Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 234

Historique | Voir | Annoter | Télécharger (79,08 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
def pobyso_dirty_find_zeros_so_so(funcSo, rangeSo):
365
    """
366
    Thin wrapper over sollya_lib_dirtyfindzeros()
367
    """
368
    return sollya_lib_dirtyfindzeros(funcSo, rangeSo)
369
# End pobys_dirty_find_zeros
370

    
371
def pobyso_float_list_so_sa(listSo):
372
    """
373
    Return a Sollya list of floating-point as a Sage list of 
374
    floating-point numbers.
375
    """
376
    listSa   = []
377
    ## Remember pobyso_get_list_elements_so_so returns more information
378
    #  than just the elements of the list (# elements, is_ellipti)
379
    listSaSo = pobyso_get_list_elements_so_so(listSo)[0]
380
    ## Search first for the maximum precision of the elements
381
    maxPrecSa = 0
382
    for floatSo in listSaSo:
383
        pobyso_autoprint(floatSo)
384
        curPrecSa =  pobyso_get_prec_of_constant_so_sa(floatSo)
385
        if curPrecSa > maxPrecSa:
386
            maxPrecSa = curPrecSa
387
    ##
388
    RF = RealField(maxPrecSa)
389
    ##
390
    for floatSo in listSaSo:
391
        listSa.append(pobyso_get_constant_as_rn_with_rf_so_sa(floatSo))
392
    return listSa
393
# End pobyso_float_list_so_sa
394

    
395
def pobyso_float_poly_sa_so(polySa, precSa = None):
396
    """
397
    Create a Sollya polynomial from a Sage RealField polynomial.
398
    """
399
    ## TODO: filter arguments.
400
    ## Precision. If a precision is given, convert the polynomial
401
    #  into the right polynomial field. If not convert it straight
402
    #  to Sollya.
403
    sollyaPrecChanged = False 
404
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
405
    if precSa is None:
406
        precSa = polySa.parent().base_ring().precision()
407
    if (precSa > initialSollyaPrecSa):
408
        assert precSa >= 2, "Precision change <2 requested"
409
        if precSa <= 2:
410
            print inspect.stack()[0][3], ": precision change <= 2 requested"
411
        precSo = pobyso_constant_from_int(precSa)
412
        pobyso_set_prec_so_so(precSo)
413
        sollya_lib_clear_obj(precSo)
414
        sollyaPrecChanged = True    
415
    ## Get exponents and coefficients.
416
    exponentsSa     = polySa.exponents()
417
    coefficientsSa  = polySa.coefficients()
418
    ## Build the polynomial.
419
    polySo = None
420
    for coefficientSa, exponentSa in zip(coefficientsSa, exponentsSa):
421
        #print coefficientSa.n(prec=precSa), exponentSa
422
        coefficientSo = \
423
            pobyso_constant_sa_so(coefficientSa)
424
        #pobyso_autoprint(coefficientSo)
425
        exponentSo = \
426
            pobyso_constant_from_int_sa_so(exponentSa)
427
        #pobyso_autoprint(exponentSo)
428
        monomialSo = sollya_lib_build_function_pow(
429
                       sollya_lib_build_function_free_variable(),
430
                       exponentSo)
431
        polyTermSo = sollya_lib_build_function_mul(coefficientSo,
432
                                                       monomialSo)
433
        if polySo is None:
434
            polySo = polyTermSo
435
        else:
436
            polySo = sollya_lib_build_function_add(polySo, polyTermSo)
437
    if sollyaPrecChanged:
438
        pobyso_set_prec_so_so(initialSollyaPrecSo)
439
        sollya_lib_clear_obj(initialSollyaPrecSo)
440
    return polySo
441
# End pobyso_float_poly_sa_so    
442

    
443
def pobyso_float_poly_so_sa(polySo, realFieldSa=None):
444
    """
445
    Convert a Sollya polynomial into a Sage floating-point polynomial.
446
    If no realField is given, a RealField corresponding to the maximum 
447
    precision of the coefficients is internally computed.
448
    The real field is not returned but can be easily retrieved from 
449
    the polynomial itself.
450
    ALGORITHM:
451
    - (optional) compute the RealField of the coefficients;
452
    - convert the Sollya expression into a Sage expression;
453
    - convert the Sage expression into a Sage polynomial
454
    """    
455
    if realFieldSa is None:
456
        expressionPrecSa = pobyso_get_max_prec_of_exp_so_sa(polySo)
457
        #print "Maximum precision of Sollya polynomial coefficients:", expressionPrecSa
458
        if expressionPrecSa < 2 or expressionPrecSa > 2147483391:
459
            print "Maximum degree of expression:", expressionPrecSa
460
        realFieldSa      = RealField(expressionPrecSa)
461
    #print "Sollya expression before...",
462
    #pobyso_autoprint(polySo)
463

    
464
    expressionSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo,
465
                                                             realFieldSa)
466
    #print "...Sollya expression after."
467
    #pobyso_autoprint(polySo)
468
    polyVariableSa = expressionSa.variables()[0]
469
    polyRingSa     = realFieldSa[str(polyVariableSa)]
470
    #print polyRingSa
471
    # Do not use the polynomial(expressionSa, ring=polyRingSa) form!
472
    polynomialSa = polyRingSa(expressionSa)
473
    polyCoeffsListSa = polynomialSa.coefficients()
474
    #for coeff in polyCoeffsListSa:
475
    #    print coeff.abs().n()
476
    return polynomialSa
477
# End pobyso_float_poly_so_sa
478

    
479
def pobyso_free_variable():
480
    """
481
    Ultra thin wrapper around the sollya_lib_function_build_free_variable function.
482
    """
483
    return sollya_lib_build_function_free_variable()
484
    
485
def pobyso_function_type_as_string(funcType):
486
    """ Legacy function. See pobyso_function_type_as_string_so_sa. """
487
    return(pobyso_function_type_as_string_so_sa(funcType))
488

    
489
def pobyso_function_type_as_string_so_sa(funcType):
490
    """
491
    Numeric Sollya function codes -> Sage mathematical function names.
492
    Notice that pow -> ^ (a la Sage, not a la Python).
493
    """
494
    if funcType == SOLLYA_BASE_FUNC_ABS:
495
        return "abs"
496
    elif funcType == SOLLYA_BASE_FUNC_ACOS:
497
        return "arccos"
498
    elif funcType == SOLLYA_BASE_FUNC_ACOSH:
499
        return "arccosh"
500
    elif funcType == SOLLYA_BASE_FUNC_ADD:
501
        return "+"
502
    elif funcType == SOLLYA_BASE_FUNC_ASIN:
503
        return "arcsin"
504
    elif funcType == SOLLYA_BASE_FUNC_ASINH:
505
        return "arcsinh"
506
    elif funcType == SOLLYA_BASE_FUNC_ATAN:
507
        return "arctan"
508
    elif funcType == SOLLYA_BASE_FUNC_ATANH:
509
        return "arctanh"
510
    elif funcType == SOLLYA_BASE_FUNC_CEIL:
511
        return "ceil"
512
    elif funcType == SOLLYA_BASE_FUNC_CONSTANT:
513
        return "cte"
514
    elif funcType == SOLLYA_BASE_FUNC_COS:
515
        return "cos"
516
    elif funcType == SOLLYA_BASE_FUNC_COSH:
517
        return "cosh"
518
    elif funcType == SOLLYA_BASE_FUNC_DIV:
519
        return "/"
520
    elif funcType == SOLLYA_BASE_FUNC_DOUBLE:
521
        return "double"
522
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEDOUBLE:
523
        return "doubleDouble"
524
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEEXTENDED:
525
        return "doubleDxtended"
526
    elif funcType == SOLLYA_BASE_FUNC_ERF:
527
        return "erf"
528
    elif funcType == SOLLYA_BASE_FUNC_ERFC:
529
        return "erfc"
530
    elif funcType == SOLLYA_BASE_FUNC_EXP:
531
        return "exp"
532
    elif funcType == SOLLYA_BASE_FUNC_EXP_M1:
533
        return "expm1"
534
    elif funcType == SOLLYA_BASE_FUNC_FLOOR:
535
        return "floor"
536
    elif funcType == SOLLYA_BASE_FUNC_FREE_VARIABLE:
537
        return "freeVariable"
538
    elif funcType == SOLLYA_BASE_FUNC_HALFPRECISION:
539
        return "halfPrecision"
540
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYCONSTANT:
541
        return "libraryConstant"
542
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYFUNCTION:
543
        return "libraryFunction"
544
    elif funcType == SOLLYA_BASE_FUNC_LOG:
545
        return "log"
546
    elif funcType == SOLLYA_BASE_FUNC_LOG_10:
547
        return "log10"
548
    elif funcType == SOLLYA_BASE_FUNC_LOG_1P:
549
        return "log1p"
550
    elif funcType == SOLLYA_BASE_FUNC_LOG_2:
551
        return "log2"
552
    elif funcType == SOLLYA_BASE_FUNC_MUL:
553
        return "*"
