Laboratoire de l'Informatique et du Parallélisme » Pobyso

"""

@file pobyso.py

Actual functions to use in Sage

ST 2012-11-13

Command line syntax:

  use from Sage (via the "load" or the "attach" commands)

pobyso functions come in five flavors: 

- the _so_so (arguments and returned objects are pointers to Sollya objects, 

  includes the void function and the no arguments function that return a 

  pointer to a Sollya object);

- the _so_sa (argument are pointers to Sollya objects, returned objects are

  Sage/Python objects or, more generally, information is transfered from the

  Sollya world to Sage/Python world; e.g. functions without arguments that

  return a Sage/Python object);

- the _sa_so (arguments are Sage/Python objects, returned objects are 

  pointers to Sollya objects);

- the sa_sa (arguments and returned objects are all Sage/Python objects);

- a catch all flavor, without any suffix, (e. g. functions that have no argument 

  nor return value).

This classification is not always very strict. Conversion functions from Sollya

to Sage/Python are sometimes decorated with Sage/Python arguments to set

the precision. These functions remain in the so_sa category.

NOTES:

Reported errors in Eclipse come from the calls to

the Sollya library

ToDo (among other things): 

 -memory management.

"""

from ctypes import *

import re

from sage.symbolic.expression_conversions import polynomial

from sage.symbolic.expression_conversions import PolynomialConverter

"""

Create the equivalent to an enum for the Sollya function types.

"""

(SOLLYA_BASE_FUNC_ABS,

SOLLYA_BASE_FUNC_ACOS,

    SOLLYA_BASE_FUNC_ACOSH,

    SOLLYA_BASE_FUNC_ADD,

    SOLLYA_BASE_FUNC_ASIN,

    SOLLYA_BASE_FUNC_ASINH,

    SOLLYA_BASE_FUNC_ATAN,

    SOLLYA_BASE_FUNC_ATANH,

    SOLLYA_BASE_FUNC_CEIL,

    SOLLYA_BASE_FUNC_CONSTANT,

    SOLLYA_BASE_FUNC_COS,

    SOLLYA_BASE_FUNC_COSH,

    SOLLYA_BASE_FUNC_DIV,

    SOLLYA_BASE_FUNC_DOUBLE,

    SOLLYA_BASE_FUNC_DOUBLEDOUBLE,

    SOLLYA_BASE_FUNC_DOUBLEEXTENDED,

    SOLLYA_BASE_FUNC_ERF,

    SOLLYA_BASE_FUNC_ERFC,

    SOLLYA_BASE_FUNC_EXP,

    SOLLYA_BASE_FUNC_EXP_M1,

    SOLLYA_BASE_FUNC_FLOOR,

    SOLLYA_BASE_FUNC_FREE_VARIABLE,

    SOLLYA_BASE_FUNC_HALFPRECISION,

    SOLLYA_BASE_FUNC_LIBRARYCONSTANT,

    SOLLYA_BASE_FUNC_LIBRARYFUNCTION,

    SOLLYA_BASE_FUNC_LOG,

    SOLLYA_BASE_FUNC_LOG_10,

    SOLLYA_BASE_FUNC_LOG_1P,

    SOLLYA_BASE_FUNC_LOG_2,

    SOLLYA_BASE_FUNC_MUL,

    SOLLYA_BASE_FUNC_NEARESTINT,

    SOLLYA_BASE_FUNC_NEG,

    SOLLYA_BASE_FUNC_PI,

    SOLLYA_BASE_FUNC_POW,

    SOLLYA_BASE_FUNC_PROCEDUREFUNCTION,

    SOLLYA_BASE_FUNC_QUAD,

    SOLLYA_BASE_FUNC_SIN,

    SOLLYA_BASE_FUNC_SINGLE,

    SOLLYA_BASE_FUNC_SINH,

    SOLLYA_BASE_FUNC_SQRT,

    SOLLYA_BASE_FUNC_SUB,

    SOLLYA_BASE_FUNC_TAN,

    SOLLYA_BASE_FUNC_TANH,

SOLLYA_BASE_FUNC_TRIPLEDOUBLE) = map(int,xrange(44))

print "\nSuperficial pobyso check..."    

print "First constant - SOLLYA_BASE_FUNC_ABS: ", SOLLYA_BASE_FUNC_ABS

print "Last constant  - SOLLYA_BASE_FUNC_TRIPLEDOUBLE: ", SOLLYA_BASE_FUNC_TRIPLEDOUBLE

pobyso_max_arity = 9

def pobyso_absolute_so_so():

    return(sollya_lib_absolute(None))

def pobyso_autoprint(arg):

    sollya_lib_autoprint(arg,None)

def pobyso_autoprint_so_so(arg):

    sollya_lib_autoprint(arg,None)

def pobyso_bounds_to_range_sa_so(rnLowerBoundSa, rnUpperBoundSa, \

                                 precisionSa=None):

    """

    Return a Sollya range from to 2 RealField Sage elements.

    The Sollya range element has a sufficient precision to hold all

    the digits of the widest of the Sage bounds.

    """

    # Sanity check.

    if rnLowerBoundSa > rnUpperBoundSa:

        return None

    # Precision stuff.

    if precisionSa is None:

        # Check for the largest precision.

        lbPrecSa = rnLowerBoundSa.parent().precision()

        ubPrecSa = rnLowerBoundSa.parent().precision()

        maxPrecSa = max(lbPrecSa, ubPrecSa)

    else:

        maxPrecSa = precisionSa

    # From Sage to Sollya bounds.

#    lowerBoundSo = sollya_lib_constant(get_rn_value(rnLowerBoundSa),

#                                       maxPrecSa)

    lowerBoundSo = pobyso_constant_sa_so(rnLowerBoundSa,

                                         maxPrecSa)

    upperBoundSo = pobyso_constant_sa_so(rnUpperBoundSa,

                                       maxPrecSa)

    # From Sollya bounds to range.

    rangeSo = sollya_lib_range(lowerBoundSo, upperBoundSo)

    # Back to original precision.

    # Clean up

    sollya_lib_clear_obj(lowerBoundSo)

    sollya_lib_clear_obj(upperBoundSo)

    return rangeSo

