Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 284

Historique | Voir | Annoter | Télécharger (98,23 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
    try:
125
        interval = theRIF(lowerBound, upperBound)
126
    except TypeError:
127
        return None
128
    else:
129
        return interval
130
# End pobyso_bounds_to_interval_sa_sa
131
    
132
def pobyso_bounds_to_range_sa_so(rnLowerBoundSa, rnUpperBoundSa, \
133
                                 precisionSa=None):
134
    """
135
    Return a Sollya range from to 2 RealField Sage elements.
136
    The Sollya range element has a sufficient precision to hold all
137
    the digits of the widest of the Sage bounds.
138
    """
139
    # Sanity check.
140
    if rnLowerBoundSa > rnUpperBoundSa:
141
        return None
142
    # Precision stuff.
143
    if precisionSa is None:
144
        # Check for the largest precision.
145
        lbPrecSa = rnLowerBoundSa.parent().precision()
146
        ubPrecSa = rnLowerBoundSa.parent().precision()
147
        maxPrecSa = max(lbPrecSa, ubPrecSa)
148
    else:
149
        maxPrecSa = precisionSa
150
    # From Sage to Sollya bounds.
151
#    lowerBoundSo = sollya_lib_constant(get_rn_value(rnLowerBoundSa),
152
#                                       maxPrecSa)
153
    lowerBoundSo = pobyso_constant_sa_so(rnLowerBoundSa,
154
                                         maxPrecSa)
155
    upperBoundSo = pobyso_constant_sa_so(rnUpperBoundSa,
156
                                       maxPrecSa)
157
    
158
    # From Sollya bounds to range.
159
    rangeSo = sollya_lib_range(lowerBoundSo, upperBoundSo)
160
    # Back to original precision.
161
    # Clean up
162
    sollya_lib_clear_obj(lowerBoundSo)
163
    sollya_lib_clear_obj(upperBoundSo)
164
    return rangeSo
165
# End pobyso_bounds_to_range_sa_so
166

    
167
def pobyso_build_end_elliptic_list_so_so(*args):
168
    """
169
    From Sollya argument objects, create a Sollya end elliptic list.
170
    Elements of the list are "eaten" (should not be cleared individually, 
171
    are cleared when the list is cleared).
172
    """
173
    if len(args) == 0:
174
        ## When called with an empty list, sollya_lib_build_end_elliptic_list,
175
        #  produces "error".
176
        return sollya_lib_build_list(None)
177
    index = 0
178
    ## One can not append elements to an elliptic list, prepend only is 
179
    #  permitted.
180
    for argument in reversed(args):
181
        if index == 0:
182
            listSo = sollya_lib_build_end_elliptic_list(argument, None)
183
        else:
184
            listSo = sollya_lib_prepend(argument, listSo)
185
        index += 1
186
    return listSo
187
        
188
# End pobyso_build_end_elliptic_list_so_so
189

    
190
def pobyso_build_function_sub_so_so(exp1So, exp2So):
191
    return sollya_lib_build_function_sub(exp1So, exp2So)
192

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

    
228
def pobyso_change_var_in_function_so_so(funcSo, chvarExpSo):
229
    """
230
    Variable change in a function.
231
    """
232
    return(sollya_lib_evaluate(funcSo,chvarExpSo))
233
# End pobyso_change_var_in_function_so_so     
234

    
235
def pobyso_chebyshevform_so_so(functionSo, degreeSo, intervalSo):
236
    resultSo = sollya_lib_chebyshevform(functionSo, degreeSo, intervalSo)
237
    return(resultSo)
238
# End pobyso_chebyshevform_so_so.
239

    
240
def pobyso_clear_full_list_elements_sa_so(objectListSaSo):
241
    """
242
    Clear the elements of list created by the 
243
    pobyso_get_list_elements_so_so function.
244
    objectListSaSo is as follows:
245
    - objectListSaSo[0]: a list of Sollya objects;
246
    - objectListSaSo[1]: the number of elements of the previous list;
247
    - objectListSaSo[2]: an integer that if != 0 states that the list is 
248
                         end-elliptic
249
    The objects to clear are the elements of the objectListSaSo[0] list. 
250
    """
251
    for index in xrange(0, objectListSaSo[1]):
252
        sollya_lib_clear_obj(objectListSaSo[0][index])
253
# End pobyso_clear_full_list_elements_sa_so
254

    
255
def pobyso_clear_list_elements_sa_so(objectListSaSo):
256
    """
257
    Clear the elements of list of references to Sollya objects
258
    """
259
    for index in xrange(0, len(objectListSaSo)):
260
        sollya_lib_clear_obj(objectListSaSo[index])
261
# End pobyso_clear_list_elements_sa_so
262

    
263
def pobyso_clear_obj(objSo):
264
    """
265
    Free a Sollya object's memory.
266
    Very thin wrapper around sollya_lib_clear_obj().
267
    """
268
    sollya_lib_clear_obj(objSo)
269
# End pobyso_clear_obj. 
270

    
271
def pobyso_clear_taylorform_sa_so(taylorFormSaSo):
272
    """
273
    This method is rapper around pobyso_clear_list_elements_sa_so.
274
    It is a legacy method left here since it may be used in existing code 
275
    where Taylor forms are used as Sage lists obtained by converting
276
    Sollya Taylor forms (a list made of:
277
    - a polynomial;
278
    - a list of intervals enclosing the errors accumulated when computing
279
      the polynomial coefficients;
280
    - a bound on the approximation error between the polynomial and the 
281
      function.)
282
    A Taylor form directly obtained from pobyso_taylorform_so_so is cleared
283
    by sollya_lib_clear_obj. 
284
    """
285
    pobyso_clear_list_elements_sa_so(taylorFormSaSo)
286
# End pobyso_clear_taylorform_sa_so 
287

    
288
def pobyso_cmp(rnArgSa, cteSo):
289
    """
290
    Deprecated, use pobyso_cmp_sa_so_sa instead.
291
    """
292
    print "Deprecated, use pobyso_cmp_sa_so_sa instead."
293
    return pobyso_cmp_sa_so_sa(rnArgSa, cteSo)
294
# End  pobyso_cmp
295
    
296
def pobyso_cmp_sa_so_sa(rnArgSa, cteSo):
297
    """
298
    Compare the MPFR value a RealNumber with that of a Sollya constant.
299
    
300
    Get the value of the Sollya constant into a RealNumber and compare
301
    using MPFR. Could be optimized by working directly with a mpfr_t
302
    for the intermediate number. 
303
    """
304
    # Get the precision of the Sollya constant to build a Sage RealNumber 
305
    # with enough precision.to hold it.
306
    precisionOfCte = c_int(0)
307
    # From the Sollya constant, create a local Sage RealNumber.
308
    sollya_lib_get_prec_of_constant(precisionOfCte, cteSo) 
309
    #print "Precision of constant: ", precisionOfCte
310
    RRRR = RealField(precisionOfCte.value)
311
    rnLocalSa = RRRR(0)
312
    sollya_lib_get_constant(get_rn_value(rnLocalSa), cteSo)
313
    #
314
    ## Compare the Sage RealNumber version of the Sollya constant with rnArg
315
    #  through a direct comparison of underlying MPFR numbers.
316
    return cmp_rn_value(rnArgSa, rnLocal)
317
# End pobyso_smp_sa_so_sa
318

    
319
def pobyso_compute_pos_function_abs_val_bounds_sa_sa(funcSa, lowerBoundSa, \
320
                                                     upperBoundSa):
321
    """
322
    TODO: completely rework and test.
323
    """
324
    pobyso = pobyso_name_free_variable_sa_so(funcSa.variables()[0])
325
    funcSo = pobyso_parse_string(funcSa._assume_str().replace('_SAGE_VAR_', ''))
326
    rangeSo = pobyso_range_sa_so(lowerBoundSa, upperBoundSa)
327
    infnormSo = pobyso_infnorm_so_so(funcSo,rangeSo)
328
    # Sollya return the infnorm as an interval.
329
    fMaxSa = pobyso_get_interval_from_range_so_sa(infnormSo)
330
    # Get the top bound and compute the binade top limit.
331
    fMaxUpperBoundSa = fMaxSa.upper()
332
    binadeTopLimitSa = 2**ceil(fMaxUpperBoundSa.log2())
333
    # Put up together the function to use to compute the lower bound.
334
    funcAuxSo = pobyso_parse_string(str(binadeTopLimitSa) +  \
335
                                    '-(' + f._assume_str().replace('_SAGE_VAR_', '') + ')')
336
    pobyso_autoprint(funcAuxSo)
337
    # Clear the Sollya range before a new call to infnorm and issue the call.
338
    sollya_lib_clear_obj(infnormSo)
339
    infnormSo = pobyso_infnorm_so_so(funcAuxSo,rangeSo)
340
    fMinSa = pobyso_get_interval_from_range_so_sa(infnormSo)
341
    sollya_lib_clear_obj(infnormSo)
342
    fMinLowerBoundSa = binadeTopLimitSa - fMinSa.lower()
343
    # Compute the maximum of the precisions of the different bounds.
344
    maxPrecSa = max([fMinLowerBoundSa.parent().precision(), \
345
                     fMaxUpperBoundSa.parent().precision()])
346
    # Create a RealIntervalField and create an interval with the "good" bounds.
347
    RRRI = RealIntervalField(maxPrecSa)
348
    imageIntervalSa = RRRI(fMinLowerBoundSa, fMaxUpperBoundSa)
349
    # Free the unneeded Sollya objects
350
    sollya_lib_clear_obj(funcSo)
351
    sollya_lib_clear_obj(funcAuxSo)
352
    sollya_lib_clear_obj(rangeSo)
353
    return(imageIntervalSa)
354
# End pobyso_compute_pos_function_abs_val_bounds_sa_sa
355

    
356
def pobyso_compute_precision_decay_ratio_function_sa_so():
357
    """
358
    Compute the precision decay ratio function for polynomial 
359
    coefficient progressive trucation.
360
    """
361
    functionText = """
362
    proc(deg, a, b, we, wq) 
363
    {
364
      k = we * (exp(x/a)-1) + wq * (b*x)^2 + (1-we-wq) * x;
365
      return k/k(d);
366
    };
367
    """
368
    return pobyso_parse_string_sa_so(functionText)
369
# End  pobyso_compute_precision_decay_ratio_function.
370

    
371

    
372
def pobyso_constant(rnArg):
373
    """ Legacy function. See pobyso_constant_sa_so. """
374
    return(pobyso_constant_sa_so(rnArg))
375
    
376
def pobyso_constant_sa_so(rnArgSa, precisionSa=None):
377
    """
378
    Create a Sollya constant from a Sage RealNumber.
379
    The sollya_lib_constant() function creates a constant
380
    with the same precision as the source.
381
    """
382
    ## Precision stuff. If one wants to change precisions,
383
    #  everything takes place in Sage. This only makes
384
    #  sense if one wants to reduce the precision.
385
    # TODO: revisit precision stuff with new technique.
386
    # Check that rnArgSa is a realFiedl element.
387
    try:
388
        rnArgSa.ulp()
389
    except AttributeError:
390
        return pobyso_error_so()
391
    # If a different precision is wanted modify it.
392
    if not precisionSa is None:
393
        RRR = RealField(precisionSa)
394
        rnArgSa = RRR(rnArgSa)
395
    #print rnArgSa, rnArgSa.precision()
396
    # Sollya constant creation takes place here.
397
    return sollya_lib_constant(get_rn_value(rnArgSa))
398
# End pobyso_constant_sa_so
399
     
400
def pobyso_constant_0_sa_so():
401
    """
402
    Obvious.
403
    """
404
    return pobyso_constant_from_int_sa_so(0)
405

    
406
def pobyso_constant_1():
407
    """
408
    Obvious.
409
    Legacy function. See pobyso_constant_so_so. 
410
    """
411
    return pobyso_constant_1_sa_so()
412

    
413
def pobyso_constant_1_sa_so():
414
    """
415
    Obvious.
416
    """
417
    return(pobyso_constant_from_int_sa_so(1))
418

    
419
def pobyso_constant_from_int(anInt):
420
    """ Legacy function. See pobyso_constant_from_int_sa_so. """
421
    return pobyso_constant_from_int_sa_so(anInt)
422

    
423
def pobyso_constant_from_int_sa_so(anInt):
424
    """
425
    Get a Sollya constant from a Sage int.
426
    """
427
    return sollya_lib_constant_from_int64(long(anInt))
428

    
429
def pobyso_constant_from_int_so_sa(constSo):
430
    """
431
    Get a Sage int from a Sollya int constant.
432
    Usefull for precision or powers in polynomials.
433
    """
434
    constSa = c_long(0)
435
    sollya_lib_get_constant_as_int64(byref(constSa), constSo)
436
    return constSa.value
437
# End pobyso_constant_from_int_so_sa
438

    
439
def pobyso_constant_from_mpq_sa_so(rationalSa):
440
    """
441
    Make a Sollya constant from Sage rational.
442
    The Sollya constant is an unevaluated expression.
443
    Hence no precision argument is needed.
444
    It is better to leave this way since Sollya has its own
445
    optimized evaluation mecanism that tries very hard to
446
    return exact values or at least faithful ones.
447
    """
448
    ratExprSo = \
449
        sollya_lib_constant_from_mpq(sgmp_get_rational_value(rationalSa))
450
    return ratExprSo
451
# End pobyso_constant_from_mpq_sa_so.
452

    
453
def pobyso_constant_sollya_prec_sa_so(rnArgSa):
454
    """
455
    Create a Sollya constant from a Sage RealNumber at the
456
    current precision in Sollya.
457
    """
458
    currentSollyaPrecSa = pobyso_get_prec_so_sa()
459
    return pobyso_constant_sa_so(rnArgSa, currentSollyaPrecSa)
460
# End pobyso_constant_sollya_prec_sa_so
461

    
462
def pobyso_end_elliptic_list_so_sa_so(objectsListSo, intCountSa):
463
    """
464
    Create a Sollya end elliptic list made of the objectListSo[0] to
465
     objectsListSo[intCountSa-1] objects.
466
    """     
467
    return sollya_lib_end_elliptic_list(objectSo, int(intCountSa))
468

    
469
def pobyso_error_so():
470
    return sollya_lib_error(None)
471
# End pobyso_error().
472

    
473
def pobyso_evaluate_so_sa(funcSo, argumentSo):
474
    """
475
    Evaluates funcSo for arguemntSo through sollya_lib_evaluate() and return
476
    the result as a Sage object
477
    """
478
    evalSo = sollya_lib_evaluate(funcSo, argumentSo)
479
    if pobyso_is_error_so_sa(evalSo):
480
        return None
481
    if pobyso_is_range_so_sa(evalSo):
482
        retVal = pobyso_range_to_interval_so_sa(evalSo)
483
    else:
484
        retVal = pobyso_get_constant_as_rn(evalSo)
485
    sollya_lib_clear_obj(evalSo)
486
    return retVal
487
# End pobyso_evaluate_so_sa.
488

    
489
def pobyso_evaluate_so_so(funcSo, argumentSo):
490
    """
491
    Evaluates funcSo for arguemntSo through sollya_lib_evaluate().
492
    """
493
    return sollya_lib_evaluate(funcSo, argumentSo)
494
# End pobyso_evaluate_so_so.
495
#
496
def pobyso_diff_so_so(funcSo):
497
    """
498
    Very thin wrapper around sollya_lib_diff.
499
    """
500
    ## TODO: add a check to make sure funcSo is a functional expression.
501
    return sollya_lib_diff(funcSo)
502

    
503
def pobyso_dirty_find_zeros_so_so(funcSo, rangeSo):
504
    """
505
    Thin wrapper over sollya_lib_dirtyfindzeros()
506
    """
507
    return sollya_lib_dirtyfindzeros(funcSo, rangeSo)
508
# End pobys_dirty_find_zeros
509

    
510
def pobyso_dirty_inf_norm_so_so(funcSo, rangeSo, preSo=None):
511
    """
512
    Thin wrapper around sollya_dirtyinfnorm().
513
    """
514
    # TODO: manage the precision change.
515
    
516
    return sollya_lib_dirtyinfnorm(funcSo, rangeSo)
517
# End pobyso_dirty_inf_norm_so_so
518

    
519
def pobyso_find_zeros_so_so(funcSo, rangeSo):
520
    """
521
    Thin wrapper over sollya_lib_findzeros()
522
    """
523
    return sollya_lib_findzeros(funcSo, rangeSo)
524
# End pobys_find_zeros
525

    
526
def pobyso_float_list_so_sa(listSo):
527
    """
528
    Return a Sollya list of floating-point numbers as a Sage list of 
529
    floating-point numbers.
530
    TODO: add a test to make sure that each element of the list is a constant.
531
    """
532
    listSa   = []
533
    ## The function returns none if the list is empty or an error has happened.
534
    retVal = pobyso_get_list_elements_so_so(listSo)
535
    if retVal is None:
536
        return listSa
537
    ## Just in case the interface is changed and an empty list is returned 
538
    #  instead of None.
539
    elif len(retVal) == 0:
540
        return listSa
541
    else:
542
        ## Remember pobyso_get_list_elements_so_so returns more information
543
        #  than just the elements of the list (# elements, is_elliptic)
544
        listSaSo, numElements, isEndElliptic = retVal
545
    ## Return an empty list.
546
    if numElements == 0:
547
        return listSa
548
    ## Search first for the maximum precision of the elements
549
    maxPrecSa = 0
550
    for floatSo in listSaSo:
551
        #pobyso_autoprint(floatSo)
552
        curPrecSa =  pobyso_get_prec_of_constant_so_sa(floatSo)
553
        if curPrecSa > maxPrecSa:
554
            maxPrecSa = curPrecSa
555
    ##
556
    RF = RealField(maxPrecSa)
557
    ##
558
    for floatSo in listSaSo:
559
        listSa.append(pobyso_get_constant_as_rn_with_rf_so_sa(floatSo))
560
    return listSa
561
# End pobyso_float_list_so_sa
562

    
563
def pobyso_float_poly_sa_so(polySa, precSa = None):
564
    """
565
    Create a Sollya polynomial from a Sage RealField polynomial.
566
    """
567
    ## TODO: filter arguments.
568
    ## Precision. If a precision is given, convert the polynomial
569
    #  into the right polynomial field. If not convert it straight
570
    #  to Sollya.
571
    sollyaPrecChanged = False 
572
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
573
    if precSa is None:
574
        precSa = polySa.parent().base_ring().precision()
575
    if (precSa > initialSollyaPrecSa):
576
        assert precSa >= 2, "Precision change <2 requested"
577
        if precSa <= 2:
578
            print inspect.stack()[0][3], ": precision change <= 2 requested"
579
        precSo = pobyso_constant_from_int(precSa)
580
        pobyso_set_prec_so_so(precSo)
581
        sollya_lib_clear_obj(precSo)
582
        sollyaPrecChanged = True
583
    ## Free variable stuff.
584
    freeVariableNameChanged = False 
585
    polyFreeVariableNameSa = \
586
        str(polySa.variables()[0])
587
    currentFreeVariableNameSa = pobyso_get_free_variable_name_so_sa() 
588
    if polyFreeVariableNameSa != currentFreeVariableNameSa:
589
        #print "Free variable names do not match.", polyFreeVariableNameSa
590
        sollya_lib_name_free_variable(polyFreeVariableNameSa)
591
        freeVariableNameChanged = True
592
    ## Get exponents and coefficients.
593
    exponentsSa     = polySa.exponents()
594
    coefficientsSa  = polySa.coefficients()
595
    ## Build the polynomial.
596
    polySo = None
597
    for coefficientSa, exponentSa in zip(coefficientsSa, exponentsSa):
598
        #print coefficientSa.n(prec=precSa), exponentSa
599
        coefficientSo = \
600
            pobyso_constant_sa_so(coefficientSa)
601
        #pobyso_autoprint(coefficientSo)
602
        exponentSo = \
603
            pobyso_constant_from_int_sa_so(exponentSa)
604
        #pobyso_autoprint(exponentSo)
605
        monomialSo = sollya_lib_build_function_pow(
606
                       sollya_lib_build_function_free_variable(),
607
                       exponentSo)
608
        polyTermSo = sollya_lib_build_function_mul(coefficientSo,
609
                                                       monomialSo)
610
        if polySo is None:
611
            polySo = polyTermSo
612
        else:
613
            polySo = sollya_lib_build_function_add(polySo, polyTermSo)
614
    if sollyaPrecChanged:
615
        pobyso_set_prec_so_so(initialSollyaPrecSo)
616
    sollya_lib_clear_obj(initialSollyaPrecSo)
617
    ## Do not set back the free variable name in Sollya to its initial value: 
618
    #  it will change it back, in the Sollya polynomial, to what it was in the 
619
    #  first place.
620
    return polySo
621
# End pobyso_float_poly_sa_so    
622

    
623
def pobyso_float_poly_so_sa(polySo, realFieldSa=None):
624
    """
625
    Convert a Sollya polynomial into a Sage floating-point polynomial.
626
    If no realField is given, a RealField corresponding to the maximum 
627
    precision of the coefficients is internally computed.
628
    The real field is not returned but can be easily retrieved from 
629
    the polynomial itself.
630
    ALGORITHM:
631
    - (optional) compute the RealField of the coefficients;
632
    - convert the Sollya expression into a Sage expression;
633
    - convert the Sage expression into a Sage polynomial
634
    """    
635
    if realFieldSa is None:
636
        expressionPrecSa = pobyso_get_max_prec_of_exp_so_sa(polySo)
637
        #print "Maximum precision of Sollya polynomial coefficients:", expressionPrecSa
638
        if expressionPrecSa < 2 or expressionPrecSa > 2147483391:
639
            print "Maximum degree of expression:", expressionPrecSa
640
        realFieldSa      = RealField(expressionPrecSa)
641
    #print "Sollya expression before...",
642
    #pobyso_autoprint(polySo)
643

    
644
    expressionSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo,
645
                                                             realFieldSa)
