Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 261

Historique | Voir | Annoter | Télécharger (87,63 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 54 storres
103 84 storres
def pobyso_bounds_to_range_sa_so(rnLowerBoundSa, rnUpperBoundSa, \
104 84 storres
                                 precisionSa=None):
105 84 storres
    """
106 84 storres
    Return a Sollya range from to 2 RealField Sage elements.
107 84 storres
    The Sollya range element has a sufficient precision to hold all
108 162 storres
    the digits of the widest of the Sage bounds.
109 84 storres
    """
110 84 storres
    # Sanity check.
111 84 storres
    if rnLowerBoundSa > rnUpperBoundSa:
112 84 storres
        return None
113 115 storres
    # Precision stuff.
114 85 storres
    if precisionSa is None:
115 84 storres
        # Check for the largest precision.
116 84 storres
        lbPrecSa = rnLowerBoundSa.parent().precision()
117 84 storres
        ubPrecSa = rnLowerBoundSa.parent().precision()
118 84 storres
        maxPrecSa = max(lbPrecSa, ubPrecSa)
119 84 storres
    else:
120 84 storres
        maxPrecSa = precisionSa
121 115 storres
    # From Sage to Sollya bounds.
122 162 storres
#    lowerBoundSo = sollya_lib_constant(get_rn_value(rnLowerBoundSa),
123 162 storres
#                                       maxPrecSa)
124 162 storres
    lowerBoundSo = pobyso_constant_sa_so(rnLowerBoundSa,
125 162 storres
                                         maxPrecSa)
126 162 storres
    upperBoundSo = pobyso_constant_sa_so(rnUpperBoundSa,
127 154 storres
                                       maxPrecSa)
128 162 storres
129 115 storres
    # From Sollya bounds to range.
130 84 storres
    rangeSo = sollya_lib_range(lowerBoundSo, upperBoundSo)
131 84 storres
    # Back to original precision.
132 84 storres
    # Clean up
133 84 storres
    sollya_lib_clear_obj(lowerBoundSo)
134 84 storres
    sollya_lib_clear_obj(upperBoundSo)
135 154 storres
    return rangeSo
136 84 storres
# End pobyso_bounds_to_range_sa_so
137 84 storres
138 215 storres
def pobyso_build_end_elliptic_list_so_so(*args):
139 215 storres
    """
140 237 storres
    From Sollya argument objects, create a Sollya end elliptic list.
141 237 storres
    Elements of the list are "eaten" (should not be cleared individually,
142 215 storres
    are cleared when the list is cleared).
143 215 storres
    """
144 215 storres
    if len(args) == 0:
145 215 storres
        ## Called with an empty list produced "error".
146 215 storres
        return sollya_lib_build_end_elliptic_list(None)
147 215 storres
    index = 0
148 215 storres
    ## One can not append elements to an elliptic list, prepend only is
149 215 storres
    #  permitted.
150 215 storres
    for argument in reversed(args):
151 215 storres
        if index == 0:
152 215 storres
            listSo = sollya_lib_build_end_elliptic_list(argument, None)
153 215 storres
        else:
154 215 storres
            listSo = sollya_lib_prepend(argument, listSo)
155 215 storres
        index += 1
156 215 storres
    return listSo
157 215 storres
158 215 storres
# End pobyso_build_end_elliptic_list_so_so
159 215 storres
160 54 storres
def pobyso_build_function_sub_so_so(exp1So, exp2So):
161 237 storres
    return sollya_lib_build_function_sub(exp1So, exp2So)
162 54 storres
163 237 storres
def pobyso_build_list_of_ints_sa_so(*args):
164 237 storres
    """
165 237 storres
    Build a Sollya list from Sage integral arguments.
166 237 storres
    """
167 237 storres
    if len(args) == 0:
168 237 storres
        return pobyso_build_list_so_so()
169 237 storres
    ## Make a Sage list of integral constants.
170 237 storres
    intsList = []
171 237 storres
    for intElem in args:
172 237 storres
        intsList.append(pobyso_constant_from_int_sa_so(intElem))
173 237 storres
    return pobyso_build_list_so_so(*intsList)
174 237 storres
175 237 storres
def pobyso_build_list_so_so(*args):
176 237 storres
    """
177 237 storres
    Make a Sollya list out of Sollya objects passed as arguments.
178 237 storres
    If one wants to call it with a list argument, should prepend a "*"
179 237 storres
    before the list variable name.
180 237 storres
    Elements of the list are "eaten" (should not be cleared individually,
181 237 storres
    are cleared when the list is cleared).
182 237 storres
    """
183 237 storres
    if len(args) == 0:
184 237 storres
        ## Called with an empty list produced "error".
185 237 storres
        return sollya_lib_build_list(None)
186 237 storres
    index = 0
187 237 storres
    ## Append the Sollya elements one by one.
188 237 storres
    for elementSo in args:
189 237 storres
        if index == 0:
190 237 storres
            listSo = sollya_lib_build_list(elementSo, None)
191 237 storres
        else:
192 237 storres
            listSo = sollya_lib_append(listSo, elementSo)
193 237 storres
        index += 1
194 237 storres
    return listSo
195 237 storres
# End pobyso_build list_so_so
196 237 storres
197 237 storres
198 85 storres
def pobyso_change_var_in_function_so_so(funcSo, chvarExpSo):
199 54 storres
    """
200 85 storres
    Variable change in a function.
201 85 storres
    """
202 85 storres
    return(sollya_lib_evaluate(funcSo,chvarExpSo))
203 85 storres
# End pobyso_change_var_in_function_so_so
204 85 storres
205 85 storres
def pobyso_chebyshevform_so_so(functionSo, degreeSo, intervalSo):
206 85 storres
    resultSo = sollya_lib_chebyshevform(functionSo, degreeSo, intervalSo)
207 85 storres
    return(resultSo)
208 85 storres
# End pobyso_chebyshevform_so_so.
209 85 storres
210 237 storres
def pobyso_clear_obj(objSo):
211 237 storres
    """
212 237 storres
    Free a Sollya object's memory.
213 237 storres
    Very thin wrapper around sollya_lib_clear_obj().
214 237 storres
    """
215 237 storres
    sollya_lib_clear_obj(objSo)
216 237 storres
# End pobyso_clear_obj.
217 237 storres
218 117 storres
def pobyso_clear_taylorform_sa_so(taylorFormSaSo):
219 117 storres
    """
220 117 storres
    This method is necessary to correctly clean up the memory from Taylor forms.
221 117 storres
    These are made of a Sollya object, a Sollya object list, a Sollya object.
222 117 storres
    For no clearly understood reason, sollya_lib_clear_object_list crashed
223 117 storres
    when applied to the object list.
224 117 storres
    Here, we decompose it into Sage list of Sollya objects references and we
225 237 storres
     clear them one at a time.
226 117 storres
    """
227 117 storres
    sollya_lib_clear_obj(taylorFormSaSo[0])
228 117 storres
    (coefficientsErrorsListSaSo, numElementsSa, isEndEllipticSa) = \
229 117 storres
        pobyso_get_list_elements_so_so(taylorFormSaSo[1])
230 117 storres
    for element in coefficientsErrorsListSaSo:
231 117 storres
        sollya_lib_clear_obj(element)
232 117 storres
    sollya_lib_clear_obj(taylorFormSaSo[1])
233 117 storres
    sollya_lib_clear_obj(taylorFormSaSo[2])
234 117 storres
# End pobyso_clear_taylorform_sa_so
235 117 storres
236 85 storres
def pobyso_cmp(rnArgSa, cteSo):
237 85 storres
    """
238 237 storres
    Deprecated, use pobyso_cmp_sa_so_sa instead.
239 237 storres
    """
240 237 storres
    print "Deprecated, use pobyso_cmp_sa_so_sa instead."
241 237 storres
    return pobyso_cmp_sa_so_sa(rnArgSa, cteSo)
242 237 storres
# End  pobyso_cmp
243 237 storres
244 237 storres
def pobyso_cmp_sa_so_sa(rnArgSa, cteSo):
245 237 storres
    """
246 54 storres
    Compare the MPFR value a RealNumber with that of a Sollya constant.
247 54 storres

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

811 118 storres
    INPUT:
812 118 storres
    - objectListSo: a Sollya list of Sollya objects.
813 118 storres

814 118 storres
    OUTPUT:
815 118 storres
    - a Sage/Python tuple made of:
816 118 storres
      - a Sage/Python list of Sollya objects,
817 118 storres
      - a Sage/Python int holding the number of elements,
818 118 storres
      - a Sage/Python int stating (!= 0) that the list is end-elliptic.
819 118 storres
    NOTE::
820 118 storres
        We recover the addresses of the Sollya object from the list of pointers
821 118 storres
        returned by sollya_lib_get_list_elements. The list itself is freed.
822 118 storres
    TODO::
823 118 storres
        Figure out what to do with numElements since the number of elements
824 118 storres
        can easily be recovered from the list itself.
825 118 storres
        Ditto for isEndElliptic.
826 51 storres
    """
827 5 storres
    listAddress = POINTER(c_longlong)()
828 5 storres
    numElements = c_int(0)
829 5 storres
    isEndElliptic = c_int(0)
830 117 storres
    listAsSageList = []
831 5 storres
    result = sollya_lib_get_list_elements(byref(listAddress),\
832 54 storres
                                          byref(numElements),\
833 54 storres
                                          byref(isEndElliptic),\
834 117 storres
                                          objectListSo)
835 5 storres
    if result == 0 :
836 5 storres
        return None
837 5 storres
    for i in xrange(0, numElements.value, 1):
838 118 storres
       #listAsSageList.append(sollya_lib_copy_obj(listAddress[i]))
839 118 storres
       listAsSageList.append(listAddress[i])
840 117 storres
       # Clear each of the elements returned by Sollya.
841 118 storres
       #sollya_lib_clear_obj(listAddress[i])
842 117 storres
    # Free the list itself.
843 117 storres
    sollya_lib_free(listAddress)
844 209 storres
    return (listAsSageList, numElements.value, isEndElliptic.value)
845 5 storres
846 38 storres
def pobyso_get_max_prec_of_exp(soExp):
847 38 storres
    """ Legacy function. See pobyso_get_max_prec_of_exp_so_sa. """
848 209 storres
    return pobyso_get_max_prec_of_exp_so_sa(soExp)
849 5 storres
850 85 storres
def pobyso_get_max_prec_of_exp_so_sa(expSo):
851 38 storres
    """
852 38 storres
    Get the maximum precision used for the numbers in a Sollya expression.
853 52 storres

854 52 storres
    Arguments:
855 52 storres
    soExp -- a Sollya expression pointer
856 52 storres
    Return value:
857 52 storres
    A Python integer
858 38 storres
    TODO:
859 38 storres
    - error management;
860 38 storres
    - correctly deal with numerical type such as DOUBLEEXTENDED.
861 38 storres
    """
862 226 storres
    if expSo is None:
863 226 storres
        print inspect.stack()[0][3], ": expSo is None."
864 226 storres
        return 0
865 5 storres
    maxPrecision = 0
866 52 storres
    minConstPrec = 0
867 52 storres
    currentConstPrec = 0
868 226 storres
    #pobyso_autoprint(expSo)
869 85 storres
    operator = pobyso_get_head_function_so_sa(expSo)
870 5 storres
    if (operator != SOLLYA_BASE_FUNC_CONSTANT) and \
871 5 storres
    (operator != SOLLYA_BASE_FUNC_FREE_VARIABLE):
872 85 storres
        (arity, subexpressions) = pobyso_get_subfunctions_so_sa(expSo)
873 5 storres
        for i in xrange(arity):
874 5 storres
            maxPrecisionCandidate = \
875 38 storres
                pobyso_get_max_prec_of_exp_so_sa(subexpressions[i])
876 5 storres
            if maxPrecisionCandidate > maxPrecision:
877 5 storres
                maxPrecision = maxPrecisionCandidate
878 209 storres
        return maxPrecision
879 5 storres
    elif operator == SOLLYA_BASE_FUNC_CONSTANT:
880 85 storres
        #minConstPrec = pobyso_get_min_prec_of_constant_so_sa(expSo)
881 52 storres
        #currentConstPrec = pobyso_get_min_prec_of_constant_so_sa(soExp)
882 52 storres
        #print minConstPrec, " - ", currentConstPrec
883 209 storres
        return pobyso_get_min_prec_of_constant_so_sa(expSo)
884 52 storres
885 5 storres
    elif operator == SOLLYA_BASE_FUNC_FREE_VARIABLE:
886 209 storres
        return 0
887 5 storres
    else:
888 38 storres
        print "pobyso_get_max_prec_of_exp_so_sa: unexepected operator."
889 209 storres
        return 0
890 5 storres
891 85 storres
def pobyso_get_min_prec_of_constant_so_sa(constExpSo):
892 52 storres
    """
893 52 storres
    Get the minimum precision necessary to represent the value of a Sollya
894 52 storres
    constant.
895 52 storres
    MPFR_MIN_PREC and powers of 2 are taken into account.
896 209 storres
    We assume that constExpSo is a pointer to a Sollay constant expression.
897 52 storres
    """
898 85 storres
    constExpAsRnSa = pobyso_get_constant_as_rn_so_sa(constExpSo)
899 85 storres
    return(min_mpfr_size(get_rn_value(constExpAsRnSa)))
900 52 storres
901 200 storres
def pobyso_get_poly_so_sa(polySo, realFieldSa=None):
902 200 storres
    """
903 200 storres
    Convert a Sollya polynomial into a Sage polynomial.
904 209 storres
    Legacy function. Use pobyso_float_poly_so_sa() instead.
905 200 storres
    """
906 213 storres
    return pobyso_float_poly_so_sa(polySo,realFieldSa)
907 200 storres
# End pobyso_get_poly_so_sa
908 200 storres
909 200 storres
def pobyso_get_prec():
910 200 storres
    """ Legacy function. See pobyso_get_prec_so_sa(). """
911 209 storres
    return pobyso_get_prec_so_sa()
912 200 storres
913 200 storres
def pobyso_get_prec_so():
914 200 storres
    """
915 200 storres
    Get the current default precision in Sollya.
916 200 storres
    The return value is a Sollya object.
917 200 storres
    Usefull when modifying the precision back and forth by avoiding
918 200 storres
    extra conversions.
919 200 storres
    """
920 209 storres
    return sollya_lib_get_prec(None)
921 200 storres
922 200 storres
def pobyso_get_prec_so_sa():
923 200 storres
    """
924 200 storres
    Get the current default precision in Sollya.
925 200 storres
    The return value is Sage/Python int.
926 200 storres
    """
927 227 storres
    precSo = sollya_lib_get_prec()
928 227 storres
    precSa = pobyso_constant_from_int_so_sa(precSo)
929 200 storres
    sollya_lib_clear_obj(precSo)
930 227 storres
    return precSa
931 200 storres
# End pobyso_get_prec_so_sa.
932 200 storres
933 209 storres
def pobyso_get_prec_so_so_sa():
934 209 storres
    """
935 209 storres
    Return the current precision both as a Sollya object and a
936 209 storres
    Sage integer as hybrid tuple.
937 209 storres
    To avoid multiple calls for precision manipulations.
938 209 storres
    """
939 227 storres
    precSo = sollya_lib_get_prec()
940 227 storres
    precSa = pobyso_constant_from_int_so_sa(precSo)
941 227 storres
    return (precSo, int(precSa))
942 200 storres
943 200 storres
def pobyso_get_prec_of_constant(ctExpSo):
944 200 storres
    """ Legacy function. See pobyso_get_prec_of_constant_so_sa. """
945 209 storres
    return pobyso_get_prec_of_constant_so_sa(ctExpSo)
946 200 storres
947 200 storres
def pobyso_get_prec_of_constant_so_sa(ctExpSo):
948 200 storres
    """
949 200 storres
    Tries to find a precision to represent ctExpSo without rounding.
950 200 storres
    If not possible, returns None.
951 200 storres
    """
952 218 storres
    #print "Entering pobyso_get_prec_of_constant_so_sa..."
953 200 storres
    prec = c_int(0)
954 200 storres
    retc = sollya_lib_get_prec_of_constant(byref(prec), ctExpSo, None)
955 200 storres
    if retc == 0:
956 218 storres
        #print "pobyso_get_prec_of_constant_so_sa failed."
957 209 storres
        return None
958 218 storres
    #print "...exiting pobyso_get_prec_of_constant_so_sa."
959 209 storres
    return int(prec.value)
960 200 storres
961 200 storres
def pobyso_get_prec_of_range_so_sa(rangeSo):
962 200 storres
    """
963 200 storres
    Returns the number of bits elements of a range are coded with.
964 200 storres
    """
965 200 storres
    prec = c_int(0)
966 200 storres
    retc = sollya_lib_get_prec_of_range(byref(prec), rangeSo, None)
967 200 storres
    if retc == 0:
968 200 storres
        return(None)
969 209 storres
    return int(prec.value)
970 200 storres
# End pobyso_get_prec_of_range_so_sa()
971 200 storres
972 85 storres
def pobyso_get_sage_exp_from_sollya_exp(sollyaExpSo, realField = RR):
973 38 storres
    """ Legacy function. See pobyso_get_sage_exp_from_sollya_exp_so_sa. """
974 209 storres
    return pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo,
975 209 storres
                                                     realField = RR)
