Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 259

Historique | Voir | Annoter | Télécharger (86,22 ko)

1
"""
2
@file pobyso.py
3
Actual functions to use in Sage
4
@author S.T. 
5
@date 2012-11-13
6

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

10
pobyso functions come in five flavors: 
11
- the _so_so (arguments and returned objects are pointers to Sollya objects, 
12
  includes the void function and the no arguments function that return a 
13
  pointer to a Sollya object);
14
- the _so_sa (argument are pointers to Sollya objects, returned objects are
15
  Sage/Python objects or, more generally, information is transfered from the
16
  Sollya world to Sage/Python world; e.g. functions without arguments that
17
  return a Sage/Python object);
18
- the _sa_so (arguments are Sage/Python objects, returned objects are 
19
  pointers to Sollya objects);
20
- the sa_sa (arguments and returned objects are all Sage/Python objects);
21
- a catch all flavor, without any suffix, (e. g. functions that have no 
22
  argument nor return value).
23

24
This classification is not always very strict. Conversion functions from Sollya
25
to Sage/Python are sometimes decorated with Sage/Python arguments to set
26
the precision. These functions remain in the so_sa category.
27

28
@note
29
Reported errors in Eclipse come from the calls to the Sollya library
30

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

    
89
pobyso_max_arity = 9
90

    
91
def pobyso_absolute_so_so():
92
    """
93
    Create an "absolute" Sollya object.
94
    """
95
    return(sollya_lib_absolute(None))
96

    
97
def pobyso_autoprint(arg):
98
    sollya_lib_autoprint(arg, None)
99

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

    
138
def pobyso_build_end_elliptic_list_so_so(*args):
139
    """
140
    From Sollya argument objects, create a Sollya end elliptic list.
141
    Elements of the list are "eaten" (should not be cleared individually, 
142
    are cleared when the list is cleared).
143
    """
144
    if len(args) == 0:
145
        ## Called with an empty list produced "error".
146
        return sollya_lib_build_end_elliptic_list(None)
147
    index = 0
148
    ## One can not append elements to an elliptic list, prepend only is 
149
    #  permitted.
150
    for argument in reversed(args):
151
        if index == 0:
152
            listSo = sollya_lib_build_end_elliptic_list(argument, None)
153
        else:
154
            listSo = sollya_lib_prepend(argument, listSo)
155
        index += 1
156
    return listSo
157
        
158
# End pobyso_build_end_elliptic_list_so_so
159

    
160
def pobyso_build_function_sub_so_so(exp1So, exp2So):
161
    return sollya_lib_build_function_sub(exp1So, exp2So)
162

    
163
def pobyso_build_list_of_ints_sa_so(*args):
164
    """
165
    Build a Sollya list from Sage integral arguments. 
166
    """
167
    if len(args) == 0:
168
        return pobyso_build_list_so_so()
169
    ## Make a Sage list of integral constants.
170
    intsList = []
171
    for intElem in args:
172
        intsList.append(pobyso_constant_from_int_sa_so(intElem))
173
    return pobyso_build_list_so_so(*intsList) 
174
    
175
def pobyso_build_list_so_so(*args):
176
    """
177
    Make a Sollya list out of Sollya objects passed as arguments.
178
    If one wants to call it with a list argument, should prepend a "*"
179
    before the list variable name.
180
    Elements of the list are "eaten" (should not be cleared individually, 
181
    are cleared when the list is cleared).
182
    """
183
    if len(args) == 0:
184
        ## Called with an empty list produced "error".
185
        return sollya_lib_build_list(None)
186
    index = 0
187
    ## Append the Sollya elements one by one.
188
    for elementSo in args:
189
        if index == 0:
190
            listSo = sollya_lib_build_list(elementSo, None)
191
        else:
192
            listSo = sollya_lib_append(listSo, elementSo)
193
        index += 1
194
    return listSo
195
# End pobyso_build list_so_so
196
    
197

    
198
def pobyso_change_var_in_function_so_so(funcSo, chvarExpSo):
199
    """
200
    Variable change in a function.
201
    """
202
    return(sollya_lib_evaluate(funcSo,chvarExpSo))
203
# End pobyso_change_var_in_function_so_so     
204

    
205
def pobyso_chebyshevform_so_so(functionSo, degreeSo, intervalSo):
206
    resultSo = sollya_lib_chebyshevform(functionSo, degreeSo, intervalSo)
207
    return(resultSo)
208
# End pobyso_chebyshevform_so_so.
209

    
210
def pobyso_clear_obj(objSo):
211
    """
212
    Free a Sollya object's memory.
213
    Very thin wrapper around sollya_lib_clear_obj().
214
    """
215
    sollya_lib_clear_obj(objSo)
216
# End pobyso_clear_obj. 
217

    
218
def pobyso_clear_taylorform_sa_so(taylorFormSaSo):
219
    """
220
    This method is necessary to correctly clean up the memory from Taylor forms.
221
    These are made of a Sollya object, a Sollya object list, a Sollya object.
222
    For no clearly understood reason, sollya_lib_clear_object_list crashed 
223
    when applied to the object list.
224
    Here, we decompose it into Sage list of Sollya objects references and we
225
     clear them one at a time. 
226
    """
227
    sollya_lib_clear_obj(taylorFormSaSo[0])
228
    (coefficientsErrorsListSaSo, numElementsSa, isEndEllipticSa) = \
229
        pobyso_get_list_elements_so_so(taylorFormSaSo[1])
230
    for element in coefficientsErrorsListSaSo:
231
        sollya_lib_clear_obj(element)
232
    sollya_lib_clear_obj(taylorFormSaSo[1])
233
    sollya_lib_clear_obj(taylorFormSaSo[2])
234
# End pobyso_clear_taylorform_sa_so 
235

    
236
def pobyso_cmp(rnArgSa, cteSo):
237
    """
238
    Deprecated, use pobyso_cmp_sa_so_sa instead.
239
    """
240
    print "Deprecated, use pobyso_cmp_sa_so_sa instead."
241
    return pobyso_cmp_sa_so_sa(rnArgSa, cteSo)
242
# End  pobyso_cmp
243
    
244
def pobyso_cmp_sa_so_sa(rnArgSa, cteSo):
245
    """
246
    Compare the MPFR value a RealNumber with that of a Sollya constant.
247
    
248
    Get the value of the Sollya constant into a RealNumber and compare
249
    using MPFR. Could be optimized by working directly with a mpfr_t
250
    for the intermediate number. 
251
    """
252
    # Get the precision of the Sollya constant to build a Sage RealNumber 
253
    # with enough precision.to hold it.
254
    precisionOfCte = c_int(0)
255
    # From the Sollya constant, create a local Sage RealNumber.
256
    sollya_lib_get_prec_of_constant(precisionOfCte, cteSo) 
257
    #print "Precision of constant: ", precisionOfCte
258
    RRRR = RealField(precisionOfCte.value)
259
    rnLocalSa = RRRR(0)
260
    sollya_lib_get_constant(get_rn_value(rnLocalSa), cteSo)
261
    #
262
    ## Compare the Sage RealNumber version of the Sollya constant with rnArg
263
    #  through a direct comparison of underlying MPFR numbers.
264
    return cmp_rn_value(rnArgSa, rnLocal)
265
# End pobyso_smp_sa_so_sa
266

    
267
def pobyso_compute_pos_function_abs_val_bounds_sa_sa(funcSa, lowerBoundSa, \
268
                                                     upperBoundSa):
269
    """
270
    TODO: completely rework and test.
271
    """
272
    pobyso = pobyso_name_free_variable_sa_so(funcSa.variables()[0])
273
    funcSo = pobyso_parse_string(funcSa._assume_str().replace('_SAGE_VAR_', ''))
274
    rangeSo = pobyso_range_sa_so(lowerBoundSa, upperBoundSa)
275
    infnormSo = pobyso_infnorm_so_so(funcSo,rangeSo)
276
    # Sollya return the infnorm as an interval.
277
    fMaxSa = pobyso_get_interval_from_range_so_sa(infnormSo)
278
    # Get the top bound and compute the binade top limit.
279
    fMaxUpperBoundSa = fMaxSa.upper()
280
    binadeTopLimitSa = 2**ceil(fMaxUpperBoundSa.log2())
281
    # Put up together the function to use to compute the lower bound.
282
    funcAuxSo = pobyso_parse_string(str(binadeTopLimitSa) +  \
283
                                    '-(' + f._assume_str().replace('_SAGE_VAR_', '') + ')')
284
    pobyso_autoprint(funcAuxSo)
285
    # Clear the Sollya range before a new call to infnorm and issue the call.
286
    sollya_lib_clear_obj(infnormSo)
287
    infnormSo = pobyso_infnorm_so_so(funcAuxSo,rangeSo)
288
    fMinSa = pobyso_get_interval_from_range_so_sa(infnormSo)
289
    sollya_lib_clear_obj(infnormSo)
290
    fMinLowerBoundSa = binadeTopLimitSa - fMinSa.lower()
291
    # Compute the maximum of the precisions of the different bounds.
292
    maxPrecSa = max([fMinLowerBoundSa.parent().precision(), \
293
                     fMaxUpperBoundSa.parent().precision()])
294
    # Create a RealIntervalField and create an interval with the "good" bounds.
295
    RRRI = RealIntervalField(maxPrecSa)
296
    imageIntervalSa = RRRI(fMinLowerBoundSa, fMaxUpperBoundSa)
297
    # Free the unneeded Sollya objects
298
    sollya_lib_clear_obj(funcSo)
299
    sollya_lib_clear_obj(funcAuxSo)
300
    sollya_lib_clear_obj(rangeSo)
301
    return(imageIntervalSa)
302
# End pobyso_compute_pos_function_abs_val_bounds_sa_sa
303

    
304
def pobyso_compute_precision_decay_ratio_function_sa_so():
305
    """
306
    Compute the precision decay ratio function for polynomial 
307
    coefficient progressive trucation.
308
    """
309
    functionText = """
310
    proc(deg, a, b, we, wq) 
311
    {
312
      k = we * (exp(x/a)-1) + wq * (b*x)^2 + (1-we-wq) * x;
313
      return k/k(d);
314
    };
315
    """
316
    return pobyso_parse_string_sa_so(functionText)
317
# End  pobyso_compute_precision_decay_ratio_function.
318

    
319

    
320
def pobyso_constant(rnArg):
321
    """ Legacy function. See pobyso_constant_sa_so. """
322
    return(pobyso_constant_sa_so(rnArg))
323
    
324
def pobyso_constant_sa_so(rnArgSa, precisionSa=None):
325
    """
326
    Create a Sollya constant from a Sage RealNumber.
327
    The sollya_lib_constant() function creates a constant
328
    with the same precision as the source.
329
    """
330
    ## Precision stuff. If one wants to change precisions,
331
    #  everything takes place in Sage. That only makes
332
    #  sense if one wants to reduce the precision.
333
    # TODO: revisit precision stuff with new technique.
334
    if not precisionSa is None:
335
        RRR = RealField(precisionSa)
336
        rnArgSa = RRR(rnArgSa)
337
    #print rnArgSa, rnArgSa.precision()
338
    # Sollya constant creation takes place here.
339
    return sollya_lib_constant(get_rn_value(rnArgSa))
340
# End pobyso_constant_sa_so
341
     
342
def pobyso_constant_0_sa_so():
343
    """
344
    Obvious.
345
    """
346
    return pobyso_constant_from_int_sa_so(0)
347

    
348
def pobyso_constant_1():
349
    """
350
    Obvious.
351
    Legacy function. See pobyso_constant_so_so. 
352
    """
353
    return pobyso_constant_1_sa_so()
354

    
355
def pobyso_constant_1_sa_so():
356
    """
357
    Obvious.
358
    """
359
    return(pobyso_constant_from_int_sa_so(1))
360

    
361
def pobyso_constant_from_int(anInt):
362
    """ Legacy function. See pobyso_constant_from_int_sa_so. """
363
    return pobyso_constant_from_int_sa_so(anInt)
364

    
365
def pobyso_constant_from_int_sa_so(anInt):
366
    """
367
    Get a Sollya constant from a Sage int.
368
    """
369
    return sollya_lib_constant_from_int64(long(anInt))
370

    
371
def pobyso_constant_from_int_so_sa(constSo):
372
    """
373
    Get a Sage int from a Sollya int constant.
374
    Usefull for precision or powers in polynomials.
375
    """
376
    constSa = c_long(0)
377
    sollya_lib_get_constant_as_int64(byref(constSa), constSo)
378
    return constSa.value
379
# End pobyso_constant_from_int_so_sa
380

    
381
def pobyso_constant_from_mpq_sa_so(rationalSa):
382
    """
