Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 307

Historique | Voir | Annoter | Télécharger (103,24 ko)

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

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

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

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

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

31
ToDo (among other things): 
32
 -memory management.
33
"""
34
from ctypes import *
35
import re
36
from sage.symbolic.expression_conversions import polynomial
37
from sage.symbolic.expression_conversions import PolynomialConverter
38
"""
39
Create the equivalent to an enum for the Sollya function types and other
40
common types (error 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,
86
    SOLLYA_ABSOLUTE,
87
    SOLLYA_RELATIVE) = map(int,xrange(46))
88
sys.stderr.write("SOLLYA_RELATIVE = " + str(SOLLYA_RELATIVE) + "\n")
89
sys.stderr.write("Superficial pobyso check...\n")
90
#print "First constant - SOLLYA_BASE_FUNC_ABS: ", SOLLYA_BASE_FUNC_ABS
91
#print "Last constant  - SOLLYA_BASE_FUNC_TRIPLEDOUBLE: ", SOLLYA_BASE_FUNC_TRIPLEDOUBLE
92

    
93
pobyso_max_arity = 9
94

    
95
def pobyso_absolute_so_so():
96
    """
97
    Create an "absolute" Sollya object.
98
    """
99
    return(sollya_lib_absolute(None))
100

    
101
def pobyso_autoprint(arg):
102
    sollya_lib_autoprint(arg, None)
103

    
104
def pobyso_autoprint_so_so(arg):
105
    sollya_lib_autoprint(arg, None)
106

    
107
def pobyso_bounds_to_interval_sa_sa(lowerBound, upperBound):
108
    """
109
    Convert a pair of bounds into an interval (an element of 
110
    a RealIntervalField).
111
    """
112
    # Minimal (not bullet-proof) check on bounds.
113
    if lowerBound > upperBound:
114
        return None
115
    # Try to get the maximum precision among the bounds.
116
    try: 
117
        preclb = parent(lowerBound).precision()
118
        precub = parent(upperBound).precision()
119
        prec = max(preclb, precub)
120
    except AttributeError:
121
        prec = 53
122
    # Create the RealIntervalField and the interval (if possible).
123
    theRIF = RealIntervalField(prec)
124
    theRR  = RealField(prec)
125
    # Try convert the lower bound to a RealField element or fail.
126
    try:
127
        lowerBoundRF = theRR(lowerBound)
128
    except TypeError:
129
        try:
130
            lowerBoundRF = theRR(str(lowerBound))
131
        except TypeError:
132
            return None
133
    # Try convert the upper bound to a RealField element or fail.
134
    try:
135
        upperBoundRF = theRR(upperBound)
136
    except TypeError:
137
        try:
138
            upperBoundRF = theRR(str(upperBound))
139
        except TypeError:
140
            return None
141
    return theRIF(lowerBoundRF, upperBoundRF)
142
# End pobyso_bounds_to_interval_sa_sa
143
    
144
def pobyso_bounds_to_range_sa_so(rnLowerBoundSa, rnUpperBoundSa, \
145
                                 precisionSa=None):
146
    """
147
    Return a Sollya range from to 2 RealField Sage elements.
148
    The Sollya range element has a sufficient precision to hold all
149
    the digits of the widest of the Sage bounds.
150
    """
151
    # Sanity check.
152
    if rnLowerBoundSa > rnUpperBoundSa:
153
        return None
154
    # Precision stuff.
155
    if precisionSa is None:
156
        # Check for the largest precision.
157
        lbPrecSa = rnLowerBoundSa.parent().precision()
158
        ubPrecSa = rnLowerBoundSa.parent().precision()
159
        maxPrecSa = max(lbPrecSa, ubPrecSa)
160
    else:
161
        maxPrecSa = precisionSa
162
    # From Sage to Sollya bounds.
163
#    lowerBoundSo = sollya_lib_constant(get_rn_value(rnLowerBoundSa),
164
#                                       maxPrecSa)
165
    lowerBoundSo = pobyso_constant_sa_so(rnLowerBoundSa,
166
                                         maxPrecSa)
167
    upperBoundSo = pobyso_constant_sa_so(rnUpperBoundSa,
168
                                       maxPrecSa)
169
    
170
    # From Sollya bounds to range.
171
    rangeSo = sollya_lib_range(lowerBoundSo, upperBoundSo)
172
    # Back to original precision.
173
    # Clean up
174
    sollya_lib_clear_obj(lowerBoundSo)
175
    sollya_lib_clear_obj(upperBoundSo)
176
    return rangeSo
177
# End pobyso_bounds_to_range_sa_so
178

    
179
def pobyso_build_end_elliptic_list_so_so(*args):
180
    """
181
    From Sollya argument objects, create a Sollya end elliptic list.
182
    Elements of the list are "eaten" (should not be cleared individually, 
183
    are cleared when the list is cleared).
184
    """
185
    if len(args) == 0:
186
        ## When called with an empty list, sollya_lib_build_end_elliptic_list,
187
        #  produces "error".
188
        return sollya_lib_build_list(None)
189
    index = 0
190
    ## One can not append elements to an elliptic list, prepend only is 
191
    #  permitted.
192
    for argument in reversed(args):
193
        if index == 0:
194
            listSo = sollya_lib_build_end_elliptic_list(argument, None)
195
        else:
196
            listSo = sollya_lib_prepend(argument, listSo)
197
        index += 1
198
    return listSo
199
        
200
# End pobyso_build_end_elliptic_list_so_so
201

    
202
def pobyso_build_function_sub_so_so(exp1So, exp2So):
203
    return sollya_lib_build_function_sub(exp1So, exp2So)
204

    
205
def pobyso_build_list_of_ints_sa_so(*args):
206
    """
207
    Build a Sollya list from Sage integral arguments. 
208
    """
209
    if len(args) == 0:
210
        return pobyso_build_list_so_so()
211
    ## Make a Sage list of integral constants.
212
    intsList = []
213
    for intElem in args:
214
        intsList.append(pobyso_constant_from_int_sa_so(intElem))
215
    return pobyso_build_list_so_so(*intsList) 
216
    
217
def pobyso_build_list_so_so(*args):
218
    """
219
    Make a Sollya list out of Sollya objects passed as arguments.
220
    If one wants to call it with a list argument, should prepend a "*"
221
    before the list variable name.
222
    Elements of the list are "eaten" (should not be cleared individually, 
223
    are cleared when the list is cleared).
224
    """
225
    if len(args) == 0:
226
        ## Called with an empty list produced "error".
227
        return sollya_lib_build_list(None)
228
    index = 0
229
    ## Append the Sollya elements one by one.
230
    for elementSo in args:
231
        if index == 0:
232
            listSo = sollya_lib_build_list(elementSo, None)
233
        else:
234
            listSo = sollya_lib_append(listSo, elementSo)
235
        index += 1
236
    return listSo
237
# End pobyso_build list_so_so
238
    
239

    
240
def pobyso_change_var_in_function_so_so(funcSo, chvarExpSo):
241
    """
242
    Variable change in a function.
243
    """
244
    return(sollya_lib_evaluate(funcSo,chvarExpSo))
245
# End pobyso_change_var_in_function_so_so     
246

    
247
def pobyso_chebyshevform_sa_sa(funcSa, degreeSa, intervalSa):
248
    """
249
    An chebyshevnorm call with Sage arguments that returns Sage objects.
250
    To the moment only the polynomial and the delta are returned.
251
    """
252
    # Check that funcSa is a function.
253
    if not \
254
     sage.symbolic.callable.is_CallableSymbolicExpressionRing(parent(funcSa)):
255
        return None
256
    # Check that funcSa is an integer.
257
    try:
258
        if degreeSa != int(degreeSa):
259
            print "degreeSa is not a valid integer."
260
            return None
261
    except ValueError:
262
        print "degreeSa is not a number."
263
        return None 
264
    # Create the Sollya degree...
265
    degreeSo = pobyso_constant_from_int_sa_so(degreeSa)
266
    # ...and check it.
267
    if pobyso_is_error_so_sa(degreeSo):
268
        print "Could not correctly convert degreeSa to Sollya."
269
        return None
270
    # Check that intervalSa is an interval.
271
    try:
272
        intervalSa.upper()
273
    except AttributeError:
274
        print "intervalSa is not an interval."
275
        return None
276
    # Convert the Sage polynomial into a Sollya polynomial.
277
    # Convert the Sage function into a Sollya function.
278
    funcAsStringSa = funcSa._assume_str().replace('_SAGE_VAR_', '')
279
    funcSo = pobyso_parse_string_sa_so(funcAsStringSa)
280
    if not pobyso_obj_is_function_so_sa(funcSo):
281
        sollya_lib_clear_obj(funcSo)
282
        print "The Sollya object created from the function is not a function."
283
        return None
284
    #pobyso_autoprint(funcSo)
285
    # Convert the Sage interval into a Sollya range.
286
    rangeSo = pobyso_interval_to_range_sa_so(intervalSa)
287
    if not pobyso_is_range_so_sa(rangeSo):
288
        sollya_lib_clear_obj(funcSo)
289
        sollya_lib_clear_obj(rangeSo)
290
        print "The Sollya object created from the interval is not a range."
291
        return None
292
    #pobyso_autoprint(rangeSo)
293
    retValSo = sollya_lib_chebyshevform(funcSo,
294
                                        degreeSo,
295
                                        rangeSo, 
296
                                        None)
297
    sollya_lib_clear_obj(funcSo)
298
    sollya_lib_clear_obj(degreeSo)
299
    sollya_lib_clear_obj(rangeSo)
300
    if pobyso_is_error_so_sa(retValSo):
301
        sollya_lib_clear_obj(retValSo)
302
        print "Could not compute the Chebyshev form."
303
        return None
304
    #pobyso_autoprint(retValSo)
305
    ## Get the list of Sollya objects returned by sollya_lib_chebyshevform
306
    retValSa = pobyso_get_list_elements_so_so(retValSo)
307
    if retValSa[1] != 4:
308
        print "Invalid number of elements in the Chebyshev form."
309
        sollya_lib_clear_obj(retValSo)
310
        return None
311
    ## Convert all the element of the retValSa[0] array of Sollya objects 
312
    #  into their Sage counterparts.
313
    polySa = pobyso_float_poly_so_sa(retValSa[0][0])
314
    if not pobyso_is_poly_sa_sa(polySa):
315
        print "Conversion of the polynomial of the Chebyshev form failed."
316
        sollya_lib_clear_obj(retValSo)
317
        return None
318
    deltaSa  = pobyso_get_interval_from_range_so_sa(retValSa[0][2])
319
    if not pobyso_is_interval_sa_sa(deltaSa):
320
        print "Conversion of the delta interval of the Chebyshev form failed."
321
        sollya_lib_clear_obj(retValSo)
322
        return None
323
    sollya_lib_clear_obj(retValSo)
324
    return(polySa,deltaSa)
325
# End pobyso_chebyshevform_sa_sa.
326

    
327
def pobyso_chebyshevform_so_so(functionSo, degreeSo, intervalSo):
328
    resultSo = sollya_lib_chebyshevform(functionSo, degreeSo, intervalSo)
329
    return(resultSo)
330
# End pobyso_chebyshevform_so_so.
331

    
332
def pobyso_clear_full_list_elements_sa_so(objectListSaSo):
333
    """
334
    Clear the elements of list created by the 
335
    pobyso_get_list_elements_so_so function.
336
    objectListSaSo is as follows:
337
    - objectListSaSo[0]: a list of Sollya objects;
338
    - objectListSaSo[1]: the number of elements of the previous list;
339
    - objectListSaSo[2]: an integer that if != 0 states that the list is 
340
                         end-elliptic
341
    The objects to clear are the elements of the objectListSaSo[0] list. 
342
    """
343
    for index in xrange(0, objectListSaSo[1]):
344
        sollya_lib_clear_obj(objectListSaSo[0][index])
345
# End pobyso_clear_full_list_elements_sa_so
346

    
347
def pobyso_clear_list_elements_sa_so(objectListSaSo):
348
    """
349
    Clear the elements of list of references to Sollya objects
350
    """
351
    for index in xrange(0, len(objectListSaSo)):
352
        sollya_lib_clear_obj(objectListSaSo[index])
353
# End pobyso_clear_list_elements_sa_so
354

    
355
def pobyso_clear_obj(objSo):
356
    """
357
    Free a Sollya object's memory.
358
    Very thin wrapper around sollya_lib_clear_obj().
359
    """
360
    sollya_lib_clear_obj(objSo)
361
# End pobyso_clear_obj. 
362

    
363
def pobyso_clear_taylorform_sa_so(taylorFormSaSo):
364
    """
365
    This method is wrapper around pobyso_clear_list_elements_sa_so.
366
    It is a legacy method left here since it may be used in existing code 
367
    where Taylor forms are used as Sage lists obtained by converting
368
    Sollya Taylor forms (a list made of:
369
    - a polynomial;
370
    - a list of intervals enclosing the errors accumulated when computing
371
      the polynomial coefficients;
372
    - a bound on the approximation error between the polynomial and the 
373
      function.)
374
    A Taylor form directly obtained from pobyso_taylorform_so_so is cleared
375
    by sollya_lib_clear_obj. 
376
    """
377
    pobyso_clear_list_elements_sa_so(taylorFormSaSo)
378
# End pobyso_clear_taylorform_sa_so 
379

    
380
def pobyso_cmp(rnArgSa, cteSo):
381
    """
382
    Deprecated, use pobyso_cmp_sa_so_sa instead.
383
    """
384
    print "Deprecated, use pobyso_cmp_sa_so_sa instead."
385
    return pobyso_cmp_sa_so_sa(rnArgSa, cteSo)
386
# End  pobyso_cmp
387
    
388
def pobyso_cmp_sa_so_sa(rnArgSa, cteSo):
389
    """
390
    Compare the MPFR value a RealNumber with that of a Sollya constant.
391
    
392
    Get the value of the Sollya constant into a RealNumber and compare
393
    using MPFR. Could be optimized by working directly with a mpfr_t
394
    for the intermediate number. 
395
    """
396
    # Get the precision of the Sollya constant to build a Sage RealNumber 
397
    # with enough precision.to hold it.
398
    precisionOfCte = c_int(0)
399
    # From the Sollya constant, create a local Sage RealNumber.
400
    sollya_lib_get_prec_of_constant(precisionOfCte, cteSo) 
401
    #print "Precision of constant: ", precisionOfCte
402
    RRRR = RealField(precisionOfCte.value)
403
    rnLocalSa = RRRR(0)
404
    sollya_lib_get_constant(get_rn_value(rnLocalSa), cteSo)
405
    #
406
    ## Compare the Sage RealNumber version of the Sollya constant with rnArg
407
    #  through a direct comparison of underlying MPFR numbers.
408
    return cmp_rn_value(rnArgSa, rnLocal)
409
# End pobyso_smp_sa_so_sa
410

    
411
def pobyso_compute_pos_function_abs_val_bounds_sa_sa(funcSa, lowerBoundSa, \
412
                                                     upperBoundSa):
413
    """
414
    TODO: completely rework and test.
415
    """
416
    pobyso = pobyso_name_free_variable_sa_so(funcSa.variables()[0])
417
    funcSo = pobyso_parse_string(funcSa._assume_str().replace('_SAGE_VAR_', ''))
418
    rangeSo = pobyso_range_sa_so(lowerBoundSa, upperBoundSa)
