Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 219

Historique | Voir | Annoter | Télécharger (68,56 ko)

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

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

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

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

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

681 118 storres
    INPUT:
682 118 storres
    - objectListSo: a Sollya list of Sollya objects.
683 118 storres

684 118 storres
    OUTPUT:
685 118 storres
    - a Sage/Python tuple made of:
686 118 storres
      - a Sage/Python list of Sollya objects,
687 118 storres
      - a Sage/Python int holding the number of elements,
688 118 storres
      - a Sage/Python int stating (!= 0) that the list is end-elliptic.
689 118 storres
    NOTE::
690 118 storres
        We recover the addresses of the Sollya object from the list of pointers
691 118 storres
        returned by sollya_lib_get_list_elements. The list itself is freed.
692 118 storres
    TODO::
693 118 storres
        Figure out what to do with numElements since the number of elements
694 118 storres
        can easily be recovered from the list itself.
695 118 storres
        Ditto for isEndElliptic.
696 51 storres
    """
697 5 storres
    listAddress = POINTER(c_longlong)()
698 5 storres
    numElements = c_int(0)
699 5 storres
    isEndElliptic = c_int(0)
700 117 storres
    listAsSageList = []
701 5 storres
    result = sollya_lib_get_list_elements(byref(listAddress),\
702 54 storres
                                          byref(numElements),\
703 54 storres
                                          byref(isEndElliptic),\
704 117 storres
                                          objectListSo)
705 5 storres
    if result == 0 :
706 5 storres
        return None
707 5 storres
    for i in xrange(0, numElements.value, 1):
708 118 storres
       #listAsSageList.append(sollya_lib_copy_obj(listAddress[i]))
709 118 storres
       listAsSageList.append(listAddress[i])
710 117 storres
       # Clear each of the elements returned by Sollya.
711 118 storres
       #sollya_lib_clear_obj(listAddress[i])
712 117 storres
    # Free the list itself.
713 117 storres
    sollya_lib_free(listAddress)
714 209 storres
    return (listAsSageList, numElements.value, isEndElliptic.value)
715 5 storres
716 38 storres
def pobyso_get_max_prec_of_exp(soExp):
717 38 storres
    """ Legacy function. See pobyso_get_max_prec_of_exp_so_sa. """
718 209 storres
    return pobyso_get_max_prec_of_exp_so_sa(soExp)
719 5 storres
720 85 storres
def pobyso_get_max_prec_of_exp_so_sa(expSo):
721 38 storres
    """
722 38 storres
    Get the maximum precision used for the numbers in a Sollya expression.
723 52 storres

724 52 storres
    Arguments:
725 52 storres
    soExp -- a Sollya expression pointer
726 52 storres
    Return value:
727 52 storres
    A Python integer
728 38 storres
    TODO:
729 38 storres
    - error management;
730 38 storres
    - correctly deal with numerical type such as DOUBLEEXTENDED.
731 38 storres
    """
732 5 storres
    maxPrecision = 0
733 52 storres
    minConstPrec = 0
734 52 storres
    currentConstPrec = 0
735 85 storres
    operator = pobyso_get_head_function_so_sa(expSo)
736 5 storres
    if (operator != SOLLYA_BASE_FUNC_CONSTANT) and \
737 5 storres
    (operator != SOLLYA_BASE_FUNC_FREE_VARIABLE):
738 85 storres
        (arity, subexpressions) = pobyso_get_subfunctions_so_sa(expSo)
739 5 storres
        for i in xrange(arity):
740 5 storres
            maxPrecisionCandidate = \
741 38 storres
                pobyso_get_max_prec_of_exp_so_sa(subexpressions[i])
742 5 storres
            if maxPrecisionCandidate > maxPrecision:
743 5 storres
                maxPrecision = maxPrecisionCandidate
744 209 storres
        return maxPrecision
745 5 storres
    elif operator == SOLLYA_BASE_FUNC_CONSTANT:
746 85 storres
        #minConstPrec = pobyso_get_min_prec_of_constant_so_sa(expSo)
747 52 storres
        #currentConstPrec = pobyso_get_min_prec_of_constant_so_sa(soExp)
748 52 storres
        #print minConstPrec, " - ", currentConstPrec
749 209 storres
        return pobyso_get_min_prec_of_constant_so_sa(expSo)
750 52 storres
751 5 storres
    elif operator == SOLLYA_BASE_FUNC_FREE_VARIABLE:
752 209 storres
        return 0
753 5 storres
    else:
754 38 storres
        print "pobyso_get_max_prec_of_exp_so_sa: unexepected operator."
755 209 storres
        return 0
756 5 storres
757 85 storres
def pobyso_get_min_prec_of_constant_so_sa(constExpSo):
758 52 storres
    """
759 52 storres
    Get the minimum precision necessary to represent the value of a Sollya
760 52 storres
    constant.
761 52 storres
    MPFR_MIN_PREC and powers of 2 are taken into account.
762 209 storres
    We assume that constExpSo is a pointer to a Sollay constant expression.
763 52 storres
    """
764 85 storres
    constExpAsRnSa = pobyso_get_constant_as_rn_so_sa(constExpSo)
765 85 storres
    return(min_mpfr_size(get_rn_value(constExpAsRnSa)))
766 52 storres
767 200 storres
def pobyso_get_poly_so_sa(polySo, realFieldSa=None):
768 200 storres
    """
769 200 storres
    Convert a Sollya polynomial into a Sage polynomial.
770 209 storres
    Legacy function. Use pobyso_float_poly_so_sa() instead.
771 200 storres
    """
772 213 storres
    return pobyso_float_poly_so_sa(polySo,realFieldSa)
773 200 storres
# End pobyso_get_poly_so_sa
774 200 storres
775 200 storres
def pobyso_get_prec():
776 200 storres
    """ Legacy function. See pobyso_get_prec_so_sa(). """
777 209 storres
    return pobyso_get_prec_so_sa()
778 200 storres
779 200 storres
def pobyso_get_prec_so():
780 200 storres
    """
781 200 storres
    Get the current default precision in Sollya.
782 200 storres
    The return value is a Sollya object.
783 200 storres
    Usefull when modifying the precision back and forth by avoiding
784 200 storres
    extra conversions.
785 200 storres
    """
786 209 storres
    return sollya_lib_get_prec(None)
787 200 storres
788 200 storres
def pobyso_get_prec_so_sa():
789 200 storres
    """
790 200 storres
    Get the current default precision in Sollya.
791 200 storres
    The return value is Sage/Python int.
792 200 storres
    """
793 200 storres
    precSo = sollya_lib_get_prec(None)
794 200 storres
    precSa = c_int(0)
795 200 storres
    sollya_lib_get_constant_as_int(byref(precSa), precSo)
796 200 storres
    sollya_lib_clear_obj(precSo)
797 200 storres
    return int(precSa.value)
798 200 storres
# End pobyso_get_prec_so_sa.
799 200 storres
800 209 storres
def pobyso_get_prec_so_so_sa():
801 209 storres
    """
802 209 storres
    Return the current precision both as a Sollya object and a
803 209 storres
    Sage integer as hybrid tuple.
804 209 storres
    To avoid multiple calls for precision manipulations.
805 209 storres
    """
806 209 storres
    precSo = sollya_lib_get_prec(None)
807 209 storres
    precSa = c_int(0)
808 209 storres
    sollya_lib_get_constant_as_int(byref(precSa), precSo)
809 209 storres
    return (precSo, precSa)
810 200 storres
811 200 storres
def pobyso_get_prec_of_constant(ctExpSo):
812 200 storres
    """ Legacy function. See pobyso_get_prec_of_constant_so_sa. """
813 209 storres
    return pobyso_get_prec_of_constant_so_sa(ctExpSo)
814 200 storres
815 200 storres
def pobyso_get_prec_of_constant_so_sa(ctExpSo):
816 200 storres
    """
817 200 storres
    Tries to find a precision to represent ctExpSo without rounding.
818 200 storres
    If not possible, returns None.
819 200 storres
    """
820 218 storres
    #print "Entering pobyso_get_prec_of_constant_so_sa..."
821 200 storres
    prec = c_int(0)
822 200 storres
    retc = sollya_lib_get_prec_of_constant(byref(prec), ctExpSo, None)
823 200 storres
    if retc == 0:
824 218 storres
        #print "pobyso_get_prec_of_constant_so_sa failed."
825 209 storres
        return None
826 218 storres
    #print "...exiting pobyso_get_prec_of_constant_so_sa."
827 209 storres
    return int(prec.value)
828 200 storres
829 200 storres
def pobyso_get_prec_of_range_so_sa(rangeSo):
830 200 storres
    """
