Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 268

Historique | Voir | Annoter | Télécharger (88,13 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_range_sa_so(rnLowerBoundSa, rnUpperBoundSa, \
104
                                 precisionSa=None):
105
    """
106
    Return a Sollya range from to 2 RealField Sage elements.
107
    The Sollya range element has a sufficient precision to hold all
108
    the digits of the widest of the Sage bounds.
109
    """
110
    # Sanity check.
111
    if rnLowerBoundSa > rnUpperBoundSa:
112
        return None
113
    # Precision stuff.
114
    if precisionSa is None:
115
        # Check for the largest precision.
116
        lbPrecSa = rnLowerBoundSa.parent().precision()
117
        ubPrecSa = rnLowerBoundSa.parent().precision()
118
        maxPrecSa = max(lbPrecSa, ubPrecSa)
119
    else:
120
        maxPrecSa = precisionSa
121
    # From Sage to Sollya bounds.
122
#    lowerBoundSo = sollya_lib_constant(get_rn_value(rnLowerBoundSa),
123
#                                       maxPrecSa)
124
    lowerBoundSo = pobyso_constant_sa_so(rnLowerBoundSa,
125
                                         maxPrecSa)
126
    upperBoundSo = pobyso_constant_sa_so(rnUpperBoundSa,
127
                                       maxPrecSa)
128
    
129
    # From Sollya bounds to range.
130
    rangeSo = sollya_lib_range(lowerBoundSo, upperBoundSo)
131
    # Back to original precision.
132
    # Clean up
133
    sollya_lib_clear_obj(lowerBoundSo)
134
    sollya_lib_clear_obj(upperBoundSo)
135
    return rangeSo
136
# End pobyso_bounds_to_range_sa_so
137

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

    
160
def pobyso_build_function_sub_so_so(exp1So, exp2So):
161
    return sollya_lib_build_function_sub(exp1So, exp2So)
162

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

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

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

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

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

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

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

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

    
319

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

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

    
355
def pobyso_constant_1_sa_so():
356
    """
357
    Obvious.
358
    """
359
    return(pobyso_constant_from_int_sa_so(1))
360

    
361
def pobyso_constant_from_int(anInt):
362
    """ Legacy function. See pobyso_constant_from_int_sa_so. """
363
    return pobyso_constant_from_int_sa_so(anInt)
364

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

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

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

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

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

    
411
def pobyso_error_so():
412
    return sollya_lib_error(None)
413
# End pobyso_error().
414

    
415
def pobyso_evaluate_so_sa(funcSo, argumentSo):
416
    """
417
    Evaluates funcSo for arguemntSo through sollya_lib_evaluate() and return
418
    the result as a Sage object
419
    """
420
    evalSo = sollya_lib_evaluate(funcSo, argumentSo)
421
    if pobyso_is_error_so_sa(evalSo):
422
        return None
423
    if pobyso_is_range_so_sa(evalSo):
424
        retVal = pobyso_range_to_interval_so_sa(evalSo)
425
    else:
426
        retVal = pobyso_get_constant_as_rn(evalSo)
427
    sollya_lib_clear_obj(evalSo)
428
    return retVal
429
# End pobyso_evaluate_so_sa.
430

    
431
def pobyso_evaluate_so_so(funcSo, argumentSo):
432
    """
433
    Evaluates funcSo for arguemntSo through sollya_lib_evaluate().
434
    """
435
    return sollya_lib_evaluate(funcSo, argumentSo)
436
# End pobyso_evaluate_so_so.
437
#
438
def pobyso_diff_so_so(funcSo):
439
    """
440
    Very thin wrapper around sollya_lib_diff.
441
    """
442
    ## TODO: add a check to make sure funcSo is a functional expression.
443
    return sollya_lib_diff(funcSo)
444

    
445
def pobyso_dirty_find_zeros_so_so(funcSo, rangeSo):
446
    """
447
    Thin wrapper over sollya_lib_dirtyfindzeros()
448
    """
449
    return sollya_lib_dirtyfindzeros(funcSo, rangeSo)
450
# End pobys_dirty_find_zeros
451

    
452
def pobyso_dirty_inf_norm_so_so(funcSo, rangeSo, preSo=None):
453
    """
454
    Thin wrapper around sollya_dirtyinfnorm().
455
    """
456
    # TODO: manage the precision change.
457
    
458
    return sollya_lib_dirtyinfnorm(funcSo, rangeSo)
459
# End pobyso_dirty_inf_norm_so_so
460

    
461
def pobyso_find_zeros_so_so(funcSo, rangeSo):
462
    """
463
    Thin wrapper over sollya_lib_findzeros()
464
    """
465
    return sollya_lib_findzeros(funcSo, rangeSo)
466
# End pobys_find_zeros
467

    
468
def pobyso_float_list_so_sa(listSo):
469
    """
470
    Return a Sollya list of floating-point numbers as a Sage list of 
471
    floating-point numbers.
472
    TODO: add a test to make sure that each element of the list is a constant.
473
    """
474
    listSa   = []
475
    ## The function returns none if the list is empty or an error has happened.
476
    retVal = pobyso_get_list_elements_so_so(listSo)
477
    if retVal is None:
478
        return listSa
479
    ## Just in case the interface is changed and an empty list is returned 
480
    #  instead of None.
481
    elif len(retVal) == 0:
482
        return listSa
483
    else:
484
        ## Remember pobyso_get_list_elements_so_so returns more information
485
        #  than just the elements of the list (# elements, is_elliptic)
486
        listSaSo, numElements, isEndElliptic = retVal
487
    ## Return an empty list.
488
    if numElements == 0:
489
        return listSa
490
    ## Search first for the maximum precision of the elements
491
    maxPrecSa = 0
492
    for floatSo in listSaSo:
493
        #pobyso_autoprint(floatSo)
494
        curPrecSa =  pobyso_get_prec_of_constant_so_sa(floatSo)
495
        if curPrecSa > maxPrecSa:
496
            maxPrecSa = curPrecSa
497
    ##
498
    RF = RealField(maxPrecSa)
499
    ##
500
    for floatSo in listSaSo:
501
        listSa.append(pobyso_get_constant_as_rn_with_rf_so_sa(floatSo))
502
    return listSa
503
# End pobyso_float_list_so_sa
504

    
505
def pobyso_float_poly_sa_so(polySa, precSa = None):
506
    """
507
    Create a Sollya polynomial from a Sage RealField polynomial.
508
    """
509
    ## TODO: filter arguments.
510
    ## Precision. If a precision is given, convert the polynomial
511
    #  into the right polynomial field. If not convert it straight
512
    #  to Sollya.
513
    sollyaPrecChanged = False 
514
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
515
    if precSa is None:
516
        precSa = polySa.parent().base_ring().precision()
517
    if (precSa > initialSollyaPrecSa):
518
        assert precSa >= 2, "Precision change <2 requested"
519
        if precSa <= 2:
520
            print inspect.stack()[0][3], ": precision change <= 2 requested"
521
        precSo = pobyso_constant_from_int(precSa)
522
        pobyso_set_prec_so_so(precSo)
523
        sollya_lib_clear_obj(precSo)
524
        sollyaPrecChanged = True    
525
    ## Get exponents and coefficients.
526
    exponentsSa     = polySa.exponents()
527
    coefficientsSa  = polySa.coefficients()
528
    ## Build the polynomial.
529
    polySo = None
530
    for coefficientSa, exponentSa in zip(coefficientsSa, exponentsSa):
531
        #print coefficientSa.n(prec=precSa), exponentSa
532
        coefficientSo = \
533
            pobyso_constant_sa_so(coefficientSa)
534
        #pobyso_autoprint(coefficientSo)
535
        exponentSo = \
536
            pobyso_constant_from_int_sa_so(exponentSa)
537
        #pobyso_autoprint(exponentSo)
538
        monomialSo = sollya_lib_build_function_pow(
539
                       sollya_lib_build_function_free_variable(),
540
                       exponentSo)
541
        polyTermSo = sollya_lib_build_function_mul(coefficientSo,
542
                                                       monomialSo)
543
        if polySo is None:
544
            polySo = polyTermSo
545
        else:
546
            polySo = sollya_lib_build_function_add(polySo, polyTermSo)
547
    if sollyaPrecChanged:
548
        pobyso_set_prec_so_so(initialSollyaPrecSo)
549
        sollya_lib_clear_obj(initialSollyaPrecSo)
550
    return polySo
551
# End pobyso_float_poly_sa_so    
552

    
553
def pobyso_float_poly_so_sa(polySo, realFieldSa=None):
554
    """
555
    Convert a Sollya polynomial into a Sage floating-point polynomial.
556
    If no realField is given, a RealField corresponding to the maximum 
557
    precision of the coefficients is internally computed.
558
    The real field is not returned but can be easily retrieved from 
559
    the polynomial itself.
560
    ALGORITHM:
561
    - (optional) compute the RealField of the coefficients;
562
    - convert the Sollya expression into a Sage expression;
563
    - convert the Sage expression into a Sage polynomial
564
    """    
565
    if realFieldSa is None:
566
        expressionPrecSa = pobyso_get_max_prec_of_exp_so_sa(polySo)
567
        #print "Maximum precision of Sollya polynomial coefficients:", expressionPrecSa
568
        if expressionPrecSa < 2 or expressionPrecSa > 2147483391:
569
            print "Maximum degree of expression:", expressionPrecSa
570
        realFieldSa      = RealField(expressionPrecSa)
571
    #print "Sollya expression before...",
572
    #pobyso_autoprint(polySo)
573

    
574
    expressionSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo,
575
                                                             realFieldSa)
576
    #print "...Sollya expression after."
577
    #pobyso_autoprint(polySo)
578
    polyVariableSa = expressionSa.variables()[0]
579
    polyRingSa     = realFieldSa[str(polyVariableSa)]
580
    #print polyRingSa
