Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 79

Historique | Voir | Annoter | Télécharger (37,12 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_autoprint(arg):
89
    sollya_lib_autoprint(arg,None)
90

    
91
def pobyso_autoprint_so_so(arg):
92
    sollya_lib_autoprint(arg,None)
93
    
94
def pobyso_absolute_so_so():
95
    return(sollya_lib_absolute(None))
96

    
97
def pobyso_build_function_sub_so_so(exp1So, exp2So):
98
    return(sollya_lib_build_function_sub(exp1So, exp2So))
99

    
100
def pobyso_cmp(rnArg, soCte):
101
    """
102
    Compare the MPFR value a RealNumber with that of a Sollya constant.
103
    
104
    Get the value of the Sollya constant into a RealNumber and compare
105
    using MPFR. Could be optimized by working directly with a mpfr_t
106
    for the intermediate number. 
107
    """
108
    precisionOfCte = c_int(0)
109
    # From the Sollya constant, create a local Sage RealNumber.
110
    sollya_lib_get_prec_of_constant(precisionOfCte, soCte) 
111
    #print "Precision of constant: ", precisionOfCte
112
    RRRR = RealField(precisionOfCte.value)
113
    rnLocal = RRRR(0)
114
    sollya_lib_get_constant(get_rn_value(rnLocal), soCte)
115
    #print "rnDummy: ", rnDummy
116
    # Compare the local Sage RealNumber with rnArg.
117
    return(cmp_rn_value(rnArg, rnLocal))
118

    
119
def pobyso_change_var_in_function_so_so(funcSo, chvarExpSo):
120
    return(sollya_lib_evaluate(funcSo,chvarExpSo))
121
    
122

    
123
def pobyso_chebyshevform_so_so(functionSo, degreeSo, intervalSo):
124
    resultSo = sollya_lib_chebyshevform(functionSo, degreeSo, intervalSo)
125
    return(resultSo)
126
        
127

    
128

    
129
def pobyso_compute_pos_function_abs_val_bounds_sa_sa(funcSa, lowerBoundSa, \
130
                                                     upperBoundSa):
131
    """
132
    TODO: set the variable name in Sollya.
133
    """
134
    funcSo = pobyso_parse_string(funcSa._assume_str())
135
    rangeSo = pobyso_range_sa_so(lowerBoundSa, upperBoundSa)
136
    infnormSo = pobyso_infnorm_so_so(funcSo,rangeSo)
137
    fMaxSa = pobyso_get_interval_from_range_so_sa(infnormSo)
138
    # Get the top bound and compute the binade top limit.
139
    fMaxUpperBoundSa = fMaxSa.upper()
140
    binadeTopLimitSa = 2**ceil(fMaxUpperBoundSa.log2())
141
    # Put up together the function to use to compute the lower bound.
142
    funcAuxSo = pobyso_parse_string(str(binadeTopLimitSa) +  \
143
                                    '-(' + f._assume_str() + ')')
144
    pobyso_autoprint(funcAuxSo)
145
    # Clear the Sollay range before a new call to infnorm and issue the call.
146
    sollya_lib_clear_obj(infnormSo)
147
    infnormSo = pobyso_infnorm_so_so(funcAuxSo,rangeSo)
148
    fMinSa = pobyso_get_interval_from_range_so_sa(infnormSo)
149
    sollya_lib_clear_obj(infnormSo)
150
    fMinLowerBoundSa = topBinadeLimit - fMinSa.lower()
151
    # Compute the maximum of the precisions of the different bounds.
152
    maxPrecSa = max([fMinLowerBoundSa.parent().precision(), \
153
                     fMaxUpperBoundSa.parent().precision()])
154
    # Create a RealIntervalField and create an interval with the "good" bounds.
155
    RRRI = RealIntervalField(maxPrecSa)
156
    imageIntervalSa = RRRI(fMinLowerBoundSa, fMaxUpperBoundSa)
157
    # Free the uneeded Sollya objects
158
    sollya_lib_clear_obj(funcSo)
159
    sollya_lib_clear_obj(funcAuxSo)
160
    sollya_lib_clear_obj(rangeSo)
161
    return(imageIntervalSa)
162
    # End pobyso_compute_function_abs_val_bounds_sa_sa
163

    
164
def pobyso_constant(rnArg):
165
    """ Legacy function. See pobyso_constant_sa_so. """
166
    return(pobyso_constant_sa_so(rnArg))
167
    
168
def pobyso_constant_sa_so(rnArg):
169
    """
170
    Create a Sollya constant from a RealNumber.
171
    """
172
    return(sollya_lib_constant(get_rn_value(rnArg)))
173
    
174
def pobyso_constant_0_sa_so():
175
    return(pobyso_constant_from_int_sa_so(0))
176

    
177
def pobyso_constant_1():
178
    """ Legacy function. See pobyso_constant_so_so. """
179
    return(pobyso_constant_1_sa_so())
180

    
181
def pobyso_constant_1_sa_so():
182
    return(pobyso_constant_from_int_sa_so(1))
183

    
184
def pobyso_constant_from_int(anInt):
185
    """ Legacy function. See pobyso_constant_from_int_sa_so. """
186
    return(pobyso_constant_from_int_sa_so(anInt))
187

    
188
def pobyso_constant_from_int_sa_so(anInt):
189
    return(sollya_lib_constant_from_int(int(anInt)))
190

    
191
def pobyso_function_type_as_string(funcType):
192
    """ Legacy function. See pobyso_function_type_as_string_so_sa. """
193
    return(pobyso_function_type_as_string_so_sa(funcType))
194

    
195
def pobyso_function_type_as_string_so_sa(funcType):
196
    """
197
    Numeric Sollya function codes -> Sage mathematical function names.
198
    Notice that pow -> ^ (a la Sage, not a la Python).
199
    """
200
    if funcType == SOLLYA_BASE_FUNC_ABS:
201
        return "abs"
202
    elif funcType == SOLLYA_BASE_FUNC_ACOS:
203
        return "arccos"
204
    elif funcType == SOLLYA_BASE_FUNC_ACOSH:
205
        return "arccosh"
206
    elif funcType == SOLLYA_BASE_FUNC_ADD:
207
        return "+"
208
    elif funcType == SOLLYA_BASE_FUNC_ASIN:
209
        return "arcsin"
210
    elif funcType == SOLLYA_BASE_FUNC_ASINH:
211
        return "arcsinh"
212
    elif funcType == SOLLYA_BASE_FUNC_ATAN:
213
        return "arctan"
214
    elif funcType == SOLLYA_BASE_FUNC_ATANH:
215
        return "arctanh"
216
    elif funcType == SOLLYA_BASE_FUNC_CEIL:
217
        return "ceil"
218
    elif funcType == SOLLYA_BASE_FUNC_CONSTANT:
219
        return "cte"
220
    elif funcType == SOLLYA_BASE_FUNC_COS:
221
        return "cos"
222
    elif funcType == SOLLYA_BASE_FUNC_COSH:
223
        return "cosh"
224
    elif funcType == SOLLYA_BASE_FUNC_DIV:
225
        return "/"
226
    elif funcType == SOLLYA_BASE_FUNC_DOUBLE:
227
        return "double"
228
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEDOUBLE:
229
        return "doubleDouble"
230
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEEXTENDED:
231
        return "doubleDxtended"
232
    elif funcType == SOLLYA_BASE_FUNC_ERF:
233
        return "erf"
234
    elif funcType == SOLLYA_BASE_FUNC_ERFC:
235
        return "erfc"
236
    elif funcType == SOLLYA_BASE_FUNC_EXP:
237
        return "exp"
238
    elif funcType == SOLLYA_BASE_FUNC_EXP_M1:
239
        return "expm1"
240
    elif funcType == SOLLYA_BASE_FUNC_FLOOR:
241
        return "floor"
242
    elif funcType == SOLLYA_BASE_FUNC_FREE_VARIABLE:
243
        return "freeVariable"
244
    elif funcType == SOLLYA_BASE_FUNC_HALFPRECISION:
245
        return "halfPrecision"
246
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYCONSTANT:
247
        return "libraryConstant"
248
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYFUNCTION:
249
        return "libraryFunction"
250
    elif funcType == SOLLYA_BASE_FUNC_LOG:
251
        return "log"
252
    elif funcType == SOLLYA_BASE_FUNC_LOG_10:
253
        return "log10"
254
    elif funcType == SOLLYA_BASE_FUNC_LOG_1P:
255
        return "log1p"
256
    elif funcType == SOLLYA_BASE_FUNC_LOG_2:
257
        return "log2"
258
    elif funcType == SOLLYA_BASE_FUNC_MUL:
259
        return "*"
260
    elif funcType == SOLLYA_BASE_FUNC_NEARESTINT:
261
        return "round"
262
    elif funcType == SOLLYA_BASE_FUNC_NEG:
263
        return "__neg__"
264
    elif funcType == SOLLYA_BASE_FUNC_PI:
265
        return "pi"
266
    elif funcType == SOLLYA_BASE_FUNC_POW:
267
        return "^"
268
    elif funcType == SOLLYA_BASE_FUNC_PROCEDUREFUNCTION:
269
        return "procedureFunction"
270
    elif funcType == SOLLYA_BASE_FUNC_QUAD:
271
        return "quad"
272
    elif funcType == SOLLYA_BASE_FUNC_SIN:
273
        return "sin"
274
    elif funcType == SOLLYA_BASE_FUNC_SINGLE:
275
        return "single"
276
    elif funcType == SOLLYA_BASE_FUNC_SINH:
277
        return "sinh"
278
    elif funcType == SOLLYA_BASE_FUNC_SQRT:
279
        return "sqrt"
280
    elif funcType == SOLLYA_BASE_FUNC_SUB:
281
        return "-"
282
    elif funcType == SOLLYA_BASE_FUNC_TAN:
283
        return "tan"
284
    elif funcType == SOLLYA_BASE_FUNC_TANH:
285
        return "tanh"
286
    elif funcType == SOLLYA_BASE_FUNC_TRIPLEDOUBLE:
287
        return "tripleDouble"
288
    else:
289
        return None
290

    
291
def pobyso_get_constant(rnArg, soConst):
292
    """ Legacy function. See pobyso_get_constant_so_sa. """
293
    return(pobyso_get_constant_so_sa(rnArg, soConst))
294

    
295
def pobyso_get_constant_so_sa(rnArg, soConst):
296
    """
297
    Set the value of rnArg to the value of soConst in MPFR_RNDN mode.
298
    rnArg must already exist and belong to some RealField.
299
    We assume that soConst points to a Sollya constant.
300
    """
301
    return(sollya_lib_get_constant(get_rn_value(rnArg), soConst))
302
    
303
def pobyso_get_constant_as_rn(ctExpSo):
304
    """ Legacy function. See pobyso_get_constant_as_rn_so_sa. """ 
305
    return(pobyso_get_constant_as_rn_so_sa(ctExpSo))
306
    
307
def pobyso_get_constant_as_rn_so_sa(constExpSo):
308
    precisionSa  = pobyso_get_prec_of_constant_so_sa(constExpSo) 
309
    RRRR = RealField(precisionSa)
310
    rnSa = RRRR(0)
311
    sollya_lib_get_constant(get_rn_value(rnSa), constExpSo)
312
    return(rnSa)
313

    
314
def pobyso_get_constant_as_rn_with_rf(ctExp, realField):
315
    """ Legacy function. See pobyso_get_constant_as_rn_with_rf_so_sa."""
316
    return(pobyso_get_constant_as_rn_with_rf_so_sa(ctExp, realField))
317
    
318
def pobyso_get_constant_as_rn_with_rf_so_sa(ctExpSo, realFieldSa = None):
319
    if realFieldSa is None:
320
        sollyaPrecSa = pobyso_get_prec_so_sa()
321
        realFieldSa = RealField(sollyaPrecSa)
322
    rnSa = realFieldSa(0)
323
    sollya_lib_get_constant(get_rn_value(rnSa), ctExpSo)
324
    return(rnSa)
325

    
326
def pobyso_get_free_variable_name():
327
    """ Legacy function. See pobyso_get_free_variable_name_so_sa."""
328
    return(pobyso_get_free_variable_name_so_sa())
329

    
330
def pobyso_get_free_variable_name_so_sa():
331
    return(sollya_lib_get_free_variable_name())
332
    
333
def pobyso_get_function_arity(expressionSo):
334
    """ Legacy function. See pobyso_get_function_arity_so_sa."""
335
    return(pobyso_get_function_arity_so_sa(expressionSo))
336

    
337
def pobyso_get_function_arity_so_sa(expressionSo):
338
    arity = c_int(0)
339
    sollya_lib_get_function_arity(byref(arity),expressionSo)
340
    return(int(arity.value))
341

    
342
def pobyso_get_head_function(expressionSo):
343
    """ Legacy function. See pobyso_get_head_function_so_sa. """
344
    return(pobyso_get_head_function_so_sa(expressionSo)) 
345

    
346
def pobyso_get_head_function_so_sa(expressionSo):
347
    functionType = c_int(0)
348
    sollya_lib_get_head_function(byref(functionType), expressionSo, None)
349
    return(int(functionType.value))
350

    
351
def pobyso_get_interval_from_range_so_sa(soRange, realIntervalFieldSa = None ):
352
    """
353
    Return the Sage interval corresponding to the Sollya range argument.
354
    If no reaInterval lField is passed as argument, the interval bounds are not
355
    rounded: they are elements of RealIntervalField of the "right" precision
356
    to hold all the digits.
357
    """
358
    prec = c_int(0)
359
    if realIntervalFieldSa is None:
360
        retval = sollya_lib_get_prec_of_range(byref(prec), soRange, None)
361
        if retval == 0:
362
            return(None)
363
        realIntervalFieldSa = RealIntervalField(prec.value)
364
    intervalSa = realIntervalFieldSa(0,0)
365
    retval = \
366
        sollya_lib_get_interval_from_range(get_interval_value(intervalSa),\
367
                                           soRange)
368
    if retval == 0:
369
        return(None)
370
    return(intervalSa)
371
# End pobyso_get_interval_from_range_so_sa
372

    
373
def pobyso_get_list_elements(soObj):
374
    """ Legacy function. See pobyso_get_list_elements_so_so. """
375
    return(pobyso_get_list_elements_so_so(soObj))
376
 
377
def pobyso_get_list_elements_so_so(soObj):
378
    """
379
    Get the list elements as a Sage/Python array of Sollya objects.
380
    The other data returned are also Sage/Python objects.
381
    """
382
    listAddress = POINTER(c_longlong)()
383
    numElements = c_int(0)
384
    isEndElliptic = c_int(0)
385
    listAsList = []
386
    result = sollya_lib_get_list_elements(byref(listAddress),\
387
                                          byref(numElements),\
388
                                          byref(isEndElliptic),\
389
                                          soObj)
390
    if result == 0 :
391
        return None
392
    for i in xrange(0, numElements.value, 1):
393
        listAsList.append(listAddress[i])
394
    return(listAsList, numElements.value, isEndElliptic.value)
395

    
396
def pobyso_get_max_prec_of_exp(soExp):
397
    """ Legacy function. See pobyso_get_max_prec_of_exp_so_sa. """
398
    return(pobyso_get_max_prec_of_exp_so_sa(soExp))
399

    
400
def pobyso_get_max_prec_of_exp_so_sa(soExp):
401
    """
402
    Get the maximum precision used for the numbers in a Sollya expression.
403
    
404
    Arguments:
405
    soExp -- a Sollya expression pointer
406
    Return value:
407
    A Python integer
408
    TODO: 
409
    - error management;
410
    - correctly deal with numerical type such as DOUBLEEXTENDED.
411
    """
412
    maxPrecision = 0
413
    minConstPrec = 0
414
    currentConstPrec = 0
415
    operator = pobyso_get_head_function_so_sa(soExp)
416
    if (operator != SOLLYA_BASE_FUNC_CONSTANT) and \
417
    (operator != SOLLYA_BASE_FUNC_FREE_VARIABLE):
418
        (arity, subexpressions) = pobyso_get_subfunctions_so_sa(soExp)
419
        for i in xrange(arity):
420
            maxPrecisionCandidate = \
421
                pobyso_get_max_prec_of_exp_so_sa(subexpressions[i])
422
            if maxPrecisionCandidate > maxPrecision:
423
                maxPrecision = maxPrecisionCandidate
424
        return(maxPrecision)
425
    elif operator == SOLLYA_BASE_FUNC_CONSTANT:
426
        minConstPrec = pobyso_get_min_prec_of_constant_so_sa(soExp)
427
        #currentConstPrec = pobyso_get_min_prec_of_constant_so_sa(soExp)
428
        #print minConstPrec, " - ", currentConstPrec 
429
        return(pobyso_get_min_prec_of_constant_so_sa(soExp))
430
    
431
    elif operator == SOLLYA_BASE_FUNC_FREE_VARIABLE:
432
        return(0)
433
    else:
434
        print "pobyso_get_max_prec_of_exp_so_sa: unexepected operator."
435
        return(0)
436

    
437
def pobyso_get_min_prec_of_constant_so_sa(soConstExp):
438
    """
439
    Get the minimum precision necessary to represent the value of a Sollya
440
    constant.
441
    MPFR_MIN_PREC and powers of 2 are taken into account.
442
    We assume that soCteExp is a point
443
    """
444
    constExpAsRn = pobyso_get_constant_as_rn_so_sa(soConstExp)
445
    return(min_mpfr_size(get_rn_value(constExpAsRn)))
446

    
447
def pobyso_get_sage_exp_from_sollya_exp(sollyaExp, realField = RR):
448
    """ Legacy function. See pobyso_get_sage_exp_from_sollya_exp_so_sa. """
449
    return(pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExp, realField = RR))
