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

def slz_compute_polynomial_and_interval(functionSo, degreeSo, lowerBoundSa, 

                                        upperBoundSa, approxPrecSa, 

                                        sollyaPrecSa=None):

    """

    Under the assumptions listed for slz_get_intervals_and_polynomials, compute

    a polynomial that approximates the function on a an interval starting

    at lowerBoundSa and finishing at a value that guarantees that the polynomial

    approximates with the expected precision.

    The interval upper bound is lowered until the expected approximation 

    precision is reached.

    The polynomial, the bounds, the center of the interval and the error

    are returned.

    """

    RRR = lowerBoundSa.parent()

    #goldenRatioSa = RRR(5.sqrt() / 2 - 1/2)

    #intervalShrinkConstFactorSa = goldenRatioSa

    intervalShrinkConstFactorSa = RRR('0.5')

    absoluteErrorTypeSo = pobyso_absolute_so_so()

    currentRangeSo = pobyso_bounds_to_range_sa_so(lowerBoundSa, upperBoundSa)

    currentUpperBoundSa = upperBoundSa

    currentLowerBoundSa = lowerBoundSa

    # What we want here is the polynomial without the variable change, 

    # since our actual variable will be x-intervalCenter defined over the 

    # domain [lowerBound-intervalCenter , upperBound-intervalCenter]. 

    (polySo, intervalCenterSo, maxErrorSo) = \

        pobyso_taylor_expansion_no_change_var_so_so(functionSo, degreeSo, 

                                                    currentRangeSo, 

                                                    absoluteErrorTypeSo)

    maxErrorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo)

    while maxErrorSa > approxPrecSa:

        sollya_lib_clear_obj(maxErrorSo)

        errorRatioSa = 1/(maxErrorSa/approxPrecSa).log2()

        #print "Error ratio: ", errorRatioSa

        if errorRatioSa < intervalShrinkConstFactorSa:

            #currentUpperBoundSa = currentLowerBoundSa + (currentUpperBoundSa - currentLowerBoundSa) * errorRatioSa

            currentUpperBoundSa = currentLowerBoundSa + \

                                  (currentUpperBoundSa - currentLowerBoundSa) * \

                                   intervalShrinkConstFactorSa

        else:

            currentUpperBoundSa = currentLowerBoundSa + \

                                  (currentUpperBoundSa - currentLowerBoundSa) * \

                                  intervalShrinkConstFactorSa

        #print lowerBoundSa, currentUpperBoundSa

        sollya_lib_clear_obj(currentRangeSo)

        sollya_lib_clear_obj(polySo)

        currentRangeSo = pobyso_bounds_to_range_sa_so(currentLowerBoundSa, 

                                                      currentUpperBoundSa)

        (polySo, intervalCenterSo, maxErrorSo) = \

            pobyso_taylor_expansion_no_change_var_so_so(functionSo, degreeSo, 

                                                        currentRangeSo, 

                                                        absoluteErrorTypeSo)

        #maxErrorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo, RRR)

        maxErrorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo)

    sollya_lib_clear_obj(absoluteErrorTypeSo)

    return((polySo, currentRangeSo, intervalCenterSo, maxErrorSo))

    # End slz_compute_polynomial_and_interval

