Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 232

Historique | Voir | Annoter | Télécharger (78,02 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 230 storres
#print "First constant - SOLLYA_BASE_FUNC_ABS: ", SOLLYA_BASE_FUNC_ABS
85 230 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 226 storres
    # TODO: revisit precision stuff with new technique.
277 209 storres
    if not precisionSa is None:
278 209 storres
        RRR = RealField(precisionSa)
279 209 storres
        rnArgSa = RRR(rnArgSa)
280 209 storres
    #print rnArgSa, rnArgSa.precision()
281 115 storres
    # Sollya constant creation takes place here.
282 209 storres
    return sollya_lib_constant(get_rn_value(rnArgSa))
283 115 storres
# End pobyso_constant_sa_so
284 115 storres
285 55 storres
def pobyso_constant_0_sa_so():
286 115 storres
    """
287 115 storres
    Obvious.
288 115 storres
    """
289 215 storres
    return pobyso_constant_from_int_sa_so(0)
290 55 storres
291 5 storres
def pobyso_constant_1():
292 115 storres
    """
293 115 storres
    Obvious.
294 115 storres
    Legacy function. See pobyso_constant_so_so.
295 115 storres
    """
296 215 storres
    return pobyso_constant_1_sa_so()
297 5 storres
298 52 storres
def pobyso_constant_1_sa_so():
299 115 storres
    """
300 115 storres
    Obvious.
301 115 storres
    """
302 38 storres
    return(pobyso_constant_from_int_sa_so(1))
303 38 storres
304 5 storres
def pobyso_constant_from_int(anInt):
305 38 storres
    """ Legacy function. See pobyso_constant_from_int_sa_so. """
306 215 storres
    return pobyso_constant_from_int_sa_so(anInt)
307 38 storres
308 38 storres
def pobyso_constant_from_int_sa_so(anInt):
309 115 storres
    """
310 115 storres
    Get a Sollya constant from a Sage int.
311 115 storres
    """
312 215 storres
    return sollya_lib_constant_from_int64(long(anInt))
313 5 storres
314 84 storres
def pobyso_constant_from_int_so_sa(constSo):
315 84 storres
    """
316 117 storres
    Get a Sage int from a Sollya int constant.
317 115 storres
    Usefull for precision or powers in polynomials.
318 84 storres
    """
319 200 storres
    constSa = c_long(0)
320 200 storres
    sollya_lib_get_constant_as_int64(byref(constSa), constSo)
321 215 storres
    return constSa.value
322 84 storres
# End pobyso_constant_from_int_so_sa
323 84 storres
324 209 storres
def pobyso_constant_from_mpq_sa_so(rationalSa):
325 200 storres
    """
326 200 storres
    Make a Sollya constant from Sage rational.
327 209 storres
    The Sollya constant is an unevaluated expression.
328 209 storres
    Hence no precision argument is needed.
329 209 storres
    It is better to leave this way since Sollya has its own
330 209 storres
    optimized evaluation mecanism that tries very hard to
331 209 storres
    return exact values or at least faithful ones.
332 200 storres
    """
333 200 storres
    ratExprSo = \
334 200 storres
        sollya_lib_constant_from_mpq(sgmp_get_rational_value(rationalSa))
335 209 storres
    return ratExprSo
336 200 storres
# End pobyso_constant_from_mpq_sa_so.
337 200 storres
338 209 storres
def pobyso_constant_sollya_prec_sa_so(rnArgSa):
339 209 storres
    """
340 209 storres
    Create a Sollya constant from a Sage RealNumber at the
341 209 storres
    current precision in Sollya.
342 209 storres
    """
343 209 storres
    currentSollyaPrecSa = pobyso_get_prec_so_sa()
344 209 storres
    return pobyso_constant_sa_so(rnArgSa, currentSollyaPrecSa)
345 209 storres
# End pobyso_constant_sollya_prec_sa_so
346 215 storres
347 215 storres
def pobyso_end_elliptic_list_so_sa_so(objectsListSo, intCountSa):
348 215 storres
    """
349 215 storres
    Create a Sollya end elliptic list made of the objectListSo[0] to
350 215 storres
     objectsListSo[intCountSa-1] objects.
351 215 storres
    """
352 215 storres
    return sollya_lib_end_elliptic_list(objectSo, int(intCountSa))
353 215 storres
354 155 storres
def pobyso_error_so():
355 155 storres
    return sollya_lib_error(None)
356 155 storres
# End pobyso_error().
357 155 storres
358 215 storres
def pobyso_evaluate_so_so(funcSo, argumentSo):
359 215 storres
    """
360 215 storres
    Evaluates funcSo for arguemntSo through sollya_lib_evaluate().
361 215 storres
    """
362 215 storres
    return sollya_lib_evaluate(funcSo, argumentSo)
363 215 storres
# End pobyso_evaluate_so_so.
364 215 storres
365 209 storres
def pobyso_float_poly_sa_so(polySa, precSa = None):
366 209 storres
    """
367 209 storres
    Create a Sollya polynomial from a Sage RealField polynomial.
368 209 storres
    """
369 209 storres
    ## TODO: filter arguments.
370 209 storres
    ## Precision. If a precision is given, convert the polynomial
371 209 storres
    #  into the right polynomial field. If not convert it straight
372 209 storres
    #  to Sollya.
373 218 storres
    sollyaPrecChanged = False
374 226 storres
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
375 218 storres
    if precSa is None:
376 209 storres
        precSa = polySa.parent().base_ring().precision()
377 226 storres
    if (precSa > initialSollyaPrecSa):
378 228 storres
        assert precSa >= 2, "Precision change <2 requested"
379 226 storres
        if precSa <= 2:
380 226 storres
            print inspect.stack()[0][3], ": precision change <= 2 requested"
381 218 storres
        precSo = pobyso_constant_from_int(precSa)
382 218 storres
        pobyso_set_prec_so_so(precSo)
383 218 storres
        sollya_lib_clear_obj(precSo)
384 218 storres
        sollyaPrecChanged = True
385 209 storres
    ## Get exponents and coefficients.
386 218 storres
    exponentsSa     = polySa.exponents()
387 218 storres
    coefficientsSa  = polySa.coefficients()
388 209 storres
    ## Build the polynomial.
389 209 storres
    polySo = None
390 213 storres
    for coefficientSa, exponentSa in zip(coefficientsSa, exponentsSa):
391 209 storres
        #print coefficientSa.n(prec=precSa), exponentSa
392 209 storres
        coefficientSo = \
393 209 storres
            pobyso_constant_sa_so(coefficientSa)
394 209 storres
        #pobyso_autoprint(coefficientSo)
395 209 storres
        exponentSo = \
396 209 storres
            pobyso_constant_from_int_sa_so(exponentSa)
397 209 storres
        #pobyso_autoprint(exponentSo)
398 209 storres
        monomialSo = sollya_lib_build_function_pow(
399 209 storres
                       sollya_lib_build_function_free_variable(),
400 209 storres
                       exponentSo)
401 218 storres
        polyTermSo = sollya_lib_build_function_mul(coefficientSo,
402 218 storres
                                                       monomialSo)
403 209 storres
        if polySo is None:
404 218 storres
            polySo = polyTermSo
405 209 storres
        else:
406 209 storres
            polySo = sollya_lib_build_function_add(polySo, polyTermSo)
407 218 storres
    if sollyaPrecChanged:
408 226 storres
        pobyso_set_prec_so_so(initialSollyaPrecSo)
409 226 storres
        sollya_lib_clear_obj(initialSollyaPrecSo)
410 209 storres
    return polySo
411 209 storres
# End pobyso_float_poly_sa_so
412 209 storres
413 209 storres
def pobyso_float_poly_so_sa(polySo, realFieldSa=None):
414 209 storres
    """
415 209 storres
    Convert a Sollya polynomial into a Sage floating-point polynomial.
416 209 storres
    If no realField is given, a RealField corresponding to the maximum
417 209 storres
    precision of the coefficients is internally computed.
418 209 storres
    The real field is not returned but can be easily retrieved from
419 209 storres
    the polynomial itself.
420 209 storres
    ALGORITHM:
421 209 storres
    - (optional) compute the RealField of the coefficients;
422 209 storres
    - convert the Sollya expression into a Sage expression;
423 209 storres
    - convert the Sage expression into a Sage polynomial
424 209 storres
    """
425 209 storres
    if realFieldSa is None:
426 209 storres
        expressionPrecSa = pobyso_get_max_prec_of_exp_so_sa(polySo)
427 218 storres
        #print "Maximum precision of Sollya polynomial coefficients:", expressionPrecSa
428 226 storres
        if expressionPrecSa < 2 or expressionPrecSa > 2147483391:
429 226 storres
            print "Maximum degree of expression:", expressionPrecSa
430 209 storres
        realFieldSa      = RealField(expressionPrecSa)
431 209 storres
    #print "Sollya expression before...",
432 209 storres
    #pobyso_autoprint(polySo)
433 209 storres
434 209 storres
    expressionSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo,
435 209 storres
                                                             realFieldSa)