# End pobyso_bounds_to_range_sa_so

def pobyso_build_end_elliptic_list_so_so(*args):

    """

    From argumrny Sollya objects, create a Sollya end elliptic list.

    Elements of the list are "eaten" (should not be cleared individualy, 

    are cleared when the list is cleared).

    """

    if len(args) == 0:

        ## Called with an empty list produced "error".

        return sollya_lib_build_end_elliptic_list(None)

    index = 0

    ## One can not append elements to an elliptic list, prepend only is 

    #  permitted.

    for argument in reversed(args):

        if index == 0:

            listSo = sollya_lib_build_end_elliptic_list(argument, None)

        else:

            listSo = sollya_lib_prepend(argument, listSo)

        index += 1

    return listSo

# End pobyso_build_end_elliptic_list_so_so

def pobyso_build_function_sub_so_so(exp1So, exp2So):

    return(sollya_lib_build_function_sub(exp1So, exp2So))

def pobyso_change_var_in_function_so_so(funcSo, chvarExpSo):

    """

    Variable change in a function.

    """

    return(sollya_lib_evaluate(funcSo,chvarExpSo))

# End pobyso_change_var_in_function_so_so     

def pobyso_chebyshevform_so_so(functionSo, degreeSo, intervalSo):

    resultSo = sollya_lib_chebyshevform(functionSo, degreeSo, intervalSo)

    return(resultSo)

# End pobyso_chebyshevform_so_so.

def pobyso_clear_taylorform_sa_so(taylorFormSaSo):

    """

    This method is necessary to correctly clean up the memory from Taylor forms.

    These are made of a Sollya object, a Sollya object list, a Sollya object.

    For no clearly understood reason, sollya_lib_clear_object_list crashed 

    when applied to the object list.

    Here, we decompose it into Sage list of Sollya objects references and we

     clear them one by one. 

    """

    sollya_lib_clear_obj(taylorFormSaSo[0])

    (coefficientsErrorsListSaSo, numElementsSa, isEndEllipticSa) = \

        pobyso_get_list_elements_so_so(taylorFormSaSo[1])

    for element in coefficientsErrorsListSaSo:

        sollya_lib_clear_obj(element)

    sollya_lib_clear_obj(taylorFormSaSo[1])

    sollya_lib_clear_obj(taylorFormSaSo[2])

# End pobyso_clear_taylorform_sa_so 

def pobyso_cmp(rnArgSa, cteSo):

    """

    Compare the MPFR value a RealNumber with that of a Sollya constant.

    Get the value of the Sollya constant into a RealNumber and compare

    using MPFR. Could be optimized by working directly with a mpfr_t

    for the intermediate number. 

    """

    # Get the precision of the Sollya constant to build a Sage RealNumber 

    # with enough precision.to hold it.

    precisionOfCte = c_int(0)

    # From the Sollya constant, create a local Sage RealNumber.

    sollya_lib_get_prec_of_constant(precisionOfCte, cteSo) 

    #print "Precision of constant: ", precisionOfCte

    RRRR = RealField(precisionOfCte.value)

    rnLocalSa = RRRR(0)

    sollya_lib_get_constant(get_rn_value(rnLocalSa), cteSo)

    #

    ## Compare the Sage RealNumber version of the Sollya constant with rnArg.

    return(cmp_rn_value(rnArgSa, rnLocal))

# End pobyso_smp

def pobyso_compute_pos_function_abs_val_bounds_sa_sa(funcSa, lowerBoundSa, \

                                                     upperBoundSa):

    """

    TODO: completely rework and test.

    """

    pobyso = pobyso_name_free_variable_sa_so(funcSa.variables()[0])

    funcSo = pobyso_parse_string(funcSa._assume_str().replace('_SAGE_VAR_', ''))

    rangeSo = pobyso_range_sa_so(lowerBoundSa, upperBoundSa)

    infnormSo = pobyso_infnorm_so_so(funcSo,rangeSo)

    # Sollya return the infnorm as an interval.

    fMaxSa = pobyso_get_interval_from_range_so_sa(infnormSo)

    # Get the top bound and compute the binade top limit.

    fMaxUpperBoundSa = fMaxSa.upper()

    binadeTopLimitSa = 2**ceil(fMaxUpperBoundSa.log2())

    # Put up together the function to use to compute the lower bound.

    funcAuxSo = pobyso_parse_string(str(binadeTopLimitSa) +  \

                                    '-(' + f._assume_str().replace('_SAGE_VAR_', '') + ')')

    pobyso_autoprint(funcAuxSo)

    # Clear the Sollya range before a new call to infnorm and issue the call.

    sollya_lib_clear_obj(infnormSo)

    infnormSo = pobyso_infnorm_so_so(funcAuxSo,rangeSo)

    fMinSa = pobyso_get_interval_from_range_so_sa(infnormSo)

    sollya_lib_clear_obj(infnormSo)

    fMinLowerBoundSa = binadeTopLimitSa - fMinSa.lower()

    # Compute the maximum of the precisions of the different bounds.

    maxPrecSa = max([fMinLowerBoundSa.parent().precision(), \

                     fMaxUpperBoundSa.parent().precision()])

    # Create a RealIntervalField and create an interval with the "good" bounds.

    RRRI = RealIntervalField(maxPrecSa)

    imageIntervalSa = RRRI(fMinLowerBoundSa, fMaxUpperBoundSa)

    # Free the unneeded Sollya objects

    sollya_lib_clear_obj(funcSo)

    sollya_lib_clear_obj(funcAuxSo)

    sollya_lib_clear_obj(rangeSo)

    return(imageIntervalSa)

# End pobyso_compute_pos_function_abs_val_bounds_sa_sa

def pobyso_compute_precision_decay_ratio_function_sa_so():

    """

    Compute the precision decay ratio function for polynomial 

    coefficient progressive trucation.

    """

    functionText = """

    proc(deg, a, b, we, wq) 

    {

      k = we * (exp(x/a)-1) + wq * (b*x)^2 + (1-we-wq) * x;

      return k/k(d);

    };

    """

    return pobyso_parse_string_sa_so(functionText)