554
    elif funcType == SOLLYA_BASE_FUNC_NEARESTINT:
555
        return "round"
556
    elif funcType == SOLLYA_BASE_FUNC_NEG:
557
        return "__neg__"
558
    elif funcType == SOLLYA_BASE_FUNC_PI:
559
        return "pi"
560
    elif funcType == SOLLYA_BASE_FUNC_POW:
561
        return "^"
562
    elif funcType == SOLLYA_BASE_FUNC_PROCEDUREFUNCTION:
563
        return "procedureFunction"
564
    elif funcType == SOLLYA_BASE_FUNC_QUAD:
565
        return "quad"
566
    elif funcType == SOLLYA_BASE_FUNC_SIN:
567
        return "sin"
568
    elif funcType == SOLLYA_BASE_FUNC_SINGLE:
569
        return "single"
570
    elif funcType == SOLLYA_BASE_FUNC_SINH:
571
        return "sinh"
572
    elif funcType == SOLLYA_BASE_FUNC_SQRT:
573
        return "sqrt"
574
    elif funcType == SOLLYA_BASE_FUNC_SUB:
575
        return "-"
576
    elif funcType == SOLLYA_BASE_FUNC_TAN:
577
        return "tan"
578
    elif funcType == SOLLYA_BASE_FUNC_TANH:
579
        return "tanh"
580
    elif funcType == SOLLYA_BASE_FUNC_TRIPLEDOUBLE:
581
        return "tripleDouble"
582
    else:
583
        return None
584

    
585
def pobyso_get_constant(rnArgSa, constSo):
586
    """ Legacy function. See pobyso_get_constant_so_sa. """
587
    return pobyso_get_constant_so_sa(rnArgSa, constSo)
588
# End pobyso_get_constant
589

    
590
def pobyso_get_constant_so_sa(rnArgSa, constSo):
591
    """
592
    Set the value of rnArgSo to the value of constSo in MPFR_RNDN mode.
593
    rnArg must already exist and belong to some RealField.
594
    We assume that constSo points to a Sollya constant.
595
    """
596
    outcome = sollya_lib_get_constant(get_rn_value(rnArgSa), constSo)
597
    if outcome == 0: # Failure because constSo is not a constant expression.
598
        return None
599
    else:
600
        return outcome
601
# End  pobyso_get_constant_so_sa
602
   
603
def pobyso_get_constant_as_rn(ctExpSo):
604
    """ 
605
    Legacy function. See pobyso_get_constant_as_rn_so_sa. 
606
    """ 
607
    return(pobyso_get_constant_as_rn_so_sa(ctExpSo))
608
    
609
def pobyso_get_constant_as_rn_so_sa(constExpSo):
610
    """
611
    Get a Sollya constant as a Sage "real number".
612
    The precision of the floating-point number returned is that of the Sollya
613
    constant.
614
    """
615
    #print "Before computing precision of variable..."
616
    #pobyso_autoprint(constExpSo)
617
    precisionSa  = pobyso_get_prec_of_constant_so_sa(constExpSo)
618
    #print "precisionSa:", precisionSa
619
    ## If the expression can not be exactly converted, None is returned.
620
    #  In this case opt for the Sollya current expression.
621
    if precisionSa is None:
622
        precisionSa = pobyso_get_prec_so_sa()
623
    RRRR = RealField(precisionSa)
624
    rnSa = RRRR(0)
625
    outcome = sollya_lib_get_constant(get_rn_value(rnSa), constExpSo)
626
    if outcome == 0:
627
        return None
628
    else:
629
        return rnSa
630
# End pobyso_get_constant_as_rn_so_sa
631

    
632
def pobyso_get_constant_as_rn_with_rf(ctExp, realField):
633
    """ 
634
    Legacy function. See pobyso_get_constant_as_rn_with_rf_so_sa.
635
    """
636
    return pobyso_get_constant_as_rn_with_rf_so_sa(ctExp, realField)
637
# End pobyso_get_constant_as_rn_with_rf
638
    
639
def pobyso_get_constant_as_rn_with_rf_so_sa(ctExpSo, realFieldSa = None):
640
    """
641
    Get a Sollya constant as a Sage "real number".
642
    If no real field is specified, the precision of the floating-point number 
643
    returned is that of the Sollya constant.
644
    Otherwise is is that of the real field. Hence rounding may happen.
645
    """
646
    if realFieldSa is None:
647
        return pobyso_get_constant_as_rn_so_sa(ctExpSo)
648
    rnSa = realFieldSa(0)
649
    outcome = sollya_lib_get_constant(get_rn_value(rnSa), ctExpSo)
650
    if outcome == 0:
651
        return None
652
    else:
653
        return rnSa
654
# End pobyso_get_constant_as_rn_with_rf_so_sa
655

    
656
def pobyso_get_free_variable_name():
657
    """ 
658
    Legacy function. See pobyso_get_free_variable_name_so_sa.
659
    """
660
    return(pobyso_get_free_variable_name_so_sa())
661

    
662
def pobyso_get_free_variable_name_so_sa():
663
    return sollya_lib_get_free_variable_name()
664
    
665
def pobyso_get_function_arity(expressionSo):
666
    """ 
667
    Legacy function. See pobyso_get_function_arity_so_sa.
668
    """
669
    return(pobyso_get_function_arity_so_sa(expressionSo))
670

    
671
def pobyso_get_function_arity_so_sa(expressionSo):
672
    arity = c_int(0)
673
    sollya_lib_get_function_arity(byref(arity),expressionSo)
674
    return int(arity.value)
675

    
676
def pobyso_get_head_function(expressionSo):
677
    """ 
678
    Legacy function. See pobyso_get_head_function_so_sa. 
679
    """
680
    return(pobyso_get_head_function_so_sa(expressionSo)) 
681

    
682
def pobyso_get_head_function_so_sa(expressionSo):
683
    functionType = c_int(0)
684
    sollya_lib_get_head_function(byref(functionType), expressionSo)
685
    return int(functionType.value)
686

    
687
def pobyso_get_interval_from_range_so_sa(soRange, realIntervalFieldSa = None ):
688
    """
689
    Return the Sage interval corresponding to the Sollya range argument.
690
    If no reaIntervalField is passed as an argument, the interval bounds are not
691
    rounded: they are elements of RealIntervalField of the "right" precision
692
    to hold all the digits.
693
    """
694
    prec = c_int(0)
695
    if realIntervalFieldSa is None:
696
        retval = sollya_lib_get_prec_of_range(byref(prec), soRange, None)
697
        if retval == 0:
698
            return None
699
        realIntervalFieldSa = RealIntervalField(prec.value)
700
    intervalSa = realIntervalFieldSa(0,0)
701
    retval = \
702
        sollya_lib_get_interval_from_range(get_interval_value(intervalSa),\
703
                                           soRange)
704
    if retval == 0:
705
        return None
706
    return intervalSa
707
# End pobyso_get_interval_from_range_so_sa
708

    
709
def pobyso_get_list_elements(soObj):
710
    """ Legacy function. See pobyso_get_list_elements_so_so. """
711
    return pobyso_get_list_elements_so_so(soObj)
712
 
713
def pobyso_get_list_elements_so_so(objectListSo):
714
    """
715
    Get the Sollya list elements as a Sage/Python array of Sollya objects.
716
    
717
    INPUT:
718
    - objectListSo: a Sollya list of Sollya objects.
719
    
720
    OUTPUT:
721
    - a Sage/Python tuple made of:
722
      - a Sage/Python list of Sollya objects,
723
      - a Sage/Python int holding the number of elements,
724
      - a Sage/Python int stating (!= 0) that the list is end-elliptic.
725
    NOTE::
726
        We recover the addresses of the Sollya object from the list of pointers
727
        returned by sollya_lib_get_list_elements. The list itself is freed.
728
    TODO::
729
        Figure out what to do with numElements since the number of elements
730
        can easily be recovered from the list itself. 
731
        Ditto for isEndElliptic.
732
    """
733
    listAddress = POINTER(c_longlong)()
734
    numElements = c_int(0)
735
    isEndElliptic = c_int(0)
736
    listAsSageList = []
737
    result = sollya_lib_get_list_elements(byref(listAddress),\
738
                                          byref(numElements),\
739
                                          byref(isEndElliptic),\
740
                                          objectListSo)
741
    if result == 0 :
742
        return None
743
    for i in xrange(0, numElements.value, 1):
744
       #listAsSageList.append(sollya_lib_copy_obj(listAddress[i]))
745
       listAsSageList.append(listAddress[i])
746
       # Clear each of the elements returned by Sollya.
747
       #sollya_lib_clear_obj(listAddress[i])
748
    # Free the list itself.   
749
    sollya_lib_free(listAddress)
750
    return (listAsSageList, numElements.value, isEndElliptic.value)
751

    
752
def pobyso_get_max_prec_of_exp(soExp):
753
    """ Legacy function. See pobyso_get_max_prec_of_exp_so_sa. """
754
    return pobyso_get_max_prec_of_exp_so_sa(soExp)
