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

"""

@file functions_for_cvp.sage

@package functions_for_cvp

Auxiliary functions used for CVP.

@author S.T.

"""

sys.stderr.write("Functions for CVP loading...\n")

#

def cvp_affine_from_chebyshev(lowerBound, upperBound):

    """

    Compute the affine parameters to map [-1, 1] to the interval given

    as argument.

    @param: lowerBound (of the target interval)

    @param: upperBound (of the target interval)

    @return the (scalingCoefficient, offset) tuple.

    """

    ## Check bounds consistency. Bounds must differ.

    if lowerBound >= upperBound:

        return None

    scalingCoefficient = (upperBound - lowerBound) / 2

    offset             = (lowerBound + upperBound) / 2

    return (scalingCoefficient, offset)

# End cvp_affine_from_chebyshev

#

def cvp_affine_to_chebyshev(lowerBound, upperBound):

    """

    Compute the affine parameters to map the interval given

    as argument to  [-1, 1].

    @param lowerBound (of the target interval)

    @param upperBound (of the target interval)

    @return the (scalingCoefficient, offset) tuple.

    """

    ## Check bounds consistency. Bounds must differ.

    if lowerBound >= upperBound:

        return None

    scalingCoefficient = 2 / (upperBound - lowerBound)

    ## If bounds are symmetrical with relation to 0, return 0 straight before

    #  attempting a division by 0.

    if lowerBound == -upperBound:

        offset         = 0

    else:

        offset         =  (lowerBound + upperBound) / (lowerBound - upperBound)

    return (scalingCoefficient, offset)

# End cvp_affine_to_chebyshev

#

def cvp_babai(readBasis, readBasisGso, vect):

    """

    Closest plane vector implementation as per Babaï.

    @param readBasis   : a lattice basis, preferably a reduced one; 

    @param readBasisGSO: the GSO of the previous basis;

    @param vect        : a vector, in coordinated in the ambient 

                        space of the lattice

    @return the closest vector to the input, in coordinates in the 

                        ambient space of the lattice.

    """

    ## A deep copy is not necessary here.

    curVect = copy(vect)

    #print "cvp_babai - Vector:", vect

    ## From index of last row down to 0.

    for vIndex in xrange(len(readBasis.rows())-1, -1, -1):

        curRBGSO = readBasisGso.row(vIndex)

        curVect = curVect - \

                  (round((curVect * curRBGSO)  / (curRBGSO * curRBGSO)) * \

                   readBasis.row(vIndex)) 

    return vect - curVect

# End cvp_babai

#

def cvp_bkz_reduce(intBasis, blockSize):

    """

    Thin and simplistic wrapper the HKZ function of the FP_LLL module.

    """

    fplllIntBasis = FP_LLL(intBasis)

    fplllIntBasis.BKZ(blockSize)

    return fplllIntBasis._sage_()

# End cvp_hkz_reduce.

#

def cvp_candidate_vectors(intVect, diffsList):

    """

    Computes a list of best approximation vectors based on the CVP

    computed vector using the differences between this vector and the 

    target function vector.

    @param intVect the vector from which candidates are built

    @param diffsList the list of the values that will be appended to the

           initial vector to form the candidates

    @return A matrix made of candidates rows or None if there is some invalid

            argument.

    @par How it works

    The idea is to create all the possible vectors based on the computed CVP

    vector and its differences to the target. Those are collected in 

    diffsList.@n

    From the vector we create all the possible combinations between the original

    value and the differences in diffsList. Some of those can be equal to 0.

    They are not taken into account. If nonZeroDiffsCount is the number of 

    diffsList element different from 0, there are 2^nonZeroDiffsCount 

    possible variants.@n

    To compute those, we build a 2^nonZeroDiffsCount rows matrix. For the

    first element in the vector corresponding to a non-zero in diffsList, we

    switch values in the matrix at every row.@n

    For the second such element (non-zero corresponding diffsList element), we

    switch values in the matrix at every second row for a string of two similar 

    values.@n

    For the third such element, change happens at every fourth row, for a 

    string of four similar values. And so on...

    @par Example

    @verbatim

    vect        =  [10, 20, 30, 40, 50]

    diffsList   =  [ 1,  0, -1,  0,  1]

    vectCandMat = [[10, 20, 30, 40, 50],

                   [11, 20, 30, 40, 50],

                   [10, 20, 29, 40, 50],

                   [11, 20, 29, 40, 50],

                   [10, 20, 30, 40, 51],

                   [11, 20, 30, 40, 51],

                   [10, 20, 29, 40, 51],

                   [11, 20, 29, 40, 51]]

    @endverbatim

"""

    ## Arguments check.

    vectLen = len(intVect)

    if vectLen != len(diffsList):

        return None

    #

    ## Compute the number of diffsList elements different from 0.

    nonZeroDiffsCount = 0

    for elem in diffsList:

        if elem != 0:

            nonZeroDiffsCount += 1

    if nonZeroDiffsCount == 0:

        return matrix(ZZ, 0, vectLen)

    numRows = 2^nonZeroDiffsCount

    vectMat = matrix(ZZ, numRows, vectLen)

    for rowIndex in xrange(0,numRows):

        effColIndex = 0

        for colIndex in xrange(0,vectLen):

            vectMat[rowIndex, colIndex] = intVect[colIndex]

            if diffsList[colIndex] != 0:

                if (rowIndex % 2^(effColIndex + 1)) >= 2^effColIndex:

                    vectMat[rowIndex, colIndex] += diffsList[colIndex]

                effColIndex += 1

    return vectMat

# End cvp_candidate_vectors

#

