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

r"""

Sage core functions needed for the implementation of SLZ.

AUTHORS:

- S.T. (2013-08): initial version

Examples:

TODO::

"""

print "sageSLZ loading..."

#

import inspect

#

def slz_compute_binade(number):

    """"

    For a given number, compute the "binade" that is integer m such that

    2^m <= number < 2^(m+1). If number == 0 return None.

    """

    # Checking the parameter.

    # The exception construction is used to detect if number is a RealNumber

    # since not all numbers have

    # the mro() method. sage.rings.real_mpfr.RealNumber do.

    try:

        classTree = [number.__class__] + number.mro()

        # If the number is not a RealNumber (or offspring thereof) try

        # to transform it.

        if not sage.rings.real_mpfr.RealNumber in classTree:

            numberAsRR = RR(number)

        else:

            numberAsRR = number

    except AttributeError:

        return None

    # Zero special case.

    if numberAsRR == 0:

        return RR(-infinity)

    else:

        realField           = numberAsRR.parent()

        numberLog2          = numberAsRR.abs().log2()

        floorNumberLog2     = floor(numberLog2)

        ## Do not get caught by rounding of log2() both ways.

        ## When numberLog2 is an integer, compare numberAsRR

        #  with 2^numberLog2.

        if floorNumberLog2 == numberLog2:

            if numberAsRR.abs() < realField(2^floorNumberLog2):

                return floorNumberLog2 - 1

            else:

                return floorNumberLog2

        else:

            return floorNumberLog2

# End slz_compute_binade

#

def slz_compute_binade_bounds(number, emin, emax=sys.maxint):

    """

    For given "real number", compute the bounds of the binade it belongs to.

    NOTE::

        When number >= 2^(emax+1), we return the "fake" binade

        [2^(emax+1), +infinity]. Ditto for number <= -2^(emax+1)

        with interval [-infinity, -2^(emax+1)]. We want to distinguish

        this case from that of "really" invalid arguments.

    """

    # Check the parameters.

    # RealNumbers or RealNumber offspring only.

    # The exception construction is necessary since not all objects have

    # the mro() method. sage.rings.real_mpfr.RealNumber do.

    try:

        classTree = [number.__class__] + number.mro()

        if not sage.rings.real_mpfr.RealNumber in classTree:

            return None

    except AttributeError:

        return None

    # Non zero negative integers only for emin.

    if emin >= 0 or int(emin) != emin:

        return None

    # Non zero positive integers only for emax.

    if emax <= 0 or int(emax) != emax:

        return None

    precision = number.precision()

    RF  = RealField(precision)

    if number == 0:

        return (RF(0),RF(2^(emin)) - RF(2^(emin-precision)))

    # A more precise RealField is needed to avoid unwanted rounding effects

    # when computing number.log2().

    RRF = RealField(max(2048, 2 * precision))

    # number = 0 special case, the binade bounds are

    # [0, 2^emin - 2^(emin-precision)]

    # Begin general case

    l2 = RRF(number).abs().log2()

    # Another special one: beyond largest representable -> "Fake" binade.

    if l2 >= emax + 1:

        if number > 0:

            return (RF(2^(emax+1)), RF(+infinity) )

        else:

            return (RF(-infinity), -RF(2^(emax+1)))

    # Regular case cont'd.

    offset = int(l2)

    # number.abs() >= 1.

    if l2 >= 0:

        if number >= 0:

            lb = RF(2^offset)

            ub = RF(2^(offset + 1) - 2^(-precision+offset+1))

        else: #number < 0

            lb = -RF(2^(offset + 1) - 2^(-precision+offset+1))

            ub = -RF(2^offset)

    else: # log2 < 0, number.abs() < 1.

        if l2 < emin: # Denormal

           # print "Denormal:", l2

            if number >= 0:

                lb = RF(0)

                ub = RF(2^(emin)) - RF(2^(emin-precision))

            else: # number <= 0

                lb = - RF(2^(emin)) + RF(2^(emin-precision))

                ub = RF(0)

        elif l2 > emin: # Normal number other than +/-2^emin.

            if number >= 0:

                if int(l2) == l2:

                    lb = RF(2^(offset))

                    ub = RF(2^(offset+1)) - RF(2^(-precision+offset+1))

                else:

                    lb = RF(2^(offset-1))

                    ub = RF(2^(offset)) - RF(2^(-precision+offset))

            else: # number < 0

                if int(l2) == l2: # Binade limit.

                    lb = -RF(2^(offset+1) - 2^(-precision+offset+1))

                    ub = -RF(2^(offset))

                else:

                    lb = -RF(2^(offset) - 2^(-precision+offset))

                    ub = -RF(2^(offset-1))

        else: # l2== emin, number == +/-2^emin

            if number >= 0:

                lb = RF(2^(offset))

                ub = RF(2^(offset+1)) - RF(2^(-precision+offset+1))

            else: # number < 0

                lb = -RF(2^(offset+1) - 2^(-precision+offset+1))

                ub = -RF(2^(offset))

    return (lb, ub)

