Laboratoire de l'Informatique et du Parallélisme » CVP

def cvp_affine_from_chebyshev(lowerBound, upperBound):

    """

    Compute the affine parameters to map [-1, 1] to the interval given

    as argument.

    @param lowerBound (of the target interval)

    @param upperBound (of the target interval)

    @return the (scalingCoefficient, offset) tuple.

    """

    ## Check bounds consistency. Bounds must differ.

    if lowerBound >= upperBound:

        return None

    scalingCoefficient = (upperBound - lowerBound) / 2

    offset             = (lowerBound + upperBound) / 2

    return (scalingCoefficient, offset)

# End cvp_affine_from_chebyshev

#

def cvp_affine_to_chebyshev(lowerBound, upperBound):

    """

    Compute the affine parameters to map the interval given

    as argument to  [-1, 1].

    @param lowerBound (of the target interval)

    @param upperBound (of the target interval)

    @return the (scalingCoefficient, offset) tuple.

    """

    ## Check bounds consistency. Bounds must differ.

    if lowerBound >= upperBound:

        return None

    scalingCoefficient = 2 / (upperBound - lowerBound)

    ## If bounds are symmetrical with relation to 0, return 0 straight before

    #  attempting a division by 0.

    if lowerBound == -upperBound:

        offset         = 0

    else:

        offset         =  (lowerBound + upperBound) / (lowerBound - upperBound)

    return (scalingCoefficient, offset)

# End cvp_affine_to_chebyshev

#

def cvp_coefficients_bounds_projection(executablePath, arguments):

    """

    Compute the min and max of the coefficients with linear

    programming projection.

    @param executablePath: the path to the binary program;

    @param arguments: a list of arguments to be givent to the binary

    @return: the min and max coefficients value arrays (in this order).

    """

    from subprocess import check_output

    commandLine = [executablePath] + arguments

    ## Get rid of output on stderr.

    devNull = open("/dev/null", "rw")

    ## Run ther program.

    otp = check_output(commandLine, stderr=devNull)

    devNull.close()

    ## Recover the results

    bounds = sage_eval(otp)

    minCoeffsExpList = []

    maxCoeffsExpList = [] 

    print "bounds:", bounds

    for boundsPair in bounds:

        minCoeffsExpList.append(boundsPair[0])

        maxCoeffsExpList.append(boundsPair[1])

    print "minCoeffsExpList:", minCoeffsExpList

    print "maxCoeffsExpList:", maxCoeffsExpList

    return (minCoeffsExpList, maxCoeffsExpList)

# End cvp_coefficients_bounds_projection

def cvp_coefficients_bounds_projection_exps(executablePath, 

                                            arguments, 

                                            realField = None):

    """

    Compute the min and max exponents of the coefficients with linear

    programming projection.

    @param executablePath: the path to the binary program;

    @param arguments: a list of arguments to be givent to the binary

    @param realField: the realField to use for floating-point conversion.

    @return: the min and max exponents arrays (in this order).

    """

    ## If no realField is givne, fall back on RR.

    if realField is None:

        realField = RR

    (minCoeffsExpList, maxCoeffsExpList) = \

        cvp_coefficients_bounds_projection(executablePath,arguments)

    for index in xrange(0, len(minCoeffsExpList)):

        minCoeffsExpList[index] = \

            realField(minCoeffsExpList[index]).log2().floor()

        maxCoeffsExpList[index] = \

            realField(maxCoeffsExpList[index]).log2().floor()

    return (minCoeffsExpList, maxCoeffsExpList)

# End cvp_coefficients_bounds_projection_exps

def cvp_chebyshev_zeros_1k(numPoints, realField, contFracMaxErr = None):

    """

    Compute the Chebyshev zeros for some polynomial degree (numPoints, for the

    zeros) scaled to the [lowerBound, upperBound] interval.

    The list is returned as row floating-point numbers is contFracMaxErr is None.

    Otherwise elements are transformed into rational numbers. 

    """

    if numPoints < 1:

        return None

    if numPoints == 0:

        return [0]

    ## Check that realField has a "prec()" member.

    try:

        realField.prec()

    except:

        return None

    #

    zerosList = []

    ## formulas: cos (((2*i - 1) * pi) / (2*numPoints)) for i = 1,...,numPoints.

    #  If i is -1-shifted, as in the following loop, formula is 

    #  cos (((2*i + 1) * pi) / (2*numPoints)) for i = 0,...,numPoints-1.

    for index in xrange(0, numPoints):

        ## When number of points is odd (then, the Chebyshev degree is odd two), 

        #  the central point is 0. */

        if (numPoints % 2 == 1) and ((2 * index + 1) == numPoints):

            if contFracMaxErr is None:

                zerosList.append(realField(0))

            else:

                zerosList.append(0)

        else:

            currentCheby = \

                ((realField(2*index+1) * realField(pi)) / 

                 realField(2*numPoints)).cos()

            #print  "Current Cheby:", currentCheby

            ## Result is negated to have an increasing order list.

            if contFracMaxErr is None:

                zerosList.append(-currentCheby)

            else:

                zerosList.append(-currentCheby.nearby_rational(contFracMaxErr))

            #print "Relative error:", (currentCheby/zerosList[index])

    return zerosList

