Statistiques
| Branche: | Révision :

cvp / src / functions_for_cvp.sage @ 04d64614

Historique | Voir | Annoter | Télécharger (25,49 ko)

1
"""@package functions_for_cvp
2
Auxiliary functions used for CVP.
3
"""
4
print "Functions for CVP loading..."
5
#
6
def cvp_affine_from_chebyshev(lowerBound, upperBound):
7
    """
8
    Compute the affine parameters to map [-1, 1] to the interval given
9
    as argument.
10
    @param lowerBound (of the target interval)
11
    @param upperBound (of the target interval)
12
    @return the (scalingCoefficient, offset) tuple.
13
    """
14
    ## Check bounds consistency. Bounds must differ.
15
    if lowerBound >= upperBound:
16
        return None
17
    scalingCoefficient = (upperBound - lowerBound) / 2
18
    offset             = (lowerBound + upperBound) / 2
19
    return (scalingCoefficient, offset)
20
# End cvp_affine_from_chebyshev
21
#
22
def cvp_affine_to_chebyshev(lowerBound, upperBound):
23
    """
24
    Compute the affine parameters to map the interval given
25
    as argument to  [-1, 1].
26
    @param lowerBound (of the target interval)
27
    @param upperBound (of the target interval)
28
    @return the (scalingCoefficient, offset) tuple.
29
    """
30
    ## Check bounds consistency. Bounds must differ.
31
    if lowerBound >= upperBound:
32
        return None
33
    scalingCoefficient = 2 / (upperBound - lowerBound)
34
    ## If bounds are symmetrical with relation to 0, return 0 straight before
35
    #  attempting a division by 0.
36
    if lowerBound == -upperBound:
37
        offset         = 0
38
    else:
39
        offset         =  (lowerBound + upperBound) / (lowerBound - upperBound)
40
    return (scalingCoefficient, offset)
41
# End cvp_affine_to_chebyshev
42
#
43
def cvp_babai(redBasis, redBasisGso, vect):
44
    """
45
    Closest plane vector implementation as per Babaï.
46
    @param redBasis   : a lattice basis, preferably a reduced one; 
47
    @param redBasisGSO: the GSO of the previous basis;
48
    @param vector     : a vector, in coordinated in the ambient 
49
                        space of the lattice
50
    @return the closest vector to the input, in coordinates in the 
51
                        ambient space of the lattice.
52
    """
53
    ## A deep copy is not necessary here.
54
    curVect = copy(vect)
55
    print "cvp_babai - Vector:", vect
56
    ## From index of last row down to 0.
57
    for vIndex in xrange(len(redBasis.rows())-1, -1, -1):
58
        curRBGSO = redBasisGso.row(vIndex)
59
        curVect = curVect - \
60
                  (round((curVect * curRBGSO)  / (curRBGSO * curRBGSO)) * \
61
                   redBasis.row(vIndex)) 
62
    return vect - curVect
63
# End cvp_babai
64
#
65
def cvp_bkz_reduce(intBasis, blockSize):
66
    """
67
    Thin and simplistic wrapper the HKZ function of the FP_LLL module.
68
    """
69
    fplllIntBasis = FP_LLL(intBasis)
70
    fplllIntBasis.BKZ(blockSize)
71
    return fplllIntBasis._sage_()
72
# End cvp_hkz_reduce.
73
#
74
def cvp_coefficients_bounds_projection(executablePath, arguments):
75
    """
76
    Compute the min and max of the coefficients with linear
77
    programming projection.
78
    @param executablePath: the path to the binary program;
79
    @param arguments: a list of arguments to be givent to the binary
80
    @return: the min and max coefficients value arrays (in this order).
81
    """
82
    from subprocess import check_output
83
    commandLine = [executablePath] + arguments
84
    ## Get rid of output on stderr.
85
    devNull = open("/dev/null", "rw")
86
    ## Run ther program.
87
    otp = check_output(commandLine, stderr=devNull)
88
    devNull.close()
89
    ## Recover the results
90
    bounds = sage_eval(otp)
91
    minCoeffsExpList = []
92
    maxCoeffsExpList = [] 
93
    print "bounds:", bounds
94
    for boundsPair in bounds:
95
        minCoeffsExpList.append(boundsPair[0])
96
        maxCoeffsExpList.append(boundsPair[1])
97
    print "minCoeffsExpList:", minCoeffsExpList
98
    print "maxCoeffsExpList:", maxCoeffsExpList
99
    return (minCoeffsExpList, maxCoeffsExpList)
100
# End cvp_coefficients_bounds_projection
101

    
102
def cvp_coefficients_bounds_projection_exps(executablePath, 
103
                                            arguments, 
104
                                            realField = None):
105
    """
106
    Compute the min and max exponents of the coefficients with linear
107
    programming projection.
108
    @param executablePath: the path to the binary program;
109
    @param arguments: a list of arguments to be givent to the binary
110
    @param realField: the realField to use for floating-point conversion.
111
    @return: the min and max exponents arrays (in this order).
112
    """
113
    ## If no realField is givne, fall back on RR.
114
    if realField is None:
115
        realField = RR
116
    (minCoeffsExpList, maxCoeffsExpList) = \
117
        cvp_coefficients_bounds_projection(executablePath,arguments)
118
    for index in xrange(0, len(minCoeffsExpList)):
119
        minCoeffsExpList[index] = \
120
            realField(minCoeffsExpList[index]).log2().floor()
121
        maxCoeffsExpList[index] = \
122
            realField(maxCoeffsExpList[index]).log2().floor()
123
    return (minCoeffsExpList, maxCoeffsExpList)
124
# End cvp_coefficients_bounds_projection_exps
125

    
126
def cvp_chebyshev_zeros_1k(numPoints, realField, contFracMaxErr = None):
127
    """
128
    Compute the Chebyshev zeros for some polynomial degree (numPoints, for the
129
    zeros) scaled to the [lowerBound, upperBound] interval.
130
    The list is returned as row floating-point numbers is contFracMaxErr is None.
131
    Otherwise elements are transformed into rational numbers. 
132
    """
133
    if numPoints < 1:
134
        return None
135
    if numPoints == 0:
136
        return [0]
137
    ## Check that realField has a "prec()" member.
138
    try:
139
        realField.prec()
140
    except:
141
        return None
142
    #
143
    zerosList = []
144
    ## formulas: cos (((2*i - 1) * pi) / (2*numPoints)) for i = 1,...,numPoints.
145
    #  If i is -1-shifted, as in the following loop, formula is 
146
    #  cos (((2*i + 1) * pi) / (2*numPoints)) for i = 0,...,numPoints-1.
147
    for index in xrange(0, numPoints):
148
        ## When number of points is odd (then, the Chebyshev degree is odd two), 
149
        #  the central point is 0. */
150
        if (numPoints % 2 == 1) and ((2 * index + 1) == numPoints):
151
            if contFracMaxErr is None:
152
                zerosList.append(realField(0))
153
            else:
154
                zerosList.append(0)
155
        else:
156
            currentCheby = \
157
                ((realField(2*index+1) * realField(pi)) / 
158
                 realField(2*numPoints)).cos()
159
            #print  "Current Cheby:", currentCheby
160
            ## Result is negated to have an increasing order list.
161
            if contFracMaxErr is None:
162
                zerosList.append(-currentCheby)
163
            else:
164
                zerosList.append(-currentCheby.nearby_rational(contFracMaxErr))
165
            #print "Relative error:", (currentCheby/zerosList[index])
166
    return zerosList
167
# End cvp_chebyshev_zeros_1k.
168
#
169
def cvp_chebyshev_zeros_for_interval(lowerBound, 
170
                                     upperBound, 
171
                                     numPoints,
172
                                     realField = None, 
173
                                     contFracMaxErr = None):
174
    """
175
    Compute the Chebyshev zeros for some polynomial degree (numPoints, for the
176
    zeros) scaled to the [lowerBound, upperBound] interval.
177
    If no contFracMaxErr argument is given, we return the list as "raw"
178
    floating-points. Otherwise, rational numbers are returned.
179
    """
180
    ## Arguments check.
181
    if lowerBound >= upperBound:
182
        return None
183
    if numPoints < 1:
184
        return None
185
    ## If no realField argument is given try to retrieve it from the 
186
    #  bounds. If impossible (e.g. rational bound), fall back on RR. 