831 200 storres
    Returns the number of bits elements of a range are coded with.
832 200 storres
    """
833 200 storres
    prec = c_int(0)
834 200 storres
    retc = sollya_lib_get_prec_of_range(byref(prec), rangeSo, None)
835 200 storres
    if retc == 0:
836 200 storres
        return(None)
837 209 storres
    return int(prec.value)
838 200 storres
# End pobyso_get_prec_of_range_so_sa()
839 200 storres
840 85 storres
def pobyso_get_sage_exp_from_sollya_exp(sollyaExpSo, realField = RR):
841 38 storres
    """ Legacy function. See pobyso_get_sage_exp_from_sollya_exp_so_sa. """
842 209 storres
    return pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo,
843 209 storres
                                                     realField = RR)
844 38 storres
845 85 storres
def pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, realFieldSa = RR):
846 5 storres
    """
847 38 storres
    Get a Sage expression from a Sollya expression.
848 38 storres
    Currently only tested with polynomials with floating-point coefficients.
849 5 storres
    Notice that, in the returned polynomial, the exponents are RealNumbers.
850 5 storres
    """
851 5 storres
    #pobyso_autoprint(sollyaExp)
852 85 storres
    operatorSa = pobyso_get_head_function_so_sa(sollyaExpSo)
853 83 storres
    sollyaLibFreeVariableName = sollya_lib_get_free_variable_name()
854 213 storres
    ## Get rid of the "_"'s in "_x_", if any.
855 213 storres
    sollyaLibFreeVariableName = re.sub('_', '', sollyaLibFreeVariableName)
856 5 storres
    # Constants and the free variable are special cases.
857 5 storres
    # All other operator are dealt with in the same way.
858 85 storres
    if (operatorSa != SOLLYA_BASE_FUNC_CONSTANT) and \
859 85 storres
       (operatorSa != SOLLYA_BASE_FUNC_FREE_VARIABLE):
860 85 storres
        (aritySa, subexpressionsSa) = pobyso_get_subfunctions_so_sa(sollyaExpSo)
861 85 storres
        if aritySa == 1:
862 85 storres
            sageExpSa = eval(pobyso_function_type_as_string_so_sa(operatorSa) + \
863 85 storres
            "(" + pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], \
864 85 storres
            realFieldSa) + ")")
865 85 storres
        elif aritySa == 2:
866 63 storres
            # We do not get through the preprocessor.
867 63 storres
            # The "^" operator is then a special case.
868 85 storres
            if operatorSa == SOLLYA_BASE_FUNC_POW:
869 85 storres
                operatorAsStringSa = "**"
870 5 storres
            else:
871 85 storres
                operatorAsStringSa = \
872 85 storres
                    pobyso_function_type_as_string_so_sa(operatorSa)
873 85 storres
            sageExpSa = \
874 85 storres
              eval("pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], realFieldSa)"\
875 85 storres
              + " " + operatorAsStringSa + " " + \
876 85 storres
                   "pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[1], realFieldSa)")
877 63 storres
        # We do not know yet how to deal with arity >= 3
878 63 storres
        # (is there any in Sollya anyway?).
879 5 storres
        else:
880 85 storres
            sageExpSa = eval('None')
881 209 storres
        return sageExpSa
882 85 storres
    elif operatorSa == SOLLYA_BASE_FUNC_CONSTANT:
883 5 storres
        #print "This is a constant"
884 85 storres
        return pobyso_get_constant_as_rn_with_rf_so_sa(sollyaExpSo, realFieldSa)
885 85 storres
    elif operatorSa == SOLLYA_BASE_FUNC_FREE_VARIABLE:
886 218 storres
        #print "This is the free variable"
887 209 storres
        return eval(sollyaLibFreeVariableName)
888 5 storres
    else:
889 5 storres
        print "Unexpected"
890 5 storres
        return eval('None')
891 185 storres
# End pobyso_get_sage_exp_from_sollya_exp_so_sa
892 73 storres
893 185 storres
894 38 storres
def pobyso_get_subfunctions(expressionSo):
895 38 storres
    """ Legacy function. See pobyso_get_subfunctions_so_sa. """
896 209 storres
    return pobyso_get_subfunctions_so_sa(expressionSo)
897 200 storres
# End pobyso_get_subfunctions.
898 200 storres
899 38 storres
def pobyso_get_subfunctions_so_sa(expressionSo):
900 38 storres
    """
901 38 storres
    Get the subfunctions of an expression.
902 38 storres
    Return the number of subfunctions and the list of subfunctions addresses.
903 55 storres
    S.T.: Could not figure out another way than that ugly list of declarations
904 83 storres
    to recover the addresses of the subfunctions.
905 83 storres
    We limit ourselves to arity 8 functions.
906 38 storres
    """
907 5 storres
    subf0 = c_int(0)
908 5 storres
    subf1 = c_int(0)
909 5 storres
    subf2 = c_int(0)
910 5 storres
    subf3 = c_int(0)
911 5 storres
    subf4 = c_int(0)
912 5 storres
    subf5 = c_int(0)
913 5 storres
    subf6 = c_int(0)
914 5 storres
    subf7 = c_int(0)
915 5 storres
    subf8 = c_int(0)
916 5 storres
    arity = c_int(0)
917 5 storres
    nullPtr = POINTER(c_int)()
918 38 storres
    sollya_lib_get_subfunctions(expressionSo, byref(arity), \
919 83 storres
      byref(subf0), byref(subf1), byref(subf2), byref(subf3), \
920 83 storres
      byref(subf4), byref(subf5),\
921 83 storres
      byref(subf6), byref(subf7), byref(subf8), nullPtr, None)
922 83 storres
#    byref(cast(subfunctions[0], POINTER(c_int))), \
923 83 storres
#    byref(cast(subfunctions[0], POINTER(c_int))), \
924 83 storres
#    byref(cast(subfunctions[2], POINTER(c_int))), \
925 83 storres
#    byref(cast(subfunctions[3], POINTER(c_int))), \
926 83 storres
#    byref(cast(subfunctions[4], POINTER(c_int))), \
927 83 storres
#    byref(cast(subfunctions[5], POINTER(c_int))), \
928 83 storres
#    byref(cast(subfunctions[6], POINTER(c_int))), \
929 83 storres
#    byref(cast(subfunctions[7], POINTER(c_int))), \
930 5 storres
#    byref(cast(subfunctions[8], POINTER(c_int))), nullPtr)
931 83 storres
    subfunctions = [subf0, subf1, subf2, subf3, subf4, subf5, subf6, subf7, \
932 83 storres
                    subf8]
933 5 storres
    subs = []
934 5 storres
    if arity.value > pobyso_max_arity:
935 38 storres
        return(0,[])
936 5 storres
    for i in xrange(arity.value):
937 5 storres
        subs.append(int(subfunctions[i].value))
938 5 storres
        #print subs[i]
939 209 storres
    return (int(arity.value), subs)
940 200 storres
# End pobyso_get_subfunctions_so_sa
941 5 storres
942 155 storres
def pobyso_guess_degree_sa_sa(functionSa, intervalSa, approxErrorSa,
943 155 storres
                              weightSa=None, degreeBoundSa=None):
944 155 storres
    """
945 155 storres
    Sa_sa variant of the solly_guessdegree function.
946 155 storres
    Return 0 if something goes wrong.
