Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 239

Historique | Voir | Annoter | Télécharger (84,32 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
    ## The function returns none if the list is empty or an error has happened.
440
    retVal = pobyso_get_list_elements_so_so(listSo)
441
    if retVal is None:
442
        return listSa
443
    ## Just in case the interface is changed and an empty list is returned.
444
    elif len(retVal) == 0:
445
        return listSa
446
    else:
447
        ## Remember pobyso_get_list_elements_so_so returns more information
448
        #  than just the elements of the list (# elements, is_ellipti)
449
        listSaSo, numElements, isEndElliptic = retVal
450
    if numElements == 0:
451
        return listSa
452
    ## Search first for the maximum precision of the elements
453
    maxPrecSa = 0
454
    for floatSo in listSaSo:
455
        pobyso_autoprint(floatSo)
456
        curPrecSa =  pobyso_get_prec_of_constant_so_sa(floatSo)
457
        if curPrecSa > maxPrecSa:
458
            maxPrecSa = curPrecSa
459
    ##
460
    RF = RealField(maxPrecSa)
461
    ##
462
    for floatSo in listSaSo:
463
        listSa.append(pobyso_get_constant_as_rn_with_rf_so_sa(floatSo))
464
    return listSa
465
# End pobyso_float_list_so_sa
466

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

    
515
def pobyso_float_poly_so_sa(polySo, realFieldSa=None):
516
    """
517
    Convert a Sollya polynomial into a Sage floating-point polynomial.
518
    If no realField is given, a RealField corresponding to the maximum 
519
    precision of the coefficients is internally computed.
520
    The real field is not returned but can be easily retrieved from 
521
    the polynomial itself.
522
    ALGORITHM:
523
    - (optional) compute the RealField of the coefficients;
524
    - convert the Sollya expression into a Sage expression;
525
    - convert the Sage expression into a Sage polynomial
526
    """    
527
    if realFieldSa is None:
528
        expressionPrecSa = pobyso_get_max_prec_of_exp_so_sa(polySo)
529
        #print "Maximum precision of Sollya polynomial coefficients:", expressionPrecSa
530
        if expressionPrecSa < 2 or expressionPrecSa > 2147483391:
531
            print "Maximum degree of expression:", expressionPrecSa
532
        realFieldSa      = RealField(expressionPrecSa)
533
    #print "Sollya expression before...",
534
    #pobyso_autoprint(polySo)
535

    
536
    expressionSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo,
537
                                                             realFieldSa)
538
    #print "...Sollya expression after."
539
    #pobyso_autoprint(polySo)
540
    polyVariableSa = expressionSa.variables()[0]
541
    polyRingSa     = realFieldSa[str(polyVariableSa)]
542
    #print polyRingSa
543
    # Do not use the polynomial(expressionSa, ring=polyRingSa) form!
544
    polynomialSa = polyRingSa(expressionSa)
545
    polyCoeffsListSa = polynomialSa.coefficients()
546
    #for coeff in polyCoeffsListSa:
547
    #    print coeff.abs().n()
548
    return polynomialSa
549
# End pobyso_float_poly_so_sa
550

    
551
def pobyso_free_variable():
552
    """
553
    Ultra thin wrapper around the sollya_lib_function_build_free_variable function.
554
    """
555
    return sollya_lib_build_function_free_variable()
556
    
557
def pobyso_function_type_as_string(funcType):
558
    """ Legacy function. See pobyso_function_type_as_string_so_sa. """
559
    return(pobyso_function_type_as_string_so_sa(funcType))
560

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

    
657
def pobyso_get_constant(rnArgSa, constSo):
658
    """ Legacy function. See pobyso_get_constant_so_sa. """
659
    return pobyso_get_constant_so_sa(rnArgSa, constSo)
660
# End pobyso_get_constant
661

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

    
704
def pobyso_get_constant_as_rn_with_rf(ctExp, realField):
705
    """ 
706
    Legacy function. See pobyso_get_constant_as_rn_with_rf_so_sa.
707
    """
708
    return pobyso_get_constant_as_rn_with_rf_so_sa(ctExp, realField)
