Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 237

Historique | Voir | Annoter | Télécharger (83,96 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_dirty_find_zeros_so_so(funcSo, rangeSo):
418
    """
419
    Thin wrapper over sollya_lib_dirtyfindzeros()
420
    """
421
    return sollya_lib_dirtyfindzeros(funcSo, rangeSo)
422
# End pobys_dirty_find_zeros
423

    
424
def pobyso_dirty_inf_norm_so_so(funcSo, rangeSo, preSo=None):
425
    """
426
    Thin wrapper around sollya_dirtyinfnorm().
427
    """
428
    # TODO: manage the precision change.
429
    
430
    return sollya_lib_dirtyinfnorm(funcSo, rangeSo)
431
# End pobyso_dirty_inf_norm_so_so
432

    
433
def pobyso_float_list_so_sa(listSo):
434
    """
435
    Return a Sollya list of floating-point as a Sage list of 
436
    floating-point numbers.
437
    """
438
    listSa   = []
439
    ## Remember pobyso_get_list_elements_so_so returns more information
440
    #  than just the elements of the list (# elements, is_ellipti)
441
    listSaSo = pobyso_get_list_elements_so_so(listSo)[0]
442
    ## Search first for the maximum precision of the elements
443
    maxPrecSa = 0
444
    for floatSo in listSaSo:
445
        pobyso_autoprint(floatSo)
446
        curPrecSa =  pobyso_get_prec_of_constant_so_sa(floatSo)
447
        if curPrecSa > maxPrecSa:
448
            maxPrecSa = curPrecSa
449
    ##
450
    RF = RealField(maxPrecSa)
451
    ##
452
    for floatSo in listSaSo:
453
        listSa.append(pobyso_get_constant_as_rn_with_rf_so_sa(floatSo))
454
    return listSa
455
# End pobyso_float_list_so_sa
456

    
457
def pobyso_float_poly_sa_so(polySa, precSa = None):
458
    """
459
    Create a Sollya polynomial from a Sage RealField polynomial.
460
    """
461
    ## TODO: filter arguments.
462
    ## Precision. If a precision is given, convert the polynomial
463
    #  into the right polynomial field. If not convert it straight
464
    #  to Sollya.
465
    sollyaPrecChanged = False 
466
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
467
    if precSa is None:
468
        precSa = polySa.parent().base_ring().precision()
469
    if (precSa > initialSollyaPrecSa):
470
        assert precSa >= 2, "Precision change <2 requested"
471
        if precSa <= 2:
472
            print inspect.stack()[0][3], ": precision change <= 2 requested"
473
        precSo = pobyso_constant_from_int(precSa)
474
        pobyso_set_prec_so_so(precSo)
475
        sollya_lib_clear_obj(precSo)
476
        sollyaPrecChanged = True    
477
    ## Get exponents and coefficients.
478
    exponentsSa     = polySa.exponents()
479
    coefficientsSa  = polySa.coefficients()
480
    ## Build the polynomial.
481
    polySo = None
482
    for coefficientSa, exponentSa in zip(coefficientsSa, exponentsSa):
483
        #print coefficientSa.n(prec=precSa), exponentSa
484
        coefficientSo = \
485
            pobyso_constant_sa_so(coefficientSa)
486
        #pobyso_autoprint(coefficientSo)
487
        exponentSo = \
488
            pobyso_constant_from_int_sa_so(exponentSa)
489
        #pobyso_autoprint(exponentSo)
490
        monomialSo = sollya_lib_build_function_pow(
491
                       sollya_lib_build_function_free_variable(),
492
                       exponentSo)
493
        polyTermSo = sollya_lib_build_function_mul(coefficientSo,
494
                                                       monomialSo)
495
        if polySo is None:
496
            polySo = polyTermSo
497
        else:
498
            polySo = sollya_lib_build_function_add(polySo, polyTermSo)
499
    if sollyaPrecChanged:
500
        pobyso_set_prec_so_so(initialSollyaPrecSo)
501
        sollya_lib_clear_obj(initialSollyaPrecSo)
502
    return polySo
503
# End pobyso_float_poly_sa_so    
504

    
505
def pobyso_float_poly_so_sa(polySo, realFieldSa=None):
506
    """
507
    Convert a Sollya polynomial into a Sage floating-point polynomial.
508
    If no realField is given, a RealField corresponding to the maximum 
509
    precision of the coefficients is internally computed.
510
    The real field is not returned but can be easily retrieved from 
511
    the polynomial itself.
512
    ALGORITHM:
513
    - (optional) compute the RealField of the coefficients;
514
    - convert the Sollya expression into a Sage expression;
515
    - convert the Sage expression into a Sage polynomial
516
    """    
517
    if realFieldSa is None:
518
        expressionPrecSa = pobyso_get_max_prec_of_exp_so_sa(polySo)
519
        #print "Maximum precision of Sollya polynomial coefficients:", expressionPrecSa
520
        if expressionPrecSa < 2 or expressionPrecSa > 2147483391:
521
            print "Maximum degree of expression:", expressionPrecSa
522
        realFieldSa      = RealField(expressionPrecSa)
523
    #print "Sollya expression before...",
524
    #pobyso_autoprint(polySo)
525

    
526
    expressionSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo,
527
                                                             realFieldSa)
528
    #print "...Sollya expression after."
529
    #pobyso_autoprint(polySo)
530
    polyVariableSa = expressionSa.variables()[0]
531
    polyRingSa     = realFieldSa[str(polyVariableSa)]
532
    #print polyRingSa
533
    # Do not use the polynomial(expressionSa, ring=polyRingSa) form!
534
    polynomialSa = polyRingSa(expressionSa)
535
    polyCoeffsListSa = polynomialSa.coefficients()
536
    #for coeff in polyCoeffsListSa:
537
    #    print coeff.abs().n()
538
    return polynomialSa
539
# End pobyso_float_poly_so_sa
540

    
541
def pobyso_free_variable():
542
    """
543
    Ultra thin wrapper around the sollya_lib_function_build_free_variable function.
544
    """
545
    return sollya_lib_build_function_free_variable()
546
    
547
def pobyso_function_type_as_string(funcType):
548
    """ Legacy function. See pobyso_function_type_as_string_so_sa. """
549
    return(pobyso_function_type_as_string_so_sa(funcType))
550

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

    
647
def pobyso_get_constant(rnArgSa, constSo):
648
    """ Legacy function. See pobyso_get_constant_so_sa. """
649
    return pobyso_get_constant_so_sa(rnArgSa, constSo)
650
# End pobyso_get_constant
651

    
652
def pobyso_get_constant_so_sa(rnArgSa, constSo):
653
    """
654
    Set the value of rnArgSo to the value of constSo in MPFR_RNDN mode.
655
    rnArg must already exist and belong to some RealField.
656
    We assume that constSo points to a Sollya constant.
657
    """
658
    outcome = sollya_lib_get_constant(get_rn_value(rnArgSa), constSo)
659
    if outcome == 0: # Failure because constSo is not a constant expression.
660
        return None