947 155 storres
    """
948 159 storres
    functionAsStringSa = functionSa._assume_str().replace('_SAGE_VAR_', '')
949 154 storres
    functionSo = pobyso_parse_string_sa_so(functionAsStringSa)
950 155 storres
    if pobyso_is_error_so_sa(functionSo):
951 155 storres
        sollya_lib_clear_obj(functionSo)
952 155 storres
        return 0
953 154 storres
    rangeSo = pobyso_interval_to_range_sa_so(intervalSa)
954 155 storres
    # The approximation error is expected to be a floating point number.
955 155 storres
    if pobyso_is_floating_point_number_sa_sa(approxErrorSa):
956 155 storres
        approxErrorSo = pobyso_constant_sa_so(approxErrorSa)
957 155 storres
    else:
958 155 storres
        approxErrorSo = pobyso_constant_sa_so(RR(approxErrorSa))
959 154 storres
    if not weightSa is None:
960 159 storres
        weightAsStringSa = weightSa._assume_str().replace('_SAGE_VAR_', '')
961 154 storres
        weightSo = pobyso_parse_string_sa_so(weightAsStringSa)
962 166 storres
        if pobyso_is_error_so_sa(weightSo):
963 155 storres
            sollya_lib_clear_obj(functionSo)
964 155 storres
            sollya_lib_clear_obj(rangeSo)
965 155 storres
            sollya_lib_clear_obj(approxErrorSo)
966 155 storres
            sollya_lib_clear_obj(weightSo)
967 155 storres
            return 0
968 154 storres
    else:
969 154 storres
        weightSo = None
970 154 storres
    if not degreeBoundSa is None:
971 154 storres
        degreeBoundSo = pobyso_constant_from_int_sa_so(degreeBoundSa)
972 154 storres
    else:
973 154 storres
        degreeBoundSo = None
974 154 storres
    guessedDegreeSa = pobyso_guess_degree_so_sa(functionSo,
975 162 storres
                                                rangeSo,
976 162 storres
                                                approxErrorSo,
977 162 storres
                                                weightSo,
978 162 storres
                                                degreeBoundSo)
979 154 storres
    sollya_lib_clear_obj(functionSo)
980 154 storres
    sollya_lib_clear_obj(rangeSo)
981 155 storres
    sollya_lib_clear_obj(approxErrorSo)
982 154 storres
    if not weightSo is None:
983 154 storres
        sollya_lib_clear_obj(weightSo)
984 154 storres
    if not degreeBoundSo is None:
985 154 storres
        sollya_lib_clear_obj(degreeBoundSo)
986 154 storres
    return guessedDegreeSa
987 154 storres
# End poyso_guess_degree_sa_sa
988 154 storres
989 153 storres
def pobyso_guess_degree_so_sa(functionSo, rangeSo, errorSo, weightSo=None, \
990 154 storres
                              degreeBoundSo=None):
991 154 storres
    """
992 154 storres
    Thin wrapper around the guessdegree function.
993 154 storres
    Nevertheless, some precision control stuff has been appended.
994 154 storres
    """
995 154 storres
    # Deal with Sollya internal precision issues: if it is too small,
996 154 storres
    # compared with the error, increases it to about twice -log2(error).
997 154 storres
    errorSa = pobyso_get_constant_as_rn_with_rf_so_sa(errorSo)
998 154 storres
    log2ErrorSa = errorSa.log2()
999 154 storres
    if log2ErrorSa < 0:
1000 154 storres
        neededPrecisionSa = int(2 * int(-log2ErrorSa) / 64) * 64
1001 154 storres
    else:
1002 154 storres
        neededPrecisionSa = int(2 * int(log2ErrorSa) / 64) * 64
1003 154 storres
    #print "Needed precision:", neededPrecisionSa
1004 154 storres
    currentPrecSa = pobyso_get_prec_so_sa()
1005 154 storres
    if neededPrecisionSa > currentPrecSa:
1006 154 storres
        currentPrecSo = pobyso_get_prec_so()
1007 154 storres
        pobyso_set_prec_sa_so(neededPrecisionSa)
1008 166 storres
    #print "Guessing degree..."
1009 153 storres
    # weightSo and degreeBoundsSo are optional arguments.
1010 162 storres
    # As declared, sollya_lib_guessdegree must take 5 arguments.
1011 153 storres
    if weightSo is None:
1012 162 storres
        degreeRangeSo = sollya_lib_guessdegree(functionSo, rangeSo, errorSo,
1013 162 storres
                                               0, 0, None)
1014 154 storres
    elif degreeBoundSo is None:
1015 153 storres
        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, \
1016 162 storres
                                                errorSo, weightSo, 0, None)
1017 153 storres
    else:
1018 153 storres
        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, errorSo, \
1019 154 storres
                                                weightSo, degreeBoundSo, None)
1020 166 storres
    #print "...degree guess done."
1021 154 storres
    # Restore internal precision, if applicable.
1022 154 storres
    if neededPrecisionSa > currentPrecSa:
1023 154 storres
        pobyso_set_prec_so_so(currentPrecSo)
1024 154 storres
        sollya_lib_clear_obj(currentPrecSo)
1025 154 storres
    degreeIntervalSa = pobyso_range_to_interval_so_sa(degreeRangeSo)
1026 154 storres
    sollya_lib_clear_obj(degreeRangeSo)
1027 154 storres
    # When ok, both bounds match.
1028 154 storres
    # When the degree bound is too low, the upper bound is the degree
1029 154 storres
    # for which the error can be honored.
1030 154 storres
    # When it really goes wrong, the upper bound is infinity.
1031 154 storres
    if degreeIntervalSa.lower() == degreeIntervalSa.upper():
1032 154 storres
        return int(degreeIntervalSa.lower())
1033 154 storres
    else:
1034 154 storres
        if degreeIntervalSa.upper().is_infinity():
1035 154 storres
            return None
1036 154 storres
        else:
1037 154 storres
            return int(degreeIntervalSa.upper())
1038 154 storres
    # End pobyso_guess_degree_so_sa
1039 153 storres
1040 215 storres
def pobyso_inf_so_so(intervalSo):
1041 215 storres
    """
1042 215 storres
    Very thin wrapper around sollya_lib_inf().
1043 215 storres
    """
1044 215 storres
    return sollya_lib_inf(intervalSo)
1045 215 storres
# End pobyso_inf_so_so.
1046 215 storres
1047 53 storres
def pobyso_infnorm_so_so(func, interval, file = None, intervalList = None):
1048 54 storres
    print "Do not use this function. User pobyso_supnorm_so_so instead."
1049 209 storres
    return None
1050 53 storres
1051 84 storres
def pobyso_interval_to_range_sa_so(intervalSa, precisionSa=None):
1052 84 storres
    if precisionSa is None:
1053 84 storres
        precisionSa = intervalSa.parent().precision()
1054 84 storres
    intervalSo = pobyso_bounds_to_range_sa_so(intervalSa.lower(),\
1055 84 storres
                                              intervalSa.upper(),\
1056 84 storres
                                              precisionSa)
1057 209 storres
    return intervalSo
1058 84 storres
# End pobyso_interval_to_range_sa_so
1059 84 storres
1060 155 storres
def pobyso_is_error_so_sa(objSo):
1061 155 storres
    """
1062 155 storres
    Thin wrapper around the sollya_lib_obj_is_error() function.
1063 155 storres
    """
1064 155 storres
    if sollya_lib_obj_is_error(objSo) != 0:
1065 155 storres
        return True
1066 155 storres
    else:
1067 155 storres
        return False
1068 155 storres
# End pobyso_is_error-so_sa
1069 155 storres
1070 155 storres
def pobyso_is_floating_point_number_sa_sa(numberSa):
1071 155 storres
    """
1072 209 storres
    Check whether a Sage number is floating point.
1073 209 storres
    Exception stuff added because numbers other than
1074 209 storres
    floating-point ones do not have the is_real() attribute.
1075 155 storres
    """
1076 209 storres
    try:
1077 209 storres
        return numberSa.is_real()
1078 209 storres
    except AttributeError:
1079 209 storres
        return False
1080 209 storres
# End pobyso_is_floating_piont_number_sa_sa
1081 155 storres
1082 37 storres
def pobyso_lib_init():
1083 37 storres
    sollya_lib_init(None)
1084 116 storres
1085 116 storres
def pobyso_lib_close():
1086 116 storres
    sollya_lib_close(None)
1087 37 storres
1088 85 storres
def pobyso_name_free_variable(freeVariableNameSa):
1089 38 storres
    """ Legacy function. See pobyso_name_free_variable_sa_so. """
1090 85 storres
    pobyso_name_free_variable_sa_so(freeVariableNameSa)
1091 38 storres
1092 85 storres
def pobyso_name_free_variable_sa_so(freeVariableNameSa):
1093 83 storres
    """
1094 83 storres
    Set the free variable name in Sollya from a Sage string.
1095 83 storres
    """
1096 85 storres
    sollya_lib_name_free_variable(freeVariableNameSa)
1097 37 storres
1098 5 storres
def pobyso_parse_string(string):
1099 38 storres
    """ Legacy function. See pobyso_parse_string_sa_so. """
1100 209 storres
    return pobyso_parse_string_sa_so(string)
1101 38 storres
1102 38 storres
def pobyso_parse_string_sa_so(string):
1103 83 storres
    """
