Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 281

Historique | Voir | Annoter | Télécharger (92,33 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
        ## When called with an empty list, sollya_lib_build_end_elliptic_list,
146
        #  produces "error".
147
        return sollya_lib_build_list(None)
148
    index = 0
149
    ## One can not append elements to an elliptic list, prepend only is 
150
    #  permitted.
151
    for argument in reversed(args):
152
        if index == 0:
153
            listSo = sollya_lib_build_end_elliptic_list(argument, None)
154
        else:
155
            listSo = sollya_lib_prepend(argument, listSo)
156
        index += 1
157
    return listSo
158
        
159
# End pobyso_build_end_elliptic_list_so_so
160

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

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

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

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

    
211
def pobyso_clear_full_list_elements_sa_so(objectListSaSo):
212
    """
213
    Clear the elements of list created by the 
214
    pobyso_get_list_elements_so_so function.
215
    objectListSaSo is as follows:
216
    - objectListSaSo[0]: a list of Sollya objects;
217
    - objectListSaSo[1]: the number of elements of the previous list;
218
    - objectListSaSo[2]: an integer that if != 0 states that the list is 
219
                         end-elliptic
220
    The objects to clear are the elements of the objectListSaSo[0] list. 
221
    """
222
    for index in xrange(0, objectListSaSo[1]):
223
        sollya_lib_clear_obj(objectListSaSo[0][index])
224
# End pobyso_clear_full_list_elements_sa_so
225

    
226
def pobyso_clear_list_elements_sa_so(objectListSaSo):
227
    """
228
    Clear the elements of list of references to Sollya objects
229
    """
230
    for index in xrange(0, len(objectListSaSo)):
231
        sollya_lib_clear_obj(objectListSaSo[index])
232
# End pobyso_clear_list_elements_sa_so
233

    
234
def pobyso_clear_obj(objSo):
235
    """
236
    Free a Sollya object's memory.
237
    Very thin wrapper around sollya_lib_clear_obj().
238
    """
239
    sollya_lib_clear_obj(objSo)
240
# End pobyso_clear_obj. 
241

    
242
def pobyso_clear_taylorform_sa_so(taylorFormSaSo):
243
    """
244
    This method is rapper around pobyso_clear_list_elements_sa_so.
245
    It is a legacy method left here since it may be used in existing code 
246
    where Taylor forms are used as Sage lists obtained by converting
247
    Sollya Taylor forms (a list made of:
248
    - a polynomial;
249
    - a list of intervals enclosing the errors accumulated when computing
250
      the polynomial coefficients;
251
    - a bound on the approximation error between the polynomial and the 
252
      function.)
253
    A Taylor form directly obtained from pobyso_taylorform_so_so is cleared
254
    by sollya_lib_clear_obj. 
255
    """
256
    pobyso_clear_list_elements_sa_so(taylorFormSaSo)
257
# End pobyso_clear_taylorform_sa_so 
258

    
259
def pobyso_cmp(rnArgSa, cteSo):
260
    """
261
    Deprecated, use pobyso_cmp_sa_so_sa instead.
262
    """
263
    print "Deprecated, use pobyso_cmp_sa_so_sa instead."
264
    return pobyso_cmp_sa_so_sa(rnArgSa, cteSo)
265
# End  pobyso_cmp
266
    
267
def pobyso_cmp_sa_so_sa(rnArgSa, cteSo):
268
    """
269
    Compare the MPFR value a RealNumber with that of a Sollya constant.
270
    
271
    Get the value of the Sollya constant into a RealNumber and compare
272
    using MPFR. Could be optimized by working directly with a mpfr_t
273
    for the intermediate number. 
274
    """
275
    # Get the precision of the Sollya constant to build a Sage RealNumber 
276
    # with enough precision.to hold it.
277
    precisionOfCte = c_int(0)
278
    # From the Sollya constant, create a local Sage RealNumber.
279
    sollya_lib_get_prec_of_constant(precisionOfCte, cteSo) 
280
    #print "Precision of constant: ", precisionOfCte
281
    RRRR = RealField(precisionOfCte.value)
282
    rnLocalSa = RRRR(0)
283
    sollya_lib_get_constant(get_rn_value(rnLocalSa), cteSo)
284
    #
285
    ## Compare the Sage RealNumber version of the Sollya constant with rnArg
286
    #  through a direct comparison of underlying MPFR numbers.
287
    return cmp_rn_value(rnArgSa, rnLocal)
288
# End pobyso_smp_sa_so_sa
289

    
290
def pobyso_compute_pos_function_abs_val_bounds_sa_sa(funcSa, lowerBoundSa, \
291
                                                     upperBoundSa):
292
    """
293
    TODO: completely rework and test.
294
    """
295
    pobyso = pobyso_name_free_variable_sa_so(funcSa.variables()[0])
296
    funcSo = pobyso_parse_string(funcSa._assume_str().replace('_SAGE_VAR_', ''))
297
    rangeSo = pobyso_range_sa_so(lowerBoundSa, upperBoundSa)
298
    infnormSo = pobyso_infnorm_so_so(funcSo,rangeSo)
299
    # Sollya return the infnorm as an interval.
300
    fMaxSa = pobyso_get_interval_from_range_so_sa(infnormSo)
301
    # Get the top bound and compute the binade top limit.
302
    fMaxUpperBoundSa = fMaxSa.upper()
303
    binadeTopLimitSa = 2**ceil(fMaxUpperBoundSa.log2())
304
    # Put up together the function to use to compute the lower bound.
305
    funcAuxSo = pobyso_parse_string(str(binadeTopLimitSa) +  \
306
                                    '-(' + f._assume_str().replace('_SAGE_VAR_', '') + ')')
307
    pobyso_autoprint(funcAuxSo)
308
    # Clear the Sollya range before a new call to infnorm and issue the call.
309
    sollya_lib_clear_obj(infnormSo)
310
    infnormSo = pobyso_infnorm_so_so(funcAuxSo,rangeSo)
311
    fMinSa = pobyso_get_interval_from_range_so_sa(infnormSo)
312
    sollya_lib_clear_obj(infnormSo)
313
    fMinLowerBoundSa = binadeTopLimitSa - fMinSa.lower()
314
    # Compute the maximum of the precisions of the different bounds.
315
    maxPrecSa = max([fMinLowerBoundSa.parent().precision(), \
316
                     fMaxUpperBoundSa.parent().precision()])
317
    # Create a RealIntervalField and create an interval with the "good" bounds.
318
    RRRI = RealIntervalField(maxPrecSa)
319
    imageIntervalSa = RRRI(fMinLowerBoundSa, fMaxUpperBoundSa)
320
    # Free the unneeded Sollya objects
321
    sollya_lib_clear_obj(funcSo)
322
    sollya_lib_clear_obj(funcAuxSo)
323
    sollya_lib_clear_obj(rangeSo)
324
    return(imageIntervalSa)
325
# End pobyso_compute_pos_function_abs_val_bounds_sa_sa
326

    
327
def pobyso_compute_precision_decay_ratio_function_sa_so():
328
    """
329
    Compute the precision decay ratio function for polynomial 
330
    coefficient progressive trucation.
331
    """
332
    functionText = """
333
    proc(deg, a, b, we, wq) 
334
    {
335
      k = we * (exp(x/a)-1) + wq * (b*x)^2 + (1-we-wq) * x;
336
      return k/k(d);
337
    };
338
    """
339
    return pobyso_parse_string_sa_so(functionText)
340
# End  pobyso_compute_precision_decay_ratio_function.
341

    
342

    
343
def pobyso_constant(rnArg):
344
    """ Legacy function. See pobyso_constant_sa_so. """
345
    return(pobyso_constant_sa_so(rnArg))
346
    
347
def pobyso_constant_sa_so(rnArgSa, precisionSa=None):
348
    """
349
    Create a Sollya constant from a Sage RealNumber.
350
    The sollya_lib_constant() function creates a constant
351
    with the same precision as the source.
352
    """
353
    ## Precision stuff. If one wants to change precisions,
354
    #  everything takes place in Sage. That only makes
355
    #  sense if one wants to reduce the precision.
356
    # TODO: revisit precision stuff with new technique.
357
    if not precisionSa is None:
358
        RRR = RealField(precisionSa)
359
        rnArgSa = RRR(rnArgSa)
360
    #print rnArgSa, rnArgSa.precision()
361
    # Sollya constant creation takes place here.
362
    return sollya_lib_constant(get_rn_value(rnArgSa))
363
# End pobyso_constant_sa_so
364
     
365
def pobyso_constant_0_sa_so():
366
    """
367
    Obvious.
368
    """
369
    return pobyso_constant_from_int_sa_so(0)
370

    
371
def pobyso_constant_1():
372
    """
373
    Obvious.
374
    Legacy function. See pobyso_constant_so_so. 
375
    """
376
    return pobyso_constant_1_sa_so()
377

    
378
def pobyso_constant_1_sa_so():
379
    """
380
    Obvious.
381
    """
382
    return(pobyso_constant_from_int_sa_so(1))
383

    
384
def pobyso_constant_from_int(anInt):
385
    """ Legacy function. See pobyso_constant_from_int_sa_so. """
386
    return pobyso_constant_from_int_sa_so(anInt)
387

    
388
def pobyso_constant_from_int_sa_so(anInt):
389
    """
390
    Get a Sollya constant from a Sage int.
391
    """
392
    return sollya_lib_constant_from_int64(long(anInt))
393

    
394
def pobyso_constant_from_int_so_sa(constSo):
395
    """
396
    Get a Sage int from a Sollya int constant.
397
    Usefull for precision or powers in polynomials.
398
    """
399
    constSa = c_long(0)
400
    sollya_lib_get_constant_as_int64(byref(constSa), constSo)
401
    return constSa.value
402
# End pobyso_constant_from_int_so_sa
403

    
404
def pobyso_constant_from_mpq_sa_so(rationalSa):
405
    """
406
    Make a Sollya constant from Sage rational.
407
    The Sollya constant is an unevaluated expression.
408
    Hence no precision argument is needed.
409
    It is better to leave this way since Sollya has its own
410
    optimized evaluation mecanism that tries very hard to
411
    return exact values or at least faithful ones.
412
    """
413
    ratExprSo = \
414
        sollya_lib_constant_from_mpq(sgmp_get_rational_value(rationalSa))
415
    return ratExprSo
416
# End pobyso_constant_from_mpq_sa_so.
417

    
418
def pobyso_constant_sollya_prec_sa_so(rnArgSa):
419
    """
420
    Create a Sollya constant from a Sage RealNumber at the
421
    current precision in Sollya.
422
    """
423
    currentSollyaPrecSa = pobyso_get_prec_so_sa()
424
    return pobyso_constant_sa_so(rnArgSa, currentSollyaPrecSa)
425
# End pobyso_constant_sollya_prec_sa_so
426

    
427
def pobyso_end_elliptic_list_so_sa_so(objectsListSo, intCountSa):
428
    """
429
    Create a Sollya end elliptic list made of the objectListSo[0] to
430
     objectsListSo[intCountSa-1] objects.
431
    """     
432
    return sollya_lib_end_elliptic_list(objectSo, int(intCountSa))
433

    
434
def pobyso_error_so():
435
    return sollya_lib_error(None)
436
# End pobyso_error().
437

    
438
def pobyso_evaluate_so_sa(funcSo, argumentSo):
439
    """
440
    Evaluates funcSo for arguemntSo through sollya_lib_evaluate() and return
441
    the result as a Sage object
442
    """
443
    evalSo = sollya_lib_evaluate(funcSo, argumentSo)
444
    if pobyso_is_error_so_sa(evalSo):
445
        return None
446
    if pobyso_is_range_so_sa(evalSo):
447
        retVal = pobyso_range_to_interval_so_sa(evalSo)
448
    else:
449
        retVal = pobyso_get_constant_as_rn(evalSo)
450
    sollya_lib_clear_obj(evalSo)
451
    return retVal
452
# End pobyso_evaluate_so_sa.
453

    
454
def pobyso_evaluate_so_so(funcSo, argumentSo):
455
    """
456
    Evaluates funcSo for arguemntSo through sollya_lib_evaluate().
457
    """
458
    return sollya_lib_evaluate(funcSo, argumentSo)
459
# End pobyso_evaluate_so_so.
460
#
461
def pobyso_diff_so_so(funcSo):
462
    """
463
    Very thin wrapper around sollya_lib_diff.
464
    """
465
    ## TODO: add a check to make sure funcSo is a functional expression.
466
    return sollya_lib_diff(funcSo)
467

    
468
def pobyso_dirty_find_zeros_so_so(funcSo, rangeSo):
469
    """
470
    Thin wrapper over sollya_lib_dirtyfindzeros()
471
    """
472
    return sollya_lib_dirtyfindzeros(funcSo, rangeSo)
473
# End pobys_dirty_find_zeros
474

    
475
def pobyso_dirty_inf_norm_so_so(funcSo, rangeSo, preSo=None):
476
    """
477
    Thin wrapper around sollya_dirtyinfnorm().
478
    """
479
    # TODO: manage the precision change.
480
    
481
    return sollya_lib_dirtyinfnorm(funcSo, rangeSo)
482
# End pobyso_dirty_inf_norm_so_so
483

    
484
def pobyso_find_zeros_so_so(funcSo, rangeSo):
485
    """
486
    Thin wrapper over sollya_lib_findzeros()
487
    """
488
    return sollya_lib_findzeros(funcSo, rangeSo)
489
# End pobys_find_zeros
490

    
491
def pobyso_float_list_so_sa(listSo):
492
    """
493
    Return a Sollya list of floating-point numbers as a Sage list of 
494
    floating-point numbers.
495
    TODO: add a test to make sure that each element of the list is a constant.
496
    """
497
    listSa   = []
498
    ## The function returns none if the list is empty or an error has happened.
499
    retVal = pobyso_get_list_elements_so_so(listSo)
500
    if retVal is None:
501
        return listSa
502
    ## Just in case the interface is changed and an empty list is returned 
503
    #  instead of None.
504
    elif len(retVal) == 0:
505
        return listSa
506
    else:
507
        ## Remember pobyso_get_list_elements_so_so returns more information
508
        #  than just the elements of the list (# elements, is_elliptic)
509
        listSaSo, numElements, isEndElliptic = retVal
510
    ## Return an empty list.
511
    if numElements == 0:
512
        return listSa
513
    ## Search first for the maximum precision of the elements
514
    maxPrecSa = 0
515
    for floatSo in listSaSo:
516
        #pobyso_autoprint(floatSo)
517
        curPrecSa =  pobyso_get_prec_of_constant_so_sa(floatSo)
518
        if curPrecSa > maxPrecSa:
519
            maxPrecSa = curPrecSa
520
    ##
521
    RF = RealField(maxPrecSa)
522
    ##
523
    for floatSo in listSaSo:
524
        listSa.append(pobyso_get_constant_as_rn_with_rf_so_sa(floatSo))
525
    return listSa
526
# End pobyso_float_list_so_sa
527

    
528
def pobyso_float_poly_sa_so(polySa, precSa = None):
529
    """
530
    Create a Sollya polynomial from a Sage RealField polynomial.
531
    """
532
    ## TODO: filter arguments.
533
    ## Precision. If a precision is given, convert the polynomial
534
    #  into the right polynomial field. If not convert it straight
535
    #  to Sollya.
536
    sollyaPrecChanged = False 
537
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
538
    if precSa is None:
539
        precSa = polySa.parent().base_ring().precision()
540
    if (precSa > initialSollyaPrecSa):
541
        assert precSa >= 2, "Precision change <2 requested"
542
        if precSa <= 2:
543
            print inspect.stack()[0][3], ": precision change <= 2 requested"
544
        precSo = pobyso_constant_from_int(precSa)
545
        pobyso_set_prec_so_so(precSo)
546
        sollya_lib_clear_obj(precSo)
547
        sollyaPrecChanged = True
548
    ## Free variable stuff.
549
    freeVariableNameChanged = False 
550
    polyFreeVariableNameSa = \
551
        str(polySa.variables()[0])
552
    currentFreeVariableNameSa = pobyso_get_free_variable_name_so_sa() 
553
    if polyFreeVariableNameSa != currentFreeVariableNameSa:
554
        #print "Free variable names do not match.", polyFreeVariableNameSa
555
        sollya_lib_name_free_variable(polyFreeVariableNameSa)
556
        freeVariableNameChanged = True
557
    ## Get exponents and coefficients.
558
    exponentsSa     = polySa.exponents()
559
    coefficientsSa  = polySa.coefficients()
560
    ## Build the polynomial.
561
    polySo = None
562
    for coefficientSa, exponentSa in zip(coefficientsSa, exponentsSa):
563
        #print coefficientSa.n(prec=precSa), exponentSa
564
        coefficientSo = \
565
            pobyso_constant_sa_so(coefficientSa)
566
        #pobyso_autoprint(coefficientSo)
567
        exponentSo = \
568
            pobyso_constant_from_int_sa_so(exponentSa)
569
        #pobyso_autoprint(exponentSo)
570
        monomialSo = sollya_lib_build_function_pow(
571
                       sollya_lib_build_function_free_variable(),
572
                       exponentSo)
573
        polyTermSo = sollya_lib_build_function_mul(coefficientSo,
574
                                                       monomialSo)
575
        if polySo is None:
576
            polySo = polyTermSo
577
        else:
578
            polySo = sollya_lib_build_function_add(polySo, polyTermSo)
579
    if sollyaPrecChanged:
580
        pobyso_set_prec_so_so(initialSollyaPrecSo)
581
    sollya_lib_clear_obj(initialSollyaPrecSo)
582
    ## Do not set back the free variable name in Sollya to its initial value: 
583
    #  it will change it back, in the Sollya polynomial, to what it was in the 
584
    #  first place.
585
    return polySo
586
# End pobyso_float_poly_sa_so    
587

    
588
def pobyso_float_poly_so_sa(polySo, realFieldSa=None):
589
    """
590
    Convert a Sollya polynomial into a Sage floating-point polynomial.
591
    If no realField is given, a RealField corresponding to the maximum 
592
    precision of the coefficients is internally computed.
593
    The real field is not returned but can be easily retrieved from 
594
    the polynomial itself.
595
    ALGORITHM:
596
    - (optional) compute the RealField of the coefficients;
597
    - convert the Sollya expression into a Sage expression;
598
    - convert the Sage expression into a Sage polynomial
599
    """    
600
    if realFieldSa is None:
601
        expressionPrecSa = pobyso_get_max_prec_of_exp_so_sa(polySo)
602
        #print "Maximum precision of Sollya polynomial coefficients:", expressionPrecSa
603
        if expressionPrecSa < 2 or expressionPrecSa > 2147483391:
604
            print "Maximum degree of expression:", expressionPrecSa
605
        realFieldSa      = RealField(expressionPrecSa)
606
    #print "Sollya expression before...",
607
    #pobyso_autoprint(polySo)
608

    
609
    expressionSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo,
610
                                                             realFieldSa)
