commit 9645ea29ee5fea35aaf72ab78e65fc95d3077a34
Author: Serge Torres <Serge.Torres@ens-lyon.fr>
Date:   Wed Dec 23 17:52:32 2015 +0100

    Added: lots of functions.

diff --git a/src/functions_for_cvp.sage b/src/functions_for_cvp.sage
index c447df1..d990a75 100644
--- a/src/functions_for_cvp.sage
+++ b/src/functions_for_cvp.sage
@@ -1,24 +1,245 @@
-def common_scaling_factor(precision, floatBasis, precisionFraction = None):
-   """
-   Compute the common scaling factor (a power of 2).
-   The exponent is made of:
+"""@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()
-       if currentNorm < minBasisVectNorm:
-           minBasisVectNorm = currentNorm
-   powerForMinNorm = floor(log(minBasisVectsNorm)/log2)
-   if powerForMinNorm < 0:
-       return 2^(ceil(precision*precisionFraction) - powerFromMinNorm)
-   else:
-       return 2^(ceil(precision*precisionFraction))
-   
\ No newline at end of file
+    """
+    ## 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."
