"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