Statistiques
| Révision :

root / pobysoPythonSage / src / sageSLZ / sageSLZ.sage @ 279

Historique | Voir | Annoter | Télécharger (130,06 ko)

1
r"""
2
Sage core functions needed for the implementation of SLZ.
3

    
4
AUTHORS:
5
- S.T. (2013-08): initial version
6

    
7
Examples:
8

    
9
TODO::
10
"""
11
sys.stderr.write("sageSLZ loading...\n")
12
#
13
import inspect
14
#
15
def slz_compute_binade(number):
16
    """"
17
    For a given number, compute the "binade" that is integer m such that
18
    2^m <= number < 2^(m+1). If number == 0 return None.
19
    """
20
    # Checking the parameter.
21
    # The exception construction is used to detect if number is a RealNumber
22
    # since not all numbers have
23
    # the mro() method. sage.rings.real_mpfr.RealNumber do.
24
    try:
25
        classTree = [number.__class__] + number.mro()
26
        # If the number is not a RealNumber (or offspring thereof) try
27
        # to transform it.
28
        if not sage.rings.real_mpfr.RealNumber in classTree:
29
            numberAsRR = RR(number)
30
        else:
31
            numberAsRR = number
32
    except AttributeError:
33
        return None
34
    # Zero special case.
35
    if numberAsRR == 0:
36
        return RR(-infinity)
37
    else:
38
        realField           = numberAsRR.parent()
39
        numberLog2          = numberAsRR.abs().log2()
40
        floorNumberLog2     = floor(numberLog2)
41
        ## Do not get caught by rounding of log2() both ways.
42
        ## When numberLog2 is an integer, compare numberAsRR
43
        #  with 2^numberLog2.
44
        if floorNumberLog2 == numberLog2:
45
            if numberAsRR.abs() < realField(2^floorNumberLog2):
46
                return floorNumberLog2 - 1
47
            else:
48
                return floorNumberLog2
49
        else:
50
            return floorNumberLog2
51
# End slz_compute_binade
52

    
53
#
54
def slz_compute_binade_bounds(number, emin, emax=sys.maxint):
55
    """
56
    For given "real number", compute the bounds of the binade it belongs to.
57
    
58
    NOTE::
59
        When number >= 2^(emax+1), we return the "fake" binade 
60
        [2^(emax+1), +infinity]. Ditto for number <= -2^(emax+1)
61
        with interval [-infinity, -2^(emax+1)]. We want to distinguish
62
        this case from that of "really" invalid arguments.
63
        
64
    """
65
    # Check the parameters.
66
    # RealNumbers or RealNumber offspring only.
67
    # The exception construction is necessary since not all objects have
68
    # the mro() method. sage.rings.real_mpfr.RealNumber do.
69
    try:
70
        classTree = [number.__class__] + number.mro()
71
        if not sage.rings.real_mpfr.RealNumber in classTree:
72
            return None
73
    except AttributeError:
74
        return None
75
    # Non zero negative integers only for emin.
76
    if emin >= 0 or int(emin) != emin:
77
        return None
78
    # Non zero positive integers only for emax.
79
    if emax <= 0 or int(emax) != emax:
80
        return None
81
    precision = number.precision()
82
    RF  = RealField(precision)
83
    if number == 0:
84
        return (RF(0),RF(2^(emin)) - RF(2^(emin-precision)))
85
    # A more precise RealField is needed to avoid unwanted rounding effects
86
    # when computing number.log2().
87
    RRF = RealField(max(2048, 2 * precision))
88
    # number = 0 special case, the binade bounds are 
89
    # [0, 2^emin - 2^(emin-precision)]
90
    # Begin general case
91
    l2 = RRF(number).abs().log2()
92
    # Another special one: beyond largest representable -> "Fake" binade.
93
    if l2 >= emax + 1:
94
        if number > 0:
95
            return (RF(2^(emax+1)), RF(+infinity) )
96
        else:
97
            return (RF(-infinity), -RF(2^(emax+1)))
98
    # Regular case cont'd.
99
    offset = int(l2)
100
    # number.abs() >= 1.
101
    if l2 >= 0:
102
        if number >= 0:
103
            lb = RF(2^offset)
104
            ub = RF(2^(offset + 1) - 2^(-precision+offset+1))
105
        else: #number < 0
106
            lb = -RF(2^(offset + 1) - 2^(-precision+offset+1))
107
            ub = -RF(2^offset)
108
    else: # log2 < 0, number.abs() < 1.
109
        if l2 < emin: # Denormal
110
           # print "Denormal:", l2
111
            if number >= 0:
112
                lb = RF(0)
113
                ub = RF(2^(emin)) - RF(2^(emin-precision))
114
            else: # number <= 0
115
                lb = - RF(2^(emin)) + RF(2^(emin-precision))
116
                ub = RF(0)
117
        elif l2 > emin: # Normal number other than +/-2^emin.
118
            if number >= 0:
119
                if int(l2) == l2:
120
                    lb = RF(2^(offset))
121
                    ub = RF(2^(offset+1)) - RF(2^(-precision+offset+1))
122
                else:
123
                    lb = RF(2^(offset-1))
124
                    ub = RF(2^(offset)) - RF(2^(-precision+offset))
125
            else: # number < 0
126
                if int(l2) == l2: # Binade limit.
127
                    lb = -RF(2^(offset+1) - 2^(-precision+offset+1))
128
                    ub = -RF(2^(offset))
129
                else:
130
                    lb = -RF(2^(offset) - 2^(-precision+offset))
131
                    ub = -RF(2^(offset-1))
132
        else: # l2== emin, number == +/-2^emin
133
            if number >= 0:
134
                lb = RF(2^(offset))
135
                ub = RF(2^(offset+1)) - RF(2^(-precision+offset+1))
136
            else: # number < 0
137
                lb = -RF(2^(offset+1) - 2^(-precision+offset+1))
138
                ub = -RF(2^(offset))
139
    return (lb, ub)
140
# End slz_compute_binade_bounds
141
#
142
def slz_compute_coppersmith_reduced_polynomials(inputPolynomial,
143
                                                 alpha,
144
                                                 N,
145
                                                 iBound,
146
                                                 tBound,
147
                                                 debug = False):
148
    """
149
    For a given set of arguments (see below), compute a list
150
    of "reduced polynomials" that could be used to compute roots
151
    of the inputPolynomial.
152
    INPUT:
153
    
154
    - "inputPolynomial" -- (no default) a bivariate integer polynomial;
155
    - "alpha" -- the alpha parameter of the Coppersmith algorithm;
156
    - "N" -- the modulus;
157
    - "iBound" -- the bound on the first variable;
158
    - "tBound" -- the bound on the second variable.
159
    
160
    OUTPUT:
161
    
162
    A list of bivariate integer polynomial obtained using the Coppersmith
163
    algorithm. The polynomials correspond to the rows of the LLL-reduce
164
    reduced base that comply with the Coppersmith condition.
165
    """
166
    # Arguments check.
167
    if iBound == 0 or tBound == 0:
168
        return None
169
    # End arguments check.
170
    nAtAlpha = N^alpha
171
    ## Building polynomials for matrix.
172
    polyRing   = inputPolynomial.parent()
173
    # Whatever the 2 variables are actually called, we call them
174
    # 'i' and 't' in all the variable names.
175
    (iVariable, tVariable) = inputPolynomial.variables()[:2]
176
    #print polyVars[0], type(polyVars[0])
177
    initialPolynomial = inputPolynomial.subs({iVariable:iVariable * iBound,
178
                                              tVariable:tVariable * tBound})
179
    if debug:
180
        polynomialsList = \
181
            spo_polynomial_to_polynomials_list_8(initialPolynomial,
182
                                                 alpha,
183
                                                 N,
184
                                                 iBound,
185
                                                 tBound,
186
                                                 20)
187
    else:
188
        polynomialsList = \
189
            spo_polynomial_to_polynomials_list_8(initialPolynomial,
190
                                                 alpha,
191
                                                 N,
192
                                                 iBound,
193
                                                 tBound,
194
                                                 0)
195
    #print "Polynomials list:", polynomialsList
196
    ## Building the proto matrix.
197
    knownMonomials = []
198
    protoMatrix    = []
199
    if debug:
200
        for poly in polynomialsList:
201
            spo_add_polynomial_coeffs_to_matrix_row(poly,
202
                                                    knownMonomials,
203
                                                    protoMatrix,
204
                                                    20)
205
    else:
206
        for poly in polynomialsList:
207
            spo_add_polynomial_coeffs_to_matrix_row(poly,
208
                                                    knownMonomials,
209
                                                    protoMatrix,
210
                                                    0)
211
    matrixToReduce = spo_proto_to_row_matrix(protoMatrix)
212
    #print matrixToReduce
213
    ## Reduction and checking.
214
    ## S.T. changed 'fp' to None as of Sage 6.6 complying to
215
    #  error message issued when previous code was used.
216
    #reducedMatrix = matrixToReduce.LLL(fp='fp')
217
    reducedMatrix = matrixToReduce.LLL(fp=None)
218
    isLLLReduced  = reducedMatrix.is_LLL_reduced()
219
    if not isLLLReduced:
220
        return None
221
    monomialsCount     = len(knownMonomials)
222
    monomialsCountSqrt = sqrt(monomialsCount)
223
    #print "Monomials count:", monomialsCount, monomialsCountSqrt.n()
224
    #print reducedMatrix
225
    ## Check the Coppersmith condition for each row and build the reduced 
226
    #  polynomials.
227
    ccReducedPolynomialsList = []
228
    for row in reducedMatrix.rows():
229
        l2Norm = row.norm(2)
230
        if (l2Norm * monomialsCountSqrt) < nAtAlpha:
231
            #print (l2Norm * monomialsCountSqrt).n()
232
            #print l2Norm.n()
233
            ccReducedPolynomial = \
234
                slz_compute_reduced_polynomial(row,
235
                                               knownMonomials,
236
                                               iVariable,
237
                                               iBound,
238
                                               tVariable,
239
                                               tBound)
240
            if not ccReducedPolynomial is None:
241
                ccReducedPolynomialsList.append(ccReducedPolynomial)
242
        else:
243
            #print l2Norm.n() , ">", nAtAlpha
244
            pass
245
    if len(ccReducedPolynomialsList) < 2:
246
        print "Less than 2 Coppersmith condition compliant vectors."
247
        return ()
248
    #print ccReducedPolynomialsList
249
    return ccReducedPolynomialsList
250
# End slz_compute_coppersmith_reduced_polynomials
251
#
252
def slz_compute_coppersmith_reduced_polynomials_gram(inputPolynomial,
253
                                                     alpha,
254
                                                     N,
255
                                                     iBound,
256
                                                     tBound,
257
                                                     debug = False):
258
    """
259
    For a given set of arguments (see below), compute a list
260
    of "reduced polynomials" that could be used to compute roots
261
    of the inputPolynomial.
262
    INPUT:
263
    
264
    - "inputPolynomial" -- (no default) a bivariate integer polynomial;
265
    - "alpha" -- the alpha parameter of the Coppersmith algorithm;
266
    - "N" -- the modulus;
267
    - "iBound" -- the bound on the first variable;
268
    - "tBound" -- the bound on the second variable.
269
    
270
    OUTPUT:
271
    
272
    A list of bivariate integer polynomial obtained using the Coppersmith
273
    algorithm. The polynomials correspond to the rows of the LLL-reduce
274
    reduced base that comply with the Coppersmith condition.
275
    """
276
    # Arguments check.
277
    if iBound == 0 or tBound == 0:
278
        return None
279
    # End arguments check.
280
    nAtAlpha = N^alpha
281
    ## Building polynomials for matrix.
282
    polyRing   = inputPolynomial.parent()
283
    # Whatever the 2 variables are actually called, we call them
284
    # 'i' and 't' in all the variable names.
285
    (iVariable, tVariable) = inputPolynomial.variables()[:2]
286
    #print polyVars[0], type(polyVars[0])
287
    initialPolynomial = inputPolynomial.subs({iVariable:iVariable * iBound,
288
                                              tVariable:tVariable * tBound})
289
    if debug:
290
        polynomialsList = \
291
            spo_polynomial_to_polynomials_list_8(initialPolynomial,
292
                                                 alpha,
293
                                                 N,
294
                                                 iBound,
295
                                                 tBound,
296
                                                 20)
297
    else:
298
        polynomialsList = \
299
            spo_polynomial_to_polynomials_list_8(initialPolynomial,
300
                                                 alpha,
301
                                                 N,
302
                                                 iBound,
303
                                                 tBound,
304
                                                 0)
305
    #print "Polynomials list:", polynomialsList
306
    ## Building the proto matrix.
307
    knownMonomials = []
308
    protoMatrix    = []
309
    if debug:
310
        for poly in polynomialsList:
311
            spo_add_polynomial_coeffs_to_matrix_row(poly,
312
                                                    knownMonomials,
313
                                                    protoMatrix,
314
                                                    20)
315
    else:
316
        for poly in polynomialsList:
317
            spo_add_polynomial_coeffs_to_matrix_row(poly,
318
                                                    knownMonomials,
319
                                                    protoMatrix,
320
                                                    0)
321
    matrixToReduce = spo_proto_to_row_matrix(protoMatrix)
322
    #print matrixToReduce
323
    ## Reduction and checking.
324
    ### In this variant we use the Pari LLL_gram reduction function.
325
    #   It works with the Gram matrix instead of the plain matrix.
326
    matrixToReduceTransposed = matrixToReduce.transpose() 
327
    matrixToReduceGram = matrixToReduce * matrixToReduceTransposed
328
    ### LLL_gram returns a unimodular transformation matrix such that:
329
    #   umt.transpose() * matrixToReduce * umt is reduced..
330
    umt = matrixToReduceGram.LLL_gram()
331
    #print "Unimodular transformation matrix:"
332
    #for row in umt.rows():
333
    #    print row
334
    ### The computed transformation matrix is transposed and applied to the
335
    #   left side of matrixToReduce. 
336
    reducedMatrix = umt.transpose() * matrixToReduce
337
    #print "Reduced matrix:"
338
    #for row in reducedMatrix.rows():
339
    #    print row
340
    isLLLReduced  = reducedMatrix.is_LLL_reduced()
341
    #if not isLLLReduced:
342
    #    return None
343
    monomialsCount     = len(knownMonomials)
344
    monomialsCountSqrt = sqrt(monomialsCount)
345
    #print "Monomials count:", monomialsCount, monomialsCountSqrt.n()
346
    #print reducedMatrix
347
    ## Check the Coppersmith condition for each row and build the reduced 
348
    #  polynomials.
349
    ccReducedPolynomialsList = []
350
    for row in reducedMatrix.rows():
351
        l2Norm = row.norm(2)
352
        if (l2Norm * monomialsCountSqrt) < nAtAlpha:
353
            #print (l2Norm * monomialsCountSqrt).n()
354
            #print l2Norm.n()
355
            ccReducedPolynomial = \
356
                slz_compute_reduced_polynomial(row,
357
                                               knownMonomials,
358
                                               iVariable,
359
                                               iBound,
360
                                               tVariable,
361
                                               tBound)
362
            if not ccReducedPolynomial is None:
363
                ccReducedPolynomialsList.append(ccReducedPolynomial)
364
        else:
365
            #print l2Norm.n() , ">", nAtAlpha
366
            pass
367
    if len(ccReducedPolynomialsList) < 2:
368
        print "Less than 2 Coppersmith condition compliant vectors."
369
        return ()
370
    #print ccReducedPolynomialsList
371
    return ccReducedPolynomialsList
372
# End slz_compute_coppersmith_reduced_polynomials_gram
373
#
374
def slz_compute_coppersmith_reduced_polynomials_proj(inputPolynomial,
375
                                                     alpha,
376
                                                     N,
377
                                                     iBound,
378
                                                     tBound,
379
                                                     debug = False):
380
    """
381
    For a given set of arguments (see below), compute a list
382
    of "reduced polynomials" that could be used to compute roots
383
    of the inputPolynomial.
384
    INPUT:
385
    
386
    - "inputPolynomial" -- (no default) a bivariate integer polynomial;
387
    - "alpha" -- the alpha parameter of the Coppersmith algorithm;
388
    - "N" -- the modulus;
389
    - "iBound" -- the bound on the first variable;
390
    - "tBound" -- the bound on the second variable.
391
    
392
    OUTPUT:
393
    
394
    A list of bivariate integer polynomial obtained using the Coppersmith
395
    algorithm. The polynomials correspond to the rows of the LLL-reduce
396
    reduced base that comply with the Coppersmith condition.
397
    """
398
    #@par Changes from runSLZ-113.sage
399
    # LLL reduction is not performed on the matrix itself but rather on the
400
    # product of the matrix with a uniform random matrix.
401
    # The reduced matrix obtained is discarded but the transformation matrix 
402
    # obtained is used to multiply the original matrix in order to reduced it.
403
    # If a sufficient level of reduction is obtained, we stop here. If not
404
    # the product matrix obtained above is LLL reduced. But as it has been
405
    # pre-reduced at the above step, reduction is supposed to be much faster.
406
    #
407
    # Arguments check.
408
    if iBound == 0 or tBound == 0:
409
        return None
410
    # End arguments check.
411
    nAtAlpha = N^alpha
412
    ## Building polynomials for matrix.
413
    polyRing   = inputPolynomial.parent()
414
    # Whatever the 2 variables are actually called, we call them
415
    # 'i' and 't' in all the variable names.
416
    (iVariable, tVariable) = inputPolynomial.variables()[:2]