661
    else:
662
        return outcome
663
# End  pobyso_get_constant_so_sa
664
   
665
def pobyso_get_constant_as_rn(ctExpSo):
666
    """ 
667
    Legacy function. See pobyso_get_constant_as_rn_so_sa. 
668
    """ 
669
    return(pobyso_get_constant_as_rn_so_sa(ctExpSo))
670
    
671
def pobyso_get_constant_as_rn_so_sa(constExpSo):
672
    """
673
    Get a Sollya constant as a Sage "real number".
674
    The precision of the floating-point number returned is that of the Sollya
675
    constant.
676
    """
677
    #print "Before computing precision of variable..."
678
    #pobyso_autoprint(constExpSo)
679
    precisionSa  = pobyso_get_prec_of_constant_so_sa(constExpSo)
680
    #print "precisionSa:", precisionSa
681
    ## If the expression can not be exactly converted, None is returned.
682
    #  In this case opt for the Sollya current expression.
683
    if precisionSa is None:
684
        precisionSa = pobyso_get_prec_so_sa()
685
    RRRR = RealField(precisionSa)
686
    rnSa = RRRR(0)
687
    outcome = sollya_lib_get_constant(get_rn_value(rnSa), constExpSo)
688
    if outcome == 0:
689
        return None
690
    else:
691
        return rnSa
692
# End pobyso_get_constant_as_rn_so_sa
693

    
694
def pobyso_get_constant_as_rn_with_rf(ctExp, realField):
695
    """ 
696
    Legacy function. See pobyso_get_constant_as_rn_with_rf_so_sa.
697
    """
698
    return pobyso_get_constant_as_rn_with_rf_so_sa(ctExp, realField)
699
# End pobyso_get_constant_as_rn_with_rf
700
    
701
def pobyso_get_constant_as_rn_with_rf_so_sa(ctExpSo, realFieldSa = None):
702
    """
703
    Get a Sollya constant as a Sage "real number".
704
    If no real field is specified, the precision of the floating-point number 
705
    returned is that of the Sollya constant.
706
    Otherwise is is that of the real field. Hence rounding may happen.
707
    """
708
    if realFieldSa is None:
709
        return pobyso_get_constant_as_rn_so_sa(ctExpSo)
710
    rnSa = realFieldSa(0)
711
    outcome = sollya_lib_get_constant(get_rn_value(rnSa), ctExpSo)
712
    if outcome == 0:
713
        return None
714
    else:
715
        return rnSa
716
# End pobyso_get_constant_as_rn_with_rf_so_sa
717

    
718
def pobyso_get_free_variable_name():
719
    """ 
720
    Legacy function. See pobyso_get_free_variable_name_so_sa.
721
    """
722
    return(pobyso_get_free_variable_name_so_sa())
723

    
724
def pobyso_get_free_variable_name_so_sa():
725
    return sollya_lib_get_free_variable_name()
726
    
727
def pobyso_get_function_arity(expressionSo):
728
    """ 
729
    Legacy function. See pobyso_get_function_arity_so_sa.
730
    """
731
    return(pobyso_get_function_arity_so_sa(expressionSo))
732

    
733
def pobyso_get_function_arity_so_sa(expressionSo):
734
    arity = c_int(0)
735
    sollya_lib_get_function_arity(byref(arity),expressionSo)
736
    return int(arity.value)
737

    
738
def pobyso_get_head_function(expressionSo):
739
    """ 
740
    Legacy function. See pobyso_get_head_function_so_sa. 
741
    """
742
    return(pobyso_get_head_function_so_sa(expressionSo)) 
743

    
744
def pobyso_get_head_function_so_sa(expressionSo):
745
    functionType = c_int(0)
746
    sollya_lib_get_head_function(byref(functionType), expressionSo)
747
    return int(functionType.value)
748

    
749
def pobyso_get_interval_from_range_so_sa(soRange, realIntervalFieldSa = None ):
750
    """
751
    Return the Sage interval corresponding to the Sollya range argument.
752
    If no reaIntervalField is passed as an argument, the interval bounds are not
753
    rounded: they are elements of RealIntervalField of the "right" precision
754
    to hold all the digits.
755
    """
756
    prec = c_int(0)
757
    if realIntervalFieldSa is None:
758
        retval = sollya_lib_get_prec_of_range(byref(prec), soRange, None)
759
        if retval == 0:
760
            return None
761
        realIntervalFieldSa = RealIntervalField(prec.value)
762
    intervalSa = realIntervalFieldSa(0,0)
763
    retval = \
764
        sollya_lib_get_interval_from_range(get_interval_value(intervalSa),\
765
                                           soRange)
766
    if retval == 0:
767
        return None
768
    return intervalSa
769
# End pobyso_get_interval_from_range_so_sa
770

    
771
def pobyso_get_list_elements(soObj):
772
    """ Legacy function. See pobyso_get_list_elements_so_so. """
773
    return pobyso_get_list_elements_so_so(soObj)
774
 
775
def pobyso_get_list_elements_so_so(objectListSo):
776
    """
777
    Get the Sollya list elements as a Sage/Python array of Sollya objects.
778
    
779
    INPUT:
780
    - objectListSo: a Sollya list of Sollya objects.
781
    
782
    OUTPUT:
783
    - a Sage/Python tuple made of:
784
      - a Sage/Python list of Sollya objects,
785
      - a Sage/Python int holding the number of elements,
786
      - a Sage/Python int stating (!= 0) that the list is end-elliptic.
787
    NOTE::
788
        We recover the addresses of the Sollya object from the list of pointers
789
        returned by sollya_lib_get_list_elements. The list itself is freed.
790
    TODO::
791
        Figure out what to do with numElements since the number of elements
792
        can easily be recovered from the list itself. 
793
        Ditto for isEndElliptic.
794
    """
795
    listAddress = POINTER(c_longlong)()
796
    numElements = c_int(0)
797
    isEndElliptic = c_int(0)
798
    listAsSageList = []
799
    result = sollya_lib_get_list_elements(byref(listAddress),\
800
                                          byref(numElements),\
801
                                          byref(isEndElliptic),\
802
                                          objectListSo)
803
    if result == 0 :
804
        return None
805
    for i in xrange(0, numElements.value, 1):
806
       #listAsSageList.append(sollya_lib_copy_obj(listAddress[i]))
807
       listAsSageList.append(listAddress[i])
808
       # Clear each of the elements returned by Sollya.
809
       #sollya_lib_clear_obj(listAddress[i])
810
    # Free the list itself.   
811
    sollya_lib_free(listAddress)
812
    return (listAsSageList, numElements.value, isEndElliptic.value)
813

    
814
def pobyso_get_max_prec_of_exp(soExp):
815
    """ Legacy function. See pobyso_get_max_prec_of_exp_so_sa. """
816
    return pobyso_get_max_prec_of_exp_so_sa(soExp)
