Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 195

Historique | Voir | Annoter | Télécharger (52,88 ko)

1
"""
2
Actual functions to use in Sage
3
ST 2012-11-13
4

5
Command line syntax:
6
  use from Sage (via the "load" or the "attach" commands)
7

8
pobyso functions come in five flavors: 
9
- the _so_so (arguments and returned objects are pointers to Sollya objects, 
10
  includes the void function and the no arguments function that return a 
11
  pointer to a Sollya object);
12
- the _so_sa (argument are pointers to Sollya objects, returned objects are
13
  Sage/Python objects or, more generally, information is transfered from the
14
  Sollya world to Sage/Python world; e.g. functions without arguments that
15
  return a Sage/Python object);
16
- the _sa_so (arguments are Sage/Python objects, returned objects are 
17
  pointers to Sollya objects);
18
- the sa_sa (arguments and returned objects are all Sage/Python objects);
19
- a catch all flavor, without any suffix, (e. g. functions that have no argument 
20
  nor return value).
21
This classification is not always very strict. Conversion functions from Sollya
22
to Sage/Python are sometimes decorated with Sage/Python arguments to set
23
the precision. These functions remain in the so_sa category.
24
NOTES:
25
Reported errors in Eclipse come from the calls to
26
the Sollya library
27

28
ToDo (among other things): 
29
 -memory management.
30
"""
31
from ctypes import *
32
import re
33
from sage.symbolic.expression_conversions import polynomial
34
from sage.symbolic.expression_conversions import PolynomialConverter
35
"""
36
Create the equivalent to an enum for the Sollya function types.
37
"""
38
(SOLLYA_BASE_FUNC_ABS,
39
SOLLYA_BASE_FUNC_ACOS,
40
    SOLLYA_BASE_FUNC_ACOSH,
41
    SOLLYA_BASE_FUNC_ADD,
42
    SOLLYA_BASE_FUNC_ASIN,
43
    SOLLYA_BASE_FUNC_ASINH,
44
    SOLLYA_BASE_FUNC_ATAN,
45
    SOLLYA_BASE_FUNC_ATANH,
46
    SOLLYA_BASE_FUNC_CEIL,
47
    SOLLYA_BASE_FUNC_CONSTANT,
48
    SOLLYA_BASE_FUNC_COS,
49
    SOLLYA_BASE_FUNC_COSH,
50
    SOLLYA_BASE_FUNC_DIV,
51
    SOLLYA_BASE_FUNC_DOUBLE,
52
    SOLLYA_BASE_FUNC_DOUBLEDOUBLE,
53
    SOLLYA_BASE_FUNC_DOUBLEEXTENDED,
54
    SOLLYA_BASE_FUNC_ERF,
55
    SOLLYA_BASE_FUNC_ERFC,
56
    SOLLYA_BASE_FUNC_EXP,
57
    SOLLYA_BASE_FUNC_EXP_M1,
58
    SOLLYA_BASE_FUNC_FLOOR,
59
    SOLLYA_BASE_FUNC_FREE_VARIABLE,
60
    SOLLYA_BASE_FUNC_HALFPRECISION,
61
    SOLLYA_BASE_FUNC_LIBRARYCONSTANT,
62
    SOLLYA_BASE_FUNC_LIBRARYFUNCTION,
63
    SOLLYA_BASE_FUNC_LOG,
64
    SOLLYA_BASE_FUNC_LOG_10,
65
    SOLLYA_BASE_FUNC_LOG_1P,
66
    SOLLYA_BASE_FUNC_LOG_2,
67
    SOLLYA_BASE_FUNC_MUL,
68
    SOLLYA_BASE_FUNC_NEARESTINT,
69
    SOLLYA_BASE_FUNC_NEG,
70
    SOLLYA_BASE_FUNC_PI,
71
    SOLLYA_BASE_FUNC_POW,
72
    SOLLYA_BASE_FUNC_PROCEDUREFUNCTION,
73
    SOLLYA_BASE_FUNC_QUAD,
74
    SOLLYA_BASE_FUNC_SIN,
75
    SOLLYA_BASE_FUNC_SINGLE,
76
    SOLLYA_BASE_FUNC_SINH,
77
    SOLLYA_BASE_FUNC_SQRT,
78
    SOLLYA_BASE_FUNC_SUB,
79
    SOLLYA_BASE_FUNC_TAN,
80
    SOLLYA_BASE_FUNC_TANH,
81
SOLLYA_BASE_FUNC_TRIPLEDOUBLE) = map(int,xrange(44))
82
print "\nSuperficial pobyso check..."    
83
print "First constant - SOLLYA_BASE_FUNC_ABS: ", SOLLYA_BASE_FUNC_ABS
84
print "Last constant  - SOLLYA_BASE_FUNC_TRIPLEDOUBLE: ", SOLLYA_BASE_FUNC_TRIPLEDOUBLE
85

    
86
pobyso_max_arity = 9
87

    
88
def pobyso_absolute_so_so():
89
    return(sollya_lib_absolute(None))
90

    
91
def pobyso_autoprint(arg):
92
    sollya_lib_autoprint(arg,None)
93

    
94
def pobyso_autoprint_so_so(arg):
95
    sollya_lib_autoprint(arg,None)
96
    
97
def pobyso_bounds_to_range_sa_so(rnLowerBoundSa, rnUpperBoundSa, \
98
                                 precisionSa=None):
99
    """
100
    Return a Sollya range from to 2 RealField Sage elements.
101
    The Sollya range element has a sufficient precision to hold all
102
    the digits of the widest of the Sage bounds.
103
    """
104
    # Sanity check.
105
    if rnLowerBoundSa > rnUpperBoundSa:
106
        return None
107
    # Precision stuff.
108
    if precisionSa is None:
109
        # Check for the largest precision.
110
        lbPrecSa = rnLowerBoundSa.parent().precision()
111
        ubPrecSa = rnLowerBoundSa.parent().precision()
112
        maxPrecSa = max(lbPrecSa, ubPrecSa)
113
    else:
114
        maxPrecSa = precisionSa
115
    # From Sage to Sollya bounds.
116
#    lowerBoundSo = sollya_lib_constant(get_rn_value(rnLowerBoundSa),
117
#                                       maxPrecSa)
118
    lowerBoundSo = pobyso_constant_sa_so(rnLowerBoundSa,
119
                                         maxPrecSa)
120
    upperBoundSo = pobyso_constant_sa_so(rnUpperBoundSa,
121
                                       maxPrecSa)
122
    
123
    # From Sollya bounds to range.
124
    rangeSo = sollya_lib_range(lowerBoundSo, upperBoundSo)
125
    # Back to original precision.
126
    # Clean up
127
    sollya_lib_clear_obj(lowerBoundSo)
128
    sollya_lib_clear_obj(upperBoundSo)
129
    return rangeSo
130
# End pobyso_bounds_to_range_sa_so
131

    
132
def pobyso_build_function_sub_so_so(exp1So, exp2So):
133
    return(sollya_lib_build_function_sub(exp1So, exp2So))
134

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

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

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

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

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

    
224
def pobyso_constant(rnArg):
225
    """ Legacy function. See pobyso_constant_sa_so. """
226
    return(pobyso_constant_sa_so(rnArg))
227
    
228
def pobyso_constant_sa_so(rnArgSa, precisionSa=None):
229
    """
230
    Create a Sollya constant from a Sage RealNumber.
231
    """
232
    # Precision stuff
233
    if precisionSa is None:
234
        precisionSa = rnArgSa.parent().precision()
235
    currentSollyaPrecisionSo = sollya_lib_get_prec()
236
    currentSollyaPrecisionSa = \
237
        pobyso_constant_from_int(currentSollyaPrecisionSo)