def cvp_coefficients_bounds_projection(executablePath, arguments):

    """

    Compute the min and max of the coefficients with linear

    programming projection.

    @param executablePath: the path to the binary program;

    @param arguments: a list of arguments to be givent to the binary

    @return: the min and max coefficients value arrays (in this order).

    """

    from subprocess import check_output

    commandLine = [executablePath] + arguments

    ## Get rid of output on stderr.

    devNull = open("/dev/null", "rw")

    ## Run ther program.

    otp = check_output(commandLine, stderr=devNull)

    devNull.close()

    ## Recover the results

    bounds = sage_eval(otp)

    minCoeffsExpList = []

    maxCoeffsExpList = [] 

    #print "bounds:", bounds

    for boundsPair in bounds:

        minCoeffsExpList.append(boundsPair[0])

        maxCoeffsExpList.append(boundsPair[1])

    print "minCoeffsExpList:", minCoeffsExpList

    print "maxCoeffsExpList:", maxCoeffsExpList

    return (minCoeffsExpList, maxCoeffsExpList)

# End cvp_coefficients_bounds_projection

def cvp_coefficients_bounds_projection_exps(executablePath, 

                                            arguments, 

                                            realField = None):

    """

    Compute the min and max exponents of the coefficients with linear

    programming projection.

    @param executablePath: the path to the binary program;

    @param arguments: a list of arguments to be givent to the binary

    @param realField: the realField to use for floating-point conversion.

    @return: the min and max exponents arrays (in this order).

    """

    ## If no realField is givne, fall back on RR.

    if realField is None:

        realField = RR

    (minCoeffsExpList, maxCoeffsExpList) = \

        cvp_coefficients_bounds_projection(executablePath,arguments)

    for index in xrange(0, len(minCoeffsExpList)):

        minCoeffsExpList[index] = \

            realField(minCoeffsExpList[index]).abs().log2().floor()

        maxCoeffsExpList[index] = \

            realField(maxCoeffsExpList[index]).abs().log2().floor()

    return (minCoeffsExpList, maxCoeffsExpList)

# End cvp_coefficients_bounds_projection_exps

def cvp_chebyshev_maxis_1k(numPoints, realField, contFracMaxErr = None):

    """

    Compute the Chebyshev maxis for some polynomial degree (numPoints, for the

    zeros, numPoints + 1 for the maxis) scaled to the [lowerBound, upperBound]

    interval.

    The list is returned as a row floating-point numbers is contFracMaxErr 

    is None.

    Otherwise elements are transformed into rational numbers. 

    """

    if numPoints < 1:

        return None

    if numPoints == 0:

        return [0]

    ## Check that realField has a "prec()" member.

    try:

        realField.prec()

    except:

        return None

    #

    maxisList = []

    ## n extrema are given by a degree (n-1) Chebyshev polynomial.

    #  formulas: cos ((i * pi) / (numPoints - 1)) for i = 0,...,numPoints-1.

    #  Source: https://en.wikipedia.org/wiki/Chebyshev_polynomials

    for index in xrange(0, numPoints):

        ## When number of points is odd (then, the Chebyshev degree is odd two), 

        #  the central point is 0. */

        if (numPoints % 2 == 1) and ((2 * index + 1) == numPoints):

            if contFracMaxErr is None:

                maxisList.append(realField(0))

            else:

                maxisList.append(0)

        else:

            currentCheby = \

                ((realField(index) * realField(pi)) / 

                 realField(numPoints-1)).cos()

            #print  "Current Cheby:", currentCheby

            ## Result is negated to have an increasing order list.

            if contFracMaxErr is None:

                maxisList.append(-currentCheby)

            else:

                maxisList.append(-currentCheby.nearby_rational(contFracMaxErr))

            #print "Relative error:", (currentCheby/maxisList[index])

    return maxisList

# End cvp_chebyshev_maxis_1k.

#

def cvp_chebyshev_maxis_for_interval(lowerBound, 

                                     upperBound, 

                                     numPoints,

                                     realField = None, 

                                     contFracMaxErr = None):

    """

    Compute the Chebyshev maxis for some polynomial degree (numPoints, for the

    maxis) scaled to the [lowerBound, upperBound] interval.

    If no contFracMaxErr argument is given, we return the list as "raw"

    floating-points. Otherwise, rational numbers are returned.

    """

    ## Arguments check.

    if lowerBound >= upperBound:

        return None

    if numPoints < 1:

        return None

    ## If no realField argument is given try to retrieve it from the 

    #  bounds. If impossible (e.g. rational bound), fall back on RR. 

    if realField is None:

        try:

            ### Force failure if parent does not have prec() member.

            lowerBound.parent().prec()

            ### If no exception is thrown, set realField.

            realField = lowerBound.parent()

        except:

            realField = RR

    #print "Real field:", realField

    ## We want "raw"floating-point nodes to only have one final rounding.

    chebyshevMaxisList = \

        cvp_chebyshev_maxis_1k(numPoints, realField)

    scalingFactor, offset = cvp_affine_from_chebyshev(lowerBound, upperBound)

    for index in xrange(0, len(chebyshevMaxisList)):

        chebyshevMaxisList[index] = \

            chebyshevMaxisList[index] * scalingFactor + offset

        if not contFracMaxErr is None:

            chebyshevMaxisList[index] = \

                chebyshevMaxisList[index].nearby_rational(contFracMaxErr)

    return chebyshevMaxisList                                                             

# End cvp_chebyshev_maxis_for_interval.

#

def cvp_chebyshev_zeros_1k(numPoints, realField, contFracMaxErr = None):

    """

    Compute the Chebyshev zeros for some polynomial degree (numPoints, for the

    zeros) scaled to the [lowerBound, upperBound] interval.

    The list is returned as row floating-point numbers is contFracMaxErr is None.

    Otherwise elements are transformed into rational numbers. 

    """

    if numPoints < 1:

        return None

    if numPoints == 0:

        return [0]

    ## Check that realField has a "prec()" member.

    try:

        realField.prec()

    except:

        return None

    #

    zerosList = []

    ## formulas: cos (((2*i - 1) * pi) / (2*numPoints)) for i = 1,...,numPoints.

    #  If i is -1-shifted, as in the following loop, formula is 

    #  cos (((2*i + 1) * pi) / (2*numPoints)) for i = 0,...,numPoints-1.

    for index in xrange(0, numPoints):

        ## When number of points is odd (then, the Chebyshev degree is odd two), 

        #  the central point is 0. */

        if (numPoints % 2 == 1) and ((2 * index + 1) == numPoints):

            if contFracMaxErr is None:

                zerosList.append(realField(0))

            else:

                zerosList.append(0)

        else:

            currentCheby = \

                ((realField(2*index+1) * realField(pi)) / 

                 realField(2*numPoints)).cos()

            #print  "Current Cheby:", currentCheby

            ## Result is negated to have an increasing order list.

            if contFracMaxErr is None:

                zerosList.append(-currentCheby)

            else:

                zerosList.append(-currentCheby.nearby_rational(contFracMaxErr))

            #print "Relative error:", (currentCheby/zerosList[index])

    return zerosList