976 38 storres
977 85 storres
def pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, realFieldSa = RR):
978 5 storres
    """
979 38 storres
    Get a Sage expression from a Sollya expression.
980 38 storres
    Currently only tested with polynomials with floating-point coefficients.
981 5 storres
    Notice that, in the returned polynomial, the exponents are RealNumbers.
982 5 storres
    """
983 5 storres
    #pobyso_autoprint(sollyaExp)
984 85 storres
    operatorSa = pobyso_get_head_function_so_sa(sollyaExpSo)
985 83 storres
    sollyaLibFreeVariableName = sollya_lib_get_free_variable_name()
986 213 storres
    ## Get rid of the "_"'s in "_x_", if any.
987 213 storres
    sollyaLibFreeVariableName = re.sub('_', '', sollyaLibFreeVariableName)
988 5 storres
    # Constants and the free variable are special cases.
989 5 storres
    # All other operator are dealt with in the same way.
990 85 storres
    if (operatorSa != SOLLYA_BASE_FUNC_CONSTANT) and \
991 85 storres
       (operatorSa != SOLLYA_BASE_FUNC_FREE_VARIABLE):
992 85 storres
        (aritySa, subexpressionsSa) = pobyso_get_subfunctions_so_sa(sollyaExpSo)
993 85 storres
        if aritySa == 1:
994 85 storres
            sageExpSa = eval(pobyso_function_type_as_string_so_sa(operatorSa) + \
995 85 storres
            "(" + pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], \
996 85 storres
            realFieldSa) + ")")
997 85 storres
        elif aritySa == 2:
998 63 storres
            # We do not get through the preprocessor.
999 63 storres
            # The "^" operator is then a special case.
1000 85 storres
            if operatorSa == SOLLYA_BASE_FUNC_POW:
1001 85 storres
                operatorAsStringSa = "**"
1002 5 storres
            else:
1003 85 storres
                operatorAsStringSa = \
1004 85 storres
                    pobyso_function_type_as_string_so_sa(operatorSa)
1005 85 storres
            sageExpSa = \
1006 85 storres
              eval("pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], realFieldSa)"\
1007 85 storres
              + " " + operatorAsStringSa + " " + \
1008 85 storres
                   "pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[1], realFieldSa)")
1009 63 storres
        # We do not know yet how to deal with arity >= 3
1010 63 storres
        # (is there any in Sollya anyway?).
1011 5 storres
        else:
1012 85 storres
            sageExpSa = eval('None')
1013 209 storres
        return sageExpSa
1014 85 storres
    elif operatorSa == SOLLYA_BASE_FUNC_CONSTANT:
1015 5 storres
        #print "This is a constant"
1016 85 storres
        return pobyso_get_constant_as_rn_with_rf_so_sa(sollyaExpSo, realFieldSa)
1017 85 storres
    elif operatorSa == SOLLYA_BASE_FUNC_FREE_VARIABLE:
1018 218 storres
        #print "This is the free variable"
1019 209 storres
        return eval(sollyaLibFreeVariableName)
1020 5 storres
    else:
1021 5 storres
        print "Unexpected"
1022 5 storres
        return eval('None')
1023 185 storres
# End pobyso_get_sage_exp_from_sollya_exp_so_sa
1024 73 storres
1025 185 storres
1026 38 storres
def pobyso_get_subfunctions(expressionSo):
1027 38 storres
    """ Legacy function. See pobyso_get_subfunctions_so_sa. """
1028 209 storres
    return pobyso_get_subfunctions_so_sa(expressionSo)
1029 200 storres
# End pobyso_get_subfunctions.
1030 200 storres
1031 38 storres
def pobyso_get_subfunctions_so_sa(expressionSo):
1032 38 storres
    """
1033 38 storres
    Get the subfunctions of an expression.
1034 38 storres
    Return the number of subfunctions and the list of subfunctions addresses.
1035 55 storres
    S.T.: Could not figure out another way than that ugly list of declarations
1036 83 storres
    to recover the addresses of the subfunctions.
1037 83 storres
    We limit ourselves to arity 8 functions.
1038 38 storres
    """
1039 5 storres
    subf0 = c_int(0)
1040 5 storres
    subf1 = c_int(0)
1041 5 storres
    subf2 = c_int(0)
1042 5 storres
    subf3 = c_int(0)
1043 5 storres
    subf4 = c_int(0)
1044 5 storres
    subf5 = c_int(0)
1045 5 storres
    subf6 = c_int(0)
1046 5 storres
    subf7 = c_int(0)
1047 5 storres
    subf8 = c_int(0)
1048 5 storres
    arity = c_int(0)
1049 5 storres
    nullPtr = POINTER(c_int)()
1050 38 storres
    sollya_lib_get_subfunctions(expressionSo, byref(arity), \
1051 83 storres
      byref(subf0), byref(subf1), byref(subf2), byref(subf3), \
1052 83 storres
      byref(subf4), byref(subf5),\
1053 83 storres
      byref(subf6), byref(subf7), byref(subf8), nullPtr, None)
1054 83 storres
#    byref(cast(subfunctions[0], POINTER(c_int))), \
1055 83 storres
#    byref(cast(subfunctions[0], POINTER(c_int))), \
1056 83 storres
#    byref(cast(subfunctions[2], POINTER(c_int))), \
1057 83 storres
#    byref(cast(subfunctions[3], POINTER(c_int))), \
1058 83 storres
#    byref(cast(subfunctions[4], POINTER(c_int))), \
1059 83 storres
#    byref(cast(subfunctions[5], POINTER(c_int))), \
1060 83 storres
#    byref(cast(subfunctions[6], POINTER(c_int))), \
1061 83 storres
#    byref(cast(subfunctions[7], POINTER(c_int))), \
1062 5 storres
#    byref(cast(subfunctions[8], POINTER(c_int))), nullPtr)
1063 83 storres
    subfunctions = [subf0, subf1, subf2, subf3, subf4, subf5, subf6, subf7, \
1064 83 storres
                    subf8]
1065 5 storres
    subs = []
1066 5 storres
    if arity.value > pobyso_max_arity:
1067 38 storres
        return(0,[])
1068 5 storres
    for i in xrange(arity.value):
1069 5 storres
        subs.append(int(subfunctions[i].value))
1070 5 storres
        #print subs[i]
1071 209 storres
    return (int(arity.value), subs)
1072 200 storres
# End pobyso_get_subfunctions_so_sa
1073 5 storres
1074 155 storres
def pobyso_guess_degree_sa_sa(functionSa, intervalSa, approxErrorSa,
1075 155 storres
                              weightSa=None, degreeBoundSa=None):
1076 155 storres
    """
1077 155 storres
    Sa_sa variant of the solly_guessdegree function.
1078 155 storres
    Return 0 if something goes wrong.
1079 155 storres
    """
