Uses Algorithm 2.7 from Pohst-Zassenhaus to compute H and U st H = HNF(A) = A*U */
| 66 | H = HNF(A) = A*U |
| 67 | */ |
| 68 | void HNF(matrix& H,matrix& U,const matrix& A) |
| 69 | { |
| 70 | int m=A.size(),n=A[0].size(),r,i,j,k; |
| 71 | |
| 72 | #ifdef VERBOSE |
| 73 | cerr << "HNF m=" << m << ", n=" << n << endl; |
| 74 | #endif |
| 75 | |
| 76 | H=A; |
| 77 | ident(U,n); |
| 78 | r=min(m,n); |
| 79 | i=0; |
| 80 | bool flag=true; |
| 81 | bigint mn,te; |
| 82 | int step=2; |
| 83 | while (flag) |
| 84 | { if (step==2) |
| 85 | { // Step 2 |
| 86 | k=-1; |
| 87 | mn=bigint(0); |
| 88 | for (j=i; j<n; j++) |
| 89 | { if (H[i][j]!=0 && (abs(H[i][j])<mn || mn == 0)) |
| 90 | { k=j; mn=abs(H[i][j]); } |
| 91 | } |
| 92 | if (k!=-1) |
| 93 | { if (k!=i) |
| 94 | { // Step 3 |
| 95 | for (j=0; j<m; j++) |
| 96 | { te=H[j][i]; H[j][i]=H[j][k]; H[j][k]=te; } |
| 97 | for (j=0; j<n; j++) |
| 98 | { te=U[j][i]; U[j][i]=U[j][k]; U[j][k]=te; } |
| 99 | } |
| 100 | // Step 4 |
| 101 | bool fl=true; |
| 102 | for (j=i+1; j<n; j++) |
| 103 | { te=H[i][j]/H[i][i]; |
| 104 | if (abs(H[i][j]%H[i][i])>abs(H[i][i]/2)) { te=te+1; } |
| 105 | /* |
| 106 | cout << i << " " << j << " : " ; |
| 107 | cout << H[i][j] << " " << H[i][i] << " " << te << endl; |
| 108 | */ |
| 109 | for (k=0; k<m; k++) { H[k][j]=H[k][j]-te*H[k][i]; } |
| 110 | for (k=0; k<n; k++) { U[k][j]=U[k][j]-te*U[k][i]; } |
| 111 | if (H[i][j]!=0) { fl=false; } |
| 112 | } |
| 113 | if (fl==true) { step=5; } |
| 114 | } |
| 115 | } |
| 116 | |
| 117 | if (step==5) |
| 118 | { // Step 5 |
| 119 | if (H[i][i]<0) |
| 120 | { for (k=0; k<m; k++) { H[k][i]=-H[k][i]; } |
| 121 | for (k=0; k<n; k++) { U[k][i]=-U[k][i]; } |
| 122 | } |
| 123 | for (j=0; j<i; j++) |
| 124 | { te=(H[i][j]/H[i][i]); |
| 125 | for (k=0; k<m; k++) { H[k][j]=H[k][j]-te*H[k][i]; } |