646
    #print "...Sollya expression after."
647
    #pobyso_autoprint(polySo)
648
    polyVariableSa = expressionSa.variables()[0]
649
    polyRingSa     = realFieldSa[str(polyVariableSa)]
650
    #print polyRingSa
651
    # Do not use the polynomial(expressionSa, ring=polyRingSa) form!
652
    polynomialSa = polyRingSa(expressionSa)
653
    polyCoeffsListSa = polynomialSa.coefficients()
654
    #for coeff in polyCoeffsListSa:
655
    #    print coeff.abs().n()
656
    return polynomialSa
657
# End pobyso_float_poly_so_sa
658

    
659
def pobyso_free_variable():
660
    """
661
    Ultra thin wrapper around the sollya_lib_function_build_free_variable function.
662
    """
663
    return sollya_lib_build_function_free_variable()
664
    
665
def pobyso_function_type_as_string(funcType):
666
    """ Legacy function. See pobyso_function_type_as_string_so_sa. """
667
    return(pobyso_function_type_as_string_so_sa(funcType))
668

    
669
def pobyso_function_type_as_string_so_sa(funcType):
670
    """
671
    Numeric Sollya function codes -> Sage mathematical function names.
672
    Notice that pow -> ^ (a la Sage, not a la Python).
673
    """
674
    if funcType == SOLLYA_BASE_FUNC_ABS:
675
        return "abs"
676
    elif funcType == SOLLYA_BASE_FUNC_ACOS:
677
        return "arccos"
678
    elif funcType == SOLLYA_BASE_FUNC_ACOSH:
679
        return "arccosh"
680
    elif funcType == SOLLYA_BASE_FUNC_ADD:
681
        return "+"
682
    elif funcType == SOLLYA_BASE_FUNC_ASIN:
683
        return "arcsin"
684
    elif funcType == SOLLYA_BASE_FUNC_ASINH:
685
        return "arcsinh"
686
    elif funcType == SOLLYA_BASE_FUNC_ATAN:
687
        return "arctan"
688
    elif funcType == SOLLYA_BASE_FUNC_ATANH:
689
        return "arctanh"
690
    elif funcType == SOLLYA_BASE_FUNC_CEIL:
691
        return "ceil"
692
    elif funcType == SOLLYA_BASE_FUNC_CONSTANT:
693
        return "cte"
694
    elif funcType == SOLLYA_BASE_FUNC_COS:
695
        return "cos"
696
    elif funcType == SOLLYA_BASE_FUNC_COSH:
697
        return "cosh"
698
    elif funcType == SOLLYA_BASE_FUNC_DIV:
699
        return "/"
700
    elif funcType == SOLLYA_BASE_FUNC_DOUBLE:
701
        return "double"
702
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEDOUBLE:
703
        return "doubleDouble"
704
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEEXTENDED:
705
        return "doubleDxtended"
706
    elif funcType == SOLLYA_BASE_FUNC_ERF:
707
        return "erf"
708
    elif funcType == SOLLYA_BASE_FUNC_ERFC:
709
        return "erfc"
710
    elif funcType == SOLLYA_BASE_FUNC_EXP:
711
        return "exp"
712
    elif funcType == SOLLYA_BASE_FUNC_EXP_M1:
713
        return "expm1"
714
    elif funcType == SOLLYA_BASE_FUNC_FLOOR:
715
        return "floor"
716
    elif funcType == SOLLYA_BASE_FUNC_FREE_VARIABLE:
717
        return "freeVariable"
718
    elif funcType == SOLLYA_BASE_FUNC_HALFPRECISION:
719
        return "halfPrecision"
720
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYCONSTANT:
721
        return "libraryConstant"
722
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYFUNCTION:
723
        return "libraryFunction"
724
    elif funcType == SOLLYA_BASE_FUNC_LOG:
725
        return "log"
726
    elif funcType == SOLLYA_BASE_FUNC_LOG_10:
727
        return "log10"
728
    elif funcType == SOLLYA_BASE_FUNC_LOG_1P:
729
        return "log1p"
730
    elif funcType == SOLLYA_BASE_FUNC_LOG_2:
731
        return "log2"
732
    elif funcType == SOLLYA_BASE_FUNC_MUL:
733
        return "*"
734
    elif funcType == SOLLYA_BASE_FUNC_NEARESTINT:
735
        return "round"
736
    elif funcType == SOLLYA_BASE_FUNC_NEG:
737
        return "__neg__"
738
    elif funcType == SOLLYA_BASE_FUNC_PI:
739
        return "pi"
740
    elif funcType == SOLLYA_BASE_FUNC_POW:
741
        return "^"
742
    elif funcType == SOLLYA_BASE_FUNC_PROCEDUREFUNCTION:
743
        return "procedureFunction"
744
    elif funcType == SOLLYA_BASE_FUNC_QUAD:
745
        return "quad"
746
    elif funcType == SOLLYA_BASE_FUNC_SIN:
747
        return "sin"
748
    elif funcType == SOLLYA_BASE_FUNC_SINGLE:
749
        return "single"
750
    elif funcType == SOLLYA_BASE_FUNC_SINH:
751
        return "sinh"
752
    elif funcType == SOLLYA_BASE_FUNC_SQRT:
753
        return "sqrt"
754
    elif funcType == SOLLYA_BASE_FUNC_SUB:
755
        return "-"
756
    elif funcType == SOLLYA_BASE_FUNC_TAN:
757
        return "tan"
758
    elif funcType == SOLLYA_BASE_FUNC_TANH:
759
        return "tanh"
760
    elif funcType == SOLLYA_BASE_FUNC_TRIPLEDOUBLE:
761
        return "tripleDouble"
762
    else:
763
        return None
764

    
765
def pobyso_get_constant(rnArgSa, constSo):
766
    """ Legacy function. See pobyso_get_constant_so_sa. """
767
    return pobyso_get_constant_so_sa(rnArgSa, constSo)
768
# End pobyso_get_constant
769

    
770
def pobyso_get_constant_so_sa(rnArgSa, constSo):
771
    """
772
    Set the value of rnArgSa to the value of constSo in MPFR_RNDN mode.
773
    rnArg must already exist and belong to some RealField.
774
    We assume that constSo points to a Sollya constant.
775
    """
776
    outcome = sollya_lib_get_constant(get_rn_value(rnArgSa), constSo)
777
    if outcome == 0: # Failure because constSo is not a constant expression.
