Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 228

Historique | Voir | Annoter | Télécharger (77,64 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 argumrny Sollya objects, create a Sollya end elliptic list.
136
    Elements of the list are "eaten" (should not be cleared individualy, 
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_change_var_in_function_so_so(funcSo, chvarExpSo):
159
    """
160
    Variable change in a function.
161
    """
162
    return(sollya_lib_evaluate(funcSo,chvarExpSo))
163
# End pobyso_change_var_in_function_so_so     
164

    
165
def pobyso_chebyshevform_so_so(functionSo, degreeSo, intervalSo):
166
    resultSo = sollya_lib_chebyshevform(functionSo, degreeSo, intervalSo)
167
    return(resultSo)
168
# End pobyso_chebyshevform_so_so.
169

    
170
def pobyso_clear_taylorform_sa_so(taylorFormSaSo):
171
    """
172
    This method is necessary to correctly clean up the memory from Taylor forms.
173
    These are made of a Sollya object, a Sollya object list, a Sollya object.
174
    For no clearly understood reason, sollya_lib_clear_object_list crashed 
175
    when applied to the object list.
176
    Here, we decompose it into Sage list of Sollya objects references and we
177
     clear them one by one. 
178
    """
179
    sollya_lib_clear_obj(taylorFormSaSo[0])
180
    (coefficientsErrorsListSaSo, numElementsSa, isEndEllipticSa) = \
181
        pobyso_get_list_elements_so_so(taylorFormSaSo[1])
182
    for element in coefficientsErrorsListSaSo:
183
        sollya_lib_clear_obj(element)
184
    sollya_lib_clear_obj(taylorFormSaSo[1])
185
    sollya_lib_clear_obj(taylorFormSaSo[2])
186
# End pobyso_clear_taylorform_sa_so 
187

    
188
def pobyso_cmp(rnArgSa, cteSo):
189
    """
190
    Compare the MPFR value a RealNumber with that of a Sollya constant.
191
    
192
    Get the value of the Sollya constant into a RealNumber and compare
193
    using MPFR. Could be optimized by working directly with a mpfr_t
194
    for the intermediate number. 
195
    """
196
    # Get the precision of the Sollya constant to build a Sage RealNumber 
197
    # with enough precision.to hold it.
198
    precisionOfCte = c_int(0)
199
    # From the Sollya constant, create a local Sage RealNumber.
200
    sollya_lib_get_prec_of_constant(precisionOfCte, cteSo) 
201
    #print "Precision of constant: ", precisionOfCte
202
    RRRR = RealField(precisionOfCte.value)
203
    rnLocalSa = RRRR(0)
204
    sollya_lib_get_constant(get_rn_value(rnLocalSa), cteSo)
205
    #
206
    ## Compare the Sage RealNumber version of the Sollya constant with rnArg.
207
    return(cmp_rn_value(rnArgSa, rnLocal))
208
# End pobyso_smp
209

    
210
def pobyso_compute_pos_function_abs_val_bounds_sa_sa(funcSa, lowerBoundSa, \
211
                                                     upperBoundSa):
212
    """
213
    TODO: completely rework and test.
214
    """
215
    pobyso = pobyso_name_free_variable_sa_so(funcSa.variables()[0])
216
    funcSo = pobyso_parse_string(funcSa._assume_str().replace('_SAGE_VAR_', ''))
217
    rangeSo = pobyso_range_sa_so(lowerBoundSa, upperBoundSa)
218
    infnormSo = pobyso_infnorm_so_so(funcSo,rangeSo)
219
    # Sollya return the infnorm as an interval.
220
    fMaxSa = pobyso_get_interval_from_range_so_sa(infnormSo)
221
    # Get the top bound and compute the binade top limit.
222
    fMaxUpperBoundSa = fMaxSa.upper()
223
    binadeTopLimitSa = 2**ceil(fMaxUpperBoundSa.log2())
224
    # Put up together the function to use to compute the lower bound.
225
    funcAuxSo = pobyso_parse_string(str(binadeTopLimitSa) +  \
226
                                    '-(' + f._assume_str().replace('_SAGE_VAR_', '') + ')')
227
    pobyso_autoprint(funcAuxSo)
228
    # Clear the Sollya range before a new call to infnorm and issue the call.
229
    sollya_lib_clear_obj(infnormSo)
230
    infnormSo = pobyso_infnorm_so_so(funcAuxSo,rangeSo)
231
    fMinSa = pobyso_get_interval_from_range_so_sa(infnormSo)
232
    sollya_lib_clear_obj(infnormSo)
233
    fMinLowerBoundSa = binadeTopLimitSa - fMinSa.lower()
234
    # Compute the maximum of the precisions of the different bounds.
235
    maxPrecSa = max([fMinLowerBoundSa.parent().precision(), \
236
                     fMaxUpperBoundSa.parent().precision()])
237
    # Create a RealIntervalField and create an interval with the "good" bounds.
238
    RRRI = RealIntervalField(maxPrecSa)
239
    imageIntervalSa = RRRI(fMinLowerBoundSa, fMaxUpperBoundSa)
240
    # Free the unneeded Sollya objects
241
    sollya_lib_clear_obj(funcSo)
242
    sollya_lib_clear_obj(funcAuxSo)
243
    sollya_lib_clear_obj(rangeSo)
244
    return(imageIntervalSa)
245
# End pobyso_compute_pos_function_abs_val_bounds_sa_sa
246

    
247
def pobyso_compute_precision_decay_ratio_function_sa_so():
248
    """
249
    Compute the precision decay ratio function for polynomial 
250
    coefficient progressive trucation.
251
    """
252
    functionText = """
253
    proc(deg, a, b, we, wq) 
254
    {
255
      k = we * (exp(x/a)-1) + wq * (b*x)^2 + (1-we-wq) * x;
256
      return k/k(d);
257
    };
258
    """
259
    return pobyso_parse_string_sa_so(functionText)
260
# End  pobyso_compute_precision_decay_ratio_function.
261

    
262

    
263
def pobyso_constant(rnArg):
264
    """ Legacy function. See pobyso_constant_sa_so. """
265
    return(pobyso_constant_sa_so(rnArg))
266
    
267
def pobyso_constant_sa_so(rnArgSa, precisionSa=None):
268
    """
269
    Create a Sollya constant from a Sage RealNumber.
270
    The sollya_lib_constant() function creates a constant
271
    with the same precision as the source.
272
    """
273
    ## Precision stuff. If one wants to change precisions,
274
    #  everything takes place in Sage. That only makes
275
    #  sense if one wants to reduce the precision.
276
    # TODO: revisit precision stuff with new technique.
277
    if not precisionSa is None:
278
        RRR = RealField(precisionSa)
279
        rnArgSa = RRR(rnArgSa)
280
    #print rnArgSa, rnArgSa.precision()
281
    # Sollya constant creation takes place here.
282
    return sollya_lib_constant(get_rn_value(rnArgSa))
283
# End pobyso_constant_sa_so
284
     
285
def pobyso_constant_0_sa_so():
286
    """
287
    Obvious.
288
    """
289
    return pobyso_constant_from_int_sa_so(0)
290

    
291
def pobyso_constant_1():
292
    """
293
    Obvious.
294
    Legacy function. See pobyso_constant_so_so. 
295
    """
296
    return pobyso_constant_1_sa_so()
297

    
298
def pobyso_constant_1_sa_so():
299
    """
300
    Obvious.
301
    """
302
    return(pobyso_constant_from_int_sa_so(1))
303

    
304
def pobyso_constant_from_int(anInt):
305
    """ Legacy function. See pobyso_constant_from_int_sa_so. """
306
    return pobyso_constant_from_int_sa_so(anInt)
307

    
308
def pobyso_constant_from_int_sa_so(anInt):
309
    """
310
    Get a Sollya constant from a Sage int.
311
    """
312
    return sollya_lib_constant_from_int64(long(anInt))
313

    
314
def pobyso_constant_from_int_so_sa(constSo):
315
    """
316
    Get a Sage int from a Sollya int constant.
317
    Usefull for precision or powers in polynomials.
318
    """
319
    constSa = c_long(0)
320
    sollya_lib_get_constant_as_int64(byref(constSa), constSo)
321
    return constSa.value
322
# End pobyso_constant_from_int_so_sa
323

    
324
def pobyso_constant_from_mpq_sa_so(rationalSa):
325
    """
326
    Make a Sollya constant from Sage rational.
327
    The Sollya constant is an unevaluated expression.
328
    Hence no precision argument is needed.
329
    It is better to leave this way since Sollya has its own
330
    optimized evaluation mecanism that tries very hard to
331
    return exact values or at least faithful ones.
332
    """
333
    ratExprSo = \
334
        sollya_lib_constant_from_mpq(sgmp_get_rational_value(rationalSa))
335
    return ratExprSo
336
# End pobyso_constant_from_mpq_sa_so.
337

    
338
def pobyso_constant_sollya_prec_sa_so(rnArgSa):
339
    """
340
    Create a Sollya constant from a Sage RealNumber at the
341
    current precision in Sollya.
342
    """
343
    currentSollyaPrecSa = pobyso_get_prec_so_sa()
344
    return pobyso_constant_sa_so(rnArgSa, currentSollyaPrecSa)
345
# End pobyso_constant_sollya_prec_sa_so
346

    
347
def pobyso_end_elliptic_list_so_sa_so(objectsListSo, intCountSa):
348
    """
349
    Create a Sollya end elliptic list made of the objectListSo[0] to
350
     objectsListSo[intCountSa-1] objects.
351
    """     
352
    return sollya_lib_end_elliptic_list(objectSo, int(intCountSa))
353

    
354
def pobyso_error_so():
355
    return sollya_lib_error(None)
356
# End pobyso_error().
357

    
358
def pobyso_evaluate_so_so(funcSo, argumentSo):
359
    """
360
    Evaluates funcSo for arguemntSo through sollya_lib_evaluate().
361
    """
362
    return sollya_lib_evaluate(funcSo, argumentSo)
363
# End pobyso_evaluate_so_so.
364

    
365
def pobyso_float_poly_sa_so(polySa, precSa = None):
366
    """
367
    Create a Sollya polynomial from a Sage RealField polynomial.
368
    """
369
    ## TODO: filter arguments.
370
    ## Precision. If a precision is given, convert the polynomial
371
    #  into the right polynomial field. If not convert it straight
372
    #  to Sollya.
373
    sollyaPrecChanged = False 
374
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
375
    if precSa is None:
376
        precSa = polySa.parent().base_ring().precision()
377
    if (precSa > initialSollyaPrecSa):
378
        assert precSa >= 2, "Precision change <2 requested"
379
        if precSa <= 2:
380
            print inspect.stack()[0][3], ": precision change <= 2 requested"
381
        precSo = pobyso_constant_from_int(precSa)
382
        pobyso_set_prec_so_so(precSo)
383
        sollya_lib_clear_obj(precSo)
384
        sollyaPrecChanged = True    
385
    ## Get exponents and coefficients.
386
    exponentsSa     = polySa.exponents()
387
    coefficientsSa  = polySa.coefficients()
388
    ## Build the polynomial.
389
    polySo = None
390
    for coefficientSa, exponentSa in zip(coefficientsSa, exponentsSa):
391
        #print coefficientSa.n(prec=precSa), exponentSa
392
        coefficientSo = \
393
            pobyso_constant_sa_so(coefficientSa)
394
        #pobyso_autoprint(coefficientSo)
395
        exponentSo = \
396
            pobyso_constant_from_int_sa_so(exponentSa)
397
        #pobyso_autoprint(exponentSo)
398
        monomialSo = sollya_lib_build_function_pow(
399
                       sollya_lib_build_function_free_variable(),
400
                       exponentSo)
401
        polyTermSo = sollya_lib_build_function_mul(coefficientSo,
402
                                                       monomialSo)
403
        if polySo is None:
404
            polySo = polyTermSo
405
        else:
406
            polySo = sollya_lib_build_function_add(polySo, polyTermSo)
407
    if sollyaPrecChanged:
408
        pobyso_set_prec_so_so(initialSollyaPrecSo)
409
        sollya_lib_clear_obj(initialSollyaPrecSo)
410
    return polySo
411
# End pobyso_float_poly_sa_so    
412

    
413
def pobyso_float_poly_so_sa(polySo, realFieldSa=None):
414
    """
415
    Convert a Sollya polynomial into a Sage floating-point polynomial.
416
    If no realField is given, a RealField corresponding to the maximum 
417
    precision of the coefficients is internally computed.
418
    The real field is not returned but can be easily retrieved from 
419
    the polynomial itself.
420
    ALGORITHM:
421
    - (optional) compute the RealField of the coefficients;
422
    - convert the Sollya expression into a Sage expression;
423
    - convert the Sage expression into a Sage polynomial
424
    """    
425
    if realFieldSa is None:
426
        expressionPrecSa = pobyso_get_max_prec_of_exp_so_sa(polySo)
427
        #print "Maximum precision of Sollya polynomial coefficients:", expressionPrecSa
428
        if expressionPrecSa < 2 or expressionPrecSa > 2147483391:
429
            print "Maximum degree of expression:", expressionPrecSa
430
        realFieldSa      = RealField(expressionPrecSa)
431
    #print "Sollya expression before...",
432
    #pobyso_autoprint(polySo)
433

    
434
    expressionSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo,
435
                                                             realFieldSa)
436
    #print "...Sollya expression after."
437
    #pobyso_autoprint(polySo)
438
    polyVariableSa = expressionSa.variables()[0]
439
    polyRingSa     = realFieldSa[str(polyVariableSa)]
440
    #print polyRingSa
441
    # Do not use the polynomial(expressionSa, ring=polyRingSa) form!
442
    polynomialSa = polyRingSa(expressionSa)
443
    polyCoeffsListSa = polynomialSa.coefficients()
444
    #for coeff in polyCoeffsListSa:
445
    #    print coeff.abs().n()
446
    return polynomialSa
447
# End pobyso_float_poly_so_sa
448

    
449
def pobyso_free_variable():
450
    """
451
    Ultra thin wrapper around the sollya_lib_function_build_free_variable function.
452
    """
453
    return sollya_lib_build_function_free_variable()
454
    
455
def pobyso_function_type_as_string(funcType):
456
    """ Legacy function. See pobyso_function_type_as_string_so_sa. """
457
    return(pobyso_function_type_as_string_so_sa(funcType))
458

    
459
def pobyso_function_type_as_string_so_sa(funcType):
460
    """
461
    Numeric Sollya function codes -> Sage mathematical function names.
462
    Notice that pow -> ^ (a la Sage, not a la Python).
463
    """
464
    if funcType == SOLLYA_BASE_FUNC_ABS:
465
        return "abs"
466
    elif funcType == SOLLYA_BASE_FUNC_ACOS:
467
        return "arccos"
468
    elif funcType == SOLLYA_BASE_FUNC_ACOSH:
469
        return "arccosh"
470
    elif funcType == SOLLYA_BASE_FUNC_ADD:
471
        return "+"
472
    elif funcType == SOLLYA_BASE_FUNC_ASIN:
473
        return "arcsin"
474
    elif funcType == SOLLYA_BASE_FUNC_ASINH:
475
        return "arcsinh"
476
    elif funcType == SOLLYA_BASE_FUNC_ATAN:
477
        return "arctan"
478
    elif funcType == SOLLYA_BASE_FUNC_ATANH:
479
        return "arctanh"
480
    elif funcType == SOLLYA_BASE_FUNC_CEIL:
481
        return "ceil"
482
    elif funcType == SOLLYA_BASE_FUNC_CONSTANT:
483
        return "cte"
484
    elif funcType == SOLLYA_BASE_FUNC_COS:
485
        return "cos"
486
    elif funcType == SOLLYA_BASE_FUNC_COSH:
487
        return "cosh"
488
    elif funcType == SOLLYA_BASE_FUNC_DIV:
489
        return "/"
490
    elif funcType == SOLLYA_BASE_FUNC_DOUBLE:
491
        return "double"
492
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEDOUBLE:
493
        return "doubleDouble"
494
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEEXTENDED:
495
        return "doubleDxtended"
496
    elif funcType == SOLLYA_BASE_FUNC_ERF:
497
        return "erf"
498
    elif funcType == SOLLYA_BASE_FUNC_ERFC:
499
        return "erfc"
500
    elif funcType == SOLLYA_BASE_FUNC_EXP:
501
        return "exp"
502
    elif funcType == SOLLYA_BASE_FUNC_EXP_M1:
503
        return "expm1"
504
    elif funcType == SOLLYA_BASE_FUNC_FLOOR:
505
        return "floor"
506
    elif funcType == SOLLYA_BASE_FUNC_FREE_VARIABLE:
507
        return "freeVariable"
508
    elif funcType == SOLLYA_BASE_FUNC_HALFPRECISION:
509
        return "halfPrecision"
510
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYCONSTANT:
511
        return "libraryConstant"
512
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYFUNCTION:
513
        return "libraryFunction"
514
    elif funcType == SOLLYA_BASE_FUNC_LOG:
515
        return "log"
516
    elif funcType == SOLLYA_BASE_FUNC_LOG_10:
517
        return "log10"
518
    elif funcType == SOLLYA_BASE_FUNC_LOG_1P:
519
        return "log1p"
520
    elif funcType == SOLLYA_BASE_FUNC_LOG_2:
521
        return "log2"
522
    elif funcType == SOLLYA_BASE_FUNC_MUL:
523
        return "*"
524
    elif funcType == SOLLYA_BASE_FUNC_NEARESTINT:
525
        return "round"
526
    elif funcType == SOLLYA_BASE_FUNC_NEG:
527
        return "__neg__"
528
    elif funcType == SOLLYA_BASE_FUNC_PI:
529
        return "pi"
530
    elif funcType == SOLLYA_BASE_FUNC_POW:
531
        return "^"
532
    elif funcType == SOLLYA_BASE_FUNC_PROCEDUREFUNCTION:
533
        return "procedureFunction"
534
    elif funcType == SOLLYA_BASE_FUNC_QUAD:
535
        return "quad"
536
    elif funcType == SOLLYA_BASE_FUNC_SIN:
537
        return "sin"
538
    elif funcType == SOLLYA_BASE_FUNC_SINGLE:
539
        return "single"
540
    elif funcType == SOLLYA_BASE_FUNC_SINH:
541
        return "sinh"
542
    elif funcType == SOLLYA_BASE_FUNC_SQRT:
543
        return "sqrt"
544
    elif funcType == SOLLYA_BASE_FUNC_SUB:
545
        return "-"
546
    elif funcType == SOLLYA_BASE_FUNC_TAN:
547
        return "tan"
548
    elif funcType == SOLLYA_BASE_FUNC_TANH:
549
        return "tanh"
550
    elif funcType == SOLLYA_BASE_FUNC_TRIPLEDOUBLE:
551
        return "tripleDouble"
552
    else:
553
        return None
554

    
555
def pobyso_get_constant(rnArgSa, constSo):
556
    """ Legacy function. See pobyso_get_constant_so_sa. """
557
    return pobyso_get_constant_so_sa(rnArgSa, constSo)
558
# End pobyso_get_constant
559

    
560
def pobyso_get_constant_so_sa(rnArgSa, constSo):
561
    """
562
    Set the value of rnArgSo to the value of constSo in MPFR_RNDN mode.
563
    rnArg must already exist and belong to some RealField.
564
    We assume that constSo points to a Sollya constant.
565
    """
566
    outcome = sollya_lib_get_constant(get_rn_value(rnArgSa), constSo)
567
    if outcome == 0: # Failure because constSo is not a constant expression.
568
        return None
569
    else:
570
        return outcome
571
# End  pobyso_get_constant_so_sa
572
   
573
def pobyso_get_constant_as_rn(ctExpSo):
574
    """ 
575
    Legacy function. See pobyso_get_constant_as_rn_so_sa. 
576
    """ 
577
    return(pobyso_get_constant_as_rn_so_sa(ctExpSo))
578
    
579
def pobyso_get_constant_as_rn_so_sa(constExpSo):
580
    """
581
    Get a Sollya constant as a Sage "real number".
582
    The precision of the floating-point number returned is that of the Sollya
583
    constant.
584
    """
585
    #print "Before computing precision of variable..."
586
    #pobyso_autoprint(constExpSo)
587
    precisionSa  = pobyso_get_prec_of_constant_so_sa(constExpSo)
588
    #print "precisionSa:", precisionSa
589
    ## If the expression can not be exactly converted, None is returned.
590
    #  In this case opt for the Sollya current expression.
591
    if precisionSa is None:
592
        precisionSa = pobyso_get_prec_so_sa()
593
    RRRR = RealField(precisionSa)
594
    rnSa = RRRR(0)
595
    outcome = sollya_lib_get_constant(get_rn_value(rnSa), constExpSo)
596
    if outcome == 0:
597
        return None
598
    else:
599
        return rnSa
600
# End pobyso_get_constant_as_rn_so_sa
601

    
602
def pobyso_get_constant_as_rn_with_rf(ctExp, realField):
603
    """ 
604
    Legacy function. See pobyso_get_constant_as_rn_with_rf_so_sa.
605
    """
606
    return pobyso_get_constant_as_rn_with_rf_so_sa(ctExp, realField)
607
# End pobyso_get_constant_as_rn_with_rf
608
    
609
def pobyso_get_constant_as_rn_with_rf_so_sa(ctExpSo, realFieldSa = None):
610
    """
611
    Get a Sollya constant as a Sage "real number".
612
    If no real field is specified, the precision of the floating-point number 
613
    returned is that of the Sollya constant.
614
    Otherwise is is that of the real field. Hence rounding may happen.
615
    """
616
    if realFieldSa is None:
617
        return pobyso_get_constant_as_rn_so_sa(ctExpSo)
618
    rnSa = realFieldSa(0)
619
    outcome = sollya_lib_get_constant(get_rn_value(rnSa), ctExpSo)
620
    if outcome == 0:
621
        return None
622
    else:
623
        return rnSa
624
# End pobyso_get_constant_as_rn_with_rf_so_sa
625

    
626
def pobyso_get_free_variable_name():
627
    """ 
628
    Legacy function. See pobyso_get_free_variable_name_so_sa.
629
    """
630
    return(pobyso_get_free_variable_name_so_sa())
631

    
632
def pobyso_get_free_variable_name_so_sa():
633
    return sollya_lib_get_free_variable_name()
634
    
635
def pobyso_get_function_arity(expressionSo):
636
    """ 
637
    Legacy function. See pobyso_get_function_arity_so_sa.
638
    """
639
    return(pobyso_get_function_arity_so_sa(expressionSo))
640

    
641
def pobyso_get_function_arity_so_sa(expressionSo):
642
    arity = c_int(0)
643
    sollya_lib_get_function_arity(byref(arity),expressionSo)
644
    return int(arity.value)
645

    
646
def pobyso_get_head_function(expressionSo):
647
    """ 
648
    Legacy function. See pobyso_get_head_function_so_sa. 
649
    """
650
    return(pobyso_get_head_function_so_sa(expressionSo)) 
651

    
652
def pobyso_get_head_function_so_sa(expressionSo):
653
    functionType = c_int(0)
654
    sollya_lib_get_head_function(byref(functionType), expressionSo)
655
    return int(functionType.value)
656

    
657
def pobyso_get_interval_from_range_so_sa(soRange, realIntervalFieldSa = None ):
658
    """
659
    Return the Sage interval corresponding to the Sollya range argument.
660
    If no reaIntervalField is passed as an argument, the interval bounds are not
661
    rounded: they are elements of RealIntervalField of the "right" precision
662
    to hold all the digits.
663
    """
664
    prec = c_int(0)
665
    if realIntervalFieldSa is None:
666
        retval = sollya_lib_get_prec_of_range(byref(prec), soRange, None)
667
        if retval == 0:
668
            return None
669
        realIntervalFieldSa = RealIntervalField(prec.value)
670
    intervalSa = realIntervalFieldSa(0,0)
671
    retval = \
672
        sollya_lib_get_interval_from_range(get_interval_value(intervalSa),\
673
                                           soRange)
674
    if retval == 0:
675
        return None
676
    return intervalSa
677
# End pobyso_get_interval_from_range_so_sa
678

    
679
def pobyso_get_list_elements(soObj):
680
    """ Legacy function. See pobyso_get_list_elements_so_so. """
681
    return pobyso_get_list_elements_so_so(soObj)
682
 
683
def pobyso_get_list_elements_so_so(objectListSo):
684
    """
685
    Get the Sollya list elements as a Sage/Python array of Sollya objects.
686
    
687
    INPUT:
688
    - objectListSo: a Sollya list of Sollya objects.
689
    
690
    OUTPUT:
691
    - a Sage/Python tuple made of:
692
      - a Sage/Python list of Sollya objects,
693
      - a Sage/Python int holding the number of elements,
694
      - a Sage/Python int stating (!= 0) that the list is end-elliptic.
695
    NOTE::
696
        We recover the addresses of the Sollya object from the list of pointers
697
        returned by sollya_lib_get_list_elements. The list itself is freed.
698
    TODO::
699
        Figure out what to do with numElements since the number of elements
700
        can easily be recovered from the list itself. 
701
        Ditto for isEndElliptic.
702
    """
703
    listAddress = POINTER(c_longlong)()
704
    numElements = c_int(0)
705
    isEndElliptic = c_int(0)
706
    listAsSageList = []
707
    result = sollya_lib_get_list_elements(byref(listAddress),\
708
                                          byref(numElements),\
709
                                          byref(isEndElliptic),\
710
                                          objectListSo)
711
    if result == 0 :
712
        return None
713
    for i in xrange(0, numElements.value, 1):
714
       #listAsSageList.append(sollya_lib_copy_obj(listAddress[i]))
715
       listAsSageList.append(listAddress[i])
716
       # Clear each of the elements returned by Sollya.
717
       #sollya_lib_clear_obj(listAddress[i])
718
    # Free the list itself.   
719
    sollya_lib_free(listAddress)
720
    return (listAsSageList, numElements.value, isEndElliptic.value)
721

    
722
def pobyso_get_max_prec_of_exp(soExp):
723
    """ Legacy function. See pobyso_get_max_prec_of_exp_so_sa. """
724
    return pobyso_get_max_prec_of_exp_so_sa(soExp)
725

    
726
def pobyso_get_max_prec_of_exp_so_sa(expSo):
727
    """
728
    Get the maximum precision used for the numbers in a Sollya expression.
729
    
730
    Arguments:
731
    soExp -- a Sollya expression pointer
732
    Return value:
733
    A Python integer
734
    TODO: 
735
    - error management;
736
    - correctly deal with numerical type such as DOUBLEEXTENDED.
737
    """
738
    if expSo is None:
739
        print inspect.stack()[0][3], ": expSo is None."
740
        return 0
741
    maxPrecision = 0
742
    minConstPrec = 0
743
    currentConstPrec = 0
744
    #pobyso_autoprint(expSo)
745
    operator = pobyso_get_head_function_so_sa(expSo)
746
    if (operator != SOLLYA_BASE_FUNC_CONSTANT) and \
747
    (operator != SOLLYA_BASE_FUNC_FREE_VARIABLE):
748
        (arity, subexpressions) = pobyso_get_subfunctions_so_sa(expSo)
749
        for i in xrange(arity):
750
            maxPrecisionCandidate = \
751
                pobyso_get_max_prec_of_exp_so_sa(subexpressions[i])
752
            if maxPrecisionCandidate > maxPrecision:
753
                maxPrecision = maxPrecisionCandidate
754
        return maxPrecision
755
    elif operator == SOLLYA_BASE_FUNC_CONSTANT:
756
        #minConstPrec = pobyso_get_min_prec_of_constant_so_sa(expSo)
757
        #currentConstPrec = pobyso_get_min_prec_of_constant_so_sa(soExp)
758
        #print minConstPrec, " - ", currentConstPrec 
759
        return pobyso_get_min_prec_of_constant_so_sa(expSo)
760
    
761
    elif operator == SOLLYA_BASE_FUNC_FREE_VARIABLE:
762
        return 0
763
    else:
764
        print "pobyso_get_max_prec_of_exp_so_sa: unexepected operator."
765
        return 0
766

    
767
def pobyso_get_min_prec_of_constant_so_sa(constExpSo):
768
    """
769
    Get the minimum precision necessary to represent the value of a Sollya
770
    constant.
771
    MPFR_MIN_PREC and powers of 2 are taken into account.
772
    We assume that constExpSo is a pointer to a Sollay constant expression.
773
    """
774
    constExpAsRnSa = pobyso_get_constant_as_rn_so_sa(constExpSo)
775
    return(min_mpfr_size(get_rn_value(constExpAsRnSa)))
776

    
777
def pobyso_get_poly_so_sa(polySo, realFieldSa=None):
778
    """
779
    Convert a Sollya polynomial into a Sage polynomial.
780
    Legacy function. Use pobyso_float_poly_so_sa() instead.
781
    """    
782
    return pobyso_float_poly_so_sa(polySo,realFieldSa)
783
# End pobyso_get_poly_so_sa
784

    
785
def pobyso_get_prec():
786
    """ Legacy function. See pobyso_get_prec_so_sa(). """
787
    return pobyso_get_prec_so_sa()
788

    
789
def pobyso_get_prec_so():
790
    """
791
    Get the current default precision in Sollya.
792
    The return value is a Sollya object.
793
    Usefull when modifying the precision back and forth by avoiding
794
    extra conversions.
795
    """
796
    return sollya_lib_get_prec(None)
797
    
798
def pobyso_get_prec_so_sa():
799
    """
800
    Get the current default precision in Sollya.
801
    The return value is Sage/Python int.
802
    """
803
    precSo = sollya_lib_get_prec()
804
    precSa = pobyso_constant_from_int_so_sa(precSo)
805
    sollya_lib_clear_obj(precSo)
806
    return precSa
807
# End pobyso_get_prec_so_sa.
808

    
809
def pobyso_get_prec_so_so_sa():
810
    """
811
    Return the current precision both as a Sollya object and a
812
    Sage integer as hybrid tuple.
813
    To avoid multiple calls for precision manipulations. 
814
    """
815
    precSo = sollya_lib_get_prec()
816
    precSa = pobyso_constant_from_int_so_sa(precSo)
817
    return (precSo, int(precSa))
818
    
819
def pobyso_get_prec_of_constant(ctExpSo):
820
    """ Legacy function. See pobyso_get_prec_of_constant_so_sa. """
821
    return pobyso_get_prec_of_constant_so_sa(ctExpSo)
822

    
823
def pobyso_get_prec_of_constant_so_sa(ctExpSo):
824
    """
825
    Tries to find a precision to represent ctExpSo without rounding.
826
    If not possible, returns None.
827
    """
828
    #print "Entering pobyso_get_prec_of_constant_so_sa..."
829
    prec = c_int(0)
830
    retc = sollya_lib_get_prec_of_constant(byref(prec), ctExpSo, None)
831
    if retc == 0:
832
        #print "pobyso_get_prec_of_constant_so_sa failed."
833
        return None
834
    #print "...exiting pobyso_get_prec_of_constant_so_sa."
835
    return int(prec.value)
836

    
837
def pobyso_get_prec_of_range_so_sa(rangeSo):
838
    """
839
    Returns the number of bits elements of a range are coded with.
840
    """
841
    prec = c_int(0)
842
    retc = sollya_lib_get_prec_of_range(byref(prec), rangeSo, None)
843
    if retc == 0:
844
        return(None)
845
    return int(prec.value)
846
# End pobyso_get_prec_of_range_so_sa()
847

    
848
def pobyso_get_sage_exp_from_sollya_exp(sollyaExpSo, realField = RR):
849
    """ Legacy function. See pobyso_get_sage_exp_from_sollya_exp_so_sa. """
850
    return pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, 
851
                                                     realField = RR)
