Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 200

Historique | Voir | Annoter | Télécharger (55,4 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 widest 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
    # From Sage to Sollya bounds.
116
#    lowerBoundSo = sollya_lib_constant(get_rn_value(rnLowerBoundSa),
117
#                                       maxPrecSa)
118
    lowerBoundSo = pobyso_constant_sa_so(rnLowerBoundSa,
119
                                         maxPrecSa)
120
    upperBoundSo = pobyso_constant_sa_so(rnUpperBoundSa,
121
                                       maxPrecSa)
122
    
123
    # From Sollya bounds to range.
124
    rangeSo = sollya_lib_range(lowerBoundSo, upperBoundSo)
125
    # Back to original precision.
126
    # Clean up
127
    sollya_lib_clear_obj(lowerBoundSo)
128
    sollya_lib_clear_obj(upperBoundSo)
129
    return rangeSo
130
# End pobyso_bounds_to_range_sa_so
131

    
132
def pobyso_build_function_sub_so_so(exp1So, exp2So):
133
    return(sollya_lib_build_function_sub(exp1So, exp2So))
134

    
135
def pobyso_change_var_in_function_so_so(funcSo, chvarExpSo):
136
    """
137
    Variable change in a function.
138
    """
139
    return(sollya_lib_evaluate(funcSo,chvarExpSo))
140
# End pobyso_change_var_in_function_so_so     
141

    
142
def pobyso_chebyshevform_so_so(functionSo, degreeSo, intervalSo):
143
    resultSo = sollya_lib_chebyshevform(functionSo, degreeSo, intervalSo)
144
    return(resultSo)
145
# End pobyso_chebyshevform_so_so.
146

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

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

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

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

    
255
def pobyso_constant_1():
256
    """
257
    Obvious.
258
    Legacy function. See pobyso_constant_so_so. 
259
    """
260
    return(pobyso_constant_1_sa_so())
261

    
262
def pobyso_constant_1_sa_so():
263
    """
264
    Obvious.
265
    """
266
    return(pobyso_constant_from_int_sa_so(1))
267

    
268
def pobyso_constant_from_int(anInt):
269
    """ Legacy function. See pobyso_constant_from_int_sa_so. """
270
    return(pobyso_constant_from_int_sa_so(anInt))
271

    
272
def pobyso_constant_from_int_sa_so(anInt):
273
    """
274
    Get a Sollya constant from a Sage int.
275
    """
276
    return(sollya_lib_constant_from_int64(long(anInt)))
277

    
278
def pobyso_constant_from_int_so_sa(constSo):
279
    """
280
    Get a Sage int from a Sollya int constant.
281
    Usefull for precision or powers in polynomials.
282
    """
283
    constSa = c_long(0)
284
    sollya_lib_get_constant_as_int64(byref(constSa), constSo)
285
    return(constSa.value)
286
# End pobyso_constant_from_int_so_sa
287

    
288
def pobyso_constant_from_mpq_sa_so(rationalSa, precision = None):
289
    """
290
    Make a Sollya constant from Sage rational.
291
    A bit convoluted to take into account precision and the fact
292
    that mpq constants in Sollya a non-evaluated expressions.
293
    Function building and evaluation is needed to make it a "real"
294
    evaluated constant.
295
    """
296
    ## Deal with precision stuff.
297
    curPrecSa = None
298
    curPrecSo = None
299
    if not precision is None:
300
        curPrecSo = pobyso_get_prec_so()
301
        curPrecSaSa = c_int(0)
302
        sollya_lib_get_constant_as_int(byref(curPrecSaSa), curPrecSo)
303
        curPrecSa = int(curPrecSaSa.value)
304
        if curPrecSa != precision:
305
            pobyso_set_prec_sa_so(precision)
306
    ## In Sollya mpq constants are non-evaluated expressions.
307
    #  We must force evaluation.
308
    zeroSo    = pobyso_constant_0_sa_so()
309
    oneSo     = pobyso_constant_1_sa_so()
310
    ratExprSo = \
311
        sollya_lib_constant_from_mpq(sgmp_get_rational_value(rationalSa))
312
    addExprSo = sollya_lib_build_function_add(zeroSo, ratExprSo)
313
    constSo   = sollya_lib_evaluate(addExprSo,oneSo)
314
    ## Clears expression and arguments, as the former was created with a
315
    #  "build" variant.
316
    sollya_lib_clear_obj(addExprSo)
317
    sollya_lib_clear_obj(oneSo)
318
    if curPrecSa != precision:
319
        pobyso_set_prec_so_so(curPrecSo)
320
        sollya_lib_clear_obj(curPrecSo)
321
    return constSo
322
# End pobyso_constant_from_mpq_sa_so.
323

    
324
def pobyso_error_so():
325
    return sollya_lib_error(None)
326
# End pobyso_error().
327

    
328
def pobyso_function_type_as_string(funcType):
329
    """ Legacy function. See pobyso_function_type_as_string_so_sa. """
330
    return(pobyso_function_type_as_string_so_sa(funcType))
331

    
332
def pobyso_function_type_as_string_so_sa(funcType):
333
    """
334
    Numeric Sollya function codes -> Sage mathematical function names.
335
    Notice that pow -> ^ (a la Sage, not a la Python).
336
    """
337
    if funcType == SOLLYA_BASE_FUNC_ABS:
338
        return "abs"
339
    elif funcType == SOLLYA_BASE_FUNC_ACOS:
340
        return "arccos"
341
    elif funcType == SOLLYA_BASE_FUNC_ACOSH:
342
        return "arccosh"
343
    elif funcType == SOLLYA_BASE_FUNC_ADD:
344
        return "+"
345
    elif funcType == SOLLYA_BASE_FUNC_ASIN:
346
        return "arcsin"
347
    elif funcType == SOLLYA_BASE_FUNC_ASINH:
348
        return "arcsinh"
349
    elif funcType == SOLLYA_BASE_FUNC_ATAN:
350
        return "arctan"
351
    elif funcType == SOLLYA_BASE_FUNC_ATANH:
352
        return "arctanh"
353
    elif funcType == SOLLYA_BASE_FUNC_CEIL:
354
        return "ceil"
355
    elif funcType == SOLLYA_BASE_FUNC_CONSTANT:
356
        return "cte"
357
    elif funcType == SOLLYA_BASE_FUNC_COS:
358
        return "cos"
359
    elif funcType == SOLLYA_BASE_FUNC_COSH:
360
        return "cosh"
361
    elif funcType == SOLLYA_BASE_FUNC_DIV:
362
        return "/"
363
    elif funcType == SOLLYA_BASE_FUNC_DOUBLE:
364
        return "double"
365
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEDOUBLE:
366
        return "doubleDouble"
367
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEEXTENDED:
368
        return "doubleDxtended"
369
    elif funcType == SOLLYA_BASE_FUNC_ERF:
370
        return "erf"
371
    elif funcType == SOLLYA_BASE_FUNC_ERFC:
372
        return "erfc"
373
    elif funcType == SOLLYA_BASE_FUNC_EXP:
374
        return "exp"
375
    elif funcType == SOLLYA_BASE_FUNC_EXP_M1:
376
        return "expm1"
377
    elif funcType == SOLLYA_BASE_FUNC_FLOOR:
378
        return "floor"
379
    elif funcType == SOLLYA_BASE_FUNC_FREE_VARIABLE:
380
        return "freeVariable"
381
    elif funcType == SOLLYA_BASE_FUNC_HALFPRECISION:
382
        return "halfPrecision"
383
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYCONSTANT:
384
        return "libraryConstant"
385
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYFUNCTION:
386
        return "libraryFunction"
387
    elif funcType == SOLLYA_BASE_FUNC_LOG:
388
        return "log"
389
    elif funcType == SOLLYA_BASE_FUNC_LOG_10:
390
        return "log10"
391
    elif funcType == SOLLYA_BASE_FUNC_LOG_1P:
392
        return "log1p"
393
    elif funcType == SOLLYA_BASE_FUNC_LOG_2:
394
        return "log2"
395
    elif funcType == SOLLYA_BASE_FUNC_MUL:
396
        return "*"
397
    elif funcType == SOLLYA_BASE_FUNC_NEARESTINT:
398
        return "round"
399
    elif funcType == SOLLYA_BASE_FUNC_NEG:
400
        return "__neg__"
401
    elif funcType == SOLLYA_BASE_FUNC_PI:
402
        return "pi"
403
    elif funcType == SOLLYA_BASE_FUNC_POW:
404
        return "^"
405
    elif funcType == SOLLYA_BASE_FUNC_PROCEDUREFUNCTION:
406
        return "procedureFunction"
407
    elif funcType == SOLLYA_BASE_FUNC_QUAD:
408
        return "quad"
409
    elif funcType == SOLLYA_BASE_FUNC_SIN:
410
        return "sin"
411
    elif funcType == SOLLYA_BASE_FUNC_SINGLE:
412
        return "single"
413
    elif funcType == SOLLYA_BASE_FUNC_SINH:
414
        return "sinh"
415
    elif funcType == SOLLYA_BASE_FUNC_SQRT:
416
        return "sqrt"
417
    elif funcType == SOLLYA_BASE_FUNC_SUB:
418
        return "-"
419
    elif funcType == SOLLYA_BASE_FUNC_TAN:
420
        return "tan"
421
    elif funcType == SOLLYA_BASE_FUNC_TANH:
422
        return "tanh"
423
    elif funcType == SOLLYA_BASE_FUNC_TRIPLEDOUBLE:
424
        return "tripleDouble"
425
    else:
426
        return None
427

    
428
def pobyso_get_constant(rnArgSa, constSo):
429
    """ Legacy function. See pobyso_get_constant_so_sa. """
430
    return(pobyso_get_constant_so_sa(rnArgSa, constSo))
431

    
432
def pobyso_get_constant_so_sa(rnArgSa, constSo):
433
    """
434
    Set the value of rnArgSo to the value of constSo in MPFR_RNDN mode.
435
    rnArg must already exist and belong to some RealField.
436
    We assume that constSo points to a Sollya constant.
437
    """
438
    return(sollya_lib_get_constant(get_rn_value(rnArgSa), constSo))
439
    
440
def pobyso_get_constant_as_rn(ctExpSo):
441
    """ 
442
    Legacy function. See pobyso_get_constant_as_rn_so_sa. 
443
    """ 
444
    return(pobyso_get_constant_as_rn_so_sa(ctExpSo))
445
    
446
def pobyso_get_constant_as_rn_so_sa(constExpSo):
447
    """
448
    Get a Sollya constant as a Sage "real number".
449
    The precision of the floating-point number returned is that of the Sollya
450
    constant.
451
    """
452
    precisionSa  = pobyso_get_prec_of_constant_so_sa(constExpSo) 
453
    RRRR = RealField(precisionSa)
454
    rnSa = RRRR(0)
455
    sollya_lib_get_constant(get_rn_value(rnSa), constExpSo)
456
    return(rnSa)
457
# End pobyso_get_constant_as_rn_so_sa
458

    
459
def pobyso_get_constant_as_rn_with_rf(ctExp, realField):
460
    """ 
461
    Legacy function. See pobyso_get_constant_as_rn_with_rf_so_sa.
462
    """
463
    return(pobyso_get_constant_as_rn_with_rf_so_sa(ctExp, realField))
464
    
465
def pobyso_get_constant_as_rn_with_rf_so_sa(ctExpSo, realFieldSa = None):
466
    """
467
    Get a Sollya constant as a Sage "real number".
468
    If no real field is specified, the precision of the floating-point number 
469
    returned is that of the Sollya constant.
470
    Otherwise is is that of the real field. Hence rounding may happen.
471
    """
472
    if realFieldSa is None:
473
        sollyaPrecSa = pobyso_get_prec_of_constant_so_sa(ctExpSo)
474
        realFieldSa = RealField(sollyaPrecSa)
475
    rnSa = realFieldSa(0)
476
    sollya_lib_get_constant(get_rn_value(rnSa), ctExpSo)
477
    return(rnSa)
478
# End pobyso_get_constant_as_rn_with_rf_so_sa
479

    
480
def pobyso_get_float_poly_sa_so(polySa, realFieldSa=None):
481
    """
482
    Create a Sollya polynomial from a Sage polynomial.
483
    """
484
    pass
485
# End get_float_poly_sa_so.
486
    
487
def pobyso_get_free_variable_name():
488
    """ 
489
    Legacy function. See pobyso_get_free_variable_name_so_sa.
490
    """
491
    return(pobyso_get_free_variable_name_so_sa())
492

    
493
def pobyso_get_free_variable_name_so_sa():
494
    return(sollya_lib_get_free_variable_name())
495
    
496
def pobyso_get_function_arity(expressionSo):
497
    """ 
498
    Legacy function. See pobyso_get_function_arity_so_sa.
499
    """
500
    return(pobyso_get_function_arity_so_sa(expressionSo))
501

    
502
def pobyso_get_function_arity_so_sa(expressionSo):
503
    arity = c_int(0)
504
    sollya_lib_get_function_arity(byref(arity),expressionSo)
505
    return(int(arity.value))
506

    
507
def pobyso_get_head_function(expressionSo):
508
    """ 
509
    Legacy function. See pobyso_get_head_function_so_sa. 
510
    """
511
    return(pobyso_get_head_function_so_sa(expressionSo)) 
512

    
513
def pobyso_get_head_function_so_sa(expressionSo):
514
    functionType = c_int(0)
515
    sollya_lib_get_head_function(byref(functionType), expressionSo, None)
516
    return(int(functionType.value))
517

    
518
def pobyso_get_interval_from_range_so_sa(soRange, realIntervalFieldSa = None ):
519
    """
520
    Return the Sage interval corresponding to the Sollya range argument.
521
    If no reaIntervalField is passed as an argument, the interval bounds are not
522
    rounded: they are elements of RealIntervalField of the "right" precision
523
    to hold all the digits.
524
    """
525
    prec = c_int(0)
526
    if realIntervalFieldSa is None:
527
        retval = sollya_lib_get_prec_of_range(byref(prec), soRange, None)
528
        if retval == 0:
529
            return(None)
530
        realIntervalFieldSa = RealIntervalField(prec.value)
531
    intervalSa = realIntervalFieldSa(0,0)
532
    retval = \
533
        sollya_lib_get_interval_from_range(get_interval_value(intervalSa),\
534
                                           soRange)
535
    if retval == 0:
536
        return(None)
537
    return(intervalSa)
538
# End pobyso_get_interval_from_range_so_sa
539

    
540
def pobyso_get_list_elements(soObj):
541
    """ Legacy function. See pobyso_get_list_elements_so_so. """
542
    return(pobyso_get_list_elements_so_so(soObj))
543
 
544
def pobyso_get_list_elements_so_so(objectListSo):
545
    """
546
    Get the Sollya list elements as a Sage/Python array of Sollya objects.
547
    
548
    INPUT:
549
    - objectListSo: a Sollya list of Sollya objects.
550
    
551
    OUTPUT:
552
    - a Sage/Python tuple made of:
553
      - a Sage/Python list of Sollya objects,
554
      - a Sage/Python int holding the number of elements,
555
      - a Sage/Python int stating (!= 0) that the list is end-elliptic.
556
    NOTE::
557
        We recover the addresses of the Sollya object from the list of pointers
558
        returned by sollya_lib_get_list_elements. The list itself is freed.
559
    TODO::
560
        Figure out what to do with numElements since the number of elements
561
        can easily be recovered from the list itself. 
562
        Ditto for isEndElliptic.
563
    """
564
    listAddress = POINTER(c_longlong)()
565
    numElements = c_int(0)
566
    isEndElliptic = c_int(0)
567
    listAsSageList = []
568
    result = sollya_lib_get_list_elements(byref(listAddress),\
569
                                          byref(numElements),\
570
                                          byref(isEndElliptic),\
571
                                          objectListSo)
572
    if result == 0 :
573
        return None
574
    for i in xrange(0, numElements.value, 1):
575
       #listAsSageList.append(sollya_lib_copy_obj(listAddress[i]))
576
       listAsSageList.append(listAddress[i])
577
       # Clear each of the elements returned by Sollya.
578
       #sollya_lib_clear_obj(listAddress[i])
579
    # Free the list itself.   
580
    sollya_lib_free(listAddress)
581
    return(listAsSageList, numElements.value, isEndElliptic.value)
582

    
583
def pobyso_get_max_prec_of_exp(soExp):
584
    """ Legacy function. See pobyso_get_max_prec_of_exp_so_sa. """
585
    return(pobyso_get_max_prec_of_exp_so_sa(soExp))
586

    
587
def pobyso_get_max_prec_of_exp_so_sa(expSo):
588
    """
589
    Get the maximum precision used for the numbers in a Sollya expression.
590
    
591
    Arguments:
592
    soExp -- a Sollya expression pointer
593
    Return value:
594
    A Python integer
595
    TODO: 
596
    - error management;
597
    - correctly deal with numerical type such as DOUBLEEXTENDED.
598
    """
599
    maxPrecision = 0
600
    minConstPrec = 0
601
    currentConstPrec = 0
602
    operator = pobyso_get_head_function_so_sa(expSo)
603
    if (operator != SOLLYA_BASE_FUNC_CONSTANT) and \
604
    (operator != SOLLYA_BASE_FUNC_FREE_VARIABLE):
605
        (arity, subexpressions) = pobyso_get_subfunctions_so_sa(expSo)
606
        for i in xrange(arity):
607
            maxPrecisionCandidate = \
608
                pobyso_get_max_prec_of_exp_so_sa(subexpressions[i])
609
            if maxPrecisionCandidate > maxPrecision:
610
                maxPrecision = maxPrecisionCandidate
611
        return(maxPrecision)
612
    elif operator == SOLLYA_BASE_FUNC_CONSTANT:
613
        #minConstPrec = pobyso_get_min_prec_of_constant_so_sa(expSo)
614
        #currentConstPrec = pobyso_get_min_prec_of_constant_so_sa(soExp)
615
        #print minConstPrec, " - ", currentConstPrec 
616
        return(pobyso_get_min_prec_of_constant_so_sa(expSo))
617
    
618
    elif operator == SOLLYA_BASE_FUNC_FREE_VARIABLE:
619
        return(0)
620
    else:
621
        print "pobyso_get_max_prec_of_exp_so_sa: unexepected operator."
622
        return(0)
623

    
624
def pobyso_get_min_prec_of_constant_so_sa(constExpSo):
625
    """
626
    Get the minimum precision necessary to represent the value of a Sollya
627
    constant.
628
    MPFR_MIN_PREC and powers of 2 are taken into account.
629
    We assume that constExpSo is a point
630
    """
631
    constExpAsRnSa = pobyso_get_constant_as_rn_so_sa(constExpSo)
632
    return(min_mpfr_size(get_rn_value(constExpAsRnSa)))
633

    
634
def pobyso_get_poly_so_sa(polySo, realFieldSa=None):
635
    """
636
    Convert a Sollya polynomial into a Sage polynomial.
637
    We assume that the polynomial is in canonical form.
638
    If no realField is given, a RealField corresponding to the maximum 
639
    precision of the coefficients is internally computed.
640
    The real field is not returned but can be easily retrieved from 
641
    the polynomial itself.
642
    ALGORITHM:
643
    - (optional) compute the RealField of the coefficients;
644
    - convert the Sollya expression into a Sage expression;
645
    - convert the Sage expression into a Sage polynomial
646
    TODO: the canonical thing for the polynomial.
647
    """    
648
    if realFieldSa is None:
649
        expressionPrecSa = pobyso_get_max_prec_of_exp_so_sa(polySo)
650
        realFieldSa = RealField(expressionPrecSa)
651
    #print "Sollya expression before...",
652
    #pobyso_autoprint(polySo)
653

    
654
    expressionSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, \
655
                                                             realFieldSa)