# End slz_compute_binade_bounds

#

def slz_compute_coppersmith_reduced_polynomials(inputPolynomial,

                                                 alpha,

                                                 N,

                                                 iBound,

                                                 tBound,

                                                 debug = False):

    """

    For a given set of arguments (see below), compute a list

    of "reduced polynomials" that could be used to compute roots

    of the inputPolynomial.

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

    if debug:

        polynomialsList = \

            spo_polynomial_to_polynomials_list_8(initialPolynomial,

                                                 alpha,

                                                 N,

                                                 iBound,

                                                 tBound,

                                                 20)

    else:

        polynomialsList = \

            spo_polynomial_to_polynomials_list_8(initialPolynomial,

                                                 alpha,

                                                 N,

                                                 iBound,

                                                 tBound,

                                                 0)

    #print "Polynomials list:", polynomialsList

    ## Building the proto matrix.

    knownMonomials = []

    protoMatrix    = []

    if debug:

        for poly in polynomialsList:

            spo_add_polynomial_coeffs_to_matrix_row(poly,

                                                    knownMonomials,

                                                    protoMatrix,

                                                    20)

    else:

        for poly in polynomialsList:

            spo_add_polynomial_coeffs_to_matrix_row(poly,

                                                    knownMonomials,

                                                    protoMatrix,

                                                    0)

    matrixToReduce = spo_proto_to_row_matrix(protoMatrix)

    #print matrixToReduce

    ## Reduction and checking.

    ## S.T. changed 'fp' to None as of Sage 6.6 complying to

    #  error message issued when previous code was used.

    #reducedMatrix = matrixToReduce.LLL(fp='fp')

    reducedMatrix = matrixToReduce.LLL(fp=None)

    isLLLReduced  = reducedMatrix.is_LLL_reduced()

    if not isLLLReduced:

        return None

    monomialsCount     = len(knownMonomials)

    monomialsCountSqrt = sqrt(monomialsCount)

    #print "Monomials count:", monomialsCount, monomialsCountSqrt.n()

    #print reducedMatrix

    ## Check the Coppersmith condition for each row and build the reduced

    #  polynomials.

    ccReducedPolynomialsList = []

    for row in reducedMatrix.rows():

        l2Norm = row.norm(2)

        if (l2Norm * monomialsCountSqrt) < nAtAlpha:

            #print (l2Norm * monomialsCountSqrt).n()

            #print l2Norm.n()

            ccReducedPolynomial = \

                slz_compute_reduced_polynomial(row,

                                               knownMonomials,

                                               iVariable,

                                               iBound,

                                               tVariable,

                                               tBound)

            if not ccReducedPolynomial is None:

                ccReducedPolynomialsList.append(ccReducedPolynomial)

        else:

            #print l2Norm.n() , ">", nAtAlpha

            pass

    if len(ccReducedPolynomialsList) < 2:

        print "Less than 2 Coppersmith condition compliant vectors."

        return ()

    #print ccReducedPolynomialsList

    return ccReducedPolynomialsList

# End slz_compute_coppersmith_reduced_polynomials

def slz_compute_coppersmith_reduced_polynomials_with_lattice_volume(inputPolynomial,

                                                                    alpha,

                                                                    N,

                                                                    iBound,

                                                                    tBound,

                                                                    debug = False):

    """

    For a given set of arguments (see below), compute a list

    of "reduced polynomials" that could be used to compute roots

    of the inputPolynomial.

    Print the volume of the initial basis as well.

    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

    if debug:

        print "N at alpha:", nAtAlpha

    ## 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,

