Index: runSLZ-01.sage
===================================================================
--- runSLZ-01.sage	(revision 188)
+++ runSLZ-01.sage	(revision 189)
@@ -1,11 +1,12 @@
 #! /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/sageMpfr.spyx")
     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")
@@ -13,25 +14,291 @@
     load("/home/storres/recherche/arithmetique/pobysoPythonSage/src/sageSLZ/sagePolynomialOperations.sage")
 
 
-def run_SLZ(function,
-            lowerBound,
-            upperBound,
+def run_SLZ(inputFunction,
+            inputLowerBound,
+            inputUpperBound,
             alpha,
             degree,
             precision,
+            emin,
+            emax,
             targetHardnessToRound,
             debug = False):
 
     if debug:
-        print "Function                :", function
-        print "Lower bound             :", lowerBound
-        print "Upper bounds            :", upperBound
+        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()
-run_SLZ(exp(x), 1/4, 1/2, 2, 2, 53, 100, True)
-
+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)
