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

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

#load(str('/home/storres/recherche/arithmetique/pobysoPythonSage/src/sageSLZ/sageMatrixOperations.sage'))

print "sagePolynomialOperations loading..."

def spo_add_polynomial_coeffs_to_matrix_row(poly, 

                                            knownMonomials, 

                                            protoMatrixRows, 

                                            columnsWidth=0):

    """

    For a given polynomial ,

    add the coefficients of the protoMatrix (a list of proto matrix rows).

    Coefficients are added to the protoMatrix row in the order imposed by the 

    monomials discovery list (the knownMonomials list) built as construction

    goes on. 

    As a bonus, data can be printed out for a visual check.

    poly           : the polynomial; in argument;

    knownMonomials : the list of the already known monomials; will determine

                     the order of the coefficients appending to a row; in-out

                     argument (new monomials may be discovered and then 

                     appended the the knowMonomials list);

    protoMatrixRows: a list of lists, each one holding the coefficients of the 

                     monomials of a polynomial; in-out argument: a new row is

                     added at each call;

    columnWith     : the width, in characters, of the displayed column ; if 0, 

                     do not display anything; in argument.

    """

    pMonomials   = poly.monomials()

    pCoefficients = poly.coefficients()

    # We have started with the smaller degrees in the first variable.

    pMonomials.reverse()

    pCoefficients.reverse()

    # New empty proto matrix row.

    protoMatrixRowCoefficients = []

    # We work according to the order of the already known monomials

    # No known monomials yet: add the pMonomials to knownMonomials 

    # and add the coefficients to the proto matrix row.

    if len(knownMonomials) == 0: 

        for pmIdx in xrange(0, len(pMonomials)):

            knownMonomials.append(pMonomials[pmIdx])

            protoMatrixRowCoefficients.append(pCoefficients[pmIdx])

            if columnsWidth != 0:

                monomialAsString = str(pCoefficients[pmIdx]) + " " + \

                                   str(pMonomials[pmIdx])

                print monomialAsString, " " * \

                      (columnsWidth - len(monomialAsString)),

    # There are some known monomials. We search for them in pMonomials and 

    # add their coefficients to the proto matrix row.

    else: 

        for knownMonomialIndex in xrange(0,len(knownMonomials)):

            # We lazily use an exception here since pMonomials.index() function

            # may fail throwing the ValueError exception.

            try:

                indexInPmonomials = \

                    pMonomials.index(knownMonomials[knownMonomialIndex])

                if columnsWidth != 0:

                    monomialAsString = str(pCoefficients[indexInPmonomials]) + \

                        " " + str(knownMonomials[knownMonomialIndex])

                    print monomialAsString, " " * \

                        (columnsWidth - len(monomialAsString)),

                # Add the coefficient to the proto matrix row and delete the 

                # known monomial from the current pMonomial list 

                # (and the corresponding coefficient as well).

                protoMatrixRowCoefficients.append(pCoefficients[indexInPmonomials])

                del pMonomials[indexInPmonomials]

                del pCoefficients[indexInPmonomials]

            # The knownMonomials element is not in pMonomials

            except ValueError:   

                protoMatrixRowCoefficients.append(0)

                if columnsWidth != 0:

                    monomialAsString = "0" + " "+ \

                        str(knownMonomials[knownMonomialIndex])

                    print monomialAsString, " " * \

                        (columnsWidth - len(monomialAsString)),

        # End for knownMonomialKey loop. 

        # We now append the remaining monomials of pMonomials to knownMonomials 

        # and the corresponding coefficients to proto matrix row.

        for pmIdx in xrange(0, len(pMonomials)):

            knownMonomials.append(pMonomials[pmIdx])

            protoMatrixRowCoefficients.append(pCoefficients[pmIdx])

            if columnsWidth != 0:

                monomialAsString = str(pCoefficients[pmIdx]) + " " \

                    + str(pMonomials[pmIdx])

                print monomialAsString, " " * \

                    (columnsWidth - len(monomialAsString)),

        # End for pmIdx loop.

    # Add the new list row elements to the proto matrix.

    protoMatrixRows.append(protoMatrixRowCoefficients)

    if columnsWidth != 0:

        print    

# End spo_add_polynomial_coeffs_to_matrix_row

def spo_float_poly_of_float_to_rat_poly_of_rat(polyOfFloat):

    """

    Create a polynomial over the rationals from a polynomial over

    a RealField.

    Important warning: default Sage behavior is to convert coefficients

    using continued fractions instead of making a simple conversion

    with a powers of two at denominators (and possible simplification).

    Hence conversion is not exact but with a relative error around 10^(-521).

    """

    ratPolynomialRing = QQ[str(polyOfFloat.variables()[0])]

    return(ratPolynomialRing(polyOfFloat))

# End spo_float_poly_of_float_to_rat_poly_of_rat.

def spo_float_poly_of_float_to_rat_poly_of_rat_pow_two(polyOfFloat):

    """

    Create a polynomial over the rationals from a polynomial over

    a RealField where all denominators are

    powers of two.

    Allows for exact conversions (and lcm computation of the coefficients

    denominator).

    """

    polyVariable = polyOfFloat.variables()[0] 

    RPR = QQ[str(polyVariable)]

    polyCoeffs      = polyOfFloat.coefficients()

    #print polyCoeffs

    polyExponents   = polyOfFloat.exponents()

    #print polyExponents

    polyDenomPtwoCoeffs = []

    for coeff in polyCoeffs:

        polyDenomPtwoCoeffs.append(sno_float_to_rat_pow_of_two_denom(coeff))

        #print "Converted coefficient:", sno_float_to_rat_pow_of_two_denom(coeff),

        #print type(sno_float_to_rat_pow_of_two_denom(coeff))

    ratPoly = RPR(0)

    #print type(ratPoly)

    ## !!! CAUTION !!! Do not use the RPR(coeff * polyVariagle^exponent)

    #  construction.

    #  The coefficient becomes plainly wrong when exponent == 0.

    #  No clue as to why.

    for coeff, exponent in zip(polyDenomPtwoCoeffs, polyExponents):

        ratPoly += coeff * RPR(polyVariable^exponent)

    return ratPoly

