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 259 : mathnf0(GEN x, long flag)
35 : {
36 259 : switch(typ(x))
37 : {
38 70 : case t_VEC:
39 70 : if (RgV_is_ZV(x))
40 : 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 21 : case 1:
46 : case 4:
47 21 : return ZV_hnfgcdext(x);
48 : }
49 35 : x = gtomat(x); break;
50 189 : case t_MAT: break;
51 0 : default: pari_err_TYPE("mathnf0",x);
52 : }
53 :
54 224 : switch(flag)
55 : {
56 168 : 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 8035855 : count(GEN mat, long row, long len, long *firstnonzero)
72 : {
73 8035855 : long j, n = 0;
74 :
75 722586859 : for (j=1; j<=len; j++)
76 : {
77 716259827 : long p = mael(mat,j,row);
78 716259827 : if (p)
79 : {
80 22822856 : if (labs(p)!=1) return -1;
81 21114033 : n++; *firstnonzero=j;
82 : }
83 : }
84 6327032 : return n;
85 : }
86 :
87 : static int
88 404551 : count2(GEN mat, long row, long len)
89 : {
90 : long j;
91 3961782 : for (j=len; j; j--)
92 3841845 : if (labs(mael(mat,j,row)) == 1) return j;
93 119937 : return 0;
94 : }
95 :
96 : static GEN
97 207114 : hnffinal(GEN matgen,GEN perm,GEN* ptdep,GEN* ptB,GEN* ptC)
98 : {
99 : GEN p1,p2,U,H,Hnew,Bnew,Cnew,diagH1;
100 207114 : GEN B = *ptB, C = *ptC, dep = *ptdep, depnew;
101 : pari_sp av;
102 : long i,j,k,s,i1,j1,zc;
103 207114 : long co = lg(C);
104 207114 : long col = lg(matgen)-1;
105 : long lnz, nlze, lig;
106 :
107 207114 : if (col == 0) return matgen;
108 207114 : lnz = nbrows(matgen); /* was called lnz-1 - nr in hnfspec */
109 207114 : nlze = nbrows(dep);
110 207114 : lig = nlze + lnz;
111 : /* H: lnz x lnz [disregarding initial 0 cols], U: col x col */
112 207114 : H = ZM_hnflll(matgen, &U, 0);
113 207113 : 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 207113 : zc = col - lnz; /* # of 0 columns, correspond to units */
117 207113 : if (nlze) { dep = ZM_mul(dep,U); dep += zc; }
118 :
119 207113 : diagH1 = new_chunk(lnz+1); /* diagH1[i] = 0 iff H[i,i] != 1 (set later) */
120 :
121 207114 : av = avma;
122 207114 : Cnew = cgetg(co, typ(C));
123 207113 : setlg(C, col+1); p1 = gmul(C,U);
124 1848032 : for (j=1; j<=col; j++) gel(Cnew,j) = gel(p1,j);
125 3058899 : for ( ; j<co ; j++) gel(Cnew,j) = gel(C,j);
126 :
127 : /* Clean up B using new H */
128 842745 : for (s=0,i=lnz; i; i--)
129 : {
130 635607 : GEN Di = gel(dep,i), Hi = gel(H,i);
131 635607 : GEN h = gel(Hi,i); /* H[i,i] */
132 635607 : if ( (diagH1[i] = is_pm1(h)) ) { h = NULL; s++; }
133 15481699 : for (j=col+1; j<co; j++)
134 : {
135 14846374 : GEN z = gel(B,j-col);
136 14846374 : p1 = gel(z,i+nlze);
137 14846374 : if (h) p1 = truedivii(p1,h);
138 14844462 : if (!signe(p1)) continue;
139 23825745 : for (k=1; k<=nlze; k++) gel(z,k) = subii(gel(z,k), mulii(p1, gel(Di,k)));
140 122408068 : for ( ; k<=lig; k++) gel(z,k) = subii(gel(z,k), mulii(p1, gel(Hi,k-nlze)));
141 8848103 : gel(Cnew,j) = gsub(gel(Cnew,j), gmul(p1, gel(Cnew,i+zc)));
142 : }
143 635325 : if (gc_needed(av,2))
144 : {
145 714 : if(DEBUGMEM>1) pari_warn(warnmem,"hnffinal, i = %ld",i);
146 714 : gerepileall(av, 2, &Cnew, &B);
147 : }
148 : }
149 207138 : p1 = cgetg(lnz+1,t_VEC); p2 = perm + nlze;
150 842785 : for (i1=0, j1=lnz-s, i=1; i<=lnz; i++) /* push the 1 rows down */
151 635674 : if (diagH1[i])
152 303262 : gel(p1,++j1) = gel(p2,i);
153 : else
154 332412 : gel(p2,++i1) = gel(p2,i);
155 510375 : 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 207111 : lig -= s; col -= s; lnz -= s;
166 207111 : Hnew = cgetg(lnz+1,t_MAT);
167 207113 : depnew = cgetg(lnz+1,t_MAT); /* only used if nlze > 0 */
168 207113 : Bnew = cgetg(co-col,t_MAT);
169 207113 : C = shallowcopy(Cnew);
170 842794 : for (j=1,i1=j1=0; j<=lnz+s; j++)
171 : {
172 635682 : GEN z = gel(H,j);
173 635682 : if (diagH1[j])
174 : { /* hit exactly s times */
175 303266 : i1++; C[i1+col] = Cnew[j+zc];
176 303266 : p1 = cgetg(lig+1,t_COL); gel(Bnew,i1) = p1;
177 525814 : for (i=1; i<=nlze; i++) gel(p1,i) = gcoeff(dep,i,j);
178 303266 : p1 += nlze;
179 : }
180 : else
181 : {
182 332416 : j1++; C[j1+zc] = Cnew[j+zc];
183 332416 : p1 = cgetg(lnz+1,t_COL); gel(Hnew,j1) = p1;
184 332411 : depnew[j1] = dep[j];
185 : }
186 2414614 : for (i=k=1; k<=lnz; i++)
187 1778937 : if (!diagH1[i]) p1[k++] = z[i];
188 : }
189 3058823 : for (j=s+1; j<co-col; j++)
190 : {
191 2851706 : GEN z = gel(B,j-s);
192 2851706 : p1 = cgetg(lig+1,t_COL); gel(Bnew,j) = p1;
193 5369938 : for (i=1; i<=nlze; i++) gel(p1,i) = gel(z,i);
194 2851711 : z += nlze; p1 += nlze;
195 11960621 : for (i=k=1; k<=lnz; i++)
196 9108910 : if (!diagH1[i]) gel(p1,k++) = gel(z,i);
197 : }
198 207117 : *ptdep = depnew;
199 207117 : *ptC = C;
200 207117 : *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 2652800 : col_dup(long l, GEN col)
217 : {
218 2652800 : GEN c = new_chunk(l);
219 2652792 : 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 207113 : ZM_rowrankprofile(GEN x, long *nlze)
225 : {
226 207113 : pari_sp av = avma;
227 : GEN d, y;
228 : long i, j, k, l, r;
229 :
230 207113 : x = shallowtrans(x); l = lg(x);
231 207114 : (void)new_chunk(l); /* HACK */
232 207114 : d = ZM_pivots(x,&r); set_avma(av);
233 207111 : *nlze = r;
234 207111 : if (!d) return identity_perm(l-1);
235 192709 : y = cgetg(l,t_VECSMALL);
236 837900 : for (i = j = 1, k = r+1; i<l; i++)
237 645191 : if (d[i]) y[k++] = i; else y[j++] = i;
238 192709 : 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 133924 : 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 133924 : const long li = lg(perm); /* = lgcols(mat0) */
262 133924 : const long CO = lg(mat0);
263 :
264 133924 : n = 0; /* -Wall */
265 :
266 133924 : C = *ptC; co = CO;
267 133924 : 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 133924 : if (DEBUGLEVEL>5)
274 : {
275 0 : err_printf("Entering hnfspec\n");
276 0 : p_mat(mat0,perm,0);
277 : }
278 133924 : matt = cgetg(co, t_MAT); /* dense part of mat (top) */
279 133923 : mat = cgetg(co, t_MAT);
280 2786729 : for (j = 1; j < co; j++)
281 : {
282 2652804 : GEN matj = col_dup(li, gel(mat0,j));
283 2652800 : p1 = cgetg(k0+1,t_COL); gel(matt,j) = p1; gel(mat,j) = matj;
284 12334509 : for (i=1; i<=k0; i++) gel(p1,i) = stoi(matj[perm[i]]);
285 : }
286 133925 : av = avma;
287 :
288 133925 : i = lig = li-1; col = co-1; lk0 = k0;
289 133925 : 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 6380455 : while (i > lk0 && col)
292 6246471 : switch( count(mat,perm[i],col,&n) )
293 : {
294 27406 : case 0: /* move zero lines between k0+1 and lk0 */
295 27406 : lk0++; lswap(perm[i], perm[lk0]);
296 27406 : i = lig; continue;
297 :
298 482358 : case 1: /* move trivial generator between lig+1 and li */
299 482358 : lswap(perm[i], perm[lig]);
300 482358 : if (T) swap(gel(T,n), gel(T,col));
301 482358 : swap(gel(mat,n), gel(mat,col)); p = gel(mat,col);
302 482358 : if (p[perm[lig]] < 0) /* = -1 */
303 : { /* convert relation -g = 0 to g = 0 */
304 10379841 : for (i=lk0+1; i<lig; i++) p[perm[i]] = -p[perm[i]];
305 435187 : if (T)
306 : {
307 435190 : p1 = gel(T,col);
308 11107447 : for (i=1; ; i++) /* T = permuted identity: single nonzero entry */
309 11107447 : if (signe(gel(p1,i))) { togglesign_safe(&gel(p1,i)); break; }
310 : }
311 : }
312 482382 : lig--; col--; i = lig; continue;
313 :
314 5736742 : default: i--;
315 : }
316 133984 : 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 133925 : s = 0;
323 776526 : while (lig > lk0 && col && s < (long)(HIGHBIT>>1))
324 : {
325 1913995 : for (i=lig; i>lk0; i--)
326 1789479 : if (count(mat,perm[i],col,&n) > 0) break;
327 767206 : if (i == lk0) break;
328 :
329 : /* only 0, +/- 1 entries, at least 2 of them nonzero */
330 642690 : lswap(perm[i], perm[lig]);
331 642690 : swap(gel(mat,n), gel(mat,col)); p = gel(mat,col);
332 642690 : if (T) swap(gel(T,n), gel(T,col));
333 642690 : if (p[perm[lig]] < 0)
334 : {
335 7520143 : for (i=lk0+1; i<=lig; i++) p[perm[i]] = -p[perm[i]];
336 414661 : if (T) ZV_togglesign(gel(T,col));
337 : }
338 20604408 : for (j=1; j<col; j++)
339 : {
340 19961807 : GEN matj = gel(mat,j);
341 : long t;
342 19961807 : if (! (t = matj[perm[lig]]) ) continue;
343 1510489 : if (t == 1) {
344 24784151 : for (i=lk0+1; i<=lig; i++) absmax(s, matj[perm[i]] -= p[perm[i]]);
345 : }
346 : else { /* t = -1 */
347 18106599 : for (i=lk0+1; i<=lig; i++) absmax(s, matj[perm[i]] += p[perm[i]]);
348 : }
349 1510489 : if (T) ZC_lincomb1_inplace(gel(T,j), gel(T,col), stoi(-t));
350 : }
351 642601 : lig--; col--;
352 642601 : if (gc_needed(av,3))
353 : {
354 0 : if(DEBUGMEM>1) pari_warn(warnmem,"hnfspec[1]");
355 0 : if (T) T = gerepilecopy(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 133841 : vmax = cgetg(co,t_VECSMALL);
361 1661727 : for (j=1; j<=col; j++)
362 : {
363 1527803 : GEN matj = gel(mat,j);
364 6952132 : for (s=0, i=lk0+1; i<=lig; i++) absmax(s, matj[i]);
365 1527803 : vmax[j] = s;
366 : }
367 418539 : while (lig > lk0 && col)
368 : {
369 430759 : for (i=lig; i>lk0; i--)
370 404551 : if ( (n = count2(mat,perm[i],col)) ) break;
371 310822 : if (i == lk0) break;
372 :
373 284616 : lswap(vmax[n], vmax[col]);
374 284616 : lswap(perm[i], perm[lig]);
375 284616 : swap(gel(mat,n), gel(mat,col)); p = gel(mat,col);
376 284616 : if (T) swap(gel(T,n), gel(T,col));
377 284616 : if (p[perm[lig]] < 0)
378 : {
379 577803 : for (i=lk0+1; i<=lig; i++) p[perm[i]] = -p[perm[i]];
380 100873 : if (T) ZV_togglesign(gel(T,col));
381 : }
382 3861266 : for (j=1; j<col; j++)
383 : {
384 3576710 : GEN matj = gel(mat,j);
385 : long t;
386 3576710 : if (! (t = matj[perm[lig]]) ) continue;
387 1793841 : if (vmax[col] && (ulong)labs(t) >= (HIGHBIT-vmax[j]) / vmax[col])
388 0 : goto END2;
389 :
390 18331421 : for (s=0, i=lk0+1; i<=lig; i++) absmax(s, matj[perm[i]] -= t*p[perm[i]]);
391 1793841 : vmax[j] = s;
392 1793841 : if (T) ZC_lincomb1_inplace(gel(T,j), gel(T,col), stoi(-t));
393 : }
394 284556 : lig--; col--;
395 284556 : if (gc_needed(av,3))
396 : {
397 0 : if(DEBUGMEM>1) pari_warn(warnmem,"hnfspec[2]");
398 0 : gerepileall(av, T? 2: 1, &vmax, &T);
399 : }
400 : }
401 :
402 107719 : END2: /* clean up mat: remove everything to the right of the 1s on diagonal */
403 : /* go multiprecision first */
404 133925 : matb = cgetg(co,t_MAT); /* bottom part (complement of matt) */
405 2786684 : for (j=1; j<co; j++)
406 : {
407 2652837 : GEN matj = gel(mat,j);
408 2652837 : p1 = cgetg(li-k0,t_COL); gel(matb,j) = p1;
409 2652853 : p1 -= k0;
410 80268976 : for (i=k0+1; i<li; i++) gel(p1,i) = stoi(matj[perm[i]]);
411 : }
412 133847 : if (DEBUGLEVEL>5)
413 : {
414 0 : err_printf(" after phase2:\n");
415 0 : p_mat(mat,perm,lk0);
416 : }
417 1410056 : for (i=li-2; i>lig; i--)
418 : {
419 1276138 : long h, i0 = i - k0, k = i + co-li;
420 1276138 : GEN Bk = gel(matb,k);
421 28181914 : for (j=k+1; j<co; j++)
422 : {
423 26906084 : GEN Bj = gel(matb,j), v = gel(Bj,i0);
424 26906084 : s = signe(v); if (!s) continue;
425 :
426 5192606 : gel(Bj,i0) = gen_0;
427 5192606 : if (is_pm1(v))
428 : {
429 3008569 : if (s > 0) /* v = 1 */
430 44145729 : { for (h=1; h<i0; h++) gel(Bj,h) = subii(gel(Bj,h), gel(Bk,h)); }
431 : else /* v = -1 */
432 37829696 : { for (h=1; h<i0; h++) gel(Bj,h) = addii(gel(Bj,h), gel(Bk,h)); }
433 : }
434 : else {
435 50899325 : for (h=1; h<i0; h++) gel(Bj,h) = subii(gel(Bj,h), mulii(v,gel(Bk,h)));
436 : }
437 5192092 : if (T) ZC_lincomb1_inplace(gel(T,j), gel(T,k), negi(v));
438 : }
439 1275830 : if (gc_needed(av,2))
440 : {
441 6 : if(DEBUGMEM>1) pari_warn(warnmem,"hnfspec[3], i = %ld", i);
442 1578 : for (h=1; h<co; h++) setlg(matb[h], i0+1); /* bottom can be forgotten */
443 6 : gerepileall(av, T? 2: 1, &matb, &T);
444 : }
445 : }
446 2786773 : for (j=1; j<co; j++) setlg(matb[j], lig-k0+1); /* bottom can be forgotten */
447 133918 : gerepileall(av, T? 2: 1, &matb, &T);
448 133924 : if (DEBUGLEVEL>5) err_printf(" matb cleaned up (using Id block)\n");
449 :
450 133924 : nlze = lk0 - k0; /* # of 0 rows */
451 133924 : lnz = lig-nlze+1; /* 1 + # of nonzero rows (!= 0...0 1 0 ... 0) */
452 133924 : if (T) matt = ZM_mul(matt,T); /* update top rows */
453 133922 : extramat = cgetg(col+1,t_MAT); /* = new C minus the 0 rows */
454 1377096 : for (j=1; j<=col; j++)
455 : {
456 1243176 : GEN z = gel(matt,j);
457 1243176 : GEN t = (gel(matb,j)) + nlze - k0;
458 1243176 : p2=cgetg(lnz,t_COL); gel(extramat,j) = p2;
459 5211839 : for (i=1; i<=k0; i++) gel(p2,i) = gel(z,i); /* top k0 rows */
460 1970556 : for ( ; i<lnz; i++) gel(p2,i) = gel(t,i); /* other nonzero rows */
461 : }
462 133920 : if (!col) {
463 0 : permpro = identity_perm(lnz);
464 0 : nr = lnz;
465 : }
466 : else
467 133920 : permpro = ZM_rowrankprofile(extramat, &nr);
468 : /* lnz = lg(permpro) */
469 133924 : if (nlze)
470 : { /* put the nlze 0 rows (trivial generators) at the top */
471 11499 : p1 = new_chunk(lk0+1);
472 38905 : for (i=1; i<=nlze; i++) p1[i] = perm[i + k0];
473 59871 : for ( ; i<=lk0; i++) p1[i] = perm[i - nlze];
474 87277 : for (i=1; i<=lk0; i++) perm[i] = p1[i];
475 : }
476 : /* sort other rows according to permpro (nr redundant generators first) */
477 133924 : p1 = new_chunk(lnz); p2 = perm + nlze;
478 586015 : for (i=1; i<lnz; i++) p1[i] = p2[permpro[i]];
479 586015 : 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 133925 : matbnew = cgetg(col+1,t_MAT); /* dense+toobig, maximal rank. For hnffinal */
495 133924 : dep = cgetg(col+1,t_MAT); /* rows dependent from the ones in matbnew */
496 1377109 : for (j=1; j<=col; j++)
497 : {
498 1243184 : GEN z = gel(extramat,j);
499 1243184 : p1 = cgetg(nlze+nr+1,t_COL); gel(dep,j) = p1;
500 1243181 : p2 = cgetg(lnz-nr,t_COL); gel(matbnew,j) = p2;
501 1586879 : for (i=1; i<=nlze; i++) gel(p1,i) = gen_0;
502 1311921 : p1 += nlze; for (i=1; i<=nr; i++) p1[i] = z[permpro[i]];
503 5870511 : 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 133925 : B = cgetg(co-col,t_MAT);
509 1543580 : for (j=col+1; j<co; j++)
510 : {
511 1409665 : GEN y = gel(matt,j);
512 1409665 : GEN z = gel(matb,j);
513 1409665 : p1=cgetg(lig+1,t_COL); gel(B,j-col) = p1;
514 3104928 : for (i=1; i<=nlze; i++) gel(p1,i) = gel(z,i);
515 1409655 : p1 += nlze; z += nlze-k0;
516 9853334 : for (k=1; k<lnz; k++)
517 : {
518 8443679 : i = permpro[k];
519 8443679 : gel(p1,k) = (i <= k0)? gel(y,i): gel(z,i);
520 : }
521 : }
522 133915 : if (T) C = typ(C)==t_MAT? RgM_ZM_mul(C,T): RgV_RgM_mul(C,T);
523 133915 : gerepileall(av, 4, &matbnew, &B, &dep, &C);
524 133925 : *ptdep = dep;
525 133925 : *ptB = B;
526 133925 : H = hnffinal(matbnew, perm, ptdep, ptB, &C);
527 133924 : 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 : gerepileall(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 133924 : *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 73189 : hnfadd_i(GEN H, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, /* cf hnfspec */
613 : GEN extramat,GEN extraC)
614 : {
615 73189 : GEN matb, extratop, Cnew, permpro, B = *ptB, C = *ptC, dep = *ptdep;
616 : long i, lH, lB, li, lig, co, col, nlze;
617 :
618 73189 : if (lg(extramat) == 1) return H;
619 73189 : co = lg(C)-1;
620 73189 : lH = lg(H)-1;
621 73189 : lB = lg(B)-1;
622 73189 : li = lg(perm)-1;
623 73189 : lig = li - lB;
624 73189 : col = co - lB;
625 73189 : nlze = lig - lH;
626 :
627 : /* col co
628 : * [ 0 |dep | ]
629 : * nlze [--------| B ]
630 : * [ 0 | H | ]
631 : * [--------|-----] lig
632 : * [ 0 | Id ]
633 : * [ | ] li */
634 73189 : extratop = zm_to_ZM( rowslicepermute(extramat, perm, 1, lig) );
635 73189 : if (li != lig)
636 : { /* zero out bottom part, using the Id block */
637 73175 : GEN A = vecslice(C, col+1, co);
638 73175 : GEN c = rowslicepermute(extramat, perm, lig+1, li);
639 73175 : extraC = gsub(extraC, typ(A)==t_MAT? RgM_zm_mul(A, c): RgV_zm_mul(A,c));
640 73175 : extratop = ZM_sub(extratop, ZM_zm_mul(B, c));
641 : }
642 :
643 73187 : extramat = shallowconcat(extratop, vconcat(dep, H));
644 73189 : Cnew = shallowconcat(extraC, vecslice(C, col-lH+1, co));
645 73189 : if (DEBUGLEVEL>5) err_printf(" 1st phase done\n");
646 73189 : permpro = ZM_rowrankprofile(extramat, &nlze);
647 73189 : extramat = rowpermute(extramat, permpro);
648 73189 : *ptB = rowpermute(B, permpro);
649 73189 : permpro = vecsmallpermute(perm, permpro);
650 266304 : for (i=1; i<=lig; i++) perm[i] = permpro[i]; /* perm o= permpro */
651 :
652 73189 : *ptdep = rowslice(extramat, 1, nlze);
653 73189 : matb = rowslice(extramat, nlze+1, lig);
654 73189 : if (DEBUGLEVEL>5) err_printf(" 2nd phase done\n");
655 73189 : H = hnffinal(matb,perm,ptdep,ptB,&Cnew);
656 73189 : *ptC = shallowconcat(vecslice(C, 1, col-lH), Cnew);
657 73189 : 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 36670043 : ZC_elem(GEN aj, GEN ak, GEN A, GEN U, long j, long k)
673 : {
674 : GEN p1,u,v,d;
675 :
676 36670043 : if (!signe(ak)) {
677 101748 : swap(gel(A,j), gel(A,k));
678 101748 : if (U) swap(gel(U,j), gel(U,k));
679 33285441 : return;
680 : }
681 36568295 : d = bezout(aj,ak,&u,&v);
682 : /* frequent special case (u,v) = (1,0) or (0,1) */
683 36560392 : if (!signe(u))
684 : { /* ak | aj */
685 17822108 : p1 = diviiexact(aj,ak); togglesign(p1);
686 17824443 : ZC_lincomb1_inplace(gel(A,j), gel(A,k), p1);
687 17833019 : if (U)
688 1791885 : ZC_lincomb1_inplace(gel(U,j), gel(U,k), p1);
689 17832134 : return;
690 : }
691 18738284 : if (!signe(v))
692 : { /* aj | ak */
693 15353440 : p1 = diviiexact(ak,aj); togglesign(p1);
694 15351627 : ZC_lincomb1_inplace(gel(A,k), gel(A,j), p1);
695 15351537 : swap(gel(A,j), gel(A,k));
696 15351537 : if (U) {
697 333979 : ZC_lincomb1_inplace(gel(U,k), gel(U,j), p1);
698 334001 : swap(gel(U,j), gel(U,k));
699 : }
700 15351559 : return;
701 : }
702 :
703 3384844 : if (!is_pm1(d)) { aj = diviiexact(aj, d); ak = diviiexact(ak, d); }
704 3384828 : p1 = gel(A,k); aj = negi(aj); /* NOT togglesign */
705 3385136 : gel(A,k) = ZC_lincomb(u,v, gel(A,j),p1);
706 3384680 : gel(A,j) = ZC_lincomb(aj,ak, p1,gel(A,j));
707 3384939 : if (U)
708 : {
709 646038 : p1 = gel(U,k);
710 646038 : gel(U,k) = ZC_lincomb(u,v, gel(U,j),p1);
711 646034 : 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;
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 : if (degpol(d)) { a = RgX_div(a, d); b = RgX_div(b, d); }
731 238 : else if (typ(gel(d,2)) == t_REAL && lg(gel(d,2)) <= 3)
732 : #if 1
733 : { /* possible accuracy problem */
734 0 : GEN D = RgX_gcd_simple(a,b);
735 0 : if (degpol(D)) {
736 0 : D = RgX_normalize(D);
737 0 : a = RgX_div(a, D);
738 0 : b = RgX_div(b, D);
739 0 : d = RgX_extgcd(a,b, pu,pv); /* retry now */
740 0 : d = RgX_mul(d, D);
741 : }
742 : }
743 : #else
744 : { /* less stable */
745 : d = RgX_extgcd_simple(a,b, pu,pv);
746 : if (degpol(d)) { a = RgX_div(a, d); b = RgX_div(b, d); }
747 : }
748 : #endif
749 882 : *pa = a;
750 882 : *pb = b; return d;
751 : }
752 : static GEN
753 2694051 : col_mul(GEN x, GEN c)
754 : {
755 2694051 : if (typ(x) == t_INT)
756 : {
757 2692109 : long s = signe(x);
758 2692109 : if (!s) return NULL;
759 2041064 : if (is_pm1(x)) return (s > 0)? c: RgC_neg(c);
760 : }
761 309887 : return RgC_Rg_mul(c, x);
762 : }
763 : static void
764 0 : do_zero(GEN x)
765 : {
766 0 : long i, lx = lg(x);
767 0 : for (i=1; i<lx; i++) gel(x,i) = gen_0;
768 0 : }
769 :
770 : /* (c1, c2) *= [u,-b; v,a] */
771 : static void
772 673507 : update(GEN u, GEN v, GEN a, GEN b, GEN *c1, GEN *c2)
773 : {
774 : GEN p1,p2;
775 :
776 673507 : u = col_mul(u,*c1);
777 673589 : v = col_mul(v,*c2);
778 673597 : if (u) p1 = v? gadd(u,v): u;
779 6436 : else p1 = v? v: NULL;
780 :
781 673593 : a = col_mul(a,*c2);
782 673611 : b = col_mul(gneg_i(b),*c1);
783 673587 : if (a) p2 = b? RgC_add(a,b): a;
784 0 : else p2 = b? b: NULL;
785 :
786 673497 : if (!p1) do_zero(*c1); else *c1 = p1;
787 673498 : if (!p2) do_zero(*c2); else *c2 = p2;
788 673498 : }
789 :
790 : /* zero aj = Aij (!= 0) using ak = Aik (maybe 0), via linear combination of
791 : * A[j] and A[k] of determinant 1. If U != NULL, likewise update its columns */
792 : static void
793 511 : RgC_elem(GEN aj, GEN ak, GEN A, GEN V, long j, long k, long li, long vx)
794 : {
795 511 : GEN u,v, d = gbezout_step(&aj, &ak, &u, &v, vx);
796 : long l;
797 : /* (A[,k], A[,j]) *= [v, -aj; u, ak ] */
798 1778 : for (l = 1; l < li; l++)
799 : {
800 1267 : GEN t = gadd(gmul(u,gcoeff(A,l,j)), gmul(v,gcoeff(A,l,k)));
801 1267 : gcoeff(A,l,j) = gsub(gmul(ak,gcoeff(A,l,j)), gmul(aj,gcoeff(A,l,k)));
802 1267 : gcoeff(A,l,k) = t;
803 : }
804 511 : gcoeff(A,li,j) = gen_0;
805 511 : gcoeff(A,li,k) = d;
806 511 : if (V) update(v,u,ak,aj,(GEN*)(V+k),(GEN*)(V+j));
807 511 : }
808 :
809 : /* reduce A[i,j] mod A[i,j0] for j=j0+1... via column operations */
810 : static void
811 4153461 : ZM_reduce(GEN A, GEN U, long i, long j0)
812 : {
813 4153461 : long j, lA = lg(A);
814 4153461 : GEN d = gcoeff(A,i,j0);
815 4153461 : if (signe(d) < 0)
816 : {
817 28234 : ZV_neg_inplace(gel(A,j0));
818 28234 : if (U) ZV_togglesign(gel(U,j0));
819 28234 : d = gcoeff(A,i,j0);
820 : }
821 8353333 : for (j=j0+1; j<lA; j++)
822 : {
823 4199839 : GEN q = truedivii(gcoeff(A,i,j), d);
824 4199824 : if (!signe(q)) continue;
825 :
826 257763 : togglesign(q);
827 257864 : ZC_lincomb1_inplace(gel(A,j), gel(A,j0), q);
828 257867 : if (U) ZC_lincomb1_inplace(gel(U,j), gel(U,j0), q);
829 : }
830 4153494 : }
831 :
832 : /* normalize T as if it were a t_POL in variable v */
833 : static GEN
834 364 : normalize_as_RgX(GEN T, long v, GEN *pd)
835 : {
836 : GEN d;
837 364 : if (!is_RgX(T,v)) { *pd = T; return gen_1; }
838 336 : d = leading_coeff(T);
839 336 : while (gequal0(d) || (typ(d) == t_REAL && lg(d) == 3
840 0 : && gexpo(T) - expo(d) > (long)BITS_IN_LONG)) {
841 14 : T = normalizepol_lg(T, lg(T)-1);
842 14 : if (!signe(T)) { *pd = gen_1; return T; }
843 0 : d = leading_coeff(T);
844 : }
845 322 : if (degpol(T)) T = RgX_Rg_div(T,d); else { d = gel(T,2); T = gen_1; }
846 322 : *pd = d; return T;
847 : }
848 : /* reduce A[i,j] mod A[i,j0] for j=j0+1... via column operations */
849 : static void
850 84 : RgM_reduce(GEN A, GEN U, long i, long j0, long vx)
851 : {
852 84 : long j, lA = lg(A);
853 84 : GEN d, T = normalize_as_RgX(gcoeff(A,i,j0), vx, &d);
854 84 : if (U && !gequal1(d)) gel(U,j0) = RgC_Rg_div(gel(U,j0), d);
855 84 : gcoeff(A,i,j0) = T;
856 :
857 196 : for (j=j0+1; j<lA; j++)
858 : {
859 112 : GEN t = gcoeff(A,i,j), q;
860 112 : if (gequal0(t)) continue;
861 14 : if (T == gen_1)
862 0 : q = t;
863 14 : else if (is_RgX(t,vx))
864 14 : q = RgX_div(t, T);
865 0 : else continue;
866 :
867 14 : if (gequal0(q)) continue;
868 7 : gel(A,j) = RgC_sub(gel(A,j), RgC_Rg_mul(gel(A,j0), q));
869 7 : if (U) gel(U,j) = RgC_sub(gel(U,j), RgC_Rg_mul(gel(U,j0), q));
870 : }
871 84 : }
872 :
873 : /* A,B square integral in upper HNF, of the same dimension > 0. Return Au
874 : * in Z^n (v in Z^n not computed), such that Au + Bv = [1, 0, ..., 0] */
875 : GEN
876 637288 : hnfmerge_get_1(GEN A, GEN B)
877 : {
878 637288 : pari_sp av = avma;
879 637288 : long j, k, l = lg(A), lb;
880 637288 : GEN b, U = cgetg(l + 1, t_MAT), C = cgetg(l + 1, t_VEC);
881 :
882 637328 : b = gcoeff(B,1,1); lb = lgefint(b);
883 1202795 : for (j = 1; j < l; j++)
884 : {
885 : GEN t;
886 1202793 : long c = j+1;
887 1202793 : gel(U,j) = col_ei(l-1, j);
888 1202974 : gel(U,c) = zerocol(l-1); /* dummy */
889 1203010 : gel(C,j) = vecslice(gel(A,j), 1,j);
890 1203010 : gel(C,c) = vecslice(gel(B,j), 1,j);
891 3278582 : for (k = j; k > 0; k--)
892 : {
893 2076121 : t = gcoeff(C,k,c);
894 2076121 : if (gequal0(t)) continue;
895 1869825 : setlg(C[c], k+1);
896 1869798 : ZC_elem(t, gcoeff(C,k,k), C, U, c, k);
897 1868716 : if (lgefint(gcoeff(C,k,k)) > lb) gel(C,k) = FpC_red(gel(C,k), b);
898 1868716 : if (j > 4)
899 : {
900 89125 : GEN u = gel(U,k);
901 : long h;
902 1332695 : for (h=1; h<l; h++)
903 1243570 : if (lgefint(gel(u,h)) > lb) gel(u,h) = remii(gel(u,h), b);
904 : }
905 : }
906 1202461 : if (j == 1)
907 637077 : t = gcoeff(C,1,1);
908 : else
909 : {
910 : GEN u;
911 565384 : t = bezout(gcoeff(C,1,1), b, &u, NULL); /* >= 0 */
912 565785 : if (signe(u) && !equali1(u)) gel(U,1) = ZC_Z_mul(gel(U,1), u);
913 565794 : gcoeff(C,1,1) = t;
914 : }
915 1202871 : if (equali1(t)) break;
916 : }
917 637237 : if (j >= l) return NULL;
918 637237 : b = lcmii(gcoeff(A,1,1),b);
919 637220 : A = FpC_red(ZM_ZC_mul(A,gel(U,1)), b);
920 637042 : return gerepileupto(av, FpC_center(A, b, shifti(b,-1)));
921 : }
922 :
923 : /* remove the first r columns */
924 : static void
925 308867 : remove_0cols(long r, GEN *pA, GEN *pB, long remove)
926 : {
927 308867 : GEN A = *pA, B = *pB;
928 308867 : long l = lg(A);
929 308867 : A += r; A[0] = evaltyp(t_MAT) | evallg(l-r);
930 308865 : if (B && remove == 2) { B += r; B[0] = A[0]; }
931 308865 : *pA = A; *pB = B;
932 308865 : }
933 :
934 : /* Inefficient compared to hnfall. 'remove' = throw away lin.dep columns */
935 : static GEN
936 54482 : hnf_i(GEN A, int remove)
937 : {
938 54482 : pari_sp av0 = avma, av;
939 : long s, n, m, j, k, li, def, ldef;
940 :
941 54482 : RgM_dimensions(A, &m, &n);
942 54482 : if (!n) return cgetg(1,t_MAT);
943 54230 : av = avma;
944 54230 : A = RgM_shallowcopy(A);
945 54230 : def = n; ldef = (m>n)? m-n: 0;
946 147584 : for (li=m; li>ldef; li--)
947 : {
948 384643 : for (j=def-1; j; j--)
949 : {
950 291289 : GEN a = gcoeff(A,li,j);
951 291289 : if (!signe(a)) continue;
952 :
953 : /* zero a = Aij using b = Aik */
954 150967 : k = (j==1)? def: j-1;
955 150967 : ZC_elem(a,gcoeff(A,li,k), A,NULL, j,k);
956 150967 : if (gc_needed(av,1))
957 : {
958 0 : if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnf[1]. li=%ld",li);
959 0 : A = gerepilecopy(av, A);
960 : }
961 : }
962 93354 : s = signe(gcoeff(A,li,def));
963 93354 : if (s)
964 : {
965 92108 : if (s < 0) ZV_neg_inplace(gel(A,def));
966 92108 : ZM_reduce(A, NULL, li,def);
967 92108 : def--;
968 : }
969 : else
970 1246 : if (ldef) ldef--;
971 93354 : if (gc_needed(av,1))
972 : {
973 0 : if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnf[2]. li=%ld",li);
974 0 : A = gerepilecopy(av, A);
975 : }
976 : }
977 : /* rank A = n - def */
978 54230 : if (remove) { GEN B = NULL; remove_0cols(def, &A, &B, remove); }
979 54230 : return gerepileupto(av0, ZM_copy(A));
980 : }
981 :
982 : GEN
983 76681 : ZM_hnf(GEN x) { return lg(x) > 8? ZM_hnfall(x, NULL, 1): hnf_i(x, 1); }
984 :
985 : /* u*z[1..k] mod p, in place */
986 : static void
987 4703834 : FpV_Fp_mul_part_ip(GEN z, GEN u, GEN p, long k)
988 : {
989 : long i;
990 4703834 : if (is_pm1(u)) {
991 191976 : if (signe(u) > 0) {
992 570385 : for (i = 1; i <= k; i++)
993 431721 : if (signe(gel(z,i))) gel(z,i) = modii(gel(z,i), p);
994 : } else {
995 152890 : for (i = 1; i <= k; i++)
996 99577 : if (signe(gel(z,i))) gel(z,i) = modii(negi(gel(z,i)), p);
997 : }
998 : }
999 : else {
1000 17835589 : for (i = 1; i <= k; i++)
1001 13324175 : if (signe(gel(z,i))) gel(z,i) = Fp_mul(u,gel(z,i), p);
1002 : }
1003 4703391 : }
1004 : static void
1005 33003933 : FpV_red_part_ipvec(GEN z, GEN p, long k)
1006 : {
1007 : long i;
1008 100990400 : for (i = 1; i <= k; i++) gel(z,i) = modii(gel(z,i), gel(p,i));
1009 32998573 : }
1010 :
1011 : /* return x * U, in echelon form (mod p^m), where (det(U),p) = 1.
1012 : * If early_abort is set, return NULL as soon as one pivot is 0 (mod p^m) */
1013 : GEN
1014 522875 : ZpM_echelon(GEN x, long early_abort, GEN p, GEN pm)
1015 : {
1016 522875 : pari_sp av0 = avma, av;
1017 : long m, li, co, i, j, k, def, ldef;
1018 :
1019 522875 : co = lg(x); if (co == 1) return cgetg(1,t_MAT);
1020 522875 : li = lgcols(x);
1021 522874 : av = avma;
1022 522874 : x = RgM_shallowcopy(x);
1023 522877 : m = Z_pval(pm, p);
1024 :
1025 522875 : ldef = (li > co)? li - co: 0;
1026 2495570 : for (def = co-1,i = li-1; i > ldef; i--)
1027 : {
1028 1973338 : long vmin = LONG_MAX, kmin = 0;
1029 1973338 : GEN umin = gen_0, pvmin, q;
1030 9700939 : for (k = 1; k <= def; k++)
1031 : {
1032 8186824 : GEN u = gcoeff(x,i,k);
1033 : long v;
1034 8186824 : if (!signe(u)) continue;
1035 3856673 : v = Z_pvalrem(u, p, &u);
1036 3856533 : if (v >= m) gcoeff(x,i,k) = gen_0;
1037 2715682 : else if (v < vmin) {
1038 1390598 : vmin = v; kmin = k; umin = u;
1039 1390598 : if (!vmin) break;
1040 : }
1041 : }
1042 1973198 : if (!kmin)
1043 : {
1044 925712 : if (early_abort) return NULL;
1045 925022 : gcoeff(x,i,def) = gen_0;
1046 925022 : ldef--;
1047 925022 : if (ldef < 0) ldef = 0;
1048 925022 : continue;
1049 : }
1050 1047486 : if (kmin != def) swap(gel(x,def), gel(x,kmin));
1051 1047486 : q = vmin? powiu(p, m-vmin): pm;
1052 : /* pivot has valuation vmin */
1053 1047473 : umin = modii(umin, q);
1054 1047569 : if (!equali1(umin))
1055 525037 : FpV_Fp_mul_part_ip(gel(x,def), Fp_inv(umin,q), pm, i-1);
1056 1047578 : gcoeff(x, i, def) = pvmin = powiu(p, vmin);
1057 7117489 : for (j = def-1; j; j--)
1058 : { /* zero x[i, 1..def-1] using x[i,def] = pvmin */
1059 6069816 : GEN t, a = gcoeff(x,i,j) = modii(gcoeff(x,i,j), pm);
1060 6069285 : if (!signe(a)) continue;
1061 :
1062 2839337 : t = diviiexact(a, pvmin); togglesign(t);
1063 2839952 : ZC_lincomb1_inplace(gel(x,j), gel(x,def), t);
1064 2839954 : if (gc_needed(av,1))
1065 : {
1066 14 : if (DEBUGMEM>1) pari_warn(warnmem,"ZpM_echelon. i=%ld",i);
1067 14 : x = gerepilecopy(av, x); pvmin = gcoeff(x,i,def);
1068 : }
1069 : }
1070 1047673 : def--;
1071 : }
1072 522232 : if (co > li)
1073 : {
1074 0 : x += co - li;
1075 0 : x[0] = evaltyp(t_MAT) | evallg(li);
1076 : }
1077 522232 : return gerepilecopy(av0, x);
1078 : }
1079 : GEN
1080 4469572 : zlm_echelon(GEN x, long early_abort, ulong p, ulong pm)
1081 : {
1082 4469572 : pari_sp av0 = avma;
1083 : long li, co, i, j, k, def, ldef;
1084 : ulong m;
1085 :
1086 4469572 : co = lg(x); if (co == 1) return cgetg(1,t_MAT);
1087 4469572 : li = lgcols(x);
1088 4469630 : x = Flm_copy(x);
1089 4469682 : m = u_lval(pm, p);
1090 :
1091 4469687 : ldef = (li > co)? li - co: 0;
1092 20012171 : for (def = co-1,i = li-1; i > ldef; i--)
1093 : {
1094 15727209 : long vmin = LONG_MAX, kmin = 0;
1095 15727209 : ulong umin = 0, pvmin, q;
1096 48217990 : for (k = 1; k <= def; k++)
1097 : {
1098 40834294 : ulong u = ucoeff(x,i,k);
1099 : long v;
1100 40834294 : if (!u) continue;
1101 22168004 : v = u_lvalrem(u, p, &u);
1102 22167934 : if (v >= (long) m) ucoeff(x,i,k) = 0;
1103 22167934 : else if (v < vmin) {
1104 15776868 : vmin = v; kmin = k; umin = u;
1105 15776868 : if (!vmin) break;
1106 : }
1107 : }
1108 15727139 : if (!kmin)
1109 : {
1110 1014373 : if (early_abort) return NULL;
1111 829982 : ucoeff(x,i,def) = 0;
1112 829982 : ldef--;
1113 829982 : if (ldef < 0) ldef = 0;
1114 829982 : continue;
1115 : }
1116 14712766 : if (kmin != def) swap(gel(x,def), gel(x,kmin));
1117 14712766 : q = vmin? upowuu(p, m-vmin): pm;
1118 : /* pivot has valuation vmin */
1119 14712828 : umin %= q;
1120 14712828 : if (umin != 1)
1121 7539309 : Flv_Fl_mul_part_inplace(gel(x,def), Fl_inv(umin,q), pm, i-1);
1122 14712774 : ucoeff(x, i, def) = pvmin = upowuu(p, vmin);
1123 66553234 : for (j = def-1; j; j--)
1124 : { /* zero x[i, 1..def-1] using x[i,def] = pvmin */
1125 51840732 : ulong t, a = ucoeff(x,i,j);
1126 51840732 : if (!a) continue;
1127 :
1128 29271348 : t = Fl_neg(a / pvmin, q);
1129 29271317 : Flc_lincomb1_inplace(gel(x,j), gel(x,def), t, pm);
1130 : }
1131 14712502 : def--;
1132 : }
1133 4284962 : if (co > li)
1134 : {
1135 0 : x += co - li;
1136 0 : x[0] = evaltyp(t_MAT) | evallg(li);
1137 : }
1138 4284962 : return gerepilecopy(av0, x);
1139 : }
1140 :
1141 : static int
1142 219589 : ZV_allequal(GEN v)
1143 : {
1144 219589 : long i, l = lg(v);
1145 219589 : if (l > 1)
1146 : {
1147 219589 : GEN x = gel(v,1);
1148 388861 : for (i = 2; i < l; i++) if (!equalii(x,gel(v,i))) return 0;
1149 : }
1150 198898 : return 1;
1151 : }
1152 : /* compute optimal D for hnfmod: x upper triangular */
1153 : static GEN
1154 4456086 : optimal_D(GEN x, GEN D)
1155 : {
1156 4456086 : long i, n = nbrows(x);
1157 4456295 : GEN C = shallowcopy(D);
1158 4457552 : gel(C,1) = gcoeff(x,1,1);
1159 4910526 : for (i = 2; i < n; i++)
1160 : {
1161 2296031 : GEN c = mulii(gel(C,i-1), gcoeff(x,i,i));
1162 2295757 : if (signe(c) < 0) togglesign(c);
1163 2295759 : if (cmpii(c, gel(D,i)) >= 0) break;
1164 452974 : gel(C,i) = c;
1165 : }
1166 4457445 : return C;
1167 : }
1168 :
1169 : /* D = multiple of det x (usually detint(x)) or vector of positive moduli
1170 : * (compute hnf(x | D))
1171 : * flag & hnf_MODID: reduce mod D * matid [ otherwise as above ].
1172 : * flag & hnf_PART: don't reduce once diagonal is known
1173 : * flag & hnf_CENTER: centermod HNF (2|x[i,j]| <] x[i,i]) */
1174 : GEN
1175 4459178 : ZM_hnfmodall_i(GEN x, GEN D, long flag)
1176 : {
1177 : pari_sp av;
1178 4459178 : const long center = (flag & hnf_CENTER);
1179 : long moddiag, modsame, nli, li, co, i, j, k, def, ldef;
1180 : GEN u, LDM;
1181 :
1182 4459178 : co = lg(x);
1183 4459178 : if (co == 1)
1184 : {
1185 1554 : if (typ(D) == t_INT || lg(D) == 1) return cgetg(1,t_MAT);
1186 105 : x = diagonal_shallow(D); /* handle flags properly */
1187 105 : co = lg(x);
1188 : }
1189 4457729 : li = lgcols(x);
1190 4457949 : if (li == 1)
1191 : {
1192 131 : if (typ(D) != t_INT && lg(D) != li) pari_err_DIM("ZM_hnfmod");
1193 131 : return cgetg(1,t_MAT);
1194 : }
1195 4457854 : nli = li - 1;
1196 4457854 : modsame = typ(D)==t_INT;
1197 4457854 : if (!modsame)
1198 : {
1199 219603 : if (lg(D) != li) pari_err_DIM("ZM_hnfmod");
1200 219589 : if (ZV_allequal(D)) { modsame = 1; D = gel(D,1); }
1201 : }
1202 4457841 : moddiag = (flag & hnf_MODID) || !modsame;
1203 : /* modsame: triangularize mod fixed d*Id;
1204 : * moddiag: modulo diagonal matrix, else modulo multiple of determinant */
1205 :
1206 4457841 : if (modsame)
1207 : {
1208 4437131 : LDM = const_vecsmall(nli, 2*lgefint(D)-2);
1209 4437240 : D = const_vec(nli,D);
1210 : }
1211 : else
1212 : {
1213 20710 : LDM = cgetg(li, t_VECSMALL);
1214 68039 : for (i=1; i<li; i++) LDM[i] = lgefint(gel(D,i));
1215 : }
1216 4458051 : av = avma;
1217 4458051 : x = RgM_shallowcopy(x);
1218 :
1219 4458164 : ldef = 0;
1220 4458164 : if (li > co)
1221 : {
1222 135054 : ldef = li - co;
1223 135054 : if (!moddiag)
1224 7 : pari_err_DOMAIN("ZM_hnfmod","nb lines",">", strtoGENstr("nb columns"), x);
1225 : }
1226 17736827 : for (def = co-1,i = nli; i > ldef; i--,def--)
1227 : {
1228 13276905 : GEN d = gel(D,i);
1229 13276905 : long add_N = modsame;
1230 58952160 : for (j = 1; j < def; j++)
1231 : {
1232 45673500 : GEN p1, p2, b, a = gcoeff(x,i,j) = remii(gcoeff(x,i,j), d);
1233 57668123 : if (!signe(a)) continue;
1234 :
1235 31209801 : k = j+1;
1236 31209801 : b = gcoeff(x,i,k) = remii(gcoeff(x,i,k), d);
1237 31216215 : if (!signe(b)) { swap(gel(x,j), gel(x,k)); continue; }
1238 19210260 : if (add_N)
1239 : { /* ensure the moving pivot on row i divides d from now on */
1240 6725688 : add_N = 0;
1241 6725688 : if (!equali1(a))
1242 : { /* x[j] *= u; after this, a = x[i,j] | d */
1243 5594066 : GEN u = Fp_invgen(a, d, &a);
1244 : long t;
1245 5595549 : p1 = gel(x,j);
1246 19946442 : for (t = 1; t < i; t++) gel(p1,t) = mulii(gel(p1,t), u);
1247 5593659 : FpV_red_part_ipvec(p1, D, i-1);
1248 5593637 : gel(p1,i) = a;
1249 5593637 : if (2*lg(a) < lg(b))
1250 : { /* reduce x[i,k] mod x[i,j]: helps ZC_elem */
1251 26164 : GEN r, q = dvmdii(b, a, &r);
1252 26164 : togglesign(q);
1253 26164 : ZC_lincomb1_inplace_i(gel(x,k), gel(x,j), q, i-1);
1254 26164 : FpV_red_part_ipvec(gel(x,k), D, i-1);
1255 26164 : gcoeff(x,i,k) = b = r;
1256 : }
1257 : }
1258 : }
1259 19209819 : ZC_elem(a,b, x, NULL, j,k);
1260 19216933 : p1 = gel(x,j);
1261 19216933 : p2 = gel(x,k);
1262 : /* prevent coeffs explosion */
1263 113084491 : for (k = 1; k < i; k++)
1264 : {
1265 93867558 : if (lgefint(gel(p1,k)) > LDM[k]) gel(p1,k) = remii(gel(p1,k),gel(D,k));
1266 93867558 : if (lgefint(gel(p2,k)) > LDM[k]) gel(p2,k) = remii(gel(p2,k),gel(D,k));
1267 : }
1268 : }
1269 13278660 : if (gc_needed(av,2))
1270 : {
1271 25 : if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnfmod[1]. i=%ld",i);
1272 25 : x = gerepilecopy(av, x);
1273 : }
1274 13278660 : if (moddiag && !signe(gcoeff(x,i,def)))
1275 : { /* missing pivot on line i, insert column */
1276 903981 : GEN a = cgetg(co + 1, t_MAT);
1277 4876704 : for (k = 1; k <= def; k++) gel(a,k) = gel(x,k);
1278 903983 : gel(a,k++) = Rg_col_ei(gel(D,i), nli, i);
1279 3931713 : for ( ; k <= co; k++) gel(a,k) = gel(x,k-1);
1280 903989 : ldef--; if (ldef < 0) ldef = 0;
1281 903989 : co++; def++; x = a;
1282 : }
1283 : }
1284 4459922 : if (co < li)
1285 : { /* implies moddiag, add missing diag(D) components */
1286 103013 : GEN a = cgetg(li+1, t_MAT);
1287 242049 : for (k = 1; k <= li-co; k++) gel(a,k) = Rg_col_ei(gel(D,k), nli, k);
1288 282857 : for (i = 1; i < co; i++) gel(a,k-1+i) = gel(x,i);
1289 103012 : gel(a,li) = zerocol(nli); x = a;
1290 : }
1291 : else
1292 : {
1293 4356909 : x += co - li;
1294 4356909 : x[0] = evaltyp(t_MAT) | evallg(li); /* kill 0 columns */
1295 4354709 : if (moddiag) x = shallowconcat(x, zerocol(nli));
1296 : }
1297 4458097 : if (moddiag)
1298 : { /* x[li]: extra column, an accumulator discarded at the end */
1299 : GEN D2;
1300 4303059 : gcoeff(x,1,1) = gcdii(gcoeff(x,1,1), gel(D,1));
1301 4301138 : D2 = optimal_D(x,D);
1302 : /* add up missing diag(D) components */
1303 17287372 : for (i = nli; i > 0; i--)
1304 : {
1305 12984947 : gcoeff(x, i, li) = gel(D,i);
1306 50629298 : for (j = i; j > 0; j--)
1307 : {
1308 37647131 : GEN a = gcoeff(x, j, li);
1309 37647131 : if (!signe(a)) continue;
1310 13704448 : ZC_elem(a, gcoeff(x,j,j), x, NULL, li,j);
1311 13702203 : FpV_red_part_ipvec(gel(x,li), D, j-1);
1312 13701879 : FpV_red_part_ipvec(gel(x,j), D, j-1);
1313 : }
1314 12982167 : if (gc_needed(av,1))
1315 : {
1316 26 : if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnfmod[2]. i=%ld", i);
1317 26 : gerepileall(av, 2, &x, &D2);
1318 : }
1319 : }
1320 4302425 : D = D2;
1321 : }
1322 : else
1323 : {
1324 155038 : GEN b = gel(D,1);
1325 582661 : for (i = nli; i > 0; i--)
1326 : {
1327 427625 : GEN d = bezout(gcoeff(x,i,i),b, &u,NULL);
1328 427634 : gcoeff(x,i,i) = d;
1329 427634 : FpV_Fp_mul_part_ip(gel(x,i), u, b, i-1);
1330 427634 : if (i > 1) b = diviiexact(b,d);
1331 : }
1332 155036 : D = optimal_D(x,D);
1333 : }
1334 4457461 : x[0] = evaltyp(t_MAT) | evallg(li); /* kill 0 columns / discard accumulator */
1335 4457483 : if (flag & hnf_PART) return x;
1336 :
1337 17869651 : for (i = nli; i > 0; i--)
1338 : {
1339 13415905 : GEN diag = gcoeff(x,i,i);
1340 13415905 : if (signe(diag) < 0) { gel(x,i) = ZC_neg(gel(x,i)); diag = gcoeff(x,i,i); }
1341 13415910 : if (i != nli)
1342 25293114 : for (j = 1; j < i; j++) gcoeff(x,j,i) = remii(gcoeff(x,j,i), gel(D,j));
1343 38705341 : for (j = i+1; j < li; j++)
1344 : {
1345 25291215 : GEN b = gcoeff(x,i,j) = remii(gcoeff(x,i,j), gel(D,i));
1346 25289326 : GEN r, q = truedvmdii(b, diag, &r);
1347 : /* ensure -diag/2 <= r < diag/2 */
1348 25289054 : if (center && signe(r) && abscmpii(shifti(r,1),diag) >= 0)
1349 626464 : { r = subii(r,diag); q = addiu(q,1); }
1350 25291327 : if (!signe(q)) continue;
1351 10408197 : togglesign(q);
1352 10408201 : ZC_lincomb1_inplace_i(gel(x,j), gel(x,i), q, i-1);
1353 10406501 : gcoeff(x,i,j) = r;
1354 : }
1355 13414126 : if (gc_needed(av,1))
1356 : {
1357 39 : if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnfmod[3]. i=%ld", i);
1358 39 : gerepileall(av, 2, &x, &D);
1359 : }
1360 : }
1361 4453746 : return x;
1362 : }
1363 : GEN
1364 4409618 : ZM_hnfmodall(GEN x, GEN dm, long flag)
1365 : {
1366 4409618 : pari_sp av = avma;
1367 4409618 : return gerepilecopy(av, ZM_hnfmodall_i(x, dm, flag));
1368 : }
1369 : GEN
1370 155031 : ZM_hnfmod(GEN x, GEN d) { return ZM_hnfmodall(x,d,0); }
1371 : GEN
1372 3856178 : ZM_hnfmodid(GEN x, GEN d)
1373 3856178 : { return ZM_hnfmodall(x,d,hnf_MODID); }
1374 : /* return the column echelon form of x with 1's as pivots,
1375 : * P contains the row indices containing the pivots in increasing order */
1376 : static GEN
1377 2754298 : FpM_echelon(GEN x, GEN *pP, GEN p)
1378 : {
1379 : pari_sp av;
1380 : long iP, li, co, i, j, k, def, ldef;
1381 : GEN P;
1382 :
1383 2754298 : co = lg(x); if (co == 1) { *pP = cgetg(1,t_VECSMALL); return cgetg(1,t_MAT); }
1384 2754298 : li = lgcols(x);
1385 2754285 : iP = 1;
1386 2754285 : *pP = P = cgetg(li, t_VECSMALL);
1387 2754287 : av = avma;
1388 2754287 : x = FpM_red(x, p);
1389 :
1390 2753541 : ldef = (li > co)? li - co: 0;
1391 10402566 : for (def = co-1,i = li-1; i > ldef; i--)
1392 : {
1393 7646600 : GEN u = NULL;
1394 11791698 : for (k = def; k; k--)
1395 : {
1396 8842988 : u = gcoeff(x,i,k);
1397 8842988 : if (signe(u)) break;
1398 : }
1399 7646600 : if (!k)
1400 : {
1401 2950626 : if (--ldef < 0) ldef = 0;
1402 2950626 : continue;
1403 : }
1404 4695974 : P[iP++] = i;
1405 4695974 : if (k != def) swap(gel(x,def), gel(x,k));
1406 4695974 : if (!equali1(u))
1407 3751514 : FpV_Fp_mul_part_ip(gel(x,def), Fp_inv(u,p), p, i-1);
1408 4697920 : gcoeff(x, i, def) = gen_1;
1409 15420437 : for (j = def-1; j; j--)
1410 : { /* zero x[i, 1..def-1] using x[i,def] = 1*/
1411 10722710 : GEN xj = gel(x,j), u = gel(xj,i);
1412 10722710 : if (!signe(u)) continue;
1413 :
1414 6759608 : ZC_lincomb1_inplace(xj, gel(x,def), negi(u));
1415 37558097 : for (k = 1; k < i; k++) gel(xj,k) = modii(gel(xj,k), p);
1416 : }
1417 4697727 : if (gc_needed(av,2))
1418 : {
1419 0 : if (DEBUGMEM>1) pari_warn(warnmem,"FpM_echelon. i=%ld",i);
1420 0 : x = gerepilecopy(av, x);
1421 : }
1422 4698399 : def--;
1423 : }
1424 : /* rank = iP - 1 */
1425 2755966 : setlg(P, iP); vecsmall_sort(P);
1426 2754339 : if (co > iP) x += co - iP;
1427 2754339 : x[0] = evaltyp(t_MAT) | evallg(iP);
1428 2754322 : return x;
1429 : }
1430 : /* given x square of maximal rank with 1 or p on diagonal from hnfmodid
1431 : * (=> a column containing p has its other entries at 0 ), return the HNF */
1432 : static GEN
1433 2754503 : FpM_hnfend(pari_sp av, GEN x, GEN p)
1434 : {
1435 2754503 : long i, l = lgcols(x);
1436 10402458 : for (i = l-1; i > 0; i--)
1437 : {
1438 7649635 : GEN diag = gcoeff(x,i,i);
1439 : long j;
1440 7649635 : if (is_pm1(diag))
1441 10156965 : for (j = i+1; j < l; j++)
1442 : {
1443 5457977 : GEN xj = gel(x,j), b = gel(xj,i);
1444 : long k;
1445 5457977 : if (!signe(b)) continue;
1446 4158711 : ZC_lincomb1_inplace(xj, gel(x,i), negi(b));
1447 20170626 : for (k=1; k<i; k++)
1448 16012132 : if (lgefint(gel(xj,k)) > 3) gel(xj,k) = remii(gel(xj,k), p);
1449 : }
1450 : else
1451 8488319 : for (j = i+1; j < l; j++) gcoeff(x,i,j) = modii(gcoeff(x,i,j), p);
1452 7648882 : if (gc_needed(av,2))
1453 : {
1454 0 : if (DEBUGMEM>1) pari_warn(warnmem,"FpM_hnfend. i=%ld",i);
1455 0 : x = gerepilecopy(av, x);
1456 : }
1457 : }
1458 2752823 : return x;
1459 : }
1460 : GEN
1461 2754296 : ZM_hnfmodprime(GEN x, GEN p)
1462 : {
1463 2754296 : pari_sp av = avma;
1464 : GEN P, y;
1465 : long l, lP, i;
1466 2754296 : if (lg(x) == 1) return cgetg(1, t_MAT);
1467 2754296 : l = lgcols(x);
1468 2754285 : x = FpM_echelon(x, &P, p);
1469 2754321 : lP = lg(P); /* rank = lP-1 */
1470 2754321 : if (lP == l) { set_avma(av); return matid(l-1); }
1471 2754321 : y = scalarmat_shallow(p, l-1);
1472 7453788 : for (i = 1; i < lP; i++) gel(y,P[i]) = gel(x,i);
1473 2754508 : return gerepilecopy(av, FpM_hnfend(av,y,p));
1474 : }
1475 :
1476 : static GEN
1477 397172 : allhnfmod(GEN x, GEN dm, int flag)
1478 : {
1479 397172 : if (typ(x)!=t_MAT) pari_err_TYPE("allhnfmod",x);
1480 397172 : RgM_check_ZM(x, "allhnfmod");
1481 397166 : if (isintzero(dm)) return ZM_hnf(x);
1482 397164 : return ZM_hnfmodall(x, dm, flag);
1483 : }
1484 : GEN
1485 14 : hnfmod(GEN x, GEN d)
1486 : {
1487 14 : if (typ(d) != t_INT) pari_err_TYPE("mathnfmod",d);
1488 14 : return allhnfmod(x, d, 0);
1489 : }
1490 : GEN
1491 397158 : hnfmodid(GEN x, GEN d)
1492 : {
1493 397158 : switch(typ(d))
1494 : {
1495 380007 : case t_INT: break;
1496 17151 : case t_VEC: case t_COL:
1497 17151 : if (RgV_is_ZV(d)) break;
1498 0 : default: pari_err_TYPE("mathnfmodid",d);
1499 : }
1500 397158 : return allhnfmod(x, d, hnf_MODID);
1501 : }
1502 :
1503 : /* M a ZM in HNF. Normalize with *centered* residues */
1504 : GEN
1505 25100 : ZM_hnfcenter(GEN M)
1506 : {
1507 25100 : long i, j, k, N = lg(M)-1;
1508 25100 : pari_sp av = avma;
1509 :
1510 88413 : for (j=N-1; j>0; j--) /* skip last line */
1511 : {
1512 63308 : GEN Mj = gel(M,j), a = gel(Mj,j);
1513 211002 : for (k = j+1; k <= N; k++)
1514 : {
1515 147689 : GEN Mk = gel(M,k), q = diviiround(gel(Mk,j), a);
1516 147697 : long s = signe(q);
1517 147697 : if (!s) continue;
1518 69021 : if (is_pm1(q))
1519 : {
1520 30922 : if (s < 0)
1521 7417 : for (i = 1; i <= j; i++) gel(Mk,i) = addii(gel(Mk,i), gel(Mj,i));
1522 : else
1523 88672 : for (i = 1; i <= j; i++) gel(Mk,i) = subii(gel(Mk,i), gel(Mj,i));
1524 : }
1525 : else
1526 212159 : for (i = 1; i <= j; i++) gel(Mk,i) = subii(gel(Mk,i), mulii(q,gel(Mj,i)));
1527 69018 : if (gc_needed(av,1))
1528 : {
1529 0 : if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnfcenter, j = %ld",j);
1530 0 : M = gerepilecopy(av, M);
1531 : }
1532 : }
1533 : }
1534 25105 : return M;
1535 : }
1536 :
1537 : /***********************************************************************/
1538 : /* */
1539 : /* HNFLLL (Havas, Majewski, Mathews) */
1540 : /* */
1541 : /***********************************************************************/
1542 :
1543 : static void
1544 2026469 : Minus(long j, GEN lambda)
1545 : {
1546 2026469 : long k, n = lg(lambda);
1547 :
1548 11822518 : for (k=1 ; k<j; k++) togglesign_safe(&gcoeff(lambda,k,j));
1549 17209904 : for (k=j+1; k<n; k++) togglesign_safe(&gcoeff(lambda,j,k));
1550 2026434 : }
1551 :
1552 : /* index of first nonzero entry */
1553 : static long
1554 110659698 : findi(GEN M)
1555 : {
1556 110659698 : long i, n = lg(M);
1557 2901975200 : for (i=1; i<n; i++)
1558 2871075733 : if (signe(gel(M,i))) return i;
1559 30899467 : return 0;
1560 : }
1561 :
1562 : static long
1563 110650730 : findi_normalize(GEN Aj, GEN B, long j, GEN lambda)
1564 : {
1565 110650730 : long r = findi(Aj);
1566 110700491 : if (r && signe(gel(Aj,r)) < 0)
1567 : {
1568 2026447 : ZV_togglesign(Aj); if (B) ZV_togglesign(gel(B,j));
1569 2026472 : Minus(j,lambda);
1570 : }
1571 110703968 : return r;
1572 : }
1573 :
1574 : static void
1575 55284771 : reduce2(GEN A, GEN B, long k, long j, long *row0, long *row1, GEN lambda, GEN D)
1576 : {
1577 : GEN q;
1578 : long i;
1579 :
1580 55284771 : *row0 = findi_normalize(gel(A,j), B,j,lambda);
1581 55357800 : *row1 = findi_normalize(gel(A,k), B,k,lambda);
1582 55373655 : if (*row0)
1583 33563012 : q = truedivii(gcoeff(A,*row0,k), gcoeff(A,*row0,j));
1584 21810643 : else if (abscmpii(shifti(gcoeff(lambda,j,k), 1), gel(D,j)) > 0)
1585 7227548 : q = diviiround(gcoeff(lambda,j,k), gel(D,j));
1586 : else
1587 14570125 : return;
1588 :
1589 40773060 : if (signe(q))
1590 : {
1591 14285891 : GEN Lk = gel(lambda,k), Lj = gel(lambda,j);
1592 14285891 : togglesign_safe(&q);
1593 14304605 : if (*row0) ZC_lincomb1_inplace(gel(A,k),gel(A,j),q);
1594 14301087 : if (B) ZC_lincomb1_inplace(gel(B,k),gel(B,j),q);
1595 14242842 : gel(Lk,j) = addmulii(gel(Lk,j), q, gel(D,j));
1596 14228730 : if (is_pm1(q))
1597 : {
1598 5379143 : if (signe(q) > 0)
1599 : {
1600 5364461 : for (i=1; i<j; i++)
1601 3314601 : if (signe(gel(Lj,i))) gel(Lk,i) = addii(gel(Lk,i), gel(Lj,i));
1602 : }
1603 : else
1604 : {
1605 11453886 : for (i=1; i<j; i++)
1606 8125387 : if (signe(gel(Lj,i))) gel(Lk,i) = subii(gel(Lk,i), gel(Lj,i));
1607 : }
1608 : }
1609 : else
1610 : {
1611 39125726 : for (i=1; i<j; i++)
1612 30232065 : if (signe(gel(Lj,i))) gel(Lk,i) = addmulii(gel(Lk,i), q, gel(Lj,i));
1613 : }
1614 : }
1615 : }
1616 :
1617 : static void
1618 8385480 : hnfswap(GEN A, GEN B, long k, GEN lambda, GEN D)
1619 : {
1620 8385480 : GEN t, p1, p2, Lk = gel(lambda,k);
1621 8385480 : long i,j,n = lg(A);
1622 :
1623 8385480 : swap(gel(A,k), gel(A,k-1));
1624 8385480 : if (B) swap(gel(B,k), gel(B,k-1));
1625 48529458 : for (j=k-2; j; j--) swap(gcoeff(lambda,j,k-1), gel(Lk,j));
1626 108391754 : for (i=k+1; i<n; i++)
1627 : {
1628 100005809 : GEN Li = gel(lambda,i);
1629 100005809 : p1 = mulii(gel(Li,k-1), gel(D,k));
1630 99899510 : p2 = mulii(gel(Li,k), gel(Lk,k-1));
1631 99924997 : t = subii(p1,p2);
1632 :
1633 99933576 : p1 = mulii(gel(Li,k), gel(D,k-2));
1634 99911666 : p2 = mulii(gel(Li,k-1), gel(Lk,k-1));
1635 99915122 : gel(Li,k-1) = diviiexact(addii(p1,p2), gel(D,k-1));
1636 99959665 : gel(Li,k) = diviiexact(t, gel(D,k-1));
1637 : }
1638 8385945 : p1 = mulii(gel(D,k-2), gel(D,k));
1639 8360970 : p2 = sqri(gel(Lk,k-1));
1640 8362438 : gel(D,k-1) = diviiexact(addii(p1,p2), gel(D,k-1));
1641 8363260 : }
1642 :
1643 : /* reverse row order in matrix A, IN PLACE */
1644 : static GEN
1645 420577 : reverse_rows(GEN A)
1646 : {
1647 420577 : long i, j, h, n = lg(A);
1648 420577 : if (n == 1) return A;
1649 420577 : h = lgcols(A);
1650 3738095 : for (j=1; j<n; j++)
1651 : {
1652 3317517 : GEN c = gel(A,j);
1653 : /* start at (h-1) >>1 : if h = 2i even, no need to swap c[i] and itself */
1654 9446733 : for (i=(h-1)>>1; i; i--) swap(gel(c,i), gel(c,h-i));
1655 : }
1656 420578 : return A;
1657 : }
1658 : /* decide whether to swap */
1659 : static int
1660 4698402 : must_swap(long k, GEN lambda, GEN D)
1661 : {
1662 4698402 : pari_sp av = avma;
1663 4698402 : GEN z = addii(mulii(gel(D,k-2),gel(D,k)), sqri(gcoeff(lambda,k-1,k)));
1664 4692736 : return gc_bool(av, cmpii(z, sqri(gel(D,k-1))) < 0);
1665 : }
1666 :
1667 : GEN
1668 210292 : ZM_hnflll(GEN A, GEN *ptB, int remove)
1669 : {
1670 210292 : pari_sp av = avma;
1671 : long n, k, kmax;
1672 : GEN B, lambda, D;
1673 :
1674 210292 : n = lg(A);
1675 210292 : A = reverse_rows(ZM_copy(A)); /* ZM_copy for in place findi_normalize() */
1676 210291 : B = ptB? matid(n-1): NULL;
1677 210292 : D = const_vec(n, gen_1) + 1;
1678 210291 : lambda = zeromatcopy(n-1,n-1);
1679 185685 : k = kmax = 2;
1680 15557566 : while (k < n)
1681 : {
1682 : long row0, row1;
1683 : int do_swap;
1684 15347281 : reduce2(A,B,k,k-1,&row0,&row1,lambda,D);
1685 15353888 : if (row0) do_swap = (!row1 || row0 <= row1);
1686 5560496 : else if (row1) do_swap = 0;
1687 4429011 : else do_swap = must_swap(k,lambda,D);
1688 15394986 : if (do_swap)
1689 : {
1690 8135550 : hnfswap(A,B,k,lambda,D);
1691 8121610 : if (k > 2) k--;
1692 : }
1693 : else
1694 : {
1695 : long i;
1696 47182058 : for (i=k-2; i; i--)
1697 : {
1698 : long row0, row1;
1699 39931787 : reduce2(A,B,k,i,&row0,&row1,lambda,D);
1700 : }
1701 7250271 : if (++k > kmax) kmax = k;
1702 : }
1703 15371881 : if (gc_needed(av,3))
1704 : {
1705 182 : GEN b = D-1;
1706 182 : if (DEBUGMEM>1) pari_warn(warnmem,"hnflll, kmax = %ld / %ld",kmax,n-1);
1707 182 : gerepileall(av, B? 4: 3, &A, &lambda, &b, &B);
1708 182 : if (gc_needed(av,1)) paristack_resize(0); /* avoid desperation GC */
1709 182 : D = b+1;
1710 : }
1711 : }
1712 : /* handle trivial case: return negative diag coefficient otherwise */
1713 210285 : if (n == 2) (void)findi_normalize(gel(A,1), B,1,lambda);
1714 210285 : A = reverse_rows(A);
1715 210291 : if (remove)
1716 : {
1717 : long i;
1718 7343 : for (i = 1; i < n; i++)
1719 7119 : if (!ZV_equal0(gel(A,i))) break;
1720 2100 : remove_0cols(i-1, &A, &B, remove);
1721 : }
1722 210291 : gerepileall(av, B? 2: 1, &A, &B);
1723 210291 : if (B) *ptB = B;
1724 210291 : return A;
1725 : }
1726 :
1727 : GEN
1728 7 : hnflll(GEN x)
1729 : {
1730 7 : GEN z = cgetg(3, t_VEC);
1731 7 : gel(z,1) = ZM_hnflll(x, &gel(z,2), 1);
1732 7 : return z;
1733 : }
1734 :
1735 : /* Variation on HNFLLL: Extended GCD */
1736 :
1737 : static void
1738 952053 : reduce1(GEN A, GEN B, long k, long j, GEN lambda, GEN D)
1739 : {
1740 : GEN q;
1741 : long i;
1742 :
1743 952053 : if (signe(gel(A,j)))
1744 172947 : q = diviiround(gel(A,k),gel(A,j));
1745 779106 : else if (abscmpii(shifti(gcoeff(lambda,j,k), 1), gel(D,j)) > 0)
1746 96291 : q = diviiround(gcoeff(lambda,j,k), gel(D,j));
1747 : else
1748 682815 : return;
1749 :
1750 269238 : if (signe(q))
1751 : {
1752 232233 : GEN Lk = gel(lambda,k), Lj = gel(lambda,j);
1753 232233 : togglesign_safe(&q);
1754 232233 : gel(A,k) = addmulii(gel(A,k), q, gel(A,j));
1755 232233 : ZC_lincomb1_inplace(gel(B,k),gel(B,j),q);
1756 232233 : gel(Lk,j) = addmulii(gel(Lk,j), q, gel(D,j));
1757 422095 : for (i=1; i<j; i++)
1758 189862 : if (signe(gel(Lj,i))) gel(Lk,i) = addmulii(gel(Lk,i), q, gel(Lj,i));
1759 : }
1760 : }
1761 :
1762 : static GEN
1763 108829 : ZV_gcdext_i(GEN A)
1764 : {
1765 108829 : long k, n = lg(A);
1766 : GEN B, lambda, D;
1767 :
1768 108829 : if (n == 1) retmkvec2(gen_1, cgetg(1,t_MAT));
1769 108829 : A = leafcopy(A);
1770 108829 : B = matid(n-1);
1771 108829 : lambda = zeromatcopy(n-1,n-1);
1772 108829 : D = const_vec(n, gen_1) + 1;
1773 108829 : k = 2;
1774 698617 : while (k < n)
1775 : {
1776 : int do_swap;
1777 :
1778 589788 : reduce1(A,B,k,k-1,lambda,D);
1779 589788 : if (signe(gel(A,k-1))) do_swap = 1;
1780 416841 : else if (signe(gel(A,k))) do_swap = 0;
1781 224153 : else do_swap = must_swap(k,lambda,D);
1782 589788 : if (do_swap)
1783 : {
1784 241535 : hnfswap(A,B,k,lambda,D);
1785 241535 : if (k > 2) k--;
1786 : }
1787 : else
1788 : {
1789 : long i;
1790 710518 : for (i=k-2; i; i--) reduce1(A,B,k,i,lambda,D);
1791 348253 : k++;
1792 : }
1793 : }
1794 108829 : if (signe(gel(A,n-1)) < 0)
1795 : {
1796 13441 : gel(A,n-1) = negi(gel(A,n-1));
1797 13441 : ZV_togglesign(gel(B,n-1));
1798 : }
1799 108829 : return mkvec2(gel(A,n-1), B);
1800 : }
1801 : GEN
1802 108815 : ZV_extgcd(GEN A)
1803 : {
1804 108815 : pari_sp av = avma;
1805 108815 : return gerepilecopy(av, ZV_gcdext_i(A));
1806 : }
1807 : /* as ZV_extgcd, transforming the gcd into a t_MAT, for mathnf0 */
1808 : static GEN
1809 21 : ZV_hnfgcdext(GEN A)
1810 : {
1811 21 : pari_sp av = avma;
1812 : GEN z;
1813 21 : if (lg(A) == 1) retmkvec2(cgetg(1,t_MAT),cgetg(1,t_MAT));
1814 14 : z = ZV_gcdext_i(A);
1815 14 : gel(z,1) = mkmat(mkcol(gel(z,1)));
1816 14 : return gerepilecopy(av, z);
1817 : }
1818 :
1819 : /* HNF with permutation. */
1820 : GEN
1821 385 : ZM_hnfperm(GEN A, GEN *ptU, GEN *ptperm)
1822 : {
1823 : GEN U, c, l, perm, d, p, q, b;
1824 385 : pari_sp av = avma, av1;
1825 : long r, t, i, j, j1, k, m, n;
1826 :
1827 385 : n = lg(A)-1;
1828 385 : if (!n)
1829 : {
1830 7 : if (ptU) *ptU = cgetg(1,t_MAT);
1831 7 : if (ptperm) *ptperm = cgetg(1,t_VEC);
1832 7 : return cgetg(1, t_MAT);
1833 : }
1834 378 : m = nbrows(A);
1835 378 : c = zero_zv(m);
1836 378 : l = zero_zv(n);
1837 378 : perm = cgetg(m+1, t_VECSMALL);
1838 378 : av1 = avma;
1839 378 : A = RgM_shallowcopy(A);
1840 378 : U = ptU? matid(n): NULL;
1841 : /* U base change matrix : A0*U = A all along */
1842 1722 : for (r=0, k=1; k <= n; k++)
1843 : {
1844 3976 : for (j=1; j<k; j++)
1845 : {
1846 2632 : if (!l[j]) continue;
1847 2478 : t=l[j]; b=gcoeff(A,t,k);
1848 2478 : if (!signe(b)) continue;
1849 :
1850 546 : ZC_elem(b,gcoeff(A,t,j), A,U,k,j);
1851 546 : d = gcoeff(A,t,j);
1852 546 : if (signe(d) < 0)
1853 : {
1854 0 : ZV_neg_inplace(gel(A,j));
1855 0 : if (U) ZV_togglesign(gel(U,j));
1856 0 : d = gcoeff(A,t,j);
1857 : }
1858 1400 : for (j1=1; j1<j; j1++)
1859 : {
1860 854 : if (!l[j1]) continue;
1861 833 : q = truedivii(gcoeff(A,t,j1),d);
1862 833 : if (!signe(q)) continue;
1863 :
1864 329 : togglesign(q);
1865 329 : ZC_lincomb1_inplace(gel(A,j1), gel(A,j), q);
1866 329 : if (U) ZC_lincomb1_inplace(gel(U,j1), gel(U,j), q);
1867 : }
1868 : }
1869 4585 : t = m; while (t && (c[t] || !signe(gcoeff(A,t,k)))) t--;
1870 1344 : if (t)
1871 : {
1872 1225 : p = gcoeff(A,t,k);
1873 19971 : for (i=t-1; i; i--)
1874 : {
1875 18746 : q = gcoeff(A,i,k);
1876 18746 : if (signe(q) && abscmpii(p,q) > 0) { p = q; t = i; }
1877 : }
1878 1225 : perm[++r] = l[k] = t; c[t] = k;
1879 1225 : if (signe(p) < 0)
1880 : {
1881 145 : ZV_neg_inplace(gel(A,k));
1882 145 : if (U) ZV_togglesign(gel(U,k));
1883 145 : p = gcoeff(A,t,k);
1884 : }
1885 : /* p > 0 */
1886 3493 : for (j=1; j<k; j++)
1887 : {
1888 2268 : if (!l[j]) continue;
1889 2233 : q = truedivii(gcoeff(A,t,j),p);
1890 2233 : if (!signe(q)) continue;
1891 :
1892 238 : togglesign(q);
1893 238 : ZC_lincomb1_inplace(gel(A,j), gel(A,k), q);
1894 238 : if (U) ZC_lincomb1_inplace(gel(U,j), gel(U,k), q);
1895 : }
1896 : }
1897 1344 : if (gc_needed(av1,1))
1898 : {
1899 0 : if (DEBUGMEM>1) pari_warn(warnmem,"hnfperm, k=%ld",k);
1900 0 : gerepileall(av1, U? 2: 1, &A, &U);
1901 : }
1902 : }
1903 378 : if (r < m)
1904 : {
1905 5131 : for (i=1,k=r; i<=m; i++)
1906 4816 : if (!c[i]) perm[++k] = i;
1907 : }
1908 :
1909 : /* We have A0*U=A, U in Gl(n,Z)
1910 : * basis for Im(A): columns of A s.t l[j]>0 (r cols)
1911 : * basis for Ker(A): columns of U s.t l[j]=0 (n-r cols) */
1912 378 : p = cgetg(r+1,t_MAT);
1913 2814 : for (i=1; i<=m/2; i++) lswap(perm[i], perm[m+1-i]);
1914 378 : if (U)
1915 : {
1916 84 : GEN u = cgetg(n+1,t_MAT);
1917 378 : for (t=1,k=r,j=1; j<=n; j++)
1918 294 : if (l[j])
1919 : {
1920 182 : u[k + n-r] = U[j];
1921 182 : gel(p,k--) = vecpermute(gel(A,j), perm);
1922 : }
1923 : else
1924 112 : u[t++] = U[j];
1925 84 : *ptU = u;
1926 84 : if (ptperm) *ptperm = perm;
1927 84 : gerepileall(av, ptperm? 3: 2, &p, ptU, ptperm);
1928 : }
1929 : else
1930 : {
1931 1344 : for (k=r,j=1; j<=n; j++)
1932 1050 : if (l[j]) gel(p,k--) = vecpermute(gel(A,j), perm);
1933 294 : if (ptperm) *ptperm = perm;
1934 294 : gerepileall(av, ptperm? 2: 1, &p, ptperm);
1935 : }
1936 378 : return p;
1937 : }
1938 :
1939 : GEN
1940 238 : ZM_hnf_knapsack(GEN x)
1941 : {
1942 238 : GEN t, perm, H = ZM_hnfperm(x,NULL,&perm);
1943 238 : long i,j, l = lg(H), h = lgcols(H);
1944 3213 : for (i=1; i<h; i++)
1945 : {
1946 3066 : int fl = 0;
1947 16690 : for (j=1; j<l; j++)
1948 : {
1949 13715 : t = gcoeff(H,i,j);
1950 13715 : if (signe(t))
1951 : {
1952 3101 : if (!is_pm1(t) || fl) return NULL;
1953 3010 : fl = 1;
1954 : }
1955 : }
1956 : }
1957 147 : return rowpermute(H, perm_inv(perm));
1958 : }
1959 :
1960 : GEN
1961 14 : hnfperm(GEN A)
1962 : {
1963 14 : GEN y = cgetg(4, t_VEC);
1964 14 : gel(y,1) = ZM_hnfperm(A, &gel(y,2), &gel(y,3));
1965 14 : return y;
1966 : }
1967 :
1968 : /* Hermite Normal Form, with base change matrix if ptB != NULL.
1969 : * If 'remove' = 1, remove 0 columns (do NOT update *ptB accordingly)
1970 : * If 'remove' = 2, remove 0 columns and update *ptB accordingly */
1971 : GEN
1972 717370 : ZM_hnfall_i(GEN A, GEN *ptB, long remove)
1973 : {
1974 : pari_sp av;
1975 : long m, n, r, i, j, k, li;
1976 : GEN B, c, h, a;
1977 :
1978 717370 : RgM_dimensions(A, &m,&n);
1979 717366 : if (!n)
1980 : {
1981 7 : if (ptB) *ptB = cgetg(1,t_MAT);
1982 7 : return cgetg(1,t_MAT);
1983 : }
1984 717359 : c = zero_zv(m);
1985 717353 : h = const_vecsmall(n, m);
1986 717356 : av = avma;
1987 717356 : A = RgM_shallowcopy(A);
1988 717370 : B = ptB? matid(n): NULL;
1989 717370 : r = n+1;
1990 2015422 : for (li=m; li; li--)
1991 : {
1992 2897860 : for (j=1; j<r; j++)
1993 : {
1994 3529234 : for (i=h[j]; i>li; i--)
1995 : {
1996 631926 : a = gcoeff(A,i,j);
1997 631926 : k = c[i];
1998 : /* zero a = Aij using Aik */
1999 631926 : if (signe(a)) ZC_elem(a,gcoeff(A,i,k), A,B,j,k);
2000 631878 : ZM_reduce(A,B, i,k); /* ensure reduced entries even if a = 0 */
2001 : }
2002 2897308 : if (gc_needed(av,1) && (j & 0x7f) == 0)
2003 : {
2004 0 : if (DEBUGMEM>1)
2005 0 : pari_warn(warnmem,"ZM_hnfall[1], li = %ld, j = %ld", li, j);
2006 0 : gerepileall(av, B? 2: 1, &A, &B);
2007 : }
2008 2897345 : if (signe( gcoeff(A,li,j) )) break;
2009 1599795 : h[j] = li-1;
2010 : }
2011 1298082 : if (j == r) continue;
2012 1297595 : r--;
2013 1297595 : if (j < r) /* A[j] != 0 */
2014 : {
2015 839547 : swap(gel(A,j), gel(A,r));
2016 839547 : if (B) swap(gel(B,j), gel(B,r));
2017 839547 : h[j] = h[r]; h[r] = li; c[li] = r;
2018 : }
2019 1297595 : if (signe(gcoeff(A,li,r)) < 0)
2020 : {
2021 300732 : ZV_neg_inplace(gel(A,r));
2022 300740 : if (B) ZV_togglesign(gel(B,r));
2023 : }
2024 1297601 : ZM_reduce(A,B, li,r);
2025 1297564 : if (gc_needed(av,1))
2026 : {
2027 0 : if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnfall[2], li = %ld", li);
2028 0 : gerepileall(av, B? 2: 1, &A, &B);
2029 : }
2030 : }
2031 :
2032 717357 : if (DEBUGLEVEL>5) err_printf("\nhnfall, final phase: ");
2033 717362 : r--; /* first r cols are in the image the n-r (independent) last ones */
2034 2462321 : for (j=1; j<=r; j++)
2035 : {
2036 3877354 : for (i=h[j]; i; i--)
2037 : {
2038 2132407 : a = gcoeff(A,i,j);
2039 2132407 : k = c[i];
2040 2132407 : if (signe(a)) ZC_elem(a,gcoeff(A,i,k), A,B, j,k);
2041 2132349 : ZM_reduce(A,B, i,k); /* ensure reduced entries, even if a = 0 */
2042 : }
2043 1744947 : if (gc_needed(av,1) && (j & 0x7f) == 0)
2044 : {
2045 0 : if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnfall[3], j = %ld", j);
2046 0 : gerepileall(av, B? 2: 1, &A, &B);
2047 : }
2048 : }
2049 717347 : if (DEBUGLEVEL>5) err_printf("\n");
2050 717347 : if (remove) remove_0cols(r, &A, &B, remove);
2051 717356 : if (ptB) *ptB = B;
2052 717356 : return A;
2053 : }
2054 : GEN
2055 23676 : ZM_hnfall(GEN A, GEN *ptB, long remove)
2056 : {
2057 23676 : pari_sp av = avma;
2058 23676 : A = ZM_hnfall_i(A, ptB, remove);
2059 23676 : return gc_all(av, ptB? 2: 1, &A, ptB);
2060 : }
2061 :
2062 : GEN
2063 14 : hnfall(GEN x)
2064 : {
2065 14 : GEN z = cgetg(3, t_VEC);
2066 14 : gel(z,1) = ZM_hnfall(x, (GEN*)(z+2), 1);
2067 14 : return z;
2068 : }
2069 : GEN
2070 14 : hnf(GEN x) { return mathnf0(x,0); }
2071 :
2072 : /* C = A^(-1)t where A and C are integral, A is upper triangular, t t_INT */
2073 : GEN
2074 236584 : hnf_invscale(GEN A, GEN t)
2075 : {
2076 236584 : long n = lg(A)-1, i,j,k;
2077 236584 : GEN m, c = cgetg(n+1,t_MAT);
2078 :
2079 236585 : if (!n) return c;
2080 877445 : for (k=1; k<=n; k++)
2081 : { /* cf hnf_divscale with B = id, thus b = e_k */
2082 640875 : GEN u = cgetg(n+1, t_COL);
2083 640880 : pari_sp av = avma;
2084 640880 : gel(c,k) = u;
2085 640880 : gel(u,n) = k == n? gerepileuptoint(av, diviiexact(t, gcoeff(A,n,n))): gen_0;
2086 2167769 : for (i=n-1; i>0; i--)
2087 : {
2088 1526909 : av = avma; m = i == k? t: gen_0;
2089 5777809 : for (j=i+1; j<=n; j++) m = subii(m, mulii(gcoeff(A,i,j),gel(u,j)));
2090 1526645 : gel(u,i) = gerepileuptoint(av, diviiexact(m, gcoeff(A,i,i)));
2091 : }
2092 : }
2093 236570 : return c;
2094 : }
2095 :
2096 : /* C = A^(-1)(tB) where A, B, C are integral, A is upper triangular, t t_INT */
2097 : GEN
2098 196857 : hnf_divscale(GEN A, GEN B, GEN t)
2099 : {
2100 196857 : long n = lg(A)-1, i,j,k;
2101 196857 : GEN m, c = cgetg(n+1,t_MAT);
2102 :
2103 196857 : if (!n) return c;
2104 887408 : for (k=1; k<=n; k++)
2105 : {
2106 690551 : GEN u = cgetg(n+1, t_COL), b = gel(B,k);
2107 690552 : pari_sp av = avma;
2108 690552 : gel(c,k) = u; m = mulii(gel(b,n),t);
2109 690536 : gel(u,n) = gerepileuptoint(av, diviiexact(m, gcoeff(A,n,n)));
2110 4103775 : for (i=n-1; i>0; i--)
2111 : {
2112 3413224 : av = avma; m = mulii(gel(b,i),t);
2113 31738089 : for (j=i+1; j<=n; j++) m = subii(m, mulii(gcoeff(A,i,j),gel(u,j)));
2114 3413083 : gel(u,i) = gerepileuptoint(av, diviiexact(m, gcoeff(A,i,i)));
2115 : }
2116 : }
2117 196857 : return c;
2118 : }
2119 :
2120 : /* A, B integral upper HNF. A^(-1) B integral ? */
2121 : int
2122 133 : hnfdivide(GEN A, GEN B)
2123 : {
2124 133 : pari_sp av = avma;
2125 133 : long n = lg(A)-1, i,j,k;
2126 : GEN u, b, m, r;
2127 :
2128 133 : if (!n) return 1;
2129 133 : if (lg(B)-1 != n) pari_err_DIM("hnfdivide");
2130 133 : u = cgetg(n+1, t_COL);
2131 483 : for (k=1; k<=n; k++)
2132 : {
2133 364 : b = gel(B,k);
2134 364 : m = gel(b,k);
2135 364 : gel(u,k) = dvmdii(m, gcoeff(A,k,k), &r);
2136 364 : if (r != gen_0) return gc_long(av, 0);
2137 826 : for (i=k-1; i>0; i--)
2138 : {
2139 476 : m = gel(b,i);
2140 1337 : for (j=i+1; j<=k; j++) m = subii(m, mulii(gcoeff(A,i,j),gel(u,j)));
2141 476 : m = dvmdii(m, gcoeff(A,i,i), &r);
2142 476 : if (r != gen_0) return gc_long(av, 0);
2143 476 : gel(u,i) = m;
2144 : }
2145 : }
2146 119 : return gc_long(av, 1);
2147 : }
2148 :
2149 : /* A upper HNF, b integral vector. Return A^(-1) b if integral,
2150 : * NULL otherwise. Assume #A[,1] = #b. */
2151 : GEN
2152 365208 : hnf_invimage(GEN A, GEN b)
2153 : {
2154 365208 : pari_sp av = avma;
2155 365208 : long n = lg(A)-1, m, i, k;
2156 : GEN u, r;
2157 :
2158 365208 : if (!n) return lg(b)==1? cgetg(1,t_COL):NULL;
2159 365166 : m = nbrows(A); /* m >= n */
2160 365167 : u = cgetg(n+1, t_COL);
2161 826023 : for (i = n, k = m; k > 0; k--)
2162 : {
2163 826023 : pari_sp av2 = avma;
2164 : long j;
2165 826023 : GEN t = gel(b,k), Aki = gcoeff(A,k,i);
2166 826023 : if (typ(t) != t_INT) pari_err_TYPE("hnf_invimage",t);
2167 1981952 : for (j=i+1; j<=n; j++) t = subii(t, mulii(gcoeff(A,k,j),gel(u,j)));
2168 825933 : if (!signe(Aki))
2169 : {
2170 0 : if (signe(t)) return gc_NULL(av);
2171 0 : set_avma(av2); gel(u,i) = gen_0; continue;
2172 : }
2173 825933 : t = dvmdii(t, Aki, &r);
2174 825870 : if (r != gen_0) return gc_NULL(av);
2175 603609 : gel(u,i) = gerepileuptoint(av2, t);
2176 603747 : if (--i == 0) break;
2177 : }
2178 : /* If there is a solution, it must be u. Check remaining equations */
2179 285755 : for (; k > 0; k--)
2180 : {
2181 142880 : pari_sp av2 = avma;
2182 : long j;
2183 142880 : GEN t = gel(b,k);
2184 142880 : if (typ(t) != t_INT) pari_err_TYPE("hnf_invimage",t);
2185 544993 : for (j=1; j<=n; j++) t = subii(t, mulii(gcoeff(A,k,j),gel(u,j)));
2186 142876 : if (signe(t)) return gc_NULL(av);
2187 142876 : set_avma(av2);
2188 : }
2189 142875 : return u;
2190 : }
2191 :
2192 : /* A upper HNF, B integral matrix or column. Return A^(-1) B if integral,
2193 : * NULL otherwise */
2194 : GEN
2195 313891 : hnf_solve(GEN A, GEN B)
2196 : {
2197 : pari_sp av;
2198 : long i, l;
2199 : GEN C;
2200 :
2201 313891 : if (typ(B) == t_COL) return hnf_invimage(A, B);
2202 251718 : av = avma; C = cgetg_copy(B, &l);
2203 350651 : for (i = 1; i < l; i++)
2204 : {
2205 265926 : GEN c = hnf_invimage(A, gel(B,i));
2206 265891 : if (!c) return gc_NULL(av);
2207 98914 : gel(C,i) = c;
2208 : }
2209 84725 : return C;
2210 : }
2211 :
2212 : /***************************************************************/
2213 : /** **/
2214 : /** SMITH NORMAL FORM REDUCTION **/
2215 : /** **/
2216 : /***************************************************************/
2217 :
2218 : static GEN
2219 0 : trivsmith(long all)
2220 : {
2221 : GEN z;
2222 0 : if (!all) return cgetg(1,t_VEC);
2223 0 : z=cgetg(4,t_VEC);
2224 0 : gel(z,1) = cgetg(1,t_MAT);
2225 0 : gel(z,2) = cgetg(1,t_MAT);
2226 0 : gel(z,3) = cgetg(1,t_MAT); return z;
2227 : }
2228 :
2229 : static void
2230 63 : snf_pile1(pari_sp av, GEN *x, GEN *U)
2231 : {
2232 : GEN *gptr[2];
2233 63 : int c = 1; gptr[0]=x;
2234 63 : if (*U) gptr[c++] = U;
2235 63 : gerepilemany(av,gptr,c);
2236 63 : }
2237 : static void
2238 971321 : snf_pile(pari_sp av, GEN *x, GEN *U, GEN *V)
2239 : {
2240 : GEN *gptr[3];
2241 971321 : int c = 1; gptr[0]=x;
2242 971321 : if (*U) gptr[c++] = U;
2243 971321 : if (*V) gptr[c++] = V;
2244 971321 : gerepilemany(av,gptr,c);
2245 971376 : }
2246 :
2247 : static GEN
2248 695722 : bezout_step(GEN *pa, GEN *pb, GEN *pu, GEN *pv)
2249 : {
2250 695722 : GEN a = *pa, b = *pb, d;
2251 695722 : if (absequalii(a,b))
2252 : {
2253 411813 : long sa = signe(a), sb = signe(b);
2254 411813 : *pv = gen_0;
2255 411813 : if (sb == sa) {
2256 403551 : *pa = *pb = gen_1;
2257 403551 : if (sa > 0) {*pu=gen_1; return a;} else {*pu=gen_m1; return absi(a);}
2258 : }
2259 8262 : if (sa > 0) { *pa = *pu = gen_1; *pb = gen_m1; return a; }
2260 1101 : *pa = *pu = gen_m1; *pb = gen_1; return b;
2261 : }
2262 283911 : d = bezout(a,b, pu,pv);
2263 283925 : *pa = diviiexact(a, d);
2264 283918 : *pb = diviiexact(b, d); return d;
2265 : }
2266 :
2267 : static int
2268 517250 : negcmpii(void *E, GEN x, GEN y) { (void)E; return -cmpii(x,y); }
2269 :
2270 : /* x square of maximal rank; does b = x[i,i] divide all entries in
2271 : * x[1..i-1, 1..i-1] ? If so, return 0; else the index of a problematic row */
2272 : static long
2273 1062928 : ZM_snf_no_divide(GEN x, long i)
2274 : {
2275 1062928 : GEN b = gcoeff(x,i,i);
2276 : long j, k;
2277 :
2278 1062928 : if (is_pm1(b)) return 0;
2279 907508 : for (k = 1; k < i; k++)
2280 2259462 : for (j = 1; j < i; j++)
2281 1784878 : if (!dvdii(gcoeff(x,k,j),b)) return k;
2282 283064 : return 0;
2283 : }
2284 :
2285 : static void
2286 1377519 : ZM_redpart(GEN x, GEN p, long I)
2287 : {
2288 1377519 : long l = lgefint(p), i, j;
2289 5877063 : for (i = 1; i <= I; i++)
2290 27108969 : for (j = 1; j <= I; j++)
2291 : {
2292 22609425 : GEN c = gcoeff(x,i,j);
2293 22609425 : if (lgefint(c) > l) gcoeff(x,i,j) = remii(c, p);
2294 : }
2295 1377519 : }
2296 : static void
2297 860881 : ZMrow_divexact_inplace(GEN M, long i, GEN c)
2298 : {
2299 860881 : long j, l = lg(M);
2300 3410142 : for (j = 1; j < l; j++) gcoeff(M,i,j) = diviiexact(gcoeff(M,i,j), c);
2301 860855 : }
2302 :
2303 : /* Return the SNF D of matrix X. If ptU/ptV non-NULL set them to U/V
2304 : * to that D = UXV */
2305 : GEN
2306 711169 : ZM_snfall_i(GEN x, GEN *ptU, GEN *ptV, long flag)
2307 : {
2308 711169 : pari_sp av0 = avma, av;
2309 711169 : const long return_vec = flag & 1;
2310 : long i, j, k, m0, m, n0, n;
2311 711169 : GEN u, v, U, V, V0, mdet, A = NULL, perm = NULL;
2312 :
2313 711169 : n0 = n = lg(x)-1;
2314 711169 : if (!n) {
2315 41007 : if (ptU) *ptU = cgetg(1,t_MAT);
2316 41007 : if (ptV) *ptV = cgetg(1,t_MAT);
2317 41007 : return cgetg(1, return_vec? t_VEC: t_MAT);
2318 : }
2319 670162 : m0 = m = nbrows(x);
2320 :
2321 670165 : U = V = V0 = NULL; /* U = TRANSPOSE of row transform matrix [act on columns]*/
2322 670165 : if (m == n && ZM_ishnf(x))
2323 : {
2324 610094 : mdet = ZM_det_triangular(x); /* != 0 */
2325 610049 : if (ptV) *ptV = matid(n);
2326 : }
2327 : else
2328 : {
2329 60091 : mdet = ZM_detmult(x);
2330 60099 : if (!signe(mdet))
2331 77 : x = ZM_hnfperm(x, ptV, ptU? &perm: NULL);
2332 : else
2333 : { /* m <= n */
2334 60022 : if (!ptV)
2335 357 : x = ZM_hnfmod(x,mdet);
2336 59665 : else if (m == n)
2337 : {
2338 59637 : GEN H = ZM_hnfmod(x,mdet);
2339 59638 : *ptV = ZM_gauss(x,H);
2340 59638 : x = H;
2341 : }
2342 : else
2343 28 : x = ZM_hnfperm(x, ptV, ptU? &perm: NULL);
2344 60023 : mdet = ZM_det_triangular(x); /* > 0 */
2345 : }
2346 60100 : n = lg(x)-1; /* n independent columns */
2347 60100 : if (ptV)
2348 : {
2349 59687 : V = *ptV;
2350 59687 : if (n != n0)
2351 : {
2352 28 : V0 = vecslice(V, 1, n0 - n); /* kernel */
2353 28 : V = vecslice(V, n0-n+1, n0);
2354 : }
2355 : }
2356 60100 : if (!signe(mdet))
2357 : {
2358 77 : if (n)
2359 : {
2360 70 : x = ZM_snfall_i(shallowtrans(x), ptV, ptU, return_vec); /* swap V,U */
2361 70 : if (!return_vec && n != m) x = shallowtrans(x);
2362 70 : if (ptV) V = ZM_mul(V, shallowtrans(*ptV));
2363 70 : if (ptU) U = *ptU; /* TRANSPOSE */
2364 : }
2365 : else /* 0 matrix */
2366 : {
2367 7 : x = cgetg(1,t_MAT);
2368 7 : if (ptV) V = cgetg(1, t_MAT);
2369 7 : if (ptU) U = matid(m);
2370 : }
2371 77 : goto THEEND;
2372 : }
2373 : }
2374 670128 : if (ptV || ptU) U = matid(n); /* we will compute V in terms of U */
2375 670129 : if (DEBUGLEVEL>7) err_printf("starting SNF loop");
2376 :
2377 : /* square, maximal rank n */
2378 670129 : A = x; x = shallowcopy(x); av = avma;
2379 1583207 : for (i = n; i > 1; i--)
2380 : {
2381 913091 : if (DEBUGLEVEL>7) err_printf("\ni = %ld: ",i);
2382 : for(;;)
2383 693659 : {
2384 1606729 : int c = 0;
2385 : GEN a, b;
2386 4622094 : for (j = i-1; j >= 1; j--)
2387 : {
2388 3015559 : b = gcoeff(x,i,j); if (!signe(b)) continue;
2389 263451 : a = gcoeff(x,i,i);
2390 263451 : ZC_elem(b, a, x,NULL, j,i);
2391 263440 : if (gc_needed(av,1))
2392 : {
2393 0 : if (DEBUGMEM>1) pari_warn(warnmem,"[1]: ZM_snfall i = %ld", i);
2394 0 : snf_pile1(av, &x,&U);
2395 : }
2396 : }
2397 1606535 : if (DEBUGLEVEL>7) err_printf("; ");
2398 4621992 : for (j=i-1; j>=1; j--)
2399 : {
2400 : GEN d;
2401 3015448 : b = gcoeff(x,j,i); if (!signe(b)) continue;
2402 695706 : a = gcoeff(x,i,i);
2403 695706 : d = bezout_step(&a, &b, &u, &v);
2404 4170557 : for (k = 1; k < i; k++)
2405 : {
2406 3474937 : GEN t = addii(mulii(u,gcoeff(x,i,k)),mulii(v,gcoeff(x,j,k)));
2407 3475020 : gcoeff(x,j,k) = subii(mulii(a,gcoeff(x,j,k)),
2408 3475027 : mulii(b,gcoeff(x,i,k)));
2409 3474855 : gcoeff(x,i,k) = t;
2410 : }
2411 695620 : gcoeff(x,j,i) = gen_0;
2412 695620 : gcoeff(x,i,i) = d;
2413 695620 : if (U) update(u,v,a,b,(GEN*)(U+i),(GEN*)(U+j));
2414 695624 : if (gc_needed(av,1))
2415 : {
2416 63 : if (DEBUGMEM>1) pari_warn(warnmem,"[2]: ZM_snfall, i = %ld", i);
2417 63 : snf_pile1(av, &x,&U);
2418 : }
2419 695626 : c = 1;
2420 : }
2421 1606544 : if (!c)
2422 : {
2423 1062930 : k = ZM_snf_no_divide(x, i);
2424 1062877 : if (!k) break;
2425 :
2426 : /* x[k,j] != 0 mod b */
2427 574453 : for (j = 1; j <= i; j++)
2428 424666 : gcoeff(x,i,j) = addii(gcoeff(x,i,j),gcoeff(x,k,j));
2429 149787 : if (U) gel(U,i) = gadd(gel(U,i),gel(U,k));
2430 : }
2431 693403 : ZM_redpart(x, mdet, i);
2432 693616 : if (U && (flag & 2)) ZM_redpart(U, mdet, n);
2433 693658 : if (gc_needed(av,1))
2434 : {
2435 0 : if (DEBUGMEM>1) pari_warn(warnmem,"[3]: ZM_snfall");
2436 0 : snf_pile1(av, &x,&U);
2437 : }
2438 : }
2439 : }
2440 670116 : if (DEBUGLEVEL>7) err_printf("\n");
2441 2252410 : for (k = n; k; k--)
2442 : {
2443 1582652 : GEN d = gcdii(gcoeff(x,k,k), mdet);
2444 1582382 : gcoeff(x,k,k) = d;
2445 1582382 : if (!is_pm1(d)) mdet = diviiexact(mdet,d);
2446 : }
2447 669758 : THEEND:
2448 669835 : if (U) U = shallowtrans(U);
2449 670015 : if (ptV && A)
2450 : { /* U A V = D => D^(-1) U A = V^(-1) */
2451 665157 : long l = lg(x);
2452 665157 : GEN W = ZM_mul(U, A);
2453 1525824 : for (i = 1; i < l; i++)
2454 : {
2455 1334108 : GEN c = gcoeff(x,i,i);
2456 1334108 : if (is_pm1(c)) break; /* only 1 from now on */
2457 860644 : ZMrow_divexact_inplace(W, i, c);
2458 : }
2459 665241 : if (flag & 2)
2460 : {
2461 641027 : W = FpM_red(W, gcoeff(x,1,1));
2462 641021 : W = matinvmod(W, gcoeff(x,1,1));
2463 : }
2464 : else
2465 24214 : W = ZM_inv(W, NULL);
2466 665190 : V = V? ZM_mul(V, W): W;
2467 : }
2468 670031 : if (return_vec)
2469 : {
2470 646022 : long l = lg(x)-1;
2471 646022 : if (typ(x) == t_MAT) x = RgM_diagonal_shallow(x);
2472 646007 : if (m0 > l) x = shallowconcat(zerovec(m0-l), x);
2473 : }
2474 :
2475 670016 : if (V0)
2476 : { /* add kernel */
2477 28 : if (!return_vec) x = shallowconcat(zeromat(m,n0-n), x);
2478 28 : if (ptV) V = shallowconcat(V0, V);
2479 : }
2480 670016 : if (perm && U) U = vecpermute(U, perm_inv(perm));
2481 670016 : snf_pile(av0, &x,&U,&V);
2482 670227 : if (ptU) *ptU = U;
2483 670227 : if (ptV) *ptV = V;
2484 670227 : return x;
2485 : }
2486 : GEN
2487 63784 : ZM_snfall(GEN x, GEN *U, GEN *V) { return ZM_snfall_i(x, U, V, 0); }
2488 : GEN
2489 1638 : ZM_snf(GEN x) { return ZM_snfall_i(x, NULL,NULL, 1); }
2490 : GEN
2491 385 : smith(GEN x) { return ZM_snfall_i(x, NULL,NULL, 1); }
2492 : GEN
2493 35 : smithall(GEN x)
2494 : {
2495 35 : GEN z = cgetg(4, t_VEC);
2496 35 : gel(z,3) = ZM_snfall_i(x, (GEN*)(z+1),(GEN*)(z+2), 0);
2497 35 : return z;
2498 : }
2499 :
2500 : GEN
2501 213962 : ZV_snfclean(GEN d)
2502 : {
2503 213962 : long c, l = lg(d);
2504 226328 : for (c = 1; c < l; c++) { GEN t = gel(d,c); if (is_pm1(t)) break; }
2505 213963 : return c == l? d: vec_shorten(d, c-1);
2506 : }
2507 : void
2508 943578 : ZM_snfclean(GEN d, GEN u, GEN v)
2509 : {
2510 943578 : long i, c, l = lg(d);
2511 :
2512 943578 : if (typ(d) == t_VEC)
2513 2423003 : for (c=1; c<l; c++) { GEN t = gel(d,c); if (is_pm1(t)) break; }
2514 : else
2515 : {
2516 0 : for (c=1; c<l; c++) { GEN t = gcoeff(d,c,c); if (is_pm1(t)) break; }
2517 0 : if (c < l) for (i = 1; i < c; i++) setlg(gel(d,i), c);
2518 : }
2519 943574 : setlg(d, c);
2520 2917443 : if (u) for (i=1; i<l; i++) setlg(gel(u,i), c);
2521 943545 : if (v) setlg(v, c);
2522 943543 : }
2523 :
2524 : /* Assume z was computed by [g]smithall(). Remove the 1s on the diagonal */
2525 : GEN
2526 693 : smithclean(GEN z)
2527 : {
2528 : long i, j, h, l, c, d;
2529 : GEN U, V, y, D, t;
2530 :
2531 693 : if (typ(z) != t_VEC) pari_err_TYPE("smithclean",z);
2532 693 : l = lg(z); if (l == 1) return cgetg(1,t_VEC);
2533 686 : U = gel(z,1);
2534 686 : if (l != 4 || typ(U) != t_MAT)
2535 : { /* assume z = vector of elementary divisors */
2536 1827 : for (c=1; c<l; c++)
2537 1519 : if (gequal1(gel(z,c))) break;
2538 665 : return gcopy_lg(z, c);
2539 : }
2540 21 : V = gel(z,2);
2541 21 : D = gel(z,3);
2542 21 : l = lg(D);
2543 21 : if (l == 1) return gcopy(z);
2544 21 : h = lgcols(D);
2545 21 : if (h > l)
2546 : { /* D = vconcat(zero matrix, diagonal matrix) */
2547 21 : for (c=1+h-l, d=1; c<h; c++,d++)
2548 21 : if (gequal1(gcoeff(D,c,d))) break;
2549 : }
2550 7 : else if (h < l)
2551 : { /* D = concat(zero matrix, diagonal matrix) */
2552 7 : for (c=1, d=1+l-h; d<l; c++,d++)
2553 7 : if (gequal1(gcoeff(D,c,d))) break;
2554 : }
2555 : else
2556 : { /* D diagonal */
2557 0 : for (c=1; c<l; c++)
2558 0 : if (gequal1(gcoeff(D,c,c))) break;
2559 0 : d = c;
2560 : }
2561 : /* U was (h-1)x(h-1), V was (l-1)x(l-1), D was (h-1)x(l-1) */
2562 21 : y = cgetg(4,t_VEC);
2563 : /* truncate U to (c-1) x (h-1) */
2564 21 : gel(y,1) = t = cgetg(h,t_MAT);
2565 77 : for (j=1; j<h; j++) gel(t,j) = gcopy_lg(gel(U,j), c);
2566 : /* truncate V to (l-1) x (d-1) */
2567 21 : gel(y,2) = gcopy_lg(V, d);
2568 21 : gel(y,3) = t = zeromatcopy(c-1, d-1);
2569 : /* truncate D to a (c-1) x (d-1) matrix */
2570 21 : if (d > 1)
2571 : {
2572 14 : if (h > l)
2573 : {
2574 14 : for (i=1+h-l, j=1; i<c; i++,j++)
2575 7 : gcoeff(t,i,j) = gcopy(gcoeff(D,i,j));
2576 : }
2577 7 : else if (h < l)
2578 : {
2579 7 : for (i=1, j=1+l-h; j<d; i++,j++)
2580 0 : gcoeff(t,i,j) = gcopy(gcoeff(D,i,j));
2581 : }
2582 : else
2583 : {
2584 0 : for (j=1; j<d; j++)
2585 0 : gcoeff(t,j,j) = gcopy(gcoeff(D,j,j));
2586 : }
2587 : }
2588 21 : return y;
2589 : }
2590 :
2591 : /* does b = x[i,i] divide all entries in x[1..i-1,1..i-1] ? If so, return 0;
2592 : * else return the index of a problematic row */
2593 : static long
2594 196 : gsnf_no_divide(GEN x, long i, long vx)
2595 : {
2596 196 : GEN b = gcoeff(x,i,i);
2597 : long j, k;
2598 :
2599 196 : if (gequal0(b))
2600 : {
2601 14 : for (k = 1; k < i; k++)
2602 14 : for (j = 1; j < i; j++)
2603 14 : if (!gequal0(gcoeff(x,k,j))) return k;
2604 0 : return 0;
2605 : }
2606 :
2607 182 : if (!is_RgX(b,vx) || degpol(b)<=0) return 0;
2608 217 : for (k = 1; k < i; k++)
2609 378 : for (j = 1; j < i; j++)
2610 : {
2611 266 : GEN z = gcoeff(x,k,j), r;
2612 266 : if (!is_RgX(z,vx)) z = scalarpol(z, vx);
2613 266 : r = RgX_rem(z, b);
2614 266 : if (signe(r) && (! isinexactreal(r) ||
2615 0 : gexpo(r) > 16 + gexpo(b) - prec2nbits(gprecision(r)))
2616 35 : ) return k;
2617 : }
2618 70 : return 0;
2619 : }
2620 :
2621 : /* Hermite Normal Form, with base change matrix if ptB != NULL.
2622 : * If 'remove' = 1, remove 0 columns (do NOT update *ptB accordingly)
2623 : * If 'remove' = 2, remove 0 columns and update *ptB accordingly */
2624 : GEN
2625 42 : RgM_hnfall(GEN A, GEN *pB, long remove)
2626 : {
2627 : pari_sp av;
2628 : long li, j, k, m, n, def, ldef;
2629 : GEN B;
2630 42 : long vx = gvar(A);
2631 :
2632 42 : n = lg(A)-1;
2633 42 : if (vx==NO_VARIABLE || !n)
2634 : {
2635 0 : RgM_check_ZM(A, "mathnf0");
2636 0 : return ZM_hnfall(A, pB, remove);
2637 : }
2638 42 : m = nbrows(A);
2639 42 : av = avma;
2640 42 : A = RgM_shallowcopy(A);
2641 42 : B = pB? matid(n): NULL;
2642 42 : def = n; ldef = (m>n)? m-n: 0;
2643 126 : for (li=m; li>ldef; li--)
2644 : {
2645 : GEN d, T;
2646 714 : for (j=def-1; j; j--)
2647 : {
2648 630 : GEN a = gcoeff(A,li,j);
2649 630 : if (gequal0(a)) continue;
2650 :
2651 511 : k = (j==1)? def: j-1;
2652 511 : RgC_elem(a,gcoeff(A,li,k), A,B, j,k, li, vx);
2653 : }
2654 84 : T = normalize_as_RgX(gcoeff(A,li,def), vx, &d);
2655 84 : if (gequal0(T))
2656 0 : { if (ldef) ldef--; }
2657 : else
2658 : {
2659 84 : if (!gequal1(d))
2660 : {
2661 49 : gel(A,def) = RgC_Rg_div(gel(A,def), d);
2662 49 : if (B) gel(B, def) = RgC_Rg_div(gel(B, def), d);
2663 : }
2664 84 : RgM_reduce(A, B, li, def, vx);
2665 84 : def--;
2666 : }
2667 84 : if (gc_needed(av,1))
2668 : {
2669 0 : if (DEBUGMEM>1) pari_warn(warnmem,"ghnfall");
2670 0 : gerepileall(av, B? 2: 1, &A, &B);
2671 : }
2672 : }
2673 : /* rank A = n - def */
2674 42 : if (remove) remove_0cols(def, &A, &B, remove);
2675 42 : gerepileall(av, B? 2: 1, &A, &B);
2676 42 : if (B) *pB = B;
2677 42 : return A;
2678 : }
2679 :
2680 : static GEN
2681 49 : RgXM_snf(GEN x,long all)
2682 : {
2683 : pari_sp av;
2684 : long i, j, k, n;
2685 : GEN z, u, v, U, V;
2686 49 : long vx = gvar(x);
2687 49 : n = lg(x)-1; if (!n) return trivsmith(all);
2688 49 : if (vx==NO_VARIABLE) pari_err_TYPE("RgXM_snf",x);
2689 49 : if (lgcols(x) != n+1) pari_err_DIM("gsmithall");
2690 49 : av = avma;
2691 49 : x = RgM_shallowcopy(x);
2692 49 : if (all) { U = matid(n); V = matid(n); }
2693 196 : for (i=n; i>=2; i--)
2694 : {
2695 : for(;;)
2696 168 : {
2697 : GEN a, b, d;
2698 315 : int c = 0;
2699 1092 : for (j=i-1; j>=1; j--)
2700 : {
2701 777 : b = gcoeff(x,i,j); if (gequal0(b)) continue;
2702 196 : a = gcoeff(x,i,i);
2703 196 : d = gbezout_step(&b, &a, &v, &u, vx);
2704 700 : for (k = 1; k < i; k++)
2705 : {
2706 504 : GEN t = gadd(gmul(u,gcoeff(x,k,i)),gmul(v,gcoeff(x,k,j)));
2707 504 : gcoeff(x,k,j) = gsub(gmul(a,gcoeff(x,k,j)),gmul(b,gcoeff(x,k,i)));
2708 504 : gcoeff(x,k,i) = t;
2709 : }
2710 196 : gcoeff(x,i,j) = gen_0;
2711 196 : gcoeff(x,i,i) = d;
2712 196 : if (all) update(u,v,a,b,(GEN*)(V+i),(GEN*)(V+j));
2713 : }
2714 1092 : for (j=i-1; j>=1; j--)
2715 : {
2716 777 : b = gcoeff(x,j,i); if (gequal0(b)) continue;
2717 175 : a = gcoeff(x,i,i);
2718 175 : d = gbezout_step(&b, &a, &v, &u, vx);
2719 651 : for (k = 1; k < i; k++)
2720 : {
2721 476 : GEN t = gadd(gmul(u,gcoeff(x,i,k)),gmul(v,gcoeff(x,j,k)));
2722 476 : gcoeff(x,j,k) = gsub(gmul(a,gcoeff(x,j,k)),gmul(b,gcoeff(x,i,k)));
2723 476 : gcoeff(x,i,k) = t;
2724 : }
2725 175 : gcoeff(x,j,i) = gen_0;
2726 175 : gcoeff(x,i,i) = d;
2727 175 : if (all) update(u,v,a,b,(GEN*)(U+i),(GEN*)(U+j));
2728 175 : c = 1;
2729 : }
2730 315 : if (!c)
2731 : {
2732 196 : k = gsnf_no_divide(x, i, vx);
2733 196 : if (!k) break;
2734 :
2735 245 : for (j=1; j<=i; j++)
2736 196 : gcoeff(x,i,j) = gadd(gcoeff(x,i,j),gcoeff(x,k,j));
2737 49 : if (all) gel(U,i) = gadd(gel(U,i),gel(U,k));
2738 : }
2739 168 : if (gc_needed(av,1))
2740 : {
2741 0 : if (DEBUGMEM>1) pari_warn(warnmem,"gsmithall");
2742 0 : gerepileall(av, all? 3: 1, &x, &U, &V);
2743 : }
2744 : }
2745 : }
2746 245 : for (k=1; k<=n; k++)
2747 : {
2748 196 : GEN d, T = normalize_as_RgX(gcoeff(x,k,k), vx, &d);
2749 196 : if (gequal0(T)) continue;
2750 182 : if (all && !gequal1(d)) gel(V,k) = RgC_Rg_div(gel(V,k), d);
2751 182 : gcoeff(x,k,k) = T;
2752 : }
2753 49 : z = all? mkvec3(shallowtrans(U), V, x): RgM_diagonal_shallow(x);
2754 49 : return gerepilecopy(av, z);
2755 : }
2756 :
2757 : GEN
2758 469 : matsnf0(GEN x,long flag)
2759 : {
2760 469 : pari_sp av = avma;
2761 469 : if (flag > 7) pari_err_FLAG("matsnf");
2762 469 : if (typ(x) == t_VEC && flag & 4) return smithclean(x);
2763 469 : if (typ(x)!=t_MAT) pari_err_TYPE("matsnf",x);
2764 469 : if (RgM_is_ZM(x)) x = flag&1 ? smithall(x): smith(x);
2765 49 : else x = RgXM_snf(x, flag&1);
2766 469 : if (flag & 4) x = gerepileupto(av, smithclean(x));
2767 469 : return x;
2768 : }
2769 : GEN
2770 0 : gsmith(GEN x) { return RgXM_snf(x,0); }
2771 : GEN
2772 0 : gsmithall(GEN x) { return RgXM_snf(x,1); }
2773 :
2774 : /* H is a relation matrix, either in HNF or a t_VEC (diagonal HNF) */
2775 : static GEN
2776 943579 : snf_group(GEN H, GEN D, GEN *newU, GEN *newUi)
2777 : {
2778 : long i, j, l;
2779 :
2780 943579 : ZM_snfclean(D, newU? *newU: NULL, newUi? *newUi: NULL);
2781 943550 : l = lg(D);
2782 943550 : if (newU) {
2783 821820 : GEN U = *newU;
2784 2120845 : for (i = 1; i < l; i++)
2785 : {
2786 1299506 : GEN d = gel(D,i), d2 = shifti(d, 1);
2787 5080597 : for (j = 1; j < lg(U); j++)
2788 3781572 : gcoeff(U,i,j) = centermodii(gcoeff(U,i,j), d, d2);
2789 : }
2790 821339 : *newU = U;
2791 : }
2792 943069 : if (newUi && l > 1)
2793 : { /* UHV=D -> U^-1 = (HV)D^-1 -> U^-1 = H(VD^-1 mod 1) mod H */
2794 : /* Ui = ZM_inv(U, NULL); setlg(Ui, l); */
2795 856153 : GEN V = *newUi, Ui;
2796 856153 : int Hvec = (typ(H) == t_VEC);
2797 2334445 : for (i = 1; i < l; i++) gel(V,i) = FpC_red(gel(V,i), gel(D,i));
2798 855952 : if (!Hvec)
2799 : {
2800 555342 : if (ZM_isdiagonal(H)) { H = RgM_diagonal_shallow(H); Hvec = 1; }
2801 : }
2802 856265 : Ui = Hvec? ZM_diag_mul(H, V): ZM_mul(H, V);
2803 2334329 : for (i = 1; i < l; i++) gel(Ui,i) = ZC_Z_divexact(gel(Ui,i), gel(D,i));
2804 855959 : *newUi = Hvec? ZM_ZV_mod(Ui, H): ZM_hnfrem(Ui, H);
2805 : }
2806 943054 : return D;
2807 : }
2808 : /* H relation matrix among row of generators g in HNF. Let URV = D its SNF,
2809 : * newU R newV = newD its clean SNF (no 1 in Dnew). Return the diagonal of
2810 : * newD, newU and newUi such that 1/U = (newUi, ?).
2811 : * Rationale: let (G,0) = g Ui be the new generators then
2812 : * 0 = G U R --> G D = 0, g = G newU and G = g newUi */
2813 : GEN
2814 642360 : ZM_snf_group(GEN H, GEN *newU, GEN *newUi)
2815 : {
2816 642360 : GEN D = ZM_snfall_i(H, newU, newUi, 1 + 2);
2817 642446 : return snf_group(H, D, newU, newUi);
2818 : }
2819 :
2820 : /* D a ZV: SNF for matdiagonal(D). Faster because we only ensure elementary
2821 : * divisors condition: d[n] | ... | d[1] and need not convert D to matrix form*/
2822 : GEN
2823 301148 : ZV_snfall(GEN D, GEN *pU, GEN *pV)
2824 : {
2825 301148 : pari_sp av = avma;
2826 301148 : long j, n = lg(D)-1;
2827 301148 : GEN U = pU? matid(n): NULL;
2828 301151 : GEN V = pV? matid(n): NULL;
2829 : GEN p;
2830 :
2831 301158 : D = leafcopy(D);
2832 960797 : for (j = n; j > 0; j--)
2833 : {
2834 659646 : GEN b = gel(D,j);
2835 659646 : if (signe(b) < 0)
2836 : {
2837 0 : gel(D,j) = negi(b);
2838 0 : if (V) ZV_togglesign(gel(V,j));
2839 : }
2840 : }
2841 : /* entries are nonnegative integers */
2842 301151 : p = gen_indexsort(D, NULL, &negcmpii);
2843 301150 : D = vecpermute(D, p);
2844 301152 : if (U) U = vecpermute(U, p);
2845 301152 : if (V) V = vecpermute(V, p);
2846 : /* entries are sorted by decreasing value */
2847 960753 : for (j = n; j > 0; j--)
2848 : {
2849 659599 : GEN b = gel(D,j);
2850 : long i;
2851 1180543 : for (i = j-1; i > 0; i--)
2852 : { /* invariant: a >= b. If au+bv = d is a Bezout relation, A=a/d and B=b/d
2853 : * we have [B,-A;u,v]*diag(a,b)*[1-u*A,1; -u*A,1]] = diag(Ab, d) */
2854 533125 : GEN a = gel(D,i), u,v, d = bezout(a,b, &u,&v), A, Wi, Wj;
2855 533115 : if (equalii(d,b)) continue;
2856 70295 : A = diviiexact(a,d);
2857 70295 : if (V)
2858 : {
2859 70253 : GEN t = mulii(u,A);
2860 70251 : Wi = ZC_lincomb(subui(1,t), negi(t), gel(V,i), gel(V,j));
2861 70253 : Wj = ZC_add(gel(V,i), gel(V,j));
2862 70255 : gel(V,i) = Wi;
2863 70255 : gel(V,j) = Wj;
2864 : }
2865 70297 : if (U)
2866 : {
2867 70297 : GEN B = diviiexact(b,d);
2868 70298 : Wi = ZC_lincomb(B, negi(A), gel(U,i), gel(U,j));
2869 70297 : Wj = ZC_lincomb(u, v, gel(U,i), gel(U,j));
2870 70297 : gel(U,i) = Wi;
2871 70297 : gel(U,j) = Wj;
2872 : }
2873 70297 : gel(D,i) = mulii(A,b); /* lcm(a,b) */
2874 70295 : gel(D,j) = d; /* gcd(a,b) */
2875 70295 : b = gel(D,j); if (equali1(b)) break;
2876 : }
2877 : }
2878 301154 : snf_pile(av, &D,&U,&V);
2879 301158 : if (U) *pU = shallowtrans(U);
2880 301158 : if (V) *pV = V;
2881 301158 : return D;
2882 : }
2883 : GEN
2884 301149 : ZV_snf_group(GEN d, GEN *newU, GEN *newUi)
2885 : {
2886 301149 : GEN D = ZV_snfall(d, newU, newUi);
2887 301157 : return snf_group(d, D, newU, newUi);
2888 : }
2889 :
2890 : /* D a vector of elementary divisors. Truncate (setlg) to leave out trivial
2891 : * entries (= 1) */
2892 : void
2893 0 : ZV_snf_trunc(GEN D)
2894 : {
2895 0 : long i, l = lg(D);
2896 0 : for (i = 1; i < l; i++)
2897 0 : if (is_pm1(gel(D,i))) { setlg(D,i); break; }
2898 0 : }
2899 :
2900 : long
2901 0 : zv_snf_rank(GEN D, ulong p)
2902 : {
2903 0 : long i, l = lg(D);
2904 0 : if (!p) return l - 1;
2905 0 : for (i = 1; i < l; i++)
2906 0 : if (D[i] % p) break;
2907 0 : return i - 1;
2908 : }
2909 : long
2910 49 : ZV_snf_rank_u(GEN D, ulong p)
2911 : {
2912 49 : long i, l = lg(D);
2913 49 : while (l > 1 && D[l-1] == 1) l--;
2914 49 : if (!p) return l - 1;
2915 49 : if (p == 2)
2916 : {
2917 49 : for (i = 1; i < l; i++)
2918 42 : if (mpodd(gel(D,i))) break;
2919 : }
2920 35 : else if (!(p & (p-1)))
2921 : { /* power of 2 */
2922 14 : long n = vals(p);
2923 28 : for (i = 1; i < l; i++)
2924 28 : if (umodi2n(gel(D,i), n)) break;
2925 : }
2926 : else
2927 : {
2928 49 : for (i = 1; i < l; i++)
2929 42 : if (umodiu(gel(D,i), p)) break;
2930 : }
2931 49 : return i - 1;
2932 : }
2933 : long
2934 91 : ZV_snf_rank(GEN D, GEN p)
2935 : {
2936 91 : long i, l = lg(D);
2937 91 : if (lgefint(p) == 3) return ZV_snf_rank_u(D, p[2]);
2938 77 : while (l > 1 && equali1(gel(D, l-1))) l--;
2939 42 : if (!signe(p)) return l - 1;
2940 77 : for (i = 1; i < l; i++)
2941 70 : if (!dvdii(gel(D,i), p)) break;
2942 21 : return i - 1;
2943 : }
2944 : long
2945 154 : snfrank(GEN D, GEN p)
2946 : {
2947 : long i, l;
2948 154 : if (typ(D) != t_VEC) pari_err_TYPE("snfrank", D);
2949 154 : if (!p) p = gen_0;
2950 154 : l = lg(D);
2951 154 : if (l == 4 && typ(gel(D,3)) == t_MAT)
2952 : { /* from matsnf(,1) */
2953 14 : pari_sp av = avma;
2954 : long z;
2955 : GEN v;
2956 14 : D = gel(D,3); l = lg(D);
2957 14 : if (l == 1) return 0;
2958 14 : z = lgcols(D) - l; /* missing columns of 0s */
2959 14 : if (z < 0) pari_err_TYPE("snfrank", D);
2960 14 : v = cgetg(l, t_VEC);
2961 35 : for (i = 1; i < l; i++) gel(v, i) = gcoeff(D, i + z, i);
2962 14 : return gc_long(av, z + snfrank(v, p));
2963 : }
2964 140 : switch(typ(p))
2965 : {
2966 98 : case t_INT:
2967 98 : if (RgV_is_ZV(D)) return ZV_snf_rank(D, p);
2968 7 : if (!signe(p)) break; /* allow p = 0 */
2969 0 : pari_err_TYPE("snfrank", D);
2970 42 : case t_POL: break;
2971 0 : default: pari_err_TYPE("snfrank", p);
2972 : }
2973 175 : while (l > 1 && gequal1(gel(D, l-1))) l--;
2974 49 : if (gequal0(p)) return l - 1;
2975 112 : for (i = 1; i < l; i++)
2976 91 : if (!gdvd(gel(D,i), p)) break;
2977 42 : return i - 1;
2978 : }
|
