Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 222

Historique | Voir | Annoter | Télécharger (69,4 ko)

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

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

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

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

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

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

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

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