611
    #print "...Sollya expression after."
612
    #pobyso_autoprint(polySo)
613
    polyVariableSa = expressionSa.variables()[0]
614
    polyRingSa     = realFieldSa[str(polyVariableSa)]
615
    #print polyRingSa
616
    # Do not use the polynomial(expressionSa, ring=polyRingSa) form!
617
    polynomialSa = polyRingSa(expressionSa)
618
    polyCoeffsListSa = polynomialSa.coefficients()
619
    #for coeff in polyCoeffsListSa:
620
    #    print coeff.abs().n()
621
    return polynomialSa
622
# End pobyso_float_poly_so_sa
623

    
624
def pobyso_free_variable():
625
    """
626
    Ultra thin wrapper around the sollya_lib_function_build_free_variable function.
627
    """
628
    return sollya_lib_build_function_free_variable()
629
    
630
def pobyso_function_type_as_string(funcType):
631
    """ Legacy function. See pobyso_function_type_as_string_so_sa. """
632
    return(pobyso_function_type_as_string_so_sa(funcType))
633

    
634
def pobyso_function_type_as_string_so_sa(funcType):
635
    """
636
    Numeric Sollya function codes -> Sage mathematical function names.
637
    Notice that pow -> ^ (a la Sage, not a la Python).
638
    """
639
    if funcType == SOLLYA_BASE_FUNC_ABS:
640
        return "abs"
641
    elif funcType == SOLLYA_BASE_FUNC_ACOS:
642
        return "arccos"
643
    elif funcType == SOLLYA_BASE_FUNC_ACOSH:
644
        return "arccosh"
645
    elif funcType == SOLLYA_BASE_FUNC_ADD:
646
        return "+"
647
    elif funcType == SOLLYA_BASE_FUNC_ASIN:
648
        return "arcsin"
649
    elif funcType == SOLLYA_BASE_FUNC_ASINH:
650
        return "arcsinh"
651
    elif funcType == SOLLYA_BASE_FUNC_ATAN:
652
        return "arctan"
653
    elif funcType == SOLLYA_BASE_FUNC_ATANH:
654
        return "arctanh"
655
    elif funcType == SOLLYA_BASE_FUNC_CEIL:
656
        return "ceil"
657
    elif funcType == SOLLYA_BASE_FUNC_CONSTANT:
658
        return "cte"
659
    elif funcType == SOLLYA_BASE_FUNC_COS:
660
        return "cos"
661
    elif funcType == SOLLYA_BASE_FUNC_COSH:
662
        return "cosh"
663
    elif funcType == SOLLYA_BASE_FUNC_DIV:
664
        return "/"
665
    elif funcType == SOLLYA_BASE_FUNC_DOUBLE:
666
        return "double"
667
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEDOUBLE:
668
        return "doubleDouble"
669
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEEXTENDED:
670
        return "doubleDxtended"
671
    elif funcType == SOLLYA_BASE_FUNC_ERF:
672
        return "erf"
673
    elif funcType == SOLLYA_BASE_FUNC_ERFC:
674
        return "erfc"
675
    elif funcType == SOLLYA_BASE_FUNC_EXP:
676
        return "exp"
677
    elif funcType == SOLLYA_BASE_FUNC_EXP_M1:
678
        return "expm1"
679
    elif funcType == SOLLYA_BASE_FUNC_FLOOR:
680
        return "floor"
681
    elif funcType == SOLLYA_BASE_FUNC_FREE_VARIABLE:
682
        return "freeVariable"
683
    elif funcType == SOLLYA_BASE_FUNC_HALFPRECISION:
684
        return "halfPrecision"
685
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYCONSTANT:
686
        return "libraryConstant"
687
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYFUNCTION:
688
        return "libraryFunction"
689
    elif funcType == SOLLYA_BASE_FUNC_LOG:
690
        return "log"
691
    elif funcType == SOLLYA_BASE_FUNC_LOG_10:
692
        return "log10"
693
    elif funcType == SOLLYA_BASE_FUNC_LOG_1P:
694
        return "log1p"
695
    elif funcType == SOLLYA_BASE_FUNC_LOG_2:
696
        return "log2"
697
    elif funcType == SOLLYA_BASE_FUNC_MUL:
698
        return "*"
699
    elif funcType == SOLLYA_BASE_FUNC_NEARESTINT:
700
        return "round"
701
    elif funcType == SOLLYA_BASE_FUNC_NEG:
702
        return "__neg__"
703
    elif funcType == SOLLYA_BASE_FUNC_PI:
704
        return "pi"
705
    elif funcType == SOLLYA_BASE_FUNC_POW:
706
        return "^"
707
    elif funcType == SOLLYA_BASE_FUNC_PROCEDUREFUNCTION:
708
        return "procedureFunction"
709
    elif funcType == SOLLYA_BASE_FUNC_QUAD:
710
        return "quad"
711
    elif funcType == SOLLYA_BASE_FUNC_SIN:
712
        return "sin"
713
    elif funcType == SOLLYA_BASE_FUNC_SINGLE:
714
        return "single"
715
    elif funcType == SOLLYA_BASE_FUNC_SINH:
716
        return "sinh"
717
    elif funcType == SOLLYA_BASE_FUNC_SQRT:
718
        return "sqrt"
719
    elif funcType == SOLLYA_BASE_FUNC_SUB:
720
        return "-"
721
    elif funcType == SOLLYA_BASE_FUNC_TAN:
722
        return "tan"
723
    elif funcType == SOLLYA_BASE_FUNC_TANH:
724
        return "tanh"
725
    elif funcType == SOLLYA_BASE_FUNC_TRIPLEDOUBLE:
726
        return "tripleDouble"
727
    else:
728
        return None
729

    
730
def pobyso_get_constant(rnArgSa, constSo):
731
    """ Legacy function. See pobyso_get_constant_so_sa. """
732
    return pobyso_get_constant_so_sa(rnArgSa, constSo)
733
# End pobyso_get_constant
734

    
735
def pobyso_get_constant_so_sa(rnArgSa, constSo):
736
    """
737
    Set the value of rnArgSo to the value of constSo in MPFR_RNDN mode.
738
    rnArg must already exist and belong to some RealField.
739
    We assume that constSo points to a Sollya constant.
740
    """
741
    outcome = sollya_lib_get_constant(get_rn_value(rnArgSa), constSo)
742
    if outcome == 0: # Failure because constSo is not a constant expression.
743
        return None
744
    else:
745
        return outcome
746
# End  pobyso_get_constant_so_sa
747
   
748
def pobyso_get_constant_as_rn(ctExpSo):
749
    """ 
750
    Legacy function. See pobyso_get_constant_as_rn_so_sa. 
751
    """ 
752
    return(pobyso_get_constant_as_rn_so_sa(ctExpSo))
753
    
754
def pobyso_get_constant_as_rn_so_sa(constExpSo):
755
    """
756
    Get a Sollya constant as a Sage "real number".
757
    The precision of the floating-point number returned is that of the Sollya
758
    constant.
759
    """
760
    #print "Before computing precision of variable..."
761
    #pobyso_autoprint(constExpSo)
