Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 268

Historique | Voir | Annoter | Télécharger (88,13 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 218 storres
        sollyaPrecChanged = True
525 209 storres
    ## Get exponents and coefficients.
526 218 storres
    exponentsSa     = polySa.exponents()
527 218 storres
    coefficientsSa  = polySa.coefficients()
528 209 storres
    ## Build the polynomial.
529 209 storres
    polySo = None
530 213 storres
    for coefficientSa, exponentSa in zip(coefficientsSa, exponentsSa):
531 209 storres
        #print coefficientSa.n(prec=precSa), exponentSa
532 209 storres
        coefficientSo = \
533 209 storres
            pobyso_constant_sa_so(coefficientSa)
534 209 storres
        #pobyso_autoprint(coefficientSo)
535 209 storres
        exponentSo = \
536 209 storres
            pobyso_constant_from_int_sa_so(exponentSa)
537 209 storres
        #pobyso_autoprint(exponentSo)
538 209 storres
        monomialSo = sollya_lib_build_function_pow(
539 209 storres
                       sollya_lib_build_function_free_variable(),
540 209 storres
                       exponentSo)
541 218 storres
        polyTermSo = sollya_lib_build_function_mul(coefficientSo,
542 218 storres
                                                       monomialSo)
543 209 storres
        if polySo is None:
544 218 storres
            polySo = polyTermSo
545 209 storres
        else:
546 209 storres
            polySo = sollya_lib_build_function_add(polySo, polyTermSo)
547 218 storres
    if sollyaPrecChanged:
548 226 storres
        pobyso_set_prec_so_so(initialSollyaPrecSo)
549 226 storres
        sollya_lib_clear_obj(initialSollyaPrecSo)
550 209 storres
    return polySo
551 209 storres
# End pobyso_float_poly_sa_so
552 209 storres
553 209 storres
def pobyso_float_poly_so_sa(polySo, realFieldSa=None):
554 209 storres
    """
555 209 storres
    Convert a Sollya polynomial into a Sage floating-point polynomial.
556 209 storres
    If no realField is given, a RealField corresponding to the maximum
557 209 storres
    precision of the coefficients is internally computed.
558 209 storres
    The real field is not returned but can be easily retrieved from
559 209 storres
    the polynomial itself.
560 209 storres
    ALGORITHM:
561 209 storres
    - (optional) compute the RealField of the coefficients;
562 209 storres
    - convert the Sollya expression into a Sage expression;
563 209 storres
    - convert the Sage expression into a Sage polynomial
564 209 storres
    """
565 209 storres
    if realFieldSa is None:
566 209 storres
        expressionPrecSa = pobyso_get_max_prec_of_exp_so_sa(polySo)
567 218 storres
        #print "Maximum precision of Sollya polynomial coefficients:", expressionPrecSa
568 226 storres
        if expressionPrecSa < 2 or expressionPrecSa > 2147483391:
569 226 storres
            print "Maximum degree of expression:", expressionPrecSa
570 209 storres
        realFieldSa      = RealField(expressionPrecSa)
571 209 storres
    #print "Sollya expression before...",
572 209 storres
    #pobyso_autoprint(polySo)
573 209 storres
574 209 storres
    expressionSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo,
575 209 storres
                                                             realFieldSa)
576 218 storres
    #print "...Sollya expression after."
577 209 storres
    #pobyso_autoprint(polySo)
578 209 storres
    polyVariableSa = expressionSa.variables()[0]
579 209 storres
    polyRingSa     = realFieldSa[str(polyVariableSa)]
580 209 storres
    #print polyRingSa
581 209 storres
    # Do not use the polynomial(expressionSa, ring=polyRingSa) form!
582 209 storres
    polynomialSa = polyRingSa(expressionSa)
583 215 storres
    polyCoeffsListSa = polynomialSa.coefficients()
584 215 storres
    #for coeff in polyCoeffsListSa:
585 215 storres
    #    print coeff.abs().n()
586 209 storres
    return polynomialSa
587 209 storres
# End pobyso_float_poly_so_sa
588 209 storres
589 215 storres
def pobyso_free_variable():
590 215 storres
    """
591 215 storres
    Ultra thin wrapper around the sollya_lib_function_build_free_variable function.
592 215 storres
    """
593 215 storres
    return sollya_lib_build_function_free_variable()
594 209 storres
595 5 storres
def pobyso_function_type_as_string(funcType):
596 38 storres
    """ Legacy function. See pobyso_function_type_as_string_so_sa. """
