Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 222

Historique | Voir | Annoter | Télécharger (69,4 ko)

1
"""
2
@file pobyso.py
3
Actual functions to use in Sage
4
ST 2012-11-13
5

6
Command line syntax:
7
  use from Sage (via the "load" or the "attach" commands)
8

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

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

    
87
pobyso_max_arity = 9
88

    
89
def pobyso_absolute_so_so():
90
    return(sollya_lib_absolute(None))
91

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

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

    
133
def pobyso_build_end_elliptic_list_so_so(*args):
134
    """
135
    From argumrny Sollya objects, create a Sollya end elliptic list.
136
    Elements of the list are "eaten" (should not be cleared individualy, 
137
    are cleared when the list is cleared).
138
    """
139
    if len(args) == 0:
140
        ## Called with an empty list produced "error".
141
        return sollya_lib_build_end_elliptic_list(None)
142
    index = 0
143
    ## One can not append elements to an elliptic list, prepend only is 
144
    #  permitted.
145
    for argument in reversed(args):
146
        if index == 0:
147
            listSo = sollya_lib_build_end_elliptic_list(argument, None)
148
        else:
149
            listSo = sollya_lib_prepend(argument, listSo)
150
        index += 1
151
    return listSo
152
        
153
# End pobyso_build_end_elliptic_list_so_so
154

    
155
def pobyso_build_function_sub_so_so(exp1So, exp2So):
156
    return(sollya_lib_build_function_sub(exp1So, exp2So))
157

    
158
def pobyso_change_var_in_function_so_so(funcSo, chvarExpSo):
159
    """
160
    Variable change in a function.
161
    """
162
    return(sollya_lib_evaluate(funcSo,chvarExpSo))
163
# End pobyso_change_var_in_function_so_so     
164

    
165
def pobyso_chebyshevform_so_so(functionSo, degreeSo, intervalSo):
166
    resultSo = sollya_lib_chebyshevform(functionSo, degreeSo, intervalSo)
167
    return(resultSo)
168
# End pobyso_chebyshevform_so_so.
169

    
170
def pobyso_clear_taylorform_sa_so(taylorFormSaSo):
171
    """
172
    This method is necessary to correctly clean up the memory from Taylor forms.
173
    These are made of a Sollya object, a Sollya object list, a Sollya object.
174
    For no clearly understood reason, sollya_lib_clear_object_list crashed 
175
    when applied to the object list.
176
    Here, we decompose it into Sage list of Sollya objects references and we
177
     clear them one by one. 
178
    """
179
    sollya_lib_clear_obj(taylorFormSaSo[0])
180
    (coefficientsErrorsListSaSo, numElementsSa, isEndEllipticSa) = \
181
        pobyso_get_list_elements_so_so(taylorFormSaSo[1])
182
    for element in coefficientsErrorsListSaSo:
183
        sollya_lib_clear_obj(element)
184
    sollya_lib_clear_obj(taylorFormSaSo[1])
185
    sollya_lib_clear_obj(taylorFormSaSo[2])
186
# End pobyso_clear_taylorform_sa_so 
187

    
188
def pobyso_cmp(rnArgSa, cteSo):
189
    """
190
    Compare the MPFR value a RealNumber with that of a Sollya constant.
191
    
192
    Get the value of the Sollya constant into a RealNumber and compare
193
    using MPFR. Could be optimized by working directly with a mpfr_t
194
    for the intermediate number. 
195
    """
196
    # Get the precision of the Sollya constant to build a Sage RealNumber 
197
    # with enough precision.to hold it.
198
    precisionOfCte = c_int(0)
199
    # From the Sollya constant, create a local Sage RealNumber.
200
    sollya_lib_get_prec_of_constant(precisionOfCte, cteSo) 
201
    #print "Precision of constant: ", precisionOfCte
202
    RRRR = RealField(precisionOfCte.value)
203
    rnLocalSa = RRRR(0)
204
    sollya_lib_get_constant(get_rn_value(rnLocalSa), cteSo)
205
    #
206
    ## Compare the Sage RealNumber version of the Sollya constant with rnArg.
207
    return(cmp_rn_value(rnArgSa, rnLocal))
208
# End pobyso_smp
209

    
210
def pobyso_compute_pos_function_abs_val_bounds_sa_sa(funcSa, lowerBoundSa, \
211
                                                     upperBoundSa):
212
    """
213
    TODO: completely rework and test.
214
    """
215
    pobyso = pobyso_name_free_variable_sa_so(funcSa.variables()[0])
216
    funcSo = pobyso_parse_string(funcSa._assume_str().replace('_SAGE_VAR_', ''))
217
    rangeSo = pobyso_range_sa_so(lowerBoundSa, upperBoundSa)
218
    infnormSo = pobyso_infnorm_so_so(funcSo,rangeSo)
219
    # Sollya return the infnorm as an interval.
220
    fMaxSa = pobyso_get_interval_from_range_so_sa(infnormSo)
221
    # Get the top bound and compute the binade top limit.
222
    fMaxUpperBoundSa = fMaxSa.upper()
223
    binadeTopLimitSa = 2**ceil(fMaxUpperBoundSa.log2())
224
    # Put up together the function to use to compute the lower bound.
225
    funcAuxSo = pobyso_parse_string(str(binadeTopLimitSa) +  \
226
                                    '-(' + f._assume_str().replace('_SAGE_VAR_', '') + ')')
227
    pobyso_autoprint(funcAuxSo)
228
    # Clear the Sollya range before a new call to infnorm and issue the call.
229
    sollya_lib_clear_obj(infnormSo)
230
    infnormSo = pobyso_infnorm_so_so(funcAuxSo,rangeSo)
231
    fMinSa = pobyso_get_interval_from_range_so_sa(infnormSo)
232
    sollya_lib_clear_obj(infnormSo)
233
    fMinLowerBoundSa = binadeTopLimitSa - fMinSa.lower()
234
    # Compute the maximum of the precisions of the different bounds.
235
    maxPrecSa = max([fMinLowerBoundSa.parent().precision(), \
236
                     fMaxUpperBoundSa.parent().precision()])
237
    # Create a RealIntervalField and create an interval with the "good" bounds.
238
    RRRI = RealIntervalField(maxPrecSa)
239
    imageIntervalSa = RRRI(fMinLowerBoundSa, fMaxUpperBoundSa)
240
    # Free the unneeded Sollya objects
241
    sollya_lib_clear_obj(funcSo)
242
    sollya_lib_clear_obj(funcAuxSo)
243
    sollya_lib_clear_obj(rangeSo)
244
    return(imageIntervalSa)
245
# End pobyso_compute_pos_function_abs_val_bounds_sa_sa
246

    
247
def pobyso_compute_precision_decay_ratio_function_sa_so():
248
    """
249
    Compute the precision decay ratio function for polynomial 
250
    coefficient progressive trucation.
251
    """
252
    functionText = """
253
    proc(deg, a, b, we, wq) 
254
    {
255
      k = we * (exp(x/a)-1) + wq * (b*x)^2 + (1-we-wq) * x;
256
      return k/k(d);
257
    };
258
    """
259
    return pobyso_parse_string_sa_so(functionText)
260
# End  pobyso_compute_precision_decay_ratio_function.
261

    
262

    
263
def pobyso_constant(rnArg):
264
    """ Legacy function. See pobyso_constant_sa_so. """
265
    return(pobyso_constant_sa_so(rnArg))
266
    
267
def pobyso_constant_sa_so(rnArgSa, precisionSa=None):
268
    """
269
    Create a Sollya constant from a Sage RealNumber.
270
    The sollya_lib_constant() function creates a constant
271
    with the same precision as the source.
272
    """
273
    ## Precision stuff. If one wants to change precisions,
274
    #  everything takes place in Sage. That only makes
275
    #  sense if one wants to reduce the precision.
276
    if not precisionSa is None:
277
        RRR = RealField(precisionSa)
278
        rnArgSa = RRR(rnArgSa)
279
    #print rnArgSa, rnArgSa.precision()
280
    # Sollya constant creation takes place here.
281
    return sollya_lib_constant(get_rn_value(rnArgSa))
282
# End pobyso_constant_sa_so
283
     
284
def pobyso_constant_0_sa_so():
285
    """
286
    Obvious.
287
    """
288
    return pobyso_constant_from_int_sa_so(0)
289

    
290
def pobyso_constant_1():
291
    """
292
    Obvious.
293
    Legacy function. See pobyso_constant_so_so. 
294
    """
295
    return pobyso_constant_1_sa_so()
296

    
297
def pobyso_constant_1_sa_so():
298
    """
299
    Obvious.
300
    """
301
    return(pobyso_constant_from_int_sa_so(1))
302

    
303
def pobyso_constant_from_int(anInt):
304
    """ Legacy function. See pobyso_constant_from_int_sa_so. """
305
    return pobyso_constant_from_int_sa_so(anInt)
306

    
307
def pobyso_constant_from_int_sa_so(anInt):
308
    """
309
    Get a Sollya constant from a Sage int.
310
    """
311
    return sollya_lib_constant_from_int64(long(anInt))
312

    
313
def pobyso_constant_from_int_so_sa(constSo):
314
    """
315
    Get a Sage int from a Sollya int constant.
316
    Usefull for precision or powers in polynomials.
317
    """
318
    constSa = c_long(0)
319
    sollya_lib_get_constant_as_int64(byref(constSa), constSo)
320
    return constSa.value
321
# End pobyso_constant_from_int_so_sa
322

    
323
def pobyso_constant_from_mpq_sa_so(rationalSa):
324
    """
325
    Make a Sollya constant from Sage rational.
326
    The Sollya constant is an unevaluated expression.
327
    Hence no precision argument is needed.
328
    It is better to leave this way since Sollya has its own
329
    optimized evaluation mecanism that tries very hard to
330
    return exact values or at least faithful ones.
331
    """
332
    ratExprSo = \
333
        sollya_lib_constant_from_mpq(sgmp_get_rational_value(rationalSa))
334
    return ratExprSo
335
# End pobyso_constant_from_mpq_sa_so.
336

    
337
def pobyso_constant_sollya_prec_sa_so(rnArgSa):
338
    """