450

    
451
def pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExp, realField = RR):
452
    """
453
    Get a Sage expression from a Sollya expression. 
454
    Currently only tested with polynomials with floating-point coefficients.
455
    Notice that, in the returned polynomial, the exponents are RealNumbers.
456
    """
457
    #pobyso_autoprint(sollyaExp)
458
    operator = pobyso_get_head_function_so_sa(sollyaExp)
459
    # Constants and the free variable are special cases.
460
    # All other operator are dealt with in the same way.
461
    if (operator != SOLLYA_BASE_FUNC_CONSTANT) and \
462
       (operator != SOLLYA_BASE_FUNC_FREE_VARIABLE):
463
        (arity, subexpressions) = pobyso_get_subfunctions_so_sa(sollyaExp)
464
        if arity == 1:
465
            sageExp = eval(pobyso_function_type_as_string_so_sa(operator) + \
466
            "(" + pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressions[0], \
467
            realField) + ")")
468
        elif arity == 2:
469
            # We do not get through the preprocessor.
470
            # The "^" operator is then a special case.
471
            if operator == SOLLYA_BASE_FUNC_POW:
472
                operatorAsString = "**"
473
            else:
474
                operatorAsString = \
475
                    pobyso_function_type_as_string_so_sa(operator)
476
            sageExp = \