597 38 storres
    return(pobyso_function_type_as_string_so_sa(funcType))
598 38 storres
599 38 storres
def pobyso_function_type_as_string_so_sa(funcType):
600 38 storres
    """
601 38 storres
    Numeric Sollya function codes -> Sage mathematical function names.
602 38 storres
    Notice that pow -> ^ (a la Sage, not a la Python).
603 38 storres
    """
604 5 storres
    if funcType == SOLLYA_BASE_FUNC_ABS:
605 5 storres
        return "abs"
606 5 storres
    elif funcType == SOLLYA_BASE_FUNC_ACOS:
607 5 storres
        return "arccos"
608 5 storres
    elif funcType == SOLLYA_BASE_FUNC_ACOSH:
609 5 storres
        return "arccosh"
610 5 storres
    elif funcType == SOLLYA_BASE_FUNC_ADD:
611 5 storres
        return "+"
612 5 storres
    elif funcType == SOLLYA_BASE_FUNC_ASIN:
613 5 storres
        return "arcsin"
614 5 storres
    elif funcType == SOLLYA_BASE_FUNC_ASINH:
615 5 storres
        return "arcsinh"
616 5 storres
    elif funcType == SOLLYA_BASE_FUNC_ATAN:
617 5 storres
        return "arctan"
618 5 storres
    elif funcType == SOLLYA_BASE_FUNC_ATANH:
619 5 storres
        return "arctanh"
620 5 storres
    elif funcType == SOLLYA_BASE_FUNC_CEIL:
621 5 storres
        return "ceil"
622 5 storres
    elif funcType == SOLLYA_BASE_FUNC_CONSTANT:
623 5 storres
        return "cte"
624 5 storres
    elif funcType == SOLLYA_BASE_FUNC_COS:
625 5 storres
        return "cos"
626 5 storres
    elif funcType == SOLLYA_BASE_FUNC_COSH:
627 5 storres
        return "cosh"
628 5 storres
    elif funcType == SOLLYA_BASE_FUNC_DIV:
629 5 storres
        return "/"
630 5 storres
    elif funcType == SOLLYA_BASE_FUNC_DOUBLE:
631 5 storres
        return "double"
632 5 storres
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEDOUBLE:
633 5 storres
        return "doubleDouble"
634 5 storres
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEEXTENDED:
635 5 storres
        return "doubleDxtended"
636 5 storres
    elif funcType == SOLLYA_BASE_FUNC_ERF:
637 5 storres
        return "erf"
638 5 storres
    elif funcType == SOLLYA_BASE_FUNC_ERFC:
639 5 storres
        return "erfc"
640 5 storres
    elif funcType == SOLLYA_BASE_FUNC_EXP:
641 5 storres
        return "exp"
642 5 storres
    elif funcType == SOLLYA_BASE_FUNC_EXP_M1:
643 5 storres
        return "expm1"
644 5 storres
    elif funcType == SOLLYA_BASE_FUNC_FLOOR:
645 5 storres
        return "floor"
646 5 storres
    elif funcType == SOLLYA_BASE_FUNC_FREE_VARIABLE:
647 5 storres
        return "freeVariable"
648 5 storres
    elif funcType == SOLLYA_BASE_FUNC_HALFPRECISION:
649 5 storres
        return "halfPrecision"
650 5 storres
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYCONSTANT:
651 5 storres
        return "libraryConstant"
652 5 storres
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYFUNCTION:
653 5 storres
        return "libraryFunction"
654 5 storres
    elif funcType == SOLLYA_BASE_FUNC_LOG:
655 5 storres
        return "log"
656 5 storres
    elif funcType == SOLLYA_BASE_FUNC_LOG_10:
657 5 storres
        return "log10"
658 5 storres
    elif funcType == SOLLYA_BASE_FUNC_LOG_1P:
659 5 storres
        return "log1p"
660 5 storres
    elif funcType == SOLLYA_BASE_FUNC_LOG_2:
661 5 storres
        return "log2"
662 5 storres
    elif funcType == SOLLYA_BASE_FUNC_MUL:
663 5 storres
        return "*"
664 5 storres
    elif funcType == SOLLYA_BASE_FUNC_NEARESTINT:
665 5 storres
        return "round"
666 5 storres
    elif funcType == SOLLYA_BASE_FUNC_NEG:
667 5 storres
        return "__neg__"
668 5 storres
    elif funcType == SOLLYA_BASE_FUNC_PI:
669 5 storres
        return "pi"
670 5 storres
    elif funcType == SOLLYA_BASE_FUNC_POW:
