Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 213

Historique | Voir | Annoter | Télécharger (58,22 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_function_sub_so_so(exp1So, exp2So):
134
    return(sollya_lib_build_function_sub(exp1So, exp2So))
135

    
136
def pobyso_change_var_in_function_so_so(funcSo, chvarExpSo):
137
    """
138
    Variable change in a function.
139
    """
140
    return(sollya_lib_evaluate(funcSo,chvarExpSo))
141
# End pobyso_change_var_in_function_so_so     
142

    
143
def pobyso_chebyshevform_so_so(functionSo, degreeSo, intervalSo):
144
    resultSo = sollya_lib_chebyshevform(functionSo, degreeSo, intervalSo)
145
    return(resultSo)
146
# End pobyso_chebyshevform_so_so.
147

    
148
def pobyso_clear_taylorform_sa_so(taylorFormSaSo):
149
    """
150
    This method is necessary to correctly clean up the memory from Taylor forms.
151
    These are made of a Sollya object, a Sollya object list, a Sollya object.
152
    For no clearly understood reason, sollya_lib_clear_object_list crashed 
153
    when applied to the object list.
154
    Here, we decompose it into Sage list of Sollya objects references and we
155
     clear them one by one. 
156
    """
157
    sollya_lib_clear_obj(taylorFormSaSo[0])
158
    (coefficientsErrorsListSaSo, numElementsSa, isEndEllipticSa) = \
159
        pobyso_get_list_elements_so_so(taylorFormSaSo[1])
160
    for element in coefficientsErrorsListSaSo:
161
        sollya_lib_clear_obj(element)
162
    sollya_lib_clear_obj(taylorFormSaSo[1])
163
    sollya_lib_clear_obj(taylorFormSaSo[2])
164
# End pobyso_clear_taylorform_sa_so 
165

    
166
def pobyso_cmp(rnArgSa, cteSo):
167
    """
168
    Compare the MPFR value a RealNumber with that of a Sollya constant.
169
    
170
    Get the value of the Sollya constant into a RealNumber and compare
171
    using MPFR. Could be optimized by working directly with a mpfr_t
172
    for the intermediate number. 
173
    """
174
    # Get the precision of the Sollya constant to build a Sage RealNumber 
175
    # with enough precision.to hold it.
176
    precisionOfCte = c_int(0)
177
    # From the Sollya constant, create a local Sage RealNumber.
178
    sollya_lib_get_prec_of_constant(precisionOfCte, cteSo) 
179
    #print "Precision of constant: ", precisionOfCte
180
    RRRR = RealField(precisionOfCte.value)
181
    rnLocalSa = RRRR(0)
182
    sollya_lib_get_constant(get_rn_value(rnLocalSa), cteSo)
183
    #
184
    ## Compare the Sage RealNumber version of the Sollya constant with rnArg.
185
    return(cmp_rn_value(rnArgSa, rnLocal))
186
# End pobyso_smp
187

    
188
def pobyso_compute_pos_function_abs_val_bounds_sa_sa(funcSa, lowerBoundSa, \
189
                                                     upperBoundSa):
190
    """
191
    TODO: completely rework and test.
192
    """
193
    pobyso = pobyso_name_free_variable_sa_so(funcSa.variables()[0])
194
    funcSo = pobyso_parse_string(funcSa._assume_str().replace('_SAGE_VAR_', ''))
195
    rangeSo = pobyso_range_sa_so(lowerBoundSa, upperBoundSa)
196
    infnormSo = pobyso_infnorm_so_so(funcSo,rangeSo)
197
    # Sollya return the infnorm as an interval.
198
    fMaxSa = pobyso_get_interval_from_range_so_sa(infnormSo)
199
    # Get the top bound and compute the binade top limit.
200
    fMaxUpperBoundSa = fMaxSa.upper()
201
    binadeTopLimitSa = 2**ceil(fMaxUpperBoundSa.log2())
202
    # Put up together the function to use to compute the lower bound.
203
    funcAuxSo = pobyso_parse_string(str(binadeTopLimitSa) +  \
204
                                    '-(' + f._assume_str().replace('_SAGE_VAR_', '') + ')')
205
    pobyso_autoprint(funcAuxSo)
206
    # Clear the Sollya range before a new call to infnorm and issue the call.
207
    sollya_lib_clear_obj(infnormSo)
208
    infnormSo = pobyso_infnorm_so_so(funcAuxSo,rangeSo)
209
    fMinSa = pobyso_get_interval_from_range_so_sa(infnormSo)
210
    sollya_lib_clear_obj(infnormSo)
211
    fMinLowerBoundSa = binadeTopLimitSa - fMinSa.lower()
212
    # Compute the maximum of the precisions of the different bounds.
213
    maxPrecSa = max([fMinLowerBoundSa.parent().precision(), \
214
                     fMaxUpperBoundSa.parent().precision()])
215
    # Create a RealIntervalField and create an interval with the "good" bounds.
216
    RRRI = RealIntervalField(maxPrecSa)
217
    imageIntervalSa = RRRI(fMinLowerBoundSa, fMaxUpperBoundSa)
218
    # Free the unneeded Sollya objects
219
    sollya_lib_clear_obj(funcSo)
220
    sollya_lib_clear_obj(funcAuxSo)
221
    sollya_lib_clear_obj(rangeSo)
222
    return(imageIntervalSa)
223
# End pobyso_compute_pos_function_abs_val_bounds_sa_sa
224

    
225
def pobyso_constant(rnArg):
226
    """ Legacy function. See pobyso_constant_sa_so. """
227
    return(pobyso_constant_sa_so(rnArg))
228
    
229
def pobyso_constant_sa_so(rnArgSa, precisionSa=None):
230
    """
231
    Create a Sollya constant from a Sage RealNumber.
232
    The sollya_lib_constant() function creates a constant
233
    with the same precision as the source.
234
    """
235
    ## Precision stuff. If one wants to change precisions,
236
    #  everything takes place in Sage. That only makes
237
    #  sense if one wants to reduce the precision.
238
    if not precisionSa is None:
239
        RRR = RealField(precisionSa)
240
        rnArgSa = RRR(rnArgSa)
241
    #print rnArgSa, rnArgSa.precision()
242
    # Sollya constant creation takes place here.
243
    return sollya_lib_constant(get_rn_value(rnArgSa))
244
# End pobyso_constant_sa_so
245
     
246
def pobyso_constant_0_sa_so():
247
    """
248
    Obvious.
249
    """
250
    return(pobyso_constant_from_int_sa_so(0))
251

    
252
def pobyso_constant_1():
253
    """
254
    Obvious.
255
    Legacy function. See pobyso_constant_so_so. 
256
    """
257
    return(pobyso_constant_1_sa_so())
258

    
259
def pobyso_constant_1_sa_so():
260
    """
261
    Obvious.
262
    """
263
    return(pobyso_constant_from_int_sa_so(1))
264

    
265
def pobyso_constant_from_int(anInt):
266
    """ Legacy function. See pobyso_constant_from_int_sa_so. """
267
    return(pobyso_constant_from_int_sa_so(anInt))
268

    
269
def pobyso_constant_from_int_sa_so(anInt):
270
    """
271
    Get a Sollya constant from a Sage int.
272
    """
273
    return(sollya_lib_constant_from_int64(long(anInt)))
274

    
275
def pobyso_constant_from_int_so_sa(constSo):
276
    """
277
    Get a Sage int from a Sollya int constant.
278
    Usefull for precision or powers in polynomials.
279
    """
280
    constSa = c_long(0)
281
    sollya_lib_get_constant_as_int64(byref(constSa), constSo)
282
    return(constSa.value)
283
# End pobyso_constant_from_int_so_sa
284

    
285
def pobyso_constant_from_mpq_sa_so(rationalSa):
286
    """
287
    Make a Sollya constant from Sage rational.
288
    The Sollya constant is an unevaluated expression.
289
    Hence no precision argument is needed.
290
    It is better to leave this way since Sollya has its own
291
    optimized evaluation mecanism that tries very hard to
292
    return exact values or at least faithful ones.
293
    """
294
    ratExprSo = \
295
        sollya_lib_constant_from_mpq(sgmp_get_rational_value(rationalSa))
296
    return ratExprSo
297
# End pobyso_constant_from_mpq_sa_so.
298

    
299
def pobyso_constant_sollya_prec_sa_so(rnArgSa):
300
    """
301
    Create a Sollya constant from a Sage RealNumber at the
302
    current precision in Sollya.
303
    """
304
    currentSollyaPrecSa = pobyso_get_prec_so_sa()
305
    return pobyso_constant_sa_so(rnArgSa, currentSollyaPrecSa)
306
# End pobyso_constant_sollya_prec_sa_so
307
     
308
def pobyso_error_so():
309
    return sollya_lib_error(None)
310
# End pobyso_error().
311

    
312
def pobyso_float_poly_sa_so(polySa, precSa = None):
313
    """
314
    Create a Sollya polynomial from a Sage RealField polynomial.
315
    """
316
    ## TODO: filter arguments.
317
    ## Precision. If a precision is given, convert the polynomial
318
    #  into the right polynomial field. If not convert it straight
319
    #  to Sollya.
320
    if not precSa is None:
321
        RRR = RealField(precSa)
322
        ## Create a Sage polynomial in the "right" precision.
323
        P_RRR = RRR[polySa.variables()[0]]
324
        polyFloatSa = P_RRR(polySa)
325
    else:
326
        polyFloatSa = polySa
327
        precSa = polySa.parent().base_ring().precision()
328
    ## Get exponents and coefficients.
329
    exponentsSa = polyFloatSa.exponents()
330
    coefficientsSa = polyFloatSa.coefficients()
331
    ## Build the polynomial.
332
    polySo = None
333
    for coefficientSa, exponentSa in zip(coefficientsSa, exponentsSa):
334
        #print coefficientSa.n(prec=precSa), exponentSa
335
        coefficientSo = \
336
            pobyso_constant_sa_so(coefficientSa)
337
        #pobyso_autoprint(coefficientSo)
338
        exponentSo = \
339
            pobyso_constant_from_int_sa_so(exponentSa)
340
        #pobyso_autoprint(exponentSo)
341
        monomialSo = sollya_lib_build_function_pow(
342
                       sollya_lib_build_function_free_variable(),
343
                       exponentSo)
344
        if polySo is None:
345
            polySo = sollya_lib_build_function_mul(coefficientSo,
346
                                                   monomialSo)
347
        else:
348
            polyTermSo = sollya_lib_build_function_mul(coefficientSo,
349
                                                       monomialSo)
350
            polySo = sollya_lib_build_function_add(polySo, polyTermSo)
351
    return polySo
352
# End pobyso_float_poly_sa_so    
353

    
354
def pobyso_float_poly_so_sa(polySo, realFieldSa=None):
355
    """
356
    Convert a Sollya polynomial into a Sage floating-point polynomial.
357
    We assume that the polynomial is in canonical form.
358
    If no realField is given, a RealField corresponding to the maximum 
359
    precision of the coefficients is internally computed.
360
    The real field is not returned but can be easily retrieved from 
361
    the polynomial itself.
362
    ALGORITHM:
363
    - (optional) compute the RealField of the coefficients;
364
    - convert the Sollya expression into a Sage expression;
365
    - convert the Sage expression into a Sage polynomial
366
    TODO: the canonical thing for the polynomial.
367
    """    
368
    if realFieldSa is None:
369
        expressionPrecSa = pobyso_get_max_prec_of_exp_so_sa(polySo)
370
        realFieldSa      = RealField(expressionPrecSa)
371
    #print "Sollya expression before...",
372
    #pobyso_autoprint(polySo)
373

    
374
    expressionSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo,
375
                                                             realFieldSa)