1080 159 storres
    functionAsStringSa = functionSa._assume_str().replace('_SAGE_VAR_', '')
1081 154 storres
    functionSo = pobyso_parse_string_sa_so(functionAsStringSa)
1082 155 storres
    if pobyso_is_error_so_sa(functionSo):
1083 155 storres
        sollya_lib_clear_obj(functionSo)
1084 155 storres
        return 0
1085 154 storres
    rangeSo = pobyso_interval_to_range_sa_so(intervalSa)
1086 155 storres
    # The approximation error is expected to be a floating point number.
1087 155 storres
    if pobyso_is_floating_point_number_sa_sa(approxErrorSa):
1088 155 storres
        approxErrorSo = pobyso_constant_sa_so(approxErrorSa)
1089 155 storres
    else:
1090 155 storres
        approxErrorSo = pobyso_constant_sa_so(RR(approxErrorSa))
1091 154 storres
    if not weightSa is None:
1092 159 storres
        weightAsStringSa = weightSa._assume_str().replace('_SAGE_VAR_', '')
1093 154 storres
        weightSo = pobyso_parse_string_sa_so(weightAsStringSa)
1094 166 storres
        if pobyso_is_error_so_sa(weightSo):
1095 155 storres
            sollya_lib_clear_obj(functionSo)
1096 155 storres
            sollya_lib_clear_obj(rangeSo)
1097 155 storres
            sollya_lib_clear_obj(approxErrorSo)
1098 155 storres
            sollya_lib_clear_obj(weightSo)
1099 155 storres
            return 0
1100 154 storres
    else:
1101 154 storres
        weightSo = None
1102 154 storres
    if not degreeBoundSa is None:
1103 154 storres
        degreeBoundSo = pobyso_constant_from_int_sa_so(degreeBoundSa)
1104 154 storres
    else:
1105 154 storres
        degreeBoundSo = None
1106 154 storres
    guessedDegreeSa = pobyso_guess_degree_so_sa(functionSo,
1107 162 storres
                                                rangeSo,
1108 162 storres
                                                approxErrorSo,
1109 162 storres
                                                weightSo,
1110 162 storres
                                                degreeBoundSo)
1111 154 storres
    sollya_lib_clear_obj(functionSo)
1112 154 storres
    sollya_lib_clear_obj(rangeSo)
1113 155 storres
    sollya_lib_clear_obj(approxErrorSo)
1114 154 storres
    if not weightSo is None:
1115 154 storres
        sollya_lib_clear_obj(weightSo)
1116 154 storres
    if not degreeBoundSo is None:
1117 154 storres
        sollya_lib_clear_obj(degreeBoundSo)
1118 154 storres
    return guessedDegreeSa
1119 154 storres
# End poyso_guess_degree_sa_sa
1120 154 storres
1121 153 storres
def pobyso_guess_degree_so_sa(functionSo, rangeSo, errorSo, weightSo=None, \
1122 154 storres
                              degreeBoundSo=None):
1123 154 storres
    """
1124 154 storres
    Thin wrapper around the guessdegree function.
1125 154 storres
    Nevertheless, some precision control stuff has been appended.
1126 154 storres
    """
1127 154 storres
    # Deal with Sollya internal precision issues: if it is too small,
1128 154 storres
    # compared with the error, increases it to about twice -log2(error).
1129 154 storres
    errorSa = pobyso_get_constant_as_rn_with_rf_so_sa(errorSo)
1130 154 storres
    log2ErrorSa = errorSa.log2()
1131 154 storres
    if log2ErrorSa < 0:
1132 154 storres
        neededPrecisionSa = int(2 * int(-log2ErrorSa) / 64) * 64
1133 154 storres
    else:
1134 154 storres
        neededPrecisionSa = int(2 * int(log2ErrorSa) / 64) * 64
1135 154 storres
    #print "Needed precision:", neededPrecisionSa
1136 226 storres
    sollyaPrecisionHasChanged = False
1137 226 storres
    (initialPrecSo, initialPrecSa) = pobyso_get_prec_so_so_sa()
1138 226 storres
    if neededPrecisionSa > initialPrecSa:
1139 226 storres
        if neededPrecisionSa <= 2:
1140 226 storres
            print inspect.stack()[0][3], ": precision change <= 2 requested."
1141 154 storres
        pobyso_set_prec_sa_so(neededPrecisionSa)
1142 226 storres
        sollyaPrecisionHasChanged = True
1143 166 storres
    #print "Guessing degree..."
1144 153 storres
    # weightSo and degreeBoundsSo are optional arguments.
1145 162 storres
    # As declared, sollya_lib_guessdegree must take 5 arguments.
1146 153 storres
    if weightSo is None:
1147 162 storres
        degreeRangeSo = sollya_lib_guessdegree(functionSo, rangeSo, errorSo,
1148 162 storres
                                               0, 0, None)
1149 154 storres
    elif degreeBoundSo is None:
1150 153 storres
        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, \
1151 162 storres
                                                errorSo, weightSo, 0, None)
1152 153 storres
    else:
1153 153 storres
        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, errorSo, \
1154 154 storres
                                                weightSo, degreeBoundSo, None)
1155 166 storres
    #print "...degree guess done."
1156 154 storres
    # Restore internal precision, if applicable.
1157 226 storres
    if sollyaPrecisionHasChanged:
1158 226 storres
        pobyso_set_prec_so_so(initialPrecSo)
1159 226 storres
        sollya_lib_clear_obj(initialPrecSo)
1160 154 storres
    degreeIntervalSa = pobyso_range_to_interval_so_sa(degreeRangeSo)
1161 154 storres
    sollya_lib_clear_obj(degreeRangeSo)
1162 154 storres
    # When ok, both bounds match.
1163 154 storres
    # When the degree bound is too low, the upper bound is the degree
1164 154 storres
    # for which the error can be honored.
1165 154 storres
    # When it really goes wrong, the upper bound is infinity.
1166 154 storres
    if degreeIntervalSa.lower() == degreeIntervalSa.upper():
1167 154 storres
        return int(degreeIntervalSa.lower())
1168 154 storres
    else:
1169 154 storres
        if degreeIntervalSa.upper().is_infinity():
1170 154 storres
            return None
1171 154 storres
        else:
1172 154 storres
            return int(degreeIntervalSa.upper())
1173 154 storres
    # End pobyso_guess_degree_so_sa
1174 153 storres
1175 215 storres
def pobyso_inf_so_so(intervalSo):
1176 215 storres
    """
1177 215 storres
    Very thin wrapper around sollya_lib_inf().
1178 215 storres
    """
1179 215 storres
    return sollya_lib_inf(intervalSo)
1180 215 storres
# End pobyso_inf_so_so.
1181 215 storres
1182 53 storres
def pobyso_infnorm_so_so(func, interval, file = None, intervalList = None):
1183 54 storres
    print "Do not use this function. User pobyso_supnorm_so_so instead."
1184 209 storres
    return None
1185 53 storres
1186 84 storres
def pobyso_interval_to_range_sa_so(intervalSa, precisionSa=None):
1187 84 storres
    if precisionSa is None:
1188 84 storres
        precisionSa = intervalSa.parent().precision()
1189 84 storres
    intervalSo = pobyso_bounds_to_range_sa_so(intervalSa.lower(),\
1190 84 storres
                                              intervalSa.upper(),\
1191 84 storres
                                              precisionSa)
1192 209 storres
    return intervalSo
1193 84 storres
# End pobyso_interval_to_range_sa_so
1194 84 storres
1195 155 storres
def pobyso_is_error_so_sa(objSo):
1196 155 storres
    """
1197 155 storres
    Thin wrapper around the sollya_lib_obj_is_error() function.
1198 155 storres
    """
1199 155 storres
    if sollya_lib_obj_is_error(objSo) != 0:
1200 155 storres
        return True
1201 155 storres
    else:
1202 155 storres
        return False
1203 155 storres
# End pobyso_is_error-so_sa
1204 155 storres
1205 155 storres
def pobyso_is_floating_point_number_sa_sa(numberSa):
1206 155 storres
    """
1207 209 storres
    Check whether a Sage number is floating point.
1208 209 storres
    Exception stuff added because numbers other than
1209 209 storres
    floating-point ones do not have the is_real() attribute.
1210 155 storres
    """
1211 209 storres
    try:
1212 209 storres
        return numberSa.is_real()
1213 209 storres
    except AttributeError:
1214 209 storres
        return False
1215 209 storres
# End pobyso_is_floating_piont_number_sa_sa
1216 155 storres
1217 226 storres
def pobyso_is_range_so_sa(rangeCandidateSo):
1218 226 storres
    """
1219 226 storres
    Thin wrapper over sollya_lib_is_range.
1220 226 storres
    """
1221 226 storres
    return sollya_lib_obj_is_range(rangeCandidateSo) != 0
1222 226 storres
1223 226 storres
# End pobyso_is_range_so_sa
1224 226 storres
1225 226 storres
1226 37 storres
def pobyso_lib_init():
1227 37 storres
    sollya_lib_init(None)
1228 116 storres
1229 116 storres
def pobyso_lib_close():
1230 116 storres
    sollya_lib_close(None)
1231 237 storres
1232 85 storres
def pobyso_name_free_variable(freeVariableNameSa):
1233 38 storres
    """ Legacy function. See pobyso_name_free_variable_sa_so. """
1234 85 storres
    pobyso_name_free_variable_sa_so(freeVariableNameSa)
1235 38 storres
1236 85 storres
def pobyso_name_free_variable_sa_so(freeVariableNameSa):
1237 83 storres
    """
1238 83 storres
    Set the free variable name in Sollya from a Sage string.
1239 83 storres
    """
1240 85 storres
    sollya_lib_name_free_variable(freeVariableNameSa)
1241 37 storres
1242 5 storres
def pobyso_parse_string(string):
1243 38 storres
    """ Legacy function. See pobyso_parse_string_sa_so. """
1244 209 storres
    return pobyso_parse_string_sa_so(string)
1245 38 storres
1246 38 storres
def pobyso_parse_string_sa_so(string):
1247 83 storres
    """
1248 155 storres
    Get the Sollya expression computed from a Sage string or
1249 155 storres
    a Sollya error object if parsing failed.
1250 83 storres
    """
1251 209 storres
    return sollya_lib_parse_string(string)
1252 5 storres
1253 200 storres
def pobyso_precision_so_sa(ctExpSo):
1254 209 storres
    """
1255 209 storres
    Computes the necessary precision to represent a number.
1256 209 storres
    If x is not zero, it can be uniquely written as x = m · 2e
1257 209 storres
    where m is an odd integer and e is an integer.
1258 209 storres
    precision(x) returns the number of bits necessary to write m
1259 209 storres
    in binary (i.e. ceil(log2(m))).
1260 209 storres
    """
1261 209 storres
    #TODO: take care of the special case: 0, @NaN@, @Inf@
1262 200 storres
    precisionSo = sollya_lib_precision(ctExpSo)
1263 200 storres
    precisionSa = pobyso_constant_from_int_so_sa(precisionSo)
1264 200 storres
    sollya_lib_clear_obj(precisionSo)
1265 200 storres
    return precisionSa
1266 200 storres
# End pobyso_precision_so_sa
1267 215 storres
1268 217 storres
def pobyso_polynomial_coefficients_progressive_round_so_so(polySo,
1269 217 storres
                                                           funcSo,
1270 217 storres
                                                           icSo,
1271 217 storres
                                                           intervalSo,
1272 217 storres
                                                           itpSo,
1273 217 storres
                                                           ftpSo,
1274 217 storres
                                                           maxPrecSo,
1275 219 storres
                                                           maxErrSo,
1276 219 storres
                                                           debug=False):
1277 219 storres
    if debug:
1278 219 storres
        print "Input arguments:"
1279 219 storres
        pobyso_autoprint(polySo)
1280 219 storres
        pobyso_autoprint(funcSo)
1281 219 storres
        pobyso_autoprint(icSo)
1282 219 storres
        pobyso_autoprint(intervalSo)