762
    precisionSa  = pobyso_get_prec_of_constant_so_sa(constExpSo)
763
    #print "precisionSa:", precisionSa
764
    ## If the expression can not be exactly converted, None is returned.
765
    #  In this case opt for the Sollya current expression.
766
    if precisionSa is None:
767
        precisionSa = pobyso_get_prec_so_sa()
768
    RRRR = RealField(precisionSa)
769
    rnSa = RRRR(0)
770
    outcome = sollya_lib_get_constant(get_rn_value(rnSa), constExpSo)
771
    if outcome == 0:
772
        return None
773
    else:
774
        return rnSa
775
# End pobyso_get_constant_as_rn_so_sa
776

    
777
def pobyso_get_constant_as_rn_with_rf(ctExp, realField):
778
    """ 
779
    Legacy function. See pobyso_get_constant_as_rn_with_rf_so_sa.
780
    """
781
    return pobyso_get_constant_as_rn_with_rf_so_sa(ctExp, realField)
782
# End pobyso_get_constant_as_rn_with_rf
783
    
784
def pobyso_get_constant_as_rn_with_rf_so_sa(ctExpSo, realFieldSa = None):
785
    """
786
    Get a Sollya constant as a Sage "real number".
787
    If no real field is specified, the precision of the floating-point number 
788
    returned is that of the Sollya constant.
789
    Otherwise is is that of the real field. Hence rounding may happen.
790
    """
791
    if realFieldSa is None:
792
        return pobyso_get_constant_as_rn_so_sa(ctExpSo)
793
    rnSa = realFieldSa(0)
794
    outcome = sollya_lib_get_constant(get_rn_value(rnSa), ctExpSo)
795
    if outcome == 0:
796
        return None
797
    else:
798
        return rnSa
799
# End pobyso_get_constant_as_rn_with_rf_so_sa
800

    
801
def pobyso_get_free_variable_name():
802
    """ 
803
    Legacy function. See pobyso_get_free_variable_name_so_sa.
804
    """
805
    return(pobyso_get_free_variable_name_so_sa())
806

    
807
def pobyso_get_free_variable_name_so_sa():
808
    return sollya_lib_get_free_variable_name()
809
    
810
def pobyso_get_function_arity(expressionSo):
811
    """ 
812
    Legacy function. See pobyso_get_function_arity_so_sa.
813
    """
814
    return(pobyso_get_function_arity_so_sa(expressionSo))
815

    
816
def pobyso_get_function_arity_so_sa(expressionSo):
817
    arity = c_int(0)
818
    sollya_lib_get_function_arity(byref(arity),expressionSo)
819
    return int(arity.value)
820

    
821
def pobyso_get_head_function(expressionSo):
822
    """ 
823
    Legacy function. See pobyso_get_head_function_so_sa. 
824
    """
825
    return(pobyso_get_head_function_so_sa(expressionSo)) 
826

    
827
def pobyso_get_head_function_so_sa(expressionSo):
828
    functionType = c_int(0)
829
    sollya_lib_get_head_function(byref(functionType), expressionSo)
830
    return int(functionType.value)
831

    
832
def pobyso_get_interval_from_range_so_sa(soRange, realIntervalFieldSa = None ):
833
    """
834
    Return the Sage interval corresponding to the Sollya range argument.
835
    If no reaIntervalField is passed as an argument, the interval bounds are not
836
    rounded: they are elements of RealIntervalField of the "right" precision
837
    to hold all the digits.
838
    """
839
    prec = c_int(0)
840
    if realIntervalFieldSa is None:
841
        retval = sollya_lib_get_prec_of_range(byref(prec), soRange, None)
842
        if retval == 0:
843
            return None
844
        realIntervalFieldSa = RealIntervalField(prec.value)
845
    intervalSa = realIntervalFieldSa(0,0)
846
    retval = \
847
        sollya_lib_get_interval_from_range(get_interval_value(intervalSa),\
848
                                           soRange)
849
    if retval == 0:
850
        return None
851
    return intervalSa
852
# End pobyso_get_interval_from_range_so_sa
853

    
854
def pobyso_get_list_elements(soObj):
855
    """ Legacy function. See pobyso_get_list_elements_so_so. """
856
    return pobyso_get_list_elements_so_so(soObj)
857
 
858
def pobyso_get_list_elements_so_so(objectListSo):
859
    """
860
    Get the Sollya list elements as a Sage/Python array of Sollya objects.
861
    
862
    INPUT:
863
    - objectListSo: a Sollya list of Sollya objects.
864
    
865
    OUTPUT:
866
    - a Sage/Python tuple made of:
867
      - a Sage/Python list of Sollya objects,
868
      - a Sage/Python int holding the number of elements,
869
      - a Sage/Python int stating (!= 0) that the list is end-elliptic.
870
    NOTE::
871
        We recover the addresses of the Sollya object from the list of pointers
872
        returned by sollya_lib_get_list_elements. The list itself is freed.
873
    TODO::
874
        Figure out what to do with numElements since the number of elements
875
        can easily be recovered from the list itself. 
876
        Ditto for isEndElliptic.
877
    """
878
    listAddress = POINTER(c_longlong)()
879
    numElements = c_int(0)
880
    isEndElliptic = c_int(0)
881
    listAsSageList = []
882
    result = sollya_lib_get_list_elements(byref(listAddress),\
883
                                          byref(numElements),\
884
                                          byref(isEndElliptic),\
885
                                          objectListSo)
886
    if result == 0 :
887
        return None
888
    for i in xrange(0, numElements.value, 1):
889
       #listAsSageList.append(sollya_lib_copy_obj(listAddress[i]))
890
       listAsSageList.append(listAddress[i])
891
       # Clear each of the elements returned by Sollya.
892
       #sollya_lib_clear_obj(listAddress[i])
893
    # Free the list itself.   
894
    sollya_lib_free(listAddress)
895
    return (listAsSageList, numElements.value, isEndElliptic.value)
896

    
897
def pobyso_get_max_prec_of_exp(soExp):
898
    """ Legacy function. See pobyso_get_max_prec_of_exp_so_sa. """
899
    return pobyso_get_max_prec_of_exp_so_sa(soExp)
900

    
901
def pobyso_get_max_prec_of_exp_so_sa(expSo):
902
    """
903
    Get the maximum precision used for the numbers in a Sollya expression.
904
    
905
    Arguments:
906
    soExp -- a Sollya expression pointer
907
    Return value:
908
    A Python integer
909
    TODO: 
910
    - error management;
911
    - correctly deal with numerical type such as DOUBLEEXTENDED.
912
    """
913
    if expSo is None:
914
        print inspect.stack()[0][3], ": expSo is None."
915
        return 0
916
    maxPrecision = 0
917
    minConstPrec = 0
918
    currentConstPrec = 0
919
    #pobyso_autoprint(expSo)
920
    operator = pobyso_get_head_function_so_sa(expSo)
921
    if (operator != SOLLYA_BASE_FUNC_CONSTANT) and \
922
    (operator != SOLLYA_BASE_FUNC_FREE_VARIABLE):
923
        (arity, subexpressions) = pobyso_get_subfunctions_so_sa(expSo)
924
        for i in xrange(arity):
925
            maxPrecisionCandidate = \
926
                pobyso_get_max_prec_of_exp_so_sa(subexpressions[i])
927
            if maxPrecisionCandidate > maxPrecision:
928
                maxPrecision = maxPrecisionCandidate
929
        return maxPrecision
930
    elif operator == SOLLYA_BASE_FUNC_CONSTANT:
931
        #minConstPrec = pobyso_get_min_prec_of_constant_so_sa(expSo)
932
        #currentConstPrec = pobyso_get_min_prec_of_constant_so_sa(soExp)
933
        #print minConstPrec, " - ", currentConstPrec 
934
        return pobyso_get_min_prec_of_constant_so_sa(expSo)
935
    
936
    elif operator == SOLLYA_BASE_FUNC_FREE_VARIABLE:
937
        return 0
938
    else:
939
        print "pobyso_get_max_prec_of_exp_so_sa: unexepected operator."
940
        return 0
941

    
942
def pobyso_get_min_prec_of_constant_so_sa(constExpSo):
943
    """
944
    Get the minimum precision necessary to represent the value of a Sollya
945
    constant.
946
    MPFR_MIN_PREC and powers of 2 are taken into account.
947
    We assume that constExpSo is a pointer to a Sollay constant expression.
948
    """
949
    constExpAsRnSa = pobyso_get_constant_as_rn_so_sa(constExpSo)
950
    return(min_mpfr_size(get_rn_value(constExpAsRnSa)))
951

    
952
def pobyso_get_poly_so_sa(polySo, realFieldSa=None):
953
    """
954
    Convert a Sollya polynomial into a Sage polynomial.
955
    Legacy function. Use pobyso_float_poly_so_sa() instead.
956
    """    
957
    return pobyso_float_poly_so_sa(polySo,realFieldSa)
958
# End pobyso_get_poly_so_sa
959

    
960
def pobyso_get_prec():
961
    """ Legacy function. See pobyso_get_prec_so_sa(). """
962
    return pobyso_get_prec_so_sa()
963

    
964
def pobyso_get_prec_so():
965
    """
966
    Get the current default precision in Sollya.
967
    The return value is a Sollya object.
968
    Usefull when modifying the precision back and forth by avoiding
969
    extra conversions.
970
    """
971
    return sollya_lib_get_prec(None)
972
    
973
def pobyso_get_prec_so_sa():
974
    """
975
    Get the current default precision in Sollya.
976
    The return value is Sage/Python int.
977
    """
978
    precSo = sollya_lib_get_prec()
979
    precSa = pobyso_constant_from_int_so_sa(precSo)
980
    sollya_lib_clear_obj(precSo)
981
    return precSa
982
# End pobyso_get_prec_so_sa.
983

    
984
def pobyso_get_prec_so_so_sa():
985
    """
986
    Return the current precision both as a Sollya object and a
987
    Sage integer as hybrid tuple.
988
    To avoid multiple calls for precision manipulations. 
989
    """
990
    precSo = sollya_lib_get_prec()
991
    precSa = pobyso_constant_from_int_so_sa(precSo)
992
    return (precSo, int(precSa))
993
    
994
def pobyso_get_prec_of_constant(ctExpSo):
995
    """ Legacy function. See pobyso_get_prec_of_constant_so_sa. """
996
    return pobyso_get_prec_of_constant_so_sa(ctExpSo)
997

    
998
def pobyso_get_prec_of_constant_so_sa(ctExpSo):
999
    """
1000
    Tries to find a precision to represent ctExpSo without rounding.
1001
    If not possible, returns None.
1002
    """
1003
    #print "Entering pobyso_get_prec_of_constant_so_sa..."
1004
    prec = c_int(0)
1005
    retc = sollya_lib_get_prec_of_constant(byref(prec), ctExpSo, None)
1006
    if retc == 0:
1007
        #print "pobyso_get_prec_of_constant_so_sa failed."
1008
        return None
1009
    #print "...exiting pobyso_get_prec_of_constant_so_sa."
1010
    return int(prec.value)
1011

    
1012
def pobyso_get_prec_of_range_so_sa(rangeSo):
1013
    """
1014
    Returns the number of bits elements of a range are coded with.
1015
    """
1016
    prec = c_int(0)
1017
    retc = sollya_lib_get_prec_of_range(byref(prec), rangeSo, None)
1018
    if retc == 0:
1019
        return(None)
1020
    return int(prec.value)
1021
# End pobyso_get_prec_of_range_so_sa()
1022

    
1023
def pobyso_get_sage_exp_from_sollya_exp(sollyaExpSo, realField = RR):
1024
    """ Legacy function. See pobyso_get_sage_exp_from_sollya_exp_so_sa. """
1025
    return pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, 
1026
                                                     realField = RR)
1027

    
1028
def pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, realFieldSa = RR):
1029
    """
1030
    Get a Sage expression from a Sollya expression. 
1031
    Currently only tested with polynomials with floating-point coefficients.
1032
    Notice that, in the returned polynomial, the exponents are RealNumbers.
1033
    """
1034
    #pobyso_autoprint(sollyaExp)
1035
    operatorSa = pobyso_get_head_function_so_sa(sollyaExpSo)
1036
    sollyaLibFreeVariableName = sollya_lib_get_free_variable_name()
1037
    ## Get rid of the "_"'s in "_x_", if any.
1038
    sollyaLibFreeVariableName = re.sub('_', '', sollyaLibFreeVariableName)
1039
    # Constants and the free variable are special cases.
1040
    # All other operator are dealt with in the same way.
1041
    if (operatorSa != SOLLYA_BASE_FUNC_CONSTANT) and \
1042
       (operatorSa != SOLLYA_BASE_FUNC_FREE_VARIABLE):
1043
        (aritySa, subexpressionsSa) = pobyso_get_subfunctions_so_sa(sollyaExpSo)
1044
        if aritySa == 1:
1045
            sageExpSa = eval(pobyso_function_type_as_string_so_sa(operatorSa) + \
1046
            "(" + pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], \
1047
            realFieldSa) + ")")
1048
        elif aritySa == 2:
1049
            # We do not get through the preprocessor.
1050
            # The "^" operator is then a special case.
1051
            if operatorSa == SOLLYA_BASE_FUNC_POW:
1052
                operatorAsStringSa = "**"
1053
            else:
1054
                operatorAsStringSa = \
1055
                    pobyso_function_type_as_string_so_sa(operatorSa)
1056
            sageExpSa = \
1057
              eval("pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], realFieldSa)"\
1058
              + " " + operatorAsStringSa + " " + \
1059
                   "pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[1], realFieldSa)")
1060
        # We do not know yet how to deal with arity >= 3 
1061
        # (is there any in Sollya anyway?).
1062
        else:
1063
            sageExpSa = eval('None')
1064
        return sageExpSa
1065
    elif operatorSa == SOLLYA_BASE_FUNC_CONSTANT:
1066
        #print "This is a constant"
1067
        return pobyso_get_constant_as_rn_with_rf_so_sa(sollyaExpSo, realFieldSa)
1068
    elif operatorSa == SOLLYA_BASE_FUNC_FREE_VARIABLE:
1069
        #print "This is the free variable"
1070
        return eval(sollyaLibFreeVariableName)
1071
    else:
1072
        print "Unexpected"
1073
        return eval('None')
1074
# End pobyso_get_sage_exp_from_sollya_exp_so_sa
1075

    
1076

    
1077
def pobyso_get_subfunctions(expressionSo):
1078
    """ Legacy function. See pobyso_get_subfunctions_so_sa. """
1079
    return pobyso_get_subfunctions_so_sa(expressionSo) 
1080
# End pobyso_get_subfunctions.
1081
 
1082
def pobyso_get_subfunctions_so_sa(expressionSo):
1083
    """
1084
    Get the subfunctions of an expression.
1085
    Return the number of subfunctions and the list of subfunctions addresses.
1086
    S.T.: Could not figure out another way than that ugly list of declarations
1087
    to recover the addresses of the subfunctions.
1088
    We limit ourselves to arity 8 functions. 
1089
    """
1090
    subf0 = c_int(0)
1091
    subf1 = c_int(0)
1092
    subf2 = c_int(0)
1093
    subf3 = c_int(0)
1094
    subf4 = c_int(0)
1095
    subf5 = c_int(0)
1096
    subf6 = c_int(0)
1097
    subf7 = c_int(0)
1098
    subf8 = c_int(0)
1099
    arity = c_int(0)
1100
    nullPtr = POINTER(c_int)()
1101
    sollya_lib_get_subfunctions(expressionSo, byref(arity), \
1102
      byref(subf0), byref(subf1), byref(subf2), byref(subf3), \
1103
      byref(subf4), byref(subf5),\
1104
      byref(subf6), byref(subf7), byref(subf8), nullPtr, None) 
1105
#    byref(cast(subfunctions[0], POINTER(c_int))), \
1106
#    byref(cast(subfunctions[0], POINTER(c_int))), \
1107
#    byref(cast(subfunctions[2], POINTER(c_int))), \
1108
#    byref(cast(subfunctions[3], POINTER(c_int))), \
1109
#    byref(cast(subfunctions[4], POINTER(c_int))), \
1110
#    byref(cast(subfunctions[5], POINTER(c_int))), \
1111
#    byref(cast(subfunctions[6], POINTER(c_int))), \
1112
#    byref(cast(subfunctions[7], POINTER(c_int))), \
1113
#    byref(cast(subfunctions[8], POINTER(c_int))), nullPtr)
1114
    subfunctions = [subf0, subf1, subf2, subf3, subf4, subf5, subf6, subf7, \
1115
                    subf8]
1116
    subs = []
1117
    if arity.value > pobyso_max_arity:
1118
        return(0,[])
1119
    for i in xrange(arity.value):
1120
        subs.append(int(subfunctions[i].value))
1121
        #print subs[i]
1122
    return (int(arity.value), subs)
1123
# End pobyso_get_subfunctions_so_sa
1124
    
1125
def pobyso_guess_degree_sa_sa(functionSa, intervalSa, approxErrorSa, 
1126
                              weightSa=None, degreeBoundSa=None):
1127
    """
1128
    Sa_sa variant of the solly_guessdegree function.
1129
    Return 0 if something goes wrong.
1130
    """
1131
    functionAsStringSa = functionSa._assume_str().replace('_SAGE_VAR_', '')
1132
    functionSo = pobyso_parse_string_sa_so(functionAsStringSa)
1133
    if pobyso_is_error_so_sa(functionSo):
1134
        sollya_lib_clear_obj(functionSo)
1135
        return 0
1136
    rangeSo = pobyso_interval_to_range_sa_so(intervalSa)
1137
    # The approximation error is expected to be a floating point number.
1138
    if pobyso_is_floating_point_number_sa_sa(approxErrorSa):
1139
        approxErrorSo = pobyso_constant_sa_so(approxErrorSa)
1140
    else:
1141
        approxErrorSo = pobyso_constant_sa_so(RR(approxErrorSa))
1142
    if not weightSa is None:
1143
        weightAsStringSa = weightSa._assume_str().replace('_SAGE_VAR_', '')
1144
        weightSo = pobyso_parse_string_sa_so(weightAsStringSa)
1145
        if pobyso_is_error_so_sa(weightSo):
1146
            sollya_lib_clear_obj(functionSo)
1147
            sollya_lib_clear_obj(rangeSo)
1148
            sollya_lib_clear_obj(approxErrorSo)   
1149
            sollya_lib_clear_obj(weightSo)
1150
            return 0   
1151
    else:
1152
        weightSo = None
1153
    if not degreeBoundSa is None:
1154
        degreeBoundSo = pobyso_constant_from_int_sa_so(degreeBoundSa)
1155
    else:
1156
        degreeBoundSo = None
1157
    guessedDegreeSa = pobyso_guess_degree_so_sa(functionSo,
1158
                                                rangeSo,
1159
                                                approxErrorSo,
1160
                                                weightSo,
1161
                                                degreeBoundSo)
1162
    sollya_lib_clear_obj(functionSo)
1163
    sollya_lib_clear_obj(rangeSo)
1164
    sollya_lib_clear_obj(approxErrorSo)
1165
    if not weightSo is None:
1166
        sollya_lib_clear_obj(weightSo)
1167
    if not degreeBoundSo is None:
1168
        sollya_lib_clear_obj(degreeBoundSo)
1169
    return guessedDegreeSa
1170
# End poyso_guess_degree_sa_sa
1171

    
1172
def pobyso_guess_degree_so_sa(functionSo, rangeSo, errorSo, weightSo=None, \
1173
                              degreeBoundSo=None):
1174
    """
1175
    Thin wrapper around the guessdegree function.
1176
    Nevertheless, some precision control stuff has been appended.
1177
    """
1178
    # Deal with Sollya internal precision issues: if it is too small, 
1179
    # compared with the error, increases it to about twice -log2(error).
1180
    errorSa = pobyso_get_constant_as_rn_with_rf_so_sa(errorSo)
1181
    log2ErrorSa = errorSa.log2()
1182
    if log2ErrorSa < 0:
1183
        neededPrecisionSa = int(2 * int(-log2ErrorSa) / 64) * 64
1184
    else:
1185
        neededPrecisionSa = int(2 * int(log2ErrorSa) / 64) * 64
1186
    #print "Needed precision:", neededPrecisionSa
1187
    sollyaPrecisionHasChanged = False
1188
    (initialPrecSo, initialPrecSa) = pobyso_get_prec_so_so_sa()
1189
    if neededPrecisionSa > initialPrecSa:
1190
        if neededPrecisionSa <= 2:
1191
            print inspect.stack()[0][3], ": precision change <= 2 requested."
1192
        pobyso_set_prec_sa_so(neededPrecisionSa)
1193
        sollyaPrecisionHasChanged = True
1194
    #print "Guessing degree..."
1195
    # weightSo and degreeBoundsSo are optional arguments.
1196
    # As declared, sollya_lib_guessdegree must take 5 arguments.
1197
    if weightSo is None:
1198
        degreeRangeSo = sollya_lib_guessdegree(functionSo, rangeSo, errorSo, 
1199
                                               0, 0, None)
1200
    elif degreeBoundSo is None:
1201
        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, \
1202
                                                errorSo, weightSo, 0, None)
1203
    else:
1204
        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, errorSo, \
1205
                                                weightSo, degreeBoundSo, None)
1206
    #print "...degree guess done."
1207
    # Restore internal precision, if applicable.
1208
    if sollyaPrecisionHasChanged:
1209
        pobyso_set_prec_so_so(initialPrecSo)
1210
        sollya_lib_clear_obj(initialPrecSo)
1211
    degreeIntervalSa = pobyso_range_to_interval_so_sa(degreeRangeSo)
1212
    sollya_lib_clear_obj(degreeRangeSo)
1213
    # When ok, both bounds match.
1214
    # When the degree bound is too low, the upper bound is the degree
1215
    # for which the error can be honored.
1216
    # When it really goes wrong, the upper bound is infinity. 
1217
    if degreeIntervalSa.lower() == degreeIntervalSa.upper():
1218
        return int(degreeIntervalSa.lower())
1219
    else:
1220
        if degreeIntervalSa.upper().is_infinity():
1221
            return None
1222
        else:
1223
            return int(degreeIntervalSa.upper())
1224
    # End pobyso_guess_degree_so_sa
1225

    
1226
def pobyso_inf_so_so(intervalSo):
1227
    """
1228
    Very thin wrapper around sollya_lib_inf().
1229
    """
1230
    return sollya_lib_inf(intervalSo)
1231
# End pobyso_inf_so_so.
1232
    
1233
def pobyso_infnorm_so_so(func, interval, file = None, intervalList = None):
1234
    print "Do not use this function. User pobyso_supnorm_so_so instead."
1235
    return None
1236

    
1237
def pobyso_interval_to_range_sa_so(intervalSa, precisionSa=None):
1238
    if precisionSa is None:
1239
        precisionSa = intervalSa.parent().precision()
1240
    intervalSo = pobyso_bounds_to_range_sa_so(intervalSa.lower(),\
1241
                                              intervalSa.upper(),\
1242
                                              precisionSa)
1243
    return intervalSo
1244
# End pobyso_interval_to_range_sa_so
1245

    
1246
def pobyso_is_error_so_sa(objSo):
1247
    """
1248
    Thin wrapper around the sollya_lib_obj_is_error() function.
1249
    """
1250
    if sollya_lib_obj_is_error(objSo) != 0:
1251
        return True
1252
    else:
1253
        return False
1254
# End pobyso_is_error-so_sa
1255

    
1256
def pobyso_is_floating_point_number_sa_sa(numberSa):
1257
    """
1258
    Check whether a Sage number is floating point.
1259
    Exception stuff added because numbers other than
1260
    floating-point ones do not have the is_real() attribute.
1261
    """
1262
    try:
1263
        return numberSa.is_real()
1264
    except AttributeError:
1265
        return False
1266
# End pobyso_is_floating_piont_number_sa_sa
1267

    
1268
def pobyso_is_range_so_sa(rangeCandidateSo):
1269
    """
1270
    Thin wrapper over sollya_lib_is_range.
1271
    """
1272
    return sollya_lib_obj_is_range(rangeCandidateSo) != 0
1273
    
1274
# End pobyso_is_range_so_sa
1275

    
1276

    
1277
def pobyso_lib_init():
1278
    sollya_lib_init(None)
1279

    
1280
def pobyso_lib_close():
1281
    sollya_lib_close(None)
1282

    
1283
def pobyso_name_free_variable(freeVariableNameSa):
1284
    """ Legacy function. See pobyso_name_free_variable_sa_so. """
1285
    pobyso_name_free_variable_sa_so(freeVariableNameSa)
1286

    
1287
def pobyso_name_free_variable_sa_so(freeVariableNameSa):
1288
    """
1289
    Set the free variable name in Sollya from a Sage string.
1290
    """
1291
    sollya_lib_name_free_variable(freeVariableNameSa)
1292

    
1293
def pobyso_obj_is_function_so_sa(objSo):
1294
    """
1295
    Check if an object is a function.
1296
    """
1297
    if sollya_lib_obj_is_function(objSo) != 0:
1298
        return True
1299
    else:
1300
        return False
1301
# End pobyso_obj_is_function_so_sa  
1302
    
1303
def pobyso_obj_is_range_so_sa(objSo):
1304
    """
1305
    Check if an object is a function.
1306
    """
1307
    if sollya_lib_obj_is_range(objSo) != 0:
1308
        return True
1309
    else:
1310
        return False
1311
# End pobyso_obj_is_range_so_sa  
1312
    
1313
def pobyso_obj_is_string_so_sa(objSo):
1314
    """
1315
    Check if an object is a function.
1316
    """
1317
    if sollya_lib_obj_is_string(objSo) != 0:
1318
        return True
1319
    else:
1320
        return False
1321
# End pobyso_obj_is_string_so_sa  
1322
    
1323
def pobyso_parse_string(string):
1324
    """ Legacy function. See pobyso_parse_string_sa_so. """
