Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 218

Historique | Voir | Annoter | Télécharger (68,2 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 217 storres
                                                           maxErrSo):
1132 215 storres
    print "Input arguments:"
1133 218 storres
    #pobyso_autoprint(polySo)
1134 218 storres
    #pobyso_autoprint(funcSo)
1135 218 storres
    #pobyso_autoprint(icSo)
1136 218 storres
    #pobyso_autoprint(intervalSo)
1137 218 storres
    #pobyso_autoprint(itpSo)
1138 218 storres
    #pobyso_autoprint(ftpSo)
1139 218 storres
    #pobyso_autoprint(maxPrecSo)
1140 218 storres
    #pobyso_autoprint(maxErrSo)
1141 218 storres
    #print "________________"
1142 200 storres
1143 217 storres
    ## Higher order function see:
1144 217 storres
    #  http://effbot.org/pyfaq/how-do-you-make-a-higher-order-function-in-python.htm
1145 217 storres
    def precision_decay_ratio_function(degreeSa):
1146 217 storres
        def outer(x):
1147 217 storres
            def inner(x):
1148 217 storres
                we = 3/8
1149 217 storres
                wq = 2/8
1150 217 storres
                a  = 2.2
1151 217 storres
                b  = 2
1152 217 storres
                return we*(exp(x/a)-1) +  wq*((b*x)**2) + (1-we-wq)*x
1153 217 storres
            return  inner(x)/inner(degreeSa)
1154 217 storres
        return outer
1155 217 storres
1156 217 storres
    #
1157 217 storres
    degreeSa        = pobyso_polynomial_degree_so_sa(polySo)
1158 218 storres
    print "degreeSa:", degreeSa
1159 217 storres
    ratio           = precision_decay_ratio_function(degreeSa)
1160 218 storres
    print "ratio:", ratio
1161 217 storres
    itpSa           = pobyso_constant_from_int_so_sa(itpSo)
1162 218 storres
    print "itpsSa:", itpSa
1163 217 storres
    ftpSa           = pobyso_constant_from_int_so_sa(ftpSo)
1164 218 storres
    print "ftpSa:", ftpSa
1165 217 storres
    maxPrecSa       = pobyso_constant_from_int_so_sa(maxPrecSo)
1166 218 storres
    print "maxPrecSa:", maxPrecSa
1167 217 storres
    maxErrSa        = pobyso_get_constant_as_rn_so_sa(maxErrSo)
1168 218 storres
    print "maxErrSa:", maxErrSa
1169 217 storres
    lastResPolySo   = None
1170 218 storres
    lastInfNormSo   = None
1171 218 storres
    print "About to enter the while loop..."
1172 217 storres
    while True:
1173 218 storres
        resPolySo   = pobyso_constant_0_sa_so()
1174 217 storres
        pDeltaSa    = ftpSa - itpSa
1175 217 storres
        for indexSa in reversed(xrange(0,degreeSa+1)):
1176 218 storres
            #print "Index:", indexSa
1177 217 storres
            indexSo = pobyso_constant_from_int_sa_so(indexSa)
1178 217 storres
            #pobyso_autoprint(indexSo)
1179 217 storres
            #print ratio(indexSa)
1180 217 storres
            ctpSa = floor(ftpSa - (pDeltaSa * ratio(indexSa)))
1181 217 storres
            ctpSo = pobyso_constant_from_int_sa_so(ctpSa)
1182 217 storres
            print "Index:", indexSa, " - Target precision:",
1183 217 storres
            pobyso_autoprint(ctpSo)
1184 217 storres
            cmonSo  = \
1185 217 storres
                sollya_lib_build_function_mul(sollya_lib_coeff(polySo, indexSo),
1186 217 storres
                                      sollya_lib_build_function_pow( \
1187 217 storres
                                          sollya_lib_build_function_free_variable(), \
1188 217 storres
                                          indexSo))
1189 217 storres
            #pobyso_autoprint(cmonSo)
1190 217 storres
            cmonrSo = pobyso_round_coefficients_single_so_so(cmonSo, ctpSo)
1191 217 storres
            sollya_lib_clear_obj(cmonSo)
1192 217 storres
            #pobyso_autoprint(cmonrSo)
1193 217 storres
            resPolySo = sollya_lib_build_function_add(resPolySo,
1194 217 storres
                                                      cmonrSo)
1195 218 storres
            #pobyso_autoprint(resPolySo)
1196 217 storres
        # End for index
1197 217 storres
        freeVarSo     = sollya_lib_build_function_free_variable()
1198 217 storres
        changeVarSo   = sollya_lib_sub(freeVarSo, icSo)
1199 217 storres
        resPolyCvSo   = sollya_lib_evaluate(resPolySo, changeVarSo)
1200 218 storres
        errFuncSo = sollya_lib_build_function_sub(sollya_lib_copy_obj(funcSo),
1201 218 storres
                                                  resPolyCvSo)
1202 218 storres
        infNormSo = sollya_lib_dirtyinfnorm(errFuncSo, intervalSo)
1203 217 storres
        cerrSa    = pobyso_get_constant_as_rn_so_sa(infNormSo)
