Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 261

Historique | Voir | Annoter | Télécharger (87,63 ko)

1
"""
2
@file pobyso.py
3
Actual functions to use in Sage
4
@author S.T. 
5
@date 2012-11-13
6

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

10
pobyso functions come in five flavors: 
11
- the _so_so (arguments and returned objects are pointers to Sollya objects, 
12
  includes the void function and the no arguments function that return a 
13
  pointer to a Sollya object);
14
- the _so_sa (argument are pointers to Sollya objects, returned objects are
15
  Sage/Python objects or, more generally, information is transfered from the
16
  Sollya world to Sage/Python world; e.g. functions without arguments that
17
  return a Sage/Python object);
18
- the _sa_so (arguments are Sage/Python objects, returned objects are 
19
  pointers to Sollya objects);
20
- the sa_sa (arguments and returned objects are all Sage/Python objects);
21
- a catch all flavor, without any suffix, (e. g. functions that have no 
22
  argument nor return value).
23

24
This classification is not always very strict. Conversion functions from Sollya
25
to Sage/Python are sometimes decorated with Sage/Python arguments to set
26
the precision. These functions remain in the so_sa category.
27

28
@note
29
Reported errors in Eclipse come from the calls to the Sollya library
30

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

    
89
pobyso_max_arity = 9
90

    
91
def pobyso_absolute_so_so():
92
    """
93
    Create an "absolute" Sollya object.
94
    """
95
    return(sollya_lib_absolute(None))
96

    
97
def pobyso_autoprint(arg):
98
    sollya_lib_autoprint(arg, None)
99

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

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

    
160
def pobyso_build_function_sub_so_so(exp1So, exp2So):
161
    return sollya_lib_build_function_sub(exp1So, exp2So)
162

    
163
def pobyso_build_list_of_ints_sa_so(*args):
164
    """
165
    Build a Sollya list from Sage integral arguments. 
166
    """
167
    if len(args) == 0:
168
        return pobyso_build_list_so_so()
169
    ## Make a Sage list of integral constants.
170
    intsList = []
171
    for intElem in args:
172
        intsList.append(pobyso_constant_from_int_sa_so(intElem))
173
    return pobyso_build_list_so_so(*intsList) 
174
    
175
def pobyso_build_list_so_so(*args):
176
    """
177
    Make a Sollya list out of Sollya objects passed as arguments.
178
    If one wants to call it with a list argument, should prepend a "*"
179
    before the list variable name.
180
    Elements of the list are "eaten" (should not be cleared individually, 
181
    are cleared when the list is cleared).
182
    """
183
    if len(args) == 0:
184
        ## Called with an empty list produced "error".
185
        return sollya_lib_build_list(None)
186
    index = 0
187
    ## Append the Sollya elements one by one.
188
    for elementSo in args:
189
        if index == 0:
190
            listSo = sollya_lib_build_list(elementSo, None)
191
        else:
192
            listSo = sollya_lib_append(listSo, elementSo)
193
        index += 1
194
    return listSo
195
# End pobyso_build list_so_so
196
    
197

    
198
def pobyso_change_var_in_function_so_so(funcSo, chvarExpSo):
199
    """
200
    Variable change in a function.
201
    """
202
    return(sollya_lib_evaluate(funcSo,chvarExpSo))
203
# End pobyso_change_var_in_function_so_so     
204

    
205
def pobyso_chebyshevform_so_so(functionSo, degreeSo, intervalSo):
206
    resultSo = sollya_lib_chebyshevform(functionSo, degreeSo, intervalSo)
207
    return(resultSo)
208
# End pobyso_chebyshevform_so_so.
209

    
210
def pobyso_clear_obj(objSo):
211
    """
212
    Free a Sollya object's memory.
213
    Very thin wrapper around sollya_lib_clear_obj().
214
    """
215
    sollya_lib_clear_obj(objSo)
216
# End pobyso_clear_obj. 
217

    
218
def pobyso_clear_taylorform_sa_so(taylorFormSaSo):
219
    """
220
    This method is necessary to correctly clean up the memory from Taylor forms.
221
    These are made of a Sollya object, a Sollya object list, a Sollya object.
222
    For no clearly understood reason, sollya_lib_clear_object_list crashed 
223
    when applied to the object list.
224
    Here, we decompose it into Sage list of Sollya objects references and we
225
     clear them one at a time. 
226
    """
227
    sollya_lib_clear_obj(taylorFormSaSo[0])
228
    (coefficientsErrorsListSaSo, numElementsSa, isEndEllipticSa) = \
229
        pobyso_get_list_elements_so_so(taylorFormSaSo[1])
230
    for element in coefficientsErrorsListSaSo:
231
        sollya_lib_clear_obj(element)
232
    sollya_lib_clear_obj(taylorFormSaSo[1])
233
    sollya_lib_clear_obj(taylorFormSaSo[2])
234
# End pobyso_clear_taylorform_sa_so 
235

    
236
def pobyso_cmp(rnArgSa, cteSo):
237
    """
238
    Deprecated, use pobyso_cmp_sa_so_sa instead.
239
    """
240
    print "Deprecated, use pobyso_cmp_sa_so_sa instead."
241
    return pobyso_cmp_sa_so_sa(rnArgSa, cteSo)
242
# End  pobyso_cmp
243
    
244
def pobyso_cmp_sa_so_sa(rnArgSa, cteSo):
245
    """
246
    Compare the MPFR value a RealNumber with that of a Sollya constant.
247
    
248
    Get the value of the Sollya constant into a RealNumber and compare
249
    using MPFR. Could be optimized by working directly with a mpfr_t
250
    for the intermediate number. 
251
    """
252
    # Get the precision of the Sollya constant to build a Sage RealNumber 
253
    # with enough precision.to hold it.
254
    precisionOfCte = c_int(0)
255
    # From the Sollya constant, create a local Sage RealNumber.
256
    sollya_lib_get_prec_of_constant(precisionOfCte, cteSo) 
257
    #print "Precision of constant: ", precisionOfCte
258
    RRRR = RealField(precisionOfCte.value)
259
    rnLocalSa = RRRR(0)
260
    sollya_lib_get_constant(get_rn_value(rnLocalSa), cteSo)
261
    #
262
    ## Compare the Sage RealNumber version of the Sollya constant with rnArg
263
    #  through a direct comparison of underlying MPFR numbers.
264
    return cmp_rn_value(rnArgSa, rnLocal)
265
# End pobyso_smp_sa_so_sa
266

    
267
def pobyso_compute_pos_function_abs_val_bounds_sa_sa(funcSa, lowerBoundSa, \
268
                                                     upperBoundSa):
269
    """
270
    TODO: completely rework and test.
271
    """
272
    pobyso = pobyso_name_free_variable_sa_so(funcSa.variables()[0])
273
    funcSo = pobyso_parse_string(funcSa._assume_str().replace('_SAGE_VAR_', ''))
274
    rangeSo = pobyso_range_sa_so(lowerBoundSa, upperBoundSa)
275
    infnormSo = pobyso_infnorm_so_so(funcSo,rangeSo)
276
    # Sollya return the infnorm as an interval.
277
    fMaxSa = pobyso_get_interval_from_range_so_sa(infnormSo)
278
    # Get the top bound and compute the binade top limit.
279
    fMaxUpperBoundSa = fMaxSa.upper()
280
    binadeTopLimitSa = 2**ceil(fMaxUpperBoundSa.log2())
281
    # Put up together the function to use to compute the lower bound.
282
    funcAuxSo = pobyso_parse_string(str(binadeTopLimitSa) +  \
283
                                    '-(' + f._assume_str().replace('_SAGE_VAR_', '') + ')')
284
    pobyso_autoprint(funcAuxSo)
285
    # Clear the Sollya range before a new call to infnorm and issue the call.
286
    sollya_lib_clear_obj(infnormSo)
287
    infnormSo = pobyso_infnorm_so_so(funcAuxSo,rangeSo)
288
    fMinSa = pobyso_get_interval_from_range_so_sa(infnormSo)
289
    sollya_lib_clear_obj(infnormSo)
290
    fMinLowerBoundSa = binadeTopLimitSa - fMinSa.lower()
291
    # Compute the maximum of the precisions of the different bounds.
292
    maxPrecSa = max([fMinLowerBoundSa.parent().precision(), \
293
                     fMaxUpperBoundSa.parent().precision()])
294
    # Create a RealIntervalField and create an interval with the "good" bounds.
295
    RRRI = RealIntervalField(maxPrecSa)
296
    imageIntervalSa = RRRI(fMinLowerBoundSa, fMaxUpperBoundSa)
297
    # Free the unneeded Sollya objects
298
    sollya_lib_clear_obj(funcSo)
299
    sollya_lib_clear_obj(funcAuxSo)
300
    sollya_lib_clear_obj(rangeSo)
301
    return(imageIntervalSa)
302
# End pobyso_compute_pos_function_abs_val_bounds_sa_sa
303

    
304
def pobyso_compute_precision_decay_ratio_function_sa_so():
305
    """
306
    Compute the precision decay ratio function for polynomial 
307
    coefficient progressive trucation.
308
    """
309
    functionText = """
310
    proc(deg, a, b, we, wq) 
311
    {
312
      k = we * (exp(x/a)-1) + wq * (b*x)^2 + (1-we-wq) * x;
313
      return k/k(d);
314
    };
315
    """
316
    return pobyso_parse_string_sa_so(functionText)
317
# End  pobyso_compute_precision_decay_ratio_function.
318

    
319

    
320
def pobyso_constant(rnArg):
321
    """ Legacy function. See pobyso_constant_sa_so. """
322
    return(pobyso_constant_sa_so(rnArg))
323
    
324
def pobyso_constant_sa_so(rnArgSa, precisionSa=None):
325
    """
326
    Create a Sollya constant from a Sage RealNumber.
327
    The sollya_lib_constant() function creates a constant
328
    with the same precision as the source.
329
    """
330
    ## Precision stuff. If one wants to change precisions,
331
    #  everything takes place in Sage. That only makes
332
    #  sense if one wants to reduce the precision.
333
    # TODO: revisit precision stuff with new technique.
334
    if not precisionSa is None:
335
        RRR = RealField(precisionSa)
336
        rnArgSa = RRR(rnArgSa)
337
    #print rnArgSa, rnArgSa.precision()
338
    # Sollya constant creation takes place here.
339
    return sollya_lib_constant(get_rn_value(rnArgSa))
340
# End pobyso_constant_sa_so
341
     
342
def pobyso_constant_0_sa_so():
343
    """
344
    Obvious.
345
    """
346
    return pobyso_constant_from_int_sa_so(0)
347

    
348
def pobyso_constant_1():
349
    """
350
    Obvious.
351
    Legacy function. See pobyso_constant_so_so. 
352
    """
353
    return pobyso_constant_1_sa_so()
354

    
355
def pobyso_constant_1_sa_so():
356
    """
357
    Obvious.
358
    """
359
    return(pobyso_constant_from_int_sa_so(1))
360

    
361
def pobyso_constant_from_int(anInt):
362
    """ Legacy function. See pobyso_constant_from_int_sa_so. """
363
    return pobyso_constant_from_int_sa_so(anInt)
364

    
365
def pobyso_constant_from_int_sa_so(anInt):
366
    """
367
    Get a Sollya constant from a Sage int.
368
    """
369
    return sollya_lib_constant_from_int64(long(anInt))
370

    
371
def pobyso_constant_from_int_so_sa(constSo):
372
    """
373
    Get a Sage int from a Sollya int constant.
374
    Usefull for precision or powers in polynomials.
375
    """
376
    constSa = c_long(0)
377
    sollya_lib_get_constant_as_int64(byref(constSa), constSo)
378
    return constSa.value
379
# End pobyso_constant_from_int_so_sa
380

    
381
def pobyso_constant_from_mpq_sa_so(rationalSa):
382
    """
383
    Make a Sollya constant from Sage rational.
384
    The Sollya constant is an unevaluated expression.
385
    Hence no precision argument is needed.
386
    It is better to leave this way since Sollya has its own
387
    optimized evaluation mecanism that tries very hard to
388
    return exact values or at least faithful ones.
389
    """
390
    ratExprSo = \
391
        sollya_lib_constant_from_mpq(sgmp_get_rational_value(rationalSa))
392
    return ratExprSo
393
# End pobyso_constant_from_mpq_sa_so.
394

    
395
def pobyso_constant_sollya_prec_sa_so(rnArgSa):
396
    """
397
    Create a Sollya constant from a Sage RealNumber at the
398
    current precision in Sollya.
399
    """
400
    currentSollyaPrecSa = pobyso_get_prec_so_sa()
401
    return pobyso_constant_sa_so(rnArgSa, currentSollyaPrecSa)
402
# End pobyso_constant_sollya_prec_sa_so
403

    
404
def pobyso_end_elliptic_list_so_sa_so(objectsListSo, intCountSa):
405
    """