477
              eval("pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressions[0], realField)"\
478
              + " " + operatorAsString + " " + \
479
                   "pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressions[1], realField)")
480
        # We do not know yet how to deal with arity >= 3 
481
        # (is there any in Sollya anyway?).
482
        else:
483
            sageExp = eval('None')
484
        return(sageExp)
485
    elif operator == SOLLYA_BASE_FUNC_CONSTANT:
486
        #print "This is a constant"
487
        return pobyso_get_constant_as_rn_with_rf_so_sa(sollyaExp, realField)
488
    elif operator == SOLLYA_BASE_FUNC_FREE_VARIABLE:
489
        #print "This is free variable"
490
        return(eval(sollya_lib_get_free_variable_name()))
491
    else:
492
        print "Unexpected"
493
        return eval('None')
494
# End pobyso_get_sage_poly_from_sollya_poly
495

    
496
def pobyso_get_poly_sa_so(polySo, realFieldSa=None):
497
    pass
498
# pobyso_get_poly_sa_so
499

    
500
def pobyso_get_poly_so_sa(polySo, realFieldSa=None):
501
    """
502
    Convert a Sollya polynomial into a Sage polynomial.
503
    We assume that the polynomial is in canonical form.
504
    If no realField is given, a RealField corresponding to the maximum precision 
505
    of the coefficients is internally computed.
506
    It is not returned but can be easily retrieved from the polynomial itself.
507
    Main steps:
508
    - (optional) compute the RealField of the coefficients;
509
    - convert the Sollya expression into a Sage expression;
510
    - convert the Sage expression into a Sage polynomial
511
    TODO: the canonical thing for the polynomial.
512
    """    
