Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 252

Historique | Voir | Annoter | Télécharger (86,15 ko)

1 5 storres
"""
2 209 storres
@file pobyso.py
3 5 storres
Actual functions to use in Sage
4 5 storres
ST 2012-11-13
5 5 storres

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

9 38 storres
pobyso functions come in five flavors:
10 57 storres
- the _so_so (arguments and returned objects are pointers to Sollya objects,
11 57 storres
  includes the void function and the no arguments function that return a
12 57 storres
  pointer to a Sollya object);
13 38 storres
- the _so_sa (argument are pointers to Sollya objects, returned objects are
14 57 storres
  Sage/Python objects or, more generally, information is transfered from the
15 57 storres
  Sollya world to Sage/Python world; e.g. functions without arguments that
16 57 storres
  return a Sage/Python object);
17 38 storres
- the _sa_so (arguments are Sage/Python objects, returned objects are
18 38 storres
  pointers to Sollya objects);
19 38 storres
- the sa_sa (arguments and returned objects are all Sage/Python objects);
20 51 storres
- a catch all flavor, without any suffix, (e. g. functions that have no argument
21 51 storres
  nor return value).
22 57 storres
This classification is not always very strict. Conversion functions from Sollya
23 57 storres
to Sage/Python are sometimes decorated with Sage/Python arguments to set
24 57 storres
the precision. These functions remain in the so_sa category.
25 5 storres
NOTES:
26 5 storres
Reported errors in Eclipse come from the calls to
27 5 storres
the Sollya library
28 5 storres

29 10 storres
ToDo (among other things):
30 10 storres
 -memory management.
31 5 storres
"""
32 5 storres
from ctypes import *
33 37 storres
import re
34 37 storres
from sage.symbolic.expression_conversions import polynomial
35 59 storres
from sage.symbolic.expression_conversions import PolynomialConverter
36 38 storres
"""
37 38 storres
Create the equivalent to an enum for the Sollya function types.
38 38 storres
"""
39 5 storres
(SOLLYA_BASE_FUNC_ABS,
40 5 storres
SOLLYA_BASE_FUNC_ACOS,
41 5 storres
    SOLLYA_BASE_FUNC_ACOSH,
42 5 storres
    SOLLYA_BASE_FUNC_ADD,
43 5 storres
    SOLLYA_BASE_FUNC_ASIN,
44 5 storres
    SOLLYA_BASE_FUNC_ASINH,
45 5 storres
    SOLLYA_BASE_FUNC_ATAN,
46 5 storres
    SOLLYA_BASE_FUNC_ATANH,
47 5 storres
    SOLLYA_BASE_FUNC_CEIL,
48 5 storres
    SOLLYA_BASE_FUNC_CONSTANT,
49 5 storres
    SOLLYA_BASE_FUNC_COS,
50 5 storres
    SOLLYA_BASE_FUNC_COSH,
51 5 storres
    SOLLYA_BASE_FUNC_DIV,
52 5 storres
    SOLLYA_BASE_FUNC_DOUBLE,
53 5 storres
    SOLLYA_BASE_FUNC_DOUBLEDOUBLE,
54 5 storres
    SOLLYA_BASE_FUNC_DOUBLEEXTENDED,
55 5 storres
    SOLLYA_BASE_FUNC_ERF,
56 5 storres
    SOLLYA_BASE_FUNC_ERFC,
57 5 storres
    SOLLYA_BASE_FUNC_EXP,
58 5 storres
    SOLLYA_BASE_FUNC_EXP_M1,
59 5 storres
    SOLLYA_BASE_FUNC_FLOOR,
60 5 storres
    SOLLYA_BASE_FUNC_FREE_VARIABLE,
61 5 storres
    SOLLYA_BASE_FUNC_HALFPRECISION,
62 5 storres
    SOLLYA_BASE_FUNC_LIBRARYCONSTANT,
63 5 storres
    SOLLYA_BASE_FUNC_LIBRARYFUNCTION,
64 5 storres
    SOLLYA_BASE_FUNC_LOG,
65 5 storres
    SOLLYA_BASE_FUNC_LOG_10,
66 5 storres
    SOLLYA_BASE_FUNC_LOG_1P,
67 5 storres
    SOLLYA_BASE_FUNC_LOG_2,
68 5 storres
    SOLLYA_BASE_FUNC_MUL,
69 5 storres
    SOLLYA_BASE_FUNC_NEARESTINT,
70 5 storres
    SOLLYA_BASE_FUNC_NEG,
71 5 storres
    SOLLYA_BASE_FUNC_PI,
72 5 storres
    SOLLYA_BASE_FUNC_POW,
73 5 storres
    SOLLYA_BASE_FUNC_PROCEDUREFUNCTION,
74 5 storres
    SOLLYA_BASE_FUNC_QUAD,
75 5 storres
    SOLLYA_BASE_FUNC_SIN,
76 5 storres
    SOLLYA_BASE_FUNC_SINGLE,
77 5 storres
    SOLLYA_BASE_FUNC_SINH,
78 5 storres
    SOLLYA_BASE_FUNC_SQRT,
79 5 storres
    SOLLYA_BASE_FUNC_SUB,
80 5 storres
    SOLLYA_BASE_FUNC_TAN,
81 5 storres
    SOLLYA_BASE_FUNC_TANH,
82 5 storres
SOLLYA_BASE_FUNC_TRIPLEDOUBLE) = map(int,xrange(44))
83 246 storres
sys.stderr.write("Superficial pobyso check...\n")
84 230 storres
#print "First constant - SOLLYA_BASE_FUNC_ABS: ", SOLLYA_BASE_FUNC_ABS
85 230 storres
#print "Last constant  - SOLLYA_BASE_FUNC_TRIPLEDOUBLE: ", SOLLYA_BASE_FUNC_TRIPLEDOUBLE
86 5 storres
87 5 storres
pobyso_max_arity = 9
88 5 storres
89 85 storres
def pobyso_absolute_so_so():
90 85 storres
    return(sollya_lib_absolute(None))
