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..."

#

def slz_check_htr_value(function, htrValue, lowerBound, upperBound, precision, \

                        degree, targetHardnessToRound, alpha):

    """

    Check an Hard-to-round value.

    TODO::

    Full rewriting: this is hardly a draft.

    """

    polyApproxPrec = targetHardnessToRound + 1

    polyTargetHardnessToRound = targetHardnessToRound + 1

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

    RRR = htrValue.parent()

    #

    ## Compute the scaled function.

    fff = slz_compute_scaled_function(f, lowerBound, upperBound, precision)[0]

    print "Scaled function:", fff

    #

    ## Compute the scaling.

    boundsIntervalRifSa = RealIntervalField(precision)

    domainBoundsInterval = boundsIntervalRifSa(lowerBound, upperBound)

    scalingExpressions = \

        slz_interval_scaling_expression(domainBoundsInterval, i)

    #

    ## Get the polynomials, bounds, etc. for all the interval.

    resultListOfTuplesOfSo = \

        slz_get_intervals_and_polynomials(f, degree, lowerBound, upperBound, \

                                          precision, internalSollyaPrec,\

                                          2^-(polyApproxPrec))

    #

    ## We only want one interval.

    if len(resultListOfTuplesOfSo) > 1:

        print "Too many intervals! Aborting!"

        exit

    #

    ## Get the first tuple of Sollya objects as Sage objects.

    firstTupleSa = \

        slz_interval_and_polynomial_to_sage(resultListOfTuplesOfSo[0])

    pobyso_set_canonical_on()

    #

    print "Floatting point polynomial:", firstTupleSa[0]

    print "with coefficients precision:", firstTupleSa[0].base_ring().prec()

    #

    ## From a polynomial over a real ring, create a polynomial over the

    #  rationals ring.

    rationalPolynomial = \

        slz_float_poly_of_float_to_rat_poly_of_rat(firstTupleSa[0])

    print "Rational polynomial:", rationalPolynomial

    #

    ## Create a polynomial over the rationals that will take integer

    # variables instead of rational.

    rationalPolynomialOfIntegers = \

        slz_rat_poly_of_rat_to_rat_poly_of_int(rationalPolynomial, precision)

    print "Type:", type(rationalPolynomialOfIntegers)

    print "Rational polynomial of integers:", rationalPolynomialOfIntegers

    #

    ## Check the rational polynomial of integers variables.

    # (check against the scaled function).

    toIntegerFactor = 2^(precision-1)

    intervalCenterAsIntegerSa = int(firstTupleSa[3] * toIntegerFactor)

    print "Interval center as integer:", intervalCenterAsIntegerSa

    lowerBoundAsIntegerSa = int(firstTupleSa[2].endpoints()[0] * \

                                toIntegerFactor) - intervalCenterAsIntegerSa

    upperBoundAsIntegerSa = int(firstTupleSa[2].endpoints()[1] * \

                                toIntegerFactor) - intervalCenterAsIntegerSa

    print "Lower bound as integer:", lowerBoundAsIntegerSa

    print "Upper bound as integer:", upperBoundAsIntegerSa

    print "Image of the lower bound by the scaled function", \

            fff(firstTupleSa[2].endpoints()[0])

    print "Image of the lower bound by the approximation polynomial of ints:", \

            RRR(rationalPolynomialOfIntegers(lowerBoundAsIntegerSa))

    print "Image of the center by the scaled function", fff(firstTupleSa[3])

    print "Image of the center by the approximation polynomial of ints:", \

            RRR(rationalPolynomialOfIntegers(0))

    print "Image of the upper bound by the scaled function", \

            fff(firstTupleSa[2].endpoints()[1])

    print "Image of the upper bound by the approximation polynomial of ints:", \

            RRR(rationalPolynomialOfIntegers(upperBoundAsIntegerSa))

