Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 230

Historique | Voir | Annoter | Télécharger (77,65 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 223 storres
        return (polySo, currentApproxErrorSo)
1565 223 storres
    degreeSa             = pobyso_polynomial_degree_so_sa(polySo)
1566 223 storres
    intervalSa           = pobyso_range_to_interval_so_sa(intervalSo)
1567 223 storres
1568 223 storres
    if debug:
1569 228 storres
        print "degreeSa              :", degreeSa
1570 228 storres
        print "intervalSa            :", intervalSa.str(style='brackets')
1571 223 storres
        print "currentApproxErrorSa  :", currentApproxErrorSa
1572 228 storres
        print "approxAccurSa         :", approxAccurSa
1573 223 storres
    ### Start with a 0 value expression.
1574 223 storres
    radiusSa = intervalSa.absolute_diameter() / 2
1575 223 storres
    if debug:
1576 224 storres
        print "log2(radius):", RR(radiusSa).log2()
1577 224 storres
    iterIndex = 0
1578 224 storres
    while True:
1579 224 storres
        resPolySo = pobyso_constant_0_sa_so()
1580 224 storres
        roundedPolyApproxAccurSa = approxAccurSa / 2
1581 224 storres
        currentRadiusPowerSa = 1
1582 224 storres
        for degree in xrange(0,degreeSa + 1):
1583 224 storres
            #### At round 0, use the agressive formula. At round 1, run the
1584 224 storres
            #    proved formula.
1585 224 storres
            if iterIndex == 0:
1586 224 storres
                roundingPowerSa = \
1587 224 storres
                    floor(((currentRadiusPowerSa/roundedPolyApproxAccurSa)*(degree+1)).log2())
1588 224 storres
            else:
1589 224 storres
                roundingPowerSa = \
1590 224 storres
                    floor(((currentRadiusPowerSa/roundedPolyApproxAccurSa)*(degreeSa+1)).log2())
1591 228 storres
            ## Under extreme conditions the above formulas can evaluate under 2, which is the
1592 228 storres
            #  minimal precision of an MPFR number.
1593 228 storres
            if roundingPowerSa < 2:
1594 228 storres
                roundingPowerSa = 2
1595 224 storres
            if debug:
1596 224 storres
                print "roundedPolyApproxAccurSa", roundedPolyApproxAccurSa
1597 224 storres
                print "currentRadiusPowerSa", currentRadiusPowerSa
1598 224 storres
                print "Current rounding exponent:", roundingPowerSa
1599 224 storres
            currentRadiusPowerSa *= radiusSa
1600 224 storres
            index1So = pobyso_constant_from_int_sa_so(degree)
1601 224 storres
            index2So = pobyso_constant_from_int_sa_so(degree)
1602 224 storres
            ### Create a monomial with:
1603 224 storres
            #   - the coefficient in the initial monomial at the current degrree;
1604 224 storres
            #   - the current exponent;
1605 224 storres
            #   - the free variable.
1606 224 storres
            cmonSo  = \
1607 224 storres
                sollya_lib_build_function_mul(sollya_lib_coeff(polySo, index1So),
1608 224 storres
                                              sollya_lib_build_function_pow( \
1609 224 storres
                                                sollya_lib_build_function_free_variable(), \
1610 224 storres
                                                index2So))
1611 224 storres
            roundingPowerSo = pobyso_constant_from_int_sa_so(roundingPowerSa)
1612 224 storres
            cmonrSo = pobyso_round_coefficients_single_so_so(cmonSo, roundingPowerSo)
1613 224 storres
            sollya_lib_clear_obj(cmonSo)
1614 224 storres
            ### Add to the result polynomial.
1615 224 storres
            resPolySo = sollya_lib_build_function_add(resPolySo,
1616 224 storres
                                                      cmonrSo)
1617 224 storres
        # End for.
1618 224 storres
        ### Check the new polynomial.
1619 224 storres
        freeVarSo     = sollya_lib_build_function_free_variable()
1620 224 storres
        changeVarSo   = sollya_lib_sub(freeVarSo, icSo)
1621 224 storres
        resPolyCvSo   = sollya_lib_evaluate(resPolySo, changeVarSo)
1622 224 storres
        errFuncSo     = sollya_lib_build_function_sub(sollya_lib_copy_obj(funcSo),
1623 224 storres
                                                      resPolyCvSo)
1624 224 storres
        infNormSo = sollya_lib_dirtyinfnorm(errFuncSo, intervalSo)
1625 224 storres
        ### This also clears resPolyCvSo.
1626 224 storres
        sollya_lib_clear_obj(errFuncSo)
1627 224 storres
        cerrSa    = pobyso_get_constant_as_rn_so_sa(infNormSo)
1628 223 storres
        if debug:
1629 224 storres
            print "Error of the new polynomial:", cerrSa
1630 224 storres
        ### If at round 1, return the initial polynomial error. This should
1631 224 storres
        #   never happen since the rounding algorithm is proved. But some
1632 224 storres
        #   circumstances may break it (e.g. internal precision of tools).
1633 224 storres
        if cerrSa > approxAccurSa:
1634 224 storres
            if iterIndex > 0: # Round 1 and beyond.
1635 224 storres
                sollya_lib_clear_obj(resPolySo)
1636 224 storres
                sollya_lib_clear_obj(infNormSo)
1637 224 storres
                return (polySo, currentApproxErrorSo)
1638 224 storres
            else: # Round 0, got round 1
1639 224 storres
                sollya_lib_clear_obj(resPolySo)
1640 224 storres
                sollya_lib_clear_obj(infNormSo)
1641 224 storres
                iterIndex += 1
1642 224 storres
                continue
1643 224 storres
        ### If get here it is because cerrSa <= approxAccurSa
1644 224 storres
        ### Approximation error of the new polynomial is acceptable.
1645 224 storres
        return (resPolySo, infNormSo)
1646 224 storres
    # End while True
1647 224 storres
# End pobyso_round_coefficients_progressive_so_so
1648 223 storres
1649 215 storres
def pobyso_round_coefficients_single_so_so(polySo, precSo):
1650 215 storres
    """
1651 215 storres
    Create a rounded coefficients polynomial from polynomial argument to
1652 215 storres
    the number of bits in size argument.
1653 215 storres
    All coefficients are set to the same precision.
1654 215 storres
    """
1655 215 storres
    ## TODO: check arguments.
1656 215 storres
    endEllipListSo = pobyso_build_end_elliptic_list_so_so(precSo)
1657 215 storres
    polySo = sollya_lib_roundcoefficients(polySo, endEllipListSo, None)
1658 217 storres
    sollya_lib_clear_obj(endEllipListSo)
1659 215 storres
    #sollya_lib_clear_obj(endEllipListSo)
1660 215 storres
    return polySo
1661 215 storres
1662 215 storres
# End pobyso_round_coefficients_single_so_so
1663 215 storres
1664 5 storres
def pobyso_set_canonical_off():
1665 5 storres
    sollya_lib_set_canonical(sollya_lib_off())
1666 5 storres
1667 5 storres
def pobyso_set_canonical_on():
1668 5 storres
    sollya_lib_set_canonical(sollya_lib_on())
1669 5 storres
1670 5 storres
def pobyso_set_prec(p):
1671 38 storres
    """ Legacy function. See pobyso_set_prec_sa_so. """
1672 85 storres
    pobyso_set_prec_sa_so(p)
1673 38 storres
1674 38 storres
def pobyso_set_prec_sa_so(p):
1675 227 storres
    #a = c_int(p)
1676 226 storres
    #precSo = c_void_p(sollya_lib_constant_from_int(a))
1677 227 storres
    #precSo = sollya_lib_constant_from_int(a)
1678 227 storres
    precSo = pobyso_constant_from_int_sa_so(p)
1679 226 storres
    sollya_lib_set_prec(precSo)
1680 226 storres
    sollya_lib_clear_obj(precSo)
1681 215 storres
# End pobyso_set_prec_sa_so.
1682 5 storres
1683 85 storres
def pobyso_set_prec_so_so(newPrecSo):
1684 226 storres
    sollya_lib_set_prec(newPrecSo)
1685 215 storres
# End pobyso_set_prec_so_so.
1686 54 storres
1687 215 storres
def pobyso_inf_so_so(intervalSo):
1688 215 storres
    """
1689 215 storres
    Very thin wrapper around sollya_lib_inf().
1690 215 storres
    """
1691 215 storres
    return sollya_lib_inf(intervalSo)
1692 215 storres
# End pobyso_inf_so_so.
1693 215 storres
1694 85 storres
def pobyso_supnorm_so_so(polySo, funcSo, intervalSo, errorTypeSo = None,\
1695 85 storres
                         accuracySo = None):
1696 58 storres
    """
1697 85 storres
    Computes the supnorm of the approximation error between the given
1698 85 storres
    polynomial and function.
1699 85 storres
    errorTypeSo defaults to "absolute".
1700 85 storres
    accuracySo defaults to 2^(-40).
1701 85 storres
    """
1702 85 storres
    if errorTypeSo is None:
1703 85 storres
        errorTypeSo = sollya_lib_absolute(None)
1704 85 storres
        errorTypeIsNone = True
1705 85 storres
    else:
1706 85 storres
        errorTypeIsNone = False
1707 85 storres
    #
1708 85 storres
    if accuracySo is None:
1709 228 storres
        # Notice the **: we are in Pythonland!
1710 85 storres
        accuracySo = pobyso_constant_sa_so(RR(2**(-40)))
1711 85 storres
        accuracyIsNone = True
1712 85 storres
    else:
1713 85 storres
        accuracyIsNone = False
1714 85 storres
    pobyso_autoprint(accuracySo)
1715 85 storres
    resultSo = \
1716 85 storres
        sollya_lib_supnorm(polySo, funcSo, intervalSo, errorTypeSo, \
1717 85 storres
                              accuracySo)
1718 85 storres
    if errorTypeIsNone:
1719 85 storres
        sollya_lib_clear_obj(errorTypeSo)
1720 85 storres
    if accuracyIsNone:
1721 85 storres
        sollya_lib_clear_obj(accuracySo)
1722 85 storres
    return resultSo
1723 85 storres
# End pobyso_supnorm_so_so
1724 85 storres
1725 162 storres
def pobyso_taylor_expansion_no_change_var_so_so(functionSo,
1726 162 storres
                                                degreeSo,
1727 162 storres
                                                rangeSo,
1728 162 storres
                                                errorTypeSo=None,
1729 162 storres
                                                sollyaPrecSo=None):
1730 85 storres
    """
1731 162 storres
    Compute the Taylor expansion without the variable change
1732 162 storres
    x -> x-intervalCenter.
1733 58 storres
    """
1734 226 storres
    # Change internal Sollya precision, if needed.
1735 227 storres
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
1736 226 storres
    sollyaPrecChanged = False
1737 226 storres
    if sollyaPrecSo is None:
1738 226 storres
        pass
1739 226 storres
    else:
1740 58 storres
        sollya_lib_set_prec(sollyaPrecSo)
1741 226 storres
        sollyaPrecChanged = True
1742 85 storres
    # Error type stuff: default to absolute.
1743 85 storres
    if errorTypeSo is None:
1744 85 storres
        errorTypeIsNone = True
1745 85 storres
        errorTypeSo = sollya_lib_absolute(None)
1746 85 storres
    else:
1747 85 storres
        errorTypeIsNone = False
1748 162 storres
    intervalCenterSo = sollya_lib_mid(rangeSo, None)
1749 162 storres
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo,
1750 162 storres
                                         intervalCenterSo,
1751 58 storres
                                         rangeSo, errorTypeSo, None)