513
    if realFieldSa is None:
514
        expressionPrecSa = pobyso_get_max_prec_of_exp_so_sa(polySo)
515
        realFieldSa = RealField(expressionPrecSa)
516
    #print "Sollya expression before...",
517
    #pobyso_autoprint(polySo)
518

    
519
    expressionSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, \
520
                                                             realFieldSa)
521
    #print "...Sollya expression after.",
522
    #pobyso_autoprint(polySo)
523
    polyVariableSa = expressionSa.variables()[0]
524
    polyRingSa = realFieldSa[str(polyVariableSa)]
525
    print polyRingSa
526
    # Do not use the polynomial(expressionSa, ring=polyRingSa) form!
527
    polynomialSa = polyRingSa(expressionSa)
528
    return(polynomialSa)
529
# End pobyso_get_sage_poly_from_sollya_poly
530

    
531
def pobyso_get_subfunctions(expressionSo):
532
    """ Legacy function. See pobyso_get_subfunctions_so_sa. """
533
    return(pobyso_get_subfunctions_so_sa(expressionSo)) 
534

    
535
def pobyso_get_subfunctions_so_sa(expressionSo):
536
    """
537
    Get the subfunctions of an expression.
538
    Return the number of subfunctions and the list of subfunctions addresses.
539
    S.T.: Could not figure out another way than that ugly list of declarations
540
    to recover the addresses of the subfunctions. 
541
    """