406
    Create a Sollya end elliptic list made of the objectListSo[0] to
407
     objectsListSo[intCountSa-1] objects.
408
    """     
409
    return sollya_lib_end_elliptic_list(objectSo, int(intCountSa))
410

    
411
def pobyso_error_so():
412
    return sollya_lib_error(None)
413
# End pobyso_error().
414

    
415
def pobyso_evaluate_so_so(funcSo, argumentSo):
416
    """
417
    Evaluates funcSo for arguemntSo through sollya_lib_evaluate().
418
    """
419
    return sollya_lib_evaluate(funcSo, argumentSo)
420
# End pobyso_evaluate_so_so.
421
#
422
def pobyso_diff_so_so(funcSo):
423
    """
424
    Very thin wrapper around sollya_lib_diff.
425
    """
426
    ## TODO: add a check to make sure funcSo is a functional expression.
427
    return sollya_lib_diff(funcSo)
428

    
429
def pobyso_dirty_find_zeros_so_so(funcSo, rangeSo):
430
    """
431
    Thin wrapper over sollya_lib_dirtyfindzeros()
432
    """
433
    return sollya_lib_dirtyfindzeros(funcSo, rangeSo)
434
# End pobys_dirty_find_zeros
435

    
436
def pobyso_dirty_inf_norm_so_so(funcSo, rangeSo, preSo=None):
437
    """
438
    Thin wrapper around sollya_dirtyinfnorm().
439
    """
440
    # TODO: manage the precision change.
441
    
442
    return sollya_lib_dirtyinfnorm(funcSo, rangeSo)
443
# End pobyso_dirty_inf_norm_so_so
444

    
445
def pobyso_find_zeros_so_so(funcSo, rangeSo):
446
    """
447
    Thin wrapper over sollya_lib_findzeros()
448
    """
449
    return sollya_lib_findzeros(funcSo, rangeSo)
450
# End pobys_find_zeros
451

    
452
def pobyso_float_list_so_sa(listSo):
453
    """
454
    Return a Sollya list of floating-point numbers as a Sage list of 
455
    floating-point numbers.
456
    TODO: add a test to make sure that each element of the list is a constant.
457
    """
458
    listSa   = []
459
    ## The function returns none if the list is empty or an error has happened.
460
    retVal = pobyso_get_list_elements_so_so(listSo)
461
    if retVal is None:
462
        return listSa
463
    ## Just in case the interface is changed and an empty list is returned 
464
    #  instead of None.
465
    elif len(retVal) == 0:
466
        return listSa
467
    else:
468
        ## Remember pobyso_get_list_elements_so_so returns more information
469
        #  than just the elements of the list (# elements, is_elliptic)
470
        listSaSo, numElements, isEndElliptic = retVal
471
    ## Return an empty list.
472
    if numElements == 0:
473
        return listSa
474
    ## Search first for the maximum precision of the elements
475
    maxPrecSa = 0
476
    for floatSo in listSaSo:
477
        #pobyso_autoprint(floatSo)
478
        curPrecSa =  pobyso_get_prec_of_constant_so_sa(floatSo)
479
        if curPrecSa > maxPrecSa:
480
            maxPrecSa = curPrecSa
481
    ##
482
    RF = RealField(maxPrecSa)
483
    ##
484
    for floatSo in listSaSo:
485
        listSa.append(pobyso_get_constant_as_rn_with_rf_so_sa(floatSo))
486
    return listSa
487
# End pobyso_float_list_so_sa
488

    
489
def pobyso_float_poly_sa_so(polySa, precSa = None):
490
    """
491
    Create a Sollya polynomial from a Sage RealField polynomial.
492
    """
493
    ## TODO: filter arguments.
494
    ## Precision. If a precision is given, convert the polynomial
495
    #  into the right polynomial field. If not convert it straight
496
    #  to Sollya.
497
    sollyaPrecChanged = False 
498
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
499
    if precSa is None:
500
        precSa = polySa.parent().base_ring().precision()
501
    if (precSa > initialSollyaPrecSa):
502
        assert precSa >= 2, "Precision change <2 requested"
503
        if precSa <= 2:
504
            print inspect.stack()[0][3], ": precision change <= 2 requested"
505
        precSo = pobyso_constant_from_int(precSa)
506
        pobyso_set_prec_so_so(precSo)
507
        sollya_lib_clear_obj(precSo)
508
        sollyaPrecChanged = True    
509
    ## Get exponents and coefficients.
510
    exponentsSa     = polySa.exponents()
511
    coefficientsSa  = polySa.coefficients()
512
    ## Build the polynomial.
513
    polySo = None
514
    for coefficientSa, exponentSa in zip(coefficientsSa, exponentsSa):
515
        #print coefficientSa.n(prec=precSa), exponentSa
516
        coefficientSo = \
517
            pobyso_constant_sa_so(coefficientSa)
518
        #pobyso_autoprint(coefficientSo)
519
        exponentSo = \
520
            pobyso_constant_from_int_sa_so(exponentSa)
521
        #pobyso_autoprint(exponentSo)
522
        monomialSo = sollya_lib_build_function_pow(
523
                       sollya_lib_build_function_free_variable(),
524
                       exponentSo)
525
        polyTermSo = sollya_lib_build_function_mul(coefficientSo,
526
                                                       monomialSo)
527
        if polySo is None:
528
            polySo = polyTermSo
529
        else:
530
            polySo = sollya_lib_build_function_add(polySo, polyTermSo)
531
    if sollyaPrecChanged:
532
        pobyso_set_prec_so_so(initialSollyaPrecSo)
533
        sollya_lib_clear_obj(initialSollyaPrecSo)
534
    return polySo
535
# End pobyso_float_poly_sa_so    
536

    
537
def pobyso_float_poly_so_sa(polySo, realFieldSa=None):
538
    """
539
    Convert a Sollya polynomial into a Sage floating-point polynomial.
540
    If no realField is given, a RealField corresponding to the maximum 
541
    precision of the coefficients is internally computed.
542
    The real field is not returned but can be easily retrieved from 
543
    the polynomial itself.
544
    ALGORITHM:
545
    - (optional) compute the RealField of the coefficients;
546
    - convert the Sollya expression into a Sage expression;
547
    - convert the Sage expression into a Sage polynomial
548
    """    
549
    if realFieldSa is None:
550
        expressionPrecSa = pobyso_get_max_prec_of_exp_so_sa(polySo)
551
        #print "Maximum precision of Sollya polynomial coefficients:", expressionPrecSa
552
        if expressionPrecSa < 2 or expressionPrecSa > 2147483391:
553
            print "Maximum degree of expression:", expressionPrecSa
554
        realFieldSa      = RealField(expressionPrecSa)
555
    #print "Sollya expression before...",
556
    #pobyso_autoprint(polySo)
557

    
558
    expressionSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo,
559
                                                             realFieldSa)
560
    #print "...Sollya expression after."
561
    #pobyso_autoprint(polySo)
562
    polyVariableSa = expressionSa.variables()[0]
563
    polyRingSa     = realFieldSa[str(polyVariableSa)]
564
    #print polyRingSa
565
    # Do not use the polynomial(expressionSa, ring=polyRingSa) form!
566
    polynomialSa = polyRingSa(expressionSa)
567
    polyCoeffsListSa = polynomialSa.coefficients()
568
    #for coeff in polyCoeffsListSa:
569
    #    print coeff.abs().n()
570
    return polynomialSa
571
# End pobyso_float_poly_so_sa
572

    
573
def pobyso_free_variable():
574
    """
575
    Ultra thin wrapper around the sollya_lib_function_build_free_variable function.
576
    """
577
    return sollya_lib_build_function_free_variable()
578
    
579
def pobyso_function_type_as_string(funcType):
580
    """ Legacy function. See pobyso_function_type_as_string_so_sa. """
581
    return(pobyso_function_type_as_string_so_sa(funcType))
582

    
583
def pobyso_function_type_as_string_so_sa(funcType):
584
    """
585
    Numeric Sollya function codes -> Sage mathematical function names.
586
    Notice that pow -> ^ (a la Sage, not a la Python).
587
    """
588
    if funcType == SOLLYA_BASE_FUNC_ABS:
589
        return "abs"
590
    elif funcType == SOLLYA_BASE_FUNC_ACOS:
591
        return "arccos"
592
    elif funcType == SOLLYA_BASE_FUNC_ACOSH:
593
        return "arccosh"
594
    elif funcType == SOLLYA_BASE_FUNC_ADD:
595
        return "+"
596
    elif funcType == SOLLYA_BASE_FUNC_ASIN:
597
        return "arcsin"
598
    elif funcType == SOLLYA_BASE_FUNC_ASINH:
599
        return "arcsinh"
600
    elif funcType == SOLLYA_BASE_FUNC_ATAN:
601
        return "arctan"
602
    elif funcType == SOLLYA_BASE_FUNC_ATANH:
603
        return "arctanh"
604
    elif funcType == SOLLYA_BASE_FUNC_CEIL:
605
        return "ceil"
606
    elif funcType == SOLLYA_BASE_FUNC_CONSTANT:
607
        return "cte"
608
    elif funcType == SOLLYA_BASE_FUNC_COS:
609
        return "cos"
610
    elif funcType == SOLLYA_BASE_FUNC_COSH:
611
        return "cosh"
612
    elif funcType == SOLLYA_BASE_FUNC_DIV:
613
        return "/"
614
    elif funcType == SOLLYA_BASE_FUNC_DOUBLE:
615
        return "double"
616
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEDOUBLE:
617
        return "doubleDouble"
618
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEEXTENDED:
619
        return "doubleDxtended"
620
    elif funcType == SOLLYA_BASE_FUNC_ERF:
621
        return "erf"
622
    elif funcType == SOLLYA_BASE_FUNC_ERFC:
623
        return "erfc"
624
    elif funcType == SOLLYA_BASE_FUNC_EXP:
625
        return "exp"
626
    elif funcType == SOLLYA_BASE_FUNC_EXP_M1:
627
        return "expm1"
628
    elif funcType == SOLLYA_BASE_FUNC_FLOOR:
629
        return "floor"
630
    elif funcType == SOLLYA_BASE_FUNC_FREE_VARIABLE:
631
        return "freeVariable"
632
    elif funcType == SOLLYA_BASE_FUNC_HALFPRECISION:
633
        return "halfPrecision"
634
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYCONSTANT:
635
        return "libraryConstant"
636
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYFUNCTION:
637
        return "libraryFunction"
638
    elif funcType == SOLLYA_BASE_FUNC_LOG:
639
        return "log"
640
    elif funcType == SOLLYA_BASE_FUNC_LOG_10:
641
        return "log10"
642
    elif funcType == SOLLYA_BASE_FUNC_LOG_1P:
643
        return "log1p"
644
    elif funcType == SOLLYA_BASE_FUNC_LOG_2:
645
        return "log2"
646
    elif funcType == SOLLYA_BASE_FUNC_MUL:
647
        return "*"
648
    elif funcType == SOLLYA_BASE_FUNC_NEARESTINT:
649
        return "round"
650
    elif funcType == SOLLYA_BASE_FUNC_NEG:
651
        return "__neg__"
652
    elif funcType == SOLLYA_BASE_FUNC_PI:
653
        return "pi"
654
    elif funcType == SOLLYA_BASE_FUNC_POW:
655
        return "^"
656
    elif funcType == SOLLYA_BASE_FUNC_PROCEDUREFUNCTION:
657
        return "procedureFunction"
658
    elif funcType == SOLLYA_BASE_FUNC_QUAD:
659
        return "quad"
660
    elif funcType == SOLLYA_BASE_FUNC_SIN:
661
        return "sin"
662
    elif funcType == SOLLYA_BASE_FUNC_SINGLE:
663
        return "single"
664
    elif funcType == SOLLYA_BASE_FUNC_SINH:
665
        return "sinh"
666
    elif funcType == SOLLYA_BASE_FUNC_SQRT:
667
        return "sqrt"
668
    elif funcType == SOLLYA_BASE_FUNC_SUB:
669
        return "-"
670
    elif funcType == SOLLYA_BASE_FUNC_TAN:
671
        return "tan"
672
    elif funcType == SOLLYA_BASE_FUNC_TANH:
673
        return "tanh"
674
    elif funcType == SOLLYA_BASE_FUNC_TRIPLEDOUBLE:
675
        return "tripleDouble"
676
    else:
677
        return None
678

    
679
def pobyso_get_constant(rnArgSa, constSo):
680
    """ Legacy function. See pobyso_get_constant_so_sa. """
681
    return pobyso_get_constant_so_sa(rnArgSa, constSo)
682
# End pobyso_get_constant
683

    
684
def pobyso_get_constant_so_sa(rnArgSa, constSo):
685
    """