376
    #print "...Sollya expression after.",
377
    #pobyso_autoprint(polySo)
378
    polyVariableSa = expressionSa.variables()[0]
379
    polyRingSa     = realFieldSa[str(polyVariableSa)]
380
    #print polyRingSa
381
    # Do not use the polynomial(expressionSa, ring=polyRingSa) form!
382
    polynomialSa = polyRingSa(expressionSa)
383
    return polynomialSa
384
# End pobyso_float_poly_so_sa
385

    
386
    
387
def pobyso_function_type_as_string(funcType):
388
    """ Legacy function. See pobyso_function_type_as_string_so_sa. """
389
    return(pobyso_function_type_as_string_so_sa(funcType))
390

    
391
def pobyso_function_type_as_string_so_sa(funcType):
392
    """
393
    Numeric Sollya function codes -> Sage mathematical function names.
394
    Notice that pow -> ^ (a la Sage, not a la Python).
395
    """
396
    if funcType == SOLLYA_BASE_FUNC_ABS:
397
        return "abs"
398
    elif funcType == SOLLYA_BASE_FUNC_ACOS:
399
        return "arccos"
400
    elif funcType == SOLLYA_BASE_FUNC_ACOSH:
401
        return "arccosh"
402
    elif funcType == SOLLYA_BASE_FUNC_ADD:
403
        return "+"
404
    elif funcType == SOLLYA_BASE_FUNC_ASIN:
405
        return "arcsin"
406
    elif funcType == SOLLYA_BASE_FUNC_ASINH:
407
        return "arcsinh"
408
    elif funcType == SOLLYA_BASE_FUNC_ATAN:
409
        return "arctan"
410
    elif funcType == SOLLYA_BASE_FUNC_ATANH:
411
        return "arctanh"
412
    elif funcType == SOLLYA_BASE_FUNC_CEIL:
413
        return "ceil"
414
    elif funcType == SOLLYA_BASE_FUNC_CONSTANT:
415
        return "cte"
416
    elif funcType == SOLLYA_BASE_FUNC_COS:
417
        return "cos"
418
    elif funcType == SOLLYA_BASE_FUNC_COSH:
419
        return "cosh"
420
    elif funcType == SOLLYA_BASE_FUNC_DIV:
421
        return "/"
422
    elif funcType == SOLLYA_BASE_FUNC_DOUBLE:
423
        return "double"
424
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEDOUBLE:
425
        return "doubleDouble"
426
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEEXTENDED:
427
        return "doubleDxtended"
428
    elif funcType == SOLLYA_BASE_FUNC_ERF:
429
        return "erf"
430
    elif funcType == SOLLYA_BASE_FUNC_ERFC:
431
        return "erfc"
432
    elif funcType == SOLLYA_BASE_FUNC_EXP:
433
        return "exp"
434
    elif funcType == SOLLYA_BASE_FUNC_EXP_M1:
435
        return "expm1"
436
    elif funcType == SOLLYA_BASE_FUNC_FLOOR:
437
        return "floor"
438
    elif funcType == SOLLYA_BASE_FUNC_FREE_VARIABLE:
439
        return "freeVariable"
440
    elif funcType == SOLLYA_BASE_FUNC_HALFPRECISION:
441
        return "halfPrecision"
442
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYCONSTANT:
443
        return "libraryConstant"
444
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYFUNCTION:
445
        return "libraryFunction"
446
    elif funcType == SOLLYA_BASE_FUNC_LOG:
447
        return "log"
448
    elif funcType == SOLLYA_BASE_FUNC_LOG_10:
449
        return "log10"
450
    elif funcType == SOLLYA_BASE_FUNC_LOG_1P:
451
        return "log1p"
452
    elif funcType == SOLLYA_BASE_FUNC_LOG_2:
453
        return "log2"
454
    elif funcType == SOLLYA_BASE_FUNC_MUL:
455
        return "*"
456
    elif funcType == SOLLYA_BASE_FUNC_NEARESTINT:
457
        return "round"
458
    elif funcType == SOLLYA_BASE_FUNC_NEG:
459
        return "__neg__"
460
    elif funcType == SOLLYA_BASE_FUNC_PI:
461
        return "pi"
462
    elif funcType == SOLLYA_BASE_FUNC_POW:
463
        return "^"
464
    elif funcType == SOLLYA_BASE_FUNC_PROCEDUREFUNCTION:
465
        return "procedureFunction"
466
    elif funcType == SOLLYA_BASE_FUNC_QUAD:
467
        return "quad"
468
    elif funcType == SOLLYA_BASE_FUNC_SIN:
469
        return "sin"
470
    elif funcType == SOLLYA_BASE_FUNC_SINGLE:
471
        return "single"
472
    elif funcType == SOLLYA_BASE_FUNC_SINH:
473
        return "sinh"
474
    elif funcType == SOLLYA_BASE_FUNC_SQRT:
475
        return "sqrt"