852

    
853
def pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, realFieldSa = RR):
854
    """
855
    Get a Sage expression from a Sollya expression. 
856
    Currently only tested with polynomials with floating-point coefficients.
857
    Notice that, in the returned polynomial, the exponents are RealNumbers.
858
    """
859
    #pobyso_autoprint(sollyaExp)
860
    operatorSa = pobyso_get_head_function_so_sa(sollyaExpSo)
861
    sollyaLibFreeVariableName = sollya_lib_get_free_variable_name()
862
    ## Get rid of the "_"'s in "_x_", if any.
863
    sollyaLibFreeVariableName = re.sub('_', '', sollyaLibFreeVariableName)
864
    # Constants and the free variable are special cases.
865
    # All other operator are dealt with in the same way.
866
    if (operatorSa != SOLLYA_BASE_FUNC_CONSTANT) and \
867
       (operatorSa != SOLLYA_BASE_FUNC_FREE_VARIABLE):
868
        (aritySa, subexpressionsSa) = pobyso_get_subfunctions_so_sa(sollyaExpSo)
869
        if aritySa == 1:
870
            sageExpSa = eval(pobyso_function_type_as_string_so_sa(operatorSa) + \
871
            "(" + pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], \
872
            realFieldSa) + ")")
873
        elif aritySa == 2:
