Ticket #174: hnfrow.sage

File hnfrow.sage, 7.2 KB (added by was, 14 years ago)
Line 
1"""
2Modular algorithm to compute Hermite normal forms of integer matrices.
3
4AUTHORS:
5    -- Clement Pernet and William Stein (2008-02-07): initial version
6
7TODO:
8
9   [ ] implement the solve using all but last row and "one dimensional affine space" trick.
10
11   [ ] make hermite mod n fast and work.
12
13   [ ] make add_row fast(er)
14
15   [ ] double det
16   
17"""
18
19fast = True
20k = 0
21
22def fastdet(A, times=3):
23   """
24   Likely to be correct fast det (no guarantee).
25   """
26   #t = verbose("start", level=0)
27   v = random_matrix(ZZ,A.nrows(), 1, x=-1,y=1)
28   #t = verbose("got random matrix", t, level=0)   
29   w = A.solve_right(v, check_rank=False)
30   #t = verbose("solved", t, level=0)
31   d = w.denominator()
32   p = 46337
33   for i in range(times):
34      #t = verbose("changing ring", t, level=0)
35      Amod = A._reduce(p)
36      #t = verbose("change ring took: ", t, level=0)
37      #t = verbose("det mod %s"%p, t, level=0)           
38      det = Amod.determinant() / d
39      #t = verbose("done", t, level=0)                 
40      m = ZZ(det)
41      if m >= p//2:
42         m -= p
43      d *= m
44      p = previous_prime(p)
45   return d
46
47def hnf(A):
48   return A.echelon_form()
49
50def hermite_form_mod (A, d):
51   """
52   Return the Hermite form of A computed modulo d (whatever that means).
53
54   INPUT:
55       A -- an m x n matrix over ZZ
56       d -- a positive integer
57
58   OUTPUT:
59       an m x n matrix over ZZ that
60   """
61   return A.echelon_form(D = d)
62
63
64def doubleDet (A, b, c):
65   """
66   Compute the determinants of the stacked matrices A.stack(b) and A.stack(c).
67
68   INPUT:
69       A -- an (n-1) x n matrix
70       b -- an 1 x n matrix
71       c -- an 1 x n matrix
72
73   OUTPUT:
74       a pair of two integers.
75   """
76   #d1 = A.stack(b).determinant()
77   #d2 = A.stack(c).determinant()
78   d1 = fastdet(A.stack(b))
79   d2 = fastdet(A.stack(c))
80   return (d1,d2)
81
82def add_column(B, H_B, a):
83   """
84   The add column procedure.
85   
86   INPUT:
87       B   -- a non-singular square matrix
88       H_B -- the Hermite normal form of B
89       a   -- a column vector
90
91   OUTPUT:
92       x   -- a vector such that H' = H_B.augment(x) is the HNF of A = B.augment(a).
93   """
94
95   if 1:
96      # We use a direct solve method without inverse.  This
97      # is more clever than what is in Allan Steel's talk and
98      # what is in that paper, in 2 ways -- (1) no inverse need
99      # to be computed, and (2) we cleverly solve a vastly easier
100      # system and recover the solution to the original system.
101     
102      t = verbose('starting optimized solve', level=0)
103     
104      # Here's how:
105      # 1. We make a copy of B but with the last *nasty* row of B replaced
106      #    by a random very nice row.
107      C = copy(B)
108      C[C.nrows()-1] = [1]*C.ncols() #[ZZ.random_element() for _ in range(C.ncols())]
109
110      # 2. Then we find the unique solution to C * x = a
111      #    (todo -- recover from bad case.)
112      x = C.solve_right(a)
113
114      # 3. We next delete the last row of B and find a basis vector k
115      #    for the 1-dimensional kernel.
116      D = B.matrix_from_rows(range(C.nrows()-1))
117      N = D._rational_kernel_iml()     
118      if N.ncols() != 1:
119         print "need to recover gracefully from rank issues with matrix."
120      k = N.matrix_from_columns([0])
121
122      # 4. The sought for solution z to B*z = a is some linear combination
123      #       z = x + alpha*k
124      # and setting w to be the last row of B, this column vector z satisfies
125      #       w * z = a'
126      # where a' is the last entry of a.  Thus
127      #       w * (x + alpha*k) = a'
128      # so    w * x + alpha*w*k = a'
129      # so    alpha*w*k  = a' - w*x.
130      w = B[-1]  # last row of B
131      a_prime = a[-1]
132      lhs = w*k
133      rhs = a_prime - w * x
134      alpha = rhs[0] / lhs[0]
135      z = x + alpha*k
136      #print "z = \n", z
137      #print "Bz = \n", B*z
138      t = verbose('start matrix-vector multiply', t, level=0)     
139      x = H_B.change_ring(QQ)*z
140      verbose('done with matrix-vector multiply', t, level=0)
141      #print "x = \n", x
142      x = x.change_ring(ZZ)
143      return x
144      #print x.list()
145
146     
147   elif 0:
148      #print B[-1]
149      t = verbose('invert', level=0)
150      H_B_inv, denom = H_B._invert_iml(check_invertible=False)
151      t = verbose('got invert', t, level=0)
152      Uinv = (B * H_B_inv / denom).change_ring(ZZ)
153      #print "Uinv = \n", Uinv
154      #print "det = ", Uinv.det()
155      t = verbose('did multiplication', t, level=0)
156      x, d = Uinv._solve_iml(a)
157      #print "d = ", d
158      #global k; k = (Uinv, a)
159      #print "x = ", x
160      t = verbose('solve right', t, level=0)
161      return x.change_ring(ZZ)
162   elif 0:   # the direct solve without inverse -- turns out to suck
163      t = verbose('fast solve right in add_column', level=0)
164      b = B.solve_right(a, check_rank=False)
165      t = verbose('solved right', t, level=0)
166      #global k; k = (B, a)
167      t = verbose('start matrix-vector multiply', t, level=0)     
168      x = H_B.change_ring(QQ)*b
169      verbose('done with matrix-vector multiply', t, level=0)
170      x = x.change_ring(ZZ)
171      #print x.list()
172      return x
173   else:
174      verbose('slow add_column', level=0)
175      A = B.augment(a)
176      H = hnf(A)
177      pivots = H.pivots()
178      x = H.matrix_from_columns([H.ncols() - 1])
179      return x
180
181def add_row(A, b, pivots):
182   """
183   The add row procedure.
184
185   INPUT:
186       A -- an n x (n-1) matrix in Hermite normal form
187       b -- an n x 1 matrix
188       pivots -- sorted list of integers; the pivot positions of A.
189
190   OUTPUT:
191       H -- the Hermite normal form of A.stack(b)
192   """
193   if fast:
194      verbose('fast add_row', level=0)
195      return A._add_row_and_maintain_echelon_form(b.row(0), pivots)
196   else:
197      verbose('slow add_row', level=0)
198      H = A.stack(b)
199      ans = hnf(H)
200      return ans, ans.pivots()
201
202def hermite_form (A):
203   """
204   INPUT:
205       an n x m matrix A over the integers.
206   OUTPUT:
207       the Hermite normal form of A.
208   """
209   
210   n = A.nrows()
211   m = A.ncols()
212   assert n == m
213   B = A.matrix_from_rows(range(m-2)).matrix_from_columns(range(n-1))
214   c = A.matrix_from_rows([m-2]).matrix_from_columns (range(n-1))
215   d = A.matrix_from_rows([m-1]).matrix_from_columns (range(n-1))
216   b = A.matrix_from_columns([n-1]).matrix_from_rows(range(m-2))
217
218   t = verbose('double det', level=0, caller_name='')
219   (d1,d2) = doubleDet (B,c,d)
220   verbose('finished', t, level=0, caller_name='')
221
222   t = verbose('xgcd', level=0, caller_name='')
223   (g,k,l) = d1._xgcd (d2, minimal=True)
224   verbose('finished (g=%s)'%g, t, level=0, caller_name='')
225
226   W = B.stack (k*c + l*d)
227
228   t = verbose('hermite mod d', level=0, caller_name='')
229   H = hermite_form_mod (W, g)
230   verbose('finished', t, level=0, caller_name='')
231
232
233   t = verbose('add_col', level=0, caller_name='')
234   x = add_column(B.stack(k*c+l*d), H, b.stack(matrix(1,1,[k*A[m-2,m-1] + l*A[m-1,m-1]])))
235   verbose('finished', t, level=0, caller_name='')
236
237   Hprime = H.augment(x)
238   pivots = range(Hprime.nrows())
239
240   t = verbose('first add row', level=0, caller_name='')
241   Hprime, pivots = add_row(Hprime, A.matrix_from_rows([m-2]), pivots)
242   verbose('finished', t, level=0, caller_name='')
243
244   t = verbose('second add row', level=0, caller_name='')
245   Hprime, pivots = add_row(Hprime, A.matrix_from_rows([m-1]), pivots)
246   verbose('finished', t, level=0, caller_name='')
247
248   H = Hprime.matrix_from_rows(range(m))
249   return H
250
251   
252   
253
254