Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 218

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

    
1266
def pobyso_polynomial_degree_so_sa(polySo):
1267
    """
1268
    Return the degree of a Sollya polynomial as a Sage int.
1269
    """
1270
    degreeSo = sollya_lib_degree(polySo)
1271
    return pobyso_constant_from_int_so_sa(degreeSo)
1272
# End pobyso_polynomial_degree_so_sa
1273

    
1274
def pobyso_polynomial_degree_so_so(polySo):
1275
    """
1276
    Thin wrapper around lib_sollya_degree().
1277
    """
1278
    return sollya_lib_degree(polySo)
1279
# End pobyso_polynomial_degree_so_so
1280
                                                                  
1281
def pobyso_range(rnLowerBound, rnUpperBound):
1282
    """ Legacy function. See pobyso_range_sa_so. """
1283
    return pobyso_range_sa_so(rnLowerBound, rnUpperBound) 
1284

    
1285

    
1286
def pobyso_range_to_interval_so_sa(rangeSo, realIntervalFieldSa = None):
1287
    """
1288
    Get a Sage interval from a Sollya range.
1289
    If no realIntervalField is given as a parameter, the Sage interval
1290
    precision is that of the Sollya range.
1291
    Otherwise, the precision is that of the realIntervalField. In this case
1292
    rounding may happen.
1293
    """
1294
    if realIntervalFieldSa is None:
1295
        precSa = pobyso_get_prec_of_range_so_sa(rangeSo)
1296
        realIntervalFieldSa = RealIntervalField(precSa)
1297
    intervalSa = \
1298
        pobyso_get_interval_from_range_so_sa(rangeSo, realIntervalFieldSa) 
1299
    return intervalSa
1300
# End pobyso_range_to_interval_so_sa
1301

    
1302
def pobyso_rat_poly_sa_so(polySa, precSa = None):
1303
    """
1304
    Create a Sollya polynomial from a Sage rational polynomial.
1305
    """
1306
    ## TODO: filter arguments.
1307
    ## Precision. If no precision is given, use the current precision
1308
    #  of Sollya.
1309
    if precSa is None:
1310
        precSa =  pobyso_get_prec_so_sa()
1311
    #print "Precision:",  precSa
1312
    RRR = RealField(precSa)
1313
    ## Create a Sage polynomial in the "right" precision.
1314
    P_RRR = RRR[polySa.variables()[0]]
1315
    polyFloatSa = P_RRR(polySa)
1316
    ## Make sure no precision is provided: pobyso_float_poly_sa_so will
1317
    #  recover it all by itself and not make an extra conversion.
1318
    return pobyso_float_poly_sa_so(polyFloatSa)
1319
    
1320
# End pobyso_rat_poly_sa_so    
1321
    
1322
def pobyso_remez_canonical_sa_sa(func, \
1323
                                 degree, \
1324
                                 lowerBound, \
1325
                                 upperBound, \
1326
                                 weight = None, \
1327
                                 quality = None):
1328
    """
1329
    All arguments are Sage/Python.
1330
    The functions (func and weight) must be passed as expressions or strings.
1331
    Otherwise the function fails. 
1332
    The return value is a Sage polynomial.
1333
    """
1334
    var('zorglub')    # Dummy variable name for type check only. Type of 
1335
    # zorglub is "symbolic expression".
1336
    polySo = pobyso_remez_canonical_sa_so(func, \
1337
                                 degree, \
1338
                                 lowerBound, \
1339
                                 upperBound, \
1340
                                 weight, \
1341
                                 quality)
1342
    # String test
1343
    if parent(func) == parent("string"):
1344
        functionSa = eval(func)
1345
    # Expression test.
1346
    elif type(func) == type(zorglub):
1347
        functionSa = func
1348
    else:
1349
        return None
1350
    #
1351
    maxPrecision = 0
1352
    if polySo is None:
1353
        return(None)
1354
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
1355
    RRRRSa = RealField(maxPrecision)
1356
    polynomialRingSa = RRRRSa[functionSa.variables()[0]]
1357
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, RRRRSa)
1358
    polySa = polynomial(expSa, polynomialRingSa)