817

    
818
def pobyso_get_max_prec_of_exp_so_sa(expSo):
819
    """
820
    Get the maximum precision used for the numbers in a Sollya expression.
821
    
822
    Arguments:
823
    soExp -- a Sollya expression pointer
824
    Return value:
825
    A Python integer
826
    TODO: 
827
    - error management;
828
    - correctly deal with numerical type such as DOUBLEEXTENDED.
829
    """
830
    if expSo is None:
831
        print inspect.stack()[0][3], ": expSo is None."
832
        return 0
833
    maxPrecision = 0
834
    minConstPrec = 0
835
    currentConstPrec = 0
836
    #pobyso_autoprint(expSo)
837
    operator = pobyso_get_head_function_so_sa(expSo)
838
    if (operator != SOLLYA_BASE_FUNC_CONSTANT) and \
839
    (operator != SOLLYA_BASE_FUNC_FREE_VARIABLE):
840
        (arity, subexpressions) = pobyso_get_subfunctions_so_sa(expSo)
841
        for i in xrange(arity):
842
            maxPrecisionCandidate = \
843
                pobyso_get_max_prec_of_exp_so_sa(subexpressions[i])
844
            if maxPrecisionCandidate > maxPrecision:
845
                maxPrecision = maxPrecisionCandidate
846
        return maxPrecision
847
    elif operator == SOLLYA_BASE_FUNC_CONSTANT:
848
        #minConstPrec = pobyso_get_min_prec_of_constant_so_sa(expSo)
849
        #currentConstPrec = pobyso_get_min_prec_of_constant_so_sa(soExp)
850
        #print minConstPrec, " - ", currentConstPrec 
851
        return pobyso_get_min_prec_of_constant_so_sa(expSo)
852
    
853
    elif operator == SOLLYA_BASE_FUNC_FREE_VARIABLE:
854
        return 0
855
    else:
856
        print "pobyso_get_max_prec_of_exp_so_sa: unexepected operator."
857
        return 0
858

    
859
def pobyso_get_min_prec_of_constant_so_sa(constExpSo):
860
    """
861
    Get the minimum precision necessary to represent the value of a Sollya
862
    constant.
863
    MPFR_MIN_PREC and powers of 2 are taken into account.
864
    We assume that constExpSo is a pointer to a Sollay constant expression.
865
    """
866
    constExpAsRnSa = pobyso_get_constant_as_rn_so_sa(constExpSo)
867
    return(min_mpfr_size(get_rn_value(constExpAsRnSa)))
868

    
869
def pobyso_get_poly_so_sa(polySo, realFieldSa=None):
870
    """
871
    Convert a Sollya polynomial into a Sage polynomial.
872
    Legacy function. Use pobyso_float_poly_so_sa() instead.
873
    """    
874
    return pobyso_float_poly_so_sa(polySo,realFieldSa)
875
# End pobyso_get_poly_so_sa
876

    
877
def pobyso_get_prec():
878
    """ Legacy function. See pobyso_get_prec_so_sa(). """
879
    return pobyso_get_prec_so_sa()
880

    
881
def pobyso_get_prec_so():
882
    """
883
    Get the current default precision in Sollya.
884
    The return value is a Sollya object.
885
    Usefull when modifying the precision back and forth by avoiding
886
    extra conversions.
887
    """
888
    return sollya_lib_get_prec(None)
889
    
890
def pobyso_get_prec_so_sa():
891
    """
892
    Get the current default precision in Sollya.
893
    The return value is Sage/Python int.
894
    """
895
    precSo = sollya_lib_get_prec()
896
    precSa = pobyso_constant_from_int_so_sa(precSo)
897
    sollya_lib_clear_obj(precSo)
898
    return precSa
899
# End pobyso_get_prec_so_sa.
900

    
901
def pobyso_get_prec_so_so_sa():
902
    """
903
    Return the current precision both as a Sollya object and a
904
    Sage integer as hybrid tuple.
905
    To avoid multiple calls for precision manipulations. 
906
    """
907
    precSo = sollya_lib_get_prec()
908
    precSa = pobyso_constant_from_int_so_sa(precSo)
909
    return (precSo, int(precSa))
910
    
911
def pobyso_get_prec_of_constant(ctExpSo):
912
    """ Legacy function. See pobyso_get_prec_of_constant_so_sa. """
913
    return pobyso_get_prec_of_constant_so_sa(ctExpSo)
914

    
915
def pobyso_get_prec_of_constant_so_sa(ctExpSo):
916
    """
917
    Tries to find a precision to represent ctExpSo without rounding.
918
    If not possible, returns None.
919
    """
920
    #print "Entering pobyso_get_prec_of_constant_so_sa..."
921
    prec = c_int(0)
922
    retc = sollya_lib_get_prec_of_constant(byref(prec), ctExpSo, None)
923
    if retc == 0:
924
        #print "pobyso_get_prec_of_constant_so_sa failed."
925
        return None
926
    #print "...exiting pobyso_get_prec_of_constant_so_sa."
927
    return int(prec.value)
928

    
929
def pobyso_get_prec_of_range_so_sa(rangeSo):
930
    """
931
    Returns the number of bits elements of a range are coded with.
932
    """
933
    prec = c_int(0)
934
    retc = sollya_lib_get_prec_of_range(byref(prec), rangeSo, None)
935
    if retc == 0:
936
        return(None)
937
    return int(prec.value)
938
# End pobyso_get_prec_of_range_so_sa()
939

    
940
def pobyso_get_sage_exp_from_sollya_exp(sollyaExpSo, realField = RR):
941
    """ Legacy function. See pobyso_get_sage_exp_from_sollya_exp_so_sa. """
942
    return pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, 
943
                                                     realField = RR)
944

    
945
def pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, realFieldSa = RR):
946
    """
947
    Get a Sage expression from a Sollya expression. 
948
    Currently only tested with polynomials with floating-point coefficients.
949
    Notice that, in the returned polynomial, the exponents are RealNumbers.
950
    """
951
    #pobyso_autoprint(sollyaExp)
952
    operatorSa = pobyso_get_head_function_so_sa(sollyaExpSo)
953
    sollyaLibFreeVariableName = sollya_lib_get_free_variable_name()
954
    ## Get rid of the "_"'s in "_x_", if any.
955
    sollyaLibFreeVariableName = re.sub('_', '', sollyaLibFreeVariableName)
956
    # Constants and the free variable are special cases.
957
    # All other operator are dealt with in the same way.
958
    if (operatorSa != SOLLYA_BASE_FUNC_CONSTANT) and \
959
       (operatorSa != SOLLYA_BASE_FUNC_FREE_VARIABLE):
960
        (aritySa, subexpressionsSa) = pobyso_get_subfunctions_so_sa(sollyaExpSo)
961
        if aritySa == 1:
962
            sageExpSa = eval(pobyso_function_type_as_string_so_sa(operatorSa) + \
963
            "(" + pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], \
964
            realFieldSa) + ")")
965
        elif aritySa == 2:
966
            # We do not get through the preprocessor.
967
            # The "^" operator is then a special case.
968
            if operatorSa == SOLLYA_BASE_FUNC_POW:
969
                operatorAsStringSa = "**"
970
            else:
971
                operatorAsStringSa = \
972
                    pobyso_function_type_as_string_so_sa(operatorSa)
973
            sageExpSa = \
974
              eval("pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], realFieldSa)"\
975
              + " " + operatorAsStringSa + " " + \
976
                   "pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[1], realFieldSa)")
977
        # We do not know yet how to deal with arity >= 3 
978
        # (is there any in Sollya anyway?).
979
        else:
980
            sageExpSa = eval('None')
981
        return sageExpSa
982
    elif operatorSa == SOLLYA_BASE_FUNC_CONSTANT:
983
        #print "This is a constant"
984
        return pobyso_get_constant_as_rn_with_rf_so_sa(sollyaExpSo, realFieldSa)
985
    elif operatorSa == SOLLYA_BASE_FUNC_FREE_VARIABLE:
986
        #print "This is the free variable"
987
        return eval(sollyaLibFreeVariableName)
988
    else:
989
        print "Unexpected"
990
        return eval('None')
991
# End pobyso_get_sage_exp_from_sollya_exp_so_sa
992

    
993

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

    
1089
def pobyso_guess_degree_so_sa(functionSo, rangeSo, errorSo, weightSo=None, \
1090
                              degreeBoundSo=None):
1091
    """
1092
    Thin wrapper around the guessdegree function.
1093
    Nevertheless, some precision control stuff has been appended.
1094
    """
1095
    # Deal with Sollya internal precision issues: if it is too small, 
1096
    # compared with the error, increases it to about twice -log2(error).
1097
    errorSa = pobyso_get_constant_as_rn_with_rf_so_sa(errorSo)
1098
    log2ErrorSa = errorSa.log2()
1099
    if log2ErrorSa < 0:
1100
        neededPrecisionSa = int(2 * int(-log2ErrorSa) / 64) * 64
1101
    else:
1102
        neededPrecisionSa = int(2 * int(log2ErrorSa) / 64) * 64
1103
    #print "Needed precision:", neededPrecisionSa
1104
    sollyaPrecisionHasChanged = False
1105
    (initialPrecSo, initialPrecSa) = pobyso_get_prec_so_so_sa()
1106
    if neededPrecisionSa > initialPrecSa:
1107
        if neededPrecisionSa <= 2:
1108
            print inspect.stack()[0][3], ": precision change <= 2 requested."
1109
        pobyso_set_prec_sa_so(neededPrecisionSa)
1110
        sollyaPrecisionHasChanged = True
1111
    #print "Guessing degree..."
1112
    # weightSo and degreeBoundsSo are optional arguments.
1113
    # As declared, sollya_lib_guessdegree must take 5 arguments.
1114
    if weightSo is None:
1115
        degreeRangeSo = sollya_lib_guessdegree(functionSo, rangeSo, errorSo, 
1116
                                               0, 0, None)
1117
    elif degreeBoundSo is None:
1118
        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, \
1119
                                                errorSo, weightSo, 0, None)
1120
    else:
1121
        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, errorSo, \
1122
                                                weightSo, degreeBoundSo, None)
1123
    #print "...degree guess done."
1124
    # Restore internal precision, if applicable.
1125
    if sollyaPrecisionHasChanged:
1126
        pobyso_set_prec_so_so(initialPrecSo)
1127
        sollya_lib_clear_obj(initialPrecSo)
1128
    degreeIntervalSa = pobyso_range_to_interval_so_sa(degreeRangeSo)
1129
    sollya_lib_clear_obj(degreeRangeSo)
1130
    # When ok, both bounds match.
1131
    # When the degree bound is too low, the upper bound is the degree
1132
    # for which the error can be honored.
1133
    # When it really goes wrong, the upper bound is infinity. 
1134
    if degreeIntervalSa.lower() == degreeIntervalSa.upper():
1135
        return int(degreeIntervalSa.lower())
1136
    else:
1137
        if degreeIntervalSa.upper().is_infinity():
1138
            return None
1139
        else:
1140
            return int(degreeIntervalSa.upper())
1141
    # End pobyso_guess_degree_so_sa
1142

    
1143
def pobyso_inf_so_so(intervalSo):
1144
    """
1145
    Very thin wrapper around sollya_lib_inf().
1146
    """
1147
    return sollya_lib_inf(intervalSo)
1148
# End pobyso_inf_so_so.
1149
    
1150
def pobyso_infnorm_so_so(func, interval, file = None, intervalList = None):
1151
    print "Do not use this function. User pobyso_supnorm_so_so instead."
1152
    return None
1153

    
1154
def pobyso_interval_to_range_sa_so(intervalSa, precisionSa=None):
1155
    if precisionSa is None:
1156
        precisionSa = intervalSa.parent().precision()
1157
    intervalSo = pobyso_bounds_to_range_sa_so(intervalSa.lower(),\
1158
                                              intervalSa.upper(),\
1159
                                              precisionSa)
1160
    return intervalSo
1161
# End pobyso_interval_to_range_sa_so
1162

    
1163
def pobyso_is_error_so_sa(objSo):
1164
    """
1165
    Thin wrapper around the sollya_lib_obj_is_error() function.
1166
    """
1167
    if sollya_lib_obj_is_error(objSo) != 0:
1168
        return True
1169
    else:
1170
        return False
1171
# End pobyso_is_error-so_sa
1172

    
1173
def pobyso_is_floating_point_number_sa_sa(numberSa):
1174
    """
1175
    Check whether a Sage number is floating point.
1176
    Exception stuff added because numbers other than
1177
    floating-point ones do not have the is_real() attribute.
1178
    """
1179
    try:
1180
        return numberSa.is_real()
1181
    except AttributeError:
1182
        return False
1183
# End pobyso_is_floating_piont_number_sa_sa
1184

    
1185
def pobyso_is_range_so_sa(rangeCandidateSo):
1186
    """
1187
    Thin wrapper over sollya_lib_is_range.
1188
    """
1189
    return sollya_lib_obj_is_range(rangeCandidateSo) != 0
1190
    
1191
# End pobyso_is_range_so_sa
1192

    
1193

    
1194
def pobyso_lib_init():
1195
    sollya_lib_init(None)
1196

    
1197
def pobyso_lib_close():
1198
    sollya_lib_close(None)
1199

    
1200
def pobyso_name_free_variable(freeVariableNameSa):
1201
    """ Legacy function. See pobyso_name_free_variable_sa_so. """
1202
    pobyso_name_free_variable_sa_so(freeVariableNameSa)
1203

    
1204
def pobyso_name_free_variable_sa_so(freeVariableNameSa):
1205
    """
1206
    Set the free variable name in Sollya from a Sage string.
1207
    """
1208
    sollya_lib_name_free_variable(freeVariableNameSa)
1209

    
1210
def pobyso_parse_string(string):
1211
    """ Legacy function. See pobyso_parse_string_sa_so. """
1212
    return pobyso_parse_string_sa_so(string)
1213
 
1214
def pobyso_parse_string_sa_so(string):
1215
    """
1216
    Get the Sollya expression computed from a Sage string or
1217
    a Sollya error object if parsing failed.
1218
    """
