Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 279

Historique | Voir | Annoter | Télécharger (89,58 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 262 storres
def pobyso_evaluate_so_sa(funcSo, argumentSo):
416 262 storres
    """
417 262 storres
    Evaluates funcSo for arguemntSo through sollya_lib_evaluate() and return
418 262 storres
    the result as a Sage object
419 262 storres
    """
420 262 storres
    evalSo = sollya_lib_evaluate(funcSo, argumentSo)
421 262 storres
    if pobyso_is_error_so_sa(evalSo):
422 262 storres
        return None
423 262 storres
    if pobyso_is_range_so_sa(evalSo):
424 262 storres
        retVal = pobyso_range_to_interval_so_sa(evalSo)
425 262 storres
    else:
426 262 storres
        retVal = pobyso_get_constant_as_rn(evalSo)
427 262 storres
    sollya_lib_clear_obj(evalSo)
428 262 storres
    return retVal
429 262 storres
# End pobyso_evaluate_so_sa.
430 262 storres
431 215 storres
def pobyso_evaluate_so_so(funcSo, argumentSo):
432 215 storres
    """
433 215 storres
    Evaluates funcSo for arguemntSo through sollya_lib_evaluate().
434 215 storres
    """
435 215 storres
    return sollya_lib_evaluate(funcSo, argumentSo)
436 215 storres
# End pobyso_evaluate_so_so.
437 240 storres
#
438 240 storres
def pobyso_diff_so_so(funcSo):
439 240 storres
    """
440 240 storres
    Very thin wrapper around sollya_lib_diff.
441 240 storres
    """
442 240 storres
    ## TODO: add a check to make sure funcSo is a functional expression.
443 240 storres
    return sollya_lib_diff(funcSo)
444 237 storres
445 234 storres
def pobyso_dirty_find_zeros_so_so(funcSo, rangeSo):
446 234 storres
    """
447 234 storres
    Thin wrapper over sollya_lib_dirtyfindzeros()
448 234 storres
    """
449 234 storres
    return sollya_lib_dirtyfindzeros(funcSo, rangeSo)
450 234 storres
# End pobys_dirty_find_zeros
451 215 storres
452 237 storres
def pobyso_dirty_inf_norm_so_so(funcSo, rangeSo, preSo=None):
453 237 storres
    """
454 237 storres
    Thin wrapper around sollya_dirtyinfnorm().
455 237 storres
    """
456 237 storres
    # TODO: manage the precision change.
457 237 storres
458 237 storres
    return sollya_lib_dirtyinfnorm(funcSo, rangeSo)
459 237 storres
# End pobyso_dirty_inf_norm_so_so
460 237 storres
461 261 storres
def pobyso_find_zeros_so_so(funcSo, rangeSo):
462 261 storres
    """
463 261 storres
    Thin wrapper over sollya_lib_findzeros()
464 261 storres
    """
465 261 storres
    return sollya_lib_findzeros(funcSo, rangeSo)
466 261 storres
# End pobys_find_zeros
467 261 storres
468 234 storres
def pobyso_float_list_so_sa(listSo):
469 234 storres
    """
470 240 storres
    Return a Sollya list of floating-point numbers as a Sage list of
471 234 storres
    floating-point numbers.
472 240 storres
    TODO: add a test to make sure that each element of the list is a constant.
473 234 storres
    """
474 234 storres
    listSa   = []
475 239 storres
    ## The function returns none if the list is empty or an error has happened.
476 239 storres
    retVal = pobyso_get_list_elements_so_so(listSo)
477 239 storres
    if retVal is None:
478 239 storres
        return listSa
479 240 storres
    ## Just in case the interface is changed and an empty list is returned
480 240 storres
    #  instead of None.
481 239 storres
    elif len(retVal) == 0:
482 239 storres
        return listSa
483 239 storres
    else:
484 239 storres
        ## Remember pobyso_get_list_elements_so_so returns more information
485 240 storres
        #  than just the elements of the list (# elements, is_elliptic)
486 239 storres
        listSaSo, numElements, isEndElliptic = retVal
487 240 storres
    ## Return an empty list.
488 239 storres
    if numElements == 0:
489 239 storres
        return listSa
490 234 storres
    ## Search first for the maximum precision of the elements
491 234 storres
    maxPrecSa = 0
492 234 storres
    for floatSo in listSaSo:
493 240 storres
        #pobyso_autoprint(floatSo)
494 234 storres
        curPrecSa =  pobyso_get_prec_of_constant_so_sa(floatSo)
495 234 storres
        if curPrecSa > maxPrecSa:
496 234 storres
            maxPrecSa = curPrecSa
497 234 storres
    ##
498 234 storres
    RF = RealField(maxPrecSa)
499 234 storres
    ##
500 234 storres
    for floatSo in listSaSo:
501 234 storres
        listSa.append(pobyso_get_constant_as_rn_with_rf_so_sa(floatSo))
502 234 storres
    return listSa
503 234 storres
# End pobyso_float_list_so_sa
504 234 storres
505 209 storres
def pobyso_float_poly_sa_so(polySa, precSa = None):
506 209 storres
    """
507 209 storres
    Create a Sollya polynomial from a Sage RealField polynomial.
508 209 storres
    """
509 209 storres
    ## TODO: filter arguments.
510 209 storres
    ## Precision. If a precision is given, convert the polynomial
511 209 storres
    #  into the right polynomial field. If not convert it straight
512 209 storres
    #  to Sollya.
513 218 storres
    sollyaPrecChanged = False
514 226 storres
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
515 218 storres
    if precSa is None:
516 209 storres
        precSa = polySa.parent().base_ring().precision()
517 226 storres
    if (precSa > initialSollyaPrecSa):
518 228 storres
        assert precSa >= 2, "Precision change <2 requested"
519 226 storres
        if precSa <= 2:
520 226 storres
            print inspect.stack()[0][3], ": precision change <= 2 requested"
521 218 storres
        precSo = pobyso_constant_from_int(precSa)
522 218 storres
        pobyso_set_prec_so_so(precSo)
523 218 storres
        sollya_lib_clear_obj(precSo)
524 272 storres
        sollyaPrecChanged = True
525 272 storres
    ## Free variable stuff.
526 272 storres
    freeVariableNameChanged = False
527 272 storres
    polyFreeVariableNameSa = \
528 272 storres
        str(polySa.variables()[0])
529 272 storres
    currentFreeVariableNameSa = pobyso_get_free_variable_name_so_sa()
530 272 storres
    if polyFreeVariableNameSa != currentFreeVariableNameSa:
531 272 storres
        #print "Free variable names do not match.", polyFreeVariableNameSa
532 272 storres
        sollya_lib_name_free_variable(polyFreeVariableNameSa)
533 272 storres
        freeVariableNameChanged = True
534 209 storres
    ## Get exponents and coefficients.
535 218 storres
    exponentsSa     = polySa.exponents()
536 218 storres
    coefficientsSa  = polySa.coefficients()
537 209 storres
    ## Build the polynomial.
538 209 storres
    polySo = None
539 213 storres
    for coefficientSa, exponentSa in zip(coefficientsSa, exponentsSa):
540 209 storres
        #print coefficientSa.n(prec=precSa), exponentSa
541 209 storres
        coefficientSo = \
542 209 storres
            pobyso_constant_sa_so(coefficientSa)
543 209 storres
        #pobyso_autoprint(coefficientSo)
544 209 storres
        exponentSo = \
545 209 storres
            pobyso_constant_from_int_sa_so(exponentSa)
546 209 storres
        #pobyso_autoprint(exponentSo)
547 209 storres
        monomialSo = sollya_lib_build_function_pow(
548 209 storres
                       sollya_lib_build_function_free_variable(),
549 209 storres
                       exponentSo)
550 218 storres
        polyTermSo = sollya_lib_build_function_mul(coefficientSo,
551 218 storres
                                                       monomialSo)
552 209 storres
        if polySo is None:
553 218 storres
            polySo = polyTermSo
554 209 storres
        else:
555 209 storres
            polySo = sollya_lib_build_function_add(polySo, polyTermSo)
556 218 storres
    if sollyaPrecChanged:
557 226 storres
        pobyso_set_prec_so_so(initialSollyaPrecSo)
558 272 storres
    sollya_lib_clear_obj(initialSollyaPrecSo)
559 272 storres
    ## Do not set back the free variable name in Sollya to its initial value:
560 272 storres
    #  it will change it back, in the Sollya polynomial, to what it was in the
561 272 storres
    #  first place.
562 209 storres
    return polySo
563 209 storres
# End pobyso_float_poly_sa_so
564 209 storres
565 209 storres
def pobyso_float_poly_so_sa(polySo, realFieldSa=None):
566 209 storres
    """
567 209 storres
    Convert a Sollya polynomial into a Sage floating-point polynomial.
568 209 storres
    If no realField is given, a RealField corresponding to the maximum
569 209 storres
    precision of the coefficients is internally computed.
570 209 storres
    The real field is not returned but can be easily retrieved from
571 209 storres
    the polynomial itself.
572 209 storres
    ALGORITHM:
573 209 storres
    - (optional) compute the RealField of the coefficients;
574 209 storres
    - convert the Sollya expression into a Sage expression;
575 209 storres
    - convert the Sage expression into a Sage polynomial
576 209 storres
    """
577 209 storres
    if realFieldSa is None:
578 209 storres
        expressionPrecSa = pobyso_get_max_prec_of_exp_so_sa(polySo)
579 218 storres
        #print "Maximum precision of Sollya polynomial coefficients:", expressionPrecSa
580 226 storres
        if expressionPrecSa < 2 or expressionPrecSa > 2147483391:
581 226 storres
            print "Maximum degree of expression:", expressionPrecSa
582 209 storres
        realFieldSa      = RealField(expressionPrecSa)
583 209 storres
    #print "Sollya expression before...",
584 209 storres
    #pobyso_autoprint(polySo)
585 209 storres
586 209 storres
    expressionSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo,
587 209 storres
                                                             realFieldSa)
