Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 283

Historique | Voir | Annoter | Télécharger (93,91 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.
40
"""
41
(SOLLYA_BASE_FUNC_ABS,
42
SOLLYA_BASE_FUNC_ACOS,
43
    SOLLYA_BASE_FUNC_ACOSH,
44
    SOLLYA_BASE_FUNC_ADD,
45
    SOLLYA_BASE_FUNC_ASIN,
46
    SOLLYA_BASE_FUNC_ASINH,
47
    SOLLYA_BASE_FUNC_ATAN,
48
    SOLLYA_BASE_FUNC_ATANH,
49
    SOLLYA_BASE_FUNC_CEIL,
50
    SOLLYA_BASE_FUNC_CONSTANT,
51
    SOLLYA_BASE_FUNC_COS,
52
    SOLLYA_BASE_FUNC_COSH,
53
    SOLLYA_BASE_FUNC_DIV,
54
    SOLLYA_BASE_FUNC_DOUBLE,
55
    SOLLYA_BASE_FUNC_DOUBLEDOUBLE,
56
    SOLLYA_BASE_FUNC_DOUBLEEXTENDED,
57
    SOLLYA_BASE_FUNC_ERF,
58
    SOLLYA_BASE_FUNC_ERFC,
59
    SOLLYA_BASE_FUNC_EXP,
60
    SOLLYA_BASE_FUNC_EXP_M1,
61
    SOLLYA_BASE_FUNC_FLOOR,
62
    SOLLYA_BASE_FUNC_FREE_VARIABLE,
63
    SOLLYA_BASE_FUNC_HALFPRECISION,
64
    SOLLYA_BASE_FUNC_LIBRARYCONSTANT,
65
    SOLLYA_BASE_FUNC_LIBRARYFUNCTION,
66
    SOLLYA_BASE_FUNC_LOG,
67
    SOLLYA_BASE_FUNC_LOG_10,
68
    SOLLYA_BASE_FUNC_LOG_1P,
69
    SOLLYA_BASE_FUNC_LOG_2,
70
    SOLLYA_BASE_FUNC_MUL,
71
    SOLLYA_BASE_FUNC_NEARESTINT,
72
    SOLLYA_BASE_FUNC_NEG,
73
    SOLLYA_BASE_FUNC_PI,
74
    SOLLYA_BASE_FUNC_POW,
75
    SOLLYA_BASE_FUNC_PROCEDUREFUNCTION,
76
    SOLLYA_BASE_FUNC_QUAD,
77
    SOLLYA_BASE_FUNC_SIN,
78
    SOLLYA_BASE_FUNC_SINGLE,
79
    SOLLYA_BASE_FUNC_SINH,
80
    SOLLYA_BASE_FUNC_SQRT,
81
    SOLLYA_BASE_FUNC_SUB,
82
    SOLLYA_BASE_FUNC_TAN,
83
    SOLLYA_BASE_FUNC_TANH,
84
SOLLYA_BASE_FUNC_TRIPLEDOUBLE) = map(int,xrange(44))
85
sys.stderr.write("Superficial pobyso check...\n")
86
#print "First constant - SOLLYA_BASE_FUNC_ABS: ", SOLLYA_BASE_FUNC_ABS
87
#print "Last constant  - SOLLYA_BASE_FUNC_TRIPLEDOUBLE: ", SOLLYA_BASE_FUNC_TRIPLEDOUBLE
88

    
89
pobyso_max_arity = 9
90

    
91
def pobyso_absolute_so_so():
92
    """
93
    Create an "absolute" Sollya object.
94
    """
95
    return(sollya_lib_absolute(None))
96

    
97
def pobyso_autoprint(arg):
98
    sollya_lib_autoprint(arg, None)
99

    
100
def pobyso_autoprint_so_so(arg):
101
    sollya_lib_autoprint(arg,None)
102

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

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

    
186
def pobyso_build_function_sub_so_so(exp1So, exp2So):
187
    return sollya_lib_build_function_sub(exp1So, exp2So)
188

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

    
224
def pobyso_change_var_in_function_so_so(funcSo, chvarExpSo):
225
    """
226
    Variable change in a function.
227
    """
228
    return(sollya_lib_evaluate(funcSo,chvarExpSo))
229
# End pobyso_change_var_in_function_so_so     
230

    
231
def pobyso_chebyshevform_so_so(functionSo, degreeSo, intervalSo):
232
    resultSo = sollya_lib_chebyshevform(functionSo, degreeSo, intervalSo)
233
    return(resultSo)
234
# End pobyso_chebyshevform_so_so.
235

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

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

    
259
def pobyso_clear_obj(objSo):
260
    """
261
    Free a Sollya object's memory.
262
    Very thin wrapper around sollya_lib_clear_obj().
263
    """
264
    sollya_lib_clear_obj(objSo)
265
# End pobyso_clear_obj. 
266

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

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

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

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

    
367

    
368
def pobyso_constant(rnArg):
369
    """ Legacy function. See pobyso_constant_sa_so. """
370
    return(pobyso_constant_sa_so(rnArg))
371
    
372
def pobyso_constant_sa_so(rnArgSa, precisionSa=None):
373
    """
374
    Create a Sollya constant from a Sage RealNumber.
375
    The sollya_lib_constant() function creates a constant
376
    with the same precision as the source.
377
    """
378
    ## Precision stuff. If one wants to change precisions,
379
    #  everything takes place in Sage. That only makes
380
    #  sense if one wants to reduce the precision.
381
    # TODO: revisit precision stuff with new technique.
382
    if not precisionSa is None:
383
        RRR = RealField(precisionSa)
384
        rnArgSa = RRR(rnArgSa)
385
    #print rnArgSa, rnArgSa.precision()
386
    # Sollya constant creation takes place here.
387
    return sollya_lib_constant(get_rn_value(rnArgSa))
388
# End pobyso_constant_sa_so
389
     
390
def pobyso_constant_0_sa_so():
391
    """
392
    Obvious.
393
    """
394
    return pobyso_constant_from_int_sa_so(0)
395

    
396
def pobyso_constant_1():
397
    """
398
    Obvious.
399
    Legacy function. See pobyso_constant_so_so. 
400
    """
401
    return pobyso_constant_1_sa_so()
402

    
403
def pobyso_constant_1_sa_so():
404
    """
405
    Obvious.
406
    """
407
    return(pobyso_constant_from_int_sa_so(1))
408

    
409
def pobyso_constant_from_int(anInt):
410
    """ Legacy function. See pobyso_constant_from_int_sa_so. """
411
    return pobyso_constant_from_int_sa_so(anInt)
412

    
413
def pobyso_constant_from_int_sa_so(anInt):
414
    """
415
    Get a Sollya constant from a Sage int.
416
    """
417
    return sollya_lib_constant_from_int64(long(anInt))
418

    
419
def pobyso_constant_from_int_so_sa(constSo):
420
    """
421
    Get a Sage int from a Sollya int constant.
422
    Usefull for precision or powers in polynomials.
423
    """
424
    constSa = c_long(0)
425
    sollya_lib_get_constant_as_int64(byref(constSa), constSo)
426
    return constSa.value
427
# End pobyso_constant_from_int_so_sa
428

    
429
def pobyso_constant_from_mpq_sa_so(rationalSa):
430
    """
431
    Make a Sollya constant from Sage rational.
432
    The Sollya constant is an unevaluated expression.
433
    Hence no precision argument is needed.
434
    It is better to leave this way since Sollya has its own
435
    optimized evaluation mecanism that tries very hard to
436
    return exact values or at least faithful ones.
437
    """
438
    ratExprSo = \
439
        sollya_lib_constant_from_mpq(sgmp_get_rational_value(rationalSa))
440
    return ratExprSo
441
# End pobyso_constant_from_mpq_sa_so.
442

    
443
def pobyso_constant_sollya_prec_sa_so(rnArgSa):
444
    """
445
    Create a Sollya constant from a Sage RealNumber at the
446
    current precision in Sollya.
447
    """
448
    currentSollyaPrecSa = pobyso_get_prec_so_sa()
449
    return pobyso_constant_sa_so(rnArgSa, currentSollyaPrecSa)
450
# End pobyso_constant_sollya_prec_sa_so
451

    
452
def pobyso_end_elliptic_list_so_sa_so(objectsListSo, intCountSa):
453
    """
454
    Create a Sollya end elliptic list made of the objectListSo[0] to
455
     objectsListSo[intCountSa-1] objects.
456
    """     
457
    return sollya_lib_end_elliptic_list(objectSo, int(intCountSa))
458

    
459
def pobyso_error_so():
460
    return sollya_lib_error(None)
461
# End pobyso_error().
462

    
463
def pobyso_evaluate_so_sa(funcSo, argumentSo):
464
    """
465
    Evaluates funcSo for arguemntSo through sollya_lib_evaluate() and return
466
    the result as a Sage object
467
    """
468
    evalSo = sollya_lib_evaluate(funcSo, argumentSo)
469
    if pobyso_is_error_so_sa(evalSo):
470
        return None
471
    if pobyso_is_range_so_sa(evalSo):
472
        retVal = pobyso_range_to_interval_so_sa(evalSo)
473
    else:
474
        retVal = pobyso_get_constant_as_rn(evalSo)
475
    sollya_lib_clear_obj(evalSo)
476
    return retVal
477
# End pobyso_evaluate_so_sa.
478

    
479
def pobyso_evaluate_so_so(funcSo, argumentSo):
480
    """
481
    Evaluates funcSo for arguemntSo through sollya_lib_evaluate().
482
    """
483
    return sollya_lib_evaluate(funcSo, argumentSo)
484
# End pobyso_evaluate_so_so.
485
#
486
def pobyso_diff_so_so(funcSo):
487
    """
488
    Very thin wrapper around sollya_lib_diff.
489
    """
490
    ## TODO: add a check to make sure funcSo is a functional expression.
491
    return sollya_lib_diff(funcSo)
492

    
493
def pobyso_dirty_find_zeros_so_so(funcSo, rangeSo):
494
    """
495
    Thin wrapper over sollya_lib_dirtyfindzeros()
496
    """
497
    return sollya_lib_dirtyfindzeros(funcSo, rangeSo)
498
# End pobys_dirty_find_zeros
499

    
500
def pobyso_dirty_inf_norm_so_so(funcSo, rangeSo, preSo=None):
501
    """
502
    Thin wrapper around sollya_dirtyinfnorm().
503
    """
504
    # TODO: manage the precision change.
505
    
506
    return sollya_lib_dirtyinfnorm(funcSo, rangeSo)
507
# End pobyso_dirty_inf_norm_so_so
508

    
509
def pobyso_find_zeros_so_so(funcSo, rangeSo):
510
    """
511
    Thin wrapper over sollya_lib_findzeros()
512
    """
513
    return sollya_lib_findzeros(funcSo, rangeSo)
514
# End pobys_find_zeros
515

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

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

    
613
def pobyso_float_poly_so_sa(polySo, realFieldSa=None):
614
    """
615
    Convert a Sollya polynomial into a Sage floating-point polynomial.
616
    If no realField is given, a RealField corresponding to the maximum 
617
    precision of the coefficients is internally computed.
618
    The real field is not returned but can be easily retrieved from 
619
    the polynomial itself.
620
    ALGORITHM:
621
    - (optional) compute the RealField of the coefficients;
622
    - convert the Sollya expression into a Sage expression;
623
    - convert the Sage expression into a Sage polynomial
624
    """    
625
    if realFieldSa is None:
626
        expressionPrecSa = pobyso_get_max_prec_of_exp_so_sa(polySo)
627
        #print "Maximum precision of Sollya polynomial coefficients:", expressionPrecSa
628
        if expressionPrecSa < 2 or expressionPrecSa > 2147483391:
629
            print "Maximum degree of expression:", expressionPrecSa
630
        realFieldSa      = RealField(expressionPrecSa)
631
    #print "Sollya expression before...",
632
    #pobyso_autoprint(polySo)
633

    
634
    expressionSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo,
635
                                                             realFieldSa)
636
    #print "...Sollya expression after."
637
    #pobyso_autoprint(polySo)
638
    polyVariableSa = expressionSa.variables()[0]
639
    polyRingSa     = realFieldSa[str(polyVariableSa)]
640
    #print polyRingSa
641
    # Do not use the polynomial(expressionSa, ring=polyRingSa) form!
642
    polynomialSa = polyRingSa(expressionSa)
643
    polyCoeffsListSa = polynomialSa.coefficients()
644
    #for coeff in polyCoeffsListSa:
645
    #    print coeff.abs().n()
646
    return polynomialSa
647
# End pobyso_float_poly_so_sa
648

    
649
def pobyso_free_variable():
650
    """
651
    Ultra thin wrapper around the sollya_lib_function_build_free_variable function.
652
    """
653
    return sollya_lib_build_function_free_variable()
654
    
655
def pobyso_function_type_as_string(funcType):
656
    """ Legacy function. See pobyso_function_type_as_string_so_sa. """
657
    return(pobyso_function_type_as_string_so_sa(funcType))
658

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

    
755
def pobyso_get_constant(rnArgSa, constSo):
756
    """ Legacy function. See pobyso_get_constant_so_sa. """
757
    return pobyso_get_constant_so_sa(rnArgSa, constSo)
758
# End pobyso_get_constant
759

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

    
802
def pobyso_get_constant_as_rn_with_rf(ctExp, realField):
803
    """ 
804
    Legacy function. See pobyso_get_constant_as_rn_with_rf_so_sa.
805
    """
806
    return pobyso_get_constant_as_rn_with_rf_so_sa(ctExp, realField)
807
# End pobyso_get_constant_as_rn_with_rf
808
    
809
def pobyso_get_constant_as_rn_with_rf_so_sa(ctExpSo, realFieldSa = None):
810
    """
811
    Get a Sollya constant as a Sage "real number".
812
    If no real field is specified, the precision of the floating-point number 
813
    returned is that of the Sollya constant.
814
    Otherwise is is that of the real field. Hence rounding may happen.
815
    """
816
    if realFieldSa is None:
817
        return pobyso_get_constant_as_rn_so_sa(ctExpSo)
818
    rnSa = realFieldSa(0)
819
    outcome = sollya_lib_get_constant(get_rn_value(rnSa), ctExpSo)
820
    if outcome == 0:
821
        return None
822
    else:
823
        return rnSa
824
# End pobyso_get_constant_as_rn_with_rf_so_sa
825

    
826
def pobyso_get_free_variable_name():
827
    """ 
828
    Legacy function. See pobyso_get_free_variable_name_so_sa.
829
    """
830
    return(pobyso_get_free_variable_name_so_sa())
831

    
832
def pobyso_get_free_variable_name_so_sa():
833
    return sollya_lib_get_free_variable_name()
834
    
835
def pobyso_get_function_arity(expressionSo):
836
    """ 
837
    Legacy function. See pobyso_get_function_arity_so_sa.
838
    """
839
    return(pobyso_get_function_arity_so_sa(expressionSo))
840

    
841
def pobyso_get_function_arity_so_sa(expressionSo):
842
    arity = c_int(0)
843
    sollya_lib_get_function_arity(byref(arity),expressionSo)
844
    return int(arity.value)
845

    
846
def pobyso_get_head_function(expressionSo):
847
    """ 
848
    Legacy function. See pobyso_get_head_function_so_sa. 
849
    """
850
    return(pobyso_get_head_function_so_sa(expressionSo)) 
851

    
852
def pobyso_get_head_function_so_sa(expressionSo):
853
    functionType = c_int(0)
854
    sollya_lib_get_head_function(byref(functionType), expressionSo)
855
    return int(functionType.value)
856

    
857
def pobyso_get_interval_from_range_so_sa(soRange, realIntervalFieldSa = None ):
858
    """
859
    Return the Sage interval corresponding to the Sollya range argument.
860
    If no reaIntervalField is passed as an argument, the interval bounds are not
861
    rounded: they are elements of RealIntervalField of the "right" precision
862
    to hold all the digits.
863
    """
864
    prec = c_int(0)
865
    if realIntervalFieldSa is None:
866
        retval = sollya_lib_get_prec_of_range(byref(prec), soRange, None)
867
        if retval == 0:
868
            return None
869
        realIntervalFieldSa = RealIntervalField(prec.value)
870
    intervalSa = realIntervalFieldSa(0,0)
871
    retval = \
872
        sollya_lib_get_interval_from_range(get_interval_value(intervalSa),\
873
                                           soRange)
874
    if retval == 0:
875
        return None
876
    return intervalSa
877
# End pobyso_get_interval_from_range_so_sa
878

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

    
922
def pobyso_get_max_prec_of_exp(soExp):
923
    """ Legacy function. See pobyso_get_max_prec_of_exp_so_sa. """
924
    return pobyso_get_max_prec_of_exp_so_sa(soExp)
925

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

    
967
def pobyso_get_min_prec_of_constant_so_sa(constExpSo):
968
    """
969
    Get the minimum precision necessary to represent the value of a Sollya
970
    constant.
971
    MPFR_MIN_PREC and powers of 2 are taken into account.
972
    We assume that constExpSo is a pointer to a Sollay constant expression.
973
    """
974
    constExpAsRnSa = pobyso_get_constant_as_rn_so_sa(constExpSo)
975
    return(min_mpfr_size(get_rn_value(constExpAsRnSa)))
976

    
977
def pobyso_get_poly_so_sa(polySo, realFieldSa=None):
978
    """
979
    Convert a Sollya polynomial into a Sage polynomial.
980
    Legacy function. Use pobyso_float_poly_so_sa() instead.
981
    """    
982
    return pobyso_float_poly_so_sa(polySo,realFieldSa)
983
# End pobyso_get_poly_so_sa
984

    
985
def pobyso_get_prec():
986
    """ Legacy function. See pobyso_get_prec_so_sa(). """
987
    return pobyso_get_prec_so_sa()
988

    
989
def pobyso_get_prec_so():
990
    """
991
    Get the current default precision in Sollya.
992
    The return value is a Sollya object.
993
    Usefull when modifying the precision back and forth by avoiding
994
    extra conversions.
995
    """
996
    return sollya_lib_get_prec(None)
997
    
998
def pobyso_get_prec_so_sa():
999
    """
1000
    Get the current default precision in Sollya.
1001
    The return value is Sage/Python int.
1002
    """
1003
    precSo = sollya_lib_get_prec()
1004
    precSa = pobyso_constant_from_int_so_sa(precSo)
1005
    sollya_lib_clear_obj(precSo)
1006
    return precSa
1007
# End pobyso_get_prec_so_sa.
1008

    
1009
def pobyso_get_prec_so_so_sa():
1010
    """
1011
    Return the current precision both as a Sollya object and a
1012
    Sage integer as hybrid tuple.
1013
    To avoid multiple calls for precision manipulations. 
1014
    """
1015
    precSo = sollya_lib_get_prec()
1016
    precSa = pobyso_constant_from_int_so_sa(precSo)
1017
    return (precSo, int(precSa))
1018
    
1019
def pobyso_get_prec_of_constant(ctExpSo):
1020
    """ Legacy function. See pobyso_get_prec_of_constant_so_sa. """
1021
    return pobyso_get_prec_of_constant_so_sa(ctExpSo)
1022

    
1023
def pobyso_get_prec_of_constant_so_sa(ctExpSo):
1024
    """
1025
    Tries to find a precision to represent ctExpSo without rounding.
1026
    If not possible, returns None.
1027
    """
1028
    #print "Entering pobyso_get_prec_of_constant_so_sa..."
1029
    prec = c_int(0)
1030
    retc = sollya_lib_get_prec_of_constant(byref(prec), ctExpSo, None)
1031
    if retc == 0:
1032
        #print "pobyso_get_prec_of_constant_so_sa failed."
1033
        return None
1034
    #print "...exiting pobyso_get_prec_of_constant_so_sa."
1035
    return int(prec.value)
1036

    
1037
def pobyso_get_prec_of_range_so_sa(rangeSo):
1038
    """
1039
    Returns the number of bits elements of a range are coded with.
1040
    """
1041
    prec = c_int(0)
1042
    retc = sollya_lib_get_prec_of_range(byref(prec), rangeSo, None)
1043
    if retc == 0:
1044
        return(None)
1045
    return int(prec.value)
1046
# End pobyso_get_prec_of_range_so_sa()
1047

    
1048
def pobyso_get_sage_exp_from_sollya_exp(sollyaExpSo, realField = RR):
1049
    """ Legacy function. See pobyso_get_sage_exp_from_sollya_exp_so_sa. """
1050
    return pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, 
1051
                                                     realField = RR)
1052

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

    
1101

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

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

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

    
1309
def pobyso_interval_to_range_sa_so(intervalSa, precisionSa=None):
1310
    if precisionSa is None:
1311
        precisionSa = intervalSa.parent().precision()
1312
    intervalSo = pobyso_bounds_to_range_sa_so(intervalSa.lower(),\
1313
                                              intervalSa.upper(),\
1314
                                              precisionSa)
1315
    return intervalSo
1316
# End pobyso_interval_to_range_sa_so
1317

    
1318
def pobyso_is_error_so_sa(objSo):
1319
    """
1320
    Thin wrapper around the sollya_lib_obj_is_error() function.
1321
    """
1322
    if sollya_lib_obj_is_error(objSo) != 0:
1323
        return True
1324
    else:
1325
        return False
1326
# End pobyso_is_error-so_sa
1327

    
1328
def pobyso_is_floating_point_number_sa_sa(numberSa):
1329
    """
1330
    Check whether a Sage number is floating point.
1331
    Exception stuff added because numbers other than
1332
    floating-point ones do not have the is_real() attribute.
1333
    """
1334
    try:
1335
        return numberSa.is_real()
1336
    except AttributeError:
1337
        return False
1338
# End pobyso_is_floating_piont_number_sa_sa
1339

    
1340
def pobyso_is_range_so_sa(rangeCandidateSo):
1341
    """
1342
    Thin wrapper over sollya_lib_is_range.
1343
    """
1344
    return sollya_lib_obj_is_range(rangeCandidateSo) != 0
1345
    
1346
# End pobyso_is_range_so_sa
1347

    
1348

    
1349
def pobyso_lib_init():
1350
    sollya_lib_init(None)
1351

    
1352
def pobyso_lib_close():
1353
    sollya_lib_close(None)
1354

    
1355
def pobyso_name_free_variable(freeVariableNameSa):
1356
    """ Legacy function. See pobyso_name_free_variable_sa_so. """
1357
    pobyso_name_free_variable_sa_so(freeVariableNameSa)
1358

    
1359
def pobyso_name_free_variable_sa_so(freeVariableNameSa):
1360
    """
1361
    Set the free variable name in Sollya from a Sage string.
1362
    """
1363
    sollya_lib_name_free_variable(freeVariableNameSa)
1364

    
1365
def pobyso_obj_is_function_so_sa(objSo):
1366
    """
1367
    Check if an object is a function.
1368
    """
1369
    if sollya_lib_obj_is_function(objSo) != 0:
1370
        return True
1371
    else:
1372
        return False
1373
# End pobyso_obj_is_function_so_sa  
1374
    
1375
def pobyso_obj_is_range_so_sa(objSo):
1376
    """
1377
    Check if an object is a function.
1378
    """
1379
    if sollya_lib_obj_is_range(objSo) != 0:
1380
        return True
1381
    else:
1382
        return False
1383
# End pobyso_obj_is_range_so_sa  
1384
    
1385
def pobyso_obj_is_string_so_sa(objSo):
1386
    """
1387
    Check if an object is a function.
1388
    """
1389
    if sollya_lib_obj_is_string(objSo) != 0:
1390
        return True
1391
    else:
1392
        return False
1393
# End pobyso_obj_is_string_so_sa  
1394
    
1395
def pobyso_parse_string(string):
1396
    """ Legacy function. See pobyso_parse_string_sa_so. """
1397
    return pobyso_parse_string_sa_so(string)
1398
 
1399
def pobyso_parse_string_sa_so(string):
1400
    """
1401
    Get the Sollya expression computed from a Sage string or
1402
    a Sollya error object if parsing failed.
1403
    """
1404
    return sollya_lib_parse_string(string)
1405

    
1406
def pobyso_precision_so_sa(ctExpSo):
1407
    """
1408
    Computes the necessary precision to represent a number.
1409
    If x is not zero, it can be uniquely written as x = m · 2e 
1410
    where m is an odd integer and e is an integer. 
1411
    precision(x) returns the number of bits necessary to write m 
1412
    in binary (i.e. ceil(log2(m))).
1413
    """
1414
    #TODO: take care of the special case: 0, @NaN@, @Inf@
1415
    precisionSo = sollya_lib_precision(ctExpSo)
1416
    precisionSa = pobyso_constant_from_int_so_sa(precisionSo)
1417
    sollya_lib_clear_obj(precisionSo)
1418
    return precisionSa
1419
# End pobyso_precision_so_sa
1420

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

    
1587
def pobyso_polynomial_degree_so_sa(polySo):
1588
    """
1589
    Return the degree of a Sollya polynomial as a Sage int.
1590
    """
1591
    degreeSo = sollya_lib_degree(polySo)
1592
    return pobyso_constant_from_int_so_sa(degreeSo)
1593
# End pobyso_polynomial_degree_so_sa
1594

    
1595
def pobyso_polynomial_degree_so_so(polySo):
1596
    """
1597
    Thin wrapper around lib_sollya_degree().
1598
    """
1599
    return sollya_lib_degree(polySo)
1600
# End pobyso_polynomial_degree_so_so
1601
                                                                  
1602
def pobyso_range(rnLowerBound, rnUpperBound):
1603
    """ Legacy function. See pobyso_range_sa_so. """
1604
    return pobyso_range_sa_so(rnLowerBound, rnUpperBound) 
1605

    
1606
def pobyso_range_from_bounds_sa_so(rnLowerBound, rnUpperBound, precSa = None):
1607
    """
1608
    Create a Sollya range from 2 Sage real numbers as bounds
1609
    """
1610
    # TODO check precision stuff.
1611
    sollyaPrecChanged = False
1612
    (initialSollyaPrecSo, initialSollyaPrecSa) = \
1613
        pobyso_get_prec_so_so_sa()
1614
    if precSa is None:
1615
        precSa = max(rnLowerBound.parent().prec(), rnUpperBound.parent().prec())
1616
    if precSa > initialSollyaPrecSa:
1617
        if precSa <= 2:
1618
            print inspect.stack()[0][3], ": precision change <= 2 requested."
1619
        pobyso_set_prec_sa_so(precSa)
1620
        sollyaPrecChanged = True
1621
    rangeSo = sollya_lib_range_from_bounds(get_rn_value(rnLowerBound),
1622
                                           get_rn_value(rnUpperBound))
1623
    if sollyaPrecChanged:
1624
        pobyso_set_prec_so_so(initialSollyaPrecSo)
1625
        sollya_lib_clear_obj(initialSollyaPrecSo)
1626
    return rangeSo
1627
# End pobyso_range_from_bounds_sa_so
1628

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

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

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

    
1851
def pobyso_remez_canonical_so_so(funcSo, \
1852
                                 degreeSo, \
1853
                                 rangeSo, \
1854
                                 weightSo = pobyso_constant_1_sa_so(),\
1855
                                 qualitySo = None):
1856
    """
1857
    All arguments are pointers to Sollya objects.
1858
    The return value is a pointer to a Sollya function.
1859
    """
1860
    if not sollya_lib_obj_is_function(funcSo):
1861
        return(None)
1862
    return(sollya_lib_remez(funcSo, degreeSo, rangeSo, weightSo, qualitySo, None))
1863
# End pobyso_remez_canonical_so_so.
1864

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

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

    
2076
def pobyso_set_canonical_off():
2077
    sollya_lib_set_canonical(sollya_lib_off())
2078

    
2079
def pobyso_set_canonical_on():
2080
    sollya_lib_set_canonical(sollya_lib_on())
2081

    
2082
def pobyso_set_prec(p):
2083
    """ Legacy function. See pobyso_set_prec_sa_so. """
2084
    pobyso_set_prec_sa_so(p)
2085

    
2086
def pobyso_set_prec_sa_so(p):
2087
    #a = c_int(p)
2088
    #precSo = c_void_p(sollya_lib_constant_from_int(a))
2089
    #precSo = sollya_lib_constant_from_int(a)
2090
    precSo = pobyso_constant_from_int_sa_so(p)
2091
    sollya_lib_set_prec(precSo)
2092
    sollya_lib_clear_obj(precSo)
2093
# End pobyso_set_prec_sa_so.
2094

    
2095
def pobyso_set_prec_so_so(newPrecSo):
2096
    sollya_lib_set_prec(newPrecSo)
2097
# End pobyso_set_prec_so_so.
2098
#   
2099
def pobyso_supnorm_sa_sa(poly):
2100
    """
2101
    Computes the supremum norm from Sage input arguments and returns a
2102
    Sage floating-point number whose precision is set by the realFieldSa
2103
    argument.
2104
    TODO: complete this stub!
2105
    """
2106
    print("This function does nothing!")
2107
    return None
2108
# End pobyso_supnorm_sa_sa
2109

    
2110
def pobyso_supnorm_so_sa(polySo, funcSo, intervalSo, errorTypeSo = None,\
2111
                         accuracySo = None, realFieldSa = None):
2112
    """
2113
    Computes the supremum norm from Sollya input arguments and returns a
2114
    Sage floating-point number whose precision is set by the only Sage argument.
2115
    
2116
    The returned value is the maximum of the absolute values of the range
2117
    elements  returned by the Sollya supnorm functions
2118
    """
2119
    supNormRangeSo = pobyso_supnorm_so_so(polySo, 
2120
                                          funcSo, 
2121
                                          intervalSo, 
2122
                                          errorTypeSo,
2123
                                          accuracySo)
2124
    supNormSo = pobyso_range_max_abs_so_so(supNormRangeSo)
2125
    sollya_lib_clear_obj(supNormRangeSo)
2126
    #pobyso_autoprint(supNormSo)
2127
    supNormSa = pobyso_get_constant_as_rn_with_rf_so_sa(supNormSo, realFieldSa)
2128
    sollya_lib_clear_obj(supNormSo)
2129
    return supNormSa 
2130
# End pobyso_supnorm_so_sa.
2131
#
2132
def pobyso_supnorm_so_so(polySo, funcSo, intervalSo, errorTypeSo = None,\
2133
                         accuracySo = None):
2134
    """
2135
    Computes the supnorm of the approximation error between the given 
2136
    polynomial and function. Attention: returns a range!
2137
    errorTypeSo defaults to "absolute".
2138
    accuracySo defaults to 2^(-40).
2139
    """
2140
    if errorTypeSo is None:
2141
        errorTypeSo = sollya_lib_absolute(None)
2142
        errorTypeIsNone = True
2143
    else:
2144
        errorTypeIsNone = False
2145
    #
2146
    if accuracySo is None:
2147
        # Notice the **: we are in Pythonland here!
2148
        accuracySo = pobyso_constant_sa_so(RR(2**(-40)))
2149
        accuracyIsNone = True
2150
    else:
2151
        accuracyIsNone = False
2152
    #pobyso_autoprint(accuracySo)
2153
    resultSo = \
2154
        sollya_lib_supnorm(polySo, funcSo, intervalSo, errorTypeSo, \
2155
                              accuracySo)
2156
    if errorTypeIsNone:
2157
        sollya_lib_clear_obj(errorTypeSo)
2158
    if accuracyIsNone:
2159
        sollya_lib_clear_obj(accuracySo)
2160
    return resultSo
2161
# End pobyso_supnorm_so_so
2162
#
2163
def pobyso_taylor_expansion_no_change_var_so_so(functionSo, 
2164
                                                degreeSo, 
2165
                                                rangeSo,
2166
                                                errorTypeSo=None, 
2167
                                                sollyaPrecSo=None):
2168
    """
2169
    Compute the Taylor expansion without the variable change
2170
    x -> x-intervalCenter.
2171
    If errorTypeSo is None, absolute is used.
2172
    If sollyaPrecSo is None, Sollya internal precision is not changed. 
2173
    """
2174
    # Change internal Sollya precision, if needed.
2175
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
2176
    sollyaPrecChanged = False
2177
    if sollyaPrecSo is None:
2178
        pass
2179
    else:
2180
        sollya_lib_set_prec(sollyaPrecSo)
2181
        sollyaPrecChanged = True
2182
    # Error type stuff: default to absolute.
2183
    if errorTypeSo is None:
2184
        errorTypeIsNone = True
2185
        errorTypeSo = sollya_lib_absolute(None)
2186
    else:
2187
        errorTypeIsNone = False
2188
    intervalCenterSo = sollya_lib_mid(rangeSo, None)
2189
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo,
2190
                                         intervalCenterSo,
2191
                                         rangeSo, errorTypeSo, None)
2192
    # Object taylorFormListSaSo is a Python list of Sollya objects references 
2193
    #  that are copies of the elements of taylorFormSo.
2194
    # pobyso_get_list_elements_so_so clears taylorFormSo.
2195
    (taylorFormListSaSo, numElementsSa, isEndEllipticSa) = \
2196
        pobyso_get_list_elements_so_so(taylorFormSo)
2197
    ## Copy needed here since polySo will be returned and taylorFormListSaSo
2198
    #  will be cleared.
2199
    polySo = sollya_lib_copy_obj(taylorFormListSaSo[0])
2200
    #print "Num elements:", numElementsSa
2201
    sollya_lib_clear_obj(taylorFormSo)
2202
    # No copy_obj needed here: a new objects are created.
2203
    maxErrorSo    = sollya_lib_sup(taylorFormListSaSo[2])
2204
    minErrorSo    = sollya_lib_inf(taylorFormListSaSo[2])
2205
    # List taylorFormListSaSo is not needed anymore.
2206
    pobyso_clear_list_elements_sa_so(taylorFormListSaSo)
2207
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
2208
    absMinErrorSo = sollya_lib_abs(minErrorSo)
2209
    sollya_lib_clear_obj(maxErrorSo)
2210
    sollya_lib_clear_obj(minErrorSo)
2211
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
2212
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
2213
    #
2214
    if errorTypeIsNone:
2215
        sollya_lib_clear_obj(errorTypeSo)
2216
    ## If changed, reset the Sollya working precision.
2217
    if sollyaPrecChanged:
2218
        sollya_lib_set_prec(initialSollyaPrecSo)
2219
    sollya_lib_clear_obj(initialSollyaPrecSo)
2220
    ## According to what error is the largest, return the errors.
2221
    if absMaxErrorSa > absMinErrorSa:
2222
        sollya_lib_clear_obj(absMinErrorSo)
2223
        return (polySo, intervalCenterSo, absMaxErrorSo)
2224
    else:
2225
        sollya_lib_clear_obj(absMaxErrorSo)
2226
        return (polySo, intervalCenterSo, absMinErrorSo)
2227
# end pobyso_taylor_expansion_no_change_var_so_so
2228

    
2229
def pobyso_taylor_expansion_with_change_var_so_so(functionSo, degreeSo, \
2230
                                                  rangeSo, \
2231
                                                  errorTypeSo=None, \
2232
                                                  sollyaPrecSo=None):
2233
    """
2234
    Compute the Taylor expansion with the variable change
2235
    x -> (x-intervalCenter) included.
2236
    """
2237
    # Change Sollya internal precision, if need.
2238
    sollyaPrecChanged = False
2239
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
2240
    if sollyaPrecSo is None:
2241
        pass
2242
    else:
2243
        sollya_lib_set_prec(sollyaPrecSo)
2244
        sollyaPrecChanged = True
2245
    #
2246
    # Error type stuff: default to absolute.
2247
    if errorTypeSo is None:
2248
        errorTypeIsNone = True
2249
        errorTypeSo = sollya_lib_absolute(None)
2250
    else:
2251
        errorTypeIsNone = False
2252
    intervalCenterSo = sollya_lib_mid(rangeSo)
2253
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo, \
2254
                                         intervalCenterSo, \
2255
                                         rangeSo, errorTypeSo, None)
2256
    # Object taylorFormListSaSo is a Python list of Sollya objects references 
2257
    # that are copies of the elements of taylorFormSo.
2258
    # pobyso_get_list_elements_so_so clears taylorFormSo.
2259
    (taylorFormListSaSo, numElements, isEndElliptic) = \
2260
        pobyso_get_list_elements_so_so(taylorFormSo)
2261
    sollya_lib_clear_obj(taylorFormSo)
2262
    polySo = taylorFormListSaSo[0]
2263
    ## Maximum error computation with taylorFormListSaSo[2], a range
2264
    #  holding the actual error. Bounds can be negative.
2265
    maxErrorSo    = sollya_lib_sup(taylorFormListSaSo[2])
2266
    minErrorSo    = sollya_lib_inf(taylorFormListSaSo[2])
2267
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
2268
    absMinErrorSo = sollya_lib_abs(minErrorSo)
2269
    sollya_lib_clear_obj(maxErrorSo)
2270
    sollya_lib_clear_obj(minErrorSo)
2271
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
2272
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
2273
    changeVarExpSo = sollya_lib_build_function_sub(\
2274
                       sollya_lib_build_function_free_variable(),\
2275
                       sollya_lib_copy_obj(intervalCenterSo))
2276
    polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpSo)
2277
    # List taylorFormListSaSo is not needed anymore.
2278
    pobyso_clear_list_elements_sa_so(taylorFormListSaSo)
2279
    sollya_lib_clear_obj(changeVarExpSo)
2280
    # If changed, reset the Sollya working precision.
2281
    if sollyaPrecChanged:
2282
        sollya_lib_set_prec(initialSollyaPrecSo)
2283
    sollya_lib_clear_obj(initialSollyaPrecSo)
2284
    if errorTypeIsNone:
2285
        sollya_lib_clear_obj(errorTypeSo)
2286
    # Do not clear maxErrorSo.
2287
    if absMaxErrorSa > absMinErrorSa:
2288
        sollya_lib_clear_obj(absMinErrorSo)
2289
        return((polyVarChangedSo, intervalCenterSo, absMaxErrorSo))
2290
    else:
2291
        sollya_lib_clear_obj(absMaxErrorSo)
2292
        return((polyVarChangedSo, intervalCenterSo, absMinErrorSo))
2293
# end pobyso_taylor_expansion_with_change_var_so_so
2294

    
2295
def pobyso_taylor(function, degree, point):
2296
    """ Legacy function. See pobysoTaylor_so_so. """
2297
    return(pobyso_taylor_so_so(function, degree, point))
2298

    
2299
def pobyso_taylor_so_so(functionSo, degreeSo, pointSo):
2300
    return(sollya_lib_taylor(functionSo, degreeSo, pointSo))
2301
    
2302
def pobyso_taylorform(function, degree, point = None, 
2303
                      interval = None, errorType=None):
2304
    """ Legacy function. See pobyso_taylorform_sa_sa;"""
2305
    
2306
def pobyso_taylorform_sa_sa(functionSa, \
2307
                            degreeSa, \
2308
                            pointSa, \
2309
                            intervalSa=None, \
2310
                            errorTypeSa=None, \
2311
                            precisionSa=None):
2312
    """
2313
    Compute the Taylor form of 'degreeSa' for 'functionSa' at 'pointSa' 
2314
    for 'intervalSa' with 'errorTypeSa' (a string) using 'precisionSa'. 
2315
    point: must be a Real or a Real interval.
2316
    return the Taylor form as an array
2317
    TODO: take care of the interval and of the point when it is an interval;
2318
          when errorType is not None;
2319
          take care of the other elements of the Taylor form (coefficients 
2320
          errors and delta.
2321
    """
2322
    # Absolute as the default error.
2323
    if errorTypeSa is None:
2324
        errorTypeSo = sollya_lib_absolute()
2325
    elif errorTypeSa == "relative":
2326
        errorTypeSo = sollya_lib_relative()
2327
    elif errortypeSa == "absolute":
2328
        errorTypeSo = sollya_lib_absolute()
2329
    else:
2330
        # No clean up needed.
2331
        return None
2332
    # Global precision stuff
2333
    sollyaPrecisionChangedSa = False
2334
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
2335
    if precisionSa is None:
2336
        precSa = initialSollyaPrecSa
2337
    else:
2338
        if precSa > initialSollyaPrecSa:
2339
            if precSa <= 2:
2340
                print inspect.stack()[0][3], ":precision change <= 2 requested."
2341
            pobyso_set_prec_sa_so(precSa)
2342
            sollyaPrecisionChangedSa = True
2343
    #        
2344
    if len(functionSa.variables()) > 0:
2345
        varSa = functionSa.variables()[0]
2346
        pobyso_name_free_variable_sa_so(str(varSa))
2347
    # In any case (point or interval) the parent of pointSa has a precision
2348
    # method.
2349
    pointPrecSa = pointSa.parent().precision()
2350
    if precSa > pointPrecSa:
2351
        pointPrecSa = precSa
2352
    # In any case (point or interval) pointSa has a base_ring() method.
2353
    pointBaseRingString = str(pointSa.base_ring())
2354
    if re.search('Interval', pointBaseRingString) is None: # Point
2355
        pointSo = pobyso_constant_sa_so(pointSa, pointPrecSa)
2356
    else: # Interval.
2357
        pointSo = pobyso_interval_to_range_sa_so(pointSa, pointPrecSa)
2358
    # Sollyafy the function.
2359
    functionSo = pobyso_parse_string_sa_so(functionSa._assume_str().replace('_SAGE_VAR_', ''))
2360
    if sollya_lib_obj_is_error(functionSo):
2361
        print "pobyso_tailorform: function string can't be parsed!"
2362
        return None
2363
    # Sollyafy the degree
2364
    degreeSo = sollya_lib_constant_from_int(int(degreeSa))
2365
    # Sollyafy the point
2366
    # Call Sollya
2367
    taylorFormSo = \
2368
        sollya_lib_taylorform(functionSo, degreeSo, pointSo, errorTypeSo,\
2369
                                         None)
2370
    sollya_lib_clear_obj(functionSo)
2371
    sollya_lib_clear_obj(degreeSo)
2372
    sollya_lib_clear_obj(pointSo)
2373
    sollya_lib_clear_obj(errorTypeSo)
2374
    (tfsAsList, numElements, isEndElliptic) = \
2375
            pobyso_get_list_elements_so_so(taylorFormSo)
2376
    polySo = tfsAsList[0]
2377
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
2378
    polyRealField = RealField(maxPrecision)
2379
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, polyRealField)
2380
    if sollyaPrecisionChangedSa:
2381
        sollya_lib_set_prec(initialSollyaPrecSo)
2382
    sollya_lib_clear_obj(initialSollyaPrecSo)
2383
    polynomialRing = polyRealField[str(varSa)]
2384
    polySa = polynomial(expSa, polynomialRing)
2385
    taylorFormSa = [polySa]
2386
    # Final clean-up.
2387
    sollya_lib_clear_obj(taylorFormSo)
2388
    return(taylorFormSa)
2389
# End pobyso_taylor_form_sa_sa
2390

    
2391
def pobyso_taylorform_so_so(functionSo, degreeSo, pointSo, intervalSo=None, \
2392
                            errorTypeSo=None):
2393
    createdErrorType = False
2394
    if errorTypeSo is None:
2395
        errorTypeSo = sollya_lib_absolute()
2396
        createdErrorType = True
2397
    else:
2398
        #TODO: deal with the other case.
2399
        pass
2400
    if intervalSo is None:
2401
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
2402
                                         errorTypeSo, None)
2403
    else:
2404
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
2405
                                         intervalSo, errorTypeSo, None)
2406
    if createdErrorType:
2407
        sollya_lib_clear_obj(errorTypeSo)
2408
    return resultSo
2409
        
2410

    
2411
def pobyso_univar_polynomial_print_reverse(polySa):
2412
    """ Legacy function. See pobyso_univar_polynomial_print_reverse_sa_sa. """
2413
    return(pobyso_univar_polynomial_print_reverse_sa_sa(polySa))
2414

    
2415
def pobyso_univar_polynomial_print_reverse_sa_sa(polySa):
2416
    """
2417
    Return the string representation of a univariate polynomial with
2418
    monomials ordered in the x^0..x^n order of the monomials.
2419
    Remember: Sage
2420
    """
2421
    polynomialRing = polySa.base_ring()
2422
    # A very expensive solution:
2423
    # -create a fake multivariate polynomial field with only one variable,
2424
    #   specifying a negative lexicographical order;
2425
    mpolynomialRing = PolynomialRing(polynomialRing.base(), \
2426
                                     polynomialRing.variable_name(), \
2427
                                     1, order='neglex')
2428
    # - convert the univariate argument polynomial into a multivariate
2429
    #   version;
2430
    p = mpolynomialRing(polySa)
2431
    # - return the string representation of the converted form.
2432
    # There is no simple str() method defined for p's class.
2433
    return(p.__str__())
2434
#
2435
#print pobyso_get_prec()  
2436
pobyso_set_prec(165)
2437
#print pobyso_get_prec()  
2438
#a=100
2439
#print type(a)
2440
#id(a)
2441
#print "Max arity: ", pobyso_max_arity
2442
#print "Function tripleDouble (43) as a string: ", pobyso_function_type_as_string(43)
2443
#print "Function None (44) as a string: ", pobyso_function_type_as_string(44)
2444
sys.stderr.write("\t...Pobyso check done.\n")