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

#! /opt/sage/sage

# @file cvp-run-02.sage

# Compute a polynomial from the CVP of the function vector.

# Return the polynomial, the error and the (one of) point(s) where this error

# is reached.

#

# @par Notes

# Change from version 01:

#   Take the approximation mode as an argument.

#   Instead of computing a single vector, compute several of them and pick

#   the one that gives the smallest error.

#   The points list is only used to complement the Chebyshev points lists

#   that are computed at each iteration.

#

scriptName = os.path.basename(__file__)

var('x')

contFracMaxErr = 10^-7

#

def initialize_env():

    """

    Load all necessary modules.

    """

    compiledSpyxDir = \

        "/home/storres/recherche/arithmetique/pobysoPythonSage/compiledSpyx"

    if compiledSpyxDir not in sys.path:

        sys.path.append(compiledSpyxDir)

    """

    As of Sage versions above or equal to 7.0, appending MPFI is not needed and

    crashes the program.

    if not 'mpfi' in sage.misc.cython.standard_libs:

        sage.misc.cython.standard_libs.append('mpfi')

    """

    load("/home/storres/recherche/arithmetique/pobysoPythonSage/src/sollya_lib.sage")

#    load("/home/storres/recherche/arithmetique/pobysoPythonSage/src/sageMpfr.spyx")

#    load("/home/storres/recherche/arithmetique/pobysoPythonSage/src/sageGMP.spyx")

    load("/home/storres/recherche/arithmetique/pobysoPythonSage/src/pobyso.py")

    load("/home/storres/recherche/arithmetique/pobysoPythonSage/src/sageSLZ/sageSLZ.sage")

    load("/home/storres/recherche/arithmetique/pobysoPythonSage/src/sageSLZ/sageNumericalOperations.sage")

    load("/home/storres/recherche/arithmetique/pobysoPythonSage/src/sageSLZ/sageRationalOperations.sage")

    # Matrix operations are loaded by polynomial operations.

    load("/home/storres/recherche/arithmetique/pobysoPythonSage/src/sageSLZ/sagePolynomialOperations.sage")

    load("/home/storres/recherche/arithmetique/pobysoPythonSage/src/sageSLZ/sageRunSLZ.sage")

    load("/home/storres/recherche/arithmetique/cvp/src/functions_for_cvp.sage")

# End initialize_env.

#

initialize_env()

from sageMpfr import *

from sageGMP  import *

import sys

#

def usage(scriptName):

    write = sys.stderr.write

    write("Usage:\n")

    write("  " + scriptName + " <func> <lowerBound> <upperBound> <mode> <coeffsExpList> \n")

    write(" " * (len(scriptName) + 3))

    write("<coeffsPrecList> <minCoeffsBoundExpList>\n")

    write(" " * (len(scriptName) + 3))

    write("<maxCoeffsBoundExpList> <pointList> <precision>\n")

    write("\nArguments:\n")

    write("  func\n")

    write("  lowerBound\n")

    write("  upperBound\n")

    write("  mode             one of \"abs\" or \"rel\"\n")

    write("  coeffsExpList\n")

    write("  coeffsPrecList\n")

    write("  minCoeffsBoundExpList\n")

    write("  maxCoeffsBoundExpList\n")

    write("  pointsList\n")

    write("  precision\n\n")

    ##

    sys.exit(1)

# End usage

#

argsCount = len(sys.argv)

#

## Arguments recovery

if argsCount < 11:

    usage(scriptName)

## Shell weird behavior may oblige to single quote arguments.

#  They should be removed before evaluation to avoid inadequate typing. 

funcAsString                    = re.sub("^'|'$", "", sys.argv[1])

lowerBoundAsString              = re.sub("^'|'$", "", sys.argv[2])

upperBoundAsString              = re.sub("^'|'$", "", sys.argv[3])

modeAsString                    = re.sub("^'|'$", "", sys.argv[4])

coeffsExpListAsString           = re.sub("^'|'$", "", sys.argv[5])

coeffsPrecListAsString          = re.sub("^'|'$", "", sys.argv[6])

minCoeffsBoundExpListAsString   = re.sub("^'|'$", "", sys.argv[7])

maxCoeffsBoundExpListAsString   = re.sub("^'|'$", "", sys.argv[8])

pointsListAsString              = re.sub("^'|'$", "", sys.argv[9])

precisionAsString               = re.sub("^'|'$", "", sys.argv[10])

#

## Arguments conversion

### Create first the rational field that will be used throughout the script.

precision               = int(precisionAsString)

realField               = RealField(precision)

#

func(x)                 = sage_eval(funcAsString, cmds="var('x')")

try:

    lowerBound          = realField(lowerBoundAsString)