419
    infnormSo = pobyso_infnorm_so_so(funcSo,rangeSo)
420
    # Sollya return the infnorm as an interval.
421
    fMaxSa = pobyso_get_interval_from_range_so_sa(infnormSo)
422
    # Get the top bound and compute the binade top limit.
423
    fMaxUpperBoundSa = fMaxSa.upper()
424
    binadeTopLimitSa = 2**ceil(fMaxUpperBoundSa.log2())
425
    # Put up together the function to use to compute the lower bound.
426
    funcAuxSo = pobyso_parse_string(str(binadeTopLimitSa) +  \
427
                                    '-(' + f._assume_str().replace('_SAGE_VAR_', '') + ')')
428
    pobyso_autoprint(funcAuxSo)
429
    # Clear the Sollya range before a new call to infnorm and issue the call.
430
    sollya_lib_clear_obj(infnormSo)
431
    infnormSo = pobyso_infnorm_so_so(funcAuxSo,rangeSo)
432
    fMinSa = pobyso_get_interval_from_range_so_sa(infnormSo)
433
    sollya_lib_clear_obj(infnormSo)
434
    fMinLowerBoundSa = binadeTopLimitSa - fMinSa.lower()
435
    # Compute the maximum of the precisions of the different bounds.
436
    maxPrecSa = max([fMinLowerBoundSa.parent().precision(), \
437
                     fMaxUpperBoundSa.parent().precision()])
438
    # Create a RealIntervalField and create an interval with the "good" bounds.
439
    RRRI = RealIntervalField(maxPrecSa)
440
    imageIntervalSa = RRRI(fMinLowerBoundSa, fMaxUpperBoundSa)
441
    # Free the unneeded Sollya objects
442
    sollya_lib_clear_obj(funcSo)
443
    sollya_lib_clear_obj(funcAuxSo)
444
    sollya_lib_clear_obj(rangeSo)
445
    return(imageIntervalSa)
446
# End pobyso_compute_pos_function_abs_val_bounds_sa_sa
447

    
448
def pobyso_compute_precision_decay_ratio_function_sa_so():
449
    """
450
    Compute the precision decay ratio function for polynomial 
451
    coefficient progressive trucation.
452
    """
453
    functionText = """
454
    proc(deg, a, b, we, wq) 
455
    {
456
      k = we * (exp(x/a)-1) + wq * (b*x)^2 + (1-we-wq) * x;
457
      return k/k(d);
458
    };
459
    """
460
    return pobyso_parse_string_sa_so(functionText)
461
# End  pobyso_compute_precision_decay_ratio_function.
462

    
463

    
464
def pobyso_constant(rnArg):
465
    """ Legacy function. See pobyso_constant_sa_so. """
466
    return(pobyso_constant_sa_so(rnArg))
467
    
468
def pobyso_constant_sa_so(rnArgSa, precisionSa=None):
469
    """
470
    Create a Sollya constant from a Sage RealNumber.
471
    The sollya_lib_constant() function creates a constant
472
    with the same precision as the source.
473
    """
474
    ## Precision stuff. If one wants to change precisions,
475
    #  everything takes place in Sage. This only makes
476
    #  sense if one wants to reduce the precision.
477
    # TODO: revisit precision stuff with new technique.
478
    # Check that rnArgSa is a realFiedl element.
479
    try:
480
        rnArgSa.ulp()
481
    except AttributeError:
482
        return pobyso_error_so()
483
    # If a different precision is wanted modify it.
484
    if not precisionSa is None:
485
        RRR = RealField(precisionSa)
486
        rnArgSa = RRR(rnArgSa)
487
    #print rnArgSa, rnArgSa.precision()
488
    # Sollya constant creation takes place here.
489
    return sollya_lib_constant(get_rn_value(rnArgSa))
490
# End pobyso_constant_sa_so
491
     
492
def pobyso_constant_0_sa_so():
493
    """
494
    Obvious.
495
    """
496
    return pobyso_constant_from_int_sa_so(0)
497

    
498
def pobyso_constant_1():
499
    """
500
    Obvious.
501
    Legacy function. See pobyso_constant_so_so. 
502
    """
503
    return pobyso_constant_1_sa_so()
504

    
505
def pobyso_constant_1_sa_so():
506
    """
507
    Obvious.
508
    """
509
    return(pobyso_constant_from_int_sa_so(1))
510

    
511
def pobyso_constant_from_int(anInt):
512
    """ Legacy function. See pobyso_constant_from_int_sa_so. """
513
    return pobyso_constant_from_int_sa_so(anInt)
514

    
515
def pobyso_constant_from_int_sa_so(anInt):
516
    """
517
    Get a Sollya constant from a Sage int.
518
    """
519
    return sollya_lib_constant_from_int64(long(anInt))
520

    
521
def pobyso_constant_from_int_so_sa(constSo):
522
    """
523
    Get a Sage int from a Sollya int constant.
524
    Usefull for precision or powers in polynomials.
525
    """
526
    constSa = c_long(0)
527
    sollya_lib_get_constant_as_int64(byref(constSa), constSo)
528
    return constSa.value
529
# End pobyso_constant_from_int_so_sa
530

    
531
def pobyso_constant_from_mpq_sa_so(rationalSa):
532
    """
533
    Make a Sollya constant from Sage rational.
534
    The Sollya constant is an unevaluated expression.
535
    Hence no precision argument is needed.
536
    It is better to leave this way since Sollya has its own
537
    optimized evaluation mecanism that tries very hard to
538
    return exact values or at least faithful ones.
539
    """
540
    ratExprSo = \
541
        sollya_lib_constant_from_mpq(sgmp_get_rational_value(rationalSa))
542
    return ratExprSo
543
# End pobyso_constant_from_mpq_sa_so.
544

    
545
def pobyso_constant_sollya_prec_sa_so(rnArgSa):
546
    """
547
    Create a Sollya constant from a Sage RealNumber at the
548
    current precision in Sollya.
549
    """
550
    currentSollyaPrecSa = pobyso_get_prec_so_sa()
551
    return pobyso_constant_sa_so(rnArgSa, currentSollyaPrecSa)
552
# End pobyso_constant_sollya_prec_sa_so
553

    
554
def pobyso_end_elliptic_list_so_sa_so(objectsListSo, intCountSa):
555
    """
556
    Create a Sollya end elliptic list made of the objectListSo[0] to
557
     objectsListSo[intCountSa-1] objects.
558
    """     
559
    return sollya_lib_end_elliptic_list(objectSo, int(intCountSa))
560

    
561
def pobyso_error_so():
562
    return sollya_lib_error(None)
563
# End pobyso_error().
564

    
565
def pobyso_evaluate_so_sa(funcSo, argumentSo):
566
    """
567
    Evaluates funcSo for arguemntSo through sollya_lib_evaluate() and return
568
    the result as a Sage object
569
    """
570
    evalSo = sollya_lib_evaluate(funcSo, argumentSo)
571
    if pobyso_is_error_so_sa(evalSo):
572
        return None
573
    if pobyso_is_range_so_sa(evalSo):
574
        retVal = pobyso_range_to_interval_so_sa(evalSo)
575
    else:
576
        retVal = pobyso_get_constant_as_rn(evalSo)
577
    sollya_lib_clear_obj(evalSo)
578
    return retVal
579
# End pobyso_evaluate_so_sa.
580

    
581
def pobyso_evaluate_so_so(funcSo, argumentSo):
582
    """
583
    Evaluates funcSo for arguemntSo through sollya_lib_evaluate().
584
    """
585
    return sollya_lib_evaluate(funcSo, argumentSo)
586
# End pobyso_evaluate_so_so.
587
#
588
def pobyso_diff_so_so(funcSo):
589
    """
590
    Very thin wrapper around sollya_lib_diff.
591
    """
592
    ## TODO: add a check to make sure funcSo is a functional expression.
593
    return sollya_lib_diff(funcSo)
594

    
595
def pobyso_dirty_find_zeros_so_so(funcSo, rangeSo):
596
    """
597
    Thin wrapper over sollya_lib_dirtyfindzeros()
598
    """
599
    return sollya_lib_dirtyfindzeros(funcSo, rangeSo)
600
# End pobys_dirty_find_zeros
601

    
602
def pobyso_dirty_inf_norm_so_so(funcSo, rangeSo, preSo=None):
603
    """
604
    Thin wrapper around sollya_dirtyinfnorm().
605
    """
606
    # TODO: manage the precision change.
607
    
608
    return sollya_lib_dirtyinfnorm(funcSo, rangeSo)
609
# End pobyso_dirty_inf_norm_so_so
610

    
611
def pobyso_find_zeros_so_so(funcSo, rangeSo):
612
    """
613
    Thin wrapper over sollya_lib_findzeros()
614
    """
615
    return sollya_lib_findzeros(funcSo, rangeSo)
616
# End pobys_find_zeros
617

    
618
def pobyso_float_list_so_sa(listSo):
619
    """
620
    Return a Sollya list of floating-point numbers as a Sage list of 
621
    floating-point numbers.
622
    TODO: add a test to make sure that each element of the list is a constant.
623
    """
624
    listSa   = []
625
    ## The function returns none if the list is empty or an error has happened.
626
    retVal = pobyso_get_list_elements_so_so(listSo)
627
    if retVal is None:
628
        return listSa
629
    ## Just in case the interface is changed and an empty list is returned 
630
    #  instead of None.
631
    elif len(retVal) == 0:
632
        return listSa
633
    else:
634
        ## Remember pobyso_get_list_elements_so_so returns more information
635
        #  than just the elements of the list (# elements, is_elliptic)
636
        listSaSo, numElements, isEndElliptic = retVal
637
    ## Return an empty list.
638
    if numElements == 0:
639
        return listSa
640
    ## Search first for the maximum precision of the elements
641
    maxPrecSa = 0
642
    for floatSo in listSaSo:
643
        #pobyso_autoprint(floatSo)
644
        curPrecSa =  pobyso_get_prec_of_constant_so_sa(floatSo)
645
        if curPrecSa > maxPrecSa:
646
            maxPrecSa = curPrecSa
647
    ##
648
    RF = RealField(maxPrecSa)
649
    ##
650
    for floatSo in listSaSo:
651
        listSa.append(pobyso_get_constant_as_rn_with_rf_so_sa(floatSo))
652
    return listSa
653
# End pobyso_float_list_so_sa
654

    
655
def pobyso_float_poly_sa_so(polySa, precSa = None):
656
    """
657
    Create a Sollya polynomial from a Sage RealField polynomial.
658
    """
659
    ## TODO: filter arguments.
660
    ## Precision. If a precision is given, convert the polynomial
661
    #  into the right polynomial field. If not convert it straight
662
    #  to Sollya.
663
    sollyaPrecChanged = False 
664
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
665
    if precSa is None:
666
        precSa = polySa.parent().base_ring().precision()
667
    if (precSa > initialSollyaPrecSa):
668
        assert precSa >= 2, "Precision change <2 requested"
669
        if precSa <= 2:
670
            print inspect.stack()[0][3], ": precision change <= 2 requested"
671
        precSo = pobyso_constant_from_int(precSa)
672
        pobyso_set_prec_so_so(precSo)
673
        sollya_lib_clear_obj(precSo)
674
        sollyaPrecChanged = True
675
    ## Free variable stuff.
676
    freeVariableNameChanged = False 
677
    polyFreeVariableNameSa = \
678
        str(polySa.variables()[0])
679
    currentFreeVariableNameSa = pobyso_get_free_variable_name_so_sa() 
680
    if polyFreeVariableNameSa != currentFreeVariableNameSa:
681
        #print "Free variable names do not match.", polyFreeVariableNameSa
682
        sollya_lib_name_free_variable(polyFreeVariableNameSa)
683
        freeVariableNameChanged = True
684
    ## Get exponents and coefficients.
685
    exponentsSa     = polySa.exponents()
686
    coefficientsSa  = polySa.coefficients()
687
    ## Build the polynomial.
688
    polySo = None
689
    for coefficientSa, exponentSa in zip(coefficientsSa, exponentsSa):
690
        #print coefficientSa.n(prec=precSa), exponentSa
691
        coefficientSo = \
692
            pobyso_constant_sa_so(coefficientSa)
693
        #pobyso_autoprint(coefficientSo)
694
        exponentSo = \
695
            pobyso_constant_from_int_sa_so(exponentSa)
696
        #pobyso_autoprint(exponentSo)
697
        monomialSo = sollya_lib_build_function_pow(
698
                       sollya_lib_build_function_free_variable(),
699
                       exponentSo)
700
        polyTermSo = sollya_lib_build_function_mul(coefficientSo,
701
                                                       monomialSo)
702
        if polySo is None:
703
            polySo = polyTermSo
704
        else:
705
            polySo = sollya_lib_build_function_add(polySo, polyTermSo)
706
    if sollyaPrecChanged:
707
        pobyso_set_prec_so_so(initialSollyaPrecSo)
708
    sollya_lib_clear_obj(initialSollyaPrecSo)
709
    ## Do not set back the free variable name in Sollya to its initial value: 
710
    #  it will change it back, in the Sollya polynomial, to what it was in the 
711
    #  first place.
712
    return polySo
713
# End pobyso_float_poly_sa_so    
714

    
715
def pobyso_float_poly_so_sa(polySo, realFieldSa=None):
716
    """
717
    Convert a Sollya polynomial into a Sage floating-point polynomial.
718
    If no realField is given, a RealField corresponding to the maximum 
719
    precision of the coefficients is internally computed.
720
    The real field is not returned but can be easily retrieved from 
721
    the polynomial itself.
722
    ALGORITHM:
723
    - (optional) compute the RealField of the coefficients;
724
    - convert the Sollya expression into a Sage expression;
725
    - convert the Sage expression into a Sage polynomial
726
    """    
727
    if realFieldSa is None:
728
        expressionPrecSa = pobyso_get_max_prec_of_exp_so_sa(polySo)
729
        #print "Maximum precision of Sollya polynomial coefficients:", expressionPrecSa
730
        if expressionPrecSa < 2 or expressionPrecSa > 2147483391:
731
            print "Maximum degree of expression:", expressionPrecSa
732
        realFieldSa      = RealField(expressionPrecSa)
733
    #print "Sollya expression before...",
734
    #pobyso_autoprint(polySo)
735

    
736
    expressionSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo,
737
                                                             realFieldSa)
738
    #print "...Sollya expression after."
739
    #pobyso_autoprint(polySo)
740
    polyVariableSa = expressionSa.variables()[0]
741
    polyRingSa     = realFieldSa[str(polyVariableSa)]
742
    #print polyRingSa
743
    # Do not use the polynomial(expressionSa, ring=polyRingSa) form!
744
    polynomialSa = polyRingSa(expressionSa)
745
    #polyCoeffsListSa = polynomialSa.coefficients()
746
    #for coeff in polyCoeffsListSa:
747
    #    print coeff.abs().n()
748
    return polynomialSa
749
# End pobyso_float_poly_so_sa
750

    
751
def pobyso_free_variable():
752
    """
753
    Ultra thin wrapper around the sollya_lib_function_build_free_variable function.
754
    """
755
    return sollya_lib_build_function_free_variable()
756
    
757
def pobyso_function_type_as_string(funcType):
758
    """ Legacy function. See pobyso_function_type_as_string_so_sa. """
759
    return(pobyso_function_type_as_string_so_sa(funcType))