588 218 storres
    #print "...Sollya expression after."
589 209 storres
    #pobyso_autoprint(polySo)
590 209 storres
    polyVariableSa = expressionSa.variables()[0]
591 209 storres
    polyRingSa     = realFieldSa[str(polyVariableSa)]
592 209 storres
    #print polyRingSa
593 209 storres
    # Do not use the polynomial(expressionSa, ring=polyRingSa) form!
594 209 storres
    polynomialSa = polyRingSa(expressionSa)
595 215 storres
    polyCoeffsListSa = polynomialSa.coefficients()
596 215 storres
    #for coeff in polyCoeffsListSa:
597 215 storres
    #    print coeff.abs().n()
598 209 storres
    return polynomialSa
599 209 storres
# End pobyso_float_poly_so_sa
600 209 storres
601 215 storres
def pobyso_free_variable():
602 215 storres
    """
603 215 storres
    Ultra thin wrapper around the sollya_lib_function_build_free_variable function.
604 215 storres
    """
605 215 storres
    return sollya_lib_build_function_free_variable()
606 209 storres
607 5 storres
def pobyso_function_type_as_string(funcType):
608 38 storres
    """ Legacy function. See pobyso_function_type_as_string_so_sa. """
609 38 storres
    return(pobyso_function_type_as_string_so_sa(funcType))