1752 117 storres
    # taylorFormListSaSo is a Python list of Sollya objects references that
1753 117 storres
    # are copies of the elements of taylorFormSo.
1754 117 storres
    # pobyso_get_list_elements_so_so clears taylorFormSo.
1755 162 storres
    (taylorFormListSaSo, numElementsSa, isEndEllipticSa) = \
1756 58 storres
        pobyso_get_list_elements_so_so(taylorFormSo)
1757 162 storres
    polySo = sollya_lib_copy_obj(taylorFormListSaSo[0])
1758 162 storres
    #print "Num elements:", numElementsSa
1759 162 storres
    sollya_lib_clear_obj(taylorFormSo)
1760 162 storres
    #polySo = taylorFormListSaSo[0]
1761 162 storres
    #errorRangeSo = sollya_lib_copy_obj(taylorFormListSaSo[2])
1762 162 storres
    errorRangeSo = taylorFormListSaSo[2]
1763 181 storres
    # No copy_obj needed here: a new objects are created.
1764 181 storres
    maxErrorSo    = sollya_lib_sup(errorRangeSo)
1765 181 storres
    minErrorSo    = sollya_lib_inf(errorRangeSo)
1766 181 storres
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
1767 181 storres
    absMinErrorSo = sollya_lib_abs(minErrorSo)
