Statistiques
| Révision :

root / pobysoPythonSage / src / sageSLZ / sagePolynomialOperations.sage @ 130

Historique | Voir | Annoter | Télécharger (42,79 ko)

1
load "/home/storres/recherche/arithmetique/pobysoPythonSage/src/sageSLZ/sageMatrixOperations.sage"
2
print "sagePolynomialOperations loading..."
3
def spo_add_polynomial_coeffs_to_matrix_row(poly, 
4
                                            knownMonomials, 
5
                                            protoMatrixRows, 
6
                                            columnsWidth=0):
7
    """
8
    For a given polynomial ,
9
    add the coefficients of the protoMatrix (a list of proto matrix rows).
10
    Coefficients are added to the protoMatrix row in the order imposed by the 
11
    monomials discovery list (the knownMonomials list) built as construction
12
    goes on. 
13
    As a bonus, data can be printed out for a visual check.
14
    poly           : the polynomial; in argument;
15
    knownMonomials : the list of the already known monomials; will determine
16
                     the order of the coefficients appending to a row; in-out
17
                     argument (new monomials may be discovered and then 
18
                     appended the the knowMonomials list);
19
    protoMatrixRows: a list of lists, each one holding the coefficients of the 
20
                     monomials of a polynomial; in-out argument: a new row is
21
                     added at each call;
22
    columnWith     : the width, in characters, of the displayed column ; if 0, 
23
                     do not display anything; in argument.
24
    """
25
    pMonomials   = poly.monomials()
26
    pCoefficients = poly.coefficients()
27
    # We have started with the smaller degrees in the first variable.
28
    pMonomials.reverse()
29
    pCoefficients.reverse()
30
    # New empty proto matrix row.
31
    protoMatrixRowCoefficients = []
32
    # We work according to the order of the already known monomials
33
    # No known monomials yet: add the pMonomials to knownMonomials 
34
    # and add the coefficients to the proto matrix row.
35
    if len(knownMonomials) == 0: 
36
        for pmIdx in xrange(0, len(pMonomials)):
37
            knownMonomials.append(pMonomials[pmIdx])
38
            protoMatrixRowCoefficients.append(pCoefficients[pmIdx])
39
            if columnsWidth != 0:
40
                monomialAsString = str(pCoefficients[pmIdx]) + " " + \
41
                                   str(pMonomials[pmIdx])
42
                print monomialAsString, " " * \
43
                      (columnsWidth - len(monomialAsString)),
44
    # There are some known monomials. We search for them in pMonomials and 
45
    # add their coefficients to the proto matrix row.
46
    else: 
47
        for knownMonomialIndex in xrange(0,len(knownMonomials)):
48
            # We lazily use an exception here since pMonomials.index() function
49
            # may fail throwing the ValueError exception.
50
            try:
51
                indexInPmonomials = \
52
                    pMonomials.index(knownMonomials[knownMonomialIndex])
53
                if columnsWidth != 0:
54
                    monomialAsString = str(pCoefficients[indexInPmonomials]) + \
55
                        " " + str(knownMonomials[knownMonomialIndex])
56
                    print monomialAsString, " " * \
57
                        (columnsWidth - len(monomialAsString)),
58
                # Add the coefficient to the proto matrix row and delete the \
59
                # known monomial from the current pMonomial list 
60
                #(and the corresponding coefficient as well).
61
                protoMatrixRowCoefficients.append(pCoefficients[indexInPmonomials])
62
                del pMonomials[indexInPmonomials]
63
                del pCoefficients[indexInPmonomials]
64
            # The knownMonomials element is not in pMonomials
65
            except ValueError:   
66
                protoMatrixRowCoefficients.append(0)
67
                if columnsWidth != 0:
68
                    monomialAsString = "0" + " "+ \
69
                        str(knownMonomials[knownMonomialIndex])
70
                    print monomialAsString, " " * \
71
                        (columnsWidth - len(monomialAsString)),
72
        # End for knownMonomialKey loop. 
73
        # We now append the remaining monomials of pMonomials to knownMonomials 
74
        # and the corresponding coefficients to proto matrix row.
75
        for pmIdx in xrange(0, len(pMonomials)):
76
            knownMonomials.append(pMonomials[pmIdx])
77
            protoMatrixRowCoefficients.append(pCoefficients[pmIdx])
78
            if columnsWidth != 0:
79
                monomialAsString = str(pCoefficients[pmIdx]) + " " \
80
                    + str(pMonomials[pmIdx])
81
                print monomialAsString, " " * \
82
                    (columnsWidth - len(monomialAsString)),
83
        # End for pmIdx loop.
84
    # Add the new list row elements to the proto matrix.
85
    protoMatrixRows.append(protoMatrixRowCoefficients)
86
    if columnsWidth != 0:
87
        print    
88
# End spo_add_polynomial_coeffs_to_matrix_row
89

    
90
def spo_get_coefficient_for_monomial(monomialsList, coefficientsList, monomial):
91
    """
92
    Get, for a polynomial, the coefficient for a given monomial.
93
    The polynomial is given as two lists (monomials and coefficients as
94
    return by the respective methods ; indexes of the two lists must match).
95
    If the monomial is not found, 0 is returned.
96
    """