# End cvp_chebyshev_zeros_1k.

#

def cvp_chebyshev_zeros_for_interval(lowerBound, 

                                     upperBound, 

                                     numPoints,

                                     realField = None, 

                                     contFracMaxErr = None):

    """

    Compute the Chebyshev zeros for some polynomial degree (numPoints, for the

    zeros) scaled to the [lowerBound, upperBound] interval.

    If no contFracMaxErr argument is given, we return the list as "raw"

    floating-points. Otherwise, rational numbers are returned.

    """

    ## Arguments check.

    if lowerBound >= upperBound:

        return None

    if numPoints < 1:

        return None

    ## If no realField argument is given try to retrieve it from the 

    #  bounds. If impossible (e.g. rational bound), fall back on RR. 

    if realField is None:

        try:

            ### Force failure if parent does not have prec() member.

            lowerBound.parent().prec()

            ### If no exception is thrown, set realField.

            realField = lowerBound.parent()

        except:

            realField = RR

    #print "Real field:", realField

    ## We want "raw"floating-point nodes to only have one final rounding.

    chebyshevZerosList = \

        cvp_chebyshev_zeros_1k(numPoints, realField)

    scalingFactor, offset = cvp_affine_from_chebyshev(lowerBound, upperBound)

    for index in xrange(0, len(chebyshevZerosList)):

        chebyshevZerosList[index] = \

            chebyshevZerosList[index] * scalingFactor + offset

        if not contFracMaxErr is None:

            chebyshevZerosList[index] = \

                chebyshevZerosList[index].nearby_rational(contFracMaxErr)

    return chebyshevZerosList                                                             

# End cvp_chebyshev_zeros_for_interval.

#

def cvp_common_scaling_exp(precision, 

                           precisionFraction = None):

    """

    Compute the common scaling factor (a power of 2).

    The exponent is a fraction of the precision;

    A extra factor is computed for the vector,

    exra factors are computed for the basis vectors, all in the corresponding

    functions.

    """

    ## Set a default value for the precisionFraction to 1/2.

    if precisionFraction is None:

        precisionFraction = 3/4

    """

    minBasisVectsNorm = oo

    currentNorm = 0

    for vect in floatBasis.rows():

        currentNorm = vect.norm()

        print "Current norm:", currentNorm

        if currentNorm < minBasisVectsNorm:

            minBasisVectsNorm = currentNorm

    currentNorm = funcVect.norm()

    print "Current norm:", currentNorm

    if currentNorm < minBasisVectsNorm:

        minBasisVectsNorm = currentNorm

    ## Check if a push is need because the smallest norm is below 0.    

    powerForMinNorm = floor(log(minBasisVectsNorm)/log2)

    print "Power for min norm :", powerForMinNorm

    print "Power for precision:", ceil(precision*precisionFraction)

    if powerForMinNorm < 0:

        return 2^(ceil(precision*precisionFraction) - powerFromMinNorm)

    else:

    """

    return ceil(precision * precisionFraction)

# End cvp_common_scaling_factor.

#

def cvp_cvp_polynomial(func, 

                       lowerBound, 

                       upperBound,

                       coeffsExpList,

                       coeffsPrecList, 

                       minCoeffsBoundExpList, 

                       maxCoeffsBoundExpList,

                       pointsList,

                       realField):

    """

    Compute a polynomial based on a CVP vector for a function approximation. 

    """

    ## Build the basis for CVP.

    ### Common scaling exponent.

    commonScalingExp = cvp_common_scaling_exp(realField.prec())

    ### Monomials list derived from exponenets list.

    monomialsList = cvp_monomials_list(coeffsExpList)

    ### Float basis first

    floatBasis = cvp_float_basis(monomialsList, pointsList, realField)

    ### Integral basis.

    extraScalingExpsList = \

        cvp_extra_basis_scaling_exps(coeffsPrecList, minCoeffsBoundExpList)

    intBasis = \

        cvp_int_basis(floatBasis, commonScalingExp, extraScalingExpsList)

    ### The reduced basis.

    redIntBasis = cvp_lll_reduce(intBasis)

    ### The reduced basis GSO.

    (redIntBasisGso, mu) = redIntBasis.gram_schmidt()

    ## Function vector.

    #### Floating-point vector.

    floatFuncVect = \

        cvp_float_vector_for_approx(func, pointsList, realField)

    extraFuncScalingExp = \

        cvp_extra_function_scaling_exp(floatBasis, 

                                       floatFuncVect, 

                                       maxCoeffsBoundExpList)

    #### Integral vector.

    intFuncVect   = \

        cvp_int_vector_for_approx(floatFuncVect, 

                                  commonScalingExp, 

                                  extraFuncScalingExp)

    ## CPV vector

    ### In ambient basis.

    cvpVectInAmbient = cvp_babai(redIntBasis, redIntBasisGso, intFuncVect)

    ### In intBasis.

    cvpVectInIntBasis = cvp_vector_in_basis(cvpVectInAmbient, intBasis)

    ## Polynomial coefficients

    cvpPolynomialCoeffs = \

        cvp_polynomial_coeffs_from_vect(cvpVectInIntBasis, 

                                        coeffsPrecList, 

                                        minCoeffsBoundExpList, 

                                        extraFuncScalingExp)

    #print "CVP polynomial coefficients:"

    #print cvpPolynomialCoeffs

    ## Polynomial in Sage.

    cvpPolySa =  \

        cvp_polynomial_from_coeffs_and_exps(cvpPolynomialCoeffs,

                                            coeffsExpList,

                                            realField[x])

    #print cvpPolySa

    return cvpPolySa