417
    #print polyVars[0], type(polyVars[0])
418
    initialPolynomial = inputPolynomial.subs({iVariable:iVariable * iBound,
419
                                              tVariable:tVariable * tBound})
420
    if debug:
421
        polynomialsList = \
422
            spo_polynomial_to_polynomials_list_8(initialPolynomial,
423
                                                 alpha,
424
                                                 N,
425
                                                 iBound,
426
                                                 tBound,
427
                                                 20)
428
    else:
429
        polynomialsList = \
430
            spo_polynomial_to_polynomials_list_8(initialPolynomial,
431
                                                 alpha,
432
                                                 N,
433
                                                 iBound,
434
                                                 tBound,
435
                                                 0)
436
    #print "Polynomials list:", polynomialsList
437
    ## Building the proto matrix.
438
    knownMonomials = []
439
    protoMatrix    = []
440
    if debug:
441
        for poly in polynomialsList:
442
            spo_add_polynomial_coeffs_to_matrix_row(poly,
443
                                                    knownMonomials,
444
                                                    protoMatrix,
445
                                                    20)
446
    else:
447
        for poly in polynomialsList:
448
            spo_add_polynomial_coeffs_to_matrix_row(poly,
449
                                                    knownMonomials,
450
                                                    protoMatrix,
451
                                                    0)
452
    matrixToReduce = spo_proto_to_row_matrix(protoMatrix)
453
    #print matrixToReduce
454
    ## Reduction and checking.
455
    ### Reduction with projection
456
    (reducedMatrixStep1, reductionMatrixStep1) = \
457
         slz_reduce_lll_proj(matrixToReduce,16)
458
    #print "Reduced matrix:"
459
    #print reducedMatrixStep1
460
    #for row in reducedMatrix.rows():
461
    #    print row
462
    monomialsCount     = len(knownMonomials)
463
    monomialsCountSqrt = sqrt(monomialsCount)
464
    #print "Monomials count:", monomialsCount, monomialsCountSqrt.n()
465
    #print reducedMatrix
466
    ## Check the Coppersmith condition for each row and build the reduced 
467
    #  polynomials.
468
    ccReducedPolynomialsList = []
469
    for row in reducedMatrixStep1.rows():
470
        l2Norm = row.norm(2)
471
        if (l2Norm * monomialsCountSqrt) < nAtAlpha:
472
            #print (l2Norm * monomialsCountSqrt).n()
473
            #print l2Norm.n()
474
            ccReducedPolynomial = \
475
                slz_compute_reduced_polynomial(row,
476
                                               knownMonomials,
477
                                               iVariable,
478
                                               iBound,
479
                                               tVariable,
480
                                               tBound)
481
            if not ccReducedPolynomial is None:
482
                ccReducedPolynomialsList.append(ccReducedPolynomial)
483
        else:
484
            #print l2Norm.n() , ">", nAtAlpha
485
            pass
486
        # End if
487
    # End for.
488
    if len(ccReducedPolynomialsList) < 2: # Insufficient reduction.
489
        print "Less than 2 Coppersmith condition compliant vectors."
490
        print "Extra reduction starting..."
491
        reducedMatrixStep2 = reducedMatrixStep1.LLL(algorithm='fpLLL:wrapper')
492
        ## Create a fresh reduced polynomials list.
493
        ccReducedPolynomialsList = []
494
        for row in reducedMatrixStep2.rows():
495
            l2Norm = row.norm(2)
496
            if (l2Norm * monomialsCountSqrt) < nAtAlpha:
497
                #print (l2Norm * monomialsCountSqrt).n()
498
                #print l2Norm.n()
499
                ccReducedPolynomial = \
500
                    slz_compute_reduced_polynomial(row,
501
                                                   knownMonomials,
502
                                                   iVariable,
503
                                                   iBound,
504
                                                   tVariable,
505
                                                   tBound)
506
                if not ccReducedPolynomial is None:
507
                    ccReducedPolynomialsList.append(ccReducedPolynomial)
508
        # End for.
509
    else: # At least 2 small norm enough vectors from reducedMatrixStep2.
510
        print "First step of reduction afforded enough vectors"
511
    return ccReducedPolynomialsList
512
    #print ccReducedPolynomialsList
513
    ## Check again the Coppersmith condition for each row and build the reduced 
514
    #  polynomials.
515
    ccReducedPolynomialsList = []
516
    for row in reducedMatrix.rows():
517
        l2Norm = row.norm(2)
518
        if (l2Norm * monomialsCountSqrt) < nAtAlpha:
519
            #print (l2Norm * monomialsCountSqrt).n()
520
            #print l2Norm.n()
521
            ccReducedPolynomial = \
522
                slz_compute_reduced_polynomial(row,
523
                                               knownMonomials,
524
                                               iVariable,
525
                                               iBound,
526
                                               tVariable,
527
                                               tBound)
528
            if not ccReducedPolynomial is None:
529
                ccReducedPolynomialsList.append(ccReducedPolynomial)
530
        else:
531
            #print l2Norm.n() , ">", nAtAlpha
532
            pass
533
    if len(ccReducedPolynomialsList) < 2: # Insufficient reduction.
534
        print "Less than 2 Coppersmith condition compliant vectors after extra reduction."
535
        return ()
536
    else:
537
        return ccReducedPolynomialsList
538
# End slz_compute_coppersmith_reduced_polynomials_proj
539
#
540
def slz_compute_weak_coppersmith_reduced_polynomials_proj_02(inputPolynomial,
541
                                                             alpha,
542
                                                             N,
543
                                                             iBound,
544
                                                             tBound,
545
                                                             debug = False):
546
    """
547
    For a given set of arguments (see below), compute a list
548
    of "reduced polynomials" that could be used to compute roots
549
    of the inputPolynomial.
550
    INPUT:
551
    
552
    - "inputPolynomial" -- (no default) a bivariate integer polynomial;
553
    - "alpha" -- the alpha parameter of the Coppersmith algorithm;
554
    - "N" -- the modulus;
555
    - "iBound" -- the bound on the first variable;
556
    - "tBound" -- the bound on the second variable.
557
    
558
    OUTPUT:
559
    
560
    A list of bivariate integer polynomial obtained using the Coppersmith
561
    algorithm. The polynomials correspond to the rows of the LLL-reduce
562
    reduced base that comply with the weak version of Coppersmith condition.
563
    """
564
    #@par Changes from runSLZ-113.sage
565
    # LLL reduction is not performed on the matrix itself but rather on the
566
    # product of the matrix with a uniform random matrix.
567
    # The reduced matrix obtained is discarded but the transformation matrix 
568
    # obtained is used to multiply the original matrix in order to reduced it.
569
    # If a sufficient level of reduction is obtained, we stop here. If not
570
    # the product matrix obtained above is LLL reduced. But as it has been
571
    # pre-reduced at the above step, reduction is supposed to be much faster.
572
    #
573
    # Arguments check.
574
    if iBound == 0 or tBound == 0:
575
        return None
576
    # End arguments check.
577
    nAtAlpha = N^alpha
578
    ## Building polynomials for matrix.
579
    polyRing   = inputPolynomial.parent()
580
    # Whatever the 2 variables are actually called, we call them
581
    # 'i' and 't' in all the variable names.
582
    (iVariable, tVariable) = inputPolynomial.variables()[:2]
583
    #print polyVars[0], type(polyVars[0])
584
    initialPolynomial = inputPolynomial.subs({iVariable:iVariable * iBound,
585
                                              tVariable:tVariable * tBound})
586
    if debug:
587
        polynomialsList = \
588
            spo_polynomial_to_polynomials_list_8(initialPolynomial,
589
                                                 alpha,
590
                                                 N,
591
                                                 iBound,
592
                                                 tBound,
593
                                                 20)
594
    else:
595
        polynomialsList = \
596
            spo_polynomial_to_polynomials_list_8(initialPolynomial,
597
                                                 alpha,
598
                                                 N,
599
                                                 iBound,
600
                                                 tBound,
601
                                                 0)
602
    #print "Polynomials list:", polynomialsList
603
    ## Building the proto matrix.
604
    knownMonomials = []
605
    protoMatrix    = []
606
    if debug:
607
        for poly in polynomialsList:
608
            spo_add_polynomial_coeffs_to_matrix_row(poly,
609
                                                    knownMonomials,
610
                                                    protoMatrix,
611
                                                    20)
612
    else:
613
        for poly in polynomialsList:
614
            spo_add_polynomial_coeffs_to_matrix_row(poly,
615
                                                    knownMonomials,
616
                                                    protoMatrix,
617
                                                    0)
618
    matrixToReduce = spo_proto_to_row_matrix(protoMatrix)
619
    #print matrixToReduce
620
    ## Reduction and checking.
621
    ### Reduction with projection
622
    (reducedMatrixStep1, reductionMatrixStep1) = \
623
         slz_reduce_lll_proj_02(matrixToReduce,16)
624
    #print "Reduced matrix:"
625
    #print reducedMatrixStep1
626
    #for row in reducedMatrix.rows():
627
    #    print row
628
    monomialsCount     = len(knownMonomials)
629
    monomialsCountSqrt = sqrt(monomialsCount)
630
    #print "Monomials count:", monomialsCount, monomialsCountSqrt.n()
631
    #print reducedMatrix
632
    ## Check the Coppersmith condition for each row and build the reduced 
633
    #  polynomials.
634
    ccReducedPolynomialsList = []
635
    for row in reducedMatrixStep1.rows():
636
        l1Norm = row.norm(1)
637
        l2Norm = row.norm(2)
638
        if l2Norm * monomialsCountSqrt < l1Norm:
639
            print "l1norm is smaller than l2norm*sqrt(w)."
640
        else:
641
            print "l1norm is NOT smaller than l2norm*sqrt(w)."
642
            print (l2Norm * monomialsCountSqrt).n(), 'vs ', l1Norm.n()
643
        if l1Norm < nAtAlpha:
644
            #print l2Norm.n()
645
            ccReducedPolynomial = \
646
                slz_compute_reduced_polynomial(row,
647
                                               knownMonomials,
648
                                               iVariable,
649
                                               iBound,
650
                                               tVariable,
651
                                               tBound)
652
            if not ccReducedPolynomial is None:
653
                ccReducedPolynomialsList.append(ccReducedPolynomial)
654
        else:
655
            #print l2Norm.n() , ">", nAtAlpha
656
            pass
657
    if len(ccReducedPolynomialsList) < 2: # Insufficient reduction.
658
        print "Less than 2 Coppersmith condition compliant vectors."
659
        print "Extra reduction starting..."
660
        reducedMatrix = reducedMatrixStep1.LLL(algorithm='fpLLL:wrapper')
661
        ### If uncommented, the following statement avoids performing
662
        #   an actual LLL reduction. This allows for demonstrating
663
        #   the behavior of our pseudo-reduction alone.
664
        #return ()
665
    else:
666
        print "First step of reduction afforded enough vectors"
667
        return ccReducedPolynomialsList
668
    #print ccReducedPolynomialsList
669
    ## Check again the Coppersmith condition for each row and build the reduced 
670
    #  polynomials.
671
    ccReducedPolynomialsList = []
672
    for row in reducedMatrix.rows():
673
        l2Norm = row.norm(2)
674
        if (l2Norm * monomialsCountSqrt) < nAtAlpha:
675
            #print (l2Norm * monomialsCountSqrt).n()
676
            #print l2Norm.n()
677
            ccReducedPolynomial = \
678
                slz_compute_reduced_polynomial(row,
679
                                               knownMonomials,
680
                                               iVariable,
681
                                               iBound,
682
                                               tVariable,
683
                                               tBound)
684
            if not ccReducedPolynomial is None:
685
                ccReducedPolynomialsList.append(ccReducedPolynomial)
686
        else:
687
            #print l2Norm.n() , ">", nAtAlpha
688
            pass
689
    if len(ccReducedPolynomialsList) < 2: # Insufficient reduction.
690
        print "Less than 2 Coppersmith condition compliant vectors after extra reduction."
691
        return ()
692
    else:
693
        return ccReducedPolynomialsList
694
# End slz_compute_coppersmith_reduced_polynomials_proj_02
695
#
696
def slz_compute_coppersmith_reduced_polynomials_proj(inputPolynomial,
697
                                                     alpha,
698
                                                     N,
699
                                                     iBound,
700
                                                     tBound,
701
                                                     debug = False):
702
    """
703
    For a given set of arguments (see below), compute a list
704
    of "reduced polynomials" that could be used to compute roots
705
    of the inputPolynomial.
706
    INPUT:
707
    
708
    - "inputPolynomial" -- (no default) a bivariate integer polynomial;
709
    - "alpha" -- the alpha parameter of the Coppersmith algorithm;
710
    - "N" -- the modulus;
711
    - "iBound" -- the bound on the first variable;
712
    - "tBound" -- the bound on the second variable.
713
    
714
    OUTPUT:
715
    
716
    A list of bivariate integer polynomial obtained using the Coppersmith
717
    algorithm. The polynomials correspond to the rows of the LLL-reduce
718
    reduced base that comply with the Coppersmith condition.
719
    """
720
    #@par Changes from runSLZ-113.sage
721
    # LLL reduction is not performed on the matrix itself but rather on the
722
    # product of the matrix with a uniform random matrix.
723
    # The reduced matrix obtained is discarded but the transformation matrix 
724
    # obtained is used to multiply the original matrix in order to reduced it.
725
    # If a sufficient level of reduction is obtained, we stop here. If not
726
    # the product matrix obtained above is LLL reduced. But as it has been
727
    # pre-reduced at the above step, reduction is supposed to be much faster.
728
    #
729
    # Arguments check.
730
    if iBound == 0 or tBound == 0:
731
        return None
732
    # End arguments check.
733
    nAtAlpha = N^alpha
734
    ## Building polynomials for matrix.
735
    polyRing   = inputPolynomial.parent()
736
    # Whatever the 2 variables are actually called, we call them
737
    # 'i' and 't' in all the variable names.
738
    (iVariable, tVariable) = inputPolynomial.variables()[:2]
739
    #print polyVars[0], type(polyVars[0])
740
    initialPolynomial = inputPolynomial.subs({iVariable:iVariable * iBound,
741
                                              tVariable:tVariable * tBound})
742
    if debug:
743
        polynomialsList = \
744
            spo_polynomial_to_polynomials_list_8(initialPolynomial,
745
                                                 alpha,
746
                                                 N,
747
                                                 iBound,
748
                                                 tBound,
749
                                                 20)
750
    else:
751
        polynomialsList = \
752
            spo_polynomial_to_polynomials_list_8(initialPolynomial,
753
                                                 alpha,
754
                                                 N,
755
                                                 iBound,
756
                                                 tBound,
757
                                                 0)
758
    #print "Polynomials list:", polynomialsList
759
    ## Building the proto matrix.
760
    knownMonomials = []
761
    protoMatrix    = []
762
    if debug:
763
        for poly in polynomialsList:
764
            spo_add_polynomial_coeffs_to_matrix_row(poly,
765
                                                    knownMonomials,
766
                                                    protoMatrix,
767
                                                    20)
768
    else:
769
        for poly in polynomialsList:
770
            spo_add_polynomial_coeffs_to_matrix_row(poly,
771
                                                    knownMonomials,
772
                                                    protoMatrix,
773
                                                    0)
774
    matrixToReduce = spo_proto_to_row_matrix(protoMatrix)
775
    #print matrixToReduce
776
    ## Reduction and checking.
777
    ### Reduction with projection
778
    (reducedMatrixStep1, reductionMatrixStep1) = \
779
         slz_reduce_lll_proj(matrixToReduce,16)
780
    #print "Reduced matrix:"
781
    #print reducedMatrixStep1
782
    #for row in reducedMatrix.rows():
783
    #    print row
784
    monomialsCount     = len(knownMonomials)
785
    monomialsCountSqrt = sqrt(monomialsCount)
786
    #print "Monomials count:", monomialsCount, monomialsCountSqrt.n()
787
    #print reducedMatrix
788
    ## Check the Coppersmith condition for each row and build the reduced 
789
    #  polynomials.
790
    ccReducedPolynomialsList = []
791
    for row in reducedMatrixStep1.rows():
792
        l2Norm = row.norm(2)
793
        if (l2Norm * monomialsCountSqrt) < nAtAlpha:
794
            #print (l2Norm * monomialsCountSqrt).n()
795
            #print l2Norm.n()
796
            ccReducedPolynomial = \
797
                slz_compute_reduced_polynomial(row,
798
                                               knownMonomials,
799
                                               iVariable,
800
                                               iBound,
801
                                               tVariable,
802
                                               tBound)
803
            if not ccReducedPolynomial is None:
804
                ccReducedPolynomialsList.append(ccReducedPolynomial)
805
        else:
806
            #print l2Norm.n() , ">", nAtAlpha
807
            pass
808
    if len(ccReducedPolynomialsList) < 2: # Insufficient reduction.
809
        print "Less than 2 Coppersmith condition compliant vectors."
810
        print "Extra reduction starting..."
811
        reducedMatrix = reducedMatrixStep1.LLL(algorithm='fpLLL:wrapper')
812
        ### If uncommented, the following statement avoids performing
813
        #   an actual LLL reduction. This allows for demonstrating
814
        #   the behavior of our pseudo-reduction alone.
815
        #return ()
816
    else:
817
        print "First step of reduction afforded enough vectors"
818
        return ccReducedPolynomialsList
819
    #print ccReducedPolynomialsList
820
    ## Check again the Coppersmith condition for each row and build the reduced 
821
    #  polynomials.