755

    
756
def pobyso_get_max_prec_of_exp_so_sa(expSo):
757
    """
758
    Get the maximum precision used for the numbers in a Sollya expression.
759
    
760
    Arguments:
761
    soExp -- a Sollya expression pointer
762
    Return value:
763
    A Python integer
764
    TODO: 
765
    - error management;
766
    - correctly deal with numerical type such as DOUBLEEXTENDED.
767
    """
768
    if expSo is None:
769
        print inspect.stack()[0][3], ": expSo is None."
770
        return 0
771
    maxPrecision = 0
772
    minConstPrec = 0
773
    currentConstPrec = 0
774
    #pobyso_autoprint(expSo)
775
    operator = pobyso_get_head_function_so_sa(expSo)
776
    if (operator != SOLLYA_BASE_FUNC_CONSTANT) and \
777
    (operator != SOLLYA_BASE_FUNC_FREE_VARIABLE):
778
        (arity, subexpressions) = pobyso_get_subfunctions_so_sa(expSo)
779
        for i in xrange(arity):
780
            maxPrecisionCandidate = \
781
                pobyso_get_max_prec_of_exp_so_sa(subexpressions[i])
782
            if maxPrecisionCandidate > maxPrecision:
783
                maxPrecision = maxPrecisionCandidate
784
        return maxPrecision
785
    elif operator == SOLLYA_BASE_FUNC_CONSTANT:
786
        #minConstPrec = pobyso_get_min_prec_of_constant_so_sa(expSo)
787
        #currentConstPrec = pobyso_get_min_prec_of_constant_so_sa(soExp)
788
        #print minConstPrec, " - ", currentConstPrec 
789
        return pobyso_get_min_prec_of_constant_so_sa(expSo)
790
    
791
    elif operator == SOLLYA_BASE_FUNC_FREE_VARIABLE:
792
        return 0
793
    else:
794
        print "pobyso_get_max_prec_of_exp_so_sa: unexepected operator."
795
        return 0
796

    
797
def pobyso_get_min_prec_of_constant_so_sa(constExpSo):
798
    """
799
    Get the minimum precision necessary to represent the value of a Sollya
800
    constant.
801
    MPFR_MIN_PREC and powers of 2 are taken into account.
802
    We assume that constExpSo is a pointer to a Sollay constant expression.
803
    """
804
    constExpAsRnSa = pobyso_get_constant_as_rn_so_sa(constExpSo)
805
    return(min_mpfr_size(get_rn_value(constExpAsRnSa)))
806

    
807
def pobyso_get_poly_so_sa(polySo, realFieldSa=None):
808
    """
809
    Convert a Sollya polynomial into a Sage polynomial.
810
    Legacy function. Use pobyso_float_poly_so_sa() instead.
811
    """    
812
    return pobyso_float_poly_so_sa(polySo,realFieldSa)
813
# End pobyso_get_poly_so_sa
814

    
815
def pobyso_get_prec():
816
    """ Legacy function. See pobyso_get_prec_so_sa(). """
817
    return pobyso_get_prec_so_sa()
818

    
819
def pobyso_get_prec_so():
820
    """
821
    Get the current default precision in Sollya.
822
    The return value is a Sollya object.
823
    Usefull when modifying the precision back and forth by avoiding
824
    extra conversions.
825
    """
826
    return sollya_lib_get_prec(None)
827
    
828
def pobyso_get_prec_so_sa():
829
    """
830
    Get the current default precision in Sollya.
831
    The return value is Sage/Python int.
832
    """
833
    precSo = sollya_lib_get_prec()
834
    precSa = pobyso_constant_from_int_so_sa(precSo)
835
    sollya_lib_clear_obj(precSo)
836
    return precSa
837
# End pobyso_get_prec_so_sa.
838

    
839
def pobyso_get_prec_so_so_sa():
840
    """
841
    Return the current precision both as a Sollya object and a
842
    Sage integer as hybrid tuple.
843
    To avoid multiple calls for precision manipulations. 
844
    """
845
    precSo = sollya_lib_get_prec()
846
    precSa = pobyso_constant_from_int_so_sa(precSo)
847
    return (precSo, int(precSa))
848
    
849
def pobyso_get_prec_of_constant(ctExpSo):
850
    """ Legacy function. See pobyso_get_prec_of_constant_so_sa. """
851
    return pobyso_get_prec_of_constant_so_sa(ctExpSo)
852

    
853
def pobyso_get_prec_of_constant_so_sa(ctExpSo):
854
    """
855
    Tries to find a precision to represent ctExpSo without rounding.
856
    If not possible, returns None.
857
    """
858
    #print "Entering pobyso_get_prec_of_constant_so_sa..."
859
    prec = c_int(0)
860
    retc = sollya_lib_get_prec_of_constant(byref(prec), ctExpSo, None)
861
    if retc == 0:
862
        #print "pobyso_get_prec_of_constant_so_sa failed."
863
        return None
864
    #print "...exiting pobyso_get_prec_of_constant_so_sa."
865
    return int(prec.value)
866

    
867
def pobyso_get_prec_of_range_so_sa(rangeSo):
868
    """
869
    Returns the number of bits elements of a range are coded with.
870
    """
871
    prec = c_int(0)
872
    retc = sollya_lib_get_prec_of_range(byref(prec), rangeSo, None)
873
    if retc == 0:
874
        return(None)
875
    return int(prec.value)
876
# End pobyso_get_prec_of_range_so_sa()
877

    
878
def pobyso_get_sage_exp_from_sollya_exp(sollyaExpSo, realField = RR):
879
    """ Legacy function. See pobyso_get_sage_exp_from_sollya_exp_so_sa. """
880
    return pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, 
881
                                                     realField = RR)
882

    
883
def pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, realFieldSa = RR):
884
    """
885
    Get a Sage expression from a Sollya expression. 
886
    Currently only tested with polynomials with floating-point coefficients.
887
    Notice that, in the returned polynomial, the exponents are RealNumbers.
888
    """
889
    #pobyso_autoprint(sollyaExp)
890
    operatorSa = pobyso_get_head_function_so_sa(sollyaExpSo)
891
    sollyaLibFreeVariableName = sollya_lib_get_free_variable_name()
892
    ## Get rid of the "_"'s in "_x_", if any.
893
    sollyaLibFreeVariableName = re.sub('_', '', sollyaLibFreeVariableName)
894
    # Constants and the free variable are special cases.
895
    # All other operator are dealt with in the same way.
896
    if (operatorSa != SOLLYA_BASE_FUNC_CONSTANT) and \
897
       (operatorSa != SOLLYA_BASE_FUNC_FREE_VARIABLE):
898
        (aritySa, subexpressionsSa) = pobyso_get_subfunctions_so_sa(sollyaExpSo)
899
        if aritySa == 1:
900
            sageExpSa = eval(pobyso_function_type_as_string_so_sa(operatorSa) + \
901
            "(" + pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], \
902
            realFieldSa) + ")")
903
        elif aritySa == 2:
904
            # We do not get through the preprocessor.
905
            # The "^" operator is then a special case.
906
            if operatorSa == SOLLYA_BASE_FUNC_POW:
907
                operatorAsStringSa = "**"
908
            else:
909
                operatorAsStringSa = \
910
                    pobyso_function_type_as_string_so_sa(operatorSa)
911
            sageExpSa = \
912
              eval("pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], realFieldSa)"\
913
              + " " + operatorAsStringSa + " " + \
914
                   "pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[1], realFieldSa)")
915
        # We do not know yet how to deal with arity >= 3 
916
        # (is there any in Sollya anyway?).
917
        else:
918
            sageExpSa = eval('None')
919
        return sageExpSa
920
    elif operatorSa == SOLLYA_BASE_FUNC_CONSTANT:
921
        #print "This is a constant"
922
        return pobyso_get_constant_as_rn_with_rf_so_sa(sollyaExpSo, realFieldSa)
923
    elif operatorSa == SOLLYA_BASE_FUNC_FREE_VARIABLE:
924
        #print "This is the free variable"
925
        return eval(sollyaLibFreeVariableName)
926
    else:
927
        print "Unexpected"
928
        return eval('None')
929
# End pobyso_get_sage_exp_from_sollya_exp_so_sa
930

    
931

    
932
def pobyso_get_subfunctions(expressionSo):
933
    """ Legacy function. See pobyso_get_subfunctions_so_sa. """
934
    return pobyso_get_subfunctions_so_sa(expressionSo) 
935
# End pobyso_get_subfunctions.
936
 
937
def pobyso_get_subfunctions_so_sa(expressionSo):
938
    """
939
    Get the subfunctions of an expression.
940
    Return the number of subfunctions and the list of subfunctions addresses.
941
    S.T.: Could not figure out another way than that ugly list of declarations
942
    to recover the addresses of the subfunctions.
943
    We limit ourselves to arity 8 functions. 
944
    """
945
    subf0 = c_int(0)
946
    subf1 = c_int(0)
947
    subf2 = c_int(0)
948
    subf3 = c_int(0)
949
    subf4 = c_int(0)
950
    subf5 = c_int(0)
