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

# @file runSLZ-113gram.sage

#@par Changes from runSLZ-113.sage

# LLL reduction is not performed on the matrix itself but on its Gram matrix as

# an attempt for a faster reduction.

#

def slz_compute_coppersmith_reduced_polynomials_proj(inputPolynomial,

                                                     alpha,

                                                     N,

                                                     iBound,

                                                     tBound,

                                                     debug = False):

    """

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

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

    of the inputPolynomial.

    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.

    """

    #@par Changes from runSLZ-113.sage

    # LLL reduction is not performed on the matrix itself but rather on the

    # product of the matrix with a uniform random matrix.

    # The reduced matrix obtained is discarded but the transformation matrix 

    # obtained is used to multiply the original matrix in order to reduced it.

    # If a sufficient level of reduction is obtained, we stop here. If not

    # the product matrix obtained above is LLL reduced. But as it has been

    # pre-reduced at the above step, reduction is supposed to be much faster.

    #

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

    if debug:

        polynomialsList = \

            spo_polynomial_to_polynomials_list_8(initialPolynomial,

                                                 alpha,

                                                 N,

                                                 iBound,

                                                 tBound,

                                                 20)

    else:

        polynomialsList = \

            spo_polynomial_to_polynomials_list_8(initialPolynomial,

                                                 alpha,

                                                 N,

                                                 iBound,

                                                 tBound,

                                                 0)

    #print "Polynomials list:", polynomialsList

    ## Building the proto matrix.

    knownMonomials = []

    protoMatrix    = []

    if debug:

        for poly in polynomialsList:

            spo_add_polynomial_coeffs_to_matrix_row(poly,

                                                    knownMonomials,

                                                    protoMatrix,

                                                    20)

    else:

        for poly in polynomialsList:

            spo_add_polynomial_coeffs_to_matrix_row(poly,

                                                    knownMonomials,

                                                    protoMatrix,

                                                    0)

    matrixToReduce = spo_proto_to_row_matrix(protoMatrix)

    #print matrixToReduce

    ## Reduction and checking.

    ### Reduction with projection

    (reducedMatrix1, reductionMatrix1) = \

         slz_reduce_lll_proj(matrixToReduce, 2^16)

    #print "Reduced matrix:"

    matrixIsReducedEnough = reducedMatrix1.is_LLL_reduced()

    #

    if not matrixIsReducedEnough:

        reducedMatrix = reducedMatrix1.LLL(algorithm='fpLLL:wrapper')

        isLLLReduced  = reducedMatrix.is_LLL_reduced()

    else:

        reducedMatrix = reducedMatrix1

        isLLLReduced  = matrixIsReducedEnough

    #for row in reducedMatrix.rows():

    #    print row

    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_proj

#

#

def slz_pm1():

    """

    Compute a uniform RV in {-1, 1}.

    """

    ## function getrandbits(N) generates a long int with N random bits.

    return getrandbits(1) * 2-1

# End slz_pm1.

#

def slz_random_proj_pm1(r, c):

    """

    r x c matrix with \pm 1 ({-1, 1}) coefficients

    """

    M = matrix(r, c)

    for i in range(0, r):

        for j in range (0, c):

            M[i,j] = slz_pm1()

    return M

# End random_proj_pm1.

#

def slz_random_proj_uniform(n, r, c):

    """

    r x c integer matrix with uniform n-bit integer coefficients

    """

    M = matrix(r, c)

    for i in range(0, r):

        for j in range (0, c):

            M[i,j] = slz_uniform(n)

    return M       

# End slz_random_proj_uniform.

#

def slz_reduce_lll_proj(matToReduce, n):

    """

    Compute the transformation matrix that realizes an LLL reduction on

    the random uniform projected matrix.

    Return both the initial matrix "reduced" by the transformation matrix and

    the transformation matrix itself.

    """

    ## Compute the projected matrix.

    matProjector = slz_random_proj_uniform(n,

                                           matToReduce.ncols(),

                                           matToReduce.nrows())

    matProjected = matToReduce * matProjector

    ## Build the argument matrix for LLL in such a way that the transformation

    #  matrix is also returned. This matrix is obtained at almost no extra

    # cost. An identity matrix must be appended to 

    #  the left of the initial matrix. The transformation matrix will

    # will be recovered at the same location from the returned matrix .

    idMat = identity_matrix(matProjected.nrows())

    augmentedMatToReduce = idMat.augment(matProjected)

    reducedProjMat = \

        augmentedMatToReduce.LLL(algorithm='fpLLL:wrapper')

    ## Recover the transformation matrix (the left part of the reduced matrix). 

    #  We discard the reduced matrix itself.

    transMat = reducedProjMat.submatrix(0,

                                        0,

                                        reducedProjMat.nrows(),

                                        reducedProjMat.nrows())

    ## Return the initial matrix "reduced" and the transformation matrix tuple.

    return (transMat * matToReduce, transMat)