# End cvp_cvp_polynomial.

#

def cvp_evaluation_points_list(funct, 

                               maxDegree, 

                               lowerBound, 

                               upperBound, 

                               realField = None):

    """

    Compute a list of points for the vector creation.

    Strategy:

    - compute the zeros of the function-remez;

    - compute the zeros of the function-rounded_remez;

    - compute the Chebyshev points.

    Merge the whole thing.

    """

# End cvp_evaluation_points_list.

#

def cvp_extra_basis_scaling_exps(precsList, minExponentsList):

    """

    Computes the exponents for the extra scaling down of the elements of the

    basis.

    This extra scaling down is necessary when there are elements of the 

    basis that have negative exponents.

    """

    extraScalingExpsList = []

    for minExp, prec in zip(minExponentsList, precsList):

        if minExp - prec < 0:

            extraScalingExpsList.append(minExp - prec)

        else:

            extraScalingExpsList.append(0)

    return extraScalingExpsList

# End cvp_extra_basis_scaling_exps.

#

def cvp_extra_function_scaling_exp(floatBasis,

                                   floatFuncVect,

                                   maxExponentsList):

    """

    Compute the extra scaling exponent for the function vector in ordre to 

    guarantee that the maximum exponent can be reached for each element of the

    basis.

    """

    ## Check arguments

    if floatBasis.nrows() == 0 or \

       floatBasis.ncols() == 0 or \

       len(floatFuncVect)      == 0:

        return 1

    if len(maxExponentsList) != floatBasis.nrows():

        return None

    ## Compute the log of the norm of the function vector.

    funcVectNormLog  = log(floatFuncVect.norm())

    ## Compute the extra scaling factor for each vector of the basis.

    #  for norms ratios an maxExponentsList. 

    extraScalingExp = 0

    for (row, maxExp) in zip(floatBasis.rows(),maxExponentsList):

        rowNormLog     = log(row.norm())

        normsRatioLog2 = floor((funcVectNormLog - rowNormLog) / log2) 

        #print "Current norms ratio log2:", normsRatioLog2

        scalingExpCandidate = normsRatioLog2 - maxExp 

        #print "Current scaling exponent candidate:", scalingExpCandidate

        if scalingExpCandidate < 0:

            if -scalingExpCandidate > extraScalingExp:

                extraScalingExp = -scalingExpCandidate

    return extraScalingExp

# End cvp_extra_function_scaling_exp.

#

def cvp_filter_out_undefined_image(initialPointList, func):

    """

    Fom the initial list and function, built a list of points which have a

    defined image.

    """

    finalPointsList = []

    for iterPoint in initialPointList:

        try:

            finalPointsList.append(func(iterPoint))

            #print "cfoui - Point", iterPoint, "done."

        except ValueError:

            pass

            #print "cfoui - Point", iterPoint, "failed."

    return finalPointsList

# End    cvp_filter_out_undefined_image

def cvp_float_basis(monomialsList, pointsList, realField):

    """

    For a given monomials list and points list, compute the floating-point

    basis matrix.

    """

    numRows    = len(monomialsList)

    numCols    = len(pointsList)

    if numRows == 0 or numCols == 0:

        return matrix(realField, 0, 0)

    #

    floatBasis = matrix(realField, numRows, numCols)

    for rIndex in xrange(0, numRows):

        for cIndex in xrange(0, numCols):

            floatBasis[rIndex,cIndex] = \

                monomialsList[rIndex](realField(pointsList[cIndex]))

    return floatBasis

# End cvp_float_basis

#

def cvp_float_vector_for_approx(func, 

                                basisPointsList,

                                realField):

    """

    Compute a floating-point vector for the function and the points list.

    """

    #

    ## Some local variables.

    basisPointsNum = len(basisPointsList)

    #

    floatVector = vector(realField, basisPointsNum)

    ##

    for vIndex in xrange(0,basisPointsNum):

        ### We assume the all for all points their image is defined.

        floatVector[vIndex] = \

                func(basisPointsList[vIndex])

    return floatVector

# End cvp_float_vector_for_approx.

#

def cvp_func_abs_max_for_points(func, pointsList, realField):

    """

    Compute the maximum absolute value of a function for a list of 

    points.

    If several points yield the same value and this value is the maximum,

    the last visited one is returned.

    @return (theMaxPoint, theMaxValue) 

    """

    evalDict = dict()

    #

    ## An empty points list is really an error. We cannot return an empty tuple.

    if len(pointsList) == 0:

        return None

    #

    for pt in pointsList:

        try:

            evalDict[abs(realField(func(pt)))] = realField(pt)

            #print "cfamfp - point:", realField(pt) , "- image:", abs(realField(func(pt)))

        except ValueError:

            pass

    maxAbs = max(evalDict.keys())

    #print "cfamfp - maxAbs:", maxAbs, evalDict[maxAbs]

    return (evalDict[maxAbs], maxAbs)

# End cvp_func_abs_max_for_points

#

def cvp_hkz_reduce(intBasis):

    """

    Thin and simplistic wrapper on the HKZ function of the FP_LLL module.

    """

    fplllIntBasis = FP_LLL(intBasis)

    fplllIntBasis.HKZ()

    return fplllIntBasis._sage_()

# End cvp_hkz_reduce.

#

def cvp_int_basis(floatBasis, 

                  commonScalingExp, 

                  extraScalingExpsList):

    """

    From the float basis, the common scaling factor and the extra exponents 

    compute the integral basis.

    """

    ## Checking arguments.

    if floatBasis.nrows() != len(extraScalingExpsList):

        return None

    intBasis = matrix(ZZ, floatBasis.nrows(), floatBasis.ncols())

    for rIndex in xrange(0, floatBasis.nrows()):

        for cIndex in xrange(0, floatBasis.ncols()):

            intBasis[rIndex, cIndex] = \

            (floatBasis[rIndex, cIndex] * \

            2^(commonScalingExp + extraScalingExpsList[rIndex])).round()

    return intBasis

# End cvp_int_basis.

#