187
    if realField is None:
188
        try:
189
            ### Force failure if parent does not have prec() member.
190
            lowerBound.parent().prec()
191
            ### If no exception is thrown, set realField.
192
            realField = lowerBound.parent()
193
        except:
194
            realField = RR
195
    #print "Real field:", realField
196
    ## We want "raw"floating-point nodes to only have one final rounding.
197
    chebyshevZerosList = \
198
        cvp_chebyshev_zeros_1k(numPoints, realField)
199
    scalingFactor, offset = cvp_affine_from_chebyshev(lowerBound, upperBound)
200
    for index in xrange(0, len(chebyshevZerosList)):
201
        chebyshevZerosList[index] = \
202
            chebyshevZerosList[index] * scalingFactor + offset
203
        if not contFracMaxErr is None:
204
            chebyshevZerosList[index] = \
205
                chebyshevZerosList[index].nearby_rational(contFracMaxErr)
206
    return chebyshevZerosList                                                             
207
# End cvp_chebyshev_zeros_for_interval.
208
#
209
def cvp_common_scaling_exp(precision, 
210
                           precisionFraction = None):
211
    """
212
    Compute the common scaling factor (a power of 2).
213
    The exponent is a fraction of the precision;
214
    A extra factor is computed for the vector,
215
    exra factors are computed for the basis vectors, all in the corresponding
216
    functions.
217
    """
218
    ## Set a default value for the precisionFraction to 1/2.
219
    if precisionFraction is None:
220
        precisionFraction = 3/4
221
    """
222
    minBasisVectsNorm = oo
223
    currentNorm = 0
224
    for vect in floatBasis.rows():
225
        currentNorm = vect.norm()
226
        print "Current norm:", currentNorm
227
        if currentNorm < minBasisVectsNorm:
228
            minBasisVectsNorm = currentNorm
229
    currentNorm = funcVect.norm()
230
    print "Current norm:", currentNorm
231
    if currentNorm < minBasisVectsNorm:
232
        minBasisVectsNorm = currentNorm
233
    ## Check if a push is need because the smallest norm is below 0.    
234
    powerForMinNorm = floor(log(minBasisVectsNorm)/log2)
235
    print "Power for min norm :", powerForMinNorm
236
    print "Power for precision:", ceil(precision*precisionFraction)
237
    if powerForMinNorm < 0:
238
        return 2^(ceil(precision*precisionFraction) - powerFromMinNorm)
239
    else:
240
    """
241
    return ceil(precision * precisionFraction)
242
# End cvp_common_scaling_factor.
243
#
244
def cvp_evaluation_points_list(funct, 
245
                               maxDegree, 
246
                               lowerBound, 
247
                               upperBound, 
248
                               realField = None):
249
    """
250
    Compute a list of points for the vector creation.
251
    Strategy:
252
    - compute the zeros of the function-remez;
253
    - compute the zeros of the function-rounded_remez;
254
    - compute the Chebyshev points.
255
    Merge the whole thing.
256
    """
257
    
258
# End cvp_evaluation_points_list.
259
#
260
def cvp_extra_basis_scaling_exps(precsList, minExponentsList):
261
    extraScalingExpsList = []
262
    for minExp, prec in zip(minExponentsList, precsList):
263
        if minExp - prec < 0:
264
            extraScalingExpsList.append(minExp - prec)
265
        else:
266
            extraScalingExpsList.append(0)
267
    return extraScalingExpsList
268
# End cvp_extra_basis_scaling_exps.
269
#
270
def cvp_extra_function_scaling_exp(floatBasis,
271
                                   floatFuncVect,
272
                                   maxExponentsList):
273
    """
274
    Compute the extra scaling exponent for the function vector in ordre to 
275
    guarantee that the maximum exponent can be reached for each element of the
276
    basis.
277
    """
278
    ## Check arguments
279
    if floatBasis.nrows() == 0 or \
280
       floatBasis.ncols() == 0 or \
281
       len(floatFuncVect)      == 0:
282
        return 1
283
    if len(maxExponentsList) != floatBasis.nrows():
284
        return None
285
    ## Compute the log of the norm of the function vector.
286
    funcVectNormLog  = log(floatFuncVect.norm())