91 85 storres
92 5 storres
def pobyso_autoprint(arg):
93 218 storres
    sollya_lib_autoprint(arg, None)
94 5 storres
95 38 storres
def pobyso_autoprint_so_so(arg):
96 38 storres
    sollya_lib_autoprint(arg,None)
97 54 storres
98 84 storres
def pobyso_bounds_to_range_sa_so(rnLowerBoundSa, rnUpperBoundSa, \
99 84 storres
                                 precisionSa=None):
100 84 storres
    """
101 84 storres
    Return a Sollya range from to 2 RealField Sage elements.
102 84 storres
    The Sollya range element has a sufficient precision to hold all
103 162 storres
    the digits of the widest of the Sage bounds.
104 84 storres
    """
105 84 storres
    # Sanity check.
106 84 storres
    if rnLowerBoundSa > rnUpperBoundSa:
107 84 storres
        return None
108 115 storres
    # Precision stuff.
109 85 storres
    if precisionSa is None:
110 84 storres
        # Check for the largest precision.
111 84 storres
        lbPrecSa = rnLowerBoundSa.parent().precision()
112 84 storres
        ubPrecSa = rnLowerBoundSa.parent().precision()
113 84 storres
        maxPrecSa = max(lbPrecSa, ubPrecSa)
114 84 storres
    else:
115 84 storres
        maxPrecSa = precisionSa
116 115 storres
    # From Sage to Sollya bounds.
117 162 storres
#    lowerBoundSo = sollya_lib_constant(get_rn_value(rnLowerBoundSa),
118 162 storres
#                                       maxPrecSa)
119 162 storres
    lowerBoundSo = pobyso_constant_sa_so(rnLowerBoundSa,
120 162 storres
                                         maxPrecSa)
121 162 storres
    upperBoundSo = pobyso_constant_sa_so(rnUpperBoundSa,
122 154 storres
                                       maxPrecSa)
123 162 storres
124 115 storres
    # From Sollya bounds to range.
125 84 storres
    rangeSo = sollya_lib_range(lowerBoundSo, upperBoundSo)
126 84 storres
    # Back to original precision.
127 84 storres
    # Clean up
128 84 storres
    sollya_lib_clear_obj(lowerBoundSo)
129 84 storres
    sollya_lib_clear_obj(upperBoundSo)
130 154 storres
    return rangeSo
131 84 storres
# End pobyso_bounds_to_range_sa_so
132 84 storres
133 215 storres
def pobyso_build_end_elliptic_list_so_so(*args):
134 215 storres
    """
135 237 storres
    From Sollya argument objects, create a Sollya end elliptic list.
136 237 storres
    Elements of the list are "eaten" (should not be cleared individually,
137 215 storres
    are cleared when the list is cleared).
138 215 storres
    """
139 215 storres
    if len(args) == 0:
140 215 storres
        ## Called with an empty list produced "error".
141 215 storres
        return sollya_lib_build_end_elliptic_list(None)
142 215 storres
    index = 0
143 215 storres
    ## One can not append elements to an elliptic list, prepend only is
144 215 storres
    #  permitted.
145 215 storres
    for argument in reversed(args):
146 215 storres
        if index == 0:
147 215 storres
            listSo = sollya_lib_build_end_elliptic_list(argument, None)
148 215 storres
        else:
149 215 storres
            listSo = sollya_lib_prepend(argument, listSo)
150 215 storres
        index += 1
151 215 storres
    return listSo
152 215 storres
153 215 storres
# End pobyso_build_end_elliptic_list_so_so
154 215 storres
155 54 storres
def pobyso_build_function_sub_so_so(exp1So, exp2So):
156 237 storres
    return sollya_lib_build_function_sub(exp1So, exp2So)