# End  pobyso_compute_precision_decay_ratio_function.

def pobyso_constant(rnArg):

    """ Legacy function. See pobyso_constant_sa_so. """

    return(pobyso_constant_sa_so(rnArg))

def pobyso_constant_sa_so(rnArgSa, precisionSa=None):

    """

    Create a Sollya constant from a Sage RealNumber.

    The sollya_lib_constant() function creates a constant

    with the same precision as the source.

    """

    ## Precision stuff. If one wants to change precisions,

    #  everything takes place in Sage. That only makes

    #  sense if one wants to reduce the precision.

    if not precisionSa is None:

        RRR = RealField(precisionSa)

        rnArgSa = RRR(rnArgSa)

    #print rnArgSa, rnArgSa.precision()

    # Sollya constant creation takes place here.

    return sollya_lib_constant(get_rn_value(rnArgSa))

# End pobyso_constant_sa_so

def pobyso_constant_0_sa_so():

    """

    Obvious.

    """

    return pobyso_constant_from_int_sa_so(0)

def pobyso_constant_1():

    """

    Obvious.

    Legacy function. See pobyso_constant_so_so. 

    """

    return pobyso_constant_1_sa_so()

def pobyso_constant_1_sa_so():

    """

    Obvious.

    """

    return(pobyso_constant_from_int_sa_so(1))

def pobyso_constant_from_int(anInt):

    """ Legacy function. See pobyso_constant_from_int_sa_so. """

    return pobyso_constant_from_int_sa_so(anInt)

def pobyso_constant_from_int_sa_so(anInt):

    """

    Get a Sollya constant from a Sage int.

    """

    return sollya_lib_constant_from_int64(long(anInt))

def pobyso_constant_from_int_so_sa(constSo):

    """

    Get a Sage int from a Sollya int constant.

    Usefull for precision or powers in polynomials.

    """

    constSa = c_long(0)

    sollya_lib_get_constant_as_int64(byref(constSa), constSo)

    return constSa.value

# End pobyso_constant_from_int_so_sa

def pobyso_constant_from_mpq_sa_so(rationalSa):

    """

    Make a Sollya constant from Sage rational.

    The Sollya constant is an unevaluated expression.

    Hence no precision argument is needed.

    It is better to leave this way since Sollya has its own

    optimized evaluation mecanism that tries very hard to

    return exact values or at least faithful ones.

    """

    ratExprSo = \

        sollya_lib_constant_from_mpq(sgmp_get_rational_value(rationalSa))

    return ratExprSo

# End pobyso_constant_from_mpq_sa_so.

def pobyso_constant_sollya_prec_sa_so(rnArgSa):

    """

    Create a Sollya constant from a Sage RealNumber at the

    current precision in Sollya.

    """

    currentSollyaPrecSa = pobyso_get_prec_so_sa()

    return pobyso_constant_sa_so(rnArgSa, currentSollyaPrecSa)

# End pobyso_constant_sollya_prec_sa_so

def pobyso_end_elliptic_list_so_sa_so(objectsListSo, intCountSa):

    """

    Create a Sollya end elliptic list made of the objectListSo[0] to

     objectsListSo[intCountSa-1] objects.

    """     

    return sollya_lib_end_elliptic_list(objectSo, int(intCountSa))

def pobyso_error_so():

    return sollya_lib_error(None)

# End pobyso_error().

def pobyso_evaluate_so_so(funcSo, argumentSo):

    """

    Evaluates funcSo for arguemntSo through sollya_lib_evaluate().

    """

    return sollya_lib_evaluate(funcSo, argumentSo)

# End pobyso_evaluate_so_so.

def pobyso_float_poly_sa_so(polySa, precSa = None):

    """

    Create a Sollya polynomial from a Sage RealField polynomial.

    """

    ## TODO: filter arguments.

    ## Precision. If a precision is given, convert the polynomial

    #  into the right polynomial field. If not convert it straight

    #  to Sollya.

    if not precSa is None:

        RRR = RealField(precSa)

        ## Create a Sage polynomial in the "right" precision.

        P_RRR = RRR[polySa.variables()[0]]

        polyFloatSa = P_RRR(polySa)

    else:

        polyFloatSa = polySa

        precSa = polySa.parent().base_ring().precision()

    ## Get exponents and coefficients.

    exponentsSa = polyFloatSa.exponents()

    coefficientsSa = polyFloatSa.coefficients()

    ## Build the polynomial.

    polySo = None

    for coefficientSa, exponentSa in zip(coefficientsSa, exponentsSa):

        #print coefficientSa.n(prec=precSa), exponentSa

        coefficientSo = \

            pobyso_constant_sa_so(coefficientSa)

        #pobyso_autoprint(coefficientSo)

        exponentSo = \

            pobyso_constant_from_int_sa_so(exponentSa)

        #pobyso_autoprint(exponentSo)

        monomialSo = sollya_lib_build_function_pow(

                       sollya_lib_build_function_free_variable(),

                       exponentSo)

        if polySo is None:

            polySo = sollya_lib_build_function_mul(coefficientSo,

                                                   monomialSo)

        else:

            polyTermSo = sollya_lib_build_function_mul(coefficientSo,

                                                       monomialSo)

            polySo = sollya_lib_build_function_add(polySo, polyTermSo)

    return polySo

# End pobyso_float_poly_sa_so    