97
    monomialIndex = 0
98
    for mono in monomialsList:
99
        if mono == monomial:
100
            return coefficientsList[monomialIndex]
101
        monomialIndex += 1
102
    return 0
103
# End spo_get_coefficient_for_monomial.
104
    
105

    
106
def spo_expression_as_string(powI, boundI, powT, boundT, powP, powN):
107
    """
108
    Computes a string version of the i^k + t^l + p^m + N^n expression for
109
    output.
110
    """
111
    expressionAsString =""
112
    if powI != 0:
113
        expressionAsString += str(iBound^powI) + " i^" + str(powI)
114
    if powT != 0:
115
        if len(expressionAsString) != 0:
116
            expressionAsString += " * "
117
        expressionAsString += str(tBound^powT) + " t^" + str(powT)
118
    if powP != 0:
119
        if len(expressionAsString) != 0:
120
            expressionAsString += " * "
121
        expressionAsString += "p^" + str(powP)
122
    if (powN) != 0 :
123
        if len(expressionAsString) != 0:
124
            expressionAsString += " * "
125
        expressionAsString += "N^" + str(powN)
126
    return(expressionAsString)
127
# End spo_expression_as_string.
128

    
129
def spo_norm(poly, p=2):
130
    """
131
    Behaves more or less (no infinity defined) as the norm for the
132
    univariate polynomials.
133
    Quoting Sage documentation:
134
    "Definition: For integer p, the p-norm of a polynomial is the pth root of 
135
    the sum of the pth powers of the absolute values of the coefficients of 
136
    the polynomial."
137
    
138
    """
139
    # TODO: check the arguments (for p see below)..
140
    norm = 0
141
    # For infinity norm.
142
    if p == Infinity:
143
        for coefficient in poly.coefficients():
144
            coefficientAbs = coefficient.abs()
145
            if coefficientAbs > norm:
146
                norm = coefficientAbs
147
        return norm
148
    # TODO: check here the value of p
149
    # p must be a positive integer >= 1.
150
    if p < 1 or (not p in ZZ):
151
        return None
152
    # For 1 norm.
153
    if p == 1:
154
        for coefficient in poly.coefficients():
155
            norm +=  coefficient.abs()
156
        return norm
157
    # For other norms
158
    for coefficient in poly.coefficients():
159
        norm +=  coefficient.abs()^p
160
    return pow(norm, 1/p)
161
# end spo_norm
162

    
163
def spo_polynomial_to_proto_matrix(p, alpha, N, columnsWidth=0):
164
    """
165
    From a (bivariate) polynomial and some other parameters build a proto 
166
    matrix (an array of "rows") to be converted into a "true" matrix and 
167
    eventually by reduced by fpLLL.
168
    The matrix is such as those found in Boneh-Durphee and Stehlé.
169
    
170
    Parameters
171
    ----------
172
    p: the (bivariate) polynomial;
173
    pRing: the ring over which p is defined;
174
    alpha:
175
    N:
176
    columsWidth: if == 0, no information is displayed, otherwise data is 
177
                 printed in colums of columnsWitdth width.
178
    """
179
    pRing = p.parent()
180
    knownMonomials = []
181
    protoMatrixRows = []
182
    polynomialsList = []
183
    pVariables = p.variables()
184
    #print "In spo...", p, p.variables()
185
    iVariable = pVariables[0]
186
    tVariable = pVariables[1]
187
    polynomialAtPower = pRing(1)
188
    currentPolynomial = pRing(1)
189
    pIdegree = p.degree(pVariables[0])
190
    pTdegree = p.degree(pVariables[1])
191
    currentIdegree = currentPolynomial.degree(iVariable)
192
    nAtAlpha = N^alpha
193
    nAtPower = nAtAlpha
194
    polExpStr = ""
195
    # We work from p^0 * N^alpha to p^alpha * N^0
196
    for pPower in xrange(0, alpha + 1):
197
        # pPower == 0 is a special case. We introduce all the monomials but one
198
        # in i and those in t necessary to be able to introduce
199
        # p. We arbitrary choose to introduce the highest degree monomial in i
200
        # with p. We also introduce all the mixed i^k * t^l monomials with
201
        # k < p.degree(i) and l <= p.degree(t).
202
        # Mixed terms introduction is necessary here before we start "i shifts"
203
        # in the next iteration.
204
        if pPower == 0:
205
            # Notice that i^pIdegree is excluded as the bound of the xrange is
206
            # pIdegree
207
            for iPower in xrange(0, pIdegree): 
208
                for tPower in xrange(0, pTdegree + 1):
209
                    if columnsWidth != 0:
210
                        polExpStr = spo_expression_as_string(iPower, 
211
                                                             tPower, 
212
                                                             pPower, 
213
                                                             alpha-pPower)
214
                        print "->", polExpStr
215
                    currentExpression = iVariable^iPower * \
216
                                        tVariable^tPower * nAtAlpha
217
                    # polynomialAtPower == 1 here. Next line should be commented
218
                    # out but it does not work! Some conversion problem?
219
                    currentPolynomial = pRing(currentExpression)
220
                    polynomialsList.append(currentPolynomial) 
221
                    pMonomials = currentPolynomial.monomials()
222
                    pCoefficients = currentPolynomial.coefficients()
