Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 227

Historique | Voir | Annoter | Télécharger (77,31 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
        if precSa <= 2:
379
            print inspect.stack()[0][3], ": precision change <= 2 requested"
380
        precSo = pobyso_constant_from_int(precSa)
381
        pobyso_set_prec_so_so(precSo)
382
        sollya_lib_clear_obj(precSo)
383
        sollyaPrecChanged = True    
384
    ## Get exponents and coefficients.
385
    exponentsSa     = polySa.exponents()
386
    coefficientsSa  = polySa.coefficients()
387
    ## Build the polynomial.
388
    polySo = None
389
    for coefficientSa, exponentSa in zip(coefficientsSa, exponentsSa):
390
        #print coefficientSa.n(prec=precSa), exponentSa
391
        coefficientSo = \
392
            pobyso_constant_sa_so(coefficientSa)
393
        #pobyso_autoprint(coefficientSo)
394
        exponentSo = \
395
            pobyso_constant_from_int_sa_so(exponentSa)
396
        #pobyso_autoprint(exponentSo)
397
        monomialSo = sollya_lib_build_function_pow(
398
                       sollya_lib_build_function_free_variable(),
399
                       exponentSo)
400
        polyTermSo = sollya_lib_build_function_mul(coefficientSo,
401
                                                       monomialSo)
402
        if polySo is None:
403
            polySo = polyTermSo
404
        else:
405
            polySo = sollya_lib_build_function_add(polySo, polyTermSo)
406
    if sollyaPrecChanged:
407
        pobyso_set_prec_so_so(initialSollyaPrecSo)
408
        sollya_lib_clear_obj(initialSollyaPrecSo)
409
    return polySo
410
# End pobyso_float_poly_sa_so    
411

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    
900

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

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

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

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

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

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

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

    
1100

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    
1659
def pobyso_set_canonical_off():
1660
    sollya_lib_set_canonical(sollya_lib_off())
1661

    
1662
def pobyso_set_canonical_on():
1663
    sollya_lib_set_canonical(sollya_lib_on())
1664

    
1665
def pobyso_set_prec(p):
1666
    """ Legacy function. See pobyso_set_prec_sa_so. """
1667
    pobyso_set_prec_sa_so(p)
1668

    
1669
def pobyso_set_prec_sa_so(p):
1670
    #a = c_int(p)
1671
    #precSo = c_void_p(sollya_lib_constant_from_int(a))
1672
    #precSo = sollya_lib_constant_from_int(a)
1673
    precSo = pobyso_constant_from_int_sa_so(p)
1674
    sollya_lib_set_prec(precSo)
1675
    sollya_lib_clear_obj(precSo)
1676
# End pobyso_set_prec_sa_so.
1677

    
1678
def pobyso_set_prec_so_so(newPrecSo):
1679
    sollya_lib_set_prec(newPrecSo)
1680
# End pobyso_set_prec_so_so.
1681

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

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

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

    
1846
def pobyso_taylor(function, degree, point):
1847
    """ Legacy function. See pobysoTaylor_so_so. """
1848
    return(pobyso_taylor_so_so(function, degree, point))
1849

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

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

    
1962
def pobyso_univar_polynomial_print_reverse(polySa):
1963
    """ Legacy function. See pobyso_univar_polynomial_print_reverse_sa_sa. """
1964
    return(pobyso_univar_polynomial_print_reverse_sa_sa(polySa))
1965

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