1104 155 storres
    Get the Sollya expression computed from a Sage string or
1105 155 storres
    a Sollya error object if parsing failed.
1106 83 storres
    """
1107 209 storres
    return sollya_lib_parse_string(string)
1108 5 storres
1109 200 storres
def pobyso_precision_so_sa(ctExpSo):
1110 209 storres
    """
1111 209 storres
    Computes the necessary precision to represent a number.
1112 209 storres
    If x is not zero, it can be uniquely written as x = m · 2e
1113 209 storres
    where m is an odd integer and e is an integer.
1114 209 storres
    precision(x) returns the number of bits necessary to write m
1115 209 storres
    in binary (i.e. ceil(log2(m))).
1116 209 storres
    """
1117 209 storres
    #TODO: take care of the special case: 0, @NaN@, @Inf@
1118 200 storres
    precisionSo = sollya_lib_precision(ctExpSo)
1119 200 storres
    precisionSa = pobyso_constant_from_int_so_sa(precisionSo)
1120 200 storres
    sollya_lib_clear_obj(precisionSo)
1121 200 storres
    return precisionSa
1122 200 storres
# End pobyso_precision_so_sa
1123 215 storres
1124 217 storres
def pobyso_polynomial_coefficients_progressive_round_so_so(polySo,
1125 217 storres
                                                           funcSo,
1126 217 storres
                                                           icSo,
1127 217 storres
                                                           intervalSo,
1128 217 storres
                                                           itpSo,
1129 217 storres
                                                           ftpSo,
1130 217 storres
                                                           maxPrecSo,
1131 219 storres
                                                           maxErrSo,
1132 219 storres
                                                           debug=False):
1133 219 storres
    if debug:
1134 219 storres
        print "Input arguments:"
1135 219 storres
        pobyso_autoprint(polySo)
1136 219 storres
        pobyso_autoprint(funcSo)
1137 219 storres
        pobyso_autoprint(icSo)
1138 219 storres
        pobyso_autoprint(intervalSo)
1139 219 storres
        pobyso_autoprint(itpSo)
1140 219 storres
        pobyso_autoprint(ftpSo)
1141 219 storres
        pobyso_autoprint(maxPrecSo)
1142 219 storres
        pobyso_autoprint(maxErrSo)
1143 219 storres
        print "________________"
1144 200 storres
1145 217 storres
    ## Higher order function see:
1146 217 storres
    #  http://effbot.org/pyfaq/how-do-you-make-a-higher-order-function-in-python.htm
1147 217 storres
    def precision_decay_ratio_function(degreeSa):
1148 217 storres
        def outer(x):
1149 217 storres
            def inner(x):
1150 217 storres
                we = 3/8
1151 217 storres
                wq = 2/8
1152 217 storres
                a  = 2.2
1153 217 storres
                b  = 2
1154 217 storres
                return we*(exp(x/a)-1) +  wq*((b*x)**2) + (1-we-wq)*x
1155 217 storres
            return  inner(x)/inner(degreeSa)
1156 217 storres
        return outer
1157 217 storres
1158 217 storres
    #
1159 217 storres
    degreeSa        = pobyso_polynomial_degree_so_sa(polySo)
1160 217 storres
    ratio           = precision_decay_ratio_function(degreeSa)
1161 217 storres
    itpSa           = pobyso_constant_from_int_so_sa(itpSo)
1162 217 storres
    ftpSa           = pobyso_constant_from_int_so_sa(ftpSo)
1163 217 storres
    maxPrecSa       = pobyso_constant_from_int_so_sa(maxPrecSo)
1164 217 storres
    maxErrSa        = pobyso_get_constant_as_rn_so_sa(maxErrSo)
1165 219 storres
    if debug:
1166 219 storres
        print "degreeSa:", degreeSa
1167 219 storres
        print "ratio:", ratio
1168 219 storres
        print "itpsSa:", itpSa
1169 219 storres
        print "ftpSa:", ftpSa
1170 219 storres
        print "maxPrecSa:", maxPrecSa
1171 219 storres
        print "maxErrSa:", maxErrSa
1172 217 storres
    lastResPolySo   = None
1173 218 storres
    lastInfNormSo   = None
1174 219 storres
    #print "About to enter the while loop..."
1175 217 storres
    while True:
1176 218 storres
        resPolySo   = pobyso_constant_0_sa_so()
1177 217 storres
        pDeltaSa    = ftpSa - itpSa
1178 217 storres
        for indexSa in reversed(xrange(0,degreeSa+1)):
1179 218 storres
            #print "Index:", indexSa
1180 217 storres
            indexSo = pobyso_constant_from_int_sa_so(indexSa)
1181 217 storres
            #pobyso_autoprint(indexSo)
1182 217 storres
            #print ratio(indexSa)
1183 217 storres
            ctpSa = floor(ftpSa - (pDeltaSa * ratio(indexSa)))
1184 217 storres
            ctpSo = pobyso_constant_from_int_sa_so(ctpSa)
1185 219 storres
            if debug:
1186 219 storres
                print "Index:", indexSa, " - Target precision:",
1187 219 storres
                pobyso_autoprint(ctpSo)
1188 217 storres
            cmonSo  = \
1189 217 storres
                sollya_lib_build_function_mul(sollya_lib_coeff(polySo, indexSo),
1190 217 storres
                                      sollya_lib_build_function_pow( \
1191 217 storres
                                          sollya_lib_build_function_free_variable(), \
1192 217 storres
                                          indexSo))
1193 217 storres
            #pobyso_autoprint(cmonSo)
1194 217 storres
            cmonrSo = pobyso_round_coefficients_single_so_so(cmonSo, ctpSo)
1195 217 storres
            sollya_lib_clear_obj(cmonSo)
1196 217 storres
            #pobyso_autoprint(cmonrSo)
1197 217 storres
            resPolySo = sollya_lib_build_function_add(resPolySo,
1198 217 storres
                                                      cmonrSo)
1199 218 storres
            #pobyso_autoprint(resPolySo)
1200 217 storres
        # End for index
1201 217 storres
        freeVarSo     = sollya_lib_build_function_free_variable()
1202 217 storres
        changeVarSo   = sollya_lib_sub(freeVarSo, icSo)
1203 217 storres
        resPolyCvSo   = sollya_lib_evaluate(resPolySo, changeVarSo)
1204 218 storres
        errFuncSo = sollya_lib_build_function_sub(sollya_lib_copy_obj(funcSo),
1205 218 storres
                                                  resPolyCvSo)
1206 218 storres
        infNormSo = sollya_lib_dirtyinfnorm(errFuncSo, intervalSo)
1207 217 storres
        cerrSa    = pobyso_get_constant_as_rn_so_sa(infNormSo)
1208 219 storres
        if debug:
1209 219 storres
            print "Infnorm (Sollya):",
1210 219 storres
            pobyso_autoprint(infNormSo)
1211 218 storres
        sollya_lib_clear_obj(errFuncSo)
1212 218 storres
        #print "Infnorm  (Sage):", cerrSa
1213 217 storres
        if (cerrSa > maxErrSa):
1214 219 storres
            if debug:
1215 219 storres
                print "Error is too large."
1216 218 storres
            if lastResPolySo is None:
1217 219 storres
                if debug:
1218 219 storres
                    print "Enlarging prec."
1219 217 storres
                ntpSa = floor(ftpSa + ftpSa/50)
1220 217 storres
                ## Can't enlarge (numerical)
1221 217 storres
                if ntpSa == ftpSa:
1222 217 storres
                    sollya_lib_clear_obj(resPolySo)
1223 217 storres
                    return None
1224 217 storres
                ## Can't enlarge (not enough precision left)
1225 217 storres
                if ntpSa > maxPrecSa:
1226 217 storres
                    sollya_lib_clear_obj(resPolySo)
1227 217 storres
                    return None
1228 217 storres
                ftpSa = ntpSa
1229 217 storres
                continue
1230 217 storres
            ## One enlargement took place.
1231 217 storres
            else:
1232 219 storres
                if debug:
1233 219 storres
                    print "Exit with the last before last polynomial."
1234 217 storres
                sollya_lib_clear_obj(resPolySo)
1235 218 storres
                sollya_lib_clear_obj(infNormSo)
1236 218 storres
                return (lastResPolySo, lastInfNormSo)
1237 218 storres
        # cerrSa <= maxErrSa: scrap more bits, possibly.
1238 217 storres
        else:
1239 219 storres
            if debug:
1240 219 storres
                print "Error is too small"
1241 219 storres
            if cerrSa <= (maxErrSa/2):
1242 219 storres
                if debug:
1243 219 storres
                    print "Shrinking prec."
1244 218 storres
                ntpSa = floor(ftpSa - ftpSa/50)
1245 218 storres
                ## Can't shrink (numerical)
1246 218 storres
                if ntpSa == ftpSa:
1247 218 storres
                    if not lastResPolySo is None:
1248 218 storres
                        sollya_lib_clear_obj(lastResPolySo)
1249 218 storres
                    if not lastInfNormSo is None:
1250 218 storres
                        sollya_lib_clear_obj(lastInfNormSo)
1251 218 storres
                    return (resPolySo, infNormSo)
1252 218 storres
                ## Can't shrink (not enough precision left)
1253 218 storres
                if ntpSa <= itpSa:
1254 218 storres
                    if not lastResPolySo is None:
1255 218 storres
                        sollya_lib_clear_obj(lastResPolySo)
1256 218 storres
                    if not lastInfNormSo is None:
1257 218 storres
                        sollya_lib_clear_obj(lastInfNormSo)
1258 218 storres
                    return (resPolySo, infNormSo)
1259 218 storres
                ftpSa = ntpSa
1260 217 storres
                if not lastResPolySo is None:
1261 217 storres
                    sollya_lib_clear_obj(lastResPolySo)
1262 218 storres
                if not lastInfNormSo is None:
1263 218 storres
                    sollya_lib_clear_obj(lastInfNormSo)
1264 218 storres
                lastResPolySo = resPolySo
1265 218 storres
                lastInfNormSo = infNormSo
1266 218 storres
                continue
1267 218 storres
            else: # Error is not that small, just return
1268 217 storres
                if not lastResPolySo is None:
1269 217 storres
                    sollya_lib_clear_obj(lastResPolySo)
1270 218 storres
                if not lastInfNormSo is None:
1271 218 storres
                    sollya_lib_clear_obj(lastInfNormSo)
1272 218 storres
                return (resPolySo, infNormSo)
1273 217 storres
    # End wile True
1274 215 storres
    return None
1275 217 storres
# End pobyso_polynomial_coefficients_progressive_truncate_so_so.
1276 217 storres
1277 217 storres
def pobyso_polynomial_degree_so_sa(polySo):
1278 217 storres
    """