383
    Make a Sollya constant from Sage rational.
384
    The Sollya constant is an unevaluated expression.
385
    Hence no precision argument is needed.
386
    It is better to leave this way since Sollya has its own
387
    optimized evaluation mecanism that tries very hard to
388
    return exact values or at least faithful ones.
389
    """
390
    ratExprSo = \
391
        sollya_lib_constant_from_mpq(sgmp_get_rational_value(rationalSa))
392
    return ratExprSo
393
# End pobyso_constant_from_mpq_sa_so.
394

    
395
def pobyso_constant_sollya_prec_sa_so(rnArgSa):
396
    """
397
    Create a Sollya constant from a Sage RealNumber at the
398
    current precision in Sollya.
399
    """
400
    currentSollyaPrecSa = pobyso_get_prec_so_sa()
401
    return pobyso_constant_sa_so(rnArgSa, currentSollyaPrecSa)
402
# End pobyso_constant_sollya_prec_sa_so
403

    
404
def pobyso_end_elliptic_list_so_sa_so(objectsListSo, intCountSa):
405
    """
406
    Create a Sollya end elliptic list made of the objectListSo[0] to
407
     objectsListSo[intCountSa-1] objects.
408
    """     
409
    return sollya_lib_end_elliptic_list(objectSo, int(intCountSa))
410

    
411
def pobyso_error_so():
412
    return sollya_lib_error(None)
413
# End pobyso_error().
414

    
415
def pobyso_evaluate_so_so(funcSo, argumentSo):
416
    """
417
    Evaluates funcSo for arguemntSo through sollya_lib_evaluate().
418
    """
419
    return sollya_lib_evaluate(funcSo, argumentSo)
420
# End pobyso_evaluate_so_so.
421
#
422
def pobyso_diff_so_so(funcSo):
423
    """
424
    Very thin wrapper around sollya_lib_diff.
425
    """
426
    ## TODO: add a check to make sure funcSo is a functional expression.
427
    return sollya_lib_diff(funcSo)
428

    
429
def pobyso_dirty_find_zeros_so_so(funcSo, rangeSo):
430
    """
431
    Thin wrapper over sollya_lib_dirtyfindzeros()
432
    """
433
    return sollya_lib_dirtyfindzeros(funcSo, rangeSo)
434
# End pobys_dirty_find_zeros
435

    
436
def pobyso_dirty_inf_norm_so_so(funcSo, rangeSo, preSo=None):
437
    """
438
    Thin wrapper around sollya_dirtyinfnorm().
439
    """
440
    # TODO: manage the precision change.
441
    
442
    return sollya_lib_dirtyinfnorm(funcSo, rangeSo)
443
# End pobyso_dirty_inf_norm_so_so
444

    
445
def pobyso_float_list_so_sa(listSo):
446
    """
447
    Return a Sollya list of floating-point numbers as a Sage list of 
448
    floating-point numbers.
449
    TODO: add a test to make sure that each element of the list is a constant.
450
    """
451
    listSa   = []
452
    ## The function returns none if the list is empty or an error has happened.
453
    retVal = pobyso_get_list_elements_so_so(listSo)
454
    if retVal is None:
455
        return listSa
456
    ## Just in case the interface is changed and an empty list is returned 
457
    #  instead of None.
458
    elif len(retVal) == 0:
459
        return listSa
460
    else:
461
        ## Remember pobyso_get_list_elements_so_so returns more information
462
        #  than just the elements of the list (# elements, is_elliptic)
463
        listSaSo, numElements, isEndElliptic = retVal
464
    ## Return an empty list.
465
    if numElements == 0:
466
        return listSa
467
    ## Search first for the maximum precision of the elements
468
    maxPrecSa = 0
469
    for floatSo in listSaSo:
470
        #pobyso_autoprint(floatSo)
471
        curPrecSa =  pobyso_get_prec_of_constant_so_sa(floatSo)
472
        if curPrecSa > maxPrecSa:
473
            maxPrecSa = curPrecSa
474
    ##
475
    RF = RealField(maxPrecSa)
476
    ##
477
    for floatSo in listSaSo:
478
        listSa.append(pobyso_get_constant_as_rn_with_rf_so_sa(floatSo))
479
    return listSa
480
# End pobyso_float_list_so_sa
481

    
482
def pobyso_float_poly_sa_so(polySa, precSa = None):
483
    """
484
    Create a Sollya polynomial from a Sage RealField polynomial.
485
    """
486
    ## TODO: filter arguments.
487
    ## Precision. If a precision is given, convert the polynomial
488
    #  into the right polynomial field. If not convert it straight
489
    #  to Sollya.
490
    sollyaPrecChanged = False 
491
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
492
    if precSa is None:
493
        precSa = polySa.parent().base_ring().precision()
494
    if (precSa > initialSollyaPrecSa):
495
        assert precSa >= 2, "Precision change <2 requested"
496
        if precSa <= 2:
497
            print inspect.stack()[0][3], ": precision change <= 2 requested"
498
        precSo = pobyso_constant_from_int(precSa)
499
        pobyso_set_prec_so_so(precSo)
500
        sollya_lib_clear_obj(precSo)
501
        sollyaPrecChanged = True    
502
    ## Get exponents and coefficients.
503
    exponentsSa     = polySa.exponents()
504
    coefficientsSa  = polySa.coefficients()
505
    ## Build the polynomial.
506
    polySo = None
507
    for coefficientSa, exponentSa in zip(coefficientsSa, exponentsSa):
508
        #print coefficientSa.n(prec=precSa), exponentSa
509
        coefficientSo = \
510
            pobyso_constant_sa_so(coefficientSa)
511
        #pobyso_autoprint(coefficientSo)
512
        exponentSo = \
513
            pobyso_constant_from_int_sa_so(exponentSa)
514
        #pobyso_autoprint(exponentSo)
515
        monomialSo = sollya_lib_build_function_pow(
516
                       sollya_lib_build_function_free_variable(),
517
                       exponentSo)
518
        polyTermSo = sollya_lib_build_function_mul(coefficientSo,
519
                                                       monomialSo)
520
        if polySo is None:
521
            polySo = polyTermSo
522
        else:
523
            polySo = sollya_lib_build_function_add(polySo, polyTermSo)
524
    if sollyaPrecChanged:
525
        pobyso_set_prec_so_so(initialSollyaPrecSo)
526
        sollya_lib_clear_obj(initialSollyaPrecSo)
527
    return polySo
528
# End pobyso_float_poly_sa_so    
529

    
530
def pobyso_float_poly_so_sa(polySo, realFieldSa=None):
531
    """
532
    Convert a Sollya polynomial into a Sage floating-point polynomial.
533
    If no realField is given, a RealField corresponding to the maximum 
534
    precision of the coefficients is internally computed.
535
    The real field is not returned but can be easily retrieved from 
536
    the polynomial itself.
537
    ALGORITHM:
538
    - (optional) compute the RealField of the coefficients;
539
    - convert the Sollya expression into a Sage expression;
540
    - convert the Sage expression into a Sage polynomial
541
    """    
542
    if realFieldSa is None:
543
        expressionPrecSa = pobyso_get_max_prec_of_exp_so_sa(polySo)
544
        #print "Maximum precision of Sollya polynomial coefficients:", expressionPrecSa
545
        if expressionPrecSa < 2 or expressionPrecSa > 2147483391:
546
            print "Maximum degree of expression:", expressionPrecSa
547
        realFieldSa      = RealField(expressionPrecSa)
548
    #print "Sollya expression before...",
549
    #pobyso_autoprint(polySo)
550

    
551
    expressionSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo,
552
                                                             realFieldSa)
553
    #print "...Sollya expression after."
554
    #pobyso_autoprint(polySo)
555
    polyVariableSa = expressionSa.variables()[0]
556
    polyRingSa     = realFieldSa[str(polyVariableSa)]
557
    #print polyRingSa
558
    # Do not use the polynomial(expressionSa, ring=polyRingSa) form!
559
    polynomialSa = polyRingSa(expressionSa)
560
    polyCoeffsListSa = polynomialSa.coefficients()
561
    #for coeff in polyCoeffsListSa:
562
    #    print coeff.abs().n()
563
    return polynomialSa
564
# End pobyso_float_poly_so_sa
565

    
566
def pobyso_free_variable():
567
    """
568
    Ultra thin wrapper around the sollya_lib_function_build_free_variable function.
569
    """
570
    return sollya_lib_build_function_free_variable()
571
    
572
def pobyso_function_type_as_string(funcType):
573
    """ Legacy function. See pobyso_function_type_as_string_so_sa. """
574
    return(pobyso_function_type_as_string_so_sa(funcType))
575

    
576
def pobyso_function_type_as_string_so_sa(funcType):
577
    """
578
    Numeric Sollya function codes -> Sage mathematical function names.
579
    Notice that pow -> ^ (a la Sage, not a la Python).
580
    """
581
    if funcType == SOLLYA_BASE_FUNC_ABS:
582
        return "abs"
583
    elif funcType == SOLLYA_BASE_FUNC_ACOS:
584
        return "arccos"
585
    elif funcType == SOLLYA_BASE_FUNC_ACOSH:
586
        return "arccosh"
587
    elif funcType == SOLLYA_BASE_FUNC_ADD:
588
        return "+"
589
    elif funcType == SOLLYA_BASE_FUNC_ASIN:
590
        return "arcsin"
591
    elif funcType == SOLLYA_BASE_FUNC_ASINH:
592
        return "arcsinh"
593
    elif funcType == SOLLYA_BASE_FUNC_ATAN:
594
        return "arctan"
595
    elif funcType == SOLLYA_BASE_FUNC_ATANH:
596
        return "arctanh"
597
    elif funcType == SOLLYA_BASE_FUNC_CEIL:
598
        return "ceil"
599
    elif funcType == SOLLYA_BASE_FUNC_CONSTANT:
600
        return "cte"
601
    elif funcType == SOLLYA_BASE_FUNC_COS:
602
        return "cos"
603
    elif funcType == SOLLYA_BASE_FUNC_COSH:
604
        return "cosh"
605
    elif funcType == SOLLYA_BASE_FUNC_DIV:
606
        return "/"
607
    elif funcType == SOLLYA_BASE_FUNC_DOUBLE:
608
        return "double"
609
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEDOUBLE:
610
        return "doubleDouble"
611
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEEXTENDED:
612
        return "doubleDxtended"
613
    elif funcType == SOLLYA_BASE_FUNC_ERF:
614
        return "erf"
615
    elif funcType == SOLLYA_BASE_FUNC_ERFC:
616
        return "erfc"
617
    elif funcType == SOLLYA_BASE_FUNC_EXP:
618
        return "exp"
619
    elif funcType == SOLLYA_BASE_FUNC_EXP_M1:
620
        return "expm1"
621
    elif funcType == SOLLYA_BASE_FUNC_FLOOR:
622
        return "floor"
623
    elif funcType == SOLLYA_BASE_FUNC_FREE_VARIABLE:
624
        return "freeVariable"
625
    elif funcType == SOLLYA_BASE_FUNC_HALFPRECISION:
626
        return "halfPrecision"
627
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYCONSTANT:
628
        return "libraryConstant"
629
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYFUNCTION:
630
        return "libraryFunction"
631
    elif funcType == SOLLYA_BASE_FUNC_LOG:
632
        return "log"
633
    elif funcType == SOLLYA_BASE_FUNC_LOG_10:
634
        return "log10"
635
    elif funcType == SOLLYA_BASE_FUNC_LOG_1P:
636
        return "log1p"
637
    elif funcType == SOLLYA_BASE_FUNC_LOG_2:
638
        return "log2"
639
    elif funcType == SOLLYA_BASE_FUNC_MUL:
640
        return "*"
641
    elif funcType == SOLLYA_BASE_FUNC_NEARESTINT:
642
        return "round"
643
    elif funcType == SOLLYA_BASE_FUNC_NEG:
644
        return "__neg__"
645
    elif funcType == SOLLYA_BASE_FUNC_PI:
646
        return "pi"
647
    elif funcType == SOLLYA_BASE_FUNC_POW:
648
        return "^"
649
    elif funcType == SOLLYA_BASE_FUNC_PROCEDUREFUNCTION:
650
        return "procedureFunction"
651
    elif funcType == SOLLYA_BASE_FUNC_QUAD:
652
        return "quad"
653
    elif funcType == SOLLYA_BASE_FUNC_SIN:
654
        return "sin"
655
    elif funcType == SOLLYA_BASE_FUNC_SINGLE:
656
        return "single"
657
    elif funcType == SOLLYA_BASE_FUNC_SINH:
658
        return "sinh"
659
    elif funcType == SOLLYA_BASE_FUNC_SQRT:
660
        return "sqrt"
661
    elif funcType == SOLLYA_BASE_FUNC_SUB:
662
        return "-"
663
    elif funcType == SOLLYA_BASE_FUNC_TAN:
664
        return "tan"
665
    elif funcType == SOLLYA_BASE_FUNC_TANH:
666
        return "tanh"
667
    elif funcType == SOLLYA_BASE_FUNC_TRIPLEDOUBLE:
668
        return "tripleDouble"
669
    else:
670
        return None
671

    
672
def pobyso_get_constant(rnArgSa, constSo):
673
    """ Legacy function. See pobyso_get_constant_so_sa. """
674
    return pobyso_get_constant_so_sa(rnArgSa, constSo)
675
# End pobyso_get_constant
676

    
677
def pobyso_get_constant_so_sa(rnArgSa, constSo):
678
    """