581
    # Do not use the polynomial(expressionSa, ring=polyRingSa) form!
582
    polynomialSa = polyRingSa(expressionSa)
583
    polyCoeffsListSa = polynomialSa.coefficients()
584
    #for coeff in polyCoeffsListSa:
585
    #    print coeff.abs().n()
586
    return polynomialSa
587
# End pobyso_float_poly_so_sa
588

    
589
def pobyso_free_variable():
590
    """
591
    Ultra thin wrapper around the sollya_lib_function_build_free_variable function.
592
    """
593
    return sollya_lib_build_function_free_variable()
594
    
595
def pobyso_function_type_as_string(funcType):
596
    """ Legacy function. See pobyso_function_type_as_string_so_sa. """
597
    return(pobyso_function_type_as_string_so_sa(funcType))
598

    
599
def pobyso_function_type_as_string_so_sa(funcType):
600
    """
601
    Numeric Sollya function codes -> Sage mathematical function names.
602
    Notice that pow -> ^ (a la Sage, not a la Python).
603
    """
604
    if funcType == SOLLYA_BASE_FUNC_ABS:
605
        return "abs"
606
    elif funcType == SOLLYA_BASE_FUNC_ACOS:
607
        return "arccos"
608
    elif funcType == SOLLYA_BASE_FUNC_ACOSH:
609
        return "arccosh"
610
    elif funcType == SOLLYA_BASE_FUNC_ADD:
611
        return "+"
612
    elif funcType == SOLLYA_BASE_FUNC_ASIN:
613
        return "arcsin"
614
    elif funcType == SOLLYA_BASE_FUNC_ASINH:
615
        return "arcsinh"
616
    elif funcType == SOLLYA_BASE_FUNC_ATAN:
617
        return "arctan"
618
    elif funcType == SOLLYA_BASE_FUNC_ATANH:
619
        return "arctanh"
620
    elif funcType == SOLLYA_BASE_FUNC_CEIL:
621
        return "ceil"
622
    elif funcType == SOLLYA_BASE_FUNC_CONSTANT:
623
        return "cte"
624
    elif funcType == SOLLYA_BASE_FUNC_COS:
625
        return "cos"
626
    elif funcType == SOLLYA_BASE_FUNC_COSH:
627
        return "cosh"
628
    elif funcType == SOLLYA_BASE_FUNC_DIV:
629
        return "/"
630
    elif funcType == SOLLYA_BASE_FUNC_DOUBLE:
631
        return "double"
632
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEDOUBLE:
633
        return "doubleDouble"
634
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEEXTENDED:
635
        return "doubleDxtended"
636
    elif funcType == SOLLYA_BASE_FUNC_ERF:
637
        return "erf"
638
    elif funcType == SOLLYA_BASE_FUNC_ERFC:
639
        return "erfc"
640
    elif funcType == SOLLYA_BASE_FUNC_EXP:
641
        return "exp"
642
    elif funcType == SOLLYA_BASE_FUNC_EXP_M1:
643
        return "expm1"
644
    elif funcType == SOLLYA_BASE_FUNC_FLOOR:
645
        return "floor"
646
    elif funcType == SOLLYA_BASE_FUNC_FREE_VARIABLE:
647
        return "freeVariable"
648
    elif funcType == SOLLYA_BASE_FUNC_HALFPRECISION:
649
        return "halfPrecision"
650
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYCONSTANT:
651
        return "libraryConstant"
652
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYFUNCTION:
653
        return "libraryFunction"
654
    elif funcType == SOLLYA_BASE_FUNC_LOG:
655
        return "log"
656
    elif funcType == SOLLYA_BASE_FUNC_LOG_10:
657
        return "log10"
658
    elif funcType == SOLLYA_BASE_FUNC_LOG_1P:
659
        return "log1p"
660
    elif funcType == SOLLYA_BASE_FUNC_LOG_2:
661
        return "log2"
662
    elif funcType == SOLLYA_BASE_FUNC_MUL:
663
        return "*"
664
    elif funcType == SOLLYA_BASE_FUNC_NEARESTINT:
665
        return "round"
666
    elif funcType == SOLLYA_BASE_FUNC_NEG:
667
        return "__neg__"
668
    elif funcType == SOLLYA_BASE_FUNC_PI:
669
        return "pi"
670
    elif funcType == SOLLYA_BASE_FUNC_POW:
671
        return "^"
672
    elif funcType == SOLLYA_BASE_FUNC_PROCEDUREFUNCTION:
673
        return "procedureFunction"
674
    elif funcType == SOLLYA_BASE_FUNC_QUAD:
675
        return "quad"
676
    elif funcType == SOLLYA_BASE_FUNC_SIN:
677
        return "sin"
678
    elif funcType == SOLLYA_BASE_FUNC_SINGLE:
679
        return "single"
680
    elif funcType == SOLLYA_BASE_FUNC_SINH:
681
        return "sinh"
682
    elif funcType == SOLLYA_BASE_FUNC_SQRT:
683
        return "sqrt"
684
    elif funcType == SOLLYA_BASE_FUNC_SUB:
685
        return "-"
686
    elif funcType == SOLLYA_BASE_FUNC_TAN:
687
        return "tan"
688
    elif funcType == SOLLYA_BASE_FUNC_TANH:
689
        return "tanh"
690
    elif funcType == SOLLYA_BASE_FUNC_TRIPLEDOUBLE:
691
        return "tripleDouble"
692
    else:
693
        return None
694

    
695
def pobyso_get_constant(rnArgSa, constSo):
696
    """ Legacy function. See pobyso_get_constant_so_sa. """
697
    return pobyso_get_constant_so_sa(rnArgSa, constSo)
698
# End pobyso_get_constant
699

    
700
def pobyso_get_constant_so_sa(rnArgSa, constSo):
701
    """
702
    Set the value of rnArgSo to the value of constSo in MPFR_RNDN mode.
703
    rnArg must already exist and belong to some RealField.
704
    We assume that constSo points to a Sollya constant.
705
    """
706
    outcome = sollya_lib_get_constant(get_rn_value(rnArgSa), constSo)
707
    if outcome == 0: # Failure because constSo is not a constant expression.
708
        return None
709
    else:
710
        return outcome
711
# End  pobyso_get_constant_so_sa
712
   
713
def pobyso_get_constant_as_rn(ctExpSo):
714
    """ 
715
    Legacy function. See pobyso_get_constant_as_rn_so_sa. 
716
    """ 
717
    return(pobyso_get_constant_as_rn_so_sa(ctExpSo))
718
    
719
def pobyso_get_constant_as_rn_so_sa(constExpSo):
720
    """
721
    Get a Sollya constant as a Sage "real number".
722
    The precision of the floating-point number returned is that of the Sollya
723
    constant.
724
    """
725
    #print "Before computing precision of variable..."
726
    #pobyso_autoprint(constExpSo)
727
    precisionSa  = pobyso_get_prec_of_constant_so_sa(constExpSo)
728
    #print "precisionSa:", precisionSa
729
    ## If the expression can not be exactly converted, None is returned.
730
    #  In this case opt for the Sollya current expression.
731
    if precisionSa is None:
732
        precisionSa = pobyso_get_prec_so_sa()
733
    RRRR = RealField(precisionSa)
734
    rnSa = RRRR(0)
735
    outcome = sollya_lib_get_constant(get_rn_value(rnSa), constExpSo)
736
    if outcome == 0:
737
        return None
738
    else:
739
        return rnSa
740
# End pobyso_get_constant_as_rn_so_sa
741

    
742
def pobyso_get_constant_as_rn_with_rf(ctExp, realField):
743
    """ 
744
    Legacy function. See pobyso_get_constant_as_rn_with_rf_so_sa.
745
    """
746
    return pobyso_get_constant_as_rn_with_rf_so_sa(ctExp, realField)
747
# End pobyso_get_constant_as_rn_with_rf
748
    
749
def pobyso_get_constant_as_rn_with_rf_so_sa(ctExpSo, realFieldSa = None):
750
    """
751
    Get a Sollya constant as a Sage "real number".
752
    If no real field is specified, the precision of the floating-point number 
753
    returned is that of the Sollya constant.
754
    Otherwise is is that of the real field. Hence rounding may happen.
755
    """
756
    if realFieldSa is None:
757
        return pobyso_get_constant_as_rn_so_sa(ctExpSo)
758
    rnSa = realFieldSa(0)
759
    outcome = sollya_lib_get_constant(get_rn_value(rnSa), ctExpSo)
760
    if outcome == 0:
761
        return None
762
    else:
763
        return rnSa
764
# End pobyso_get_constant_as_rn_with_rf_so_sa
765

    
766
def pobyso_get_free_variable_name():
767
    """ 
768
    Legacy function. See pobyso_get_free_variable_name_so_sa.
769
    """
770
    return(pobyso_get_free_variable_name_so_sa())
771

    
772
def pobyso_get_free_variable_name_so_sa():
773
    return sollya_lib_get_free_variable_name()
774
    
775
def pobyso_get_function_arity(expressionSo):
776
    """ 
777
    Legacy function. See pobyso_get_function_arity_so_sa.
778
    """
779
    return(pobyso_get_function_arity_so_sa(expressionSo))
780

    
781
def pobyso_get_function_arity_so_sa(expressionSo):
782
    arity = c_int(0)
783
    sollya_lib_get_function_arity(byref(arity),expressionSo)
784
    return int(arity.value)
785

    
786
def pobyso_get_head_function(expressionSo):
787
    """ 
788
    Legacy function. See pobyso_get_head_function_so_sa. 
789
    """
790
    return(pobyso_get_head_function_so_sa(expressionSo)) 
791

    
792
def pobyso_get_head_function_so_sa(expressionSo):
793
    functionType = c_int(0)
794
    sollya_lib_get_head_function(byref(functionType), expressionSo)
