root / pobysoPythonSage / src / pobyso.py @ 35
Historique | Voir | Annoter | Télécharger (14,31 ko)
1 |
"""
|
---|---|
2 |
Actual functions to use in Sage
|
3 |
ST 2012-11-13
|
4 |
|
5 |
Command line syntax:
|
6 |
use from Sage (via the "load" or the "attach" commands)
|
7 |
|
8 |
NOTES:
|
9 |
Reported errors in Eclipse come from the calls to
|
10 |
the Sollya library
|
11 |
|
12 |
ToDo (among other things):
|
13 |
-memory management.
|
14 |
"""
|
15 |
from ctypes import * |
16 |
|
17 |
(SOLLYA_BASE_FUNC_ABS, |
18 |
SOLLYA_BASE_FUNC_ACOS, |
19 |
SOLLYA_BASE_FUNC_ACOSH, |
20 |
SOLLYA_BASE_FUNC_ADD, |
21 |
SOLLYA_BASE_FUNC_ASIN, |
22 |
SOLLYA_BASE_FUNC_ASINH, |
23 |
SOLLYA_BASE_FUNC_ATAN, |
24 |
SOLLYA_BASE_FUNC_ATANH, |
25 |
SOLLYA_BASE_FUNC_CEIL, |
26 |
SOLLYA_BASE_FUNC_CONSTANT, |
27 |
SOLLYA_BASE_FUNC_COS, |
28 |
SOLLYA_BASE_FUNC_COSH, |
29 |
SOLLYA_BASE_FUNC_DIV, |
30 |
SOLLYA_BASE_FUNC_DOUBLE, |
31 |
SOLLYA_BASE_FUNC_DOUBLEDOUBLE, |
32 |
SOLLYA_BASE_FUNC_DOUBLEEXTENDED, |
33 |
SOLLYA_BASE_FUNC_ERF, |
34 |
SOLLYA_BASE_FUNC_ERFC, |
35 |
SOLLYA_BASE_FUNC_EXP, |
36 |
SOLLYA_BASE_FUNC_EXP_M1, |
37 |
SOLLYA_BASE_FUNC_FLOOR, |
38 |
SOLLYA_BASE_FUNC_FREE_VARIABLE, |
39 |
SOLLYA_BASE_FUNC_HALFPRECISION, |
40 |
SOLLYA_BASE_FUNC_LIBRARYCONSTANT, |
41 |
SOLLYA_BASE_FUNC_LIBRARYFUNCTION, |
42 |
SOLLYA_BASE_FUNC_LOG, |
43 |
SOLLYA_BASE_FUNC_LOG_10, |
44 |
SOLLYA_BASE_FUNC_LOG_1P, |
45 |
SOLLYA_BASE_FUNC_LOG_2, |
46 |
SOLLYA_BASE_FUNC_MUL, |
47 |
SOLLYA_BASE_FUNC_NEARESTINT, |
48 |
SOLLYA_BASE_FUNC_NEG, |
49 |
SOLLYA_BASE_FUNC_PI, |
50 |
SOLLYA_BASE_FUNC_POW, |
51 |
SOLLYA_BASE_FUNC_PROCEDUREFUNCTION, |
52 |
SOLLYA_BASE_FUNC_QUAD, |
53 |
SOLLYA_BASE_FUNC_SIN, |
54 |
SOLLYA_BASE_FUNC_SINGLE, |
55 |
SOLLYA_BASE_FUNC_SINH, |
56 |
SOLLYA_BASE_FUNC_SQRT, |
57 |
SOLLYA_BASE_FUNC_SUB, |
58 |
SOLLYA_BASE_FUNC_TAN, |
59 |
SOLLYA_BASE_FUNC_TANH, |
60 |
SOLLYA_BASE_FUNC_TRIPLEDOUBLE) = map(int,xrange(44)) |
61 |
print "First constant - SOLLYA_BASE_FUNC_ABS: ", SOLLYA_BASE_FUNC_ABS |
62 |
print "Last constant - SOLLYA_BASE_FUNC_TRIPLEDOUBLE: ", SOLLYA_BASE_FUNC_TRIPLEDOUBLE |
63 |
|
64 |
pobyso_max_arity = 9
|
65 |
|
66 |
def pobyso_autoprint(arg): |
67 |
sollya_lib_autoprint(arg,None)
|
68 |
|
69 |
def pobyso_cmp(rnArg, soCte): |
70 |
precisionOfCte = c_int(0)
|
71 |
# From the Sollya constant, create a local Sage RealNumber.
|
72 |
sollya_lib_get_prec_of_constant(precisionOfCte, soCte) |
73 |
#print "Precision of constant: ", precisionOfCte
|
74 |
RRRR = RealField(precisionOfCte.value) |
75 |
rnLocal = RRRR(0)
|
76 |
sollya_lib_get_constant(get_rn_value(rnLocal), soCte) |
77 |
#print "rnDummy: ", rnDummy
|
78 |
# Compare the local Sage RealNumber with rnArg.
|
79 |
return(cmp_rn_value(rnArg, rnLocal))
|
80 |
|
81 |
def pobyso_constant(rnArg): |
82 |
return (sollya_lib_constant(get_rn_value(rnArg)))
|
83 |
|
84 |
def pobyso_constant_1(): |
85 |
return(pobyso_constant_from_int(1)) |
86 |
|
87 |
def pobyso_constant_from_int(anInt): |
88 |
return(sollya_lib_constant_from_int(int(anInt))) |
89 |
|
90 |
# Numeric Sollya function codes -> Sage mathematical function names
|
91 |
def pobyso_function_type_as_string(funcType): |
92 |
if funcType == SOLLYA_BASE_FUNC_ABS:
|
93 |
return "abs" |
94 |
elif funcType == SOLLYA_BASE_FUNC_ACOS:
|
95 |
return "arccos" |
96 |
elif funcType == SOLLYA_BASE_FUNC_ACOSH:
|
97 |
return "arccosh" |
98 |
elif funcType == SOLLYA_BASE_FUNC_ADD:
|
99 |
return "+" |
100 |
elif funcType == SOLLYA_BASE_FUNC_ASIN:
|
101 |
return "arcsin" |
102 |
elif funcType == SOLLYA_BASE_FUNC_ASINH:
|
103 |
return "arcsinh" |
104 |
elif funcType == SOLLYA_BASE_FUNC_ATAN:
|
105 |
return "arctan" |
106 |
elif funcType == SOLLYA_BASE_FUNC_ATANH:
|
107 |
return "arctanh" |
108 |
elif funcType == SOLLYA_BASE_FUNC_CEIL:
|
109 |
return "ceil" |
110 |
elif funcType == SOLLYA_BASE_FUNC_CONSTANT:
|
111 |
return "cte" |
112 |
elif funcType == SOLLYA_BASE_FUNC_COS:
|
113 |
return "cos" |
114 |
elif funcType == SOLLYA_BASE_FUNC_COSH:
|
115 |
return "cosh" |
116 |
elif funcType == SOLLYA_BASE_FUNC_DIV:
|
117 |
return "/" |
118 |
elif funcType == SOLLYA_BASE_FUNC_DOUBLE:
|
119 |
return "double" |
120 |
elif funcType == SOLLYA_BASE_FUNC_DOUBLEDOUBLE:
|
121 |
return "doubleDouble" |
122 |
elif funcType == SOLLYA_BASE_FUNC_DOUBLEEXTENDED:
|
123 |
return "doubleDxtended" |
124 |
elif funcType == SOLLYA_BASE_FUNC_ERF:
|
125 |
return "erf" |
126 |
elif funcType == SOLLYA_BASE_FUNC_ERFC:
|
127 |
return "erfc" |
128 |
elif funcType == SOLLYA_BASE_FUNC_EXP:
|
129 |
return "exp" |
130 |
elif funcType == SOLLYA_BASE_FUNC_EXP_M1:
|
131 |
return "expm1" |
132 |
elif funcType == SOLLYA_BASE_FUNC_FLOOR:
|
133 |
return "floor" |
134 |
elif funcType == SOLLYA_BASE_FUNC_FREE_VARIABLE:
|
135 |
return "freeVariable" |
136 |
elif funcType == SOLLYA_BASE_FUNC_HALFPRECISION:
|
137 |
return "halfPrecision" |
138 |
elif funcType == SOLLYA_BASE_FUNC_LIBRARYCONSTANT:
|
139 |
return "libraryConstant" |
140 |
elif funcType == SOLLYA_BASE_FUNC_LIBRARYFUNCTION:
|
141 |
return "libraryFunction" |
142 |
elif funcType == SOLLYA_BASE_FUNC_LOG:
|
143 |
return "log" |
144 |
elif funcType == SOLLYA_BASE_FUNC_LOG_10:
|
145 |
return "log10" |
146 |
elif funcType == SOLLYA_BASE_FUNC_LOG_1P:
|
147 |
return "log1p" |
148 |
elif funcType == SOLLYA_BASE_FUNC_LOG_2:
|
149 |
return "log2" |
150 |
elif funcType == SOLLYA_BASE_FUNC_MUL:
|
151 |
return "*" |
152 |
elif funcType == SOLLYA_BASE_FUNC_NEARESTINT:
|
153 |
return "round" |
154 |
elif funcType == SOLLYA_BASE_FUNC_NEG:
|
155 |
return "__neg__" |
156 |
elif funcType == SOLLYA_BASE_FUNC_PI:
|
157 |
return "pi" |
158 |
elif funcType == SOLLYA_BASE_FUNC_POW:
|
159 |
return "^" |
160 |
elif funcType == SOLLYA_BASE_FUNC_PROCEDUREFUNCTION:
|
161 |
return "procedureFunction" |
162 |
elif funcType == SOLLYA_BASE_FUNC_QUAD:
|
163 |
return "quad" |
164 |
elif funcType == SOLLYA_BASE_FUNC_SIN:
|
165 |
return "sin" |
166 |
elif funcType == SOLLYA_BASE_FUNC_SINGLE:
|
167 |
return "single" |
168 |
elif funcType == SOLLYA_BASE_FUNC_SINH:
|
169 |
return "sinh" |
170 |
elif funcType == SOLLYA_BASE_FUNC_SQRT:
|
171 |
return "sqrt" |
172 |
elif funcType == SOLLYA_BASE_FUNC_SUB:
|
173 |
return "-" |
174 |
elif funcType == SOLLYA_BASE_FUNC_TAN:
|
175 |
return "tan" |
176 |
elif funcType == SOLLYA_BASE_FUNC_TANH:
|
177 |
return "tanh" |
178 |
elif funcType == SOLLYA_BASE_FUNC_TRIPLEDOUBLE:
|
179 |
return "tripleDouble" |
180 |
else:
|
181 |
return None |
182 |
|
183 |
def pobyso_get_constant(rnArg, soConst): |
184 |
set_rn_value(rnArg, soConst) |
185 |
|
186 |
def pobyso_get_constant_as_rn(ctExp): |
187 |
precision = pobyso_get_prec_of_constant(ctExp) |
188 |
RRRR = RealField(precision) |
189 |
rn = RRRR(0)
|
190 |
sollya_lib_get_constant(get_rn_value(rn), ctExp) |
191 |
return(rn)
|
192 |
|
193 |
def pobyso_get_constant_as_rn_with_rf(ctExp, realField): |
194 |
rn = realField(0)
|
195 |
sollya_lib_get_constant(get_rn_value(rn), ctExp) |
196 |
return(rn)
|
197 |
def pobyso_get_free_variable_name(): |
198 |
return(sollya_lib_get_free_variable_name())
|
199 |
|
200 |
def pobyso_get_function_arity(expression): |
201 |
arity = c_int(0)
|
202 |
sollya_lib_get_function_arity(byref(arity),expression) |
203 |
return(int(arity.value)) |
204 |
|
205 |
def pobyso_get_head_function(expression): |
206 |
functionType = c_int(0)
|
207 |
sollya_lib_get_head_function(byref(functionType), expression, None)
|
208 |
return(int(functionType.value)) |
209 |
|
210 |
def pobyso_get_list_elements(soObj): |
211 |
# Type for array of pointers to sollya_obj_t
|
212 |
listAddress = POINTER(c_longlong)() |
213 |
numElements = c_int(0)
|
214 |
isEndElliptic = c_int(0)
|
215 |
listAsList = [] |
216 |
result = sollya_lib_get_list_elements(byref(listAddress),\ |
217 |
byref(numElements),\ |
218 |
byref(isEndElliptic),\ |
219 |
soObj) |
220 |
if result == 0 : |
221 |
return None |
222 |
for i in xrange(0, numElements.value, 1): |
223 |
print "address ", i, " ->", listAddress[i] |
224 |
listAsList.append(listAddress[i]) |
225 |
return(listAsList, numElements.value, isEndElliptic.value)
|
226 |
|
227 |
|
228 |
# Get the maximum precision used for the numbers in a
|
229 |
# Sollya expression.
|
230 |
# ToDo:
|
231 |
# - error management;
|
232 |
# - correctly deal with numerical type such as DOUBLEEXTENDED.
|
233 |
def pobyso_get_max_prec_of_exp(soExp): |
234 |
maxPrecision = 0
|
235 |
operator = pobyso_get_head_function(soExp) |
236 |
if (operator != SOLLYA_BASE_FUNC_CONSTANT) and \ |
237 |
(operator != SOLLYA_BASE_FUNC_FREE_VARIABLE): |
238 |
(arity, subexpressions) = pobyso_get_subfunctions(soExp) |
239 |
for i in xrange(arity): |
240 |
maxPrecisionCandidate = \ |
241 |
pobyso_get_max_prec_of_exp(subexpressions[i]) |
242 |
if maxPrecisionCandidate > maxPrecision:
|
243 |
maxPrecision = maxPrecisionCandidate |
244 |
return(maxPrecision)
|
245 |
elif operator == SOLLYA_BASE_FUNC_CONSTANT:
|
246 |
#print pobyso_get_prec_of_constant(soExp)
|
247 |
return(pobyso_get_prec_of_constant(soExp))
|
248 |
elif operator == SOLLYA_BASE_FUNC_FREE_VARIABLE:
|
249 |
return(0) |
250 |
else:
|
251 |
print "pobyso_get_max_prec_of_exp: unexepected operator." |
252 |
return(0) |
253 |
|
254 |
def pobyso_get_sage_exp_from_sollya_exp(sollyaExp, realField = RR): |
255 |
"""
|
256 |
Get a Sage expression from a Sollya expression, currently only tested
|
257 |
with polynomials with floating-point coefficients.
|
258 |
Notice that, in the returned polynomial, the exponents are RealNumbers.
|
259 |
"""
|
260 |
#pobyso_autoprint(sollyaExp)
|
261 |
operator = pobyso_get_head_function(sollyaExp) |
262 |
# Constants and the free variable are special cases.
|
263 |
# All other operator are dealt with in the same way.
|
264 |
if (operator != SOLLYA_BASE_FUNC_CONSTANT) and \ |
265 |
(operator != SOLLYA_BASE_FUNC_FREE_VARIABLE): |
266 |
(arity, subexpressions) = pobyso_get_subfunctions(sollyaExp) |
267 |
if arity == 1: |
268 |
sageExp = eval(pobyso_function_type_as_string(operator) + \
|
269 |
"(" + pobyso_get_sage_exp_from_sollya_exp(subexpressions[0], realField)\ |
270 |
+ ")")
|
271 |
elif arity == 2: |
272 |
if operator == SOLLYA_BASE_FUNC_POW:
|
273 |
operatorAsString = "**"
|
274 |
else:
|
275 |
operatorAsString = pobyso_function_type_as_string(operator) |
276 |
sageExp = \ |
277 |
eval("pobyso_get_sage_exp_from_sollya_exp(subexpressions[0], realField)"\ |
278 |
+ " " + operatorAsString + " " + \ |
279 |
"pobyso_get_sage_exp_from_sollya_exp(subexpressions[1], realField)")
|
280 |
# We do not know yet how to deal with arity > 3 (is there any in Sollya anyway?).
|
281 |
else:
|
282 |
sageExp = eval('None') |
283 |
return(sageExp)
|
284 |
elif operator == SOLLYA_BASE_FUNC_CONSTANT:
|
285 |
#print "This is a constant"
|
286 |
return pobyso_get_constant_as_rn_with_rf(sollyaExp, realField)
|
287 |
elif operator == SOLLYA_BASE_FUNC_FREE_VARIABLE:
|
288 |
#print "This is free variable"
|
289 |
return(eval(sollya_lib_get_free_variable_name())) |
290 |
else:
|
291 |
print "Unexpected" |
292 |
return eval('None') |
293 |
# End pobyso_get_sage_poly_from_sollya_poly
|
294 |
|
295 |
def pobyso_get_subfunctions(expression): |
296 |
subf0 = c_int(0)
|
297 |
subf1 = c_int(0)
|
298 |
subf2 = c_int(0)
|
299 |
subf3 = c_int(0)
|
300 |
subf4 = c_int(0)
|
301 |
subf5 = c_int(0)
|
302 |
subf6 = c_int(0)
|
303 |
subf7 = c_int(0)
|
304 |
subf8 = c_int(0)
|
305 |
arity = c_int(0)
|
306 |
nullPtr = POINTER(c_int)() |
307 |
sollya_lib_get_subfunctions(expression, byref(arity), \ |
308 |
byref(subf0), byref(subf1), byref(subf2), byref(subf3), byref(subf4), byref(subf5),\ |
309 |
byref(subf6), byref(subf7), byref(subf8), nullPtr, None)
|
310 |
# byref(cast(subfunctions[0], POINTER(c_int))), byref(cast(subfunctions[0], POINTER(c_int))), \
|
311 |
# byref(cast(subfunctions[2], POINTER(c_int))), byref(cast(subfunctions[3], POINTER(c_int))), \
|
312 |
# byref(cast(subfunctions[4], POINTER(c_int))), byref(cast(subfunctions[5], POINTER(c_int))), \
|
313 |
# byref(cast(subfunctions[6], POINTER(c_int))), byref(cast(subfunctions[7], POINTER(c_int))), \
|
314 |
# byref(cast(subfunctions[8], POINTER(c_int))), nullPtr)
|
315 |
subfunctions = [subf0, subf1, subf2, subf3, subf4, subf5, subf6, subf7, subf8] |
316 |
subs = [] |
317 |
if arity.value > pobyso_max_arity:
|
318 |
return(None,None) |
319 |
for i in xrange(arity.value): |
320 |
subs.append(int(subfunctions[i].value))
|
321 |
#print subs[i]
|
322 |
return(int(arity.value), subs) |
323 |
|
324 |
def pobyso_get_prec(): |
325 |
retc = sollya_lib_get_prec(None)
|
326 |
a = c_int(0)
|
327 |
sollya_lib_get_constant_as_int(byref(a), retc) |
328 |
return(int(a.value)) |
329 |
|
330 |
def pobyso_get_prec_of_constant(ctExp): |
331 |
prec = c_int(0)
|
332 |
retc = sollya_lib_get_prec_of_constant(byref(prec), ctExp, None)
|
333 |
return(int(prec.value)) |
334 |
|
335 |
def pobyso_parse_string(string): |
336 |
return(sollya_lib_parse_string(string))
|
337 |
|
338 |
def pobyso_univar_polynomial_print_reverse(polySa): |
339 |
"""
|
340 |
Return the string representation of a univariate polynomial with
|
341 |
monomial ordered in the x^0..x^n order of the monomials.
|
342 |
Remember: Sage
|
343 |
"""
|
344 |
polynomialRing = polySa.base_ring() |
345 |
# A very expensive solution:
|
346 |
# -create a fake multivariate polynomial field with only one variable,
|
347 |
# specifying a negative lexicographical order;
|
348 |
mpolynomialRing = PolynomialRing(polynomialRing.base(), \ |
349 |
polynomialRing.variable_name(), \ |
350 |
1, order='neglex') |
351 |
# - convert the univariate argument polynomial into a multivariate
|
352 |
# version;
|
353 |
p = mpolynomialRing(polySa) |
354 |
# - return the string representation of the converted form.
|
355 |
# There is no simple str() method defined for p's class.
|
356 |
return(p.__str__())
|
357 |
|
358 |
def pobyso_range(rnLowerBound, rnUpperBound): |
359 |
lowerBoundSo = sollya_lib_constant(get_rn_value(rnLowerBound)) |
360 |
upperBoundSo = sollya_lib_constant(get_rn_value(rnUpperBound)) |
361 |
rangeSo = sollya_lib_range(lowerBoundSo, upperBoundSo) |
362 |
return(rangeSo)
|
363 |
|
364 |
def pobyso_remez_canonical(function, \ |
365 |
degree, \ |
366 |
lowerBound, \ |
367 |
upperBound, \ |
368 |
weightSo = pobyso_constant_1(), |
369 |
quality = None):
|
370 |
if parent(function) == parent("string"): |
371 |
functionSo = sollya_lib_parse_string(function) |
372 |
# print "Is string!"
|
373 |
elif sollya_lib_obj_is_function(function):
|
374 |
functionSo = function |
375 |
# print "Is Function!"
|
376 |
degreeSo = pobyso_constant_from_int(degree) |
377 |
rangeSo = pobyso_range(lowerBound, upperBound) |
378 |
return(sollya_lib_remez(functionSo, degreeSo, rangeSo, quality, None)) |
379 |
|
380 |
def pobyso_set_canonical_off(): |
381 |
sollya_lib_set_canonical(sollya_lib_off()) |
382 |
|
383 |
def pobyso_set_canonical_on(): |
384 |
sollya_lib_set_canonical(sollya_lib_on()) |
385 |
|
386 |
def pobyso_set_prec(p): |
387 |
a = c_int(p) |
388 |
precSo = c_void_p(sollya_lib_constant_from_int(a)) |
389 |
sollya_lib_set_prec(precSo) |
390 |
|
391 |
def pobyso_taylor(function, degree, point): |
392 |
return(sollya_lib_taylor(function, degree, point))
|
393 |
|
394 |
def pobyso_taylorform(function, degree, point = None, interval = None, errorType=None): |
395 |
if errorType is None: |
396 |
errorType = sollya_lib_absolute() |
397 |
return(sollya_lib_taylorform(function, degree, point, errorType, None)) |
398 |
#
|
399 |
print "Superficial test of pobyso:" |
400 |
print pobyso_get_prec()
|
401 |
pobyso_set_prec(165)
|
402 |
print pobyso_get_prec()
|
403 |
a=100
|
404 |
print type(a) |
405 |
id(a)
|
406 |
print "Max arity: ", pobyso_max_arity |
407 |
print "Function tripleDouble (43) as a string: ", pobyso_function_type_as_string(43) |
408 |
print "Function None (44) as a string: ", pobyso_function_type_as_string(44) |