874
            # We do not get through the preprocessor.
875
            # The "^" operator is then a special case.
876
            if operatorSa == SOLLYA_BASE_FUNC_POW:
877
                operatorAsStringSa = "**"
878
            else:
879
                operatorAsStringSa = \
880
                    pobyso_function_type_as_string_so_sa(operatorSa)
881
            sageExpSa = \
882
              eval("pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], realFieldSa)"\
883
              + " " + operatorAsStringSa + " " + \
884
                   "pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[1], realFieldSa)")
885
        # We do not know yet how to deal with arity >= 3 
886
        # (is there any in Sollya anyway?).
887
        else:
888
            sageExpSa = eval('None')
889
        return sageExpSa
890
    elif operatorSa == SOLLYA_BASE_FUNC_CONSTANT:
891
        #print "This is a constant"
892
        return pobyso_get_constant_as_rn_with_rf_so_sa(sollyaExpSo, realFieldSa)
893
    elif operatorSa == SOLLYA_BASE_FUNC_FREE_VARIABLE:
894
        #print "This is the free variable"
895
        return eval(sollyaLibFreeVariableName)
896
    else:
897
        print "Unexpected"
898
        return eval('None')
899
# End pobyso_get_sage_exp_from_sollya_exp_so_sa
900

    
901

    
902
def pobyso_get_subfunctions(expressionSo):
903
    """ Legacy function. See pobyso_get_subfunctions_so_sa. """