822
    ccReducedPolynomialsList = []
823
    for row in reducedMatrix.rows():
824
        l2Norm = row.norm(2)
825
        if (l2Norm * monomialsCountSqrt) < nAtAlpha:
826
            #print (l2Norm * monomialsCountSqrt).n()
827
            #print l2Norm.n()
828
            ccReducedPolynomial = \
829
                slz_compute_reduced_polynomial(row,
830
                                               knownMonomials,
831
                                               iVariable,
832
                                               iBound,
833
                                               tVariable,
834
                                               tBound)
835
            if not ccReducedPolynomial is None:
836
                ccReducedPolynomialsList.append(ccReducedPolynomial)
837
        else:
838
            #print l2Norm.n() , ">", nAtAlpha
839
            pass
840
    if len(ccReducedPolynomialsList) < 2: # Insufficient reduction.
841
        print "Less than 2 Coppersmith condition compliant vectors after extra reduction."
842
        return ()
843
    else:
844
        return ccReducedPolynomialsList
845
# End slz_compute_coppersmith_reduced_polynomials_proj
846
def slz_compute_weak_coppersmith_reduced_polynomials_proj(inputPolynomial,
847
                                                          alpha,
848
                                                          N,
849
                                                          iBound,
850
                                                          tBound,
851
                                                          debug = False):
852
    """
853
    For a given set of arguments (see below), compute a list
854
    of "reduced polynomials" that could be used to compute roots
855
    of the inputPolynomial.
856
    INPUT:
857
    
858
    - "inputPolynomial" -- (no default) a bivariate integer polynomial;
859
    - "alpha" -- the alpha parameter of the Coppersmith algorithm;
860
    - "N" -- the modulus;
861
    - "iBound" -- the bound on the first variable;
862
    - "tBound" -- the bound on the second variable.
863
    
864
    OUTPUT:
865
    
866
    A list of bivariate integer polynomial obtained using the Coppersmith
867
    algorithm. The polynomials correspond to the rows of the LLL-reduce
868
    reduced base that comply with the weak version of Coppersmith condition.
869
    """
870
    #@par Changes from runSLZ-113.sage
871
    # LLL reduction is not performed on the matrix itself but rather on the
872
    # product of the matrix with a uniform random matrix.
873
    # The reduced matrix obtained is discarded but the transformation matrix 
874
    # obtained is used to multiply the original matrix in order to reduced it.
875
    # If a sufficient level of reduction is obtained, we stop here. If not
876
    # the product matrix obtained above is LLL reduced. But as it has been
877
    # pre-reduced at the above step, reduction is supposed to be much faster.
878
    #
879
    # Arguments check.
880
    if iBound == 0 or tBound == 0:
881
        return None
882
    # End arguments check.
883
    nAtAlpha = N^alpha
884
    ## Building polynomials for matrix.
885
    polyRing   = inputPolynomial.parent()
886
    # Whatever the 2 variables are actually called, we call them
887
    # 'i' and 't' in all the variable names.
888
    (iVariable, tVariable) = inputPolynomial.variables()[:2]
889
    #print polyVars[0], type(polyVars[0])
890
    initialPolynomial = inputPolynomial.subs({iVariable:iVariable * iBound,
891
                                              tVariable:tVariable * tBound})
892
    if debug:
893
        polynomialsList = \
894
            spo_polynomial_to_polynomials_list_8(initialPolynomial,
895
                                                 alpha,
896
                                                 N,
897
                                                 iBound,
898
                                                 tBound,
899
                                                 20)
900
    else:
901
        polynomialsList = \
902
            spo_polynomial_to_polynomials_list_8(initialPolynomial,
903
                                                 alpha,
904
                                                 N,
905
                                                 iBound,
906
                                                 tBound,
907
                                                 0)
908
    #print "Polynomials list:", polynomialsList
909
    ## Building the proto matrix.
910
    knownMonomials = []
911
    protoMatrix    = []
912
    if debug:
913
        for poly in polynomialsList:
914
            spo_add_polynomial_coeffs_to_matrix_row(poly,
915
                                                    knownMonomials,
916
                                                    protoMatrix,
917
                                                    20)
918
    else:
919
        for poly in polynomialsList:
920
            spo_add_polynomial_coeffs_to_matrix_row(poly,
921
                                                    knownMonomials,
922
                                                    protoMatrix,
923
                                                    0)
924
    matrixToReduce = spo_proto_to_row_matrix(protoMatrix)
925
    #print matrixToReduce
926
    ## Reduction and checking.
927
    ### Reduction with projection
928
    (reducedMatrixStep1, reductionMatrixStep1) = \
929
         slz_reduce_lll_proj(matrixToReduce,16)
930
    #print "Reduced matrix:"
931
    #print reducedMatrixStep1
932
    #for row in reducedMatrix.rows():
933
    #    print row
934
    monomialsCount     = len(knownMonomials)
935
    monomialsCountSqrt = sqrt(monomialsCount)
936
    #print "Monomials count:", monomialsCount, monomialsCountSqrt.n()
937
    #print reducedMatrix
938
    ## Check the Coppersmith condition for each row and build the reduced 
939
    #  polynomials.
940
    ccReducedPolynomialsList = []
941
    for row in reducedMatrixStep1.rows():
942
        l1Norm = row.norm(1)
943
        l2Norm = row.norm(2)
944
        if l2Norm * monomialsCountSqrt < l1Norm:
945
            print "l1norm is smaller than l2norm*sqrt(w)."
946
        else:
947
            print "l1norm is NOT smaller than l2norm*sqrt(w)."
948
            print (l2Norm * monomialsCountSqrt).n(), 'vs ', l1Norm.n()
949
        if l1Norm < nAtAlpha:
950
            #print l2Norm.n()
951
            ccReducedPolynomial = \
952
                slz_compute_reduced_polynomial(row,
953
                                               knownMonomials,
954
                                               iVariable,
955
                                               iBound,
956
                                               tVariable,
957
                                               tBound)
958
            if not ccReducedPolynomial is None:
959
                ccReducedPolynomialsList.append(ccReducedPolynomial)
960
        else:
961
            #print l2Norm.n() , ">", nAtAlpha
962
            pass
963
    if len(ccReducedPolynomialsList) < 2: # Insufficient reduction.
964
        print "Less than 2 Coppersmith condition compliant vectors."
965
        print "Extra reduction starting..."
966
        reducedMatrix = reducedMatrixStep1.LLL(algorithm='fpLLL:wrapper')
967
        ### If uncommented, the following statement avoids performing
968
        #   an actual LLL reduction. This allows for demonstrating
969
        #   the behavior of our pseudo-reduction alone.
970
        #return ()
971
    else:
972
        print "First step of reduction afforded enough vectors"
973
        return ccReducedPolynomialsList
974
    #print ccReducedPolynomialsList
975
    ## Check again the Coppersmith condition for each row and build the reduced 
976
    #  polynomials.
977
    ccReducedPolynomialsList = []
978
    for row in reducedMatrix.rows():
979
        l2Norm = row.norm(2)
980
        if (l2Norm * monomialsCountSqrt) < nAtAlpha:
981
            #print (l2Norm * monomialsCountSqrt).n()
982
            #print l2Norm.n()
983
            ccReducedPolynomial = \
984
                slz_compute_reduced_polynomial(row,
985
                                               knownMonomials,
986
                                               iVariable,
987
                                               iBound,
988
                                               tVariable,
989
                                               tBound)
990
            if not ccReducedPolynomial is None:
991
                ccReducedPolynomialsList.append(ccReducedPolynomial)
992
        else:
993
            #print l2Norm.n() , ">", nAtAlpha
994
            pass
995
    if len(ccReducedPolynomialsList) < 2: # Insufficient reduction.
996
        print "Less than 2 Coppersmith condition compliant vectors after extra reduction."
997
        return ()
998
    else:
999
        return ccReducedPolynomialsList
1000
# End slz_compute_coppersmith_reduced_polynomials_proj2
1001
#
1002
def slz_compute_coppersmith_reduced_polynomials_with_lattice_volume(inputPolynomial,
1003
                                                                    alpha,
1004
                                                                    N,
1005
                                                                    iBound,
1006
                                                                    tBound,
1007
                                                                    debug = False):
1008
    """
1009
    For a given set of arguments (see below), compute a list
1010
    of "reduced polynomials" that could be used to compute roots
1011
    of the inputPolynomial.
1012
    Print the volume of the initial basis as well.
1013
    INPUT:
1014
    
1015
    - "inputPolynomial" -- (no default) a bivariate integer polynomial;
1016
    - "alpha" -- the alpha parameter of the Coppersmith algorithm;
1017
    - "N" -- the modulus;
1018
    - "iBound" -- the bound on the first variable;
1019
    - "tBound" -- the bound on the second variable.
1020
    
1021
    OUTPUT:
1022
    
1023
    A list of bivariate integer polynomial obtained using the Coppersmith
1024
    algorithm. The polynomials correspond to the rows of the LLL-reduce
1025
    reduced base that comply with the Coppersmith condition.
1026
    """
1027
    # Arguments check.
1028
    if iBound == 0 or tBound == 0:
1029
        return None
1030
    # End arguments check.
1031
    nAtAlpha = N^alpha
1032
    if debug:
1033
        print "N at alpha:", nAtAlpha
1034
    ## Building polynomials for matrix.
1035
    polyRing   = inputPolynomial.parent()
1036
    # Whatever the 2 variables are actually called, we call them
1037
    # 'i' and 't' in all the variable names.
1038
    (iVariable, tVariable) = inputPolynomial.variables()[:2]
1039
    #print polyVars[0], type(polyVars[0])
1040
    initialPolynomial = inputPolynomial.subs({iVariable:iVariable * iBound,
1041
                                              tVariable:tVariable * tBound})
1042
##    polynomialsList = \
1043
##        spo_polynomial_to_polynomials_list_8(initialPolynomial,
1044
##        spo_polynomial_to_polynomials_list_5(initialPolynomial,
1045
    polynomialsList = \
1046
        spo_polynomial_to_polynomials_list_5(initialPolynomial,
1047
                                             alpha,
1048
                                             N,
1049
                                             iBound,
1050
                                             tBound,
1051
                                             0)
1052
    #print "Polynomials list:", polynomialsList
1053
    ## Building the proto matrix.
1054
    knownMonomials = []
1055
    protoMatrix    = []
1056
    for poly in polynomialsList:
1057
        spo_add_polynomial_coeffs_to_matrix_row(poly,
1058
                                                knownMonomials,
1059
                                                protoMatrix,
1060
                                                0)
1061
    matrixToReduce = spo_proto_to_row_matrix(protoMatrix)
1062
    matrixToReduceTranspose = matrixToReduce.transpose()
1063
    squareMatrix = matrixToReduce * matrixToReduceTranspose
1064
    squareMatDet = det(squareMatrix)
1065
    latticeVolume = sqrt(squareMatDet)
1066
    print "Lattice volume:", latticeVolume.n()
1067
    print "Lattice volume / N:", (latticeVolume/N).n()
1068
    #print matrixToReduce
1069
    ## Reduction and checking.
1070
    ## S.T. changed 'fp' to None as of Sage 6.6 complying to
1071
    #  error message issued when previous code was used.
1072
    #reducedMatrix = matrixToReduce.LLL(fp='fp')
1073
    reductionTimeStart = cputime() 
1074
    reducedMatrix = matrixToReduce.LLL(fp=None)
1075
    reductionTime = cputime(reductionTimeStart)
1076
    print "Reduction time:", reductionTime
1077
    isLLLReduced  = reducedMatrix.is_LLL_reduced()
1078
    if not isLLLReduced:
1079
       return None
1080
    #
1081
    if debug:
1082
        matrixFile = file('/tmp/reducedMatrix.txt', 'w')
1083
        for row in reducedMatrix.rows():
1084
            matrixFile.write(str(row) + "\n")
1085
        matrixFile.close()
1086
    #
1087
    monomialsCount     = len(knownMonomials)
1088
    monomialsCountSqrt = sqrt(monomialsCount)
1089
    #print "Monomials count:", monomialsCount, monomialsCountSqrt.n()
1090
    #print reducedMatrix
1091
    ## Check the Coppersmith condition for each row and build the reduced 
1092
    #  polynomials.
1093
    ccVectorsCount = 0
1094
    ccReducedPolynomialsList = []
1095
    for row in reducedMatrix.rows():
1096
        l2Norm = row.norm(2)
1097
        if (l2Norm * monomialsCountSqrt) < nAtAlpha:
1098
            #print (l2Norm * monomialsCountSqrt).n()
1099
            #print l2Norm.n()
1100
            ccVectorsCount +=1            
1101
            ccReducedPolynomial = \
1102
                slz_compute_reduced_polynomial(row,
1103
                                               knownMonomials,
1104
                                               iVariable,
1105
                                               iBound,
1106
                                               tVariable,
1107
                                               tBound)
1108
            if not ccReducedPolynomial is None:
1109
                ccReducedPolynomialsList.append(ccReducedPolynomial)
1110
        else:
1111
            #print l2Norm.n() , ">", nAtAlpha
1112
            pass
1113
    if debug:
1114
        print ccVectorsCount, "out of ", len(ccReducedPolynomialsList), 
1115
        print "took Coppersmith text."
1116
    if len(ccReducedPolynomialsList) < 2:
1117
        print "Less than 2 Coppersmith condition compliant vectors."
1118
        return ()
1119
    if debug:
1120
        print "Reduced and Coppersmith compliant polynomials list", ccReducedPolynomialsList
1121
    return ccReducedPolynomialsList
1122
# End slz_compute_coppersmith_reduced_polynomials_with_lattice volume
1123

    
1124
def slz_compute_initial_lattice_matrix(inputPolynomial,
1125
                                       alpha,
1126
                                       N,
1127
                                       iBound,
1128
                                       tBound,
1129
                                       debug = False):
1130
    """
1131
    For a given set of arguments (see below), compute the initial lattice
1132
    that could be reduced. 
1133
    INPUT:
1134
    
1135
    - "inputPolynomial" -- (no default) a bivariate integer polynomial;
1136
    - "alpha" -- the alpha parameter of the Coppersmith algorithm;
1137
    - "N" -- the modulus;
1138
    - "iBound" -- the bound on the first variable;
1139
    - "tBound" -- the bound on the second variable.
1140
    
1141
    OUTPUT:
1142
    
1143
    A list of bivariate integer polynomial obtained using the Coppersmith
1144
    algorithm. The polynomials correspond to the rows of the LLL-reduce
1145
    reduced base that comply with the Coppersmith condition.
1146
    """
1147
    # Arguments check.
1148
    if iBound == 0 or tBound == 0:
1149
        return None
1150
    # End arguments check.
1151
    nAtAlpha = N^alpha
1152
    ## Building polynomials for matrix.
1153
    polyRing   = inputPolynomial.parent()
1154
    # Whatever the 2 variables are actually called, we call them
1155
    # 'i' and 't' in all the variable names.
1156
    (iVariable, tVariable) = inputPolynomial.variables()[:2]
1157
    #print polyVars[0], type(polyVars[0])
1158
    initialPolynomial = inputPolynomial.subs({iVariable:iVariable * iBound,
1159
                                              tVariable:tVariable * tBound})
1160
    polynomialsList = \
1161
        spo_polynomial_to_polynomials_list_8(initialPolynomial,
1162
                                             alpha,
1163
                                             N,
1164
                                             iBound,
1165
                                             tBound,
1166
                                             0)
1167
    #print "Polynomials list:", polynomialsList
1168
    ## Building the proto matrix.
1169
    knownMonomials = []
1170
    protoMatrix    = []
1171
    for poly in polynomialsList:
1172
        spo_add_polynomial_coeffs_to_matrix_row(poly,
1173
                                                knownMonomials,
1174
                                                protoMatrix,
1175
                                                0)
1176
    matrixToReduce = spo_proto_to_row_matrix(protoMatrix)
1177
    if debug:
1178
        print "Initial basis polynomials"
1179
        for poly in polynomialsList:
1180
            print poly
1181
    return matrixToReduce
1182
# End slz_compute_initial_lattice_matrix.
1183

    
1184
def slz_compute_integer_polynomial_modular_roots(inputPolynomial,
1185
                                                 alpha,
1186
                                                 N,
1187
                                                 iBound,
1188
                                                 tBound):
1189
    """
1190
    For a given set of arguments (see below), compute the polynomial modular 
1191
    roots, if any.
1192
    
1193
    """
1194
    # Arguments check.
1195
    if iBound == 0 or tBound == 0:
1196
        return set()
1197
    # End arguments check.
1198
    nAtAlpha = N^alpha
1199
    ## Building polynomials for matrix.
1200
    polyRing   = inputPolynomial.parent()
1201
    # Whatever the 2 variables are actually called, we call them
1202
    # 'i' and 't' in all the variable names.
1203
    (iVariable, tVariable) = inputPolynomial.variables()[:2]
1204
    ccReducedPolynomialsList = \
1205
        slz_compute_coppersmith_reduced_polynomials (inputPolynomial,
1206
                                                     alpha,
1207
                                                     N,
1208
                                                     iBound,
1209
                                                     tBound)
1210
    if len(ccReducedPolynomialsList) == 0:
1211
        return set()   
1212
    ## Create the valid (poly1 and poly2 are algebraically independent) 
1213
    # resultant tuples (poly1, poly2, resultant(poly1, poly2)).
1214
    # Try to mix and match all the polynomial pairs built from the 
1215
    # ccReducedPolynomialsList to obtain non zero resultants.
1216
    resultantsInITuplesList = []
