Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 258

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

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

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

    
90
pobyso_max_arity = 9
91

    
92
def pobyso_absolute_so_so():
93
    return(sollya_lib_absolute(None))
94

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

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

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

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

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

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

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

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

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

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

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

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

    
317

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    
1016

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

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

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

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

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

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

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

    
1216

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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