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

"""@package functions_for_cvp

Auxiliary functions used for CVP.

"""

print "Functions for CVP loading..."

#

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(redBasis, redBasisGso, vect):

    """

    Closest plane vector implementation as per Babaï.

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

    @param redBasisGSO: the GSO of the previous basis;

    @param vector     : 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(redBasis.rows())-1, -1, -1):

        curRBGSO = redBasisGso.row(vIndex)

        curVect = curVect - \

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

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

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

    #

    ##

    vVect = vector(realField, basisPointsNum)

    for vIndex in xrange(0,basisPointsNum):

        vVect[vIndex] = \

            func(basisPointsList[vIndex])

    return vVect

# End cvp_float_vector_for_approx.

#

def cvp_hkz_reduce(intBasis):

    """

    Thin and simplistic wrapper 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):

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

        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_maxi(funct, 

                                        maxDegree, 

                                        lowerBound, 

                                        upperBound,

                                        precsList, 

                                        realField = None,

                                        contFracMaxErr = None):

    pass

    """

    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.

    """

# End cvp_remez_all_poly_error_func_maxi.

#

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 the (f-truncRemez_d(f)) function has very 

    few zeros, compute the zeros of the derivative instead.

    TODO: change the final behaviour! Now!

    """

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

    ## It may happen that the error function has only one or two zeros.

    #  In this case, we are interested in the variations and we consider the

    #  derivative

    if maxDegree > 3:

        if len(errorFuncTruncZerosSa) <= (maxDegree / 2):

            errorFuncDiffSo = pobyso_diff_so_so(errorFuncSo)

            errorFuncZerosSo = pobyso_dirty_find_zeros_so_so(errorFuncDiffSo,

                                                             intervalSo)

            errorFuncTruncZerosSa = pobyso_float_list_so_sa(errorFuncZerosSo)

            pobyso_clear_obj(errorFuncZerosSo)

            pobyso_clear_obj(errorFuncDiffSo) # Clears funcSo and pTruncSo as well.

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

    ## 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 dictionnary is needed

    #  to recover the values of the solutions.

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

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

#

print "... functions for CVP loaded."