656
    #print "...Sollya expression after.",
657
    #pobyso_autoprint(polySo)
658
    polyVariableSa = expressionSa.variables()[0]
659
    polyRingSa = realFieldSa[str(polyVariableSa)]
660
    #print polyRingSa
661
    # Do not use the polynomial(expressionSa, ring=polyRingSa) form!
662
    polynomialSa = polyRingSa(expressionSa)
663
    return(polynomialSa)
664
# End pobyso_get_poly_so_sa
665

    
666
def pobyso_get_prec():
667
    """ Legacy function. See pobyso_get_prec_so_sa(). """
668
    return(pobyso_get_prec_so_sa())
669

    
670
def pobyso_get_prec_so():
671
    """
672
    Get the current default precision in Sollya.
673
    The return value is a Sollya object.
674
    Usefull when modifying the precision back and forth by avoiding
675
    extra conversions.
676
    """
677
    return(sollya_lib_get_prec(None))
678
    
679
def pobyso_get_prec_so_sa():
680
    """
681
    Get the current default precision in Sollya.
682
    The return value is Sage/Python int.
683
    """
684
    precSo = sollya_lib_get_prec(None)
685
    precSa = c_int(0)
686
    sollya_lib_get_constant_as_int(byref(precSa), precSo)
687
    sollya_lib_clear_obj(precSo)