760

    
761
def pobyso_function_type_as_string_so_sa(funcType):
762
    """
763
    Numeric Sollya function codes -> Sage mathematical function names.
764
    Notice that pow -> ^ (a la Sage, not a la Python).
765
    """
766
    if funcType == SOLLYA_BASE_FUNC_ABS:
767
        return "abs"
768
    elif funcType == SOLLYA_BASE_FUNC_ACOS:
769
        return "arccos"
770
    elif funcType == SOLLYA_BASE_FUNC_ACOSH:
771
        return "arccosh"
772
    elif funcType == SOLLYA_BASE_FUNC_ADD:
773
        return "+"
774
    elif funcType == SOLLYA_BASE_FUNC_ASIN:
775
        return "arcsin"
776
    elif funcType == SOLLYA_BASE_FUNC_ASINH:
777
        return "arcsinh"
778
    elif funcType == SOLLYA_BASE_FUNC_ATAN:
779
        return "arctan"
780
    elif funcType == SOLLYA_BASE_FUNC_ATANH:
781
        return "arctanh"
782
    elif funcType == SOLLYA_BASE_FUNC_CEIL:
783
        return "ceil"
784
    elif funcType == SOLLYA_BASE_FUNC_CONSTANT:
785
        return "cte"
786
    elif funcType == SOLLYA_BASE_FUNC_COS:
787
        return "cos"
788
    elif funcType == SOLLYA_BASE_FUNC_COSH:
789
        return "cosh"
790
    elif funcType == SOLLYA_BASE_FUNC_DIV:
791
        return "/"
792
    elif funcType == SOLLYA_BASE_FUNC_DOUBLE:
793
        return "double"
794
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEDOUBLE:
795
        return "doubleDouble"
796
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEEXTENDED:
797
        return "doubleDxtended"
798
    elif funcType == SOLLYA_BASE_FUNC_ERF:
799
        return "erf"
800
    elif funcType == SOLLYA_BASE_FUNC_ERFC:
801
        return "erfc"
802
    elif funcType == SOLLYA_BASE_FUNC_EXP:
803
        return "exp"
804
    elif funcType == SOLLYA_BASE_FUNC_EXP_M1:
805
        return "expm1"
806
    elif funcType == SOLLYA_BASE_FUNC_FLOOR:
807
        return "floor"
808
    elif funcType == SOLLYA_BASE_FUNC_FREE_VARIABLE:
809
        return "freeVariable"
810
    elif funcType == SOLLYA_BASE_FUNC_HALFPRECISION:
811
        return "halfPrecision"
812
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYCONSTANT:
813
        return "libraryConstant"
814
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYFUNCTION:
815
        return "libraryFunction"
816
    elif funcType == SOLLYA_BASE_FUNC_LOG:
817
        return "log"
818
    elif funcType == SOLLYA_BASE_FUNC_LOG_10:
819
        return "log10"
820
    elif funcType == SOLLYA_BASE_FUNC_LOG_1P:
821
        return "log1p"
822
    elif funcType == SOLLYA_BASE_FUNC_LOG_2:
823
        return "log2"
824
    elif funcType == SOLLYA_BASE_FUNC_MUL:
825
        return "*"
826
    elif funcType == SOLLYA_BASE_FUNC_NEARESTINT:
827
        return "round"
828
    elif funcType == SOLLYA_BASE_FUNC_NEG:
829
        return "__neg__"
830
    elif funcType == SOLLYA_BASE_FUNC_PI:
831
        return "pi"
832
    elif funcType == SOLLYA_BASE_FUNC_POW:
833
        return "^"
834
    elif funcType == SOLLYA_BASE_FUNC_PROCEDUREFUNCTION:
835
        return "procedureFunction"
836
    elif funcType == SOLLYA_BASE_FUNC_QUAD:
837
        return "quad"
838
    elif funcType == SOLLYA_BASE_FUNC_SIN:
839
        return "sin"
840
    elif funcType == SOLLYA_BASE_FUNC_SINGLE:
841
        return "single"
842
    elif funcType == SOLLYA_BASE_FUNC_SINH:
843
        return "sinh"
844
    elif funcType == SOLLYA_BASE_FUNC_SQRT:
845
        return "sqrt"
846
    elif funcType == SOLLYA_BASE_FUNC_SUB:
847
        return "-"
848
    elif funcType == SOLLYA_BASE_FUNC_TAN:
849
        return "tan"
850
    elif funcType == SOLLYA_BASE_FUNC_TANH:
851
        return "tanh"
852
    elif funcType == SOLLYA_BASE_FUNC_TRIPLEDOUBLE:
853
        return "tripleDouble"
854
    else:
855
        return None
856

    
857
def pobyso_get_constant(rnArgSa, constSo):
858
    """ Legacy function. See pobyso_get_constant_so_sa. """
859
    return pobyso_get_constant_so_sa(rnArgSa, constSo)
860
# End pobyso_get_constant
861

    
862
def pobyso_get_constant_so_sa(rnArgSa, constSo):
863
    """
864
    Set the value of rnArgSa to the value of constSo in MPFR_RNDN mode.
865
    rnArg must already exist and belong to some RealField.
866
    We assume that constSo points to a Sollya constant.
867
    """
868
    outcome = sollya_lib_get_constant(get_rn_value(rnArgSa), constSo)
869
    if outcome == 0: # Failure because constSo is not a constant expression.
870
        return None
871
    else:
872
        return outcome
873
# End  pobyso_get_constant_so_sa
874
   
875
def pobyso_get_constant_as_rn(ctExpSo):
876
    """ 
877
    Legacy function. See pobyso_get_constant_as_rn_so_sa. 
878
    """ 
879
    return(pobyso_get_constant_as_rn_so_sa(ctExpSo))
880
    
881
def pobyso_get_constant_as_rn_so_sa(constExpSo):
882
    """
883
    Get a Sollya constant as a Sage "real number".
884
    The precision of the floating-point number returned is that of the Sollya
885
    constant.
886
    """
887
    #print "Before computing precision of variable..."
888
    #pobyso_autoprint(constExpSo)
889
    precisionSa  = pobyso_get_prec_of_constant_so_sa(constExpSo)
890
    #print "precisionSa:", precisionSa
891
    ## If the expression can not be exactly converted, None is returned.
892
    #  In this case opt for the Sollya current expression.
893
    if precisionSa is None:
894
        precisionSa = pobyso_get_prec_so_sa()
895
    RRRR = RealField(precisionSa)
896
    rnSa = RRRR(0)
897
    outcome = sollya_lib_get_constant(get_rn_value(rnSa), constExpSo)
898
    if outcome == 0:
899
        return None
900
    else:
901
        return rnSa
902
# End pobyso_get_constant_as_rn_so_sa
903

    
904
def pobyso_get_constant_as_rn_with_rf(ctExp, realField):
905
    """ 
906
    Legacy function. See pobyso_get_constant_as_rn_with_rf_so_sa.
907
    """
908
    return pobyso_get_constant_as_rn_with_rf_so_sa(ctExp, realField)
909
# End pobyso_get_constant_as_rn_with_rf
910
    
911
def pobyso_get_constant_as_rn_with_rf_so_sa(ctExpSo, realFieldSa = None):
912
    """
913
    Get a Sollya constant as a Sage "real number".
914
    If no real field is specified, the precision of the floating-point number 
915
    returned is that of the Sollya constant.
916
    Otherwise is is that of the real field. Hence rounding may happen.
917
    """
918
    if realFieldSa is None:
919
        return pobyso_get_constant_as_rn_so_sa(ctExpSo)
920
    rnSa = realFieldSa(0)
921
    outcome = sollya_lib_get_constant(get_rn_value(rnSa), ctExpSo)
922
    if outcome == 0:
923
        return None
924
    else:
925
        return rnSa
926
# End pobyso_get_constant_as_rn_with_rf_so_sa
927

    
928
def pobyso_get_free_variable_name():
929
    """ 
930
    Legacy function. See pobyso_get_free_variable_name_so_sa.
931
    """
932
    return(pobyso_get_free_variable_name_so_sa())
933

    
934
def pobyso_get_free_variable_name_so_sa():
935
    return sollya_lib_get_free_variable_name()
936
    
937
def pobyso_get_function_arity(expressionSo):
938
    """ 
939
    Legacy function. See pobyso_get_function_arity_so_sa.
940
    """
941
    return(pobyso_get_function_arity_so_sa(expressionSo))
942

    
943
def pobyso_get_function_arity_so_sa(expressionSo):
944
    arity = c_int(0)
945
    sollya_lib_get_function_arity(byref(arity),expressionSo)
946
    return int(arity.value)
947

    
948
def pobyso_get_head_function(expressionSo):
949
    """ 
950
    Legacy function. See pobyso_get_head_function_so_sa. 
951
    """
952
    return(pobyso_get_head_function_so_sa(expressionSo)) 
953

    
954
def pobyso_get_head_function_so_sa(expressionSo):
955
    functionType = c_int(0)
956
    sollya_lib_get_head_function(byref(functionType), expressionSo)
957
    return int(functionType.value)
958

    
959
def pobyso_get_interval_from_range_so_sa(soRange, realIntervalFieldSa = None ):
960
    """
961
    Return the Sage interval corresponding to the Sollya range argument.
962
    If no reaIntervalField is passed as an argument, the interval bounds are not
963
    rounded: they are elements of RealIntervalField of the "right" precision
964
    to hold all the digits.
965
    """
966
    prec = c_int(0)
967
    if realIntervalFieldSa is None:
968
        retval = sollya_lib_get_prec_of_range(byref(prec), soRange, None)
969
        if retval == 0:
970
            return None
971
        realIntervalFieldSa = RealIntervalField(prec.value)
972
    intervalSa = realIntervalFieldSa(0,0)
973
    retval = \
974
        sollya_lib_get_interval_from_range(get_interval_value(intervalSa),\
975
                                           soRange)
976
    if retval == 0:
977
        return None
978
    return intervalSa
979
# End pobyso_get_interval_from_range_so_sa
980

    
981
def pobyso_get_list_elements(soObj):
982
    """ Legacy function. See pobyso_get_list_elements_so_so. """
983
    return pobyso_get_list_elements_so_so(soObj)
984
 
985
def pobyso_get_list_elements_so_so(objectListSo):
986
    """
987
    Get the Sollya list elements as a Sage/Python array of Sollya objects.
988
    
989
    INPUT:
990
    - objectListSo: a Sollya list of Sollya objects.
991
    
992
    OUTPUT:
993
    - a Sage/Python tuple made of:
994
      - a Sage/Python list of Sollya objects,
995
      - a Sage/Python int holding the number of elements,
996
      - a Sage/Python int stating (!= 0) that the list is end-elliptic.
997
    NOTE::
998
        We recover the addresses of the Sollya object from the list of pointers
999
        returned by sollya_lib_get_list_elements. The list itself is freed.
1000
    TODO::
1001
        Figure out what to do with numElements since the number of elements
1002
        can easily be recovered from the list itself. 
1003
        Ditto for isEndElliptic.
1004
    """
1005
    listAddress = POINTER(c_longlong)()
1006
    numElements = c_int(0)
1007
    isEndElliptic = c_int(0)
1008
    listAsSageList = []
1009
    result = sollya_lib_get_list_elements(byref(listAddress),\
1010
                                          byref(numElements),\
1011
                                          byref(isEndElliptic),\
1012
                                          objectListSo)
1013
    if result == 0 :
1014
        return None
1015
    for i in xrange(0, numElements.value, 1):
1016
       #listAsSageList.append(sollya_lib_copy_obj(listAddress[i]))
1017
       listAsSageList.append(listAddress[i])
1018
       # Clear each of the elements returned by Sollya.
1019
       #sollya_lib_clear_obj(listAddress[i])
1020
    # Free the list itself.   
1021
    sollya_lib_free(listAddress)
1022
    return (listAsSageList, numElements.value, isEndElliptic.value)
1023

    
1024
def pobyso_get_max_prec_of_exp(soExp):
1025
    """ Legacy function. See pobyso_get_max_prec_of_exp_so_sa. """
1026
    return pobyso_get_max_prec_of_exp_so_sa(soExp)
1027

    
1028
def pobyso_get_max_prec_of_exp_so_sa(expSo):
1029
    """
1030
    Get the maximum precision used for the numbers in a Sollya expression.
1031
    
1032
    Arguments:
1033
    soExp -- a Sollya expression pointer
1034
    Return value:
1035
    A Python integer
1036
    TODO: 
1037
    - error management;
1038
    - correctly deal with numerical type such as DOUBLEEXTENDED.
1039
    """
1040
    if expSo is None:
1041
        print inspect.stack()[0][3], ": expSo is None."
1042
        return 0
1043
    maxPrecision = 0
1044
    minConstPrec = 0
1045
    currentConstPrec = 0
1046
    #pobyso_autoprint(expSo)
1047
    operator = pobyso_get_head_function_so_sa(expSo)
1048
    if (operator != SOLLYA_BASE_FUNC_CONSTANT) and \
1049
    (operator != SOLLYA_BASE_FUNC_FREE_VARIABLE):
1050
        (arity, subexpressions) = pobyso_get_subfunctions_so_sa(expSo)
1051
        for i in xrange(arity):
1052
            maxPrecisionCandidate = \
1053
                pobyso_get_max_prec_of_exp_so_sa(subexpressions[i])
1054
            if maxPrecisionCandidate > maxPrecision:
1055
                maxPrecision = maxPrecisionCandidate
1056
        return maxPrecision
1057
    elif operator == SOLLYA_BASE_FUNC_CONSTANT:
1058
        #minConstPrec = pobyso_get_min_prec_of_constant_so_sa(expSo)
1059
        #currentConstPrec = pobyso_get_min_prec_of_constant_so_sa(soExp)
1060
        #print minConstPrec, " - ", currentConstPrec 
1061
        return pobyso_get_min_prec_of_constant_so_sa(expSo)
1062
    
1063
    elif operator == SOLLYA_BASE_FUNC_FREE_VARIABLE:
1064
        return 0
1065
    else:
1066
        print "pobyso_get_max_prec_of_exp_so_sa: unexepected operator."
1067
        return 0
1068

    
1069
def pobyso_get_min_prec_of_constant_so_sa(constExpSo):
1070
    """
1071
    Get the minimum precision necessary to represent the value of a Sollya
1072
    constant.
1073
    MPFR_MIN_PREC and powers of 2 are taken into account.
1074
    We assume that constExpSo is a pointer to a Sollay constant expression.
1075
    """
1076
    constExpAsRnSa = pobyso_get_constant_as_rn_so_sa(constExpSo)
1077
    return(min_mpfr_size(get_rn_value(constExpAsRnSa)))
1078

    
1079
def pobyso_get_poly_so_sa(polySo, realFieldSa=None):
1080
    """
1081
    Convert a Sollya polynomial into a Sage polynomial.
1082
    Legacy function. Use pobyso_float_poly_so_sa() instead.
1083
    """
1084
    print "Deprecated: use pobyso_float_poly_so_sa() instead."    
1085
    return pobyso_float_poly_so_sa(polySo,realFieldSa)
1086
# End pobyso_get_poly_so_sa
1087

    
1088
def pobyso_get_prec():
1089
    """ Legacy function. See pobyso_get_prec_so_sa(). """
1090
    return pobyso_get_prec_so_sa()
