Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 257

Historique | Voir | Annoter | Télécharger (86,16 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
NOTES:
28
Reported errors in Eclipse come from the calls to
29
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
    return(sollya_lib_absolute(None))
93

    
94
def pobyso_autoprint(arg):
95
    sollya_lib_autoprint(arg, None)
96

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

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

    
157
def pobyso_build_function_sub_so_so(exp1So, exp2So):
158
    return sollya_lib_build_function_sub(exp1So, exp2So)
159

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

    
195
def pobyso_change_var_in_function_so_so(funcSo, chvarExpSo):
196
    """
197
    Variable change in a function.
198
    """
199
    return(sollya_lib_evaluate(funcSo,chvarExpSo))
200
# End pobyso_change_var_in_function_so_so     
201

    
202
def pobyso_chebyshevform_so_so(functionSo, degreeSo, intervalSo):
203
    resultSo = sollya_lib_chebyshevform(functionSo, degreeSo, intervalSo)
204
    return(resultSo)
205
# End pobyso_chebyshevform_so_so.
206

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

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

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

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

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

    
316

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

    
345
def pobyso_constant_1():
346
    """
347
    Obvious.
348
    Legacy function. See pobyso_constant_so_so. 
349
    """
350
    return pobyso_constant_1_sa_so()
351

    
352
def pobyso_constant_1_sa_so():
353
    """
354
    Obvious.
355
    """
356
    return(pobyso_constant_from_int_sa_so(1))
357

    
358
def pobyso_constant_from_int(anInt):
359
    """ Legacy function. See pobyso_constant_from_int_sa_so. """
360
    return pobyso_constant_from_int_sa_so(anInt)
361

    
362
def pobyso_constant_from_int_sa_so(anInt):
363
    """
364
    Get a Sollya constant from a Sage int.
365
    """
366
    return sollya_lib_constant_from_int64(long(anInt))
367

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

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

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

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

    
408
def pobyso_error_so():
409
    return sollya_lib_error(None)
410
# End pobyso_error().
411

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

    
426
def pobyso_dirty_find_zeros_so_so(funcSo, rangeSo):
427
    """
428
    Thin wrapper over sollya_lib_dirtyfindzeros()
429
    """
430
    return sollya_lib_dirtyfindzeros(funcSo, rangeSo)
431
# End pobys_dirty_find_zeros
432

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

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

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

    
527
def pobyso_float_poly_so_sa(polySo, realFieldSa=None):
528
    """
529
    Convert a Sollya polynomial into a Sage floating-point polynomial.
530
    If no realField is given, a RealField corresponding to the maximum 
531
    precision of the coefficients is internally computed.
532
    The real field is not returned but can be easily retrieved from 
533
    the polynomial itself.
534
    ALGORITHM:
535
    - (optional) compute the RealField of the coefficients;
536
    - convert the Sollya expression into a Sage expression;
537
    - convert the Sage expression into a Sage polynomial
538
    """    
539
    if realFieldSa is None:
540
        expressionPrecSa = pobyso_get_max_prec_of_exp_so_sa(polySo)
541
        #print "Maximum precision of Sollya polynomial coefficients:", expressionPrecSa
542
        if expressionPrecSa < 2 or expressionPrecSa > 2147483391:
543
            print "Maximum degree of expression:", expressionPrecSa
544
        realFieldSa      = RealField(expressionPrecSa)
545
    #print "Sollya expression before...",
546
    #pobyso_autoprint(polySo)
547

    
548
    expressionSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo,
549
                                                             realFieldSa)
550
    #print "...Sollya expression after."
551
    #pobyso_autoprint(polySo)
552
    polyVariableSa = expressionSa.variables()[0]
553
    polyRingSa     = realFieldSa[str(polyVariableSa)]
554
    #print polyRingSa
555
    # Do not use the polynomial(expressionSa, ring=polyRingSa) form!
556
    polynomialSa = polyRingSa(expressionSa)
557
    polyCoeffsListSa = polynomialSa.coefficients()
558
    #for coeff in polyCoeffsListSa:
559
    #    print coeff.abs().n()
560
    return polynomialSa
561
# End pobyso_float_poly_so_sa
562

    
563
def pobyso_free_variable():
564
    """
565
    Ultra thin wrapper around the sollya_lib_function_build_free_variable function.
566
    """
567
    return sollya_lib_build_function_free_variable()
568
    
569
def pobyso_function_type_as_string(funcType):
570
    """ Legacy function. See pobyso_function_type_as_string_so_sa. """
571
    return(pobyso_function_type_as_string_so_sa(funcType))
572

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

    
669
def pobyso_get_constant(rnArgSa, constSo):
670
    """ Legacy function. See pobyso_get_constant_so_sa. """
671
    return pobyso_get_constant_so_sa(rnArgSa, constSo)
672
# End pobyso_get_constant
673

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

    
716
def pobyso_get_constant_as_rn_with_rf(ctExp, realField):
717
    """ 
718
    Legacy function. See pobyso_get_constant_as_rn_with_rf_so_sa.
719
    """
720
    return pobyso_get_constant_as_rn_with_rf_so_sa(ctExp, realField)
721
# End pobyso_get_constant_as_rn_with_rf
722
    
723
def pobyso_get_constant_as_rn_with_rf_so_sa(ctExpSo, realFieldSa = None):
724
    """
725
    Get a Sollya constant as a Sage "real number".
726
    If no real field is specified, the precision of the floating-point number 
727
    returned is that of the Sollya constant.
728
    Otherwise is is that of the real field. Hence rounding may happen.
729
    """
730
    if realFieldSa is None:
731
        return pobyso_get_constant_as_rn_so_sa(ctExpSo)
732
    rnSa = realFieldSa(0)
733
    outcome = sollya_lib_get_constant(get_rn_value(rnSa), ctExpSo)
734
    if outcome == 0:
735
        return None
736
    else:
737
        return rnSa
738
# End pobyso_get_constant_as_rn_with_rf_so_sa
739

    
740
def pobyso_get_free_variable_name():
741
    """ 
742
    Legacy function. See pobyso_get_free_variable_name_so_sa.
743
    """
744
    return(pobyso_get_free_variable_name_so_sa())
745

    
746
def pobyso_get_free_variable_name_so_sa():
747
    return sollya_lib_get_free_variable_name()
748
    
749
def pobyso_get_function_arity(expressionSo):
750
    """ 
751
    Legacy function. See pobyso_get_function_arity_so_sa.
752
    """
753
    return(pobyso_get_function_arity_so_sa(expressionSo))
754

    
755
def pobyso_get_function_arity_so_sa(expressionSo):
756
    arity = c_int(0)
757
    sollya_lib_get_function_arity(byref(arity),expressionSo)
758
    return int(arity.value)
759

    
760
def pobyso_get_head_function(expressionSo):
761
    """ 
762
    Legacy function. See pobyso_get_head_function_so_sa. 
763
    """
764
    return(pobyso_get_head_function_so_sa(expressionSo)) 