223
                    spo_add_polynomial_coeffs_to_matrix_row(pMonomials, 
224
                                                            pCoefficients, 
225
                                                            knownMonomials, 
226
                                                            protoMatrixRows, 
227
                                                            columnsWidth)
228
                # End tPower.
229
            # End for iPower.
230
        else: # pPower > 0: (p^1..p^alpha)
231
            # This where we introduce the p^pPower * N^(alpha-pPower)
232
            # polynomial.
233
            # This step could technically be fused as the first iteration
234
            # of the next loop (with iPower starting at 0).
235
            # We set it apart for clarity.
236
            if columnsWidth != 0:
237
                polExpStr = spo_expression_as_string(0, 0, pPower, alpha-pPower)
238
                print "->", polExpStr
239
            currentPolynomial = polynomialAtPower * nAtPower
240
            polynomialsList.append(currentPolynomial) 
241
            pMonomials = currentPolynomial.monomials()
242
            pCoefficients = currentPolynomial.coefficients()
243
            spo_add_polynomial_coeffs_to_matrix_row(pMonomials, 
244
                                                    pCoefficients, 
245
                                                    knownMonomials, 
246
                                                    protoMatrixRows, 
247
                                                    columnsWidth)
248
            
249
            # The i^iPower * p^pPower polynomials: they add i^k monomials to  
250
            # p^pPower up to k < pIdegree * pPower. This only introduces i^k 
251
            # monomials since mixed terms (that were introduced at a previous
252
            # stage) are only shifted to already existing 
253
            # ones. p^pPower is "shifted" to higher degrees in i as far as 
254
            # possible, one step short of the degree in i of p^(pPower+1) .
255
            # These "pure" i^k monomials can only show up with i multiplications.
256
            for iPower in xrange(1, pIdegree):
257
                if columnsWidth != 0:
258
                    polExpStr = spo_expression_as_string(iPower, \
259
                                                         0,      \
260
                                                         pPower, \
261
                                                         alpha)
262
                    print "->", polExpStr
263
                currentExpression = i^iPower * nAtPower
264
                currentPolynomial = pRing(currentExpression) * polynomialAtPower
265
                polynomialsList.append(currentPolynomial) 
266
                pMonomials = currentPolynomial.monomials()
267
                pCoefficients = currentPolynomial.coefficients()
268
                spo_add_polynomial_coeffs_to_matrix_row(pMonomials,      \
269
                                                        pCoefficients,   \
270
                                                        knownMonomials,  \
271
                                                        protoMatrixRows, \
272
                                                        columnsWidth)
273
            # End for iPower
274
            # We want now to introduce a t * p^pPower polynomial. But before
275
            # that we must introduce some mixed monomials.
276
            # This loop is no triggered before pPower == 2.
277
            # It introduces a first set of high i degree mixed monomials.
278
            for iPower in xrange(1, pPower):
279
                tPower = pPower - iPower + 1
280
                if columnsWidth != 0:
281
                    polExpStr = spo_expression_as_string(iPower * pIdegree, 
282
                                                         tPower, 
283
                                                         0,  
284
                                                         alpha)
285
                    print "->", polExpStr
286
                currentExpression = i^(iPower * pIdegree) * t^tPower * nAtAlpha
287
                currentPolynomial = pRing(currentExpression)
288
                polynomialsList.append(currentPolynomial) 
289
                pMonomials = currentPolynomial.monomials()
290
                pCoefficients = currentPolynomial.coefficients()
291
                spo_add_polynomial_coeffs_to_matrix_row(pMonomials, 
292
                                                        pCoefficients, 
293
                                                        knownMonomials, 
294
                                                        protoMatrixRows, 
295
                                                        columnsWidth)
296
            # End for iPower
297
            #
298
            # This is the mixed monomials main loop. It introduces:
299
            # - the missing mixed monomials needed before the 
300
            #   t^l * p^pPower * N^(alpha-pPower) polynomial;
301
            # - the t^l * p^pPower * N^(alpha-pPower) itself;
302
            # - for each of i^k * t^l * p^pPower * N^(alpha-pPower) polynomials:
303
            #   - the the missing mixed monomials needed  polynomials,
304
            #   - the i^k * t^l * p^pPower * N^(alpha-pPower) itself.
305
            # The t^l * p^pPower * N^(alpha-pPower) is introduced when
306
            # 
307
            for iShift in xrange(0, pIdegree):
308
                # When pTdegree == 1, the following loop only introduces 
309
                # a single new monomial.
310
                #print "++++++++++"
311
                for outerTpower in xrange(1, pTdegree + 1):
312
                    # First one high i degree mixed monomial.
313
                    iPower = iShift + pPower * pIdegree
314
                    if columnsWidth != 0:
315
                        polExpStr = spo_expression_as_string(iPower, 
316
                                                             outerTpower,
317
                                                             0,
318
                                                             alpha)
319
                        print "->", polExpStr
320
                    currentExpression = i^iPower * t^outerTpower * nAtAlpha
321
                    currentPolynomial = pRing(currentExpression)
322
                    polynomialsList.append(currentPolynomial) 
323
                    pMonomials = currentPolynomial.monomials()
324
                    pCoefficients = currentPolynomial.coefficients()
