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

                                        upperBoundSa, approxAccurSa, 

                                        precSa=None):

      output polynomial ; this error <= approxAccurSaS.

    #print"In slz_compute_polynomial_and_interval..."

    ## Change Sollya precision, if requested.

    sollyaPrecChanged = False

    (sollyaPrecSo, sollyaPrecSa) = pobyso_get_prec_so_so_sa()

    if precSa is None:

        precSa = ceil((RR('1.5') * abs(RR(approxAccurSa).log2())) / 64) * 64

        #print "Computed internal precision:", precSa

        if precSa < 192:

            precSa = 192

    if precSa != sollyaPrecSa:

        precSo = pobyso_constant_from_int_sa_so(precSa)

        pobyso_set_prec_so_so(precSo)

        sollya_lib_clear_obj(precSo)

        sollyaPrecChanged = True

    while maxErrorSa > approxAccurSa:

        shrinkFactorSa = 1 /  (maxErrorSa/approxAccurSa).log2().abs()

        #errorRatioSa = approxAccurSa/maxErrorSa

            if sollyaPrecChanged:

                pobyso_set_prec_so_so(sollyaPrecSo)

                sollya_lib_clear_obj(sollyaPrecSo)

            return None

    # End while

    if sollyaPrecChanged:

        pobyso_set_prec_so_so(sollyaPrecSo)

        sollya_lib_clear_obj(sollyaPrecSo)