765

    
766
def pobyso_get_head_function_so_sa(expressionSo):
767
    functionType = c_int(0)
768
    sollya_lib_get_head_function(byref(functionType), expressionSo)
769
    return int(functionType.value)
770

    
771
def pobyso_get_interval_from_range_so_sa(soRange, realIntervalFieldSa = None ):
772
    """
773
    Return the Sage interval corresponding to the Sollya range argument.
774
    If no reaIntervalField is passed as an argument, the interval bounds are not
775
    rounded: they are elements of RealIntervalField of the "right" precision
776
    to hold all the digits.
777
    """
778
    prec = c_int(0)
779
    if realIntervalFieldSa is None:
780
        retval = sollya_lib_get_prec_of_range(byref(prec), soRange, None)
781
        if retval == 0:
782
            return None
783
        realIntervalFieldSa = RealIntervalField(prec.value)
784
    intervalSa = realIntervalFieldSa(0,0)
785
    retval = \
786
        sollya_lib_get_interval_from_range(get_interval_value(intervalSa),\
787
                                           soRange)
788
    if retval == 0:
789
        return None
790
    return intervalSa
791
# End pobyso_get_interval_from_range_so_sa
792

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

    
836
def pobyso_get_max_prec_of_exp(soExp):
837
    """ Legacy function. See pobyso_get_max_prec_of_exp_so_sa. """
838
    return pobyso_get_max_prec_of_exp_so_sa(soExp)
839

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

    
881
def pobyso_get_min_prec_of_constant_so_sa(constExpSo):
882
    """
883
    Get the minimum precision necessary to represent the value of a Sollya
884
    constant.
885
    MPFR_MIN_PREC and powers of 2 are taken into account.
886
    We assume that constExpSo is a pointer to a Sollay constant expression.
887
    """
888
    constExpAsRnSa = pobyso_get_constant_as_rn_so_sa(constExpSo)
889
    return(min_mpfr_size(get_rn_value(constExpAsRnSa)))
890

    
891
def pobyso_get_poly_so_sa(polySo, realFieldSa=None):
892
    """
893
    Convert a Sollya polynomial into a Sage polynomial.
894
    Legacy function. Use pobyso_float_poly_so_sa() instead.
895
    """    
896
    return pobyso_float_poly_so_sa(polySo,realFieldSa)
897
# End pobyso_get_poly_so_sa
898

    
899
def pobyso_get_prec():
900
    """ Legacy function. See pobyso_get_prec_so_sa(). """
901
    return pobyso_get_prec_so_sa()
902

    
903
def pobyso_get_prec_so():
904
    """
905
    Get the current default precision in Sollya.
906
    The return value is a Sollya object.
907
    Usefull when modifying the precision back and forth by avoiding
908
    extra conversions.
909
    """
910
    return sollya_lib_get_prec(None)
911
    
912
def pobyso_get_prec_so_sa():
913
    """
914
    Get the current default precision in Sollya.
915
    The return value is Sage/Python int.
916
    """
917
    precSo = sollya_lib_get_prec()
918
    precSa = pobyso_constant_from_int_so_sa(precSo)
919
    sollya_lib_clear_obj(precSo)
920
    return precSa
921
# End pobyso_get_prec_so_sa.
922

    
923
def pobyso_get_prec_so_so_sa():
924
    """
925
    Return the current precision both as a Sollya object and a
926
    Sage integer as hybrid tuple.
927
    To avoid multiple calls for precision manipulations. 
928
    """
929
    precSo = sollya_lib_get_prec()
930
    precSa = pobyso_constant_from_int_so_sa(precSo)
931
    return (precSo, int(precSa))
932
    
933
def pobyso_get_prec_of_constant(ctExpSo):
934
    """ Legacy function. See pobyso_get_prec_of_constant_so_sa. """
935
    return pobyso_get_prec_of_constant_so_sa(ctExpSo)
936

    
937
def pobyso_get_prec_of_constant_so_sa(ctExpSo):
938
    """
939
    Tries to find a precision to represent ctExpSo without rounding.
940
    If not possible, returns None.
941
    """
942
    #print "Entering pobyso_get_prec_of_constant_so_sa..."
943
    prec = c_int(0)
944
    retc = sollya_lib_get_prec_of_constant(byref(prec), ctExpSo, None)
945
    if retc == 0:
946
        #print "pobyso_get_prec_of_constant_so_sa failed."
947
        return None
948
    #print "...exiting pobyso_get_prec_of_constant_so_sa."
949
    return int(prec.value)
950

    
951
def pobyso_get_prec_of_range_so_sa(rangeSo):
952
    """
953
    Returns the number of bits elements of a range are coded with.
954
    """
955
    prec = c_int(0)
956
    retc = sollya_lib_get_prec_of_range(byref(prec), rangeSo, None)
957
    if retc == 0:
958
        return(None)
959
    return int(prec.value)
960
# End pobyso_get_prec_of_range_so_sa()
961

    
962
def pobyso_get_sage_exp_from_sollya_exp(sollyaExpSo, realField = RR):
963
    """ Legacy function. See pobyso_get_sage_exp_from_sollya_exp_so_sa. """
964
    return pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, 
965
                                                     realField = RR)
966

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

    
1015

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

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

    
1165
def pobyso_inf_so_so(intervalSo):
1166
    """
1167
    Very thin wrapper around sollya_lib_inf().
1168
    """
1169
    return sollya_lib_inf(intervalSo)
1170
# End pobyso_inf_so_so.
1171
    
1172
def pobyso_infnorm_so_so(func, interval, file = None, intervalList = None):
1173
    print "Do not use this function. User pobyso_supnorm_so_so instead."
1174
    return None
1175

    
1176
def pobyso_interval_to_range_sa_so(intervalSa, precisionSa=None):
1177
    if precisionSa is None:
1178
        precisionSa = intervalSa.parent().precision()
1179
    intervalSo = pobyso_bounds_to_range_sa_so(intervalSa.lower(),\
1180
                                              intervalSa.upper(),\
1181
                                              precisionSa)
1182
    return intervalSo
1183
# End pobyso_interval_to_range_sa_so
1184

    
1185
def pobyso_is_error_so_sa(objSo):
1186
    """
1187
    Thin wrapper around the sollya_lib_obj_is_error() function.
1188
    """
1189
    if sollya_lib_obj_is_error(objSo) != 0:
1190
        return True
1191
    else:
1192
        return False
1193
# End pobyso_is_error-so_sa
1194

    
1195
def pobyso_is_floating_point_number_sa_sa(numberSa):
1196
    """
1197
    Check whether a Sage number is floating point.
1198
    Exception stuff added because numbers other than
1199
    floating-point ones do not have the is_real() attribute.
1200
    """
1201
    try:
1202
        return numberSa.is_real()
1203
    except AttributeError:
1204
        return False