def cvp_int_vector_for_approx(floatVect, commonScalingExp, extraScalingExp):

    """

    Compute the integral version of the function vector.

    """

    totalScalingFactor = 2^(commonScalingExp + extraScalingExp)

    intVect = vector(ZZ, len(floatVect))

    for index in xrange(0, len(floatVect)):

        intVect[index] = (floatVect[index] * totalScalingFactor).round()

    return intVect

# End cvp_int_vector_for_approx.

#

def cvp_lll_reduce(intBasis):

    """

    Thin and simplistic wrapper around the LLL function

    """

    return intBasis.LLL()

# End cvp_lll_reduce.

#

def cvp_monomials_list(exponentsList, varName = None):

    """

    Create a list of monomials corresponding to the given exponentsList.

    Monomials are defined as functions.

    """

    monomialsList = []

    for exponent in exponentsList:

        if varName is None:

            monomialsList.append((x^exponent).function(x))

        else:

            monomialsList.append((varName^exponent).function(varName))

    return monomialsList

# End cvp_monomials_list.

#

def cvp_polynomial_coeffs_from_vect(cvpVectInBasis,

                                    coeffsPrecList, 

                                    minCoeffsExpList,

                                    extraFuncScalingExp):

    """

    Computes polynomial coefficients out of the elements of the CVP vector.

    @todo

    Rewrite the code with a single exponentiation, at the very end.

    """

    numElements = len(cvpVectInBasis)

    ## Arguments check.

    if numElements != len(minCoeffsExpList) or \

       numElements != len(coeffsPrecList):

        return None

    polynomialCoeffsList = []

    for index in xrange(0, numElements):

        currentCoeff = cvp_round_to_bits(cvpVectInBasis[index],

                                         coeffsPrecList[index])

        #print "cvp_polynomial_coeffs - current coefficient rounded:", currentCoeff

        currentCoeff *= 2^(minCoeffsExpList[index])

        #print "cvp_polynomial_coeffs - current coefficient scaled :", currentCoeff

        currentCoeff *= 2^(-coeffsPrecList[index])

        #print "cvp_polynomial_coeffs - current coefficient prec corrected :", currentCoeff

        currentCoeff *= 2^(-extraFuncScalingExp)

        #print "cvp_polynomial_coeffs - current coefficient extra func scaled :", currentCoeff

        polynomialCoeffsList.append(currentCoeff)

    return polynomialCoeffsList

# En polynomial_coeffs_from_vect.

#