436 218 storres
    #print "...Sollya expression after."
437 209 storres
    #pobyso_autoprint(polySo)
438 209 storres
    polyVariableSa = expressionSa.variables()[0]
439 209 storres
    polyRingSa     = realFieldSa[str(polyVariableSa)]
440 209 storres
    #print polyRingSa
441 209 storres
    # Do not use the polynomial(expressionSa, ring=polyRingSa) form!
442 209 storres
    polynomialSa = polyRingSa(expressionSa)
443 215 storres
    polyCoeffsListSa = polynomialSa.coefficients()
444 215 storres
    #for coeff in polyCoeffsListSa:
445 215 storres
    #    print coeff.abs().n()
446 209 storres
    return polynomialSa
447 209 storres
# End pobyso_float_poly_so_sa
448 209 storres
449 215 storres
def pobyso_free_variable():
450 215 storres
    """
451 215 storres
    Ultra thin wrapper around the sollya_lib_function_build_free_variable function.
452 215 storres
    """
453 215 storres
    return sollya_lib_build_function_free_variable()
454 209 storres
455 5 storres
def pobyso_function_type_as_string(funcType):
456 38 storres
    """ Legacy function. See pobyso_function_type_as_string_so_sa. """
457 38 storres
    return(pobyso_function_type_as_string_so_sa(funcType))