679
    Set the value of rnArgSo to the value of constSo in MPFR_RNDN mode.
680
    rnArg must already exist and belong to some RealField.
681
    We assume that constSo points to a Sollya constant.
682
    """
683
    outcome = sollya_lib_get_constant(get_rn_value(rnArgSa), constSo)
684
    if outcome == 0: # Failure because constSo is not a constant expression.
685
        return None
686
    else:
687
        return outcome
688
# End  pobyso_get_constant_so_sa
689
   
690
def pobyso_get_constant_as_rn(ctExpSo):
691
    """ 
692
    Legacy function. See pobyso_get_constant_as_rn_so_sa. 
693
    """ 
694
    return(pobyso_get_constant_as_rn_so_sa(ctExpSo))
695
    
696
def pobyso_get_constant_as_rn_so_sa(constExpSo):
697
    """
698
    Get a Sollya constant as a Sage "real number".
699
    The precision of the floating-point number returned is that of the Sollya
700
    constant.
701
    """
702
    #print "Before computing precision of variable..."
703
    #pobyso_autoprint(constExpSo)
704
    precisionSa  = pobyso_get_prec_of_constant_so_sa(constExpSo)
705
    #print "precisionSa:", precisionSa
706
    ## If the expression can not be exactly converted, None is returned.
707
    #  In this case opt for the Sollya current expression.
708
    if precisionSa is None:
709
        precisionSa = pobyso_get_prec_so_sa()
710
    RRRR = RealField(precisionSa)
711
    rnSa = RRRR(0)
712
    outcome = sollya_lib_get_constant(get_rn_value(rnSa), constExpSo)
713
    if outcome == 0:
714
        return None
715
    else:
716
        return rnSa
717
# End pobyso_get_constant_as_rn_so_sa
718

    
719
def pobyso_get_constant_as_rn_with_rf(ctExp, realField):
720
    """ 
721
    Legacy function. See pobyso_get_constant_as_rn_with_rf_so_sa.
722
    """
723
    return pobyso_get_constant_as_rn_with_rf_so_sa(ctExp, realField)
724
# End pobyso_get_constant_as_rn_with_rf
725
    
726
def pobyso_get_constant_as_rn_with_rf_so_sa(ctExpSo, realFieldSa = None):
727
    """
728
    Get a Sollya constant as a Sage "real number".
729
    If no real field is specified, the precision of the floating-point number 
730
    returned is that of the Sollya constant.
731
    Otherwise is is that of the real field. Hence rounding may happen.
732
    """
733
    if realFieldSa is None:
734
        return pobyso_get_constant_as_rn_so_sa(ctExpSo)
735
    rnSa = realFieldSa(0)
736
    outcome = sollya_lib_get_constant(get_rn_value(rnSa), ctExpSo)
737
    if outcome == 0:
738
        return None
739
    else:
740
        return rnSa
741
# End pobyso_get_constant_as_rn_with_rf_so_sa
742

    
743
def pobyso_get_free_variable_name():
744
    """ 
745
    Legacy function. See pobyso_get_free_variable_name_so_sa.
746
    """
747
    return(pobyso_get_free_variable_name_so_sa())
748

    
749
def pobyso_get_free_variable_name_so_sa():
750
    return sollya_lib_get_free_variable_name()
751
    
752
def pobyso_get_function_arity(expressionSo):
753
    """ 
754
    Legacy function. See pobyso_get_function_arity_so_sa.
755
    """
756
    return(pobyso_get_function_arity_so_sa(expressionSo))
757

    
758
def pobyso_get_function_arity_so_sa(expressionSo):
759
    arity = c_int(0)
760
    sollya_lib_get_function_arity(byref(arity),expressionSo)
761
    return int(arity.value)
762

    
763
def pobyso_get_head_function(expressionSo):
764
    """ 
765
    Legacy function. See pobyso_get_head_function_so_sa. 
766
    """
767
    return(pobyso_get_head_function_so_sa(expressionSo)) 
768

    
769
def pobyso_get_head_function_so_sa(expressionSo):
770
    functionType = c_int(0)
771
    sollya_lib_get_head_function(byref(functionType), expressionSo)
772
    return int(functionType.value)
773

    
774
def pobyso_get_interval_from_range_so_sa(soRange, realIntervalFieldSa = None ):
775
    """
776
    Return the Sage interval corresponding to the Sollya range argument.
777
    If no reaIntervalField is passed as an argument, the interval bounds are not
778
    rounded: they are elements of RealIntervalField of the "right" precision
779
    to hold all the digits.
780
    """
781
    prec = c_int(0)
782
    if realIntervalFieldSa is None:
783
        retval = sollya_lib_get_prec_of_range(byref(prec), soRange, None)
784
        if retval == 0:
785
            return None
786
        realIntervalFieldSa = RealIntervalField(prec.value)
787
    intervalSa = realIntervalFieldSa(0,0)
788
    retval = \
789
        sollya_lib_get_interval_from_range(get_interval_value(intervalSa),\
790
                                           soRange)
791
    if retval == 0:
792
        return None
793
    return intervalSa
794
# End pobyso_get_interval_from_range_so_sa
795

    
796
def pobyso_get_list_elements(soObj):
797
    """ Legacy function. See pobyso_get_list_elements_so_so. """
798
    return pobyso_get_list_elements_so_so(soObj)
799
 
800
def pobyso_get_list_elements_so_so(objectListSo):
801
    """
802
    Get the Sollya list elements as a Sage/Python array of Sollya objects.
803
    
804
    INPUT:
805
    - objectListSo: a Sollya list of Sollya objects.
806
    
807
    OUTPUT:
808
    - a Sage/Python tuple made of:
809
      - a Sage/Python list of Sollya objects,
810
      - a Sage/Python int holding the number of elements,
811
      - a Sage/Python int stating (!= 0) that the list is end-elliptic.
812
    NOTE::
813
        We recover the addresses of the Sollya object from the list of pointers
814
        returned by sollya_lib_get_list_elements. The list itself is freed.
815
    TODO::
816
        Figure out what to do with numElements since the number of elements
817
        can easily be recovered from the list itself. 
818
        Ditto for isEndElliptic.
819
    """
820
    listAddress = POINTER(c_longlong)()
821
    numElements = c_int(0)
822
    isEndElliptic = c_int(0)
823
    listAsSageList = []
824
    result = sollya_lib_get_list_elements(byref(listAddress),\
825
                                          byref(numElements),\
826
                                          byref(isEndElliptic),\
827
                                          objectListSo)
828
    if result == 0 :
829
        return None
830
    for i in xrange(0, numElements.value, 1):
831
       #listAsSageList.append(sollya_lib_copy_obj(listAddress[i]))
832
       listAsSageList.append(listAddress[i])
833
       # Clear each of the elements returned by Sollya.
834
       #sollya_lib_clear_obj(listAddress[i])
835
    # Free the list itself.   
836
    sollya_lib_free(listAddress)
837
    return (listAsSageList, numElements.value, isEndElliptic.value)
838

    
839
def pobyso_get_max_prec_of_exp(soExp):
840
    """ Legacy function. See pobyso_get_max_prec_of_exp_so_sa. """
841
    return pobyso_get_max_prec_of_exp_so_sa(soExp)
842

    
843
def pobyso_get_max_prec_of_exp_so_sa(expSo):
844
    """
845
    Get the maximum precision used for the numbers in a Sollya expression.
846
    
847
    Arguments:
848
    soExp -- a Sollya expression pointer
849
    Return value:
850
    A Python integer
851
    TODO: 
852
    - error management;
853
    - correctly deal with numerical type such as DOUBLEEXTENDED.
854
    """
855
    if expSo is None:
856
        print inspect.stack()[0][3], ": expSo is None."
857
        return 0
858
    maxPrecision = 0
859
    minConstPrec = 0
860
    currentConstPrec = 0
861
    #pobyso_autoprint(expSo)
862
    operator = pobyso_get_head_function_so_sa(expSo)
863
    if (operator != SOLLYA_BASE_FUNC_CONSTANT) and \
864
    (operator != SOLLYA_BASE_FUNC_FREE_VARIABLE):
865
        (arity, subexpressions) = pobyso_get_subfunctions_so_sa(expSo)
866
        for i in xrange(arity):
867
            maxPrecisionCandidate = \
868
                pobyso_get_max_prec_of_exp_so_sa(subexpressions[i])
869
            if maxPrecisionCandidate > maxPrecision:
870
                maxPrecision = maxPrecisionCandidate
871
        return maxPrecision
872
    elif operator == SOLLYA_BASE_FUNC_CONSTANT:
873
        #minConstPrec = pobyso_get_min_prec_of_constant_so_sa(expSo)
874
        #currentConstPrec = pobyso_get_min_prec_of_constant_so_sa(soExp)
875
        #print minConstPrec, " - ", currentConstPrec 
876
        return pobyso_get_min_prec_of_constant_so_sa(expSo)
877
    
878
    elif operator == SOLLYA_BASE_FUNC_FREE_VARIABLE:
879
        return 0
880
    else:
881
        print "pobyso_get_max_prec_of_exp_so_sa: unexepected operator."
882
        return 0
883

    
884
def pobyso_get_min_prec_of_constant_so_sa(constExpSo):
885
    """
886
    Get the minimum precision necessary to represent the value of a Sollya
887
    constant.
888
    MPFR_MIN_PREC and powers of 2 are taken into account.
889
    We assume that constExpSo is a pointer to a Sollay constant expression.
890
    """
891
    constExpAsRnSa = pobyso_get_constant_as_rn_so_sa(constExpSo)
892
    return(min_mpfr_size(get_rn_value(constExpAsRnSa)))
893

    
894
def pobyso_get_poly_so_sa(polySo, realFieldSa=None):
895
    """
896
    Convert a Sollya polynomial into a Sage polynomial.
897
    Legacy function. Use pobyso_float_poly_so_sa() instead.
898
    """    
899
    return pobyso_float_poly_so_sa(polySo,realFieldSa)
900
# End pobyso_get_poly_so_sa
901

    
902
def pobyso_get_prec():
903
    """ Legacy function. See pobyso_get_prec_so_sa(). """
904
    return pobyso_get_prec_so_sa()
905

    
906
def pobyso_get_prec_so():
907
    """
908
    Get the current default precision in Sollya.
909
    The return value is a Sollya object.
910
    Usefull when modifying the precision back and forth by avoiding
911
    extra conversions.
912
    """
913
    return sollya_lib_get_prec(None)
914
    
915
def pobyso_get_prec_so_sa():
916
    """
917
    Get the current default precision in Sollya.
918
    The return value is Sage/Python int.
919
    """
920
    precSo = sollya_lib_get_prec()
921
    precSa = pobyso_constant_from_int_so_sa(precSo)
922
    sollya_lib_clear_obj(precSo)
923
    return precSa
924
# End pobyso_get_prec_so_sa.
925

    
926
def pobyso_get_prec_so_so_sa():
927
    """
928
    Return the current precision both as a Sollya object and a
929
    Sage integer as hybrid tuple.
930
    To avoid multiple calls for precision manipulations. 
931
    """
932
    precSo = sollya_lib_get_prec()
933
    precSa = pobyso_constant_from_int_so_sa(precSo)
934
    return (precSo, int(precSa))
935
    
936
def pobyso_get_prec_of_constant(ctExpSo):
937
    """ Legacy function. See pobyso_get_prec_of_constant_so_sa. """
938
    return pobyso_get_prec_of_constant_so_sa(ctExpSo)
939

    
940
def pobyso_get_prec_of_constant_so_sa(ctExpSo):
941
    """
