Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 59

Historique | Voir | Annoter | Télécharger (36,77 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
pobyso functions come in five flavors: 
9
- the _so_so (arguments and returned objects are pointers to Sollya objects, 
10
  includes the void function and the no arguments function that return a 
11
  pointer to a Sollya object);
12
- the _so_sa (argument are pointers to Sollya objects, returned objects are
13
  Sage/Python objects or, more generally, information is transfered from the
14
  Sollya world to Sage/Python world; e.g. functions without arguments that
15
  return a Sage/Python object);
16
- the _sa_so (arguments are Sage/Python objects, returned objects are 
17
  pointers to Sollya objects);
18
- the sa_sa (arguments and returned objects are all Sage/Python objects);
19
- a catch all flavor, without any suffix, (e. g. functions that have no argument 
20
  nor return value).
21
This classification is not always very strict. Conversion functions from Sollya
22
to Sage/Python are sometimes decorated with Sage/Python arguments to set
23
the precision. These functions remain in the so_sa category.
24
NOTES:
25
Reported errors in Eclipse come from the calls to
26
the Sollya library
27

28
ToDo (among other things): 
29
 -memory management.
30
"""
31
from ctypes import *
32
import re
33
from sage.symbolic.expression_conversions import polynomial
34
from sage.symbolic.expression_conversions import PolynomialConverter
35
"""
36
Create the equivalent to an enum for the Sollya function types.
37
"""
38
(SOLLYA_BASE_FUNC_ABS,
39
SOLLYA_BASE_FUNC_ACOS,
40
    SOLLYA_BASE_FUNC_ACOSH,
41
    SOLLYA_BASE_FUNC_ADD,
42
    SOLLYA_BASE_FUNC_ASIN,
43
    SOLLYA_BASE_FUNC_ASINH,
44
    SOLLYA_BASE_FUNC_ATAN,
45
    SOLLYA_BASE_FUNC_ATANH,
46
    SOLLYA_BASE_FUNC_CEIL,
47
    SOLLYA_BASE_FUNC_CONSTANT,
48
    SOLLYA_BASE_FUNC_COS,
49
    SOLLYA_BASE_FUNC_COSH,
50
    SOLLYA_BASE_FUNC_DIV,
51
    SOLLYA_BASE_FUNC_DOUBLE,
52
    SOLLYA_BASE_FUNC_DOUBLEDOUBLE,
53
    SOLLYA_BASE_FUNC_DOUBLEEXTENDED,
54
    SOLLYA_BASE_FUNC_ERF,
55
    SOLLYA_BASE_FUNC_ERFC,
56
    SOLLYA_BASE_FUNC_EXP,
57
    SOLLYA_BASE_FUNC_EXP_M1,
58
    SOLLYA_BASE_FUNC_FLOOR,
59
    SOLLYA_BASE_FUNC_FREE_VARIABLE,
60
    SOLLYA_BASE_FUNC_HALFPRECISION,
61
    SOLLYA_BASE_FUNC_LIBRARYCONSTANT,
62
    SOLLYA_BASE_FUNC_LIBRARYFUNCTION,
63
    SOLLYA_BASE_FUNC_LOG,
64
    SOLLYA_BASE_FUNC_LOG_10,
65
    SOLLYA_BASE_FUNC_LOG_1P,
66
    SOLLYA_BASE_FUNC_LOG_2,
67
    SOLLYA_BASE_FUNC_MUL,
68
    SOLLYA_BASE_FUNC_NEARESTINT,
69
    SOLLYA_BASE_FUNC_NEG,
70
    SOLLYA_BASE_FUNC_PI,
71
    SOLLYA_BASE_FUNC_POW,
72
    SOLLYA_BASE_FUNC_PROCEDUREFUNCTION,
73
    SOLLYA_BASE_FUNC_QUAD,
74
    SOLLYA_BASE_FUNC_SIN,
75
    SOLLYA_BASE_FUNC_SINGLE,
76
    SOLLYA_BASE_FUNC_SINH,
77
    SOLLYA_BASE_FUNC_SQRT,
78
    SOLLYA_BASE_FUNC_SUB,
79
    SOLLYA_BASE_FUNC_TAN,
80
    SOLLYA_BASE_FUNC_TANH,
81
SOLLYA_BASE_FUNC_TRIPLEDOUBLE) = map(int,xrange(44))
82
print "\nSuperficial pobyso check..."    
83
print "First constant - SOLLYA_BASE_FUNC_ABS: ", SOLLYA_BASE_FUNC_ABS
84
print "Last constant  - SOLLYA_BASE_FUNC_TRIPLEDOUBLE: ", SOLLYA_BASE_FUNC_TRIPLEDOUBLE
85

    
86
pobyso_max_arity = 9
87

    
88
def pobyso_autoprint(arg):
89
    sollya_lib_autoprint(arg,None)
90

    
91
def pobyso_autoprint_so_so(arg):
92
    sollya_lib_autoprint(arg,None)
93
    
94
def pobyso_absolute_so_so():
95
    return(sollya_lib_absolute(None))
96

    
97
def pobyso_build_function_sub_so_so(exp1So, exp2So):
98
    return(sollya_lib_build_function_sub(exp1So, exp2So))
99

    
100
def pobyso_cmp(rnArg, soCte):
101
    """
102
    Compare the MPFR value a RealNumber with that of a Sollya constant.
103
    
104
    Get the value of the Sollya constant into a RealNumber and compare
105
    using MPFR. Could be optimized by working directly with a mpfr_t
106
    for the intermediate number. 
107
    """
108
    precisionOfCte = c_int(0)
109
    # From the Sollya constant, create a local Sage RealNumber.
110
    sollya_lib_get_prec_of_constant(precisionOfCte, soCte) 
111
    #print "Precision of constant: ", precisionOfCte
112
    RRRR = RealField(precisionOfCte.value)
113
    rnLocal = RRRR(0)
114
    sollya_lib_get_constant(get_rn_value(rnLocal), soCte)
115
    #print "rnDummy: ", rnDummy
116
    # Compare the local Sage RealNumber with rnArg.
117
    return(cmp_rn_value(rnArg, rnLocal))
118

    
119
def pobyso_change_var_in_function_so_so(funcSo, chvarExpSo):
120
    return(sollya_lib_evaluate(funcSo,chvarExpSo))
121
    
122

    
123
def pobyso_chebyshevform_so_so(functionSo, degreeSo, intervalSo):
124
    resultSo = sollya_lib_chebyshevform(functionSo, degreeSo, intervalSo)
125
    return(resultSo)
126
        
127

    
128

    
129
def pobyso_compute_pos_function_abs_val_bounds_sa_sa(funcSa, lowerBoundSa, \
130
                                                     upperBoundSa):
131
    """
132
    TODO: set the variable name in Sollya.
133
    """
134
    funcSo = pobyso_parse_string(funcSa._assume_str())
135
    rangeSo = pobyso_range_sa_so(lowerBoundSa, upperBoundSa)
136
    infnormSo = pobyso_infnorm_so_so(funcSo,rangeSo)
137
    fMaxSa = pobyso_get_interval_from_range_so_sa(infnormSo)
138
    # Get the top bound and compute the binade top limit.
139
    fMaxUpperBoundSa = fMaxSa.upper()
140
    binadeTopLimitSa = 2**ceil(fMaxUpperBoundSa.log2())
141
    # Put up together the function to use to compute the lower bound.
142
    funcAuxSo = pobyso_parse_string(str(binadeTopLimitSa) +  \
143
                                    '-(' + f._assume_str() + ')')
144
    pobyso_autoprint(funcAuxSo)
145
    # Clear the Sollay range before a new call to infnorm and issue the call.
146
    sollya_lib_clear_obj(infnormSo)
147
    infnormSo = pobyso_infnorm_so_so(funcAuxSo,rangeSo)
148
    fMinSa = pobyso_get_interval_from_range_so_sa(infnormSo)
149
    sollya_lib_clear_obj(infnormSo)
150
    fMinLowerBoundSa = topBinadeLimit - fMinSa.lower()
151
    # Compute the maximum of the precisions of the different bounds.
152
    maxPrecSa = max([fMinLowerBoundSa.parent().precision(), \
153
                     fMaxUpperBoundSa.parent().precision()])
154
    # Create a RealIntervalField and create an interval with the "good" bounds.
155
    RRRI = RealIntervalField(maxPrecSa)
156
    imageIntervalSa = RRRI(fMinLowerBoundSa, fMaxUpperBoundSa)
157
    # Free the uneeded Sollya objects
158
    sollya_lib_clear_obj(funcSo)
159
    sollya_lib_clear_obj(funcAuxSo)
160
    sollya_lib_clear_obj(rangeSo)
161
    return(imageIntervalSa)
162
    # End pobyso_compute_function_abs_val_bounds_sa_sa
163

    
164
def pobyso_constant(rnArg):
165
    """ Legacy function. See pobyso_constant_sa_so. """
166
    return(pobyso_constant_sa_so(rnArg))
167
    
168
def pobyso_constant_sa_so(rnArg):
169
    """
170
    Create a Sollya constant from a RealNumber.
171
    """
172
    return(sollya_lib_constant(get_rn_value(rnArg)))
173
    
174
def pobyso_constant_0_sa_so():
175
    return(pobyso_constant_from_int_sa_so(0))
176

    
177
def pobyso_constant_1():
178
    """ Legacy function. See pobyso_constant_so_so. """
179
    return(pobyso_constant_1_sa_so())
180

    
181
def pobyso_constant_1_sa_so():
182
    return(pobyso_constant_from_int_sa_so(1))
183

    
184
def pobyso_constant_from_int(anInt):
185
    """ Legacy function. See pobyso_constant_from_int_sa_so. """
186
    return(pobyso_constant_from_int_sa_so(anInt))
187

    
188
def pobyso_constant_from_int_sa_so(anInt):
189
    return(sollya_lib_constant_from_int(int(anInt)))
190

    
191
def pobyso_function_type_as_string(funcType):
192
    """ Legacy function. See pobyso_function_type_as_string_so_sa. """
193
    return(pobyso_function_type_as_string_so_sa(funcType))
194

    
195
def pobyso_function_type_as_string_so_sa(funcType):
196
    """
197
    Numeric Sollya function codes -> Sage mathematical function names.
198
    Notice that pow -> ^ (a la Sage, not a la Python).
199
    """
200
    if funcType == SOLLYA_BASE_FUNC_ABS:
201
        return "abs"
202
    elif funcType == SOLLYA_BASE_FUNC_ACOS:
203
        return "arccos"
204
    elif funcType == SOLLYA_BASE_FUNC_ACOSH:
205
        return "arccosh"
206
    elif funcType == SOLLYA_BASE_FUNC_ADD:
207
        return "+"
208
    elif funcType == SOLLYA_BASE_FUNC_ASIN:
209
        return "arcsin"
210
    elif funcType == SOLLYA_BASE_FUNC_ASINH:
211
        return "arcsinh"
212
    elif funcType == SOLLYA_BASE_FUNC_ATAN:
213
        return "arctan"
214
    elif funcType == SOLLYA_BASE_FUNC_ATANH:
215
        return "arctanh"
216
    elif funcType == SOLLYA_BASE_FUNC_CEIL:
217
        return "ceil"
218
    elif funcType == SOLLYA_BASE_FUNC_CONSTANT:
219
        return "cte"
220
    elif funcType == SOLLYA_BASE_FUNC_COS:
221
        return "cos"
222
    elif funcType == SOLLYA_BASE_FUNC_COSH:
223
        return "cosh"
224
    elif funcType == SOLLYA_BASE_FUNC_DIV:
225
        return "/"
226
    elif funcType == SOLLYA_BASE_FUNC_DOUBLE:
227
        return "double"
228
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEDOUBLE:
229
        return "doubleDouble"
230
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEEXTENDED:
231
        return "doubleDxtended"
232
    elif funcType == SOLLYA_BASE_FUNC_ERF:
233
        return "erf"
234
    elif funcType == SOLLYA_BASE_FUNC_ERFC:
235
        return "erfc"
236
    elif funcType == SOLLYA_BASE_FUNC_EXP:
237
        return "exp"
238
    elif funcType == SOLLYA_BASE_FUNC_EXP_M1:
239
        return "expm1"
240
    elif funcType == SOLLYA_BASE_FUNC_FLOOR:
241
        return "floor"
242
    elif funcType == SOLLYA_BASE_FUNC_FREE_VARIABLE:
243
        return "freeVariable"
244
    elif funcType == SOLLYA_BASE_FUNC_HALFPRECISION:
245
        return "halfPrecision"
246
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYCONSTANT:
247
        return "libraryConstant"
248
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYFUNCTION:
249
        return "libraryFunction"
250
    elif funcType == SOLLYA_BASE_FUNC_LOG:
251
        return "log"
252
    elif funcType == SOLLYA_BASE_FUNC_LOG_10:
253
        return "log10"
254
    elif funcType == SOLLYA_BASE_FUNC_LOG_1P:
255
        return "log1p"
256
    elif funcType == SOLLYA_BASE_FUNC_LOG_2:
257
        return "log2"
258
    elif funcType == SOLLYA_BASE_FUNC_MUL:
259
        return "*"
260
    elif funcType == SOLLYA_BASE_FUNC_NEARESTINT:
261
        return "round"
262
    elif funcType == SOLLYA_BASE_FUNC_NEG:
263
        return "__neg__"
264
    elif funcType == SOLLYA_BASE_FUNC_PI:
265
        return "pi"
266
    elif funcType == SOLLYA_BASE_FUNC_POW:
267
        return "^"
268
    elif funcType == SOLLYA_BASE_FUNC_PROCEDUREFUNCTION:
269
        return "procedureFunction"
270
    elif funcType == SOLLYA_BASE_FUNC_QUAD:
271
        return "quad"
272
    elif funcType == SOLLYA_BASE_FUNC_SIN:
273
        return "sin"
274
    elif funcType == SOLLYA_BASE_FUNC_SINGLE:
275
        return "single"
276
    elif funcType == SOLLYA_BASE_FUNC_SINH:
277
        return "sinh"
278
    elif funcType == SOLLYA_BASE_FUNC_SQRT:
279
        return "sqrt"
280
    elif funcType == SOLLYA_BASE_FUNC_SUB:
281
        return "-"
282
    elif funcType == SOLLYA_BASE_FUNC_TAN:
283
        return "tan"
284
    elif funcType == SOLLYA_BASE_FUNC_TANH:
285
        return "tanh"
286
    elif funcType == SOLLYA_BASE_FUNC_TRIPLEDOUBLE:
287
        return "tripleDouble"
288
    else:
289
        return None
290

    
291
def pobyso_get_constant(rnArg, soConst):
292
    """ Legacy function. See pobyso_get_constant_so_sa. """
293
    return(pobyso_get_constant_so_sa(rnArg, soConst))
294

    
295
def pobyso_get_constant_so_sa(rnArg, soConst):
296
    """
297
    Set the value of rnArg to the value of soConst in MPFR_RNDN mode.
298
    rnArg must already exist and belong to some RealField.
299
    We assume that soConst points to a Sollya constant.
300
    """
301
    return(sollya_lib_get_constant(get_rn_value(rnArg), soConst))
302
    
303
def pobyso_get_constant_as_rn(ctExpSo):
304
    """ Legacy function. See pobyso_get_constant_as_rn_so_sa. """ 
305
    return(pobyso_get_constant_as_rn_so_sa(ctExpSo))
306
    
307
def pobyso_get_constant_as_rn_so_sa(constExpSo):
308
    precisionSa  = pobyso_get_prec_of_constant_so_sa(constExpSo) 
309
    RRRR = RealField(precisionSa)
310
    rnSa = RRRR(0)
311
    sollya_lib_get_constant(get_rn_value(rnSa), constExpSo)
312
    return(rnSa)
313

    
314
def pobyso_get_constant_as_rn_with_rf(ctExp, realField):
315
    """ Legacy function. See pobyso_get_constant_as_rn_with_rf_so_sa."""
316
    return(pobyso_get_constant_as_rn_with_rf_so_sa(ctExp, realField))
317
    
318
def pobyso_get_constant_as_rn_with_rf_so_sa(ctExpSo, realFieldSa = None):
319
    if realFieldSa is None:
320
        sollyaPrecSa = pobyso_get_prec_so_sa()
321
        realFieldSa = RealField(sollyaPrecSa)
322
    rnSa = realFieldSa(0)
323
    sollya_lib_get_constant(get_rn_value(rnSa), ctExpSo)
324
    return(rnSa)
325

    
326
def pobyso_get_free_variable_name():
327
    """ Legacy function. See pobyso_get_free_variable_name_so_sa."""
328
    return(pobyso_get_free_variable_name_so_sa())
329

    
330
def pobyso_get_free_variable_name_so_sa():
331
    return(sollya_lib_get_free_variable_name())
332
    
333
def pobyso_get_function_arity(expressionSo):
334
    """ Legacy function. See pobyso_get_function_arity_so_sa."""
335
    return(pobyso_get_function_arity_so_sa(expressionSo))
336

    
337
def pobyso_get_function_arity_so_sa(expressionSo):
338
    arity = c_int(0)
339
    sollya_lib_get_function_arity(byref(arity),expressionSo)
340
    return(int(arity.value))
341

    
342
def pobyso_get_head_function(expressionSo):
343
    """ Legacy function. See pobyso_get_head_function_so_sa. """
344
    return(pobyso_get_head_function_so_sa(expressionSo)) 
345

    
346
def pobyso_get_head_function_so_sa(expressionSo):
347
    functionType = c_int(0)
348
    sollya_lib_get_head_function(byref(functionType), expressionSo, None)
349
    return(int(functionType.value))
350

    
351
def pobyso_get_interval_from_range_so_sa(soRange, realIntervalFieldSa = None ):
352
    """
353
    Return the Sage interval corresponding to the Sollya range argument.
354
    If no reaInterval lField is passed as argument, the interval bounds are not
355
    rounded: they are elements of RealIntervalField of the "right" precision
356
    to hold all the digits.
357
    """
358
    prec = c_int(0)
359
    if realIntervalFieldSa is None:
360
        retval = sollya_lib_get_prec_of_range(byref(prec), soRange, None)
361
        if retval == 0:
362
            return(None)
363
        realIntervalFieldSa = RealIntervalField(prec.value)
364
    intervalSa = realIntervalFieldSa(0,0)
365
    retval = \
366
        sollya_lib_get_interval_from_range(get_interval_value(intervalSa),\
367
                                           soRange)
368
    if retval == 0:
369
        return(None)
370
    return(intervalSa)
371
# End pobyso_get_interval_from_range_so_sa
372

    
373
def pobyso_get_list_elements(soObj):
374
    """ Legacy function. See pobyso_get_list_elements_so_so. """
375
    return(pobyso_get_list_elements_so_so(soObj))
376
 
377
def pobyso_get_list_elements_so_so(soObj):
378
    """
379
    Get the list elements as a Sage/Python array of Sollya objects.
380
    The other data returned are also Sage/Python objects.
381
    """
382
    listAddress = POINTER(c_longlong)()
383
    numElements = c_int(0)
384
    isEndElliptic = c_int(0)
385
    listAsList = []
386
    result = sollya_lib_get_list_elements(byref(listAddress),\
387
                                          byref(numElements),\
388
                                          byref(isEndElliptic),\
389
                                          soObj)
390
    if result == 0 :
391
        return None
392
    for i in xrange(0, numElements.value, 1):
393
        listAsList.append(listAddress[i])
394
    return(listAsList, numElements.value, isEndElliptic.value)
395

    
396
def pobyso_get_max_prec_of_exp(soExp):
397
    """ Legacy function. See pobyso_get_max_prec_of_exp_so_sa. """
398
    return(pobyso_get_max_prec_of_exp_so_sa(soExp))
399

    
400
def pobyso_get_max_prec_of_exp_so_sa(soExp):
401
    """
402
    Get the maximum precision used for the numbers in a Sollya expression.
403
    
404
    Arguments:
405
    soExp -- a Sollya expression pointer
406
    Return value:
407
    A Python integer
408
    TODO: 
409
    - error management;
410
    - correctly deal with numerical type such as DOUBLEEXTENDED.
411
    """
412
    maxPrecision = 0
413
    minConstPrec = 0
414
    currentConstPrec = 0
415
    operator = pobyso_get_head_function_so_sa(soExp)
416
    if (operator != SOLLYA_BASE_FUNC_CONSTANT) and \
417
    (operator != SOLLYA_BASE_FUNC_FREE_VARIABLE):
418
        (arity, subexpressions) = pobyso_get_subfunctions_so_sa(soExp)
419
        for i in xrange(arity):
420
            maxPrecisionCandidate = \
421
                pobyso_get_max_prec_of_exp_so_sa(subexpressions[i])
422
            if maxPrecisionCandidate > maxPrecision:
423
                maxPrecision = maxPrecisionCandidate
424
        return(maxPrecision)
425
    elif operator == SOLLYA_BASE_FUNC_CONSTANT:
426
        minConstPrec = pobyso_get_min_prec_of_constant_so_sa(soExp)
427
        #currentConstPrec = pobyso_get_min_prec_of_constant_so_sa(soExp)
428
        #print minConstPrec, " - ", currentConstPrec 
429
        return(pobyso_get_min_prec_of_constant_so_sa(soExp))
430
    
431
    elif operator == SOLLYA_BASE_FUNC_FREE_VARIABLE:
432
        return(0)
433
    else:
434
        print "pobyso_get_max_prec_of_exp_so_sa: unexepected operator."
435
        return(0)
436

    
437
def pobyso_get_min_prec_of_constant_so_sa(soConstExp):
438
    """
439
    Get the minimum precision necessary to represent the value of a Sollya
440
    constant.
441
    MPFR_MIN_PREC and powers of 2 are taken into account.
442
    We assume that soCteExp is a point
443
    """
444
    constExpAsRn = pobyso_get_constant_as_rn_so_sa(soConstExp)
445
    return(min_mpfr_size(get_rn_value(constExpAsRn)))
446

    
447
def pobyso_get_sage_exp_from_sollya_exp(sollyaExp, realField = RR):
448
    """ Legacy function. See pobyso_get_sage_exp_from_sollya_exp_so_sa. """
449
    return(pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExp, realField = RR))
450

    
451
def pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExp, realField = RR):
452
    """
453
    Get a Sage expression from a Sollya expression. 
454
    Currently only tested with polynomials with floating-point coefficients.
455
    Notice that, in the returned polynomial, the exponents are RealNumbers.
456
    """
457
    #pobyso_autoprint(sollyaExp)
458
    operator = pobyso_get_head_function_so_sa(sollyaExp)
459
    # Constants and the free variable are special cases.
460
    # All other operator are dealt with in the same way.
461
    if (operator != SOLLYA_BASE_FUNC_CONSTANT) and \
462
       (operator != SOLLYA_BASE_FUNC_FREE_VARIABLE):
463
        (arity, subexpressions) = pobyso_get_subfunctions_so_sa(sollyaExp)
464
        if arity == 1:
465
            sageExp = eval(pobyso_function_type_as_string_so_sa(operator) + \
466
            "(" + pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressions[0], \
467
            realField) + ")")
468
        elif arity == 2:
469
            if operator == SOLLYA_BASE_FUNC_POW:
470
                operatorAsString = "**"
471
            else:
472
                operatorAsString = \
473
                    pobyso_function_type_as_string_so_sa(operator)
474
            sageExp = \
475
              eval("pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressions[0], realField)"\
476
              + " " + operatorAsString + " " + \
477
                   "pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressions[1], realField)")
478
        # We do not know yet how to deal with arity > 3 (is there any in Sollya anyway?).
479
        else:
480
            sageExp = eval('None')
481
        return(sageExp)
482
    elif operator == SOLLYA_BASE_FUNC_CONSTANT:
483
        #print "This is a constant"
484
        return pobyso_get_constant_as_rn_with_rf_so_sa(sollyaExp, realField)
485
    elif operator == SOLLYA_BASE_FUNC_FREE_VARIABLE:
486
        #print "This is free variable"
487
        return(eval(sollya_lib_get_free_variable_name()))
488
    else:
489
        print "Unexpected"
490
        return eval('None')
491
# End pobyso_get_sage_poly_from_sollya_poly
492
def pobyso_get_sage_poly_from_sollya_poly(polySo, realFieldSa=None):
493
    """
494
    Convert a Sollya polynomial into a Sage polynomial.
495
    If no realField is given, a RealField corresponding to the maximum precision 
496
    of the coefficients is internally computed.
497
    It is not returned but can be easily retrieved from the polynomial itself.
498
    Main steps:
499
    - (optional) compute the RealField of the coefficients;
500
    - convert the Sollya expression into a Sage expression;
501
    - convert the Sage expression into a Sage polynomial
502
    TODO: the canonical thing for the polynomial.
503
    """    
504
    if realFieldSa is None:
505
        expressionPrecSa = pobyso_get_max_prec_of_exp_so_sa(polySo)
506
        realFieldSa = RealField(expressionPrecSa)
507
    
508
    expressionSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, \
509
                                                             realFieldSa)
510
    print "ssssss: ", expressionSa
511
    polyVariableSa = expressionSa.variables()[0]
512
    print type(str(polyVariableSa))
513
    polyRingSa = realFieldSa[str(polyVariableSa)]
514
    print polyRingSa
515
    polynomialSa = polynomial(polyVariableSa, ring=polyRingSa)
516
    print polyRingSa.base()
517
    print "tttttt: ", polynomialSa
518
    return(polynomialSa)
519
# End pobyso_get_sage_poly_from_sollya_poly
520

    
521
def pobyso_get_subfunctions(expressionSo):
522
    """ Legacy function. See pobyso_get_subfunctions_so_sa. """
523
    return(pobyso_get_subfunctions_so_sa(expressionSo)) 
524

    
525
def pobyso_get_subfunctions_so_sa(expressionSo):
526
    """
527
    Get the subfunctions of an expression.
528
    Return the number of subfunctions and the list of subfunctions addresses.
529
    S.T.: Could not figure out another way than that ugly list of declarations
530
    to recover the addresses of the subfunctions. 
531
    """
532
    subf0 = c_int(0)
533
    subf1 = c_int(0)
534
    subf2 = c_int(0)
535
    subf3 = c_int(0)
536
    subf4 = c_int(0)
537
    subf5 = c_int(0)
538
    subf6 = c_int(0)
539
    subf7 = c_int(0)
540
    subf8 = c_int(0)
541
    arity = c_int(0)
542
    nullPtr = POINTER(c_int)()
543
    sollya_lib_get_subfunctions(expressionSo, byref(arity), \
544
    byref(subf0), byref(subf1), byref(subf2), byref(subf3), byref(subf4), byref(subf5),\
545
     byref(subf6), byref(subf7), byref(subf8), nullPtr, None) 
546
#    byref(cast(subfunctions[0], POINTER(c_int))), byref(cast(subfunctions[0], POINTER(c_int))), \
547
#    byref(cast(subfunctions[2], POINTER(c_int))), byref(cast(subfunctions[3], POINTER(c_int))), \
548
#    byref(cast(subfunctions[4], POINTER(c_int))), byref(cast(subfunctions[5], POINTER(c_int))), \
549
#    byref(cast(subfunctions[6], POINTER(c_int))), byref(cast(subfunctions[7], POINTER(c_int))), \
550
#    byref(cast(subfunctions[8], POINTER(c_int))), nullPtr)
551
    subfunctions = [subf0, subf1, subf2, subf3, subf4, subf5, subf6, subf7, subf8]
552
    subs = []
553
    if arity.value > pobyso_max_arity:
554
        return(0,[])
555
    for i in xrange(arity.value):
556
        subs.append(int(subfunctions[i].value))
557
        #print subs[i]
558
    return(int(arity.value), subs)
559
    
560
def pobyso_get_prec():
561
    """ Legacy function. See pobyso_get_prec_so_sa(). """
562
    return(pobyso_get_prec_so_sa())
563

    
564
def pobyso_get_prec_so_sa():
565
    """
566
    Get the current default precision in Sollya.
567
    The return value is Sage/Python int.
568
    """
569
    precSo = sollya_lib_get_prec(None)
570
    precSa = c_int(0)
571
    sollya_lib_get_constant_as_int(byref(precSa), precSo)
572
    sollya_lib_clear_obj(precSo)
573
    return(int(precSa.value))
574

    
575
def pobyso_get_prec_of_constant(ctExpSo):
576
    """ Legacy function. See pobyso_get_prec_of_constant_so_sa. """
577
    return(pobyso_get_prec_of_constant_so_sa(ctExpSo))
578

    
579
def pobyso_get_prec_of_constant_so_sa(ctExpSo):
580
    prec = c_int(0)
581
    retc = sollya_lib_get_prec_of_constant(byref(prec), ctExpSo, None)
582
    if retc == 0:
583
        return(None)
584
    return(int(prec.value))
585

    
586
def pobyso_get_prec_of_range_so_sa(rangeSo):
587
    prec = c_int(0)
588
    retc = sollya_lib_get_prec_of_range(byref(prec), rangeSo, None)
589
    if retc == 0:
590
        return(None)
591
    return(int(prec.value))
592

    
593
def pobyso_infnorm_so_so(func, interval, file = None, intervalList = None):
594
    print "Do not use this function. User pobyso_supnorm_so_so instead."
595
    return(None)
596

    
597
def pobyso_lib_init():
598
    sollya_lib_init(None)
599
    
600
def pobyso_name_free_variable(freeVariableName):
601
    """ Legacy function. See pobyso_name_free_variable_sa_so. """
602
    pobyso_name_free_variable_sa_so(freeVariableName)
603

    
604
def pobyso_name_free_variable_sa_so(freeVariableName):
605
    sollya_lib_name_free_variable(freeVariableName)
606

    
607
def pobyso_parse_string(string):
608
    """ Legacy function. See pobyso_parse_string_sa_so. """
609
    return(pobyso_parse_string_sa_so(string))
610
 
611
def pobyso_parse_string_sa_so(string):
612
    return(sollya_lib_parse_string(string))
613

    
614
def pobyso_range(rnLowerBound, rnUpperBound):
615
    """ Legacy function. See pobyso_range_sa_so. """
616
    return(pobyso_range_sa_so(rnLowerBound, rnUpperBound)) 
617

    
618
def pobyso_range_sa_so(rnLowerBoundSa, rnUpperBoundSa):
619
    """
620
    Return a Sollya range from to 2 RealField elements.
621
    The Sollya range element has a sufficient precision to hold all
622
    the digits of the bounds.
623
    """
624
    # TODO: check the bounds.
625
    lbPrec = rnLowerBoundSa.parent().precision()
626
    ubPrec = rnLowerBoundSa.parent().precision()
627
    currentSollyaPrecSa = pobyso_get_prec_so_sa()
628
    maxPrecSa = max(lbPrec, ubPrec, currentSollyaPrecSa)
629
    # Change the current Sollya precision only if necessary.
630
    if maxPrecSa > currentSollyaPrecSa:
631
        currentPrecSo = sollya_lib_get_prec(None)
632
        newPrecSo = solly_lib_constant_from_uint64(maxPrecSa)
633
        sollya_lib_set_prec(newPrecSo)
634
    lowerBoundSo = sollya_lib_constant(get_rn_value(rnLowerBoundSa))
635
    upperBoundSo = sollya_lib_constant(get_rn_value(rnUpperBoundSa))
636
    rangeSo = sollya_lib_range(lowerBoundSo, upperBoundSo)
637
    currentPrecSo = sollya_lib_get_prec(None)
638
    if maxPrecSa > currentSollyaPrecSa:
639
        sollya_lib_set_prec(currentPrecSo)
640
        sollya_lib_clear_obj(currentPrecSo)
641
        sollya_lib_clear_obj(newPrecSo)
642
    sollya_lib_clear_obj(lowerBoundSo)
643
    sollya_lib_clear_obj(upperBoundSo)
644
    return(rangeSo)
645

    
646
def pobyso_range_so_sa(rangeSo, realIntervalField = None):
647
    if realIntervalField is None:
648
        precSa = pobyso_get_prec_of_range_so_sa(rangeSo)
649
        realIntervalField = RealIntervalField(precSa)
650
    intervalSa = \
651
        pobyso_get_interval_from_range_so_sa(rangeSo, realIntervalField) 
652
    return(intervalSa)
653

    
654
def pobyso_remez_canonical_sa_sa(func, \
655
                                 degree, \
656
                                 lowerBound, \
657
                                 upperBound, \
658
                                 weight = None, \
659
                                 quality = None):
660
    """
661
    All arguments are Sage/Python.
662
    The functions (func and weight) must be passed as expressions or strings.
663
    Otherwise the function fails. 
664
    The return value is a pointer is a Sage polynomial.
665
    """
666
    var('zorglub')    # Dummy variable name for type check only.
667
    polySo = pobyso_remez_canonical_sa_so(func, \
668
                                 degree, \
669
                                 lowerBound, \
670
                                 upperBound, \
671
                                 weight = None, \
672
                                 quality = None)
673
    if parent(func) == parent("string"):
674
        functionSa = eval(func)
675
    # Expression test.
676
    elif type(func) == type(zorglub):
677
        functionSa = func
678
    maxPrecision = 0
679
    if polySo is None:
680
        return(None)
681
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
682
    RRRR = RealField(maxPrecision)
683
    polynomialRing = RRRR[functionSa.variables()[0]]
684
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, RRRR)
685
    polySa = polynomial(expSa, polynomialRing)
686
    return(polySa)
687
    
688
def pobyso_remez_canonical(func, \
689
                           degree, \
690
                           lowerBound, \
691
                           upperBound, \
692
                           weight = "1", \
693
                           quality = None):
694
    """ Legacy function. See pobyso_remez_canonical_sa_so. """
695
    return(pobyso_remez_canonical_sa_so(func, \
696
                                        degree, \
697
                                        lowerBound, \
698
                                        upperBound, \
699
                                        weight, \
700
                                        quality))
701
def pobyso_remez_canonical_sa_so(func, \
702
                                 degree, \
703
                                 lowerBound, \
704
                                 upperBound, \
705
                                 weight = None, \
706
                                 quality = None):
707
    """
708
    All arguments are Sage/Python.
709
    The functions (func and weight) must be passed as expressions or strings.
710
    Otherwise the function fails. 
711
    The return value is a pointer to a Sollya function.
712
    """
713
    var('zorglub')    # Dummy variable name for type check only.
714
    currentVariableName = None
715
    # The func argument can be of different types (string, 
716
    # symbolic expression...)
717
    if parent(func) == parent("string"):
718
        functionSo = sollya_lib_parse_string(func)
719
    # Expression test.
720
    elif type(func) == type(zorglub):
721
        # Until we are able to translate Sage expressions into Sollya 
722
        # expressions : parse the string version.
723
        currentVariableName = func.variables()[0]
724
        sollya_lib_name_free_variable(str(currentVariableName))
725
        functionSo = sollya_lib_parse_string(func._assume_str())
726
    else:
727
        return(None)
728
    if weight is None:
729
        weightSo = pobyso_constant_1_sa_so()
730
    elif parent(weight) == parent("string"):
731
        weightSo = sollya_lib_parse_string(func)
732
    elif type(weight) == type(zorglub): 
733
        functionSo = sollya_lib_parse_string_sa_so(weight._assume_str())
734
    else:
735
        return(None)
736
    degreeSo = pobyso_constant_from_int(degree)
737
    rangeSo = pobyso_range_sa_so(lowerBound, upperBound)
738
    if not quality is None:
739
        qualitySo= pobyso_constant_sa_so(quality)
740
    else:
741
        qualitySo = None
742
    return(sollya_lib_remez(functionSo, \
743
                            degreeSo, \
744
                            rangeSo, \
745
                            weightSo, \
746
                            qualitySo, \
747
                            None))
748
    
749
def pobyso_remez_canonical_so_so(funcSo, \
750
                                 degreeSo, \
751
                                 rangeSo, \
752
                                 weightSo = pobyso_constant_1_sa_so(),\
753
                                 qualitySo = None):
754
    """
755
    All arguments are pointers to Sollya objects.
756
    The return value is a pointer to a Sollya function.
757
    """
758
    if not sollya_lib_obj_is_function(funcSo):
759
        return(None)
760
    return(sollya_lib_remez(funcSo, degreeSo, rangeSo, weightSo, qualitySo, None))
761
    
762
def pobyso_set_canonical_off():
763
    sollya_lib_set_canonical(sollya_lib_off())
764

    
765
def pobyso_set_canonical_on():
766
    sollya_lib_set_canonical(sollya_lib_on())
767

    
768
def pobyso_set_prec(p):
769
    """ Legacy function. See pobyso_set_prec_sa_so. """
770
    return( pobyso_set_prec_sa_so(p))
771

    
772
def pobyso_set_prec_sa_so(p):
773
    a = c_int(p)
774
    precSo = c_void_p(sollya_lib_constant_from_int(a))
775
    sollya_lib_set_prec(precSo)
776

    
777
def pobyso_supnorm_so_so(polySo, funcSo, intervalSo, errorTypeSo, accuracySo):
778
    return(sollya_lib_supnorm(polySo, funcSo, intervalSo, errorTypeSo, \
779
                              accuracySo))
780

    
781
def pobyso_taylor_expansion_with_change_var_so_so(functionSo, degreeSo, rangeSo, \
782
                                                errorTypeSo, \
783
                                                sollyaPrecSo=None):
784
    """
785
    Compute the Taylor expansion with the variable change
786
    x -> (x-intervalCenter) included.
787
    """
788
    # No global change of the working precision.
789
    if not sollyaPrecSo is None:
790
        initialPrecSo = sollya_lib_get_prec(None)
791
        sollya_lib_set_prec(sollyaPrecSo)
792
    #
793
    intervalCenterSo = sollya_lib_mid(rangeSo)
794
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo, \
795
                                         intervalCenterSo, \
796
                                         rangeSo, errorTypeSo, None)
797
    (taylorFormListSo, numElements, isEndElliptic) = \
798
        pobyso_get_list_elements_so_so(taylorFormSo)
799
    polySo = taylorFormListSo[0]
800
    errorRangeSo = taylorFormListSo[2]
801
    maxErrorSo = sollya_lib_sup(errorRangeSo)
802
    changeVarExpSo = sollya_lib_build_function_sub(\
803
                       sollya_lib_build_function_free_variable(),\
804
                       sollya_lib_copy_obj(intervalCenterSo))
805
    polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpSo) 
806
    # If changed, reset the Sollya working precision.
807
    if not sollyaPrecSo is None:
808
        sollya_lib_set_prec(initialPrecSo)
809
        sollya_lib_clear_obj(initailPrecSo)
810
    return([polyVarChangedSo, intervalCenterSo, maxErrorSo])
811
# end pobyso_taylor_expansion_with_change_var_so_so
812

    
813
def pobyso_taylor_expansion_no_change_var_so_so(functionSo, degreeSo, rangeSo, \
814
                                                errorTypeSo, \
815
                                                sollyaPrecSo=None):
816
    """
817
    Compute the Taylor expansion without the variable change
818
    x -> x-intervalCenter.
819
    """
820
    # No global change of the working precision.
821
    if not sollyaPrecSo is None:
822
        initialPrecSo = sollya_lib_get_prec(None)
823
        sollya_lib_set_prec(sollyaPrecSo)
824
    #
825
    intervalCenterSo = sollya_lib_mid(rangeSo)
826
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo, \
827
                                         intervalCenterSo, \
828
                                         rangeSo, errorTypeSo, None)
829
    (taylorFormListSo, numElements, isEndElliptic) = \
830
        pobyso_get_list_elements_so_so(taylorFormSo)
831
    polySo = taylorFormListSo[0]
832
    errorRangeSo = taylorFormListSo[2]
833
    maxErrorSo = sollya_lib_sup(errorRangeSo)
834
    # If changed, reset the Sollya working precision.
835
    if not sollyaPrecSo is None:
836
        sollya_lib_set_prec(initialPrecSo)
837
        sollya_lib_clear_obj(initailPrecSo)
838
    return([polySo, intervalCenterSo, maxErrorSo])
839
# end pobyso_taylor_expansion_no_change_var_so_so
840

    
841
def pobyso_taylor(function, degree, point):
842
    """ Legacy function. See pobysoTaylor_so_so. """
843
    return(pobyso_taylor_so_so(function, degree, point))
844

    
845
def pobyso_taylor_so_so(functionSo, degreeSo, pointSo):
846
    return(sollya_lib_taylor(functionSo, degreeSo, pointSo))
847
    
848
def pobyso_taylorform(function, degree, point = None, interval = None, errorType=None):
849
    """ Legacy function. See ;"""
850
    
851
def pobyso_taylorform_sa_sa(functionSa, \
852
                            degree, \
853
                            point, \
854
                            precision, \
855
                            interval=None, \
856
                            errorType=None):
857
    """
858
    Compute the Taylor form of 'degree' for 'functionSa' at 'point' 
859
    for 'interval' with 'errorType'. 
860
    point: must be a Real or a Real interval.
861
    return the Taylor form as an array
862
    TODO: take care of the interval and of point when it is an interval;
863
          when errorType is not None;
864
          take care of the other elements of the Taylor form (coefficients errors and
865
          delta.
866
    """
867
    # Absolute as the default error.
868
    if errorType is None:
869
        errorTypeSo = sollya_lib_absolute()
870
    else:
871
        #TODO: deal with the other case.
872
        pass
873
    varSa = functionSa.variables()[0]
874
    pointBaseRingString = str(point.base_ring())
875
    if not re.search('Real', pointBaseRingString):
876
        return None
877
    # Call Sollya but first "sollyafy" the arguments.
878
    pobyso_name_free_variable_sa_so(str(varSa))
879
    #pobyso_set_prec_sa_so(300)
880
    # Sollyafy the function.
881
    functionSo = pobyso_parse_string_sa_so(functionSa._assume_str())
882
    if sollya_lib_obj_is_error(functionSo):
883
        print "pobyso_tailorform: function string can't be parsed!"
884
        return None
885
    # Sollyafy the degree
886
    degreeSo = sollya_lib_constant_from_int(int(degree))
887
    # Sollyafy the point
888
    if not re.search('Interval', pointBaseRingString):
889
        pointSo  = pobyso_constant_sa_so(point)
890
    else:
891
        # TODO: deal with the interval case.
892
        pass
893
    # Call Sollya
894
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, errorTypeSo,\
895
                                         None)
896
    (tfsAsList, numElements, isEndElliptic) = \
897
            pobyso_get_list_elements_so_so(taylorFormSo)
898
    polySo = tfsAsList[0]
899
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
900
    polyRealField = RealField(maxPrecision)
901
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, polyRealField)
902
    sollya_lib_close()
903
    polynomialRing = polyRealField[str(varSa)]
904
    polySa = polynomial(expSa, polynomialRing)
905
    taylorFormSa = [polySa]
906
    return(taylorFormSa)
907
# End pobyso_taylor_form_sa_sa
908

    
909
def pobyso_taylorform_so_so(functionSo, degreeSo, pointSo, intervalSo=None, \
910
                            errorTypeSo=None):
911
    createdErrorType = False
912
    if errorTypeSo is None:
913
        errorTypeSo = sollya_lib_absolute()
914
        createdErrorType = True
915
    else:
916
        #TODO: deal with the other case.
917
        pass
918
    if intervalSo is None:
919
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
920
                                         errorTypeSo, None)
921
    else:
922
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
923
                                         intervalSo, errorTypeSo, None)
924
    if createdErrorType:
925
        sollya_lib_clear_obj(errorTypeSo)
926
    return(resultSo)
927
        
928

    
929
def pobyso_univar_polynomial_print_reverse(polySa):
930
    """ Legacy function. See pobyso_univar_polynomial_print_reverse_sa_sa. """
931
    return(pobyso_univar_polynomial_print_reverse_sa_sa(polySa))
932

    
933
def pobyso_univar_polynomial_print_reverse_sa_sa(polySa):
934
    """
935
    Return the string representation of a univariate polynomial with
936
    monomials ordered in the x^0..x^n order of the monomials.
937
    Remember: Sage
938
    """
939
    polynomialRing = polySa.base_ring()
940
    # A very expensive solution:
941
    # -create a fake multivariate polynomial field with only one variable,
942
    #   specifying a negative lexicographical order;
943
    mpolynomialRing = PolynomialRing(polynomialRing.base(), \
944
                                     polynomialRing.variable_name(), \
945
                                     1, order='neglex')
946
    # - convert the univariate argument polynomial into a multivariate
947
    #   version;
948
    p = mpolynomialRing(polySa)
949
    # - return the string representation of the converted form.
950
    # There is no simple str() method defined for p's class.
951
    return(p.__str__())
952
#
953
print pobyso_get_prec()  
954
pobyso_set_prec(165)
955
print pobyso_get_prec()  
956
a=100
957
print type(a)
958
id(a)
959
print "Max arity: ", pobyso_max_arity
960
print "Function tripleDouble (43) as a string: ", pobyso_function_type_as_string(43)
961
print "Function None (44) as a string: ", pobyso_function_type_as_string(44)
962
print "...Pobyso check done"