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

def slz_compute_coppersmith_reduced_polynomials_with_lattice_volume(inputPolynomial,

                                                                    alpha,

                                                                    N,

                                                                    iBound,

                                                                    tBound):

    """

    For a given set of arguments (see below), compute a list

    of "reduced polynomials" that could be used to compute roots

    of the inputPolynomial.

    Print the volume of the initial basis as well.

    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,

##        spo_polynomial_to_polynomials_list_5(initialPolynomial,

    polynomialsList = \

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

    matrixToReduceTranspose = matrixToReduce.transpose()

    squareMatrix = matrixToReduce * matrixToReduceTranspose

    squareMatDet = det(squareMatrix)

    latticeVolume = sqrt(squareMatDet)

    print "Lattice volume:", latticeVolume.n()

    print "Lattice volume / N:", (latticeVolume/N).n()

    #print matrixToReduce

    ## Reduction and checking.

    ## S.T. changed 'fp' to None as of Sage 6.6 complying to

    #  error message issued when previous code was used.

    #reducedMatrix = matrixToReduce.LLL(fp='fp')

    reductionTimeStart = cputime() 

    reducedMatrix = matrixToReduce.LLL(fp=None)

    reductionTime = cputime(reductionTimeStart)

    print "Reduction time:", reductionTime

    isLLLReduced  = reducedMatrix.is_LLL_reduced()

    if not isLLLReduced:

        return None

    monomialsCount     = len(knownMonomials)

    monomialsCountSqrt = sqrt(monomialsCount)

    #print "Monomials count:", monomialsCount, monomialsCountSqrt.n()

    #print reducedMatrix

    ## Check the Coppersmith condition for each row and build the reduced 

    #  polynomials.

    ccReducedPolynomialsList = []

    for row in reducedMatrix.rows():

        l2Norm = row.norm(2)

        if (l2Norm * monomialsCountSqrt) < nAtAlpha:

            #print (l2Norm * monomialsCountSqrt).n()

            #print l2Norm.n()

            ccReducedPolynomial = \

                slz_compute_reduced_polynomial(row,

                                               knownMonomials,

                                               iVariable,

                                               iBound,

                                               tVariable,

                                               tBound)

            if not ccReducedPolynomial is None:

                ccReducedPolynomialsList.append(ccReducedPolynomial)

        else:

            #print l2Norm.n() , ">", nAtAlpha

            pass

    if len(ccReducedPolynomialsList) < 2:

        print "Less than 2 Coppersmith condition compliant vectors."

        return ()

    #print ccReducedPolynomialsList

    return ccReducedPolynomialsList

# End slz_compute_coppersmith_reduced_polynomials_with_lattice volume

intervalRadiusInUlp = 2^49 + 2^45          

srs_run_SLZ_v03(inputFunction=func, 

                inputLowerBound = RRR(1) * 2^-31, 

                inputUpperBound = RRR(1) * 2^-30 - icUlp, 

#! /opt/sage/sage

# @file runSLZ-03.sage

#

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")

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

print "Running SLZ..."

initialize_env()

x = var('x')

func(x) = exp(x)

precision = 53

RRR = RealField(precision)

intervalCenter      = RRR("1.9E9CBBFD6080B",16)  * 2^-31

icUlp               = intervalCenter.ulp()

intervalRadiusInUlp = 2^49 + 2^45          

srs_compute_lattice_volume(inputFunction=func, 

                           inputLowerBound = RRR(1) * 2^-31, 

                           inputUpperBound = RRR(1) * 2^-30 - icUlp, 

                           alpha           = 4, 

                           degree          = 4, 

                           precision       = 53, 

                           emin            = -1022, 

                           emax            = 1023, 

                           targetHardnessToRound =  precision+50, 

                           debug           = True)

#

"""

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

"""

"""

srs_run_SLZ_v01(inputFunction=func, 

                inputLowerBound = 402653184/1073741824, 

                inputUpperBound = 402653185/1073741824, 

                alpha = 2, 

                degree = 10, 

                precision = 53, 

                emin = -1022, 

                emax = 1023, 

                targetHardnessToRound =  precision+50, 

                debug = True)

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

RR("1.9E9CBBFD6080B",16)  * 2^-31]

