Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py

Historique | Voir | Annoter | Télécharger (103,24 ko)

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

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

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

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

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

31 10 storres
ToDo (among other things):
32 10 storres
 -memory management.
33 5 storres
"""
34 5 storres
from ctypes import *
35 37 storres
import re
36 37 storres
from sage.symbolic.expression_conversions import polynomial
37 59 storres
from sage.symbolic.expression_conversions import PolynomialConverter
38 38 storres
"""
39 284 storres
Create the equivalent to an enum for the Sollya function types and other
40 284 storres
common types (error types...).
41 38 storres
"""
42 5 storres
(SOLLYA_BASE_FUNC_ABS,
43 284 storres
    SOLLYA_BASE_FUNC_ACOS,
44 5 storres
    SOLLYA_BASE_FUNC_ACOSH,
45 5 storres
    SOLLYA_BASE_FUNC_ADD,
46 5 storres
    SOLLYA_BASE_FUNC_ASIN,
47 5 storres
    SOLLYA_BASE_FUNC_ASINH,
48 5 storres
    SOLLYA_BASE_FUNC_ATAN,
49 5 storres
    SOLLYA_BASE_FUNC_ATANH,
50 5 storres
    SOLLYA_BASE_FUNC_CEIL,
51 5 storres
    SOLLYA_BASE_FUNC_CONSTANT,
52 5 storres
    SOLLYA_BASE_FUNC_COS,
53 5 storres
    SOLLYA_BASE_FUNC_COSH,
54 5 storres
    SOLLYA_BASE_FUNC_DIV,
55 5 storres
    SOLLYA_BASE_FUNC_DOUBLE,
56 5 storres
    SOLLYA_BASE_FUNC_DOUBLEDOUBLE,
57 5 storres
    SOLLYA_BASE_FUNC_DOUBLEEXTENDED,
58 5 storres
    SOLLYA_BASE_FUNC_ERF,
59 5 storres
    SOLLYA_BASE_FUNC_ERFC,
60 5 storres
    SOLLYA_BASE_FUNC_EXP,
61 5 storres
    SOLLYA_BASE_FUNC_EXP_M1,
62 5 storres
    SOLLYA_BASE_FUNC_FLOOR,
63 5 storres
    SOLLYA_BASE_FUNC_FREE_VARIABLE,
64 5 storres
    SOLLYA_BASE_FUNC_HALFPRECISION,
65 5 storres
    SOLLYA_BASE_FUNC_LIBRARYCONSTANT,
66 5 storres
    SOLLYA_BASE_FUNC_LIBRARYFUNCTION,
67 5 storres
    SOLLYA_BASE_FUNC_LOG,
68 5 storres
    SOLLYA_BASE_FUNC_LOG_10,
69 5 storres
    SOLLYA_BASE_FUNC_LOG_1P,
70 5 storres
    SOLLYA_BASE_FUNC_LOG_2,
71 5 storres
    SOLLYA_BASE_FUNC_MUL,
72 5 storres
    SOLLYA_BASE_FUNC_NEARESTINT,
73 5 storres
    SOLLYA_BASE_FUNC_NEG,
74 5 storres
    SOLLYA_BASE_FUNC_PI,
75 5 storres
    SOLLYA_BASE_FUNC_POW,
76 5 storres
    SOLLYA_BASE_FUNC_PROCEDUREFUNCTION,
77 5 storres
    SOLLYA_BASE_FUNC_QUAD,
78 5 storres
    SOLLYA_BASE_FUNC_SIN,
79 5 storres
    SOLLYA_BASE_FUNC_SINGLE,
80 5 storres
    SOLLYA_BASE_FUNC_SINH,
81 5 storres
    SOLLYA_BASE_FUNC_SQRT,
82 5 storres
    SOLLYA_BASE_FUNC_SUB,
83 5 storres
    SOLLYA_BASE_FUNC_TAN,
84 5 storres
    SOLLYA_BASE_FUNC_TANH,
85 284 storres
    SOLLYA_BASE_FUNC_TRIPLEDOUBLE,
86 284 storres
    SOLLYA_ABSOLUTE,
87 284 storres
    SOLLYA_RELATIVE) = map(int,xrange(46))
88 284 storres
sys.stderr.write("SOLLYA_RELATIVE = " + str(SOLLYA_RELATIVE) + "\n")
89 246 storres
sys.stderr.write("Superficial pobyso check...\n")
90 230 storres
#print "First constant - SOLLYA_BASE_FUNC_ABS: ", SOLLYA_BASE_FUNC_ABS
91 230 storres
#print "Last constant  - SOLLYA_BASE_FUNC_TRIPLEDOUBLE: ", SOLLYA_BASE_FUNC_TRIPLEDOUBLE
92 5 storres
93 5 storres
pobyso_max_arity = 9
94 5 storres
95 85 storres
def pobyso_absolute_so_so():
96 259 storres
    """
97 259 storres
    Create an "absolute" Sollya object.
98 259 storres
    """
99 85 storres
    return(sollya_lib_absolute(None))
100 85 storres
101 5 storres
def pobyso_autoprint(arg):
102 218 storres
    sollya_lib_autoprint(arg, None)
103 5 storres
104 38 storres
def pobyso_autoprint_so_so(arg):
105 285 storres
    sollya_lib_autoprint(arg, None)
106 282 storres
107 282 storres
def pobyso_bounds_to_interval_sa_sa(lowerBound, upperBound):
108 282 storres
    """
109 282 storres
    Convert a pair of bounds into an interval (an element of
110 282 storres
    a RealIntervalField).
111 282 storres
    """
112 282 storres
    # Minimal (not bullet-proof) check on bounds.
113 282 storres
    if lowerBound > upperBound:
114 282 storres
        return None
115 282 storres
    # Try to get the maximum precision among the bounds.
116 282 storres
    try:
117 282 storres
        preclb = parent(lowerBound).precision()
118 282 storres
        precub = parent(upperBound).precision()
119 282 storres
        prec = max(preclb, precub)
120 282 storres
    except AttributeError:
121 282 storres
        prec = 53
122 282 storres
    # Create the RealIntervalField and the interval (if possible).
123 282 storres
    theRIF = RealIntervalField(prec)
124 285 storres
    theRR  = RealField(prec)
125 285 storres
    # Try convert the lower bound to a RealField element or fail.
126 282 storres
    try:
127 285 storres
        lowerBoundRF = theRR(lowerBound)
128 282 storres
    except TypeError:
129 285 storres
        try:
130 285 storres
            lowerBoundRF = theRR(str(lowerBound))
131 285 storres
        except TypeError:
132 285 storres
            return None
133 285 storres
    # Try convert the upper bound to a RealField element or fail.
134 285 storres
    try:
135 285 storres
        upperBoundRF = theRR(upperBound)
136 285 storres
    except TypeError:
137 285 storres
        try:
138 285 storres
            upperBoundRF = theRR(str(upperBound))
139 285 storres
        except TypeError:
140 285 storres
            return None
141 285 storres
    return theRIF(lowerBoundRF, upperBoundRF)
142 282 storres
# End pobyso_bounds_to_interval_sa_sa
143 54 storres
144 84 storres
def pobyso_bounds_to_range_sa_so(rnLowerBoundSa, rnUpperBoundSa, \
145 84 storres
                                 precisionSa=None):
146 84 storres
    """
147 84 storres
    Return a Sollya range from to 2 RealField Sage elements.
148 84 storres
    The Sollya range element has a sufficient precision to hold all
149 162 storres
    the digits of the widest of the Sage bounds.
150 84 storres
    """
151 84 storres
    # Sanity check.
152 84 storres
    if rnLowerBoundSa > rnUpperBoundSa:
153 84 storres
        return None
154 115 storres
    # Precision stuff.
155 85 storres
    if precisionSa is None:
156 84 storres
        # Check for the largest precision.
157 84 storres
        lbPrecSa = rnLowerBoundSa.parent().precision()
158 84 storres
        ubPrecSa = rnLowerBoundSa.parent().precision()
159 84 storres
        maxPrecSa = max(lbPrecSa, ubPrecSa)
160 84 storres
    else:
161 84 storres
        maxPrecSa = precisionSa
162 115 storres
    # From Sage to Sollya bounds.
163 162 storres
#    lowerBoundSo = sollya_lib_constant(get_rn_value(rnLowerBoundSa),
164 162 storres
#                                       maxPrecSa)
165 162 storres
    lowerBoundSo = pobyso_constant_sa_so(rnLowerBoundSa,
166 162 storres
                                         maxPrecSa)
167 162 storres
    upperBoundSo = pobyso_constant_sa_so(rnUpperBoundSa,
168 154 storres
                                       maxPrecSa)
169 162 storres
170 115 storres
    # From Sollya bounds to range.
171 84 storres
    rangeSo = sollya_lib_range(lowerBoundSo, upperBoundSo)
172 84 storres
    # Back to original precision.
173 84 storres
    # Clean up
174 84 storres
    sollya_lib_clear_obj(lowerBoundSo)
175 84 storres
    sollya_lib_clear_obj(upperBoundSo)
176 154 storres
    return rangeSo
177 84 storres
# End pobyso_bounds_to_range_sa_so
178 84 storres
179 215 storres
def pobyso_build_end_elliptic_list_so_so(*args):
180 215 storres
    """
181 237 storres
    From Sollya argument objects, create a Sollya end elliptic list.
182 237 storres
    Elements of the list are "eaten" (should not be cleared individually,
183 215 storres
    are cleared when the list is cleared).
184 215 storres
    """
185 215 storres
    if len(args) == 0:
186 280 storres
        ## When called with an empty list, sollya_lib_build_end_elliptic_list,
187 280 storres
        #  produces "error".
188 280 storres
        return sollya_lib_build_list(None)
189 215 storres
    index = 0
190 215 storres
    ## One can not append elements to an elliptic list, prepend only is
191 215 storres
    #  permitted.
192 215 storres
    for argument in reversed(args):
193 215 storres
        if index == 0:
194 215 storres
            listSo = sollya_lib_build_end_elliptic_list(argument, None)
195 215 storres
        else:
196 215 storres
            listSo = sollya_lib_prepend(argument, listSo)
197 215 storres
        index += 1
198 215 storres
    return listSo
199 215 storres
200 215 storres
# End pobyso_build_end_elliptic_list_so_so
201 215 storres
202 54 storres
def pobyso_build_function_sub_so_so(exp1So, exp2So):
203 237 storres
    return sollya_lib_build_function_sub(exp1So, exp2So)
204 54 storres
205 237 storres
def pobyso_build_list_of_ints_sa_so(*args):
206 237 storres
    """
207 237 storres
    Build a Sollya list from Sage integral arguments.
208 237 storres
    """
209 237 storres
    if len(args) == 0:
210 237 storres
        return pobyso_build_list_so_so()
211 237 storres
    ## Make a Sage list of integral constants.
212 237 storres
    intsList = []
213 237 storres
    for intElem in args:
214 237 storres
        intsList.append(pobyso_constant_from_int_sa_so(intElem))
215 237 storres
    return pobyso_build_list_so_so(*intsList)
216 237 storres
217 237 storres
def pobyso_build_list_so_so(*args):
218 237 storres
    """
219 237 storres
    Make a Sollya list out of Sollya objects passed as arguments.
220 237 storres
    If one wants to call it with a list argument, should prepend a "*"
221 237 storres
    before the list variable name.
222 237 storres
    Elements of the list are "eaten" (should not be cleared individually,
223 237 storres
    are cleared when the list is cleared).
224 237 storres
    """
225 237 storres
    if len(args) == 0:
226 237 storres
        ## Called with an empty list produced "error".
227 237 storres
        return sollya_lib_build_list(None)
228 237 storres
    index = 0
229 237 storres
    ## Append the Sollya elements one by one.
230 237 storres
    for elementSo in args:
231 237 storres
        if index == 0:
232 237 storres
            listSo = sollya_lib_build_list(elementSo, None)
233 237 storres
        else:
234 237 storres
            listSo = sollya_lib_append(listSo, elementSo)
235 237 storres
        index += 1
236 237 storres
    return listSo
237 237 storres
# End pobyso_build list_so_so
238 237 storres
239 237 storres
240 85 storres
def pobyso_change_var_in_function_so_so(funcSo, chvarExpSo):
241 54 storres
    """
242 85 storres
    Variable change in a function.
243 85 storres
    """
244 85 storres
    return(sollya_lib_evaluate(funcSo,chvarExpSo))
245 85 storres
# End pobyso_change_var_in_function_so_so
246 85 storres
247 285 storres
def pobyso_chebyshevform_sa_sa(funcSa, degreeSa, intervalSa):
248 285 storres
    """
249 285 storres
    An chebyshevnorm call with Sage arguments that returns Sage objects.
250 285 storres
    To the moment only the polynomial and the delta are returned.
251 285 storres
    """
252 285 storres
    # Check that funcSa is a function.
253 285 storres
    if not \
254 285 storres
     sage.symbolic.callable.is_CallableSymbolicExpressionRing(parent(funcSa)):
255 285 storres
        return None
256 285 storres
    # Check that funcSa is an integer.
257 285 storres
    try:
258 285 storres
        if degreeSa != int(degreeSa):
259 285 storres
            print "degreeSa is not a valid integer."
260 285 storres
            return None
261 285 storres
    except ValueError:
262 285 storres
        print "degreeSa is not a number."
263 285 storres
        return None
264 285 storres
    # Create the Sollya degree...
265 285 storres
    degreeSo = pobyso_constant_from_int_sa_so(degreeSa)
266 285 storres
    # ...and check it.
267 285 storres
    if pobyso_is_error_so_sa(degreeSo):
268 285 storres
        print "Could not correctly convert degreeSa to Sollya."
269 285 storres
        return None
270 285 storres
    # Check that intervalSa is an interval.
271 285 storres
    try:
272 285 storres
        intervalSa.upper()
273 285 storres
    except AttributeError:
274 285 storres
        print "intervalSa is not an interval."
275 285 storres
        return None
276 285 storres
    # Convert the Sage polynomial into a Sollya polynomial.
277 285 storres
    # Convert the Sage function into a Sollya function.
278 285 storres
    funcAsStringSa = funcSa._assume_str().replace('_SAGE_VAR_', '')
279 285 storres
    funcSo = pobyso_parse_string_sa_so(funcAsStringSa)
280 285 storres
    if not pobyso_obj_is_function_so_sa(funcSo):
281 285 storres
        sollya_lib_clear_obj(funcSo)
282 285 storres
        print "The Sollya object created from the function is not a function."
283 285 storres
        return None
284 285 storres
    #pobyso_autoprint(funcSo)
285 285 storres
    # Convert the Sage interval into a Sollya range.
286 285 storres
    rangeSo = pobyso_interval_to_range_sa_so(intervalSa)
287 285 storres
    if not pobyso_is_range_so_sa(rangeSo):
288 285 storres
        sollya_lib_clear_obj(funcSo)
289 285 storres
        sollya_lib_clear_obj(rangeSo)
290 285 storres
        print "The Sollya object created from the interval is not a range."
291 285 storres
        return None
292 285 storres
    #pobyso_autoprint(rangeSo)
293 285 storres
    retValSo = sollya_lib_chebyshevform(funcSo,
294 285 storres
                                        degreeSo,
295 285 storres
                                        rangeSo,
296 285 storres
                                        None)
297 285 storres
    sollya_lib_clear_obj(funcSo)
298 285 storres
    sollya_lib_clear_obj(degreeSo)
299 285 storres
    sollya_lib_clear_obj(rangeSo)
300 285 storres
    if pobyso_is_error_so_sa(retValSo):
301 285 storres
        sollya_lib_clear_obj(retValSo)
302 285 storres
        print "Could not compute the Chebyshev form."
303 285 storres
        return None
304 285 storres
    #pobyso_autoprint(retValSo)
305 285 storres
    ## Get the list of Sollya objects returned by sollya_lib_chebyshevform
306 285 storres
    retValSa = pobyso_get_list_elements_so_so(retValSo)
307 285 storres
    if retValSa[1] != 4:
308 285 storres
        print "Invalid number of elements in the Chebyshev form."
309 285 storres
        sollya_lib_clear_obj(retValSo)
310 285 storres
        return None
311 285 storres
    ## Convert all the element of the retValSa[0] array of Sollya objects
312 285 storres
    #  into their Sage counterparts.
313 285 storres
    polySa = pobyso_float_poly_so_sa(retValSa[0][0])
314 285 storres
    if not pobyso_is_poly_sa_sa(polySa):
315 285 storres
        print "Conversion of the polynomial of the Chebyshev form failed."
316 285 storres
        sollya_lib_clear_obj(retValSo)
317 285 storres
        return None
318 285 storres
    deltaSa  = pobyso_get_interval_from_range_so_sa(retValSa[0][2])
319 285 storres
    if not pobyso_is_interval_sa_sa(deltaSa):
320 285 storres
        print "Conversion of the delta interval of the Chebyshev form failed."
321 285 storres
        sollya_lib_clear_obj(retValSo)
322 285 storres
        return None
323 285 storres
    sollya_lib_clear_obj(retValSo)
324 285 storres
    return(polySa,deltaSa)
325 285 storres
# End pobyso_chebyshevform_sa_sa.
326 285 storres
327 85 storres
def pobyso_chebyshevform_so_so(functionSo, degreeSo, intervalSo):
328 85 storres
    resultSo = sollya_lib_chebyshevform(functionSo, degreeSo, intervalSo)
329 85 storres
    return(resultSo)
330 85 storres
# End pobyso_chebyshevform_so_so.
331 85 storres
332 280 storres
def pobyso_clear_full_list_elements_sa_so(objectListSaSo):
333 280 storres
    """
334 280 storres
    Clear the elements of list created by the
335 280 storres
    pobyso_get_list_elements_so_so function.
336 280 storres
    objectListSaSo is as follows:
337 280 storres
    - objectListSaSo[0]: a list of Sollya objects;
338 280 storres
    - objectListSaSo[1]: the number of elements of the previous list;
339 280 storres
    - objectListSaSo[2]: an integer that if != 0 states that the list is
340 280 storres
                         end-elliptic
341 280 storres
    The objects to clear are the elements of the objectListSaSo[0] list.
342 280 storres
    """
343 280 storres
    for index in xrange(0, objectListSaSo[1]):
344 280 storres
        sollya_lib_clear_obj(objectListSaSo[0][index])
345 280 storres
# End pobyso_clear_full_list_elements_sa_so
346 280 storres
347 280 storres
def pobyso_clear_list_elements_sa_so(objectListSaSo):
348 280 storres
    """
349 280 storres
    Clear the elements of list of references to Sollya objects
350 280 storres
    """
351 280 storres
    for index in xrange(0, len(objectListSaSo)):
352 280 storres
        sollya_lib_clear_obj(objectListSaSo[index])
353 280 storres
# End pobyso_clear_list_elements_sa_so
354 280 storres
355 237 storres
def pobyso_clear_obj(objSo):
356 237 storres
    """
357 237 storres
    Free a Sollya object's memory.
358 237 storres
    Very thin wrapper around sollya_lib_clear_obj().
359 237 storres
    """
360 237 storres
    sollya_lib_clear_obj(objSo)
361 237 storres
# End pobyso_clear_obj.
362 237 storres
363 117 storres
def pobyso_clear_taylorform_sa_so(taylorFormSaSo):
364 117 storres
    """
365 285 storres
    This method is wrapper around pobyso_clear_list_elements_sa_so.
366 280 storres
    It is a legacy method left here since it may be used in existing code
367 280 storres
    where Taylor forms are used as Sage lists obtained by converting