942
    Tries to find a precision to represent ctExpSo without rounding.
943
    If not possible, returns None.
944
    """
945
    #print "Entering pobyso_get_prec_of_constant_so_sa..."
946
    prec = c_int(0)
947
    retc = sollya_lib_get_prec_of_constant(byref(prec), ctExpSo, None)
948
    if retc == 0:
949
        #print "pobyso_get_prec_of_constant_so_sa failed."
950
        return None
951
    #print "...exiting pobyso_get_prec_of_constant_so_sa."
952
    return int(prec.value)
953

    
954
def pobyso_get_prec_of_range_so_sa(rangeSo):
955
    """
956
    Returns the number of bits elements of a range are coded with.
957
    """
958
    prec = c_int(0)
959
    retc = sollya_lib_get_prec_of_range(byref(prec), rangeSo, None)
960
    if retc == 0:
961
        return(None)
962
    return int(prec.value)
963
# End pobyso_get_prec_of_range_so_sa()
964

    
965
def pobyso_get_sage_exp_from_sollya_exp(sollyaExpSo, realField = RR):
966
    """ Legacy function. See pobyso_get_sage_exp_from_sollya_exp_so_sa. """
967
    return pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, 
968
                                                     realField = RR)
969

    
970
def pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, realFieldSa = RR):
971
    """
972
    Get a Sage expression from a Sollya expression. 
973
    Currently only tested with polynomials with floating-point coefficients.
974
    Notice that, in the returned polynomial, the exponents are RealNumbers.
975
    """
976
    #pobyso_autoprint(sollyaExp)
977
    operatorSa = pobyso_get_head_function_so_sa(sollyaExpSo)
978
    sollyaLibFreeVariableName = sollya_lib_get_free_variable_name()
979
    ## Get rid of the "_"'s in "_x_", if any.
980
    sollyaLibFreeVariableName = re.sub('_', '', sollyaLibFreeVariableName)
981
    # Constants and the free variable are special cases.
982
    # All other operator are dealt with in the same way.
983
    if (operatorSa != SOLLYA_BASE_FUNC_CONSTANT) and \
984
       (operatorSa != SOLLYA_BASE_FUNC_FREE_VARIABLE):
985
        (aritySa, subexpressionsSa) = pobyso_get_subfunctions_so_sa(sollyaExpSo)
986
        if aritySa == 1:
987
            sageExpSa = eval(pobyso_function_type_as_string_so_sa(operatorSa) + \
988
            "(" + pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], \
989
            realFieldSa) + ")")
990
        elif aritySa == 2:
991
            # We do not get through the preprocessor.
992
            # The "^" operator is then a special case.
993
            if operatorSa == SOLLYA_BASE_FUNC_POW:
994
                operatorAsStringSa = "**"
995
            else:
996
                operatorAsStringSa = \
997
                    pobyso_function_type_as_string_so_sa(operatorSa)
998
            sageExpSa = \
999
              eval("pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], realFieldSa)"\
1000
              + " " + operatorAsStringSa + " " + \
1001
                   "pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[1], realFieldSa)")
1002
        # We do not know yet how to deal with arity >= 3 
1003
        # (is there any in Sollya anyway?).
1004
        else:
1005
            sageExpSa = eval('None')
1006
        return sageExpSa
1007
    elif operatorSa == SOLLYA_BASE_FUNC_CONSTANT:
1008
        #print "This is a constant"
1009
        return pobyso_get_constant_as_rn_with_rf_so_sa(sollyaExpSo, realFieldSa)
1010
    elif operatorSa == SOLLYA_BASE_FUNC_FREE_VARIABLE:
1011
        #print "This is the free variable"
1012
        return eval(sollyaLibFreeVariableName)
1013
    else:
1014
        print "Unexpected"
1015
        return eval('None')
1016
# End pobyso_get_sage_exp_from_sollya_exp_so_sa
1017

    
1018

    
1019
def pobyso_get_subfunctions(expressionSo):
1020
    """ Legacy function. See pobyso_get_subfunctions_so_sa. """
1021
    return pobyso_get_subfunctions_so_sa(expressionSo) 
1022
# End pobyso_get_subfunctions.
1023
 
1024
def pobyso_get_subfunctions_so_sa(expressionSo):
1025
    """
1026
    Get the subfunctions of an expression.
1027
    Return the number of subfunctions and the list of subfunctions addresses.
1028
    S.T.: Could not figure out another way than that ugly list of declarations
1029
    to recover the addresses of the subfunctions.
1030
    We limit ourselves to arity 8 functions. 
1031
    """
1032
    subf0 = c_int(0)
1033
    subf1 = c_int(0)
1034
    subf2 = c_int(0)
1035
    subf3 = c_int(0)
1036
    subf4 = c_int(0)
1037
    subf5 = c_int(0)
1038
    subf6 = c_int(0)
1039
    subf7 = c_int(0)
1040
    subf8 = c_int(0)
1041
    arity = c_int(0)
1042
    nullPtr = POINTER(c_int)()
1043
    sollya_lib_get_subfunctions(expressionSo, byref(arity), \
1044
      byref(subf0), byref(subf1), byref(subf2), byref(subf3), \
1045
      byref(subf4), byref(subf5),\
1046
      byref(subf6), byref(subf7), byref(subf8), nullPtr, None) 
1047
#    byref(cast(subfunctions[0], POINTER(c_int))), \
1048
#    byref(cast(subfunctions[0], POINTER(c_int))), \
1049
#    byref(cast(subfunctions[2], POINTER(c_int))), \
1050
#    byref(cast(subfunctions[3], POINTER(c_int))), \
1051
#    byref(cast(subfunctions[4], POINTER(c_int))), \
1052
#    byref(cast(subfunctions[5], POINTER(c_int))), \
1053
#    byref(cast(subfunctions[6], POINTER(c_int))), \
1054
#    byref(cast(subfunctions[7], POINTER(c_int))), \
1055
#    byref(cast(subfunctions[8], POINTER(c_int))), nullPtr)
1056
    subfunctions = [subf0, subf1, subf2, subf3, subf4, subf5, subf6, subf7, \
1057
                    subf8]
1058
    subs = []
1059
    if arity.value > pobyso_max_arity:
1060
        return(0,[])
1061
    for i in xrange(arity.value):
1062
        subs.append(int(subfunctions[i].value))
1063
        #print subs[i]
1064
    return (int(arity.value), subs)
1065
# End pobyso_get_subfunctions_so_sa
1066
    
1067
def pobyso_guess_degree_sa_sa(functionSa, intervalSa, approxErrorSa, 
1068
                              weightSa=None, degreeBoundSa=None):
1069
    """
1070
    Sa_sa variant of the solly_guessdegree function.
1071
    Return 0 if something goes wrong.
1072
    """
1073
    functionAsStringSa = functionSa._assume_str().replace('_SAGE_VAR_', '')
1074
    functionSo = pobyso_parse_string_sa_so(functionAsStringSa)
1075
    if pobyso_is_error_so_sa(functionSo):
1076
        sollya_lib_clear_obj(functionSo)
1077
        return 0
1078
    rangeSo = pobyso_interval_to_range_sa_so(intervalSa)
1079
    # The approximation error is expected to be a floating point number.
1080
    if pobyso_is_floating_point_number_sa_sa(approxErrorSa):
1081
        approxErrorSo = pobyso_constant_sa_so(approxErrorSa)
1082
    else:
1083
        approxErrorSo = pobyso_constant_sa_so(RR(approxErrorSa))
1084
    if not weightSa is None:
1085
        weightAsStringSa = weightSa._assume_str().replace('_SAGE_VAR_', '')
1086
        weightSo = pobyso_parse_string_sa_so(weightAsStringSa)
1087
        if pobyso_is_error_so_sa(weightSo):
1088
            sollya_lib_clear_obj(functionSo)
1089
            sollya_lib_clear_obj(rangeSo)
1090
            sollya_lib_clear_obj(approxErrorSo)   
1091
            sollya_lib_clear_obj(weightSo)
1092
            return 0   
1093
    else:
1094
        weightSo = None
1095
    if not degreeBoundSa is None:
1096
        degreeBoundSo = pobyso_constant_from_int_sa_so(degreeBoundSa)
1097
    else:
1098
        degreeBoundSo = None
1099
    guessedDegreeSa = pobyso_guess_degree_so_sa(functionSo,
1100
                                                rangeSo,
1101
                                                approxErrorSo,
1102
                                                weightSo,
1103
                                                degreeBoundSo)
1104
    sollya_lib_clear_obj(functionSo)
1105
    sollya_lib_clear_obj(rangeSo)
1106
    sollya_lib_clear_obj(approxErrorSo)
1107
    if not weightSo is None:
1108
        sollya_lib_clear_obj(weightSo)
1109
    if not degreeBoundSo is None:
1110
        sollya_lib_clear_obj(degreeBoundSo)
1111
    return guessedDegreeSa
1112
# End poyso_guess_degree_sa_sa
1113

    
1114
def pobyso_guess_degree_so_sa(functionSo, rangeSo, errorSo, weightSo=None, \
1115
                              degreeBoundSo=None):
1116
    """
1117
    Thin wrapper around the guessdegree function.
1118
    Nevertheless, some precision control stuff has been appended.
1119
    """
1120
    # Deal with Sollya internal precision issues: if it is too small, 
1121
    # compared with the error, increases it to about twice -log2(error).
1122
    errorSa = pobyso_get_constant_as_rn_with_rf_so_sa(errorSo)
1123
    log2ErrorSa = errorSa.log2()
1124
    if log2ErrorSa < 0:
1125
        neededPrecisionSa = int(2 * int(-log2ErrorSa) / 64) * 64
1126
    else:
1127
        neededPrecisionSa = int(2 * int(log2ErrorSa) / 64) * 64
1128
    #print "Needed precision:", neededPrecisionSa
1129
    sollyaPrecisionHasChanged = False
1130
    (initialPrecSo, initialPrecSa) = pobyso_get_prec_so_so_sa()
1131
    if neededPrecisionSa > initialPrecSa:
1132
        if neededPrecisionSa <= 2:
1133
            print inspect.stack()[0][3], ": precision change <= 2 requested."
1134
        pobyso_set_prec_sa_so(neededPrecisionSa)
1135
        sollyaPrecisionHasChanged = True
1136
    #print "Guessing degree..."
1137
    # weightSo and degreeBoundsSo are optional arguments.
1138
    # As declared, sollya_lib_guessdegree must take 5 arguments.
1139
    if weightSo is None:
1140
        degreeRangeSo = sollya_lib_guessdegree(functionSo, rangeSo, errorSo, 
1141
                                               0, 0, None)
1142
    elif degreeBoundSo is None:
1143
        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, \
1144
                                                errorSo, weightSo, 0, None)
1145
    else:
1146
        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, errorSo, \
1147
                                                weightSo, degreeBoundSo, None)
1148
    #print "...degree guess done."
1149
    # Restore internal precision, if applicable.
1150
    if sollyaPrecisionHasChanged:
1151
        pobyso_set_prec_so_so(initialPrecSo)
1152
        sollya_lib_clear_obj(initialPrecSo)
1153
    degreeIntervalSa = pobyso_range_to_interval_so_sa(degreeRangeSo)
1154
    sollya_lib_clear_obj(degreeRangeSo)
1155
    # When ok, both bounds match.
1156
    # When the degree bound is too low, the upper bound is the degree
1157
    # for which the error can be honored.
1158
    # When it really goes wrong, the upper bound is infinity. 
1159
    if degreeIntervalSa.lower() == degreeIntervalSa.upper():
1160
        return int(degreeIntervalSa.lower())
1161
    else:
1162
        if degreeIntervalSa.upper().is_infinity():
1163
            return None
1164
        else:
1165
            return int(degreeIntervalSa.upper())
1166
    # End pobyso_guess_degree_so_sa
1167

    
1168
def pobyso_inf_so_so(intervalSo):
1169
    """
1170
    Very thin wrapper around sollya_lib_inf().
1171
    """
1172
    return sollya_lib_inf(intervalSo)
1173
# End pobyso_inf_so_so.
1174
    
1175
def pobyso_infnorm_so_so(func, interval, file = None, intervalList = None):
1176
    print "Do not use this function. User pobyso_supnorm_so_so instead."
1177
    return None
1178

    
1179
def pobyso_interval_to_range_sa_so(intervalSa, precisionSa=None):
1180
    if precisionSa is None:
1181
        precisionSa = intervalSa.parent().precision()
1182
    intervalSo = pobyso_bounds_to_range_sa_so(intervalSa.lower(),\
1183
                                              intervalSa.upper(),\
1184
                                              precisionSa)
1185
    return intervalSo
1186
# End pobyso_interval_to_range_sa_so
1187

    
1188
def pobyso_is_error_so_sa(objSo):
1189
    """