1279 217 storres
    Return the degree of a Sollya polynomial as a Sage int.
1280 217 storres
    """
1281 217 storres
    degreeSo = sollya_lib_degree(polySo)
1282 217 storres
    return pobyso_constant_from_int_so_sa(degreeSo)
1283 217 storres
# End pobyso_polynomial_degree_so_sa
1284 217 storres
1285 217 storres
def pobyso_polynomial_degree_so_so(polySo):
1286 217 storres
    """
1287 217 storres
    Thin wrapper around lib_sollya_degree().
1288 217 storres
    """
1289 217 storres
    return sollya_lib_degree(polySo)
1290 217 storres
# End pobyso_polynomial_degree_so_so
1291 217 storres
1292 5 storres
def pobyso_range(rnLowerBound, rnUpperBound):
1293 38 storres
    """ Legacy function. See pobyso_range_sa_so. """
1294 209 storres
    return pobyso_range_sa_so(rnLowerBound, rnUpperBound)
1295 38 storres
1296 5 storres
1297 85 storres
def pobyso_range_to_interval_so_sa(rangeSo, realIntervalFieldSa = None):
1298 83 storres
    """
1299 83 storres
    Get a Sage interval from a Sollya range.
1300 83 storres
    If no realIntervalField is given as a parameter, the Sage interval
1301 83 storres
    precision is that of the Sollya range.
1302 85 storres
    Otherwise, the precision is that of the realIntervalField. In this case
1303 85 storres
    rounding may happen.
1304 83 storres
    """
1305 85 storres
    if realIntervalFieldSa is None:
1306 56 storres
        precSa = pobyso_get_prec_of_range_so_sa(rangeSo)
1307 85 storres
        realIntervalFieldSa = RealIntervalField(precSa)
1308 56 storres
    intervalSa = \
1309 85 storres
        pobyso_get_interval_from_range_so_sa(rangeSo, realIntervalFieldSa)
1310 209 storres
    return intervalSa
1311 209 storres
# End pobyso_range_to_interval_so_sa
1312 56 storres
1313 209 storres
def pobyso_rat_poly_sa_so(polySa, precSa = None):
1314 209 storres
    """
1315 209 storres
    Create a Sollya polynomial from a Sage rational polynomial.
1316 209 storres
    """
1317 209 storres
    ## TODO: filter arguments.
1318 209 storres
    ## Precision. If no precision is given, use the current precision
1319 209 storres
    #  of Sollya.
1320 209 storres
    if precSa is None:
1321 209 storres
        precSa =  pobyso_get_prec_so_sa()
1322 209 storres
    #print "Precision:",  precSa
1323 209 storres
    RRR = RealField(precSa)
1324 209 storres
    ## Create a Sage polynomial in the "right" precision.
1325 209 storres
    P_RRR = RRR[polySa.variables()[0]]
1326 209 storres
    polyFloatSa = P_RRR(polySa)
1327 213 storres
    ## Make sure no precision is provided: pobyso_float_poly_sa_so will
1328 213 storres
    #  recover it all by itself and not make an extra conversion.
1329 209 storres
    return pobyso_float_poly_sa_so(polyFloatSa)
1330 209 storres
1331 209 storres
# End pobyso_rat_poly_sa_so
1332 209 storres
1333 52 storres
def pobyso_remez_canonical_sa_sa(func, \
1334 52 storres
                                 degree, \
1335 52 storres
                                 lowerBound, \
1336 52 storres
                                 upperBound, \
1337 52 storres
                                 weight = None, \
1338 52 storres
                                 quality = None):
1339 52 storres
    """
1340 52 storres
    All arguments are Sage/Python.
1341 52 storres
    The functions (func and weight) must be passed as expressions or strings.
1342 52 storres
    Otherwise the function fails.
1343 83 storres
    The return value is a Sage polynomial.
1344 52 storres
    """
1345 83 storres
    var('zorglub')    # Dummy variable name for type check only. Type of
1346 83 storres
    # zorglub is "symbolic expression".
1347 52 storres
    polySo = pobyso_remez_canonical_sa_so(func, \
1348 52 storres
                                 degree, \
1349 52 storres
                                 lowerBound, \
1350 52 storres
                                 upperBound, \
1351 85 storres
                                 weight, \
1352 85 storres
                                 quality)
1353 83 storres
    # String test
1354 52 storres
    if parent(func) == parent("string"):
1355 52 storres
        functionSa = eval(func)
1356 52 storres
    # Expression test.
1357 52 storres
    elif type(func) == type(zorglub):
1358 52 storres
        functionSa = func
1359 83 storres
    else:
1360 83 storres
        return None
1361 83 storres
    #
1362 52 storres
    maxPrecision = 0
1363 52 storres
    if polySo is None:
1364 52 storres
        return(None)
1365 52 storres
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
1366 85 storres
    RRRRSa = RealField(maxPrecision)
1367 85 storres
    polynomialRingSa = RRRRSa[functionSa.variables()[0]]
1368 85 storres
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, RRRRSa)
1369 85 storres
    polySa = polynomial(expSa, polynomialRingSa)
1370 83 storres
    sollya_lib_clear_obj(polySo)
1371 52 storres
    return(polySa)
1372 85 storres
# End pobyso_remez_canonical_sa_sa
1373 52 storres
1374 38 storres
def pobyso_remez_canonical(func, \
1375 5 storres
                           degree, \
1376 5 storres
                           lowerBound, \
1377 5 storres
                           upperBound, \
1378 38 storres
                           weight = "1", \
1379 5 storres
                           quality = None):
1380 38 storres
    """ Legacy function. See pobyso_remez_canonical_sa_so. """
1381 51 storres
    return(pobyso_remez_canonical_sa_so(func, \
1382 51 storres
                                        degree, \
1383 51 storres
                                        lowerBound, \
1384 51 storres
                                        upperBound, \
1385 51 storres
                                        weight, \
1386 51 storres
                                        quality))
1387 200 storres
# End pobyso_remez_canonical.
1388 200 storres
1389 38 storres
def pobyso_remez_canonical_sa_so(func, \
1390 38 storres
                                 degree, \
1391 38 storres
                                 lowerBound, \
1392 38 storres
                                 upperBound, \
1393 52 storres
                                 weight = None, \
1394 38 storres
                                 quality = None):
1395 38 storres
    """
