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

                                                           maxErrSo,

                                                           debug=False):

    if debug:

        print "Input arguments:"

        pobyso_autoprint(polySo)

        pobyso_autoprint(funcSo)

        pobyso_autoprint(icSo)

        pobyso_autoprint(intervalSo)

        pobyso_autoprint(itpSo)

        pobyso_autoprint(ftpSo)

        pobyso_autoprint(maxPrecSo)

        pobyso_autoprint(maxErrSo)

        print "________________"

    if debug:

        print "degreeSa:", degreeSa

        print "ratio:", ratio

        print "itpsSa:", itpSa

        print "ftpSa:", ftpSa

        print "maxPrecSa:", maxPrecSa

        print "maxErrSa:", maxErrSa

    #print "About to enter the while loop..."

            if debug:

                print "Index:", indexSa, " - Target precision:", 

                pobyso_autoprint(ctpSo)

        if debug:

            print "Infnorm (Sollya):", 

            pobyso_autoprint(infNormSo)

            if debug:

                print "Error is too large."

                if debug:

                    print "Enlarging prec."

                if debug:

                    print "Exit with the last before last polynomial."

            if debug:

                print "Error is too small"

            if cerrSa <= (maxErrSa/2):

                if debug: 

                    print "Shrinking prec."

            - roots in i are searched for;