709
# End pobyso_get_constant_as_rn_with_rf
710
    
711
def pobyso_get_constant_as_rn_with_rf_so_sa(ctExpSo, realFieldSa = None):
712
    """
713
    Get a Sollya constant as a Sage "real number".
714
    If no real field is specified, the precision of the floating-point number 
715
    returned is that of the Sollya constant.
716
    Otherwise is is that of the real field. Hence rounding may happen.
717
    """
718
    if realFieldSa is None:
719
        return pobyso_get_constant_as_rn_so_sa(ctExpSo)
720
    rnSa = realFieldSa(0)
721
    outcome = sollya_lib_get_constant(get_rn_value(rnSa), ctExpSo)
722
    if outcome == 0:
723
        return None
724
    else:
725
        return rnSa
726
# End pobyso_get_constant_as_rn_with_rf_so_sa
727

    
728
def pobyso_get_free_variable_name():
729
    """ 
730
    Legacy function. See pobyso_get_free_variable_name_so_sa.
731
    """
732
    return(pobyso_get_free_variable_name_so_sa())
733

    
734
def pobyso_get_free_variable_name_so_sa():
735
    return sollya_lib_get_free_variable_name()
736
    
737
def pobyso_get_function_arity(expressionSo):
738
    """ 
739
    Legacy function. See pobyso_get_function_arity_so_sa.
740
    """
741
    return(pobyso_get_function_arity_so_sa(expressionSo))
742

    
743
def pobyso_get_function_arity_so_sa(expressionSo):
744
    arity = c_int(0)
745
    sollya_lib_get_function_arity(byref(arity),expressionSo)
746
    return int(arity.value)
747

    
748
def pobyso_get_head_function(expressionSo):
749
    """ 
750
    Legacy function. See pobyso_get_head_function_so_sa. 
751
    """
752
    return(pobyso_get_head_function_so_sa(expressionSo)) 
753

    
754
def pobyso_get_head_function_so_sa(expressionSo):
755
    functionType = c_int(0)
756
    sollya_lib_get_head_function(byref(functionType), expressionSo)
757
    return int(functionType.value)
758

    
759
def pobyso_get_interval_from_range_so_sa(soRange, realIntervalFieldSa = None ):
760
    """
761
    Return the Sage interval corresponding to the Sollya range argument.
762
    If no reaIntervalField is passed as an argument, the interval bounds are not
763
    rounded: they are elements of RealIntervalField of the "right" precision
764
    to hold all the digits.
765
    """
766
    prec = c_int(0)
767
    if realIntervalFieldSa is None:
768
        retval = sollya_lib_get_prec_of_range(byref(prec), soRange, None)
769
        if retval == 0:
770
            return None
771
        realIntervalFieldSa = RealIntervalField(prec.value)
772
    intervalSa = realIntervalFieldSa(0,0)
773
    retval = \
774
        sollya_lib_get_interval_from_range(get_interval_value(intervalSa),\
775
                                           soRange)
776
    if retval == 0:
777
        return None
778
    return intervalSa
779
# End pobyso_get_interval_from_range_so_sa
780

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

    
824
def pobyso_get_max_prec_of_exp(soExp):
825
    """ Legacy function. See pobyso_get_max_prec_of_exp_so_sa. """
826
    return pobyso_get_max_prec_of_exp_so_sa(soExp)
827

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

    
869
def pobyso_get_min_prec_of_constant_so_sa(constExpSo):
870
    """
871
    Get the minimum precision necessary to represent the value of a Sollya
872
    constant.
873
    MPFR_MIN_PREC and powers of 2 are taken into account.
874
    We assume that constExpSo is a pointer to a Sollay constant expression.
875
    """
876
    constExpAsRnSa = pobyso_get_constant_as_rn_so_sa(constExpSo)
877
    return(min_mpfr_size(get_rn_value(constExpAsRnSa)))
878

    
879
def pobyso_get_poly_so_sa(polySo, realFieldSa=None):
880
    """
881
    Convert a Sollya polynomial into a Sage polynomial.
882
    Legacy function. Use pobyso_float_poly_so_sa() instead.
883
    """    