795
    return int(functionType.value)
796

    
797
def pobyso_get_interval_from_range_so_sa(soRange, realIntervalFieldSa = None ):
798
    """
799
    Return the Sage interval corresponding to the Sollya range argument.
800
    If no reaIntervalField is passed as an argument, the interval bounds are not
801
    rounded: they are elements of RealIntervalField of the "right" precision
802
    to hold all the digits.
803
    """
804
    prec = c_int(0)
805
    if realIntervalFieldSa is None:
806
        retval = sollya_lib_get_prec_of_range(byref(prec), soRange, None)
807
        if retval == 0:
808
            return None
809
        realIntervalFieldSa = RealIntervalField(prec.value)
810
    intervalSa = realIntervalFieldSa(0,0)
811
    retval = \
812
        sollya_lib_get_interval_from_range(get_interval_value(intervalSa),\
813
                                           soRange)
814
    if retval == 0:
815
        return None
816
    return intervalSa
817
# End pobyso_get_interval_from_range_so_sa
818

    
819
def pobyso_get_list_elements(soObj):
820
    """ Legacy function. See pobyso_get_list_elements_so_so. """
821
    return pobyso_get_list_elements_so_so(soObj)
822
 
823
def pobyso_get_list_elements_so_so(objectListSo):
824
    """
825
    Get the Sollya list elements as a Sage/Python array of Sollya objects.
826
    
827
    INPUT:
828
    - objectListSo: a Sollya list of Sollya objects.
829
    
830
    OUTPUT:
831
    - a Sage/Python tuple made of:
832
      - a Sage/Python list of Sollya objects,
833
      - a Sage/Python int holding the number of elements,
834
      - a Sage/Python int stating (!= 0) that the list is end-elliptic.
835
    NOTE::
836
        We recover the addresses of the Sollya object from the list of pointers
837
        returned by sollya_lib_get_list_elements. The list itself is freed.
838
    TODO::
839
        Figure out what to do with numElements since the number of elements
840
        can easily be recovered from the list itself. 
841
        Ditto for isEndElliptic.
842
    """
843
    listAddress = POINTER(c_longlong)()
844
    numElements = c_int(0)
845
    isEndElliptic = c_int(0)
846
    listAsSageList = []
847
    result = sollya_lib_get_list_elements(byref(listAddress),\
848
                                          byref(numElements),\
849
                                          byref(isEndElliptic),\
850
                                          objectListSo)
851
    if result == 0 :
852
        return None
853
    for i in xrange(0, numElements.value, 1):
854
       #listAsSageList.append(sollya_lib_copy_obj(listAddress[i]))
855
       listAsSageList.append(listAddress[i])
856
       # Clear each of the elements returned by Sollya.
857
       #sollya_lib_clear_obj(listAddress[i])
858
    # Free the list itself.   
859
    sollya_lib_free(listAddress)
860
    return (listAsSageList, numElements.value, isEndElliptic.value)
861

    
862
def pobyso_get_max_prec_of_exp(soExp):
863
    """ Legacy function. See pobyso_get_max_prec_of_exp_so_sa. """
864
    return pobyso_get_max_prec_of_exp_so_sa(soExp)
865

    
866
def pobyso_get_max_prec_of_exp_so_sa(expSo):
867
    """
868
    Get the maximum precision used for the numbers in a Sollya expression.
869
    
870
    Arguments:
871
    soExp -- a Sollya expression pointer
872
    Return value:
873
    A Python integer
874
    TODO: 
875
    - error management;
876
    - correctly deal with numerical type such as DOUBLEEXTENDED.
877
    """
878
    if expSo is None:
879
        print inspect.stack()[0][3], ": expSo is None."
880
        return 0
881
    maxPrecision = 0
882
    minConstPrec = 0
883
    currentConstPrec = 0
884
    #pobyso_autoprint(expSo)
885
    operator = pobyso_get_head_function_so_sa(expSo)
886
    if (operator != SOLLYA_BASE_FUNC_CONSTANT) and \
887
    (operator != SOLLYA_BASE_FUNC_FREE_VARIABLE):
888
        (arity, subexpressions) = pobyso_get_subfunctions_so_sa(expSo)
889
        for i in xrange(arity):
890
            maxPrecisionCandidate = \
891
                pobyso_get_max_prec_of_exp_so_sa(subexpressions[i])
892
            if maxPrecisionCandidate > maxPrecision:
893
                maxPrecision = maxPrecisionCandidate
894
        return maxPrecision
895
    elif operator == SOLLYA_BASE_FUNC_CONSTANT:
896
        #minConstPrec = pobyso_get_min_prec_of_constant_so_sa(expSo)
897
        #currentConstPrec = pobyso_get_min_prec_of_constant_so_sa(soExp)
898
        #print minConstPrec, " - ", currentConstPrec 
899
        return pobyso_get_min_prec_of_constant_so_sa(expSo)
900
    
901
    elif operator == SOLLYA_BASE_FUNC_FREE_VARIABLE:
902
        return 0
903
    else:
904
        print "pobyso_get_max_prec_of_exp_so_sa: unexepected operator."
905
        return 0
906

    
907
def pobyso_get_min_prec_of_constant_so_sa(constExpSo):
908
    """
909
    Get the minimum precision necessary to represent the value of a Sollya
910
    constant.
911
    MPFR_MIN_PREC and powers of 2 are taken into account.
912
    We assume that constExpSo is a pointer to a Sollay constant expression.
913
    """
914
    constExpAsRnSa = pobyso_get_constant_as_rn_so_sa(constExpSo)
915
    return(min_mpfr_size(get_rn_value(constExpAsRnSa)))
916

    
917
def pobyso_get_poly_so_sa(polySo, realFieldSa=None):
918
    """
919
    Convert a Sollya polynomial into a Sage polynomial.
920
    Legacy function. Use pobyso_float_poly_so_sa() instead.
921
    """    
922
    return pobyso_float_poly_so_sa(polySo,realFieldSa)
923
# End pobyso_get_poly_so_sa
924

    
925
def pobyso_get_prec():
926
    """ Legacy function. See pobyso_get_prec_so_sa(). """
927
    return pobyso_get_prec_so_sa()
928

    
929
def pobyso_get_prec_so():
930
    """
931
    Get the current default precision in Sollya.
932
    The return value is a Sollya object.
933
    Usefull when modifying the precision back and forth by avoiding
934
    extra conversions.
935
    """
936
    return sollya_lib_get_prec(None)
937
    
938
def pobyso_get_prec_so_sa():
939
    """
940
    Get the current default precision in Sollya.
941
    The return value is Sage/Python int.
942
    """
943
    precSo = sollya_lib_get_prec()
944
    precSa = pobyso_constant_from_int_so_sa(precSo)
945
    sollya_lib_clear_obj(precSo)
946
    return precSa
947
# End pobyso_get_prec_so_sa.
948

    
949
def pobyso_get_prec_so_so_sa():
950
    """
951
    Return the current precision both as a Sollya object and a
952
    Sage integer as hybrid tuple.
953
    To avoid multiple calls for precision manipulations. 
954
    """
955
    precSo = sollya_lib_get_prec()
956
    precSa = pobyso_constant_from_int_so_sa(precSo)
957
    return (precSo, int(precSa))
958
    
959
def pobyso_get_prec_of_constant(ctExpSo):
960
    """ Legacy function. See pobyso_get_prec_of_constant_so_sa. """
961
    return pobyso_get_prec_of_constant_so_sa(ctExpSo)
962

    
963
def pobyso_get_prec_of_constant_so_sa(ctExpSo):
964
    """
965
    Tries to find a precision to represent ctExpSo without rounding.
966
    If not possible, returns None.
967
    """
968
    #print "Entering pobyso_get_prec_of_constant_so_sa..."
969
    prec = c_int(0)
970
    retc = sollya_lib_get_prec_of_constant(byref(prec), ctExpSo, None)
971
    if retc == 0:
972
        #print "pobyso_get_prec_of_constant_so_sa failed."
973
        return None
974
    #print "...exiting pobyso_get_prec_of_constant_so_sa."
975
    return int(prec.value)
976

    
977
def pobyso_get_prec_of_range_so_sa(rangeSo):
978
    """
979
    Returns the number of bits elements of a range are coded with.
980
    """
981
    prec = c_int(0)
982
    retc = sollya_lib_get_prec_of_range(byref(prec), rangeSo, None)
983
    if retc == 0:
984
        return(None)
985
    return int(prec.value)
986
# End pobyso_get_prec_of_range_so_sa()
987

    
988
def pobyso_get_sage_exp_from_sollya_exp(sollyaExpSo, realField = RR):
989
    """ Legacy function. See pobyso_get_sage_exp_from_sollya_exp_so_sa. """
990
    return pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, 
991
                                                     realField = RR)
992

    
993
def pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, realFieldSa = RR):
994
    """
995
    Get a Sage expression from a Sollya expression. 
996
    Currently only tested with polynomials with floating-point coefficients.
997
    Notice that, in the returned polynomial, the exponents are RealNumbers.
998
    """
999
    #pobyso_autoprint(sollyaExp)
1000
    operatorSa = pobyso_get_head_function_so_sa(sollyaExpSo)
1001
    sollyaLibFreeVariableName = sollya_lib_get_free_variable_name()
1002
    ## Get rid of the "_"'s in "_x_", if any.
1003
    sollyaLibFreeVariableName = re.sub('_', '', sollyaLibFreeVariableName)
1004
    # Constants and the free variable are special cases.
1005
    # All other operator are dealt with in the same way.
1006
    if (operatorSa != SOLLYA_BASE_FUNC_CONSTANT) and \
1007
       (operatorSa != SOLLYA_BASE_FUNC_FREE_VARIABLE):
1008
        (aritySa, subexpressionsSa) = pobyso_get_subfunctions_so_sa(sollyaExpSo)