610 38 storres
611 38 storres
def pobyso_function_type_as_string_so_sa(funcType):
612 38 storres
    """
613 38 storres
    Numeric Sollya function codes -> Sage mathematical function names.
614 38 storres
    Notice that pow -> ^ (a la Sage, not a la Python).
615 38 storres
    """
616 5 storres
    if funcType == SOLLYA_BASE_FUNC_ABS:
617 5 storres
        return "abs"
618 5 storres
    elif funcType == SOLLYA_BASE_FUNC_ACOS:
619 5 storres
        return "arccos"
620 5 storres
    elif funcType == SOLLYA_BASE_FUNC_ACOSH:
621 5 storres
        return "arccosh"
622 5 storres
    elif funcType == SOLLYA_BASE_FUNC_ADD:
623 5 storres
        return "+"
624 5 storres
    elif funcType == SOLLYA_BASE_FUNC_ASIN:
625 5 storres
        return "arcsin"
626 5 storres
    elif funcType == SOLLYA_BASE_FUNC_ASINH:
627 5 storres
        return "arcsinh"
628 5 storres
    elif funcType == SOLLYA_BASE_FUNC_ATAN:
629 5 storres
        return "arctan"
630 5 storres
    elif funcType == SOLLYA_BASE_FUNC_ATANH:
631 5 storres
        return "arctanh"
632 5 storres
    elif funcType == SOLLYA_BASE_FUNC_CEIL:
633 5 storres
        return "ceil"
634 5 storres
    elif funcType == SOLLYA_BASE_FUNC_CONSTANT:
635 5 storres
        return "cte"
636 5 storres
    elif funcType == SOLLYA_BASE_FUNC_COS:
637 5 storres
        return "cos"
638 5 storres
    elif funcType == SOLLYA_BASE_FUNC_COSH:
639 5 storres
        return "cosh"
640 5 storres
    elif funcType == SOLLYA_BASE_FUNC_DIV:
641 5 storres
        return "/"
642 5 storres
    elif funcType == SOLLYA_BASE_FUNC_DOUBLE:
643 5 storres
        return "double"
644 5 storres
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEDOUBLE:
645 5 storres
        return "doubleDouble"
646 5 storres
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEEXTENDED:
647 5 storres
        return "doubleDxtended"
648 5 storres
    elif funcType == SOLLYA_BASE_FUNC_ERF:
649 5 storres
        return "erf"
650 5 storres
    elif funcType == SOLLYA_BASE_FUNC_ERFC:
651 5 storres
        return "erfc"
652 5 storres
    elif funcType == SOLLYA_BASE_FUNC_EXP:
653 5 storres
        return "exp"
654 5 storres
    elif funcType == SOLLYA_BASE_FUNC_EXP_M1:
655 5 storres
        return "expm1"
656 5 storres
    elif funcType == SOLLYA_BASE_FUNC_FLOOR:
657 5 storres
        return "floor"
658 5 storres
    elif funcType == SOLLYA_BASE_FUNC_FREE_VARIABLE:
659 5 storres
        return "freeVariable"
660 5 storres
    elif funcType == SOLLYA_BASE_FUNC_HALFPRECISION:
661 5 storres
        return "halfPrecision"
662 5 storres
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYCONSTANT:
663 5 storres
        return "libraryConstant"
664 5 storres
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYFUNCTION:
665 5 storres
        return "libraryFunction"
666 5 storres
    elif funcType == SOLLYA_BASE_FUNC_LOG:
667 5 storres
        return "log"
668 5 storres
    elif funcType == SOLLYA_BASE_FUNC_LOG_10:
669 5 storres
        return "log10"
670 5 storres
    elif funcType == SOLLYA_BASE_FUNC_LOG_1P:
671 5 storres
        return "log1p"
672 5 storres
    elif funcType == SOLLYA_BASE_FUNC_LOG_2:
673 5 storres
        return "log2"
674 5 storres
    elif funcType == SOLLYA_BASE_FUNC_MUL:
675 5 storres
        return "*"
676 5 storres
    elif funcType == SOLLYA_BASE_FUNC_NEARESTINT:
677 5 storres
        return "round"
678 5 storres
    elif funcType == SOLLYA_BASE_FUNC_NEG:
679 5 storres
        return "__neg__"
680 5 storres
    elif funcType == SOLLYA_BASE_FUNC_PI:
681 5 storres
        return "pi"
682 5 storres
    elif funcType == SOLLYA_BASE_FUNC_POW:
683 5 storres
        return "^"
684 5 storres
    elif funcType == SOLLYA_BASE_FUNC_PROCEDUREFUNCTION:
685 5 storres
        return "procedureFunction"
686 5 storres
    elif funcType == SOLLYA_BASE_FUNC_QUAD:
687 5 storres
        return "quad"
688 5 storres
    elif funcType == SOLLYA_BASE_FUNC_SIN:
689 5 storres
        return "sin"
690 5 storres
    elif funcType == SOLLYA_BASE_FUNC_SINGLE:
691 5 storres
        return "single"
692 5 storres
    elif funcType == SOLLYA_BASE_FUNC_SINH:
693 5 storres
        return "sinh"
694 5 storres
    elif funcType == SOLLYA_BASE_FUNC_SQRT:
695 5 storres
        return "sqrt"
696 5 storres
    elif funcType == SOLLYA_BASE_FUNC_SUB:
697 5 storres
        return "-"
698 5 storres
    elif funcType == SOLLYA_BASE_FUNC_TAN:
699 5 storres
        return "tan"
700 5 storres
    elif funcType == SOLLYA_BASE_FUNC_TANH:
701 5 storres
        return "tanh"
702 5 storres
    elif funcType == SOLLYA_BASE_FUNC_TRIPLEDOUBLE:
703 5 storres
        return "tripleDouble"
704 5 storres
    else:
705 5 storres
        return None
706 5 storres
707 85 storres
def pobyso_get_constant(rnArgSa, constSo):
708 38 storres
    """ Legacy function. See pobyso_get_constant_so_sa. """
709 209 storres
    return pobyso_get_constant_so_sa(rnArgSa, constSo)
710 209 storres
# End pobyso_get_constant
711 38 storres
712 84 storres
def pobyso_get_constant_so_sa(rnArgSa, constSo):
713 52 storres
    """
714 85 storres
    Set the value of rnArgSo to the value of constSo in MPFR_RNDN mode.
715 52 storres
    rnArg must already exist and belong to some RealField.
716 85 storres
    We assume that constSo points to a Sollya constant.
717 52 storres
    """
718 209 storres
    outcome = sollya_lib_get_constant(get_rn_value(rnArgSa), constSo)
719 209 storres
    if outcome == 0: # Failure because constSo is not a constant expression.
720 209 storres
        return None
721 209 storres
    else:
722 209 storres
        return outcome
723 209 storres
# End  pobyso_get_constant_so_sa
724 209 storres
725 57 storres
def pobyso_get_constant_as_rn(ctExpSo):
726 83 storres
    """
727 83 storres
    Legacy function. See pobyso_get_constant_as_rn_so_sa.
728 83 storres
    """
729 57 storres
    return(pobyso_get_constant_as_rn_so_sa(ctExpSo))
730 38 storres
731 56 storres
def pobyso_get_constant_as_rn_so_sa(constExpSo):
732 83 storres
    """
733 83 storres
    Get a Sollya constant as a Sage "real number".
734 83 storres
    The precision of the floating-point number returned is that of the Sollya
735 83 storres
    constant.
736 83 storres
    """
737 218 storres
    #print "Before computing precision of variable..."
738 218 storres
    #pobyso_autoprint(constExpSo)
739 209 storres
    precisionSa  = pobyso_get_prec_of_constant_so_sa(constExpSo)
740 218 storres
    #print "precisionSa:", precisionSa
741 209 storres
    ## If the expression can not be exactly converted, None is returned.
742 209 storres
    #  In this case opt for the Sollya current expression.