1217
    for polyOuterIndex in xrange(0, len(ccReducedPolynomialsList)-1):
1218
        for polyInnerIndex in xrange(polyOuterIndex+1, 
1219
                                     len(ccReducedPolynomialsList)):
1220
            # Compute the resultant in resultants in the
1221
            # first variable (is it the optimal choice?).
1222
            resultantInI = \
1223
               ccReducedPolynomialsList[polyOuterIndex].resultant(ccReducedPolynomialsList[polyInnerIndex],
1224
                                                                  ccReducedPolynomialsList[0].parent(str(iVariable)))
1225
            #print "Resultant", resultantInI
1226
            # Test algebraic independence.
1227
            if not resultantInI.is_zero():
1228
                resultantsInITuplesList.append((ccReducedPolynomialsList[polyOuterIndex],
1229
                                                 ccReducedPolynomialsList[polyInnerIndex],
1230
                                                 resultantInI))
1231
    # If no non zero resultant was found: we can't get no algebraically 
1232
    # independent polynomials pair. Give up!
1233
    if len(resultantsInITuplesList) == 0:
1234
        return set()
1235
    #print resultantsInITuplesList
1236
    # Compute the roots.
1237
    Zi = ZZ[str(iVariable)]
1238
    Zt = ZZ[str(tVariable)]
1239
    polynomialRootsSet = set()
1240
    # First, solve in the second variable since resultants are in the first
1241
    # variable.
1242
    for resultantInITuple in resultantsInITuplesList:
1243
        tRootsList = Zt(resultantInITuple[2]).roots()
1244
        # For each tRoot, compute the corresponding iRoots and check 
1245
        # them in the input polynomial.
1246
        for tRoot in tRootsList:
1247
            #print "tRoot:", tRoot
1248
            # Roots returned by root() are (value, multiplicity) tuples.
1249
            iRootsList = \
1250
                Zi(resultantInITuple[0].subs({resultantInITuple[0].variables()[1]:tRoot[0]})).roots()
1251
            print iRootsList
1252
            # The iRootsList can be empty, hence the test.
1253
            if len(iRootsList) != 0:
1254
                for iRoot in iRootsList:
1255
                    polyEvalModN = inputPolynomial(iRoot[0], tRoot[0]) / N
1256
                    # polyEvalModN must be an integer.
1257
                    if polyEvalModN == int(polyEvalModN):
1258
                        polynomialRootsSet.add((iRoot[0],tRoot[0]))
1259
    return polynomialRootsSet
1260
# End slz_compute_integer_polynomial_modular_roots.
1261
#
1262
def slz_compute_integer_polynomial_modular_roots_2(inputPolynomial,
1263
                                                 alpha,
1264
                                                 N,
1265
                                                 iBound,
1266
                                                 tBound):
1267
    """
1268
    For a given set of arguments (see below), compute the polynomial modular 
1269
    roots, if any.
1270
    This version differs in the way resultants are computed.   
1271
    """
1272
    # Arguments check.
1273
    if iBound == 0 or tBound == 0:
1274
        return set()
1275
    # End arguments check.
1276
    nAtAlpha = N^alpha
1277
    ## Building polynomials for matrix.
1278
    polyRing   = inputPolynomial.parent()
1279
    # Whatever the 2 variables are actually called, we call them
1280
    # 'i' and 't' in all the variable names.
1281
    (iVariable, tVariable) = inputPolynomial.variables()[:2]
1282
    #print polyVars[0], type(polyVars[0])
1283
    ccReducedPolynomialsList = \
1284
        slz_compute_coppersmith_reduced_polynomials (inputPolynomial,
1285
                                                     alpha,
1286
                                                     N,
1287
                                                     iBound,
1288
                                                     tBound)
1289
    if len(ccReducedPolynomialsList) == 0:
1290
        return set()   
1291
    ## Create the valid (poly1 and poly2 are algebraically independent) 
1292
    # resultant tuples (poly1, poly2, resultant(poly1, poly2)).
1293
    # Try to mix and match all the polynomial pairs built from the 
1294
    # ccReducedPolynomialsList to obtain non zero resultants.
1295
    resultantsInTTuplesList = []
1296
    for polyOuterIndex in xrange(0, len(ccReducedPolynomialsList)-1):
1297
        for polyInnerIndex in xrange(polyOuterIndex+1, 
1298
                                     len(ccReducedPolynomialsList)):
1299
            # Compute the resultant in resultants in the
1300
            # first variable (is it the optimal choice?).
1301
            resultantInT = \
1302
               ccReducedPolynomialsList[polyOuterIndex].resultant(ccReducedPolynomialsList[polyInnerIndex],
1303
                                                                  ccReducedPolynomialsList[0].parent(str(tVariable)))
1304
            #print "Resultant", resultantInT
1305
            # Test algebraic independence.
1306
            if not resultantInT.is_zero():
1307
                resultantsInTTuplesList.append((ccReducedPolynomialsList[polyOuterIndex],
1308
                                                 ccReducedPolynomialsList[polyInnerIndex],
1309
                                                 resultantInT))
1310
    # If no non zero resultant was found: we can't get no algebraically 
1311
    # independent polynomials pair. Give up!
1312
    if len(resultantsInTTuplesList) == 0:
1313
        return set()
1314
    #print resultantsInITuplesList
1315
    # Compute the roots.
1316
    Zi = ZZ[str(iVariable)]
1317
    Zt = ZZ[str(tVariable)]
1318
    polynomialRootsSet = set()
1319
    # First, solve in the second variable since resultants are in the first
1320
    # variable.
1321
    for resultantInTTuple in resultantsInTTuplesList:
1322
        iRootsList = Zi(resultantInTTuple[2]).roots()
1323
        # For each iRoot, compute the corresponding tRoots and check 
1324
        # them in the input polynomial.
1325
        for iRoot in iRootsList:
1326
            #print "iRoot:", iRoot
1327
            # Roots returned by root() are (value, multiplicity) tuples.
1328
            tRootsList = \
1329
                Zt(resultantInTTuple[0].subs({resultantInTTuple[0].variables()[0]:iRoot[0]})).roots()
1330
            print tRootsList
1331
            # The tRootsList can be empty, hence the test.
1332
            if len(tRootsList) != 0:
1333
                for tRoot in tRootsList:
1334
                    polyEvalModN = inputPolynomial(iRoot[0],tRoot[0]) / N
1335
                    # polyEvalModN must be an integer.
1336
                    if polyEvalModN == int(polyEvalModN):
1337
                        polynomialRootsSet.add((iRoot[0],tRoot[0]))
1338
    return polynomialRootsSet
1339
# End slz_compute_integer_polynomial_modular_roots_2.
1340
#
1341
def slz_compute_polynomial_and_interval(functionSo, degreeSo, lowerBoundSa, 
1342
                                        upperBoundSa, approxAccurSa, 
1343
                                        precSa=None, debug=False):
1344
    """
1345
    Under the assumptions listed for slz_get_intervals_and_polynomials, compute
1346
    a polynomial that approximates the function on a an interval starting
1347
    at lowerBoundSa and finishing at a value that guarantees that the polynomial
1348
    approximates with the expected precision.
1349
    The interval upper bound is lowered until the expected approximation 
1350
    precision is reached.
1351
    The polynomial, the bounds, the center of the interval and the error
1352
    are returned.
1353
    Argument debug is not used.
1354
    OUTPUT:
1355
    A tuple made of 4 Sollya objects:
1356
    - a polynomial;
1357
    - an range (an interval, not in the sense of number given as an interval);
1358
    - the center of the interval;
1359
    - the maximum error in the approximation of the input functionSo by the
1360
      output polynomial ; this error <= approxAccurSaS.
1361
    
1362
    """
1363
    #print"In slz_compute_polynomial_and_interval..."
1364
    ## Superficial argument check.
1365
    if lowerBoundSa > upperBoundSa:
1366
        return None
1367
    ## Change Sollya precision, if requested.
1368
    if precSa is None:
1369
        precSa = ceil((RR('1.5') * abs(RR(approxAccurSa).log2())) / 64) * 64
1370
        #print "Computed internal precision:", precSa
1371
        if precSa < 192:
1372
            precSa = 192
1373
    sollyaPrecChanged = False
1374
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
1375
    if precSa > initialSollyaPrecSa:
1376
        if precSa <= 2:
1377
            print inspect.stack()[0][3], ": precision change <=2 requested."
1378
        pobyso_set_prec_sa_so(precSa)
1379
        sollyaPrecChanged = True
1380
    RRR = lowerBoundSa.parent()
1381
    intervalShrinkConstFactorSa = RRR('0.9')
1382
    absoluteErrorTypeSo = pobyso_absolute_so_so()
1383
    currentRangeSo = pobyso_bounds_to_range_sa_so(lowerBoundSa, upperBoundSa)
1384
    currentUpperBoundSa = upperBoundSa
1385
    currentLowerBoundSa = lowerBoundSa
1386
    # What we want here is the polynomial without the variable change, 
1387
    # since our actual variable will be x-intervalCenter defined over the 
1388
    # domain [lowerBound-intervalCenter , upperBound-intervalCenter]. 
1389
    (polySo, intervalCenterSo, maxErrorSo) = \
1390
        pobyso_taylor_expansion_no_change_var_so_so(functionSo, degreeSo, 
1391
                                                    currentRangeSo, 
1392
                                                    absoluteErrorTypeSo)
1393
    maxErrorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo)
1394
    while maxErrorSa > approxAccurSa:
1395
        print "++Approximation error:", maxErrorSa.n()
1396
        sollya_lib_clear_obj(polySo)
1397
        sollya_lib_clear_obj(intervalCenterSo)
1398
        sollya_lib_clear_obj(maxErrorSo)
1399
        # Very empirical shrinking factor.
1400
        shrinkFactorSa = 1 /  (maxErrorSa/approxAccurSa).log2().abs()
1401
        print "Shrink factor:", \
1402
              shrinkFactorSa.n(), \
1403
              intervalShrinkConstFactorSa
1404
        print
1405
        #errorRatioSa = approxAccurSa/maxErrorSa
1406
        #print "Error ratio: ", errorRatioSa
1407
        # Make sure interval shrinks.
1408
        if shrinkFactorSa > intervalShrinkConstFactorSa:
1409
            actualShrinkFactorSa = intervalShrinkConstFactorSa
1410
            #print "Fixed"
1411
        else:
1412
            actualShrinkFactorSa = shrinkFactorSa
1413
            #print "Computed",shrinkFactorSa,maxErrorSa
1414
            #print shrinkFactorSa, maxErrorSa
1415
        #print "Shrink factor", actualShrinkFactorSa
1416
        currentUpperBoundSa = currentLowerBoundSa + \
1417
                                (currentUpperBoundSa - currentLowerBoundSa) * \
1418
                                actualShrinkFactorSa
1419
        #print "Current upper bound:", currentUpperBoundSa
1420
        sollya_lib_clear_obj(currentRangeSo)
1421
        # Check what is left with the bounds.
1422
        if currentUpperBoundSa <= currentLowerBoundSa or \
1423
              currentUpperBoundSa == currentLowerBoundSa.nextabove():
1424
            sollya_lib_clear_obj(absoluteErrorTypeSo)
1425
            print "Can't find an interval."
1426
            print "Use either or both a higher polynomial degree or a higher",
1427
            print "internal precision."
1428
            print "Aborting!"
1429
            if sollyaPrecChanged:
1430
                pobyso_set_prec_so_so(initialSollyaPrecSo)
1431
            sollya_lib_clear_obj(initialSollyaPrecSo)
1432
            return None
1433
        currentRangeSo = pobyso_bounds_to_range_sa_so(currentLowerBoundSa, 
1434
                                                      currentUpperBoundSa)
1435
        # print "New interval:",
1436
        # pobyso_autoprint(currentRangeSo)
1437
        #print "Second Taylor expansion call."
1438
        (polySo, intervalCenterSo, maxErrorSo) = \
1439
            pobyso_taylor_expansion_no_change_var_so_so(functionSo, degreeSo, 
1440
                                                        currentRangeSo, 
1441
                                                        absoluteErrorTypeSo)
1442
        #maxErrorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo, RRR)
1443
        #print "Max errorSo:",
1444
        #pobyso_autoprint(maxErrorSo)
1445
        maxErrorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo)
1446
        #print "Max errorSa:", maxErrorSa
1447
        #print "Sollya prec:",
1448
        #pobyso_autoprint(sollya_lib_get_prec(None))
1449
    # End while
1450
    sollya_lib_clear_obj(absoluteErrorTypeSo)
1451
    if sollyaPrecChanged:
1452
        pobyso_set_prec_so_so(initialSollyaPrecSo)
1453
    sollya_lib_clear_obj(initialSollyaPrecSo)
1454
    return (polySo, currentRangeSo, intervalCenterSo, maxErrorSo)
1455
# End slz_compute_polynomial_and_interval
1456

    
1457
def slz_compute_polynomial_and_interval_01(functionSo, degreeSo, lowerBoundSa, 
1458
                                        upperBoundSa, approxAccurSa, 
1459
                                        sollyaPrecSa=None, debug=False):
1460
    """
1461
    Under the assumptions listed for slz_get_intervals_and_polynomials, compute
1462
    a polynomial that approximates the function on a an interval starting
1463
    at lowerBoundSa and finishing at a value that guarantees that the polynomial
1464
    approximates with the expected precision.
1465
    The interval upper bound is lowered until the expected approximation 
1466
    precision is reached.
1467
    The polynomial, the bounds, the center of the interval and the error
1468
    are returned.
1469
    OUTPUT:
1470
    A tuple made of 4 Sollya objects:
1471
    - a polynomial;
1472
    - an range (an interval, not in the sense of number given as an interval);
1473
    - the center of the interval;
1474
    - the maximum error in the approximation of the input functionSo by the
1475
      output polynomial ; this error <= approxAccurSaS.
1476
      
1477
    Changes from version 0 (no number):
1478
    - make use of debug argument;
1479
    - polynomial coefficients are "shaved".
1480
    
1481
    """
1482
    #print"In slz_compute_polynomial_and_interval..."
1483
    ## Superficial argument check.
1484
    if lowerBoundSa > upperBoundSa:
1485
        print inspect.stack()[0][3], ": lower bound is larger than upper bound. "
1486
        return None
1487
    ## Change Sollya precision, if requested.
1488
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
1489
    sollyaPrecChangedSa = False
1490
    if sollyaPrecSa is None:
1491
        sollyaPrecSa = initialSollyaPrecSa
1492
    else:
1493
        if sollyaPrecSa > initialSollyaPrecSa:
1494
            if sollyaPrecSa <= 2:
1495
                print inspect.stack()[0][3], ": precision change <= 2 requested."
1496
            pobyso_set_prec_sa_so(sollyaPrecSa)
1497
            sollyaPrecChangedSa = True
1498
    ## Other initializations and data recovery.
1499
    RRR = lowerBoundSa.parent()
1500
    intervalShrinkConstFactorSa = RRR('0.9')
1501
    absoluteErrorTypeSo = pobyso_absolute_so_so()
1502
    currentRangeSo = pobyso_bounds_to_range_sa_so(lowerBoundSa, upperBoundSa)
1503
    currentUpperBoundSa = upperBoundSa
1504
    currentLowerBoundSa = lowerBoundSa
1505
    # What we want here is the polynomial without the variable change, 
1506
    # since our actual variable will be x-intervalCenter defined over the 
1507
    # domain [lowerBound-intervalCenter , upperBound-intervalCenter]. 
1508
    (polySo, intervalCenterSo, maxErrorSo) = \
1509
        pobyso_taylor_expansion_no_change_var_so_so(functionSo, degreeSo, 
1510
                                                    currentRangeSo, 
1511
                                                    absoluteErrorTypeSo)
1512
    maxErrorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo)
1513
    ## Expression "approxAccurSa/2" (instead of just approxAccurSa) is added 
1514
    #  since we use polynomial shaving. In order for it to work, we must start 
1515
    #  with a better accuracy.
1516
    while maxErrorSa > approxAccurSa/2:
1517
        print "++Approximation error:", maxErrorSa.n()
1518
        sollya_lib_clear_obj(polySo)
1519
        sollya_lib_clear_obj(intervalCenterSo)
1520
        sollya_lib_clear_obj(maxErrorSo)
1521
        # Very empirical shrinking factor.
1522
        shrinkFactorSa = 1 /  (maxErrorSa/approxAccurSa).log2().abs()
1523
        print "Shrink factor:", \
1524
              shrinkFactorSa.n(), \
1525
              intervalShrinkConstFactorSa
1526
        print
1527
        #errorRatioSa = approxAccurSa/maxErrorSa
1528
        #print "Error ratio: ", errorRatioSa
1529
        # Make sure interval shrinks.
1530
        if shrinkFactorSa > intervalShrinkConstFactorSa:
1531
            actualShrinkFactorSa = intervalShrinkConstFactorSa
1532
            #print "Fixed"
1533
        else:
1534
            actualShrinkFactorSa = shrinkFactorSa
1535
            #print "Computed",shrinkFactorSa,maxErrorSa
1536
            #print shrinkFactorSa, maxErrorSa
1537
        #print "Shrink factor", actualShrinkFactorSa
1538
        currentUpperBoundSa = currentLowerBoundSa + \
1539
                                (currentUpperBoundSa - currentLowerBoundSa) * \
1540
                                actualShrinkFactorSa
1541
        #print "Current upper bound:", currentUpperBoundSa
1542
        sollya_lib_clear_obj(currentRangeSo)
1543
        # Check what is left with the bounds.