671 5 storres
        return "^"
672 5 storres
    elif funcType == SOLLYA_BASE_FUNC_PROCEDUREFUNCTION:
673 5 storres
        return "procedureFunction"
674 5 storres
    elif funcType == SOLLYA_BASE_FUNC_QUAD:
675 5 storres
        return "quad"
676 5 storres
    elif funcType == SOLLYA_BASE_FUNC_SIN:
677 5 storres
        return "sin"
678 5 storres
    elif funcType == SOLLYA_BASE_FUNC_SINGLE:
679 5 storres
        return "single"
680 5 storres
    elif funcType == SOLLYA_BASE_FUNC_SINH:
681 5 storres
        return "sinh"
682 5 storres
    elif funcType == SOLLYA_BASE_FUNC_SQRT:
683 5 storres
        return "sqrt"
684 5 storres
    elif funcType == SOLLYA_BASE_FUNC_SUB:
685 5 storres
        return "-"
686 5 storres
    elif funcType == SOLLYA_BASE_FUNC_TAN:
687 5 storres
        return "tan"
688 5 storres
    elif funcType == SOLLYA_BASE_FUNC_TANH:
689 5 storres
        return "tanh"
690 5 storres
    elif funcType == SOLLYA_BASE_FUNC_TRIPLEDOUBLE:
691 5 storres
        return "tripleDouble"
692 5 storres
    else:
693 5 storres
        return None
694 5 storres
695 85 storres
def pobyso_get_constant(rnArgSa, constSo):
696 38 storres
    """ Legacy function. See pobyso_get_constant_so_sa. """
697 209 storres
    return pobyso_get_constant_so_sa(rnArgSa, constSo)
698 209 storres
# End pobyso_get_constant
699 38 storres
700 84 storres
def pobyso_get_constant_so_sa(rnArgSa, constSo):
701 52 storres
    """
702 85 storres
    Set the value of rnArgSo to the value of constSo in MPFR_RNDN mode.
703 52 storres
    rnArg must already exist and belong to some RealField.
704 85 storres
    We assume that constSo points to a Sollya constant.
705 52 storres
    """
706 209 storres
    outcome = sollya_lib_get_constant(get_rn_value(rnArgSa), constSo)
707 209 storres
    if outcome == 0: # Failure because constSo is not a constant expression.
708 209 storres
        return None
709 209 storres
    else:
710 209 storres
        return outcome
711 209 storres
# End  pobyso_get_constant_so_sa
712 209 storres
713 57 storres
def pobyso_get_constant_as_rn(ctExpSo):
714 83 storres
    """
715 83 storres
    Legacy function. See pobyso_get_constant_as_rn_so_sa.
716 83 storres
    """
717 57 storres
    return(pobyso_get_constant_as_rn_so_sa(ctExpSo))
718 38 storres
719 56 storres
def pobyso_get_constant_as_rn_so_sa(constExpSo):
720 83 storres
    """
721 83 storres
    Get a Sollya constant as a Sage "real number".
722 83 storres
    The precision of the floating-point number returned is that of the Sollya
723 83 storres
    constant.
724 83 storres
    """
725 218 storres
    #print "Before computing precision of variable..."
726 218 storres
    #pobyso_autoprint(constExpSo)
727 209 storres
    precisionSa  = pobyso_get_prec_of_constant_so_sa(constExpSo)
728 218 storres
    #print "precisionSa:", precisionSa
729 209 storres
    ## If the expression can not be exactly converted, None is returned.
730 209 storres
    #  In this case opt for the Sollya current expression.
731 209 storres
    if precisionSa is None:
732 209 storres
        precisionSa = pobyso_get_prec_so_sa()
733 56 storres
    RRRR = RealField(precisionSa)
734 56 storres
    rnSa = RRRR(0)
735 209 storres
    outcome = sollya_lib_get_constant(get_rn_value(rnSa), constExpSo)
736 209 storres
    if outcome == 0:
737 209 storres
        return None
738 209 storres
    else:
739 209 storres
        return rnSa
740 83 storres
# End pobyso_get_constant_as_rn_so_sa
741 38 storres
742 38 storres
def pobyso_get_constant_as_rn_with_rf(ctExp, realField):
743 83 storres
    """
744 83 storres
    Legacy function. See pobyso_get_constant_as_rn_with_rf_so_sa.
745 83 storres
    """
746 209 storres
    return pobyso_get_constant_as_rn_with_rf_so_sa(ctExp, realField)