def pobyso_float_poly_so_sa(polySo, realFieldSa=None):

    """

    Convert a Sollya polynomial into a Sage floating-point polynomial.

    We assume that the polynomial is in canonical form.

    If no realField is given, a RealField corresponding to the maximum 

    precision of the coefficients is internally computed.

    The real field is not returned but can be easily retrieved from 

    the polynomial itself.

    ALGORITHM:

    - (optional) compute the RealField of the coefficients;

    - convert the Sollya expression into a Sage expression;

    - convert the Sage expression into a Sage polynomial

    TODO: the canonical thing for the polynomial.

    """    

    if realFieldSa is None:

        expressionPrecSa = pobyso_get_max_prec_of_exp_so_sa(polySo)

        print "Maximum precision of Sollya polynomial coefficients:", expressionPrecSa

        realFieldSa      = RealField(expressionPrecSa)

    #print "Sollya expression before...",

    #pobyso_autoprint(polySo)

    expressionSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo,

                                                             realFieldSa)

    #print "...Sollya expression after.",

    #pobyso_autoprint(polySo)

    polyVariableSa = expressionSa.variables()[0]

    polyRingSa     = realFieldSa[str(polyVariableSa)]

    #print polyRingSa

    # Do not use the polynomial(expressionSa, ring=polyRingSa) form!

    polynomialSa = polyRingSa(expressionSa)

    polyCoeffsListSa = polynomialSa.coefficients()

    #for coeff in polyCoeffsListSa:

    #    print coeff.abs().n()

    return polynomialSa

# End pobyso_float_poly_so_sa

def pobyso_free_variable():

    """

    Ultra thin wrapper around the sollya_lib_function_build_free_variable function.

    """

    return sollya_lib_build_function_free_variable()

def pobyso_function_type_as_string(funcType):

    """ Legacy function. See pobyso_function_type_as_string_so_sa. """

    return(pobyso_function_type_as_string_so_sa(funcType))

def pobyso_function_type_as_string_so_sa(funcType):

    """

    Numeric Sollya function codes -> Sage mathematical function names.

    Notice that pow -> ^ (a la Sage, not a la Python).

    """

    if funcType == SOLLYA_BASE_FUNC_ABS:

        return "abs"

    elif funcType == SOLLYA_BASE_FUNC_ACOS:

        return "arccos"

    elif funcType == SOLLYA_BASE_FUNC_ACOSH:

        return "arccosh"

    elif funcType == SOLLYA_BASE_FUNC_ADD:

        return "+"

    elif funcType == SOLLYA_BASE_FUNC_ASIN:

        return "arcsin"

    elif funcType == SOLLYA_BASE_FUNC_ASINH:

        return "arcsinh"

    elif funcType == SOLLYA_BASE_FUNC_ATAN:

        return "arctan"

    elif funcType == SOLLYA_BASE_FUNC_ATANH:

        return "arctanh"

    elif funcType == SOLLYA_BASE_FUNC_CEIL:

        return "ceil"

    elif funcType == SOLLYA_BASE_FUNC_CONSTANT:

        return "cte"

    elif funcType == SOLLYA_BASE_FUNC_COS:

        return "cos"

    elif funcType == SOLLYA_BASE_FUNC_COSH:

        return "cosh"

    elif funcType == SOLLYA_BASE_FUNC_DIV:

        return "/"

    elif funcType == SOLLYA_BASE_FUNC_DOUBLE:

        return "double"

    elif funcType == SOLLYA_BASE_FUNC_DOUBLEDOUBLE:

        return "doubleDouble"

    elif funcType == SOLLYA_BASE_FUNC_DOUBLEEXTENDED:

        return "doubleDxtended"

    elif funcType == SOLLYA_BASE_FUNC_ERF:

        return "erf"

    elif funcType == SOLLYA_BASE_FUNC_ERFC:

        return "erfc"

    elif funcType == SOLLYA_BASE_FUNC_EXP:

        return "exp"

    elif funcType == SOLLYA_BASE_FUNC_EXP_M1:

        return "expm1"

    elif funcType == SOLLYA_BASE_FUNC_FLOOR:

        return "floor"

    elif funcType == SOLLYA_BASE_FUNC_FREE_VARIABLE:

        return "freeVariable"

    elif funcType == SOLLYA_BASE_FUNC_HALFPRECISION:

        return "halfPrecision"

    elif funcType == SOLLYA_BASE_FUNC_LIBRARYCONSTANT:

        return "libraryConstant"

    elif funcType == SOLLYA_BASE_FUNC_LIBRARYFUNCTION:

        return "libraryFunction"

    elif funcType == SOLLYA_BASE_FUNC_LOG:

        return "log"

    elif funcType == SOLLYA_BASE_FUNC_LOG_10:

        return "log10"

    elif funcType == SOLLYA_BASE_FUNC_LOG_1P:

        return "log1p"

    elif funcType == SOLLYA_BASE_FUNC_LOG_2:

        return "log2"

    elif funcType == SOLLYA_BASE_FUNC_MUL:

        return "*"

    elif funcType == SOLLYA_BASE_FUNC_NEARESTINT:

        return "round"

    elif funcType == SOLLYA_BASE_FUNC_NEG:

        return "__neg__"

    elif funcType == SOLLYA_BASE_FUNC_PI:

        return "pi"

    elif funcType == SOLLYA_BASE_FUNC_POW:

        return "^"

    elif funcType == SOLLYA_BASE_FUNC_PROCEDUREFUNCTION:

        return "procedureFunction"

    elif funcType == SOLLYA_BASE_FUNC_QUAD:

        return "quad"

    elif funcType == SOLLYA_BASE_FUNC_SIN:

        return "sin"

    elif funcType == SOLLYA_BASE_FUNC_SINGLE:

        return "single"

    elif funcType == SOLLYA_BASE_FUNC_SINH:

        return "sinh"

    elif funcType == SOLLYA_BASE_FUNC_SQRT:

        return "sqrt"

    elif funcType == SOLLYA_BASE_FUNC_SUB:

        return "-"

    elif funcType == SOLLYA_BASE_FUNC_TAN:

        return "tan"

    elif funcType == SOLLYA_BASE_FUNC_TANH:

        return "tanh"

    elif funcType == SOLLYA_BASE_FUNC_TRIPLEDOUBLE:

        return "tripleDouble"

    else:

        return None

def pobyso_get_constant(rnArgSa, constSo):

    """ Legacy function. See pobyso_get_constant_so_sa. """

    return pobyso_get_constant_so_sa(rnArgSa, constSo)

# End pobyso_get_constant