951
    subf6 = c_int(0)
952
    subf7 = c_int(0)
953
    subf8 = c_int(0)
954
    arity = c_int(0)
955
    nullPtr = POINTER(c_int)()
956
    sollya_lib_get_subfunctions(expressionSo, byref(arity), \
957
      byref(subf0), byref(subf1), byref(subf2), byref(subf3), \
958
      byref(subf4), byref(subf5),\
959
      byref(subf6), byref(subf7), byref(subf8), nullPtr, None) 
960
#    byref(cast(subfunctions[0], POINTER(c_int))), \
961
#    byref(cast(subfunctions[0], POINTER(c_int))), \
962
#    byref(cast(subfunctions[2], POINTER(c_int))), \
963
#    byref(cast(subfunctions[3], POINTER(c_int))), \
964
#    byref(cast(subfunctions[4], POINTER(c_int))), \
965
#    byref(cast(subfunctions[5], POINTER(c_int))), \
966
#    byref(cast(subfunctions[6], POINTER(c_int))), \
967
#    byref(cast(subfunctions[7], POINTER(c_int))), \
968
#    byref(cast(subfunctions[8], POINTER(c_int))), nullPtr)
969
    subfunctions = [subf0, subf1, subf2, subf3, subf4, subf5, subf6, subf7, \
970
                    subf8]
971
    subs = []
972
    if arity.value > pobyso_max_arity:
973
        return(0,[])
974
    for i in xrange(arity.value):
975
        subs.append(int(subfunctions[i].value))
976
        #print subs[i]
977
    return (int(arity.value), subs)
978
# End pobyso_get_subfunctions_so_sa
979
    
980
def pobyso_guess_degree_sa_sa(functionSa, intervalSa, approxErrorSa, 
981
                              weightSa=None, degreeBoundSa=None):
982
    """
983
    Sa_sa variant of the solly_guessdegree function.
984
    Return 0 if something goes wrong.
985
    """
986
    functionAsStringSa = functionSa._assume_str().replace('_SAGE_VAR_', '')
987
    functionSo = pobyso_parse_string_sa_so(functionAsStringSa)
988
    if pobyso_is_error_so_sa(functionSo):
989
        sollya_lib_clear_obj(functionSo)
990
        return 0
991
    rangeSo = pobyso_interval_to_range_sa_so(intervalSa)
992
    # The approximation error is expected to be a floating point number.
993
    if pobyso_is_floating_point_number_sa_sa(approxErrorSa):
994
        approxErrorSo = pobyso_constant_sa_so(approxErrorSa)
995
    else:
996
        approxErrorSo = pobyso_constant_sa_so(RR(approxErrorSa))
997
    if not weightSa is None:
998
        weightAsStringSa = weightSa._assume_str().replace('_SAGE_VAR_', '')
999
        weightSo = pobyso_parse_string_sa_so(weightAsStringSa)
1000
        if pobyso_is_error_so_sa(weightSo):
1001
            sollya_lib_clear_obj(functionSo)
1002
            sollya_lib_clear_obj(rangeSo)
1003
            sollya_lib_clear_obj(approxErrorSo)   
1004
            sollya_lib_clear_obj(weightSo)
1005
            return 0   
1006
    else:
1007
        weightSo = None
1008
    if not degreeBoundSa is None:
1009
        degreeBoundSo = pobyso_constant_from_int_sa_so(degreeBoundSa)
1010
    else:
1011
        degreeBoundSo = None
1012
    guessedDegreeSa = pobyso_guess_degree_so_sa(functionSo,
1013
                                                rangeSo,
1014
                                                approxErrorSo,
1015
                                                weightSo,
1016
                                                degreeBoundSo)
1017
    sollya_lib_clear_obj(functionSo)
1018
    sollya_lib_clear_obj(rangeSo)
1019
    sollya_lib_clear_obj(approxErrorSo)
1020
    if not weightSo is None:
1021
        sollya_lib_clear_obj(weightSo)
1022
    if not degreeBoundSo is None:
1023
        sollya_lib_clear_obj(degreeBoundSo)
1024
    return guessedDegreeSa
1025
# End poyso_guess_degree_sa_sa
1026

    
1027
def pobyso_guess_degree_so_sa(functionSo, rangeSo, errorSo, weightSo=None, \
1028
                              degreeBoundSo=None):
1029
    """
1030
    Thin wrapper around the guessdegree function.
1031
    Nevertheless, some precision control stuff has been appended.
1032
    """
1033
    # Deal with Sollya internal precision issues: if it is too small, 
1034
    # compared with the error, increases it to about twice -log2(error).
1035
    errorSa = pobyso_get_constant_as_rn_with_rf_so_sa(errorSo)
1036
    log2ErrorSa = errorSa.log2()
1037
    if log2ErrorSa < 0:
1038
        neededPrecisionSa = int(2 * int(-log2ErrorSa) / 64) * 64
1039
    else:
1040
        neededPrecisionSa = int(2 * int(log2ErrorSa) / 64) * 64
1041
    #print "Needed precision:", neededPrecisionSa
1042
    sollyaPrecisionHasChanged = False
1043
    (initialPrecSo, initialPrecSa) = pobyso_get_prec_so_so_sa()
1044
    if neededPrecisionSa > initialPrecSa:
1045
        if neededPrecisionSa <= 2:
1046
            print inspect.stack()[0][3], ": precision change <= 2 requested."
1047
        pobyso_set_prec_sa_so(neededPrecisionSa)
1048
        sollyaPrecisionHasChanged = True
1049
    #print "Guessing degree..."
1050
    # weightSo and degreeBoundsSo are optional arguments.
1051
    # As declared, sollya_lib_guessdegree must take 5 arguments.
1052
    if weightSo is None:
1053
        degreeRangeSo = sollya_lib_guessdegree(functionSo, rangeSo, errorSo, 
1054
                                               0, 0, None)
1055
    elif degreeBoundSo is None:
1056
        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, \
1057
                                                errorSo, weightSo, 0, None)
1058
    else:
1059
        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, errorSo, \
1060
                                                weightSo, degreeBoundSo, None)
1061
    #print "...degree guess done."
1062
    # Restore internal precision, if applicable.
1063
    if sollyaPrecisionHasChanged:
1064
        pobyso_set_prec_so_so(initialPrecSo)
1065
        sollya_lib_clear_obj(initialPrecSo)
1066
    degreeIntervalSa = pobyso_range_to_interval_so_sa(degreeRangeSo)
1067
    sollya_lib_clear_obj(degreeRangeSo)
1068
    # When ok, both bounds match.
1069
    # When the degree bound is too low, the upper bound is the degree
1070
    # for which the error can be honored.
1071
    # When it really goes wrong, the upper bound is infinity. 
1072
    if degreeIntervalSa.lower() == degreeIntervalSa.upper():
1073
        return int(degreeIntervalSa.lower())
1074
    else:
1075
        if degreeIntervalSa.upper().is_infinity():
1076
            return None
1077
        else:
1078
            return int(degreeIntervalSa.upper())
1079
    # End pobyso_guess_degree_so_sa
1080

    
1081
def pobyso_inf_so_so(intervalSo):
1082
    """
1083
    Very thin wrapper around sollya_lib_inf().
1084
    """
1085
    return sollya_lib_inf(intervalSo)
1086
# End pobyso_inf_so_so.
1087
    
1088
def pobyso_infnorm_so_so(func, interval, file = None, intervalList = None):
1089
    print "Do not use this function. User pobyso_supnorm_so_so instead."
1090
    return None
1091

    
1092
def pobyso_interval_to_range_sa_so(intervalSa, precisionSa=None):
1093
    if precisionSa is None:
1094
        precisionSa = intervalSa.parent().precision()
1095
    intervalSo = pobyso_bounds_to_range_sa_so(intervalSa.lower(),\
1096
                                              intervalSa.upper(),\
1097
                                              precisionSa)
1098
    return intervalSo
1099
# End pobyso_interval_to_range_sa_so
1100

    
1101
def pobyso_is_error_so_sa(objSo):
1102
    """
1103
    Thin wrapper around the sollya_lib_obj_is_error() function.
1104
    """
1105
    if sollya_lib_obj_is_error(objSo) != 0:
1106
        return True
1107
    else:
1108
        return False
1109
# End pobyso_is_error-so_sa
1110

    
1111
def pobyso_is_floating_point_number_sa_sa(numberSa):
1112
    """
1113
    Check whether a Sage number is floating point.
1114
    Exception stuff added because numbers other than
1115
    floating-point ones do not have the is_real() attribute.
1116
    """
1117
    try:
1118
        return numberSa.is_real()
1119
    except AttributeError:
1120
        return False
1121
# End pobyso_is_floating_piont_number_sa_sa
1122

    
1123
def pobyso_is_range_so_sa(rangeCandidateSo):
1124
    """
1125
    Thin wrapper over sollya_lib_is_range.
1126
    """
1127
    return sollya_lib_obj_is_range(rangeCandidateSo) != 0
1128
    