238
    # Sollya constant creation takes place here.
239
    if precisionSa > currentSollyaPrecisionSa:
240
        pobyso_set_prec_sa_so(precisionSa)
241
        constantSo = sollya_lib_constant(get_rn_value(rnArgSa))
242
        pobyso_set_prec_so_so(currentSollyaPrecisionSo)
243
    else:
244
        constantSo = sollya_lib_constant(get_rn_value(rnArgSa))
245
    sollya_lib_clear_obj(currentSollyaPrecisionSo)
246
    return constantSo
247
# End pobyso_constant_sa_so
248
     
249
def pobyso_constant_0_sa_so():
250
    """
251
    Obvious.
252
    """
253
    return(pobyso_constant_from_int_sa_so(0))
254

    
255
def pobyso_constant_1():
256
    """
257
    Obvious.
258
    Legacy function. See pobyso_constant_so_so. 
259
    """
260
    return(pobyso_constant_1_sa_so())
261

    
262
def pobyso_constant_1_sa_so():
263
    """
264
    Obvious.
265
    """
266
    return(pobyso_constant_from_int_sa_so(1))
267

    
268
def pobyso_constant_from_int(anInt):
269
    """ Legacy function. See pobyso_constant_from_int_sa_so. """
270
    return(pobyso_constant_from_int_sa_so(anInt))
271

    
272
def pobyso_constant_from_int_sa_so(anInt):
273
    """
274
    Get a Sollya constant from a Sage int.
275
    """
276
    return(sollya_lib_constant_from_int(int(anInt)))
277

    
278
def pobyso_constant_from_int_so_sa(constSo):
279
    """
280
    Get a Sage int from a Sollya int constant.
281
    Usefull for precision or powers in polynomials.
282
    """
283
    constSa = c_int(0)
284
    sollya_lib_get_constant_as_int(byref(constSa), constSo)
285
    return(constSa.value)
286
# End pobyso_constant_from_int_so_sa
287

    
288
def pobyso_error_so():
289
    return sollya_lib_error(None)
290
# End pobyso_error().
291

    
292
def pobyso_function_type_as_string(funcType):
293
    """ Legacy function. See pobyso_function_type_as_string_so_sa. """
294
    return(pobyso_function_type_as_string_so_sa(funcType))
295

    
296
def pobyso_function_type_as_string_so_sa(funcType):
297
    """
298
    Numeric Sollya function codes -> Sage mathematical function names.
299
    Notice that pow -> ^ (a la Sage, not a la Python).
300
    """
301
    if funcType == SOLLYA_BASE_FUNC_ABS:
302
        return "abs"
303
    elif funcType == SOLLYA_BASE_FUNC_ACOS:
304
        return "arccos"
305
    elif funcType == SOLLYA_BASE_FUNC_ACOSH:
306
        return "arccosh"
307
    elif funcType == SOLLYA_BASE_FUNC_ADD:
308
        return "+"
309
    elif funcType == SOLLYA_BASE_FUNC_ASIN:
310
        return "arcsin"
311
    elif funcType == SOLLYA_BASE_FUNC_ASINH:
312
        return "arcsinh"
313
    elif funcType == SOLLYA_BASE_FUNC_ATAN:
314
        return "arctan"
315
    elif funcType == SOLLYA_BASE_FUNC_ATANH:
316
        return "arctanh"
317
    elif funcType == SOLLYA_BASE_FUNC_CEIL:
318
        return "ceil"
319
    elif funcType == SOLLYA_BASE_FUNC_CONSTANT:
320
        return "cte"
321
    elif funcType == SOLLYA_BASE_FUNC_COS:
322
        return "cos"
323
    elif funcType == SOLLYA_BASE_FUNC_COSH:
324
        return "cosh"
325
    elif funcType == SOLLYA_BASE_FUNC_DIV:
326
        return "/"
327
    elif funcType == SOLLYA_BASE_FUNC_DOUBLE:
328
        return "double"
329
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEDOUBLE:
330
        return "doubleDouble"
331
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEEXTENDED:
332
        return "doubleDxtended"
333
    elif funcType == SOLLYA_BASE_FUNC_ERF:
334
        return "erf"
335
    elif funcType == SOLLYA_BASE_FUNC_ERFC:
336
        return "erfc"
337
    elif funcType == SOLLYA_BASE_FUNC_EXP:
338
        return "exp"
339
    elif funcType == SOLLYA_BASE_FUNC_EXP_M1:
340
        return "expm1"
341
    elif funcType == SOLLYA_BASE_FUNC_FLOOR:
342
        return "floor"
343
    elif funcType == SOLLYA_BASE_FUNC_FREE_VARIABLE:
344
        return "freeVariable"
345
    elif funcType == SOLLYA_BASE_FUNC_HALFPRECISION:
346
        return "halfPrecision"
347
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYCONSTANT:
348
        return "libraryConstant"
349
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYFUNCTION:
350
        return "libraryFunction"
351
    elif funcType == SOLLYA_BASE_FUNC_LOG:
352
        return "log"
353
    elif funcType == SOLLYA_BASE_FUNC_LOG_10:
354
        return "log10"
355
    elif funcType == SOLLYA_BASE_FUNC_LOG_1P:
356
        return "log1p"
357
    elif funcType == SOLLYA_BASE_FUNC_LOG_2:
358
        return "log2"
359
    elif funcType == SOLLYA_BASE_FUNC_MUL:
360
        return "*"
361
    elif funcType == SOLLYA_BASE_FUNC_NEARESTINT:
362
        return "round"
363
    elif funcType == SOLLYA_BASE_FUNC_NEG:
364
        return "__neg__"
365
    elif funcType == SOLLYA_BASE_FUNC_PI:
366
        return "pi"
367
    elif funcType == SOLLYA_BASE_FUNC_POW:
368
        return "^"
369
    elif funcType == SOLLYA_BASE_FUNC_PROCEDUREFUNCTION:
370
        return "procedureFunction"
371
    elif funcType == SOLLYA_BASE_FUNC_QUAD:
372
        return "quad"
373
    elif funcType == SOLLYA_BASE_FUNC_SIN:
374
        return "sin"
375
    elif funcType == SOLLYA_BASE_FUNC_SINGLE:
376
        return "single"
377
    elif funcType == SOLLYA_BASE_FUNC_SINH:
378
        return "sinh"
379
    elif funcType == SOLLYA_BASE_FUNC_SQRT:
380
        return "sqrt"
381
    elif funcType == SOLLYA_BASE_FUNC_SUB:
382
        return "-"
383
    elif funcType == SOLLYA_BASE_FUNC_TAN:
384
        return "tan"
385
    elif funcType == SOLLYA_BASE_FUNC_TANH:
386
        return "tanh"
387
    elif funcType == SOLLYA_BASE_FUNC_TRIPLEDOUBLE:
388
        return "tripleDouble"
389
    else:
390
        return None
391

    
392
def pobyso_get_constant(rnArgSa, constSo):
393
    """ Legacy function. See pobyso_get_constant_so_sa. """
394
    return(pobyso_get_constant_so_sa(rnArgSa, constSo))
395

    
396
def pobyso_get_constant_so_sa(rnArgSa, constSo):
397
    """
398
    Set the value of rnArgSo to the value of constSo in MPFR_RNDN mode.
399
    rnArg must already exist and belong to some RealField.
400
    We assume that constSo points to a Sollya constant.
401
    """
402
    return(sollya_lib_get_constant(get_rn_value(rnArgSa), constSo))
403
    
404
def pobyso_get_constant_as_rn(ctExpSo):
405
    """ 
406
    Legacy function. See pobyso_get_constant_as_rn_so_sa. 
407
    """ 
