Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 285

Historique | Voir | Annoter | Télécharger (102,42 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
    """
1457
    try:
1458
        coefficientsListSa = objectSa
1459
        if len(objectSa.coefficients()) < 0:
1460
            return False
1461
    except AttributeError:
1462
        return False
1463
    return True
1464
# End pobyso_is_poly_sa_sa
1465

    
1466
def pobyso_is_range_so_sa(rangeCandidateSo):
1467
    """
1468
    Thin wrapper over sollya_lib_is_range.
1469
    """
1470
    return sollya_lib_obj_is_range(rangeCandidateSo) != 0
1471
    
1472
# End pobyso_is_range_so_sa
1473

    
1474

    
1475
def pobyso_lib_init():
1476
    sollya_lib_init(None)
1477

    
1478
def pobyso_lib_close():
1479
    sollya_lib_close(None)
1480

    
1481
def pobyso_name_free_variable(freeVariableNameSa):
1482
    """ Legacy function. See pobyso_name_free_variable_sa_so. """
1483
    pobyso_name_free_variable_sa_so(freeVariableNameSa)
1484

    
1485
def pobyso_name_free_variable_sa_so(freeVariableNameSa):
1486
    """
1487
    Set the free variable name in Sollya from a Sage string.
1488
    """
1489
    sollya_lib_name_free_variable(freeVariableNameSa)
1490

    
1491
def pobyso_obj_is_function_so_sa(objSo):
1492
    """
1493
    Check if an object is a function.
1494
    """
1495
    if sollya_lib_obj_is_function(objSo) != 0:
1496
        return True
1497
    else:
1498
        return False
1499
# End pobyso_obj_is_function_so_sa  
1500
    
1501
def pobyso_obj_is_range_so_sa(objSo):
1502
    """
1503
    Check if an object is a function.
1504
    """
1505
    if sollya_lib_obj_is_range(objSo) != 0:
1506
        return True
1507
    else:
1508
        return False
1509
# End pobyso_obj_is_range_so_sa  
1510
    
1511
def pobyso_obj_is_string_so_sa(objSo):
1512
    """
1513
    Check if an object is a function.
1514
    """
1515
    if sollya_lib_obj_is_string(objSo) != 0:
1516
        return True
1517
    else:
1518
        return False
1519
# End pobyso_obj_is_string_so_sa  
1520
    
1521
def pobyso_parse_string(string):
1522
    """ Legacy function. See pobyso_parse_string_sa_so. """
1523
    return pobyso_parse_string_sa_so(string)
1524
 
1525
def pobyso_parse_string_sa_so(string):
1526
    """
1527
    Get the Sollya expression computed from a Sage string or
1528
    a Sollya error object if parsing failed.
1529
    """
1530
    return sollya_lib_parse_string(string)
1531

    
1532
def pobyso_precision_so_sa(ctExpSo):
1533
    """
1534
    Computes the necessary precision to represent a number.
1535
    If x is not zero, it can be uniquely written as x = m · 2e 
1536
    where m is an odd integer and e is an integer. 
1537
    precision(x) returns the number of bits necessary to write m 
1538
    in binary (i.e. ceil(log2(m))).
1539
    """
1540
    #TODO: take care of the special case: 0, @NaN@, @Inf@
1541
    precisionSo = sollya_lib_precision(ctExpSo)
1542
    precisionSa = pobyso_constant_from_int_so_sa(precisionSo)
1543
    sollya_lib_clear_obj(precisionSo)
1544
    return precisionSa
1545
# End pobyso_precision_so_sa
1546

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

    
1713
def pobyso_polynomial_degree_so_sa(polySo):
1714
    """
1715
    Return the degree of a Sollya polynomial as a Sage int.
1716
    """
1717
    degreeSo = sollya_lib_degree(polySo)
1718
    return pobyso_constant_from_int_so_sa(degreeSo)
1719
# End pobyso_polynomial_degree_so_sa
1720

    
1721
def pobyso_polynomial_degree_so_so(polySo):
1722
    """
1723
    Thin wrapper around lib_sollya_degree().
1724
    """
1725
    return sollya_lib_degree(polySo)
1726
# End pobyso_polynomial_degree_so_so
1727
                                                                  
1728
def pobyso_range(rnLowerBound, rnUpperBound):
1729
    """ Legacy function. See pobyso_range_sa_so. """
1730
    return pobyso_range_sa_so(rnLowerBound, rnUpperBound) 
1731

    
1732
def pobyso_range_from_bounds_sa_so(rnLowerBound, rnUpperBound, precSa = None):
1733
    """
1734
    Create a Sollya range from 2 Sage real numbers as bounds
1735
    """
1736
    # TODO check precision stuff.
1737
    sollyaPrecChanged = False
1738
    (initialSollyaPrecSo, initialSollyaPrecSa) = \
1739
        pobyso_get_prec_so_so_sa()
1740
    if precSa is None:
1741
        precSa = max(rnLowerBound.parent().prec(), rnUpperBound.parent().prec())
1742
    if precSa > initialSollyaPrecSa:
1743
        if precSa <= 2:
1744
            print inspect.stack()[0][3], ": precision change <= 2 requested."
1745
        pobyso_set_prec_sa_so(precSa)
1746
        sollyaPrecChanged = True
1747
    rangeSo = sollya_lib_range_from_bounds(get_rn_value(rnLowerBound),
1748
                                           get_rn_value(rnUpperBound))
1749
    if sollyaPrecChanged:
1750
        pobyso_set_prec_so_so(initialSollyaPrecSo)
1751
        sollya_lib_clear_obj(initialSollyaPrecSo)
1752
    return rangeSo
1753
# End pobyso_range_from_bounds_sa_so
1754

    
1755
def pobyso_range_list_so_sa(listSo):
1756
    """
1757
    Return a Sollya list of ranges as a Sage list of 
1758
    floating-point intervals.
1759
    """
1760
    listSa   = []
1761
    ## The function returns none if the list is empty or an error has happened.
1762
    retVal = pobyso_get_list_elements_so_so(listSo)
1763
    if retVal is None:
1764
        return listSa
1765
    ## Just in case the interface is changed and an empty list is returned 
1766
    #  instead of None.
1767
    elif len(retVal) == 0:
1768
        return listSa
1769
    else:
1770
        ## Remember pobyso_get_list_elements_so_so returns more information
1771
        #  than just the elements of the list (# elements, is_elliptic)
1772
        listSaSo, numElements, isEndElliptic = retVal
1773
    ## Return an empty list.
1774
    if numElements == 0:
1775
        return listSa
1776
    ## Search first for the maximum precision of the elements
1777
    maxPrecSa = 0
1778
    for rangeSo in listSaSo:
1779
        #pobyso_autoprint(floatSo)
1780
        curPrecSa =  pobyso_get_prec_of_range_so_sa(rangeSo)
1781
        if curPrecSa > maxPrecSa:
1782
            maxPrecSa = curPrecSa
1783
    ##
1784
    intervalField = RealIntervalField(maxPrecSa)
1785
    ##
1786
    for rangeSo in listSaSo:
1787
        listSa.append(pobyso_range_to_interval_so_sa(rangeSo, intervalField))
1788
    return listSa
1789
# End pobyso_range_list_so_sa
1790

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

    
1910
def pobyso_remez_canonical_sa_so(func, \
1911
                                 degree, \
1912
                                 lowerBound, \
1913
                                 upperBound, \
1914
                                 weight = None, \
1915
                                 quality = None):
1916
    """
1917
    All arguments are Sage/Python.
1918
    The functions (func and weight) must be passed as expressions or strings.
1919
    Otherwise the function fails. 
1920
    The return value is a pointer to a Sollya function.
1921
    lowerBound and upperBound mus be reals.
1922
    """
1923
    var('zorglub')    # Dummy variable name for type check only. Type of
1924
    # zorglub is "symbolic expression".
1925
    currentVariableNameSa = None
1926
    # The func argument can be of different types (string, 
1927
    # symbolic expression...)
1928
    if parent(func) == parent("string"):
1929
        localFuncSa = sage_eval(func,globals())
1930
        if len(localFuncSa.variables()) > 0:
1931
            currentVariableNameSa = localFuncSa.variables()[0]
1932
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1933
            functionSo = \
1934
              sollya_lib_parse_string(localFuncSa._assume_str().replace('_SAGE_VAR_', ''))
1935
    # Expression test.
1936
    elif type(func) == type(zorglub):
1937
        # Until we are able to translate Sage expressions into Sollya 
1938
        # expressions : parse the string version.
1939
        if len(func.variables()) > 0:
1940
            currentVariableNameSa = func.variables()[0]
1941
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1942
            functionSo = \
1943
              sollya_lib_parse_string(func._assume_str().replace('_SAGE_VAR_', ''))
1944
    else:
1945
        return(None)
1946
    if weight is None: # No weight given -> 1.
1947
        weightSo = pobyso_constant_1_sa_so()
1948
    elif parent(weight) == parent("string"): # Weight given as string: parse it.
1949
        weightSo = sollya_lib_parse_string(func)
1950
    elif type(weight) == type(zorglub): # Weight given as symbolice expression.
1951
        functionSo = \
1952
          sollya_lib_parse_string_sa_so(weight._assume_str().replace('_SAGE_VAR_', ''))
1953
    else:
1954
        return(None)
1955
    degreeSo = pobyso_constant_from_int(degree)
1956
    rangeSo = pobyso_bounds_to_range_sa_so(lowerBound, upperBound)
1957
    if not quality is None:
1958
        qualitySo= pobyso_constant_sa_so(quality)
1959
    else:
1960
        qualitySo = None
1961
        
1962
    remezPolySo = sollya_lib_remez(functionSo, \
1963
                                   degreeSo, \
1964
                                   rangeSo, \
1965
                                   weightSo, \
1966
                                   qualitySo, \
1967
                                   None)
1968
    sollya_lib_clear_obj(functionSo)
1969
    sollya_lib_clear_obj(degreeSo)
1970
    sollya_lib_clear_obj(rangeSo)
1971
    sollya_lib_clear_obj(weightSo)
1972
    if not qualitySo is None:
1973
        sollya_lib_clear_obj(qualitySo)
1974
    return(remezPolySo)
1975
# End pobyso_remez_canonical_sa_so
1976

    
1977
def pobyso_remez_canonical_so_so(funcSo, \
1978
                                 degreeSo, \
1979
                                 rangeSo, \
1980
                                 weightSo = pobyso_constant_1_sa_so(),\
1981
                                 qualitySo = None):
1982
    """
1983
    All arguments are pointers to Sollya objects.
1984
    The return value is a pointer to a Sollya function.
1985
    """
1986
    if not sollya_lib_obj_is_function(funcSo):
1987
        return(None)
1988
    return(sollya_lib_remez(funcSo, degreeSo, rangeSo, weightSo, qualitySo, None))
1989
# End pobyso_remez_canonical_so_so.
1990

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

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

    
2202
def pobyso_set_canonical_off():
2203
    sollya_lib_set_canonical(sollya_lib_off())
2204

    
2205
def pobyso_set_canonical_on():
2206
    sollya_lib_set_canonical(sollya_lib_on())
2207

    
2208
def pobyso_set_prec(p):
2209
    """ Legacy function. See pobyso_set_prec_sa_so. """
2210
    pobyso_set_prec_sa_so(p)
2211

    
2212
def pobyso_set_prec_sa_so(p):
2213
    #a = c_int(p)
2214
    #precSo = c_void_p(sollya_lib_constant_from_int(a))
2215
    #precSo = sollya_lib_constant_from_int(a)
2216
    precSo = pobyso_constant_from_int_sa_so(p)
2217
    sollya_lib_set_prec(precSo)
2218
    sollya_lib_clear_obj(precSo)
2219
# End pobyso_set_prec_sa_so.
2220

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

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

    
2455
def pobyso_taylor_expansion_with_change_var_so_so(functionSo, degreeSo, \
2456
                                                  rangeSo, \
2457
                                                  errorTypeSo=None, \
2458
                                                  sollyaPrecSo=None):
2459
    """
2460
    Compute the Taylor expansion with the variable change
2461
    x -> (x-intervalCenter) included.
2462
    """
2463
    # Change Sollya internal precision, if need.
2464
    sollyaPrecChanged = False
2465
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
2466
    if sollyaPrecSo is None:
2467
        pass
2468
    else:
2469
        sollya_lib_set_prec(sollyaPrecSo)
2470
        sollyaPrecChanged = True
2471
    #
2472
    # Error type stuff: default to absolute.
2473
    if errorTypeSo is None:
2474
        errorTypeIsNone = True
2475
        errorTypeSo = sollya_lib_absolute(None)
2476
    else:
2477
        errorTypeIsNone = False
2478
    intervalCenterSo = sollya_lib_mid(rangeSo)
2479
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo, \
2480
                                         intervalCenterSo, \
2481
                                         rangeSo, errorTypeSo, None)
2482
    # Object taylorFormListSaSo is a Python list of Sollya objects references 
2483
    # that are copies of the elements of taylorFormSo.
2484
    # pobyso_get_list_elements_so_so clears taylorFormSo.
2485
    (taylorFormListSaSo, numElements, isEndElliptic) = \
2486
        pobyso_get_list_elements_so_so(taylorFormSo)
2487
    sollya_lib_clear_obj(taylorFormSo)
2488
    polySo = taylorFormListSaSo[0]
2489
    ## Maximum error computation with taylorFormListSaSo[2], a range
2490
    #  holding the actual error. Bounds can be negative.
2491
    maxErrorSo    = sollya_lib_sup(taylorFormListSaSo[2])
2492
    minErrorSo    = sollya_lib_inf(taylorFormListSaSo[2])
2493
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
2494
    absMinErrorSo = sollya_lib_abs(minErrorSo)
2495
    sollya_lib_clear_obj(maxErrorSo)
2496
    sollya_lib_clear_obj(minErrorSo)
2497
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
2498
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
2499
    changeVarExpSo = sollya_lib_build_function_sub(\
2500
                       sollya_lib_build_function_free_variable(),\
2501
                       sollya_lib_copy_obj(intervalCenterSo))
2502
    polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpSo)
2503
    # List taylorFormListSaSo is not needed anymore.
2504
    pobyso_clear_list_elements_sa_so(taylorFormListSaSo)
2505
    sollya_lib_clear_obj(changeVarExpSo)
2506
    # If changed, reset the Sollya working precision.
2507
    if sollyaPrecChanged:
2508
        sollya_lib_set_prec(initialSollyaPrecSo)
2509
    sollya_lib_clear_obj(initialSollyaPrecSo)
2510
    if errorTypeIsNone:
2511
        sollya_lib_clear_obj(errorTypeSo)
2512
    # Do not clear maxErrorSo.
2513
    if absMaxErrorSa > absMinErrorSa:
2514
        sollya_lib_clear_obj(absMinErrorSo)
2515
        return((polyVarChangedSo, intervalCenterSo, absMaxErrorSo))
2516
    else:
2517
        sollya_lib_clear_obj(absMaxErrorSo)
2518
        return((polyVarChangedSo, intervalCenterSo, absMinErrorSo))
2519
# end pobyso_taylor_expansion_with_change_var_so_so
2520

    
2521
def pobyso_taylor(function, degree, point):
2522
    """ Legacy function. See pobysoTaylor_so_so. """
2523
    return(pobyso_taylor_so_so(function, degree, point))
2524

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

    
2617
def pobyso_taylorform_so_so(functionSo, degreeSo, pointSo, intervalSo=None, \
2618
                            errorTypeSo=None):
2619
    createdErrorType = False
2620
    if errorTypeSo is None:
2621
        errorTypeSo = sollya_lib_absolute()
2622
        createdErrorType = True
2623
    else:
2624
        #TODO: deal with the other case.
2625
        pass
2626
    if intervalSo is None:
2627
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
2628
                                         errorTypeSo, None)
2629
    else:
2630
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
2631
                                         intervalSo, errorTypeSo, None)
2632
    if createdErrorType:
2633
        sollya_lib_clear_obj(errorTypeSo)
2634
    return resultSo
2635
        
2636

    
2637
def pobyso_univar_polynomial_print_reverse(polySa):
2638
    """ Legacy function. See pobyso_univar_polynomial_print_reverse_sa_sa. """
2639
    return(pobyso_univar_polynomial_print_reverse_sa_sa(polySa))
2640

    
2641
def pobyso_univar_polynomial_print_reverse_sa_sa(polySa):
2642
    """
2643
    Return the string representation of a univariate polynomial with
2644
    monomials ordered in the x^0..x^n order of the monomials.
2645
    Remember: Sage
2646
    """
2647
    polynomialRing = polySa.base_ring()
2648
    # A very expensive solution:
2649
    # -create a fake multivariate polynomial field with only one variable,
2650
    #   specifying a negative lexicographical order;
2651
    mpolynomialRing = PolynomialRing(polynomialRing.base(), \
2652
                                     polynomialRing.variable_name(), \
2653
                                     1, order='neglex')
2654
    # - convert the univariate argument polynomial into a multivariate
2655
    #   version;
2656
    p = mpolynomialRing(polySa)
2657
    # - return the string representation of the converted form.
2658
    # There is no simple str() method defined for p's class.
2659
    return(p.__str__())
2660
#
2661
#print pobyso_get_prec()  
2662
pobyso_set_prec(165)
2663
#print pobyso_get_prec()  
2664
#a=100
2665
#print type(a)
2666
#id(a)
2667
#print "Max arity: ", pobyso_max_arity
2668
#print "Function tripleDouble (43) as a string: ", pobyso_function_type_as_string(43)
2669
#print "Function None (44) as a string: ", pobyso_function_type_as_string(44)
2670
sys.stderr.write("\t...Pobyso check done.\n")