1129
# End pobyso_is_range_so_sa
1130

    
1131

    
1132
def pobyso_lib_init():
1133
    sollya_lib_init(None)
1134

    
1135
def pobyso_lib_close():
1136
    sollya_lib_close(None)
1137
    
1138
def pobyso_name_free_variable(freeVariableNameSa):
1139
    """ Legacy function. See pobyso_name_free_variable_sa_so. """
1140
    pobyso_name_free_variable_sa_so(freeVariableNameSa)
1141

    
1142
def pobyso_name_free_variable_sa_so(freeVariableNameSa):
1143
    """
1144
    Set the free variable name in Sollya from a Sage string.
1145
    """
1146
    sollya_lib_name_free_variable(freeVariableNameSa)
1147

    
1148
def pobyso_parse_string(string):
1149
    """ Legacy function. See pobyso_parse_string_sa_so. """
1150
    return pobyso_parse_string_sa_so(string)
1151
 
1152
def pobyso_parse_string_sa_so(string):
1153
    """
1154
    Get the Sollya expression computed from a Sage string or
1155
    a Sollya error object if parsing failed.
1156
    """
1157
    return sollya_lib_parse_string(string)
1158

    
1159
def pobyso_precision_so_sa(ctExpSo):
1160
    """
1161
    Computes the necessary precision to represent a number.
1162
    If x is not zero, it can be uniquely written as x = m · 2e 
1163
    where m is an odd integer and e is an integer. 
1164
    precision(x) returns the number of bits necessary to write m 
1165
    in binary (i.e. ceil(log2(m))).
1166
    """
1167
    #TODO: take care of the special case: 0, @NaN@, @Inf@
1168
    precisionSo = sollya_lib_precision(ctExpSo)
1169
    precisionSa = pobyso_constant_from_int_so_sa(precisionSo)
1170
    sollya_lib_clear_obj(precisionSo)
1171
    return precisionSa
1172
# End pobyso_precision_so_sa
1173

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

    
1340
def pobyso_polynomial_degree_so_sa(polySo):
1341
    """
1342
    Return the degree of a Sollya polynomial as a Sage int.
1343
    """
1344
    degreeSo = sollya_lib_degree(polySo)
1345
    return pobyso_constant_from_int_so_sa(degreeSo)
1346
# End pobyso_polynomial_degree_so_sa
1347

    
1348
def pobyso_polynomial_degree_so_so(polySo):
1349
    """
1350
    Thin wrapper around lib_sollya_degree().
1351
    """
1352
    return sollya_lib_degree(polySo)
1353
# End pobyso_polynomial_degree_so_so
1354
                                                                  
1355
def pobyso_range(rnLowerBound, rnUpperBound):
1356
    """ Legacy function. See pobyso_range_sa_so. """
1357
    return pobyso_range_sa_so(rnLowerBound, rnUpperBound) 
1358

    
1359
def pobyso_range_from_bounds_sa_so(rnLowerBound, rnUpperBound, precSa = None):
1360
    """
1361
    Create a Sollya range from 2 Sage real numbers as bounds
1362
    """
1363
    # TODO check precision stuff.
1364
    sollyaPrecChanged = False
1365
    (initialSollyaPrecSo, initialSollyaPrecSa) = \
1366
        pobyso_get_prec_so_so_sa()
1367
    if precSa is None:
1368
        precSa = max(rnLowerBound.parent().prec(), rnUpperBound.parent().prec())
1369
    if precSa > initialSollyaPrecSa:
1370
        if precSa <= 2:
1371
            print inspect.stack()[0][3], ": precision change <= 2 requested."
1372
        pobyso_set_prec_sa_so(precSa)
1373
        sollyaPrecChanged = True
1374
    rangeSo = sollya_lib_range_from_bounds(get_rn_value(rnLowerBound),
1375
                                           get_rn_value(rnUpperBound))
1376
    if sollyaPrecChanged:
1377
        pobyso_set_prec_so_so(initialSollyaPrecSo)
1378
        sollya_lib_clear_obj(initialSollyaPrecSo)
1379
    return rangeSo
1380
# End pobyso_range_from_bounds_sa_so
1381

    
1382
def pobyso_range_max_abs_so_so(rangeSo):
1383
    """
1384
    Return, as Sollya constant, the maximum absolute value of a Sollay range.
1385
    """
1386
    lowerBoundSo = sollya_lib_inf(rangeSo)
1387
    upperBoundSo = sollya_lib_sup(rangeSo)
1388
    #
1389
    lowerBoundSo = sollya_lib_build_function_abs(lowerBoundSo)
1390
    upperBoundSo = sollya_lib_build_function_abs(upperBoundSo)
1391
    #pobyso_autoprint(lowerBoundSo)
1392
    #pobyso_autoprint(upperBoundSo)
1393
    #
1394
    maxAbsSo = sollya_lib_max(lowerBoundSo, upperBoundSo, None)
1395
    #sollya_lib_clear_obj(lowerBoundSo)
1396
    #sollya_lib_clear_obj(upperBoundSo)
1397
    return maxAbsSo
1398
# End pobyso_range_max_abs_so_so
1399
 
1400
def pobyso_range_to_interval_so_sa(rangeSo, realIntervalFieldSa = None):
1401
    """
1402
    Get a Sage interval from a Sollya range.
1403
    If no realIntervalField is given as a parameter, the Sage interval
1404
    precision is that of the Sollya range.
1405
    Otherwise, the precision is that of the realIntervalField. In this case
1406
    rounding may happen.
1407
    """
1408
    if realIntervalFieldSa is None:
1409
        precSa = pobyso_get_prec_of_range_so_sa(rangeSo)
1410
        realIntervalFieldSa = RealIntervalField(precSa)
1411
    intervalSa = \
1412
        pobyso_get_interval_from_range_so_sa(rangeSo, realIntervalFieldSa) 
1413
    return intervalSa
1414
# End pobyso_range_to_interval_so_sa
1415

    
1416
def pobyso_rat_poly_sa_so(polySa, precSa = None):
1417
    """
1418
    Create a Sollya polynomial from a Sage rational polynomial.
1419
    """
1420
    ## TODO: filter arguments.
1421
    ## Precision. If no precision is given, use the current precision
1422
    #  of Sollya.
1423
    if precSa is None:
1424
        precSa =  pobyso_get_prec_so_sa()
1425
    #print "Precision:",  precSa
1426
    RRR = RealField(precSa)
1427
    ## Create a Sage polynomial in the "right" precision.
1428
    P_RRR = RRR[polySa.variables()[0]]
1429
    polyFloatSa = P_RRR(polySa)
1430
    ## Make sure no precision is provided: pobyso_float_poly_sa_so will
1431
    #  recover it all by itself and not make an extra conversion.
1432
    return pobyso_float_poly_sa_so(polyFloatSa)
1433
    
1434
# End pobyso_rat_poly_sa_so    
1435
    
1436
def pobyso_remez_canonical_sa_sa(func, \
1437
                                 degree, \
1438
                                 lowerBound, \
1439
                                 upperBound, \
1440
                                 weight = None, \
1441
                                 quality = None):
1442
    """
1443
    All arguments are Sage/Python.
1444
    The functions (func and weight) must be passed as expressions or strings.
1445
    Otherwise the function fails. 
1446
    The return value is a Sage polynomial.
1447
    """
1448
    var('zorglub')    # Dummy variable name for type check only. Type of 
1449
    # zorglub is "symbolic expression".
1450
    polySo = pobyso_remez_canonical_sa_so(func, \
1451
                                 degree, \
1452
                                 lowerBound, \
1453
                                 upperBound, \
1454
                                 weight, \
1455
                                 quality)
1456
    # String test
1457
    if parent(func) == parent("string"):
1458
        functionSa = eval(func)
1459
    # Expression test.
1460
    elif type(func) == type(zorglub):
1461
        functionSa = func
1462
    else:
1463
        return None
1464
    #
1465
    maxPrecision = 0
1466
    if polySo is None:
1467
        return(None)
1468
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
1469
    RRRRSa = RealField(maxPrecision)
1470
    polynomialRingSa = RRRRSa[functionSa.variables()[0]]
1471
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, RRRRSa)
1472
    polySa = polynomial(expSa, polynomialRingSa)
1473
    sollya_lib_clear_obj(polySo)
1474
    return(polySa)
1475
# End pobyso_remez_canonical_sa_sa
1476
    
1477
def pobyso_remez_canonical(func, \
1478
                           degree, \
1479
                           lowerBound, \
1480
                           upperBound, \
1481
                           weight = "1", \
1482
                           quality = None):
1483
    """ Legacy function. See pobyso_remez_canonical_sa_so. """
1484
    return(pobyso_remez_canonical_sa_so(func, \
1485
                                        degree, \
1486
                                        lowerBound, \
1487
                                        upperBound, \
1488
                                        weight, \
1489
                                        quality))
1490
# End pobyso_remez_canonical.
1491

    
1492
def pobyso_remez_canonical_sa_so(func, \
1493
                                 degree, \
1494
                                 lowerBound, \
1495
                                 upperBound, \
1496
                                 weight = None, \
1497
                                 quality = None):