686
    Set the value of rnArgSo to the value of constSo in MPFR_RNDN mode.
687
    rnArg must already exist and belong to some RealField.
688
    We assume that constSo points to a Sollya constant.
689
    """
690
    outcome = sollya_lib_get_constant(get_rn_value(rnArgSa), constSo)
691
    if outcome == 0: # Failure because constSo is not a constant expression.
692
        return None
693
    else:
694
        return outcome
695
# End  pobyso_get_constant_so_sa
696
   
697
def pobyso_get_constant_as_rn(ctExpSo):
698
    """ 
699
    Legacy function. See pobyso_get_constant_as_rn_so_sa. 
700
    """ 
701
    return(pobyso_get_constant_as_rn_so_sa(ctExpSo))
702
    
703
def pobyso_get_constant_as_rn_so_sa(constExpSo):
704
    """
705
    Get a Sollya constant as a Sage "real number".
706
    The precision of the floating-point number returned is that of the Sollya
707
    constant.
708
    """
709
    #print "Before computing precision of variable..."
710
    #pobyso_autoprint(constExpSo)
711
    precisionSa  = pobyso_get_prec_of_constant_so_sa(constExpSo)
712
    #print "precisionSa:", precisionSa
713
    ## If the expression can not be exactly converted, None is returned.
714
    #  In this case opt for the Sollya current expression.
715
    if precisionSa is None:
716
        precisionSa = pobyso_get_prec_so_sa()
717
    RRRR = RealField(precisionSa)
718
    rnSa = RRRR(0)
719
    outcome = sollya_lib_get_constant(get_rn_value(rnSa), constExpSo)
720
    if outcome == 0:
721
        return None
722
    else:
723
        return rnSa
724
# End pobyso_get_constant_as_rn_so_sa
725

    
726
def pobyso_get_constant_as_rn_with_rf(ctExp, realField):
727
    """ 
728
    Legacy function. See pobyso_get_constant_as_rn_with_rf_so_sa.
729
    """
730
    return pobyso_get_constant_as_rn_with_rf_so_sa(ctExp, realField)
731
# End pobyso_get_constant_as_rn_with_rf
732
    
733
def pobyso_get_constant_as_rn_with_rf_so_sa(ctExpSo, realFieldSa = None):
734
    """
735
    Get a Sollya constant as a Sage "real number".
736
    If no real field is specified, the precision of the floating-point number 
737
    returned is that of the Sollya constant.
738
    Otherwise is is that of the real field. Hence rounding may happen.
739
    """
740
    if realFieldSa is None:
741
        return pobyso_get_constant_as_rn_so_sa(ctExpSo)
742
    rnSa = realFieldSa(0)
743
    outcome = sollya_lib_get_constant(get_rn_value(rnSa), ctExpSo)
744
    if outcome == 0:
745
        return None
746
    else:
747
        return rnSa
748
# End pobyso_get_constant_as_rn_with_rf_so_sa
749

    
750
def pobyso_get_free_variable_name():
751
    """ 
752
    Legacy function. See pobyso_get_free_variable_name_so_sa.
753
    """
754
    return(pobyso_get_free_variable_name_so_sa())
755

    
756
def pobyso_get_free_variable_name_so_sa():
757
    return sollya_lib_get_free_variable_name()
758
    
759
def pobyso_get_function_arity(expressionSo):
760
    """ 
761
    Legacy function. See pobyso_get_function_arity_so_sa.
762
    """
763
    return(pobyso_get_function_arity_so_sa(expressionSo))
764

    
765
def pobyso_get_function_arity_so_sa(expressionSo):
766
    arity = c_int(0)
767
    sollya_lib_get_function_arity(byref(arity),expressionSo)
768
    return int(arity.value)
769

    
770
def pobyso_get_head_function(expressionSo):
771
    """ 
772
    Legacy function. See pobyso_get_head_function_so_sa. 
773
    """
774
    return(pobyso_get_head_function_so_sa(expressionSo)) 
775

    
776
def pobyso_get_head_function_so_sa(expressionSo):
777
    functionType = c_int(0)
778
    sollya_lib_get_head_function(byref(functionType), expressionSo)
779
    return int(functionType.value)
780

    
781
def pobyso_get_interval_from_range_so_sa(soRange, realIntervalFieldSa = None ):
782
    """
783
    Return the Sage interval corresponding to the Sollya range argument.
784
    If no reaIntervalField is passed as an argument, the interval bounds are not
785
    rounded: they are elements of RealIntervalField of the "right" precision
786
    to hold all the digits.
787
    """
788
    prec = c_int(0)
789
    if realIntervalFieldSa is None:
790
        retval = sollya_lib_get_prec_of_range(byref(prec), soRange, None)
791
        if retval == 0:
792
            return None
793
        realIntervalFieldSa = RealIntervalField(prec.value)
794
    intervalSa = realIntervalFieldSa(0,0)
795
    retval = \
796
        sollya_lib_get_interval_from_range(get_interval_value(intervalSa),\
797
                                           soRange)
798
    if retval == 0:
799
        return None
800
    return intervalSa
801
# End pobyso_get_interval_from_range_so_sa
802

    
803
def pobyso_get_list_elements(soObj):
804
    """ Legacy function. See pobyso_get_list_elements_so_so. """
805
    return pobyso_get_list_elements_so_so(soObj)
806
 
807
def pobyso_get_list_elements_so_so(objectListSo):
808
    """
809
    Get the Sollya list elements as a Sage/Python array of Sollya objects.
810
    
811
    INPUT:
812
    - objectListSo: a Sollya list of Sollya objects.
813
    
814
    OUTPUT:
815
    - a Sage/Python tuple made of:
816
      - a Sage/Python list of Sollya objects,
817
      - a Sage/Python int holding the number of elements,
818
      - a Sage/Python int stating (!= 0) that the list is end-elliptic.
819
    NOTE::
820
        We recover the addresses of the Sollya object from the list of pointers
821
        returned by sollya_lib_get_list_elements. The list itself is freed.
822
    TODO::
823
        Figure out what to do with numElements since the number of elements
824
        can easily be recovered from the list itself. 
825
        Ditto for isEndElliptic.
826
    """
827
    listAddress = POINTER(c_longlong)()
828
    numElements = c_int(0)
829
    isEndElliptic = c_int(0)
830
    listAsSageList = []
831
    result = sollya_lib_get_list_elements(byref(listAddress),\
832
                                          byref(numElements),\
833
                                          byref(isEndElliptic),\
834
                                          objectListSo)
835
    if result == 0 :
836
        return None
837
    for i in xrange(0, numElements.value, 1):
838
       #listAsSageList.append(sollya_lib_copy_obj(listAddress[i]))
839
       listAsSageList.append(listAddress[i])
840
       # Clear each of the elements returned by Sollya.
841
       #sollya_lib_clear_obj(listAddress[i])
842
    # Free the list itself.   
843
    sollya_lib_free(listAddress)
844
    return (listAsSageList, numElements.value, isEndElliptic.value)
845

    
846
def pobyso_get_max_prec_of_exp(soExp):
847
    """ Legacy function. See pobyso_get_max_prec_of_exp_so_sa. """
848
    return pobyso_get_max_prec_of_exp_so_sa(soExp)
849

    
850
def pobyso_get_max_prec_of_exp_so_sa(expSo):
851
    """
852
    Get the maximum precision used for the numbers in a Sollya expression.
853
    
854
    Arguments:
855
    soExp -- a Sollya expression pointer
856
    Return value:
857
    A Python integer
858
    TODO: 
859
    - error management;
860
    - correctly deal with numerical type such as DOUBLEEXTENDED.
861
    """
862
    if expSo is None:
863
        print inspect.stack()[0][3], ": expSo is None."
864
        return 0
865
    maxPrecision = 0
866
    minConstPrec = 0
867
    currentConstPrec = 0
868
    #pobyso_autoprint(expSo)
869
    operator = pobyso_get_head_function_so_sa(expSo)
870
    if (operator != SOLLYA_BASE_FUNC_CONSTANT) and \
871
    (operator != SOLLYA_BASE_FUNC_FREE_VARIABLE):
872
        (arity, subexpressions) = pobyso_get_subfunctions_so_sa(expSo)
873
        for i in xrange(arity):
874
            maxPrecisionCandidate = \
875
                pobyso_get_max_prec_of_exp_so_sa(subexpressions[i])
876
            if maxPrecisionCandidate > maxPrecision:
877
                maxPrecision = maxPrecisionCandidate
878
        return maxPrecision
879
    elif operator == SOLLYA_BASE_FUNC_CONSTANT:
880
        #minConstPrec = pobyso_get_min_prec_of_constant_so_sa(expSo)
881
        #currentConstPrec = pobyso_get_min_prec_of_constant_so_sa(soExp)
882
        #print minConstPrec, " - ", currentConstPrec 
883
        return pobyso_get_min_prec_of_constant_so_sa(expSo)
884
    
885
    elif operator == SOLLYA_BASE_FUNC_FREE_VARIABLE:
886
        return 0
887
    else:
888
        print "pobyso_get_max_prec_of_exp_so_sa: unexepected operator."
889
        return 0
890

    
891
def pobyso_get_min_prec_of_constant_so_sa(constExpSo):
892
    """
893
    Get the minimum precision necessary to represent the value of a Sollya
894
    constant.
895
    MPFR_MIN_PREC and powers of 2 are taken into account.
896
    We assume that constExpSo is a pointer to a Sollay constant expression.
897
    """
898
    constExpAsRnSa = pobyso_get_constant_as_rn_so_sa(constExpSo)
899
    return(min_mpfr_size(get_rn_value(constExpAsRnSa)))
900

    
901
def pobyso_get_poly_so_sa(polySo, realFieldSa=None):
902
    """
903
    Convert a Sollya polynomial into a Sage polynomial.
904
    Legacy function. Use pobyso_float_poly_so_sa() instead.
905
    """    
906
    return pobyso_float_poly_so_sa(polySo,realFieldSa)
907
# End pobyso_get_poly_so_sa
908

    
909
def pobyso_get_prec():
910
    """ Legacy function. See pobyso_get_prec_so_sa(). """
911
    return pobyso_get_prec_so_sa()
912

    
913
def pobyso_get_prec_so():
914
    """
915
    Get the current default precision in Sollya.
916
    The return value is a Sollya object.
917
    Usefull when modifying the precision back and forth by avoiding
918
    extra conversions.
919
    """
920
    return sollya_lib_get_prec(None)
921
    
922
def pobyso_get_prec_so_sa():
923
    """
924
    Get the current default precision in Sollya.
925
    The return value is Sage/Python int.
926
    """
927
    precSo = sollya_lib_get_prec()
928
    precSa = pobyso_constant_from_int_so_sa(precSo)
929
    sollya_lib_clear_obj(precSo)
930
    return precSa
931
# End pobyso_get_prec_so_sa.
932

    
933
def pobyso_get_prec_so_so_sa():
934
    """
935
    Return the current precision both as a Sollya object and a
936
    Sage integer as hybrid tuple.
937
    To avoid multiple calls for precision manipulations. 
938
    """
939
    precSo = sollya_lib_get_prec()
940
    precSa = pobyso_constant_from_int_so_sa(precSo)
941
    return (precSo, int(precSa))
942
    
943
def pobyso_get_prec_of_constant(ctExpSo):
944
    """ Legacy function. See pobyso_get_prec_of_constant_so_sa. """
945
    return pobyso_get_prec_of_constant_so_sa(ctExpSo)
946

    
947
def pobyso_get_prec_of_constant_so_sa(ctExpSo):
948
    """
949
    Tries to find a precision to represent ctExpSo without rounding.
950
    If not possible, returns None.
951
    """
952
    #print "Entering pobyso_get_prec_of_constant_so_sa..."
953
    prec = c_int(0)
954
    retc = sollya_lib_get_prec_of_constant(byref(prec), ctExpSo, None)
955
    if retc == 0:
956
        #print "pobyso_get_prec_of_constant_so_sa failed."
957
        return None
958
    #print "...exiting pobyso_get_prec_of_constant_so_sa."
959
    return int(prec.value)
960

    
961
def pobyso_get_prec_of_range_so_sa(rangeSo):
962
    """
963
    Returns the number of bits elements of a range are coded with.