# End  slz_reduce_lll_proj.

#

def slz_uniform(n):

    """

    Compute a uniform RV in {-1, 1}.

    """

    ## function getrandbits(n) generates a long int with n random bits.

    return getrandbits(n) - 2^(n-1)

# End slz_uniform.

#

#! /opt/sage/sage

# @file runSLZ-113proj.sage

#

#@par Changes from runSLZ-113.sage

# LLL reduction is not performed on the matrix itself but rather on the

# product of the matrix with a uniform random matrix.

# The reduced matrix obtained is discarded but the transformation matrix 

# obtained is used to multiply the original matrix in order to reduced it.

# If a sufficient level of reduction is obtained, we stop here. If not

# the product matrix obtained above is LLL reduced. But as it has been

# pre-reduced at the above step, reduction is supposed to be much faster.

#

# Both reductions combined should hopefully be faster than a straight single

# reduction.

#

# Run SLZ for p=113

#from scipy.constants.codata import precision

def initialize_env():

    """

    Load all necessary modules.

    """

    compiledSpyxDir = "/home/storres/recherche/arithmetique/pobysoPythonSage/compiledSpyx"

    if compiledSpyxDir not in sys.path:

        sys.path.append(compiledSpyxDir)

    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/sageGMP.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()

from sageMpfr import *

from sageGMP  import *

import sys

from subprocess import call

#

## Main variables and parameters.

x         = var('x')

func(x)   = exp(x)

precision = 113

emin      = -16382

emax      = 16383 

RRR = RealField(precision)

degree              = 0

alpha               = 0

htrn                = 0

intervalCenter      = 0

intervalRadius      = 0

debugMode           = False          

## Local functions

#

def usage():

    write = sys.stderr.write

    write("\nUsage:\n")

    write("  " + scriptName + " <degree> <alpha> <htrn> <intervalCenter>\n")

    write("               <numberOfNumbers> [debug]\n")

    write("\nArguments:\n")

    write("  degree          the degree of the polynomial (integer)\n")

    write("  alpha           alpha (integer)\n")

    write("  htrn            hardness-to-round - a number of bits (integer)\n")

    write("  intervalCenter  the interval center (a floating-point number)\n")

    write("  numberOfNumbers the number of floating-point numbers in the interval\n")

    write("                  as a positive integral expression\n")

    write("  debug           debug mode (\"debug\", in any case)\n\n")

    sys.exit(2)

# End usage.

#

argsCount = len(sys.argv)

scriptName = os.path.basename(__file__)

if argsCount < 5:

    usage()

for index in xrange(1,argsCount):

    if index == 1:

        degree = int(sys.argv[index])

    elif index == 2:

        alpha = int(sys.argv[index])

    elif index == 3:

        htrn = int(eval(sys.argv[index]))

    elif index == 4:

        try:

            intervalCenter = QQ(sage_eval(sys.argv[index]))

        except:

            intervalCenter = RRR(sys.argv[index])

        intervalCenter = RRR(intervalCenter)

    elif index == 5:

        ## Can be read as rational number but must end up as an integer.

        numberOfNumbers = QQ(sage_eval(sys.argv[index]))

        if numberOfNumbers != numberOfNumbers.round():

            raise Exception("Invalid number of numbers: " + sys.argv[index] + ".")

        numberOfNumbers = numberOfNumbers.round()

        ## The number must be strictly positive.

        if numberOfNumbers <= 0:

            raise Exception("Invalid number of numbers: " + sys.argv[index] + ".")

    elif index == 6:

        debugMode = sys.argv[index].upper()

        debugMode = (debugMode == "DEBUG")

# Done with command line arguments collection.

#

## Debug printing

print "degree         :", degree

print "alpha          :", alpha

print "htrn           :", htrn

print "interval center:", intervalCenter.n(prec=10).str(truncate=False)

print "num of nums    :", RR(numberOfNumbers).log2().n(prec=10).str(truncate=False)

print "debug mode     :", debugMode

print

#

## Set the terminal window title.

terminalWindowTitle = ['stt', str(degree), str(alpha), str(htrn), 

                       intervalCenter.n(prec=10).str(truncate=False), 

                       RRR(numberOfNumbers).log2().n(prec=10).str(truncate=False)]