339
    Create a Sollya constant from a Sage RealNumber at the
340
    current precision in Sollya.
341
    """
342
    currentSollyaPrecSa = pobyso_get_prec_so_sa()
343
    return pobyso_constant_sa_so(rnArgSa, currentSollyaPrecSa)
344
# End pobyso_constant_sollya_prec_sa_so
345

    
346
def pobyso_end_elliptic_list_so_sa_so(objectsListSo, intCountSa):
347
    """
348
    Create a Sollya end elliptic list made of the objectListSo[0] to
349
     objectsListSo[intCountSa-1] objects.
350
    """     
351
    return sollya_lib_end_elliptic_list(objectSo, int(intCountSa))
352

    
353
def pobyso_error_so():
354
    return sollya_lib_error(None)
355
# End pobyso_error().
356

    
357
def pobyso_evaluate_so_so(funcSo, argumentSo):
358
    """
359
    Evaluates funcSo for arguemntSo through sollya_lib_evaluate().
360
    """
361
    return sollya_lib_evaluate(funcSo, argumentSo)
362
# End pobyso_evaluate_so_so.
363

    
364
def pobyso_float_poly_sa_so(polySa, precSa = None):
365
    """
366
    Create a Sollya polynomial from a Sage RealField polynomial.
367
    """
368
    ## TODO: filter arguments.
369
    ## Precision. If a precision is given, convert the polynomial
370
    #  into the right polynomial field. If not convert it straight
371
    #  to Sollya.
372
    sollyaPrecChanged = False 
373
    (curSollyaPrecSo, curSollyaPrecSa) = pobyso_get_prec_so_so_sa()
374
    if precSa is None:
375
        precSa = polySa.parent().base_ring().precision()
376
    if (precSa != curSollyaPrecSa):
377
        precSo = pobyso_constant_from_int(precSa)
378
        pobyso_set_prec_so_so(precSo)
379
        sollya_lib_clear_obj(precSo)
380
        sollyaPrecChanged = True    
381
    ## Get exponents and coefficients.
382
    exponentsSa     = polySa.exponents()
383
    coefficientsSa  = polySa.coefficients()
384
    ## Build the polynomial.
385
    polySo = None
386
    for coefficientSa, exponentSa in zip(coefficientsSa, exponentsSa):
387
        #print coefficientSa.n(prec=precSa), exponentSa
388
        coefficientSo = \
389
            pobyso_constant_sa_so(coefficientSa)
390
        #pobyso_autoprint(coefficientSo)
391
        exponentSo = \
392
            pobyso_constant_from_int_sa_so(exponentSa)
393
        #pobyso_autoprint(exponentSo)
394
        monomialSo = sollya_lib_build_function_pow(
395
                       sollya_lib_build_function_free_variable(),
396
                       exponentSo)
397
        polyTermSo = sollya_lib_build_function_mul(coefficientSo,
398
                                                       monomialSo)
399
        if polySo is None:
400
            polySo = polyTermSo
401
        else:
402
            polySo = sollya_lib_build_function_add(polySo, polyTermSo)
403
    if sollyaPrecChanged:
404
        pobyso_set_prec_so_so(curSollyaPrecSo)
405
        sollya_lib_clear_obj(curSollyaPrecSo)
406
    return polySo
407
# End pobyso_float_poly_sa_so    
408

    
409
def pobyso_float_poly_so_sa(polySo, realFieldSa=None):
410
    """
411
    Convert a Sollya polynomial into a Sage floating-point polynomial.
412
    If no realField is given, a RealField corresponding to the maximum 
413
    precision of the coefficients is internally computed.
414
    The real field is not returned but can be easily retrieved from 
415
    the polynomial itself.
416
    ALGORITHM:
417
    - (optional) compute the RealField of the coefficients;
418
    - convert the Sollya expression into a Sage expression;
419
    - convert the Sage expression into a Sage polynomial
420
    """    
421
    if realFieldSa is None:
422
        expressionPrecSa = pobyso_get_max_prec_of_exp_so_sa(polySo)
423
        #print "Maximum precision of Sollya polynomial coefficients:", expressionPrecSa
424
        realFieldSa      = RealField(expressionPrecSa)
425
    #print "Sollya expression before...",
426
    #pobyso_autoprint(polySo)
427

    
428
    expressionSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo,
429
                                                             realFieldSa)
430
    #print "...Sollya expression after."
431
    #pobyso_autoprint(polySo)
432
    polyVariableSa = expressionSa.variables()[0]
433
    polyRingSa     = realFieldSa[str(polyVariableSa)]
434
    #print polyRingSa
435
    # Do not use the polynomial(expressionSa, ring=polyRingSa) form!
436
    polynomialSa = polyRingSa(expressionSa)
437
    polyCoeffsListSa = polynomialSa.coefficients()
438
    #for coeff in polyCoeffsListSa:
439
    #    print coeff.abs().n()
440
    return polynomialSa
441
# End pobyso_float_poly_so_sa
442

    
443
def pobyso_free_variable():
444
    """
445
    Ultra thin wrapper around the sollya_lib_function_build_free_variable function.
446
    """
447
    return sollya_lib_build_function_free_variable()
448
    
449
def pobyso_function_type_as_string(funcType):
450
    """ Legacy function. See pobyso_function_type_as_string_so_sa. """
451
    return(pobyso_function_type_as_string_so_sa(funcType))
452

    
453
def pobyso_function_type_as_string_so_sa(funcType):
454
    """
455
    Numeric Sollya function codes -> Sage mathematical function names.
456
    Notice that pow -> ^ (a la Sage, not a la Python).
457
    """
458
    if funcType == SOLLYA_BASE_FUNC_ABS:
459
        return "abs"
460
    elif funcType == SOLLYA_BASE_FUNC_ACOS:
461
        return "arccos"
462
    elif funcType == SOLLYA_BASE_FUNC_ACOSH:
463
        return "arccosh"
464
    elif funcType == SOLLYA_BASE_FUNC_ADD:
465
        return "+"
466
    elif funcType == SOLLYA_BASE_FUNC_ASIN:
467
        return "arcsin"
468
    elif funcType == SOLLYA_BASE_FUNC_ASINH:
469
        return "arcsinh"
470
    elif funcType == SOLLYA_BASE_FUNC_ATAN:
471
        return "arctan"
472
    elif funcType == SOLLYA_BASE_FUNC_ATANH:
473
        return "arctanh"
474
    elif funcType == SOLLYA_BASE_FUNC_CEIL:
475
        return "ceil"
476
    elif funcType == SOLLYA_BASE_FUNC_CONSTANT:
477
        return "cte"
478
    elif funcType == SOLLYA_BASE_FUNC_COS:
479
        return "cos"
480
    elif funcType == SOLLYA_BASE_FUNC_COSH:
481
        return "cosh"
482
    elif funcType == SOLLYA_BASE_FUNC_DIV:
483
        return "/"
484
    elif funcType == SOLLYA_BASE_FUNC_DOUBLE:
485
        return "double"
486
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEDOUBLE:
487
        return "doubleDouble"
488
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEEXTENDED:
489
        return "doubleDxtended"
490
    elif funcType == SOLLYA_BASE_FUNC_ERF:
491
        return "erf"
492
    elif funcType == SOLLYA_BASE_FUNC_ERFC:
493
        return "erfc"
494
    elif funcType == SOLLYA_BASE_FUNC_EXP:
495
        return "exp"
496
    elif funcType == SOLLYA_BASE_FUNC_EXP_M1:
497
        return "expm1"
498
    elif funcType == SOLLYA_BASE_FUNC_FLOOR:
499
        return "floor"
500
    elif funcType == SOLLYA_BASE_FUNC_FREE_VARIABLE:
501
        return "freeVariable"
502
    elif funcType == SOLLYA_BASE_FUNC_HALFPRECISION:
503
        return "halfPrecision"
504
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYCONSTANT:
505
        return "libraryConstant"
506
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYFUNCTION:
507
        return "libraryFunction"
508
    elif funcType == SOLLYA_BASE_FUNC_LOG:
509
        return "log"
510
    elif funcType == SOLLYA_BASE_FUNC_LOG_10:
511
        return "log10"
512
    elif funcType == SOLLYA_BASE_FUNC_LOG_1P:
513
        return "log1p"
514
    elif funcType == SOLLYA_BASE_FUNC_LOG_2:
515
        return "log2"
516
    elif funcType == SOLLYA_BASE_FUNC_MUL:
517
        return "*"
518
    elif funcType == SOLLYA_BASE_FUNC_NEARESTINT:
519
        return "round"
520
    elif funcType == SOLLYA_BASE_FUNC_NEG:
521
        return "__neg__"
522
    elif funcType == SOLLYA_BASE_FUNC_PI:
523
        return "pi"
524
    elif funcType == SOLLYA_BASE_FUNC_POW:
525
        return "^"
526
    elif funcType == SOLLYA_BASE_FUNC_PROCEDUREFUNCTION:
527
        return "procedureFunction"
528
    elif funcType == SOLLYA_BASE_FUNC_QUAD:
529
        return "quad"
530
    elif funcType == SOLLYA_BASE_FUNC_SIN:
531
        return "sin"
532
    elif funcType == SOLLYA_BASE_FUNC_SINGLE:
533
        return "single"
534
    elif funcType == SOLLYA_BASE_FUNC_SINH:
535
        return "sinh"
536
    elif funcType == SOLLYA_BASE_FUNC_SQRT:
537
        return "sqrt"
538
    elif funcType == SOLLYA_BASE_FUNC_SUB:
539
        return "-"
540
    elif funcType == SOLLYA_BASE_FUNC_TAN:
541
        return "tan"
542
    elif funcType == SOLLYA_BASE_FUNC_TANH:
543
        return "tanh"
544
    elif funcType == SOLLYA_BASE_FUNC_TRIPLEDOUBLE:
545
        return "tripleDouble"
546
    else:
547
        return None
548

    
549
def pobyso_get_constant(rnArgSa, constSo):
550
    """ Legacy function. See pobyso_get_constant_so_sa. """
551
    return pobyso_get_constant_so_sa(rnArgSa, constSo)
552
# End pobyso_get_constant
553

    
554
def pobyso_get_constant_so_sa(rnArgSa, constSo):
555
    """