408
    return(pobyso_get_constant_as_rn_so_sa(ctExpSo))
409
    
410
def pobyso_get_constant_as_rn_so_sa(constExpSo):
411
    """
412
    Get a Sollya constant as a Sage "real number".
413
    The precision of the floating-point number returned is that of the Sollya
414
    constant.
415
    """
416
    precisionSa  = pobyso_get_prec_of_constant_so_sa(constExpSo) 
417
    RRRR = RealField(precisionSa)
418
    rnSa = RRRR(0)
419
    sollya_lib_get_constant(get_rn_value(rnSa), constExpSo)
420
    return(rnSa)
421
# End pobyso_get_constant_as_rn_so_sa
422

    
423
def pobyso_get_constant_as_rn_with_rf(ctExp, realField):
424
    """ 
425
    Legacy function. See pobyso_get_constant_as_rn_with_rf_so_sa.
426
    """
427
    return(pobyso_get_constant_as_rn_with_rf_so_sa(ctExp, realField))
428
    
429
def pobyso_get_constant_as_rn_with_rf_so_sa(ctExpSo, realFieldSa = None):
430
    """
431
    Get a Sollya constant as a Sage "real number".
432
    If no real field is specified, the precision of the floating-point number 
433
    returned is that of the Sollya constant.
434
    Otherwise is is that of the real field. Hence rounding may happen.
435
    """
436
    if realFieldSa is None:
437
        sollyaPrecSa = pobyso_get_prec_of_constant_so_sa(ctExpSo)
438
        realFieldSa = RealField(sollyaPrecSa)
439
    rnSa = realFieldSa(0)
440
    sollya_lib_get_constant(get_rn_value(rnSa), ctExpSo)
441
    return(rnSa)
442
# End pobyso_get_constant_as_rn_with_rf_so_sa
443

    
444
def pobyso_get_free_variable_name():
445
    """ 
446
    Legacy function. See pobyso_get_free_variable_name_so_sa.
447
    """
448
    return(pobyso_get_free_variable_name_so_sa())
449

    
450
def pobyso_get_free_variable_name_so_sa():
451
    return(sollya_lib_get_free_variable_name())
452
    
453
def pobyso_get_function_arity(expressionSo):
454
    """ 
455
    Legacy function. See pobyso_get_function_arity_so_sa.
456
    """
457
    return(pobyso_get_function_arity_so_sa(expressionSo))
458

    
459
def pobyso_get_function_arity_so_sa(expressionSo):
460
    arity = c_int(0)
461
    sollya_lib_get_function_arity(byref(arity),expressionSo)
462
    return(int(arity.value))
463

    
464
def pobyso_get_head_function(expressionSo):
465
    """ 
466
    Legacy function. See pobyso_get_head_function_so_sa. 
467
    """
468
    return(pobyso_get_head_function_so_sa(expressionSo)) 
469

    
470
def pobyso_get_head_function_so_sa(expressionSo):
471
    functionType = c_int(0)
472
    sollya_lib_get_head_function(byref(functionType), expressionSo, None)
473
    return(int(functionType.value))
474

    
475
def pobyso_get_interval_from_range_so_sa(soRange, realIntervalFieldSa = None ):
476
    """
477
    Return the Sage interval corresponding to the Sollya range argument.
478
    If no reaIntervalField is passed as an argument, the interval bounds are not
479
    rounded: they are elements of RealIntervalField of the "right" precision
480
    to hold all the digits.
481
    """
482
    prec = c_int(0)
483
    if realIntervalFieldSa is None:
484
        retval = sollya_lib_get_prec_of_range(byref(prec), soRange, None)
485
        if retval == 0:
486
            return(None)
487
        realIntervalFieldSa = RealIntervalField(prec.value)
488
    intervalSa = realIntervalFieldSa(0,0)
489
    retval = \
490
        sollya_lib_get_interval_from_range(get_interval_value(intervalSa),\
491
                                           soRange)
492
    if retval == 0:
493
        return(None)
494
    return(intervalSa)
495
# End pobyso_get_interval_from_range_so_sa
496

    
497
def pobyso_get_list_elements(soObj):
498
    """ Legacy function. See pobyso_get_list_elements_so_so. """
499
    return(pobyso_get_list_elements_so_so(soObj))
500
 
501
def pobyso_get_list_elements_so_so(objectListSo):
502
    """
503
    Get the Sollya list elements as a Sage/Python array of Sollya objects.
504
    
505
    INPUT:
506
    - objectListSo: a Sollya list of Sollya objects.
507
    
508
    OUTPUT:
509
    - a Sage/Python tuple made of:
510
      - a Sage/Python list of Sollya objects,
511
      - a Sage/Python int holding the number of elements,
512
      - a Sage/Python int stating (!= 0) that the list is end-elliptic.
513
    NOTE::
514
        We recover the addresses of the Sollya object from the list of pointers
515
        returned by sollya_lib_get_list_elements. The list itself is freed.
516
    TODO::
517
        Figure out what to do with numElements since the number of elements
518
        can easily be recovered from the list itself. 
519
        Ditto for isEndElliptic.
520
    """
521
    listAddress = POINTER(c_longlong)()
522
    numElements = c_int(0)
523
    isEndElliptic = c_int(0)
524
    listAsSageList = []
525
    result = sollya_lib_get_list_elements(byref(listAddress),\
526
                                          byref(numElements),\
527
                                          byref(isEndElliptic),\
528
                                          objectListSo)
529
    if result == 0 :
530
        return None
531
    for i in xrange(0, numElements.value, 1):
532
       #listAsSageList.append(sollya_lib_copy_obj(listAddress[i]))
533
       listAsSageList.append(listAddress[i])
534
       # Clear each of the elements returned by Sollya.
535
       #sollya_lib_clear_obj(listAddress[i])
536
    # Free the list itself.   
537
    sollya_lib_free(listAddress)
538
    return(listAsSageList, numElements.value, isEndElliptic.value)
539

    
540
def pobyso_get_max_prec_of_exp(soExp):
541
    """ Legacy function. See pobyso_get_max_prec_of_exp_so_sa. """
542
    return(pobyso_get_max_prec_of_exp_so_sa(soExp))
543

    
544
def pobyso_get_max_prec_of_exp_so_sa(expSo):
545
    """
546
    Get the maximum precision used for the numbers in a Sollya expression.
547
    
548
    Arguments:
549
    soExp -- a Sollya expression pointer
550
    Return value:
551
    A Python integer
552
    TODO: 
553
    - error management;
554
    - correctly deal with numerical type such as DOUBLEEXTENDED.
555
    """
556
    maxPrecision = 0
557
    minConstPrec = 0
558
    currentConstPrec = 0
559
    operator = pobyso_get_head_function_so_sa(expSo)
560
    if (operator != SOLLYA_BASE_FUNC_CONSTANT) and \
561
    (operator != SOLLYA_BASE_FUNC_FREE_VARIABLE):
562
        (arity, subexpressions) = pobyso_get_subfunctions_so_sa(expSo)
563
        for i in xrange(arity):
564
            maxPrecisionCandidate = \
565
                pobyso_get_max_prec_of_exp_so_sa(subexpressions[i])
566
            if maxPrecisionCandidate > maxPrecision:
567
                maxPrecision = maxPrecisionCandidate
568
        return(maxPrecision)
569
    elif operator == SOLLYA_BASE_FUNC_CONSTANT:
570
        #minConstPrec = pobyso_get_min_prec_of_constant_so_sa(expSo)
571
        #currentConstPrec = pobyso_get_min_prec_of_constant_so_sa(soExp)
572
        #print minConstPrec, " - ", currentConstPrec 
573
        return(pobyso_get_min_prec_of_constant_so_sa(expSo))