476
    elif funcType == SOLLYA_BASE_FUNC_SUB:
477
        return "-"
478
    elif funcType == SOLLYA_BASE_FUNC_TAN:
479
        return "tan"
480
    elif funcType == SOLLYA_BASE_FUNC_TANH:
481
        return "tanh"
482
    elif funcType == SOLLYA_BASE_FUNC_TRIPLEDOUBLE:
483
        return "tripleDouble"
484
    else:
485
        return None
486

    
487
def pobyso_get_constant(rnArgSa, constSo):
488
    """ Legacy function. See pobyso_get_constant_so_sa. """
489
    return pobyso_get_constant_so_sa(rnArgSa, constSo)
490
# End pobyso_get_constant
491

    
492
def pobyso_get_constant_so_sa(rnArgSa, constSo):
493
    """
494
    Set the value of rnArgSo to the value of constSo in MPFR_RNDN mode.
495
    rnArg must already exist and belong to some RealField.
496
    We assume that constSo points to a Sollya constant.
497
    """
498
    outcome = sollya_lib_get_constant(get_rn_value(rnArgSa), constSo)
499
    if outcome == 0: # Failure because constSo is not a constant expression.
500
        return None
501
    else:
502
        return outcome
503
# End  pobyso_get_constant_so_sa
504
   
505
def pobyso_get_constant_as_rn(ctExpSo):
506
    """ 
507
    Legacy function. See pobyso_get_constant_as_rn_so_sa. 
508
    """ 
509
    return(pobyso_get_constant_as_rn_so_sa(ctExpSo))
510
    
511
def pobyso_get_constant_as_rn_so_sa(constExpSo):
512
    """
513
    Get a Sollya constant as a Sage "real number".
514
    The precision of the floating-point number returned is that of the Sollya
515
    constant.
516
    """
517
    precisionSa  = pobyso_get_prec_of_constant_so_sa(constExpSo)
518
    ## If the expression can not be exactly converted, None is returned.
519
    #  In this case opt for the Sollya current expression.
520
    if precisionSa is None:
521
        precisionSa = pobyso_get_prec_so_sa()
522
    RRRR = RealField(precisionSa)
523
    rnSa = RRRR(0)
524
    outcome = sollya_lib_get_constant(get_rn_value(rnSa), constExpSo)
525
    if outcome == 0:
526
        return None
527
    else:
528
        return rnSa
529
# End pobyso_get_constant_as_rn_so_sa
530

    
531
def pobyso_get_constant_as_rn_with_rf(ctExp, realField):
532
    """ 
533
    Legacy function. See pobyso_get_constant_as_rn_with_rf_so_sa.
534
    """
535
    return pobyso_get_constant_as_rn_with_rf_so_sa(ctExp, realField)
536
# End pobyso_get_constant_as_rn_with_rf
537
    
538
def pobyso_get_constant_as_rn_with_rf_so_sa(ctExpSo, realFieldSa = None):
539
    """
540
    Get a Sollya constant as a Sage "real number".
541
    If no real field is specified, the precision of the floating-point number 
542
    returned is that of the Sollya constant.
543
    Otherwise is is that of the real field. Hence rounding may happen.
544
    """
545
    if realFieldSa is None:
546
        return pobyso_get_constant_as_rn_so_sa(ctExpSo)
547
    rnSa = realFieldSa(0)
548
    outcome = sollya_lib_get_constant(get_rn_value(rnSa), ctExpSo)
549
    if outcome == 0:
550
        return None
551
    else:
552
        return rnSa
553
# End pobyso_get_constant_as_rn_with_rf_so_sa
554

    
555
def pobyso_get_free_variable_name():
556
    """ 
557
    Legacy function. See pobyso_get_free_variable_name_so_sa.
558
    """
559
    return(pobyso_get_free_variable_name_so_sa())
560

    
561
def pobyso_get_free_variable_name_so_sa():
562
    return sollya_lib_get_free_variable_name()
563
    
564
def pobyso_get_function_arity(expressionSo):
565
    """ 
566
    Legacy function. See pobyso_get_function_arity_so_sa.
567
    """
568
    return(pobyso_get_function_arity_so_sa(expressionSo))
569

    
570
def pobyso_get_function_arity_so_sa(expressionSo):
571
    arity = c_int(0)
572
    sollya_lib_get_function_arity(byref(arity),expressionSo)
573
    return int(arity.value)
574

    
575
def pobyso_get_head_function(expressionSo):
576
    """ 
577
    Legacy function. See pobyso_get_head_function_so_sa. 
578
    """
579
    return(pobyso_get_head_function_so_sa(expressionSo)) 
580

    
581
def pobyso_get_head_function_so_sa(expressionSo):
582
    functionType = c_int(0)
583
    sollya_lib_get_head_function(byref(functionType), expressionSo, None)
584
    return int(functionType.value)
585

    
586
def pobyso_get_interval_from_range_so_sa(soRange, realIntervalFieldSa = None ):
587
    """
588
    Return the Sage interval corresponding to the Sollya range argument.
589
    If no reaIntervalField is passed as an argument, the interval bounds are not
590
    rounded: they are elements of RealIntervalField of the "right" precision
591
    to hold all the digits.
592
    """
593
    prec = c_int(0)
594
    if realIntervalFieldSa is None:
595
        retval = sollya_lib_get_prec_of_range(byref(prec), soRange, None)
596
        if retval == 0:
597
            return None
598
        realIntervalFieldSa = RealIntervalField(prec.value)
599
    intervalSa = realIntervalFieldSa(0,0)
600
    retval = \
601
        sollya_lib_get_interval_from_range(get_interval_value(intervalSa),\
602
                                           soRange)
603
    if retval == 0:
604
        return None
605
    return intervalSa
606
# End pobyso_get_interval_from_range_so_sa
607

    
608
def pobyso_get_list_elements(soObj):
609
    """ Legacy function. See pobyso_get_list_elements_so_so. """
610
    return pobyso_get_list_elements_so_so(soObj)
611
 
612
def pobyso_get_list_elements_so_so(objectListSo):
613
    """
614
    Get the Sollya list elements as a Sage/Python array of Sollya objects.
615
    
616
    INPUT:
617
    - objectListSo: a Sollya list of Sollya objects.
618
    
619
    OUTPUT:
620
    - a Sage/Python tuple made of:
621
      - a Sage/Python list of Sollya objects,
622
      - a Sage/Python int holding the number of elements,
623
      - a Sage/Python int stating (!= 0) that the list is end-elliptic.
624
    NOTE::
625
        We recover the addresses of the Sollya object from the list of pointers
626
        returned by sollya_lib_get_list_elements. The list itself is freed.
627
    TODO::
628
        Figure out what to do with numElements since the number of elements
629
        can easily be recovered from the list itself. 
630
        Ditto for isEndElliptic.
631
    """
632
    listAddress = POINTER(c_longlong)()
633
    numElements = c_int(0)
634
    isEndElliptic = c_int(0)
635
    listAsSageList = []
636
    result = sollya_lib_get_list_elements(byref(listAddress),\
637
                                          byref(numElements),\
638
                                          byref(isEndElliptic),\
639
                                          objectListSo)
640
    if result == 0 :
641
        return None
642
    for i in xrange(0, numElements.value, 1):
643
       #listAsSageList.append(sollya_lib_copy_obj(listAddress[i]))
644
       listAsSageList.append(listAddress[i])
645
       # Clear each of the elements returned by Sollya.
646
       #sollya_lib_clear_obj(listAddress[i])
647
    # Free the list itself.   
648
    sollya_lib_free(listAddress)
649
    return (listAsSageList, numElements.value, isEndElliptic.value)
650

    
651
def pobyso_get_max_prec_of_exp(soExp):
652
    """ Legacy function. See pobyso_get_max_prec_of_exp_so_sa. """
653
    return pobyso_get_max_prec_of_exp_so_sa(soExp)
654

    
655
def pobyso_get_max_prec_of_exp_so_sa(expSo):
656
    """
657
    Get the maximum precision used for the numbers in a Sollya expression.
658
    
659
    Arguments:
660
    soExp -- a Sollya expression pointer
661
    Return value:
662
    A Python integer
663
    TODO: 
664
    - error management;
665
    - correctly deal with numerical type such as DOUBLEEXTENDED.
666
    """
667
    maxPrecision = 0