1204 218 storres
        print "Infnorm (Sollya):", pobyso_autoprint(infNormSo)
1205 218 storres
        sollya_lib_clear_obj(errFuncSo)
1206 218 storres
        #print "Infnorm  (Sage):", cerrSa
1207 217 storres
        if (cerrSa > maxErrSa):
1208 217 storres
            print "Error is too large."
1209 218 storres
            if lastResPolySo is None:
1210 217 storres
                print "Enlarging prec."
1211 217 storres
                ntpSa = floor(ftpSa + ftpSa/50)
1212 217 storres
                ## Can't enlarge (numerical)
1213 217 storres
                if ntpSa == ftpSa:
1214 217 storres
                    sollya_lib_clear_obj(resPolySo)
1215 217 storres
                    return None
1216 217 storres
                ## Can't enlarge (not enough precision left)
1217 217 storres
                if ntpSa > maxPrecSa:
1218 217 storres
                    sollya_lib_clear_obj(resPolySo)
1219 217 storres
                    return None
1220 217 storres
                ftpSa = ntpSa
1221 217 storres
                continue
1222 217 storres
            ## One enlargement took place.
1223 217 storres
            else:
1224 218 storres
                print "Exit with the last before last polynomial."
1225 217 storres
                sollya_lib_clear_obj(resPolySo)
1226 218 storres
                sollya_lib_clear_obj(infNormSo)
1227 218 storres
                return (lastResPolySo, lastInfNormSo)
1228 218 storres
        # cerrSa <= maxErrSa: scrap more bits, possibly.
1229 217 storres
        else:
1230 218 storres
            print "Error is too small"
1231 218 storres
            if cerrSa <= (maxErrSa/2):
1232 218 storres
                print "Shrinking prec."
1233 218 storres
                ntpSa = floor(ftpSa - ftpSa/50)
1234 218 storres
                ## Can't shrink (numerical)
1235 218 storres
                if ntpSa == ftpSa:
1236 218 storres
                    if not lastResPolySo is None:
1237 218 storres
                        sollya_lib_clear_obj(lastResPolySo)
1238 218 storres
                    if not lastInfNormSo is None:
1239 218 storres
                        sollya_lib_clear_obj(lastInfNormSo)
1240 218 storres
                    return (resPolySo, infNormSo)
1241 218 storres
                ## Can't shrink (not enough precision left)
1242 218 storres
                if ntpSa <= itpSa:
1243 218 storres
                    if not lastResPolySo is None:
1244 218 storres
                        sollya_lib_clear_obj(lastResPolySo)
1245 218 storres
                    if not lastInfNormSo is None:
1246 218 storres
                        sollya_lib_clear_obj(lastInfNormSo)
1247 218 storres
                    return (resPolySo, infNormSo)
1248 218 storres
                ftpSa = ntpSa
1249 217 storres
                if not lastResPolySo is None:
1250 217 storres
                    sollya_lib_clear_obj(lastResPolySo)
1251 218 storres
                if not lastInfNormSo is None:
1252 218 storres
                    sollya_lib_clear_obj(lastInfNormSo)
1253 218 storres
                lastResPolySo = resPolySo
1254 218 storres
                lastInfNormSo = infNormSo
1255 218 storres
                continue
1256 218 storres
            else: # Error is not that small, just return
1257 217 storres
                if not lastResPolySo is None:
1258 217 storres
                    sollya_lib_clear_obj(lastResPolySo)
1259 218 storres
                if not lastInfNormSo is None:
1260 218 storres
                    sollya_lib_clear_obj(lastInfNormSo)
1261 218 storres
                return (resPolySo, infNormSo)
1262 217 storres
    # End wile True
1263 215 storres
    return None
1264 217 storres
# End pobyso_polynomial_coefficients_progressive_truncate_so_so.
1265 217 storres
1266 217 storres
def pobyso_polynomial_degree_so_sa(polySo):
1267 217 storres
    """
1268 217 storres
    Return the degree of a Sollya polynomial as a Sage int.
1269 217 storres
    """
1270 217 storres
    degreeSo = sollya_lib_degree(polySo)
1271 217 storres
    return pobyso_constant_from_int_so_sa(degreeSo)
1272 217 storres
# End pobyso_polynomial_degree_so_sa
1273 217 storres
1274 217 storres
def pobyso_polynomial_degree_so_so(polySo):
1275 217 storres
    """
1276 217 storres
    Thin wrapper around lib_sollya_degree().
1277 217 storres
    """
1278 217 storres
    return sollya_lib_degree(polySo)
1279 217 storres
# End pobyso_polynomial_degree_so_so
1280 217 storres
1281 5 storres
def pobyso_range(rnLowerBound, rnUpperBound):
1282 38 storres
    """ Legacy function. See pobyso_range_sa_so. """
1283 209 storres
    return pobyso_range_sa_so(rnLowerBound, rnUpperBound)
