Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 283

Historique | Voir | Annoter | Télécharger (93,91 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 38 storres
Create the equivalent to an enum for the Sollya function types.
40 38 storres
"""
41 5 storres
(SOLLYA_BASE_FUNC_ABS,
42 5 storres
SOLLYA_BASE_FUNC_ACOS,
43 5 storres
    SOLLYA_BASE_FUNC_ACOSH,
44 5 storres
    SOLLYA_BASE_FUNC_ADD,
45 5 storres
    SOLLYA_BASE_FUNC_ASIN,
46 5 storres
    SOLLYA_BASE_FUNC_ASINH,
47 5 storres
    SOLLYA_BASE_FUNC_ATAN,
48 5 storres
    SOLLYA_BASE_FUNC_ATANH,
49 5 storres
    SOLLYA_BASE_FUNC_CEIL,
50 5 storres
    SOLLYA_BASE_FUNC_CONSTANT,
51 5 storres
    SOLLYA_BASE_FUNC_COS,
52 5 storres
    SOLLYA_BASE_FUNC_COSH,
53 5 storres
    SOLLYA_BASE_FUNC_DIV,
54 5 storres
    SOLLYA_BASE_FUNC_DOUBLE,
55 5 storres
    SOLLYA_BASE_FUNC_DOUBLEDOUBLE,
56 5 storres
    SOLLYA_BASE_FUNC_DOUBLEEXTENDED,
57 5 storres
    SOLLYA_BASE_FUNC_ERF,
58 5 storres
    SOLLYA_BASE_FUNC_ERFC,
59 5 storres
    SOLLYA_BASE_FUNC_EXP,
60 5 storres
    SOLLYA_BASE_FUNC_EXP_M1,
61 5 storres
    SOLLYA_BASE_FUNC_FLOOR,
62 5 storres
    SOLLYA_BASE_FUNC_FREE_VARIABLE,
63 5 storres
    SOLLYA_BASE_FUNC_HALFPRECISION,
64 5 storres
    SOLLYA_BASE_FUNC_LIBRARYCONSTANT,
65 5 storres
    SOLLYA_BASE_FUNC_LIBRARYFUNCTION,
66 5 storres
    SOLLYA_BASE_FUNC_LOG,
67 5 storres
    SOLLYA_BASE_FUNC_LOG_10,
68 5 storres
    SOLLYA_BASE_FUNC_LOG_1P,
69 5 storres
    SOLLYA_BASE_FUNC_LOG_2,
70 5 storres
    SOLLYA_BASE_FUNC_MUL,
71 5 storres
    SOLLYA_BASE_FUNC_NEARESTINT,
72 5 storres
    SOLLYA_BASE_FUNC_NEG,
73 5 storres
    SOLLYA_BASE_FUNC_PI,
74 5 storres
    SOLLYA_BASE_FUNC_POW,
75 5 storres
    SOLLYA_BASE_FUNC_PROCEDUREFUNCTION,
76 5 storres
    SOLLYA_BASE_FUNC_QUAD,
77 5 storres
    SOLLYA_BASE_FUNC_SIN,
78 5 storres
    SOLLYA_BASE_FUNC_SINGLE,
79 5 storres
    SOLLYA_BASE_FUNC_SINH,
80 5 storres
    SOLLYA_BASE_FUNC_SQRT,
81 5 storres
    SOLLYA_BASE_FUNC_SUB,
82 5 storres
    SOLLYA_BASE_FUNC_TAN,
83 5 storres
    SOLLYA_BASE_FUNC_TANH,
84 5 storres
SOLLYA_BASE_FUNC_TRIPLEDOUBLE) = map(int,xrange(44))
85 246 storres
sys.stderr.write("Superficial pobyso check...\n")
86 230 storres
#print "First constant - SOLLYA_BASE_FUNC_ABS: ", SOLLYA_BASE_FUNC_ABS
87 230 storres
#print "Last constant  - SOLLYA_BASE_FUNC_TRIPLEDOUBLE: ", SOLLYA_BASE_FUNC_TRIPLEDOUBLE
88 5 storres
89 5 storres
pobyso_max_arity = 9
90 5 storres
91 85 storres
def pobyso_absolute_so_so():
92 259 storres
    """
93 259 storres
    Create an "absolute" Sollya object.
94 259 storres
    """
95 85 storres
    return(sollya_lib_absolute(None))
96 85 storres
97 5 storres
def pobyso_autoprint(arg):
98 218 storres
    sollya_lib_autoprint(arg, None)
99 5 storres
100 38 storres
def pobyso_autoprint_so_so(arg):
101 38 storres
    sollya_lib_autoprint(arg,None)
102 282 storres
103 282 storres
def pobyso_bounds_to_interval_sa_sa(lowerBound, upperBound):
104 282 storres
    """
105 282 storres
    Convert a pair of bounds into an interval (an element of
106 282 storres
    a RealIntervalField).
107 282 storres
    """
108 282 storres
    # Minimal (not bullet-proof) check on bounds.
109 282 storres
    if lowerBound > upperBound:
110 282 storres
        return None
111 282 storres
    # Try to get the maximum precision among the bounds.
112 282 storres
    try:
113 282 storres
        preclb = parent(lowerBound).precision()
114 282 storres
        precub = parent(upperBound).precision()
115 282 storres
        prec = max(preclb, precub)
116 282 storres
    except AttributeError:
117 282 storres
        prec = 53
118 282 storres
    # Create the RealIntervalField and the interval (if possible).
119 282 storres
    theRIF = RealIntervalField(prec)
120 282 storres
    try:
121 282 storres
        interval = theRIF(lowerBound, upperBound)
122 282 storres
    except TypeError:
123 282 storres
        return None
124 282 storres
    else:
125 282 storres
        return interval
126 282 storres
# End pobyso_bounds_to_interval_sa_sa
127 54 storres
128 84 storres
def pobyso_bounds_to_range_sa_so(rnLowerBoundSa, rnUpperBoundSa, \
129 84 storres
                                 precisionSa=None):
130 84 storres
    """
131 84 storres
    Return a Sollya range from to 2 RealField Sage elements.
132 84 storres
    The Sollya range element has a sufficient precision to hold all
133 162 storres
    the digits of the widest of the Sage bounds.
134 84 storres
    """
135 84 storres
    # Sanity check.
136 84 storres
    if rnLowerBoundSa > rnUpperBoundSa:
137 84 storres
        return None
138 115 storres
    # Precision stuff.
139 85 storres
    if precisionSa is None:
140 84 storres
        # Check for the largest precision.
141 84 storres
        lbPrecSa = rnLowerBoundSa.parent().precision()
142 84 storres
        ubPrecSa = rnLowerBoundSa.parent().precision()
143 84 storres
        maxPrecSa = max(lbPrecSa, ubPrecSa)
144 84 storres
    else:
145 84 storres
        maxPrecSa = precisionSa
146 115 storres
    # From Sage to Sollya bounds.
147 162 storres
#    lowerBoundSo = sollya_lib_constant(get_rn_value(rnLowerBoundSa),
148 162 storres
#                                       maxPrecSa)
149 162 storres
    lowerBoundSo = pobyso_constant_sa_so(rnLowerBoundSa,
150 162 storres
                                         maxPrecSa)
151 162 storres
    upperBoundSo = pobyso_constant_sa_so(rnUpperBoundSa,
152 154 storres
                                       maxPrecSa)
153 162 storres
154 115 storres
    # From Sollya bounds to range.
155 84 storres
    rangeSo = sollya_lib_range(lowerBoundSo, upperBoundSo)
156 84 storres
    # Back to original precision.
157 84 storres
    # Clean up
158 84 storres
    sollya_lib_clear_obj(lowerBoundSo)
159 84 storres
    sollya_lib_clear_obj(upperBoundSo)
160 154 storres
    return rangeSo
161 84 storres
# End pobyso_bounds_to_range_sa_so
162 84 storres
163 215 storres
def pobyso_build_end_elliptic_list_so_so(*args):
164 215 storres
    """
165 237 storres
    From Sollya argument objects, create a Sollya end elliptic list.
166 237 storres
    Elements of the list are "eaten" (should not be cleared individually,
167 215 storres
    are cleared when the list is cleared).
168 215 storres
    """
169 215 storres
    if len(args) == 0:
170 280 storres
        ## When called with an empty list, sollya_lib_build_end_elliptic_list,
171 280 storres
        #  produces "error".
172 280 storres
        return sollya_lib_build_list(None)
173 215 storres
    index = 0
174 215 storres
    ## One can not append elements to an elliptic list, prepend only is
175 215 storres
    #  permitted.
176 215 storres
    for argument in reversed(args):
177 215 storres
        if index == 0:
178 215 storres
            listSo = sollya_lib_build_end_elliptic_list(argument, None)
179 215 storres
        else:
180 215 storres
            listSo = sollya_lib_prepend(argument, listSo)
181 215 storres
        index += 1
182 215 storres
    return listSo
183 215 storres
184 215 storres
# End pobyso_build_end_elliptic_list_so_so
185 215 storres
186 54 storres
def pobyso_build_function_sub_so_so(exp1So, exp2So):
187 237 storres
    return sollya_lib_build_function_sub(exp1So, exp2So)
188 54 storres
189 237 storres
def pobyso_build_list_of_ints_sa_so(*args):
190 237 storres
    """
191 237 storres
    Build a Sollya list from Sage integral arguments.
192 237 storres
    """
193 237 storres
    if len(args) == 0:
194 237 storres
        return pobyso_build_list_so_so()
195 237 storres
    ## Make a Sage list of integral constants.
196 237 storres
    intsList = []
197 237 storres
    for intElem in args:
198 237 storres
        intsList.append(pobyso_constant_from_int_sa_so(intElem))
199 237 storres
    return pobyso_build_list_so_so(*intsList)
200 237 storres
201 237 storres
def pobyso_build_list_so_so(*args):
202 237 storres
    """
203 237 storres
    Make a Sollya list out of Sollya objects passed as arguments.
204 237 storres
    If one wants to call it with a list argument, should prepend a "*"
205 237 storres
    before the list variable name.
206 237 storres
    Elements of the list are "eaten" (should not be cleared individually,
207 237 storres
    are cleared when the list is cleared).
208 237 storres
    """
209 237 storres
    if len(args) == 0:
210 237 storres
        ## Called with an empty list produced "error".
211 237 storres
        return sollya_lib_build_list(None)
212 237 storres
    index = 0
213 237 storres
    ## Append the Sollya elements one by one.
214 237 storres
    for elementSo in args:
215 237 storres
        if index == 0:
216 237 storres
            listSo = sollya_lib_build_list(elementSo, None)
217 237 storres
        else:
218 237 storres
            listSo = sollya_lib_append(listSo, elementSo)
219 237 storres
        index += 1
220 237 storres
    return listSo
221 237 storres
# End pobyso_build list_so_so
222 237 storres
223 237 storres
224 85 storres
def pobyso_change_var_in_function_so_so(funcSo, chvarExpSo):
225 54 storres
    """
226 85 storres
    Variable change in a function.
227 85 storres
    """
228 85 storres
    return(sollya_lib_evaluate(funcSo,chvarExpSo))
229 85 storres
# End pobyso_change_var_in_function_so_so
230 85 storres
231 85 storres
def pobyso_chebyshevform_so_so(functionSo, degreeSo, intervalSo):
232 85 storres
    resultSo = sollya_lib_chebyshevform(functionSo, degreeSo, intervalSo)
233 85 storres
    return(resultSo)
234 85 storres
# End pobyso_chebyshevform_so_so.
235 85 storres
236 280 storres
def pobyso_clear_full_list_elements_sa_so(objectListSaSo):
237 280 storres
    """
238 280 storres
    Clear the elements of list created by the
239 280 storres
    pobyso_get_list_elements_so_so function.
240 280 storres
    objectListSaSo is as follows:
241 280 storres
    - objectListSaSo[0]: a list of Sollya objects;
242 280 storres
    - objectListSaSo[1]: the number of elements of the previous list;
243 280 storres
    - objectListSaSo[2]: an integer that if != 0 states that the list is
244 280 storres
                         end-elliptic
245 280 storres
    The objects to clear are the elements of the objectListSaSo[0] list.
246 280 storres
    """
247 280 storres
    for index in xrange(0, objectListSaSo[1]):
248 280 storres
        sollya_lib_clear_obj(objectListSaSo[0][index])
249 280 storres
# End pobyso_clear_full_list_elements_sa_so
250 280 storres
251 280 storres
def pobyso_clear_list_elements_sa_so(objectListSaSo):
252 280 storres
    """
253 280 storres
    Clear the elements of list of references to Sollya objects
254 280 storres
    """
255 280 storres
    for index in xrange(0, len(objectListSaSo)):
256 280 storres
        sollya_lib_clear_obj(objectListSaSo[index])
257 280 storres
# End pobyso_clear_list_elements_sa_so
258 280 storres
259 237 storres
def pobyso_clear_obj(objSo):
260 237 storres
    """
261 237 storres
    Free a Sollya object's memory.
262 237 storres
    Very thin wrapper around sollya_lib_clear_obj().
263 237 storres
    """
264 237 storres
    sollya_lib_clear_obj(objSo)
265 237 storres
# End pobyso_clear_obj.
266 237 storres
267 117 storres
def pobyso_clear_taylorform_sa_so(taylorFormSaSo):
268 117 storres
    """
269 280 storres
    This method is rapper around pobyso_clear_list_elements_sa_so.
270 280 storres
    It is a legacy method left here since it may be used in existing code
271 280 storres
    where Taylor forms are used as Sage lists obtained by converting
272 280 storres
    Sollya Taylor forms (a list made of:
273 280 storres
    - a polynomial;
274 280 storres
    - a list of intervals enclosing the errors accumulated when computing
275 280 storres
      the polynomial coefficients;
276 280 storres
    - a bound on the approximation error between the polynomial and the
277 280 storres
      function.)
278 280 storres
    A Taylor form directly obtained from pobyso_taylorform_so_so is cleared
279 280 storres
    by sollya_lib_clear_obj.
280 117 storres
    """
281 280 storres
    pobyso_clear_list_elements_sa_so(taylorFormSaSo)
282 117 storres
# End pobyso_clear_taylorform_sa_so
283 117 storres
284 85 storres
def pobyso_cmp(rnArgSa, cteSo):
285 85 storres
    """
286 237 storres
    Deprecated, use pobyso_cmp_sa_so_sa instead.
287 237 storres
    """
288 237 storres
    print "Deprecated, use pobyso_cmp_sa_so_sa instead."
289 237 storres
    return pobyso_cmp_sa_so_sa(rnArgSa, cteSo)
290 237 storres
# End  pobyso_cmp
291 237 storres
292 237 storres
def pobyso_cmp_sa_so_sa(rnArgSa, cteSo):
293 237 storres
    """
294 54 storres
    Compare the MPFR value a RealNumber with that of a Sollya constant.
295 54 storres

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

887 118 storres
    INPUT:
888 118 storres
    - objectListSo: a Sollya list of Sollya objects.
889 118 storres

890 118 storres
    OUTPUT:
891 118 storres
    - a Sage/Python tuple made of:
892 118 storres
      - a Sage/Python list of Sollya objects,
893 118 storres
      - a Sage/Python int holding the number of elements,
894 118 storres
      - a Sage/Python int stating (!= 0) that the list is end-elliptic.
895 118 storres
    NOTE::
896 118 storres
        We recover the addresses of the Sollya object from the list of pointers
897 118 storres
        returned by sollya_lib_get_list_elements. The list itself is freed.
898 118 storres
    TODO::
899 118 storres
        Figure out what to do with numElements since the number of elements
900 118 storres
        can easily be recovered from the list itself.
901 118 storres
        Ditto for isEndElliptic.
902 51 storres
    """
903 5 storres
    listAddress = POINTER(c_longlong)()
904 5 storres
    numElements = c_int(0)
905 5 storres
    isEndElliptic = c_int(0)
906 117 storres
    listAsSageList = []
907 5 storres
    result = sollya_lib_get_list_elements(byref(listAddress),\
908 54 storres
                                          byref(numElements),\
909 54 storres
                                          byref(isEndElliptic),\
910 117 storres
                                          objectListSo)
911 5 storres
    if result == 0 :
912 5 storres
        return None
913 5 storres
    for i in xrange(0, numElements.value, 1):
914 118 storres
       #listAsSageList.append(sollya_lib_copy_obj(listAddress[i]))
915 118 storres
       listAsSageList.append(listAddress[i])
916 117 storres
       # Clear each of the elements returned by Sollya.
917 118 storres
       #sollya_lib_clear_obj(listAddress[i])
918 117 storres
    # Free the list itself.
919 117 storres
    sollya_lib_free(listAddress)
920 209 storres
    return (listAsSageList, numElements.value, isEndElliptic.value)
921 5 storres
922 38 storres
def pobyso_get_max_prec_of_exp(soExp):
923 38 storres
    """ Legacy function. See pobyso_get_max_prec_of_exp_so_sa. """
924 209 storres
    return pobyso_get_max_prec_of_exp_so_sa(soExp)
925 5 storres
926 85 storres
def pobyso_get_max_prec_of_exp_so_sa(expSo):
927 38 storres
    """
928 38 storres
    Get the maximum precision used for the numbers in a Sollya expression.
929 52 storres

930 52 storres
    Arguments:
931 52 storres
    soExp -- a Sollya expression pointer
932 52 storres
    Return value:
933 52 storres
    A Python integer
934 38 storres
    TODO:
935 38 storres
    - error management;
936 38 storres
    - correctly deal with numerical type such as DOUBLEEXTENDED.
937 38 storres
    """
938 226 storres
    if expSo is None:
939 226 storres
        print inspect.stack()[0][3], ": expSo is None."
940 226 storres
        return 0
941 5 storres
    maxPrecision = 0
942 52 storres
    minConstPrec = 0
943 52 storres
    currentConstPrec = 0
944 226 storres
    #pobyso_autoprint(expSo)
945 85 storres
    operator = pobyso_get_head_function_so_sa(expSo)
946 5 storres
    if (operator != SOLLYA_BASE_FUNC_CONSTANT) and \
947 5 storres
    (operator != SOLLYA_BASE_FUNC_FREE_VARIABLE):
948 85 storres
        (arity, subexpressions) = pobyso_get_subfunctions_so_sa(expSo)
949 5 storres
        for i in xrange(arity):
950 5 storres
            maxPrecisionCandidate = \
951 38 storres
                pobyso_get_max_prec_of_exp_so_sa(subexpressions[i])
952 5 storres
            if maxPrecisionCandidate > maxPrecision:
953 5 storres
                maxPrecision = maxPrecisionCandidate
954 209 storres
        return maxPrecision
955 5 storres
    elif operator == SOLLYA_BASE_FUNC_CONSTANT:
956 85 storres
        #minConstPrec = pobyso_get_min_prec_of_constant_so_sa(expSo)
957 52 storres
        #currentConstPrec = pobyso_get_min_prec_of_constant_so_sa(soExp)
958 52 storres
        #print minConstPrec, " - ", currentConstPrec
959 209 storres
        return pobyso_get_min_prec_of_constant_so_sa(expSo)
960 52 storres
961 5 storres
    elif operator == SOLLYA_BASE_FUNC_FREE_VARIABLE:
962 209 storres
        return 0
963 5 storres
    else:
964 38 storres
        print "pobyso_get_max_prec_of_exp_so_sa: unexepected operator."
965 209 storres
        return 0
966 5 storres
967 85 storres
def pobyso_get_min_prec_of_constant_so_sa(constExpSo):
968 52 storres
    """
969 52 storres
    Get the minimum precision necessary to represent the value of a Sollya
970 52 storres
    constant.
971 52 storres
    MPFR_MIN_PREC and powers of 2 are taken into account.
972 209 storres
    We assume that constExpSo is a pointer to a Sollay constant expression.
973 52 storres
    """
974 85 storres
    constExpAsRnSa = pobyso_get_constant_as_rn_so_sa(constExpSo)
975 85 storres
    return(min_mpfr_size(get_rn_value(constExpAsRnSa)))
976 52 storres
977 200 storres
def pobyso_get_poly_so_sa(polySo, realFieldSa=None):
978 200 storres
    """
979 200 storres
    Convert a Sollya polynomial into a Sage polynomial.
980 209 storres
    Legacy function. Use pobyso_float_poly_so_sa() instead.
981 200 storres
    """
982 213 storres
    return pobyso_float_poly_so_sa(polySo,realFieldSa)
983 200 storres
# End pobyso_get_poly_so_sa
984 200 storres
985 200 storres
def pobyso_get_prec():
986 200 storres
    """ Legacy function. See pobyso_get_prec_so_sa(). """
987 209 storres
    return pobyso_get_prec_so_sa()
988 200 storres
989 200 storres
def pobyso_get_prec_so():
990 200 storres
    """
991 200 storres
    Get the current default precision in Sollya.
992 200 storres
    The return value is a Sollya object.
993 200 storres
    Usefull when modifying the precision back and forth by avoiding
994 200 storres
    extra conversions.
995 200 storres
    """
996 209 storres
    return sollya_lib_get_prec(None)
997 200 storres
998 200 storres
def pobyso_get_prec_so_sa():
999 200 storres
    """
1000 200 storres
    Get the current default precision in Sollya.
1001 200 storres
    The return value is Sage/Python int.
1002 200 storres
    """
1003 227 storres
    precSo = sollya_lib_get_prec()
1004 227 storres
    precSa = pobyso_constant_from_int_so_sa(precSo)
1005 200 storres
    sollya_lib_clear_obj(precSo)
1006 227 storres
    return precSa
1007 200 storres
# End pobyso_get_prec_so_sa.
1008 200 storres
1009 209 storres
def pobyso_get_prec_so_so_sa():
1010 209 storres
    """
1011 209 storres
    Return the current precision both as a Sollya object and a
1012 209 storres
    Sage integer as hybrid tuple.
1013 209 storres
    To avoid multiple calls for precision manipulations.
1014 209 storres
    """
1015 227 storres
    precSo = sollya_lib_get_prec()
1016 227 storres
    precSa = pobyso_constant_from_int_so_sa(precSo)
1017 227 storres
    return (precSo, int(precSa))
1018 200 storres
1019 200 storres
def pobyso_get_prec_of_constant(ctExpSo):
1020 200 storres
    """ Legacy function. See pobyso_get_prec_of_constant_so_sa. """
1021 209 storres
    return pobyso_get_prec_of_constant_so_sa(ctExpSo)
1022 200 storres
1023 200 storres
def pobyso_get_prec_of_constant_so_sa(ctExpSo):
1024 200 storres
    """
1025 200 storres
    Tries to find a precision to represent ctExpSo without rounding.
1026 200 storres
    If not possible, returns None.
1027 200 storres
    """
1028 218 storres
    #print "Entering pobyso_get_prec_of_constant_so_sa..."
1029 200 storres
    prec = c_int(0)
1030 200 storres
    retc = sollya_lib_get_prec_of_constant(byref(prec), ctExpSo, None)
1031 200 storres
    if retc == 0:
1032 218 storres
        #print "pobyso_get_prec_of_constant_so_sa failed."
1033 209 storres
        return None
1034 218 storres
    #print "...exiting pobyso_get_prec_of_constant_so_sa."
1035 209 storres
    return int(prec.value)
1036 200 storres
1037 200 storres
def pobyso_get_prec_of_range_so_sa(rangeSo):
1038 200 storres
    """
1039 200 storres
    Returns the number of bits elements of a range are coded with.
1040 200 storres
    """
1041 200 storres
    prec = c_int(0)
1042 200 storres
    retc = sollya_lib_get_prec_of_range(byref(prec), rangeSo, None)
1043 200 storres
    if retc == 0:
1044 200 storres
        return(None)
1045 209 storres
    return int(prec.value)
1046 200 storres
# End pobyso_get_prec_of_range_so_sa()
1047 200 storres
1048 85 storres
def pobyso_get_sage_exp_from_sollya_exp(sollyaExpSo, realField = RR):
1049 38 storres
    """ Legacy function. See pobyso_get_sage_exp_from_sollya_exp_so_sa. """
1050 209 storres
    return pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo,
1051 209 storres
                                                     realField = RR)
1052 38 storres
1053 85 storres
def pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, realFieldSa = RR):
1054 5 storres
    """
1055 38 storres
    Get a Sage expression from a Sollya expression.
1056 38 storres
    Currently only tested with polynomials with floating-point coefficients.
1057 5 storres
    Notice that, in the returned polynomial, the exponents are RealNumbers.
1058 5 storres
    """
1059 5 storres
    #pobyso_autoprint(sollyaExp)
1060 85 storres
    operatorSa = pobyso_get_head_function_so_sa(sollyaExpSo)
1061 83 storres
    sollyaLibFreeVariableName = sollya_lib_get_free_variable_name()
1062 213 storres
    ## Get rid of the "_"'s in "_x_", if any.
1063 213 storres
    sollyaLibFreeVariableName = re.sub('_', '', sollyaLibFreeVariableName)
1064 5 storres
    # Constants and the free variable are special cases.
1065 5 storres
    # All other operator are dealt with in the same way.
1066 85 storres
    if (operatorSa != SOLLYA_BASE_FUNC_CONSTANT) and \
1067 85 storres
       (operatorSa != SOLLYA_BASE_FUNC_FREE_VARIABLE):
1068 85 storres
        (aritySa, subexpressionsSa) = pobyso_get_subfunctions_so_sa(sollyaExpSo)
1069 85 storres
        if aritySa == 1:
1070 85 storres
            sageExpSa = eval(pobyso_function_type_as_string_so_sa(operatorSa) + \
1071 85 storres
            "(" + pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], \
1072 85 storres
            realFieldSa) + ")")
1073 85 storres
        elif aritySa == 2:
1074 63 storres
            # We do not get through the preprocessor.
1075 63 storres
            # The "^" operator is then a special case.
1076 85 storres
            if operatorSa == SOLLYA_BASE_FUNC_POW:
1077 85 storres
                operatorAsStringSa = "**"
1078 5 storres
            else:
1079 85 storres
                operatorAsStringSa = \
1080 85 storres
                    pobyso_function_type_as_string_so_sa(operatorSa)
1081 85 storres
            sageExpSa = \
1082 85 storres
              eval("pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], realFieldSa)"\
1083 85 storres
              + " " + operatorAsStringSa + " " + \
1084 85 storres
                   "pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[1], realFieldSa)")
1085 63 storres
        # We do not know yet how to deal with arity >= 3
1086 63 storres
        # (is there any in Sollya anyway?).
1087 5 storres
        else:
1088 85 storres
            sageExpSa = eval('None')
1089 209 storres
        return sageExpSa
1090 85 storres
    elif operatorSa == SOLLYA_BASE_FUNC_CONSTANT:
1091 5 storres
        #print "This is a constant"
1092 85 storres
        return pobyso_get_constant_as_rn_with_rf_so_sa(sollyaExpSo, realFieldSa)
1093 85 storres
    elif operatorSa == SOLLYA_BASE_FUNC_FREE_VARIABLE:
1094 218 storres
        #print "This is the free variable"
1095 209 storres
        return eval(sollyaLibFreeVariableName)
1096 5 storres
    else:
1097 5 storres
        print "Unexpected"
1098 5 storres
        return eval('None')
1099 185 storres
# End pobyso_get_sage_exp_from_sollya_exp_so_sa
1100 73 storres
1101 185 storres
1102 38 storres
def pobyso_get_subfunctions(expressionSo):
1103 38 storres
    """ Legacy function. See pobyso_get_subfunctions_so_sa. """
1104 209 storres
    return pobyso_get_subfunctions_so_sa(expressionSo)
1105 200 storres
# End pobyso_get_subfunctions.
1106 200 storres
1107 38 storres
def pobyso_get_subfunctions_so_sa(expressionSo):
1108 38 storres
    """
1109 38 storres
    Get the subfunctions of an expression.
1110 38 storres
    Return the number of subfunctions and the list of subfunctions addresses.
1111 55 storres
    S.T.: Could not figure out another way than that ugly list of declarations
1112 83 storres
    to recover the addresses of the subfunctions.
1113 83 storres
    We limit ourselves to arity 8 functions.
1114 38 storres
    """
1115 5 storres
    subf0 = c_int(0)
1116 5 storres
    subf1 = c_int(0)
1117 5 storres
    subf2 = c_int(0)
1118 5 storres
    subf3 = c_int(0)
1119 5 storres
    subf4 = c_int(0)
1120 5 storres
    subf5 = c_int(0)
1121 5 storres
    subf6 = c_int(0)
1122 5 storres
    subf7 = c_int(0)
1123 5 storres
    subf8 = c_int(0)
1124 5 storres
    arity = c_int(0)
1125 5 storres
    nullPtr = POINTER(c_int)()
1126 38 storres
    sollya_lib_get_subfunctions(expressionSo, byref(arity), \
1127 83 storres
      byref(subf0), byref(subf1), byref(subf2), byref(subf3), \
1128 83 storres
      byref(subf4), byref(subf5),\
1129 83 storres
      byref(subf6), byref(subf7), byref(subf8), nullPtr, None)
1130 83 storres
#    byref(cast(subfunctions[0], POINTER(c_int))), \
1131 83 storres
#    byref(cast(subfunctions[0], POINTER(c_int))), \
1132 83 storres
#    byref(cast(subfunctions[2], POINTER(c_int))), \
1133 83 storres
#    byref(cast(subfunctions[3], POINTER(c_int))), \
1134 83 storres
#    byref(cast(subfunctions[4], POINTER(c_int))), \
1135 83 storres
#    byref(cast(subfunctions[5], POINTER(c_int))), \
1136 83 storres
#    byref(cast(subfunctions[6], POINTER(c_int))), \
1137 83 storres
#    byref(cast(subfunctions[7], POINTER(c_int))), \
1138 5 storres
#    byref(cast(subfunctions[8], POINTER(c_int))), nullPtr)
1139 83 storres
    subfunctions = [subf0, subf1, subf2, subf3, subf4, subf5, subf6, subf7, \
1140 83 storres
                    subf8]
1141 5 storres
    subs = []
1142 5 storres
    if arity.value > pobyso_max_arity:
1143 38 storres
        return(0,[])
1144 5 storres
    for i in xrange(arity.value):
1145 5 storres
        subs.append(int(subfunctions[i].value))
1146 5 storres
        #print subs[i]
1147 209 storres
    return (int(arity.value), subs)
1148 200 storres
# End pobyso_get_subfunctions_so_sa
1149 5 storres
1150 155 storres
def pobyso_guess_degree_sa_sa(functionSa, intervalSa, approxErrorSa,
1151 155 storres
                              weightSa=None, degreeBoundSa=None):
1152 155 storres
    """
1153 155 storres
    Sa_sa variant of the solly_guessdegree function.
1154 155 storres
    Return 0 if something goes wrong.
1155 155 storres
    """
1156 159 storres
    functionAsStringSa = functionSa._assume_str().replace('_SAGE_VAR_', '')
1157 154 storres
    functionSo = pobyso_parse_string_sa_so(functionAsStringSa)
1158 155 storres
    if pobyso_is_error_so_sa(functionSo):
1159 155 storres
        sollya_lib_clear_obj(functionSo)
1160 155 storres
        return 0
1161 154 storres
    rangeSo = pobyso_interval_to_range_sa_so(intervalSa)
1162 155 storres
    # The approximation error is expected to be a floating point number.
1163 155 storres
    if pobyso_is_floating_point_number_sa_sa(approxErrorSa):
1164 155 storres
        approxErrorSo = pobyso_constant_sa_so(approxErrorSa)
1165 155 storres
    else:
1166 155 storres
        approxErrorSo = pobyso_constant_sa_so(RR(approxErrorSa))
1167 154 storres
    if not weightSa is None:
1168 159 storres
        weightAsStringSa = weightSa._assume_str().replace('_SAGE_VAR_', '')
1169 154 storres
        weightSo = pobyso_parse_string_sa_so(weightAsStringSa)
1170 166 storres
        if pobyso_is_error_so_sa(weightSo):
1171 155 storres
            sollya_lib_clear_obj(functionSo)
1172 155 storres
            sollya_lib_clear_obj(rangeSo)
1173 155 storres
            sollya_lib_clear_obj(approxErrorSo)
1174 155 storres
            sollya_lib_clear_obj(weightSo)
1175 155 storres
            return 0
1176 154 storres
    else:
1177 154 storres
        weightSo = None
1178 154 storres
    if not degreeBoundSa is None:
1179 154 storres
        degreeBoundSo = pobyso_constant_from_int_sa_so(degreeBoundSa)
1180 154 storres
    else:
1181 154 storres
        degreeBoundSo = None
1182 154 storres
    guessedDegreeSa = pobyso_guess_degree_so_sa(functionSo,
1183 162 storres
                                                rangeSo,
1184 162 storres
                                                approxErrorSo,
1185 162 storres
                                                weightSo,
1186 162 storres
                                                degreeBoundSo)
1187 154 storres
    sollya_lib_clear_obj(functionSo)
1188 154 storres
    sollya_lib_clear_obj(rangeSo)
1189 155 storres
    sollya_lib_clear_obj(approxErrorSo)
1190 154 storres
    if not weightSo is None:
1191 154 storres
        sollya_lib_clear_obj(weightSo)
1192 154 storres
    if not degreeBoundSo is None:
1193 154 storres
        sollya_lib_clear_obj(degreeBoundSo)
1194 154 storres
    return guessedDegreeSa
1195 154 storres
# End poyso_guess_degree_sa_sa
1196 154 storres
1197 153 storres
def pobyso_guess_degree_so_sa(functionSo, rangeSo, errorSo, weightSo=None, \
1198 154 storres
                              degreeBoundSo=None):
1199 154 storres
    """
1200 154 storres
    Thin wrapper around the guessdegree function.
1201 154 storres
    Nevertheless, some precision control stuff has been appended.
1202 154 storres
    """
1203 154 storres
    # Deal with Sollya internal precision issues: if it is too small,
1204 154 storres
    # compared with the error, increases it to about twice -log2(error).
1205 154 storres
    errorSa = pobyso_get_constant_as_rn_with_rf_so_sa(errorSo)
1206 154 storres
    log2ErrorSa = errorSa.log2()
1207 154 storres
    if log2ErrorSa < 0:
1208 154 storres
        neededPrecisionSa = int(2 * int(-log2ErrorSa) / 64) * 64
1209 154 storres
    else:
1210 154 storres
        neededPrecisionSa = int(2 * int(log2ErrorSa) / 64) * 64
1211 154 storres
    #print "Needed precision:", neededPrecisionSa
1212 226 storres
    sollyaPrecisionHasChanged = False
1213 226 storres
    (initialPrecSo, initialPrecSa) = pobyso_get_prec_so_so_sa()
1214 226 storres
    if neededPrecisionSa > initialPrecSa:
1215 226 storres
        if neededPrecisionSa <= 2:
1216 226 storres
            print inspect.stack()[0][3], ": precision change <= 2 requested."
1217 154 storres
        pobyso_set_prec_sa_so(neededPrecisionSa)
1218 226 storres
        sollyaPrecisionHasChanged = True
1219 166 storres
    #print "Guessing degree..."
1220 153 storres
    # weightSo and degreeBoundsSo are optional arguments.
1221 162 storres
    # As declared, sollya_lib_guessdegree must take 5 arguments.
1222 153 storres
    if weightSo is None:
1223 162 storres
        degreeRangeSo = sollya_lib_guessdegree(functionSo, rangeSo, errorSo,
1224 162 storres
                                               0, 0, None)
1225 154 storres
    elif degreeBoundSo is None:
1226 153 storres
        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, \
1227 162 storres
                                                errorSo, weightSo, 0, None)
1228 153 storres
    else:
1229 153 storres
        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, errorSo, \
1230 154 storres
                                                weightSo, degreeBoundSo, None)
1231 166 storres
    #print "...degree guess done."
1232 154 storres
    # Restore internal precision, if applicable.
1233 226 storres
    if sollyaPrecisionHasChanged:
1234 226 storres
        pobyso_set_prec_so_so(initialPrecSo)
1235 226 storres
        sollya_lib_clear_obj(initialPrecSo)
1236 154 storres
    degreeIntervalSa = pobyso_range_to_interval_so_sa(degreeRangeSo)
1237 154 storres
    sollya_lib_clear_obj(degreeRangeSo)
1238 154 storres
    # When ok, both bounds match.
1239 154 storres
    # When the degree bound is too low, the upper bound is the degree
1240 154 storres
    # for which the error can be honored.
1241 154 storres
    # When it really goes wrong, the upper bound is infinity.
1242 154 storres
    if degreeIntervalSa.lower() == degreeIntervalSa.upper():
1243 154 storres
        return int(degreeIntervalSa.lower())
1244 154 storres
    else:
1245 154 storres
        if degreeIntervalSa.upper().is_infinity():
1246 154 storres
            return None
1247 154 storres
        else:
1248 154 storres
            return int(degreeIntervalSa.upper())
1249 154 storres
    # End pobyso_guess_degree_so_sa
1250 153 storres
1251 215 storres
def pobyso_inf_so_so(intervalSo):
1252 215 storres
    """
1253 215 storres
    Very thin wrapper around sollya_lib_inf().
1254 215 storres
    """
1255 215 storres
    return sollya_lib_inf(intervalSo)
1256 215 storres
# End pobyso_inf_so_so.
1257 282 storres
#
1258 282 storres
def pobyso_infnorm_sa_sa(funcSa, intervalSa):
1259 282 storres
    """
1260 282 storres
    An infnorm call with Sage arguments.
1261 282 storres
    We only take into account the 2 first arguments (the function and
1262 282 storres
    the interval (a range). Managing the other arguments (the file for
1263 282 storres
    the proof and the exclusion intervals list) will be performed later
1264 282 storres
    Changes will be needed in sollya_lib.py file too.
1265 282 storres
    """
1266 282 storres
    # Check that funcSa is a function.
1267 282 storres
    if not \
1268 282 storres
     sage.symbolic.callable.is_CallableSymbolicExpressionRing(parent(funcSa)):
1269 282 storres
        return None
1270 282 storres
    # Check that intervalSa is an interval.
1271 282 storres
    try:
1272 282 storres
        intervalSa.upper()
1273 282 storres
    except AttributeError:
1274 282 storres
        return None
1275 282 storres
    # Convert the Sage function into a Sollya function.
1276 282 storres
    funcAsStringSa = funcSa._assume_str().replace('_SAGE_VAR_', '')
1277 282 storres
    funcSo = pobyso_parse_string_sa_so(funcAsStringSa)
1278 282 storres
    if not pobyso_obj_is_function_so_sa(funcSo):
1279 282 storres
        sollya_lib_clear_obj(funcSo)
1280 282 storres
        return None
1281 282 storres
    # Convert the Sage interval into a Sollya range.
1282 282 storres
    rangeSo = pobyso_interval_to_range_sa_so(intervalSa)
1283 282 storres
    retValSo = sollya_lib_infnorm(funcSo, rangeSo, None)
1284 282 storres
    sollya_lib_clear_obj(funcSo)
1285 282 storres
    sollya_lib_clear_obj(rangeSo)
1286 282 storres
    if pobyso_is_error_so_sa(retValSo):
1287 282 storres
        sollya_lib_clear_obj(retValSo)
1288 282 storres
        return None
1289 282 storres
    retValSa = pobyso_range_to_interval_so_sa(retValSo)
1290 282 storres
    sollya_lib_clear_obj(retValSo)
1291 282 storres
    return retValSa
1292 282 storres
# End pobyso_infnorm_so_so.
1293 282 storres
#
1294 282 storres
def pobyso_infnorm_so_so(funcSo, rangeSo):
1295 282 storres
    """
1296 282 storres
    Very thin wrapper around sollya_lib_infnorm().
1297 282 storres
    We only take into account the 2 first arguments (the function and
1298 282 storres
    the interval (a range). Managing the other arguments (the file for
1299 282 storres
    the proof and the exclusion intervals list) will be performed later
1300 282 storres
    Changes will be needed in sollya_lib.py file too.
1301 215 storres

1302 282 storres
    As per Sollya manual, this function should not be used anymore and
1303 282 storres
    supnorm should be called instead. Nevertheless, supnorm breaks
1304 282 storres
    sometimes whereas infnorm still returns a satisfactory answer.
1305 282 storres
    """
1306 282 storres
    return sollya_lib_infnorm(funcSo, rangeSo, None)
1307 282 storres
# End pobyso_infnorm_so_so.
1308 53 storres
1309 84 storres
def pobyso_interval_to_range_sa_so(intervalSa, precisionSa=None):
1310 84 storres
    if precisionSa is None:
1311 84 storres
        precisionSa = intervalSa.parent().precision()
1312 84 storres
    intervalSo = pobyso_bounds_to_range_sa_so(intervalSa.lower(),\
1313 84 storres
                                              intervalSa.upper(),\
1314 84 storres
                                              precisionSa)
1315 209 storres
    return intervalSo
1316 84 storres
# End pobyso_interval_to_range_sa_so
1317 84 storres
1318 155 storres
def pobyso_is_error_so_sa(objSo):
1319 155 storres
    """
1320 155 storres
    Thin wrapper around the sollya_lib_obj_is_error() function.
1321 155 storres
    """
1322 155 storres
    if sollya_lib_obj_is_error(objSo) != 0:
1323 155 storres
        return True
1324 155 storres
    else:
1325 155 storres
        return False
1326 155 storres
# End pobyso_is_error-so_sa
1327 155 storres
1328 155 storres
def pobyso_is_floating_point_number_sa_sa(numberSa):
1329 155 storres
    """
1330 209 storres
    Check whether a Sage number is floating point.
1331 209 storres
    Exception stuff added because numbers other than
1332 209 storres
    floating-point ones do not have the is_real() attribute.
1333 155 storres
    """
1334 209 storres
    try:
1335 209 storres
        return numberSa.is_real()
1336 209 storres
    except AttributeError:
1337 209 storres
        return False
1338 209 storres
# End pobyso_is_floating_piont_number_sa_sa
1339 155 storres
1340 226 storres
def pobyso_is_range_so_sa(rangeCandidateSo):
1341 226 storres
    """
1342 226 storres
    Thin wrapper over sollya_lib_is_range.
1343 226 storres
    """
1344 226 storres
    return sollya_lib_obj_is_range(rangeCandidateSo) != 0
1345 226 storres
1346 226 storres
# End pobyso_is_range_so_sa
1347 226 storres
1348 226 storres
1349 37 storres
def pobyso_lib_init():
1350 37 storres
    sollya_lib_init(None)
1351 116 storres
1352 116 storres
def pobyso_lib_close():
1353 116 storres
    sollya_lib_close(None)
1354 237 storres
1355 85 storres
def pobyso_name_free_variable(freeVariableNameSa):
1356 38 storres
    """ Legacy function. See pobyso_name_free_variable_sa_so. """
1357 85 storres
    pobyso_name_free_variable_sa_so(freeVariableNameSa)
1358 38 storres
1359 85 storres
def pobyso_name_free_variable_sa_so(freeVariableNameSa):
1360 83 storres
    """
1361 83 storres
    Set the free variable name in Sollya from a Sage string.
1362 83 storres
    """
1363 85 storres
    sollya_lib_name_free_variable(freeVariableNameSa)
1364 37 storres
1365 281 storres
def pobyso_obj_is_function_so_sa(objSo):
1366 281 storres
    """
1367 281 storres
    Check if an object is a function.
1368 281 storres
    """
1369 281 storres
    if sollya_lib_obj_is_function(objSo) != 0:
1370 281 storres
        return True
1371 281 storres
    else:
1372 281 storres
        return False
1373 281 storres
# End pobyso_obj_is_function_so_sa
1374 281 storres
1375 281 storres
def pobyso_obj_is_range_so_sa(objSo):
1376 281 storres
    """
1377 281 storres
    Check if an object is a function.
1378 281 storres
    """
1379 281 storres
    if sollya_lib_obj_is_range(objSo) != 0:
1380 281 storres
        return True
1381 281 storres
    else:
1382 281 storres
        return False
1383 281 storres
# End pobyso_obj_is_range_so_sa
1384 281 storres
1385 281 storres
def pobyso_obj_is_string_so_sa(objSo):
1386 281 storres
    """
1387 281 storres
    Check if an object is a function.
1388 281 storres
    """
1389 281 storres
    if sollya_lib_obj_is_string(objSo) != 0:
1390 281 storres
        return True
1391 281 storres
    else:
1392 281 storres
        return False
1393 281 storres
# End pobyso_obj_is_string_so_sa
1394 281 storres
1395 5 storres
def pobyso_parse_string(string):
1396 38 storres
    """ Legacy function. See pobyso_parse_string_sa_so. """
1397 209 storres
    return pobyso_parse_string_sa_so(string)
1398 38 storres
1399 38 storres
def pobyso_parse_string_sa_so(string):
1400 83 storres
    """
1401 155 storres
    Get the Sollya expression computed from a Sage string or
1402 155 storres
    a Sollya error object if parsing failed.
1403 83 storres
    """
1404 209 storres
    return sollya_lib_parse_string(string)
1405 5 storres
1406 200 storres
def pobyso_precision_so_sa(ctExpSo):
1407 209 storres
    """
1408 209 storres
    Computes the necessary precision to represent a number.
1409 209 storres
    If x is not zero, it can be uniquely written as x = m · 2e
1410 209 storres
    where m is an odd integer and e is an integer.
1411 209 storres
    precision(x) returns the number of bits necessary to write m
1412 209 storres
    in binary (i.e. ceil(log2(m))).
1413 209 storres
    """
1414 209 storres
    #TODO: take care of the special case: 0, @NaN@, @Inf@
1415 200 storres
    precisionSo = sollya_lib_precision(ctExpSo)
1416 200 storres
    precisionSa = pobyso_constant_from_int_so_sa(precisionSo)
1417 200 storres
    sollya_lib_clear_obj(precisionSo)
1418 200 storres
    return precisionSa
1419 200 storres
# End pobyso_precision_so_sa
1420 215 storres
1421 217 storres
def pobyso_polynomial_coefficients_progressive_round_so_so(polySo,
1422 217 storres
                                                           funcSo,
1423 217 storres
                                                           icSo,
1424 217 storres
                                                           intervalSo,
1425 217 storres
                                                           itpSo,
1426 217 storres
                                                           ftpSo,
1427 217 storres
                                                           maxPrecSo,
1428 219 storres
                                                           maxErrSo,
1429 219 storres
                                                           debug=False):
1430 219 storres
    if debug:
1431 219 storres
        print "Input arguments:"
1432 219 storres
        pobyso_autoprint(polySo)
1433 219 storres
        pobyso_autoprint(funcSo)
1434 219 storres
        pobyso_autoprint(icSo)
1435 219 storres
        pobyso_autoprint(intervalSo)
1436 219 storres
        pobyso_autoprint(itpSo)
1437 219 storres
        pobyso_autoprint(ftpSo)
1438 219 storres
        pobyso_autoprint(maxPrecSo)
1439 219 storres
        pobyso_autoprint(maxErrSo)
1440 219 storres
        print "________________"
1441 200 storres
1442 217 storres
    ## Higher order function see:
1443 217 storres
    #  http://effbot.org/pyfaq/how-do-you-make-a-higher-order-function-in-python.htm
1444 217 storres
    def precision_decay_ratio_function(degreeSa):
1445 217 storres
        def outer(x):
1446 217 storres
            def inner(x):
1447 217 storres
                we = 3/8
1448 217 storres
                wq = 2/8
1449 217 storres
                a  = 2.2
1450 217 storres
                b  = 2
1451 217 storres
                return we*(exp(x/a)-1) +  wq*((b*x)**2) + (1-we-wq)*x
1452 217 storres
            return  inner(x)/inner(degreeSa)
1453 217 storres
        return outer
1454 217 storres
1455 217 storres
    #
1456 217 storres
    degreeSa        = pobyso_polynomial_degree_so_sa(polySo)
1457 217 storres
    ratio           = precision_decay_ratio_function(degreeSa)
1458 217 storres
    itpSa           = pobyso_constant_from_int_so_sa(itpSo)
1459 217 storres
    ftpSa           = pobyso_constant_from_int_so_sa(ftpSo)
1460 217 storres
    maxPrecSa       = pobyso_constant_from_int_so_sa(maxPrecSo)
1461 217 storres
    maxErrSa        = pobyso_get_constant_as_rn_so_sa(maxErrSo)
1462 219 storres
    if debug:
1463 219 storres
        print "degreeSa:", degreeSa
1464 219 storres
        print "ratio:", ratio
1465 219 storres
        print "itpsSa:", itpSa
1466 219 storres
        print "ftpSa:", ftpSa
1467 219 storres
        print "maxPrecSa:", maxPrecSa
1468 219 storres
        print "maxErrSa:", maxErrSa
1469 217 storres
    lastResPolySo   = None
1470 218 storres
    lastInfNormSo   = None
1471 219 storres
    #print "About to enter the while loop..."
1472 217 storres
    while True:
1473 218 storres
        resPolySo   = pobyso_constant_0_sa_so()
1474 217 storres
        pDeltaSa    = ftpSa - itpSa
1475 217 storres
        for indexSa in reversed(xrange(0,degreeSa+1)):
1476 218 storres
            #print "Index:", indexSa
1477 217 storres
            indexSo = pobyso_constant_from_int_sa_so(indexSa)
1478 217 storres
            #pobyso_autoprint(indexSo)
1479 217 storres
            #print ratio(indexSa)
1480 217 storres
            ctpSa = floor(ftpSa - (pDeltaSa * ratio(indexSa)))
1481 217 storres
            ctpSo = pobyso_constant_from_int_sa_so(ctpSa)
1482 219 storres
            if debug:
1483 219 storres
                print "Index:", indexSa, " - Target precision:",
1484 219 storres
                pobyso_autoprint(ctpSo)
1485 217 storres
            cmonSo  = \
1486 217 storres
                sollya_lib_build_function_mul(sollya_lib_coeff(polySo, indexSo),
1487 217 storres
                                      sollya_lib_build_function_pow( \
1488 217 storres
                                          sollya_lib_build_function_free_variable(), \
1489 217 storres
                                          indexSo))
1490 217 storres
            #pobyso_autoprint(cmonSo)
1491 217 storres
            cmonrSo = pobyso_round_coefficients_single_so_so(cmonSo, ctpSo)
1492 217 storres
            sollya_lib_clear_obj(cmonSo)
1493 217 storres
            #pobyso_autoprint(cmonrSo)
1494 217 storres
            resPolySo = sollya_lib_build_function_add(resPolySo,
1495 217 storres
                                                      cmonrSo)
1496 218 storres
            #pobyso_autoprint(resPolySo)
1497 217 storres
        # End for index
1498 217 storres
        freeVarSo     = sollya_lib_build_function_free_variable()
1499 217 storres
        changeVarSo   = sollya_lib_sub(freeVarSo, icSo)
1500 217 storres
        resPolyCvSo   = sollya_lib_evaluate(resPolySo, changeVarSo)
1501 218 storres
        errFuncSo = sollya_lib_build_function_sub(sollya_lib_copy_obj(funcSo),
1502 218 storres
                                                  resPolyCvSo)
1503 218 storres
        infNormSo = sollya_lib_dirtyinfnorm(errFuncSo, intervalSo)
1504 217 storres
        cerrSa    = pobyso_get_constant_as_rn_so_sa(infNormSo)
1505 219 storres
        if debug:
1506 219 storres
            print "Infnorm (Sollya):",
1507 219 storres
            pobyso_autoprint(infNormSo)
1508 218 storres
        sollya_lib_clear_obj(errFuncSo)
1509 218 storres
        #print "Infnorm  (Sage):", cerrSa
1510 217 storres
        if (cerrSa > maxErrSa):
1511 219 storres
            if debug:
1512 219 storres
                print "Error is too large."
1513 218 storres
            if lastResPolySo is None:
1514 219 storres
                if debug:
1515 219 storres
                    print "Enlarging prec."
1516 217 storres
                ntpSa = floor(ftpSa + ftpSa/50)
1517 217 storres
                ## Can't enlarge (numerical)
1518 217 storres
                if ntpSa == ftpSa:
1519 217 storres
                    sollya_lib_clear_obj(resPolySo)
1520 217 storres
                    return None
1521 217 storres
                ## Can't enlarge (not enough precision left)
1522 217 storres
                if ntpSa > maxPrecSa:
1523 217 storres
                    sollya_lib_clear_obj(resPolySo)
1524 217 storres
                    return None
1525 217 storres
                ftpSa = ntpSa
1526 217 storres
                continue
1527 217 storres
            ## One enlargement took place.
1528 217 storres
            else:
1529 219 storres
                if debug:
1530 219 storres
                    print "Exit with the last before last polynomial."
1531 222 storres
                    print "Precision of highest degree monomial:", itpSa
1532 222 storres
                    print "Precision of constant term          :", ftpSa
1533 217 storres
                sollya_lib_clear_obj(resPolySo)
1534 218 storres
                sollya_lib_clear_obj(infNormSo)
1535 218 storres
                return (lastResPolySo, lastInfNormSo)
1536 218 storres
        # cerrSa <= maxErrSa: scrap more bits, possibly.
1537 217 storres
        else:
1538 219 storres
            if debug:
1539 219 storres
                print "Error is too small"
1540 219 storres
            if cerrSa <= (maxErrSa/2):
1541 219 storres
                if debug:
1542 219 storres
                    print "Shrinking prec."
1543 218 storres
                ntpSa = floor(ftpSa - ftpSa/50)
1544 218 storres
                ## Can't shrink (numerical)
1545 218 storres
                if ntpSa == ftpSa:
1546 218 storres
                    if not lastResPolySo is None:
1547 218 storres
                        sollya_lib_clear_obj(lastResPolySo)
1548 218 storres
                    if not lastInfNormSo is None:
1549 218 storres
                        sollya_lib_clear_obj(lastInfNormSo)
1550 222 storres
                    if debug:
1551 222 storres
                        print "Exit because can't shrink anymore (numerically)."
1552 222 storres
                        print "Precision of highest degree monomial:", itpSa
1553 222 storres
                        print "Precision of constant term          :", ftpSa
1554 218 storres
                    return (resPolySo, infNormSo)
1555 218 storres
                ## Can't shrink (not enough precision left)
1556 218 storres
                if ntpSa <= itpSa:
1557 218 storres
                    if not lastResPolySo is None:
1558 218 storres
                        sollya_lib_clear_obj(lastResPolySo)
1559 218 storres
                    if not lastInfNormSo is None:
1560 218 storres
                        sollya_lib_clear_obj(lastInfNormSo)
1561 222 storres
                        print "Exit because can't shrink anymore (no bits left)."
1562 222 storres
                        print "Precision of highest degree monomial:", itpSa
1563 222 storres
                        print "Precision of constant term          :", ftpSa
1564 218 storres
                    return (resPolySo, infNormSo)
1565 218 storres
                ftpSa = ntpSa
1566 217 storres
                if not lastResPolySo is None:
1567 217 storres
                    sollya_lib_clear_obj(lastResPolySo)
1568 218 storres
                if not lastInfNormSo is None:
1569 218 storres
                    sollya_lib_clear_obj(lastInfNormSo)
1570 218 storres
                lastResPolySo = resPolySo
1571 218 storres
                lastInfNormSo = infNormSo
1572 218 storres
                continue
1573 218 storres
            else: # Error is not that small, just return
1574 217 storres
                if not lastResPolySo is None:
1575 217 storres
                    sollya_lib_clear_obj(lastResPolySo)
1576 218 storres
                if not lastInfNormSo is None:
1577 218 storres
                    sollya_lib_clear_obj(lastInfNormSo)
1578 222 storres
                if debug:
1579 222 storres
                    print "Exit normally."
1580 222 storres
                    print "Precision of highest degree monomial:", itpSa
1581 222 storres
                    print "Precision of constant term          :", ftpSa
1582 218 storres
                return (resPolySo, infNormSo)
1583 217 storres
    # End wile True
1584 215 storres
    return None
1585 217 storres
# End pobyso_polynomial_coefficients_progressive_truncate_so_so.
1586 217 storres
1587 217 storres
def pobyso_polynomial_degree_so_sa(polySo):
1588 217 storres
    """
1589 217 storres
    Return the degree of a Sollya polynomial as a Sage int.
1590 217 storres
    """
1591 217 storres
    degreeSo = sollya_lib_degree(polySo)
1592 217 storres
    return pobyso_constant_from_int_so_sa(degreeSo)
1593 217 storres
# End pobyso_polynomial_degree_so_sa
1594 217 storres
1595 217 storres
def pobyso_polynomial_degree_so_so(polySo):
1596 217 storres
    """
1597 217 storres
    Thin wrapper around lib_sollya_degree().
1598 217 storres
    """
1599 217 storres
    return sollya_lib_degree(polySo)
1600 217 storres
# End pobyso_polynomial_degree_so_so
1601 217 storres
1602 5 storres
def pobyso_range(rnLowerBound, rnUpperBound):
1603 38 storres
    """ Legacy function. See pobyso_range_sa_so. """
1604 209 storres
    return pobyso_range_sa_so(rnLowerBound, rnUpperBound)
1605 38 storres
1606 226 storres
def pobyso_range_from_bounds_sa_so(rnLowerBound, rnUpperBound, precSa = None):
1607 226 storres
    """
1608 226 storres
    Create a Sollya range from 2 Sage real numbers as bounds
1609 226 storres
    """
1610 226 storres
    # TODO check precision stuff.
1611 226 storres
    sollyaPrecChanged = False
1612 226 storres
    (initialSollyaPrecSo, initialSollyaPrecSa) = \
1613 226 storres
        pobyso_get_prec_so_so_sa()
1614 226 storres
    if precSa is None:
1615 226 storres
        precSa = max(rnLowerBound.parent().prec(), rnUpperBound.parent().prec())
1616 226 storres
    if precSa > initialSollyaPrecSa:
1617 226 storres
        if precSa <= 2:
1618 226 storres
            print inspect.stack()[0][3], ": precision change <= 2 requested."
1619 226 storres
        pobyso_set_prec_sa_so(precSa)
1620 226 storres
        sollyaPrecChanged = True
1621 226 storres
    rangeSo = sollya_lib_range_from_bounds(get_rn_value(rnLowerBound),
1622 226 storres
                                           get_rn_value(rnUpperBound))
1623 226 storres
    if sollyaPrecChanged:
1624 226 storres
        pobyso_set_prec_so_so(initialSollyaPrecSo)
1625 234 storres
        sollya_lib_clear_obj(initialSollyaPrecSo)
1626 226 storres
    return rangeSo
1627 226 storres
# End pobyso_range_from_bounds_sa_so
1628 5 storres
1629 261 storres
def pobyso_range_list_so_sa(listSo):
1630 261 storres
    """
1631 261 storres
    Return a Sollya list of ranges as a Sage list of
1632 261 storres
    floating-point intervals.
1633 261 storres
    """
1634 261 storres
    listSa   = []
1635 261 storres
    ## The function returns none if the list is empty or an error has happened.
1636 261 storres
    retVal = pobyso_get_list_elements_so_so(listSo)
1637 261 storres
    if retVal is None:
1638 261 storres
        return listSa
1639 261 storres
    ## Just in case the interface is changed and an empty list is returned
1640 261 storres
    #  instead of None.
1641 261 storres
    elif len(retVal) == 0:
1642 261 storres
        return listSa
1643 261 storres
    else:
1644 261 storres
        ## Remember pobyso_get_list_elements_so_so returns more information
1645 261 storres
        #  than just the elements of the list (# elements, is_elliptic)
1646 261 storres
        listSaSo, numElements, isEndElliptic = retVal
1647 261 storres
    ## Return an empty list.
1648 261 storres
    if numElements == 0:
1649 261 storres
        return listSa
1650 261 storres
    ## Search first for the maximum precision of the elements
1651 261 storres
    maxPrecSa = 0
1652 261 storres
    for rangeSo in listSaSo:
1653 261 storres
        #pobyso_autoprint(floatSo)
1654 261 storres
        curPrecSa =  pobyso_get_prec_of_range_so_sa(rangeSo)
1655 261 storres
        if curPrecSa > maxPrecSa:
1656 261 storres
            maxPrecSa = curPrecSa
1657 261 storres
    ##
1658 261 storres
    intervalField = RealIntervalField(maxPrecSa)
1659 261 storres
    ##
1660 261 storres
    for rangeSo in listSaSo:
1661 261 storres
        listSa.append(pobyso_range_to_interval_so_sa(rangeSo, intervalField))
1662 261 storres
    return listSa
1663 261 storres
# End pobyso_range_list_so_sa
1664 261 storres
1665 226 storres
def pobyso_range_max_abs_so_so(rangeSo):
1666 226 storres
    """
1667 226 storres
    Return, as Sollya constant, the maximum absolute value of a Sollay range.
1668 226 storres
    """
1669 226 storres
    lowerBoundSo = sollya_lib_inf(rangeSo)
1670 226 storres
    upperBoundSo = sollya_lib_sup(rangeSo)
1671 226 storres
    #
1672 226 storres
    lowerBoundSo = sollya_lib_build_function_abs(lowerBoundSo)
1673 226 storres
    upperBoundSo = sollya_lib_build_function_abs(upperBoundSo)
1674 226 storres
    #pobyso_autoprint(lowerBoundSo)
1675 226 storres
    #pobyso_autoprint(upperBoundSo)
1676 226 storres
    #
1677 226 storres
    maxAbsSo = sollya_lib_max(lowerBoundSo, upperBoundSo, None)
1678 226 storres
    #sollya_lib_clear_obj(lowerBoundSo)
1679 226 storres
    #sollya_lib_clear_obj(upperBoundSo)
1680 226 storres
    return maxAbsSo
1681 226 storres
# End pobyso_range_max_abs_so_so
1682 226 storres
1683 85 storres
def pobyso_range_to_interval_so_sa(rangeSo, realIntervalFieldSa = None):
1684 83 storres
    """
1685 83 storres
    Get a Sage interval from a Sollya range.
1686 83 storres
    If no realIntervalField is given as a parameter, the Sage interval
1687 83 storres
    precision is that of the Sollya range.
1688 85 storres
    Otherwise, the precision is that of the realIntervalField. In this case
1689 85 storres
    rounding may happen.
1690 83 storres
    """
1691 85 storres
    if realIntervalFieldSa is None:
1692 56 storres
        precSa = pobyso_get_prec_of_range_so_sa(rangeSo)
1693 85 storres
        realIntervalFieldSa = RealIntervalField(precSa)
1694 56 storres
    intervalSa = \
1695 85 storres
        pobyso_get_interval_from_range_so_sa(rangeSo, realIntervalFieldSa)
1696 209 storres
    return intervalSa
1697 209 storres
# End pobyso_range_to_interval_so_sa
1698 240 storres
#
1699 240 storres
def pobyso_relative_so_so():
1700 240 storres
    """
1701 240 storres
    Very thin wrapper around the sollya_lib_relative function.
1702 240 storres
    """
1703 240 storres
    return sollya_lib_relative()
1704 240 storres
# End pobyso_relative_so_so
1705 240 storres
#
1706 209 storres
def pobyso_rat_poly_sa_so(polySa, precSa = None):
1707 209 storres
    """
1708 209 storres
    Create a Sollya polynomial from a Sage rational polynomial.
1709 272 storres
    We first convert the rational polynomial into a floating-point
1710 272 storres
    polynomial.
1711 209 storres
    """
1712 209 storres
    ## TODO: filter arguments.
1713 209 storres
    ## Precision. If no precision is given, use the current precision
1714 209 storres
    #  of Sollya.
1715 209 storres
    if precSa is None:
1716 209 storres
        precSa =  pobyso_get_prec_so_sa()
1717 209 storres
    #print "Precision:",  precSa
1718 209 storres
    RRR = RealField(precSa)
1719 209 storres
    ## Create a Sage polynomial in the "right" precision.
1720 209 storres
    P_RRR = RRR[polySa.variables()[0]]
1721 209 storres
    polyFloatSa = P_RRR(polySa)
1722 213 storres
    ## Make sure no precision is provided: pobyso_float_poly_sa_so will
1723 272 storres
    #  recover it all by itself and will not make any extra conversion.
1724 209 storres
    return pobyso_float_poly_sa_so(polyFloatSa)
1725 209 storres
1726 209 storres
# End pobyso_rat_poly_sa_so
1727 209 storres
1728 52 storres
def pobyso_remez_canonical_sa_sa(func, \
1729 52 storres
                                 degree, \
1730 52 storres
                                 lowerBound, \
1731 52 storres
                                 upperBound, \
1732 52 storres
                                 weight = None, \
1733 52 storres
                                 quality = None):
1734 52 storres
    """
1735 52 storres
    All arguments are Sage/Python.
1736 52 storres
    The functions (func and weight) must be passed as expressions or strings.
1737 52 storres
    Otherwise the function fails.
1738 83 storres
    The return value is a Sage polynomial.
1739 52 storres
    """
1740 83 storres
    var('zorglub')    # Dummy variable name for type check only. Type of
1741 83 storres
    # zorglub is "symbolic expression".
1742 52 storres
    polySo = pobyso_remez_canonical_sa_so(func, \
1743 52 storres
                                 degree, \
1744 52 storres
                                 lowerBound, \
1745 52 storres
                                 upperBound, \
1746 85 storres
                                 weight, \
1747 85 storres
                                 quality)
1748 83 storres
    # String test
1749 52 storres
    if parent(func) == parent("string"):
1750 52 storres
        functionSa = eval(func)
1751 52 storres
    # Expression test.
1752 52 storres
    elif type(func) == type(zorglub):
1753 52 storres
        functionSa = func
1754 83 storres
    else:
1755 83 storres
        return None
1756 83 storres
    #
1757 52 storres
    maxPrecision = 0
1758 52 storres
    if polySo is None:
1759 52 storres
        return(None)
1760 52 storres
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
1761 85 storres
    RRRRSa = RealField(maxPrecision)
1762 85 storres
    polynomialRingSa = RRRRSa[functionSa.variables()[0]]
1763 85 storres
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, RRRRSa)
1764 85 storres
    polySa = polynomial(expSa, polynomialRingSa)
1765 83 storres
    sollya_lib_clear_obj(polySo)
1766 52 storres
    return(polySa)
1767 85 storres
# End pobyso_remez_canonical_sa_sa
1768 52 storres
1769 38 storres
def pobyso_remez_canonical(func, \
1770 5 storres
                           degree, \
1771 5 storres
                           lowerBound, \
1772 5 storres
                           upperBound, \
1773 38 storres
                           weight = "1", \
1774 5 storres
                           quality = None):
1775 38 storres
    """ Legacy function. See pobyso_remez_canonical_sa_so. """
1776 51 storres
    return(pobyso_remez_canonical_sa_so(func, \
1777 51 storres
                                        degree, \
1778 51 storres
                                        lowerBound, \
1779 51 storres
                                        upperBound, \
1780 51 storres
                                        weight, \
1781 51 storres
                                        quality))
1782 200 storres
# End pobyso_remez_canonical.
1783 200 storres
1784 38 storres
def pobyso_remez_canonical_sa_so(func, \
1785 38 storres
                                 degree, \
1786 38 storres
                                 lowerBound, \
1787 38 storres
                                 upperBound, \
1788 52 storres
                                 weight = None, \
1789 38 storres
                                 quality = None):
1790 38 storres
    """
1791 38 storres
    All arguments are Sage/Python.
1792 51 storres
    The functions (func and weight) must be passed as expressions or strings.
1793 51 storres
    Otherwise the function fails.
1794 38 storres
    The return value is a pointer to a Sollya function.
1795 234 storres
    lowerBound and upperBound mus be reals.
1796 38 storres
    """
1797 83 storres
    var('zorglub')    # Dummy variable name for type check only. Type of
1798 83 storres
    # zorglub is "symbolic expression".
1799 85 storres
    currentVariableNameSa = None
1800 52 storres
    # The func argument can be of different types (string,
1801 52 storres
    # symbolic expression...)
1802 38 storres
    if parent(func) == parent("string"):
1803 234 storres
        localFuncSa = sage_eval(func,globals())
1804 85 storres
        if len(localFuncSa.variables()) > 0:
1805 85 storres
            currentVariableNameSa = localFuncSa.variables()[0]
1806 85 storres
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1807 159 storres
            functionSo = \
1808 159 storres
              sollya_lib_parse_string(localFuncSa._assume_str().replace('_SAGE_VAR_', ''))
1809 51 storres
    # Expression test.
1810 52 storres
    elif type(func) == type(zorglub):
1811 52 storres
        # Until we are able to translate Sage expressions into Sollya
1812 52 storres
        # expressions : parse the string version.
1813 85 storres
        if len(func.variables()) > 0:
1814 85 storres
            currentVariableNameSa = func.variables()[0]
1815 85 storres
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1816 159 storres
            functionSo = \
1817 159 storres
              sollya_lib_parse_string(func._assume_str().replace('_SAGE_VAR_', ''))
1818 38 storres
    else:
1819 38 storres
        return(None)
1820 85 storres
    if weight is None: # No weight given -> 1.
1821 52 storres
        weightSo = pobyso_constant_1_sa_so()
1822 85 storres
    elif parent(weight) == parent("string"): # Weight given as string: parse it.
1823 51 storres
        weightSo = sollya_lib_parse_string(func)
1824 85 storres
    elif type(weight) == type(zorglub): # Weight given as symbolice expression.
1825 159 storres
        functionSo = \
1826 159 storres
          sollya_lib_parse_string_sa_so(weight._assume_str().replace('_SAGE_VAR_', ''))
1827 51 storres
    else:
1828 51 storres
        return(None)
1829 5 storres
    degreeSo = pobyso_constant_from_int(degree)
1830 85 storres
    rangeSo = pobyso_bounds_to_range_sa_so(lowerBound, upperBound)
1831 38 storres
    if not quality is None:
1832 38 storres
        qualitySo= pobyso_constant_sa_so(quality)
1833 52 storres
    else:
1834 52 storres
        qualitySo = None
1835 83 storres
1836 83 storres
    remezPolySo = sollya_lib_remez(functionSo, \
1837 83 storres
                                   degreeSo, \
1838 83 storres
                                   rangeSo, \
1839 83 storres
                                   weightSo, \
1840 83 storres
                                   qualitySo, \
1841 83 storres
                                   None)
1842 83 storres
    sollya_lib_clear_obj(functionSo)
1843 83 storres
    sollya_lib_clear_obj(degreeSo)
1844 83 storres
    sollya_lib_clear_obj(rangeSo)
1845 83 storres
    sollya_lib_clear_obj(weightSo)
1846 83 storres
    if not qualitySo is None:
1847 85 storres
        sollya_lib_clear_obj(qualitySo)
1848 83 storres
    return(remezPolySo)
1849 83 storres
# End pobyso_remez_canonical_sa_so
1850 83 storres
1851 38 storres
def pobyso_remez_canonical_so_so(funcSo, \
1852 38 storres
                                 degreeSo, \
1853 38 storres
                                 rangeSo, \
1854 52 storres
                                 weightSo = pobyso_constant_1_sa_so(),\
1855 38 storres
                                 qualitySo = None):
1856 38 storres
    """
1857 38 storres
    All arguments are pointers to Sollya objects.
1858 38 storres
    The return value is a pointer to a Sollya function.
1859 38 storres
    """
1860 38 storres
    if not sollya_lib_obj_is_function(funcSo):
1861 38 storres
        return(None)
1862 38 storres
    return(sollya_lib_remez(funcSo, degreeSo, rangeSo, weightSo, qualitySo, None))
1863 200 storres
# End pobyso_remez_canonical_so_so.
1864 200 storres
1865 237 storres
def pobyso_remez_exponents_list_sa_so(func,     \
1866 237 storres
                                 exponentsList, \
1867 237 storres
                                 lowerBound,    \
1868 237 storres
                                 upperBound,    \
1869 237 storres
                                 weight = None, \
1870 237 storres
                                 quality = None):
1871 237 storres
    """
1872 237 storres
    All arguments are Sage/Python.
1873 237 storres
    The functions (func and weight) must be passed as expressions or strings.
1874 237 storres
    Otherwise the function fails.
1875 237 storres
    The return value is a pointer to a Sollya function.
1876 237 storres
    lowerBound and upperBound mus be reals.
1877 237 storres
    """
1878 237 storres
    var('zorglub')    # Dummy variable name for type check only. Type of
1879 237 storres
    # zorglub is "symbolic expression".
1880 237 storres
    currentVariableNameSa = None
1881 237 storres
    # The func argument can be of different types (string,
1882 237 storres
    # symbolic expression...)
1883 237 storres
    if parent(func) == parent("string"):
1884 237 storres
        localFuncSa = sage_eval(func,globals())
1885 237 storres
        if len(localFuncSa.variables()) > 0:
1886 237 storres
            currentVariableNameSa = localFuncSa.variables()[0]
1887 237 storres
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1888 237 storres
            functionSo = \
1889 237 storres
              sollya_lib_parse_string(localFuncSa._assume_str().replace('_SAGE_VAR_', ''))
1890 237 storres
    # Expression test.
1891 237 storres
    elif type(func) == type(zorglub):
1892 237 storres
        # Until we are able to translate Sage expressions into Sollya
1893 237 storres
        # expressions : parse the string version.
1894 237 storres
        if len(func.variables()) > 0:
1895 237 storres
            currentVariableNameSa = func.variables()[0]
1896 237 storres
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1897 237 storres
            functionSo = \
1898 237 storres
              sollya_lib_parse_string(func._assume_str().replace('_SAGE_VAR_', ''))
1899 237 storres
    else:
1900 237 storres
        return(None)
1901 237 storres
    ## Deal with the weight, much in the same way as with the function.
1902 237 storres
    if weight is None: # No weight given -> 1.
1903 237 storres
        weightSo = pobyso_constant_1_sa_so()
1904 237 storres
    elif parent(weight) == parent("string"): # Weight given as string: parse it.
1905 237 storres
        weightSo = sollya_lib_parse_string(func)
1906 237 storres
    elif type(weight) == type(zorglub): # Weight given as symbolice expression.
1907 237 storres
        functionSo = \
1908 237 storres
          sollya_lib_parse_string(weight._assume_str().replace('_SAGE_VAR_', '',100))
1909 237 storres
    else:
1910 237 storres
        return(None)
1911 237 storres
    rangeSo = pobyso_bounds_to_range_sa_so(lowerBound, upperBound)
1912 237 storres
    if not quality is None:
1913 237 storres
        qualitySo= pobyso_constant_sa_so(quality)
1914 237 storres
    else:
1915 237 storres
        qualitySo = None
1916 237 storres
    #
1917 237 storres
    ## Tranform the Sage list of exponents into a Sollya list.
1918 237 storres
    exponentsListSo = pobyso_build_list_of_ints_sa_so(*exponentsList)
1919 278 storres
    remezPolySo = sollya_lib_remez(functionSo,      \
1920 237 storres
                                   exponentsListSo, \
1921 278 storres
                                   rangeSo,         \
1922 278 storres
                                   weightSo,        \
1923 278 storres
                                   qualitySo,       \
1924 237 storres
                                   None)
1925 237 storres
    sollya_lib_clear_obj(functionSo)
1926 237 storres
    sollya_lib_clear_obj(exponentsListSo)
1927 237 storres
    sollya_lib_clear_obj(rangeSo)
1928 237 storres
    sollya_lib_clear_obj(weightSo)
1929 237 storres
    if not qualitySo is None:
1930 237 storres
        sollya_lib_clear_obj(qualitySo)
1931 237 storres
    return(remezPolySo)
1932 237 storres
# End pobyso_remez_exponentsList_sa_so
1933 240 storres
#
1934 240 storres
def pobyso_round_coefficients_so_so(polySo, truncFormatListSo):
1935 240 storres
    """
1936 240 storres
    A wrapper around the "classical" sollya_lib_roundcoefficients: a Sollya
1937 240 storres
    polynomial and a Sollya list are given as arguments.
1938 240 storres
    """
1939 240 storres
    return sollya_lib_roundcoefficients(polySo, truncFormatListSo)
1940 237 storres
1941 224 storres
def pobyso_round_coefficients_progressive_so_so(polySo,
1942 224 storres
                                                funcSo,
1943 224 storres
                                                precSo,
1944 224 storres
                                                intervalSo,
1945 224 storres
                                                icSo,
1946 224 storres
                                                currentApproxErrorSo,
1947 224 storres
                                                approxAccurSo,
1948 224 storres
                                                debug=False):
1949 272 storres
    """
1950 272 storres
    From an input approximation polynomial, compute an output one with
1951 272 storres
    smaller coefficients and yet yields a sufficient approximation accuracy.
1952 272 storres
    """
1953 223 storres
    if debug:
1954 223 storres
        print "Input arguments:"
1955 228 storres
        print "Polynomial: ", ; pobyso_autoprint(polySo)
1956 228 storres
        print "Function: ", ; pobyso_autoprint(funcSo)
1957 228 storres
        print "Internal precision: ", ; pobyso_autoprint(precSo)
1958 228 storres
        print "Interval: ", ; pobyso_autoprint(intervalSo)
1959 228 storres
        print "Current approximation error: ", ; pobyso_autoprint(currentApproxErrorSo)
1960 228 storres
        print "Requested approxiation error: ", ; pobyso_autoprint(approxAccurSo)
1961 223 storres
        print "________________"
1962 223 storres
    approxAccurSa        = pobyso_get_constant_as_rn_so_sa(approxAccurSo)
1963 223 storres
    currentApproxErrorSa = pobyso_get_constant_as_rn_so_sa(currentApproxErrorSo)
1964 223 storres
    ## If the current approximation error is too close to the target, there is
1965 223 storres
    #  no possible gain.
1966 223 storres
    if currentApproxErrorSa >= approxAccurSa / 2:
1967 232 storres
        #### Do not return the initial argument but copies: caller may free
1968 272 storres
        #    the former as inutile after call.
1969 232 storres
        return (sollya_lib_copy_obj(polySo),
1970 232 storres
                sollya_lib_copy_obj(currentApproxErrorSo))
1971 278 storres
    #
1972 278 storres
    ## Try to round the coefficients.
1973 223 storres
    degreeSa             = pobyso_polynomial_degree_so_sa(polySo)
1974 223 storres
    intervalSa           = pobyso_range_to_interval_so_sa(intervalSo)
1975 223 storres
1976 223 storres
    if debug:
1977 228 storres
        print "degreeSa              :", degreeSa
1978 228 storres
        print "intervalSa            :", intervalSa.str(style='brackets')
1979 223 storres
        print "currentApproxErrorSa  :", currentApproxErrorSa
1980 228 storres
        print "approxAccurSa         :", approxAccurSa
1981 223 storres
    radiusSa = intervalSa.absolute_diameter() / 2
1982 223 storres
    if debug:
1983 224 storres
        print "log2(radius):", RR(radiusSa).log2()
1984 224 storres
    iterIndex = 0
1985 272 storres
    ## Build the "shaved" polynomial.
1986 224 storres
    while True:
1987 272 storres
        ### Start with a 0 value expression.
1988 224 storres
        resPolySo = pobyso_constant_0_sa_so()
1989 224 storres
        roundedPolyApproxAccurSa = approxAccurSa / 2
1990 224 storres
        currentRadiusPowerSa = 1
1991 224 storres
        for degree in xrange(0,degreeSa + 1):
1992 224 storres
            #### At round 0, use the agressive formula. At round 1, run the
1993 224 storres
            #    proved formula.
1994 224 storres
            if iterIndex == 0:
1995 224 storres
                roundingPowerSa = \
1996 224 storres
                    floor(((currentRadiusPowerSa/roundedPolyApproxAccurSa)*(degree+1)).log2())
1997 224 storres
            else:
1998 224 storres
                roundingPowerSa = \
1999 224 storres
                    floor(((currentRadiusPowerSa/roundedPolyApproxAccurSa)*(degreeSa+1)).log2())
2000 272 storres
            ## Under extreme conditions the above formulas can evaluate under 2,
2001 272 storres
            #   which is the minimal precision of an MPFR number.
2002 228 storres
            if roundingPowerSa < 2:
2003 228 storres
                roundingPowerSa = 2
2004 224 storres
            if debug:
2005 224 storres
                print "roundedPolyApproxAccurSa", roundedPolyApproxAccurSa
2006 224 storres
                print "currentRadiusPowerSa", currentRadiusPowerSa
2007 224 storres
                print "Current rounding exponent:", roundingPowerSa
2008 224 storres
            currentRadiusPowerSa *= radiusSa
2009 224 storres
            index1So = pobyso_constant_from_int_sa_so(degree)
2010 224 storres
            index2So = pobyso_constant_from_int_sa_so(degree)
2011 224 storres
            ### Create a monomial with:
2012 224 storres
            #   - the coefficient in the initial monomial at the current degrree;
2013 224 storres
            #   - the current exponent;
2014 224 storres
            #   - the free variable.
2015 224 storres
            cmonSo  = \
2016 224 storres
                sollya_lib_build_function_mul(sollya_lib_coeff(polySo, index1So),
2017 224 storres
                                              sollya_lib_build_function_pow( \
2018 224 storres
                                                sollya_lib_build_function_free_variable(), \
2019 224 storres
                                                index2So))
2020 224 storres
            roundingPowerSo = pobyso_constant_from_int_sa_so(roundingPowerSa)
2021 224 storres
            cmonrSo = pobyso_round_coefficients_single_so_so(cmonSo, roundingPowerSo)
2022 224 storres
            sollya_lib_clear_obj(cmonSo)
2023 224 storres
            ### Add to the result polynomial.
2024 224 storres
            resPolySo = sollya_lib_build_function_add(resPolySo,
2025 224 storres
                                                      cmonrSo)
2026 224 storres
        # End for.
2027 224 storres
        ### Check the new polynomial.
2028 224 storres
        freeVarSo     = sollya_lib_build_function_free_variable()
2029 224 storres
        changeVarSo   = sollya_lib_sub(freeVarSo, icSo)
2030 224 storres
        resPolyCvSo   = sollya_lib_evaluate(resPolySo, changeVarSo)
2031 224 storres
        errFuncSo     = sollya_lib_build_function_sub(sollya_lib_copy_obj(funcSo),
2032 224 storres
                                                      resPolyCvSo)
2033 224 storres
        infNormSo = sollya_lib_dirtyinfnorm(errFuncSo, intervalSo)
2034 224 storres
        ### This also clears resPolyCvSo.
2035 224 storres
        sollya_lib_clear_obj(errFuncSo)
2036 224 storres
        cerrSa    = pobyso_get_constant_as_rn_so_sa(infNormSo)
2037 223 storres
        if debug:
2038 224 storres
            print "Error of the new polynomial:", cerrSa
2039 224 storres
        ### If at round 1, return the initial polynomial error. This should
2040 224 storres
        #   never happen since the rounding algorithm is proved. But some
2041 224 storres
        #   circumstances may break it (e.g. internal precision of tools).
2042 224 storres
        if cerrSa > approxAccurSa:
2043 278 storres
            if iterIndex == 0: # Round 0 is agressive rounding, got round 1 (proved rounding)
2044 224 storres
                sollya_lib_clear_obj(resPolySo)
2045 224 storres
                sollya_lib_clear_obj(infNormSo)
2046 278 storres
                iterIndex += 1
2047 278 storres
                continue
2048 278 storres
            else: # Round 1 and beyond : just return the oroginal polynomial.
2049 278 storres
                sollya_lib_clear_obj(resPolySo)
2050 278 storres
                sollya_lib_clear_obj(infNormSo)
2051 232 storres
                #### Do not return the arguments but copies: the caller may free
2052 232 storres
                #    free the former as inutile after call.
2053 232 storres
                return (sollya_lib_copy_obj(polySo),
2054 232 storres
                        sollya_lib_copy_obj(currentApproxErrorSo))
2055 224 storres
        ### If get here it is because cerrSa <= approxAccurSa
2056 224 storres
        ### Approximation error of the new polynomial is acceptable.
2057 224 storres
        return (resPolySo, infNormSo)
2058 224 storres
    # End while True
2059 224 storres
# End pobyso_round_coefficients_progressive_so_so
2060 223 storres
2061 240 storres
def pobyso_round_coefficients_single_so_so(polySo, commonPrecSo):
2062 215 storres
    """
2063 215 storres
    Create a rounded coefficients polynomial from polynomial argument to
2064 215 storres
    the number of bits in size argument.
2065 215 storres
    All coefficients are set to the same precision.
2066 215 storres
    """
2067 215 storres
    ## TODO: check arguments.
2068 240 storres
    endEllipListSo = pobyso_build_end_elliptic_list_so_so(commonPrecSo)
2069 215 storres
    polySo = sollya_lib_roundcoefficients(polySo, endEllipListSo, None)
2070 217 storres
    sollya_lib_clear_obj(endEllipListSo)
2071 215 storres
    #sollya_lib_clear_obj(endEllipListSo)
2072 215 storres
    return polySo
2073 215 storres
2074 215 storres
# End pobyso_round_coefficients_single_so_so
2075 215 storres
2076 5 storres
def pobyso_set_canonical_off():
2077 5 storres
    sollya_lib_set_canonical(sollya_lib_off())
2078 5 storres
2079 5 storres
def pobyso_set_canonical_on():
2080 5 storres
    sollya_lib_set_canonical(sollya_lib_on())
2081 5 storres
2082 5 storres
def pobyso_set_prec(p):
2083 38 storres
    """ Legacy function. See pobyso_set_prec_sa_so. """
2084 85 storres
    pobyso_set_prec_sa_so(p)
2085 38 storres
2086 38 storres
def pobyso_set_prec_sa_so(p):
2087 227 storres
    #a = c_int(p)
2088 226 storres
    #precSo = c_void_p(sollya_lib_constant_from_int(a))
2089 227 storres
    #precSo = sollya_lib_constant_from_int(a)
2090 227 storres
    precSo = pobyso_constant_from_int_sa_so(p)
2091 226 storres
    sollya_lib_set_prec(precSo)
2092 226 storres
    sollya_lib_clear_obj(precSo)
2093 215 storres
# End pobyso_set_prec_sa_so.
2094 5 storres
2095 85 storres
def pobyso_set_prec_so_so(newPrecSo):
2096 226 storres
    sollya_lib_set_prec(newPrecSo)
2097 215 storres
# End pobyso_set_prec_so_so.
2098 242 storres
#
2099 281 storres
def pobyso_supnorm_sa_sa(poly):
2100 281 storres
    """
2101 281 storres
    Computes the supremum norm from Sage input arguments and returns a
2102 281 storres
    Sage floating-point number whose precision is set by the realFieldSa
2103 281 storres
    argument.
2104 281 storres
    TODO: complete this stub!
2105 281 storres
    """
2106 281 storres
    print("This function does nothing!")
2107 281 storres
    return None
2108 281 storres
# End pobyso_supnorm_sa_sa
2109 281 storres
2110 242 storres
def pobyso_supnorm_so_sa(polySo, funcSo, intervalSo, errorTypeSo = None,\
2111 242 storres
                         accuracySo = None, realFieldSa = None):
2112 242 storres
    """
2113 281 storres
    Computes the supremum norm from Sollya input arguments and returns a
2114 242 storres
    Sage floating-point number whose precision is set by the only Sage argument.
2115 215 storres

2116 242 storres
    The returned value is the maximum of the absolute values of the range
2117 242 storres
    elements  returned by the Sollya supnorm functions
2118 242 storres
    """
2119 242 storres
    supNormRangeSo = pobyso_supnorm_so_so(polySo,
2120 242 storres
                                          funcSo,
2121 242 storres
                                          intervalSo,
2122 242 storres
                                          errorTypeSo,
2123 242 storres
                                          accuracySo)
2124 242 storres
    supNormSo = pobyso_range_max_abs_so_so(supNormRangeSo)
2125 242 storres
    sollya_lib_clear_obj(supNormRangeSo)
2126 242 storres
    #pobyso_autoprint(supNormSo)
2127 242 storres
    supNormSa = pobyso_get_constant_as_rn_with_rf_so_sa(supNormSo, realFieldSa)
2128 242 storres
    sollya_lib_clear_obj(supNormSo)
2129 242 storres
    return supNormSa
2130 242 storres
# End pobyso_supnorm_so_sa.
2131 242 storres
#
2132 85 storres
def pobyso_supnorm_so_so(polySo, funcSo, intervalSo, errorTypeSo = None,\
2133 85 storres
                         accuracySo = None):
2134 58 storres
    """
2135 85 storres
    Computes the supnorm of the approximation error between the given
2136 242 storres
    polynomial and function. Attention: returns a range!
2137 85 storres
    errorTypeSo defaults to "absolute".
2138 85 storres
    accuracySo defaults to 2^(-40).
2139 85 storres
    """
2140 85 storres
    if errorTypeSo is None:
2141 85 storres
        errorTypeSo = sollya_lib_absolute(None)
2142 85 storres
        errorTypeIsNone = True
2143 85 storres
    else:
2144 85 storres
        errorTypeIsNone = False
2145 85 storres
    #
2146 85 storres
    if accuracySo is None:
2147 237 storres
        # Notice the **: we are in Pythonland here!
2148 85 storres
        accuracySo = pobyso_constant_sa_so(RR(2**(-40)))
2149 85 storres
        accuracyIsNone = True
2150 85 storres
    else:
2151 85 storres
        accuracyIsNone = False
2152 238 storres
    #pobyso_autoprint(accuracySo)
2153 85 storres
    resultSo = \
2154 85 storres
        sollya_lib_supnorm(polySo, funcSo, intervalSo, errorTypeSo, \
2155 85 storres
                              accuracySo)
2156 85 storres
    if errorTypeIsNone:
2157 85 storres
        sollya_lib_clear_obj(errorTypeSo)
2158 85 storres
    if accuracyIsNone:
2159 85 storres
        sollya_lib_clear_obj(accuracySo)
2160 85 storres
    return resultSo
2161 85 storres
# End pobyso_supnorm_so_so
2162 242 storres
#
2163 162 storres
def pobyso_taylor_expansion_no_change_var_so_so(functionSo,
2164 162 storres
                                                degreeSo,
2165 162 storres
                                                rangeSo,
2166 162 storres
                                                errorTypeSo=None,
2167 162 storres
                                                sollyaPrecSo=None):
2168 85 storres
    """
2169 162 storres
    Compute the Taylor expansion without the variable change
2170 162 storres
    x -> x-intervalCenter.
2171 272 storres
    If errorTypeSo is None, absolute is used.
2172 272 storres
    If sollyaPrecSo is None, Sollya internal precision is not changed.
2173 58 storres
    """
2174 226 storres
    # Change internal Sollya precision, if needed.
2175 227 storres
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
2176 226 storres
    sollyaPrecChanged = False
2177 226 storres
    if sollyaPrecSo is None:
2178 226 storres
        pass
2179 226 storres
    else:
2180 58 storres
        sollya_lib_set_prec(sollyaPrecSo)
2181 226 storres
        sollyaPrecChanged = True
2182 85 storres
    # Error type stuff: default to absolute.
2183 85 storres
    if errorTypeSo is None:
2184 85 storres
        errorTypeIsNone = True
2185 85 storres
        errorTypeSo = sollya_lib_absolute(None)
2186 85 storres
    else:
2187 85 storres
        errorTypeIsNone = False
2188 162 storres
    intervalCenterSo = sollya_lib_mid(rangeSo, None)
2189 162 storres
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo,
2190 162 storres
                                         intervalCenterSo,
2191 58 storres
                                         rangeSo, errorTypeSo, None)
2192 272 storres
    # Object taylorFormListSaSo is a Python list of Sollya objects references
2193 272 storres
    #  that are copies of the elements of taylorFormSo.
2194 117 storres
    # pobyso_get_list_elements_so_so clears taylorFormSo.
2195 162 storres
    (taylorFormListSaSo, numElementsSa, isEndEllipticSa) = \
2196 58 storres
        pobyso_get_list_elements_so_so(taylorFormSo)
2197 272 storres
    ## Copy needed here since polySo will be returned and taylorFormListSaSo
2198 272 storres
    #  will be cleared.
2199 162 storres
    polySo = sollya_lib_copy_obj(taylorFormListSaSo[0])
2200 162 storres
    #print "Num elements:", numElementsSa
2201 162 storres
    sollya_lib_clear_obj(taylorFormSo)
2202 181 storres
    # No copy_obj needed here: a new objects are created.
2203 272 storres
    maxErrorSo    = sollya_lib_sup(taylorFormListSaSo[2])
2204 272 storres
    minErrorSo    = sollya_lib_inf(taylorFormListSaSo[2])
2205 272 storres
    # List taylorFormListSaSo is not needed anymore.
2206 280 storres
    pobyso_clear_list_elements_sa_so(taylorFormListSaSo)
2207 181 storres
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
2208 181 storres
    absMinErrorSo = sollya_lib_abs(minErrorSo)
2209 181 storres
    sollya_lib_clear_obj(maxErrorSo)
2210 181 storres
    sollya_lib_clear_obj(minErrorSo)
2211 181 storres
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
2212 181 storres
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
2213 272 storres
    #
2214 272 storres
    if errorTypeIsNone:
2215 272 storres
        sollya_lib_clear_obj(errorTypeSo)
2216 272 storres
    ## If changed, reset the Sollya working precision.
2217 226 storres
    if sollyaPrecChanged:
2218 226 storres
        sollya_lib_set_prec(initialSollyaPrecSo)
2219 226 storres
    sollya_lib_clear_obj(initialSollyaPrecSo)
2220 272 storres
    ## According to what error is the largest, return the errors.
2221 181 storres
    if absMaxErrorSa > absMinErrorSa:
2222 181 storres
        sollya_lib_clear_obj(absMinErrorSo)
2223 227 storres
        return (polySo, intervalCenterSo, absMaxErrorSo)
2224 181 storres
    else:
2225 181 storres
        sollya_lib_clear_obj(absMaxErrorSo)
2226 227 storres
        return (polySo, intervalCenterSo, absMinErrorSo)
2227 162 storres
# end pobyso_taylor_expansion_no_change_var_so_so
2228 58 storres
2229 162 storres
def pobyso_taylor_expansion_with_change_var_so_so(functionSo, degreeSo, \
2230 162 storres
                                                  rangeSo, \
2231 162 storres
                                                  errorTypeSo=None, \
2232 162 storres
                                                  sollyaPrecSo=None):
2233 58 storres
    """
2234 162 storres
    Compute the Taylor expansion with the variable change
2235 162 storres
    x -> (x-intervalCenter) included.
2236 58 storres
    """
2237 226 storres
    # Change Sollya internal precision, if need.
2238 226 storres
    sollyaPrecChanged = False
2239 271 storres
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
2240 226 storres
    if sollyaPrecSo is None:
2241 226 storres
        pass
2242 226 storres
    else:
2243 56 storres
        sollya_lib_set_prec(sollyaPrecSo)
2244 226 storres
        sollyaPrecChanged = True
2245 162 storres
    #
2246 85 storres
    # Error type stuff: default to absolute.
2247 85 storres
    if errorTypeSo is None:
2248 85 storres
        errorTypeIsNone = True
2249 85 storres
        errorTypeSo = sollya_lib_absolute(None)
2250 85 storres
    else:
2251 85 storres
        errorTypeIsNone = False
2252 162 storres
    intervalCenterSo = sollya_lib_mid(rangeSo)
2253 162 storres
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo, \
2254 162 storres
                                         intervalCenterSo, \
2255 56 storres
                                         rangeSo, errorTypeSo, None)
2256 272 storres
    # Object taylorFormListSaSo is a Python list of Sollya objects references
2257 272 storres
    # that are copies of the elements of taylorFormSo.
2258 116 storres
    # pobyso_get_list_elements_so_so clears taylorFormSo.
2259 272 storres
    (taylorFormListSaSo, numElements, isEndElliptic) = \
2260 56 storres
        pobyso_get_list_elements_so_so(taylorFormSo)
2261 272 storres
    sollya_lib_clear_obj(taylorFormSo)
2262 272 storres
    polySo = taylorFormListSaSo[0]
2263 272 storres
    ## Maximum error computation with taylorFormListSaSo[2], a range
2264 272 storres
    #  holding the actual error. Bounds can be negative.
2265 272 storres
    maxErrorSo    = sollya_lib_sup(taylorFormListSaSo[2])
2266 272 storres
    minErrorSo    = sollya_lib_inf(taylorFormListSaSo[2])
2267 181 storres
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
2268 181 storres
    absMinErrorSo = sollya_lib_abs(minErrorSo)
2269 181 storres
    sollya_lib_clear_obj(maxErrorSo)
2270 181 storres
    sollya_lib_clear_obj(minErrorSo)
2271 181 storres
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
2272 181 storres
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
2273 162 storres
    changeVarExpSo = sollya_lib_build_function_sub(\
2274 162 storres
                       sollya_lib_build_function_free_variable(),\
2275 162 storres
                       sollya_lib_copy_obj(intervalCenterSo))
2276 181 storres
    polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpSo)
2277 272 storres
    # List taylorFormListSaSo is not needed anymore.
2278 280 storres
    pobyso_clear_list_elements_sa_so(taylorFormListSaSo)
2279 162 storres
    sollya_lib_clear_obj(changeVarExpSo)
2280 56 storres
    # If changed, reset the Sollya working precision.
2281 226 storres
    if sollyaPrecChanged:
2282 226 storres
        sollya_lib_set_prec(initialSollyaPrecSo)
2283 226 storres
    sollya_lib_clear_obj(initialSollyaPrecSo)
2284 85 storres
    if errorTypeIsNone:
2285 85 storres
        sollya_lib_clear_obj(errorTypeSo)
2286 162 storres
    # Do not clear maxErrorSo.
2287 181 storres
    if absMaxErrorSa > absMinErrorSa:
2288 181 storres
        sollya_lib_clear_obj(absMinErrorSo)
2289 181 storres
        return((polyVarChangedSo, intervalCenterSo, absMaxErrorSo))
2290 181 storres
    else:
2291 181 storres
        sollya_lib_clear_obj(absMaxErrorSo)
2292 181 storres
        return((polyVarChangedSo, intervalCenterSo, absMinErrorSo))
2293 162 storres
# end pobyso_taylor_expansion_with_change_var_so_so
2294 56 storres
2295 5 storres
def pobyso_taylor(function, degree, point):
2296 38 storres
    """ Legacy function. See pobysoTaylor_so_so. """
2297 38 storres
    return(pobyso_taylor_so_so(function, degree, point))
2298 38 storres
2299 56 storres
def pobyso_taylor_so_so(functionSo, degreeSo, pointSo):
2300 56 storres
    return(sollya_lib_taylor(functionSo, degreeSo, pointSo))
2301 5 storres
2302 85 storres
def pobyso_taylorform(function, degree, point = None,
2303 85 storres
                      interval = None, errorType=None):
2304 85 storres
    """ Legacy function. See pobyso_taylorform_sa_sa;"""
2305 38 storres
2306 38 storres
def pobyso_taylorform_sa_sa(functionSa, \
2307 84 storres
                            degreeSa, \
2308 84 storres
                            pointSa, \
2309 84 storres
                            intervalSa=None, \
2310 84 storres
                            errorTypeSa=None, \
2311 84 storres
                            precisionSa=None):
2312 37 storres
    """
2313 85 storres
    Compute the Taylor form of 'degreeSa' for 'functionSa' at 'pointSa'
2314 85 storres
    for 'intervalSa' with 'errorTypeSa' (a string) using 'precisionSa'.
2315 37 storres
    point: must be a Real or a Real interval.
2316 37 storres
    return the Taylor form as an array
2317 83 storres
    TODO: take care of the interval and of the point when it is an interval;
2318 38 storres
          when errorType is not None;
2319 83 storres
          take care of the other elements of the Taylor form (coefficients
2320 83 storres
          errors and delta.
2321 37 storres
    """
2322 37 storres
    # Absolute as the default error.
2323 84 storres
    if errorTypeSa is None:
2324 37 storres
        errorTypeSo = sollya_lib_absolute()
2325 84 storres
    elif errorTypeSa == "relative":
2326 84 storres
        errorTypeSo = sollya_lib_relative()
2327 84 storres
    elif errortypeSa == "absolute":
2328 84 storres
        errorTypeSo = sollya_lib_absolute()
2329 37 storres
    else:
2330 84 storres
        # No clean up needed.
2331 84 storres
        return None
2332 84 storres
    # Global precision stuff
2333 226 storres
    sollyaPrecisionChangedSa = False
2334 226 storres
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
2335 226 storres
    if precisionSa is None:
2336 226 storres
        precSa = initialSollyaPrecSa
2337 226 storres
    else:
2338 226 storres
        if precSa > initialSollyaPrecSa:
2339 226 storres
            if precSa <= 2:
2340 226 storres
                print inspect.stack()[0][3], ":precision change <= 2 requested."
2341 226 storres
            pobyso_set_prec_sa_so(precSa)
2342 226 storres
            sollyaPrecisionChangedSa = True
2343 226 storres
    #
2344 85 storres
    if len(functionSa.variables()) > 0:
2345 85 storres
        varSa = functionSa.variables()[0]
2346 85 storres
        pobyso_name_free_variable_sa_so(str(varSa))
2347 84 storres
    # In any case (point or interval) the parent of pointSa has a precision
2348 84 storres
    # method.
2349 84 storres
    pointPrecSa = pointSa.parent().precision()
2350 226 storres
    if precSa > pointPrecSa:
2351 226 storres
        pointPrecSa = precSa
2352 84 storres
    # In any case (point or interval) pointSa has a base_ring() method.
2353 84 storres
    pointBaseRingString = str(pointSa.base_ring())
2354 84 storres
    if re.search('Interval', pointBaseRingString) is None: # Point
2355 84 storres
        pointSo = pobyso_constant_sa_so(pointSa, pointPrecSa)
2356 84 storres
    else: # Interval.
2357 84 storres
        pointSo = pobyso_interval_to_range_sa_so(pointSa, pointPrecSa)
2358 37 storres
    # Sollyafy the function.
2359 159 storres
    functionSo = pobyso_parse_string_sa_so(functionSa._assume_str().replace('_SAGE_VAR_', ''))
2360 37 storres
    if sollya_lib_obj_is_error(functionSo):
2361 37 storres
        print "pobyso_tailorform: function string can't be parsed!"
2362 37 storres
        return None
2363 37 storres
    # Sollyafy the degree
2364 84 storres
    degreeSo = sollya_lib_constant_from_int(int(degreeSa))
2365 37 storres
    # Sollyafy the point
2366 37 storres
    # Call Sollya
2367 83 storres
    taylorFormSo = \
2368 83 storres
        sollya_lib_taylorform(functionSo, degreeSo, pointSo, errorTypeSo,\
2369 37 storres
                                         None)
2370 85 storres
    sollya_lib_clear_obj(functionSo)
2371 85 storres
    sollya_lib_clear_obj(degreeSo)
2372 85 storres
    sollya_lib_clear_obj(pointSo)
2373 85 storres
    sollya_lib_clear_obj(errorTypeSo)
2374 38 storres
    (tfsAsList, numElements, isEndElliptic) = \
2375 38 storres
            pobyso_get_list_elements_so_so(taylorFormSo)
2376 37 storres
    polySo = tfsAsList[0]
2377 38 storres
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
2378 37 storres
    polyRealField = RealField(maxPrecision)
2379 38 storres
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, polyRealField)
2380 226 storres
    if sollyaPrecisionChangedSa:
2381 226 storres
        sollya_lib_set_prec(initialSollyaPrecSo)
2382 226 storres
    sollya_lib_clear_obj(initialSollyaPrecSo)
2383 37 storres
    polynomialRing = polyRealField[str(varSa)]
2384 37 storres
    polySa = polynomial(expSa, polynomialRing)
2385 37 storres
    taylorFormSa = [polySa]
2386 85 storres
    # Final clean-up.
2387 85 storres
    sollya_lib_clear_obj(taylorFormSo)
2388 51 storres
    return(taylorFormSa)
2389 51 storres
# End pobyso_taylor_form_sa_sa
2390 54 storres
2391 54 storres
def pobyso_taylorform_so_so(functionSo, degreeSo, pointSo, intervalSo=None, \
2392 54 storres
                            errorTypeSo=None):
2393 54 storres
    createdErrorType = False
2394 51 storres
    if errorTypeSo is None:
2395 51 storres
        errorTypeSo = sollya_lib_absolute()
2396 54 storres
        createdErrorType = True
2397 51 storres
    else:
2398 51 storres
        #TODO: deal with the other case.
2399 51 storres
        pass
2400 51 storres
    if intervalSo is None:
2401 54 storres
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
2402 54 storres
                                         errorTypeSo, None)
2403 51 storres
    else:
2404 54 storres
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
2405 54 storres
                                         intervalSo, errorTypeSo, None)
2406 54 storres
    if createdErrorType:
2407 54 storres
        sollya_lib_clear_obj(errorTypeSo)
2408 215 storres
    return resultSo
2409 51 storres
2410 37 storres
2411 37 storres
def pobyso_univar_polynomial_print_reverse(polySa):
2412 51 storres
    """ Legacy function. See pobyso_univar_polynomial_print_reverse_sa_sa. """
2413 51 storres
    return(pobyso_univar_polynomial_print_reverse_sa_sa(polySa))
2414 38 storres
2415 51 storres
def pobyso_univar_polynomial_print_reverse_sa_sa(polySa):
2416 37 storres
    """
2417 37 storres
    Return the string representation of a univariate polynomial with
2418 38 storres
    monomials ordered in the x^0..x^n order of the monomials.
2419 37 storres
    Remember: Sage
2420 37 storres
    """
2421 37 storres
    polynomialRing = polySa.base_ring()
2422 37 storres
    # A very expensive solution:
2423 37 storres
    # -create a fake multivariate polynomial field with only one variable,
2424 37 storres
    #   specifying a negative lexicographical order;
2425 37 storres
    mpolynomialRing = PolynomialRing(polynomialRing.base(), \
2426 37 storres
                                     polynomialRing.variable_name(), \
2427 37 storres
                                     1, order='neglex')
2428 37 storres
    # - convert the univariate argument polynomial into a multivariate
2429 37 storres
    #   version;
2430 37 storres
    p = mpolynomialRing(polySa)
2431 37 storres
    # - return the string representation of the converted form.
2432 37 storres
    # There is no simple str() method defined for p's class.
2433 37 storres
    return(p.__str__())
2434 5 storres
#
2435 230 storres
#print pobyso_get_prec()
2436 5 storres
pobyso_set_prec(165)
2437 230 storres
#print pobyso_get_prec()
2438 230 storres
#a=100
2439 230 storres
#print type(a)
2440 230 storres
#id(a)
2441 230 storres
#print "Max arity: ", pobyso_max_arity
2442 230 storres
#print "Function tripleDouble (43) as a string: ", pobyso_function_type_as_string(43)
2443 230 storres
#print "Function None (44) as a string: ", pobyso_function_type_as_string(44)
2444 246 storres
sys.stderr.write("\t...Pobyso check done.\n")