574
    
575
    elif operator == SOLLYA_BASE_FUNC_FREE_VARIABLE:
576
        return(0)
577
    else:
578
        print "pobyso_get_max_prec_of_exp_so_sa: unexepected operator."
579
        return(0)
580

    
581
def pobyso_get_min_prec_of_constant_so_sa(constExpSo):
582
    """
583
    Get the minimum precision necessary to represent the value of a Sollya
584
    constant.
585
    MPFR_MIN_PREC and powers of 2 are taken into account.
586
    We assume that constExpSo is a point
587
    """
588
    constExpAsRnSa = pobyso_get_constant_as_rn_so_sa(constExpSo)
589
    return(min_mpfr_size(get_rn_value(constExpAsRnSa)))
590

    
591
def pobyso_get_sage_exp_from_sollya_exp(sollyaExpSo, realField = RR):
592
    """ Legacy function. See pobyso_get_sage_exp_from_sollya_exp_so_sa. """
593
    return(pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, \
594
                                                     realField = RR))
595

    
596
def pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, realFieldSa = RR):
597
    """
598
    Get a Sage expression from a Sollya expression. 
599
    Currently only tested with polynomials with floating-point coefficients.
600
    Notice that, in the returned polynomial, the exponents are RealNumbers.
601
    """
602
    #pobyso_autoprint(sollyaExp)
603
    operatorSa = pobyso_get_head_function_so_sa(sollyaExpSo)
604
    sollyaLibFreeVariableName = sollya_lib_get_free_variable_name()
605
    # Constants and the free variable are special cases.
606
    # All other operator are dealt with in the same way.
607
    if (operatorSa != SOLLYA_BASE_FUNC_CONSTANT) and \
608
       (operatorSa != SOLLYA_BASE_FUNC_FREE_VARIABLE):
609
        (aritySa, subexpressionsSa) = pobyso_get_subfunctions_so_sa(sollyaExpSo)
610
        if aritySa == 1:
611
            sageExpSa = eval(pobyso_function_type_as_string_so_sa(operatorSa) + \
612
            "(" + pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], \
613
            realFieldSa) + ")")
614
        elif aritySa == 2:
615
            # We do not get through the preprocessor.
616
            # The "^" operator is then a special case.
617
            if operatorSa == SOLLYA_BASE_FUNC_POW:
618
                operatorAsStringSa = "**"
619
            else:
620
                operatorAsStringSa = \
621
                    pobyso_function_type_as_string_so_sa(operatorSa)
622
            sageExpSa = \
623
              eval("pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], realFieldSa)"\
624
              + " " + operatorAsStringSa + " " + \
625
                   "pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[1], realFieldSa)")
626
        # We do not know yet how to deal with arity >= 3 
627
        # (is there any in Sollya anyway?).
628
        else:
629
            sageExpSa = eval('None')
630
        return(sageExpSa)
631
    elif operatorSa == SOLLYA_BASE_FUNC_CONSTANT:
632
        #print "This is a constant"
633
        return pobyso_get_constant_as_rn_with_rf_so_sa(sollyaExpSo, realFieldSa)
634
    elif operatorSa == SOLLYA_BASE_FUNC_FREE_VARIABLE:
635
        #print "This is free variable"
636
        return(eval(sollyaLibFreeVariableName))
637
    else:
638
        print "Unexpected"
639
        return eval('None')
640
# End pobyso_get_sage_exp_from_sollya_exp_so_sa
641

    
642

    
643
def pobyso_get_poly_sa_so(polySo, realFieldSa=None):
644
    """
645
    Create a Sollya polynomial from a Sage polynomial.
646
    """
647
    pass
648
# pobyso_get_poly_sa_so
649

    
650
def pobyso_get_poly_so_sa(polySo, realFieldSa=None):
651
    """
652
    Convert a Sollya polynomial into a Sage polynomial.
653
    We assume that the polynomial is in canonical form.
654
    If no realField is given, a RealField corresponding to the maximum 
655
    precision of the coefficients is internally computed.
656
    The real field is not returned but can be easily retrieved from 
657
    the polynomial itself.
658
    ALGORITHM:
659
    - (optional) compute the RealField of the coefficients;
660
    - convert the Sollya expression into a Sage expression;
661
    - convert the Sage expression into a Sage polynomial
662
    TODO: the canonical thing for the polynomial.
663
    """    
664
    if realFieldSa is None:
665
        expressionPrecSa = pobyso_get_max_prec_of_exp_so_sa(polySo)
666
        realFieldSa = RealField(expressionPrecSa)
667
    #print "Sollya expression before...",
668
    #pobyso_autoprint(polySo)
669

    
670
    expressionSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, \
671
                                                             realFieldSa)
672
    #print "...Sollya expression after.",
673
    #pobyso_autoprint(polySo)
674
    polyVariableSa = expressionSa.variables()[0]
675
    polyRingSa = realFieldSa[str(polyVariableSa)]
676
    #print polyRingSa
677
    # Do not use the polynomial(expressionSa, ring=polyRingSa) form!
678
    polynomialSa = polyRingSa(expressionSa)
679
    return(polynomialSa)
680
# End pobyso_get_sage_poly_from_sollya_poly
681

    
682
def pobyso_get_subfunctions(expressionSo):
683
    """ Legacy function. See pobyso_get_subfunctions_so_sa. """
684
    return(pobyso_get_subfunctions_so_sa(expressionSo)) 
685

    
686
def pobyso_get_subfunctions_so_sa(expressionSo):
687
    """
688
    Get the subfunctions of an expression.
689
    Return the number of subfunctions and the list of subfunctions addresses.
690
    S.T.: Could not figure out another way than that ugly list of declarations
691
    to recover the addresses of the subfunctions.
692
    We limit ourselves to arity 8 functions. 
693
    """
694
    subf0 = c_int(0)
695
    subf1 = c_int(0)
696
    subf2 = c_int(0)
697
    subf3 = c_int(0)
698
    subf4 = c_int(0)
699
    subf5 = c_int(0)
700
    subf6 = c_int(0)
701
    subf7 = c_int(0)
702
    subf8 = c_int(0)
703
    arity = c_int(0)
704
    nullPtr = POINTER(c_int)()
705
    sollya_lib_get_subfunctions(expressionSo, byref(arity), \
706
      byref(subf0), byref(subf1), byref(subf2), byref(subf3), \
707
      byref(subf4), byref(subf5),\
708
      byref(subf6), byref(subf7), byref(subf8), nullPtr, None) 
709
#    byref(cast(subfunctions[0], POINTER(c_int))), \
710
#    byref(cast(subfunctions[0], POINTER(c_int))), \
711
#    byref(cast(subfunctions[2], POINTER(c_int))), \
712
#    byref(cast(subfunctions[3], POINTER(c_int))), \
713
#    byref(cast(subfunctions[4], POINTER(c_int))), \
714
#    byref(cast(subfunctions[5], POINTER(c_int))), \
715
#    byref(cast(subfunctions[6], POINTER(c_int))), \
716
#    byref(cast(subfunctions[7], POINTER(c_int))), \
717
#    byref(cast(subfunctions[8], POINTER(c_int))), nullPtr)
718
    subfunctions = [subf0, subf1, subf2, subf3, subf4, subf5, subf6, subf7, \
719
                    subf8]
720
    subs = []
721
    if arity.value > pobyso_max_arity:
722
        return(0,[])
723
    for i in xrange(arity.value):
724
        subs.append(int(subfunctions[i].value))
725
        #print subs[i]
726
    return(int(arity.value), subs)
727
    
728
def pobyso_get_prec():
729
    """ Legacy function. See pobyso_get_prec_so_sa(). """