1325
    return pobyso_parse_string_sa_so(string)
1326
 
1327
def pobyso_parse_string_sa_so(string):
1328
    """
1329
    Get the Sollya expression computed from a Sage string or
1330
    a Sollya error object if parsing failed.
1331
    """
1332
    return sollya_lib_parse_string(string)
1333

    
1334
def pobyso_precision_so_sa(ctExpSo):
1335
    """
1336
    Computes the necessary precision to represent a number.
1337
    If x is not zero, it can be uniquely written as x = m · 2e 
1338
    where m is an odd integer and e is an integer. 
1339
    precision(x) returns the number of bits necessary to write m 
1340
    in binary (i.e. ceil(log2(m))).
1341
    """
1342
    #TODO: take care of the special case: 0, @NaN@, @Inf@
1343
    precisionSo = sollya_lib_precision(ctExpSo)
1344
    precisionSa = pobyso_constant_from_int_so_sa(precisionSo)
1345
    sollya_lib_clear_obj(precisionSo)
1346
    return precisionSa
1347
# End pobyso_precision_so_sa
1348

    
1349
def pobyso_polynomial_coefficients_progressive_round_so_so(polySo,
1350
                                                           funcSo,
1351
                                                           icSo,
1352
                                                           intervalSo,
1353
                                                           itpSo,
1354
                                                           ftpSo,
1355
                                                           maxPrecSo,
1356
                                                           maxErrSo,
1357
                                                           debug=False):
1358
    if debug:
1359
        print "Input arguments:"
1360
        pobyso_autoprint(polySo)
1361
        pobyso_autoprint(funcSo)
1362
        pobyso_autoprint(icSo)
1363
        pobyso_autoprint(intervalSo)
1364
        pobyso_autoprint(itpSo)
1365
        pobyso_autoprint(ftpSo)
1366
        pobyso_autoprint(maxPrecSo)
1367
        pobyso_autoprint(maxErrSo)
1368
        print "________________"
1369
    
1370
    ## Higher order function see: 
1371
    #  http://effbot.org/pyfaq/how-do-you-make-a-higher-order-function-in-python.htm
1372
    def precision_decay_ratio_function(degreeSa):
1373
        def outer(x):
1374
            def inner(x):
1375
                we = 3/8
1376
                wq = 2/8
1377
                a  = 2.2
1378
                b  = 2 
1379
                return we*(exp(x/a)-1) +  wq*((b*x)**2) + (1-we-wq)*x
1380
            return  inner(x)/inner(degreeSa)
1381
        return outer
1382
    
1383
    #
1384
    degreeSa        = pobyso_polynomial_degree_so_sa(polySo)
1385
    ratio           = precision_decay_ratio_function(degreeSa)
1386
    itpSa           = pobyso_constant_from_int_so_sa(itpSo)
1387
    ftpSa           = pobyso_constant_from_int_so_sa(ftpSo)
1388
    maxPrecSa       = pobyso_constant_from_int_so_sa(maxPrecSo)
1389
    maxErrSa        = pobyso_get_constant_as_rn_so_sa(maxErrSo)
1390
    if debug:
1391
        print "degreeSa:", degreeSa
1392
        print "ratio:", ratio
1393
        print "itpsSa:", itpSa
1394
        print "ftpSa:", ftpSa
1395
        print "maxPrecSa:", maxPrecSa
1396
        print "maxErrSa:", maxErrSa
1397
    lastResPolySo   = None
1398
    lastInfNormSo   = None
1399
    #print "About to enter the while loop..."
1400
    while True: 
1401
        resPolySo   = pobyso_constant_0_sa_so()
1402
        pDeltaSa    = ftpSa - itpSa
1403
        for indexSa in reversed(xrange(0,degreeSa+1)):
1404
            #print "Index:", indexSa
1405
            indexSo = pobyso_constant_from_int_sa_so(indexSa)
1406
            #pobyso_autoprint(indexSo)
1407
            #print ratio(indexSa)
1408
            ctpSa = floor(ftpSa - (pDeltaSa * ratio(indexSa)))
1409
            ctpSo = pobyso_constant_from_int_sa_so(ctpSa)
1410
            if debug:
1411
                print "Index:", indexSa, " - Target precision:", 
1412
                pobyso_autoprint(ctpSo)
1413
            cmonSo  = \
1414
                sollya_lib_build_function_mul(sollya_lib_coeff(polySo, indexSo),
1415
                                      sollya_lib_build_function_pow( \
1416
                                          sollya_lib_build_function_free_variable(), \
1417
                                          indexSo))
1418
            #pobyso_autoprint(cmonSo)
1419
            cmonrSo = pobyso_round_coefficients_single_so_so(cmonSo, ctpSo)
1420
            sollya_lib_clear_obj(cmonSo)
1421
            #pobyso_autoprint(cmonrSo)
1422
            resPolySo = sollya_lib_build_function_add(resPolySo,
1423
                                                      cmonrSo)
1424
            #pobyso_autoprint(resPolySo)
1425
        # End for index
1426
        freeVarSo     = sollya_lib_build_function_free_variable()
1427
        changeVarSo   = sollya_lib_sub(freeVarSo, icSo)
1428
        resPolyCvSo   = sollya_lib_evaluate(resPolySo, changeVarSo)
1429
        errFuncSo = sollya_lib_build_function_sub(sollya_lib_copy_obj(funcSo), 
1430
                                                  resPolyCvSo)
1431
        infNormSo = sollya_lib_dirtyinfnorm(errFuncSo, intervalSo)
1432
        cerrSa    = pobyso_get_constant_as_rn_so_sa(infNormSo)
1433
        if debug:
1434
            print "Infnorm (Sollya):", 
1435
            pobyso_autoprint(infNormSo)
1436
        sollya_lib_clear_obj(errFuncSo)
1437
        #print "Infnorm  (Sage):", cerrSa
1438
        if (cerrSa > maxErrSa):
1439
            if debug:
1440
                print "Error is too large."
1441
            if lastResPolySo is None:
1442
                if debug:
1443
                    print "Enlarging prec."
1444
                ntpSa = floor(ftpSa + ftpSa/50)
1445
                ## Can't enlarge (numerical)
1446
                if ntpSa == ftpSa:
1447
                    sollya_lib_clear_obj(resPolySo)
1448
                    return None
1449
                ## Can't enlarge (not enough precision left)
1450
                if ntpSa > maxPrecSa:
1451
                    sollya_lib_clear_obj(resPolySo)
1452
                    return None
1453
                ftpSa = ntpSa
1454
                continue
1455
            ## One enlargement took place.
1456
            else:
1457
                if debug:
1458
                    print "Exit with the last before last polynomial."
1459
                    print "Precision of highest degree monomial:", itpSa
1460
                    print "Precision of constant term          :", ftpSa
1461
                sollya_lib_clear_obj(resPolySo)
1462
                sollya_lib_clear_obj(infNormSo)
1463
                return (lastResPolySo, lastInfNormSo)
1464
        # cerrSa <= maxErrSa: scrap more bits, possibly.
1465
        else:
1466
            if debug:
1467
                print "Error is too small"
1468
            if cerrSa <= (maxErrSa/2):
1469
                if debug: 
1470
                    print "Shrinking prec."
1471
                ntpSa = floor(ftpSa - ftpSa/50)
1472
                ## Can't shrink (numerical)
1473
                if ntpSa == ftpSa:
1474
                    if not lastResPolySo is None:
1475
                        sollya_lib_clear_obj(lastResPolySo)
1476
                    if not lastInfNormSo is None:
1477
                        sollya_lib_clear_obj(lastInfNormSo)
1478
                    if debug:
1479
                        print "Exit because can't shrink anymore (numerically)."
1480
                        print "Precision of highest degree monomial:", itpSa
1481
                        print "Precision of constant term          :", ftpSa
1482
                    return (resPolySo, infNormSo)
1483
                ## Can't shrink (not enough precision left)
1484
                if ntpSa <= itpSa:
1485
                    if not lastResPolySo is None:
1486
                        sollya_lib_clear_obj(lastResPolySo)
1487
                    if not lastInfNormSo is None:
1488
                        sollya_lib_clear_obj(lastInfNormSo)
1489
                        print "Exit because can't shrink anymore (no bits left)."
1490
                        print "Precision of highest degree monomial:", itpSa
1491
                        print "Precision of constant term          :", ftpSa
1492
                    return (resPolySo, infNormSo)
1493
                ftpSa = ntpSa
1494
                if not lastResPolySo is None:
1495
                    sollya_lib_clear_obj(lastResPolySo)
1496
                if not lastInfNormSo is None:
1497
                    sollya_lib_clear_obj(lastInfNormSo)
1498
                lastResPolySo = resPolySo
1499
                lastInfNormSo = infNormSo
1500
                continue
1501
            else: # Error is not that small, just return
1502
                if not lastResPolySo is None:
1503
                    sollya_lib_clear_obj(lastResPolySo)
1504
                if not lastInfNormSo is None:
1505
                    sollya_lib_clear_obj(lastInfNormSo)
1506
                if debug:
1507
                    print "Exit normally."
1508
                    print "Precision of highest degree monomial:", itpSa
1509
                    print "Precision of constant term          :", ftpSa
1510
                return (resPolySo, infNormSo)                
1511
    # End wile True
1512
    return None
1513
# End pobyso_polynomial_coefficients_progressive_truncate_so_so.
1514

    
1515
def pobyso_polynomial_degree_so_sa(polySo):
1516
    """
1517
    Return the degree of a Sollya polynomial as a Sage int.
1518
    """
1519
    degreeSo = sollya_lib_degree(polySo)
1520
    return pobyso_constant_from_int_so_sa(degreeSo)
1521
# End pobyso_polynomial_degree_so_sa
1522

    
1523
def pobyso_polynomial_degree_so_so(polySo):
1524
    """
1525
    Thin wrapper around lib_sollya_degree().
1526
    """
1527
    return sollya_lib_degree(polySo)
1528
# End pobyso_polynomial_degree_so_so
1529
                                                                  
1530
def pobyso_range(rnLowerBound, rnUpperBound):
1531
    """ Legacy function. See pobyso_range_sa_so. """
1532
    return pobyso_range_sa_so(rnLowerBound, rnUpperBound) 
1533

    
1534
def pobyso_range_from_bounds_sa_so(rnLowerBound, rnUpperBound, precSa = None):
1535
    """
1536
    Create a Sollya range from 2 Sage real numbers as bounds
1537
    """
1538
    # TODO check precision stuff.
1539
    sollyaPrecChanged = False
1540
    (initialSollyaPrecSo, initialSollyaPrecSa) = \
1541
        pobyso_get_prec_so_so_sa()
1542
    if precSa is None:
1543
        precSa = max(rnLowerBound.parent().prec(), rnUpperBound.parent().prec())
1544
    if precSa > initialSollyaPrecSa:
1545
        if precSa <= 2:
1546
            print inspect.stack()[0][3], ": precision change <= 2 requested."
1547
        pobyso_set_prec_sa_so(precSa)
1548
        sollyaPrecChanged = True
1549
    rangeSo = sollya_lib_range_from_bounds(get_rn_value(rnLowerBound),
1550
                                           get_rn_value(rnUpperBound))
1551
    if sollyaPrecChanged:
1552
        pobyso_set_prec_so_so(initialSollyaPrecSo)
1553
        sollya_lib_clear_obj(initialSollyaPrecSo)
1554
    return rangeSo
1555
# End pobyso_range_from_bounds_sa_so
1556

    
1557
def pobyso_range_list_so_sa(listSo):
1558
    """
1559
    Return a Sollya list of ranges as a Sage list of 
1560
    floating-point intervals.
1561
    """
1562
    listSa   = []
1563
    ## The function returns none if the list is empty or an error has happened.
1564
    retVal = pobyso_get_list_elements_so_so(listSo)
1565
    if retVal is None:
1566
        return listSa
1567
    ## Just in case the interface is changed and an empty list is returned 
1568
    #  instead of None.
1569
    elif len(retVal) == 0:
1570
        return listSa
1571
    else:
1572
        ## Remember pobyso_get_list_elements_so_so returns more information
1573
        #  than just the elements of the list (# elements, is_elliptic)
1574
        listSaSo, numElements, isEndElliptic = retVal
1575
    ## Return an empty list.
1576
    if numElements == 0:
1577
        return listSa
1578
    ## Search first for the maximum precision of the elements
1579
    maxPrecSa = 0
1580
    for rangeSo in listSaSo:
1581
        #pobyso_autoprint(floatSo)
1582
        curPrecSa =  pobyso_get_prec_of_range_so_sa(rangeSo)
1583
        if curPrecSa > maxPrecSa:
1584
            maxPrecSa = curPrecSa
1585
    ##
1586
    intervalField = RealIntervalField(maxPrecSa)
1587
    ##
1588
    for rangeSo in listSaSo:
1589
        listSa.append(pobyso_range_to_interval_so_sa(rangeSo, intervalField))
1590
    return listSa