1009
        if aritySa == 1:
1010
            sageExpSa = eval(pobyso_function_type_as_string_so_sa(operatorSa) + \
1011
            "(" + pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], \
1012
            realFieldSa) + ")")
1013
        elif aritySa == 2:
1014
            # We do not get through the preprocessor.
1015
            # The "^" operator is then a special case.
1016
            if operatorSa == SOLLYA_BASE_FUNC_POW:
1017
                operatorAsStringSa = "**"
1018
            else:
1019
                operatorAsStringSa = \
1020
                    pobyso_function_type_as_string_so_sa(operatorSa)
1021
            sageExpSa = \
1022
              eval("pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], realFieldSa)"\
1023
              + " " + operatorAsStringSa + " " + \
1024
                   "pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[1], realFieldSa)")
1025
        # We do not know yet how to deal with arity >= 3 
1026
        # (is there any in Sollya anyway?).
1027
        else:
1028
            sageExpSa = eval('None')
1029
        return sageExpSa
1030
    elif operatorSa == SOLLYA_BASE_FUNC_CONSTANT:
1031
        #print "This is a constant"
1032
        return pobyso_get_constant_as_rn_with_rf_so_sa(sollyaExpSo, realFieldSa)
1033
    elif operatorSa == SOLLYA_BASE_FUNC_FREE_VARIABLE:
1034
        #print "This is the free variable"
1035
        return eval(sollyaLibFreeVariableName)
1036
    else:
1037
        print "Unexpected"
1038
        return eval('None')
1039
# End pobyso_get_sage_exp_from_sollya_exp_so_sa
1040

    
1041

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

    
1137
def pobyso_guess_degree_so_sa(functionSo, rangeSo, errorSo, weightSo=None, \
1138
                              degreeBoundSo=None):
1139
    """
1140
    Thin wrapper around the guessdegree function.
1141
    Nevertheless, some precision control stuff has been appended.
1142
    """
1143
    # Deal with Sollya internal precision issues: if it is too small, 
1144
    # compared with the error, increases it to about twice -log2(error).
1145
    errorSa = pobyso_get_constant_as_rn_with_rf_so_sa(errorSo)
1146
    log2ErrorSa = errorSa.log2()
1147
    if log2ErrorSa < 0:
1148
        neededPrecisionSa = int(2 * int(-log2ErrorSa) / 64) * 64
1149
    else:
1150
        neededPrecisionSa = int(2 * int(log2ErrorSa) / 64) * 64
1151
    #print "Needed precision:", neededPrecisionSa
1152
    sollyaPrecisionHasChanged = False
1153
    (initialPrecSo, initialPrecSa) = pobyso_get_prec_so_so_sa()
1154
    if neededPrecisionSa > initialPrecSa:
1155
        if neededPrecisionSa <= 2:
1156
            print inspect.stack()[0][3], ": precision change <= 2 requested."
1157
        pobyso_set_prec_sa_so(neededPrecisionSa)
1158
        sollyaPrecisionHasChanged = True
1159
    #print "Guessing degree..."
1160
    # weightSo and degreeBoundsSo are optional arguments.
1161
    # As declared, sollya_lib_guessdegree must take 5 arguments.
1162
    if weightSo is None:
1163
        degreeRangeSo = sollya_lib_guessdegree(functionSo, rangeSo, errorSo, 
1164
                                               0, 0, None)
1165
    elif degreeBoundSo is None:
1166
        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, \
1167
                                                errorSo, weightSo, 0, None)
1168
    else:
1169
        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, errorSo, \
1170
                                                weightSo, degreeBoundSo, None)
1171
    #print "...degree guess done."
1172
    # Restore internal precision, if applicable.
1173
    if sollyaPrecisionHasChanged:
1174
        pobyso_set_prec_so_so(initialPrecSo)
1175
        sollya_lib_clear_obj(initialPrecSo)
1176
    degreeIntervalSa = pobyso_range_to_interval_so_sa(degreeRangeSo)
1177
    sollya_lib_clear_obj(degreeRangeSo)
1178
    # When ok, both bounds match.
1179
    # When the degree bound is too low, the upper bound is the degree
1180
    # for which the error can be honored.
1181
    # When it really goes wrong, the upper bound is infinity. 
1182
    if degreeIntervalSa.lower() == degreeIntervalSa.upper():
1183
        return int(degreeIntervalSa.lower())
1184
    else:
1185
        if degreeIntervalSa.upper().is_infinity():
1186
            return None
1187
        else:
1188
            return int(degreeIntervalSa.upper())
1189
    # End pobyso_guess_degree_so_sa
1190

    
1191
def pobyso_inf_so_so(intervalSo):
1192
    """
1193
    Very thin wrapper around sollya_lib_inf().
1194
    """
1195
    return sollya_lib_inf(intervalSo)
1196
# End pobyso_inf_so_so.
1197
    
1198
def pobyso_infnorm_so_so(func, interval, file = None, intervalList = None):
1199
    print "Do not use this function. User pobyso_supnorm_so_so instead."
1200
    return None
1201

    
1202
def pobyso_interval_to_range_sa_so(intervalSa, precisionSa=None):
1203
    if precisionSa is None:
1204
        precisionSa = intervalSa.parent().precision()
1205
    intervalSo = pobyso_bounds_to_range_sa_so(intervalSa.lower(),\
1206
                                              intervalSa.upper(),\
1207
                                              precisionSa)
1208
    return intervalSo
1209
# End pobyso_interval_to_range_sa_so
1210

    
1211
def pobyso_is_error_so_sa(objSo):
1212
    """
1213
    Thin wrapper around the sollya_lib_obj_is_error() function.
1214
    """
1215
    if sollya_lib_obj_is_error(objSo) != 0:
1216
        return True
1217
    else:
1218
        return False
1219
# End pobyso_is_error-so_sa
1220

    
1221
def pobyso_is_floating_point_number_sa_sa(numberSa):
1222
    """
1223
    Check whether a Sage number is floating point.
1224
    Exception stuff added because numbers other than
1225
    floating-point ones do not have the is_real() attribute.
1226
    """
1227
    try:
1228
        return numberSa.is_real()
1229
    except AttributeError:
1230
        return False
1231
# End pobyso_is_floating_piont_number_sa_sa
1232

    
1233
def pobyso_is_range_so_sa(rangeCandidateSo):
1234
    """
1235
    Thin wrapper over sollya_lib_is_range.
1236
    """
1237
    return sollya_lib_obj_is_range(rangeCandidateSo) != 0
1238
    
1239
# End pobyso_is_range_so_sa
1240

    
1241

    
1242
def pobyso_lib_init():
1243
    sollya_lib_init(None)
1244

    
1245
def pobyso_lib_close():
1246
    sollya_lib_close(None)
1247

    
1248
def pobyso_name_free_variable(freeVariableNameSa):
1249
    """ Legacy function. See pobyso_name_free_variable_sa_so. """
1250
    pobyso_name_free_variable_sa_so(freeVariableNameSa)
1251

    
1252
def pobyso_name_free_variable_sa_so(freeVariableNameSa):
1253
    """
1254
    Set the free variable name in Sollya from a Sage string.
1255
    """
1256
    sollya_lib_name_free_variable(freeVariableNameSa)
1257

    
1258
def pobyso_parse_string(string):
1259
    """ Legacy function. See pobyso_parse_string_sa_so. """
1260
    return pobyso_parse_string_sa_so(string)
1261
 
1262
def pobyso_parse_string_sa_so(string):
1263
    """
1264
    Get the Sollya expression computed from a Sage string or
1265
    a Sollya error object if parsing failed.
1266
    """
1267
    return sollya_lib_parse_string(string)
1268

    
1269
def pobyso_precision_so_sa(ctExpSo):
1270
    """
1271
    Computes the necessary precision to represent a number.
1272
    If x is not zero, it can be uniquely written as x = m · 2e 
1273
    where m is an odd integer and e is an integer. 
1274
    precision(x) returns the number of bits necessary to write m 
1275
    in binary (i.e. ceil(log2(m))).
1276
    """
1277
    #TODO: take care of the special case: 0, @NaN@, @Inf@
1278
    precisionSo = sollya_lib_precision(ctExpSo)
1279
    precisionSa = pobyso_constant_from_int_so_sa(precisionSo)
1280
    sollya_lib_clear_obj(precisionSo)
1281
    return precisionSa
1282
# End pobyso_precision_so_sa
1283

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

    
1450
def pobyso_polynomial_degree_so_sa(polySo):
1451
    """
1452
    Return the degree of a Sollya polynomial as a Sage int.
1453
    """
1454
    degreeSo = sollya_lib_degree(polySo)
1455
    return pobyso_constant_from_int_so_sa(degreeSo)
1456
# End pobyso_polynomial_degree_so_sa
1457

    
1458
def pobyso_polynomial_degree_so_so(polySo):
1459
    """
1460
    Thin wrapper around lib_sollya_degree().
1461
    """
1462
    return sollya_lib_degree(polySo)
1463
# End pobyso_polynomial_degree_so_so
1464
                                                                  
1465
def pobyso_range(rnLowerBound, rnUpperBound):
1466
    """ Legacy function. See pobyso_range_sa_so. """
1467
    return pobyso_range_sa_so(rnLowerBound, rnUpperBound) 
1468

    
1469
def pobyso_range_from_bounds_sa_so(rnLowerBound, rnUpperBound, precSa = None):
1470
    """
1471
    Create a Sollya range from 2 Sage real numbers as bounds
1472
    """
1473
    # TODO check precision stuff.
1474
    sollyaPrecChanged = False
1475
    (initialSollyaPrecSo, initialSollyaPrecSa) = \
1476
        pobyso_get_prec_so_so_sa()
1477
    if precSa is None:
1478
        precSa = max(rnLowerBound.parent().prec(), rnUpperBound.parent().prec())