688
    return int(precSa.value)
689
# End pobyso_get_prec_so_sa.
690

    
691
    
692
def pobyso_get_prec_of_constant(ctExpSo):
693
    """ Legacy function. See pobyso_get_prec_of_constant_so_sa. """
694
    return(pobyso_get_prec_of_constant_so_sa(ctExpSo))
695

    
696
def pobyso_get_prec_of_constant_so_sa(ctExpSo):
697
    """
698
    Tries to find a precision to represent ctExpSo without rounding.
699
    If not possible, returns None.
700
    """
701
    prec = c_int(0)
702
    retc = sollya_lib_get_prec_of_constant(byref(prec), ctExpSo, None)
703
    if retc == 0:
704
        return(None)
705
    return(int(prec.value))
706

    
707
def pobyso_get_prec_of_range_so_sa(rangeSo):
708
    """
709
    Returns the number of bits elements of a range are coded with.
710
    """
711
    prec = c_int(0)
712
    retc = sollya_lib_get_prec_of_range(byref(prec), rangeSo, None)
713
    if retc == 0:
714
        return(None)
715
    return(int(prec.value))
716
# End pobyso_get_prec_of_range_so_sa()
717

    
718
def pobyso_get_rat_poly_sa_so(polySa, precSa = None):
719
    """
720
    Create a Sollya polynomial from a Sage rational polynomial.
721
    """