1219
    return sollya_lib_parse_string(string)
1220

    
1221
def pobyso_precision_so_sa(ctExpSo):
1222
    """
1223
    Computes the necessary precision to represent a number.
1224
    If x is not zero, it can be uniquely written as x = m · 2e 
1225
    where m is an odd integer and e is an integer. 
1226
    precision(x) returns the number of bits necessary to write m 
1227
    in binary (i.e. ceil(log2(m))).
1228
    """
1229
    #TODO: take care of the special case: 0, @NaN@, @Inf@
1230
    precisionSo = sollya_lib_precision(ctExpSo)
1231
    precisionSa = pobyso_constant_from_int_so_sa(precisionSo)
1232
    sollya_lib_clear_obj(precisionSo)
1233
    return precisionSa
1234
# End pobyso_precision_so_sa
1235

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

    
1402
def pobyso_polynomial_degree_so_sa(polySo):
1403
    """
1404
    Return the degree of a Sollya polynomial as a Sage int.
1405
    """
1406
    degreeSo = sollya_lib_degree(polySo)
1407
    return pobyso_constant_from_int_so_sa(degreeSo)
1408
# End pobyso_polynomial_degree_so_sa
1409

    
1410
def pobyso_polynomial_degree_so_so(polySo):
1411
    """
1412
    Thin wrapper around lib_sollya_degree().
1413
    """
1414
    return sollya_lib_degree(polySo)
1415
# End pobyso_polynomial_degree_so_so
1416
                                                                  
1417
def pobyso_range(rnLowerBound, rnUpperBound):
1418
    """ Legacy function. See pobyso_range_sa_so. """
1419
    return pobyso_range_sa_so(rnLowerBound, rnUpperBound) 
1420

    
1421
def pobyso_range_from_bounds_sa_so(rnLowerBound, rnUpperBound, precSa = None):
1422
    """
1423
    Create a Sollya range from 2 Sage real numbers as bounds
1424
    """
1425
    # TODO check precision stuff.
1426
    sollyaPrecChanged = False
1427
    (initialSollyaPrecSo, initialSollyaPrecSa) = \
1428
        pobyso_get_prec_so_so_sa()
1429
    if precSa is None:
1430
        precSa = max(rnLowerBound.parent().prec(), rnUpperBound.parent().prec())
1431
    if precSa > initialSollyaPrecSa:
1432
        if precSa <= 2:
1433
            print inspect.stack()[0][3], ": precision change <= 2 requested."
1434
        pobyso_set_prec_sa_so(precSa)
1435
        sollyaPrecChanged = True
1436
    rangeSo = sollya_lib_range_from_bounds(get_rn_value(rnLowerBound),
1437
                                           get_rn_value(rnUpperBound))
1438
    if sollyaPrecChanged:
1439
        pobyso_set_prec_so_so(initialSollyaPrecSo)
1440
        sollya_lib_clear_obj(initialSollyaPrecSo)
1441
    return rangeSo
1442
# End pobyso_range_from_bounds_sa_so
1443

    
1444
def pobyso_range_max_abs_so_so(rangeSo):
1445
    """
1446
    Return, as Sollya constant, the maximum absolute value of a Sollay range.
1447
    """
1448
    lowerBoundSo = sollya_lib_inf(rangeSo)
1449
    upperBoundSo = sollya_lib_sup(rangeSo)
1450
    #
1451
    lowerBoundSo = sollya_lib_build_function_abs(lowerBoundSo)
1452
    upperBoundSo = sollya_lib_build_function_abs(upperBoundSo)
1453
    #pobyso_autoprint(lowerBoundSo)
1454
    #pobyso_autoprint(upperBoundSo)
1455
    #
1456
    maxAbsSo = sollya_lib_max(lowerBoundSo, upperBoundSo, None)
1457
    #sollya_lib_clear_obj(lowerBoundSo)
1458
    #sollya_lib_clear_obj(upperBoundSo)
1459
    return maxAbsSo
1460
# End pobyso_range_max_abs_so_so
1461
 
1462
def pobyso_range_to_interval_so_sa(rangeSo, realIntervalFieldSa = None):
1463
    """
1464
    Get a Sage interval from a Sollya range.
1465
    If no realIntervalField is given as a parameter, the Sage interval
1466
    precision is that of the Sollya range.
1467
    Otherwise, the precision is that of the realIntervalField. In this case
1468
    rounding may happen.
1469
    """
1470
    if realIntervalFieldSa is None:
1471
        precSa = pobyso_get_prec_of_range_so_sa(rangeSo)
1472
        realIntervalFieldSa = RealIntervalField(precSa)
1473
    intervalSa = \
1474
        pobyso_get_interval_from_range_so_sa(rangeSo, realIntervalFieldSa) 
1475
    return intervalSa
1476
# End pobyso_range_to_interval_so_sa
1477

    
1478
def pobyso_rat_poly_sa_so(polySa, precSa = None):
1479
    """
1480
    Create a Sollya polynomial from a Sage rational polynomial.
1481
    """
1482
    ## TODO: filter arguments.
1483
    ## Precision. If no precision is given, use the current precision
1484
    #  of Sollya.
1485
    if precSa is None:
1486
        precSa =  pobyso_get_prec_so_sa()
1487
    #print "Precision:",  precSa
1488
    RRR = RealField(precSa)
1489
    ## Create a Sage polynomial in the "right" precision.
1490
    P_RRR = RRR[polySa.variables()[0]]
1491
    polyFloatSa = P_RRR(polySa)
1492
    ## Make sure no precision is provided: pobyso_float_poly_sa_so will
1493
    #  recover it all by itself and not make an extra conversion.
1494
    return pobyso_float_poly_sa_so(polyFloatSa)
1495
    
1496
# End pobyso_rat_poly_sa_so    
1497
    
1498
def pobyso_remez_canonical_sa_sa(func, \
1499
                                 degree, \
1500
                                 lowerBound, \
1501
                                 upperBound, \
1502
                                 weight = None, \
1503
                                 quality = None):
1504
    """
1505
    All arguments are Sage/Python.
1506
    The functions (func and weight) must be passed as expressions or strings.
1507
    Otherwise the function fails. 
1508
    The return value is a Sage polynomial.
1509
    """
1510
    var('zorglub')    # Dummy variable name for type check only. Type of 
1511
    # zorglub is "symbolic expression".
1512
    polySo = pobyso_remez_canonical_sa_so(func, \
1513
                                 degree, \
1514
                                 lowerBound, \
1515
                                 upperBound, \
1516
                                 weight, \
1517
                                 quality)
1518
    # String test
1519
    if parent(func) == parent("string"):
1520
        functionSa = eval(func)
1521
    # Expression test.
1522
    elif type(func) == type(zorglub):
1523
        functionSa = func
1524
    else:
1525
        return None
1526
    #
1527
    maxPrecision = 0
1528
    if polySo is None:
1529
        return(None)
1530
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
1531
    RRRRSa = RealField(maxPrecision)
