Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 124

Historique | Voir | Annoter | Télécharger (46,78 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_absolute_so_so():
89
    return(sollya_lib_absolute(None))
90

    
91
def pobyso_autoprint(arg):
92
    sollya_lib_autoprint(arg,None)
93

    
94
def pobyso_autoprint_so_so(arg):
95
    sollya_lib_autoprint(arg,None)
96
    
97
def pobyso_bounds_to_range_sa_so(rnLowerBoundSa, rnUpperBoundSa, \
98
                                 precisionSa=None):
99
    """
100
    Return a Sollya range from to 2 RealField Sage elements.
101
    The Sollya range element has a sufficient precision to hold all
102
    the digits of the Sage bounds.
103
    """
104
    # Sanity check.
105
    if rnLowerBoundSa > rnUpperBoundSa:
106
        return None
107
    # Precision stuff.
108
    if precisionSa is None:
109
        # Check for the largest precision.
110
        lbPrecSa = rnLowerBoundSa.parent().precision()
111
        ubPrecSa = rnLowerBoundSa.parent().precision()
112
        maxPrecSa = max(lbPrecSa, ubPrecSa)
113
    else:
114
        maxPrecSa = precisionSa
115
    sollyaCurrentPrecSo = pobyso_get_prec_so()
116
    sollyaCurrentPrecSa = pobyso_constant_from_int_so_sa(sollyaCurrentPrecSo)
117
    # Change the current Sollya precision only if necessary.
118
    if maxPrecSa > sollyaCurrentPrecSa:
119
        pobyso_set_prec_sa_so(maxPrecSa)
120
    # From Sage to Sollya bounds.
121
    lowerBoundSo = sollya_lib_constant(get_rn_value(rnLowerBoundSa))
122
    upperBoundSo = sollya_lib_constant(get_rn_value(rnUpperBoundSa))
123
    # From Sollya bounds to range.
124
    rangeSo = sollya_lib_range(lowerBoundSo, upperBoundSo)
125
    # Back to original precision.
126
    if maxPrecSa > sollyaCurrentPrecSa:
127
        sollya_lib_set_prec(sollyaCurrentPrecSo)
128
    # Clean up
129
    sollya_lib_clear_obj(sollyaCurrentPrecSo)
130
    sollya_lib_clear_obj(lowerBoundSo)
131
    sollya_lib_clear_obj(upperBoundSo)
132
    return(rangeSo)
133
# End pobyso_bounds_to_range_sa_so
134

    
135
def pobyso_build_function_sub_so_so(exp1So, exp2So):
136
    return(sollya_lib_build_function_sub(exp1So, exp2So))
137

    
138
def pobyso_change_var_in_function_so_so(funcSo, chvarExpSo):
139
    """
140
    Variable change in a function.
141
    """
142
    return(sollya_lib_evaluate(funcSo,chvarExpSo))
143
# End pobyso_change_var_in_function_so_so     
144

    
145
def pobyso_chebyshevform_so_so(functionSo, degreeSo, intervalSo):
146
    resultSo = sollya_lib_chebyshevform(functionSo, degreeSo, intervalSo)
147
    return(resultSo)
148
# End pobyso_chebyshevform_so_so.
149

    
150
def pobyso_clear_taylorform_sa_so(taylorFormSaSo):
151
    """
152
    This method is necessary to correctly clean up the memory from Taylor forms.
153
    These are made of a Sollya object, a Sollya object list, a Sollya object.
154
    For no clearly understood reason, sollya_lib_clear_object_list crashed 
155
    when applied to the object list.
156
    Here, we decompose it into Sage list of Sollya objects references and we
157
     clear them one by one. 
158
    """
159
    sollya_lib_clear_obj(taylorFormSaSo[0])
160
    (coefficientsErrorsListSaSo, numElementsSa, isEndEllipticSa) = \
161
        pobyso_get_list_elements_so_so(taylorFormSaSo[1])
162
    for element in coefficientsErrorsListSaSo:
163
        sollya_lib_clear_obj(element)
164
    sollya_lib_clear_obj(taylorFormSaSo[1])
165
    sollya_lib_clear_obj(taylorFormSaSo[2])
166
# End pobyso_clear_taylorform_sa_so 
167

    
168
def pobyso_cmp(rnArgSa, cteSo):
169
    """
170
    Compare the MPFR value a RealNumber with that of a Sollya constant.
171
    
172
    Get the value of the Sollya constant into a RealNumber and compare
173
    using MPFR. Could be optimized by working directly with a mpfr_t
174
    for the intermediate number. 
175
    """
176
    # Get the precision of the Sollya constant to build a Sage RealNumber 
177
    # with enough precision.to hold it.
178
    precisionOfCte = c_int(0)
179
    # From the Sollya constant, create a local Sage RealNumber.
180
    sollya_lib_get_prec_of_constant(precisionOfCte, cteSo) 
181
    #print "Precision of constant: ", precisionOfCte
182
    RRRR = RealField(precisionOfCte.value)
183
    rnLocalSa = RRRR(0)
184
    sollya_lib_get_constant(get_rn_value(rnLocalSa), cteSo)
185
    #
186
    ## Compare the Sage RealNumber version of the Sollya constant with rnArg.
187
    return(cmp_rn_value(rnArgSa, rnLocal))
188
# End pobyso_smp
189

    
190
def pobyso_compute_pos_function_abs_val_bounds_sa_sa(funcSa, lowerBoundSa, \
191
                                                     upperBoundSa):
192
    """
193
    TODO: completely rework and test.
194
    """
195
    pobyso = pobyso_name_free_variable_sa_so(funcSa.variables()[0])
196
    funcSo = pobyso_parse_string(funcSa._assume_str())
197
    rangeSo = pobyso_range_sa_so(lowerBoundSa, upperBoundSa)
198
    infnormSo = pobyso_infnorm_so_so(funcSo,rangeSo)
199
    # Sollya return the infnorm as an interval.
200
    fMaxSa = pobyso_get_interval_from_range_so_sa(infnormSo)
201
    # Get the top bound and compute the binade top limit.
202
    fMaxUpperBoundSa = fMaxSa.upper()
203
    binadeTopLimitSa = 2**ceil(fMaxUpperBoundSa.log2())
204
    # Put up together the function to use to compute the lower bound.
205
    funcAuxSo = pobyso_parse_string(str(binadeTopLimitSa) +  \
206
                                    '-(' + f._assume_str() + ')')
207
    pobyso_autoprint(funcAuxSo)
208
    # Clear the Sollya range before a new call to infnorm and issue the call.
209
    sollya_lib_clear_obj(infnormSo)
210
    infnormSo = pobyso_infnorm_so_so(funcAuxSo,rangeSo)
211
    fMinSa = pobyso_get_interval_from_range_so_sa(infnormSo)
212
    sollya_lib_clear_obj(infnormSo)
213
    fMinLowerBoundSa = binadeTopLimitSa - fMinSa.lower()
214
    # Compute the maximum of the precisions of the different bounds.
215
    maxPrecSa = max([fMinLowerBoundSa.parent().precision(), \
216
                     fMaxUpperBoundSa.parent().precision()])
217
    # Create a RealIntervalField and create an interval with the "good" bounds.
218
    RRRI = RealIntervalField(maxPrecSa)
219
    imageIntervalSa = RRRI(fMinLowerBoundSa, fMaxUpperBoundSa)
220
    # Free the unneeded Sollya objects
221
    sollya_lib_clear_obj(funcSo)
222
    sollya_lib_clear_obj(funcAuxSo)
223
    sollya_lib_clear_obj(rangeSo)
224
    return(imageIntervalSa)
225
# End pobyso_compute_pos_function_abs_val_bounds_sa_sa
226

    
227
def pobyso_constant(rnArg):
228
    """ Legacy function. See pobyso_constant_sa_so. """
229
    return(pobyso_constant_sa_so(rnArg))
230
    
231
def pobyso_constant_sa_so(rnArgSa, precisionSa=None):
232
    """
233
    Create a Sollya constant from a Sage RealNumber.
234
    """
235
    # Precision stuff
236
    if precisionSa is None:
237
        precisionSa = rnArgSa.parent().precision()
238
    currentSollyaPrecisionSo = sollya_lib_get_prec()
239
    currentSollyaPrecisionSa = \
240
        pobyso_constant_from_int(currentSollyaPrecisionSo)
241
    # Sollya constant creation takes place here.
242
    if precisionSa > currentSollyaPrecisionSa:
243
        pobyso_set_prec_sa_so(precisionSa)
244
        constantSo = sollya_lib_constant(get_rn_value(rnArgSa))
245
        pobyso_set_prec_sa_so(currentSollyaPrecision)
246
    else:
247
        constantSo = sollya_lib_constant(get_rn_value(rnArgSa))
248
    sollya_lib_clear_obj(currentSollyaPrecisionSo)
249
    return(constantSo)
250
# End pobyso_constant_sa_so
251
     
252
def pobyso_constant_0_sa_so():
253
    """
254
    Obvious.
255
    """
256
    return(pobyso_constant_from_int_sa_so(0))
257

    
258
def pobyso_constant_1():
259
    """
260
    Obvious.
261
    Legacy function. See pobyso_constant_so_so. 
262
    """
263
    return(pobyso_constant_1_sa_so())
264

    
265
def pobyso_constant_1_sa_so():
266
    """
267
    Obvious.
268
    """
269
    return(pobyso_constant_from_int_sa_so(1))
270

    
271
def pobyso_constant_from_int(anInt):
272
    """ Legacy function. See pobyso_constant_from_int_sa_so. """
273
    return(pobyso_constant_from_int_sa_so(anInt))
274

    
275
def pobyso_constant_from_int_sa_so(anInt):
276
    """
277
    Get a Sollya constant from a Sage int.
278
    """
279
    return(sollya_lib_constant_from_int(int(anInt)))
280

    
281
def pobyso_constant_from_int_so_sa(constSo):
282
    """
283
    Get a Sage int from a Sollya int constant.
284
    Usefull for precision or powers in polynomials.
285
    """
286
    constSa = c_int(0)
287
    sollya_lib_get_constant_as_int(byref(constSa), constSo)
288
    return(constSa.value)
289
# End pobyso_constant_from_int_so_sa
290

    
291
def pobyso_function_type_as_string(funcType):
292
    """ Legacy function. See pobyso_function_type_as_string_so_sa. """
293
    return(pobyso_function_type_as_string_so_sa(funcType))
294

    
295
def pobyso_function_type_as_string_so_sa(funcType):
296
    """
297
    Numeric Sollya function codes -> Sage mathematical function names.
298
    Notice that pow -> ^ (a la Sage, not a la Python).
299
    """
300
    if funcType == SOLLYA_BASE_FUNC_ABS:
301
        return "abs"
302
    elif funcType == SOLLYA_BASE_FUNC_ACOS:
303
        return "arccos"
304
    elif funcType == SOLLYA_BASE_FUNC_ACOSH:
305
        return "arccosh"
306
    elif funcType == SOLLYA_BASE_FUNC_ADD:
307
        return "+"
308
    elif funcType == SOLLYA_BASE_FUNC_ASIN:
309
        return "arcsin"
310
    elif funcType == SOLLYA_BASE_FUNC_ASINH:
311
        return "arcsinh"
312
    elif funcType == SOLLYA_BASE_FUNC_ATAN:
313
        return "arctan"
314
    elif funcType == SOLLYA_BASE_FUNC_ATANH:
315
        return "arctanh"
316
    elif funcType == SOLLYA_BASE_FUNC_CEIL:
317
        return "ceil"
318
    elif funcType == SOLLYA_BASE_FUNC_CONSTANT:
319
        return "cte"
320
    elif funcType == SOLLYA_BASE_FUNC_COS:
321
        return "cos"
322
    elif funcType == SOLLYA_BASE_FUNC_COSH:
323
        return "cosh"
324
    elif funcType == SOLLYA_BASE_FUNC_DIV:
325
        return "/"
326
    elif funcType == SOLLYA_BASE_FUNC_DOUBLE:
327
        return "double"
328
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEDOUBLE:
329
        return "doubleDouble"
330
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEEXTENDED:
331
        return "doubleDxtended"
332
    elif funcType == SOLLYA_BASE_FUNC_ERF:
333
        return "erf"
334
    elif funcType == SOLLYA_BASE_FUNC_ERFC:
335
        return "erfc"
336
    elif funcType == SOLLYA_BASE_FUNC_EXP:
337
        return "exp"
338
    elif funcType == SOLLYA_BASE_FUNC_EXP_M1:
339
        return "expm1"
340
    elif funcType == SOLLYA_BASE_FUNC_FLOOR:
341
        return "floor"
342
    elif funcType == SOLLYA_BASE_FUNC_FREE_VARIABLE:
343
        return "freeVariable"
344
    elif funcType == SOLLYA_BASE_FUNC_HALFPRECISION:
345
        return "halfPrecision"
346
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYCONSTANT:
347
        return "libraryConstant"
348
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYFUNCTION:
349
        return "libraryFunction"
350
    elif funcType == SOLLYA_BASE_FUNC_LOG:
351
        return "log"
352
    elif funcType == SOLLYA_BASE_FUNC_LOG_10:
353
        return "log10"
354
    elif funcType == SOLLYA_BASE_FUNC_LOG_1P:
355
        return "log1p"
356
    elif funcType == SOLLYA_BASE_FUNC_LOG_2:
357
        return "log2"
358
    elif funcType == SOLLYA_BASE_FUNC_MUL:
359
        return "*"
360
    elif funcType == SOLLYA_BASE_FUNC_NEARESTINT:
361
        return "round"
362
    elif funcType == SOLLYA_BASE_FUNC_NEG:
363
        return "__neg__"
364
    elif funcType == SOLLYA_BASE_FUNC_PI:
365
        return "pi"
366
    elif funcType == SOLLYA_BASE_FUNC_POW:
367
        return "^"
368
    elif funcType == SOLLYA_BASE_FUNC_PROCEDUREFUNCTION:
369
        return "procedureFunction"
370
    elif funcType == SOLLYA_BASE_FUNC_QUAD:
371
        return "quad"
372
    elif funcType == SOLLYA_BASE_FUNC_SIN:
373
        return "sin"
374
    elif funcType == SOLLYA_BASE_FUNC_SINGLE:
375
        return "single"
376
    elif funcType == SOLLYA_BASE_FUNC_SINH:
377
        return "sinh"
378
    elif funcType == SOLLYA_BASE_FUNC_SQRT:
379
        return "sqrt"
380
    elif funcType == SOLLYA_BASE_FUNC_SUB:
381
        return "-"
382
    elif funcType == SOLLYA_BASE_FUNC_TAN:
383
        return "tan"
384
    elif funcType == SOLLYA_BASE_FUNC_TANH:
385
        return "tanh"
386
    elif funcType == SOLLYA_BASE_FUNC_TRIPLEDOUBLE:
387
        return "tripleDouble"
388
    else:
389
        return None
390

    
391
def pobyso_get_constant(rnArgSa, constSo):
392
    """ Legacy function. See pobyso_get_constant_so_sa. """
393
    return(pobyso_get_constant_so_sa(rnArgSa, constSo))
394

    
395
def pobyso_get_constant_so_sa(rnArgSa, constSo):
396
    """
397
    Set the value of rnArgSo to the value of constSo in MPFR_RNDN mode.
398
    rnArg must already exist and belong to some RealField.
399
    We assume that constSo points to a Sollya constant.
400
    """
401
    return(sollya_lib_get_constant(get_rn_value(rnArgSa), constSo))
402
    
403
def pobyso_get_constant_as_rn(ctExpSo):
404
    """ 
405
    Legacy function. See pobyso_get_constant_as_rn_so_sa. 
406
    """ 
407
    return(pobyso_get_constant_as_rn_so_sa(ctExpSo))
408
    
409
def pobyso_get_constant_as_rn_so_sa(constExpSo):
410
    """
411
    Get a Sollya constant as a Sage "real number".
412
    The precision of the floating-point number returned is that of the Sollya
413
    constant.
414
    """
415
    precisionSa  = pobyso_get_prec_of_constant_so_sa(constExpSo) 
416
    RRRR = RealField(precisionSa)
417
    rnSa = RRRR(0)
418
    sollya_lib_get_constant(get_rn_value(rnSa), constExpSo)
419
    return(rnSa)
420
# End pobyso_get_constant_as_rn_so_sa
421

    
422
def pobyso_get_constant_as_rn_with_rf(ctExp, realField):
423
    """ 
424
    Legacy function. See pobyso_get_constant_as_rn_with_rf_so_sa.
425
    """
426
    return(pobyso_get_constant_as_rn_with_rf_so_sa(ctExp, realField))
427
    
428
def pobyso_get_constant_as_rn_with_rf_so_sa(ctExpSo, realFieldSa = None):
429
    """
430
    Get a Sollya constant as a Sage "real number".
431
    If no real field is specified, the precision of the floating-point number 
432
    returned is that of the Sollya constant.
433
    Otherwise is is that of the real field. Hence rounding may happen.
434
    """
435
    if realFieldSa is None:
436
        sollyaPrecSa = pobyso_get_prec_of_constant_so_sa(ctExpSo)
437
        realFieldSa = RealField(sollyaPrecSa)
438
    rnSa = realFieldSa(0)
439
    sollya_lib_get_constant(get_rn_value(rnSa), ctExpSo)
440
    return(rnSa)
441
# End pobyso_get_constant_as_rn_with_rf_so_sa
442

    
443
def pobyso_get_free_variable_name():
444
    """ 
445
    Legacy function. See pobyso_get_free_variable_name_so_sa.
446
    """
447
    return(pobyso_get_free_variable_name_so_sa())
448

    
449
def pobyso_get_free_variable_name_so_sa():
450
    return(sollya_lib_get_free_variable_name())
451
    
452
def pobyso_get_function_arity(expressionSo):
453
    """ 
454
    Legacy function. See pobyso_get_function_arity_so_sa.
455
    """
456
    return(pobyso_get_function_arity_so_sa(expressionSo))
457

    
458
def pobyso_get_function_arity_so_sa(expressionSo):
459
    arity = c_int(0)
460
    sollya_lib_get_function_arity(byref(arity),expressionSo)
461
    return(int(arity.value))
462

    
463
def pobyso_get_head_function(expressionSo):
464
    """ 
465
    Legacy function. See pobyso_get_head_function_so_sa. 
466
    """
467
    return(pobyso_get_head_function_so_sa(expressionSo)) 
468

    
469
def pobyso_get_head_function_so_sa(expressionSo):
470
    functionType = c_int(0)
471
    sollya_lib_get_head_function(byref(functionType), expressionSo, None)
472
    return(int(functionType.value))
473

    
474
def pobyso_get_interval_from_range_so_sa(soRange, realIntervalFieldSa = None ):
475
    """
476
    Return the Sage interval corresponding to the Sollya range argument.
477
    If no reaIntervalField is passed as an argument, the interval bounds are not
478
    rounded: they are elements of RealIntervalField of the "right" precision
479
    to hold all the digits.
480
    """
481
    prec = c_int(0)
482
    if realIntervalFieldSa is None:
483
        retval = sollya_lib_get_prec_of_range(byref(prec), soRange, None)
484
        if retval == 0:
485
            return(None)
486
        realIntervalFieldSa = RealIntervalField(prec.value)
487
    intervalSa = realIntervalFieldSa(0,0)
488
    retval = \
489
        sollya_lib_get_interval_from_range(get_interval_value(intervalSa),\
490
                                           soRange)
491
    if retval == 0:
492
        return(None)
493
    return(intervalSa)
494
# End pobyso_get_interval_from_range_so_sa
495

    
496
def pobyso_get_list_elements(soObj):
497
    """ Legacy function. See pobyso_get_list_elements_so_so. """
498
    return(pobyso_get_list_elements_so_so(soObj))
499
 
500
def pobyso_get_list_elements_so_so(objectListSo):
501
    """
502
    Get the Sollya list elements as a Sage/Python array of Sollya objects.
503
    
504
    INPUT:
505
    - objectListSo: a Sollya list of Sollya objects.
506
    
507
    OUTPUT:
508
    - a Sage/Python tuple made of:
509
      - a Sage/Python list of Sollya objects,
510
      - a Sage/Python int holding the number of elements,
511
      - a Sage/Python int stating (!= 0) that the list is end-elliptic.
512
    NOTE::
513
        We recover the addresses of the Sollya object from the list of pointers
514
        returned by sollya_lib_get_list_elements. The list itself is freed.
515
    TODO::
516
        Figure out what to do with numElements since the number of elements
517
        can easily be recovered from the list itself. 
518
        Ditto for isEndElliptic.
519
    """
520
    listAddress = POINTER(c_longlong)()
521
    numElements = c_int(0)
522
    isEndElliptic = c_int(0)
523
    listAsSageList = []
524
    result = sollya_lib_get_list_elements(byref(listAddress),\
525
                                          byref(numElements),\
526
                                          byref(isEndElliptic),\
527
                                          objectListSo)
528
    if result == 0 :
529
        return None
530
    for i in xrange(0, numElements.value, 1):
531
       #listAsSageList.append(sollya_lib_copy_obj(listAddress[i]))
532
       listAsSageList.append(listAddress[i])
533
       # Clear each of the elements returned by Sollya.
534
       #sollya_lib_clear_obj(listAddress[i])
535
    # Free the list itself.   
536
    sollya_lib_free(listAddress)
537
    return(listAsSageList, numElements.value, isEndElliptic.value)
538

    
539
def pobyso_get_max_prec_of_exp(soExp):
540
    """ Legacy function. See pobyso_get_max_prec_of_exp_so_sa. """
541
    return(pobyso_get_max_prec_of_exp_so_sa(soExp))
542

    
543
def pobyso_get_max_prec_of_exp_so_sa(expSo):
544
    """
545
    Get the maximum precision used for the numbers in a Sollya expression.
546
    
547
    Arguments:
548
    soExp -- a Sollya expression pointer
549
    Return value:
550
    A Python integer
551
    TODO: 
552
    - error management;
553
    - correctly deal with numerical type such as DOUBLEEXTENDED.
554
    """
555
    maxPrecision = 0
556
    minConstPrec = 0
557
    currentConstPrec = 0
558
    operator = pobyso_get_head_function_so_sa(expSo)
559
    if (operator != SOLLYA_BASE_FUNC_CONSTANT) and \
560
    (operator != SOLLYA_BASE_FUNC_FREE_VARIABLE):
561
        (arity, subexpressions) = pobyso_get_subfunctions_so_sa(expSo)
562
        for i in xrange(arity):
563
            maxPrecisionCandidate = \
564
                pobyso_get_max_prec_of_exp_so_sa(subexpressions[i])
565
            if maxPrecisionCandidate > maxPrecision:
566
                maxPrecision = maxPrecisionCandidate
567
        return(maxPrecision)
568
    elif operator == SOLLYA_BASE_FUNC_CONSTANT:
569
        #minConstPrec = pobyso_get_min_prec_of_constant_so_sa(expSo)
570
        #currentConstPrec = pobyso_get_min_prec_of_constant_so_sa(soExp)
571
        #print minConstPrec, " - ", currentConstPrec 
572
        return(pobyso_get_min_prec_of_constant_so_sa(expSo))
573
    
574
    elif operator == SOLLYA_BASE_FUNC_FREE_VARIABLE:
575
        return(0)
576
    else:
577
        print "pobyso_get_max_prec_of_exp_so_sa: unexepected operator."
578
        return(0)
579

    
580
def pobyso_get_min_prec_of_constant_so_sa(constExpSo):
581
    """
582
    Get the minimum precision necessary to represent the value of a Sollya
583
    constant.
584
    MPFR_MIN_PREC and powers of 2 are taken into account.
585
    We assume that constExpSo is a point
586
    """
587
    constExpAsRnSa = pobyso_get_constant_as_rn_so_sa(constExpSo)
588
    return(min_mpfr_size(get_rn_value(constExpAsRnSa)))
589

    
590
def pobyso_get_sage_exp_from_sollya_exp(sollyaExpSo, realField = RR):
591
    """ Legacy function. See pobyso_get_sage_exp_from_sollya_exp_so_sa. """
592
    return(pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, \
593
                                                     realField = RR))