except:

    lowerBoundAsRat     = QQ(lowerBoundAsString)

    lowerBound          = realField(lowerBoundAsRat)

try:

    upperBound          = realField(upperBoundAsString)

except:

    upperBoundAsRat     = QQ(upperBoundAsString)

    upperBound          = realField(upperBoundAsRat)

approxMode              = modeAsString.lower()

coeffsExpList           = sage_eval(coeffsExpListAsString)

coeffsPrecList          = sage_eval(coeffsPrecListAsString)

minCoeffsBoundExpList   = sage_eval(minCoeffsBoundExpListAsString)

maxCoeffsBoundExpList   = sage_eval(maxCoeffsBoundExpListAsString)

pointsList              = sage_eval(pointsListAsString)

#

if approxMode != "abs" and approxMode != "rel":

    sys.stderr.write("\n" +  scriptName + ": " + "\"" + approxMode + "\"")

    sys.stderr.write("is an invalide mode.\n")

    sys.stderr.write("Mode must be either \"abs\" or \"rel\". Aborting!\n")

    sys.exit(int(1))

maxExponent             = max(coeffsExpList)

## Debug printing.

sys.stderr.write("func: " + func._assume_str().replace("_SAGE_VAR_","",100) + \

                 "\n")

sys.stderr.write("precision: " + str(precision) + "\n")

sys.stderr.write("lowerBound: " + str(lowerBound) + "\n")

sys.stderr.write("upperBound: " + str(upperBound) + "\n")

sys.stderr.write("mode      :" + approxMode + "\n")

sys.stderr.write("coeffsExpList: " + str(coeffsExpList) + "\n")

sys.stderr.write("coeffsPrecList: " + str(coeffsPrecList) + "\n")

sys.stderr.write("minCoeffsBoundExpList: " + str(minCoeffsBoundExpList) + "\n")

sys.stderr.write("maxCoeffsBoundExpList: " + str(maxCoeffsBoundExpList) + "\n")

sys.stderr.write("pointsList: " + str(pointsList) + "\n")

sys.stderr.write("maxExponent: " + str(maxExponent) + "\n")

#

##

chebyshevPointsList = []

minMaxErrPoint  = oo

minMaxErr       = oo

for pointsNum in xrange(maxExponent+1, round(maxExponent * 8/2)+1):

    ## Compute the Chebyshev points and append the user given points list.

    chebyshevPointsList =  \

        cvp_chebyshev_maxis_for_interval(lowerBound, 

                                         upperBound, 

                                         pointsNum,

                                         realField, 

                                         contFracMaxErr)

    #print "cr-02 - Chebyshe points list computed."

    chebyshevPointsList = chebyshevPointsList + pointsList

    #

    ## Filter out those points that do not have a defined image.

    chebyshevPointsList = cvp_filter_out_undefined_image(chebyshevPointsList, func)

    #

    ##

    poly = cvp_cvp_polynomial(func, 

                              lowerBound, 

                              upperBound,

                              coeffsExpList,

                              coeffsPrecList, 

                              minCoeffsBoundExpList, 

                              maxCoeffsBoundExpList,

                              chebyshevPointsList,

                              realField)

    #sys.stderr.write(str(poly) + "\n")

    (maxisList,maxErr)  = cvp_polynomial_error_func_maxis(func,

                                                          poly, 

                                                          lowerBound, 

                                                          upperBound,

                                                          realField,

                                                          contFracMaxErr)

    ## In some cases

    #print "cr-02 - mcr-02 maxis list:", maxisList

    #print "cr-02 - maxErr:", maxErr

    #polyRing = PolynomialRing(realField, func.variables()[0])

    #errFunc(x) = func(x) - polyRing(poly)

    if approxMode == "abs":

        errFunc(x) = func(x) - poly(x)

    else:

        errFunc(x) = 1 - poly(x) / func(x)

    if len(maxisList) != 0:

        (maxErrPoint, maxErrPointErr) = \

            cvp_func_abs_max_for_points(errFunc,

                                        maxisList,

                                        realField)

        sys.stderr.write("# points: " + str(pointsNum))

        sys.stderr.write(" - MaxErr: " + str(maxErr.n(prec=53)))

        sys.stderr.write(" - MaxErrPoint: " + str(maxErrPoint) + "\n")

        if maxErr < minMaxErr:

            minMaxErr      = maxErr

            minMaxErrPoint = maxErrPoint

# End for

sys.stdout.write(str(minMaxErrPoint) + " " + str(minMaxErr.n(prec=53)) + "\n")

## Use int(0) as argument for sys.exit() in Sage because otherwise

# sys.exit(0) -> sys.exit(Integer(0)) (where Integer(0) is a Sage

# object, not the integer value 0 expected by Python) hence Python 

# returns 1!

sys.exit(int(0))