# End slz_float_poly_of_float_to_rat_poly_of_rat.

def spo_get_coefficient_for_monomial(monomialsList, coefficientsList, monomial):

    """

    Get, for a polynomial, the coefficient for a given monomial.

    The polynomial is given as two lists (monomials and coefficients as

    return by the respective methods ; indexes of the two lists must match).

    If the monomial is not found, 0 is returned.

    """

    monomialIndex = 0

    for mono in monomialsList:

        if mono == monomial:

            return coefficientsList[monomialIndex]

        monomialIndex += 1

    return 0

# End spo_get_coefficient_for_monomial.

def spo_expression_as_string(powI, boundI, powT, boundT, powP, powN):

    """

    Computes a string version of the i^k + t^l + p^m + N^n expression for

    output.

    """

    expressionAsString =""

    if powI != 0:

        expressionAsString += str(iBound^powI) + " i^" + str(powI)

    if powT != 0:

        if len(expressionAsString) != 0:

            expressionAsString += " * "

        expressionAsString += str(tBound^powT) + " t^" + str(powT)

    if powP != 0:

        if len(expressionAsString) != 0:

            expressionAsString += " * "

        expressionAsString += "p^" + str(powP)

    if (powN) != 0 :

        if len(expressionAsString) != 0:

            expressionAsString += " * "

        expressionAsString += "N^" + str(powN)

    return(expressionAsString)

# End spo_expression_as_string.

def spo_norm(poly, p=2):

    """

    Behaves more or less (no infinity defined) as the norm for the

    univariate polynomials.

    Quoting Sage documentation:

    "Definition: For integer p, the p-norm of a polynomial is the pth root of 

    the sum of the pth powers of the absolute values of the coefficients of 

    the polynomial."

    """

    # TODO: check the arguments (for p see below)..

    norm = 0

    # For infinity norm.

    if p == Infinity:

        for coefficient in poly.coefficients():

            coefficientAbs = coefficient.abs()

            if coefficientAbs > norm:

                norm = coefficientAbs

        return norm

    # TODO: check here the value of p

    # p must be a positive integer >= 1.

    if p < 1 or (not p in ZZ):

        return None

    # For 1 norm.

    if p == 1:

        for coefficient in poly.coefficients():

            norm +=  coefficient.abs()

        return norm

    # For other norms

    for coefficient in poly.coefficients():

        norm +=  coefficient.abs()^p

    return pow(norm, 1/p)

# end spo_norm