594

    
595
def pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, realFieldSa = RR):
596
    """
597
    Get a Sage expression from a Sollya expression. 
598
    Currently only tested with polynomials with floating-point coefficients.
599
    Notice that, in the returned polynomial, the exponents are RealNumbers.
600
    """
601
    #pobyso_autoprint(sollyaExp)
602
    operatorSa = pobyso_get_head_function_so_sa(sollyaExpSo)
603
    sollyaLibFreeVariableName = sollya_lib_get_free_variable_name()
604
    # Constants and the free variable are special cases.
605
    # All other operator are dealt with in the same way.
606
    if (operatorSa != SOLLYA_BASE_FUNC_CONSTANT) and \
607
       (operatorSa != SOLLYA_BASE_FUNC_FREE_VARIABLE):
608
        (aritySa, subexpressionsSa) = pobyso_get_subfunctions_so_sa(sollyaExpSo)
609
        if aritySa == 1:
610
            sageExpSa = eval(pobyso_function_type_as_string_so_sa(operatorSa) + \
611
            "(" + pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], \
612
            realFieldSa) + ")")
613
        elif aritySa == 2:
614
            # We do not get through the preprocessor.
615
            # The "^" operator is then a special case.
616
            if operatorSa == SOLLYA_BASE_FUNC_POW:
617
                operatorAsStringSa = "**"
618
            else:
619
                operatorAsStringSa = \
620
                    pobyso_function_type_as_string_so_sa(operatorSa)
621
            sageExpSa = \
622
              eval("pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], realFieldSa)"\
623
              + " " + operatorAsStringSa + " " + \
624
                   "pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[1], realFieldSa)")
625
        # We do not know yet how to deal with arity >= 3 
626
        # (is there any in Sollya anyway?).
627
        else:
628
            sageExpSa = eval('None')
629
        return(sageExpSa)
630
    elif operatorSa == SOLLYA_BASE_FUNC_CONSTANT:
631
        #print "This is a constant"
632
        return pobyso_get_constant_as_rn_with_rf_so_sa(sollyaExpSo, realFieldSa)
633
    elif operatorSa == SOLLYA_BASE_FUNC_FREE_VARIABLE:
634
        #print "This is free variable"
635
        return(eval(sollyaLibFreeVariableName))
636
    else:
637
        print "Unexpected"
638
        return eval('None')
639
# End pobyso_get_sage_poly_from_sollya_poly
640

    
641
def pobyso_get_poly_sa_so(polySo, realFieldSa=None):
642
    """
643
    Create a Sollya polynomial from a Sage polynomial.
644
    """
645
    pass