1479
    if precSa > initialSollyaPrecSa:
1480
        if precSa <= 2:
1481
            print inspect.stack()[0][3], ": precision change <= 2 requested."
1482
        pobyso_set_prec_sa_so(precSa)
1483
        sollyaPrecChanged = True
1484
    rangeSo = sollya_lib_range_from_bounds(get_rn_value(rnLowerBound),
1485
                                           get_rn_value(rnUpperBound))
1486
    if sollyaPrecChanged:
1487
        pobyso_set_prec_so_so(initialSollyaPrecSo)
1488
        sollya_lib_clear_obj(initialSollyaPrecSo)
1489
    return rangeSo
1490
# End pobyso_range_from_bounds_sa_so
1491

    
1492
def pobyso_range_list_so_sa(listSo):
1493
    """
1494
    Return a Sollya list of ranges as a Sage list of 
1495
    floating-point intervals.
1496
    """
1497
    listSa   = []
1498
    ## The function returns none if the list is empty or an error has happened.
1499
    retVal = pobyso_get_list_elements_so_so(listSo)
1500
    if retVal is None:
1501
        return listSa
1502
    ## Just in case the interface is changed and an empty list is returned 
1503
    #  instead of None.
1504
    elif len(retVal) == 0:
1505
        return listSa
1506
    else:
1507
        ## Remember pobyso_get_list_elements_so_so returns more information
1508
        #  than just the elements of the list (# elements, is_elliptic)
1509
        listSaSo, numElements, isEndElliptic = retVal
1510
    ## Return an empty list.
1511
    if numElements == 0:
1512
        return listSa
1513
    ## Search first for the maximum precision of the elements
1514
    maxPrecSa = 0
1515
    for rangeSo in listSaSo:
1516
        #pobyso_autoprint(floatSo)
1517
        curPrecSa =  pobyso_get_prec_of_range_so_sa(rangeSo)
1518
        if curPrecSa > maxPrecSa:
1519
            maxPrecSa = curPrecSa
1520
    ##
1521
    intervalField = RealIntervalField(maxPrecSa)
1522
    ##
1523
    for rangeSo in listSaSo:
1524
        listSa.append(pobyso_range_to_interval_so_sa(rangeSo, intervalField))
1525
    return listSa
1526
# End pobyso_range_list_so_sa
1527

    
1528
def pobyso_range_max_abs_so_so(rangeSo):
1529
    """
1530
    Return, as Sollya constant, the maximum absolute value of a Sollay range.
1531
    """
1532
    lowerBoundSo = sollya_lib_inf(rangeSo)
1533
    upperBoundSo = sollya_lib_sup(rangeSo)
1534
    #
1535
    lowerBoundSo = sollya_lib_build_function_abs(lowerBoundSo)
1536
    upperBoundSo = sollya_lib_build_function_abs(upperBoundSo)
1537
    #pobyso_autoprint(lowerBoundSo)
1538
    #pobyso_autoprint(upperBoundSo)
1539
    #
1540
    maxAbsSo = sollya_lib_max(lowerBoundSo, upperBoundSo, None)
1541
    #sollya_lib_clear_obj(lowerBoundSo)
1542
    #sollya_lib_clear_obj(upperBoundSo)
1543
    return maxAbsSo
1544
# End pobyso_range_max_abs_so_so
1545
 
1546
def pobyso_range_to_interval_so_sa(rangeSo, realIntervalFieldSa = None):
1547
    """
1548
    Get a Sage interval from a Sollya range.
1549
    If no realIntervalField is given as a parameter, the Sage interval
1550
    precision is that of the Sollya range.
1551
    Otherwise, the precision is that of the realIntervalField. In this case
1552
    rounding may happen.
1553
    """
1554
    if realIntervalFieldSa is None:
1555
        precSa = pobyso_get_prec_of_range_so_sa(rangeSo)
1556
        realIntervalFieldSa = RealIntervalField(precSa)
1557
    intervalSa = \
1558
        pobyso_get_interval_from_range_so_sa(rangeSo, realIntervalFieldSa) 
1559
    return intervalSa
1560
# End pobyso_range_to_interval_so_sa
1561
#
1562
def pobyso_relative_so_so():
1563
    """
1564
    Very thin wrapper around the sollya_lib_relative function.
1565
    """
1566
    return sollya_lib_relative()
1567
# End pobyso_relative_so_so
1568
#
1569
def pobyso_rat_poly_sa_so(polySa, precSa = None):
1570
    """
1571
    Create a Sollya polynomial from a Sage rational polynomial.
1572
    """
1573
    ## TODO: filter arguments.
1574
    ## Precision. If no precision is given, use the current precision
1575
    #  of Sollya.
1576
    if precSa is None:
1577
        precSa =  pobyso_get_prec_so_sa()
1578
    #print "Precision:",  precSa
1579
    RRR = RealField(precSa)
1580
    ## Create a Sage polynomial in the "right" precision.
1581
    P_RRR = RRR[polySa.variables()[0]]
1582
    polyFloatSa = P_RRR(polySa)
1583
    ## Make sure no precision is provided: pobyso_float_poly_sa_so will
1584
    #  recover it all by itself and not make an extra conversion.
1585
    return pobyso_float_poly_sa_so(polyFloatSa)
1586
    
1587
# End pobyso_rat_poly_sa_so    
1588
    
1589
def pobyso_remez_canonical_sa_sa(func, \
1590
                                 degree, \
1591
                                 lowerBound, \
1592
                                 upperBound, \
1593
                                 weight = None, \
1594
                                 quality = None):
1595
    """
1596
    All arguments are Sage/Python.
1597
    The functions (func and weight) must be passed as expressions or strings.
1598
    Otherwise the function fails. 
1599
    The return value is a Sage polynomial.
1600
    """
1601
    var('zorglub')    # Dummy variable name for type check only. Type of 
1602
    # zorglub is "symbolic expression".
1603
    polySo = pobyso_remez_canonical_sa_so(func, \
1604
                                 degree, \
1605
                                 lowerBound, \
1606
                                 upperBound, \
1607
                                 weight, \
1608
                                 quality)
1609
    # String test
1610
    if parent(func) == parent("string"):
1611
        functionSa = eval(func)
1612
    # Expression test.
1613
    elif type(func) == type(zorglub):
1614
        functionSa = func
1615
    else:
1616
        return None
1617
    #
1618
    maxPrecision = 0
1619
    if polySo is None:
1620
        return(None)
1621
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
1622
    RRRRSa = RealField(maxPrecision)
1623
    polynomialRingSa = RRRRSa[functionSa.variables()[0]]
1624
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, RRRRSa)
1625
    polySa = polynomial(expSa, polynomialRingSa)
1626
    sollya_lib_clear_obj(polySo)
1627
    return(polySa)
1628
# End pobyso_remez_canonical_sa_sa
1629
    
1630
def pobyso_remez_canonical(func, \
1631
                           degree, \
1632
                           lowerBound, \
1633
                           upperBound, \
1634
                           weight = "1", \
1635
                           quality = None):
1636
    """ Legacy function. See pobyso_remez_canonical_sa_so. """
1637
    return(pobyso_remez_canonical_sa_so(func, \
1638
                                        degree, \
1639
                                        lowerBound, \
1640
                                        upperBound, \
1641
                                        weight, \
1642
                                        quality))
1643
# End pobyso_remez_canonical.
1644

    
1645
def pobyso_remez_canonical_sa_so(func, \
1646
                                 degree, \
1647
                                 lowerBound, \
1648
                                 upperBound, \
1649
                                 weight = None, \
1650
                                 quality = None):
1651
    """
1652
    All arguments are Sage/Python.
1653
    The functions (func and weight) must be passed as expressions or strings.
1654
    Otherwise the function fails. 
1655
    The return value is a pointer to a Sollya function.
1656
    lowerBound and upperBound mus be reals.
1657
    """
1658
    var('zorglub')    # Dummy variable name for type check only. Type of
1659
    # zorglub is "symbolic expression".
1660
    currentVariableNameSa = None
1661
    # The func argument can be of different types (string, 
1662
    # symbolic expression...)
1663
    if parent(func) == parent("string"):
1664
        localFuncSa = sage_eval(func,globals())
1665
        if len(localFuncSa.variables()) > 0:
1666
            currentVariableNameSa = localFuncSa.variables()[0]
1667
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1668
            functionSo = \
1669
              sollya_lib_parse_string(localFuncSa._assume_str().replace('_SAGE_VAR_', ''))
1670
    # Expression test.
1671
    elif type(func) == type(zorglub):
1672
        # Until we are able to translate Sage expressions into Sollya 
1673
        # expressions : parse the string version.
1674
        if len(func.variables()) > 0:
1675
            currentVariableNameSa = func.variables()[0]
1676
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1677
            functionSo = \
1678
              sollya_lib_parse_string(func._assume_str().replace('_SAGE_VAR_', ''))
1679
    else:
1680
        return(None)
1681
    if weight is None: # No weight given -> 1.
1682
        weightSo = pobyso_constant_1_sa_so()
1683
    elif parent(weight) == parent("string"): # Weight given as string: parse it.
1684
        weightSo = sollya_lib_parse_string(func)
1685
    elif type(weight) == type(zorglub): # Weight given as symbolice expression.
1686
        functionSo = \
1687
          sollya_lib_parse_string_sa_so(weight._assume_str().replace('_SAGE_VAR_', ''))
1688
    else:
1689
        return(None)
1690
    degreeSo = pobyso_constant_from_int(degree)
1691
    rangeSo = pobyso_bounds_to_range_sa_so(lowerBound, upperBound)
1692
    if not quality is None:
1693
        qualitySo= pobyso_constant_sa_so(quality)
1694
    else:
1695
        qualitySo = None
1696
        
1697
    remezPolySo = sollya_lib_remez(functionSo, \
1698
                                   degreeSo, \
1699
                                   rangeSo, \
1700
                                   weightSo, \
1701
                                   qualitySo, \
1702
                                   None)
