Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 240

Historique | Voir | Annoter | Télécharger (85,11 ko)

1
"""
2
@file pobyso.py
3
Actual functions to use in Sage
4
ST 2012-11-13
5

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

9
pobyso functions come in five flavors: 
10
- the _so_so (arguments and returned objects are pointers to Sollya objects, 
11
  includes the void function and the no arguments function that return a 
12
  pointer to a Sollya object);
13
- the _so_sa (argument are pointers to Sollya objects, returned objects are
14
  Sage/Python objects or, more generally, information is transfered from the
15
  Sollya world to Sage/Python world; e.g. functions without arguments that
16
  return a Sage/Python object);
17
- the _sa_so (arguments are Sage/Python objects, returned objects are 
18
  pointers to Sollya objects);
19
- the sa_sa (arguments and returned objects are all Sage/Python objects);
20
- a catch all flavor, without any suffix, (e. g. functions that have no argument 
21
  nor return value).
22
This classification is not always very strict. Conversion functions from Sollya
23
to Sage/Python are sometimes decorated with Sage/Python arguments to set
24
the precision. These functions remain in the so_sa category.
25
NOTES:
26
Reported errors in Eclipse come from the calls to
27
the Sollya library
28

29
ToDo (among other things): 
30
 -memory management.
31
"""
32
from ctypes import *
33
import re
34
from sage.symbolic.expression_conversions import polynomial
35
from sage.symbolic.expression_conversions import PolynomialConverter
36
"""
37
Create the equivalent to an enum for the Sollya function types.
38
"""
39
(SOLLYA_BASE_FUNC_ABS,
40
SOLLYA_BASE_FUNC_ACOS,
41
    SOLLYA_BASE_FUNC_ACOSH,
42
    SOLLYA_BASE_FUNC_ADD,
43
    SOLLYA_BASE_FUNC_ASIN,
44
    SOLLYA_BASE_FUNC_ASINH,
45
    SOLLYA_BASE_FUNC_ATAN,
46
    SOLLYA_BASE_FUNC_ATANH,
47
    SOLLYA_BASE_FUNC_CEIL,
48
    SOLLYA_BASE_FUNC_CONSTANT,
49
    SOLLYA_BASE_FUNC_COS,
50
    SOLLYA_BASE_FUNC_COSH,
51
    SOLLYA_BASE_FUNC_DIV,
52
    SOLLYA_BASE_FUNC_DOUBLE,
53
    SOLLYA_BASE_FUNC_DOUBLEDOUBLE,
54
    SOLLYA_BASE_FUNC_DOUBLEEXTENDED,
55
    SOLLYA_BASE_FUNC_ERF,
56
    SOLLYA_BASE_FUNC_ERFC,
57
    SOLLYA_BASE_FUNC_EXP,
58
    SOLLYA_BASE_FUNC_EXP_M1,
59
    SOLLYA_BASE_FUNC_FLOOR,
60
    SOLLYA_BASE_FUNC_FREE_VARIABLE,
61
    SOLLYA_BASE_FUNC_HALFPRECISION,
62
    SOLLYA_BASE_FUNC_LIBRARYCONSTANT,
63
    SOLLYA_BASE_FUNC_LIBRARYFUNCTION,
64
    SOLLYA_BASE_FUNC_LOG,
65
    SOLLYA_BASE_FUNC_LOG_10,
66
    SOLLYA_BASE_FUNC_LOG_1P,
67
    SOLLYA_BASE_FUNC_LOG_2,
68
    SOLLYA_BASE_FUNC_MUL,
69
    SOLLYA_BASE_FUNC_NEARESTINT,
70
    SOLLYA_BASE_FUNC_NEG,
71
    SOLLYA_BASE_FUNC_PI,
72
    SOLLYA_BASE_FUNC_POW,
73
    SOLLYA_BASE_FUNC_PROCEDUREFUNCTION,
74
    SOLLYA_BASE_FUNC_QUAD,
75
    SOLLYA_BASE_FUNC_SIN,
76
    SOLLYA_BASE_FUNC_SINGLE,
77
    SOLLYA_BASE_FUNC_SINH,
78
    SOLLYA_BASE_FUNC_SQRT,
79
    SOLLYA_BASE_FUNC_SUB,
80
    SOLLYA_BASE_FUNC_TAN,
81
    SOLLYA_BASE_FUNC_TANH,
82
SOLLYA_BASE_FUNC_TRIPLEDOUBLE) = map(int,xrange(44))
83
print "\nSuperficial pobyso check..."    
84
#print "First constant - SOLLYA_BASE_FUNC_ABS: ", SOLLYA_BASE_FUNC_ABS
85
#print "Last constant  - SOLLYA_BASE_FUNC_TRIPLEDOUBLE: ", SOLLYA_BASE_FUNC_TRIPLEDOUBLE
86

    
87
pobyso_max_arity = 9
88

    
89
def pobyso_absolute_so_so():
90
    return(sollya_lib_absolute(None))
91

    
92
def pobyso_autoprint(arg):
93
    sollya_lib_autoprint(arg, None)
94

    
95
def pobyso_autoprint_so_so(arg):
96
    sollya_lib_autoprint(arg,None)
97
    
98
def pobyso_bounds_to_range_sa_so(rnLowerBoundSa, rnUpperBoundSa, \
99
                                 precisionSa=None):
100
    """
101
    Return a Sollya range from to 2 RealField Sage elements.
102
    The Sollya range element has a sufficient precision to hold all
103
    the digits of the widest of the Sage bounds.
104
    """
105
    # Sanity check.
106
    if rnLowerBoundSa > rnUpperBoundSa:
107
        return None
108
    # Precision stuff.
109
    if precisionSa is None:
110
        # Check for the largest precision.
111
        lbPrecSa = rnLowerBoundSa.parent().precision()
112
        ubPrecSa = rnLowerBoundSa.parent().precision()
113
        maxPrecSa = max(lbPrecSa, ubPrecSa)
114
    else:
115
        maxPrecSa = precisionSa
116
    # From Sage to Sollya bounds.
117
#    lowerBoundSo = sollya_lib_constant(get_rn_value(rnLowerBoundSa),
118
#                                       maxPrecSa)
119
    lowerBoundSo = pobyso_constant_sa_so(rnLowerBoundSa,
120
                                         maxPrecSa)
121
    upperBoundSo = pobyso_constant_sa_so(rnUpperBoundSa,
122
                                       maxPrecSa)
123
    
124
    # From Sollya bounds to range.
125
    rangeSo = sollya_lib_range(lowerBoundSo, upperBoundSo)
126
    # Back to original precision.
127
    # Clean up
128
    sollya_lib_clear_obj(lowerBoundSo)
129
    sollya_lib_clear_obj(upperBoundSo)
130
    return rangeSo
131
# End pobyso_bounds_to_range_sa_so
132

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

    
155
def pobyso_build_function_sub_so_so(exp1So, exp2So):
156
    return sollya_lib_build_function_sub(exp1So, exp2So)
157

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

    
193
def pobyso_change_var_in_function_so_so(funcSo, chvarExpSo):
194
    """
195
    Variable change in a function.
196
    """
197
    return(sollya_lib_evaluate(funcSo,chvarExpSo))
198
# End pobyso_change_var_in_function_so_so     
199

    
200
def pobyso_chebyshevform_so_so(functionSo, degreeSo, intervalSo):
201
    resultSo = sollya_lib_chebyshevform(functionSo, degreeSo, intervalSo)
202
    return(resultSo)
203
# End pobyso_chebyshevform_so_so.
204

    
205
def pobyso_clear_obj(objSo):
206
    """
207
    Free a Sollya object's memory.
208
    Very thin wrapper around sollya_lib_clear_obj().
209
    """
210
    sollya_lib_clear_obj(objSo)
211
# End pobyso_clear_obj. 
212

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

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

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

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

    
314

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

    
343
def pobyso_constant_1():
344
    """
345
    Obvious.
346
    Legacy function. See pobyso_constant_so_so. 
347
    """
348
    return pobyso_constant_1_sa_so()
349

    
350
def pobyso_constant_1_sa_so():
351
    """
352
    Obvious.
353
    """
354
    return(pobyso_constant_from_int_sa_so(1))
355

    
356
def pobyso_constant_from_int(anInt):
357
    """ Legacy function. See pobyso_constant_from_int_sa_so. """
358
    return pobyso_constant_from_int_sa_so(anInt)
359

    
360
def pobyso_constant_from_int_sa_so(anInt):
361
    """
362
    Get a Sollya constant from a Sage int.
363
    """
364
    return sollya_lib_constant_from_int64(long(anInt))
365

    
366
def pobyso_constant_from_int_so_sa(constSo):
367
    """
368
    Get a Sage int from a Sollya int constant.
369
    Usefull for precision or powers in polynomials.
370
    """
371
    constSa = c_long(0)
372
    sollya_lib_get_constant_as_int64(byref(constSa), constSo)
373
    return constSa.value
374
# End pobyso_constant_from_int_so_sa
375

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

    
390
def pobyso_constant_sollya_prec_sa_so(rnArgSa):
391
    """
392
    Create a Sollya constant from a Sage RealNumber at the
393
    current precision in Sollya.
394
    """
395
    currentSollyaPrecSa = pobyso_get_prec_so_sa()
396
    return pobyso_constant_sa_so(rnArgSa, currentSollyaPrecSa)
397
# End pobyso_constant_sollya_prec_sa_so
398

    
399
def pobyso_end_elliptic_list_so_sa_so(objectsListSo, intCountSa):
400
    """
401
    Create a Sollya end elliptic list made of the objectListSo[0] to
402
     objectsListSo[intCountSa-1] objects.
403
    """     
404
    return sollya_lib_end_elliptic_list(objectSo, int(intCountSa))
405

    
406
def pobyso_error_so():
407
    return sollya_lib_error(None)
408
# End pobyso_error().
409

    
410
def pobyso_evaluate_so_so(funcSo, argumentSo):
411
    """
412
    Evaluates funcSo for arguemntSo through sollya_lib_evaluate().
413
    """
414
    return sollya_lib_evaluate(funcSo, argumentSo)
415
# End pobyso_evaluate_so_so.
416
#
417
def pobyso_diff_so_so(funcSo):
418
    """
419
    Very thin wrapper around sollya_lib_diff.
420
    """
421
    ## TODO: add a check to make sure funcSo is a functional expression.
422
    return sollya_lib_diff(funcSo)
423

    
424
def pobyso_dirty_find_zeros_so_so(funcSo, rangeSo):
425
    """
426
    Thin wrapper over sollya_lib_dirtyfindzeros()
427
    """
428
    return sollya_lib_dirtyfindzeros(funcSo, rangeSo)
429
# End pobys_dirty_find_zeros
430

    
431
def pobyso_dirty_inf_norm_so_so(funcSo, rangeSo, preSo=None):
432
    """
433
    Thin wrapper around sollya_dirtyinfnorm().
434
    """