1091

    
1092
def pobyso_get_prec_so():
1093
    """
1094
    Get the current default precision in Sollya.
1095
    The return value is a Sollya object.
1096
    Usefull when modifying the precision back and forth by avoiding
1097
    extra conversions.
1098
    """
1099
    return sollya_lib_get_prec(None)
1100
    
1101
def pobyso_get_prec_so_sa():
1102
    """
1103
    Get the current default precision in Sollya.
1104
    The return value is Sage/Python int.
1105
    """
1106
    precSo = sollya_lib_get_prec()
1107
    precSa = pobyso_constant_from_int_so_sa(precSo)
1108
    sollya_lib_clear_obj(precSo)
1109
    return precSa
1110
# End pobyso_get_prec_so_sa.
1111

    
1112
def pobyso_get_prec_so_so_sa():
1113
    """
1114
    Return the current precision both as a Sollya object and a
1115
    Sage integer as hybrid tuple.
1116
    To avoid multiple calls for precision manipulations. 
1117
    """
1118
    precSo = sollya_lib_get_prec()
1119
    precSa = pobyso_constant_from_int_so_sa(precSo)
1120
    return (precSo, int(precSa))
1121
    
1122
def pobyso_get_prec_of_constant(ctExpSo):
1123
    """ Legacy function. See pobyso_get_prec_of_constant_so_sa. """
1124
    return pobyso_get_prec_of_constant_so_sa(ctExpSo)
1125

    
1126
def pobyso_get_prec_of_constant_so_sa(ctExpSo):
1127
    """
1128
    Tries to find a precision to represent ctExpSo without rounding.
1129
    If not possible, returns None.
1130
    """
1131
    #print "Entering pobyso_get_prec_of_constant_so_sa..."
1132
    prec = c_int(0)
1133
    retc = sollya_lib_get_prec_of_constant(byref(prec), ctExpSo, None)
1134
    if retc == 0:
1135
        #print "pobyso_get_prec_of_constant_so_sa failed."
1136
        return None
1137
    #print "...exiting pobyso_get_prec_of_constant_so_sa."
1138
    return int(prec.value)
1139

    
1140
def pobyso_get_prec_of_range_so_sa(rangeSo):
1141
    """
1142
    Returns the number of bits elements of a range are coded with.
1143
    """
1144
    prec = c_int(0)
1145
    retc = sollya_lib_get_prec_of_range(byref(prec), rangeSo, None)
1146
    if retc == 0:
1147
        return(None)
1148
    return int(prec.value)
1149
# End pobyso_get_prec_of_range_so_sa()
1150

    
1151
def pobyso_get_sage_exp_from_sollya_exp(sollyaExpSo, realField = RR):
1152
    """ Legacy function. See pobyso_get_sage_exp_from_sollya_exp_so_sa. """
1153
    return pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, 
1154
                                                     realField = RR)
1155

    
1156
def pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, realFieldSa = RR):
1157
    """
1158
    Get a Sage expression from a Sollya expression. 
1159
    Currently only tested with polynomials with floating-point coefficients.
1160
    Notice that, in the returned polynomial, the exponents are RealNumbers.
1161
    """
1162
    #pobyso_autoprint(sollyaExp)
1163
    operatorSa = pobyso_get_head_function_so_sa(sollyaExpSo)
1164
    sollyaLibFreeVariableName = sollya_lib_get_free_variable_name()
1165
    ## Get rid of the "_"'s in "_x_", if any.
1166
    sollyaLibFreeVariableName = re.sub('_', '', sollyaLibFreeVariableName)
1167
    # Constants and the free variable are special cases.
1168
    # All other operator are dealt with in the same way.
1169
    if (operatorSa != SOLLYA_BASE_FUNC_CONSTANT) and \
1170
       (operatorSa != SOLLYA_BASE_FUNC_FREE_VARIABLE):
1171
        (aritySa, subexpressionsSa) = pobyso_get_subfunctions_so_sa(sollyaExpSo)
1172
        if aritySa == 1:
1173
            sageExpSa = eval(pobyso_function_type_as_string_so_sa(operatorSa) + \
1174
            "(" + pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], \
1175
            realFieldSa) + ")")
1176
        elif aritySa == 2:
1177
            # We do not get through the preprocessor.
1178
            # The "^" operator is then a special case.
1179
            if operatorSa == SOLLYA_BASE_FUNC_POW:
1180
                operatorAsStringSa = "**"
1181
            else:
1182
                operatorAsStringSa = \
1183
                    pobyso_function_type_as_string_so_sa(operatorSa)
1184
            sageExpSa = \
1185
              eval("pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], realFieldSa)"\
1186
              + " " + operatorAsStringSa + " " + \
1187
                   "pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[1], realFieldSa)")
1188
        # We do not know yet how to deal with arity >= 3 
1189
        # (is there any in Sollya anyway?).
1190
        else:
1191
            sageExpSa = eval('None')
1192
        return sageExpSa
1193
    elif operatorSa == SOLLYA_BASE_FUNC_CONSTANT:
1194
        #print "This is a constant"
1195
        return pobyso_get_constant_as_rn_with_rf_so_sa(sollyaExpSo, realFieldSa)
1196
    elif operatorSa == SOLLYA_BASE_FUNC_FREE_VARIABLE:
1197
        #print "This is the free variable"
1198
        return eval(sollyaLibFreeVariableName)
1199
    else:
1200
        print "Unexpected"
1201
        return eval('None')
1202
# End pobyso_get_sage_exp_from_sollya_exp_so_sa
1203

    
1204

    
1205
def pobyso_get_subfunctions(expressionSo):
1206
    """ Legacy function. See pobyso_get_subfunctions_so_sa. """
1207
    return pobyso_get_subfunctions_so_sa(expressionSo) 
1208
# End pobyso_get_subfunctions.
1209
 
1210
def pobyso_get_subfunctions_so_sa(expressionSo):
1211
    """
1212
    Get the subfunctions of an expression.
1213
    Return the number of subfunctions and the list of subfunctions addresses.
1214
    S.T.: Could not figure out another way than that ugly list of declarations
1215
    to recover the addresses of the subfunctions.
1216
    We limit ourselves to arity 8 functions. 
1217
    """
1218
    subf0 = c_int(0)
1219
    subf1 = c_int(0)
1220
    subf2 = c_int(0)
1221
    subf3 = c_int(0)
1222
    subf4 = c_int(0)
1223
    subf5 = c_int(0)
1224
    subf6 = c_int(0)
1225
    subf7 = c_int(0)
1226
    subf8 = c_int(0)
1227
    arity = c_int(0)
1228
    nullPtr = POINTER(c_int)()
1229
    sollya_lib_get_subfunctions(expressionSo, byref(arity), \
1230
      byref(subf0), byref(subf1), byref(subf2), byref(subf3), \
1231
      byref(subf4), byref(subf5),\
1232
      byref(subf6), byref(subf7), byref(subf8), nullPtr, None) 
1233
#    byref(cast(subfunctions[0], POINTER(c_int))), \
1234
#    byref(cast(subfunctions[0], POINTER(c_int))), \
1235
#    byref(cast(subfunctions[2], POINTER(c_int))), \
1236
#    byref(cast(subfunctions[3], POINTER(c_int))), \
1237
#    byref(cast(subfunctions[4], POINTER(c_int))), \
1238
#    byref(cast(subfunctions[5], POINTER(c_int))), \
1239
#    byref(cast(subfunctions[6], POINTER(c_int))), \
1240
#    byref(cast(subfunctions[7], POINTER(c_int))), \
1241
#    byref(cast(subfunctions[8], POINTER(c_int))), nullPtr)
1242
    subfunctions = [subf0, subf1, subf2, subf3, subf4, subf5, subf6, subf7, \
1243
                    subf8]
1244
    subs = []
1245
    if arity.value > pobyso_max_arity:
1246
        return(0,[])
1247
    for i in xrange(arity.value):
1248
        subs.append(int(subfunctions[i].value))
1249
        #print subs[i]
1250
    return (int(arity.value), subs)
1251
# End pobyso_get_subfunctions_so_sa
1252
    
1253
def pobyso_guess_degree_sa_sa(functionSa, intervalSa, approxErrorSa, 
1254
                              weightSa=None, degreeBoundSa=None):
1255
    """
1256
    Sa_sa variant of the solly_guessdegree function.
1257
    Return 0 if something goes wrong.
1258
    """
1259
    functionAsStringSa = functionSa._assume_str().replace('_SAGE_VAR_', '')
1260
    functionSo = pobyso_parse_string_sa_so(functionAsStringSa)
1261
    if pobyso_is_error_so_sa(functionSo):
1262
        sollya_lib_clear_obj(functionSo)
1263
        return 0
1264
    rangeSo = pobyso_interval_to_range_sa_so(intervalSa)
1265
    # The approximation error is expected to be a floating point number.
1266
    if pobyso_is_floating_point_number_sa_sa(approxErrorSa):
1267
        approxErrorSo = pobyso_constant_sa_so(approxErrorSa)
1268
    else:
1269
        approxErrorSo = pobyso_constant_sa_so(RR(approxErrorSa))
1270
    if not weightSa is None:
1271
        weightAsStringSa = weightSa._assume_str().replace('_SAGE_VAR_', '')
1272
        weightSo = pobyso_parse_string_sa_so(weightAsStringSa)
1273
        if pobyso_is_error_so_sa(weightSo):
1274
            sollya_lib_clear_obj(functionSo)
1275
            sollya_lib_clear_obj(rangeSo)
1276
            sollya_lib_clear_obj(approxErrorSo)   
1277
            sollya_lib_clear_obj(weightSo)
1278
            return 0   
1279
    else:
1280
        weightSo = None
1281
    if not degreeBoundSa is None:
1282
        degreeBoundSo = pobyso_constant_from_int_sa_so(degreeBoundSa)
1283
    else:
1284
        degreeBoundSo = None
1285
    guessedDegreeSa = pobyso_guess_degree_so_sa(functionSo,
1286
                                                rangeSo,
1287
                                                approxErrorSo,
1288
                                                weightSo,
1289
                                                degreeBoundSo)
1290
    sollya_lib_clear_obj(functionSo)
1291
    sollya_lib_clear_obj(rangeSo)
1292
    sollya_lib_clear_obj(approxErrorSo)
1293
    if not weightSo is None:
1294
        sollya_lib_clear_obj(weightSo)
1295
    if not degreeBoundSo is None:
1296
        sollya_lib_clear_obj(degreeBoundSo)
1297
    return guessedDegreeSa
1298
# End poyso_guess_degree_sa_sa
1299

    
1300
def pobyso_guess_degree_so_sa(functionSo, rangeSo, errorSo, weightSo=None, \
1301
                              degreeBoundSo=None):
1302
    """
1303
    Thin wrapper around the guessdegree function.
1304
    Nevertheless, some precision control stuff has been appended.
1305
    """
1306
    # Deal with Sollya internal precision issues: if it is too small, 
1307
    # compared with the error, increases it to about twice -log2(error).
1308
    errorSa = pobyso_get_constant_as_rn_with_rf_so_sa(errorSo)
1309
    log2ErrorSa = errorSa.log2()
1310
    if log2ErrorSa < 0:
1311
        neededPrecisionSa = int(2 * int(-log2ErrorSa) / 64) * 64
1312
    else:
1313
        neededPrecisionSa = int(2 * int(log2ErrorSa) / 64) * 64
1314
    #print "Needed precision:", neededPrecisionSa
1315
    sollyaPrecisionHasChanged = False
1316
    (initialPrecSo, initialPrecSa) = pobyso_get_prec_so_so_sa()
1317
    if neededPrecisionSa > initialPrecSa:
1318
        if neededPrecisionSa <= 2:
1319
            print inspect.stack()[0][3], ": precision change <= 2 requested."
1320
        pobyso_set_prec_sa_so(neededPrecisionSa)
1321
        sollyaPrecisionHasChanged = True
1322
    #print "Guessing degree..."
1323
    # weightSo and degreeBoundsSo are optional arguments.
1324
    # As declared, sollya_lib_guessdegree must take 5 arguments.
1325
    if weightSo is None:
1326
        degreeRangeSo = sollya_lib_guessdegree(functionSo, rangeSo, errorSo, 
1327
                                               0, 0, None)
1328
    elif degreeBoundSo is None:
1329
        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, \
1330
                                                errorSo, weightSo, 0, None)
1331
    else:
1332
        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, errorSo, \
1333
                                                weightSo, degreeBoundSo, None)
1334
    #print "...degree guess done."
1335
    # Restore internal precision, if applicable.
1336
    if sollyaPrecisionHasChanged:
1337
        pobyso_set_prec_so_so(initialPrecSo)
1338
        sollya_lib_clear_obj(initialPrecSo)
1339
    degreeIntervalSa = pobyso_range_to_interval_so_sa(degreeRangeSo)
1340
    sollya_lib_clear_obj(degreeRangeSo)
1341
    # When ok, both bounds match.
1342
    # When the degree bound is too low, the upper bound is the degree
1343
    # for which the error can be honored.
1344
    # When it really goes wrong, the upper bound is infinity. 
1345
    if degreeIntervalSa.lower() == degreeIntervalSa.upper():
1346
        return int(degreeIntervalSa.lower())
1347
    else:
1348
        if degreeIntervalSa.upper().is_infinity():
1349
            return None
1350
        else:
1351
            return int(degreeIntervalSa.upper())
1352
    # End pobyso_guess_degree_so_sa
1353

    
1354
def pobyso_inf_so_so(intervalSo):
1355
    """
1356
    Very thin wrapper around sollya_lib_inf().
1357
    """
1358
    return sollya_lib_inf(intervalSo)
1359
# End pobyso_inf_so_so.
1360
#    
1361
def pobyso_infnorm_sa_sa(funcSa, intervalSa):
1362
    """
1363
    An infnorm call with Sage arguments.
1364
    We only take into account the 2 first arguments (the function and
1365
    the interval (a range). Managing the other arguments (the file for
1366
    the proof and the exclusion intervals list) will be performed later
1367
    Changes will be needed in sollya_lib.py file too.
1368
    """
1369
    # Check that funcSa is a function.
1370
    if not \
1371
     sage.symbolic.callable.is_CallableSymbolicExpressionRing(parent(funcSa)):
1372
        return None
1373
    # Check that intervalSa is an interval.
1374
    try:
1375
        intervalSa.upper()
1376
    except AttributeError:
1377
        return None
1378
    # Convert the Sage function into a Sollya function.
1379
    funcAsStringSa = funcSa._assume_str().replace('_SAGE_VAR_', '')
1380
    funcSo = pobyso_parse_string_sa_so(funcAsStringSa)
1381
    if not pobyso_obj_is_function_so_sa(funcSo):
1382
        sollya_lib_clear_obj(funcSo)
1383
        return None
1384
    # Convert the Sage interval into a Sollya range.
1385
    rangeSo = pobyso_interval_to_range_sa_so(intervalSa)
1386
    retValSo = sollya_lib_infnorm(funcSo, rangeSo, None)
1387
    sollya_lib_clear_obj(funcSo)
1388
    sollya_lib_clear_obj(rangeSo)
1389
    if pobyso_is_error_so_sa(retValSo):
1390
        sollya_lib_clear_obj(retValSo)
1391
        return None
1392
    retValSa = pobyso_range_to_interval_so_sa(retValSo)
1393
    sollya_lib_clear_obj(retValSo)
1394
    return retValSa