1544
        if currentUpperBoundSa <= currentLowerBoundSa or \
1545
              currentUpperBoundSa == currentLowerBoundSa.nextabove():
1546
            sollya_lib_clear_obj(absoluteErrorTypeSo)
1547
            print "Can't find an interval."
1548
            print "Use either or both a higher polynomial degree or a higher",
1549
            print "internal precision."
1550
            print "Aborting!"
1551
            if sollyaPrecChangedSa:
1552
                pobyso_set_prec_so_so(initialSollyaPrecSo)
1553
            sollya_lib_clear_obj(initialSollyaPrecSo)
1554
            return None
1555
        currentRangeSo = pobyso_bounds_to_range_sa_so(currentLowerBoundSa, 
1556
                                                      currentUpperBoundSa)
1557
        # print "New interval:",
1558
        # pobyso_autoprint(currentRangeSo)
1559
        #print "Second Taylor expansion call."
1560
        (polySo, intervalCenterSo, maxErrorSo) = \
1561
            pobyso_taylor_expansion_no_change_var_so_so(functionSo, degreeSo, 
1562
                                                        currentRangeSo, 
1563
                                                        absoluteErrorTypeSo)
1564
        #maxErrorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo, RRR)
1565
        #print "Max errorSo:",
1566
        #pobyso_autoprint(maxErrorSo)
1567
        maxErrorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo)
1568
        #print "Max errorSa:", maxErrorSa
1569
        #print "Sollya prec:",
1570
        #pobyso_autoprint(sollya_lib_get_prec(None))
1571
    # End while
1572
    sollya_lib_clear_obj(absoluteErrorTypeSo)
1573
    itpSo           = pobyso_constant_from_int_sa_so(floor(sollyaPrecSa/3))
1574
    ftpSo           = pobyso_constant_from_int_sa_so(floor(2*sollyaPrecSa/3))
1575
    maxPrecSo       = pobyso_constant_from_int_sa_so(sollyaPrecSa)
1576
    approxAccurSo   = pobyso_constant_sa_so(RR(approxAccurSa))
1577
    if debug:
1578
        print inspect.stack()[0][3], "SollyaPrecSa:", sollyaPrecSa
1579
        print "About to call polynomial rounding with:"
1580
        print "polySo: ",           ; pobyso_autoprint(polySo)
1581
        print "functionSo: ",       ; pobyso_autoprint(functionSo)
1582
        print "intervalCenterSo: ", ; pobyso_autoprint(intervalCenterSo)
1583
        print "currentRangeSo: ",   ; pobyso_autoprint(currentRangeSo)
1584
        print "itpSo: ",            ; pobyso_autoprint(itpSo)
1585
        print "ftpSo: ",            ; pobyso_autoprint(ftpSo)
1586
        print "maxPrecSo: ",        ; pobyso_autoprint(maxPrecSo)
1587
        print "approxAccurSo: ",    ; pobyso_autoprint(approxAccurSo)
1588
    """
1589
    # Naive rounding.
1590
    (roundedPolySo, roundedPolyMaxErrSo) = \
1591
    pobyso_polynomial_coefficients_progressive_round_so_so(polySo,
1592
                                                           functionSo,
1593
                                                           intervalCenterSo,
1594
                                                           currentRangeSo,
1595
                                                           itpSo,
1596
                                                           ftpSo,
1597
                                                           maxPrecSo,
1598
                                                           approxAccurSo)
1599
    """
1600
    # Proved rounding.
1601
    (roundedPolySo, roundedPolyMaxErrSo) = \
1602
        pobyso_round_coefficients_progressive_so_so(polySo, 
1603
                                                    functionSo,
1604
                                                    maxPrecSo,
1605
                                                    currentRangeSo,
1606
                                                    intervalCenterSo,
1607
                                                    maxErrorSo,
1608
                                                    approxAccurSo,
1609
                                                    debug=False)
1610
    #### Comment out the two next lines when polynomial rounding is activated.
1611
    #roundedPolySo = sollya_lib_copy_obj(polySo)
1612
    #roundedPolyMaxErrSo =  sollya_lib_copy_obj(maxErrorSo)
1613
    sollya_lib_clear_obj(polySo)
1614
    sollya_lib_clear_obj(maxErrorSo)
1615
    sollya_lib_clear_obj(itpSo)
1616
    sollya_lib_clear_obj(ftpSo)
1617
    sollya_lib_clear_obj(approxAccurSo)
1618
    if sollyaPrecChangedSa:
1619
        pobyso_set_prec_so_so(initialSollyaPrecSo)
1620
    sollya_lib_clear_obj(initialSollyaPrecSo)
1621
    if debug:
1622
        print "1: ", ; pobyso_autoprint(roundedPolySo)
1623
        print "2: ", ; pobyso_autoprint(currentRangeSo)
1624
        print "3: ", ; pobyso_autoprint(intervalCenterSo)
1625
        print "4: ", ; pobyso_autoprint(roundedPolyMaxErrSo)
1626
    return (roundedPolySo, currentRangeSo, intervalCenterSo, roundedPolyMaxErrSo)
1627
# End slz_compute_polynomial_and_interval_01
1628

    
1629
def slz_compute_polynomial_and_interval_02(functionSo, degreeSo, lowerBoundSa, 
1630
                                        upperBoundSa, approxAccurSa, 
1631
                                        sollyaPrecSa=None, debug=True ):
1632
    """
1633
    Under the assumptions listed for slz_get_intervals_and_polynomials, compute
1634
    a polynomial that approximates the function on a an interval starting
1635
    at lowerBoundSa and finishing at a value that guarantees that the polynomial
1636
    approximates with the expected precision.
1637
    The interval upper bound is lowered until the expected approximation 
1638
    precision is reached.
1639
    The polynomial, the bounds, the center of the interval and the error
1640
    are returned.
1641
    OUTPUT:
1642
    A tuple made of 4 Sollya objects:
1643
    - a polynomial;
1644
    - an range (an interval, not in the sense of number given as an interval);
1645
    - the center of the interval;
1646
    - the maximum error in the approximation of the input functionSo by the
1647
      output polynomial ; this error <= approxAccurSaS.
1648
      
1649
    Changes fom v 01:
1650
        extra verbose.
1651
    """
1652
    print"In slz_compute_polynomial_and_interval..."
1653
    ## Superficial argument check.
1654
    if lowerBoundSa > upperBoundSa:
1655
        return None
1656
    ## Change Sollya precision, if requested.
1657
    sollyaPrecChanged = False
1658
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
1659
    #print "Initial Sollya prec:", initialSollyaPrecSa, type(initialSollyaPrecSa)
1660
    if sollyaPrecSa is None:
1661
        sollyaPrecSa = initialSollyaPrecSa
1662
    else:
1663
        if sollyaPrecSa <= 2:
1664
            print inspect.stack()[0][3], ": precision change <=2 requested."
1665
        pobyso_set_prec_sa_so(sollyaPrecSa)
1666
        sollyaPrecChanged = True
1667
    RRR = lowerBoundSa.parent()
1668
    intervalShrinkConstFactorSa = RRR('0.9')
1669
    absoluteErrorTypeSo = pobyso_absolute_so_so()
1670
    currentRangeSo = pobyso_bounds_to_range_sa_so(lowerBoundSa, upperBoundSa)
1671
    currentUpperBoundSa = upperBoundSa
1672
    currentLowerBoundSa = lowerBoundSa
1673
    #pobyso_autoprint(functionSo)
1674
    #pobyso_autoprint(degreeSo)
1675
    #pobyso_autoprint(currentRangeSo)
1676
    #pobyso_autoprint(absoluteErrorTypeSo)
1677
    ## What we want here is the polynomial without the variable change, 
1678
    #  since our actual variable will be x-intervalCenter defined over the 
1679
    #  domain [lowerBound-intervalCenter , upperBound-intervalCenter].
1680
    (polySo, intervalCenterSo, maxErrorSo) = \
1681
        pobyso_taylor_expansion_no_change_var_so_so(functionSo, degreeSo, 
1682
                                                    currentRangeSo, 
1683
                                                    absoluteErrorTypeSo)
1684
    maxErrorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo)
1685
    print "...after Taylor expansion."
1686
    ## Expression "approxAccurSa/2" (instead of just approxAccurSa) is added 
1687
    #  since we use polynomial shaving. In order for it to work, we must start 
1688
    #  with a better accuracy.
1689
    while maxErrorSa > approxAccurSa/2:
1690
        print "++Approximation error:", maxErrorSa.n()
1691
        sollya_lib_clear_obj(polySo)
1692
        sollya_lib_clear_obj(intervalCenterSo)
1693
        sollya_lib_clear_obj(maxErrorSo)
1694
        # Very empirical shrinking factor.
1695
        shrinkFactorSa = 1 /  (maxErrorSa/approxAccurSa).log2().abs()
1696
        print "Shrink factor:", \
1697
              shrinkFactorSa.n(), \
1698
              intervalShrinkConstFactorSa
1699
        print
1700
        #errorRatioSa = approxAccurSa/maxErrorSa
1701
        #print "Error ratio: ", errorRatioSa
1702
        # Make sure interval shrinks.
1703
        if shrinkFactorSa > intervalShrinkConstFactorSa:
1704
            actualShrinkFactorSa = intervalShrinkConstFactorSa
1705
            #print "Fixed"
1706
        else:
1707
            actualShrinkFactorSa = shrinkFactorSa
1708
            #print "Computed",shrinkFactorSa,maxErrorSa
1709
            #print shrinkFactorSa, maxErrorSa
1710
        #print "Shrink factor", actualShrinkFactorSa
1711
        currentUpperBoundSa = currentLowerBoundSa + \
1712
                                (currentUpperBoundSa - currentLowerBoundSa) * \
1713
                                actualShrinkFactorSa
1714
        #print "Current upper bound:", currentUpperBoundSa
1715
        sollya_lib_clear_obj(currentRangeSo)
1716
        # Check what is left with the bounds.
1717
        if currentUpperBoundSa <= currentLowerBoundSa or \
1718
              currentUpperBoundSa == currentLowerBoundSa.nextabove():
1719
            sollya_lib_clear_obj(absoluteErrorTypeSo)
1720
            print "Can't find an interval."
1721
            print "Use either or both a higher polynomial degree or a higher",
1722
            print "internal precision."
1723
            print "Aborting!"
1724
            if sollyaPrecChanged:
1725
                pobyso_set_prec_so_so(initialSollyaPrecSo)
1726
            sollya_lib_clear_obj(initialSollyaPrecSo)
1727
            return None
1728
        currentRangeSo = pobyso_bounds_to_range_sa_so(currentLowerBoundSa, 
1729
                                                      currentUpperBoundSa)
1730
        # print "New interval:",
1731
        # pobyso_autoprint(currentRangeSo)
1732
        #print "Second Taylor expansion call."
1733
        (polySo, intervalCenterSo, maxErrorSo) = \
1734
            pobyso_taylor_expansion_no_change_var_so_so(functionSo, degreeSo, 
1735
                                                        currentRangeSo, 
1736
                                                        absoluteErrorTypeSo)
1737
        #maxErrorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo, RRR)
1738
        #print "Max errorSo:",
1739
        #pobyso_autoprint(maxErrorSo)
1740
        maxErrorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo)
1741
        #print "Max errorSa:", maxErrorSa
1742
        #print "Sollya prec:",
1743
        #pobyso_autoprint(sollya_lib_get_prec(None))
1744
    # End while
1745
    sollya_lib_clear_obj(absoluteErrorTypeSo)
1746
    itpSo           = pobyso_constant_from_int_sa_so(floor(sollyaPrecSa/3))
1747
    ftpSo           = pobyso_constant_from_int_sa_so(floor(2*sollyaPrecSa/3))
1748
    maxPrecSo       = pobyso_constant_from_int_sa_so(sollyaPrecSa)
1749
    approxAccurSo   = pobyso_constant_sa_so(RR(approxAccurSa))
1750
    print "About to call polynomial rounding with:"
1751
    print "polySo: ",           ; pobyso_autoprint(polySo)
1752
    print "functionSo: ",       ; pobyso_autoprint(functionSo)
1753
    print "intervalCenterSo: ", ; pobyso_autoprint(intervalCenterSo)
1754
    print "currentRangeSo: ",   ; pobyso_autoprint(currentRangeSo)
1755
    print "itpSo: ",            ; pobyso_autoprint(itpSo)
1756
    print "ftpSo: ",            ; pobyso_autoprint(ftpSo)
1757
    print "maxPrecSo: ",        ; pobyso_autoprint(maxPrecSo)
1758
    print "approxAccurSo: ",    ; pobyso_autoprint(approxAccurSo)
1759
    (roundedPolySo, roundedPolyMaxErrSo) = \
1760
    pobyso_round_coefficients_progressive_so_so(polySo,
1761
                                                functionSo,
1762
                                                maxPrecSo,
1763
                                                currentRangeSo,
1764
                                                intervalCenterSo,
1765
                                                maxErrorSo,
1766
                                                approxAccurSo,
1767
                                                debug = True)
1768
    
1769
    sollya_lib_clear_obj(polySo)
1770
    sollya_lib_clear_obj(maxErrorSo)
1771
    sollya_lib_clear_obj(itpSo)
1772
    sollya_lib_clear_obj(ftpSo)
1773
    sollya_lib_clear_obj(approxAccurSo)
1774
    if sollyaPrecChanged:
1775
        pobyso_set_prec_so_so(initialSollyaPrecSo)
1776
    sollya_lib_clear_obj(initialSollyaPrecSo)
1777
    print "1: ", ; pobyso_autoprint(roundedPolySo)
1778
    print "2: ", ; pobyso_autoprint(currentRangeSo)
1779
    print "3: ", ; pobyso_autoprint(intervalCenterSo)
1780
    print "4: ", ; pobyso_autoprint(roundedPolyMaxErrSo)
1781
    return (roundedPolySo, currentRangeSo, intervalCenterSo, roundedPolyMaxErrSo)
1782
# End slz_compute_polynomial_and_interval_02
1783

    
1784
def slz_compute_reduced_polynomial(matrixRow,
1785
                                    knownMonomials,
1786
                                    var1,
1787
                                    var1Bound,
1788
                                    var2,
1789
                                    var2Bound):
1790
    """
1791
    Compute a polynomial from a single reduced matrix row.
1792
    This function was introduced in order to avoid the computation of the
1793
    all the polynomials  from the full matrix (even those built from rows 
1794
    that do no verify the Coppersmith condition) as this may involves 
1795
    expensive operations over (large) integers.
1796
    """
1797
    ## Check arguments.
1798
    if len(knownMonomials) == 0:
1799
        return None
1800
    # varNounds can be zero since 0^0 returns 1.
1801
    if (var1Bound < 0) or (var2Bound < 0):
1802
        return None
1803
    ## Initialisations. 
1804
    polynomialRing = knownMonomials[0].parent() 
1805
    currentPolynomial = polynomialRing(0)
1806
    # TODO: use zip instead of indices.
1807
    for colIndex in xrange(0, len(knownMonomials)):
1808
        currentCoefficient = matrixRow[colIndex]
1809
        if currentCoefficient != 0:
1810
            #print "Current coefficient:", currentCoefficient
1811
            currentMonomial = knownMonomials[colIndex]
1812
            #print "Monomial as multivariate polynomial:", \
1813
            #currentMonomial, type(currentMonomial)
1814
            degreeInVar1 = currentMonomial.degree(var1)
1815
            #print "Degree in var1", var1, ":", degreeInVar1
1816
            degreeInVar2 = currentMonomial.degree(var2)
1817
            #print "Degree in var2", var2, ":", degreeInVar2
1818
            if degreeInVar1 > 0:
1819
                currentCoefficient = \
1820
                    currentCoefficient / (var1Bound^degreeInVar1)
1821
                #print "varBound1 in degree:", var1Bound^degreeInVar1
1822
                #print "Current coefficient(1)", currentCoefficient
1823
            if degreeInVar2 > 0:
1824
                currentCoefficient = \
1825
                    currentCoefficient / (var2Bound^degreeInVar2)
1826
                #print "Current coefficient(2)", currentCoefficient
1827
            #print "Current reduced monomial:", (currentCoefficient * \
1828
            #                                    currentMonomial)
1829
            currentPolynomial += (currentCoefficient * currentMonomial)
1830
            #print "Current polynomial:", currentPolynomial
1831
        # End if
1832
    # End for colIndex.
1833
    #print "Type of the current polynomial:", type(currentPolynomial)
1834
    return(currentPolynomial)
1835
# End slz_compute_reduced_polynomial
1836
#
1837
def slz_compute_reduced_polynomials(reducedMatrix,
1838
                                        knownMonomials,
1839
                                        var1,
1840
                                        var1Bound,
1841
                                        var2,
1842
                                        var2Bound):
1843
    """
1844
    Legacy function, use slz_compute_reduced_polynomials_list
1845
    """
1846
    return(slz_compute_reduced_polynomials_list(reducedMatrix,
1847
                                                knownMonomials,
1848
                                                var1,
1849
                                                var1Bound,
1850
                                                var2,
1851
                                                var2Bound)
1852
    )
1853
#
1854
def slz_compute_reduced_polynomials_list(reducedMatrix,
1855
                                         knownMonomials,
1856
                                         var1,
1857
                                         var1Bound,
1858
                                         var2,
1859
                                         var2Bound):