435
    # TODO: manage the precision change.
436
    
437
    return sollya_lib_dirtyinfnorm(funcSo, rangeSo)
438
# End pobyso_dirty_inf_norm_so_so
439

    
440
def pobyso_float_list_so_sa(listSo):
441
    """
442
    Return a Sollya list of floating-point numbers as a Sage list of 
443
    floating-point numbers.
444
    TODO: add a test to make sure that each element of the list is a constant.
445
    """
446
    listSa   = []
447
    ## The function returns none if the list is empty or an error has happened.
448
    retVal = pobyso_get_list_elements_so_so(listSo)
449
    if retVal is None:
450
        return listSa
451
    ## Just in case the interface is changed and an empty list is returned 
452
    #  instead of None.
453
    elif len(retVal) == 0:
454
        return listSa
455
    else:
456
        ## Remember pobyso_get_list_elements_so_so returns more information
457
        #  than just the elements of the list (# elements, is_elliptic)
458
        listSaSo, numElements, isEndElliptic = retVal
459
    ## Return an empty list.
460
    if numElements == 0:
461
        return listSa
462
    ## Search first for the maximum precision of the elements
463
    maxPrecSa = 0
464
    for floatSo in listSaSo:
465
        #pobyso_autoprint(floatSo)
466
        curPrecSa =  pobyso_get_prec_of_constant_so_sa(floatSo)
467
        if curPrecSa > maxPrecSa:
468
            maxPrecSa = curPrecSa
469
    ##
470
    RF = RealField(maxPrecSa)
471
    ##
472
    for floatSo in listSaSo:
473
        listSa.append(pobyso_get_constant_as_rn_with_rf_so_sa(floatSo))
474
    return listSa
475
# End pobyso_float_list_so_sa
476

    
477
def pobyso_float_poly_sa_so(polySa, precSa = None):
478
    """
479
    Create a Sollya polynomial from a Sage RealField polynomial.
480
    """
481
    ## TODO: filter arguments.
482
    ## Precision. If a precision is given, convert the polynomial
483
    #  into the right polynomial field. If not convert it straight
484
    #  to Sollya.
485
    sollyaPrecChanged = False 
486
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
487
    if precSa is None:
488
        precSa = polySa.parent().base_ring().precision()
489
    if (precSa > initialSollyaPrecSa):
490
        assert precSa >= 2, "Precision change <2 requested"
491
        if precSa <= 2:
492
            print inspect.stack()[0][3], ": precision change <= 2 requested"
493
        precSo = pobyso_constant_from_int(precSa)
494
        pobyso_set_prec_so_so(precSo)
495
        sollya_lib_clear_obj(precSo)
496
        sollyaPrecChanged = True    
497
    ## Get exponents and coefficients.
498
    exponentsSa     = polySa.exponents()
499
    coefficientsSa  = polySa.coefficients()
500
    ## Build the polynomial.
501
    polySo = None
502
    for coefficientSa, exponentSa in zip(coefficientsSa, exponentsSa):
503
        #print coefficientSa.n(prec=precSa), exponentSa
504
        coefficientSo = \
505
            pobyso_constant_sa_so(coefficientSa)
506
        #pobyso_autoprint(coefficientSo)
507
        exponentSo = \
508
            pobyso_constant_from_int_sa_so(exponentSa)
509
        #pobyso_autoprint(exponentSo)
510
        monomialSo = sollya_lib_build_function_pow(
511
                       sollya_lib_build_function_free_variable(),
512
                       exponentSo)
513
        polyTermSo = sollya_lib_build_function_mul(coefficientSo,
514
                                                       monomialSo)
515
        if polySo is None:
516
            polySo = polyTermSo
517
        else:
518
            polySo = sollya_lib_build_function_add(polySo, polyTermSo)
519
    if sollyaPrecChanged:
520
        pobyso_set_prec_so_so(initialSollyaPrecSo)
521
        sollya_lib_clear_obj(initialSollyaPrecSo)
522
    return polySo
523
# End pobyso_float_poly_sa_so    
524

    
525
def pobyso_float_poly_so_sa(polySo, realFieldSa=None):
526
    """
527
    Convert a Sollya polynomial into a Sage floating-point polynomial.
528
    If no realField is given, a RealField corresponding to the maximum 
529
    precision of the coefficients is internally computed.
530
    The real field is not returned but can be easily retrieved from 
531
    the polynomial itself.
532
    ALGORITHM:
533
    - (optional) compute the RealField of the coefficients;
534
    - convert the Sollya expression into a Sage expression;
535
    - convert the Sage expression into a Sage polynomial
536
    """    
537
    if realFieldSa is None:
538
        expressionPrecSa = pobyso_get_max_prec_of_exp_so_sa(polySo)
539
        #print "Maximum precision of Sollya polynomial coefficients:", expressionPrecSa
540
        if expressionPrecSa < 2 or expressionPrecSa > 2147483391:
541
            print "Maximum degree of expression:", expressionPrecSa
542
        realFieldSa      = RealField(expressionPrecSa)
543
    #print "Sollya expression before...",
544
    #pobyso_autoprint(polySo)
545

    
546
    expressionSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo,
547
                                                             realFieldSa)
548
    #print "...Sollya expression after."
549
    #pobyso_autoprint(polySo)
550
    polyVariableSa = expressionSa.variables()[0]
551
    polyRingSa     = realFieldSa[str(polyVariableSa)]
552
    #print polyRingSa
553
    # Do not use the polynomial(expressionSa, ring=polyRingSa) form!
554
    polynomialSa = polyRingSa(expressionSa)
555
    polyCoeffsListSa = polynomialSa.coefficients()
556
    #for coeff in polyCoeffsListSa:
557
    #    print coeff.abs().n()
558
    return polynomialSa
559
# End pobyso_float_poly_so_sa
560

    
561
def pobyso_free_variable():
562
    """
563
    Ultra thin wrapper around the sollya_lib_function_build_free_variable function.
564
    """
565
    return sollya_lib_build_function_free_variable()
566
    
567
def pobyso_function_type_as_string(funcType):
568
    """ Legacy function. See pobyso_function_type_as_string_so_sa. """
569
    return(pobyso_function_type_as_string_so_sa(funcType))
570

    
571
def pobyso_function_type_as_string_so_sa(funcType):
572
    """
573
    Numeric Sollya function codes -> Sage mathematical function names.
574
    Notice that pow -> ^ (a la Sage, not a la Python).
575
    """
576
    if funcType == SOLLYA_BASE_FUNC_ABS:
577
        return "abs"
578
    elif funcType == SOLLYA_BASE_FUNC_ACOS:
579
        return "arccos"
580
    elif funcType == SOLLYA_BASE_FUNC_ACOSH:
581
        return "arccosh"
582
    elif funcType == SOLLYA_BASE_FUNC_ADD:
583
        return "+"
584
    elif funcType == SOLLYA_BASE_FUNC_ASIN:
585
        return "arcsin"
586
    elif funcType == SOLLYA_BASE_FUNC_ASINH:
587
        return "arcsinh"
588
    elif funcType == SOLLYA_BASE_FUNC_ATAN:
589
        return "arctan"
590
    elif funcType == SOLLYA_BASE_FUNC_ATANH:
591
        return "arctanh"
592
    elif funcType == SOLLYA_BASE_FUNC_CEIL:
593
        return "ceil"
594
    elif funcType == SOLLYA_BASE_FUNC_CONSTANT:
595
        return "cte"
596
    elif funcType == SOLLYA_BASE_FUNC_COS:
597
        return "cos"
598
    elif funcType == SOLLYA_BASE_FUNC_COSH:
599
        return "cosh"
600
    elif funcType == SOLLYA_BASE_FUNC_DIV:
601
        return "/"
602
    elif funcType == SOLLYA_BASE_FUNC_DOUBLE:
603
        return "double"
604
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEDOUBLE:
605
        return "doubleDouble"
606
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEEXTENDED:
607
        return "doubleDxtended"
608
    elif funcType == SOLLYA_BASE_FUNC_ERF:
609
        return "erf"
610
    elif funcType == SOLLYA_BASE_FUNC_ERFC:
611
        return "erfc"
612
    elif funcType == SOLLYA_BASE_FUNC_EXP:
613
        return "exp"
614
    elif funcType == SOLLYA_BASE_FUNC_EXP_M1:
615
        return "expm1"
616
    elif funcType == SOLLYA_BASE_FUNC_FLOOR:
617
        return "floor"
618
    elif funcType == SOLLYA_BASE_FUNC_FREE_VARIABLE:
619
        return "freeVariable"
620
    elif funcType == SOLLYA_BASE_FUNC_HALFPRECISION:
621
        return "halfPrecision"
622
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYCONSTANT:
623
        return "libraryConstant"
624
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYFUNCTION:
625
        return "libraryFunction"
626
    elif funcType == SOLLYA_BASE_FUNC_LOG:
627
        return "log"
628
    elif funcType == SOLLYA_BASE_FUNC_LOG_10:
629
        return "log10"
630
    elif funcType == SOLLYA_BASE_FUNC_LOG_1P:
631
        return "log1p"
632
    elif funcType == SOLLYA_BASE_FUNC_LOG_2:
633
        return "log2"
634
    elif funcType == SOLLYA_BASE_FUNC_MUL:
635
        return "*"
636
    elif funcType == SOLLYA_BASE_FUNC_NEARESTINT:
637
        return "round"
638
    elif funcType == SOLLYA_BASE_FUNC_NEG:
639
        return "__neg__"
640
    elif funcType == SOLLYA_BASE_FUNC_PI:
641
        return "pi"
642
    elif funcType == SOLLYA_BASE_FUNC_POW:
643
        return "^"
644
    elif funcType == SOLLYA_BASE_FUNC_PROCEDUREFUNCTION:
645
        return "procedureFunction"
646
    elif funcType == SOLLYA_BASE_FUNC_QUAD:
647
        return "quad"
648
    elif funcType == SOLLYA_BASE_FUNC_SIN:
649
        return "sin"
650
    elif funcType == SOLLYA_BASE_FUNC_SINGLE:
651
        return "single"
652
    elif funcType == SOLLYA_BASE_FUNC_SINH:
653
        return "sinh"
654
    elif funcType == SOLLYA_BASE_FUNC_SQRT:
655
        return "sqrt"
656
    elif funcType == SOLLYA_BASE_FUNC_SUB:
657
        return "-"
658
    elif funcType == SOLLYA_BASE_FUNC_TAN:
659
        return "tan"
660
    elif funcType == SOLLYA_BASE_FUNC_TANH:
661
        return "tanh"
662
    elif funcType == SOLLYA_BASE_FUNC_TRIPLEDOUBLE:
663
        return "tripleDouble"
664
    else:
665
        return None
666

    
667
def pobyso_get_constant(rnArgSa, constSo):
668
    """ Legacy function. See pobyso_get_constant_so_sa. """
