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

#! /opt/sage/sage

from scipy.constants.codata import precision

def initialize_env():

    #Load all necessary modules.

    if not 'mpfi' in sage.misc.cython.standard_libs:

        sage.misc.cython.standard_libs.append('mpfi')

    load("/home/storres/recherche/arithmetique/pobysoPythonSage/src/sollya_lib.sage")

    load("/home/storres/recherche/arithmetique/pobysoPythonSage/src/sageMpfr.spyx")

    load("/home/storres/recherche/arithmetique/pobysoPythonSage/src/pobyso.py")

    load("/home/storres/recherche/arithmetique/pobysoPythonSage/src/sageSLZ/sageSLZ.sage")

    load("/home/storres/recherche/arithmetique/pobysoPythonSage/src/sageSLZ/sageNumericalOperations.sage")

    load("/home/storres/recherche/arithmetique/pobysoPythonSage/src/sageSLZ/sageRationalOperations.sage")

    # Matrix operations are loaded by polynomial operations.

    load("/home/storres/recherche/arithmetique/pobysoPythonSage/src/sageSLZ/sagePolynomialOperations.sage")

def run_SLZ_v01(inputFunction,

                inputLowerBound,

                inputUpperBound,

                alpha,

                degree,

                precision,

                emin,

                emax,

                targetHardnessToRound,

                debug = False):

    if debug:

        print "Function                :", inputFunction

        print "Lower bound             :", inputLowerBound

        print "Upper bounds            :", inputUpperBound

        print "Alpha                   :", alpha

        print "Degree                  :", degree

        print "Precision               :", precision

        print "Emin                    :", emin

        print "Emax                    :", emax

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

        print

    ## Structures.

    RRR         = RealField(precision)

    RRIF        = RealIntervalField(precision)

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

    lowerBound = RRR(inputLowerBound)

    upperBound = RRR(inputUpperBound)

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

    print inputFunction(1).n()

    (lb,ub) = slz_fix_bounds_for_binades(lowerBound, upperBound, inputFunction)

    if lb != lowerBound or ub != upperBound:

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

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

        lowerBound, upperBound, "Images:", \

        inputFunction(lowerBound), inputFunction(upperBound)

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

    #

    ## Progam initialization

    ### Approximation polynomial accuracy and hardness to round.

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

    polyTargetHardnessToRound = targetHardnessToRound + 1

    ### Significand to integer conversion ratio.

    toIntegerFactor           = 2^(precision-1)

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

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

    i=var('i')

    t=var('t')

    inputFunctionVariable = inputFunction.variables()[0]

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

    #  Polynomial Rings over the integers, for root finding.

    Zi = ZZ[i]

    Zt = ZZ[t]

    Zit = ZZ[i,t]

    ## Number of iterations limit.

    maxIter = 100000

    #

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

    #  once for all.

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

        slz_compute_scaled_function(function, lowerBound, upperBound, precision)

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

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

    degreeSo  = pobyso_constant_from_int_sa_so(degree)

    #

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

    domainBoundsInterval   = RRIF(lowerBound, upperBound)

    (unscalingFunction, scalingFunction) = \

        slz_interval_scaling_expression(domainBoundsInterval, i)

    #print scalingFunction, unscalingFunction

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

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

    if internalSollyaPrec < 192:

        internalSollyaPrec = 192

    pobyso_set_prec_sa_so(internalSollyaPrec)

    print "Sollya internal precision:", internalSollyaPrec

    ## Some variables.

    ### General variables

    lb                  = sdlb

    ub                  = sdub

    nbw                 = 0

    intervalUlp         = ub.ulp()

    #### Will be set by slz_interval_and_polynomila_to_sage.

    ic                  = 0 

    icAsInt             = 0    # Set from ic.

    solutionsSet        = set()

    tsErrorWidth        = []

    csErrorVectors      = []

    csVectorsResultants = []

    floatP              = 0  # Taylor polynomial.

    floatPcv            = 0  # Ditto with variable change.

    intvl               = "" # Taylor interval

    terr                = 0  # Taylor error.

    iterCount = 0

    htrnSet = set()

    ### Timers and counters.

    wallTimeStart                   = 0

    cpuTimeStart                    = 0

    taylCondFailedCount             = 0

    coppCondFailedCount             = 0

    resultCondFailedCount           = 0

    coppCondFailed                  = False

    resultCondFailed                = False

    globalResultsList               = []

    basisConstructionsCount         = 0

    basisConstructionsFullTime      = 0

    basisConstructionTime           = 0

    reductionsCount                 = 0

    reductionsFullTime              = 0

    reductionTime                   = 0

    resultantsComputationsCount     = 0

    resultantsComputationsFullTime  = 0

    resultantsComputationTime       = 0

    rootsComputationsCount          = 0

    rootsComputationsFullTime       = 0

    rootsComputationTime            = 0

    print

    ## Global times are started here.

    wallTimeStart                   = walltime()

    cpuTimeStart                    = cputime()

    ## Main loop.

    while True:

        if lb >= sdub:

            print "Lower bound reached upper bound."

            break

        if iterCount == maxIter:

            print "Reached maxIter. Aborting"

            break

        iterCount += 1

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

            "log2(numbers)." 

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

        prceSo = slz_compute_polynomial_and_interval(scaledfSo,

                                                     degreeSo,

                                                     lb,

                                                     ub,

                                                     polyApproxAccur)

        ### Convert back the data into Sage space.                         

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

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

                                                 prceSo[1], prceSo[2], 

                                                 prceSo[3]))

        intvl = RRIF(intvl)

        ## Clean-up Sollya stuff.

        for elem in prceSo:

            sollya_lib_clear_obj(elem)

        #print  floatP, floatPcv, intvl, ic, terr

        #print floatP

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

        #      ic.n(),

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

        ### Check returned data.

        #### Is approximation error OK?

        if terr > polyApproxAccur:

            exceptionErrorMess  = \

                "Approximation failed - computed error:" + \

                str(terr) + " - target error: "

            exceptionErrorMess += \

                str(polyApproxAccur) + ". Aborting!"

            raise Exception(exceptionErrorMess)

        #### Is lower bound OK?

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

            exceptionErrorMess =  "Wrong lower bound:" + \

                               str(lb) + ". Aborting!"

            raise Exception(exceptionErrorMess)

        #### Set upper bound.

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

            ub = intvl.endpoints()[1]

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

            "log2(numbers)." 

            taylCondFailedCount += 1

        #### Is interval not degenerate?

        if lb >= ub:

            exceptionErrorMess = "Degenerate interval: " + \

                                "lowerBound(" + str(lb) +\

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

                                "). Aborting!"

            raise Exception(exceptionErrorMess)

        #### Is interval center ok?

        if ic <= lb or ic >= ub:

            exceptionErrorMess =  "Invalid interval center for " + \

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

                                str(ub) + ". Aborting!"

            raise Exception(exceptionErrorMess)

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

        bw = ub - lb

        nbw = 0

        icAsInt = int(ic * toIntegerFactor)

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

        #  enlarge the interval

        curTaylErrRat = polyApproxAccur / terr

        ## Make the  integral transformations.

        ### First for interval center and bounds.

        intIc = int(ic * toIntegerFactor)

        intLb = int(lb * toIntegerFactor) - intIc

        intUb = int(ub * toIntegerFactor) - intIc

        #

        #### Loop flesh

        #### For polynomials

        basisConstructionTime         = cputime()

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

        ratRatP = slz_float_poly_of_float_to_rat_poly_of_rat_pow_two(floatP)

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

        ratIntP = \

            slz_rat_poly_of_rat_to_rat_poly_of_int(ratRatP, precision)

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

        #      arguments.

        coppersmithTuple = \

            slz_rat_poly_of_int_to_poly_for_coppersmith(ratIntP, 

                                                        precision, 

                                                        targetHardnessToRound, 

                                                        i, t)

        #### Recover Coppersmith information.

        intIntP = coppersmithTuple[0]

        N = coppersmithTuple[1]

        nAtAlpha = N^alpha

        tBound = coppersmithTuple[2]

        leastCommonMultiple = coppersmithTuple[3]

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

        basisConstructionsFullTime        += cputime(basisConstructionTime)

        basisConstructionsCount           += 1

        reductionTime                     = cputime()

        # Compute the reduced polynomials.

        ccReducedPolynomialsList =  \

            slz_compute_coppersmith_reduced_polynomials(intIntP, 

                                                        alpha, 

                                                        N, 

                                                        iBound, 

                                                        tBound)

        if ccReducedPolynomialsList is None:

            raise Exception("Reduction failed.")

        reductionsFullTime    += cputime(reductionTime)

        reductionsCount       += 1

        if len(ccReducedPolynomialsList) < 2:

            print "Nothing to form resultants with."

            print

            coppCondFailedCount += 1

            coppCondFailed = True

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

            #  the number of complient polynomials.

            if len(ccReducedPolynomialsList) == 0:

                ub  = lb + bw / 4

            else: # At least one complient polynomial.

                ub = lb + bw / 2

            if ub > sdub:

                ub = sdub

            if lb == ub:

                raise Exception("Cant shrink interval \

                anymore to get Coppersmith condition.")

            nbw = 0

            continue

        #### We have at least two polynomials.

        #   Let us try to compute resultants.

        resultantsComputationTime      = cputime()

        resultantsInTTuplesList = []

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

            for polyInnerIndex in xrange(polyOuterIndex+1, 

                                         len(ccReducedPolynomialsList)):

                resultantTuple = \

                slz_resultant_tuple(ccReducedPolynomialsList[polyOuterIndex],

                                ccReducedPolynomialsList[polyInnerIndex],

                                t)

                if len(resultantTuple) > 2:                    

                    #print resultantTuple[2]

                    resultantsInTTuplesList.append(resultantTuple)

                else:

                    print "No non nul resultant"

        print len(resultantsInTTuplesList), "resultant(s) in t tuple(s) created."

        resultantsComputationsFullTime += cputime(resultantsComputationTime)

        resultantsComputationsCount    += 1

        if len(resultantsInTTuplesList) == 0:

            print "Only null resultants, shrinking interval."

            resultCondFailed      =  True

            resultCondFailedCount += 1

            ### Shrink interval for next iteration.

            ub = lb + bw / 2

            if ub > sdub:

                ub = sdub

            nbw = 0

            continue

        #### Compute roots.

        reducedPolynomialsRootsSet = set()

        ##### Solve in the second variable since resultants are in the first

        #     variable.

        for resultantInTTuple in resultantsInTTuplesList:

            currentResultant = resultantInTTuple[2]

            ##### If the resultant degree is not at least 1, there are no roots.

            if currentResultant.degree() < 1:

                print "Resultant is constant:", currentResultant

                continue # Next resultantInTTuple

            ##### Compute i roots

            iRootsList = Zi(currentResultant).roots()

            ##### For each iRoot, compute the corresponding tRoots and check 

            #     them in the input polynomial.

            for iRoot in iRootsList:

                ####### Roots returned by roots() are (value, multiplicity) 

                #       tuples.

                #print "iRoot:", iRoot

                ###### Use the tRoot against each polynomial, alternatively.

                for indexInTuple in range(0,2):

                    currentPolynomial = resultantInTTuple[indexInTuple]

                    ####### If the polynomial is univariate, just drop it.

                    if len(currentPolynomial.variables()) < 2:

                        print "    Current polynomial is not in two variables."

                        continue # Next indexInTuple

                    tRootsList = \

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

                    ####### The tRootsList can be empty, hence the test.

                    if len(tRootsList) == 0:

                        print "    No t root."

                        continue # Next indexInTuple

                    for tRoot in tRootsList:

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

        # End of roots computation

    ## Prepare for results.

    intervalResultsList = []

    intervalResultsList.append((lb, ub))

    ## Check roots.

    rootsResultsList = []

    rootsComputationTime        = cputime()

    for root in reducedPolynomialsRootsSet:

        specificRootResultsList = []

        failingBounds = []

        intIntPdivN = intIntP(root[0], root[1]) / N

        if int(intIntPdivN) != intIntPdivN:

            continue # Next root

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

        hardToRoundCaseAsFloat = RRR((icAsInt + root[0]) / toIntegerFactor)

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

        #print scalingFunction

        scaledHardToRoundCaseAsFloat = scalingFunction(hardToRoundCaseAsFloat) 

        print "Candidate HTRNc at x =", scaledHardToRoundCaseAsFloat.n().str(base=2),

        if slz_is_htrn(scaledHardToRoundCaseAsFloat,

                       f,

                       2^-(targetHardnessToRound),

                       RRR):

            print hardToRoundCaseAsFloat, "is HTRN case."

            if lb <= hardToRoundCaseAsFloat and hardToRoundCaseAsFloat <= ub:

                print "Found in interval."

            else:

                print "Found out of interval."

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

            # Check the root is in the bounds

            if abs(root[0]) > iBound or abs(root[1]) > tBound:

                print "Root", root, "is out of bounds."

                if abs(root[0]) > iBound:

                    print "root[0]:", root[0]

                    print "i bound:", iBound

                    failingBounds.append('i')

                    failingBounds.append(root[0])

                    failingBounds.append(iBound)

                if abs(root[1]) > tBound:

                    print "root[1]:", root[1]

                    print "t bound:", tBound

                    failingBounds.append('t')

                    failingBounds.append(root[1])

                    failingBounds.append(tBound)

            if len(failingBounds) > 0:

                specificRootResultsList.append(failingBounds)

        else: # From slz_is_htrn...

            print "is not an HTRN case."

        if len(specificRootResultsList) > 0:

            rootsResultsList.append(specificRootResultsList)

    rootsComputationsFullTime       =   cputime(rootsComputationTime)

    rootsComputationsCount          +=  1

    if len(rootsResultsList) > 0:

        intervalResultsList.append(rootsResultsList)

    ## An intervalResultsList has at least the bounds.

    globalResultsList.append(intervalResultsList)   

        #### End loop flesh.

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

        #    if not Coppersmith condition nor resultant condition

        #    failed at the previous run. 

        if not coppCondFailed and not resultCondFailed:

            nbw = (1 + 2^(-5)) * bw

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

            #     again if needed.

        else:

            nbw = bw

        coppCondFailed   = False

        resultCondFailed = False

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

        #print "nbw:", nbw

        lb  = ub

        ub += nbw     

        if ub > sdub:

            ub = sdub

        print

    # End while True

    ## Main loop just ended.

    globalWallTime = walltime(wallTimeStart)

    globalCpuTime  = cputime(cpuTimeStart)

    ## Output results

    print ; print "Intervals and HTRNs"

    for intervalResultsList in globalResultsList:

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

        if len(intervalResultsList) > 1:

            rootsResultsList = intervalResultsList[1]

            for specificRootResultsList in rootsResultsList:

                print "\t", specificRootResultsList[0],

                if len(specificRootResultsList) > 1:

                    print specificRootResultsList[1],

        print ; print

    #print globalResultsList

    #

    print "Timers and counters"

    print

    print "Number of iterations:", iterCount

    print "Taylor condition failures:", taylCondFailedCount

    print "Coppersmith condition failures:", coppCondFailedCount

    print "Resultant condition failures:", resultCondFailedCount

    print "Iterations count: ", iterCount

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

    print "Number of basis constructions:", basisConstructionsCount 

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

        basisConstructionsFullTime

    if basisConstructionsCount != 0:

        print "Average basis construction CPU time:", \

            basisConstructionsFullTime/basisConstructionsCount

    print "Number of reductions:", reductionsCount

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

    if reductionsCount != 0:

        print "Average reduction CPU time:", \

            reductionsFullTime/reductionsCount

    print "Number of resultants computation rounds:", \

        resultantsComputationsCount

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

        resultantsComputationsFullTime

    if resultantsComputationsCount != 0:

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

            resultantsComputationsFullTime/resultantsComputationsCount

    print "Number of root finding rounds:", rootsComputationsCount

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

        rootsComputationsFullTime

    if rootsComputationsCount != 0:

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

            rootsComputationsFullTime/rootsComputationsCount

    print "Global Wall time:", globalWallTime

    print "Global CPU time:", globalCpuTime

    ## Output counters

# End runSLZ-v01

print "Running SLZ..."

initialize_env()

x = var('x')

func(x) = exp(x)

precision = 53

RRR = RealField(precision)

run_SLZ_v01(inputFunction=func, 

            inputLowerBound = 1/4, 

            inputUpperBound = RRR(1/2) - RRR(1/4).ulp(), 

            alpha = 2, 

            degree = 10, 

            precision = 53, 

            emin = -1022, 

            emax = 1023, 

            targetHardnessToRound =  precision+50, 

            debug = True)