Statistiques
| Révision :

root / pobysoPythonSage / src / pobyso.py @ 258

Historique | Voir | Annoter | Télécharger (86,16 ko)

1 5 storres
"""
2 209 storres
@file pobyso.py
3 5 storres
Actual functions to use in Sage
4 256 storres
@author S.T.
5 256 storres
@date 2012-11-13
6 5 storres

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

10 38 storres
pobyso functions come in five flavors:
11 57 storres
- the _so_so (arguments and returned objects are pointers to Sollya objects,
12 57 storres
  includes the void function and the no arguments function that return a
13 57 storres
  pointer to a Sollya object);
14 38 storres
- the _so_sa (argument are pointers to Sollya objects, returned objects are
15 57 storres
  Sage/Python objects or, more generally, information is transfered from the
16 57 storres
  Sollya world to Sage/Python world; e.g. functions without arguments that
17 57 storres
  return a Sage/Python object);
18 38 storres
- the _sa_so (arguments are Sage/Python objects, returned objects are
19 38 storres
  pointers to Sollya objects);
20 38 storres
- the sa_sa (arguments and returned objects are all Sage/Python objects);
21 257 storres
- a catch all flavor, without any suffix, (e. g. functions that have no
22 257 storres
  argument nor return value).
23 257 storres

24 57 storres
This classification is not always very strict. Conversion functions from Sollya
25 57 storres
to Sage/Python are sometimes decorated with Sage/Python arguments to set
26 57 storres
the precision. These functions remain in the so_sa category.
27 258 storres

28 258 storres
@note
29 5 storres
Reported errors in Eclipse come from the calls to
30 5 storres
the Sollya library
31 5 storres

32 10 storres
ToDo (among other things):
33 10 storres
 -memory management.
34 5 storres
"""
35 5 storres
from ctypes import *
36 37 storres
import re
37 37 storres
from sage.symbolic.expression_conversions import polynomial
38 59 storres
from sage.symbolic.expression_conversions import PolynomialConverter
39 38 storres
"""
40 38 storres
Create the equivalent to an enum for the Sollya function types.
41 38 storres
"""
42 5 storres
(SOLLYA_BASE_FUNC_ABS,
43 5 storres
SOLLYA_BASE_FUNC_ACOS,
44 5 storres
    SOLLYA_BASE_FUNC_ACOSH,
45 5 storres
    SOLLYA_BASE_FUNC_ADD,
46 5 storres
    SOLLYA_BASE_FUNC_ASIN,
47 5 storres
    SOLLYA_BASE_FUNC_ASINH,
48 5 storres
    SOLLYA_BASE_FUNC_ATAN,
49 5 storres
    SOLLYA_BASE_FUNC_ATANH,
50 5 storres
    SOLLYA_BASE_FUNC_CEIL,
51 5 storres
    SOLLYA_BASE_FUNC_CONSTANT,
52 5 storres
    SOLLYA_BASE_FUNC_COS,
53 5 storres
    SOLLYA_BASE_FUNC_COSH,
54 5 storres
    SOLLYA_BASE_FUNC_DIV,
55 5 storres
    SOLLYA_BASE_FUNC_DOUBLE,
56 5 storres
    SOLLYA_BASE_FUNC_DOUBLEDOUBLE,
57 5 storres
    SOLLYA_BASE_FUNC_DOUBLEEXTENDED,
58 5 storres
    SOLLYA_BASE_FUNC_ERF,
59 5 storres
    SOLLYA_BASE_FUNC_ERFC,
60 5 storres
    SOLLYA_BASE_FUNC_EXP,
61 5 storres
    SOLLYA_BASE_FUNC_EXP_M1,
62 5 storres
    SOLLYA_BASE_FUNC_FLOOR,
63 5 storres
    SOLLYA_BASE_FUNC_FREE_VARIABLE,
64 5 storres
    SOLLYA_BASE_FUNC_HALFPRECISION,
65 5 storres
    SOLLYA_BASE_FUNC_LIBRARYCONSTANT,
66 5 storres
    SOLLYA_BASE_FUNC_LIBRARYFUNCTION,
67 5 storres
    SOLLYA_BASE_FUNC_LOG,
68 5 storres
    SOLLYA_BASE_FUNC_LOG_10,
69 5 storres
    SOLLYA_BASE_FUNC_LOG_1P,
70 5 storres
    SOLLYA_BASE_FUNC_LOG_2,
71 5 storres
    SOLLYA_BASE_FUNC_MUL,
72 5 storres
    SOLLYA_BASE_FUNC_NEARESTINT,
73 5 storres
    SOLLYA_BASE_FUNC_NEG,
74 5 storres
    SOLLYA_BASE_FUNC_PI,
75 5 storres
    SOLLYA_BASE_FUNC_POW,
76 5 storres
    SOLLYA_BASE_FUNC_PROCEDUREFUNCTION,
77 5 storres
    SOLLYA_BASE_FUNC_QUAD,
78 5 storres
    SOLLYA_BASE_FUNC_SIN,
79 5 storres
    SOLLYA_BASE_FUNC_SINGLE,
80 5 storres
    SOLLYA_BASE_FUNC_SINH,
81 5 storres
    SOLLYA_BASE_FUNC_SQRT,
82 5 storres
    SOLLYA_BASE_FUNC_SUB,
83 5 storres
    SOLLYA_BASE_FUNC_TAN,
84 5 storres
    SOLLYA_BASE_FUNC_TANH,
85 5 storres
SOLLYA_BASE_FUNC_TRIPLEDOUBLE) = map(int,xrange(44))
86 246 storres
sys.stderr.write("Superficial pobyso check...\n")
87 230 storres
#print "First constant - SOLLYA_BASE_FUNC_ABS: ", SOLLYA_BASE_FUNC_ABS
88 230 storres
#print "Last constant  - SOLLYA_BASE_FUNC_TRIPLEDOUBLE: ", SOLLYA_BASE_FUNC_TRIPLEDOUBLE
89 5 storres
90 5 storres
pobyso_max_arity = 9
91 5 storres
92 85 storres
def pobyso_absolute_so_so():
93 85 storres
    return(sollya_lib_absolute(None))