1498
    """
1499
    All arguments are Sage/Python.
1500
    The functions (func and weight) must be passed as expressions or strings.
1501
    Otherwise the function fails. 
1502
    The return value is a pointer to a Sollya function.
1503
    lowerBound and upperBound mus be reals.
1504
    """
1505
    var('zorglub')    # Dummy variable name for type check only. Type of
1506
    # zorglub is "symbolic expression".
1507
    currentVariableNameSa = None
1508
    # The func argument can be of different types (string, 
1509
    # symbolic expression...)
1510
    if parent(func) == parent("string"):
1511
        localFuncSa = sage_eval(func,globals())
1512
        if len(localFuncSa.variables()) > 0:
1513
            currentVariableNameSa = localFuncSa.variables()[0]
1514
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1515
            functionSo = \
1516
              sollya_lib_parse_string(localFuncSa._assume_str().replace('_SAGE_VAR_', ''))
1517
    # Expression test.
1518
    elif type(func) == type(zorglub):
1519
        # Until we are able to translate Sage expressions into Sollya 
1520
        # expressions : parse the string version.
1521
        if len(func.variables()) > 0:
1522
            currentVariableNameSa = func.variables()[0]
1523
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1524
            functionSo = \
1525
              sollya_lib_parse_string(func._assume_str().replace('_SAGE_VAR_', ''))
1526
    else:
1527
        return(None)
1528
    if weight is None: # No weight given -> 1.
1529
        weightSo = pobyso_constant_1_sa_so()
1530
    elif parent(weight) == parent("string"): # Weight given as string: parse it.
1531
        weightSo = sollya_lib_parse_string(func)
1532
    elif type(weight) == type(zorglub): # Weight given as symbolice expression.
1533
        functionSo = \
1534
          sollya_lib_parse_string_sa_so(weight._assume_str().replace('_SAGE_VAR_', ''))
1535
    else:
1536
        return(None)
1537
    degreeSo = pobyso_constant_from_int(degree)
1538
    rangeSo = pobyso_bounds_to_range_sa_so(lowerBound, upperBound)
1539
    if not quality is None:
1540
        qualitySo= pobyso_constant_sa_so(quality)
1541
    else:
1542
        qualitySo = None
1543
        
1544
    remezPolySo = sollya_lib_remez(functionSo, \
1545
                                   degreeSo, \
1546
                                   rangeSo, \
1547
                                   weightSo, \
1548
                                   qualitySo, \
1549
                                   None)
1550
    sollya_lib_clear_obj(functionSo)
1551
    sollya_lib_clear_obj(degreeSo)
1552
    sollya_lib_clear_obj(rangeSo)
1553
    sollya_lib_clear_obj(weightSo)
1554
    if not qualitySo is None:
1555
        sollya_lib_clear_obj(qualitySo)
1556
    return(remezPolySo)
1557
# End pobyso_remez_canonical_sa_so
1558

    
1559
def pobyso_remez_canonical_so_so(funcSo, \
1560
                                 degreeSo, \
1561
                                 rangeSo, \
1562
                                 weightSo = pobyso_constant_1_sa_so(),\
1563
                                 qualitySo = None):
1564
    """
1565
    All arguments are pointers to Sollya objects.
1566
    The return value is a pointer to a Sollya function.
1567
    """
1568
    if not sollya_lib_obj_is_function(funcSo):
1569
        return(None)
1570
    return(sollya_lib_remez(funcSo, degreeSo, rangeSo, weightSo, qualitySo, None))
1571
# End pobyso_remez_canonical_so_so.
1572

    
1573
def pobyso_round_coefficients_progressive_so_so(polySo, 
1574
                                                funcSo,
1575
                                                precSo,
1576
                                                intervalSo,
1577
                                                icSo,
1578
                                                currentApproxErrorSo,
1579
                                                approxAccurSo,
1580
                                                debug=False):
1581
    if debug:
1582
        print "Input arguments:"
1583
        print "Polynomial: ", ; pobyso_autoprint(polySo)
1584
        print "Function: ", ; pobyso_autoprint(funcSo)
1585
        print "Internal precision: ", ; pobyso_autoprint(precSo)
1586
        print "Interval: ", ; pobyso_autoprint(intervalSo)
1587
        print "Current approximation error: ", ; pobyso_autoprint(currentApproxErrorSo)
1588
        print "Requested approxiation error: ", ; pobyso_autoprint(approxAccurSo)
1589
        print "________________"
1590
    approxAccurSa        = pobyso_get_constant_as_rn_so_sa(approxAccurSo)
1591
    currentApproxErrorSa = pobyso_get_constant_as_rn_so_sa(currentApproxErrorSo)
1592
    ## If the current approximation error is too close to the target, there is
1593
    #  no possible gain.
1594
    if currentApproxErrorSa >= approxAccurSa / 2:
1595
        #### Do not return the initial argument but copies: caller may free 
1596
        #    as inutile after call.
1597
        return (sollya_lib_copy_obj(polySo),
1598
                sollya_lib_copy_obj(currentApproxErrorSo))
1599
    degreeSa             = pobyso_polynomial_degree_so_sa(polySo)
1600
    intervalSa           = pobyso_range_to_interval_so_sa(intervalSo)
1601
    
1602
    if debug:
1603
        print "degreeSa              :", degreeSa
1604
        print "intervalSa            :", intervalSa.str(style='brackets')
1605
        print "currentApproxErrorSa  :", currentApproxErrorSa 
1606
        print "approxAccurSa         :", approxAccurSa 
1607
    ### Start with a 0 value expression.
1608
    radiusSa = intervalSa.absolute_diameter() / 2
1609
    if debug:
1610
        print "log2(radius):", RR(radiusSa).log2()
1611
    iterIndex = 0
1612
    while True: 
1613
        resPolySo = pobyso_constant_0_sa_so()
1614
        roundedPolyApproxAccurSa = approxAccurSa / 2
1615
        currentRadiusPowerSa = 1 
1616
        for degree in xrange(0,degreeSa + 1):
1617
            #### At round 0, use the agressive formula. At round 1, run the
1618
            #    proved formula.
1619
            if iterIndex == 0:
1620
                roundingPowerSa = \
1621
                    floor(((currentRadiusPowerSa/roundedPolyApproxAccurSa)*(degree+1)).log2())
1622
            else:
1623
                roundingPowerSa = \
1624
                    floor(((currentRadiusPowerSa/roundedPolyApproxAccurSa)*(degreeSa+1)).log2())
1625
            ## Under extreme conditions the above formulas can evaluate under 2, which is the
1626
            #  minimal precision of an MPFR number.
1627
            if roundingPowerSa < 2:
1628
                roundingPowerSa = 2
1629
            if debug:
1630
                print "roundedPolyApproxAccurSa", roundedPolyApproxAccurSa
1631
                print "currentRadiusPowerSa", currentRadiusPowerSa
1632
                print "Current rounding exponent:", roundingPowerSa
1633
            currentRadiusPowerSa *= radiusSa
1634
            index1So = pobyso_constant_from_int_sa_so(degree)
1635
            index2So = pobyso_constant_from_int_sa_so(degree)
1636
            ### Create a monomial with:
1637
            #   - the coefficient in the initial monomial at the current degrree;
1638
            #   - the current exponent;
1639
            #   - the free variable.
1640
            cmonSo  = \
1641
                sollya_lib_build_function_mul(sollya_lib_coeff(polySo, index1So),
1642
                                              sollya_lib_build_function_pow( \
1643
                                                sollya_lib_build_function_free_variable(), \
1644
                                                index2So))
1645
            roundingPowerSo = pobyso_constant_from_int_sa_so(roundingPowerSa)
1646
            cmonrSo = pobyso_round_coefficients_single_so_so(cmonSo, roundingPowerSo)
1647
            sollya_lib_clear_obj(cmonSo)
1648
            ### Add to the result polynomial.
1649
            resPolySo = sollya_lib_build_function_add(resPolySo,
1650
                                                      cmonrSo)
1651
        # End for.
1652
        ### Check the new polynomial.
1653
        freeVarSo     = sollya_lib_build_function_free_variable()
1654
        changeVarSo   = sollya_lib_sub(freeVarSo, icSo)
1655
        resPolyCvSo   = sollya_lib_evaluate(resPolySo, changeVarSo)
1656
        errFuncSo     = sollya_lib_build_function_sub(sollya_lib_copy_obj(funcSo), 
1657
                                                      resPolyCvSo)
1658
        infNormSo = sollya_lib_dirtyinfnorm(errFuncSo, intervalSo)
1659
        ### This also clears resPolyCvSo.
1660
        sollya_lib_clear_obj(errFuncSo)
1661
        cerrSa    = pobyso_get_constant_as_rn_so_sa(infNormSo)
1662
        if debug:
1663
            print "Error of the new polynomial:", cerrSa
1664
        ### If at round 1, return the initial polynomial error. This should
1665
        #   never happen since the rounding algorithm is proved. But some 
