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

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 bounds for which the polynomial approximates the function 

      (a Sage object);

    - the center of the interval;

    - the approximation error (a Sage object). 

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

    The initial interval is, possibly, splitted into a list of smaller interval

    each associated to an approximation polynomial and a 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.

    functionSo = pobyso_parse_string_sa_so(functionSa._assume_str())

    degreeSo = pobyso_constant_from_int_sa_so(degreeSa)

    scaledBoundsSo = pobyso_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_so_sa(boundsSo, realIntervalField)

    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_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_compute_polynomial_and_interval(functionSo, degreeSo, lowerBoundSa, 

                                        upperBoundSa, approxPrecSa, 

                                        sollyaPrecSa=None):

    RRR = lowerBoundSa.parent()

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

    #intervalShrinkConstFactorSa = goldenRatioSa

    intervalShrinkConstFactorSa = RRR('0.5')

    absoluteErrorTypeSo = pobyso_absolute_so_so()

    currentRangeSo = pobyso_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_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_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