747 209 storres
# End pobyso_get_constant_as_rn_with_rf
748 5 storres
749 56 storres
def pobyso_get_constant_as_rn_with_rf_so_sa(ctExpSo, realFieldSa = None):
750 83 storres
    """
751 83 storres
    Get a Sollya constant as a Sage "real number".
752 83 storres
    If no real field is specified, the precision of the floating-point number
753 85 storres
    returned is that of the Sollya constant.
754 83 storres
    Otherwise is is that of the real field. Hence rounding may happen.
755 83 storres
    """
756 56 storres
    if realFieldSa is None:
757 209 storres
        return pobyso_get_constant_as_rn_so_sa(ctExpSo)
758 56 storres
    rnSa = realFieldSa(0)
759 209 storres
    outcome = sollya_lib_get_constant(get_rn_value(rnSa), ctExpSo)
760 209 storres
    if outcome == 0:
761 209 storres
        return None
762 209 storres
    else:
763 209 storres
        return rnSa
764 83 storres
# End pobyso_get_constant_as_rn_with_rf_so_sa
765 38 storres
766 5 storres
def pobyso_get_free_variable_name():
767 83 storres
    """
768 83 storres
    Legacy function. See pobyso_get_free_variable_name_so_sa.
769 83 storres
    """
770 38 storres
    return(pobyso_get_free_variable_name_so_sa())
771 38 storres
772 38 storres
def pobyso_get_free_variable_name_so_sa():
773 209 storres
    return sollya_lib_get_free_variable_name()
774 5 storres
775 38 storres
def pobyso_get_function_arity(expressionSo):
776 83 storres
    """
777 83 storres
    Legacy function. See pobyso_get_function_arity_so_sa.
778 83 storres
    """
779 38 storres
    return(pobyso_get_function_arity_so_sa(expressionSo))
780 38 storres
781 38 storres
def pobyso_get_function_arity_so_sa(expressionSo):
782 5 storres
    arity = c_int(0)
783 38 storres
    sollya_lib_get_function_arity(byref(arity),expressionSo)
784 209 storres
    return int(arity.value)
785 5 storres
786 38 storres
def pobyso_get_head_function(expressionSo):
787 83 storres
    """
788 83 storres
    Legacy function. See pobyso_get_head_function_so_sa.
789 83 storres
    """
790 38 storres
    return(pobyso_get_head_function_so_sa(expressionSo))
791 38 storres
792 38 storres
def pobyso_get_head_function_so_sa(expressionSo):
793 5 storres
    functionType = c_int(0)
794 218 storres
    sollya_lib_get_head_function(byref(functionType), expressionSo)
795 209 storres
    return int(functionType.value)
796 5 storres
797 56 storres
def pobyso_get_interval_from_range_so_sa(soRange, realIntervalFieldSa = None ):
798 53 storres
    """
799 53 storres
    Return the Sage interval corresponding to the Sollya range argument.
800 83 storres
    If no reaIntervalField is passed as an argument, the interval bounds are not
801 56 storres
    rounded: they are elements of RealIntervalField of the "right" precision
802 56 storres
    to hold all the digits.
803 53 storres
    """
804 53 storres
    prec = c_int(0)
805 56 storres
    if realIntervalFieldSa is None:
806 56 storres
        retval = sollya_lib_get_prec_of_range(byref(prec), soRange, None)
807 56 storres
        if retval == 0:
808 209 storres
            return None
809 56 storres
        realIntervalFieldSa = RealIntervalField(prec.value)
810 56 storres
    intervalSa = realIntervalFieldSa(0,0)
811 53 storres
    retval = \
812 53 storres
        sollya_lib_get_interval_from_range(get_interval_value(intervalSa),\
813 53 storres
                                           soRange)
814 53 storres
    if retval == 0:
815 209 storres
        return None
816 209 storres
    return intervalSa
817 56 storres
# End pobyso_get_interval_from_range_so_sa
818 56 storres
819 5 storres
def pobyso_get_list_elements(soObj):
820 38 storres
    """ Legacy function. See pobyso_get_list_elements_so_so. """
821 209 storres
    return pobyso_get_list_elements_so_so(soObj)
822 38 storres
823 117 storres
def pobyso_get_list_elements_so_so(objectListSo):
824 51 storres
    """
825 118 storres
    Get the Sollya list elements as a Sage/Python array of Sollya objects.
826 118 storres

827 118 storres
    INPUT:
828 118 storres
    - objectListSo: a Sollya list of Sollya objects.
829 118 storres

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

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

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