def pobyso_get_constant_so_sa(rnArgSa, constSo):

    """

    Set the value of rnArgSo to the value of constSo in MPFR_RNDN mode.

    rnArg must already exist and belong to some RealField.

    We assume that constSo points to a Sollya constant.

    """

    outcome = sollya_lib_get_constant(get_rn_value(rnArgSa), constSo)

    if outcome == 0: # Failure because constSo is not a constant expression.

        return None

    else:

        return outcome

# End  pobyso_get_constant_so_sa

def pobyso_get_constant_as_rn(ctExpSo):

    """ 

    Legacy function. See pobyso_get_constant_as_rn_so_sa. 

    """ 

    return(pobyso_get_constant_as_rn_so_sa(ctExpSo))

def pobyso_get_constant_as_rn_so_sa(constExpSo):

    """

    Get a Sollya constant as a Sage "real number".

    The precision of the floating-point number returned is that of the Sollya

    constant.

    """

    precisionSa  = pobyso_get_prec_of_constant_so_sa(constExpSo)

    ## If the expression can not be exactly converted, None is returned.

    #  In this case opt for the Sollya current expression.

    if precisionSa is None:

        precisionSa = pobyso_get_prec_so_sa()

    RRRR = RealField(precisionSa)

    rnSa = RRRR(0)

    outcome = sollya_lib_get_constant(get_rn_value(rnSa), constExpSo)

    if outcome == 0:

        return None

    else:

        return rnSa

# End pobyso_get_constant_as_rn_so_sa

def pobyso_get_constant_as_rn_with_rf(ctExp, realField):

    """ 

    Legacy function. See pobyso_get_constant_as_rn_with_rf_so_sa.

    """

    return pobyso_get_constant_as_rn_with_rf_so_sa(ctExp, realField)

# End pobyso_get_constant_as_rn_with_rf

def pobyso_get_constant_as_rn_with_rf_so_sa(ctExpSo, realFieldSa = None):

    """

    Get a Sollya constant as a Sage "real number".

    If no real field is specified, the precision of the floating-point number 

    returned is that of the Sollya constant.

    Otherwise is is that of the real field. Hence rounding may happen.

    """

    if realFieldSa is None:

        return pobyso_get_constant_as_rn_so_sa(ctExpSo)

    rnSa = realFieldSa(0)

    outcome = sollya_lib_get_constant(get_rn_value(rnSa), ctExpSo)

    if outcome == 0:

        return None

    else:

        return rnSa

# End pobyso_get_constant_as_rn_with_rf_so_sa

def pobyso_get_free_variable_name():

    """ 

    Legacy function. See pobyso_get_free_variable_name_so_sa.

    """

    return(pobyso_get_free_variable_name_so_sa())

def pobyso_get_free_variable_name_so_sa():

    return sollya_lib_get_free_variable_name()

def pobyso_get_function_arity(expressionSo):

    """ 

    Legacy function. See pobyso_get_function_arity_so_sa.

    """

    return(pobyso_get_function_arity_so_sa(expressionSo))

def pobyso_get_function_arity_so_sa(expressionSo):

    arity = c_int(0)

    sollya_lib_get_function_arity(byref(arity),expressionSo)

    return int(arity.value)

def pobyso_get_head_function(expressionSo):

    """ 

    Legacy function. See pobyso_get_head_function_so_sa. 

    """

    return(pobyso_get_head_function_so_sa(expressionSo)) 

def pobyso_get_head_function_so_sa(expressionSo):

    functionType = c_int(0)

    sollya_lib_get_head_function(byref(functionType), expressionSo, None)

    return int(functionType.value)

def pobyso_get_interval_from_range_so_sa(soRange, realIntervalFieldSa = None ):

    """

    Return the Sage interval corresponding to the Sollya range argument.

    If no reaIntervalField is passed as an argument, the interval bounds are not

    rounded: they are elements of RealIntervalField of the "right" precision

    to hold all the digits.

    """

    prec = c_int(0)

    if realIntervalFieldSa is None:

        retval = sollya_lib_get_prec_of_range(byref(prec), soRange, None)

        if retval == 0:

            return None

        realIntervalFieldSa = RealIntervalField(prec.value)

    intervalSa = realIntervalFieldSa(0,0)

    retval = \

        sollya_lib_get_interval_from_range(get_interval_value(intervalSa),\

                                           soRange)

    if retval == 0:

        return None

    return intervalSa

# End pobyso_get_interval_from_range_so_sa

def pobyso_get_list_elements(soObj):

    """ Legacy function. See pobyso_get_list_elements_so_so. """

    return pobyso_get_list_elements_so_so(soObj)

def pobyso_get_list_elements_so_so(objectListSo):

    """

    Get the Sollya list elements as a Sage/Python array of Sollya objects.

    INPUT:

    - objectListSo: a Sollya list of Sollya objects.

    OUTPUT:

    - a Sage/Python tuple made of:

      - a Sage/Python list of Sollya objects,

      - a Sage/Python int holding the number of elements,

      - a Sage/Python int stating (!= 0) that the list is end-elliptic.

    NOTE::

        We recover the addresses of the Sollya object from the list of pointers

        returned by sollya_lib_get_list_elements. The list itself is freed.

    TODO::

        Figure out what to do with numElements since the number of elements

        can easily be recovered from the list itself. 

        Ditto for isEndElliptic.

    """

    listAddress = POINTER(c_longlong)()

    numElements = c_int(0)

    isEndElliptic = c_int(0)

    listAsSageList = []

    result = sollya_lib_get_list_elements(byref(listAddress),\

                                          byref(numElements),\

                                          byref(isEndElliptic),\

                                          objectListSo)

    if result == 0 :

        return None

    for i in xrange(0, numElements.value, 1):

       #listAsSageList.append(sollya_lib_copy_obj(listAddress[i]))

       listAsSageList.append(listAddress[i])

       # Clear each of the elements returned by Sollya.

       #sollya_lib_clear_obj(listAddress[i])

    # Free the list itself.   

    sollya_lib_free(listAddress)

    return (listAsSageList, numElements.value, isEndElliptic.value)

def pobyso_get_max_prec_of_exp(soExp):

    """ Legacy function. See pobyso_get_max_prec_of_exp_so_sa. """

    return pobyso_get_max_prec_of_exp_so_sa(soExp)

