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

#

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_factor(floatBasis, 

                              funcVect, 

                              precision , 

                              precisionFraction = None):

    """

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

    The exponent is made of:

       - a fraction of the precision;

       - an integer depending on the smallest norm of the vectors of the basis.

    """

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

    if precisionFraction is None:

        precisionFraction = 1/2

    ## Compute the norms of the vectors and get the smallest one.

    #  Start with  "oo" (+Infinty)

    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 2^(ceil(precision*precisionFraction))

# End cvp_common_scaling_factor.

#

def cvp_extra_scaling_factor(floatBasis,

                             funcVect,

                             coeffsPrecList,

                             minExponentsList):

    """

    Compute the extra scaling factor for each element of the basis.

    """

    if floatBasis.nrows() == 0 or \

       floatBasis.ncols() == 0 or \

       len(funcVect) == 0:

        return 0

    ## Compute the norm of the function vector.

    funcVectNorm = funcVect.norm()

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

    ###  Pass #1 for norms ratios an minExponentsList 

    scalingPowers = []

    normsRatio = 0

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

        normsRatio =  funcVectNorm / row.norm() 

        normsRatioLog2 = normsRatio.log2().floor()

        print "Current norms ratio:", normsRatio, normsRatioLog2

        scalingPower = normsRatioLog2 + maxExp

        print "Current scaling power:", scalingPower

        scalingPowers.append(scalingPower)

    ### Pass #2 for coefficients precision.

    minScalingPower = oo

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

        scalingPowers[index] -= coeffsPrecList[index]

        if scalingPowers[index] < minScalingPower:

            minScalingPower = scalingPowers[index]

            print "Current scaling power:", scalingPowers[index]

    if minScalingPower < 0:

        scalingFactor = 2^(-minScalingPower)

    else:

        scalingFactor = 1

    return (scalingFactor, scalingPowers)

# End cvp_extra_scaling_factor.

#

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, commonScalingFactor):

    """

    From the float basis and the common scaling factor, compute the

    integral basis.

    """

    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] * commonScalingFactor

    return intBasis

# End cvp_int_basis.

#

def cvp_int_vector_for_approx(floatVect, commonScalingFactor, extraScalingFactor):

    totalScalingFactor = commonScalingFactor * extraScalingFactor

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