Line data Source code
1 : /* Copyright (C) 2000 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation; either version 2 of the License, or (at your option) any later
8 : version. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 :
15 : #include "pari.h"
16 : #include "paripriv.h"
17 :
18 : #define DEBUGLEVEL DEBUGLEVEL_mathnf
19 :
20 : /**************************************************************/
21 : /** **/
22 : /** HERMITE NORMAL FORM REDUCTION **/
23 : /** **/
24 : /**************************************************************/
25 : static GEN ZV_hnfgcdext(GEN A);
26 : static GEN
27 21 : hnfallgen(GEN x)
28 : {
29 21 : GEN z = cgetg(3, t_VEC);
30 21 : gel(z,1) = RgM_hnfall(x, (GEN*)(z+2), 1);
31 21 : return z;
32 : }
33 : GEN
34 1428 : mathnf0(GEN x, long flag)
35 : {
36 1428 : switch(typ(x))
37 : {
38 1211 : case t_VEC:
39 1211 : if (RgV_is_ZV(x))
40 1183 : switch (flag)
41 : {
42 14 : case 0:
43 14 : if (lg(x) == 1) return cgetg(1, t_MAT);
44 7 : retmkmat(mkcol(ZV_content(x)));
45 1162 : case 1:
46 : case 4:
47 1162 : return ZV_hnfgcdext(x);
48 : }
49 35 : x = gtomat(x); break;
50 217 : case t_MAT: break;
51 0 : default: pari_err_TYPE("mathnf0",x);
52 : }
53 :
54 252 : switch(flag)
55 : {
56 196 : case 0: case 2: return RgM_is_ZM(x)? ZM_hnf(x): RgM_hnfall(x,NULL,1);
57 35 : case 1: case 3: return RgM_is_ZM(x)? hnfall(x): hnfallgen(x);
58 7 : case 4: RgM_check_ZM(x, "mathnf0"); return hnflll(x);
59 14 : case 5: RgM_check_ZM(x, "mathnf0"); return hnfperm(x);
60 0 : default: pari_err_FLAG("mathnf");
61 : }
62 : return NULL; /* LCOV_EXCL_LINE */
63 : }
64 :
65 : /*******************************************************************/
66 : /* */
67 : /* SPECIAL HNF (FOR INTERNAL USE !!!) */
68 : /* */
69 : /*******************************************************************/
70 : static int
71 10218390 : count(GEN mat, long row, long len, long *firstnonzero)
72 : {
73 10218390 : long j, n = 0;
74 :
75 776282317 : for (j=1; j<=len; j++)
76 : {
77 768260317 : long p = mael(mat,j,row);
78 768260317 : if (p)
79 : {
80 27307028 : if (labs(p)!=1) return -1;
81 25110638 : n++; *firstnonzero=j;
82 : }
83 : }
84 8022000 : return n;
85 : }
86 :
87 : static int
88 555329 : count2(GEN mat, long row, long len)
89 : {
90 : long j;
91 5218497 : for (j=len; j; j--)
92 5047682 : if (labs(mael(mat,j,row)) == 1) return j;
93 170815 : return 0;
94 : }
95 :
96 : static GEN
97 303844 : hnffinal(GEN matgen,GEN perm,GEN* ptdep,GEN* ptB,GEN* ptC)
98 : {
99 303844 : GEN p1, p2, H, Hnew, Bnew, Cnew, diagH1, U = NULL; /* -Wall */
100 303844 : GEN B = *ptB, C = *ptC, dep = *ptdep, depnew;
101 : pari_sp av;
102 : long i,j,k,s,i1,j1,zc;
103 303844 : long co = lg(C);
104 303844 : long col = lg(matgen)-1;
105 : long lnz, nlze, lig;
106 :
107 303844 : if (col == 0) return matgen;
108 303844 : lnz = nbrows(matgen); /* was called lnz-1 - nr in hnfspec */
109 303844 : nlze = nbrows(dep);
110 303844 : lig = nlze + lnz;
111 : /* H: lnz x lnz [disregarding initial 0 cols], U: col x col */
112 303844 : H = ZM_hnflll(matgen, &U, 0);
113 303844 : H += lg(H)-1 - lnz; H[0] = evaltyp(t_MAT) | _evallg(lnz+1);
114 : /* Only keep the part above the H (above the 0s is 0 since the dep rows
115 : * are dependent from the ones in matgen) */
116 303844 : zc = col - lnz; /* # of 0 columns, correspond to units */
117 303844 : if (nlze) { dep = ZM_mul(dep,U); dep += zc; }
118 :
119 303844 : diagH1 = new_chunk(lnz+1); /* diagH1[i] = 0 iff H[i,i] != 1 (set later) */
120 :
121 303844 : av = avma;
122 303844 : Cnew = cgetg(co, typ(C));
123 303844 : setlg(C, col+1); p1 = gmul(C,U);
124 2621880 : for (j=1; j<=col; j++) gel(Cnew,j) = gel(p1,j);
125 4369205 : for ( ; j<co ; j++) gel(Cnew,j) = gel(C,j);
126 :
127 : /* Clean up B using new H */
128 1267826 : for (s=0,i=lnz; i; i--)
129 : {
130 963982 : GEN Di = gel(dep,i), Hi = gel(H,i);
131 963982 : GEN h = gel(Hi,i); /* H[i,i] */
132 963982 : if ( (diagH1[i] = is_pm1(h)) ) { h = NULL; s++; }
133 20855998 : for (j=col+1; j<co; j++)
134 : {
135 19892016 : GEN z = gel(B,j-col);
136 19892016 : p1 = gel(z,i+nlze);
137 19892016 : if (h) p1 = truedivii(p1,h);
138 19892016 : if (!signe(p1)) continue;
139 30803638 : for (k=1; k<=nlze; k++) gel(z,k) = subii(gel(z,k), mulii(p1, gel(Di,k)));
140 148662591 : for ( ; k<=lig; k++) gel(z,k) = subii(gel(z,k), mulii(p1, gel(Hi,k-nlze)));
141 12084347 : gel(Cnew,j) = gsub(gel(Cnew,j), gmul(p1, gel(Cnew,i+zc)));
142 : }
143 963982 : if (gc_needed(av,2))
144 : {
145 746 : if(DEBUGMEM>1) pari_warn(warnmem,"hnffinal, i = %ld",i);
146 746 : (void)gc_all(av, 2, &Cnew, &B);
147 : }
148 : }
149 303844 : p1 = cgetg(lnz+1,t_VEC); p2 = perm + nlze;
150 1267826 : for (i1=0, j1=lnz-s, i=1; i<=lnz; i++) /* push the 1 rows down */
151 963982 : if (diagH1[i])
152 407423 : gel(p1,++j1) = gel(p2,i);
153 : else
154 556559 : gel(p2,++i1) = gel(p2,i);
155 711267 : for (i=i1+1; i<=lnz; i++) gel(p2,i) = gel(p1,i);
156 :
157 : /* s = # extra redundant generators taken from H
158 : * zc col-s co zc = col - lnz
159 : * [ 0 |dep | ] i = nlze + lnz - s = lig - s
160 : * nlze [--------| B' ]
161 : * [ 0 | H' | ] H' = H minus the s rows with a 1 on diagonal
162 : * i [--------|-----] lig-s (= "1-rows")
163 : * [ 0 | Id ]
164 : * [ | ] li */
165 303844 : lig -= s; col -= s; lnz -= s;
166 303844 : Hnew = cgetg(lnz+1,t_MAT);
167 303844 : depnew = cgetg(lnz+1,t_MAT); /* only used if nlze > 0 */
168 303844 : Bnew = cgetg(co-col,t_MAT);
169 303844 : C = shallowcopy(Cnew);
170 1267826 : for (j=1,i1=j1=0; j<=lnz+s; j++)
171 : {
172 963982 : GEN z = gel(H,j);
173 963982 : if (diagH1[j])
174 : { /* hit exactly s times */
175 407423 : i1++; C[i1+col] = Cnew[j+zc];
176 407423 : p1 = cgetg(lig+1,t_COL); gel(Bnew,i1) = p1;
177 651367 : for (i=1; i<=nlze; i++) gel(p1,i) = gcoeff(dep,i,j);
178 407423 : p1 += nlze;
179 : }
180 : else
181 : {
182 556559 : j1++; C[j1+zc] = Cnew[j+zc];
183 556559 : p1 = cgetg(lnz+1,t_COL); gel(Hnew,j1) = p1;
184 556559 : depnew[j1] = dep[j];
185 : }
186 3774401 : for (i=k=1; k<=lnz; i++)
187 2810419 : if (!diagH1[i]) p1[k++] = z[i];
188 : }
189 4369205 : for (j=s+1; j<co-col; j++)
190 : {
191 4065361 : GEN z = gel(B,j-s);
192 4065361 : p1 = cgetg(lig+1,t_COL); gel(Bnew,j) = p1;
193 6465943 : for (i=1; i<=nlze; i++) gel(p1,i) = gel(z,i);
194 4065361 : z += nlze; p1 += nlze;
195 17275853 : for (i=k=1; k<=lnz; i++)
196 13210492 : if (!diagH1[i]) gel(p1,k++) = gel(z,i);
197 : }
198 303844 : *ptdep = depnew;
199 303844 : *ptC = C;
200 303844 : *ptB = Bnew; return Hnew;
201 : }
202 :
203 : /* for debugging */
204 : static void
205 0 : p_mat(GEN mat, GEN perm, long k)
206 : {
207 0 : pari_sp av = avma;
208 0 : perm = vecslice(perm, k+1, lg(perm)-1);
209 0 : err_printf("Permutation: %Ps\n",perm);
210 0 : if (DEBUGLEVEL > 6)
211 0 : err_printf("matgen = %Ps\n", zm_to_ZM( rowpermute(mat, perm) ));
212 0 : set_avma(av);
213 0 : }
214 :
215 : static GEN
216 3771826 : col_dup(long l, GEN col)
217 : {
218 3771826 : GEN c = new_chunk(l);
219 3771826 : memcpy(c,col,l * sizeof(long)); return c;
220 : }
221 :
222 : /* permutation giving imagecompl(x') | image(x'), x' = transpose of x */
223 : static GEN
224 303844 : ZM_rowrankprofile(GEN x, long *nlze)
225 : {
226 303844 : pari_sp av = avma;
227 : GEN d, y;
228 : long i, j, k, l, r;
229 :
230 303844 : x = shallowtrans(x); l = lg(x);
231 303844 : (void)new_chunk(l); /* HACK */
232 303844 : d = ZM_pivots(x,&r); set_avma(av);
233 303844 : *nlze = r;
234 303844 : if (!d) return identity_perm(l-1);
235 293370 : y = cgetg(l,t_VECSMALL);
236 1266583 : for (i = j = 1, k = r+1; i<l; i++)
237 973213 : if (d[i]) y[k++] = i; else y[j++] = i;
238 293370 : return y;
239 : }
240 :
241 : /* HNF reduce a relation matrix (column operations + row permutation)
242 : ** Input:
243 : ** mat = (li-1) x (co-1) matrix of long
244 : ** C = r x (co-1) matrix of GEN
245 : ** perm= permutation vector (length li-1), indexing the rows of mat: easier
246 : ** to maintain perm than to copy rows. For columns we can do it directly
247 : ** using e.g. swap(mat[i], mat[j])
248 : ** k0 = integer. The k0 first lines of mat are dense, the others are sparse.
249 : ** Output: cf ASCII art in the function body
250 : **
251 : ** row permutations applied to perm
252 : ** column operations applied to C. IN PLACE
253 : **/
254 : GEN
255 186566 : hnfspec_i(GEN mat0, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, long k0)
256 : {
257 : pari_sp av;
258 : long co, n, s, nlze, lnz, nr, i, j, k, lk0, col, lig, *p;
259 : GEN mat;
260 : GEN p1, p2, matb, matbnew, vmax, matt, T, extramat, B, C, H, dep, permpro;
261 186566 : const long li = lg(perm); /* = lgcols(mat0) */
262 186566 : const long CO = lg(mat0);
263 :
264 186566 : n = 0; /* -Wall */
265 :
266 186566 : C = *ptC; co = CO;
267 186566 : if (co > 300 && co > 1.5 * li)
268 : { /* treat the rest at the end */
269 0 : co = (long)(1.2 * li);
270 0 : setlg(C, co);
271 : }
272 :
273 186566 : if (DEBUGLEVEL>5)
274 : {
275 0 : err_printf("Entering hnfspec\n");
276 0 : p_mat(mat0,perm,0);
277 : }
278 186566 : matt = cgetg(co, t_MAT); /* dense part of mat (top) */
279 186566 : mat = cgetg(co, t_MAT);
280 3958392 : for (j = 1; j < co; j++)
281 : {
282 3771826 : GEN matj = col_dup(li, gel(mat0,j));
283 3771826 : p1 = cgetg(k0+1,t_COL); gel(matt,j) = p1; gel(mat,j) = matj;
284 17528364 : for (i=1; i<=k0; i++) gel(p1,i) = stoi(matj[perm[i]]);
285 : }
286 186566 : av = avma;
287 :
288 186566 : i = lig = li-1; col = co-1; lk0 = k0;
289 186566 : T = (k0 || (lg(C) > 1 && lgcols(C) > 1))? matid(col): NULL;
290 : /* Look for lines with a single nonzero entry, equal to 1 in absolute value */
291 8023320 : while (i > lk0 && col)
292 7836754 : switch( count(mat,perm[i],col,&n) )
293 : {
294 34628 : case 0: /* move zero lines between k0+1 and lk0 */
295 34628 : lk0++; lswap(perm[i], perm[lk0]);
296 34628 : i = lig; continue;
297 :
298 783845 : case 1: /* move trivial generator between lig+1 and li */
299 783845 : lswap(perm[i], perm[lig]);
300 783845 : if (T) swap(gel(T,n), gel(T,col));
301 783845 : swap(gel(mat,n), gel(mat,col)); p = gel(mat,col);
302 783845 : if (p[perm[lig]] < 0) /* = -1 */
303 : { /* convert relation -g = 0 to g = 0 */
304 14526112 : for (i=lk0+1; i<lig; i++) p[perm[i]] = -p[perm[i]];
305 730659 : if (T)
306 : {
307 730659 : p1 = gel(T,col);
308 16035202 : for (i=1; ; i++) /* T = permuted identity: single nonzero entry */
309 16035202 : if (signe(gel(p1,i))) { togglesign_safe(&gel(p1,i)); break; }
310 : }
311 : }
312 783845 : lig--; col--; i = lig; continue;
313 :
314 7018281 : default: i--;
315 : }
316 186566 : if (DEBUGLEVEL>5) { err_printf(" after phase1:\n"); p_mat(mat,perm,0); }
317 :
318 : #define absmax(s,z) {long _z; _z = labs(z); if (_z > s) s = _z;}
319 : /* Get rid of all lines containing only 0 and +/- 1, keeping track of column
320 : * operations in T. Leave the rows 1..lk0 alone [up to k0, coefficient
321 : * explosion, between k0+1 and lk0, row is 0] */
322 186566 : s = 0;
323 1123156 : while (lig > lk0 && col && s < (long)(HIGHBIT>>1))
324 : {
325 2554908 : for (i=lig; i>lk0; i--)
326 2381636 : if (count(mat,perm[i],col,&n) > 0) break;
327 1109862 : if (i == lk0) break;
328 :
329 : /* only 0, +/- 1 entries, at least 2 of them nonzero */
330 936590 : lswap(perm[i], perm[lig]);
331 936590 : swap(gel(mat,n), gel(mat,col)); p = gel(mat,col);
332 936590 : if (T) swap(gel(T,n), gel(T,col));
333 936590 : if (p[perm[lig]] < 0)
334 : {
335 9526089 : for (i=lk0+1; i<=lig; i++) p[perm[i]] = -p[perm[i]];
336 645314 : if (T) ZV_togglesign(gel(T,col));
337 : }
338 25888069 : for (j=1; j<col; j++)
339 : {
340 24951479 : GEN matj = gel(mat,j);
341 : long t;
342 24951479 : if (! (t = matj[perm[lig]]) ) continue;
343 1977651 : if (t == 1) {
344 27641234 : for (i=lk0+1; i<=lig; i++) absmax(s, matj[perm[i]] -= p[perm[i]]);
345 : }
346 : else { /* t = -1 */
347 21429030 : for (i=lk0+1; i<=lig; i++) absmax(s, matj[perm[i]] += p[perm[i]]);
348 : }
349 1977651 : if (T) ZC_lincomb1_inplace(gel(T,j), gel(T,col), stoi(-t));
350 : }
351 936590 : lig--; col--;
352 936590 : if (gc_needed(av,3))
353 : {
354 0 : if(DEBUGMEM>1) pari_warn(warnmem,"hnfspec[1]");
355 0 : if (T) T = gc_GEN(av, T); else set_avma(av);
356 : }
357 : }
358 : /* As above with lines containing a +/- 1 (no other assumption).
359 : * Stop when single precision becomes dangerous */
360 186566 : vmax = cgetg(co,t_VECSMALL);
361 2237957 : for (j=1; j<=col; j++)
362 : {
363 2051391 : GEN matj = gel(mat,j);
364 8966912 : for (s=0, i=lk0+1; i<=lig; i++) absmax(s, matj[i]);
365 2051391 : vmax[j] = s;
366 : }
367 571080 : while (lig > lk0 && col)
368 : {
369 599091 : for (i=lig; i>lk0; i--)
370 555329 : if ( (n = count2(mat,perm[i],col)) ) break;
371 428276 : if (i == lk0) break;
372 :
373 384514 : lswap(vmax[n], vmax[col]);
374 384514 : lswap(perm[i], perm[lig]);
375 384514 : swap(gel(mat,n), gel(mat,col)); p = gel(mat,col);
376 384514 : if (T) swap(gel(T,n), gel(T,col));
377 384514 : if (p[perm[lig]] < 0)
378 : {
379 750360 : for (i=lk0+1; i<=lig; i++) p[perm[i]] = -p[perm[i]];
380 144162 : if (T) ZV_togglesign(gel(T,col));
381 : }
382 4888203 : for (j=1; j<col; j++)
383 : {
384 4503689 : GEN matj = gel(mat,j);
385 : long t;
386 4503689 : if (! (t = matj[perm[lig]]) ) continue;
387 2037514 : if (vmax[col] && (ulong)labs(t) >= (HIGHBIT-vmax[j]) / vmax[col])
388 0 : goto END2;
389 :
390 20910882 : for (s=0, i=lk0+1; i<=lig; i++) absmax(s, matj[perm[i]] -= t*p[perm[i]]);
391 2037514 : vmax[j] = s;
392 2037514 : if (T) ZC_lincomb1_inplace(gel(T,j), gel(T,col), stoi(-t));
393 : }
394 384514 : lig--; col--;
395 384514 : if (gc_needed(av,3))
396 : {
397 0 : if(DEBUGMEM>1) pari_warn(warnmem,"hnfspec[2]");
398 0 : (void)gc_all(av, T? 2: 1, &vmax, &T);
399 : }
400 : }
401 :
402 142804 : END2: /* clean up mat: remove everything to the right of the 1s on diagonal */
403 : /* go multiprecision first */
404 186566 : matb = cgetg(co,t_MAT); /* bottom part (complement of matt) */
405 3958392 : for (j=1; j<co; j++)
406 : {
407 3771826 : GEN matj = gel(mat,j);
408 3771826 : p1 = cgetg(li-k0,t_COL); gel(matb,j) = p1;
409 3771826 : p1 -= k0;
410 101702782 : for (i=k0+1; i<li; i++) gel(p1,i) = stoi(matj[perm[i]]);
411 : }
412 186566 : if (DEBUGLEVEL>5)
413 : {
414 0 : err_printf(" after phase2:\n");
415 0 : p_mat(mat,perm,lk0);
416 : }
417 2105406 : for (i=li-2; i>lig; i--)
418 : {
419 1918840 : long h, i0 = i - k0, k = i + co-li;
420 1918840 : GEN Bk = gel(matb,k);
421 34896401 : for (j=k+1; j<co; j++)
422 : {
423 32977561 : GEN Bj = gel(matb,j), v = gel(Bj,i0);
424 32977561 : s = signe(v); if (!s) continue;
425 :
426 6217114 : gel(Bj,i0) = gen_0;
427 6217114 : if (is_pm1(v))
428 : {
429 3843126 : if (s > 0) /* v = 1 */
430 48960369 : { for (h=1; h<i0; h++) gel(Bj,h) = subii(gel(Bj,h), gel(Bk,h)); }
431 : else /* v = -1 */
432 42472717 : { for (h=1; h<i0; h++) gel(Bj,h) = addii(gel(Bj,h), gel(Bk,h)); }
433 : }
434 : else {
435 55776503 : for (h=1; h<i0; h++) gel(Bj,h) = subii(gel(Bj,h), mulii(v,gel(Bk,h)));
436 : }
437 6217114 : if (T) ZC_lincomb1_inplace(gel(T,j), gel(T,k), negi(v));
438 : }
439 1918840 : if (gc_needed(av,2))
440 : {
441 6 : if(DEBUGMEM>1) pari_warn(warnmem,"hnfspec[3], i = %ld", i);
442 1590 : for (h=1; h<co; h++) setlg(matb[h], i0+1); /* bottom can be forgotten */
443 6 : (void)gc_all(av, T? 2: 1, &matb, &T);
444 : }
445 : }
446 3958392 : for (j=1; j<co; j++) setlg(matb[j], lig-k0+1); /* bottom can be forgotten */
447 186566 : (void)gc_all(av, T? 2: 1, &matb, &T);
448 186566 : if (DEBUGLEVEL>5) err_printf(" matb cleaned up (using Id block)\n");
449 :
450 186566 : nlze = lk0 - k0; /* # of 0 rows */
451 186566 : lnz = lig-nlze+1; /* 1 + # of nonzero rows (!= 0...0 1 0 ... 0) */
452 186566 : if (T) matt = ZM_mul(matt,T); /* update top rows */
453 186566 : extramat = cgetg(col+1,t_MAT); /* = new C minus the 0 rows */
454 1853443 : for (j=1; j<=col; j++)
455 : {
456 1666877 : GEN z = gel(matt,j);
457 1666877 : GEN t = (gel(matb,j)) + nlze - k0;
458 1666877 : p2=cgetg(lnz,t_COL); gel(extramat,j) = p2;
459 7246201 : for (i=1; i<=k0; i++) gel(p2,i) = gel(z,i); /* top k0 rows */
460 2685412 : for ( ; i<lnz; i++) gel(p2,i) = gel(t,i); /* other nonzero rows */
461 : }
462 186566 : if (!col) {
463 0 : permpro = identity_perm(lnz);
464 0 : nr = lnz;
465 : }
466 : else
467 186566 : permpro = ZM_rowrankprofile(extramat, &nr);
468 : /* lnz = lg(permpro) */
469 186566 : if (nlze)
470 : { /* put the nlze 0 rows (trivial generators) at the top */
471 19004 : p1 = new_chunk(lk0+1);
472 53632 : for (i=1; i<=nlze; i++) p1[i] = perm[i + k0];
473 93155 : for ( ; i<=lk0; i++) p1[i] = perm[i - nlze];
474 127783 : for (i=1; i<=lk0; i++) perm[i] = p1[i];
475 : }
476 : /* sort other rows according to permpro (nr redundant generators first) */
477 186566 : p1 = new_chunk(lnz); p2 = perm + nlze;
478 846838 : for (i=1; i<lnz; i++) p1[i] = p2[permpro[i]];
479 846838 : for (i=1; i<lnz; i++) p2[i] = p1[i];
480 : /* perm indexes the rows of mat
481 : * |_0__|__redund__|__dense__|__too big__|_____done______|
482 : * 0 nlze lig li
483 : * \___nr___/ \___k0__/
484 : * \____________lnz ______________/
485 : *
486 : * col co
487 : * [dep | ]
488 : * i0 [--------| B ] (i0 = nlze + nr)
489 : * [matbnew | ] matbnew has maximal rank = lnz-1 - nr
490 : * mat = [--------|-----] lig
491 : * [ 0 | Id ]
492 : * [ | ] li */
493 :
494 186566 : matbnew = cgetg(col+1,t_MAT); /* dense+toobig, maximal rank. For hnffinal */
495 186566 : dep = cgetg(col+1,t_MAT); /* rows dependent from the ones in matbnew */
496 1853443 : for (j=1; j<=col; j++)
497 : {
498 1666877 : GEN z = gel(extramat,j);
499 1666877 : p1 = cgetg(nlze+nr+1,t_COL); gel(dep,j) = p1;
500 1666877 : p2 = cgetg(lnz-nr,t_COL); gel(matbnew,j) = p2;
501 2109943 : for (i=1; i<=nlze; i++) gel(p1,i) = gen_0;
502 1751243 : p1 += nlze; for (i=1; i<=nr; i++) p1[i] = z[permpro[i]];
503 8180370 : p2 -= nr; for ( ; i<lnz; i++) p2[i] = z[permpro[i]];
504 : }
505 :
506 : /* redundant generators in terms of the genuine generators
507 : * (x_i) = - (g_i) B */
508 186566 : B = cgetg(co-col,t_MAT);
509 2291515 : for (j=col+1; j<co; j++)
510 : {
511 2104949 : GEN y = gel(matt,j);
512 2104949 : GEN z = gel(matb,j);
513 2104949 : p1=cgetg(lig+1,t_COL); gel(B,j-col) = p1;
514 3799393 : for (i=1; i<=nlze; i++) gel(p1,i) = gel(z,i);
515 2104949 : p1 += nlze; z += nlze-k0;
516 13537321 : for (k=1; k<lnz; k++)
517 : {
518 11432372 : i = permpro[k];
519 11432372 : gel(p1,k) = (i <= k0)? gel(y,i): gel(z,i);
520 : }
521 : }
522 186566 : if (T) C = typ(C)==t_MAT? RgM_ZM_mul(C,T): RgV_RgM_mul(C,T);
523 186566 : (void)gc_all(av, 4, &matbnew, &B, &dep, &C);
524 186566 : *ptdep = dep;
525 186566 : *ptB = B;
526 186566 : H = hnffinal(matbnew, perm, ptdep, ptB, &C);
527 186566 : if (CO > co)
528 : { /* treat the rest, N cols at a time (hnflll slow otherwise) */
529 0 : const long N = 300;
530 0 : long a, L = CO - co, l = minss(L, N); /* L columns to add */
531 0 : GEN CC = *ptC, m0 = mat0;
532 0 : setlg(CC, CO); /* restore */
533 0 : CC += co-1;
534 0 : m0 += co-1;
535 0 : for (a = l;;)
536 0 : {
537 : GEN MAT, emb;
538 0 : (void)gc_all(av, 4, &H,&C,ptB,ptdep);
539 0 : MAT = cgetg(l + 1, t_MAT);
540 0 : emb = cgetg(l + 1, typ(C));
541 0 : for (j = 1 ; j <= l; j++)
542 : {
543 0 : gel(MAT,j) = gel(m0,j);
544 0 : emb[j] = CC[j];
545 : }
546 0 : H = hnfadd_i(H, perm, ptdep, ptB, &C, MAT, emb);
547 0 : if (a == L) break;
548 0 : CC += l;
549 0 : m0 += l;
550 0 : a += l; if (a > L) { l = L - (a - l); a = L; }
551 : }
552 : }
553 186566 : *ptC = C; return H;
554 : }
555 :
556 : GEN
557 0 : hnfspec(GEN mat, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, long k0)
558 : {
559 0 : pari_sp av = avma;
560 0 : GEN H = hnfspec_i(mat, perm, ptdep, ptB, ptC, k0);
561 0 : return gc_all(av, 4, &H, ptC, ptdep, ptB);
562 : }
563 :
564 : /* HNF reduce x, apply same transforms to C */
565 : GEN
566 0 : mathnfspec(GEN x, GEN *pperm, GEN *pdep, GEN *pB, GEN *pC)
567 : {
568 0 : long i, j, k, l, n, ly, lx = lg(x);
569 : GEN z, v1, perm;
570 0 : if (lx == 1) return cgetg(1, t_MAT);
571 0 : ly = lgcols(x);
572 0 : *pperm = perm = identity_perm(ly-1);
573 0 : z = cgetg(lx,t_MAT);
574 0 : for (i=1; i<lx; i++)
575 : {
576 0 : GEN C = cgetg(ly,t_COL), D = gel(x,i);
577 0 : gel(z,i) = C;
578 0 : for (j=1; j<ly; j++)
579 : {
580 0 : GEN d = gel(D,j);
581 0 : if (is_bigint(d)) goto TOOLARGE;
582 0 : C[j] = itos(d);
583 : }
584 : }
585 : /* [ dep | ]
586 : * [-----| B ]
587 : * [ H | ]
588 : * [-----|-----]
589 : * [ 0 | Id ] */
590 0 : return hnfspec(z,perm, pdep, pB, pC, 0);
591 :
592 0 : TOOLARGE:
593 0 : if (lg(*pC) > 1 && lgcols(*pC) > 1)
594 0 : pari_err_IMPL("mathnfspec with large entries");
595 0 : x = ZM_hnf(x); lx = lg(x);
596 0 : v1 = cgetg(ly, t_VECSMALL);
597 0 : n = lx - ly;
598 0 : for (i = k = l = 1; i < ly; i++)
599 0 : if (equali1(gcoeff(x,i,i + n))) v1[l++] = i; else perm[k++] = i;
600 0 : setlg(perm, k);
601 0 : setlg(v1, l);
602 0 : x = rowpermute(x, perm); /* upper part */
603 0 : *pperm = vecsmall_concat(perm, v1);
604 0 : *pB = vecslice(x, k+n, lx-1);
605 0 : setlg(x, k);
606 0 : *pdep = rowslice(x, 1, n);
607 0 : return n? rowslice(x, n+1, k-1): x; /* H */
608 : }
609 :
610 : /* add new relations to a matrix treated by hnfspec (extramat / extraC) */
611 : GEN
612 117278 : hnfadd_i(GEN H, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, /* cf hnfspec */
613 : GEN extramat,GEN extraC)
614 : {
615 117278 : GEN matb, extratop, Cnew, permpro, B = *ptB, C = *ptC, dep = *ptdep;
616 : long i, lH, lB, li, lig, co, col, nlze;
617 :
618 117278 : if (lg(extramat) == 1) return H;
619 117278 : co = lg(C)-1;
620 117278 : lH = lg(H)-1;
621 117278 : lB = lg(B)-1;
622 117278 : li = lg(perm)-1;
623 117278 : lig = li - lB;
624 117278 : col = co - lB;
625 117278 : nlze = lig - lH;
626 :
627 : /* col co
628 : * [ 0 |dep | ]
629 : * nlze [--------| B ]
630 : * [ 0 | H | ]
631 : * [--------|-----] lig
632 : * [ 0 | Id ]
633 : * [ | ] li */
634 117278 : extratop = zm_to_ZM( rowslicepermute(extramat, perm, 1, lig) );
635 117278 : if (li != lig)
636 : { /* zero out bottom part, using the Id block */
637 117110 : GEN A = vecslice(C, col+1, co);
638 117110 : GEN c = rowslicepermute(extramat, perm, lig+1, li);
639 117110 : extraC = gsub(extraC, typ(A)==t_MAT? RgM_zm_mul(A, c): RgV_zm_mul(A,c));
640 117110 : extratop = ZM_sub(extratop, ZM_zm_mul(B, c));
641 : }
642 :
643 117278 : extramat = shallowconcat(extratop, vconcat(dep, H));
644 117278 : Cnew = shallowconcat(extraC, vecslice(C, col-lH+1, co));
645 117278 : if (DEBUGLEVEL>5) err_printf(" 1st phase done\n");
646 117278 : permpro = ZM_rowrankprofile(extramat, &nlze);
647 117278 : extramat = rowpermute(extramat, permpro);
648 117278 : *ptB = rowpermute(B, permpro);
649 117278 : permpro = vecsmallpermute(perm, permpro);
650 430219 : for (i=1; i<=lig; i++) perm[i] = permpro[i]; /* perm o= permpro */
651 :
652 117278 : *ptdep = rowslice(extramat, 1, nlze);
653 117278 : matb = rowslice(extramat, nlze+1, lig);
654 117278 : if (DEBUGLEVEL>5) err_printf(" 2nd phase done\n");
655 117278 : H = hnffinal(matb,perm,ptdep,ptB,&Cnew);
656 117278 : *ptC = shallowconcat(vecslice(C, 1, col-lH), Cnew);
657 117278 : return H;
658 : }
659 :
660 : GEN
661 0 : hnfadd(GEN H, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, /* cf hnfspec */
662 : GEN extramat,GEN extraC)
663 : {
664 0 : pari_sp av = avma;
665 0 : H = hnfadd_i(H, perm, ptdep, ptB, ptC, ZM_to_zm(extramat), extraC);
666 0 : return gc_all(av, 4, &H, ptC, ptdep, ptB);
667 : }
668 :
669 : /* zero aj = Aij (!= 0) using ak = Aik (maybe 0), via linear combination of
670 : * A[j] and A[k] of determinant 1. If U != NULL, likewise update its columns */
671 : static void
672 41407713 : ZC_elem(GEN aj, GEN ak, GEN A, GEN U, long j, long k)
673 : {
674 : GEN p1,u,v,d;
675 :
676 41407713 : if (!signe(ak)) {
677 111138 : swap(gel(A,j), gel(A,k));
678 111138 : if (U) swap(gel(U,j), gel(U,k));
679 37583260 : return;
680 : }
681 41296575 : d = bezout(aj,ak,&u,&v);
682 : /* frequent special case (u,v) = (1,0) or (0,1) */
683 41296575 : if (!signe(u))
684 : { /* ak | aj */
685 20627298 : p1 = diviiexact(aj,ak); togglesign(p1);
686 20627298 : ZC_lincomb1_inplace(gel(A,j), gel(A,k), p1);
687 20627298 : if (U)
688 2445283 : ZC_lincomb1_inplace(gel(U,j), gel(U,k), p1);
689 20627298 : return;
690 : }
691 20669277 : if (!signe(v))
692 : { /* aj | ak */
693 16844824 : p1 = diviiexact(ak,aj); togglesign(p1);
694 16844824 : ZC_lincomb1_inplace(gel(A,k), gel(A,j), p1);
695 16844824 : swap(gel(A,j), gel(A,k));
696 16844824 : if (U) {
697 440902 : ZC_lincomb1_inplace(gel(U,k), gel(U,j), p1);
698 440902 : swap(gel(U,j), gel(U,k));
699 : }
700 16844824 : return;
701 : }
702 :
703 3824453 : if (!is_pm1(d)) { aj = diviiexact(aj, d); ak = diviiexact(ak, d); }
704 3824453 : p1 = gel(A,k); aj = negi(aj); /* NOT togglesign */
705 3824453 : gel(A,k) = ZC_lincomb(u,v, gel(A,j),p1);
706 3824453 : gel(A,j) = ZC_lincomb(aj,ak, p1,gel(A,j));
707 3824453 : if (U)
708 : {
709 798398 : p1 = gel(U,k);
710 798398 : gel(U,k) = ZC_lincomb(u,v, gel(U,j),p1);
711 798398 : gel(U,j) = ZC_lincomb(aj,ak, p1,gel(U,j));
712 : }
713 : }
714 :
715 : INLINE int
716 2590 : is_RgX(GEN a, long v) { return typ(a) == t_POL && varn(a)==v; }
717 : /* set u,v such that au + bv = gcd(a,b), divide a,b by the gcd */
718 : static GEN
719 882 : gbezout_step(GEN *pa, GEN *pb, GEN *pu, GEN *pv, long vx)
720 : {
721 882 : GEN a = *pa, b = *pb, d, l;
722 882 : if (gequal0(a))
723 : {
724 0 : *pa = gen_0; *pu = gen_0;
725 0 : *pb = gen_1; *pv = gen_1; return b;
726 : }
727 882 : a = is_RgX(a,vx)? RgX_renormalize(a): scalarpol(a, vx);
728 882 : b = is_RgX(b,vx)? RgX_renormalize(b): scalarpol(b, vx);
729 882 : d = RgX_extgcd(a,b, pu,pv);
730 882 : l = pollead(d,vx);
731 882 : d = RgX_Rg_div(d,l);
732 882 : *pu = RgX_Rg_div(*pu,l);
733 882 : *pv = RgX_Rg_div(*pv,l);
734 882 : if (degpol(d)) { a = RgX_div(a, d); b = RgX_div(b, d); }
735 238 : else if (typ(gel(d,2)) == t_REAL && lg(gel(d,2)) <= 3)
736 : #if 1
737 : { /* possible accuracy problem */
738 0 : GEN D = RgX_gcd_simple(a,b);
739 0 : if (degpol(D)) {
740 0 : D = RgX_normalize(D);
741 0 : a = RgX_div(a, D);
742 0 : b = RgX_div(b, D);
743 0 : d = RgX_extgcd(a,b, pu,pv); /* retry now */
744 0 : d = RgX_mul(d, D);
745 : }
746 : }
747 : #else
748 : { /* less stable */
749 : d = RgX_extgcd_simple(a,b, pu,pv);
750 : if (degpol(d)) { a = RgX_div(a, d); b = RgX_div(b, d); }
751 : }
752 : #endif
753 882 : *pa = a;
754 882 : *pb = b; return d;
755 : }
756 : static GEN
757 3369760 : col_mul(GEN x, GEN c)
758 : {
759 3369760 : if (typ(x) == t_INT)
760 : {
761 3367744 : long s = signe(x);
762 3367744 : if (!s) return NULL;
763 2563924 : if (is_pm1(x)) return (s > 0)? c: RgC_neg(c);
764 : }
765 410168 : return RgC_Rg_mul(c, x);
766 : }
767 : static void
768 0 : do_zero(GEN x)
769 : {
770 0 : long i, lx = lg(x);
771 0 : for (i=1; i<lx; i++) gel(x,i) = gen_0;
772 0 : }
773 :
774 : /* (c1, c2) *= [u,-b; v,a] */
775 : static void
776 842440 : update(GEN u, GEN v, GEN a, GEN b, GEN *c1, GEN *c2)
777 : {
778 : GEN p1,p2;
779 :
780 842440 : u = col_mul(u,*c1);
781 842440 : v = col_mul(v,*c2);
782 842440 : if (u) p1 = v? gadd(u,v): u;
783 9768 : else p1 = v? v: NULL;
784 :
785 842440 : a = col_mul(a,*c2);
786 842440 : b = col_mul(gneg_i(b),*c1);
787 842440 : if (a) p2 = b? RgC_add(a,b): a;
788 0 : else p2 = b? b: NULL;
789 :
790 842440 : if (!p1) do_zero(*c1); else *c1 = p1;
791 842440 : if (!p2) do_zero(*c2); else *c2 = p2;
792 842440 : }
793 :
794 : /* zero aj = Aij (!= 0) using ak = Aik (maybe 0), via linear combination of
795 : * A[j] and A[k] of determinant 1. If U != NULL, likewise update its columns */
796 : static void
797 511 : RgC_elem(GEN aj, GEN ak, GEN A, GEN V, long j, long k, long li, long vx)
798 : {
799 511 : GEN u,v, d = gbezout_step(&aj, &ak, &u, &v, vx);
800 : long l;
801 : /* (A[,k], A[,j]) *= [v, -aj; u, ak ] */
802 1778 : for (l = 1; l < li; l++)
803 : {
804 1267 : GEN t = gadd(gmul(u,gcoeff(A,l,j)), gmul(v,gcoeff(A,l,k)));
805 1267 : gcoeff(A,l,j) = gsub(gmul(ak,gcoeff(A,l,j)), gmul(aj,gcoeff(A,l,k)));
806 1267 : gcoeff(A,l,k) = t;
807 : }
808 511 : gcoeff(A,li,j) = gen_0;
809 511 : gcoeff(A,li,k) = d;
810 511 : if (V) update(v,u,ak,aj,(GEN*)(V+k),(GEN*)(V+j));
811 511 : }
812 :
813 : /* reduce A[i,j] mod A[i,j0] for j=j0+1... via column operations */
814 : static void
815 4454982 : ZM_reduce(GEN A, GEN U, long i, long j0)
816 : {
817 4454982 : long j, lA = lg(A);
818 4454982 : GEN d = gcoeff(A,i,j0);
819 4454982 : if (signe(d) < 0)
820 : {
821 38172 : ZV_neg_inplace(gel(A,j0));
822 38172 : if (U) ZV_togglesign(gel(U,j0));
823 38172 : d = gcoeff(A,i,j0);
824 : }
825 9093666 : for (j=j0+1; j<lA; j++)
826 : {
827 4638684 : GEN q = truedivii(gcoeff(A,i,j), d);
828 4638684 : if (!signe(q)) continue;
829 :
830 303115 : togglesign(q);
831 303115 : ZC_lincomb1_inplace(gel(A,j), gel(A,j0), q);
832 303115 : if (U) ZC_lincomb1_inplace(gel(U,j), gel(U,j0), q);
833 : }
834 4454982 : }
835 :
836 : /* normalize T as if it were a t_POL in variable v */
837 : static GEN
838 364 : normalize_as_RgX(GEN T, long v, GEN *pd)
839 : {
840 : GEN d;
841 364 : if (!is_RgX(T,v)) { *pd = T; return gen_1; }
842 336 : d = leading_coeff(T);
843 336 : while (gequal0(d) || (typ(d) == t_REAL && lg(d) == 3
844 0 : && gexpo(T) - expo(d) > (long)BITS_IN_LONG)) {
845 14 : T = normalizepol_lg(T, lg(T)-1);
846 14 : if (!signe(T)) { *pd = gen_1; return T; }
847 0 : d = leading_coeff(T);
848 : }
849 322 : if (degpol(T)) T = RgX_Rg_div(T,d); else { d = gel(T,2); T = gen_1; }
850 322 : *pd = d; return T;
851 : }
852 : /* reduce A[i,j] mod A[i,j0] for j=j0+1... via column operations */
853 : static void
854 84 : RgM_reduce(GEN A, GEN U, long i, long j0, long vx)
855 : {
856 84 : long j, lA = lg(A);
857 84 : GEN d, T = normalize_as_RgX(gcoeff(A,i,j0), vx, &d);
858 84 : if (U && !gequal1(d)) gel(U,j0) = RgC_Rg_div(gel(U,j0), d);
859 84 : gcoeff(A,i,j0) = T;
860 :
861 196 : for (j=j0+1; j<lA; j++)
862 : {
863 112 : GEN t = gcoeff(A,i,j), q;
864 112 : if (gequal0(t)) continue;
865 14 : if (T == gen_1)
866 0 : q = t;
867 14 : else if (is_RgX(t,vx))
868 14 : q = RgX_div(t, T);
869 0 : else continue;
870 :
871 14 : if (gequal0(q)) continue;
872 7 : gel(A,j) = RgC_sub(gel(A,j), RgC_Rg_mul(gel(A,j0), q));
873 7 : if (U) gel(U,j) = RgC_sub(gel(U,j), RgC_Rg_mul(gel(U,j0), q));
874 : }
875 84 : }
876 :
877 : /* A,B square integral in upper HNF, of the same dimension > 0. Return Au
878 : * in Z^n (v in Z^n not computed), such that Au + Bv = [1, 0, ..., 0] */
879 : GEN
880 816035 : hnfmerge_get_1(GEN A, GEN B)
881 : {
882 816035 : pari_sp av = avma;
883 816035 : long j, k, l = lg(A), lb;
884 816035 : GEN b, U = cgetg(l + 1, t_MAT), C = cgetg(l + 1, t_VEC);
885 :
886 816035 : b = gcoeff(B,1,1); lb = lgefint(b);
887 1649669 : for (j = 1; j < l; j++)
888 : {
889 : GEN t;
890 1649669 : long c = j+1;
891 1649669 : gel(U,j) = col_ei(l-1, j);
892 1649669 : gel(U,c) = zerocol(l-1); /* dummy */
893 1649669 : gel(C,j) = vecslice(gel(A,j), 1,j);
894 1649669 : gel(C,c) = vecslice(gel(B,j), 1,j);
895 4552041 : for (k = j; k > 0; k--)
896 : {
897 2902372 : t = gcoeff(C,k,c);
898 2902372 : if (gequal0(t)) continue;
899 2685741 : setlg(C[c], k+1);
900 2685741 : ZC_elem(t, gcoeff(C,k,k), C, U, c, k);
901 2685741 : if (lgefint(gcoeff(C,k,k)) > lb) gel(C,k) = FpC_red(gel(C,k), b);
902 2685741 : if (j > 4)
903 : {
904 94877 : GEN u = gel(U,k);
905 : long h;
906 1405304 : for (h=1; h<l; h++)
907 1310427 : if (lgefint(gel(u,h)) > lb) gel(u,h) = remii(gel(u,h), b);
908 : }
909 : }
910 1649669 : if (j == 1)
911 816035 : t = gcoeff(C,1,1);
912 : else
913 : {
914 : GEN u;
915 833634 : t = bezout(gcoeff(C,1,1), b, &u, NULL); /* >= 0 */
916 833634 : if (signe(u) && !equali1(u)) gel(U,1) = ZC_Z_mul(gel(U,1), u);
917 833634 : gcoeff(C,1,1) = t;
918 : }
919 1649669 : if (equali1(t)) break;
920 : }
921 816035 : if (j >= l) return NULL;
922 816035 : b = lcmii(gcoeff(A,1,1),b);
923 816035 : A = FpC_red(ZM_ZC_mul(A,gel(U,1)), b);
924 816035 : return gc_upto(av, FpC_center(A, b, shifti(b,-1)));
925 : }
926 :
927 : /* remove the first r columns */
928 : static void
929 336479 : remove_0cols(long r, GEN *pA, GEN *pB, long remove)
930 : {
931 336479 : GEN A = *pA, B = *pB;
932 336479 : long l = lg(A);
933 336479 : A += r; A[0] = evaltyp(t_MAT) | _evallg(l-r);
934 336479 : if (B && remove == 2) { B += r; B[0] = A[0]; }
935 336479 : *pA = A; *pB = B;
936 336479 : }
937 :
938 : /* Inefficient compared to hnfall. 'remove' = throw away lin.dep columns */
939 : static GEN
940 62582 : hnf_i(GEN A, int remove)
941 : {
942 62582 : pari_sp av0 = avma, av;
943 : long s, n, m, j, k, li, def, ldef;
944 :
945 62582 : RgM_dimensions(A, &m, &n);
946 62582 : if (!n) return cgetg(1,t_MAT);
947 62323 : av = avma;
948 62323 : A = RgM_shallowcopy(A);
949 62323 : def = n; ldef = (m>n)? m-n: 0;
950 167135 : for (li=m; li>ldef; li--)
951 : {
952 415786 : for (j=def-1; j; j--)
953 : {
954 310974 : GEN a = gcoeff(A,li,j);
955 310974 : if (!signe(a)) continue;
956 :
957 : /* zero a = Aij using b = Aik */
958 160723 : k = (j==1)? def: j-1;
959 160723 : ZC_elem(a,gcoeff(A,li,k), A,NULL, j,k);
960 160723 : if (gc_needed(av,1))
961 : {
962 0 : if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnf[1]. li=%ld",li);
963 0 : A = gc_GEN(av, A);
964 : }
965 : }
966 104812 : s = signe(gcoeff(A,li,def));
967 104812 : if (s)
968 : {
969 102928 : if (s < 0) ZV_neg_inplace(gel(A,def));
970 102928 : ZM_reduce(A, NULL, li,def);
971 102928 : def--;
972 : }
973 : else
974 1884 : if (ldef) ldef--;
975 104812 : if (gc_needed(av,1))
976 : {
977 0 : if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnf[2]. li=%ld",li);
978 0 : A = gc_GEN(av, A);
979 : }
980 : }
981 : /* rank A = n - def */
982 62323 : if (remove) { GEN B = NULL; remove_0cols(def, &A, &B, remove); }
983 62323 : return gc_upto(av0, ZM_copy(A));
984 : }
985 :
986 : GEN
987 90247 : ZM_hnf(GEN x) { return lg(x) > 8? ZM_hnfall(x, NULL, 1): hnf_i(x, 1); }
988 :
989 : /* u*z[1..k] mod p, in place */
990 : static void
991 5420872 : FpV_Fp_mul_part_ip(GEN z, GEN u, GEN p, long k)
992 : {
993 : long i;
994 5420872 : if (is_pm1(u)) {
995 196236 : if (signe(u) > 0) {
996 581948 : for (i = 1; i <= k; i++)
997 441035 : if (signe(gel(z,i))) gel(z,i) = modii(gel(z,i), p);
998 : } else {
999 170371 : for (i = 1; i <= k; i++)
1000 115048 : if (signe(gel(z,i))) gel(z,i) = modii(negi(gel(z,i)), p);
1001 : }
1002 : }
1003 : else {
1004 19714337 : for (i = 1; i <= k; i++)
1005 14489701 : if (signe(gel(z,i))) gel(z,i) = Fp_mul(u,gel(z,i), p);
1006 : }
1007 5420872 : }
1008 : static void
1009 37830162 : FpV_red_part_ipvec(GEN z, GEN p, long k)
1010 : {
1011 : long i;
1012 114789335 : for (i = 1; i <= k; i++) gel(z,i) = modii(gel(z,i), gel(p,i));
1013 37830162 : }
1014 :
1015 : /* return x * U, in echelon form (mod p^m), where (det(U),p) = 1.
1016 : * If early_abort is set, return NULL as soon as one pivot is 0 (mod p^m) */
1017 : GEN
1018 529499 : ZpM_echelon(GEN x, long early_abort, GEN p, GEN pm)
1019 : {
1020 529499 : pari_sp av0 = avma, av;
1021 : long m, li, co, i, j, k, def, ldef;
1022 :
1023 529499 : co = lg(x); if (co == 1) return cgetg(1,t_MAT);
1024 529499 : li = lgcols(x);
1025 529499 : av = avma;
1026 529499 : x = RgM_shallowcopy(x);
1027 529499 : m = Z_pval(pm, p);
1028 :
1029 529499 : ldef = (li > co)? li - co: 0;
1030 2553612 : for (def = co-1,i = li-1; i > ldef; i--)
1031 : {
1032 2024926 : long vmin = LONG_MAX, kmin = 0;
1033 2024926 : GEN umin = gen_0, pvmin, q;
1034 10642048 : for (k = 1; k <= def; k++)
1035 : {
1036 9086691 : GEN u = gcoeff(x,i,k);
1037 : long v;
1038 9086691 : if (!signe(u)) continue;
1039 4222942 : v = Z_pvalrem(u, p, &u);
1040 4222942 : if (v >= m) gcoeff(x,i,k) = gen_0;
1041 2790746 : else if (v < vmin) {
1042 1424726 : vmin = v; kmin = k; umin = u;
1043 1424726 : if (!vmin) break;
1044 : }
1045 : }
1046 2024926 : if (!kmin)
1047 : {
1048 951267 : if (early_abort) return NULL;
1049 950454 : gcoeff(x,i,def) = gen_0;
1050 950454 : ldef--;
1051 950454 : if (ldef < 0) ldef = 0;
1052 950454 : continue;
1053 : }
1054 1073659 : if (kmin != def) swap(gel(x,def), gel(x,kmin));
1055 1073659 : q = vmin? powiu(p, m-vmin): pm;
1056 : /* pivot has valuation vmin */
1057 1073659 : umin = modii(umin, q);
1058 1073659 : if (!equali1(umin))
1059 543681 : FpV_Fp_mul_part_ip(gel(x,def), Fp_inv(umin,q), pm, i-1);
1060 1073659 : gcoeff(x, i, def) = pvmin = powiu(p, vmin);
1061 7332171 : for (j = def-1; j; j--)
1062 : { /* zero x[i, 1..def-1] using x[i,def] = pvmin */
1063 6258512 : GEN t, a = gcoeff(x,i,j) = modii(gcoeff(x,i,j), pm);
1064 6258512 : if (!signe(a)) continue;
1065 :
1066 2940414 : t = diviiexact(a, pvmin); togglesign(t);
1067 2940414 : ZC_lincomb1_inplace(gel(x,j), gel(x,def), t);
1068 2940414 : if (gc_needed(av,1))
1069 : {
1070 14 : if (DEBUGMEM>1) pari_warn(warnmem,"ZpM_echelon. i=%ld",i);
1071 14 : x = gc_GEN(av, x); pvmin = gcoeff(x,i,def);
1072 : }
1073 : }
1074 1073659 : def--;
1075 : }
1076 528686 : if (co > li)
1077 : {
1078 0 : x += co - li;
1079 0 : x[0] = evaltyp(t_MAT) | _evallg(li);
1080 : }
1081 528686 : return gc_GEN(av0, x);
1082 : }
1083 : GEN
1084 4808736 : zlm_echelon(GEN x, long early_abort, ulong p, ulong pm)
1085 : {
1086 4808736 : pari_sp av0 = avma;
1087 : long li, co, i, j, k, def, ldef;
1088 : ulong m;
1089 :
1090 4808736 : co = lg(x); if (co == 1) return cgetg(1,t_MAT);
1091 4808736 : li = lgcols(x);
1092 4808736 : x = Flm_copy(x);
1093 4808736 : m = u_lval(pm, p);
1094 :
1095 4808736 : ldef = (li > co)? li - co: 0;
1096 21403136 : for (def = co-1,i = li-1; i > ldef; i--)
1097 : {
1098 16796992 : long vmin = LONG_MAX, kmin = 0;
1099 16796992 : ulong umin = 0, pvmin, q;
1100 50738129 : for (k = 1; k <= def; k++)
1101 : {
1102 42897701 : ulong u = ucoeff(x,i,k);
1103 : long v;
1104 42897701 : if (!u) continue;
1105 23480861 : v = u_lvalrem(u, p, &u);
1106 23480861 : if (v >= (long) m) ucoeff(x,i,k) = 0;
1107 23480861 : else if (v < vmin) {
1108 16838085 : vmin = v; kmin = k; umin = u;
1109 16838085 : if (!vmin) break;
1110 : }
1111 : }
1112 16796992 : if (!kmin)
1113 : {
1114 1052593 : if (early_abort) return NULL;
1115 850001 : ucoeff(x,i,def) = 0;
1116 850001 : ldef--;
1117 850001 : if (ldef < 0) ldef = 0;
1118 850001 : continue;
1119 : }
1120 15744399 : if (kmin != def) swap(gel(x,def), gel(x,kmin));
1121 15744399 : q = vmin? upowuu(p, m-vmin): pm;
1122 : /* pivot has valuation vmin */
1123 15744399 : umin %= q;
1124 15744399 : if (umin != 1)
1125 7984732 : Flv_Fl_mul_part_inplace(gel(x,def), Fl_inv(umin,q), pm, i-1);
1126 15744399 : ucoeff(x, i, def) = pvmin = upowuu(p, vmin);
1127 69946898 : for (j = def-1; j; j--)
1128 : { /* zero x[i, 1..def-1] using x[i,def] = pvmin */
1129 54202499 : ulong t, a = ucoeff(x,i,j);
1130 54202499 : if (!a) continue;
1131 :
1132 31171171 : t = Fl_neg(a / pvmin, q);
1133 31171171 : Flc_lincomb1_inplace(gel(x,j), gel(x,def), t, pm);
1134 : }
1135 15744399 : def--;
1136 : }
1137 4606144 : if (co > li)
1138 : {
1139 0 : x += co - li;
1140 0 : x[0] = evaltyp(t_MAT) | _evallg(li);
1141 : }
1142 4606144 : return gc_GEN(av0, x);
1143 : }
1144 :
1145 : static int
1146 322729 : ZV_allequal(GEN v)
1147 : {
1148 322729 : long i, l = lg(v);
1149 322729 : if (l > 1)
1150 : {
1151 322729 : GEN x = gel(v,1);
1152 555792 : for (i = 2; i < l; i++) if (!equalii(x,gel(v,i))) return 0;
1153 : }
1154 257719 : return 1;
1155 : }
1156 : /* compute optimal D for hnfmod: x upper triangular */
1157 : static GEN
1158 5195213 : optimal_D(GEN x, GEN D)
1159 : {
1160 5195213 : long i, n = nbrows(x);
1161 5195213 : GEN C = shallowcopy(D);
1162 5195213 : gel(C,1) = gcoeff(x,1,1);
1163 5608653 : for (i = 2; i < n; i++)
1164 : {
1165 2627834 : GEN c = mulii(gel(C,i-1), gcoeff(x,i,i));
1166 2627834 : if (signe(c) < 0) togglesign(c);
1167 2627834 : if (cmpii(c, gel(D,i)) >= 0) break;
1168 413440 : gel(C,i) = c;
1169 : }
1170 5195213 : return C;
1171 : }
1172 :
1173 : /* D = multiple of det x (usually detint(x)) or vector of positive moduli
1174 : * (compute hnf(x | D))
1175 : * flag & hnf_MODID: reduce mod D * matid [ otherwise as above ].
1176 : * flag & hnf_PART: don't reduce once diagonal is known
1177 : * flag & hnf_CENTER: centermod HNF (2|x[i,j]| <] x[i,i]) */
1178 : GEN
1179 5196851 : ZM_hnfmodall_i(GEN x, GEN D, long flag)
1180 : {
1181 : pari_sp av;
1182 5196851 : const long center = (flag & hnf_CENTER);
1183 : long moddiag, modsame, nli, li, co, i, j, k, def, ldef;
1184 : GEN u, LDM;
1185 :
1186 5196851 : co = lg(x);
1187 5196851 : if (co == 1)
1188 : {
1189 1533 : if (typ(D) == t_INT || lg(D) == 1) return cgetg(1,t_MAT);
1190 105 : x = diagonal_shallow(D); /* handle flags properly */
1191 105 : co = lg(x);
1192 : }
1193 5195423 : li = lgcols(x);
1194 5195423 : if (li == 1)
1195 : {
1196 189 : if (typ(D) != t_INT && lg(D) != li) pari_err_DIM("ZM_hnfmod");
1197 189 : return cgetg(1,t_MAT);
1198 : }
1199 5195234 : nli = li - 1;
1200 5195234 : modsame = typ(D)==t_INT;
1201 5195234 : if (!modsame)
1202 : {
1203 322743 : if (lg(D) != li) pari_err_DIM("ZM_hnfmod");
1204 322729 : if (ZV_allequal(D)) { modsame = 1; D = gel(D,1); }
1205 : }
1206 5195220 : moddiag = (flag & hnf_MODID) || !modsame;
1207 : /* modsame: triangularize mod fixed d*Id;
1208 : * moddiag: modulo diagonal matrix, else modulo multiple of determinant */
1209 :
1210 5195220 : if (modsame)
1211 : {
1212 5130210 : LDM = const_vecsmall(nli, 2*lgefint(D)-2);
1213 5130210 : D = const_vec(nli,D);
1214 : }
1215 : else
1216 : {
1217 65010 : LDM = cgetg(li, t_VECSMALL);
1218 276107 : for (i=1; i<li; i++) LDM[i] = lgefint(gel(D,i));
1219 : }
1220 5195220 : av = avma;
1221 5195220 : x = RgM_shallowcopy(x);
1222 :
1223 5195220 : ldef = 0;
1224 5195220 : if (li > co)
1225 : {
1226 146989 : ldef = li - co;
1227 146989 : if (!moddiag)
1228 7 : pari_err_DOMAIN("ZM_hnfmod","nb lines",">", strtoGENstr("nb columns"), x);
1229 : }
1230 20491737 : for (def = co-1,i = nli; i > ldef; i--,def--)
1231 : {
1232 15296524 : GEN d = gel(D,i);
1233 15296524 : long add_N = modsame;
1234 66006465 : for (j = 1; j < def; j++)
1235 : {
1236 50709941 : GEN p1, p2, b, a = gcoeff(x,i,j) = remii(gcoeff(x,i,j), d);
1237 63591799 : if (!signe(a)) continue;
1238 :
1239 33740548 : k = j+1;
1240 33740548 : b = gcoeff(x,i,k) = remii(gcoeff(x,i,k), d);
1241 33740548 : if (!signe(b)) { swap(gel(x,j), gel(x,k)); continue; }
1242 20858690 : if (add_N)
1243 : { /* ensure the moving pivot on row i divides d from now on */
1244 7507330 : add_N = 0;
1245 7507330 : if (!equali1(a))
1246 : { /* x[j] *= u; after this, a = x[i,j] | d */
1247 6110397 : GEN u = Fp_invgen(a, d, &a);
1248 : long t;
1249 6110397 : p1 = gel(x,j);
1250 21522337 : for (t = 1; t < i; t++) gel(p1,t) = mulii(gel(p1,t), u);
1251 6110397 : FpV_red_part_ipvec(p1, D, i-1);
1252 6110397 : gel(p1,i) = a;
1253 6110397 : if (2*lg(a) < lg(b))
1254 : { /* reduce x[i,k] mod x[i,j]: helps ZC_elem */
1255 29053 : GEN r, q = dvmdii(b, a, &r);
1256 29053 : togglesign(q);
1257 29053 : ZC_lincomb1_inplace_i(gel(x,k), gel(x,j), q, i-1);
1258 29053 : FpV_red_part_ipvec(gel(x,k), D, i-1);
1259 29053 : gcoeff(x,i,k) = b = r;
1260 : }
1261 : }
1262 : }
1263 20858690 : ZC_elem(a,b, x, NULL, j,k);
1264 20858690 : p1 = gel(x,j);
1265 20858690 : p2 = gel(x,k);
1266 : /* prevent coeffs explosion */
1267 124780784 : for (k = 1; k < i; k++)
1268 : {
1269 103922094 : if (lgefint(gel(p1,k)) > LDM[k]) gel(p1,k) = remii(gel(p1,k),gel(D,k));
1270 103922094 : if (lgefint(gel(p2,k)) > LDM[k]) gel(p2,k) = remii(gel(p2,k),gel(D,k));
1271 : }
1272 : }
1273 15296524 : if (gc_needed(av,2))
1274 : {
1275 25 : if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnfmod[1]. i=%ld",i);
1276 25 : x = gc_GEN(av, x);
1277 : }
1278 15296524 : if (moddiag && !signe(gcoeff(x,i,def)))
1279 : { /* missing pivot on line i, insert column */
1280 1092533 : GEN a = cgetg(co + 1, t_MAT);
1281 5634343 : for (k = 1; k <= def; k++) gel(a,k) = gel(x,k);
1282 1092533 : gel(a,k++) = Rg_col_ei(gel(D,i), nli, i);
1283 4573183 : for ( ; k <= co; k++) gel(a,k) = gel(x,k-1);
1284 1092533 : ldef--; if (ldef < 0) ldef = 0;
1285 1092533 : co++; def++; x = a;
1286 : }
1287 : }
1288 5195213 : if (co < li)
1289 : { /* implies moddiag, add missing diag(D) components */
1290 113335 : GEN a = cgetg(li+1, t_MAT);
1291 279736 : for (k = 1; k <= li-co; k++) gel(a,k) = Rg_col_ei(gel(D,k), nli, k);
1292 345080 : for (i = 1; i < co; i++) gel(a,k-1+i) = gel(x,i);
1293 113335 : gel(a,li) = zerocol(nli); x = a;
1294 : }
1295 : else
1296 : {
1297 5081878 : x += co - li;
1298 5081878 : x[0] = evaltyp(t_MAT) | _evallg(li); /* kill 0 columns */
1299 5081878 : if (moddiag) x = shallowconcat(x, zerocol(nli));
1300 : }
1301 5195213 : if (moddiag)
1302 : { /* x[li]: extra column, an accumulator discarded at the end */
1303 : GEN D2;
1304 5032550 : gcoeff(x,1,1) = gcdii(gcoeff(x,1,1), gel(D,1));
1305 5032550 : D2 = optimal_D(x,D);
1306 : /* add up missing diag(D) components */
1307 20050553 : for (i = nli; i > 0; i--)
1308 : {
1309 15018003 : gcoeff(x, i, li) = gel(D,i);
1310 57929068 : for (j = i; j > 0; j--)
1311 : {
1312 42911065 : GEN a = gcoeff(x, j, li);
1313 42911065 : if (!signe(a)) continue;
1314 15845356 : ZC_elem(a, gcoeff(x,j,j), x, NULL, li,j);
1315 15845356 : FpV_red_part_ipvec(gel(x,li), D, j-1);
1316 15845356 : FpV_red_part_ipvec(gel(x,j), D, j-1);
1317 : }
1318 15018003 : if (gc_needed(av,1))
1319 : {
1320 35 : if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnfmod[2]. i=%ld", i);
1321 35 : (void)gc_all(av, 2, &x, &D2);
1322 : }
1323 : }
1324 5032550 : D = D2;
1325 : }
1326 : else
1327 : {
1328 162663 : GEN b = gel(D,1);
1329 607585 : for (i = nli; i > 0; i--)
1330 : {
1331 444922 : GEN d = bezout(gcoeff(x,i,i),b, &u,NULL);
1332 444922 : gcoeff(x,i,i) = d;
1333 444922 : FpV_Fp_mul_part_ip(gel(x,i), u, b, i-1);
1334 444922 : if (i > 1) b = diviiexact(b,d);
1335 : }
1336 162663 : D = optimal_D(x,D);
1337 : }
1338 5195213 : x[0] = evaltyp(t_MAT) | _evallg(li); /* kill 0 columns, discard accumulator */
1339 5195213 : if (flag & hnf_PART) return x;
1340 :
1341 20658138 : for (i = nli; i > 0; i--)
1342 : {
1343 15462925 : GEN diag = gcoeff(x,i,i);
1344 15462925 : if (signe(diag) < 0) { gel(x,i) = ZC_neg(gel(x,i)); diag = gcoeff(x,i,i); }
1345 15462925 : if (i != nli)
1346 28553662 : for (j = 1; j < i; j++) gcoeff(x,j,i) = remii(gcoeff(x,j,i), gel(D,j));
1347 44016587 : for (j = i+1; j < li; j++)
1348 : {
1349 28553662 : GEN b = gcoeff(x,i,j) = remii(gcoeff(x,i,j), gel(D,i));
1350 28553662 : GEN r, q = truedvmdii(b, diag, &r);
1351 : /* ensure -diag/2 <= r < diag/2 */
1352 28553662 : if (center && signe(r) && abscmpii(shifti(r,1),diag) >= 0)
1353 672579 : { r = subii(r,diag); q = addiu(q,1); }
1354 28553662 : if (!signe(q)) continue;
1355 11152410 : togglesign(q);
1356 11152410 : ZC_lincomb1_inplace_i(gel(x,j), gel(x,i), q, i-1);
1357 11152410 : gcoeff(x,i,j) = r;
1358 : }
1359 15462925 : if (gc_needed(av,1))
1360 : {
1361 44 : if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnfmod[3]. i=%ld", i);
1362 44 : (void)gc_all(av, 2, &x, &D);
1363 : }
1364 : }
1365 5195213 : return x;
1366 : }
1367 : GEN
1368 5145621 : ZM_hnfmodall(GEN x, GEN dm, long flag)
1369 : {
1370 5145621 : pari_sp av = avma;
1371 5145621 : return gc_GEN(av, ZM_hnfmodall_i(x, dm, flag));
1372 : }
1373 : GEN
1374 162656 : ZM_hnfmod(GEN x, GEN d) { return ZM_hnfmodall(x,d,0); }
1375 : GEN
1376 4542452 : ZM_hnfmodid(GEN x, GEN d)
1377 4542452 : { return ZM_hnfmodall(x,d,hnf_MODID); }
1378 : /* return the column echelon form of x with 1's as pivots,
1379 : * P contains the row indices containing the pivots in increasing order */
1380 : static GEN
1381 3366410 : FpM_echelon(GEN x, GEN *pP, GEN p)
1382 : {
1383 : pari_sp av;
1384 : long iP, li, co, i, j, k, def, ldef;
1385 : GEN P;
1386 :
1387 3366410 : co = lg(x); if (co == 1) { *pP = cgetg(1,t_VECSMALL); return cgetg(1,t_MAT); }
1388 3366410 : li = lgcols(x);
1389 3366410 : iP = 1;
1390 3366410 : *pP = P = cgetg(li, t_VECSMALL);
1391 3366410 : av = avma;
1392 3366410 : x = FpM_red(x, p);
1393 :
1394 3366410 : ldef = (li > co)? li - co: 0;
1395 12687042 : for (def = co-1,i = li-1; i > ldef; i--)
1396 : {
1397 9320632 : GEN u = NULL;
1398 14435837 : for (k = def; k; k--)
1399 : {
1400 10799044 : u = gcoeff(x,i,k);
1401 10799044 : if (signe(u)) break;
1402 : }
1403 9320632 : if (!k)
1404 : {
1405 3636793 : if (--ldef < 0) ldef = 0;
1406 3636793 : continue;
1407 : }
1408 5683839 : P[iP++] = i;
1409 5683839 : if (k != def) swap(gel(x,def), gel(x,k));
1410 5683839 : if (!equali1(u))
1411 4432269 : FpV_Fp_mul_part_ip(gel(x,def), Fp_inv(u,p), p, i-1);
1412 5683839 : gcoeff(x, i, def) = gen_1;
1413 17782981 : for (j = def-1; j; j--)
1414 : { /* zero x[i, 1..def-1] using x[i,def] = 1*/
1415 12099142 : GEN xj = gel(x,j), u = gel(xj,i);
1416 12099142 : if (!signe(u)) continue;
1417 :
1418 7666106 : ZC_lincomb1_inplace(xj, gel(x,def), negi(u));
1419 39083312 : for (k = 1; k < i; k++) gel(xj,k) = modii(gel(xj,k), p);
1420 : }
1421 5683839 : if (gc_needed(av,2))
1422 : {
1423 0 : if (DEBUGMEM>1) pari_warn(warnmem,"FpM_echelon. i=%ld",i);
1424 0 : x = gc_GEN(av, x);
1425 : }
1426 5683839 : def--;
1427 : }
1428 : /* rank = iP - 1 */
1429 3366410 : setlg(P, iP); vecsmall_sort(P);
1430 3366410 : if (co > iP) x += co - iP;
1431 3366410 : x[0] = evaltyp(t_MAT) | _evallg(iP);
1432 3366410 : return x;
1433 : }
1434 : /* given x square of maximal rank with 1 or p on diagonal from hnfmodid
1435 : * (=> a column containing p has its other entries at 0 ), return the HNF */
1436 : static GEN
1437 3366410 : FpM_hnfend(pari_sp av, GEN x, GEN p)
1438 : {
1439 3366410 : long i, l = lgcols(x);
1440 12687042 : for (i = l-1; i > 0; i--)
1441 : {
1442 9320632 : GEN diag = gcoeff(x,i,i);
1443 : long j;
1444 9320632 : if (is_pm1(diag))
1445 11464096 : for (j = i+1; j < l; j++)
1446 : {
1447 5780257 : GEN xj = gel(x,j), b = gel(xj,i);
1448 : long k;
1449 5780257 : if (!signe(b)) continue;
1450 4067225 : ZC_lincomb1_inplace(xj, gel(x,i), negi(b));
1451 19687805 : for (k=1; k<i; k++)
1452 15620580 : if (lgefint(gel(xj,k)) > 3) gel(xj,k) = remii(gel(xj,k), p);
1453 : }
1454 : else
1455 10302396 : for (j = i+1; j < l; j++) gcoeff(x,i,j) = modii(gcoeff(x,i,j), p);
1456 9320632 : if (gc_needed(av,2))
1457 : {
1458 0 : if (DEBUGMEM>1) pari_warn(warnmem,"FpM_hnfend. i=%ld",i);
1459 0 : x = gc_GEN(av, x);
1460 : }
1461 : }
1462 3366410 : return x;
1463 : }
1464 : GEN
1465 3366410 : ZM_hnfmodprime(GEN x, GEN p)
1466 : {
1467 3366410 : pari_sp av = avma;
1468 : GEN P, y;
1469 : long l, lP, i;
1470 3366410 : if (lg(x) == 1) return cgetg(1, t_MAT);
1471 3366410 : l = lgcols(x);
1472 3366410 : x = FpM_echelon(x, &P, p);
1473 3366410 : lP = lg(P); /* rank = lP-1 */
1474 3366410 : if (lP == l) { set_avma(av); return matid(l-1); }
1475 3366410 : y = scalarmat_shallow(p, l-1);
1476 9050249 : for (i = 1; i < lP; i++) gel(y,P[i]) = gel(x,i);
1477 3366410 : return gc_GEN(av, FpM_hnfend(av,y,p));
1478 : }
1479 :
1480 : static GEN
1481 438580 : allhnfmod(GEN x, GEN dm, int flag)
1482 : {
1483 438580 : if (typ(x)!=t_MAT) pari_err_TYPE("allhnfmod",x);
1484 438580 : RgM_check_ZM(x, "allhnfmod");
1485 438573 : if (isintzero(dm)) return ZM_hnf(x);
1486 438573 : return ZM_hnfmodall(x, dm, flag);
1487 : }
1488 : GEN
1489 14 : hnfmod(GEN x, GEN d)
1490 : {
1491 14 : if (typ(d) != t_INT) pari_err_TYPE("mathnfmod",d);
1492 14 : return allhnfmod(x, d, 0);
1493 : }
1494 : GEN
1495 438566 : hnfmodid(GEN x, GEN d)
1496 : {
1497 438566 : switch(typ(d))
1498 : {
1499 317192 : case t_INT: break;
1500 121374 : case t_VEC: case t_COL:
1501 121374 : if (RgV_is_ZV(d)) break;
1502 0 : default: pari_err_TYPE("mathnfmodid",d);
1503 : }
1504 438566 : return allhnfmod(x, d, hnf_MODID);
1505 : }
1506 :
1507 : /* M a ZM in HNF. Normalize with *centered* residues */
1508 : GEN
1509 25387 : ZM_hnfcenter(GEN M)
1510 : {
1511 25387 : long i, j, k, N = lg(M)-1;
1512 25387 : pari_sp av = avma;
1513 :
1514 89869 : for (j=N-1; j>0; j--) /* skip last line */
1515 : {
1516 64482 : GEN Mj = gel(M,j), a = gel(Mj,j);
1517 215828 : for (k = j+1; k <= N; k++)
1518 : {
1519 151346 : GEN Mk = gel(M,k), q = diviiround(gel(Mk,j), a);
1520 151346 : long s = signe(q);
1521 151346 : if (!s) continue;
1522 69265 : if (is_pm1(q))
1523 : {
1524 32863 : if (s < 0)
1525 9221 : for (i = 1; i <= j; i++) gel(Mk,i) = addii(gel(Mk,i), gel(Mj,i));
1526 : else
1527 95617 : for (i = 1; i <= j; i++) gel(Mk,i) = subii(gel(Mk,i), gel(Mj,i));
1528 : }
1529 : else
1530 204575 : for (i = 1; i <= j; i++) gel(Mk,i) = subii(gel(Mk,i), mulii(q,gel(Mj,i)));
1531 69265 : if (gc_needed(av,1))
1532 : {
1533 0 : if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnfcenter, j = %ld",j);
1534 0 : M = gc_GEN(av, M);
1535 : }
1536 : }
1537 : }
1538 25387 : return M;
1539 : }
1540 :
1541 : /***********************************************************************/
1542 : /* */
1543 : /* HNFLLL (Havas, Majewski, Mathews) */
1544 : /* */
1545 : /***********************************************************************/
1546 :
1547 : static void
1548 2997278 : Minus(long j, GEN lambda)
1549 : {
1550 2997278 : long k, n = lg(lambda);
1551 :
1552 15683027 : for (k=1 ; k<j; k++) togglesign_safe(&gcoeff(lambda,k,j));
1553 22849526 : for (k=j+1; k<n; k++) togglesign_safe(&gcoeff(lambda,j,k));
1554 2997278 : }
1555 :
1556 : /* index of first nonzero entry */
1557 : static long
1558 129171008 : findi(GEN M)
1559 : {
1560 129171008 : long i, n = lg(M);
1561 2701401932 : for (i=1; i<n; i++)
1562 2661122121 : if (signe(gel(M,i))) return i;
1563 40279811 : return 0;
1564 : }
1565 :
1566 : static long
1567 129171008 : findi_normalize(GEN Aj, GEN B, long j, GEN lambda)
1568 : {
1569 129171008 : long r = findi(Aj);
1570 129171008 : if (r && signe(gel(Aj,r)) < 0)
1571 : {
1572 2997278 : ZV_togglesign(Aj); if (B) ZV_togglesign(gel(B,j));
1573 2997278 : Minus(j,lambda);
1574 : }
1575 129171008 : return r;
1576 : }
1577 :
1578 : static void
1579 7929808 : subzi(GEN *a, GEN b)
1580 : {
1581 7929808 : pari_sp av = avma;
1582 7929808 : b = subii(*a, b);
1583 7929808 : if (lgefint(b)<=lg(*a) && isonstack(*a)) { affii(b,*a); set_avma(av); }
1584 1559566 : else *a = b;
1585 7929808 : }
1586 :
1587 : static void
1588 283906255 : addzi(GEN *a, GEN b)
1589 : {
1590 283906255 : pari_sp av = avma;
1591 283906255 : b = addii(*a, b);
1592 283906255 : if (lgefint(b)<=lg(*a) && isonstack(*a)) { affii(b,*a); set_avma(av); }
1593 23660044 : else *a = b;
1594 283906255 : }
1595 :
1596 : /* x + y*z */
1597 : static void
1598 592164488 : addmulzii(GEN *x, GEN y, GEN z)
1599 : {
1600 592164488 : long lx = lgefint(*x);
1601 : pari_sp av;
1602 : GEN t;
1603 592164488 : long ly = lgefint(y), lz;
1604 592164488 : if (ly == 2) return;
1605 280462450 : lz = lgefint(z);
1606 280462450 : av = avma; (void)new_chunk(lx+ly+lz); /*HACK*/
1607 280462450 : t = mulii(z, y);
1608 280462450 : set_avma(av); return addzi(x,t);
1609 : }
1610 :
1611 : static void
1612 30080471 : ZC_lincomb1z_i(GEN X, GEN Y, GEN v, long n)
1613 : {
1614 : long i;
1615 572972709 : for (i = n; i; i--) addmulzii(&gel(X,i), gel(Y,i), v);
1616 30080471 : }
1617 : /* X <- X + v Y (elementary col operation) */
1618 : static void
1619 30080471 : ZC_lincomb1z(GEN X, GEN Y, GEN v)
1620 : {
1621 30080471 : if (lgefint(v) != 2) return ZC_lincomb1z_i(X, Y, v, lg(X)-1);
1622 : }
1623 :
1624 : static void
1625 64584596 : reduce2(GEN A, GEN B, long k, long j, long *row0, long *row1, GEN lambda, GEN D)
1626 : {
1627 : GEN q;
1628 : long i;
1629 :
1630 64584596 : *row0 = findi_normalize(gel(A,j), B,j,lambda);
1631 64584596 : *row1 = findi_normalize(gel(A,k), B,k,lambda);
1632 64584596 : if (*row0)
1633 36006564 : q = truedivii(gcoeff(A,*row0,k), gcoeff(A,*row0,j));
1634 : else
1635 : {
1636 28578032 : pari_sp btop = avma;
1637 28578032 : int c = abscmpii(shifti(gcoeff(lambda,j,k), 1), gel(D,j));
1638 28578032 : set_avma(btop);
1639 28578032 : if (c > 0)
1640 10335054 : q = diviiround(gcoeff(lambda,j,k), gel(D,j));
1641 : else
1642 18242978 : return;
1643 : }
1644 :
1645 46341618 : if (signe(q))
1646 : {
1647 20208137 : GEN Lk = gel(lambda,k), Lj = gel(lambda,j);
1648 20208137 : togglesign_safe(&q);
1649 20208137 : if (*row0) ZC_lincomb1z(gel(A,k),gel(A,j),q);
1650 20208137 : if (B) ZC_lincomb1z(gel(B,k),gel(B,j),q);
1651 20208137 : addmulzii(&gel(Lk,j), q, gel(D,j));
1652 20208137 : if (is_pm1(q))
1653 : {
1654 7668906 : if (signe(q) > 0)
1655 : {
1656 7405312 : for (i=1; i<j; i++)
1657 4394031 : if (signe(gel(Lj,i))) addzi(&gel(Lk,i), gel(Lj,i));
1658 : }
1659 : else
1660 : {
1661 14162629 : for (i=1; i<j; i++)
1662 9505004 : if (signe(gel(Lj,i))) subzi(&gel(Lk,i), gel(Lj,i));
1663 : }
1664 : }
1665 : else
1666 : {
1667 49099343 : for (i=1; i<j; i++)
1668 36560112 : if (signe(gel(Lj,i))) addmulzii(&gel(Lk,i), q, gel(Lj,i));
1669 : }
1670 : }
1671 : }
1672 : static void
1673 54526249 : affii2_or_copy_gc(pari_sp av, GEN x, GEN *y, GEN x1, GEN *y1)
1674 : {
1675 54526249 : long l = lg(*y), l1 = lg(*y1);
1676 54526249 : if(x==*y && x1==*y1) return;
1677 54526249 : if (lgefint(x) <= l && isonstack(*y) && lgefint(x1) <= l1 && isonstack(*y1))
1678 : {
1679 52979086 : affii(x,*y);
1680 52979086 : affii(x1,*y1);
1681 52979086 : set_avma(av);
1682 : }
1683 1547163 : else if (lgefint(x) <= l && isonstack(*y))
1684 : {
1685 1109234 : affii(x,*y);
1686 1109234 : *y1 = gc_GEN(av, x1);
1687 : }
1688 437929 : else if (lgefint(x1) <= l1 && isonstack(*y1))
1689 : {
1690 385025 : affii(x1,*y1);
1691 385025 : *y = gc_GEN(av, x);
1692 : }
1693 : else
1694 : {
1695 52904 : *y = cgeti(2*lg(x));
1696 52904 : *y1 = cgeti(2*lg(x1));
1697 52904 : affii(x,*y);
1698 52904 : affii(x1,*y1);
1699 52904 : (void)gc_all(av,2,y,y1);
1700 : }
1701 : }
1702 : static void
1703 11659440 : hnfswap(GEN A, GEN B, long k, GEN lambda, GEN D)
1704 : {
1705 : pari_sp av;
1706 11659440 : GEN t, p1, p2, Lk = gel(lambda,k);
1707 11659440 : long i,j,n = lg(A);
1708 :
1709 11659440 : swap(gel(A,k), gel(A,k-1));
1710 11659440 : if (B) swap(gel(B,k), gel(B,k-1));
1711 55413320 : for (j=k-2; j; j--) swap(gcoeff(lambda,j,k-1), gel(Lk,j));
1712 125360197 : for (i=k+1; i<n; i++)
1713 : {
1714 113700757 : pari_sp btop = avma;
1715 113700757 : GEN Li = gel(lambda,i);
1716 113700757 : if (signe(gel(Li,k-1))==0 && signe(gel(Li,k))==0) continue;
1717 54526249 : p1 = mulii(gel(Li,k-1), gel(D,k));
1718 54526249 : p2 = mulii(gel(Li,k), gel(Lk,k-1));
1719 54526249 : t = subii(p1,p2);
1720 54526249 : p1 = mulii(gel(Li,k), gel(D,k-2));
1721 54526249 : p2 = mulii(gel(Li,k-1), gel(Lk,k-1));
1722 54526249 : affii2_or_copy_gc(btop, diviiexact(addii(p1,p2), gel(D,k-1)),
1723 54526249 : &gel(Li,k-1), diviiexact(t, gel(D,k-1)), &gel(Li,k));
1724 : }
1725 11659440 : av = avma;
1726 11659440 : p1 = mulii(gel(D,k-2), gel(D,k));
1727 11659440 : p2 = sqri(gel(Lk,k-1));
1728 11659440 : gc_INT_affii(av, diviiexact(addii(p1,p2), gel(D,k-1)), &gel(D,k-1));
1729 11659440 : }
1730 :
1731 : /* reverse row order in matrix A, IN PLACE */
1732 : static GEN
1733 618720 : reverse_rows(GEN A)
1734 : {
1735 618720 : long i, j, h, n = lg(A);
1736 618720 : if (n == 1) return A;
1737 618720 : h = lgcols(A);
1738 5300866 : for (j=1; j<n; j++)
1739 : {
1740 4682146 : GEN c = gel(A,j);
1741 : /* start at (h-1) >>1 : if h = 2i even, no need to swap c[i] and itself */
1742 13081868 : for (i=(h-1)>>1; i; i--) swap(gel(c,i), gel(c,h-i));
1743 : }
1744 618720 : return A;
1745 : }
1746 : /* decide whether to swap */
1747 : static int
1748 6205727 : must_swap(long k, GEN lambda, GEN D)
1749 : {
1750 6205727 : pari_sp av = avma;
1751 6205727 : GEN z = addii(mulii(gel(D,k-2),gel(D,k)), sqri(gcoeff(lambda,k-1,k)));
1752 6205727 : return gc_bool(av, cmpii(z, sqri(gel(D,k-1))) < 0);
1753 : }
1754 :
1755 : GEN
1756 309360 : ZM_hnflll(GEN A, GEN *ptB, int remove)
1757 : {
1758 309360 : pari_sp av = avma;
1759 : long n, k, kmax;
1760 : GEN B, lambda, D;
1761 :
1762 309360 : n = lg(A);
1763 309360 : A = reverse_rows(ZM_copy(A)); /* ZM_copy for in place findi_normalize() */
1764 309360 : B = ptB? matid(n-1): NULL;
1765 309360 : D = const_vec(n, gen_1) + 1;
1766 309360 : lambda = zeromatcopy(n-1,n-1);
1767 309360 : k = kmax = 2;
1768 21672565 : while (k < n)
1769 : {
1770 : long row0, row1;
1771 : int do_swap;
1772 21363205 : reduce2(A,B,k,k-1,&row0,&row1,lambda,D);
1773 21363205 : if (row0) do_swap = (!row1 || row0 <= row1);
1774 7562050 : else if (row1) do_swap = 0;
1775 5965068 : else do_swap = must_swap(k,lambda,D);
1776 21363205 : if (do_swap)
1777 : {
1778 11401546 : hnfswap(A,B,k,lambda,D);
1779 11401546 : if (k > 2) k--;
1780 : }
1781 : else
1782 : {
1783 : long i;
1784 53183050 : for (i=k-2; i; i--)
1785 : {
1786 : long row0, row1;
1787 43221391 : reduce2(A,B,k,i,&row0,&row1,lambda,D);
1788 : }
1789 9961659 : if (++k > kmax) kmax = k;
1790 : }
1791 21363205 : if (gc_needed(av,3))
1792 : {
1793 0 : GEN b = D-1;
1794 0 : if (DEBUGMEM>1) pari_warn(warnmem,"hnflll, kmax = %ld / %ld",kmax,n-1);
1795 0 : (void)gc_all(av, B? 4: 3, &A, &lambda, &b, &B);
1796 0 : if (_gc_needed(av,1)) paristack_resize(0); /* avoid desperation GC */
1797 0 : D = b+1;
1798 : }
1799 : }
1800 : /* handle trivial case: return negative diag coefficient otherwise */
1801 309360 : if (n == 2) (void)findi_normalize(gel(A,1), B,1,lambda);
1802 309360 : A = reverse_rows(A);
1803 309360 : if (remove)
1804 : {
1805 : long i;
1806 7525 : for (i = 1; i < n; i++)
1807 7301 : if (!ZV_equal0(gel(A,i))) break;
1808 2275 : remove_0cols(i-1, &A, &B, remove);
1809 : }
1810 309360 : (void)gc_all(av, B? 2: 1, &A, &B);
1811 309360 : if (B) *ptB = B;
1812 309360 : return A;
1813 : }
1814 :
1815 : GEN
1816 7 : hnflll(GEN x)
1817 : {
1818 7 : GEN z = cgetg(3, t_VEC);
1819 7 : gel(z,1) = ZM_hnflll(x, &gel(z,2), 1);
1820 7 : return z;
1821 : }
1822 :
1823 : /* Variation on HNFLLL: Extended GCD */
1824 :
1825 : static void
1826 1016504 : reduce1(GEN A, GEN B, long k, long j, GEN lambda, GEN D)
1827 : {
1828 : GEN q;
1829 : long i;
1830 :
1831 1016504 : if (signe(gel(A,j)))
1832 182669 : q = diviiround(gel(A,k),gel(A,j));
1833 833835 : else if (abscmpii(shifti(gcoeff(lambda,j,k), 1), gel(D,j)) > 0)
1834 114994 : q = diviiround(gcoeff(lambda,j,k), gel(D,j));
1835 : else
1836 718841 : return;
1837 :
1838 297663 : if (signe(q))
1839 : {
1840 260161 : GEN Lk = gel(lambda,k), Lj = gel(lambda,j);
1841 260161 : togglesign_safe(&q);
1842 260161 : gel(A,k) = addmulii(gel(A,k), q, gel(A,j));
1843 260161 : ZC_lincomb1_inplace(gel(B,k),gel(B,j),q);
1844 260161 : gel(Lk,j) = addmulii(gel(Lk,j), q, gel(D,j));
1845 491772 : for (i=1; i<j; i++)
1846 231611 : if (signe(gel(Lj,i))) gel(Lk,i) = addmulii(gel(Lk,i), q, gel(Lj,i));
1847 : }
1848 : }
1849 :
1850 : static GEN
1851 111734 : ZV_gcdext_i(GEN A)
1852 : {
1853 111734 : long k, n = lg(A);
1854 : GEN B, lambda, D;
1855 :
1856 111734 : if (n == 1) retmkvec2(gen_1, cgetg(1,t_MAT));
1857 111734 : A = leafcopy(A);
1858 111734 : B = matid(n-1);
1859 111734 : lambda = zeromatcopy(n-1,n-1);
1860 111734 : D = const_vec(n, gen_1) + 1;
1861 111734 : k = 2;
1862 735891 : while (k < n)
1863 : {
1864 : int do_swap;
1865 :
1866 624157 : reduce1(A,B,k,k-1,lambda,D);
1867 624157 : if (signe(gel(A,k-1))) do_swap = 1;
1868 441488 : else if (signe(gel(A,k))) do_swap = 0;
1869 240659 : else do_swap = must_swap(k,lambda,D);
1870 624157 : if (do_swap)
1871 : {
1872 257894 : hnfswap(A,B,k,lambda,D);
1873 257894 : if (k > 2) k--;
1874 : }
1875 : else
1876 : {
1877 : long i;
1878 758610 : for (i=k-2; i; i--) reduce1(A,B,k,i,lambda,D);
1879 366263 : k++;
1880 : }
1881 : }
1882 111734 : if (signe(gel(A,n-1)) < 0)
1883 : {
1884 14724 : gel(A,n-1) = negi(gel(A,n-1));
1885 14724 : ZV_togglesign(gel(B,n-1));
1886 : }
1887 111734 : return mkvec2(gel(A,n-1), B);
1888 : }
1889 : GEN
1890 110579 : ZV_extgcd(GEN A)
1891 : {
1892 110579 : pari_sp av = avma;
1893 110579 : return gc_GEN(av, ZV_gcdext_i(A));
1894 : }
1895 : /* as ZV_extgcd, transforming the gcd into a t_MAT, for mathnf0 */
1896 : static GEN
1897 1162 : ZV_hnfgcdext(GEN A)
1898 : {
1899 1162 : pari_sp av = avma;
1900 : GEN z;
1901 1162 : if (lg(A) == 1) retmkvec2(cgetg(1,t_MAT),cgetg(1,t_MAT));
1902 1155 : z = ZV_gcdext_i(A);
1903 1155 : gel(z,1) = mkmat(mkcol(gel(z,1)));
1904 1155 : return gc_GEN(av, z);
1905 : }
1906 :
1907 : /* HNF with permutation. */
1908 : GEN
1909 385 : ZM_hnfperm(GEN A, GEN *ptU, GEN *ptperm)
1910 : {
1911 : GEN U, c, l, perm, d, p, q, b;
1912 385 : pari_sp av = avma, av1;
1913 : long r, t, i, j, j1, k, m, n;
1914 :
1915 385 : n = lg(A)-1;
1916 385 : if (!n)
1917 : {
1918 7 : if (ptU) *ptU = cgetg(1,t_MAT);
1919 7 : if (ptperm) *ptperm = cgetg(1,t_VEC);
1920 7 : return cgetg(1, t_MAT);
1921 : }
1922 378 : m = nbrows(A);
1923 378 : c = zero_zv(m);
1924 378 : l = zero_zv(n);
1925 378 : perm = cgetg(m+1, t_VECSMALL);
1926 378 : av1 = avma;
1927 378 : A = RgM_shallowcopy(A);
1928 378 : U = ptU? matid(n): NULL;
1929 : /* U base change matrix : A0*U = A all along */
1930 1736 : for (r=0, k=1; k <= n; k++)
1931 : {
1932 4018 : for (j=1; j<k; j++)
1933 : {
1934 2660 : if (!l[j]) continue;
1935 2506 : t=l[j]; b=gcoeff(A,t,k);
1936 2506 : if (!signe(b)) continue;
1937 :
1938 560 : ZC_elem(b,gcoeff(A,t,j), A,U,k,j);
1939 560 : d = gcoeff(A,t,j);
1940 560 : if (signe(d) < 0)
1941 : {
1942 0 : ZV_neg_inplace(gel(A,j));
1943 0 : if (U) ZV_togglesign(gel(U,j));
1944 0 : d = gcoeff(A,t,j);
1945 : }
1946 1456 : for (j1=1; j1<j; j1++)
1947 : {
1948 896 : if (!l[j1]) continue;
1949 875 : q = truedivii(gcoeff(A,t,j1),d);
1950 875 : if (!signe(q)) continue;
1951 :
1952 329 : togglesign(q);
1953 329 : ZC_lincomb1_inplace(gel(A,j1), gel(A,j), q);
1954 329 : if (U) ZC_lincomb1_inplace(gel(U,j1), gel(U,j), q);
1955 : }
1956 : }
1957 4725 : t = m; while (t && (c[t] || !signe(gcoeff(A,t,k)))) t--;
1958 1358 : if (t)
1959 : {
1960 1239 : p = gcoeff(A,t,k);
1961 20041 : for (i=t-1; i; i--)
1962 : {
1963 18802 : q = gcoeff(A,i,k);
1964 18802 : if (signe(q) && abscmpii(p,q) > 0) { p = q; t = i; }
1965 : }
1966 1239 : perm[++r] = l[k] = t; c[t] = k;
1967 1239 : if (signe(p) < 0)
1968 : {
1969 317 : ZV_neg_inplace(gel(A,k));
1970 317 : if (U) ZV_togglesign(gel(U,k));
1971 317 : p = gcoeff(A,t,k);
1972 : }
1973 : /* p > 0 */
1974 3535 : for (j=1; j<k; j++)
1975 : {
1976 2296 : if (!l[j]) continue;
1977 2261 : q = truedivii(gcoeff(A,t,j),p);
1978 2261 : if (!signe(q)) continue;
1979 :
1980 224 : togglesign(q);
1981 224 : ZC_lincomb1_inplace(gel(A,j), gel(A,k), q);
1982 224 : if (U) ZC_lincomb1_inplace(gel(U,j), gel(U,k), q);
1983 : }
1984 : }
1985 1358 : if (gc_needed(av1,1))
1986 : {
1987 0 : if (DEBUGMEM>1) pari_warn(warnmem,"hnfperm, k=%ld",k);
1988 0 : (void)gc_all(av1, U? 2: 1, &A, &U);
1989 : }
1990 : }
1991 378 : if (r < m)
1992 : {
1993 5145 : for (i=1,k=r; i<=m; i++)
1994 4830 : if (!c[i]) perm[++k] = i;
1995 : }
1996 :
1997 : /* We have A0*U=A, U in Gl(n,Z)
1998 : * basis for Im(A): columns of A s.t l[j]>0 (r cols)
1999 : * basis for Ker(A): columns of U s.t l[j]=0 (n-r cols) */
2000 378 : p = cgetg(r+1,t_MAT);
2001 2828 : for (i=1; i<=m/2; i++) lswap(perm[i], perm[m+1-i]);
2002 378 : if (U)
2003 : {
2004 84 : GEN u = cgetg(n+1,t_MAT);
2005 378 : for (t=1,k=r,j=1; j<=n; j++)
2006 294 : if (l[j])
2007 : {
2008 182 : u[k + n-r] = U[j];
2009 182 : gel(p,k--) = vecpermute(gel(A,j), perm);
2010 : }
2011 : else
2012 112 : u[t++] = U[j];
2013 84 : *ptU = u;
2014 84 : if (ptperm) *ptperm = perm;
2015 84 : (void)gc_all(av, ptperm? 3: 2, &p, ptU, ptperm);
2016 : }
2017 : else
2018 : {
2019 1358 : for (k=r,j=1; j<=n; j++)
2020 1064 : if (l[j]) gel(p,k--) = vecpermute(gel(A,j), perm);
2021 294 : if (ptperm) *ptperm = perm;
2022 294 : (void)gc_all(av, ptperm? 2: 1, &p, ptperm);
2023 : }
2024 378 : return p;
2025 : }
2026 :
2027 : GEN
2028 238 : ZM_hnf_knapsack(GEN x)
2029 : {
2030 238 : GEN t, perm, H = ZM_hnfperm(x,NULL,&perm);
2031 238 : long i,j, l = lg(H), h = lgcols(H);
2032 3311 : for (i=1; i<h; i++)
2033 : {
2034 3164 : int fl = 0;
2035 17101 : for (j=1; j<l; j++)
2036 : {
2037 14028 : t = gcoeff(H,i,j);
2038 14028 : if (signe(t))
2039 : {
2040 3213 : if (!is_pm1(t) || fl) return NULL;
2041 3122 : fl = 1;
2042 : }
2043 : }
2044 : }
2045 147 : return rowpermute(H, perm_inv(perm));
2046 : }
2047 :
2048 : GEN
2049 14 : hnfperm(GEN A)
2050 : {
2051 14 : GEN y = cgetg(4, t_VEC);
2052 14 : gel(y,1) = ZM_hnfperm(A, &gel(y,2), &gel(y,3));
2053 14 : return y;
2054 : }
2055 :
2056 : /* Hermite Normal Form, with base change matrix if ptB != NULL.
2057 : * If 'remove' = 1, remove 0 columns (do NOT update *ptB accordingly)
2058 : * If 'remove' = 2, remove 0 columns and update *ptB accordingly */
2059 : GEN
2060 752457 : ZM_hnfall_i(GEN A, GEN *ptB, long remove)
2061 : {
2062 : pari_sp av;
2063 : long m, n, r, i, j, k, li;
2064 : GEN B, c, h, a;
2065 :
2066 752457 : RgM_dimensions(A, &m,&n);
2067 752457 : if (!n)
2068 : {
2069 7 : if (ptB) *ptB = cgetg(1,t_MAT);
2070 7 : return cgetg(1,t_MAT);
2071 : }
2072 752450 : c = zero_zv(m);
2073 752450 : h = const_vecsmall(n, m);
2074 752450 : av = avma;
2075 752450 : A = RgM_shallowcopy(A);
2076 752450 : B = ptB? matid(n): NULL;
2077 752450 : r = n+1;
2078 2151216 : for (li=m; li; li--)
2079 : {
2080 3086654 : for (j=1; j<r; j++)
2081 : {
2082 3885839 : for (i=h[j]; i>li; i--)
2083 : {
2084 799507 : a = gcoeff(A,i,j);
2085 799507 : k = c[i];
2086 : /* zero a = Aij using Aik */
2087 799507 : if (signe(a)) ZC_elem(a,gcoeff(A,i,k), A,B,j,k);
2088 799507 : ZM_reduce(A,B, i,k); /* ensure reduced entries even if a = 0 */
2089 : }
2090 3086332 : if (gc_needed(av,1) && (j & 0x7f) == 0)
2091 : {
2092 0 : if (DEBUGMEM>1)
2093 0 : pari_warn(warnmem,"ZM_hnfall[1], li = %ld, j = %ld", li, j);
2094 0 : (void)gc_all(av, B? 2: 1, &A, &B);
2095 : }
2096 3086332 : if (signe( gcoeff(A,li,j) )) break;
2097 1687888 : h[j] = li-1;
2098 : }
2099 1398766 : if (j == r) continue;
2100 1398444 : r--;
2101 1398444 : if (j < r) /* A[j] != 0 */
2102 : {
2103 906134 : swap(gel(A,j), gel(A,r));
2104 906134 : if (B) swap(gel(B,j), gel(B,r));
2105 906134 : h[j] = h[r]; h[r] = li; c[li] = r;
2106 : }
2107 1398444 : if (signe(gcoeff(A,li,r)) < 0)
2108 : {
2109 340538 : ZV_neg_inplace(gel(A,r));
2110 340538 : if (B) ZV_togglesign(gel(B,r));
2111 : }
2112 1398444 : ZM_reduce(A,B, li,r);
2113 1398444 : if (gc_needed(av,1))
2114 : {
2115 0 : if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnfall[2], li = %ld", li);
2116 0 : (void)gc_all(av, B? 2: 1, &A, &B);
2117 : }
2118 : }
2119 :
2120 752450 : if (DEBUGLEVEL>5) err_printf("\nhnfall, final phase: ");
2121 752450 : r--; /* first r cols are in the image the n-r (independent) last ones */
2122 2568149 : for (j=1; j<=r; j++)
2123 : {
2124 3969802 : for (i=h[j]; i; i--)
2125 : {
2126 2154103 : a = gcoeff(A,i,j);
2127 2154103 : k = c[i];
2128 2154103 : if (signe(a)) ZC_elem(a,gcoeff(A,i,k), A,B, j,k);
2129 2154103 : ZM_reduce(A,B, i,k); /* ensure reduced entries, even if a = 0 */
2130 : }
2131 1815699 : if (gc_needed(av,1) && (j & 0x7f) == 0)
2132 : {
2133 0 : if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnfall[3], j = %ld", j);
2134 0 : (void)gc_all(av, B? 2: 1, &A, &B);
2135 : }
2136 : }
2137 752450 : if (DEBUGLEVEL>5) err_printf("\n");
2138 752450 : if (remove) remove_0cols(r, &A, &B, remove);
2139 752450 : if (ptB) *ptB = B;
2140 752450 : return A;
2141 : }
2142 : GEN
2143 29170 : ZM_hnfall(GEN A, GEN *ptB, long remove)
2144 : {
2145 29170 : pari_sp av = avma;
2146 29170 : A = ZM_hnfall_i(A, ptB, remove);
2147 29170 : return gc_all(av, ptB? 2: 1, &A, ptB);
2148 : }
2149 :
2150 : GEN
2151 14 : hnfall(GEN x)
2152 : {
2153 14 : GEN z = cgetg(3, t_VEC);
2154 14 : gel(z,1) = ZM_hnfall(x, (GEN*)(z+2), 1);
2155 14 : return z;
2156 : }
2157 : GEN
2158 14 : hnf(GEN x) { return mathnf0(x,0); }
2159 :
2160 : /* C = A^(-1)t where A and C are integral, A is upper triangular, t t_INT */
2161 : GEN
2162 275331 : hnf_invscale(GEN A, GEN t)
2163 : {
2164 275331 : long n = lg(A)-1, i,j,k;
2165 275331 : GEN m, c = cgetg(n+1,t_MAT);
2166 :
2167 275331 : if (!n) return c;
2168 996275 : for (k=1; k<=n; k++)
2169 : { /* cf hnf_divscale with B = id, thus b = e_k */
2170 720944 : GEN u = cgetg(n+1, t_COL);
2171 720944 : pari_sp av = avma;
2172 720944 : gel(c,k) = u;
2173 720944 : gel(u,n) = k == n? gc_INT(av, diviiexact(t, gcoeff(A,n,n))): gen_0;
2174 2341038 : for (i=n-1; i>0; i--)
2175 : {
2176 1620094 : av = avma; m = i == k? t: gen_0;
2177 5990117 : for (j=i+1; j<=n; j++) m = subii(m, mulii(gcoeff(A,i,j),gel(u,j)));
2178 1620094 : gel(u,i) = gc_INT(av, diviiexact(m, gcoeff(A,i,i)));
2179 : }
2180 : }
2181 275331 : return c;
2182 : }
2183 :
2184 : /* C = A^(-1)(tB) where A, B, C are integral, A is upper triangular, t t_INT */
2185 : GEN
2186 197801 : hnf_divscale(GEN A, GEN B, GEN t)
2187 : {
2188 197801 : long n = lg(A)-1, i,j,k;
2189 197801 : GEN m, c = cgetg(n+1,t_MAT);
2190 :
2191 197801 : if (!n) return c;
2192 879139 : for (k=1; k<=n; k++)
2193 : {
2194 681338 : GEN u = cgetg(n+1, t_COL), b = gel(B,k);
2195 681338 : pari_sp av = avma;
2196 681338 : gel(c,k) = u; m = mulii(gel(b,n),t);
2197 681338 : gel(u,n) = gc_INT(av, diviiexact(m, gcoeff(A,n,n)));
2198 3947530 : for (i=n-1; i>0; i--)
2199 : {
2200 3266192 : av = avma; m = mulii(gel(b,i),t);
2201 30942472 : for (j=i+1; j<=n; j++) m = subii(m, mulii(gcoeff(A,i,j),gel(u,j)));
2202 3266192 : gel(u,i) = gc_INT(av, diviiexact(m, gcoeff(A,i,i)));
2203 : }
2204 : }
2205 197801 : return c;
2206 : }
2207 :
2208 : /* A, B integral upper HNF. A^(-1) B integral ? */
2209 : int
2210 133 : hnfdivide(GEN A, GEN B)
2211 : {
2212 133 : pari_sp av = avma;
2213 133 : long n = lg(A)-1, i,j,k;
2214 : GEN u, b, m, r;
2215 :
2216 133 : if (!n) return 1;
2217 133 : if (lg(B)-1 != n) pari_err_DIM("hnfdivide");
2218 133 : u = cgetg(n+1, t_COL);
2219 483 : for (k=1; k<=n; k++)
2220 : {
2221 364 : b = gel(B,k);
2222 364 : m = gel(b,k);
2223 364 : gel(u,k) = dvmdii(m, gcoeff(A,k,k), &r);
2224 364 : if (r != gen_0) return gc_long(av, 0);
2225 826 : for (i=k-1; i>0; i--)
2226 : {
2227 476 : m = gel(b,i);
2228 1337 : for (j=i+1; j<=k; j++) m = subii(m, mulii(gcoeff(A,i,j),gel(u,j)));
2229 476 : m = dvmdii(m, gcoeff(A,i,i), &r);
2230 476 : if (r != gen_0) return gc_long(av, 0);
2231 476 : gel(u,i) = m;
2232 : }
2233 : }
2234 119 : return gc_long(av, 1);
2235 : }
2236 :
2237 : /* A upper HNF, b integral vector. Return A^(-1) b if integral,
2238 : * NULL otherwise. Assume #A[,1] = #b. */
2239 : GEN
2240 547578 : hnf_invimage(GEN A, GEN b)
2241 : {
2242 547578 : pari_sp av = avma;
2243 547578 : long n = lg(A)-1, m, i, k;
2244 : GEN u, r;
2245 :
2246 547578 : if (!n) return lg(b)==1? cgetg(1,t_COL):NULL;
2247 547536 : m = nbrows(A); /* m >= n */
2248 547536 : u = cgetg(n+1, t_COL);
2249 1505933 : for (i = n, k = m; k > 0; k--)
2250 : {
2251 1505933 : pari_sp av2 = avma;
2252 : long j;
2253 1505933 : GEN t = gel(b,k), Aki = gcoeff(A,k,i);
2254 1505933 : if (typ(t) != t_INT) pari_err_TYPE("hnf_invimage",t);
2255 4110428 : for (j=i+1; j<=n; j++) t = subii(t, mulii(gcoeff(A,k,j),gel(u,j)));
2256 1505933 : if (!signe(Aki))
2257 : {
2258 0 : if (signe(t)) return gc_NULL(av);
2259 0 : set_avma(av2); gel(u,i) = gen_0; continue;
2260 : }
2261 1505933 : t = dvmdii(t, Aki, &r);
2262 1505933 : if (r != gen_0) return gc_NULL(av);
2263 1164365 : gel(u,i) = gc_INT(av2, t);
2264 1164365 : if (--i == 0) break;
2265 : }
2266 : /* If there is a solution, it must be u. Check remaining equations */
2267 411936 : for (; k > 0; k--)
2268 : {
2269 205968 : pari_sp av2 = avma;
2270 : long j;
2271 205968 : GEN t = gel(b,k);
2272 205968 : if (typ(t) != t_INT) pari_err_TYPE("hnf_invimage",t);
2273 924952 : for (j=1; j<=n; j++) t = subii(t, mulii(gcoeff(A,k,j),gel(u,j)));
2274 205968 : if (signe(t)) return gc_NULL(av);
2275 205968 : set_avma(av2);
2276 : }
2277 205968 : return u;
2278 : }
2279 :
2280 : /* A upper HNF, B integral matrix or column. Return A^(-1) B if integral,
2281 : * NULL otherwise */
2282 : GEN
2283 461209 : hnf_solve(GEN A, GEN B)
2284 : {
2285 : pari_sp av;
2286 : long i, l;
2287 : GEN C;
2288 :
2289 461209 : if (typ(B) == t_COL) return hnf_invimage(A, B);
2290 357175 : av = avma; C = cgetg_copy(B, &l);
2291 484645 : for (i = 1; i < l; i++)
2292 : {
2293 372113 : GEN c = hnf_invimage(A, gel(B,i));
2294 372113 : if (!c) return gc_NULL(av);
2295 127470 : gel(C,i) = c;
2296 : }
2297 112532 : return C;
2298 : }
2299 :
2300 : /***************************************************************/
2301 : /** **/
2302 : /** SMITH NORMAL FORM REDUCTION **/
2303 : /** **/
2304 : /***************************************************************/
2305 :
2306 : static GEN
2307 0 : trivsmith(long all)
2308 : {
2309 : GEN z;
2310 0 : if (!all) return cgetg(1,t_VEC);
2311 0 : z=cgetg(4,t_VEC);
2312 0 : gel(z,1) = cgetg(1,t_MAT);
2313 0 : gel(z,2) = cgetg(1,t_MAT);
2314 0 : gel(z,3) = cgetg(1,t_MAT); return z;
2315 : }
2316 :
2317 : static void
2318 48 : snf_pile1(pari_sp av, GEN *x, GEN *U)
2319 48 : { (void)gc_all(av, *U? 2: 1, x, U); }
2320 : static void
2321 1124029 : snf_pile(pari_sp av, GEN *x, GEN *U, GEN *V)
2322 : {
2323 1124029 : int c = 1;
2324 1124029 : if (*U) c++;
2325 1124029 : if (*V) c++;
2326 1124029 : (void)gc_all(av, c, x, U, V);
2327 1124029 : }
2328 :
2329 : static GEN
2330 891006 : bezout_step(GEN *pa, GEN *pb, GEN *pu, GEN *pv)
2331 : {
2332 891006 : GEN a = *pa, b = *pb, d;
2333 891006 : if (absequalii(a,b))
2334 : {
2335 536476 : long sa = signe(a), sb = signe(b);
2336 536476 : *pv = gen_0;
2337 536476 : if (sb == sa) {
2338 523115 : *pa = *pb = gen_1;
2339 523115 : if (sa > 0) {*pu=gen_1; return a;} else {*pu=gen_m1; return absi(a);}
2340 : }
2341 13361 : if (sa > 0) { *pa = *pu = gen_1; *pb = gen_m1; return a; }
2342 2484 : *pa = *pu = gen_m1; *pb = gen_1; return b;
2343 : }
2344 354530 : d = bezout(a,b, pu,pv);
2345 354530 : *pa = diviiexact(a, d);
2346 354530 : *pb = diviiexact(b, d); return d;
2347 : }
2348 :
2349 : static int
2350 602537 : negcmpii(void *E, GEN x, GEN y) { (void)E; return -cmpii(x,y); }
2351 :
2352 : /* x square of maximal rank; does b = x[i,i] divide all entries in
2353 : * x[1..i-1, 1..i-1] ? If so, return 0; else the index of a problematic row */
2354 : static long
2355 1302063 : ZM_snf_no_divide(GEN x, long i)
2356 : {
2357 1302063 : GEN b = gcoeff(x,i,i);
2358 : long j, k;
2359 :
2360 1302063 : if (is_pm1(b)) return 0;
2361 1163802 : for (k = 1; k < i; k++)
2362 2717285 : for (j = 1; j < i; j++)
2363 2102064 : if (!dvdii(gcoeff(x,k,j),b)) return k;
2364 375969 : return 0;
2365 : }
2366 :
2367 : static void
2368 1690757 : ZM_redpart(GEN x, GEN p, long I)
2369 : {
2370 1690757 : long l = lgefint(p), i, j;
2371 7229121 : for (i = 1; i <= I; i++)
2372 32400646 : for (j = 1; j <= I; j++)
2373 : {
2374 26862282 : GEN c = gcoeff(x,i,j);
2375 26862282 : if (lgefint(c) > l) gcoeff(x,i,j) = remii(c, p);
2376 : }
2377 1690757 : }
2378 : static void
2379 1071268 : ZMrow_divexact_inplace(GEN M, long i, GEN c)
2380 : {
2381 1071268 : long j, l = lg(M);
2382 4263293 : for (j = 1; j < l; j++) gcoeff(M,i,j) = diviiexact(gcoeff(M,i,j), c);
2383 1071268 : }
2384 :
2385 : /* Return the SNF D of matrix X. If ptU/ptV non-NULL set them to U/V
2386 : * to that D = UXV */
2387 : GEN
2388 848676 : ZM_snfall_i(GEN x, GEN *ptU, GEN *ptV, long flag)
2389 : {
2390 848676 : pari_sp av0 = avma, av;
2391 848676 : const long return_vec = flag & 1;
2392 : long i, j, k, m0, m, n0, n;
2393 848676 : GEN u, v, U, V, V0, mdet, A = NULL, perm = NULL;
2394 :
2395 848676 : n0 = n = lg(x)-1;
2396 848676 : if (!n) {
2397 41594 : if (ptU) *ptU = cgetg(1,t_MAT);
2398 41594 : if (ptV) *ptV = cgetg(1,t_MAT);
2399 41594 : return cgetg(1, return_vec? t_VEC: t_MAT);
2400 : }
2401 807082 : m0 = m = nbrows(x);
2402 :
2403 807082 : U = V = V0 = NULL; /* U = TRANSPOSE of row transform matrix [act on columns]*/
2404 807082 : if (m == n && ZM_ishnf(x))
2405 : {
2406 737338 : mdet = ZM_det_triangular(x); /* != 0 */
2407 737338 : if (ptV) *ptV = matid(n);
2408 : }
2409 : else
2410 : {
2411 69744 : mdet = ZM_detmult(x);
2412 69744 : if (!signe(mdet))
2413 77 : x = ZM_hnfperm(x, ptV, ptU? &perm: NULL);
2414 : else
2415 : { /* m <= n */
2416 69667 : if (!ptV)
2417 9576 : x = ZM_hnfmod(x,mdet);
2418 60091 : else if (m == n)
2419 : {
2420 60063 : GEN H = ZM_hnfmod(x,mdet);
2421 60063 : *ptV = ZM_gauss(x,H);
2422 60063 : x = H;
2423 : }
2424 : else
2425 28 : x = ZM_hnfperm(x, ptV, ptU? &perm: NULL);
2426 69667 : mdet = ZM_det_triangular(x); /* > 0 */
2427 : }
2428 69744 : n = lg(x)-1; /* n independent columns */
2429 69744 : if (ptV)
2430 : {
2431 60112 : V = *ptV;
2432 60112 : if (n != n0)
2433 : {
2434 28 : V0 = vecslice(V, 1, n0 - n); /* kernel */
2435 28 : V = vecslice(V, n0-n+1, n0);
2436 : }
2437 : }
2438 69744 : if (!signe(mdet))
2439 : {
2440 77 : if (n)
2441 : {
2442 70 : x = ZM_snfall_i(shallowtrans(x), ptV, ptU, return_vec); /* swap V,U */
2443 70 : if (!return_vec && n != m) x = shallowtrans(x);
2444 70 : if (ptV) V = ZM_mul(V, shallowtrans(*ptV));
2445 70 : if (ptU) U = *ptU; /* TRANSPOSE */
2446 : }
2447 : else /* 0 matrix */
2448 : {
2449 7 : x = cgetg(1,t_MAT);
2450 7 : if (ptV) V = cgetg(1, t_MAT);
2451 7 : if (ptU) U = matid(m);
2452 : }
2453 77 : goto THEEND;
2454 : }
2455 : }
2456 807005 : if (ptV || ptU) U = matid(n); /* we will compute V in terms of U */
2457 807005 : if (DEBUGLEVEL>7) err_printf("starting SNF loop");
2458 :
2459 : /* square, maximal rank n */
2460 807005 : A = x; x = shallowcopy(x); av = avma;
2461 1936456 : for (i = n; i > 1; i--)
2462 : {
2463 1129451 : if (DEBUGLEVEL>7) err_printf("\ni = %ld: ",i);
2464 : for(;;)
2465 864358 : {
2466 1993809 : int c = 0;
2467 : GEN a, b;
2468 5780341 : for (j = i-1; j >= 1; j--)
2469 : {
2470 3786532 : b = gcoeff(x,i,j); if (!signe(b)) continue;
2471 332935 : a = gcoeff(x,i,i);
2472 332935 : ZC_elem(b, a, x,NULL, j,i);
2473 332935 : if (gc_needed(av,1))
2474 : {
2475 0 : if (DEBUGMEM>1) pari_warn(warnmem,"[1]: ZM_snfall i = %ld", i);
2476 0 : snf_pile1(av, &x,&U);
2477 : }
2478 : }
2479 1993809 : if (DEBUGLEVEL>7) err_printf("; ");
2480 5780341 : for (j=i-1; j>=1; j--)
2481 : {
2482 : GEN d;
2483 3786532 : b = gcoeff(x,j,i); if (!signe(b)) continue;
2484 891006 : a = gcoeff(x,i,i);
2485 891006 : d = bezout_step(&a, &b, &u, &v);
2486 4869984 : for (k = 1; k < i; k++)
2487 : {
2488 3978978 : GEN t = addii(mulii(u,gcoeff(x,i,k)),mulii(v,gcoeff(x,j,k)));
2489 3978978 : gcoeff(x,j,k) = subii(mulii(a,gcoeff(x,j,k)),
2490 3978978 : mulii(b,gcoeff(x,i,k)));
2491 3978978 : gcoeff(x,i,k) = t;
2492 : }
2493 891006 : gcoeff(x,j,i) = gen_0;
2494 891006 : gcoeff(x,i,i) = d;
2495 891006 : if (U) update(u,v,a,b,(GEN*)(U+i),(GEN*)(U+j));
2496 891006 : if (gc_needed(av,1))
2497 : {
2498 48 : if (DEBUGMEM>1) pari_warn(warnmem,"[2]: ZM_snfall, i = %ld", i);
2499 48 : snf_pile1(av, &x,&U);
2500 : }
2501 891006 : c = 1;
2502 : }
2503 1993809 : if (!c)
2504 : {
2505 1302063 : k = ZM_snf_no_divide(x, i);
2506 1302063 : if (!k) break;
2507 :
2508 : /* x[k,j] != 0 mod b */
2509 670294 : for (j = 1; j <= i; j++)
2510 497682 : gcoeff(x,i,j) = addii(gcoeff(x,i,j),gcoeff(x,k,j));
2511 172612 : if (U) gel(U,i) = gadd(gel(U,i),gel(U,k));
2512 : }
2513 864358 : ZM_redpart(x, mdet, i);
2514 864358 : if (U && (flag & 2)) ZM_redpart(U, mdet, n);
2515 864358 : if (gc_needed(av,1))
2516 : {
2517 0 : if (DEBUGMEM>1) pari_warn(warnmem,"[3]: ZM_snfall");
2518 0 : snf_pile1(av, &x,&U);
2519 : }
2520 : }
2521 : }
2522 807005 : if (DEBUGLEVEL>7) err_printf("\n");
2523 2743461 : for (k = n; k; k--)
2524 : {
2525 1936456 : GEN d = gcdii(gcoeff(x,k,k), mdet);
2526 1936456 : gcoeff(x,k,k) = d;
2527 1936456 : if (!is_pm1(d)) mdet = diviiexact(mdet,d);
2528 : }
2529 807005 : THEEND:
2530 807082 : if (U) U = shallowtrans(U);
2531 807082 : if (ptV && A)
2532 : { /* U A V = D => D^(-1) U A = V^(-1) */
2533 787363 : long l = lg(x);
2534 787363 : GEN W = ZM_mul(U, A);
2535 1858631 : for (i = 1; i < l; i++)
2536 : {
2537 1616100 : GEN c = gcoeff(x,i,i);
2538 1616100 : if (is_pm1(c)) break; /* only 1 from now on */
2539 1071268 : ZMrow_divexact_inplace(W, i, c);
2540 : }
2541 787363 : if (flag & 2)
2542 : {
2543 763346 : W = FpM_red(W, gcoeff(x,1,1));
2544 763346 : W = matinvmod(W, gcoeff(x,1,1));
2545 : }
2546 : else
2547 24017 : W = ZM_inv(W, NULL);
2548 787363 : V = V? ZM_mul(V, W): W;
2549 : }
2550 807082 : if (return_vec)
2551 : {
2552 782876 : long l = lg(x)-1;
2553 782876 : if (typ(x) == t_MAT) x = RgM_diagonal_shallow(x);
2554 782876 : if (m0 > l) x = shallowconcat(zerovec(m0-l), x);
2555 : }
2556 :
2557 807082 : if (V0)
2558 : { /* add kernel */
2559 28 : if (!return_vec) x = shallowconcat(zeromat(m,n0-n), x);
2560 28 : if (ptV) V = shallowconcat(V0, V);
2561 : }
2562 807082 : if (perm && U) U = vecpermute(U, perm_inv(perm));
2563 807082 : snf_pile(av0, &x,&U,&V);
2564 807082 : if (ptU) *ptU = U;
2565 807082 : if (ptV) *ptV = V;
2566 807082 : return x;
2567 : }
2568 : GEN
2569 64246 : ZM_snfall(GEN x, GEN *U, GEN *V) { return ZM_snfall_i(x, U, V, 0); }
2570 : GEN
2571 16184 : ZM_snf(GEN x) { return ZM_snfall_i(x, NULL,NULL, 1); }
2572 : GEN
2573 385 : smith(GEN x) { return ZM_snfall_i(x, NULL,NULL, 1); }
2574 : GEN
2575 35 : smithall(GEN x)
2576 : {
2577 35 : GEN z = cgetg(4, t_VEC);
2578 35 : gel(z,3) = ZM_snfall_i(x, (GEN*)(z+1),(GEN*)(z+2), 0);
2579 35 : return z;
2580 : }
2581 :
2582 : GEN
2583 227101 : ZV_snfclean(GEN d)
2584 : {
2585 227101 : long c, l = lg(d);
2586 239470 : for (c = 1; c < l; c++) { GEN t = gel(d,c); if (is_pm1(t)) break; }
2587 227101 : return c == l? d: vec_shorten(d, c-1);
2588 : }
2589 : void
2590 1081644 : ZM_snfclean(GEN d, GEN u, GEN v)
2591 : {
2592 1081644 : long i, c, l = lg(d);
2593 :
2594 1081644 : if (typ(d) == t_VEC)
2595 2837867 : for (c=1; c<l; c++) { GEN t = gel(d,c); if (is_pm1(t)) break; }
2596 : else
2597 : {
2598 0 : for (c=1; c<l; c++) { GEN t = gcoeff(d,c,c); if (is_pm1(t)) break; }
2599 0 : if (c < l) for (i = 1; i < c; i++) setlg(gel(d,i), c);
2600 : }
2601 1081644 : setlg(d, c);
2602 3293978 : if (u) for (i=1; i<l; i++) setlg(gel(u,i), c);
2603 1081644 : if (v) setlg(v, c);
2604 1081644 : }
2605 :
2606 : /* Assume z was computed by [g]smithall(). Remove the 1s on the diagonal */
2607 : GEN
2608 952 : smithclean(GEN z)
2609 : {
2610 : long i, j, h, l, c, d;
2611 : GEN U, V, y, D, t;
2612 :
2613 952 : if (typ(z) != t_VEC) pari_err_TYPE("smithclean",z);
2614 952 : l = lg(z); if (l == 1) return cgetg(1,t_VEC);
2615 938 : U = gel(z,1);
2616 938 : if (l != 4 || typ(U) != t_MAT)
2617 : { /* assume z = vector of elementary divisors */
2618 2366 : for (c=1; c<l; c++)
2619 1806 : if (gequal1(gel(z,c))) break;
2620 917 : return gcopy_lg(z, c);
2621 : }
2622 21 : V = gel(z,2);
2623 21 : D = gel(z,3);
2624 21 : l = lg(D);
2625 21 : if (l == 1) return gcopy(z);
2626 21 : h = lgcols(D);
2627 21 : if (h > l)
2628 : { /* D = vconcat(zero matrix, diagonal matrix) */
2629 21 : for (c=1+h-l, d=1; c<h; c++,d++)
2630 21 : if (gequal1(gcoeff(D,c,d))) break;
2631 : }
2632 7 : else if (h < l)
2633 : { /* D = concat(zero matrix, diagonal matrix) */
2634 7 : for (c=1, d=1+l-h; d<l; c++,d++)
2635 7 : if (gequal1(gcoeff(D,c,d))) break;
2636 : }
2637 : else
2638 : { /* D diagonal */
2639 0 : for (c=1; c<l; c++)
2640 0 : if (gequal1(gcoeff(D,c,c))) break;
2641 0 : d = c;
2642 : }
2643 : /* U was (h-1)x(h-1), V was (l-1)x(l-1), D was (h-1)x(l-1) */
2644 21 : y = cgetg(4,t_VEC);
2645 : /* truncate U to (c-1) x (h-1) */
2646 21 : gel(y,1) = t = cgetg(h,t_MAT);
2647 77 : for (j=1; j<h; j++) gel(t,j) = gcopy_lg(gel(U,j), c);
2648 : /* truncate V to (l-1) x (d-1) */
2649 21 : gel(y,2) = gcopy_lg(V, d);
2650 21 : gel(y,3) = t = zeromatcopy(c-1, d-1);
2651 : /* truncate D to a (c-1) x (d-1) matrix */
2652 21 : if (d > 1)
2653 : {
2654 14 : if (h > l)
2655 : {
2656 14 : for (i=1+h-l, j=1; i<c; i++,j++)
2657 7 : gcoeff(t,i,j) = gcopy(gcoeff(D,i,j));
2658 : }
2659 7 : else if (h < l)
2660 : {
2661 7 : for (i=1, j=1+l-h; j<d; i++,j++)
2662 0 : gcoeff(t,i,j) = gcopy(gcoeff(D,i,j));
2663 : }
2664 : else
2665 : {
2666 0 : for (j=1; j<d; j++)
2667 0 : gcoeff(t,j,j) = gcopy(gcoeff(D,j,j));
2668 : }
2669 : }
2670 21 : return y;
2671 : }
2672 :
2673 : /* does b = x[i,i] divide all entries in x[1..i-1,1..i-1] ? If so, return 0;
2674 : * else return the index of a problematic row */
2675 : static long
2676 196 : gsnf_no_divide(GEN x, long i, long vx)
2677 : {
2678 196 : GEN b = gcoeff(x,i,i);
2679 : long j, k;
2680 :
2681 196 : if (gequal0(b))
2682 : {
2683 14 : for (k = 1; k < i; k++)
2684 14 : for (j = 1; j < i; j++)
2685 14 : if (!gequal0(gcoeff(x,k,j))) return k;
2686 0 : return 0;
2687 : }
2688 :
2689 182 : if (!is_RgX(b,vx) || degpol(b)<=0) return 0;
2690 217 : for (k = 1; k < i; k++)
2691 378 : for (j = 1; j < i; j++)
2692 : {
2693 266 : GEN z = gcoeff(x,k,j), r;
2694 266 : if (!is_RgX(z,vx)) z = scalarpol(z, vx);
2695 266 : r = RgX_rem(z, b);
2696 266 : if (signe(r) && (! isinexactreal(r) ||
2697 0 : gexpo(r) > 16 + gexpo(b) - gprecision(r))
2698 35 : ) return k;
2699 : }
2700 70 : return 0;
2701 : }
2702 :
2703 : /* Hermite Normal Form, with base change matrix if ptB != NULL.
2704 : * If 'remove' = 1, remove 0 columns (do NOT update *ptB accordingly)
2705 : * If 'remove' = 2, remove 0 columns and update *ptB accordingly */
2706 : GEN
2707 42 : RgM_hnfall(GEN A, GEN *pB, long remove)
2708 : {
2709 : pari_sp av;
2710 : long li, j, k, m, n, def, ldef;
2711 : GEN B;
2712 42 : long vx = gvar(A);
2713 :
2714 42 : n = lg(A)-1;
2715 42 : if (vx==NO_VARIABLE || !n)
2716 : {
2717 0 : RgM_check_ZM(A, "mathnf0");
2718 0 : return ZM_hnfall(A, pB, remove);
2719 : }
2720 42 : m = nbrows(A);
2721 42 : av = avma;
2722 42 : A = RgM_shallowcopy(A);
2723 42 : B = pB? matid(n): NULL;
2724 42 : def = n; ldef = (m>n)? m-n: 0;
2725 126 : for (li=m; li>ldef; li--)
2726 : {
2727 : GEN d, T;
2728 714 : for (j=def-1; j; j--)
2729 : {
2730 630 : GEN a = gcoeff(A,li,j);
2731 630 : if (gequal0(a)) continue;
2732 :
2733 511 : k = (j==1)? def: j-1;
2734 511 : RgC_elem(a,gcoeff(A,li,k), A,B, j,k, li, vx);
2735 : }
2736 84 : T = normalize_as_RgX(gcoeff(A,li,def), vx, &d);
2737 84 : if (gequal0(T))
2738 0 : { if (ldef) ldef--; }
2739 : else
2740 : {
2741 84 : if (!gequal1(d))
2742 : {
2743 7 : gel(A,def) = RgC_Rg_div(gel(A,def), d);
2744 7 : if (B) gel(B, def) = RgC_Rg_div(gel(B, def), d);
2745 : }
2746 84 : RgM_reduce(A, B, li, def, vx);
2747 84 : def--;
2748 : }
2749 84 : if (gc_needed(av,1))
2750 : {
2751 0 : if (DEBUGMEM>1) pari_warn(warnmem,"ghnfall");
2752 0 : (void)gc_all(av, B? 2: 1, &A, &B);
2753 : }
2754 : }
2755 : /* rank A = n - def */
2756 42 : if (remove) remove_0cols(def, &A, &B, remove);
2757 42 : (void)gc_all(av, B? 2: 1, &A, &B);
2758 42 : if (B) *pB = B;
2759 42 : return A;
2760 : }
2761 :
2762 : static GEN
2763 49 : RgXM_snf(GEN x,long all)
2764 : {
2765 : pari_sp av;
2766 : long i, j, k, n;
2767 : GEN z, u, v, U, V;
2768 49 : long vx = gvar(x);
2769 49 : n = lg(x)-1; if (!n) return trivsmith(all);
2770 49 : if (vx==NO_VARIABLE) pari_err_TYPE("RgXM_snf",x);
2771 49 : if (lgcols(x) != n+1) pari_err_DIM("gsmithall");
2772 49 : av = avma;
2773 49 : x = RgM_shallowcopy(x);
2774 49 : if (all) { U = matid(n); V = matid(n); }
2775 196 : for (i=n; i>=2; i--)
2776 : {
2777 : for(;;)
2778 168 : {
2779 : GEN a, b, d;
2780 315 : int c = 0;
2781 1092 : for (j=i-1; j>=1; j--)
2782 : {
2783 777 : b = gcoeff(x,i,j); if (gequal0(b)) continue;
2784 196 : a = gcoeff(x,i,i);
2785 196 : d = gbezout_step(&b, &a, &v, &u, vx);
2786 700 : for (k = 1; k < i; k++)
2787 : {
2788 504 : GEN t = gadd(gmul(u,gcoeff(x,k,i)),gmul(v,gcoeff(x,k,j)));
2789 504 : gcoeff(x,k,j) = gsub(gmul(a,gcoeff(x,k,j)),gmul(b,gcoeff(x,k,i)));
2790 504 : gcoeff(x,k,i) = t;
2791 : }
2792 196 : gcoeff(x,i,j) = gen_0;
2793 196 : gcoeff(x,i,i) = d;
2794 196 : if (all) update(u,v,a,b,(GEN*)(V+i),(GEN*)(V+j));
2795 : }
2796 1092 : for (j=i-1; j>=1; j--)
2797 : {
2798 777 : b = gcoeff(x,j,i); if (gequal0(b)) continue;
2799 175 : a = gcoeff(x,i,i);
2800 175 : d = gbezout_step(&b, &a, &v, &u, vx);
2801 651 : for (k = 1; k < i; k++)
2802 : {
2803 476 : GEN t = gadd(gmul(u,gcoeff(x,i,k)),gmul(v,gcoeff(x,j,k)));
2804 476 : gcoeff(x,j,k) = gsub(gmul(a,gcoeff(x,j,k)),gmul(b,gcoeff(x,i,k)));
2805 476 : gcoeff(x,i,k) = t;
2806 : }
2807 175 : gcoeff(x,j,i) = gen_0;
2808 175 : gcoeff(x,i,i) = d;
2809 175 : if (all) update(u,v,a,b,(GEN*)(U+i),(GEN*)(U+j));
2810 175 : c = 1;
2811 : }
2812 315 : if (!c)
2813 : {
2814 196 : k = gsnf_no_divide(x, i, vx);
2815 196 : if (!k) break;
2816 :
2817 245 : for (j=1; j<=i; j++)
2818 196 : gcoeff(x,i,j) = gadd(gcoeff(x,i,j),gcoeff(x,k,j));
2819 49 : if (all) gel(U,i) = gadd(gel(U,i),gel(U,k));
2820 : }
2821 168 : if (gc_needed(av,1))
2822 : {
2823 0 : if (DEBUGMEM>1) pari_warn(warnmem,"gsmithall");
2824 0 : (void)gc_all(av, all? 3: 1, &x, &U, &V);
2825 : }
2826 : }
2827 : }
2828 245 : for (k=1; k<=n; k++)
2829 : {
2830 196 : GEN d, T = normalize_as_RgX(gcoeff(x,k,k), vx, &d);
2831 196 : if (gequal0(T)) continue;
2832 182 : if (all && !gequal1(d)) gel(V,k) = RgC_Rg_div(gel(V,k), d);
2833 182 : gcoeff(x,k,k) = T;
2834 : }
2835 49 : z = all? mkvec3(shallowtrans(U), V, x): RgM_diagonal_shallow(x);
2836 49 : return gc_GEN(av, z);
2837 : }
2838 :
2839 : GEN
2840 469 : matsnf0(GEN x,long flag)
2841 : {
2842 469 : pari_sp av = avma;
2843 469 : if (flag > 7) pari_err_FLAG("matsnf");
2844 469 : if (typ(x) == t_VEC && flag & 4) return smithclean(x);
2845 469 : if (typ(x)!=t_MAT) pari_err_TYPE("matsnf",x);
2846 469 : if (RgM_is_ZM(x)) x = flag&1 ? smithall(x): smith(x);
2847 49 : else x = RgXM_snf(x, flag&1);
2848 469 : if (flag & 4) x = gc_upto(av, smithclean(x));
2849 469 : return x;
2850 : }
2851 : GEN
2852 0 : gsmith(GEN x) { return RgXM_snf(x,0); }
2853 : GEN
2854 0 : gsmithall(GEN x) { return RgXM_snf(x,1); }
2855 :
2856 : /* H is a relation matrix, either in HNF or a t_VEC (diagonal HNF) */
2857 : static GEN
2858 1081644 : snf_group(GEN H, GEN D, GEN *newU, GEN *newUi)
2859 : {
2860 : long i, j, l;
2861 :
2862 1081644 : ZM_snfclean(D, newU? *newU: NULL, newUi? *newUi: NULL);
2863 1081644 : l = lg(D);
2864 1081644 : if (newU) {
2865 902966 : GEN U = *newU;
2866 2375057 : for (i = 1; i < l; i++)
2867 : {
2868 1472091 : GEN d = gel(D,i), d2 = shifti(d, 1);
2869 5975630 : for (j = 1; j < lg(U); j++)
2870 4503539 : gcoeff(U,i,j) = centermodii(gcoeff(U,i,j), d, d2);
2871 : }
2872 902966 : *newU = U;
2873 : }
2874 1081644 : if (newUi && l > 1)
2875 : { /* UHV=D -> U^-1 = (HV)D^-1 -> U^-1 = H(VD^-1 mod 1) mod H */
2876 : /* Ui = ZM_inv(U, NULL); setlg(Ui, l); */
2877 994116 : GEN V = *newUi, Ui;
2878 994116 : int Hvec = (typ(H) == t_VEC);
2879 2749779 : for (i = 1; i < l; i++) gel(V,i) = FpC_red(gel(V,i), gel(D,i));
2880 994116 : if (!Hvec)
2881 : {
2882 677442 : if (ZM_isdiagonal(H)) { H = RgM_diagonal_shallow(H); Hvec = 1; }
2883 : }
2884 994116 : Ui = Hvec? ZM_diag_mul(H, V): ZM_mul(H, V);
2885 2749779 : for (i = 1; i < l; i++) gel(Ui,i) = ZC_Z_divexact(gel(Ui,i), gel(D,i));
2886 994116 : *newUi = Hvec? ZM_ZV_mod(Ui, H): ZM_hnfrem(Ui, H);
2887 : }
2888 1081644 : return D;
2889 : }
2890 : /* H relation matrix among row of generators g in HNF. Let URV = D its SNF,
2891 : * newU R newV = newD its clean SNF (no 1 in Dnew). Return the diagonal of
2892 : * newD, newU and newUi such that 1/U = (newUi, ?).
2893 : * Rationale: let (G,0) = g Ui be the new generators then
2894 : * 0 = G U R --> G D = 0, g = G newU and G = g newUi */
2895 : GEN
2896 764697 : ZM_snf_group(GEN H, GEN *newU, GEN *newUi)
2897 : {
2898 764697 : GEN D = ZM_snfall_i(H, newU, newUi, 1 + 2);
2899 764697 : return snf_group(H, D, newU, newUi);
2900 : }
2901 :
2902 : /* D a ZV: SNF for matdiagonal(D). Faster because we only ensure elementary
2903 : * divisors condition: d[n] | ... | d[1] and need not convert D to matrix form*/
2904 : GEN
2905 316947 : ZV_snfall(GEN D, GEN *pU, GEN *pV)
2906 : {
2907 316947 : pari_sp av = avma;
2908 316947 : long j, n = lg(D)-1;
2909 316947 : GEN U = pU? matid(n): NULL;
2910 316947 : GEN V = pV? matid(n): NULL;
2911 : GEN p;
2912 :
2913 316947 : D = leafcopy(D);
2914 1043167 : for (j = n; j > 0; j--)
2915 : {
2916 726220 : GEN b = gel(D,j);
2917 726220 : if (signe(b) < 0)
2918 : {
2919 0 : gel(D,j) = negi(b);
2920 0 : if (V) ZV_togglesign(gel(V,j));
2921 : }
2922 : }
2923 : /* entries are nonnegative integers */
2924 316947 : p = gen_indexsort(D, NULL, &negcmpii);
2925 316947 : D = vecpermute(D, p);
2926 316947 : if (U) U = vecpermute(U, p);
2927 316947 : if (V) V = vecpermute(V, p);
2928 : /* entries are sorted by decreasing value */
2929 1043167 : for (j = n; j > 0; j--)
2930 : {
2931 726220 : GEN b = gel(D,j);
2932 : long i;
2933 1390800 : for (i = j-1; i > 0; i--)
2934 : { /* invariant: a >= b. If au+bv = d is a Bezout relation, A=a/d and B=b/d
2935 : * we have [B,-A;u,v]*diag(a,b)*[1-u*A,1; -u*A,1]] = diag(Ab, d) */
2936 676774 : GEN a = gel(D,i), u,v, d = bezout(a,b, &u,&v), A, Wi, Wj;
2937 676774 : if (equalii(d,b)) continue;
2938 70522 : A = diviiexact(a,d);
2939 70522 : if (V)
2940 : {
2941 70466 : GEN t = mulii(u,A);
2942 70466 : Wi = ZC_lincomb(subui(1,t), negi(t), gel(V,i), gel(V,j));
2943 70466 : Wj = ZC_add(gel(V,i), gel(V,j));
2944 70466 : gel(V,i) = Wi;
2945 70466 : gel(V,j) = Wj;
2946 : }
2947 70522 : if (U)
2948 : {
2949 70522 : GEN B = diviiexact(b,d);
2950 70522 : Wi = ZC_lincomb(B, negi(A), gel(U,i), gel(U,j));
2951 70522 : Wj = ZC_lincomb(u, v, gel(U,i), gel(U,j));
2952 70522 : gel(U,i) = Wi;
2953 70522 : gel(U,j) = Wj;
2954 : }
2955 70522 : gel(D,i) = mulii(A,b); /* lcm(a,b) */
2956 70522 : gel(D,j) = d; /* gcd(a,b) */
2957 70522 : b = gel(D,j); if (equali1(b)) break;
2958 : }
2959 : }
2960 316947 : snf_pile(av, &D,&U,&V);
2961 316947 : if (U) *pU = shallowtrans(U);
2962 316947 : if (V) *pV = V;
2963 316947 : return D;
2964 : }
2965 : GEN
2966 316947 : ZV_snf_group(GEN d, GEN *newU, GEN *newUi)
2967 : {
2968 316947 : GEN D = ZV_snfall(d, newU, newUi);
2969 316947 : return snf_group(d, D, newU, newUi);
2970 : }
2971 :
2972 : /* D a vector of elementary divisors. Truncate (setlg) to leave out trivial
2973 : * entries (= 1) */
2974 : void
2975 0 : ZV_snf_trunc(GEN D)
2976 : {
2977 0 : long i, l = lg(D);
2978 0 : for (i = 1; i < l; i++)
2979 0 : if (is_pm1(gel(D,i))) { setlg(D,i); break; }
2980 0 : }
2981 :
2982 : long
2983 0 : zv_snf_rank(GEN D, ulong p)
2984 : {
2985 0 : long i, l = lg(D);
2986 0 : if (!p) return l - 1;
2987 0 : for (i = 1; i < l; i++)
2988 0 : if (D[i] % p) break;
2989 0 : return i - 1;
2990 : }
2991 : long
2992 49 : ZV_snf_rank_u(GEN D, ulong p)
2993 : {
2994 49 : long i, l = lg(D);
2995 49 : while (l > 1 && D[l-1] == 1) l--;
2996 49 : if (!p) return l - 1;
2997 49 : if (p == 2)
2998 : {
2999 49 : for (i = 1; i < l; i++)
3000 42 : if (mpodd(gel(D,i))) break;
3001 : }
3002 35 : else if (!(p & (p-1)))
3003 : { /* power of 2 */
3004 14 : long n = vals(p);
3005 28 : for (i = 1; i < l; i++)
3006 28 : if (umodi2n(gel(D,i), n)) break;
3007 : }
3008 : else
3009 : {
3010 49 : for (i = 1; i < l; i++)
3011 42 : if (umodiu(gel(D,i), p)) break;
3012 : }
3013 49 : return i - 1;
3014 : }
3015 : long
3016 91 : ZV_snf_rank(GEN D, GEN p)
3017 : {
3018 91 : long i, l = lg(D);
3019 91 : if (lgefint(p) == 3) return ZV_snf_rank_u(D, p[2]);
3020 77 : while (l > 1 && equali1(gel(D, l-1))) l--;
3021 42 : if (!signe(p)) return l - 1;
3022 77 : for (i = 1; i < l; i++)
3023 70 : if (!dvdii(gel(D,i), p)) break;
3024 21 : return i - 1;
3025 : }
3026 : long
3027 154 : snfrank(GEN D, GEN p)
3028 : {
3029 : long i, l;
3030 154 : if (typ(D) != t_VEC) pari_err_TYPE("snfrank", D);
3031 154 : if (!p) p = gen_0;
3032 154 : l = lg(D);
3033 154 : if (l == 4 && typ(gel(D,3)) == t_MAT)
3034 : { /* from matsnf(,1) */
3035 14 : pari_sp av = avma;
3036 : long z;
3037 : GEN v;
3038 14 : D = gel(D,3); l = lg(D);
3039 14 : if (l == 1) return 0;
3040 14 : z = lgcols(D) - l; /* missing columns of 0s */
3041 14 : if (z < 0) pari_err_TYPE("snfrank", D);
3042 14 : v = cgetg(l, t_VEC);
3043 35 : for (i = 1; i < l; i++) gel(v, i) = gcoeff(D, i + z, i);
3044 14 : return gc_long(av, z + snfrank(v, p));
3045 : }
3046 140 : switch(typ(p))
3047 : {
3048 98 : case t_INT:
3049 98 : if (RgV_is_ZV(D)) return ZV_snf_rank(D, p);
3050 7 : if (!signe(p)) break; /* allow p = 0 */
3051 0 : pari_err_TYPE("snfrank", D);
3052 42 : case t_POL: break;
3053 0 : default: pari_err_TYPE("snfrank", p);
3054 : }
3055 175 : while (l > 1 && gequal1(gel(D, l-1))) l--;
3056 49 : if (gequal0(p)) return l - 1;
3057 112 : for (i = 1; i < l; i++)
3058 91 : if (!gdvd(gel(D,i), p)) break;
3059 42 : return i - 1;
3060 : }
3061 :
3062 : /* Assume L1/L2 are of maximal rank */
3063 : GEN
3064 28 : ZM_hnfintersectmod(GEN L1, GEN L2, GEN N)
3065 : {
3066 28 : pari_sp av = avma;
3067 28 : long n = lg(L1)-1;
3068 28 : GEN L1N = FpM_red(L1, N);
3069 28 : GEN K = matkermod(shallowconcat(L1N, L2), N, NULL);
3070 28 : GEN H = ZM_hnfmodid(FpM_mul(L1N, rowslice(K, 1, n), N), N);
3071 28 : return gc_upto(av, H);
3072 : }
|