def slz_compute_polynomial_and_interval_01(functionSo, degreeSo, lowerBoundSa, 

                                        upperBoundSa, approxAccurSa, 

                                        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.

    OUTPUT:

    A tuple made of 4 Sollya objects:

    - a polynomial;

    - an range (an interval, not in the sense of number given as an interval);

    - the center of the interval;

    - the maximum error in the approximation of the input functionSo by the

      output polynomial ; this error <= approxAccurSaS.

    """

    print"In slz_compute_polynomial_and_interval..."

    ## Superficial argument check.

    if lowerBoundSa > upperBoundSa:

        return None

    ## Change Sollya precision, if requested.

    if not sollyaPrecSa is None:

        sollyaPrecSo = pobyso_get_prec_so()

        pobyso_set_prec_sa_so(sollyaPrecSa)

    else:

        sollyaPrecSa = pobyso_get_prec_so_sa()

        sollyaPrecSo = None

    RRR = lowerBoundSa.parent()

    intervalShrinkConstFactorSa = RRR('0.9')

    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 > approxAccurSa:

        print "++Approximation error:", maxErrorSa.n()

        sollya_lib_clear_obj(polySo)

        sollya_lib_clear_obj(intervalCenterSo)

        sollya_lib_clear_obj(maxErrorSo)

        # Very empirical shrinking factor.

        shrinkFactorSa = 1 /  (maxErrorSa/approxAccurSa).log2().abs()

        print "Shrink factor:", \

              shrinkFactorSa.n(), \

              intervalShrinkConstFactorSa

        print

        #errorRatioSa = approxAccurSa/maxErrorSa

        #print "Error ratio: ", errorRatioSa

        # Make sure interval shrinks.

        if shrinkFactorSa > intervalShrinkConstFactorSa:

            actualShrinkFactorSa = intervalShrinkConstFactorSa

            #print "Fixed"

        else:

            actualShrinkFactorSa = shrinkFactorSa

            #print "Computed",shrinkFactorSa,maxErrorSa

            #print shrinkFactorSa, maxErrorSa

        #print "Shrink factor", actualShrinkFactorSa

        currentUpperBoundSa = currentLowerBoundSa + \

                                (currentUpperBoundSa - currentLowerBoundSa) * \

                                actualShrinkFactorSa

        #print "Current upper bound:", currentUpperBoundSa

        sollya_lib_clear_obj(currentRangeSo)

        # Check what is left with the bounds.

        if currentUpperBoundSa <= currentLowerBoundSa or \

              currentUpperBoundSa == currentLowerBoundSa.nextabove():

            sollya_lib_clear_obj(absoluteErrorTypeSo)

            print "Can't find an interval."

            print "Use either or both a higher polynomial degree or a higher",

            print "internal precision."

            print "Aborting!"

            if not sollyaPrecSo is None:

                pobyso_set_prec_so_so(sollyaPrecSo)

                sollya_lib_clear_obj(sollyaPrecSo)

            return None

        currentRangeSo = pobyso_bounds_to_range_sa_so(currentLowerBoundSa, 

                                                      currentUpperBoundSa)

        # print "New interval:",

        # pobyso_autoprint(currentRangeSo)

        #print "Second Taylor expansion call."

        (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)

        #print "Max errorSo:",

        #pobyso_autoprint(maxErrorSo)

        maxErrorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo)

        #print "Max errorSa:", maxErrorSa

        #print "Sollya prec:",

        #pobyso_autoprint(sollya_lib_get_prec(None))

    # End while

    sollya_lib_clear_obj(absoluteErrorTypeSo)

    itpSo           = pobyso_constant_from_int_sa_so(floor(sollyaPrecSa/3))

    ftpSo           = pobyso_constant_from_int_sa_so(floor(2*sollyaPrecSa/3))

    maxPrecSo       = pobyso_constant_from_int_sa_so(sollyaPrecSa)

    approxAccurSo   = pobyso_constant_sa_so(RR(approxAccurSa))

    print "About to call polynomial rounding with:"

    print "polySo: ",           ; pobyso_autoprint(polySo)

    print "functionSo: ",       ; pobyso_autoprint(functionSo)

    print "intervalCenterSo: ", ; pobyso_autoprint(intervalCenterSo)

    print "currentRangeSo: ",   ; pobyso_autoprint(currentRangeSo)

    print "itpSo: ",            ; pobyso_autoprint(itpSo)

    print "ftpSo: ",            ; pobyso_autoprint(ftpSo)

    print "maxPrecSo: ",        ; pobyso_autoprint(maxPrecSo)

    print "approxAccurSo: ",    ; pobyso_autoprint(approxAccurSo)

    (roundedPolySo, roundedPolyMaxErrSo) = \

    pobyso_polynomial_coefficients_progressive_round_so_so(polySo,

                                                           functionSo,

                                                           intervalCenterSo,

                                                           currentRangeSo,

                                                           itpSo,

                                                           ftpSo,

                                                           maxPrecSo,

                                                           approxAccurSo)

    sollya_lib_clear_obj(polySo)

    sollya_lib_clear_obj(maxErrorSo)

    sollya_lib_clear_obj(itpSo)

    sollya_lib_clear_obj(ftpSo)

    sollya_lib_clear_obj(approxAccurSo)

    if not sollyaPrecSo is None:

        pobyso_set_prec_so_so(sollyaPrecSo)

        sollya_lib_clear_obj(sollyaPrecSo)

    print "1: ", ; pobyso_autoprint(roundedPolySo)

    print "2: ", ; pobyso_autoprint(currentRangeSo)

    print "3: ", ; pobyso_autoprint(intervalCenterSo)

    print "4: ", ; pobyso_autoprint(roundedPolyMaxErrSo)

    return (roundedPolySo, currentRangeSo, intervalCenterSo, roundedPolyMaxErrSo)

# End slz_compute_polynomial_and_interval_01

    print "Approximation accuracy: ", RR(approxAccurSa)

                                        approxAccurSa, internalSollyaPrecSa)

            (errorSa <= approxAccurSa):

        currentErrorRatio = approxAccurSa / errorSa

    The function returns a 5 elements tuple: formed with all the 

    polynomialSa = pobyso_float_poly_so_sa(polyRangeCenterErrorSo[0])

    #print "Polynomial after first conversion: ", pobyso_autoprint(polyRangeCenterErrorSo[1])

    polynomialChangedVarSa = pobyso_float_poly_so_sa(polyRangeCenterErrorSo[1])