368 280 storres
    Sollya Taylor forms (a list made of:
369 280 storres
    - a polynomial;
370 280 storres
    - a list of intervals enclosing the errors accumulated when computing
371 280 storres
      the polynomial coefficients;
372 280 storres
    - a bound on the approximation error between the polynomial and the
373 280 storres
      function.)
374 280 storres
    A Taylor form directly obtained from pobyso_taylorform_so_so is cleared
375 280 storres
    by sollya_lib_clear_obj.
376 117 storres
    """
377 280 storres
    pobyso_clear_list_elements_sa_so(taylorFormSaSo)
378 117 storres
# End pobyso_clear_taylorform_sa_so
379 117 storres
380 85 storres
def pobyso_cmp(rnArgSa, cteSo):
381 85 storres
    """
382 237 storres
    Deprecated, use pobyso_cmp_sa_so_sa instead.
383 237 storres
    """
384 237 storres
    print "Deprecated, use pobyso_cmp_sa_so_sa instead."
385 237 storres
    return pobyso_cmp_sa_so_sa(rnArgSa, cteSo)
386 237 storres
# End  pobyso_cmp
387 237 storres
388 237 storres
def pobyso_cmp_sa_so_sa(rnArgSa, cteSo):
389 237 storres
    """
390 54 storres
    Compare the MPFR value a RealNumber with that of a Sollya constant.
391 54 storres

392 54 storres
    Get the value of the Sollya constant into a RealNumber and compare
393 54 storres
    using MPFR. Could be optimized by working directly with a mpfr_t
394 54 storres
    for the intermediate number.
395 54 storres
    """
396 115 storres
    # Get the precision of the Sollya constant to build a Sage RealNumber
397 115 storres
    # with enough precision.to hold it.
398 5 storres
    precisionOfCte = c_int(0)
399 5 storres
    # From the Sollya constant, create a local Sage RealNumber.
400 85 storres
    sollya_lib_get_prec_of_constant(precisionOfCte, cteSo)
401 5 storres
    #print "Precision of constant: ", precisionOfCte
402 5 storres
    RRRR = RealField(precisionOfCte.value)
403 85 storres
    rnLocalSa = RRRR(0)
404 85 storres
    sollya_lib_get_constant(get_rn_value(rnLocalSa), cteSo)
405 115 storres
    #
406 237 storres
    ## Compare the Sage RealNumber version of the Sollya constant with rnArg
407 237 storres
    #  through a direct comparison of underlying MPFR numbers.
408 237 storres
    return cmp_rn_value(rnArgSa, rnLocal)
409 237 storres
# End pobyso_smp_sa_so_sa
410 5 storres
411 54 storres
def pobyso_compute_pos_function_abs_val_bounds_sa_sa(funcSa, lowerBoundSa, \
412 54 storres
                                                     upperBoundSa):
413 54 storres
    """
414 85 storres
    TODO: completely rework and test.
415 54 storres
    """
416 85 storres
    pobyso = pobyso_name_free_variable_sa_so(funcSa.variables()[0])
417 159 storres
    funcSo = pobyso_parse_string(funcSa._assume_str().replace('_SAGE_VAR_', ''))
418 54 storres
    rangeSo = pobyso_range_sa_so(lowerBoundSa, upperBoundSa)
419 54 storres
    infnormSo = pobyso_infnorm_so_so(funcSo,rangeSo)
420 83 storres
    # Sollya return the infnorm as an interval.
421 54 storres
    fMaxSa = pobyso_get_interval_from_range_so_sa(infnormSo)
422 54 storres
    # Get the top bound and compute the binade top limit.
423 54 storres
    fMaxUpperBoundSa = fMaxSa.upper()
424 54 storres
    binadeTopLimitSa = 2**ceil(fMaxUpperBoundSa.log2())
425 54 storres
    # Put up together the function to use to compute the lower bound.
426 54 storres
    funcAuxSo = pobyso_parse_string(str(binadeTopLimitSa) +  \
427 159 storres
                                    '-(' + f._assume_str().replace('_SAGE_VAR_', '') + ')')
428 54 storres
    pobyso_autoprint(funcAuxSo)
429 83 storres
    # Clear the Sollya range before a new call to infnorm and issue the call.
430 54 storres
    sollya_lib_clear_obj(infnormSo)
431 54 storres
    infnormSo = pobyso_infnorm_so_so(funcAuxSo,rangeSo)
432 54 storres
    fMinSa = pobyso_get_interval_from_range_so_sa(infnormSo)
433 54 storres
    sollya_lib_clear_obj(infnormSo)
434 85 storres
    fMinLowerBoundSa = binadeTopLimitSa - fMinSa.lower()
435 54 storres
    # Compute the maximum of the precisions of the different bounds.
436 54 storres
    maxPrecSa = max([fMinLowerBoundSa.parent().precision(), \
437 54 storres
                     fMaxUpperBoundSa.parent().precision()])
438 54 storres
    # Create a RealIntervalField and create an interval with the "good" bounds.
439 54 storres
    RRRI = RealIntervalField(maxPrecSa)
440 54 storres
    imageIntervalSa = RRRI(fMinLowerBoundSa, fMaxUpperBoundSa)
441 83 storres
    # Free the unneeded Sollya objects
442 54 storres
    sollya_lib_clear_obj(funcSo)
443 54 storres
    sollya_lib_clear_obj(funcAuxSo)
444 54 storres
    sollya_lib_clear_obj(rangeSo)
445 54 storres
    return(imageIntervalSa)
446 83 storres
# End pobyso_compute_pos_function_abs_val_bounds_sa_sa
447 54 storres
448 215 storres
def pobyso_compute_precision_decay_ratio_function_sa_so():
449 215 storres
    """
450 215 storres
    Compute the precision decay ratio function for polynomial
451 215 storres
    coefficient progressive trucation.
452 215 storres
    """
453 215 storres
    functionText = """
454 215 storres
    proc(deg, a, b, we, wq)
455 215 storres
    {
456 215 storres
      k = we * (exp(x/a)-1) + wq * (b*x)^2 + (1-we-wq) * x;
457 215 storres
      return k/k(d);
458 215 storres
    };
459 215 storres
    """
460 215 storres
    return pobyso_parse_string_sa_so(functionText)
461 215 storres
# End  pobyso_compute_precision_decay_ratio_function.
462 215 storres
463 215 storres
464 5 storres
def pobyso_constant(rnArg):
465 38 storres
    """ Legacy function. See pobyso_constant_sa_so. """
466 38 storres
    return(pobyso_constant_sa_so(rnArg))
467 5 storres
468 84 storres
def pobyso_constant_sa_so(rnArgSa, precisionSa=None):
469 52 storres
    """
470 115 storres
    Create a Sollya constant from a Sage RealNumber.
471 209 storres
    The sollya_lib_constant() function creates a constant
472 209 storres
    with the same precision as the source.
473 52 storres
    """
474 209 storres
    ## Precision stuff. If one wants to change precisions,
475 284 storres
    #  everything takes place in Sage. This only makes
476 209 storres
    #  sense if one wants to reduce the precision.
477 226 storres
    # TODO: revisit precision stuff with new technique.
478 284 storres
    # Check that rnArgSa is a realFiedl element.
479 284 storres
    try:
480 284 storres
        rnArgSa.ulp()
481 284 storres
    except AttributeError:
482 284 storres
        return pobyso_error_so()
483 284 storres
    # If a different precision is wanted modify it.
484 209 storres
    if not precisionSa is None:
485 209 storres
        RRR = RealField(precisionSa)
486 209 storres
        rnArgSa = RRR(rnArgSa)
487 209 storres
    #print rnArgSa, rnArgSa.precision()
488 115 storres
    # Sollya constant creation takes place here.
489 209 storres
    return sollya_lib_constant(get_rn_value(rnArgSa))
490 115 storres
# End pobyso_constant_sa_so
491 115 storres
492 55 storres
def pobyso_constant_0_sa_so():
493 115 storres
    """
494 115 storres
    Obvious.
495 115 storres
    """
496 215 storres
    return pobyso_constant_from_int_sa_so(0)
497 55 storres
498 5 storres
def pobyso_constant_1():
499 115 storres
    """
500 115 storres
    Obvious.
501 115 storres
    Legacy function. See pobyso_constant_so_so.
502 115 storres
    """
503 215 storres
    return pobyso_constant_1_sa_so()
504 5 storres
505 52 storres
def pobyso_constant_1_sa_so():
506 115 storres
    """
507 115 storres
    Obvious.
508 115 storres
    """
509 38 storres
    return(pobyso_constant_from_int_sa_so(1))
510 38 storres
511 5 storres
def pobyso_constant_from_int(anInt):
512 38 storres
    """ Legacy function. See pobyso_constant_from_int_sa_so. """
513 215 storres
    return pobyso_constant_from_int_sa_so(anInt)
514 38 storres
515 38 storres
def pobyso_constant_from_int_sa_so(anInt):
516 115 storres
    """
517 115 storres
    Get a Sollya constant from a Sage int.
518 115 storres
    """
519 215 storres
    return sollya_lib_constant_from_int64(long(anInt))
520 5 storres
521 84 storres
def pobyso_constant_from_int_so_sa(constSo):
522 84 storres
    """
523 117 storres
    Get a Sage int from a Sollya int constant.
524 115 storres
    Usefull for precision or powers in polynomials.
525 84 storres
    """
526 200 storres
    constSa = c_long(0)
527 200 storres
    sollya_lib_get_constant_as_int64(byref(constSa), constSo)
528 215 storres
    return constSa.value
529 84 storres
# End pobyso_constant_from_int_so_sa
530 84 storres
531 209 storres
def pobyso_constant_from_mpq_sa_so(rationalSa):
532 200 storres
    """
533 200 storres
    Make a Sollya constant from Sage rational.
534 209 storres
    The Sollya constant is an unevaluated expression.
535 209 storres
    Hence no precision argument is needed.
536 209 storres
    It is better to leave this way since Sollya has its own
537 209 storres
    optimized evaluation mecanism that tries very hard to
538 209 storres
    return exact values or at least faithful ones.
539 200 storres
    """
540 200 storres
    ratExprSo = \
541 200 storres
        sollya_lib_constant_from_mpq(sgmp_get_rational_value(rationalSa))
542 209 storres
    return ratExprSo
543 200 storres
# End pobyso_constant_from_mpq_sa_so.
544 200 storres
545 209 storres
def pobyso_constant_sollya_prec_sa_so(rnArgSa):
546 209 storres
    """
547 209 storres
    Create a Sollya constant from a Sage RealNumber at the
548 209 storres
    current precision in Sollya.
549 209 storres
    """
550 209 storres
    currentSollyaPrecSa = pobyso_get_prec_so_sa()
551 209 storres
    return pobyso_constant_sa_so(rnArgSa, currentSollyaPrecSa)
552 209 storres
# End pobyso_constant_sollya_prec_sa_so
553 215 storres
554 215 storres
def pobyso_end_elliptic_list_so_sa_so(objectsListSo, intCountSa):
555 215 storres
    """
556 215 storres
    Create a Sollya end elliptic list made of the objectListSo[0] to
557 215 storres
     objectsListSo[intCountSa-1] objects.
558 215 storres
    """
559 215 storres
    return sollya_lib_end_elliptic_list(objectSo, int(intCountSa))
560 215 storres
561 155 storres
def pobyso_error_so():
562 155 storres
    return sollya_lib_error(None)
563 155 storres
# End pobyso_error().
564 155 storres
565 262 storres
def pobyso_evaluate_so_sa(funcSo, argumentSo):
566 262 storres
    """
567 262 storres
    Evaluates funcSo for arguemntSo through sollya_lib_evaluate() and return
568 262 storres
    the result as a Sage object
569 262 storres
    """
570 262 storres
    evalSo = sollya_lib_evaluate(funcSo, argumentSo)
571 262 storres
    if pobyso_is_error_so_sa(evalSo):
572 262 storres
        return None
573 262 storres
    if pobyso_is_range_so_sa(evalSo):
574 262 storres
        retVal = pobyso_range_to_interval_so_sa(evalSo)
575 262 storres
    else:
576 262 storres
        retVal = pobyso_get_constant_as_rn(evalSo)
577 262 storres
    sollya_lib_clear_obj(evalSo)
578 262 storres
    return retVal
579 262 storres
# End pobyso_evaluate_so_sa.
580 262 storres
581 215 storres
def pobyso_evaluate_so_so(funcSo, argumentSo):
582 215 storres
    """
583 215 storres
    Evaluates funcSo for arguemntSo through sollya_lib_evaluate().
584 215 storres
    """
585 215 storres
    return sollya_lib_evaluate(funcSo, argumentSo)
586 215 storres
# End pobyso_evaluate_so_so.
587 240 storres
#
588 240 storres
def pobyso_diff_so_so(funcSo):
589 240 storres
    """
590 240 storres
    Very thin wrapper around sollya_lib_diff.
591 240 storres
    """
592 240 storres
    ## TODO: add a check to make sure funcSo is a functional expression.
593 240 storres
    return sollya_lib_diff(funcSo)
594 237 storres
595 234 storres
def pobyso_dirty_find_zeros_so_so(funcSo, rangeSo):
596 234 storres
    """
597 234 storres
    Thin wrapper over sollya_lib_dirtyfindzeros()
598 234 storres
    """
599 234 storres
    return sollya_lib_dirtyfindzeros(funcSo, rangeSo)
600 234 storres
# End pobys_dirty_find_zeros
601 215 storres
602 237 storres
def pobyso_dirty_inf_norm_so_so(funcSo, rangeSo, preSo=None):
603 237 storres
    """
604 237 storres
    Thin wrapper around sollya_dirtyinfnorm().
605 237 storres
    """
606 237 storres
    # TODO: manage the precision change.
607 237 storres
608 237 storres
    return sollya_lib_dirtyinfnorm(funcSo, rangeSo)
609 237 storres
# End pobyso_dirty_inf_norm_so_so
610 237 storres
611 261 storres
def pobyso_find_zeros_so_so(funcSo, rangeSo):
612 261 storres
    """
613 261 storres
    Thin wrapper over sollya_lib_findzeros()
614 261 storres
    """
615 261 storres
    return sollya_lib_findzeros(funcSo, rangeSo)
616 261 storres
# End pobys_find_zeros
617 261 storres
618 234 storres
def pobyso_float_list_so_sa(listSo):
619 234 storres
    """
620 240 storres
    Return a Sollya list of floating-point numbers as a Sage list of
621 234 storres
    floating-point numbers.
622 240 storres
    TODO: add a test to make sure that each element of the list is a constant.
623 234 storres
    """
624 234 storres
    listSa   = []
625 239 storres
    ## The function returns none if the list is empty or an error has happened.
626 239 storres
    retVal = pobyso_get_list_elements_so_so(listSo)
627 239 storres
    if retVal is None:
628 239 storres
        return listSa
629 240 storres
    ## Just in case the interface is changed and an empty list is returned
630 240 storres
    #  instead of None.
631 239 storres
    elif len(retVal) == 0:
632 239 storres
        return listSa
633 239 storres
    else:
634 239 storres
        ## Remember pobyso_get_list_elements_so_so returns more information
635 240 storres
        #  than just the elements of the list (# elements, is_elliptic)
636 239 storres
        listSaSo, numElements, isEndElliptic = retVal
637 240 storres
    ## Return an empty list.
638 239 storres
    if numElements == 0:
639 239 storres
        return listSa
640 234 storres
    ## Search first for the maximum precision of the elements
641 234 storres
    maxPrecSa = 0
642 234 storres
    for floatSo in listSaSo:
643 240 storres
        #pobyso_autoprint(floatSo)
644 234 storres
        curPrecSa =  pobyso_get_prec_of_constant_so_sa(floatSo)
645 234 storres
        if curPrecSa > maxPrecSa:
646 234 storres
            maxPrecSa = curPrecSa
647 234 storres
    ##
648 234 storres
    RF = RealField(maxPrecSa)
649 234 storres
    ##
650 234 storres
    for floatSo in listSaSo:
651 234 storres
        listSa.append(pobyso_get_constant_as_rn_with_rf_so_sa(floatSo))
652 234 storres
    return listSa
653 234 storres
# End pobyso_float_list_so_sa
654 234 storres
655 209 storres
def pobyso_float_poly_sa_so(polySa, precSa = None):
656 209 storres
    """
657 209 storres
    Create a Sollya polynomial from a Sage RealField polynomial.
658 209 storres
    """
659 209 storres
    ## TODO: filter arguments.
660 209 storres
    ## Precision. If a precision is given, convert the polynomial
661 209 storres
    #  into the right polynomial field. If not convert it straight
662 209 storres
    #  to Sollya.
663 218 storres
    sollyaPrecChanged = False
664 226 storres
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
665 218 storres
    if precSa is None:
666 209 storres
        precSa = polySa.parent().base_ring().precision()
667 226 storres
    if (precSa > initialSollyaPrecSa):
668 228 storres
        assert precSa >= 2, "Precision change <2 requested"
669 226 storres
        if precSa <= 2:
670 226 storres
            print inspect.stack()[0][3], ": precision change <= 2 requested"
671 218 storres
        precSo = pobyso_constant_from_int(precSa)
672 218 storres
        pobyso_set_prec_so_so(precSo)
673 218 storres
        sollya_lib_clear_obj(precSo)
674 272 storres
        sollyaPrecChanged = True
675 272 storres
    ## Free variable stuff.
676 272 storres
    freeVariableNameChanged = False
677 272 storres
    polyFreeVariableNameSa = \
678 272 storres
        str(polySa.variables()[0])
679 272 storres
    currentFreeVariableNameSa = pobyso_get_free_variable_name_so_sa()
680 272 storres
    if polyFreeVariableNameSa != currentFreeVariableNameSa:
681 272 storres
        #print "Free variable names do not match.", polyFreeVariableNameSa
682 272 storres
        sollya_lib_name_free_variable(polyFreeVariableNameSa)
683 272 storres
        freeVariableNameChanged = True
684 209 storres
    ## Get exponents and coefficients.
685 218 storres
    exponentsSa     = polySa.exponents()
686 218 storres
    coefficientsSa  = polySa.coefficients()
687 209 storres
    ## Build the polynomial.
688 209 storres
    polySo = None
689 213 storres
    for coefficientSa, exponentSa in zip(coefficientsSa, exponentsSa):
690 209 storres
        #print coefficientSa.n(prec=precSa), exponentSa
691 209 storres
        coefficientSo = \
692 209 storres
            pobyso_constant_sa_so(coefficientSa)
693 209 storres
        #pobyso_autoprint(coefficientSo)
694 209 storres
        exponentSo = \
695 209 storres
            pobyso_constant_from_int_sa_so(exponentSa)
696 209 storres
        #pobyso_autoprint(exponentSo)
697 209 storres
        monomialSo = sollya_lib_build_function_pow(
698 209 storres
                       sollya_lib_build_function_free_variable(),
699 209 storres
                       exponentSo)
700 218 storres
        polyTermSo = sollya_lib_build_function_mul(coefficientSo,
701 218 storres
                                                       monomialSo)
702 209 storres
        if polySo is None:
703 218 storres
            polySo = polyTermSo
704 209 storres
        else:
705 209 storres
            polySo = sollya_lib_build_function_add(polySo, polyTermSo)
706 218 storres
    if sollyaPrecChanged:
707 226 storres
        pobyso_set_prec_so_so(initialSollyaPrecSo)
708 272 storres
    sollya_lib_clear_obj(initialSollyaPrecSo)
709 272 storres
    ## Do not set back the free variable name in Sollya to its initial value:
710 272 storres
    #  it will change it back, in the Sollya polynomial, to what it was in the
711 272 storres
    #  first place.
712 209 storres
    return polySo
713 209 storres
# End pobyso_float_poly_sa_so
714 209 storres
715 209 storres
def pobyso_float_poly_so_sa(polySo, realFieldSa=None):
716 209 storres
    """