1205
# End pobyso_is_floating_piont_number_sa_sa
1206

    
1207
def pobyso_is_range_so_sa(rangeCandidateSo):
1208
    """
1209
    Thin wrapper over sollya_lib_is_range.
1210
    """
1211
    return sollya_lib_obj_is_range(rangeCandidateSo) != 0
1212
    
1213
# End pobyso_is_range_so_sa
1214

    
1215

    
1216
def pobyso_lib_init():
1217
    sollya_lib_init(None)
1218

    
1219
def pobyso_lib_close():
1220
    sollya_lib_close(None)
1221

    
1222
def pobyso_name_free_variable(freeVariableNameSa):
1223
    """ Legacy function. See pobyso_name_free_variable_sa_so. """
1224
    pobyso_name_free_variable_sa_so(freeVariableNameSa)
1225

    
1226
def pobyso_name_free_variable_sa_so(freeVariableNameSa):
1227
    """
1228
    Set the free variable name in Sollya from a Sage string.
1229
    """
1230
    sollya_lib_name_free_variable(freeVariableNameSa)
1231

    
1232
def pobyso_parse_string(string):
1233
    """ Legacy function. See pobyso_parse_string_sa_so. """
1234
    return pobyso_parse_string_sa_so(string)
1235
 
1236
def pobyso_parse_string_sa_so(string):
1237
    """
1238
    Get the Sollya expression computed from a Sage string or
1239
    a Sollya error object if parsing failed.
1240
    """
1241
    return sollya_lib_parse_string(string)
1242

    
1243
def pobyso_precision_so_sa(ctExpSo):
1244
    """
1245
    Computes the necessary precision to represent a number.
1246
    If x is not zero, it can be uniquely written as x = m · 2e 
1247
    where m is an odd integer and e is an integer. 
1248
    precision(x) returns the number of bits necessary to write m 
1249
    in binary (i.e. ceil(log2(m))).
1250
    """
1251
    #TODO: take care of the special case: 0, @NaN@, @Inf@
1252
    precisionSo = sollya_lib_precision(ctExpSo)
1253
    precisionSa = pobyso_constant_from_int_so_sa(precisionSo)
1254
    sollya_lib_clear_obj(precisionSo)
1255
    return precisionSa
1256
# End pobyso_precision_so_sa
1257

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

    
1424
def pobyso_polynomial_degree_so_sa(polySo):
1425
    """
1426
    Return the degree of a Sollya polynomial as a Sage int.
1427
    """
1428
    degreeSo = sollya_lib_degree(polySo)
1429
    return pobyso_constant_from_int_so_sa(degreeSo)
1430
# End pobyso_polynomial_degree_so_sa
1431

    
1432
def pobyso_polynomial_degree_so_so(polySo):
1433
    """
1434
    Thin wrapper around lib_sollya_degree().
1435
    """
1436
    return sollya_lib_degree(polySo)
1437
# End pobyso_polynomial_degree_so_so
1438
                                                                  
1439
def pobyso_range(rnLowerBound, rnUpperBound):
1440
    """ Legacy function. See pobyso_range_sa_so. """
1441
    return pobyso_range_sa_so(rnLowerBound, rnUpperBound) 
1442

    
1443
def pobyso_range_from_bounds_sa_so(rnLowerBound, rnUpperBound, precSa = None):
1444
    """
1445
    Create a Sollya range from 2 Sage real numbers as bounds
1446
    """
1447
    # TODO check precision stuff.
1448
    sollyaPrecChanged = False
1449
    (initialSollyaPrecSo, initialSollyaPrecSa) = \
1450
        pobyso_get_prec_so_so_sa()
1451
    if precSa is None:
1452
        precSa = max(rnLowerBound.parent().prec(), rnUpperBound.parent().prec())
1453
    if precSa > initialSollyaPrecSa:
1454
        if precSa <= 2:
1455
            print inspect.stack()[0][3], ": precision change <= 2 requested."
1456
        pobyso_set_prec_sa_so(precSa)
1457
        sollyaPrecChanged = True
1458
    rangeSo = sollya_lib_range_from_bounds(get_rn_value(rnLowerBound),
1459
                                           get_rn_value(rnUpperBound))
1460
    if sollyaPrecChanged:
1461
        pobyso_set_prec_so_so(initialSollyaPrecSo)
1462
        sollya_lib_clear_obj(initialSollyaPrecSo)
1463
    return rangeSo
1464
# End pobyso_range_from_bounds_sa_so
1465

    
1466
def pobyso_range_max_abs_so_so(rangeSo):
1467
    """
1468
    Return, as Sollya constant, the maximum absolute value of a Sollay range.
1469
    """
1470
    lowerBoundSo = sollya_lib_inf(rangeSo)
1471
    upperBoundSo = sollya_lib_sup(rangeSo)
1472
    #
1473
    lowerBoundSo = sollya_lib_build_function_abs(lowerBoundSo)
1474
    upperBoundSo = sollya_lib_build_function_abs(upperBoundSo)
1475
    #pobyso_autoprint(lowerBoundSo)
1476
    #pobyso_autoprint(upperBoundSo)
1477
    #
1478
    maxAbsSo = sollya_lib_max(lowerBoundSo, upperBoundSo, None)
1479
    #sollya_lib_clear_obj(lowerBoundSo)
1480
    #sollya_lib_clear_obj(upperBoundSo)
1481
    return maxAbsSo
1482
# End pobyso_range_max_abs_so_so
1483
 
1484
def pobyso_range_to_interval_so_sa(rangeSo, realIntervalFieldSa = None):
1485
    """
1486
    Get a Sage interval from a Sollya range.
1487
    If no realIntervalField is given as a parameter, the Sage interval
1488
    precision is that of the Sollya range.
1489
    Otherwise, the precision is that of the realIntervalField. In this case
1490
    rounding may happen.
1491
    """
1492
    if realIntervalFieldSa is None:
1493
        precSa = pobyso_get_prec_of_range_so_sa(rangeSo)
1494
        realIntervalFieldSa = RealIntervalField(precSa)
1495
    intervalSa = \
1496
        pobyso_get_interval_from_range_so_sa(rangeSo, realIntervalFieldSa) 
1497
    return intervalSa
1498
# End pobyso_range_to_interval_so_sa
1499
#
1500
def pobyso_relative_so_so():
1501
    """
1502
    Very thin wrapper around the sollya_lib_relative function.
1503
    """
1504
    return sollya_lib_relative()
1505
# End pobyso_relative_so_so
1506
#
1507
def pobyso_rat_poly_sa_so(polySa, precSa = None):
1508
    """
1509
    Create a Sollya polynomial from a Sage rational polynomial.
1510
    """
1511
    ## TODO: filter arguments.
1512
    ## Precision. If no precision is given, use the current precision
1513
    #  of Sollya.
1514
    if precSa is None:
1515
        precSa =  pobyso_get_prec_so_sa()
1516
    #print "Precision:",  precSa
1517
    RRR = RealField(precSa)
1518
    ## Create a Sage polynomial in the "right" precision.