def pobyso_get_max_prec_of_exp_so_sa(expSo):

    """

    Get the maximum precision used for the numbers in a Sollya expression.

    Arguments:

    soExp -- a Sollya expression pointer

    Return value:

    A Python integer

    TODO: 

    - error management;

    - correctly deal with numerical type such as DOUBLEEXTENDED.

    """

    maxPrecision = 0

    minConstPrec = 0

    currentConstPrec = 0

    operator = pobyso_get_head_function_so_sa(expSo)

    if (operator != SOLLYA_BASE_FUNC_CONSTANT) and \

    (operator != SOLLYA_BASE_FUNC_FREE_VARIABLE):

        (arity, subexpressions) = pobyso_get_subfunctions_so_sa(expSo)

        for i in xrange(arity):

            maxPrecisionCandidate = \

                pobyso_get_max_prec_of_exp_so_sa(subexpressions[i])

            if maxPrecisionCandidate > maxPrecision:

                maxPrecision = maxPrecisionCandidate

        return maxPrecision

    elif operator == SOLLYA_BASE_FUNC_CONSTANT:

        #minConstPrec = pobyso_get_min_prec_of_constant_so_sa(expSo)

        #currentConstPrec = pobyso_get_min_prec_of_constant_so_sa(soExp)

        #print minConstPrec, " - ", currentConstPrec 

        return pobyso_get_min_prec_of_constant_so_sa(expSo)

    elif operator == SOLLYA_BASE_FUNC_FREE_VARIABLE:

        return 0

    else:

        print "pobyso_get_max_prec_of_exp_so_sa: unexepected operator."

        return 0

def pobyso_get_min_prec_of_constant_so_sa(constExpSo):

    """

    Get the minimum precision necessary to represent the value of a Sollya

    constant.

    MPFR_MIN_PREC and powers of 2 are taken into account.

    We assume that constExpSo is a pointer to a Sollay constant expression.

    """

    constExpAsRnSa = pobyso_get_constant_as_rn_so_sa(constExpSo)

    return(min_mpfr_size(get_rn_value(constExpAsRnSa)))

def pobyso_get_poly_so_sa(polySo, realFieldSa=None):

    """

    Convert a Sollya polynomial into a Sage polynomial.

    Legacy function. Use pobyso_float_poly_so_sa() instead.

    """    

    return pobyso_float_poly_so_sa(polySo,realFieldSa)

# End pobyso_get_poly_so_sa

def pobyso_get_prec():

    """ Legacy function. See pobyso_get_prec_so_sa(). """

    return pobyso_get_prec_so_sa()

def pobyso_get_prec_so():

    """

    Get the current default precision in Sollya.

    The return value is a Sollya object.

    Usefull when modifying the precision back and forth by avoiding

    extra conversions.

    """

    return sollya_lib_get_prec(None)

def pobyso_get_prec_so_sa():

    """

    Get the current default precision in Sollya.

    The return value is Sage/Python int.

    """

    precSo = sollya_lib_get_prec(None)

    precSa = c_int(0)

    sollya_lib_get_constant_as_int(byref(precSa), precSo)

    sollya_lib_clear_obj(precSo)

    return int(precSa.value)

# End pobyso_get_prec_so_sa.

def pobyso_get_prec_so_so_sa():

    """

    Return the current precision both as a Sollya object and a

    Sage integer as hybrid tuple.

    To avoid multiple calls for precision manipulations. 

    """

    precSo = sollya_lib_get_prec(None)

    precSa = c_int(0)

    sollya_lib_get_constant_as_int(byref(precSa), precSo)

    return (precSo, precSa)

def pobyso_get_prec_of_constant(ctExpSo):

    """ Legacy function. See pobyso_get_prec_of_constant_so_sa. """

    return pobyso_get_prec_of_constant_so_sa(ctExpSo)

def pobyso_get_prec_of_constant_so_sa(ctExpSo):

    """

    Tries to find a precision to represent ctExpSo without rounding.

    If not possible, returns None.

    """

    prec = c_int(0)

    retc = sollya_lib_get_prec_of_constant(byref(prec), ctExpSo, None)

    if retc == 0:

        return None

    return int(prec.value)

def pobyso_get_prec_of_range_so_sa(rangeSo):

    """

    Returns the number of bits elements of a range are coded with.

    """

    prec = c_int(0)

    retc = sollya_lib_get_prec_of_range(byref(prec), rangeSo, None)

    if retc == 0:

        return(None)

    return int(prec.value)

# End pobyso_get_prec_of_range_so_sa()

def pobyso_get_sage_exp_from_sollya_exp(sollyaExpSo, realField = RR):

    """ Legacy function. See pobyso_get_sage_exp_from_sollya_exp_so_sa. """

    return pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, 

                                                     realField = RR)

def pobyso_get_sage_exp_from_sollya_exp_so_sa(sollyaExpSo, realFieldSa = RR):

    """

    Get a Sage expression from a Sollya expression. 

    Currently only tested with polynomials with floating-point coefficients.

    Notice that, in the returned polynomial, the exponents are RealNumbers.

    """

    #pobyso_autoprint(sollyaExp)

    operatorSa = pobyso_get_head_function_so_sa(sollyaExpSo)

    sollyaLibFreeVariableName = sollya_lib_get_free_variable_name()

    ## Get rid of the "_"'s in "_x_", if any.

    sollyaLibFreeVariableName = re.sub('_', '', sollyaLibFreeVariableName)

    # Constants and the free variable are special cases.

    # All other operator are dealt with in the same way.

    if (operatorSa != SOLLYA_BASE_FUNC_CONSTANT) and \

       (operatorSa != SOLLYA_BASE_FUNC_FREE_VARIABLE):

        (aritySa, subexpressionsSa) = pobyso_get_subfunctions_so_sa(sollyaExpSo)

        if aritySa == 1:

            sageExpSa = eval(pobyso_function_type_as_string_so_sa(operatorSa) + \

            "(" + pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], \

            realFieldSa) + ")")

        elif aritySa == 2:

            # We do not get through the preprocessor.

            # The "^" operator is then a special case.

            if operatorSa == SOLLYA_BASE_FUNC_POW:

                operatorAsStringSa = "**"

            else:

                operatorAsStringSa = \

                    pobyso_function_type_as_string_so_sa(operatorSa)

            sageExpSa = \

              eval("pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[0], realFieldSa)"\

              + " " + operatorAsStringSa + " " + \

                   "pobyso_get_sage_exp_from_sollya_exp_so_sa(subexpressionsSa[1], realFieldSa)")

        # We do not know yet how to deal with arity >= 3 

        # (is there any in Sollya anyway?).

        else:

            sageExpSa = eval('None')

        return sageExpSa

    elif operatorSa == SOLLYA_BASE_FUNC_CONSTANT:

        #print "This is a constant"

        return pobyso_get_constant_as_rn_with_rf_so_sa(sollyaExpSo, realFieldSa)

    elif operatorSa == SOLLYA_BASE_FUNC_FREE_VARIABLE:

        #print "This is free variable"

        return eval(sollyaLibFreeVariableName)

    else:

        print "Unexpected"

        return eval('None')