556
    Set the value of rnArgSo to the value of constSo in MPFR_RNDN mode.
557
    rnArg must already exist and belong to some RealField.
558
    We assume that constSo points to a Sollya constant.
559
    """
560
    outcome = sollya_lib_get_constant(get_rn_value(rnArgSa), constSo)
561
    if outcome == 0: # Failure because constSo is not a constant expression.
562
        return None
563
    else:
564
        return outcome
565
# End  pobyso_get_constant_so_sa
566
   
567
def pobyso_get_constant_as_rn(ctExpSo):
568
    """ 
569
    Legacy function. See pobyso_get_constant_as_rn_so_sa. 
570
    """ 
571
    return(pobyso_get_constant_as_rn_so_sa(ctExpSo))
572
    
573
def pobyso_get_constant_as_rn_so_sa(constExpSo):
574
    """
575
    Get a Sollya constant as a Sage "real number".
576
    The precision of the floating-point number returned is that of the Sollya
577
    constant.
578
    """
579
    #print "Before computing precision of variable..."
580
    #pobyso_autoprint(constExpSo)
581
    precisionSa  = pobyso_get_prec_of_constant_so_sa(constExpSo)
582
    #print "precisionSa:", precisionSa
583
    ## If the expression can not be exactly converted, None is returned.
584
    #  In this case opt for the Sollya current expression.
585
    if precisionSa is None:
586
        precisionSa = pobyso_get_prec_so_sa()
587
    RRRR = RealField(precisionSa)
588
    rnSa = RRRR(0)
589
    outcome = sollya_lib_get_constant(get_rn_value(rnSa), constExpSo)
590
    if outcome == 0:
591
        return None
592
    else:
593
        return rnSa
594
# End pobyso_get_constant_as_rn_so_sa
595

    
596
def pobyso_get_constant_as_rn_with_rf(ctExp, realField):
597
    """ 
598
    Legacy function. See pobyso_get_constant_as_rn_with_rf_so_sa.
599
    """
600
    return pobyso_get_constant_as_rn_with_rf_so_sa(ctExp, realField)
601
# End pobyso_get_constant_as_rn_with_rf
602
    
603
def pobyso_get_constant_as_rn_with_rf_so_sa(ctExpSo, realFieldSa = None):
604
    """
605
    Get a Sollya constant as a Sage "real number".
606
    If no real field is specified, the precision of the floating-point number 
607
    returned is that of the Sollya constant.
608
    Otherwise is is that of the real field. Hence rounding may happen.
609
    """
610
    if realFieldSa is None:
611
        return pobyso_get_constant_as_rn_so_sa(ctExpSo)
612
    rnSa = realFieldSa(0)
613
    outcome = sollya_lib_get_constant(get_rn_value(rnSa), ctExpSo)
614
    if outcome == 0:
615
        return None
616
    else:
617
        return rnSa
618
# End pobyso_get_constant_as_rn_with_rf_so_sa
619

    
620
def pobyso_get_free_variable_name():
621
    """ 
622
    Legacy function. See pobyso_get_free_variable_name_so_sa.
623
    """
624
    return(pobyso_get_free_variable_name_so_sa())
625

    
626
def pobyso_get_free_variable_name_so_sa():
627
    return sollya_lib_get_free_variable_name()
628
    
629
def pobyso_get_function_arity(expressionSo):
630
    """ 
631
    Legacy function. See pobyso_get_function_arity_so_sa.
632
    """
633
    return(pobyso_get_function_arity_so_sa(expressionSo))
634

    
635
def pobyso_get_function_arity_so_sa(expressionSo):
636
    arity = c_int(0)
637
    sollya_lib_get_function_arity(byref(arity),expressionSo)
638
    return int(arity.value)
639

    
640
def pobyso_get_head_function(expressionSo):
641
    """ 
642
    Legacy function. See pobyso_get_head_function_so_sa. 
643
    """
644
    return(pobyso_get_head_function_so_sa(expressionSo)) 
645

    
646
def pobyso_get_head_function_so_sa(expressionSo):
647
    functionType = c_int(0)
648
    sollya_lib_get_head_function(byref(functionType), expressionSo)
649
    return int(functionType.value)
650

    
651
def pobyso_get_interval_from_range_so_sa(soRange, realIntervalFieldSa = None ):
652
    """
653
    Return the Sage interval corresponding to the Sollya range argument.
654
    If no reaIntervalField is passed as an argument, the interval bounds are not
655
    rounded: they are elements of RealIntervalField of the "right" precision
656
    to hold all the digits.
657
    """
658
    prec = c_int(0)
659
    if realIntervalFieldSa is None:
660
        retval = sollya_lib_get_prec_of_range(byref(prec), soRange, None)
661
        if retval == 0:
662
            return None
663
        realIntervalFieldSa = RealIntervalField(prec.value)
664
    intervalSa = realIntervalFieldSa(0,0)
665
    retval = \
666
        sollya_lib_get_interval_from_range(get_interval_value(intervalSa),\
667
                                           soRange)
668
    if retval == 0:
669
        return None
670
    return intervalSa
671
# End pobyso_get_interval_from_range_so_sa
672

    
673
def pobyso_get_list_elements(soObj):
674
    """ Legacy function. See pobyso_get_list_elements_so_so. """
675
    return pobyso_get_list_elements_so_so(soObj)
676
 
677
def pobyso_get_list_elements_so_so(objectListSo):
678
    """
679
    Get the Sollya list elements as a Sage/Python array of Sollya objects.
680
    
681
    INPUT:
682
    - objectListSo: a Sollya list of Sollya objects.
683
    
684
    OUTPUT:
685
    - a Sage/Python tuple made of:
686
      - a Sage/Python list of Sollya objects,
687
      - a Sage/Python int holding the number of elements,
688
      - a Sage/Python int stating (!= 0) that the list is end-elliptic.
689
    NOTE::
690
        We recover the addresses of the Sollya object from the list of pointers
691
        returned by sollya_lib_get_list_elements. The list itself is freed.
692
    TODO::
693
        Figure out what to do with numElements since the number of elements
694
        can easily be recovered from the list itself. 
695
        Ditto for isEndElliptic.
696
    """
697
    listAddress = POINTER(c_longlong)()
698
    numElements = c_int(0)
699
    isEndElliptic = c_int(0)
700
    listAsSageList = []
701
    result = sollya_lib_get_list_elements(byref(listAddress),\
702
                                          byref(numElements),\
703
                                          byref(isEndElliptic),\
704
                                          objectListSo)
705
    if result == 0 :
706
        return None
707
    for i in xrange(0, numElements.value, 1):
708
       #listAsSageList.append(sollya_lib_copy_obj(listAddress[i]))
709
       listAsSageList.append(listAddress[i])
710
       # Clear each of the elements returned by Sollya.
711
       #sollya_lib_clear_obj(listAddress[i])
712
    # Free the list itself.   
713
    sollya_lib_free(listAddress)
714
    return (listAsSageList, numElements.value, isEndElliptic.value)
715

    
716
def pobyso_get_max_prec_of_exp(soExp):
717
    """ Legacy function. See pobyso_get_max_prec_of_exp_so_sa. """
718
    return pobyso_get_max_prec_of_exp_so_sa(soExp)
719

    
720
def pobyso_get_max_prec_of_exp_so_sa(expSo):
721
    """
722
    Get the maximum precision used for the numbers in a Sollya expression.
723
    
724
    Arguments:
725
    soExp -- a Sollya expression pointer
726
    Return value:
727
    A Python integer
728
    TODO: 
729
    - error management;
730
    - correctly deal with numerical type such as DOUBLEEXTENDED.
731
    """
732
    maxPrecision = 0
733
    minConstPrec = 0
734
    currentConstPrec = 0
735
    operator = pobyso_get_head_function_so_sa(expSo)
736
    if (operator != SOLLYA_BASE_FUNC_CONSTANT) and \
737
    (operator != SOLLYA_BASE_FUNC_FREE_VARIABLE):
738
        (arity, subexpressions) = pobyso_get_subfunctions_so_sa(expSo)
739
        for i in xrange(arity):
740
            maxPrecisionCandidate = \
741
                pobyso_get_max_prec_of_exp_so_sa(subexpressions[i])
742
            if maxPrecisionCandidate > maxPrecision:
743
                maxPrecision = maxPrecisionCandidate
744
        return maxPrecision
745
    elif operator == SOLLYA_BASE_FUNC_CONSTANT:
746
        #minConstPrec = pobyso_get_min_prec_of_constant_so_sa(expSo)
747
        #currentConstPrec = pobyso_get_min_prec_of_constant_so_sa(soExp)
748
        #print minConstPrec, " - ", currentConstPrec 
749
        return pobyso_get_min_prec_of_constant_so_sa(expSo)
750
    
751
    elif operator == SOLLYA_BASE_FUNC_FREE_VARIABLE:
752
        return 0
753
    else:
754
        print "pobyso_get_max_prec_of_exp_so_sa: unexepected operator."
755
        return 0
756

    
757
def pobyso_get_min_prec_of_constant_so_sa(constExpSo):
758
    """
759
    Get the minimum precision necessary to represent the value of a Sollya
760
    constant.
761
    MPFR_MIN_PREC and powers of 2 are taken into account.
762
    We assume that constExpSo is a pointer to a Sollay constant expression.
763
    """
764
    constExpAsRnSa = pobyso_get_constant_as_rn_so_sa(constExpSo)
765
    return(min_mpfr_size(get_rn_value(constExpAsRnSa)))
766

    
767
def pobyso_get_poly_so_sa(polySo, realFieldSa=None):
768
    """
769
    Convert a Sollya polynomial into a Sage polynomial.
