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

print "sageSLZ loading..."

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.

    """

    RRR = lowerBoundSa.parent()

    intervalShrinkConstFactorSa = RRR('0.5')

    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:

        sollya_lib_clear_obj(maxErrorSo)

        sollya_lib_clear_obj(polySo)

        sollya_lib_clear_obj(intervalCenterSo)

        shrinkFactorSa = RRR('5.0')/(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

        currentUpperBoundSa = currentLowerBoundSa + \

                                  (currentUpperBoundSa - currentLowerBoundSa) * \

                                   actualShrinkFactorSa

        #print "Current upper bound:", currentUpperBoundSa

        sollya_lib_clear_obj(currentRangeSo)

        sollya_lib_clear_obj(polySo)

        if currentUpperBoundSa <= currentLowerBoundSa:

            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)

        (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_scaled_function(functionSa, \

                                      lowerBoundSa, \

                                      upperBoundSa, \

                                      floatingPointPrecSa):

    """

    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)

    (domainScalingExpressionSa, invDomainScalingExpressionSa) = \

        slz_interval_scaling_expression(domainBoundsIntervalSa, x)

    print "domainScalingExpression for argument :", domainScalingExpressionSa

    print "f: ", f

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

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

    domainScalingFunction(x) = invDomainScalingExpressionSa

    scaledLowerBoundSa = domainScalingFunction(lowerBoundSa).n()

    scaledUpperBoundSa = domainScalingFunction(upperBoundSa).n()

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

    #

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

    flbSa = f(lowerBoundSa).n()

    fubSa = f(upperBoundSa).n()

    if flbSa <= fubSa: # Increasing

        imageBinadeBottomSa = floor(flbSa.log2())

    else: # Decreasing

        imageBinadeBottomSa = floor(fubSa.log2())

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

    imageBoundsIntervalSa = boundsIntervalRifSa(flbSa, fubSa)

    (imageScalingExpressionSa, invImageScalingExpressionSa) = \

        slz_interval_scaling_expression(imageBoundsIntervalSa, x)

    iis = invImageScalingExpressionSa.function(x)

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

    print "fff:", fff,

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

    return([fff, scaledLowerBoundSa, scaledUpperBoundSa, \

            fff(scaledLowerBoundSa), fff(scaledUpperBoundSa)])

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.

    """

    currentSollyaPrecSo = pobyso_get_prec_so()

    currentSollyaPrecSa = pobyso_constant_from_int_so_sa(currentSollyaPrecSo)

    if internalSollyaPrecSa > currentSollyaPrecSa:

        pobyso_set_prec_sa_so(internalSollyaPrecSa)

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

    (fff, scaledLowerBoundSa, scaledUpperBoundSa, \

            scaledLowerBoundImageSa, scaledUpperBoundImageSa) = \

                slz_compute_scaled_function(functionSa, \

                                            lowerBoundSa, \

                                            upperBoundSa, \

                                            floatingPointPrecSa)

    #

    resultArray = []

    #

    print "Approximation precision: ", RR(approxPrecSa)

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

    functionSo = pobyso_parse_string_sa_so(fff._assume_str())

    degreeSo = pobyso_constant_from_int_sa_so(degreeSa)

    scaledBoundsSo = pobyso_bounds_to_range_sa_so(scaledLowerBoundSa, 

                                                  scaledUpperBoundSa)

    # Compute the first Taylor expansion.

    (polySo, boundsSo, intervalCenterSo, maxErrorSo) = \

        slz_compute_polynomial_and_interval(functionSo, degreeSo,

                                        scaledLowerBoundSa, scaledUpperBoundSa,

                                        approxPrecSa, internalSollyaPrecSa)

    if polySo is None:

        print "Aborting"

        return None

    # Change variable stuff

    changeVarExpressionSo = sollya_lib_build_function_sub(

                              sollya_lib_build_function_free_variable(), 

                              sollya_lib_copy_obj(intervalCenterSo))

    polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpressionSo) 

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

                         maxErrorSo))

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

                                              upperBoundSa.parent().precision()))

    boundsSa = pobyso_range_to_interval_so_sa(boundsSo, realIntervalField)

    # Compute the next upper bound.

    # If the error of approximation is more than half of the target,

    # use the same interval.

    # If it is less, increase it a bit.

    errorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo)

    currentErrorRatio = approxPrecSa / errorSa

    currentScaledUpperBoundSa = boundsSa.endpoints()[1]

    if currentErrorRatio  < 2 :

        currentScaledUpperBoundSa += \

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

    else:

        currentScaledUpperBoundSa += \

            (boundsSa.endpoints()[1] - boundsSa.endpoints()[0]) \

                         * currentErrorRatio.log2() * 2

    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)

        # Change variable stuff

        changeVarExpressionSo = sollya_lib_build_function_sub(

                                  sollya_lib_build_function_free_variable(), 

                                  sollya_lib_copy_obj(intervalCenterSo))

        polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpressionSo) 

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

                            intervalCenterSo, maxErrorSo))

        boundsSa = pobyso_range_to_interval_so_sa(boundsSo, realIntervalField)

        # Compute the next upper bound.

        # If the error of approximation is more than half of the target,

        # use the same interval.

        # If it is less, increase it a bit.

        errorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo)

        currentErrorRatio = approxPrecSa / errorSa

        if currentErrorRatio  < RR('1.5') :

            currentScaledUpperBoundSa = \

                            boundsSa.endpoints()[1] + \

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

        elif currentErrorRatio < 2:

            currentScaledUpperBoundSa = \

                            boundsSa.endpoints()[1] + \

                            (boundsSa.endpoints()[1] - boundsSa.endpoints()[0]) \

                             * currentErrorRatio.log2()

        else:

            currentScaledUpperBoundSa = \

                            boundsSa.endpoints()[1] + \

                            (boundsSa.endpoints()[1] - boundsSa.endpoints()[0]) \

                             * currentErrorRatio.log2() * 2

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

            sollya_lib_clear_obj(functionSo)

            sollya_lib_clear_obj(degreeSo)

            sollya_lib_clear_obj(scaledBoundsSo)

            return None

        if currentScaledUpperBoundSa > scaledUpperBoundSa:

            currentScaledUpperBoundSa = scaledUpperBoundSa

    sollya_lib_clear_obj(functionSo)

    sollya_lib_clear_obj(degreeSo)

    sollya_lib_clear_obj(scaledBoundsSo)

    if internalSollyaPrecSa > currentSollyaPrecSa:

        pobyso_set_prec_so_so(currentSollyaPrecSo)

    return(resultArray)

# End slz_get_intervals_and_polynomials

def slz_interval_scaling_expression(boundsInterval, expVar):

    """

    Compute the scaling expression to map an interval that span only

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

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

    """

    # The scaling offset is only used for negative numbers.

    if abs(boundsInterval.endpoints()[0]) < 1:

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

            scalingCoeff = 2^floor(boundsInterval.endpoints()[0].log2())

            invScalingCoeff = 1/scalingCoeff

            return((scalingCoeff * expVar, 

                    invScalingCoeff * expVar))

        else:

            scalingCoeff = \

                2^(floor((-boundsInterval.endpoints()[0]).log2()) - 1)

            scalingOffset = -3 * scalingCoeff

            return((scalingCoeff * expVar + scalingOffset,

                    1/scalingCoeff * expVar + 3))

    else: 

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

            scalingCoeff = 2^floor(boundsInterval.endpoints()[0].log2())

            scalingOffset = 0

            return((scalingCoeff * expVar, 

                    1/scalingCoeff * expVar))

        else:

            scalingCoeff = \

                2^(floor((-boundsInterval.endpoints()[1]).log2()))

            scalingOffset =  -3 * scalingCoeff

            #scalingOffset = 0

            return((scalingCoeff * expVar + scalingOffset,

                    1/scalingCoeff * expVar + 3))

def slz_interval_and_polynomial_to_sage(polyRangeCenterErrorSo):

    """

    Compute the Sage version of the Taylor polynomial and it's

    companion data (interval, center...)

    The input parameter is a five elements tuple:

    - [0]: the polyomial (without variable change), as polynomial over a

           real ring;

    - [1]: the polyomial (with variable change done in Sollya), as polynomial

           over a real ring;

    - [2]: the interval (as Sollya range);

    - [3]: the interval center;

    - [4]: the approximation error. 

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

    input elements converted into their Sollya counterpart.

    """

    polynomialSa = pobyso_get_poly_so_sa(polyRangeCenterErrorSo[0])

    polynomialChangedVarSa = pobyso_get_poly_so_sa(polyRangeCenterErrorSo[1])

    intervalSa = \

            pobyso_get_interval_from_range_so_sa(polyRangeCenterErrorSo[2])

    centerSa = \

            pobyso_get_constant_as_rn_with_rf_so_sa(polyRangeCenterErrorSo[3])

    errorSa = \

            pobyso_get_constant_as_rn_with_rf_so_sa(polyRangeCenterErrorSo[4])

    return((polynomialSa, polynomialChangedVarSa, intervalSa, centerSa, errorSa))

    # End slz_interval_and_polynomial_to_sage

def slz_rat_poly_of_int_to_poly_for_coppersmith(ratPolyOfInt, 

                                                precision,

                                                targetHardnessToRound,

                                                variable1,

                                                variable2):

    """

    Creates a new polynomial with integer coefficients for use with the

    Coppersmith method.

    A the same time it computes :

    - 2^K (N);

    - 2^k

    - lcm

    """

    # Create a new integer polynomial ring.

    IP = PolynomialRing(ZZ, name=str(variable1) + "," + str(variable2))

    # Coefficients are issued in the increasing power order.

    ratPolyCoefficients = ratPolyOfInt.coefficients()

    # Build the list of number we compute the lcmm of.

    coefficientDenominators = sro_denominators(ratPolyCoefficients)

    coefficientDenominators.append(2^precision)

    coefficientDenominators.append(2^(targetHardnessToRound + 1))

    # Compute the lcm

    leastCommonMultiple = lcm(coefficientDenominators)

    # Compute the expression corresponding to the new polynomial

    coefficientNumerators =  sro_numerators(ratPolyCoefficients)

    print coefficientNumerators

    polynomialExpression = 0

    power = 0

    # Iterate over two lists at the same time, stop when the shorter is

    # exhausted.

    for numerator, denominator in \

                            zip(coefficientNumerators, coefficientDenominators):

        multiplicator = leastCommonMultiple / denominator 

        newCoefficient = numerator * multiplicator

        polynomialExpression += newCoefficient * variable1^power

        power +=1

    polynomialExpression += - variable2

    return (IP(polynomialExpression),

            leastCommonMultiple / 2^precision, # 2^K or N.

            leastCommonMultiple / 2 ^(targetHardnessToRound + 1), # tBound

            leastCommonMultiple)

# End slz_ratPoly_of_int_to_poly_for_coppersmith

def slz_rat_poly_of_rat_to_rat_poly_of_int(ratPolyOfRat, 

                                          precision):

    """

    Makes a variable substitution into the input polynomial so that the output

    polynomial can take integer arguments.

    All variables of the input polynomial "have precision p". That is to say

    that they are rationals with denominator == 2^precision: x = y/2^precision

    We "incorporate" these denominators into the coefficients with, 

    respectively, the "right" power.

    """

    polynomialField = ratPolyOfRat.parent()

    polynomialVariable = rationalPolynomial.variables()[0]

    print "The polynomial field is:", polynomialField

    return \

        polynomialField(rationalPolynomial.subs({polynomialVariable : \

                                   polynomialVariable/2^(precision-1)}))

    # Return a tuple:

    # - the bivariate integer polynomial in (i,j);

    # - 2^K

# End slz_rat_poly_of_rat_to_rat_poly_of_int

print "\t...sageSLZ loaded"