94 85 storres
95 5 storres
def pobyso_autoprint(arg):
96 218 storres
    sollya_lib_autoprint(arg, None)
97 5 storres
98 38 storres
def pobyso_autoprint_so_so(arg):
99 38 storres
    sollya_lib_autoprint(arg,None)
100 54 storres
101 84 storres
def pobyso_bounds_to_range_sa_so(rnLowerBoundSa, rnUpperBoundSa, \
102 84 storres
                                 precisionSa=None):
103 84 storres
    """
104 84 storres
    Return a Sollya range from to 2 RealField Sage elements.
105 84 storres
    The Sollya range element has a sufficient precision to hold all
106 162 storres
    the digits of the widest of the Sage bounds.
107 84 storres
    """
108 84 storres
    # Sanity check.
109 84 storres
    if rnLowerBoundSa > rnUpperBoundSa:
110 84 storres
        return None
111 115 storres
    # Precision stuff.
112 85 storres
    if precisionSa is None:
113 84 storres
        # Check for the largest precision.
114 84 storres
        lbPrecSa = rnLowerBoundSa.parent().precision()
115 84 storres
        ubPrecSa = rnLowerBoundSa.parent().precision()
116 84 storres
        maxPrecSa = max(lbPrecSa, ubPrecSa)
117 84 storres
    else:
118 84 storres
        maxPrecSa = precisionSa
119 115 storres
    # From Sage to Sollya bounds.
120 162 storres
#    lowerBoundSo = sollya_lib_constant(get_rn_value(rnLowerBoundSa),
121 162 storres
#                                       maxPrecSa)
122 162 storres
    lowerBoundSo = pobyso_constant_sa_so(rnLowerBoundSa,
123 162 storres
                                         maxPrecSa)
124 162 storres
    upperBoundSo = pobyso_constant_sa_so(rnUpperBoundSa,
125 154 storres
                                       maxPrecSa)
126 162 storres
127 115 storres
    # From Sollya bounds to range.
128 84 storres
    rangeSo = sollya_lib_range(lowerBoundSo, upperBoundSo)
129 84 storres
    # Back to original precision.
130 84 storres
    # Clean up
131 84 storres
    sollya_lib_clear_obj(lowerBoundSo)
132 84 storres
    sollya_lib_clear_obj(upperBoundSo)
133 154 storres
    return rangeSo
134 84 storres
# End pobyso_bounds_to_range_sa_so
135 84 storres
136 215 storres
def pobyso_build_end_elliptic_list_so_so(*args):
137 215 storres
    """
138 237 storres
    From Sollya argument objects, create a Sollya end elliptic list.
139 237 storres
    Elements of the list are "eaten" (should not be cleared individually,
140 215 storres
    are cleared when the list is cleared).
141 215 storres
    """
142 215 storres
    if len(args) == 0:
143 215 storres
        ## Called with an empty list produced "error".
144 215 storres
        return sollya_lib_build_end_elliptic_list(None)
145 215 storres
    index = 0
146 215 storres
    ## One can not append elements to an elliptic list, prepend only is
147 215 storres
    #  permitted.
148 215 storres
    for argument in reversed(args):
149 215 storres
        if index == 0:
150 215 storres
            listSo = sollya_lib_build_end_elliptic_list(argument, None)
151 215 storres
        else:
152 215 storres
            listSo = sollya_lib_prepend(argument, listSo)
153 215 storres
        index += 1
154 215 storres
    return listSo
155 215 storres
156 215 storres
# End pobyso_build_end_elliptic_list_so_so
157 215 storres
158 54 storres
def pobyso_build_function_sub_so_so(exp1So, exp2So):
159 237 storres
    return sollya_lib_build_function_sub(exp1So, exp2So)
160 54 storres
161 237 storres
def pobyso_build_list_of_ints_sa_so(*args):
162 237 storres
    """
163 237 storres
    Build a Sollya list from Sage integral arguments.
164 237 storres
    """
165 237 storres
    if len(args) == 0:
166 237 storres
        return pobyso_build_list_so_so()
167 237 storres
    ## Make a Sage list of integral constants.
168 237 storres
    intsList = []
169 237 storres
    for intElem in args:
170 237 storres
        intsList.append(pobyso_constant_from_int_sa_so(intElem))
171 237 storres
    return pobyso_build_list_so_so(*intsList)
172 237 storres
173 237 storres
def pobyso_build_list_so_so(*args):
174 237 storres
    """
175 237 storres
    Make a Sollya list out of Sollya objects passed as arguments.
176 237 storres
    If one wants to call it with a list argument, should prepend a "*"
177 237 storres
    before the list variable name.
178 237 storres
    Elements of the list are "eaten" (should not be cleared individually,
179 237 storres
    are cleared when the list is cleared).
180 237 storres
    """
181 237 storres
    if len(args) == 0:
182 237 storres
        ## Called with an empty list produced "error".
183 237 storres
        return sollya_lib_build_list(None)
184 237 storres
    index = 0
185 237 storres
    ## Append the Sollya elements one by one.
186 237 storres
    for elementSo in args:
187 237 storres
        if index == 0:
188 237 storres
            listSo = sollya_lib_build_list(elementSo, None)
189 237 storres
        else:
190 237 storres
            listSo = sollya_lib_append(listSo, elementSo)
191 237 storres
        index += 1
192 237 storres
    return listSo
193 237 storres
# End pobyso_build list_so_so
194 237 storres
195 237 storres
196 85 storres
def pobyso_change_var_in_function_so_so(funcSo, chvarExpSo):
197 54 storres
    """
198 85 storres
    Variable change in a function.
199 85 storres
    """
200 85 storres
    return(sollya_lib_evaluate(funcSo,chvarExpSo))
201 85 storres
# End pobyso_change_var_in_function_so_so
202 85 storres
203 85 storres
def pobyso_chebyshevform_so_so(functionSo, degreeSo, intervalSo):
204 85 storres
    resultSo = sollya_lib_chebyshevform(functionSo, degreeSo, intervalSo)
205 85 storres
    return(resultSo)
206 85 storres
# End pobyso_chebyshevform_so_so.
207 85 storres
208 237 storres
def pobyso_clear_obj(objSo):
209 237 storres
    """
210 237 storres
    Free a Sollya object's memory.
211 237 storres
    Very thin wrapper around sollya_lib_clear_obj().
212 237 storres
    """
213 237 storres
    sollya_lib_clear_obj(objSo)
214 237 storres
# End pobyso_clear_obj.
215 237 storres
216 117 storres
def pobyso_clear_taylorform_sa_so(taylorFormSaSo):
217 117 storres
    """
218 117 storres
    This method is necessary to correctly clean up the memory from Taylor forms.
219 117 storres
    These are made of a Sollya object, a Sollya object list, a Sollya object.
220 117 storres
    For no clearly understood reason, sollya_lib_clear_object_list crashed
221 117 storres
    when applied to the object list.
222 117 storres
    Here, we decompose it into Sage list of Sollya objects references and we
223 237 storres
     clear them one at a time.
224 117 storres
    """
225 117 storres
    sollya_lib_clear_obj(taylorFormSaSo[0])
226 117 storres
    (coefficientsErrorsListSaSo, numElementsSa, isEndEllipticSa) = \
227 117 storres
        pobyso_get_list_elements_so_so(taylorFormSaSo[1])
228 117 storres
    for element in coefficientsErrorsListSaSo:
229 117 storres
        sollya_lib_clear_obj(element)
230 117 storres
    sollya_lib_clear_obj(taylorFormSaSo[1])
231 117 storres
    sollya_lib_clear_obj(taylorFormSaSo[2])
232 117 storres
# End pobyso_clear_taylorform_sa_so
233 117 storres
234 85 storres
def pobyso_cmp(rnArgSa, cteSo):
235 85 storres
    """
236 237 storres
    Deprecated, use pobyso_cmp_sa_so_sa instead.
237 237 storres
    """