# End pobyso_get_sage_exp_from_sollya_exp_so_sa

def pobyso_get_subfunctions(expressionSo):

    """ Legacy function. See pobyso_get_subfunctions_so_sa. """

    return pobyso_get_subfunctions_so_sa(expressionSo) 

# End pobyso_get_subfunctions.

def pobyso_get_subfunctions_so_sa(expressionSo):

    """

    Get the subfunctions of an expression.

    Return the number of subfunctions and the list of subfunctions addresses.

    S.T.: Could not figure out another way than that ugly list of declarations

    to recover the addresses of the subfunctions.

    We limit ourselves to arity 8 functions. 

    """

    subf0 = c_int(0)

    subf1 = c_int(0)

    subf2 = c_int(0)

    subf3 = c_int(0)

    subf4 = c_int(0)

    subf5 = c_int(0)

    subf6 = c_int(0)

    subf7 = c_int(0)

    subf8 = c_int(0)

    arity = c_int(0)

    nullPtr = POINTER(c_int)()

    sollya_lib_get_subfunctions(expressionSo, byref(arity), \

      byref(subf0), byref(subf1), byref(subf2), byref(subf3), \

      byref(subf4), byref(subf5),\

      byref(subf6), byref(subf7), byref(subf8), nullPtr, None) 

#    byref(cast(subfunctions[0], POINTER(c_int))), \

#    byref(cast(subfunctions[0], POINTER(c_int))), \

#    byref(cast(subfunctions[2], POINTER(c_int))), \

#    byref(cast(subfunctions[3], POINTER(c_int))), \

#    byref(cast(subfunctions[4], POINTER(c_int))), \

#    byref(cast(subfunctions[5], POINTER(c_int))), \

#    byref(cast(subfunctions[6], POINTER(c_int))), \

#    byref(cast(subfunctions[7], POINTER(c_int))), \

#    byref(cast(subfunctions[8], POINTER(c_int))), nullPtr)

    subfunctions = [subf0, subf1, subf2, subf3, subf4, subf5, subf6, subf7, \

                    subf8]

    subs = []

    if arity.value > pobyso_max_arity:

        return(0,[])

    for i in xrange(arity.value):

        subs.append(int(subfunctions[i].value))

        #print subs[i]

    return (int(arity.value), subs)

# End pobyso_get_subfunctions_so_sa

def pobyso_guess_degree_sa_sa(functionSa, intervalSa, approxErrorSa, 

                              weightSa=None, degreeBoundSa=None):

    """

    Sa_sa variant of the solly_guessdegree function.

    Return 0 if something goes wrong.

    """

    functionAsStringSa = functionSa._assume_str().replace('_SAGE_VAR_', '')

    functionSo = pobyso_parse_string_sa_so(functionAsStringSa)

    if pobyso_is_error_so_sa(functionSo):

        sollya_lib_clear_obj(functionSo)

        return 0

    rangeSo = pobyso_interval_to_range_sa_so(intervalSa)

    # The approximation error is expected to be a floating point number.

    if pobyso_is_floating_point_number_sa_sa(approxErrorSa):

        approxErrorSo = pobyso_constant_sa_so(approxErrorSa)

    else:

        approxErrorSo = pobyso_constant_sa_so(RR(approxErrorSa))

    if not weightSa is None:

        weightAsStringSa = weightSa._assume_str().replace('_SAGE_VAR_', '')

        weightSo = pobyso_parse_string_sa_so(weightAsStringSa)

        if pobyso_is_error_so_sa(weightSo):

            sollya_lib_clear_obj(functionSo)

            sollya_lib_clear_obj(rangeSo)

            sollya_lib_clear_obj(approxErrorSo)   

            sollya_lib_clear_obj(weightSo)

            return 0   

    else:

        weightSo = None

    if not degreeBoundSa is None:

        degreeBoundSo = pobyso_constant_from_int_sa_so(degreeBoundSa)

    else:

        degreeBoundSo = None

    guessedDegreeSa = pobyso_guess_degree_so_sa(functionSo,

                                                rangeSo,

                                                approxErrorSo,

                                                weightSo,

                                                degreeBoundSo)

    sollya_lib_clear_obj(functionSo)

    sollya_lib_clear_obj(rangeSo)

    sollya_lib_clear_obj(approxErrorSo)

    if not weightSo is None:

        sollya_lib_clear_obj(weightSo)

    if not degreeBoundSo is None:

        sollya_lib_clear_obj(degreeBoundSo)

    return guessedDegreeSa

# End poyso_guess_degree_sa_sa