1396 38 storres
    All arguments are Sage/Python.
1397 51 storres
    The functions (func and weight) must be passed as expressions or strings.
1398 51 storres
    Otherwise the function fails.
1399 38 storres
    The return value is a pointer to a Sollya function.
1400 38 storres
    """
1401 83 storres
    var('zorglub')    # Dummy variable name for type check only. Type of
1402 83 storres
    # zorglub is "symbolic expression".
1403 85 storres
    currentVariableNameSa = None
1404 52 storres
    # The func argument can be of different types (string,
1405 52 storres
    # symbolic expression...)
1406 38 storres
    if parent(func) == parent("string"):
1407 85 storres
        localFuncSa = eval(func)
1408 85 storres
        if len(localFuncSa.variables()) > 0:
1409 85 storres
            currentVariableNameSa = localFuncSa.variables()[0]
1410 85 storres
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1411 159 storres
            functionSo = \
1412 159 storres
              sollya_lib_parse_string(localFuncSa._assume_str().replace('_SAGE_VAR_', ''))
1413 51 storres
    # Expression test.
1414 52 storres
    elif type(func) == type(zorglub):
1415 52 storres
        # Until we are able to translate Sage expressions into Sollya
1416 52 storres
        # expressions : parse the string version.
1417 85 storres
        if len(func.variables()) > 0:
1418 85 storres
            currentVariableNameSa = func.variables()[0]
1419 85 storres
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1420 159 storres
            functionSo = \
1421 159 storres
              sollya_lib_parse_string(func._assume_str().replace('_SAGE_VAR_', ''))
1422 38 storres
    else:
1423 38 storres
        return(None)
1424 85 storres
    if weight is None: # No weight given -> 1.
1425 52 storres
        weightSo = pobyso_constant_1_sa_so()
1426 85 storres
    elif parent(weight) == parent("string"): # Weight given as string: parse it.
1427 51 storres
        weightSo = sollya_lib_parse_string(func)
1428 85 storres
    elif type(weight) == type(zorglub): # Weight given as symbolice expression.
1429 159 storres
        functionSo = \
1430 159 storres
          sollya_lib_parse_string_sa_so(weight._assume_str().replace('_SAGE_VAR_', ''))
1431 51 storres
    else:
1432 51 storres
        return(None)
1433 5 storres
    degreeSo = pobyso_constant_from_int(degree)
1434 85 storres
    rangeSo = pobyso_bounds_to_range_sa_so(lowerBound, upperBound)
1435 38 storres
    if not quality is None:
1436 38 storres
        qualitySo= pobyso_constant_sa_so(quality)
1437 52 storres
    else:
1438 52 storres
        qualitySo = None
1439 83 storres
1440 83 storres
    remezPolySo = sollya_lib_remez(functionSo, \
1441 83 storres
                                   degreeSo, \
1442 83 storres
                                   rangeSo, \
1443 83 storres
                                   weightSo, \
1444 83 storres
                                   qualitySo, \
1445 83 storres
                                   None)
1446 83 storres
    sollya_lib_clear_obj(functionSo)
1447 83 storres
    sollya_lib_clear_obj(degreeSo)
1448 83 storres
    sollya_lib_clear_obj(rangeSo)
1449 83 storres
    sollya_lib_clear_obj(weightSo)
1450 83 storres
    if not qualitySo is None:
1451 85 storres
        sollya_lib_clear_obj(qualitySo)
1452 83 storres
    return(remezPolySo)
1453 83 storres
# End pobyso_remez_canonical_sa_so
1454 83 storres
1455 38 storres
def pobyso_remez_canonical_so_so(funcSo, \
1456 38 storres
                                 degreeSo, \
1457 38 storres
                                 rangeSo, \
1458 52 storres
                                 weightSo = pobyso_constant_1_sa_so(),\
1459 38 storres
                                 qualitySo = None):
1460 38 storres
    """
1461 38 storres
    All arguments are pointers to Sollya objects.
1462 38 storres
    The return value is a pointer to a Sollya function.
1463 38 storres
    """
1464 38 storres
    if not sollya_lib_obj_is_function(funcSo):
1465 38 storres
        return(None)
1466 38 storres
    return(sollya_lib_remez(funcSo, degreeSo, rangeSo, weightSo, qualitySo, None))
1467 200 storres
# End pobyso_remez_canonical_so_so.
1468 200 storres
1469 215 storres
def pobyso_round_coefficients_single_so_so(polySo, precSo):
1470 215 storres
    """
1471 215 storres
    Create a rounded coefficients polynomial from polynomial argument to
1472 215 storres
    the number of bits in size argument.
1473 215 storres
    All coefficients are set to the same precision.
1474 215 storres
    """
1475 215 storres
    ## TODO: check arguments.
1476 215 storres
    endEllipListSo = pobyso_build_end_elliptic_list_so_so(precSo)
1477 215 storres
    polySo = sollya_lib_roundcoefficients(polySo, endEllipListSo, None)
1478 217 storres
    sollya_lib_clear_obj(endEllipListSo)
1479 215 storres
    #sollya_lib_clear_obj(endEllipListSo)
1480 215 storres
    return polySo
1481 215 storres
1482 215 storres
# End pobyso_round_coefficients_single_so_so
1483 215 storres
1484 5 storres
def pobyso_set_canonical_off():
1485 5 storres
    sollya_lib_set_canonical(sollya_lib_off())
1486 5 storres
1487 5 storres
def pobyso_set_canonical_on():
1488 5 storres
    sollya_lib_set_canonical(sollya_lib_on())
1489 5 storres
1490 5 storres
def pobyso_set_prec(p):
1491 38 storres
    """ Legacy function. See pobyso_set_prec_sa_so. """
1492 85 storres
    pobyso_set_prec_sa_so(p)
1493 38 storres
1494 38 storres
def pobyso_set_prec_sa_so(p):
1495 5 storres
    a = c_int(p)
1496 5 storres
    precSo = c_void_p(sollya_lib_constant_from_int(a))
1497 85 storres
    sollya_lib_set_prec(precSo, None)
1498 215 storres
# End pobyso_set_prec_sa_so.
1499 5 storres
1500 85 storres
def pobyso_set_prec_so_so(newPrecSo):
1501 85 storres
    sollya_lib_set_prec(newPrecSo, None)
1502 215 storres
# End pobyso_set_prec_so_so.
1503 54 storres
1504 215 storres
def pobyso_inf_so_so(intervalSo):
1505 215 storres
    """
1506 215 storres
    Very thin wrapper around sollya_lib_inf().
1507 215 storres
    """
1508 215 storres
    return sollya_lib_inf(intervalSo)
1509 215 storres
# End pobyso_inf_so_so.
1510 215 storres
1511 85 storres
def pobyso_supnorm_so_so(polySo, funcSo, intervalSo, errorTypeSo = None,\
1512 85 storres
                         accuracySo = None):
1513 58 storres
    """
1514 85 storres
    Computes the supnorm of the approximation error between the given
1515 85 storres
    polynomial and function.
1516 85 storres
    errorTypeSo defaults to "absolute".
1517 85 storres
    accuracySo defaults to 2^(-40).
1518 85 storres
    """
1519 85 storres
    if errorTypeSo is None:
1520 85 storres
        errorTypeSo = sollya_lib_absolute(None)
1521 85 storres
        errorTypeIsNone = True
1522 85 storres
    else:
1523 85 storres
        errorTypeIsNone = False
1524 85 storres
    #
1525 85 storres
    if accuracySo is None:
1526 85 storres
        # Notice the **!
1527 85 storres
        accuracySo = pobyso_constant_sa_so(RR(2**(-40)))
1528 85 storres
        accuracyIsNone = True
1529 85 storres
    else:
1530 85 storres
        accuracyIsNone = False
1531 85 storres
    pobyso_autoprint(accuracySo)
1532 85 storres
    resultSo = \
1533 85 storres
        sollya_lib_supnorm(polySo, funcSo, intervalSo, errorTypeSo, \
1534 85 storres
                              accuracySo)
1535 85 storres
    if errorTypeIsNone:
1536 85 storres
        sollya_lib_clear_obj(errorTypeSo)
1537 85 storres
    if accuracyIsNone:
1538 85 storres
        sollya_lib_clear_obj(accuracySo)
1539 85 storres
    return resultSo
1540 85 storres
# End pobyso_supnorm_so_so
1541 85 storres
1542 162 storres
def pobyso_taylor_expansion_no_change_var_so_so(functionSo,
1543 162 storres
                                                degreeSo,
1544 162 storres
                                                rangeSo,
1545 162 storres
                                                errorTypeSo=None,
1546 162 storres
                                                sollyaPrecSo=None):
1547 85 storres
    """
