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

"""

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_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_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.

    """

    # Precision stuff

    if precisionSa is None:

        precisionSa = rnArgSa.parent().precision()

    currentSollyaPrecisionSo = sollya_lib_get_prec()

    currentSollyaPrecisionSa = \

        pobyso_constant_from_int(currentSollyaPrecisionSo)

    # Sollya constant creation takes place here.

    if precisionSa > currentSollyaPrecisionSa:

        pobyso_set_prec_sa_so(precisionSa)

        constantSo = sollya_lib_constant(get_rn_value(rnArgSa))

        pobyso_set_prec_so_so(currentSollyaPrecisionSo)

    else:

        constantSo = sollya_lib_constant(get_rn_value(rnArgSa))

    sollya_lib_clear_obj(currentSollyaPrecisionSo)

    return constantSo

# 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_int(int(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_int(0)

    sollya_lib_get_constant_as_int(byref(constSa), constSo)

    return(constSa.value)

# End pobyso_constant_from_int_so_sa

def pobyso_error_so():

    return sollya_lib_error(None)

# End pobyso_error().

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))

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.

    """

    return(sollya_lib_get_constant(get_rn_value(rnArgSa), constSo))

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) 

    RRRR = RealField(precisionSa)

    rnSa = RRRR(0)

    sollya_lib_get_constant(get_rn_value(rnSa), constExpSo)

    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))

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:

        sollyaPrecSa = pobyso_get_prec_of_constant_so_sa(ctExpSo)

        realFieldSa = RealField(sollyaPrecSa)

    rnSa = realFieldSa(0)

    sollya_lib_get_constant(get_rn_value(rnSa), ctExpSo)

    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 point

    """

    constExpAsRnSa = pobyso_get_constant_as_rn_so_sa(constExpSo)

    return(min_mpfr_size(get_rn_value(constExpAsRnSa)))

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()

    # 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_poly_from_sollya_poly

def pobyso_get_poly_sa_so(polySo, realFieldSa=None):

    """

    Create a Sollya polynomial from a Sage polynomial.

    """

    pass

# pobyso_get_poly_sa_so

def pobyso_get_poly_so_sa(polySo, realFieldSa=None):

    """

    Convert a Sollya polynomial into a Sage 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)

        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)

    return(polynomialSa)

# End pobyso_get_sage_poly_from_sollya_poly

def pobyso_get_subfunctions(expressionSo):

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

    return(pobyso_get_subfunctions_so_sa(expressionSo)) 

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)

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_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):

    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):

    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_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_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

    """

    return numberSa.parent().__class__ is RR.__class__

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_range(rnLowerBound, rnUpperBound):

    """ Legacy function. See pobyso_range_sa_so. """

    return(pobyso_range_sa_so(rnLowerBound, rnUpperBound)) 

def pobyso_range_to_interval_so_sa(rangeSo, realIntervalFieldSa = None):

    """

    Get a Sage interval from a Sollya range.

    If no realIntervalField is given as a parameter, the Sage interval

    precision is that of the Sollya range.

    Otherwise, the precision is that of the realIntervalField. In this case

    rounding may happen.

    """

    if realIntervalFieldSa is None:

        precSa = pobyso_get_prec_of_range_so_sa(rangeSo)

        realIntervalFieldSa = RealIntervalField(precSa)

    intervalSa = \

        pobyso_get_interval_from_range_so_sa(rangeSo, realIntervalFieldSa) 

    return(intervalSa)