904
    return pobyso_get_subfunctions_so_sa(expressionSo) 
905
# End pobyso_get_subfunctions.
906
 
907
def pobyso_get_subfunctions_so_sa(expressionSo):
908
    """
909
    Get the subfunctions of an expression.
910
    Return the number of subfunctions and the list of subfunctions addresses.
911
    S.T.: Could not figure out another way than that ugly list of declarations
912
    to recover the addresses of the subfunctions.
913
    We limit ourselves to arity 8 functions. 
914
    """
915
    subf0 = c_int(0)
916
    subf1 = c_int(0)
917
    subf2 = c_int(0)
918
    subf3 = c_int(0)
919
    subf4 = c_int(0)
920
    subf5 = c_int(0)
921
    subf6 = c_int(0)
922
    subf7 = c_int(0)
923
    subf8 = c_int(0)
924
    arity = c_int(0)
925
    nullPtr = POINTER(c_int)()
926
    sollya_lib_get_subfunctions(expressionSo, byref(arity), \
927
      byref(subf0), byref(subf1), byref(subf2), byref(subf3), \
928
      byref(subf4), byref(subf5),\
929
      byref(subf6), byref(subf7), byref(subf8), nullPtr, None) 
930
#    byref(cast(subfunctions[0], POINTER(c_int))), \
931
#    byref(cast(subfunctions[0], POINTER(c_int))), \
932
#    byref(cast(subfunctions[2], POINTER(c_int))), \
933
#    byref(cast(subfunctions[3], POINTER(c_int))), \
934
#    byref(cast(subfunctions[4], POINTER(c_int))), \
935
#    byref(cast(subfunctions[5], POINTER(c_int))), \
936
#    byref(cast(subfunctions[6], POINTER(c_int))), \
937
#    byref(cast(subfunctions[7], POINTER(c_int))), \
938
#    byref(cast(subfunctions[8], POINTER(c_int))), nullPtr)
939
    subfunctions = [subf0, subf1, subf2, subf3, subf4, subf5, subf6, subf7, \
940
                    subf8]
941
    subs = []
942
    if arity.value > pobyso_max_arity:
943
        return(0,[])
944
    for i in xrange(arity.value):
945
        subs.append(int(subfunctions[i].value))
946
        #print subs[i]
947
    return (int(arity.value), subs)
948
# End pobyso_get_subfunctions_so_sa
949
    
950
def pobyso_guess_degree_sa_sa(functionSa, intervalSa, approxErrorSa, 
951
                              weightSa=None, degreeBoundSa=None):
952
    """
953
    Sa_sa variant of the solly_guessdegree function.
954
    Return 0 if something goes wrong.
955
    """
956
    functionAsStringSa = functionSa._assume_str().replace('_SAGE_VAR_', '')
957
    functionSo = pobyso_parse_string_sa_so(functionAsStringSa)
958
    if pobyso_is_error_so_sa(functionSo):
959
        sollya_lib_clear_obj(functionSo)
960
        return 0
961
    rangeSo = pobyso_interval_to_range_sa_so(intervalSa)
962
    # The approximation error is expected to be a floating point number.
963
    if pobyso_is_floating_point_number_sa_sa(approxErrorSa):
964
        approxErrorSo = pobyso_constant_sa_so(approxErrorSa)
965
    else:
966
        approxErrorSo = pobyso_constant_sa_so(RR(approxErrorSa))
967
    if not weightSa is None:
968
        weightAsStringSa = weightSa._assume_str().replace('_SAGE_VAR_', '')
969
        weightSo = pobyso_parse_string_sa_so(weightAsStringSa)
970
        if pobyso_is_error_so_sa(weightSo):
971
            sollya_lib_clear_obj(functionSo)
972
            sollya_lib_clear_obj(rangeSo)
973
            sollya_lib_clear_obj(approxErrorSo)   
974
            sollya_lib_clear_obj(weightSo)
975
            return 0   
976
    else:
977
        weightSo = None
978
    if not degreeBoundSa is None:
979
        degreeBoundSo = pobyso_constant_from_int_sa_so(degreeBoundSa)
980
    else:
981
        degreeBoundSo = None
982
    guessedDegreeSa = pobyso_guess_degree_so_sa(functionSo,
983
                                                rangeSo,
984
                                                approxErrorSo,
985
                                                weightSo,
986
                                                degreeBoundSo)
987
    sollya_lib_clear_obj(functionSo)
988
    sollya_lib_clear_obj(rangeSo)
989
    sollya_lib_clear_obj(approxErrorSo)
990
    if not weightSo is None:
991
        sollya_lib_clear_obj(weightSo)
992
    if not degreeBoundSo is None:
993
        sollya_lib_clear_obj(degreeBoundSo)
994
    return guessedDegreeSa
995
# End poyso_guess_degree_sa_sa
996

    
997
def pobyso_guess_degree_so_sa(functionSo, rangeSo, errorSo, weightSo=None, \
998
                              degreeBoundSo=None):
999
    """
1000
    Thin wrapper around the guessdegree function.
1001
    Nevertheless, some precision control stuff has been appended.
1002
    """
1003
    # Deal with Sollya internal precision issues: if it is too small, 
1004
    # compared with the error, increases it to about twice -log2(error).
1005
    errorSa = pobyso_get_constant_as_rn_with_rf_so_sa(errorSo)
1006
    log2ErrorSa = errorSa.log2()
1007
    if log2ErrorSa < 0:
1008
        neededPrecisionSa = int(2 * int(-log2ErrorSa) / 64) * 64
1009
    else:
1010
        neededPrecisionSa = int(2 * int(log2ErrorSa) / 64) * 64
1011
    #print "Needed precision:", neededPrecisionSa
1012
    sollyaPrecisionHasChanged = False
1013
    (initialPrecSo, initialPrecSa) = pobyso_get_prec_so_so_sa()
1014
    if neededPrecisionSa > initialPrecSa:
1015
        if neededPrecisionSa <= 2:
1016
            print inspect.stack()[0][3], ": precision change <= 2 requested."
1017
        pobyso_set_prec_sa_so(neededPrecisionSa)
1018
        sollyaPrecisionHasChanged = True
1019
    #print "Guessing degree..."
1020
    # weightSo and degreeBoundsSo are optional arguments.
1021
    # As declared, sollya_lib_guessdegree must take 5 arguments.
1022
    if weightSo is None:
1023
        degreeRangeSo = sollya_lib_guessdegree(functionSo, rangeSo, errorSo, 
1024
                                               0, 0, None)
1025
    elif degreeBoundSo is None:
1026
        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, \
1027
                                                errorSo, weightSo, 0, None)
1028
    else:
1029
        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, errorSo, \
1030
                                                weightSo, degreeBoundSo, None)
1031
    #print "...degree guess done."
1032
    # Restore internal precision, if applicable.
1033
    if sollyaPrecisionHasChanged:
1034
        pobyso_set_prec_so_so(initialPrecSo)
1035
        sollya_lib_clear_obj(initialPrecSo)
1036
    degreeIntervalSa = pobyso_range_to_interval_so_sa(degreeRangeSo)
1037
    sollya_lib_clear_obj(degreeRangeSo)
1038
    # When ok, both bounds match.
1039
    # When the degree bound is too low, the upper bound is the degree
1040
    # for which the error can be honored.
1041
    # When it really goes wrong, the upper bound is infinity. 
1042
    if degreeIntervalSa.lower() == degreeIntervalSa.upper():
1043
        return int(degreeIntervalSa.lower())
1044
    else:
1045
        if degreeIntervalSa.upper().is_infinity():
1046
            return None
1047
        else:
1048
            return int(degreeIntervalSa.upper())
1049
    # End pobyso_guess_degree_so_sa
1050

    
1051
def pobyso_inf_so_so(intervalSo):
1052
    """
1053
    Very thin wrapper around sollya_lib_inf().
1054
    """
1055
    return sollya_lib_inf(intervalSo)
1056
# End pobyso_inf_so_so.
1057
    
1058
def pobyso_infnorm_so_so(func, interval, file = None, intervalList = None):
1059
    print "Do not use this function. User pobyso_supnorm_so_so instead."
1060
    return None
1061

    
1062
def pobyso_interval_to_range_sa_so(intervalSa, precisionSa=None):
1063
    if precisionSa is None:
1064
        precisionSa = intervalSa.parent().precision()
1065
    intervalSo = pobyso_bounds_to_range_sa_so(intervalSa.lower(),\
1066
                                              intervalSa.upper(),\
1067
                                              precisionSa)
1068
    return intervalSo
1069
# End pobyso_interval_to_range_sa_so
1070

    
1071
def pobyso_is_error_so_sa(objSo):
1072
    """
1073
    Thin wrapper around the sollya_lib_obj_is_error() function.
1074
    """
1075
    if sollya_lib_obj_is_error(objSo) != 0:
1076
        return True
1077
    else:
1078
        return False
1079
# End pobyso_is_error-so_sa
1080

    
1081
def pobyso_is_floating_point_number_sa_sa(numberSa):
1082
    """
1083
    Check whether a Sage number is floating point.
1084
    Exception stuff added because numbers other than
1085
    floating-point ones do not have the is_real() attribute.
1086
    """
1087
    try:
1088
        return numberSa.is_real()
1089
    except AttributeError:
1090
        return False
