Index: pobysoPythonSage/src/pobyso.py
===================================================================
--- pobysoPythonSage/src/pobyso.py	(revision 218)
+++ pobysoPythonSage/src/pobyso.py	(revision 219)
@@ -1128,17 +1128,19 @@
                                                            itpSo,
                                                            ftpSo,
                                                            maxPrecSo,
-                                                           maxErrSo):
-    print "Input arguments:"
-    #pobyso_autoprint(polySo)
-    #pobyso_autoprint(funcSo)
-    #pobyso_autoprint(icSo)
-    #pobyso_autoprint(intervalSo)
-    #pobyso_autoprint(itpSo)
-    #pobyso_autoprint(ftpSo)
-    #pobyso_autoprint(maxPrecSo)
-    #pobyso_autoprint(maxErrSo)
-    #print "________________"
+                                                           maxErrSo,
+                                                           debug=False):
+    if debug:
+        print "Input arguments:"
+        pobyso_autoprint(polySo)
+        pobyso_autoprint(funcSo)
+        pobyso_autoprint(icSo)
+        pobyso_autoprint(intervalSo)
+        pobyso_autoprint(itpSo)
+        pobyso_autoprint(ftpSo)
+        pobyso_autoprint(maxPrecSo)
+        pobyso_autoprint(maxErrSo)
+        print "________________"
     
     ## Higher order function see: 
     #  http://effbot.org/pyfaq/how-do-you-make-a-higher-order-function-in-python.htm
@@ -1155,20 +1157,21 @@
     
     #
     degreeSa        = pobyso_polynomial_degree_so_sa(polySo)
-    print "degreeSa:", degreeSa
     ratio           = precision_decay_ratio_function(degreeSa)
-    print "ratio:", ratio
     itpSa           = pobyso_constant_from_int_so_sa(itpSo)
-    print "itpsSa:", itpSa
     ftpSa           = pobyso_constant_from_int_so_sa(ftpSo)
-    print "ftpSa:", ftpSa
     maxPrecSa       = pobyso_constant_from_int_so_sa(maxPrecSo)
-    print "maxPrecSa:", maxPrecSa
     maxErrSa        = pobyso_get_constant_as_rn_so_sa(maxErrSo)
-    print "maxErrSa:", maxErrSa
+    if debug:
+        print "degreeSa:", degreeSa
+        print "ratio:", ratio
+        print "itpsSa:", itpSa
+        print "ftpSa:", ftpSa
+        print "maxPrecSa:", maxPrecSa
+        print "maxErrSa:", maxErrSa
     lastResPolySo   = None
     lastInfNormSo   = None
-    print "About to enter the while loop..."
+    #print "About to enter the while loop..."
     while True: 
         resPolySo   = pobyso_constant_0_sa_so()
         pDeltaSa    = ftpSa - itpSa
@@ -1179,8 +1182,9 @@
             #print ratio(indexSa)
             ctpSa = floor(ftpSa - (pDeltaSa * ratio(indexSa)))
             ctpSo = pobyso_constant_from_int_sa_so(ctpSa)