1283 219 storres
        pobyso_autoprint(itpSo)
1284 219 storres
        pobyso_autoprint(ftpSo)
1285 219 storres
        pobyso_autoprint(maxPrecSo)
1286 219 storres
        pobyso_autoprint(maxErrSo)
1287 219 storres
        print "________________"
1288 200 storres
1289 217 storres
    ## Higher order function see:
1290 217 storres
    #  http://effbot.org/pyfaq/how-do-you-make-a-higher-order-function-in-python.htm
1291 217 storres
    def precision_decay_ratio_function(degreeSa):
1292 217 storres
        def outer(x):
1293 217 storres
            def inner(x):
1294 217 storres
                we = 3/8
1295 217 storres
                wq = 2/8
1296 217 storres
                a  = 2.2
1297 217 storres
                b  = 2
1298 217 storres
                return we*(exp(x/a)-1) +  wq*((b*x)**2) + (1-we-wq)*x
1299 217 storres
            return  inner(x)/inner(degreeSa)
1300 217 storres
        return outer
1301 217 storres
1302 217 storres
    #
1303 217 storres
    degreeSa        = pobyso_polynomial_degree_so_sa(polySo)
1304 217 storres
    ratio           = precision_decay_ratio_function(degreeSa)
1305 217 storres
    itpSa           = pobyso_constant_from_int_so_sa(itpSo)
1306 217 storres
    ftpSa           = pobyso_constant_from_int_so_sa(ftpSo)
1307 217 storres
    maxPrecSa       = pobyso_constant_from_int_so_sa(maxPrecSo)
1308 217 storres
    maxErrSa        = pobyso_get_constant_as_rn_so_sa(maxErrSo)
1309 219 storres
    if debug:
1310 219 storres
        print "degreeSa:", degreeSa
1311 219 storres
        print "ratio:", ratio
1312 219 storres
        print "itpsSa:", itpSa
1313 219 storres
        print "ftpSa:", ftpSa
1314 219 storres
        print "maxPrecSa:", maxPrecSa
1315 219 storres
        print "maxErrSa:", maxErrSa
1316 217 storres
    lastResPolySo   = None
1317 218 storres
    lastInfNormSo   = None
1318 219 storres
    #print "About to enter the while loop..."
1319 217 storres
    while True:
1320 218 storres
        resPolySo   = pobyso_constant_0_sa_so()
1321 217 storres
        pDeltaSa    = ftpSa - itpSa
1322 217 storres
        for indexSa in reversed(xrange(0,degreeSa+1)):
1323 218 storres
            #print "Index:", indexSa
1324 217 storres
            indexSo = pobyso_constant_from_int_sa_so(indexSa)
1325 217 storres
            #pobyso_autoprint(indexSo)
1326 217 storres
            #print ratio(indexSa)
1327 217 storres
            ctpSa = floor(ftpSa - (pDeltaSa * ratio(indexSa)))
1328 217 storres
            ctpSo = pobyso_constant_from_int_sa_so(ctpSa)
1329 219 storres
            if debug:
1330 219 storres
                print "Index:", indexSa, " - Target precision:",
1331 219 storres
                pobyso_autoprint(ctpSo)
1332 217 storres
            cmonSo  = \
1333 217 storres
                sollya_lib_build_function_mul(sollya_lib_coeff(polySo, indexSo),
1334 217 storres
                                      sollya_lib_build_function_pow( \
1335 217 storres
                                          sollya_lib_build_function_free_variable(), \
1336 217 storres
                                          indexSo))
1337 217 storres
            #pobyso_autoprint(cmonSo)
1338 217 storres
            cmonrSo = pobyso_round_coefficients_single_so_so(cmonSo, ctpSo)
1339 217 storres
            sollya_lib_clear_obj(cmonSo)
1340 217 storres
            #pobyso_autoprint(cmonrSo)
1341 217 storres
            resPolySo = sollya_lib_build_function_add(resPolySo,
1342 217 storres
                                                      cmonrSo)
1343 218 storres
            #pobyso_autoprint(resPolySo)
1344 217 storres
        # End for index
1345 217 storres
        freeVarSo     = sollya_lib_build_function_free_variable()
1346 217 storres
        changeVarSo   = sollya_lib_sub(freeVarSo, icSo)
1347 217 storres
        resPolyCvSo   = sollya_lib_evaluate(resPolySo, changeVarSo)
1348 218 storres
        errFuncSo = sollya_lib_build_function_sub(sollya_lib_copy_obj(funcSo),
1349 218 storres
                                                  resPolyCvSo)
1350 218 storres
        infNormSo = sollya_lib_dirtyinfnorm(errFuncSo, intervalSo)
1351 217 storres
        cerrSa    = pobyso_get_constant_as_rn_so_sa(infNormSo)
1352 219 storres
        if debug:
1353 219 storres
            print "Infnorm (Sollya):",
1354 219 storres
            pobyso_autoprint(infNormSo)
1355 218 storres
        sollya_lib_clear_obj(errFuncSo)
1356 218 storres
        #print "Infnorm  (Sage):", cerrSa
1357 217 storres
        if (cerrSa > maxErrSa):
1358 219 storres
            if debug:
1359 219 storres
                print "Error is too large."
1360 218 storres
            if lastResPolySo is None:
1361 219 storres
                if debug:
1362 219 storres
                    print "Enlarging prec."
1363 217 storres
                ntpSa = floor(ftpSa + ftpSa/50)
1364 217 storres
                ## Can't enlarge (numerical)
1365 217 storres
                if ntpSa == ftpSa:
1366 217 storres
                    sollya_lib_clear_obj(resPolySo)
1367 217 storres
                    return None
1368 217 storres
                ## Can't enlarge (not enough precision left)
1369 217 storres
                if ntpSa > maxPrecSa:
1370 217 storres
                    sollya_lib_clear_obj(resPolySo)
1371 217 storres
                    return None
1372 217 storres
                ftpSa = ntpSa
1373 217 storres
                continue
1374 217 storres
            ## One enlargement took place.
1375 217 storres
            else:
1376 219 storres
                if debug:
1377 219 storres
                    print "Exit with the last before last polynomial."
1378 222 storres
                    print "Precision of highest degree monomial:", itpSa
1379 222 storres
                    print "Precision of constant term          :", ftpSa
1380 217 storres
                sollya_lib_clear_obj(resPolySo)
1381 218 storres
                sollya_lib_clear_obj(infNormSo)
1382 218 storres
                return (lastResPolySo, lastInfNormSo)
1383 218 storres
        # cerrSa <= maxErrSa: scrap more bits, possibly.
1384 217 storres
        else:
1385 219 storres
            if debug:
1386 219 storres
                print "Error is too small"
1387 219 storres
            if cerrSa <= (maxErrSa/2):
1388 219 storres
                if debug:
1389 219 storres
                    print "Shrinking prec."
1390 218 storres
                ntpSa = floor(ftpSa - ftpSa/50)
1391 218 storres
                ## Can't shrink (numerical)
1392 218 storres
                if ntpSa == ftpSa:
1393 218 storres
                    if not lastResPolySo is None:
1394 218 storres
                        sollya_lib_clear_obj(lastResPolySo)
1395 218 storres
                    if not lastInfNormSo is None:
1396 218 storres
                        sollya_lib_clear_obj(lastInfNormSo)
1397 222 storres
                    if debug:
1398 222 storres
                        print "Exit because can't shrink anymore (numerically)."
1399 222 storres
                        print "Precision of highest degree monomial:", itpSa
1400 222 storres
                        print "Precision of constant term          :", ftpSa
1401 218 storres
                    return (resPolySo, infNormSo)
1402 218 storres
                ## Can't shrink (not enough precision left)
1403 218 storres
                if ntpSa <= itpSa:
1404 218 storres
                    if not lastResPolySo is None:
1405 218 storres
                        sollya_lib_clear_obj(lastResPolySo)
1406 218 storres
                    if not lastInfNormSo is None:
1407 218 storres
                        sollya_lib_clear_obj(lastInfNormSo)
1408 222 storres
                        print "Exit because can't shrink anymore (no bits left)."
1409 222 storres
                        print "Precision of highest degree monomial:", itpSa
1410 222 storres
                        print "Precision of constant term          :", ftpSa
1411 218 storres
                    return (resPolySo, infNormSo)
1412 218 storres
                ftpSa = ntpSa
1413 217 storres
                if not lastResPolySo is None:
1414 217 storres
                    sollya_lib_clear_obj(lastResPolySo)
1415 218 storres
                if not lastInfNormSo is None:
1416 218 storres
                    sollya_lib_clear_obj(lastInfNormSo)
1417 218 storres
                lastResPolySo = resPolySo
1418 218 storres
                lastInfNormSo = infNormSo
1419 218 storres
                continue
1420 218 storres
            else: # Error is not that small, just return
1421 217 storres
                if not lastResPolySo is None:
1422 217 storres
                    sollya_lib_clear_obj(lastResPolySo)
1423 218 storres
                if not lastInfNormSo is None:
1424 218 storres
                    sollya_lib_clear_obj(lastInfNormSo)
1425 222 storres
                if debug:
1426 222 storres
                    print "Exit normally."
1427 222 storres
                    print "Precision of highest degree monomial:", itpSa
1428 222 storres
                    print "Precision of constant term          :", ftpSa
1429 218 storres
                return (resPolySo, infNormSo)
1430 217 storres
    # End wile True
1431 215 storres
    return None
1432 217 storres
# End pobyso_polynomial_coefficients_progressive_truncate_so_so.
1433 217 storres
1434 217 storres
def pobyso_polynomial_degree_so_sa(polySo):
1435 217 storres
    """
1436 217 storres
    Return the degree of a Sollya polynomial as a Sage int.
1437 217 storres
    """
1438 217 storres
    degreeSo = sollya_lib_degree(polySo)
1439 217 storres
    return pobyso_constant_from_int_so_sa(degreeSo)
1440 217 storres
# End pobyso_polynomial_degree_so_sa
1441 217 storres
1442 217 storres
def pobyso_polynomial_degree_so_so(polySo):
1443 217 storres
    """
1444 217 storres
    Thin wrapper around lib_sollya_degree().
1445 217 storres
    """
1446 217 storres
    return sollya_lib_degree(polySo)
1447 217 storres
# End pobyso_polynomial_degree_so_so
1448 217 storres
1449 5 storres
def pobyso_range(rnLowerBound, rnUpperBound):
1450 38 storres
    """ Legacy function. See pobyso_range_sa_so. """
1451 209 storres
    return pobyso_range_sa_so(rnLowerBound, rnUpperBound)
1452 38 storres
1453 226 storres
def pobyso_range_from_bounds_sa_so(rnLowerBound, rnUpperBound, precSa = None):
1454 226 storres
    """
1455 226 storres
    Create a Sollya range from 2 Sage real numbers as bounds
1456 226 storres
    """
1457 226 storres
    # TODO check precision stuff.
1458 226 storres
    sollyaPrecChanged = False
1459 226 storres
    (initialSollyaPrecSo, initialSollyaPrecSa) = \
1460 226 storres
        pobyso_get_prec_so_so_sa()
1461 226 storres
    if precSa is None:
1462 226 storres
        precSa = max(rnLowerBound.parent().prec(), rnUpperBound.parent().prec())
1463 226 storres
    if precSa > initialSollyaPrecSa:
1464 226 storres
        if precSa <= 2:
1465 226 storres
            print inspect.stack()[0][3], ": precision change <= 2 requested."
1466 226 storres
        pobyso_set_prec_sa_so(precSa)
1467 226 storres
        sollyaPrecChanged = True
1468 226 storres
    rangeSo = sollya_lib_range_from_bounds(get_rn_value(rnLowerBound),
1469 226 storres
                                           get_rn_value(rnUpperBound))
1470 226 storres
    if sollyaPrecChanged:
1471 226 storres
        pobyso_set_prec_so_so(initialSollyaPrecSo)
1472 234 storres
        sollya_lib_clear_obj(initialSollyaPrecSo)
1473 226 storres
    return rangeSo