542
    subf0 = c_int(0)
543
    subf1 = c_int(0)
544
    subf2 = c_int(0)
545
    subf3 = c_int(0)
546
    subf4 = c_int(0)
547
    subf5 = c_int(0)
548
    subf6 = c_int(0)
549
    subf7 = c_int(0)
550
    subf8 = c_int(0)
551
    arity = c_int(0)
552
    nullPtr = POINTER(c_int)()
553
    sollya_lib_get_subfunctions(expressionSo, byref(arity), \
554
    byref(subf0), byref(subf1), byref(subf2), byref(subf3), byref(subf4), byref(subf5),\
555
     byref(subf6), byref(subf7), byref(subf8), nullPtr, None) 
556
#    byref(cast(subfunctions[0], POINTER(c_int))), byref(cast(subfunctions[0], POINTER(c_int))), \
557
#    byref(cast(subfunctions[2], POINTER(c_int))), byref(cast(subfunctions[3], POINTER(c_int))), \
558
#    byref(cast(subfunctions[4], POINTER(c_int))), byref(cast(subfunctions[5], POINTER(c_int))), \
559
#    byref(cast(subfunctions[6], POINTER(c_int))), byref(cast(subfunctions[7], POINTER(c_int))), \
560
#    byref(cast(subfunctions[8], POINTER(c_int))), nullPtr)
561
    subfunctions = [subf0, subf1, subf2, subf3, subf4, subf5, subf6, subf7, subf8]
562
    subs = []
563
    if arity.value > pobyso_max_arity:
564
        return(0,[])
565
    for i in xrange(arity.value):
566
        subs.append(int(subfunctions[i].value))
567
        #print subs[i]
568
    return(int(arity.value), subs)
569
    
570
def pobyso_get_prec():
571
    """ Legacy function. See pobyso_get_prec_so_sa(). """
572
    return(pobyso_get_prec_so_sa())
573

    
574
def pobyso_get_prec_so_sa():
575
    """
576
    Get the current default precision in Sollya.
577
    The return value is Sage/Python int.
578
    """
579
    precSo = sollya_lib_get_prec(None)
580
    precSa = c_int(0)
581
    sollya_lib_get_constant_as_int(byref(precSa), precSo)
582
    sollya_lib_clear_obj(precSo)
583
    return(int(precSa.value))
584

    
585
def pobyso_get_prec_of_constant(ctExpSo):
586
    """ Legacy function. See pobyso_get_prec_of_constant_so_sa. """
587
    return(pobyso_get_prec_of_constant_so_sa(ctExpSo))
588

    
589
def pobyso_get_prec_of_constant_so_sa(ctExpSo):
590
    prec = c_int(0)
591
    retc = sollya_lib_get_prec_of_constant(byref(prec), ctExpSo, None)
592
    if retc == 0:
593
        return(None)
594
    return(int(prec.value))
595

    
596
def pobyso_get_prec_of_range_so_sa(rangeSo):
597
    prec = c_int(0)
598
    retc = sollya_lib_get_prec_of_range(byref(prec), rangeSo, None)
599
    if retc == 0:
600
        return(None)
601
    return(int(prec.value))
602

    
603
def pobyso_infnorm_so_so(func, interval, file = None, intervalList = None):
604
    print "Do not use this function. User pobyso_supnorm_so_so instead."
605
    return(None)
606

    
607
def pobyso_lib_init():
608
    sollya_lib_init(None)
609
    
610
def pobyso_name_free_variable(freeVariableName):
611
    """ Legacy function. See pobyso_name_free_variable_sa_so. """
612
    pobyso_name_free_variable_sa_so(freeVariableName)
613

    
614
def pobyso_name_free_variable_sa_so(freeVariableName):
615
    sollya_lib_name_free_variable(freeVariableName)
616

    
617
def pobyso_parse_string(string):
618
    """ Legacy function. See pobyso_parse_string_sa_so. """