1395
# End pobyso_infnorm_so_so.
1396
#   
1397
def pobyso_infnorm_so_so(funcSo, rangeSo):
1398
    """
1399
    Very thin wrapper around sollya_lib_infnorm().
1400
    We only take into account the 2 first arguments (the function and
1401
    the interval (a range). Managing the other arguments (the file for
1402
    the proof and the exclusion intervals list) will be performed later
1403
    Changes will be needed in sollya_lib.py file too.
1404
    
1405
    As per Sollya manual, this function should not be used anymore and
1406
    supnorm should be called instead. Nevertheless, supnorm breaks 
1407
    sometimes whereas infnorm still returns a satisfactory answer.
1408
    """
1409
    return sollya_lib_infnorm(funcSo, rangeSo, None)
1410
# End pobyso_infnorm_so_so.
1411

    
1412
def pobyso_interval_to_range_sa_so(intervalSa, precisionSa=None):
1413
    if precisionSa is None:
1414
        precisionSa = intervalSa.parent().precision()
1415
    intervalSo = pobyso_bounds_to_range_sa_so(intervalSa.lower(),\
1416
                                              intervalSa.upper(),\
1417
                                              precisionSa)
1418
    return intervalSo
1419
# End pobyso_interval_to_range_sa_so
1420

    
1421
def pobyso_is_error_so_sa(objSo):
1422
    """
1423
    Thin wrapper around the sollya_lib_obj_is_error() function.
1424
    """
1425
    if sollya_lib_obj_is_error(objSo) != 0:
1426
        return True
1427
    else:
1428
        return False
1429
# End pobyso_is_error-so_sa
1430

    
1431
def pobyso_is_floating_point_number_sa_sa(numberSa):
1432
    """
1433
    Check whether a Sage number is floating point.
1434
    Exception stuff added because numbers other than
1435
    floating-point ones do not have the is_real() attribute.
1436
    """
1437
    try:
1438
        return numberSa.is_real()
1439
    except AttributeError:
1440
        return False
1441
# End pobyso_is_floating_piont_number_sa_sa
1442

    
1443
def pobyso_is_interval_sa_sa(objectSa):
1444
    """
1445
    Check that objectSa is an interval.
1446
    """
1447
    try:
1448
        objectSa.upper()
1449
    except AttributeError:
1450
        return False
1451
    return True
1452

    
1453
def pobyso_is_poly_sa_sa(objectSa):
1454
    """
1455
    Check if an object is a polynomial.
1456
    We consider 3 kinds of polynomials:
1457
    - symbolicExpressions that are polynomials;
1458
    - callable symbolicExpression that are polynomials;
1459
    - polynomials over a ring.
1460
    """
1461
    # We consider a constant as a polynomial. All 3 kinds of polynomials
1462
    # have the is_constant() method.
1463
    try:
1464
        if objectSa.is_constant():
1465
            return True
1466
    except AttributeError:
1467
        return False
1468
    # All 3 kinds of polynomials have the args() method.
1469
    # A polynomial has at least one argument.
1470
    # We catch the TypeError in case the return value of args() has
1471
    # no __getitem__ method.
1472
    try:
1473
        argument = objectSa.args()[0]
1474
    except (AttributeError, TypeError):
1475
        return False
1476
    # The 2 first kinds have a is_polynomial() method that returns True
1477
    try:
1478
        if not objectSa.is_polynomial(argument):
1479
            return False
1480
    except AttributeError:
1481
        try:
1482
            objectSa.polynomial(argument)
1483
        except AttributeError:
1484
            return False
1485
    return True
1486
# End pobyso_is_poly_sa_sa
1487

    
1488
def pobyso_is_range_so_sa(rangeCandidateSo):
1489
    """
1490
    Thin wrapper over sollya_lib_is_range.
1491
    """
1492
    return sollya_lib_obj_is_range(rangeCandidateSo) != 0
1493
    
1494
# End pobyso_is_range_so_sa
1495

    
1496

    
1497
def pobyso_lib_init():
1498
    sollya_lib_init(None)
1499

    
1500
def pobyso_lib_close():
1501
    sollya_lib_close(None)
1502

    
1503
def pobyso_name_free_variable(freeVariableNameSa):
1504
    """ Legacy function. See pobyso_name_free_variable_sa_so. """
1505
    pobyso_name_free_variable_sa_so(freeVariableNameSa)
1506

    
1507
def pobyso_name_free_variable_sa_so(freeVariableNameSa):
1508
    """
1509
    Set the free variable name in Sollya from a Sage string.
1510
    """
1511
    sollya_lib_name_free_variable(freeVariableNameSa)
1512

    
1513
def pobyso_obj_is_function_so_sa(objSo):
1514
    """
1515
    Check if an object is a function.
1516
    """
1517
    if sollya_lib_obj_is_function(objSo) != 0:
1518
        return True
1519
    else:
1520
        return False
1521
# End pobyso_obj_is_function_so_sa  
1522
    
1523
def pobyso_obj_is_range_so_sa(objSo):
1524
    """
1525
    Check if an object is a function.
1526
    """
1527
    if sollya_lib_obj_is_range(objSo) != 0:
1528
        return True
1529
    else:
1530
        return False
1531
# End pobyso_obj_is_range_so_sa  
1532
    
1533
def pobyso_obj_is_string_so_sa(objSo):
1534
    """
1535
    Check if an object is a function.
1536
    """
1537
    if sollya_lib_obj_is_string(objSo) != 0:
1538
        return True
1539
    else:
1540
        return False
1541
# End pobyso_obj_is_string_so_sa  
1542
    
1543
def pobyso_parse_string(string):
1544
    """ Legacy function. See pobyso_parse_string_sa_so. """
1545
    return pobyso_parse_string_sa_so(string)
1546
 
1547
def pobyso_parse_string_sa_so(string):
1548
    """
1549
    Get the Sollya expression computed from a Sage string or
1550
    a Sollya error object if parsing failed.
1551
    """
1552
    return sollya_lib_parse_string(string)
1553

    
1554
def pobyso_precision_so_sa(ctExpSo):
1555
    """
1556
    Computes the necessary precision to represent a number.
1557
    If x is not zero, it can be uniquely written as x = m · 2e 
1558
    where m is an odd integer and e is an integer. 
1559
    precision(x) returns the number of bits necessary to write m 
1560
    in binary (i.e. ceil(log2(m))).
1561
    """
1562
    #TODO: take care of the special case: 0, @NaN@, @Inf@
1563
    precisionSo = sollya_lib_precision(ctExpSo)
1564
    precisionSa = pobyso_constant_from_int_so_sa(precisionSo)
1565
    sollya_lib_clear_obj(precisionSo)
1566
    return precisionSa
1567
# End pobyso_precision_so_sa
1568

    
1569
def pobyso_polynomial_coefficients_progressive_round_so_so(polySo,
1570
                                                           funcSo,
1571
                                                           icSo,
1572
                                                           intervalSo,
1573
                                                           itpSo,
1574
                                                           ftpSo,
1575
                                                           maxPrecSo,
1576
                                                           maxErrSo,
1577
                                                           debug=False):
1578
    if debug:
1579
        print "Input arguments:"
1580
        pobyso_autoprint(polySo)
1581
        pobyso_autoprint(funcSo)
1582
        pobyso_autoprint(icSo)
1583
        pobyso_autoprint(intervalSo)
1584
        pobyso_autoprint(itpSo)
1585
        pobyso_autoprint(ftpSo)
1586
        pobyso_autoprint(maxPrecSo)
1587
        pobyso_autoprint(maxErrSo)
1588
        print "________________"
1589
    
1590
    ## Higher order function see: 
1591
    #  http://effbot.org/pyfaq/how-do-you-make-a-higher-order-function-in-python.htm
1592
    def precision_decay_ratio_function(degreeSa):
1593
        def outer(x):
1594
            def inner(x):
1595
                we = 3/8
1596
                wq = 2/8
1597
                a  = 2.2
1598
                b  = 2 
1599
                return we*(exp(x/a)-1) +  wq*((b*x)**2) + (1-we-wq)*x
1600
            return  inner(x)/inner(degreeSa)
1601
        return outer
1602
    
1603
    #
1604
    degreeSa        = pobyso_polynomial_degree_so_sa(polySo)
1605
    ratio           = precision_decay_ratio_function(degreeSa)
1606
    itpSa           = pobyso_constant_from_int_so_sa(itpSo)
1607
    ftpSa           = pobyso_constant_from_int_so_sa(ftpSo)
1608
    maxPrecSa       = pobyso_constant_from_int_so_sa(maxPrecSo)
1609
    maxErrSa        = pobyso_get_constant_as_rn_so_sa(maxErrSo)
1610
    if debug:
1611
        print "degreeSa:", degreeSa
1612
        print "ratio:", ratio
1613
        print "itpsSa:", itpSa
1614
        print "ftpSa:", ftpSa
1615
        print "maxPrecSa:", maxPrecSa
1616
        print "maxErrSa:", maxErrSa
1617
    lastResPolySo   = None
1618
    lastInfNormSo   = None
1619
    #print "About to enter the while loop..."
1620
    while True: 
1621
        resPolySo   = pobyso_constant_0_sa_so()
1622
        pDeltaSa    = ftpSa - itpSa
1623
        for indexSa in reversed(xrange(0,degreeSa+1)):
1624
            #print "Index:", indexSa
1625
            indexSo = pobyso_constant_from_int_sa_so(indexSa)
1626
            #pobyso_autoprint(indexSo)
1627
            #print ratio(indexSa)
1628
            ctpSa = floor(ftpSa - (pDeltaSa * ratio(indexSa)))
1629
            ctpSo = pobyso_constant_from_int_sa_so(ctpSa)
1630
            if debug:
1631
                print "Index:", indexSa, " - Target precision:", 
1632
                pobyso_autoprint(ctpSo)
1633
            cmonSo  = \
1634
                sollya_lib_build_function_mul(sollya_lib_coeff(polySo, indexSo),
1635
                                      sollya_lib_build_function_pow( \
1636
                                          sollya_lib_build_function_free_variable(), \
1637
                                          indexSo))
1638
            #pobyso_autoprint(cmonSo)
1639
            cmonrSo = pobyso_round_coefficients_single_so_so(cmonSo, ctpSo)
1640
            sollya_lib_clear_obj(cmonSo)
1641
            #pobyso_autoprint(cmonrSo)
1642
            resPolySo = sollya_lib_build_function_add(resPolySo,
1643
                                                      cmonrSo)
1644
            #pobyso_autoprint(resPolySo)
1645
        # End for index
1646
        freeVarSo     = sollya_lib_build_function_free_variable()
1647
        changeVarSo   = sollya_lib_sub(freeVarSo, icSo)
1648
        resPolyCvSo   = sollya_lib_evaluate(resPolySo, changeVarSo)
1649
        errFuncSo = sollya_lib_build_function_sub(sollya_lib_copy_obj(funcSo), 
1650
                                                  resPolyCvSo)
1651
        infNormSo = sollya_lib_dirtyinfnorm(errFuncSo, intervalSo)
1652
        cerrSa    = pobyso_get_constant_as_rn_so_sa(infNormSo)
1653
        if debug:
1654
            print "Infnorm (Sollya):", 
1655
            pobyso_autoprint(infNormSo)
1656
        sollya_lib_clear_obj(errFuncSo)
1657
        #print "Infnorm  (Sage):", cerrSa
1658
        if (cerrSa > maxErrSa):
1659
            if debug:
1660
                print "Error is too large."
1661
            if lastResPolySo is None:
1662
                if debug:
1663
                    print "Enlarging prec."
1664
                ntpSa = floor(ftpSa + ftpSa/50)
1665
                ## Can't enlarge (numerical)
1666
                if ntpSa == ftpSa:
1667
                    sollya_lib_clear_obj(resPolySo)
1668
                    return None
1669
                ## Can't enlarge (not enough precision left)
1670
                if ntpSa > maxPrecSa:
1671
                    sollya_lib_clear_obj(resPolySo)
1672
                    return None
1673
                ftpSa = ntpSa
1674
                continue
1675
            ## One enlargement took place.
1676
            else:
1677
                if debug:
1678
                    print "Exit with the last before last polynomial."
1679
                    print "Precision of highest degree monomial:", itpSa
1680
                    print "Precision of constant term          :", ftpSa
1681
                sollya_lib_clear_obj(resPolySo)
1682
                sollya_lib_clear_obj(infNormSo)
1683
                return (lastResPolySo, lastInfNormSo)
1684
        # cerrSa <= maxErrSa: scrap more bits, possibly.
1685
        else:
1686
            if debug:
1687
                print "Error is too small"
1688
            if cerrSa <= (maxErrSa/2):
1689
                if debug: 
1690
                    print "Shrinking prec."
1691
                ntpSa = floor(ftpSa - ftpSa/50)
1692
                ## Can't shrink (numerical)
1693
                if ntpSa == ftpSa:
1694
                    if not lastResPolySo is None:
1695
                        sollya_lib_clear_obj(lastResPolySo)
1696
                    if not lastInfNormSo is None:
1697
                        sollya_lib_clear_obj(lastInfNormSo)
1698
                    if debug:
1699
                        print "Exit because can't shrink anymore (numerically)."
1700
                        print "Precision of highest degree monomial:", itpSa
1701
                        print "Precision of constant term          :", ftpSa
1702
                    return (resPolySo, infNormSo)
1703
                ## Can't shrink (not enough precision left)
1704
                if ntpSa <= itpSa:
1705
                    if not lastResPolySo is None:
1706
                        sollya_lib_clear_obj(lastResPolySo)
1707
                    if not lastInfNormSo is None:
1708
                        sollya_lib_clear_obj(lastInfNormSo)
1709
                        print "Exit because can't shrink anymore (no bits left)."
1710
                        print "Precision of highest degree monomial:", itpSa
1711
                        print "Precision of constant term          :", ftpSa
1712
                    return (resPolySo, infNormSo)
1713
                ftpSa = ntpSa
1714
                if not lastResPolySo is None:
1715
                    sollya_lib_clear_obj(lastResPolySo)
1716
                if not lastInfNormSo is None:
1717
                    sollya_lib_clear_obj(lastInfNormSo)
1718
                lastResPolySo = resPolySo
1719
                lastInfNormSo = infNormSo
1720
                continue
1721
            else: # Error is not that small, just return
1722
                if not lastResPolySo is None:
1723
                    sollya_lib_clear_obj(lastResPolySo)
1724
                if not lastInfNormSo is None:
1725
                    sollya_lib_clear_obj(lastInfNormSo)
1726
                if debug:
1727
                    print "Exit normally."
1728
                    print "Precision of highest degree monomial:", itpSa
1729
                    print "Precision of constant term          :", ftpSa
1730
                return (resPolySo, infNormSo)                
1731
    # End wile True
1732
    return None
1733
# End pobyso_polynomial_coefficients_progressive_truncate_so_so.
1734

    
1735
def pobyso_polynomial_degree_so_sa(polySo):
1736
    """
1737
    Return the degree of a Sollya polynomial as a Sage int.
1738
    """
1739
    degreeSo = sollya_lib_degree(polySo)
1740
    return pobyso_constant_from_int_so_sa(degreeSo)
1741
# End pobyso_polynomial_degree_so_sa
1742

    
1743
def pobyso_polynomial_degree_so_so(polySo):
1744
    """
1745
    Thin wrapper around lib_sollya_degree().
1746
    """