1768 181 storres
    sollya_lib_clear_obj(maxErrorSo)
1769 181 storres
    sollya_lib_clear_obj(minErrorSo)
1770 181 storres
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
1771 181 storres
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
1772 58 storres
    # If changed, reset the Sollya working precision.
1773 226 storres
    if sollyaPrecChanged:
1774 226 storres
        sollya_lib_set_prec(initialSollyaPrecSo)
1775 226 storres
    sollya_lib_clear_obj(initialSollyaPrecSo)
1776 85 storres
    if errorTypeIsNone:
1777 85 storres
        sollya_lib_clear_obj(errorTypeSo)
1778 162 storres
    pobyso_clear_taylorform_sa_so(taylorFormListSaSo)
1779 181 storres
    if absMaxErrorSa > absMinErrorSa:
1780 181 storres
        sollya_lib_clear_obj(absMinErrorSo)
1781 227 storres
        return (polySo, intervalCenterSo, absMaxErrorSo)
1782 181 storres
    else:
1783 181 storres
        sollya_lib_clear_obj(absMaxErrorSo)
1784 227 storres
        return (polySo, intervalCenterSo, absMinErrorSo)
1785 162 storres
# end pobyso_taylor_expansion_no_change_var_so_so
1786 58 storres
1787 162 storres
def pobyso_taylor_expansion_with_change_var_so_so(functionSo, degreeSo, \
1788 162 storres
                                                  rangeSo, \
1789 162 storres
                                                  errorTypeSo=None, \
1790 162 storres
                                                  sollyaPrecSo=None):