619
    return(pobyso_parse_string_sa_so(string))
620
 
621
def pobyso_parse_string_sa_so(string):
622
    return(sollya_lib_parse_string(string))
623

    
624
def pobyso_range(rnLowerBound, rnUpperBound):
625
    """ Legacy function. See pobyso_range_sa_so. """
626
    return(pobyso_range_sa_so(rnLowerBound, rnUpperBound)) 
627

    
628
def pobyso_bounds_to_range_sa_so(rnLowerBoundSa, rnUpperBoundSa):
629
    """
630
    Return a Sollya range from to 2 RealField elements.
631
    The Sollya range element has a sufficient precision to hold all
632
    the digits of the bounds.
633
    """
634
    if rnLowerBoundSa > rnUpperBoundSa:
635
        return None
636
    lbPrec = rnLowerBoundSa.parent().precision()
637
    ubPrec = rnLowerBoundSa.parent().precision()
638
    currentSollyaPrecSa = pobyso_get_prec_so_sa()
639
    maxPrecSa = max(lbPrec, ubPrec, currentSollyaPrecSa)
640
    # Change the current Sollya precision only if necessary.
641
    if maxPrecSa > currentSollyaPrecSa:
642
        currentPrecSo = sollya_lib_get_prec(None)
643
        newPrecSo = solly_lib_constant_from_uint64(maxPrecSa)
644
        sollya_lib_set_prec(newPrecSo)
645
    lowerBoundSo = sollya_lib_constant(get_rn_value(rnLowerBoundSa))
646
    upperBoundSo = sollya_lib_constant(get_rn_value(rnUpperBoundSa))
647
    rangeSo = sollya_lib_range(lowerBoundSo, upperBoundSo)
648
    currentPrecSo = sollya_lib_get_prec(None)
649
    if maxPrecSa > currentSollyaPrecSa:
650
        sollya_lib_set_prec(currentPrecSo)
651
        sollya_lib_clear_obj(currentPrecSo)
652
        sollya_lib_clear_obj(newPrecSo)
653
    sollya_lib_clear_obj(lowerBoundSo)
654
    sollya_lib_clear_obj(upperBoundSo)
655
    return(rangeSo)
656

    
657
def pobyso_range_to_interval_so_sa(rangeSo, realIntervalField = None):
658
    if realIntervalField is None:
659
        precSa = pobyso_get_prec_of_range_so_sa(rangeSo)
660
        realIntervalField = RealIntervalField(precSa)
661
    intervalSa = \
662
        pobyso_get_interval_from_range_so_sa(rangeSo, realIntervalField) 
663
    return(intervalSa)
664

    
665
def pobyso_remez_canonical_sa_sa(func, \
666
                                 degree, \
667
                                 lowerBound, \
668
                                 upperBound, \
669
                                 weight = None, \
670
                                 quality = None):
671
    """
672
    All arguments are Sage/Python.
673
    The functions (func and weight) must be passed as expressions or strings.
674
    Otherwise the function fails. 
675
    The return value is a pointer is a Sage polynomial.
676
    """
677
    var('zorglub')    # Dummy variable name for type check only.
678
    polySo = pobyso_remez_canonical_sa_so(func, \
679
                                 degree, \
680
                                 lowerBound, \
681
                                 upperBound, \
682
                                 weight = None, \
683
                                 quality = None)
684
    if parent(func) == parent("string"):
685
        functionSa = eval(func)
686
    # Expression test.
687
    elif type(func) == type(zorglub):
688
        functionSa = func
689
    maxPrecision = 0
690
    if polySo is None:
691
        return(None)
692
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
693
    RRRR = RealField(maxPrecision)
694
    polynomialRing = RRRR[functionSa.variables()[0]]
695
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, RRRR)
696
    polySa = polynomial(expSa, polynomialRing)
697
    return(polySa)
698
    
699
def pobyso_remez_canonical(func, \
700
                           degree, \
701
                           lowerBound, \
702
                           upperBound, \
703
                           weight = "1", \
704
                           quality = None):
705
    """ Legacy function. See pobyso_remez_canonical_sa_so. """
706
    return(pobyso_remez_canonical_sa_so(func, \
707
                                        degree, \
708
                                        lowerBound, \
709
                                        upperBound, \
710
                                        weight, \
711
                                        quality))
712
def pobyso_remez_canonical_sa_so(func, \
713
                                 degree, \
714
                                 lowerBound, \
715
                                 upperBound, \
716
                                 weight = None, \
717
                                 quality = None):
718
    """
719
    All arguments are Sage/Python.
720
    The functions (func and weight) must be passed as expressions or strings.
721
    Otherwise the function fails. 
722
    The return value is a pointer to a Sollya function.
723
    """
724
    var('zorglub')    # Dummy variable name for type check only.
725
    currentVariableName = None
726
    # The func argument can be of different types (string, 
727
    # symbolic expression...)
728
    if parent(func) == parent("string"):
729
        functionSo = sollya_lib_parse_string(func)
730
    # Expression test.
731
    elif type(func) == type(zorglub):
732
        # Until we are able to translate Sage expressions into Sollya 
733
        # expressions : parse the string version.
734
        currentVariableName = func.variables()[0]
735
        sollya_lib_name_free_variable(str(currentVariableName))