287
    ## Compute the extra scaling factor for each vector of the basis.
288
    #  for norms ratios an maxExponentsList. 
289
    extraScalingExp = 0
290
    for (row, maxExp) in zip(floatBasis.rows(),maxExponentsList):
291
        rowNormLog     = log(row.norm())
292
        normsRatioLog2 = floor((funcVectNormLog - rowNormLog) / log2) 
293
        #print "Current norms ratio log2:", normsRatioLog2
294
        scalingExpCandidate = normsRatioLog2 - maxExp 
295
        #print "Current scaling exponent candidate:", scalingExpCandidate
296
        if scalingExpCandidate < 0:
297
            if -scalingExpCandidate > extraScalingExp:
298
                extraScalingExp = -scalingExpCandidate
299
    return extraScalingExp
300
# End cvp_extra_function_scaling_exp.
301
#
302
def cvp_float_basis(monomialsList, pointsList, realField):
303
    """
304
    For a given monomials list and points list, compute the floating-point
305
    basis matrix.
306
    """
307
    numRows    = len(monomialsList)
308
    numCols    = len(pointsList)
309
    if numRows == 0 or numCols == 0:
310
        return matrix(realField, 0, 0)
311
    #
312
    floatBasis = matrix(realField, numRows, numCols)
313
    for rIndex in xrange(0, numRows):
314
        for cIndex in xrange(0, numCols):
315
            floatBasis[rIndex,cIndex] = \
316
                monomialsList[rIndex](realField(pointsList[cIndex]))
317
    return floatBasis
318
# End cvp_float_basis
319
#
320
def cvp_float_vector_for_approx(func, 
321
                                basisPointsList,
322
                                realField):
323
    """
324
    Compute a floating-point vector for the function and the points list.
325
    """
326
    #
327
    ## Some local variables.
328
    basisPointsNum = len(basisPointsList)
329
    #
330
    ##
331
    vVect = vector(realField, basisPointsNum)
332
    for vIndex in xrange(0,basisPointsNum):
333
        vVect[vIndex] = \
334
            func(basisPointsList[vIndex])
335
    return vVect
336
# End cvp_float_vector_for_approx.
337
#
338
def cvp_hkz_reduce(intBasis):
339
    """
340
    Thin and simplistic wrapper the HKZ function of the FP_LLL module.
341
    """
342
    fplllIntBasis = FP_LLL(intBasis)
343
    fplllIntBasis.HKZ()
344
    return fplllIntBasis._sage_()
345
# End cvp_hkz_reduce.
346
#
347
def cvp_int_basis(floatBasis, 
348
                  commonScalingExp, 
349
                  extraScalingExpsList):
350
    """
351
    From the float basis, the common scaling factor and the extra exponents 
352
    compute the integral basis.
353
    """
354
    ## Checking arguments.
355
    if floatBasis.nrows() != len(extraScalingExpsList):
356
        return None
357
    intBasis = matrix(ZZ, floatBasis.nrows(), floatBasis.ncols())
358
    for rIndex in xrange(0, floatBasis.nrows()):
359
        for cIndex in xrange(0, floatBasis.ncols()):
360
            intBasis[rIndex, cIndex] = \
361
            floatBasis[rIndex, cIndex] * \
362
            2^(commonScalingExp + extraScalingExpsList[rIndex])
363
    return intBasis
364
# End cvp_int_basis.
365
#
366
def cvp_int_vector_for_approx(floatVect, commonScalingExp, extraScalingExp):
367
    totalScalingFactor = 2^(commonScalingExp + extraScalingExp)
368
    intVect = vector(ZZ, len(floatVect))
369
    for index in xrange(0, len(floatVect)):
370
        intVect[index] = (floatVect[index] * totalScalingFactor).round()
371
    return intVect
372
# End cvp_int_vector_for_approx.
373
#
374
def cvp_lll_reduce(intBasis):
375
    """
376
    Thin and simplistic wrapper around the LLL function
377
    """
378
    return intBasis.LLL()
379
# End cvp_lll_reduce.
380
#
381
def cvp_monomials_list(exponentsList, varName = None):
382
    """
383
    Create a list of monomials corresponding to the given exponentsList.
384
    Monomials are defined as functions.
385
    """