325
                    spo_add_polynomial_coeffs_to_matrix_row(pMonomials, 
326
                                                            pCoefficients, 
327
                                                            knownMonomials, 
328
                                                            protoMatrixRows, 
329
                                                            columnsWidth)
330
                    #print "+++++"
331
                    # At iShift == 0, the following innerTpower loop adds  
332
                    # duplicate monomials, since no extra i^l * t^k is needed 
333
                    # before introducing the  
334
                    # i^iShift * t^outerPpower * p^pPower * N^(alpha-pPower) 
335
                    # polynomial.
336
                    # It introduces smaller i degree monomials than the
337
                    # one(s) added previously (no pPower multiplication).
338
                    # Here the exponent of t decreases as that of i increases.
339
                    # This conditional is not entered before pPower == 1.
340
                    # The innerTpower loop does not produce anything before
341
                    # pPower == 2. We keep it anyway for other configuration of
342
                    # p.
343
                    if iShift > 0:
344
                        iPower = pIdegree + iShift
345
                        for innerTpower in xrange(pPower, 1, -1):
346
                            if columnsWidth != 0:
347
                                polExpStr = spo_expression_as_string(iPower, 
348
                                                                     innerTpower,
349
                                                                     0,
350
                                                                     alpha) 
351
                            currentExpression = \
352
                                    i^(iPower) * t^(innerTpower) * nAtAlpha
353
                            currentPolynomial = pRing(currentExpression)
354
                            polynomialsList.append(currentPolynomial) 
355
                            pMonomials = currentPolynomial.monomials()
356
                            pCoefficients = currentPolynomial.coefficients()
357
                            spo_add_polynomial_coeffs_to_matrix_row(pMonomials, 
358
                                                                pCoefficients, 
359
                                                                knownMonomials, 
360
                                                                protoMatrixRows, 
361
                                                                columnsWidth)
362
                            iPower += pIdegree
363
                        # End for innerTpower
364
                    # End of if iShift > 0
365
                    # When iShift == 0, just after each of the  
366
                    # p^pPower * N^(alpha-pPower) polynomials has 
367
                    # been introduced (followed by a string of 
368
                    # i^k * p^pPower * N^(alpha-pPower) polynomials) a
369
                    # t^l *  p^pPower * N^(alpha-pPower) is introduced here.
370
                    # 
371
                    # Eventually, the following section introduces the 
372
                    # i^iShift * t^outerTpower * p^iPower * N^(alpha-pPower) 
373
                    # polynomials.
374
                    if columnsWidth != 0:
375
                        polExpStr = spo_expression_as_string(iShift, 
376
                                                             outerTpower,
377
                                                             pPower,
378
                                                             alpha-pPower)
379
                        print "->", polExpStr
380
                    currentExpression = i^iShift * t^outerTpower * nAtPower
381
                    currentPolynomial = pRing(currentExpression) * \
382
                                            polynomialAtPower
383
                    polynomialsList.append(currentPolynomial) 
384
                    pMonomials = currentPolynomial.monomials()
385
                    pCoefficients = currentPolynomial.coefficients()
386
                    spo_add_polynomial_coeffs_to_matrix_row(pMonomials, 
387
                                                            pCoefficients, 
388
                                                            knownMonomials, 
389
                                                            protoMatrixRows, 
390
                                                            columnsWidth)
391
                # End for outerTpower 
392
                #print "++++++++++"
393
            # End for iShift
394
        polynomialAtPower *= p  
395
        nAtPower /= N
396
    # End for pPower loop
397
    return ((protoMatrixRows, knownMonomials, polynomialsList))
398
# End spo_polynomial_to_proto_matrix
399

    
400
def spo_polynomial_to_polynomials_list_2(p, alpha, N, iBound, tBound,
401
                                         columnsWidth=0):
402
    """
403
    Badly out of sync code: check with versions 3 or 4.
404
    
405
    From p, alpha, N build a list of polynomials...
406
    TODO: clean up the comments below!
407
    
408
    From a (bivariate) polynomial and some other parameters build a proto 
409
    matrix (an array of "rows") to be converted into a "true" matrix and 
410
    eventually by reduced by fpLLL.
411
    The matrix is based on a list of polynomials that are built in a way
412
    that one and only monomial is added at each new polynomial. Among the many
413
    possible ways to build this list we pick one strongly dependent on the 
414
    structure of the polynomial and of the problem.
415
    We consider here the polynomials of the form: 
416
    a_k*i^k + a_(k-1)*i^(k-1) + ... + a_1*i + a_0 - t 
417
    The values of i and t are bounded and we eventually look for (i_0,t_0) 
418
    pairs such that:
419
    a_k*i_0^k + a_(k-1)*i_0^(k-1) + ... + a_1*i_0 + a_0 = t_0
420
    Hence, departing from the procedure in described in Boneh-Durfee, we will 
421
    not use "t-shifts" but only "i-shifts".
422
        
423
    Parameters
424
    ----------
425
    p: the (bivariate) polynomial;
426
    pRing: the ring over which p is defined;
427
    alpha:
428
    N:
429
    columsWidth: if == 0, no information is displayed, otherwise data is 
430
                 printed in colums of columnsWitdth width.
431
    """