1284 38 storres
1285 5 storres
1286 85 storres
def pobyso_range_to_interval_so_sa(rangeSo, realIntervalFieldSa = None):
1287 83 storres
    """
1288 83 storres
    Get a Sage interval from a Sollya range.
1289 83 storres
    If no realIntervalField is given as a parameter, the Sage interval
1290 83 storres
    precision is that of the Sollya range.
1291 85 storres
    Otherwise, the precision is that of the realIntervalField. In this case
1292 85 storres
    rounding may happen.
1293 83 storres
    """
1294 85 storres
    if realIntervalFieldSa is None:
1295 56 storres
        precSa = pobyso_get_prec_of_range_so_sa(rangeSo)
1296 85 storres
        realIntervalFieldSa = RealIntervalField(precSa)
1297 56 storres
    intervalSa = \
1298 85 storres
        pobyso_get_interval_from_range_so_sa(rangeSo, realIntervalFieldSa)
1299 209 storres
    return intervalSa
1300 209 storres
# End pobyso_range_to_interval_so_sa
1301 56 storres
1302 209 storres
def pobyso_rat_poly_sa_so(polySa, precSa = None):
1303 209 storres
    """
1304 209 storres
    Create a Sollya polynomial from a Sage rational polynomial.
1305 209 storres
    """
1306 209 storres
    ## TODO: filter arguments.
1307 209 storres
    ## Precision. If no precision is given, use the current precision
1308 209 storres
    #  of Sollya.
1309 209 storres
    if precSa is None:
1310 209 storres
        precSa =  pobyso_get_prec_so_sa()
1311 209 storres
    #print "Precision:",  precSa
1312 209 storres
    RRR = RealField(precSa)
1313 209 storres
    ## Create a Sage polynomial in the "right" precision.
1314 209 storres
    P_RRR = RRR[polySa.variables()[0]]
1315 209 storres
    polyFloatSa = P_RRR(polySa)
1316 213 storres
    ## Make sure no precision is provided: pobyso_float_poly_sa_so will
1317 213 storres
    #  recover it all by itself and not make an extra conversion.
1318 209 storres
    return pobyso_float_poly_sa_so(polyFloatSa)
1319 209 storres
1320 209 storres
# End pobyso_rat_poly_sa_so
1321 209 storres
1322 52 storres
def pobyso_remez_canonical_sa_sa(func, \
1323 52 storres
                                 degree, \
1324 52 storres
                                 lowerBound, \
1325 52 storres
                                 upperBound, \
1326 52 storres
                                 weight = None, \
1327 52 storres
                                 quality = None):
1328 52 storres
    """
1329 52 storres
    All arguments are Sage/Python.
1330 52 storres
    The functions (func and weight) must be passed as expressions or strings.
1331 52 storres
    Otherwise the function fails.
1332 83 storres
    The return value is a Sage polynomial.
1333 52 storres
    """
1334 83 storres
    var('zorglub')    # Dummy variable name for type check only. Type of
1335 83 storres
    # zorglub is "symbolic expression".
1336 52 storres
    polySo = pobyso_remez_canonical_sa_so(func, \
1337 52 storres
                                 degree, \
1338 52 storres
                                 lowerBound, \
1339 52 storres
                                 upperBound, \
1340 85 storres
                                 weight, \
1341 85 storres
                                 quality)
1342 83 storres
    # String test
1343 52 storres
    if parent(func) == parent("string"):
1344 52 storres
        functionSa = eval(func)
1345 52 storres
    # Expression test.
1346 52 storres
    elif type(func) == type(zorglub):
1347 52 storres
        functionSa = func
1348 83 storres
    else:
1349 83 storres
        return None
1350 83 storres
    #
1351 52 storres
    maxPrecision = 0
1352 52 storres
    if polySo is None:
1353 52 storres
        return(None)
1354 52 storres
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
1355 85 storres
    RRRRSa = RealField(maxPrecision)
1356 85 storres
    polynomialRingSa = RRRRSa[functionSa.variables()[0]]
1357 85 storres
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, RRRRSa)
1358 85 storres
    polySa = polynomial(expSa, polynomialRingSa)
1359 83 storres
    sollya_lib_clear_obj(polySo)
1360 52 storres
    return(polySa)
1361 85 storres
# End pobyso_remez_canonical_sa_sa
1362 52 storres
1363 38 storres
def pobyso_remez_canonical(func, \
1364 5 storres
                           degree, \
1365 5 storres
                           lowerBound, \
1366 5 storres
                           upperBound, \
1367 38 storres
                           weight = "1", \
1368 5 storres
                           quality = None):
1369 38 storres
    """ Legacy function. See pobyso_remez_canonical_sa_so. """
1370 51 storres
    return(pobyso_remez_canonical_sa_so(func, \
1371 51 storres
                                        degree, \
1372 51 storres
                                        lowerBound, \
1373 51 storres
                                        upperBound, \
1374 51 storres
                                        weight, \
1375 51 storres
                                        quality))
1376 200 storres
# End pobyso_remez_canonical.
1377 200 storres
1378 38 storres
def pobyso_remez_canonical_sa_so(func, \
1379 38 storres
                                 degree, \
1380 38 storres
                                 lowerBound, \
1381 38 storres
                                 upperBound, \
1382 52 storres
                                 weight = None, \
1383 38 storres
                                 quality = None):