646
# pobyso_get_poly_sa_so
647

    
648
def pobyso_get_poly_so_sa(polySo, realFieldSa=None):
649
    """
650
    Convert a Sollya polynomial into a Sage polynomial.
651
    We assume that the polynomial is in canonical form.
652
    If no realField is given, a RealField corresponding to the maximum 
653
    precision of the coefficients is internally computed.
654
    The real field is not returned but can be easily retrieved from 
655
    the polynomial itself.
656
    ALGORITHM:
657
    - (optional) compute the RealField of the coefficients;
658
    - convert the Sollya expression into a Sage expression;
659
    - convert the Sage expression into a Sage polynomial
660
    TODO: the canonical thing for the polynomial.
661
    """    
662
    if realFieldSa is None:
663
        expressionPrecSa = pobyso_get_max_prec_of_exp_so_sa(polySo)
664
        realFieldSa = RealField(expressionPrecSa)
665
    #print "Sollya expression before...",
666
    #pobyso_autoprint(polySo)
667

    
668
    expressionSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, \
669
                                                             realFieldSa)
670
    #print "...Sollya expression after.",
671
    #pobyso_autoprint(polySo)
672
    polyVariableSa = expressionSa.variables()[0]
673
    polyRingSa = realFieldSa[str(polyVariableSa)]