1190
    Thin wrapper around the sollya_lib_obj_is_error() function.
1191
    """
1192
    if sollya_lib_obj_is_error(objSo) != 0:
1193
        return True
1194
    else:
1195
        return False
1196
# End pobyso_is_error-so_sa
1197

    
1198
def pobyso_is_floating_point_number_sa_sa(numberSa):
1199
    """
1200
    Check whether a Sage number is floating point.
1201
    Exception stuff added because numbers other than
1202
    floating-point ones do not have the is_real() attribute.
1203
    """
1204
    try:
1205
        return numberSa.is_real()
1206
    except AttributeError:
1207
        return False
1208
# End pobyso_is_floating_piont_number_sa_sa
1209

    
1210
def pobyso_is_range_so_sa(rangeCandidateSo):
1211
    """
1212
    Thin wrapper over sollya_lib_is_range.
1213
    """
1214
    return sollya_lib_obj_is_range(rangeCandidateSo) != 0
1215
    
1216
# End pobyso_is_range_so_sa
1217

    
1218

    
1219
def pobyso_lib_init():
1220
    sollya_lib_init(None)
1221

    
1222
def pobyso_lib_close():
1223
    sollya_lib_close(None)
1224

    
1225
def pobyso_name_free_variable(freeVariableNameSa):
1226
    """ Legacy function. See pobyso_name_free_variable_sa_so. """
1227
    pobyso_name_free_variable_sa_so(freeVariableNameSa)
1228

    
1229
def pobyso_name_free_variable_sa_so(freeVariableNameSa):
1230
    """
1231
    Set the free variable name in Sollya from a Sage string.
1232
    """
1233
    sollya_lib_name_free_variable(freeVariableNameSa)
1234

    
1235
def pobyso_parse_string(string):
1236
    """ Legacy function. See pobyso_parse_string_sa_so. """
1237
    return pobyso_parse_string_sa_so(string)
1238
 
1239
def pobyso_parse_string_sa_so(string):
1240
    """
1241
    Get the Sollya expression computed from a Sage string or
1242
    a Sollya error object if parsing failed.
1243
    """
1244
    return sollya_lib_parse_string(string)
1245

    
1246
def pobyso_precision_so_sa(ctExpSo):
1247
    """
1248
    Computes the necessary precision to represent a number.
1249
    If x is not zero, it can be uniquely written as x = m · 2e 
1250
    where m is an odd integer and e is an integer. 
1251
    precision(x) returns the number of bits necessary to write m 
1252
    in binary (i.e. ceil(log2(m))).
1253
    """
1254
    #TODO: take care of the special case: 0, @NaN@, @Inf@
1255
    precisionSo = sollya_lib_precision(ctExpSo)
1256
    precisionSa = pobyso_constant_from_int_so_sa(precisionSo)
1257
    sollya_lib_clear_obj(precisionSo)
1258
    return precisionSa
1259
# End pobyso_precision_so_sa
1260

    
1261
def pobyso_polynomial_coefficients_progressive_round_so_so(polySo,
1262
                                                           funcSo,
1263
                                                           icSo,
1264
                                                           intervalSo,
1265
                                                           itpSo,
1266
                                                           ftpSo,
1267
                                                           maxPrecSo,
1268
                                                           maxErrSo,
1269
                                                           debug=False):
1270
    if debug:
1271
        print "Input arguments:"
1272
        pobyso_autoprint(polySo)
1273
        pobyso_autoprint(funcSo)
1274
        pobyso_autoprint(icSo)
1275
        pobyso_autoprint(intervalSo)
1276
        pobyso_autoprint(itpSo)
1277
        pobyso_autoprint(ftpSo)
1278
        pobyso_autoprint(maxPrecSo)
1279
        pobyso_autoprint(maxErrSo)
1280
        print "________________"
1281
    
1282
    ## Higher order function see: 
1283
    #  http://effbot.org/pyfaq/how-do-you-make-a-higher-order-function-in-python.htm
1284
    def precision_decay_ratio_function(degreeSa):
1285
        def outer(x):
1286
            def inner(x):
1287
                we = 3/8
1288
                wq = 2/8
1289
                a  = 2.2
1290
                b  = 2 
1291
                return we*(exp(x/a)-1) +  wq*((b*x)**2) + (1-we-wq)*x
1292
            return  inner(x)/inner(degreeSa)
1293
        return outer
1294
    
1295
    #
1296
    degreeSa        = pobyso_polynomial_degree_so_sa(polySo)
1297
    ratio           = precision_decay_ratio_function(degreeSa)
1298
    itpSa           = pobyso_constant_from_int_so_sa(itpSo)
1299
    ftpSa           = pobyso_constant_from_int_so_sa(ftpSo)
1300
    maxPrecSa       = pobyso_constant_from_int_so_sa(maxPrecSo)
1301
    maxErrSa        = pobyso_get_constant_as_rn_so_sa(maxErrSo)
1302
    if debug:
1303
        print "degreeSa:", degreeSa
1304
        print "ratio:", ratio
1305
        print "itpsSa:", itpSa
1306
        print "ftpSa:", ftpSa
1307
        print "maxPrecSa:", maxPrecSa
1308
        print "maxErrSa:", maxErrSa
1309
    lastResPolySo   = None
1310
    lastInfNormSo   = None
1311
    #print "About to enter the while loop..."
1312
    while True: 
1313
        resPolySo   = pobyso_constant_0_sa_so()
1314
        pDeltaSa    = ftpSa - itpSa
1315
        for indexSa in reversed(xrange(0,degreeSa+1)):
1316
            #print "Index:", indexSa
1317
            indexSo = pobyso_constant_from_int_sa_so(indexSa)
1318
            #pobyso_autoprint(indexSo)
1319
            #print ratio(indexSa)
1320
            ctpSa = floor(ftpSa - (pDeltaSa * ratio(indexSa)))
1321
            ctpSo = pobyso_constant_from_int_sa_so(ctpSa)
1322
            if debug:
1323
                print "Index:", indexSa, " - Target precision:", 
1324
                pobyso_autoprint(ctpSo)
1325
            cmonSo  = \
1326
                sollya_lib_build_function_mul(sollya_lib_coeff(polySo, indexSo),
1327
                                      sollya_lib_build_function_pow( \
1328
                                          sollya_lib_build_function_free_variable(), \
1329
                                          indexSo))
1330
            #pobyso_autoprint(cmonSo)
1331
            cmonrSo = pobyso_round_coefficients_single_so_so(cmonSo, ctpSo)
1332
            sollya_lib_clear_obj(cmonSo)
1333
            #pobyso_autoprint(cmonrSo)
1334
            resPolySo = sollya_lib_build_function_add(resPolySo,
1335
                                                      cmonrSo)
1336
            #pobyso_autoprint(resPolySo)
1337
        # End for index
1338
        freeVarSo     = sollya_lib_build_function_free_variable()
1339
        changeVarSo   = sollya_lib_sub(freeVarSo, icSo)
1340
        resPolyCvSo   = sollya_lib_evaluate(resPolySo, changeVarSo)
1341
        errFuncSo = sollya_lib_build_function_sub(sollya_lib_copy_obj(funcSo), 
1342
                                                  resPolyCvSo)
1343
        infNormSo = sollya_lib_dirtyinfnorm(errFuncSo, intervalSo)
1344
        cerrSa    = pobyso_get_constant_as_rn_so_sa(infNormSo)
1345
        if debug:
1346
            print "Infnorm (Sollya):", 
1347
            pobyso_autoprint(infNormSo)
1348
        sollya_lib_clear_obj(errFuncSo)
1349
        #print "Infnorm  (Sage):", cerrSa
1350
        if (cerrSa > maxErrSa):
1351
            if debug:
1352
                print "Error is too large."
1353
            if lastResPolySo is None:
1354
                if debug:
1355
                    print "Enlarging prec."
1356
                ntpSa = floor(ftpSa + ftpSa/50)
1357
                ## Can't enlarge (numerical)
1358
                if ntpSa == ftpSa:
1359
                    sollya_lib_clear_obj(resPolySo)
1360
                    return None
1361
                ## Can't enlarge (not enough precision left)
1362
                if ntpSa > maxPrecSa:
1363
                    sollya_lib_clear_obj(resPolySo)
1364
                    return None
1365
                ftpSa = ntpSa
1366
                continue
1367
            ## One enlargement took place.
1368
            else:
1369
                if debug:
1370
                    print "Exit with the last before last polynomial."
1371
                    print "Precision of highest degree monomial:", itpSa
1372
                    print "Precision of constant term          :", ftpSa
1373
                sollya_lib_clear_obj(resPolySo)
1374
                sollya_lib_clear_obj(infNormSo)
1375
                return (lastResPolySo, lastInfNormSo)
1376
        # cerrSa <= maxErrSa: scrap more bits, possibly.
1377
        else:
1378
            if debug:
1379
                print "Error is too small"
1380
            if cerrSa <= (maxErrSa/2):
1381
                if debug: 
1382
                    print "Shrinking prec."
1383
                ntpSa = floor(ftpSa - ftpSa/50)
1384
                ## Can't shrink (numerical)
1385
                if ntpSa == ftpSa:
1386
                    if not lastResPolySo is None:
1387
                        sollya_lib_clear_obj(lastResPolySo)
1388
                    if not lastInfNormSo is None:
1389
                        sollya_lib_clear_obj(lastInfNormSo)
1390
                    if debug:
1391
                        print "Exit because can't shrink anymore (numerically)."
1392
                        print "Precision of highest degree monomial:", itpSa
1393
                        print "Precision of constant term          :", ftpSa
1394
                    return (resPolySo, infNormSo)
1395
                ## Can't shrink (not enough precision left)
1396
                if ntpSa <= itpSa:
1397
                    if not lastResPolySo is None:
1398
                        sollya_lib_clear_obj(lastResPolySo)
1399
                    if not lastInfNormSo is None:
1400
                        sollya_lib_clear_obj(lastInfNormSo)
1401
                        print "Exit because can't shrink anymore (no bits left)."
1402
                        print "Precision of highest degree monomial:", itpSa
1403
                        print "Precision of constant term          :", ftpSa
1404
                    return (resPolySo, infNormSo)
1405
                ftpSa = ntpSa
1406
                if not lastResPolySo is None:
1407
                    sollya_lib_clear_obj(lastResPolySo)
1408
                if not lastInfNormSo is None:
1409
                    sollya_lib_clear_obj(lastInfNormSo)
1410
                lastResPolySo = resPolySo
1411
                lastInfNormSo = infNormSo
1412
                continue
1413
            else: # Error is not that small, just return
1414
                if not lastResPolySo is None:
1415
                    sollya_lib_clear_obj(lastResPolySo)
1416
                if not lastInfNormSo is None:
1417
                    sollya_lib_clear_obj(lastInfNormSo)
1418
                if debug:
1419
                    print "Exit normally."
1420
                    print "Precision of highest degree monomial:", itpSa
1421
                    print "Precision of constant term          :", ftpSa
1422
                return (resPolySo, infNormSo)                
1423
    # End wile True
1424
    return None
1425
# End pobyso_polynomial_coefficients_progressive_truncate_so_so.
1426

    
1427
def pobyso_polynomial_degree_so_sa(polySo):
1428
    """
1429
    Return the degree of a Sollya polynomial as a Sage int.
1430
    """
1431
    degreeSo = sollya_lib_degree(polySo)
1432
    return pobyso_constant_from_int_so_sa(degreeSo)
1433
# End pobyso_polynomial_degree_so_sa
1434

    
1435
def pobyso_polynomial_degree_so_so(polySo):
1436
    """
1437
    Thin wrapper around lib_sollya_degree().
1438
    """
1439
    return sollya_lib_degree(polySo)
1440
# End pobyso_polynomial_degree_so_so
1441
                                                                  
1442
def pobyso_range(rnLowerBound, rnUpperBound):
1443
    """ Legacy function. See pobyso_range_sa_so. """
1444
    return pobyso_range_sa_so(rnLowerBound, rnUpperBound) 
1445

    
1446
def pobyso_range_from_bounds_sa_so(rnLowerBound, rnUpperBound, precSa = None):
1447
    """
1448
    Create a Sollya range from 2 Sage real numbers as bounds