def spo_polynomial_to_proto_matrix(p, alpha, N, columnsWidth=0):

    """

    From a (bivariate) polynomial and some other parameters build a proto 

    matrix (an array of "rows") to be converted into a "true" matrix and 

    eventually by reduced by fpLLL.

    The matrix is such as those found in Boneh-Durphee and Stehlé.

    Parameters

    ----------

    p: the (bivariate) polynomial;

    pRing: the ring over which p is defined;

    alpha:

    N:

    columsWidth: if == 0, no information is displayed, otherwise data is 

                 printed in colums of columnsWitdth width.

    """

    pRing = p.parent()

    knownMonomials = []

    protoMatrixRows = []

    polynomialsList = []

    pVariables = p.variables()

    #print "In spo...", p, p.variables()

    iVariable = pVariables[0]

    tVariable = pVariables[1]

    polynomialAtPower = pRing(1)

    currentPolynomial = pRing(1)

    pIdegree = p.degree(pVariables[0])

    pTdegree = p.degree(pVariables[1])

    currentIdegree = currentPolynomial.degree(iVariable)

    nAtAlpha = N^alpha

    nAtPower = nAtAlpha

    polExpStr = ""

    # We work from p^0 * N^alpha to p^alpha * N^0

    for pPower in xrange(0, alpha + 1):

        # pPower == 0 is a special case. We introduce all the monomials but one

        # in i and those in t necessary to be able to introduce

        # p. We arbitrary choose to introduce the highest degree monomial in i

        # with p. We also introduce all the mixed i^k * t^l monomials with

        # k < p.degree(i) and l <= p.degree(t).

        # Mixed terms introduction is necessary here before we start "i shifts"

        # in the next iteration.

        if pPower == 0:

            # Notice that i^pIdegree is excluded as the bound of the xrange is

            # pIdegree

            for iPower in xrange(0, pIdegree): 

                for tPower in xrange(0, pTdegree + 1):

                    if columnsWidth != 0:

                        polExpStr = spo_expression_as_string(iPower, 

                                                             tPower, 

                                                             pPower, 

                                                             alpha-pPower)

                        print "->", polExpStr

                    currentExpression = iVariable^iPower * \

                                        tVariable^tPower * nAtAlpha

                    # polynomialAtPower == 1 here. Next line should be commented

                    # out but it does not work! Some conversion problem?

                    currentPolynomial = pRing(currentExpression)

                    polynomialsList.append(currentPolynomial) 

                    pMonomials = currentPolynomial.monomials()

                    pCoefficients = currentPolynomial.coefficients()

                    spo_add_polynomial_coeffs_to_matrix_row(pMonomials, 

                                                            pCoefficients, 

                                                            knownMonomials, 

                                                            protoMatrixRows, 

                                                            columnsWidth)

                # End tPower.

            # End for iPower.

        else: # pPower > 0: (p^1..p^alpha)

            # This where we introduce the p^pPower * N^(alpha-pPower)

            # polynomial.

            # This step could technically be fused as the first iteration

            # of the next loop (with iPower starting at 0).

            # We set it apart for clarity.

            if columnsWidth != 0:

                polExpStr = spo_expression_as_string(0, 0, pPower, alpha-pPower)

                print "->", polExpStr

            currentPolynomial = polynomialAtPower * nAtPower

            polynomialsList.append(currentPolynomial) 

            pMonomials = currentPolynomial.monomials()

            pCoefficients = currentPolynomial.coefficients()

            spo_add_polynomial_coeffs_to_matrix_row(pMonomials, 

                                                    pCoefficients, 

                                                    knownMonomials, 

                                                    protoMatrixRows, 

                                                    columnsWidth)

            # The i^iPower * p^pPower polynomials: they add i^k monomials to  

            # p^pPower up to k < pIdegree * pPower. This only introduces i^k 

            # monomials since mixed terms (that were introduced at a previous

            # stage) are only shifted to already existing 

            # ones. p^pPower is "shifted" to higher degrees in i as far as 

            # possible, one step short of the degree in i of p^(pPower+1) .

            # These "pure" i^k monomials can only show up with i multiplications.

            for iPower in xrange(1, pIdegree):

                if columnsWidth != 0:

                    polExpStr = spo_expression_as_string(iPower, \

                                                         0,      \

                                                         pPower, \

                                                         alpha)

                    print "->", polExpStr

                currentExpression = i^iPower * nAtPower

                currentPolynomial = pRing(currentExpression) * polynomialAtPower

                polynomialsList.append(currentPolynomial) 

                pMonomials = currentPolynomial.monomials()

                pCoefficients = currentPolynomial.coefficients()

                spo_add_polynomial_coeffs_to_matrix_row(pMonomials,      \

                                                        pCoefficients,   \

                                                        knownMonomials,  \

                                                        protoMatrixRows, \

                                                        columnsWidth)

            # End for iPower

            # We want now to introduce a t * p^pPower polynomial. But before

            # that we must introduce some mixed monomials.

            # This loop is no triggered before pPower == 2.

            # It introduces a first set of high i degree mixed monomials.

            for iPower in xrange(1, pPower):

                tPower = pPower - iPower + 1

                if columnsWidth != 0:

                    polExpStr = spo_expression_as_string(iPower * pIdegree, 

                                                         tPower, 

                                                         0,  

                                                         alpha)

                    print "->", polExpStr

                currentExpression = i^(iPower * pIdegree) * t^tPower * nAtAlpha

                currentPolynomial = pRing(currentExpression)

                polynomialsList.append(currentPolynomial) 

                pMonomials = currentPolynomial.monomials()

                pCoefficients = currentPolynomial.coefficients()

                spo_add_polynomial_coeffs_to_matrix_row(pMonomials, 

                                                        pCoefficients, 

                                                        knownMonomials, 

                                                        protoMatrixRows, 

                                                        columnsWidth)

            # End for iPower

            #

            # This is the mixed monomials main loop. It introduces:

            # - the missing mixed monomials needed before the 

            #   t^l * p^pPower * N^(alpha-pPower) polynomial;

            # - the t^l * p^pPower * N^(alpha-pPower) itself;

            # - for each of i^k * t^l * p^pPower * N^(alpha-pPower) polynomials:

            #   - the the missing mixed monomials needed  polynomials,

            #   - the i^k * t^l * p^pPower * N^(alpha-pPower) itself.

            # The t^l * p^pPower * N^(alpha-pPower) is introduced when

            # 

            for iShift in xrange(0, pIdegree):

                # When pTdegree == 1, the following loop only introduces 

                # a single new monomial.

                #print "++++++++++"

                for outerTpower in xrange(1, pTdegree + 1):

                    # First one high i degree mixed monomial.

                    iPower = iShift + pPower * pIdegree

                    if columnsWidth != 0:

                        polExpStr = spo_expression_as_string(iPower, 

                                                             outerTpower,

                                                             0,

                                                             alpha)

                        print "->", polExpStr

                    currentExpression = i^iPower * t^outerTpower * nAtAlpha

                    currentPolynomial = pRing(currentExpression)

                    polynomialsList.append(currentPolynomial) 

                    pMonomials = currentPolynomial.monomials()

                    pCoefficients = currentPolynomial.coefficients()

                    spo_add_polynomial_coeffs_to_matrix_row(pMonomials, 

                                                            pCoefficients, 

                                                            knownMonomials, 

                                                            protoMatrixRows, 

                                                            columnsWidth)

                    #print "+++++"

                    # At iShift == 0, the following innerTpower loop adds  

                    # duplicate monomials, since no extra i^l * t^k is needed 

                    # before introducing the  

                    # i^iShift * t^outerPpower * p^pPower * N^(alpha-pPower) 

                    # polynomial.

                    # It introduces smaller i degree monomials than the

                    # one(s) added previously (no pPower multiplication).

                    # Here the exponent of t decreases as that of i increases.

                    # This conditional is not entered before pPower == 1.

                    # The innerTpower loop does not produce anything before

                    # pPower == 2. We keep it anyway for other configuration of

                    # p.

                    if iShift > 0:

                        iPower = pIdegree + iShift

                        for innerTpower in xrange(pPower, 1, -1):

                            if columnsWidth != 0:

                                polExpStr = spo_expression_as_string(iPower, 

                                                                     innerTpower,

                                                                     0,

                                                                     alpha) 

                            currentExpression = \

                                    i^(iPower) * t^(innerTpower) * nAtAlpha

                            currentPolynomial = pRing(currentExpression)

                            polynomialsList.append(currentPolynomial) 

                            pMonomials = currentPolynomial.monomials()

                            pCoefficients = currentPolynomial.coefficients()

                            spo_add_polynomial_coeffs_to_matrix_row(pMonomials, 

                                                                pCoefficients, 

                                                                knownMonomials, 

                                                                protoMatrixRows, 

                                                                columnsWidth)

                            iPower += pIdegree

                        # End for innerTpower

                    # End of if iShift > 0

                    # When iShift == 0, just after each of the  

                    # p^pPower * N^(alpha-pPower) polynomials has 

                    # been introduced (followed by a string of 

                    # i^k * p^pPower * N^(alpha-pPower) polynomials) a

                    # t^l *  p^pPower * N^(alpha-pPower) is introduced here.

                    # 

                    # Eventually, the following section introduces the 

                    # i^iShift * t^outerTpower * p^iPower * N^(alpha-pPower) 

                    # polynomials.

                    if columnsWidth != 0:

                        polExpStr = spo_expression_as_string(iShift, 

                                                             outerTpower,

                                                             pPower,

                                                             alpha-pPower)

                        print "->", polExpStr

                    currentExpression = i^iShift * t^outerTpower * nAtPower

                    currentPolynomial = pRing(currentExpression) * \

                                            polynomialAtPower

                    polynomialsList.append(currentPolynomial) 

                    pMonomials = currentPolynomial.monomials()

                    pCoefficients = currentPolynomial.coefficients()

                    spo_add_polynomial_coeffs_to_matrix_row(pMonomials, 

                                                            pCoefficients, 

                                                            knownMonomials, 

                                                            protoMatrixRows, 

                                                            columnsWidth)

                # End for outerTpower 

                #print "++++++++++"

            # End for iShift

        polynomialAtPower *= p  

        nAtPower /= N

    # End for pPower loop

    return ((protoMatrixRows, knownMonomials, polynomialsList))