##        spo_polynomial_to_polynomials_list_5(initialPolynomial,

    polynomialsList = \

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

    matrixToReduceTranspose = matrixToReduce.transpose()

    squareMatrix = matrixToReduce * matrixToReduceTranspose

    squareMatDet = det(squareMatrix)

    latticeVolume = sqrt(squareMatDet)

    print "Lattice volume:", latticeVolume.n()

    print "Lattice volume / N:", (latticeVolume/N).n()

    #print matrixToReduce

    ## Reduction and checking.

    ## S.T. changed 'fp' to None as of Sage 6.6 complying to

    #  error message issued when previous code was used.

    #reducedMatrix = matrixToReduce.LLL(fp='fp')

    reductionTimeStart = cputime()

    reducedMatrix = matrixToReduce.LLL(fp=None)

    reductionTime = cputime(reductionTimeStart)

    print "Reduction time:", reductionTime

    isLLLReduced  = reducedMatrix.is_LLL_reduced()

    if not isLLLReduced:

       return None

    #

    if debug:

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

        for row in reducedMatrix.rows():

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

        matrixFile.close()

    #

    monomialsCount     = len(knownMonomials)

    monomialsCountSqrt = sqrt(monomialsCount)

    #print "Monomials count:", monomialsCount, monomialsCountSqrt.n()

    #print reducedMatrix

    ## Check the Coppersmith condition for each row and build the reduced

    #  polynomials.

    ccVectorsCount = 0

    ccReducedPolynomialsList = []

    for row in reducedMatrix.rows():

        l2Norm = row.norm(2)

        if (l2Norm * monomialsCountSqrt) < nAtAlpha:

            #print (l2Norm * monomialsCountSqrt).n()

            #print l2Norm.n()

            ccVectorsCount +=1

            ccReducedPolynomial = \

                slz_compute_reduced_polynomial(row,

                                               knownMonomials,

                                               iVariable,

                                               iBound,

                                               tVariable,

                                               tBound)

            if not ccReducedPolynomial is None:

                ccReducedPolynomialsList.append(ccReducedPolynomial)

        else:

            #print l2Norm.n() , ">", nAtAlpha

            pass

    if debug:

        print ccVectorsCount, "out of ", len(ccReducedPolynomialsList),

        print "took Coppersmith text."

    if len(ccReducedPolynomialsList) < 2:

        print "Less than 2 Coppersmith condition compliant vectors."

        return ()

    if debug:

        print "Reduced and Coppersmith compliant polynomials list", ccReducedPolynomialsList

    return ccReducedPolynomialsList

# End slz_compute_coppersmith_reduced_polynomials_with_lattice volume