1449
    """
1450
    # TODO check precision stuff.
1451
    sollyaPrecChanged = False
1452
    (initialSollyaPrecSo, initialSollyaPrecSa) = \
1453
        pobyso_get_prec_so_so_sa()
1454
    if precSa is None:
1455
        precSa = max(rnLowerBound.parent().prec(), rnUpperBound.parent().prec())
1456
    if precSa > initialSollyaPrecSa:
1457
        if precSa <= 2:
1458
            print inspect.stack()[0][3], ": precision change <= 2 requested."
1459
        pobyso_set_prec_sa_so(precSa)
1460
        sollyaPrecChanged = True
1461
    rangeSo = sollya_lib_range_from_bounds(get_rn_value(rnLowerBound),
1462
                                           get_rn_value(rnUpperBound))
1463
    if sollyaPrecChanged:
1464
        pobyso_set_prec_so_so(initialSollyaPrecSo)
1465
        sollya_lib_clear_obj(initialSollyaPrecSo)
1466
    return rangeSo
1467
# End pobyso_range_from_bounds_sa_so
1468

    
1469
def pobyso_range_max_abs_so_so(rangeSo):
1470
    """
1471
    Return, as Sollya constant, the maximum absolute value of a Sollay range.
1472
    """
1473
    lowerBoundSo = sollya_lib_inf(rangeSo)
1474
    upperBoundSo = sollya_lib_sup(rangeSo)
1475
    #
1476
    lowerBoundSo = sollya_lib_build_function_abs(lowerBoundSo)
1477
    upperBoundSo = sollya_lib_build_function_abs(upperBoundSo)
1478
    #pobyso_autoprint(lowerBoundSo)
1479
    #pobyso_autoprint(upperBoundSo)
1480
    #
1481
    maxAbsSo = sollya_lib_max(lowerBoundSo, upperBoundSo, None)
1482
    #sollya_lib_clear_obj(lowerBoundSo)
1483
    #sollya_lib_clear_obj(upperBoundSo)
1484
    return maxAbsSo
1485
# End pobyso_range_max_abs_so_so
1486
 
1487
def pobyso_range_to_interval_so_sa(rangeSo, realIntervalFieldSa = None):
1488
    """
1489
    Get a Sage interval from a Sollya range.
1490
    If no realIntervalField is given as a parameter, the Sage interval
1491
    precision is that of the Sollya range.
1492
    Otherwise, the precision is that of the realIntervalField. In this case
1493
    rounding may happen.
1494
    """
1495
    if realIntervalFieldSa is None:
1496
        precSa = pobyso_get_prec_of_range_so_sa(rangeSo)
1497
        realIntervalFieldSa = RealIntervalField(precSa)
1498
    intervalSa = \
1499
        pobyso_get_interval_from_range_so_sa(rangeSo, realIntervalFieldSa) 
1500
    return intervalSa
1501
# End pobyso_range_to_interval_so_sa
1502
#
1503
def pobyso_relative_so_so():
1504
    """
1505
    Very thin wrapper around the sollya_lib_relative function.
1506
    """
1507
    return sollya_lib_relative()
1508
# End pobyso_relative_so_so
1509
#
1510
def pobyso_rat_poly_sa_so(polySa, precSa = None):
1511
    """
1512
    Create a Sollya polynomial from a Sage rational polynomial.
1513
    """
1514
    ## TODO: filter arguments.
1515
    ## Precision. If no precision is given, use the current precision
1516
    #  of Sollya.
1517
    if precSa is None:
1518
        precSa =  pobyso_get_prec_so_sa()
1519
    #print "Precision:",  precSa
1520
    RRR = RealField(precSa)
1521
    ## Create a Sage polynomial in the "right" precision.
1522
    P_RRR = RRR[polySa.variables()[0]]
1523
    polyFloatSa = P_RRR(polySa)
1524
    ## Make sure no precision is provided: pobyso_float_poly_sa_so will
1525
    #  recover it all by itself and not make an extra conversion.
1526
    return pobyso_float_poly_sa_so(polyFloatSa)
1527
    
1528
# End pobyso_rat_poly_sa_so    
1529
    
1530
def pobyso_remez_canonical_sa_sa(func, \
1531
                                 degree, \
1532
                                 lowerBound, \
1533
                                 upperBound, \
1534
                                 weight = None, \
1535
                                 quality = None):
1536
    """
1537
    All arguments are Sage/Python.
1538
    The functions (func and weight) must be passed as expressions or strings.
1539
    Otherwise the function fails. 
1540
    The return value is a Sage polynomial.
1541
    """
1542
    var('zorglub')    # Dummy variable name for type check only. Type of 
1543
    # zorglub is "symbolic expression".
1544
    polySo = pobyso_remez_canonical_sa_so(func, \
1545
                                 degree, \
1546
                                 lowerBound, \
1547
                                 upperBound, \
1548
                                 weight, \
1549
                                 quality)
1550
    # String test
1551
    if parent(func) == parent("string"):
1552
        functionSa = eval(func)
1553
    # Expression test.
1554
    elif type(func) == type(zorglub):
1555
        functionSa = func
1556
    else:
1557
        return None
1558
    #
1559
    maxPrecision = 0
1560
    if polySo is None:
1561
        return(None)
1562
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
1563
    RRRRSa = RealField(maxPrecision)
1564
    polynomialRingSa = RRRRSa[functionSa.variables()[0]]
1565
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, RRRRSa)
1566
    polySa = polynomial(expSa, polynomialRingSa)
1567
    sollya_lib_clear_obj(polySo)
1568
    return(polySa)
1569
# End pobyso_remez_canonical_sa_sa
1570
    
1571
def pobyso_remez_canonical(func, \
1572
                           degree, \
1573
                           lowerBound, \
1574
                           upperBound, \
1575
                           weight = "1", \
1576
                           quality = None):
1577
    """ Legacy function. See pobyso_remez_canonical_sa_so. """
1578
    return(pobyso_remez_canonical_sa_so(func, \
1579
                                        degree, \
1580
                                        lowerBound, \
1581
                                        upperBound, \
1582
                                        weight, \
1583
                                        quality))
1584
# End pobyso_remez_canonical.
1585

    
1586
def pobyso_remez_canonical_sa_so(func, \
1587
                                 degree, \
1588
                                 lowerBound, \
1589
                                 upperBound, \
1590
                                 weight = None, \
1591
                                 quality = None):
1592
    """
1593
    All arguments are Sage/Python.
1594
    The functions (func and weight) must be passed as expressions or strings.
1595
    Otherwise the function fails. 
1596
    The return value is a pointer to a Sollya function.
1597
    lowerBound and upperBound mus be reals.
1598
    """
1599
    var('zorglub')    # Dummy variable name for type check only. Type of
1600
    # zorglub is "symbolic expression".
1601
    currentVariableNameSa = None
1602
    # The func argument can be of different types (string, 
1603
    # symbolic expression...)
1604
    if parent(func) == parent("string"):
1605
        localFuncSa = sage_eval(func,globals())
1606
        if len(localFuncSa.variables()) > 0:
1607
            currentVariableNameSa = localFuncSa.variables()[0]
1608
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1609
            functionSo = \
1610
              sollya_lib_parse_string(localFuncSa._assume_str().replace('_SAGE_VAR_', ''))
1611
    # Expression test.
1612
    elif type(func) == type(zorglub):
1613
        # Until we are able to translate Sage expressions into Sollya 
1614
        # expressions : parse the string version.
1615
        if len(func.variables()) > 0:
1616
            currentVariableNameSa = func.variables()[0]
1617
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1618
            functionSo = \
1619
              sollya_lib_parse_string(func._assume_str().replace('_SAGE_VAR_', ''))
1620
    else:
1621
        return(None)
1622
    if weight is None: # No weight given -> 1.
1623
        weightSo = pobyso_constant_1_sa_so()
1624
    elif parent(weight) == parent("string"): # Weight given as string: parse it.
1625
        weightSo = sollya_lib_parse_string(func)
1626
    elif type(weight) == type(zorglub): # Weight given as symbolice expression.
1627
        functionSo = \
1628
          sollya_lib_parse_string_sa_so(weight._assume_str().replace('_SAGE_VAR_', ''))
1629
    else:
1630
        return(None)
1631
    degreeSo = pobyso_constant_from_int(degree)
1632
    rangeSo = pobyso_bounds_to_range_sa_so(lowerBound, upperBound)
1633
    if not quality is None:
1634
        qualitySo= pobyso_constant_sa_so(quality)
1635
    else:
1636
        qualitySo = None
1637
        
1638
    remezPolySo = sollya_lib_remez(functionSo, \
1639
                                   degreeSo, \
1640
                                   rangeSo, \
1641
                                   weightSo, \
1642
                                   qualitySo, \
1643
                                   None)
1644
    sollya_lib_clear_obj(functionSo)
1645
    sollya_lib_clear_obj(degreeSo)
1646
    sollya_lib_clear_obj(rangeSo)
1647
    sollya_lib_clear_obj(weightSo)
1648
    if not qualitySo is None:
1649
        sollya_lib_clear_obj(qualitySo)
1650
    return(remezPolySo)
1651
# End pobyso_remez_canonical_sa_so
1652

    
1653
def pobyso_remez_canonical_so_so(funcSo, \
1654
                                 degreeSo, \
1655
                                 rangeSo, \
1656
                                 weightSo = pobyso_constant_1_sa_so(),\
1657
                                 qualitySo = None):
1658
    """
1659
    All arguments are pointers to Sollya objects.
1660
    The return value is a pointer to a Sollya function.
1661
    """
1662
    if not sollya_lib_obj_is_function(funcSo):
1663
        return(None)
1664
    return(sollya_lib_remez(funcSo, degreeSo, rangeSo, weightSo, qualitySo, None))
1665
# End pobyso_remez_canonical_so_so.
1666

    
1667
def pobyso_remez_exponents_list_sa_so(func,     \
1668
                                 exponentsList, \
1669
                                 lowerBound,    \
1670
                                 upperBound,    \
1671
                                 weight = None, \
1672
                                 quality = None):
1673
    """
1674
    All arguments are Sage/Python.
1675
    The functions (func and weight) must be passed as expressions or strings.
1676
    Otherwise the function fails. 
1677
    The return value is a pointer to a Sollya function.
1678
    lowerBound and upperBound mus be reals.
1679
    """
1680
    var('zorglub')    # Dummy variable name for type check only. Type of
1681
    # zorglub is "symbolic expression".
1682
    currentVariableNameSa = None
1683
    # The func argument can be of different types (string, 
1684
    # symbolic expression...)
1685
    if parent(func) == parent("string"):
1686
        localFuncSa = sage_eval(func,globals())
1687
        if len(localFuncSa.variables()) > 0:
1688
            currentVariableNameSa = localFuncSa.variables()[0]
1689
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1690
            functionSo = \
1691
              sollya_lib_parse_string(localFuncSa._assume_str().replace('_SAGE_VAR_', ''))
1692
    # Expression test.
1693
    elif type(func) == type(zorglub):
1694
        # Until we are able to translate Sage expressions into Sollya 
1695
        # expressions : parse the string version.
1696
        if len(func.variables()) > 0:
1697
            currentVariableNameSa = func.variables()[0]
1698
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1699
            functionSo = \
1700
              sollya_lib_parse_string(func._assume_str().replace('_SAGE_VAR_', ''))
1701
    else:
1702
        return(None)
1703
    ## Deal with the weight, much in the same way as with the function.
1704
    if weight is None: # No weight given -> 1.
1705
        weightSo = pobyso_constant_1_sa_so()
1706
    elif parent(weight) == parent("string"): # Weight given as string: parse it.
1707
        weightSo = sollya_lib_parse_string(func)
1708
    elif type(weight) == type(zorglub): # Weight given as symbolice expression.
1709
        functionSo = \
1710
          sollya_lib_parse_string(weight._assume_str().replace('_SAGE_VAR_', '',100))
1711
    else:
1712
        return(None)
1713
    rangeSo = pobyso_bounds_to_range_sa_so(lowerBound, upperBound)
1714
    if not quality is None:
1715
        qualitySo= pobyso_constant_sa_so(quality)
1716
    else:
1717
        qualitySo = None
1718
    #
1719
    ## Tranform the Sage list of exponents into a Sollya list.
1720
    exponentsListSo = pobyso_build_list_of_ints_sa_so(*exponentsList)
1721
    remezPolySo = sollya_lib_remez(functionSo, \
1722
                                   exponentsListSo, \
1723
                                   rangeSo, \
1724
                                   weightSo, \
1725
                                   qualitySo, \
1726
                                   None)
1727
    sollya_lib_clear_obj(functionSo)
1728
    sollya_lib_clear_obj(exponentsListSo)
1729
    sollya_lib_clear_obj(rangeSo)
1730
    sollya_lib_clear_obj(weightSo)
1731
    if not qualitySo is None:
1732
        sollya_lib_clear_obj(qualitySo)
1733
    return(remezPolySo)
1734
# End pobyso_remez_exponentsList_sa_so
1735
#
1736
def pobyso_round_coefficients_so_so(polySo, truncFormatListSo):
1737
    """