1591
# End pobyso_range_list_so_sa
1592

    
1593
def pobyso_range_max_abs_so_so(rangeSo):
1594
    """
1595
    Return, as Sollya constant, the maximum absolute value of a Sollay range.
1596
    """
1597
    lowerBoundSo = sollya_lib_inf(rangeSo)
1598
    upperBoundSo = sollya_lib_sup(rangeSo)
1599
    #
1600
    lowerBoundSo = sollya_lib_build_function_abs(lowerBoundSo)
1601
    upperBoundSo = sollya_lib_build_function_abs(upperBoundSo)
1602
    #pobyso_autoprint(lowerBoundSo)
1603
    #pobyso_autoprint(upperBoundSo)
1604
    #
1605
    maxAbsSo = sollya_lib_max(lowerBoundSo, upperBoundSo, None)
1606
    #sollya_lib_clear_obj(lowerBoundSo)
1607
    #sollya_lib_clear_obj(upperBoundSo)
1608
    return maxAbsSo
1609
# End pobyso_range_max_abs_so_so
1610
 
1611
def pobyso_range_to_interval_so_sa(rangeSo, realIntervalFieldSa = None):
1612
    """
1613
    Get a Sage interval from a Sollya range.
1614
    If no realIntervalField is given as a parameter, the Sage interval
1615
    precision is that of the Sollya range.
1616
    Otherwise, the precision is that of the realIntervalField. In this case
1617
    rounding may happen.
1618
    """
1619
    if realIntervalFieldSa is None:
1620
        precSa = pobyso_get_prec_of_range_so_sa(rangeSo)
1621
        realIntervalFieldSa = RealIntervalField(precSa)
1622
    intervalSa = \
1623
        pobyso_get_interval_from_range_so_sa(rangeSo, realIntervalFieldSa) 
1624
    return intervalSa
1625
# End pobyso_range_to_interval_so_sa
1626
#
1627
def pobyso_relative_so_so():
1628
    """
1629
    Very thin wrapper around the sollya_lib_relative function.
1630
    """
1631
    return sollya_lib_relative()
1632
# End pobyso_relative_so_so
1633
#
1634
def pobyso_rat_poly_sa_so(polySa, precSa = None):
1635
    """
1636
    Create a Sollya polynomial from a Sage rational polynomial.
1637
    We first convert the rational polynomial into a floating-point 
1638
    polynomial.
1639
    """
1640
    ## TODO: filter arguments.
1641
    ## Precision. If no precision is given, use the current precision
1642
    #  of Sollya.
1643
    if precSa is None:
1644
        precSa =  pobyso_get_prec_so_sa()
1645
    #print "Precision:",  precSa
1646
    RRR = RealField(precSa)
1647
    ## Create a Sage polynomial in the "right" precision.
1648
    P_RRR = RRR[polySa.variables()[0]]
1649
    polyFloatSa = P_RRR(polySa)
1650
    ## Make sure no precision is provided: pobyso_float_poly_sa_so will
1651
    #  recover it all by itself and will not make any extra conversion.
1652
    return pobyso_float_poly_sa_so(polyFloatSa)
1653
    
1654
# End pobyso_rat_poly_sa_so    
1655
    
1656
def pobyso_remez_canonical_sa_sa(func, \
1657
                                 degree, \
1658
                                 lowerBound, \
1659
                                 upperBound, \
1660
                                 weight = None, \
1661
                                 quality = None):
1662
    """
1663
    All arguments are Sage/Python.
1664
    The functions (func and weight) must be passed as expressions or strings.
1665
    Otherwise the function fails. 
1666
    The return value is a Sage polynomial.
1667
    """
1668
    var('zorglub')    # Dummy variable name for type check only. Type of 
1669
    # zorglub is "symbolic expression".
1670
    polySo = pobyso_remez_canonical_sa_so(func, \
1671
                                 degree, \
1672
                                 lowerBound, \
1673
                                 upperBound, \
1674
                                 weight, \
1675
                                 quality)
1676
    # String test
1677
    if parent(func) == parent("string"):
1678
        functionSa = eval(func)
1679
    # Expression test.
1680
    elif type(func) == type(zorglub):
1681
        functionSa = func
1682
    else:
1683
        return None
1684
    #
1685
    maxPrecision = 0
1686
    if polySo is None:
1687
        return(None)
1688
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
1689
    RRRRSa = RealField(maxPrecision)
1690
    polynomialRingSa = RRRRSa[functionSa.variables()[0]]
1691
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, RRRRSa)
1692
    polySa = polynomial(expSa, polynomialRingSa)
1693
    sollya_lib_clear_obj(polySo)
1694
    return(polySa)
1695
# End pobyso_remez_canonical_sa_sa
1696
    
1697
def pobyso_remez_canonical(func, \
1698
                           degree, \
1699
                           lowerBound, \
1700
                           upperBound, \
1701
                           weight = "1", \
1702
                           quality = None):
1703
    """ Legacy function. See pobyso_remez_canonical_sa_so. """
1704
    return(pobyso_remez_canonical_sa_so(func, \
1705
                                        degree, \
1706
                                        lowerBound, \
1707
                                        upperBound, \
1708
                                        weight, \
1709
                                        quality))
1710
# End pobyso_remez_canonical.
1711

    
1712
def pobyso_remez_canonical_sa_so(func, \
1713
                                 degree, \
1714
                                 lowerBound, \
1715
                                 upperBound, \
1716
                                 weight = None, \
1717
                                 quality = None):
1718
    """
1719
    All arguments are Sage/Python.
1720
    The functions (func and weight) must be passed as expressions or strings.
1721
    Otherwise the function fails. 
1722
    The return value is a pointer to a Sollya function.
1723
    lowerBound and upperBound mus be reals.
1724
    """
1725
    var('zorglub')    # Dummy variable name for type check only. Type of
1726
    # zorglub is "symbolic expression".
1727
    currentVariableNameSa = None
1728
    # The func argument can be of different types (string, 
1729
    # symbolic expression...)
1730
    if parent(func) == parent("string"):
1731
        localFuncSa = sage_eval(func,globals())
1732
        if len(localFuncSa.variables()) > 0:
1733
            currentVariableNameSa = localFuncSa.variables()[0]
1734
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1735
            functionSo = \
1736
              sollya_lib_parse_string(localFuncSa._assume_str().replace('_SAGE_VAR_', ''))
1737
    # Expression test.
1738
    elif type(func) == type(zorglub):
1739
        # Until we are able to translate Sage expressions into Sollya 
1740
        # expressions : parse the string version.
1741
        if len(func.variables()) > 0:
1742
            currentVariableNameSa = func.variables()[0]
1743
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1744
            functionSo = \
1745
              sollya_lib_parse_string(func._assume_str().replace('_SAGE_VAR_', ''))
1746
    else:
1747
        return(None)
1748
    if weight is None: # No weight given -> 1.
1749
        weightSo = pobyso_constant_1_sa_so()
1750
    elif parent(weight) == parent("string"): # Weight given as string: parse it.
1751
        weightSo = sollya_lib_parse_string(func)
1752
    elif type(weight) == type(zorglub): # Weight given as symbolice expression.
1753
        functionSo = \
1754
          sollya_lib_parse_string_sa_so(weight._assume_str().replace('_SAGE_VAR_', ''))
1755
    else:
1756
        return(None)
1757
    degreeSo = pobyso_constant_from_int(degree)
1758
    rangeSo = pobyso_bounds_to_range_sa_so(lowerBound, upperBound)
1759
    if not quality is None:
1760
        qualitySo= pobyso_constant_sa_so(quality)
1761
    else:
1762
        qualitySo = None
1763
        
1764
    remezPolySo = sollya_lib_remez(functionSo, \
1765
                                   degreeSo, \
1766
                                   rangeSo, \
1767
                                   weightSo, \
1768
                                   qualitySo, \
1769
                                   None)
1770
    sollya_lib_clear_obj(functionSo)
1771
    sollya_lib_clear_obj(degreeSo)
1772
    sollya_lib_clear_obj(rangeSo)
1773
    sollya_lib_clear_obj(weightSo)
1774
    if not qualitySo is None:
1775
        sollya_lib_clear_obj(qualitySo)
1776
    return(remezPolySo)
1777
# End pobyso_remez_canonical_sa_so
1778

    
1779
def pobyso_remez_canonical_so_so(funcSo, \
1780
                                 degreeSo, \
1781
                                 rangeSo, \
1782
                                 weightSo = pobyso_constant_1_sa_so(),\
1783
                                 qualitySo = None):
1784
    """
1785
    All arguments are pointers to Sollya objects.
1786
    The return value is a pointer to a Sollya function.
1787
    """
1788
    if not sollya_lib_obj_is_function(funcSo):
1789
        return(None)
1790
    return(sollya_lib_remez(funcSo, degreeSo, rangeSo, weightSo, qualitySo, None))
1791
# End pobyso_remez_canonical_so_so.
1792

    
1793
def pobyso_remez_exponents_list_sa_so(func,     \
1794
                                 exponentsList, \
1795
                                 lowerBound,    \
1796
                                 upperBound,    \
1797
                                 weight = None, \
1798
                                 quality = None):
1799
    """
1800
    All arguments are Sage/Python.
1801
    The functions (func and weight) must be passed as expressions or strings.
1802
    Otherwise the function fails. 
1803
    The return value is a pointer to a Sollya function.
1804
    lowerBound and upperBound mus be reals.
1805
    """
1806
    var('zorglub')    # Dummy variable name for type check only. Type of
1807
    # zorglub is "symbolic expression".
1808
    currentVariableNameSa = None
1809
    # The func argument can be of different types (string, 
1810
    # symbolic expression...)
1811
    if parent(func) == parent("string"):
1812
        localFuncSa = sage_eval(func,globals())
1813
        if len(localFuncSa.variables()) > 0:
1814
            currentVariableNameSa = localFuncSa.variables()[0]
1815
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1816
            functionSo = \
1817
              sollya_lib_parse_string(localFuncSa._assume_str().replace('_SAGE_VAR_', ''))
1818
    # Expression test.
1819
    elif type(func) == type(zorglub):
1820
        # Until we are able to translate Sage expressions into Sollya 
1821
        # expressions : parse the string version.
1822
        if len(func.variables()) > 0:
1823
            currentVariableNameSa = func.variables()[0]
1824
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1825
            functionSo = \
1826
              sollya_lib_parse_string(func._assume_str().replace('_SAGE_VAR_', ''))
1827
    else:
1828
        return(None)
1829
    ## Deal with the weight, much in the same way as with the function.
1830
    if weight is None: # No weight given -> 1.
1831
        weightSo = pobyso_constant_1_sa_so()
1832
    elif parent(weight) == parent("string"): # Weight given as string: parse it.
1833
        weightSo = sollya_lib_parse_string(func)
1834
    elif type(weight) == type(zorglub): # Weight given as symbolice expression.
1835
        functionSo = \
1836
          sollya_lib_parse_string(weight._assume_str().replace('_SAGE_VAR_', '',100))
1837
    else:
1838
        return(None)
1839
    rangeSo = pobyso_bounds_to_range_sa_so(lowerBound, upperBound)
1840
    if not quality is None:
1841
        qualitySo= pobyso_constant_sa_so(quality)
1842
    else:
1843
        qualitySo = None
1844
    #
1845
    ## Tranform the Sage list of exponents into a Sollya list.
1846
    exponentsListSo = pobyso_build_list_of_ints_sa_so(*exponentsList)
1847
    remezPolySo = sollya_lib_remez(functionSo,      \
1848
                                   exponentsListSo, \
1849
                                   rangeSo,         \
1850
                                   weightSo,        \
1851
                                   qualitySo,       \
1852
                                   None)
1853
    sollya_lib_clear_obj(functionSo)
1854
    sollya_lib_clear_obj(exponentsListSo)
1855
    sollya_lib_clear_obj(rangeSo)
1856
    sollya_lib_clear_obj(weightSo)
1857
    if not qualitySo is None:
1858
        sollya_lib_clear_obj(qualitySo)
1859
    return(remezPolySo)
1860
# End pobyso_remez_exponentsList_sa_so
1861
#
1862
def pobyso_round_coefficients_so_so(polySo, truncFormatListSo):
1863
    """
1864
    A wrapper around the "classical" sollya_lib_roundcoefficients: a Sollya
1865
    polynomial and a Sollya list are given as arguments.
1866
    """
1867
    return sollya_lib_roundcoefficients(polySo, truncFormatListSo)
1868

    
1869
def pobyso_round_coefficients_progressive_so_so(polySo, 
1870
                                                funcSo,
1871
                                                precSo,
1872
                                                intervalSo,
1873
                                                icSo,
1874
                                                currentApproxErrorSo,
1875
                                                approxAccurSo,
1876
                                                debug=False):
1877
    """
1878
    From an input approximation polynomial, compute an output one with 
1879
    smaller coefficients and yet yields a sufficient approximation accuracy.
1880
    """
1881
    if debug:
1882
        print "Input arguments:"
1883
        print "Polynomial: ", ; pobyso_autoprint(polySo)
1884
        print "Function: ", ; pobyso_autoprint(funcSo)