668
    minConstPrec = 0
669
    currentConstPrec = 0
670
    operator = pobyso_get_head_function_so_sa(expSo)
671
    if (operator != SOLLYA_BASE_FUNC_CONSTANT) and \
672
    (operator != SOLLYA_BASE_FUNC_FREE_VARIABLE):
673
        (arity, subexpressions) = pobyso_get_subfunctions_so_sa(expSo)
674
        for i in xrange(arity):
675
            maxPrecisionCandidate = \
676
                pobyso_get_max_prec_of_exp_so_sa(subexpressions[i])
677
            if maxPrecisionCandidate > maxPrecision:
678
                maxPrecision = maxPrecisionCandidate
679
        return maxPrecision
680
    elif operator == SOLLYA_BASE_FUNC_CONSTANT:
681
        #minConstPrec = pobyso_get_min_prec_of_constant_so_sa(expSo)
682
        #currentConstPrec = pobyso_get_min_prec_of_constant_so_sa(soExp)
683
        #print minConstPrec, " - ", currentConstPrec 
684
        return pobyso_get_min_prec_of_constant_so_sa(expSo)
685
    
686
    elif operator == SOLLYA_BASE_FUNC_FREE_VARIABLE:
687
        return 0
688
    else:
689
        print "pobyso_get_max_prec_of_exp_so_sa: unexepected operator."
690
        return 0
691

    
692
def pobyso_get_min_prec_of_constant_so_sa(constExpSo):
693
    """
694
    Get the minimum precision necessary to represent the value of a Sollya
695
    constant.
696
    MPFR_MIN_PREC and powers of 2 are taken into account.
697
    We assume that constExpSo is a pointer to a Sollay constant expression.
698
    """
699
    constExpAsRnSa = pobyso_get_constant_as_rn_so_sa(constExpSo)
700
    return(min_mpfr_size(get_rn_value(constExpAsRnSa)))
701

    
702
def pobyso_get_poly_so_sa(polySo, realFieldSa=None):
703
    """
704
    Convert a Sollya polynomial into a Sage polynomial.
705
    Legacy function. Use pobyso_float_poly_so_sa() instead.
706
    """    
707
    return pobyso_float_poly_so_sa(polySo,realFieldSa)
708
# End pobyso_get_poly_so_sa
709

    
710
def pobyso_get_prec():
711
    """ Legacy function. See pobyso_get_prec_so_sa(). """
712
    return pobyso_get_prec_so_sa()
713

    
714
def pobyso_get_prec_so():
715
    """
716
    Get the current default precision in Sollya.
717
    The return value is a Sollya object.
718
    Usefull when modifying the precision back and forth by avoiding
719
    extra conversions.
720
    """
721
    return sollya_lib_get_prec(None)
722
    
723
def pobyso_get_prec_so_sa():
724
    """
725
    Get the current default precision in Sollya.
726
    The return value is Sage/Python int.
727
    """
728
    precSo = sollya_lib_get_prec(None)
729
    precSa = c_int(0)
730
    sollya_lib_get_constant_as_int(byref(precSa), precSo)
731
    sollya_lib_clear_obj(precSo)
732
    return int(precSa.value)
733
# End pobyso_get_prec_so_sa.
734

    
735
def pobyso_get_prec_so_so_sa():
736
    """
737
    Return the current precision both as a Sollya object and a
738
    Sage integer as hybrid tuple.
739
    To avoid multiple calls for precision manipulations. 
740
    """
741
    precSo = sollya_lib_get_prec(None)
742
    precSa = c_int(0)
743
    sollya_lib_get_constant_as_int(byref(precSa), precSo)
744
    return (precSo, precSa)
745
    
746
def pobyso_get_prec_of_constant(ctExpSo):
747
    """ Legacy function. See pobyso_get_prec_of_constant_so_sa. """
748
    return pobyso_get_prec_of_constant_so_sa(ctExpSo)
749

    
750
def pobyso_get_prec_of_constant_so_sa(ctExpSo):
751
    """
752
    Tries to find a precision to represent ctExpSo without rounding.
753
    If not possible, returns None.
754
    """
755
    prec = c_int(0)
756
    retc = sollya_lib_get_prec_of_constant(byref(prec), ctExpSo, None)
757
    if retc == 0:
758
        return None
759
    return int(prec.value)
760

    
761
def pobyso_get_prec_of_range_so_sa(rangeSo):
762
    """
763
    Returns the number of bits elements of a range are coded with.
764
    """
765
    prec = c_int(0)
766
    retc = sollya_lib_get_prec_of_range(byref(prec), rangeSo, None)
767
    if retc == 0:
768
        return(None)
769
    return int(prec.value)
770
# End pobyso_get_prec_of_range_so_sa()
771

    
772
def pobyso_get_sage_exp_from_sollya_exp(sollyaExpSo, realField = RR):
773
    """ Legacy function. See pobyso_get_sage_exp_from_sollya_exp_so_sa. """
774
    return pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, 
775
                                                     realField = RR)
776

    
777
def pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, realFieldSa = RR):
778
    """
779
    Get a Sage expression from a Sollya expression. 
780
    Currently only tested with polynomials with floating-point coefficients.
781
    Notice that, in the returned polynomial, the exponents are RealNumbers.
782
    """
783
    #pobyso_autoprint(sollyaExp)
784
    operatorSa = pobyso_get_head_function_so_sa(sollyaExpSo)
785
    sollyaLibFreeVariableName = sollya_lib_get_free_variable_name()
786
    ## Get rid of the "_"'s in "_x_", if any.
787
    sollyaLibFreeVariableName = re.sub('_', '', sollyaLibFreeVariableName)
788
    # Constants and the free variable are special cases.
789
    # All other operator are dealt with in the same way.
790
    if (operatorSa != SOLLYA_BASE_FUNC_CONSTANT) and \
791
       (operatorSa != SOLLYA_BASE_FUNC_FREE_VARIABLE):
792
        (aritySa, subexpressionsSa) = pobyso_get_subfunctions_so_sa(sollyaExpSo)
793
        if aritySa == 1:
794
            sageExpSa = eval(pobyso_function_type_as_string_so_sa(operatorSa) + \
795
            "(" + pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], \
796
            realFieldSa) + ")")
797
        elif aritySa == 2:
798
            # We do not get through the preprocessor.
799
            # The "^" operator is then a special case.
800
            if operatorSa == SOLLYA_BASE_FUNC_POW:
801
                operatorAsStringSa = "**"
802
            else:
803
                operatorAsStringSa = \
804
                    pobyso_function_type_as_string_so_sa(operatorSa)
805
            sageExpSa = \
806
              eval("pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], realFieldSa)"\
807
              + " " + operatorAsStringSa + " " + \
808
                   "pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[1], realFieldSa)")
809
        # We do not know yet how to deal with arity >= 3 
810
        # (is there any in Sollya anyway?).
811
        else:
812
            sageExpSa = eval('None')
813
        return sageExpSa
814
    elif operatorSa == SOLLYA_BASE_FUNC_CONSTANT:
815
        #print "This is a constant"
816
        return pobyso_get_constant_as_rn_with_rf_so_sa(sollyaExpSo, realFieldSa)
817
    elif operatorSa == SOLLYA_BASE_FUNC_FREE_VARIABLE:
818
        #print "This is free variable"
819
        return eval(sollyaLibFreeVariableName)
820
    else:
821
        print "Unexpected"
822
        return eval('None')
823
# End pobyso_get_sage_exp_from_sollya_exp_so_sa
824

    
825

    
826
def pobyso_get_subfunctions(expressionSo):
827
    """ Legacy function. See pobyso_get_subfunctions_so_sa. """
828
    return pobyso_get_subfunctions_so_sa(expressionSo) 
829
# End pobyso_get_subfunctions.
830
 
831
def pobyso_get_subfunctions_so_sa(expressionSo):
832
    """
833
    Get the subfunctions of an expression.
834
    Return the number of subfunctions and the list of subfunctions addresses.
835
    S.T.: Could not figure out another way than that ugly list of declarations
836
    to recover the addresses of the subfunctions.
837
    We limit ourselves to arity 8 functions. 
838
    """
839
    subf0 = c_int(0)
840
    subf1 = c_int(0)
841
    subf2 = c_int(0)
842
    subf3 = c_int(0)