1738
    A wrapper around the "classical" sollya_lib_roundcoefficients: a Sollya
1739
    polynomial and a Sollya list are given as arguments.
1740
    """
1741
    return sollya_lib_roundcoefficients(polySo, truncFormatListSo)
1742

    
1743
def pobyso_round_coefficients_progressive_so_so(polySo, 
1744
                                                funcSo,
1745
                                                precSo,
1746
                                                intervalSo,
1747
                                                icSo,
1748
                                                currentApproxErrorSo,
1749
                                                approxAccurSo,
1750
                                                debug=False):
1751
    if debug:
1752
        print "Input arguments:"
1753
        print "Polynomial: ", ; pobyso_autoprint(polySo)
1754
        print "Function: ", ; pobyso_autoprint(funcSo)
1755
        print "Internal precision: ", ; pobyso_autoprint(precSo)
1756
        print "Interval: ", ; pobyso_autoprint(intervalSo)
1757
        print "Current approximation error: ", ; pobyso_autoprint(currentApproxErrorSo)
1758
        print "Requested approxiation error: ", ; pobyso_autoprint(approxAccurSo)
1759
        print "________________"
1760
    approxAccurSa        = pobyso_get_constant_as_rn_so_sa(approxAccurSo)
1761
    currentApproxErrorSa = pobyso_get_constant_as_rn_so_sa(currentApproxErrorSo)
1762
    ## If the current approximation error is too close to the target, there is
1763
    #  no possible gain.
1764
    if currentApproxErrorSa >= approxAccurSa / 2:
1765
        #### Do not return the initial argument but copies: caller may free 
1766
        #    as inutile after call.
1767
        return (sollya_lib_copy_obj(polySo),
1768
                sollya_lib_copy_obj(currentApproxErrorSo))
1769
    degreeSa             = pobyso_polynomial_degree_so_sa(polySo)
1770
    intervalSa           = pobyso_range_to_interval_so_sa(intervalSo)
1771
    
1772
    if debug:
1773
        print "degreeSa              :", degreeSa
1774
        print "intervalSa            :", intervalSa.str(style='brackets')
1775
        print "currentApproxErrorSa  :", currentApproxErrorSa 
1776
        print "approxAccurSa         :", approxAccurSa 
1777
    ### Start with a 0 value expression.
1778
    radiusSa = intervalSa.absolute_diameter() / 2
1779
    if debug:
1780
        print "log2(radius):", RR(radiusSa).log2()
1781
    iterIndex = 0
1782
    while True: 
1783
        resPolySo = pobyso_constant_0_sa_so()
1784
        roundedPolyApproxAccurSa = approxAccurSa / 2
1785
        currentRadiusPowerSa = 1 
1786
        for degree in xrange(0,degreeSa + 1):
1787
            #### At round 0, use the agressive formula. At round 1, run the
1788
            #    proved formula.
1789
            if iterIndex == 0:
1790
                roundingPowerSa = \
1791
                    floor(((currentRadiusPowerSa/roundedPolyApproxAccurSa)*(degree+1)).log2())
1792
            else:
1793
                roundingPowerSa = \
1794
                    floor(((currentRadiusPowerSa/roundedPolyApproxAccurSa)*(degreeSa+1)).log2())
1795
            ## Under extreme conditions the above formulas can evaluate under 2, which is the
1796
            #  minimal precision of an MPFR number.
1797
            if roundingPowerSa < 2:
1798
                roundingPowerSa = 2
1799
            if debug:
1800
                print "roundedPolyApproxAccurSa", roundedPolyApproxAccurSa
1801
                print "currentRadiusPowerSa", currentRadiusPowerSa
1802
                print "Current rounding exponent:", roundingPowerSa
1803
            currentRadiusPowerSa *= radiusSa
1804
            index1So = pobyso_constant_from_int_sa_so(degree)
1805
            index2So = pobyso_constant_from_int_sa_so(degree)
1806
            ### Create a monomial with:
1807
            #   - the coefficient in the initial monomial at the current degrree;
1808
            #   - the current exponent;
1809
            #   - the free variable.
1810
            cmonSo  = \
1811
                sollya_lib_build_function_mul(sollya_lib_coeff(polySo, index1So),
1812
                                              sollya_lib_build_function_pow( \
1813
                                                sollya_lib_build_function_free_variable(), \
1814
                                                index2So))
1815
            roundingPowerSo = pobyso_constant_from_int_sa_so(roundingPowerSa)
1816
            cmonrSo = pobyso_round_coefficients_single_so_so(cmonSo, roundingPowerSo)
1817
            sollya_lib_clear_obj(cmonSo)
1818
            ### Add to the result polynomial.
1819
            resPolySo = sollya_lib_build_function_add(resPolySo,
1820
                                                      cmonrSo)
1821
        # End for.
1822
        ### Check the new polynomial.
1823
        freeVarSo     = sollya_lib_build_function_free_variable()
1824
        changeVarSo   = sollya_lib_sub(freeVarSo, icSo)
1825
        resPolyCvSo   = sollya_lib_evaluate(resPolySo, changeVarSo)
1826
        errFuncSo     = sollya_lib_build_function_sub(sollya_lib_copy_obj(funcSo), 
1827
                                                      resPolyCvSo)
1828
        infNormSo = sollya_lib_dirtyinfnorm(errFuncSo, intervalSo)
1829
        ### This also clears resPolyCvSo.
1830
        sollya_lib_clear_obj(errFuncSo)
1831
        cerrSa    = pobyso_get_constant_as_rn_so_sa(infNormSo)
1832
        if debug:
1833
            print "Error of the new polynomial:", cerrSa
1834
        ### If at round 1, return the initial polynomial error. This should
1835
        #   never happen since the rounding algorithm is proved. But some 
1836
        #   circumstances may break it (e.g. internal precision of tools).
1837
        if cerrSa > approxAccurSa:
1838
            if iterIndex > 0: # Round 1 and beyond.
1839
                sollya_lib_clear_obj(resPolySo)
1840
                sollya_lib_clear_obj(infNormSo)
1841
                #### Do not return the arguments but copies: the caller may free
1842
                #    free the former as inutile after call.
1843
                return (sollya_lib_copy_obj(polySo), 
1844
                        sollya_lib_copy_obj(currentApproxErrorSo))
1845
            else: # Round 0, got round 1
1846
                sollya_lib_clear_obj(resPolySo)
1847
                sollya_lib_clear_obj(infNormSo)
1848
                iterIndex += 1
1849
                continue
1850
        ### If get here it is because cerrSa <= approxAccurSa
1851
        ### Approximation error of the new polynomial is acceptable.
1852
        return (resPolySo, infNormSo)
1853
    # End while True
1854
# End pobyso_round_coefficients_progressive_so_so
1855
    
1856
def pobyso_round_coefficients_single_so_so(polySo, commonPrecSo):
1857
    """
1858
    Create a rounded coefficients polynomial from polynomial argument to
1859
    the number of bits in size argument.
1860
    All coefficients are set to the same precision.
1861
    """
1862
    ## TODO: check arguments.
1863
    endEllipListSo = pobyso_build_end_elliptic_list_so_so(commonPrecSo)
1864
    polySo = sollya_lib_roundcoefficients(polySo, endEllipListSo, None)
1865
    sollya_lib_clear_obj(endEllipListSo)
1866
    #sollya_lib_clear_obj(endEllipListSo)
1867
    return polySo
1868
    
1869
# End pobyso_round_coefficients_single_so_so
1870

    
1871
def pobyso_set_canonical_off():
1872
    sollya_lib_set_canonical(sollya_lib_off())
1873

    
1874
def pobyso_set_canonical_on():
1875
    sollya_lib_set_canonical(sollya_lib_on())
1876

    
1877
def pobyso_set_prec(p):
1878
    """ Legacy function. See pobyso_set_prec_sa_so. """
1879
    pobyso_set_prec_sa_so(p)
1880

    
1881
def pobyso_set_prec_sa_so(p):
1882
    #a = c_int(p)
1883
    #precSo = c_void_p(sollya_lib_constant_from_int(a))
1884
    #precSo = sollya_lib_constant_from_int(a)
1885
    precSo = pobyso_constant_from_int_sa_so(p)
1886
    sollya_lib_set_prec(precSo)
1887
    sollya_lib_clear_obj(precSo)
1888
# End pobyso_set_prec_sa_so.
1889

    
1890
def pobyso_set_prec_so_so(newPrecSo):
1891
    sollya_lib_set_prec(newPrecSo)
1892
# End pobyso_set_prec_so_so.
1893

    
1894
def pobyso_inf_so_so(intervalSo):
1895
    """
1896
    Very thin wrapper around sollya_lib_inf().
1897
    """
1898
    return sollya_lib_inf(intervalSo)
1899
# End pobyso_inf_so_so.
1900
#   
1901
def pobyso_supnorm_so_sa(polySo, funcSo, intervalSo, errorTypeSo = None,\
1902
                         accuracySo = None, realFieldSa = None):
1903
    """
1904
    Computes the supremum norm from Solly input arguments and returns a
1905
    Sage floating-point number whose precision is set by the only Sage argument.
1906
    
1907
    The returned value is the maximum of the absolute values of the range
1908
    elements  returned by the Sollya supnorm functions
1909
    """
1910
    supNormRangeSo = pobyso_supnorm_so_so(polySo, 
1911
                                          funcSo, 
1912
                                          intervalSo, 
1913
                                          errorTypeSo,
1914
                                          accuracySo)
1915
    supNormSo = pobyso_range_max_abs_so_so(supNormRangeSo)
1916
    sollya_lib_clear_obj(supNormRangeSo)
1917
    #pobyso_autoprint(supNormSo)
1918
    supNormSa = pobyso_get_constant_as_rn_with_rf_so_sa(supNormSo, realFieldSa)
1919
    sollya_lib_clear_obj(supNormSo)
1920
    return supNormSa 
1921
# End pobyso_supnorm_so_sa.
1922
#
1923
def pobyso_supnorm_so_so(polySo, funcSo, intervalSo, errorTypeSo = None,\
1924
                         accuracySo = None):
1925
    """
1926
    Computes the supnorm of the approximation error between the given 
1927
    polynomial and function. Attention: returns a range!
1928
    errorTypeSo defaults to "absolute".
1929
    accuracySo defaults to 2^(-40).
1930
    """
1931
    if errorTypeSo is None:
1932
        errorTypeSo = sollya_lib_absolute(None)
1933
        errorTypeIsNone = True
1934
    else:
1935
        errorTypeIsNone = False
1936
    #
1937
    if accuracySo is None:
1938
        # Notice the **: we are in Pythonland here!
1939
        accuracySo = pobyso_constant_sa_so(RR(2**(-40)))
1940
        accuracyIsNone = True
1941
    else:
1942
        accuracyIsNone = False
1943
    #pobyso_autoprint(accuracySo)
1944
    resultSo = \
1945
        sollya_lib_supnorm(polySo, funcSo, intervalSo, errorTypeSo, \
1946
                              accuracySo)
1947
    if errorTypeIsNone:
1948
        sollya_lib_clear_obj(errorTypeSo)
1949
    if accuracyIsNone:
1950
        sollya_lib_clear_obj(accuracySo)
1951
    return resultSo
1952
# End pobyso_supnorm_so_so
1953
#
1954
def pobyso_taylor_expansion_no_change_var_so_so(functionSo, 
1955
                                                degreeSo, 
1956
                                                rangeSo,
1957
                                                errorTypeSo=None, 
1958
                                                sollyaPrecSo=None):
1959
    """
1960
    Compute the Taylor expansion without the variable change
1961
    x -> x-intervalCenter.