1384 38 storres
    """
1385 38 storres
    All arguments are Sage/Python.
1386 51 storres
    The functions (func and weight) must be passed as expressions or strings.
1387 51 storres
    Otherwise the function fails.
1388 38 storres
    The return value is a pointer to a Sollya function.
1389 38 storres
    """
1390 83 storres
    var('zorglub')    # Dummy variable name for type check only. Type of
1391 83 storres
    # zorglub is "symbolic expression".
1392 85 storres
    currentVariableNameSa = None
1393 52 storres
    # The func argument can be of different types (string,
1394 52 storres
    # symbolic expression...)
1395 38 storres
    if parent(func) == parent("string"):
1396 85 storres
        localFuncSa = eval(func)
1397 85 storres
        if len(localFuncSa.variables()) > 0:
1398 85 storres
            currentVariableNameSa = localFuncSa.variables()[0]
1399 85 storres
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1400 159 storres
            functionSo = \
1401 159 storres
              sollya_lib_parse_string(localFuncSa._assume_str().replace('_SAGE_VAR_', ''))
1402 51 storres
    # Expression test.
1403 52 storres
    elif type(func) == type(zorglub):
1404 52 storres
        # Until we are able to translate Sage expressions into Sollya
1405 52 storres
        # expressions : parse the string version.
1406 85 storres
        if len(func.variables()) > 0:
1407 85 storres
            currentVariableNameSa = func.variables()[0]
1408 85 storres
            sollya_lib_name_free_variable(str(currentVariableNameSa))
1409 159 storres
            functionSo = \
1410 159 storres
              sollya_lib_parse_string(func._assume_str().replace('_SAGE_VAR_', ''))
1411 38 storres
    else:
1412 38 storres
        return(None)
1413 85 storres
    if weight is None: # No weight given -> 1.
1414 52 storres
        weightSo = pobyso_constant_1_sa_so()
1415 85 storres
    elif parent(weight) == parent("string"): # Weight given as string: parse it.
1416 51 storres
        weightSo = sollya_lib_parse_string(func)
1417 85 storres
    elif type(weight) == type(zorglub): # Weight given as symbolice expression.
1418 159 storres
        functionSo = \
1419 159 storres
          sollya_lib_parse_string_sa_so(weight._assume_str().replace('_SAGE_VAR_', ''))
1420 51 storres
    else:
1421 51 storres
        return(None)
1422 5 storres
    degreeSo = pobyso_constant_from_int(degree)
1423 85 storres
    rangeSo = pobyso_bounds_to_range_sa_so(lowerBound, upperBound)
1424 38 storres
    if not quality is None:
1425 38 storres
        qualitySo= pobyso_constant_sa_so(quality)
1426 52 storres
    else:
1427 52 storres
        qualitySo = None
1428 83 storres
1429 83 storres
    remezPolySo = sollya_lib_remez(functionSo, \
1430 83 storres
                                   degreeSo, \
1431 83 storres
                                   rangeSo, \
1432 83 storres
                                   weightSo, \
1433 83 storres
                                   qualitySo, \
1434 83 storres
                                   None)
1435 83 storres
    sollya_lib_clear_obj(functionSo)
1436 83 storres
    sollya_lib_clear_obj(degreeSo)
1437 83 storres
    sollya_lib_clear_obj(rangeSo)
1438 83 storres
    sollya_lib_clear_obj(weightSo)
1439 83 storres
    if not qualitySo is None:
1440 85 storres
        sollya_lib_clear_obj(qualitySo)
1441 83 storres
    return(remezPolySo)
1442 83 storres
# End pobyso_remez_canonical_sa_so
1443 83 storres
1444 38 storres
def pobyso_remez_canonical_so_so(funcSo, \
1445 38 storres
                                 degreeSo, \
1446 38 storres
                                 rangeSo, \
1447 52 storres
                                 weightSo = pobyso_constant_1_sa_so(),\
1448 38 storres
                                 qualitySo = None):
1449 38 storres
    """
1450 38 storres
    All arguments are pointers to Sollya objects.
1451 38 storres
    The return value is a pointer to a Sollya function.
1452 38 storres
    """
1453 38 storres
    if not sollya_lib_obj_is_function(funcSo):
1454 38 storres
        return(None)
1455 38 storres
    return(sollya_lib_remez(funcSo, degreeSo, rangeSo, weightSo, qualitySo, None))
1456 200 storres
# End pobyso_remez_canonical_so_so.
1457 200 storres
1458 215 storres
def pobyso_round_coefficients_single_so_so(polySo, precSo):
1459 215 storres
    """
1460 215 storres
    Create a rounded coefficients polynomial from polynomial argument to
1461 215 storres
    the number of bits in size argument.
1462 215 storres
    All coefficients are set to the same precision.
1463 215 storres
    """
1464 215 storres
    ## TODO: check arguments.
1465 215 storres
    endEllipListSo = pobyso_build_end_elliptic_list_so_so(precSo)
1466 215 storres
    polySo = sollya_lib_roundcoefficients(polySo, endEllipListSo, None)