717 209 storres
    Convert a Sollya polynomial into a Sage floating-point polynomial.
718 209 storres
    If no realField is given, a RealField corresponding to the maximum
719 209 storres
    precision of the coefficients is internally computed.
720 209 storres
    The real field is not returned but can be easily retrieved from
721 209 storres
    the polynomial itself.
722 209 storres
    ALGORITHM:
723 209 storres
    - (optional) compute the RealField of the coefficients;
724 209 storres
    - convert the Sollya expression into a Sage expression;
725 209 storres
    - convert the Sage expression into a Sage polynomial
726 209 storres
    """
727 209 storres
    if realFieldSa is None:
728 209 storres
        expressionPrecSa = pobyso_get_max_prec_of_exp_so_sa(polySo)
729 218 storres
        #print "Maximum precision of Sollya polynomial coefficients:", expressionPrecSa
730 226 storres
        if expressionPrecSa < 2 or expressionPrecSa > 2147483391:
731 226 storres
            print "Maximum degree of expression:", expressionPrecSa
732 209 storres
        realFieldSa      = RealField(expressionPrecSa)
733 209 storres
    #print "Sollya expression before...",
734 209 storres
    #pobyso_autoprint(polySo)
735 209 storres
736 209 storres
    expressionSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo,
737 209 storres
                                                             realFieldSa)
738 218 storres
    #print "...Sollya expression after."
739 209 storres
    #pobyso_autoprint(polySo)
740 209 storres
    polyVariableSa = expressionSa.variables()[0]
741 209 storres
    polyRingSa     = realFieldSa[str(polyVariableSa)]
742 209 storres
    #print polyRingSa
743 209 storres
    # Do not use the polynomial(expressionSa, ring=polyRingSa) form!
744 209 storres
    polynomialSa = polyRingSa(expressionSa)
745 285 storres
    #polyCoeffsListSa = polynomialSa.coefficients()
746 215 storres
    #for coeff in polyCoeffsListSa:
747 215 storres
    #    print coeff.abs().n()
748 209 storres
    return polynomialSa
749 209 storres
# End pobyso_float_poly_so_sa
750 209 storres
751 215 storres
def pobyso_free_variable():
752 215 storres
    """
753 215 storres
    Ultra thin wrapper around the sollya_lib_function_build_free_variable function.
754 215 storres
    """
755 215 storres
    return sollya_lib_build_function_free_variable()
756 209 storres
757 5 storres
def pobyso_function_type_as_string(funcType):
758 38 storres
    """ Legacy function. See pobyso_function_type_as_string_so_sa. """
759 38 storres
    return(pobyso_function_type_as_string_so_sa(funcType))
760 38 storres
761 38 storres
def pobyso_function_type_as_string_so_sa(funcType):
762 38 storres
    """
763 38 storres
    Numeric Sollya function codes -> Sage mathematical function names.
764 38 storres
    Notice that pow -> ^ (a la Sage, not a la Python).
765 38 storres
    """
766 5 storres
    if funcType == SOLLYA_BASE_FUNC_ABS:
767 5 storres
        return "abs"
768 5 storres
    elif funcType == SOLLYA_BASE_FUNC_ACOS:
769 5 storres
        return "arccos"
770 5 storres
    elif funcType == SOLLYA_BASE_FUNC_ACOSH:
771 5 storres
        return "arccosh"
772 5 storres
    elif funcType == SOLLYA_BASE_FUNC_ADD:
773 5 storres
        return "+"
774 5 storres
    elif funcType == SOLLYA_BASE_FUNC_ASIN:
775 5 storres
        return "arcsin"
776 5 storres
    elif funcType == SOLLYA_BASE_FUNC_ASINH:
777 5 storres
        return "arcsinh"
778 5 storres
    elif funcType == SOLLYA_BASE_FUNC_ATAN:
779 5 storres
        return "arctan"
780 5 storres
    elif funcType == SOLLYA_BASE_FUNC_ATANH:
781 5 storres
        return "arctanh"
782 5 storres
    elif funcType == SOLLYA_BASE_FUNC_CEIL:
783 5 storres
        return "ceil"
784 5 storres
    elif funcType == SOLLYA_BASE_FUNC_CONSTANT:
785 5 storres
        return "cte"
786 5 storres
    elif funcType == SOLLYA_BASE_FUNC_COS:
787 5 storres
        return "cos"
788 5 storres
    elif funcType == SOLLYA_BASE_FUNC_COSH:
789 5 storres
        return "cosh"
790 5 storres
    elif funcType == SOLLYA_BASE_FUNC_DIV:
791 5 storres
        return "/"
792 5 storres
    elif funcType == SOLLYA_BASE_FUNC_DOUBLE:
793 5 storres
        return "double"
794 5 storres
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEDOUBLE:
795 5 storres
        return "doubleDouble"
796 5 storres
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEEXTENDED:
797 5 storres
        return "doubleDxtended"
798 5 storres
    elif funcType == SOLLYA_BASE_FUNC_ERF:
799 5 storres
        return "erf"
800 5 storres
    elif funcType == SOLLYA_BASE_FUNC_ERFC:
801 5 storres
        return "erfc"
802 5 storres
    elif funcType == SOLLYA_BASE_FUNC_EXP:
803 5 storres
        return "exp"
804 5 storres
    elif funcType == SOLLYA_BASE_FUNC_EXP_M1:
805 5 storres
        return "expm1"
806 5 storres
    elif funcType == SOLLYA_BASE_FUNC_FLOOR:
807 5 storres
        return "floor"
808 5 storres
    elif funcType == SOLLYA_BASE_FUNC_FREE_VARIABLE:
809 5 storres
        return "freeVariable"
810 5 storres
    elif funcType == SOLLYA_BASE_FUNC_HALFPRECISION:
811 5 storres
        return "halfPrecision"
812 5 storres
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYCONSTANT:
813 5 storres
        return "libraryConstant"
814 5 storres
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYFUNCTION:
815 5 storres
        return "libraryFunction"
816 5 storres
    elif funcType == SOLLYA_BASE_FUNC_LOG:
817 5 storres
        return "log"
818 5 storres
    elif funcType == SOLLYA_BASE_FUNC_LOG_10:
819 5 storres
        return "log10"
820 5 storres
    elif funcType == SOLLYA_BASE_FUNC_LOG_1P:
821 5 storres
        return "log1p"
822 5 storres
    elif funcType == SOLLYA_BASE_FUNC_LOG_2:
823 5 storres
        return "log2"
824 5 storres
    elif funcType == SOLLYA_BASE_FUNC_MUL:
825 5 storres
        return "*"
826 5 storres
    elif funcType == SOLLYA_BASE_FUNC_NEARESTINT:
827 5 storres
        return "round"
828 5 storres
    elif funcType == SOLLYA_BASE_FUNC_NEG:
829 5 storres
        return "__neg__"
830 5 storres
    elif funcType == SOLLYA_BASE_FUNC_PI:
831 5 storres
        return "pi"
832 5 storres
    elif funcType == SOLLYA_BASE_FUNC_POW:
833 5 storres
        return "^"
834 5 storres
    elif funcType == SOLLYA_BASE_FUNC_PROCEDUREFUNCTION:
835 5 storres
        return "procedureFunction"
836 5 storres
    elif funcType == SOLLYA_BASE_FUNC_QUAD:
837 5 storres
        return "quad"
838 5 storres
    elif funcType == SOLLYA_BASE_FUNC_SIN:
839 5 storres
        return "sin"
840 5 storres
    elif funcType == SOLLYA_BASE_FUNC_SINGLE:
841 5 storres
        return "single"
842 5 storres
    elif funcType == SOLLYA_BASE_FUNC_SINH:
843 5 storres
        return "sinh"
844 5 storres
    elif funcType == SOLLYA_BASE_FUNC_SQRT:
845 5 storres
        return "sqrt"
846 5 storres
    elif funcType == SOLLYA_BASE_FUNC_SUB:
847 5 storres
        return "-"
848 5 storres
    elif funcType == SOLLYA_BASE_FUNC_TAN:
849 5 storres
        return "tan"
850 5 storres
    elif funcType == SOLLYA_BASE_FUNC_TANH:
851 5 storres
        return "tanh"
852 5 storres
    elif funcType == SOLLYA_BASE_FUNC_TRIPLEDOUBLE:
853 5 storres
        return "tripleDouble"
854 5 storres
    else:
855 5 storres
        return None
856 5 storres
857 85 storres
def pobyso_get_constant(rnArgSa, constSo):
858 38 storres
    """ Legacy function. See pobyso_get_constant_so_sa. """
859 209 storres
    return pobyso_get_constant_so_sa(rnArgSa, constSo)
860 209 storres
# End pobyso_get_constant
861 38 storres
862 84 storres
def pobyso_get_constant_so_sa(rnArgSa, constSo):
863 52 storres
    """
864 284 storres
    Set the value of rnArgSa to the value of constSo in MPFR_RNDN mode.
865 52 storres
    rnArg must already exist and belong to some RealField.
866 85 storres
    We assume that constSo points to a Sollya constant.
867 52 storres
    """
868 209 storres
    outcome = sollya_lib_get_constant(get_rn_value(rnArgSa), constSo)
869 209 storres
    if outcome == 0: # Failure because constSo is not a constant expression.
870 209 storres
        return None
871 209 storres
    else:
872 209 storres
        return outcome
873 209 storres
# End  pobyso_get_constant_so_sa
874 209 storres
875 57 storres
def pobyso_get_constant_as_rn(ctExpSo):
876 83 storres
    """
877 83 storres
    Legacy function. See pobyso_get_constant_as_rn_so_sa.
878 83 storres
    """
879 57 storres
    return(pobyso_get_constant_as_rn_so_sa(ctExpSo))
880 38 storres
881 56 storres
def pobyso_get_constant_as_rn_so_sa(constExpSo):
882 83 storres
    """
883 83 storres
    Get a Sollya constant as a Sage "real number".
884 83 storres
    The precision of the floating-point number returned is that of the Sollya
885 83 storres
    constant.
886 83 storres
    """
887 218 storres
    #print "Before computing precision of variable..."
888 218 storres
    #pobyso_autoprint(constExpSo)
889 209 storres
    precisionSa  = pobyso_get_prec_of_constant_so_sa(constExpSo)
890 218 storres
    #print "precisionSa:", precisionSa
891 209 storres
    ## If the expression can not be exactly converted, None is returned.
892 209 storres
    #  In this case opt for the Sollya current expression.
893 209 storres
    if precisionSa is None:
894 209 storres
        precisionSa = pobyso_get_prec_so_sa()
895 56 storres
    RRRR = RealField(precisionSa)
896 56 storres
    rnSa = RRRR(0)
897 209 storres
    outcome = sollya_lib_get_constant(get_rn_value(rnSa), constExpSo)
898 209 storres
    if outcome == 0:
899 209 storres
        return None
900 209 storres
    else:
901 209 storres
        return rnSa
902 83 storres
# End pobyso_get_constant_as_rn_so_sa
903 38 storres
904 38 storres
def pobyso_get_constant_as_rn_with_rf(ctExp, realField):
905 83 storres
    """
906 83 storres
    Legacy function. See pobyso_get_constant_as_rn_with_rf_so_sa.
907 83 storres
    """
908 209 storres
    return pobyso_get_constant_as_rn_with_rf_so_sa(ctExp, realField)
909 209 storres
# End pobyso_get_constant_as_rn_with_rf
910 5 storres
911 56 storres
def pobyso_get_constant_as_rn_with_rf_so_sa(ctExpSo, realFieldSa = None):
912 83 storres
    """
913 83 storres
    Get a Sollya constant as a Sage "real number".
914 83 storres
    If no real field is specified, the precision of the floating-point number
915 85 storres
    returned is that of the Sollya constant.
916 83 storres
    Otherwise is is that of the real field. Hence rounding may happen.
917 83 storres
    """
918 56 storres
    if realFieldSa is None:
919 209 storres
        return pobyso_get_constant_as_rn_so_sa(ctExpSo)
920 56 storres
    rnSa = realFieldSa(0)
921 209 storres
    outcome = sollya_lib_get_constant(get_rn_value(rnSa), ctExpSo)
922 209 storres
    if outcome == 0:
923 209 storres
        return None
924 209 storres
    else:
925 209 storres
        return rnSa
926 83 storres
# End pobyso_get_constant_as_rn_with_rf_so_sa
927 38 storres
928 5 storres
def pobyso_get_free_variable_name():
929 83 storres
    """
930 83 storres
    Legacy function. See pobyso_get_free_variable_name_so_sa.
931 83 storres
    """
932 38 storres
    return(pobyso_get_free_variable_name_so_sa())
933 38 storres
934 38 storres
def pobyso_get_free_variable_name_so_sa():
935 209 storres
    return sollya_lib_get_free_variable_name()
936 5 storres
937 38 storres
def pobyso_get_function_arity(expressionSo):
938 83 storres
    """
939 83 storres
    Legacy function. See pobyso_get_function_arity_so_sa.
940 83 storres
    """
941 38 storres
    return(pobyso_get_function_arity_so_sa(expressionSo))
942 38 storres
943 38 storres
def pobyso_get_function_arity_so_sa(expressionSo):
944 5 storres
    arity = c_int(0)
945 38 storres
    sollya_lib_get_function_arity(byref(arity),expressionSo)
946 209 storres
    return int(arity.value)
947 5 storres
948 38 storres
def pobyso_get_head_function(expressionSo):
949 83 storres
    """
950 83 storres
    Legacy function. See pobyso_get_head_function_so_sa.
951 83 storres
    """
952 38 storres
    return(pobyso_get_head_function_so_sa(expressionSo))
953 38 storres
954 38 storres
def pobyso_get_head_function_so_sa(expressionSo):
955 5 storres
    functionType = c_int(0)
956 218 storres
    sollya_lib_get_head_function(byref(functionType), expressionSo)
957 209 storres
    return int(functionType.value)
958 5 storres
959 56 storres
def pobyso_get_interval_from_range_so_sa(soRange, realIntervalFieldSa = None ):
960 53 storres
    """
961 53 storres
    Return the Sage interval corresponding to the Sollya range argument.
962 83 storres
    If no reaIntervalField is passed as an argument, the interval bounds are not
963 56 storres
    rounded: they are elements of RealIntervalField of the "right" precision
964 56 storres
    to hold all the digits.
965 53 storres
    """
966 53 storres
    prec = c_int(0)
967 56 storres
    if realIntervalFieldSa is None:
968 56 storres
        retval = sollya_lib_get_prec_of_range(byref(prec), soRange, None)
969 56 storres
        if retval == 0:
970 209 storres
            return None
971 56 storres
        realIntervalFieldSa = RealIntervalField(prec.value)
972 56 storres
    intervalSa = realIntervalFieldSa(0,0)
973 53 storres
    retval = \
974 53 storres
        sollya_lib_get_interval_from_range(get_interval_value(intervalSa),\
975 53 storres
                                           soRange)
976 53 storres
    if retval == 0:
977 209 storres
        return None
978 209 storres
    return intervalSa
979 56 storres
# End pobyso_get_interval_from_range_so_sa
980 56 storres
981 5 storres
def pobyso_get_list_elements(soObj):
982 38 storres
    """ Legacy function. See pobyso_get_list_elements_so_so. """
983 209 storres
    return pobyso_get_list_elements_so_so(soObj)
984 38 storres
985 117 storres
def pobyso_get_list_elements_so_so(objectListSo):
986 51 storres
    """
987 118 storres
    Get the Sollya list elements as a Sage/Python array of Sollya objects.
988 118 storres

989 118 storres
    INPUT:
990 118 storres
    - objectListSo: a Sollya list of Sollya objects.
991 118 storres

992 118 storres
    OUTPUT:
993 118 storres
    - a Sage/Python tuple made of:
994 118 storres
      - a Sage/Python list of Sollya objects,
995 118 storres
      - a Sage/Python int holding the number of elements,
996 118 storres
      - a Sage/Python int stating (!= 0) that the list is end-elliptic.
997 118 storres
    NOTE::
998 118 storres
        We recover the addresses of the Sollya object from the list of pointers
999 118 storres
        returned by sollya_lib_get_list_elements. The list itself is freed.
1000 118 storres
    TODO::
1001 118 storres
        Figure out what to do with numElements since the number of elements
1002 118 storres
        can easily be recovered from the list itself.
1003 118 storres
        Ditto for isEndElliptic.
1004 51 storres
    """
1005 5 storres
    listAddress = POINTER(c_longlong)()
1006 5 storres
    numElements = c_int(0)
1007 5 storres
    isEndElliptic = c_int(0)
1008 117 storres
    listAsSageList = []
1009 5 storres
    result = sollya_lib_get_list_elements(byref(listAddress),\
1010 54 storres
                                          byref(numElements),\
1011 54 storres
                                          byref(isEndElliptic),\
1012 117 storres
                                          objectListSo)
1013 5 storres
    if result == 0 :
1014 5 storres
        return None
1015 5 storres
    for i in xrange(0, numElements.value, 1):
1016 118 storres
       #listAsSageList.append(sollya_lib_copy_obj(listAddress[i]))
1017 118 storres
       listAsSageList.append(listAddress[i])
1018 117 storres
       # Clear each of the elements returned by Sollya.
1019 118 storres
       #sollya_lib_clear_obj(listAddress[i])
1020 117 storres
    # Free the list itself.
1021 117 storres
    sollya_lib_free(listAddress)
1022 209 storres
    return (listAsSageList, numElements.value, isEndElliptic.value)
1023 5 storres
1024 38 storres
def pobyso_get_max_prec_of_exp(soExp):
1025 38 storres
    """ Legacy function. See pobyso_get_max_prec_of_exp_so_sa. """
1026 209 storres
    return pobyso_get_max_prec_of_exp_so_sa(soExp)
1027 5 storres
1028 85 storres
def pobyso_get_max_prec_of_exp_so_sa(expSo):
1029 38 storres
    """
1030 38 storres
    Get the maximum precision used for the numbers in a Sollya expression.
1031 52 storres

1032 52 storres
    Arguments:
1033 52 storres
    soExp -- a Sollya expression pointer
1034 52 storres
    Return value:
1035 52 storres
    A Python integer
1036 38 storres
    TODO:
1037 38 storres
    - error management;
1038 38 storres
    - correctly deal with numerical type such as DOUBLEEXTENDED.