674
    #print polyRingSa
675
    # Do not use the polynomial(expressionSa, ring=polyRingSa) form!
676
    polynomialSa = polyRingSa(expressionSa)
677
    return(polynomialSa)
678
# End pobyso_get_sage_poly_from_sollya_poly
679

    
680
def pobyso_get_subfunctions(expressionSo):
681
    """ Legacy function. See pobyso_get_subfunctions_so_sa. """
682
    return(pobyso_get_subfunctions_so_sa(expressionSo)) 
683

    
684
def pobyso_get_subfunctions_so_sa(expressionSo):
685
    """
686
    Get the subfunctions of an expression.
687
    Return the number of subfunctions and the list of subfunctions addresses.
688
    S.T.: Could not figure out another way than that ugly list of declarations
689
    to recover the addresses of the subfunctions.
690
    We limit ourselves to arity 8 functions. 
691
    """
692
    subf0 = c_int(0)
693
    subf1 = c_int(0)
694
    subf2 = c_int(0)
695
    subf3 = c_int(0)
696
    subf4 = c_int(0)
697
    subf5 = c_int(0)
698
    subf6 = c_int(0)
699
    subf7 = c_int(0)
700
    subf8 = c_int(0)
701
    arity = c_int(0)
702
    nullPtr = POINTER(c_int)()
703
    sollya_lib_get_subfunctions(expressionSo, byref(arity), \
704
      byref(subf0), byref(subf1), byref(subf2), byref(subf3), \
705
      byref(subf4), byref(subf5),\
706
      byref(subf6), byref(subf7), byref(subf8), nullPtr, None) 
707
#    byref(cast(subfunctions[0], POINTER(c_int))), \
708
#    byref(cast(subfunctions[0], POINTER(c_int))), \
709
#    byref(cast(subfunctions[2], POINTER(c_int))), \
710
#    byref(cast(subfunctions[3], POINTER(c_int))), \
711
#    byref(cast(subfunctions[4], POINTER(c_int))), \
712
#    byref(cast(subfunctions[5], POINTER(c_int))), \
713
#    byref(cast(subfunctions[6], POINTER(c_int))), \
714
#    byref(cast(subfunctions[7], POINTER(c_int))), \
715
#    byref(cast(subfunctions[8], POINTER(c_int))), nullPtr)
716
    subfunctions = [subf0, subf1, subf2, subf3, subf4, subf5, subf6, subf7, \
717
                    subf8]
718
    subs = []
719
    if arity.value > pobyso_max_arity:
720
        return(0,[])
721
    for i in xrange(arity.value):
722
        subs.append(int(subfunctions[i].value))
723
        #print subs[i]
724
    return(int(arity.value), subs)
725
    
726
def pobyso_get_prec():
727
    """ Legacy function. See pobyso_get_prec_so_sa(). """
728
    return(pobyso_get_prec_so_sa())
729

    
730
def pobyso_get_prec_so():
731
    """
732
    Get the current default precision in Sollya.
733
    The return value is a Sollya object.
734
    Usefull when modifying the precision back and forth by avoiding
735
    extra conversions.
736
    """
737
    return(sollya_lib_get_prec(None))
738
    
739
def pobyso_get_prec_so_sa():
740
    """
741
    Get the current default precision in Sollya.
742
    The return value is Sage/Python int.
743
    """
744
    precSo = sollya_lib_get_prec(None)
745
    precSa = c_int(0)
746
    sollya_lib_get_constant_as_int(byref(precSa), precSo)
747
    sollya_lib_clear_obj(precSo)
748
    return(int(precSa.value))
749
# End pobyso_get_prec_so_sa.
750

    
751
def pobyso_get_prec_of_constant(ctExpSo):
752
    """ Legacy function. See pobyso_get_prec_of_constant_so_sa. """
753
    return(pobyso_get_prec_of_constant_so_sa(ctExpSo))
754

    
755
def pobyso_get_prec_of_constant_so_sa(ctExpSo):
756
    prec = c_int(0)
757
    retc = sollya_lib_get_prec_of_constant(byref(prec), ctExpSo, None)
758
    if retc == 0:
759
        return(None)
760
    return(int(prec.value))
761

    
762
def pobyso_get_prec_of_range_so_sa(rangeSo):
763
    prec = c_int(0)
764
    retc = sollya_lib_get_prec_of_range(byref(prec), rangeSo, None)
765
    if retc == 0:
766
        return(None)
767
    return(int(prec.value))
768

    
769
def pobyso_infnorm_so_so(func, interval, file = None, intervalList = None):
770
    print "Do not use this function. User pobyso_supnorm_so_so instead."
771
    return(None)
772

    
773
def pobyso_interval_to_range_sa_so(intervalSa, precisionSa=None):
774
    if precisionSa is None:
775
        precisionSa = intervalSa.parent().precision()
776
    intervalSo = pobyso_bounds_to_range_sa_so(intervalSa.lower(),\
777
                                              intervalSa.upper(),\
778
                                              precisionSa)
779
    return(intervalSo)
780
# End pobyso_interval_to_range_sa_so
781

    
782
def pobyso_lib_init():
783
    sollya_lib_init(None)
784

    
785
def pobyso_lib_close():
786
    sollya_lib_close(None)
787
    
788
def pobyso_name_free_variable(freeVariableNameSa):
789
    """ Legacy function. See pobyso_name_free_variable_sa_so. """
790
    pobyso_name_free_variable_sa_so(freeVariableNameSa)
791

    
792
def pobyso_name_free_variable_sa_so(freeVariableNameSa):
793
    """
794
    Set the free variable name in Sollya from a Sage string.
795
    """
796
    sollya_lib_name_free_variable(freeVariableNameSa)
797

    
798
def pobyso_parse_string(string):
799
    """ Legacy function. See pobyso_parse_string_sa_so. """