1519
    P_RRR = RRR[polySa.variables()[0]]
1520
    polyFloatSa = P_RRR(polySa)
1521
    ## Make sure no precision is provided: pobyso_float_poly_sa_so will
1522
    #  recover it all by itself and not make an extra conversion.
1523
    return pobyso_float_poly_sa_so(polyFloatSa)
1524
    
1525
# End pobyso_rat_poly_sa_so    
1526
    
1527
def pobyso_remez_canonical_sa_sa(func, \
1528
                                 degree, \
1529
                                 lowerBound, \
1530
                                 upperBound, \
1531
                                 weight = None, \
1532
                                 quality = None):
1533
    """
1534
    All arguments are Sage/Python.
1535
    The functions (func and weight) must be passed as expressions or strings.
1536
    Otherwise the function fails. 
1537
    The return value is a Sage polynomial.
1538
    """
1539
    var('zorglub')    # Dummy variable name for type check only. Type of 
1540
    # zorglub is "symbolic expression".
1541
    polySo = pobyso_remez_canonical_sa_so(func, \
1542
                                 degree, \
1543
                                 lowerBound, \
1544
                                 upperBound, \
1545
                                 weight, \
1546
                                 quality)
1547
    # String test
1548
    if parent(func) == parent("string"):
1549
        functionSa = eval(func)
1550
    # Expression test.
1551
    elif type(func) == type(zorglub):
1552
        functionSa = func
1553
    else:
1554
        return None
1555
    #
1556
    maxPrecision = 0
1557
    if polySo is None:
1558
        return(None)
1559
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
1560
    RRRRSa = RealField(maxPrecision)
1561
    polynomialRingSa = RRRRSa[functionSa.variables()[0]]
1562
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, RRRRSa)
1563
    polySa = polynomial(expSa, polynomialRingSa)
1564
    sollya_lib_clear_obj(polySo)
1565
    return(polySa)
1566
# End pobyso_remez_canonical_sa_sa
1567
    
1568
def pobyso_remez_canonical(func, \
1569
                           degree, \
1570
                           lowerBound, \
1571
                           upperBound, \
1572
                           weight = "1", \
1573
                           quality = None):
1574
    """ Legacy function. See pobyso_remez_canonical_sa_so. """
1575
    return(pobyso_remez_canonical_sa_so(func, \
1576
                                        degree, \
1577
                                        lowerBound, \
1578
                                        upperBound, \
1579
                                        weight, \
1580
                                        quality))
1581
# End pobyso_remez_canonical.
1582

    
1583
def pobyso_remez_canonical_sa_so(func, \
1584
                                 degree, \
1585
                                 lowerBound, \
1586
                                 upperBound, \
1587
                                 weight = None, \
1588
                                 quality = None):
1589
    """
1590
    All arguments are Sage/Python.
1591
    The functions (func and weight) must be passed as expressions or strings.
1592
    Otherwise the function fails. 
1593
    The return value is a pointer to a Sollya function.
1594
    lowerBound and upperBound mus be reals.
1595
    """
1596
    var('zorglub')    # Dummy variable name for type check only. Type of
1597
    # zorglub is "symbolic expression".
1598
    currentVariableNameSa = None
1599
    # The func argument can be of different types (string, 
1600
    # symbolic expression...)
1601
    if parent(func) == parent("string"):
1602
        localFuncSa = sage_eval(func,globals())
1603
        if len(localFuncSa.variables()) > 0:
1604
            currentVariableNameSa = localFuncSa.variables()[0]
1605
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1606
            functionSo = \
1607
              sollya_lib_parse_string(localFuncSa._assume_str().replace('_SAGE_VAR_', ''))
1608
    # Expression test.
1609
    elif type(func) == type(zorglub):
1610
        # Until we are able to translate Sage expressions into Sollya 
1611
        # expressions : parse the string version.
1612
        if len(func.variables()) > 0:
1613
            currentVariableNameSa = func.variables()[0]
1614
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1615
            functionSo = \
1616
              sollya_lib_parse_string(func._assume_str().replace('_SAGE_VAR_', ''))
1617
    else:
1618
        return(None)
1619
    if weight is None: # No weight given -> 1.
1620
        weightSo = pobyso_constant_1_sa_so()
1621
    elif parent(weight) == parent("string"): # Weight given as string: parse it.
1622
        weightSo = sollya_lib_parse_string(func)
1623
    elif type(weight) == type(zorglub): # Weight given as symbolice expression.
1624
        functionSo = \
1625
          sollya_lib_parse_string_sa_so(weight._assume_str().replace('_SAGE_VAR_', ''))
1626
    else:
1627
        return(None)
1628
    degreeSo = pobyso_constant_from_int(degree)
1629
    rangeSo = pobyso_bounds_to_range_sa_so(lowerBound, upperBound)
1630
    if not quality is None:
1631
        qualitySo= pobyso_constant_sa_so(quality)
1632
    else:
1633
        qualitySo = None
1634
        
1635
    remezPolySo = sollya_lib_remez(functionSo, \
1636
                                   degreeSo, \
1637
                                   rangeSo, \
1638
                                   weightSo, \
1639
                                   qualitySo, \
1640
                                   None)
1641
    sollya_lib_clear_obj(functionSo)
1642
    sollya_lib_clear_obj(degreeSo)
1643
    sollya_lib_clear_obj(rangeSo)
1644
    sollya_lib_clear_obj(weightSo)
1645
    if not qualitySo is None:
1646
        sollya_lib_clear_obj(qualitySo)
1647
    return(remezPolySo)
1648
# End pobyso_remez_canonical_sa_so
1649

    
1650
def pobyso_remez_canonical_so_so(funcSo, \
1651
                                 degreeSo, \
1652
                                 rangeSo, \
1653
                                 weightSo = pobyso_constant_1_sa_so(),\
1654
                                 qualitySo = None):
1655
    """
1656
    All arguments are pointers to Sollya objects.
1657
    The return value is a pointer to a Sollya function.
1658
    """
1659
    if not sollya_lib_obj_is_function(funcSo):
1660
        return(None)
1661
    return(sollya_lib_remez(funcSo, degreeSo, rangeSo, weightSo, qualitySo, None))
1662
# End pobyso_remez_canonical_so_so.
1663

    
1664
def pobyso_remez_exponents_list_sa_so(func,     \
1665
                                 exponentsList, \
1666
                                 lowerBound,    \
1667
                                 upperBound,    \
1668
                                 weight = None, \
1669
                                 quality = None):
1670
    """
1671
    All arguments are Sage/Python.
1672
    The functions (func and weight) must be passed as expressions or strings.
1673
    Otherwise the function fails. 
1674
    The return value is a pointer to a Sollya function.
1675
    lowerBound and upperBound mus be reals.
1676
    """
1677
    var('zorglub')    # Dummy variable name for type check only. Type of
1678
    # zorglub is "symbolic expression".
1679
    currentVariableNameSa = None
1680
    # The func argument can be of different types (string, 