# End spo_polynomial_to_proto_matrix

def spo_polynomial_to_polynomials_list_2(p, alpha, N, iBound, tBound,

                                         columnsWidth=0):

    """

    Badly out of sync code: check with versions 3 or 4.

    From p, alpha, N build a list of polynomials...

    TODO: clean up the comments below!

    From a (bivariate) polynomial and some other parameters build a proto 

    matrix (an array of "rows") to be converted into a "true" matrix and 

    eventually by reduced by fpLLL.

    The matrix is based on a list of polynomials that are built in a way

    that one and only monomial is added at each new polynomial. Among the many

    possible ways to build this list we pick one strongly dependent on the 

    structure of the polynomial and of the problem.

    We consider here the polynomials of the form: 

    a_k*i^k + a_(k-1)*i^(k-1) + ... + a_1*i + a_0 - t 

    The values of i and t are bounded and we eventually look for (i_0,t_0) 

    pairs such that:

    a_k*i_0^k + a_(k-1)*i_0^(k-1) + ... + a_1*i_0 + a_0 = t_0

    Hence, departing from the procedure in described in Boneh-Durfee, we will 

    not use "t-shifts" but only "i-shifts".

    Parameters

    ----------

    p: the (bivariate) polynomial;

    pRing: the ring over which p is defined;

    alpha:

    N:

    columsWidth: if == 0, no information is displayed, otherwise data is 

                 printed in colums of columnsWitdth width.

    """

    pRing = p.parent()

    polynomialsList = []

    pVariables = p.variables()

    iVariable = pVariables[0]

    tVariable = pVariables[1]

    polynomialAtPower = pRing(1)

    currentPolynomial = pRing(1)

    pIdegree = p.degree(iVariable)

    pTdegree = p.degree(tVariable)

    currentIdegree = currentPolynomial.degree(iVariable)

    nAtAlpha = N^alpha

    nAtPower = nAtAlpha

    polExpStr = ""

    # We work from p^0 * N^alpha to p^alpha * N^0

    for pPower in xrange(0, alpha + 1):

        # pPower == 0 is a special case. We introduce all the monomials in i 

        # up to i^pIdegree.

        if pPower == 0:

            # Notice who iPower runs up to i^pIdegree.

            for iPower in xrange(0, pIdegree + 1): 

                # No t power is taken into account as we limit our selves to

                # degree 1 in t and make no "t-shifts".

                if columnsWidth != 0:

                    polExpStr = spo_expression_as_string(iPower,

                                                         iBound, 

                                                         0, 

                                                         tBound,

                                                         0, 

                                                         alpha)

                    print "->", polExpStr

                currentExpression = iVariable^iPower * nAtAlpha * iBound^iPower

                # polynomialAtPower == 1 here. Next line should be commented

                # out but it does not work! Some conversion problem?

                currentPolynomial = pRing(currentExpression)

                polynomialsList.append(currentPolynomial) 

            # End for iPower.

        else: # pPower > 0: (p^1..p^alpha)

            # This where we introduce the p^pPower * N^(alpha-pPower)

            # polynomial. This is also where the t^pPower monomials shows up for

            # the first time.

            if columnsWidth != 0:

                polExpStr = spo_expression_as_string(0, iBound, 0, tBound, \

                                                     pPower, alpha-pPower)

                print "->", polExpStr

            currentPolynomial = polynomialAtPower * nAtPower

            polynomialsList.append(currentPolynomial) 

            # Exit when pPower == alpha

            if pPower == alpha:

                return polynomialsList

            # This is where the "i-shifts" take place. Mixed terms, i^k * t^l

            # (that were introduced at a previous

            # stage or are introduced now) are only shifted to already existing 

            # ones with the notable exception of i^iPower * t^pPower, which

            # must be manually introduced.

            # p^pPower is "shifted" to higher degrees in i as far as 

            # possible, up to of the degree in i of p^(pPower+1).

            # These "pure" i^k monomials can only show up with i multiplications.

            for iPower in xrange(1, pIdegree + 1):

                # The i^iPower * t^pPower monomial. Notice the alpha exponent

                # for N.

                internalIpower = iPower

                for tPower in xrange(pPower,0,-1):

                    if columnsWidth != 0:

                        polExpStr = spo_expression_as_string(internalIpower, 

                                                             iBound,

                                                             tPower,

                                                             tBound,

                                                             0,

                                                             alpha)

                        print "->", polExpStr

                    currentExpression = i^internalIpower * t^tPower * \

                                        nAtAlpha * iBound^internalIpower * \

                                        tBound^tPower

                    currentPolynomial = pRing(currentExpression)

                    polynomialsList.append(currentPolynomial) 

                    internalIpower += pIdegree

                # End for tPower

                # The i^iPower * p^pPower * N^(alpha-pPower) i-shift.

                if columnsWidth != 0:

                    polExpStr = spo_expression_as_string(iPower, 

                                                         iBound,

                                                         0,      

                                                         tBound,

                                                         pPower, 

                                                         alpha-pPower)

                    print "->", polExpStr

                currentExpression = i^iPower * nAtPower * iBound^iPower

                currentPolynomial = pRing(currentExpression) * polynomialAtPower

                polynomialsList.append(currentPolynomial) 

            # End for iPower

        polynomialAtPower *= p  

        nAtPower /= N

    # End for pPower loop

    return polynomialsList