730
    return(pobyso_get_prec_so_sa())
731

    
732
def pobyso_get_prec_so():
733
    """
734
    Get the current default precision in Sollya.
735
    The return value is a Sollya object.
736
    Usefull when modifying the precision back and forth by avoiding
737
    extra conversions.
738
    """
739
    return(sollya_lib_get_prec(None))
740
    
741
def pobyso_get_prec_so_sa():
742
    """
743
    Get the current default precision in Sollya.
744
    The return value is Sage/Python int.
745
    """
746
    precSo = sollya_lib_get_prec(None)
747
    precSa = c_int(0)
748
    sollya_lib_get_constant_as_int(byref(precSa), precSo)
749
    sollya_lib_clear_obj(precSo)
750
    return int(precSa.value)
751
# End pobyso_get_prec_so_sa.
752

    
753
def pobyso_get_prec_of_constant(ctExpSo):
754
    """ Legacy function. See pobyso_get_prec_of_constant_so_sa. """
755
    return(pobyso_get_prec_of_constant_so_sa(ctExpSo))
756

    
757
def pobyso_get_prec_of_constant_so_sa(ctExpSo):
758
    prec = c_int(0)
759
    retc = sollya_lib_get_prec_of_constant(byref(prec), ctExpSo, None)
760
    if retc == 0:
761
        return(None)
762
    return(int(prec.value))
763

    
764
def pobyso_get_prec_of_range_so_sa(rangeSo):
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_guess_degree_sa_sa(functionSa, intervalSa, approxErrorSa, 
773
                              weightSa=None, degreeBoundSa=None):
774
    """
775
    Sa_sa variant of the solly_guessdegree function.
776
    Return 0 if something goes wrong.
777
    """
778
    functionAsStringSa = functionSa._assume_str().replace('_SAGE_VAR_', '')
779
    functionSo = pobyso_parse_string_sa_so(functionAsStringSa)
780
    if pobyso_is_error_so_sa(functionSo):
781
        sollya_lib_clear_obj(functionSo)
782
        return 0
783
    rangeSo = pobyso_interval_to_range_sa_so(intervalSa)
784
    # The approximation error is expected to be a floating point number.
785
    if pobyso_is_floating_point_number_sa_sa(approxErrorSa):
786
        approxErrorSo = pobyso_constant_sa_so(approxErrorSa)
787
    else:
788
        approxErrorSo = pobyso_constant_sa_so(RR(approxErrorSa))
789
    if not weightSa is None:
790
        weightAsStringSa = weightSa._assume_str().replace('_SAGE_VAR_', '')
791
        weightSo = pobyso_parse_string_sa_so(weightAsStringSa)
792
        if pobyso_is_error_so_sa(weightSo):
793
            sollya_lib_clear_obj(functionSo)
794
            sollya_lib_clear_obj(rangeSo)
795
            sollya_lib_clear_obj(approxErrorSo)   
796
            sollya_lib_clear_obj(weightSo)
797
            return 0   
798
    else:
799
        weightSo = None
800
    if not degreeBoundSa is None:
801
        degreeBoundSo = pobyso_constant_from_int_sa_so(degreeBoundSa)
802
    else:
803
        degreeBoundSo = None
804
    guessedDegreeSa = pobyso_guess_degree_so_sa(functionSo,
805
                                                rangeSo,
806
                                                approxErrorSo,
807
                                                weightSo,
808
                                                degreeBoundSo)
809
    sollya_lib_clear_obj(functionSo)
810
    sollya_lib_clear_obj(rangeSo)
811
    sollya_lib_clear_obj(approxErrorSo)
812
    if not weightSo is None:
813
        sollya_lib_clear_obj(weightSo)
814
    if not degreeBoundSo is None:
815
        sollya_lib_clear_obj(degreeBoundSo)
816
    return guessedDegreeSa
817
# End poyso_guess_degree_sa_sa
818

    
819
def pobyso_guess_degree_so_sa(functionSo, rangeSo, errorSo, weightSo=None, \
820
                              degreeBoundSo=None):
821
    """
822
    Thin wrapper around the guessdegree function.
823
    Nevertheless, some precision control stuff has been appended.
824
    """
825
    # Deal with Sollya internal precision issues: if it is too small, 
826
    # compared with the error, increases it to about twice -log2(error).
827
    errorSa = pobyso_get_constant_as_rn_with_rf_so_sa(errorSo)
828
    log2ErrorSa = errorSa.log2()
829
    if log2ErrorSa < 0:
830
        neededPrecisionSa = int(2 * int(-log2ErrorSa) / 64) * 64
831
    else:
832
        neededPrecisionSa = int(2 * int(log2ErrorSa) / 64) * 64
833
    #print "Needed precision:", neededPrecisionSa
834
    currentPrecSa = pobyso_get_prec_so_sa()
835
    if neededPrecisionSa > currentPrecSa:
836
        currentPrecSo = pobyso_get_prec_so()
837
        pobyso_set_prec_sa_so(neededPrecisionSa)
838
    #print "Guessing degree..."
839
    # weightSo and degreeBoundsSo are optional arguments.
840
    # As declared, sollya_lib_guessdegree must take 5 arguments.
841
    if weightSo is None:
842
        degreeRangeSo = sollya_lib_guessdegree(functionSo, rangeSo, errorSo, 
843
                                               0, 0, None)
844
    elif degreeBoundSo is None:
845
        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, \
846
                                                errorSo, weightSo, 0, None)
847
    else:
848
        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, errorSo, \
849
                                                weightSo, degreeBoundSo, None)
850
    #print "...degree guess done."
851
    # Restore internal precision, if applicable.
852
    if neededPrecisionSa > currentPrecSa:
853
        pobyso_set_prec_so_so(currentPrecSo)
854
        sollya_lib_clear_obj(currentPrecSo)
855
    degreeIntervalSa = pobyso_range_to_interval_so_sa(degreeRangeSo)
856
    sollya_lib_clear_obj(degreeRangeSo)
857
    # When ok, both bounds match.
858
    # When the degree bound is too low, the upper bound is the degree
859
    # for which the error can be honored.
860
    # When it really goes wrong, the upper bound is infinity. 
861
    if degreeIntervalSa.lower() == degreeIntervalSa.upper():
862
        return int(degreeIntervalSa.lower())
863
    else:
864
        if degreeIntervalSa.upper().is_infinity():
865
            return None
866
        else:
867
            return int(degreeIntervalSa.upper())
868
    # End pobyso_guess_degree_so_sa
869

    
870
def pobyso_infnorm_so_so(func, interval, file = None, intervalList = None):
871
    print "Do not use this function. User pobyso_supnorm_so_so instead."
872
    return(None)
873

    
874
def pobyso_interval_to_range_sa_so(intervalSa, precisionSa=None):
875
    if precisionSa is None:
876
        precisionSa = intervalSa.parent().precision()
877
    intervalSo = pobyso_bounds_to_range_sa_so(intervalSa.lower(),\
878
                                              intervalSa.upper(),\
879
                                              precisionSa)
880
    return(intervalSo)
881
# End pobyso_interval_to_range_sa_so
882

    
883
def pobyso_is_error_so_sa(objSo):
884
    """
885
    Thin wrapper around the sollya_lib_obj_is_error() function.
886
    """
887
    if sollya_lib_obj_is_error(objSo) != 0:
888
        return True
889
    else:
890
        return False
891
# End pobyso_is_error-so_sa
892

    
893
def pobyso_is_floating_point_number_sa_sa(numberSa):
894
    """
895
    Check whether a Sage number is floating point
896
    """
897
    return numberSa.parent().__class__ is RR.__class__
898

    
899
def pobyso_lib_init():
900
    sollya_lib_init(None)