843
    subf4 = c_int(0)
844
    subf5 = c_int(0)
845
    subf6 = c_int(0)
846
    subf7 = c_int(0)
847
    subf8 = c_int(0)
848
    arity = c_int(0)
849
    nullPtr = POINTER(c_int)()
850
    sollya_lib_get_subfunctions(expressionSo, byref(arity), \
851
      byref(subf0), byref(subf1), byref(subf2), byref(subf3), \
852
      byref(subf4), byref(subf5),\
853
      byref(subf6), byref(subf7), byref(subf8), nullPtr, None) 
854
#    byref(cast(subfunctions[0], POINTER(c_int))), \
855
#    byref(cast(subfunctions[0], POINTER(c_int))), \
856
#    byref(cast(subfunctions[2], POINTER(c_int))), \
857
#    byref(cast(subfunctions[3], POINTER(c_int))), \
858
#    byref(cast(subfunctions[4], POINTER(c_int))), \
859
#    byref(cast(subfunctions[5], POINTER(c_int))), \
860
#    byref(cast(subfunctions[6], POINTER(c_int))), \
861
#    byref(cast(subfunctions[7], POINTER(c_int))), \
862
#    byref(cast(subfunctions[8], POINTER(c_int))), nullPtr)
863
    subfunctions = [subf0, subf1, subf2, subf3, subf4, subf5, subf6, subf7, \
864
                    subf8]
865
    subs = []
866
    if arity.value > pobyso_max_arity:
867
        return(0,[])
868
    for i in xrange(arity.value):
869
        subs.append(int(subfunctions[i].value))
870
        #print subs[i]
871
    return (int(arity.value), subs)
872
# End pobyso_get_subfunctions_so_sa
873
    
874
def pobyso_guess_degree_sa_sa(functionSa, intervalSa, approxErrorSa, 
875
                              weightSa=None, degreeBoundSa=None):
876
    """
877
    Sa_sa variant of the solly_guessdegree function.
878
    Return 0 if something goes wrong.
879
    """
880
    functionAsStringSa = functionSa._assume_str().replace('_SAGE_VAR_', '')
881
    functionSo = pobyso_parse_string_sa_so(functionAsStringSa)
882
    if pobyso_is_error_so_sa(functionSo):
883
        sollya_lib_clear_obj(functionSo)
884
        return 0
885
    rangeSo = pobyso_interval_to_range_sa_so(intervalSa)
886
    # The approximation error is expected to be a floating point number.
887
    if pobyso_is_floating_point_number_sa_sa(approxErrorSa):
888
        approxErrorSo = pobyso_constant_sa_so(approxErrorSa)
889
    else:
890
        approxErrorSo = pobyso_constant_sa_so(RR(approxErrorSa))
891
    if not weightSa is None:
892
        weightAsStringSa = weightSa._assume_str().replace('_SAGE_VAR_', '')
893
        weightSo = pobyso_parse_string_sa_so(weightAsStringSa)
894
        if pobyso_is_error_so_sa(weightSo):
895
            sollya_lib_clear_obj(functionSo)
896
            sollya_lib_clear_obj(rangeSo)
897
            sollya_lib_clear_obj(approxErrorSo)   
898
            sollya_lib_clear_obj(weightSo)
899
            return 0   
900
    else:
901
        weightSo = None
902
    if not degreeBoundSa is None:
903
        degreeBoundSo = pobyso_constant_from_int_sa_so(degreeBoundSa)
904
    else:
905
        degreeBoundSo = None
906
    guessedDegreeSa = pobyso_guess_degree_so_sa(functionSo,
907
                                                rangeSo,
908
                                                approxErrorSo,
909
                                                weightSo,
910
                                                degreeBoundSo)
911
    sollya_lib_clear_obj(functionSo)
912
    sollya_lib_clear_obj(rangeSo)
913
    sollya_lib_clear_obj(approxErrorSo)
914
    if not weightSo is None:
915
        sollya_lib_clear_obj(weightSo)
916
    if not degreeBoundSo is None:
917
        sollya_lib_clear_obj(degreeBoundSo)
918
    return guessedDegreeSa
919
# End poyso_guess_degree_sa_sa
920

    
921
def pobyso_guess_degree_so_sa(functionSo, rangeSo, errorSo, weightSo=None, \
922
                              degreeBoundSo=None):
923
    """
924
    Thin wrapper around the guessdegree function.
925
    Nevertheless, some precision control stuff has been appended.
926
    """
927
    # Deal with Sollya internal precision issues: if it is too small, 
928
    # compared with the error, increases it to about twice -log2(error).
929
    errorSa = pobyso_get_constant_as_rn_with_rf_so_sa(errorSo)
930
    log2ErrorSa = errorSa.log2()
931
    if log2ErrorSa < 0:
932
        neededPrecisionSa = int(2 * int(-log2ErrorSa) / 64) * 64
933
    else:
934
        neededPrecisionSa = int(2 * int(log2ErrorSa) / 64) * 64
935
    #print "Needed precision:", neededPrecisionSa
936
    currentPrecSa = pobyso_get_prec_so_sa()
937
    if neededPrecisionSa > currentPrecSa:
938
        currentPrecSo = pobyso_get_prec_so()
939
        pobyso_set_prec_sa_so(neededPrecisionSa)
940
    #print "Guessing degree..."
941
    # weightSo and degreeBoundsSo are optional arguments.
942
    # As declared, sollya_lib_guessdegree must take 5 arguments.
943
    if weightSo is None:
944
        degreeRangeSo = sollya_lib_guessdegree(functionSo, rangeSo, errorSo, 
945
                                               0, 0, None)
946
    elif degreeBoundSo is None:
947
        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, \
948
                                                errorSo, weightSo, 0, None)
949
    else:
950
        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, errorSo, \
951
                                                weightSo, degreeBoundSo, None)
952
    #print "...degree guess done."
953
    # Restore internal precision, if applicable.
954
    if neededPrecisionSa > currentPrecSa:
955
        pobyso_set_prec_so_so(currentPrecSo)
956
        sollya_lib_clear_obj(currentPrecSo)
957
    degreeIntervalSa = pobyso_range_to_interval_so_sa(degreeRangeSo)
958
    sollya_lib_clear_obj(degreeRangeSo)
959
    # When ok, both bounds match.
960
    # When the degree bound is too low, the upper bound is the degree
961
    # for which the error can be honored.
962
    # When it really goes wrong, the upper bound is infinity. 
963
    if degreeIntervalSa.lower() == degreeIntervalSa.upper():
964
        return int(degreeIntervalSa.lower())
965
    else:
966
        if degreeIntervalSa.upper().is_infinity():
967
            return None
968
        else:
969
            return int(degreeIntervalSa.upper())
970
    # End pobyso_guess_degree_so_sa
971

    
972
def pobyso_infnorm_so_so(func, interval, file = None, intervalList = None):
973
    print "Do not use this function. User pobyso_supnorm_so_so instead."
974
    return None
975

    
976
def pobyso_interval_to_range_sa_so(intervalSa, precisionSa=None):
977
    if precisionSa is None:
978
        precisionSa = intervalSa.parent().precision()
979
    intervalSo = pobyso_bounds_to_range_sa_so(intervalSa.lower(),\
980
                                              intervalSa.upper(),\
981
                                              precisionSa)
982
    return intervalSo
983
# End pobyso_interval_to_range_sa_so
984

    
985
def pobyso_is_error_so_sa(objSo):
986
    """
987
    Thin wrapper around the sollya_lib_obj_is_error() function.
988
    """
989
    if sollya_lib_obj_is_error(objSo) != 0:
990
        return True
991
    else:
992
        return False
993
# End pobyso_is_error-so_sa
994

    
995
def pobyso_is_floating_point_number_sa_sa(numberSa):
996
    """
997
    Check whether a Sage number is floating point.
998
    Exception stuff added because numbers other than
999
    floating-point ones do not have the is_real() attribute.
1000
    """
1001
    try:
1002
        return numberSa.is_real()
1003
    except AttributeError:
1004
        return False
1005
# End pobyso_is_floating_piont_number_sa_sa
1006

    
1007
def pobyso_lib_init():
1008
    sollya_lib_init(None)
1009

    
1010
def pobyso_lib_close():
1011
    sollya_lib_close(None)
1012
    