238 237 storres
    print "Deprecated, use pobyso_cmp_sa_so_sa instead."
239 237 storres
    return pobyso_cmp_sa_so_sa(rnArgSa, cteSo)
240 237 storres
# End  pobyso_cmp
241 237 storres
242 237 storres
def pobyso_cmp_sa_so_sa(rnArgSa, cteSo):
243 237 storres
    """
244 54 storres
    Compare the MPFR value a RealNumber with that of a Sollya constant.
245 54 storres

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

802 118 storres
    INPUT:
803 118 storres
    - objectListSo: a Sollya list of Sollya objects.
804 118 storres

805 118 storres
    OUTPUT:
806 118 storres
    - a Sage/Python tuple made of:
807 118 storres
      - a Sage/Python list of Sollya objects,
808 118 storres
      - a Sage/Python int holding the number of elements,
809 118 storres
      - a Sage/Python int stating (!= 0) that the list is end-elliptic.
810 118 storres
    NOTE::
811 118 storres
        We recover the addresses of the Sollya object from the list of pointers
812 118 storres
        returned by sollya_lib_get_list_elements. The list itself is freed.
813 118 storres
    TODO::
814 118 storres
        Figure out what to do with numElements since the number of elements
815 118 storres
        can easily be recovered from the list itself.
816 118 storres
        Ditto for isEndElliptic.
817 51 storres
    """
818 5 storres
    listAddress = POINTER(c_longlong)()
819 5 storres
    numElements = c_int(0)
820 5 storres
    isEndElliptic = c_int(0)
821 117 storres
    listAsSageList = []
822 5 storres
    result = sollya_lib_get_list_elements(byref(listAddress),\
823 54 storres
                                          byref(numElements),\
824 54 storres
                                          byref(isEndElliptic),\
825 117 storres
                                          objectListSo)
826 5 storres
    if result == 0 :
827 5 storres
        return None
828 5 storres
    for i in xrange(0, numElements.value, 1):
829 118 storres
       #listAsSageList.append(sollya_lib_copy_obj(listAddress[i]))
830 118 storres
       listAsSageList.append(listAddress[i])
831 117 storres
       # Clear each of the elements returned by Sollya.
832 118 storres
       #sollya_lib_clear_obj(listAddress[i])
833 117 storres
    # Free the list itself.
834 117 storres
    sollya_lib_free(listAddress)
835 209 storres
    return (listAsSageList, numElements.value, isEndElliptic.value)
836 5 storres
837 38 storres
def pobyso_get_max_prec_of_exp(soExp):
838 38 storres
    """ Legacy function. See pobyso_get_max_prec_of_exp_so_sa. """
839 209 storres
    return pobyso_get_max_prec_of_exp_so_sa(soExp)
840 5 storres
841 85 storres
def pobyso_get_max_prec_of_exp_so_sa(expSo):
842 38 storres
    """
843 38 storres
    Get the maximum precision used for the numbers in a Sollya expression.
844 52 storres

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

1905 242 storres
    The returned value is the maximum of the absolute values of the range
1906 242 storres
    elements  returned by the Sollya supnorm functions
1907 242 storres
    """
1908 242 storres
    supNormRangeSo = pobyso_supnorm_so_so(polySo,
1909 242 storres
                                          funcSo,
1910 242 storres
                                          intervalSo,
1911 242 storres
                                          errorTypeSo,
1912 242 storres
                                          accuracySo)
1913 242 storres
    supNormSo = pobyso_range_max_abs_so_so(supNormRangeSo)
1914 242 storres
    sollya_lib_clear_obj(supNormRangeSo)
1915 242 storres
    #pobyso_autoprint(supNormSo)
1916 242 storres
    supNormSa = pobyso_get_constant_as_rn_with_rf_so_sa(supNormSo, realFieldSa)
1917 242 storres
    sollya_lib_clear_obj(supNormSo)
1918 242 storres
    return supNormSa
1919 242 storres
# End pobyso_supnorm_so_sa.
1920 242 storres
#
1921 85 storres
def pobyso_supnorm_so_so(polySo, funcSo, intervalSo, errorTypeSo = None,\
1922 85 storres
                         accuracySo = None):
1923 58 storres
    """
1924 85 storres
    Computes the supnorm of the approximation error between the given
1925 242 storres
    polynomial and function. Attention: returns a range!
1926 85 storres
    errorTypeSo defaults to "absolute".
1927 85 storres
    accuracySo defaults to 2^(-40).
1928 85 storres
    """
1929 85 storres
    if errorTypeSo is None:
1930 85 storres
        errorTypeSo = sollya_lib_absolute(None)
