"Fossies" - the Fresh Open Source Software Archive

Member "NZMATH-1.2.0/nzmath/lattice.py" (19 Nov 2012, 6653 Bytes) of package /linux/misc/old/NZMATH-1.2.0.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Python source code syntax highlighting (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file. For more information about "lattice.py" see the Fossies "Dox" file reference documentation.

    1 from __future__ import division
    2 from math import floor
    3 from nzmath.matrix import Matrix
    4 from nzmath.matrix import VectorsNotIndependent
    5 from nzmath.vector import *
    6 from nzmath.arith1 import *
    7 from nzmath.gcd import *
    8 #import nzmath.matrix as matrix
    9 
   10 class Lattice:
   11     """
   12     A class of lattice.
   13     """
   14     def __init__(self, basis, quadraticForm):
   15         self.basis = basis.copy()  # in form of Matrix
   16         self.quadraticForm = quadraticForm.copy()  # in form of Matrix
   17         if self.basis.determinant() == 0:
   18             raise ValueError
   19 
   20     def createElement(self, compo):
   21         return LatticeElement(self, compo)
   22 
   23     def bilinearForm(self, v1, v2):
   24         return v2.transpose() * self.quadraticForm * v1
   25 
   26     def isCyclic(self):
   27         """
   28         Check whether given lattice is a cyclic lattice. 
   29         """
   30         n = self.basis.column
   31         def rot(x):
   32             x_list = []
   33             for i in range(n):
   34                 x_list.append(x[(n-1+i)%n])
   35             return Vector(x_list)
   36         Rot = []
   37         for i in range(n):
   38             X_list = []
   39             for j in range(n):
   40                 X_list.append(self.basis.compo[j][i])
   41             Rot.append(rot(X_list))
   42         T = self.basis.inverse()*Matrix(n, n, Rot)
   43         for i in range(n):
   44             for j in range(n):
   45                 if T.compo[i][j].denominator != 1:
   46                     return False
   47         return True
   48 
   49     def isIdeal(self):
   50         """
   51         Check whether given lattice is a ideal lattice. 
   52         """
   53         n = self.basis.column
   54         Compo = []
   55         for i in range(n):
   56             for j in range(n):
   57                 if i == j + 1:
   58                     Compo.append(1)
   59                 else:
   60                     Compo.append(0)
   61         M = Matrix(n, n, Compo)
   62         d = self.basis.determinant()
   63         B = self.basis.hermiteNormalForm()
   64         z = B.compo[n-1][n-1]
   65         A = B.adjugateMatrix()
   66         z = B.compo[n-1][n-1]
   67         P = A*M*B
   68         for i in range(n):
   69             for j in range(n):
   70                 P.compo[i][j] = P.compo[i][j]%d
   71         sum = 0
   72         for i in range(n):
   73             for j in range(n-1):
   74                 sum = sum + P.compo[i][j]
   75         if sum == 0:
   76             c_compo = []
   77             for i in range(n):
   78                 c_compo.append(P.compo[i][n-1])
   79         else:
   80             return False
   81         c_sum = 0
   82         for i in range(n):
   83             c_sum = c_sum + c_compo[i]
   84         if c_sum == 0:
   85             c_compo = []
   86             for i in range(n):
   87                 c_compo.append(d)
   88         c = Vector(c_compo)
   89 
   90         if z == 1:
   91             qstar = c
   92             poly = B*qstar
   93         elif gcd(z, d//z) != 1:
   94             qstar = c
   95             poly = B*qstar
   96         else:
   97             sum0 = 0
   98             for i in range(n):
   99                 sum0 = sum0 + c.compo[i]%z
  100             if sum0 == 0:
  101                 qstar_compo = []
  102                 for j in range(n):
  103                     qstar_compo.append(CRT([(c.compo[j]//z, d//z), (0, z)]))
  104                 qstar = Vector(qstar_compo)
  105                 poly = B*qstar
  106             else:
  107                 return False
  108 
  109         sum1 = 0
  110         for i in range(n):
  111             sum1 = sum1 + poly.compo[i]%(d//z)
  112         if sum1 == 0:
  113             q_compo = []
  114             for j in range(n):
  115                 q_compo.append(poly.compo[j]//d)
  116             q = Vector(q_compo)
  117             return True, q
  118         else:
  119             return False
  120 
  121 def LLL(_basis, quadraticForm=None ,delta = 0.75):
  122     """
  123     LLL returns LLL-reduced basis
  124     """
  125     basis = _basis.copy()
  126     k=2 
  127     kmax = 1
  128     mu = []
  129     for i in range(basis.column + 1):
  130         mu.append( [0] * i)
  131 
  132     def _innerProduct(v1, v2):
  133         if quadraticForm:
  134             val =  ((v1 * quadraticForm) * v2.toMatrix().transpose())
  135             return val.compo[0]
  136         else:
  137             return innerProduct(v1, v2)
  138     
  139     bstar = [0] * ( basis.column + 1)
  140     bstar[1] = basis[1].copy()
  141     B = [0] * (basis.column + 1)
  142     B[1] = _innerProduct(basis[1], basis[1])
  143     H = basis.getRing().unitMatrix()
  144     
  145     def _RED(k, l):
  146         if 2 * abs(mu[k][l]) <= 1:
  147             return
  148         q = int( floor(0.5 + mu[k][l]) )
  149         basis[k] -= q * basis[l]
  150         H[k] -= q * H[l] 
  151         mu[k][l] -= q
  152         for i in range(1, l):
  153             mu[k][i] -= q * mu[l][i]
  154         return
  155 
  156     def _SWAP(k):
  157         basis.swapColumn(k, k-1)
  158         H.swapColumn(k, k-1)
  159         if k > 2:
  160             for j in range(1, k-1):
  161                 mu[k][j], mu[k-1][j] = mu[k-1][j], mu[k][j]
  162         _mu = mu[k][k-1]
  163         _B = B[k] + _mu * _mu * B[k-1]
  164         
  165         if abs(_B) < (2**(-30)):
  166             B[k], B[k-1] = B[k-1], B[k]
  167             bstar[k], bstar[k-1] = bstar[k-1], bstar[k]
  168             for i in range(k+1, kmax+1):
  169                 mu[i][k], mu[i][k-1] = mu[i][k-1], mu[i][k]
  170         elif abs(B[k]) < (2**(-30)) and _mu != 0:
  171             B[k-1] = _B
  172             bstar[k-1] = _mu * bstar[k-1]
  173             mu[k][k-1] = 1 / _mu
  174             for i in range(k+1, kmax+1):
  175                 mu[i][k-1] = mu[i][k-1] / _mu
  176         elif B[k] != 0:
  177             t = B[k-1] / _B
  178             mu[k][k-1] = _mu * t
  179             _b = bstar[k-1].copy()
  180             bstar[k-1] = bstar[k] + _mu * _b
  181             bstar[k] = -mu[k][k-1]*bstar[k] + (B[k]/_B) * _b
  182             B[k] = B[k] * t
  183             B[k-1] = _B
  184             for i in range(k+1, kmax+1):
  185                 t = mu[i][k]
  186                 mu[i][k] = mu[i][k-1] - _mu * t
  187                 mu[i][k-1]  = t + mu[k][k-1] * mu[i][k]
  188         return
  189     
  190     #step2
  191     while k <= basis.column :
  192         if (k > kmax):
  193             kmax = k
  194             bstar[k] = basis[k].copy()
  195             for j in range(1, k):
  196                 if abs(B[j]) < (2**(-30)):
  197                     mu[k][j] = 0
  198                 else:
  199                     mu[k][j] = _innerProduct(basis[k], bstar[j]) / B[j]
  200                 bstar[k] -= mu[k][j] * bstar[j]
  201             B[k] = _innerProduct(bstar[k], bstar[k])
  202         #step3
  203         while 1:
  204             #print basis
  205             _RED(k,k-1)
  206             if B[k] < (delta - mu[k][k-1] ** 2) * B[k-1]:
  207                 _SWAP(k)
  208                 k = max([2,k-1])
  209             else:
  210                 for l in range(k-2, 0, -1):
  211                     _RED(k, l)
  212                 k += 1
  213                 break
  214     return basis, H
  215 
  216 class LatticeElement(Matrix):
  217 
  218     def __init__(self, lattice, compo):
  219         self.lattice = lattice
  220         self.row = len(compo)
  221         self.column = 1
  222         self.compo = []
  223         for x in compo:
  224             self.compo.append([x])
  225 
  226     def getLattice(self):
  227         return self.lattice
  228