1747
    return sollya_lib_degree(polySo)
1748
# End pobyso_polynomial_degree_so_so
1749
                                                                  
1750
def pobyso_range(rnLowerBound, rnUpperBound):
1751
    """ Legacy function. See pobyso_range_sa_so. """
1752
    return pobyso_range_sa_so(rnLowerBound, rnUpperBound) 
1753

    
1754
def pobyso_range_from_bounds_sa_so(rnLowerBound, rnUpperBound, precSa = None):
1755
    """
1756
    Create a Sollya range from 2 Sage real numbers as bounds
1757
    """
1758
    # TODO check precision stuff.
1759
    sollyaPrecChanged = False
1760
    (initialSollyaPrecSo, initialSollyaPrecSa) = \
1761
        pobyso_get_prec_so_so_sa()
1762
    if precSa is None:
1763
        precSa = max(rnLowerBound.parent().prec(), rnUpperBound.parent().prec())
1764
    if precSa > initialSollyaPrecSa:
1765
        if precSa <= 2:
1766
            print inspect.stack()[0][3], ": precision change <= 2 requested."
1767
        pobyso_set_prec_sa_so(precSa)
1768
        sollyaPrecChanged = True
1769
    rangeSo = sollya_lib_range_from_bounds(get_rn_value(rnLowerBound),
1770
                                           get_rn_value(rnUpperBound))
1771
    if sollyaPrecChanged:
1772
        pobyso_set_prec_so_so(initialSollyaPrecSo)
1773
        sollya_lib_clear_obj(initialSollyaPrecSo)
1774
    return rangeSo
1775
# End pobyso_range_from_bounds_sa_so
1776

    
1777
def pobyso_range_list_so_sa(listSo):
1778
    """
1779
    Return a Sollya list of ranges as a Sage list of 
1780
    floating-point intervals.
1781
    """
1782
    listSa   = []
1783
    ## The function returns none if the list is empty or an error has happened.
1784
    retVal = pobyso_get_list_elements_so_so(listSo)
1785
    if retVal is None:
1786
        return listSa
1787
    ## Just in case the interface is changed and an empty list is returned 
1788
    #  instead of None.
1789
    elif len(retVal) == 0:
1790
        return listSa
1791
    else:
1792
        ## Remember pobyso_get_list_elements_so_so returns more information
1793
        #  than just the elements of the list (# elements, is_elliptic)
1794
        listSaSo, numElements, isEndElliptic = retVal
1795
    ## Return an empty list.
1796
    if numElements == 0:
1797
        return listSa
1798
    ## Search first for the maximum precision of the elements
1799
    maxPrecSa = 0
1800
    for rangeSo in listSaSo:
1801
        #pobyso_autoprint(floatSo)
1802
        curPrecSa =  pobyso_get_prec_of_range_so_sa(rangeSo)
1803
        if curPrecSa > maxPrecSa:
1804
            maxPrecSa = curPrecSa
1805
    ##
1806
    intervalField = RealIntervalField(maxPrecSa)
1807
    ##
1808
    for rangeSo in listSaSo:
1809
        listSa.append(pobyso_range_to_interval_so_sa(rangeSo, intervalField))
1810
    return listSa
1811
# End pobyso_range_list_so_sa
1812

    
1813
def pobyso_range_max_abs_so_so(rangeSo):
1814
    """
1815
    Return, as Sollya constant, the maximum absolute value of a Sollay range.
1816
    """
1817
    lowerBoundSo = sollya_lib_inf(rangeSo)
1818
    upperBoundSo = sollya_lib_sup(rangeSo)
1819
    #
1820
    lowerBoundSo = sollya_lib_build_function_abs(lowerBoundSo)
1821
    upperBoundSo = sollya_lib_build_function_abs(upperBoundSo)
1822
    #pobyso_autoprint(lowerBoundSo)
1823
    #pobyso_autoprint(upperBoundSo)
1824
    #
1825
    maxAbsSo = sollya_lib_max(lowerBoundSo, upperBoundSo, None)
1826
    #sollya_lib_clear_obj(lowerBoundSo)
1827
    #sollya_lib_clear_obj(upperBoundSo)
1828
    return maxAbsSo
1829
# End pobyso_range_max_abs_so_so
1830
 
1831
def pobyso_range_to_interval_so_sa(rangeSo, realIntervalFieldSa = None):
1832
    """
1833
    Get a Sage interval from a Sollya range.
1834
    If no realIntervalField is given as a parameter, the Sage interval
1835
    precision is that of the Sollya range.
1836
    Otherwise, the precision is that of the realIntervalField. In this case
1837
    rounding may happen.
1838
    """
1839
    if realIntervalFieldSa is None:
1840
        precSa = pobyso_get_prec_of_range_so_sa(rangeSo)
1841
        realIntervalFieldSa = RealIntervalField(precSa)
1842
    intervalSa = \
1843
        pobyso_get_interval_from_range_so_sa(rangeSo, realIntervalFieldSa) 
1844
    return intervalSa
1845
# End pobyso_range_to_interval_so_sa
1846
#
1847
def pobyso_relative_so_so():
1848
    """
1849
    Very thin wrapper around the sollya_lib_relative function.
1850
    """
1851
    return sollya_lib_relative()
1852
# End pobyso_relative_so_so
1853
#
1854
def pobyso_rat_poly_sa_so(polySa, precSa = None):
1855
    """
1856
    Create a Sollya polynomial from a Sage rational polynomial.
1857
    We first convert the rational polynomial into a floating-point 
1858
    polynomial.
1859
    """
1860
    ## TODO: filter arguments.
1861
    ## Precision. If no precision is given, use the current precision
1862
    #  of Sollya.
1863
    if precSa is None:
1864
        precSa =  pobyso_get_prec_so_sa()
1865
    #print "Precision:",  precSa
1866
    RRR = RealField(precSa)
1867
    ## Create a Sage polynomial in the "right" precision.
1868
    P_RRR = RRR[polySa.variables()[0]]
1869
    polyFloatSa = P_RRR(polySa)
1870
    ## Make sure no precision is provided: pobyso_float_poly_sa_so will
1871
    #  recover it all by itself and will not make any extra conversion.
1872
    return pobyso_float_poly_sa_so(polyFloatSa)
1873
    
1874
# End pobyso_rat_poly_sa_so    
1875
    
1876
def pobyso_remez_canonical_sa_sa(func, \
1877
                                 degree, \
1878
                                 lowerBound, \
1879
                                 upperBound, \
1880
                                 weight = None, \
1881
                                 quality = None):
1882
    """
1883
    All arguments are Sage/Python.
1884
    The functions (func and weight) must be passed as expressions or strings.
1885
    Otherwise the function fails. 
1886
    The return value is a Sage polynomial.
1887
    """
1888
    var('zorglub')    # Dummy variable name for type check only. Type of 
1889
    # zorglub is "symbolic expression".
1890
    polySo = pobyso_remez_canonical_sa_so(func, \
1891
                                 degree, \
1892
                                 lowerBound, \
1893
                                 upperBound, \
1894
                                 weight, \
1895
                                 quality)
1896
    # String test
1897
    if parent(func) == parent("string"):
1898
        functionSa = eval(func)
1899
    # Expression test.
1900
    elif type(func) == type(zorglub):
1901
        functionSa = func
1902
    else:
1903
        return None
1904
    #
1905
    maxPrecision = 0
1906
    if polySo is None:
1907
        return(None)
1908
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
1909
    RRRRSa = RealField(maxPrecision)
1910
    polynomialRingSa = RRRRSa[functionSa.variables()[0]]
1911
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, RRRRSa)
1912
    polySa = polynomial(expSa, polynomialRingSa)
1913
    sollya_lib_clear_obj(polySo)
1914
    return(polySa)
1915
# End pobyso_remez_canonical_sa_sa
1916
    
1917
def pobyso_remez_canonical(func, \
1918
                           degree, \
1919
                           lowerBound, \
1920
                           upperBound, \
1921
                           weight = "1", \
1922
                           quality = None):
1923
    """ Legacy function. See pobyso_remez_canonical_sa_so. """
1924
    return(pobyso_remez_canonical_sa_so(func, \
1925
                                        degree, \
1926
                                        lowerBound, \
1927
                                        upperBound, \
1928
                                        weight, \
1929
                                        quality))
1930
# End pobyso_remez_canonical.
1931

    
1932
def pobyso_remez_canonical_sa_so(func, \
1933
                                 degree, \
1934
                                 lowerBound, \
1935
                                 upperBound, \
1936
                                 weight = None, \
1937
                                 quality = None):
1938
    """
1939
    All arguments are Sage/Python.
1940
    The functions (func and weight) must be passed as expressions or strings.
1941
    Otherwise the function fails. 
1942
    The return value is a pointer to a Sollya function.
1943
    lowerBound and upperBound mus be reals.
1944
    """
1945
    var('zorglub')    # Dummy variable name for type check only. Type of
1946
    # zorglub is "symbolic expression".
1947
    currentVariableNameSa = None
1948
    # The func argument can be of different types (string, 
1949
    # symbolic expression...)
1950
    if parent(func) == parent("string"):
1951
        localFuncSa = sage_eval(func,globals())
1952
        if len(localFuncSa.variables()) > 0:
1953
            currentVariableNameSa = localFuncSa.variables()[0]
1954
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1955
            functionSo = \
1956
              sollya_lib_parse_string(localFuncSa._assume_str().replace('_SAGE_VAR_', ''))
1957
    # Expression test.
1958
    elif type(func) == type(zorglub):
1959
        # Until we are able to translate Sage expressions into Sollya 
1960
        # expressions : parse the string version.
1961
        if len(func.variables()) > 0:
1962
            currentVariableNameSa = func.variables()[0]
1963
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1964
            functionSo = \
1965
              sollya_lib_parse_string(func._assume_str().replace('_SAGE_VAR_', ''))
1966
    else:
1967
        return(None)
1968
    if weight is None: # No weight given -> 1.
1969
        weightSo = pobyso_constant_1_sa_so()
1970
    elif parent(weight) == parent("string"): # Weight given as string: parse it.
1971
        weightSo = sollya_lib_parse_string(func)
1972
    elif type(weight) == type(zorglub): # Weight given as symbolice expression.
1973
        functionSo = \
1974
          sollya_lib_parse_string_sa_so(weight._assume_str().replace('_SAGE_VAR_', ''))
1975
    else:
1976
        return(None)
1977
    degreeSo = pobyso_constant_from_int(degree)
1978
    rangeSo = pobyso_bounds_to_range_sa_so(lowerBound, upperBound)
1979
    if not quality is None:
1980
        qualitySo= pobyso_constant_sa_so(quality)
1981
    else:
1982
        qualitySo = None
1983
        
1984
    remezPolySo = sollya_lib_remez(functionSo, \
1985
                                   degreeSo, \
1986
                                   rangeSo, \
1987
                                   weightSo, \
1988
                                   qualitySo, \
1989
                                   None)
1990
    sollya_lib_clear_obj(functionSo)
1991
    sollya_lib_clear_obj(degreeSo)
1992
    sollya_lib_clear_obj(rangeSo)
1993
    sollya_lib_clear_obj(weightSo)
1994
    if not qualitySo is None:
1995
        sollya_lib_clear_obj(qualitySo)
1996
    return(remezPolySo)
1997
# End pobyso_remez_canonical_sa_so
1998

    
1999
def pobyso_remez_canonical_so_so(funcSo, \
2000
                                 degreeSo, \
2001
                                 rangeSo, \
2002
                                 weightSo = pobyso_constant_1_sa_so(),\
2003
                                 qualitySo = None):
2004
    """
2005
    All arguments are pointers to Sollya objects.
2006
    The return value is a pointer to a Sollya function.
2007
    """
2008
    if not sollya_lib_obj_is_function(funcSo):
2009
        return(None)
2010
    return(sollya_lib_remez(funcSo, degreeSo, rangeSo, weightSo, qualitySo, None))
2011
# End pobyso_remez_canonical_so_so.
2012

    
2013
def pobyso_remez_exponents_list_sa_so(func,     \
2014
                                 exponentsList, \
2015
                                 lowerBound,    \
2016
                                 upperBound,    \
2017
                                 weight = None, \
2018
                                 quality = None):
2019
    """
2020
    All arguments are Sage/Python.
2021
    The functions (func and weight) must be passed as expressions or strings.
2022
    Otherwise the function fails. 
2023
    The return value is a pointer to a Sollya function.
2024
    lowerBound and upperBound mus be reals.
2025
    """
2026
    var('zorglub')    # Dummy variable name for type check only. Type of
2027
    # zorglub is "symbolic expression".
2028
    currentVariableNameSa = None
2029
    # The func argument can be of different types (string, 
2030
    # symbolic expression...)
2031
    if parent(func) == parent("string"):
2032
        localFuncSa = sage_eval(func,globals())
2033
        if len(localFuncSa.variables()) > 0:
2034
            currentVariableNameSa = localFuncSa.variables()[0]
2035
            sollya_lib_name_free_variable(str(currentVariableNameSa))
2036
            functionSo = \
2037
              sollya_lib_parse_string(localFuncSa._assume_str().replace('_SAGE_VAR_', ''))
2038
    # Expression test.
2039
    elif type(func) == type(zorglub):
2040
        # Until we are able to translate Sage expressions into Sollya 
2041
        # expressions : parse the string version.
2042
        if len(func.variables()) > 0:
2043
            currentVariableNameSa = func.variables()[0]
2044
            sollya_lib_name_free_variable(str(currentVariableNameSa))
2045
            functionSo = \
2046
              sollya_lib_parse_string(func._assume_str().replace('_SAGE_VAR_', ''))
2047
    else:
2048
        return(None)
2049
    ## Deal with the weight, much in the same way as with the function.
2050
    if weight is None: # No weight given -> 1.
2051
        weightSo = pobyso_constant_1_sa_so()
2052
    elif parent(weight) == parent("string"): # Weight given as string: parse it.
2053
        weightSo = sollya_lib_parse_string(func)
2054
    elif type(weight) == type(zorglub): # Weight given as symbolice expression.
2055
        functionSo = \
2056
          sollya_lib_parse_string(weight._assume_str().replace('_SAGE_VAR_', '',100))
2057
    else:
2058
        return(None)
2059
    rangeSo = pobyso_bounds_to_range_sa_so(lowerBound, upperBound)
2060
    if not quality is None:
2061
        qualitySo= pobyso_constant_sa_so(quality)
2062
    else:
2063
        qualitySo = None
2064
    #
2065
    ## Tranform the Sage list of exponents into a Sollya list.
2066
    exponentsListSo = pobyso_build_list_of_ints_sa_so(*exponentsList)
2067
    remezPolySo = sollya_lib_remez(functionSo,      \
2068
                                   exponentsListSo, \
2069
                                   rangeSo,         \
2070
                                   weightSo,        \
2071
                                   qualitySo,       \
2072
                                   None)
2073
    sollya_lib_clear_obj(functionSo)
2074
    sollya_lib_clear_obj(exponentsListSo)
2075
    sollya_lib_clear_obj(rangeSo)
2076
    sollya_lib_clear_obj(weightSo)
2077
    if not qualitySo is None:
2078
        sollya_lib_clear_obj(qualitySo)
2079
    return(remezPolySo)