-            print "Index:", indexSa, " - Target precision:", 
-            pobyso_autoprint(ctpSo)
+            if debug:
+                print "Index:", indexSa, " - Target precision:", 
+                pobyso_autoprint(ctpSo)
             cmonSo  = \
                 sollya_lib_build_function_mul(sollya_lib_coeff(polySo, indexSo),
                                       sollya_lib_build_function_pow( \
@@ -1201,13 +1205,17 @@
                                                   resPolyCvSo)
         infNormSo = sollya_lib_dirtyinfnorm(errFuncSo, intervalSo)
         cerrSa    = pobyso_get_constant_as_rn_so_sa(infNormSo)
-        print "Infnorm (Sollya):", pobyso_autoprint(infNormSo)
+        if debug:
+            print "Infnorm (Sollya):", 
+            pobyso_autoprint(infNormSo)
         sollya_lib_clear_obj(errFuncSo)
         #print "Infnorm  (Sage):", cerrSa
         if (cerrSa > maxErrSa):
-            print "Error is too large."
+            if debug:
+                print "Error is too large."
             if lastResPolySo is None:
-                print "Enlarging prec."
+                if debug:
+                    print "Enlarging prec."
                 ntpSa = floor(ftpSa + ftpSa/50)
                 ## Can't enlarge (numerical)
                 if ntpSa == ftpSa:
@@ -1221,15 +1229,18 @@
                 continue
             ## One enlargement took place.
             else:
-                print "Exit with the last before last polynomial."
+                if debug:
+                    print "Exit with the last before last polynomial."
                 sollya_lib_clear_obj(resPolySo)
                 sollya_lib_clear_obj(infNormSo)
                 return (lastResPolySo, lastInfNormSo)
         # cerrSa <= maxErrSa: scrap more bits, possibly.
         else:
-            print "Error is too small"
-            if cerrSa <= (maxErrSa/2): 
-                print "Shrinking prec."
+            if debug:
+                print "Error is too small"
+            if cerrSa <= (maxErrSa/2):
+                if debug: 
+                    print "Shrinking prec."
                 ntpSa = floor(ftpSa - ftpSa/50)
                 ## Can't shrink (numerical)
                 if ntpSa == ftpSa:
Index: pobysoPythonSage/src/sageSLZ/sageRunSLZ.sage
===================================================================
--- pobysoPythonSage/src/sageSLZ/sageRunSLZ.sage	(revision 218)
+++ pobysoPythonSage/src/sageSLZ/sageRunSLZ.sage	(revision 219)
@@ -1950,7 +1950,7 @@
     Changes from V3:
         Root search is changed again:
             - only resultants in i are computed;
-            - root are searched for;
+            - roots in i are searched for;
             - if any, they are tested for hardness-to-round.
     Changes from V2:
         Root search is changed:
@@ -2396,3 +2396,477 @@
     ## Output counters
 # End srs_runSLZ-v04
 
+def srs_run_SLZ_v05(inputFunction,
+                    inputLowerBound,
+                    inputUpperBound,
+                    alpha,
+                    degree,
+                    precision,
+                    emin,
+                    emax,
+                    targetHardnessToRound,
+                    debug = False):
+    """
+    Changes from V4:
+        Approximation polynomial has coefficients rounded.
+    Changes from V3:
+        Root search is changed again:
+            - only resultants in i are computed;
+            - roots in i are searched for;
+            - if any, they are tested for hardness-to-round.
+    Changes from V2:
+        Root search is changed:
+            - we compute the resultants in i and in t;
+            - we compute the roots set of each of these resultants;
+            - we combine all the possible pairs between the two sets;
+            - we check these pairs in polynomials for correctness.
+    Changes from V1: 
+        1- check for roots as soon as a resultant is computed;
+        2- once a non null resultant is found, check for roots;
+        3- constant resultant == no root.
+    """
+
+    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
+    ## Important constants.
+    ### Stretch the interval if no error happens.
+    noErrorIntervalStretch = 1 + 2^(-5)
+    ### If no vector validates the Coppersmith condition, shrink the interval
+    #   by the following factor.
+    noCoppersmithIntervalShrink = 1/2
+    ### If only (or at least) one vector validates the Coppersmith condition, 
+    #   shrink the interval by the following factor.
+    oneCoppersmithIntervalShrink = 3/4
+    #### If only null resultants are found, shrink the interval by the 
+    #    following factor.
+    onlyNullResultantsShrink     = 3/4
+    ## 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()
+    output = slz_fix_bounds_for_binades(lowerBound, upperBound, inputFunction)
+    if output is None:
+        print "Invalid domain/image binades. Domain:",\
+        lowerBound, upperBound, "Images:", \
+        inputFunction(lowerBound), inputFunction(upperBound)
+        raise Exception("Invalid domain/image binades.")
+    lb = output[0] ; ub = output[1]
+    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_01(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.
+        #### Bounds and interval center.
+        intIc = int(ic * toIntegerFactor)
+        intLb = int(lb * toIntegerFactor) - intIc
+        intUb = int(ub * toIntegerFactor) - intIc
+        #
+        #### 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 multivariate 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
+        """
+        #### Compute the matrix to reduce.
+        matrixToReduce = slz_compute_initial_lattice_matrix(intIntP,
+                                                            alpha,
+                                                            N,
+                                                            iBound,
+                                                            tBound)
+        matrixFile = file('/tmp/matrixToReduce.txt', 'w')
+        for row in matrixToReduce.rows():
+            matrixFile.write(str(row) + "\n")
+        matrixFile.close()
+        raise Exception("Deliberate stop here.")
+        """
+        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 compliant polynomials.
+            if len(ccReducedPolynomialsList) == 0:
+                ub = lb + bw * noCoppersmithIntervalShrink
+            else: # At least one compliant polynomial.
+                ub = lb + bw * oneCoppersmithIntervalShrink
+            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.
+        #    For each resultant computed, go for the solutions.
+        ##### Build the pairs list.
+        polyPairsList           = []
+        for polyOuterIndex in xrange(0, len(ccReducedPolynomialsList) - 1):
+            for polyInnerIndex in xrange(polyOuterIndex+1, 
+                                         len(ccReducedPolynomialsList)):
+                polyPairsList.append((ccReducedPolynomialsList[polyOuterIndex],
+                                      ccReducedPolynomialsList[polyInnerIndex]))
+        #### Actual root search.
+        iRootsSet           = set()
+        hasNonNullResultant = False
+        for polyPair in polyPairsList:
+            resultantsComputationTime   = cputime()
+            currentResultantI = \
+                slz_resultant(polyPair[0],
+                              polyPair[1],
+                              t)
+            resultantsComputationsCount    += 1
+            resultantsComputationsFullTime +=  \
+                cputime(resultantsComputationTime)
+            #### Function slz_resultant returns None both for None and O
+            #    resultants.
+            if currentResultantI is None:
+                print "Nul resultant"
+                continue # Next polyPair.
+            ## We deleted the currentResultantI computation.
+            #### We have a non null resultant. From now on, whatever this
+            #    root search yields, no extra root search is necessary.
+            hasNonNullResultant = True
+            #### A constant resultant leads to no root. Root search is done.
+            if currentResultantI.degree() < 1:
+                print "Resultant is constant:", currentResultantI
+                break # There is no root.
+            #### Actual iroots computation.
+            rootsComputationTime        = cputime()
+            iRootsList = Zi(currentResultantI).roots()
+            rootsComputationsCount      +=  1
+            rootsComputationsFullTime   =   cputime(rootsComputationTime)
+            if len(iRootsList) == 0:
+                print "No roots in \"i\"."
+                break # No roots in i.
+            else:
+                for iRoot in iRootsList:
+                    # A root is given as a (value, multiplicity) tuple.
+                    iRootsSet.add(iRoot[0])
+        # End loop for polyPair in polyParsList. We only loop again if a 
+        # None or zero resultant is found.
+        #### Prepare for results for the current interval..
+        intervalResultsList = []
+        intervalResultsList.append((lb, ub))
+        #### Check roots.
+        rootsResultsList = []
+        for iRoot in iRootsSet:
+            specificRootResultsList = []
+            failingBounds           = []
+            # Root qualifies for modular equation, test it for hardness to round.
+            hardToRoundCaseAsFloat = RRR((icAsInt + iRoot) / 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,
+                           function,
+                           2^-(targetHardnessToRound),
+                           RRR):
+                print hardToRoundCaseAsFloat, "is HTRN case."
+                specificRootResultsList.append(hardToRoundCaseAsFloat.n().str(base=2))
+                if lb <= hardToRoundCaseAsFloat and hardToRoundCaseAsFloat <= ub:
+                    print "Found in interval."
+                else:
+                    print "Found out of interval."
+                # Check the i root is within the i bound.
+                if abs(iRoot) > iBound:
+                    print "IRoot", iRoot, "is out of bounds for modular equation."
+                    print "i bound:", iBound
+                    failingBounds.append('i')
+                    failingBounds.append(iRoot)
+                    failingBounds.append(iBound)
+                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)
+        if len(rootsResultsList) > 0:
+            intervalResultsList.append(rootsResultsList)
+        ### Check if a non null resultant was found. If not shrink the interval.
+        if not hasNonNullResultant:
+            print "Only null resultants for this reduction, shrinking interval."
+            resultCondFailed      =  True
+            resultCondFailedCount += 1
+            ### Shrink interval for next iteration.
+            ub = lb + bw * onlyNullResultantsShrink
+            if ub > sdub:
+                ub = sdub
+            nbw = 0
+            continue
+        #### An intervalResultsList has at least the bounds.
+        globalResultsList.append(intervalResultsList)   
+        #### 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 = noErrorIntervalStretch * bw
+        else:
+            nbw = bw
+        ##### Reset the failure flags. They will be raised
+        #     again if needed.
+        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" ; print
+    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 srs_runSLZ-v05
Index: pobysoPythonSage/src/sageSLZ/runSLZ-04.sage
===================================================================
--- pobysoPythonSage/src/sageSLZ/runSLZ-04.sage	(revision 218)
+++ pobysoPythonSage/src/sageSLZ/runSLZ-04.sage	(revision 219)
@@ -22,21 +22,22 @@
 initialize_env()
 x = var('x')
 func(x) = exp(x)
-precision = 53
+precision = 113
 RRR = RealField(precision)
-intervalCenter      = RRR("1.9E9CBBFD6080B",16)  * 2^-31
+intervalCenter      = RRR(3/8)
 icUlp               = intervalCenter.ulp()
-intervalRadiusInUlp = 2^49 + 2^45          
+intervalRadiusInUlp = 2^49 + 2^45
+intervalRadius      = RRR(2^(-35))          
 srs_run_SLZ_v04(inputFunction=func, 
-                inputLowerBound = RRR(1) * 2^-31, 
-                inputUpperBound = RRR(1) * 2^-30 - icUlp, 
-                alpha           = 2, 
-                degree          = 2, 
-                precision       = 53, 
-                emin            = -1022, 
-                emax            = 1023, 
-                targetHardnessToRound =  precision+50, 
-                debug           = True)
+                inputLowerBound         = intervalCenter - intervalRadius, 
+                inputUpperBound         = intervalCenter + intervalRadius, 
+                alpha                   = 6, 
+                degree                  = 18, 
+                precision               = precision, 
+                emin                    = -16382, 
+                emax                    = 16383, 
+                targetHardnessToRound   = precision * 6, 
+                debug                   = True)
 #
 """
 srs_run_SLZ_v01(inputFunction=func, 
Index: pobysoPythonSage/src/sageSLZ/sageSLZ.sage
===================================================================
--- pobysoPythonSage/src/sageSLZ/sageSLZ.sage	(revision 218)
+++ pobysoPythonSage/src/sageSLZ/sageSLZ.sage	(revision 219)
@@ -336,6 +336,61 @@
     return ccReducedPolynomialsList
 # End slz_compute_coppersmith_reduced_polynomials_with_lattice volume
 
+def slz_compute_initial_lattice_matrix(inputPolynomial,
+                                       alpha,
+                                       N,
+                                       iBound,
+                                       tBound):
+    """
+    For a given set of arguments (see below), compute the initial lattice
+    that could be reduced. 
+    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 None
+    # 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_8(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)
+    return matrixToReduce
+# End slz_compute_initial_lattice_matrix.
+
 def slz_compute_integer_polynomial_modular_roots(inputPolynomial,
                                                  alpha,
                                                  N,
@@ -610,6 +665,149 @@
 
 def slz_compute_polynomial_and_interval_01(functionSo, degreeSo, lowerBoundSa, 
                                         upperBoundSa, approxAccurSa, 
+                                        sollyaPrecSa=None, debug=False):
+    """
+    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 <= approxAccurSaS.
+    
+    """
+    #print"In slz_compute_polynomial_and_interval..."
+    ## Superficial argument check.
+    if lowerBoundSa > upperBoundSa:
+        return None
+    ## Change Sollya precision, if requested.
+    if not sollyaPrecSa is None:
+        sollyaPrecSo = pobyso_get_prec_so()
+        pobyso_set_prec_sa_so(sollyaPrecSa)
+    else:
+        sollyaPrecSa = pobyso_get_prec_so_sa()
+        sollyaPrecSo = None
+    RRR = lowerBoundSa.parent()
+    intervalShrinkConstFactorSa = RRR('0.9')
+    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 > approxAccurSa:
+        print "++Approximation error:", maxErrorSa.n()
+        sollya_lib_clear_obj(polySo)
+        sollya_lib_clear_obj(intervalCenterSo)
+        sollya_lib_clear_obj(maxErrorSo)
+        # Very empirical shrinking factor.
+        shrinkFactorSa = 1 /  (maxErrorSa/approxAccurSa).log2().abs()
+        print "Shrink factor:", \
+              shrinkFactorSa.n(), \
+              intervalShrinkConstFactorSa
+        print
+        #errorRatioSa = approxAccurSa/maxErrorSa
+        #print "Error ratio: ", errorRatioSa
+        # Make sure interval shrinks.
+        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)
+        # Check what is left with the bounds.
+        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!"
+            if not sollyaPrecSo is None:
+                pobyso_set_prec_so_so(sollyaPrecSo)
+                sollya_lib_clear_obj(sollyaPrecSo)
+            return 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))
+    # End while
+    sollya_lib_clear_obj(absoluteErrorTypeSo)
+    itpSo           = pobyso_constant_from_int_sa_so(floor(sollyaPrecSa/3))
+    ftpSo           = pobyso_constant_from_int_sa_so(floor(2*sollyaPrecSa/3))
+    maxPrecSo       = pobyso_constant_from_int_sa_so(sollyaPrecSa)
+    approxAccurSo   = pobyso_constant_sa_so(RR(approxAccurSa))
+    if debug:
+        print "About to call polynomial rounding with:"
+        print "polySo: ",           ; pobyso_autoprint(polySo)
+        print "functionSo: ",       ; pobyso_autoprint(functionSo)
+        print "intervalCenterSo: ", ; pobyso_autoprint(intervalCenterSo)
+        print "currentRangeSo: ",   ; pobyso_autoprint(currentRangeSo)
+        print "itpSo: ",            ; pobyso_autoprint(itpSo)
+        print "ftpSo: ",            ; pobyso_autoprint(ftpSo)
+        print "maxPrecSo: ",        ; pobyso_autoprint(maxPrecSo)
+        print "approxAccurSo: ",    ; pobyso_autoprint(approxAccurSo)
+    (roundedPolySo, roundedPolyMaxErrSo) = \
+    pobyso_polynomial_coefficients_progressive_round_so_so(polySo,
+                                                           functionSo,
+                                                           intervalCenterSo,
+                                                           currentRangeSo,
+                                                           itpSo,
+                                                           ftpSo,
+                                                           maxPrecSo,
+                                                           approxAccurSo)
+    
+    sollya_lib_clear_obj(polySo)
+    sollya_lib_clear_obj(maxErrorSo)
+    sollya_lib_clear_obj(itpSo)
+    sollya_lib_clear_obj(ftpSo)
+    sollya_lib_clear_obj(approxAccurSo)
+    if not sollyaPrecSo is None:
+        pobyso_set_prec_so_so(sollyaPrecSo)
+        sollya_lib_clear_obj(sollyaPrecSo)
+    if debug:
+        print "1: ", ; pobyso_autoprint(roundedPolySo)
+        print "2: ", ; pobyso_autoprint(currentRangeSo)
+        print "3: ", ; pobyso_autoprint(intervalCenterSo)
+        print "4: ", ; pobyso_autoprint(roundedPolyMaxErrSo)
+    return (roundedPolySo, currentRangeSo, intervalCenterSo, roundedPolyMaxErrSo)
+# End slz_compute_polynomial_and_interval_01
+
+def slz_compute_polynomial_and_interval_02(functionSo, degreeSo, lowerBoundSa, 
+                                        upperBoundSa, approxAccurSa, 
                                         sollyaPrecSa=None):
     """
     Under the assumptions listed for slz_get_intervals_and_polynomials, compute
@@ -747,7 +945,7 @@
     print "3: ", ; pobyso_autoprint(intervalCenterSo)
     print "4: ", ; pobyso_autoprint(roundedPolyMaxErrSo)
     return (roundedPolySo, currentRangeSo, intervalCenterSo, roundedPolyMaxErrSo)
-# End slz_compute_polynomial_and_interval_01
+# End slz_compute_polynomial_and_interval_02
 
 def slz_compute_reduced_polynomial(matrixRow,
                                     knownMonomials,