669
    return pobyso_get_constant_so_sa(rnArgSa, constSo)
670
# End pobyso_get_constant
671

    
672
def pobyso_get_constant_so_sa(rnArgSa, constSo):
673
    """
674
    Set the value of rnArgSo to the value of constSo in MPFR_RNDN mode.
675
    rnArg must already exist and belong to some RealField.
676
    We assume that constSo points to a Sollya constant.
677
    """
678
    outcome = sollya_lib_get_constant(get_rn_value(rnArgSa), constSo)
679
    if outcome == 0: # Failure because constSo is not a constant expression.
680
        return None
681
    else:
682
        return outcome
683
# End  pobyso_get_constant_so_sa
684
   
685
def pobyso_get_constant_as_rn(ctExpSo):
686
    """ 
687
    Legacy function. See pobyso_get_constant_as_rn_so_sa. 
688
    """ 
689
    return(pobyso_get_constant_as_rn_so_sa(ctExpSo))
690
    
691
def pobyso_get_constant_as_rn_so_sa(constExpSo):
692
    """
693
    Get a Sollya constant as a Sage "real number".
694
    The precision of the floating-point number returned is that of the Sollya
695
    constant.
696
    """
697
    #print "Before computing precision of variable..."
698
    #pobyso_autoprint(constExpSo)
699
    precisionSa  = pobyso_get_prec_of_constant_so_sa(constExpSo)
700
    #print "precisionSa:", precisionSa
701
    ## If the expression can not be exactly converted, None is returned.
702
    #  In this case opt for the Sollya current expression.
703
    if precisionSa is None:
704
        precisionSa = pobyso_get_prec_so_sa()
705
    RRRR = RealField(precisionSa)
706
    rnSa = RRRR(0)
707
    outcome = sollya_lib_get_constant(get_rn_value(rnSa), constExpSo)
708
    if outcome == 0:
709
        return None
710
    else:
711
        return rnSa
712
# End pobyso_get_constant_as_rn_so_sa
713

    
714
def pobyso_get_constant_as_rn_with_rf(ctExp, realField):
715
    """ 
716
    Legacy function. See pobyso_get_constant_as_rn_with_rf_so_sa.
717
    """
718
    return pobyso_get_constant_as_rn_with_rf_so_sa(ctExp, realField)
719
# End pobyso_get_constant_as_rn_with_rf
720
    
721
def pobyso_get_constant_as_rn_with_rf_so_sa(ctExpSo, realFieldSa = None):
722
    """
723
    Get a Sollya constant as a Sage "real number".
724
    If no real field is specified, the precision of the floating-point number 
725
    returned is that of the Sollya constant.
726
    Otherwise is is that of the real field. Hence rounding may happen.
727
    """
728
    if realFieldSa is None:
729
        return pobyso_get_constant_as_rn_so_sa(ctExpSo)
730
    rnSa = realFieldSa(0)
731
    outcome = sollya_lib_get_constant(get_rn_value(rnSa), ctExpSo)
732
    if outcome == 0:
733
        return None
734
    else:
735
        return rnSa
736
# End pobyso_get_constant_as_rn_with_rf_so_sa
737

    
738
def pobyso_get_free_variable_name():
739
    """ 
740
    Legacy function. See pobyso_get_free_variable_name_so_sa.
741
    """
742
    return(pobyso_get_free_variable_name_so_sa())
743

    
744
def pobyso_get_free_variable_name_so_sa():
745
    return sollya_lib_get_free_variable_name()
746
    
747
def pobyso_get_function_arity(expressionSo):
748
    """ 
749
    Legacy function. See pobyso_get_function_arity_so_sa.
750
    """
751
    return(pobyso_get_function_arity_so_sa(expressionSo))
752

    
753
def pobyso_get_function_arity_so_sa(expressionSo):
754
    arity = c_int(0)
755
    sollya_lib_get_function_arity(byref(arity),expressionSo)
756
    return int(arity.value)
757

    
758
def pobyso_get_head_function(expressionSo):
759
    """ 
760
    Legacy function. See pobyso_get_head_function_so_sa. 
761
    """
762
    return(pobyso_get_head_function_so_sa(expressionSo)) 
763

    
764
def pobyso_get_head_function_so_sa(expressionSo):
765
    functionType = c_int(0)
766
    sollya_lib_get_head_function(byref(functionType), expressionSo)
767
    return int(functionType.value)
768

    
769
def pobyso_get_interval_from_range_so_sa(soRange, realIntervalFieldSa = None ):
770
    """
771
    Return the Sage interval corresponding to the Sollya range argument.
772
    If no reaIntervalField is passed as an argument, the interval bounds are not
773
    rounded: they are elements of RealIntervalField of the "right" precision
774
    to hold all the digits.
775
    """
776
    prec = c_int(0)
777
    if realIntervalFieldSa is None:
778
        retval = sollya_lib_get_prec_of_range(byref(prec), soRange, None)
779
        if retval == 0:
780
            return None
781
        realIntervalFieldSa = RealIntervalField(prec.value)
782
    intervalSa = realIntervalFieldSa(0,0)
783
    retval = \
784
        sollya_lib_get_interval_from_range(get_interval_value(intervalSa),\
785
                                           soRange)
786
    if retval == 0:
787
        return None
788
    return intervalSa
789
# End pobyso_get_interval_from_range_so_sa
790

    
791
def pobyso_get_list_elements(soObj):
792
    """ Legacy function. See pobyso_get_list_elements_so_so. """
793
    return pobyso_get_list_elements_so_so(soObj)
794
 
795
def pobyso_get_list_elements_so_so(objectListSo):
796
    """
797
    Get the Sollya list elements as a Sage/Python array of Sollya objects.
798
    
799
    INPUT:
800
    - objectListSo: a Sollya list of Sollya objects.
801
    
802
    OUTPUT:
803
    - a Sage/Python tuple made of:
804
      - a Sage/Python list of Sollya objects,
805
      - a Sage/Python int holding the number of elements,
806
      - a Sage/Python int stating (!= 0) that the list is end-elliptic.
807
    NOTE::
808
        We recover the addresses of the Sollya object from the list of pointers
809
        returned by sollya_lib_get_list_elements. The list itself is freed.
810
    TODO::
811
        Figure out what to do with numElements since the number of elements
812
        can easily be recovered from the list itself. 
813
        Ditto for isEndElliptic.
814
    """
815
    listAddress = POINTER(c_longlong)()
816
    numElements = c_int(0)
817
    isEndElliptic = c_int(0)
818
    listAsSageList = []
819
    result = sollya_lib_get_list_elements(byref(listAddress),\
820
                                          byref(numElements),\
821
                                          byref(isEndElliptic),\
822
                                          objectListSo)
823
    if result == 0 :
824
        return None
825
    for i in xrange(0, numElements.value, 1):
826
       #listAsSageList.append(sollya_lib_copy_obj(listAddress[i]))
827
       listAsSageList.append(listAddress[i])
828
       # Clear each of the elements returned by Sollya.
829
       #sollya_lib_clear_obj(listAddress[i])
830
    # Free the list itself.   
831
    sollya_lib_free(listAddress)
832
    return (listAsSageList, numElements.value, isEndElliptic.value)
833

    
834
def pobyso_get_max_prec_of_exp(soExp):
835
    """ Legacy function. See pobyso_get_max_prec_of_exp_so_sa. """
836
    return pobyso_get_max_prec_of_exp_so_sa(soExp)
837

    
838
def pobyso_get_max_prec_of_exp_so_sa(expSo):
839
    """
840
    Get the maximum precision used for the numbers in a Sollya expression.
841
    
842
    Arguments:
843
    soExp -- a Sollya expression pointer
844
    Return value:
845
    A Python integer
846
    TODO: 
847
    - error management;
848
    - correctly deal with numerical type such as DOUBLEEXTENDED.
849
    """
850
    if expSo is None:
851
        print inspect.stack()[0][3], ": expSo is None."
852
        return 0
853
    maxPrecision = 0
854
    minConstPrec = 0
855
    currentConstPrec = 0
856
    #pobyso_autoprint(expSo)
857
    operator = pobyso_get_head_function_so_sa(expSo)
858
    if (operator != SOLLYA_BASE_FUNC_CONSTANT) and \
859
    (operator != SOLLYA_BASE_FUNC_FREE_VARIABLE):
860
        (arity, subexpressions) = pobyso_get_subfunctions_so_sa(expSo)
861
        for i in xrange(arity):
862
            maxPrecisionCandidate = \
863
                pobyso_get_max_prec_of_exp_so_sa(subexpressions[i])
864
            if maxPrecisionCandidate > maxPrecision:
865
                maxPrecision = maxPrecisionCandidate
866
        return maxPrecision
867
    elif operator == SOLLYA_BASE_FUNC_CONSTANT:
868
        #minConstPrec = pobyso_get_min_prec_of_constant_so_sa(expSo)
869
        #currentConstPrec = pobyso_get_min_prec_of_constant_so_sa(soExp)
870
        #print minConstPrec, " - ", currentConstPrec 
871
        return pobyso_get_min_prec_of_constant_so_sa(expSo)
872
    
873
    elif operator == SOLLYA_BASE_FUNC_FREE_VARIABLE:
874
        return 0
875
    else:
876
        print "pobyso_get_max_prec_of_exp_so_sa: unexepected operator."
877
        return 0
878

    
879
def pobyso_get_min_prec_of_constant_so_sa(constExpSo):
880
    """
881
    Get the minimum precision necessary to represent the value of a Sollya
882
    constant.
883
    MPFR_MIN_PREC and powers of 2 are taken into account.
884
    We assume that constExpSo is a pointer to a Sollay constant expression.
885
    """
886
    constExpAsRnSa = pobyso_get_constant_as_rn_so_sa(constExpSo)
887
    return(min_mpfr_size(get_rn_value(constExpAsRnSa)))
888

    
889
def pobyso_get_poly_so_sa(polySo, realFieldSa=None):
890
    """
891
    Convert a Sollya polynomial into a Sage polynomial.
892
    Legacy function. Use pobyso_float_poly_so_sa() instead.
893
    """    
894
    return pobyso_float_poly_so_sa(polySo,realFieldSa)
895
# End pobyso_get_poly_so_sa
896

    
897
def pobyso_get_prec():
898
    """ Legacy function. See pobyso_get_prec_so_sa(). """
899
    return pobyso_get_prec_so_sa()
900

    
901
def pobyso_get_prec_so():
902
    """
903
    Get the current default precision in Sollya.
904
    The return value is a Sollya object.
905
    Usefull when modifying the precision back and forth by avoiding
906
    extra conversions.
907
    """
908
    return sollya_lib_get_prec(None)
909
    
910
def pobyso_get_prec_so_sa():
911
    """
912
    Get the current default precision in Sollya.
913
    The return value is Sage/Python int.
914
    """
915
    precSo = sollya_lib_get_prec()
916
    precSa = pobyso_constant_from_int_so_sa(precSo)
