Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 260

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

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

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

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

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

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

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

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

804 118 storres
    INPUT:
805 118 storres
    - objectListSo: a Sollya list of Sollya objects.
806 118 storres

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

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

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