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

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

def spo_polynomial_to_matrix(p, pRing, alpha, N, columnsWidth=0):

    """

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

    to be reduced by fpLLL.

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

    p: the (bivariate) polynomial

    alpha:

    N:

    """

    knownMonomials = []

    protoMatrixRows = []

    pVariables = p.variables()

    iVariable = pVariables[0]

    tVariable = pVariables[1]

    polynomialAtPower = P(1)

    currentPolynomial = P(1)

    pIdegree = p.degree(pVariables[0])

    pTdegree = p.degree(pVariables[1])

    currentIdegree = currentPolynomial.degree(i)

    nAtPower = N^alpha

    # 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, 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 before we start "i shifts" in

        # the next iteration.

        if pPower == 0:

            # iter1: power of i

            # Notice how i^pIdegree is excluded.

            for iPower in xrange(0, pIdegree):

                # iter5: power of t.

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

                    if columnsWidth != 0:

                        print "->", spo_expression_as_string(iPower,

                                                             tPower,

                                                             pPower,

                                                             alpha)

                    currentExpression = iVariable^iPower * \

                                        tVariable^tPower * nAtPower

                    # polynomialAtPower == 1 here. Next line should be commented out but it does not work!

                    # Some convertion problem?

                    currentPolynomial = pRing(currentExpression) * \

                                        polynomialAtPower

                    pMonomials = currentPolynomial.monomials()

                    pCoefficients = currentPolynomial.coefficients()

                    spo_add_polynomial_coeffs_to_matrix(pMonomials,

                                                        pCoefficients,

                                                        knownMonomials,

                                                        protoMatrixRows,

                                                        columnsWidth)

                # End iter5.

            # End for iter1.

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

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

            # 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:

                print "->", spo_expression_as_string(0, 0, pPower, alpha)

            currentPolynomial = polynomialAtPower * nAtPower

            pMonomials = currentPolynomial.monomials()

            pCoefficients = currentPolynomial.coefficients()

            spo_add_polynomial_coeffs_to_matrix(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):

                print "->", spo_expression_as_string(iPower, 0, pPower, alpha)

                currentExpression = i^iPower * nAtPower

                currentPolynomial = P(currentExpression) * polynomialAtPower

                pMonomials = currentPolynomial.monomials()

                pCoefficients = currentPolynomial.coefficients()

                spo_add_polynomial_coeffs_to_matrix(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 mixed monomials.

            for iPower in xrange(1, pPower):

                tPower = pPower - iPower + 1

                if columnsWidth != 0:

                    print "->", spo_expression_as_string(iPower * pIdegree,

                                                         tPower,

                                                         0,

                                                         alpha)

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

                currentPolynomial = P(currentExpression)

                pMonomials = currentPolynomial.monomials()

                pCoefficients = currentPolynomial.coefficients()

                spo_add_polynomial_coeffs_to_matrix(pMonomials,

                                                    pCoefficients,

                                                    knownMonomials,

                                                    protoMatrixRows,

                                                    columnsWidth)

            # End for iPower

            # This is the main loop. It introduces:

            # - the missing mixed monomials needed before the

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

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

            # - the missing mixed monomials needed before each of the

            #   i^k * t^l * p^pPower * N^(alpha-pPower) polynomials;

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

            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:

                        print "->", spo_expression_as_string(iPower,

                                                             outerTpower,

                                                             0,

                                                             alpha)

                    currentExpression = i^iPower * t^outerTpower * nAtPower

                    currentPolynomial = P(currentExpression)

                    pMonomials = currentPolynomial.monomials()

                    pCoefficients = currentPolynomial.coefficients()

                    spo_add_polynomial_coeffs_to_matrix(pMonomials,

                                                        pCoefficients,

                                                        knownMonomials,

                                                        protoMatrixRows,

                                                        columnsWidth)

                    #print "+++++"

                    # At iShift == 0, the following 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.

                    if iShift > 0:

                        iPower = pIdegree + iShift

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

                            if columnsWidth != 0:

                                print "->", spo_expression_as_string(iPower,

                                                                     innerTpower,

                                                                     0,

                                                                     alpha)

                            currentExpression = \

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

                            currentPolynomial = P(currentExpression)

                            pMonomials = currentPolynomial.monomials()

                            pCoefficients = currentPolynomial.coefficients()

                            spo_add_polynomial_coeffs_to_matrix(pMonomials,

                                                                pCoefficients,

                                                                knownMonomials,

                                                                protoMatrixRows,

                                                                columnsWidth)

                            iPower += pIdegree

                        # End for innerTpower

                    # End of if iShift > 0

                    # Eventually, the following section introduces the

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

                    # polynomials.

                    if columnsWidth != 0:

                        print "->", spo_expression_as_string(iShift,

                                                             outerTpower,

                                                             pPower,

                                                             alpha)

                    currentExpression = i^iShift * t^outerTpower * nAtPower

                    currentPolynomial = P(currentExpression) * polynomialAtPower

                    pMonomials = currentPolynomial.monomials()

                    pCoefficients = currentPolynomial.coefficients()

                    spo_add_polynomial_coeffs_to_matrix(pMonomials,

                                                        pCoefficients,

                                                        knownMonomials,

                                                        protoMatrixRows,

                                                        columnsWidth)

                # End for outerTpower

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

            # End for iShift

        polynomialAtPower *= p

        nAtPower /= N

    # End for pPower loop

    return protoMatrixRows

# End spo_polynomial_to_matrix

def spo_add_polynomial_coeffs_to_matrix(pMonomials,

                                        pCoefficients,

                                        knownMonomials,

                                        protoMatrixRows,

                                        columnsWidth=0):

    """

    For a given polynomial (under the form of monomials and coefficents lists),

    add the coefficients of the protoMatrix (a list of proto 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.

    pMonomials     : the list of the monomials coming form some polynomial;

    pCoefficients  : the list of the corresponding coefficients to add to

                     the protoMatrix in the exact same order as the monomials;

    knownMonomials : the list of the already knonw monomials;

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

                     monomials

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

                     do display anything.

    """

    # 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

def spo_expression_as_string(powI, powT, powP, alpha):

    """

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

    output.

    """

    expressionAsString =""

    if powI != 0:

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

    if powT != 0:

        if len(expressionAsString) != 0:

            expressionAsString += " * "

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

    if powP != 0:

        if len(expressionAsString) != 0:

            expressionAsString += " * "

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

    if (alpha - powP) != 0 :

        if len(expressionAsString) != 0:

            expressionAsString += " * "

        expressionAsString += "N^" + str(alpha - powP)

    return(expressionAsString)

# End spo_expression_as_string.