770
    Legacy function. Use pobyso_float_poly_so_sa() instead.
771
    """    
772
    return pobyso_float_poly_so_sa(polySo,realFieldSa)
773
# End pobyso_get_poly_so_sa
774

    
775
def pobyso_get_prec():
776
    """ Legacy function. See pobyso_get_prec_so_sa(). """
777
    return pobyso_get_prec_so_sa()
778

    
779
def pobyso_get_prec_so():
780
    """
781
    Get the current default precision in Sollya.
782
    The return value is a Sollya object.
783
    Usefull when modifying the precision back and forth by avoiding
784
    extra conversions.
785
    """
786
    return sollya_lib_get_prec(None)
787
    
788
def pobyso_get_prec_so_sa():
789
    """
790
    Get the current default precision in Sollya.
791
    The return value is Sage/Python int.
792
    """
793
    precSo = sollya_lib_get_prec(None)
794
    precSa = c_int(0)
795
    sollya_lib_get_constant_as_int(byref(precSa), precSo)
796
    sollya_lib_clear_obj(precSo)
797
    return int(precSa.value)
798
# End pobyso_get_prec_so_sa.
799

    
800
def pobyso_get_prec_so_so_sa():
801
    """
802
    Return the current precision both as a Sollya object and a
803
    Sage integer as hybrid tuple.
804
    To avoid multiple calls for precision manipulations. 
805
    """
806
    precSo = sollya_lib_get_prec(None)
807
    precSa = c_int(0)
808
    sollya_lib_get_constant_as_int(byref(precSa), precSo)
809
    return (precSo, precSa)
810
    
811
def pobyso_get_prec_of_constant(ctExpSo):
812
    """ Legacy function. See pobyso_get_prec_of_constant_so_sa. """
813
    return pobyso_get_prec_of_constant_so_sa(ctExpSo)
814

    
815
def pobyso_get_prec_of_constant_so_sa(ctExpSo):
816
    """
817
    Tries to find a precision to represent ctExpSo without rounding.
818
    If not possible, returns None.
819
    """
820
    #print "Entering pobyso_get_prec_of_constant_so_sa..."
821
    prec = c_int(0)
822
    retc = sollya_lib_get_prec_of_constant(byref(prec), ctExpSo, None)
823
    if retc == 0:
824
        #print "pobyso_get_prec_of_constant_so_sa failed."
825
        return None
826
    #print "...exiting pobyso_get_prec_of_constant_so_sa."
827
    return int(prec.value)
828

    
829
def pobyso_get_prec_of_range_so_sa(rangeSo):
830
    """
831
    Returns the number of bits elements of a range are coded with.
832
    """
833
    prec = c_int(0)
834
    retc = sollya_lib_get_prec_of_range(byref(prec), rangeSo, None)
835
    if retc == 0:
836
        return(None)
837
    return int(prec.value)
838
# End pobyso_get_prec_of_range_so_sa()
839

    
840
def pobyso_get_sage_exp_from_sollya_exp(sollyaExpSo, realField = RR):
841
    """ Legacy function. See pobyso_get_sage_exp_from_sollya_exp_so_sa. """
842
    return pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, 
843
                                                     realField = RR)
844

    
845
def pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, realFieldSa = RR):
846
    """
847
    Get a Sage expression from a Sollya expression. 
848
    Currently only tested with polynomials with floating-point coefficients.
849
    Notice that, in the returned polynomial, the exponents are RealNumbers.
850
    """
851
    #pobyso_autoprint(sollyaExp)
852
    operatorSa = pobyso_get_head_function_so_sa(sollyaExpSo)
853
    sollyaLibFreeVariableName = sollya_lib_get_free_variable_name()
854
    ## Get rid of the "_"'s in "_x_", if any.
855
    sollyaLibFreeVariableName = re.sub('_', '', sollyaLibFreeVariableName)
856
    # Constants and the free variable are special cases.
857
    # All other operator are dealt with in the same way.
858
    if (operatorSa != SOLLYA_BASE_FUNC_CONSTANT) and \
859
       (operatorSa != SOLLYA_BASE_FUNC_FREE_VARIABLE):
860
        (aritySa, subexpressionsSa) = pobyso_get_subfunctions_so_sa(sollyaExpSo)
861
        if aritySa == 1:
862
            sageExpSa = eval(pobyso_function_type_as_string_so_sa(operatorSa) + \
863
            "(" + pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], \
864
            realFieldSa) + ")")
865
        elif aritySa == 2:
866
            # We do not get through the preprocessor.
867
            # The "^" operator is then a special case.
868
            if operatorSa == SOLLYA_BASE_FUNC_POW:
869
                operatorAsStringSa = "**"
870
            else:
871
                operatorAsStringSa = \
872
                    pobyso_function_type_as_string_so_sa(operatorSa)
873
            sageExpSa = \
874
              eval("pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], realFieldSa)"\
875
              + " " + operatorAsStringSa + " " + \
876
                   "pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[1], realFieldSa)")
877
        # We do not know yet how to deal with arity >= 3 
878
        # (is there any in Sollya anyway?).
879
        else:
880
            sageExpSa = eval('None')
881
        return sageExpSa
882
    elif operatorSa == SOLLYA_BASE_FUNC_CONSTANT:
883
        #print "This is a constant"
884
        return pobyso_get_constant_as_rn_with_rf_so_sa(sollyaExpSo, realFieldSa)
885
    elif operatorSa == SOLLYA_BASE_FUNC_FREE_VARIABLE:
886
        #print "This is the free variable"
887
        return eval(sollyaLibFreeVariableName)
888
    else:
889
        print "Unexpected"
890
        return eval('None')
891
# End pobyso_get_sage_exp_from_sollya_exp_so_sa
892

    
893

    
894
def pobyso_get_subfunctions(expressionSo):
895
    """ Legacy function. See pobyso_get_subfunctions_so_sa. """
896
    return pobyso_get_subfunctions_so_sa(expressionSo) 
897
# End pobyso_get_subfunctions.
898
 
899
def pobyso_get_subfunctions_so_sa(expressionSo):
900
    """
901
    Get the subfunctions of an expression.
902
    Return the number of subfunctions and the list of subfunctions addresses.
903
    S.T.: Could not figure out another way than that ugly list of declarations
904
    to recover the addresses of the subfunctions.
905
    We limit ourselves to arity 8 functions. 
906
    """
907
    subf0 = c_int(0)
908
    subf1 = c_int(0)
909
    subf2 = c_int(0)
910
    subf3 = c_int(0)
911
    subf4 = c_int(0)
912
    subf5 = c_int(0)
913
    subf6 = c_int(0)
914
    subf7 = c_int(0)
915
    subf8 = c_int(0)
916
    arity = c_int(0)
917
    nullPtr = POINTER(c_int)()
918
    sollya_lib_get_subfunctions(expressionSo, byref(arity), \
919
      byref(subf0), byref(subf1), byref(subf2), byref(subf3), \
920
      byref(subf4), byref(subf5),\
921
      byref(subf6), byref(subf7), byref(subf8), nullPtr, None) 
922
#    byref(cast(subfunctions[0], POINTER(c_int))), \
923
#    byref(cast(subfunctions[0], POINTER(c_int))), \
924
#    byref(cast(subfunctions[2], POINTER(c_int))), \
925
#    byref(cast(subfunctions[3], POINTER(c_int))), \
926
#    byref(cast(subfunctions[4], POINTER(c_int))), \
927
#    byref(cast(subfunctions[5], POINTER(c_int))), \
928
#    byref(cast(subfunctions[6], POINTER(c_int))), \
929
#    byref(cast(subfunctions[7], POINTER(c_int))), \
930
#    byref(cast(subfunctions[8], POINTER(c_int))), nullPtr)
931
    subfunctions = [subf0, subf1, subf2, subf3, subf4, subf5, subf6, subf7, \
932
                    subf8]
933
    subs = []
934
    if arity.value > pobyso_max_arity:
935
        return(0,[])
936
    for i in xrange(arity.value):
937
        subs.append(int(subfunctions[i].value))
938
        #print subs[i]
939
    return (int(arity.value), subs)
940
# End pobyso_get_subfunctions_so_sa
941
    
942
def pobyso_guess_degree_sa_sa(functionSa, intervalSa, approxErrorSa, 
943
                              weightSa=None, degreeBoundSa=None):
944
    """
945
    Sa_sa variant of the solly_guessdegree function.
946
    Return 0 if something goes wrong.
947
    """
948
    functionAsStringSa = functionSa._assume_str().replace('_SAGE_VAR_', '')
949
    functionSo = pobyso_parse_string_sa_so(functionAsStringSa)
950
    if pobyso_is_error_so_sa(functionSo):
951
        sollya_lib_clear_obj(functionSo)
952
        return 0
953
    rangeSo = pobyso_interval_to_range_sa_so(intervalSa)
954
    # The approximation error is expected to be a floating point number.
955
    if pobyso_is_floating_point_number_sa_sa(approxErrorSa):
956
        approxErrorSo = pobyso_constant_sa_so(approxErrorSa)
957
    else:
958
        approxErrorSo = pobyso_constant_sa_so(RR(approxErrorSa))
959
    if not weightSa is None:
960
        weightAsStringSa = weightSa._assume_str().replace('_SAGE_VAR_', '')
961
        weightSo = pobyso_parse_string_sa_so(weightAsStringSa)
962
        if pobyso_is_error_so_sa(weightSo):
963
            sollya_lib_clear_obj(functionSo)
964
            sollya_lib_clear_obj(rangeSo)
965
            sollya_lib_clear_obj(approxErrorSo)   
966
            sollya_lib_clear_obj(weightSo)
967
            return 0   
968
    else:
969
        weightSo = None
970
    if not degreeBoundSa is None:
971
        degreeBoundSo = pobyso_constant_from_int_sa_so(degreeBoundSa)
972
    else:
973
        degreeBoundSo = None
974
    guessedDegreeSa = pobyso_guess_degree_so_sa(functionSo,
975
                                                rangeSo,
976
                                                approxErrorSo,
977
                                                weightSo,
978
                                                degreeBoundSo)
979
    sollya_lib_clear_obj(functionSo)
980
    sollya_lib_clear_obj(rangeSo)
981
    sollya_lib_clear_obj(approxErrorSo)
982
    if not weightSo is None:
983
        sollya_lib_clear_obj(weightSo)
984
    if not degreeBoundSo is None:
985
        sollya_lib_clear_obj(degreeBoundSo)
986
    return guessedDegreeSa
987
# End poyso_guess_degree_sa_sa
988

    
989
def pobyso_guess_degree_so_sa(functionSo, rangeSo, errorSo, weightSo=None, \
990
                              degreeBoundSo=None):
991
    """