778
        return None
779
    else:
780
        return outcome
781
# End  pobyso_get_constant_so_sa
782
   
783
def pobyso_get_constant_as_rn(ctExpSo):
784
    """ 
785
    Legacy function. See pobyso_get_constant_as_rn_so_sa. 
786
    """ 
787
    return(pobyso_get_constant_as_rn_so_sa(ctExpSo))
788
    
789
def pobyso_get_constant_as_rn_so_sa(constExpSo):
790
    """
791
    Get a Sollya constant as a Sage "real number".
792
    The precision of the floating-point number returned is that of the Sollya
793
    constant.
794
    """
795
    #print "Before computing precision of variable..."
796
    #pobyso_autoprint(constExpSo)
797
    precisionSa  = pobyso_get_prec_of_constant_so_sa(constExpSo)
798
    #print "precisionSa:", precisionSa
799
    ## If the expression can not be exactly converted, None is returned.
800
    #  In this case opt for the Sollya current expression.
801
    if precisionSa is None:
802
        precisionSa = pobyso_get_prec_so_sa()
803
    RRRR = RealField(precisionSa)
804
    rnSa = RRRR(0)
805
    outcome = sollya_lib_get_constant(get_rn_value(rnSa), constExpSo)
806
    if outcome == 0:
807
        return None
808
    else:
809
        return rnSa
810
# End pobyso_get_constant_as_rn_so_sa
811

    
812
def pobyso_get_constant_as_rn_with_rf(ctExp, realField):
813
    """ 
814
    Legacy function. See pobyso_get_constant_as_rn_with_rf_so_sa.
815
    """
816
    return pobyso_get_constant_as_rn_with_rf_so_sa(ctExp, realField)
817
# End pobyso_get_constant_as_rn_with_rf
818
    
819
def pobyso_get_constant_as_rn_with_rf_so_sa(ctExpSo, realFieldSa = None):
820
    """
821
    Get a Sollya constant as a Sage "real number".
822
    If no real field is specified, the precision of the floating-point number 
823
    returned is that of the Sollya constant.
824
    Otherwise is is that of the real field. Hence rounding may happen.
825
    """
826
    if realFieldSa is None:
827
        return pobyso_get_constant_as_rn_so_sa(ctExpSo)
828
    rnSa = realFieldSa(0)
829
    outcome = sollya_lib_get_constant(get_rn_value(rnSa), ctExpSo)
830
    if outcome == 0:
831
        return None
832
    else:
833
        return rnSa
834
# End pobyso_get_constant_as_rn_with_rf_so_sa
835

    
836
def pobyso_get_free_variable_name():
837
    """ 
838
    Legacy function. See pobyso_get_free_variable_name_so_sa.
839
    """
840
    return(pobyso_get_free_variable_name_so_sa())
841

    
842
def pobyso_get_free_variable_name_so_sa():
843
    return sollya_lib_get_free_variable_name()
844
    
845
def pobyso_get_function_arity(expressionSo):
846
    """ 
847
    Legacy function. See pobyso_get_function_arity_so_sa.
848
    """
849
    return(pobyso_get_function_arity_so_sa(expressionSo))
850

    
851
def pobyso_get_function_arity_so_sa(expressionSo):
852
    arity = c_int(0)
853
    sollya_lib_get_function_arity(byref(arity),expressionSo)
854
    return int(arity.value)
855

    
856
def pobyso_get_head_function(expressionSo):
857
    """ 
858
    Legacy function. See pobyso_get_head_function_so_sa. 
859
    """
860
    return(pobyso_get_head_function_so_sa(expressionSo)) 
861

    
862
def pobyso_get_head_function_so_sa(expressionSo):
863
    functionType = c_int(0)
864
    sollya_lib_get_head_function(byref(functionType), expressionSo)
865
    return int(functionType.value)
866

    
867
def pobyso_get_interval_from_range_so_sa(soRange, realIntervalFieldSa = None ):
868
    """
869
    Return the Sage interval corresponding to the Sollya range argument.
870
    If no reaIntervalField is passed as an argument, the interval bounds are not
871
    rounded: they are elements of RealIntervalField of the "right" precision
872
    to hold all the digits.
873
    """
874
    prec = c_int(0)
875
    if realIntervalFieldSa is None:
876
        retval = sollya_lib_get_prec_of_range(byref(prec), soRange, None)
877
        if retval == 0:
878
            return None
879
        realIntervalFieldSa = RealIntervalField(prec.value)
880
    intervalSa = realIntervalFieldSa(0,0)
881
    retval = \
882
        sollya_lib_get_interval_from_range(get_interval_value(intervalSa),\
883
                                           soRange)
884
    if retval == 0:
885
        return None
886
    return intervalSa
887
# End pobyso_get_interval_from_range_so_sa
888

    
889
def pobyso_get_list_elements(soObj):
890
    """ Legacy function. See pobyso_get_list_elements_so_so. """
891
    return pobyso_get_list_elements_so_so(soObj)
892
 
893
def pobyso_get_list_elements_so_so(objectListSo):
894
    """
895
    Get the Sollya list elements as a Sage/Python array of Sollya objects.
896
    
897
    INPUT:
898
    - objectListSo: a Sollya list of Sollya objects.
899
    
900
    OUTPUT:
901
    - a Sage/Python tuple made of:
902
      - a Sage/Python list of Sollya objects,
903
      - a Sage/Python int holding the number of elements,
904
      - a Sage/Python int stating (!= 0) that the list is end-elliptic.
905
    NOTE::
906
        We recover the addresses of the Sollya object from the list of pointers
907
        returned by sollya_lib_get_list_elements. The list itself is freed.
908
    TODO::
909
        Figure out what to do with numElements since the number of elements
910
        can easily be recovered from the list itself. 
911
        Ditto for isEndElliptic.
912
    """
913
    listAddress = POINTER(c_longlong)()
914
    numElements = c_int(0)
915
    isEndElliptic = c_int(0)
916
    listAsSageList = []
917
    result = sollya_lib_get_list_elements(byref(listAddress),\
918
                                          byref(numElements),\
919
                                          byref(isEndElliptic),\
920
                                          objectListSo)
921
    if result == 0 :
922
        return None
923
    for i in xrange(0, numElements.value, 1):
924
       #listAsSageList.append(sollya_lib_copy_obj(listAddress[i]))
925
       listAsSageList.append(listAddress[i])
926
       # Clear each of the elements returned by Sollya.
927
       #sollya_lib_clear_obj(listAddress[i])
928
    # Free the list itself.   
929
    sollya_lib_free(listAddress)
930
    return (listAsSageList, numElements.value, isEndElliptic.value)
931

    
932
def pobyso_get_max_prec_of_exp(soExp):
933
    """ Legacy function. See pobyso_get_max_prec_of_exp_so_sa. """
934
    return pobyso_get_max_prec_of_exp_so_sa(soExp)
935

    
936
def pobyso_get_max_prec_of_exp_so_sa(expSo):
937
    """
938
    Get the maximum precision used for the numbers in a Sollya expression.
939
    
940
    Arguments:
941
    soExp -- a Sollya expression pointer
942
    Return value:
943
    A Python integer
944
    TODO: 
945
    - error management;
946
    - correctly deal with numerical type such as DOUBLEEXTENDED.
947
    """
948
    if expSo is None:
949
        print inspect.stack()[0][3], ": expSo is None."
950
        return 0
951
    maxPrecision = 0
952
    minConstPrec = 0
953
    currentConstPrec = 0
954
    #pobyso_autoprint(expSo)
955
    operator = pobyso_get_head_function_so_sa(expSo)
956
    if (operator != SOLLYA_BASE_FUNC_CONSTANT) and \
957
    (operator != SOLLYA_BASE_FUNC_FREE_VARIABLE):
958
        (arity, subexpressions) = pobyso_get_subfunctions_so_sa(expSo)
959
        for i in xrange(arity):
960
            maxPrecisionCandidate = \
961
                pobyso_get_max_prec_of_exp_so_sa(subexpressions[i])
962
            if maxPrecisionCandidate > maxPrecision:
963
                maxPrecision = maxPrecisionCandidate
964
        return maxPrecision
965
    elif operator == SOLLYA_BASE_FUNC_CONSTANT:
966
        #minConstPrec = pobyso_get_min_prec_of_constant_so_sa(expSo)
967
        #currentConstPrec = pobyso_get_min_prec_of_constant_so_sa(soExp)
968
        #print minConstPrec, " - ", currentConstPrec 
969
        return pobyso_get_min_prec_of_constant_so_sa(expSo)
970
    
971
    elif operator == SOLLYA_BASE_FUNC_FREE_VARIABLE:
972
        return 0
973
    else:
974
        print "pobyso_get_max_prec_of_exp_so_sa: unexepected operator."
975
        return 0
976

    
977
def pobyso_get_min_prec_of_constant_so_sa(constExpSo):
978
    """
979
    Get the minimum precision necessary to represent the value of a Sollya
980
    constant.
981
    MPFR_MIN_PREC and powers of 2 are taken into account.
982
    We assume that constExpSo is a pointer to a Sollay constant expression.
983
    """
984
    constExpAsRnSa = pobyso_get_constant_as_rn_so_sa(constExpSo)
985
    return(min_mpfr_size(get_rn_value(constExpAsRnSa)))
986

    
987
def pobyso_get_poly_so_sa(polySo, realFieldSa=None):
988
    """
989
    Convert a Sollya polynomial into a Sage polynomial.
990
    Legacy function. Use pobyso_float_poly_so_sa() instead.
991
    """    
992
    return pobyso_float_poly_so_sa(polySo,realFieldSa)
993
# End pobyso_get_poly_so_sa
994

    
995
def pobyso_get_prec():
996
    """ Legacy function. See pobyso_get_prec_so_sa(). """
997
    return pobyso_get_prec_so_sa()
998

    
999
def pobyso_get_prec_so():
1000
    """
1001
    Get the current default precision in Sollya.
1002
    The return value is a Sollya object.
1003
    Usefull when modifying the precision back and forth by avoiding
1004
    extra conversions.
1005
    """
1006
    return sollya_lib_get_prec(None)
1007
    
1008
def pobyso_get_prec_so_sa():
1009
    """
1010
    Get the current default precision in Sollya.
1011
    The return value is Sage/Python int.
1012
    """
1013
    precSo = sollya_lib_get_prec()
1014
    precSa = pobyso_constant_from_int_so_sa(precSo)
1015
    sollya_lib_clear_obj(precSo)
1016
    return precSa
1017
# End pobyso_get_prec_so_sa.
1018

    
1019
def pobyso_get_prec_so_so_sa():
1020
    """
1021
    Return the current precision both as a Sollya object and a
1022
    Sage integer as hybrid tuple.
1023
    To avoid multiple calls for precision manipulations. 
1024
    """
1025
    precSo = sollya_lib_get_prec()
1026
    precSa = pobyso_constant_from_int_so_sa(precSo)
1027
    return (precSo, int(precSa))
1028
    
1029
def pobyso_get_prec_of_constant(ctExpSo):
1030
    """ Legacy function. See pobyso_get_prec_of_constant_so_sa. """
1031
    return pobyso_get_prec_of_constant_so_sa(ctExpSo)
1032

    
1033
def pobyso_get_prec_of_constant_so_sa(ctExpSo):
1034
    """
1035
    Tries to find a precision to represent ctExpSo without rounding.
1036
    If not possible, returns None.
1037
    """
1038
    #print "Entering pobyso_get_prec_of_constant_so_sa..."
1039
    prec = c_int(0)
1040
    retc = sollya_lib_get_prec_of_constant(byref(prec), ctExpSo, None)
1041
    if retc == 0:
1042
        #print "pobyso_get_prec_of_constant_so_sa failed."
1043
        return None
1044
    #print "...exiting pobyso_get_prec_of_constant_so_sa."
1045
    return int(prec.value)
1046

    
1047
def pobyso_get_prec_of_range_so_sa(rangeSo):
1048
    """
1049
    Returns the number of bits elements of a range are coded with.
1050
    """
1051
    prec = c_int(0)
1052
    retc = sollya_lib_get_prec_of_range(byref(prec), rangeSo, None)
1053
    if retc == 0:
1054
        return(None)
1055
    return int(prec.value)
1056
# End pobyso_get_prec_of_range_so_sa()
1057

    
1058
def pobyso_get_sage_exp_from_sollya_exp(sollyaExpSo, realField = RR):
1059
    """ Legacy function. See pobyso_get_sage_exp_from_sollya_exp_so_sa. """
1060
    return pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, 
1061
                                                     realField = RR)
1062

    
1063
def pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, realFieldSa = RR):
1064
    """
1065
    Get a Sage expression from a Sollya expression. 
1066
    Currently only tested with polynomials with floating-point coefficients.
1067
    Notice that, in the returned polynomial, the exponents are RealNumbers.
1068
    """
1069
    #pobyso_autoprint(sollyaExp)
1070
    operatorSa = pobyso_get_head_function_so_sa(sollyaExpSo)
1071
    sollyaLibFreeVariableName = sollya_lib_get_free_variable_name()
1072
    ## Get rid of the "_"'s in "_x_", if any.
1073
    sollyaLibFreeVariableName = re.sub('_', '', sollyaLibFreeVariableName)
1074
    # Constants and the free variable are special cases.
1075
    # All other operator are dealt with in the same way.
1076
    if (operatorSa != SOLLYA_BASE_FUNC_CONSTANT) and \
1077
       (operatorSa != SOLLYA_BASE_FUNC_FREE_VARIABLE):
1078
        (aritySa, subexpressionsSa) = pobyso_get_subfunctions_so_sa(sollyaExpSo)
1079
        if aritySa == 1:
1080
            sageExpSa = eval(pobyso_function_type_as_string_so_sa(operatorSa) + \
1081
            "(" + pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], \
1082
            realFieldSa) + ")")
1083
        elif aritySa == 2:
1084
            # We do not get through the preprocessor.
1085
            # The "^" operator is then a special case.
1086
            if operatorSa == SOLLYA_BASE_FUNC_POW:
1087
                operatorAsStringSa = "**"
1088
            else:
1089
                operatorAsStringSa = \