1467 217 storres
    sollya_lib_clear_obj(endEllipListSo)
1468 215 storres
    #sollya_lib_clear_obj(endEllipListSo)
1469 215 storres
    return polySo
1470 215 storres
1471 215 storres
# End pobyso_round_coefficients_single_so_so
1472 215 storres
1473 5 storres
def pobyso_set_canonical_off():
1474 5 storres
    sollya_lib_set_canonical(sollya_lib_off())
1475 5 storres
1476 5 storres
def pobyso_set_canonical_on():
1477 5 storres
    sollya_lib_set_canonical(sollya_lib_on())
1478 5 storres
1479 5 storres
def pobyso_set_prec(p):
1480 38 storres
    """ Legacy function. See pobyso_set_prec_sa_so. """
1481 85 storres
    pobyso_set_prec_sa_so(p)
1482 38 storres
1483 38 storres
def pobyso_set_prec_sa_so(p):
1484 5 storres
    a = c_int(p)
1485 5 storres
    precSo = c_void_p(sollya_lib_constant_from_int(a))
1486 85 storres
    sollya_lib_set_prec(precSo, None)
1487 215 storres
# End pobyso_set_prec_sa_so.
1488 5 storres
1489 85 storres
def pobyso_set_prec_so_so(newPrecSo):
1490 85 storres
    sollya_lib_set_prec(newPrecSo, None)
1491 215 storres
# End pobyso_set_prec_so_so.
1492 54 storres
1493 215 storres
def pobyso_inf_so_so(intervalSo):
1494 215 storres
    """
1495 215 storres
    Very thin wrapper around sollya_lib_inf().
1496 215 storres
    """
1497 215 storres
    return sollya_lib_inf(intervalSo)
1498 215 storres
# End pobyso_inf_so_so.
1499 215 storres
1500 85 storres
def pobyso_supnorm_so_so(polySo, funcSo, intervalSo, errorTypeSo = None,\
1501 85 storres
                         accuracySo = None):
1502 58 storres
    """
1503 85 storres
    Computes the supnorm of the approximation error between the given
1504 85 storres
    polynomial and function.
1505 85 storres
    errorTypeSo defaults to "absolute".
1506 85 storres
    accuracySo defaults to 2^(-40).
1507 85 storres
    """
1508 85 storres
    if errorTypeSo is None:
1509 85 storres
        errorTypeSo = sollya_lib_absolute(None)
1510 85 storres
        errorTypeIsNone = True
1511 85 storres
    else:
1512 85 storres
        errorTypeIsNone = False
1513 85 storres
    #
1514 85 storres
    if accuracySo is None:
1515 85 storres
        # Notice the **!
1516 85 storres
        accuracySo = pobyso_constant_sa_so(RR(2**(-40)))
1517 85 storres
        accuracyIsNone = True
1518 85 storres
    else:
1519 85 storres
        accuracyIsNone = False
1520 85 storres
    pobyso_autoprint(accuracySo)
1521 85 storres
    resultSo = \
1522 85 storres
        sollya_lib_supnorm(polySo, funcSo, intervalSo, errorTypeSo, \
1523 85 storres
                              accuracySo)
1524 85 storres
    if errorTypeIsNone:
1525 85 storres
        sollya_lib_clear_obj(errorTypeSo)
1526 85 storres
    if accuracyIsNone:
1527 85 storres
        sollya_lib_clear_obj(accuracySo)
1528 85 storres
    return resultSo
1529 85 storres
# End pobyso_supnorm_so_so
1530 85 storres
1531 162 storres
def pobyso_taylor_expansion_no_change_var_so_so(functionSo,
1532 162 storres
                                                degreeSo,
1533 162 storres
                                                rangeSo,
1534 162 storres
                                                errorTypeSo=None,
1535 162 storres
                                                sollyaPrecSo=None):
1536 85 storres
    """
1537 162 storres
    Compute the Taylor expansion without the variable change
1538 162 storres
    x -> x-intervalCenter.
1539 58 storres
    """
1540 58 storres
    # No global change of the working precision.
1541 58 storres
    if not sollyaPrecSo is None:
1542 58 storres
        initialPrecSo = sollya_lib_get_prec(None)
1543 58 storres
        sollya_lib_set_prec(sollyaPrecSo)
1544 85 storres
    # Error type stuff: default to absolute.
1545 85 storres
    if errorTypeSo is None:
1546 85 storres
        errorTypeIsNone = True
1547 85 storres
        errorTypeSo = sollya_lib_absolute(None)
1548 85 storres
    else:
1549 85 storres
        errorTypeIsNone = False
1550 162 storres
    intervalCenterSo = sollya_lib_mid(rangeSo, None)
1551 162 storres
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo,
1552 162 storres
                                         intervalCenterSo,
1553 58 storres
                                         rangeSo, errorTypeSo, None)
1554 117 storres
    # taylorFormListSaSo is a Python list of Sollya objects references that
1555 117 storres
    # are copies of the elements of taylorFormSo.
1556 117 storres
    # pobyso_get_list_elements_so_so clears taylorFormSo.