992
    Thin wrapper around the guessdegree function.
993
    Nevertheless, some precision control stuff has been appended.
994
    """
995
    # Deal with Sollya internal precision issues: if it is too small, 
996
    # compared with the error, increases it to about twice -log2(error).
997
    errorSa = pobyso_get_constant_as_rn_with_rf_so_sa(errorSo)
998
    log2ErrorSa = errorSa.log2()
999
    if log2ErrorSa < 0:
1000
        neededPrecisionSa = int(2 * int(-log2ErrorSa) / 64) * 64
1001
    else:
1002
        neededPrecisionSa = int(2 * int(log2ErrorSa) / 64) * 64
1003
    #print "Needed precision:", neededPrecisionSa
1004
    currentPrecSa = pobyso_get_prec_so_sa()
1005
    if neededPrecisionSa > currentPrecSa:
1006
        currentPrecSo = pobyso_get_prec_so()
1007
        pobyso_set_prec_sa_so(neededPrecisionSa)
1008
    #print "Guessing degree..."
1009
    # weightSo and degreeBoundsSo are optional arguments.
1010
    # As declared, sollya_lib_guessdegree must take 5 arguments.
1011
    if weightSo is None:
1012
        degreeRangeSo = sollya_lib_guessdegree(functionSo, rangeSo, errorSo, 
1013
                                               0, 0, None)
1014
    elif degreeBoundSo is None:
1015
        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, \
1016
                                                errorSo, weightSo, 0, None)
1017
    else:
1018
        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, errorSo, \
1019
                                                weightSo, degreeBoundSo, None)
1020
    #print "...degree guess done."
1021
    # Restore internal precision, if applicable.
1022
    if neededPrecisionSa > currentPrecSa:
1023
        pobyso_set_prec_so_so(currentPrecSo)
1024
        sollya_lib_clear_obj(currentPrecSo)
1025
    degreeIntervalSa = pobyso_range_to_interval_so_sa(degreeRangeSo)
1026
    sollya_lib_clear_obj(degreeRangeSo)
1027
    # When ok, both bounds match.
1028
    # When the degree bound is too low, the upper bound is the degree
1029
    # for which the error can be honored.
1030
    # When it really goes wrong, the upper bound is infinity. 
1031
    if degreeIntervalSa.lower() == degreeIntervalSa.upper():
1032
        return int(degreeIntervalSa.lower())
1033
    else:
1034
        if degreeIntervalSa.upper().is_infinity():
1035
            return None
1036
        else:
1037
            return int(degreeIntervalSa.upper())
1038
    # End pobyso_guess_degree_so_sa
1039

    
1040
def pobyso_inf_so_so(intervalSo):
1041
    """
1042
    Very thin wrapper around sollya_lib_inf().
1043
    """
1044
    return sollya_lib_inf(intervalSo)
1045
# End pobyso_inf_so_so.
1046
    
1047
def pobyso_infnorm_so_so(func, interval, file = None, intervalList = None):
1048
    print "Do not use this function. User pobyso_supnorm_so_so instead."
1049
    return None
1050

    
1051
def pobyso_interval_to_range_sa_so(intervalSa, precisionSa=None):
1052
    if precisionSa is None:
1053
        precisionSa = intervalSa.parent().precision()
1054
    intervalSo = pobyso_bounds_to_range_sa_so(intervalSa.lower(),\
1055
                                              intervalSa.upper(),\
1056
                                              precisionSa)
1057
    return intervalSo
1058
# End pobyso_interval_to_range_sa_so
1059

    
1060
def pobyso_is_error_so_sa(objSo):
1061
    """
1062
    Thin wrapper around the sollya_lib_obj_is_error() function.
1063
    """
1064
    if sollya_lib_obj_is_error(objSo) != 0:
1065
        return True
1066
    else:
1067
        return False
1068
# End pobyso_is_error-so_sa
1069

    
1070
def pobyso_is_floating_point_number_sa_sa(numberSa):
1071
    """
1072
    Check whether a Sage number is floating point.
1073
    Exception stuff added because numbers other than
1074
    floating-point ones do not have the is_real() attribute.
1075
    """
1076
    try:
1077
        return numberSa.is_real()
1078
    except AttributeError:
1079
        return False
1080
# End pobyso_is_floating_piont_number_sa_sa
1081

    
1082
def pobyso_lib_init():
1083
    sollya_lib_init(None)
1084

    
1085
def pobyso_lib_close():
1086
    sollya_lib_close(None)
1087
    
1088
def pobyso_name_free_variable(freeVariableNameSa):
1089
    """ Legacy function. See pobyso_name_free_variable_sa_so. """
1090
    pobyso_name_free_variable_sa_so(freeVariableNameSa)
1091

    
1092
def pobyso_name_free_variable_sa_so(freeVariableNameSa):
1093
    """
1094
    Set the free variable name in Sollya from a Sage string.
1095
    """
1096
    sollya_lib_name_free_variable(freeVariableNameSa)
1097

    
1098
def pobyso_parse_string(string):
1099
    """ Legacy function. See pobyso_parse_string_sa_so. """
1100
    return pobyso_parse_string_sa_so(string)
1101
 
1102
def pobyso_parse_string_sa_so(string):
1103
    """
1104
    Get the Sollya expression computed from a Sage string or
1105
    a Sollya error object if parsing failed.
1106
    """
1107
    return sollya_lib_parse_string(string)
1108

    
1109
def pobyso_precision_so_sa(ctExpSo):
1110
    """
1111
    Computes the necessary precision to represent a number.
1112
    If x is not zero, it can be uniquely written as x = m · 2e 
1113
    where m is an odd integer and e is an integer. 
1114
    precision(x) returns the number of bits necessary to write m 
1115
    in binary (i.e. ceil(log2(m))).