1532
    polynomialRingSa = RRRRSa[functionSa.variables()[0]]
1533
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, RRRRSa)
1534
    polySa = polynomial(expSa, polynomialRingSa)
1535
    sollya_lib_clear_obj(polySo)
1536
    return(polySa)
1537
# End pobyso_remez_canonical_sa_sa
1538
    
1539
def pobyso_remez_canonical(func, \
1540
                           degree, \
1541
                           lowerBound, \
1542
                           upperBound, \
1543
                           weight = "1", \
1544
                           quality = None):
1545
    """ Legacy function. See pobyso_remez_canonical_sa_so. """
1546
    return(pobyso_remez_canonical_sa_so(func, \
1547
                                        degree, \
1548
                                        lowerBound, \
1549
                                        upperBound, \
1550
                                        weight, \
1551
                                        quality))
1552
# End pobyso_remez_canonical.
1553

    
1554
def pobyso_remez_canonical_sa_so(func, \
1555
                                 degree, \
1556
                                 lowerBound, \
1557
                                 upperBound, \
1558
                                 weight = None, \
1559
                                 quality = None):
1560
    """
1561
    All arguments are Sage/Python.
1562
    The functions (func and weight) must be passed as expressions or strings.
1563
    Otherwise the function fails. 
1564
    The return value is a pointer to a Sollya function.
1565
    lowerBound and upperBound mus be reals.
1566
    """
1567
    var('zorglub')    # Dummy variable name for type check only. Type of
1568
    # zorglub is "symbolic expression".
1569
    currentVariableNameSa = None
1570
    # The func argument can be of different types (string, 
1571
    # symbolic expression...)
1572
    if parent(func) == parent("string"):
1573
        localFuncSa = sage_eval(func,globals())
1574
        if len(localFuncSa.variables()) > 0:
1575
            currentVariableNameSa = localFuncSa.variables()[0]
1576
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1577
            functionSo = \
1578
              sollya_lib_parse_string(localFuncSa._assume_str().replace('_SAGE_VAR_', ''))
1579
    # Expression test.
1580
    elif type(func) == type(zorglub):
1581
        # Until we are able to translate Sage expressions into Sollya 
1582
        # expressions : parse the string version.
1583
        if len(func.variables()) > 0:
1584
            currentVariableNameSa = func.variables()[0]
1585
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1586
            functionSo = \
1587
              sollya_lib_parse_string(func._assume_str().replace('_SAGE_VAR_', ''))
1588
    else:
1589
        return(None)
1590
    if weight is None: # No weight given -> 1.
1591
        weightSo = pobyso_constant_1_sa_so()
1592
    elif parent(weight) == parent("string"): # Weight given as string: parse it.
1593
        weightSo = sollya_lib_parse_string(func)
1594
    elif type(weight) == type(zorglub): # Weight given as symbolice expression.
1595
        functionSo = \
1596
          sollya_lib_parse_string_sa_so(weight._assume_str().replace('_SAGE_VAR_', ''))
1597
    else:
1598
        return(None)
1599
    degreeSo = pobyso_constant_from_int(degree)
1600
    rangeSo = pobyso_bounds_to_range_sa_so(lowerBound, upperBound)
1601
    if not quality is None:
1602
        qualitySo= pobyso_constant_sa_so(quality)
1603
    else:
1604
        qualitySo = None
1605
        
1606
    remezPolySo = sollya_lib_remez(functionSo, \
1607
                                   degreeSo, \
1608
                                   rangeSo, \
1609
                                   weightSo, \
1610
                                   qualitySo, \
1611
                                   None)
1612
    sollya_lib_clear_obj(functionSo)
1613
    sollya_lib_clear_obj(degreeSo)
1614
    sollya_lib_clear_obj(rangeSo)
1615
    sollya_lib_clear_obj(weightSo)
1616
    if not qualitySo is None:
1617
        sollya_lib_clear_obj(qualitySo)
1618
    return(remezPolySo)
1619
# End pobyso_remez_canonical_sa_so
1620

    
1621
def pobyso_remez_canonical_so_so(funcSo, \
1622
                                 degreeSo, \
1623
                                 rangeSo, \
1624
                                 weightSo = pobyso_constant_1_sa_so(),\
1625
                                 qualitySo = None):
1626
    """
1627
    All arguments are pointers to Sollya objects.
1628
    The return value is a pointer to a Sollya function.
1629
    """
1630
    if not sollya_lib_obj_is_function(funcSo):
1631
        return(None)
1632
    return(sollya_lib_remez(funcSo, degreeSo, rangeSo, weightSo, qualitySo, None))
1633
# End pobyso_remez_canonical_so_so.
1634

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

    
1704

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

    
1833
def pobyso_set_canonical_off():
1834
    sollya_lib_set_canonical(sollya_lib_off())
1835

    
1836
def pobyso_set_canonical_on():
1837
    sollya_lib_set_canonical(sollya_lib_on())
1838

    
1839
def pobyso_set_prec(p):
1840
    """ Legacy function. See pobyso_set_prec_sa_so. """
1841
    pobyso_set_prec_sa_so(p)
1842

    
1843
def pobyso_set_prec_sa_so(p):
1844
    #a = c_int(p)
1845
    #precSo = c_void_p(sollya_lib_constant_from_int(a))
1846
    #precSo = sollya_lib_constant_from_int(a)
1847
    precSo = pobyso_constant_from_int_sa_so(p)
1848
    sollya_lib_set_prec(precSo)
1849
    sollya_lib_clear_obj(precSo)
1850
# End pobyso_set_prec_sa_so.
1851

    
1852
def pobyso_set_prec_so_so(newPrecSo):
1853
    sollya_lib_set_prec(newPrecSo)
1854
# End pobyso_set_prec_so_so.
1855

    
1856
def pobyso_inf_so_so(intervalSo):
1857
    """
1858
    Very thin wrapper around sollya_lib_inf().
1859
    """
1860
    return sollya_lib_inf(intervalSo)
1861
# End pobyso_inf_so_so.
1862
    
1863
def pobyso_supnorm_so_so(polySo, funcSo, intervalSo, errorTypeSo = None,\
1864
                         accuracySo = None):
1865
    """
1866
    Computes the supnorm of the approximation error between the given 
1867
    polynomial and function.
1868
    errorTypeSo defaults to "absolute".
1869
    accuracySo defaults to 2^(-40).
1870
    """
1871
    if errorTypeSo is None:
1872
        errorTypeSo = sollya_lib_absolute(None)
1873
        errorTypeIsNone = True
1874
    else:
1875
        errorTypeIsNone = False
1876
    #
1877
    if accuracySo is None:
1878
        # Notice the **: we are in Pythonland here!
1879
        accuracySo = pobyso_constant_sa_so(RR(2**(-40)))
1880
        accuracyIsNone = True
1881
    else:
1882
        accuracyIsNone = False
1883
    pobyso_autoprint(accuracySo)
1884
    resultSo = \
1885
        sollya_lib_supnorm(polySo, funcSo, intervalSo, errorTypeSo, \
1886
                              accuracySo)