1557 162 storres
    (taylorFormListSaSo, numElementsSa, isEndEllipticSa) = \
1558 58 storres
        pobyso_get_list_elements_so_so(taylorFormSo)
1559 162 storres
    polySo = sollya_lib_copy_obj(taylorFormListSaSo[0])
1560 162 storres
    #print "Num elements:", numElementsSa
1561 162 storres
    sollya_lib_clear_obj(taylorFormSo)
1562 162 storres
    #polySo = taylorFormListSaSo[0]
1563 162 storres
    #errorRangeSo = sollya_lib_copy_obj(taylorFormListSaSo[2])
1564 162 storres
    errorRangeSo = taylorFormListSaSo[2]
1565 181 storres
    # No copy_obj needed here: a new objects are created.
1566 181 storres
    maxErrorSo    = sollya_lib_sup(errorRangeSo)
1567 181 storres
    minErrorSo    = sollya_lib_inf(errorRangeSo)
1568 181 storres
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
1569 181 storres
    absMinErrorSo = sollya_lib_abs(minErrorSo)
1570 181 storres
    sollya_lib_clear_obj(maxErrorSo)
1571 181 storres
    sollya_lib_clear_obj(minErrorSo)
1572 181 storres
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
1573 181 storres
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
1574 58 storres
    # If changed, reset the Sollya working precision.
1575 58 storres
    if not sollyaPrecSo is None:
1576 58 storres
        sollya_lib_set_prec(initialPrecSo)
1577 83 storres
        sollya_lib_clear_obj(initialPrecSo)
1578 85 storres
    if errorTypeIsNone:
1579 85 storres
        sollya_lib_clear_obj(errorTypeSo)
1580 162 storres
    pobyso_clear_taylorform_sa_so(taylorFormListSaSo)
1581 181 storres
    if absMaxErrorSa > absMinErrorSa:
1582 181 storres
        sollya_lib_clear_obj(absMinErrorSo)
1583 181 storres
        return((polySo, intervalCenterSo, absMaxErrorSo))
1584 181 storres
    else:
1585 181 storres
        sollya_lib_clear_obj(absMaxErrorSo)
1586 181 storres
        return((polySo, intervalCenterSo, absMinErrorSo))
1587 162 storres
# end pobyso_taylor_expansion_no_change_var_so_so
1588 58 storres
1589 162 storres
def pobyso_taylor_expansion_with_change_var_so_so(functionSo, degreeSo, \
1590 162 storres
                                                  rangeSo, \
1591 162 storres
                                                  errorTypeSo=None, \
1592 162 storres
                                                  sollyaPrecSo=None):
1593 58 storres
    """
1594 162 storres
    Compute the Taylor expansion with the variable change
1595 162 storres
    x -> (x-intervalCenter) included.
1596 58 storres
    """
1597 56 storres
    # No global change of the working precision.
1598 56 storres
    if not sollyaPrecSo is None:
1599 56 storres
        initialPrecSo = sollya_lib_get_prec(None)
1600 56 storres
        sollya_lib_set_prec(sollyaPrecSo)
1601 162 storres
    #
1602 85 storres
    # Error type stuff: default to absolute.
1603 85 storres
    if errorTypeSo is None:
1604 85 storres
        errorTypeIsNone = True
1605 85 storres
        errorTypeSo = sollya_lib_absolute(None)
1606 85 storres
    else:
1607 85 storres
        errorTypeIsNone = False
1608 162 storres
    intervalCenterSo = sollya_lib_mid(rangeSo)
1609 162 storres
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo, \
1610 162 storres
                                         intervalCenterSo, \
1611 56 storres
                                         rangeSo, errorTypeSo, None)
1612 116 storres
    # taylorFormListSaSo is a Python list of Sollya objects references that
1613 116 storres
    # are copies of the elements of taylorFormSo.
1614 116 storres
    # pobyso_get_list_elements_so_so clears taylorFormSo.
1615 162 storres
    (taylorFormListSo, numElements, isEndElliptic) = \
1616 56 storres
        pobyso_get_list_elements_so_so(taylorFormSo)
1617 162 storres
    polySo = taylorFormListSo[0]
1618 162 storres
    errorRangeSo = taylorFormListSo[2]
1619 181 storres
    maxErrorSo    = sollya_lib_sup(errorRangeSo)
1620 181 storres
    minErrorSo    = sollya_lib_inf(errorRangeSo)
1621 181 storres
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
1622 181 storres
    absMinErrorSo = sollya_lib_abs(minErrorSo)
1623 181 storres
    sollya_lib_clear_obj(maxErrorSo)
1624 181 storres
    sollya_lib_clear_obj(minErrorSo)
1625 181 storres
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
1626 181 storres
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
1627 162 storres
    changeVarExpSo = sollya_lib_build_function_sub(\
1628 162 storres
                       sollya_lib_build_function_free_variable(),\
1629 162 storres
                       sollya_lib_copy_obj(intervalCenterSo))
1630 181 storres
    polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpSo)
1631 181 storres
    sollya_lib_clear_obj(polySo)