1703
    sollya_lib_clear_obj(functionSo)
1704
    sollya_lib_clear_obj(degreeSo)
1705
    sollya_lib_clear_obj(rangeSo)
1706
    sollya_lib_clear_obj(weightSo)
1707
    if not qualitySo is None:
1708
        sollya_lib_clear_obj(qualitySo)
1709
    return(remezPolySo)
1710
# End pobyso_remez_canonical_sa_so
1711

    
1712
def pobyso_remez_canonical_so_so(funcSo, \
1713
                                 degreeSo, \
1714
                                 rangeSo, \
1715
                                 weightSo = pobyso_constant_1_sa_so(),\
1716
                                 qualitySo = None):
1717
    """
1718
    All arguments are pointers to Sollya objects.
1719
    The return value is a pointer to a Sollya function.
1720
    """
1721
    if not sollya_lib_obj_is_function(funcSo):
1722
        return(None)
1723
    return(sollya_lib_remez(funcSo, degreeSo, rangeSo, weightSo, qualitySo, None))
1724
# End pobyso_remez_canonical_so_so.
1725

    
1726
def pobyso_remez_exponents_list_sa_so(func,     \
1727
                                 exponentsList, \
1728
                                 lowerBound,    \
1729
                                 upperBound,    \
1730
                                 weight = None, \
1731
                                 quality = None):
1732
    """
1733
    All arguments are Sage/Python.
1734
    The functions (func and weight) must be passed as expressions or strings.
1735
    Otherwise the function fails. 
1736
    The return value is a pointer to a Sollya function.
1737
    lowerBound and upperBound mus be reals.
1738
    """
1739
    var('zorglub')    # Dummy variable name for type check only. Type of
1740
    # zorglub is "symbolic expression".
1741
    currentVariableNameSa = None
1742
    # The func argument can be of different types (string, 
1743
    # symbolic expression...)
1744
    if parent(func) == parent("string"):
1745
        localFuncSa = sage_eval(func,globals())
1746
        if len(localFuncSa.variables()) > 0:
1747
            currentVariableNameSa = localFuncSa.variables()[0]
1748
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1749
            functionSo = \
1750
              sollya_lib_parse_string(localFuncSa._assume_str().replace('_SAGE_VAR_', ''))
1751
    # Expression test.
1752
    elif type(func) == type(zorglub):
1753
        # Until we are able to translate Sage expressions into Sollya 
1754
        # expressions : parse the string version.
1755
        if len(func.variables()) > 0:
1756
            currentVariableNameSa = func.variables()[0]
1757
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1758
            functionSo = \
1759
              sollya_lib_parse_string(func._assume_str().replace('_SAGE_VAR_', ''))
1760
    else:
1761
        return(None)
1762
    ## Deal with the weight, much in the same way as with the function.
1763
    if weight is None: # No weight given -> 1.
1764
        weightSo = pobyso_constant_1_sa_so()
1765
    elif parent(weight) == parent("string"): # Weight given as string: parse it.
1766
        weightSo = sollya_lib_parse_string(func)
1767
    elif type(weight) == type(zorglub): # Weight given as symbolice expression.
1768
        functionSo = \
1769
          sollya_lib_parse_string(weight._assume_str().replace('_SAGE_VAR_', '',100))
1770
    else:
1771
        return(None)
1772
    rangeSo = pobyso_bounds_to_range_sa_so(lowerBound, upperBound)
1773
    if not quality is None:
1774
        qualitySo= pobyso_constant_sa_so(quality)
1775
    else:
1776
        qualitySo = None
1777
    #
1778
    ## Tranform the Sage list of exponents into a Sollya list.
1779
    exponentsListSo = pobyso_build_list_of_ints_sa_so(*exponentsList)
1780
    remezPolySo = sollya_lib_remez(functionSo, \
1781
                                   exponentsListSo, \
1782
                                   rangeSo, \
1783
                                   weightSo, \
1784
                                   qualitySo, \
1785
                                   None)
1786
    sollya_lib_clear_obj(functionSo)
1787
    sollya_lib_clear_obj(exponentsListSo)
1788
    sollya_lib_clear_obj(rangeSo)
1789
    sollya_lib_clear_obj(weightSo)
1790
    if not qualitySo is None:
1791
        sollya_lib_clear_obj(qualitySo)
1792
    return(remezPolySo)
1793
# End pobyso_remez_exponentsList_sa_so
1794
#
1795
def pobyso_round_coefficients_so_so(polySo, truncFormatListSo):
1796
    """
1797
    A wrapper around the "classical" sollya_lib_roundcoefficients: a Sollya
1798
    polynomial and a Sollya list are given as arguments.
1799
    """
1800
    return sollya_lib_roundcoefficients(polySo, truncFormatListSo)
1801

    
1802
def pobyso_round_coefficients_progressive_so_so(polySo, 
1803
                                                funcSo,
1804
                                                precSo,
1805
                                                intervalSo,
1806
                                                icSo,
1807
                                                currentApproxErrorSo,
1808
                                                approxAccurSo,
1809
                                                debug=False):
1810
    if debug:
1811
        print "Input arguments:"
1812
        print "Polynomial: ", ; pobyso_autoprint(polySo)
1813
        print "Function: ", ; pobyso_autoprint(funcSo)
1814
        print "Internal precision: ", ; pobyso_autoprint(precSo)
1815
        print "Interval: ", ; pobyso_autoprint(intervalSo)
1816
        print "Current approximation error: ", ; pobyso_autoprint(currentApproxErrorSo)
1817
        print "Requested approxiation error: ", ; pobyso_autoprint(approxAccurSo)
1818
        print "________________"
1819
    approxAccurSa        = pobyso_get_constant_as_rn_so_sa(approxAccurSo)
1820
    currentApproxErrorSa = pobyso_get_constant_as_rn_so_sa(currentApproxErrorSo)
1821
    ## If the current approximation error is too close to the target, there is
1822
    #  no possible gain.
1823
    if currentApproxErrorSa >= approxAccurSa / 2:
1824
        #### Do not return the initial argument but copies: caller may free 
1825
        #    as inutile after call.
1826
        return (sollya_lib_copy_obj(polySo),
1827
                sollya_lib_copy_obj(currentApproxErrorSo))
1828
    degreeSa             = pobyso_polynomial_degree_so_sa(polySo)
1829
    intervalSa           = pobyso_range_to_interval_so_sa(intervalSo)
1830
    
1831
    if debug:
1832
        print "degreeSa              :", degreeSa
1833
        print "intervalSa            :", intervalSa.str(style='brackets')
1834
        print "currentApproxErrorSa  :", currentApproxErrorSa 
1835
        print "approxAccurSa         :", approxAccurSa 
1836
    ### Start with a 0 value expression.
1837
    radiusSa = intervalSa.absolute_diameter() / 2
1838
    if debug:
1839
        print "log2(radius):", RR(radiusSa).log2()
1840
    iterIndex = 0
1841
    while True: 
1842
        resPolySo = pobyso_constant_0_sa_so()
1843
        roundedPolyApproxAccurSa = approxAccurSa / 2
1844
        currentRadiusPowerSa = 1 
1845
        for degree in xrange(0,degreeSa + 1):
1846
            #### At round 0, use the agressive formula. At round 1, run the
1847
            #    proved formula.
1848
            if iterIndex == 0:
1849
                roundingPowerSa = \
1850
                    floor(((currentRadiusPowerSa/roundedPolyApproxAccurSa)*(degree+1)).log2())
1851
            else:
1852
                roundingPowerSa = \
1853
                    floor(((currentRadiusPowerSa/roundedPolyApproxAccurSa)*(degreeSa+1)).log2())
1854
            ## Under extreme conditions the above formulas can evaluate under 2, which is the
1855
            #  minimal precision of an MPFR number.
1856
            if roundingPowerSa < 2:
1857
                roundingPowerSa = 2
1858
            if debug:
1859
                print "roundedPolyApproxAccurSa", roundedPolyApproxAccurSa
1860
                print "currentRadiusPowerSa", currentRadiusPowerSa
1861
                print "Current rounding exponent:", roundingPowerSa
1862
            currentRadiusPowerSa *= radiusSa
1863
            index1So = pobyso_constant_from_int_sa_so(degree)
1864
            index2So = pobyso_constant_from_int_sa_so(degree)
1865
            ### Create a monomial with:
1866
            #   - the coefficient in the initial monomial at the current degrree;
1867
            #   - the current exponent;
1868
            #   - the free variable.
1869
            cmonSo  = \
1870
                sollya_lib_build_function_mul(sollya_lib_coeff(polySo, index1So),
1871
                                              sollya_lib_build_function_pow( \
1872
                                                sollya_lib_build_function_free_variable(), \
1873
                                                index2So))
1874
            roundingPowerSo = pobyso_constant_from_int_sa_so(roundingPowerSa)
1875
            cmonrSo = pobyso_round_coefficients_single_so_so(cmonSo, roundingPowerSo)
1876
            sollya_lib_clear_obj(cmonSo)
1877
            ### Add to the result polynomial.
1878
            resPolySo = sollya_lib_build_function_add(resPolySo,
1879
                                                      cmonrSo)
1880
        # End for.
1881
        ### Check the new polynomial.
1882
        freeVarSo     = sollya_lib_build_function_free_variable()
1883
        changeVarSo   = sollya_lib_sub(freeVarSo, icSo)
1884
        resPolyCvSo   = sollya_lib_evaluate(resPolySo, changeVarSo)
1885
        errFuncSo     = sollya_lib_build_function_sub(sollya_lib_copy_obj(funcSo), 
1886
                                                      resPolyCvSo)
1887
        infNormSo = sollya_lib_dirtyinfnorm(errFuncSo, intervalSo)