1860
    """
1861
    From a reduced matrix, holding the coefficients, from a monomials list,
1862
    from the bounds of each variable, compute the corresponding polynomials
1863
    scaled back by dividing by the "right" powers of the variables bounds.
1864
    
1865
    The elements in knownMonomials must be of the "right" polynomial type.
1866
    They set the polynomial type of the output polynomials in list.
1867
    @param reducedMatrix:  the reduced matrix as output from LLL;
1868
    @param kwnonMonomials: the ordered list of the monomials used to
1869
                           build the polynomials;
1870
    @param var1:           the first variable (of the "right" type);
1871
    @param var1Bound:      the first variable bound;
1872
    @param var2:           the second variable (of the "right" type);
1873
    @param var2Bound:      the second variable bound.
1874
    @return: a list of polynomials obtained with the reduced coefficients
1875
             and scaled down with the bounds
1876
    """
1877
    
1878
    # TODO: check input arguments.
1879
    reducedPolynomials = []
1880
    #print "type var1:", type(var1), " - type var2:", type(var2)
1881
    for matrixRow in reducedMatrix.rows():
1882
        currentPolynomial = 0
1883
        for colIndex in xrange(0, len(knownMonomials)):
1884
            currentCoefficient = matrixRow[colIndex]
1885
            if currentCoefficient != 0:
1886
                #print "Current coefficient:", currentCoefficient
1887
                currentMonomial = knownMonomials[colIndex]
1888
                parentRing = currentMonomial.parent()
1889
                #print "Monomial as multivariate polynomial:", \
1890
                #currentMonomial, type(currentMonomial)
1891
                degreeInVar1 = currentMonomial.degree(parentRing(var1))
1892
                #print "Degree in var", var1, ":", degreeInVar1
1893
                degreeInVar2 = currentMonomial.degree(parentRing(var2))
1894
                #print "Degree in var", var2, ":", degreeInVar2
1895
                if degreeInVar1 > 0:
1896
                    currentCoefficient /= var1Bound^degreeInVar1
1897
                    #print "varBound1 in degree:", var1Bound^degreeInVar1
1898
                    #print "Current coefficient(1)", currentCoefficient
1899
                if degreeInVar2 > 0:
1900
                    currentCoefficient /= var2Bound^degreeInVar2
1901
                    #print "Current coefficient(2)", currentCoefficient
1902
                #print "Current reduced monomial:", (currentCoefficient * \
1903
                #                                    currentMonomial)
1904
                currentPolynomial += (currentCoefficient * currentMonomial)
1905
                #if degreeInVar1 == 0 and degreeInVar2 == 0:
1906
                    #print "!!!! constant term !!!!"
1907
                #print "Current polynomial:", currentPolynomial
1908
            # End if
1909
        # End for colIndex.
1910
        #print "Type of the current polynomial:", type(currentPolynomial)
1911
        reducedPolynomials.append(currentPolynomial)
1912
    return reducedPolynomials 
1913
# End slz_compute_reduced_polynomials_list.
1914

    
1915
def slz_compute_reduced_polynomials_list_from_rows(rowsList,
1916
                                                   knownMonomials,
1917
                                                   var1,
1918
                                                   var1Bound,
1919
                                                   var2,
1920
                                                   var2Bound):
1921
    """
1922
    From a list of rows, holding the coefficients, from a monomials list,
1923
    from the bounds of each variable, compute the corresponding polynomials
1924
    scaled back by dividing by the "right" powers of the variables bounds.
1925
    
1926
    The elements in knownMonomials must be of the "right" polynomial type.
1927
    They set the polynomial type of the output polynomials in list.
1928
    @param rowsList:       the reduced matrix as output from LLL;
1929
    @param kwnonMonomials: the ordered list of the monomials used to
1930
                           build the polynomials;
1931
    @param var1:           the first variable (of the "right" type);
1932
    @param var1Bound:      the first variable bound;
1933
    @param var2:           the second variable (of the "right" type);
1934
    @param var2Bound:      the second variable bound.
1935
    @return: a list of polynomials obtained with the reduced coefficients
1936
             and scaled down with the bounds
1937
    """
1938
    
1939
    # TODO: check input arguments.
1940
    reducedPolynomials = []
1941
    #print "type var1:", type(var1), " - type var2:", type(var2)
1942
    for matrixRow in rowsList:
1943
        currentPolynomial = 0
1944
        for colIndex in xrange(0, len(knownMonomials)):
1945
            currentCoefficient = matrixRow[colIndex]
1946
            if currentCoefficient != 0:
1947
                #print "Current coefficient:", currentCoefficient
1948
                currentMonomial = knownMonomials[colIndex]
1949
                parentRing = currentMonomial.parent()
1950
                #print "Monomial as multivariate polynomial:", \
1951
                #currentMonomial, type(currentMonomial)
1952
                degreeInVar1 = currentMonomial.degree(parentRing(var1))
1953
                #print "Degree in var", var1, ":", degreeInVar1
1954
                degreeInVar2 = currentMonomial.degree(parentRing(var2))
1955
                #print "Degree in var", var2, ":", degreeInVar2
1956
                if degreeInVar1 > 0:
1957
                    currentCoefficient /= var1Bound^degreeInVar1
1958
                    #print "varBound1 in degree:", var1Bound^degreeInVar1
1959
                    #print "Current coefficient(1)", currentCoefficient
1960
                if degreeInVar2 > 0:
1961
                    currentCoefficient /= var2Bound^degreeInVar2
1962
                    #print "Current coefficient(2)", currentCoefficient
1963
                #print "Current reduced monomial:", (currentCoefficient * \
1964
                #                                    currentMonomial)
1965
                currentPolynomial += (currentCoefficient * currentMonomial)
1966
                #if degreeInVar1 == 0 and degreeInVar2 == 0:
1967
                    #print "!!!! constant term !!!!"
1968
                #print "Current polynomial:", currentPolynomial
1969
            # End if
1970
        # End for colIndex.
1971
        #print "Type of the current polynomial:", type(currentPolynomial)
1972
        reducedPolynomials.append(currentPolynomial)
1973
    return reducedPolynomials 
1974
# End slz_compute_reduced_polynomials_list_from_rows.
1975
#
1976
def slz_compute_scaled_function(functionSa,
1977
                                lowerBoundSa,
1978
                                upperBoundSa,
1979
                                floatingPointPrecSa,
1980
                                debug=False):
1981
    """
1982
    From a function, compute the scaled function whose domain
1983
    is included in [1, 2) and whose image is also included in [1,2).
1984
    Return a tuple: 
1985
    [0]: the scaled function
1986
    [1]: the scaled domain lower bound
1987
    [2]: the scaled domain upper bound
1988
    [3]: the scaled image lower bound
1989
    [4]: the scaled image upper bound
1990
    """
1991
    ## The variable can be called anything.
1992
    x = functionSa.variables()[0]
1993
    # Scalling the domain -> [1,2[.
1994
    boundsIntervalRifSa = RealIntervalField(floatingPointPrecSa)
1995
    domainBoundsIntervalSa = boundsIntervalRifSa(lowerBoundSa, upperBoundSa)
1996
    (invDomainScalingExpressionSa, domainScalingExpressionSa) = \
1997
        slz_interval_scaling_expression(domainBoundsIntervalSa, x)
1998
    if debug:
1999
        print "domainScalingExpression for argument :", \
2000
        invDomainScalingExpressionSa
2001
        print "function: ", functionSa
2002
    ff = functionSa.subs({x : domainScalingExpressionSa})
2003
    ## Bless expression as a function.
2004
    ff = ff.function(x)
2005
    #ff = f.subs_expr(x==domainScalingExpressionSa)
2006
    #domainScalingFunction(x) = invDomainScalingExpressionSa
2007
    domainScalingFunction = invDomainScalingExpressionSa.function(x)
2008
    scaledLowerBoundSa = \
2009
        domainScalingFunction(lowerBoundSa).n(prec=floatingPointPrecSa)
2010
    scaledUpperBoundSa = \
2011
        domainScalingFunction(upperBoundSa).n(prec=floatingPointPrecSa)
2012
    if debug:
2013
        print 'ff:', ff, "- Domain:", scaledLowerBoundSa, \
2014
              scaledUpperBoundSa
2015
    ## If both bounds are negative, swap them.
2016
    if lowerBoundSa < 0 and upperBoundSa < 0:
2017
        scaledLowerBoundSa, scaledUpperBoundSa = \
2018
        scaledUpperBoundSa,scaledLowerBoundSa
2019
    #
2020
    # Scalling the image -> [1,2[.
2021
    flbSa = ff(scaledLowerBoundSa).n(prec=floatingPointPrecSa)
2022
    fubSa = ff(scaledUpperBoundSa).n(prec=floatingPointPrecSa)
2023
    if flbSa <= fubSa: # Increasing
2024
        imageBinadeBottomSa = floor(flbSa.log2())
2025
    else: # Decreasing
2026
        imageBinadeBottomSa = floor(fubSa.log2())
2027
    if debug:
2028
        print 'ff:', ff, '- Image:', flbSa, fubSa, imageBinadeBottomSa
2029
    imageBoundsIntervalSa = boundsIntervalRifSa(flbSa, fubSa)
2030
    (invImageScalingExpressionSa,imageScalingExpressionSa) = \
2031
        slz_interval_scaling_expression(imageBoundsIntervalSa, x)
2032
    if debug:
2033
        print "imageScalingExpression for argument :", \
2034
              invImageScalingExpressionSa
2035
    iis = invImageScalingExpressionSa.function(x)
2036
    fff = iis.subs({x:ff})
2037
    if debug:
2038
        print "fff:", fff,
2039
        print " - Image:", fff(scaledLowerBoundSa), fff(scaledUpperBoundSa)
2040
    return([fff, scaledLowerBoundSa, scaledUpperBoundSa, \
2041
            fff(scaledLowerBoundSa), fff(scaledUpperBoundSa)])
2042
# End slz_compute_scaled_function
2043

    
2044
def slz_fix_bounds_for_binades(lowerBound, 
2045
                               upperBound, 
2046
                               func = None, 
2047
                               domainDirection = -1,
2048
                               imageDirection  = -1):
2049
    """
2050
    Assuming the function is increasing or decreasing over the 
2051
    [lowerBound, upperBound] interval, return a lower bound lb and 
2052
    an upper bound ub such that:
2053
    - lb and ub belong to the same binade;
2054
    - func(lb) and func(ub) belong to the same binade.
2055
    domainDirection indicate how bounds move to fit into the same binade:
2056
    - a negative value move the upper bound down;
2057
    - a positive value move the lower bound up.
2058
    imageDirection indicate how bounds move in order to have their image
2059
    fit into the same binade, variation of the function is also condidered.
2060
    For an increasing function:
2061
    - negative value moves the upper bound down (and its image value as well);
2062
    - a positive values moves the lower bound up (and its image value as well);
2063
    For a decreasing function it is the other way round.
2064
    """
2065
    ## Arguments check
2066
    if lowerBound > upperBound:
2067
        return None
2068
    if lowerBound == upperBound:
2069
        return (lowerBound, upperBound)
2070
    if func is None:
2071
        return None
2072
    #
2073
    #varFunc = func.variables()[0]
2074
    lb = lowerBound
2075
    ub = upperBound
2076
    lbBinade = slz_compute_binade(lb) 
2077
    ubBinade = slz_compute_binade(ub)
2078
    ## Domain binade.
2079
    while lbBinade != ubBinade:
2080
        newIntervalWidth = (ub - lb) / 2
2081
        if domainDirection < 0:
2082
            ub = lb + newIntervalWidth
2083
            ubBinade = slz_compute_binade(ub)
2084
        else:
2085
            lb = lb + newIntervalWidth
2086
            lbBinade = slz_compute_binade(lb) 
2087
    #print "sfbfb: lower bound:", lb.str(truncate=False)
2088
    #print "sfbfb: upper bound:", ub.str(truncate=False)
2089
    ## At this point, both bounds belond to the same binade.
2090
    ## Image binade.
2091
    if lb == ub:
2092
        return (lb, ub)
2093
    lbImg = func(lb)
2094
    ubImg = func(ub)
2095
    funcIsInc = ((ubImg - lbImg) >= 0)
2096
    lbImgBinade = slz_compute_binade(lbImg)
2097
    ubImgBinade = slz_compute_binade(ubImg)
2098
    while lbImgBinade != ubImgBinade:
2099
        newIntervalWidth = (ub - lb) / 2
2100
        if imageDirection < 0:
2101
            if funcIsInc:
2102
                ub = lb + newIntervalWidth
2103
                ubImgBinade = slz_compute_binade(func(ub))
2104
                #print ubImgBinade
2105
            else:
2106
                lb = lb + newIntervalWidth
2107
                lbImgBinade = slz_compute_binade(func(lb))
2108
                #print lbImgBinade
2109
        else:
2110
            if funcIsInc:
2111
                lb = lb + newIntervalWidth
2112
                lbImgBinade = slz_compute_binade(func(lb))
2113
                #print lbImgBinade
2114
            else:
2115
                ub = lb + newIntervalWidth
2116
                ubImgBinade = slz_compute_binade(func(ub))
2117
                #print ubImgBinade
2118
    # End while lbImgBinade != ubImgBinade:
2119
    return (lb, ub)
2120
# End slz_fix_bounds_for_binades.
2121

    
2122
def slz_float_poly_of_float_to_rat_poly_of_rat(polyOfFloat):
2123
    # Create a polynomial over the rationals.
2124
    ratPolynomialRing = QQ[str(polyOfFloat.variables()[0])]
2125
    return(ratPolynomialRing(polyOfFloat))
2126
# End slz_float_poly_of_float_to_rat_poly_of_rat.
2127

    
2128
def slz_float_poly_of_float_to_rat_poly_of_rat_pow_two(polyOfFloat):
2129
    """
2130
     Create a polynomial over the rationals where all denominators are
2131
     powers of two.
2132
    """
2133
    polyVariable = polyOfFloat.variables()[0] 
2134
    RPR = QQ[str(polyVariable)]
2135
    polyCoeffs      = polyOfFloat.coefficients()
2136
    #print polyCoeffs
2137
    polyExponents   = polyOfFloat.exponents()
2138
    #print polyExponents
2139
    polyDenomPtwoCoeffs = []
2140
    for coeff in polyCoeffs:
2141
        polyDenomPtwoCoeffs.append(sno_float_to_rat_pow_of_two_denom(coeff))
2142
        #print "Converted coefficient:", sno_float_to_rat_pow_of_two_denom(coeff),
2143
        #print type(sno_float_to_rat_pow_of_two_denom(coeff))
2144
    ratPoly = RPR(0)
2145
    #print type(ratPoly)
2146
    ## !!! CAUTION !!! Do not use the RPR(coeff * polyVariagle^exponent)
2147
    #  The coefficient becomes plainly wrong when exponent == 0.
2148
    #  No clue as to why.
2149
    for coeff, exponent in zip(polyDenomPtwoCoeffs, polyExponents):
2150
        ratPoly += coeff * RPR(polyVariable^exponent)
2151
    return ratPoly
2152
# End slz_float_poly_of_float_to_rat_poly_of_rat.
2153

    
2154
def slz_get_intervals_and_polynomials(functionSa, degreeSa, 
2155
                                      lowerBoundSa, 
2156
                                      upperBoundSa, 
2157
                                      floatingPointPrecSa, 
2158
                                      internalSollyaPrecSa, 
2159
                                      approxAccurSa):
2160
    """
2161
    Under the assumption that:
2162
    - functionSa is monotonic on the [lowerBoundSa, upperBoundSa] interval;
2163
    - lowerBound and upperBound belong to the same binade.
2164
    from a:
2165
    - function;
2166
    - a degree
2167
    - a pair of bounds;
2168
    - the floating-point precision we work on;
2169
    - the internal Sollya precision;
2170
    - the requested approximation error
2171
    The initial interval is, possibly, splitted into smaller intervals.
2172
    It return a list of tuples, each made of:
2173
    - a first polynomial (without the changed variable f(x) = p(x-x0));
2174
    - a second polynomial (with a changed variable f(x) = q(x))
2175
    - the approximation interval;
2176
    - the center, x0, of the interval;
2177
    - the corresponding approximation error.
2178
    TODO: fix endless looping for some parameters sets.
2179
    """
2180
    resultArray = []
2181
    # Set Sollya to the necessary internal precision.
2182
    sollyaPrecChangedSa = False
2183
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
2184
    if internalSollyaPrecSa > initialSollyaPrecSa:
2185
        if internalSollyaPrecSa <= 2:
2186
            print inspect.stack()[0][3], ": precision change <=2 requested."
2187
        pobyso_set_prec_sa_so(internalSollyaPrecSa)
2188
        sollyaPrecChangedSa = True
2189
    #
2190
    x = functionSa.variables()[0] # Actual variable name can be anything.
2191
    # Scaled function: [1=,2] -> [1,2].
2192
    (fff, scaledLowerBoundSa, scaledUpperBoundSa,                             \
2193
            scaledLowerBoundImageSa, scaledUpperBoundImageSa) =               \
2194
                slz_compute_scaled_function(functionSa,                       \
2195
                                            lowerBoundSa,                     \
2196
                                            upperBoundSa,                     \
2197
                                            floatingPointPrecSa)
2198
    # In case bounds were in the negative real one may need to
2199
    # switch scaled bounds.
2200
    if scaledLowerBoundSa > scaledUpperBoundSa:
2201
        scaledLowerBoundSa, scaledUpperBoundSa = \
2202
            scaledUpperBoundSa, scaledLowerBoundSa
2203
        #print "Switching!"
2204
    print "Approximation accuracy: ", RR(approxAccurSa)
2205
    # Prepare the arguments for the Taylor expansion computation with Sollya.