964
    """
965
    prec = c_int(0)
966
    retc = sollya_lib_get_prec_of_range(byref(prec), rangeSo, None)
967
    if retc == 0:
968
        return(None)
969
    return int(prec.value)
970
# End pobyso_get_prec_of_range_so_sa()
971

    
972
def pobyso_get_sage_exp_from_sollya_exp(sollyaExpSo, realField = RR):
973
    """ Legacy function. See pobyso_get_sage_exp_from_sollya_exp_so_sa. """
974
    return pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, 
975
                                                     realField = RR)
976

    
977
def pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, realFieldSa = RR):
978
    """
979
    Get a Sage expression from a Sollya expression. 
980
    Currently only tested with polynomials with floating-point coefficients.
981
    Notice that, in the returned polynomial, the exponents are RealNumbers.
982
    """
983
    #pobyso_autoprint(sollyaExp)
984
    operatorSa = pobyso_get_head_function_so_sa(sollyaExpSo)
985
    sollyaLibFreeVariableName = sollya_lib_get_free_variable_name()
986
    ## Get rid of the "_"'s in "_x_", if any.
987
    sollyaLibFreeVariableName = re.sub('_', '', sollyaLibFreeVariableName)
988
    # Constants and the free variable are special cases.
989
    # All other operator are dealt with in the same way.
990
    if (operatorSa != SOLLYA_BASE_FUNC_CONSTANT) and \
991
       (operatorSa != SOLLYA_BASE_FUNC_FREE_VARIABLE):
992
        (aritySa, subexpressionsSa) = pobyso_get_subfunctions_so_sa(sollyaExpSo)
993
        if aritySa == 1:
994
            sageExpSa = eval(pobyso_function_type_as_string_so_sa(operatorSa) + \
995
            "(" + pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], \
996
            realFieldSa) + ")")
997
        elif aritySa == 2:
998
            # We do not get through the preprocessor.
999
            # The "^" operator is then a special case.
1000
            if operatorSa == SOLLYA_BASE_FUNC_POW:
1001
                operatorAsStringSa = "**"
1002
            else:
1003
                operatorAsStringSa = \
1004
                    pobyso_function_type_as_string_so_sa(operatorSa)
1005
            sageExpSa = \
1006
              eval("pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], realFieldSa)"\
1007
              + " " + operatorAsStringSa + " " + \
1008
                   "pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[1], realFieldSa)")
1009
        # We do not know yet how to deal with arity >= 3 
1010
        # (is there any in Sollya anyway?).
1011
        else:
1012
            sageExpSa = eval('None')
1013
        return sageExpSa
1014
    elif operatorSa == SOLLYA_BASE_FUNC_CONSTANT:
1015
        #print "This is a constant"
1016
        return pobyso_get_constant_as_rn_with_rf_so_sa(sollyaExpSo, realFieldSa)
1017
    elif operatorSa == SOLLYA_BASE_FUNC_FREE_VARIABLE:
1018
        #print "This is the free variable"
1019
        return eval(sollyaLibFreeVariableName)
1020
    else:
1021
        print "Unexpected"
1022
        return eval('None')
1023
# End pobyso_get_sage_exp_from_sollya_exp_so_sa
1024

    
1025

    
1026
def pobyso_get_subfunctions(expressionSo):
1027
    """ Legacy function. See pobyso_get_subfunctions_so_sa. """
1028
    return pobyso_get_subfunctions_so_sa(expressionSo) 
1029
# End pobyso_get_subfunctions.
1030
 
1031
def pobyso_get_subfunctions_so_sa(expressionSo):
1032
    """
1033
    Get the subfunctions of an expression.
1034
    Return the number of subfunctions and the list of subfunctions addresses.
1035
    S.T.: Could not figure out another way than that ugly list of declarations
1036
    to recover the addresses of the subfunctions.
1037
    We limit ourselves to arity 8 functions. 
1038
    """
1039
    subf0 = c_int(0)
1040
    subf1 = c_int(0)
1041
    subf2 = c_int(0)
1042
    subf3 = c_int(0)
1043
    subf4 = c_int(0)
1044
    subf5 = c_int(0)
1045
    subf6 = c_int(0)
1046
    subf7 = c_int(0)
1047
    subf8 = c_int(0)
1048
    arity = c_int(0)
1049
    nullPtr = POINTER(c_int)()
1050
    sollya_lib_get_subfunctions(expressionSo, byref(arity), \
1051
      byref(subf0), byref(subf1), byref(subf2), byref(subf3), \
1052
      byref(subf4), byref(subf5),\
1053
      byref(subf6), byref(subf7), byref(subf8), nullPtr, None) 
1054
#    byref(cast(subfunctions[0], POINTER(c_int))), \
1055
#    byref(cast(subfunctions[0], POINTER(c_int))), \
1056
#    byref(cast(subfunctions[2], POINTER(c_int))), \
1057
#    byref(cast(subfunctions[3], POINTER(c_int))), \
1058
#    byref(cast(subfunctions[4], POINTER(c_int))), \
1059
#    byref(cast(subfunctions[5], POINTER(c_int))), \
1060
#    byref(cast(subfunctions[6], POINTER(c_int))), \
1061
#    byref(cast(subfunctions[7], POINTER(c_int))), \
1062
#    byref(cast(subfunctions[8], POINTER(c_int))), nullPtr)
1063
    subfunctions = [subf0, subf1, subf2, subf3, subf4, subf5, subf6, subf7, \
1064
                    subf8]
1065
    subs = []
1066
    if arity.value > pobyso_max_arity:
1067
        return(0,[])
1068
    for i in xrange(arity.value):
1069
        subs.append(int(subfunctions[i].value))
1070
        #print subs[i]
1071
    return (int(arity.value), subs)
1072
# End pobyso_get_subfunctions_so_sa
1073
    
1074
def pobyso_guess_degree_sa_sa(functionSa, intervalSa, approxErrorSa, 
1075
                              weightSa=None, degreeBoundSa=None):
1076
    """
1077
    Sa_sa variant of the solly_guessdegree function.
1078
    Return 0 if something goes wrong.
1079
    """
1080
    functionAsStringSa = functionSa._assume_str().replace('_SAGE_VAR_', '')
1081
    functionSo = pobyso_parse_string_sa_so(functionAsStringSa)
1082
    if pobyso_is_error_so_sa(functionSo):
1083
        sollya_lib_clear_obj(functionSo)
1084
        return 0
1085
    rangeSo = pobyso_interval_to_range_sa_so(intervalSa)
1086
    # The approximation error is expected to be a floating point number.
1087
    if pobyso_is_floating_point_number_sa_sa(approxErrorSa):
1088
        approxErrorSo = pobyso_constant_sa_so(approxErrorSa)
1089
    else:
1090
        approxErrorSo = pobyso_constant_sa_so(RR(approxErrorSa))
1091
    if not weightSa is None:
1092
        weightAsStringSa = weightSa._assume_str().replace('_SAGE_VAR_', '')
1093
        weightSo = pobyso_parse_string_sa_so(weightAsStringSa)
1094
        if pobyso_is_error_so_sa(weightSo):
1095
            sollya_lib_clear_obj(functionSo)
1096
            sollya_lib_clear_obj(rangeSo)
1097
            sollya_lib_clear_obj(approxErrorSo)   
1098
            sollya_lib_clear_obj(weightSo)
1099
            return 0   
1100
    else:
1101
        weightSo = None
1102
    if not degreeBoundSa is None:
1103
        degreeBoundSo = pobyso_constant_from_int_sa_so(degreeBoundSa)
1104
    else:
1105
        degreeBoundSo = None
1106
    guessedDegreeSa = pobyso_guess_degree_so_sa(functionSo,
1107
                                                rangeSo,
1108
                                                approxErrorSo,
1109
                                                weightSo,
1110
                                                degreeBoundSo)
1111
    sollya_lib_clear_obj(functionSo)
1112
    sollya_lib_clear_obj(rangeSo)
1113
    sollya_lib_clear_obj(approxErrorSo)
1114
    if not weightSo is None:
1115
        sollya_lib_clear_obj(weightSo)
1116
    if not degreeBoundSo is None:
1117
        sollya_lib_clear_obj(degreeBoundSo)
1118
    return guessedDegreeSa
1119
# End poyso_guess_degree_sa_sa
1120

    
1121
def pobyso_guess_degree_so_sa(functionSo, rangeSo, errorSo, weightSo=None, \
1122
                              degreeBoundSo=None):
1123
    """
1124
    Thin wrapper around the guessdegree function.
1125
    Nevertheless, some precision control stuff has been appended.
1126
    """
1127
    # Deal with Sollya internal precision issues: if it is too small, 
1128
    # compared with the error, increases it to about twice -log2(error).
1129
    errorSa = pobyso_get_constant_as_rn_with_rf_so_sa(errorSo)
1130
    log2ErrorSa = errorSa.log2()
1131
    if log2ErrorSa < 0:
1132
        neededPrecisionSa = int(2 * int(-log2ErrorSa) / 64) * 64
1133
    else:
1134
        neededPrecisionSa = int(2 * int(log2ErrorSa) / 64) * 64
1135
    #print "Needed precision:", neededPrecisionSa
1136
    sollyaPrecisionHasChanged = False
1137
    (initialPrecSo, initialPrecSa) = pobyso_get_prec_so_so_sa()
1138
    if neededPrecisionSa > initialPrecSa:
1139
        if neededPrecisionSa <= 2:
1140
            print inspect.stack()[0][3], ": precision change <= 2 requested."
1141
        pobyso_set_prec_sa_so(neededPrecisionSa)
1142
        sollyaPrecisionHasChanged = True
1143
    #print "Guessing degree..."
1144
    # weightSo and degreeBoundsSo are optional arguments.
1145
    # As declared, sollya_lib_guessdegree must take 5 arguments.
1146
    if weightSo is None:
1147
        degreeRangeSo = sollya_lib_guessdegree(functionSo, rangeSo, errorSo, 
1148
                                               0, 0, None)
1149
    elif degreeBoundSo is None:
1150
        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, \
1151
                                                errorSo, weightSo, 0, None)
1152
    else:
1153
        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, errorSo, \
1154
                                                weightSo, degreeBoundSo, None)
1155
    #print "...degree guess done."
1156
    # Restore internal precision, if applicable.
1157
    if sollyaPrecisionHasChanged:
1158
        pobyso_set_prec_so_so(initialPrecSo)
1159
        sollya_lib_clear_obj(initialPrecSo)
1160
    degreeIntervalSa = pobyso_range_to_interval_so_sa(degreeRangeSo)
1161
    sollya_lib_clear_obj(degreeRangeSo)
1162
    # When ok, both bounds match.
1163
    # When the degree bound is too low, the upper bound is the degree
1164
    # for which the error can be honored.
1165
    # When it really goes wrong, the upper bound is infinity. 
1166
    if degreeIntervalSa.lower() == degreeIntervalSa.upper():
1167
        return int(degreeIntervalSa.lower())
1168
    else:
1169
        if degreeIntervalSa.upper().is_infinity():
1170
            return None
1171
        else:
1172
            return int(degreeIntervalSa.upper())
1173
    # End pobyso_guess_degree_so_sa
1174

    
1175
def pobyso_inf_so_so(intervalSo):
1176
    """
1177
    Very thin wrapper around sollya_lib_inf().
1178
    """
1179
    return sollya_lib_inf(intervalSo)
1180
# End pobyso_inf_so_so.
1181
    
1182
def pobyso_infnorm_so_so(func, interval, file = None, intervalList = None):
1183
    print "Do not use this function. User pobyso_supnorm_so_so instead."
1184
    return None
1185

    
1186
def pobyso_interval_to_range_sa_so(intervalSa, precisionSa=None):
1187
    if precisionSa is None:
1188
        precisionSa = intervalSa.parent().precision()
1189
    intervalSo = pobyso_bounds_to_range_sa_so(intervalSa.lower(),\
1190
                                              intervalSa.upper(),\
1191
                                              precisionSa)
1192
    return intervalSo
1193
# End pobyso_interval_to_range_sa_so
1194

    
1195
def pobyso_is_error_so_sa(objSo):
1196
    """
1197
    Thin wrapper around the sollya_lib_obj_is_error() function.
1198
    """
1199
    if sollya_lib_obj_is_error(objSo) != 0:
1200
        return True
1201
    else:
1202
        return False
1203
# End pobyso_is_error-so_sa
1204

    
1205
def pobyso_is_floating_point_number_sa_sa(numberSa):
1206
    """
1207
    Check whether a Sage number is floating point.
1208
    Exception stuff added because numbers other than
1209
    floating-point ones do not have the is_real() attribute.
1210
    """
1211
    try:
1212
        return numberSa.is_real()
1213
    except AttributeError:
1214
        return False
1215
# End pobyso_is_floating_piont_number_sa_sa
1216

    
1217
def pobyso_is_range_so_sa(rangeCandidateSo):
1218
    """
1219
    Thin wrapper over sollya_lib_is_range.