736
        functionSo = sollya_lib_parse_string(func._assume_str())
737
    else:
738
        return(None)
739
    if weight is None:
740
        weightSo = pobyso_constant_1_sa_so()
741
    elif parent(weight) == parent("string"):
742
        weightSo = sollya_lib_parse_string(func)
743
    elif type(weight) == type(zorglub): 
744
        functionSo = sollya_lib_parse_string_sa_so(weight._assume_str())
745
    else:
746
        return(None)
747
    degreeSo = pobyso_constant_from_int(degree)
748
    rangeSo = pobyso_range_sa_so(lowerBound, upperBound)
749
    if not quality is None:
750
        qualitySo= pobyso_constant_sa_so(quality)
751
    else:
752
        qualitySo = None
753
    return(sollya_lib_remez(functionSo, \
754
                            degreeSo, \
755
                            rangeSo, \
756
                            weightSo, \
757
                            qualitySo, \
758
                            None))
759
    
760
def pobyso_remez_canonical_so_so(funcSo, \
761
                                 degreeSo, \
762
                                 rangeSo, \
763
                                 weightSo = pobyso_constant_1_sa_so(),\
764
                                 qualitySo = None):
765
    """
766
    All arguments are pointers to Sollya objects.
767
    The return value is a pointer to a Sollya function.
768
    """
769
    if not sollya_lib_obj_is_function(funcSo):
770
        return(None)
771
    return(sollya_lib_remez(funcSo, degreeSo, rangeSo, weightSo, qualitySo, None))
772
    
773
def pobyso_set_canonical_off():
774
    sollya_lib_set_canonical(sollya_lib_off())
775

    
776
def pobyso_set_canonical_on():
777
    sollya_lib_set_canonical(sollya_lib_on())
778

    
779
def pobyso_set_prec(p):
780
    """ Legacy function. See pobyso_set_prec_sa_so. """
781
    return( pobyso_set_prec_sa_so(p))
782

    
783
def pobyso_set_prec_sa_so(p):
784
    a = c_int(p)
785
    precSo = c_void_p(sollya_lib_constant_from_int(a))
786
    sollya_lib_set_prec(precSo)
787

    
788
def pobyso_supnorm_so_so(polySo, funcSo, intervalSo, errorTypeSo, accuracySo):
789
    return(sollya_lib_supnorm(polySo, funcSo, intervalSo, errorTypeSo, \
790
                              accuracySo))
791

    
792
def pobyso_taylor_expansion_with_change_var_so_so(functionSo, degreeSo, rangeSo, \
793
                                                errorTypeSo, \
794
                                                sollyaPrecSo=None):
795
    """
796
    Compute the Taylor expansion with the variable change
797
    x -> (x-intervalCenter) included.
798
    """
799
    # No global change of the working precision.
800
    if not sollyaPrecSo is None:
801
        initialPrecSo = sollya_lib_get_prec(None)
802
        sollya_lib_set_prec(sollyaPrecSo)
803
    #
804
    intervalCenterSo = sollya_lib_mid(rangeSo)
805
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo, \
806
                                         intervalCenterSo, \
807
                                         rangeSo, errorTypeSo, None)
808
    (taylorFormListSo, numElements, isEndElliptic) = \
809
        pobyso_get_list_elements_so_so(taylorFormSo)
810
    polySo = taylorFormListSo[0]
811
    errorRangeSo = taylorFormListSo[2]
812
    maxErrorSo = sollya_lib_sup(errorRangeSo)
813
    changeVarExpSo = sollya_lib_build_function_sub(\
814
                       sollya_lib_build_function_free_variable(),\
815
                       sollya_lib_copy_obj(intervalCenterSo))
816
    polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpSo) 
817
    # If changed, reset the Sollya working precision.
818
    if not sollyaPrecSo is None:
819
        sollya_lib_set_prec(initialPrecSo)
820
        sollya_lib_clear_obj(initailPrecSo)
821
    return((polyVarChangedSo, intervalCenterSo, maxErrorSo))
822
# end pobyso_taylor_expansion_with_change_var_so_so
823

    
824
def pobyso_taylor_expansion_no_change_var_so_so(functionSo, degreeSo, rangeSo, \
825
                                                errorTypeSo, \
826
                                                sollyaPrecSo=None):
827
    """
828
    Compute the Taylor expansion without the variable change
829
    x -> x-intervalCenter.
830
    """
831
    # No global change of the working precision.
832
    if not sollyaPrecSo is None:
833
        initialPrecSo = sollya_lib_get_prec(None)
834
        sollya_lib_set_prec(sollyaPrecSo)
835
    #
836
    intervalCenterSo = sollya_lib_mid(rangeSo)
837
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo, \
838
                                         intervalCenterSo, \
839
                                         rangeSo, errorTypeSo, None)
840
    (taylorFormListSo, numElements, isEndElliptic) = \
841
        pobyso_get_list_elements_so_so(taylorFormSo)
842
    polySo = taylorFormListSo[0]
843
    errorRangeSo = taylorFormListSo[2]
844
    maxErrorSo = sollya_lib_sup(errorRangeSo)
845
    # If changed, reset the Sollya working precision.
846
    if not sollyaPrecSo is None:
847
        sollya_lib_set_prec(initialPrecSo)
848
        sollya_lib_clear_obj(initialPrecSo)