1474 226 storres
# End pobyso_range_from_bounds_sa_so
1475 5 storres
1476 261 storres
def pobyso_range_list_so_sa(listSo):
1477 261 storres
    """
1478 261 storres
    Return a Sollya list of ranges as a Sage list of
1479 261 storres
    floating-point intervals.
1480 261 storres
    """
1481 261 storres
    listSa   = []
1482 261 storres
    ## The function returns none if the list is empty or an error has happened.
1483 261 storres
    retVal = pobyso_get_list_elements_so_so(listSo)
1484 261 storres
    if retVal is None:
1485 261 storres
        return listSa
1486 261 storres
    ## Just in case the interface is changed and an empty list is returned
1487 261 storres
    #  instead of None.
1488 261 storres
    elif len(retVal) == 0:
1489 261 storres
        return listSa
1490 261 storres
    else:
1491 261 storres
        ## Remember pobyso_get_list_elements_so_so returns more information
1492 261 storres
        #  than just the elements of the list (# elements, is_elliptic)
1493 261 storres
        listSaSo, numElements, isEndElliptic = retVal
1494 261 storres
    ## Return an empty list.
1495 261 storres
    if numElements == 0:
1496 261 storres
        return listSa
1497 261 storres
    ## Search first for the maximum precision of the elements
1498 261 storres
    maxPrecSa = 0
1499 261 storres
    for rangeSo in listSaSo:
1500 261 storres
        #pobyso_autoprint(floatSo)
1501 261 storres
        curPrecSa =  pobyso_get_prec_of_range_so_sa(rangeSo)
1502 261 storres
        if curPrecSa > maxPrecSa:
1503 261 storres
            maxPrecSa = curPrecSa
1504 261 storres
    ##
1505 261 storres
    intervalField = RealIntervalField(maxPrecSa)
1506 261 storres
    ##
1507 261 storres
    for rangeSo in listSaSo:
1508 261 storres
        listSa.append(pobyso_range_to_interval_so_sa(rangeSo, intervalField))
1509 261 storres
    return listSa
1510 261 storres
# End pobyso_range_list_so_sa
1511 261 storres
1512 226 storres
def pobyso_range_max_abs_so_so(rangeSo):
1513 226 storres
    """
1514 226 storres
    Return, as Sollya constant, the maximum absolute value of a Sollay range.
1515 226 storres
    """
1516 226 storres
    lowerBoundSo = sollya_lib_inf(rangeSo)
1517 226 storres
    upperBoundSo = sollya_lib_sup(rangeSo)
1518 226 storres
    #
1519 226 storres
    lowerBoundSo = sollya_lib_build_function_abs(lowerBoundSo)
1520 226 storres
    upperBoundSo = sollya_lib_build_function_abs(upperBoundSo)
1521 226 storres
    #pobyso_autoprint(lowerBoundSo)
1522 226 storres
    #pobyso_autoprint(upperBoundSo)
1523 226 storres
    #
1524 226 storres
    maxAbsSo = sollya_lib_max(lowerBoundSo, upperBoundSo, None)
1525 226 storres
    #sollya_lib_clear_obj(lowerBoundSo)
1526 226 storres
    #sollya_lib_clear_obj(upperBoundSo)
1527 226 storres
    return maxAbsSo
1528 226 storres
# End pobyso_range_max_abs_so_so
1529 226 storres
1530 85 storres
def pobyso_range_to_interval_so_sa(rangeSo, realIntervalFieldSa = None):
1531 83 storres
    """
1532 83 storres
    Get a Sage interval from a Sollya range.
1533 83 storres
    If no realIntervalField is given as a parameter, the Sage interval
1534 83 storres
    precision is that of the Sollya range.
1535 85 storres
    Otherwise, the precision is that of the realIntervalField. In this case
1536 85 storres
    rounding may happen.
1537 83 storres
    """
1538 85 storres
    if realIntervalFieldSa is None:
1539 56 storres
        precSa = pobyso_get_prec_of_range_so_sa(rangeSo)
1540 85 storres
        realIntervalFieldSa = RealIntervalField(precSa)
1541 56 storres
    intervalSa = \
1542 85 storres
        pobyso_get_interval_from_range_so_sa(rangeSo, realIntervalFieldSa)
1543 209 storres
    return intervalSa
1544 209 storres
# End pobyso_range_to_interval_so_sa
1545 240 storres
#
1546 240 storres
def pobyso_relative_so_so():
1547 240 storres
    """
1548 240 storres
    Very thin wrapper around the sollya_lib_relative function.
1549 240 storres
    """
1550 240 storres
    return sollya_lib_relative()
1551 240 storres
# End pobyso_relative_so_so
1552 240 storres
#
1553 209 storres
def pobyso_rat_poly_sa_so(polySa, precSa = None):
1554 209 storres
    """
1555 209 storres
    Create a Sollya polynomial from a Sage rational polynomial.
1556 209 storres
    """
1557 209 storres
    ## TODO: filter arguments.
1558 209 storres
    ## Precision. If no precision is given, use the current precision
1559 209 storres
    #  of Sollya.
1560 209 storres
    if precSa is None:
1561 209 storres
        precSa =  pobyso_get_prec_so_sa()
1562 209 storres
    #print "Precision:",  precSa
1563 209 storres
    RRR = RealField(precSa)
1564 209 storres
    ## Create a Sage polynomial in the "right" precision.
1565 209 storres
    P_RRR = RRR[polySa.variables()[0]]
1566 209 storres
    polyFloatSa = P_RRR(polySa)
1567 213 storres
    ## Make sure no precision is provided: pobyso_float_poly_sa_so will
1568 213 storres
    #  recover it all by itself and not make an extra conversion.
1569 209 storres
    return pobyso_float_poly_sa_so(polyFloatSa)
1570 209 storres
1571 209 storres
# End pobyso_rat_poly_sa_so
1572 209 storres
1573 52 storres
def pobyso_remez_canonical_sa_sa(func, \
1574 52 storres
                                 degree, \
1575 52 storres
                                 lowerBound, \
1576 52 storres
                                 upperBound, \
1577 52 storres
                                 weight = None, \
1578 52 storres
                                 quality = None):
1579 52 storres
    """
1580 52 storres
    All arguments are Sage/Python.
1581 52 storres
    The functions (func and weight) must be passed as expressions or strings.
1582 52 storres
    Otherwise the function fails.
1583 83 storres
    The return value is a Sage polynomial.
1584 52 storres
    """
1585 83 storres
    var('zorglub')    # Dummy variable name for type check only. Type of
1586 83 storres
    # zorglub is "symbolic expression".
1587 52 storres
    polySo = pobyso_remez_canonical_sa_so(func, \
1588 52 storres
                                 degree, \
1589 52 storres
                                 lowerBound, \
1590 52 storres
                                 upperBound, \
1591 85 storres
                                 weight, \
1592 85 storres
                                 quality)
1593 83 storres
    # String test
1594 52 storres
    if parent(func) == parent("string"):
1595 52 storres
        functionSa = eval(func)
1596 52 storres
    # Expression test.
1597 52 storres
    elif type(func) == type(zorglub):
1598 52 storres
        functionSa = func
1599 83 storres
    else:
1600 83 storres
        return None
1601 83 storres
    #
1602 52 storres
    maxPrecision = 0
1603 52 storres
    if polySo is None:
1604 52 storres
        return(None)
1605 52 storres
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
1606 85 storres
    RRRRSa = RealField(maxPrecision)
1607 85 storres
    polynomialRingSa = RRRRSa[functionSa.variables()[0]]
1608 85 storres
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, RRRRSa)
1609 85 storres
    polySa = polynomial(expSa, polynomialRingSa)
1610 83 storres
    sollya_lib_clear_obj(polySo)
1611 52 storres
    return(polySa)
1612 85 storres
# End pobyso_remez_canonical_sa_sa
1613 52 storres
1614 38 storres
def pobyso_remez_canonical(func, \
1615 5 storres
                           degree, \
1616 5 storres
                           lowerBound, \
1617 5 storres
                           upperBound, \
1618 38 storres
                           weight = "1", \
1619 5 storres
                           quality = None):
1620 38 storres
    """ Legacy function. See pobyso_remez_canonical_sa_so. """
1621 51 storres
    return(pobyso_remez_canonical_sa_so(func, \
1622 51 storres
                                        degree, \
1623 51 storres
                                        lowerBound, \
1624 51 storres
                                        upperBound, \
1625 51 storres
                                        weight, \
1626 51 storres
                                        quality))
1627 200 storres
# End pobyso_remez_canonical.
1628 200 storres
1629 38 storres
def pobyso_remez_canonical_sa_so(func, \
1630 38 storres
                                 degree, \
1631 38 storres
                                 lowerBound, \
1632 38 storres
                                 upperBound, \
1633 52 storres
                                 weight = None, \
1634 38 storres
                                 quality = None):
1635 38 storres
    """
1636 38 storres
    All arguments are Sage/Python.
1637 51 storres
    The functions (func and weight) must be passed as expressions or strings.
1638 51 storres
    Otherwise the function fails.
1639 38 storres
    The return value is a pointer to a Sollya function.
1640 234 storres
    lowerBound and upperBound mus be reals.
1641 38 storres
    """
1642 83 storres
    var('zorglub')    # Dummy variable name for type check only. Type of
1643 83 storres
    # zorglub is "symbolic expression".
1644 85 storres
    currentVariableNameSa = None
1645 52 storres
    # The func argument can be of different types (string,
1646 52 storres
    # symbolic expression...)
1647 38 storres
    if parent(func) == parent("string"):
1648 234 storres
        localFuncSa = sage_eval(func,globals())
1649 85 storres
        if len(localFuncSa.variables()) > 0:
1650 85 storres
            currentVariableNameSa = localFuncSa.variables()[0]
1651 85 storres
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1652 159 storres
            functionSo = \
1653 159 storres
              sollya_lib_parse_string(localFuncSa._assume_str().replace('_SAGE_VAR_', ''))
1654 51 storres
    # Expression test.
1655 52 storres
    elif type(func) == type(zorglub):
1656 52 storres
        # Until we are able to translate Sage expressions into Sollya
1657 52 storres
        # expressions : parse the string version.
1658 85 storres
        if len(func.variables()) > 0:
1659 85 storres
            currentVariableNameSa = func.variables()[0]
1660 85 storres
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1661 159 storres
            functionSo = \
1662 159 storres
              sollya_lib_parse_string(func._assume_str().replace('_SAGE_VAR_', ''))
1663 38 storres
    else:
1664 38 storres
        return(None)
1665 85 storres
    if weight is None: # No weight given -> 1.
1666 52 storres
        weightSo = pobyso_constant_1_sa_so()
1667 85 storres
    elif parent(weight) == parent("string"): # Weight given as string: parse it.
1668 51 storres
        weightSo = sollya_lib_parse_string(func)
1669 85 storres
    elif type(weight) == type(zorglub): # Weight given as symbolice expression.
1670 159 storres
        functionSo = \
1671 159 storres
          sollya_lib_parse_string_sa_so(weight._assume_str().replace('_SAGE_VAR_', ''))
1672 51 storres
    else:
1673 51 storres
        return(None)
1674 5 storres
    degreeSo = pobyso_constant_from_int(degree)
1675 85 storres
    rangeSo = pobyso_bounds_to_range_sa_so(lowerBound, upperBound)
1676 38 storres
    if not quality is None:
1677 38 storres
        qualitySo= pobyso_constant_sa_so(quality)
1678 52 storres
    else:
1679 52 storres
        qualitySo = None
1680 83 storres
1681 83 storres
    remezPolySo = sollya_lib_remez(functionSo, \
1682 83 storres
                                   degreeSo, \
1683 83 storres
                                   rangeSo, \
1684 83 storres
                                   weightSo, \
1685 83 storres
                                   qualitySo, \
1686 83 storres
                                   None)
1687 83 storres
    sollya_lib_clear_obj(functionSo)
1688 83 storres
    sollya_lib_clear_obj(degreeSo)
1689 83 storres
    sollya_lib_clear_obj(rangeSo)