def slz_compute_initial_lattice_matrix(inputPolynomial,

                                       alpha,

                                       N,

                                       iBound,

                                       tBound,

                                       debug = False):

    """

    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)

    if debug:

        print "Initial basis polynomials"

        for poly in polynomialsList:

            print poly

    return matrixToReduce

# End slz_compute_initial_lattice_matrix.

def slz_compute_integer_polynomial_modular_roots(inputPolynomial,

                                                 alpha,

                                                 N,

                                                 iBound,

                                                 tBound):

    """

    For a given set of arguments (see below), compute the polynomial modular

    roots, if any.

    """

    # Arguments check.

    if iBound == 0 or tBound == 0:

        return set()

    # 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]

    ccReducedPolynomialsList = \

        slz_compute_coppersmith_reduced_polynomials (inputPolynomial,

                                                     alpha,

                                                     N,

                                                     iBound,

                                                     tBound)

    if len(ccReducedPolynomialsList) == 0:

        return set()

    ## Create the valid (poly1 and poly2 are algebraically independent)

    # resultant tuples (poly1, poly2, resultant(poly1, poly2)).

    # Try to mix and match all the polynomial pairs built from the

    # ccReducedPolynomialsList to obtain non zero resultants.

    resultantsInITuplesList = []

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

        for polyInnerIndex in xrange(polyOuterIndex+1,

                                     len(ccReducedPolynomialsList)):

            # Compute the resultant in resultants in the

            # first variable (is it the optimal choice?).

            resultantInI = \

               ccReducedPolynomialsList[polyOuterIndex].resultant(ccReducedPolynomialsList[polyInnerIndex],

                                                                  ccReducedPolynomialsList[0].parent(str(iVariable)))

            #print "Resultant", resultantInI

            # Test algebraic independence.

            if not resultantInI.is_zero():

                resultantsInITuplesList.append((ccReducedPolynomialsList[polyOuterIndex],

                                                 ccReducedPolynomialsList[polyInnerIndex],

                                                 resultantInI))

    # If no non zero resultant was found: we can't get no algebraically

    # independent polynomials pair. Give up!

    if len(resultantsInITuplesList) == 0:

        return set()

    #print resultantsInITuplesList

    # Compute the roots.

    Zi = ZZ[str(iVariable)]

    Zt = ZZ[str(tVariable)]

    polynomialRootsSet = set()

    # First, solve in the second variable since resultants are in the first

    # variable.

    for resultantInITuple in resultantsInITuplesList:

        tRootsList = Zt(resultantInITuple[2]).roots()

        # For each tRoot, compute the corresponding iRoots and check

        # them in the input polynomial.

        for tRoot in tRootsList:

            #print "tRoot:", tRoot

            # Roots returned by root() are (value, multiplicity) tuples.

            iRootsList = \

                Zi(resultantInITuple[0].subs({resultantInITuple[0].variables()[1]:tRoot[0]})).roots()

            print iRootsList

            # The iRootsList can be empty, hence the test.

            if len(iRootsList) != 0:

                for iRoot in iRootsList:

                    polyEvalModN = inputPolynomial(iRoot[0], tRoot[0]) / N

                    # polyEvalModN must be an integer.

                    if polyEvalModN == int(polyEvalModN):

                        polynomialRootsSet.add((iRoot[0],tRoot[0]))

    return polynomialRootsSet

# End slz_compute_integer_polynomial_modular_roots.

#

def slz_compute_integer_polynomial_modular_roots_2(inputPolynomial,

                                                 alpha,

                                                 N,

                                                 iBound,

                                                 tBound):

    """

    For a given set of arguments (see below), compute the polynomial modular

    roots, if any.

    This version differs in the way resultants are computed.

    """

    # Arguments check.

    if iBound == 0 or tBound == 0:

        return set()

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

    ccReducedPolynomialsList = \

        slz_compute_coppersmith_reduced_polynomials (inputPolynomial,

                                                     alpha,

                                                     N,

                                                     iBound,

                                                     tBound)

    if len(ccReducedPolynomialsList) == 0:

        return set()

    ## Create the valid (poly1 and poly2 are algebraically independent)

    # resultant tuples (poly1, poly2, resultant(poly1, poly2)).

    # Try to mix and match all the polynomial pairs built from the

    # ccReducedPolynomialsList to obtain non zero resultants.

    resultantsInTTuplesList = []

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

        for polyInnerIndex in xrange(polyOuterIndex+1,

                                     len(ccReducedPolynomialsList)):

            # Compute the resultant in resultants in the

            # first variable (is it the optimal choice?).

            resultantInT = \

               ccReducedPolynomialsList[polyOuterIndex].resultant(ccReducedPolynomialsList[polyInnerIndex],

                                                                  ccReducedPolynomialsList[0].parent(str(tVariable)))

            #print "Resultant", resultantInT

            # Test algebraic independence.

            if not resultantInT.is_zero():

                resultantsInTTuplesList.append((ccReducedPolynomialsList[polyOuterIndex],

                                                 ccReducedPolynomialsList[polyInnerIndex],

                                                 resultantInT))

    # If no non zero resultant was found: we can't get no algebraically

    # independent polynomials pair. Give up!

    if len(resultantsInTTuplesList) == 0:

        return set()

    #print resultantsInITuplesList

    # Compute the roots.

    Zi = ZZ[str(iVariable)]

    Zt = ZZ[str(tVariable)]

    polynomialRootsSet = set()

    # First, solve in the second variable since resultants are in the first

    # variable.

    for resultantInTTuple in resultantsInTTuplesList:

        iRootsList = Zi(resultantInTTuple[2]).roots()

        # For each iRoot, compute the corresponding tRoots and check

        # them in the input polynomial.

        for iRoot in iRootsList:

            #print "iRoot:", iRoot

            # Roots returned by root() are (value, multiplicity) tuples.

            tRootsList = \

                Zt(resultantInTTuple[0].subs({resultantInTTuple[0].variables()[0]:iRoot[0]})).roots()

            print tRootsList

            # The tRootsList can be empty, hence the test.

            if len(tRootsList) != 0:

                for tRoot in tRootsList:

                    polyEvalModN = inputPolynomial(iRoot[0],tRoot[0]) / N

                    # polyEvalModN must be an integer.

                    if polyEvalModN == int(polyEvalModN):

                        polynomialRootsSet.add((iRoot[0],tRoot[0]))

    return polynomialRootsSet

# End slz_compute_integer_polynomial_modular_roots_2.

#

def slz_compute_polynomial_and_interval(functionSo, degreeSo, lowerBoundSa,

                                        upperBoundSa, approxAccurSa,

                                        precSa=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 precSa is None:

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

        #print "Computed internal precision:", precSa

        if precSa < 192:

            precSa = 192

    sollyaPrecChanged = False

    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()

    if precSa > initialSollyaPrecSa:

        if precSa <= 2:

            print inspect.stack()[0][3], ": precision change <=2 requested."

        pobyso_set_prec_sa_so(precSa)

        sollyaPrecChanged = True

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

                pobyso_set_prec_so_so(initialSollyaPrecSo)

            sollya_lib_clear_obj(initialSollyaPrecSo)

            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)

    if sollyaPrecChanged:

        pobyso_set_prec_so_so(initialSollyaPrecSo)

    sollya_lib_clear_obj(initialSollyaPrecSo)

    return (polySo, currentRangeSo, intervalCenterSo, maxErrorSo)

# End slz_compute_polynomial_and_interval

def slz_compute_polynomial_and_interval_01(functionSo, degreeSo, lowerBoundSa,

                                        upperBoundSa, approxAccurSa,

                                        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:

        print inspect.stack()[0][3], ": lower bound is larger than upper bound. "

        return None

    ## Change Sollya precision, if requested.

    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()

    sollyaPrecChangedSa = False

    if sollyaPrecSa is None:

        sollyaPrecSa = initialSollyaPrecSa

    else:

        if sollyaPrecSa > initialSollyaPrecSa:

            if sollyaPrecSa <= 2:

                print inspect.stack()[0][3], ": precision change <= 2 requested."

            pobyso_set_prec_sa_so(sollyaPrecSa)

            sollyaPrecChangedSa = True

    ## Other initializations and data recovery.

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

                pobyso_set_prec_so_so(initialSollyaPrecSo)

            sollya_lib_clear_obj(initialSollyaPrecSo)

            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 inspect.stack()[0][3], "SollyaPrecSa:", sollyaPrecSa

        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)

    """

    # Naive rounding.

    (roundedPolySo, roundedPolyMaxErrSo) = \

    pobyso_polynomial_coefficients_progressive_round_so_so(polySo,

                                                           functionSo,

                                                           intervalCenterSo,

                                                           currentRangeSo,

                                                           itpSo,

                                                           ftpSo,

                                                           maxPrecSo,

                                                           approxAccurSo)

    """

    # Proved rounding.

    (roundedPolySo, roundedPolyMaxErrSo) = \

        pobyso_round_coefficients_progressive_so_so(polySo,

                                                    functionSo,

                                                    maxPrecSo,

                                                    currentRangeSo,

                                                    intervalCenterSo,

                                                    maxErrorSo,

                                                    approxAccurSo,

                                                    debug=False)

    #### Comment out the two next lines when polynomial rounding is activated.

    #roundedPolySo = sollya_lib_copy_obj(polySo)

    #roundedPolyMaxErrSo =  sollya_lib_copy_obj(maxErrorSo)

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

        pobyso_set_prec_so_so(initialSollyaPrecSo)

    sollya_lib_clear_obj(initialSollyaPrecSo)

    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,

                                        sollyaPrecSa=None, debug=True ):

    """

    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.

    Changes fom v 01:

        extra verbose.

    """

    print"In slz_compute_polynomial_and_interval..."

    ## Superficial argument check.

    if lowerBoundSa > upperBoundSa:

        return None

    ## Change Sollya precision, if requested.

    sollyaPrecChanged = False

    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()

    #print "Initial Sollya prec:", initialSollyaPrecSa, type(initialSollyaPrecSa)

    if sollyaPrecSa is None:

        sollyaPrecSa = initialSollyaPrecSa

    else:

        if sollyaPrecSa <= 2:

            print inspect.stack()[0][3], ": precision change <=2 requested."

        pobyso_set_prec_sa_so(sollyaPrecSa)

        sollyaPrecChanged = True

    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

    #pobyso_autoprint(functionSo)

    #pobyso_autoprint(degreeSo)

    #pobyso_autoprint(currentRangeSo)

    #pobyso_autoprint(absoluteErrorTypeSo)

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

    print "...after Taylor expansion."

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

                pobyso_set_prec_so_so(initialSollyaPrecSo)

            sollya_lib_clear_obj(initialSollyaPrecSo)

            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)