1888
        ### This also clears resPolyCvSo.
1889
        sollya_lib_clear_obj(errFuncSo)
1890
        cerrSa    = pobyso_get_constant_as_rn_so_sa(infNormSo)
1891
        if debug:
1892
            print "Error of the new polynomial:", cerrSa
1893
        ### If at round 1, return the initial polynomial error. This should
1894
        #   never happen since the rounding algorithm is proved. But some 
1895
        #   circumstances may break it (e.g. internal precision of tools).
1896
        if cerrSa > approxAccurSa:
1897
            if iterIndex > 0: # Round 1 and beyond.
1898
                sollya_lib_clear_obj(resPolySo)
1899
                sollya_lib_clear_obj(infNormSo)
1900
                #### Do not return the arguments but copies: the caller may free
1901
                #    free the former as inutile after call.
1902
                return (sollya_lib_copy_obj(polySo), 
1903
                        sollya_lib_copy_obj(currentApproxErrorSo))
1904
            else: # Round 0, got round 1
1905
                sollya_lib_clear_obj(resPolySo)
1906
                sollya_lib_clear_obj(infNormSo)
1907
                iterIndex += 1
1908
                continue
1909
        ### If get here it is because cerrSa <= approxAccurSa
1910
        ### Approximation error of the new polynomial is acceptable.
1911
        return (resPolySo, infNormSo)
1912
    # End while True
1913
# End pobyso_round_coefficients_progressive_so_so
1914
    
1915
def pobyso_round_coefficients_single_so_so(polySo, commonPrecSo):
1916
    """
1917
    Create a rounded coefficients polynomial from polynomial argument to
1918
    the number of bits in size argument.
1919
    All coefficients are set to the same precision.
1920
    """
1921
    ## TODO: check arguments.
1922
    endEllipListSo = pobyso_build_end_elliptic_list_so_so(commonPrecSo)
1923
    polySo = sollya_lib_roundcoefficients(polySo, endEllipListSo, None)
1924
    sollya_lib_clear_obj(endEllipListSo)
1925
    #sollya_lib_clear_obj(endEllipListSo)
1926
    return polySo
1927
    
1928
# End pobyso_round_coefficients_single_so_so
1929

    
1930
def pobyso_set_canonical_off():
1931
    sollya_lib_set_canonical(sollya_lib_off())
1932

    
1933
def pobyso_set_canonical_on():
1934
    sollya_lib_set_canonical(sollya_lib_on())
1935

    
1936
def pobyso_set_prec(p):
1937
    """ Legacy function. See pobyso_set_prec_sa_so. """
1938
    pobyso_set_prec_sa_so(p)
1939

    
1940
def pobyso_set_prec_sa_so(p):
1941
    #a = c_int(p)
1942
    #precSo = c_void_p(sollya_lib_constant_from_int(a))
1943
    #precSo = sollya_lib_constant_from_int(a)
1944
    precSo = pobyso_constant_from_int_sa_so(p)
1945
    sollya_lib_set_prec(precSo)
1946
    sollya_lib_clear_obj(precSo)
1947
# End pobyso_set_prec_sa_so.
1948

    
1949
def pobyso_set_prec_so_so(newPrecSo):
1950
    sollya_lib_set_prec(newPrecSo)
1951
# End pobyso_set_prec_so_so.
1952

    
1953
def pobyso_inf_so_so(intervalSo):
1954
    """
1955
    Very thin wrapper around sollya_lib_inf().
1956
    """
1957
    return sollya_lib_inf(intervalSo)
1958
# End pobyso_inf_so_so.
1959
#   
1960
def pobyso_supnorm_so_sa(polySo, funcSo, intervalSo, errorTypeSo = None,\
1961
                         accuracySo = None, realFieldSa = None):
1962
    """
1963
    Computes the supremum norm from Solly input arguments and returns a
1964
    Sage floating-point number whose precision is set by the only Sage argument.
1965
    
1966
    The returned value is the maximum of the absolute values of the range
1967
    elements  returned by the Sollya supnorm functions
1968
    """
1969
    supNormRangeSo = pobyso_supnorm_so_so(polySo, 
1970
                                          funcSo, 
1971
                                          intervalSo, 
1972
                                          errorTypeSo,
1973
                                          accuracySo)
1974
    supNormSo = pobyso_range_max_abs_so_so(supNormRangeSo)
1975
    sollya_lib_clear_obj(supNormRangeSo)
1976
    #pobyso_autoprint(supNormSo)
1977
    supNormSa = pobyso_get_constant_as_rn_with_rf_so_sa(supNormSo, realFieldSa)
1978
    sollya_lib_clear_obj(supNormSo)
1979
    return supNormSa 
1980
# End pobyso_supnorm_so_sa.
1981
#
1982
def pobyso_supnorm_so_so(polySo, funcSo, intervalSo, errorTypeSo = None,\
1983
                         accuracySo = None):
1984
    """
1985
    Computes the supnorm of the approximation error between the given 
1986
    polynomial and function. Attention: returns a range!
1987
    errorTypeSo defaults to "absolute".
1988
    accuracySo defaults to 2^(-40).
1989
    """
1990
    if errorTypeSo is None:
1991
        errorTypeSo = sollya_lib_absolute(None)
1992
        errorTypeIsNone = True
1993
    else:
1994
        errorTypeIsNone = False
1995
    #
1996
    if accuracySo is None:
1997
        # Notice the **: we are in Pythonland here!
1998
        accuracySo = pobyso_constant_sa_so(RR(2**(-40)))
1999
        accuracyIsNone = True
2000
    else:
2001
        accuracyIsNone = False
2002
    #pobyso_autoprint(accuracySo)
2003
    resultSo = \
2004
        sollya_lib_supnorm(polySo, funcSo, intervalSo, errorTypeSo, \
2005
                              accuracySo)
2006
    if errorTypeIsNone:
2007
        sollya_lib_clear_obj(errorTypeSo)
2008
    if accuracyIsNone:
2009
        sollya_lib_clear_obj(accuracySo)
2010
    return resultSo
2011
# End pobyso_supnorm_so_so
2012
#
2013
def pobyso_taylor_expansion_no_change_var_so_so(functionSo, 
2014
                                                degreeSo, 
2015
                                                rangeSo,
2016
                                                errorTypeSo=None, 
2017
                                                sollyaPrecSo=None):
2018
    """
2019
    Compute the Taylor expansion without the variable change
2020
    x -> x-intervalCenter.
2021
    """
2022
    # Change internal Sollya precision, if needed.
2023
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
2024
    sollyaPrecChanged = False
2025
    if sollyaPrecSo is None:
2026
        pass
2027
    else:
2028
        sollya_lib_set_prec(sollyaPrecSo)
2029
        sollyaPrecChanged = True
2030
    # Error type stuff: default to absolute.
2031
    if errorTypeSo is None:
2032
        errorTypeIsNone = True
2033
        errorTypeSo = sollya_lib_absolute(None)
2034
    else:
2035
        errorTypeIsNone = False
2036
    intervalCenterSo = sollya_lib_mid(rangeSo, None)
2037
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo,
2038
                                         intervalCenterSo,
2039
                                         rangeSo, errorTypeSo, None)
2040
    # taylorFormListSaSo is a Python list of Sollya objects references that 
2041
    # are copies of the elements of taylorFormSo.
2042
    # pobyso_get_list_elements_so_so clears taylorFormSo.
2043
    (taylorFormListSaSo, numElementsSa, isEndEllipticSa) = \
2044
        pobyso_get_list_elements_so_so(taylorFormSo)
2045
    polySo = sollya_lib_copy_obj(taylorFormListSaSo[0])
2046
    #print "Num elements:", numElementsSa
2047
    sollya_lib_clear_obj(taylorFormSo)
2048
    #polySo = taylorFormListSaSo[0]
2049
    #errorRangeSo = sollya_lib_copy_obj(taylorFormListSaSo[2])
2050
    errorRangeSo = taylorFormListSaSo[2]
2051
    # No copy_obj needed here: a new objects are created.
2052
    maxErrorSo    = sollya_lib_sup(errorRangeSo)
2053
    minErrorSo    = sollya_lib_inf(errorRangeSo)
2054
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
2055
    absMinErrorSo = sollya_lib_abs(minErrorSo)
2056
    sollya_lib_clear_obj(maxErrorSo)
2057
    sollya_lib_clear_obj(minErrorSo)
2058
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
2059
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
2060
    # If changed, reset the Sollya working precision.
2061
    if sollyaPrecChanged:
2062
        sollya_lib_set_prec(initialSollyaPrecSo)
2063
    sollya_lib_clear_obj(initialSollyaPrecSo)
2064
    if errorTypeIsNone:
2065
        sollya_lib_clear_obj(errorTypeSo)
2066
    pobyso_clear_taylorform_sa_so(taylorFormListSaSo)
2067
    if absMaxErrorSa > absMinErrorSa:
2068
        sollya_lib_clear_obj(absMinErrorSo)
2069
        return (polySo, intervalCenterSo, absMaxErrorSo)
2070
    else:
2071
        sollya_lib_clear_obj(absMaxErrorSo)
2072
        return (polySo, intervalCenterSo, absMinErrorSo)
2073
# end pobyso_taylor_expansion_no_change_var_so_so
2074

    
2075
def pobyso_taylor_expansion_with_change_var_so_so(functionSo, degreeSo, \
2076
                                                  rangeSo, \
2077
                                                  errorTypeSo=None, \
2078
                                                  sollyaPrecSo=None):
2079
    """
2080
    Compute the Taylor expansion with the variable change
2081
    x -> (x-intervalCenter) included.
2082
    """
2083
    # Change Sollya internal precision, if need.
2084
    sollyaPrecChanged = False
2085
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_sos_sa()
2086
    if sollyaPrecSo is None:
2087
        pass
2088
    else:
2089
        sollya_lib_set_prec(sollyaPrecSo)
2090
        sollyaPrecChanged = True
2091
    #
2092
    # Error type stuff: default to absolute.
2093
    if errorTypeSo is None:
2094
        errorTypeIsNone = True
2095
        errorTypeSo = sollya_lib_absolute(None)
2096
    else:
2097
        errorTypeIsNone = False
2098
    intervalCenterSo = sollya_lib_mid(rangeSo)
2099
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo, \
2100
                                         intervalCenterSo, \
2101
                                         rangeSo, errorTypeSo, None)
2102
    # taylorFormListSaSo is a Python list of Sollya objects references that 
2103
    # are copies of the elements of taylorFormSo.
2104
    # pobyso_get_list_elements_so_so clears taylorFormSo.
2105
    (taylorFormListSo, numElements, isEndElliptic) = \
2106
        pobyso_get_list_elements_so_so(taylorFormSo)
2107
    polySo = taylorFormListSo[0]
2108
    errorRangeSo = taylorFormListSo[2]
2109
    maxErrorSo    = sollya_lib_sup(errorRangeSo)
2110
    minErrorSo    = sollya_lib_inf(errorRangeSo)
2111
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
2112
    absMinErrorSo = sollya_lib_abs(minErrorSo)
2113
    sollya_lib_clear_obj(maxErrorSo)
2114
    sollya_lib_clear_obj(minErrorSo)
2115
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
2116
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
2117
    changeVarExpSo = sollya_lib_build_function_sub(\
2118
                       sollya_lib_build_function_free_variable(),\
2119
                       sollya_lib_copy_obj(intervalCenterSo))
2120
    polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpSo)
2121
    sollya_lib_clear_obj(polySo) 
2122
    sollya_lib_clear_obj(changeVarExpSo)
2123
    # If changed, reset the Sollya working precision.
2124
    if sollyaPrecChanged:
2125
        sollya_lib_set_prec(initialSollyaPrecSo)
2126
    sollya_lib_clear_obj(initialSollyaPrecSo)
2127
    if errorTypeIsNone:
2128
        sollya_lib_clear_obj(errorTypeSo)
2129
    sollya_lib_clear_obj(taylorFormSo)
2130
    # Do not clear maxErrorSo.
2131
    if absMaxErrorSa > absMinErrorSa:
2132
        sollya_lib_clear_obj(absMinErrorSo)
2133
        return((polyVarChangedSo, intervalCenterSo, absMaxErrorSo))
2134
    else:
2135
        sollya_lib_clear_obj(absMaxErrorSo)
2136
        return((polyVarChangedSo, intervalCenterSo, absMinErrorSo))
2137
# end pobyso_taylor_expansion_with_change_var_so_so
2138

    
2139
def pobyso_taylor(function, degree, point):
2140
    """ Legacy function. See pobysoTaylor_so_so. """
2141
    return(pobyso_taylor_so_so(function, degree, point))
2142

    
2143
def pobyso_taylor_so_so(functionSo, degreeSo, pointSo):
2144
    return(sollya_lib_taylor(functionSo, degreeSo, pointSo))
2145
    
2146
def pobyso_taylorform(function, degree, point = None, 
2147
                      interval = None, errorType=None):
2148
    """ Legacy function. See pobyso_taylorform_sa_sa;"""
2149
    
2150
def pobyso_taylorform_sa_sa(functionSa, \
2151
                            degreeSa, \
2152
                            pointSa, \
2153
                            intervalSa=None, \
2154
                            errorTypeSa=None, \
2155
                            precisionSa=None):
2156
    """
2157
    Compute the Taylor form of 'degreeSa' for 'functionSa' at 'pointSa' 
2158
    for 'intervalSa' with 'errorTypeSa' (a string) using 'precisionSa'. 
2159
    point: must be a Real or a Real interval.
2160
    return the Taylor form as an array
2161
    TODO: take care of the interval and of the point when it is an interval;
2162
          when errorType is not None;
2163
          take care of the other elements of the Taylor form (coefficients 
2164
          errors and delta.
2165
    """
2166
    # Absolute as the default error.
2167
    if errorTypeSa is None:
2168
        errorTypeSo = sollya_lib_absolute()
2169
    elif errorTypeSa == "relative":
2170
        errorTypeSo = sollya_lib_relative()
2171
    elif errortypeSa == "absolute":
2172
        errorTypeSo = sollya_lib_absolute()
2173
    else:
2174
        # No clean up needed.
2175
        return None
2176
    # Global precision stuff
2177
    sollyaPrecisionChangedSa = False
2178
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
2179
    if precisionSa is None:
2180
        precSa = initialSollyaPrecSa
2181
    else:
2182
        if precSa > initialSollyaPrecSa:
2183
            if precSa <= 2:
2184
                print inspect.stack()[0][3], ":precision change <= 2 requested."
2185
            pobyso_set_prec_sa_so(precSa)
2186
            sollyaPrecisionChangedSa = True
2187
    #        
2188
    if len(functionSa.variables()) > 0:
2189
        varSa = functionSa.variables()[0]
2190
        pobyso_name_free_variable_sa_so(str(varSa))
2191
    # In any case (point or interval) the parent of pointSa has a precision
2192
    # method.
2193
    pointPrecSa = pointSa.parent().precision()
2194
    if precSa > pointPrecSa:
2195
        pointPrecSa = precSa
2196
    # In any case (point or interval) pointSa has a base_ring() method.
2197
    pointBaseRingString = str(pointSa.base_ring())
2198
    if re.search('Interval', pointBaseRingString) is None: # Point
2199
        pointSo = pobyso_constant_sa_so(pointSa, pointPrecSa)
2200
    else: # Interval.
2201
        pointSo = pobyso_interval_to_range_sa_so(pointSa, pointPrecSa)
2202
    # Sollyafy the function.
2203
    functionSo = pobyso_parse_string_sa_so(functionSa._assume_str().replace('_SAGE_VAR_', ''))
2204
    if sollya_lib_obj_is_error(functionSo):
2205
        print "pobyso_tailorform: function string can't be parsed!"
2206
        return None
2207
    # Sollyafy the degree
2208
    degreeSo = sollya_lib_constant_from_int(int(degreeSa))
2209
    # Sollyafy the point
2210
    # Call Sollya
2211
    taylorFormSo = \
2212
        sollya_lib_taylorform(functionSo, degreeSo, pointSo, errorTypeSo,\
2213
                                         None)
2214
    sollya_lib_clear_obj(functionSo)
2215
    sollya_lib_clear_obj(degreeSo)
2216
    sollya_lib_clear_obj(pointSo)
2217
    sollya_lib_clear_obj(errorTypeSo)
2218
    (tfsAsList, numElements, isEndElliptic) = \
2219
            pobyso_get_list_elements_so_so(taylorFormSo)
2220
    polySo = tfsAsList[0]
2221
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
2222
    polyRealField = RealField(maxPrecision)
2223
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, polyRealField)
2224
    if sollyaPrecisionChangedSa:
2225
        sollya_lib_set_prec(initialSollyaPrecSo)
2226
    sollya_lib_clear_obj(initialSollyaPrecSo)
2227
    polynomialRing = polyRealField[str(varSa)]
2228
    polySa = polynomial(expSa, polynomialRing)
2229
    taylorFormSa = [polySa]
2230
    # Final clean-up.
2231
    sollya_lib_clear_obj(taylorFormSo)
2232
    return(taylorFormSa)
2233
# End pobyso_taylor_form_sa_sa
2234

    
2235
def pobyso_taylorform_so_so(functionSo, degreeSo, pointSo, intervalSo=None, \
2236
                            errorTypeSo=None):
2237
    createdErrorType = False
2238
    if errorTypeSo is None:
2239
        errorTypeSo = sollya_lib_absolute()
2240
        createdErrorType = True
2241
    else:
2242
        #TODO: deal with the other case.
2243
        pass
2244
    if intervalSo is None:
2245
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
2246
                                         errorTypeSo, None)
2247
    else:
2248
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
2249
                                         intervalSo, errorTypeSo, None)
2250
    if createdErrorType:
2251
        sollya_lib_clear_obj(errorTypeSo)
2252
    return resultSo
2253
        
2254

    
2255
def pobyso_univar_polynomial_print_reverse(polySa):
2256
    """ Legacy function. See pobyso_univar_polynomial_print_reverse_sa_sa. """
2257
    return(pobyso_univar_polynomial_print_reverse_sa_sa(polySa))
2258

    
2259
def pobyso_univar_polynomial_print_reverse_sa_sa(polySa):
2260
    """
2261
    Return the string representation of a univariate polynomial with
2262
    monomials ordered in the x^0..x^n order of the monomials.
2263
    Remember: Sage
2264
    """
2265
    polynomialRing = polySa.base_ring()
2266
    # A very expensive solution:
2267
    # -create a fake multivariate polynomial field with only one variable,
2268
    #   specifying a negative lexicographical order;
2269
    mpolynomialRing = PolynomialRing(polynomialRing.base(), \
2270
                                     polynomialRing.variable_name(), \
2271
                                     1, order='neglex')
2272
    # - convert the univariate argument polynomial into a multivariate
2273
    #   version;
2274
    p = mpolynomialRing(polySa)
2275
    # - return the string representation of the converted form.
2276
    # There is no simple str() method defined for p's class.
2277
    return(p.__str__())
2278
#
2279
#print pobyso_get_prec()  
2280
pobyso_set_prec(165)
2281
#print pobyso_get_prec()  
2282
#a=100
2283
#print type(a)
2284
#id(a)
2285
#print "Max arity: ", pobyso_max_arity
2286
#print "Function tripleDouble (43) as a string: ", pobyso_function_type_as_string(43)
2287
#print "Function None (44) as a string: ", pobyso_function_type_as_string(44)
2288
sys.stderr.write("\t...Pobyso check done.\n")