1039 38 storres
    """
1040 226 storres
    if expSo is None:
1041 226 storres
        print inspect.stack()[0][3], ": expSo is None."
1042 226 storres
        return 0
1043 5 storres
    maxPrecision = 0
1044 52 storres
    minConstPrec = 0
1045 52 storres
    currentConstPrec = 0
1046 226 storres
    #pobyso_autoprint(expSo)
1047 85 storres
    operator = pobyso_get_head_function_so_sa(expSo)
1048 5 storres
    if (operator != SOLLYA_BASE_FUNC_CONSTANT) and \
1049 5 storres
    (operator != SOLLYA_BASE_FUNC_FREE_VARIABLE):
1050 85 storres
        (arity, subexpressions) = pobyso_get_subfunctions_so_sa(expSo)
1051 5 storres
        for i in xrange(arity):
1052 5 storres
            maxPrecisionCandidate = \
1053 38 storres
                pobyso_get_max_prec_of_exp_so_sa(subexpressions[i])
1054 5 storres
            if maxPrecisionCandidate > maxPrecision:
1055 5 storres
                maxPrecision = maxPrecisionCandidate
1056 209 storres
        return maxPrecision
1057 5 storres
    elif operator == SOLLYA_BASE_FUNC_CONSTANT:
1058 85 storres
        #minConstPrec = pobyso_get_min_prec_of_constant_so_sa(expSo)
1059 52 storres
        #currentConstPrec = pobyso_get_min_prec_of_constant_so_sa(soExp)
1060 52 storres
        #print minConstPrec, " - ", currentConstPrec
1061 209 storres
        return pobyso_get_min_prec_of_constant_so_sa(expSo)
1062 52 storres
1063 5 storres
    elif operator == SOLLYA_BASE_FUNC_FREE_VARIABLE:
1064 209 storres
        return 0
1065 5 storres
    else:
1066 38 storres
        print "pobyso_get_max_prec_of_exp_so_sa: unexepected operator."
1067 209 storres
        return 0
1068 5 storres
1069 85 storres
def pobyso_get_min_prec_of_constant_so_sa(constExpSo):
1070 52 storres
    """
1071 52 storres
    Get the minimum precision necessary to represent the value of a Sollya
1072 52 storres
    constant.
1073 52 storres
    MPFR_MIN_PREC and powers of 2 are taken into account.
1074 209 storres
    We assume that constExpSo is a pointer to a Sollay constant expression.
1075 52 storres
    """
1076 85 storres
    constExpAsRnSa = pobyso_get_constant_as_rn_so_sa(constExpSo)
1077 85 storres
    return(min_mpfr_size(get_rn_value(constExpAsRnSa)))
1078 52 storres
1079 200 storres
def pobyso_get_poly_so_sa(polySo, realFieldSa=None):
1080 200 storres
    """
1081 200 storres
    Convert a Sollya polynomial into a Sage polynomial.
1082 209 storres
    Legacy function. Use pobyso_float_poly_so_sa() instead.
1083 285 storres
    """
1084 285 storres
    print "Deprecated: use pobyso_float_poly_so_sa() instead."
1085 213 storres
    return pobyso_float_poly_so_sa(polySo,realFieldSa)
1086 200 storres
# End pobyso_get_poly_so_sa
1087 200 storres
1088 200 storres
def pobyso_get_prec():
1089 200 storres
    """ Legacy function. See pobyso_get_prec_so_sa(). """
1090 209 storres
    return pobyso_get_prec_so_sa()
1091 200 storres
1092 200 storres
def pobyso_get_prec_so():
1093 200 storres
    """
1094 200 storres
    Get the current default precision in Sollya.
1095 200 storres
    The return value is a Sollya object.
1096 200 storres
    Usefull when modifying the precision back and forth by avoiding
1097 200 storres
    extra conversions.
1098 200 storres
    """
1099 209 storres
    return sollya_lib_get_prec(None)
1100 200 storres
1101 200 storres
def pobyso_get_prec_so_sa():
1102 200 storres
    """
1103 200 storres
    Get the current default precision in Sollya.
1104 200 storres
    The return value is Sage/Python int.
1105 200 storres
    """
1106 227 storres
    precSo = sollya_lib_get_prec()
1107 227 storres
    precSa = pobyso_constant_from_int_so_sa(precSo)
1108 200 storres
    sollya_lib_clear_obj(precSo)
1109 227 storres
    return precSa
1110 200 storres
# End pobyso_get_prec_so_sa.
1111 200 storres
1112 209 storres
def pobyso_get_prec_so_so_sa():
1113 209 storres
    """
1114 209 storres
    Return the current precision both as a Sollya object and a
1115 209 storres
    Sage integer as hybrid tuple.
1116 209 storres
    To avoid multiple calls for precision manipulations.
1117 209 storres
    """
1118 227 storres
    precSo = sollya_lib_get_prec()
1119 227 storres
    precSa = pobyso_constant_from_int_so_sa(precSo)
1120 227 storres
    return (precSo, int(precSa))
1121 200 storres
1122 200 storres
def pobyso_get_prec_of_constant(ctExpSo):
1123 200 storres
    """ Legacy function. See pobyso_get_prec_of_constant_so_sa. """
1124 209 storres
    return pobyso_get_prec_of_constant_so_sa(ctExpSo)
1125 200 storres
1126 200 storres
def pobyso_get_prec_of_constant_so_sa(ctExpSo):
1127 200 storres
    """
1128 200 storres
    Tries to find a precision to represent ctExpSo without rounding.
1129 200 storres
    If not possible, returns None.
1130 200 storres
    """
1131 218 storres
    #print "Entering pobyso_get_prec_of_constant_so_sa..."
1132 200 storres
    prec = c_int(0)
1133 200 storres
    retc = sollya_lib_get_prec_of_constant(byref(prec), ctExpSo, None)
1134 200 storres
    if retc == 0:
1135 218 storres
        #print "pobyso_get_prec_of_constant_so_sa failed."
1136 209 storres
        return None
1137 218 storres
    #print "...exiting pobyso_get_prec_of_constant_so_sa."
1138 209 storres
    return int(prec.value)
1139 200 storres
1140 200 storres
def pobyso_get_prec_of_range_so_sa(rangeSo):
1141 200 storres
    """
1142 200 storres
    Returns the number of bits elements of a range are coded with.
1143 200 storres
    """
1144 200 storres
    prec = c_int(0)
1145 200 storres
    retc = sollya_lib_get_prec_of_range(byref(prec), rangeSo, None)
1146 200 storres
    if retc == 0:
1147 200 storres
        return(None)
1148 209 storres
    return int(prec.value)
1149 200 storres
# End pobyso_get_prec_of_range_so_sa()
1150 200 storres
1151 85 storres
def pobyso_get_sage_exp_from_sollya_exp(sollyaExpSo, realField = RR):
1152 38 storres
    """ Legacy function. See pobyso_get_sage_exp_from_sollya_exp_so_sa. """
1153 209 storres
    return pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo,
1154 209 storres
                                                     realField = RR)
1155 38 storres
1156 85 storres
def pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, realFieldSa = RR):
1157 5 storres
    """
1158 38 storres
    Get a Sage expression from a Sollya expression.
1159 38 storres
    Currently only tested with polynomials with floating-point coefficients.
1160 5 storres
    Notice that, in the returned polynomial, the exponents are RealNumbers.
1161 5 storres
    """
1162 5 storres
    #pobyso_autoprint(sollyaExp)
1163 85 storres
    operatorSa = pobyso_get_head_function_so_sa(sollyaExpSo)
1164 83 storres
    sollyaLibFreeVariableName = sollya_lib_get_free_variable_name()
1165 213 storres
    ## Get rid of the "_"'s in "_x_", if any.
1166 213 storres
    sollyaLibFreeVariableName = re.sub('_', '', sollyaLibFreeVariableName)
1167 5 storres
    # Constants and the free variable are special cases.
1168 5 storres
    # All other operator are dealt with in the same way.
1169 85 storres
    if (operatorSa != SOLLYA_BASE_FUNC_CONSTANT) and \
1170 85 storres
       (operatorSa != SOLLYA_BASE_FUNC_FREE_VARIABLE):
1171 85 storres
        (aritySa, subexpressionsSa) = pobyso_get_subfunctions_so_sa(sollyaExpSo)
1172 85 storres
        if aritySa == 1:
1173 85 storres
            sageExpSa = eval(pobyso_function_type_as_string_so_sa(operatorSa) + \
1174 85 storres
            "(" + pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], \
1175 85 storres
            realFieldSa) + ")")
1176 85 storres
        elif aritySa == 2:
1177 63 storres
            # We do not get through the preprocessor.
1178 63 storres
            # The "^" operator is then a special case.
1179 85 storres
            if operatorSa == SOLLYA_BASE_FUNC_POW:
1180 85 storres
                operatorAsStringSa = "**"
1181 5 storres
            else:
1182 85 storres
                operatorAsStringSa = \
1183 85 storres
                    pobyso_function_type_as_string_so_sa(operatorSa)
1184 85 storres
            sageExpSa = \
1185 85 storres
              eval("pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], realFieldSa)"\
1186 85 storres
              + " " + operatorAsStringSa + " " + \
1187 85 storres
                   "pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[1], realFieldSa)")
1188 63 storres
        # We do not know yet how to deal with arity >= 3
1189 63 storres
        # (is there any in Sollya anyway?).
1190 5 storres
        else:
1191 85 storres
            sageExpSa = eval('None')
1192 209 storres
        return sageExpSa
1193 85 storres
    elif operatorSa == SOLLYA_BASE_FUNC_CONSTANT:
1194 5 storres
        #print "This is a constant"
1195 85 storres
        return pobyso_get_constant_as_rn_with_rf_so_sa(sollyaExpSo, realFieldSa)
1196 85 storres
    elif operatorSa == SOLLYA_BASE_FUNC_FREE_VARIABLE:
1197 218 storres
        #print "This is the free variable"
1198 209 storres
        return eval(sollyaLibFreeVariableName)
1199 5 storres
    else:
1200 5 storres
        print "Unexpected"
1201 5 storres
        return eval('None')
1202 185 storres
# End pobyso_get_sage_exp_from_sollya_exp_so_sa
1203 73 storres
1204 185 storres
1205 38 storres
def pobyso_get_subfunctions(expressionSo):
1206 38 storres
    """ Legacy function. See pobyso_get_subfunctions_so_sa. """
1207 209 storres
    return pobyso_get_subfunctions_so_sa(expressionSo)
1208 200 storres
# End pobyso_get_subfunctions.
1209 200 storres
1210 38 storres
def pobyso_get_subfunctions_so_sa(expressionSo):
1211 38 storres
    """
1212 38 storres
    Get the subfunctions of an expression.
1213 38 storres
    Return the number of subfunctions and the list of subfunctions addresses.
1214 55 storres
    S.T.: Could not figure out another way than that ugly list of declarations
1215 83 storres
    to recover the addresses of the subfunctions.
1216 83 storres
    We limit ourselves to arity 8 functions.
1217 38 storres
    """
1218 5 storres
    subf0 = c_int(0)
1219 5 storres
    subf1 = c_int(0)
1220 5 storres
    subf2 = c_int(0)
1221 5 storres
    subf3 = c_int(0)
1222 5 storres
    subf4 = c_int(0)
1223 5 storres
    subf5 = c_int(0)
1224 5 storres
    subf6 = c_int(0)
1225 5 storres
    subf7 = c_int(0)
1226 5 storres
    subf8 = c_int(0)
1227 5 storres
    arity = c_int(0)
1228 5 storres
    nullPtr = POINTER(c_int)()
1229 38 storres
    sollya_lib_get_subfunctions(expressionSo, byref(arity), \
1230 83 storres
      byref(subf0), byref(subf1), byref(subf2), byref(subf3), \
1231 83 storres
      byref(subf4), byref(subf5),\
1232 83 storres
      byref(subf6), byref(subf7), byref(subf8), nullPtr, None)
1233 83 storres
#    byref(cast(subfunctions[0], POINTER(c_int))), \
1234 83 storres
#    byref(cast(subfunctions[0], POINTER(c_int))), \
1235 83 storres
#    byref(cast(subfunctions[2], POINTER(c_int))), \
1236 83 storres
#    byref(cast(subfunctions[3], POINTER(c_int))), \
1237 83 storres
#    byref(cast(subfunctions[4], POINTER(c_int))), \
1238 83 storres
#    byref(cast(subfunctions[5], POINTER(c_int))), \
1239 83 storres
#    byref(cast(subfunctions[6], POINTER(c_int))), \
1240 83 storres
#    byref(cast(subfunctions[7], POINTER(c_int))), \
1241 5 storres
#    byref(cast(subfunctions[8], POINTER(c_int))), nullPtr)
1242 83 storres
    subfunctions = [subf0, subf1, subf2, subf3, subf4, subf5, subf6, subf7, \
1243 83 storres
                    subf8]
1244 5 storres
    subs = []
1245 5 storres
    if arity.value > pobyso_max_arity:
1246 38 storres
        return(0,[])
1247 5 storres
    for i in xrange(arity.value):
1248 5 storres
        subs.append(int(subfunctions[i].value))
1249 5 storres
        #print subs[i]
1250 209 storres
    return (int(arity.value), subs)
1251 200 storres
# End pobyso_get_subfunctions_so_sa
1252 5 storres
1253 155 storres
def pobyso_guess_degree_sa_sa(functionSa, intervalSa, approxErrorSa,
1254 155 storres
                              weightSa=None, degreeBoundSa=None):
1255 155 storres
    """
1256 155 storres
    Sa_sa variant of the solly_guessdegree function.
1257 155 storres
    Return 0 if something goes wrong.
1258 155 storres
    """
1259 159 storres
    functionAsStringSa = functionSa._assume_str().replace('_SAGE_VAR_', '')
1260 154 storres
    functionSo = pobyso_parse_string_sa_so(functionAsStringSa)
1261 155 storres
    if pobyso_is_error_so_sa(functionSo):
1262 155 storres
        sollya_lib_clear_obj(functionSo)
1263 155 storres
        return 0
1264 154 storres
    rangeSo = pobyso_interval_to_range_sa_so(intervalSa)
1265 155 storres
    # The approximation error is expected to be a floating point number.
1266 155 storres
    if pobyso_is_floating_point_number_sa_sa(approxErrorSa):
1267 155 storres
        approxErrorSo = pobyso_constant_sa_so(approxErrorSa)
1268 155 storres
    else:
1269 155 storres
        approxErrorSo = pobyso_constant_sa_so(RR(approxErrorSa))
1270 154 storres
    if not weightSa is None:
1271 159 storres
        weightAsStringSa = weightSa._assume_str().replace('_SAGE_VAR_', '')
1272 154 storres
        weightSo = pobyso_parse_string_sa_so(weightAsStringSa)
1273 166 storres
        if pobyso_is_error_so_sa(weightSo):
1274 155 storres
            sollya_lib_clear_obj(functionSo)
1275 155 storres
            sollya_lib_clear_obj(rangeSo)
1276 155 storres
            sollya_lib_clear_obj(approxErrorSo)
1277 155 storres
            sollya_lib_clear_obj(weightSo)
1278 155 storres
            return 0
1279 154 storres
    else:
1280 154 storres
        weightSo = None
1281 154 storres
    if not degreeBoundSa is None:
1282 154 storres
        degreeBoundSo = pobyso_constant_from_int_sa_so(degreeBoundSa)
1283 154 storres
    else:
1284 154 storres
        degreeBoundSo = None
1285 154 storres
    guessedDegreeSa = pobyso_guess_degree_so_sa(functionSo,
1286 162 storres
                                                rangeSo,
1287 162 storres
                                                approxErrorSo,
1288 162 storres
                                                weightSo,
1289 162 storres
                                                degreeBoundSo)
1290 154 storres
    sollya_lib_clear_obj(functionSo)
1291 154 storres
    sollya_lib_clear_obj(rangeSo)
1292 155 storres
    sollya_lib_clear_obj(approxErrorSo)
1293 154 storres
    if not weightSo is None:
1294 154 storres
        sollya_lib_clear_obj(weightSo)
1295 154 storres
    if not degreeBoundSo is None:
1296 154 storres
        sollya_lib_clear_obj(degreeBoundSo)
1297 154 storres
    return guessedDegreeSa
1298 154 storres
# End poyso_guess_degree_sa_sa
1299 154 storres
1300 153 storres
def pobyso_guess_degree_so_sa(functionSo, rangeSo, errorSo, weightSo=None, \
1301 154 storres
                              degreeBoundSo=None):
1302 154 storres
    """
1303 154 storres
    Thin wrapper around the guessdegree function.
1304 154 storres
    Nevertheless, some precision control stuff has been appended.
1305 154 storres
    """
1306 154 storres
    # Deal with Sollya internal precision issues: if it is too small,
1307 154 storres
    # compared with the error, increases it to about twice -log2(error).
1308 154 storres
    errorSa = pobyso_get_constant_as_rn_with_rf_so_sa(errorSo)
1309 154 storres
    log2ErrorSa = errorSa.log2()
1310 154 storres
    if log2ErrorSa < 0:
1311 154 storres
        neededPrecisionSa = int(2 * int(-log2ErrorSa) / 64) * 64
1312 154 storres
    else:
1313 154 storres
        neededPrecisionSa = int(2 * int(log2ErrorSa) / 64) * 64
1314 154 storres
    #print "Needed precision:", neededPrecisionSa
1315 226 storres
    sollyaPrecisionHasChanged = False
1316 226 storres
    (initialPrecSo, initialPrecSa) = pobyso_get_prec_so_so_sa()
1317 226 storres
    if neededPrecisionSa > initialPrecSa:
1318 226 storres
        if neededPrecisionSa <= 2:
1319 226 storres
            print inspect.stack()[0][3], ": precision change <= 2 requested."
1320 154 storres
        pobyso_set_prec_sa_so(neededPrecisionSa)
1321 226 storres
        sollyaPrecisionHasChanged = True
1322 166 storres
    #print "Guessing degree..."
1323 153 storres
    # weightSo and degreeBoundsSo are optional arguments.
1324 162 storres
    # As declared, sollya_lib_guessdegree must take 5 arguments.
1325 153 storres
    if weightSo is None:
1326 162 storres
        degreeRangeSo = sollya_lib_guessdegree(functionSo, rangeSo, errorSo,
1327 162 storres
                                               0, 0, None)
1328 154 storres
    elif degreeBoundSo is None:
1329 153 storres
        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, \
1330 162 storres
                                                errorSo, weightSo, 0, None)
1331 153 storres
    else:
1332 153 storres
        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, errorSo, \
1333 154 storres
                                                weightSo, degreeBoundSo, None)
1334 166 storres
    #print "...degree guess done."
1335 154 storres
    # Restore internal precision, if applicable.
1336 226 storres
    if sollyaPrecisionHasChanged:
1337 226 storres
        pobyso_set_prec_so_so(initialPrecSo)
1338 226 storres
        sollya_lib_clear_obj(initialPrecSo)
1339 154 storres
    degreeIntervalSa = pobyso_range_to_interval_so_sa(degreeRangeSo)
1340 154 storres
    sollya_lib_clear_obj(degreeRangeSo)
1341 154 storres
    # When ok, both bounds match.
1342 154 storres
    # When the degree bound is too low, the upper bound is the degree
1343 154 storres
    # for which the error can be honored.
1344 154 storres
    # When it really goes wrong, the upper bound is infinity.
1345 154 storres
    if degreeIntervalSa.lower() == degreeIntervalSa.upper():
1346 154 storres
        return int(degreeIntervalSa.lower())
1347 154 storres
    else:
1348 154 storres
        if degreeIntervalSa.upper().is_infinity():
1349 154 storres
            return None
1350 154 storres
        else:
1351 154 storres
            return int(degreeIntervalSa.upper())
1352 154 storres
    # End pobyso_guess_degree_so_sa
1353 153 storres
1354 215 storres
def pobyso_inf_so_so(intervalSo):
1355 215 storres
    """
1356 215 storres
    Very thin wrapper around sollya_lib_inf().
1357 215 storres
    """
1358 215 storres
    return sollya_lib_inf(intervalSo)
1359 215 storres
# End pobyso_inf_so_so.
1360 282 storres
#
1361 282 storres
def pobyso_infnorm_sa_sa(funcSa, intervalSa):
1362 282 storres
    """
1363 282 storres
    An infnorm call with Sage arguments.