1931 85 storres
        errorTypeIsNone = True
1932 85 storres
    else:
1933 85 storres
        errorTypeIsNone = False
1934 85 storres
    #
1935 85 storres
    if accuracySo is None:
1936 237 storres
        # Notice the **: we are in Pythonland here!
1937 85 storres
        accuracySo = pobyso_constant_sa_so(RR(2**(-40)))
1938 85 storres
        accuracyIsNone = True
1939 85 storres
    else:
1940 85 storres
        accuracyIsNone = False
1941 238 storres
    #pobyso_autoprint(accuracySo)
1942 85 storres
    resultSo = \
1943 85 storres
        sollya_lib_supnorm(polySo, funcSo, intervalSo, errorTypeSo, \
1944 85 storres
                              accuracySo)
1945 85 storres
    if errorTypeIsNone:
1946 85 storres
        sollya_lib_clear_obj(errorTypeSo)
1947 85 storres
    if accuracyIsNone:
1948 85 storres
        sollya_lib_clear_obj(accuracySo)
1949 85 storres
    return resultSo
1950 85 storres
# End pobyso_supnorm_so_so
1951 242 storres
#
1952 162 storres
def pobyso_taylor_expansion_no_change_var_so_so(functionSo,
1953 162 storres
                                                degreeSo,
1954 162 storres
                                                rangeSo,
1955 162 storres
                                                errorTypeSo=None,
1956 162 storres
                                                sollyaPrecSo=None):
1957 85 storres
    """
1958 162 storres
    Compute the Taylor expansion without the variable change
1959 162 storres
    x -> x-intervalCenter.
1960 58 storres
    """
1961 226 storres
    # Change internal Sollya precision, if needed.
1962 227 storres
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
1963 226 storres
    sollyaPrecChanged = False
1964 226 storres
    if sollyaPrecSo is None:
1965 226 storres
        pass
1966 226 storres
    else:
1967 58 storres
        sollya_lib_set_prec(sollyaPrecSo)
1968 226 storres
        sollyaPrecChanged = True
1969 85 storres
    # Error type stuff: default to absolute.
1970 85 storres
    if errorTypeSo is None:
1971 85 storres
        errorTypeIsNone = True
1972 85 storres
        errorTypeSo = sollya_lib_absolute(None)
1973 85 storres
    else:
1974 85 storres
        errorTypeIsNone = False
1975 162 storres
    intervalCenterSo = sollya_lib_mid(rangeSo, None)
1976 162 storres
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo,
1977 162 storres
                                         intervalCenterSo,
1978 58 storres
                                         rangeSo, errorTypeSo, None)
1979 117 storres
    # taylorFormListSaSo is a Python list of Sollya objects references that
1980 117 storres
    # are copies of the elements of taylorFormSo.
1981 117 storres
    # pobyso_get_list_elements_so_so clears taylorFormSo.
1982 162 storres
    (taylorFormListSaSo, numElementsSa, isEndEllipticSa) = \
1983 58 storres
        pobyso_get_list_elements_so_so(taylorFormSo)
1984 162 storres
    polySo = sollya_lib_copy_obj(taylorFormListSaSo[0])
1985 162 storres
    #print "Num elements:", numElementsSa
1986 162 storres
    sollya_lib_clear_obj(taylorFormSo)
1987 162 storres
    #polySo = taylorFormListSaSo[0]
1988 162 storres
    #errorRangeSo = sollya_lib_copy_obj(taylorFormListSaSo[2])
1989 162 storres
    errorRangeSo = taylorFormListSaSo[2]
1990 181 storres
    # No copy_obj needed here: a new objects are created.
1991 181 storres
    maxErrorSo    = sollya_lib_sup(errorRangeSo)
1992 181 storres
    minErrorSo    = sollya_lib_inf(errorRangeSo)
1993 181 storres
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
1994 181 storres
    absMinErrorSo = sollya_lib_abs(minErrorSo)
1995 181 storres
    sollya_lib_clear_obj(maxErrorSo)
1996 181 storres
    sollya_lib_clear_obj(minErrorSo)
1997 181 storres
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
1998 181 storres
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
1999 58 storres
    # If changed, reset the Sollya working precision.
2000 226 storres
    if sollyaPrecChanged:
2001 226 storres
        sollya_lib_set_prec(initialSollyaPrecSo)
2002 226 storres
    sollya_lib_clear_obj(initialSollyaPrecSo)
2003 85 storres
    if errorTypeIsNone:
2004 85 storres
        sollya_lib_clear_obj(errorTypeSo)
2005 162 storres
    pobyso_clear_taylorform_sa_so(taylorFormListSaSo)
2006 181 storres
    if absMaxErrorSa > absMinErrorSa:
2007 181 storres
        sollya_lib_clear_obj(absMinErrorSo)
2008 227 storres
        return (polySo, intervalCenterSo, absMaxErrorSo)
2009 181 storres
    else:
2010 181 storres
        sollya_lib_clear_obj(absMaxErrorSo)
2011 227 storres
        return (polySo, intervalCenterSo, absMinErrorSo)
2012 162 storres
# end pobyso_taylor_expansion_no_change_var_so_so
2013 58 storres
2014 162 storres
def pobyso_taylor_expansion_with_change_var_so_so(functionSo, degreeSo, \
2015 162 storres
                                                  rangeSo, \
2016 162 storres
                                                  errorTypeSo=None, \
2017 162 storres
                                                  sollyaPrecSo=None):
2018 58 storres
    """
2019 162 storres
    Compute the Taylor expansion with the variable change
2020 162 storres
    x -> (x-intervalCenter) included.
2021 58 storres
    """
2022 226 storres
    # Change Sollya internal precision, if need.
2023 226 storres
    sollyaPrecChanged = False
2024 226 storres
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_sos_sa()
2025 226 storres
    if sollyaPrecSo is None:
2026 226 storres
        pass
2027 226 storres
    else:
2028 56 storres
        sollya_lib_set_prec(sollyaPrecSo)
2029 226 storres
        sollyaPrecChanged = True
2030 162 storres
    #
2031 85 storres
    # Error type stuff: default to absolute.
2032 85 storres
    if errorTypeSo is None:
2033 85 storres
        errorTypeIsNone = True
2034 85 storres
        errorTypeSo = sollya_lib_absolute(None)
2035 85 storres
    else:
2036 85 storres
        errorTypeIsNone = False
2037 162 storres
    intervalCenterSo = sollya_lib_mid(rangeSo)
2038 162 storres
    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo, \
2039 162 storres
                                         intervalCenterSo, \
2040 56 storres
                                         rangeSo, errorTypeSo, None)
2041 116 storres
    # taylorFormListSaSo is a Python list of Sollya objects references that
2042 116 storres
    # are copies of the elements of taylorFormSo.
2043 116 storres
    # pobyso_get_list_elements_so_so clears taylorFormSo.
2044 162 storres
    (taylorFormListSo, numElements, isEndElliptic) = \
2045 56 storres
        pobyso_get_list_elements_so_so(taylorFormSo)
2046 162 storres
    polySo = taylorFormListSo[0]
2047 162 storres
    errorRangeSo = taylorFormListSo[2]
2048 181 storres
    maxErrorSo    = sollya_lib_sup(errorRangeSo)
2049 181 storres
    minErrorSo    = sollya_lib_inf(errorRangeSo)
2050 181 storres
    absMaxErrorSo = sollya_lib_abs(maxErrorSo)
2051 181 storres
    absMinErrorSo = sollya_lib_abs(minErrorSo)
2052 181 storres
    sollya_lib_clear_obj(maxErrorSo)
2053 181 storres
    sollya_lib_clear_obj(minErrorSo)
2054 181 storres
    absMaxErrorSa = pobyso_get_constant_as_rn_so_sa(absMaxErrorSo)
2055 181 storres
    absMinErrorSa = pobyso_get_constant_as_rn_so_sa(absMinErrorSo)
2056 162 storres
    changeVarExpSo = sollya_lib_build_function_sub(\
2057 162 storres
                       sollya_lib_build_function_free_variable(),\
2058 162 storres
                       sollya_lib_copy_obj(intervalCenterSo))
2059 181 storres
    polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpSo)
2060 181 storres
    sollya_lib_clear_obj(polySo)
2061 162 storres
    sollya_lib_clear_obj(changeVarExpSo)
2062 56 storres
    # If changed, reset the Sollya working precision.
2063 226 storres
    if sollyaPrecChanged:
2064 226 storres
        sollya_lib_set_prec(initialSollyaPrecSo)
2065 226 storres
    sollya_lib_clear_obj(initialSollyaPrecSo)
2066 85 storres
    if errorTypeIsNone:
2067 85 storres
        sollya_lib_clear_obj(errorTypeSo)
2068 162 storres
    sollya_lib_clear_obj(taylorFormSo)
2069 162 storres
    # Do not clear maxErrorSo.
2070 181 storres
    if absMaxErrorSa > absMinErrorSa:
2071 181 storres
        sollya_lib_clear_obj(absMinErrorSo)
2072 181 storres
        return((polyVarChangedSo, intervalCenterSo, absMaxErrorSo))
2073 181 storres
    else:
2074 181 storres
        sollya_lib_clear_obj(absMaxErrorSo)