722
    ## TODO: filter arguments.
723
    ##
724
    if not precSa is None:
725
        initialPrecSo = pobyso_get_prec_so() 
726
    coefficients = polySa.coefficients()
727
    exponents    = polySa.exponents()
728
    ## TODO: deal with variables names.
729
    
730
    ## If necessary, return Sollya to its initial default precision.
731
    if not precSa is None:
732
        pobyso_set_prec_so_so(initialPrecSo)
733
# End pobyso_get_rat_poly_sa_so    
734
    
735
def pobyso_get_sage_exp_from_sollya_exp(sollyaExpSo, realField = RR):
736
    """ Legacy function. See pobyso_get_sage_exp_from_sollya_exp_so_sa. """
737
    return(pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, \
738
                                                     realField = RR))
739

    
740
def pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, realFieldSa = RR):
741
    """
742
    Get a Sage expression from a Sollya expression. 
743
    Currently only tested with polynomials with floating-point coefficients.
744
    Notice that, in the returned polynomial, the exponents are RealNumbers.
745
    """
746
    #pobyso_autoprint(sollyaExp)
747
    operatorSa = pobyso_get_head_function_so_sa(sollyaExpSo)
748
    sollyaLibFreeVariableName = sollya_lib_get_free_variable_name()
749
    # Constants and the free variable are special cases.
750
    # All other operator are dealt with in the same way.
751
    if (operatorSa != SOLLYA_BASE_FUNC_CONSTANT) and \
752
       (operatorSa != SOLLYA_BASE_FUNC_FREE_VARIABLE):
753
        (aritySa, subexpressionsSa) = pobyso_get_subfunctions_so_sa(sollyaExpSo)
754
        if aritySa == 1:
755
            sageExpSa = eval(pobyso_function_type_as_string_so_sa(operatorSa) + \
756
            "(" + pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], \
757
            realFieldSa) + ")")
758
        elif aritySa == 2:
759
            # We do not get through the preprocessor.
760
            # The "^" operator is then a special case.
761
            if operatorSa == SOLLYA_BASE_FUNC_POW:
762
                operatorAsStringSa = "**"
763
            else:
764
                operatorAsStringSa = \
765
                    pobyso_function_type_as_string_so_sa(operatorSa)
766
            sageExpSa = \
767
              eval("pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], realFieldSa)"\
768
              + " " + operatorAsStringSa + " " + \
769
                   "pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[1], realFieldSa)")
770
        # We do not know yet how to deal with arity >= 3 
771
        # (is there any in Sollya anyway?).
772
        else:
773
            sageExpSa = eval('None')
774
        return(sageExpSa)
775
    elif operatorSa == SOLLYA_BASE_FUNC_CONSTANT:
776
        #print "This is a constant"
777
        return pobyso_get_constant_as_rn_with_rf_so_sa(sollyaExpSo, realFieldSa)
778
    elif operatorSa == SOLLYA_BASE_FUNC_FREE_VARIABLE:
779
        #print "This is free variable"
780
        return(eval(sollyaLibFreeVariableName))
781
    else:
782
        print "Unexpected"
783
        return eval('None')
784
# End pobyso_get_sage_exp_from_sollya_exp_so_sa
785

    
786

    
787
def pobyso_get_subfunctions(expressionSo):
788
    """ Legacy function. See pobyso_get_subfunctions_so_sa. """
789
    return(pobyso_get_subfunctions_so_sa(expressionSo)) 
790
# End pobyso_get_subfunctions.
791
 
792
def pobyso_get_subfunctions_so_sa(expressionSo):
793
    """
794
    Get the subfunctions of an expression.
795
    Return the number of subfunctions and the list of subfunctions addresses.
796
    S.T.: Could not figure out another way than that ugly list of declarations
797
    to recover the addresses of the subfunctions.
798
    We limit ourselves to arity 8 functions. 
799
    """
800
    subf0 = c_int(0)
801
    subf1 = c_int(0)
802
    subf2 = c_int(0)
803
    subf3 = c_int(0)
804
    subf4 = c_int(0)
805
    subf5 = c_int(0)
806
    subf6 = c_int(0)
807
    subf7 = c_int(0)
808
    subf8 = c_int(0)
809
    arity = c_int(0)
810
    nullPtr = POINTER(c_int)()
811
    sollya_lib_get_subfunctions(expressionSo, byref(arity), \
812
      byref(subf0), byref(subf1), byref(subf2), byref(subf3), \
813
      byref(subf4), byref(subf5),\
814
      byref(subf6), byref(subf7), byref(subf8), nullPtr, None) 
815
#    byref(cast(subfunctions[0], POINTER(c_int))), \
816
#    byref(cast(subfunctions[0], POINTER(c_int))), \
817
#    byref(cast(subfunctions[2], POINTER(c_int))), \
818
#    byref(cast(subfunctions[3], POINTER(c_int))), \
819
#    byref(cast(subfunctions[4], POINTER(c_int))), \
820
#    byref(cast(subfunctions[5], POINTER(c_int))), \
821
#    byref(cast(subfunctions[6], POINTER(c_int))), \
822
#    byref(cast(subfunctions[7], POINTER(c_int))), \
823
#    byref(cast(subfunctions[8], POINTER(c_int))), nullPtr)
824
    subfunctions = [subf0, subf1, subf2, subf3, subf4, subf5, subf6, subf7, \
825
                    subf8]
826
    subs = []
827
    if arity.value > pobyso_max_arity:
828
        return(0,[])
829
    for i in xrange(arity.value):
830
        subs.append(int(subfunctions[i].value))
831
        #print subs[i]
832
    return(int(arity.value), subs)
833
# End pobyso_get_subfunctions_so_sa
834
    
835
def pobyso_guess_degree_sa_sa(functionSa, intervalSa, approxErrorSa, 
836
                              weightSa=None, degreeBoundSa=None):
837
    """
838
    Sa_sa variant of the solly_guessdegree function.
839
    Return 0 if something goes wrong.
840
    """
841
    functionAsStringSa = functionSa._assume_str().replace('_SAGE_VAR_', '')
842
    functionSo = pobyso_parse_string_sa_so(functionAsStringSa)
843
    if pobyso_is_error_so_sa(functionSo):
844
        sollya_lib_clear_obj(functionSo)
845
        return 0
846
    rangeSo = pobyso_interval_to_range_sa_so(intervalSa)
847
    # The approximation error is expected to be a floating point number.