1791 58 storres
    """
1792 162 storres
    Compute the Taylor expansion with the variable change
1793 162 storres
    x -> (x-intervalCenter) included.
1794 58 storres
    """
1795 226 storres
    # Change Sollya internal precision, if need.
1796 226 storres
    sollyaPrecChanged = False
1797 226 storres
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_sos_sa()
1798 226 storres
    if sollyaPrecSo is None:
1799 226 storres
        pass
1800 226 storres
    else:
1801 56 storres
        sollya_lib_set_prec(sollyaPrecSo)
1802 226 storres
        sollyaPrecChanged = True
1803 162 storres
    #
1804 85 storres
    # Error type stuff: default to absolute.
1805 85 storres
    if errorTypeSo is None:
1806 85 storres
        errorTypeIsNone = True
1807 85 storres
        errorTypeSo = sollya_lib_absolute(None)
1808 85 storres
    else:
1809 85 storres
        errorTypeIsNone = False
1810 162 storres
    intervalCenterSo = sollya_lib_mid(rangeSo)
1811 162 storres
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo, \
1812 162 storres
                                         intervalCenterSo, \
1813 56 storres
                                         rangeSo, errorTypeSo, None)
1814 116 storres
    # taylorFormListSaSo is a Python list of Sollya objects references that