917
    sollya_lib_clear_obj(precSo)
918
    return precSa
919
# End pobyso_get_prec_so_sa.
920

    
921
def pobyso_get_prec_so_so_sa():
922
    """
923
    Return the current precision both as a Sollya object and a
924
    Sage integer as hybrid tuple.
925
    To avoid multiple calls for precision manipulations. 
926
    """
927
    precSo = sollya_lib_get_prec()
928
    precSa = pobyso_constant_from_int_so_sa(precSo)
929
    return (precSo, int(precSa))
930
    
931
def pobyso_get_prec_of_constant(ctExpSo):
932
    """ Legacy function. See pobyso_get_prec_of_constant_so_sa. """
933
    return pobyso_get_prec_of_constant_so_sa(ctExpSo)
934

    
935
def pobyso_get_prec_of_constant_so_sa(ctExpSo):
936
    """
937
    Tries to find a precision to represent ctExpSo without rounding.
938
    If not possible, returns None.
939
    """
940
    #print "Entering pobyso_get_prec_of_constant_so_sa..."
941
    prec = c_int(0)
942
    retc = sollya_lib_get_prec_of_constant(byref(prec), ctExpSo, None)
943
    if retc == 0:
944
        #print "pobyso_get_prec_of_constant_so_sa failed."
945
        return None
946
    #print "...exiting pobyso_get_prec_of_constant_so_sa."
947
    return int(prec.value)
948

    
949
def pobyso_get_prec_of_range_so_sa(rangeSo):
950
    """
951
    Returns the number of bits elements of a range are coded with.
952
    """
953
    prec = c_int(0)
954
    retc = sollya_lib_get_prec_of_range(byref(prec), rangeSo, None)
955
    if retc == 0:
956
        return(None)
957
    return int(prec.value)
958
# End pobyso_get_prec_of_range_so_sa()
959

    
960
def pobyso_get_sage_exp_from_sollya_exp(sollyaExpSo, realField = RR):
961
    """ Legacy function. See pobyso_get_sage_exp_from_sollya_exp_so_sa. """
962
    return pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, 
963
                                                     realField = RR)
964

    
965
def pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, realFieldSa = RR):
966
    """
967
    Get a Sage expression from a Sollya expression. 
968
    Currently only tested with polynomials with floating-point coefficients.
969
    Notice that, in the returned polynomial, the exponents are RealNumbers.
970
    """
971
    #pobyso_autoprint(sollyaExp)
972
    operatorSa = pobyso_get_head_function_so_sa(sollyaExpSo)
973
    sollyaLibFreeVariableName = sollya_lib_get_free_variable_name()
974
    ## Get rid of the "_"'s in "_x_", if any.
975
    sollyaLibFreeVariableName = re.sub('_', '', sollyaLibFreeVariableName)
976
    # Constants and the free variable are special cases.
977
    # All other operator are dealt with in the same way.
978
    if (operatorSa != SOLLYA_BASE_FUNC_CONSTANT) and \
979
       (operatorSa != SOLLYA_BASE_FUNC_FREE_VARIABLE):
980
        (aritySa, subexpressionsSa) = pobyso_get_subfunctions_so_sa(sollyaExpSo)
981
        if aritySa == 1:
982
            sageExpSa = eval(pobyso_function_type_as_string_so_sa(operatorSa) + \
983
            "(" + pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], \
984
            realFieldSa) + ")")
985
        elif aritySa == 2:
986
            # We do not get through the preprocessor.
987
            # The "^" operator is then a special case.
988
            if operatorSa == SOLLYA_BASE_FUNC_POW:
989
                operatorAsStringSa = "**"
990
            else:
991
                operatorAsStringSa = \
992
                    pobyso_function_type_as_string_so_sa(operatorSa)
993
            sageExpSa = \
994
              eval("pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], realFieldSa)"\
995
              + " " + operatorAsStringSa + " " + \
996
                   "pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[1], realFieldSa)")
997
        # We do not know yet how to deal with arity >= 3 
998
        # (is there any in Sollya anyway?).
999
        else:
1000
            sageExpSa = eval('None')
1001
        return sageExpSa
1002
    elif operatorSa == SOLLYA_BASE_FUNC_CONSTANT:
1003
        #print "This is a constant"
1004
        return pobyso_get_constant_as_rn_with_rf_so_sa(sollyaExpSo, realFieldSa)
1005
    elif operatorSa == SOLLYA_BASE_FUNC_FREE_VARIABLE:
1006
        #print "This is the free variable"
1007
        return eval(sollyaLibFreeVariableName)
1008
    else:
1009
        print "Unexpected"
1010
        return eval('None')
1011
# End pobyso_get_sage_exp_from_sollya_exp_so_sa
1012

    
1013

    
1014
def pobyso_get_subfunctions(expressionSo):
1015
    """ Legacy function. See pobyso_get_subfunctions_so_sa. """
1016
    return pobyso_get_subfunctions_so_sa(expressionSo) 
1017
# End pobyso_get_subfunctions.
1018
 
1019
def pobyso_get_subfunctions_so_sa(expressionSo):
1020
    """
1021
    Get the subfunctions of an expression.
1022
    Return the number of subfunctions and the list of subfunctions addresses.
1023
    S.T.: Could not figure out another way than that ugly list of declarations
1024
    to recover the addresses of the subfunctions.
1025
    We limit ourselves to arity 8 functions. 
1026
    """
1027
    subf0 = c_int(0)
1028
    subf1 = c_int(0)
1029
    subf2 = c_int(0)
1030
    subf3 = c_int(0)
1031
    subf4 = c_int(0)
1032
    subf5 = c_int(0)
1033
    subf6 = c_int(0)
1034
    subf7 = c_int(0)
1035
    subf8 = c_int(0)
1036
    arity = c_int(0)
1037
    nullPtr = POINTER(c_int)()
1038
    sollya_lib_get_subfunctions(expressionSo, byref(arity), \
1039
      byref(subf0), byref(subf1), byref(subf2), byref(subf3), \
1040
      byref(subf4), byref(subf5),\
1041
      byref(subf6), byref(subf7), byref(subf8), nullPtr, None) 
1042
#    byref(cast(subfunctions[0], POINTER(c_int))), \
1043
#    byref(cast(subfunctions[0], POINTER(c_int))), \
1044
#    byref(cast(subfunctions[2], POINTER(c_int))), \
1045
#    byref(cast(subfunctions[3], POINTER(c_int))), \
1046
#    byref(cast(subfunctions[4], POINTER(c_int))), \
1047
#    byref(cast(subfunctions[5], POINTER(c_int))), \
1048
#    byref(cast(subfunctions[6], POINTER(c_int))), \
1049
#    byref(cast(subfunctions[7], POINTER(c_int))), \
1050
#    byref(cast(subfunctions[8], POINTER(c_int))), nullPtr)
1051
    subfunctions = [subf0, subf1, subf2, subf3, subf4, subf5, subf6, subf7, \
1052
                    subf8]
1053
    subs = []
1054
    if arity.value > pobyso_max_arity:
1055
        return(0,[])
1056
    for i in xrange(arity.value):
1057
        subs.append(int(subfunctions[i].value))
1058
        #print subs[i]
1059
    return (int(arity.value), subs)
1060
# End pobyso_get_subfunctions_so_sa
1061
    
1062
def pobyso_guess_degree_sa_sa(functionSa, intervalSa, approxErrorSa, 
1063
                              weightSa=None, degreeBoundSa=None):
1064
    """
1065
    Sa_sa variant of the solly_guessdegree function.
1066
    Return 0 if something goes wrong.
1067
    """
1068
    functionAsStringSa = functionSa._assume_str().replace('_SAGE_VAR_', '')
1069
    functionSo = pobyso_parse_string_sa_so(functionAsStringSa)
1070
    if pobyso_is_error_so_sa(functionSo):
1071
        sollya_lib_clear_obj(functionSo)
1072
        return 0
1073
    rangeSo = pobyso_interval_to_range_sa_so(intervalSa)
1074
    # The approximation error is expected to be a floating point number.
1075
    if pobyso_is_floating_point_number_sa_sa(approxErrorSa):
1076
        approxErrorSo = pobyso_constant_sa_so(approxErrorSa)
1077
    else:
1078
        approxErrorSo = pobyso_constant_sa_so(RR(approxErrorSa))
1079
    if not weightSa is None:
1080
        weightAsStringSa = weightSa._assume_str().replace('_SAGE_VAR_', '')
1081
        weightSo = pobyso_parse_string_sa_so(weightAsStringSa)
1082
        if pobyso_is_error_so_sa(weightSo):
1083
            sollya_lib_clear_obj(functionSo)
1084
            sollya_lib_clear_obj(rangeSo)
1085
            sollya_lib_clear_obj(approxErrorSo)   
1086
            sollya_lib_clear_obj(weightSo)
1087
            return 0   
1088
    else:
1089
        weightSo = None
1090
    if not degreeBoundSa is None:
1091
        degreeBoundSo = pobyso_constant_from_int_sa_so(degreeBoundSa)
1092
    else:
1093
        degreeBoundSo = None
1094
    guessedDegreeSa = pobyso_guess_degree_so_sa(functionSo,
1095
                                                rangeSo,
1096
                                                approxErrorSo,
1097
                                                weightSo,
1098
                                                degreeBoundSo)
1099
    sollya_lib_clear_obj(functionSo)
1100
    sollya_lib_clear_obj(rangeSo)
1101
    sollya_lib_clear_obj(approxErrorSo)
1102
    if not weightSo is None:
1103
        sollya_lib_clear_obj(weightSo)
1104
    if not degreeBoundSo is None:
1105
        sollya_lib_clear_obj(degreeBoundSo)
1106
    return guessedDegreeSa
1107
# End poyso_guess_degree_sa_sa
1108

    
1109
def pobyso_guess_degree_so_sa(functionSo, rangeSo, errorSo, weightSo=None, \
1110
                              degreeBoundSo=None):
1111
    """
1112
    Thin wrapper around the guessdegree function.
1113
    Nevertheless, some precision control stuff has been appended.
1114
    """
1115
    # Deal with Sollya internal precision issues: if it is too small, 
1116
    # compared with the error, increases it to about twice -log2(error).
1117
    errorSa = pobyso_get_constant_as_rn_with_rf_so_sa(errorSo)
1118
    log2ErrorSa = errorSa.log2()
1119
    if log2ErrorSa < 0:
1120
        neededPrecisionSa = int(2 * int(-log2ErrorSa) / 64) * 64
1121
    else:
1122
        neededPrecisionSa = int(2 * int(log2ErrorSa) / 64) * 64
1123
    #print "Needed precision:", neededPrecisionSa
1124
    sollyaPrecisionHasChanged = False
1125
    (initialPrecSo, initialPrecSa) = pobyso_get_prec_so_so_sa()
1126
    if neededPrecisionSa > initialPrecSa:
1127
        if neededPrecisionSa <= 2:
1128
            print inspect.stack()[0][3], ": precision change <= 2 requested."
1129
        pobyso_set_prec_sa_so(neededPrecisionSa)
1130
        sollyaPrecisionHasChanged = True
1131
    #print "Guessing degree..."
1132
    # weightSo and degreeBoundsSo are optional arguments.
1133
    # As declared, sollya_lib_guessdegree must take 5 arguments.
1134
    if weightSo is None:
1135
        degreeRangeSo = sollya_lib_guessdegree(functionSo, rangeSo, errorSo, 
1136
                                               0, 0, None)
1137
    elif degreeBoundSo is None:
1138
        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, \
1139
                                                errorSo, weightSo, 0, None)
1140
    else:
1141
        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, errorSo, \
1142
                                                weightSo, degreeBoundSo, None)
1143
    #print "...degree guess done."
1144
    # Restore internal precision, if applicable.
1145
    if sollyaPrecisionHasChanged:
1146
        pobyso_set_prec_so_so(initialPrecSo)
1147
        sollya_lib_clear_obj(initialPrecSo)
1148
    degreeIntervalSa = pobyso_range_to_interval_so_sa(degreeRangeSo)
1149
    sollya_lib_clear_obj(degreeRangeSo)
1150
    # When ok, both bounds match.
1151
    # When the degree bound is too low, the upper bound is the degree
1152
    # for which the error can be honored.
1153
    # When it really goes wrong, the upper bound is infinity. 
1154
    if degreeIntervalSa.lower() == degreeIntervalSa.upper():
1155
        return int(degreeIntervalSa.lower())
1156
    else:
1157
        if degreeIntervalSa.upper().is_infinity():
1158
            return None
1159
        else:
1160
            return int(degreeIntervalSa.upper())
1161
    # End pobyso_guess_degree_so_sa
1162

    
1163
def pobyso_inf_so_so(intervalSo):
1164
    """
1165
    Very thin wrapper around sollya_lib_inf().
1166
    """
1167
    return sollya_lib_inf(intervalSo)
1168
# End pobyso_inf_so_so.
1169
    
1170
def pobyso_infnorm_so_so(func, interval, file = None, intervalList = None):
1171
    print "Do not use this function. User pobyso_supnorm_so_so instead."
1172
    return None
1173

    
1174
def pobyso_interval_to_range_sa_so(intervalSa, precisionSa=None):
1175
    if precisionSa is None:
1176
        precisionSa = intervalSa.parent().precision()
1177
    intervalSo = pobyso_bounds_to_range_sa_so(intervalSa.lower(),\
1178
                                              intervalSa.upper(),\
1179
                                              precisionSa)
1180
    return intervalSo
1181
# End pobyso_interval_to_range_sa_so
1182

    
1183
def pobyso_is_error_so_sa(objSo):
1184
    """
1185
    Thin wrapper around the sollya_lib_obj_is_error() function.
1186
    """
1187
    if sollya_lib_obj_is_error(objSo) != 0:
1188
        return True
1189
    else:
1190
        return False
1191
# End pobyso_is_error-so_sa
1192

    
1193
def pobyso_is_floating_point_number_sa_sa(numberSa):
1194
    """
1195
    Check whether a Sage number is floating point.
1196
    Exception stuff added because numbers other than
1197
    floating-point ones do not have the is_real() attribute.
1198
    """
1199
    try:
1200
        return numberSa.is_real()
1201
    except AttributeError:
1202
        return False
1203
# End pobyso_is_floating_piont_number_sa_sa
1204

    
1205
def pobyso_is_range_so_sa(rangeCandidateSo):
1206
    """
1207
    Thin wrapper over sollya_lib_is_range.
1208
    """
1209
    return sollya_lib_obj_is_range(rangeCandidateSo) != 0
1210
    
1211
# End pobyso_is_range_so_sa
1212

    
1213

    
1214
def pobyso_lib_init():
1215
    sollya_lib_init(None)
1216

    
1217
def pobyso_lib_close():
1218
    sollya_lib_close(None)
1219

    
1220
def pobyso_name_free_variable(freeVariableNameSa):
1221
    """ Legacy function. See pobyso_name_free_variable_sa_so. """
1222
    pobyso_name_free_variable_sa_so(freeVariableNameSa)
1223

    
1224
def pobyso_name_free_variable_sa_so(freeVariableNameSa):
1225
    """
1226
    Set the free variable name in Sollya from a Sage string.
1227
    """
1228
    sollya_lib_name_free_variable(freeVariableNameSa)
1229

    
1230
def pobyso_parse_string(string):
1231
    """ Legacy function. See pobyso_parse_string_sa_so. """
1232
    return pobyso_parse_string_sa_so(string)
1233
 
1234
def pobyso_parse_string_sa_so(string):
1235
    """
1236
    Get the Sollya expression computed from a Sage string or
1237
    a Sollya error object if parsing failed.
1238
    """
1239
    return sollya_lib_parse_string(string)
1240

    
1241
def pobyso_precision_so_sa(ctExpSo):
1242
    """
1243
    Computes the necessary precision to represent a number.
1244
    If x is not zero, it can be uniquely written as x = m · 2e 
1245
    where m is an odd integer and e is an integer. 
1246
    precision(x) returns the number of bits necessary to write m 
1247
    in binary (i.e. ceil(log2(m))).
1248
    """
1249
    #TODO: take care of the special case: 0, @NaN@, @Inf@
1250
    precisionSo = sollya_lib_precision(ctExpSo)
1251
    precisionSa = pobyso_constant_from_int_so_sa(precisionSo)
1252
    sollya_lib_clear_obj(precisionSo)
1253
    return precisionSa
1254
# End pobyso_precision_so_sa
1255

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

    
1422
def pobyso_polynomial_degree_so_sa(polySo):
1423
    """
1424
    Return the degree of a Sollya polynomial as a Sage int.
1425
    """
1426
    degreeSo = sollya_lib_degree(polySo)
1427
    return pobyso_constant_from_int_so_sa(degreeSo)
1428
# End pobyso_polynomial_degree_so_sa
1429

    
1430
def pobyso_polynomial_degree_so_so(polySo):
1431
    """
1432
    Thin wrapper around lib_sollya_degree().
1433
    """
1434
    return sollya_lib_degree(polySo)
1435
# End pobyso_polynomial_degree_so_so
1436
                                                                  
1437
def pobyso_range(rnLowerBound, rnUpperBound):
1438
    """ Legacy function. See pobyso_range_sa_so. """
1439
    return pobyso_range_sa_so(rnLowerBound, rnUpperBound) 
1440

    
1441
def pobyso_range_from_bounds_sa_so(rnLowerBound, rnUpperBound, precSa = None):
1442
    """
1443
    Create a Sollya range from 2 Sage real numbers as bounds
1444
    """
1445
    # TODO check precision stuff.
1446
    sollyaPrecChanged = False
1447
    (initialSollyaPrecSo, initialSollyaPrecSa) = \
1448
        pobyso_get_prec_so_so_sa()
1449
    if precSa is None:
1450
        precSa = max(rnLowerBound.parent().prec(), rnUpperBound.parent().prec())
1451
    if precSa > initialSollyaPrecSa:
1452
        if precSa <= 2:
1453
            print inspect.stack()[0][3], ": precision change <= 2 requested."
1454
        pobyso_set_prec_sa_so(precSa)
1455
        sollyaPrecChanged = True
1456
    rangeSo = sollya_lib_range_from_bounds(get_rn_value(rnLowerBound),
1457
                                           get_rn_value(rnUpperBound))
1458
    if sollyaPrecChanged:
1459
        pobyso_set_prec_so_so(initialSollyaPrecSo)
1460
        sollya_lib_clear_obj(initialSollyaPrecSo)
1461
    return rangeSo
1462
# End pobyso_range_from_bounds_sa_so
1463

    
1464
def pobyso_range_max_abs_so_so(rangeSo):
1465
    """
1466
    Return, as Sollya constant, the maximum absolute value of a Sollay range.
1467
    """
1468
    lowerBoundSo = sollya_lib_inf(rangeSo)
1469
    upperBoundSo = sollya_lib_sup(rangeSo)
1470
    #
1471
    lowerBoundSo = sollya_lib_build_function_abs(lowerBoundSo)
1472
    upperBoundSo = sollya_lib_build_function_abs(upperBoundSo)
1473
    #pobyso_autoprint(lowerBoundSo)
1474
    #pobyso_autoprint(upperBoundSo)
1475
    #
1476
    maxAbsSo = sollya_lib_max(lowerBoundSo, upperBoundSo, None)
1477
    #sollya_lib_clear_obj(lowerBoundSo)
1478
    #sollya_lib_clear_obj(upperBoundSo)
1479
    return maxAbsSo
1480
# End pobyso_range_max_abs_so_so
1481
 
1482
def pobyso_range_to_interval_so_sa(rangeSo, realIntervalFieldSa = None):
1483
    """
1484
    Get a Sage interval from a Sollya range.
1485
    If no realIntervalField is given as a parameter, the Sage interval
1486
    precision is that of the Sollya range.
1487
    Otherwise, the precision is that of the realIntervalField. In this case
1488
    rounding may happen.
1489
    """
1490
    if realIntervalFieldSa is None:
1491
        precSa = pobyso_get_prec_of_range_so_sa(rangeSo)
1492
        realIntervalFieldSa = RealIntervalField(precSa)
1493
    intervalSa = \
1494
        pobyso_get_interval_from_range_so_sa(rangeSo, realIntervalFieldSa) 
1495
    return intervalSa
1496
# End pobyso_range_to_interval_so_sa
1497
#
1498
def pobyso_relative_so_so():
1499
    """
1500
    Very thin wrapper around the sollya_lib_relative function.
1501
    """
1502
    return sollya_lib_relative()
1503
# End pobyso_relative_so_so
1504
#
1505
def pobyso_rat_poly_sa_so(polySa, precSa = None):
1506
    """
1507
    Create a Sollya polynomial from a Sage rational polynomial.
1508
    """
1509
    ## TODO: filter arguments.
1510
    ## Precision. If no precision is given, use the current precision
1511
    #  of Sollya.
1512
    if precSa is None:
1513
        precSa =  pobyso_get_prec_so_sa()
1514
    #print "Precision:",  precSa
1515
    RRR = RealField(precSa)
1516
    ## Create a Sage polynomial in the "right" precision.
1517
    P_RRR = RRR[polySa.variables()[0]]
1518
    polyFloatSa = P_RRR(polySa)
1519
    ## Make sure no precision is provided: pobyso_float_poly_sa_so will
1520
    #  recover it all by itself and not make an extra conversion.
1521
    return pobyso_float_poly_sa_so(polyFloatSa)
1522
    
1523
# End pobyso_rat_poly_sa_so    
1524
    