848
    if pobyso_is_floating_point_number_sa_sa(approxErrorSa):
849
        approxErrorSo = pobyso_constant_sa_so(approxErrorSa)
850
    else:
851
        approxErrorSo = pobyso_constant_sa_so(RR(approxErrorSa))
852
    if not weightSa is None:
853
        weightAsStringSa = weightSa._assume_str().replace('_SAGE_VAR_', '')
854
        weightSo = pobyso_parse_string_sa_so(weightAsStringSa)
855
        if pobyso_is_error_so_sa(weightSo):
856
            sollya_lib_clear_obj(functionSo)
857
            sollya_lib_clear_obj(rangeSo)
858
            sollya_lib_clear_obj(approxErrorSo)   
859
            sollya_lib_clear_obj(weightSo)
860
            return 0   
861
    else:
862
        weightSo = None
863
    if not degreeBoundSa is None:
864
        degreeBoundSo = pobyso_constant_from_int_sa_so(degreeBoundSa)
865
    else:
866
        degreeBoundSo = None
867
    guessedDegreeSa = pobyso_guess_degree_so_sa(functionSo,
868
                                                rangeSo,
869
                                                approxErrorSo,
870
                                                weightSo,
871
                                                degreeBoundSo)
872
    sollya_lib_clear_obj(functionSo)
873
    sollya_lib_clear_obj(rangeSo)
874
    sollya_lib_clear_obj(approxErrorSo)
875
    if not weightSo is None:
876
        sollya_lib_clear_obj(weightSo)
877
    if not degreeBoundSo is None:
878
        sollya_lib_clear_obj(degreeBoundSo)
879
    return guessedDegreeSa
880
# End poyso_guess_degree_sa_sa
881

    
882
def pobyso_guess_degree_so_sa(functionSo, rangeSo, errorSo, weightSo=None, \
883
                              degreeBoundSo=None):
884
    """
885
    Thin wrapper around the guessdegree function.
886
    Nevertheless, some precision control stuff has been appended.
887
    """
888
    # Deal with Sollya internal precision issues: if it is too small, 
889
    # compared with the error, increases it to about twice -log2(error).
890
    errorSa = pobyso_get_constant_as_rn_with_rf_so_sa(errorSo)
891
    log2ErrorSa = errorSa.log2()
892
    if log2ErrorSa < 0:
893
        neededPrecisionSa = int(2 * int(-log2ErrorSa) / 64) * 64
894
    else:
895
        neededPrecisionSa = int(2 * int(log2ErrorSa) / 64) * 64
896
    #print "Needed precision:", neededPrecisionSa
897
    currentPrecSa = pobyso_get_prec_so_sa()
898
    if neededPrecisionSa > currentPrecSa:
899
        currentPrecSo = pobyso_get_prec_so()
900
        pobyso_set_prec_sa_so(neededPrecisionSa)
901
    #print "Guessing degree..."
902
    # weightSo and degreeBoundsSo are optional arguments.
903
    # As declared, sollya_lib_guessdegree must take 5 arguments.
904
    if weightSo is None:
905
        degreeRangeSo = sollya_lib_guessdegree(functionSo, rangeSo, errorSo, 
906
                                               0, 0, None)
907
    elif degreeBoundSo is None:
908
        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, \
909
                                                errorSo, weightSo, 0, None)
910
    else:
911
        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, errorSo, \
912
                                                weightSo, degreeBoundSo, None)
913
    #print "...degree guess done."
914
    # Restore internal precision, if applicable.
915
    if neededPrecisionSa > currentPrecSa:
916
        pobyso_set_prec_so_so(currentPrecSo)
917
        sollya_lib_clear_obj(currentPrecSo)
918
    degreeIntervalSa = pobyso_range_to_interval_so_sa(degreeRangeSo)
919
    sollya_lib_clear_obj(degreeRangeSo)
920
    # When ok, both bounds match.
921
    # When the degree bound is too low, the upper bound is the degree
922
    # for which the error can be honored.
923
    # When it really goes wrong, the upper bound is infinity. 
924
    if degreeIntervalSa.lower() == degreeIntervalSa.upper():
925
        return int(degreeIntervalSa.lower())
926
    else:
927
        if degreeIntervalSa.upper().is_infinity():
928
            return None
929
        else:
930
            return int(degreeIntervalSa.upper())
931
    # End pobyso_guess_degree_so_sa
932

    
933
def pobyso_infnorm_so_so(func, interval, file = None, intervalList = None):
934
    print "Do not use this function. User pobyso_supnorm_so_so instead."
935
    return(None)
936

    
937
def pobyso_interval_to_range_sa_so(intervalSa, precisionSa=None):
938
    if precisionSa is None:
939
        precisionSa = intervalSa.parent().precision()
940
    intervalSo = pobyso_bounds_to_range_sa_so(intervalSa.lower(),\
941
                                              intervalSa.upper(),\
942
                                              precisionSa)
943
    return(intervalSo)
944
# End pobyso_interval_to_range_sa_so
945

    
946
def pobyso_is_error_so_sa(objSo):
947
    """
948
    Thin wrapper around the sollya_lib_obj_is_error() function.
949
    """
950
    if sollya_lib_obj_is_error(objSo) != 0:
951
        return True
952
    else:
953
        return False
954
# End pobyso_is_error-so_sa
955

    
956
def pobyso_is_floating_point_number_sa_sa(numberSa):
957
    """
958
    Check whether a Sage number is floating point
959
    """
960
    return numberSa.parent().__class__ is RR.__class__
961

    
962
def pobyso_lib_init():
963
    sollya_lib_init(None)
964

    
965
def pobyso_lib_close():
966
    sollya_lib_close(None)
967
    
968
def pobyso_name_free_variable(freeVariableNameSa):
969
    """ Legacy function. See pobyso_name_free_variable_sa_so. """
970
    pobyso_name_free_variable_sa_so(freeVariableNameSa)
971

    
972
def pobyso_name_free_variable_sa_so(freeVariableNameSa):
973
    """
974
    Set the free variable name in Sollya from a Sage string.
975
    """
976
    sollya_lib_name_free_variable(freeVariableNameSa)
977

    
978
def pobyso_parse_string(string):
979
    """ Legacy function. See pobyso_parse_string_sa_so. """
980
    return(pobyso_parse_string_sa_so(string))
981
 
982
def pobyso_parse_string_sa_so(string):
983
    """
984
    Get the Sollya expression computed from a Sage string or
985
    a Sollya error object if parsing failed.
986
    """
987
    return(sollya_lib_parse_string(string))
988

    
989
def pobyso_precision_so_sa(ctExpSo):
990
    precisionSo = sollya_lib_precision(ctExpSo)
991
    precisionSa = pobyso_constant_from_int_so_sa(precisionSo)
992
    sollya_lib_clear_obj(precisionSo)
993
    return precisionSa
994
# End pobyso_precision_so_sa
995
    
996
def pobyso_range(rnLowerBound, rnUpperBound):
997
    """ Legacy function. See pobyso_range_sa_so. """
998
    return(pobyso_range_sa_so(rnLowerBound, rnUpperBound)) 
999

    
1000

    
1001
def pobyso_range_to_interval_so_sa(rangeSo, realIntervalFieldSa = None):
1002
    """
1003
    Get a Sage interval from a Sollya range.
1004
    If no realIntervalField is given as a parameter, the Sage interval
1005
    precision is that of the Sollya range.
1006
    Otherwise, the precision is that of the realIntervalField. In this case
1007
    rounding may happen.
1008
    """