1091
# End pobyso_is_floating_piont_number_sa_sa
1092

    
1093
def pobyso_is_range_so_sa(rangeCandidateSo):
1094
    """
1095
    Thin wrapper over sollya_lib_is_range.
1096
    """
1097
    return sollya_lib_obj_is_range(rangeCandidateSo) != 0
1098
    
1099
# End pobyso_is_range_so_sa
1100

    
1101

    
1102
def pobyso_lib_init():
1103
    sollya_lib_init(None)
1104

    
1105
def pobyso_lib_close():
1106
    sollya_lib_close(None)
1107
    
1108
def pobyso_name_free_variable(freeVariableNameSa):
1109
    """ Legacy function. See pobyso_name_free_variable_sa_so. """
1110
    pobyso_name_free_variable_sa_so(freeVariableNameSa)
1111

    
1112
def pobyso_name_free_variable_sa_so(freeVariableNameSa):
1113
    """
1114
    Set the free variable name in Sollya from a Sage string.
1115
    """
1116
    sollya_lib_name_free_variable(freeVariableNameSa)
1117

    
1118
def pobyso_parse_string(string):
1119
    """ Legacy function. See pobyso_parse_string_sa_so. """
1120
    return pobyso_parse_string_sa_so(string)
1121
 
1122
def pobyso_parse_string_sa_so(string):
1123
    """
1124
    Get the Sollya expression computed from a Sage string or
1125
    a Sollya error object if parsing failed.
1126
    """
1127
    return sollya_lib_parse_string(string)
1128

    
1129
def pobyso_precision_so_sa(ctExpSo):
1130
    """
1131
    Computes the necessary precision to represent a number.
1132
    If x is not zero, it can be uniquely written as x = m · 2e 
1133
    where m is an odd integer and e is an integer. 
1134
    precision(x) returns the number of bits necessary to write m 
1135
    in binary (i.e. ceil(log2(m))).
1136
    """
1137
    #TODO: take care of the special case: 0, @NaN@, @Inf@
1138
    precisionSo = sollya_lib_precision(ctExpSo)
1139
    precisionSa = pobyso_constant_from_int_so_sa(precisionSo)
1140
    sollya_lib_clear_obj(precisionSo)
1141
    return precisionSa
1142
# End pobyso_precision_so_sa
1143

    
1144
def pobyso_polynomial_coefficients_progressive_round_so_so(polySo,
1145
                                                           funcSo,
1146
                                                           icSo,
1147
                                                           intervalSo,
1148
                                                           itpSo,
1149
                                                           ftpSo,
1150
                                                           maxPrecSo,
1151
                                                           maxErrSo,
1152
                                                           debug=False):
1153
    if debug:
1154
        print "Input arguments:"
1155
        pobyso_autoprint(polySo)
1156
        pobyso_autoprint(funcSo)
1157
        pobyso_autoprint(icSo)
1158
        pobyso_autoprint(intervalSo)
1159
        pobyso_autoprint(itpSo)
1160
        pobyso_autoprint(ftpSo)
1161
        pobyso_autoprint(maxPrecSo)
1162
        pobyso_autoprint(maxErrSo)
1163
        print "________________"
1164
    
1165
    ## Higher order function see: 
1166
    #  http://effbot.org/pyfaq/how-do-you-make-a-higher-order-function-in-python.htm
1167
    def precision_decay_ratio_function(degreeSa):
1168
        def outer(x):
1169
            def inner(x):
1170
                we = 3/8
1171
                wq = 2/8
1172
                a  = 2.2
1173
                b  = 2 
1174
                return we*(exp(x/a)-1) +  wq*((b*x)**2) + (1-we-wq)*x
1175
            return  inner(x)/inner(degreeSa)
1176
        return outer
1177
    
1178
    #
1179
    degreeSa        = pobyso_polynomial_degree_so_sa(polySo)
1180
    ratio           = precision_decay_ratio_function(degreeSa)
1181
    itpSa           = pobyso_constant_from_int_so_sa(itpSo)
1182
    ftpSa           = pobyso_constant_from_int_so_sa(ftpSo)
1183
    maxPrecSa       = pobyso_constant_from_int_so_sa(maxPrecSo)
1184
    maxErrSa        = pobyso_get_constant_as_rn_so_sa(maxErrSo)
1185
    if debug:
1186
        print "degreeSa:", degreeSa
1187
        print "ratio:", ratio
1188
        print "itpsSa:", itpSa
1189
        print "ftpSa:", ftpSa
1190
        print "maxPrecSa:", maxPrecSa
1191
        print "maxErrSa:", maxErrSa
1192
    lastResPolySo   = None
1193
    lastInfNormSo   = None
1194
    #print "About to enter the while loop..."
1195
    while True: 
1196
        resPolySo   = pobyso_constant_0_sa_so()
1197
        pDeltaSa    = ftpSa - itpSa
1198
        for indexSa in reversed(xrange(0,degreeSa+1)):
1199
            #print "Index:", indexSa
1200
            indexSo = pobyso_constant_from_int_sa_so(indexSa)
1201
            #pobyso_autoprint(indexSo)
1202
            #print ratio(indexSa)
1203
            ctpSa = floor(ftpSa - (pDeltaSa * ratio(indexSa)))
1204
            ctpSo = pobyso_constant_from_int_sa_so(ctpSa)
1205
            if debug:
1206
                print "Index:", indexSa, " - Target precision:", 
1207
                pobyso_autoprint(ctpSo)
1208
            cmonSo  = \
1209
                sollya_lib_build_function_mul(sollya_lib_coeff(polySo, indexSo),
1210
                                      sollya_lib_build_function_pow( \
1211
                                          sollya_lib_build_function_free_variable(), \
1212
                                          indexSo))
1213
            #pobyso_autoprint(cmonSo)
1214
            cmonrSo = pobyso_round_coefficients_single_so_so(cmonSo, ctpSo)
1215
            sollya_lib_clear_obj(cmonSo)
1216
            #pobyso_autoprint(cmonrSo)
1217
            resPolySo = sollya_lib_build_function_add(resPolySo,
1218
                                                      cmonrSo)
1219
            #pobyso_autoprint(resPolySo)
1220
        # End for index
1221
        freeVarSo     = sollya_lib_build_function_free_variable()
1222
        changeVarSo   = sollya_lib_sub(freeVarSo, icSo)
1223
        resPolyCvSo   = sollya_lib_evaluate(resPolySo, changeVarSo)
1224
        errFuncSo = sollya_lib_build_function_sub(sollya_lib_copy_obj(funcSo), 
1225
                                                  resPolyCvSo)
1226
        infNormSo = sollya_lib_dirtyinfnorm(errFuncSo, intervalSo)
1227
        cerrSa    = pobyso_get_constant_as_rn_so_sa(infNormSo)
1228
        if debug:
1229
            print "Infnorm (Sollya):", 
1230
            pobyso_autoprint(infNormSo)
1231
        sollya_lib_clear_obj(errFuncSo)
1232
        #print "Infnorm  (Sage):", cerrSa
1233
        if (cerrSa > maxErrSa):
1234
            if debug:
1235
                print "Error is too large."
1236
            if lastResPolySo is None:
1237
                if debug:
1238
                    print "Enlarging prec."
1239
                ntpSa = floor(ftpSa + ftpSa/50)
1240
                ## Can't enlarge (numerical)
1241
                if ntpSa == ftpSa:
1242
                    sollya_lib_clear_obj(resPolySo)
1243
                    return None
1244
                ## Can't enlarge (not enough precision left)
1245
                if ntpSa > maxPrecSa:
1246
                    sollya_lib_clear_obj(resPolySo)
1247
                    return None
1248
                ftpSa = ntpSa
1249
                continue
1250
            ## One enlargement took place.
1251
            else:
1252
                if debug:
1253
                    print "Exit with the last before last polynomial."
1254
                    print "Precision of highest degree monomial:", itpSa
1255
                    print "Precision of constant term          :", ftpSa
1256
                sollya_lib_clear_obj(resPolySo)
1257
                sollya_lib_clear_obj(infNormSo)
1258
                return (lastResPolySo, lastInfNormSo)
1259
        # cerrSa <= maxErrSa: scrap more bits, possibly.
1260
        else:
1261
            if debug:
1262
                print "Error is too small"
1263
            if cerrSa <= (maxErrSa/2):
1264
                if debug: 
1265
                    print "Shrinking prec."
1266
                ntpSa = floor(ftpSa - ftpSa/50)
1267
                ## Can't shrink (numerical)
1268
                if ntpSa == ftpSa:
1269
                    if not lastResPolySo is None:
1270
                        sollya_lib_clear_obj(lastResPolySo)
1271
                    if not lastInfNormSo is None:
1272
                        sollya_lib_clear_obj(lastInfNormSo)
1273
                    if debug:
1274
                        print "Exit because can't shrink anymore (numerically)."
1275
                        print "Precision of highest degree monomial:", itpSa
1276
                        print "Precision of constant term          :", ftpSa
1277
                    return (resPolySo, infNormSo)
1278
                ## Can't shrink (not enough precision left)
1279
                if ntpSa <= itpSa:
1280
                    if not lastResPolySo is None:
1281
                        sollya_lib_clear_obj(lastResPolySo)
1282
                    if not lastInfNormSo is None:
1283
                        sollya_lib_clear_obj(lastInfNormSo)
1284
                        print "Exit because can't shrink anymore (no bits left)."
1285
                        print "Precision of highest degree monomial:", itpSa
1286
                        print "Precision of constant term          :", ftpSa
1287
                    return (resPolySo, infNormSo)
1288
                ftpSa = ntpSa
1289
                if not lastResPolySo is None:
1290
                    sollya_lib_clear_obj(lastResPolySo)
1291
                if not lastInfNormSo is None:
1292
                    sollya_lib_clear_obj(lastInfNormSo)
1293
                lastResPolySo = resPolySo
1294
                lastInfNormSo = infNormSo
1295
                continue
1296
            else: # Error is not that small, just return
1297
                if not lastResPolySo is None:
1298
                    sollya_lib_clear_obj(lastResPolySo)
1299
                if not lastInfNormSo is None:
1300
                    sollya_lib_clear_obj(lastInfNormSo)
1301
                if debug:
1302
                    print "Exit normally."
1303
                    print "Precision of highest degree monomial:", itpSa
1304
                    print "Precision of constant term          :", ftpSa
1305
                return (resPolySo, infNormSo)                
1306
    # End wile True
1307
    return None
1308
# End pobyso_polynomial_coefficients_progressive_truncate_so_so.
1309

    
1310
def pobyso_polynomial_degree_so_sa(polySo):
1311
    """
1312
    Return the degree of a Sollya polynomial as a Sage int.
1313
    """
1314
    degreeSo = sollya_lib_degree(polySo)
1315
    return pobyso_constant_from_int_so_sa(degreeSo)
1316
# End pobyso_polynomial_degree_so_sa
1317

    
1318
def pobyso_polynomial_degree_so_so(polySo):
1319
    """
1320
    Thin wrapper around lib_sollya_degree().
1321
    """