2206
    functionSo = \
2207
      pobyso_parse_string_sa_so(fff._assume_str().replace('_SAGE_VAR_', ''))
2208
    degreeSo = pobyso_constant_from_int_sa_so(degreeSa)
2209
    scaledBoundsSo = pobyso_bounds_to_range_sa_so(scaledLowerBoundSa, 
2210
                                                  scaledUpperBoundSa)
2211

    
2212
    realIntervalField = RealIntervalField(max(lowerBoundSa.parent().precision(),
2213
                                              upperBoundSa.parent().precision()))
2214
    currentScaledLowerBoundSa = scaledLowerBoundSa
2215
    currentScaledUpperBoundSa = scaledUpperBoundSa
2216
    while True:
2217
        ## Compute the first Taylor expansion.
2218
        print "Computing a Taylor expansion..."
2219
        (polySo, boundsSo, intervalCenterSo, maxErrorSo) = \
2220
            slz_compute_polynomial_and_interval(functionSo, degreeSo,
2221
                                        currentScaledLowerBoundSa, 
2222
                                        currentScaledUpperBoundSa,
2223
                                        approxAccurSa, internalSollyaPrecSa)
2224
        print "...done."
2225
        ## If slz_compute_polynomial_and_interval fails, it returns None.
2226
        #  This value goes to the first variable: polySo. Other variables are
2227
        #  not assigned and should not be tested.
2228
        if polySo is None:
2229
            print "slz_get_intervals_and_polynomials: Aborting and returning None!"
2230
            if precChangedSa:
2231
                pobyso_set_prec_so_so(initialSollyaPrecSo)
2232
            sollya_lib_clear_obj(initialSollyaPrecSo)
2233
            sollya_lib_clear_obj(functionSo)
2234
            sollya_lib_clear_obj(degreeSo)
2235
            sollya_lib_clear_obj(scaledBoundsSo)
2236
            return None
2237
        ## Add to the result array.
2238
        ### Change variable stuff in Sollya x -> x0-x.
2239
        changeVarExpressionSo = \
2240
            sollya_lib_build_function_sub( \
2241
                              sollya_lib_build_function_free_variable(), 
2242
                              sollya_lib_copy_obj(intervalCenterSo))
2243
        polyVarChangedSo = \
2244
                    sollya_lib_evaluate(polySo, changeVarExpressionSo)
2245
        #### Get rid of the variable change Sollya stuff. 
2246
        sollya_lib_clear_obj(changeVarExpressionSo)
2247
        resultArray.append((polySo, polyVarChangedSo, boundsSo,
2248
                            intervalCenterSo, maxErrorSo))
2249
        boundsSa = pobyso_range_to_interval_so_sa(boundsSo, realIntervalField)
2250
        errorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo)
2251
        print "Computed approximation error:", errorSa.n(digits=10)
2252
        # If the error and interval are OK a the first try, just return.
2253
        if (boundsSa.endpoints()[1] >= scaledUpperBoundSa) and \
2254
            (errorSa <= approxAccurSa):
2255
            if sollyaPrecChangedSa:
2256
                pobyso_set_prec_so_so(initialSollyaPrecSo)
2257
            sollya_lib_clear_obj(initialSollyaPrecSo)
2258
            sollya_lib_clear_obj(functionSo)
2259
            sollya_lib_clear_obj(degreeSo)
2260
            sollya_lib_clear_obj(scaledBoundsSo)
2261
            #print "Approximation error:", errorSa
2262
            return resultArray
2263
        ## The returned interval upper bound does not reach the requested upper
2264
        #  upper bound: compute the next upper bound.
2265
        ## The following ratio is always >= 1. If errorSa, we may want to
2266
        #  enlarge the interval
2267
        currentErrorRatio = approxAccurSa / errorSa
2268
        ## --|--------------------------------------------------------------|--
2269
        #  --|--------------------|--------------------------------------------
2270
        #  --|----------------------------|------------------------------------
2271
        ## Starting point for the next upper bound.
2272
        boundsWidthSa = boundsSa.endpoints()[1] - boundsSa.endpoints()[0]
2273
        # Compute the increment. 
2274
        newBoundsWidthSa = \
2275
                        ((currentErrorRatio.log() / 10) + 1) * boundsWidthSa
2276
        currentScaledLowerBoundSa = boundsSa.endpoints()[1]
2277
        currentScaledUpperBoundSa = boundsSa.endpoints()[1] + newBoundsWidthSa
2278
        # Take into account the original interval upper bound.                     
2279
        if currentScaledUpperBoundSa > scaledUpperBoundSa:
2280
            currentScaledUpperBoundSa = scaledUpperBoundSa
2281
        if currentScaledUpperBoundSa == currentScaledLowerBoundSa:
2282
            print "Can't shrink the interval anymore!"
2283
            print "You should consider increasing the Sollya internal precision"
2284
            print "or the polynomial degree."
2285
            print "Giving up!"
2286
            if precChangedSa:
2287
                pobyso_set_prec_so_so(initialSollyaPrecSo)
2288
            sollya_lib_clear_obj(initialSollyaPrecSo)
2289
            sollya_lib_clear_obj(functionSo)
2290
            sollya_lib_clear_obj(degreeSo)
2291
            sollya_lib_clear_obj(scaledBoundsSo)
2292
            return None
2293
        # Compute the other expansions.
2294
        # Test for insufficient precision.
2295
# End slz_get_intervals_and_polynomials
2296

    
2297
def slz_interval_scaling_expression(boundsInterval, expVar):
2298
    """
2299
    Compute the scaling expression to map an interval that spans at most
2300
    a single binade into [1, 2) and the inverse expression as well.
2301
    If the interval spans more than one binade, result is wrong since
2302
    scaling is only based on the lower bound.
2303
    Not very sure that the transformation makes sense for negative numbers.
2304
    """
2305
    # The "one of the bounds is 0" special case. Here we consider
2306
    # the interval as the subnormals binade.
2307
    if boundsInterval.endpoints()[0] == 0 or boundsInterval.endpoints()[1] == 0:
2308
        # The upper bound is (or should be) positive.
2309
        if boundsInterval.endpoints()[0] == 0:
2310
            if boundsInterval.endpoints()[1] == 0:
2311
                return None
2312
            binade = slz_compute_binade(boundsInterval.endpoints()[1])
2313
            l2     = boundsInterval.endpoints()[1].abs().log2()
2314
            # The upper bound is a power of two
2315
            if binade == l2:
2316
                scalingCoeff    = 2^(-binade)
2317
                invScalingCoeff = 2^(binade)
2318
                scalingOffset   = 1
2319
                return \
2320
                  ((scalingCoeff * expVar + scalingOffset).function(expVar),
2321
                  ((expVar - scalingOffset) * invScalingCoeff).function(expVar))
2322
            else:
2323
                scalingCoeff    = 2^(-binade-1)
2324
                invScalingCoeff = 2^(binade+1)
2325
                scalingOffset   = 1
2326
                return((scalingCoeff * expVar + scalingOffset),\
2327
                    ((expVar - scalingOffset) * invScalingCoeff))
2328
        # The lower bound is (or should be) negative.
2329
        if boundsInterval.endpoints()[1] == 0:
2330
            if boundsInterval.endpoints()[0] == 0:
2331
                return None
2332
            binade = slz_compute_binade(boundsInterval.endpoints()[0])
2333
            l2     = boundsInterval.endpoints()[0].abs().log2()
2334
            # The upper bound is a power of two
2335
            if binade == l2:
2336
                scalingCoeff    = -2^(-binade)
2337
                invScalingCoeff = -2^(binade)
2338
                scalingOffset   = 1
2339
                return((scalingCoeff * expVar + scalingOffset),\
2340
                    ((expVar - scalingOffset) * invScalingCoeff))
2341
            else:
2342
                scalingCoeff    = -2^(-binade-1)
2343
                invScalingCoeff = -2^(binade+1)
2344
                scalingOffset   = 1
2345
                return((scalingCoeff * expVar + scalingOffset),\
2346
                   ((expVar - scalingOffset) * invScalingCoeff))
2347
    # General case
2348
    lbBinade = slz_compute_binade(boundsInterval.endpoints()[0])
2349
    ubBinade = slz_compute_binade(boundsInterval.endpoints()[1])
2350
    # We allow for a single exception in binade spanning is when the
2351
    # "outward bound" is a power of 2 and its binade is that of the
2352
    # "inner bound" + 1.
2353
    if boundsInterval.endpoints()[0] > 0:
2354
        ubL2 = boundsInterval.endpoints()[1].abs().log2()
2355
        if lbBinade != ubBinade:
2356
            print "Different binades."
2357
            if ubL2 != ubBinade:
2358
                print "Not a power of 2."
2359
                return None
2360
            elif abs(ubBinade - lbBinade) > 1:
2361
                print "Too large span:", abs(ubBinade - lbBinade)
2362
                return None
2363
    else:
2364
        lbL2 = boundsInterval.endpoints()[0].abs().log2()
2365
        if lbBinade != ubBinade:
2366
            print "Different binades."
2367
            if lbL2 != lbBinade:
2368
                print "Not a power of 2."
2369
                return None
2370
            elif abs(ubBinade - lbBinade) > 1:
2371
                print "Too large span:", abs(ubBinade - lbBinade)
2372
                return None
2373
    #print "Lower bound binade:", binade
2374
    if boundsInterval.endpoints()[0] > 0:
2375
        return \
2376
            ((2^(-lbBinade) * expVar).function(expVar),
2377
             (2^(lbBinade) * expVar).function(expVar))
2378
    else:
2379
        return \
2380
            ((-(2^(-ubBinade)) * expVar).function(expVar),
2381
             (-(2^(ubBinade)) * expVar).function(expVar))
2382
"""
2383
    # Code sent to attic. Should be the base for a 
2384
    # "slz_interval_translate_expression" rather than scale.
2385
    # Extra control and special cases code  added in  
2386
    # slz_interval_scaling_expression could (should ?) be added to
2387
    # this new function.
2388
    # The scaling offset is only used for negative numbers.
2389
    # When the absolute value of the lower bound is < 0.
2390
    if abs(boundsInterval.endpoints()[0]) < 1:
2391
        if boundsInterval.endpoints()[0] >= 0:
2392
            scalingCoeff = 2^floor(boundsInterval.endpoints()[0].log2())
2393
            invScalingCoeff = 1/scalingCoeff
2394
            return((scalingCoeff * expVar, 
2395
                    invScalingCoeff * expVar))
2396
        else:
2397
            scalingCoeff = \
2398
                2^(floor((-boundsInterval.endpoints()[0]).log2()) - 1)
2399
            scalingOffset = -3 * scalingCoeff
2400
            return((scalingCoeff * expVar + scalingOffset,
2401
                    1/scalingCoeff * expVar + 3))
2402
    else: 
2403
        if boundsInterval.endpoints()[0] >= 0:
2404
            scalingCoeff = 2^floor(boundsInterval.endpoints()[0].log2())
2405
            scalingOffset = 0
2406
            return((scalingCoeff * expVar, 
2407
                    1/scalingCoeff * expVar))
2408
        else:
2409
            scalingCoeff = \
2410
                2^(floor((-boundsInterval.endpoints()[1]).log2()))
2411
            scalingOffset =  -3 * scalingCoeff
2412
            #scalingOffset = 0
2413
            return((scalingCoeff * expVar + scalingOffset,
2414
                    1/scalingCoeff * expVar + 3))
2415
"""
2416
# End slz_interval_scaling_expression
2417
   
2418
def slz_interval_and_polynomial_to_sage(polyRangeCenterErrorSo):
2419
    """
2420
    Compute the Sage version of the Taylor polynomial and it's
2421
    companion data (interval, center...)
2422
    The input parameter is a five elements tuple:
2423
    - [0]: the polyomial (without variable change), as polynomial over a
2424
           real ring;
2425
    - [1]: the polyomial (with variable change done in Sollya), as polynomial
2426
           over a real ring;
2427
    - [2]: the interval (as Sollya range);
2428
    - [3]: the interval center;
2429
    - [4]: the approximation error. 
2430
    
2431
    The function returns a 5 elements tuple: formed with all the 
2432
    input elements converted into their Sollya counterpart.
2433
    """
2434
    #print "Sollya polynomial to convert:", 
2435
    #pobyso_autoprint(polyRangeCenterErrorSo)
2436
    polynomialSa = pobyso_float_poly_so_sa(polyRangeCenterErrorSo[0])
2437
    #print "Polynomial after first conversion: ", pobyso_autoprint(polyRangeCenterErrorSo[1])
2438
    polynomialChangedVarSa = pobyso_float_poly_so_sa(polyRangeCenterErrorSo[1])
2439
    intervalSa = \
2440
            pobyso_get_interval_from_range_so_sa(polyRangeCenterErrorSo[2])
2441
    centerSa = \
2442
            pobyso_get_constant_as_rn_with_rf_so_sa(polyRangeCenterErrorSo[3])
2443
    errorSa = \
2444
            pobyso_get_constant_as_rn_with_rf_so_sa(polyRangeCenterErrorSo[4])
2445
    return((polynomialSa, polynomialChangedVarSa, intervalSa, centerSa, errorSa))
2446
    # End slz_interval_and_polynomial_to_sage
2447

    
2448
def slz_is_htrn(argument, 
2449
                function, 
2450
                targetAccuracy, 
2451
                targetRF = None, 
2452
                targetPlusOnePrecRF = None, 
2453
                quasiExactRF = None):
2454
    """
2455
      Check if an element (argument) of the domain of function (function)
2456
      yields a HTRN case (rounding to next) for the target precision 
2457
      (as impersonated by targetRF) for a given accuracy (targetAccuracy). 
2458
      
2459
      The strategy is: 
2460
      - compute the image at high (quasi-exact) precision;
2461
      - round it to nearest to precision;
2462
      - round it to exactly to (precision+1), the computed number has two
2463
        midpoint neighbors;
2464
      - check the distance between these neighbors and the quasi-exact value;
2465
        - if none of them is closer than the targetAccuracy, return False,
2466
        - else return true.
2467
      - Powers of two are a special case when comparing the midpoint in
2468
        the 0 direction..
2469
    """
2470
    ## Arguments filtering. 
2471
    ## TODO: filter the first argument for consistence.
2472
    if targetRF is None:
2473
        targetRF = argument.parent()
2474
    ## Ditto for the real field holding the midpoints.
2475
    if targetPlusOnePrecRF is None:
2476
        targetPlusOnePrecRF = RealField(targetRF.prec()+1)
2477
    ## If no quasiExactField is provided, create a high accuracy RealField.
2478
    if quasiExactRF is None:
2479
        quasiExactRF = RealField(1536)
2480
    function                      = function.function(function.variables()[0])
2481
    exactArgument                 = quasiExactRF(argument)
2482
    quasiExactValue               = function(exactArgument)
2483
    roundedValue                  = targetRF(quasiExactValue)
2484
    roundedValuePrecPlusOne       = targetPlusOnePrecRF(roundedValue)
2485
    # Upper midpoint.
2486
    roundedValuePrecPlusOneNext   = roundedValuePrecPlusOne.nextabove()
2487
    # Lower midpoint.
2488
    roundedValuePrecPlusOnePrev   = roundedValuePrecPlusOne.nextbelow()
2489
    binade                        = slz_compute_binade(roundedValue)
2490
    binadeCorrectedTargetAccuracy = targetAccuracy * 2^binade
2491
    #print "Argument:", argument
2492
    #print "f(x):", roundedValue, binade, floor(binade), ceil(binade)
2493
    #print "Binade:", binade
2494
    #print "binadeCorrectedTargetAccuracy:", \
2495
    #binadeCorrectedTargetAccuracy.n(prec=100)
2496
    #print "binadeCorrectedTargetAccuracy:", \
2497
    #  binadeCorrectedTargetAccuracy.n(prec=100).str(base=2)
2498
    #print "Upper midpoint:", \
2499
    #  roundedValuePrecPlusOneNext.n(prec=targetPlusOnePrecRF.prec()).str(base=2)
2500
    #print "Rounded value :", \
2501
    #  roundedValuePrecPlusOne.n(prec=targetPlusOnePrecRF.prec()).str(base=2), \
2502
    #  roundedValuePrecPlusOne.ulp().n(prec=2).str(base=2)
2503
    #print "QuasiEx value :", quasiExactValue.n(prec=250).str(base=2)
2504
    #print "Lower midpoint:", \
2505
    #  roundedValuePrecPlusOnePrev.n(prec=targetPlusOnePrecRF.prec()).str(base=2)
2506
    ## Make quasiExactValue = 0 a special case to move it out of the way.
2507
    if quasiExactValue == 0:
2508
        return False
2509
    ## Begining of the general case : check with the midpoint of 
2510
    #  greatest absolute value.
2511
    if quasiExactValue > 0:
2512
        if abs(quasiExactRF(roundedValuePrecPlusOneNext) - quasiExactValue) <\
2513
           binadeCorrectedTargetAccuracy:
2514
            #print "Too close to the upper midpoint: ", \
2515
            #abs(quasiExactRF(roundedValuePrecPlusOneNext) - quasiExactValue).n(prec=100)
2516
            #print argument.n().str(base=16, \
2517
            #  no_sci=False)
2518
            #print "Lower midpoint:", \
2519
            #  roundedValuePrecPlusOnePrev.n(prec=targetPlusOnePrecRF.prec()).str(base=2)
2520
            #print "Value         :", \
2521
            # quasiExactValue.n(prec=200).str(base=2)
2522
            #print "Upper midpoint:", \