1525
def pobyso_remez_canonical_sa_sa(func, \
1526
                                 degree, \
1527
                                 lowerBound, \
1528
                                 upperBound, \
1529
                                 weight = None, \
1530
                                 quality = None):
1531
    """
1532
    All arguments are Sage/Python.
1533
    The functions (func and weight) must be passed as expressions or strings.
1534
    Otherwise the function fails. 
1535
    The return value is a Sage polynomial.
1536
    """
1537
    var('zorglub')    # Dummy variable name for type check only. Type of 
1538
    # zorglub is "symbolic expression".
1539
    polySo = pobyso_remez_canonical_sa_so(func, \
1540
                                 degree, \
1541
                                 lowerBound, \
1542
                                 upperBound, \
1543
                                 weight, \
1544
                                 quality)
1545
    # String test
1546
    if parent(func) == parent("string"):
1547
        functionSa = eval(func)
1548
    # Expression test.
1549
    elif type(func) == type(zorglub):
1550
        functionSa = func
1551
    else:
1552
        return None
1553
    #
1554
    maxPrecision = 0
1555
    if polySo is None:
1556
        return(None)
1557
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
1558
    RRRRSa = RealField(maxPrecision)
1559
    polynomialRingSa = RRRRSa[functionSa.variables()[0]]
1560
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, RRRRSa)
1561
    polySa = polynomial(expSa, polynomialRingSa)
1562
    sollya_lib_clear_obj(polySo)
1563
    return(polySa)
1564
# End pobyso_remez_canonical_sa_sa
1565
    
1566
def pobyso_remez_canonical(func, \
1567
                           degree, \
1568
                           lowerBound, \
1569
                           upperBound, \
1570
                           weight = "1", \
1571
                           quality = None):
1572
    """ Legacy function. See pobyso_remez_canonical_sa_so. """
1573
    return(pobyso_remez_canonical_sa_so(func, \
1574
                                        degree, \
1575
                                        lowerBound, \
1576
                                        upperBound, \
1577
                                        weight, \
1578
                                        quality))
1579
# End pobyso_remez_canonical.
1580

    
1581
def pobyso_remez_canonical_sa_so(func, \
1582
                                 degree, \
1583
                                 lowerBound, \
1584
                                 upperBound, \
1585
                                 weight = None, \
1586
                                 quality = None):
1587
    """
1588
    All arguments are Sage/Python.
1589
    The functions (func and weight) must be passed as expressions or strings.
1590
    Otherwise the function fails. 
1591
    The return value is a pointer to a Sollya function.
1592
    lowerBound and upperBound mus be reals.
1593
    """
1594
    var('zorglub')    # Dummy variable name for type check only. Type of
1595
    # zorglub is "symbolic expression".
1596
    currentVariableNameSa = None
1597
    # The func argument can be of different types (string, 
1598
    # symbolic expression...)
1599
    if parent(func) == parent("string"):
1600
        localFuncSa = sage_eval(func,globals())
1601
        if len(localFuncSa.variables()) > 0:
1602
            currentVariableNameSa = localFuncSa.variables()[0]
1603
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1604
            functionSo = \
1605
              sollya_lib_parse_string(localFuncSa._assume_str().replace('_SAGE_VAR_', ''))
1606
    # Expression test.
1607
    elif type(func) == type(zorglub):
1608
        # Until we are able to translate Sage expressions into Sollya 
1609
        # expressions : parse the string version.
1610
        if len(func.variables()) > 0:
1611
            currentVariableNameSa = func.variables()[0]
1612
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1613
            functionSo = \
1614
              sollya_lib_parse_string(func._assume_str().replace('_SAGE_VAR_', ''))
1615
    else:
1616
        return(None)
1617
    if weight is None: # No weight given -> 1.
1618
        weightSo = pobyso_constant_1_sa_so()
1619
    elif parent(weight) == parent("string"): # Weight given as string: parse it.
1620
        weightSo = sollya_lib_parse_string(func)
1621
    elif type(weight) == type(zorglub): # Weight given as symbolice expression.
1622
        functionSo = \
1623
          sollya_lib_parse_string_sa_so(weight._assume_str().replace('_SAGE_VAR_', ''))
1624
    else:
1625
        return(None)
1626
    degreeSo = pobyso_constant_from_int(degree)
1627
    rangeSo = pobyso_bounds_to_range_sa_so(lowerBound, upperBound)
1628
    if not quality is None:
1629
        qualitySo= pobyso_constant_sa_so(quality)
1630
    else:
1631
        qualitySo = None
1632
        
1633
    remezPolySo = sollya_lib_remez(functionSo, \
1634
                                   degreeSo, \
1635
                                   rangeSo, \
1636
                                   weightSo, \
1637
                                   qualitySo, \
1638
                                   None)
1639
    sollya_lib_clear_obj(functionSo)
1640
    sollya_lib_clear_obj(degreeSo)
1641
    sollya_lib_clear_obj(rangeSo)
1642
    sollya_lib_clear_obj(weightSo)
1643
    if not qualitySo is None:
1644
        sollya_lib_clear_obj(qualitySo)
1645
    return(remezPolySo)
1646
# End pobyso_remez_canonical_sa_so
1647

    
1648
def pobyso_remez_canonical_so_so(funcSo, \
1649
                                 degreeSo, \
1650
                                 rangeSo, \
1651
                                 weightSo = pobyso_constant_1_sa_so(),\
1652
                                 qualitySo = None):
1653
    """
1654
    All arguments are pointers to Sollya objects.
1655
    The return value is a pointer to a Sollya function.
1656
    """
1657
    if not sollya_lib_obj_is_function(funcSo):
1658
        return(None)
1659
    return(sollya_lib_remez(funcSo, degreeSo, rangeSo, weightSo, qualitySo, None))
1660
# End pobyso_remez_canonical_so_so.
1661

    
1662
def pobyso_remez_exponents_list_sa_so(func,     \
1663
                                 exponentsList, \
1664
                                 lowerBound,    \
1665
                                 upperBound,    \
1666
                                 weight = None, \
1667
                                 quality = None):
1668
    """
1669
    All arguments are Sage/Python.
1670
    The functions (func and weight) must be passed as expressions or strings.
1671
    Otherwise the function fails. 
1672
    The return value is a pointer to a Sollya function.
1673
    lowerBound and upperBound mus be reals.
1674
    """
1675
    var('zorglub')    # Dummy variable name for type check only. Type of
1676
    # zorglub is "symbolic expression".
1677
    currentVariableNameSa = None
1678
    # The func argument can be of different types (string, 
1679
    # symbolic expression...)
1680
    if parent(func) == parent("string"):
1681
        localFuncSa = sage_eval(func,globals())
1682
        if len(localFuncSa.variables()) > 0:
1683
            currentVariableNameSa = localFuncSa.variables()[0]
1684
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1685
            functionSo = \
1686
              sollya_lib_parse_string(localFuncSa._assume_str().replace('_SAGE_VAR_', ''))
1687
    # Expression test.
1688
    elif type(func) == type(zorglub):
1689
        # Until we are able to translate Sage expressions into Sollya 
1690
        # expressions : parse the string version.
1691
        if len(func.variables()) > 0:
1692
            currentVariableNameSa = func.variables()[0]
1693
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1694
            functionSo = \
1695
              sollya_lib_parse_string(func._assume_str().replace('_SAGE_VAR_', ''))
1696
    else:
1697
        return(None)
1698
    ## Deal with the weight, much in the same way as with the function.
1699
    if weight is None: # No weight given -> 1.
1700
        weightSo = pobyso_constant_1_sa_so()
1701
    elif parent(weight) == parent("string"): # Weight given as string: parse it.
1702
        weightSo = sollya_lib_parse_string(func)
1703
    elif type(weight) == type(zorglub): # Weight given as symbolice expression.
1704
        functionSo = \
1705
          sollya_lib_parse_string(weight._assume_str().replace('_SAGE_VAR_', '',100))
1706
    else:
1707
        return(None)
1708
    rangeSo = pobyso_bounds_to_range_sa_so(lowerBound, upperBound)
1709
    if not quality is None:
1710
        qualitySo= pobyso_constant_sa_so(quality)
1711
    else:
1712
        qualitySo = None
1713
    #
1714
    ## Tranform the Sage list of exponents into a Sollya list.
1715
    exponentsListSo = pobyso_build_list_of_ints_sa_so(*exponentsList)
1716
    remezPolySo = sollya_lib_remez(functionSo, \
1717
                                   exponentsListSo, \
1718
                                   rangeSo, \
1719
                                   weightSo, \
1720
                                   qualitySo, \
1721
                                   None)
1722
    sollya_lib_clear_obj(functionSo)
1723
    sollya_lib_clear_obj(exponentsListSo)
1724
    sollya_lib_clear_obj(rangeSo)
1725
    sollya_lib_clear_obj(weightSo)
1726
    if not qualitySo is None:
1727
        sollya_lib_clear_obj(qualitySo)
1728
    return(remezPolySo)
1729
# End pobyso_remez_exponentsList_sa_so
1730
#
1731
def pobyso_round_coefficients_so_so(polySo, truncFormatListSo):
1732
    """
1733
    A wrapper around the "classical" sollya_lib_roundcoefficients: a Sollya
1734
    polynomial and a Sollya list are given as arguments.
1735
    """
1736
    return sollya_lib_roundcoefficients(polySo, truncFormatListSo)
1737

    
1738
def pobyso_round_coefficients_progressive_so_so(polySo, 
1739
                                                funcSo,
1740
                                                precSo,
1741
                                                intervalSo,
1742
                                                icSo,
1743
                                                currentApproxErrorSo,
1744
                                                approxAccurSo,
1745
                                                debug=False):
1746
    if debug:
1747
        print "Input arguments:"
1748
        print "Polynomial: ", ; pobyso_autoprint(polySo)
1749
        print "Function: ", ; pobyso_autoprint(funcSo)
1750
        print "Internal precision: ", ; pobyso_autoprint(precSo)
1751
        print "Interval: ", ; pobyso_autoprint(intervalSo)
1752
        print "Current approximation error: ", ; pobyso_autoprint(currentApproxErrorSo)
1753
        print "Requested approxiation error: ", ; pobyso_autoprint(approxAccurSo)
1754
        print "________________"
1755
    approxAccurSa        = pobyso_get_constant_as_rn_so_sa(approxAccurSo)
1756
    currentApproxErrorSa = pobyso_get_constant_as_rn_so_sa(currentApproxErrorSo)
1757
    ## If the current approximation error is too close to the target, there is
1758
    #  no possible gain.
1759
    if currentApproxErrorSa >= approxAccurSa / 2:
1760
        #### Do not return the initial argument but copies: caller may free 
1761
        #    as inutile after call.
1762
        return (sollya_lib_copy_obj(polySo),
1763
                sollya_lib_copy_obj(currentApproxErrorSo))
1764
    degreeSa             = pobyso_polynomial_degree_so_sa(polySo)
1765
    intervalSa           = pobyso_range_to_interval_so_sa(intervalSo)