1116
    """
1117
    #TODO: take care of the special case: 0, @NaN@, @Inf@
1118
    precisionSo = sollya_lib_precision(ctExpSo)
1119
    precisionSa = pobyso_constant_from_int_so_sa(precisionSo)
1120
    sollya_lib_clear_obj(precisionSo)
1121
    return precisionSa
1122
# End pobyso_precision_so_sa
1123

    
1124
def pobyso_polynomial_coefficients_progressive_round_so_so(polySo,
1125
                                                           funcSo,
1126
                                                           icSo,
1127
                                                           intervalSo,
1128
                                                           itpSo,
1129
                                                           ftpSo,
1130
                                                           maxPrecSo,
1131
                                                           maxErrSo,
1132
                                                           debug=False):
1133
    if debug:
1134
        print "Input arguments:"
1135
        pobyso_autoprint(polySo)
1136
        pobyso_autoprint(funcSo)
1137
        pobyso_autoprint(icSo)
1138
        pobyso_autoprint(intervalSo)
1139
        pobyso_autoprint(itpSo)
1140
        pobyso_autoprint(ftpSo)
1141
        pobyso_autoprint(maxPrecSo)
1142
        pobyso_autoprint(maxErrSo)
1143
        print "________________"
1144
    
1145
    ## Higher order function see: 
1146
    #  http://effbot.org/pyfaq/how-do-you-make-a-higher-order-function-in-python.htm
1147
    def precision_decay_ratio_function(degreeSa):
1148
        def outer(x):
1149
            def inner(x):
1150
                we = 3/8
1151
                wq = 2/8
1152
                a  = 2.2
1153
                b  = 2 
1154
                return we*(exp(x/a)-1) +  wq*((b*x)**2) + (1-we-wq)*x
1155
            return  inner(x)/inner(degreeSa)
1156
        return outer
1157
    
1158
    #
1159
    degreeSa        = pobyso_polynomial_degree_so_sa(polySo)
1160
    ratio           = precision_decay_ratio_function(degreeSa)
1161
    itpSa           = pobyso_constant_from_int_so_sa(itpSo)
1162
    ftpSa           = pobyso_constant_from_int_so_sa(ftpSo)
1163
    maxPrecSa       = pobyso_constant_from_int_so_sa(maxPrecSo)
1164
    maxErrSa        = pobyso_get_constant_as_rn_so_sa(maxErrSo)
1165
    if debug:
1166
        print "degreeSa:", degreeSa
1167
        print "ratio:", ratio
1168
        print "itpsSa:", itpSa
1169
        print "ftpSa:", ftpSa
1170
        print "maxPrecSa:", maxPrecSa
1171
        print "maxErrSa:", maxErrSa
1172
    lastResPolySo   = None
1173
    lastInfNormSo   = None
1174
    #print "About to enter the while loop..."
1175
    while True: 
1176
        resPolySo   = pobyso_constant_0_sa_so()
1177
        pDeltaSa    = ftpSa - itpSa
1178
        for indexSa in reversed(xrange(0,degreeSa+1)):
1179
            #print "Index:", indexSa
1180
            indexSo = pobyso_constant_from_int_sa_so(indexSa)
1181
            #pobyso_autoprint(indexSo)
1182
            #print ratio(indexSa)
1183
            ctpSa = floor(ftpSa - (pDeltaSa * ratio(indexSa)))
1184
            ctpSo = pobyso_constant_from_int_sa_so(ctpSa)
1185
            if debug:
1186
                print "Index:", indexSa, " - Target precision:", 
1187
                pobyso_autoprint(ctpSo)
1188
            cmonSo  = \
1189
                sollya_lib_build_function_mul(sollya_lib_coeff(polySo, indexSo),
1190
                                      sollya_lib_build_function_pow( \
1191
                                          sollya_lib_build_function_free_variable(), \
1192
                                          indexSo))
1193
            #pobyso_autoprint(cmonSo)
1194
            cmonrSo = pobyso_round_coefficients_single_so_so(cmonSo, ctpSo)
1195
            sollya_lib_clear_obj(cmonSo)
1196
            #pobyso_autoprint(cmonrSo)
1197
            resPolySo = sollya_lib_build_function_add(resPolySo,
1198
                                                      cmonrSo)
1199
            #pobyso_autoprint(resPolySo)
1200
        # End for index
1201
        freeVarSo     = sollya_lib_build_function_free_variable()
1202
        changeVarSo   = sollya_lib_sub(freeVarSo, icSo)
1203
        resPolyCvSo   = sollya_lib_evaluate(resPolySo, changeVarSo)
1204
        errFuncSo = sollya_lib_build_function_sub(sollya_lib_copy_obj(funcSo), 
1205
                                                  resPolyCvSo)
1206
        infNormSo = sollya_lib_dirtyinfnorm(errFuncSo, intervalSo)
1207
        cerrSa    = pobyso_get_constant_as_rn_so_sa(infNormSo)
1208
        if debug:
1209
            print "Infnorm (Sollya):", 
1210
            pobyso_autoprint(infNormSo)
1211
        sollya_lib_clear_obj(errFuncSo)
1212
        #print "Infnorm  (Sage):", cerrSa
1213
        if (cerrSa > maxErrSa):
1214
            if debug:
1215
                print "Error is too large."
1216
            if lastResPolySo is None:
1217
                if debug:
1218
                    print "Enlarging prec."
1219
                ntpSa = floor(ftpSa + ftpSa/50)
1220
                ## Can't enlarge (numerical)
1221
                if ntpSa == ftpSa:
1222
                    sollya_lib_clear_obj(resPolySo)
1223
                    return None
1224
                ## Can't enlarge (not enough precision left)
1225
                if ntpSa > maxPrecSa:
1226
                    sollya_lib_clear_obj(resPolySo)
1227
                    return None
1228
                ftpSa = ntpSa
1229
                continue
1230
            ## One enlargement took place.
1231
            else:
1232
                if debug:
1233
                    print "Exit with the last before last polynomial."
1234
                    print "Precision of highest degree monomial:", itpSa
1235
                    print "Precision of constant term          :", ftpSa
1236
                sollya_lib_clear_obj(resPolySo)
1237
                sollya_lib_clear_obj(infNormSo)
1238
                return (lastResPolySo, lastInfNormSo)
1239
        # cerrSa <= maxErrSa: scrap more bits, possibly.
1240
        else:
1241
            if debug:
1242
                print "Error is too small"
1243
            if cerrSa <= (maxErrSa/2):
1244
                if debug: 
1245
                    print "Shrinking prec."
1246
                ntpSa = floor(ftpSa - ftpSa/50)
1247
                ## Can't shrink (numerical)
1248
                if ntpSa == ftpSa:
1249
                    if not lastResPolySo is None:
1250
                        sollya_lib_clear_obj(lastResPolySo)
1251
                    if not lastInfNormSo is None:
1252
                        sollya_lib_clear_obj(lastInfNormSo)
1253
                    if debug:
1254
                        print "Exit because can't shrink anymore (numerically)."
1255
                        print "Precision of highest degree monomial:", itpSa
1256
                        print "Precision of constant term          :", ftpSa
1257
                    return (resPolySo, infNormSo)
1258
                ## Can't shrink (not enough precision left)
1259
                if ntpSa <= itpSa:
1260
                    if not lastResPolySo is None:
1261
                        sollya_lib_clear_obj(lastResPolySo)
1262
                    if not lastInfNormSo is None:
1263
                        sollya_lib_clear_obj(lastInfNormSo)
1264
                        print "Exit because can't shrink anymore (no bits left)."
1265
                        print "Precision of highest degree monomial:", itpSa
1266
                        print "Precision of constant term          :", ftpSa
1267
                    return (resPolySo, infNormSo)
1268
                ftpSa = ntpSa
1269
                if not lastResPolySo is None:
1270
                    sollya_lib_clear_obj(lastResPolySo)
1271
                if not lastInfNormSo is None:
1272
                    sollya_lib_clear_obj(lastInfNormSo)
1273
                lastResPolySo = resPolySo
1274
                lastInfNormSo = infNormSo
1275
                continue
1276
            else: # Error is not that small, just return
1277
                if not lastResPolySo is None:
1278
                    sollya_lib_clear_obj(lastResPolySo)
1279
                if not lastInfNormSo is None:
1280
                    sollya_lib_clear_obj(lastInfNormSo)
1281
                if debug:
1282
                    print "Exit normally."
1283
                    print "Precision of highest degree monomial:", itpSa
1284
                    print "Precision of constant term          :", ftpSa
1285
                return (resPolySo, infNormSo)                
1286
    # End wile True
1287
    return None
1288
# End pobyso_polynomial_coefficients_progressive_truncate_so_so.
1289

    
1290
def pobyso_polynomial_degree_so_sa(polySo):
1291
    """
1292
    Return the degree of a Sollya polynomial as a Sage int.
1293
    """
1294
    degreeSo = sollya_lib_degree(polySo)
1295
    return pobyso_constant_from_int_so_sa(degreeSo)
1296
# End pobyso_polynomial_degree_so_sa
1297

    
1298
def pobyso_polynomial_degree_so_so(polySo):
1299
    """
1300
    Thin wrapper around lib_sollya_degree().
1301
    """
1302
    return sollya_lib_degree(polySo)
1303
# End pobyso_polynomial_degree_so_so
1304
                                                                  
1305
def pobyso_range(rnLowerBound, rnUpperBound):
1306
    """ Legacy function. See pobyso_range_sa_so. """
1307
    return pobyso_range_sa_so(rnLowerBound, rnUpperBound) 
1308

    
1309

    
1310
def pobyso_range_to_interval_so_sa(rangeSo, realIntervalFieldSa = None):
1311
    """
1312
    Get a Sage interval from a Sollya range.
1313
    If no realIntervalField is given as a parameter, the Sage interval
1314
    precision is that of the Sollya range.
1315
    Otherwise, the precision is that of the realIntervalField. In this case
1316
    rounding may happen.
1317
    """
1318
    if realIntervalFieldSa is None:
1319
        precSa = pobyso_get_prec_of_range_so_sa(rangeSo)
1320
        realIntervalFieldSa = RealIntervalField(precSa)
1321
    intervalSa = \
1322
        pobyso_get_interval_from_range_so_sa(rangeSo, realIntervalFieldSa) 
1323
    return intervalSa
1324
# End pobyso_range_to_interval_so_sa
1325

    
1326
def pobyso_rat_poly_sa_so(polySa, precSa = None):
1327
    """
1328
    Create a Sollya polynomial from a Sage rational polynomial.
1329
    """
1330
    ## TODO: filter arguments.
1331
    ## Precision. If no precision is given, use the current precision
1332
    #  of Sollya.
1333
    if precSa is None:
1334
        precSa =  pobyso_get_prec_so_sa()
1335
    #print "Precision:",  precSa
1336
    RRR = RealField(precSa)
1337
    ## Create a Sage polynomial in the "right" precision.
1338
    P_RRR = RRR[polySa.variables()[0]]
1339
    polyFloatSa = P_RRR(polySa)
1340
    ## Make sure no precision is provided: pobyso_float_poly_sa_so will
1341
    #  recover it all by itself and not make an extra conversion.
1342
    return pobyso_float_poly_sa_so(polyFloatSa)
1343
    
1344
# End pobyso_rat_poly_sa_so    
1345
    
1346
def pobyso_remez_canonical_sa_sa(func, \
1347
                                 degree, \
1348
                                 lowerBound, \
1349
                                 upperBound, \
1350
                                 weight = None, \
1351
                                 quality = None):
1352
    """
1353
    All arguments are Sage/Python.
1354
    The functions (func and weight) must be passed as expressions or strings.
1355
    Otherwise the function fails. 
1356
    The return value is a Sage polynomial.
1357
    """
1358
    var('zorglub')    # Dummy variable name for type check only. Type of 
1359
    # zorglub is "symbolic expression".
1360
    polySo = pobyso_remez_canonical_sa_so(func, \
1361
                                 degree, \
1362
                                 lowerBound, \
1363
                                 upperBound, \
1364
                                 weight, \
1365
                                 quality)
1366
    # String test
1367
    if parent(func) == parent("string"):
1368
        functionSa = eval(func)
1369
    # Expression test.
1370
    elif type(func) == type(zorglub):
1371
        functionSa = func
1372
    else:
1373
        return None
1374
    #
1375
    maxPrecision = 0
1376
    if polySo is None:
1377
        return(None)
1378
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
1379
    RRRRSa = RealField(maxPrecision)
1380
    polynomialRingSa = RRRRSa[functionSa.variables()[0]]
1381
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, RRRRSa)
1382
    polySa = polynomial(expSa, polynomialRingSa)
1383
    sollya_lib_clear_obj(polySo)
1384
    return(polySa)
1385
# End pobyso_remez_canonical_sa_sa
1386
    
1387
def pobyso_remez_canonical(func, \
1388
                           degree, \
1389
                           lowerBound, \
1390
                           upperBound, \
1391
                           weight = "1", \
1392
                           quality = None):
1393
    """ Legacy function. See pobyso_remez_canonical_sa_so. """
1394
    return(pobyso_remez_canonical_sa_so(func, \
1395
                                        degree, \
1396
                                        lowerBound, \
1397
                                        upperBound, \
1398
                                        weight, \
1399
                                        quality))
1400
# End pobyso_remez_canonical.
1401

    
1402
def pobyso_remez_canonical_sa_so(func, \
1403
                                 degree, \
1404
                                 lowerBound, \
1405
                                 upperBound, \
1406
                                 weight = None, \
1407
                                 quality = None):
1408
    """