1220
    """
1221
    return sollya_lib_obj_is_range(rangeCandidateSo) != 0
1222
    
1223
# End pobyso_is_range_so_sa
1224

    
1225

    
1226
def pobyso_lib_init():
1227
    sollya_lib_init(None)
1228

    
1229
def pobyso_lib_close():
1230
    sollya_lib_close(None)
1231

    
1232
def pobyso_name_free_variable(freeVariableNameSa):
1233
    """ Legacy function. See pobyso_name_free_variable_sa_so. """
1234
    pobyso_name_free_variable_sa_so(freeVariableNameSa)
1235

    
1236
def pobyso_name_free_variable_sa_so(freeVariableNameSa):
1237
    """
1238
    Set the free variable name in Sollya from a Sage string.
1239
    """
1240
    sollya_lib_name_free_variable(freeVariableNameSa)
1241

    
1242
def pobyso_parse_string(string):
1243
    """ Legacy function. See pobyso_parse_string_sa_so. """
1244
    return pobyso_parse_string_sa_so(string)
1245
 
1246
def pobyso_parse_string_sa_so(string):
1247
    """
1248
    Get the Sollya expression computed from a Sage string or
1249
    a Sollya error object if parsing failed.
1250
    """
1251
    return sollya_lib_parse_string(string)
1252

    
1253
def pobyso_precision_so_sa(ctExpSo):
1254
    """
1255
    Computes the necessary precision to represent a number.
1256
    If x is not zero, it can be uniquely written as x = m · 2e 
1257
    where m is an odd integer and e is an integer. 
1258
    precision(x) returns the number of bits necessary to write m 
1259
    in binary (i.e. ceil(log2(m))).
1260
    """
1261
    #TODO: take care of the special case: 0, @NaN@, @Inf@
1262
    precisionSo = sollya_lib_precision(ctExpSo)
1263
    precisionSa = pobyso_constant_from_int_so_sa(precisionSo)
1264
    sollya_lib_clear_obj(precisionSo)
1265
    return precisionSa
1266
# End pobyso_precision_so_sa
1267

    
1268
def pobyso_polynomial_coefficients_progressive_round_so_so(polySo,
1269
                                                           funcSo,
1270
                                                           icSo,
1271
                                                           intervalSo,
1272
                                                           itpSo,
1273
                                                           ftpSo,
1274
                                                           maxPrecSo,
1275
                                                           maxErrSo,
1276
                                                           debug=False):
1277
    if debug:
1278
        print "Input arguments:"
1279
        pobyso_autoprint(polySo)
1280
        pobyso_autoprint(funcSo)
1281
        pobyso_autoprint(icSo)
1282
        pobyso_autoprint(intervalSo)
1283
        pobyso_autoprint(itpSo)
1284
        pobyso_autoprint(ftpSo)
1285
        pobyso_autoprint(maxPrecSo)
1286
        pobyso_autoprint(maxErrSo)
1287
        print "________________"
1288
    
1289
    ## Higher order function see: 
1290
    #  http://effbot.org/pyfaq/how-do-you-make-a-higher-order-function-in-python.htm
1291
    def precision_decay_ratio_function(degreeSa):
1292
        def outer(x):
1293
            def inner(x):
1294
                we = 3/8
1295
                wq = 2/8
1296
                a  = 2.2
1297
                b  = 2 
1298
                return we*(exp(x/a)-1) +  wq*((b*x)**2) + (1-we-wq)*x
1299
            return  inner(x)/inner(degreeSa)
1300
        return outer
1301
    
1302
    #
1303
    degreeSa        = pobyso_polynomial_degree_so_sa(polySo)
1304
    ratio           = precision_decay_ratio_function(degreeSa)
1305
    itpSa           = pobyso_constant_from_int_so_sa(itpSo)
1306
    ftpSa           = pobyso_constant_from_int_so_sa(ftpSo)
1307
    maxPrecSa       = pobyso_constant_from_int_so_sa(maxPrecSo)
1308
    maxErrSa        = pobyso_get_constant_as_rn_so_sa(maxErrSo)
1309
    if debug:
1310
        print "degreeSa:", degreeSa
1311
        print "ratio:", ratio
1312
        print "itpsSa:", itpSa
1313
        print "ftpSa:", ftpSa
1314
        print "maxPrecSa:", maxPrecSa
1315
        print "maxErrSa:", maxErrSa
1316
    lastResPolySo   = None
1317
    lastInfNormSo   = None
1318
    #print "About to enter the while loop..."
1319
    while True: 
1320
        resPolySo   = pobyso_constant_0_sa_so()
1321
        pDeltaSa    = ftpSa - itpSa
1322
        for indexSa in reversed(xrange(0,degreeSa+1)):
1323
            #print "Index:", indexSa
1324
            indexSo = pobyso_constant_from_int_sa_so(indexSa)
1325
            #pobyso_autoprint(indexSo)
1326
            #print ratio(indexSa)
1327
            ctpSa = floor(ftpSa - (pDeltaSa * ratio(indexSa)))
1328
            ctpSo = pobyso_constant_from_int_sa_so(ctpSa)
1329
            if debug:
1330
                print "Index:", indexSa, " - Target precision:", 
1331
                pobyso_autoprint(ctpSo)
1332
            cmonSo  = \
1333
                sollya_lib_build_function_mul(sollya_lib_coeff(polySo, indexSo),
1334
                                      sollya_lib_build_function_pow( \
1335
                                          sollya_lib_build_function_free_variable(), \
1336
                                          indexSo))
1337
            #pobyso_autoprint(cmonSo)
1338
            cmonrSo = pobyso_round_coefficients_single_so_so(cmonSo, ctpSo)
1339
            sollya_lib_clear_obj(cmonSo)
1340
            #pobyso_autoprint(cmonrSo)
1341
            resPolySo = sollya_lib_build_function_add(resPolySo,
1342
                                                      cmonrSo)
1343
            #pobyso_autoprint(resPolySo)
1344
        # End for index
1345
        freeVarSo     = sollya_lib_build_function_free_variable()
1346
        changeVarSo   = sollya_lib_sub(freeVarSo, icSo)
1347
        resPolyCvSo   = sollya_lib_evaluate(resPolySo, changeVarSo)
1348
        errFuncSo = sollya_lib_build_function_sub(sollya_lib_copy_obj(funcSo), 
1349
                                                  resPolyCvSo)
1350
        infNormSo = sollya_lib_dirtyinfnorm(errFuncSo, intervalSo)
1351
        cerrSa    = pobyso_get_constant_as_rn_so_sa(infNormSo)
1352
        if debug:
1353
            print "Infnorm (Sollya):", 
1354
            pobyso_autoprint(infNormSo)
1355
        sollya_lib_clear_obj(errFuncSo)
1356
        #print "Infnorm  (Sage):", cerrSa
1357
        if (cerrSa > maxErrSa):
1358
            if debug:
1359
                print "Error is too large."
1360
            if lastResPolySo is None:
1361
                if debug:
1362
                    print "Enlarging prec."
1363
                ntpSa = floor(ftpSa + ftpSa/50)
1364
                ## Can't enlarge (numerical)
1365
                if ntpSa == ftpSa:
1366
                    sollya_lib_clear_obj(resPolySo)
1367
                    return None
1368
                ## Can't enlarge (not enough precision left)
1369
                if ntpSa > maxPrecSa:
1370
                    sollya_lib_clear_obj(resPolySo)
1371
                    return None
1372
                ftpSa = ntpSa
1373
                continue
1374
            ## One enlargement took place.
1375
            else:
1376
                if debug:
1377
                    print "Exit with the last before last polynomial."
1378
                    print "Precision of highest degree monomial:", itpSa
1379
                    print "Precision of constant term          :", ftpSa
1380
                sollya_lib_clear_obj(resPolySo)
1381
                sollya_lib_clear_obj(infNormSo)
1382
                return (lastResPolySo, lastInfNormSo)
1383
        # cerrSa <= maxErrSa: scrap more bits, possibly.
1384
        else:
1385
            if debug:
1386
                print "Error is too small"
1387
            if cerrSa <= (maxErrSa/2):
1388
                if debug: 
1389
                    print "Shrinking prec."
1390
                ntpSa = floor(ftpSa - ftpSa/50)
1391
                ## Can't shrink (numerical)
1392
                if ntpSa == ftpSa:
1393
                    if not lastResPolySo is None:
1394
                        sollya_lib_clear_obj(lastResPolySo)
1395
                    if not lastInfNormSo is None:
1396
                        sollya_lib_clear_obj(lastInfNormSo)
1397
                    if debug:
1398
                        print "Exit because can't shrink anymore (numerically)."
1399
                        print "Precision of highest degree monomial:", itpSa
1400
                        print "Precision of constant term          :", ftpSa
1401
                    return (resPolySo, infNormSo)
1402
                ## Can't shrink (not enough precision left)
1403
                if ntpSa <= itpSa:
1404
                    if not lastResPolySo is None:
1405
                        sollya_lib_clear_obj(lastResPolySo)
1406
                    if not lastInfNormSo is None:
1407
                        sollya_lib_clear_obj(lastInfNormSo)
1408
                        print "Exit because can't shrink anymore (no bits left)."
1409
                        print "Precision of highest degree monomial:", itpSa
1410
                        print "Precision of constant term          :", ftpSa
1411
                    return (resPolySo, infNormSo)
1412
                ftpSa = ntpSa
1413
                if not lastResPolySo is None:
1414
                    sollya_lib_clear_obj(lastResPolySo)
1415
                if not lastInfNormSo is None:
1416
                    sollya_lib_clear_obj(lastInfNormSo)
1417
                lastResPolySo = resPolySo
1418
                lastInfNormSo = infNormSo
1419
                continue
1420
            else: # Error is not that small, just return
1421
                if not lastResPolySo is None:
1422
                    sollya_lib_clear_obj(lastResPolySo)
1423
                if not lastInfNormSo is None:
1424
                    sollya_lib_clear_obj(lastInfNormSo)
1425
                if debug:
1426
                    print "Exit normally."
1427
                    print "Precision of highest degree monomial:", itpSa
1428
                    print "Precision of constant term          :", ftpSa
1429
                return (resPolySo, infNormSo)                
1430
    # End wile True
1431
    return None
1432
# End pobyso_polynomial_coefficients_progressive_truncate_so_so.
1433

    
1434
def pobyso_polynomial_degree_so_sa(polySo):
1435
    """
1436
    Return the degree of a Sollya polynomial as a Sage int.
1437
    """
1438
    degreeSo = sollya_lib_degree(polySo)
1439
    return pobyso_constant_from_int_so_sa(degreeSo)
1440
# End pobyso_polynomial_degree_so_sa
1441

    
1442
def pobyso_polynomial_degree_so_so(polySo):
1443
    """
1444
    Thin wrapper around lib_sollya_degree().
1445
    """
1446
    return sollya_lib_degree(polySo)
1447
# End pobyso_polynomial_degree_so_so
1448
                                                                  
1449
def pobyso_range(rnLowerBound, rnUpperBound):
1450
    """ Legacy function. See pobyso_range_sa_so. """
1451
    return pobyso_range_sa_so(rnLowerBound, rnUpperBound) 
1452

    
1453
def pobyso_range_from_bounds_sa_so(rnLowerBound, rnUpperBound, precSa = None):
1454
    """
1455
    Create a Sollya range from 2 Sage real numbers as bounds
1456
    """
1457
    # TODO check precision stuff.
1458
    sollyaPrecChanged = False
1459
    (initialSollyaPrecSo, initialSollyaPrecSa) = \
1460
        pobyso_get_prec_so_so_sa()
1461
    if precSa is None:
1462
        precSa = max(rnLowerBound.parent().prec(), rnUpperBound.parent().prec())
1463
    if precSa > initialSollyaPrecSa:
1464
        if precSa <= 2:
1465
            print inspect.stack()[0][3], ": precision change <= 2 requested."
1466
        pobyso_set_prec_sa_so(precSa)
1467
        sollyaPrecChanged = True
1468
    rangeSo = sollya_lib_range_from_bounds(get_rn_value(rnLowerBound),
1469
                                           get_rn_value(rnUpperBound))
1470
    if sollyaPrecChanged:
1471
        pobyso_set_prec_so_so(initialSollyaPrecSo)
1472
        sollya_lib_clear_obj(initialSollyaPrecSo)
1473
    return rangeSo
1474
# End pobyso_range_from_bounds_sa_so
1475

    
1476
def pobyso_range_list_so_sa(listSo):
1477
    """
1478
    Return a Sollya list of ranges as a Sage list of 
1479
    floating-point intervals.