800
    return(pobyso_parse_string_sa_so(string))
801
 
802
def pobyso_parse_string_sa_so(string):
803
    """
804
    Get the Sollya expression computed from a Sage string.
805
    """
806
    return(sollya_lib_parse_string(string))
807

    
808
def pobyso_range(rnLowerBound, rnUpperBound):
809
    """ Legacy function. See pobyso_range_sa_so. """
810
    return(pobyso_range_sa_so(rnLowerBound, rnUpperBound)) 
811

    
812

    
813
def pobyso_range_to_interval_so_sa(rangeSo, realIntervalFieldSa = None):
814
    """
815
    Get a Sage interval from a Sollya range.
816
    If no realIntervalField is given as a parameter, the Sage interval
817
    precision is that of the Sollya range.
818
    Otherwise, the precision is that of the realIntervalField. In this case
819
    rounding may happen.
820
    """
821
    if realIntervalFieldSa is None:
822
        precSa = pobyso_get_prec_of_range_so_sa(rangeSo)
823
        realIntervalFieldSa = RealIntervalField(precSa)
824
    intervalSa = \
825
        pobyso_get_interval_from_range_so_sa(rangeSo, realIntervalFieldSa) 
826
    return(intervalSa)
827

    
828
def pobyso_remez_canonical_sa_sa(func, \
829
                                 degree, \
830
                                 lowerBound, \
831
                                 upperBound, \
832
                                 weight = None, \
833
                                 quality = None):
834
    """
835
    All arguments are Sage/Python.
836
    The functions (func and weight) must be passed as expressions or strings.
837
    Otherwise the function fails. 
838
    The return value is a Sage polynomial.
839
    """
840
    var('zorglub')    # Dummy variable name for type check only. Type of 
841
    # zorglub is "symbolic expression".
842
    polySo = pobyso_remez_canonical_sa_so(func, \
843
                                 degree, \
844
                                 lowerBound, \
845
                                 upperBound, \
846
                                 weight, \
847
                                 quality)
848
    # String test
849
    if parent(func) == parent("string"):
850
        functionSa = eval(func)
851
    # Expression test.
852
    elif type(func) == type(zorglub):
853
        functionSa = func
854
    else:
855
        return None
856
    #
857
    maxPrecision = 0
858
    if polySo is None:
859
        return(None)
860
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
861
    RRRRSa = RealField(maxPrecision)
862
    polynomialRingSa = RRRRSa[functionSa.variables()[0]]
863
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, RRRRSa)
864
    polySa = polynomial(expSa, polynomialRingSa)
865
    sollya_lib_clear_obj(polySo)
866
    return(polySa)
867
# End pobyso_remez_canonical_sa_sa
868
    
869
def pobyso_remez_canonical(func, \
870
                           degree, \
871
                           lowerBound, \
872
                           upperBound, \
873
                           weight = "1", \
874
                           quality = None):
875
    """ Legacy function. See pobyso_remez_canonical_sa_so. """
876
    return(pobyso_remez_canonical_sa_so(func, \
877
                                        degree, \
878
                                        lowerBound, \
879
                                        upperBound, \
880
                                        weight, \
881
                                        quality))
882
def pobyso_remez_canonical_sa_so(func, \
883
                                 degree, \
884
                                 lowerBound, \
885
                                 upperBound, \
886
                                 weight = None, \
887
                                 quality = None):
888
    """
889
    All arguments are Sage/Python.
890
    The functions (func and weight) must be passed as expressions or strings.
891
    Otherwise the function fails. 
892
    The return value is a pointer to a Sollya function.
893
    """
894
    var('zorglub')    # Dummy variable name for type check only. Type of
895
    # zorglub is "symbolic expression".
896
    currentVariableNameSa = None
897
    # The func argument can be of different types (string, 
898
    # symbolic expression...)
899
    if parent(func) == parent("string"):
900
        localFuncSa = eval(func)
901
        if len(localFuncSa.variables()) > 0:
902
            currentVariableNameSa = localFuncSa.variables()[0]
903
            sollya_lib_name_free_variable(str(currentVariableNameSa))
904
            functionSo = sollya_lib_parse_string(localFuncSa._assume_str())
905
    # Expression test.
906
    elif type(func) == type(zorglub):
907
        # Until we are able to translate Sage expressions into Sollya 
908
        # expressions : parse the string version.
909
        if len(func.variables()) > 0:
910
            currentVariableNameSa = func.variables()[0]
911
            sollya_lib_name_free_variable(str(currentVariableNameSa))
912
            functionSo = sollya_lib_parse_string(func._assume_str())
913
    else:
914
        return(None)
915
    if weight is None: # No weight given -> 1.
916
        weightSo = pobyso_constant_1_sa_so()
917
    elif parent(weight) == parent("string"): # Weight given as string: parse it.
918
        weightSo = sollya_lib_parse_string(func)
919
    elif type(weight) == type(zorglub): # Weight given as symbolice expression.
920
        functionSo = sollya_lib_parse_string_sa_so(weight._assume_str())
921
    else:
922
        return(None)
923
    degreeSo = pobyso_constant_from_int(degree)
924
    rangeSo = pobyso_bounds_to_range_sa_so(lowerBound, upperBound)
925
    if not quality is None:
926
        qualitySo= pobyso_constant_sa_so(quality)
927
    else:
928
        qualitySo = None
929
        
930
    remezPolySo = sollya_lib_remez(functionSo, \
931
                                   degreeSo, \
932
                                   rangeSo, \
933
                                   weightSo, \
934
                                   qualitySo, \
935
                                   None)
936
    sollya_lib_clear_obj(functionSo)
937
    sollya_lib_clear_obj(degreeSo)
938
    sollya_lib_clear_obj(rangeSo)
939
    sollya_lib_clear_obj(weightSo)
940
    if not qualitySo is None:
941
        sollya_lib_clear_obj(qualitySo)
942
    return(remezPolySo)
943
# End pobyso_remez_canonical_sa_so
944

    
945
def pobyso_remez_canonical_so_so(funcSo, \
946
                                 degreeSo, \
947
                                 rangeSo, \
948
                                 weightSo = pobyso_constant_1_sa_so(),\
949
                                 qualitySo = None):
950
    """
951
    All arguments are pointers to Sollya objects.
952
    The return value is a pointer to a Sollya function.
953
    """
954
    if not sollya_lib_obj_is_function(funcSo):
955
        return(None)
956
    return(sollya_lib_remez(funcSo, degreeSo, rangeSo, weightSo, qualitySo, None))
957
    
958
def pobyso_set_canonical_off():
959
    sollya_lib_set_canonical(sollya_lib_off())
960

    
961
def pobyso_set_canonical_on():
962
    sollya_lib_set_canonical(sollya_lib_on())
963

    
964
def pobyso_set_prec(p):
965
    """ Legacy function. See pobyso_set_prec_sa_so. """