901

    
902
def pobyso_lib_close():
903
    sollya_lib_close(None)
904
    
905
def pobyso_name_free_variable(freeVariableNameSa):
906
    """ Legacy function. See pobyso_name_free_variable_sa_so. """
907
    pobyso_name_free_variable_sa_so(freeVariableNameSa)
908

    
909
def pobyso_name_free_variable_sa_so(freeVariableNameSa):
910
    """
911
    Set the free variable name in Sollya from a Sage string.
912
    """
913
    sollya_lib_name_free_variable(freeVariableNameSa)
914

    
915
def pobyso_parse_string(string):
916
    """ Legacy function. See pobyso_parse_string_sa_so. """
917
    return(pobyso_parse_string_sa_so(string))
918
 
919
def pobyso_parse_string_sa_so(string):
920
    """
921
    Get the Sollya expression computed from a Sage string or
922
    a Sollya error object if parsing failed.
923
    """
924
    return(sollya_lib_parse_string(string))
925

    
926
def pobyso_range(rnLowerBound, rnUpperBound):
927
    """ Legacy function. See pobyso_range_sa_so. """
928
    return(pobyso_range_sa_so(rnLowerBound, rnUpperBound)) 
929

    
930

    
931
def pobyso_range_to_interval_so_sa(rangeSo, realIntervalFieldSa = None):
932
    """
933
    Get a Sage interval from a Sollya range.
934
    If no realIntervalField is given as a parameter, the Sage interval
935
    precision is that of the Sollya range.
936
    Otherwise, the precision is that of the realIntervalField. In this case
937
    rounding may happen.
938
    """
939
    if realIntervalFieldSa is None:
940
        precSa = pobyso_get_prec_of_range_so_sa(rangeSo)
941
        realIntervalFieldSa = RealIntervalField(precSa)
942
    intervalSa = \
943
        pobyso_get_interval_from_range_so_sa(rangeSo, realIntervalFieldSa) 
944
    return(intervalSa)
945

    
946
def pobyso_remez_canonical_sa_sa(func, \
947
                                 degree, \
948
                                 lowerBound, \
949
                                 upperBound, \
950
                                 weight = None, \
951
                                 quality = None):
952
    """
953
    All arguments are Sage/Python.
954
    The functions (func and weight) must be passed as expressions or strings.
955
    Otherwise the function fails. 
956
    The return value is a Sage polynomial.
957
    """
958
    var('zorglub')    # Dummy variable name for type check only. Type of 
959
    # zorglub is "symbolic expression".
960
    polySo = pobyso_remez_canonical_sa_so(func, \
961
                                 degree, \
962
                                 lowerBound, \
963
                                 upperBound, \
964
                                 weight, \
965
                                 quality)
966
    # String test
967
    if parent(func) == parent("string"):
968
        functionSa = eval(func)
969
    # Expression test.
970
    elif type(func) == type(zorglub):
971
        functionSa = func
972
    else:
973
        return None
974
    #
975
    maxPrecision = 0
976
    if polySo is None:
977
        return(None)
978
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
979
    RRRRSa = RealField(maxPrecision)
980
    polynomialRingSa = RRRRSa[functionSa.variables()[0]]
981
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, RRRRSa)
982
    polySa = polynomial(expSa, polynomialRingSa)
983
    sollya_lib_clear_obj(polySo)
984
    return(polySa)
985
# End pobyso_remez_canonical_sa_sa
986
    
987
def pobyso_remez_canonical(func, \
988
                           degree, \
989
                           lowerBound, \
990
                           upperBound, \
991
                           weight = "1", \
992
                           quality = None):
993
    """ Legacy function. See pobyso_remez_canonical_sa_so. """
994
    return(pobyso_remez_canonical_sa_so(func, \
995
                                        degree, \
996
                                        lowerBound, \
997
                                        upperBound, \
998
                                        weight, \
999
                                        quality))
1000
def pobyso_remez_canonical_sa_so(func, \
1001
                                 degree, \
1002
                                 lowerBound, \
1003
                                 upperBound, \
1004
                                 weight = None, \
1005
                                 quality = None):
1006
    """
1007
    All arguments are Sage/Python.
1008
    The functions (func and weight) must be passed as expressions or strings.
1009
    Otherwise the function fails. 
1010
    The return value is a pointer to a Sollya function.
1011
    """
1012
    var('zorglub')    # Dummy variable name for type check only. Type of
1013
    # zorglub is "symbolic expression".
1014
    currentVariableNameSa = None
1015
    # The func argument can be of different types (string, 
1016
    # symbolic expression...)
1017
    if parent(func) == parent("string"):
1018
        localFuncSa = eval(func)
1019
        if len(localFuncSa.variables()) > 0:
1020
            currentVariableNameSa = localFuncSa.variables()[0]
1021
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1022
            functionSo = \
1023
              sollya_lib_parse_string(localFuncSa._assume_str().replace('_SAGE_VAR_', ''))
1024
    # Expression test.
1025
    elif type(func) == type(zorglub):
1026
        # Until we are able to translate Sage expressions into Sollya 
1027
        # expressions : parse the string version.
1028
        if len(func.variables()) > 0:
1029
            currentVariableNameSa = func.variables()[0]
1030
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1031
            functionSo = \
1032
              sollya_lib_parse_string(func._assume_str().replace('_SAGE_VAR_', ''))
1033
    else:
1034
        return(None)
1035
    if weight is None: # No weight given -> 1.
1036
        weightSo = pobyso_constant_1_sa_so()
1037
    elif parent(weight) == parent("string"): # Weight given as string: parse it.
1038
        weightSo = sollya_lib_parse_string(func)
1039
    elif type(weight) == type(zorglub): # Weight given as symbolice expression.
1040
        functionSo = \
1041
          sollya_lib_parse_string_sa_so(weight._assume_str().replace('_SAGE_VAR_', ''))
1042
    else:
1043
        return(None)
1044
    degreeSo = pobyso_constant_from_int(degree)
1045
    rangeSo = pobyso_bounds_to_range_sa_so(lowerBound, upperBound)
1046
    if not quality is None:
1047
        qualitySo= pobyso_constant_sa_so(quality)
1048
    else:
1049
        qualitySo = None
1050
        
1051
    remezPolySo = sollya_lib_remez(functionSo, \
1052
                                   degreeSo, \
1053
                                   rangeSo, \
1054
                                   weightSo, \
1055
                                   qualitySo, \
1056
                                   None)
1057
    sollya_lib_clear_obj(functionSo)
1058
    sollya_lib_clear_obj(degreeSo)
1059
    sollya_lib_clear_obj(rangeSo)
1060
    sollya_lib_clear_obj(weightSo)
1061
    if not qualitySo is None:
1062
        sollya_lib_clear_obj(qualitySo)
1063
    return(remezPolySo)
1064
# End pobyso_remez_canonical_sa_so
1065

    
1066
def pobyso_remez_canonical_so_so(funcSo, \
1067
                                 degreeSo, \
1068
                                 rangeSo, \
1069
                                 weightSo = pobyso_constant_1_sa_so(),\
1070
                                 qualitySo = None):
1071
    """
1072
    All arguments are pointers to Sollya objects.
1073
    The return value is a pointer to a Sollya function.
1074
    """
1075
    if not sollya_lib_obj_is_function(funcSo):
1076
        return(None)
1077
    return(sollya_lib_remez(funcSo, degreeSo, rangeSo, weightSo, qualitySo, None))
1078
    
1079
def pobyso_set_canonical_off():
1080
    sollya_lib_set_canonical(sollya_lib_off())
1081

    
1082
def pobyso_set_canonical_on():
1083
    sollya_lib_set_canonical(sollya_lib_on())