1322
    return sollya_lib_degree(polySo)
1323
# End pobyso_polynomial_degree_so_so
1324
                                                                  
1325
def pobyso_range(rnLowerBound, rnUpperBound):
1326
    """ Legacy function. See pobyso_range_sa_so. """
1327
    return pobyso_range_sa_so(rnLowerBound, rnUpperBound) 
1328

    
1329
def pobyso_range_from_bounds_sa_so(rnLowerBound, rnUpperBound, precSa = None):
1330
    """
1331
    Create a Sollya range from 2 Sage real numbers as bounds
1332
    """
1333
    # TODO check precision stuff.
1334
    sollyaPrecChanged = False
1335
    (initialSollyaPrecSo, initialSollyaPrecSa) = \
1336
        pobyso_get_prec_so_so_sa()
1337
    if precSa is None:
1338
        precSa = max(rnLowerBound.parent().prec(), rnUpperBound.parent().prec())
1339
    if precSa > initialSollyaPrecSa:
1340
        if precSa <= 2:
1341
            print inspect.stack()[0][3], ": precision change <= 2 requested."
1342
        pobyso_set_prec_sa_so(precSa)
1343
        sollyaPrecChanged = True
1344
    rangeSo = sollya_lib_range_from_bounds(get_rn_value(rnLowerBound),
1345
                                           get_rn_value(rnUpperBound))
1346
    if sollyaPrecChanged:
1347
        pobyso_set_prec_so_so(initialSollyaPrecSo)
1348
        pobyso_lib_clear_obj(initialSollyaPrecSo)
1349
    return rangeSo
1350
# End pobyso_range_from_bounds_sa_so
1351

    
1352
def pobyso_range_max_abs_so_so(rangeSo):
1353
    """
1354
    Return, as Sollya constant, the maximum absolute value of a Sollay range.
1355
    """
1356
    lowerBoundSo = sollya_lib_inf(rangeSo)
1357
    upperBoundSo = sollya_lib_sup(rangeSo)
1358
    #
1359
    lowerBoundSo = sollya_lib_build_function_abs(lowerBoundSo)
1360
    upperBoundSo = sollya_lib_build_function_abs(upperBoundSo)
1361
    #pobyso_autoprint(lowerBoundSo)
1362
    #pobyso_autoprint(upperBoundSo)
1363
    #
1364
    maxAbsSo = sollya_lib_max(lowerBoundSo, upperBoundSo, None)
1365
    #sollya_lib_clear_obj(lowerBoundSo)
1366
    #sollya_lib_clear_obj(upperBoundSo)
1367
    return maxAbsSo
1368
# End pobyso_range_max_abs_so_so
1369
 
1370
def pobyso_range_to_interval_so_sa(rangeSo, realIntervalFieldSa = None):
1371
    """
1372
    Get a Sage interval from a Sollya range.
1373
    If no realIntervalField is given as a parameter, the Sage interval
1374
    precision is that of the Sollya range.
1375
    Otherwise, the precision is that of the realIntervalField. In this case
1376
    rounding may happen.
1377
    """
1378
    if realIntervalFieldSa is None:
1379
        precSa = pobyso_get_prec_of_range_so_sa(rangeSo)
1380
        realIntervalFieldSa = RealIntervalField(precSa)
1381
    intervalSa = \
1382
        pobyso_get_interval_from_range_so_sa(rangeSo, realIntervalFieldSa) 
1383
    return intervalSa
1384
# End pobyso_range_to_interval_so_sa
1385

    
1386
def pobyso_rat_poly_sa_so(polySa, precSa = None):
1387
    """
1388
    Create a Sollya polynomial from a Sage rational polynomial.
1389
    """
1390
    ## TODO: filter arguments.
1391
    ## Precision. If no precision is given, use the current precision
1392
    #  of Sollya.
1393
    if precSa is None:
1394
        precSa =  pobyso_get_prec_so_sa()
1395
    #print "Precision:",  precSa
1396
    RRR = RealField(precSa)
1397
    ## Create a Sage polynomial in the "right" precision.
1398
    P_RRR = RRR[polySa.variables()[0]]
1399
    polyFloatSa = P_RRR(polySa)
1400
    ## Make sure no precision is provided: pobyso_float_poly_sa_so will
1401
    #  recover it all by itself and not make an extra conversion.
1402
    return pobyso_float_poly_sa_so(polyFloatSa)
1403
    
1404
# End pobyso_rat_poly_sa_so    
1405
    
1406
def pobyso_remez_canonical_sa_sa(func, \
1407
                                 degree, \
1408
                                 lowerBound, \
1409
                                 upperBound, \
1410
                                 weight = None, \
1411
                                 quality = None):
1412
    """
1413
    All arguments are Sage/Python.
1414
    The functions (func and weight) must be passed as expressions or strings.
1415
    Otherwise the function fails. 
1416
    The return value is a Sage polynomial.
1417
    """
1418
    var('zorglub')    # Dummy variable name for type check only. Type of 
1419
    # zorglub is "symbolic expression".
1420
    polySo = pobyso_remez_canonical_sa_so(func, \
1421
                                 degree, \
1422
                                 lowerBound, \
1423
                                 upperBound, \
1424
                                 weight, \
1425
                                 quality)
1426
    # String test
1427
    if parent(func) == parent("string"):
1428
        functionSa = eval(func)
1429
    # Expression test.
1430
    elif type(func) == type(zorglub):
1431
        functionSa = func
1432
    else:
1433
        return None
1434
    #
1435
    maxPrecision = 0
1436
    if polySo is None:
1437
        return(None)
1438
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
1439
    RRRRSa = RealField(maxPrecision)
1440
    polynomialRingSa = RRRRSa[functionSa.variables()[0]]
1441
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, RRRRSa)
1442
    polySa = polynomial(expSa, polynomialRingSa)
1443
    sollya_lib_clear_obj(polySo)
1444
    return(polySa)
1445
# End pobyso_remez_canonical_sa_sa
1446
    
1447
def pobyso_remez_canonical(func, \
1448
                           degree, \
1449
                           lowerBound, \
1450
                           upperBound, \
1451
                           weight = "1", \
1452
                           quality = None):
1453
    """ Legacy function. See pobyso_remez_canonical_sa_so. """
1454
    return(pobyso_remez_canonical_sa_so(func, \
1455
                                        degree, \
1456
                                        lowerBound, \
1457
                                        upperBound, \
1458
                                        weight, \
1459
                                        quality))
1460
# End pobyso_remez_canonical.
1461

    
1462
def pobyso_remez_canonical_sa_so(func, \
1463
                                 degree, \
1464
                                 lowerBound, \
1465
                                 upperBound, \
1466
                                 weight = None, \
1467
                                 quality = None):
1468
    """
1469
    All arguments are Sage/Python.
1470
    The functions (func and weight) must be passed as expressions or strings.
1471
    Otherwise the function fails. 
1472
    The return value is a pointer to a Sollya function.
1473
    """
1474
    var('zorglub')    # Dummy variable name for type check only. Type of
1475
    # zorglub is "symbolic expression".
1476
    currentVariableNameSa = None
1477
    # The func argument can be of different types (string, 
1478
    # symbolic expression...)
1479
    if parent(func) == parent("string"):
1480
        localFuncSa = eval(func)
1481
        if len(localFuncSa.variables()) > 0:
1482
            currentVariableNameSa = localFuncSa.variables()[0]
1483
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1484
            functionSo = \
1485
              sollya_lib_parse_string(localFuncSa._assume_str().replace('_SAGE_VAR_', ''))
1486
    # Expression test.
1487
    elif type(func) == type(zorglub):
1488
        # Until we are able to translate Sage expressions into Sollya 
1489
        # expressions : parse the string version.
1490
        if len(func.variables()) > 0:
1491
            currentVariableNameSa = func.variables()[0]
1492
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1493
            functionSo = \
1494
              sollya_lib_parse_string(func._assume_str().replace('_SAGE_VAR_', ''))
1495
    else:
1496
        return(None)
1497
    if weight is None: # No weight given -> 1.
1498
        weightSo = pobyso_constant_1_sa_so()
1499
    elif parent(weight) == parent("string"): # Weight given as string: parse it.
1500
        weightSo = sollya_lib_parse_string(func)
1501
    elif type(weight) == type(zorglub): # Weight given as symbolice expression.
1502
        functionSo = \
1503
          sollya_lib_parse_string_sa_so(weight._assume_str().replace('_SAGE_VAR_', ''))
1504
    else:
1505
        return(None)
1506
    degreeSo = pobyso_constant_from_int(degree)
1507
    rangeSo = pobyso_bounds_to_range_sa_so(lowerBound, upperBound)
1508
    if not quality is None:
1509
        qualitySo= pobyso_constant_sa_so(quality)
1510
    else:
1511
        qualitySo = None
1512
        
1513
    remezPolySo = sollya_lib_remez(functionSo, \
1514
                                   degreeSo, \
1515
                                   rangeSo, \
1516
                                   weightSo, \
1517
                                   qualitySo, \
1518
                                   None)
1519
    sollya_lib_clear_obj(functionSo)
1520
    sollya_lib_clear_obj(degreeSo)
1521
    sollya_lib_clear_obj(rangeSo)
1522
    sollya_lib_clear_obj(weightSo)
1523
    if not qualitySo is None:
1524
        sollya_lib_clear_obj(qualitySo)
1525
    return(remezPolySo)
1526
# End pobyso_remez_canonical_sa_so
1527

    
1528
def pobyso_remez_canonical_so_so(funcSo, \
1529
                                 degreeSo, \
1530
                                 rangeSo, \
1531
                                 weightSo = pobyso_constant_1_sa_so(),\
1532
                                 qualitySo = None):
1533
    """
1534
    All arguments are pointers to Sollya objects.
1535
    The return value is a pointer to a Sollya function.
1536
    """
1537
    if not sollya_lib_obj_is_function(funcSo):
1538
        return(None)
1539
    return(sollya_lib_remez(funcSo, degreeSo, rangeSo, weightSo, qualitySo, None))
1540
# End pobyso_remez_canonical_so_so.
1541

    
1542
def pobyso_round_coefficients_progressive_so_so(polySo, 
1543
                                                funcSo,
1544
                                                precSo,
1545
                                                intervalSo,
1546
                                                icSo,
1547
                                                currentApproxErrorSo,
1548
                                                approxAccurSo,
1549
                                                debug=False):
1550
    if debug:
1551
        print "Input arguments:"
1552
        print "Polynomial: ", ; pobyso_autoprint(polySo)
1553
        print "Function: ", ; pobyso_autoprint(funcSo)
1554
        print "Internal precision: ", ; pobyso_autoprint(precSo)
1555
        print "Interval: ", ; pobyso_autoprint(intervalSo)
1556
        print "Current approximation error: ", ; pobyso_autoprint(currentApproxErrorSo)
1557
        print "Requested approxiation error: ", ; pobyso_autoprint(approxAccurSo)
1558
        print "________________"
1559
    approxAccurSa        = pobyso_get_constant_as_rn_so_sa(approxAccurSo)