2075 181 storres
        return((polyVarChangedSo, intervalCenterSo, absMinErrorSo))
2076 162 storres
# end pobyso_taylor_expansion_with_change_var_so_so
2077 56 storres
2078 5 storres
def pobyso_taylor(function, degree, point):
2079 38 storres
    """ Legacy function. See pobysoTaylor_so_so. """
2080 38 storres
    return(pobyso_taylor_so_so(function, degree, point))
2081 38 storres
2082 56 storres
def pobyso_taylor_so_so(functionSo, degreeSo, pointSo):
2083 56 storres
    return(sollya_lib_taylor(functionSo, degreeSo, pointSo))
2084 5 storres
2085 85 storres
def pobyso_taylorform(function, degree, point = None,
2086 85 storres
                      interval = None, errorType=None):
2087 85 storres
    """ Legacy function. See pobyso_taylorform_sa_sa;"""
2088 38 storres
2089 38 storres
def pobyso_taylorform_sa_sa(functionSa, \
2090 84 storres
                            degreeSa, \
2091 84 storres
                            pointSa, \
2092 84 storres
                            intervalSa=None, \
2093 84 storres
                            errorTypeSa=None, \
2094 84 storres
                            precisionSa=None):
2095 37 storres
    """
2096 85 storres
    Compute the Taylor form of 'degreeSa' for 'functionSa' at 'pointSa'
2097 85 storres
    for 'intervalSa' with 'errorTypeSa' (a string) using 'precisionSa'.
2098 37 storres
    point: must be a Real or a Real interval.
2099 37 storres
    return the Taylor form as an array
2100 83 storres
    TODO: take care of the interval and of the point when it is an interval;
2101 38 storres
          when errorType is not None;
2102 83 storres
          take care of the other elements of the Taylor form (coefficients
2103 83 storres
          errors and delta.
2104 37 storres
    """
2105 37 storres
    # Absolute as the default error.
2106 84 storres
    if errorTypeSa is None:
2107 37 storres
        errorTypeSo = sollya_lib_absolute()
2108 84 storres
    elif errorTypeSa == "relative":
2109 84 storres
        errorTypeSo = sollya_lib_relative()
2110 84 storres
    elif errortypeSa == "absolute":
2111 84 storres
        errorTypeSo = sollya_lib_absolute()
2112 37 storres
    else:
2113 84 storres
        # No clean up needed.
2114 84 storres
        return None
2115 84 storres
    # Global precision stuff
2116 226 storres
    sollyaPrecisionChangedSa = False
2117 226 storres
    (initialSollyaPrecSo, initialSollyaPrecSa) = pobyso_get_prec_so_so_sa()
2118 226 storres
    if precisionSa is None:
2119 226 storres
        precSa = initialSollyaPrecSa
2120 226 storres
    else:
2121 226 storres
        if precSa > initialSollyaPrecSa:
2122 226 storres
            if precSa <= 2:
2123 226 storres
                print inspect.stack()[0][3], ":precision change <= 2 requested."
2124 226 storres
            pobyso_set_prec_sa_so(precSa)
2125 226 storres
            sollyaPrecisionChangedSa = True
2126 226 storres
    #
2127 85 storres
    if len(functionSa.variables()) > 0:
2128 85 storres
        varSa = functionSa.variables()[0]
2129 85 storres
        pobyso_name_free_variable_sa_so(str(varSa))
2130 84 storres
    # In any case (point or interval) the parent of pointSa has a precision
2131 84 storres
    # method.
2132 84 storres
    pointPrecSa = pointSa.parent().precision()
2133 226 storres
    if precSa > pointPrecSa:
2134 226 storres
        pointPrecSa = precSa
2135 84 storres
    # In any case (point or interval) pointSa has a base_ring() method.
2136 84 storres
    pointBaseRingString = str(pointSa.base_ring())
2137 84 storres
    if re.search('Interval', pointBaseRingString) is None: # Point
2138 84 storres
        pointSo = pobyso_constant_sa_so(pointSa, pointPrecSa)
2139 84 storres
    else: # Interval.
2140 84 storres
        pointSo = pobyso_interval_to_range_sa_so(pointSa, pointPrecSa)
2141 37 storres
    # Sollyafy the function.
2142 159 storres
    functionSo = pobyso_parse_string_sa_so(functionSa._assume_str().replace('_SAGE_VAR_', ''))
2143 37 storres
    if sollya_lib_obj_is_error(functionSo):
2144 37 storres
        print "pobyso_tailorform: function string can't be parsed!"
2145 37 storres
        return None
2146 37 storres
    # Sollyafy the degree
2147 84 storres
    degreeSo = sollya_lib_constant_from_int(int(degreeSa))
2148 37 storres
    # Sollyafy the point
