1 | |
---|
2 | #include <NTL/HNF.h> |
---|
3 | |
---|
4 | #include <NTL/new.h> |
---|
5 | |
---|
6 | NTL_START_IMPL |
---|
7 | |
---|
8 | |
---|
9 | // This implements a variation of an algorithm in |
---|
10 | // [P. Domich, R. Kannan and L. Trotter, Math. Oper. Research 12:50-59, 1987]. |
---|
11 | // I started with the description in Henri Cohen's book, but had to modify |
---|
12 | // that because Cohen does not actually keep the numbers reduced modulo |
---|
13 | // the determinant, which leads to larger than necessary numbers. |
---|
14 | // This modifiaction was put in place in v3.9b. |
---|
15 | |
---|
16 | static |
---|
17 | void EuclUpdate(vec_ZZ& u, vec_ZZ& v, |
---|
18 | const ZZ& a, const ZZ& b, const ZZ& c, const ZZ& d, |
---|
19 | const ZZ& M) |
---|
20 | |
---|
21 | { |
---|
22 | long m = u.length(); |
---|
23 | long i; |
---|
24 | |
---|
25 | ZZ M1; |
---|
26 | RightShift(M1, M, 1); |
---|
27 | |
---|
28 | ZZ t1, t2, t3; |
---|
29 | |
---|
30 | for (i = 1; i <= m; i++) { |
---|
31 | mul(t1, u(i), a); |
---|
32 | mul(t2, v(i), b); |
---|
33 | add(t1, t1, t2); |
---|
34 | rem(t1, t1, M); |
---|
35 | if (t1 > M1) |
---|
36 | sub(t1, t1, M); |
---|
37 | |
---|
38 | t3 = t1; |
---|
39 | |
---|
40 | mul(t1, u(i), c); |
---|
41 | mul(t2, v(i), d); |
---|
42 | add(t1, t1, t2); |
---|
43 | rem(t1, t1, M); |
---|
44 | if (t1 > M1) |
---|
45 | sub(t1, t1, M); |
---|
46 | |
---|
47 | u(i) = t3; |
---|
48 | v(i) = t1; |
---|
49 | } |
---|
50 | } |
---|
51 | |
---|
52 | |
---|
53 | static |
---|
54 | void FixDiag(vec_ZZ& u, const ZZ& a, const vec_ZZ& v, const ZZ& M, long m) |
---|
55 | { |
---|
56 | long i; |
---|
57 | ZZ t1; |
---|
58 | |
---|
59 | for (i = 1; i <= m; i++) { |
---|
60 | mul(t1, a, v(i)); |
---|
61 | rem(u(i), t1, M); |
---|
62 | } |
---|
63 | } |
---|
64 | |
---|
65 | |
---|
66 | static |
---|
67 | void ReduceW(vec_ZZ& u, const ZZ& a, const vec_ZZ& v, const ZZ& M, long m) |
---|
68 | { |
---|
69 | long i; |
---|
70 | ZZ t1, t2; |
---|
71 | |
---|
72 | for (i = 1; i <= m; i++) { |
---|
73 | mul(t1, a, v(i)); |
---|
74 | sub(t2, u(i), t1); |
---|
75 | rem(u(i), t2, M); |
---|
76 | } |
---|
77 | } |
---|
78 | |
---|
79 | |
---|
80 | |
---|
81 | void HNF(mat_ZZ& W, const mat_ZZ& A_in, const ZZ& D_in) |
---|
82 | { |
---|
83 | mat_ZZ A = A_in; |
---|
84 | |
---|
85 | long n = A.NumRows(); |
---|
86 | long m = A.NumCols(); |
---|
87 | |
---|
88 | ZZ D = D_in; |
---|
89 | if (D < 0) |
---|
90 | negate(D, D); |
---|
91 | |
---|
92 | if (n == 0 || m == 0 || D == 0) |
---|
93 | Error("HNF: bad input"); |
---|
94 | |
---|
95 | W.SetDims(m, m); |
---|
96 | clear(W); |
---|
97 | |
---|
98 | long i, j, k; |
---|
99 | ZZ d, u, v, c1, c2; |
---|
100 | |
---|
101 | k = n; |
---|
102 | |
---|
103 | for (i = m; i >= 1; i--) { |
---|
104 | for (j = k-1; j >= 1; j--) { |
---|
105 | if (A(j, i) != 0) { |
---|
106 | XGCD(d, u, v, A(k, i), A(j, i)); |
---|
107 | div(c1, A(k, i), d); |
---|
108 | div(c2, A(j, i), d); |
---|
109 | negate(c2, c2); |
---|
110 | EuclUpdate(A(j), A(k), c1, c2, v, u, D); |
---|
111 | } |
---|
112 | } |
---|
113 | |
---|
114 | XGCD(d, u, v, A(k, i), D); |
---|
115 | FixDiag(W(i), u, A(k), D, i); |
---|
116 | if (W(i, i) == 0) W(i, i) = D; |
---|
117 | |
---|
118 | for (j = i+1; j <= m; j++) { |
---|
119 | div(c1, W(j, i), W(i, i)); |
---|
120 | ReduceW(W(j), c1, W(i), D, i); |
---|
121 | } |
---|
122 | |
---|
123 | div(D, D, d); |
---|
124 | k--; |
---|
125 | } |
---|
126 | } |
---|
127 | |
---|
128 | NTL_END_IMPL |
---|