2523
            #  roundedValuePrecPlusOneNext.n(prec=targetPlusOnePrecRF.prec()).str(base=2)
2524
            return True
2525
    else: # quasiExactValue < 0, the 0 case has been previously filtered out.
2526
        if abs(quasiExactRF(roundedValuePrecPlusOnePrev) - quasiExactValue) < \
2527
           binadeCorrectedTargetAccuracy:
2528
            #print "Too close to the upper midpoint: ", \
2529
            #  abs(quasiExactRF(roundedValuePrecPlusOneNext) - quasiExactValue).n(prec=100)
2530
            #print argument.n().str(base=16, \
2531
            #  no_sci=False)
2532
            #print "Lower midpoint:", \
2533
            #  roundedValuePrecPlusOnePrev.n(prec=targetPlusOnePrecRF.prec()).str(base=2)
2534
            #print "Value         :", \
2535
            #  quasiExactValue.n(prec=200).str(base=2)
2536
            #print "Upper midpoint:", \
2537
            #  roundedValuePrecPlusOneNext.n(prec=targetPlusOnePrecRF.prec()).str(base=2)
2538
            #print
2539
            return True
2540
#2345678901234567890123456789012345678901234567890123456789012345678901234567890
2541
    ## For the midpoint of smaller absolute value, 
2542
    #  split cases with the powers of 2.
2543
    if sno_abs_is_power_of_two(roundedValue):
2544
        if quasiExactValue > 0:        
2545
            if abs(quasiExactRF(roundedValuePrecPlusOnePrev) - quasiExactValue) <\
2546
               binadeCorrectedTargetAccuracy / 2:
2547
                #print "Lower midpoint:", \
2548
                #  roundedValuePrecPlusOnePrev.n(prec=targetPlusOnePrecRF.prec()).str(base=2)
2549
                #print "Value         :", \
2550
                #  quasiExactValue.n(prec=200).str(base=2)
2551
                #print "Upper midpoint:", \
2552
                #  roundedValuePrecPlusOneNext.n(prec=targetPlusOnePrecRF.prec()).str(base=2)
2553
                #print
2554
                return True
2555
        else:
2556
            if abs(quasiExactRF(roundedValuePrecPlusOneNext) - quasiExactValue) < \
2557
              binadeCorrectedTargetAccuracy / 2:
2558
                #print "Lower midpoint:", \
2559
                #  roundedValuePrecPlusOnePrev.n(prec=targetPlusOnePrecRF.prec()).str(base=2)
2560
                #print "Value         :", 
2561
                #  quasiExactValue.n(prec=200).str(base=2)
2562
                #print "Upper midpoint:", 
2563
                #  roundedValuePrecPlusOneNext.n(prec=targetPlusOnePrecRF.prec()).str(base=2)
2564
                #print
2565
                return True
2566
#2345678901234567890123456789012345678901234567890123456789012345678901234567890
2567
    else: ## Not a power of 2, regular comparison with the midpoint of 
2568
          #  smaller absolute value.
2569
        if quasiExactValue > 0:        
2570
            if abs(quasiExactRF(roundedValuePrecPlusOnePrev) - quasiExactValue) < \
2571
              binadeCorrectedTargetAccuracy:
2572
                #print "Lower midpoint:", \
2573
                #  roundedValuePrecPlusOnePrev.n(prec=targetPlusOnePrecRF.prec()).str(base=2)
2574
                #print "Value         :", \
2575
                #  quasiExactValue.n(prec=200).str(base=2)
2576
                #print "Upper midpoint:", \
2577
                #  roundedValuePrecPlusOneNext.n(prec=targetPlusOnePrecRF.prec()).str(base=2)
2578
                #print
2579
                return True
2580
        else: # quasiExactValue <= 0
2581
            if abs(quasiExactRF(roundedValuePrecPlusOneNext) - quasiExactValue) < \
2582
              binadeCorrectedTargetAccuracy:
2583
                #print "Lower midpoint:", \
2584
                #  roundedValuePrecPlusOnePrev.n(prec=targetPlusOnePrecRF.prec()).str(base=2)
2585
                #print "Value         :", \
2586
                #  quasiExactValue.n(prec=200).str(base=2)
2587
                #print "Upper midpoint:", \
2588
                #  roundedValuePrecPlusOneNext.n(prec=targetPlusOnePrecRF.prec()).str(base=2)
2589
                #print
2590
                return True
2591
    #print "distance to the upper midpoint:", \
2592
    #  abs(quasiExactRF(roundedValuePrecPlusOneNext) - quasiExactValue).n(prec=100).str(base=2)
2593
    #print "distance to the lower midpoint:", \
2594
    #  abs(quasiExactRF(roundedValuePrecPlusOnePrev) - quasiExactValue).n(prec=100).str(base=2)
2595
    return False
2596
# End slz_is_htrn
2597
#
2598
def slz_pm1():
2599
    """
2600
    Compute a uniform RV in {-1, 1} (not zero).
2601
    """
2602
    ## function getrandbits(N) generates a long int with N random bits.
2603
    #  Here it generates either 0 or 1. The multiplication by 2 and the
2604
    #  subtraction of 1 make that:
2605
    #  if getranbits(1) == 0 -> O * 2 - 1 = -1
2606
    #  else                  -> 1 * 1 - 1 =  1. 
2607
    return getrandbits(1) * 2-1
2608
# End slz_pm1.
2609
#
2610
def slz_random_proj_pm1(r, c):
2611
    """
2612
    r x c matrix with \pm 1 ({-1, 1}) coefficients
2613
    """
2614
    M = matrix(r, c)
2615
    for i in range(0, r):
2616
        for j in range (0, c):
2617
            M[i,j] = slz_pm1()
2618
    return M
2619
# End random_proj_pm1.
2620
#
2621
def slz_random_proj_uniform(n, r, c):
2622
    """
2623
    r x c integer matrix with uniform n-bit integer elements.
2624
    """
2625
    M = matrix(r, c)
2626
    for i in range(0, r):
2627
        for j in range (0, c):
2628
            M[i,j] = slz_uniform(n)
2629
    return M       
2630
# End slz_random_proj_uniform.
2631
#
2632
def slz_rat_poly_of_int_to_poly_for_coppersmith(ratPolyOfInt, 
2633
                                                precision,
2634
                                                targetHardnessToRound,
2635
                                                variable1,
2636
                                                variable2):
2637
    """
2638
    Creates a new multivariate polynomial with integer coefficients for use
2639
     with the Coppersmith method.
2640
    A the same time it computes :
2641
    - 2^K (N);
2642
    - 2^k (bound on the second variable)
2643
    - lcm
2644

    
2645
    :param ratPolyOfInt: a polynomial with rational coefficients and integer
2646
                         variables.
2647
    :param precision: the precision of the floating-point coefficients.
2648
    :param targetHardnessToRound: the hardness to round we want to check.
2649
    :param variable1: the first variable of the polynomial (an expression).
2650
    :param variable2: the second variable of the polynomial (an expression).
2651
    
2652
    :returns: a 4 elements tuple:
2653
                - the polynomial;
2654
                - the modulus (N);
2655
                - the t bound;
2656
                - the lcm used to compute the integral coefficients and the 
2657
                  module.
2658
    """
2659
    # Create a new integer polynomial ring.
2660
    IP = PolynomialRing(ZZ, name=str(variable1) + "," + str(variable2))
2661
    # Coefficients are issued in the increasing power order.
2662
    ratPolyCoefficients = ratPolyOfInt.coefficients()
2663
    # Print the reversed list for debugging.
2664
    #print
2665
    #print "Rational polynomial coefficients:", ratPolyCoefficients[::-1]
2666
    # Build the list of number we compute the lcm of.
2667
    coefficientDenominators = sro_denominators(ratPolyCoefficients)
2668
    #print "Coefficient denominators:", coefficientDenominators
2669
    coefficientDenominators.append(2^precision)
2670
    coefficientDenominators.append(2^(targetHardnessToRound))
2671
    leastCommonMultiple = lcm(coefficientDenominators)
2672
    # Compute the expression corresponding to the new polynomial
2673
    coefficientNumerators =  sro_numerators(ratPolyCoefficients)
2674
    #print coefficientNumerators
2675
    polynomialExpression = 0
2676
    power = 0
2677
    # Iterate over two lists at the same time, stop when the shorter
2678
    # (is this case coefficientsNumerators) is 
2679
    # exhausted. Both lists are ordered in the order of increasing powers
2680
    # of variable1.
2681
    for numerator, denominator in \
2682
                        zip(coefficientNumerators, coefficientDenominators):
2683
        multiplicator = leastCommonMultiple / denominator 
2684
        newCoefficient = numerator * multiplicator
2685
        polynomialExpression += newCoefficient * variable1^power
2686
        power +=1
2687
    polynomialExpression += - variable2
2688
    return (IP(polynomialExpression),
2689
            leastCommonMultiple / 2^precision, # 2^K aka N.
2690
            #leastCommonMultiple / 2^(targetHardnessToRound + 1), # tBound
2691
            leastCommonMultiple / 2^(targetHardnessToRound), # tBound
2692
            leastCommonMultiple) # If we want to make test computations.
2693
        
2694
# End slz_rat_poly_of_int_to_poly_for_coppersmith
2695

    
2696
def slz_rat_poly_of_rat_to_rat_poly_of_int(ratPolyOfRat, 
2697
                                          precision):
2698
    """
2699
    Makes a variable substitution into the input polynomial so that the output
2700
    polynomial can take integer arguments.
2701
    All variables of the input polynomial "have precision p". That is to say
2702
    that they are rationals with denominator == 2^(precision - 1): 
2703
    x = y/2^(precision - 1).
2704
    We "incorporate" these denominators into the coefficients with, 
2705
    respectively, the "right" power.
2706
    """
2707
    polynomialField = ratPolyOfRat.parent()
2708
    polynomialVariable = ratPolyOfRat.variables()[0]
2709
    #print "The polynomial field is:", polynomialField
2710
    return \
2711
        polynomialField(ratPolyOfRat.subs({polynomialVariable : \
2712
                                   polynomialVariable/2^(precision-1)}))
2713
    
2714
# End slz_rat_poly_of_rat_to_rat_poly_of_int
2715

    
2716
def slz_reduce_and_test_base(matrixToReduce,
2717
                             nAtAlpha,
2718
                             monomialsCountSqrt):
2719
    """
2720
    Reduces the basis, tests the Coppersmith condition and returns
2721
    a list of rows that comply with the condition.
2722
    """
2723
    ccComplientRowsList = []
2724
    reducedMatrix = matrixToReduce.LLL(None)
2725
    if not reducedMatrix.is_LLL_reduced():
2726
        raise Exception("reducedMatrix is not actually reduced. Aborting!")
2727
    else:
2728
        print "reducedMatrix is actually reduced."
2729
    print "N^alpha:", nAtAlpha.n()
2730
    rowIndex = 0
2731
    for row in reducedMatrix.rows():
2732
        l2Norm = row.norm(2)
2733
        print "L_2 norm for vector # ", rowIndex, "= ", RR(l2Norm), "*", \
2734
                monomialsCountSqrt.n(), ". Is vector OK?", 
2735
        if (l2Norm * monomialsCountSqrt < nAtAlpha):
2736
            ccComplientRowsList.append(row)
2737
            print "True"
2738
        else:
2739
            print "False"
2740
    # End for
2741
    return ccComplientRowsList
2742
# End slz_reduce_and_test_base
2743

    
2744
def slz_reduce_lll_proj(matToReduce, n):
2745
    """
2746
    Compute the transformation matrix that realizes an LLL reduction on
2747
    the random uniform projected matrix.
2748
    Return both the initial matrix "reduced" by the transformation matrix and
2749
    the transformation matrix itself.
2750
    """
2751
    ## Compute the projected matrix.
2752
    """
2753
    # Random matrix elements {-2^(n-1),...,0,...,2^(n-1)-1}.
2754
    matProjector = slz_random_proj_uniform(n,
2755
                                           matToReduce.ncols(),
2756
                                           matToReduce.nrows())
2757
    """
2758
    reductionOK = False
2759
    while not reductionOK:
2760
        # Random matrix elements in {-1,0,1}. 
2761
        matProjector = slz_random_proj_pm1(matToReduce.ncols(),
2762
                                           matToReduce.nrows())
2763
        matProjected = matToReduce * matProjector
2764
        ## Build the argument matrix for LLL in such a way that the
2765
        #  transformationmatrix is also returned. This matrix is obtained
2766
        #  at almost no extra cost. An identity matrix must be appended to 
2767
        #  the left of the initial matrix. The transformation matrix will
2768
        # will be recovered at the same location from the returned matrix .
2769
        idMat = identity_matrix(matProjected.nrows())
2770
        augmentedMatToReduce = idMat.augment(matProjected)
2771
        try:
2772
            reducedProjMat = \
2773
                augmentedMatToReduce.LLL(algorithm='fpLLL:wrapper')
2774
            reductionOK = True
2775
        except ValueError:
2776
            print "Problem with projected matrix." 
2777
    # Enb while not reduction OK
2778
    #
2779
    ## Recover the transformation matrix (the left part of the reduced matrix). 
2780
    #  We discard the reduced matrix itself.
2781
    transMat = reducedProjMat.submatrix(0,
2782
                                        0,
2783
                                        reducedProjMat.nrows(),
2784
                                        reducedProjMat.nrows())
2785
    ## Return the initial matrix "reduced" and the transformation matrix tuple.
2786
    return (transMat * matToReduce, transMat)
2787
# End  slz_reduce_lll_proj.
2788
#
2789
def slz_reduce_lll_proj_02(matToReduce, n):
2790
    """
2791
    Compute the transformation matrix that realizes an LLL reduction on
2792
    the random uniform projected matrix.
2793
    Return both the initial matrix "reduced" by the transformation matrix and
2794
    the transformation matrix itself.
2795
    """
2796
    ## Compute the projected matrix.
2797
    """
2798
    # Random matrix elements {-2^(n-1),...,0,...,2^(n-1)-1}.
2799
    matProjector = slz_random_proj_uniform(n,
2800
                                           matToReduce.ncols(),
2801
                                           matToReduce.nrows())
2802
    """
2803
    # Random matrix elements in {-8,0,7} -> 4. 
2804
    matProjector = slz_random_proj_uniform(matToReduce.ncols(),
2805
                                           matToReduce.nrows(),
2806
                                           4)
2807
    matProjected = matToReduce * matProjector
2808
    ## Build the argument matrix for LLL in such a way that the transformation
2809
    #  matrix is also returned. This matrix is obtained at almost no extra
2810
    # cost. An identity matrix must be appended to 
2811
    #  the left of the initial matrix. The transformation matrix will
2812
    # will be recovered at the same location from the returned matrix .
2813
    idMat = identity_matrix(matProjected.nrows())
2814
    augmentedMatToReduce = idMat.augment(matProjected)
2815
    reducedProjMat = \
2816
        augmentedMatToReduce.LLL(algorithm='fpLLL:wrapper')
2817
    ## Recover the transformation matrix (the left part of the reduced matrix). 
2818
    #  We discard the reduced matrix itself.
2819
    transMat = reducedProjMat.submatrix(0,
2820
                                        0,
2821
                                        reducedProjMat.nrows(),
2822
                                        reducedProjMat.nrows())
2823
    ## Return the initial matrix "reduced" and the transformation matrix tuple.
2824
    return (transMat * matToReduce, transMat)
2825
# End  slz_reduce_lll_proj_02.
2826
#
2827
def slz_resultant(poly1, poly2, elimVar, debug = False):
2828
    """
2829
    Compute the resultant for two polynomials for a given variable
2830
    and return the (poly1, poly2, resultant) if the resultant
2831
    is not 0.
2832
    Return () otherwise.
2833
    """
2834
    polynomialRing0    = poly1.parent()
2835
    resultantInElimVar = poly1.resultant(poly2,polynomialRing0(elimVar))
2836
    if resultantInElimVar is None:
2837
        if debug:
2838
            print poly1
2839
            print poly2
2840
            print "have GCD = ", poly1.gcd(poly2)
2841
        return None
2842
    if resultantInElimVar.is_zero():
2843
        if debug:
2844
            print poly1
2845
            print poly2
2846
            print "have GCD = ", poly1.gcd(poly2)
2847
        return None
2848
    else:
2849
        if debug:
2850
            print poly1
2851
            print poly2
2852
            print "have resultant in t = ", resultantInElimVar
2853
        return resultantInElimVar
2854
# End slz_resultant.
2855
#
2856
def slz_resultant_tuple(poly1, poly2, elimVar):
2857
    """
2858
    Compute the resultant for two polynomials for a given variable
2859
    and return the (poly1, poly2, resultant) if the resultant
2860
    is not 0.
2861
    Return () otherwise.
2862
    """
2863
    polynomialRing0    = poly1.parent()
2864
    resultantInElimVar = poly1.resultant(poly2,polynomialRing0(elimVar))
2865
    if resultantInElimVar.is_zero():
2866
        return ()
2867
    else:
2868
        return (poly1, poly2, resultantInElimVar)
2869
# End slz_resultant_tuple.
2870
#
2871
def slz_uniform(n):
2872
    """
2873
    Compute a uniform RV in [-2^(n-1), 2^(n-1)-1] (zero is included).
2874
    """
2875
    ## function getrandbits(n) generates a long int with n random bits.
2876
    return getrandbits(n) - 2^(n-1)
2877
# End slz_uniform.
2878
#
2879
sys.stderr.write("\t...sageSLZ loaded\n")
2880