1090
                    pobyso_function_type_as_string_so_sa(operatorSa)
1091
            sageExpSa = \
1092
              eval("pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], realFieldSa)"\
1093
              + " " + operatorAsStringSa + " " + \
1094
                   "pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[1], realFieldSa)")
1095
        # We do not know yet how to deal with arity >= 3 
1096
        # (is there any in Sollya anyway?).
1097
        else:
1098
            sageExpSa = eval('None')
1099
        return sageExpSa
1100
    elif operatorSa == SOLLYA_BASE_FUNC_CONSTANT:
1101
        #print "This is a constant"
1102
        return pobyso_get_constant_as_rn_with_rf_so_sa(sollyaExpSo, realFieldSa)
1103
    elif operatorSa == SOLLYA_BASE_FUNC_FREE_VARIABLE:
1104
        #print "This is the free variable"
1105
        return eval(sollyaLibFreeVariableName)
1106
    else:
1107
        print "Unexpected"
1108
        return eval('None')
1109
# End pobyso_get_sage_exp_from_sollya_exp_so_sa
1110

    
1111

    
1112
def pobyso_get_subfunctions(expressionSo):
1113
    """ Legacy function. See pobyso_get_subfunctions_so_sa. """
1114
    return pobyso_get_subfunctions_so_sa(expressionSo) 
1115
# End pobyso_get_subfunctions.
1116
 
1117
def pobyso_get_subfunctions_so_sa(expressionSo):
1118
    """
1119
    Get the subfunctions of an expression.
1120
    Return the number of subfunctions and the list of subfunctions addresses.
1121
    S.T.: Could not figure out another way than that ugly list of declarations
1122
    to recover the addresses of the subfunctions.
1123
    We limit ourselves to arity 8 functions. 
1124
    """
1125
    subf0 = c_int(0)
1126
    subf1 = c_int(0)
1127
    subf2 = c_int(0)
1128
    subf3 = c_int(0)
1129
    subf4 = c_int(0)
1130
    subf5 = c_int(0)
1131
    subf6 = c_int(0)
1132
    subf7 = c_int(0)
1133
    subf8 = c_int(0)
1134
    arity = c_int(0)
1135
    nullPtr = POINTER(c_int)()
1136
    sollya_lib_get_subfunctions(expressionSo, byref(arity), \
1137
      byref(subf0), byref(subf1), byref(subf2), byref(subf3), \
1138
      byref(subf4), byref(subf5),\
1139
      byref(subf6), byref(subf7), byref(subf8), nullPtr, None) 
1140
#    byref(cast(subfunctions[0], POINTER(c_int))), \
1141
#    byref(cast(subfunctions[0], POINTER(c_int))), \
1142
#    byref(cast(subfunctions[2], POINTER(c_int))), \
1143
#    byref(cast(subfunctions[3], POINTER(c_int))), \
1144
#    byref(cast(subfunctions[4], POINTER(c_int))), \
1145
#    byref(cast(subfunctions[5], POINTER(c_int))), \
1146
#    byref(cast(subfunctions[6], POINTER(c_int))), \
1147
#    byref(cast(subfunctions[7], POINTER(c_int))), \
1148
#    byref(cast(subfunctions[8], POINTER(c_int))), nullPtr)
1149
    subfunctions = [subf0, subf1, subf2, subf3, subf4, subf5, subf6, subf7, \
1150
                    subf8]
1151
    subs = []
1152
    if arity.value > pobyso_max_arity:
1153
        return(0,[])
1154
    for i in xrange(arity.value):
1155
        subs.append(int(subfunctions[i].value))
1156
        #print subs[i]
1157
    return (int(arity.value), subs)
1158
# End pobyso_get_subfunctions_so_sa
1159
    
1160
def pobyso_guess_degree_sa_sa(functionSa, intervalSa, approxErrorSa, 
1161
                              weightSa=None, degreeBoundSa=None):
1162
    """
1163
    Sa_sa variant of the solly_guessdegree function.
1164
    Return 0 if something goes wrong.
1165
    """
1166
    functionAsStringSa = functionSa._assume_str().replace('_SAGE_VAR_', '')
1167
    functionSo = pobyso_parse_string_sa_so(functionAsStringSa)
1168
    if pobyso_is_error_so_sa(functionSo):
1169
        sollya_lib_clear_obj(functionSo)
1170
        return 0
1171
    rangeSo = pobyso_interval_to_range_sa_so(intervalSa)
1172
    # The approximation error is expected to be a floating point number.
1173
    if pobyso_is_floating_point_number_sa_sa(approxErrorSa):
1174
        approxErrorSo = pobyso_constant_sa_so(approxErrorSa)
1175
    else:
1176
        approxErrorSo = pobyso_constant_sa_so(RR(approxErrorSa))
1177
    if not weightSa is None:
1178
        weightAsStringSa = weightSa._assume_str().replace('_SAGE_VAR_', '')
1179
        weightSo = pobyso_parse_string_sa_so(weightAsStringSa)
1180
        if pobyso_is_error_so_sa(weightSo):
1181
            sollya_lib_clear_obj(functionSo)
1182
            sollya_lib_clear_obj(rangeSo)
1183
            sollya_lib_clear_obj(approxErrorSo)   
1184
            sollya_lib_clear_obj(weightSo)
1185
            return 0   
1186
    else:
1187
        weightSo = None
1188
    if not degreeBoundSa is None:
1189
        degreeBoundSo = pobyso_constant_from_int_sa_so(degreeBoundSa)
1190
    else:
1191
        degreeBoundSo = None
1192
    guessedDegreeSa = pobyso_guess_degree_so_sa(functionSo,
1193
                                                rangeSo,
1194
                                                approxErrorSo,
1195
                                                weightSo,
1196
                                                degreeBoundSo)
1197
    sollya_lib_clear_obj(functionSo)
1198
    sollya_lib_clear_obj(rangeSo)
1199
    sollya_lib_clear_obj(approxErrorSo)
1200
    if not weightSo is None:
1201
        sollya_lib_clear_obj(weightSo)
1202
    if not degreeBoundSo is None:
1203
        sollya_lib_clear_obj(degreeBoundSo)
1204
    return guessedDegreeSa
1205
# End poyso_guess_degree_sa_sa
1206

    
1207
def pobyso_guess_degree_so_sa(functionSo, rangeSo, errorSo, weightSo=None, \
1208
                              degreeBoundSo=None):
1209
    """
1210
    Thin wrapper around the guessdegree function.
1211
    Nevertheless, some precision control stuff has been appended.
1212
    """
1213
    # Deal with Sollya internal precision issues: if it is too small, 
1214
    # compared with the error, increases it to about twice -log2(error).
1215
    errorSa = pobyso_get_constant_as_rn_with_rf_so_sa(errorSo)
1216
    log2ErrorSa = errorSa.log2()
1217
    if log2ErrorSa < 0:
1218
        neededPrecisionSa = int(2 * int(-log2ErrorSa) / 64) * 64
1219
    else:
1220
        neededPrecisionSa = int(2 * int(log2ErrorSa) / 64) * 64
1221
    #print "Needed precision:", neededPrecisionSa
1222
    sollyaPrecisionHasChanged = False
1223
    (initialPrecSo, initialPrecSa) = pobyso_get_prec_so_so_sa()
1224
    if neededPrecisionSa > initialPrecSa:
1225
        if neededPrecisionSa <= 2:
1226
            print inspect.stack()[0][3], ": precision change <= 2 requested."
1227
        pobyso_set_prec_sa_so(neededPrecisionSa)
1228
        sollyaPrecisionHasChanged = True
1229
    #print "Guessing degree..."
1230
    # weightSo and degreeBoundsSo are optional arguments.
1231
    # As declared, sollya_lib_guessdegree must take 5 arguments.
1232
    if weightSo is None:
1233
        degreeRangeSo = sollya_lib_guessdegree(functionSo, rangeSo, errorSo, 
1234
                                               0, 0, None)
1235
    elif degreeBoundSo is None:
1236
        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, \
1237
                                                errorSo, weightSo, 0, None)
1238
    else:
1239
        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, errorSo, \
1240
                                                weightSo, degreeBoundSo, None)
1241
    #print "...degree guess done."
1242
    # Restore internal precision, if applicable.
1243
    if sollyaPrecisionHasChanged:
1244
        pobyso_set_prec_so_so(initialPrecSo)
1245
        sollya_lib_clear_obj(initialPrecSo)
1246
    degreeIntervalSa = pobyso_range_to_interval_so_sa(degreeRangeSo)
1247
    sollya_lib_clear_obj(degreeRangeSo)
1248
    # When ok, both bounds match.
1249
    # When the degree bound is too low, the upper bound is the degree
1250
    # for which the error can be honored.
1251
    # When it really goes wrong, the upper bound is infinity. 
1252
    if degreeIntervalSa.lower() == degreeIntervalSa.upper():
1253
        return int(degreeIntervalSa.lower())
1254
    else:
1255
        if degreeIntervalSa.upper().is_infinity():
1256
            return None
1257
        else:
1258
            return int(degreeIntervalSa.upper())
1259
    # End pobyso_guess_degree_so_sa
1260

    
1261
def pobyso_inf_so_so(intervalSo):
1262
    """
1263
    Very thin wrapper around sollya_lib_inf().
1264
    """
1265
    return sollya_lib_inf(intervalSo)
1266
# End pobyso_inf_so_so.
1267
#    
1268
def pobyso_infnorm_sa_sa(funcSa, intervalSa):
1269
    """
1270
    An infnorm call with Sage arguments.
1271
    We only take into account the 2 first arguments (the function and
1272
    the interval (a range). Managing the other arguments (the file for
1273
    the proof and the exclusion intervals list) will be performed later
1274
    Changes will be needed in sollya_lib.py file too.
1275
    """
1276
    # Check that funcSa is a function.
1277
    if not \
1278
     sage.symbolic.callable.is_CallableSymbolicExpressionRing(parent(funcSa)):
1279
        return None
1280
    # Check that intervalSa is an interval.
1281
    try:
1282
        intervalSa.upper()
1283
    except AttributeError:
1284
        return None
1285
    # Convert the Sage function into a Sollya function.
1286
    funcAsStringSa = funcSa._assume_str().replace('_SAGE_VAR_', '')
1287
    funcSo = pobyso_parse_string_sa_so(funcAsStringSa)
1288
    if not pobyso_obj_is_function_so_sa(funcSo):
1289
        sollya_lib_clear_obj(funcSo)
1290
        return None
1291
    # Convert the Sage interval into a Sollya range.
1292
    rangeSo = pobyso_interval_to_range_sa_so(intervalSa)
1293
    retValSo = sollya_lib_infnorm(funcSo, rangeSo, None)
1294
    sollya_lib_clear_obj(funcSo)
1295
    sollya_lib_clear_obj(rangeSo)
1296
    if pobyso_is_error_so_sa(retValSo):
1297
        sollya_lib_clear_obj(retValSo)
1298
        return None
1299
    retValSa = pobyso_range_to_interval_so_sa(retValSo)
1300
    sollya_lib_clear_obj(retValSo)
1301
    return retValSa
1302
# End pobyso_infnorm_so_so.
1303
#   
1304
def pobyso_infnorm_so_so(funcSo, rangeSo):
1305
    """
1306
    Very thin wrapper around sollya_lib_infnorm().
1307
    We only take into account the 2 first arguments (the function and
1308
    the interval (a range). Managing the other arguments (the file for
1309
    the proof and the exclusion intervals list) will be performed later
1310
    Changes will be needed in sollya_lib.py file too.
1311
    
1312
    As per Sollya manual, this function should not be used anymore and
1313
    supnorm should be called instead. Nevertheless, supnorm breaks 
1314
    sometimes whereas infnorm still returns a satisfactory answer.
1315
    """
1316
    return sollya_lib_infnorm(funcSo, rangeSo, None)
1317
# End pobyso_infnorm_so_so.
1318

    
1319
def pobyso_interval_to_range_sa_so(intervalSa, precisionSa=None):
1320
    if precisionSa is None:
1321
        precisionSa = intervalSa.parent().precision()
1322
    intervalSo = pobyso_bounds_to_range_sa_so(intervalSa.lower(),\
1323
                                              intervalSa.upper(),\
1324
                                              precisionSa)
1325
    return intervalSo
1326
# End pobyso_interval_to_range_sa_so
1327

    
1328
def pobyso_is_error_so_sa(objSo):
1329
    """
1330
    Thin wrapper around the sollya_lib_obj_is_error() function.
1331
    """
1332
    if sollya_lib_obj_is_error(objSo) != 0:
1333
        return True
1334
    else:
1335
        return False
1336
# End pobyso_is_error-so_sa
1337

    
1338
def pobyso_is_floating_point_number_sa_sa(numberSa):
1339
    """
1340
    Check whether a Sage number is floating point.
1341
    Exception stuff added because numbers other than
1342
    floating-point ones do not have the is_real() attribute.
1343
    """
1344
    try:
1345
        return numberSa.is_real()
1346
    except AttributeError:
1347
        return False
1348
# End pobyso_is_floating_piont_number_sa_sa
1349

    
1350
def pobyso_is_range_so_sa(rangeCandidateSo):
1351
    """
1352
    Thin wrapper over sollya_lib_is_range.
1353
    """
1354
    return sollya_lib_obj_is_range(rangeCandidateSo) != 0
1355
    
1356
# End pobyso_is_range_so_sa
1357

    
1358

    
1359
def pobyso_lib_init():
1360
    sollya_lib_init(None)
1361

    
1362
def pobyso_lib_close():
1363
    sollya_lib_close(None)
1364

    
1365
def pobyso_name_free_variable(freeVariableNameSa):
1366
    """ Legacy function. See pobyso_name_free_variable_sa_so. """
1367
    pobyso_name_free_variable_sa_so(freeVariableNameSa)
1368

    
1369
def pobyso_name_free_variable_sa_so(freeVariableNameSa):
1370
    """
1371
    Set the free variable name in Sollya from a Sage string.
1372
    """
1373
    sollya_lib_name_free_variable(freeVariableNameSa)
1374

    
1375
def pobyso_obj_is_function_so_sa(objSo):
1376
    """
1377
    Check if an object is a function.
1378
    """