"""

                slz_compute_coppersmith_reduced_polynomials(intIntP, 

                                                            alpha, 

                                                            N, 

                                                            iBound, 

                                                            tBound)

                slz_compute_coppersmith_reduced_polynomials(intIntP, 

                                                            alpha, 

                                                            N, 

                                                            iBound, 

                                                            tBound)

def srs_run_SLZ_v03(inputFunction,

                    inputLowerBound,

                    inputUpperBound,

                    alpha,

                    degree,

                    precision,

                    emin,

                    emax,

                    targetHardnessToRound,

                    debug = False):

    """

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

        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.

        rootsSet            = set()

        hasNonNullResultant = False

        for polyPair in polyPairsList:

            if hasNonNullResultant:

                break

            resultantsComputationTime   = cputime()

            currentResultantI = \

                slz_resultant(polyPair[0],

                              polyPair[1],

                              t)

            resultantsComputationsCount    += 1

            if currentResultantI is None:

                resultantsComputationsFullTime +=  \

                    cputime(resultantsComputationTime)

                print "Nul resultant"

                continue # Next polyPair.

            currentResultantT = \

                slz_resultant(polyPair[0],

                              polyPair[1],

                              i)

            resultantsComputationsFullTime += cputime(resultantsComputationTime)

            resultantsComputationsCount    += 1

            if currentResultantT is None:                    

                print "Nul resultant"

                continue # Next polyPair.

            else:

                hasNonNullResultant = True

            #### We have a non null resultants pair. From now on, whatever the

            #    root search yields, no extra root search is necessary.

            #### A constant resultant leads to no root. Root search is done.

            if currentResultantI.degree() < 1:

                print "Resultant is constant:", currentResultantI

                break # Next polyPair and should break.

            if currentResultantT.degree() < 1:

                print "Resultant is constant:", currentResultantT

                break # Next polyPair and should break.

            #### Actual roots computation.

            rootsComputationTime       = cputime()

            ##### Compute i roots

            iRootsList = Zi(currentResultantI).roots()

            rootsComputationsCount      +=  1

            if len(iRootsList) == 0:

                rootsComputationsFullTime   =   cputime(rootsComputationTime)

                print "No roots in \"i\"."

                break # No roots in i.

            tRootsList = Zt(currentResultantT).roots()

            rootsComputationsFullTime   =   cputime(rootsComputationTime)

            rootsComputationsCount      +=  1

            if len(tRootsList) == 0:

                print "No roots in \"t\"."

                break # No roots in i.

            ##### For each iRoot, get a tRoot and check against the polynomials.

            for iRoot in iRootsList:

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

                #       tuples.

                #print "iRoot:", iRoot

                for tRoot in tRootsList:

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

                    if polyPair[0](iRoot[0],tRoot[0]) != 0:

                        continue

                    if polyPair[1](iRoot[0],tRoot[0]) != 0:

                        continue

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

            # End of roots computation.

        # End loop for polyPair in polyParsList.  Will break at next iteration.

        # since a non null resultant was found.

        #### Prepare for results for the current interval..

        intervalResultsList = []

        intervalResultsList.append((lb, ub))

        #### Check roots.

        rootsResultsList = []

        for root in rootsSet:

            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,

                           function,

                           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 for modular equation."

                    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)

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

def srs_compute_lattice_volume(inputFunction,

                    inputLowerBound,

                    inputUpperBound,

                    alpha,

                    degree,

                    precision,

                    emin,

                    emax,

                    targetHardnessToRound,

                    debug = False):

    """

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

        reductionTime                     = cputime()

        #### Compute the reduced polynomials.

        ccReducedPolynomialsList =  \

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

        rootsSet            = set()

        hasNonNullResultant = False

        for polyPair in polyPairsList:

            if hasNonNullResultant:

                break

            resultantsComputationTime   = cputime()

            currentResultantI = \

                slz_resultant(polyPair[0],

                              polyPair[1],

                              t)

            resultantsComputationsCount    += 1

            if currentResultantI is None:

                resultantsComputationsFullTime +=  \

                    cputime(resultantsComputationTime)

                print "Nul resultant"

                continue # Next polyPair.

            currentResultantT = \

                slz_resultant(polyPair[0],

                              polyPair[1],

                              i)

            resultantsComputationsFullTime += cputime(resultantsComputationTime)

            resultantsComputationsCount    += 1

            if currentResultantT is None:                    

                print "Nul resultant"

                continue # Next polyPair.

            else:

                hasNonNullResultant = True

            #### We have a non null resultants pair. From now on, whatever the

            #    root search yields, no extra root search is necessary.

            #### A constant resultant leads to no root. Root search is done.

            if currentResultantI.degree() < 1:

                print "Resultant is constant:", currentResultantI

                break # Next polyPair and should break.

            if currentResultantT.degree() < 1:

                print "Resultant is constant:", currentResultantT

                break # Next polyPair and should break.

            #### Actual roots computation.

            rootsComputationTime       = cputime()

            ##### Compute i roots

            iRootsList = Zi(currentResultantI).roots()

            rootsComputationsCount      +=  1

            if len(iRootsList) == 0:

                rootsComputationsFullTime   =   cputime(rootsComputationTime)

                print "No roots in \"i\"."

                break # No roots in i.

            tRootsList = Zt(currentResultantT).roots()

            rootsComputationsFullTime   =   cputime(rootsComputationTime)

            rootsComputationsCount      +=  1

            if len(tRootsList) == 0:

                print "No roots in \"t\"."

                break # No roots in i.

            ##### For each iRoot, get a tRoot and check against the polynomials.

            for iRoot in iRootsList:

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

                #       tuples.

                #print "iRoot:", iRoot

                for tRoot in tRootsList:

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

                    if polyPair[0](iRoot[0],tRoot[0]) != 0:

                        continue

                    if polyPair[1](iRoot[0],tRoot[0]) != 0:

                        continue

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

            # End of roots computation.

        # End loop for polyPair in polyParsList.  Will break at next iteration.

        # since a non null resultant was found.

        #### Prepare for results for the current interval..

        intervalResultsList = []

        intervalResultsList.append((lb, ub))

        #### Check roots.

        rootsResultsList = []

        for root in rootsSet:

            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,

                           function,

                           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 for modular equation."

                    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)

        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_compute_lattice_volume