157 54 storres
158 237 storres
def pobyso_build_list_of_ints_sa_so(*args):
159 237 storres
    """
160 237 storres
    Build a Sollya list from Sage integral arguments.
161 237 storres
    """
162 237 storres
    if len(args) == 0:
163 237 storres
        return pobyso_build_list_so_so()
164 237 storres
    ## Make a Sage list of integral constants.
165 237 storres
    intsList = []
166 237 storres
    for intElem in args:
167 237 storres
        intsList.append(pobyso_constant_from_int_sa_so(intElem))
168 237 storres
    return pobyso_build_list_so_so(*intsList)
169 237 storres
170 237 storres
def pobyso_build_list_so_so(*args):
171 237 storres
    """
172 237 storres
    Make a Sollya list out of Sollya objects passed as arguments.
173 237 storres
    If one wants to call it with a list argument, should prepend a "*"
174 237 storres
    before the list variable name.
175 237 storres
    Elements of the list are "eaten" (should not be cleared individually,
176 237 storres
    are cleared when the list is cleared).
177 237 storres
    """
178 237 storres
    if len(args) == 0:
179 237 storres
        ## Called with an empty list produced "error".
180 237 storres
        return sollya_lib_build_list(None)
181 237 storres
    index = 0
182 237 storres
    ## Append the Sollya elements one by one.
183 237 storres
    for elementSo in args:
184 237 storres
        if index == 0:
185 237 storres
            listSo = sollya_lib_build_list(elementSo, None)
186 237 storres
        else:
187 237 storres
            listSo = sollya_lib_append(listSo, elementSo)
188 237 storres
        index += 1
189 237 storres
    return listSo
190 237 storres
# End pobyso_build list_so_so
191 237 storres
192 237 storres
193 85 storres
def pobyso_change_var_in_function_so_so(funcSo, chvarExpSo):
194 54 storres
    """
195 85 storres
    Variable change in a function.
196 85 storres
    """
197 85 storres
    return(sollya_lib_evaluate(funcSo,chvarExpSo))
198 85 storres
# End pobyso_change_var_in_function_so_so
199 85 storres
200 85 storres
def pobyso_chebyshevform_so_so(functionSo, degreeSo, intervalSo):
201 85 storres
    resultSo = sollya_lib_chebyshevform(functionSo, degreeSo, intervalSo)
202 85 storres
    return(resultSo)
203 85 storres
# End pobyso_chebyshevform_so_so.
204 85 storres
205 237 storres
def pobyso_clear_obj(objSo):
206 237 storres
    """
207 237 storres
    Free a Sollya object's memory.
208 237 storres
    Very thin wrapper around sollya_lib_clear_obj().
209 237 storres
    """
210 237 storres
    sollya_lib_clear_obj(objSo)
211 237 storres
# End pobyso_clear_obj.
212 237 storres
213 117 storres
def pobyso_clear_taylorform_sa_so(taylorFormSaSo):
214 117 storres
    """
215 117 storres
    This method is necessary to correctly clean up the memory from Taylor forms.
216 117 storres
    These are made of a Sollya object, a Sollya object list, a Sollya object.
217 117 storres
    For no clearly understood reason, sollya_lib_clear_object_list crashed
218 117 storres
    when applied to the object list.
219 117 storres
    Here, we decompose it into Sage list of Sollya objects references and we
220 237 storres
     clear them one at a time.
221 117 storres
    """
222 117 storres
    sollya_lib_clear_obj(taylorFormSaSo[0])
223 117 storres
    (coefficientsErrorsListSaSo, numElementsSa, isEndEllipticSa) = \
224 117 storres
        pobyso_get_list_elements_so_so(taylorFormSaSo[1])
225 117 storres
    for element in coefficientsErrorsListSaSo:
226 117 storres
        sollya_lib_clear_obj(element)
227 117 storres
    sollya_lib_clear_obj(taylorFormSaSo[1])
228 117 storres
    sollya_lib_clear_obj(taylorFormSaSo[2])
229 117 storres
# End pobyso_clear_taylorform_sa_so
230 117 storres
231 85 storres
def pobyso_cmp(rnArgSa, cteSo):
232 85 storres
    """
233 237 storres
    Deprecated, use pobyso_cmp_sa_so_sa instead.
234 237 storres
    """
235 237 storres
    print "Deprecated, use pobyso_cmp_sa_so_sa instead."
236 237 storres
    return pobyso_cmp_sa_so_sa(rnArgSa, cteSo)
237 237 storres
# End  pobyso_cmp
238 237 storres
239 237 storres
def pobyso_cmp_sa_so_sa(rnArgSa, cteSo):
240 237 storres
    """
241 54 storres
    Compare the MPFR value a RealNumber with that of a Sollya constant.
242 54 storres

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

799 118 storres
    INPUT:
800 118 storres
    - objectListSo: a Sollya list of Sollya objects.
801 118 storres

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

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

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