def srs_run_SLZ_v05(inputFunction,

                    inputLowerBound,

                    inputUpperBound,

                    alpha,

                    degree,

                    precision,

                    emin,

                    emax,

                    targetHardnessToRound,

                    debug = False):

    """

    Changes from V4:

        Approximation polynomial has coefficients rounded.

    Changes from V3:

        Root search is changed again:

            - only resultants in i are computed;

            - roots in i are searched for;

            - if any, they are tested for hardness-to-round.

    Changes from V2:

        Root search is changed:

            - we compute the resultants in i and in t;

            - we compute the roots set of each of these resultants;

            - we combine all the possible pairs between the two sets;

            - we check these pairs in polynomials for correctness.

    Changes from V1: 

        1- check for roots as soon as a resultant is computed;

        2- once a non null resultant is found, check for roots;

        3- constant resultant == no root.

    """

    if debug:

        print "Function                :", inputFunction

        print "Lower bound             :", inputLowerBound

        print "Upper bounds            :", inputUpperBound

        print "Alpha                   :", alpha

        print "Degree                  :", degree

        print "Precision               :", precision

        print "Emin                    :", emin

        print "Emax                    :", emax

        print "Target hardness-to-round:", targetHardnessToRound

        print

    ## Important constants.

    ### Stretch the interval if no error happens.

    noErrorIntervalStretch = 1 + 2^(-5)

    ### If no vector validates the Coppersmith condition, shrink the interval

    #   by the following factor.

    noCoppersmithIntervalShrink = 1/2

    ### If only (or at least) one vector validates the Coppersmith condition, 

    #   shrink the interval by the following factor.

    oneCoppersmithIntervalShrink = 3/4

    #### If only null resultants are found, shrink the interval by the 

    #    following factor.

    onlyNullResultantsShrink     = 3/4

    ## Structures.

    RRR         = RealField(precision)

    RRIF        = RealIntervalField(precision)

    ## Converting input bound into the "right" field.

    lowerBound = RRR(inputLowerBound)

    upperBound = RRR(inputUpperBound)

    ## Before going any further, check domain and image binade conditions.

    print inputFunction(1).n()

    output = slz_fix_bounds_for_binades(lowerBound, upperBound, inputFunction)

    if output is None:

        print "Invalid domain/image binades. Domain:",\

        lowerBound, upperBound, "Images:", \

        inputFunction(lowerBound), inputFunction(upperBound)

        raise Exception("Invalid domain/image binades.")

    lb = output[0] ; ub = output[1]

    if lb != lowerBound or ub != upperBound:

        print "lb:", lb, " - ub:", ub

        print "Invalid domain/image binades. Domain:",\

        lowerBound, upperBound, "Images:", \

        inputFunction(lowerBound), inputFunction(upperBound)

        raise Exception("Invalid domain/image binades.")

    #

    ## Progam initialization

    ### Approximation polynomial accuracy and hardness to round.

    polyApproxAccur           = 2^(-(targetHardnessToRound + 1))

    polyTargetHardnessToRound = targetHardnessToRound + 1

    ### Significand to integer conversion ratio.

    toIntegerFactor           = 2^(precision-1)

    print "Polynomial approximation required accuracy:", polyApproxAccur.n()

    ### Variables and rings for polynomials and root searching.

    i=var('i')

    t=var('t')

    inputFunctionVariable = inputFunction.variables()[0]

    function = inputFunction.subs({inputFunctionVariable:i})

    #  Polynomial Rings over the integers, for root finding.

    Zi = ZZ[i]

    Zt = ZZ[t]

    Zit = ZZ[i,t]

    ## Number of iterations limit.

    maxIter = 100000

    #

    ## Compute the scaled function and the degree, in their Sollya version 

    #  once for all.

    (scaledf, sdlb, sdub, silb, siub) = \

        slz_compute_scaled_function(function, lowerBound, upperBound, precision)

    print "Scaled function:", scaledf._assume_str().replace('_SAGE_VAR_', '')

    scaledfSo = sollya_lib_parse_string(scaledf._assume_str().replace('_SAGE_VAR_', ''))

    degreeSo  = pobyso_constant_from_int_sa_so(degree)

    #

    ## Compute the scaling. boundsIntervalRifSa defined out of the loops.

    domainBoundsInterval   = RRIF(lowerBound, upperBound)

    (unscalingFunction, scalingFunction) = \

        slz_interval_scaling_expression(domainBoundsInterval, i)

    #print scalingFunction, unscalingFunction

    ## Set the Sollya internal precision (with an arbitrary minimum of 192).

    internalSollyaPrec = ceil((RR('1.5') * targetHardnessToRound) / 64) * 64

    if internalSollyaPrec < 192:

        internalSollyaPrec = 192

    pobyso_set_prec_sa_so(internalSollyaPrec)

    print "Sollya internal precision:", internalSollyaPrec

    ## Some variables.

    ### General variables

    lb                  = sdlb

    ub                  = sdub

    nbw                 = 0

    intervalUlp         = ub.ulp()

    #### Will be set by slz_interval_and_polynomila_to_sage.

    ic                  = 0 

    icAsInt             = 0    # Set from ic.

    solutionsSet        = set()

    tsErrorWidth        = []

    csErrorVectors      = []

    csVectorsResultants = []

    floatP              = 0  # Taylor polynomial.

    floatPcv            = 0  # Ditto with variable change.

    intvl               = "" # Taylor interval

    terr                = 0  # Taylor error.

    iterCount = 0

    htrnSet = set()

    ### Timers and counters.

    wallTimeStart                   = 0

    cpuTimeStart                    = 0

    taylCondFailedCount             = 0

    coppCondFailedCount             = 0

    resultCondFailedCount           = 0

    coppCondFailed                  = False

    resultCondFailed                = False

    globalResultsList               = []

    basisConstructionsCount         = 0

    basisConstructionsFullTime      = 0

    basisConstructionTime           = 0

    reductionsCount                 = 0

    reductionsFullTime              = 0

    reductionTime                   = 0

    resultantsComputationsCount     = 0

    resultantsComputationsFullTime  = 0

    resultantsComputationTime       = 0

    rootsComputationsCount          = 0

    rootsComputationsFullTime       = 0

    rootsComputationTime            = 0

    print

    ## Global times are started here.

    wallTimeStart                   = walltime()

    cpuTimeStart                    = cputime()

    ## Main loop.

    while True:

        if lb >= sdub:

            print "Lower bound reached upper bound."

            break

        if iterCount == maxIter:

            print "Reached maxIter. Aborting"

            break

        iterCount += 1

        print "[", lb, ",", ub, "]", ((ub - lb) / intervalUlp).log2().n(), \

            "log2(numbers)." 

        ### Compute a Sollya polynomial that will honor the Taylor condition.

        prceSo = slz_compute_polynomial_and_interval_01(scaledfSo,

                                                        degreeSo,

                                                        lb,

                                                        ub,

                                                        polyApproxAccur)

        ### Convert back the data into Sage space.                         

        (floatP, floatPcv, intvl, ic, terr) = \

            slz_interval_and_polynomial_to_sage((prceSo[0], prceSo[0],

                                                 prceSo[1], prceSo[2], 

                                                 prceSo[3]))

        intvl = RRIF(intvl)

        ## Clean-up Sollya stuff.

        for elem in prceSo:

            sollya_lib_clear_obj(elem)

        #print  floatP, floatPcv, intvl, ic, terr

        #print floatP

        #print intvl.endpoints()[0].n(), \

        #      ic.n(),

        #intvl.endpoints()[1].n()

        ### Check returned data.

        #### Is approximation error OK?

        if terr > polyApproxAccur:

            exceptionErrorMess  = \

                "Approximation failed - computed error:" + \

                str(terr) + " - target error: "

            exceptionErrorMess += \

                str(polyApproxAccur) + ". Aborting!"

            raise Exception(exceptionErrorMess)

        #### Is lower bound OK?

        if lb != intvl.endpoints()[0]:

            exceptionErrorMess =  "Wrong lower bound:" + \

                               str(lb) + ". Aborting!"

            raise Exception(exceptionErrorMess)

        #### Set upper bound.

        if ub > intvl.endpoints()[1]:

            ub = intvl.endpoints()[1]

            print "[", lb, ",", ub, "]", ((ub - lb) / intervalUlp).log2().n(), \

            "log2(numbers)." 

            taylCondFailedCount += 1

        #### Is interval not degenerate?

        if lb >= ub:

            exceptionErrorMess = "Degenerate interval: " + \

                                "lowerBound(" + str(lb) +\

                                ")>= upperBound(" + str(ub) + \

                                "). Aborting!"

            raise Exception(exceptionErrorMess)

        #### Is interval center ok?

        if ic <= lb or ic >= ub:

            exceptionErrorMess =  "Invalid interval center for " + \

                                str(lb) + ',' + str(ic) + ',' +  \

                                str(ub) + ". Aborting!"

            raise Exception(exceptionErrorMess)

        ##### Current interval width and reset future interval width.

        bw = ub - lb

        nbw = 0

        icAsInt = int(ic * toIntegerFactor)

        #### The following ratio is always >= 1. In case we may want to

        #    enlarge the interval

        curTaylErrRat = polyApproxAccur / terr

        ### Make the  integral transformations.

        #### Bounds and interval center.

        intIc = int(ic * toIntegerFactor)

        intLb = int(lb * toIntegerFactor) - intIc

        intUb = int(ub * toIntegerFactor) - intIc

        #

        #### Polynomials

        basisConstructionTime         = cputime()

        ##### To a polynomial with rational coefficients with rational arguments

        ratRatP = slz_float_poly_of_float_to_rat_poly_of_rat_pow_two(floatP)

        ##### To a polynomial with rational coefficients with integer arguments

        ratIntP = \

            slz_rat_poly_of_rat_to_rat_poly_of_int(ratRatP, precision)

        #####  Ultimately a multivariate polynomial with integer coefficients  

        #      with integer arguments.

        coppersmithTuple = \

            slz_rat_poly_of_int_to_poly_for_coppersmith(ratIntP, 

                                                        precision, 

                                                        targetHardnessToRound, 

                                                        i, t)

        #### Recover Coppersmith information.

        intIntP = coppersmithTuple[0]

        N = coppersmithTuple[1]

        nAtAlpha = N^alpha

        tBound = coppersmithTuple[2]

        leastCommonMultiple = coppersmithTuple[3]

        iBound = max(abs(intLb),abs(intUb))

        basisConstructionsFullTime        += cputime(basisConstructionTime)

        basisConstructionsCount           += 1

        """

        #### Compute the matrix to reduce.

        matrixToReduce = slz_compute_initial_lattice_matrix(intIntP,

                                                            alpha,

                                                            N,

                                                            iBound,

                                                            tBound)

        matrixFile = file('/tmp/matrixToReduce.txt', 'w')

        for row in matrixToReduce.rows():

            matrixFile.write(str(row) + "\n")

        matrixFile.close()

        raise Exception("Deliberate stop here.")

        """

        reductionTime                     = cputime()

        #### Compute the reduced polynomials.

        ccReducedPolynomialsList =  \

                slz_compute_coppersmith_reduced_polynomials(intIntP, 

                                                            alpha, 

                                                            N, 

                                                            iBound, 

                                                            tBound)

        if ccReducedPolynomialsList is None:

            raise Exception("Reduction failed.")

        reductionsFullTime    += cputime(reductionTime)

        reductionsCount       += 1

        if len(ccReducedPolynomialsList) < 2:

            print "Nothing to form resultants with."

            print

            coppCondFailedCount += 1

            coppCondFailed = True

            ##### Apply a different shrink factor according to 

            #  the number of compliant polynomials.

            if len(ccReducedPolynomialsList) == 0:

                ub = lb + bw * noCoppersmithIntervalShrink

            else: # At least one compliant polynomial.

                ub = lb + bw * oneCoppersmithIntervalShrink

            if ub > sdub:

                ub = sdub

            if lb == ub:

                raise Exception("Cant shrink interval \

                anymore to get Coppersmith condition.")

            nbw = 0

            continue

        #### We have at least two polynomials.

        #    Let us try to compute resultants.

        #    For each resultant computed, go for the solutions.

        ##### Build the pairs list.

        polyPairsList           = []

        for polyOuterIndex in xrange(0, len(ccReducedPolynomialsList) - 1):

            for polyInnerIndex in xrange(polyOuterIndex+1, 

                                         len(ccReducedPolynomialsList)):

                polyPairsList.append((ccReducedPolynomialsList[polyOuterIndex],

                                      ccReducedPolynomialsList[polyInnerIndex]))

        #### Actual root search.

        iRootsSet           = set()

        hasNonNullResultant = False

        for polyPair in polyPairsList:

            resultantsComputationTime   = cputime()

            currentResultantI = \

                slz_resultant(polyPair[0],

                              polyPair[1],

                              t)

            resultantsComputationsCount    += 1

            resultantsComputationsFullTime +=  \

                cputime(resultantsComputationTime)

            #### Function slz_resultant returns None both for None and O

            #    resultants.

            if currentResultantI is None:

                print "Nul resultant"

                continue # Next polyPair.

            ## We deleted the currentResultantI computation.

            #### We have a non null resultant. From now on, whatever this

            #    root search yields, no extra root search is necessary.

            hasNonNullResultant = True

            #### A constant resultant leads to no root. Root search is done.

            if currentResultantI.degree() < 1:

                print "Resultant is constant:", currentResultantI

                break # There is no root.

            #### Actual iroots computation.

            rootsComputationTime        = cputime()

            iRootsList = Zi(currentResultantI).roots()

            rootsComputationsCount      +=  1

            rootsComputationsFullTime   =   cputime(rootsComputationTime)

            if len(iRootsList) == 0:

                print "No roots in \"i\"."

                break # No roots in i.

            else:

                for iRoot in iRootsList:

                    # A root is given as a (value, multiplicity) tuple.

                    iRootsSet.add(iRoot[0])

        # End loop for polyPair in polyParsList. We only loop again if a 

        # None or zero resultant is found.

        #### Prepare for results for the current interval..

        intervalResultsList = []

        intervalResultsList.append((lb, ub))

        #### Check roots.

        rootsResultsList = []

        for iRoot in iRootsSet:

            specificRootResultsList = []

            failingBounds           = []

            # Root qualifies for modular equation, test it for hardness to round.

            hardToRoundCaseAsFloat = RRR((icAsInt + iRoot) / toIntegerFactor)

            #print "Before unscaling:", hardToRoundCaseAsFloat.n(prec=precision)

            #print scalingFunction

            scaledHardToRoundCaseAsFloat = \

                scalingFunction(hardToRoundCaseAsFloat) 

            print "Candidate HTRNc at x =", \

                scaledHardToRoundCaseAsFloat.n().str(base=2),

            if slz_is_htrn(scaledHardToRoundCaseAsFloat,

                           function,

                           2^-(targetHardnessToRound),

                           RRR):

                print hardToRoundCaseAsFloat, "is HTRN case."

                specificRootResultsList.append(hardToRoundCaseAsFloat.n().str(base=2))

                if lb <= hardToRoundCaseAsFloat and hardToRoundCaseAsFloat <= ub:

                    print "Found in interval."

                else:

                    print "Found out of interval."

                # Check the i root is within the i bound.

                if abs(iRoot) > iBound:

                    print "IRoot", iRoot, "is out of bounds for modular equation."

                    print "i bound:", iBound

                    failingBounds.append('i')

                    failingBounds.append(iRoot)

                    failingBounds.append(iBound)

                if len(failingBounds) > 0:

                    specificRootResultsList.append(failingBounds)

            else: # From slz_is_htrn...

                print "is not an HTRN case."

            if len(specificRootResultsList) > 0:

                rootsResultsList.append(specificRootResultsList)

        if len(rootsResultsList) > 0:

            intervalResultsList.append(rootsResultsList)

        ### Check if a non null resultant was found. If not shrink the interval.

        if not hasNonNullResultant:

            print "Only null resultants for this reduction, shrinking interval."

            resultCondFailed      =  True

            resultCondFailedCount += 1

            ### Shrink interval for next iteration.

            ub = lb + bw * onlyNullResultantsShrink

            if ub > sdub:

                ub = sdub

            nbw = 0

            continue

        #### An intervalResultsList has at least the bounds.

        globalResultsList.append(intervalResultsList)   

        #### Compute an incremented width for next upper bound, only

        #    if not Coppersmith condition nor resultant condition

        #    failed at the previous run. 

        if not coppCondFailed and not resultCondFailed:

            nbw = noErrorIntervalStretch * bw

        else:

            nbw = bw

        ##### Reset the failure flags. They will be raised

        #     again if needed.

        coppCondFailed   = False

        resultCondFailed = False

        #### For next iteration (at end of loop)

        #print "nbw:", nbw

        lb  = ub

        ub += nbw     

        if ub > sdub:

            ub = sdub

        print

    # End while True

    ## Main loop just ended.

    globalWallTime = walltime(wallTimeStart)

    globalCpuTime  = cputime(cpuTimeStart)

    ## Output results

    print ; print "Intervals and HTRNs" ; print

    for intervalResultsList in globalResultsList:

        print "[", intervalResultsList[0][0], ",",intervalResultsList[0][1], "]",

        if len(intervalResultsList) > 1:

            rootsResultsList = intervalResultsList[1]

            for specificRootResultsList in rootsResultsList:

                print "\t", specificRootResultsList[0],

                if len(specificRootResultsList) > 1:

                    print specificRootResultsList[1],

        print ; print

    #print globalResultsList

    #

    print "Timers and counters"

    print

    print "Number of iterations:", iterCount

    print "Taylor condition failures:", taylCondFailedCount

    print "Coppersmith condition failures:", coppCondFailedCount

    print "Resultant condition failures:", resultCondFailedCount

    print "Iterations count: ", iterCount

    print "Number of intervals:", len(globalResultsList)

    print "Number of basis constructions:", basisConstructionsCount 

    print "Total CPU time spent in basis constructions:", \

        basisConstructionsFullTime

    if basisConstructionsCount != 0:

        print "Average basis construction CPU time:", \

            basisConstructionsFullTime/basisConstructionsCount

    print "Number of reductions:", reductionsCount

    print "Total CPU time spent in reductions:", reductionsFullTime

    if reductionsCount != 0:

        print "Average reduction CPU time:", \

            reductionsFullTime/reductionsCount

    print "Number of resultants computation rounds:", \

        resultantsComputationsCount

    print "Total CPU time spent in resultants computation rounds:", \

        resultantsComputationsFullTime

    if resultantsComputationsCount != 0:

        print "Average resultants computation round CPU time:", \

            resultantsComputationsFullTime/resultantsComputationsCount

    print "Number of root finding rounds:", rootsComputationsCount

    print "Total CPU time spent in roots finding rounds:", \

        rootsComputationsFullTime

    if rootsComputationsCount != 0:

        print "Average roots finding round CPU time:", \

            rootsComputationsFullTime/rootsComputationsCount

    print "Global Wall time:", globalWallTime

    print "Global CPU time:", globalCpuTime

    ## Output counters

# End srs_runSLZ-v05

precision = 113

intervalCenter      = RRR(3/8)

intervalRadiusInUlp = 2^49 + 2^45

intervalRadius      = RRR(2^(-35))          

                inputLowerBound         = intervalCenter - intervalRadius, 

                inputUpperBound         = intervalCenter + intervalRadius, 

                alpha                   = 6, 

                degree                  = 18, 

                precision               = precision, 

                emin                    = -16382, 

                emax                    = 16383, 

                targetHardnessToRound   = precision * 6, 

                debug                   = True)