1480
    """
1481
    listSa   = []
1482
    ## The function returns none if the list is empty or an error has happened.
1483
    retVal = pobyso_get_list_elements_so_so(listSo)
1484
    if retVal is None:
1485
        return listSa
1486
    ## Just in case the interface is changed and an empty list is returned 
1487
    #  instead of None.
1488
    elif len(retVal) == 0:
1489
        return listSa
1490
    else:
1491
        ## Remember pobyso_get_list_elements_so_so returns more information
1492
        #  than just the elements of the list (# elements, is_elliptic)
1493
        listSaSo, numElements, isEndElliptic = retVal
1494
    ## Return an empty list.
1495
    if numElements == 0:
1496
        return listSa
1497
    ## Search first for the maximum precision of the elements
1498
    maxPrecSa = 0
1499
    for rangeSo in listSaSo:
1500
        #pobyso_autoprint(floatSo)
1501
        curPrecSa =  pobyso_get_prec_of_range_so_sa(rangeSo)
1502
        if curPrecSa > maxPrecSa:
1503
            maxPrecSa = curPrecSa
1504
    ##
1505
    intervalField = RealIntervalField(maxPrecSa)
1506
    ##
1507
    for rangeSo in listSaSo:
1508
        listSa.append(pobyso_range_to_interval_so_sa(rangeSo, intervalField))
1509
    return listSa
1510
# End pobyso_range_list_so_sa
1511

    
1512
def pobyso_range_max_abs_so_so(rangeSo):
1513
    """
1514
    Return, as Sollya constant, the maximum absolute value of a Sollay range.
1515
    """
1516
    lowerBoundSo = sollya_lib_inf(rangeSo)
1517
    upperBoundSo = sollya_lib_sup(rangeSo)
1518
    #
1519
    lowerBoundSo = sollya_lib_build_function_abs(lowerBoundSo)
1520
    upperBoundSo = sollya_lib_build_function_abs(upperBoundSo)
1521
    #pobyso_autoprint(lowerBoundSo)
1522
    #pobyso_autoprint(upperBoundSo)
1523
    #
1524
    maxAbsSo = sollya_lib_max(lowerBoundSo, upperBoundSo, None)
1525
    #sollya_lib_clear_obj(lowerBoundSo)
1526
    #sollya_lib_clear_obj(upperBoundSo)
1527
    return maxAbsSo
1528
# End pobyso_range_max_abs_so_so
1529
 
1530
def pobyso_range_to_interval_so_sa(rangeSo, realIntervalFieldSa = None):
1531
    """
1532
    Get a Sage interval from a Sollya range.
1533
    If no realIntervalField is given as a parameter, the Sage interval
1534
    precision is that of the Sollya range.
1535
    Otherwise, the precision is that of the realIntervalField. In this case
1536
    rounding may happen.
1537
    """
1538
    if realIntervalFieldSa is None:
1539
        precSa = pobyso_get_prec_of_range_so_sa(rangeSo)
1540
        realIntervalFieldSa = RealIntervalField(precSa)
1541
    intervalSa = \
1542
        pobyso_get_interval_from_range_so_sa(rangeSo, realIntervalFieldSa) 
1543
    return intervalSa
1544
# End pobyso_range_to_interval_so_sa
1545
#
1546
def pobyso_relative_so_so():
1547
    """
1548
    Very thin wrapper around the sollya_lib_relative function.
1549
    """
1550
    return sollya_lib_relative()
1551
# End pobyso_relative_so_so
1552
#
1553
def pobyso_rat_poly_sa_so(polySa, precSa = None):
1554
    """
1555
    Create a Sollya polynomial from a Sage rational polynomial.
1556
    """
1557
    ## TODO: filter arguments.
1558
    ## Precision. If no precision is given, use the current precision
1559
    #  of Sollya.
1560
    if precSa is None:
1561
        precSa =  pobyso_get_prec_so_sa()
1562
    #print "Precision:",  precSa
1563
    RRR = RealField(precSa)
1564
    ## Create a Sage polynomial in the "right" precision.
1565
    P_RRR = RRR[polySa.variables()[0]]
1566
    polyFloatSa = P_RRR(polySa)
1567
    ## Make sure no precision is provided: pobyso_float_poly_sa_so will
1568
    #  recover it all by itself and not make an extra conversion.
1569
    return pobyso_float_poly_sa_so(polyFloatSa)
1570
    
1571
# End pobyso_rat_poly_sa_so    
1572
    
1573
def pobyso_remez_canonical_sa_sa(func, \
1574
                                 degree, \
1575
                                 lowerBound, \
1576
                                 upperBound, \
1577
                                 weight = None, \
1578
                                 quality = None):
1579
    """
1580
    All arguments are Sage/Python.
1581
    The functions (func and weight) must be passed as expressions or strings.
1582
    Otherwise the function fails. 
1583
    The return value is a Sage polynomial.
1584
    """
1585
    var('zorglub')    # Dummy variable name for type check only. Type of 
1586
    # zorglub is "symbolic expression".
1587
    polySo = pobyso_remez_canonical_sa_so(func, \
1588
                                 degree, \
1589
                                 lowerBound, \
1590
                                 upperBound, \
1591
                                 weight, \
1592
                                 quality)
1593
    # String test
1594
    if parent(func) == parent("string"):
1595
        functionSa = eval(func)
1596
    # Expression test.
1597
    elif type(func) == type(zorglub):
1598
        functionSa = func
1599
    else:
1600
        return None
1601
    #
1602
    maxPrecision = 0
1603
    if polySo is None:
1604
        return(None)
1605
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
1606
    RRRRSa = RealField(maxPrecision)
1607
    polynomialRingSa = RRRRSa[functionSa.variables()[0]]
1608
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, RRRRSa)
1609
    polySa = polynomial(expSa, polynomialRingSa)
1610
    sollya_lib_clear_obj(polySo)
1611
    return(polySa)
1612
# End pobyso_remez_canonical_sa_sa
1613
    
1614
def pobyso_remez_canonical(func, \
1615
                           degree, \
1616
                           lowerBound, \
1617
                           upperBound, \
1618
                           weight = "1", \
1619
                           quality = None):
1620
    """ Legacy function. See pobyso_remez_canonical_sa_so. """
1621
    return(pobyso_remez_canonical_sa_so(func, \
1622
                                        degree, \
1623
                                        lowerBound, \
1624
                                        upperBound, \
1625
                                        weight, \
1626
                                        quality))
1627
# End pobyso_remez_canonical.
1628

    
1629
def pobyso_remez_canonical_sa_so(func, \
1630
                                 degree, \
1631
                                 lowerBound, \
1632
                                 upperBound, \
1633
                                 weight = None, \
1634
                                 quality = None):
1635
    """
1636
    All arguments are Sage/Python.
1637
    The functions (func and weight) must be passed as expressions or strings.
1638
    Otherwise the function fails. 
1639
    The return value is a pointer to a Sollya function.
1640
    lowerBound and upperBound mus be reals.
1641
    """
1642
    var('zorglub')    # Dummy variable name for type check only. Type of
1643
    # zorglub is "symbolic expression".
1644
    currentVariableNameSa = None
1645
    # The func argument can be of different types (string, 
1646
    # symbolic expression...)
1647
    if parent(func) == parent("string"):
1648
        localFuncSa = sage_eval(func,globals())
1649
        if len(localFuncSa.variables()) > 0:
1650
            currentVariableNameSa = localFuncSa.variables()[0]
1651
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1652
            functionSo = \
1653
              sollya_lib_parse_string(localFuncSa._assume_str().replace('_SAGE_VAR_', ''))
1654
    # Expression test.
1655
    elif type(func) == type(zorglub):
1656
        # Until we are able to translate Sage expressions into Sollya 
1657
        # expressions : parse the string version.
1658
        if len(func.variables()) > 0:
1659
            currentVariableNameSa = func.variables()[0]
1660
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1661
            functionSo = \
1662
              sollya_lib_parse_string(func._assume_str().replace('_SAGE_VAR_', ''))
1663
    else:
1664
        return(None)
1665
    if weight is None: # No weight given -> 1.
1666
        weightSo = pobyso_constant_1_sa_so()
1667
    elif parent(weight) == parent("string"): # Weight given as string: parse it.
1668
        weightSo = sollya_lib_parse_string(func)
1669
    elif type(weight) == type(zorglub): # Weight given as symbolice expression.
1670
        functionSo = \
1671
          sollya_lib_parse_string_sa_so(weight._assume_str().replace('_SAGE_VAR_', ''))
1672
    else:
1673
        return(None)
1674
    degreeSo = pobyso_constant_from_int(degree)
1675
    rangeSo = pobyso_bounds_to_range_sa_so(lowerBound, upperBound)
1676
    if not quality is None:
1677
        qualitySo= pobyso_constant_sa_so(quality)
1678
    else:
1679
        qualitySo = None
1680
        
1681
    remezPolySo = sollya_lib_remez(functionSo, \
1682
                                   degreeSo, \
1683
                                   rangeSo, \
1684
                                   weightSo, \
1685
                                   qualitySo, \
1686
                                   None)
1687
    sollya_lib_clear_obj(functionSo)
1688
    sollya_lib_clear_obj(degreeSo)
1689
    sollya_lib_clear_obj(rangeSo)
1690
    sollya_lib_clear_obj(weightSo)
1691
    if not qualitySo is None:
1692
        sollya_lib_clear_obj(qualitySo)
1693
    return(remezPolySo)
1694
# End pobyso_remez_canonical_sa_so
1695

    
1696
def pobyso_remez_canonical_so_so(funcSo, \
1697
                                 degreeSo, \
1698
                                 rangeSo, \
1699
                                 weightSo = pobyso_constant_1_sa_so(),\
1700
                                 qualitySo = None):
1701
    """
1702
    All arguments are pointers to Sollya objects.
1703
    The return value is a pointer to a Sollya function.
1704
    """
1705
    if not sollya_lib_obj_is_function(funcSo):
1706
        return(None)
1707
    return(sollya_lib_remez(funcSo, degreeSo, rangeSo, weightSo, qualitySo, None))
1708
# End pobyso_remez_canonical_so_so.
1709

    
1710
def pobyso_remez_exponents_list_sa_so(func,     \
1711
                                 exponentsList, \
1712
                                 lowerBound,    \
1713
                                 upperBound,    \
1714
                                 weight = None, \
1715
                                 quality = None):
1716
    """
1717
    All arguments are Sage/Python.
1718
    The functions (func and weight) must be passed as expressions or strings.
1719
    Otherwise the function fails. 
1720
    The return value is a pointer to a Sollya function.
1721
    lowerBound and upperBound mus be reals.
1722
    """
1723
    var('zorglub')    # Dummy variable name for type check only. Type of
1724
    # zorglub is "symbolic expression".
1725
    currentVariableNameSa = None
1726
    # The func argument can be of different types (string, 
1727
    # symbolic expression...)
1728
    if parent(func) == parent("string"):
1729
        localFuncSa = sage_eval(func,globals())
1730
        if len(localFuncSa.variables()) > 0:
1731
            currentVariableNameSa = localFuncSa.variables()[0]
1732
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1733
            functionSo = \
1734
              sollya_lib_parse_string(localFuncSa._assume_str().replace('_SAGE_VAR_', ''))
1735
    # Expression test.
1736
    elif type(func) == type(zorglub):
1737
        # Until we are able to translate Sage expressions into Sollya 
1738
        # expressions : parse the string version.
1739
        if len(func.variables()) > 0:
1740
            currentVariableNameSa = func.variables()[0]
1741
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1742
            functionSo = \
1743
              sollya_lib_parse_string(func._assume_str().replace('_SAGE_VAR_', ''))
1744
    else:
1745
        return(None)
1746
    ## Deal with the weight, much in the same way as with the function.
1747
    if weight is None: # No weight given -> 1.
1748
        weightSo = pobyso_constant_1_sa_so()
1749
    elif parent(weight) == parent("string"): # Weight given as string: parse it.
1750
        weightSo = sollya_lib_parse_string(func)
1751
    elif type(weight) == type(zorglub): # Weight given as symbolice expression.
1752
        functionSo = \
1753
          sollya_lib_parse_string(weight._assume_str().replace('_SAGE_VAR_', '',100))
1754
    else:
1755
        return(None)
1756
    rangeSo = pobyso_bounds_to_range_sa_so(lowerBound, upperBound)
1757
    if not quality is None:
1758
        qualitySo= pobyso_constant_sa_so(quality)
1759
    else:
1760
        qualitySo = None
1761
    #
1762
    ## Tranform the Sage list of exponents into a Sollya list.
1763
    exponentsListSo = pobyso_build_list_of_ints_sa_so(*exponentsList)
1764
    remezPolySo = sollya_lib_remez(functionSo, \
1765
                                   exponentsListSo, \
1766
                                   rangeSo, \
1767
                                   weightSo, \
1768
                                   qualitySo, \
1769
                                   None)
1770
    sollya_lib_clear_obj(functionSo)
1771
    sollya_lib_clear_obj(exponentsListSo)
1772
    sollya_lib_clear_obj(rangeSo)
1773
    sollya_lib_clear_obj(weightSo)
1774
    if not qualitySo is None:
1775
        sollya_lib_clear_obj(qualitySo)
1776
    return(remezPolySo)
1777
# End pobyso_remez_exponentsList_sa_so
1778
#
1779
def pobyso_round_coefficients_so_so(polySo, truncFormatListSo):
1780
    """