1887
    if errorTypeIsNone:
1888
        sollya_lib_clear_obj(errorTypeSo)
1889
    if accuracyIsNone:
1890
        sollya_lib_clear_obj(accuracySo)
1891
    return resultSo
1892
# End pobyso_supnorm_so_so
1893

    
1894
def pobyso_taylor_expansion_no_change_var_so_so(functionSo, 
1895
                                                degreeSo, 
1896
                                                rangeSo,
1897
                                                errorTypeSo=None, 
1898
                                                sollyaPrecSo=None):
1899
    """
1900
    Compute the Taylor expansion without the variable change
1901
    x -> x-intervalCenter.
1902
    """
1903
    # Change internal Sollya precision, if needed.
1904
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
1905
    sollyaPrecChanged = False
1906
    if sollyaPrecSo is None:
1907
        pass
1908
    else:
1909
        sollya_lib_set_prec(sollyaPrecSo)
1910
        sollyaPrecChanged = True
1911
    # Error type stuff: default to absolute.
1912
    if errorTypeSo is None:
1913
        errorTypeIsNone = True
1914
        errorTypeSo = sollya_lib_absolute(None)
1915
    else:
1916
        errorTypeIsNone = False
1917
    intervalCenterSo = sollya_lib_mid(rangeSo, None)
1918
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo,
1919
                                         intervalCenterSo,
1920
                                         rangeSo, errorTypeSo, None)
1921
    # taylorFormListSaSo is a Python list of Sollya objects references that 
1922
    # are copies of the elements of taylorFormSo.
1923
    # pobyso_get_list_elements_so_so clears taylorFormSo.
1924
    (taylorFormListSaSo, numElementsSa, isEndEllipticSa) = \
1925
        pobyso_get_list_elements_so_so(taylorFormSo)
1926
    polySo = sollya_lib_copy_obj(taylorFormListSaSo[0])
1927
    #print "Num elements:", numElementsSa
1928
    sollya_lib_clear_obj(taylorFormSo)
1929
    #polySo = taylorFormListSaSo[0]
1930
    #errorRangeSo = sollya_lib_copy_obj(taylorFormListSaSo[2])
1931
    errorRangeSo = taylorFormListSaSo[2]
1932
    # No copy_obj needed here: a new objects are created.
1933
    maxErrorSo    = sollya_lib_sup(errorRangeSo)
1934
    minErrorSo    = sollya_lib_inf(errorRangeSo)
1935
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
1936
    absMinErrorSo = sollya_lib_abs(minErrorSo)
1937
    sollya_lib_clear_obj(maxErrorSo)
1938
    sollya_lib_clear_obj(minErrorSo)
1939
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
1940
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
1941
    # If changed, reset the Sollya working precision.
1942
    if sollyaPrecChanged:
1943
        sollya_lib_set_prec(initialSollyaPrecSo)
1944
    sollya_lib_clear_obj(initialSollyaPrecSo)
1945
    if errorTypeIsNone:
1946
        sollya_lib_clear_obj(errorTypeSo)
1947
    pobyso_clear_taylorform_sa_so(taylorFormListSaSo)
1948
    if absMaxErrorSa > absMinErrorSa:
1949
        sollya_lib_clear_obj(absMinErrorSo)
1950
        return (polySo, intervalCenterSo, absMaxErrorSo)
1951
    else:
1952
        sollya_lib_clear_obj(absMaxErrorSo)
1953
        return (polySo, intervalCenterSo, absMinErrorSo)
1954
# end pobyso_taylor_expansion_no_change_var_so_so
1955

    
1956
def pobyso_taylor_expansion_with_change_var_so_so(functionSo, degreeSo, \
1957
                                                  rangeSo, \
1958
                                                  errorTypeSo=None, \
1959
                                                  sollyaPrecSo=None):
1960
    """
1961
    Compute the Taylor expansion with the variable change
1962
    x -> (x-intervalCenter) included.
1963
    """
1964
    # Change Sollya internal precision, if need.
1965
    sollyaPrecChanged = False
1966
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_sos_sa()
1967
    if sollyaPrecSo is None:
1968
        pass
1969
    else:
1970
        sollya_lib_set_prec(sollyaPrecSo)
1971
        sollyaPrecChanged = True
1972
    #
1973
    # Error type stuff: default to absolute.
1974
    if errorTypeSo is None:
1975
        errorTypeIsNone = True
1976
        errorTypeSo = sollya_lib_absolute(None)
1977
    else:
1978
        errorTypeIsNone = False
1979
    intervalCenterSo = sollya_lib_mid(rangeSo)
1980
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo, \
1981
                                         intervalCenterSo, \
1982
                                         rangeSo, errorTypeSo, None)
1983
    # taylorFormListSaSo is a Python list of Sollya objects references that 
1984
    # are copies of the elements of taylorFormSo.
1985
    # pobyso_get_list_elements_so_so clears taylorFormSo.
1986
    (taylorFormListSo, numElements, isEndElliptic) = \
1987
        pobyso_get_list_elements_so_so(taylorFormSo)
1988
    polySo = taylorFormListSo[0]
1989
    errorRangeSo = taylorFormListSo[2]
1990
    maxErrorSo    = sollya_lib_sup(errorRangeSo)
1991
    minErrorSo    = sollya_lib_inf(errorRangeSo)
1992
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
1993
    absMinErrorSo = sollya_lib_abs(minErrorSo)
1994
    sollya_lib_clear_obj(maxErrorSo)
1995
    sollya_lib_clear_obj(minErrorSo)
1996
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
1997
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
1998
    changeVarExpSo = sollya_lib_build_function_sub(\
1999
                       sollya_lib_build_function_free_variable(),\
2000
                       sollya_lib_copy_obj(intervalCenterSo))
2001
    polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpSo)
2002
    sollya_lib_clear_obj(polySo) 
2003
    sollya_lib_clear_obj(changeVarExpSo)
2004
    # If changed, reset the Sollya working precision.
2005
    if sollyaPrecChanged:
2006
        sollya_lib_set_prec(initialSollyaPrecSo)
2007
    sollya_lib_clear_obj(initialSollyaPrecSo)
2008
    if errorTypeIsNone:
2009
        sollya_lib_clear_obj(errorTypeSo)
2010
    sollya_lib_clear_obj(taylorFormSo)
2011
    # Do not clear maxErrorSo.
2012
    if absMaxErrorSa > absMinErrorSa:
2013
        sollya_lib_clear_obj(absMinErrorSo)
2014
        return((polyVarChangedSo, intervalCenterSo, absMaxErrorSo))
2015
    else:
2016
        sollya_lib_clear_obj(absMaxErrorSo)
2017
        return((polyVarChangedSo, intervalCenterSo, absMinErrorSo))
2018
# end pobyso_taylor_expansion_with_change_var_so_so
2019

    
2020
def pobyso_taylor(function, degree, point):
2021
    """ Legacy function. See pobysoTaylor_so_so. """
2022
    return(pobyso_taylor_so_so(function, degree, point))