call(terminalWindowTitle)

#

intervalCenterBinade = slz_compute_binade(intervalCenter)

intervalRadius       = \

    2^intervalCenterBinade * 2^(-precision + 1) * numberOfNumbers / 2

srs_run_SLZ_v05_gram(inputFunction=func, 

                     inputLowerBound         = intervalCenter - intervalRadius, 

                     inputUpperBound         = intervalCenter + intervalRadius, 

                     alpha                   = alpha, 

                     degree                  = degree, 

                     precision               = precision, 

                     emin                    = emin, 

                     emax                    = emax, 

                     targetHardnessToRound   = htrn, 

                     debug                   = debugMode)

# This file was *autogenerated* from the file ./runSLZ-113gram.sage

from sage.all_cmdline import *   # import sage library

_sage_const_3 = Integer(3); _sage_const_2 = Integer(2); _sage_const_1 = Integer(1); _sage_const_0 = Integer(0); _sage_const_6 = Integer(6); _sage_const_5 = Integer(5); _sage_const_4 = Integer(4); _sage_const_113 = Integer(113); _sage_const_10 = Integer(10); _sage_const_16383 = Integer(16383); _sage_const_16382 = Integer(16382)#! /opt/sage/sage

# @file runSLZ-113.sage

#

# Run SLZ for p=113

#from scipy.constants.codata import precision

def initialize_env():

    """

    Load all necessary modules.

    """

    compiledSpyxDir = "/home/storres/recherche/arithmetique/pobysoPythonSage/compiledSpyx"

    if compiledSpyxDir not in sys.path:

        sys.path.append(compiledSpyxDir)

    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/sageGMP.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()

from sageMpfr import *

from sageGMP  import *

import sys

from subprocess import call

#

## Main variables and parameters.

x         = var('x')

__tmp__=var("x"); func = symbolic_expression(exp(x)).function(x)

precision = _sage_const_113 

emin      = -_sage_const_16382 

emax      = _sage_const_16383  

RRR = RealField(precision)

degree              = _sage_const_0 

alpha               = _sage_const_0 

htrn                = _sage_const_0 

intervalCenter      = _sage_const_0 

intervalRadius      = _sage_const_0 

debugMode           = False          

## Local functions

#

def usage():

    write = sys.stderr.write

    write("\nUsage:\n")

    write("  " + scriptName + " <degree> <alpha> <htrn> <intervalCenter>\n")

    write("               <numberOfNumbers> [debug]\n")

    write("\nArguments:\n")

    write("  degree          the degree of the polynomial (integer)\n")

    write("  alpha           alpha (integer)\n")

    write("  htrn            hardness-to-round - a number of bits (integer)\n")

    write("  intervalCenter  the interval center (a floating-point number)\n")

    write("  numberOfNumbers the number of floating-point numbers in the interval\n")

    write("                  as a positive integral expression\n")

    write("  debug           debug mode (\"debug\", in any case)\n\n")

    sys.exit(_sage_const_2 )

# End usage.

#

argsCount = len(sys.argv)

scriptName = os.path.basename(__file__)

if argsCount < _sage_const_5 :

    usage()

for index in xrange(_sage_const_1 ,argsCount):

    if index == _sage_const_1 :

        degree = int(sys.argv[index])

    elif index == _sage_const_2 :

        alpha = int(sys.argv[index])

    elif index == _sage_const_3 :

        htrn = int(eval(sys.argv[index]))

    elif index == _sage_const_4 :

        try:

            intervalCenter = QQ(sage_eval(sys.argv[index]))

        except:

            intervalCenter = RRR(sys.argv[index])

        intervalCenter = RRR(intervalCenter)

    elif index == _sage_const_5 :

        ## Can be read as rational number but must end up as an integer.

        numberOfNumbers = QQ(sage_eval(sys.argv[index]))

        if numberOfNumbers != numberOfNumbers.round():

            raise Exception("Invalid number of numbers: " + sys.argv[index] + ".")

        numberOfNumbers = numberOfNumbers.round()

        ## The number must be strictly positive.

        if numberOfNumbers <= _sage_const_0 :

            raise Exception("Invalid number of numbers: " + sys.argv[index] + ".")

    elif index == _sage_const_6 :

        debugMode = sys.argv[index].upper()

        debugMode = (debugMode == "DEBUG")

# Done with command line arguments collection.

#

## Debug printing

print "degree         :", degree

print "alpha          :", alpha

print "htrn           :", htrn