1379
    if sollya_lib_obj_is_function(objSo) != 0:
1380
        return True
1381
    else:
1382
        return False
1383
# End pobyso_obj_is_function_so_sa  
1384
    
1385
def pobyso_obj_is_range_so_sa(objSo):
1386
    """
1387
    Check if an object is a function.
1388
    """
1389
    if sollya_lib_obj_is_range(objSo) != 0:
1390
        return True
1391
    else:
1392
        return False
1393
# End pobyso_obj_is_range_so_sa  
1394
    
1395
def pobyso_obj_is_string_so_sa(objSo):
1396
    """
1397
    Check if an object is a function.
1398
    """
1399
    if sollya_lib_obj_is_string(objSo) != 0:
1400
        return True
1401
    else:
1402
        return False
1403
# End pobyso_obj_is_string_so_sa  
1404
    
1405
def pobyso_parse_string(string):
1406
    """ Legacy function. See pobyso_parse_string_sa_so. """
1407
    return pobyso_parse_string_sa_so(string)
1408
 
1409
def pobyso_parse_string_sa_so(string):
1410
    """
1411
    Get the Sollya expression computed from a Sage string or
1412
    a Sollya error object if parsing failed.
1413
    """
1414
    return sollya_lib_parse_string(string)
1415

    
1416
def pobyso_precision_so_sa(ctExpSo):
1417
    """
1418
    Computes the necessary precision to represent a number.
1419
    If x is not zero, it can be uniquely written as x = m · 2e 
1420
    where m is an odd integer and e is an integer. 
1421
    precision(x) returns the number of bits necessary to write m 
1422
    in binary (i.e. ceil(log2(m))).
1423
    """
1424
    #TODO: take care of the special case: 0, @NaN@, @Inf@
1425
    precisionSo = sollya_lib_precision(ctExpSo)
1426
    precisionSa = pobyso_constant_from_int_so_sa(precisionSo)
1427
    sollya_lib_clear_obj(precisionSo)
1428
    return precisionSa
1429
# End pobyso_precision_so_sa
1430

    
1431
def pobyso_polynomial_coefficients_progressive_round_so_so(polySo,
1432
                                                           funcSo,
1433
                                                           icSo,
1434
                                                           intervalSo,
1435
                                                           itpSo,
1436
                                                           ftpSo,
1437
                                                           maxPrecSo,
1438
                                                           maxErrSo,
1439
                                                           debug=False):
1440
    if debug:
1441
        print "Input arguments:"
1442
        pobyso_autoprint(polySo)
1443
        pobyso_autoprint(funcSo)
1444
        pobyso_autoprint(icSo)
1445
        pobyso_autoprint(intervalSo)
1446
        pobyso_autoprint(itpSo)
1447
        pobyso_autoprint(ftpSo)
1448
        pobyso_autoprint(maxPrecSo)
1449
        pobyso_autoprint(maxErrSo)
1450
        print "________________"
1451
    
1452
    ## Higher order function see: 
1453
    #  http://effbot.org/pyfaq/how-do-you-make-a-higher-order-function-in-python.htm
1454
    def precision_decay_ratio_function(degreeSa):
1455
        def outer(x):
1456
            def inner(x):
1457
                we = 3/8
1458
                wq = 2/8
1459
                a  = 2.2
1460
                b  = 2 
1461
                return we*(exp(x/a)-1) +  wq*((b*x)**2) + (1-we-wq)*x
1462
            return  inner(x)/inner(degreeSa)
1463
        return outer
1464
    
1465
    #
1466
    degreeSa        = pobyso_polynomial_degree_so_sa(polySo)
1467
    ratio           = precision_decay_ratio_function(degreeSa)
1468
    itpSa           = pobyso_constant_from_int_so_sa(itpSo)
1469
    ftpSa           = pobyso_constant_from_int_so_sa(ftpSo)
1470
    maxPrecSa       = pobyso_constant_from_int_so_sa(maxPrecSo)
1471
    maxErrSa        = pobyso_get_constant_as_rn_so_sa(maxErrSo)
1472
    if debug:
1473
        print "degreeSa:", degreeSa
1474
        print "ratio:", ratio
1475
        print "itpsSa:", itpSa
1476
        print "ftpSa:", ftpSa
1477
        print "maxPrecSa:", maxPrecSa
1478
        print "maxErrSa:", maxErrSa
1479
    lastResPolySo   = None
1480
    lastInfNormSo   = None
1481
    #print "About to enter the while loop..."
1482
    while True: 
1483
        resPolySo   = pobyso_constant_0_sa_so()
1484
        pDeltaSa    = ftpSa - itpSa
1485
        for indexSa in reversed(xrange(0,degreeSa+1)):
1486
            #print "Index:", indexSa
1487
            indexSo = pobyso_constant_from_int_sa_so(indexSa)
1488
            #pobyso_autoprint(indexSo)
1489
            #print ratio(indexSa)
1490
            ctpSa = floor(ftpSa - (pDeltaSa * ratio(indexSa)))
1491
            ctpSo = pobyso_constant_from_int_sa_so(ctpSa)
1492
            if debug:
1493
                print "Index:", indexSa, " - Target precision:", 
1494
                pobyso_autoprint(ctpSo)
1495
            cmonSo  = \
1496
                sollya_lib_build_function_mul(sollya_lib_coeff(polySo, indexSo),
1497
                                      sollya_lib_build_function_pow( \
1498
                                          sollya_lib_build_function_free_variable(), \
1499
                                          indexSo))
1500
            #pobyso_autoprint(cmonSo)
1501
            cmonrSo = pobyso_round_coefficients_single_so_so(cmonSo, ctpSo)
1502
            sollya_lib_clear_obj(cmonSo)
1503
            #pobyso_autoprint(cmonrSo)
1504
            resPolySo = sollya_lib_build_function_add(resPolySo,
1505
                                                      cmonrSo)
1506
            #pobyso_autoprint(resPolySo)
1507
        # End for index
1508
        freeVarSo     = sollya_lib_build_function_free_variable()
1509
        changeVarSo   = sollya_lib_sub(freeVarSo, icSo)
1510
        resPolyCvSo   = sollya_lib_evaluate(resPolySo, changeVarSo)
1511
        errFuncSo = sollya_lib_build_function_sub(sollya_lib_copy_obj(funcSo), 
1512
                                                  resPolyCvSo)
1513
        infNormSo = sollya_lib_dirtyinfnorm(errFuncSo, intervalSo)
1514
        cerrSa    = pobyso_get_constant_as_rn_so_sa(infNormSo)
1515
        if debug:
1516
            print "Infnorm (Sollya):", 
1517
            pobyso_autoprint(infNormSo)
1518
        sollya_lib_clear_obj(errFuncSo)
1519
        #print "Infnorm  (Sage):", cerrSa
1520
        if (cerrSa > maxErrSa):
1521
            if debug:
1522
                print "Error is too large."
1523
            if lastResPolySo is None:
1524
                if debug:
1525
                    print "Enlarging prec."
1526
                ntpSa = floor(ftpSa + ftpSa/50)
1527
                ## Can't enlarge (numerical)
1528
                if ntpSa == ftpSa:
1529
                    sollya_lib_clear_obj(resPolySo)
1530
                    return None
1531
                ## Can't enlarge (not enough precision left)
1532
                if ntpSa > maxPrecSa:
1533
                    sollya_lib_clear_obj(resPolySo)
1534
                    return None
1535
                ftpSa = ntpSa
1536
                continue
1537
            ## One enlargement took place.
1538
            else:
1539
                if debug:
1540
                    print "Exit with the last before last polynomial."
1541
                    print "Precision of highest degree monomial:", itpSa
1542
                    print "Precision of constant term          :", ftpSa
1543
                sollya_lib_clear_obj(resPolySo)
1544
                sollya_lib_clear_obj(infNormSo)
1545
                return (lastResPolySo, lastInfNormSo)
1546
        # cerrSa <= maxErrSa: scrap more bits, possibly.
1547
        else:
1548
            if debug:
1549
                print "Error is too small"
1550
            if cerrSa <= (maxErrSa/2):
1551
                if debug: 
1552
                    print "Shrinking prec."
1553
                ntpSa = floor(ftpSa - ftpSa/50)
1554
                ## Can't shrink (numerical)
1555
                if ntpSa == ftpSa:
1556
                    if not lastResPolySo is None:
1557
                        sollya_lib_clear_obj(lastResPolySo)
1558
                    if not lastInfNormSo is None:
1559
                        sollya_lib_clear_obj(lastInfNormSo)
1560
                    if debug:
1561
                        print "Exit because can't shrink anymore (numerically)."
1562
                        print "Precision of highest degree monomial:", itpSa
1563
                        print "Precision of constant term          :", ftpSa
1564
                    return (resPolySo, infNormSo)
1565
                ## Can't shrink (not enough precision left)
1566
                if ntpSa <= itpSa:
1567
                    if not lastResPolySo is None:
1568
                        sollya_lib_clear_obj(lastResPolySo)
1569
                    if not lastInfNormSo is None:
1570
                        sollya_lib_clear_obj(lastInfNormSo)
1571
                        print "Exit because can't shrink anymore (no bits left)."
1572
                        print "Precision of highest degree monomial:", itpSa
1573
                        print "Precision of constant term          :", ftpSa
1574
                    return (resPolySo, infNormSo)
1575
                ftpSa = ntpSa
1576
                if not lastResPolySo is None:
1577
                    sollya_lib_clear_obj(lastResPolySo)
1578
                if not lastInfNormSo is None:
1579
                    sollya_lib_clear_obj(lastInfNormSo)
1580
                lastResPolySo = resPolySo
1581
                lastInfNormSo = infNormSo
1582
                continue
1583
            else: # Error is not that small, just return
1584
                if not lastResPolySo is None:
1585
                    sollya_lib_clear_obj(lastResPolySo)
1586
                if not lastInfNormSo is None:
1587
                    sollya_lib_clear_obj(lastInfNormSo)
1588
                if debug:
1589
                    print "Exit normally."
1590
                    print "Precision of highest degree monomial:", itpSa
1591
                    print "Precision of constant term          :", ftpSa
1592
                return (resPolySo, infNormSo)                
1593
    # End wile True
1594
    return None
1595
# End pobyso_polynomial_coefficients_progressive_truncate_so_so.
1596

    
1597
def pobyso_polynomial_degree_so_sa(polySo):
1598
    """
1599
    Return the degree of a Sollya polynomial as a Sage int.
1600
    """
1601
    degreeSo = sollya_lib_degree(polySo)
1602
    return pobyso_constant_from_int_so_sa(degreeSo)
1603
# End pobyso_polynomial_degree_so_sa
1604

    
1605
def pobyso_polynomial_degree_so_so(polySo):
1606
    """
1607
    Thin wrapper around lib_sollya_degree().
1608
    """
1609
    return sollya_lib_degree(polySo)
1610
# End pobyso_polynomial_degree_so_so
1611
                                                                  
1612
def pobyso_range(rnLowerBound, rnUpperBound):
1613
    """ Legacy function. See pobyso_range_sa_so. """
1614
    return pobyso_range_sa_so(rnLowerBound, rnUpperBound) 
1615

    
1616
def pobyso_range_from_bounds_sa_so(rnLowerBound, rnUpperBound, precSa = None):
1617
    """
1618
    Create a Sollya range from 2 Sage real numbers as bounds
1619
    """
1620
    # TODO check precision stuff.
1621
    sollyaPrecChanged = False
1622
    (initialSollyaPrecSo, initialSollyaPrecSa) = \
1623
        pobyso_get_prec_so_so_sa()
1624
    if precSa is None:
1625
        precSa = max(rnLowerBound.parent().prec(), rnUpperBound.parent().prec())
1626
    if precSa > initialSollyaPrecSa:
1627
        if precSa <= 2:
1628
            print inspect.stack()[0][3], ": precision change <= 2 requested."
1629
        pobyso_set_prec_sa_so(precSa)
1630
        sollyaPrecChanged = True
1631
    rangeSo = sollya_lib_range_from_bounds(get_rn_value(rnLowerBound),
1632
                                           get_rn_value(rnUpperBound))
1633
    if sollyaPrecChanged:
1634
        pobyso_set_prec_so_so(initialSollyaPrecSo)
1635
        sollya_lib_clear_obj(initialSollyaPrecSo)
1636
    return rangeSo
1637
# End pobyso_range_from_bounds_sa_so
1638

    
1639
def pobyso_range_list_so_sa(listSo):
1640
    """
1641
    Return a Sollya list of ranges as a Sage list of 
1642
    floating-point intervals.
1643
    """
1644
    listSa   = []
1645
    ## The function returns none if the list is empty or an error has happened.
1646
    retVal = pobyso_get_list_elements_so_so(listSo)
1647
    if retVal is None:
1648
        return listSa
1649
    ## Just in case the interface is changed and an empty list is returned 
1650
    #  instead of None.
1651
    elif len(retVal) == 0:
1652
        return listSa
1653
    else:
1654
        ## Remember pobyso_get_list_elements_so_so returns more information
1655
        #  than just the elements of the list (# elements, is_elliptic)
1656
        listSaSo, numElements, isEndElliptic = retVal
1657
    ## Return an empty list.
1658
    if numElements == 0:
1659
        return listSa
1660
    ## Search first for the maximum precision of the elements
1661
    maxPrecSa = 0
1662
    for rangeSo in listSaSo:
1663
        #pobyso_autoprint(floatSo)
1664
        curPrecSa =  pobyso_get_prec_of_range_so_sa(rangeSo)
1665
        if curPrecSa > maxPrecSa:
1666
            maxPrecSa = curPrecSa
1667
    ##
1668
    intervalField = RealIntervalField(maxPrecSa)
1669
    ##
1670
    for rangeSo in listSaSo:
1671
        listSa.append(pobyso_range_to_interval_so_sa(rangeSo, intervalField))
1672
    return listSa
1673
# End pobyso_range_list_so_sa
1674

    
1675
def pobyso_range_max_abs_so_so(rangeSo):
1676
    """
1677
    Return, as Sollya constant, the maximum absolute value of a Sollay range.
1678
    """
1679
    lowerBoundSo = sollya_lib_inf(rangeSo)