1690 83 storres
    sollya_lib_clear_obj(weightSo)
1691 83 storres
    if not qualitySo is None:
1692 85 storres
        sollya_lib_clear_obj(qualitySo)
1693 83 storres
    return(remezPolySo)
1694 83 storres
# End pobyso_remez_canonical_sa_so
1695 83 storres
1696 38 storres
def pobyso_remez_canonical_so_so(funcSo, \
1697 38 storres
                                 degreeSo, \
1698 38 storres
                                 rangeSo, \
1699 52 storres
                                 weightSo = pobyso_constant_1_sa_so(),\
1700 38 storres
                                 qualitySo = None):
1701 38 storres
    """
1702 38 storres
    All arguments are pointers to Sollya objects.
1703 38 storres
    The return value is a pointer to a Sollya function.
1704 38 storres
    """
1705 38 storres
    if not sollya_lib_obj_is_function(funcSo):
1706 38 storres
        return(None)
1707 38 storres
    return(sollya_lib_remez(funcSo, degreeSo, rangeSo, weightSo, qualitySo, None))
1708 200 storres
# End pobyso_remez_canonical_so_so.
1709 200 storres
1710 237 storres
def pobyso_remez_exponents_list_sa_so(func,     \
1711 237 storres
                                 exponentsList, \
1712 237 storres
                                 lowerBound,    \
1713 237 storres
                                 upperBound,    \
1714 237 storres
                                 weight = None, \
1715 237 storres
                                 quality = None):
1716 237 storres
    """
1717 237 storres
    All arguments are Sage/Python.
1718 237 storres
    The functions (func and weight) must be passed as expressions or strings.
1719 237 storres
    Otherwise the function fails.
1720 237 storres
    The return value is a pointer to a Sollya function.
1721 237 storres
    lowerBound and upperBound mus be reals.
1722 237 storres
    """
1723 237 storres
    var('zorglub')    # Dummy variable name for type check only. Type of
1724 237 storres
    # zorglub is "symbolic expression".
1725 237 storres
    currentVariableNameSa = None
1726 237 storres
    # The func argument can be of different types (string,
1727 237 storres
    # symbolic expression...)
1728 237 storres
    if parent(func) == parent("string"):
1729 237 storres
        localFuncSa = sage_eval(func,globals())
1730 237 storres
        if len(localFuncSa.variables()) > 0:
1731 237 storres
            currentVariableNameSa = localFuncSa.variables()[0]
1732 237 storres
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1733 237 storres
            functionSo = \
1734 237 storres
              sollya_lib_parse_string(localFuncSa._assume_str().replace('_SAGE_VAR_', ''))
1735 237 storres
    # Expression test.
1736 237 storres
    elif type(func) == type(zorglub):
1737 237 storres
        # Until we are able to translate Sage expressions into Sollya
1738 237 storres
        # expressions : parse the string version.
1739 237 storres
        if len(func.variables()) > 0:
1740 237 storres
            currentVariableNameSa = func.variables()[0]
1741 237 storres
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1742 237 storres
            functionSo = \
1743 237 storres
              sollya_lib_parse_string(func._assume_str().replace('_SAGE_VAR_', ''))
1744 237 storres
    else:
1745 237 storres
        return(None)
1746 237 storres
    ## Deal with the weight, much in the same way as with the function.
1747 237 storres
    if weight is None: # No weight given -> 1.
1748 237 storres
        weightSo = pobyso_constant_1_sa_so()
1749 237 storres
    elif parent(weight) == parent("string"): # Weight given as string: parse it.
1750 237 storres
        weightSo = sollya_lib_parse_string(func)
1751 237 storres
    elif type(weight) == type(zorglub): # Weight given as symbolice expression.
1752 237 storres
        functionSo = \
1753 237 storres
          sollya_lib_parse_string(weight._assume_str().replace('_SAGE_VAR_', '',100))
1754 237 storres
    else:
1755 237 storres
        return(None)
1756 237 storres
    rangeSo = pobyso_bounds_to_range_sa_so(lowerBound, upperBound)
1757 237 storres
    if not quality is None:
1758 237 storres
        qualitySo= pobyso_constant_sa_so(quality)
1759 237 storres
    else:
1760 237 storres
        qualitySo = None
1761 237 storres
    #
1762 237 storres
    ## Tranform the Sage list of exponents into a Sollya list.
1763 237 storres
    exponentsListSo = pobyso_build_list_of_ints_sa_so(*exponentsList)
1764 237 storres
    remezPolySo = sollya_lib_remez(functionSo, \
1765 237 storres
                                   exponentsListSo, \
1766 237 storres
                                   rangeSo, \
1767 237 storres
                                   weightSo, \
1768 237 storres
                                   qualitySo, \
1769 237 storres
                                   None)
1770 237 storres
    sollya_lib_clear_obj(functionSo)
1771 237 storres
    sollya_lib_clear_obj(exponentsListSo)
1772 237 storres
    sollya_lib_clear_obj(rangeSo)
1773 237 storres
    sollya_lib_clear_obj(weightSo)
1774 237 storres
    if not qualitySo is None:
1775 237 storres
        sollya_lib_clear_obj(qualitySo)
1776 237 storres
    return(remezPolySo)
1777 237 storres
# End pobyso_remez_exponentsList_sa_so
1778 240 storres
#
1779 240 storres
def pobyso_round_coefficients_so_so(polySo, truncFormatListSo):
1780 240 storres
    """
1781 240 storres
    A wrapper around the "classical" sollya_lib_roundcoefficients: a Sollya
1782 240 storres
    polynomial and a Sollya list are given as arguments.
1783 240 storres
    """
1784 240 storres
    return sollya_lib_roundcoefficients(polySo, truncFormatListSo)
1785 237 storres
1786 224 storres
def pobyso_round_coefficients_progressive_so_so(polySo,
1787 224 storres
                                                funcSo,
1788 224 storres
                                                precSo,
1789 224 storres
                                                intervalSo,
1790 224 storres
                                                icSo,
1791 224 storres
                                                currentApproxErrorSo,
1792 224 storres
                                                approxAccurSo,
1793 224 storres
                                                debug=False):
1794 223 storres
    if debug:
1795 223 storres
        print "Input arguments:"
1796 228 storres
        print "Polynomial: ", ; pobyso_autoprint(polySo)
1797 228 storres
        print "Function: ", ; pobyso_autoprint(funcSo)
1798 228 storres
        print "Internal precision: ", ; pobyso_autoprint(precSo)
1799 228 storres
        print "Interval: ", ; pobyso_autoprint(intervalSo)
1800 228 storres
        print "Current approximation error: ", ; pobyso_autoprint(currentApproxErrorSo)
1801 228 storres
        print "Requested approxiation error: ", ; pobyso_autoprint(approxAccurSo)
1802 223 storres
        print "________________"
1803 223 storres
    approxAccurSa        = pobyso_get_constant_as_rn_so_sa(approxAccurSo)
1804 223 storres
    currentApproxErrorSa = pobyso_get_constant_as_rn_so_sa(currentApproxErrorSo)
1805 223 storres
    ## If the current approximation error is too close to the target, there is
1806 223 storres
    #  no possible gain.
1807 223 storres
    if currentApproxErrorSa >= approxAccurSa / 2:
1808 232 storres
        #### Do not return the initial argument but copies: caller may free
1809 232 storres
        #    as inutile after call.
1810 232 storres
        return (sollya_lib_copy_obj(polySo),
1811 232 storres
                sollya_lib_copy_obj(currentApproxErrorSo))
1812 223 storres
    degreeSa             = pobyso_polynomial_degree_so_sa(polySo)
1813 223 storres
    intervalSa           = pobyso_range_to_interval_so_sa(intervalSo)
1814 223 storres
1815 223 storres
    if debug:
1816 228 storres
        print "degreeSa              :", degreeSa
1817 228 storres
        print "intervalSa            :", intervalSa.str(style='brackets')
1818 223 storres
        print "currentApproxErrorSa  :", currentApproxErrorSa
1819 228 storres
        print "approxAccurSa         :", approxAccurSa
1820 223 storres
    ### Start with a 0 value expression.
1821 223 storres
    radiusSa = intervalSa.absolute_diameter() / 2
1822 223 storres
    if debug:
1823 224 storres
        print "log2(radius):", RR(radiusSa).log2()
1824 224 storres
    iterIndex = 0
1825 224 storres
    while True:
1826 224 storres
        resPolySo = pobyso_constant_0_sa_so()
1827 224 storres
        roundedPolyApproxAccurSa = approxAccurSa / 2
1828 224 storres
        currentRadiusPowerSa = 1
1829 224 storres
        for degree in xrange(0,degreeSa + 1):
1830 224 storres
            #### At round 0, use the agressive formula. At round 1, run the
1831 224 storres
            #    proved formula.
1832 224 storres
            if iterIndex == 0:
1833 224 storres
                roundingPowerSa = \
1834 224 storres
                    floor(((currentRadiusPowerSa/roundedPolyApproxAccurSa)*(degree+1)).log2())
1835 224 storres
            else:
1836 224 storres
                roundingPowerSa = \
1837 224 storres
                    floor(((currentRadiusPowerSa/roundedPolyApproxAccurSa)*(degreeSa+1)).log2())
1838 228 storres
            ## Under extreme conditions the above formulas can evaluate under 2, which is the
1839 228 storres
            #  minimal precision of an MPFR number.
1840 228 storres
            if roundingPowerSa < 2:
1841 228 storres
                roundingPowerSa = 2
1842 224 storres
            if debug:
1843 224 storres
                print "roundedPolyApproxAccurSa", roundedPolyApproxAccurSa
1844 224 storres
                print "currentRadiusPowerSa", currentRadiusPowerSa
1845 224 storres
                print "Current rounding exponent:", roundingPowerSa
1846 224 storres
            currentRadiusPowerSa *= radiusSa
1847 224 storres
            index1So = pobyso_constant_from_int_sa_so(degree)
1848 224 storres
            index2So = pobyso_constant_from_int_sa_so(degree)
1849 224 storres
            ### Create a monomial with:
1850 224 storres
            #   - the coefficient in the initial monomial at the current degrree;
1851 224 storres
            #   - the current exponent;
1852 224 storres
            #   - the free variable.
1853 224 storres
            cmonSo  = \
1854 224 storres
                sollya_lib_build_function_mul(sollya_lib_coeff(polySo, index1So),
1855 224 storres
                                              sollya_lib_build_function_pow( \
1856 224 storres
                                                sollya_lib_build_function_free_variable(), \
1857 224 storres
                                                index2So))
1858 224 storres
            roundingPowerSo = pobyso_constant_from_int_sa_so(roundingPowerSa)
1859 224 storres
            cmonrSo = pobyso_round_coefficients_single_so_so(cmonSo, roundingPowerSo)
1860 224 storres
            sollya_lib_clear_obj(cmonSo)
1861 224 storres
            ### Add to the result polynomial.
1862 224 storres
            resPolySo = sollya_lib_build_function_add(resPolySo,
1863 224 storres
                                                      cmonrSo)
1864 224 storres
        # End for.
1865 224 storres
        ### Check the new polynomial.
1866 224 storres
        freeVarSo     = sollya_lib_build_function_free_variable()
1867 224 storres
        changeVarSo   = sollya_lib_sub(freeVarSo, icSo)
1868 224 storres
        resPolyCvSo   = sollya_lib_evaluate(resPolySo, changeVarSo)
1869 224 storres
        errFuncSo     = sollya_lib_build_function_sub(sollya_lib_copy_obj(funcSo),
1870 224 storres
                                                      resPolyCvSo)
1871 224 storres
        infNormSo = sollya_lib_dirtyinfnorm(errFuncSo, intervalSo)
1872 224 storres
        ### This also clears resPolyCvSo.
1873 224 storres
        sollya_lib_clear_obj(errFuncSo)
1874 224 storres
        cerrSa    = pobyso_get_constant_as_rn_so_sa(infNormSo)
1875 223 storres
        if debug:
1876 224 storres
            print "Error of the new polynomial:", cerrSa