1962
    """
1963
    # Change internal Sollya precision, if needed.
1964
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
1965
    sollyaPrecChanged = False
1966
    if sollyaPrecSo is None:
1967
        pass
1968
    else:
1969
        sollya_lib_set_prec(sollyaPrecSo)
1970
        sollyaPrecChanged = True
1971
    # Error type stuff: default to absolute.
1972
    if errorTypeSo is None:
1973
        errorTypeIsNone = True
1974
        errorTypeSo = sollya_lib_absolute(None)
1975
    else:
1976
        errorTypeIsNone = False
1977
    intervalCenterSo = sollya_lib_mid(rangeSo, None)
1978
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo,
1979
                                         intervalCenterSo,
1980
                                         rangeSo, errorTypeSo, None)
1981
    # taylorFormListSaSo is a Python list of Sollya objects references that 
1982
    # are copies of the elements of taylorFormSo.
1983
    # pobyso_get_list_elements_so_so clears taylorFormSo.
1984
    (taylorFormListSaSo, numElementsSa, isEndEllipticSa) = \
1985
        pobyso_get_list_elements_so_so(taylorFormSo)
1986
    polySo = sollya_lib_copy_obj(taylorFormListSaSo[0])
1987
    #print "Num elements:", numElementsSa
1988
    sollya_lib_clear_obj(taylorFormSo)
1989
    #polySo = taylorFormListSaSo[0]
1990
    #errorRangeSo = sollya_lib_copy_obj(taylorFormListSaSo[2])
1991
    errorRangeSo = taylorFormListSaSo[2]
1992
    # No copy_obj needed here: a new objects are created.
1993
    maxErrorSo    = sollya_lib_sup(errorRangeSo)
1994
    minErrorSo    = sollya_lib_inf(errorRangeSo)
1995
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
1996
    absMinErrorSo = sollya_lib_abs(minErrorSo)
1997
    sollya_lib_clear_obj(maxErrorSo)
1998
    sollya_lib_clear_obj(minErrorSo)
1999
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
2000
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
2001
    # If changed, reset the Sollya working precision.
2002
    if sollyaPrecChanged:
2003
        sollya_lib_set_prec(initialSollyaPrecSo)
2004
    sollya_lib_clear_obj(initialSollyaPrecSo)
2005
    if errorTypeIsNone:
2006
        sollya_lib_clear_obj(errorTypeSo)
2007
    pobyso_clear_taylorform_sa_so(taylorFormListSaSo)
2008
    if absMaxErrorSa > absMinErrorSa:
2009
        sollya_lib_clear_obj(absMinErrorSo)
2010
        return (polySo, intervalCenterSo, absMaxErrorSo)
2011
    else:
2012
        sollya_lib_clear_obj(absMaxErrorSo)
2013
        return (polySo, intervalCenterSo, absMinErrorSo)
2014
# end pobyso_taylor_expansion_no_change_var_so_so
2015

    
2016
def pobyso_taylor_expansion_with_change_var_so_so(functionSo, degreeSo, \
2017
                                                  rangeSo, \
2018
                                                  errorTypeSo=None, \
2019
                                                  sollyaPrecSo=None):
2020
    """
2021
    Compute the Taylor expansion with the variable change
2022
    x -> (x-intervalCenter) included.
2023
    """
2024
    # Change Sollya internal precision, if need.
2025
    sollyaPrecChanged = False
2026
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_sos_sa()
2027
    if sollyaPrecSo is None:
2028
        pass
2029
    else:
2030
        sollya_lib_set_prec(sollyaPrecSo)
2031
        sollyaPrecChanged = True
2032
    #
2033
    # Error type stuff: default to absolute.
2034
    if errorTypeSo is None:
2035
        errorTypeIsNone = True
2036
        errorTypeSo = sollya_lib_absolute(None)
2037
    else:
2038
        errorTypeIsNone = False
2039
    intervalCenterSo = sollya_lib_mid(rangeSo)
2040
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo, \
2041
                                         intervalCenterSo, \
2042
                                         rangeSo, errorTypeSo, None)
2043
    # taylorFormListSaSo is a Python list of Sollya objects references that 
2044
    # are copies of the elements of taylorFormSo.
2045
    # pobyso_get_list_elements_so_so clears taylorFormSo.
2046
    (taylorFormListSo, numElements, isEndElliptic) = \
2047
        pobyso_get_list_elements_so_so(taylorFormSo)
2048
    polySo = taylorFormListSo[0]
2049
    errorRangeSo = taylorFormListSo[2]
2050
    maxErrorSo    = sollya_lib_sup(errorRangeSo)
2051
    minErrorSo    = sollya_lib_inf(errorRangeSo)
2052
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
2053
    absMinErrorSo = sollya_lib_abs(minErrorSo)
2054
    sollya_lib_clear_obj(maxErrorSo)
2055
    sollya_lib_clear_obj(minErrorSo)
2056
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
2057
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
2058
    changeVarExpSo = sollya_lib_build_function_sub(\
2059
                       sollya_lib_build_function_free_variable(),\
2060
                       sollya_lib_copy_obj(intervalCenterSo))
2061
    polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpSo)
2062
    sollya_lib_clear_obj(polySo) 
2063
    sollya_lib_clear_obj(changeVarExpSo)
2064
    # If changed, reset the Sollya working precision.
2065
    if sollyaPrecChanged:
2066
        sollya_lib_set_prec(initialSollyaPrecSo)
2067
    sollya_lib_clear_obj(initialSollyaPrecSo)
2068
    if errorTypeIsNone:
2069
        sollya_lib_clear_obj(errorTypeSo)
2070
    sollya_lib_clear_obj(taylorFormSo)
2071
    # Do not clear maxErrorSo.
2072
    if absMaxErrorSa > absMinErrorSa:
2073
        sollya_lib_clear_obj(absMinErrorSo)
2074
        return((polyVarChangedSo, intervalCenterSo, absMaxErrorSo))
2075
    else:
2076
        sollya_lib_clear_obj(absMaxErrorSo)
2077
        return((polyVarChangedSo, intervalCenterSo, absMinErrorSo))
2078
# end pobyso_taylor_expansion_with_change_var_so_so
2079

    
2080
def pobyso_taylor(function, degree, point):
2081
    """ Legacy function. See pobysoTaylor_so_so. """
2082
    return(pobyso_taylor_so_so(function, degree, point))
2083

    
2084
def pobyso_taylor_so_so(functionSo, degreeSo, pointSo):
2085
    return(sollya_lib_taylor(functionSo, degreeSo, pointSo))
2086
    
2087
def pobyso_taylorform(function, degree, point = None, 
2088
                      interval = None, errorType=None):
2089
    """ Legacy function. See pobyso_taylorform_sa_sa;"""
2090
    
2091
def pobyso_taylorform_sa_sa(functionSa, \
2092
                            degreeSa, \
2093
                            pointSa, \
2094
                            intervalSa=None, \
2095
                            errorTypeSa=None, \
2096
                            precisionSa=None):
2097
    """
2098
    Compute the Taylor form of 'degreeSa' for 'functionSa' at 'pointSa' 
2099
    for 'intervalSa' with 'errorTypeSa' (a string) using 'precisionSa'. 
2100
    point: must be a Real or a Real interval.
2101
    return the Taylor form as an array
2102
    TODO: take care of the interval and of the point when it is an interval;
2103
          when errorType is not None;
2104
          take care of the other elements of the Taylor form (coefficients 
2105
          errors and delta.
2106
    """
2107
    # Absolute as the default error.
2108
    if errorTypeSa is None:
2109
        errorTypeSo = sollya_lib_absolute()
2110
    elif errorTypeSa == "relative":
2111
        errorTypeSo = sollya_lib_relative()
2112
    elif errortypeSa == "absolute":
2113
        errorTypeSo = sollya_lib_absolute()
2114
    else:
2115
        # No clean up needed.
2116
        return None
2117
    # Global precision stuff
2118
    sollyaPrecisionChangedSa = False
2119
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
2120
    if precisionSa is None:
2121
        precSa = initialSollyaPrecSa
2122
    else:
2123
        if precSa > initialSollyaPrecSa:
2124
            if precSa <= 2:
2125
                print inspect.stack()[0][3], ":precision change <= 2 requested."
2126
            pobyso_set_prec_sa_so(precSa)
2127
            sollyaPrecisionChangedSa = True
2128
    #        
2129
    if len(functionSa.variables()) > 0:
2130
        varSa = functionSa.variables()[0]
2131
        pobyso_name_free_variable_sa_so(str(varSa))
2132
    # In any case (point or interval) the parent of pointSa has a precision
2133
    # method.
2134
    pointPrecSa = pointSa.parent().precision()
2135
    if precSa > pointPrecSa:
2136
        pointPrecSa = precSa
2137
    # In any case (point or interval) pointSa has a base_ring() method.
2138
    pointBaseRingString = str(pointSa.base_ring())
2139
    if re.search('Interval', pointBaseRingString) is None: # Point
2140
        pointSo = pobyso_constant_sa_so(pointSa, pointPrecSa)
2141
    else: # Interval.
2142
        pointSo = pobyso_interval_to_range_sa_so(pointSa, pointPrecSa)
2143
    # Sollyafy the function.
2144
    functionSo = pobyso_parse_string_sa_so(functionSa._assume_str().replace('_SAGE_VAR_', ''))
2145
    if sollya_lib_obj_is_error(functionSo):
2146
        print "pobyso_tailorform: function string can't be parsed!"
2147
        return None
2148
    # Sollyafy the degree
2149
    degreeSo = sollya_lib_constant_from_int(int(degreeSa))
2150
    # Sollyafy the point
2151
    # Call Sollya
2152
    taylorFormSo = \
2153
        sollya_lib_taylorform(functionSo, degreeSo, pointSo, errorTypeSo,\
2154
                                         None)
2155
    sollya_lib_clear_obj(functionSo)
2156
    sollya_lib_clear_obj(degreeSo)
2157
    sollya_lib_clear_obj(pointSo)
2158
    sollya_lib_clear_obj(errorTypeSo)
2159
    (tfsAsList, numElements, isEndElliptic) = \
2160
            pobyso_get_list_elements_so_so(taylorFormSo)
2161
    polySo = tfsAsList[0]
2162
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
2163
    polyRealField = RealField(maxPrecision)
2164
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, polyRealField)
2165
    if sollyaPrecisionChangedSa:
2166
        sollya_lib_set_prec(initialSollyaPrecSo)
2167
    sollya_lib_clear_obj(initialSollyaPrecSo)
2168
    polynomialRing = polyRealField[str(varSa)]
2169
    polySa = polynomial(expSa, polynomialRing)
2170
    taylorFormSa = [polySa]
2171
    # Final clean-up.
2172
    sollya_lib_clear_obj(taylorFormSo)
2173
    return(taylorFormSa)
2174
# End pobyso_taylor_form_sa_sa
2175

    
2176
def pobyso_taylorform_so_so(functionSo, degreeSo, pointSo, intervalSo=None, \
2177
                            errorTypeSo=None):
2178
    createdErrorType = False
2179
    if errorTypeSo is None:
2180
        errorTypeSo = sollya_lib_absolute()
2181
        createdErrorType = True
2182
    else:
2183
        #TODO: deal with the other case.
2184
        pass
2185
    if intervalSo is None:
2186
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
2187
                                         errorTypeSo, None)
2188
    else:
2189
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
2190
                                         intervalSo, errorTypeSo, None)
2191
    if createdErrorType:
2192
        sollya_lib_clear_obj(errorTypeSo)
2193
    return resultSo
2194
        
2195

    
2196
def pobyso_univar_polynomial_print_reverse(polySa):
2197
    """ Legacy function. See pobyso_univar_polynomial_print_reverse_sa_sa. """
2198
    return(pobyso_univar_polynomial_print_reverse_sa_sa(polySa))
2199

    
2200
def pobyso_univar_polynomial_print_reverse_sa_sa(polySa):
2201
    """
2202
    Return the string representation of a univariate polynomial with
2203
    monomials ordered in the x^0..x^n order of the monomials.
2204
    Remember: Sage
2205
    """
2206
    polynomialRing = polySa.base_ring()
2207
    # A very expensive solution:
2208
    # -create a fake multivariate polynomial field with only one variable,
2209
    #   specifying a negative lexicographical order;
2210
    mpolynomialRing = PolynomialRing(polynomialRing.base(), \
2211
                                     polynomialRing.variable_name(), \
2212
                                     1, order='neglex')
2213
    # - convert the univariate argument polynomial into a multivariate
2214
    #   version;
2215
    p = mpolynomialRing(polySa)
2216
    # - return the string representation of the converted form.
2217
    # There is no simple str() method defined for p's class.
2218
    return(p.__str__())
2219
#
2220
#print pobyso_get_prec()  
2221
pobyso_set_prec(165)
2222
#print pobyso_get_prec()  
2223
#a=100
2224
#print type(a)
2225
#id(a)
2226
#print "Max arity: ", pobyso_max_arity
2227
#print "Function tripleDouble (43) as a string: ", pobyso_function_type_as_string(43)
2228
#print "Function None (44) as a string: ", pobyso_function_type_as_string(44)
2229
sys.stderr.write("\t...Pobyso check done.\n")