1680
    upperBoundSo = sollya_lib_sup(rangeSo)
1681
    #
1682
    lowerBoundSo = sollya_lib_build_function_abs(lowerBoundSo)
1683
    upperBoundSo = sollya_lib_build_function_abs(upperBoundSo)
1684
    #pobyso_autoprint(lowerBoundSo)
1685
    #pobyso_autoprint(upperBoundSo)
1686
    #
1687
    maxAbsSo = sollya_lib_max(lowerBoundSo, upperBoundSo, None)
1688
    #sollya_lib_clear_obj(lowerBoundSo)
1689
    #sollya_lib_clear_obj(upperBoundSo)
1690
    return maxAbsSo
1691
# End pobyso_range_max_abs_so_so
1692
 
1693
def pobyso_range_to_interval_so_sa(rangeSo, realIntervalFieldSa = None):
1694
    """
1695
    Get a Sage interval from a Sollya range.
1696
    If no realIntervalField is given as a parameter, the Sage interval
1697
    precision is that of the Sollya range.
1698
    Otherwise, the precision is that of the realIntervalField. In this case
1699
    rounding may happen.
1700
    """
1701
    if realIntervalFieldSa is None:
1702
        precSa = pobyso_get_prec_of_range_so_sa(rangeSo)
1703
        realIntervalFieldSa = RealIntervalField(precSa)
1704
    intervalSa = \
1705
        pobyso_get_interval_from_range_so_sa(rangeSo, realIntervalFieldSa) 
1706
    return intervalSa
1707
# End pobyso_range_to_interval_so_sa
1708
#
1709
def pobyso_relative_so_so():
1710
    """
1711
    Very thin wrapper around the sollya_lib_relative function.
1712
    """
1713
    return sollya_lib_relative()
1714
# End pobyso_relative_so_so
1715
#
1716
def pobyso_rat_poly_sa_so(polySa, precSa = None):
1717
    """
1718
    Create a Sollya polynomial from a Sage rational polynomial.
1719
    We first convert the rational polynomial into a floating-point 
1720
    polynomial.
1721
    """
1722
    ## TODO: filter arguments.
1723
    ## Precision. If no precision is given, use the current precision
1724
    #  of Sollya.
1725
    if precSa is None:
1726
        precSa =  pobyso_get_prec_so_sa()
1727
    #print "Precision:",  precSa
1728
    RRR = RealField(precSa)
1729
    ## Create a Sage polynomial in the "right" precision.
1730
    P_RRR = RRR[polySa.variables()[0]]
1731
    polyFloatSa = P_RRR(polySa)
1732
    ## Make sure no precision is provided: pobyso_float_poly_sa_so will
1733
    #  recover it all by itself and will not make any extra conversion.
1734
    return pobyso_float_poly_sa_so(polyFloatSa)
1735
    
1736
# End pobyso_rat_poly_sa_so    
1737
    
1738
def pobyso_remez_canonical_sa_sa(func, \
1739
                                 degree, \
1740
                                 lowerBound, \
1741
                                 upperBound, \
1742
                                 weight = None, \
1743
                                 quality = None):
1744
    """
1745
    All arguments are Sage/Python.
1746
    The functions (func and weight) must be passed as expressions or strings.
1747
    Otherwise the function fails. 
1748
    The return value is a Sage polynomial.
1749
    """
1750
    var('zorglub')    # Dummy variable name for type check only. Type of 
1751
    # zorglub is "symbolic expression".
1752
    polySo = pobyso_remez_canonical_sa_so(func, \
1753
                                 degree, \
1754
                                 lowerBound, \
1755
                                 upperBound, \
1756
                                 weight, \
1757
                                 quality)
1758
    # String test
1759
    if parent(func) == parent("string"):
1760
        functionSa = eval(func)
1761
    # Expression test.
1762
    elif type(func) == type(zorglub):
1763
        functionSa = func
1764
    else:
1765
        return None
1766
    #
1767
    maxPrecision = 0
1768
    if polySo is None:
1769
        return(None)
1770
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
1771
    RRRRSa = RealField(maxPrecision)
1772
    polynomialRingSa = RRRRSa[functionSa.variables()[0]]
1773
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, RRRRSa)
1774
    polySa = polynomial(expSa, polynomialRingSa)
1775
    sollya_lib_clear_obj(polySo)
1776
    return(polySa)
1777
# End pobyso_remez_canonical_sa_sa
1778
    
1779
def pobyso_remez_canonical(func, \
1780
                           degree, \
1781
                           lowerBound, \
1782
                           upperBound, \
1783
                           weight = "1", \
1784
                           quality = None):
1785
    """ Legacy function. See pobyso_remez_canonical_sa_so. """
1786
    return(pobyso_remez_canonical_sa_so(func, \
1787
                                        degree, \
1788
                                        lowerBound, \
1789
                                        upperBound, \
1790
                                        weight, \
1791
                                        quality))
1792
# End pobyso_remez_canonical.
1793

    
1794
def pobyso_remez_canonical_sa_so(func, \
1795
                                 degree, \
1796
                                 lowerBound, \
1797
                                 upperBound, \
1798
                                 weight = None, \
1799
                                 quality = None):
1800
    """
1801
    All arguments are Sage/Python.
1802
    The functions (func and weight) must be passed as expressions or strings.
1803
    Otherwise the function fails. 
1804
    The return value is a pointer to a Sollya function.
1805
    lowerBound and upperBound mus be reals.
1806
    """
1807
    var('zorglub')    # Dummy variable name for type check only. Type of
1808
    # zorglub is "symbolic expression".
1809
    currentVariableNameSa = None
1810
    # The func argument can be of different types (string, 
1811
    # symbolic expression...)
1812
    if parent(func) == parent("string"):
1813
        localFuncSa = sage_eval(func,globals())
1814
        if len(localFuncSa.variables()) > 0:
1815
            currentVariableNameSa = localFuncSa.variables()[0]
1816
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1817
            functionSo = \
1818
              sollya_lib_parse_string(localFuncSa._assume_str().replace('_SAGE_VAR_', ''))
1819
    # Expression test.
1820
    elif type(func) == type(zorglub):
1821
        # Until we are able to translate Sage expressions into Sollya 
1822
        # expressions : parse the string version.
1823
        if len(func.variables()) > 0:
1824
            currentVariableNameSa = func.variables()[0]
1825
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1826
            functionSo = \
1827
              sollya_lib_parse_string(func._assume_str().replace('_SAGE_VAR_', ''))
1828
    else:
1829
        return(None)
1830
    if weight is None: # No weight given -> 1.
1831
        weightSo = pobyso_constant_1_sa_so()
1832
    elif parent(weight) == parent("string"): # Weight given as string: parse it.
1833
        weightSo = sollya_lib_parse_string(func)
1834
    elif type(weight) == type(zorglub): # Weight given as symbolice expression.
1835
        functionSo = \
1836
          sollya_lib_parse_string_sa_so(weight._assume_str().replace('_SAGE_VAR_', ''))
1837
    else:
1838
        return(None)
1839
    degreeSo = pobyso_constant_from_int(degree)
1840
    rangeSo = pobyso_bounds_to_range_sa_so(lowerBound, upperBound)
1841
    if not quality is None:
1842
        qualitySo= pobyso_constant_sa_so(quality)
1843
    else:
1844
        qualitySo = None
1845
        
1846
    remezPolySo = sollya_lib_remez(functionSo, \
1847
                                   degreeSo, \
1848
                                   rangeSo, \
1849
                                   weightSo, \
1850
                                   qualitySo, \
1851
                                   None)
1852
    sollya_lib_clear_obj(functionSo)
1853
    sollya_lib_clear_obj(degreeSo)
1854
    sollya_lib_clear_obj(rangeSo)
1855
    sollya_lib_clear_obj(weightSo)
1856
    if not qualitySo is None:
1857
        sollya_lib_clear_obj(qualitySo)
1858
    return(remezPolySo)
1859
# End pobyso_remez_canonical_sa_so
1860

    
1861
def pobyso_remez_canonical_so_so(funcSo, \
1862
                                 degreeSo, \
1863
                                 rangeSo, \
1864
                                 weightSo = pobyso_constant_1_sa_so(),\
1865
                                 qualitySo = None):
1866
    """
1867
    All arguments are pointers to Sollya objects.
1868
    The return value is a pointer to a Sollya function.
1869
    """
1870
    if not sollya_lib_obj_is_function(funcSo):
1871
        return(None)
1872
    return(sollya_lib_remez(funcSo, degreeSo, rangeSo, weightSo, qualitySo, None))
1873
# End pobyso_remez_canonical_so_so.
1874

    
1875
def pobyso_remez_exponents_list_sa_so(func,     \
1876
                                 exponentsList, \
1877
                                 lowerBound,    \
1878
                                 upperBound,    \
1879
                                 weight = None, \
1880
                                 quality = None):
1881
    """
1882
    All arguments are Sage/Python.
1883
    The functions (func and weight) must be passed as expressions or strings.
1884
    Otherwise the function fails. 
1885
    The return value is a pointer to a Sollya function.
1886
    lowerBound and upperBound mus be reals.
1887
    """
1888
    var('zorglub')    # Dummy variable name for type check only. Type of
1889
    # zorglub is "symbolic expression".
1890
    currentVariableNameSa = None
1891
    # The func argument can be of different types (string, 
1892
    # symbolic expression...)
1893
    if parent(func) == parent("string"):
1894
        localFuncSa = sage_eval(func,globals())
1895
        if len(localFuncSa.variables()) > 0:
1896
            currentVariableNameSa = localFuncSa.variables()[0]
1897
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1898
            functionSo = \
1899
              sollya_lib_parse_string(localFuncSa._assume_str().replace('_SAGE_VAR_', ''))
1900
    # Expression test.
1901
    elif type(func) == type(zorglub):
1902
        # Until we are able to translate Sage expressions into Sollya 
1903
        # expressions : parse the string version.
1904
        if len(func.variables()) > 0:
1905
            currentVariableNameSa = func.variables()[0]
1906
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1907
            functionSo = \
1908
              sollya_lib_parse_string(func._assume_str().replace('_SAGE_VAR_', ''))
1909
    else:
1910
        return(None)
1911
    ## Deal with the weight, much in the same way as with the function.
1912
    if weight is None: # No weight given -> 1.
1913
        weightSo = pobyso_constant_1_sa_so()
1914
    elif parent(weight) == parent("string"): # Weight given as string: parse it.
1915
        weightSo = sollya_lib_parse_string(func)
1916
    elif type(weight) == type(zorglub): # Weight given as symbolice expression.
1917
        functionSo = \
1918
          sollya_lib_parse_string(weight._assume_str().replace('_SAGE_VAR_', '',100))
1919
    else:
1920
        return(None)
1921
    rangeSo = pobyso_bounds_to_range_sa_so(lowerBound, upperBound)
1922
    if not quality is None:
1923
        qualitySo= pobyso_constant_sa_so(quality)
1924
    else:
1925
        qualitySo = None
1926
    #
1927
    ## Tranform the Sage list of exponents into a Sollya list.
1928
    exponentsListSo = pobyso_build_list_of_ints_sa_so(*exponentsList)
1929
    remezPolySo = sollya_lib_remez(functionSo,      \
1930
                                   exponentsListSo, \
1931
                                   rangeSo,         \
1932
                                   weightSo,        \
1933
                                   qualitySo,       \
1934
                                   None)
1935
    sollya_lib_clear_obj(functionSo)
1936
    sollya_lib_clear_obj(exponentsListSo)
1937
    sollya_lib_clear_obj(rangeSo)
1938
    sollya_lib_clear_obj(weightSo)
1939
    if not qualitySo is None:
1940
        sollya_lib_clear_obj(qualitySo)
1941
    return(remezPolySo)
1942
# End pobyso_remez_exponentsList_sa_so
1943
#
1944
def pobyso_round_coefficients_so_so(polySo, truncFormatListSo):
1945
    """
1946
    A wrapper around the "classical" sollya_lib_roundcoefficients: a Sollya
1947
    polynomial and a Sollya list are given as arguments.
1948
    """
1949
    return sollya_lib_roundcoefficients(polySo, truncFormatListSo)
1950

    
1951
def pobyso_round_coefficients_progressive_so_so(polySo, 
1952
                                                funcSo,
1953
                                                precSo,
1954
                                                intervalSo,
1955
                                                icSo,
1956
                                                currentApproxErrorSo,
1957
                                                approxAccurSo,
1958
                                                debug=False):
1959
    """
1960
    From an input approximation polynomial, compute an output one with 
1961
    smaller coefficients and yet yields a sufficient approximation accuracy.
1962
    """
1963
    if debug:
1964
        print "Input arguments:"
1965
        print "Polynomial: ", ; pobyso_autoprint(polySo)
1966
        print "Function: ", ; pobyso_autoprint(funcSo)
1967
        print "Internal precision: ", ; pobyso_autoprint(precSo)
1968
        print "Interval: ", ; pobyso_autoprint(intervalSo)
1969
        print "Current approximation error: ", ; pobyso_autoprint(currentApproxErrorSo)
1970
        print "Requested approxiation error: ", ; pobyso_autoprint(approxAccurSo)
1971
        print "________________"
1972
    approxAccurSa        = pobyso_get_constant_as_rn_so_sa(approxAccurSo)
1973
    currentApproxErrorSa = pobyso_get_constant_as_rn_so_sa(currentApproxErrorSo)
1974
    ## If the current approximation error is too close to the target, there is
1975
    #  no possible gain.
1976
    if currentApproxErrorSa >= approxAccurSa / 2:
1977
        #### Do not return the initial argument but copies: caller may free 
1978
        #    the former as inutile after call.
1979
        return (sollya_lib_copy_obj(polySo),
1980
                sollya_lib_copy_obj(currentApproxErrorSo))
1981
    #
1982
    ## Try to round the coefficients.
1983
    degreeSa             = pobyso_polynomial_degree_so_sa(polySo)
1984
    intervalSa           = pobyso_range_to_interval_so_sa(intervalSo)
1985
    
1986
    if debug:
1987
        print "degreeSa              :", degreeSa
1988
        print "intervalSa            :", intervalSa.str(style='brackets')