1409
    All arguments are Sage/Python.
1410
    The functions (func and weight) must be passed as expressions or strings.
1411
    Otherwise the function fails. 
1412
    The return value is a pointer to a Sollya function.
1413
    """
1414
    var('zorglub')    # Dummy variable name for type check only. Type of
1415
    # zorglub is "symbolic expression".
1416
    currentVariableNameSa = None
1417
    # The func argument can be of different types (string, 
1418
    # symbolic expression...)
1419
    if parent(func) == parent("string"):
1420
        localFuncSa = eval(func)
1421
        if len(localFuncSa.variables()) > 0:
1422
            currentVariableNameSa = localFuncSa.variables()[0]
1423
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1424
            functionSo = \
1425
              sollya_lib_parse_string(localFuncSa._assume_str().replace('_SAGE_VAR_', ''))
1426
    # Expression test.
1427
    elif type(func) == type(zorglub):
1428
        # Until we are able to translate Sage expressions into Sollya 
1429
        # expressions : parse the string version.
1430
        if len(func.variables()) > 0:
1431
            currentVariableNameSa = func.variables()[0]
1432
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1433
            functionSo = \
1434
              sollya_lib_parse_string(func._assume_str().replace('_SAGE_VAR_', ''))
1435
    else:
1436
        return(None)
1437
    if weight is None: # No weight given -> 1.
1438
        weightSo = pobyso_constant_1_sa_so()
1439
    elif parent(weight) == parent("string"): # Weight given as string: parse it.
1440
        weightSo = sollya_lib_parse_string(func)
1441
    elif type(weight) == type(zorglub): # Weight given as symbolice expression.
1442
        functionSo = \
1443
          sollya_lib_parse_string_sa_so(weight._assume_str().replace('_SAGE_VAR_', ''))
1444
    else:
1445
        return(None)
1446
    degreeSo = pobyso_constant_from_int(degree)
1447
    rangeSo = pobyso_bounds_to_range_sa_so(lowerBound, upperBound)
1448
    if not quality is None:
1449
        qualitySo= pobyso_constant_sa_so(quality)
1450
    else:
1451
        qualitySo = None
1452
        
1453
    remezPolySo = sollya_lib_remez(functionSo, \
1454
                                   degreeSo, \
1455
                                   rangeSo, \
1456
                                   weightSo, \
1457
                                   qualitySo, \
1458
                                   None)
1459
    sollya_lib_clear_obj(functionSo)
1460
    sollya_lib_clear_obj(degreeSo)
1461
    sollya_lib_clear_obj(rangeSo)
1462
    sollya_lib_clear_obj(weightSo)
1463
    if not qualitySo is None:
1464
        sollya_lib_clear_obj(qualitySo)
1465
    return(remezPolySo)
1466
# End pobyso_remez_canonical_sa_so
1467

    
1468
def pobyso_remez_canonical_so_so(funcSo, \
1469
                                 degreeSo, \
1470
                                 rangeSo, \
1471
                                 weightSo = pobyso_constant_1_sa_so(),\
1472
                                 qualitySo = None):
1473
    """
1474
    All arguments are pointers to Sollya objects.
1475
    The return value is a pointer to a Sollya function.
1476
    """
1477
    if not sollya_lib_obj_is_function(funcSo):
1478
        return(None)
1479
    return(sollya_lib_remez(funcSo, degreeSo, rangeSo, weightSo, qualitySo, None))
1480
# End pobyso_remez_canonical_so_so.
1481

    
1482
def pobyso_round_coefficients_single_so_so(polySo, precSo):
1483
    """
1484
    Create a rounded coefficients polynomial from polynomial argument to
1485
    the number of bits in size argument.
1486
    All coefficients are set to the same precision.
1487
    """
1488
    ## TODO: check arguments.
1489
    endEllipListSo = pobyso_build_end_elliptic_list_so_so(precSo)
1490
    polySo = sollya_lib_roundcoefficients(polySo, endEllipListSo, None)
1491
    sollya_lib_clear_obj(endEllipListSo)
1492
    #sollya_lib_clear_obj(endEllipListSo)
1493
    return polySo
1494
    
1495
# End pobyso_round_coefficients_single_so_so
1496

    
1497
def pobyso_set_canonical_off():
1498
    sollya_lib_set_canonical(sollya_lib_off())
1499

    
1500
def pobyso_set_canonical_on():
1501
    sollya_lib_set_canonical(sollya_lib_on())
1502

    
1503
def pobyso_set_prec(p):
1504
    """ Legacy function. See pobyso_set_prec_sa_so. """
1505
    pobyso_set_prec_sa_so(p)
1506

    
1507
def pobyso_set_prec_sa_so(p):
1508
    a = c_int(p)
1509
    precSo = c_void_p(sollya_lib_constant_from_int(a))
1510
    sollya_lib_set_prec(precSo, None)
1511
# End pobyso_set_prec_sa_so.
1512

    
1513
def pobyso_set_prec_so_so(newPrecSo):
1514
    sollya_lib_set_prec(newPrecSo, None)
1515
# End pobyso_set_prec_so_so.
1516

    
1517
def pobyso_inf_so_so(intervalSo):
1518
    """
1519
    Very thin wrapper around sollya_lib_inf().
1520
    """
1521
    return sollya_lib_inf(intervalSo)
1522
# End pobyso_inf_so_so.
1523
    
1524
def pobyso_supnorm_so_so(polySo, funcSo, intervalSo, errorTypeSo = None,\
1525
                         accuracySo = None):
1526
    """
1527
    Computes the supnorm of the approximation error between the given 
1528
    polynomial and function.
1529
    errorTypeSo defaults to "absolute".
1530
    accuracySo defaults to 2^(-40).
1531
    """
1532
    if errorTypeSo is None:
1533
        errorTypeSo = sollya_lib_absolute(None)
1534
        errorTypeIsNone = True
1535
    else:
1536
        errorTypeIsNone = False
1537
    #
1538
    if accuracySo is None:
1539
        # Notice the **!
1540
        accuracySo = pobyso_constant_sa_so(RR(2**(-40)))
1541
        accuracyIsNone = True
1542
    else:
1543
        accuracyIsNone = False
1544
    pobyso_autoprint(accuracySo)
1545
    resultSo = \
1546
        sollya_lib_supnorm(polySo, funcSo, intervalSo, errorTypeSo, \
1547
                              accuracySo)
1548
    if errorTypeIsNone:
1549
        sollya_lib_clear_obj(errorTypeSo)
1550
    if accuracyIsNone:
1551
        sollya_lib_clear_obj(accuracySo)
1552
    return resultSo
1553
# End pobyso_supnorm_so_so
1554

    
1555
def pobyso_taylor_expansion_no_change_var_so_so(functionSo, 
1556
                                                degreeSo, 
1557
                                                rangeSo,
1558
                                                errorTypeSo=None, 
1559
                                                sollyaPrecSo=None):
1560
    """
1561
    Compute the Taylor expansion without the variable change
1562
    x -> x-intervalCenter.
1563
    """
1564
    # No global change of the working precision.
1565
    if not sollyaPrecSo is None:
1566
        initialPrecSo = sollya_lib_get_prec(None)
1567
        sollya_lib_set_prec(sollyaPrecSo)
1568
    # Error type stuff: default to absolute.
1569
    if errorTypeSo is None:
1570
        errorTypeIsNone = True
1571
        errorTypeSo = sollya_lib_absolute(None)
1572
    else:
1573
        errorTypeIsNone = False
1574
    intervalCenterSo = sollya_lib_mid(rangeSo, None)
1575
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo,
1576
                                         intervalCenterSo,
1577
                                         rangeSo, errorTypeSo, None)
1578
    # taylorFormListSaSo is a Python list of Sollya objects references that 
1579
    # are copies of the elements of taylorFormSo.
1580
    # pobyso_get_list_elements_so_so clears taylorFormSo.
1581
    (taylorFormListSaSo, numElementsSa, isEndEllipticSa) = \
1582
        pobyso_get_list_elements_so_so(taylorFormSo)
1583
    polySo = sollya_lib_copy_obj(taylorFormListSaSo[0])
1584
    #print "Num elements:", numElementsSa
1585
    sollya_lib_clear_obj(taylorFormSo)
1586
    #polySo = taylorFormListSaSo[0]
1587
    #errorRangeSo = sollya_lib_copy_obj(taylorFormListSaSo[2])
1588
    errorRangeSo = taylorFormListSaSo[2]
1589
    # No copy_obj needed here: a new objects are created.
1590
    maxErrorSo    = sollya_lib_sup(errorRangeSo)
1591
    minErrorSo    = sollya_lib_inf(errorRangeSo)
1592
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
1593
    absMinErrorSo = sollya_lib_abs(minErrorSo)
1594
    sollya_lib_clear_obj(maxErrorSo)
1595
    sollya_lib_clear_obj(minErrorSo)
1596
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
1597
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
1598
    # If changed, reset the Sollya working precision.
1599
    if not sollyaPrecSo is None:
1600
        sollya_lib_set_prec(initialPrecSo)
1601
        sollya_lib_clear_obj(initialPrecSo)
1602
    if errorTypeIsNone:
1603
        sollya_lib_clear_obj(errorTypeSo)
1604
    pobyso_clear_taylorform_sa_so(taylorFormListSaSo)
1605
    if absMaxErrorSa > absMinErrorSa:
1606
        sollya_lib_clear_obj(absMinErrorSo)
1607
        return((polySo, intervalCenterSo, absMaxErrorSo))
1608
    else:
1609
        sollya_lib_clear_obj(absMaxErrorSo)
1610
        return((polySo, intervalCenterSo, absMinErrorSo))
1611
# end pobyso_taylor_expansion_no_change_var_so_so
1612

    
1613
def pobyso_taylor_expansion_with_change_var_so_so(functionSo, degreeSo, \
1614
                                                  rangeSo, \
1615
                                                  errorTypeSo=None, \
1616
                                                  sollyaPrecSo=None):
1617
    """