1681
    # symbolic expression...)
1682
    if parent(func) == parent("string"):
1683
        localFuncSa = sage_eval(func,globals())
1684
        if len(localFuncSa.variables()) > 0:
1685
            currentVariableNameSa = localFuncSa.variables()[0]
1686
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1687
            functionSo = \
1688
              sollya_lib_parse_string(localFuncSa._assume_str().replace('_SAGE_VAR_', ''))
1689
    # Expression test.
1690
    elif type(func) == type(zorglub):
1691
        # Until we are able to translate Sage expressions into Sollya 
1692
        # expressions : parse the string version.
1693
        if len(func.variables()) > 0:
1694
            currentVariableNameSa = func.variables()[0]
1695
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1696
            functionSo = \
1697
              sollya_lib_parse_string(func._assume_str().replace('_SAGE_VAR_', ''))
1698
    else:
1699
        return(None)
1700
    ## Deal with the weight, much in the same way as with the function.
1701
    if weight is None: # No weight given -> 1.
1702
        weightSo = pobyso_constant_1_sa_so()
1703
    elif parent(weight) == parent("string"): # Weight given as string: parse it.
1704
        weightSo = sollya_lib_parse_string(func)
1705
    elif type(weight) == type(zorglub): # Weight given as symbolice expression.
1706
        functionSo = \
1707
          sollya_lib_parse_string(weight._assume_str().replace('_SAGE_VAR_', '',100))
1708
    else:
1709
        return(None)
1710
    rangeSo = pobyso_bounds_to_range_sa_so(lowerBound, upperBound)
1711
    if not quality is None:
1712
        qualitySo= pobyso_constant_sa_so(quality)
1713
    else:
1714
        qualitySo = None
1715
    #
1716
    ## Tranform the Sage list of exponents into a Sollya list.
1717
    exponentsListSo = pobyso_build_list_of_ints_sa_so(*exponentsList)
1718
    remezPolySo = sollya_lib_remez(functionSo, \
1719
                                   exponentsListSo, \
1720
                                   rangeSo, \
1721
                                   weightSo, \
1722
                                   qualitySo, \
1723
                                   None)
1724
    sollya_lib_clear_obj(functionSo)
1725
    sollya_lib_clear_obj(exponentsListSo)
1726
    sollya_lib_clear_obj(rangeSo)
1727
    sollya_lib_clear_obj(weightSo)
1728
    if not qualitySo is None:
1729
        sollya_lib_clear_obj(qualitySo)
1730
    return(remezPolySo)
1731
# End pobyso_remez_exponentsList_sa_so
1732
#
1733
def pobyso_round_coefficients_so_so(polySo, truncFormatListSo):
1734
    """
1735
    A wrapper around the "classical" sollya_lib_roundcoefficients: a Sollya
1736
    polynomial and a Sollya list are given as arguments.
1737
    """
1738
    return sollya_lib_roundcoefficients(polySo, truncFormatListSo)
1739

    
1740
def pobyso_round_coefficients_progressive_so_so(polySo, 
1741
                                                funcSo,
1742
                                                precSo,
1743
                                                intervalSo,
1744
                                                icSo,
1745
                                                currentApproxErrorSo,
1746
                                                approxAccurSo,
1747
                                                debug=False):
1748
    if debug:
1749
        print "Input arguments:"
1750
        print "Polynomial: ", ; pobyso_autoprint(polySo)
1751
        print "Function: ", ; pobyso_autoprint(funcSo)
1752
        print "Internal precision: ", ; pobyso_autoprint(precSo)
1753
        print "Interval: ", ; pobyso_autoprint(intervalSo)
1754
        print "Current approximation error: ", ; pobyso_autoprint(currentApproxErrorSo)
1755
        print "Requested approxiation error: ", ; pobyso_autoprint(approxAccurSo)
1756
        print "________________"
1757
    approxAccurSa        = pobyso_get_constant_as_rn_so_sa(approxAccurSo)
1758
    currentApproxErrorSa = pobyso_get_constant_as_rn_so_sa(currentApproxErrorSo)
1759
    ## If the current approximation error is too close to the target, there is
1760
    #  no possible gain.
1761
    if currentApproxErrorSa >= approxAccurSa / 2:
1762
        #### Do not return the initial argument but copies: caller may free 
1763
        #    as inutile after call.
1764
        return (sollya_lib_copy_obj(polySo),
1765
                sollya_lib_copy_obj(currentApproxErrorSo))
1766
    degreeSa             = pobyso_polynomial_degree_so_sa(polySo)
1767
    intervalSa           = pobyso_range_to_interval_so_sa(intervalSo)
1768
    
1769
    if debug:
1770
        print "degreeSa              :", degreeSa
1771
        print "intervalSa            :", intervalSa.str(style='brackets')
1772
        print "currentApproxErrorSa  :", currentApproxErrorSa 
1773
        print "approxAccurSa         :", approxAccurSa 
1774
    ### Start with a 0 value expression.
1775
    radiusSa = intervalSa.absolute_diameter() / 2
1776
    if debug:
1777
        print "log2(radius):", RR(radiusSa).log2()
1778
    iterIndex = 0
1779
    while True: 
1780
        resPolySo = pobyso_constant_0_sa_so()
1781
        roundedPolyApproxAccurSa = approxAccurSa / 2
1782
        currentRadiusPowerSa = 1 
1783
        for degree in xrange(0,degreeSa + 1):
1784
            #### At round 0, use the agressive formula. At round 1, run the
1785
            #    proved formula.
1786
            if iterIndex == 0:
1787
                roundingPowerSa = \
1788
                    floor(((currentRadiusPowerSa/roundedPolyApproxAccurSa)*(degree+1)).log2())
1789
            else:
1790
                roundingPowerSa = \
1791
                    floor(((currentRadiusPowerSa/roundedPolyApproxAccurSa)*(degreeSa+1)).log2())
1792
            ## Under extreme conditions the above formulas can evaluate under 2, which is the
1793
            #  minimal precision of an MPFR number.
1794
            if roundingPowerSa < 2:
1795
                roundingPowerSa = 2
1796
            if debug:
1797
                print "roundedPolyApproxAccurSa", roundedPolyApproxAccurSa
1798
                print "currentRadiusPowerSa", currentRadiusPowerSa
1799
                print "Current rounding exponent:", roundingPowerSa
1800
            currentRadiusPowerSa *= radiusSa
1801
            index1So = pobyso_constant_from_int_sa_so(degree)
1802
            index2So = pobyso_constant_from_int_sa_so(degree)
1803
            ### Create a monomial with:
1804
            #   - the coefficient in the initial monomial at the current degrree;
1805
            #   - the current exponent;
1806
            #   - the free variable.
1807
            cmonSo  = \
1808
                sollya_lib_build_function_mul(sollya_lib_coeff(polySo, index1So),
1809
                                              sollya_lib_build_function_pow( \
1810
                                                sollya_lib_build_function_free_variable(), \
1811
                                                index2So))