1359
    sollya_lib_clear_obj(polySo)
1360
    return(polySa)
1361
# End pobyso_remez_canonical_sa_sa
1362
    
1363
def pobyso_remez_canonical(func, \
1364
                           degree, \
1365
                           lowerBound, \
1366
                           upperBound, \
1367
                           weight = "1", \
1368
                           quality = None):
1369
    """ Legacy function. See pobyso_remez_canonical_sa_so. """
1370
    return(pobyso_remez_canonical_sa_so(func, \
1371
                                        degree, \
1372
                                        lowerBound, \
1373
                                        upperBound, \
1374
                                        weight, \
1375
                                        quality))
1376
# End pobyso_remez_canonical.
1377

    
1378
def pobyso_remez_canonical_sa_so(func, \
1379
                                 degree, \
1380
                                 lowerBound, \
1381
                                 upperBound, \
1382
                                 weight = None, \
1383
                                 quality = None):
1384
    """
1385
    All arguments are Sage/Python.
1386
    The functions (func and weight) must be passed as expressions or strings.
1387
    Otherwise the function fails. 
1388
    The return value is a pointer to a Sollya function.
1389
    """
1390
    var('zorglub')    # Dummy variable name for type check only. Type of
1391
    # zorglub is "symbolic expression".
1392
    currentVariableNameSa = None
1393
    # The func argument can be of different types (string, 
1394
    # symbolic expression...)
1395
    if parent(func) == parent("string"):
1396
        localFuncSa = eval(func)
1397
        if len(localFuncSa.variables()) > 0:
1398
            currentVariableNameSa = localFuncSa.variables()[0]
1399
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1400
            functionSo = \
1401
              sollya_lib_parse_string(localFuncSa._assume_str().replace('_SAGE_VAR_', ''))
1402
    # Expression test.
1403
    elif type(func) == type(zorglub):
1404
        # Until we are able to translate Sage expressions into Sollya 
1405
        # expressions : parse the string version.
1406
        if len(func.variables()) > 0:
1407
            currentVariableNameSa = func.variables()[0]
1408
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1409
            functionSo = \
1410
              sollya_lib_parse_string(func._assume_str().replace('_SAGE_VAR_', ''))
1411
    else:
1412
        return(None)
1413
    if weight is None: # No weight given -> 1.
1414
        weightSo = pobyso_constant_1_sa_so()
1415
    elif parent(weight) == parent("string"): # Weight given as string: parse it.
1416
        weightSo = sollya_lib_parse_string(func)
1417
    elif type(weight) == type(zorglub): # Weight given as symbolice expression.
1418
        functionSo = \
1419
          sollya_lib_parse_string_sa_so(weight._assume_str().replace('_SAGE_VAR_', ''))
1420
    else:
1421
        return(None)
1422
    degreeSo = pobyso_constant_from_int(degree)
1423
    rangeSo = pobyso_bounds_to_range_sa_so(lowerBound, upperBound)
1424
    if not quality is None:
1425
        qualitySo= pobyso_constant_sa_so(quality)
1426
    else:
1427
        qualitySo = None
1428
        
1429
    remezPolySo = sollya_lib_remez(functionSo, \
1430
                                   degreeSo, \
1431
                                   rangeSo, \
1432
                                   weightSo, \
1433
                                   qualitySo, \
1434
                                   None)
1435
    sollya_lib_clear_obj(functionSo)
1436
    sollya_lib_clear_obj(degreeSo)
1437
    sollya_lib_clear_obj(rangeSo)
1438
    sollya_lib_clear_obj(weightSo)
1439
    if not qualitySo is None:
1440
        sollya_lib_clear_obj(qualitySo)
1441
    return(remezPolySo)
1442
# End pobyso_remez_canonical_sa_so
1443

    
1444
def pobyso_remez_canonical_so_so(funcSo, \
1445
                                 degreeSo, \
1446
                                 rangeSo, \
1447
                                 weightSo = pobyso_constant_1_sa_so(),\
1448
                                 qualitySo = None):
1449
    """
1450
    All arguments are pointers to Sollya objects.
1451
    The return value is a pointer to a Sollya function.
1452
    """
1453
    if not sollya_lib_obj_is_function(funcSo):
