1 /*
2 Copyright (C) 2014, 2015 Alex J. Best
3
4 This file is part of FLINT.
5
6 FLINT is free software: you can redistribute it and/or modify it under
7 the terms of the GNU Lesser General Public License (LGPL) as published
8 by the Free Software Foundation; either version 2.1 of the License, or
9 (at your option) any later version. See <http://www.gnu.org/licenses/>.
10 */
11
12 #include "fmpz_mat.h"
13 #include "fmpq_mat.h"
14 #include "perm.h"
15
16 static void
add_columns(fmpz_mat_t H,const fmpz_mat_t B,const fmpz_mat_t H1,flint_rand_t state)17 add_columns(fmpz_mat_t H, const fmpz_mat_t B, const fmpz_mat_t H1, flint_rand_t state)
18 {
19 int neg;
20 slong i, j, n, bits;
21 fmpz_t den, tmp, one;
22 fmpq_t num, alpha;
23 fmpz_mat_t Bu, B1, cols, k;
24 fmpq_mat_t H1_q, cols_q, x;
25
26 n = B->r;
27
28 fmpz_mat_init(Bu, n, n);
29 fmpz_mat_init(B1, n - 1, n);
30 fmpz_mat_init(cols, n, B->c - n);
31 fmpz_mat_init(k, n, 1);
32 fmpq_mat_init(x, n, B->c - n);
33 fmpq_mat_init(cols_q, n, B->c - n);
34 fmpq_mat_init(H1_q, n, n);
35
36 for (i = 0; i < n; i++)
37 for (j = 0; j < cols->c; j++)
38 fmpz_set(fmpz_mat_entry(cols, i, j), fmpz_mat_entry(B, i, n + j));
39 for (i = 0; i < n - 1; i++)
40 {
41 for (j = 0; j < n; j++)
42 {
43 fmpz_set(fmpz_mat_entry(Bu, i, j), fmpz_mat_entry(B, i, j));
44 fmpz_set(fmpz_mat_entry(B1, i, j), fmpz_mat_entry(B, i, j));
45 }
46 }
47
48 /* find kernel basis vector */
49 if (fmpz_mat_nullspace(k, B1) != 1)
50 {
51 flint_printf("Exception (fmpz_mat_hnf_pernet_stein). "
52 "Nullspace was not dimension one.\n");
53 flint_abort();
54 }
55
56 bits = fmpz_mat_max_bits(B1);
57 if (bits < 0)
58 bits = -bits;
59
60 fmpz_mat_clear(B1);
61
62 fmpz_init(tmp);
63
64 /* set the last row of Bu to be random, such that Bu is nonsingular */
65 while (fmpz_is_zero(tmp))
66 {
67 _fmpz_vec_randtest(Bu->rows[n - 1], state, n, bits);
68
69 fmpz_zero(tmp);
70 for (j = 0; j < n; j++)
71 fmpz_addmul(tmp, fmpz_mat_entry(Bu, n - 1, j),
72 fmpz_mat_entry(k, j, 0));
73 }
74 fmpz_clear(tmp);
75
76 /* solve Bu*x = cols */
77 if (!fmpq_mat_solve_fmpz_mat(x, Bu, cols))
78 {
79 flint_printf("Exception (fmpz_mat_hnf_pernet_stein). "
80 "Singular input matrix for solve.");
81 flint_abort();
82 }
83
84 /* fix final row */
85 fmpq_init(num);
86 fmpz_init(den);
87 fmpq_init(alpha);
88 fmpz_init(one);
89 fmpz_one(one);
90
91 /* compute denominator */
92 for (i = 0; i < n; i++)
93 fmpz_addmul(den, fmpz_mat_entry(B, n - 1, i), fmpz_mat_entry(k, i, 0));
94 neg = (fmpz_sgn(den) < 0);
95 if (neg)
96 fmpz_neg(den, den);
97
98 for (j = 0; j < B->c - H1->c; j++)
99 {
100 fmpq_zero(num);
101 for (i = 0; i < n; i++)
102 {
103 _fmpq_addmul(fmpq_numref(num), fmpq_denref(num),
104 fmpz_mat_entry(B, n - 1, i), one,
105 fmpq_mat_entry_num(x, i, j),
106 fmpq_mat_entry_den(x, i, j));
107 }
108 _fmpq_sub(fmpq_numref(alpha), fmpq_denref(alpha),
109 fmpz_mat_entry(B, n - 1, n + j), one,
110 fmpq_numref(num), fmpq_denref(num));
111
112 _fmpq_mul(fmpq_numref(alpha), fmpq_denref(alpha),
113 fmpq_numref(alpha), fmpq_denref(alpha), one, den);
114 if (neg)
115 fmpq_neg(alpha, alpha);
116
117 /* x_i += alpha*k */
118 for (i = 0; i < n; i++)
119 {
120 _fmpq_addmul(fmpq_mat_entry_num(x, i, j),
121 fmpq_mat_entry_den(x, i, j), fmpq_numref(alpha),
122 fmpq_denref(alpha), fmpz_mat_entry(k, i, 0), one);
123 }
124 }
125
126 fmpq_clear(num);
127 fmpz_clear(den);
128 fmpz_clear(one);
129 fmpq_clear(alpha);
130
131 /* set cols = H1*x and place in position in H */
132 fmpq_mat_set_fmpz_mat(H1_q, H1);
133 fmpq_mat_mul(cols_q, H1_q, x);
134 fmpq_mat_get_fmpz_mat(cols, cols_q);
135 for (i = 0; i < n; i++)
136 {
137 for (j = 0; j < n; j++)
138 fmpz_set(fmpz_mat_entry(H, i, j), fmpz_mat_entry(H1, i, j));
139 for (j = n; j < H->c; j++)
140 fmpz_set(fmpz_mat_entry(H, i, j), fmpz_mat_entry(cols, i, j - n));
141 }
142
143 fmpq_mat_clear(H1_q);
144 fmpq_mat_clear(x);
145 fmpq_mat_clear(cols_q);
146 fmpz_mat_clear(k);
147 fmpz_mat_clear(cols);
148 fmpz_mat_clear(Bu);
149 }
150
151 /* takes input matrix H with rows 0 to start_row - 1 in HNF to a HNF matrix */
152 static void
add_rows(fmpz_mat_t H,slong start_row,slong * pivots,slong num_pivots)153 add_rows(fmpz_mat_t H, slong start_row, slong *pivots, slong num_pivots)
154 {
155 slong i, i2, j, j2, new_row, row;
156 fmpz_t b, d, u, v, r1d, r2d, q;
157
158 fmpz_init(b);
159 fmpz_init(d);
160 fmpz_init(u);
161 fmpz_init(v);
162 fmpz_init(r1d);
163 fmpz_init(r2d);
164 fmpz_init(q);
165
166 for (row = start_row; row < H->r; row++)
167 {
168 /* reduce row to be added with existing */
169 for (i = j = 0; i < num_pivots; i++)
170 {
171 /* check if added row can still be reduced */
172 for (; j < pivots[i]; j++)
173 if (!fmpz_is_zero(fmpz_mat_entry(H, row, j)))
174 break;
175 if (j < pivots[i])
176 break;
177 if (fmpz_is_zero(fmpz_mat_entry(H, row, j)))
178 continue;
179 fmpz_xgcd(d, u, v, fmpz_mat_entry(H, i, j),
180 fmpz_mat_entry(H, row, j));
181 fmpz_divexact(r1d, fmpz_mat_entry(H, i, j), d);
182 fmpz_divexact(r2d, fmpz_mat_entry(H, row, j), d);
183 for (j2 = j; j2 < H->c; j2++)
184 {
185 fmpz_mul(b, u, fmpz_mat_entry(H, i, j2));
186 fmpz_addmul(b, v, fmpz_mat_entry(H, row, j2));
187 fmpz_mul(fmpz_mat_entry(H, row, j2), r1d,
188 fmpz_mat_entry(H, row, j2));
189 fmpz_submul(fmpz_mat_entry(H, row, j2), r2d,
190 fmpz_mat_entry(H, i, j2));
191 fmpz_set(fmpz_mat_entry(H, i, j2), b);
192 }
193 }
194
195 /* find first non-zero entry of the added row */
196 for (j = 0; j < H->c && fmpz_is_zero(fmpz_mat_entry(H, row, j)); j++) ;
197 new_row = row;
198 if (j != H->c) /* last row non-zero, move to correct position */
199 {
200 if (fmpz_sgn(fmpz_mat_entry(H, row, j)) < 0)
201 {
202 for (j2 = j; j2 < H->c; j2++)
203 {
204 fmpz_neg(fmpz_mat_entry(H, row, j2),
205 fmpz_mat_entry(H, row, j2));
206 }
207 }
208 do
209 {
210 if (new_row < row)
211 fmpz_mat_swap_rows(H, NULL, new_row, new_row + 1);
212 if (new_row == 0)
213 break;
214 new_row--;
215 for (j2 = 0; j2 < H->c &&
216 fmpz_is_zero(fmpz_mat_entry(H, new_row, j2)); j2++) ;
217 }
218 while (j2 > j);
219 }
220
221 /* recompute pivots */
222 for (i = new_row, j = 0; i <= row && i < H->c; i++, j++)
223 {
224 for (; j < H->c && fmpz_is_zero(fmpz_mat_entry(H, i, j)); j++) ;
225 if (j == H->c)
226 break;
227 pivots[i] = j;
228 num_pivots = i + 1;
229 }
230
231 /* reduce above pivot entries */
232 for (i = 0; i < num_pivots; i++)
233 {
234 for (i2 = 0; i2 < i; i2++)
235 {
236 fmpz_fdiv_q(q, fmpz_mat_entry(H, i2, pivots[i]),
237 fmpz_mat_entry(H, i, pivots[i]));
238 for (j2 = pivots[i]; j2 < H->c; j2++)
239 {
240 fmpz_submul(fmpz_mat_entry(H, i2, j2), q,
241 fmpz_mat_entry(H, i, j2));
242 }
243 }
244 }
245 }
246
247 fmpz_clear(q);
248 fmpz_clear(r2d);
249 fmpz_clear(r1d);
250 fmpz_clear(v);
251 fmpz_clear(u);
252 fmpz_clear(d);
253 fmpz_clear(b);
254 }
255
256 static void
double_det(fmpz_t d1,fmpz_t d2,const fmpz_mat_t B,const fmpz_mat_t c,const fmpz_mat_t d)257 double_det(fmpz_t d1, fmpz_t d2, const fmpz_mat_t B, const fmpz_mat_t c,
258 const fmpz_mat_t d)
259 {
260 slong i, j, n;
261 slong *P;
262 mp_limb_t p, u1mod, u2mod, v1mod, v2mod;
263 fmpz_t bound, prod, s1, s2, t, u1, u2, v1, v2;
264 fmpz_mat_t dt, Bt;
265 fmpq_t tmpq;
266 fmpq_mat_t x;
267 nmod_mat_t Btmod;
268
269 n = B->c;
270
271 fmpz_mat_init(dt, n, 1);
272 fmpz_mat_init(Bt, n, n);
273 fmpq_mat_init(x, n, 1);
274
275 for (i = 0; i < n; i++)
276 {
277 for (j = 0; j < n - 1; j++)
278 fmpz_set(fmpz_mat_entry(Bt, i, j), fmpz_mat_entry(B, j, i));
279 fmpz_set(fmpz_mat_entry(Bt, i, n - 1), fmpz_mat_entry(c, 0, i));
280 }
281
282 /* solve B^Tx = d^T */
283 fmpz_mat_transpose(dt, d);
284 fmpq_mat_solve_fmpz_mat(x, Bt, dt);
285
286 if (!fmpq_is_zero(fmpq_mat_entry(x, n - 1, 0)))
287 {
288 fmpz_init(bound);
289 fmpz_init(prod);
290 fmpz_init(t);
291 fmpz_init(s1);
292 fmpz_init(s2);
293 fmpz_init(u1);
294 fmpz_init(u2);
295 fmpz_init(v1);
296 fmpz_init(v2);
297
298 /* compute lcm of denominators of vectors x and y */
299 fmpq_init(tmpq);
300 fmpz_one(u1);
301 fmpz_one(u2);
302 for (i = 0; i < n - 1; i++)
303 {
304 fmpz_lcm(u1, u1, fmpq_mat_entry_den(x, i, 0));
305 fmpq_div(tmpq, fmpq_mat_entry(x, i, 0),
306 fmpq_mat_entry(x, n - 1, 0));
307 fmpz_lcm(u2, u2, fmpq_denref(tmpq));
308 }
309 fmpz_lcm(u1, u1, fmpq_mat_entry_den(x, n - 1, 0));
310 fmpq_inv(tmpq, fmpq_mat_entry(x, n - 1, 0));
311 fmpz_lcm(u2, u2, fmpq_denref(tmpq));
312 fmpq_clear(tmpq);
313
314 /* compute Hadamard bounds */
315 fmpz_one(bound);
316 for (j = 0; j < n - 1; j++)
317 {
318 fmpz_zero(s1);
319 for (i = 0; i < n; i++)
320 fmpz_addmul(s1, fmpz_mat_entry(Bt, i, j),
321 fmpz_mat_entry(Bt, i, j));
322 fmpz_sqrtrem(s1, t, s1);
323 if (!fmpz_is_zero(t))
324 fmpz_add_ui(s1, s1, UWORD(1));
325 fmpz_mul(bound, bound, s1);
326 }
327 fmpz_zero(s1);
328 fmpz_zero(s2);
329 for (j = 0; j < n; j++)
330 {
331 fmpz_addmul(s1, fmpz_mat_entry(c, 0, j), fmpz_mat_entry(c, 0, j));
332 fmpz_addmul(s2, fmpz_mat_entry(d, 0, j), fmpz_mat_entry(d, 0, j));
333 }
334 fmpz_sqrtrem(s1, t, s1);
335 if (!fmpz_is_zero(t))
336 fmpz_add_ui(s1, s1, UWORD(1));
337 fmpz_sqrtrem(s2, t, s2);
338 if (!fmpz_is_zero(t))
339 fmpz_add_ui(s2, s2, UWORD(1));
340 fmpz_mul(s1, s1, bound);
341 fmpz_mul(s2, s2, bound);
342 fmpz_cdiv_q(s1, s1, u1);
343 fmpz_cdiv_q(s2, s2, u2);
344 if (fmpz_cmp(s1, s2) > 0)
345 fmpz_set(bound, s1);
346 else
347 fmpz_set(bound, s2);
348 fmpz_mul_ui(bound, bound, UWORD(2));
349
350 fmpz_one(prod);
351 P = _perm_init(n);
352 nmod_mat_init(Btmod, n, n, 2);
353 p = UWORD(1) << NMOD_MAT_OPTIMAL_MODULUS_BITS;
354 /* compute determinants divided by u1 and u2 */
355 while (fmpz_cmp(prod, bound) <= 0)
356 {
357 p = n_nextprime(p, 0);
358 u1mod = fmpz_fdiv_ui(u1, p);
359 u2mod = fmpz_fdiv_ui(u2, p);
360 if (!(u1mod || u2mod))
361 continue;
362 _nmod_mat_set_mod(Btmod, p);
363 for (i = 0; i < n; i++)
364 {
365 for (j = 0; j < n - 1; j++)
366 nmod_mat_entry(Btmod, i, j) =
367 fmpz_fdiv_ui(fmpz_mat_entry(B, j, i), p);
368 nmod_mat_entry(Btmod, i, n - 1) =
369 fmpz_fdiv_ui(fmpz_mat_entry(c, 0, i), p);
370 }
371 nmod_mat_lu(P, Btmod, 0);
372 v1mod = UWORD(1);
373 for (i = 0; i < n; i++)
374 v1mod = n_mulmod2_preinv(v1mod, nmod_mat_entry(Btmod, i, i), p,
375 Btmod->mod.ninv);
376 if (_perm_parity(P, n) == 1)
377 v1mod = nmod_neg(v1mod, Btmod->mod);
378
379 for (i = 0; i < n; i++)
380 {
381 for (j = 0; j < n - 1; j++)
382 nmod_mat_entry(Btmod, i, j) =
383 fmpz_fdiv_ui(fmpz_mat_entry(B, j, i), p);
384 nmod_mat_entry(Btmod, i, n - 1) =
385 fmpz_fdiv_ui(fmpz_mat_entry(d, 0, i), p);
386 }
387 nmod_mat_lu(P, Btmod, 0);
388 v2mod = UWORD(1);
389 for (i = 0; i < n; i++)
390 v2mod = n_mulmod2_preinv(v2mod, nmod_mat_entry(Btmod, i, i), p,
391 Btmod->mod.ninv);
392 if (_perm_parity(P, n) == 1)
393 v2mod = nmod_neg(v2mod, Btmod->mod);
394
395 v1mod = n_mulmod2_preinv(v1mod, n_invmod(u1mod, p), p,
396 Btmod->mod.ninv);
397 v2mod = n_mulmod2_preinv(v2mod, n_invmod(u2mod, p), p,
398 Btmod->mod.ninv);
399 fmpz_CRT_ui(v1, v1, prod, v1mod, p, 1);
400 fmpz_CRT_ui(v2, v2, prod, v2mod, p, 1);
401 fmpz_mul_ui(prod, prod, p);
402 }
403
404 fmpz_mul(d1, u1, v1);
405 fmpz_mul(d2, u2, v2);
406
407 fmpz_clear(bound);
408 fmpz_clear(prod);
409 fmpz_clear(s1);
410 fmpz_clear(s2);
411 fmpz_clear(u1);
412 fmpz_clear(u2);
413 fmpz_clear(v1);
414 fmpz_clear(v2);
415 fmpz_clear(t);
416 _perm_clear(P);
417 nmod_mat_clear(Btmod);
418 }
419 else /* can't use the clever method above so naively compute both dets */
420 {
421 fmpz_mat_det(d1, Bt);
422 for (j = 0; j < n; j++)
423 fmpz_set(fmpz_mat_entry(Bt, j, n - 1), fmpz_mat_entry(d, 0, j));
424 fmpz_mat_det(d2, Bt);
425 }
426
427 fmpz_mat_clear(dt);
428 fmpz_mat_clear(Bt);
429 fmpq_mat_clear(x);
430 }
431
432 void
fmpz_mat_hnf_pernet_stein(fmpz_mat_t H,const fmpz_mat_t A,flint_rand_t state)433 fmpz_mat_hnf_pernet_stein(fmpz_mat_t H, const fmpz_mat_t A, flint_rand_t state)
434 {
435 slong i, j, m, n, p, r, *P, *pivots, finished;
436 fmpz_t d1, d2, g, s, t;
437 fmpz_mat_t c, d, B, C, H1, H2, H3;
438 nmod_mat_t Amod;
439
440 m = fmpz_mat_nrows(A);
441 n = fmpz_mat_ncols(A);
442
443 if (m == 0 || n == 0)
444 return;
445
446 /* find permutation so we can ensure first rows of H are nonsingular */
447 P = _perm_init(m);
448 pivots = _perm_init(n);
449
450 finished = 0;
451
452 while (!finished)
453 {
454 p = n_randprime(state, NMOD_MAT_OPTIMAL_MODULUS_BITS, 1);
455 nmod_mat_init(Amod, m, n, p);
456
457 fmpz_mat_get_nmod_mat(Amod, A);
458 r = _nmod_mat_rref(Amod, pivots, P);
459
460 nmod_mat_clear(Amod);
461
462 /* rank is zero so matrix is possibly zero too */
463 if (r == 0)
464 {
465 if (fmpz_mat_is_zero(A))
466 {
467 fmpz_mat_zero(H);
468 _perm_clear(P);
469 _perm_clear(pivots);
470 return;
471 }
472 continue;
473 }
474
475 /* if A has full column rank we might wish to use minors based hnf */
476 if (r == n && n < 52)
477 {
478 slong b = fmpz_mat_max_bits(A), cutoff = 52;
479 if (b < 0)
480 b = -b;
481
482 if (b <= 8)
483 cutoff = 35;
484 else if (b <= 32)
485 cutoff = 44;
486 else if (b <= 256)
487 cutoff = 48;
488
489 if (n < cutoff)
490 {
491 fmpz_mat_hnf_minors(H, A);
492 _perm_clear(P);
493 _perm_clear(pivots);
494 return;
495 }
496 }
497
498 fmpz_mat_init(c, 1, r - 1);
499 fmpz_mat_init(d, 1, r - 1);
500 fmpz_mat_init(B, FLINT_MAX(r - 2, 0), r - 1);
501 fmpz_mat_init(C, r - 1, r - 1);
502
503 for (i = 0; i < r - 2; i++)
504 {
505 for (j = 0; j < r - 1; j++)
506 {
507 fmpz_set(fmpz_mat_entry(B, i, j),
508 fmpz_mat_entry(A, P[i], pivots[j]));
509 fmpz_set(fmpz_mat_entry(C, i, j),
510 fmpz_mat_entry(A, P[i], pivots[j]));
511 }
512 fmpz_set(fmpz_mat_entry(C, i, r - 1),
513 fmpz_mat_entry(A, P[i], pivots[r - 1]));
514 }
515 for (j = 0; j < r - 1; j++)
516 {
517 fmpz_set(fmpz_mat_entry(c, 0, j),
518 fmpz_mat_entry(A, P[r - 2], pivots[j]));
519 fmpz_set(fmpz_mat_entry(d, 0, j),
520 fmpz_mat_entry(A, P[r - 1], pivots[j]));
521 }
522
523 fmpz_init(g);
524 fmpz_init(s);
525 fmpz_init(t);
526
527 /* if rank is too low leave g = 0 so we don't try to decompose later */
528 if (r > 2)
529 {
530 fmpz_init(d1);
531 fmpz_init(d2);
532
533 double_det(d1, d2, B, c, d);
534 fmpz_xgcd(g, s, t, d1, d2);
535
536 for (j = 0; j < r - 1; j++)
537 {
538 fmpz_mul(fmpz_mat_entry(C, r - 2, j), s,
539 fmpz_mat_entry(A, P[r - 2], pivots[j]));
540 fmpz_addmul(fmpz_mat_entry(C, r - 2, j), t,
541 fmpz_mat_entry(A, P[r - 1], pivots[j]));
542 }
543
544 fmpz_clear(d2);
545 fmpz_clear(d1);
546 }
547
548 if (!fmpz_is_zero(g)) /* chosen matrix invertible */
549 {
550 fmpz_mat_init(H1, r - 1, r - 1);
551
552 if (COEFF_IS_MPZ(*g) && C->r > 3) /* if g is too big, recurse */
553 fmpz_mat_hnf_pernet_stein(H1, C, state);
554 else /* use modulo determinant algorithm to compute HNF of C */
555 fmpz_mat_hnf_modular(H1, C, g);
556
557 fmpz_mat_clear(B);
558 fmpz_mat_init(B, r - 1, n);
559
560 for (j = 0; j < n; j++)
561 {
562 for (i = 0; i < r - 2; i++)
563 fmpz_set(fmpz_mat_entry(B, i, j),
564 fmpz_mat_entry(A, P[i], pivots[j]));
565 fmpz_mul(fmpz_mat_entry(B, r - 2, j), s,
566 fmpz_mat_entry(A, P[r - 2], pivots[j]));
567 fmpz_addmul(fmpz_mat_entry(B, r - 2, j), t,
568 fmpz_mat_entry(A, P[r - 1], pivots[j]));
569 }
570
571 fmpz_mat_init(H2, r - 1, n);
572 fmpz_mat_init(H3, m + 1, n);
573
574 add_columns(H2, B, H1, state);
575
576 for (i = 0; i < r - 1; i++)
577 for (j = 0; j < n; j++)
578 fmpz_set(fmpz_mat_entry(H3, i, pivots[j]),
579 fmpz_mat_entry(H2, i, j));
580
581 for (i = 1; i <= m - r + 2; i++)
582 for (j = 0; j < n; j++)
583 fmpz_set(fmpz_mat_entry(H3, H3->r - i, j),
584 fmpz_mat_entry(A, P[m - i], j));
585
586 /* check the pivots of H3 are as expected */
587 for (i = 0; i < r - 1; i++)
588 {
589 for (j = 0; j < H3->c && fmpz_is_zero(fmpz_mat_entry(H3, i, j)); j++);
590 if (pivots[i] != j)
591 break;
592 }
593
594 /* the pivots were as expected so our choice of prime was ok */
595 if (i == r - 1)
596 {
597 /* add final rows in */
598 add_rows(H3, r - 1, pivots, r - 1);
599
600 /* fill H with HNF */
601 for (i = 0; i < m; i++)
602 for (j = 0; j < n; j++)
603 fmpz_set(fmpz_mat_entry(H, i, j), fmpz_mat_entry(H3, i, j));
604
605 finished = 1;
606 }
607 /* otherwise we must restart as our random prime gave us incorrect pivots */
608
609 fmpz_mat_clear(H1);
610 fmpz_mat_clear(H2);
611 fmpz_mat_clear(H3);
612 }
613 else
614 {
615 if (r == n) /* if A has full column rank we can use minors based hnf */
616 fmpz_mat_hnf_minors(H, A);
617 else
618 fmpz_mat_hnf_classical(H, A);
619 finished = 1;
620 }
621
622 fmpz_clear(t);
623 fmpz_clear(s);
624 fmpz_clear(g);
625 fmpz_mat_clear(C);
626 fmpz_mat_clear(B);
627 fmpz_mat_clear(c);
628 fmpz_mat_clear(d);
629 }
630
631 _perm_clear(P);
632 _perm_clear(pivots);
633 }
634
635