1877 224 storres
        ### If at round 1, return the initial polynomial error. This should
1878 224 storres
        #   never happen since the rounding algorithm is proved. But some
1879 224 storres
        #   circumstances may break it (e.g. internal precision of tools).
1880 224 storres
        if cerrSa > approxAccurSa:
1881 224 storres
            if iterIndex > 0: # Round 1 and beyond.
1882 224 storres
                sollya_lib_clear_obj(resPolySo)
1883 224 storres
                sollya_lib_clear_obj(infNormSo)
1884 232 storres
                #### Do not return the arguments but copies: the caller may free
1885 232 storres
                #    free the former as inutile after call.
1886 232 storres
                return (sollya_lib_copy_obj(polySo),
1887 232 storres
                        sollya_lib_copy_obj(currentApproxErrorSo))
1888 224 storres
            else: # Round 0, got round 1
1889 224 storres
                sollya_lib_clear_obj(resPolySo)
1890 224 storres
                sollya_lib_clear_obj(infNormSo)
1891 224 storres
                iterIndex += 1
1892 224 storres
                continue
1893 224 storres
        ### If get here it is because cerrSa <= approxAccurSa
1894 224 storres
        ### Approximation error of the new polynomial is acceptable.
1895 224 storres
        return (resPolySo, infNormSo)
1896 224 storres
    # End while True
1897 224 storres
# End pobyso_round_coefficients_progressive_so_so
1898 223 storres
1899 240 storres
def pobyso_round_coefficients_single_so_so(polySo, commonPrecSo):
1900 215 storres
    """
1901 215 storres
    Create a rounded coefficients polynomial from polynomial argument to
1902 215 storres
    the number of bits in size argument.
1903 215 storres
    All coefficients are set to the same precision.
1904 215 storres
    """
1905 215 storres
    ## TODO: check arguments.
1906 240 storres
    endEllipListSo = pobyso_build_end_elliptic_list_so_so(commonPrecSo)
1907 215 storres
    polySo = sollya_lib_roundcoefficients(polySo, endEllipListSo, None)
1908 217 storres
    sollya_lib_clear_obj(endEllipListSo)
1909 215 storres
    #sollya_lib_clear_obj(endEllipListSo)
1910 215 storres
    return polySo
1911 215 storres
1912 215 storres
# End pobyso_round_coefficients_single_so_so
1913 215 storres
1914 5 storres
def pobyso_set_canonical_off():
1915 5 storres
    sollya_lib_set_canonical(sollya_lib_off())
1916 5 storres
1917 5 storres
def pobyso_set_canonical_on():
1918 5 storres
    sollya_lib_set_canonical(sollya_lib_on())
1919 5 storres
1920 5 storres
def pobyso_set_prec(p):
1921 38 storres
    """ Legacy function. See pobyso_set_prec_sa_so. """
1922 85 storres
    pobyso_set_prec_sa_so(p)
1923 38 storres
1924 38 storres
def pobyso_set_prec_sa_so(p):
1925 227 storres
    #a = c_int(p)
1926 226 storres
    #precSo = c_void_p(sollya_lib_constant_from_int(a))
1927 227 storres
    #precSo = sollya_lib_constant_from_int(a)
1928 227 storres
    precSo = pobyso_constant_from_int_sa_so(p)
1929 226 storres
    sollya_lib_set_prec(precSo)
1930 226 storres
    sollya_lib_clear_obj(precSo)
1931 215 storres
# End pobyso_set_prec_sa_so.
1932 5 storres
1933 85 storres
def pobyso_set_prec_so_so(newPrecSo):
1934 226 storres
    sollya_lib_set_prec(newPrecSo)
1935 215 storres
# End pobyso_set_prec_so_so.
1936 54 storres
1937 215 storres
def pobyso_inf_so_so(intervalSo):
1938 215 storres
    """
1939 215 storres
    Very thin wrapper around sollya_lib_inf().
1940 215 storres
    """
1941 215 storres
    return sollya_lib_inf(intervalSo)
1942 215 storres
# End pobyso_inf_so_so.
1943 242 storres
#
1944 242 storres
def pobyso_supnorm_so_sa(polySo, funcSo, intervalSo, errorTypeSo = None,\
1945 242 storres
                         accuracySo = None, realFieldSa = None):
1946 242 storres
    """
1947 242 storres
    Computes the supremum norm from Solly input arguments and returns a
1948 242 storres
    Sage floating-point number whose precision is set by the only Sage argument.
1949 215 storres

1950 242 storres
    The returned value is the maximum of the absolute values of the range
1951 242 storres
    elements  returned by the Sollya supnorm functions
1952 242 storres
    """
1953 242 storres
    supNormRangeSo = pobyso_supnorm_so_so(polySo,
1954 242 storres
                                          funcSo,
1955 242 storres
                                          intervalSo,
1956 242 storres
                                          errorTypeSo,
1957 242 storres
                                          accuracySo)
1958 242 storres
    supNormSo = pobyso_range_max_abs_so_so(supNormRangeSo)
1959 242 storres
    sollya_lib_clear_obj(supNormRangeSo)
1960 242 storres
    #pobyso_autoprint(supNormSo)
1961 242 storres
    supNormSa = pobyso_get_constant_as_rn_with_rf_so_sa(supNormSo, realFieldSa)
1962 242 storres
    sollya_lib_clear_obj(supNormSo)
1963 242 storres
    return supNormSa
1964 242 storres
# End pobyso_supnorm_so_sa.
1965 242 storres
#
1966 85 storres
def pobyso_supnorm_so_so(polySo, funcSo, intervalSo, errorTypeSo = None,\
1967 85 storres
                         accuracySo = None):
1968 58 storres
    """
1969 85 storres
    Computes the supnorm of the approximation error between the given
1970 242 storres
    polynomial and function. Attention: returns a range!
1971 85 storres
    errorTypeSo defaults to "absolute".
1972 85 storres
    accuracySo defaults to 2^(-40).
1973 85 storres
    """
1974 85 storres
    if errorTypeSo is None:
1975 85 storres
        errorTypeSo = sollya_lib_absolute(None)
1976 85 storres
        errorTypeIsNone = True
1977 85 storres
    else:
1978 85 storres
        errorTypeIsNone = False
1979 85 storres
    #
1980 85 storres
    if accuracySo is None:
1981 237 storres
        # Notice the **: we are in Pythonland here!
1982 85 storres
        accuracySo = pobyso_constant_sa_so(RR(2**(-40)))
1983 85 storres
        accuracyIsNone = True
1984 85 storres
    else:
1985 85 storres
        accuracyIsNone = False
1986 238 storres
    #pobyso_autoprint(accuracySo)
1987 85 storres
    resultSo = \
1988 85 storres
        sollya_lib_supnorm(polySo, funcSo, intervalSo, errorTypeSo, \
1989 85 storres
                              accuracySo)
1990 85 storres
    if errorTypeIsNone:
1991 85 storres
        sollya_lib_clear_obj(errorTypeSo)
1992 85 storres
    if accuracyIsNone:
1993 85 storres
        sollya_lib_clear_obj(accuracySo)
1994 85 storres
    return resultSo
1995 85 storres
# End pobyso_supnorm_so_so
1996 242 storres
#
1997 162 storres
def pobyso_taylor_expansion_no_change_var_so_so(functionSo,
1998 162 storres
                                                degreeSo,
1999 162 storres
                                                rangeSo,
2000 162 storres
                                                errorTypeSo=None,
2001 162 storres
                                                sollyaPrecSo=None):
2002 85 storres
    """
2003 162 storres
    Compute the Taylor expansion without the variable change
2004 162 storres
    x -> x-intervalCenter.
2005 58 storres
    """
2006 226 storres
    # Change internal Sollya precision, if needed.
2007 227 storres
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
2008 226 storres
    sollyaPrecChanged = False
2009 226 storres
    if sollyaPrecSo is None:
2010 226 storres
        pass
2011 226 storres
    else:
2012 58 storres
        sollya_lib_set_prec(sollyaPrecSo)
2013 226 storres
        sollyaPrecChanged = True
2014 85 storres
    # Error type stuff: default to absolute.
2015 85 storres
    if errorTypeSo is None:
2016 85 storres
        errorTypeIsNone = True
2017 85 storres
        errorTypeSo = sollya_lib_absolute(None)
2018 85 storres
    else:
2019 85 storres
        errorTypeIsNone = False
2020 162 storres
    intervalCenterSo = sollya_lib_mid(rangeSo, None)
2021 162 storres
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo,
2022 162 storres
                                         intervalCenterSo,
2023 58 storres
                                         rangeSo, errorTypeSo, None)
2024 117 storres
    # taylorFormListSaSo is a Python list of Sollya objects references that
2025 117 storres
    # are copies of the elements of taylorFormSo.
2026 117 storres
    # pobyso_get_list_elements_so_so clears taylorFormSo.
2027 162 storres
    (taylorFormListSaSo, numElementsSa, isEndEllipticSa) = \
2028 58 storres
        pobyso_get_list_elements_so_so(taylorFormSo)
2029 162 storres
    polySo = sollya_lib_copy_obj(taylorFormListSaSo[0])
2030 162 storres
    #print "Num elements:", numElementsSa
2031 162 storres
    sollya_lib_clear_obj(taylorFormSo)
2032 162 storres
    #polySo = taylorFormListSaSo[0]
2033 162 storres
    #errorRangeSo = sollya_lib_copy_obj(taylorFormListSaSo[2])
2034 162 storres
    errorRangeSo = taylorFormListSaSo[2]
2035 181 storres
    # No copy_obj needed here: a new objects are created.
2036 181 storres
    maxErrorSo    = sollya_lib_sup(errorRangeSo)
2037 181 storres
    minErrorSo    = sollya_lib_inf(errorRangeSo)
2038 181 storres
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
2039 181 storres
    absMinErrorSo = sollya_lib_abs(minErrorSo)
2040 181 storres
    sollya_lib_clear_obj(maxErrorSo)
2041 181 storres
    sollya_lib_clear_obj(minErrorSo)
2042 181 storres
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
2043 181 storres
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
2044 58 storres
    # If changed, reset the Sollya working precision.
2045 226 storres
    if sollyaPrecChanged:
2046 226 storres
        sollya_lib_set_prec(initialSollyaPrecSo)
2047 226 storres
    sollya_lib_clear_obj(initialSollyaPrecSo)
2048 85 storres
    if errorTypeIsNone:
2049 85 storres
        sollya_lib_clear_obj(errorTypeSo)
2050 162 storres
    pobyso_clear_taylorform_sa_so(taylorFormListSaSo)
2051 181 storres
    if absMaxErrorSa > absMinErrorSa:
2052 181 storres
        sollya_lib_clear_obj(absMinErrorSo)
2053 227 storres
        return (polySo, intervalCenterSo, absMaxErrorSo)
2054 181 storres
    else:
2055 181 storres
        sollya_lib_clear_obj(absMaxErrorSo)
2056 227 storres
        return (polySo, intervalCenterSo, absMinErrorSo)
2057 162 storres
# end pobyso_taylor_expansion_no_change_var_so_so
2058 58 storres
2059 162 storres
def pobyso_taylor_expansion_with_change_var_so_so(functionSo, degreeSo, \
2060 162 storres
                                                  rangeSo, \
2061 162 storres
                                                  errorTypeSo=None, \
2062 162 storres
                                                  sollyaPrecSo=None):
2063 58 storres
    """
2064 162 storres
    Compute the Taylor expansion with the variable change
2065 162 storres
    x -> (x-intervalCenter) included.
2066 58 storres
    """
2067 226 storres
    # Change Sollya internal precision, if need.
2068 226 storres
    sollyaPrecChanged = False
2069 226 storres
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_sos_sa()
2070 226 storres
    if sollyaPrecSo is None:
2071 226 storres
        pass
2072 226 storres
    else:
2073 56 storres
        sollya_lib_set_prec(sollyaPrecSo)
2074 226 storres
        sollyaPrecChanged = True
2075 162 storres
    #