1364 282 storres
    We only take into account the 2 first arguments (the function and
1365 282 storres
    the interval (a range). Managing the other arguments (the file for
1366 282 storres
    the proof and the exclusion intervals list) will be performed later
1367 282 storres
    Changes will be needed in sollya_lib.py file too.
1368 282 storres
    """
1369 282 storres
    # Check that funcSa is a function.
1370 282 storres
    if not \
1371 282 storres
     sage.symbolic.callable.is_CallableSymbolicExpressionRing(parent(funcSa)):
1372 282 storres
        return None
1373 282 storres
    # Check that intervalSa is an interval.
1374 282 storres
    try:
1375 282 storres
        intervalSa.upper()
1376 282 storres
    except AttributeError:
1377 282 storres
        return None
1378 282 storres
    # Convert the Sage function into a Sollya function.
1379 282 storres
    funcAsStringSa = funcSa._assume_str().replace('_SAGE_VAR_', '')
1380 282 storres
    funcSo = pobyso_parse_string_sa_so(funcAsStringSa)
1381 282 storres
    if not pobyso_obj_is_function_so_sa(funcSo):
1382 282 storres
        sollya_lib_clear_obj(funcSo)
1383 282 storres
        return None
1384 282 storres
    # Convert the Sage interval into a Sollya range.
1385 282 storres
    rangeSo = pobyso_interval_to_range_sa_so(intervalSa)
1386 282 storres
    retValSo = sollya_lib_infnorm(funcSo, rangeSo, None)
1387 282 storres
    sollya_lib_clear_obj(funcSo)
1388 282 storres
    sollya_lib_clear_obj(rangeSo)
1389 282 storres
    if pobyso_is_error_so_sa(retValSo):
1390 282 storres
        sollya_lib_clear_obj(retValSo)
1391 282 storres
        return None
1392 282 storres
    retValSa = pobyso_range_to_interval_so_sa(retValSo)
1393 282 storres
    sollya_lib_clear_obj(retValSo)
1394 282 storres
    return retValSa
1395 282 storres
# End pobyso_infnorm_so_so.
1396 282 storres
#
1397 282 storres
def pobyso_infnorm_so_so(funcSo, rangeSo):
1398 282 storres
    """
1399 282 storres
    Very thin wrapper around sollya_lib_infnorm().
1400 282 storres
    We only take into account the 2 first arguments (the function and
1401 282 storres
    the interval (a range). Managing the other arguments (the file for
1402 282 storres
    the proof and the exclusion intervals list) will be performed later
1403 282 storres
    Changes will be needed in sollya_lib.py file too.
1404 215 storres

1405 282 storres
    As per Sollya manual, this function should not be used anymore and
1406 282 storres
    supnorm should be called instead. Nevertheless, supnorm breaks
1407 282 storres
    sometimes whereas infnorm still returns a satisfactory answer.
1408 282 storres
    """
1409 282 storres
    return sollya_lib_infnorm(funcSo, rangeSo, None)
1410 282 storres
# End pobyso_infnorm_so_so.
1411 53 storres
1412 84 storres
def pobyso_interval_to_range_sa_so(intervalSa, precisionSa=None):
1413 84 storres
    if precisionSa is None:
1414 84 storres
        precisionSa = intervalSa.parent().precision()
1415 84 storres
    intervalSo = pobyso_bounds_to_range_sa_so(intervalSa.lower(),\
1416 84 storres
                                              intervalSa.upper(),\
1417 84 storres
                                              precisionSa)
1418 209 storres
    return intervalSo
1419 84 storres
# End pobyso_interval_to_range_sa_so
1420 84 storres
1421 155 storres
def pobyso_is_error_so_sa(objSo):
1422 155 storres
    """
1423 155 storres
    Thin wrapper around the sollya_lib_obj_is_error() function.
1424 155 storres
    """
1425 155 storres
    if sollya_lib_obj_is_error(objSo) != 0:
1426 155 storres
        return True
1427 155 storres
    else:
1428 155 storres
        return False
1429 155 storres
# End pobyso_is_error-so_sa
1430 155 storres
1431 155 storres
def pobyso_is_floating_point_number_sa_sa(numberSa):
1432 155 storres
    """
1433 209 storres
    Check whether a Sage number is floating point.
1434 209 storres
    Exception stuff added because numbers other than
1435 209 storres
    floating-point ones do not have the is_real() attribute.
1436 155 storres
    """
1437 209 storres
    try:
1438 209 storres
        return numberSa.is_real()
1439 209 storres
    except AttributeError:
1440 209 storres
        return False
1441 209 storres
# End pobyso_is_floating_piont_number_sa_sa
1442 155 storres
1443 285 storres
def pobyso_is_interval_sa_sa(objectSa):
1444 285 storres
    """
1445 285 storres
    Check that objectSa is an interval.
1446 285 storres
    """
1447 285 storres
    try:
1448 285 storres
        objectSa.upper()
1449 285 storres
    except AttributeError:
1450 285 storres
        return False
1451 285 storres
    return True
1452 285 storres
1453 285 storres
def pobyso_is_poly_sa_sa(objectSa):
1454 285 storres
    """
1455 285 storres
    Check if an object is a polynomial.
1456 307 storres
    We consider 3 kinds of polynomials:
1457 307 storres
    - symbolicExpressions that are polynomials;
1458 307 storres
    - callable symbolicExpression that are polynomials;
1459 307 storres
    - polynomials over a ring.
1460 285 storres
    """
1461 307 storres
    # We consider a constant as a polynomial. All 3 kinds of polynomials
1462 307 storres
    # have the is_constant() method.
1463 285 storres
    try:
1464 307 storres
        if objectSa.is_constant():
1465 307 storres
            return True
1466 307 storres
    except AttributeError:
1467 307 storres
        return False
1468 307 storres
    # All 3 kinds of polynomials have the args() method.
1469 307 storres
    # A polynomial has at least one argument.
1470 307 storres
    # We catch the TypeError in case the return value of args() has
1471 307 storres
    # no __getitem__ method.
1472 307 storres
    try:
1473 307 storres
        argument = objectSa.args()[0]
1474 307 storres
    except (AttributeError, TypeError):
1475 307 storres
        return False
1476 307 storres
    # The 2 first kinds have a is_polynomial() method that returns True
1477 307 storres
    try:
1478 307 storres
        if not objectSa.is_polynomial(argument):
1479 285 storres
            return False
1480 285 storres
    except AttributeError:
1481 307 storres
        try:
1482 307 storres
            objectSa.polynomial(argument)
1483 307 storres
        except AttributeError:
1484 307 storres
            return False
1485 285 storres
    return True
1486 285 storres
# End pobyso_is_poly_sa_sa
1487 285 storres
1488 226 storres
def pobyso_is_range_so_sa(rangeCandidateSo):
1489 226 storres
    """
1490 226 storres
    Thin wrapper over sollya_lib_is_range.
1491 226 storres
    """
1492 226 storres
    return sollya_lib_obj_is_range(rangeCandidateSo) != 0
1493 226 storres
1494 226 storres
# End pobyso_is_range_so_sa
1495 226 storres
1496 226 storres
1497 37 storres
def pobyso_lib_init():
1498 37 storres
    sollya_lib_init(None)
1499 116 storres
1500 116 storres
def pobyso_lib_close():
1501 116 storres
    sollya_lib_close(None)
1502 237 storres
1503 85 storres
def pobyso_name_free_variable(freeVariableNameSa):
1504 38 storres
    """ Legacy function. See pobyso_name_free_variable_sa_so. """
1505 85 storres
    pobyso_name_free_variable_sa_so(freeVariableNameSa)
1506 38 storres
1507 85 storres
def pobyso_name_free_variable_sa_so(freeVariableNameSa):
1508 83 storres
    """
1509 83 storres
    Set the free variable name in Sollya from a Sage string.
1510 83 storres
    """
1511 85 storres
    sollya_lib_name_free_variable(freeVariableNameSa)
1512 37 storres
1513 281 storres
def pobyso_obj_is_function_so_sa(objSo):
1514 281 storres
    """
1515 281 storres
    Check if an object is a function.
1516 281 storres
    """
1517 281 storres
    if sollya_lib_obj_is_function(objSo) != 0:
1518 281 storres
        return True
1519 281 storres
    else:
1520 281 storres
        return False
1521 281 storres
# End pobyso_obj_is_function_so_sa
1522 281 storres
1523 281 storres
def pobyso_obj_is_range_so_sa(objSo):
1524 281 storres
    """
1525 281 storres
    Check if an object is a function.
1526 281 storres
    """
1527 281 storres
    if sollya_lib_obj_is_range(objSo) != 0:
1528 281 storres
        return True
1529 281 storres
    else:
1530 281 storres
        return False
1531 281 storres
# End pobyso_obj_is_range_so_sa
1532 281 storres
1533 281 storres
def pobyso_obj_is_string_so_sa(objSo):
1534 281 storres
    """
1535 281 storres
    Check if an object is a function.
1536 281 storres
    """
1537 281 storres
    if sollya_lib_obj_is_string(objSo) != 0:
1538 281 storres
        return True
1539 281 storres
    else:
1540 281 storres
        return False
1541 281 storres
# End pobyso_obj_is_string_so_sa
1542 281 storres
1543 5 storres
def pobyso_parse_string(string):
1544 38 storres
    """ Legacy function. See pobyso_parse_string_sa_so. """
1545 209 storres
    return pobyso_parse_string_sa_so(string)
1546 38 storres
1547 38 storres
def pobyso_parse_string_sa_so(string):
1548 83 storres
    """
1549 155 storres
    Get the Sollya expression computed from a Sage string or
1550 155 storres
    a Sollya error object if parsing failed.
1551 83 storres
    """
1552 209 storres
    return sollya_lib_parse_string(string)
1553 5 storres
1554 200 storres
def pobyso_precision_so_sa(ctExpSo):
1555 209 storres
    """
1556 209 storres
    Computes the necessary precision to represent a number.
1557 209 storres
    If x is not zero, it can be uniquely written as x = m · 2e
1558 209 storres
    where m is an odd integer and e is an integer.
1559 209 storres
    precision(x) returns the number of bits necessary to write m
1560 209 storres
    in binary (i.e. ceil(log2(m))).
1561 209 storres
    """
1562 209 storres
    #TODO: take care of the special case: 0, @NaN@, @Inf@
1563 200 storres
    precisionSo = sollya_lib_precision(ctExpSo)
1564 200 storres
    precisionSa = pobyso_constant_from_int_so_sa(precisionSo)
1565 200 storres
    sollya_lib_clear_obj(precisionSo)
1566 200 storres
    return precisionSa
1567 200 storres
# End pobyso_precision_so_sa
1568 215 storres
1569 217 storres
def pobyso_polynomial_coefficients_progressive_round_so_so(polySo,
1570 217 storres
                                                           funcSo,
1571 217 storres
                                                           icSo,
1572 217 storres
                                                           intervalSo,
1573 217 storres
                                                           itpSo,
1574 217 storres
                                                           ftpSo,
1575 217 storres
                                                           maxPrecSo,
1576 219 storres
                                                           maxErrSo,
1577 219 storres
                                                           debug=False):
1578 219 storres
    if debug:
1579 219 storres
        print "Input arguments:"
1580 219 storres
        pobyso_autoprint(polySo)
1581 219 storres
        pobyso_autoprint(funcSo)
1582 219 storres
        pobyso_autoprint(icSo)
1583 219 storres
        pobyso_autoprint(intervalSo)
1584 219 storres
        pobyso_autoprint(itpSo)
1585 219 storres
        pobyso_autoprint(ftpSo)
1586 219 storres
        pobyso_autoprint(maxPrecSo)
1587 219 storres
        pobyso_autoprint(maxErrSo)
1588 219 storres
        print "________________"
1589 200 storres
1590 217 storres
    ## Higher order function see:
1591 217 storres
    #  http://effbot.org/pyfaq/how-do-you-make-a-higher-order-function-in-python.htm
1592 217 storres
    def precision_decay_ratio_function(degreeSa):
1593 217 storres
        def outer(x):
1594 217 storres
            def inner(x):
1595 217 storres
                we = 3/8
1596 217 storres
                wq = 2/8
1597 217 storres
                a  = 2.2
1598 217 storres
                b  = 2
1599 217 storres
                return we*(exp(x/a)-1) +  wq*((b*x)**2) + (1-we-wq)*x
1600 217 storres
            return  inner(x)/inner(degreeSa)
1601 217 storres
        return outer
1602 217 storres
1603 217 storres
    #
1604 217 storres
    degreeSa        = pobyso_polynomial_degree_so_sa(polySo)
1605 217 storres
    ratio           = precision_decay_ratio_function(degreeSa)
1606 217 storres
    itpSa           = pobyso_constant_from_int_so_sa(itpSo)
1607 217 storres
    ftpSa           = pobyso_constant_from_int_so_sa(ftpSo)
1608 217 storres
    maxPrecSa       = pobyso_constant_from_int_so_sa(maxPrecSo)
1609 217 storres
    maxErrSa        = pobyso_get_constant_as_rn_so_sa(maxErrSo)
1610 219 storres
    if debug:
1611 219 storres
        print "degreeSa:", degreeSa
1612 219 storres
        print "ratio:", ratio
1613 219 storres
        print "itpsSa:", itpSa
1614 219 storres
        print "ftpSa:", ftpSa
1615 219 storres
        print "maxPrecSa:", maxPrecSa
1616 219 storres
        print "maxErrSa:", maxErrSa
1617 217 storres
    lastResPolySo   = None
1618 218 storres
    lastInfNormSo   = None
1619 219 storres
    #print "About to enter the while loop..."
1620 217 storres
    while True:
1621 218 storres
        resPolySo   = pobyso_constant_0_sa_so()
1622 217 storres
        pDeltaSa    = ftpSa - itpSa
1623 217 storres
        for indexSa in reversed(xrange(0,degreeSa+1)):
1624 218 storres
            #print "Index:", indexSa
1625 217 storres
            indexSo = pobyso_constant_from_int_sa_so(indexSa)
1626 217 storres
            #pobyso_autoprint(indexSo)
1627 217 storres
            #print ratio(indexSa)
1628 217 storres
            ctpSa = floor(ftpSa - (pDeltaSa * ratio(indexSa)))
1629 217 storres
            ctpSo = pobyso_constant_from_int_sa_so(ctpSa)
1630 219 storres
            if debug:
1631 219 storres
                print "Index:", indexSa, " - Target precision:",
1632 219 storres
                pobyso_autoprint(ctpSo)
1633 217 storres
            cmonSo  = \
1634 217 storres
                sollya_lib_build_function_mul(sollya_lib_coeff(polySo, indexSo),
1635 217 storres
                                      sollya_lib_build_function_pow( \
1636 217 storres
                                          sollya_lib_build_function_free_variable(), \
1637 217 storres
                                          indexSo))
1638 217 storres
            #pobyso_autoprint(cmonSo)
1639 217 storres
            cmonrSo = pobyso_round_coefficients_single_so_so(cmonSo, ctpSo)
1640 217 storres
            sollya_lib_clear_obj(cmonSo)
1641 217 storres
            #pobyso_autoprint(cmonrSo)
1642 217 storres
            resPolySo = sollya_lib_build_function_add(resPolySo,
1643 217 storres
                                                      cmonrSo)
1644 218 storres
            #pobyso_autoprint(resPolySo)
1645 217 storres
        # End for index
1646 217 storres
        freeVarSo     = sollya_lib_build_function_free_variable()
1647 217 storres
        changeVarSo   = sollya_lib_sub(freeVarSo, icSo)
1648 217 storres
        resPolyCvSo   = sollya_lib_evaluate(resPolySo, changeVarSo)
1649 218 storres
        errFuncSo = sollya_lib_build_function_sub(sollya_lib_copy_obj(funcSo),
1650 218 storres
                                                  resPolyCvSo)
1651 218 storres
        infNormSo = sollya_lib_dirtyinfnorm(errFuncSo, intervalSo)
1652 217 storres
        cerrSa    = pobyso_get_constant_as_rn_so_sa(infNormSo)
1653 219 storres
        if debug:
1654 219 storres
            print "Infnorm (Sollya):",
1655 219 storres
            pobyso_autoprint(infNormSo)
1656 218 storres
        sollya_lib_clear_obj(errFuncSo)
1657 218 storres
        #print "Infnorm  (Sage):", cerrSa
1658 217 storres
        if (cerrSa > maxErrSa):
1659 219 storres
            if debug:
1660 219 storres
                print "Error is too large."
1661 218 storres
            if lastResPolySo is None:
1662 219 storres
                if debug:
1663 219 storres
                    print "Enlarging prec."
1664 217 storres
                ntpSa = floor(ftpSa + ftpSa/50)
1665 217 storres
                ## Can't enlarge (numerical)
1666 217 storres
                if ntpSa == ftpSa:
1667 217 storres
                    sollya_lib_clear_obj(resPolySo)
1668 217 storres
                    return None
1669 217 storres
                ## Can't enlarge (not enough precision left)
1670 217 storres
                if ntpSa > maxPrecSa:
1671 217 storres
                    sollya_lib_clear_obj(resPolySo)
1672 217 storres
                    return None
1673 217 storres
                ftpSa = ntpSa
1674 217 storres
                continue
1675 217 storres
            ## One enlargement took place.
1676 217 storres
            else:
1677 219 storres
                if debug:
1678 219 storres
                    print "Exit with the last before last polynomial."
1679 222 storres
                    print "Precision of highest degree monomial:", itpSa
1680 222 storres
                    print "Precision of constant term          :", ftpSa
1681 217 storres
                sollya_lib_clear_obj(resPolySo)
1682 218 storres
                sollya_lib_clear_obj(infNormSo)
1683 218 storres
                return (lastResPolySo, lastInfNormSo)
1684 218 storres
        # cerrSa <= maxErrSa: scrap more bits, possibly.
1685 217 storres
        else:
1686 219 storres
            if debug:
1687 219 storres
                print "Error is too small"
1688 219 storres
            if cerrSa <= (maxErrSa/2):
1689 219 storres
                if debug:
1690 219 storres
                    print "Shrinking prec."
1691 218 storres
                ntpSa = floor(ftpSa - ftpSa/50)
1692 218 storres
                ## Can't shrink (numerical)
1693 218 storres
                if ntpSa == ftpSa:
1694 218 storres
                    if not lastResPolySo is None:
1695 218 storres
                        sollya_lib_clear_obj(lastResPolySo)
1696 218 storres
                    if not lastInfNormSo is None:
1697 218 storres
                        sollya_lib_clear_obj(lastInfNormSo)
1698 222 storres
                    if debug:
1699 222 storres
                        print "Exit because can't shrink anymore (numerically)."
1700 222 storres
                        print "Precision of highest degree monomial:", itpSa
1701 222 storres
                        print "Precision of constant term          :", ftpSa
1702 218 storres
                    return (resPolySo, infNormSo)
1703 218 storres
                ## Can't shrink (not enough precision left)
1704 218 storres
                if ntpSa <= itpSa:
1705 218 storres
                    if not lastResPolySo is None:
1706 218 storres
                        sollya_lib_clear_obj(lastResPolySo)
1707 218 storres
                    if not lastInfNormSo is None:
1708 218 storres
                        sollya_lib_clear_obj(lastInfNormSo)
1709 222 storres
                        print "Exit because can't shrink anymore (no bits left)."
1710 222 storres
                        print "Precision of highest degree monomial:", itpSa
1711 222 storres
                        print "Precision of constant term          :", ftpSa
1712 218 storres
                    return (resPolySo, infNormSo)
1713 218 storres
                ftpSa = ntpSa
1714 217 storres
                if not lastResPolySo is None:
1715 217 storres
                    sollya_lib_clear_obj(lastResPolySo)
1716 218 storres
                if not lastInfNormSo is None:
1717 218 storres
                    sollya_lib_clear_obj(lastInfNormSo)
1718 218 storres
                lastResPolySo = resPolySo
1719 218 storres
                lastInfNormSo = infNormSo
1720 218 storres
                continue
1721 218 storres
            else: # Error is not that small, just return
1722 217 storres
                if not lastResPolySo is None:
1723 217 storres
                    sollya_lib_clear_obj(lastResPolySo)
1724 218 storres
                if not lastInfNormSo is None:
1725 218 storres
                    sollya_lib_clear_obj(lastInfNormSo)
1726 222 storres
                if debug:
1727 222 storres
                    print "Exit normally."
1728 222 storres
                    print "Precision of highest degree monomial:", itpSa
1729 222 storres
                    print "Precision of constant term          :", ftpSa
1730 218 storres
                return (resPolySo, infNormSo)
1731 217 storres
    # End wile True
1732 215 storres
    return None
1733 217 storres
# End pobyso_polynomial_coefficients_progressive_truncate_so_so.
1734 217 storres
1735 217 storres
def pobyso_polynomial_degree_so_sa(polySo):
1736 217 storres
    """
1737 217 storres
    Return the degree of a Sollya polynomial as a Sage int.
1738 217 storres
    """
1739 217 storres
    degreeSo = sollya_lib_degree(polySo)
1740 217 storres
    return pobyso_constant_from_int_so_sa(degreeSo)
1741 217 storres
# End pobyso_polynomial_degree_so_sa
1742 217 storres
1743 217 storres
def pobyso_polynomial_degree_so_so(polySo):
1744 217 storres
    """
1745 217 storres
    Thin wrapper around lib_sollya_degree().
1746 217 storres
    """
1747 217 storres
    return sollya_lib_degree(polySo)
1748 217 storres
# End pobyso_polynomial_degree_so_so
1749 217 storres
1750 5 storres
def pobyso_range(rnLowerBound, rnUpperBound):
1751 38 storres
    """ Legacy function. See pobyso_range_sa_so. """
1752 209 storres
    return pobyso_range_sa_so(rnLowerBound, rnUpperBound)
1753 38 storres
1754 226 storres
def pobyso_range_from_bounds_sa_so(rnLowerBound, rnUpperBound, precSa = None):
1755 226 storres
    """
1756 226 storres
    Create a Sollya range from 2 Sage real numbers as bounds
1757 226 storres
    """
1758 226 storres
    # TODO check precision stuff.
1759 226 storres
    sollyaPrecChanged = False
1760 226 storres
    (initialSollyaPrecSo, initialSollyaPrecSa) = \
1761 226 storres
        pobyso_get_prec_so_so_sa()
1762 226 storres
    if precSa is None:
1763 226 storres
        precSa = max(rnLowerBound.parent().prec(), rnUpperBound.parent().prec())
1764 226 storres
    if precSa > initialSollyaPrecSa:
1765 226 storres
        if precSa <= 2:
1766 226 storres
            print inspect.stack()[0][3], ": precision change <= 2 requested."
1767 226 storres
        pobyso_set_prec_sa_so(precSa)
1768 226 storres
        sollyaPrecChanged = True
1769 226 storres
    rangeSo = sollya_lib_range_from_bounds(get_rn_value(rnLowerBound),
1770 226 storres
                                           get_rn_value(rnUpperBound))
1771 226 storres
    if sollyaPrecChanged:
1772 226 storres
        pobyso_set_prec_so_so(initialSollyaPrecSo)
1773 234 storres
        sollya_lib_clear_obj(initialSollyaPrecSo)
1774 226 storres
    return rangeSo
1775 226 storres
# End pobyso_range_from_bounds_sa_so
1776 5 storres
1777 261 storres
def pobyso_range_list_so_sa(listSo):
1778 261 storres
    """
1779 261 storres
    Return a Sollya list of ranges as a Sage list of
1780 261 storres
    floating-point intervals.
1781 261 storres
    """
1782 261 storres
    listSa   = []
1783 261 storres
    ## The function returns none if the list is empty or an error has happened.
1784 261 storres
    retVal = pobyso_get_list_elements_so_so(listSo)
1785 261 storres
    if retVal is None:
1786 261 storres
        return listSa
1787 261 storres
    ## Just in case the interface is changed and an empty list is returned
1788 261 storres
    #  instead of None.
1789 261 storres
    elif len(retVal) == 0:
1790 261 storres
        return listSa
1791 261 storres
    else:
1792 261 storres
        ## Remember pobyso_get_list_elements_so_so returns more information
1793 261 storres
        #  than just the elements of the list (# elements, is_elliptic)
1794 261 storres
        listSaSo, numElements, isEndElliptic = retVal
1795 261 storres
    ## Return an empty list.
1796 261 storres
    if numElements == 0:
1797 261 storres
        return listSa
1798 261 storres
    ## Search first for the maximum precision of the elements
1799 261 storres
    maxPrecSa = 0
1800 261 storres
    for rangeSo in listSaSo:
1801 261 storres
        #pobyso_autoprint(floatSo)
1802 261 storres
        curPrecSa =  pobyso_get_prec_of_range_so_sa(rangeSo)
1803 261 storres
        if curPrecSa > maxPrecSa:
1804 261 storres
            maxPrecSa = curPrecSa
1805 261 storres
    ##
1806 261 storres
    intervalField = RealIntervalField(maxPrecSa)
1807 261 storres
    ##
1808 261 storres
    for rangeSo in listSaSo:
1809 261 storres
        listSa.append(pobyso_range_to_interval_so_sa(rangeSo, intervalField))
1810 261 storres
    return listSa
1811 261 storres
# End pobyso_range_list_so_sa
1812 261 storres
1813 226 storres
def pobyso_range_max_abs_so_so(rangeSo):
1814 226 storres
    """
1815 226 storres
    Return, as Sollya constant, the maximum absolute value of a Sollay range.
1816 226 storres
    """
1817 226 storres
    lowerBoundSo = sollya_lib_inf(rangeSo)
1818 226 storres
    upperBoundSo = sollya_lib_sup(rangeSo)
1819 226 storres
    #
1820 226 storres
    lowerBoundSo = sollya_lib_build_function_abs(lowerBoundSo)
1821 226 storres
    upperBoundSo = sollya_lib_build_function_abs(upperBoundSo)
1822 226 storres
    #pobyso_autoprint(lowerBoundSo)
1823 226 storres
    #pobyso_autoprint(upperBoundSo)
1824 226 storres
    #
1825 226 storres
    maxAbsSo = sollya_lib_max(lowerBoundSo, upperBoundSo, None)
1826 226 storres
    #sollya_lib_clear_obj(lowerBoundSo)
1827 226 storres
    #sollya_lib_clear_obj(upperBoundSo)
1828 226 storres
    return maxAbsSo
1829 226 storres
# End pobyso_range_max_abs_so_so
1830 226 storres
1831 85 storres
def pobyso_range_to_interval_so_sa(rangeSo, realIntervalFieldSa = None):
1832 83 storres
    """
1833 83 storres
    Get a Sage interval from a Sollya range.
1834 83 storres
    If no realIntervalField is given as a parameter, the Sage interval
1835 83 storres
    precision is that of the Sollya range.
1836 85 storres
    Otherwise, the precision is that of the realIntervalField. In this case
1837 85 storres
    rounding may happen.
1838 83 storres
    """
1839 85 storres
    if realIntervalFieldSa is None:
1840 56 storres
        precSa = pobyso_get_prec_of_range_so_sa(rangeSo)
1841 85 storres
        realIntervalFieldSa = RealIntervalField(precSa)
1842 56 storres
    intervalSa = \
1843 85 storres
        pobyso_get_interval_from_range_so_sa(rangeSo, realIntervalFieldSa)
1844 209 storres
    return intervalSa
1845 209 storres
# End pobyso_range_to_interval_so_sa
1846 240 storres
#
1847 240 storres
def pobyso_relative_so_so():
1848 240 storres
    """
1849 240 storres
    Very thin wrapper around the sollya_lib_relative function.
1850 240 storres
    """
1851 240 storres
    return sollya_lib_relative()
1852 240 storres
# End pobyso_relative_so_so
1853 240 storres
#
1854 209 storres
def pobyso_rat_poly_sa_so(polySa, precSa = None):
1855 209 storres
    """
1856 209 storres
    Create a Sollya polynomial from a Sage rational polynomial.
1857 272 storres
    We first convert the rational polynomial into a floating-point
1858 272 storres
    polynomial.
1859 209 storres
    """
1860 209 storres
    ## TODO: filter arguments.
1861 209 storres
    ## Precision. If no precision is given, use the current precision
1862 209 storres
    #  of Sollya.
1863 209 storres
    if precSa is None:
1864 209 storres
        precSa =  pobyso_get_prec_so_sa()
1865 209 storres
    #print "Precision:",  precSa
1866 209 storres
    RRR = RealField(precSa)
1867 209 storres
    ## Create a Sage polynomial in the "right" precision.
1868 209 storres
    P_RRR = RRR[polySa.variables()[0]]
1869 209 storres
    polyFloatSa = P_RRR(polySa)
1870 213 storres
    ## Make sure no precision is provided: pobyso_float_poly_sa_so will
1871 272 storres
    #  recover it all by itself and will not make any extra conversion.
1872 209 storres
    return pobyso_float_poly_sa_so(polyFloatSa)
1873 209 storres
1874 209 storres
# End pobyso_rat_poly_sa_so
1875 209 storres
1876 52 storres
def pobyso_remez_canonical_sa_sa(func, \
1877 52 storres
                                 degree, \
1878 52 storres
                                 lowerBound, \
1879 52 storres
                                 upperBound, \
1880 52 storres
                                 weight = None, \
1881 52 storres
                                 quality = None):
1882 52 storres
    """
1883 52 storres
    All arguments are Sage/Python.
1884 52 storres
    The functions (func and weight) must be passed as expressions or strings.
1885 52 storres
    Otherwise the function fails.
1886 83 storres
    The return value is a Sage polynomial.
1887 52 storres
    """
1888 83 storres
    var('zorglub')    # Dummy variable name for type check only. Type of
1889 83 storres
    # zorglub is "symbolic expression".
1890 52 storres
    polySo = pobyso_remez_canonical_sa_so(func, \
1891 52 storres
                                 degree, \
1892 52 storres
                                 lowerBound, \
1893 52 storres
                                 upperBound, \
1894 85 storres
                                 weight, \
1895 85 storres
                                 quality)
1896 83 storres
    # String test
1897 52 storres
    if parent(func) == parent("string"):
1898 52 storres
        functionSa = eval(func)
1899 52 storres
    # Expression test.
1900 52 storres
    elif type(func) == type(zorglub):
1901 52 storres
        functionSa = func
1902 83 storres
    else:
1903 83 storres
        return None
1904 83 storres
    #
1905 52 storres
    maxPrecision = 0
1906 52 storres
    if polySo is None:
1907 52 storres
        return(None)
1908 52 storres
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
1909 85 storres
    RRRRSa = RealField(maxPrecision)
1910 85 storres
    polynomialRingSa = RRRRSa[functionSa.variables()[0]]
1911 85 storres
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, RRRRSa)
1912 85 storres
    polySa = polynomial(expSa, polynomialRingSa)
1913 83 storres
    sollya_lib_clear_obj(polySo)
1914 52 storres
    return(polySa)
1915 85 storres
# End pobyso_remez_canonical_sa_sa
1916 52 storres
1917 38 storres
def pobyso_remez_canonical(func, \
1918 5 storres
                           degree, \
1919 5 storres
                           lowerBound, \
1920 5 storres
                           upperBound, \
1921 38 storres
                           weight = "1", \
1922 5 storres
                           quality = None):
1923 38 storres
    """ Legacy function. See pobyso_remez_canonical_sa_so. """
1924 51 storres
    return(pobyso_remez_canonical_sa_so(func, \
1925 51 storres
                                        degree, \
1926 51 storres
                                        lowerBound, \
1927 51 storres
                                        upperBound, \
1928 51 storres
                                        weight, \
1929 51 storres
                                        quality))
1930 200 storres
# End pobyso_remez_canonical.
1931 200 storres
1932 38 storres
def pobyso_remez_canonical_sa_so(func, \
1933 38 storres
                                 degree, \
1934 38 storres
                                 lowerBound, \
1935 38 storres
                                 upperBound, \
1936 52 storres
                                 weight = None, \
1937 38 storres
                                 quality = None):
1938 38 storres
    """
1939 38 storres
    All arguments are Sage/Python.
1940 51 storres
    The functions (func and weight) must be passed as expressions or strings.
1941 51 storres
    Otherwise the function fails.
1942 38 storres
    The return value is a pointer to a Sollya function.
1943 234 storres
    lowerBound and upperBound mus be reals.
1944 38 storres
    """
1945 83 storres
    var('zorglub')    # Dummy variable name for type check only. Type of
1946 83 storres
    # zorglub is "symbolic expression".
1947 85 storres
    currentVariableNameSa = None
1948 52 storres
    # The func argument can be of different types (string,
1949 52 storres
    # symbolic expression...)
1950 38 storres
    if parent(func) == parent("string"):
1951 234 storres
        localFuncSa = sage_eval(func,globals())
1952 85 storres
        if len(localFuncSa.variables()) > 0:
1953 85 storres
            currentVariableNameSa = localFuncSa.variables()[0]
1954 85 storres
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1955 159 storres
            functionSo = \
1956 159 storres
              sollya_lib_parse_string(localFuncSa._assume_str().replace('_SAGE_VAR_', ''))
1957 51 storres
    # Expression test.
1958 52 storres
    elif type(func) == type(zorglub):
1959 52 storres
        # Until we are able to translate Sage expressions into Sollya
1960 52 storres
        # expressions : parse the string version.
1961 85 storres
        if len(func.variables()) > 0:
1962 85 storres
            currentVariableNameSa = func.variables()[0]
1963 85 storres
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1964 159 storres
            functionSo = \
1965 159 storres
              sollya_lib_parse_string(func._assume_str().replace('_SAGE_VAR_', ''))
1966 38 storres
    else:
1967 38 storres
        return(None)
1968 85 storres
    if weight is None: # No weight given -> 1.
1969 52 storres
        weightSo = pobyso_constant_1_sa_so()
1970 85 storres
    elif parent(weight) == parent("string"): # Weight given as string: parse it.
1971 51 storres
        weightSo = sollya_lib_parse_string(func)
1972 85 storres
    elif type(weight) == type(zorglub): # Weight given as symbolice expression.
1973 159 storres
        functionSo = \
1974 159 storres
          sollya_lib_parse_string_sa_so(weight._assume_str().replace('_SAGE_VAR_', ''))
1975 51 storres
    else:
1976 51 storres
        return(None)
1977 5 storres
    degreeSo = pobyso_constant_from_int(degree)
1978 85 storres
    rangeSo = pobyso_bounds_to_range_sa_so(lowerBound, upperBound)
1979 38 storres
    if not quality is None:
1980 38 storres
        qualitySo= pobyso_constant_sa_so(quality)
1981 52 storres
    else:
1982 52 storres
        qualitySo = None
1983 83 storres
1984 83 storres
    remezPolySo = sollya_lib_remez(functionSo, \
1985 83 storres
                                   degreeSo, \
1986 83 storres
                                   rangeSo, \
1987 83 storres
                                   weightSo, \
1988 83 storres
                                   qualitySo, \
1989 83 storres
                                   None)
1990 83 storres
    sollya_lib_clear_obj(functionSo)
1991 83 storres
    sollya_lib_clear_obj(degreeSo)
1992 83 storres
    sollya_lib_clear_obj(rangeSo)
1993 83 storres
    sollya_lib_clear_obj(weightSo)
1994 83 storres
    if not qualitySo is None:
1995 85 storres
        sollya_lib_clear_obj(qualitySo)
1996 83 storres
    return(remezPolySo)
1997 83 storres
# End pobyso_remez_canonical_sa_so
1998 83 storres
1999 38 storres
def pobyso_remez_canonical_so_so(funcSo, \
2000 38 storres
                                 degreeSo, \
2001 38 storres
                                 rangeSo, \
2002 52 storres
                                 weightSo = pobyso_constant_1_sa_so(),\
2003 38 storres
                                 qualitySo = None):
2004 38 storres
    """
2005 38 storres
    All arguments are pointers to Sollya objects.
2006 38 storres
    The return value is a pointer to a Sollya function.
2007 38 storres
    """
2008 38 storres
    if not sollya_lib_obj_is_function(funcSo):
2009 38 storres
        return(None)
2010 38 storres
    return(sollya_lib_remez(funcSo, degreeSo, rangeSo, weightSo, qualitySo, None))
2011 200 storres
# End pobyso_remez_canonical_so_so.
2012 200 storres
2013 237 storres
def pobyso_remez_exponents_list_sa_so(func,     \
2014 237 storres
                                 exponentsList, \
2015 237 storres
                                 lowerBound,    \
2016 237 storres
                                 upperBound,    \
2017 237 storres
                                 weight = None, \
2018 237 storres
                                 quality = None):
2019 237 storres
    """
2020 237 storres
    All arguments are Sage/Python.
2021 237 storres
    The functions (func and weight) must be passed as expressions or strings.
2022 237 storres
    Otherwise the function fails.
2023 237 storres
    The return value is a pointer to a Sollya function.
2024 237 storres
    lowerBound and upperBound mus be reals.
2025 237 storres
    """
2026 237 storres
    var('zorglub')    # Dummy variable name for type check only. Type of
2027 237 storres
    # zorglub is "symbolic expression".
2028 237 storres
    currentVariableNameSa = None
2029 237 storres
    # The func argument can be of different types (string,
2030 237 storres
    # symbolic expression...)
2031 237 storres
    if parent(func) == parent("string"):
2032 237 storres
        localFuncSa = sage_eval(func,globals())
2033 237 storres
        if len(localFuncSa.variables()) > 0:
2034 237 storres
            currentVariableNameSa = localFuncSa.variables()[0]
2035 237 storres
            sollya_lib_name_free_variable(str(currentVariableNameSa))
2036 237 storres
            functionSo = \
2037 237 storres
              sollya_lib_parse_string(localFuncSa._assume_str().replace('_SAGE_VAR_', ''))
2038 237 storres
    # Expression test.
2039 237 storres
    elif type(func) == type(zorglub):
2040 237 storres
        # Until we are able to translate Sage expressions into Sollya
2041 237 storres
        # expressions : parse the string version.
2042 237 storres
        if len(func.variables()) > 0:
2043 237 storres
            currentVariableNameSa = func.variables()[0]
2044 237 storres
            sollya_lib_name_free_variable(str(currentVariableNameSa))
2045 237 storres
            functionSo = \
2046 237 storres
              sollya_lib_parse_string(func._assume_str().replace('_SAGE_VAR_', ''))
2047 237 storres
    else:
2048 237 storres
        return(None)
2049 237 storres
    ## Deal with the weight, much in the same way as with the function.
2050 237 storres
    if weight is None: # No weight given -> 1.
2051 237 storres
        weightSo = pobyso_constant_1_sa_so()
2052 237 storres
    elif parent(weight) == parent("string"): # Weight given as string: parse it.
2053 237 storres
        weightSo = sollya_lib_parse_string(func)
2054 237 storres
    elif type(weight) == type(zorglub): # Weight given as symbolice expression.
2055 237 storres
        functionSo = \
2056 237 storres
          sollya_lib_parse_string(weight._assume_str().replace('_SAGE_VAR_', '',100))
2057 237 storres
    else:
2058 237 storres
        return(None)
2059 237 storres
    rangeSo = pobyso_bounds_to_range_sa_so(lowerBound, upperBound)
2060 237 storres
    if not quality is None:
2061 237 storres
        qualitySo= pobyso_constant_sa_so(quality)
2062 237 storres
    else:
2063 237 storres
        qualitySo = None
2064 237 storres
    #
2065 237 storres
    ## Tranform the Sage list of exponents into a Sollya list.
2066 237 storres
    exponentsListSo = pobyso_build_list_of_ints_sa_so(*exponentsList)
2067 278 storres
    remezPolySo = sollya_lib_remez(functionSo,      \
2068 237 storres
                                   exponentsListSo, \
2069 278 storres
                                   rangeSo,         \
2070 278 storres
                                   weightSo,        \
2071 278 storres
                                   qualitySo,       \
2072 237 storres
                                   None)
2073 237 storres
    sollya_lib_clear_obj(functionSo)
2074 237 storres
    sollya_lib_clear_obj(exponentsListSo)
2075 237 storres
    sollya_lib_clear_obj(rangeSo)
2076 237 storres
    sollya_lib_clear_obj(weightSo)
2077 237 storres
    if not qualitySo is None:
2078 237 storres
        sollya_lib_clear_obj(qualitySo)
2079 237 storres
    return(remezPolySo)
2080 237 storres
# End pobyso_remez_exponentsList_sa_so
2081 240 storres
#
2082 240 storres
def pobyso_round_coefficients_so_so(polySo, truncFormatListSo):
2083 240 storres
    """
2084 240 storres
    A wrapper around the "classical" sollya_lib_roundcoefficients: a Sollya
2085 240 storres
    polynomial and a Sollya list are given as arguments.
2086 240 storres
    """
2087 240 storres
    return sollya_lib_roundcoefficients(polySo, truncFormatListSo)
2088 237 storres
2089 224 storres
def pobyso_round_coefficients_progressive_so_so(polySo,
2090 224 storres
                                                funcSo,
2091 224 storres
                                                precSo,
2092 224 storres
                                                intervalSo,
2093 224 storres
                                                icSo,
2094 224 storres
                                                currentApproxErrorSo,
2095 224 storres
                                                approxAccurSo,
2096 224 storres
                                                debug=False):
2097 272 storres
    """
2098 272 storres
    From an input approximation polynomial, compute an output one with
2099 272 storres
    smaller coefficients and yet yields a sufficient approximation accuracy.
2100 272 storres
    """
2101 223 storres
    if debug:
2102 223 storres
        print "Input arguments:"
2103 228 storres
        print "Polynomial: ", ; pobyso_autoprint(polySo)
2104 228 storres
        print "Function: ", ; pobyso_autoprint(funcSo)
2105 228 storres
        print "Internal precision: ", ; pobyso_autoprint(precSo)
2106 228 storres
        print "Interval: ", ; pobyso_autoprint(intervalSo)
2107 228 storres
        print "Current approximation error: ", ; pobyso_autoprint(currentApproxErrorSo)
2108 228 storres
        print "Requested approxiation error: ", ; pobyso_autoprint(approxAccurSo)
2109 223 storres
        print "________________"
2110 223 storres
    approxAccurSa        = pobyso_get_constant_as_rn_so_sa(approxAccurSo)
2111 223 storres
    currentApproxErrorSa = pobyso_get_constant_as_rn_so_sa(currentApproxErrorSo)
2112 223 storres
    ## If the current approximation error is too close to the target, there is
2113 223 storres
    #  no possible gain.
2114 223 storres
    if currentApproxErrorSa >= approxAccurSa / 2:
2115 232 storres
        #### Do not return the initial argument but copies: caller may free
2116 272 storres
        #    the former as inutile after call.
2117 232 storres
        return (sollya_lib_copy_obj(polySo),
2118 232 storres
                sollya_lib_copy_obj(currentApproxErrorSo))
2119 278 storres
    #
2120 278 storres
    ## Try to round the coefficients.
2121 223 storres
    degreeSa             = pobyso_polynomial_degree_so_sa(polySo)
2122 223 storres
    intervalSa           = pobyso_range_to_interval_so_sa(intervalSo)
2123 223 storres
2124 223 storres
    if debug:
2125 228 storres
        print "degreeSa              :", degreeSa
2126 228 storres
        print "intervalSa            :", intervalSa.str(style='brackets')
2127 223 storres
        print "currentApproxErrorSa  :", currentApproxErrorSa
2128 228 storres
        print "approxAccurSa         :", approxAccurSa
2129 223 storres
    radiusSa = intervalSa.absolute_diameter() / 2
2130 223 storres
    if debug:
2131 224 storres
        print "log2(radius):", RR(radiusSa).log2()
2132 224 storres
    iterIndex = 0
2133 272 storres
    ## Build the "shaved" polynomial.
2134 224 storres
    while True:
2135 272 storres
        ### Start with a 0 value expression.
2136 224 storres
        resPolySo = pobyso_constant_0_sa_so()
2137 224 storres
        roundedPolyApproxAccurSa = approxAccurSa / 2
2138 224 storres
        currentRadiusPowerSa = 1
2139 224 storres
        for degree in xrange(0,degreeSa + 1):
2140 224 storres
            #### At round 0, use the agressive formula. At round 1, run the
2141 224 storres
            #    proved formula.
2142 224 storres
            if iterIndex == 0:
2143 224 storres
                roundingPowerSa = \
2144 224 storres
                    floor(((currentRadiusPowerSa/roundedPolyApproxAccurSa)*(degree+1)).log2())
2145 224 storres
            else:
2146 224 storres
                roundingPowerSa = \
2147 224 storres
                    floor(((currentRadiusPowerSa/roundedPolyApproxAccurSa)*(degreeSa+1)).log2())
2148 272 storres
            ## Under extreme conditions the above formulas can evaluate under 2,
2149 272 storres
            #   which is the minimal precision of an MPFR number.
2150 228 storres
            if roundingPowerSa < 2:
2151 228 storres
                roundingPowerSa = 2
2152 224 storres
            if debug:
2153 224 storres
                print "roundedPolyApproxAccurSa", roundedPolyApproxAccurSa
2154 224 storres
                print "currentRadiusPowerSa", currentRadiusPowerSa
2155 224 storres
                print "Current rounding exponent:", roundingPowerSa
2156 224 storres
            currentRadiusPowerSa *= radiusSa
2157 224 storres
            index1So = pobyso_constant_from_int_sa_so(degree)
2158 224 storres
            index2So = pobyso_constant_from_int_sa_so(degree)
2159 224 storres
            ### Create a monomial with:
2160 224 storres
            #   - the coefficient in the initial monomial at the current degrree;
2161 224 storres
            #   - the current exponent;
2162 224 storres
            #   - the free variable.
2163 224 storres
            cmonSo  = \
2164 224 storres
                sollya_lib_build_function_mul(sollya_lib_coeff(polySo, index1So),
2165 224 storres
                                              sollya_lib_build_function_pow( \
2166 224 storres
                                                sollya_lib_build_function_free_variable(), \
2167 224 storres
                                                index2So))
2168 224 storres
            roundingPowerSo = pobyso_constant_from_int_sa_so(roundingPowerSa)
2169 224 storres
            cmonrSo = pobyso_round_coefficients_single_so_so(cmonSo, roundingPowerSo)
2170 224 storres
            sollya_lib_clear_obj(cmonSo)
2171 224 storres
            ### Add to the result polynomial.
2172 224 storres
            resPolySo = sollya_lib_build_function_add(resPolySo,
2173 224 storres
                                                      cmonrSo)
2174 224 storres
        # End for.
2175 224 storres
        ### Check the new polynomial.
2176 224 storres
        freeVarSo     = sollya_lib_build_function_free_variable()
2177 224 storres
        changeVarSo   = sollya_lib_sub(freeVarSo, icSo)
2178 224 storres
        resPolyCvSo   = sollya_lib_evaluate(resPolySo, changeVarSo)
2179 224 storres
        errFuncSo     = sollya_lib_build_function_sub(sollya_lib_copy_obj(funcSo),
2180 224 storres
                                                      resPolyCvSo)
2181 224 storres
        infNormSo = sollya_lib_dirtyinfnorm(errFuncSo, intervalSo)
2182 224 storres
        ### This also clears resPolyCvSo.
2183 224 storres
        sollya_lib_clear_obj(errFuncSo)
2184 224 storres
        cerrSa    = pobyso_get_constant_as_rn_so_sa(infNormSo)
2185 223 storres
        if debug:
2186 224 storres
            print "Error of the new polynomial:", cerrSa
2187 224 storres
        ### If at round 1, return the initial polynomial error. This should
2188 224 storres
        #   never happen since the rounding algorithm is proved. But some
2189 224 storres
        #   circumstances may break it (e.g. internal precision of tools).
2190 224 storres
        if cerrSa > approxAccurSa:
2191 278 storres
            if iterIndex == 0: # Round 0 is agressive rounding, got round 1 (proved rounding)
2192 224 storres
                sollya_lib_clear_obj(resPolySo)
2193 224 storres
                sollya_lib_clear_obj(infNormSo)
2194 278 storres
                iterIndex += 1
2195 278 storres
                continue
2196 278 storres
            else: # Round 1 and beyond : just return the oroginal polynomial.
2197 278 storres
                sollya_lib_clear_obj(resPolySo)
2198 278 storres
                sollya_lib_clear_obj(infNormSo)
2199 232 storres
                #### Do not return the arguments but copies: the caller may free
2200 232 storres
                #    free the former as inutile after call.
2201 232 storres
                return (sollya_lib_copy_obj(polySo),
2202 232 storres
                        sollya_lib_copy_obj(currentApproxErrorSo))
2203 224 storres
        ### If get here it is because cerrSa <= approxAccurSa
2204 224 storres
        ### Approximation error of the new polynomial is acceptable.
2205 224 storres
        return (resPolySo, infNormSo)
2206 224 storres
    # End while True
2207 224 storres
# End pobyso_round_coefficients_progressive_so_so
2208 223 storres
2209 240 storres
def pobyso_round_coefficients_single_so_so(polySo, commonPrecSo):
2210 215 storres
    """
2211 215 storres
    Create a rounded coefficients polynomial from polynomial argument to
2212 215 storres
    the number of bits in size argument.
2213 215 storres
    All coefficients are set to the same precision.
2214 215 storres
    """
2215 215 storres
    ## TODO: check arguments.
2216 240 storres
    endEllipListSo = pobyso_build_end_elliptic_list_so_so(commonPrecSo)
2217 215 storres
    polySo = sollya_lib_roundcoefficients(polySo, endEllipListSo, None)
2218 217 storres
    sollya_lib_clear_obj(endEllipListSo)
2219 215 storres
    #sollya_lib_clear_obj(endEllipListSo)
2220 215 storres
    return polySo
2221 215 storres
2222 215 storres
# End pobyso_round_coefficients_single_so_so
2223 215 storres
2224 5 storres
def pobyso_set_canonical_off():
2225 5 storres
    sollya_lib_set_canonical(sollya_lib_off())
2226 5 storres
2227 5 storres
def pobyso_set_canonical_on():
2228 5 storres
    sollya_lib_set_canonical(sollya_lib_on())
2229 5 storres
2230 5 storres
def pobyso_set_prec(p):
2231 38 storres
    """ Legacy function. See pobyso_set_prec_sa_so. """
2232 85 storres
    pobyso_set_prec_sa_so(p)
2233 38 storres
2234 38 storres
def pobyso_set_prec_sa_so(p):
2235 227 storres
    #a = c_int(p)
2236 226 storres
    #precSo = c_void_p(sollya_lib_constant_from_int(a))
2237 227 storres
    #precSo = sollya_lib_constant_from_int(a)
2238 227 storres
    precSo = pobyso_constant_from_int_sa_so(p)
2239 226 storres
    sollya_lib_set_prec(precSo)
2240 226 storres
    sollya_lib_clear_obj(precSo)
2241 215 storres
# End pobyso_set_prec_sa_so.
2242 5 storres
2243 85 storres
def pobyso_set_prec_so_so(newPrecSo):
2244 226 storres
    sollya_lib_set_prec(newPrecSo)
2245 215 storres
# End pobyso_set_prec_so_so.
2246 242 storres
#
2247 284 storres
def pobyso_supnorm_sa_sa(polySa,
2248 284 storres
                         funcSa,
2249 284 storres
                         intervalSa,
2250 284 storres
                         errorTypeSa=SOLLYA_ABSOLUTE,
2251 285 storres
                         accuracySa=RR(2**-40)):
2252 281 storres
    """
2253 284 storres
    An supnorm call with Sage arguments.
2254 281 storres
    """
2255 284 storres
    # Check that polySa is a univariate polynomial. We only check here at
2256 284 storres
    # expression level, not at polynomial ring level
2257 284 storres
    ## Make sure it is not a multivariate polynomial. The number of variables
2258 284 storres
    #  can be zero in the case of a constant.
2259 284 storres
    try:
2260 284 storres
        if len(polySa.free_variables()) > 1 :
2261 284 storres
            print "Invalid number of free variables in polySa:", polySa,
2262 284 storres
            " -", polySa.free_variables()
2263 284 storres
            return None
2264 284 storres
    except AttributeError:
2265 284 storres
        print "polySa is not a SymbolicExpression."
2266 284 storres
        return None
2267 284 storres
    ## Make sure it is a polynomial.
2268 284 storres
    if not polySa.is_polynomial(polySa.default_variable()):
2269 284 storres
        print "polySa is not a polynomila expression."
2270 284 storres
        return None
2271 284 storres
    # Check that funcSa is a function.
2272 284 storres
    if not \
2273 284 storres
     sage.symbolic.callable.is_CallableSymbolicExpressionRing(parent(funcSa)):
2274 284 storres
        return None
2275 284 storres
    # Check that intervalSa is an interval.
2276 284 storres
    try:
2277 284 storres
        intervalSa.upper()
2278 284 storres
    except AttributeError:
2279 284 storres
        print "intervalSa is not an interval."
2280 284 storres
        return None
2281 284 storres
    # Convert the Sage polynomial into a Sollya polynomial.
2282 284 storres
    polyAsStringSa = polySa._assume_str().replace('_SAGE_VAR_', '')
2283 284 storres
    polySo = pobyso_parse_string_sa_so(polyAsStringSa)
2284 284 storres
    if not pobyso_obj_is_function_so_sa(polySo):
2285 284 storres
        sollya_lib_clear_obj(polySo)
2286 284 storres
        print "The Sollya object created from the polynomial is not a function."
2287 284 storres
        return None
2288 284 storres
    #pobyso_autoprint(polySo)
2289 284 storres
    # Convert the Sage function into a Sollya function.
2290 284 storres
    funcAsStringSa = funcSa._assume_str().replace('_SAGE_VAR_', '')
2291 284 storres
    funcSo = pobyso_parse_string_sa_so(funcAsStringSa)
2292 284 storres
    if not pobyso_obj_is_function_so_sa(funcSo):
2293 284 storres
        sollya_lib_clear_obj(polySo)
2294 284 storres
        sollya_lib_clear_obj(funcSo)
2295 284 storres
        print "The Sollya object created from the function is not a function."
2296 284 storres
        return None
2297 284 storres
    #pobyso_autoprint(funcSo)
2298 284 storres
    # Convert the Sage interval into a Sollya range.
2299 284 storres
    rangeSo = pobyso_interval_to_range_sa_so(intervalSa)
2300 284 storres
    if not pobyso_is_range_so_sa(rangeSo):
2301 284 storres
        sollya_lib_clear_obj(polySo)
2302 284 storres
        sollya_lib_clear_obj(funcSo)
2303 284 storres
        sollya_lib_clear_obj(rangeSo)
2304 284 storres
        print "The Sollya object created from the interval is not a range."
2305 284 storres
        return None
2306 284 storres
    #pobyso_autoprint(rangeSo)
2307 284 storres
    # Check the error type. We do not check the returned object since we
2308 284 storres
    # trust our code and Sollya on this one.
2309 284 storres
    if errorTypeSa != SOLLYA_ABSOLUTE and errorTypeSa != SOLLYA_RELATIVE:
2310 284 storres
        sollya_lib_clear_obj(polySo)
2311 284 storres
        sollya_lib_clear_obj(funcSo)
2312 284 storres
        sollya_lib_clear_obj(rangeSo)
2313 284 storres
        return None
2314 284 storres
    if errorTypeSa == SOLLYA_ABSOLUTE:
2315 284 storres
        errorTypeSo = pobyso_absolute_so_so()
2316 284 storres
    else:
2317 284 storres
        errorTypeSo = pobyso_relative_so_so()
2318 284 storres
    #pobyso_autoprint(errorTypeSo)
2319 284 storres
    # Check if accuracySa is an element of a RealField. If not, try to
2320 284 storres
    # convert it.
2321 284 storres
    try:
2322 284 storres
        accuracySa.ulp()
2323 284 storres
    except AttributeError:
2324 284 storres
        try:
2325 284 storres
            accuracySa = RR(accuracySa)
2326 284 storres
        except TypeError:
2327 284 storres
            accuracySa = RR(str(accuracySa))
2328 284 storres
    # Create the Sollya object
2329 284 storres
    accuracySo = pobyso_constant_sa_so(accuracySa)
2330 284 storres
    #pobyso_autoprint(accuracySo)
2331 284 storres
    if pobyso_is_error_so_sa(accuracySo):
2332 284 storres
        sollya_lib_clear_obj(polySo)
2333 284 storres
        sollya_lib_clear_obj(funcSo)
2334 284 storres
        sollya_lib_clear_obj(rangeSo)
2335 284 storres
        sollya_lib_clear_obj(errorTypeSo)
2336 284 storres
        sollya_lib_clear_obj(accuracySo)
2337 284 storres
        return None
2338 284 storres
    retValSo = sollya_lib_supnorm(polySo,
2339 284 storres
                                  funcSo,
2340 284 storres
                                  rangeSo,
2341 284 storres
                                  errorTypeSo,
2342 284 storres
                                  accuracySo,
2343 284 storres
                                  None)
2344 284 storres
    sollya_lib_clear_obj(polySo)
2345 284 storres
    sollya_lib_clear_obj(funcSo)
2346 284 storres
    sollya_lib_clear_obj(rangeSo)
2347 284 storres
    sollya_lib_clear_obj(errorTypeSo)
2348 284 storres
    sollya_lib_clear_obj(accuracySo)
2349 284 storres
    if pobyso_is_error_so_sa(retValSo):
2350 284 storres
        sollya_lib_clear_obj(retValSo)
2351 284 storres
        return None
2352 284 storres
    #pobyso_autoprint(retValSo)
2353 284 storres
    retValSa = pobyso_range_to_interval_so_sa(retValSo)
2354 284 storres
    sollya_lib_clear_obj(retValSo)
2355 284 storres
    return retValSa
2356 281 storres
# End pobyso_supnorm_sa_sa
2357 281 storres
2358 242 storres
def pobyso_supnorm_so_sa(polySo, funcSo, intervalSo, errorTypeSo = None,\
2359 242 storres
                         accuracySo = None, realFieldSa = None):
2360 242 storres
    """
2361 281 storres
    Computes the supremum norm from Sollya input arguments and returns a
2362 242 storres
    Sage floating-point number whose precision is set by the only Sage argument.
2363 215 storres

2364 242 storres
    The returned value is the maximum of the absolute values of the range
2365 242 storres
    elements  returned by the Sollya supnorm functions
2366 242 storres
    """
2367 242 storres
    supNormRangeSo = pobyso_supnorm_so_so(polySo,
2368 242 storres
                                          funcSo,
2369 242 storres
                                          intervalSo,
2370 242 storres
                                          errorTypeSo,
2371 242 storres
                                          accuracySo)
2372 242 storres
    supNormSo = pobyso_range_max_abs_so_so(supNormRangeSo)
2373 242 storres
    sollya_lib_clear_obj(supNormRangeSo)
2374 242 storres
    #pobyso_autoprint(supNormSo)
2375 242 storres
    supNormSa = pobyso_get_constant_as_rn_with_rf_so_sa(supNormSo, realFieldSa)
2376 242 storres
    sollya_lib_clear_obj(supNormSo)
2377 242 storres
    return supNormSa
2378 242 storres
# End pobyso_supnorm_so_sa.
2379 242 storres
#
2380 85 storres
def pobyso_supnorm_so_so(polySo, funcSo, intervalSo, errorTypeSo = None,\
2381 85 storres
                         accuracySo = None):
2382 58 storres
    """
2383 85 storres
    Computes the supnorm of the approximation error between the given
2384 242 storres
    polynomial and function. Attention: returns a range!
2385 85 storres
    errorTypeSo defaults to "absolute".
2386 85 storres
    accuracySo defaults to 2^(-40).
2387 85 storres
    """
2388 85 storres
    if errorTypeSo is None:
2389 85 storres
        errorTypeSo = sollya_lib_absolute(None)
2390 85 storres
        errorTypeIsNone = True
2391 85 storres
    else:
2392 85 storres
        errorTypeIsNone = False
2393 85 storres
    #
2394 85 storres
    if accuracySo is None:
2395 237 storres
        # Notice the **: we are in Pythonland here!
2396 85 storres
        accuracySo = pobyso_constant_sa_so(RR(2**(-40)))
2397 85 storres
        accuracyIsNone = True
2398 85 storres
    else:
2399 85 storres
        accuracyIsNone = False
2400 238 storres
    #pobyso_autoprint(accuracySo)
2401 85 storres
    resultSo = \
2402 85 storres
        sollya_lib_supnorm(polySo, funcSo, intervalSo, errorTypeSo, \
2403 85 storres
                              accuracySo)
2404 85 storres
    if errorTypeIsNone:
2405 85 storres
        sollya_lib_clear_obj(errorTypeSo)
2406 85 storres
    if accuracyIsNone:
2407 85 storres
        sollya_lib_clear_obj(accuracySo)
2408 85 storres
    return resultSo
2409 85 storres
# End pobyso_supnorm_so_so
2410 242 storres
#
2411 162 storres
def pobyso_taylor_expansion_no_change_var_so_so(functionSo,
2412 162 storres
                                                degreeSo,
2413 162 storres
                                                rangeSo,
2414 162 storres
                                                errorTypeSo=None,
2415 162 storres
                                                sollyaPrecSo=None):
2416 85 storres
    """
2417 162 storres
    Compute the Taylor expansion without the variable change
2418 162 storres
    x -> x-intervalCenter.
2419 272 storres
    If errorTypeSo is None, absolute is used.
2420 272 storres
    If sollyaPrecSo is None, Sollya internal precision is not changed.
2421 58 storres
    """
2422 226 storres
    # Change internal Sollya precision, if needed.
2423 227 storres
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
2424 226 storres
    sollyaPrecChanged = False
2425 226 storres
    if sollyaPrecSo is None:
2426 226 storres
        pass
2427 226 storres
    else:
2428 58 storres
        sollya_lib_set_prec(sollyaPrecSo)
2429 226 storres
        sollyaPrecChanged = True
2430 85 storres
    # Error type stuff: default to absolute.
2431 85 storres
    if errorTypeSo is None:
2432 85 storres
        errorTypeIsNone = True
2433 85 storres
        errorTypeSo = sollya_lib_absolute(None)
2434 85 storres
    else:
2435 85 storres
        errorTypeIsNone = False
2436 162 storres
    intervalCenterSo = sollya_lib_mid(rangeSo, None)
2437 162 storres
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo,
2438 162 storres
                                         intervalCenterSo,
2439 58 storres
                                         rangeSo, errorTypeSo, None)
2440 272 storres
    # Object taylorFormListSaSo is a Python list of Sollya objects references
2441 272 storres
    #  that are copies of the elements of taylorFormSo.
2442 117 storres
    # pobyso_get_list_elements_so_so clears taylorFormSo.
2443 162 storres
    (taylorFormListSaSo, numElementsSa, isEndEllipticSa) = \
2444 58 storres
        pobyso_get_list_elements_so_so(taylorFormSo)
2445 272 storres
    ## Copy needed here since polySo will be returned and taylorFormListSaSo
2446 272 storres
    #  will be cleared.
2447 162 storres
    polySo = sollya_lib_copy_obj(taylorFormListSaSo[0])
2448 162 storres
    #print "Num elements:", numElementsSa
2449 162 storres
    sollya_lib_clear_obj(taylorFormSo)
2450 181 storres
    # No copy_obj needed here: a new objects are created.
2451 272 storres
    maxErrorSo    = sollya_lib_sup(taylorFormListSaSo[2])
2452 272 storres
    minErrorSo    = sollya_lib_inf(taylorFormListSaSo[2])
2453 272 storres
    # List taylorFormListSaSo is not needed anymore.
2454 280 storres
    pobyso_clear_list_elements_sa_so(taylorFormListSaSo)
2455 181 storres
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
2456 181 storres
    absMinErrorSo = sollya_lib_abs(minErrorSo)
2457 181 storres
    sollya_lib_clear_obj(maxErrorSo)
2458 181 storres
    sollya_lib_clear_obj(minErrorSo)
2459 181 storres
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
2460 181 storres
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
2461 272 storres
    #
2462 272 storres
    if errorTypeIsNone:
2463 272 storres
        sollya_lib_clear_obj(errorTypeSo)
2464 272 storres
    ## If changed, reset the Sollya working precision.
2465 226 storres
    if sollyaPrecChanged:
2466 226 storres
        sollya_lib_set_prec(initialSollyaPrecSo)
2467 226 storres
    sollya_lib_clear_obj(initialSollyaPrecSo)
2468 272 storres
    ## According to what error is the largest, return the errors.
2469 181 storres
    if absMaxErrorSa > absMinErrorSa:
2470 181 storres
        sollya_lib_clear_obj(absMinErrorSo)
2471 227 storres
        return (polySo, intervalCenterSo, absMaxErrorSo)
2472 181 storres
    else:
2473 181 storres
        sollya_lib_clear_obj(absMaxErrorSo)
2474 227 storres
        return (polySo, intervalCenterSo, absMinErrorSo)
2475 162 storres
# end pobyso_taylor_expansion_no_change_var_so_so
2476 58 storres
2477 162 storres
def pobyso_taylor_expansion_with_change_var_so_so(functionSo, degreeSo, \
2478 162 storres
                                                  rangeSo, \
2479 162 storres
                                                  errorTypeSo=None, \
2480 162 storres
                                                  sollyaPrecSo=None):
2481 58 storres
    """
2482 162 storres
    Compute the Taylor expansion with the variable change
2483 162 storres
    x -> (x-intervalCenter) included.
2484 58 storres
    """
2485 226 storres
    # Change Sollya internal precision, if need.
2486 226 storres
    sollyaPrecChanged = False
2487 271 storres
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
2488 226 storres
    if sollyaPrecSo is None:
2489 226 storres
        pass
2490 226 storres
    else:
2491 56 storres
        sollya_lib_set_prec(sollyaPrecSo)
2492 226 storres
        sollyaPrecChanged = True
2493 162 storres
    #
2494 85 storres
    # Error type stuff: default to absolute.
2495 85 storres
    if errorTypeSo is None:
2496 85 storres
        errorTypeIsNone = True
2497 85 storres
        errorTypeSo = sollya_lib_absolute(None)
2498 85 storres
    else:
2499 85 storres
        errorTypeIsNone = False
2500 162 storres
    intervalCenterSo = sollya_lib_mid(rangeSo)
2501 162 storres
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo, \
2502 162 storres
                                         intervalCenterSo, \
2503 56 storres
                                         rangeSo, errorTypeSo, None)
2504 272 storres
    # Object taylorFormListSaSo is a Python list of Sollya objects references
2505 272 storres
    # that are copies of the elements of taylorFormSo.
2506 116 storres
    # pobyso_get_list_elements_so_so clears taylorFormSo.
2507 272 storres
    (taylorFormListSaSo, numElements, isEndElliptic) = \
2508 56 storres
        pobyso_get_list_elements_so_so(taylorFormSo)
2509 272 storres
    sollya_lib_clear_obj(taylorFormSo)
2510 272 storres
    polySo = taylorFormListSaSo[0]
2511 272 storres
    ## Maximum error computation with taylorFormListSaSo[2], a range
2512 272 storres
    #  holding the actual error. Bounds can be negative.
2513 272 storres
    maxErrorSo    = sollya_lib_sup(taylorFormListSaSo[2])
2514 272 storres
    minErrorSo    = sollya_lib_inf(taylorFormListSaSo[2])
2515 181 storres
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
2516 181 storres
    absMinErrorSo = sollya_lib_abs(minErrorSo)
2517 181 storres
    sollya_lib_clear_obj(maxErrorSo)
2518 181 storres
    sollya_lib_clear_obj(minErrorSo)
2519 181 storres
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
2520 181 storres
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
2521 162 storres
    changeVarExpSo = sollya_lib_build_function_sub(\
2522 162 storres
                       sollya_lib_build_function_free_variable(),\
2523 162 storres
                       sollya_lib_copy_obj(intervalCenterSo))
2524 181 storres
    polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpSo)
2525 272 storres
    # List taylorFormListSaSo is not needed anymore.
2526 280 storres
    pobyso_clear_list_elements_sa_so(taylorFormListSaSo)
2527 162 storres
    sollya_lib_clear_obj(changeVarExpSo)
2528 56 storres
    # If changed, reset the Sollya working precision.
2529 226 storres
    if sollyaPrecChanged:
2530 226 storres
        sollya_lib_set_prec(initialSollyaPrecSo)
2531 226 storres
    sollya_lib_clear_obj(initialSollyaPrecSo)
2532 85 storres
    if errorTypeIsNone:
2533 85 storres
        sollya_lib_clear_obj(errorTypeSo)
2534 162 storres
    # Do not clear maxErrorSo.
2535 181 storres
    if absMaxErrorSa > absMinErrorSa:
2536 181 storres
        sollya_lib_clear_obj(absMinErrorSo)
2537 181 storres
        return((polyVarChangedSo, intervalCenterSo, absMaxErrorSo))
2538 181 storres
    else:
2539 181 storres
        sollya_lib_clear_obj(absMaxErrorSo)
2540 181 storres
        return((polyVarChangedSo, intervalCenterSo, absMinErrorSo))
2541 162 storres
# end pobyso_taylor_expansion_with_change_var_so_so
2542 56 storres
2543 5 storres
def pobyso_taylor(function, degree, point):
2544 38 storres
    """ Legacy function. See pobysoTaylor_so_so. """
2545 38 storres
    return(pobyso_taylor_so_so(function, degree, point))
2546 38 storres
2547 56 storres
def pobyso_taylor_so_so(functionSo, degreeSo, pointSo):
2548 56 storres
    return(sollya_lib_taylor(functionSo, degreeSo, pointSo))
2549 5 storres
2550 85 storres
def pobyso_taylorform(function, degree, point = None,
2551 85 storres
                      interval = None, errorType=None):
2552 85 storres
    """ Legacy function. See pobyso_taylorform_sa_sa;"""
2553 38 storres
2554 38 storres
def pobyso_taylorform_sa_sa(functionSa, \
2555 84 storres
                            degreeSa, \
2556 84 storres
                            pointSa, \
2557 84 storres
                            intervalSa=None, \
2558 84 storres
                            errorTypeSa=None, \
2559 84 storres
                            precisionSa=None):
2560 37 storres
    """
2561 85 storres
    Compute the Taylor form of 'degreeSa' for 'functionSa' at 'pointSa'
2562 85 storres
    for 'intervalSa' with 'errorTypeSa' (a string) using 'precisionSa'.
2563 37 storres
    point: must be a Real or a Real interval.
2564 37 storres
    return the Taylor form as an array
2565 83 storres
    TODO: take care of the interval and of the point when it is an interval;
2566 38 storres
          when errorType is not None;
2567 83 storres
          take care of the other elements of the Taylor form (coefficients
2568 83 storres
          errors and delta.
2569 37 storres
    """
2570 37 storres
    # Absolute as the default error.
2571 84 storres
    if errorTypeSa is None:
2572 37 storres
        errorTypeSo = sollya_lib_absolute()
2573 84 storres
    elif errorTypeSa == "relative":
2574 84 storres
        errorTypeSo = sollya_lib_relative()
2575 84 storres
    elif errortypeSa == "absolute":
2576 84 storres
        errorTypeSo = sollya_lib_absolute()
2577 37 storres
    else:
2578 84 storres
        # No clean up needed.
2579 84 storres
        return None
2580 84 storres
    # Global precision stuff
2581 226 storres
    sollyaPrecisionChangedSa = False
2582 226 storres
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
2583 226 storres
    if precisionSa is None:
2584 226 storres
        precSa = initialSollyaPrecSa
2585 226 storres
    else:
2586 226 storres
        if precSa > initialSollyaPrecSa:
2587 226 storres
            if precSa <= 2:
2588 226 storres
                print inspect.stack()[0][3], ":precision change <= 2 requested."
2589 226 storres
            pobyso_set_prec_sa_so(precSa)
2590 226 storres
            sollyaPrecisionChangedSa = True
2591 226 storres
    #
2592 85 storres
    if len(functionSa.variables()) > 0:
2593 85 storres
        varSa = functionSa.variables()[0]
2594 85 storres
        pobyso_name_free_variable_sa_so(str(varSa))
2595 84 storres
    # In any case (point or interval) the parent of pointSa has a precision
2596 84 storres
    # method.
2597 84 storres
    pointPrecSa = pointSa.parent().precision()
2598 226 storres
    if precSa > pointPrecSa:
2599 226 storres
        pointPrecSa = precSa
2600 84 storres
    # In any case (point or interval) pointSa has a base_ring() method.
2601 84 storres
    pointBaseRingString = str(pointSa.base_ring())
2602 84 storres
    if re.search('Interval', pointBaseRingString) is None: # Point
2603 84 storres
        pointSo = pobyso_constant_sa_so(pointSa, pointPrecSa)
2604 84 storres
    else: # Interval.
2605 84 storres
        pointSo = pobyso_interval_to_range_sa_so(pointSa, pointPrecSa)
2606 37 storres
    # Sollyafy the function.
2607 159 storres
    functionSo = pobyso_parse_string_sa_so(functionSa._assume_str().replace('_SAGE_VAR_', ''))
2608 37 storres
    if sollya_lib_obj_is_error(functionSo):
2609 37 storres
        print "pobyso_tailorform: function string can't be parsed!"
2610 37 storres
        return None
2611 37 storres
    # Sollyafy the degree
2612 84 storres
    degreeSo = sollya_lib_constant_from_int(int(degreeSa))
2613 37 storres
    # Sollyafy the point
2614 37 storres
    # Call Sollya
2615 83 storres
    taylorFormSo = \
2616 83 storres
        sollya_lib_taylorform(functionSo, degreeSo, pointSo, errorTypeSo,\
2617 37 storres
                                         None)
2618 85 storres
    sollya_lib_clear_obj(functionSo)
2619 85 storres
    sollya_lib_clear_obj(degreeSo)
2620 85 storres
    sollya_lib_clear_obj(pointSo)
2621 85 storres
    sollya_lib_clear_obj(errorTypeSo)
2622 38 storres
    (tfsAsList, numElements, isEndElliptic) = \
2623 38 storres
            pobyso_get_list_elements_so_so(taylorFormSo)
2624 37 storres
    polySo = tfsAsList[0]
2625 38 storres
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
2626 37 storres
    polyRealField = RealField(maxPrecision)
2627 38 storres
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, polyRealField)
2628 226 storres
    if sollyaPrecisionChangedSa:
2629 226 storres
        sollya_lib_set_prec(initialSollyaPrecSo)
2630 226 storres
    sollya_lib_clear_obj(initialSollyaPrecSo)
2631 37 storres
    polynomialRing = polyRealField[str(varSa)]
2632 37 storres
    polySa = polynomial(expSa, polynomialRing)
2633 37 storres
    taylorFormSa = [polySa]
2634 85 storres
    # Final clean-up.
2635 85 storres
    sollya_lib_clear_obj(taylorFormSo)
2636 51 storres
    return(taylorFormSa)
2637 51 storres
# End pobyso_taylor_form_sa_sa
2638 54 storres
2639 54 storres
def pobyso_taylorform_so_so(functionSo, degreeSo, pointSo, intervalSo=None, \
2640 54 storres
                            errorTypeSo=None):
2641 54 storres
    createdErrorType = False
2642 51 storres
    if errorTypeSo is None:
2643 51 storres
        errorTypeSo = sollya_lib_absolute()
2644 54 storres
        createdErrorType = True
2645 51 storres
    else:
2646 51 storres
        #TODO: deal with the other case.
2647 51 storres
        pass
2648 51 storres
    if intervalSo is None:
2649 54 storres
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
2650 54 storres
                                         errorTypeSo, None)
2651 51 storres
    else:
2652 54 storres
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
2653 54 storres
                                         intervalSo, errorTypeSo, None)
2654 54 storres
    if createdErrorType:
2655 54 storres
        sollya_lib_clear_obj(errorTypeSo)
2656 215 storres
    return resultSo
2657 51 storres
2658 37 storres
2659 37 storres
def pobyso_univar_polynomial_print_reverse(polySa):
2660 51 storres
    """ Legacy function. See pobyso_univar_polynomial_print_reverse_sa_sa. """
2661 51 storres
    return(pobyso_univar_polynomial_print_reverse_sa_sa(polySa))
2662 38 storres
2663 51 storres
def pobyso_univar_polynomial_print_reverse_sa_sa(polySa):
2664 37 storres
    """
2665 37 storres
    Return the string representation of a univariate polynomial with
2666 38 storres
    monomials ordered in the x^0..x^n order of the monomials.
2667 37 storres
    Remember: Sage
2668 37 storres
    """
2669 37 storres
    polynomialRing = polySa.base_ring()
2670 37 storres
    # A very expensive solution:
2671 37 storres
    # -create a fake multivariate polynomial field with only one variable,
2672 37 storres
    #   specifying a negative lexicographical order;
2673 37 storres
    mpolynomialRing = PolynomialRing(polynomialRing.base(), \
2674 37 storres
                                     polynomialRing.variable_name(), \
2675 37 storres
                                     1, order='neglex')
2676 37 storres
    # - convert the univariate argument polynomial into a multivariate
2677 37 storres
    #   version;
2678 37 storres
    p = mpolynomialRing(polySa)
2679 37 storres
    # - return the string representation of the converted form.
2680 37 storres
    # There is no simple str() method defined for p's class.
2681 37 storres
    return(p.__str__())
2682 5 storres
#
2683 230 storres
#print pobyso_get_prec()
2684 5 storres
pobyso_set_prec(165)
2685 230 storres
#print pobyso_get_prec()
2686 230 storres
#a=100
2687 230 storres
#print type(a)
2688 230 storres
#id(a)
2689 230 storres
#print "Max arity: ", pobyso_max_arity
2690 230 storres
#print "Function tripleDouble (43) as a string: ", pobyso_function_type_as_string(43)
2691 230 storres
#print "Function None (44) as a string: ", pobyso_function_type_as_string(44)
2692 246 storres
sys.stderr.write("\t...Pobyso check done.\n")