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)

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

(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 "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_autoprint(arg):

    sollya_lib_autoprint(arg,None)

def pobyso_cmp(rnArg, soCte):

    precisionOfCte = c_int(0)

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

    sollya_lib_get_prec_of_constant(precisionOfCte, soCte) 

    #print "Precision of constant: ", precisionOfCte

    RRRR = RealField(precisionOfCte.value)

    rnLocal = RRRR(0)

    sollya_lib_get_constant(get_rn_value(rnLocal), soCte)

    #print "rnDummy: ", rnDummy

    # Compare the local Sage RealNumber with rnArg.

    return(cmp_rn_value(rnArg, rnLocal))

def pobyso_constant(rnArg):

    return (sollya_lib_constant(get_rn_value(rnArg)))

def pobyso_constant_1():

    return(pobyso_constant_from_int(1))

def pobyso_constant_from_int(anInt):

    return(sollya_lib_constant_from_int(int(anInt)))

# Numeric Sollya function codes -> Sage mathematical function names

def pobyso_function_type_as_string(funcType):

    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(rnArg, soConst):

    set_rn_value(rnArg, soConst)

def pobyso_get_constant_as_rn(ctExp):

    precision  = pobyso_get_prec_of_constant(ctExp) 

    RRRR = RealField(precision)

    rn = RRRR(0)

    sollya_lib_get_constant(get_rn_value(rn), ctExp)

    return(rn)

def pobyso_get_constant_as_rn_with_rf(ctExp, realField):

    rn = realField(0)

    sollya_lib_get_constant(get_rn_value(rn), ctExp)

    return(rn)

def pobyso_get_free_variable_name():

    return(sollya_lib_get_free_variable_name())

def pobyso_get_function_arity(expression):

    arity = c_int(0)

    sollya_lib_get_function_arity(byref(arity),expression)

    return(int(arity.value))

def pobyso_get_head_function(expression):

    functionType = c_int(0)

    sollya_lib_get_head_function(byref(functionType), expression, None)

    return(int(functionType.value))

def pobyso_get_list_elements(soObj):

    # Type for array of pointers to sollya_obj_t

    listAddress = POINTER(c_longlong)()

    numElements = c_int(0)

    isEndElliptic = c_int(0)

    listAsList = []

    result = sollya_lib_get_list_elements(byref(listAddress),\

                                        byref(numElements),\

                                        byref(isEndElliptic),\

                                        soObj)

    if result == 0 :

        return None

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

        print "address ", i, " ->", listAddress[i]

        listAsList.append(listAddress[i])

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

# Get the maximum precision used for the numbers in a 

# Sollya expression.

# ToDo: 

# - error management;

# - correctly deal with numerical type such as DOUBLEEXTENDED.

def pobyso_get_max_prec_of_exp(soExp):

    maxPrecision = 0

    operator = pobyso_get_head_function(soExp)

    if (operator != SOLLYA_BASE_FUNC_CONSTANT) and \

    (operator != SOLLYA_BASE_FUNC_FREE_VARIABLE):

        (arity, subexpressions) = pobyso_get_subfunctions(soExp)

        for i in xrange(arity):

            maxPrecisionCandidate = \

            pobyso_get_max_prec_of_exp(subexpressions[i])

            if maxPrecisionCandidate > maxPrecision:

                maxPrecision = maxPrecisionCandidate

        return(maxPrecision)

    elif operator == SOLLYA_BASE_FUNC_CONSTANT:

        #print pobyso_get_prec_of_constant(soExp)

        return(pobyso_get_prec_of_constant(soExp))

    elif operator == SOLLYA_BASE_FUNC_FREE_VARIABLE:

        return(0)

    else:

        print "pobyso_get_max_prec_of_exp: unexepected operator."

        return(0)

def pobyso_get_sage_exp_from_sollya_exp(sollyaExp, realField = 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)

    operator = pobyso_get_head_function(sollyaExp)

    # Constants and the free variable are special cases.

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

    if (operator != SOLLYA_BASE_FUNC_CONSTANT) and \

       (operator != SOLLYA_BASE_FUNC_FREE_VARIABLE):

        (arity, subexpressions) = pobyso_get_subfunctions(sollyaExp)

        if arity == 1:

            sageExp = eval(pobyso_function_type_as_string(operator) + \

            "(" + pobyso_get_sage_exp_from_sollya_exp(subexpressions[0], realField)\

             + ")")

        elif arity == 2:

            if operator == SOLLYA_BASE_FUNC_POW:

                operatorAsString = "**"

            else:

                operatorAsString = pobyso_function_type_as_string(operator)

            sageExp = \

              eval("pobyso_get_sage_exp_from_sollya_exp(subexpressions[0], realField)"\

              + " " + operatorAsString + " " + \

                   "pobyso_get_sage_exp_from_sollya_exp(subexpressions[1], realField)")

        # We do not know yet how to deal with arity > 3 (is there any in Sollya anyway?).

        else:

            sageExp = eval('None')

        return(sageExp)

    elif operator == SOLLYA_BASE_FUNC_CONSTANT:

        #print "This is a constant"

        return pobyso_get_constant_as_rn_with_rf(sollyaExp, realField)

    elif operator == SOLLYA_BASE_FUNC_FREE_VARIABLE:

        #print "This is free variable"

        return(eval(sollya_lib_get_free_variable_name()))

    else:

        print "Unexpected"

        return eval('None')

# End pobyso_get_sage_poly_from_sollya_poly

def pobyso_get_subfunctions(expression):

    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(expression, 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(None,None)

    for i in xrange(arity.value):

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

        #print subs[i]

    return(int(arity.value), subs)

def pobyso_get_prec():

    retc = sollya_lib_get_prec(None)

    a = c_int(0)

    sollya_lib_get_constant_as_int(byref(a), retc)

    return(int(a.value))

def pobyso_get_prec_of_constant(ctExp):

    prec = c_int(0)

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

    return(int(prec.value))

def pobyso_lib_init():

    sollya_lib_init(None)

def pobyso_name_free_variable(freeVariableName):

    sollya_lib_name_free_variable(freeVariableName)

def pobyso_parse_string(string):

    return(sollya_lib_parse_string(string))

def pobyso_range(rnLowerBound, rnUpperBound):

    lowerBoundSo = sollya_lib_constant(get_rn_value(rnLowerBound))

    upperBoundSo = sollya_lib_constant(get_rn_value(rnUpperBound))

    rangeSo = sollya_lib_range(lowerBoundSo, upperBoundSo)

    return(rangeSo)

def pobyso_remez_canonical(function, \

                           degree, \

                           lowerBound, \

                           upperBound, \

                           weightSo = pobyso_constant_1(),

                           quality = None):

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

        functionSo = sollya_lib_parse_string(function)

    #    print "Is string!"

    elif sollya_lib_obj_is_function(function):

        functionSo = function

    #    print "Is Function!"

    degreeSo = pobyso_constant_from_int(degree)

    rangeSo = pobyso_range(lowerBound, upperBound)

    return(sollya_lib_remez(functionSo, degreeSo, rangeSo, quality, 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):

    a = c_int(p)

    precSo = c_void_p(sollya_lib_constant_from_int(a))

    sollya_lib_set_prec(precSo)

def pobyso_taylor(function, degree, point):

    return(sollya_lib_taylor(function, degree, point))

# TODO: take care of the interval and of point when it is an interval;

#       when errorType is not None;

#       take care of the other elements of the Taylor form (coefficients errors and

#       delta.

def pobyso_taylorform(function, degree, point = None, interval = None, errorType=None):

    """

    Compute the Taylor form of 'degre' for 'function' at 'point' 

    for 'interval' with 'errorType'. 

    point: must be a Real or a Real interval.

    return the Taylor form as an array

    """

    # Absolute as the default error.

    if errorType is None:

        errorTypeSo = sollya_lib_absolute()

    else:

        #TODO: deal with the other case.

        pass

    varSa = function.variables()[0]

    pointBaseRingString = str(point.base_ring())

    if not re.search('Real', pointBaseRingString):

        return None

    # Call Sollya but first "sollyafy" the arguments.

    sollya_lib_init(None)

    pobyso_name_free_variable(str(varSa))

    # Sollyafy the function.

    functionSo = pobyso_parse_string(function._assume_str())

    if sollya_lib_obj_is_error(functionSo):

        print "pobyso_tailorform: function string can't be parsed!"

        return None

    # Sollyafy the degree

    degreeSo = sollya_lib_constant_from_int(int(degree))

    # Sollyafy the point

    if not re.search('Interval', pointBaseRingString):

        pointSo  = pobyso_constant(point)

    else:

        # TODO: deal with the interval case.

        pass

    # Call Sollya

    taylorFormSo = sollya_lib_taylorform(functionSo, degreeSo, pointSo, errorTypeSo,\

                                         None)

    (tfsAsList, numElements, isEndElliptic) = pobyso_get_list_elements(taylorFormSo)

    polySo = tfsAsList[0]

    maxPrecision = pobyso_get_max_prec_of_exp(polySo)

    polyRealField = RealField(maxPrecision)

    expSa = pobyso_get_sage_exp_from_sollya_exp(polySo, polyRealField)

    sollya_lib_close()

    polynomialRing = polyRealField[str(varSa)]

    polySa = polynomial(expSa, polynomialRing)

    taylorFormSa = [polySa]

    return(taylorFormSa)                    

def pobyso_univar_polynomial_print_reverse(polySa):

    """

    Return the string representation of a univariate polynomial with

    monomial ordered in the x^0..x^n order of the monomials.

    Remember: Sage

    """

    polynomialRing = polySa.base_ring()

    # A very expensive solution:

    # -create a fake multivariate polynomial field with only one variable,

    #   specifying a negative lexicographical order;

    mpolynomialRing = PolynomialRing(polynomialRing.base(), \

                                     polynomialRing.variable_name(), \

                                     1, order='neglex')

    # - convert the univariate argument polynomial into a multivariate

    #   version;

    p = mpolynomialRing(polySa)

    # - return the string representation of the converted form.

    # There is no simple str() method defined for p's class.

    return(p.__str__())

#

print "Superficial test of pobyso:"    

print pobyso_get_prec()  

pobyso_set_prec(165)

print pobyso_get_prec()  

a=100

print type(a)

id(a)

print "Max arity: ", pobyso_max_arity

print "Function tripleDouble (43) as a string: ", pobyso_function_type_as_string(43)

print "Function None (44) as a string: ", pobyso_function_type_as_string(44)