386
    monomialsList = []
387
    for exponent in exponentsList:
388
        if varName is None:
389
            monomialsList.append((x^exponent).function(x))
390
        else:
391
            monomialsList.append((varName^exponent).function(varName))
392
    return monomialsList
393
# End cvp_monomials_list.
394
#
395
def cvp_polynomial_coeffs_from_vect(cvpVectInBasis,
396
                                   coeffsPrecList, 
397
                                   minCoeffsExpList,
398
                                   extraFuncScalingExp):
399
    numElements = len(cvpVectInBasis)
400
    ## Arguments check.
401
    if numElements != len(minCoeffsExpList) or \
402
       numElements != len(coeffsPrecList):
403
        return None
404
    polynomialCoeffsList = []
405
    for index in xrange(0, numElements):
406
        currentCoeff = cvp_round_to_bits(cvpVectInBasis[index],
407
                                         coeffsPrecList[index])
408
        print "cvp_polynomial_coeffs - current coefficient rounded:", currentCoeff
409
        currentCoeff *= 2^(minCoeffsExpList[index])
410
        print "cvp_polynomial_coeffs - current coefficient scaled :", currentCoeff
411
        currentCoeff *= 2^(-coeffsPrecList[index])
412
        print "cvp_polynomial_coeffs - current coefficient prec corrected :", currentCoeff
413
        currentCoeff *= 2^(-extraFuncScalingExp)
414
        print "cvp_polynomial_coeffs - current coefficient extra func scaled :", currentCoeff
415
        polynomialCoeffsList.append(currentCoeff)
416
    return polynomialCoeffsList
417
# En polynomial_coeffs_from_vect.
418
#
419
def cvp_remez_all_poly_error_func_zeros(funct, 
420
                                        maxDegree, 
421
                                        lowerBound, 
422
                                        upperBound,
423
                                        precsList, 
424
                                        realField = None,
425
                                        contFracMaxErr = None):
426
    """
427
    For a given function f, a degree and an interval, compute the 
428
    zeros of the (f-remez_d(f)) function and those of the (f-truncRemez_d(f)) 
429
    function over the interval. If the (f-truncRemez_d(f)) function has very 
430
    few zeros, compute the zeros of the derivative instead.
431
    """
432
    ## If no realField argument is given try to retrieve it from the 
433
    #  bounds. If impossible (e.g. rational bound), fall back on RR. 
434
    if realField is None:
435
        try:
436
            ### Force failure if parent does not have prec() member.
437
            lowerBound.parent().prec()
438
            ### If no exception is thrown, set realField.
439
            realField = lowerBound.parent()
440
        except:
441
            realField = RR
442
    #print "Real field:", realField
443
    funcSa         = func
444
    funcAsStringSa = func._assume_str().replace("_SAGE_VAR_","",100)
445
    print "Function as a string:", funcAsStringSa
446
    ## Create the Sollya version of the function and compute the
447
    #  Remez approximant
448
    funcSo  = pobyso_parse_string(funcAsStringSa)
449
    pStarSo = pobyso_remez_canonical_sa_so(funcAsStringSa, 
450
                                           maxDegree, 
451
                                           lowerBound,
452
                                           upperBound)
453
    ## Compute the zeroes of the error function.
454
    ### First create the error function with copies since they are needed later.
455
    errorFuncSo = sollya_lib_build_function_sub(sollya_lib_copy_obj(pStarSo), 
456
                                                sollya_lib_copy_obj(funcSo))
457
    intervalSo = pobyso_range_from_bounds_sa_so(lowerBound, upperBound)
458
    errorFuncZerosSo = pobyso_dirty_find_zeros_so_so(errorFuncSo, intervalSo)
459
    pobyso_clear_obj(errorFuncSo) # Clears funcSo and pStarSo as well.
460
    #print "Zeroes of the error function from Sollya:"
461
    #pobyso_autoprint(errorFuncZerosSo)
462
    errorFuncZerosSa = pobyso_float_list_so_sa(errorFuncZerosSo)
463
    pobyso_clear_obj(errorFuncZerosSo)
464
    #print "\nZeroes of the error function from Sage:"
465
    #print errorFuncZerosSa
