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

    return intBasis

# End cvp_int_basis.

#

def cvp_int_vector_for_approx(floatVect, commonScalingExp, extraScalingExp):

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