1013
def pobyso_name_free_variable(freeVariableNameSa):
1014
    """ Legacy function. See pobyso_name_free_variable_sa_so. """
1015
    pobyso_name_free_variable_sa_so(freeVariableNameSa)
1016

    
1017
def pobyso_name_free_variable_sa_so(freeVariableNameSa):
1018
    """
1019
    Set the free variable name in Sollya from a Sage string.
1020
    """
1021
    sollya_lib_name_free_variable(freeVariableNameSa)
1022

    
1023
def pobyso_parse_string(string):
1024
    """ Legacy function. See pobyso_parse_string_sa_so. """
1025
    return pobyso_parse_string_sa_so(string)
1026
 
1027
def pobyso_parse_string_sa_so(string):
1028
    """
1029
    Get the Sollya expression computed from a Sage string or
1030
    a Sollya error object if parsing failed.
1031
    """
1032
    return sollya_lib_parse_string(string)
1033

    
1034
def pobyso_precision_so_sa(ctExpSo):
1035
    """
1036
    Computes the necessary precision to represent a number.
1037
    If x is not zero, it can be uniquely written as x = m · 2e 
1038
    where m is an odd integer and e is an integer. 
1039
    precision(x) returns the number of bits necessary to write m 
1040
    in binary (i.e. ceil(log2(m))).
1041
    """
1042
    #TODO: take care of the special case: 0, @NaN@, @Inf@
1043
    precisionSo = sollya_lib_precision(ctExpSo)
1044
    precisionSa = pobyso_constant_from_int_so_sa(precisionSo)
1045
    sollya_lib_clear_obj(precisionSo)
1046
    return precisionSa
1047
# End pobyso_precision_so_sa
1048
    
1049
def pobyso_range(rnLowerBound, rnUpperBound):
1050
    """ Legacy function. See pobyso_range_sa_so. """
1051
    return pobyso_range_sa_so(rnLowerBound, rnUpperBound) 
1052

    
1053

    
1054
def pobyso_range_to_interval_so_sa(rangeSo, realIntervalFieldSa = None):
1055
    """
1056
    Get a Sage interval from a Sollya range.
1057
    If no realIntervalField is given as a parameter, the Sage interval
1058
    precision is that of the Sollya range.
1059
    Otherwise, the precision is that of the realIntervalField. In this case
1060
    rounding may happen.
1061
    """
1062
    if realIntervalFieldSa is None:
1063
        precSa = pobyso_get_prec_of_range_so_sa(rangeSo)
1064
        realIntervalFieldSa = RealIntervalField(precSa)
1065
    intervalSa = \
1066
        pobyso_get_interval_from_range_so_sa(rangeSo, realIntervalFieldSa) 
1067
    return intervalSa
1068
# End pobyso_range_to_interval_so_sa
1069

    
1070
def pobyso_rat_poly_sa_so(polySa, precSa = None):
1071
    """
1072
    Create a Sollya polynomial from a Sage rational polynomial.
1073
    """
1074
    ## TODO: filter arguments.
1075
    ## Precision. If no precision is given, use the current precision
1076
    #  of Sollya.
1077
    if precSa is None:
1078
        precSa =  pobyso_get_prec_so_sa()
1079
    #print "Precision:",  precSa
1080
    RRR = RealField(precSa)
1081
    ## Create a Sage polynomial in the "right" precision.
1082
    P_RRR = RRR[polySa.variables()[0]]
1083
    polyFloatSa = P_RRR(polySa)
1084
    ## Make sure no precision is provided: pobyso_float_poly_sa_so will
1085
    #  recover it all by itself and not make an extra conversion.
1086
    return pobyso_float_poly_sa_so(polyFloatSa)
1087
    
1088
# End pobyso_rat_poly_sa_so    
1089
    
1090
def pobyso_remez_canonical_sa_sa(func, \
1091
                                 degree, \
1092
                                 lowerBound, \
1093
                                 upperBound, \
1094
                                 weight = None, \
1095
                                 quality = None):
1096
    """
1097
    All arguments are Sage/Python.
1098
    The functions (func and weight) must be passed as expressions or strings.
1099
    Otherwise the function fails. 
1100
    The return value is a Sage polynomial.
1101
    """
1102
    var('zorglub')    # Dummy variable name for type check only. Type of 
1103
    # zorglub is "symbolic expression".
1104
    polySo = pobyso_remez_canonical_sa_so(func, \
1105
                                 degree, \
1106
                                 lowerBound, \
1107
                                 upperBound, \
1108
                                 weight, \
1109
                                 quality)
1110
    # String test
1111
    if parent(func) == parent("string"):
1112
        functionSa = eval(func)
1113
    # Expression test.
1114
    elif type(func) == type(zorglub):
1115
        functionSa = func
1116
    else:
1117
        return None
1118
    #
1119
    maxPrecision = 0
1120
    if polySo is None:
1121
        return(None)
1122
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
1123
    RRRRSa = RealField(maxPrecision)
1124
    polynomialRingSa = RRRRSa[functionSa.variables()[0]]
1125
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, RRRRSa)
1126
    polySa = polynomial(expSa, polynomialRingSa)
1127
    sollya_lib_clear_obj(polySo)
1128
    return(polySa)
1129
# End pobyso_remez_canonical_sa_sa
1130
    
1131
def pobyso_remez_canonical(func, \
1132
                           degree, \
1133
                           lowerBound, \
1134
                           upperBound, \
1135
                           weight = "1", \
1136
                           quality = None):
1137
    """ Legacy function. See pobyso_remez_canonical_sa_so. """
1138
    return(pobyso_remez_canonical_sa_so(func, \
1139
                                        degree, \
1140
                                        lowerBound, \
1141
                                        upperBound, \
1142
                                        weight, \
1143
                                        quality))
1144
# End pobyso_remez_canonical.
1145

    
1146
def pobyso_remez_canonical_sa_so(func, \
1147
                                 degree, \
1148
                                 lowerBound, \
1149
                                 upperBound, \
1150
                                 weight = None, \
1151
                                 quality = None):
1152
    """
1153
    All arguments are Sage/Python.
1154
    The functions (func and weight) must be passed as expressions or strings.
1155
    Otherwise the function fails. 
1156
    The return value is a pointer to a Sollya function.
1157
    """
1158
    var('zorglub')    # Dummy variable name for type check only. Type of
1159
    # zorglub is "symbolic expression".
1160
    currentVariableNameSa = None
1161
    # The func argument can be of different types (string, 
1162
    # symbolic expression...)
1163
    if parent(func) == parent("string"):
1164
        localFuncSa = eval(func)
1165
        if len(localFuncSa.variables()) > 0:
1166
            currentVariableNameSa = localFuncSa.variables()[0]
1167
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1168
            functionSo = \
1169
              sollya_lib_parse_string(localFuncSa._assume_str().replace('_SAGE_VAR_', ''))
1170
    # Expression test.
1171
    elif type(func) == type(zorglub):
1172
        # Until we are able to translate Sage expressions into Sollya 
1173
        # expressions : parse the string version.
1174
        if len(func.variables()) > 0:
1175
            currentVariableNameSa = func.variables()[0]
1176
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1177
            functionSo = \
1178
              sollya_lib_parse_string(func._assume_str().replace('_SAGE_VAR_', ''))
1179
    else:
1180
        return(None)
1181
    if weight is None: # No weight given -> 1.
1182
        weightSo = pobyso_constant_1_sa_so()
1183
    elif parent(weight) == parent("string"): # Weight given as string: parse it.
1184
        weightSo = sollya_lib_parse_string(func)
1185
    elif type(weight) == type(zorglub): # Weight given as symbolice expression.
1186
        functionSo = \
1187
          sollya_lib_parse_string_sa_so(weight._assume_str().replace('_SAGE_VAR_', ''))
1188
    else:
1189
        return(None)
1190
    degreeSo = pobyso_constant_from_int(degree)
1191
    rangeSo = pobyso_bounds_to_range_sa_so(lowerBound, upperBound)
1192
    if not quality is None:
1193
        qualitySo= pobyso_constant_sa_so(quality)
1194
    else:
1195
        qualitySo = None
1196
        
1197
    remezPolySo = sollya_lib_remez(functionSo, \
1198
                                   degreeSo, \
1199
                                   rangeSo, \
1200
                                   weightSo, \
1201
                                   qualitySo, \
1202
                                   None)