1618
    Compute the Taylor expansion with the variable change
1619
    x -> (x-intervalCenter) included.
1620
    """
1621
    # No global change of the working precision.
1622
    if not sollyaPrecSo is None:
1623
        initialPrecSo = sollya_lib_get_prec(None)
1624
        sollya_lib_set_prec(sollyaPrecSo)
1625
    #
1626
    # Error type stuff: default to absolute.
1627
    if errorTypeSo is None:
1628
        errorTypeIsNone = True
1629
        errorTypeSo = sollya_lib_absolute(None)
1630
    else:
1631
        errorTypeIsNone = False
1632
    intervalCenterSo = sollya_lib_mid(rangeSo)
1633
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo, \
1634
                                         intervalCenterSo, \
1635
                                         rangeSo, errorTypeSo, None)
1636
    # taylorFormListSaSo is a Python list of Sollya objects references that 
1637
    # are copies of the elements of taylorFormSo.
1638
    # pobyso_get_list_elements_so_so clears taylorFormSo.
1639
    (taylorFormListSo, numElements, isEndElliptic) = \
1640
        pobyso_get_list_elements_so_so(taylorFormSo)
1641
    polySo = taylorFormListSo[0]
1642
    errorRangeSo = taylorFormListSo[2]
1643
    maxErrorSo    = sollya_lib_sup(errorRangeSo)
1644
    minErrorSo    = sollya_lib_inf(errorRangeSo)
1645
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
1646
    absMinErrorSo = sollya_lib_abs(minErrorSo)
1647
    sollya_lib_clear_obj(maxErrorSo)
1648
    sollya_lib_clear_obj(minErrorSo)
1649
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
1650
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
1651
    changeVarExpSo = sollya_lib_build_function_sub(\
1652
                       sollya_lib_build_function_free_variable(),\
1653
                       sollya_lib_copy_obj(intervalCenterSo))
1654
    polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpSo)
1655
    sollya_lib_clear_obj(polySo) 
1656
    sollya_lib_clear_obj(changeVarExpSo)
1657
    # If changed, reset the Sollya working precision.
1658
    if not sollyaPrecSo is None:
1659
        sollya_lib_set_prec(initialPrecSo)
1660
        sollya_lib_clear_obj(initialPrecSo)
1661
    if errorTypeIsNone:
1662
        sollya_lib_clear_obj(errorTypeSo)
1663
    sollya_lib_clear_obj(taylorFormSo)
1664
    # Do not clear maxErrorSo.
1665
    if absMaxErrorSa > absMinErrorSa:
1666
        sollya_lib_clear_obj(absMinErrorSo)
1667
        return((polyVarChangedSo, intervalCenterSo, absMaxErrorSo))
1668
    else:
1669
        sollya_lib_clear_obj(absMaxErrorSo)
1670
        return((polyVarChangedSo, intervalCenterSo, absMinErrorSo))
1671
# end pobyso_taylor_expansion_with_change_var_so_so
1672

    
1673
def pobyso_taylor(function, degree, point):
1674
    """ Legacy function. See pobysoTaylor_so_so. """
1675
    return(pobyso_taylor_so_so(function, degree, point))
1676

    
1677
def pobyso_taylor_so_so(functionSo, degreeSo, pointSo):
1678
    return(sollya_lib_taylor(functionSo, degreeSo, pointSo))
1679
    
1680
def pobyso_taylorform(function, degree, point = None, 
1681
                      interval = None, errorType=None):
1682
    """ Legacy function. See pobyso_taylorform_sa_sa;"""
1683
    
1684
def pobyso_taylorform_sa_sa(functionSa, \
1685
                            degreeSa, \
1686
                            pointSa, \
1687
                            intervalSa=None, \
1688
                            errorTypeSa=None, \
1689
                            precisionSa=None):
1690
    """
1691
    Compute the Taylor form of 'degreeSa' for 'functionSa' at 'pointSa' 
1692
    for 'intervalSa' with 'errorTypeSa' (a string) using 'precisionSa'. 
1693
    point: must be a Real or a Real interval.
1694
    return the Taylor form as an array
1695
    TODO: take care of the interval and of the point when it is an interval;
1696
          when errorType is not None;
1697
          take care of the other elements of the Taylor form (coefficients 
1698
          errors and delta.
1699
    """
1700
    # Absolute as the default error.
1701
    if errorTypeSa is None:
1702
        errorTypeSo = sollya_lib_absolute()
1703
    elif errorTypeSa == "relative":
1704
        errorTypeSo = sollya_lib_relative()
1705
    elif errortypeSa == "absolute":
1706
        errorTypeSo = sollya_lib_absolute()
1707
    else:
1708
        # No clean up needed.
1709
        return None
1710
    # Global precision stuff
1711
    precisionChangedSa = False
1712
    currentSollyaPrecSo = pobyso_get_prec_so()
1713
    currentSollyaPrecSa = pobyso_constant_from_int_so_sa(currentSollyaPrecSo)
1714
    if not precisionSa is None:
1715
        if precisionSa > currentSollyaPrecSa:
1716
            pobyso_set_prec_sa_so(precisionSa)
1717
            precisionChangedSa = True
1718
            
1719
    if len(functionSa.variables()) > 0:
1720
        varSa = functionSa.variables()[0]
1721
        pobyso_name_free_variable_sa_so(str(varSa))
1722
    # In any case (point or interval) the parent of pointSa has a precision
1723
    # method.
1724
    pointPrecSa = pointSa.parent().precision()
1725
    if precisionSa > pointPrecSa:
1726
        pointPrecSa = precisionSa
1727
    # In any case (point or interval) pointSa has a base_ring() method.
1728
    pointBaseRingString = str(pointSa.base_ring())
1729
    if re.search('Interval', pointBaseRingString) is None: # Point
1730
        pointSo = pobyso_constant_sa_so(pointSa, pointPrecSa)
1731
    else: # Interval.
1732
        pointSo = pobyso_interval_to_range_sa_so(pointSa, pointPrecSa)
1733
    # Sollyafy the function.
1734
    functionSo = pobyso_parse_string_sa_so(functionSa._assume_str().replace('_SAGE_VAR_', ''))
1735
    if sollya_lib_obj_is_error(functionSo):
1736
        print "pobyso_tailorform: function string can't be parsed!"
1737
        return None
1738
    # Sollyafy the degree
1739
    degreeSo = sollya_lib_constant_from_int(int(degreeSa))
1740
    # Sollyafy the point
1741
    # Call Sollya
1742
    taylorFormSo = \
1743
        sollya_lib_taylorform(functionSo, degreeSo, pointSo, errorTypeSo,\
1744
                                         None)
1745
    sollya_lib_clear_obj(functionSo)
1746
    sollya_lib_clear_obj(degreeSo)
1747
    sollya_lib_clear_obj(pointSo)
1748
    sollya_lib_clear_obj(errorTypeSo)
1749
    (tfsAsList, numElements, isEndElliptic) = \
1750
            pobyso_get_list_elements_so_so(taylorFormSo)
1751
    polySo = tfsAsList[0]
1752
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
1753
    polyRealField = RealField(maxPrecision)
1754
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, polyRealField)
1755
    if precisionChangedSa:
1756
        sollya_lib_set_prec(currentSollyaPrecSo)
1757
        sollya_lib_clear_obj(currentSollyaPrecSo)
1758
    polynomialRing = polyRealField[str(varSa)]
1759
    polySa = polynomial(expSa, polynomialRing)
1760
    taylorFormSa = [polySa]
1761
    # Final clean-up.
1762
    sollya_lib_clear_obj(taylorFormSo)
1763
    return(taylorFormSa)
1764
# End pobyso_taylor_form_sa_sa
1765

    
1766
def pobyso_taylorform_so_so(functionSo, degreeSo, pointSo, intervalSo=None, \
1767
                            errorTypeSo=None):
1768
    createdErrorType = False
1769
    if errorTypeSo is None:
1770
        errorTypeSo = sollya_lib_absolute()
1771
        createdErrorType = True
1772
    else:
1773
        #TODO: deal with the other case.
1774
        pass
1775
    if intervalSo is None:
1776
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
1777
                                         errorTypeSo, None)
1778
    else:
1779
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
1780
                                         intervalSo, errorTypeSo, None)
1781
    if createdErrorType:
1782
        sollya_lib_clear_obj(errorTypeSo)
1783
    return resultSo
1784
        
1785

    
1786
def pobyso_univar_polynomial_print_reverse(polySa):
1787
    """ Legacy function. See pobyso_univar_polynomial_print_reverse_sa_sa. """
1788
    return(pobyso_univar_polynomial_print_reverse_sa_sa(polySa))
1789

    
1790
def pobyso_univar_polynomial_print_reverse_sa_sa(polySa):
1791
    """
1792
    Return the string representation of a univariate polynomial with
1793
    monomials ordered in the x^0..x^n order of the monomials.
1794
    Remember: Sage
1795
    """
1796
    polynomialRing = polySa.base_ring()
1797
    # A very expensive solution:
1798
    # -create a fake multivariate polynomial field with only one variable,
1799
    #   specifying a negative lexicographical order;
1800
    mpolynomialRing = PolynomialRing(polynomialRing.base(), \
1801
                                     polynomialRing.variable_name(), \
1802
                                     1, order='neglex')
1803
    # - convert the univariate argument polynomial into a multivariate
1804
    #   version;
1805
    p = mpolynomialRing(polySa)
1806
    # - return the string representation of the converted form.
1807
    # There is no simple str() method defined for p's class.
1808
    return(p.__str__())
1809
#
1810
print pobyso_get_prec()  
1811
pobyso_set_prec(165)
1812
print pobyso_get_prec()  
1813
a=100
1814
print type(a)
1815
id(a)
1816
print "Max arity: ", pobyso_max_arity
1817
print "Function tripleDouble (43) as a string: ", pobyso_function_type_as_string(43)
1818
print "Function None (44) as a string: ", pobyso_function_type_as_string(44)
1819
print "...Pobyso check done"