1632 162 storres
    sollya_lib_clear_obj(changeVarExpSo)
1633 56 storres
    # If changed, reset the Sollya working precision.
1634 56 storres
    if not sollyaPrecSo is None:
1635 56 storres
        sollya_lib_set_prec(initialPrecSo)
1636 63 storres
        sollya_lib_clear_obj(initialPrecSo)
1637 85 storres
    if errorTypeIsNone:
1638 85 storres
        sollya_lib_clear_obj(errorTypeSo)
1639 162 storres
    sollya_lib_clear_obj(taylorFormSo)
1640 162 storres
    # Do not clear maxErrorSo.
1641 181 storres
    if absMaxErrorSa > absMinErrorSa:
1642 181 storres
        sollya_lib_clear_obj(absMinErrorSo)
1643 181 storres
        return((polyVarChangedSo, intervalCenterSo, absMaxErrorSo))
1644 181 storres
    else:
1645 181 storres
        sollya_lib_clear_obj(absMaxErrorSo)
1646 181 storres
        return((polyVarChangedSo, intervalCenterSo, absMinErrorSo))
1647 162 storres
# end pobyso_taylor_expansion_with_change_var_so_so
1648 56 storres
1649 5 storres
def pobyso_taylor(function, degree, point):
1650 38 storres
    """ Legacy function. See pobysoTaylor_so_so. """
1651 38 storres
    return(pobyso_taylor_so_so(function, degree, point))
1652 38 storres
1653 56 storres
def pobyso_taylor_so_so(functionSo, degreeSo, pointSo):
1654 56 storres
    return(sollya_lib_taylor(functionSo, degreeSo, pointSo))
1655 5 storres
1656 85 storres
def pobyso_taylorform(function, degree, point = None,
1657 85 storres
                      interval = None, errorType=None):
1658 85 storres
    """ Legacy function. See pobyso_taylorform_sa_sa;"""
1659 38 storres
1660 38 storres
def pobyso_taylorform_sa_sa(functionSa, \
1661 84 storres
                            degreeSa, \
1662 84 storres
                            pointSa, \
1663 84 storres
                            intervalSa=None, \
1664 84 storres
                            errorTypeSa=None, \
1665 84 storres
                            precisionSa=None):
1666 37 storres
    """
1667 85 storres
    Compute the Taylor form of 'degreeSa' for 'functionSa' at 'pointSa'
1668 85 storres
    for 'intervalSa' with 'errorTypeSa' (a string) using 'precisionSa'.
1669 37 storres
    point: must be a Real or a Real interval.
1670 37 storres
    return the Taylor form as an array
1671 83 storres
    TODO: take care of the interval and of the point when it is an interval;
1672 38 storres
          when errorType is not None;
1673 83 storres
          take care of the other elements of the Taylor form (coefficients
1674 83 storres
          errors and delta.
1675 37 storres
    """
1676 37 storres
    # Absolute as the default error.
1677 84 storres
    if errorTypeSa is None:
1678 37 storres
        errorTypeSo = sollya_lib_absolute()
1679 84 storres
    elif errorTypeSa == "relative":
1680 84 storres
        errorTypeSo = sollya_lib_relative()
1681 84 storres
    elif errortypeSa == "absolute":
1682 84 storres
        errorTypeSo = sollya_lib_absolute()
1683 37 storres
    else:
1684 84 storres
        # No clean up needed.
1685 84 storres
        return None
1686 84 storres
    # Global precision stuff
1687 84 storres
    precisionChangedSa = False
1688 84 storres
    currentSollyaPrecSo = pobyso_get_prec_so()
1689 84 storres
    currentSollyaPrecSa = pobyso_constant_from_int_so_sa(currentSollyaPrecSo)
1690 84 storres
    if not precisionSa is None:
1691 84 storres
        if precisionSa > currentSollyaPrecSa:
1692 84 storres
            pobyso_set_prec_sa_so(precisionSa)
1693 84 storres
            precisionChangedSa = True
1694 84 storres
1695 85 storres
    if len(functionSa.variables()) > 0:
1696 85 storres
        varSa = functionSa.variables()[0]
1697 85 storres
        pobyso_name_free_variable_sa_so(str(varSa))
1698 84 storres
    # In any case (point or interval) the parent of pointSa has a precision
1699 84 storres
    # method.
1700 84 storres
    pointPrecSa = pointSa.parent().precision()
1701 84 storres
    if precisionSa > pointPrecSa:
1702 84 storres
        pointPrecSa = precisionSa
1703 84 storres
    # In any case (point or interval) pointSa has a base_ring() method.
1704 84 storres
    pointBaseRingString = str(pointSa.base_ring())
1705 84 storres
    if re.search('Interval', pointBaseRingString) is None: # Point
1706 84 storres
        pointSo = pobyso_constant_sa_so(pointSa, pointPrecSa)
1707 84 storres
    else: # Interval.
1708 84 storres
        pointSo = pobyso_interval_to_range_sa_so(pointSa, pointPrecSa)
1709 37 storres
    # Sollyafy the function.