# End spo_polynomial_to_proto_matrix_2

def spo_polynomial_to_polynomials_list_3(p, alpha, N, iBound, tBound,

                                         columnsWidth=0):

    """

    From p, alpha, N build a list of polynomials...

    TODO: more in depth rationale...

    Our goal is to introduce each monomial with the smallest coefficient.

    Parameters

    ----------

    p: the (bivariate) polynomial;

    pRing: the ring over which p is defined;

    alpha:

    N:

    columsWidth: if == 0, no information is displayed, otherwise data is 

                 printed in colums of columnsWitdth width.

    """

    pRing = p.parent()

    polynomialsList = []

    pVariables = p.variables()

    iVariable = pVariables[0]

    tVariable = pVariables[1]

    polynomialAtPower = pRing(1)

    currentPolynomial = pRing(1)

    pIdegree = p.degree(iVariable)

    pTdegree = p.degree(tVariable)

    currentIdegree = currentPolynomial.degree(iVariable)

    nAtAlpha = N^alpha

    nAtPower = nAtAlpha

    polExpStr = ""

    # We work from p^0 * N^alpha to p^alpha * N^0

    for pPower in xrange(0, alpha + 1):

        # pPower == 0 is a special case. We introduce all the monomials in i 

        # up to i^pIdegree.

        if pPower == 0:

            # Notice who iPower runs up to i^pIdegree.

            for iPower in xrange(0, pIdegree + 1): 

                # No t power is taken into account as we limit our selves to

                # degree 1 in t and make no "t-shifts".

                if columnsWidth != 0:

                    polExpStr = spo_expression_as_string(iPower, 

                                                         iBound,

                                                         0, 

                                                         tBound,

                                                         0, 

                                                         alpha)

                    print "->", polExpStr

                currentExpression = iVariable^iPower * nAtAlpha * iBound^iPower

                # polynomialAtPower == 1 here. Next line should be commented

                # out but it does not work! Some conversion problem?

                currentPolynomial = pRing(currentExpression)

                polynomialsList.append(currentPolynomial) 

            # End for iPower.

        else: # pPower > 0: (p^1..p^alpha)

            # This where we introduce the p^pPower * N^(alpha-pPower)

            # polynomial. This is also where the t^pPower monomials shows up for

            # the first time. It app

            if columnsWidth != 0:

                polExpStr = spo_expression_as_string(0, iBound, 

                                                     0, tBound,

                                                     pPower, alpha-pPower)

                print "->", polExpStr

            currentPolynomial = polynomialAtPower * nAtPower

            polynomialsList.append(currentPolynomial) 

            # Exit when pPower == alpha

            if pPower == alpha:

                return polynomialsList

            # This is where the "i-shifts" take place. Mixed terms, i^k * t^l

            # (that were introduced at a previous

            # stage or are introduced now) are only shifted to already existing 

            # ones with the notable exception of i^iPower * t^pPower, which

            # must be manually introduced.

            # p^pPower is "shifted" to higher degrees in i as far as 

            # possible, up to of the degree in i of p^(pPower+1).

            # These "pure" i^k monomials can only show up with i multiplications.

            for iPower in xrange(1, pIdegree + 1):

                # The i^iPower * t^pPower monomial. Notice the alpha exponent

                # for N.

                internalIpower = iPower

                for tPower in xrange(pPower,0,-1):

                    if columnsWidth != 0:

                        polExpStr = spo_expression_as_string(internalIpower,

                                                             iBound, 

                                                             tPower,  

                                                             tBound,

                                                             0,

                                                             alpha)

                        print "->", polExpStr

                    currentExpression = i^internalIpower * t^tPower * nAtAlpha * \

                                        iBound^internalIpower * tBound^tPower

                    currentPolynomial = pRing(currentExpression)

                    polynomialsList.append(currentPolynomial) 

                    internalIpower += pIdegree

                # End for tPower

                # Here we have to choose between a 

                # i^iPower * p^pPower * N^(alpha-pPower) i-shift and

                # i^iPower * i^(d_i(p) * pPower) * N^alpha, depending on which 

                # coefficient is smallest.

                IcurrentExponent = iPower + \

                                (pPower * polynomialAtPower.degree(iVariable))

                currentMonomial = pRing(iVariable^IcurrentExponent)

                currentPolynomial = pRing(iVariable^iPower * nAtPower * \

                                          iBound^iPower) * \

                                          polynomialAtPower

                currMonomials = currentPolynomial.monomials()

                currCoefficients = currentPolynomial.coefficients()

                currentCoefficient = spo_get_coefficient_for_monomial( \

                                                    currMonomials, 

                                                    currCoefficients, 

                                                    currentMonomial)

                print "Current coefficient:", currentCoefficient

                alterCoefficient = iBound^IcurrentExponent * nAtAlpha 

                print "N^alpha * ibound^", IcurrentExponent, ":", \

                         alterCoefficient

                if currentCoefficient > alterCoefficient :

                    if columnsWidth != 0:

                        polExpStr = spo_expression_as_string(IcurrentExponent,

                                                             iBound,

                                                             0,

                                                             tBound,

                                                             0,

                                                             alpha)

                        print "->", polExpStr

                    polynomialsList.append(currentMonomial * \

                                           alterCoefficient)

                else:

                    if columnsWidth != 0:

                        polExpStr = spo_expression_as_string(iPower, iBound,

                                                             0, tBound,

                                                             pPower, 

                                                             alpha-pPower)

                        print "->", polExpStr

                    polynomialsList.append(currentPolynomial) 

            # End for iPower

        polynomialAtPower *= p  

        nAtPower /= N

    # End for pPower loop

    return polynomialsList 