1989
        print "currentApproxErrorSa  :", currentApproxErrorSa 
1990
        print "approxAccurSa         :", approxAccurSa 
1991
    radiusSa = intervalSa.absolute_diameter() / 2
1992
    if debug:
1993
        print "log2(radius):", RR(radiusSa).log2()
1994
    iterIndex = 0
1995
    ## Build the "shaved" polynomial.
1996
    while True: 
1997
        ### Start with a 0 value expression.
1998
        resPolySo = pobyso_constant_0_sa_so()
1999
        roundedPolyApproxAccurSa = approxAccurSa / 2
2000
        currentRadiusPowerSa = 1 
2001
        for degree in xrange(0,degreeSa + 1):
2002
            #### At round 0, use the agressive formula. At round 1, run the
2003
            #    proved formula.
2004
            if iterIndex == 0:
2005
                roundingPowerSa = \
2006
                    floor(((currentRadiusPowerSa/roundedPolyApproxAccurSa)*(degree+1)).log2())
2007
            else:
2008
                roundingPowerSa = \
2009
                    floor(((currentRadiusPowerSa/roundedPolyApproxAccurSa)*(degreeSa+1)).log2())
2010
            ## Under extreme conditions the above formulas can evaluate under 2,
2011
            #   which is the minimal precision of an MPFR number.
2012
            if roundingPowerSa < 2:
2013
                roundingPowerSa = 2
2014
            if debug:
2015
                print "roundedPolyApproxAccurSa", roundedPolyApproxAccurSa
2016
                print "currentRadiusPowerSa", currentRadiusPowerSa
2017
                print "Current rounding exponent:", roundingPowerSa
2018
            currentRadiusPowerSa *= radiusSa
2019
            index1So = pobyso_constant_from_int_sa_so(degree)
2020
            index2So = pobyso_constant_from_int_sa_so(degree)
2021
            ### Create a monomial with:
2022
            #   - the coefficient in the initial monomial at the current degrree;
2023
            #   - the current exponent;
2024
            #   - the free variable.
2025
            cmonSo  = \
2026
                sollya_lib_build_function_mul(sollya_lib_coeff(polySo, index1So),
2027
                                              sollya_lib_build_function_pow( \
2028
                                                sollya_lib_build_function_free_variable(), \
2029
                                                index2So))
2030
            roundingPowerSo = pobyso_constant_from_int_sa_so(roundingPowerSa)
2031
            cmonrSo = pobyso_round_coefficients_single_so_so(cmonSo, roundingPowerSo)
2032
            sollya_lib_clear_obj(cmonSo)
2033
            ### Add to the result polynomial.
2034
            resPolySo = sollya_lib_build_function_add(resPolySo,
2035
                                                      cmonrSo)
2036
        # End for.
2037
        ### Check the new polynomial.
2038
        freeVarSo     = sollya_lib_build_function_free_variable()
2039
        changeVarSo   = sollya_lib_sub(freeVarSo, icSo)
2040
        resPolyCvSo   = sollya_lib_evaluate(resPolySo, changeVarSo)
2041
        errFuncSo     = sollya_lib_build_function_sub(sollya_lib_copy_obj(funcSo), 
2042
                                                      resPolyCvSo)
2043
        infNormSo = sollya_lib_dirtyinfnorm(errFuncSo, intervalSo)
2044
        ### This also clears resPolyCvSo.
2045
        sollya_lib_clear_obj(errFuncSo)
2046
        cerrSa    = pobyso_get_constant_as_rn_so_sa(infNormSo)
2047
        if debug:
2048
            print "Error of the new polynomial:", cerrSa
2049
        ### If at round 1, return the initial polynomial error. This should
2050
        #   never happen since the rounding algorithm is proved. But some 
2051
        #   circumstances may break it (e.g. internal precision of tools).
2052
        if cerrSa > approxAccurSa:
2053
            if iterIndex == 0: # Round 0 is agressive rounding, got round 1 (proved rounding)
2054
                sollya_lib_clear_obj(resPolySo)
2055
                sollya_lib_clear_obj(infNormSo)
2056
                iterIndex += 1
2057
                continue
2058
            else: # Round 1 and beyond : just return the oroginal polynomial.
2059
                sollya_lib_clear_obj(resPolySo)
2060
                sollya_lib_clear_obj(infNormSo)
2061
                #### Do not return the arguments but copies: the caller may free
2062
                #    free the former as inutile after call.
2063
                return (sollya_lib_copy_obj(polySo), 
2064
                        sollya_lib_copy_obj(currentApproxErrorSo))
2065
        ### If get here it is because cerrSa <= approxAccurSa
2066
        ### Approximation error of the new polynomial is acceptable.
2067
        return (resPolySo, infNormSo)
2068
    # End while True
2069
# End pobyso_round_coefficients_progressive_so_so
2070
    
2071
def pobyso_round_coefficients_single_so_so(polySo, commonPrecSo):
2072
    """
2073
    Create a rounded coefficients polynomial from polynomial argument to
2074
    the number of bits in size argument.
2075
    All coefficients are set to the same precision.
2076
    """
2077
    ## TODO: check arguments.
2078
    endEllipListSo = pobyso_build_end_elliptic_list_so_so(commonPrecSo)
2079
    polySo = sollya_lib_roundcoefficients(polySo, endEllipListSo, None)
2080
    sollya_lib_clear_obj(endEllipListSo)
2081
    #sollya_lib_clear_obj(endEllipListSo)
2082
    return polySo
2083
    
2084
# End pobyso_round_coefficients_single_so_so
2085

    
2086
def pobyso_set_canonical_off():
2087
    sollya_lib_set_canonical(sollya_lib_off())
2088

    
2089
def pobyso_set_canonical_on():
2090
    sollya_lib_set_canonical(sollya_lib_on())
2091

    
2092
def pobyso_set_prec(p):
2093
    """ Legacy function. See pobyso_set_prec_sa_so. """
2094
    pobyso_set_prec_sa_so(p)
2095

    
2096
def pobyso_set_prec_sa_so(p):
2097
    #a = c_int(p)
2098
    #precSo = c_void_p(sollya_lib_constant_from_int(a))
2099
    #precSo = sollya_lib_constant_from_int(a)
2100
    precSo = pobyso_constant_from_int_sa_so(p)
2101
    sollya_lib_set_prec(precSo)
2102
    sollya_lib_clear_obj(precSo)
2103
# End pobyso_set_prec_sa_so.
2104

    
2105
def pobyso_set_prec_so_so(newPrecSo):
2106
    sollya_lib_set_prec(newPrecSo)
2107
# End pobyso_set_prec_so_so.
2108
#   
2109
def pobyso_supnorm_sa_sa(polySa, 
2110
                         funcSa, 
2111
                         intervalSa, 
2112
                         errorTypeSa=SOLLYA_ABSOLUTE, 
2113
                         accuracySa=RR(2^-40)):
2114
    """
2115
    An supnorm call with Sage arguments.
2116
    """
2117
    # Check that polySa is a univariate polynomial. We only check here at
2118
    # expression level, not at polynomial ring level
2119
    ## Make sure it is not a multivariate polynomial. The number of variables
2120
    #  can be zero in the case of a constant.
2121
    try:
2122
        if len(polySa.free_variables()) > 1 :
2123
            print "Invalid number of free variables in polySa:", polySa,
2124
            " -", polySa.free_variables()
2125
            return None
2126
    except AttributeError:
2127
        print "polySa is not a SymbolicExpression."
2128
        return None 
2129
    ## Make sure it is a polynomial.
2130
    if not polySa.is_polynomial(polySa.default_variable()):
2131
        print "polySa is not a polynomila expression."
2132
        return None
2133
    # Check that funcSa is a function.
2134
    if not \
2135
     sage.symbolic.callable.is_CallableSymbolicExpressionRing(parent(funcSa)):
2136
        return None
2137
    # Check that intervalSa is an interval.
2138
    try:
2139
        intervalSa.upper()
2140
    except AttributeError:
2141
        print "intervalSa is not an interval."
2142
        return None
2143
    # Convert the Sage polynomial into a Sollya polynomial.
2144
    polyAsStringSa = polySa._assume_str().replace('_SAGE_VAR_', '')
2145
    polySo = pobyso_parse_string_sa_so(polyAsStringSa)
2146
    if not pobyso_obj_is_function_so_sa(polySo):
2147
        sollya_lib_clear_obj(polySo)
2148
        print "The Sollya object created from the polynomial is not a function."
2149
        return None
2150
    #pobyso_autoprint(polySo)
2151
    # Convert the Sage function into a Sollya function.
2152
    funcAsStringSa = funcSa._assume_str().replace('_SAGE_VAR_', '')
2153
    funcSo = pobyso_parse_string_sa_so(funcAsStringSa)
2154
    if not pobyso_obj_is_function_so_sa(funcSo):
2155
        sollya_lib_clear_obj(polySo)
2156
        sollya_lib_clear_obj(funcSo)
2157
        print "The Sollya object created from the function is not a function."
2158
        return None
2159
    #pobyso_autoprint(funcSo)
2160
    # Convert the Sage interval into a Sollya range.
2161
    rangeSo = pobyso_interval_to_range_sa_so(intervalSa)
2162
    if not pobyso_is_range_so_sa(rangeSo):
2163
        sollya_lib_clear_obj(polySo)
2164
        sollya_lib_clear_obj(funcSo)
2165
        sollya_lib_clear_obj(rangeSo)
2166
        print "The Sollya object created from the interval is not a range."
2167
        return None
2168
    #pobyso_autoprint(rangeSo)
2169
    # Check the error type. We do not check the returned object since we
2170
    # trust our code and Sollya on this one.
2171
    if errorTypeSa != SOLLYA_ABSOLUTE and errorTypeSa != SOLLYA_RELATIVE:
2172
        sollya_lib_clear_obj(polySo)
2173
        sollya_lib_clear_obj(funcSo)
2174
        sollya_lib_clear_obj(rangeSo)
2175
        return None
2176
    if errorTypeSa == SOLLYA_ABSOLUTE:
2177
        errorTypeSo = pobyso_absolute_so_so()
2178
    else:
2179
        errorTypeSo = pobyso_relative_so_so()
2180
    #pobyso_autoprint(errorTypeSo)
2181
    # Check if accuracySa is an element of a RealField. If not, try to 
2182
    # convert it.
2183
    try:
2184
        accuracySa.ulp()
2185
    except AttributeError:
2186
        try:
2187
            accuracySa = RR(accuracySa)
2188
        except TypeError:
2189
            accuracySa = RR(str(accuracySa))
2190
    # Create the Sollya object
2191
    accuracySo = pobyso_constant_sa_so(accuracySa)
2192
    #pobyso_autoprint(accuracySo)
2193
    if pobyso_is_error_so_sa(accuracySo):
2194
        sollya_lib_clear_obj(polySo)
2195
        sollya_lib_clear_obj(funcSo)
2196
        sollya_lib_clear_obj(rangeSo)
2197
        sollya_lib_clear_obj(errorTypeSo)
2198
        sollya_lib_clear_obj(accuracySo)
2199
        return None
2200
    retValSo = sollya_lib_supnorm(polySo, 
2201
                                  funcSo, 
2202
                                  rangeSo, 
2203
                                  errorTypeSo, 
2204
                                  accuracySo, 
2205
                                  None)
2206
    sollya_lib_clear_obj(polySo)
2207
    sollya_lib_clear_obj(funcSo)
2208
    sollya_lib_clear_obj(rangeSo)
2209
    sollya_lib_clear_obj(errorTypeSo)
2210
    sollya_lib_clear_obj(accuracySo)
2211
    if pobyso_is_error_so_sa(retValSo):
2212
        sollya_lib_clear_obj(retValSo)
2213
        return None
2214
    #pobyso_autoprint(retValSo)
2215
    retValSa = pobyso_range_to_interval_so_sa(retValSo)
2216
    sollya_lib_clear_obj(retValSo)
2217
    return retValSa
2218
# End pobyso_supnorm_sa_sa
2219

    
2220
def pobyso_supnorm_so_sa(polySo, funcSo, intervalSo, errorTypeSo = None,\
2221
                         accuracySo = None, realFieldSa = None):
2222
    """
2223
    Computes the supremum norm from Sollya input arguments and returns a
2224
    Sage floating-point number whose precision is set by the only Sage argument.
2225
    
2226
    The returned value is the maximum of the absolute values of the range
2227
    elements  returned by the Sollya supnorm functions
2228
    """
2229
    supNormRangeSo = pobyso_supnorm_so_so(polySo, 
2230
                                          funcSo, 
2231
                                          intervalSo, 
2232
                                          errorTypeSo,
2233
                                          accuracySo)
2234
    supNormSo = pobyso_range_max_abs_so_so(supNormRangeSo)
2235
    sollya_lib_clear_obj(supNormRangeSo)
2236
    #pobyso_autoprint(supNormSo)
2237
    supNormSa = pobyso_get_constant_as_rn_with_rf_so_sa(supNormSo, realFieldSa)
2238
    sollya_lib_clear_obj(supNormSo)
2239
    return supNormSa 
2240
# End pobyso_supnorm_so_sa.
2241
#
2242
def pobyso_supnorm_so_so(polySo, funcSo, intervalSo, errorTypeSo = None,\
2243
                         accuracySo = None):
2244
    """
2245
    Computes the supnorm of the approximation error between the given 
2246
    polynomial and function. Attention: returns a range!
2247
    errorTypeSo defaults to "absolute".
2248
    accuracySo defaults to 2^(-40).
2249
    """
2250
    if errorTypeSo is None:
2251
        errorTypeSo = sollya_lib_absolute(None)
2252
        errorTypeIsNone = True
2253
    else:
2254
        errorTypeIsNone = False
2255
    #
2256
    if accuracySo is None:
2257
        # Notice the **: we are in Pythonland here!
2258
        accuracySo = pobyso_constant_sa_so(RR(2**(-40)))
2259
        accuracyIsNone = True
2260
    else:
2261
        accuracyIsNone = False
2262
    #pobyso_autoprint(accuracySo)
2263
    resultSo = \
2264
        sollya_lib_supnorm(polySo, funcSo, intervalSo, errorTypeSo, \
2265
                              accuracySo)