1548 162 storres
    Compute the Taylor expansion without the variable change
1549 162 storres
    x -> x-intervalCenter.
1550 58 storres
    """
1551 58 storres
    # No global change of the working precision.
1552 58 storres
    if not sollyaPrecSo is None:
1553 58 storres
        initialPrecSo = sollya_lib_get_prec(None)
1554 58 storres
        sollya_lib_set_prec(sollyaPrecSo)
1555 85 storres
    # Error type stuff: default to absolute.
1556 85 storres
    if errorTypeSo is None:
1557 85 storres
        errorTypeIsNone = True
1558 85 storres
        errorTypeSo = sollya_lib_absolute(None)
1559 85 storres
    else:
1560 85 storres
        errorTypeIsNone = False
1561 162 storres
    intervalCenterSo = sollya_lib_mid(rangeSo, None)
1562 162 storres
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo,
1563 162 storres
                                         intervalCenterSo,
1564 58 storres
                                         rangeSo, errorTypeSo, None)
1565 117 storres
    # taylorFormListSaSo is a Python list of Sollya objects references that
1566 117 storres
    # are copies of the elements of taylorFormSo.
1567 117 storres
    # pobyso_get_list_elements_so_so clears taylorFormSo.
1568 162 storres
    (taylorFormListSaSo, numElementsSa, isEndEllipticSa) = \
1569 58 storres
        pobyso_get_list_elements_so_so(taylorFormSo)
1570 162 storres
    polySo = sollya_lib_copy_obj(taylorFormListSaSo[0])
1571 162 storres
    #print "Num elements:", numElementsSa
1572 162 storres
    sollya_lib_clear_obj(taylorFormSo)
1573 162 storres
    #polySo = taylorFormListSaSo[0]
1574 162 storres
    #errorRangeSo = sollya_lib_copy_obj(taylorFormListSaSo[2])
1575 162 storres
    errorRangeSo = taylorFormListSaSo[2]
1576 181 storres
    # No copy_obj needed here: a new objects are created.
1577 181 storres
    maxErrorSo    = sollya_lib_sup(errorRangeSo)
1578 181 storres
    minErrorSo    = sollya_lib_inf(errorRangeSo)
1579 181 storres
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
1580 181 storres
    absMinErrorSo = sollya_lib_abs(minErrorSo)
1581 181 storres
    sollya_lib_clear_obj(maxErrorSo)
1582 181 storres
    sollya_lib_clear_obj(minErrorSo)
1583 181 storres
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
1584 181 storres
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
1585 58 storres
    # If changed, reset the Sollya working precision.
1586 58 storres
    if not sollyaPrecSo is None:
1587 58 storres
        sollya_lib_set_prec(initialPrecSo)
1588 83 storres
        sollya_lib_clear_obj(initialPrecSo)
1589 85 storres
    if errorTypeIsNone:
1590 85 storres
        sollya_lib_clear_obj(errorTypeSo)
1591 162 storres
    pobyso_clear_taylorform_sa_so(taylorFormListSaSo)
1592 181 storres
    if absMaxErrorSa > absMinErrorSa:
1593 181 storres
        sollya_lib_clear_obj(absMinErrorSo)
1594 181 storres
        return((polySo, intervalCenterSo, absMaxErrorSo))
1595 181 storres
    else:
1596 181 storres
        sollya_lib_clear_obj(absMaxErrorSo)
1597 181 storres
        return((polySo, intervalCenterSo, absMinErrorSo))
1598 162 storres
# end pobyso_taylor_expansion_no_change_var_so_so
1599 58 storres
1600 162 storres
def pobyso_taylor_expansion_with_change_var_so_so(functionSo, degreeSo, \
1601 162 storres
                                                  rangeSo, \
1602 162 storres
                                                  errorTypeSo=None, \
1603 162 storres
                                                  sollyaPrecSo=None):
1604 58 storres
    """
1605 162 storres
    Compute the Taylor expansion with the variable change
1606 162 storres
    x -> (x-intervalCenter) included.
1607 58 storres
    """
1608 56 storres
    # No global change of the working precision.
1609 56 storres
    if not sollyaPrecSo is None:
1610 56 storres
        initialPrecSo = sollya_lib_get_prec(None)
1611 56 storres
        sollya_lib_set_prec(sollyaPrecSo)
1612 162 storres
    #
1613 85 storres
    # Error type stuff: default to absolute.
1614 85 storres
    if errorTypeSo is None:
1615 85 storres
        errorTypeIsNone = True
1616 85 storres
        errorTypeSo = sollya_lib_absolute(None)
1617 85 storres
    else:
1618 85 storres
        errorTypeIsNone = False
1619 162 storres
    intervalCenterSo = sollya_lib_mid(rangeSo)
1620 162 storres
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo, \
1621 162 storres
                                         intervalCenterSo, \
1622 56 storres
                                         rangeSo, errorTypeSo, None)
1623 116 storres
    # taylorFormListSaSo is a Python list of Sollya objects references that
1624 116 storres
    # are copies of the elements of taylorFormSo.
1625 116 storres
    # pobyso_get_list_elements_so_so clears taylorFormSo.
1626 162 storres
    (taylorFormListSo, numElements, isEndElliptic) = \
1627 56 storres
        pobyso_get_list_elements_so_so(taylorFormSo)
1628 162 storres
    polySo = taylorFormListSo[0]
1629 162 storres
    errorRangeSo = taylorFormListSo[2]
1630 181 storres
    maxErrorSo    = sollya_lib_sup(errorRangeSo)
1631 181 storres
    minErrorSo    = sollya_lib_inf(errorRangeSo)
1632 181 storres
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
1633 181 storres
    absMinErrorSo = sollya_lib_abs(minErrorSo)
1634 181 storres
    sollya_lib_clear_obj(maxErrorSo)
1635 181 storres
    sollya_lib_clear_obj(minErrorSo)
1636 181 storres
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
1637 181 storres
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
1638 162 storres
    changeVarExpSo = sollya_lib_build_function_sub(\
1639 162 storres
                       sollya_lib_build_function_free_variable(),\
1640 162 storres
                       sollya_lib_copy_obj(intervalCenterSo))
1641 181 storres
    polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpSo)
1642 181 storres
    sollya_lib_clear_obj(polySo)
1643 162 storres
    sollya_lib_clear_obj(changeVarExpSo)
1644 56 storres
    # If changed, reset the Sollya working precision.
1645 56 storres
    if not sollyaPrecSo is None:
1646 56 storres
        sollya_lib_set_prec(initialPrecSo)
1647 63 storres
        sollya_lib_clear_obj(initialPrecSo)
1648 85 storres
    if errorTypeIsNone:
1649 85 storres
        sollya_lib_clear_obj(errorTypeSo)
1650 162 storres
    sollya_lib_clear_obj(taylorFormSo)
1651 162 storres
    # Do not clear maxErrorSo.
1652 181 storres
    if absMaxErrorSa > absMinErrorSa:
1653 181 storres
        sollya_lib_clear_obj(absMinErrorSo)
1654 181 storres
        return((polyVarChangedSo, intervalCenterSo, absMaxErrorSo))
1655 181 storres
    else:
1656 181 storres
        sollya_lib_clear_obj(absMaxErrorSo)
1657 181 storres
        return((polyVarChangedSo, intervalCenterSo, absMinErrorSo))
1658 162 storres
# end pobyso_taylor_expansion_with_change_var_so_so
1659 56 storres
1660 5 storres
def pobyso_taylor(function, degree, point):
1661 38 storres
    """ Legacy function. See pobysoTaylor_so_so. """
1662 38 storres
    return(pobyso_taylor_so_so(function, degree, point))
1663 38 storres
1664 56 storres
def pobyso_taylor_so_so(functionSo, degreeSo, pointSo):
1665 56 storres
    return(sollya_lib_taylor(functionSo, degreeSo, pointSo))
1666 5 storres
1667 85 storres
def pobyso_taylorform(function, degree, point = None,
1668 85 storres
                      interval = None, errorType=None):
1669 85 storres
    """ Legacy function. See pobyso_taylorform_sa_sa;"""
1670 38 storres
1671 38 storres
def pobyso_taylorform_sa_sa(functionSa, \
1672 84 storres
                            degreeSa, \
1673 84 storres
                            pointSa, \
1674 84 storres
                            intervalSa=None, \
1675 84 storres
                            errorTypeSa=None, \
1676 84 storres
                            precisionSa=None):
1677 37 storres
    """