1710 159 storres
    functionSo = pobyso_parse_string_sa_so(functionSa._assume_str().replace('_SAGE_VAR_', ''))
1711 37 storres
    if sollya_lib_obj_is_error(functionSo):
1712 37 storres
        print "pobyso_tailorform: function string can't be parsed!"
1713 37 storres
        return None
1714 37 storres
    # Sollyafy the degree
1715 84 storres
    degreeSo = sollya_lib_constant_from_int(int(degreeSa))
1716 37 storres
    # Sollyafy the point
1717 37 storres
    # Call Sollya
1718 83 storres
    taylorFormSo = \
1719 83 storres
        sollya_lib_taylorform(functionSo, degreeSo, pointSo, errorTypeSo,\
1720 37 storres
                                         None)
1721 85 storres
    sollya_lib_clear_obj(functionSo)
1722 85 storres
    sollya_lib_clear_obj(degreeSo)
1723 85 storres
    sollya_lib_clear_obj(pointSo)
1724 85 storres
    sollya_lib_clear_obj(errorTypeSo)
1725 38 storres
    (tfsAsList, numElements, isEndElliptic) = \
1726 38 storres
            pobyso_get_list_elements_so_so(taylorFormSo)
1727 37 storres
    polySo = tfsAsList[0]
1728 38 storres
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
1729 37 storres
    polyRealField = RealField(maxPrecision)
1730 38 storres
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, polyRealField)
1731 84 storres
    if precisionChangedSa:
1732 84 storres
        sollya_lib_set_prec(currentSollyaPrecSo)
1733 84 storres
        sollya_lib_clear_obj(currentSollyaPrecSo)
1734 37 storres
    polynomialRing = polyRealField[str(varSa)]
1735 37 storres
    polySa = polynomial(expSa, polynomialRing)
1736 37 storres
    taylorFormSa = [polySa]
1737 85 storres
    # Final clean-up.
1738 85 storres
    sollya_lib_clear_obj(taylorFormSo)
1739 51 storres
    return(taylorFormSa)
1740 51 storres
# End pobyso_taylor_form_sa_sa
1741 54 storres
1742 54 storres
def pobyso_taylorform_so_so(functionSo, degreeSo, pointSo, intervalSo=None, \
1743 54 storres
                            errorTypeSo=None):
1744 54 storres
    createdErrorType = False
1745 51 storres
    if errorTypeSo is None:
1746 51 storres
        errorTypeSo = sollya_lib_absolute()
1747 54 storres
        createdErrorType = True
1748 51 storres
    else:
1749 51 storres
        #TODO: deal with the other case.
1750 51 storres
        pass
1751 51 storres
    if intervalSo is None:
1752 54 storres
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
1753 54 storres
                                         errorTypeSo, None)
1754 51 storres
    else:
1755 54 storres
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
1756 54 storres
                                         intervalSo, errorTypeSo, None)
1757 54 storres
    if createdErrorType:
1758 54 storres
        sollya_lib_clear_obj(errorTypeSo)
1759 215 storres
    return resultSo
1760 51 storres
1761 37 storres
1762 37 storres
def pobyso_univar_polynomial_print_reverse(polySa):
1763 51 storres
    """ Legacy function. See pobyso_univar_polynomial_print_reverse_sa_sa. """
1764 51 storres
    return(pobyso_univar_polynomial_print_reverse_sa_sa(polySa))
1765 38 storres
1766 51 storres
def pobyso_univar_polynomial_print_reverse_sa_sa(polySa):
1767 37 storres
    """
1768 37 storres
    Return the string representation of a univariate polynomial with
1769 38 storres
    monomials ordered in the x^0..x^n order of the monomials.
1770 37 storres
    Remember: Sage
1771 37 storres
    """
1772 37 storres
    polynomialRing = polySa.base_ring()
1773 37 storres
    # A very expensive solution:
1774 37 storres
    # -create a fake multivariate polynomial field with only one variable,
1775 37 storres
    #   specifying a negative lexicographical order;
1776 37 storres
    mpolynomialRing = PolynomialRing(polynomialRing.base(), \
1777 37 storres
                                     polynomialRing.variable_name(), \
1778 37 storres
                                     1, order='neglex')
1779 37 storres
    # - convert the univariate argument polynomial into a multivariate
1780 37 storres
    #   version;
1781 37 storres
    p = mpolynomialRing(polySa)
1782 37 storres
    # - return the string representation of the converted form.
1783 37 storres
    # There is no simple str() method defined for p's class.
1784 37 storres
    return(p.__str__())
1785 5 storres
#
1786 5 storres
print pobyso_get_prec()
1787 5 storres
pobyso_set_prec(165)
1788 5 storres
print pobyso_get_prec()
1789 5 storres
a=100
1790 5 storres
print type(a)
1791 5 storres
id(a)
1792 5 storres
print "Max arity: ", pobyso_max_arity
1793 5 storres
print "Function tripleDouble (43) as a string: ", pobyso_function_type_as_string(43)
1794 56 storres
print "Function None (44) as a string: ", pobyso_function_type_as_string(44)
1795 56 storres
print "...Pobyso check done"