1454
        return(None)
1455
    return(sollya_lib_remez(funcSo, degreeSo, rangeSo, weightSo, qualitySo, None))
1456
# End pobyso_remez_canonical_so_so.
1457

    
1458
def pobyso_round_coefficients_single_so_so(polySo, precSo):
1459
    """
1460
    Create a rounded coefficients polynomial from polynomial argument to
1461
    the number of bits in size argument.
1462
    All coefficients are set to the same precision.
1463
    """
1464
    ## TODO: check arguments.
1465
    endEllipListSo = pobyso_build_end_elliptic_list_so_so(precSo)
1466
    polySo = sollya_lib_roundcoefficients(polySo, endEllipListSo, None)
1467
    sollya_lib_clear_obj(endEllipListSo)
1468
    #sollya_lib_clear_obj(endEllipListSo)
1469
    return polySo
1470
    
1471
# End pobyso_round_coefficients_single_so_so
1472

    
1473
def pobyso_set_canonical_off():
1474
    sollya_lib_set_canonical(sollya_lib_off())
1475

    
1476
def pobyso_set_canonical_on():
1477
    sollya_lib_set_canonical(sollya_lib_on())
1478

    
1479
def pobyso_set_prec(p):
1480
    """ Legacy function. See pobyso_set_prec_sa_so. """
1481
    pobyso_set_prec_sa_so(p)
1482

    
1483
def pobyso_set_prec_sa_so(p):
1484
    a = c_int(p)
1485
    precSo = c_void_p(sollya_lib_constant_from_int(a))
1486
    sollya_lib_set_prec(precSo, None)
1487
# End pobyso_set_prec_sa_so.
1488

    
1489
def pobyso_set_prec_so_so(newPrecSo):
1490
    sollya_lib_set_prec(newPrecSo, None)
1491
# End pobyso_set_prec_so_so.
1492

    
1493
def pobyso_inf_so_so(intervalSo):
1494
    """
1495
    Very thin wrapper around sollya_lib_inf().
1496
    """
1497
    return sollya_lib_inf(intervalSo)
1498
# End pobyso_inf_so_so.
1499
    
1500
def pobyso_supnorm_so_so(polySo, funcSo, intervalSo, errorTypeSo = None,\
1501
                         accuracySo = None):
1502
    """
1503
    Computes the supnorm of the approximation error between the given 
1504
    polynomial and function.
1505
    errorTypeSo defaults to "absolute".
1506
    accuracySo defaults to 2^(-40).
1507
    """
1508
    if errorTypeSo is None:
1509
        errorTypeSo = sollya_lib_absolute(None)
1510
        errorTypeIsNone = True
1511
    else:
1512
        errorTypeIsNone = False
1513
    #
1514
    if accuracySo is None:
1515
        # Notice the **!
1516
        accuracySo = pobyso_constant_sa_so(RR(2**(-40)))
1517
        accuracyIsNone = True
1518
    else:
1519
        accuracyIsNone = False
1520
    pobyso_autoprint(accuracySo)
1521
    resultSo = \
1522
        sollya_lib_supnorm(polySo, funcSo, intervalSo, errorTypeSo, \
1523
                              accuracySo)
1524
    if errorTypeIsNone:
1525
        sollya_lib_clear_obj(errorTypeSo)
1526
    if accuracyIsNone:
1527
        sollya_lib_clear_obj(accuracySo)
1528
    return resultSo
1529
# End pobyso_supnorm_so_so
1530

    
1531
def pobyso_taylor_expansion_no_change_var_so_so(functionSo, 
1532
                                                degreeSo, 
1533
                                                rangeSo,
1534
                                                errorTypeSo=None, 
1535
                                                sollyaPrecSo=None):
1536
    """
1537
    Compute the Taylor expansion without the variable change
1538
    x -> x-intervalCenter.
1539
    """
1540
    # No global change of the working precision.
1541
    if not sollyaPrecSo is None:
1542
        initialPrecSo = sollya_lib_get_prec(None)
1543
        sollya_lib_set_prec(sollyaPrecSo)
1544
    # Error type stuff: default to absolute.
1545
    if errorTypeSo is None:
1546
        errorTypeIsNone = True
1547
        errorTypeSo = sollya_lib_absolute(None)
1548
    else:
1549
        errorTypeIsNone = False
1550
    intervalCenterSo = sollya_lib_mid(rangeSo, None)
1551
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo,
1552
                                         intervalCenterSo,
1553
                                         rangeSo, errorTypeSo, None)
1554
    # taylorFormListSaSo is a Python list of Sollya objects references that 
1555
    # are copies of the elements of taylorFormSo.
1556
    # pobyso_get_list_elements_so_so clears taylorFormSo.
1557
    (taylorFormListSaSo, numElementsSa, isEndEllipticSa) = \
1558
        pobyso_get_list_elements_so_so(taylorFormSo)
1559
    polySo = sollya_lib_copy_obj(taylorFormListSaSo[0])
1560
    #print "Num elements:", numElementsSa
1561
    sollya_lib_clear_obj(taylorFormSo)
1562
    #polySo = taylorFormListSaSo[0]
1563
    #errorRangeSo = sollya_lib_copy_obj(taylorFormListSaSo[2])
1564
    errorRangeSo = taylorFormListSaSo[2]
1565
    # No copy_obj needed here: a new objects are created.
1566
    maxErrorSo    = sollya_lib_sup(errorRangeSo)
1567
    minErrorSo    = sollya_lib_inf(errorRangeSo)
1568
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
1569
    absMinErrorSo = sollya_lib_abs(minErrorSo)
1570
    sollya_lib_clear_obj(maxErrorSo)
1571
    sollya_lib_clear_obj(minErrorSo)
1572
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
1573
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
1574
    # If changed, reset the Sollya working precision.
1575
    if not sollyaPrecSo is None:
1576
        sollya_lib_set_prec(initialPrecSo)
1577
        sollya_lib_clear_obj(initialPrecSo)
1578
    if errorTypeIsNone:
1579
        sollya_lib_clear_obj(errorTypeSo)
1580
    pobyso_clear_taylorform_sa_so(taylorFormListSaSo)
1581
    if absMaxErrorSa > absMinErrorSa:
1582
        sollya_lib_clear_obj(absMinErrorSo)
1583
        return((polySo, intervalCenterSo, absMaxErrorSo))
1584
    else:
1585
        sollya_lib_clear_obj(absMaxErrorSo)
1586
        return((polySo, intervalCenterSo, absMinErrorSo))
1587
# end pobyso_taylor_expansion_no_change_var_so_so
1588

    
1589
def pobyso_taylor_expansion_with_change_var_so_so(functionSo, degreeSo, \
1590
                                                  rangeSo, \
1591
                                                  errorTypeSo=None, \
1592
                                                  sollyaPrecSo=None):
1593
    """
1594
    Compute the Taylor expansion with the variable change
1595
    x -> (x-intervalCenter) included.
1596
    """
1597
    # No global change of the working precision.
1598
    if not sollyaPrecSo is None:
1599
        initialPrecSo = sollya_lib_get_prec(None)
1600
        sollya_lib_set_prec(sollyaPrecSo)
1601
    #
1602
    # Error type stuff: default to absolute.
1603
    if errorTypeSo is None:
1604
        errorTypeIsNone = True
1605
        errorTypeSo = sollya_lib_absolute(None)
1606
    else:
1607
        errorTypeIsNone = False
1608
    intervalCenterSo = sollya_lib_mid(rangeSo)
1609
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo, \
1610
                                         intervalCenterSo, \
1611
                                         rangeSo, errorTypeSo, None)
1612
    # taylorFormListSaSo is a Python list of Sollya objects references that 
1613
    # are copies of the elements of taylorFormSo.
1614
    # pobyso_get_list_elements_so_so clears taylorFormSo.
1615
    (taylorFormListSo, numElements, isEndElliptic) = \
1616
        pobyso_get_list_elements_so_so(taylorFormSo)
1617
    polySo = taylorFormListSo[0]
1618
    errorRangeSo = taylorFormListSo[2]
1619
    maxErrorSo    = sollya_lib_sup(errorRangeSo)
