Statistiques
| Révision :

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)