966
    pobyso_set_prec_sa_so(p)
967

    
968
def pobyso_set_prec_sa_so(p):
969
    a = c_int(p)
970
    precSo = c_void_p(sollya_lib_constant_from_int(a))
971
    sollya_lib_set_prec(precSo, None)
972

    
973
def pobyso_set_prec_so_so(newPrecSo):
974
    sollya_lib_set_prec(newPrecSo, None)
975

    
976
def pobyso_supnorm_so_so(polySo, funcSo, intervalSo, errorTypeSo = None,\
977
                         accuracySo = None):
978
    """
979
    Computes the supnorm of the approximation error between the given 
980
    polynomial and function.
981
    errorTypeSo defaults to "absolute".
982
    accuracySo defaults to 2^(-40).
983
    """
984
    if errorTypeSo is None:
985
        errorTypeSo = sollya_lib_absolute(None)
986
        errorTypeIsNone = True
987
    else:
988
        errorTypeIsNone = False
989
    #
990
    if accuracySo is None:
991
        # Notice the **!
992
        accuracySo = pobyso_constant_sa_so(RR(2**(-40)))
993
        accuracyIsNone = True
994
    else:
995
        accuracyIsNone = False
996
    pobyso_autoprint(accuracySo)
997
    resultSo = \
998
        sollya_lib_supnorm(polySo, funcSo, intervalSo, errorTypeSo, \
999
                              accuracySo)
1000
    if errorTypeIsNone:
1001
        sollya_lib_clear_obj(errorTypeSo)
1002
    if accuracyIsNone:
1003
        sollya_lib_clear_obj(accuracySo)
1004
    return resultSo
1005
# End pobyso_supnorm_so_so
1006

    
1007
def pobyso_taylor_expansion_with_change_var_so_so(functionSo, degreeSo, \
1008
                                                  rangeSo, \
1009
                                                  errorTypeSo=None, \
1010
                                                  sollyaPrecSo=None):
1011
    """
1012
    Compute the Taylor expansion with the variable change
1013
    x -> (x-intervalCenter) included.
1014
    """
1015
    # No global change of the working precision.
1016
    if not sollyaPrecSo is None:
1017
        initialPrecSo = sollya_lib_get_prec(None)
1018
        sollya_lib_set_prec(sollyaPrecSo)
1019
    #
1020
    # Error type stuff: default to absolute.
1021
    if errorTypeSo is None:
1022
        errorTypeIsNone = True
1023
        errorTypeSo = sollya_lib_absolute(None)
1024
    else:
1025
        errorTypeIsNone = False
1026
    intervalCenterSo = sollya_lib_mid(rangeSo)
1027
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo, \
1028
                                         intervalCenterSo, \
1029
                                         rangeSo, errorTypeSo, None)
1030
    # taylorFormListSaSo is a Python list of Sollya objects references that 
1031
    # are copies of the elements of taylorFormSo.
1032
    # pobyso_get_list_elements_so_so clears taylorFormSo.
1033
    (taylorFormListSo, numElements, isEndElliptic) = \
1034
        pobyso_get_list_elements_so_so(taylorFormSo)
1035
    polySo = taylorFormListSo[0]
1036
    errorRangeSo = taylorFormListSo[2]
1037
    maxErrorSo = sollya_lib_sup(errorRangeSo)
1038
    changeVarExpSo = sollya_lib_build_function_sub(\
1039
                       sollya_lib_build_function_free_variable(),\
1040
                       sollya_lib_copy_obj(intervalCenterSo))
1041
    polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpSo) 
1042
    sollya_lib_clear_obj(changeVarExpSo)
1043
    # If changed, reset the Sollya working precision.
1044
    if not sollyaPrecSo is None:
1045
        sollya_lib_set_prec(initialPrecSo)
1046
        sollya_lib_clear_obj(initialPrecSo)
1047
    if errorTypeIsNone:
1048
        sollya_lib_clear_obj(errorTypeSo)
1049
    sollya_lib_clear_obj(taylorFormSo)
1050
    # Do not clear maxErrorSo.
1051
    return((polyVarChangedSo, intervalCenterSo, maxErrorSo))
1052
# end pobyso_taylor_expansion_with_change_var_so_so
1053

    
1054
def pobyso_taylor_expansion_no_change_var_so_so(functionSo, degreeSo, rangeSo,
1055
                                                errorTypeSo=None, 
1056
                                                sollyaPrecSo=None):
1057
    """
1058
    Compute the Taylor expansion without the variable change
1059
    x -> x-intervalCenter.
1060
    """
1061
    # No global change of the working precision.
1062
    if not sollyaPrecSo is None:
1063
        initialPrecSo = sollya_lib_get_prec(None)
1064
        sollya_lib_set_prec(sollyaPrecSo)
1065
    # Error type stuff: default to absolute.
1066
    if errorTypeSo is None:
1067
        errorTypeIsNone = True
1068
        errorTypeSo = sollya_lib_absolute(None)
1069
    else:
1070
        errorTypeIsNone = False
1071
    intervalCenterSo = sollya_lib_mid(rangeSo, None)
1072
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo,
1073
                                         intervalCenterSo,
1074
                                         rangeSo, errorTypeSo, None)
1075
    # taylorFormListSaSo is a Python list of Sollya objects references that 
1076
    # are copies of the elements of taylorFormSo.
1077
    # pobyso_get_list_elements_so_so clears taylorFormSo.
1078
    (taylorFormListSaSo, numElementsSa, isEndEllipticSa) = \
1079
        pobyso_get_list_elements_so_so(taylorFormSo)
1080
    polySo = sollya_lib_copy_obj(taylorFormListSaSo[0])
1081
    #print "Num elements:", numElementsSa
1082
    sollya_lib_clear_obj(taylorFormSo)
1083
    #polySo = taylorFormListSaSo[0]
1084
    #errorRangeSo = sollya_lib_copy_obj(taylorFormListSaSo[2])
1085
    errorRangeSo = taylorFormListSaSo[2]
1086
    # No copy_obj needed here: a new object is created.
1087
    maxErrorSo = sollya_lib_sup(errorRangeSo)
1088
    # If changed, reset the Sollya working precision.
1089
    if not sollyaPrecSo is None:
1090
        sollya_lib_set_prec(initialPrecSo)
1091
        sollya_lib_clear_obj(initialPrecSo)
1092
    if errorTypeIsNone:
1093
        sollya_lib_clear_obj(errorTypeSo)
1094
    pobyso_clear_taylorform_sa_so(taylorFormListSaSo)
1095
    return((polySo, intervalCenterSo, maxErrorSo))
1096
# end pobyso_taylor_expansion_no_change_var_so_so
1097

    
1098
def pobyso_taylor(function, degree, point):
1099
    """ Legacy function. See pobysoTaylor_so_so. """
1100
    return(pobyso_taylor_so_so(function, degree, point))
