Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 303

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

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