1009
    if realIntervalFieldSa is None:
1010
        precSa = pobyso_get_prec_of_range_so_sa(rangeSo)
1011
        realIntervalFieldSa = RealIntervalField(precSa)
1012
    intervalSa = \
1013
        pobyso_get_interval_from_range_so_sa(rangeSo, realIntervalFieldSa) 
1014
    return(intervalSa)
1015

    
1016
def pobyso_remez_canonical_sa_sa(func, \
1017
                                 degree, \
1018
                                 lowerBound, \
1019
                                 upperBound, \
1020
                                 weight = None, \
1021
                                 quality = None):
1022
    """
1023
    All arguments are Sage/Python.
1024
    The functions (func and weight) must be passed as expressions or strings.
1025
    Otherwise the function fails. 
1026
    The return value is a Sage polynomial.
1027
    """
1028
    var('zorglub')    # Dummy variable name for type check only. Type of 
1029
    # zorglub is "symbolic expression".
1030
    polySo = pobyso_remez_canonical_sa_so(func, \
1031
                                 degree, \
1032
                                 lowerBound, \
1033
                                 upperBound, \
1034
                                 weight, \
1035
                                 quality)
1036
    # String test
1037
    if parent(func) == parent("string"):
1038
        functionSa = eval(func)
1039
    # Expression test.
1040
    elif type(func) == type(zorglub):
1041
        functionSa = func
1042
    else:
1043
        return None
1044
    #
1045
    maxPrecision = 0
1046
    if polySo is None:
1047
        return(None)
1048
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
1049
    RRRRSa = RealField(maxPrecision)
1050
    polynomialRingSa = RRRRSa[functionSa.variables()[0]]
1051
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, RRRRSa)
1052
    polySa = polynomial(expSa, polynomialRingSa)
1053
    sollya_lib_clear_obj(polySo)
1054
    return(polySa)
1055
# End pobyso_remez_canonical_sa_sa
1056
    
1057
def pobyso_remez_canonical(func, \
1058
                           degree, \
1059
                           lowerBound, \
1060
                           upperBound, \
1061
                           weight = "1", \
1062
                           quality = None):
1063
    """ Legacy function. See pobyso_remez_canonical_sa_so. """
1064
    return(pobyso_remez_canonical_sa_so(func, \
1065
                                        degree, \
1066
                                        lowerBound, \
1067
                                        upperBound, \
1068
                                        weight, \
1069
                                        quality))
1070
# End pobyso_remez_canonical.
1071

    
1072
def pobyso_remez_canonical_sa_so(func, \
1073
                                 degree, \
1074
                                 lowerBound, \
1075
                                 upperBound, \
1076
                                 weight = None, \
1077
                                 quality = None):
1078
    """
1079
    All arguments are Sage/Python.
1080
    The functions (func and weight) must be passed as expressions or strings.
1081
    Otherwise the function fails. 
1082
    The return value is a pointer to a Sollya function.
1083
    """
1084
    var('zorglub')    # Dummy variable name for type check only. Type of
1085
    # zorglub is "symbolic expression".
1086
    currentVariableNameSa = None
1087
    # The func argument can be of different types (string, 
1088
    # symbolic expression...)
1089
    if parent(func) == parent("string"):
1090
        localFuncSa = eval(func)
1091
        if len(localFuncSa.variables()) > 0:
1092
            currentVariableNameSa = localFuncSa.variables()[0]
1093
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1094
            functionSo = \
1095
              sollya_lib_parse_string(localFuncSa._assume_str().replace('_SAGE_VAR_', ''))
1096
    # Expression test.
1097
    elif type(func) == type(zorglub):
1098
        # Until we are able to translate Sage expressions into Sollya 
1099
        # expressions : parse the string version.
1100
        if len(func.variables()) > 0:
1101
            currentVariableNameSa = func.variables()[0]
1102
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1103
            functionSo = \
1104
              sollya_lib_parse_string(func._assume_str().replace('_SAGE_VAR_', ''))
1105
    else:
1106
        return(None)
1107
    if weight is None: # No weight given -> 1.
1108
        weightSo = pobyso_constant_1_sa_so()
1109
    elif parent(weight) == parent("string"): # Weight given as string: parse it.
1110
        weightSo = sollya_lib_parse_string(func)
1111
    elif type(weight) == type(zorglub): # Weight given as symbolice expression.
1112
        functionSo = \
1113
          sollya_lib_parse_string_sa_so(weight._assume_str().replace('_SAGE_VAR_', ''))
1114
    else:
1115
        return(None)
1116
    degreeSo = pobyso_constant_from_int(degree)
1117
    rangeSo = pobyso_bounds_to_range_sa_so(lowerBound, upperBound)
1118
    if not quality is None:
1119
        qualitySo= pobyso_constant_sa_so(quality)
1120
    else:
1121
        qualitySo = None
1122
        
1123
    remezPolySo = sollya_lib_remez(functionSo, \
1124
                                   degreeSo, \
1125
                                   rangeSo, \
1126
                                   weightSo, \
1127
                                   qualitySo, \
1128
                                   None)
1129
    sollya_lib_clear_obj(functionSo)
1130
    sollya_lib_clear_obj(degreeSo)
1131
    sollya_lib_clear_obj(rangeSo)
1132
    sollya_lib_clear_obj(weightSo)
1133
    if not qualitySo is None:
1134
        sollya_lib_clear_obj(qualitySo)
1135
    return(remezPolySo)
1136
# End pobyso_remez_canonical_sa_so
1137

    
1138
def pobyso_remez_canonical_so_so(funcSo, \
1139
                                 degreeSo, \
1140
                                 rangeSo, \
1141
                                 weightSo = pobyso_constant_1_sa_so(),\
1142
                                 qualitySo = None):
1143
    """
1144
    All arguments are pointers to Sollya objects.
1145
    The return value is a pointer to a Sollya function.
1146
    """
1147
    if not sollya_lib_obj_is_function(funcSo):
1148
        return(None)
1149
    return(sollya_lib_remez(funcSo, degreeSo, rangeSo, weightSo, qualitySo, None))
1150
# End pobyso_remez_canonical_so_so.
1151

    
1152
def pobyso_set_canonical_off():
1153
    sollya_lib_set_canonical(sollya_lib_off())
1154

    
1155
def pobyso_set_canonical_on():
1156
    sollya_lib_set_canonical(sollya_lib_on())
1157

    
1158
def pobyso_set_prec(p):
1159
    """ Legacy function. See pobyso_set_prec_sa_so. """
1160
    pobyso_set_prec_sa_so(p)
1161

    
1162
def pobyso_set_prec_sa_so(p):
1163
    a = c_int(p)
1164
    precSo = c_void_p(sollya_lib_constant_from_int(a))
1165
    sollya_lib_set_prec(precSo, None)
1166

    
1167
def pobyso_set_prec_so_so(newPrecSo):
1168
    sollya_lib_set_prec(newPrecSo, None)
1169

    
1170
def pobyso_supnorm_so_so(polySo, funcSo, intervalSo, errorTypeSo = None,\
1171
                         accuracySo = None):