1812
            roundingPowerSo = pobyso_constant_from_int_sa_so(roundingPowerSa)
1813
            cmonrSo = pobyso_round_coefficients_single_so_so(cmonSo, roundingPowerSo)
1814
            sollya_lib_clear_obj(cmonSo)
1815
            ### Add to the result polynomial.
1816
            resPolySo = sollya_lib_build_function_add(resPolySo,
1817
                                                      cmonrSo)
1818
        # End for.
1819
        ### Check the new polynomial.
1820
        freeVarSo     = sollya_lib_build_function_free_variable()
1821
        changeVarSo   = sollya_lib_sub(freeVarSo, icSo)
1822
        resPolyCvSo   = sollya_lib_evaluate(resPolySo, changeVarSo)
1823
        errFuncSo     = sollya_lib_build_function_sub(sollya_lib_copy_obj(funcSo), 
1824
                                                      resPolyCvSo)
1825
        infNormSo = sollya_lib_dirtyinfnorm(errFuncSo, intervalSo)
1826
        ### This also clears resPolyCvSo.
1827
        sollya_lib_clear_obj(errFuncSo)
1828
        cerrSa    = pobyso_get_constant_as_rn_so_sa(infNormSo)
1829
        if debug:
1830
            print "Error of the new polynomial:", cerrSa
1831
        ### If at round 1, return the initial polynomial error. This should
1832
        #   never happen since the rounding algorithm is proved. But some 
1833
        #   circumstances may break it (e.g. internal precision of tools).
1834
        if cerrSa > approxAccurSa:
1835
            if iterIndex > 0: # Round 1 and beyond.
1836
                sollya_lib_clear_obj(resPolySo)
1837
                sollya_lib_clear_obj(infNormSo)
1838
                #### Do not return the arguments but copies: the caller may free
1839
                #    free the former as inutile after call.
1840
                return (sollya_lib_copy_obj(polySo), 
1841
                        sollya_lib_copy_obj(currentApproxErrorSo))
1842
            else: # Round 0, got round 1
1843
                sollya_lib_clear_obj(resPolySo)
1844
                sollya_lib_clear_obj(infNormSo)
1845
                iterIndex += 1
1846
                continue
1847
        ### If get here it is because cerrSa <= approxAccurSa
1848
        ### Approximation error of the new polynomial is acceptable.
1849
        return (resPolySo, infNormSo)
1850
    # End while True
1851
# End pobyso_round_coefficients_progressive_so_so
1852
    
1853
def pobyso_round_coefficients_single_so_so(polySo, commonPrecSo):
1854
    """
1855
    Create a rounded coefficients polynomial from polynomial argument to
1856
    the number of bits in size argument.
1857
    All coefficients are set to the same precision.
1858
    """
1859
    ## TODO: check arguments.
1860
    endEllipListSo = pobyso_build_end_elliptic_list_so_so(commonPrecSo)
1861
    polySo = sollya_lib_roundcoefficients(polySo, endEllipListSo, None)
1862
    sollya_lib_clear_obj(endEllipListSo)
1863
    #sollya_lib_clear_obj(endEllipListSo)
1864
    return polySo
1865
    
1866
# End pobyso_round_coefficients_single_so_so
1867

    
1868
def pobyso_set_canonical_off():
1869
    sollya_lib_set_canonical(sollya_lib_off())
1870

    
1871
def pobyso_set_canonical_on():
1872
    sollya_lib_set_canonical(sollya_lib_on())
1873

    
1874
def pobyso_set_prec(p):
1875
    """ Legacy function. See pobyso_set_prec_sa_so. """
1876
    pobyso_set_prec_sa_so(p)
1877

    
1878
def pobyso_set_prec_sa_so(p):
1879
    #a = c_int(p)
1880
    #precSo = c_void_p(sollya_lib_constant_from_int(a))
1881
    #precSo = sollya_lib_constant_from_int(a)
1882
    precSo = pobyso_constant_from_int_sa_so(p)
1883
    sollya_lib_set_prec(precSo)
1884
    sollya_lib_clear_obj(precSo)
1885
# End pobyso_set_prec_sa_so.
1886

    
1887
def pobyso_set_prec_so_so(newPrecSo):
1888
    sollya_lib_set_prec(newPrecSo)
1889
# End pobyso_set_prec_so_so.
1890

    
1891
def pobyso_inf_so_so(intervalSo):
1892
    """
1893
    Very thin wrapper around sollya_lib_inf().
1894
    """
1895
    return sollya_lib_inf(intervalSo)
1896
# End pobyso_inf_so_so.
1897
#   
1898
def pobyso_supnorm_so_sa(polySo, funcSo, intervalSo, errorTypeSo = None,\
1899
                         accuracySo = None, realFieldSa = None):
1900
    """
1901
    Computes the supremum norm from Solly input arguments and returns a
1902
    Sage floating-point number whose precision is set by the only Sage argument.
1903
    
1904
    The returned value is the maximum of the absolute values of the range
1905
    elements  returned by the Sollya supnorm functions
1906
    """
1907
    supNormRangeSo = pobyso_supnorm_so_so(polySo, 
1908
                                          funcSo, 
1909
                                          intervalSo, 
1910
                                          errorTypeSo,
1911
                                          accuracySo)
1912
    supNormSo = pobyso_range_max_abs_so_so(supNormRangeSo)
1913
    sollya_lib_clear_obj(supNormRangeSo)
1914
    #pobyso_autoprint(supNormSo)
1915
    supNormSa = pobyso_get_constant_as_rn_with_rf_so_sa(supNormSo, realFieldSa)
1916
    sollya_lib_clear_obj(supNormSo)
1917
    return supNormSa 
1918
# End pobyso_supnorm_so_sa.
1919
#
1920
def pobyso_supnorm_so_so(polySo, funcSo, intervalSo, errorTypeSo = None,\
1921
                         accuracySo = None):
1922
    """
1923
    Computes the supnorm of the approximation error between the given 
1924
    polynomial and function. Attention: returns a range!
1925
    errorTypeSo defaults to "absolute".
1926
    accuracySo defaults to 2^(-40).
1927
    """
1928
    if errorTypeSo is None:
1929
        errorTypeSo = sollya_lib_absolute(None)
1930
        errorTypeIsNone = True
1931
    else:
1932
        errorTypeIsNone = False
1933
    #
1934
    if accuracySo is None:
1935
        # Notice the **: we are in Pythonland here!
1936
        accuracySo = pobyso_constant_sa_so(RR(2**(-40)))
1937
        accuracyIsNone = True
1938
    else:
1939
        accuracyIsNone = False
1940
    #pobyso_autoprint(accuracySo)
1941
    resultSo = \
1942
        sollya_lib_supnorm(polySo, funcSo, intervalSo, errorTypeSo, \
1943
                              accuracySo)
1944
    if errorTypeIsNone:
1945
        sollya_lib_clear_obj(errorTypeSo)
1946
    if accuracyIsNone:
1947
        sollya_lib_clear_obj(accuracySo)
1948
    return resultSo
1949
# End pobyso_supnorm_so_so
1950
#
1951
def pobyso_taylor_expansion_no_change_var_so_so(functionSo, 
1952
                                                degreeSo, 
1953
                                                rangeSo,
1954
                                                errorTypeSo=None, 
1955
                                                sollyaPrecSo=None):
1956
    """
1957
    Compute the Taylor expansion without the variable change
1958
    x -> x-intervalCenter.
1959
    """
1960
    # Change internal Sollya precision, if needed.
1961
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
1962
    sollyaPrecChanged = False
1963
    if sollyaPrecSo is None:
1964
        pass
1965
    else:
1966
        sollya_lib_set_prec(sollyaPrecSo)
1967
        sollyaPrecChanged = True
1968
    # Error type stuff: default to absolute.
1969
    if errorTypeSo is None:
1970
        errorTypeIsNone = True
1971
        errorTypeSo = sollya_lib_absolute(None)
1972
    else:
1973
        errorTypeIsNone = False
1974
    intervalCenterSo = sollya_lib_mid(rangeSo, None)
1975
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo,
1976
                                         intervalCenterSo,
1977
                                         rangeSo, errorTypeSo, None)
1978
    # taylorFormListSaSo is a Python list of Sollya objects references that 
1979
    # are copies of the elements of taylorFormSo.
1980
    # pobyso_get_list_elements_so_so clears taylorFormSo.
1981
    (taylorFormListSaSo, numElementsSa, isEndEllipticSa) = \
1982
        pobyso_get_list_elements_so_so(taylorFormSo)
1983
    polySo = sollya_lib_copy_obj(taylorFormListSaSo[0])
1984
    #print "Num elements:", numElementsSa
1985
    sollya_lib_clear_obj(taylorFormSo)
1986
    #polySo = taylorFormListSaSo[0]
1987
    #errorRangeSo = sollya_lib_copy_obj(taylorFormListSaSo[2])
1988
    errorRangeSo = taylorFormListSaSo[2]
1989
    # No copy_obj needed here: a new objects are created.
1990
    maxErrorSo    = sollya_lib_sup(errorRangeSo)
1991
    minErrorSo    = sollya_lib_inf(errorRangeSo)
1992
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
1993
    absMinErrorSo = sollya_lib_abs(minErrorSo)
1994
    sollya_lib_clear_obj(maxErrorSo)
1995
    sollya_lib_clear_obj(minErrorSo)
1996
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
1997
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
1998
    # If changed, reset the Sollya working precision.
1999
    if sollyaPrecChanged:
2000
        sollya_lib_set_prec(initialSollyaPrecSo)
2001
    sollya_lib_clear_obj(initialSollyaPrecSo)
2002
    if errorTypeIsNone:
2003
        sollya_lib_clear_obj(errorTypeSo)
2004
    pobyso_clear_taylorform_sa_so(taylorFormListSaSo)
2005
    if absMaxErrorSa > absMinErrorSa:
2006
        sollya_lib_clear_obj(absMinErrorSo)
2007
        return (polySo, intervalCenterSo, absMaxErrorSo)
2008
    else:
2009
        sollya_lib_clear_obj(absMaxErrorSo)
2010
        return (polySo, intervalCenterSo, absMinErrorSo)
2011
# end pobyso_taylor_expansion_no_change_var_so_so
2012

    
2013
def pobyso_taylor_expansion_with_change_var_so_so(functionSo, degreeSo, \
2014
                                                  rangeSo, \
2015
                                                  errorTypeSo=None, \
2016
                                                  sollyaPrecSo=None):
2017
    """
2018
    Compute the Taylor expansion with the variable change
2019
    x -> (x-intervalCenter) included.
2020
    """
2021
    # Change Sollya internal precision, if need.
2022
    sollyaPrecChanged = False
2023
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_sos_sa()
2024
    if sollyaPrecSo is None:
2025
        pass
2026
    else:
2027
        sollya_lib_set_prec(sollyaPrecSo)
2028
        sollyaPrecChanged = True
2029
    #
2030
    # Error type stuff: default to absolute.
2031
    if errorTypeSo is None:
2032
        errorTypeIsNone = True
2033
        errorTypeSo = sollya_lib_absolute(None)
2034
    else:
2035
        errorTypeIsNone = False
2036
    intervalCenterSo = sollya_lib_mid(rangeSo)
2037
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo, \
2038
                                         intervalCenterSo, \
2039
                                         rangeSo, errorTypeSo, None)
2040
    # taylorFormListSaSo is a Python list of Sollya objects references that 
2041
    # are copies of the elements of taylorFormSo.
2042
    # pobyso_get_list_elements_so_so clears taylorFormSo.
2043
    (taylorFormListSo, numElements, isEndElliptic) = \
2044
        pobyso_get_list_elements_so_so(taylorFormSo)
2045
    polySo = taylorFormListSo[0]
2046
    errorRangeSo = taylorFormListSo[2]
2047
    maxErrorSo    = sollya_lib_sup(errorRangeSo)
2048
    minErrorSo    = sollya_lib_inf(errorRangeSo)
2049
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
2050
    absMinErrorSo = sollya_lib_abs(minErrorSo)
2051
    sollya_lib_clear_obj(maxErrorSo)
2052
    sollya_lib_clear_obj(minErrorSo)
2053
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
2054
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
2055
    changeVarExpSo = sollya_lib_build_function_sub(\
2056
                       sollya_lib_build_function_free_variable(),\
2057
                       sollya_lib_copy_obj(intervalCenterSo))
2058
    polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpSo)
2059
    sollya_lib_clear_obj(polySo) 
2060
    sollya_lib_clear_obj(changeVarExpSo)
2061
    # If changed, reset the Sollya working precision.
2062
    if sollyaPrecChanged:
2063
        sollya_lib_set_prec(initialSollyaPrecSo)
2064
    sollya_lib_clear_obj(initialSollyaPrecSo)
2065
    if errorTypeIsNone:
2066
        sollya_lib_clear_obj(errorTypeSo)
2067
    sollya_lib_clear_obj(taylorFormSo)
2068
    # Do not clear maxErrorSo.
2069
    if absMaxErrorSa > absMinErrorSa:
2070
        sollya_lib_clear_obj(absMinErrorSo)
2071
        return((polyVarChangedSo, intervalCenterSo, absMaxErrorSo))
2072
    else:
2073
        sollya_lib_clear_obj(absMaxErrorSo)
2074
        return((polyVarChangedSo, intervalCenterSo, absMinErrorSo))
2075
# end pobyso_taylor_expansion_with_change_var_so_so
2076

    
2077
def pobyso_taylor(function, degree, point):
2078
    """ Legacy function. See pobysoTaylor_so_so. """
2079
    return(pobyso_taylor_so_so(function, degree, point))