def cvp_polynomial_error_func_maxis(funcSa,

                                    polySa, 

                                    lowerBoundSa, 

                                    upperBoundSa,

                                    realField = None,

                                    contFracMaxErr = None):

    """

    For a given function, approximation polynomial and interval (given 

    as bounds) compute the maxima of the functSa - polySo on the interval.

    Also computes the infnorm of the error function.

    """

    ## Arguments check.

    # @todo.

    ## If no realField argument is given try to retrieve it from the 

    #  bounds. If impossible (e.g. rational bound), fall back on RR. 

    if realField is None:

        try:

            ### Force failure if parent does not have prec() member.

            lowerBound.parent().prec()

            ### If no exception is thrown, set realField.

            realField = lowerBound.parent()

        except:

            realField = RR

    #print "Real field:", realField

    ## Compute the Sollya version of the function.

    funcAsStringSa = funcSa._assume_str().replace("_SAGE_VAR_","",100)

    #print "cpefa:", funcAsStringSa

    funcSo  = pobyso_parse_string(funcAsStringSa)

    if pobyso_is_error_so_sa(funcSo):

        pobyso_clear_obj(funcSo)

        return None

    #print "cpefa: function conversion OK."

    ## Compute the Sollya version of the polynomial.

    ## The conversion is made from a floating-point coefficients polynomial.

    try:

        polySa.base_ring().prec()

        convFloatPolySa = polySa

    except:

        convFloatPolySa = realField[polySa.variables()[0]](polySa)

    polySo = pobyso_float_poly_sa_so(convFloatPolySa)

    if pobyso_is_error_so_sa(funcSo):

        pobyso_clear_obj(funcSo)

        pobyso_clear_obj(polySo)

        return None

    #print "cpefa: polynomial conversion to Sollya OK."

    ## Copy both funcSo and polySo as they are needed later for the infnorm..

    errorFuncSo = \

        sollya_lib_build_function_sub(sollya_lib_copy_obj(funcSo), 

                                      sollya_lib_copy_obj(polySo))

    if pobyso_is_error_so_sa(errorFuncSo):

        pobyso_clear_obj(errorFuncSo)

        return None

    #print "cpefa: error function in Sollya OK."

    ## Evaluate the error function at both bounds, it may be usefull, if

    #  the error function is monotonous

    ## Lower bound stuff.

    lowerBoundSo = pobyso_constant_sa_so(lowerBoundSa)

    if pobyso_is_error_so_sa(lowerBoundSo):

        pobyso_clear_obj(errorFuncSo)

        pobyso_clear_obj(lowerBoundSo)

        return None

    errFuncAtLbSa = pobyso_evaluate_so_sa(errorFuncSo, lowerBoundSo)

    pobyso_clear_obj(lowerBoundSo)

    if errFuncAtLbSa is None:

        pobyso_clear_obj(errorFuncSo)

        return None

    ## The result of the evaluation can be an interval.

    try:

        errFuncAtLbSa = errFuncAtLbSa.center()

    except:

        errFuncAtLbSa = errFuncAtLbSa

    #print "cpefa: errFuncAtLbSa:", errFuncAtLbSa

    ## Upper bound stuff.

    upperBoundSo = pobyso_constant_sa_so(upperBoundSa)

    if pobyso_is_error_so_sa(upperBoundSo):

        pobyso_clear_obj(errorFuncSo)

        pobyso_clear_obj(upperBoundSo)

        return None

    errFuncAtUbSa = pobyso_evaluate_so_sa(errorFuncSo, upperBoundSo)

    pobyso_clear_obj(upperBoundSo)

    if errFuncAtUbSa is None:

        return None

    ## The result of the evaluation can be an interval.

    try:

        errFuncAtUbSa = errFuncAtUbSa.center()

    except:

        errFuncAtUbSa = errFuncAtUbSa

    #print "cpefa: errFuncAtUbSa:", errFuncAtUbSa

    ## Compute the derivative.

    diffErrorFuncSo = pobyso_diff_so_so(errorFuncSo)

    pobyso_clear_obj(errorFuncSo)

    if pobyso_is_error_so_sa(diffErrorFuncSo):

        pobyso_clear_obj(diffErrorFuncSo)

        return None

        #print "cpefa: polynomial conversion to Sollya OK."

    #print "cpefa: error function derivative in Sollya OK."

    ## Compute the interval.

    intervalSo = pobyso_range_from_bounds_sa_so(lowerBound, upperBound)

    if pobyso_is_error_so_sa(intervalSo):

        pobyso_clear_obj(diffErrorFuncSo)

        pobyso_clear_obj(intervalSo)

        return None

    ## Compute the infnorm of the error function.

    errFuncSupNormSa = pobyso_supnorm_so_sa(polySo, 

                                            funcSo, 

                                            intervalSo, 

                                            realFieldSa = realField)

    pobyso_clear_obj(polySo)

    pobyso_clear_obj(funcSo)

    ## Compute the zeros of the derivative.

    ### Give it a first try with dirtyfindzeros but it may

    #   fail.

    errorFuncMaxisSo = pobyso_dirty_find_zeros_so_so(diffErrorFuncSo, 

                                                     intervalSo)

    if pobyso_is_error_so_sa(errorFuncMaxisSo):

        pobyso_clear_obj(diffErrorFuncSo)

        pobyso_clear_obj(intervalSo)

        pobyso_clear_obj(errorFuncMaxisSo)

        return None

    #print "cpefa: error function maxis in Sollya OK."

    errorFuncMaxisSa = pobyso_float_list_so_sa(errorFuncMaxisSo)

    pobyso_clear_obj(errorFuncMaxisSo)

    #### If the list is empty try the more sophisticated

    #    findzeros.

    if len(errorFuncMaxisSa) == 0:

        errorFuncMaxisSo = \

            pobyso_find_zeros_so_so(diffErrorFuncSo, 

                                    intervalSo)

        if pobyso_is_error_so_sa(errorFuncMaxisSo):

            pobyso_clear_obj(diffErrorFuncSo)

            pobyso_clear_obj(intervalSo)

            pobyso_clear_obj(errorFuncMaxisSo)

            return None

        #print "cpefa: error function maxis in Sollya OK."

        #print "cpefa:", ; pobyso_autoprint(errorFuncMaxisSo)

        ##### The findzeros functions returns intervals (ranges).

        errorFuncRangeMaxisSa = pobyso_range_list_so_sa(errorFuncMaxisSo)

        #print "cpefa:", errorFuncRangeMaxisSa

        pobyso_clear_obj(errorFuncMaxisSo)

        errorFuncMaxisSa = []

        if len(errorFuncRangeMaxisSa) != 0:

            for interval in errorFuncRangeMaxisSa:

                errorFuncMaxisSa.append(interval.center())

            #print "cpefa:", errorFuncMaxisSa

        else: # The error function is monotonous. It reaches its maximum

              # at either of the bounds.

              errFuncMaxAtBounds = max(abs(errFuncAtLbSa), abs(errFuncAtUbSa))

              if errFuncMaxAtBounds == abs(errFuncAtLbSa):

                  errorFuncMaxisSa.append(lowerBoundSa)

              else:

                  errorFuncMaxisSa.append(upperBoundSa)

    pobyso_clear_obj(diffErrorFuncSo)

    pobyso_clear_obj(intervalSo)

    ## If required, convert the numbers to rational numbers.

    #print "cpefa - errorFuncMaxisSa:", errorFuncMaxisSa

    if not contFracMaxErr is None:

        for index in xrange(0, len(errorFuncMaxisSa)):

            errorFuncMaxisSa[index] = \

                errorFuncMaxisSa[index].nearby_rational(contFracMaxErr)

    return (errorFuncMaxisSa, errFuncSupNormSa)

# End cvp_polynomial_error_func_maxis

#

def cvp_polynomial_from_coeffs_and_exps(coeffsList, 

                                        expsList, 

                                        polyField = None,

                                        theVar = None):

    """

    Build a polynomial in the classical monomials (possibly lacunary) basis 

    from a list of coefficients and a list of exponents. The polynomial is in

    the polynomial field given as argument.

    """

    if len(coeffsList) != len(expsList):

        return None

    ## If no variable is given, "x" is used.

    ## If no polyField is given, we fall back on a polynomial field on RR.

    if theVar is None:

        if polyField is None:

            theVar = x

            polyField = RR[theVar]

        else:

            theVar = var(polyField.variable_name())

    else: ### theVar is set.

        ### No polyField, fallback on RR[theVar].

        if polyField is None:

            polyField = RR[theVar]

        else: ### Both the polyFiled and theVar are set: create a new polyField

              #   with theVar. The original polyField is not affected by the 

              #   change.

            polyField = polyField.change_var(theVar)  

    ## Seed the polynomial.

    poly = polyField(0)

    ## Append the terms.

    for coeff, exp in zip(coeffsList, expsList):

        poly += polyField(coeff * theVar^exp)

    return poly

# End cvp_polynomial_from_coeffs_and_exps.

#