743 209 storres
    if precisionSa is None:
744 209 storres
        precisionSa = pobyso_get_prec_so_sa()
745 56 storres
    RRRR = RealField(precisionSa)
746 56 storres
    rnSa = RRRR(0)
747 209 storres
    outcome = sollya_lib_get_constant(get_rn_value(rnSa), constExpSo)
748 209 storres
    if outcome == 0:
749 209 storres
        return None
750 209 storres
    else:
751 209 storres
        return rnSa
752 83 storres
# End pobyso_get_constant_as_rn_so_sa
753 38 storres
754 38 storres
def pobyso_get_constant_as_rn_with_rf(ctExp, realField):
755 83 storres
    """
756 83 storres
    Legacy function. See pobyso_get_constant_as_rn_with_rf_so_sa.
757 83 storres
    """
758 209 storres
    return pobyso_get_constant_as_rn_with_rf_so_sa(ctExp, realField)
759 209 storres
# End pobyso_get_constant_as_rn_with_rf
760 5 storres
761 56 storres
def pobyso_get_constant_as_rn_with_rf_so_sa(ctExpSo, realFieldSa = None):
762 83 storres
    """
763 83 storres
    Get a Sollya constant as a Sage "real number".
764 83 storres
    If no real field is specified, the precision of the floating-point number
765 85 storres
    returned is that of the Sollya constant.
766 83 storres
    Otherwise is is that of the real field. Hence rounding may happen.
767 83 storres
    """
768 56 storres
    if realFieldSa is None:
769 209 storres
        return pobyso_get_constant_as_rn_so_sa(ctExpSo)
770 56 storres
    rnSa = realFieldSa(0)
771 209 storres
    outcome = sollya_lib_get_constant(get_rn_value(rnSa), ctExpSo)
772 209 storres
    if outcome == 0:
773 209 storres
        return None
774 209 storres
    else:
775 209 storres
        return rnSa
776 83 storres
# End pobyso_get_constant_as_rn_with_rf_so_sa
777 38 storres
778 5 storres
def pobyso_get_free_variable_name():
779 83 storres
    """
780 83 storres
    Legacy function. See pobyso_get_free_variable_name_so_sa.
781 83 storres
    """
782 38 storres
    return(pobyso_get_free_variable_name_so_sa())
783 38 storres
784 38 storres
def pobyso_get_free_variable_name_so_sa():
785 209 storres
    return sollya_lib_get_free_variable_name()
786 5 storres
787 38 storres
def pobyso_get_function_arity(expressionSo):
788 83 storres
    """
789 83 storres
    Legacy function. See pobyso_get_function_arity_so_sa.
790 83 storres
    """
791 38 storres
    return(pobyso_get_function_arity_so_sa(expressionSo))
792 38 storres
793 38 storres
def pobyso_get_function_arity_so_sa(expressionSo):
794 5 storres
    arity = c_int(0)
795 38 storres
    sollya_lib_get_function_arity(byref(arity),expressionSo)
796 209 storres
    return int(arity.value)
797 5 storres
798 38 storres
def pobyso_get_head_function(expressionSo):
799 83 storres
    """
800 83 storres
    Legacy function. See pobyso_get_head_function_so_sa.
801 83 storres
    """
802 38 storres
    return(pobyso_get_head_function_so_sa(expressionSo))
803 38 storres
804 38 storres
def pobyso_get_head_function_so_sa(expressionSo):
805 5 storres
    functionType = c_int(0)
806 218 storres
    sollya_lib_get_head_function(byref(functionType), expressionSo)
807 209 storres
    return int(functionType.value)
808 5 storres
809 56 storres
def pobyso_get_interval_from_range_so_sa(soRange, realIntervalFieldSa = None ):
810 53 storres
    """
811 53 storres
    Return the Sage interval corresponding to the Sollya range argument.
812 83 storres
    If no reaIntervalField is passed as an argument, the interval bounds are not
813 56 storres
    rounded: they are elements of RealIntervalField of the "right" precision
814 56 storres
    to hold all the digits.
815 53 storres
    """