1766
    
1767
    if debug:
1768
        print "degreeSa              :", degreeSa
1769
        print "intervalSa            :", intervalSa.str(style='brackets')
1770
        print "currentApproxErrorSa  :", currentApproxErrorSa 
1771
        print "approxAccurSa         :", approxAccurSa 
1772
    ### Start with a 0 value expression.
1773
    radiusSa = intervalSa.absolute_diameter() / 2
1774
    if debug:
1775
        print "log2(radius):", RR(radiusSa).log2()
1776
    iterIndex = 0
1777
    while True: 
1778
        resPolySo = pobyso_constant_0_sa_so()
1779
        roundedPolyApproxAccurSa = approxAccurSa / 2
1780
        currentRadiusPowerSa = 1 
1781
        for degree in xrange(0,degreeSa + 1):
1782
            #### At round 0, use the agressive formula. At round 1, run the
1783
            #    proved formula.
1784
            if iterIndex == 0:
1785
                roundingPowerSa = \
1786
                    floor(((currentRadiusPowerSa/roundedPolyApproxAccurSa)*(degree+1)).log2())
1787
            else:
1788
                roundingPowerSa = \
1789
                    floor(((currentRadiusPowerSa/roundedPolyApproxAccurSa)*(degreeSa+1)).log2())
1790
            ## Under extreme conditions the above formulas can evaluate under 2, which is the
1791
            #  minimal precision of an MPFR number.
1792
            if roundingPowerSa < 2:
1793
                roundingPowerSa = 2
1794
            if debug:
1795
                print "roundedPolyApproxAccurSa", roundedPolyApproxAccurSa
1796
                print "currentRadiusPowerSa", currentRadiusPowerSa
1797
                print "Current rounding exponent:", roundingPowerSa
1798
            currentRadiusPowerSa *= radiusSa
1799
            index1So = pobyso_constant_from_int_sa_so(degree)
1800
            index2So = pobyso_constant_from_int_sa_so(degree)
1801
            ### Create a monomial with:
1802
            #   - the coefficient in the initial monomial at the current degrree;
1803
            #   - the current exponent;
1804
            #   - the free variable.
1805
            cmonSo  = \
1806
                sollya_lib_build_function_mul(sollya_lib_coeff(polySo, index1So),
1807
                                              sollya_lib_build_function_pow( \
1808
                                                sollya_lib_build_function_free_variable(), \
1809
                                                index2So))
1810
            roundingPowerSo = pobyso_constant_from_int_sa_so(roundingPowerSa)
1811
            cmonrSo = pobyso_round_coefficients_single_so_so(cmonSo, roundingPowerSo)
1812
            sollya_lib_clear_obj(cmonSo)
1813
            ### Add to the result polynomial.
1814
            resPolySo = sollya_lib_build_function_add(resPolySo,
1815
                                                      cmonrSo)
1816
        # End for.
1817
        ### Check the new polynomial.
1818
        freeVarSo     = sollya_lib_build_function_free_variable()
1819
        changeVarSo   = sollya_lib_sub(freeVarSo, icSo)
1820
        resPolyCvSo   = sollya_lib_evaluate(resPolySo, changeVarSo)
1821
        errFuncSo     = sollya_lib_build_function_sub(sollya_lib_copy_obj(funcSo), 
1822
                                                      resPolyCvSo)
1823
        infNormSo = sollya_lib_dirtyinfnorm(errFuncSo, intervalSo)
1824
        ### This also clears resPolyCvSo.
1825
        sollya_lib_clear_obj(errFuncSo)
1826
        cerrSa    = pobyso_get_constant_as_rn_so_sa(infNormSo)
1827
        if debug:
1828
            print "Error of the new polynomial:", cerrSa
1829
        ### If at round 1, return the initial polynomial error. This should
1830
        #   never happen since the rounding algorithm is proved. But some 
1831
        #   circumstances may break it (e.g. internal precision of tools).
1832
        if cerrSa > approxAccurSa:
1833
            if iterIndex > 0: # Round 1 and beyond.
1834
                sollya_lib_clear_obj(resPolySo)
1835
                sollya_lib_clear_obj(infNormSo)
1836
                #### Do not return the arguments but copies: the caller may free
1837
                #    free the former as inutile after call.
1838
                return (sollya_lib_copy_obj(polySo), 
1839
                        sollya_lib_copy_obj(currentApproxErrorSo))
1840
            else: # Round 0, got round 1
1841
                sollya_lib_clear_obj(resPolySo)
1842
                sollya_lib_clear_obj(infNormSo)
1843
                iterIndex += 1
1844
                continue
1845
        ### If get here it is because cerrSa <= approxAccurSa
1846
        ### Approximation error of the new polynomial is acceptable.
1847
        return (resPolySo, infNormSo)
1848
    # End while True
1849
# End pobyso_round_coefficients_progressive_so_so
1850
    
1851
def pobyso_round_coefficients_single_so_so(polySo, commonPrecSo):
1852
    """
1853
    Create a rounded coefficients polynomial from polynomial argument to
1854
    the number of bits in size argument.
1855
    All coefficients are set to the same precision.
1856
    """
1857
    ## TODO: check arguments.
1858
    endEllipListSo = pobyso_build_end_elliptic_list_so_so(commonPrecSo)
1859
    polySo = sollya_lib_roundcoefficients(polySo, endEllipListSo, None)
1860
    sollya_lib_clear_obj(endEllipListSo)
1861
    #sollya_lib_clear_obj(endEllipListSo)
1862
    return polySo
1863
    
1864
# End pobyso_round_coefficients_single_so_so
1865

    
1866
def pobyso_set_canonical_off():
1867
    sollya_lib_set_canonical(sollya_lib_off())
1868

    
1869
def pobyso_set_canonical_on():
1870
    sollya_lib_set_canonical(sollya_lib_on())
1871

    
1872
def pobyso_set_prec(p):
1873
    """ Legacy function. See pobyso_set_prec_sa_so. """
1874
    pobyso_set_prec_sa_so(p)
1875

    
1876
def pobyso_set_prec_sa_so(p):
1877
    #a = c_int(p)
1878
    #precSo = c_void_p(sollya_lib_constant_from_int(a))
1879
    #precSo = sollya_lib_constant_from_int(a)
1880
    precSo = pobyso_constant_from_int_sa_so(p)
1881
    sollya_lib_set_prec(precSo)
1882
    sollya_lib_clear_obj(precSo)
1883
# End pobyso_set_prec_sa_so.
1884

    
1885
def pobyso_set_prec_so_so(newPrecSo):
1886
    sollya_lib_set_prec(newPrecSo)
1887
# End pobyso_set_prec_so_so.
1888

    
1889
def pobyso_inf_so_so(intervalSo):
1890
    """
1891
    Very thin wrapper around sollya_lib_inf().
1892
    """
1893
    return sollya_lib_inf(intervalSo)
1894
# End pobyso_inf_so_so.
1895
    
1896
def pobyso_supnorm_so_so(polySo, funcSo, intervalSo, errorTypeSo = None,\
1897
                         accuracySo = None):
1898
    """
1899
    Computes the supnorm of the approximation error between the given 
1900
    polynomial and function.
1901
    errorTypeSo defaults to "absolute".
1902
    accuracySo defaults to 2^(-40).
1903
    """
1904
    if errorTypeSo is None:
1905
        errorTypeSo = sollya_lib_absolute(None)
1906
        errorTypeIsNone = True
1907
    else:
1908
        errorTypeIsNone = False
1909
    #
1910
    if accuracySo is None:
1911
        # Notice the **: we are in Pythonland here!
1912
        accuracySo = pobyso_constant_sa_so(RR(2**(-40)))
1913
        accuracyIsNone = True
1914
    else:
1915
        accuracyIsNone = False
1916
    #pobyso_autoprint(accuracySo)
1917
    resultSo = \
1918
        sollya_lib_supnorm(polySo, funcSo, intervalSo, errorTypeSo, \
1919
                              accuracySo)
1920
    if errorTypeIsNone:
1921
        sollya_lib_clear_obj(errorTypeSo)
1922
    if accuracyIsNone:
1923
        sollya_lib_clear_obj(accuracySo)
1924
    return resultSo
1925
# End pobyso_supnorm_so_so
1926

    
1927
def pobyso_taylor_expansion_no_change_var_so_so(functionSo, 
1928
                                                degreeSo, 
1929
                                                rangeSo,
1930
                                                errorTypeSo=None, 
1931
                                                sollyaPrecSo=None):
1932
    """
1933
    Compute the Taylor expansion without the variable change
1934
    x -> x-intervalCenter.
1935
    """
1936
    # Change internal Sollya precision, if needed.
1937
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
1938
    sollyaPrecChanged = False
1939
    if sollyaPrecSo is None:
1940
        pass
1941
    else:
1942
        sollya_lib_set_prec(sollyaPrecSo)
1943
        sollyaPrecChanged = True
1944
    # Error type stuff: default to absolute.
1945
    if errorTypeSo is None:
1946
        errorTypeIsNone = True
1947
        errorTypeSo = sollya_lib_absolute(None)
1948
    else:
1949
        errorTypeIsNone = False
1950
    intervalCenterSo = sollya_lib_mid(rangeSo, None)
1951
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo,
1952
                                         intervalCenterSo,
1953
                                         rangeSo, errorTypeSo, None)
1954
    # taylorFormListSaSo is a Python list of Sollya objects references that 
1955
    # are copies of the elements of taylorFormSo.
1956
    # pobyso_get_list_elements_so_so clears taylorFormSo.
1957
    (taylorFormListSaSo, numElementsSa, isEndEllipticSa) = \
1958
        pobyso_get_list_elements_so_so(taylorFormSo)
1959
    polySo = sollya_lib_copy_obj(taylorFormListSaSo[0])
1960
    #print "Num elements:", numElementsSa
1961
    sollya_lib_clear_obj(taylorFormSo)
1962
    #polySo = taylorFormListSaSo[0]
1963
    #errorRangeSo = sollya_lib_copy_obj(taylorFormListSaSo[2])
1964
    errorRangeSo = taylorFormListSaSo[2]
1965
    # No copy_obj needed here: a new objects are created.
1966
    maxErrorSo    = sollya_lib_sup(errorRangeSo)
1967
    minErrorSo    = sollya_lib_inf(errorRangeSo)
1968
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
1969
    absMinErrorSo = sollya_lib_abs(minErrorSo)
1970
    sollya_lib_clear_obj(maxErrorSo)
1971
    sollya_lib_clear_obj(minErrorSo)
1972
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
1973
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
1974
    # If changed, reset the Sollya working precision.
1975
    if sollyaPrecChanged:
1976
        sollya_lib_set_prec(initialSollyaPrecSo)
1977
    sollya_lib_clear_obj(initialSollyaPrecSo)
1978
    if errorTypeIsNone:
1979
        sollya_lib_clear_obj(errorTypeSo)