def pobyso_remez_canonical_sa_sa(func, \

                                 degree, \

                                 lowerBound, \

                                 upperBound, \

                                 weight = None, \

                                 quality = None):

    """

    All arguments are Sage/Python.

    The functions (func and weight) must be passed as expressions or strings.

    Otherwise the function fails. 

    The return value is a Sage polynomial.

    """

    var('zorglub')    # Dummy variable name for type check only. Type of 

    # zorglub is "symbolic expression".

    polySo = pobyso_remez_canonical_sa_so(func, \

                                 degree, \

                                 lowerBound, \

                                 upperBound, \

                                 weight, \

                                 quality)

    # String test

    if parent(func) == parent("string"):

        functionSa = eval(func)

    # Expression test.

    elif type(func) == type(zorglub):

        functionSa = func

    else:

        return None

    #

    maxPrecision = 0

    if polySo is None:

        return(None)

    maxPrecision = pobyso_get_max_prec_of_exp_so_sa(polySo)

    RRRRSa = RealField(maxPrecision)

    polynomialRingSa = RRRRSa[functionSa.variables()[0]]

    expSa = pobyso_get_sage_exp_from_sollya_exp_so_sa(polySo, RRRRSa)

    polySa = polynomial(expSa, polynomialRingSa)

    sollya_lib_clear_obj(polySo)

    return(polySa)

# End pobyso_remez_canonical_sa_sa

def pobyso_remez_canonical(func, \

                           degree, \

                           lowerBound, \

                           upperBound, \

                           weight = "1", \

                           quality = None):

    """ Legacy function. See pobyso_remez_canonical_sa_so. """

    return(pobyso_remez_canonical_sa_so(func, \

                                        degree, \

                                        lowerBound, \

                                        upperBound, \

                                        weight, \

                                        quality))

def pobyso_remez_canonical_sa_so(func, \

                                 degree, \

                                 lowerBound, \

                                 upperBound, \

                                 weight = None, \

                                 quality = None):

    """

    All arguments are Sage/Python.

    The functions (func and weight) must be passed as expressions or strings.

    Otherwise the function fails. 

    The return value is a pointer to a Sollya function.

    """

    var('zorglub')    # Dummy variable name for type check only. Type of

    # zorglub is "symbolic expression".

    currentVariableNameSa = None

    # The func argument can be of different types (string, 

    # symbolic expression...)

    if parent(func) == parent("string"):

        localFuncSa = eval(func)

        if len(localFuncSa.variables()) > 0:

            currentVariableNameSa = localFuncSa.variables()[0]

            sollya_lib_name_free_variable(str(currentVariableNameSa))

            functionSo = \

              sollya_lib_parse_string(localFuncSa._assume_str().replace('_SAGE_VAR_', ''))

    # Expression test.

    elif type(func) == type(zorglub):

        # Until we are able to translate Sage expressions into Sollya 

        # expressions : parse the string version.

        if len(func.variables()) > 0:

            currentVariableNameSa = func.variables()[0]

            sollya_lib_name_free_variable(str(currentVariableNameSa))

            functionSo = \

              sollya_lib_parse_string(func._assume_str().replace('_SAGE_VAR_', ''))

    else:

        return(None)

    if weight is None: # No weight given -> 1.

        weightSo = pobyso_constant_1_sa_so()

    elif parent(weight) == parent("string"): # Weight given as string: parse it.

        weightSo = sollya_lib_parse_string(func)

    elif type(weight) == type(zorglub): # Weight given as symbolice expression.

        functionSo = \

          sollya_lib_parse_string_sa_so(weight._assume_str().replace('_SAGE_VAR_', ''))

    else:

        return(None)

    degreeSo = pobyso_constant_from_int(degree)

    rangeSo = pobyso_bounds_to_range_sa_so(lowerBound, upperBound)

    if not quality is None:

        qualitySo= pobyso_constant_sa_so(quality)

    else:

        qualitySo = None

    remezPolySo = sollya_lib_remez(functionSo, \

                                   degreeSo, \

                                   rangeSo, \

                                   weightSo, \

                                   qualitySo, \

                                   None)

    sollya_lib_clear_obj(functionSo)

    sollya_lib_clear_obj(degreeSo)

    sollya_lib_clear_obj(rangeSo)

    sollya_lib_clear_obj(weightSo)

    if not qualitySo is None:

        sollya_lib_clear_obj(qualitySo)

    return(remezPolySo)

# End pobyso_remez_canonical_sa_so