466
    #
467
    ## Deal with the truncated polynomial now.
468
    ### Create the formats list. Notice the "*" before the list variable name.
469
    truncFormatListSo = pobyso_build_list_of_ints_sa_so(*precsList)
470
    print "Precisions list as Sollya object:", 
471
    pobyso_autoprint(truncFormatListSo)
472
    pTruncSo = pobyso_round_coefficients_so_so(pStarSo, truncFormatListSo)
473
    print "Truncated polynomial:", ; pobyso_autoprint(pTruncSo)
474
    ### Compute the error function. This time we consume both the function
475
    #   and the polynomial.
476
    errorFuncSo = sollya_lib_build_function_sub(pTruncSo, 
477
                                                funcSo)
478
    errorFuncZerosSo = pobyso_dirty_find_zeros_so_so(errorFuncSo, intervalSo)
479
    pobyso_clear_obj(pStarSo)
480
    print "Zeroes of the error function for the truncated polynomial from Sollya:"
481
    pobyso_autoprint(errorFuncZerosSo)
482
    errorFuncTruncZerosSa = pobyso_float_list_so_sa(errorFuncZerosSo)
483
    pobyso_clear_obj(errorFuncZerosSo)
484
    ## It may happen that the error function has only one or two zeros.
485
    #  In this case, we are interested in the variations and we consider the
486
    #  derivative
487
    if maxDegree > 3:
488
        if len(errorFuncTruncZerosSa) <= 2:
489
            errorFuncDiffSo = pobyso_diff_so_so(errorFuncSo)
490
            errorFuncZerosSo = pobyso_dirty_find_zeros_so_so(errorFuncDiffSo,
491
                                                             intervalSo)
492
            errorFuncTruncZerosSa = pobyso_float_list_so_sa(errorFuncZerosSo)
493
            pobyso_clear_obj(errorFuncZerosSo)
494
            pobyso_clear_obj(errorFuncDiffSo) # Clears funcSo and pTruncSo as well.
495
    pobyso_clear_obj(errorFuncSo) # Clears funcSo and pTruncSo as well.
496
    ## Sollya objects clean up.
497
    pobyso_clear_obj(intervalSo)
498
    errorFuncZerosSa += errorFuncTruncZerosSa
499
    errorFuncZerosSa.sort()
500
    ## If required, convert the numbers to rational numbers.
501
    if not contFracMaxErr is None:
502
        for index in xrange(0, len(errorFuncZerosSa)):
503
            errorFuncZerosSa[index] = \
504
                errorFuncZerosSa[index].nearby_rational(contFracMaxErr)
505
    return errorFuncZerosSa
506
    #
507
# End remez_all_poly_error_func_zeros
508
#
509
def cvp_remez_poly_error_func_zeros(funct, 
510
                                    maxDegree, 
511
                                    lowerBound, 
512
                                    upperBound,
513
                                    realField = None,
514
                                    contFracMaxErr = None):
515
    """
516
    For a given function f, a degree and an interval, compute the 
517
    zeros of the (f-remez_d(f)) function.
518
    """
519
    ## If no realField argument is given try to retrieve it from the 
520
    #  bounds. If impossible (e.g. rational bound), fall back on RR. 
521
    if realField is None:
522
        try:
523
            ### Force failure if parent does not have prec() member.
524
            lowerBound.parent().prec()
525
            ### If no exception is thrown, set realField.
526
            realField = lowerBound.parent()
527
        except:
528
            realField = RR
529
    #print "Real field:", realField
530
    funcSa         = func
531
    funcAsStringSa = func._assume_str().replace("_SAGE_VAR_","",100)
532
    #print "Function as a string:", funcAsStringSa
533
    ## Create the Sollya version of the function and compute the
534
    #  Remez approximant
535
    funcSo  = pobyso_parse_string(funcAsStringSa)
536
    pStarSo = pobyso_remez_canonical_sa_so(funcAsStringSa, 
537
                                           maxDegree, 
538
                                           lowerBound,
539
                                           upperBound)
540
    ## Compute the zeroes of the error function.
541
    ### Error function creation consumes both pStarSo and funcSo.
542
    errorFuncSo = sollya_lib_build_function_sub(pStarSo, funcSo)