# End spo_polynomial_to_proto_matrix_3

def spo_polynomial_to_polynomials_list_4(p, alpha, N, iBound, tBound,

                                         columnsWidth=0):

    """

    From p, alpha, N build a list of polynomials...

    TODO: more in depth rationale...

    Our goal is to introduce each monomial with the smallest coefficient.

    Parameters

    ----------

    p: the (bivariate) polynomial;

    pRing: the ring over which p is defined;

    alpha:

    N:

    columsWidth: if == 0, no information is displayed, otherwise data is 

                 printed in colums of columnsWitdth width.

    """

    pRing = p.parent()

    polynomialsList = []

    pVariables = p.variables()

    iVariable = pVariables[0]

    tVariable = pVariables[1]

    polynomialAtPower = copy(p)

    currentPolynomial = pRing(1)

    pIdegree = p.degree(iVariable)

    pTdegree = p.degree(tVariable)

    maxIdegree = pIdegree * alpha

    currentIdegree = currentPolynomial.degree(iVariable)

    nAtAlpha = N^alpha

    nAtPower = nAtAlpha

    polExpStr = ""

    # We first introduce all the monomials in i alone multiplied by N^alpha.

    for iPower in xrange(0, maxIdegree + 1):

        if columnsWidth !=0:

            polExpStr = spo_expression_as_string(iPower, iBound,

                                                 0, tBound, 

                                                 0, alpha)

            print "->", polExpStr

        currentExpression = iVariable^iPower * nAtAlpha * iBound^iPower

        currentPolynomial = pRing(currentExpression)

        polynomialsList.append(currentPolynomial)

    # End for iPower

    # We work from p^1 * N^alpha-1 to p^alpha * N^0

    for pPower in xrange(1, alpha + 1):

        # First of all the p^pPower * N^(alpha-pPower) polynomial.

        nAtPower /= N

        if columnsWidth !=0:

            polExpStr = spo_expression_as_string(0, iBound,

                                                 0, tBound,

                                                 pPower, alpha-pPower)

            print "->", polExpStr

        currentPolynomial = polynomialAtPower * nAtPower

        polynomialsList.append(currentPolynomial)

        # Exit when pPower == alpha

        if pPower == alpha:

            return polynomialsList

        # We now introduce the mixed i^k * t^l monomials by i^m * p^n * N^(alpha-n)

        for iPower in xrange(1, pIdegree + 1): 

            if columnsWidth != 0:

                polExpStr = spo_expression_as_string(iPower, iBound,

                                                     0, tBound, 

                                                     pPower, alpha-pPower)

                print "->", polExpStr

            currentExpression = i^iPower * iBound^iPower * nAtPower

            currentPolynomial = pRing(currentExpression) * polynomialAtPower

            polynomialsList.append(currentPolynomial) 

        # End for iPower

        polynomialAtPower *= p  

    # End for pPower loop

    return polynomialsList 

# End spo_polynomial_to_proto_matrix_4