1620
    minErrorSo    = sollya_lib_inf(errorRangeSo)
1621
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
1622
    absMinErrorSo = sollya_lib_abs(minErrorSo)
1623
    sollya_lib_clear_obj(maxErrorSo)
1624
    sollya_lib_clear_obj(minErrorSo)
1625
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
1626
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
1627
    changeVarExpSo = sollya_lib_build_function_sub(\
1628
                       sollya_lib_build_function_free_variable(),\
1629
                       sollya_lib_copy_obj(intervalCenterSo))
1630
    polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpSo)
1631
    sollya_lib_clear_obj(polySo) 
1632
    sollya_lib_clear_obj(changeVarExpSo)
1633
    # If changed, reset the Sollya working precision.
1634
    if not sollyaPrecSo is None:
1635
        sollya_lib_set_prec(initialPrecSo)
1636
        sollya_lib_clear_obj(initialPrecSo)
1637
    if errorTypeIsNone:
1638
        sollya_lib_clear_obj(errorTypeSo)
1639
    sollya_lib_clear_obj(taylorFormSo)
1640
    # Do not clear maxErrorSo.
1641
    if absMaxErrorSa > absMinErrorSa:
1642
        sollya_lib_clear_obj(absMinErrorSo)
1643
        return((polyVarChangedSo, intervalCenterSo, absMaxErrorSo))
1644
    else:
1645
        sollya_lib_clear_obj(absMaxErrorSo)
1646
        return((polyVarChangedSo, intervalCenterSo, absMinErrorSo))
1647
# end pobyso_taylor_expansion_with_change_var_so_so
1648

    
1649
def pobyso_taylor(function, degree, point):
1650
    """ Legacy function. See pobysoTaylor_so_so. """
1651
    return(pobyso_taylor_so_so(function, degree, point))
1652

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

    
1742
def pobyso_taylorform_so_so(functionSo, degreeSo, pointSo, intervalSo=None, \
1743
                            errorTypeSo=None):
1744
    createdErrorType = False
1745
    if errorTypeSo is None:
1746
        errorTypeSo = sollya_lib_absolute()
1747
        createdErrorType = True
1748
    else:
1749
        #TODO: deal with the other case.
1750
        pass
1751
    if intervalSo is None:
1752
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
1753
                                         errorTypeSo, None)
1754
    else:
1755
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
1756
                                         intervalSo, errorTypeSo, None)
1757
    if createdErrorType:
1758
        sollya_lib_clear_obj(errorTypeSo)
1759
    return resultSo
1760
        
1761

    
1762
def pobyso_univar_polynomial_print_reverse(polySa):
1763
    """ Legacy function. See pobyso_univar_polynomial_print_reverse_sa_sa. """
1764
    return(pobyso_univar_polynomial_print_reverse_sa_sa(polySa))
1765

    
1766
def pobyso_univar_polynomial_print_reverse_sa_sa(polySa):
1767
    """
1768
    Return the string representation of a univariate polynomial with
1769
    monomials ordered in the x^0..x^n order of the monomials.
1770
    Remember: Sage
1771
    """
1772
    polynomialRing = polySa.base_ring()
1773
    # A very expensive solution:
1774
    # -create a fake multivariate polynomial field with only one variable,
1775
    #   specifying a negative lexicographical order;
1776
    mpolynomialRing = PolynomialRing(polynomialRing.base(), \
1777
                                     polynomialRing.variable_name(), \
1778
                                     1, order='neglex')
1779
    # - convert the univariate argument polynomial into a multivariate
1780
    #   version;
1781
    p = mpolynomialRing(polySa)
1782
    # - return the string representation of the converted form.
1783
    # There is no simple str() method defined for p's class.
1784
    return(p.__str__())
1785
#
1786
print pobyso_get_prec()  
1787
pobyso_set_prec(165)
1788
print pobyso_get_prec()  
1789
a=100
1790
print type(a)
1791
id(a)
1792
print "Max arity: ", pobyso_max_arity
1793
print "Function tripleDouble (43) as a string: ", pobyso_function_type_as_string(43)
1794
print "Function None (44) as a string: ", pobyso_function_type_as_string(44)
1795
print "...Pobyso check done"