def pobyso_guess_degree_so_sa(functionSo, rangeSo, errorSo, weightSo=None, \

                              degreeBoundSo=None):

    """

    Thin wrapper around the guessdegree function.

    Nevertheless, some precision control stuff has been appended.

    """

    # Deal with Sollya internal precision issues: if it is too small, 

    # compared with the error, increases it to about twice -log2(error).

    errorSa = pobyso_get_constant_as_rn_with_rf_so_sa(errorSo)

    log2ErrorSa = errorSa.log2()

    if log2ErrorSa < 0:

        neededPrecisionSa = int(2 * int(-log2ErrorSa) / 64) * 64

    else:

        neededPrecisionSa = int(2 * int(log2ErrorSa) / 64) * 64

    #print "Needed precision:", neededPrecisionSa

    currentPrecSa = pobyso_get_prec_so_sa()

    if neededPrecisionSa > currentPrecSa:

        currentPrecSo = pobyso_get_prec_so()

        pobyso_set_prec_sa_so(neededPrecisionSa)

    #print "Guessing degree..."

    # weightSo and degreeBoundsSo are optional arguments.

    # As declared, sollya_lib_guessdegree must take 5 arguments.

    if weightSo is None:

        degreeRangeSo = sollya_lib_guessdegree(functionSo, rangeSo, errorSo, 

                                               0, 0, None)

    elif degreeBoundSo is None:

        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, \

                                                errorSo, weightSo, 0, None)

    else:

        degreeRangeSo =  sollya_lib_guessdegree(functionSo, rangeSo, errorSo, \

                                                weightSo, degreeBoundSo, None)

    #print "...degree guess done."

    # Restore internal precision, if applicable.

    if neededPrecisionSa > currentPrecSa:

        pobyso_set_prec_so_so(currentPrecSo)

        sollya_lib_clear_obj(currentPrecSo)

    degreeIntervalSa = pobyso_range_to_interval_so_sa(degreeRangeSo)

    sollya_lib_clear_obj(degreeRangeSo)

    # When ok, both bounds match.

    # When the degree bound is too low, the upper bound is the degree

    # for which the error can be honored.

    # When it really goes wrong, the upper bound is infinity. 

    if degreeIntervalSa.lower() == degreeIntervalSa.upper():

        return int(degreeIntervalSa.lower())

    else:

        if degreeIntervalSa.upper().is_infinity():

            return None

        else:

            return int(degreeIntervalSa.upper())

    # End pobyso_guess_degree_so_sa

def pobyso_inf_so_so(intervalSo):

    """

    Very thin wrapper around sollya_lib_inf().

    """

    return sollya_lib_inf(intervalSo)

# End pobyso_inf_so_so.

def pobyso_infnorm_so_so(func, interval, file = None, intervalList = None):

    print "Do not use this function. User pobyso_supnorm_so_so instead."

    return None

def pobyso_interval_to_range_sa_so(intervalSa, precisionSa=None):

    if precisionSa is None:

        precisionSa = intervalSa.parent().precision()

    intervalSo = pobyso_bounds_to_range_sa_so(intervalSa.lower(),\

                                              intervalSa.upper(),\

                                              precisionSa)

    return intervalSo

# End pobyso_interval_to_range_sa_so

def pobyso_is_error_so_sa(objSo):

    """

    Thin wrapper around the sollya_lib_obj_is_error() function.

    """

    if sollya_lib_obj_is_error(objSo) != 0:

        return True

    else:

        return False

# End pobyso_is_error-so_sa

def pobyso_is_floating_point_number_sa_sa(numberSa):

    """

    Check whether a Sage number is floating point.

    Exception stuff added because numbers other than

    floating-point ones do not have the is_real() attribute.

    """

    try:

        return numberSa.is_real()

    except AttributeError:

        return False

# End pobyso_is_floating_piont_number_sa_sa

def pobyso_lib_init():

    sollya_lib_init(None)

def pobyso_lib_close():

    sollya_lib_close(None)

def pobyso_name_free_variable(freeVariableNameSa):

    """ Legacy function. See pobyso_name_free_variable_sa_so. """

    pobyso_name_free_variable_sa_so(freeVariableNameSa)

def pobyso_name_free_variable_sa_so(freeVariableNameSa):

    """

    Set the free variable name in Sollya from a Sage string.

    """

    sollya_lib_name_free_variable(freeVariableNameSa)

def pobyso_parse_string(string):

    """ Legacy function. See pobyso_parse_string_sa_so. """

    return pobyso_parse_string_sa_so(string)

def pobyso_parse_string_sa_so(string):

    """

    Get the Sollya expression computed from a Sage string or

    a Sollya error object if parsing failed.

    """

    return sollya_lib_parse_string(string)

def pobyso_precision_so_sa(ctExpSo):

    """

    Computes the necessary precision to represent a number.

    If x is not zero, it can be uniquely written as x = m · 2e 

    where m is an odd integer and e is an integer. 

    precision(x) returns the number of bits necessary to write m 

    in binary (i.e. ceil(log2(m))).

    """

    #TODO: take care of the special case: 0, @NaN@, @Inf@

    precisionSo = sollya_lib_precision(ctExpSo)

    precisionSa = pobyso_constant_from_int_so_sa(precisionSo)

    sollya_lib_clear_obj(precisionSo)

    return precisionSa

# End pobyso_precision_so_sa

def pobyso_polynomial_coefficients_progressive_round_so_so(polySo,

                                                           funcSo,

                                                           icSo,

                                                           intervalSo,

                                                           itpSo,

                                                           ftpSo,

                                                           maxPrecSo,

                                                           maxErrSo):

    print "Input arguments:"

    pobyso_autoprint(polySo)

    pobyso_autoprint(funcSo)

    pobyso_autoprint(icSo)

    pobyso_autoprint(intervalSo)

    pobyso_autoprint(itpSo)

    pobyso_autoprint(ftpSo)

    pobyso_autoprint(maxPrecSo)

    pobyso_autoprint(maxErrSo)

    print "________________"

    ## Higher order function see: 

    #  http://effbot.org/pyfaq/how-do-you-make-a-higher-order-function-in-python.htm

    def precision_decay_ratio_function(degreeSa):

        def outer(x):

            def inner(x):

                we = 3/8

                wq = 2/8

                a  = 2.2

                b  = 2 

                return we*(exp(x/a)-1) +  wq*((b*x)**2) + (1-we-wq)*x

            return  inner(x)/inner(degreeSa)

        return outer