2023

    
2024
def pobyso_taylor_so_so(functionSo, degreeSo, pointSo):
2025
    return(sollya_lib_taylor(functionSo, degreeSo, pointSo))
2026
    
2027
def pobyso_taylorform(function, degree, point = None, 
2028
                      interval = None, errorType=None):
2029
    """ Legacy function. See pobyso_taylorform_sa_sa;"""
2030
    
2031
def pobyso_taylorform_sa_sa(functionSa, \
2032
                            degreeSa, \
2033
                            pointSa, \
2034
                            intervalSa=None, \
2035
                            errorTypeSa=None, \
2036
                            precisionSa=None):
2037
    """
2038
    Compute the Taylor form of 'degreeSa' for 'functionSa' at 'pointSa' 
2039
    for 'intervalSa' with 'errorTypeSa' (a string) using 'precisionSa'. 
2040
    point: must be a Real or a Real interval.
2041
    return the Taylor form as an array
2042
    TODO: take care of the interval and of the point when it is an interval;
2043
          when errorType is not None;
2044
          take care of the other elements of the Taylor form (coefficients 
2045
          errors and delta.
2046
    """
2047
    # Absolute as the default error.
2048
    if errorTypeSa is None:
2049
        errorTypeSo = sollya_lib_absolute()
2050
    elif errorTypeSa == "relative":
2051
        errorTypeSo = sollya_lib_relative()
2052
    elif errortypeSa == "absolute":
2053
        errorTypeSo = sollya_lib_absolute()
2054
    else:
2055
        # No clean up needed.
2056
        return None
2057
    # Global precision stuff
2058
    sollyaPrecisionChangedSa = False
2059
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
2060
    if precisionSa is None:
2061
        precSa = initialSollyaPrecSa
2062
    else:
2063
        if precSa > initialSollyaPrecSa:
2064
            if precSa <= 2:
2065
                print inspect.stack()[0][3], ":precision change <= 2 requested."
2066
            pobyso_set_prec_sa_so(precSa)
2067
            sollyaPrecisionChangedSa = True
2068
    #        
2069
    if len(functionSa.variables()) > 0:
2070
        varSa = functionSa.variables()[0]
2071
        pobyso_name_free_variable_sa_so(str(varSa))
2072
    # In any case (point or interval) the parent of pointSa has a precision
2073
    # method.
2074
    pointPrecSa = pointSa.parent().precision()
2075
    if precSa > pointPrecSa:
2076
        pointPrecSa = precSa
2077
    # In any case (point or interval) pointSa has a base_ring() method.
2078
    pointBaseRingString = str(pointSa.base_ring())
2079
    if re.search('Interval', pointBaseRingString) is None: # Point
2080
        pointSo = pobyso_constant_sa_so(pointSa, pointPrecSa)
2081
    else: # Interval.
2082
        pointSo = pobyso_interval_to_range_sa_so(pointSa, pointPrecSa)
2083
    # Sollyafy the function.
2084
    functionSo = pobyso_parse_string_sa_so(functionSa._assume_str().replace('_SAGE_VAR_', ''))
2085
    if sollya_lib_obj_is_error(functionSo):
2086
        print "pobyso_tailorform: function string can't be parsed!"
2087
        return None
2088
    # Sollyafy the degree
2089
    degreeSo = sollya_lib_constant_from_int(int(degreeSa))
2090
    # Sollyafy the point
2091
    # Call Sollya
2092
    taylorFormSo = \
2093
        sollya_lib_taylorform(functionSo, degreeSo, pointSo, errorTypeSo,\
2094
                                         None)
2095
    sollya_lib_clear_obj(functionSo)
2096
    sollya_lib_clear_obj(degreeSo)
2097
    sollya_lib_clear_obj(pointSo)
2098
    sollya_lib_clear_obj(errorTypeSo)
2099
    (tfsAsList, numElements, isEndElliptic) = \
2100
            pobyso_get_list_elements_so_so(taylorFormSo)
2101
    polySo = tfsAsList[0]
2102
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
2103
    polyRealField = RealField(maxPrecision)
2104
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, polyRealField)
2105
    if sollyaPrecisionChangedSa:
2106
        sollya_lib_set_prec(initialSollyaPrecSo)
2107
    sollya_lib_clear_obj(initialSollyaPrecSo)
2108
    polynomialRing = polyRealField[str(varSa)]
2109
    polySa = polynomial(expSa, polynomialRing)
2110
    taylorFormSa = [polySa]
2111
    # Final clean-up.
2112
    sollya_lib_clear_obj(taylorFormSo)
2113
    return(taylorFormSa)
2114
# End pobyso_taylor_form_sa_sa
2115

    
2116
def pobyso_taylorform_so_so(functionSo, degreeSo, pointSo, intervalSo=None, \
2117
                            errorTypeSo=None):
2118
    createdErrorType = False
2119
    if errorTypeSo is None:
2120
        errorTypeSo = sollya_lib_absolute()
2121
        createdErrorType = True
2122
    else:
2123
        #TODO: deal with the other case.
2124
        pass
2125
    if intervalSo is None:
2126
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
2127
                                         errorTypeSo, None)
2128
    else:
2129
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
2130
                                         intervalSo, errorTypeSo, None)
2131
    if createdErrorType:
2132
        sollya_lib_clear_obj(errorTypeSo)
2133
    return resultSo
2134
        
2135

    
2136
def pobyso_univar_polynomial_print_reverse(polySa):
2137
    """ Legacy function. See pobyso_univar_polynomial_print_reverse_sa_sa. """
2138
    return(pobyso_univar_polynomial_print_reverse_sa_sa(polySa))
2139

    
2140
def pobyso_univar_polynomial_print_reverse_sa_sa(polySa):
2141
    """
2142
    Return the string representation of a univariate polynomial with
2143
    monomials ordered in the x^0..x^n order of the monomials.
2144
    Remember: Sage
2145
    """
2146
    polynomialRing = polySa.base_ring()
2147
    # A very expensive solution:
2148
    # -create a fake multivariate polynomial field with only one variable,
2149
    #   specifying a negative lexicographical order;
2150
    mpolynomialRing = PolynomialRing(polynomialRing.base(), \
2151
                                     polynomialRing.variable_name(), \
2152
                                     1, order='neglex')
2153
    # - convert the univariate argument polynomial into a multivariate
2154
    #   version;
2155
    p = mpolynomialRing(polySa)
2156
    # - return the string representation of the converted form.
2157
    # There is no simple str() method defined for p's class.
2158
    return(p.__str__())
2159
#
2160
#print pobyso_get_prec()  
2161
pobyso_set_prec(165)
2162
#print pobyso_get_prec()  
2163
#a=100
2164
#print type(a)
2165
#id(a)
2166
#print "Max arity: ", pobyso_max_arity
2167
#print "Function tripleDouble (43) as a string: ", pobyso_function_type_as_string(43)
2168
#print "Function None (44) as a string: ", pobyso_function_type_as_string(44)
2169
print "...Pobyso check done"