432
    pRing = p.parent()
433
    polynomialsList = []
434
    pVariables = p.variables()
435
    iVariable = pVariables[0]
436
    tVariable = pVariables[1]
437
    polynomialAtPower = pRing(1)
438
    currentPolynomial = pRing(1)
439
    pIdegree = p.degree(iVariable)
440
    pTdegree = p.degree(tVariable)
441
    currentIdegree = currentPolynomial.degree(iVariable)
442
    nAtAlpha = N^alpha
443
    nAtPower = nAtAlpha
444
    polExpStr = ""
445
    # We work from p^0 * N^alpha to p^alpha * N^0
446
    for pPower in xrange(0, alpha + 1):
447
        # pPower == 0 is a special case. We introduce all the monomials in i 
448
        # up to i^pIdegree.
449
        if pPower == 0:
450
            # Notice who iPower runs up to i^pIdegree.
451
            for iPower in xrange(0, pIdegree + 1): 
452
                # No t power is taken into account as we limit our selves to
453
                # degree 1 in t and make no "t-shifts".
454
                if columnsWidth != 0:
455
                    polExpStr = spo_expression_as_string(iPower,
456
                                                         iBound, 
457
                                                         0, 
458
                                                         tBound,
459
                                                         0, 
460
                                                         alpha)
461
                    print "->", polExpStr
462
                currentExpression = iVariable^iPower * nAtAlpha * iBound^iPower
463
                # polynomialAtPower == 1 here. Next line should be commented
464
                # out but it does not work! Some conversion problem?
465
                currentPolynomial = pRing(currentExpression)
466
                polynomialsList.append(currentPolynomial) 
467
            # End for iPower.
468
        else: # pPower > 0: (p^1..p^alpha)
469
            # This where we introduce the p^pPower * N^(alpha-pPower)
470
            # polynomial. This is also where the t^pPower monomials shows up for
471
            # the first time.
472
            if columnsWidth != 0:
473
                polExpStr = spo_expression_as_string(0, iBound, 0, tBound, \
474
                                                     pPower, alpha-pPower)
475
                print "->", polExpStr
476
            currentPolynomial = polynomialAtPower * nAtPower
477
            polynomialsList.append(currentPolynomial) 
478
            # Exit when pPower == alpha
479
            if pPower == alpha:
480
                return polynomialsList
481
            # This is where the "i-shifts" take place. Mixed terms, i^k * t^l
482
            # (that were introduced at a previous
483
            # stage or are introduced now) are only shifted to already existing 
484
            # ones with the notable exception of i^iPower * t^pPower, which
485
            # must be manually introduced.
486
            # p^pPower is "shifted" to higher degrees in i as far as 
487
            # possible, up to of the degree in i of p^(pPower+1).
488
            # These "pure" i^k monomials can only show up with i multiplications.
489
            for iPower in xrange(1, pIdegree + 1):
490
                # The i^iPower * t^pPower monomial. Notice the alpha exponent
491
                # for N.
492
                internalIpower = iPower
493
                for tPower in xrange(pPower,0,-1):
494
                    if columnsWidth != 0:
495
                        polExpStr = spo_expression_as_string(internalIpower, 
496
                                                             iBound,
497
                                                             tPower,
498
                                                             tBound,
499
                                                             0,
500
                                                             alpha)
501
                        print "->", polExpStr
502
                    currentExpression = i^internalIpower * t^tPower * \
503
                                        nAtAlpha * iBound^internalIpower * \
504
                                        tBound^tPower
505
                    
506
                    currentPolynomial = pRing(currentExpression)
507
                    polynomialsList.append(currentPolynomial) 
508
                    internalIpower += pIdegree
509
                # End for tPower
510
                # The i^iPower * p^pPower * N^(alpha-pPower) i-shift.
511
                if columnsWidth != 0:
512
                    polExpStr = spo_expression_as_string(iPower, 
513
                                                         iBound,
514
                                                         0,      
515
                                                         tBound,
516
                                                         pPower, 
517
                                                         alpha-pPower)
518
                    print "->", polExpStr
519
                currentExpression = i^iPower * nAtPower * iBound^iPower
520
                currentPolynomial = pRing(currentExpression) * polynomialAtPower
521
                polynomialsList.append(currentPolynomial) 
522
            # End for iPower
523
        polynomialAtPower *= p  
524
        nAtPower /= N
525
    # End for pPower loop
526
    return polynomialsList
527
# End spo_polynomial_to_proto_matrix_2
528

    
529
def spo_polynomial_to_polynomials_list_3(p, alpha, N, iBound, tBound,
530
                                         columnsWidth=0):
531
    """
532
    From p, alpha, N build a list of polynomials...
533
    TODO: more in depth rationale...
534
    
535
    Our goal is to introduce each monomial with the smallest coefficient.
536
     
537
    
538
        
539
    Parameters
540
    ----------
541
    p: the (bivariate) polynomial;
542
    pRing: the ring over which p is defined;
543
    alpha:
544
    N:
545
    columsWidth: if == 0, no information is displayed, otherwise data is 
546
                 printed in colums of columnsWitdth width.
547
    """
548
    pRing = p.parent()
549
    polynomialsList = []
550
    pVariables = p.variables()
551
    iVariable = pVariables[0]