1678 85 storres
    Compute the Taylor form of 'degreeSa' for 'functionSa' at 'pointSa'
1679 85 storres
    for 'intervalSa' with 'errorTypeSa' (a string) using 'precisionSa'.
1680 37 storres
    point: must be a Real or a Real interval.
1681 37 storres
    return the Taylor form as an array
1682 83 storres
    TODO: take care of the interval and of the point when it is an interval;
1683 38 storres
          when errorType is not None;
1684 83 storres
          take care of the other elements of the Taylor form (coefficients
1685 83 storres
          errors and delta.
1686 37 storres
    """
1687 37 storres
    # Absolute as the default error.
1688 84 storres
    if errorTypeSa is None:
1689 37 storres
        errorTypeSo = sollya_lib_absolute()
1690 84 storres
    elif errorTypeSa == "relative":
1691 84 storres
        errorTypeSo = sollya_lib_relative()
1692 84 storres
    elif errortypeSa == "absolute":
1693 84 storres
        errorTypeSo = sollya_lib_absolute()
1694 37 storres
    else:
1695 84 storres
        # No clean up needed.
1696 84 storres
        return None
1697 84 storres
    # Global precision stuff
1698 84 storres
    precisionChangedSa = False
1699 84 storres
    currentSollyaPrecSo = pobyso_get_prec_so()
1700 84 storres
    currentSollyaPrecSa = pobyso_constant_from_int_so_sa(currentSollyaPrecSo)
1701 84 storres
    if not precisionSa is None:
1702 84 storres
        if precisionSa > currentSollyaPrecSa:
1703 84 storres
            pobyso_set_prec_sa_so(precisionSa)
1704 84 storres
            precisionChangedSa = True
1705 84 storres
1706 85 storres
    if len(functionSa.variables()) > 0:
1707 85 storres
        varSa = functionSa.variables()[0]
1708 85 storres
        pobyso_name_free_variable_sa_so(str(varSa))
1709 84 storres
    # In any case (point or interval) the parent of pointSa has a precision
1710 84 storres
    # method.
1711 84 storres
    pointPrecSa = pointSa.parent().precision()
1712 84 storres
    if precisionSa > pointPrecSa:
1713 84 storres
        pointPrecSa = precisionSa
1714 84 storres
    # In any case (point or interval) pointSa has a base_ring() method.
1715 84 storres
    pointBaseRingString = str(pointSa.base_ring())
1716 84 storres
    if re.search('Interval', pointBaseRingString) is None: # Point
1717 84 storres
        pointSo = pobyso_constant_sa_so(pointSa, pointPrecSa)
1718 84 storres
    else: # Interval.
1719 84 storres
        pointSo = pobyso_interval_to_range_sa_so(pointSa, pointPrecSa)
1720 37 storres
    # Sollyafy the function.
1721 159 storres
    functionSo = pobyso_parse_string_sa_so(functionSa._assume_str().replace('_SAGE_VAR_', ''))
1722 37 storres
    if sollya_lib_obj_is_error(functionSo):
1723 37 storres
        print "pobyso_tailorform: function string can't be parsed!"
1724 37 storres
        return None
1725 37 storres
    # Sollyafy the degree
1726 84 storres
    degreeSo = sollya_lib_constant_from_int(int(degreeSa))
1727 37 storres
    # Sollyafy the point
1728 37 storres
    # Call Sollya
1729 83 storres
    taylorFormSo = \
1730 83 storres
        sollya_lib_taylorform(functionSo, degreeSo, pointSo, errorTypeSo,\
1731 37 storres
                                         None)
1732 85 storres
    sollya_lib_clear_obj(functionSo)
1733 85 storres
    sollya_lib_clear_obj(degreeSo)
1734 85 storres
    sollya_lib_clear_obj(pointSo)
1735 85 storres
    sollya_lib_clear_obj(errorTypeSo)
1736 38 storres
    (tfsAsList, numElements, isEndElliptic) = \
1737 38 storres
            pobyso_get_list_elements_so_so(taylorFormSo)
1738 37 storres
    polySo = tfsAsList[0]
1739 38 storres
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
1740 37 storres
    polyRealField = RealField(maxPrecision)
1741 38 storres
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, polyRealField)
1742 84 storres
    if precisionChangedSa:
1743 84 storres
        sollya_lib_set_prec(currentSollyaPrecSo)
1744 84 storres
        sollya_lib_clear_obj(currentSollyaPrecSo)
1745 37 storres
    polynomialRing = polyRealField[str(varSa)]
1746 37 storres
    polySa = polynomial(expSa, polynomialRing)
1747 37 storres
    taylorFormSa = [polySa]
1748 85 storres
    # Final clean-up.
1749 85 storres
    sollya_lib_clear_obj(taylorFormSo)
1750 51 storres
    return(taylorFormSa)
1751 51 storres
# End pobyso_taylor_form_sa_sa
1752 54 storres
1753 54 storres
def pobyso_taylorform_so_so(functionSo, degreeSo, pointSo, intervalSo=None, \
1754 54 storres
                            errorTypeSo=None):
1755 54 storres
    createdErrorType = False
1756 51 storres
    if errorTypeSo is None:
1757 51 storres
        errorTypeSo = sollya_lib_absolute()
1758 54 storres
        createdErrorType = True
1759 51 storres
    else:
1760 51 storres
        #TODO: deal with the other case.
1761 51 storres
        pass
1762 51 storres
    if intervalSo is None:
1763 54 storres
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
1764 54 storres
                                         errorTypeSo, None)
1765 51 storres
    else:
1766 54 storres
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
1767 54 storres
                                         intervalSo, errorTypeSo, None)
1768 54 storres
    if createdErrorType:
1769 54 storres
        sollya_lib_clear_obj(errorTypeSo)
1770 215 storres
    return resultSo
1771 51 storres
1772 37 storres
1773 37 storres
def pobyso_univar_polynomial_print_reverse(polySa):
1774 51 storres
    """ Legacy function. See pobyso_univar_polynomial_print_reverse_sa_sa. """
1775 51 storres
    return(pobyso_univar_polynomial_print_reverse_sa_sa(polySa))
1776 38 storres
1777 51 storres
def pobyso_univar_polynomial_print_reverse_sa_sa(polySa):
1778 37 storres
    """
1779 37 storres
    Return the string representation of a univariate polynomial with
1780 38 storres
    monomials ordered in the x^0..x^n order of the monomials.
1781 37 storres
    Remember: Sage
1782 37 storres
    """
1783 37 storres
    polynomialRing = polySa.base_ring()
1784 37 storres
    # A very expensive solution:
1785 37 storres
    # -create a fake multivariate polynomial field with only one variable,
1786 37 storres
    #   specifying a negative lexicographical order;
1787 37 storres
    mpolynomialRing = PolynomialRing(polynomialRing.base(), \
1788 37 storres
                                     polynomialRing.variable_name(), \
1789 37 storres
                                     1, order='neglex')
1790 37 storres
    # - convert the univariate argument polynomial into a multivariate
1791 37 storres
    #   version;
1792 37 storres
    p = mpolynomialRing(polySa)
1793 37 storres
    # - return the string representation of the converted form.
1794 37 storres
    # There is no simple str() method defined for p's class.
1795 37 storres
    return(p.__str__())
1796 5 storres
#
1797 5 storres
print pobyso_get_prec()
1798 5 storres
pobyso_set_prec(165)
1799 5 storres
print pobyso_get_prec()
1800 5 storres
a=100
1801 5 storres
print type(a)
1802 5 storres
id(a)
1803 5 storres
print "Max arity: ", pobyso_max_arity
1804 5 storres
print "Function tripleDouble (43) as a string: ", pobyso_function_type_as_string(43)
1805 56 storres
print "Function None (44) as a string: ", pobyso_function_type_as_string(44)
1806 56 storres
print "...Pobyso check done"