1815 116 storres
    # are copies of the elements of taylorFormSo.
1816 116 storres
    # pobyso_get_list_elements_so_so clears taylorFormSo.
1817 162 storres
    (taylorFormListSo, numElements, isEndElliptic) = \
1818 56 storres
        pobyso_get_list_elements_so_so(taylorFormSo)
1819 162 storres
    polySo = taylorFormListSo[0]
1820 162 storres
    errorRangeSo = taylorFormListSo[2]
1821 181 storres
    maxErrorSo    = sollya_lib_sup(errorRangeSo)
1822 181 storres
    minErrorSo    = sollya_lib_inf(errorRangeSo)
1823 181 storres
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
1824 181 storres
    absMinErrorSo = sollya_lib_abs(minErrorSo)
1825 181 storres
    sollya_lib_clear_obj(maxErrorSo)
1826 181 storres
    sollya_lib_clear_obj(minErrorSo)
1827 181 storres
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
1828 181 storres
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
1829 162 storres
    changeVarExpSo = sollya_lib_build_function_sub(\
1830 162 storres
                       sollya_lib_build_function_free_variable(),\
1831 162 storres
                       sollya_lib_copy_obj(intervalCenterSo))
1832 181 storres
    polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpSo)
1833 181 storres
    sollya_lib_clear_obj(polySo)
1834 162 storres
    sollya_lib_clear_obj(changeVarExpSo)
1835 56 storres
    # If changed, reset the Sollya working precision.
1836 226 storres
    if sollyaPrecChanged:
1837 226 storres
        sollya_lib_set_prec(initialSollyaPrecSo)
1838 226 storres
    sollya_lib_clear_obj(initialSollyaPrecSo)
1839 85 storres
    if errorTypeIsNone:
1840 85 storres
        sollya_lib_clear_obj(errorTypeSo)
1841 162 storres
    sollya_lib_clear_obj(taylorFormSo)
1842 162 storres
    # Do not clear maxErrorSo.
1843 181 storres
    if absMaxErrorSa > absMinErrorSa:
1844 181 storres
        sollya_lib_clear_obj(absMinErrorSo)
1845 181 storres
        return((polyVarChangedSo, intervalCenterSo, absMaxErrorSo))
1846 181 storres
    else:
1847 181 storres
        sollya_lib_clear_obj(absMaxErrorSo)
1848 181 storres
        return((polyVarChangedSo, intervalCenterSo, absMinErrorSo))
1849 162 storres
# end pobyso_taylor_expansion_with_change_var_so_so
1850 56 storres
1851 5 storres
def pobyso_taylor(function, degree, point):
1852 38 storres
    """ Legacy function. See pobysoTaylor_so_so. """
1853 38 storres
    return(pobyso_taylor_so_so(function, degree, point))
1854 38 storres
1855 56 storres
def pobyso_taylor_so_so(functionSo, degreeSo, pointSo):
1856 56 storres
    return(sollya_lib_taylor(functionSo, degreeSo, pointSo))
1857 5 storres
1858 85 storres
def pobyso_taylorform(function, degree, point = None,
1859 85 storres
                      interval = None, errorType=None):
1860 85 storres
    """ Legacy function. See pobyso_taylorform_sa_sa;"""
1861 38 storres
1862 38 storres
def pobyso_taylorform_sa_sa(functionSa, \
1863 84 storres
                            degreeSa, \
1864 84 storres
                            pointSa, \
1865 84 storres
                            intervalSa=None, \
1866 84 storres
                            errorTypeSa=None, \
1867 84 storres
                            precisionSa=None):