1084

    
1085
def pobyso_set_prec(p):
1086
    """ Legacy function. See pobyso_set_prec_sa_so. """
1087
    pobyso_set_prec_sa_so(p)
1088

    
1089
def pobyso_set_prec_sa_so(p):
1090
    a = c_int(p)
1091
    precSo = c_void_p(sollya_lib_constant_from_int(a))
1092
    sollya_lib_set_prec(precSo, None)
1093

    
1094
def pobyso_set_prec_so_so(newPrecSo):
1095
    sollya_lib_set_prec(newPrecSo, None)
1096

    
1097
def pobyso_supnorm_so_so(polySo, funcSo, intervalSo, errorTypeSo = None,\
1098
                         accuracySo = None):
1099
    """
1100
    Computes the supnorm of the approximation error between the given 
1101
    polynomial and function.
1102
    errorTypeSo defaults to "absolute".
1103
    accuracySo defaults to 2^(-40).
1104
    """
1105
    if errorTypeSo is None:
1106
        errorTypeSo = sollya_lib_absolute(None)
1107
        errorTypeIsNone = True
1108
    else:
1109
        errorTypeIsNone = False
1110
    #
1111
    if accuracySo is None:
1112
        # Notice the **!
1113
        accuracySo = pobyso_constant_sa_so(RR(2**(-40)))
1114
        accuracyIsNone = True
1115
    else:
1116
        accuracyIsNone = False
1117
    pobyso_autoprint(accuracySo)
1118
    resultSo = \
1119
        sollya_lib_supnorm(polySo, funcSo, intervalSo, errorTypeSo, \
1120
                              accuracySo)
1121
    if errorTypeIsNone:
1122
        sollya_lib_clear_obj(errorTypeSo)
1123
    if accuracyIsNone:
1124
        sollya_lib_clear_obj(accuracySo)
1125
    return resultSo
1126
# End pobyso_supnorm_so_so
1127

    
1128
def pobyso_taylor_expansion_no_change_var_so_so(functionSo, 
1129
                                                degreeSo, 
1130
                                                rangeSo,
1131
                                                errorTypeSo=None, 
1132
                                                sollyaPrecSo=None):
1133
    """
1134
    Compute the Taylor expansion without the variable change
1135
    x -> x-intervalCenter.
1136
    """
1137
    # No global change of the working precision.
1138
    if not sollyaPrecSo is None:
1139
        initialPrecSo = sollya_lib_get_prec(None)
1140
        sollya_lib_set_prec(sollyaPrecSo)
1141
    # Error type stuff: default to absolute.
1142
    if errorTypeSo is None:
1143
        errorTypeIsNone = True
1144
        errorTypeSo = sollya_lib_absolute(None)
1145
    else:
1146
        errorTypeIsNone = False
1147
    intervalCenterSo = sollya_lib_mid(rangeSo, None)
1148
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo,
1149
                                         intervalCenterSo,
1150
                                         rangeSo, errorTypeSo, None)
1151
    # taylorFormListSaSo is a Python list of Sollya objects references that 
1152
    # are copies of the elements of taylorFormSo.
1153
    # pobyso_get_list_elements_so_so clears taylorFormSo.
1154
    (taylorFormListSaSo, numElementsSa, isEndEllipticSa) = \
1155
        pobyso_get_list_elements_so_so(taylorFormSo)
1156
    polySo = sollya_lib_copy_obj(taylorFormListSaSo[0])
1157
    #print "Num elements:", numElementsSa
1158
    sollya_lib_clear_obj(taylorFormSo)
1159
    #polySo = taylorFormListSaSo[0]
1160
    #errorRangeSo = sollya_lib_copy_obj(taylorFormListSaSo[2])
1161
    errorRangeSo = taylorFormListSaSo[2]
1162
    # No copy_obj needed here: a new objects are created.
1163
    maxErrorSo    = sollya_lib_sup(errorRangeSo)
1164
    minErrorSo    = sollya_lib_inf(errorRangeSo)
1165
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
1166
    absMinErrorSo = sollya_lib_abs(minErrorSo)
1167
    sollya_lib_clear_obj(maxErrorSo)
1168
    sollya_lib_clear_obj(minErrorSo)
1169
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
1170
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
1171
    # If changed, reset the Sollya working precision.
1172
    if not sollyaPrecSo is None:
1173
        sollya_lib_set_prec(initialPrecSo)
1174
        sollya_lib_clear_obj(initialPrecSo)
1175
    if errorTypeIsNone:
1176
        sollya_lib_clear_obj(errorTypeSo)
1177
    pobyso_clear_taylorform_sa_so(taylorFormListSaSo)
1178
    if absMaxErrorSa > absMinErrorSa:
1179
        sollya_lib_clear_obj(absMinErrorSo)
1180
        return((polySo, intervalCenterSo, absMaxErrorSo))
1181
    else:
1182
        sollya_lib_clear_obj(absMaxErrorSo)
1183
        return((polySo, intervalCenterSo, absMinErrorSo))
1184
# end pobyso_taylor_expansion_no_change_var_so_so
1185

    
1186
def pobyso_taylor_expansion_with_change_var_so_so(functionSo, degreeSo, \
1187
                                                  rangeSo, \
1188
                                                  errorTypeSo=None, \
1189
                                                  sollyaPrecSo=None):
1190
    """
1191
    Compute the Taylor expansion with the variable change
1192
    x -> (x-intervalCenter) included.
1193
    """
1194
    # No global change of the working precision.
1195
    if not sollyaPrecSo is None:
1196
        initialPrecSo = sollya_lib_get_prec(None)
1197
        sollya_lib_set_prec(sollyaPrecSo)
1198
    #
1199
    # Error type stuff: default to absolute.
1200
    if errorTypeSo is None:
1201
        errorTypeIsNone = True
1202
        errorTypeSo = sollya_lib_absolute(None)
1203
    else:
1204
        errorTypeIsNone = False
1205
    intervalCenterSo = sollya_lib_mid(rangeSo)
1206
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo, \
1207
                                         intervalCenterSo, \
1208
                                         rangeSo, errorTypeSo, None)
1209
    # taylorFormListSaSo is a Python list of Sollya objects references that 
1210
    # are copies of the elements of taylorFormSo.
1211
    # pobyso_get_list_elements_so_so clears taylorFormSo.
1212
    (taylorFormListSo, numElements, isEndElliptic) = \
1213
        pobyso_get_list_elements_so_so(taylorFormSo)
1214
    polySo = taylorFormListSo[0]
1215
    errorRangeSo = taylorFormListSo[2]
1216
    maxErrorSo    = sollya_lib_sup(errorRangeSo)
1217
    minErrorSo    = sollya_lib_inf(errorRangeSo)
1218
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
1219
    absMinErrorSo = sollya_lib_abs(minErrorSo)
1220
    sollya_lib_clear_obj(maxErrorSo)
1221
    sollya_lib_clear_obj(minErrorSo)
1222
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
1223
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
1224
    changeVarExpSo = sollya_lib_build_function_sub(\
1225
                       sollya_lib_build_function_free_variable(),\
1226
                       sollya_lib_copy_obj(intervalCenterSo))
1227
    polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpSo)
1228
    sollya_lib_clear_obj(polySo) 
1229
    sollya_lib_clear_obj(changeVarExpSo)
1230
    # If changed, reset the Sollya working precision.
1231
    if not sollyaPrecSo is None:
1232
        sollya_lib_set_prec(initialPrecSo)
1233
        sollya_lib_clear_obj(initialPrecSo)
1234
    if errorTypeIsNone:
1235
        sollya_lib_clear_obj(errorTypeSo)
