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." |