884
    return pobyso_float_poly_so_sa(polySo,realFieldSa)
885
# End pobyso_get_poly_so_sa
886

    
887
def pobyso_get_prec():
888
    """ Legacy function. See pobyso_get_prec_so_sa(). """
889
    return pobyso_get_prec_so_sa()
890

    
891
def pobyso_get_prec_so():
892
    """
893
    Get the current default precision in Sollya.
894
    The return value is a Sollya object.
895
    Usefull when modifying the precision back and forth by avoiding
896
    extra conversions.
897
    """
898
    return sollya_lib_get_prec(None)
899
    
900
def pobyso_get_prec_so_sa():
901
    """
902
    Get the current default precision in Sollya.
903
    The return value is Sage/Python int.
904
    """
905
    precSo = sollya_lib_get_prec()
906
    precSa = pobyso_constant_from_int_so_sa(precSo)
907
    sollya_lib_clear_obj(precSo)
908
    return precSa
909
# End pobyso_get_prec_so_sa.
910

    
911
def pobyso_get_prec_so_so_sa():
912
    """
913
    Return the current precision both as a Sollya object and a
914
    Sage integer as hybrid tuple.
915
    To avoid multiple calls for precision manipulations. 
916
    """
917
    precSo = sollya_lib_get_prec()
918
    precSa = pobyso_constant_from_int_so_sa(precSo)
919
    return (precSo, int(precSa))
920
    
921
def pobyso_get_prec_of_constant(ctExpSo):
922
    """ Legacy function. See pobyso_get_prec_of_constant_so_sa. """
923
    return pobyso_get_prec_of_constant_so_sa(ctExpSo)
924

    
925
def pobyso_get_prec_of_constant_so_sa(ctExpSo):
926
    """
927
    Tries to find a precision to represent ctExpSo without rounding.
928
    If not possible, returns None.
929
    """
930
    #print "Entering pobyso_get_prec_of_constant_so_sa..."
931
    prec = c_int(0)
932
    retc = sollya_lib_get_prec_of_constant(byref(prec), ctExpSo, None)
933
    if retc == 0:
934
        #print "pobyso_get_prec_of_constant_so_sa failed."
935
        return None
936
    #print "...exiting pobyso_get_prec_of_constant_so_sa."
937
    return int(prec.value)
938

    
939
def pobyso_get_prec_of_range_so_sa(rangeSo):
940
    """
941
    Returns the number of bits elements of a range are coded with.
942
    """
943
    prec = c_int(0)
944
    retc = sollya_lib_get_prec_of_range(byref(prec), rangeSo, None)
945
    if retc == 0:
946
        return(None)
947
    return int(prec.value)
948
# End pobyso_get_prec_of_range_so_sa()
949

    
950
def pobyso_get_sage_exp_from_sollya_exp(sollyaExpSo, realField = RR):
951
    """ Legacy function. See pobyso_get_sage_exp_from_sollya_exp_so_sa. """
952
    return pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, 
953
                                                     realField = RR)
954

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

    
1003

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

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

    
1153
def pobyso_inf_so_so(intervalSo):
1154
    """
1155
    Very thin wrapper around sollya_lib_inf().
1156
    """
1157
    return sollya_lib_inf(intervalSo)
1158
# End pobyso_inf_so_so.
1159
    
1160
def pobyso_infnorm_so_so(func, interval, file = None, intervalList = None):
1161
    print "Do not use this function. User pobyso_supnorm_so_so instead."
1162
    return None
1163

    
1164
def pobyso_interval_to_range_sa_so(intervalSa, precisionSa=None):
1165
    if precisionSa is None:
1166
        precisionSa = intervalSa.parent().precision()
1167
    intervalSo = pobyso_bounds_to_range_sa_so(intervalSa.lower(),\
1168
                                              intervalSa.upper(),\
1169
                                              precisionSa)
1170
    return intervalSo
1171
# End pobyso_interval_to_range_sa_so
1172

    
1173
def pobyso_is_error_so_sa(objSo):
1174
    """
1175
    Thin wrapper around the sollya_lib_obj_is_error() function.
1176
    """