def slz_compute_initial_lattice_matrix(inputPolynomial,

                                       alpha,

                                       N,

                                       iBound,

                                       tBound):

    """

    For a given set of arguments (see below), compute the initial lattice

    that could be reduced. 

    INPUT:

    - "inputPolynomial" -- (no default) a bivariate integer polynomial;

    - "alpha" -- the alpha parameter of the Coppersmith algorithm;

    - "N" -- the modulus;

    - "iBound" -- the bound on the first variable;

    - "tBound" -- the bound on the second variable.

    OUTPUT:

    A list of bivariate integer polynomial obtained using the Coppersmith

    algorithm. The polynomials correspond to the rows of the LLL-reduce

    reduced base that comply with the Coppersmith condition.

    """

    # Arguments check.

    if iBound == 0 or tBound == 0:

        return None

    # End arguments check.

    nAtAlpha = N^alpha

    ## Building polynomials for matrix.

    polyRing   = inputPolynomial.parent()

    # Whatever the 2 variables are actually called, we call them

    # 'i' and 't' in all the variable names.

    (iVariable, tVariable) = inputPolynomial.variables()[:2]

    #print polyVars[0], type(polyVars[0])

    initialPolynomial = inputPolynomial.subs({iVariable:iVariable * iBound,

                                              tVariable:tVariable * tBound})

    polynomialsList = \

        spo_polynomial_to_polynomials_list_8(initialPolynomial,

                                             alpha,

                                             N,

                                             iBound,

                                             tBound,

                                             0)

    #print "Polynomials list:", polynomialsList

    ## Building the proto matrix.

    knownMonomials = []

    protoMatrix    = []

    for poly in polynomialsList:

        spo_add_polynomial_coeffs_to_matrix_row(poly,

                                                knownMonomials,

                                                protoMatrix,

                                                0)

    matrixToReduce = spo_proto_to_row_matrix(protoMatrix)

    return matrixToReduce

# End slz_compute_initial_lattice_matrix.

                                        sollyaPrecSa=None, debug=False):

    """

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

    if debug:

        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)

    if debug:

        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

def slz_compute_polynomial_and_interval_02(functionSo, degreeSo, lowerBoundSa, 

                                        upperBoundSa, approxAccurSa, 

# End slz_compute_polynomial_and_interval_02