458 38 storres
459 38 storres
def pobyso_function_type_as_string_so_sa(funcType):
460 38 storres
    """
461 38 storres
    Numeric Sollya function codes -> Sage mathematical function names.
462 38 storres
    Notice that pow -> ^ (a la Sage, not a la Python).
463 38 storres
    """
464 5 storres
    if funcType == SOLLYA_BASE_FUNC_ABS:
465 5 storres
        return "abs"
466 5 storres
    elif funcType == SOLLYA_BASE_FUNC_ACOS:
467 5 storres
        return "arccos"
468 5 storres
    elif funcType == SOLLYA_BASE_FUNC_ACOSH:
469 5 storres
        return "arccosh"
470 5 storres
    elif funcType == SOLLYA_BASE_FUNC_ADD:
471 5 storres
        return "+"
472 5 storres
    elif funcType == SOLLYA_BASE_FUNC_ASIN:
473 5 storres
        return "arcsin"
474 5 storres
    elif funcType == SOLLYA_BASE_FUNC_ASINH:
475 5 storres
        return "arcsinh"
476 5 storres
    elif funcType == SOLLYA_BASE_FUNC_ATAN:
477 5 storres
        return "arctan"
478 5 storres
    elif funcType == SOLLYA_BASE_FUNC_ATANH:
479 5 storres
        return "arctanh"
480 5 storres
    elif funcType == SOLLYA_BASE_FUNC_CEIL:
481 5 storres
        return "ceil"
482 5 storres
    elif funcType == SOLLYA_BASE_FUNC_CONSTANT:
483 5 storres
        return "cte"
484 5 storres
    elif funcType == SOLLYA_BASE_FUNC_COS:
485 5 storres
        return "cos"
486 5 storres
    elif funcType == SOLLYA_BASE_FUNC_COSH:
487 5 storres
        return "cosh"
488 5 storres
    elif funcType == SOLLYA_BASE_FUNC_DIV:
489 5 storres
        return "/"
490 5 storres
    elif funcType == SOLLYA_BASE_FUNC_DOUBLE:
491 5 storres
        return "double"
492 5 storres
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEDOUBLE:
493 5 storres
        return "doubleDouble"
494 5 storres
    elif funcType == SOLLYA_BASE_FUNC_DOUBLEEXTENDED:
495 5 storres
        return "doubleDxtended"
496 5 storres
    elif funcType == SOLLYA_BASE_FUNC_ERF:
497 5 storres
        return "erf"
498 5 storres
    elif funcType == SOLLYA_BASE_FUNC_ERFC:
499 5 storres
        return "erfc"
500 5 storres
    elif funcType == SOLLYA_BASE_FUNC_EXP:
501 5 storres
        return "exp"
502 5 storres
    elif funcType == SOLLYA_BASE_FUNC_EXP_M1:
503 5 storres
        return "expm1"
504 5 storres
    elif funcType == SOLLYA_BASE_FUNC_FLOOR:
505 5 storres
        return "floor"
506 5 storres
    elif funcType == SOLLYA_BASE_FUNC_FREE_VARIABLE:
507 5 storres
        return "freeVariable"
508 5 storres
    elif funcType == SOLLYA_BASE_FUNC_HALFPRECISION:
509 5 storres
        return "halfPrecision"
510 5 storres
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYCONSTANT:
511 5 storres
        return "libraryConstant"
512 5 storres
    elif funcType == SOLLYA_BASE_FUNC_LIBRARYFUNCTION:
513 5 storres
        return "libraryFunction"
514 5 storres
    elif funcType == SOLLYA_BASE_FUNC_LOG:
515 5 storres
        return "log"
516 5 storres
    elif funcType == SOLLYA_BASE_FUNC_LOG_10:
517 5 storres
        return "log10"
518 5 storres
    elif funcType == SOLLYA_BASE_FUNC_LOG_1P:
519 5 storres
        return "log1p"
520 5 storres
    elif funcType == SOLLYA_BASE_FUNC_LOG_2:
521 5 storres
        return "log2"
522 5 storres
    elif funcType == SOLLYA_BASE_FUNC_MUL:
523 5 storres
        return "*"
524 5 storres
    elif funcType == SOLLYA_BASE_FUNC_NEARESTINT:
525 5 storres
        return "round"
526 5 storres
    elif funcType == SOLLYA_BASE_FUNC_NEG:
527 5 storres
        return "__neg__"
528 5 storres
    elif funcType == SOLLYA_BASE_FUNC_PI:
529 5 storres
        return "pi"
530 5 storres
    elif funcType == SOLLYA_BASE_FUNC_POW:
531 5 storres
        return "^"
532 5 storres
    elif funcType == SOLLYA_BASE_FUNC_PROCEDUREFUNCTION:
533 5 storres
        return "procedureFunction"
534 5 storres
    elif funcType == SOLLYA_BASE_FUNC_QUAD:
535 5 storres
        return "quad"
536 5 storres
    elif funcType == SOLLYA_BASE_FUNC_SIN:
537 5 storres
        return "sin"
538 5 storres
    elif funcType == SOLLYA_BASE_FUNC_SINGLE:
539 5 storres
        return "single"
540 5 storres
    elif funcType == SOLLYA_BASE_FUNC_SINH:
541 5 storres
        return "sinh"
542 5 storres
    elif funcType == SOLLYA_BASE_FUNC_SQRT:
543 5 storres
        return "sqrt"
544 5 storres
    elif funcType == SOLLYA_BASE_FUNC_SUB:
545 5 storres
        return "-"