1781
    A wrapper around the "classical" sollya_lib_roundcoefficients: a Sollya
1782
    polynomial and a Sollya list are given as arguments.
1783
    """
1784
    return sollya_lib_roundcoefficients(polySo, truncFormatListSo)
1785

    
1786
def pobyso_round_coefficients_progressive_so_so(polySo, 
1787
                                                funcSo,
1788
                                                precSo,
1789
                                                intervalSo,
1790
                                                icSo,
1791
                                                currentApproxErrorSo,
1792
                                                approxAccurSo,
1793
                                                debug=False):
1794
    if debug:
1795
        print "Input arguments:"
1796
        print "Polynomial: ", ; pobyso_autoprint(polySo)
1797
        print "Function: ", ; pobyso_autoprint(funcSo)
1798
        print "Internal precision: ", ; pobyso_autoprint(precSo)
1799
        print "Interval: ", ; pobyso_autoprint(intervalSo)
1800
        print "Current approximation error: ", ; pobyso_autoprint(currentApproxErrorSo)
1801
        print "Requested approxiation error: ", ; pobyso_autoprint(approxAccurSo)
1802
        print "________________"
1803
    approxAccurSa        = pobyso_get_constant_as_rn_so_sa(approxAccurSo)
1804
    currentApproxErrorSa = pobyso_get_constant_as_rn_so_sa(currentApproxErrorSo)
1805
    ## If the current approximation error is too close to the target, there is
1806
    #  no possible gain.
1807
    if currentApproxErrorSa >= approxAccurSa / 2:
1808
        #### Do not return the initial argument but copies: caller may free 
1809
        #    as inutile after call.
1810
        return (sollya_lib_copy_obj(polySo),
1811
                sollya_lib_copy_obj(currentApproxErrorSo))
1812
    degreeSa             = pobyso_polynomial_degree_so_sa(polySo)
1813
    intervalSa           = pobyso_range_to_interval_so_sa(intervalSo)
1814
    
1815
    if debug:
1816
        print "degreeSa              :", degreeSa
1817
        print "intervalSa            :", intervalSa.str(style='brackets')
1818
        print "currentApproxErrorSa  :", currentApproxErrorSa 
1819
        print "approxAccurSa         :", approxAccurSa 
1820
    ### Start with a 0 value expression.
1821
    radiusSa = intervalSa.absolute_diameter() / 2
1822
    if debug:
1823
        print "log2(radius):", RR(radiusSa).log2()
1824
    iterIndex = 0
1825
    while True: 
1826
        resPolySo = pobyso_constant_0_sa_so()
1827
        roundedPolyApproxAccurSa = approxAccurSa / 2
1828
        currentRadiusPowerSa = 1 
1829
        for degree in xrange(0,degreeSa + 1):
1830
            #### At round 0, use the agressive formula. At round 1, run the
1831
            #    proved formula.
1832
            if iterIndex == 0:
1833
                roundingPowerSa = \
1834
                    floor(((currentRadiusPowerSa/roundedPolyApproxAccurSa)*(degree+1)).log2())
1835
            else:
1836
                roundingPowerSa = \
1837
                    floor(((currentRadiusPowerSa/roundedPolyApproxAccurSa)*(degreeSa+1)).log2())
1838
            ## Under extreme conditions the above formulas can evaluate under 2, which is the
1839
            #  minimal precision of an MPFR number.
1840
            if roundingPowerSa < 2:
1841
                roundingPowerSa = 2
1842
            if debug:
1843
                print "roundedPolyApproxAccurSa", roundedPolyApproxAccurSa
1844
                print "currentRadiusPowerSa", currentRadiusPowerSa
1845
                print "Current rounding exponent:", roundingPowerSa
1846
            currentRadiusPowerSa *= radiusSa
1847
            index1So = pobyso_constant_from_int_sa_so(degree)
1848
            index2So = pobyso_constant_from_int_sa_so(degree)
1849
            ### Create a monomial with:
1850
            #   - the coefficient in the initial monomial at the current degrree;
1851
            #   - the current exponent;
1852
            #   - the free variable.
1853
            cmonSo  = \
1854
                sollya_lib_build_function_mul(sollya_lib_coeff(polySo, index1So),
1855
                                              sollya_lib_build_function_pow( \
1856
                                                sollya_lib_build_function_free_variable(), \
1857
                                                index2So))
1858
            roundingPowerSo = pobyso_constant_from_int_sa_so(roundingPowerSa)
1859
            cmonrSo = pobyso_round_coefficients_single_so_so(cmonSo, roundingPowerSo)
1860
            sollya_lib_clear_obj(cmonSo)
1861
            ### Add to the result polynomial.
1862
            resPolySo = sollya_lib_build_function_add(resPolySo,
1863
                                                      cmonrSo)
1864
        # End for.
1865
        ### Check the new polynomial.
1866
        freeVarSo     = sollya_lib_build_function_free_variable()
1867
        changeVarSo   = sollya_lib_sub(freeVarSo, icSo)
1868
        resPolyCvSo   = sollya_lib_evaluate(resPolySo, changeVarSo)
1869
        errFuncSo     = sollya_lib_build_function_sub(sollya_lib_copy_obj(funcSo), 
1870
                                                      resPolyCvSo)
1871
        infNormSo = sollya_lib_dirtyinfnorm(errFuncSo, intervalSo)
1872
        ### This also clears resPolyCvSo.
1873
        sollya_lib_clear_obj(errFuncSo)
1874
        cerrSa    = pobyso_get_constant_as_rn_so_sa(infNormSo)
1875
        if debug:
1876
            print "Error of the new polynomial:", cerrSa
1877
        ### If at round 1, return the initial polynomial error. This should
1878
        #   never happen since the rounding algorithm is proved. But some 
1879
        #   circumstances may break it (e.g. internal precision of tools).
1880
        if cerrSa > approxAccurSa:
1881
            if iterIndex > 0: # Round 1 and beyond.
1882
                sollya_lib_clear_obj(resPolySo)
1883
                sollya_lib_clear_obj(infNormSo)
1884
                #### Do not return the arguments but copies: the caller may free
1885
                #    free the former as inutile after call.
1886
                return (sollya_lib_copy_obj(polySo), 
1887
                        sollya_lib_copy_obj(currentApproxErrorSo))
1888
            else: # Round 0, got round 1
1889
                sollya_lib_clear_obj(resPolySo)
1890
                sollya_lib_clear_obj(infNormSo)
1891
                iterIndex += 1
1892
                continue
1893
        ### If get here it is because cerrSa <= approxAccurSa
1894
        ### Approximation error of the new polynomial is acceptable.
1895
        return (resPolySo, infNormSo)
1896
    # End while True
1897
# End pobyso_round_coefficients_progressive_so_so
1898
    
1899
def pobyso_round_coefficients_single_so_so(polySo, commonPrecSo):
1900
    """
1901
    Create a rounded coefficients polynomial from polynomial argument to
1902
    the number of bits in size argument.
1903
    All coefficients are set to the same precision.
1904
    """
1905
    ## TODO: check arguments.
1906
    endEllipListSo = pobyso_build_end_elliptic_list_so_so(commonPrecSo)
1907
    polySo = sollya_lib_roundcoefficients(polySo, endEllipListSo, None)
1908
    sollya_lib_clear_obj(endEllipListSo)
1909
    #sollya_lib_clear_obj(endEllipListSo)
1910
    return polySo
1911
    
1912
# End pobyso_round_coefficients_single_so_so
1913

    
1914
def pobyso_set_canonical_off():
1915
    sollya_lib_set_canonical(sollya_lib_off())
1916

    
1917
def pobyso_set_canonical_on():
1918
    sollya_lib_set_canonical(sollya_lib_on())
1919

    
1920
def pobyso_set_prec(p):
1921
    """ Legacy function. See pobyso_set_prec_sa_so. """
1922
    pobyso_set_prec_sa_so(p)
1923

    
1924
def pobyso_set_prec_sa_so(p):
1925
    #a = c_int(p)
1926
    #precSo = c_void_p(sollya_lib_constant_from_int(a))
1927
    #precSo = sollya_lib_constant_from_int(a)
1928
    precSo = pobyso_constant_from_int_sa_so(p)
1929
    sollya_lib_set_prec(precSo)
1930
    sollya_lib_clear_obj(precSo)
1931
# End pobyso_set_prec_sa_so.
1932

    
1933
def pobyso_set_prec_so_so(newPrecSo):
1934
    sollya_lib_set_prec(newPrecSo)
1935
# End pobyso_set_prec_so_so.
1936

    
1937
def pobyso_inf_so_so(intervalSo):
1938
    """
1939
    Very thin wrapper around sollya_lib_inf().
1940
    """
1941
    return sollya_lib_inf(intervalSo)
1942
# End pobyso_inf_so_so.
1943
#   
1944
def pobyso_supnorm_so_sa(polySo, funcSo, intervalSo, errorTypeSo = None,\
1945
                         accuracySo = None, realFieldSa = None):
1946
    """
1947
    Computes the supremum norm from Solly input arguments and returns a
1948
    Sage floating-point number whose precision is set by the only Sage argument.
1949
    
1950
    The returned value is the maximum of the absolute values of the range
1951
    elements  returned by the Sollya supnorm functions
1952
    """
1953
    supNormRangeSo = pobyso_supnorm_so_so(polySo, 
1954
                                          funcSo, 
1955
                                          intervalSo, 
1956
                                          errorTypeSo,
1957
                                          accuracySo)
1958
    supNormSo = pobyso_range_max_abs_so_so(supNormRangeSo)
1959
    sollya_lib_clear_obj(supNormRangeSo)
1960
    #pobyso_autoprint(supNormSo)
1961
    supNormSa = pobyso_get_constant_as_rn_with_rf_so_sa(supNormSo, realFieldSa)
1962
    sollya_lib_clear_obj(supNormSo)
1963
    return supNormSa 
1964
# End pobyso_supnorm_so_sa.
1965
#
1966
def pobyso_supnorm_so_so(polySo, funcSo, intervalSo, errorTypeSo = None,\
1967
                         accuracySo = None):
1968
    """
1969
    Computes the supnorm of the approximation error between the given 
1970
    polynomial and function. Attention: returns a range!
1971
    errorTypeSo defaults to "absolute".
1972
    accuracySo defaults to 2^(-40).
1973
    """
1974
    if errorTypeSo is None:
1975
        errorTypeSo = sollya_lib_absolute(None)
1976
        errorTypeIsNone = True
1977
    else:
1978
        errorTypeIsNone = False
1979
    #
1980
    if accuracySo is None:
1981
        # Notice the **: we are in Pythonland here!
1982
        accuracySo = pobyso_constant_sa_so(RR(2**(-40)))
1983
        accuracyIsNone = True
1984
    else:
1985
        accuracyIsNone = False
1986
    #pobyso_autoprint(accuracySo)
1987
    resultSo = \
1988
        sollya_lib_supnorm(polySo, funcSo, intervalSo, errorTypeSo, \
1989
                              accuracySo)
1990
    if errorTypeIsNone:
1991
        sollya_lib_clear_obj(errorTypeSo)
1992
    if accuracyIsNone:
1993
        sollya_lib_clear_obj(accuracySo)
1994
    return resultSo
1995
# End pobyso_supnorm_so_so
1996
#
1997
def pobyso_taylor_expansion_no_change_var_so_so(functionSo, 
1998
                                                degreeSo, 
1999
                                                rangeSo,
2000
                                                errorTypeSo=None, 
2001
                                                sollyaPrecSo=None):
2002
    """
2003
    Compute the Taylor expansion without the variable change
2004
    x -> x-intervalCenter.