1560
    currentApproxErrorSa = pobyso_get_constant_as_rn_so_sa(currentApproxErrorSo)
1561
    ## If the current approximation error is too close to the target, there is
1562
    #  no possible gain.
1563
    if currentApproxErrorSa >= approxAccurSa / 2:
1564
        return (polySo, currentApproxErrorSo)
1565
    degreeSa             = pobyso_polynomial_degree_so_sa(polySo)
1566
    intervalSa           = pobyso_range_to_interval_so_sa(intervalSo)
1567
    
1568
    if debug:
1569
        print "degreeSa              :", degreeSa
1570
        print "intervalSa            :", intervalSa.str(style='brackets')
1571
        print "currentApproxErrorSa  :", currentApproxErrorSa 
1572
        print "approxAccurSa         :", approxAccurSa 
1573
    ### Start with a 0 value expression.
1574
    radiusSa = intervalSa.absolute_diameter() / 2
1575
    if debug:
1576
        print "log2(radius):", RR(radiusSa).log2()
1577
    iterIndex = 0
1578
    while True: 
1579
        resPolySo = pobyso_constant_0_sa_so()
1580
        roundedPolyApproxAccurSa = approxAccurSa / 2
1581
        currentRadiusPowerSa = 1 
1582
        for degree in xrange(0,degreeSa + 1):
1583
            #### At round 0, use the agressive formula. At round 1, run the
1584
            #    proved formula.
1585
            if iterIndex == 0:
1586
                roundingPowerSa = \
1587
                    floor(((currentRadiusPowerSa/roundedPolyApproxAccurSa)*(degree+1)).log2())
1588
            else:
1589
                roundingPowerSa = \
1590
                    floor(((currentRadiusPowerSa/roundedPolyApproxAccurSa)*(degreeSa+1)).log2())
1591
            ## Under extreme conditions the above formulas can evaluate under 2, which is the
1592
            #  minimal precision of an MPFR number.
1593
            if roundingPowerSa < 2:
1594
                roundingPowerSa = 2
1595
            if debug:
1596
                print "roundedPolyApproxAccurSa", roundedPolyApproxAccurSa
1597
                print "currentRadiusPowerSa", currentRadiusPowerSa
1598
                print "Current rounding exponent:", roundingPowerSa
1599
            currentRadiusPowerSa *= radiusSa
1600
            index1So = pobyso_constant_from_int_sa_so(degree)
1601
            index2So = pobyso_constant_from_int_sa_so(degree)
1602
            ### Create a monomial with:
1603
            #   - the coefficient in the initial monomial at the current degrree;
1604
            #   - the current exponent;
1605
            #   - the free variable.
1606
            cmonSo  = \
1607
                sollya_lib_build_function_mul(sollya_lib_coeff(polySo, index1So),
1608
                                              sollya_lib_build_function_pow( \
1609
                                                sollya_lib_build_function_free_variable(), \
1610
                                                index2So))
1611
            roundingPowerSo = pobyso_constant_from_int_sa_so(roundingPowerSa)
1612
            cmonrSo = pobyso_round_coefficients_single_so_so(cmonSo, roundingPowerSo)
1613
            sollya_lib_clear_obj(cmonSo)
1614
            ### Add to the result polynomial.
1615
            resPolySo = sollya_lib_build_function_add(resPolySo,
1616
                                                      cmonrSo)
1617
        # End for.
1618
        ### Check the new polynomial.
1619
        freeVarSo     = sollya_lib_build_function_free_variable()
1620
        changeVarSo   = sollya_lib_sub(freeVarSo, icSo)
1621
        resPolyCvSo   = sollya_lib_evaluate(resPolySo, changeVarSo)
1622
        errFuncSo     = sollya_lib_build_function_sub(sollya_lib_copy_obj(funcSo), 
1623
                                                      resPolyCvSo)
1624
        infNormSo = sollya_lib_dirtyinfnorm(errFuncSo, intervalSo)
1625
        ### This also clears resPolyCvSo.
1626
        sollya_lib_clear_obj(errFuncSo)
1627
        cerrSa    = pobyso_get_constant_as_rn_so_sa(infNormSo)
1628
        if debug:
1629
            print "Error of the new polynomial:", cerrSa
1630
        ### If at round 1, return the initial polynomial error. This should
1631
        #   never happen since the rounding algorithm is proved. But some 
1632
        #   circumstances may break it (e.g. internal precision of tools).
1633
        if cerrSa > approxAccurSa:
1634
            if iterIndex > 0: # Round 1 and beyond.
1635
                sollya_lib_clear_obj(resPolySo)
1636
                sollya_lib_clear_obj(infNormSo)
1637
                return (polySo, currentApproxErrorSo)
1638
            else: # Round 0, got round 1
1639
                sollya_lib_clear_obj(resPolySo)
1640
                sollya_lib_clear_obj(infNormSo)
1641
                iterIndex += 1
1642
                continue
1643
        ### If get here it is because cerrSa <= approxAccurSa
1644
        ### Approximation error of the new polynomial is acceptable.
1645
        return (resPolySo, infNormSo)
1646
    # End while True
1647
# End pobyso_round_coefficients_progressive_so_so
1648
    
1649
def pobyso_round_coefficients_single_so_so(polySo, precSo):
1650
    """
1651
    Create a rounded coefficients polynomial from polynomial argument to
1652
    the number of bits in size argument.
1653
    All coefficients are set to the same precision.
1654
    """
1655
    ## TODO: check arguments.
1656
    endEllipListSo = pobyso_build_end_elliptic_list_so_so(precSo)
1657
    polySo = sollya_lib_roundcoefficients(polySo, endEllipListSo, None)
1658
    sollya_lib_clear_obj(endEllipListSo)
1659
    #sollya_lib_clear_obj(endEllipListSo)
1660
    return polySo
1661
    
1662
# End pobyso_round_coefficients_single_so_so
1663

    
1664
def pobyso_set_canonical_off():
1665
    sollya_lib_set_canonical(sollya_lib_off())
1666

    
1667
def pobyso_set_canonical_on():
1668
    sollya_lib_set_canonical(sollya_lib_on())
1669

    
1670
def pobyso_set_prec(p):
1671
    """ Legacy function. See pobyso_set_prec_sa_so. """
1672
    pobyso_set_prec_sa_so(p)
1673

    
1674
def pobyso_set_prec_sa_so(p):
1675
    #a = c_int(p)
1676
    #precSo = c_void_p(sollya_lib_constant_from_int(a))
1677
    #precSo = sollya_lib_constant_from_int(a)
1678
    precSo = pobyso_constant_from_int_sa_so(p)
1679
    sollya_lib_set_prec(precSo)
1680
    sollya_lib_clear_obj(precSo)
1681
# End pobyso_set_prec_sa_so.
1682

    
1683
def pobyso_set_prec_so_so(newPrecSo):
1684
    sollya_lib_set_prec(newPrecSo)
1685
# End pobyso_set_prec_so_so.
1686

    
1687
def pobyso_inf_so_so(intervalSo):
1688
    """
1689
    Very thin wrapper around sollya_lib_inf().
1690
    """
1691
    return sollya_lib_inf(intervalSo)
1692
# End pobyso_inf_so_so.
1693
    
1694
def pobyso_supnorm_so_so(polySo, funcSo, intervalSo, errorTypeSo = None,\
1695
                         accuracySo = None):
1696
    """
1697
    Computes the supnorm of the approximation error between the given 
1698
    polynomial and function.
1699
    errorTypeSo defaults to "absolute".
1700
    accuracySo defaults to 2^(-40).
1701
    """
1702
    if errorTypeSo is None:
1703
        errorTypeSo = sollya_lib_absolute(None)
1704
        errorTypeIsNone = True
1705
    else:
1706
        errorTypeIsNone = False
1707
    #
1708
    if accuracySo is None:
1709
        # Notice the **: we are in Pythonland!
1710
        accuracySo = pobyso_constant_sa_so(RR(2**(-40)))
1711
        accuracyIsNone = True
1712
    else:
1713
        accuracyIsNone = False
1714
    pobyso_autoprint(accuracySo)
1715
    resultSo = \
1716
        sollya_lib_supnorm(polySo, funcSo, intervalSo, errorTypeSo, \
1717
                              accuracySo)
1718
    if errorTypeIsNone:
1719
        sollya_lib_clear_obj(errorTypeSo)
1720
    if accuracyIsNone:
1721
        sollya_lib_clear_obj(accuracySo)
1722
    return resultSo
1723
# End pobyso_supnorm_so_so
1724

    
1725
def pobyso_taylor_expansion_no_change_var_so_so(functionSo, 
1726
                                                degreeSo, 
1727
                                                rangeSo,
1728
                                                errorTypeSo=None, 
1729
                                                sollyaPrecSo=None):
1730
    """
1731
    Compute the Taylor expansion without the variable change
1732
    x -> x-intervalCenter.
1733
    """
1734
    # Change internal Sollya precision, if needed.
1735
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
1736
    sollyaPrecChanged = False
1737
    if sollyaPrecSo is None:
1738
        pass
1739
    else:
1740
        sollya_lib_set_prec(sollyaPrecSo)
1741
        sollyaPrecChanged = True
1742
    # Error type stuff: default to absolute.
1743
    if errorTypeSo is None:
1744
        errorTypeIsNone = True
1745
        errorTypeSo = sollya_lib_absolute(None)
1746
    else:
1747
        errorTypeIsNone = False
1748
    intervalCenterSo = sollya_lib_mid(rangeSo, None)
1749
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo,
1750
                                         intervalCenterSo,
1751
                                         rangeSo, errorTypeSo, None)
1752
    # taylorFormListSaSo is a Python list of Sollya objects references that 
1753
    # are copies of the elements of taylorFormSo.
1754
    # pobyso_get_list_elements_so_so clears taylorFormSo.
1755
    (taylorFormListSaSo, numElementsSa, isEndEllipticSa) = \
1756
        pobyso_get_list_elements_so_so(taylorFormSo)
1757
    polySo = sollya_lib_copy_obj(taylorFormListSaSo[0])
1758
    #print "Num elements:", numElementsSa
1759
    sollya_lib_clear_obj(taylorFormSo)
1760
    #polySo = taylorFormListSaSo[0]
1761
    #errorRangeSo = sollya_lib_copy_obj(taylorFormListSaSo[2])
1762
    errorRangeSo = taylorFormListSaSo[2]
1763
    # No copy_obj needed here: a new objects are created.
1764
    maxErrorSo    = sollya_lib_sup(errorRangeSo)
1765
    minErrorSo    = sollya_lib_inf(errorRangeSo)
1766
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
1767
    absMinErrorSo = sollya_lib_abs(minErrorSo)
1768
    sollya_lib_clear_obj(maxErrorSo)
1769
    sollya_lib_clear_obj(minErrorSo)
1770
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
1771
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
1772
    # If changed, reset the Sollya working precision.
1773
    if sollyaPrecChanged:
1774
        sollya_lib_set_prec(initialSollyaPrecSo)