2076 85 storres
    # Error type stuff: default to absolute.
2077 85 storres
    if errorTypeSo is None:
2078 85 storres
        errorTypeIsNone = True
2079 85 storres
        errorTypeSo = sollya_lib_absolute(None)
2080 85 storres
    else:
2081 85 storres
        errorTypeIsNone = False
2082 162 storres
    intervalCenterSo = sollya_lib_mid(rangeSo)
2083 162 storres
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo, \
2084 162 storres
                                         intervalCenterSo, \
2085 56 storres
                                         rangeSo, errorTypeSo, None)
2086 116 storres
    # taylorFormListSaSo is a Python list of Sollya objects references that
2087 116 storres
    # are copies of the elements of taylorFormSo.
2088 116 storres
    # pobyso_get_list_elements_so_so clears taylorFormSo.
2089 162 storres
    (taylorFormListSo, numElements, isEndElliptic) = \
2090 56 storres
        pobyso_get_list_elements_so_so(taylorFormSo)
2091 162 storres
    polySo = taylorFormListSo[0]
2092 162 storres
    errorRangeSo = taylorFormListSo[2]
2093 181 storres
    maxErrorSo    = sollya_lib_sup(errorRangeSo)
2094 181 storres
    minErrorSo    = sollya_lib_inf(errorRangeSo)
2095 181 storres
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
2096 181 storres
    absMinErrorSo = sollya_lib_abs(minErrorSo)
2097 181 storres
    sollya_lib_clear_obj(maxErrorSo)
2098 181 storres
    sollya_lib_clear_obj(minErrorSo)
2099 181 storres
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
2100 181 storres
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
2101 162 storres
    changeVarExpSo = sollya_lib_build_function_sub(\
2102 162 storres
                       sollya_lib_build_function_free_variable(),\
2103 162 storres
                       sollya_lib_copy_obj(intervalCenterSo))
2104 181 storres
    polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpSo)
2105 181 storres
    sollya_lib_clear_obj(polySo)
2106 162 storres
    sollya_lib_clear_obj(changeVarExpSo)
2107 56 storres
    # If changed, reset the Sollya working precision.
2108 226 storres
    if sollyaPrecChanged:
2109 226 storres
        sollya_lib_set_prec(initialSollyaPrecSo)
2110 226 storres
    sollya_lib_clear_obj(initialSollyaPrecSo)
2111 85 storres
    if errorTypeIsNone:
2112 85 storres
        sollya_lib_clear_obj(errorTypeSo)
2113 162 storres
    sollya_lib_clear_obj(taylorFormSo)
2114 162 storres
    # Do not clear maxErrorSo.
2115 181 storres
    if absMaxErrorSa > absMinErrorSa:
2116 181 storres
        sollya_lib_clear_obj(absMinErrorSo)
2117 181 storres
        return((polyVarChangedSo, intervalCenterSo, absMaxErrorSo))
2118 181 storres
    else:
2119 181 storres
        sollya_lib_clear_obj(absMaxErrorSo)
2120 181 storres
        return((polyVarChangedSo, intervalCenterSo, absMinErrorSo))
2121 162 storres
# end pobyso_taylor_expansion_with_change_var_so_so
2122 56 storres
2123 5 storres
def pobyso_taylor(function, degree, point):
2124 38 storres
    """ Legacy function. See pobysoTaylor_so_so. """
2125 38 storres
    return(pobyso_taylor_so_so(function, degree, point))
2126 38 storres
2127 56 storres
def pobyso_taylor_so_so(functionSo, degreeSo, pointSo):
2128 56 storres
    return(sollya_lib_taylor(functionSo, degreeSo, pointSo))
2129 5 storres
2130 85 storres
def pobyso_taylorform(function, degree, point = None,
2131 85 storres
                      interval = None, errorType=None):
2132 85 storres
    """ Legacy function. See pobyso_taylorform_sa_sa;"""
2133 38 storres
2134 38 storres
def pobyso_taylorform_sa_sa(functionSa, \
2135 84 storres
                            degreeSa, \
2136 84 storres
                            pointSa, \
2137 84 storres
                            intervalSa=None, \
2138 84 storres
                            errorTypeSa=None, \
2139 84 storres
                            precisionSa=None):
2140 37 storres
    """
2141 85 storres
    Compute the Taylor form of 'degreeSa' for 'functionSa' at 'pointSa'
2142 85 storres
    for 'intervalSa' with 'errorTypeSa' (a string) using 'precisionSa'.
2143 37 storres
    point: must be a Real or a Real interval.
2144 37 storres
    return the Taylor form as an array
2145 83 storres
    TODO: take care of the interval and of the point when it is an interval;
2146 38 storres
          when errorType is not None;
2147 83 storres
          take care of the other elements of the Taylor form (coefficients
2148 83 storres
          errors and delta.
2149 37 storres
    """
2150 37 storres
    # Absolute as the default error.
2151 84 storres
    if errorTypeSa is None:
2152 37 storres
        errorTypeSo = sollya_lib_absolute()
2153 84 storres
    elif errorTypeSa == "relative":
2154 84 storres
        errorTypeSo = sollya_lib_relative()
2155 84 storres
    elif errortypeSa == "absolute":
2156 84 storres
        errorTypeSo = sollya_lib_absolute()
2157 37 storres
    else:
2158 84 storres
        # No clean up needed.
2159 84 storres
        return None
2160 84 storres
    # Global precision stuff
2161 226 storres
    sollyaPrecisionChangedSa = False
2162 226 storres
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
2163 226 storres
    if precisionSa is None:
2164 226 storres
        precSa = initialSollyaPrecSa
2165 226 storres
    else:
2166 226 storres
        if precSa > initialSollyaPrecSa:
2167 226 storres
            if precSa <= 2:
2168 226 storres
                print inspect.stack()[0][3], ":precision change <= 2 requested."
2169 226 storres
            pobyso_set_prec_sa_so(precSa)
2170 226 storres
            sollyaPrecisionChangedSa = True
2171 226 storres
    #
2172 85 storres
    if len(functionSa.variables()) > 0:
2173 85 storres
        varSa = functionSa.variables()[0]
2174 85 storres
        pobyso_name_free_variable_sa_so(str(varSa))
2175 84 storres
    # In any case (point or interval) the parent of pointSa has a precision
2176 84 storres
    # method.
2177 84 storres
    pointPrecSa = pointSa.parent().precision()
2178 226 storres
    if precSa > pointPrecSa:
2179 226 storres
        pointPrecSa = precSa
2180 84 storres
    # In any case (point or interval) pointSa has a base_ring() method.
2181 84 storres
    pointBaseRingString = str(pointSa.base_ring())
2182 84 storres
    if re.search('Interval', pointBaseRingString) is None: # Point
2183 84 storres
        pointSo = pobyso_constant_sa_so(pointSa, pointPrecSa)
2184 84 storres
    else: # Interval.
2185 84 storres
        pointSo = pobyso_interval_to_range_sa_so(pointSa, pointPrecSa)
2186 37 storres
    # Sollyafy the function.
2187 159 storres
    functionSo = pobyso_parse_string_sa_so(functionSa._assume_str().replace('_SAGE_VAR_', ''))
2188 37 storres
    if sollya_lib_obj_is_error(functionSo):
2189 37 storres
        print "pobyso_tailorform: function string can't be parsed!"
2190 37 storres
        return None
2191 37 storres
    # Sollyafy the degree
2192 84 storres
    degreeSo = sollya_lib_constant_from_int(int(degreeSa))
2193 37 storres
    # Sollyafy the point
2194 37 storres
    # Call Sollya
2195 83 storres
    taylorFormSo = \
2196 83 storres
        sollya_lib_taylorform(functionSo, degreeSo, pointSo, errorTypeSo,\
2197 37 storres
                                         None)
2198 85 storres
    sollya_lib_clear_obj(functionSo)
2199 85 storres
    sollya_lib_clear_obj(degreeSo)
2200 85 storres
    sollya_lib_clear_obj(pointSo)
2201 85 storres
    sollya_lib_clear_obj(errorTypeSo)
2202 38 storres
    (tfsAsList, numElements, isEndElliptic) = \
2203 38 storres
            pobyso_get_list_elements_so_so(taylorFormSo)
2204 37 storres
    polySo = tfsAsList[0]
2205 38 storres
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
2206 37 storres
    polyRealField = RealField(maxPrecision)
2207 38 storres
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, polyRealField)
2208 226 storres
    if sollyaPrecisionChangedSa:
2209 226 storres
        sollya_lib_set_prec(initialSollyaPrecSo)
2210 226 storres
    sollya_lib_clear_obj(initialSollyaPrecSo)
2211 37 storres
    polynomialRing = polyRealField[str(varSa)]
2212 37 storres
    polySa = polynomial(expSa, polynomialRing)
2213 37 storres
    taylorFormSa = [polySa]
2214 85 storres
    # Final clean-up.
2215 85 storres
    sollya_lib_clear_obj(taylorFormSo)
2216 51 storres
    return(taylorFormSa)
2217 51 storres
# End pobyso_taylor_form_sa_sa
2218 54 storres
2219 54 storres
def pobyso_taylorform_so_so(functionSo, degreeSo, pointSo, intervalSo=None, \
2220 54 storres
                            errorTypeSo=None):
2221 54 storres
    createdErrorType = False
2222 51 storres
    if errorTypeSo is None:
2223 51 storres
        errorTypeSo = sollya_lib_absolute()
2224 54 storres
        createdErrorType = True
2225 51 storres
    else:
2226 51 storres
        #TODO: deal with the other case.
2227 51 storres
        pass
2228 51 storres
    if intervalSo is None:
2229 54 storres
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
2230 54 storres
                                         errorTypeSo, None)
2231 51 storres
    else:
2232 54 storres
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
2233 54 storres
                                         intervalSo, errorTypeSo, None)
2234 54 storres
    if createdErrorType:
2235 54 storres
        sollya_lib_clear_obj(errorTypeSo)
2236 215 storres
    return resultSo
2237 51 storres
2238 37 storres
2239 37 storres
def pobyso_univar_polynomial_print_reverse(polySa):
2240 51 storres
    """ Legacy function. See pobyso_univar_polynomial_print_reverse_sa_sa. """
2241 51 storres
    return(pobyso_univar_polynomial_print_reverse_sa_sa(polySa))
2242 38 storres
2243 51 storres
def pobyso_univar_polynomial_print_reverse_sa_sa(polySa):
2244 37 storres
    """
2245 37 storres
    Return the string representation of a univariate polynomial with
2246 38 storres
    monomials ordered in the x^0..x^n order of the monomials.
2247 37 storres
    Remember: Sage
2248 37 storres
    """
2249 37 storres
    polynomialRing = polySa.base_ring()
2250 37 storres
    # A very expensive solution:
2251 37 storres
    # -create a fake multivariate polynomial field with only one variable,
2252 37 storres
    #   specifying a negative lexicographical order;
2253 37 storres
    mpolynomialRing = PolynomialRing(polynomialRing.base(), \
2254 37 storres
                                     polynomialRing.variable_name(), \
2255 37 storres
                                     1, order='neglex')
2256 37 storres
    # - convert the univariate argument polynomial into a multivariate
2257 37 storres
    #   version;
2258 37 storres
    p = mpolynomialRing(polySa)
2259 37 storres
    # - return the string representation of the converted form.
2260 37 storres
    # There is no simple str() method defined for p's class.
2261 37 storres
    return(p.__str__())
2262 5 storres
#
2263 230 storres
#print pobyso_get_prec()
2264 5 storres
pobyso_set_prec(165)
2265 230 storres
#print pobyso_get_prec()
2266 230 storres
#a=100
2267 230 storres
#print type(a)
2268 230 storres
#id(a)
2269 230 storres
#print "Max arity: ", pobyso_max_arity
2270 230 storres
#print "Function tripleDouble (43) as a string: ", pobyso_function_type_as_string(43)
2271 230 storres
#print "Function None (44) as a string: ", pobyso_function_type_as_string(44)
2272 246 storres
sys.stderr.write("\t...Pobyso check done.\n")