1172
    """
1173
    Computes the supnorm of the approximation error between the given 
1174
    polynomial and function.
1175
    errorTypeSo defaults to "absolute".
1176
    accuracySo defaults to 2^(-40).
1177
    """
1178
    if errorTypeSo is None:
1179
        errorTypeSo = sollya_lib_absolute(None)
1180
        errorTypeIsNone = True
1181
    else:
1182
        errorTypeIsNone = False
1183
    #
1184
    if accuracySo is None:
1185
        # Notice the **!
1186
        accuracySo = pobyso_constant_sa_so(RR(2**(-40)))
1187
        accuracyIsNone = True
1188
    else:
1189
        accuracyIsNone = False
1190
    pobyso_autoprint(accuracySo)
1191
    resultSo = \
1192
        sollya_lib_supnorm(polySo, funcSo, intervalSo, errorTypeSo, \
1193
                              accuracySo)
1194
    if errorTypeIsNone:
1195
        sollya_lib_clear_obj(errorTypeSo)
1196
    if accuracyIsNone:
1197
        sollya_lib_clear_obj(accuracySo)
1198
    return resultSo
1199
# End pobyso_supnorm_so_so
1200

    
1201
def pobyso_taylor_expansion_no_change_var_so_so(functionSo, 
1202
                                                degreeSo, 
1203
                                                rangeSo,
1204
                                                errorTypeSo=None, 
1205
                                                sollyaPrecSo=None):
1206
    """
1207
    Compute the Taylor expansion without the variable change
1208
    x -> x-intervalCenter.
1209
    """
1210
    # No global change of the working precision.
1211
    if not sollyaPrecSo is None:
1212
        initialPrecSo = sollya_lib_get_prec(None)
1213
        sollya_lib_set_prec(sollyaPrecSo)
1214
    # Error type stuff: default to absolute.
1215
    if errorTypeSo is None:
1216
        errorTypeIsNone = True
1217
        errorTypeSo = sollya_lib_absolute(None)
1218
    else:
1219
        errorTypeIsNone = False
1220
    intervalCenterSo = sollya_lib_mid(rangeSo, None)
1221
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo,
1222
                                         intervalCenterSo,
1223
                                         rangeSo, errorTypeSo, None)
1224
    # taylorFormListSaSo is a Python list of Sollya objects references that 
1225
    # are copies of the elements of taylorFormSo.
1226
    # pobyso_get_list_elements_so_so clears taylorFormSo.
1227
    (taylorFormListSaSo, numElementsSa, isEndEllipticSa) = \
1228
        pobyso_get_list_elements_so_so(taylorFormSo)
1229
    polySo = sollya_lib_copy_obj(taylorFormListSaSo[0])
1230
    #print "Num elements:", numElementsSa
1231
    sollya_lib_clear_obj(taylorFormSo)
1232
    #polySo = taylorFormListSaSo[0]
1233
    #errorRangeSo = sollya_lib_copy_obj(taylorFormListSaSo[2])
1234
    errorRangeSo = taylorFormListSaSo[2]
1235
    # No copy_obj needed here: a new objects are created.
1236
    maxErrorSo    = sollya_lib_sup(errorRangeSo)
1237
    minErrorSo    = sollya_lib_inf(errorRangeSo)
1238
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
1239
    absMinErrorSo = sollya_lib_abs(minErrorSo)
1240
    sollya_lib_clear_obj(maxErrorSo)
1241
    sollya_lib_clear_obj(minErrorSo)
1242
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
1243
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
1244
    # If changed, reset the Sollya working precision.
1245
    if not sollyaPrecSo is None:
1246
        sollya_lib_set_prec(initialPrecSo)
1247
        sollya_lib_clear_obj(initialPrecSo)
1248
    if errorTypeIsNone:
1249
        sollya_lib_clear_obj(errorTypeSo)
1250
    pobyso_clear_taylorform_sa_so(taylorFormListSaSo)
1251
    if absMaxErrorSa > absMinErrorSa:
1252
        sollya_lib_clear_obj(absMinErrorSo)
1253
        return((polySo, intervalCenterSo, absMaxErrorSo))
1254
    else:
1255
        sollya_lib_clear_obj(absMaxErrorSo)
1256
        return((polySo, intervalCenterSo, absMinErrorSo))
1257
# end pobyso_taylor_expansion_no_change_var_so_so
1258

    
1259
def pobyso_taylor_expansion_with_change_var_so_so(functionSo, degreeSo, \
1260
                                                  rangeSo, \
1261
                                                  errorTypeSo=None, \
1262
                                                  sollyaPrecSo=None):
1263
    """
1264
    Compute the Taylor expansion with the variable change
1265
    x -> (x-intervalCenter) included.
1266
    """
1267
    # No global change of the working precision.
1268
    if not sollyaPrecSo is None:
1269
        initialPrecSo = sollya_lib_get_prec(None)
1270
        sollya_lib_set_prec(sollyaPrecSo)
1271
    #
1272
    # Error type stuff: default to absolute.
1273
    if errorTypeSo is None:
1274
        errorTypeIsNone = True
1275
        errorTypeSo = sollya_lib_absolute(None)
1276
    else:
1277
        errorTypeIsNone = False
1278
    intervalCenterSo = sollya_lib_mid(rangeSo)
1279
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo, \
1280
                                         intervalCenterSo, \
1281
                                         rangeSo, errorTypeSo, None)
1282
    # taylorFormListSaSo is a Python list of Sollya objects references that 
1283
    # are copies of the elements of taylorFormSo.
1284
    # pobyso_get_list_elements_so_so clears taylorFormSo.
1285
    (taylorFormListSo, numElements, isEndElliptic) = \
1286
        pobyso_get_list_elements_so_so(taylorFormSo)
1287
    polySo = taylorFormListSo[0]
1288
    errorRangeSo = taylorFormListSo[2]
1289
    maxErrorSo    = sollya_lib_sup(errorRangeSo)
1290
    minErrorSo    = sollya_lib_inf(errorRangeSo)
1291
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
1292
    absMinErrorSo = sollya_lib_abs(minErrorSo)
1293
    sollya_lib_clear_obj(maxErrorSo)
1294
    sollya_lib_clear_obj(minErrorSo)
1295
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
1296
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
1297
    changeVarExpSo = sollya_lib_build_function_sub(\
1298
                       sollya_lib_build_function_free_variable(),\
1299
                       sollya_lib_copy_obj(intervalCenterSo))
1300
    polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpSo)
1301
    sollya_lib_clear_obj(polySo) 
1302
    sollya_lib_clear_obj(changeVarExpSo)
1303
    # If changed, reset the Sollya working precision.
1304
    if not sollyaPrecSo is None:
1305
        sollya_lib_set_prec(initialPrecSo)
1306
        sollya_lib_clear_obj(initialPrecSo)
1307
    if errorTypeIsNone:
1308
        sollya_lib_clear_obj(errorTypeSo)
1309
    sollya_lib_clear_obj(taylorFormSo)
1310
    # Do not clear maxErrorSo.
1311
    if absMaxErrorSa > absMinErrorSa:
1312
        sollya_lib_clear_obj(absMinErrorSo)
1313
        return((polyVarChangedSo, intervalCenterSo, absMaxErrorSo))
1314
    else:
1315
        sollya_lib_clear_obj(absMaxErrorSo)