1203
    sollya_lib_clear_obj(functionSo)
1204
    sollya_lib_clear_obj(degreeSo)
1205
    sollya_lib_clear_obj(rangeSo)
1206
    sollya_lib_clear_obj(weightSo)
1207
    if not qualitySo is None:
1208
        sollya_lib_clear_obj(qualitySo)
1209
    return(remezPolySo)
1210
# End pobyso_remez_canonical_sa_so
1211

    
1212
def pobyso_remez_canonical_so_so(funcSo, \
1213
                                 degreeSo, \
1214
                                 rangeSo, \
1215
                                 weightSo = pobyso_constant_1_sa_so(),\
1216
                                 qualitySo = None):
1217
    """
1218
    All arguments are pointers to Sollya objects.
1219
    The return value is a pointer to a Sollya function.
1220
    """
1221
    if not sollya_lib_obj_is_function(funcSo):
1222
        return(None)
1223
    return(sollya_lib_remez(funcSo, degreeSo, rangeSo, weightSo, qualitySo, None))
1224
# End pobyso_remez_canonical_so_so.
1225

    
1226
def pobyso_set_canonical_off():
1227
    sollya_lib_set_canonical(sollya_lib_off())
1228

    
1229
def pobyso_set_canonical_on():
1230
    sollya_lib_set_canonical(sollya_lib_on())
1231

    
1232
def pobyso_set_prec(p):
1233
    """ Legacy function. See pobyso_set_prec_sa_so. """
1234
    pobyso_set_prec_sa_so(p)
1235

    
1236
def pobyso_set_prec_sa_so(p):
1237
    a = c_int(p)
1238
    precSo = c_void_p(sollya_lib_constant_from_int(a))
1239
    sollya_lib_set_prec(precSo, None)
1240

    
1241
def pobyso_set_prec_so_so(newPrecSo):
1242
    sollya_lib_set_prec(newPrecSo, None)
1243

    
1244
def pobyso_supnorm_so_so(polySo, funcSo, intervalSo, errorTypeSo = None,\
1245
                         accuracySo = None):
1246
    """
1247
    Computes the supnorm of the approximation error between the given 
1248
    polynomial and function.
1249
    errorTypeSo defaults to "absolute".
1250
    accuracySo defaults to 2^(-40).
1251
    """
1252
    if errorTypeSo is None:
1253
        errorTypeSo = sollya_lib_absolute(None)
1254
        errorTypeIsNone = True
1255
    else:
1256
        errorTypeIsNone = False
1257
    #
1258
    if accuracySo is None:
1259
        # Notice the **!
1260
        accuracySo = pobyso_constant_sa_so(RR(2**(-40)))
1261
        accuracyIsNone = True
1262
    else:
1263
        accuracyIsNone = False
1264
    pobyso_autoprint(accuracySo)
1265
    resultSo = \
1266
        sollya_lib_supnorm(polySo, funcSo, intervalSo, errorTypeSo, \
1267
                              accuracySo)
1268
    if errorTypeIsNone:
1269
        sollya_lib_clear_obj(errorTypeSo)
1270
    if accuracyIsNone:
1271
        sollya_lib_clear_obj(accuracySo)
1272
    return resultSo
1273
# End pobyso_supnorm_so_so
1274

    
1275
def pobyso_taylor_expansion_no_change_var_so_so(functionSo, 
1276
                                                degreeSo, 
1277
                                                rangeSo,
1278
                                                errorTypeSo=None, 
1279
                                                sollyaPrecSo=None):
1280
    """
1281
    Compute the Taylor expansion without the variable change
1282
    x -> x-intervalCenter.
1283
    """
1284
    # No global change of the working precision.
1285
    if not sollyaPrecSo is None:
1286
        initialPrecSo = sollya_lib_get_prec(None)
1287
        sollya_lib_set_prec(sollyaPrecSo)
1288
    # Error type stuff: default to absolute.
1289
    if errorTypeSo is None:
1290
        errorTypeIsNone = True
1291
        errorTypeSo = sollya_lib_absolute(None)
1292
    else:
1293
        errorTypeIsNone = False
1294
    intervalCenterSo = sollya_lib_mid(rangeSo, None)
1295
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo,
1296
                                         intervalCenterSo,
1297
                                         rangeSo, errorTypeSo, None)
1298
    # taylorFormListSaSo is a Python list of Sollya objects references that 
1299
    # are copies of the elements of taylorFormSo.
1300
    # pobyso_get_list_elements_so_so clears taylorFormSo.
1301
    (taylorFormListSaSo, numElementsSa, isEndEllipticSa) = \
1302
        pobyso_get_list_elements_so_so(taylorFormSo)
1303
    polySo = sollya_lib_copy_obj(taylorFormListSaSo[0])
1304
    #print "Num elements:", numElementsSa
1305
    sollya_lib_clear_obj(taylorFormSo)
1306
    #polySo = taylorFormListSaSo[0]
1307
    #errorRangeSo = sollya_lib_copy_obj(taylorFormListSaSo[2])
1308
    errorRangeSo = taylorFormListSaSo[2]
1309
    # No copy_obj needed here: a new objects are created.
1310
    maxErrorSo    = sollya_lib_sup(errorRangeSo)
1311
    minErrorSo    = sollya_lib_inf(errorRangeSo)
1312
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
1313
    absMinErrorSo = sollya_lib_abs(minErrorSo)
1314
    sollya_lib_clear_obj(maxErrorSo)
1315
    sollya_lib_clear_obj(minErrorSo)
1316
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
1317
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
1318
    # If changed, reset the Sollya working precision.
1319
    if not sollyaPrecSo is None:
1320
        sollya_lib_set_prec(initialPrecSo)
1321
        sollya_lib_clear_obj(initialPrecSo)
1322
    if errorTypeIsNone:
1323
        sollya_lib_clear_obj(errorTypeSo)
1324
    pobyso_clear_taylorform_sa_so(taylorFormListSaSo)
1325
    if absMaxErrorSa > absMinErrorSa:
1326
        sollya_lib_clear_obj(absMinErrorSo)
1327
        return((polySo, intervalCenterSo, absMaxErrorSo))
1328
    else:
1329
        sollya_lib_clear_obj(absMaxErrorSo)
1330
        return((polySo, intervalCenterSo, absMinErrorSo))
1331
# end pobyso_taylor_expansion_no_change_var_so_so
1332

    
1333
def pobyso_taylor_expansion_with_change_var_so_so(functionSo, degreeSo, \
1334
                                                  rangeSo, \
1335
                                                  errorTypeSo=None, \
1336
                                                  sollyaPrecSo=None):
1337
    """
1338
    Compute the Taylor expansion with the variable change
1339
    x -> (x-intervalCenter) included.
1340
    """
1341
    # No global change of the working precision.
1342
    if not sollyaPrecSo is None:
1343
        initialPrecSo = sollya_lib_get_prec(None)
1344
        sollya_lib_set_prec(sollyaPrecSo)
1345
    #
1346
    # Error type stuff: default to absolute.
1347
    if errorTypeSo is None:
1348
        errorTypeIsNone = True
1349
        errorTypeSo = sollya_lib_absolute(None)
1350
    else:
1351
        errorTypeIsNone = False
1352
    intervalCenterSo = sollya_lib_mid(rangeSo)
1353
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo, \
1354
                                         intervalCenterSo, \
1355
                                         rangeSo, errorTypeSo, None)
1356
    # taylorFormListSaSo is a Python list of Sollya objects references that 
1357
    # are copies of the elements of taylorFormSo.
1358
    # pobyso_get_list_elements_so_so clears taylorFormSo.
1359
    (taylorFormListSo, numElements, isEndElliptic) = \
1360
        pobyso_get_list_elements_so_so(taylorFormSo)
1361
    polySo = taylorFormListSo[0]
1362
    errorRangeSo = taylorFormListSo[2]
1363
    maxErrorSo    = sollya_lib_sup(errorRangeSo)
1364
    minErrorSo    = sollya_lib_inf(errorRangeSo)
1365
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
1366
    absMinErrorSo = sollya_lib_abs(minErrorSo)
1367
    sollya_lib_clear_obj(maxErrorSo)
1368
    sollya_lib_clear_obj(minErrorSo)