546 5 storres
    elif funcType == SOLLYA_BASE_FUNC_TAN:
547 5 storres
        return "tan"
548 5 storres
    elif funcType == SOLLYA_BASE_FUNC_TANH:
549 5 storres
        return "tanh"
550 5 storres
    elif funcType == SOLLYA_BASE_FUNC_TRIPLEDOUBLE:
551 5 storres
        return "tripleDouble"
552 5 storres
    else:
553 5 storres
        return None
554 5 storres
555 85 storres
def pobyso_get_constant(rnArgSa, constSo):
556 38 storres
    """ Legacy function. See pobyso_get_constant_so_sa. """
557 209 storres
    return pobyso_get_constant_so_sa(rnArgSa, constSo)
558 209 storres
# End pobyso_get_constant
559 38 storres
560 84 storres
def pobyso_get_constant_so_sa(rnArgSa, constSo):
561 52 storres
    """
562 85 storres
    Set the value of rnArgSo to the value of constSo in MPFR_RNDN mode.
563 52 storres
    rnArg must already exist and belong to some RealField.
564 85 storres
    We assume that constSo points to a Sollya constant.
565 52 storres
    """
566 209 storres
    outcome = sollya_lib_get_constant(get_rn_value(rnArgSa), constSo)
567 209 storres
    if outcome == 0: # Failure because constSo is not a constant expression.
568 209 storres
        return None
569 209 storres
    else:
570 209 storres
        return outcome
571 209 storres
# End  pobyso_get_constant_so_sa
572 209 storres
573 57 storres
def pobyso_get_constant_as_rn(ctExpSo):
574 83 storres
    """
575 83 storres
    Legacy function. See pobyso_get_constant_as_rn_so_sa.
576 83 storres
    """
577 57 storres
    return(pobyso_get_constant_as_rn_so_sa(ctExpSo))
578 38 storres
579 56 storres
def pobyso_get_constant_as_rn_so_sa(constExpSo):
580 83 storres
    """
581 83 storres
    Get a Sollya constant as a Sage "real number".
582 83 storres
    The precision of the floating-point number returned is that of the Sollya
583 83 storres
    constant.
584 83 storres
    """
585 218 storres
    #print "Before computing precision of variable..."
586 218 storres
    #pobyso_autoprint(constExpSo)
587 209 storres
    precisionSa  = pobyso_get_prec_of_constant_so_sa(constExpSo)
588 218 storres
    #print "precisionSa:", precisionSa
589 209 storres
    ## If the expression can not be exactly converted, None is returned.
590 209 storres
    #  In this case opt for the Sollya current expression.
591 209 storres
    if precisionSa is None:
592 209 storres
        precisionSa = pobyso_get_prec_so_sa()
593 56 storres
    RRRR = RealField(precisionSa)
594 56 storres
    rnSa = RRRR(0)
595 209 storres
    outcome = sollya_lib_get_constant(get_rn_value(rnSa), constExpSo)
596 209 storres
    if outcome == 0:
597 209 storres
        return None
598 209 storres
    else:
599 209 storres
        return rnSa
600 83 storres
# End pobyso_get_constant_as_rn_so_sa
601 38 storres
602 38 storres
def pobyso_get_constant_as_rn_with_rf(ctExp, realField):
603 83 storres
    """
604 83 storres
    Legacy function. See pobyso_get_constant_as_rn_with_rf_so_sa.
605 83 storres
    """
606 209 storres
    return pobyso_get_constant_as_rn_with_rf_so_sa(ctExp, realField)
607 209 storres
# End pobyso_get_constant_as_rn_with_rf
608 5 storres
609 56 storres
def pobyso_get_constant_as_rn_with_rf_so_sa(ctExpSo, realFieldSa = None):
610 83 storres
    """
611 83 storres
    Get a Sollya constant as a Sage "real number".
612 83 storres
    If no real field is specified, the precision of the floating-point number
613 85 storres
    returned is that of the Sollya constant.
614 83 storres
    Otherwise is is that of the real field. Hence rounding may happen.
615 83 storres
    """