# End cvp_chebyshev_zeros_1k.

#

def cvp_chebyshev_zeros_for_interval(lowerBound, 

                                     upperBound, 

                                     numPoints,

                                     realField = None, 

                                     contFracMaxErr = None):

    """

    Compute the Chebyshev zeros for some polynomial degree (numPoints, for the

    zeros) scaled to the [lowerBound, upperBound] interval.

    If no contFracMaxErr argument is given, we return the list as "raw"

    floating-points. Otherwise, rational numbers are returned.

    """

    ## Arguments check.

    if lowerBound >= upperBound:

        return None

    if numPoints < 1:

        return None

    ## If no realField argument is given try to retrieve it from the 

    #  bounds. If impossible (e.g. rational bound), fall back on RR. 

    if realField is None:

        try:

            ### Force failure if parent does not have prec() member.

            lowerBound.parent().prec()

            ### If no exception is thrown, set realField.

            realField = lowerBound.parent()

        except:

            realField = RR

    #print "Real field:", realField

    ## We want "raw"floating-point nodes to only have one final rounding.

    chebyshevZerosList = \

        cvp_chebyshev_zeros_1k(numPoints, realField)

    scalingFactor, offset = cvp_affine_from_chebyshev(lowerBound, upperBound)

    for index in xrange(0, len(chebyshevZerosList)):

        chebyshevZerosList[index] = \

            chebyshevZerosList[index] * scalingFactor + offset

        if not contFracMaxErr is None:

            chebyshevZerosList[index] = \

                chebyshevZerosList[index].nearby_rational(contFracMaxErr)

    return chebyshevZerosList                                                             

# End cvp_chebyshev_zeros_for_interval.

#

                                        realField = None,

                                        contFracMaxErr = None):

    zeros of the (f-remez_d(f)) function and those of the (f-truncRemez_d(f)) 

    function over the interval. If the (f-truncRemez_d(f)) function has very 

    few zeros, compute the zeros of the derivative instead.

    ## If no realField argument is given try to retrieve it from the 

    #  bounds. If impossible (e.g. rational bound), fall back on RR. 

    if realField is None:

        try:

            ### Force failure if parent does not have prec() member.

            lowerBound.parent().prec()

            ### If no exception is thrown, set realField.

            realField = lowerBound.parent()

        except:

            realField = RR

    #print "Real field:", realField

    intervalSo = pobyso_range_from_bounds_sa_so(lowerBound, upperBound)

    ## If required, convert the numbers to rational numbers.

    if not contFracMaxErr is None:

        for index in xrange(0, len(errorFuncZerosSa)):

            errorFuncZerosSa[index] = \

                errorFuncZerosSa[index].nearby_rational(contFracMaxErr)

    return errorFuncZerosSa

    #

# End remez_all_poly_error_func_zeros

#

def cvp_remez_poly_error_func_zeros(funct, 

                                    maxDegree, 

                                    lowerBound, 

                                    upperBound,

                                    realField = None,

                                    contFracMaxErr = None):

    """

    For a given function f, a degree and an interval, compute the 

    zeros of the (f-remez_d(f)) function.

    """

    ## If no realField argument is given try to retrieve it from the 

    #  bounds. If impossible (e.g. rational bound), fall back on RR. 

    if realField is None:

        try:

            ### Force failure if parent does not have prec() member.

            lowerBound.parent().prec()

            ### If no exception is thrown, set realField.

            realField = lowerBound.parent()

        except:

            realField = RR

    #print "Real field:", realField

    funcSa         = func

    funcAsStringSa = func._assume_str().replace("_SAGE_VAR_","",100)

    #print "Function as a string:", funcAsStringSa

    ## Create the Sollya version of the function and compute the

    #  Remez approximant

    funcSo  = pobyso_parse_string(funcAsStringSa)

    pStarSo = pobyso_remez_canonical_sa_so(funcAsStringSa, 

                                           maxDegree, 

                                           lowerBound,

                                           upperBound)

    ## Compute the zeroes of the error function.

    ### Error function creation consumes both pStarSo and funcSo.

    errorFuncSo = sollya_lib_build_function_sub(pStarSo, funcSo)

    intervalSo = pobyso_range_from_bounds_sa_so(lowerBound, upperBound)

    errorFuncZerosSo = pobyso_dirty_find_zeros_so_so(errorFuncSo, intervalSo)

    pobyso_clear_obj(errorFuncSo) # Clears funcSo and pStarSo as well.

    #print "Zeroes of the error function from Sollya:"

    #pobyso_autoprint(errorFuncZerosSo)

    errorFuncZerosSa = pobyso_float_list_so_sa(errorFuncZerosSo)

    pobyso_clear_obj(errorFuncZerosSo)

    ## Sollya objects clean up.

    pobyso_clear_obj(intervalSo)

    ## Converting to rational numbers, if required.

    if not contFracMaxErr is None:

        for index in xrange(0, len(errorFuncZerosSa)):

            errorFuncZerosSa[index] = \

                errorFuncZerosSa[index].nearby_rational(contFracMaxErr)