# End slz_check_htr_value.

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:

        return floor(numberAsRR.abs().log2())

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

    """

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

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

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

    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_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, approxPrecSa,

                                        sollyaPrecSa=None):

    """

    Under the assumptions listed for slz_get_intervals_and_polynomials, compute

    a polynomial that approximates the function on a an interval starting

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

    approximates with the expected precision.

    The interval upper bound is lowered until the expected approximation

    precision is reached.

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

    are returned.

    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 <= approxPrecSaS.

    """

    ## Superficial argument check.

    if lowerBoundSa > upperBoundSa:

        return None

    RRR = lowerBoundSa.parent()

    intervalShrinkConstFactorSa = RRR('0.8')

    absoluteErrorTypeSo = pobyso_absolute_so_so()

    currentRangeSo = pobyso_bounds_to_range_sa_so(lowerBoundSa, upperBoundSa)

    currentUpperBoundSa = upperBoundSa

    currentLowerBoundSa = lowerBoundSa

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

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

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

    (polySo, intervalCenterSo, maxErrorSo) = \

        pobyso_taylor_expansion_no_change_var_so_so(functionSo, degreeSo,

                                                    currentRangeSo,

                                                    absoluteErrorTypeSo)

    maxErrorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo)

    while maxErrorSa > approxPrecSa:

        print "++Approximation error:", maxErrorSa

        sollya_lib_clear_obj(polySo)

        sollya_lib_clear_obj(intervalCenterSo)

        sollya_lib_clear_obj(maxErrorSo)

        shrinkFactorSa = RRR('5')/(maxErrorSa/approxPrecSa).log2().abs()

        #shrinkFactorSa = 1.5/(maxErrorSa/approxPrecSa)

        #errorRatioSa = approxPrecSa/maxErrorSa

        #print "Error ratio: ", errorRatioSa

        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)

        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!"

            return (None, None, None, 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))

    sollya_lib_clear_obj(absoluteErrorTypeSo)

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

# End slz_compute_polynomial_and_interval

def slz_compute_reduced_polynomial(matrixRow,

                                    knownMonomials,

                                    var1,

                                    var1Bound,

                                    var2,

                                    var2Bound):

    """

    Compute a polynomial from a single reduced matrix row.

    This function was introduced in order to avoid the computation of the

    all the polynomials  from the full matrix (even those built from rows

    that do no verify the Coppersmith condition) as this may involves

    expensive operations over (large) integers.

    """

    ## Check arguments.

    if len(knownMonomials) == 0:

        return None

    # varNounds can be zero since 0^0 returns 1.

    if (var1Bound < 0) or (var2Bound < 0):

        return None

    ## Initialisations.

    polynomialRing = knownMonomials[0].parent()

    currentPolynomial = polynomialRing(0)

    # TODO: use zip instead of indices.

    for colIndex in xrange(0, len(knownMonomials)):

        currentCoefficient = matrixRow[colIndex]

        if currentCoefficient != 0:

            #print "Current coefficient:", currentCoefficient

            currentMonomial = knownMonomials[colIndex]

            #print "Monomial as multivariate polynomial:", \

            #currentMonomial, type(currentMonomial)

            degreeInVar1 = currentMonomial.degree(var1)

            #print "Degree in var1", var1, ":", degreeInVar1

            degreeInVar2 = currentMonomial.degree(var2)

            #print "Degree in var2", var2, ":", degreeInVar2

            if degreeInVar1 > 0:

                currentCoefficient = \

                    currentCoefficient / (var1Bound^degreeInVar1)

                #print "varBound1 in degree:", var1Bound^degreeInVar1

                #print "Current coefficient(1)", currentCoefficient

            if degreeInVar2 > 0:

                currentCoefficient = \

                    currentCoefficient / (var2Bound^degreeInVar2)

                #print "Current coefficient(2)", currentCoefficient

            #print "Current reduced monomial:", (currentCoefficient * \

            #                                    currentMonomial)

            currentPolynomial += (currentCoefficient * currentMonomial)

            #print "Current polynomial:", currentPolynomial

        # End if

    # End for colIndex.

    #print "Type of the current polynomial:", type(currentPolynomial)

    return(currentPolynomial)

# End slz_compute_reduced_polynomial

#

def slz_compute_reduced_polynomials(reducedMatrix,

                                        knownMonomials,

                                        var1,

                                        var1Bound,

                                        var2,

                                        var2Bound):

    """

    Legacy function, use slz_compute_reduced_polynomials_list

    """

    return(slz_compute_reduced_polynomials_list(reducedMatrix,

                                                knownMonomials,

                                                var1,

                                                var1Bound,

                                                var2,

                                                var2Bound)

    )

def slz_compute_reduced_polynomials_list(reducedMatrix,

                                         knownMonomials,

                                         var1,

                                         var1Bound,

                                         var2,

                                         var2Bound):

    """

    From a reduced matrix, holding the coefficients, from a monomials list,

    from the bounds of each variable, compute the corresponding polynomials

    scaled back by dividing by the "right" powers of the variables bounds.

    The elements in knownMonomials must be of the "right" polynomial type.

    They set the polynomial type of the output polynomials in list.

    @param reducedMatrix:  the reduced matrix as output from LLL;

    @param kwnonMonomials: the ordered list of the monomials used to

                           build the polynomials;

    @param var1:           the first variable (of the "right" type);

    @param var1Bound:      the first variable bound;

    @param var2:           the second variable (of the "right" type);

    @param var2Bound:      the second variable bound.

    @return: a list of polynomials obtained with the reduced coefficients

             and scaled down with the bounds

    """

    # TODO: check input arguments.

    reducedPolynomials = []

    #print "type var1:", type(var1), " - type var2:", type(var2)

    for matrixRow in reducedMatrix.rows():

        currentPolynomial = 0

        for colIndex in xrange(0, len(knownMonomials)):

            currentCoefficient = matrixRow[colIndex]

            if currentCoefficient != 0:

                #print "Current coefficient:", currentCoefficient

                currentMonomial = knownMonomials[colIndex]

                parentRing = currentMonomial.parent()

                #print "Monomial as multivariate polynomial:", \

                #currentMonomial, type(currentMonomial)

                degreeInVar1 = currentMonomial.degree(parentRing(var1))

                #print "Degree in var", var1, ":", degreeInVar1

                degreeInVar2 = currentMonomial.degree(parentRing(var2))

                #print "Degree in var", var2, ":", degreeInVar2

                if degreeInVar1 > 0:

                    currentCoefficient /= var1Bound^degreeInVar1

                    #print "varBound1 in degree:", var1Bound^degreeInVar1

                    #print "Current coefficient(1)", currentCoefficient

                if degreeInVar2 > 0:

                    currentCoefficient /= var2Bound^degreeInVar2

                    #print "Current coefficient(2)", currentCoefficient

                #print "Current reduced monomial:", (currentCoefficient * \

                #                                    currentMonomial)

                currentPolynomial += (currentCoefficient * currentMonomial)

                #if degreeInVar1 == 0 and degreeInVar2 == 0:

                    #print "!!!! constant term !!!!"

                #print "Current polynomial:", currentPolynomial

            # End if

        # End for colIndex.

        #print "Type of the current polynomial:", type(currentPolynomial)

        reducedPolynomials.append(currentPolynomial)

    return reducedPolynomials

# End slz_compute_reduced_polynomials.

def slz_compute_scaled_function(functionSa,

                                lowerBoundSa,

                                upperBoundSa,

                                floatingPointPrecSa,

                                debug=False):

    """

    From a function, compute the scaled function whose domain

    is included in [1, 2) and whose image is also included in [1,2).

    Return a tuple:

    [0]: the scaled function

    [1]: the scaled domain lower bound

    [2]: the scaled domain upper bound

    [3]: the scaled image lower bound

    [4]: the scaled image upper bound

    """

    x = functionSa.variables()[0]

    # Reassert f as a function (an not a mere expression).

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

    boundsIntervalRifSa = RealIntervalField(floatingPointPrecSa)

    domainBoundsIntervalSa = boundsIntervalRifSa(lowerBoundSa, upperBoundSa)

    (invDomainScalingExpressionSa, domainScalingExpressionSa) = \

        slz_interval_scaling_expression(domainBoundsIntervalSa, x)

    if debug:

        print "domainScalingExpression for argument :", \

        invDomainScalingExpressionSa

        print "f: ", f

    ff = f.subs({x : domainScalingExpressionSa})

    #ff = f.subs_expr(x==domainScalingExpressionSa)

    domainScalingFunction(x) = invDomainScalingExpressionSa

    scaledLowerBoundSa = \

        domainScalingFunction(lowerBoundSa).n(prec=floatingPointPrecSa)

    scaledUpperBoundSa = \

        domainScalingFunction(upperBoundSa).n(prec=floatingPointPrecSa)

    if debug:

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

              scaledUpperBoundSa

    #

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

    flbSa = ff(scaledLowerBoundSa).n(prec=floatingPointPrecSa)

    fubSa = ff(scaledUpperBoundSa).n(prec=floatingPointPrecSa)

    if flbSa <= fubSa: # Increasing

        imageBinadeBottomSa = floor(flbSa.log2())

    else: # Decreasing

        imageBinadeBottomSa = floor(fubSa.log2())

    if debug:

        print 'ff:', ff, '- Image:', flbSa, fubSa, imageBinadeBottomSa

    imageBoundsIntervalSa = boundsIntervalRifSa(flbSa, fubSa)

    (invImageScalingExpressionSa,imageScalingExpressionSa) = \

        slz_interval_scaling_expression(imageBoundsIntervalSa, x)

    if debug:

        print "imageScalingExpression for argument :", \

              invImageScalingExpressionSa

    iis = invImageScalingExpressionSa.function(x)

    fff = iis.subs({x:ff})

    if debug:

        print "fff:", fff,

        print " - Image:", fff(scaledLowerBoundSa), fff(scaledUpperBoundSa)

    return([fff, scaledLowerBoundSa, scaledUpperBoundSa, \

            fff(scaledLowerBoundSa), fff(scaledUpperBoundSa)])

# End slz_compute_scaled_function

def slz_float_poly_of_float_to_rat_poly_of_rat(polyOfFloat):

    # Create a polynomial over the rationals.

    polynomialRing = QQ[str(polyOfFloat.variables()[0])]

    return(polynomialRing(polyOfFloat))

# End slz_float_poly_of_float_to_rat_poly_of_rat.

def slz_get_intervals_and_polynomials(functionSa, degreeSa,

                                      lowerBoundSa,

                                      upperBoundSa, floatingPointPrecSa,

                                      internalSollyaPrecSa, approxPrecSa):

    """

    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;

    - the requested approximation error

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

    It return a list of tuples, each made of:

    - a first polynomial (without the changed variable f(x) = p(x-x0));

    - a second polynomial (with a changed variable f(x) = q(x))

    - the approximation interval;

    - the center, x0, of the interval;

    - the corresponding approximation error.

    TODO: fix endless looping for some parameters sets.

    """

    resultArray = []

    # Set Sollya to the necessary internal precision.

    precChangedSa = False

    currentSollyaPrecSo = pobyso_get_prec_so()

    currentSollyaPrecSa = pobyso_constant_from_int_so_sa(currentSollyaPrecSo)

    if internalSollyaPrecSa > currentSollyaPrecSa:

        pobyso_set_prec_sa_so(internalSollyaPrecSa)

        precChangedSa = True

    #

    x = functionSa.variables()[0] # Actual variable name can be anything.

    # Scaled function: [1=,2] -> [1,2].

    (fff, scaledLowerBoundSa, scaledUpperBoundSa,                             \

            scaledLowerBoundImageSa, scaledUpperBoundImageSa) =               \

                slz_compute_scaled_function(functionSa,                       \

                                            lowerBoundSa,                     \

                                            upperBoundSa,                     \

                                            floatingPointPrecSa)

    # In case bounds were in the negative real one may need to

    # switch scaled bounds.

    if scaledLowerBoundSa > scaledUpperBoundSa:

        scaledLowerBoundSa, scaledUpperBoundSa = \

            scaledUpperBoundSa, scaledLowerBoundSa

        #print "Switching!"

    print "Approximation precision: ", RR(approxPrecSa)

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

    functionSo = \

      pobyso_parse_string_sa_so(fff._assume_str().replace('_SAGE_VAR_', ''))

    degreeSo = pobyso_constant_from_int_sa_so(degreeSa)

    scaledBoundsSo = pobyso_bounds_to_range_sa_so(scaledLowerBoundSa,

                                                  scaledUpperBoundSa)

    # Compute the first Taylor expansion.

    print "Computing first Taylor expansion..."

    (polySo, boundsSo, intervalCenterSo, maxErrorSo) = \

        slz_compute_polynomial_and_interval(functionSo, degreeSo,

                                        scaledLowerBoundSa, scaledUpperBoundSa,

                                        approxPrecSa, internalSollyaPrecSa)

    print "...done."

    if polySo is None:

        print "slz_get_intervals_and_polynomials: Aborting and returning None!"

        if precChangedSa:

            pobyso_set_prec_so_so(currentSollyaPrecSo)

            sollya_lib_clear_obj(currentSollyaPrecSo)

        sollya_lib_clear_obj(functionSo)

        sollya_lib_clear_obj(degreeSo)

        sollya_lib_clear_obj(scaledBoundsSo)

        return None

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

                                              upperBoundSa.parent().precision()))

    boundsSa = pobyso_range_to_interval_so_sa(boundsSo, realIntervalField)

    errorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo)

    print "First approximation error:", errorSa.n(digits=10)

    # If the error and interval are OK a the first try, just return.

    if boundsSa.endpoints()[1] >= scaledUpperBoundSa:

        # Change variable stuff in Sollya x -> x0-x.

        changeVarExpressionSo = sollya_lib_build_function_sub( \

                              sollya_lib_build_function_free_variable(), \

                              sollya_lib_copy_obj(intervalCenterSo))

        polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpressionSo)

        sollya_lib_clear_obj(changeVarExpressionSo)

        resultArray.append((polySo, polyVarChangedSo, boundsSo, \

                            intervalCenterSo, maxErrorSo))

        if internalSollyaPrecSa != currentSollyaPrecSa:

            pobyso_set_prec_sa_so(currentSollyaPrecSa)

        sollya_lib_clear_obj(currentSollyaPrecSo)

        sollya_lib_clear_obj(functionSo)

        sollya_lib_clear_obj(degreeSo)

        sollya_lib_clear_obj(scaledBoundsSo)

        #print "Approximation error:", errorSa

        return resultArray

    # The returned interval upper bound does not reach the requested upper

    # upper bound: compute the next upper bound.

    # The following ratio is always >= 1

    currentErrorRatio = approxPrecSa / errorSa

    # Starting point for the next upper bound.

    currentScaledUpperBoundSa = boundsSa.endpoints()[1]

    boundsWidthSa = boundsSa.endpoints()[1] - boundsSa.endpoints()[0]

    # Compute the increment.

    if currentErrorRatio > RR('1000'): # ]1.5, infinity[

        currentScaledUpperBoundSa += \

                        currentErrorRatio * boundsWidthSa * 2

    else:  # [1, 1.5]

        currentScaledUpperBoundSa += \

                        (RR('1.0') + currentErrorRatio.log() / 500) * boundsWidthSa

    # Take into account the original interval upper bound.

    if currentScaledUpperBoundSa > scaledUpperBoundSa:

        currentScaledUpperBoundSa = scaledUpperBoundSa

    # Compute the other expansions.

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

        currentScaledLowerBoundSa = boundsSa.endpoints()[1]

        (polySo, boundsSo, intervalCenterSo, maxErrorSo) = \

            slz_compute_polynomial_and_interval(functionSo, degreeSo,

                                            currentScaledLowerBoundSa,

                                            currentScaledUpperBoundSa,

                                            approxPrecSa,

                                            internalSollyaPrecSa)

        errorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo)

        if  errorSa < approxPrecSa:

            # Change variable stuff

            #print "Approximation error:", errorSa

            changeVarExpressionSo = sollya_lib_build_function_sub(

                                    sollya_lib_build_function_free_variable(),

                                    sollya_lib_copy_obj(intervalCenterSo))

            polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpressionSo)

            sollya_lib_clear_obj(changeVarExpressionSo)

            resultArray.append((polySo, polyVarChangedSo, boundsSo, \

                                intervalCenterSo, maxErrorSo))

        boundsSa = pobyso_range_to_interval_so_sa(boundsSo, realIntervalField)

        # Compute the next upper bound.

        # The following ratio is always >= 1

        currentErrorRatio = approxPrecSa / errorSa

        # Starting point for the next upper bound.

        currentScaledUpperBoundSa = boundsSa.endpoints()[1]

        boundsWidthSa = boundsSa.endpoints()[1] - boundsSa.endpoints()[0]

        # Compute the increment.

        if currentErrorRatio > RR('1000'): # ]1.5, infinity[

            currentScaledUpperBoundSa += \

                            currentErrorRatio * boundsWidthSa * 2

        else:  # [1, 1.5]

            currentScaledUpperBoundSa += \

                            (RR('1.0') + currentErrorRatio.log()/500) * boundsWidthSa

        #print "currentErrorRatio:", currentErrorRatio

        #print "currentScaledUpperBoundSa", currentScaledUpperBoundSa

        # Test for insufficient precision.

        if currentScaledUpperBoundSa == scaledLowerBoundSa:

            print "Can't shrink the interval anymore!"

            print "You should consider increasing the Sollya internal precision"

            print "or the polynomial degree."

            print "Giving up!"

            if internalSollyaPrecSa != currentSollyaPrecSa:

                pobyso_set_prec_sa_so(currentSollyaPrecSa)

            sollya_lib_clear_obj(currentSollyaPrecSo)

            sollya_lib_clear_obj(functionSo)

            sollya_lib_clear_obj(degreeSo)

            sollya_lib_clear_obj(scaledBoundsSo)

            return None

        if currentScaledUpperBoundSa > scaledUpperBoundSa:

            currentScaledUpperBoundSa = scaledUpperBoundSa

    if internalSollyaPrecSa > currentSollyaPrecSa:

        pobyso_set_prec_so_so(currentSollyaPrecSo)

    sollya_lib_clear_obj(currentSollyaPrecSo)

    sollya_lib_clear_obj(functionSo)

    sollya_lib_clear_obj(degreeSo)

    sollya_lib_clear_obj(scaledBoundsSo)

    return(resultArray)

# End slz_get_intervals_and_polynomials

def slz_interval_scaling_expression(boundsInterval, expVar):

    """

    Compute the scaling expression to map an interval that spans at most

    a single binade into [1, 2) and the inverse expression as well.

    If the interval spans more than one binade, result is wrong since

    scaling is only based on the lower bound.

    Not very sure that the transformation makes sense for negative numbers.

    """

    # The "one of the bounds is 0" special case. Here we consider

    # the interval as the subnormals binade.

    if boundsInterval.endpoints()[0] == 0 or boundsInterval.endpoints()[1] == 0:

        # The upper bound is (or should be) positive.

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

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

                return None

            binade = slz_compute_binade(boundsInterval.endpoints()[1])

            l2     = boundsInterval.endpoints()[1].abs().log2()

            # The upper bound is a power of two

            if binade == l2:

                scalingCoeff    = 2^(-binade)

                invScalingCoeff = 2^(binade)

                scalingOffset   = 1

                return((scalingCoeff * expVar + scalingOffset),\

                       ((expVar - scalingOffset) * invScalingCoeff))

            else:

                scalingCoeff    = 2^(-binade-1)

                invScalingCoeff = 2^(binade+1)

                scalingOffset   = 1

                return((scalingCoeff * expVar + scalingOffset),\

                    ((expVar - scalingOffset) * invScalingCoeff))

        # The lower bound is (or should be) negative.

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

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

                return None

            binade = slz_compute_binade(boundsInterval.endpoints()[0])

            l2     = boundsInterval.endpoints()[0].abs().log2()

            # The upper bound is a power of two

            if binade == l2:

                scalingCoeff    = -2^(-binade)

                invScalingCoeff = -2^(binade)

                scalingOffset   = 1

                return((scalingCoeff * expVar + scalingOffset),\

                    ((expVar - scalingOffset) * invScalingCoeff))

            else:

                scalingCoeff    = -2^(-binade-1)

                invScalingCoeff = -2^(binade+1)

                scalingOffset   = 1

                return((scalingCoeff * expVar + scalingOffset),\

                   ((expVar - scalingOffset) * invScalingCoeff))

    # General case

    lbBinade = slz_compute_binade(boundsInterval.endpoints()[0])

    ubBinade = slz_compute_binade(boundsInterval.endpoints()[1])

    # We allow for a single exception in binade spanning is when the

    # "outward bound" is a power of 2 and its binade is that of the

    # "inner bound" + 1.

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

        ubL2 = boundsInterval.endpoints()[1].abs().log2()

        if lbBinade != ubBinade:

            print "Different binades."

            if ubL2 != ubBinade:

                print "Not a power of 2."

                return None

            elif abs(ubBinade - lbBinade) > 1:

                print "Too large span:", abs(ubBinade - lbBinade)

                return None

    else:

        lbL2 = boundsInterval.endpoints()[0].abs().log2()

        if lbBinade != ubBinade:

            print "Different binades."

            if lbL2 != lbBinade:

                print "Not a power of 2."

                return None

            elif abs(ubBinade - lbBinade) > 1:

                print "Too large span:", abs(ubBinade - lbBinade)

                return None

    #print "Lower bound binade:", binade

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

        return((2^(-lbBinade) * expVar),(2^(lbBinade) * expVar))