2080
# End pobyso_remez_exponentsList_sa_so
2081
#
2082
def pobyso_round_coefficients_so_so(polySo, truncFormatListSo):
2083
    """
2084
    A wrapper around the "classical" sollya_lib_roundcoefficients: a Sollya
2085
    polynomial and a Sollya list are given as arguments.
2086
    """
2087
    return sollya_lib_roundcoefficients(polySo, truncFormatListSo)
2088

    
2089
def pobyso_round_coefficients_progressive_so_so(polySo, 
2090
                                                funcSo,
2091
                                                precSo,
2092
                                                intervalSo,
2093
                                                icSo,
2094
                                                currentApproxErrorSo,
2095
                                                approxAccurSo,
2096
                                                debug=False):
2097
    """
2098
    From an input approximation polynomial, compute an output one with 
2099
    smaller coefficients and yet yields a sufficient approximation accuracy.
2100
    """
2101
    if debug:
2102
        print "Input arguments:"
2103
        print "Polynomial: ", ; pobyso_autoprint(polySo)
2104
        print "Function: ", ; pobyso_autoprint(funcSo)
2105
        print "Internal precision: ", ; pobyso_autoprint(precSo)
2106
        print "Interval: ", ; pobyso_autoprint(intervalSo)
2107
        print "Current approximation error: ", ; pobyso_autoprint(currentApproxErrorSo)
2108
        print "Requested approxiation error: ", ; pobyso_autoprint(approxAccurSo)
2109
        print "________________"
2110
    approxAccurSa        = pobyso_get_constant_as_rn_so_sa(approxAccurSo)
2111
    currentApproxErrorSa = pobyso_get_constant_as_rn_so_sa(currentApproxErrorSo)
2112
    ## If the current approximation error is too close to the target, there is
2113
    #  no possible gain.
2114
    if currentApproxErrorSa >= approxAccurSa / 2:
2115
        #### Do not return the initial argument but copies: caller may free 
2116
        #    the former as inutile after call.
2117
        return (sollya_lib_copy_obj(polySo),
2118
                sollya_lib_copy_obj(currentApproxErrorSo))
2119
    #
2120
    ## Try to round the coefficients.
2121
    degreeSa             = pobyso_polynomial_degree_so_sa(polySo)
2122
    intervalSa           = pobyso_range_to_interval_so_sa(intervalSo)
2123
    
2124
    if debug:
2125
        print "degreeSa              :", degreeSa
2126
        print "intervalSa            :", intervalSa.str(style='brackets')
2127
        print "currentApproxErrorSa  :", currentApproxErrorSa 
2128
        print "approxAccurSa         :", approxAccurSa 
2129
    radiusSa = intervalSa.absolute_diameter() / 2
2130
    if debug:
2131
        print "log2(radius):", RR(radiusSa).log2()
2132
    iterIndex = 0
2133
    ## Build the "shaved" polynomial.
2134
    while True: 
2135
        ### Start with a 0 value expression.
2136
        resPolySo = pobyso_constant_0_sa_so()
2137
        roundedPolyApproxAccurSa = approxAccurSa / 2
2138
        currentRadiusPowerSa = 1 
2139
        for degree in xrange(0,degreeSa + 1):
2140
            #### At round 0, use the agressive formula. At round 1, run the
2141
            #    proved formula.
2142
            if iterIndex == 0:
2143
                roundingPowerSa = \
2144
                    floor(((currentRadiusPowerSa/roundedPolyApproxAccurSa)*(degree+1)).log2())
2145
            else:
2146
                roundingPowerSa = \
2147
                    floor(((currentRadiusPowerSa/roundedPolyApproxAccurSa)*(degreeSa+1)).log2())
2148
            ## Under extreme conditions the above formulas can evaluate under 2,
2149
            #   which is the minimal precision of an MPFR number.
2150
            if roundingPowerSa < 2:
2151
                roundingPowerSa = 2
2152
            if debug:
2153
                print "roundedPolyApproxAccurSa", roundedPolyApproxAccurSa
2154
                print "currentRadiusPowerSa", currentRadiusPowerSa
2155
                print "Current rounding exponent:", roundingPowerSa
2156
            currentRadiusPowerSa *= radiusSa
2157
            index1So = pobyso_constant_from_int_sa_so(degree)
2158
            index2So = pobyso_constant_from_int_sa_so(degree)
2159
            ### Create a monomial with:
2160
            #   - the coefficient in the initial monomial at the current degrree;
2161
            #   - the current exponent;
2162
            #   - the free variable.
2163
            cmonSo  = \
2164
                sollya_lib_build_function_mul(sollya_lib_coeff(polySo, index1So),
2165
                                              sollya_lib_build_function_pow( \
2166
                                                sollya_lib_build_function_free_variable(), \
2167
                                                index2So))
2168
            roundingPowerSo = pobyso_constant_from_int_sa_so(roundingPowerSa)
2169
            cmonrSo = pobyso_round_coefficients_single_so_so(cmonSo, roundingPowerSo)
2170
            sollya_lib_clear_obj(cmonSo)
2171
            ### Add to the result polynomial.
2172
            resPolySo = sollya_lib_build_function_add(resPolySo,
2173
                                                      cmonrSo)
2174
        # End for.
2175
        ### Check the new polynomial.
2176
        freeVarSo     = sollya_lib_build_function_free_variable()
2177
        changeVarSo   = sollya_lib_sub(freeVarSo, icSo)
2178
        resPolyCvSo   = sollya_lib_evaluate(resPolySo, changeVarSo)
2179
        errFuncSo     = sollya_lib_build_function_sub(sollya_lib_copy_obj(funcSo), 
2180
                                                      resPolyCvSo)
2181
        infNormSo = sollya_lib_dirtyinfnorm(errFuncSo, intervalSo)
2182
        ### This also clears resPolyCvSo.
2183
        sollya_lib_clear_obj(errFuncSo)
2184
        cerrSa    = pobyso_get_constant_as_rn_so_sa(infNormSo)
2185
        if debug:
2186
            print "Error of the new polynomial:", cerrSa
2187
        ### If at round 1, return the initial polynomial error. This should
2188
        #   never happen since the rounding algorithm is proved. But some 
2189
        #   circumstances may break it (e.g. internal precision of tools).
2190
        if cerrSa > approxAccurSa:
2191
            if iterIndex == 0: # Round 0 is agressive rounding, got round 1 (proved rounding)
2192
                sollya_lib_clear_obj(resPolySo)
2193
                sollya_lib_clear_obj(infNormSo)
2194
                iterIndex += 1
2195
                continue
2196
            else: # Round 1 and beyond : just return the oroginal polynomial.
2197
                sollya_lib_clear_obj(resPolySo)
2198
                sollya_lib_clear_obj(infNormSo)
2199
                #### Do not return the arguments but copies: the caller may free
2200
                #    free the former as inutile after call.
2201
                return (sollya_lib_copy_obj(polySo), 
2202
                        sollya_lib_copy_obj(currentApproxErrorSo))
2203
        ### If get here it is because cerrSa <= approxAccurSa
2204
        ### Approximation error of the new polynomial is acceptable.
2205
        return (resPolySo, infNormSo)
2206
    # End while True
2207
# End pobyso_round_coefficients_progressive_so_so
2208
    
2209
def pobyso_round_coefficients_single_so_so(polySo, commonPrecSo):
2210
    """
2211
    Create a rounded coefficients polynomial from polynomial argument to
2212
    the number of bits in size argument.
2213
    All coefficients are set to the same precision.
2214
    """
2215
    ## TODO: check arguments.
2216
    endEllipListSo = pobyso_build_end_elliptic_list_so_so(commonPrecSo)
2217
    polySo = sollya_lib_roundcoefficients(polySo, endEllipListSo, None)
2218
    sollya_lib_clear_obj(endEllipListSo)
2219
    #sollya_lib_clear_obj(endEllipListSo)
2220
    return polySo
2221
    
2222
# End pobyso_round_coefficients_single_so_so
2223

    
2224
def pobyso_set_canonical_off():
2225
    sollya_lib_set_canonical(sollya_lib_off())
2226

    
2227
def pobyso_set_canonical_on():
2228
    sollya_lib_set_canonical(sollya_lib_on())
2229

    
2230
def pobyso_set_prec(p):
2231
    """ Legacy function. See pobyso_set_prec_sa_so. """
2232
    pobyso_set_prec_sa_so(p)
2233

    
2234
def pobyso_set_prec_sa_so(p):
2235
    #a = c_int(p)
2236
    #precSo = c_void_p(sollya_lib_constant_from_int(a))
2237
    #precSo = sollya_lib_constant_from_int(a)
2238
    precSo = pobyso_constant_from_int_sa_so(p)
2239
    sollya_lib_set_prec(precSo)
2240
    sollya_lib_clear_obj(precSo)
2241
# End pobyso_set_prec_sa_so.
2242

    
2243
def pobyso_set_prec_so_so(newPrecSo):
2244
    sollya_lib_set_prec(newPrecSo)
2245
# End pobyso_set_prec_so_so.
2246
#   
2247
def pobyso_supnorm_sa_sa(polySa, 
2248
                         funcSa, 
2249
                         intervalSa, 
2250
                         errorTypeSa=SOLLYA_ABSOLUTE, 
2251
                         accuracySa=RR(2**-40)):
2252
    """
2253
    An supnorm call with Sage arguments.
2254
    """
2255
    # Check that polySa is a univariate polynomial. We only check here at
2256
    # expression level, not at polynomial ring level
2257
    ## Make sure it is not a multivariate polynomial. The number of variables
2258
    #  can be zero in the case of a constant.
2259
    try:
2260
        if len(polySa.free_variables()) > 1 :
2261
            print "Invalid number of free variables in polySa:", polySa,
2262
            " -", polySa.free_variables()
2263
            return None
2264
    except AttributeError:
2265
        print "polySa is not a SymbolicExpression."
2266
        return None 
2267
    ## Make sure it is a polynomial.
2268
    if not polySa.is_polynomial(polySa.default_variable()):
2269
        print "polySa is not a polynomila expression."
2270
        return None
2271
    # Check that funcSa is a function.
2272
    if not \
2273
     sage.symbolic.callable.is_CallableSymbolicExpressionRing(parent(funcSa)):
2274
        return None
2275
    # Check that intervalSa is an interval.
2276
    try:
2277
        intervalSa.upper()
2278
    except AttributeError:
2279
        print "intervalSa is not an interval."
2280
        return None
2281
    # Convert the Sage polynomial into a Sollya polynomial.
2282
    polyAsStringSa = polySa._assume_str().replace('_SAGE_VAR_', '')
2283
    polySo = pobyso_parse_string_sa_so(polyAsStringSa)
2284
    if not pobyso_obj_is_function_so_sa(polySo):
2285
        sollya_lib_clear_obj(polySo)
2286
        print "The Sollya object created from the polynomial is not a function."
2287
        return None
2288
    #pobyso_autoprint(polySo)
2289
    # Convert the Sage function into a Sollya function.
2290
    funcAsStringSa = funcSa._assume_str().replace('_SAGE_VAR_', '')
2291
    funcSo = pobyso_parse_string_sa_so(funcAsStringSa)
2292
    if not pobyso_obj_is_function_so_sa(funcSo):
2293
        sollya_lib_clear_obj(polySo)
2294
        sollya_lib_clear_obj(funcSo)
2295
        print "The Sollya object created from the function is not a function."
2296
        return None
2297
    #pobyso_autoprint(funcSo)
2298
    # Convert the Sage interval into a Sollya range.
2299
    rangeSo = pobyso_interval_to_range_sa_so(intervalSa)
2300
    if not pobyso_is_range_so_sa(rangeSo):
2301
        sollya_lib_clear_obj(polySo)
2302
        sollya_lib_clear_obj(funcSo)
2303
        sollya_lib_clear_obj(rangeSo)
2304
        print "The Sollya object created from the interval is not a range."
2305
        return None
2306
    #pobyso_autoprint(rangeSo)
2307
    # Check the error type. We do not check the returned object since we
2308
    # trust our code and Sollya on this one.
2309
    if errorTypeSa != SOLLYA_ABSOLUTE and errorTypeSa != SOLLYA_RELATIVE:
2310
        sollya_lib_clear_obj(polySo)
2311
        sollya_lib_clear_obj(funcSo)
2312
        sollya_lib_clear_obj(rangeSo)
2313
        return None
2314
    if errorTypeSa == SOLLYA_ABSOLUTE:
2315
        errorTypeSo = pobyso_absolute_so_so()
2316
    else:
2317
        errorTypeSo = pobyso_relative_so_so()
2318
    #pobyso_autoprint(errorTypeSo)
2319
    # Check if accuracySa is an element of a RealField. If not, try to 
2320
    # convert it.
2321
    try:
2322
        accuracySa.ulp()
2323
    except AttributeError:
2324
        try:
2325
            accuracySa = RR(accuracySa)
2326
        except TypeError:
2327
            accuracySa = RR(str(accuracySa))
2328
    # Create the Sollya object
2329
    accuracySo = pobyso_constant_sa_so(accuracySa)
2330
    #pobyso_autoprint(accuracySo)
2331
    if pobyso_is_error_so_sa(accuracySo):
2332
        sollya_lib_clear_obj(polySo)
2333
        sollya_lib_clear_obj(funcSo)
2334
        sollya_lib_clear_obj(rangeSo)
2335
        sollya_lib_clear_obj(errorTypeSo)
2336
        sollya_lib_clear_obj(accuracySo)
2337
        return None
2338
    retValSo = sollya_lib_supnorm(polySo, 
2339
                                  funcSo, 
2340
                                  rangeSo, 
2341
                                  errorTypeSo, 
2342
                                  accuracySo, 
2343
                                  None)
2344
    sollya_lib_clear_obj(polySo)
2345
    sollya_lib_clear_obj(funcSo)
2346
    sollya_lib_clear_obj(rangeSo)
2347
    sollya_lib_clear_obj(errorTypeSo)
2348
    sollya_lib_clear_obj(accuracySo)
2349
    if pobyso_is_error_so_sa(retValSo):
2350
        sollya_lib_clear_obj(retValSo)
2351
        return None
2352
    #pobyso_autoprint(retValSo)
2353
    retValSa = pobyso_range_to_interval_so_sa(retValSo)
2354
    sollya_lib_clear_obj(retValSo)
2355
    return retValSa
2356
# End pobyso_supnorm_sa_sa
2357

    
2358
def pobyso_supnorm_so_sa(polySo, funcSo, intervalSo, errorTypeSo = None,\
2359
                         accuracySo = None, realFieldSa = None):
2360
    """
2361
    Computes the supremum norm from Sollya input arguments and returns a
2362
    Sage floating-point number whose precision is set by the only Sage argument.
2363
    
2364
    The returned value is the maximum of the absolute values of the range
2365
    elements  returned by the Sollya supnorm functions
2366
    """
2367
    supNormRangeSo = pobyso_supnorm_so_so(polySo, 
2368
                                          funcSo, 
2369
                                          intervalSo, 
2370
                                          errorTypeSo,
2371
                                          accuracySo)
2372
    supNormSo = pobyso_range_max_abs_so_so(supNormRangeSo)
2373
    sollya_lib_clear_obj(supNormRangeSo)
2374
    #pobyso_autoprint(supNormSo)
2375
    supNormSa = pobyso_get_constant_as_rn_with_rf_so_sa(supNormSo, realFieldSa)
2376
    sollya_lib_clear_obj(supNormSo)