def spo_polynomial_to_polynomials_list_5(p, alpha, N, iBound, tBound,

                                         columnsWidth=0):

    """

    From p, alpha, N build a list of polynomials use to create a base

    that will eventually be reduced with LLL.

    The bounds are computed for the coefficients that will be used to

    form the base.

    We try to introduce only one new monomial at a time, to obtain a

    triangular matrix (it is easy to compute the volume of the underlining

    latice if the matrix is triangular).

    There are many possibilities to introduce the monomials: our goal is also 

    to introduce each of them on the diagonal with the smallest coefficient.

    The method depends on the structure of the polynomial. Here it is adapted

    to the a_n*i^n + ... + a_1 * i - t + b form.     

    Parameters

    ----------

    p: the (bivariate) polynomial;

    alpha:

    N:

    iBound:

    tBound:

    columsWidth: if == 0, no information is displayed, otherwise data is 

                 printed in colums of columnsWitdth width.

    """

    pRing = p.parent()

    polynomialsList = []

    pVariables = p.variables()

    iVariable = pVariables[0]

    tVariable = pVariables[1]

    polynomialAtPower = copy(p)

    currentPolynomial = pRing(1)

    pIdegree = p.degree(iVariable)

    pTdegree = p.degree(tVariable)

    maxIdegree = pIdegree * alpha

    currentIdegree = currentPolynomial.degree(iVariable)

    nAtAlpha = N^alpha

    nAtPower = nAtAlpha

    polExpStr = ""

    # We first introduce all the monomials in i alone multiplied by N^alpha.

    for iPower in xrange(0, maxIdegree + 1):

        if columnsWidth !=0:

            polExpStr = spo_expression_as_string(iPower, iBound,

                                                 0, tBound, 

                                                 0, alpha)

            print "->", polExpStr

        currentExpression = iVariable^iPower * nAtAlpha * iBound^iPower

        currentPolynomial = pRing(currentExpression)

        polynomialsList.append(currentPolynomial)

    # End for iPower

    # We work from p^1 * N^alpha-1 to p^alpha * N^0

    for pPower in xrange(1, alpha + 1):

        # First of all the p^pPower * N^(alpha-pPower) polynomial.

        nAtPower /= N

        if columnsWidth !=0:

            polExpStr = spo_expression_as_string(0, iBound,

                                                 0, tBound,

                                                 pPower, alpha-pPower)

            print "->", polExpStr

        currentPolynomial = polynomialAtPower * nAtPower

        polynomialsList.append(currentPolynomial)

        # Exit when pPower == alpha

        if pPower == alpha:

            return polynomialsList

        for iPower in xrange(1, pIdegree + 1):

            iCurrentPower = pIdegree + iPower

            for tPower in xrange(pPower-1, 0, -1):

                #print "tPower:", tPower

                if columnsWidth != 0:

                    polExpStr = spo_expression_as_string(iCurrentPower, iBound,

                                                         tPower, tBound, 

                                                         0, alpha)

                    print "->", polExpStr

                currentExpression = i^iCurrentPower * iBound^iCurrentPower * t^tPower * tBound^tPower *nAtAlpha 

                currentPolynomial = pRing(currentExpression)

                polynomialsList.append(currentPolynomial)

                iCurrentPower += pIdegree 

            # End for tPower 

        # We now introduce the mixed i^k * t^l monomials by i^m * p^n * N^(alpha-n)

            if columnsWidth != 0:

                polExpStr = spo_expression_as_string(iPower, iBound,

                                                     0, tBound, 

                                                     pPower, alpha-pPower)

                print "->", polExpStr

            currentExpression = i^iPower * iBound^iPower * nAtPower

            currentPolynomial = pRing(currentExpression) * polynomialAtPower

            polynomialsList.append(currentPolynomial) 

        # End for iPower

        polynomialAtPower *= p  

    # End for pPower loop

    return polynomialsList 

# End spo_polynomial_to_proto_matrix_5

def spo_polynomial_to_polynomials_list_6(p, alpha, N, iBound, tBound,

                                         columnsWidth=0):

    """

    From p, alpha, N build a list of polynomials use to create a base

    that will eventually be reduced with LLL.

    The bounds are computed for the coefficients that will be used to

    form the base.

    We try to introduce only one new monomial at a time, whithout trying to

    obtain a triangular matrix.

    There are many possibilities to introduce the monomials: our goal is also 

    to introduce each of them on the diagonal with the smallest coefficient.

    The method depends on the structure of the polynomial. Here it is adapted

    to the a_n*i^n + ... + a_1 * i - t + b form.     

    Parameters

    ----------

    p: the (bivariate) polynomial;

    alpha:

    N:

    iBound:

    tBound:

    columsWidth: if == 0, no information is displayed, otherwise data is 

                 printed in colums of columnsWitdth width.

    """

    pRing = p.parent()

    polynomialsList = []

    pVariables = p.variables()

    iVariable = pVariables[0]

    tVariable = pVariables[1]

    polynomialAtPower = copy(p)

    currentPolynomial = pRing(1)     # Constant term.

    pIdegree = p.degree(iVariable)

    pTdegree = p.degree(tVariable)

    maxIdegree = pIdegree * alpha

    currentIdegree = currentPolynomial.degree(iVariable)

    nAtAlpha = N^alpha

    nAtPower = nAtAlpha

    polExpStr = ""

    #

    """

    ## Bound for iPower + pIdegree*tPower <= alpha*pIdegree

    print "degree in i:", pIdegree

    powersRangeUpperBound = alpha * pIdegree + 1 # +1 for the range.

    for iPower in xrange(0, powersRangeUpperBound):

        tPower = 0

        while (iPower + tPower * pIdegree) < powersRangeUpperBound:

                print "iPower:", iPower, " tPower:", tPower

                q = pRing(iVariable * iBound)^iPower * ((p * N)^tPower)

                print "q monomials:", q.monomials()

                polynomialsList.append(q)

                tPower += 1

    """

    """

    Start from iExp = 0 since starting from 1 does not allow for 

    resultants != 0.

    """

    for iExp in xrange(0, alpha+1):

        tExp = 0

        while iExp + tExp <= alpha:

            q = pRing(iVariable * iBound)^iExp * ((p * N)^tExp)

            sys.stdout.write("q " + str(iExp) + "," + str(tExp) + ": ")

            print q

            polynomialsList.append(q)

            tExp += 1

    return polynomialsList

    """

    # We first introduce all the monomials in i alone multiplied by N^alpha.

    for iPower in xrange(0, maxIdegree + 1):

        if columnsWidth !=0:

            polExpStr = spo_expression_as_string(iPower, iBound,

                                                 0, tBound, 

                                                 0, alpha)

            print "->", polExpStr

        currentExpression = iVariable^iPower * nAtAlpha * iBound^iPower

        currentPolynomial = pRing(currentExpression)

        polynomialsList.append(currentPolynomial)

    # End for iPower

    # We work from p^1 * N^alpha-1 to p^alpha * N^0

    for pPower in xrange(1, alpha + 1):

        # First of all the p^pPower * N^(alpha-pPower) polynomial.

        nAtPower /= N

        if columnsWidth !=0:

            polExpStr = spo_expression_as_string(0, iBound,

                                                 0, tBound,

                                                 pPower, alpha-pPower)

            print "->", polExpStr

        currentPolynomial = polynomialAtPower * nAtPower

        polynomialsList.append(currentPolynomial)

        # Exit when pPower == alpha

        if pPower == alpha:

            return polynomialsList

        for iPower in xrange(1, pIdegree + 1):

            iCurrentPower = pIdegree + iPower

            for tPower in xrange(pPower-1, 0, -1):

                #print "tPower:", tPower

                if columnsWidth != 0:

                    polExpStr = spo_expression_as_string(iCurrentPower, iBound,

                                                         tPower, tBound, 

                                                         0, alpha)

                    print "->", polExpStr

                currentExpression = i^iCurrentPower * iBound^iCurrentPower * t^tPower * tBound^tPower *nAtAlpha 

                currentPolynomial = pRing(currentExpression)

                polynomialsList.append(currentPolynomial)

                iCurrentPower += pIdegree 

            # End for tPower 

        # We now introduce the mixed i^k * t^l monomials by i^m * p^n * N^(alpha-n)

            if columnsWidth != 0:

                polExpStr = spo_expression_as_string(iPower, iBound,

                                                     0, tBound, 

                                                     pPower, alpha-pPower)

                print "->", polExpStr

            currentExpression = i^iPower * iBound^iPower * nAtPower

            currentPolynomial = pRing(currentExpression) * polynomialAtPower

            polynomialsList.append(currentPolynomial) 

        # End for iPower

        polynomialAtPower *= p  

    # End for pPower loop

    """

    return polynomialsList 