1369
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
1370
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
1371
    changeVarExpSo = sollya_lib_build_function_sub(\
1372
                       sollya_lib_build_function_free_variable(),\
1373
                       sollya_lib_copy_obj(intervalCenterSo))
1374
    polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpSo)
1375
    sollya_lib_clear_obj(polySo) 
1376
    sollya_lib_clear_obj(changeVarExpSo)
1377
    # If changed, reset the Sollya working precision.
1378
    if not sollyaPrecSo is None:
1379
        sollya_lib_set_prec(initialPrecSo)
1380
        sollya_lib_clear_obj(initialPrecSo)
1381
    if errorTypeIsNone:
1382
        sollya_lib_clear_obj(errorTypeSo)
1383
    sollya_lib_clear_obj(taylorFormSo)
1384
    # Do not clear maxErrorSo.
1385
    if absMaxErrorSa > absMinErrorSa:
1386
        sollya_lib_clear_obj(absMinErrorSo)
1387
        return((polyVarChangedSo, intervalCenterSo, absMaxErrorSo))
1388
    else:
1389
        sollya_lib_clear_obj(absMaxErrorSo)
1390
        return((polyVarChangedSo, intervalCenterSo, absMinErrorSo))
1391
# end pobyso_taylor_expansion_with_change_var_so_so
1392

    
1393
def pobyso_taylor(function, degree, point):
1394
    """ Legacy function. See pobysoTaylor_so_so. """
1395
    return(pobyso_taylor_so_so(function, degree, point))
1396

    
1397
def pobyso_taylor_so_so(functionSo, degreeSo, pointSo):
1398
    return(sollya_lib_taylor(functionSo, degreeSo, pointSo))
1399
    
1400
def pobyso_taylorform(function, degree, point = None, 
1401
                      interval = None, errorType=None):
1402
    """ Legacy function. See pobyso_taylorform_sa_sa;"""
1403
    
1404
def pobyso_taylorform_sa_sa(functionSa, \
1405
                            degreeSa, \
1406
                            pointSa, \
1407
                            intervalSa=None, \
1408
                            errorTypeSa=None, \
1409
                            precisionSa=None):
1410
    """
1411
    Compute the Taylor form of 'degreeSa' for 'functionSa' at 'pointSa' 
1412
    for 'intervalSa' with 'errorTypeSa' (a string) using 'precisionSa'. 
1413
    point: must be a Real or a Real interval.
1414
    return the Taylor form as an array
1415
    TODO: take care of the interval and of the point when it is an interval;
1416
          when errorType is not None;
1417
          take care of the other elements of the Taylor form (coefficients 
1418
          errors and delta.
1419
    """
1420
    # Absolute as the default error.
1421
    if errorTypeSa is None:
1422
        errorTypeSo = sollya_lib_absolute()
1423
    elif errorTypeSa == "relative":
1424
        errorTypeSo = sollya_lib_relative()
1425
    elif errortypeSa == "absolute":
1426
        errorTypeSo = sollya_lib_absolute()
1427
    else:
1428
        # No clean up needed.
1429
        return None
1430
    # Global precision stuff
1431
    precisionChangedSa = False
1432
    currentSollyaPrecSo = pobyso_get_prec_so()
1433
    currentSollyaPrecSa = pobyso_constant_from_int_so_sa(currentSollyaPrecSo)
1434
    if not precisionSa is None:
1435
        if precisionSa > currentSollyaPrecSa:
1436
            pobyso_set_prec_sa_so(precisionSa)
1437
            precisionChangedSa = True
1438
            
1439
    if len(functionSa.variables()) > 0:
1440
        varSa = functionSa.variables()[0]
1441
        pobyso_name_free_variable_sa_so(str(varSa))
1442
    # In any case (point or interval) the parent of pointSa has a precision
1443
    # method.
1444
    pointPrecSa = pointSa.parent().precision()
1445
    if precisionSa > pointPrecSa:
1446
        pointPrecSa = precisionSa
1447
    # In any case (point or interval) pointSa has a base_ring() method.
1448
    pointBaseRingString = str(pointSa.base_ring())
1449
    if re.search('Interval', pointBaseRingString) is None: # Point
1450
        pointSo = pobyso_constant_sa_so(pointSa, pointPrecSa)
1451
    else: # Interval.
1452
        pointSo = pobyso_interval_to_range_sa_so(pointSa, pointPrecSa)
1453
    # Sollyafy the function.
1454
    functionSo = pobyso_parse_string_sa_so(functionSa._assume_str().replace('_SAGE_VAR_', ''))
1455
    if sollya_lib_obj_is_error(functionSo):
1456
        print "pobyso_tailorform: function string can't be parsed!"
1457
        return None
1458
    # Sollyafy the degree
1459
    degreeSo = sollya_lib_constant_from_int(int(degreeSa))
1460
    # Sollyafy the point
1461
    # Call Sollya
1462
    taylorFormSo = \
1463
        sollya_lib_taylorform(functionSo, degreeSo, pointSo, errorTypeSo,\
1464
                                         None)
1465
    sollya_lib_clear_obj(functionSo)
1466
    sollya_lib_clear_obj(degreeSo)
1467
    sollya_lib_clear_obj(pointSo)
1468
    sollya_lib_clear_obj(errorTypeSo)
1469
    (tfsAsList, numElements, isEndElliptic) = \
1470
            pobyso_get_list_elements_so_so(taylorFormSo)
1471
    polySo = tfsAsList[0]
1472
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
1473
    polyRealField = RealField(maxPrecision)
1474
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, polyRealField)
1475
    if precisionChangedSa:
1476
        sollya_lib_set_prec(currentSollyaPrecSo)
1477
        sollya_lib_clear_obj(currentSollyaPrecSo)
1478
    polynomialRing = polyRealField[str(varSa)]
1479
    polySa = polynomial(expSa, polynomialRing)
1480
    taylorFormSa = [polySa]
1481
    # Final clean-up.
1482
    sollya_lib_clear_obj(taylorFormSo)
1483
    return(taylorFormSa)
1484
# End pobyso_taylor_form_sa_sa
1485

    
1486
def pobyso_taylorform_so_so(functionSo, degreeSo, pointSo, intervalSo=None, \
1487
                            errorTypeSo=None):
1488
    createdErrorType = False
1489
    if errorTypeSo is None:
1490
        errorTypeSo = sollya_lib_absolute()
1491
        createdErrorType = True
1492
    else:
1493
        #TODO: deal with the other case.
1494
        pass
1495
    if intervalSo is None:
1496
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
1497
                                         errorTypeSo, None)
1498
    else:
1499
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
1500
                                         intervalSo, errorTypeSo, None)
1501
    if createdErrorType:
1502
        sollya_lib_clear_obj(errorTypeSo)
1503
    return(resultSo)
1504
        
1505

    
1506
def pobyso_univar_polynomial_print_reverse(polySa):
1507
    """ Legacy function. See pobyso_univar_polynomial_print_reverse_sa_sa. """
1508
    return(pobyso_univar_polynomial_print_reverse_sa_sa(polySa))
1509

    
1510
def pobyso_univar_polynomial_print_reverse_sa_sa(polySa):
1511
    """
1512
    Return the string representation of a univariate polynomial with
1513
    monomials ordered in the x^0..x^n order of the monomials.
1514
    Remember: Sage
1515
    """
1516
    polynomialRing = polySa.base_ring()
1517
    # A very expensive solution:
1518
    # -create a fake multivariate polynomial field with only one variable,
1519
    #   specifying a negative lexicographical order;
1520
    mpolynomialRing = PolynomialRing(polynomialRing.base(), \
1521
                                     polynomialRing.variable_name(), \
1522
                                     1, order='neglex')
1523
    # - convert the univariate argument polynomial into a multivariate
1524
    #   version;
1525
    p = mpolynomialRing(polySa)
1526
    # - return the string representation of the converted form.
1527
    # There is no simple str() method defined for p's class.
1528
    return(p.__str__())
1529
#
1530
print pobyso_get_prec()  
1531
pobyso_set_prec(165)
1532
print pobyso_get_prec()  
1533
a=100
1534
print type(a)
1535
id(a)
1536
print "Max arity: ", pobyso_max_arity
1537
print "Function tripleDouble (43) as a string: ", pobyso_function_type_as_string(43)
1538
print "Function None (44) as a string: ", pobyso_function_type_as_string(44)
1539
print "...Pobyso check done"