1236
    sollya_lib_clear_obj(taylorFormSo)
1237
    # Do not clear maxErrorSo.
1238
    if absMaxErrorSa > absMinErrorSa:
1239
        sollya_lib_clear_obj(absMinErrorSo)
1240
        return((polyVarChangedSo, intervalCenterSo, absMaxErrorSo))
1241
    else:
1242
        sollya_lib_clear_obj(absMaxErrorSo)
1243
        return((polyVarChangedSo, intervalCenterSo, absMinErrorSo))
1244
# end pobyso_taylor_expansion_with_change_var_so_so
1245

    
1246
def pobyso_taylor(function, degree, point):
1247
    """ Legacy function. See pobysoTaylor_so_so. """
1248
    return(pobyso_taylor_so_so(function, degree, point))
1249

    
1250
def pobyso_taylor_so_so(functionSo, degreeSo, pointSo):
1251
    return(sollya_lib_taylor(functionSo, degreeSo, pointSo))
1252
    
1253
def pobyso_taylorform(function, degree, point = None, 
1254
                      interval = None, errorType=None):
1255
    """ Legacy function. See pobyso_taylorform_sa_sa;"""
1256
    
1257
def pobyso_taylorform_sa_sa(functionSa, \
1258
                            degreeSa, \
1259
                            pointSa, \
1260
                            intervalSa=None, \
1261
                            errorTypeSa=None, \
1262
                            precisionSa=None):
1263
    """
1264
    Compute the Taylor form of 'degreeSa' for 'functionSa' at 'pointSa' 
1265
    for 'intervalSa' with 'errorTypeSa' (a string) using 'precisionSa'. 
1266
    point: must be a Real or a Real interval.
1267
    return the Taylor form as an array
1268
    TODO: take care of the interval and of the point when it is an interval;
1269
          when errorType is not None;
1270
          take care of the other elements of the Taylor form (coefficients 
1271
          errors and delta.
1272
    """
1273
    # Absolute as the default error.
1274
    if errorTypeSa is None:
1275
        errorTypeSo = sollya_lib_absolute()
1276
    elif errorTypeSa == "relative":
1277
        errorTypeSo = sollya_lib_relative()
1278
    elif errortypeSa == "absolute":
1279
        errorTypeSo = sollya_lib_absolute()
1280
    else:
1281
        # No clean up needed.
1282
        return None
1283
    # Global precision stuff
1284
    precisionChangedSa = False
1285
    currentSollyaPrecSo = pobyso_get_prec_so()
1286
    currentSollyaPrecSa = pobyso_constant_from_int_so_sa(currentSollyaPrecSo)
1287
    if not precisionSa is None:
1288
        if precisionSa > currentSollyaPrecSa:
1289
            pobyso_set_prec_sa_so(precisionSa)
1290
            precisionChangedSa = True
1291
            
1292
    if len(functionSa.variables()) > 0:
1293
        varSa = functionSa.variables()[0]
1294
        pobyso_name_free_variable_sa_so(str(varSa))
1295
    # In any case (point or interval) the parent of pointSa has a precision
1296
    # method.
1297
    pointPrecSa = pointSa.parent().precision()
1298
    if precisionSa > pointPrecSa:
1299
        pointPrecSa = precisionSa
1300
    # In any case (point or interval) pointSa has a base_ring() method.
1301
    pointBaseRingString = str(pointSa.base_ring())
1302
    if re.search('Interval', pointBaseRingString) is None: # Point
1303
        pointSo = pobyso_constant_sa_so(pointSa, pointPrecSa)
1304
    else: # Interval.
1305
        pointSo = pobyso_interval_to_range_sa_so(pointSa, pointPrecSa)
1306
    # Sollyafy the function.
1307
    functionSo = pobyso_parse_string_sa_so(functionSa._assume_str().replace('_SAGE_VAR_', ''))
1308
    if sollya_lib_obj_is_error(functionSo):
1309
        print "pobyso_tailorform: function string can't be parsed!"
1310
        return None
1311
    # Sollyafy the degree
1312
    degreeSo = sollya_lib_constant_from_int(int(degreeSa))
1313
    # Sollyafy the point
1314
    # Call Sollya
1315
    taylorFormSo = \
1316
        sollya_lib_taylorform(functionSo, degreeSo, pointSo, errorTypeSo,\
1317
                                         None)
1318
    sollya_lib_clear_obj(functionSo)
1319
    sollya_lib_clear_obj(degreeSo)
1320
    sollya_lib_clear_obj(pointSo)
1321
    sollya_lib_clear_obj(errorTypeSo)
1322
    (tfsAsList, numElements, isEndElliptic) = \
1323
            pobyso_get_list_elements_so_so(taylorFormSo)
1324
    polySo = tfsAsList[0]
1325
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
1326
    polyRealField = RealField(maxPrecision)
1327
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, polyRealField)
1328
    if precisionChangedSa:
1329
        sollya_lib_set_prec(currentSollyaPrecSo)
1330
        sollya_lib_clear_obj(currentSollyaPrecSo)
1331
    polynomialRing = polyRealField[str(varSa)]
1332
    polySa = polynomial(expSa, polynomialRing)
1333
    taylorFormSa = [polySa]
1334
    # Final clean-up.
1335
    sollya_lib_clear_obj(taylorFormSo)
1336
    return(taylorFormSa)
1337
# End pobyso_taylor_form_sa_sa
1338

    
1339
def pobyso_taylorform_so_so(functionSo, degreeSo, pointSo, intervalSo=None, \
1340
                            errorTypeSo=None):
1341
    createdErrorType = False
1342
    if errorTypeSo is None:
1343
        errorTypeSo = sollya_lib_absolute()
1344
        createdErrorType = True
1345
    else:
1346
        #TODO: deal with the other case.
1347
        pass
1348
    if intervalSo is None:
1349
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
1350
                                         errorTypeSo, None)
1351
    else:
1352
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
1353
                                         intervalSo, errorTypeSo, None)
1354
    if createdErrorType:
1355
        sollya_lib_clear_obj(errorTypeSo)
1356
    return(resultSo)
1357
        
1358

    
1359
def pobyso_univar_polynomial_print_reverse(polySa):
1360
    """ Legacy function. See pobyso_univar_polynomial_print_reverse_sa_sa. """
1361
    return(pobyso_univar_polynomial_print_reverse_sa_sa(polySa))
1362

    
1363
def pobyso_univar_polynomial_print_reverse_sa_sa(polySa):
1364
    """
1365
    Return the string representation of a univariate polynomial with
1366
    monomials ordered in the x^0..x^n order of the monomials.
1367
    Remember: Sage
1368
    """
1369
    polynomialRing = polySa.base_ring()
1370
    # A very expensive solution:
1371
    # -create a fake multivariate polynomial field with only one variable,
1372
    #   specifying a negative lexicographical order;
1373
    mpolynomialRing = PolynomialRing(polynomialRing.base(), \
1374
                                     polynomialRing.variable_name(), \
1375
                                     1, order='neglex')
1376
    # - convert the univariate argument polynomial into a multivariate
1377
    #   version;
1378
    p = mpolynomialRing(polySa)
1379
    # - return the string representation of the converted form.
1380
    # There is no simple str() method defined for p's class.
1381
    return(p.__str__())
1382
#
1383
print pobyso_get_prec()  
1384
pobyso_set_prec(165)
1385
print pobyso_get_prec()  
1386
a=100
1387
print type(a)
1388
id(a)
1389
print "Max arity: ", pobyso_max_arity
1390
print "Function tripleDouble (43) as a string: ", pobyso_function_type_as_string(43)
1391
print "Function None (44) as a string: ", pobyso_function_type_as_string(44)
1392
print "...Pobyso check done"