2266
    if errorTypeIsNone:
2267
        sollya_lib_clear_obj(errorTypeSo)
2268
    if accuracyIsNone:
2269
        sollya_lib_clear_obj(accuracySo)
2270
    return resultSo
2271
# End pobyso_supnorm_so_so
2272
#
2273
def pobyso_taylor_expansion_no_change_var_so_so(functionSo, 
2274
                                                degreeSo, 
2275
                                                rangeSo,
2276
                                                errorTypeSo=None, 
2277
                                                sollyaPrecSo=None):
2278
    """
2279
    Compute the Taylor expansion without the variable change
2280
    x -> x-intervalCenter.
2281
    If errorTypeSo is None, absolute is used.
2282
    If sollyaPrecSo is None, Sollya internal precision is not changed. 
2283
    """
2284
    # Change internal Sollya precision, if needed.
2285
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
2286
    sollyaPrecChanged = False
2287
    if sollyaPrecSo is None:
2288
        pass
2289
    else:
2290
        sollya_lib_set_prec(sollyaPrecSo)
2291
        sollyaPrecChanged = True
2292
    # Error type stuff: default to absolute.
2293
    if errorTypeSo is None:
2294
        errorTypeIsNone = True
2295
        errorTypeSo = sollya_lib_absolute(None)
2296
    else:
2297
        errorTypeIsNone = False
2298
    intervalCenterSo = sollya_lib_mid(rangeSo, None)
2299
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo,
2300
                                         intervalCenterSo,
2301
                                         rangeSo, errorTypeSo, None)
2302
    # Object taylorFormListSaSo is a Python list of Sollya objects references 
2303
    #  that are copies of the elements of taylorFormSo.
2304
    # pobyso_get_list_elements_so_so clears taylorFormSo.
2305
    (taylorFormListSaSo, numElementsSa, isEndEllipticSa) = \
2306
        pobyso_get_list_elements_so_so(taylorFormSo)
2307
    ## Copy needed here since polySo will be returned and taylorFormListSaSo
2308
    #  will be cleared.
2309
    polySo = sollya_lib_copy_obj(taylorFormListSaSo[0])
2310
    #print "Num elements:", numElementsSa
2311
    sollya_lib_clear_obj(taylorFormSo)
2312
    # No copy_obj needed here: a new objects are created.
2313
    maxErrorSo    = sollya_lib_sup(taylorFormListSaSo[2])
2314
    minErrorSo    = sollya_lib_inf(taylorFormListSaSo[2])
2315
    # List taylorFormListSaSo is not needed anymore.
2316
    pobyso_clear_list_elements_sa_so(taylorFormListSaSo)
2317
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
2318
    absMinErrorSo = sollya_lib_abs(minErrorSo)
2319
    sollya_lib_clear_obj(maxErrorSo)
2320
    sollya_lib_clear_obj(minErrorSo)
2321
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
2322
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
2323
    #
2324
    if errorTypeIsNone:
2325
        sollya_lib_clear_obj(errorTypeSo)
2326
    ## If changed, reset the Sollya working precision.
2327
    if sollyaPrecChanged:
2328
        sollya_lib_set_prec(initialSollyaPrecSo)
2329
    sollya_lib_clear_obj(initialSollyaPrecSo)
2330
    ## According to what error is the largest, return the errors.
2331
    if absMaxErrorSa > absMinErrorSa:
2332
        sollya_lib_clear_obj(absMinErrorSo)
2333
        return (polySo, intervalCenterSo, absMaxErrorSo)
2334
    else:
2335
        sollya_lib_clear_obj(absMaxErrorSo)
2336
        return (polySo, intervalCenterSo, absMinErrorSo)
2337
# end pobyso_taylor_expansion_no_change_var_so_so
2338

    
2339
def pobyso_taylor_expansion_with_change_var_so_so(functionSo, degreeSo, \
2340
                                                  rangeSo, \
2341
                                                  errorTypeSo=None, \
2342
                                                  sollyaPrecSo=None):
2343
    """
2344
    Compute the Taylor expansion with the variable change
2345
    x -> (x-intervalCenter) included.
2346
    """
2347
    # Change Sollya internal precision, if need.
2348
    sollyaPrecChanged = False
2349
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
2350
    if sollyaPrecSo is None:
2351
        pass
2352
    else:
2353
        sollya_lib_set_prec(sollyaPrecSo)
2354
        sollyaPrecChanged = True
2355
    #
2356
    # Error type stuff: default to absolute.
2357
    if errorTypeSo is None:
2358
        errorTypeIsNone = True
2359
        errorTypeSo = sollya_lib_absolute(None)
2360
    else:
2361
        errorTypeIsNone = False
2362
    intervalCenterSo = sollya_lib_mid(rangeSo)
2363
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo, \
2364
                                         intervalCenterSo, \
2365
                                         rangeSo, errorTypeSo, None)
2366
    # Object taylorFormListSaSo is a Python list of Sollya objects references 
2367
    # that are copies of the elements of taylorFormSo.
2368
    # pobyso_get_list_elements_so_so clears taylorFormSo.
2369
    (taylorFormListSaSo, numElements, isEndElliptic) = \
2370
        pobyso_get_list_elements_so_so(taylorFormSo)
2371
    sollya_lib_clear_obj(taylorFormSo)
2372
    polySo = taylorFormListSaSo[0]
2373
    ## Maximum error computation with taylorFormListSaSo[2], a range
2374
    #  holding the actual error. Bounds can be negative.
2375
    maxErrorSo    = sollya_lib_sup(taylorFormListSaSo[2])
2376
    minErrorSo    = sollya_lib_inf(taylorFormListSaSo[2])
2377
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
2378
    absMinErrorSo = sollya_lib_abs(minErrorSo)
2379
    sollya_lib_clear_obj(maxErrorSo)
2380
    sollya_lib_clear_obj(minErrorSo)
2381
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
2382
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
2383
    changeVarExpSo = sollya_lib_build_function_sub(\
2384
                       sollya_lib_build_function_free_variable(),\
2385
                       sollya_lib_copy_obj(intervalCenterSo))
2386
    polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpSo)
2387
    # List taylorFormListSaSo is not needed anymore.
2388
    pobyso_clear_list_elements_sa_so(taylorFormListSaSo)
2389
    sollya_lib_clear_obj(changeVarExpSo)
2390
    # If changed, reset the Sollya working precision.
2391
    if sollyaPrecChanged:
2392
        sollya_lib_set_prec(initialSollyaPrecSo)
2393
    sollya_lib_clear_obj(initialSollyaPrecSo)
2394
    if errorTypeIsNone:
2395
        sollya_lib_clear_obj(errorTypeSo)
2396
    # Do not clear maxErrorSo.
2397
    if absMaxErrorSa > absMinErrorSa:
2398
        sollya_lib_clear_obj(absMinErrorSo)
2399
        return((polyVarChangedSo, intervalCenterSo, absMaxErrorSo))
2400
    else:
2401
        sollya_lib_clear_obj(absMaxErrorSo)
2402
        return((polyVarChangedSo, intervalCenterSo, absMinErrorSo))
2403
# end pobyso_taylor_expansion_with_change_var_so_so
2404

    
2405
def pobyso_taylor(function, degree, point):
2406
    """ Legacy function. See pobysoTaylor_so_so. """
2407
    return(pobyso_taylor_so_so(function, degree, point))
2408

    
2409
def pobyso_taylor_so_so(functionSo, degreeSo, pointSo):
2410
    return(sollya_lib_taylor(functionSo, degreeSo, pointSo))
2411
    
2412
def pobyso_taylorform(function, degree, point = None, 
2413
                      interval = None, errorType=None):
2414
    """ Legacy function. See pobyso_taylorform_sa_sa;"""
2415
    
2416
def pobyso_taylorform_sa_sa(functionSa, \
2417
                            degreeSa, \
2418
                            pointSa, \
2419
                            intervalSa=None, \
2420
                            errorTypeSa=None, \
2421
                            precisionSa=None):
2422
    """
2423
    Compute the Taylor form of 'degreeSa' for 'functionSa' at 'pointSa' 
2424
    for 'intervalSa' with 'errorTypeSa' (a string) using 'precisionSa'. 
2425
    point: must be a Real or a Real interval.
2426
    return the Taylor form as an array
2427
    TODO: take care of the interval and of the point when it is an interval;
2428
          when errorType is not None;
2429
          take care of the other elements of the Taylor form (coefficients 
2430
          errors and delta.
2431
    """
2432
    # Absolute as the default error.
2433
    if errorTypeSa is None:
2434
        errorTypeSo = sollya_lib_absolute()
2435
    elif errorTypeSa == "relative":
2436
        errorTypeSo = sollya_lib_relative()
2437
    elif errortypeSa == "absolute":
2438
        errorTypeSo = sollya_lib_absolute()
2439
    else:
2440
        # No clean up needed.
2441
        return None
2442
    # Global precision stuff
2443
    sollyaPrecisionChangedSa = False
2444
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
2445
    if precisionSa is None:
2446
        precSa = initialSollyaPrecSa
2447
    else:
2448
        if precSa > initialSollyaPrecSa:
2449
            if precSa <= 2:
2450
                print inspect.stack()[0][3], ":precision change <= 2 requested."
2451
            pobyso_set_prec_sa_so(precSa)
2452
            sollyaPrecisionChangedSa = True
2453
    #        
2454
    if len(functionSa.variables()) > 0:
2455
        varSa = functionSa.variables()[0]
2456
        pobyso_name_free_variable_sa_so(str(varSa))
2457
    # In any case (point or interval) the parent of pointSa has a precision
2458
    # method.
2459
    pointPrecSa = pointSa.parent().precision()
2460
    if precSa > pointPrecSa:
2461
        pointPrecSa = precSa
2462
    # In any case (point or interval) pointSa has a base_ring() method.
2463
    pointBaseRingString = str(pointSa.base_ring())
2464
    if re.search('Interval', pointBaseRingString) is None: # Point
2465
        pointSo = pobyso_constant_sa_so(pointSa, pointPrecSa)
2466
    else: # Interval.
2467
        pointSo = pobyso_interval_to_range_sa_so(pointSa, pointPrecSa)
2468
    # Sollyafy the function.
2469
    functionSo = pobyso_parse_string_sa_so(functionSa._assume_str().replace('_SAGE_VAR_', ''))
2470
    if sollya_lib_obj_is_error(functionSo):
2471
        print "pobyso_tailorform: function string can't be parsed!"
2472
        return None
2473
    # Sollyafy the degree
2474
    degreeSo = sollya_lib_constant_from_int(int(degreeSa))
2475
    # Sollyafy the point
2476
    # Call Sollya
2477
    taylorFormSo = \
2478
        sollya_lib_taylorform(functionSo, degreeSo, pointSo, errorTypeSo,\
2479
                                         None)
2480
    sollya_lib_clear_obj(functionSo)
2481
    sollya_lib_clear_obj(degreeSo)
2482
    sollya_lib_clear_obj(pointSo)
2483
    sollya_lib_clear_obj(errorTypeSo)
2484
    (tfsAsList, numElements, isEndElliptic) = \
2485
            pobyso_get_list_elements_so_so(taylorFormSo)
2486
    polySo = tfsAsList[0]
2487
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
2488
    polyRealField = RealField(maxPrecision)
2489
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, polyRealField)
2490
    if sollyaPrecisionChangedSa:
2491
        sollya_lib_set_prec(initialSollyaPrecSo)
2492
    sollya_lib_clear_obj(initialSollyaPrecSo)
2493
    polynomialRing = polyRealField[str(varSa)]
2494
    polySa = polynomial(expSa, polynomialRing)
2495
    taylorFormSa = [polySa]
2496
    # Final clean-up.
2497
    sollya_lib_clear_obj(taylorFormSo)
2498
    return(taylorFormSa)
2499
# End pobyso_taylor_form_sa_sa
2500

    
2501
def pobyso_taylorform_so_so(functionSo, degreeSo, pointSo, intervalSo=None, \
2502
                            errorTypeSo=None):
2503
    createdErrorType = False
2504
    if errorTypeSo is None:
2505
        errorTypeSo = sollya_lib_absolute()
2506
        createdErrorType = True
2507
    else:
2508
        #TODO: deal with the other case.
2509
        pass
2510
    if intervalSo is None:
2511
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
2512
                                         errorTypeSo, None)
2513
    else:
2514
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
2515
                                         intervalSo, errorTypeSo, None)
2516
    if createdErrorType:
2517
        sollya_lib_clear_obj(errorTypeSo)
2518
    return resultSo
2519
        
2520

    
2521
def pobyso_univar_polynomial_print_reverse(polySa):
2522
    """ Legacy function. See pobyso_univar_polynomial_print_reverse_sa_sa. """
2523
    return(pobyso_univar_polynomial_print_reverse_sa_sa(polySa))
2524

    
2525
def pobyso_univar_polynomial_print_reverse_sa_sa(polySa):
2526
    """
2527
    Return the string representation of a univariate polynomial with
2528
    monomials ordered in the x^0..x^n order of the monomials.
2529
    Remember: Sage
2530
    """
2531
    polynomialRing = polySa.base_ring()
2532
    # A very expensive solution:
2533
    # -create a fake multivariate polynomial field with only one variable,
2534
    #   specifying a negative lexicographical order;
2535
    mpolynomialRing = PolynomialRing(polynomialRing.base(), \
2536
                                     polynomialRing.variable_name(), \
2537
                                     1, order='neglex')
2538
    # - convert the univariate argument polynomial into a multivariate
2539
    #   version;
2540
    p = mpolynomialRing(polySa)
2541
    # - return the string representation of the converted form.
2542
    # There is no simple str() method defined for p's class.
2543
    return(p.__str__())
2544
#
2545
#print pobyso_get_prec()  
2546
pobyso_set_prec(165)
2547
#print pobyso_get_prec()  
2548
#a=100
2549
#print type(a)
2550
#id(a)
2551
#print "Max arity: ", pobyso_max_arity
2552
#print "Function tripleDouble (43) as a string: ", pobyso_function_type_as_string(43)
2553
#print "Function None (44) as a string: ", pobyso_function_type_as_string(44)
2554
sys.stderr.write("\t...Pobyso check done.\n")