Statistiques
| Branche: | Révision :

cvp / src / functions_for_cvp.sage @ 04d64614

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

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