# End spo_polynomial_to_proto_matrix_6

def spo_polynomial_to_polynomials_list_7(p, alpha, N, iBound, tBound,

                                         columnsWidth=0):

    """

    As per Random Bits... direct loops nesting.

    """

    pRing = p.parent()

    polynomialsList = []

    pVariables = p.variables()

    iVariable = pVariables[0]

    tVariable = pVariables[1]

    polynomialAtPower = copy(p)

    currentPolynomial = pRing(1)     # Constant term.

    for iExp in xrange(0, alpha+1):

        pExp = 0

        while (iExp + pExp) <= alpha:

            print "iExp:", iExp, \

                "- pExp:", pExp, \

                "- alpha-pExp:", alpha-pExp

            q = pRing(iVariable * iBound)^iExp * p^pExp * N^(alpha-pExp)

            print q.monomials()

            polynomialsList.append(q)

            pExp += 1

    return polynomialsList

# End spo_polynomial_to_polynomials_list_7

def spo_polynomial_to_polynomials_list_8(p, alpha, N, iBound, tBound,

                                         columnsWidth=0):

    """

    As per Random Bits... (reversed loop nesting)

    """

    pRing = p.parent()

    polynomialsList = []

    pVariables = p.variables()

    iVariable = pVariables[0]

    tVariable = pVariables[1]

    polynomialAtPower = copy(p)

    currentPolynomial = pRing(1)     # Constant term.

    for pExp in xrange(0, alpha+1):

        iExp = 0

        while (iExp + pExp) <= alpha:

            #print "iExp:", iExp, \

            #    "- pExp:", pExp, \

            #    "- alpha-pExp:", alpha-pExp

            q = pRing(iVariable * iBound)^iExp * p^pExp * N^(alpha-pExp)

            #print q.monomials()

            polynomialsList.append(q)

            iExp += 1

    return polynomialsList

# End spo_polynomial_to_polynomials_list_8

def spo_proto_to_column_matrix(protoMatrixColumns):

    """

    Create a column (each row holds the coefficients for one monomial) matrix.

    Parameters

    ----------

    protoMatrixColumns: a list of coefficient lists.

    """

    numColumns = len(protoMatrixColumns)

    if numColumns == 0:

        return None

    # The last column holds has the maximum length. 

    numRows = len(protoMatrixColumns[numColumns-1])

    if numColumns == 0:

        return None

    baseMatrix = matrix(ZZ, numRows, numColumns)

    for colIndex in xrange(0, numColumns):

        for rowIndex in xrange(0, len(protoMatrixColumns[colIndex])):

            if protoMatrixColumns[colIndex][rowIndex] != 0:

                baseMatrix[rowIndex, colIndex] = \

                    protoMatrixColumns[colIndex][rowIndex]

    return baseMatrix

# End spo_proto_to_column_matrix.

#

def spo_proto_to_row_matrix(protoMatrixRows):

    """

    Create a row (each column holds the coefficients corresponding to one 

    monomial) matrix from the protoMatrixRows list.

    Parameters

    ----------

    protoMatrixRows: a list of coefficient lists.

    """

    numRows = len(protoMatrixRows)

    if numRows == 0:

        return None

    # Search for the longest row to get the number of columns.

    numColumns = 0

    for row in protoMatrixRows:

        rowLength = len(row)

        if numColumns < rowLength:

            numColumns = rowLength

    if numColumns == 0:

        return None

    baseMatrix = matrix(ZZ, numRows, numColumns)

    for rowIndex in xrange(0, numRows):

        for colIndex in xrange(0, len(protoMatrixRows[rowIndex])):

            if protoMatrixRows[rowIndex][colIndex] !=  0:

                baseMatrix[rowIndex, colIndex] = \

                    protoMatrixRows[rowIndex][colIndex]

            #print rowIndex, colIndex,

            #print protoMatrixRows[rowIndex][colIndex],

            #print knownMonomialsList[colIndex](boundVar1,boundVar2)

    return baseMatrix

# End spo_proto_to_row_matrix.

#

print "\t...sagePolynomialOperations loaded"