543
    intervalSo = pobyso_range_from_bounds_sa_so(lowerBound, upperBound)
544
    errorFuncZerosSo = pobyso_dirty_find_zeros_so_so(errorFuncSo, intervalSo)
545
    pobyso_clear_obj(errorFuncSo) # Clears funcSo and pStarSo as well.
546
    #print "Zeroes of the error function from Sollya:"
547
    #pobyso_autoprint(errorFuncZerosSo)
548
    errorFuncZerosSa = pobyso_float_list_so_sa(errorFuncZerosSo)
549
    pobyso_clear_obj(errorFuncZerosSo)
550
    ## Sollya objects clean up.
551
    pobyso_clear_obj(intervalSo)
552
    ## Converting to rational numbers, if required.
553
    if not contFracMaxErr is None:
554
        for index in xrange(0, len(errorFuncZerosSa)):
555
            errorFuncZerosSa[index] = \
556
                errorFuncZerosSa[index].nearby_rational(contFracMaxErr)
557
    return errorFuncZerosSa
558
    #
559
# End remez_poly_error_func_zeros
560
#
561
def cvp_round_to_bits(num, numBits):
562
    """
563
    Round "num" to "numBits" bits.
564
    @param num    : the number to round;
565
    @param numBits: the number of bits to round to.
566
    """
567
    if num == 0:
568
        return num
569
    log2ofNum  = 0
570
    numRounded = 0
571
    ## Avoid conversion to floating-point if not necessary.
572
    try:
573
        log2OfNum = num.abs().log2()
574
    except:
575
        log2OfNum = RR(num).abs().log2()
576
    ## Works equally well for num >= 1 and num < 1.
577
    log2OfNum = ceil(log2OfNum)
578
    # print "cvp_round_to_bits - Log2OfNum:", log2OfNum
579
    divMult = log2OfNum - numBits
580
    #print "cvp_round_to_bits - DivMult:", divMult
581
    if divMult != 0:
582
        numRounded = round(num/2^divMult) * 2^divMult
583
    else:
584
        numRounded = num
585
    return numRounded
586
# End cvp_round_to_bits
587
#
588
def cvp_vector_in_basis(vect, basis):
589
    """
590
    Compute the coordinates of "vect" in "basis" by
591
    solving a linear system.
592
    @param vect : the vector we want to compute the coordinates of
593
                  in coordinates of the ambient space;
594
    @param basis: the basis we want to compute the coordinates in
595
                  as a matrix relative to the ambient space.
596
    """
597
    ## Create the variables for the linear equations.
598
    varDeclString = ""
599
    basisIterationsRange = xrange(0, basis.nrows())
600
    ### Building variables declaration string.
601
    for index in basisIterationsRange:
602
        varName = "a" + str(index)
603
        if varDeclString == "":
604
            varDeclString += varName
605
        else:
606
            varDeclString += "," + varName
607
    ### Variables declaration
608
    varsList = var(varDeclString)
609
    sage_eval("var('" + varDeclString + "')")
610
    ## Create the equations
611
    equationString = ""
612
    equationsList = []
613
    for vIndex in xrange(0, len(vect)):
614
        equationString = ""
615
        for bIndex in basisIterationsRange:
616
            if equationString != "":
617
                equationString += "+"
618
            equationString += str(varsList[bIndex]) + "*" + str(basis[bIndex,vIndex])
619
        equationString += " == " + str(vect[vIndex])
620
        equationsList.append(sage_eval(equationString))
621
    ## Solve the equations system. The solution dictionnary is needed
622
    #  to recover the values of the solutions.
623
    solutions = solve(equationsList,varsList, solution_dict=True)
624
    #print "Solutions:", s
625
    ## Recover solutions in rational components vector.
626
    vectInBasis = vector(QQ, basis.nrows())
627
    ### There can be several solution, in the general case.
628
    #   Here, there is only one. For each solution, each 
629
    #   variable has to be searched for.
630
    for sol in solutions:
631
        for bIndex in basisIterationsRange:
632
            vectInBasis[bIndex] = sol[varsList[bIndex]]
633
    return vectInBasis
634
# End cpv_vector_in_basis.
635
#
636
print "... functions for CVP loaded."