root / pobysoPythonSage / src / pobyso.py @ 284
Historique | Voir | Annoter | Télécharger (98,23 ko)
1 |
"""
|
---|---|
2 |
@file pobyso.py
|
3 |
Actual functions to use in Sage
|
4 |
@author S.T.
|
5 |
@date 2012-11-13
|
6 |
|
7 |
Command line syntax:
|
8 |
use from Sage (via the "load" or the "attach" commands)
|
9 |
|
10 |
pobyso functions come in five flavors:
|
11 |
- the _so_so (arguments and returned objects are pointers to Sollya objects,
|
12 |
includes the void function and the no arguments function that return a
|
13 |
pointer to a Sollya object);
|
14 |
- the _so_sa (argument are pointers to Sollya objects, returned objects are
|
15 |
Sage/Python objects or, more generally, information is transfered from the
|
16 |
Sollya world to Sage/Python world; e.g. functions without arguments that
|
17 |
return a Sage/Python object);
|
18 |
- the _sa_so (arguments are Sage/Python objects, returned objects are
|
19 |
pointers to Sollya objects);
|
20 |
- the sa_sa (arguments and returned objects are all Sage/Python objects);
|
21 |
- a catch all flavor, without any suffix, (e. g. functions that have no
|
22 |
argument nor return value).
|
23 |
|
24 |
This classification is not always very strict. Conversion functions from Sollya
|
25 |
to Sage/Python are sometimes decorated with Sage/Python arguments to set
|
26 |
the precision. These functions remain in the so_sa category.
|
27 |
|
28 |
@note
|
29 |
Reported errors in Eclipse come from the calls to the Sollya library
|
30 |
|
31 |
ToDo (among other things):
|
32 |
-memory management.
|
33 |
"""
|
34 |
from ctypes import * |
35 |
import re |
36 |
from sage.symbolic.expression_conversions import polynomial |
37 |
from sage.symbolic.expression_conversions import PolynomialConverter |
38 |
"""
|
39 |
Create the equivalent to an enum for the Sollya function types and other
|
40 |
common types (error types...).
|
41 |
"""
|
42 |
(SOLLYA_BASE_FUNC_ABS, |
43 |
SOLLYA_BASE_FUNC_ACOS, |
44 |
SOLLYA_BASE_FUNC_ACOSH, |
45 |
SOLLYA_BASE_FUNC_ADD, |
46 |
SOLLYA_BASE_FUNC_ASIN, |
47 |
SOLLYA_BASE_FUNC_ASINH, |
48 |
SOLLYA_BASE_FUNC_ATAN, |
49 |
SOLLYA_BASE_FUNC_ATANH, |
50 |
SOLLYA_BASE_FUNC_CEIL, |
51 |
SOLLYA_BASE_FUNC_CONSTANT, |
52 |
SOLLYA_BASE_FUNC_COS, |
53 |
SOLLYA_BASE_FUNC_COSH, |
54 |
SOLLYA_BASE_FUNC_DIV, |
55 |
SOLLYA_BASE_FUNC_DOUBLE, |
56 |
SOLLYA_BASE_FUNC_DOUBLEDOUBLE, |
57 |
SOLLYA_BASE_FUNC_DOUBLEEXTENDED, |
58 |
SOLLYA_BASE_FUNC_ERF, |
59 |
SOLLYA_BASE_FUNC_ERFC, |
60 |
SOLLYA_BASE_FUNC_EXP, |
61 |
SOLLYA_BASE_FUNC_EXP_M1, |
62 |
SOLLYA_BASE_FUNC_FLOOR, |
63 |
SOLLYA_BASE_FUNC_FREE_VARIABLE, |
64 |
SOLLYA_BASE_FUNC_HALFPRECISION, |
65 |
SOLLYA_BASE_FUNC_LIBRARYCONSTANT, |
66 |
SOLLYA_BASE_FUNC_LIBRARYFUNCTION, |
67 |
SOLLYA_BASE_FUNC_LOG, |
68 |
SOLLYA_BASE_FUNC_LOG_10, |
69 |
SOLLYA_BASE_FUNC_LOG_1P, |
70 |
SOLLYA_BASE_FUNC_LOG_2, |
71 |
SOLLYA_BASE_FUNC_MUL, |
72 |
SOLLYA_BASE_FUNC_NEARESTINT, |
73 |
SOLLYA_BASE_FUNC_NEG, |
74 |
SOLLYA_BASE_FUNC_PI, |
75 |
SOLLYA_BASE_FUNC_POW, |
76 |
SOLLYA_BASE_FUNC_PROCEDUREFUNCTION, |
77 |
SOLLYA_BASE_FUNC_QUAD, |
78 |
SOLLYA_BASE_FUNC_SIN, |
79 |
SOLLYA_BASE_FUNC_SINGLE, |
80 |
SOLLYA_BASE_FUNC_SINH, |
81 |
SOLLYA_BASE_FUNC_SQRT, |
82 |
SOLLYA_BASE_FUNC_SUB, |
83 |
SOLLYA_BASE_FUNC_TAN, |
84 |
SOLLYA_BASE_FUNC_TANH, |
85 |
SOLLYA_BASE_FUNC_TRIPLEDOUBLE, |
86 |
SOLLYA_ABSOLUTE, |
87 |
SOLLYA_RELATIVE) = map(int,xrange(46)) |
88 |
sys.stderr.write("SOLLYA_RELATIVE = " + str(SOLLYA_RELATIVE) + "\n") |
89 |
sys.stderr.write("Superficial pobyso check...\n")
|
90 |
#print "First constant - SOLLYA_BASE_FUNC_ABS: ", SOLLYA_BASE_FUNC_ABS
|
91 |
#print "Last constant - SOLLYA_BASE_FUNC_TRIPLEDOUBLE: ", SOLLYA_BASE_FUNC_TRIPLEDOUBLE
|
92 |
|
93 |
pobyso_max_arity = 9
|
94 |
|
95 |
def pobyso_absolute_so_so(): |
96 |
"""
|
97 |
Create an "absolute" Sollya object.
|
98 |
"""
|
99 |
return(sollya_lib_absolute(None)) |
100 |
|
101 |
def pobyso_autoprint(arg): |
102 |
sollya_lib_autoprint(arg, None)
|
103 |
|
104 |
def pobyso_autoprint_so_so(arg): |
105 |
sollya_lib_autoprint(arg,None)
|
106 |
|
107 |
def pobyso_bounds_to_interval_sa_sa(lowerBound, upperBound): |
108 |
"""
|
109 |
Convert a pair of bounds into an interval (an element of
|
110 |
a RealIntervalField).
|
111 |
"""
|
112 |
# Minimal (not bullet-proof) check on bounds.
|
113 |
if lowerBound > upperBound:
|
114 |
return None |
115 |
# Try to get the maximum precision among the bounds.
|
116 |
try:
|
117 |
preclb = parent(lowerBound).precision() |
118 |
precub = parent(upperBound).precision() |
119 |
prec = max(preclb, precub)
|
120 |
except AttributeError: |
121 |
prec = 53
|
122 |
# Create the RealIntervalField and the interval (if possible).
|
123 |
theRIF = RealIntervalField(prec) |
124 |
try:
|
125 |
interval = theRIF(lowerBound, upperBound) |
126 |
except TypeError: |
127 |
return None |
128 |
else:
|
129 |
return interval
|
130 |
# End pobyso_bounds_to_interval_sa_sa
|
131 |
|
132 |
def pobyso_bounds_to_range_sa_so(rnLowerBoundSa, rnUpperBoundSa, \ |
133 |
precisionSa=None):
|
134 |
"""
|
135 |
Return a Sollya range from to 2 RealField Sage elements.
|
136 |
The Sollya range element has a sufficient precision to hold all
|
137 |
the digits of the widest of the Sage bounds.
|
138 |
"""
|
139 |
# Sanity check.
|
140 |
if rnLowerBoundSa > rnUpperBoundSa:
|
141 |
return None |
142 |
# Precision stuff.
|
143 |
if precisionSa is None: |
144 |
# Check for the largest precision.
|
145 |
lbPrecSa = rnLowerBoundSa.parent().precision() |
146 |
ubPrecSa = rnLowerBoundSa.parent().precision() |
147 |
maxPrecSa = max(lbPrecSa, ubPrecSa)
|
148 |
else:
|
149 |
maxPrecSa = precisionSa |
150 |
# From Sage to Sollya bounds.
|
151 |
# lowerBoundSo = sollya_lib_constant(get_rn_value(rnLowerBoundSa),
|
152 |
# maxPrecSa)
|
153 |
lowerBoundSo = pobyso_constant_sa_so(rnLowerBoundSa, |
154 |
maxPrecSa) |
155 |
upperBoundSo = pobyso_constant_sa_so(rnUpperBoundSa, |
156 |
maxPrecSa) |
157 |
|
158 |
# From Sollya bounds to range.
|
159 |
rangeSo = sollya_lib_range(lowerBoundSo, upperBoundSo) |
160 |
# Back to original precision.
|
161 |
# Clean up
|
162 |
sollya_lib_clear_obj(lowerBoundSo) |
163 |
sollya_lib_clear_obj(upperBoundSo) |
164 |
return rangeSo
|
165 |
# End pobyso_bounds_to_range_sa_so
|
166 |
|
167 |
def pobyso_build_end_elliptic_list_so_so(*args): |
168 |
"""
|
169 |
From Sollya argument objects, create a Sollya end elliptic list.
|
170 |
Elements of the list are "eaten" (should not be cleared individually,
|
171 |
are cleared when the list is cleared).
|
172 |
"""
|
173 |
if len(args) == 0: |
174 |
## When called with an empty list, sollya_lib_build_end_elliptic_list,
|
175 |
# produces "error".
|
176 |
return sollya_lib_build_list(None) |
177 |
index = 0
|
178 |
## One can not append elements to an elliptic list, prepend only is
|
179 |
# permitted.
|
180 |
for argument in reversed(args): |
181 |
if index == 0: |
182 |
listSo = sollya_lib_build_end_elliptic_list(argument, None)
|
183 |
else:
|
184 |
listSo = sollya_lib_prepend(argument, listSo) |
185 |
index += 1
|
186 |
return listSo
|
187 |
|
188 |
# End pobyso_build_end_elliptic_list_so_so
|
189 |
|
190 |
def pobyso_build_function_sub_so_so(exp1So, exp2So): |
191 |
return sollya_lib_build_function_sub(exp1So, exp2So)
|
192 |
|
193 |
def pobyso_build_list_of_ints_sa_so(*args): |
194 |
"""
|
195 |
Build a Sollya list from Sage integral arguments.
|
196 |
"""
|
197 |
if len(args) == 0: |
198 |
return pobyso_build_list_so_so()
|
199 |
## Make a Sage list of integral constants.
|
200 |
intsList = [] |
201 |
for intElem in args: |
202 |
intsList.append(pobyso_constant_from_int_sa_so(intElem)) |
203 |
return pobyso_build_list_so_so(*intsList)
|
204 |
|
205 |
def pobyso_build_list_so_so(*args): |
206 |
"""
|
207 |
Make a Sollya list out of Sollya objects passed as arguments.
|
208 |
If one wants to call it with a list argument, should prepend a "*"
|
209 |
before the list variable name.
|
210 |
Elements of the list are "eaten" (should not be cleared individually,
|
211 |
are cleared when the list is cleared).
|
212 |
"""
|
213 |
if len(args) == 0: |
214 |
## Called with an empty list produced "error".
|
215 |
return sollya_lib_build_list(None) |
216 |
index = 0
|
217 |
## Append the Sollya elements one by one.
|
218 |
for elementSo in args: |
219 |
if index == 0: |
220 |
listSo = sollya_lib_build_list(elementSo, None)
|
221 |
else:
|
222 |
listSo = sollya_lib_append(listSo, elementSo) |
223 |
index += 1
|
224 |
return listSo
|
225 |
# End pobyso_build list_so_so
|
226 |
|
227 |
|
228 |
def pobyso_change_var_in_function_so_so(funcSo, chvarExpSo): |
229 |
"""
|
230 |
Variable change in a function.
|
231 |
"""
|
232 |
return(sollya_lib_evaluate(funcSo,chvarExpSo))
|
233 |
# End pobyso_change_var_in_function_so_so
|
234 |
|
235 |
def pobyso_chebyshevform_so_so(functionSo, degreeSo, intervalSo): |
236 |
resultSo = sollya_lib_chebyshevform(functionSo, degreeSo, intervalSo) |
237 |
return(resultSo)
|
238 |
# End pobyso_chebyshevform_so_so.
|
239 |
|
240 |
def pobyso_clear_full_list_elements_sa_so(objectListSaSo): |
241 |
"""
|
242 |
Clear the elements of list created by the
|
243 |
pobyso_get_list_elements_so_so function.
|
244 |
objectListSaSo is as follows:
|
245 |
- objectListSaSo[0]: a list of Sollya objects;
|
246 |
- objectListSaSo[1]: the number of elements of the previous list;
|
247 |
- objectListSaSo[2]: an integer that if != 0 states that the list is
|
248 |
end-elliptic
|
249 |
The objects to clear are the elements of the objectListSaSo[0] list.
|
250 |
"""
|
251 |
for index in xrange(0, objectListSaSo[1]): |
252 |
sollya_lib_clear_obj(objectListSaSo[0][index])
|
253 |
# End pobyso_clear_full_list_elements_sa_so
|
254 |
|
255 |
def pobyso_clear_list_elements_sa_so(objectListSaSo): |
256 |
"""
|
257 |
Clear the elements of list of references to Sollya objects
|
258 |
"""
|
259 |
for index in xrange(0, len(objectListSaSo)): |
260 |
sollya_lib_clear_obj(objectListSaSo[index]) |
261 |
# End pobyso_clear_list_elements_sa_so
|
262 |
|
263 |
def pobyso_clear_obj(objSo): |
264 |
"""
|
265 |
Free a Sollya object's memory.
|
266 |
Very thin wrapper around sollya_lib_clear_obj().
|
267 |
"""
|
268 |
sollya_lib_clear_obj(objSo) |
269 |
# End pobyso_clear_obj.
|
270 |
|
271 |
def pobyso_clear_taylorform_sa_so(taylorFormSaSo): |
272 |
"""
|
273 |
This method is rapper around pobyso_clear_list_elements_sa_so.
|
274 |
It is a legacy method left here since it may be used in existing code
|
275 |
where Taylor forms are used as Sage lists obtained by converting
|
276 |
Sollya Taylor forms (a list made of:
|
277 |
- a polynomial;
|
278 |
- a list of intervals enclosing the errors accumulated when computing
|
279 |
the polynomial coefficients;
|
280 |
- a bound on the approximation error between the polynomial and the
|
281 |
function.)
|
282 |
A Taylor form directly obtained from pobyso_taylorform_so_so is cleared
|
283 |
by sollya_lib_clear_obj.
|
284 |
"""
|
285 |
pobyso_clear_list_elements_sa_so(taylorFormSaSo) |
286 |
# End pobyso_clear_taylorform_sa_so
|
287 |
|
288 |
def pobyso_cmp(rnArgSa, cteSo): |
289 |
"""
|
290 |
Deprecated, use pobyso_cmp_sa_so_sa instead.
|
291 |
"""
|
292 |
print "Deprecated, use pobyso_cmp_sa_so_sa instead." |
293 |
return pobyso_cmp_sa_so_sa(rnArgSa, cteSo)
|
294 |
# End pobyso_cmp
|
295 |
|
296 |
def pobyso_cmp_sa_so_sa(rnArgSa, cteSo): |
297 |
"""
|
298 |
Compare the MPFR value a RealNumber with that of a Sollya constant.
|
299 |
|
300 |
Get the value of the Sollya constant into a RealNumber and compare
|
301 |
using MPFR. Could be optimized by working directly with a mpfr_t
|
302 |
for the intermediate number.
|
303 |
"""
|
304 |
# Get the precision of the Sollya constant to build a Sage RealNumber
|
305 |
# with enough precision.to hold it.
|
306 |
precisionOfCte = c_int(0)
|
307 |
# From the Sollya constant, create a local Sage RealNumber.
|
308 |
sollya_lib_get_prec_of_constant(precisionOfCte, cteSo) |
309 |
#print "Precision of constant: ", precisionOfCte
|
310 |
RRRR = RealField(precisionOfCte.value) |
311 |
rnLocalSa = RRRR(0)
|
312 |
sollya_lib_get_constant(get_rn_value(rnLocalSa), cteSo) |
313 |
#
|
314 |
## Compare the Sage RealNumber version of the Sollya constant with rnArg
|
315 |
# through a direct comparison of underlying MPFR numbers.
|
316 |
return cmp_rn_value(rnArgSa, rnLocal)
|
317 |
# End pobyso_smp_sa_so_sa
|
318 |
|
319 |
def pobyso_compute_pos_function_abs_val_bounds_sa_sa(funcSa, lowerBoundSa, \ |
320 |
upperBoundSa): |
321 |
"""
|
322 |
TODO: completely rework and test.
|
323 |
"""
|
324 |
pobyso = pobyso_name_free_variable_sa_so(funcSa.variables()[0])
|
325 |
funcSo = pobyso_parse_string(funcSa._assume_str().replace('_SAGE_VAR_', '')) |
326 |
rangeSo = pobyso_range_sa_so(lowerBoundSa, upperBoundSa) |
327 |
infnormSo = pobyso_infnorm_so_so(funcSo,rangeSo) |
328 |
# Sollya return the infnorm as an interval.
|
329 |
fMaxSa = pobyso_get_interval_from_range_so_sa(infnormSo) |
330 |
# Get the top bound and compute the binade top limit.
|
331 |
fMaxUpperBoundSa = fMaxSa.upper() |
332 |
binadeTopLimitSa = 2**ceil(fMaxUpperBoundSa.log2())
|
333 |
# Put up together the function to use to compute the lower bound.
|
334 |
funcAuxSo = pobyso_parse_string(str(binadeTopLimitSa) + \
|
335 |
'-(' + f._assume_str().replace('_SAGE_VAR_', '') + ')') |
336 |
pobyso_autoprint(funcAuxSo) |
337 |
# Clear the Sollya range before a new call to infnorm and issue the call.
|
338 |
sollya_lib_clear_obj(infnormSo) |
339 |
infnormSo = pobyso_infnorm_so_so(funcAuxSo,rangeSo) |
340 |
fMinSa = pobyso_get_interval_from_range_so_sa(infnormSo) |
341 |
sollya_lib_clear_obj(infnormSo) |
342 |
fMinLowerBoundSa = binadeTopLimitSa - fMinSa.lower() |
343 |
# Compute the maximum of the precisions of the different bounds.
|
344 |
maxPrecSa = max([fMinLowerBoundSa.parent().precision(), \
|
345 |
fMaxUpperBoundSa.parent().precision()]) |
346 |
# Create a RealIntervalField and create an interval with the "good" bounds.
|
347 |
RRRI = RealIntervalField(maxPrecSa) |
348 |
imageIntervalSa = RRRI(fMinLowerBoundSa, fMaxUpperBoundSa) |
349 |
# Free the unneeded Sollya objects
|
350 |
sollya_lib_clear_obj(funcSo) |
351 |
sollya_lib_clear_obj(funcAuxSo) |
352 |
sollya_lib_clear_obj(rangeSo) |
353 |
return(imageIntervalSa)
|
354 |
# End pobyso_compute_pos_function_abs_val_bounds_sa_sa
|
355 |
|
356 |
def pobyso_compute_precision_decay_ratio_function_sa_so(): |
357 |
"""
|
358 |
Compute the precision decay ratio function for polynomial
|
359 |
coefficient progressive trucation.
|
360 |
"""
|
361 |
functionText = """
|
362 |
proc(deg, a, b, we, wq)
|
363 |
{
|
364 |
k = we * (exp(x/a)-1) + wq * (b*x)^2 + (1-we-wq) * x;
|
365 |
return k/k(d);
|
366 |
};
|
367 |
"""
|
368 |
return pobyso_parse_string_sa_so(functionText)
|
369 |
# End pobyso_compute_precision_decay_ratio_function.
|
370 |
|
371 |
|
372 |
def pobyso_constant(rnArg): |
373 |
""" Legacy function. See pobyso_constant_sa_so. """
|
374 |
return(pobyso_constant_sa_so(rnArg))
|
375 |
|
376 |
def pobyso_constant_sa_so(rnArgSa, precisionSa=None): |
377 |
"""
|
378 |
Create a Sollya constant from a Sage RealNumber.
|
379 |
The sollya_lib_constant() function creates a constant
|
380 |
with the same precision as the source.
|
381 |
"""
|
382 |
## Precision stuff. If one wants to change precisions,
|
383 |
# everything takes place in Sage. This only makes
|
384 |
# sense if one wants to reduce the precision.
|
385 |
# TODO: revisit precision stuff with new technique.
|
386 |
# Check that rnArgSa is a realFiedl element.
|
387 |
try:
|
388 |
rnArgSa.ulp() |
389 |
except AttributeError: |
390 |
return pobyso_error_so()
|
391 |
# If a different precision is wanted modify it.
|
392 |
if not precisionSa is None: |
393 |
RRR = RealField(precisionSa) |
394 |
rnArgSa = RRR(rnArgSa) |
395 |
#print rnArgSa, rnArgSa.precision()
|
396 |
# Sollya constant creation takes place here.
|
397 |
return sollya_lib_constant(get_rn_value(rnArgSa))
|
398 |
# End pobyso_constant_sa_so
|
399 |
|
400 |
def pobyso_constant_0_sa_so(): |
401 |
"""
|
402 |
Obvious.
|
403 |
"""
|
404 |
return pobyso_constant_from_int_sa_so(0) |
405 |
|
406 |
def pobyso_constant_1(): |
407 |
"""
|
408 |
Obvious.
|
409 |
Legacy function. See pobyso_constant_so_so.
|
410 |
"""
|
411 |
return pobyso_constant_1_sa_so()
|
412 |
|
413 |
def pobyso_constant_1_sa_so(): |
414 |
"""
|
415 |
Obvious.
|
416 |
"""
|
417 |
return(pobyso_constant_from_int_sa_so(1)) |
418 |
|
419 |
def pobyso_constant_from_int(anInt): |
420 |
""" Legacy function. See pobyso_constant_from_int_sa_so. """
|
421 |
return pobyso_constant_from_int_sa_so(anInt)
|
422 |
|
423 |
def pobyso_constant_from_int_sa_so(anInt): |
424 |
"""
|
425 |
Get a Sollya constant from a Sage int.
|
426 |
"""
|
427 |
return sollya_lib_constant_from_int64(long(anInt)) |
428 |
|
429 |
def pobyso_constant_from_int_so_sa(constSo): |
430 |
"""
|
431 |
Get a Sage int from a Sollya int constant.
|
432 |
Usefull for precision or powers in polynomials.
|
433 |
"""
|
434 |
constSa = c_long(0)
|
435 |
sollya_lib_get_constant_as_int64(byref(constSa), constSo) |
436 |
return constSa.value
|
437 |
# End pobyso_constant_from_int_so_sa
|
438 |
|
439 |
def pobyso_constant_from_mpq_sa_so(rationalSa): |
440 |
"""
|
441 |
Make a Sollya constant from Sage rational.
|
442 |
The Sollya constant is an unevaluated expression.
|
443 |
Hence no precision argument is needed.
|
444 |
It is better to leave this way since Sollya has its own
|
445 |
optimized evaluation mecanism that tries very hard to
|
446 |
return exact values or at least faithful ones.
|
447 |
"""
|
448 |
ratExprSo = \ |
449 |
sollya_lib_constant_from_mpq(sgmp_get_rational_value(rationalSa)) |
450 |
return ratExprSo
|
451 |
# End pobyso_constant_from_mpq_sa_so.
|
452 |
|
453 |
def pobyso_constant_sollya_prec_sa_so(rnArgSa): |
454 |
"""
|
455 |
Create a Sollya constant from a Sage RealNumber at the
|
456 |
current precision in Sollya.
|
457 |
"""
|
458 |
currentSollyaPrecSa = pobyso_get_prec_so_sa() |
459 |
return pobyso_constant_sa_so(rnArgSa, currentSollyaPrecSa)
|
460 |
# End pobyso_constant_sollya_prec_sa_so
|
461 |
|
462 |
def pobyso_end_elliptic_list_so_sa_so(objectsListSo, intCountSa): |
463 |
"""
|
464 |
Create a Sollya end elliptic list made of the objectListSo[0] to
|
465 |
objectsListSo[intCountSa-1] objects.
|
466 |
"""
|
467 |
return sollya_lib_end_elliptic_list(objectSo, int(intCountSa)) |
468 |
|
469 |
def pobyso_error_so(): |
470 |
return sollya_lib_error(None) |
471 |
# End pobyso_error().
|
472 |
|
473 |
def pobyso_evaluate_so_sa(funcSo, argumentSo): |
474 |
"""
|
475 |
Evaluates funcSo for arguemntSo through sollya_lib_evaluate() and return
|
476 |
the result as a Sage object
|
477 |
"""
|
478 |
evalSo = sollya_lib_evaluate(funcSo, argumentSo) |
479 |
if pobyso_is_error_so_sa(evalSo):
|
480 |
return None |
481 |
if pobyso_is_range_so_sa(evalSo):
|
482 |
retVal = pobyso_range_to_interval_so_sa(evalSo) |
483 |
else:
|
484 |
retVal = pobyso_get_constant_as_rn(evalSo) |
485 |
sollya_lib_clear_obj(evalSo) |
486 |
return retVal
|
487 |
# End pobyso_evaluate_so_sa.
|
488 |
|
489 |
def pobyso_evaluate_so_so(funcSo, argumentSo): |
490 |
"""
|
491 |
Evaluates funcSo for arguemntSo through sollya_lib_evaluate().
|
492 |
"""
|
493 |
return sollya_lib_evaluate(funcSo, argumentSo)
|
494 |
# End pobyso_evaluate_so_so.
|
495 |
#
|
496 |
def pobyso_diff_so_so(funcSo): |
497 |
"""
|
498 |
Very thin wrapper around sollya_lib_diff.
|
499 |
"""
|
500 |
## TODO: add a check to make sure funcSo is a functional expression.
|
501 |
return sollya_lib_diff(funcSo)
|
502 |
|
503 |
def pobyso_dirty_find_zeros_so_so(funcSo, rangeSo): |
504 |
"""
|
505 |
Thin wrapper over sollya_lib_dirtyfindzeros()
|
506 |
"""
|
507 |
return sollya_lib_dirtyfindzeros(funcSo, rangeSo)
|
508 |
# End pobys_dirty_find_zeros
|
509 |
|
510 |
def pobyso_dirty_inf_norm_so_so(funcSo, rangeSo, preSo=None): |
511 |
"""
|
512 |
Thin wrapper around sollya_dirtyinfnorm().
|
513 |
"""
|
514 |
# TODO: manage the precision change.
|
515 |
|
516 |
return sollya_lib_dirtyinfnorm(funcSo, rangeSo)
|
517 |
# End pobyso_dirty_inf_norm_so_so
|
518 |
|
519 |
def pobyso_find_zeros_so_so(funcSo, rangeSo): |
520 |
"""
|
521 |
Thin wrapper over sollya_lib_findzeros()
|
522 |
"""
|
523 |
return sollya_lib_findzeros(funcSo, rangeSo)
|
524 |
# End pobys_find_zeros
|
525 |
|
526 |
def pobyso_float_list_so_sa(listSo): |
527 |
"""
|
528 |
Return a Sollya list of floating-point numbers as a Sage list of
|
529 |
floating-point numbers.
|
530 |
TODO: add a test to make sure that each element of the list is a constant.
|
531 |
"""
|
532 |
listSa = [] |
533 |
## The function returns none if the list is empty or an error has happened.
|
534 |
retVal = pobyso_get_list_elements_so_so(listSo) |
535 |
if retVal is None: |
536 |
return listSa
|
537 |
## Just in case the interface is changed and an empty list is returned
|
538 |
# instead of None.
|
539 |
elif len(retVal) == 0: |
540 |
return listSa
|
541 |
else:
|
542 |
## Remember pobyso_get_list_elements_so_so returns more information
|
543 |
# than just the elements of the list (# elements, is_elliptic)
|
544 |
listSaSo, numElements, isEndElliptic = retVal |
545 |
## Return an empty list.
|
546 |
if numElements == 0: |
547 |
return listSa
|
548 |
## Search first for the maximum precision of the elements
|
549 |
maxPrecSa = 0
|
550 |
for floatSo in listSaSo: |
551 |
#pobyso_autoprint(floatSo)
|
552 |
curPrecSa = pobyso_get_prec_of_constant_so_sa(floatSo) |
553 |
if curPrecSa > maxPrecSa:
|
554 |
maxPrecSa = curPrecSa |
555 |
##
|
556 |
RF = RealField(maxPrecSa) |
557 |
##
|
558 |
for floatSo in listSaSo: |
559 |
listSa.append(pobyso_get_constant_as_rn_with_rf_so_sa(floatSo)) |
560 |
return listSa
|
561 |
# End pobyso_float_list_so_sa
|
562 |
|
563 |
def pobyso_float_poly_sa_so(polySa, precSa = None): |
564 |
"""
|
565 |
Create a Sollya polynomial from a Sage RealField polynomial.
|
566 |
"""
|
567 |
## TODO: filter arguments.
|
568 |
## Precision. If a precision is given, convert the polynomial
|
569 |
# into the right polynomial field. If not convert it straight
|
570 |
# to Sollya.
|
571 |
sollyaPrecChanged = False
|
572 |
(initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa() |
573 |
if precSa is None: |
574 |
precSa = polySa.parent().base_ring().precision() |
575 |
if (precSa > initialSollyaPrecSa):
|
576 |
assert precSa >= 2, "Precision change <2 requested" |
577 |
if precSa <= 2: |
578 |
print inspect.stack()[0][3], ": precision change <= 2 requested" |
579 |
precSo = pobyso_constant_from_int(precSa) |
580 |
pobyso_set_prec_so_so(precSo) |
581 |
sollya_lib_clear_obj(precSo) |
582 |
sollyaPrecChanged = True
|
583 |
## Free variable stuff.
|
584 |
freeVariableNameChanged = False
|
585 |
polyFreeVariableNameSa = \ |
586 |
str(polySa.variables()[0]) |
587 |
currentFreeVariableNameSa = pobyso_get_free_variable_name_so_sa() |
588 |
if polyFreeVariableNameSa != currentFreeVariableNameSa:
|
589 |
#print "Free variable names do not match.", polyFreeVariableNameSa
|
590 |
sollya_lib_name_free_variable(polyFreeVariableNameSa) |
591 |
freeVariableNameChanged = True
|
592 |
## Get exponents and coefficients.
|
593 |
exponentsSa = polySa.exponents() |
594 |
coefficientsSa = polySa.coefficients() |
595 |
## Build the polynomial.
|
596 |
polySo = None
|
597 |
for coefficientSa, exponentSa in zip(coefficientsSa, exponentsSa): |
598 |
#print coefficientSa.n(prec=precSa), exponentSa
|
599 |
coefficientSo = \ |
600 |
pobyso_constant_sa_so(coefficientSa) |
601 |
#pobyso_autoprint(coefficientSo)
|
602 |
exponentSo = \ |
603 |
pobyso_constant_from_int_sa_so(exponentSa) |
604 |
#pobyso_autoprint(exponentSo)
|
605 |
monomialSo = sollya_lib_build_function_pow( |
606 |
sollya_lib_build_function_free_variable(), |
607 |
exponentSo) |
608 |
polyTermSo = sollya_lib_build_function_mul(coefficientSo, |
609 |
monomialSo) |
610 |
if polySo is None: |
611 |
polySo = polyTermSo |
612 |
else:
|
613 |
polySo = sollya_lib_build_function_add(polySo, polyTermSo) |
614 |
if sollyaPrecChanged:
|
615 |
pobyso_set_prec_so_so(initialSollyaPrecSo) |
616 |
sollya_lib_clear_obj(initialSollyaPrecSo) |
617 |
## Do not set back the free variable name in Sollya to its initial value:
|
618 |
# it will change it back, in the Sollya polynomial, to what it was in the
|
619 |
# first place.
|
620 |
return polySo
|
621 |
# End pobyso_float_poly_sa_so
|
622 |
|
623 |
def pobyso_float_poly_so_sa(polySo, realFieldSa=None): |
624 |
"""
|
625 |
Convert a Sollya polynomial into a Sage floating-point polynomial.
|
626 |
If no realField is given, a RealField corresponding to the maximum
|
627 |
precision of the coefficients is internally computed.
|
628 |
The real field is not returned but can be easily retrieved from
|
629 |
the polynomial itself.
|
630 |
ALGORITHM:
|
631 |
- (optional) compute the RealField of the coefficients;
|
632 |
- convert the Sollya expression into a Sage expression;
|
633 |
- convert the Sage expression into a Sage polynomial
|
634 |
"""
|
635 |
if realFieldSa is None: |
636 |
expressionPrecSa = pobyso_get_max_prec_of_exp_so_sa(polySo) |
637 |
#print "Maximum precision of Sollya polynomial coefficients:", expressionPrecSa
|
638 |
if expressionPrecSa < 2 or expressionPrecSa > 2147483391: |
639 |
print "Maximum degree of expression:", expressionPrecSa |
640 |
realFieldSa = RealField(expressionPrecSa) |
641 |
#print "Sollya expression before...",
|
642 |
#pobyso_autoprint(polySo)
|
643 |
|
644 |
expressionSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, |
645 |
realFieldSa) |
646 |
#print "...Sollya expression after."
|
647 |
#pobyso_autoprint(polySo)
|
648 |
polyVariableSa = expressionSa.variables()[0]
|
649 |
polyRingSa = realFieldSa[str(polyVariableSa)]
|
650 |
#print polyRingSa
|
651 |
# Do not use the polynomial(expressionSa, ring=polyRingSa) form!
|
652 |
polynomialSa = polyRingSa(expressionSa) |
653 |
polyCoeffsListSa = polynomialSa.coefficients() |
654 |
#for coeff in polyCoeffsListSa:
|
655 |
# print coeff.abs().n()
|
656 |
return polynomialSa
|
657 |
# End pobyso_float_poly_so_sa
|
658 |
|
659 |
def pobyso_free_variable(): |
660 |
"""
|
661 |
Ultra thin wrapper around the sollya_lib_function_build_free_variable function.
|
662 |
"""
|
663 |
return sollya_lib_build_function_free_variable()
|
664 |
|
665 |
def pobyso_function_type_as_string(funcType): |
666 |
""" Legacy function. See pobyso_function_type_as_string_so_sa. """
|
667 |
return(pobyso_function_type_as_string_so_sa(funcType))
|
668 |
|
669 |
def pobyso_function_type_as_string_so_sa(funcType): |
670 |
"""
|
671 |
Numeric Sollya function codes -> Sage mathematical function names.
|
672 |
Notice that pow -> ^ (a la Sage, not a la Python).
|
673 |
"""
|
674 |
if funcType == SOLLYA_BASE_FUNC_ABS:
|
675 |
return "abs" |
676 |
elif funcType == SOLLYA_BASE_FUNC_ACOS:
|
677 |
return "arccos" |
678 |
elif funcType == SOLLYA_BASE_FUNC_ACOSH:
|
679 |
return "arccosh" |
680 |
elif funcType == SOLLYA_BASE_FUNC_ADD:
|
681 |
return "+" |
682 |
elif funcType == SOLLYA_BASE_FUNC_ASIN:
|
683 |
return "arcsin" |
684 |
elif funcType == SOLLYA_BASE_FUNC_ASINH:
|
685 |
return "arcsinh" |
686 |
elif funcType == SOLLYA_BASE_FUNC_ATAN:
|
687 |
return "arctan" |
688 |
elif funcType == SOLLYA_BASE_FUNC_ATANH:
|
689 |
return "arctanh" |
690 |
elif funcType == SOLLYA_BASE_FUNC_CEIL:
|
691 |
return "ceil" |
692 |
elif funcType == SOLLYA_BASE_FUNC_CONSTANT:
|
693 |
return "cte" |
694 |
elif funcType == SOLLYA_BASE_FUNC_COS:
|
695 |
return "cos" |
696 |
elif funcType == SOLLYA_BASE_FUNC_COSH:
|
697 |
return "cosh" |
698 |
elif funcType == SOLLYA_BASE_FUNC_DIV:
|
699 |
return "/" |
700 |
elif funcType == SOLLYA_BASE_FUNC_DOUBLE:
|
701 |
return "double" |
702 |
elif funcType == SOLLYA_BASE_FUNC_DOUBLEDOUBLE:
|
703 |
return "doubleDouble" |
704 |
elif funcType == SOLLYA_BASE_FUNC_DOUBLEEXTENDED:
|
705 |
return "doubleDxtended" |
706 |
elif funcType == SOLLYA_BASE_FUNC_ERF:
|
707 |
return "erf" |
708 |
elif funcType == SOLLYA_BASE_FUNC_ERFC:
|
709 |
return "erfc" |
710 |
elif funcType == SOLLYA_BASE_FUNC_EXP:
|
711 |
return "exp" |
712 |
elif funcType == SOLLYA_BASE_FUNC_EXP_M1:
|
713 |
return "expm1" |
714 |
elif funcType == SOLLYA_BASE_FUNC_FLOOR:
|
715 |
return "floor" |
716 |
elif funcType == SOLLYA_BASE_FUNC_FREE_VARIABLE:
|
717 |
return "freeVariable" |
718 |
elif funcType == SOLLYA_BASE_FUNC_HALFPRECISION:
|
719 |
return "halfPrecision" |
720 |
elif funcType == SOLLYA_BASE_FUNC_LIBRARYCONSTANT:
|
721 |
return "libraryConstant" |
722 |
elif funcType == SOLLYA_BASE_FUNC_LIBRARYFUNCTION:
|
723 |
return "libraryFunction" |
724 |
elif funcType == SOLLYA_BASE_FUNC_LOG:
|
725 |
return "log" |
726 |
elif funcType == SOLLYA_BASE_FUNC_LOG_10:
|
727 |
return "log10" |
728 |
elif funcType == SOLLYA_BASE_FUNC_LOG_1P:
|
729 |
return "log1p" |
730 |
elif funcType == SOLLYA_BASE_FUNC_LOG_2:
|
731 |
return "log2" |
732 |
elif funcType == SOLLYA_BASE_FUNC_MUL:
|
733 |
return "*" |
734 |
elif funcType == SOLLYA_BASE_FUNC_NEARESTINT:
|
735 |
return "round" |
736 |
elif funcType == SOLLYA_BASE_FUNC_NEG:
|
737 |
return "__neg__" |
738 |
elif funcType == SOLLYA_BASE_FUNC_PI:
|
739 |
return "pi" |
740 |
elif funcType == SOLLYA_BASE_FUNC_POW:
|
741 |
return "^" |
742 |
elif funcType == SOLLYA_BASE_FUNC_PROCEDUREFUNCTION:
|
743 |
return "procedureFunction" |
744 |
elif funcType == SOLLYA_BASE_FUNC_QUAD:
|
745 |
return "quad" |
746 |
elif funcType == SOLLYA_BASE_FUNC_SIN:
|
747 |
return "sin" |
748 |
elif funcType == SOLLYA_BASE_FUNC_SINGLE:
|
749 |
return "single" |
750 |
elif funcType == SOLLYA_BASE_FUNC_SINH:
|
751 |
return "sinh" |
752 |
elif funcType == SOLLYA_BASE_FUNC_SQRT:
|
753 |
return "sqrt" |
754 |
elif funcType == SOLLYA_BASE_FUNC_SUB:
|
755 |
return "-" |
756 |
elif funcType == SOLLYA_BASE_FUNC_TAN:
|
757 |
return "tan" |
758 |
elif funcType == SOLLYA_BASE_FUNC_TANH:
|
759 |
return "tanh" |
760 |
elif funcType == SOLLYA_BASE_FUNC_TRIPLEDOUBLE:
|
761 |
return "tripleDouble" |
762 |
else:
|
763 |
return None |
764 |
|
765 |
def pobyso_get_constant(rnArgSa, constSo): |
766 |
""" Legacy function. See pobyso_get_constant_so_sa. """
|
767 |
return pobyso_get_constant_so_sa(rnArgSa, constSo)
|
768 |
# End pobyso_get_constant
|
769 |
|
770 |
def pobyso_get_constant_so_sa(rnArgSa, constSo): |
771 |
"""
|
772 |
Set the value of rnArgSa to the value of constSo in MPFR_RNDN mode.
|
773 |
rnArg must already exist and belong to some RealField.
|
774 |
We assume that constSo points to a Sollya constant.
|
775 |
"""
|
776 |
outcome = sollya_lib_get_constant(get_rn_value(rnArgSa), constSo) |
777 |
if outcome == 0: # Failure because constSo is not a constant expression. |
778 |
return None |
779 |
else:
|
780 |
return outcome
|
781 |
# End pobyso_get_constant_so_sa
|
782 |
|
783 |
def pobyso_get_constant_as_rn(ctExpSo): |
784 |
"""
|
785 |
Legacy function. See pobyso_get_constant_as_rn_so_sa.
|
786 |
"""
|
787 |
return(pobyso_get_constant_as_rn_so_sa(ctExpSo))
|
788 |
|
789 |
def pobyso_get_constant_as_rn_so_sa(constExpSo): |
790 |
"""
|
791 |
Get a Sollya constant as a Sage "real number".
|
792 |
The precision of the floating-point number returned is that of the Sollya
|
793 |
constant.
|
794 |
"""
|
795 |
#print "Before computing precision of variable..."
|
796 |
#pobyso_autoprint(constExpSo)
|
797 |
precisionSa = pobyso_get_prec_of_constant_so_sa(constExpSo) |
798 |
#print "precisionSa:", precisionSa
|
799 |
## If the expression can not be exactly converted, None is returned.
|
800 |
# In this case opt for the Sollya current expression.
|
801 |
if precisionSa is None: |
802 |
precisionSa = pobyso_get_prec_so_sa() |
803 |
RRRR = RealField(precisionSa) |
804 |
rnSa = RRRR(0)
|
805 |
outcome = sollya_lib_get_constant(get_rn_value(rnSa), constExpSo) |
806 |
if outcome == 0: |
807 |
return None |
808 |
else:
|
809 |
return rnSa
|
810 |
# End pobyso_get_constant_as_rn_so_sa
|
811 |
|
812 |
def pobyso_get_constant_as_rn_with_rf(ctExp, realField): |
813 |
"""
|
814 |
Legacy function. See pobyso_get_constant_as_rn_with_rf_so_sa.
|
815 |
"""
|
816 |
return pobyso_get_constant_as_rn_with_rf_so_sa(ctExp, realField)
|
817 |
# End pobyso_get_constant_as_rn_with_rf
|
818 |
|
819 |
def pobyso_get_constant_as_rn_with_rf_so_sa(ctExpSo, realFieldSa = None): |
820 |
"""
|
821 |
Get a Sollya constant as a Sage "real number".
|
822 |
If no real field is specified, the precision of the floating-point number
|
823 |
returned is that of the Sollya constant.
|
824 |
Otherwise is is that of the real field. Hence rounding may happen.
|
825 |
"""
|
826 |
if realFieldSa is None: |
827 |
return pobyso_get_constant_as_rn_so_sa(ctExpSo)
|
828 |
rnSa = realFieldSa(0)
|
829 |
outcome = sollya_lib_get_constant(get_rn_value(rnSa), ctExpSo) |
830 |
if outcome == 0: |
831 |
return None |
832 |
else:
|
833 |
return rnSa
|
834 |
# End pobyso_get_constant_as_rn_with_rf_so_sa
|
835 |
|
836 |
def pobyso_get_free_variable_name(): |
837 |
"""
|
838 |
Legacy function. See pobyso_get_free_variable_name_so_sa.
|
839 |
"""
|
840 |
return(pobyso_get_free_variable_name_so_sa())
|
841 |
|
842 |
def pobyso_get_free_variable_name_so_sa(): |
843 |
return sollya_lib_get_free_variable_name()
|
844 |
|
845 |
def pobyso_get_function_arity(expressionSo): |
846 |
"""
|
847 |
Legacy function. See pobyso_get_function_arity_so_sa.
|
848 |
"""
|
849 |
return(pobyso_get_function_arity_so_sa(expressionSo))
|
850 |
|
851 |
def pobyso_get_function_arity_so_sa(expressionSo): |
852 |
arity = c_int(0)
|
853 |
sollya_lib_get_function_arity(byref(arity),expressionSo) |
854 |
return int(arity.value) |
855 |
|
856 |
def pobyso_get_head_function(expressionSo): |
857 |
"""
|
858 |
Legacy function. See pobyso_get_head_function_so_sa.
|
859 |
"""
|
860 |
return(pobyso_get_head_function_so_sa(expressionSo))
|
861 |
|
862 |
def pobyso_get_head_function_so_sa(expressionSo): |
863 |
functionType = c_int(0)
|
864 |
sollya_lib_get_head_function(byref(functionType), expressionSo) |
865 |
return int(functionType.value) |
866 |
|
867 |
def pobyso_get_interval_from_range_so_sa(soRange, realIntervalFieldSa = None ): |
868 |
"""
|
869 |
Return the Sage interval corresponding to the Sollya range argument.
|
870 |
If no reaIntervalField is passed as an argument, the interval bounds are not
|
871 |
rounded: they are elements of RealIntervalField of the "right" precision
|
872 |
to hold all the digits.
|
873 |
"""
|
874 |
prec = c_int(0)
|
875 |
if realIntervalFieldSa is None: |
876 |
retval = sollya_lib_get_prec_of_range(byref(prec), soRange, None)
|
877 |
if retval == 0: |
878 |
return None |
879 |
realIntervalFieldSa = RealIntervalField(prec.value) |
880 |
intervalSa = realIntervalFieldSa(0,0) |
881 |
retval = \ |
882 |
sollya_lib_get_interval_from_range(get_interval_value(intervalSa),\ |
883 |
soRange) |
884 |
if retval == 0: |
885 |
return None |
886 |
return intervalSa
|
887 |
# End pobyso_get_interval_from_range_so_sa
|
888 |
|
889 |
def pobyso_get_list_elements(soObj): |
890 |
""" Legacy function. See pobyso_get_list_elements_so_so. """
|
891 |
return pobyso_get_list_elements_so_so(soObj)
|
892 |
|
893 |
def pobyso_get_list_elements_so_so(objectListSo): |
894 |
"""
|
895 |
Get the Sollya list elements as a Sage/Python array of Sollya objects.
|
896 |
|
897 |
INPUT:
|
898 |
- objectListSo: a Sollya list of Sollya objects.
|
899 |
|
900 |
OUTPUT:
|
901 |
- a Sage/Python tuple made of:
|
902 |
- a Sage/Python list of Sollya objects,
|
903 |
- a Sage/Python int holding the number of elements,
|
904 |
- a Sage/Python int stating (!= 0) that the list is end-elliptic.
|
905 |
NOTE::
|
906 |
We recover the addresses of the Sollya object from the list of pointers
|
907 |
returned by sollya_lib_get_list_elements. The list itself is freed.
|
908 |
TODO::
|
909 |
Figure out what to do with numElements since the number of elements
|
910 |
can easily be recovered from the list itself.
|
911 |
Ditto for isEndElliptic.
|
912 |
"""
|
913 |
listAddress = POINTER(c_longlong)() |
914 |
numElements = c_int(0)
|
915 |
isEndElliptic = c_int(0)
|
916 |
listAsSageList = [] |
917 |
result = sollya_lib_get_list_elements(byref(listAddress),\ |
918 |
byref(numElements),\ |
919 |
byref(isEndElliptic),\ |
920 |
objectListSo) |
921 |
if result == 0 : |
922 |
return None |
923 |
for i in xrange(0, numElements.value, 1): |
924 |
#listAsSageList.append(sollya_lib_copy_obj(listAddress[i]))
|
925 |
listAsSageList.append(listAddress[i]) |
926 |
# Clear each of the elements returned by Sollya.
|
927 |
#sollya_lib_clear_obj(listAddress[i])
|
928 |
# Free the list itself.
|
929 |
sollya_lib_free(listAddress) |
930 |
return (listAsSageList, numElements.value, isEndElliptic.value)
|
931 |
|
932 |
def pobyso_get_max_prec_of_exp(soExp): |
933 |
""" Legacy function. See pobyso_get_max_prec_of_exp_so_sa. """
|
934 |
return pobyso_get_max_prec_of_exp_so_sa(soExp)
|
935 |
|
936 |
def pobyso_get_max_prec_of_exp_so_sa(expSo): |
937 |
"""
|
938 |
Get the maximum precision used for the numbers in a Sollya expression.
|
939 |
|
940 |
Arguments:
|
941 |
soExp -- a Sollya expression pointer
|
942 |
Return value:
|
943 |
A Python integer
|
944 |
TODO:
|
945 |
- error management;
|
946 |
- correctly deal with numerical type such as DOUBLEEXTENDED.
|
947 |
"""
|
948 |
if expSo is None: |
949 |
print inspect.stack()[0][3], ": expSo is None." |
950 |
return 0 |
951 |
maxPrecision = 0
|
952 |
minConstPrec = 0
|
953 |
currentConstPrec = 0
|
954 |
#pobyso_autoprint(expSo)
|
955 |
operator = pobyso_get_head_function_so_sa(expSo) |
956 |
if (operator != SOLLYA_BASE_FUNC_CONSTANT) and \ |
957 |
(operator != SOLLYA_BASE_FUNC_FREE_VARIABLE): |
958 |
(arity, subexpressions) = pobyso_get_subfunctions_so_sa(expSo) |
959 |
for i in xrange(arity): |
960 |
maxPrecisionCandidate = \ |
961 |
pobyso_get_max_prec_of_exp_so_sa(subexpressions[i]) |
962 |
if maxPrecisionCandidate > maxPrecision:
|
963 |
maxPrecision = maxPrecisionCandidate |
964 |
return maxPrecision
|
965 |
elif operator == SOLLYA_BASE_FUNC_CONSTANT:
|
966 |
#minConstPrec = pobyso_get_min_prec_of_constant_so_sa(expSo)
|
967 |
#currentConstPrec = pobyso_get_min_prec_of_constant_so_sa(soExp)
|
968 |
#print minConstPrec, " - ", currentConstPrec
|
969 |
return pobyso_get_min_prec_of_constant_so_sa(expSo)
|
970 |
|
971 |
elif operator == SOLLYA_BASE_FUNC_FREE_VARIABLE:
|
972 |
return 0 |
973 |
else:
|
974 |
print "pobyso_get_max_prec_of_exp_so_sa: unexepected operator." |
975 |
return 0 |
976 |
|
977 |
def pobyso_get_min_prec_of_constant_so_sa(constExpSo): |
978 |
"""
|
979 |
Get the minimum precision necessary to represent the value of a Sollya
|
980 |
constant.
|
981 |
MPFR_MIN_PREC and powers of 2 are taken into account.
|
982 |
We assume that constExpSo is a pointer to a Sollay constant expression.
|
983 |
"""
|
984 |
constExpAsRnSa = pobyso_get_constant_as_rn_so_sa(constExpSo) |
985 |
return(min_mpfr_size(get_rn_value(constExpAsRnSa)))
|
986 |
|
987 |
def pobyso_get_poly_so_sa(polySo, realFieldSa=None): |
988 |
"""
|
989 |
Convert a Sollya polynomial into a Sage polynomial.
|
990 |
Legacy function. Use pobyso_float_poly_so_sa() instead.
|
991 |
"""
|
992 |
return pobyso_float_poly_so_sa(polySo,realFieldSa)
|
993 |
# End pobyso_get_poly_so_sa
|
994 |
|
995 |
def pobyso_get_prec(): |
996 |
""" Legacy function. See pobyso_get_prec_so_sa(). """
|
997 |
return pobyso_get_prec_so_sa()
|
998 |
|
999 |
def pobyso_get_prec_so(): |
1000 |
"""
|
1001 |
Get the current default precision in Sollya.
|
1002 |
The return value is a Sollya object.
|
1003 |
Usefull when modifying the precision back and forth by avoiding
|
1004 |
extra conversions.
|
1005 |
"""
|
1006 |
return sollya_lib_get_prec(None) |
1007 |
|
1008 |
def pobyso_get_prec_so_sa(): |
1009 |
"""
|
1010 |
Get the current default precision in Sollya.
|
1011 |
The return value is Sage/Python int.
|
1012 |
"""
|
1013 |
precSo = sollya_lib_get_prec() |
1014 |
precSa = pobyso_constant_from_int_so_sa(precSo) |
1015 |
sollya_lib_clear_obj(precSo) |
1016 |
return precSa
|
1017 |
# End pobyso_get_prec_so_sa.
|
1018 |
|
1019 |
def pobyso_get_prec_so_so_sa(): |
1020 |
"""
|
1021 |
Return the current precision both as a Sollya object and a
|
1022 |
Sage integer as hybrid tuple.
|
1023 |
To avoid multiple calls for precision manipulations.
|
1024 |
"""
|
1025 |
precSo = sollya_lib_get_prec() |
1026 |
precSa = pobyso_constant_from_int_so_sa(precSo) |
1027 |
return (precSo, int(precSa)) |
1028 |
|
1029 |
def pobyso_get_prec_of_constant(ctExpSo): |
1030 |
""" Legacy function. See pobyso_get_prec_of_constant_so_sa. """
|
1031 |
return pobyso_get_prec_of_constant_so_sa(ctExpSo)
|
1032 |
|
1033 |
def pobyso_get_prec_of_constant_so_sa(ctExpSo): |
1034 |
"""
|
1035 |
Tries to find a precision to represent ctExpSo without rounding.
|
1036 |
If not possible, returns None.
|
1037 |
"""
|
1038 |
#print "Entering pobyso_get_prec_of_constant_so_sa..."
|
1039 |
prec = c_int(0)
|
1040 |
retc = sollya_lib_get_prec_of_constant(byref(prec), ctExpSo, None)
|
1041 |
if retc == 0: |
1042 |
#print "pobyso_get_prec_of_constant_so_sa failed."
|
1043 |
return None |
1044 |
#print "...exiting pobyso_get_prec_of_constant_so_sa."
|
1045 |
return int(prec.value) |
1046 |
|
1047 |
def pobyso_get_prec_of_range_so_sa(rangeSo): |
1048 |
"""
|
1049 |
Returns the number of bits elements of a range are coded with.
|
1050 |
"""
|
1051 |
prec = c_int(0)
|
1052 |
retc = sollya_lib_get_prec_of_range(byref(prec), rangeSo, None)
|
1053 |
if retc == 0: |
1054 |
return(None) |
1055 |
return int(prec.value) |
1056 |
# End pobyso_get_prec_of_range_so_sa()
|
1057 |
|
1058 |
def pobyso_get_sage_exp_from_sollya_exp(sollyaExpSo, realField = RR): |
1059 |
""" Legacy function. See pobyso_get_sage_exp_from_sollya_exp_so_sa. """
|
1060 |
return pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo,
|
1061 |
realField = RR) |
1062 |
|
1063 |
def pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, realFieldSa = RR): |
1064 |
"""
|
1065 |
Get a Sage expression from a Sollya expression.
|
1066 |
Currently only tested with polynomials with floating-point coefficients.
|
1067 |
Notice that, in the returned polynomial, the exponents are RealNumbers.
|
1068 |
"""
|
1069 |
#pobyso_autoprint(sollyaExp)
|
1070 |
operatorSa = pobyso_get_head_function_so_sa(sollyaExpSo) |
1071 |
sollyaLibFreeVariableName = sollya_lib_get_free_variable_name() |
1072 |
## Get rid of the "_"'s in "_x_", if any.
|
1073 |
sollyaLibFreeVariableName = re.sub('_', '', sollyaLibFreeVariableName) |
1074 |
# Constants and the free variable are special cases.
|
1075 |
# All other operator are dealt with in the same way.
|
1076 |
if (operatorSa != SOLLYA_BASE_FUNC_CONSTANT) and \ |
1077 |
(operatorSa != SOLLYA_BASE_FUNC_FREE_VARIABLE): |
1078 |
(aritySa, subexpressionsSa) = pobyso_get_subfunctions_so_sa(sollyaExpSo) |
1079 |
if aritySa == 1: |
1080 |
sageExpSa = eval(pobyso_function_type_as_string_so_sa(operatorSa) + \
|
1081 |
"(" + pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], \ |
1082 |
realFieldSa) + ")")
|
1083 |
elif aritySa == 2: |
1084 |
# We do not get through the preprocessor.
|
1085 |
# The "^" operator is then a special case.
|
1086 |
if operatorSa == SOLLYA_BASE_FUNC_POW:
|
1087 |
operatorAsStringSa = "**"
|
1088 |
else:
|
1089 |
operatorAsStringSa = \ |
1090 |
pobyso_function_type_as_string_so_sa(operatorSa) |
1091 |
sageExpSa = \ |
1092 |
eval("pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], realFieldSa)"\ |
1093 |
+ " " + operatorAsStringSa + " " + \ |
1094 |
"pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[1], realFieldSa)")
|
1095 |
# We do not know yet how to deal with arity >= 3
|
1096 |
# (is there any in Sollya anyway?).
|
1097 |
else:
|
1098 |
sageExpSa = eval('None') |
1099 |
return sageExpSa
|
1100 |
elif operatorSa == SOLLYA_BASE_FUNC_CONSTANT:
|
1101 |
#print "This is a constant"
|
1102 |
return pobyso_get_constant_as_rn_with_rf_so_sa(sollyaExpSo, realFieldSa)
|
1103 |
elif operatorSa == SOLLYA_BASE_FUNC_FREE_VARIABLE:
|
1104 |
#print "This is the free variable"
|
1105 |
return eval(sollyaLibFreeVariableName) |
1106 |
else:
|
1107 |
print "Unexpected" |
1108 |
return eval('None') |
1109 |
# End pobyso_get_sage_exp_from_sollya_exp_so_sa
|
1110 |
|
1111 |
|
1112 |
def pobyso_get_subfunctions(expressionSo): |
1113 |
""" Legacy function. See pobyso_get_subfunctions_so_sa. """
|
1114 |
return pobyso_get_subfunctions_so_sa(expressionSo)
|
1115 |
# End pobyso_get_subfunctions.
|
1116 |
|
1117 |
def pobyso_get_subfunctions_so_sa(expressionSo): |
1118 |
"""
|
1119 |
Get the subfunctions of an expression.
|
1120 |
Return the number of subfunctions and the list of subfunctions addresses.
|
1121 |
S.T.: Could not figure out another way than that ugly list of declarations
|
1122 |
to recover the addresses of the subfunctions.
|
1123 |
We limit ourselves to arity 8 functions.
|
1124 |
"""
|
1125 |
subf0 = c_int(0)
|
1126 |
subf1 = c_int(0)
|
1127 |
subf2 = c_int(0)
|
1128 |
subf3 = c_int(0)
|
1129 |
subf4 = c_int(0)
|
1130 |
subf5 = c_int(0)
|
1131 |
subf6 = c_int(0)
|
1132 |
subf7 = c_int(0)
|
1133 |
subf8 = c_int(0)
|
1134 |
arity = c_int(0)
|
1135 |
nullPtr = POINTER(c_int)() |
1136 |
sollya_lib_get_subfunctions(expressionSo, byref(arity), \ |
1137 |
byref(subf0), byref(subf1), byref(subf2), byref(subf3), \ |
1138 |
byref(subf4), byref(subf5),\ |
1139 |
byref(subf6), byref(subf7), byref(subf8), nullPtr, None)
|
1140 |
# byref(cast(subfunctions[0], POINTER(c_int))), \
|
1141 |
# byref(cast(subfunctions[0], POINTER(c_int))), \
|
1142 |
# byref(cast(subfunctions[2], POINTER(c_int))), \
|
1143 |
# byref(cast(subfunctions[3], POINTER(c_int))), \
|
1144 |
# byref(cast(subfunctions[4], POINTER(c_int))), \
|
1145 |
# byref(cast(subfunctions[5], POINTER(c_int))), \
|
1146 |
# byref(cast(subfunctions[6], POINTER(c_int))), \
|
1147 |
# byref(cast(subfunctions[7], POINTER(c_int))), \
|
1148 |
# byref(cast(subfunctions[8], POINTER(c_int))), nullPtr)
|
1149 |
subfunctions = [subf0, subf1, subf2, subf3, subf4, subf5, subf6, subf7, \ |
1150 |
subf8] |
1151 |
subs = [] |
1152 |
if arity.value > pobyso_max_arity:
|
1153 |
return(0,[]) |
1154 |
for i in xrange(arity.value): |
1155 |
subs.append(int(subfunctions[i].value))
|
1156 |
#print subs[i]
|
1157 |
return (int(arity.value), subs) |
1158 |
# End pobyso_get_subfunctions_so_sa
|
1159 |
|
1160 |
def pobyso_guess_degree_sa_sa(functionSa, intervalSa, approxErrorSa, |
1161 |
weightSa=None, degreeBoundSa=None): |
1162 |
"""
|
1163 |
Sa_sa variant of the solly_guessdegree function.
|
1164 |
Return 0 if something goes wrong.
|
1165 |
"""
|
1166 |
functionAsStringSa = functionSa._assume_str().replace('_SAGE_VAR_', '') |
1167 |
functionSo = pobyso_parse_string_sa_so(functionAsStringSa) |
1168 |
if pobyso_is_error_so_sa(functionSo):
|
1169 |
sollya_lib_clear_obj(functionSo) |
1170 |
return 0 |
1171 |
rangeSo = pobyso_interval_to_range_sa_so(intervalSa) |
1172 |
# The approximation error is expected to be a floating point number.
|
1173 |
if pobyso_is_floating_point_number_sa_sa(approxErrorSa):
|
1174 |
approxErrorSo = pobyso_constant_sa_so(approxErrorSa) |
1175 |
else:
|
1176 |
approxErrorSo = pobyso_constant_sa_so(RR(approxErrorSa)) |
1177 |
if not weightSa is None: |
1178 |
weightAsStringSa = weightSa._assume_str().replace('_SAGE_VAR_', '') |
1179 |
weightSo = pobyso_parse_string_sa_so(weightAsStringSa) |
1180 |
if pobyso_is_error_so_sa(weightSo):
|
1181 |
sollya_lib_clear_obj(functionSo) |
1182 |
sollya_lib_clear_obj(rangeSo) |
1183 |
sollya_lib_clear_obj(approxErrorSo) |
1184 |
sollya_lib_clear_obj(weightSo) |
1185 |
return 0 |
1186 |
else:
|
1187 |
weightSo = None
|
1188 |
if not degreeBoundSa is None: |
1189 |
degreeBoundSo = pobyso_constant_from_int_sa_so(degreeBoundSa) |
1190 |
else:
|
1191 |
degreeBoundSo = None
|
1192 |
guessedDegreeSa = pobyso_guess_degree_so_sa(functionSo, |
1193 |
rangeSo, |
1194 |
approxErrorSo, |
1195 |
weightSo, |
1196 |
degreeBoundSo) |
1197 |
sollya_lib_clear_obj(functionSo) |
1198 |
sollya_lib_clear_obj(rangeSo) |
1199 |
sollya_lib_clear_obj(approxErrorSo) |
1200 |
if not weightSo is None: |
1201 |
sollya_lib_clear_obj(weightSo) |
1202 |
if not degreeBoundSo is None: |
1203 |
sollya_lib_clear_obj(degreeBoundSo) |
1204 |
return guessedDegreeSa
|
1205 |
# End poyso_guess_degree_sa_sa
|
1206 |
|
1207 |
def pobyso_guess_degree_so_sa(functionSo, rangeSo, errorSo, weightSo=None, \ |
1208 |
degreeBoundSo=None):
|
1209 |
"""
|
1210 |
Thin wrapper around the guessdegree function.
|
1211 |
Nevertheless, some precision control stuff has been appended.
|
1212 |
"""
|
1213 |
# Deal with Sollya internal precision issues: if it is too small,
|
1214 |
# compared with the error, increases it to about twice -log2(error).
|
1215 |
errorSa = pobyso_get_constant_as_rn_with_rf_so_sa(errorSo) |
1216 |
log2ErrorSa = errorSa.log2() |
1217 |
if log2ErrorSa < 0: |
1218 |
neededPrecisionSa = int(2 * int(-log2ErrorSa) / 64) * 64 |
1219 |
else:
|
1220 |
neededPrecisionSa = int(2 * int(log2ErrorSa) / 64) * 64 |
1221 |
#print "Needed precision:", neededPrecisionSa
|
1222 |
sollyaPrecisionHasChanged = False
|
1223 |
(initialPrecSo, initialPrecSa) = pobyso_get_prec_so_so_sa() |
1224 |
if neededPrecisionSa > initialPrecSa:
|
1225 |
if neededPrecisionSa <= 2: |
1226 |
print inspect.stack()[0][3], ": precision change <= 2 requested." |
1227 |
pobyso_set_prec_sa_so(neededPrecisionSa) |
1228 |
sollyaPrecisionHasChanged = True
|
1229 |
#print "Guessing degree..."
|
1230 |
# weightSo and degreeBoundsSo are optional arguments.
|
1231 |
# As declared, sollya_lib_guessdegree must take 5 arguments.
|
1232 |
if weightSo is None: |
1233 |
degreeRangeSo = sollya_lib_guessdegree(functionSo, rangeSo, errorSo, |
1234 |
0, 0, None) |
1235 |
elif degreeBoundSo is None: |
1236 |
degreeRangeSo = sollya_lib_guessdegree(functionSo, rangeSo, \ |
1237 |
errorSo, weightSo, 0, None) |
1238 |
else:
|
1239 |
degreeRangeSo = sollya_lib_guessdegree(functionSo, rangeSo, errorSo, \ |
1240 |
weightSo, degreeBoundSo, None)
|
1241 |
#print "...degree guess done."
|
1242 |
# Restore internal precision, if applicable.
|
1243 |
if sollyaPrecisionHasChanged:
|
1244 |
pobyso_set_prec_so_so(initialPrecSo) |
1245 |
sollya_lib_clear_obj(initialPrecSo) |
1246 |
degreeIntervalSa = pobyso_range_to_interval_so_sa(degreeRangeSo) |
1247 |
sollya_lib_clear_obj(degreeRangeSo) |
1248 |
# When ok, both bounds match.
|
1249 |
# When the degree bound is too low, the upper bound is the degree
|
1250 |
# for which the error can be honored.
|
1251 |
# When it really goes wrong, the upper bound is infinity.
|
1252 |
if degreeIntervalSa.lower() == degreeIntervalSa.upper():
|
1253 |
return int(degreeIntervalSa.lower()) |
1254 |
else:
|
1255 |
if degreeIntervalSa.upper().is_infinity():
|
1256 |
return None |
1257 |
else:
|
1258 |
return int(degreeIntervalSa.upper()) |
1259 |
# End pobyso_guess_degree_so_sa
|
1260 |
|
1261 |
def pobyso_inf_so_so(intervalSo): |
1262 |
"""
|
1263 |
Very thin wrapper around sollya_lib_inf().
|
1264 |
"""
|
1265 |
return sollya_lib_inf(intervalSo)
|
1266 |
# End pobyso_inf_so_so.
|
1267 |
#
|
1268 |
def pobyso_infnorm_sa_sa(funcSa, intervalSa): |
1269 |
"""
|
1270 |
An infnorm call with Sage arguments.
|
1271 |
We only take into account the 2 first arguments (the function and
|
1272 |
the interval (a range). Managing the other arguments (the file for
|
1273 |
the proof and the exclusion intervals list) will be performed later
|
1274 |
Changes will be needed in sollya_lib.py file too.
|
1275 |
"""
|
1276 |
# Check that funcSa is a function.
|
1277 |
if not \ |
1278 |
sage.symbolic.callable.is_CallableSymbolicExpressionRing(parent(funcSa)): |
1279 |
return None |
1280 |
# Check that intervalSa is an interval.
|
1281 |
try:
|
1282 |
intervalSa.upper() |
1283 |
except AttributeError: |
1284 |
return None |
1285 |
# Convert the Sage function into a Sollya function.
|
1286 |
funcAsStringSa = funcSa._assume_str().replace('_SAGE_VAR_', '') |
1287 |
funcSo = pobyso_parse_string_sa_so(funcAsStringSa) |
1288 |
if not pobyso_obj_is_function_so_sa(funcSo): |
1289 |
sollya_lib_clear_obj(funcSo) |
1290 |
return None |
1291 |
# Convert the Sage interval into a Sollya range.
|
1292 |
rangeSo = pobyso_interval_to_range_sa_so(intervalSa) |
1293 |
retValSo = sollya_lib_infnorm(funcSo, rangeSo, None)
|
1294 |
sollya_lib_clear_obj(funcSo) |
1295 |
sollya_lib_clear_obj(rangeSo) |
1296 |
if pobyso_is_error_so_sa(retValSo):
|
1297 |
sollya_lib_clear_obj(retValSo) |
1298 |
return None |
1299 |
retValSa = pobyso_range_to_interval_so_sa(retValSo) |
1300 |
sollya_lib_clear_obj(retValSo) |
1301 |
return retValSa
|
1302 |
# End pobyso_infnorm_so_so.
|
1303 |
#
|
1304 |
def pobyso_infnorm_so_so(funcSo, rangeSo): |
1305 |
"""
|
1306 |
Very thin wrapper around sollya_lib_infnorm().
|
1307 |
We only take into account the 2 first arguments (the function and
|
1308 |
the interval (a range). Managing the other arguments (the file for
|
1309 |
the proof and the exclusion intervals list) will be performed later
|
1310 |
Changes will be needed in sollya_lib.py file too.
|
1311 |
|
1312 |
As per Sollya manual, this function should not be used anymore and
|
1313 |
supnorm should be called instead. Nevertheless, supnorm breaks
|
1314 |
sometimes whereas infnorm still returns a satisfactory answer.
|
1315 |
"""
|
1316 |
return sollya_lib_infnorm(funcSo, rangeSo, None) |
1317 |
# End pobyso_infnorm_so_so.
|
1318 |
|
1319 |
def pobyso_interval_to_range_sa_so(intervalSa, precisionSa=None): |
1320 |
if precisionSa is None: |
1321 |
precisionSa = intervalSa.parent().precision() |
1322 |
intervalSo = pobyso_bounds_to_range_sa_so(intervalSa.lower(),\ |
1323 |
intervalSa.upper(),\ |
1324 |
precisionSa) |
1325 |
return intervalSo
|
1326 |
# End pobyso_interval_to_range_sa_so
|
1327 |
|
1328 |
def pobyso_is_error_so_sa(objSo): |
1329 |
"""
|
1330 |
Thin wrapper around the sollya_lib_obj_is_error() function.
|
1331 |
"""
|
1332 |
if sollya_lib_obj_is_error(objSo) != 0: |
1333 |
return True |
1334 |
else:
|
1335 |
return False |
1336 |
# End pobyso_is_error-so_sa
|
1337 |
|
1338 |
def pobyso_is_floating_point_number_sa_sa(numberSa): |
1339 |
"""
|
1340 |
Check whether a Sage number is floating point.
|
1341 |
Exception stuff added because numbers other than
|
1342 |
floating-point ones do not have the is_real() attribute.
|
1343 |
"""
|
1344 |
try:
|
1345 |
return numberSa.is_real()
|
1346 |
except AttributeError: |
1347 |
return False |
1348 |
# End pobyso_is_floating_piont_number_sa_sa
|
1349 |
|
1350 |
def pobyso_is_range_so_sa(rangeCandidateSo): |
1351 |
"""
|
1352 |
Thin wrapper over sollya_lib_is_range.
|
1353 |
"""
|
1354 |
return sollya_lib_obj_is_range(rangeCandidateSo) != 0 |
1355 |
|
1356 |
# End pobyso_is_range_so_sa
|
1357 |
|
1358 |
|
1359 |
def pobyso_lib_init(): |
1360 |
sollya_lib_init(None)
|
1361 |
|
1362 |
def pobyso_lib_close(): |
1363 |
sollya_lib_close(None)
|
1364 |
|
1365 |
def pobyso_name_free_variable(freeVariableNameSa): |
1366 |
""" Legacy function. See pobyso_name_free_variable_sa_so. """
|
1367 |
pobyso_name_free_variable_sa_so(freeVariableNameSa) |
1368 |
|
1369 |
def pobyso_name_free_variable_sa_so(freeVariableNameSa): |
1370 |
"""
|
1371 |
Set the free variable name in Sollya from a Sage string.
|
1372 |
"""
|
1373 |
sollya_lib_name_free_variable(freeVariableNameSa) |
1374 |
|
1375 |
def pobyso_obj_is_function_so_sa(objSo): |
1376 |
"""
|
1377 |
Check if an object is a function.
|
1378 |
"""
|
1379 |
if sollya_lib_obj_is_function(objSo) != 0: |
1380 |
return True |
1381 |
else:
|
1382 |
return False |
1383 |
# End pobyso_obj_is_function_so_sa
|
1384 |
|
1385 |
def pobyso_obj_is_range_so_sa(objSo): |
1386 |
"""
|
1387 |
Check if an object is a function.
|
1388 |
"""
|
1389 |
if sollya_lib_obj_is_range(objSo) != 0: |
1390 |
return True |
1391 |
else:
|
1392 |
return False |
1393 |
# End pobyso_obj_is_range_so_sa
|
1394 |
|
1395 |
def pobyso_obj_is_string_so_sa(objSo): |
1396 |
"""
|
1397 |
Check if an object is a function.
|
1398 |
"""
|
1399 |
if sollya_lib_obj_is_string(objSo) != 0: |
1400 |
return True |
1401 |
else:
|
1402 |
return False |
1403 |
# End pobyso_obj_is_string_so_sa
|
1404 |
|
1405 |
def pobyso_parse_string(string): |
1406 |
""" Legacy function. See pobyso_parse_string_sa_so. """
|
1407 |
return pobyso_parse_string_sa_so(string)
|
1408 |
|
1409 |
def pobyso_parse_string_sa_so(string): |
1410 |
"""
|
1411 |
Get the Sollya expression computed from a Sage string or
|
1412 |
a Sollya error object if parsing failed.
|
1413 |
"""
|
1414 |
return sollya_lib_parse_string(string)
|
1415 |
|
1416 |
def pobyso_precision_so_sa(ctExpSo): |
1417 |
"""
|
1418 |
Computes the necessary precision to represent a number.
|
1419 |
If x is not zero, it can be uniquely written as x = m · 2e
|
1420 |
where m is an odd integer and e is an integer.
|
1421 |
precision(x) returns the number of bits necessary to write m
|
1422 |
in binary (i.e. ceil(log2(m))).
|
1423 |
"""
|
1424 |
#TODO: take care of the special case: 0, @NaN@, @Inf@
|
1425 |
precisionSo = sollya_lib_precision(ctExpSo) |
1426 |
precisionSa = pobyso_constant_from_int_so_sa(precisionSo) |
1427 |
sollya_lib_clear_obj(precisionSo) |
1428 |
return precisionSa
|
1429 |
# End pobyso_precision_so_sa
|
1430 |
|
1431 |
def pobyso_polynomial_coefficients_progressive_round_so_so(polySo, |
1432 |
funcSo, |
1433 |
icSo, |
1434 |
intervalSo, |
1435 |
itpSo, |
1436 |
ftpSo, |
1437 |
maxPrecSo, |
1438 |
maxErrSo, |
1439 |
debug=False):
|
1440 |
if debug:
|
1441 |
print "Input arguments:" |
1442 |
pobyso_autoprint(polySo) |
1443 |
pobyso_autoprint(funcSo) |
1444 |
pobyso_autoprint(icSo) |
1445 |
pobyso_autoprint(intervalSo) |
1446 |
pobyso_autoprint(itpSo) |
1447 |
pobyso_autoprint(ftpSo) |
1448 |
pobyso_autoprint(maxPrecSo) |
1449 |
pobyso_autoprint(maxErrSo) |
1450 |
print "________________" |
1451 |
|
1452 |
## Higher order function see:
|
1453 |
# http://effbot.org/pyfaq/how-do-you-make-a-higher-order-function-in-python.htm
|
1454 |
def precision_decay_ratio_function(degreeSa): |
1455 |
def outer(x): |
1456 |
def inner(x): |
1457 |
we = 3/8 |
1458 |
wq = 2/8 |
1459 |
a = 2.2
|
1460 |
b = 2
|
1461 |
return we*(exp(x/a)-1) + wq*((b*x)**2) + (1-we-wq)*x |
1462 |
return inner(x)/inner(degreeSa)
|
1463 |
return outer
|
1464 |
|
1465 |
#
|
1466 |
degreeSa = pobyso_polynomial_degree_so_sa(polySo) |
1467 |
ratio = precision_decay_ratio_function(degreeSa) |
1468 |
itpSa = pobyso_constant_from_int_so_sa(itpSo) |
1469 |
ftpSa = pobyso_constant_from_int_so_sa(ftpSo) |
1470 |
maxPrecSa = pobyso_constant_from_int_so_sa(maxPrecSo) |
1471 |
maxErrSa = pobyso_get_constant_as_rn_so_sa(maxErrSo) |
1472 |
if debug:
|
1473 |
print "degreeSa:", degreeSa |
1474 |
print "ratio:", ratio |
1475 |
print "itpsSa:", itpSa |
1476 |
print "ftpSa:", ftpSa |
1477 |
print "maxPrecSa:", maxPrecSa |
1478 |
print "maxErrSa:", maxErrSa |
1479 |
lastResPolySo = None
|
1480 |
lastInfNormSo = None
|
1481 |
#print "About to enter the while loop..."
|
1482 |
while True: |
1483 |
resPolySo = pobyso_constant_0_sa_so() |
1484 |
pDeltaSa = ftpSa - itpSa |
1485 |
for indexSa in reversed(xrange(0,degreeSa+1)): |
1486 |
#print "Index:", indexSa
|
1487 |
indexSo = pobyso_constant_from_int_sa_so(indexSa) |
1488 |
#pobyso_autoprint(indexSo)
|
1489 |
#print ratio(indexSa)
|
1490 |
ctpSa = floor(ftpSa - (pDeltaSa * ratio(indexSa))) |
1491 |
ctpSo = pobyso_constant_from_int_sa_so(ctpSa) |
1492 |
if debug:
|
1493 |
print "Index:", indexSa, " - Target precision:", |
1494 |
pobyso_autoprint(ctpSo) |
1495 |
cmonSo = \ |
1496 |
sollya_lib_build_function_mul(sollya_lib_coeff(polySo, indexSo), |
1497 |
sollya_lib_build_function_pow( \ |
1498 |
sollya_lib_build_function_free_variable(), \ |
1499 |
indexSo)) |
1500 |
#pobyso_autoprint(cmonSo)
|
1501 |
cmonrSo = pobyso_round_coefficients_single_so_so(cmonSo, ctpSo) |
1502 |
sollya_lib_clear_obj(cmonSo) |
1503 |
#pobyso_autoprint(cmonrSo)
|
1504 |
resPolySo = sollya_lib_build_function_add(resPolySo, |
1505 |
cmonrSo) |
1506 |
#pobyso_autoprint(resPolySo)
|
1507 |
# End for index
|
1508 |
freeVarSo = sollya_lib_build_function_free_variable() |
1509 |
changeVarSo = sollya_lib_sub(freeVarSo, icSo) |
1510 |
resPolyCvSo = sollya_lib_evaluate(resPolySo, changeVarSo) |
1511 |
errFuncSo = sollya_lib_build_function_sub(sollya_lib_copy_obj(funcSo), |
1512 |
resPolyCvSo) |
1513 |
infNormSo = sollya_lib_dirtyinfnorm(errFuncSo, intervalSo) |
1514 |
cerrSa = pobyso_get_constant_as_rn_so_sa(infNormSo) |
1515 |
if debug:
|
1516 |
print "Infnorm (Sollya):", |
1517 |
pobyso_autoprint(infNormSo) |
1518 |
sollya_lib_clear_obj(errFuncSo) |
1519 |
#print "Infnorm (Sage):", cerrSa
|
1520 |
if (cerrSa > maxErrSa):
|
1521 |
if debug:
|
1522 |
print "Error is too large." |
1523 |
if lastResPolySo is None: |
1524 |
if debug:
|
1525 |
print "Enlarging prec." |
1526 |
ntpSa = floor(ftpSa + ftpSa/50)
|
1527 |
## Can't enlarge (numerical)
|
1528 |
if ntpSa == ftpSa:
|
1529 |
sollya_lib_clear_obj(resPolySo) |
1530 |
return None |
1531 |
## Can't enlarge (not enough precision left)
|
1532 |
if ntpSa > maxPrecSa:
|
1533 |
sollya_lib_clear_obj(resPolySo) |
1534 |
return None |
1535 |
ftpSa = ntpSa |
1536 |
continue
|
1537 |
## One enlargement took place.
|
1538 |
else:
|
1539 |
if debug:
|
1540 |
print "Exit with the last before last polynomial." |
1541 |
print "Precision of highest degree monomial:", itpSa |
1542 |
print "Precision of constant term :", ftpSa |
1543 |
sollya_lib_clear_obj(resPolySo) |
1544 |
sollya_lib_clear_obj(infNormSo) |
1545 |
return (lastResPolySo, lastInfNormSo)
|
1546 |
# cerrSa <= maxErrSa: scrap more bits, possibly.
|
1547 |
else:
|
1548 |
if debug:
|
1549 |
print "Error is too small" |
1550 |
if cerrSa <= (maxErrSa/2): |
1551 |
if debug:
|
1552 |
print "Shrinking prec." |
1553 |
ntpSa = floor(ftpSa - ftpSa/50)
|
1554 |
## Can't shrink (numerical)
|
1555 |
if ntpSa == ftpSa:
|
1556 |
if not lastResPolySo is None: |
1557 |
sollya_lib_clear_obj(lastResPolySo) |
1558 |
if not lastInfNormSo is None: |
1559 |
sollya_lib_clear_obj(lastInfNormSo) |
1560 |
if debug:
|
1561 |
print "Exit because can't shrink anymore (numerically)." |
1562 |
print "Precision of highest degree monomial:", itpSa |
1563 |
print "Precision of constant term :", ftpSa |
1564 |
return (resPolySo, infNormSo)
|
1565 |
## Can't shrink (not enough precision left)
|
1566 |
if ntpSa <= itpSa:
|
1567 |
if not lastResPolySo is None: |
1568 |
sollya_lib_clear_obj(lastResPolySo) |
1569 |
if not lastInfNormSo is None: |
1570 |
sollya_lib_clear_obj(lastInfNormSo) |
1571 |
print "Exit because can't shrink anymore (no bits left)." |
1572 |
print "Precision of highest degree monomial:", itpSa |
1573 |
print "Precision of constant term :", ftpSa |
1574 |
return (resPolySo, infNormSo)
|
1575 |
ftpSa = ntpSa |
1576 |
if not lastResPolySo is None: |
1577 |
sollya_lib_clear_obj(lastResPolySo) |
1578 |
if not lastInfNormSo is None: |
1579 |
sollya_lib_clear_obj(lastInfNormSo) |
1580 |
lastResPolySo = resPolySo |
1581 |
lastInfNormSo = infNormSo |
1582 |
continue
|
1583 |
else: # Error is not that small, just return |
1584 |
if not lastResPolySo is None: |
1585 |
sollya_lib_clear_obj(lastResPolySo) |
1586 |
if not lastInfNormSo is None: |
1587 |
sollya_lib_clear_obj(lastInfNormSo) |
1588 |
if debug:
|
1589 |
print "Exit normally." |
1590 |
print "Precision of highest degree monomial:", itpSa |
1591 |
print "Precision of constant term :", ftpSa |
1592 |
return (resPolySo, infNormSo)
|
1593 |
# End wile True
|
1594 |
return None |
1595 |
# End pobyso_polynomial_coefficients_progressive_truncate_so_so.
|
1596 |
|
1597 |
def pobyso_polynomial_degree_so_sa(polySo): |
1598 |
"""
|
1599 |
Return the degree of a Sollya polynomial as a Sage int.
|
1600 |
"""
|
1601 |
degreeSo = sollya_lib_degree(polySo) |
1602 |
return pobyso_constant_from_int_so_sa(degreeSo)
|
1603 |
# End pobyso_polynomial_degree_so_sa
|
1604 |
|
1605 |
def pobyso_polynomial_degree_so_so(polySo): |
1606 |
"""
|
1607 |
Thin wrapper around lib_sollya_degree().
|
1608 |
"""
|
1609 |
return sollya_lib_degree(polySo)
|
1610 |
# End pobyso_polynomial_degree_so_so
|
1611 |
|
1612 |
def pobyso_range(rnLowerBound, rnUpperBound): |
1613 |
""" Legacy function. See pobyso_range_sa_so. """
|
1614 |
return pobyso_range_sa_so(rnLowerBound, rnUpperBound)
|
1615 |
|
1616 |
def pobyso_range_from_bounds_sa_so(rnLowerBound, rnUpperBound, precSa = None): |
1617 |
"""
|
1618 |
Create a Sollya range from 2 Sage real numbers as bounds
|
1619 |
"""
|
1620 |
# TODO check precision stuff.
|
1621 |
sollyaPrecChanged = False
|
1622 |
(initialSollyaPrecSo, initialSollyaPrecSa) = \ |
1623 |
pobyso_get_prec_so_so_sa() |
1624 |
if precSa is None: |
1625 |
precSa = max(rnLowerBound.parent().prec(), rnUpperBound.parent().prec())
|
1626 |
if precSa > initialSollyaPrecSa:
|
1627 |
if precSa <= 2: |
1628 |
print inspect.stack()[0][3], ": precision change <= 2 requested." |
1629 |
pobyso_set_prec_sa_so(precSa) |
1630 |
sollyaPrecChanged = True
|
1631 |
rangeSo = sollya_lib_range_from_bounds(get_rn_value(rnLowerBound), |
1632 |
get_rn_value(rnUpperBound)) |
1633 |
if sollyaPrecChanged:
|
1634 |
pobyso_set_prec_so_so(initialSollyaPrecSo) |
1635 |
sollya_lib_clear_obj(initialSollyaPrecSo) |
1636 |
return rangeSo
|
1637 |
# End pobyso_range_from_bounds_sa_so
|
1638 |
|
1639 |
def pobyso_range_list_so_sa(listSo): |
1640 |
"""
|
1641 |
Return a Sollya list of ranges as a Sage list of
|
1642 |
floating-point intervals.
|
1643 |
"""
|
1644 |
listSa = [] |
1645 |
## The function returns none if the list is empty or an error has happened.
|
1646 |
retVal = pobyso_get_list_elements_so_so(listSo) |
1647 |
if retVal is None: |
1648 |
return listSa
|
1649 |
## Just in case the interface is changed and an empty list is returned
|
1650 |
# instead of None.
|
1651 |
elif len(retVal) == 0: |
1652 |
return listSa
|
1653 |
else:
|
1654 |
## Remember pobyso_get_list_elements_so_so returns more information
|
1655 |
# than just the elements of the list (# elements, is_elliptic)
|
1656 |
listSaSo, numElements, isEndElliptic = retVal |
1657 |
## Return an empty list.
|
1658 |
if numElements == 0: |
1659 |
return listSa
|
1660 |
## Search first for the maximum precision of the elements
|
1661 |
maxPrecSa = 0
|
1662 |
for rangeSo in listSaSo: |
1663 |
#pobyso_autoprint(floatSo)
|
1664 |
curPrecSa = pobyso_get_prec_of_range_so_sa(rangeSo) |
1665 |
if curPrecSa > maxPrecSa:
|
1666 |
maxPrecSa = curPrecSa |
1667 |
##
|
1668 |
intervalField = RealIntervalField(maxPrecSa) |
1669 |
##
|
1670 |
for rangeSo in listSaSo: |
1671 |
listSa.append(pobyso_range_to_interval_so_sa(rangeSo, intervalField)) |
1672 |
return listSa
|
1673 |
# End pobyso_range_list_so_sa
|
1674 |
|
1675 |
def pobyso_range_max_abs_so_so(rangeSo): |
1676 |
"""
|
1677 |
Return, as Sollya constant, the maximum absolute value of a Sollay range.
|
1678 |
"""
|
1679 |
lowerBoundSo = sollya_lib_inf(rangeSo) |
1680 |
upperBoundSo = sollya_lib_sup(rangeSo) |
1681 |
#
|
1682 |
lowerBoundSo = sollya_lib_build_function_abs(lowerBoundSo) |
1683 |
upperBoundSo = sollya_lib_build_function_abs(upperBoundSo) |
1684 |
#pobyso_autoprint(lowerBoundSo)
|
1685 |
#pobyso_autoprint(upperBoundSo)
|
1686 |
#
|
1687 |
maxAbsSo = sollya_lib_max(lowerBoundSo, upperBoundSo, None)
|
1688 |
#sollya_lib_clear_obj(lowerBoundSo)
|
1689 |
#sollya_lib_clear_obj(upperBoundSo)
|
1690 |
return maxAbsSo
|
1691 |
# End pobyso_range_max_abs_so_so
|
1692 |
|
1693 |
def pobyso_range_to_interval_so_sa(rangeSo, realIntervalFieldSa = None): |
1694 |
"""
|
1695 |
Get a Sage interval from a Sollya range.
|
1696 |
If no realIntervalField is given as a parameter, the Sage interval
|
1697 |
precision is that of the Sollya range.
|
1698 |
Otherwise, the precision is that of the realIntervalField. In this case
|
1699 |
rounding may happen.
|
1700 |
"""
|
1701 |
if realIntervalFieldSa is None: |
1702 |
precSa = pobyso_get_prec_of_range_so_sa(rangeSo) |
1703 |
realIntervalFieldSa = RealIntervalField(precSa) |
1704 |
intervalSa = \ |
1705 |
pobyso_get_interval_from_range_so_sa(rangeSo, realIntervalFieldSa) |
1706 |
return intervalSa
|
1707 |
# End pobyso_range_to_interval_so_sa
|
1708 |
#
|
1709 |
def pobyso_relative_so_so(): |
1710 |
"""
|
1711 |
Very thin wrapper around the sollya_lib_relative function.
|
1712 |
"""
|
1713 |
return sollya_lib_relative()
|
1714 |
# End pobyso_relative_so_so
|
1715 |
#
|
1716 |
def pobyso_rat_poly_sa_so(polySa, precSa = None): |
1717 |
"""
|
1718 |
Create a Sollya polynomial from a Sage rational polynomial.
|
1719 |
We first convert the rational polynomial into a floating-point
|
1720 |
polynomial.
|
1721 |
"""
|
1722 |
## TODO: filter arguments.
|
1723 |
## Precision. If no precision is given, use the current precision
|
1724 |
# of Sollya.
|
1725 |
if precSa is None: |
1726 |
precSa = pobyso_get_prec_so_sa() |
1727 |
#print "Precision:", precSa
|
1728 |
RRR = RealField(precSa) |
1729 |
## Create a Sage polynomial in the "right" precision.
|
1730 |
P_RRR = RRR[polySa.variables()[0]]
|
1731 |
polyFloatSa = P_RRR(polySa) |
1732 |
## Make sure no precision is provided: pobyso_float_poly_sa_so will
|
1733 |
# recover it all by itself and will not make any extra conversion.
|
1734 |
return pobyso_float_poly_sa_so(polyFloatSa)
|
1735 |
|
1736 |
# End pobyso_rat_poly_sa_so
|
1737 |
|
1738 |
def pobyso_remez_canonical_sa_sa(func, \ |
1739 |
degree, \ |
1740 |
lowerBound, \ |
1741 |
upperBound, \ |
1742 |
weight = None, \
|
1743 |
quality = None):
|
1744 |
"""
|
1745 |
All arguments are Sage/Python.
|
1746 |
The functions (func and weight) must be passed as expressions or strings.
|
1747 |
Otherwise the function fails.
|
1748 |
The return value is a Sage polynomial.
|
1749 |
"""
|
1750 |
var('zorglub') # Dummy variable name for type check only. Type of |
1751 |
# zorglub is "symbolic expression".
|
1752 |
polySo = pobyso_remez_canonical_sa_so(func, \ |
1753 |
degree, \ |
1754 |
lowerBound, \ |
1755 |
upperBound, \ |
1756 |
weight, \ |
1757 |
quality) |
1758 |
# String test
|
1759 |
if parent(func) == parent("string"): |
1760 |
functionSa = eval(func)
|
1761 |
# Expression test.
|
1762 |
elif type(func) == type(zorglub): |
1763 |
functionSa = func |
1764 |
else:
|
1765 |
return None |
1766 |
#
|
1767 |
maxPrecision = 0
|
1768 |
if polySo is None: |
1769 |
return(None) |
1770 |
maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo) |
1771 |
RRRRSa = RealField(maxPrecision) |
1772 |
polynomialRingSa = RRRRSa[functionSa.variables()[0]]
|
1773 |
expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, RRRRSa) |
1774 |
polySa = polynomial(expSa, polynomialRingSa) |
1775 |
sollya_lib_clear_obj(polySo) |
1776 |
return(polySa)
|
1777 |
# End pobyso_remez_canonical_sa_sa
|
1778 |
|
1779 |
def pobyso_remez_canonical(func, \ |
1780 |
degree, \ |
1781 |
lowerBound, \ |
1782 |
upperBound, \ |
1783 |
weight = "1", \
|
1784 |
quality = None):
|
1785 |
""" Legacy function. See pobyso_remez_canonical_sa_so. """
|
1786 |
return(pobyso_remez_canonical_sa_so(func, \
|
1787 |
degree, \ |
1788 |
lowerBound, \ |
1789 |
upperBound, \ |
1790 |
weight, \ |
1791 |
quality)) |
1792 |
# End pobyso_remez_canonical.
|
1793 |
|
1794 |
def pobyso_remez_canonical_sa_so(func, \ |
1795 |
degree, \ |
1796 |
lowerBound, \ |
1797 |
upperBound, \ |
1798 |
weight = None, \
|
1799 |
quality = None):
|
1800 |
"""
|
1801 |
All arguments are Sage/Python.
|
1802 |
The functions (func and weight) must be passed as expressions or strings.
|
1803 |
Otherwise the function fails.
|
1804 |
The return value is a pointer to a Sollya function.
|
1805 |
lowerBound and upperBound mus be reals.
|
1806 |
"""
|
1807 |
var('zorglub') # Dummy variable name for type check only. Type of |
1808 |
# zorglub is "symbolic expression".
|
1809 |
currentVariableNameSa = None
|
1810 |
# The func argument can be of different types (string,
|
1811 |
# symbolic expression...)
|
1812 |
if parent(func) == parent("string"): |
1813 |
localFuncSa = sage_eval(func,globals())
|
1814 |
if len(localFuncSa.variables()) > 0: |
1815 |
currentVariableNameSa = localFuncSa.variables()[0]
|
1816 |
sollya_lib_name_free_variable(str(currentVariableNameSa))
|
1817 |
functionSo = \ |
1818 |
sollya_lib_parse_string(localFuncSa._assume_str().replace('_SAGE_VAR_', '')) |
1819 |
# Expression test.
|
1820 |
elif type(func) == type(zorglub): |
1821 |
# Until we are able to translate Sage expressions into Sollya
|
1822 |
# expressions : parse the string version.
|
1823 |
if len(func.variables()) > 0: |
1824 |
currentVariableNameSa = func.variables()[0]
|
1825 |
sollya_lib_name_free_variable(str(currentVariableNameSa))
|
1826 |
functionSo = \ |
1827 |
sollya_lib_parse_string(func._assume_str().replace('_SAGE_VAR_', '')) |
1828 |
else:
|
1829 |
return(None) |
1830 |
if weight is None: # No weight given -> 1. |
1831 |
weightSo = pobyso_constant_1_sa_so() |
1832 |
elif parent(weight) == parent("string"): # Weight given as string: parse it. |
1833 |
weightSo = sollya_lib_parse_string(func) |
1834 |
elif type(weight) == type(zorglub): # Weight given as symbolice expression. |
1835 |
functionSo = \ |
1836 |
sollya_lib_parse_string_sa_so(weight._assume_str().replace('_SAGE_VAR_', '')) |
1837 |
else:
|
1838 |
return(None) |
1839 |
degreeSo = pobyso_constant_from_int(degree) |
1840 |
rangeSo = pobyso_bounds_to_range_sa_so(lowerBound, upperBound) |
1841 |
if not quality is None: |
1842 |
qualitySo= pobyso_constant_sa_so(quality) |
1843 |
else:
|
1844 |
qualitySo = None
|
1845 |
|
1846 |
remezPolySo = sollya_lib_remez(functionSo, \ |
1847 |
degreeSo, \ |
1848 |
rangeSo, \ |
1849 |
weightSo, \ |
1850 |
qualitySo, \ |
1851 |
None)
|
1852 |
sollya_lib_clear_obj(functionSo) |
1853 |
sollya_lib_clear_obj(degreeSo) |
1854 |
sollya_lib_clear_obj(rangeSo) |
1855 |
sollya_lib_clear_obj(weightSo) |
1856 |
if not qualitySo is None: |
1857 |
sollya_lib_clear_obj(qualitySo) |
1858 |
return(remezPolySo)
|
1859 |
# End pobyso_remez_canonical_sa_so
|
1860 |
|
1861 |
def pobyso_remez_canonical_so_so(funcSo, \ |
1862 |
degreeSo, \ |
1863 |
rangeSo, \ |
1864 |
weightSo = pobyso_constant_1_sa_so(),\ |
1865 |
qualitySo = None):
|
1866 |
"""
|
1867 |
All arguments are pointers to Sollya objects.
|
1868 |
The return value is a pointer to a Sollya function.
|
1869 |
"""
|
1870 |
if not sollya_lib_obj_is_function(funcSo): |
1871 |
return(None) |
1872 |
return(sollya_lib_remez(funcSo, degreeSo, rangeSo, weightSo, qualitySo, None)) |
1873 |
# End pobyso_remez_canonical_so_so.
|
1874 |
|
1875 |
def pobyso_remez_exponents_list_sa_so(func, \ |
1876 |
exponentsList, \ |
1877 |
lowerBound, \ |
1878 |
upperBound, \ |
1879 |
weight = None, \
|
1880 |
quality = None):
|
1881 |
"""
|
1882 |
All arguments are Sage/Python.
|
1883 |
The functions (func and weight) must be passed as expressions or strings.
|
1884 |
Otherwise the function fails.
|
1885 |
The return value is a pointer to a Sollya function.
|
1886 |
lowerBound and upperBound mus be reals.
|
1887 |
"""
|
1888 |
var('zorglub') # Dummy variable name for type check only. Type of |
1889 |
# zorglub is "symbolic expression".
|
1890 |
currentVariableNameSa = None
|
1891 |
# The func argument can be of different types (string,
|
1892 |
# symbolic expression...)
|
1893 |
if parent(func) == parent("string"): |
1894 |
localFuncSa = sage_eval(func,globals())
|
1895 |
if len(localFuncSa.variables()) > 0: |
1896 |
currentVariableNameSa = localFuncSa.variables()[0]
|
1897 |
sollya_lib_name_free_variable(str(currentVariableNameSa))
|
1898 |
functionSo = \ |
1899 |
sollya_lib_parse_string(localFuncSa._assume_str().replace('_SAGE_VAR_', '')) |
1900 |
# Expression test.
|
1901 |
elif type(func) == type(zorglub): |
1902 |
# Until we are able to translate Sage expressions into Sollya
|
1903 |
# expressions : parse the string version.
|
1904 |
if len(func.variables()) > 0: |
1905 |
currentVariableNameSa = func.variables()[0]
|
1906 |
sollya_lib_name_free_variable(str(currentVariableNameSa))
|
1907 |
functionSo = \ |
1908 |
sollya_lib_parse_string(func._assume_str().replace('_SAGE_VAR_', '')) |
1909 |
else:
|
1910 |
return(None) |
1911 |
## Deal with the weight, much in the same way as with the function.
|
1912 |
if weight is None: # No weight given -> 1. |
1913 |
weightSo = pobyso_constant_1_sa_so() |
1914 |
elif parent(weight) == parent("string"): # Weight given as string: parse it. |
1915 |
weightSo = sollya_lib_parse_string(func) |
1916 |
elif type(weight) == type(zorglub): # Weight given as symbolice expression. |
1917 |
functionSo = \ |
1918 |
sollya_lib_parse_string(weight._assume_str().replace('_SAGE_VAR_', '',100)) |
1919 |
else:
|
1920 |
return(None) |
1921 |
rangeSo = pobyso_bounds_to_range_sa_so(lowerBound, upperBound) |
1922 |
if not quality is None: |
1923 |
qualitySo= pobyso_constant_sa_so(quality) |
1924 |
else:
|
1925 |
qualitySo = None
|
1926 |
#
|
1927 |
## Tranform the Sage list of exponents into a Sollya list.
|
1928 |
exponentsListSo = pobyso_build_list_of_ints_sa_so(*exponentsList) |
1929 |
remezPolySo = sollya_lib_remez(functionSo, \ |
1930 |
exponentsListSo, \ |
1931 |
rangeSo, \ |
1932 |
weightSo, \ |
1933 |
qualitySo, \ |
1934 |
None)
|
1935 |
sollya_lib_clear_obj(functionSo) |
1936 |
sollya_lib_clear_obj(exponentsListSo) |
1937 |
sollya_lib_clear_obj(rangeSo) |
1938 |
sollya_lib_clear_obj(weightSo) |
1939 |
if not qualitySo is None: |
1940 |
sollya_lib_clear_obj(qualitySo) |
1941 |
return(remezPolySo)
|
1942 |
# End pobyso_remez_exponentsList_sa_so
|
1943 |
#
|
1944 |
def pobyso_round_coefficients_so_so(polySo, truncFormatListSo): |
1945 |
"""
|
1946 |
A wrapper around the "classical" sollya_lib_roundcoefficients: a Sollya
|
1947 |
polynomial and a Sollya list are given as arguments.
|
1948 |
"""
|
1949 |
return sollya_lib_roundcoefficients(polySo, truncFormatListSo)
|
1950 |
|
1951 |
def pobyso_round_coefficients_progressive_so_so(polySo, |
1952 |
funcSo, |
1953 |
precSo, |
1954 |
intervalSo, |
1955 |
icSo, |
1956 |
currentApproxErrorSo, |
1957 |
approxAccurSo, |
1958 |
debug=False):
|
1959 |
"""
|
1960 |
From an input approximation polynomial, compute an output one with
|
1961 |
smaller coefficients and yet yields a sufficient approximation accuracy.
|
1962 |
"""
|
1963 |
if debug:
|
1964 |
print "Input arguments:" |
1965 |
print "Polynomial: ", ; pobyso_autoprint(polySo) |
1966 |
print "Function: ", ; pobyso_autoprint(funcSo) |
1967 |
print "Internal precision: ", ; pobyso_autoprint(precSo) |
1968 |
print "Interval: ", ; pobyso_autoprint(intervalSo) |
1969 |
print "Current approximation error: ", ; pobyso_autoprint(currentApproxErrorSo) |
1970 |
print "Requested approxiation error: ", ; pobyso_autoprint(approxAccurSo) |
1971 |
print "________________" |
1972 |
approxAccurSa = pobyso_get_constant_as_rn_so_sa(approxAccurSo) |
1973 |
currentApproxErrorSa = pobyso_get_constant_as_rn_so_sa(currentApproxErrorSo) |
1974 |
## If the current approximation error is too close to the target, there is
|
1975 |
# no possible gain.
|
1976 |
if currentApproxErrorSa >= approxAccurSa / 2: |
1977 |
#### Do not return the initial argument but copies: caller may free
|
1978 |
# the former as inutile after call.
|
1979 |
return (sollya_lib_copy_obj(polySo),
|
1980 |
sollya_lib_copy_obj(currentApproxErrorSo)) |
1981 |
#
|
1982 |
## Try to round the coefficients.
|
1983 |
degreeSa = pobyso_polynomial_degree_so_sa(polySo) |
1984 |
intervalSa = pobyso_range_to_interval_so_sa(intervalSo) |
1985 |
|
1986 |
if debug:
|
1987 |
print "degreeSa :", degreeSa |
1988 |
print "intervalSa :", intervalSa.str(style='brackets') |
1989 |
print "currentApproxErrorSa :", currentApproxErrorSa |
1990 |
print "approxAccurSa :", approxAccurSa |
1991 |
radiusSa = intervalSa.absolute_diameter() / 2
|
1992 |
if debug:
|
1993 |
print "log2(radius):", RR(radiusSa).log2() |
1994 |
iterIndex = 0
|
1995 |
## Build the "shaved" polynomial.
|
1996 |
while True: |
1997 |
### Start with a 0 value expression.
|
1998 |
resPolySo = pobyso_constant_0_sa_so() |
1999 |
roundedPolyApproxAccurSa = approxAccurSa / 2
|
2000 |
currentRadiusPowerSa = 1
|
2001 |
for degree in xrange(0,degreeSa + 1): |
2002 |
#### At round 0, use the agressive formula. At round 1, run the
|
2003 |
# proved formula.
|
2004 |
if iterIndex == 0: |
2005 |
roundingPowerSa = \ |
2006 |
floor(((currentRadiusPowerSa/roundedPolyApproxAccurSa)*(degree+1)).log2())
|
2007 |
else:
|
2008 |
roundingPowerSa = \ |
2009 |
floor(((currentRadiusPowerSa/roundedPolyApproxAccurSa)*(degreeSa+1)).log2())
|
2010 |
## Under extreme conditions the above formulas can evaluate under 2,
|
2011 |
# which is the minimal precision of an MPFR number.
|
2012 |
if roundingPowerSa < 2: |
2013 |
roundingPowerSa = 2
|
2014 |
if debug:
|
2015 |
print "roundedPolyApproxAccurSa", roundedPolyApproxAccurSa |
2016 |
print "currentRadiusPowerSa", currentRadiusPowerSa |
2017 |
print "Current rounding exponent:", roundingPowerSa |
2018 |
currentRadiusPowerSa *= radiusSa |
2019 |
index1So = pobyso_constant_from_int_sa_so(degree) |
2020 |
index2So = pobyso_constant_from_int_sa_so(degree) |
2021 |
### Create a monomial with:
|
2022 |
# - the coefficient in the initial monomial at the current degrree;
|
2023 |
# - the current exponent;
|
2024 |
# - the free variable.
|
2025 |
cmonSo = \ |
2026 |
sollya_lib_build_function_mul(sollya_lib_coeff(polySo, index1So), |
2027 |
sollya_lib_build_function_pow( \ |
2028 |
sollya_lib_build_function_free_variable(), \ |
2029 |
index2So)) |
2030 |
roundingPowerSo = pobyso_constant_from_int_sa_so(roundingPowerSa) |
2031 |
cmonrSo = pobyso_round_coefficients_single_so_so(cmonSo, roundingPowerSo) |
2032 |
sollya_lib_clear_obj(cmonSo) |
2033 |
### Add to the result polynomial.
|
2034 |
resPolySo = sollya_lib_build_function_add(resPolySo, |
2035 |
cmonrSo) |
2036 |
# End for.
|
2037 |
### Check the new polynomial.
|
2038 |
freeVarSo = sollya_lib_build_function_free_variable() |
2039 |
changeVarSo = sollya_lib_sub(freeVarSo, icSo) |
2040 |
resPolyCvSo = sollya_lib_evaluate(resPolySo, changeVarSo) |
2041 |
errFuncSo = sollya_lib_build_function_sub(sollya_lib_copy_obj(funcSo), |
2042 |
resPolyCvSo) |
2043 |
infNormSo = sollya_lib_dirtyinfnorm(errFuncSo, intervalSo) |
2044 |
### This also clears resPolyCvSo.
|
2045 |
sollya_lib_clear_obj(errFuncSo) |
2046 |
cerrSa = pobyso_get_constant_as_rn_so_sa(infNormSo) |
2047 |
if debug:
|
2048 |
print "Error of the new polynomial:", cerrSa |
2049 |
### If at round 1, return the initial polynomial error. This should
|
2050 |
# never happen since the rounding algorithm is proved. But some
|
2051 |
# circumstances may break it (e.g. internal precision of tools).
|
2052 |
if cerrSa > approxAccurSa:
|
2053 |
if iterIndex == 0: # Round 0 is agressive rounding, got round 1 (proved rounding) |
2054 |
sollya_lib_clear_obj(resPolySo) |
2055 |
sollya_lib_clear_obj(infNormSo) |
2056 |
iterIndex += 1
|
2057 |
continue
|
2058 |
else: # Round 1 and beyond : just return the oroginal polynomial. |
2059 |
sollya_lib_clear_obj(resPolySo) |
2060 |
sollya_lib_clear_obj(infNormSo) |
2061 |
#### Do not return the arguments but copies: the caller may free
|
2062 |
# free the former as inutile after call.
|
2063 |
return (sollya_lib_copy_obj(polySo),
|
2064 |
sollya_lib_copy_obj(currentApproxErrorSo)) |
2065 |
### If get here it is because cerrSa <= approxAccurSa
|
2066 |
### Approximation error of the new polynomial is acceptable.
|
2067 |
return (resPolySo, infNormSo)
|
2068 |
# End while True
|
2069 |
# End pobyso_round_coefficients_progressive_so_so
|
2070 |
|
2071 |
def pobyso_round_coefficients_single_so_so(polySo, commonPrecSo): |
2072 |
"""
|
2073 |
Create a rounded coefficients polynomial from polynomial argument to
|
2074 |
the number of bits in size argument.
|
2075 |
All coefficients are set to the same precision.
|
2076 |
"""
|
2077 |
## TODO: check arguments.
|
2078 |
endEllipListSo = pobyso_build_end_elliptic_list_so_so(commonPrecSo) |
2079 |
polySo = sollya_lib_roundcoefficients(polySo, endEllipListSo, None)
|
2080 |
sollya_lib_clear_obj(endEllipListSo) |
2081 |
#sollya_lib_clear_obj(endEllipListSo)
|
2082 |
return polySo
|
2083 |
|
2084 |
# End pobyso_round_coefficients_single_so_so
|
2085 |
|
2086 |
def pobyso_set_canonical_off(): |
2087 |
sollya_lib_set_canonical(sollya_lib_off()) |
2088 |
|
2089 |
def pobyso_set_canonical_on(): |
2090 |
sollya_lib_set_canonical(sollya_lib_on()) |
2091 |
|
2092 |
def pobyso_set_prec(p): |
2093 |
""" Legacy function. See pobyso_set_prec_sa_so. """
|
2094 |
pobyso_set_prec_sa_so(p) |
2095 |
|
2096 |
def pobyso_set_prec_sa_so(p): |
2097 |
#a = c_int(p)
|
2098 |
#precSo = c_void_p(sollya_lib_constant_from_int(a))
|
2099 |
#precSo = sollya_lib_constant_from_int(a)
|
2100 |
precSo = pobyso_constant_from_int_sa_so(p) |
2101 |
sollya_lib_set_prec(precSo) |
2102 |
sollya_lib_clear_obj(precSo) |
2103 |
# End pobyso_set_prec_sa_so.
|
2104 |
|
2105 |
def pobyso_set_prec_so_so(newPrecSo): |
2106 |
sollya_lib_set_prec(newPrecSo) |
2107 |
# End pobyso_set_prec_so_so.
|
2108 |
#
|
2109 |
def pobyso_supnorm_sa_sa(polySa, |
2110 |
funcSa, |
2111 |
intervalSa, |
2112 |
errorTypeSa=SOLLYA_ABSOLUTE, |
2113 |
accuracySa=RR(2^-40)): |
2114 |
"""
|
2115 |
An supnorm call with Sage arguments.
|
2116 |
"""
|
2117 |
# Check that polySa is a univariate polynomial. We only check here at
|
2118 |
# expression level, not at polynomial ring level
|
2119 |
## Make sure it is not a multivariate polynomial. The number of variables
|
2120 |
# can be zero in the case of a constant.
|
2121 |
try:
|
2122 |
if len(polySa.free_variables()) > 1 : |
2123 |
print "Invalid number of free variables in polySa:", polySa, |
2124 |
" -", polySa.free_variables()
|
2125 |
return None |
2126 |
except AttributeError: |
2127 |
print "polySa is not a SymbolicExpression." |
2128 |
return None |
2129 |
## Make sure it is a polynomial.
|
2130 |
if not polySa.is_polynomial(polySa.default_variable()): |
2131 |
print "polySa is not a polynomila expression." |
2132 |
return None |
2133 |
# Check that funcSa is a function.
|
2134 |
if not \ |
2135 |
sage.symbolic.callable.is_CallableSymbolicExpressionRing(parent(funcSa)): |
2136 |
return None |
2137 |
# Check that intervalSa is an interval.
|
2138 |
try:
|
2139 |
intervalSa.upper() |
2140 |
except AttributeError: |
2141 |
print "intervalSa is not an interval." |
2142 |
return None |
2143 |
# Convert the Sage polynomial into a Sollya polynomial.
|
2144 |
polyAsStringSa = polySa._assume_str().replace('_SAGE_VAR_', '') |
2145 |
polySo = pobyso_parse_string_sa_so(polyAsStringSa) |
2146 |
if not pobyso_obj_is_function_so_sa(polySo): |
2147 |
sollya_lib_clear_obj(polySo) |
2148 |
print "The Sollya object created from the polynomial is not a function." |
2149 |
return None |
2150 |
#pobyso_autoprint(polySo)
|
2151 |
# Convert the Sage function into a Sollya function.
|
2152 |
funcAsStringSa = funcSa._assume_str().replace('_SAGE_VAR_', '') |
2153 |
funcSo = pobyso_parse_string_sa_so(funcAsStringSa) |
2154 |
if not pobyso_obj_is_function_so_sa(funcSo): |
2155 |
sollya_lib_clear_obj(polySo) |
2156 |
sollya_lib_clear_obj(funcSo) |
2157 |
print "The Sollya object created from the function is not a function." |
2158 |
return None |
2159 |
#pobyso_autoprint(funcSo)
|
2160 |
# Convert the Sage interval into a Sollya range.
|
2161 |
rangeSo = pobyso_interval_to_range_sa_so(intervalSa) |
2162 |
if not pobyso_is_range_so_sa(rangeSo): |
2163 |
sollya_lib_clear_obj(polySo) |
2164 |
sollya_lib_clear_obj(funcSo) |
2165 |
sollya_lib_clear_obj(rangeSo) |
2166 |
print "The Sollya object created from the interval is not a range." |
2167 |
return None |
2168 |
#pobyso_autoprint(rangeSo)
|
2169 |
# Check the error type. We do not check the returned object since we
|
2170 |
# trust our code and Sollya on this one.
|
2171 |
if errorTypeSa != SOLLYA_ABSOLUTE and errorTypeSa != SOLLYA_RELATIVE: |
2172 |
sollya_lib_clear_obj(polySo) |
2173 |
sollya_lib_clear_obj(funcSo) |
2174 |
sollya_lib_clear_obj(rangeSo) |
2175 |
return None |
2176 |
if errorTypeSa == SOLLYA_ABSOLUTE:
|
2177 |
errorTypeSo = pobyso_absolute_so_so() |
2178 |
else:
|
2179 |
errorTypeSo = pobyso_relative_so_so() |
2180 |
#pobyso_autoprint(errorTypeSo)
|
2181 |
# Check if accuracySa is an element of a RealField. If not, try to
|
2182 |
# convert it.
|
2183 |
try:
|
2184 |
accuracySa.ulp() |
2185 |
except AttributeError: |
2186 |
try:
|
2187 |
accuracySa = RR(accuracySa) |
2188 |
except TypeError: |
2189 |
accuracySa = RR(str(accuracySa))
|
2190 |
# Create the Sollya object
|
2191 |
accuracySo = pobyso_constant_sa_so(accuracySa) |
2192 |
#pobyso_autoprint(accuracySo)
|
2193 |
if pobyso_is_error_so_sa(accuracySo):
|
2194 |
sollya_lib_clear_obj(polySo) |
2195 |
sollya_lib_clear_obj(funcSo) |
2196 |
sollya_lib_clear_obj(rangeSo) |
2197 |
sollya_lib_clear_obj(errorTypeSo) |
2198 |
sollya_lib_clear_obj(accuracySo) |
2199 |
return None |
2200 |
retValSo = sollya_lib_supnorm(polySo, |
2201 |
funcSo, |
2202 |
rangeSo, |
2203 |
errorTypeSo, |
2204 |
accuracySo, |
2205 |
None)
|
2206 |
sollya_lib_clear_obj(polySo) |
2207 |
sollya_lib_clear_obj(funcSo) |
2208 |
sollya_lib_clear_obj(rangeSo) |
2209 |
sollya_lib_clear_obj(errorTypeSo) |
2210 |
sollya_lib_clear_obj(accuracySo) |
2211 |
if pobyso_is_error_so_sa(retValSo):
|
2212 |
sollya_lib_clear_obj(retValSo) |
2213 |
return None |
2214 |
#pobyso_autoprint(retValSo)
|
2215 |
retValSa = pobyso_range_to_interval_so_sa(retValSo) |
2216 |
sollya_lib_clear_obj(retValSo) |
2217 |
return retValSa
|
2218 |
# End pobyso_supnorm_sa_sa
|
2219 |
|
2220 |
def pobyso_supnorm_so_sa(polySo, funcSo, intervalSo, errorTypeSo = None,\ |
2221 |
accuracySo = None, realFieldSa = None): |
2222 |
"""
|
2223 |
Computes the supremum norm from Sollya input arguments and returns a
|
2224 |
Sage floating-point number whose precision is set by the only Sage argument.
|
2225 |
|
2226 |
The returned value is the maximum of the absolute values of the range
|
2227 |
elements returned by the Sollya supnorm functions
|
2228 |
"""
|
2229 |
supNormRangeSo = pobyso_supnorm_so_so(polySo, |
2230 |
funcSo, |
2231 |
intervalSo, |
2232 |
errorTypeSo, |
2233 |
accuracySo) |
2234 |
supNormSo = pobyso_range_max_abs_so_so(supNormRangeSo) |
2235 |
sollya_lib_clear_obj(supNormRangeSo) |
2236 |
#pobyso_autoprint(supNormSo)
|
2237 |
supNormSa = pobyso_get_constant_as_rn_with_rf_so_sa(supNormSo, realFieldSa) |
2238 |
sollya_lib_clear_obj(supNormSo) |
2239 |
return supNormSa
|
2240 |
# End pobyso_supnorm_so_sa.
|
2241 |
#
|
2242 |
def pobyso_supnorm_so_so(polySo, funcSo, intervalSo, errorTypeSo = None,\ |
2243 |
accuracySo = None):
|
2244 |
"""
|
2245 |
Computes the supnorm of the approximation error between the given
|
2246 |
polynomial and function. Attention: returns a range!
|
2247 |
errorTypeSo defaults to "absolute".
|
2248 |
accuracySo defaults to 2^(-40).
|
2249 |
"""
|
2250 |
if errorTypeSo is None: |
2251 |
errorTypeSo = sollya_lib_absolute(None)
|
2252 |
errorTypeIsNone = True
|
2253 |
else:
|
2254 |
errorTypeIsNone = False
|
2255 |
#
|
2256 |
if accuracySo is None: |
2257 |
# Notice the **: we are in Pythonland here!
|
2258 |
accuracySo = pobyso_constant_sa_so(RR(2**(-40))) |
2259 |
accuracyIsNone = True
|
2260 |
else:
|
2261 |
accuracyIsNone = False
|
2262 |
#pobyso_autoprint(accuracySo)
|
2263 |
resultSo = \ |
2264 |
sollya_lib_supnorm(polySo, funcSo, intervalSo, errorTypeSo, \ |
2265 |
accuracySo) |
2266 |
if errorTypeIsNone:
|
2267 |
sollya_lib_clear_obj(errorTypeSo) |
2268 |
if accuracyIsNone:
|
2269 |
sollya_lib_clear_obj(accuracySo) |
2270 |
return resultSo
|
2271 |
# End pobyso_supnorm_so_so
|
2272 |
#
|
2273 |
def pobyso_taylor_expansion_no_change_var_so_so(functionSo, |
2274 |
degreeSo, |
2275 |
rangeSo, |
2276 |
errorTypeSo=None,
|
2277 |
sollyaPrecSo=None):
|
2278 |
"""
|
2279 |
Compute the Taylor expansion without the variable change
|
2280 |
x -> x-intervalCenter.
|
2281 |
If errorTypeSo is None, absolute is used.
|
2282 |
If sollyaPrecSo is None, Sollya internal precision is not changed.
|
2283 |
"""
|
2284 |
# Change internal Sollya precision, if needed.
|
2285 |
(initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa() |
2286 |
sollyaPrecChanged = False
|
2287 |
if sollyaPrecSo is None: |
2288 |
pass
|
2289 |
else:
|
2290 |
sollya_lib_set_prec(sollyaPrecSo) |
2291 |
sollyaPrecChanged = True
|
2292 |
# Error type stuff: default to absolute.
|
2293 |
if errorTypeSo is None: |
2294 |
errorTypeIsNone = True
|
2295 |
errorTypeSo = sollya_lib_absolute(None)
|
2296 |
else:
|
2297 |
errorTypeIsNone = False
|
2298 |
intervalCenterSo = sollya_lib_mid(rangeSo, None)
|
2299 |
taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo, |
2300 |
intervalCenterSo, |
2301 |
rangeSo, errorTypeSo, None)
|
2302 |
# Object taylorFormListSaSo is a Python list of Sollya objects references
|
2303 |
# that are copies of the elements of taylorFormSo.
|
2304 |
# pobyso_get_list_elements_so_so clears taylorFormSo.
|
2305 |
(taylorFormListSaSo, numElementsSa, isEndEllipticSa) = \ |
2306 |
pobyso_get_list_elements_so_so(taylorFormSo) |
2307 |
## Copy needed here since polySo will be returned and taylorFormListSaSo
|
2308 |
# will be cleared.
|
2309 |
polySo = sollya_lib_copy_obj(taylorFormListSaSo[0])
|
2310 |
#print "Num elements:", numElementsSa
|
2311 |
sollya_lib_clear_obj(taylorFormSo) |
2312 |
# No copy_obj needed here: a new objects are created.
|
2313 |
maxErrorSo = sollya_lib_sup(taylorFormListSaSo[2])
|
2314 |
minErrorSo = sollya_lib_inf(taylorFormListSaSo[2])
|
2315 |
# List taylorFormListSaSo is not needed anymore.
|
2316 |
pobyso_clear_list_elements_sa_so(taylorFormListSaSo) |
2317 |
absMaxErrorSo = sollya_lib_abs(maxErrorSo) |
2318 |
absMinErrorSo = sollya_lib_abs(minErrorSo) |
2319 |
sollya_lib_clear_obj(maxErrorSo) |
2320 |
sollya_lib_clear_obj(minErrorSo) |
2321 |
absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo) |
2322 |
absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo) |
2323 |
#
|
2324 |
if errorTypeIsNone:
|
2325 |
sollya_lib_clear_obj(errorTypeSo) |
2326 |
## If changed, reset the Sollya working precision.
|
2327 |
if sollyaPrecChanged:
|
2328 |
sollya_lib_set_prec(initialSollyaPrecSo) |
2329 |
sollya_lib_clear_obj(initialSollyaPrecSo) |
2330 |
## According to what error is the largest, return the errors.
|
2331 |
if absMaxErrorSa > absMinErrorSa:
|
2332 |
sollya_lib_clear_obj(absMinErrorSo) |
2333 |
return (polySo, intervalCenterSo, absMaxErrorSo)
|
2334 |
else:
|
2335 |
sollya_lib_clear_obj(absMaxErrorSo) |
2336 |
return (polySo, intervalCenterSo, absMinErrorSo)
|
2337 |
# end pobyso_taylor_expansion_no_change_var_so_so
|
2338 |
|
2339 |
def pobyso_taylor_expansion_with_change_var_so_so(functionSo, degreeSo, \ |
2340 |
rangeSo, \ |
2341 |
errorTypeSo=None, \
|
2342 |
sollyaPrecSo=None):
|
2343 |
"""
|
2344 |
Compute the Taylor expansion with the variable change
|
2345 |
x -> (x-intervalCenter) included.
|
2346 |
"""
|
2347 |
# Change Sollya internal precision, if need.
|
2348 |
sollyaPrecChanged = False
|
2349 |
(initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa() |
2350 |
if sollyaPrecSo is None: |
2351 |
pass
|
2352 |
else:
|
2353 |
sollya_lib_set_prec(sollyaPrecSo) |
2354 |
sollyaPrecChanged = True
|
2355 |
#
|
2356 |
# Error type stuff: default to absolute.
|
2357 |
if errorTypeSo is None: |
2358 |
errorTypeIsNone = True
|
2359 |
errorTypeSo = sollya_lib_absolute(None)
|
2360 |
else:
|
2361 |
errorTypeIsNone = False
|
2362 |
intervalCenterSo = sollya_lib_mid(rangeSo) |
2363 |
taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo, \ |
2364 |
intervalCenterSo, \ |
2365 |
rangeSo, errorTypeSo, None)
|
2366 |
# Object taylorFormListSaSo is a Python list of Sollya objects references
|
2367 |
# that are copies of the elements of taylorFormSo.
|
2368 |
# pobyso_get_list_elements_so_so clears taylorFormSo.
|
2369 |
(taylorFormListSaSo, numElements, isEndElliptic) = \ |
2370 |
pobyso_get_list_elements_so_so(taylorFormSo) |
2371 |
sollya_lib_clear_obj(taylorFormSo) |
2372 |
polySo = taylorFormListSaSo[0]
|
2373 |
## Maximum error computation with taylorFormListSaSo[2], a range
|
2374 |
# holding the actual error. Bounds can be negative.
|
2375 |
maxErrorSo = sollya_lib_sup(taylorFormListSaSo[2])
|
2376 |
minErrorSo = sollya_lib_inf(taylorFormListSaSo[2])
|
2377 |
absMaxErrorSo = sollya_lib_abs(maxErrorSo) |
2378 |
absMinErrorSo = sollya_lib_abs(minErrorSo) |
2379 |
sollya_lib_clear_obj(maxErrorSo) |
2380 |
sollya_lib_clear_obj(minErrorSo) |
2381 |
absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo) |
2382 |
absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo) |
2383 |
changeVarExpSo = sollya_lib_build_function_sub(\ |
2384 |
sollya_lib_build_function_free_variable(),\ |
2385 |
sollya_lib_copy_obj(intervalCenterSo)) |
2386 |
polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpSo) |
2387 |
# List taylorFormListSaSo is not needed anymore.
|
2388 |
pobyso_clear_list_elements_sa_so(taylorFormListSaSo) |
2389 |
sollya_lib_clear_obj(changeVarExpSo) |
2390 |
# If changed, reset the Sollya working precision.
|
2391 |
if sollyaPrecChanged:
|
2392 |
sollya_lib_set_prec(initialSollyaPrecSo) |
2393 |
sollya_lib_clear_obj(initialSollyaPrecSo) |
2394 |
if errorTypeIsNone:
|
2395 |
sollya_lib_clear_obj(errorTypeSo) |
2396 |
# Do not clear maxErrorSo.
|
2397 |
if absMaxErrorSa > absMinErrorSa:
|
2398 |
sollya_lib_clear_obj(absMinErrorSo) |
2399 |
return((polyVarChangedSo, intervalCenterSo, absMaxErrorSo))
|
2400 |
else:
|
2401 |
sollya_lib_clear_obj(absMaxErrorSo) |
2402 |
return((polyVarChangedSo, intervalCenterSo, absMinErrorSo))
|
2403 |
# end pobyso_taylor_expansion_with_change_var_so_so
|
2404 |
|
2405 |
def pobyso_taylor(function, degree, point): |
2406 |
""" Legacy function. See pobysoTaylor_so_so. """
|
2407 |
return(pobyso_taylor_so_so(function, degree, point))
|
2408 |
|
2409 |
def pobyso_taylor_so_so(functionSo, degreeSo, pointSo): |
2410 |
return(sollya_lib_taylor(functionSo, degreeSo, pointSo))
|
2411 |
|
2412 |
def pobyso_taylorform(function, degree, point = None, |
2413 |
interval = None, errorType=None): |
2414 |
""" Legacy function. See pobyso_taylorform_sa_sa;"""
|
2415 |
|
2416 |
def pobyso_taylorform_sa_sa(functionSa, \ |
2417 |
degreeSa, \ |
2418 |
pointSa, \ |
2419 |
intervalSa=None, \
|
2420 |
errorTypeSa=None, \
|
2421 |
precisionSa=None):
|
2422 |
"""
|
2423 |
Compute the Taylor form of 'degreeSa' for 'functionSa' at 'pointSa'
|
2424 |
for 'intervalSa' with 'errorTypeSa' (a string) using 'precisionSa'.
|
2425 |
point: must be a Real or a Real interval.
|
2426 |
return the Taylor form as an array
|
2427 |
TODO: take care of the interval and of the point when it is an interval;
|
2428 |
when errorType is not None;
|
2429 |
take care of the other elements of the Taylor form (coefficients
|
2430 |
errors and delta.
|
2431 |
"""
|
2432 |
# Absolute as the default error.
|
2433 |
if errorTypeSa is None: |
2434 |
errorTypeSo = sollya_lib_absolute() |
2435 |
elif errorTypeSa == "relative": |
2436 |
errorTypeSo = sollya_lib_relative() |
2437 |
elif errortypeSa == "absolute": |
2438 |
errorTypeSo = sollya_lib_absolute() |
2439 |
else:
|
2440 |
# No clean up needed.
|
2441 |
return None |
2442 |
# Global precision stuff
|
2443 |
sollyaPrecisionChangedSa = False
|
2444 |
(initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa() |
2445 |
if precisionSa is None: |
2446 |
precSa = initialSollyaPrecSa |
2447 |
else:
|
2448 |
if precSa > initialSollyaPrecSa:
|
2449 |
if precSa <= 2: |
2450 |
print inspect.stack()[0][3], ":precision change <= 2 requested." |
2451 |
pobyso_set_prec_sa_so(precSa) |
2452 |
sollyaPrecisionChangedSa = True
|
2453 |
#
|
2454 |
if len(functionSa.variables()) > 0: |
2455 |
varSa = functionSa.variables()[0]
|
2456 |
pobyso_name_free_variable_sa_so(str(varSa))
|
2457 |
# In any case (point or interval) the parent of pointSa has a precision
|
2458 |
# method.
|
2459 |
pointPrecSa = pointSa.parent().precision() |
2460 |
if precSa > pointPrecSa:
|
2461 |
pointPrecSa = precSa |
2462 |
# In any case (point or interval) pointSa has a base_ring() method.
|
2463 |
pointBaseRingString = str(pointSa.base_ring())
|
2464 |
if re.search('Interval', pointBaseRingString) is None: # Point |
2465 |
pointSo = pobyso_constant_sa_so(pointSa, pointPrecSa) |
2466 |
else: # Interval. |
2467 |
pointSo = pobyso_interval_to_range_sa_so(pointSa, pointPrecSa) |
2468 |
# Sollyafy the function.
|
2469 |
functionSo = pobyso_parse_string_sa_so(functionSa._assume_str().replace('_SAGE_VAR_', '')) |
2470 |
if sollya_lib_obj_is_error(functionSo):
|
2471 |
print "pobyso_tailorform: function string can't be parsed!" |
2472 |
return None |
2473 |
# Sollyafy the degree
|
2474 |
degreeSo = sollya_lib_constant_from_int(int(degreeSa))
|
2475 |
# Sollyafy the point
|
2476 |
# Call Sollya
|
2477 |
taylorFormSo = \ |
2478 |
sollya_lib_taylorform(functionSo, degreeSo, pointSo, errorTypeSo,\ |
2479 |
None)
|
2480 |
sollya_lib_clear_obj(functionSo) |
2481 |
sollya_lib_clear_obj(degreeSo) |
2482 |
sollya_lib_clear_obj(pointSo) |
2483 |
sollya_lib_clear_obj(errorTypeSo) |
2484 |
(tfsAsList, numElements, isEndElliptic) = \ |
2485 |
pobyso_get_list_elements_so_so(taylorFormSo) |
2486 |
polySo = tfsAsList[0]
|
2487 |
maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo) |
2488 |
polyRealField = RealField(maxPrecision) |
2489 |
expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, polyRealField) |
2490 |
if sollyaPrecisionChangedSa:
|
2491 |
sollya_lib_set_prec(initialSollyaPrecSo) |
2492 |
sollya_lib_clear_obj(initialSollyaPrecSo) |
2493 |
polynomialRing = polyRealField[str(varSa)]
|
2494 |
polySa = polynomial(expSa, polynomialRing) |
2495 |
taylorFormSa = [polySa] |
2496 |
# Final clean-up.
|
2497 |
sollya_lib_clear_obj(taylorFormSo) |
2498 |
return(taylorFormSa)
|
2499 |
# End pobyso_taylor_form_sa_sa
|
2500 |
|
2501 |
def pobyso_taylorform_so_so(functionSo, degreeSo, pointSo, intervalSo=None, \ |
2502 |
errorTypeSo=None):
|
2503 |
createdErrorType = False
|
2504 |
if errorTypeSo is None: |
2505 |
errorTypeSo = sollya_lib_absolute() |
2506 |
createdErrorType = True
|
2507 |
else:
|
2508 |
#TODO: deal with the other case.
|
2509 |
pass
|
2510 |
if intervalSo is None: |
2511 |
resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \ |
2512 |
errorTypeSo, None)
|
2513 |
else:
|
2514 |
resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \ |
2515 |
intervalSo, errorTypeSo, None)
|
2516 |
if createdErrorType:
|
2517 |
sollya_lib_clear_obj(errorTypeSo) |
2518 |
return resultSo
|
2519 |
|
2520 |
|
2521 |
def pobyso_univar_polynomial_print_reverse(polySa): |
2522 |
""" Legacy function. See pobyso_univar_polynomial_print_reverse_sa_sa. """
|
2523 |
return(pobyso_univar_polynomial_print_reverse_sa_sa(polySa))
|
2524 |
|
2525 |
def pobyso_univar_polynomial_print_reverse_sa_sa(polySa): |
2526 |
"""
|
2527 |
Return the string representation of a univariate polynomial with
|
2528 |
monomials ordered in the x^0..x^n order of the monomials.
|
2529 |
Remember: Sage
|
2530 |
"""
|
2531 |
polynomialRing = polySa.base_ring() |
2532 |
# A very expensive solution:
|
2533 |
# -create a fake multivariate polynomial field with only one variable,
|
2534 |
# specifying a negative lexicographical order;
|
2535 |
mpolynomialRing = PolynomialRing(polynomialRing.base(), \ |
2536 |
polynomialRing.variable_name(), \ |
2537 |
1, order='neglex') |
2538 |
# - convert the univariate argument polynomial into a multivariate
|
2539 |
# version;
|
2540 |
p = mpolynomialRing(polySa) |
2541 |
# - return the string representation of the converted form.
|
2542 |
# There is no simple str() method defined for p's class.
|
2543 |
return(p.__str__())
|
2544 |
#
|
2545 |
#print pobyso_get_prec()
|
2546 |
pobyso_set_prec(165)
|
2547 |
#print pobyso_get_prec()
|
2548 |
#a=100
|
2549 |
#print type(a)
|
2550 |
#id(a)
|
2551 |
#print "Max arity: ", pobyso_max_arity
|
2552 |
#print "Function tripleDouble (43) as a string: ", pobyso_function_type_as_string(43)
|
2553 |
#print "Function None (44) as a string: ", pobyso_function_type_as_string(44)
|
2554 |
sys.stderr.write("\t...Pobyso check done.\n")
|