Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 221

Historique | Voir | Annoter | Télécharger (68,56 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
                sollya_lib_clear_obj(resPolySo)
1235
                sollya_lib_clear_obj(infNormSo)
1236
                return (lastResPolySo, lastInfNormSo)
1237
        # cerrSa <= maxErrSa: scrap more bits, possibly.
1238
        else:
1239
            if debug:
1240
                print "Error is too small"
1241
            if cerrSa <= (maxErrSa/2):
1242
                if debug: 
1243
                    print "Shrinking prec."
1244
                ntpSa = floor(ftpSa - ftpSa/50)
1245
                ## Can't shrink (numerical)
1246
                if ntpSa == ftpSa:
1247
                    if not lastResPolySo is None:
1248
                        sollya_lib_clear_obj(lastResPolySo)
1249
                    if not lastInfNormSo is None:
1250
                        sollya_lib_clear_obj(lastInfNormSo)
1251
                    return (resPolySo, infNormSo)
1252
                ## Can't shrink (not enough precision left)
1253
                if ntpSa <= itpSa:
1254
                    if not lastResPolySo is None:
1255
                        sollya_lib_clear_obj(lastResPolySo)
1256
                    if not lastInfNormSo is None:
1257
                        sollya_lib_clear_obj(lastInfNormSo)
1258
                    return (resPolySo, infNormSo)
1259
                ftpSa = ntpSa
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
                lastResPolySo = resPolySo
1265
                lastInfNormSo = infNormSo
1266
                continue
1267
            else: # Error is not that small, just return
1268
                if not lastResPolySo is None:
1269
                    sollya_lib_clear_obj(lastResPolySo)
1270
                if not lastInfNormSo is None:
1271
                    sollya_lib_clear_obj(lastInfNormSo)
1272
                return (resPolySo, infNormSo)                
1273
    # End wile True
1274
    return None
1275
# End pobyso_polynomial_coefficients_progressive_truncate_so_so.
1276

    
1277
def pobyso_polynomial_degree_so_sa(polySo):
1278
    """
1279
    Return the degree of a Sollya polynomial as a Sage int.
1280
    """
1281
    degreeSo = sollya_lib_degree(polySo)
1282
    return pobyso_constant_from_int_so_sa(degreeSo)
1283
# End pobyso_polynomial_degree_so_sa
1284

    
1285
def pobyso_polynomial_degree_so_so(polySo):
1286
    """
1287
    Thin wrapper around lib_sollya_degree().
1288
    """
1289
    return sollya_lib_degree(polySo)
1290
# End pobyso_polynomial_degree_so_so
1291
                                                                  
1292
def pobyso_range(rnLowerBound, rnUpperBound):
1293
    """ Legacy function. See pobyso_range_sa_so. """
1294
    return pobyso_range_sa_so(rnLowerBound, rnUpperBound) 
1295

    
1296

    
1297
def pobyso_range_to_interval_so_sa(rangeSo, realIntervalFieldSa = None):
1298
    """
1299
    Get a Sage interval from a Sollya range.
1300
    If no realIntervalField is given as a parameter, the Sage interval
1301
    precision is that of the Sollya range.
1302
    Otherwise, the precision is that of the realIntervalField. In this case
1303
    rounding may happen.
1304
    """
1305
    if realIntervalFieldSa is None:
1306
        precSa = pobyso_get_prec_of_range_so_sa(rangeSo)
1307
        realIntervalFieldSa = RealIntervalField(precSa)
1308
    intervalSa = \
1309
        pobyso_get_interval_from_range_so_sa(rangeSo, realIntervalFieldSa) 
1310
    return intervalSa
1311
# End pobyso_range_to_interval_so_sa
1312

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

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

    
1455
def pobyso_remez_canonical_so_so(funcSo, \
1456
                                 degreeSo, \
1457
                                 rangeSo, \
1458
                                 weightSo = pobyso_constant_1_sa_so(),\
1459
                                 qualitySo = None):
1460
    """
1461
    All arguments are pointers to Sollya objects.
1462
    The return value is a pointer to a Sollya function.
1463
    """
1464
    if not sollya_lib_obj_is_function(funcSo):
1465
        return(None)
1466
    return(sollya_lib_remez(funcSo, degreeSo, rangeSo, weightSo, qualitySo, None))
1467
# End pobyso_remez_canonical_so_so.
1468

    
1469
def pobyso_round_coefficients_single_so_so(polySo, precSo):
1470
    """
1471
    Create a rounded coefficients polynomial from polynomial argument to
1472
    the number of bits in size argument.
1473
    All coefficients are set to the same precision.
1474
    """
1475
    ## TODO: check arguments.
1476
    endEllipListSo = pobyso_build_end_elliptic_list_so_so(precSo)
1477
    polySo = sollya_lib_roundcoefficients(polySo, endEllipListSo, None)
1478
    sollya_lib_clear_obj(endEllipListSo)
1479
    #sollya_lib_clear_obj(endEllipListSo)
1480
    return polySo
1481
    
1482
# End pobyso_round_coefficients_single_so_so
1483

    
1484
def pobyso_set_canonical_off():
1485
    sollya_lib_set_canonical(sollya_lib_off())
1486

    
1487
def pobyso_set_canonical_on():
1488
    sollya_lib_set_canonical(sollya_lib_on())
1489

    
1490
def pobyso_set_prec(p):
1491
    """ Legacy function. See pobyso_set_prec_sa_so. """
1492
    pobyso_set_prec_sa_so(p)
1493

    
1494
def pobyso_set_prec_sa_so(p):
1495
    a = c_int(p)
1496
    precSo = c_void_p(sollya_lib_constant_from_int(a))
1497
    sollya_lib_set_prec(precSo, None)
1498
# End pobyso_set_prec_sa_so.
1499

    
1500
def pobyso_set_prec_so_so(newPrecSo):
1501
    sollya_lib_set_prec(newPrecSo, None)
1502
# End pobyso_set_prec_so_so.
1503

    
1504
def pobyso_inf_so_so(intervalSo):
1505
    """
1506
    Very thin wrapper around sollya_lib_inf().
1507
    """
1508
    return sollya_lib_inf(intervalSo)
1509
# End pobyso_inf_so_so.
1510
    
1511
def pobyso_supnorm_so_so(polySo, funcSo, intervalSo, errorTypeSo = None,\
1512
                         accuracySo = None):
1513
    """
1514
    Computes the supnorm of the approximation error between the given 
1515
    polynomial and function.
1516
    errorTypeSo defaults to "absolute".
1517
    accuracySo defaults to 2^(-40).
1518
    """
1519
    if errorTypeSo is None:
1520
        errorTypeSo = sollya_lib_absolute(None)
1521
        errorTypeIsNone = True
1522
    else:
1523
        errorTypeIsNone = False
1524
    #
1525
    if accuracySo is None:
1526
        # Notice the **!
1527
        accuracySo = pobyso_constant_sa_so(RR(2**(-40)))
1528
        accuracyIsNone = True
1529
    else:
1530
        accuracyIsNone = False
1531
    pobyso_autoprint(accuracySo)
1532
    resultSo = \
1533
        sollya_lib_supnorm(polySo, funcSo, intervalSo, errorTypeSo, \
1534
                              accuracySo)
1535
    if errorTypeIsNone:
1536
        sollya_lib_clear_obj(errorTypeSo)
1537
    if accuracyIsNone:
1538
        sollya_lib_clear_obj(accuracySo)
1539
    return resultSo
1540
# End pobyso_supnorm_so_so
1541

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

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

    
1660
def pobyso_taylor(function, degree, point):
1661
    """ Legacy function. See pobysoTaylor_so_so. """
1662
    return(pobyso_taylor_so_so(function, degree, point))
1663

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

    
1753
def pobyso_taylorform_so_so(functionSo, degreeSo, pointSo, intervalSo=None, \
1754
                            errorTypeSo=None):
1755
    createdErrorType = False
1756
    if errorTypeSo is None:
1757
        errorTypeSo = sollya_lib_absolute()
1758
        createdErrorType = True
1759
    else:
1760
        #TODO: deal with the other case.
1761
        pass
1762
    if intervalSo is None:
1763
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
1764
                                         errorTypeSo, None)
1765
    else:
1766
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
1767
                                         intervalSo, errorTypeSo, None)
1768
    if createdErrorType:
1769
        sollya_lib_clear_obj(errorTypeSo)
1770
    return resultSo
1771
        
1772

    
1773
def pobyso_univar_polynomial_print_reverse(polySa):
1774
    """ Legacy function. See pobyso_univar_polynomial_print_reverse_sa_sa. """
1775
    return(pobyso_univar_polynomial_print_reverse_sa_sa(polySa))
1776

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