1666
        #   circumstances may break it (e.g. internal precision of tools).
1667
        if cerrSa > approxAccurSa:
1668
            if iterIndex > 0: # Round 1 and beyond.
1669
                sollya_lib_clear_obj(resPolySo)
1670
                sollya_lib_clear_obj(infNormSo)
1671
                #### Do not return the arguments but copies: the caller may free
1672
                #    free the former as inutile after call.
1673
                return (sollya_lib_copy_obj(polySo), 
1674
                        sollya_lib_copy_obj(currentApproxErrorSo))
1675
            else: # Round 0, got round 1
1676
                sollya_lib_clear_obj(resPolySo)
1677
                sollya_lib_clear_obj(infNormSo)
1678
                iterIndex += 1
1679
                continue
1680
        ### If get here it is because cerrSa <= approxAccurSa
1681
        ### Approximation error of the new polynomial is acceptable.
1682
        return (resPolySo, infNormSo)
1683
    # End while True
1684
# End pobyso_round_coefficients_progressive_so_so
1685
    
1686
def pobyso_round_coefficients_single_so_so(polySo, precSo):
1687
    """
1688
    Create a rounded coefficients polynomial from polynomial argument to
1689
    the number of bits in size argument.
1690
    All coefficients are set to the same precision.
1691
    """
1692
    ## TODO: check arguments.
1693
    endEllipListSo = pobyso_build_end_elliptic_list_so_so(precSo)
1694
    polySo = sollya_lib_roundcoefficients(polySo, endEllipListSo, None)
1695
    sollya_lib_clear_obj(endEllipListSo)
1696
    #sollya_lib_clear_obj(endEllipListSo)
1697
    return polySo
1698
    
1699
# End pobyso_round_coefficients_single_so_so
1700

    
1701
def pobyso_set_canonical_off():
1702
    sollya_lib_set_canonical(sollya_lib_off())
1703

    
1704
def pobyso_set_canonical_on():
1705
    sollya_lib_set_canonical(sollya_lib_on())
1706

    
1707
def pobyso_set_prec(p):
1708
    """ Legacy function. See pobyso_set_prec_sa_so. """
1709
    pobyso_set_prec_sa_so(p)
1710

    
1711
def pobyso_set_prec_sa_so(p):
1712
    #a = c_int(p)
1713
    #precSo = c_void_p(sollya_lib_constant_from_int(a))
1714
    #precSo = sollya_lib_constant_from_int(a)
1715
    precSo = pobyso_constant_from_int_sa_so(p)
1716
    sollya_lib_set_prec(precSo)
1717
    sollya_lib_clear_obj(precSo)
1718
# End pobyso_set_prec_sa_so.
1719

    
1720
def pobyso_set_prec_so_so(newPrecSo):
1721
    sollya_lib_set_prec(newPrecSo)
1722
# End pobyso_set_prec_so_so.
1723

    
1724
def pobyso_inf_so_so(intervalSo):
1725
    """
1726
    Very thin wrapper around sollya_lib_inf().
1727
    """
1728
    return sollya_lib_inf(intervalSo)
1729
# End pobyso_inf_so_so.
1730
    
1731
def pobyso_supnorm_so_so(polySo, funcSo, intervalSo, errorTypeSo = None,\
1732
                         accuracySo = None):
1733
    """
1734
    Computes the supnorm of the approximation error between the given 
1735
    polynomial and function.
1736
    errorTypeSo defaults to "absolute".
1737
    accuracySo defaults to 2^(-40).
1738
    """
1739
    if errorTypeSo is None:
1740
        errorTypeSo = sollya_lib_absolute(None)
1741
        errorTypeIsNone = True
1742
    else:
1743
        errorTypeIsNone = False
1744
    #
1745
    if accuracySo is None:
1746
        # Notice the **: we are in Pythonland!
1747
        accuracySo = pobyso_constant_sa_so(RR(2**(-40)))
1748
        accuracyIsNone = True
1749
    else:
1750
        accuracyIsNone = False
1751
    pobyso_autoprint(accuracySo)
1752
    resultSo = \
1753
        sollya_lib_supnorm(polySo, funcSo, intervalSo, errorTypeSo, \
1754
                              accuracySo)
1755
    if errorTypeIsNone:
1756
        sollya_lib_clear_obj(errorTypeSo)
1757
    if accuracyIsNone:
1758
        sollya_lib_clear_obj(accuracySo)
1759
    return resultSo
1760
# End pobyso_supnorm_so_so
1761

    
1762
def pobyso_taylor_expansion_no_change_var_so_so(functionSo, 
1763
                                                degreeSo, 
1764
                                                rangeSo,
1765
                                                errorTypeSo=None, 
1766
                                                sollyaPrecSo=None):
1767
    """
1768
    Compute the Taylor expansion without the variable change
1769
    x -> x-intervalCenter.
1770
    """
1771
    # Change internal Sollya precision, if needed.
1772
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
1773
    sollyaPrecChanged = False
1774
    if sollyaPrecSo is None:
1775
        pass
1776
    else:
1777
        sollya_lib_set_prec(sollyaPrecSo)
1778
        sollyaPrecChanged = True
1779
    # Error type stuff: default to absolute.
1780
    if errorTypeSo is None:
1781
        errorTypeIsNone = True
1782
        errorTypeSo = sollya_lib_absolute(None)
1783
    else:
1784
        errorTypeIsNone = False
1785
    intervalCenterSo = sollya_lib_mid(rangeSo, None)
1786
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo,
1787
                                         intervalCenterSo,
1788
                                         rangeSo, errorTypeSo, None)
1789
    # taylorFormListSaSo is a Python list of Sollya objects references that 
1790
    # are copies of the elements of taylorFormSo.
1791
    # pobyso_get_list_elements_so_so clears taylorFormSo.
1792
    (taylorFormListSaSo, numElementsSa, isEndEllipticSa) = \
1793
        pobyso_get_list_elements_so_so(taylorFormSo)
1794
    polySo = sollya_lib_copy_obj(taylorFormListSaSo[0])
1795
    #print "Num elements:", numElementsSa
1796
    sollya_lib_clear_obj(taylorFormSo)
1797
    #polySo = taylorFormListSaSo[0]
1798
    #errorRangeSo = sollya_lib_copy_obj(taylorFormListSaSo[2])
1799
    errorRangeSo = taylorFormListSaSo[2]
1800
    # No copy_obj needed here: a new objects are created.
1801
    maxErrorSo    = sollya_lib_sup(errorRangeSo)
1802
    minErrorSo    = sollya_lib_inf(errorRangeSo)
1803
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
1804
    absMinErrorSo = sollya_lib_abs(minErrorSo)
1805
    sollya_lib_clear_obj(maxErrorSo)
1806
    sollya_lib_clear_obj(minErrorSo)
1807
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
1808
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
1809
    # If changed, reset the Sollya working precision.
1810
    if sollyaPrecChanged:
1811
        sollya_lib_set_prec(initialSollyaPrecSo)
1812
    sollya_lib_clear_obj(initialSollyaPrecSo)
1813
    if errorTypeIsNone:
1814
        sollya_lib_clear_obj(errorTypeSo)
1815
    pobyso_clear_taylorform_sa_so(taylorFormListSaSo)
1816
    if absMaxErrorSa > absMinErrorSa:
1817
        sollya_lib_clear_obj(absMinErrorSo)
1818
        return (polySo, intervalCenterSo, absMaxErrorSo)
1819
    else:
1820
        sollya_lib_clear_obj(absMaxErrorSo)
1821
        return (polySo, intervalCenterSo, absMinErrorSo)
1822
# end pobyso_taylor_expansion_no_change_var_so_so
1823

    
1824
def pobyso_taylor_expansion_with_change_var_so_so(functionSo, degreeSo, \
1825
                                                  rangeSo, \
1826
                                                  errorTypeSo=None, \
1827
                                                  sollyaPrecSo=None):
1828
    """
1829
    Compute the Taylor expansion with the variable change
1830
    x -> (x-intervalCenter) included.
1831
    """
1832
    # Change Sollya internal precision, if need.
1833
    sollyaPrecChanged = False
1834
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_sos_sa()
1835
    if sollyaPrecSo is None:
1836
        pass
1837
    else:
1838
        sollya_lib_set_prec(sollyaPrecSo)
1839
        sollyaPrecChanged = True
1840
    #
1841
    # Error type stuff: default to absolute.
1842
    if errorTypeSo is None:
1843
        errorTypeIsNone = True
1844
        errorTypeSo = sollya_lib_absolute(None)
1845
    else:
1846
        errorTypeIsNone = False
1847
    intervalCenterSo = sollya_lib_mid(rangeSo)
1848
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo, \
1849
                                         intervalCenterSo, \
1850
                                         rangeSo, errorTypeSo, None)
1851
    # taylorFormListSaSo is a Python list of Sollya objects references that 
1852
    # are copies of the elements of taylorFormSo.
1853
    # pobyso_get_list_elements_so_so clears taylorFormSo.
1854
    (taylorFormListSo, numElements, isEndElliptic) = \