2080

    
2081
def pobyso_taylor_so_so(functionSo, degreeSo, pointSo):
2082
    return(sollya_lib_taylor(functionSo, degreeSo, pointSo))
2083
    
2084
def pobyso_taylorform(function, degree, point = None, 
2085
                      interval = None, errorType=None):
2086
    """ Legacy function. See pobyso_taylorform_sa_sa;"""
2087
    
2088
def pobyso_taylorform_sa_sa(functionSa, \
2089
                            degreeSa, \
2090
                            pointSa, \
2091
                            intervalSa=None, \
2092
                            errorTypeSa=None, \
2093
                            precisionSa=None):
2094
    """
2095
    Compute the Taylor form of 'degreeSa' for 'functionSa' at 'pointSa' 
2096
    for 'intervalSa' with 'errorTypeSa' (a string) using 'precisionSa'. 
2097
    point: must be a Real or a Real interval.
2098
    return the Taylor form as an array
2099
    TODO: take care of the interval and of the point when it is an interval;
2100
          when errorType is not None;
2101
          take care of the other elements of the Taylor form (coefficients 
2102
          errors and delta.
2103
    """
2104
    # Absolute as the default error.
2105
    if errorTypeSa is None:
2106
        errorTypeSo = sollya_lib_absolute()
2107
    elif errorTypeSa == "relative":
2108
        errorTypeSo = sollya_lib_relative()
2109
    elif errortypeSa == "absolute":
2110
        errorTypeSo = sollya_lib_absolute()
2111
    else:
2112
        # No clean up needed.
2113
        return None
2114
    # Global precision stuff
2115
    sollyaPrecisionChangedSa = False
2116
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
2117
    if precisionSa is None:
2118
        precSa = initialSollyaPrecSa
2119
    else:
2120
        if precSa > initialSollyaPrecSa:
2121
            if precSa <= 2:
2122
                print inspect.stack()[0][3], ":precision change <= 2 requested."
2123
            pobyso_set_prec_sa_so(precSa)
2124
            sollyaPrecisionChangedSa = True
2125
    #        
2126
    if len(functionSa.variables()) > 0:
2127
        varSa = functionSa.variables()[0]
2128
        pobyso_name_free_variable_sa_so(str(varSa))
2129
    # In any case (point or interval) the parent of pointSa has a precision
2130
    # method.
2131
    pointPrecSa = pointSa.parent().precision()
2132
    if precSa > pointPrecSa:
2133
        pointPrecSa = precSa
2134
    # In any case (point or interval) pointSa has a base_ring() method.
2135
    pointBaseRingString = str(pointSa.base_ring())
2136
    if re.search('Interval', pointBaseRingString) is None: # Point
2137
        pointSo = pobyso_constant_sa_so(pointSa, pointPrecSa)
2138
    else: # Interval.
2139
        pointSo = pobyso_interval_to_range_sa_so(pointSa, pointPrecSa)
2140
    # Sollyafy the function.
2141
    functionSo = pobyso_parse_string_sa_so(functionSa._assume_str().replace('_SAGE_VAR_', ''))
2142
    if sollya_lib_obj_is_error(functionSo):
2143
        print "pobyso_tailorform: function string can't be parsed!"
2144
        return None
2145
    # Sollyafy the degree
2146
    degreeSo = sollya_lib_constant_from_int(int(degreeSa))
2147
    # Sollyafy the point
2148
    # Call Sollya
2149
    taylorFormSo = \
2150
        sollya_lib_taylorform(functionSo, degreeSo, pointSo, errorTypeSo,\
2151
                                         None)
2152
    sollya_lib_clear_obj(functionSo)
2153
    sollya_lib_clear_obj(degreeSo)
2154
    sollya_lib_clear_obj(pointSo)
2155
    sollya_lib_clear_obj(errorTypeSo)
2156
    (tfsAsList, numElements, isEndElliptic) = \
2157
            pobyso_get_list_elements_so_so(taylorFormSo)
2158
    polySo = tfsAsList[0]
2159
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
2160
    polyRealField = RealField(maxPrecision)
2161
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, polyRealField)
2162
    if sollyaPrecisionChangedSa:
2163
        sollya_lib_set_prec(initialSollyaPrecSo)
2164
    sollya_lib_clear_obj(initialSollyaPrecSo)
2165
    polynomialRing = polyRealField[str(varSa)]
2166
    polySa = polynomial(expSa, polynomialRing)
2167
    taylorFormSa = [polySa]
2168
    # Final clean-up.
2169
    sollya_lib_clear_obj(taylorFormSo)
2170
    return(taylorFormSa)
2171
# End pobyso_taylor_form_sa_sa
2172

    
2173
def pobyso_taylorform_so_so(functionSo, degreeSo, pointSo, intervalSo=None, \
2174
                            errorTypeSo=None):
2175
    createdErrorType = False
2176
    if errorTypeSo is None:
2177
        errorTypeSo = sollya_lib_absolute()
2178
        createdErrorType = True
2179
    else:
2180
        #TODO: deal with the other case.
2181
        pass
2182
    if intervalSo is None:
2183
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
2184
                                         errorTypeSo, None)
2185
    else:
2186
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
2187
                                         intervalSo, errorTypeSo, None)
2188
    if createdErrorType:
2189
        sollya_lib_clear_obj(errorTypeSo)
2190
    return resultSo
2191
        
2192

    
2193
def pobyso_univar_polynomial_print_reverse(polySa):
2194
    """ Legacy function. See pobyso_univar_polynomial_print_reverse_sa_sa. """
2195
    return(pobyso_univar_polynomial_print_reverse_sa_sa(polySa))
2196

    
2197
def pobyso_univar_polynomial_print_reverse_sa_sa(polySa):
2198
    """
2199
    Return the string representation of a univariate polynomial with
2200
    monomials ordered in the x^0..x^n order of the monomials.
2201
    Remember: Sage
2202
    """
2203
    polynomialRing = polySa.base_ring()
2204
    # A very expensive solution:
2205
    # -create a fake multivariate polynomial field with only one variable,
2206
    #   specifying a negative lexicographical order;
2207
    mpolynomialRing = PolynomialRing(polynomialRing.base(), \
2208
                                     polynomialRing.variable_name(), \
2209
                                     1, order='neglex')
2210
    # - convert the univariate argument polynomial into a multivariate
2211
    #   version;
2212
    p = mpolynomialRing(polySa)
2213
    # - return the string representation of the converted form.
2214
    # There is no simple str() method defined for p's class.
2215
    return(p.__str__())
2216
#
2217
#print pobyso_get_prec()  
2218
pobyso_set_prec(165)
2219
#print pobyso_get_prec()  
2220
#a=100
2221
#print type(a)
2222
#id(a)
2223
#print "Max arity: ", pobyso_max_arity
2224
#print "Function tripleDouble (43) as a string: ", pobyso_function_type_as_string(43)
2225
#print "Function None (44) as a string: ", pobyso_function_type_as_string(44)
2226
sys.stderr.write("\t...Pobyso check done.\n")