1868 37 storres
    """
1869 85 storres
    Compute the Taylor form of 'degreeSa' for 'functionSa' at 'pointSa'
1870 85 storres
    for 'intervalSa' with 'errorTypeSa' (a string) using 'precisionSa'.
1871 37 storres
    point: must be a Real or a Real interval.
1872 37 storres
    return the Taylor form as an array
1873 83 storres
    TODO: take care of the interval and of the point when it is an interval;
1874 38 storres
          when errorType is not None;
1875 83 storres
          take care of the other elements of the Taylor form (coefficients
1876 83 storres
          errors and delta.
1877 37 storres
    """
1878 37 storres
    # Absolute as the default error.
1879 84 storres
    if errorTypeSa is None:
1880 37 storres
        errorTypeSo = sollya_lib_absolute()
1881 84 storres
    elif errorTypeSa == "relative":
1882 84 storres
        errorTypeSo = sollya_lib_relative()
1883 84 storres
    elif errortypeSa == "absolute":
1884 84 storres
        errorTypeSo = sollya_lib_absolute()
1885 37 storres
    else:
1886 84 storres
        # No clean up needed.
1887 84 storres
        return None
1888 84 storres
    # Global precision stuff
1889 226 storres
    sollyaPrecisionChangedSa = False
1890 226 storres
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
1891 226 storres
    if precisionSa is None:
1892 226 storres
        precSa = initialSollyaPrecSa
1893 226 storres
    else:
1894 226 storres
        if precSa > initialSollyaPrecSa:
1895 226 storres
            if precSa <= 2:
1896 226 storres
                print inspect.stack()[0][3], ":precision change <= 2 requested."
1897 226 storres
            pobyso_set_prec_sa_so(precSa)
1898 226 storres
            sollyaPrecisionChangedSa = True
1899 226 storres
    #
1900 85 storres
    if len(functionSa.variables()) > 0:
1901 85 storres
        varSa = functionSa.variables()[0]
1902 85 storres
        pobyso_name_free_variable_sa_so(str(varSa))
1903 84 storres
    # In any case (point or interval) the parent of pointSa has a precision
1904 84 storres
    # method.
1905 84 storres
    pointPrecSa = pointSa.parent().precision()
1906 226 storres
    if precSa > pointPrecSa:
1907 226 storres
        pointPrecSa = precSa
1908 84 storres
    # In any case (point or interval) pointSa has a base_ring() method.
1909 84 storres
    pointBaseRingString = str(pointSa.base_ring())
1910 84 storres
    if re.search('Interval', pointBaseRingString) is None: # Point
1911 84 storres
        pointSo = pobyso_constant_sa_so(pointSa, pointPrecSa)
1912 84 storres
    else: # Interval.
1913 84 storres
        pointSo = pobyso_interval_to_range_sa_so(pointSa, pointPrecSa)
1914 37 storres
    # Sollyafy the function.
1915 159 storres
    functionSo = pobyso_parse_string_sa_so(functionSa._assume_str().replace('_SAGE_VAR_', ''))
1916 37 storres
    if sollya_lib_obj_is_error(functionSo):
1917 37 storres
        print "pobyso_tailorform: function string can't be parsed!"
1918 37 storres
        return None
1919 37 storres
    # Sollyafy the degree
1920 84 storres
    degreeSo = sollya_lib_constant_from_int(int(degreeSa))
1921 37 storres
    # Sollyafy the point
1922 37 storres
    # Call Sollya
1923 83 storres
    taylorFormSo = \
1924 83 storres
        sollya_lib_taylorform(functionSo, degreeSo, pointSo, errorTypeSo,\
1925 37 storres
                                         None)
1926 85 storres
    sollya_lib_clear_obj(functionSo)
1927 85 storres
    sollya_lib_clear_obj(degreeSo)
1928 85 storres
    sollya_lib_clear_obj(pointSo)
1929 85 storres
    sollya_lib_clear_obj(errorTypeSo)
1930 38 storres
    (tfsAsList, numElements, isEndElliptic) = \
1931 38 storres
            pobyso_get_list_elements_so_so(taylorFormSo)