552
    tVariable = pVariables[1]
553
    polynomialAtPower = pRing(1)
554
    currentPolynomial = pRing(1)
555
    pIdegree = p.degree(iVariable)
556
    pTdegree = p.degree(tVariable)
557
    currentIdegree = currentPolynomial.degree(iVariable)
558
    nAtAlpha = N^alpha
559
    nAtPower = nAtAlpha
560
    polExpStr = ""
561
    # We work from p^0 * N^alpha to p^alpha * N^0
562
    for pPower in xrange(0, alpha + 1):
563
        # pPower == 0 is a special case. We introduce all the monomials in i 
564
        # up to i^pIdegree.
565
        if pPower == 0:
566
            # Notice who iPower runs up to i^pIdegree.
567
            for iPower in xrange(0, pIdegree + 1): 
568
                # No t power is taken into account as we limit our selves to
569
                # degree 1 in t and make no "t-shifts".
570
                if columnsWidth != 0:
571
                    polExpStr = spo_expression_as_string(iPower, 
572
                                                         iBound,
573
                                                         0, 
574
                                                         tBound,
575
                                                         0, 
576
                                                         alpha)
577
                    print "->", polExpStr
578
                currentExpression = iVariable^iPower * nAtAlpha * iBound^iPower
579
                # polynomialAtPower == 1 here. Next line should be commented
580
                # out but it does not work! Some conversion problem?
581
                currentPolynomial = pRing(currentExpression)
582
                polynomialsList.append(currentPolynomial) 
583
            # End for iPower.
584
        else: # pPower > 0: (p^1..p^alpha)
585
            # This where we introduce the p^pPower * N^(alpha-pPower)
586
            # polynomial. This is also where the t^pPower monomials shows up for
587
            # the first time. It app
588
            if columnsWidth != 0:
589
                polExpStr = spo_expression_as_string(0, iBound, 
590
                                                     0, tBound,
591
                                                     pPower, alpha-pPower)
592
                print "->", polExpStr
593
            currentPolynomial = polynomialAtPower * nAtPower
594
            polynomialsList.append(currentPolynomial) 
595
            # Exit when pPower == alpha
596
            if pPower == alpha:
597
                return polynomialsList
598
            # This is where the "i-shifts" take place. Mixed terms, i^k * t^l
599
            # (that were introduced at a previous
600
            # stage or are introduced now) are only shifted to already existing 
601
            # ones with the notable exception of i^iPower * t^pPower, which
602
            # must be manually introduced.
603
            # p^pPower is "shifted" to higher degrees in i as far as 
604
            # possible, up to of the degree in i of p^(pPower+1).
605
            # These "pure" i^k monomials can only show up with i multiplications.
606
            for iPower in xrange(1, pIdegree + 1):
607
                # The i^iPower * t^pPower monomial. Notice the alpha exponent
608
                # for N.
609
                internalIpower = iPower
610
                for tPower in xrange(pPower,0,-1):
611
                    if columnsWidth != 0:
612
                        polExpStr = spo_expression_as_string(internalIpower,
613
                                                             iBound, 
614
                                                             tPower,  
615
                                                             tBound,
616
                                                             0,
617
                                                             alpha)
618
                        print "->", polExpStr
619
                    currentExpression = i^internalIpower * t^tPower * nAtAlpha * \
620
                                        iBound^internalIpower * tBound^tPower
621
                    currentPolynomial = pRing(currentExpression)
622
                    polynomialsList.append(currentPolynomial) 
623
                    internalIpower += pIdegree
624
                # End for tPower
625
                # Here we have to choose between a 
626
                # i^iPower * p^pPower * N^(alpha-pPower) i-shift and
627
                # i^iPower * i^(d_i(p) * pPower) * N^alpha, depending on which 
628
                # coefficient is smallest.
629
                IcurrentExponent = iPower + \
630
                                (pPower * polynomialAtPower.degree(iVariable))
631
                currentMonomial = pRing(iVariable^IcurrentExponent)
632
                currentPolynomial = pRing(iVariable^iPower * nAtPower * \
633
                                          iBound^iPower) * \
634
                                          polynomialAtPower
635
                currMonomials = currentPolynomial.monomials()
636
                currCoefficients = currentPolynomial.coefficients()
637
                currentCoefficient = spo_get_coefficient_for_monomial( \
638
                                                    currMonomials, 
639
                                                    currCoefficients, 
640
                                                    currentMonomial)
641
                print "Current coefficient:", currentCoefficient
642
                alterCoefficient = iBound^IcurrentExponent * nAtAlpha 
643
                print "N^alpha * ibound^", IcurrentExponent, ":", \
644
                         alterCoefficient
645
                if currentCoefficient > alterCoefficient :
646
                    if columnsWidth != 0:
647
                        polExpStr = spo_expression_as_string(IcurrentExponent,
648
                                                             iBound,
649
                                                             0,
650
                                                             tBound,
651
                                                             0,
652
                                                             alpha)
653
                        print "->", polExpStr
654
                    polynomialsList.append(currentMonomial * \
655
                                           alterCoefficient)
656
                else:
657
                    if columnsWidth != 0:
658
                        polExpStr = spo_expression_as_string(iPower, iBound,
659
                                                             0, tBound,
660
                                                             pPower, 
661
                                                             alpha-pPower)
662
                        print "->", polExpStr
663
                    polynomialsList.append(currentPolynomial) 
664
            # End for iPower
665
        polynomialAtPower *= p  
666
        nAtPower /= N
667
    # End for pPower loop
668
    return polynomialsList 
669
# End spo_polynomial_to_proto_matrix_3
670

    
671
def spo_polynomial_to_polynomials_list_4(p, alpha, N, iBound, tBound,
672
                                         columnsWidth=0):
673
    """
674
    From p, alpha, N build a list of polynomials...
675
    TODO: more in depth rationale...
676
    
677
    Our goal is to introduce each monomial with the smallest coefficient.
678
     
679
    
680
        
681
    Parameters
682
    ----------
683
    p: the (bivariate) polynomial;
684
    pRing: the ring over which p is defined;
685
    alpha:
686
    N:
687
    columsWidth: if == 0, no information is displayed, otherwise data is 
688
                 printed in colums of columnsWitdth width.
689
    """
690
    pRing = p.parent()
691
    polynomialsList = []
692
    pVariables = p.variables()
693
    iVariable = pVariables[0]
694
    tVariable = pVariables[1]
695
    polynomialAtPower = copy(p)
696
    currentPolynomial = pRing(1)
697
    pIdegree = p.degree(iVariable)
698
    pTdegree = p.degree(tVariable)
699
    maxIdegree = pIdegree * alpha
700
    currentIdegree = currentPolynomial.degree(iVariable)
701
    nAtAlpha = N^alpha
702
    nAtPower = nAtAlpha
703
    polExpStr = ""
704
    # We first introduce all the monomials in i alone multiplied by N^alpha.
705
    for iPower in xrange(0, maxIdegree + 1):
706
        if columnsWidth !=0:
707
            polExpStr = spo_expression_as_string(iPower, iBound,
708
                                                 0, tBound, 
709
                                                 0, alpha)
710
            print "->", polExpStr
711
        currentExpression = iVariable^iPower * nAtAlpha * iBound^iPower
712
        currentPolynomial = pRing(currentExpression)
713
        polynomialsList.append(currentPolynomial)
714
    # End for iPower
715
    # We work from p^1 * N^alpha-1 to p^alpha * N^0
716
    for pPower in xrange(1, alpha + 1):
717
        # First of all the p^pPower * N^(alpha-pPower) polynomial.
718
        nAtPower /= N
719
        if columnsWidth !=0:
720
            polExpStr = spo_expression_as_string(0, iBound,
721
                                                 0, tBound,
722
                                                 pPower, alpha-pPower)
723
            print "->", polExpStr
724
        currentPolynomial = polynomialAtPower * nAtPower
725
        polynomialsList.append(currentPolynomial)
726
        # Exit when pPower == alpha
727
        if pPower == alpha:
728
            return polynomialsList
729
        # We now introduce the mixed i^k * t^l monomials by i^m * p^n * N^(alpha-n)
730
        for iPower in xrange(1, pIdegree + 1): 
731
            if columnsWidth != 0:
732
                polExpStr = spo_expression_as_string(iPower, iBound,
733
                                                     0, tBound, 
734
                                                     pPower, alpha-pPower)
735
                print "->", polExpStr
736
            currentExpression = i^iPower * iBound^iPower * nAtPower
737
            currentPolynomial = pRing(currentExpression) * polynomialAtPower
738
            polynomialsList.append(currentPolynomial) 
739
        # End for iPower
740
        polynomialAtPower *= p  
741
    # End for pPower loop
742
    return polynomialsList 
743
# End spo_polynomial_to_proto_matrix_4
744

    
745
def spo_polynomial_to_polynomials_list_5(p, alpha, N, iBound, tBound,
746
                                         columnsWidth=0):
747
    """
748
    From p, alpha, N build a list of polynomials use to create a base
749
    that will eventually be reduced with LLL.
750
    
751
    The bounds are computed for the coefficients that will be used to
752
    form the base.
753
    
754
    We try to introduce only one new monomial at a time, to obtain a
755
    triangular matrix (it is easy to compute the volume of the underlining
756
    latice if the matrix is triangular).
757

    
758
    There are many possibilities to introduce the monomials: our goal is also 
759
    to introduce each of them on the diagonal with the smallest coefficient.
760

    
761
    The method depends on the structure of the polynomial. Here it is adapted
762
    to the a_n*i^n + ... + a_1 * i - t + b form.     
763
        
764
    Parameters
765
    ----------
766
    p: the (bivariate) polynomial;
767
    alpha:
768
    N:
769
    iBound:
770
    tBound:
771
    columsWidth: if == 0, no information is displayed, otherwise data is 
772
                 printed in colums of columnsWitdth width.
773
    """
774
    pRing = p.parent()
775
    polynomialsList = []
776
    pVariables = p.variables()
777
    iVariable = pVariables[0]
778
    tVariable = pVariables[1]
779
    polynomialAtPower = copy(p)
780
    currentPolynomial = pRing(1)
781
    pIdegree = p.degree(iVariable)
782
    pTdegree = p.degree(tVariable)
783
    maxIdegree = pIdegree * alpha