1980
    pobyso_clear_taylorform_sa_so(taylorFormListSaSo)
1981
    if absMaxErrorSa > absMinErrorSa:
1982
        sollya_lib_clear_obj(absMinErrorSo)
1983
        return (polySo, intervalCenterSo, absMaxErrorSo)
1984
    else:
1985
        sollya_lib_clear_obj(absMaxErrorSo)
1986
        return (polySo, intervalCenterSo, absMinErrorSo)
1987
# end pobyso_taylor_expansion_no_change_var_so_so
1988

    
1989
def pobyso_taylor_expansion_with_change_var_so_so(functionSo, degreeSo, \
1990
                                                  rangeSo, \
1991
                                                  errorTypeSo=None, \
1992
                                                  sollyaPrecSo=None):
1993
    """
1994
    Compute the Taylor expansion with the variable change
1995
    x -> (x-intervalCenter) included.
1996
    """
1997
    # Change Sollya internal precision, if need.
1998
    sollyaPrecChanged = False
1999
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_sos_sa()
2000
    if sollyaPrecSo is None:
2001
        pass
2002
    else:
2003
        sollya_lib_set_prec(sollyaPrecSo)
2004
        sollyaPrecChanged = True
2005
    #
2006
    # Error type stuff: default to absolute.
2007
    if errorTypeSo is None:
2008
        errorTypeIsNone = True
2009
        errorTypeSo = sollya_lib_absolute(None)
2010
    else:
2011
        errorTypeIsNone = False
2012
    intervalCenterSo = sollya_lib_mid(rangeSo)
2013
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo, \
2014
                                         intervalCenterSo, \
2015
                                         rangeSo, errorTypeSo, None)
2016
    # taylorFormListSaSo is a Python list of Sollya objects references that 
2017
    # are copies of the elements of taylorFormSo.
2018
    # pobyso_get_list_elements_so_so clears taylorFormSo.
2019
    (taylorFormListSo, numElements, isEndElliptic) = \
2020
        pobyso_get_list_elements_so_so(taylorFormSo)
2021
    polySo = taylorFormListSo[0]
2022
    errorRangeSo = taylorFormListSo[2]
2023
    maxErrorSo    = sollya_lib_sup(errorRangeSo)
2024
    minErrorSo    = sollya_lib_inf(errorRangeSo)
2025
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
2026
    absMinErrorSo = sollya_lib_abs(minErrorSo)
2027
    sollya_lib_clear_obj(maxErrorSo)
2028
    sollya_lib_clear_obj(minErrorSo)
2029
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
2030
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
2031
    changeVarExpSo = sollya_lib_build_function_sub(\
2032
                       sollya_lib_build_function_free_variable(),\
2033
                       sollya_lib_copy_obj(intervalCenterSo))
2034
    polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpSo)
2035
    sollya_lib_clear_obj(polySo) 
2036
    sollya_lib_clear_obj(changeVarExpSo)
2037
    # If changed, reset the Sollya working precision.
2038
    if sollyaPrecChanged:
2039
        sollya_lib_set_prec(initialSollyaPrecSo)
2040
    sollya_lib_clear_obj(initialSollyaPrecSo)
2041
    if errorTypeIsNone:
2042
        sollya_lib_clear_obj(errorTypeSo)
2043
    sollya_lib_clear_obj(taylorFormSo)
2044
    # Do not clear maxErrorSo.
2045
    if absMaxErrorSa > absMinErrorSa:
2046
        sollya_lib_clear_obj(absMinErrorSo)
2047
        return((polyVarChangedSo, intervalCenterSo, absMaxErrorSo))
2048
    else:
2049
        sollya_lib_clear_obj(absMaxErrorSo)
2050
        return((polyVarChangedSo, intervalCenterSo, absMinErrorSo))
2051
# end pobyso_taylor_expansion_with_change_var_so_so
2052

    
2053
def pobyso_taylor(function, degree, point):
2054
    """ Legacy function. See pobysoTaylor_so_so. """
2055
    return(pobyso_taylor_so_so(function, degree, point))
2056

    
2057
def pobyso_taylor_so_so(functionSo, degreeSo, pointSo):
2058
    return(sollya_lib_taylor(functionSo, degreeSo, pointSo))
2059
    
2060
def pobyso_taylorform(function, degree, point = None, 
2061
                      interval = None, errorType=None):
2062
    """ Legacy function. See pobyso_taylorform_sa_sa;"""
2063
    
2064
def pobyso_taylorform_sa_sa(functionSa, \
2065
                            degreeSa, \
2066
                            pointSa, \
2067
                            intervalSa=None, \
2068
                            errorTypeSa=None, \
2069
                            precisionSa=None):
2070
    """
2071
    Compute the Taylor form of 'degreeSa' for 'functionSa' at 'pointSa' 
2072
    for 'intervalSa' with 'errorTypeSa' (a string) using 'precisionSa'. 
2073
    point: must be a Real or a Real interval.
2074
    return the Taylor form as an array
2075
    TODO: take care of the interval and of the point when it is an interval;
2076
          when errorType is not None;
2077
          take care of the other elements of the Taylor form (coefficients 
2078
          errors and delta.
2079
    """
2080
    # Absolute as the default error.
2081
    if errorTypeSa is None:
2082
        errorTypeSo = sollya_lib_absolute()
2083
    elif errorTypeSa == "relative":
2084
        errorTypeSo = sollya_lib_relative()
2085
    elif errortypeSa == "absolute":
2086
        errorTypeSo = sollya_lib_absolute()
2087
    else:
2088
        # No clean up needed.
2089
        return None
2090
    # Global precision stuff
2091
    sollyaPrecisionChangedSa = False
2092
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
2093
    if precisionSa is None:
2094
        precSa = initialSollyaPrecSa
2095
    else:
2096
        if precSa > initialSollyaPrecSa:
2097
            if precSa <= 2:
2098
                print inspect.stack()[0][3], ":precision change <= 2 requested."
2099
            pobyso_set_prec_sa_so(precSa)
2100
            sollyaPrecisionChangedSa = True
2101
    #        
2102
    if len(functionSa.variables()) > 0:
2103
        varSa = functionSa.variables()[0]
2104
        pobyso_name_free_variable_sa_so(str(varSa))
2105
    # In any case (point or interval) the parent of pointSa has a precision
2106
    # method.
2107
    pointPrecSa = pointSa.parent().precision()
2108
    if precSa > pointPrecSa:
2109
        pointPrecSa = precSa
2110
    # In any case (point or interval) pointSa has a base_ring() method.
2111
    pointBaseRingString = str(pointSa.base_ring())
2112
    if re.search('Interval', pointBaseRingString) is None: # Point
2113
        pointSo = pobyso_constant_sa_so(pointSa, pointPrecSa)
2114
    else: # Interval.
2115
        pointSo = pobyso_interval_to_range_sa_so(pointSa, pointPrecSa)
2116
    # Sollyafy the function.
2117
    functionSo = pobyso_parse_string_sa_so(functionSa._assume_str().replace('_SAGE_VAR_', ''))
2118
    if sollya_lib_obj_is_error(functionSo):
2119
        print "pobyso_tailorform: function string can't be parsed!"
2120
        return None
2121
    # Sollyafy the degree
2122
    degreeSo = sollya_lib_constant_from_int(int(degreeSa))
2123
    # Sollyafy the point
2124
    # Call Sollya
2125
    taylorFormSo = \
2126
        sollya_lib_taylorform(functionSo, degreeSo, pointSo, errorTypeSo,\
2127
                                         None)
2128
    sollya_lib_clear_obj(functionSo)
2129
    sollya_lib_clear_obj(degreeSo)
2130
    sollya_lib_clear_obj(pointSo)
2131
    sollya_lib_clear_obj(errorTypeSo)
2132
    (tfsAsList, numElements, isEndElliptic) = \
2133
            pobyso_get_list_elements_so_so(taylorFormSo)
2134
    polySo = tfsAsList[0]
2135
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
2136
    polyRealField = RealField(maxPrecision)
2137
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, polyRealField)
2138
    if sollyaPrecisionChangedSa:
2139
        sollya_lib_set_prec(initialSollyaPrecSo)
2140
    sollya_lib_clear_obj(initialSollyaPrecSo)
2141
    polynomialRing = polyRealField[str(varSa)]
2142
    polySa = polynomial(expSa, polynomialRing)
2143
    taylorFormSa = [polySa]
2144
    # Final clean-up.
2145
    sollya_lib_clear_obj(taylorFormSo)
2146
    return(taylorFormSa)
2147
# End pobyso_taylor_form_sa_sa
2148

    
2149
def pobyso_taylorform_so_so(functionSo, degreeSo, pointSo, intervalSo=None, \
2150
                            errorTypeSo=None):
2151
    createdErrorType = False
2152
    if errorTypeSo is None:
2153
        errorTypeSo = sollya_lib_absolute()
2154
        createdErrorType = True
2155
    else:
2156
        #TODO: deal with the other case.
2157
        pass
2158
    if intervalSo is None:
2159
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
2160
                                         errorTypeSo, None)
2161
    else:
2162
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
2163
                                         intervalSo, errorTypeSo, None)
2164
    if createdErrorType:
2165
        sollya_lib_clear_obj(errorTypeSo)
2166
    return resultSo
2167
        
2168

    
2169
def pobyso_univar_polynomial_print_reverse(polySa):
2170
    """ Legacy function. See pobyso_univar_polynomial_print_reverse_sa_sa. """
2171
    return(pobyso_univar_polynomial_print_reverse_sa_sa(polySa))
2172

    
2173
def pobyso_univar_polynomial_print_reverse_sa_sa(polySa):
2174
    """
2175
    Return the string representation of a univariate polynomial with
2176
    monomials ordered in the x^0..x^n order of the monomials.
2177
    Remember: Sage
2178
    """
2179
    polynomialRing = polySa.base_ring()
2180
    # A very expensive solution:
2181
    # -create a fake multivariate polynomial field with only one variable,
2182
    #   specifying a negative lexicographical order;
2183
    mpolynomialRing = PolynomialRing(polynomialRing.base(), \
2184
                                     polynomialRing.variable_name(), \
2185
                                     1, order='neglex')
2186
    # - convert the univariate argument polynomial into a multivariate
2187
    #   version;
2188
    p = mpolynomialRing(polySa)
2189
    # - return the string representation of the converted form.
2190
    # There is no simple str() method defined for p's class.
2191
    return(p.__str__())
2192
#
2193
#print pobyso_get_prec()  
2194
pobyso_set_prec(165)
2195
#print pobyso_get_prec()  
2196
#a=100
2197
#print type(a)
2198
#id(a)
2199
#print "Max arity: ", pobyso_max_arity
2200
#print "Function tripleDouble (43) as a string: ", pobyso_function_type_as_string(43)
2201
#print "Function None (44) as a string: ", pobyso_function_type_as_string(44)
2202
print "...Pobyso check done"