def pobyso_remez_canonical_so_so(funcSo, \

                                 degreeSo, \

                                 rangeSo, \

                                 weightSo = pobyso_constant_1_sa_so(),\

                                 qualitySo = None):

    """

    All arguments are pointers to Sollya objects.

    The return value is a pointer to a Sollya function.

    """

    if not sollya_lib_obj_is_function(funcSo):

        return(None)

    return(sollya_lib_remez(funcSo, degreeSo, rangeSo, weightSo, qualitySo, None))

def pobyso_set_canonical_off():

    sollya_lib_set_canonical(sollya_lib_off())

def pobyso_set_canonical_on():

    sollya_lib_set_canonical(sollya_lib_on())

def pobyso_set_prec(p):

    """ Legacy function. See pobyso_set_prec_sa_so. """

    pobyso_set_prec_sa_so(p)

def pobyso_set_prec_sa_so(p):

    a = c_int(p)

    precSo = c_void_p(sollya_lib_constant_from_int(a))

    sollya_lib_set_prec(precSo, None)

def pobyso_set_prec_so_so(newPrecSo):

    sollya_lib_set_prec(newPrecSo, None)

def pobyso_supnorm_so_so(polySo, funcSo, intervalSo, errorTypeSo = None,\

                         accuracySo = None):

    """

    Computes the supnorm of the approximation error between the given 

    polynomial and function.

    errorTypeSo defaults to "absolute".

    accuracySo defaults to 2^(-40).

    """

    if errorTypeSo is None:

        errorTypeSo = sollya_lib_absolute(None)

        errorTypeIsNone = True

    else:

        errorTypeIsNone = False

    #

    if accuracySo is None:

        # Notice the **!

        accuracySo = pobyso_constant_sa_so(RR(2**(-40)))

        accuracyIsNone = True

    else:

        accuracyIsNone = False

    pobyso_autoprint(accuracySo)

    resultSo = \

        sollya_lib_supnorm(polySo, funcSo, intervalSo, errorTypeSo, \

                              accuracySo)

    if errorTypeIsNone:

        sollya_lib_clear_obj(errorTypeSo)

    if accuracyIsNone:

        sollya_lib_clear_obj(accuracySo)

    return resultSo

# End pobyso_supnorm_so_so

def pobyso_taylor_expansion_no_change_var_so_so(functionSo, 

                                                degreeSo, 

                                                rangeSo,

                                                errorTypeSo=None, 

                                                sollyaPrecSo=None):

    """

    Compute the Taylor expansion without the variable change

    x -> x-intervalCenter.

    """

    # No global change of the working precision.

    if not sollyaPrecSo is None:

        initialPrecSo = sollya_lib_get_prec(None)

        sollya_lib_set_prec(sollyaPrecSo)

    # Error type stuff: default to absolute.

    if errorTypeSo is None:

        errorTypeIsNone = True

        errorTypeSo = sollya_lib_absolute(None)

    else:

        errorTypeIsNone = False

    intervalCenterSo = sollya_lib_mid(rangeSo, None)

    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo,

                                         intervalCenterSo,

                                         rangeSo, errorTypeSo, None)

    # taylorFormListSaSo is a Python list of Sollya objects references that 

    # are copies of the elements of taylorFormSo.

    # pobyso_get_list_elements_so_so clears taylorFormSo.

    (taylorFormListSaSo, numElementsSa, isEndEllipticSa) = \

        pobyso_get_list_elements_so_so(taylorFormSo)

    polySo = sollya_lib_copy_obj(taylorFormListSaSo[0])

    #print "Num elements:", numElementsSa

    sollya_lib_clear_obj(taylorFormSo)

    #polySo = taylorFormListSaSo[0]

    #errorRangeSo = sollya_lib_copy_obj(taylorFormListSaSo[2])

    errorRangeSo = taylorFormListSaSo[2]

    # No copy_obj needed here: a new object is created.

    maxErrorSo = sollya_lib_sup(errorRangeSo)