1932 37 storres
    polySo = tfsAsList[0]
1933 38 storres
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
1934 37 storres
    polyRealField = RealField(maxPrecision)
1935 38 storres
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, polyRealField)
1936 226 storres
    if sollyaPrecisionChangedSa:
1937 226 storres
        sollya_lib_set_prec(initialSollyaPrecSo)
1938 226 storres
    sollya_lib_clear_obj(initialSollyaPrecSo)
1939 37 storres
    polynomialRing = polyRealField[str(varSa)]
1940 37 storres
    polySa = polynomial(expSa, polynomialRing)
1941 37 storres
    taylorFormSa = [polySa]
1942 85 storres
    # Final clean-up.
1943 85 storres
    sollya_lib_clear_obj(taylorFormSo)
1944 51 storres
    return(taylorFormSa)
1945 51 storres
# End pobyso_taylor_form_sa_sa
1946 54 storres
1947 54 storres
def pobyso_taylorform_so_so(functionSo, degreeSo, pointSo, intervalSo=None, \
1948 54 storres
                            errorTypeSo=None):
1949 54 storres
    createdErrorType = False
1950 51 storres
    if errorTypeSo is None:
1951 51 storres
        errorTypeSo = sollya_lib_absolute()
1952 54 storres
        createdErrorType = True
1953 51 storres
    else:
1954 51 storres
        #TODO: deal with the other case.
1955 51 storres
        pass
1956 51 storres
    if intervalSo is None:
1957 54 storres
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
1958 54 storres
                                         errorTypeSo, None)
1959 51 storres
    else:
1960 54 storres
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
1961 54 storres
                                         intervalSo, errorTypeSo, None)
1962 54 storres
    if createdErrorType:
1963 54 storres
        sollya_lib_clear_obj(errorTypeSo)
1964 215 storres
    return resultSo
1965 51 storres
1966 37 storres
1967 37 storres
def pobyso_univar_polynomial_print_reverse(polySa):
1968 51 storres
    """ Legacy function. See pobyso_univar_polynomial_print_reverse_sa_sa. """
1969 51 storres
    return(pobyso_univar_polynomial_print_reverse_sa_sa(polySa))
1970 38 storres
1971 51 storres
def pobyso_univar_polynomial_print_reverse_sa_sa(polySa):
1972 37 storres
    """
1973 37 storres
    Return the string representation of a univariate polynomial with
1974 38 storres
    monomials ordered in the x^0..x^n order of the monomials.
1975 37 storres
    Remember: Sage
1976 37 storres
    """
1977 37 storres
    polynomialRing = polySa.base_ring()
1978 37 storres
    # A very expensive solution:
1979 37 storres
    # -create a fake multivariate polynomial field with only one variable,
1980 37 storres
    #   specifying a negative lexicographical order;
1981 37 storres
    mpolynomialRing = PolynomialRing(polynomialRing.base(), \
1982 37 storres
                                     polynomialRing.variable_name(), \
1983 37 storres
                                     1, order='neglex')
1984 37 storres
    # - convert the univariate argument polynomial into a multivariate
1985 37 storres
    #   version;
1986 37 storres
    p = mpolynomialRing(polySa)
1987 37 storres
    # - return the string representation of the converted form.
1988 37 storres
    # There is no simple str() method defined for p's class.
1989 37 storres
    return(p.__str__())
1990 5 storres
#
1991 230 storres
#print pobyso_get_prec()
1992 5 storres
pobyso_set_prec(165)
1993 230 storres
#print pobyso_get_prec()
1994 230 storres
#a=100
1995 230 storres
#print type(a)
1996 230 storres
#id(a)
1997 230 storres
#print "Max arity: ", pobyso_max_arity
1998 230 storres
#print "Function tripleDouble (43) as a string: ", pobyso_function_type_as_string(43)
1999 230 storres
#print "Function None (44) as a string: ", pobyso_function_type_as_string(44)
2000 56 storres
print "...Pobyso check done"