1101

    
1102
def pobyso_taylor_so_so(functionSo, degreeSo, pointSo):
1103
    return(sollya_lib_taylor(functionSo, degreeSo, pointSo))
1104
    
1105
def pobyso_taylorform(function, degree, point = None, 
1106
                      interval = None, errorType=None):
1107
    """ Legacy function. See pobyso_taylorform_sa_sa;"""
1108
    
1109
def pobyso_taylorform_sa_sa(functionSa, \
1110
                            degreeSa, \
1111
                            pointSa, \
1112
                            intervalSa=None, \
1113
                            errorTypeSa=None, \
1114
                            precisionSa=None):
1115
    """
1116
    Compute the Taylor form of 'degreeSa' for 'functionSa' at 'pointSa' 
1117
    for 'intervalSa' with 'errorTypeSa' (a string) using 'precisionSa'. 
1118
    point: must be a Real or a Real interval.
1119
    return the Taylor form as an array
1120
    TODO: take care of the interval and of the point when it is an interval;
1121
          when errorType is not None;
1122
          take care of the other elements of the Taylor form (coefficients 
1123
          errors and delta.
1124
    """
1125
    # Absolute as the default error.
1126
    if errorTypeSa is None:
1127
        errorTypeSo = sollya_lib_absolute()
1128
    elif errorTypeSa == "relative":
1129
        errorTypeSo = sollya_lib_relative()
1130
    elif errortypeSa == "absolute":
1131
        errorTypeSo = sollya_lib_absolute()
1132
    else:
1133
        # No clean up needed.
1134
        return None
1135
    # Global precision stuff
1136
    precisionChangedSa = False
1137
    currentSollyaPrecSo = pobyso_get_prec_so()
1138
    currentSollyaPrecSa = pobyso_constant_from_int_so_sa(currentSollyaPrecSo)
1139
    if not precisionSa is None:
1140
        if precisionSa > currentSollyaPrecSa:
1141
            pobyso_set_prec_sa_so(precisionSa)
1142
            precisionChangedSa = True
1143
            
1144
    if len(functionSa.variables()) > 0:
1145
        varSa = functionSa.variables()[0]
1146
        pobyso_name_free_variable_sa_so(str(varSa))
1147
    # In any case (point or interval) the parent of pointSa has a precision
1148
    # method.
1149
    pointPrecSa = pointSa.parent().precision()
1150
    if precisionSa > pointPrecSa:
1151
        pointPrecSa = precisionSa
1152
    # In any case (point or interval) pointSa has a base_ring() method.
1153
    pointBaseRingString = str(pointSa.base_ring())
1154
    if re.search('Interval', pointBaseRingString) is None: # Point
1155
        pointSo = pobyso_constant_sa_so(pointSa, pointPrecSa)
1156
    else: # Interval.
1157
        pointSo = pobyso_interval_to_range_sa_so(pointSa, pointPrecSa)
1158
    # Sollyafy the function.
1159
    functionSo = pobyso_parse_string_sa_so(functionSa._assume_str())
1160
    if sollya_lib_obj_is_error(functionSo):
1161
        print "pobyso_tailorform: function string can't be parsed!"
1162
        return None
1163
    # Sollyafy the degree
1164
    degreeSo = sollya_lib_constant_from_int(int(degreeSa))
1165
    # Sollyafy the point
1166
    # Call Sollya
1167
    taylorFormSo = \
1168
        sollya_lib_taylorform(functionSo, degreeSo, pointSo, errorTypeSo,\
1169
                                         None)
1170
    sollya_lib_clear_obj(functionSo)
1171
    sollya_lib_clear_obj(degreeSo)
1172
    sollya_lib_clear_obj(pointSo)
1173
    sollya_lib_clear_obj(errorTypeSo)
1174
    (tfsAsList, numElements, isEndElliptic) = \
1175
            pobyso_get_list_elements_so_so(taylorFormSo)
1176
    polySo = tfsAsList[0]
1177
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
1178
    polyRealField = RealField(maxPrecision)
1179
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, polyRealField)
1180
    if precisionChangedSa:
1181
        sollya_lib_set_prec(currentSollyaPrecSo)
1182
        sollya_lib_clear_obj(currentSollyaPrecSo)
1183
    polynomialRing = polyRealField[str(varSa)]
1184
    polySa = polynomial(expSa, polynomialRing)
1185
    taylorFormSa = [polySa]
1186
    # Final clean-up.
1187
    sollya_lib_clear_obj(taylorFormSo)
1188
    return(taylorFormSa)
1189
# End pobyso_taylor_form_sa_sa
1190

    
1191
def pobyso_taylorform_so_so(functionSo, degreeSo, pointSo, intervalSo=None, \
1192
                            errorTypeSo=None):
1193
    createdErrorType = False
1194
    if errorTypeSo is None:
1195
        errorTypeSo = sollya_lib_absolute()
1196
        createdErrorType = True
1197
    else:
1198
        #TODO: deal with the other case.
1199
        pass
1200
    if intervalSo is None:
1201
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
1202
                                         errorTypeSo, None)
1203
    else:
1204
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
1205
                                         intervalSo, errorTypeSo, None)
1206
    if createdErrorType:
1207
        sollya_lib_clear_obj(errorTypeSo)
1208
    return(resultSo)
1209
        
1210

    
1211
def pobyso_univar_polynomial_print_reverse(polySa):
1212
    """ Legacy function. See pobyso_univar_polynomial_print_reverse_sa_sa. """
1213
    return(pobyso_univar_polynomial_print_reverse_sa_sa(polySa))
1214

    
1215
def pobyso_univar_polynomial_print_reverse_sa_sa(polySa):
1216
    """
1217
    Return the string representation of a univariate polynomial with
1218
    monomials ordered in the x^0..x^n order of the monomials.
1219
    Remember: Sage
1220
    """
1221
    polynomialRing = polySa.base_ring()
1222
    # A very expensive solution:
1223
    # -create a fake multivariate polynomial field with only one variable,
1224
    #   specifying a negative lexicographical order;
1225
    mpolynomialRing = PolynomialRing(polynomialRing.base(), \
1226
                                     polynomialRing.variable_name(), \
1227
                                     1, order='neglex')
1228
    # - convert the univariate argument polynomial into a multivariate
1229
    #   version;
1230
    p = mpolynomialRing(polySa)
1231
    # - return the string representation of the converted form.
1232
    # There is no simple str() method defined for p's class.
1233
    return(p.__str__())
1234
#
1235
print pobyso_get_prec()  
1236
pobyso_set_prec(165)
1237
print pobyso_get_prec()  
1238
a=100
1239
print type(a)
1240
id(a)
1241
print "Max arity: ", pobyso_max_arity
1242
print "Function tripleDouble (43) as a string: ", pobyso_function_type_as_string(43)
1243
print "Function None (44) as a string: ", pobyso_function_type_as_string(44)
1244
print "...Pobyso check done"