784
    currentIdegree = currentPolynomial.degree(iVariable)
785
    nAtAlpha = N^alpha
786
    nAtPower = nAtAlpha
787
    polExpStr = ""
788
    # We first introduce all the monomials in i alone multiplied by N^alpha.
789
    for iPower in xrange(0, maxIdegree + 1):
790
        if columnsWidth !=0:
791
            polExpStr = spo_expression_as_string(iPower, iBound,
792
                                                 0, tBound, 
793
                                                 0, alpha)
794
            print "->", polExpStr
795
        currentExpression = iVariable^iPower * nAtAlpha * iBound^iPower
796
        currentPolynomial = pRing(currentExpression)
797
        polynomialsList.append(currentPolynomial)
798
    # End for iPower
799
    # We work from p^1 * N^alpha-1 to p^alpha * N^0
800
    for pPower in xrange(1, alpha + 1):
801
        # First of all the p^pPower * N^(alpha-pPower) polynomial.
802
        nAtPower /= N
803
        if columnsWidth !=0:
804
            polExpStr = spo_expression_as_string(0, iBound,
805
                                                 0, tBound,
806
                                                 pPower, alpha-pPower)
807
            print "->", polExpStr
808
        currentPolynomial = polynomialAtPower * nAtPower
809
        polynomialsList.append(currentPolynomial)
810
        # Exit when pPower == alpha
811
        if pPower == alpha:
812
            return polynomialsList
813
        for iPower in xrange(1, pIdegree + 1):
814
            iCurrentPower = pIdegree + iPower
815
            for tPower in xrange(pPower-1, 0, -1):
816
                #print "tPower:", tPower
817
                if columnsWidth != 0:
818
                    polExpStr = spo_expression_as_string(iCurrentPower, iBound,
819
                                                         tPower, tBound, 
820
                                                         0, alpha)
821
                    print "->", polExpStr
822
                currentExpression = i^iCurrentPower * iBound^iCurrentPower * t^tPower * tBound^tPower *nAtAlpha 
823
                currentPolynomial = pRing(currentExpression)
824
                polynomialsList.append(currentPolynomial)
825
                iCurrentPower += pIdegree 
826
            # End for tPower 
827
        # We now introduce the mixed i^k * t^l monomials by i^m * p^n * N^(alpha-n)
828
            if columnsWidth != 0:
829
                polExpStr = spo_expression_as_string(iPower, iBound,
830
                                                     0, tBound, 
831
                                                     pPower, alpha-pPower)
832
                print "->", polExpStr
833
            currentExpression = i^iPower * iBound^iPower * nAtPower
834
            currentPolynomial = pRing(currentExpression) * polynomialAtPower
835
            polynomialsList.append(currentPolynomial) 
836
        # End for iPower
837
        polynomialAtPower *= p  
838
    # End for pPower loop
839
    return polynomialsList 
840
# End spo_polynomial_to_proto_matrix_5
841

    
842
def spo_proto_to_column_matrix(protoMatrixColumns):
843
    """
844
    Create a column (each row holds the coefficients for one monomial) matrix.
845
    
846
    Parameters
847
    ----------
848
    protoMatrixColumns: a list of coefficient lists.
849
    """
850
    numColumns = len(protoMatrixColumns)
851
    if numColumns == 0:
852
        return None
853
    # The last column holds has the maximum length. 
854
    numRows = len(protoMatrixColumns[numColumns-1])
855
    if numColumns == 0:
856
        return None
857
    baseMatrix = matrix(ZZ, numRows, numColumns)
858
    for colIndex in xrange(0, numColumns):
859
        for rowIndex in xrange(0, len(protoMatrixColumns[colIndex])):
860
            if protoMatrixColumns[colIndex][rowIndex] != 0:
861
                baseMatrix[rowIndex, colIndex] = \
862
                    protoMatrixColumns[colIndex][rowIndex]
863
    return baseMatrix
864
# End spo_proto_to_column_matrix.
865
#
866
def spo_proto_to_row_matrix(protoMatrixRows):
867
    """
868
    Create a row (each column holds the coefficients corresponding to one 
869
    monomial) matrix from the protoMatrixRows list.
870
    
871
    Parameters
872
    ----------
873
    protoMatrixRows: a list of coefficient lists.
874
    """
875
    numRows = len(protoMatrixRows)
876
    if numRows == 0:
877
        return None
878
    # The last row is the longest one.
879
    numColumns = len(protoMatrixRows[numRows-1])
880
    if numColumns == 0:
881
        return None
882
    baseMatrix = matrix(ZZ, numRows, numColumns)
883
    for rowIndex in xrange(0, numRows):
884
        for colIndex in xrange(0, len(protoMatrixRows[rowIndex])):
885
            if protoMatrixRows[rowIndex][colIndex] !=  0:
886
                baseMatrix[rowIndex, colIndex] = \
887
                    protoMatrixRows[rowIndex][colIndex]
888
            #print rowIndex, colIndex,
889
            #print protoMatrixRows[rowIndex][colIndex],
890
            #print knownMonomialsList[colIndex](boundVar1,boundVar2)
891
    return baseMatrix
892
# End spo_proto_to_row_matrix.
893
#
894
print "\t...sagePolynomialOperations loaded"