616 56 storres
    if realFieldSa is None:
617 209 storres
        return pobyso_get_constant_as_rn_so_sa(ctExpSo)
618 56 storres
    rnSa = realFieldSa(0)
619 209 storres
    outcome = sollya_lib_get_constant(get_rn_value(rnSa), ctExpSo)
620 209 storres
    if outcome == 0:
621 209 storres
        return None
622 209 storres
    else:
623 209 storres
        return rnSa
624 83 storres
# End pobyso_get_constant_as_rn_with_rf_so_sa
625 38 storres
626 5 storres
def pobyso_get_free_variable_name():
627 83 storres
    """
628 83 storres
    Legacy function. See pobyso_get_free_variable_name_so_sa.
629 83 storres
    """
630 38 storres
    return(pobyso_get_free_variable_name_so_sa())
631 38 storres
632 38 storres
def pobyso_get_free_variable_name_so_sa():
633 209 storres
    return sollya_lib_get_free_variable_name()
634 5 storres
635 38 storres
def pobyso_get_function_arity(expressionSo):
636 83 storres
    """
637 83 storres
    Legacy function. See pobyso_get_function_arity_so_sa.
638 83 storres
    """
639 38 storres
    return(pobyso_get_function_arity_so_sa(expressionSo))
640 38 storres
641 38 storres
def pobyso_get_function_arity_so_sa(expressionSo):
642 5 storres
    arity = c_int(0)
643 38 storres
    sollya_lib_get_function_arity(byref(arity),expressionSo)
644 209 storres
    return int(arity.value)
645 5 storres
646 38 storres
def pobyso_get_head_function(expressionSo):
647 83 storres
    """
648 83 storres
    Legacy function. See pobyso_get_head_function_so_sa.
649 83 storres
    """
650 38 storres
    return(pobyso_get_head_function_so_sa(expressionSo))
651 38 storres
652 38 storres
def pobyso_get_head_function_so_sa(expressionSo):
653 5 storres
    functionType = c_int(0)
654 218 storres
    sollya_lib_get_head_function(byref(functionType), expressionSo)
655 209 storres
    return int(functionType.value)
656 5 storres
657 56 storres
def pobyso_get_interval_from_range_so_sa(soRange, realIntervalFieldSa = None ):
658 53 storres
    """
659 53 storres
    Return the Sage interval corresponding to the Sollya range argument.
660 83 storres
    If no reaIntervalField is passed as an argument, the interval bounds are not
661 56 storres
    rounded: they are elements of RealIntervalField of the "right" precision
662 56 storres
    to hold all the digits.
663 53 storres
    """
664 53 storres
    prec = c_int(0)
665 56 storres
    if realIntervalFieldSa is None:
666 56 storres
        retval = sollya_lib_get_prec_of_range(byref(prec), soRange, None)
667 56 storres
        if retval == 0:
668 209 storres
            return None
669 56 storres
        realIntervalFieldSa = RealIntervalField(prec.value)
670 56 storres
    intervalSa = realIntervalFieldSa(0,0)
671 53 storres
    retval = \
672 53 storres
        sollya_lib_get_interval_from_range(get_interval_value(intervalSa),\
673 53 storres
                                           soRange)
674 53 storres
    if retval == 0:
675 209 storres
        return None
676 209 storres
    return intervalSa
677 56 storres
# End pobyso_get_interval_from_range_so_sa
678 56 storres
679 5 storres
def pobyso_get_list_elements(soObj):
680 38 storres
    """ Legacy function. See pobyso_get_list_elements_so_so. """
681 209 storres
    return pobyso_get_list_elements_so_so(soObj)
682 38 storres
683 117 storres
def pobyso_get_list_elements_so_so(objectListSo):
684 51 storres
    """
685 118 storres
    Get the Sollya list elements as a Sage/Python array of Sollya objects.
686 118 storres

687 118 storres
    INPUT:
688 118 storres
    - objectListSo: a Sollya list of Sollya objects.
689 118 storres

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

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