1855
        pobyso_get_list_elements_so_so(taylorFormSo)
1856
    polySo = taylorFormListSo[0]
1857
    errorRangeSo = taylorFormListSo[2]
1858
    maxErrorSo    = sollya_lib_sup(errorRangeSo)
1859
    minErrorSo    = sollya_lib_inf(errorRangeSo)
1860
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
1861
    absMinErrorSo = sollya_lib_abs(minErrorSo)
1862
    sollya_lib_clear_obj(maxErrorSo)
1863
    sollya_lib_clear_obj(minErrorSo)
1864
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
1865
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
1866
    changeVarExpSo = sollya_lib_build_function_sub(\
1867
                       sollya_lib_build_function_free_variable(),\
1868
                       sollya_lib_copy_obj(intervalCenterSo))
1869
    polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpSo)
1870
    sollya_lib_clear_obj(polySo) 
1871
    sollya_lib_clear_obj(changeVarExpSo)
1872
    # If changed, reset the Sollya working precision.
1873
    if sollyaPrecChanged:
1874
        sollya_lib_set_prec(initialSollyaPrecSo)
1875
    sollya_lib_clear_obj(initialSollyaPrecSo)
1876
    if errorTypeIsNone:
1877
        sollya_lib_clear_obj(errorTypeSo)
1878
    sollya_lib_clear_obj(taylorFormSo)
1879
    # Do not clear maxErrorSo.
1880
    if absMaxErrorSa > absMinErrorSa:
1881
        sollya_lib_clear_obj(absMinErrorSo)
1882
        return((polyVarChangedSo, intervalCenterSo, absMaxErrorSo))
1883
    else:
1884
        sollya_lib_clear_obj(absMaxErrorSo)
1885
        return((polyVarChangedSo, intervalCenterSo, absMinErrorSo))
1886
# end pobyso_taylor_expansion_with_change_var_so_so
1887

    
1888
def pobyso_taylor(function, degree, point):
1889
    """ Legacy function. See pobysoTaylor_so_so. """
1890
    return(pobyso_taylor_so_so(function, degree, point))
1891

    
1892
def pobyso_taylor_so_so(functionSo, degreeSo, pointSo):
1893
    return(sollya_lib_taylor(functionSo, degreeSo, pointSo))
1894
    
1895
def pobyso_taylorform(function, degree, point = None, 
1896
                      interval = None, errorType=None):
1897
    """ Legacy function. See pobyso_taylorform_sa_sa;"""
1898
    
1899
def pobyso_taylorform_sa_sa(functionSa, \
1900
                            degreeSa, \
1901
                            pointSa, \
1902
                            intervalSa=None, \
1903
                            errorTypeSa=None, \
1904
                            precisionSa=None):
1905
    """
1906
    Compute the Taylor form of 'degreeSa' for 'functionSa' at 'pointSa' 
1907
    for 'intervalSa' with 'errorTypeSa' (a string) using 'precisionSa'. 
1908
    point: must be a Real or a Real interval.
1909
    return the Taylor form as an array
1910
    TODO: take care of the interval and of the point when it is an interval;
1911
          when errorType is not None;
1912
          take care of the other elements of the Taylor form (coefficients 
1913
          errors and delta.
1914
    """
1915
    # Absolute as the default error.
1916
    if errorTypeSa is None:
1917
        errorTypeSo = sollya_lib_absolute()
1918
    elif errorTypeSa == "relative":
1919
        errorTypeSo = sollya_lib_relative()
1920
    elif errortypeSa == "absolute":
1921
        errorTypeSo = sollya_lib_absolute()
1922
    else:
1923
        # No clean up needed.
1924
        return None
1925
    # Global precision stuff
1926
    sollyaPrecisionChangedSa = False
1927
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
1928
    if precisionSa is None:
1929
        precSa = initialSollyaPrecSa
1930
    else:
1931
        if precSa > initialSollyaPrecSa:
1932
            if precSa <= 2:
1933
                print inspect.stack()[0][3], ":precision change <= 2 requested."
1934
            pobyso_set_prec_sa_so(precSa)
1935
            sollyaPrecisionChangedSa = True
1936
    #        
1937
    if len(functionSa.variables()) > 0:
1938
        varSa = functionSa.variables()[0]
1939
        pobyso_name_free_variable_sa_so(str(varSa))
1940
    # In any case (point or interval) the parent of pointSa has a precision
1941
    # method.
1942
    pointPrecSa = pointSa.parent().precision()
1943
    if precSa > pointPrecSa:
1944
        pointPrecSa = precSa
1945
    # In any case (point or interval) pointSa has a base_ring() method.
1946
    pointBaseRingString = str(pointSa.base_ring())
1947
    if re.search('Interval', pointBaseRingString) is None: # Point
1948
        pointSo = pobyso_constant_sa_so(pointSa, pointPrecSa)
1949
    else: # Interval.
1950
        pointSo = pobyso_interval_to_range_sa_so(pointSa, pointPrecSa)
1951
    # Sollyafy the function.
1952
    functionSo = pobyso_parse_string_sa_so(functionSa._assume_str().replace('_SAGE_VAR_', ''))
1953
    if sollya_lib_obj_is_error(functionSo):
1954
        print "pobyso_tailorform: function string can't be parsed!"
1955
        return None
1956
    # Sollyafy the degree
1957
    degreeSo = sollya_lib_constant_from_int(int(degreeSa))
1958
    # Sollyafy the point
1959
    # Call Sollya
1960
    taylorFormSo = \
1961
        sollya_lib_taylorform(functionSo, degreeSo, pointSo, errorTypeSo,\
1962
                                         None)
1963
    sollya_lib_clear_obj(functionSo)
1964
    sollya_lib_clear_obj(degreeSo)
1965
    sollya_lib_clear_obj(pointSo)
1966
    sollya_lib_clear_obj(errorTypeSo)
1967
    (tfsAsList, numElements, isEndElliptic) = \
1968
            pobyso_get_list_elements_so_so(taylorFormSo)
1969
    polySo = tfsAsList[0]
1970
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
1971
    polyRealField = RealField(maxPrecision)
1972
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, polyRealField)
1973
    if sollyaPrecisionChangedSa:
1974
        sollya_lib_set_prec(initialSollyaPrecSo)
1975
    sollya_lib_clear_obj(initialSollyaPrecSo)
1976
    polynomialRing = polyRealField[str(varSa)]
1977
    polySa = polynomial(expSa, polynomialRing)
1978
    taylorFormSa = [polySa]
1979
    # Final clean-up.
1980
    sollya_lib_clear_obj(taylorFormSo)
1981
    return(taylorFormSa)
1982
# End pobyso_taylor_form_sa_sa
1983

    
1984
def pobyso_taylorform_so_so(functionSo, degreeSo, pointSo, intervalSo=None, \
1985
                            errorTypeSo=None):
1986
    createdErrorType = False
1987
    if errorTypeSo is None:
1988
        errorTypeSo = sollya_lib_absolute()
1989
        createdErrorType = True
1990
    else:
1991
        #TODO: deal with the other case.
1992
        pass
1993
    if intervalSo is None:
1994
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
1995
                                         errorTypeSo, None)
1996
    else:
1997
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
1998
                                         intervalSo, errorTypeSo, None)
1999
    if createdErrorType:
2000
        sollya_lib_clear_obj(errorTypeSo)
2001
    return resultSo
2002
        
2003

    
2004
def pobyso_univar_polynomial_print_reverse(polySa):
2005
    """ Legacy function. See pobyso_univar_polynomial_print_reverse_sa_sa. """
2006
    return(pobyso_univar_polynomial_print_reverse_sa_sa(polySa))
2007

    
2008
def pobyso_univar_polynomial_print_reverse_sa_sa(polySa):
2009
    """
2010
    Return the string representation of a univariate polynomial with
2011
    monomials ordered in the x^0..x^n order of the monomials.
2012
    Remember: Sage
2013
    """
2014
    polynomialRing = polySa.base_ring()
2015
    # A very expensive solution:
2016
    # -create a fake multivariate polynomial field with only one variable,
2017
    #   specifying a negative lexicographical order;
2018
    mpolynomialRing = PolynomialRing(polynomialRing.base(), \
2019
                                     polynomialRing.variable_name(), \
2020
                                     1, order='neglex')
2021
    # - convert the univariate argument polynomial into a multivariate
2022
    #   version;
2023
    p = mpolynomialRing(polySa)
2024
    # - return the string representation of the converted form.
2025
    # There is no simple str() method defined for p's class.
2026
    return(p.__str__())
2027
#
2028
#print pobyso_get_prec()  
2029
pobyso_set_prec(165)
2030
#print pobyso_get_prec()  
2031
#a=100
2032
#print type(a)
2033
#id(a)
2034
#print "Max arity: ", pobyso_max_arity
2035
#print "Function tripleDouble (43) as a string: ", pobyso_function_type_as_string(43)
2036
#print "Function None (44) as a string: ", pobyso_function_type_as_string(44)
2037
print "...Pobyso check done"