def cvp_remez_all_poly_error_func_maxis(funct, 

                                        maxDegree, 

                                        lowerBound, 

                                        upperBound,

                                        precsList, 

                                        realField = None,

                                        contFracMaxErr = None):

    """

    For a given function f, a degree and an interval, compute the 

    maxima of the (f-remez_d(f)) function and those of the (f-truncRemez_d(f)) 

    function over the interval.

    """

    ## If no realField argument is given try to retrieve it from the 

    #  bounds. If impossible (e.g. rational bound), fall back on RR. 

    if realField is None:

        try:

            ### Force failure if parent does not have prec() member.

            lowerBound.parent().prec()

            ### If no exception is thrown, set realField.

            realField = lowerBound.parent()

        except:

            realField = RR

    #print "Real field:", realField

    funcSa         = func

    funcAsStringSa = func._assume_str().replace("_SAGE_VAR_","",100)

    print "Function as a string:", funcAsStringSa

    ## Create the Sollya version of the function and compute the

    #  Remez approximant

    funcSo  = pobyso_parse_string(funcAsStringSa)

    pStarSo = pobyso_remez_canonical_sa_so(funcAsStringSa, 

                                           maxDegree, 

                                           lowerBound,

                                           upperBound)

    ## Compute the error function and its derivative

    ### First create the error function with copies since they are needed later.

    errorFuncSo = sollya_lib_build_function_sub(sollya_lib_copy_obj(pStarSo), 

                                                sollya_lib_copy_obj(funcSo))

    intervalSo = pobyso_range_from_bounds_sa_so(lowerBound, upperBound)

    diffErrorFuncSo = pobyso_diff_so_so(errorFuncSo)

    pobyso_clear_obj(errorFuncSo)

    ## Compute the zeros of the derivative.

    errorFuncMaxisSo = pobyso_dirty_find_zeros_so_so(diffErrorFuncSo, intervalSo)

    pobyso_clear_obj(diffErrorFuncSo)

    errorFuncMaxisSa = pobyso_float_list_so_sa(errorFuncMaxisSo)

    pobyso_clear_obj(errorFuncMaxisSo)

    ## Compute the truncated Remez polynomial and the error function.

    truncFormatListSo = pobyso_build_list_of_ints_sa_so(*precsList)

    pTruncSo = pobyso_round_coefficients_so_so(pStarSo, truncFormatListSo)

    pobyso_clear_obj(pStarSo)

    ### Compute the error function. This time we consume both the function

    #   and the polynomial.

    errorFuncSo = sollya_lib_build_function_sub(pTruncSo, funcSo)

    ## Compute the derivative.

    diffErrorFuncSo = pobyso_diff_so_so(errorFuncSo)

    pobyso_clear_obj(errorFuncSo)

    ## Compute the zeros of the derivative.

    errorFuncMaxisSo = pobyso_dirty_find_zeros_so_so(diffErrorFuncSo, intervalSo)

    pobyso_clear_obj(diffErrorFuncSo)

    pobyso_clear_obj(intervalSo)

    errorFuncTruncMaxisSa = pobyso_float_list_so_sa(errorFuncMaxisSo)

    pobyso_clear_obj(errorFuncMaxisSo)

    ## Merge with the first list, removing duplicates if any.

    errorFuncMaxisSa.extend(elem for elem in errorFuncTruncMaxisSa \

                            if elem not in errorFuncMaxisSa)

    errorFuncMaxisSa.sort()

    ## If required, convert the numbers to rational numbers.

    if not contFracMaxErr is None:

        for index in xrange(0, len(errorFuncMaxisSa)):

            errorFuncMaxisSa[index] = \

                errorFuncMaxisSa[index].nearby_rational(contFracMaxErr)

    return errorFuncMaxisSa

# End cvp_remez_all_poly_error_func_maxis.

#

def cvp_remez_all_poly_error_func_zeros(funct, 

                                        maxDegree, 

                                        lowerBound, 

                                        upperBound,

                                        precsList, 

                                        realField = None,

                                        contFracMaxErr = None):

    """

    For a given function f, a degree and an interval, compute the 

    zeros of the (f-remez_d(f)) function and those of the (f-truncRemez_d(f)) 

    function over the interval.

    """

    ## If no realField argument is given try to retrieve it from the 

    #  bounds. If impossible (e.g. rational bound), fall back on RR. 

    if realField is None:

        try:

            ### Force failure if parent does not have prec() member.

            lowerBound.parent().prec()

            ### If no exception is thrown, set realField.

            realField = lowerBound.parent()

        except:

            realField = RR

    #print "Real field:", realField

    funcSa         = func

    funcAsStringSa = func._assume_str().replace("_SAGE_VAR_","",100)

    print "Function as a string:", funcAsStringSa

    ## Create the Sollya version of the function and compute the

    #  Remez approximant

    funcSo  = pobyso_parse_string(funcAsStringSa)

    pStarSo = pobyso_remez_canonical_sa_so(funcAsStringSa, 

                                           maxDegree, 

                                           lowerBound,

                                           upperBound)

    ## Compute the zeroes of the error function.

    ### First create the error function with copies since they are needed later.

    errorFuncSo = sollya_lib_build_function_sub(sollya_lib_copy_obj(pStarSo), 

                                                sollya_lib_copy_obj(funcSo))

    intervalSo = pobyso_range_from_bounds_sa_so(lowerBound, upperBound)

    errorFuncZerosSo = pobyso_dirty_find_zeros_so_so(errorFuncSo, intervalSo)

    pobyso_clear_obj(errorFuncSo) # Clears funcSo and pStarSo as well.

    #print "Zeroes of the error function from Sollya:"

    #pobyso_autoprint(errorFuncZerosSo)

    errorFuncZerosSa = pobyso_float_list_so_sa(errorFuncZerosSo)

    pobyso_clear_obj(errorFuncZerosSo)

    #print "\nZeroes of the error function from Sage:"

    #print errorFuncZerosSa

    #

    ## Deal with the truncated polynomial now.

    ### Create the formats list. Notice the "*" before the list variable name.

    truncFormatListSo = pobyso_build_list_of_ints_sa_so(*precsList)

    #print "Precisions list as Sollya object:", 

    #pobyso_autoprint(truncFormatListSo)

    pTruncSo = pobyso_round_coefficients_so_so(pStarSo, truncFormatListSo)

    #print "Truncated polynomial:", ; pobyso_autoprint(pTruncSo)

    ### Compute the error function. This time we consume both the function

    #   and the polynomial.

    errorFuncSo = sollya_lib_build_function_sub(pTruncSo, 

                                                funcSo)

    errorFuncZerosSo = pobyso_dirty_find_zeros_so_so(errorFuncSo, intervalSo)

    pobyso_clear_obj(pStarSo)

    #print "Zeroes of the error function for the truncated polynomial from Sollya:"

    #pobyso_autoprint(errorFuncZerosSo)

    errorFuncTruncZerosSa = pobyso_float_list_so_sa(errorFuncZerosSo)

    pobyso_clear_obj(errorFuncZerosSo)

    ## Sollya objects clean up.

    pobyso_clear_obj(intervalSo)

    errorFuncZerosSa += errorFuncTruncZerosSa

    errorFuncZerosSa.sort()

    ## If required, convert the numbers to rational numbers.

    if not contFracMaxErr is None:

        for index in xrange(0, len(errorFuncZerosSa)):

            errorFuncZerosSa[index] = \

                errorFuncZerosSa[index].nearby_rational(contFracMaxErr)

    return errorFuncZerosSa

    #