816 53 storres
    prec = c_int(0)
817 56 storres
    if realIntervalFieldSa is None:
818 56 storres
        retval = sollya_lib_get_prec_of_range(byref(prec), soRange, None)
819 56 storres
        if retval == 0:
820 209 storres
            return None
821 56 storres
        realIntervalFieldSa = RealIntervalField(prec.value)
822 56 storres
    intervalSa = realIntervalFieldSa(0,0)
823 53 storres
    retval = \
824 53 storres
        sollya_lib_get_interval_from_range(get_interval_value(intervalSa),\
825 53 storres
                                           soRange)
826 53 storres
    if retval == 0:
827 209 storres
        return None
828 209 storres
    return intervalSa
829 56 storres
# End pobyso_get_interval_from_range_so_sa
830 56 storres
831 5 storres
def pobyso_get_list_elements(soObj):
832 38 storres
    """ Legacy function. See pobyso_get_list_elements_so_so. """
833 209 storres
    return pobyso_get_list_elements_so_so(soObj)
834 38 storres
835 117 storres
def pobyso_get_list_elements_so_so(objectListSo):
836 51 storres
    """
837 118 storres
    Get the Sollya list elements as a Sage/Python array of Sollya objects.
838 118 storres

839 118 storres
    INPUT:
840 118 storres
    - objectListSo: a Sollya list of Sollya objects.
841 118 storres

842 118 storres
    OUTPUT:
843 118 storres
    - a Sage/Python tuple made of:
844 118 storres
      - a Sage/Python list of Sollya objects,
845 118 storres
      - a Sage/Python int holding the number of elements,
846 118 storres
      - a Sage/Python int stating (!= 0) that the list is end-elliptic.
847 118 storres
    NOTE::
848 118 storres
        We recover the addresses of the Sollya object from the list of pointers
849 118 storres
        returned by sollya_lib_get_list_elements. The list itself is freed.
850 118 storres
    TODO::
851 118 storres
        Figure out what to do with numElements since the number of elements
852 118 storres
        can easily be recovered from the list itself.
853 118 storres
        Ditto for isEndElliptic.
854 51 storres
    """
855 5 storres
    listAddress = POINTER(c_longlong)()
856 5 storres
    numElements = c_int(0)
857 5 storres
    isEndElliptic = c_int(0)
858 117 storres
    listAsSageList = []
859 5 storres
    result = sollya_lib_get_list_elements(byref(listAddress),\
860 54 storres
                                          byref(numElements),\
861 54 storres
                                          byref(isEndElliptic),\
862 117 storres
                                          objectListSo)
863 5 storres
    if result == 0 :
864 5 storres
        return None
865 5 storres
    for i in xrange(0, numElements.value, 1):
866 118 storres
       #listAsSageList.append(sollya_lib_copy_obj(listAddress[i]))
867 118 storres
       listAsSageList.append(listAddress[i])
868 117 storres
       # Clear each of the elements returned by Sollya.
869 118 storres
       #sollya_lib_clear_obj(listAddress[i])
870 117 storres
    # Free the list itself.
871 117 storres
    sollya_lib_free(listAddress)
872 209 storres
    return (listAsSageList, numElements.value, isEndElliptic.value)
873 5 storres
874 38 storres
def pobyso_get_max_prec_of_exp(soExp):
875 38 storres
    """ Legacy function. See pobyso_get_max_prec_of_exp_so_sa. """
876 209 storres
    return pobyso_get_max_prec_of_exp_so_sa(soExp)
877 5 storres
878 85 storres
def pobyso_get_max_prec_of_exp_so_sa(expSo):
879 38 storres
    """
880 38 storres
    Get the maximum precision used for the numbers in a Sollya expression.
881 52 storres

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

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