def slz_get_intervals_and_polynomials(functionSa, degreeSa, lowerBoundSa, 

                                      upperBoundSa, floatingPointPrecSa, 

                                      internalSollyaPrecSa):

    """

    Under the assumption that:

    - functionSa is monotonic on the [lowerBoundSa, upperBoundSa] interval;

    - lowerBound and upperBound belong to the same binade.

    from a:

    - function;

    - a degree

    - a pair of bounds;

    - the floating-point precision we work on;

    - the internal Sollya precision;

    compute a list of tuples made of:

    - a polynomial approximating the function (a Sollya object);

    - the range for which the polynomial approximates the function 

      (a Sollya object);

    - the center of the interval (a Sollya object);

    - the approximation error (a Sage object). 

      with the error given as the last element (a Sage object);

    The initial interval is, possibly, splitted into smaller intervals.

    It return a list of tuples, each made of:

    - a polynomial;

    - the approximation interval;

    - the center, x0,  of the interval (the polynomial is defined as p(x-x0));

    - the corresponding approximation error.

    """

    # Scalling the domain -> [1,2[.

    # Notice the clumsy notation for log2.

    domainScalingFactorSa = floor(lowerBound.log2()) + 1

    print "domainScalingFactor for argument :", domainScalingFactorSa.n()

    ff(x) = f(x * domainScalingFactorSa)

    scaledLowerBoundSa = lowerBoundSa/domainScalingFactorSa

    scaledUpperBoundSa = upperBoundSa/domainScalingFactorSa

    print 'ff:', ff, "- Domain:", scaledLowerBoundSa, scaledUpperBoundSa

    #

    # Scalling the image -> [1,2[.

    flb = f(lowerBoundSa).n()

    fub = f(upperBoundSa).n()

    if flb <= fub: # Increasing

        imageBinadeBottom = floor(flb.log2())

    else: # Decreasing

        imageBinadeBottom = floor(fub.log2())

    print 'ff:', ff, '- Image:', flb, fub, imageBinadeBottom

    #

    resultArray = []

    #

    approxPrecSa = 1/(2^(floatingPointPrecSa + 1))

    print "Approximation precision: ", RR(approxPrecSa)

    # Prepare the arguments for the Taylor expansion computation with Sollya.

    functionSo = pobyso_parse_string_sa_so(functionSa._assume_str())

    degreeSo = pobyso_constant_from_int_sa_so(degreeSa)

    scaledBoundsSo = pobyso_bounds_to_range_sa_so(scaledLowerBoundSa, 

                                                  scaledUpperBoundSa)

    absoluteErrorTypeSo = pobyso_absolute_so_so()

    # Compute the first Taylor expansion.

    (polySo, boundsSo, intervalCenterSo, maxErrorSo) = \

        slz_compute_polynomial_and_interval(functionSo, degreeSo,

                                        scaledLowerBoundSa, scaledUpperBoundSa,

                                        approxPrecSa, internalSollyaPrecSa)

    resultArray.append((polySo, boundsSo, intervalCenterSo, maxErrorSo))

    realIntervalField = RealIntervalField(max(lowerBoundSa.parent().precision(),

                                              upperBoundSa.parent().precision()))

    boundsSa = pobyso_range_to_interval_so_sa(boundsSo, realIntervalField)

    # Compute the other expansions.

    while boundsSa.endpoints()[1] < scaledUpperBoundSa:

        currentScaledLowerBoundSa = boundsSa.endpoints()[1]

        (polySo, boundsSo, intervalCenterSo, maxErrorSo) = \

            slz_compute_polynomial_and_interval(functionSo, degreeSo,

                                            currentScaledLowerBoundSa, 

                                            scaledUpperBoundSa, approxPrecSa, 

                                            internalSollyaPrecSa)

        resultArray.append((polySo, boundsSo, intervalCenterSo, maxErrorSo))

        boundsSa = pobyso_range_to_interval_so_sa(boundsSo, realIntervalField)

    sollya_lib_clear_obj(functionSo)

    sollya_lib_clear_obj(degreeSo)

    sollya_lib_clear_obj(scaledBoundsSo)

    sollya_lib_clear_obj(absoluteErrorTypeSo)

    return(resultArray)

    # End slz_get_intervals_and_polynomials

def slz_interval_scaling_expression(boundsInterval, varName):

    """

    Compute the scaling expression to map an interval that span only

    a binade to [1, 2)

    """

    if abs(boundsInterval.endpoints()[0]) < 1:

        if boundsInterval.endpoints()[0] >= 0:

            scalingCoeff = 2^(-floor(boundsInterval.endpoints()[0].log2()))

            return(scalingCoeff * eval(varName))

        else:

            scalingCoeff = 2^(-floor((-boundsInterval.endpoints()[1]).log2()))

            scalingOffset = -ceil(scalingCoeff * boundsInterval.endpoints()[0])

            return(scalingCoeff * eval(varName) + scalingOffset)

    else: 

        if boundsInterval.endpoints()[0] >= 0:

            scalingCoeff = 2^(-floor(boundsInterval.endpoints()[0].log2()))

            scalingOffset = 0

            return(scalingCoeff * eval(varName))

        else:

            scalingCoeff = 2^(-floor((-boundsInterval.endpoints()[1]).log2()))

            scalingOffset = floor(-(scalingCoeff * boundsInterval.endpoints()[1]) + 2)

            return(scalingCoeff * eval(varName) + scalingOffset)

def slz_polynomial_and_interval_to_sage(polyRangeCenterErrorSo):

    polynomialSa = pobyso_get_poly_so_sa(polyRangeCenterErrorSo[0])

    intervalSa = \

            pobyso_get_interval_from_range_so_sa(polyRangeCenterErrorSo[1])

    centerSa = \

            pobyso_get_constant_as_rn_with_rf_so_sa(polyRangeCenterErrorSo[2])

    errorSa = \

            pobyso_get_constant_as_rn_with_rf_so_sa(polyRangeCenterErrorSo[3])

    return((polynomialSa, intervalSa, centerSa, errorSa))

    # End slz_polynomial_and_interval_to_sage