# End remez_all_poly_error_func_zeros

#

def cvp_remez_poly_error_func_zeros(funct, 

                                    maxDegree, 

                                    lowerBound, 

                                    upperBound,

                                    realField = None,

                                    contFracMaxErr = None):

    """

    For a given function f, a degree and an interval, compute the 

    zeros of the (f-remez_d(f)) function.

    """

    ## If no realField argument is given try to retrieve it from the 

    #  bounds. If impossible (e.g. rational bound), fall back on RR. 

    if realField is None:

        try:

            ### Force failure if parent does not have prec() member.

            lowerBound.parent().prec()

            ### If no exception is thrown, set realField.

            realField = lowerBound.parent()

        except:

            realField = RR

    #print "Real field:", realField

    funcSa         = func

    funcAsStringSa = func._assume_str().replace("_SAGE_VAR_","",100)

    #print "Function as a string:", funcAsStringSa

    ## Create the Sollya version of the function and compute the

    #  Remez approximant

    funcSo  = pobyso_parse_string(funcAsStringSa)

    pStarSo = pobyso_remez_canonical_sa_so(funcAsStringSa, 

                                           maxDegree, 

                                           lowerBound,

                                           upperBound)

    ## Compute the zeroes of the error function.

    ### Error function creation consumes both pStarSo and funcSo.

    errorFuncSo = sollya_lib_build_function_sub(pStarSo, funcSo)

    intervalSo = pobyso_range_from_bounds_sa_so(lowerBound, upperBound)

    errorFuncZerosSo = pobyso_dirty_find_zeros_so_so(errorFuncSo, intervalSo)

    pobyso_clear_obj(errorFuncSo) # Clears funcSo and pStarSo as well.

    #print "Zeroes of the error function from Sollya:"

    #pobyso_autoprint(errorFuncZerosSo)

    errorFuncZerosSa = pobyso_float_list_so_sa(errorFuncZerosSo)

    pobyso_clear_obj(errorFuncZerosSo)

    ## Sollya objects clean up.

    pobyso_clear_obj(intervalSo)

    ## Converting to rational numbers, if required.

    if not contFracMaxErr is None:

        for index in xrange(0, len(errorFuncZerosSa)):

            errorFuncZerosSa[index] = \

                errorFuncZerosSa[index].nearby_rational(contFracMaxErr)

    return errorFuncZerosSa

    #

# End remez_poly_error_func_zeros

#

def cvp_round_to_bits(num, numBits):

    """

    Round "num" to "numBits" bits.

    @param num    : the number to round;

    @param numBits: the number of bits to round to.

    """

    if num == 0:

        return num

    log2ofNum  = 0

    numRounded = 0

    ## Avoid conversion to floating-point if not necessary.

    try:

        log2OfNum = num.abs().log2()

    except:

        log2OfNum = RR(num).abs().log2()

    ## Works equally well for num >= 1 and num < 1.

    log2OfNum = ceil(log2OfNum)

    # print "cvp_round_to_bits - Log2OfNum:", log2OfNum

    divMult = log2OfNum - numBits

    #print "cvp_round_to_bits - DivMult:", divMult

    if divMult != 0:

        numRounded = round(num/2^divMult) * 2^divMult

    else:

        numRounded = num

    return numRounded

# End cvp_round_to_bits

#

def cvp_vector_in_basis(vect, basis):

    """

    Compute the coordinates of "vect" in "basis" by

    solving a linear system.

    @param vect:   the vector we want to compute the coordinates of

                   in coordinates of the ambient space;

    @param basis:  the basis we want to compute the coordinates in

                   as a matrix relative to the ambient space.

    """

    ## Create the variables for the linear equations.

    varDeclString = ""

    basisIterationsRange = xrange(0, basis.nrows())

    ### Building variables declaration string.

    for index in basisIterationsRange:

        varName = "a" + str(index)

        if varDeclString == "":

            varDeclString += varName

        else:

            varDeclString += "," + varName

    ### Variables declaration

    varsList = var(varDeclString)

    sage_eval("var('" + varDeclString + "')")

    ## Create the equations

    equationString = ""

    equationsList = []

    for vIndex in xrange(0, len(vect)):

        equationString = ""

        for bIndex in basisIterationsRange:

            if equationString != "":

                equationString += "+"

            equationString += str(varsList[bIndex]) + "*" + str(basis[bIndex,vIndex])

        equationString += " == " + str(vect[vIndex])

        equationsList.append(sage_eval(equationString))

    ## Solve the equations system. The solution dictionary is needed

    #  to recover the values of the solutions.

    solutions = solve(equationsList,varsList, solution_dict=True)

    ## This code is deactivated: did not solve the problem.

    """

    if len(solutions) == 0:

        print "Trying Maxima to_poly_sove."

        solutions = solve(equationsList,varsList, solution_dict=True, to_poly_solve = 'force')

    """

    #print "Solutions:", s

    ## Recover solutions in rational components vector.

    vectInBasis = vector(QQ, basis.nrows())

    ### There can be several solution, in the general case.

    #   Here, there is only one (or none). For each solution, each 

    #   variable has to be searched for.

    for sol in solutions:

        for bIndex in basisIterationsRange:

            vectInBasis[bIndex] = sol[varsList[bIndex]]

    return vectInBasis

# End cpv_vector_in_basis.

#

sys.stderr.write("... functions for CVP loaded.\n")