2377
    return supNormSa 
2378
# End pobyso_supnorm_so_sa.
2379
#
2380
def pobyso_supnorm_so_so(polySo, funcSo, intervalSo, errorTypeSo = None,\
2381
                         accuracySo = None):
2382
    """
2383
    Computes the supnorm of the approximation error between the given 
2384
    polynomial and function. Attention: returns a range!
2385
    errorTypeSo defaults to "absolute".
2386
    accuracySo defaults to 2^(-40).
2387
    """
2388
    if errorTypeSo is None:
2389
        errorTypeSo = sollya_lib_absolute(None)
2390
        errorTypeIsNone = True
2391
    else:
2392
        errorTypeIsNone = False
2393
    #
2394
    if accuracySo is None:
2395
        # Notice the **: we are in Pythonland here!
2396
        accuracySo = pobyso_constant_sa_so(RR(2**(-40)))
2397
        accuracyIsNone = True
2398
    else:
2399
        accuracyIsNone = False
2400
    #pobyso_autoprint(accuracySo)
2401
    resultSo = \
2402
        sollya_lib_supnorm(polySo, funcSo, intervalSo, errorTypeSo, \
2403
                              accuracySo)
2404
    if errorTypeIsNone:
2405
        sollya_lib_clear_obj(errorTypeSo)
2406
    if accuracyIsNone:
2407
        sollya_lib_clear_obj(accuracySo)
2408
    return resultSo
2409
# End pobyso_supnorm_so_so
2410
#
2411
def pobyso_taylor_expansion_no_change_var_so_so(functionSo, 
2412
                                                degreeSo, 
2413
                                                rangeSo,
2414
                                                errorTypeSo=None, 
2415
                                                sollyaPrecSo=None):
2416
    """
2417
    Compute the Taylor expansion without the variable change
2418
    x -> x-intervalCenter.
2419
    If errorTypeSo is None, absolute is used.
2420
    If sollyaPrecSo is None, Sollya internal precision is not changed. 
2421
    """
2422
    # Change internal Sollya precision, if needed.
2423
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
2424
    sollyaPrecChanged = False
2425
    if sollyaPrecSo is None:
2426
        pass
2427
    else:
2428
        sollya_lib_set_prec(sollyaPrecSo)
2429
        sollyaPrecChanged = True
2430
    # Error type stuff: default to absolute.
2431
    if errorTypeSo is None:
2432
        errorTypeIsNone = True
2433
        errorTypeSo = sollya_lib_absolute(None)
2434
    else:
2435
        errorTypeIsNone = False
2436
    intervalCenterSo = sollya_lib_mid(rangeSo, None)
2437
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo,
2438
                                         intervalCenterSo,
2439
                                         rangeSo, errorTypeSo, None)
2440
    # Object taylorFormListSaSo is a Python list of Sollya objects references 
2441
    #  that are copies of the elements of taylorFormSo.
2442
    # pobyso_get_list_elements_so_so clears taylorFormSo.
2443
    (taylorFormListSaSo, numElementsSa, isEndEllipticSa) = \
2444
        pobyso_get_list_elements_so_so(taylorFormSo)
2445
    ## Copy needed here since polySo will be returned and taylorFormListSaSo
2446
    #  will be cleared.
2447
    polySo = sollya_lib_copy_obj(taylorFormListSaSo[0])
2448
    #print "Num elements:", numElementsSa
2449
    sollya_lib_clear_obj(taylorFormSo)
2450
    # No copy_obj needed here: a new objects are created.
2451
    maxErrorSo    = sollya_lib_sup(taylorFormListSaSo[2])
2452
    minErrorSo    = sollya_lib_inf(taylorFormListSaSo[2])
2453
    # List taylorFormListSaSo is not needed anymore.
2454
    pobyso_clear_list_elements_sa_so(taylorFormListSaSo)
2455
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
2456
    absMinErrorSo = sollya_lib_abs(minErrorSo)
2457
    sollya_lib_clear_obj(maxErrorSo)
2458
    sollya_lib_clear_obj(minErrorSo)
2459
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
2460
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
2461
    #
2462
    if errorTypeIsNone:
2463
        sollya_lib_clear_obj(errorTypeSo)
2464
    ## If changed, reset the Sollya working precision.
2465
    if sollyaPrecChanged:
2466
        sollya_lib_set_prec(initialSollyaPrecSo)
2467
    sollya_lib_clear_obj(initialSollyaPrecSo)
2468
    ## According to what error is the largest, return the errors.
2469
    if absMaxErrorSa > absMinErrorSa:
2470
        sollya_lib_clear_obj(absMinErrorSo)
2471
        return (polySo, intervalCenterSo, absMaxErrorSo)
2472
    else:
2473
        sollya_lib_clear_obj(absMaxErrorSo)
2474
        return (polySo, intervalCenterSo, absMinErrorSo)
2475
# end pobyso_taylor_expansion_no_change_var_so_so
2476

    
2477
def pobyso_taylor_expansion_with_change_var_so_so(functionSo, degreeSo, \
2478
                                                  rangeSo, \
2479
                                                  errorTypeSo=None, \
2480
                                                  sollyaPrecSo=None):
2481
    """
2482
    Compute the Taylor expansion with the variable change
2483
    x -> (x-intervalCenter) included.
2484
    """
2485
    # Change Sollya internal precision, if need.
2486
    sollyaPrecChanged = False
2487
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
2488
    if sollyaPrecSo is None:
2489
        pass
2490
    else:
2491
        sollya_lib_set_prec(sollyaPrecSo)
2492
        sollyaPrecChanged = True
2493
    #
2494
    # Error type stuff: default to absolute.
2495
    if errorTypeSo is None:
2496
        errorTypeIsNone = True
2497
        errorTypeSo = sollya_lib_absolute(None)
2498
    else:
2499
        errorTypeIsNone = False
2500
    intervalCenterSo = sollya_lib_mid(rangeSo)
2501
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo, \
2502
                                         intervalCenterSo, \
2503
                                         rangeSo, errorTypeSo, None)
2504
    # Object taylorFormListSaSo is a Python list of Sollya objects references 
2505
    # that are copies of the elements of taylorFormSo.
2506
    # pobyso_get_list_elements_so_so clears taylorFormSo.
2507
    (taylorFormListSaSo, numElements, isEndElliptic) = \
2508
        pobyso_get_list_elements_so_so(taylorFormSo)
2509
    sollya_lib_clear_obj(taylorFormSo)
2510
    polySo = taylorFormListSaSo[0]
2511
    ## Maximum error computation with taylorFormListSaSo[2], a range
2512
    #  holding the actual error. Bounds can be negative.
2513
    maxErrorSo    = sollya_lib_sup(taylorFormListSaSo[2])
2514
    minErrorSo    = sollya_lib_inf(taylorFormListSaSo[2])
2515
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
2516
    absMinErrorSo = sollya_lib_abs(minErrorSo)
2517
    sollya_lib_clear_obj(maxErrorSo)
2518
    sollya_lib_clear_obj(minErrorSo)
2519
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
2520
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
2521
    changeVarExpSo = sollya_lib_build_function_sub(\
2522
                       sollya_lib_build_function_free_variable(),\
2523
                       sollya_lib_copy_obj(intervalCenterSo))
2524
    polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpSo)
2525
    # List taylorFormListSaSo is not needed anymore.
2526
    pobyso_clear_list_elements_sa_so(taylorFormListSaSo)
2527
    sollya_lib_clear_obj(changeVarExpSo)
2528
    # If changed, reset the Sollya working precision.
2529
    if sollyaPrecChanged:
2530
        sollya_lib_set_prec(initialSollyaPrecSo)
2531
    sollya_lib_clear_obj(initialSollyaPrecSo)
2532
    if errorTypeIsNone:
2533
        sollya_lib_clear_obj(errorTypeSo)
2534
    # Do not clear maxErrorSo.
2535
    if absMaxErrorSa > absMinErrorSa:
2536
        sollya_lib_clear_obj(absMinErrorSo)
2537
        return((polyVarChangedSo, intervalCenterSo, absMaxErrorSo))
2538
    else:
2539
        sollya_lib_clear_obj(absMaxErrorSo)
2540
        return((polyVarChangedSo, intervalCenterSo, absMinErrorSo))
2541
# end pobyso_taylor_expansion_with_change_var_so_so
2542

    
2543
def pobyso_taylor(function, degree, point):
2544
    """ Legacy function. See pobysoTaylor_so_so. """
2545
    return(pobyso_taylor_so_so(function, degree, point))
2546

    
2547
def pobyso_taylor_so_so(functionSo, degreeSo, pointSo):
2548
    return(sollya_lib_taylor(functionSo, degreeSo, pointSo))
2549
    
2550
def pobyso_taylorform(function, degree, point = None, 
2551
                      interval = None, errorType=None):
2552
    """ Legacy function. See pobyso_taylorform_sa_sa;"""
2553
    
2554
def pobyso_taylorform_sa_sa(functionSa, \
2555
                            degreeSa, \
2556
                            pointSa, \
2557
                            intervalSa=None, \
2558
                            errorTypeSa=None, \
2559
                            precisionSa=None):
2560
    """
2561
    Compute the Taylor form of 'degreeSa' for 'functionSa' at 'pointSa' 
2562
    for 'intervalSa' with 'errorTypeSa' (a string) using 'precisionSa'. 
2563
    point: must be a Real or a Real interval.
2564
    return the Taylor form as an array
2565
    TODO: take care of the interval and of the point when it is an interval;
2566
          when errorType is not None;
2567
          take care of the other elements of the Taylor form (coefficients 
2568
          errors and delta.
2569
    """
2570
    # Absolute as the default error.
2571
    if errorTypeSa is None:
2572
        errorTypeSo = sollya_lib_absolute()
2573
    elif errorTypeSa == "relative":
2574
        errorTypeSo = sollya_lib_relative()
2575
    elif errortypeSa == "absolute":
2576
        errorTypeSo = sollya_lib_absolute()
2577
    else:
2578
        # No clean up needed.
2579
        return None
2580
    # Global precision stuff
2581
    sollyaPrecisionChangedSa = False
2582
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
2583
    if precisionSa is None:
2584
        precSa = initialSollyaPrecSa
2585
    else:
2586
        if precSa > initialSollyaPrecSa:
2587
            if precSa <= 2:
2588
                print inspect.stack()[0][3], ":precision change <= 2 requested."
2589
            pobyso_set_prec_sa_so(precSa)
2590
            sollyaPrecisionChangedSa = True
2591
    #        
2592
    if len(functionSa.variables()) > 0:
2593
        varSa = functionSa.variables()[0]
2594
        pobyso_name_free_variable_sa_so(str(varSa))
2595
    # In any case (point or interval) the parent of pointSa has a precision
2596
    # method.
2597
    pointPrecSa = pointSa.parent().precision()
2598
    if precSa > pointPrecSa:
2599
        pointPrecSa = precSa
2600
    # In any case (point or interval) pointSa has a base_ring() method.
2601
    pointBaseRingString = str(pointSa.base_ring())
2602
    if re.search('Interval', pointBaseRingString) is None: # Point
2603
        pointSo = pobyso_constant_sa_so(pointSa, pointPrecSa)
2604
    else: # Interval.
2605
        pointSo = pobyso_interval_to_range_sa_so(pointSa, pointPrecSa)
2606
    # Sollyafy the function.
2607
    functionSo = pobyso_parse_string_sa_so(functionSa._assume_str().replace('_SAGE_VAR_', ''))
2608
    if sollya_lib_obj_is_error(functionSo):
2609
        print "pobyso_tailorform: function string can't be parsed!"
2610
        return None
2611
    # Sollyafy the degree
2612
    degreeSo = sollya_lib_constant_from_int(int(degreeSa))
2613
    # Sollyafy the point
2614
    # Call Sollya
2615
    taylorFormSo = \
2616
        sollya_lib_taylorform(functionSo, degreeSo, pointSo, errorTypeSo,\
2617
                                         None)
2618
    sollya_lib_clear_obj(functionSo)
2619
    sollya_lib_clear_obj(degreeSo)
2620
    sollya_lib_clear_obj(pointSo)
2621
    sollya_lib_clear_obj(errorTypeSo)
2622
    (tfsAsList, numElements, isEndElliptic) = \
2623
            pobyso_get_list_elements_so_so(taylorFormSo)
2624
    polySo = tfsAsList[0]
2625
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
2626
    polyRealField = RealField(maxPrecision)
2627
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, polyRealField)
2628
    if sollyaPrecisionChangedSa:
2629
        sollya_lib_set_prec(initialSollyaPrecSo)
2630
    sollya_lib_clear_obj(initialSollyaPrecSo)
2631
    polynomialRing = polyRealField[str(varSa)]
2632
    polySa = polynomial(expSa, polynomialRing)
2633
    taylorFormSa = [polySa]
2634
    # Final clean-up.
2635
    sollya_lib_clear_obj(taylorFormSo)
2636
    return(taylorFormSa)
2637
# End pobyso_taylor_form_sa_sa
2638

    
2639
def pobyso_taylorform_so_so(functionSo, degreeSo, pointSo, intervalSo=None, \
2640
                            errorTypeSo=None):
2641
    createdErrorType = False
2642
    if errorTypeSo is None:
2643
        errorTypeSo = sollya_lib_absolute()
2644
        createdErrorType = True
2645
    else:
2646
        #TODO: deal with the other case.
2647
        pass
2648
    if intervalSo is None:
2649
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
2650
                                         errorTypeSo, None)
2651
    else:
2652
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
2653
                                         intervalSo, errorTypeSo, None)
2654
    if createdErrorType:
2655
        sollya_lib_clear_obj(errorTypeSo)
2656
    return resultSo
2657
        
2658

    
2659
def pobyso_univar_polynomial_print_reverse(polySa):
2660
    """ Legacy function. See pobyso_univar_polynomial_print_reverse_sa_sa. """
2661
    return(pobyso_univar_polynomial_print_reverse_sa_sa(polySa))
2662

    
2663
def pobyso_univar_polynomial_print_reverse_sa_sa(polySa):
2664
    """
2665
    Return the string representation of a univariate polynomial with
2666
    monomials ordered in the x^0..x^n order of the monomials.
2667
    Remember: Sage
2668
    """
2669
    polynomialRing = polySa.base_ring()
2670
    # A very expensive solution:
2671
    # -create a fake multivariate polynomial field with only one variable,
2672
    #   specifying a negative lexicographical order;
2673
    mpolynomialRing = PolynomialRing(polynomialRing.base(), \
2674
                                     polynomialRing.variable_name(), \
2675
                                     1, order='neglex')
2676
    # - convert the univariate argument polynomial into a multivariate
2677
    #   version;
2678
    p = mpolynomialRing(polySa)
2679
    # - return the string representation of the converted form.
2680
    # There is no simple str() method defined for p's class.
2681
    return(p.__str__())
2682
#
2683
#print pobyso_get_prec()  
2684
pobyso_set_prec(165)
2685
#print pobyso_get_prec()  
2686
#a=100
2687
#print type(a)
2688
#id(a)
2689
#print "Max arity: ", pobyso_max_arity
2690
#print "Function tripleDouble (43) as a string: ", pobyso_function_type_as_string(43)
2691
#print "Function None (44) as a string: ", pobyso_function_type_as_string(44)
2692
sys.stderr.write("\t...Pobyso check done.\n")