1885
        print "Internal precision: ", ; pobyso_autoprint(precSo)
1886
        print "Interval: ", ; pobyso_autoprint(intervalSo)
1887
        print "Current approximation error: ", ; pobyso_autoprint(currentApproxErrorSo)
1888
        print "Requested approxiation error: ", ; pobyso_autoprint(approxAccurSo)
1889
        print "________________"
1890
    approxAccurSa        = pobyso_get_constant_as_rn_so_sa(approxAccurSo)
1891
    currentApproxErrorSa = pobyso_get_constant_as_rn_so_sa(currentApproxErrorSo)
1892
    ## If the current approximation error is too close to the target, there is
1893
    #  no possible gain.
1894
    if currentApproxErrorSa >= approxAccurSa / 2:
1895
        #### Do not return the initial argument but copies: caller may free 
1896
        #    the former as inutile after call.
1897
        return (sollya_lib_copy_obj(polySo),
1898
                sollya_lib_copy_obj(currentApproxErrorSo))
1899
    #
1900
    ## Try to round the coefficients.
1901
    degreeSa             = pobyso_polynomial_degree_so_sa(polySo)
1902
    intervalSa           = pobyso_range_to_interval_so_sa(intervalSo)
1903
    
1904
    if debug:
1905
        print "degreeSa              :", degreeSa
1906
        print "intervalSa            :", intervalSa.str(style='brackets')
1907
        print "currentApproxErrorSa  :", currentApproxErrorSa 
1908
        print "approxAccurSa         :", approxAccurSa 
1909
    radiusSa = intervalSa.absolute_diameter() / 2
1910
    if debug:
1911
        print "log2(radius):", RR(radiusSa).log2()
1912
    iterIndex = 0
1913
    ## Build the "shaved" polynomial.
1914
    while True: 
1915
        ### Start with a 0 value expression.
1916
        resPolySo = pobyso_constant_0_sa_so()
1917
        roundedPolyApproxAccurSa = approxAccurSa / 2
1918
        currentRadiusPowerSa = 1 
1919
        for degree in xrange(0,degreeSa + 1):
1920
            #### At round 0, use the agressive formula. At round 1, run the
1921
            #    proved formula.
1922
            if iterIndex == 0:
1923
                roundingPowerSa = \
1924
                    floor(((currentRadiusPowerSa/roundedPolyApproxAccurSa)*(degree+1)).log2())
1925
            else:
1926
                roundingPowerSa = \
1927
                    floor(((currentRadiusPowerSa/roundedPolyApproxAccurSa)*(degreeSa+1)).log2())
1928
            ## Under extreme conditions the above formulas can evaluate under 2,
1929
            #   which is the minimal precision of an MPFR number.
1930
            if roundingPowerSa < 2:
1931
                roundingPowerSa = 2
1932
            if debug:
1933
                print "roundedPolyApproxAccurSa", roundedPolyApproxAccurSa
1934
                print "currentRadiusPowerSa", currentRadiusPowerSa
1935
                print "Current rounding exponent:", roundingPowerSa
1936
            currentRadiusPowerSa *= radiusSa
1937
            index1So = pobyso_constant_from_int_sa_so(degree)
1938
            index2So = pobyso_constant_from_int_sa_so(degree)
1939
            ### Create a monomial with:
1940
            #   - the coefficient in the initial monomial at the current degrree;
1941
            #   - the current exponent;
1942
            #   - the free variable.
1943
            cmonSo  = \
1944
                sollya_lib_build_function_mul(sollya_lib_coeff(polySo, index1So),
1945
                                              sollya_lib_build_function_pow( \
1946
                                                sollya_lib_build_function_free_variable(), \
1947
                                                index2So))
1948
            roundingPowerSo = pobyso_constant_from_int_sa_so(roundingPowerSa)
1949
            cmonrSo = pobyso_round_coefficients_single_so_so(cmonSo, roundingPowerSo)
1950
            sollya_lib_clear_obj(cmonSo)
1951
            ### Add to the result polynomial.
1952
            resPolySo = sollya_lib_build_function_add(resPolySo,
1953
                                                      cmonrSo)
1954
        # End for.
1955
        ### Check the new polynomial.
1956
        freeVarSo     = sollya_lib_build_function_free_variable()
1957
        changeVarSo   = sollya_lib_sub(freeVarSo, icSo)
1958
        resPolyCvSo   = sollya_lib_evaluate(resPolySo, changeVarSo)
1959
        errFuncSo     = sollya_lib_build_function_sub(sollya_lib_copy_obj(funcSo), 
1960
                                                      resPolyCvSo)
1961
        infNormSo = sollya_lib_dirtyinfnorm(errFuncSo, intervalSo)
1962
        ### This also clears resPolyCvSo.
1963
        sollya_lib_clear_obj(errFuncSo)
1964
        cerrSa    = pobyso_get_constant_as_rn_so_sa(infNormSo)
1965
        if debug:
1966
            print "Error of the new polynomial:", cerrSa
1967
        ### If at round 1, return the initial polynomial error. This should
1968
        #   never happen since the rounding algorithm is proved. But some 
1969
        #   circumstances may break it (e.g. internal precision of tools).
1970
        if cerrSa > approxAccurSa:
1971
            if iterIndex == 0: # Round 0 is agressive rounding, got round 1 (proved rounding)
1972
                sollya_lib_clear_obj(resPolySo)
1973
                sollya_lib_clear_obj(infNormSo)
1974
                iterIndex += 1
1975
                continue
1976
            else: # Round 1 and beyond : just return the oroginal polynomial.
1977
                sollya_lib_clear_obj(resPolySo)
1978
                sollya_lib_clear_obj(infNormSo)
1979
                #### Do not return the arguments but copies: the caller may free
1980
                #    free the former as inutile after call.
1981
                return (sollya_lib_copy_obj(polySo), 
1982
                        sollya_lib_copy_obj(currentApproxErrorSo))
1983
        ### If get here it is because cerrSa <= approxAccurSa
1984
        ### Approximation error of the new polynomial is acceptable.
1985
        return (resPolySo, infNormSo)
1986
    # End while True
1987
# End pobyso_round_coefficients_progressive_so_so
1988
    
1989
def pobyso_round_coefficients_single_so_so(polySo, commonPrecSo):
1990
    """
1991
    Create a rounded coefficients polynomial from polynomial argument to
1992
    the number of bits in size argument.
1993
    All coefficients are set to the same precision.
1994
    """
1995
    ## TODO: check arguments.
1996
    endEllipListSo = pobyso_build_end_elliptic_list_so_so(commonPrecSo)
1997
    polySo = sollya_lib_roundcoefficients(polySo, endEllipListSo, None)
1998
    sollya_lib_clear_obj(endEllipListSo)
1999
    #sollya_lib_clear_obj(endEllipListSo)
2000
    return polySo
2001
    
2002
# End pobyso_round_coefficients_single_so_so
2003

    
2004
def pobyso_set_canonical_off():
2005
    sollya_lib_set_canonical(sollya_lib_off())
2006

    
2007
def pobyso_set_canonical_on():
2008
    sollya_lib_set_canonical(sollya_lib_on())
2009

    
2010
def pobyso_set_prec(p):
2011
    """ Legacy function. See pobyso_set_prec_sa_so. """
2012
    pobyso_set_prec_sa_so(p)
2013

    
2014
def pobyso_set_prec_sa_so(p):
2015
    #a = c_int(p)
2016
    #precSo = c_void_p(sollya_lib_constant_from_int(a))
2017
    #precSo = sollya_lib_constant_from_int(a)
2018
    precSo = pobyso_constant_from_int_sa_so(p)
2019
    sollya_lib_set_prec(precSo)
2020
    sollya_lib_clear_obj(precSo)
2021
# End pobyso_set_prec_sa_so.
2022

    
2023
def pobyso_set_prec_so_so(newPrecSo):
2024
    sollya_lib_set_prec(newPrecSo)
2025
# End pobyso_set_prec_so_so.
2026

    
2027
def pobyso_inf_so_so(rangeSo):
2028
    """
2029
    Very thin wrapper around sollya_lib_inf().
2030
    """
2031
    return sollya_lib_inf(rangeSo)
2032
# End pobyso_inf_so_so.
2033
#   
2034
def pobyso_infnorm_sa_sa(funcSa, interevalSa):
2035
    """
2036
    Very thin wrapper around sollya_lib_infnorm().
2037
    We only take into account the 2 first arguments (the function and
2038
    the interval (a range). Managing the other arguments (the file for
2039
    the proof and the exclusion intervals list) will be performed later
2040
    Changes will be needed in sollya_lib.py file too.
2041
    """
2042
    return sollya_lib_infnorm(funcSo, rangeSo, None)
2043
# End pobyso_infnorm_so_so.
2044
#   
2045
def pobyso_infnorm_so_so(funcSo, rangeSo):
2046
    """
2047
    Very thin wrapper around sollya_lib_infnorm().
2048
    We only take into account the 2 first arguments (the function and
2049
    the interval (a range). Managing the other arguments (the file for
2050
    the proof and the exclusion intervals list) will be performed later
2051
    Changes will be needed in sollya_lib.py file too.
2052
    """
2053
    return sollya_lib_infnorm(funcSo, rangeSo, None)
2054
# End pobyso_infnorm_so_so.
2055
#   
2056
def pobyso_supnorm_sa_sa(poly):
2057
    """
2058
    Computes the supremum norm from Sage input arguments and returns a
2059
    Sage floating-point number whose precision is set by the realFieldSa
2060
    argument.
2061
    TODO: complete this stub!
2062
    """
2063
    print("This function does nothing!")
2064
    return None
2065
# End pobyso_supnorm_sa_sa
2066

    
2067
def pobyso_supnorm_so_sa(polySo, funcSo, intervalSo, errorTypeSo = None,\
2068
                         accuracySo = None, realFieldSa = None):
2069
    """
2070
    Computes the supremum norm from Sollya input arguments and returns a
2071
    Sage floating-point number whose precision is set by the only Sage argument.
2072
    
2073
    The returned value is the maximum of the absolute values of the range
2074
    elements  returned by the Sollya supnorm functions
2075
    """
2076
    supNormRangeSo = pobyso_supnorm_so_so(polySo, 
2077
                                          funcSo, 
2078
                                          intervalSo, 
2079
                                          errorTypeSo,
2080
                                          accuracySo)
2081
    supNormSo = pobyso_range_max_abs_so_so(supNormRangeSo)
2082
    sollya_lib_clear_obj(supNormRangeSo)
2083
    #pobyso_autoprint(supNormSo)
2084
    supNormSa = pobyso_get_constant_as_rn_with_rf_so_sa(supNormSo, realFieldSa)
2085
    sollya_lib_clear_obj(supNormSo)
2086
    return supNormSa 
2087
# End pobyso_supnorm_so_sa.
2088
#
2089
def pobyso_supnorm_so_so(polySo, funcSo, intervalSo, errorTypeSo = None,\
2090
                         accuracySo = None):
2091
    """
2092
    Computes the supnorm of the approximation error between the given 
2093
    polynomial and function. Attention: returns a range!
2094
    errorTypeSo defaults to "absolute".
2095
    accuracySo defaults to 2^(-40).
2096
    """
2097
    if errorTypeSo is None:
2098
        errorTypeSo = sollya_lib_absolute(None)
2099
        errorTypeIsNone = True
2100
    else:
2101
        errorTypeIsNone = False
2102
    #
2103
    if accuracySo is None:
2104
        # Notice the **: we are in Pythonland here!
2105
        accuracySo = pobyso_constant_sa_so(RR(2**(-40)))
2106
        accuracyIsNone = True
2107
    else:
2108
        accuracyIsNone = False
2109
    #pobyso_autoprint(accuracySo)
2110
    resultSo = \
2111
        sollya_lib_supnorm(polySo, funcSo, intervalSo, errorTypeSo, \
2112
                              accuracySo)
2113
    if errorTypeIsNone:
2114
        sollya_lib_clear_obj(errorTypeSo)
2115
    if accuracyIsNone:
2116
        sollya_lib_clear_obj(accuracySo)
2117
    return resultSo
2118
# End pobyso_supnorm_so_so
2119
#
2120
def pobyso_taylor_expansion_no_change_var_so_so(functionSo, 
2121
                                                degreeSo, 
2122
                                                rangeSo,
2123
                                                errorTypeSo=None, 
2124
                                                sollyaPrecSo=None):
2125
    """
2126
    Compute the Taylor expansion without the variable change
2127
    x -> x-intervalCenter.
2128
    If errorTypeSo is None, absolute is used.
2129
    If sollyaPrecSo is None, Sollya internal precision is not changed. 
2130
    """
2131
    # Change internal Sollya precision, if needed.
2132
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
2133
    sollyaPrecChanged = False
2134
    if sollyaPrecSo is None:
2135
        pass
2136
    else:
2137
        sollya_lib_set_prec(sollyaPrecSo)
2138
        sollyaPrecChanged = True
2139
    # Error type stuff: default to absolute.
2140
    if errorTypeSo is None:
2141
        errorTypeIsNone = True
2142
        errorTypeSo = sollya_lib_absolute(None)
2143
    else:
2144
        errorTypeIsNone = False
2145
    intervalCenterSo = sollya_lib_mid(rangeSo, None)
2146
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo,
2147
                                         intervalCenterSo,
2148
                                         rangeSo, errorTypeSo, None)
2149
    # Object taylorFormListSaSo is a Python list of Sollya objects references 
2150
    #  that are copies of the elements of taylorFormSo.
2151
    # pobyso_get_list_elements_so_so clears taylorFormSo.
2152
    (taylorFormListSaSo, numElementsSa, isEndEllipticSa) = \
2153
        pobyso_get_list_elements_so_so(taylorFormSo)
2154
    ## Copy needed here since polySo will be returned and taylorFormListSaSo
2155
    #  will be cleared.
2156
    polySo = sollya_lib_copy_obj(taylorFormListSaSo[0])
2157
    #print "Num elements:", numElementsSa
2158
    sollya_lib_clear_obj(taylorFormSo)
2159
    # No copy_obj needed here: a new objects are created.
2160
    maxErrorSo    = sollya_lib_sup(taylorFormListSaSo[2])
2161
    minErrorSo    = sollya_lib_inf(taylorFormListSaSo[2])
2162
    # List taylorFormListSaSo is not needed anymore.
2163
    pobyso_clear_list_elements_sa_so(taylorFormListSaSo)
2164
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
2165
    absMinErrorSo = sollya_lib_abs(minErrorSo)
2166
    sollya_lib_clear_obj(maxErrorSo)
2167
    sollya_lib_clear_obj(minErrorSo)
2168
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
2169
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
2170
    #
2171
    if errorTypeIsNone:
2172
        sollya_lib_clear_obj(errorTypeSo)
2173
    ## If changed, reset the Sollya working precision.
2174
    if sollyaPrecChanged:
2175
        sollya_lib_set_prec(initialSollyaPrecSo)
2176
    sollya_lib_clear_obj(initialSollyaPrecSo)
2177
    ## According to what error is the largest, return the errors.
2178
    if absMaxErrorSa > absMinErrorSa:
2179
        sollya_lib_clear_obj(absMinErrorSo)
2180
        return (polySo, intervalCenterSo, absMaxErrorSo)
2181
    else:
2182
        sollya_lib_clear_obj(absMaxErrorSo)
2183
        return (polySo, intervalCenterSo, absMinErrorSo)
2184
# end pobyso_taylor_expansion_no_change_var_so_so
2185

    
2186
def pobyso_taylor_expansion_with_change_var_so_so(functionSo, degreeSo, \
2187
                                                  rangeSo, \
2188
                                                  errorTypeSo=None, \
2189
                                                  sollyaPrecSo=None):
2190
    """
2191
    Compute the Taylor expansion with the variable change
2192
    x -> (x-intervalCenter) included.
2193
    """
2194
    # Change Sollya internal precision, if need.
2195
    sollyaPrecChanged = False
2196
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
2197
    if sollyaPrecSo is None:
2198
        pass
2199
    else:
2200
        sollya_lib_set_prec(sollyaPrecSo)
2201
        sollyaPrecChanged = True
2202
    #
2203
    # Error type stuff: default to absolute.
2204
    if errorTypeSo is None:
2205
        errorTypeIsNone = True
2206
        errorTypeSo = sollya_lib_absolute(None)
2207
    else:
2208
        errorTypeIsNone = False
2209
    intervalCenterSo = sollya_lib_mid(rangeSo)
2210
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo, \
2211
                                         intervalCenterSo, \
2212
                                         rangeSo, errorTypeSo, None)
2213
    # Object taylorFormListSaSo is a Python list of Sollya objects references 
2214
    # that are copies of the elements of taylorFormSo.
2215
    # pobyso_get_list_elements_so_so clears taylorFormSo.
2216
    (taylorFormListSaSo, numElements, isEndElliptic) = \
2217
        pobyso_get_list_elements_so_so(taylorFormSo)
2218
    sollya_lib_clear_obj(taylorFormSo)
2219
    polySo = taylorFormListSaSo[0]
2220
    ## Maximum error computation with taylorFormListSaSo[2], a range
2221
    #  holding the actual error. Bounds can be negative.
2222
    maxErrorSo    = sollya_lib_sup(taylorFormListSaSo[2])
2223
    minErrorSo    = sollya_lib_inf(taylorFormListSaSo[2])
2224
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
2225
    absMinErrorSo = sollya_lib_abs(minErrorSo)
2226
    sollya_lib_clear_obj(maxErrorSo)
2227
    sollya_lib_clear_obj(minErrorSo)
2228
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
2229
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
2230
    changeVarExpSo = sollya_lib_build_function_sub(\
2231
                       sollya_lib_build_function_free_variable(),\
2232
                       sollya_lib_copy_obj(intervalCenterSo))
2233
    polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpSo)
2234
    # List taylorFormListSaSo is not needed anymore.
2235
    pobyso_clear_list_elements_sa_so(taylorFormListSaSo)
2236
    sollya_lib_clear_obj(changeVarExpSo)
2237
    # If changed, reset the Sollya working precision.
2238
    if sollyaPrecChanged:
2239
        sollya_lib_set_prec(initialSollyaPrecSo)
2240
    sollya_lib_clear_obj(initialSollyaPrecSo)
2241
    if errorTypeIsNone:
2242
        sollya_lib_clear_obj(errorTypeSo)
2243
    # Do not clear maxErrorSo.
2244
    if absMaxErrorSa > absMinErrorSa:
2245
        sollya_lib_clear_obj(absMinErrorSo)
2246
        return((polyVarChangedSo, intervalCenterSo, absMaxErrorSo))
2247
    else:
2248
        sollya_lib_clear_obj(absMaxErrorSo)
2249
        return((polyVarChangedSo, intervalCenterSo, absMinErrorSo))
2250
# end pobyso_taylor_expansion_with_change_var_so_so
2251

    
2252
def pobyso_taylor(function, degree, point):
2253
    """ Legacy function. See pobysoTaylor_so_so. """
2254
    return(pobyso_taylor_so_so(function, degree, point))
2255

    
2256
def pobyso_taylor_so_so(functionSo, degreeSo, pointSo):
2257
    return(sollya_lib_taylor(functionSo, degreeSo, pointSo))
2258
    
2259
def pobyso_taylorform(function, degree, point = None, 
2260
                      interval = None, errorType=None):
2261
    """ Legacy function. See pobyso_taylorform_sa_sa;"""
2262
    
2263
def pobyso_taylorform_sa_sa(functionSa, \
2264
                            degreeSa, \
2265
                            pointSa, \
2266
                            intervalSa=None, \
2267
                            errorTypeSa=None, \
2268
                            precisionSa=None):
2269
    """
2270
    Compute the Taylor form of 'degreeSa' for 'functionSa' at 'pointSa' 
2271
    for 'intervalSa' with 'errorTypeSa' (a string) using 'precisionSa'. 
2272
    point: must be a Real or a Real interval.
2273
    return the Taylor form as an array
2274
    TODO: take care of the interval and of the point when it is an interval;
2275
          when errorType is not None;
2276
          take care of the other elements of the Taylor form (coefficients 
2277
          errors and delta.
2278
    """
2279
    # Absolute as the default error.
2280
    if errorTypeSa is None:
2281
        errorTypeSo = sollya_lib_absolute()
2282
    elif errorTypeSa == "relative":
2283
        errorTypeSo = sollya_lib_relative()
2284
    elif errortypeSa == "absolute":
2285
        errorTypeSo = sollya_lib_absolute()
2286
    else:
2287
        # No clean up needed.
2288
        return None
2289
    # Global precision stuff
2290
    sollyaPrecisionChangedSa = False
2291
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
2292
    if precisionSa is None:
2293
        precSa = initialSollyaPrecSa
2294
    else:
2295
        if precSa > initialSollyaPrecSa:
2296
            if precSa <= 2:
2297
                print inspect.stack()[0][3], ":precision change <= 2 requested."
2298
            pobyso_set_prec_sa_so(precSa)
2299
            sollyaPrecisionChangedSa = True
2300
    #        
2301
    if len(functionSa.variables()) > 0:
2302
        varSa = functionSa.variables()[0]
2303
        pobyso_name_free_variable_sa_so(str(varSa))
2304
    # In any case (point or interval) the parent of pointSa has a precision
2305
    # method.
2306
    pointPrecSa = pointSa.parent().precision()
2307
    if precSa > pointPrecSa:
2308
        pointPrecSa = precSa
2309
    # In any case (point or interval) pointSa has a base_ring() method.
2310
    pointBaseRingString = str(pointSa.base_ring())
2311
    if re.search('Interval', pointBaseRingString) is None: # Point
2312
        pointSo = pobyso_constant_sa_so(pointSa, pointPrecSa)
2313
    else: # Interval.
2314
        pointSo = pobyso_interval_to_range_sa_so(pointSa, pointPrecSa)
2315
    # Sollyafy the function.
2316
    functionSo = pobyso_parse_string_sa_so(functionSa._assume_str().replace('_SAGE_VAR_', ''))
2317
    if sollya_lib_obj_is_error(functionSo):
2318
        print "pobyso_tailorform: function string can't be parsed!"
2319
        return None
2320
    # Sollyafy the degree
2321
    degreeSo = sollya_lib_constant_from_int(int(degreeSa))
2322
    # Sollyafy the point
2323
    # Call Sollya
2324
    taylorFormSo = \
2325
        sollya_lib_taylorform(functionSo, degreeSo, pointSo, errorTypeSo,\
2326
                                         None)
2327
    sollya_lib_clear_obj(functionSo)
2328
    sollya_lib_clear_obj(degreeSo)
2329
    sollya_lib_clear_obj(pointSo)
2330
    sollya_lib_clear_obj(errorTypeSo)
2331
    (tfsAsList, numElements, isEndElliptic) = \
2332
            pobyso_get_list_elements_so_so(taylorFormSo)
2333
    polySo = tfsAsList[0]
2334
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
2335
    polyRealField = RealField(maxPrecision)
2336
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, polyRealField)
2337
    if sollyaPrecisionChangedSa:
2338
        sollya_lib_set_prec(initialSollyaPrecSo)
2339
    sollya_lib_clear_obj(initialSollyaPrecSo)
2340
    polynomialRing = polyRealField[str(varSa)]
2341
    polySa = polynomial(expSa, polynomialRing)
2342
    taylorFormSa = [polySa]
2343
    # Final clean-up.
2344
    sollya_lib_clear_obj(taylorFormSo)
2345
    return(taylorFormSa)
2346
# End pobyso_taylor_form_sa_sa
2347

    
2348
def pobyso_taylorform_so_so(functionSo, degreeSo, pointSo, intervalSo=None, \
2349
                            errorTypeSo=None):
2350
    createdErrorType = False
2351
    if errorTypeSo is None:
2352
        errorTypeSo = sollya_lib_absolute()
2353
        createdErrorType = True
2354
    else:
2355
        #TODO: deal with the other case.
2356
        pass
2357
    if intervalSo is None:
2358
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
2359
                                         errorTypeSo, None)
2360
    else:
2361
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
2362
                                         intervalSo, errorTypeSo, None)
2363
    if createdErrorType:
2364
        sollya_lib_clear_obj(errorTypeSo)
2365
    return resultSo
2366
        
2367

    
2368
def pobyso_univar_polynomial_print_reverse(polySa):
2369
    """ Legacy function. See pobyso_univar_polynomial_print_reverse_sa_sa. """
2370
    return(pobyso_univar_polynomial_print_reverse_sa_sa(polySa))
2371

    
2372
def pobyso_univar_polynomial_print_reverse_sa_sa(polySa):
2373
    """
2374
    Return the string representation of a univariate polynomial with
2375
    monomials ordered in the x^0..x^n order of the monomials.
2376
    Remember: Sage
2377
    """
2378
    polynomialRing = polySa.base_ring()
2379
    # A very expensive solution:
2380
    # -create a fake multivariate polynomial field with only one variable,
2381
    #   specifying a negative lexicographical order;
2382
    mpolynomialRing = PolynomialRing(polynomialRing.base(), \
2383
                                     polynomialRing.variable_name(), \
2384
                                     1, order='neglex')
2385
    # - convert the univariate argument polynomial into a multivariate
2386
    #   version;
2387
    p = mpolynomialRing(polySa)
2388
    # - return the string representation of the converted form.
2389
    # There is no simple str() method defined for p's class.
2390
    return(p.__str__())
2391
#
2392
#print pobyso_get_prec()  
2393
pobyso_set_prec(165)
2394
#print pobyso_get_prec()  
2395
#a=100
2396
#print type(a)
2397
#id(a)
2398
#print "Max arity: ", pobyso_max_arity
2399
#print "Function tripleDouble (43) as a string: ", pobyso_function_type_as_string(43)
2400
#print "Function None (44) as a string: ", pobyso_function_type_as_string(44)
2401
sys.stderr.write("\t...Pobyso check done.\n")