1316
        return((polyVarChangedSo, intervalCenterSo, absMinErrorSo))
1317
# end pobyso_taylor_expansion_with_change_var_so_so
1318

    
1319
def pobyso_taylor(function, degree, point):
1320
    """ Legacy function. See pobysoTaylor_so_so. """
1321
    return(pobyso_taylor_so_so(function, degree, point))
1322

    
1323
def pobyso_taylor_so_so(functionSo, degreeSo, pointSo):
1324
    return(sollya_lib_taylor(functionSo, degreeSo, pointSo))
1325
    
1326
def pobyso_taylorform(function, degree, point = None, 
1327
                      interval = None, errorType=None):
1328
    """ Legacy function. See pobyso_taylorform_sa_sa;"""
1329
    
1330
def pobyso_taylorform_sa_sa(functionSa, \
1331
                            degreeSa, \
1332
                            pointSa, \
1333
                            intervalSa=None, \
1334
                            errorTypeSa=None, \
1335
                            precisionSa=None):
1336
    """
1337
    Compute the Taylor form of 'degreeSa' for 'functionSa' at 'pointSa' 
1338
    for 'intervalSa' with 'errorTypeSa' (a string) using 'precisionSa'. 
1339
    point: must be a Real or a Real interval.
1340
    return the Taylor form as an array
1341
    TODO: take care of the interval and of the point when it is an interval;
1342
          when errorType is not None;
1343
          take care of the other elements of the Taylor form (coefficients 
1344
          errors and delta.
1345
    """
1346
    # Absolute as the default error.
1347
    if errorTypeSa is None:
1348
        errorTypeSo = sollya_lib_absolute()
1349
    elif errorTypeSa == "relative":
1350
        errorTypeSo = sollya_lib_relative()
1351
    elif errortypeSa == "absolute":
1352
        errorTypeSo = sollya_lib_absolute()
1353
    else:
1354
        # No clean up needed.
1355
        return None
1356
    # Global precision stuff
1357
    precisionChangedSa = False
1358
    currentSollyaPrecSo = pobyso_get_prec_so()
1359
    currentSollyaPrecSa = pobyso_constant_from_int_so_sa(currentSollyaPrecSo)
1360
    if not precisionSa is None:
1361
        if precisionSa > currentSollyaPrecSa:
1362
            pobyso_set_prec_sa_so(precisionSa)
1363
            precisionChangedSa = True
1364
            
1365
    if len(functionSa.variables()) > 0:
1366
        varSa = functionSa.variables()[0]
1367
        pobyso_name_free_variable_sa_so(str(varSa))
1368
    # In any case (point or interval) the parent of pointSa has a precision
1369
    # method.
1370
    pointPrecSa = pointSa.parent().precision()
1371
    if precisionSa > pointPrecSa:
1372
        pointPrecSa = precisionSa
1373
    # In any case (point or interval) pointSa has a base_ring() method.
1374
    pointBaseRingString = str(pointSa.base_ring())
1375
    if re.search('Interval', pointBaseRingString) is None: # Point
1376
        pointSo = pobyso_constant_sa_so(pointSa, pointPrecSa)
1377
    else: # Interval.
1378
        pointSo = pobyso_interval_to_range_sa_so(pointSa, pointPrecSa)
1379
    # Sollyafy the function.
1380
    functionSo = pobyso_parse_string_sa_so(functionSa._assume_str().replace('_SAGE_VAR_', ''))
1381
    if sollya_lib_obj_is_error(functionSo):
1382
        print "pobyso_tailorform: function string can't be parsed!"
1383
        return None
1384
    # Sollyafy the degree
1385
    degreeSo = sollya_lib_constant_from_int(int(degreeSa))
1386
    # Sollyafy the point
1387
    # Call Sollya
1388
    taylorFormSo = \
1389
        sollya_lib_taylorform(functionSo, degreeSo, pointSo, errorTypeSo,\
1390
                                         None)
1391
    sollya_lib_clear_obj(functionSo)
1392
    sollya_lib_clear_obj(degreeSo)
1393
    sollya_lib_clear_obj(pointSo)
1394
    sollya_lib_clear_obj(errorTypeSo)
1395
    (tfsAsList, numElements, isEndElliptic) = \
1396
            pobyso_get_list_elements_so_so(taylorFormSo)
1397
    polySo = tfsAsList[0]
1398
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
1399
    polyRealField = RealField(maxPrecision)
1400
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, polyRealField)
1401
    if precisionChangedSa:
1402
        sollya_lib_set_prec(currentSollyaPrecSo)
1403
        sollya_lib_clear_obj(currentSollyaPrecSo)
1404
    polynomialRing = polyRealField[str(varSa)]
1405
    polySa = polynomial(expSa, polynomialRing)
1406
    taylorFormSa = [polySa]
1407
    # Final clean-up.
1408
    sollya_lib_clear_obj(taylorFormSo)
1409
    return(taylorFormSa)
1410
# End pobyso_taylor_form_sa_sa
1411

    
1412
def pobyso_taylorform_so_so(functionSo, degreeSo, pointSo, intervalSo=None, \
1413
                            errorTypeSo=None):
1414
    createdErrorType = False
1415
    if errorTypeSo is None:
1416
        errorTypeSo = sollya_lib_absolute()
1417
        createdErrorType = True
1418
    else:
1419
        #TODO: deal with the other case.
1420
        pass
1421
    if intervalSo is None:
1422
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
1423
                                         errorTypeSo, None)
1424
    else:
1425
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
1426
                                         intervalSo, errorTypeSo, None)
1427
    if createdErrorType:
1428
        sollya_lib_clear_obj(errorTypeSo)
1429
    return(resultSo)
1430
        
1431

    
1432
def pobyso_univar_polynomial_print_reverse(polySa):
1433
    """ Legacy function. See pobyso_univar_polynomial_print_reverse_sa_sa. """
1434
    return(pobyso_univar_polynomial_print_reverse_sa_sa(polySa))
1435

    
1436
def pobyso_univar_polynomial_print_reverse_sa_sa(polySa):
1437
    """
1438
    Return the string representation of a univariate polynomial with
1439
    monomials ordered in the x^0..x^n order of the monomials.
1440
    Remember: Sage
1441
    """
1442
    polynomialRing = polySa.base_ring()
1443
    # A very expensive solution:
1444
    # -create a fake multivariate polynomial field with only one variable,
1445
    #   specifying a negative lexicographical order;
1446
    mpolynomialRing = PolynomialRing(polynomialRing.base(), \
1447
                                     polynomialRing.variable_name(), \
1448
                                     1, order='neglex')
1449
    # - convert the univariate argument polynomial into a multivariate
1450
    #   version;
1451
    p = mpolynomialRing(polySa)
1452
    # - return the string representation of the converted form.
1453
    # There is no simple str() method defined for p's class.
1454
    return(p.__str__())
1455
#
1456
print pobyso_get_prec()  
1457
pobyso_set_prec(165)
1458
print pobyso_get_prec()  
1459
a=100
1460
print type(a)
1461
id(a)
1462
print "Max arity: ", pobyso_max_arity
1463
print "Function tripleDouble (43) as a string: ", pobyso_function_type_as_string(43)
1464
print "Function None (44) as a string: ", pobyso_function_type_as_string(44)
1465
print "...Pobyso check done"