2149 37 storres
    # Call Sollya
2150 83 storres
    taylorFormSo = \
2151 83 storres
        sollya_lib_taylorform(functionSo, degreeSo, pointSo, errorTypeSo,\
2152 37 storres
                                         None)
2153 85 storres
    sollya_lib_clear_obj(functionSo)
2154 85 storres
    sollya_lib_clear_obj(degreeSo)
2155 85 storres
    sollya_lib_clear_obj(pointSo)
2156 85 storres
    sollya_lib_clear_obj(errorTypeSo)
2157 38 storres
    (tfsAsList, numElements, isEndElliptic) = \
2158 38 storres
            pobyso_get_list_elements_so_so(taylorFormSo)
2159 37 storres
    polySo = tfsAsList[0]
2160 38 storres
    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)
2161 37 storres
    polyRealField = RealField(maxPrecision)
2162 38 storres
    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, polyRealField)
2163 226 storres
    if sollyaPrecisionChangedSa:
2164 226 storres
        sollya_lib_set_prec(initialSollyaPrecSo)
2165 226 storres
    sollya_lib_clear_obj(initialSollyaPrecSo)
2166 37 storres
    polynomialRing = polyRealField[str(varSa)]
2167 37 storres
    polySa = polynomial(expSa, polynomialRing)
2168 37 storres
    taylorFormSa = [polySa]
2169 85 storres
    # Final clean-up.
2170 85 storres
    sollya_lib_clear_obj(taylorFormSo)
2171 51 storres
    return(taylorFormSa)
2172 51 storres
# End pobyso_taylor_form_sa_sa
2173 54 storres
2174 54 storres
def pobyso_taylorform_so_so(functionSo, degreeSo, pointSo, intervalSo=None, \
2175 54 storres
                            errorTypeSo=None):
2176 54 storres
    createdErrorType = False
2177 51 storres
    if errorTypeSo is None:
2178 51 storres
        errorTypeSo = sollya_lib_absolute()
2179 54 storres
        createdErrorType = True
2180 51 storres
    else:
2181 51 storres
        #TODO: deal with the other case.
2182 51 storres
        pass
2183 51 storres
    if intervalSo is None:
2184 54 storres
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
2185 54 storres
                                         errorTypeSo, None)
2186 51 storres
    else:
2187 54 storres
        resultSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, \
2188 54 storres
                                         intervalSo, errorTypeSo, None)
2189 54 storres
    if createdErrorType:
2190 54 storres
        sollya_lib_clear_obj(errorTypeSo)
2191 215 storres
    return resultSo
2192 51 storres
2193 37 storres
2194 37 storres
def pobyso_univar_polynomial_print_reverse(polySa):
2195 51 storres
    """ Legacy function. See pobyso_univar_polynomial_print_reverse_sa_sa. """
2196 51 storres
    return(pobyso_univar_polynomial_print_reverse_sa_sa(polySa))
2197 38 storres
2198 51 storres
def pobyso_univar_polynomial_print_reverse_sa_sa(polySa):
2199 37 storres
    """
2200 37 storres
    Return the string representation of a univariate polynomial with
2201 38 storres
    monomials ordered in the x^0..x^n order of the monomials.
2202 37 storres
    Remember: Sage
2203 37 storres
    """
2204 37 storres
    polynomialRing = polySa.base_ring()
2205 37 storres
    # A very expensive solution:
2206 37 storres
    # -create a fake multivariate polynomial field with only one variable,
2207 37 storres
    #   specifying a negative lexicographical order;
2208 37 storres
    mpolynomialRing = PolynomialRing(polynomialRing.base(), \
2209 37 storres
                                     polynomialRing.variable_name(), \
2210 37 storres
                                     1, order='neglex')
2211 37 storres
    # - convert the univariate argument polynomial into a multivariate
2212 37 storres
    #   version;
2213 37 storres
    p = mpolynomialRing(polySa)
2214 37 storres
    # - return the string representation of the converted form.
2215 37 storres
    # There is no simple str() method defined for p's class.
2216 37 storres
    return(p.__str__())
2217 5 storres
#
2218 230 storres
#print pobyso_get_prec()
2219 5 storres
pobyso_set_prec(165)
2220 230 storres
#print pobyso_get_prec()
2221 230 storres
#a=100
2222 230 storres
#print type(a)
2223 230 storres
#id(a)
2224 230 storres
#print "Max arity: ", pobyso_max_arity
2225 230 storres
#print "Function tripleDouble (43) as a string: ", pobyso_function_type_as_string(43)
2226 230 storres
#print "Function None (44) as a string: ", pobyso_function_type_as_string(44)
2227 246 storres
sys.stderr.write("\t...Pobyso check done.\n")