print "interval center:", intervalCenter.n(prec=_sage_const_10 ).str(truncate=False)

print "num of nums    :", RR(numberOfNumbers).log2().n(prec=_sage_const_10 ).str(truncate=False)

print "debug mode     :", debugMode

print

#

## Set the terminal window title.

terminalWindowTitle = ['stt', str(degree), str(alpha), str(htrn), 

                       intervalCenter.n(prec=_sage_const_10 ).str(truncate=False), 

                       RRR(numberOfNumbers).log2().n(prec=_sage_const_10 ).str(truncate=False)]

call(terminalWindowTitle)

#

intervalCenterBinade = slz_compute_binade(intervalCenter)

intervalRadius       = \

    _sage_const_2 **intervalCenterBinade * _sage_const_2 **(-precision + _sage_const_1 ) * numberOfNumbers / _sage_const_2 

srs_run_SLZ_v05_gram(inputFunction=func, 

                     inputLowerBound         = intervalCenter - intervalRadius, 

                     inputUpperBound         = intervalCenter + intervalRadius, 

                     alpha                   = alpha, 

                     degree                  = degree, 

                     precision               = precision, 

                     emin                    = emin, 

                     emax                    = emax, 

                     targetHardnessToRound   = htrn, 

                     debug                   = debugMode)

def srs_run_SLZ_v05_proj(inputFunction,

                         inputLowerBound,

                         inputUpperBound,

                         alpha,

                         degree,

                         precision,

                         emin,

                         emax,

                         targetHardnessToRound,

                         debug = False):

    """

    changes from plain V5:

       LLL reduction is not performed on the matrix itself but rather on the

       product of the matrix with a uniform random matrix.

       The reduced matrix obtained is discarded but the transformation matrix 

       obtained is used to multiply the original matrix in order to reduced it.

       If a sufficient level of reduction is obtained, we stop here. If not

       the product matrix obtained above is LLL reduced. But as it has been

       pre-reduced at the above step, reduction is supposed to be much fastet.

       Both reductions combined should hopefully be faster than a straight

       single reduction.

    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

    #

    ## Set the variable name in Sollya.

    pobyso_name_free_variable_sa_so(str(function.variables()[0]))

    ## 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_', '')

    #print "Scaled bounds:", sdlb, sdub

    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)

        if debug:

            print "Approximation polynomial computed."

        if prceSo is None:

            raise Exception("Could not compute an approximation polynomial.")

        ### 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 for debug purpose. Otherwise

        #    slz_compute_coppersmith_reduced_polynomials does the job

        #    invisibly.

        if debug:

            matrixToReduce = slz_compute_initial_lattice_matrix(intIntP,

                                                                alpha,

                                                                N,

                                                                iBound,

                                                                tBound)

            maxNorm     = 0

            latticeSize = 0

            matrixFile = file('/tmp/matrixToReduce.txt', 'w')

            for row in matrixToReduce.rows():

                currentNorm = row.norm()

                if currentNorm > maxNorm:

                    maxNorm = currentNorm

                latticeSize += 1

                for elem in row:

                    matrixFile.write(elem.str(base=2) + ",")

                matrixFile.write("\n")

            #matrixFile.write(matrixToReduce.str(radix="2") + "\n")

            matrixFile.close()

            #### We use here binary length as defined in LLL princepts.

            binaryLength = latticeSize * log(maxNorm)

            print "Binary length:", binaryLength.n()

            raise Exception("Deliberate stop here.")

        # End if debug

        reductionTime                     = cputime()

        #### Compute the reduced polynomials.

        print "Starting reduction..."

        ccReducedPolynomialsList =  \

                slz_compute_coppersmith_reduced_polynomials_proj(intIntP, 

                                                                 alpha, 

                                                                 N, 

                                                                 iBound, 

                                                                 tBound)

        print "...reduction accomplished in", cputime(reductionTime), "s."

        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:

        intervalResultString = "[" + str(intervalResultsList[0][0]) +\

                               "," + str(intervalResultsList[0][1])  + "]"

        print intervalResultString,

        if len(intervalResultsList) > 1:

            rootsResultsList = intervalResultsList[1]

            specificRootResultIndex = 0

            for specificRootResultsList in rootsResultsList:

                if specificRootResultIndex == 0:

                    print "\t", specificRootResultsList[0],

                else:

                    print " " * len(intervalResultString), "\t", \

                          specificRootResultsList[0],

                if len(specificRootResultsList) > 1:

                    print specificRootResultsList[1]

                specificRootResultIndex += 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_proj

#