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

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

        rootsComputationsFullTime       =   cputime(rootsComputationTime)

        rootsComputationsCount          +=  1

        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

print "Running SLZ..."

initialize_env()

x = var('x')

func(x) = exp(x)

precision = 53

RRR = RealField(precision)

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