2005
    """
2006
    # Change internal Sollya precision, if needed.
2007
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
2008
    sollyaPrecChanged = False
2009
    if sollyaPrecSo is None:
2010
        pass
2011
    else:
2012
        sollya_lib_set_prec(sollyaPrecSo)
2013
        sollyaPrecChanged = True
2014
    # Error type stuff: default to absolute.
2015
    if errorTypeSo is None:
2016
        errorTypeIsNone = True
2017
        errorTypeSo = sollya_lib_absolute(None)
2018
    else:
2019
        errorTypeIsNone = False
2020
    intervalCenterSo = sollya_lib_mid(rangeSo, None)
2021
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo,
2022
                                         intervalCenterSo,
2023
                                         rangeSo, errorTypeSo, None)
2024
    # taylorFormListSaSo is a Python list of Sollya objects references that 
2025
    # are copies of the elements of taylorFormSo.
2026
    # pobyso_get_list_elements_so_so clears taylorFormSo.
2027
    (taylorFormListSaSo, numElementsSa, isEndEllipticSa) = \
2028
        pobyso_get_list_elements_so_so(taylorFormSo)
2029
    polySo = sollya_lib_copy_obj(taylorFormListSaSo[0])
2030
    #print "Num elements:", numElementsSa
2031
    sollya_lib_clear_obj(taylorFormSo)
2032
    #polySo = taylorFormListSaSo[0]
2033
    #errorRangeSo = sollya_lib_copy_obj(taylorFormListSaSo[2])
2034
    errorRangeSo = taylorFormListSaSo[2]
2035
    # No copy_obj needed here: a new objects are created.
2036
    maxErrorSo    = sollya_lib_sup(errorRangeSo)
2037
    minErrorSo    = sollya_lib_inf(errorRangeSo)
2038
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
2039
    absMinErrorSo = sollya_lib_abs(minErrorSo)
2040
    sollya_lib_clear_obj(maxErrorSo)
2041
    sollya_lib_clear_obj(minErrorSo)
2042
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
2043
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
2044
    # If changed, reset the Sollya working precision.
2045
    if sollyaPrecChanged:
2046
        sollya_lib_set_prec(initialSollyaPrecSo)
2047
    sollya_lib_clear_obj(initialSollyaPrecSo)
2048
    if errorTypeIsNone:
2049
        sollya_lib_clear_obj(errorTypeSo)
2050
    pobyso_clear_taylorform_sa_so(taylorFormListSaSo)
2051
    if absMaxErrorSa > absMinErrorSa:
2052
        sollya_lib_clear_obj(absMinErrorSo)
2053
        return (polySo, intervalCenterSo, absMaxErrorSo)
2054
    else:
2055
        sollya_lib_clear_obj(absMaxErrorSo)
2056
        return (polySo, intervalCenterSo, absMinErrorSo)
2057
# end pobyso_taylor_expansion_no_change_var_so_so
2058

    
2059
def pobyso_taylor_expansion_with_change_var_so_so(functionSo, degreeSo, \
2060
                                                  rangeSo, \
2061
                                                  errorTypeSo=None, \
2062
                                                  sollyaPrecSo=None):
2063
    """
2064
    Compute the Taylor expansion with the variable change
2065
    x -> (x-intervalCenter) included.
2066
    """
2067
    # Change Sollya internal precision, if need.
2068
    sollyaPrecChanged = False
2069
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_sos_sa()
2070
    if sollyaPrecSo is None:
2071
        pass
2072
    else:
2073
        sollya_lib_set_prec(sollyaPrecSo)
2074
        sollyaPrecChanged = True
2075
    #
2076
    # Error type stuff: default to absolute.
2077
    if errorTypeSo is None:
2078
        errorTypeIsNone = True
2079
        errorTypeSo = sollya_lib_absolute(None)
2080
    else:
2081
        errorTypeIsNone = False
2082
    intervalCenterSo = sollya_lib_mid(rangeSo)
2083
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo, \
2084
                                         intervalCenterSo, \
2085
                                         rangeSo, errorTypeSo, None)
2086
    # taylorFormListSaSo is a Python list of Sollya objects references that 
2087
    # are copies of the elements of taylorFormSo.
2088
    # pobyso_get_list_elements_so_so clears taylorFormSo.
2089
    (taylorFormListSo, numElements, isEndElliptic) = \
2090
        pobyso_get_list_elements_so_so(taylorFormSo)
2091
    polySo = taylorFormListSo[0]
2092
    errorRangeSo = taylorFormListSo[2]
2093
    maxErrorSo    = sollya_lib_sup(errorRangeSo)
2094
    minErrorSo    = sollya_lib_inf(errorRangeSo)
2095
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
2096
    absMinErrorSo = sollya_lib_abs(minErrorSo)
2097
    sollya_lib_clear_obj(maxErrorSo)
2098
    sollya_lib_clear_obj(minErrorSo)
2099
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
2100
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
2101
    changeVarExpSo = sollya_lib_build_function_sub(\
2102
                       sollya_lib_build_function_free_variable(),\
2103
                       sollya_lib_copy_obj(intervalCenterSo))
2104
    polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpSo)
2105
    sollya_lib_clear_obj(polySo) 
2106
    sollya_lib_clear_obj(changeVarExpSo)
2107
    # If changed, reset the Sollya working precision.
2108
    if sollyaPrecChanged:
2109
        sollya_lib_set_prec(initialSollyaPrecSo)
2110
    sollya_lib_clear_obj(initialSollyaPrecSo)
2111
    if errorTypeIsNone:
2112
        sollya_lib_clear_obj(errorTypeSo)
2113
    sollya_lib_clear_obj(taylorFormSo)
2114
    # Do not clear maxErrorSo.
2115
    if absMaxErrorSa > absMinErrorSa:
2116
        sollya_lib_clear_obj(absMinErrorSo)
2117
        return((polyVarChangedSo, intervalCenterSo, absMaxErrorSo))
2118
    else:
2119
        sollya_lib_clear_obj(absMaxErrorSo)
2120
        return((polyVarChangedSo, intervalCenterSo, absMinErrorSo))
2121
# end pobyso_taylor_expansion_with_change_var_so_so
2122

    
2123
def pobyso_taylor(function, degree, point):
2124
    """ Legacy function. See pobysoTaylor_so_so. """
2125
    return(pobyso_taylor_so_so(function, degree, point))
2126

    
2127
def pobyso_taylor_so_so(functionSo, degreeSo, pointSo):
2128
    return(sollya_lib_taylor(functionSo, degreeSo, pointSo))
2129
    
2130
def pobyso_taylorform(function, degree, point = None, 
2131
                      interval = None, errorType=None):
2132
    """ Legacy function. See pobyso_taylorform_sa_sa;"""
2133
    
2134
def pobyso_taylorform_sa_sa(functionSa, \
2135
                            degreeSa, \
2136
                            pointSa, \
2137
                            intervalSa=None, \
2138
                            errorTypeSa=None, \
2139
                            precisionSa=None):
2140
    """
2141
    Compute the Taylor form of 'degreeSa' for 'functionSa' at 'pointSa' 
2142
    for 'intervalSa' with 'errorTypeSa' (a string) using 'precisionSa'. 
2143
    point: must be a Real or a Real interval.
2144
    return the Taylor form as an array
2145
    TODO: take care of the interval and of the point when it is an interval;
2146
          when errorType is not None;
2147
          take care of the other elements of the Taylor form (coefficients 
2148
          errors and delta.
2149
    """
2150
    # Absolute as the default error.
2151
    if errorTypeSa is None:
2152
        errorTypeSo = sollya_lib_absolute()
2153
    elif errorTypeSa == "relative":
2154
        errorTypeSo = sollya_lib_relative()
2155
    elif errortypeSa == "absolute":
2156
        errorTypeSo = sollya_lib_absolute()
2157
    else:
2158
        # No clean up needed.
2159
        return None
2160
    # Global precision stuff
2161
    sollyaPrecisionChangedSa = False
2162
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
2163
    if precisionSa is None:
2164
        precSa = initialSollyaPrecSa
2165
    else:
2166
        if precSa > initialSollyaPrecSa:
2167
            if precSa <= 2:
2168
                print inspect.stack()[0][3], ":precision change <= 2 requested."
2169
            pobyso_set_prec_sa_so(precSa)
2170
            sollyaPrecisionChangedSa = True
2171
    #        
2172
    if len(functionSa.variables()) > 0:
2173
        varSa = functionSa.variables()[0]
2174
        pobyso_name_free_variable_sa_so(str(varSa))
2175
    # In any case (point or interval) the parent of pointSa has a precision
2176
    # method.
2177
    pointPrecSa = pointSa.parent().precision()
2178
    if precSa > pointPrecSa:
2179
        pointPrecSa = precSa
2180
    # In any case (point or interval) pointSa has a base_ring() method.
2181
    pointBaseRingString = str(pointSa.base_ring())
2182
    if re.search('Interval', pointBaseRingString) is None: # Point
2183
        pointSo = pobyso_constant_sa_so(pointSa, pointPrecSa)
2184
    else: # Interval.
2185
        pointSo = pobyso_interval_to_range_sa_so(pointSa, pointPrecSa)
2186
    # Sollyafy the function.
2187
    functionSo = pobyso_parse_string_sa_so(functionSa._assume_str().replace('_SAGE_VAR_', ''))
2188
    if sollya_lib_obj_is_error(functionSo):
2189
        print "pobyso_tailorform: function string can't be parsed!"
2190
        return None
2191
    # Sollyafy the degree
2192
    degreeSo = sollya_lib_constant_from_int(int(degreeSa))
2193
    # Sollyafy the point
2194
    # Call Sollya
2195
    taylorFormSo = \
2196
        sollya_lib_taylorform(functionSo, degreeSo, pointSo, errorTypeSo,\
2197
                                         None)
2198
    sollya_lib_clear_obj(functionSo)
2199
    sollya_lib_clear_obj(degreeSo)
2200
    sollya_lib_clear_obj(pointSo)
2201
    sollya_lib_clear_obj(errorTypeSo)
2202
    (tfsAsList, numElements, isEndElliptic) = \
2203
            pobyso_get_list_elements_so_so(taylorFormSo)
2204
    polySo = tfsAsList[0]
2205
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
2206
    polyRealField = RealField(maxPrecision)
2207
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, polyRealField)
2208
    if sollyaPrecisionChangedSa:
2209
        sollya_lib_set_prec(initialSollyaPrecSo)
2210
    sollya_lib_clear_obj(initialSollyaPrecSo)
2211
    polynomialRing = polyRealField[str(varSa)]
2212
    polySa = polynomial(expSa, polynomialRing)
2213
    taylorFormSa = [polySa]
2214
    # Final clean-up.
2215
    sollya_lib_clear_obj(taylorFormSo)
2216
    return(taylorFormSa)
2217
# End pobyso_taylor_form_sa_sa
2218

    
2219
def pobyso_taylorform_so_so(functionSo, degreeSo, pointSo, intervalSo=None, \
2220
                            errorTypeSo=None):
2221
    createdErrorType = False
2222
    if errorTypeSo is None:
2223
        errorTypeSo = sollya_lib_absolute()
2224
        createdErrorType = True
2225
    else:
2226
        #TODO: deal with the other case.
2227
        pass
2228
    if intervalSo is None:
2229
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
2230
                                         errorTypeSo, None)
2231
    else:
2232
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
2233
                                         intervalSo, errorTypeSo, None)
2234
    if createdErrorType:
2235
        sollya_lib_clear_obj(errorTypeSo)
2236
    return resultSo
2237
        
2238

    
2239
def pobyso_univar_polynomial_print_reverse(polySa):
2240
    """ Legacy function. See pobyso_univar_polynomial_print_reverse_sa_sa. """
2241
    return(pobyso_univar_polynomial_print_reverse_sa_sa(polySa))
2242

    
2243
def pobyso_univar_polynomial_print_reverse_sa_sa(polySa):
2244
    """
2245
    Return the string representation of a univariate polynomial with
2246
    monomials ordered in the x^0..x^n order of the monomials.
2247
    Remember: Sage
2248
    """
2249
    polynomialRing = polySa.base_ring()
2250
    # A very expensive solution:
2251
    # -create a fake multivariate polynomial field with only one variable,
2252
    #   specifying a negative lexicographical order;
2253
    mpolynomialRing = PolynomialRing(polynomialRing.base(), \
2254
                                     polynomialRing.variable_name(), \
2255
                                     1, order='neglex')
2256
    # - convert the univariate argument polynomial into a multivariate
2257
    #   version;
2258
    p = mpolynomialRing(polySa)
2259
    # - return the string representation of the converted form.
2260
    # There is no simple str() method defined for p's class.
2261
    return(p.__str__())
2262
#
2263
#print pobyso_get_prec()  
2264
pobyso_set_prec(165)
2265
#print pobyso_get_prec()  
2266
#a=100
2267
#print type(a)
2268
#id(a)
2269
#print "Max arity: ", pobyso_max_arity
2270
#print "Function tripleDouble (43) as a string: ", pobyso_function_type_as_string(43)
2271
#print "Function None (44) as a string: ", pobyso_function_type_as_string(44)
2272
sys.stderr.write("\t...Pobyso check done.\n")