Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 234

Historique | Voir | Annoter | Télécharger (79,08 ko)

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

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

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

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

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

717 118 storres
    INPUT:
718 118 storres
    - objectListSo: a Sollya list of Sollya objects.
719 118 storres

720 118 storres
    OUTPUT:
721 118 storres
    - a Sage/Python tuple made of:
722 118 storres
      - a Sage/Python list of Sollya objects,
723 118 storres
      - a Sage/Python int holding the number of elements,
724 118 storres
      - a Sage/Python int stating (!= 0) that the list is end-elliptic.
725 118 storres
    NOTE::
726 118 storres
        We recover the addresses of the Sollya object from the list of pointers
727 118 storres
        returned by sollya_lib_get_list_elements. The list itself is freed.
728 118 storres
    TODO::
729 118 storres
        Figure out what to do with numElements since the number of elements
730 118 storres
        can easily be recovered from the list itself.
731 118 storres
        Ditto for isEndElliptic.
732 51 storres
    """
733 5 storres
    listAddress = POINTER(c_longlong)()
734 5 storres
    numElements = c_int(0)
735 5 storres
    isEndElliptic = c_int(0)
736 117 storres
    listAsSageList = []
737 5 storres
    result = sollya_lib_get_list_elements(byref(listAddress),\
738 54 storres
                                          byref(numElements),\
739 54 storres
                                          byref(isEndElliptic),\
740 117 storres
                                          objectListSo)
741 5 storres
    if result == 0 :
742 5 storres
        return None
743 5 storres
    for i in xrange(0, numElements.value, 1):
744 118 storres
       #listAsSageList.append(sollya_lib_copy_obj(listAddress[i]))
745 118 storres
       listAsSageList.append(listAddress[i])
746 117 storres
       # Clear each of the elements returned by Sollya.
747 118 storres
       #sollya_lib_clear_obj(listAddress[i])
748 117 storres
    # Free the list itself.
749 117 storres
    sollya_lib_free(listAddress)
750 209 storres
    return (listAsSageList, numElements.value, isEndElliptic.value)
751 5 storres
752 38 storres
def pobyso_get_max_prec_of_exp(soExp):
753 38 storres
    """ Legacy function. See pobyso_get_max_prec_of_exp_so_sa. """
754 209 storres
    return pobyso_get_max_prec_of_exp_so_sa(soExp)
755 5 storres
756 85 storres
def pobyso_get_max_prec_of_exp_so_sa(expSo):
757 38 storres
    """
758 38 storres
    Get the maximum precision used for the numbers in a Sollya expression.
759 52 storres

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