1177
    if sollya_lib_obj_is_error(objSo) != 0:
1178
        return True
1179
    else:
1180
        return False
1181
# End pobyso_is_error-so_sa
1182

    
1183
def pobyso_is_floating_point_number_sa_sa(numberSa):
1184
    """
1185
    Check whether a Sage number is floating point.
1186
    Exception stuff added because numbers other than
1187
    floating-point ones do not have the is_real() attribute.
1188
    """
1189
    try:
1190
        return numberSa.is_real()
1191
    except AttributeError:
1192
        return False
1193
# End pobyso_is_floating_piont_number_sa_sa
1194

    
1195
def pobyso_is_range_so_sa(rangeCandidateSo):
1196
    """
1197
    Thin wrapper over sollya_lib_is_range.
1198
    """
1199
    return sollya_lib_obj_is_range(rangeCandidateSo) != 0
1200
    
1201
# End pobyso_is_range_so_sa
1202

    
1203

    
1204
def pobyso_lib_init():
1205
    sollya_lib_init(None)
1206

    
1207
def pobyso_lib_close():
1208
    sollya_lib_close(None)
1209

    
1210
def pobyso_name_free_variable(freeVariableNameSa):
1211
    """ Legacy function. See pobyso_name_free_variable_sa_so. """
1212
    pobyso_name_free_variable_sa_so(freeVariableNameSa)
1213

    
1214
def pobyso_name_free_variable_sa_so(freeVariableNameSa):
1215
    """
1216
    Set the free variable name in Sollya from a Sage string.
1217
    """
1218
    sollya_lib_name_free_variable(freeVariableNameSa)
1219

    
1220
def pobyso_parse_string(string):
1221
    """ Legacy function. See pobyso_parse_string_sa_so. """
1222
    return pobyso_parse_string_sa_so(string)
1223
 
1224
def pobyso_parse_string_sa_so(string):
1225
    """
1226
    Get the Sollya expression computed from a Sage string or
1227
    a Sollya error object if parsing failed.
1228
    """
1229
    return sollya_lib_parse_string(string)
1230

    
1231
def pobyso_precision_so_sa(ctExpSo):
1232
    """
1233
    Computes the necessary precision to represent a number.
1234
    If x is not zero, it can be uniquely written as x = m · 2e 
1235
    where m is an odd integer and e is an integer. 
1236
    precision(x) returns the number of bits necessary to write m 
1237
    in binary (i.e. ceil(log2(m))).
1238
    """
1239
    #TODO: take care of the special case: 0, @NaN@, @Inf@
1240
    precisionSo = sollya_lib_precision(ctExpSo)
1241
    precisionSa = pobyso_constant_from_int_so_sa(precisionSo)
1242
    sollya_lib_clear_obj(precisionSo)
1243
    return precisionSa
1244
# End pobyso_precision_so_sa
1245

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

    
1412
def pobyso_polynomial_degree_so_sa(polySo):
1413
    """
1414
    Return the degree of a Sollya polynomial as a Sage int.
1415
    """
1416
    degreeSo = sollya_lib_degree(polySo)
1417
    return pobyso_constant_from_int_so_sa(degreeSo)
1418
# End pobyso_polynomial_degree_so_sa
1419

    
1420
def pobyso_polynomial_degree_so_so(polySo):
1421
    """
1422
    Thin wrapper around lib_sollya_degree().
1423
    """
1424
    return sollya_lib_degree(polySo)
1425
# End pobyso_polynomial_degree_so_so
1426
                                                                  
1427
def pobyso_range(rnLowerBound, rnUpperBound):
1428
    """ Legacy function. See pobyso_range_sa_so. """
1429
    return pobyso_range_sa_so(rnLowerBound, rnUpperBound) 
1430

    
1431
def pobyso_range_from_bounds_sa_so(rnLowerBound, rnUpperBound, precSa = None):
1432
    """
1433
    Create a Sollya range from 2 Sage real numbers as bounds
1434
    """
1435
    # TODO check precision stuff.