849
    return((polySo, intervalCenterSo, maxErrorSo))
850
# end pobyso_taylor_expansion_no_change_var_so_so
851

    
852
def pobyso_taylor(function, degree, point):
853
    """ Legacy function. See pobysoTaylor_so_so. """
854
    return(pobyso_taylor_so_so(function, degree, point))
855

    
856
def pobyso_taylor_so_so(functionSo, degreeSo, pointSo):
857
    return(sollya_lib_taylor(functionSo, degreeSo, pointSo))
858
    
859
def pobyso_taylorform(function, degree, point = None, interval = None, errorType=None):
860
    """ Legacy function. See ;"""
861
    
862
def pobyso_taylorform_sa_sa(functionSa, \
863
                            degree, \
864
                            point, \
865
                            precision, \
866
                            interval=None, \
867
                            errorType=None):
868
    """
869
    Compute the Taylor form of 'degree' for 'functionSa' at 'point' 
870
    for 'interval' with 'errorType'. 
871
    point: must be a Real or a Real interval.
872
    return the Taylor form as an array
873
    TODO: take care of the interval and of point when it is an interval;
874
          when errorType is not None;
875
          take care of the other elements of the Taylor form (coefficients errors and
876
          delta.
877
    """
878
    # Absolute as the default error.
879
    if errorType is None:
880
        errorTypeSo = sollya_lib_absolute()
881
    else:
882
        #TODO: deal with the other case.
883
        pass
884
    varSa = functionSa.variables()[0]
885
    pointBaseRingString = str(point.base_ring())
886
    if not re.search('Real', pointBaseRingString):
887
        return None
888
    # Call Sollya but first "sollyafy" the arguments.
889
    pobyso_name_free_variable_sa_so(str(varSa))
890
    #pobyso_set_prec_sa_so(300)
891
    # Sollyafy the function.
892
    functionSo = pobyso_parse_string_sa_so(functionSa._assume_str())
893
    if sollya_lib_obj_is_error(functionSo):
894
        print "pobyso_tailorform: function string can't be parsed!"
895
        return None
896
    # Sollyafy the degree
897
    degreeSo = sollya_lib_constant_from_int(int(degree))
898
    # Sollyafy the point
899
    if not re.search('Interval', pointBaseRingString):
900
        pointSo  = pobyso_constant_sa_so(point)
901
    else:
902
        # TODO: deal with the interval case.
903
        pass
904
    # Call Sollya
905
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, errorTypeSo,\
906
                                         None)
907
    (tfsAsList, numElements, isEndElliptic) = \
908
            pobyso_get_list_elements_so_so(taylorFormSo)
909
    polySo = tfsAsList[0]
910
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
911
    polyRealField = RealField(maxPrecision)
912
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, polyRealField)
913
    sollya_lib_close()
914
    polynomialRing = polyRealField[str(varSa)]
915
    polySa = polynomial(expSa, polynomialRing)
916
    taylorFormSa = [polySa]
917
    return(taylorFormSa)
918
# End pobyso_taylor_form_sa_sa
919

    
920
def pobyso_taylorform_so_so(functionSo, degreeSo, pointSo, intervalSo=None, \
921
                            errorTypeSo=None):
922
    createdErrorType = False
923
    if errorTypeSo is None:
924
        errorTypeSo = sollya_lib_absolute()
925
        createdErrorType = True
926
    else:
927
        #TODO: deal with the other case.
928
        pass
929
    if intervalSo is None:
930
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
931
                                         errorTypeSo, None)
932
    else:
933
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
934
                                         intervalSo, errorTypeSo, None)
935
    if createdErrorType:
936
        sollya_lib_clear_obj(errorTypeSo)
937
    return(resultSo)
938
        
939

    
940
def pobyso_univar_polynomial_print_reverse(polySa):
941
    """ Legacy function. See pobyso_univar_polynomial_print_reverse_sa_sa. """
942
    return(pobyso_univar_polynomial_print_reverse_sa_sa(polySa))
943

    
944
def pobyso_univar_polynomial_print_reverse_sa_sa(polySa):
945
    """
946
    Return the string representation of a univariate polynomial with
947
    monomials ordered in the x^0..x^n order of the monomials.
948
    Remember: Sage
949
    """
950
    polynomialRing = polySa.base_ring()
951
    # A very expensive solution:
952
    # -create a fake multivariate polynomial field with only one variable,
953
    #   specifying a negative lexicographical order;
954
    mpolynomialRing = PolynomialRing(polynomialRing.base(), \
955
                                     polynomialRing.variable_name(), \
956
                                     1, order='neglex')
957
    # - convert the univariate argument polynomial into a multivariate
958
    #   version;
959
    p = mpolynomialRing(polySa)
960
    # - return the string representation of the converted form.
961
    # There is no simple str() method defined for p's class.
962
    return(p.__str__())
963
#
964
print pobyso_get_prec()  
965
pobyso_set_prec(165)
966
print pobyso_get_prec()  
967
a=100
968
print type(a)
969
id(a)
970
print "Max arity: ", pobyso_max_arity
971
print "Function tripleDouble (43) as a string: ", pobyso_function_type_as_string(43)
972
print "Function None (44) as a string: ", pobyso_function_type_as_string(44)
973
print "...Pobyso check done"