1775
    sollya_lib_clear_obj(initialSollyaPrecSo)
1776
    if errorTypeIsNone:
1777
        sollya_lib_clear_obj(errorTypeSo)
1778
    pobyso_clear_taylorform_sa_so(taylorFormListSaSo)
1779
    if absMaxErrorSa > absMinErrorSa:
1780
        sollya_lib_clear_obj(absMinErrorSo)
1781
        return (polySo, intervalCenterSo, absMaxErrorSo)
1782
    else:
1783
        sollya_lib_clear_obj(absMaxErrorSo)
1784
        return (polySo, intervalCenterSo, absMinErrorSo)
1785
# end pobyso_taylor_expansion_no_change_var_so_so
1786

    
1787
def pobyso_taylor_expansion_with_change_var_so_so(functionSo, degreeSo, \
1788
                                                  rangeSo, \
1789
                                                  errorTypeSo=None, \
1790
                                                  sollyaPrecSo=None):
1791
    """
1792
    Compute the Taylor expansion with the variable change
1793
    x -> (x-intervalCenter) included.
1794
    """
1795
    # Change Sollya internal precision, if need.
1796
    sollyaPrecChanged = False
1797
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_sos_sa()
1798
    if sollyaPrecSo is None:
1799
        pass
1800
    else:
1801
        sollya_lib_set_prec(sollyaPrecSo)
1802
        sollyaPrecChanged = True
1803
    #
1804
    # Error type stuff: default to absolute.
1805
    if errorTypeSo is None:
1806
        errorTypeIsNone = True
1807
        errorTypeSo = sollya_lib_absolute(None)
1808
    else:
1809
        errorTypeIsNone = False
1810
    intervalCenterSo = sollya_lib_mid(rangeSo)
1811
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo, \
1812
                                         intervalCenterSo, \
1813
                                         rangeSo, errorTypeSo, None)
1814
    # taylorFormListSaSo is a Python list of Sollya objects references that 
1815
    # are copies of the elements of taylorFormSo.
1816
    # pobyso_get_list_elements_so_so clears taylorFormSo.
1817
    (taylorFormListSo, numElements, isEndElliptic) = \
1818
        pobyso_get_list_elements_so_so(taylorFormSo)
1819
    polySo = taylorFormListSo[0]
1820
    errorRangeSo = taylorFormListSo[2]
1821
    maxErrorSo    = sollya_lib_sup(errorRangeSo)
1822
    minErrorSo    = sollya_lib_inf(errorRangeSo)
1823
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
1824
    absMinErrorSo = sollya_lib_abs(minErrorSo)
1825
    sollya_lib_clear_obj(maxErrorSo)
1826
    sollya_lib_clear_obj(minErrorSo)
1827
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
1828
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
1829
    changeVarExpSo = sollya_lib_build_function_sub(\
1830
                       sollya_lib_build_function_free_variable(),\
1831
                       sollya_lib_copy_obj(intervalCenterSo))
1832
    polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpSo)
1833
    sollya_lib_clear_obj(polySo) 
1834
    sollya_lib_clear_obj(changeVarExpSo)
1835
    # If changed, reset the Sollya working precision.
1836
    if sollyaPrecChanged:
1837
        sollya_lib_set_prec(initialSollyaPrecSo)
1838
    sollya_lib_clear_obj(initialSollyaPrecSo)
1839
    if errorTypeIsNone:
1840
        sollya_lib_clear_obj(errorTypeSo)
1841
    sollya_lib_clear_obj(taylorFormSo)
1842
    # Do not clear maxErrorSo.
1843
    if absMaxErrorSa > absMinErrorSa:
1844
        sollya_lib_clear_obj(absMinErrorSo)
1845
        return((polyVarChangedSo, intervalCenterSo, absMaxErrorSo))
1846
    else:
1847
        sollya_lib_clear_obj(absMaxErrorSo)
1848
        return((polyVarChangedSo, intervalCenterSo, absMinErrorSo))
1849
# end pobyso_taylor_expansion_with_change_var_so_so
1850

    
1851
def pobyso_taylor(function, degree, point):
1852
    """ Legacy function. See pobysoTaylor_so_so. """
1853
    return(pobyso_taylor_so_so(function, degree, point))
1854

    
1855
def pobyso_taylor_so_so(functionSo, degreeSo, pointSo):
1856
    return(sollya_lib_taylor(functionSo, degreeSo, pointSo))
1857
    
1858
def pobyso_taylorform(function, degree, point = None, 
1859
                      interval = None, errorType=None):
1860
    """ Legacy function. See pobyso_taylorform_sa_sa;"""
1861
    
1862
def pobyso_taylorform_sa_sa(functionSa, \
1863
                            degreeSa, \
1864
                            pointSa, \
1865
                            intervalSa=None, \
1866
                            errorTypeSa=None, \
1867
                            precisionSa=None):
1868
    """
1869
    Compute the Taylor form of 'degreeSa' for 'functionSa' at 'pointSa' 
1870
    for 'intervalSa' with 'errorTypeSa' (a string) using 'precisionSa'. 
1871
    point: must be a Real or a Real interval.
1872
    return the Taylor form as an array
1873
    TODO: take care of the interval and of the point when it is an interval;
1874
          when errorType is not None;
1875
          take care of the other elements of the Taylor form (coefficients 
1876
          errors and delta.
1877
    """
1878
    # Absolute as the default error.
1879
    if errorTypeSa is None:
1880
        errorTypeSo = sollya_lib_absolute()
1881
    elif errorTypeSa == "relative":
1882
        errorTypeSo = sollya_lib_relative()
1883
    elif errortypeSa == "absolute":
1884
        errorTypeSo = sollya_lib_absolute()
1885
    else:
1886
        # No clean up needed.
1887
        return None
1888
    # Global precision stuff
1889
    sollyaPrecisionChangedSa = False
1890
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
1891
    if precisionSa is None:
1892
        precSa = initialSollyaPrecSa
1893
    else:
1894
        if precSa > initialSollyaPrecSa:
1895
            if precSa <= 2:
1896
                print inspect.stack()[0][3], ":precision change <= 2 requested."
1897
            pobyso_set_prec_sa_so(precSa)
1898
            sollyaPrecisionChangedSa = True
1899
    #        
1900
    if len(functionSa.variables()) > 0:
1901
        varSa = functionSa.variables()[0]
1902
        pobyso_name_free_variable_sa_so(str(varSa))
1903
    # In any case (point or interval) the parent of pointSa has a precision
1904
    # method.
1905
    pointPrecSa = pointSa.parent().precision()
1906
    if precSa > pointPrecSa:
1907
        pointPrecSa = precSa
1908
    # In any case (point or interval) pointSa has a base_ring() method.
1909
    pointBaseRingString = str(pointSa.base_ring())
1910
    if re.search('Interval', pointBaseRingString) is None: # Point
1911
        pointSo = pobyso_constant_sa_so(pointSa, pointPrecSa)
1912
    else: # Interval.
1913
        pointSo = pobyso_interval_to_range_sa_so(pointSa, pointPrecSa)
1914
    # Sollyafy the function.
1915
    functionSo = pobyso_parse_string_sa_so(functionSa._assume_str().replace('_SAGE_VAR_', ''))
1916
    if sollya_lib_obj_is_error(functionSo):
1917
        print "pobyso_tailorform: function string can't be parsed!"
1918
        return None
1919
    # Sollyafy the degree
1920
    degreeSo = sollya_lib_constant_from_int(int(degreeSa))
1921
    # Sollyafy the point
1922
    # Call Sollya
1923
    taylorFormSo = \
1924
        sollya_lib_taylorform(functionSo, degreeSo, pointSo, errorTypeSo,\
1925
                                         None)
1926
    sollya_lib_clear_obj(functionSo)
1927
    sollya_lib_clear_obj(degreeSo)
1928
    sollya_lib_clear_obj(pointSo)
1929
    sollya_lib_clear_obj(errorTypeSo)
1930
    (tfsAsList, numElements, isEndElliptic) = \
1931
            pobyso_get_list_elements_so_so(taylorFormSo)
1932
    polySo = tfsAsList[0]
1933
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
1934
    polyRealField = RealField(maxPrecision)
1935
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, polyRealField)
1936
    if sollyaPrecisionChangedSa:
1937
        sollya_lib_set_prec(initialSollyaPrecSo)
1938
    sollya_lib_clear_obj(initialSollyaPrecSo)
1939
    polynomialRing = polyRealField[str(varSa)]
1940
    polySa = polynomial(expSa, polynomialRing)
1941
    taylorFormSa = [polySa]
1942
    # Final clean-up.
1943
    sollya_lib_clear_obj(taylorFormSo)
1944
    return(taylorFormSa)
1945
# End pobyso_taylor_form_sa_sa
1946

    
1947
def pobyso_taylorform_so_so(functionSo, degreeSo, pointSo, intervalSo=None, \
1948
                            errorTypeSo=None):
1949
    createdErrorType = False
1950
    if errorTypeSo is None:
1951
        errorTypeSo = sollya_lib_absolute()
1952
        createdErrorType = True
1953
    else:
1954
        #TODO: deal with the other case.
1955
        pass
1956
    if intervalSo is None:
1957
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
1958
                                         errorTypeSo, None)
1959
    else:
1960
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
1961
                                         intervalSo, errorTypeSo, None)
1962
    if createdErrorType:
1963
        sollya_lib_clear_obj(errorTypeSo)
1964
    return resultSo
1965
        
1966

    
1967
def pobyso_univar_polynomial_print_reverse(polySa):
1968
    """ Legacy function. See pobyso_univar_polynomial_print_reverse_sa_sa. """
1969
    return(pobyso_univar_polynomial_print_reverse_sa_sa(polySa))
1970

    
1971
def pobyso_univar_polynomial_print_reverse_sa_sa(polySa):
1972
    """
1973
    Return the string representation of a univariate polynomial with
1974
    monomials ordered in the x^0..x^n order of the monomials.
1975
    Remember: Sage
1976
    """
1977
    polynomialRing = polySa.base_ring()
1978
    # A very expensive solution:
1979
    # -create a fake multivariate polynomial field with only one variable,
1980
    #   specifying a negative lexicographical order;
1981
    mpolynomialRing = PolynomialRing(polynomialRing.base(), \
1982
                                     polynomialRing.variable_name(), \
1983
                                     1, order='neglex')
1984
    # - convert the univariate argument polynomial into a multivariate
1985
    #   version;
1986
    p = mpolynomialRing(polySa)
1987
    # - return the string representation of the converted form.
1988
    # There is no simple str() method defined for p's class.
1989
    return(p.__str__())
1990
#
1991
print pobyso_get_prec()  
1992
pobyso_set_prec(165)
1993
print pobyso_get_prec()  
1994
a=100
1995
print type(a)
1996
id(a)
1997
print "Max arity: ", pobyso_max_arity
1998
print "Function tripleDouble (43) as a string: ", pobyso_function_type_as_string(43)
1999
print "Function None (44) as a string: ", pobyso_function_type_as_string(44)
2000
print "...Pobyso check done"