1436
    sollyaPrecChanged = False
1437
    (initialSollyaPrecSo, initialSollyaPrecSa) = \
1438
        pobyso_get_prec_so_so_sa()
1439
    if precSa is None:
1440
        precSa = max(rnLowerBound.parent().prec(), rnUpperBound.parent().prec())
1441
    if precSa > initialSollyaPrecSa:
1442
        if precSa <= 2:
1443
            print inspect.stack()[0][3], ": precision change <= 2 requested."
1444
        pobyso_set_prec_sa_so(precSa)
1445
        sollyaPrecChanged = True
1446
    rangeSo = sollya_lib_range_from_bounds(get_rn_value(rnLowerBound),
1447
                                           get_rn_value(rnUpperBound))
1448
    if sollyaPrecChanged:
1449
        pobyso_set_prec_so_so(initialSollyaPrecSo)
1450
        sollya_lib_clear_obj(initialSollyaPrecSo)
1451
    return rangeSo
1452
# End pobyso_range_from_bounds_sa_so
1453

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

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

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

    
1631
def pobyso_remez_canonical_so_so(funcSo, \
1632
                                 degreeSo, \
1633
                                 rangeSo, \
1634
                                 weightSo = pobyso_constant_1_sa_so(),\
1635
                                 qualitySo = None):
1636
    """
1637
    All arguments are pointers to Sollya objects.
1638
    The return value is a pointer to a Sollya function.
1639
    """
1640
    if not sollya_lib_obj_is_function(funcSo):
1641
        return(None)
1642
    return(sollya_lib_remez(funcSo, degreeSo, rangeSo, weightSo, qualitySo, None))
1643
# End pobyso_remez_canonical_so_so.
1644

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

    
1714

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

    
1843
def pobyso_set_canonical_off():
1844
    sollya_lib_set_canonical(sollya_lib_off())
1845

    
1846
def pobyso_set_canonical_on():
1847
    sollya_lib_set_canonical(sollya_lib_on())
1848

    
1849
def pobyso_set_prec(p):
1850
    """ Legacy function. See pobyso_set_prec_sa_so. """
1851
    pobyso_set_prec_sa_so(p)
1852

    
1853
def pobyso_set_prec_sa_so(p):
1854
    #a = c_int(p)
1855
    #precSo = c_void_p(sollya_lib_constant_from_int(a))
1856
    #precSo = sollya_lib_constant_from_int(a)
1857
    precSo = pobyso_constant_from_int_sa_so(p)
1858
    sollya_lib_set_prec(precSo)
1859
    sollya_lib_clear_obj(precSo)
1860
# End pobyso_set_prec_sa_so.
1861

    
1862
def pobyso_set_prec_so_so(newPrecSo):
1863
    sollya_lib_set_prec(newPrecSo)
1864
# End pobyso_set_prec_so_so.
1865

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

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

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

    
2030
def pobyso_taylor(function, degree, point):
2031
    """ Legacy function. See pobysoTaylor_so_so. """
2032
    return(pobyso_taylor_so_so(function, degree, point))
2033

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

    
2126
def pobyso_taylorform_so_so(functionSo, degreeSo, pointSo, intervalSo=None, \
2127
                            errorTypeSo=None):
2128
    createdErrorType = False
2129
    if errorTypeSo is None:
2130
        errorTypeSo = sollya_lib_absolute()
2131
        createdErrorType = True
2132
    else:
2133
        #TODO: deal with the other case.
2134
        pass
2135
    if intervalSo is None:
2136
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
2137
                                         errorTypeSo, None)
2138
    else:
2139
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
2140
                                         intervalSo, errorTypeSo, None)
2141
    if createdErrorType:
2142
        sollya_lib_clear_obj(errorTypeSo)
2143
    return resultSo
2144
        
2145

    
2146
def pobyso_univar_polynomial_print_reverse(polySa):
2147
    """ Legacy function. See pobyso_univar_polynomial_print_reverse_sa_sa. """
2148
    return(pobyso_univar_polynomial_print_reverse_sa_sa(polySa))
2149

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