1 /* $Id: alglin2.c 11336 2008-11-22 17:33:45Z kb $
2
3 Copyright (C) 2000 The PARI group.
4
5 This file is part of the PARI/GP package.
6
7 PARI/GP is free software; you can redistribute it and/or modify it under the
8 terms of the GNU General Public License as published by the Free Software
9 Foundation. It is distributed in the hope that it will be useful, but WITHOUT
10 ANY WARRANTY WHATSOEVER.
11
12 Check the License for details. You should have received a copy of it, along
13 with the package; see the file 'COPYING'. If not, write to the Free Software
14 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
15
16 /********************************************************************/
17 /** **/
18 /** LINEAR ALGEBRA **/
19 /** (second part) **/
20 /** **/
21 /********************************************************************/
22 #include "pari.h"
23 #include "paripriv.h"
24 /*******************************************************************/
25 /* */
26 /* CHARACTERISTIC POLYNOMIAL */
27 /* */
28 /*******************************************************************/
29
30 GEN
charpoly0(GEN x,long v,long flag)31 charpoly0(GEN x, long v, long flag)
32 {
33 if (v<0) v = 0;
34 switch(flag)
35 {
36 case 0: return caradj(x,v,NULL);
37 case 1: return caract(x,v);
38 case 2: return carhess(x,v);
39 }
40 pari_err(flagerr,"charpoly"); return NULL; /* not reached */
41 }
42
43 /* (v - x)^d */
44 static GEN
caract_const(GEN x,long v,long d)45 caract_const(GEN x, long v, long d)
46 {
47 pari_sp av = avma;
48 return gerepileupto(av, gpowgs(gadd(pol_x[v], gneg_i(x)), d));
49 }
50
51 static GEN
caract2_i(GEN p,GEN x,long v,GEN (subres_f)(GEN,GEN,GEN *))52 caract2_i(GEN p, GEN x, long v, GEN (subres_f)(GEN,GEN,GEN*))
53 {
54 pari_sp av = avma;
55 long d = degpol(p), dx;
56 GEN ch, L;
57
58 if (typ(x) != t_POL) return caract_const(x, v, d);
59 dx = degpol(x);
60 if (dx <= 0) return dx? monomial(gen_1, d, v): caract_const(gel(x,2), v, d);
61
62 x = gneg_i(x);
63 if (varn(x) == MAXVARN) { setvarn(x, 0); p = shallowcopy(p); setvarn(p, 0); }
64 gel(x,2) = gadd(gel(x,2), pol_x[MAXVARN]);
65 ch = subres_f(p, x, NULL);
66 if (v != MAXVARN)
67 {
68 if (typ(ch) == t_POL && varn(ch) == MAXVARN)
69 setvarn(ch, v);
70 else
71 ch = gsubst(ch, MAXVARN, pol_x[v]);
72 }
73 L = leading_term(ch);
74 if (!gcmp1(L)) ch = gdiv(ch, L);
75 return gerepileupto(av, ch);
76 }
77
78 /* return caract(Mod(x,p)) in variable v */
79 GEN
caract2(GEN p,GEN x,long v)80 caract2(GEN p, GEN x, long v)
81 {
82 return caract2_i(p,x,v, subresall);
83 }
84 GEN
caractducos(GEN p,GEN x,long v)85 caractducos(GEN p, GEN x, long v)
86 {
87 return caract2_i(p,x,v, (GEN (*)(GEN,GEN,GEN*))resultantducos);
88 }
89
90 /* characteristic pol. Easy cases. Return NULL in case it's not so easy. */
91 static GEN
easychar(GEN x,long v,GEN * py)92 easychar(GEN x, long v, GEN *py)
93 {
94 pari_sp av;
95 long lx;
96 GEN p1;
97
98 switch(typ(x))
99 {
100 case t_INT: case t_REAL: case t_INTMOD:
101 case t_FRAC: case t_PADIC:
102 p1=cgetg(4,t_POL);
103 p1[1]=evalsigne(1) | evalvarn(v);
104 gel(p1,2) = gneg(x); gel(p1,3) = gen_1;
105 if (py) *py = mkmat(mkcolcopy(x));
106 return p1;
107
108 case t_COMPLEX: case t_QUAD:
109 if (py) pari_err(typeer,"easychar");
110 p1 = cgetg(5,t_POL);
111 p1[1] = evalsigne(1) | evalvarn(v);
112 gel(p1,2) = gnorm(x); av = avma;
113 gel(p1,3) = gerepileupto(av, gneg(gtrace(x)));
114 gel(p1,4) = gen_1; return p1;
115
116 case t_POLMOD:
117 if (py) pari_err(typeer,"easychar");
118 return caract2(gel(x,1), gel(x,2), v);
119
120 case t_MAT:
121 lx=lg(x);
122 if (lx==1)
123 {
124 if (py) *py = cgetg(1,t_MAT);
125 return pol_1[v];
126 }
127 if (lg(x[1]) != lx) break;
128 return NULL;
129 }
130 pari_err(mattype1,"easychar");
131 return NULL; /* not reached */
132 }
133
134 GEN
caract(GEN x,long v)135 caract(GEN x, long v)
136 {
137 long k, n;
138 pari_sp av=avma;
139 GEN p1, p2, p3, x_k, Q;
140
141 if ((p1 = easychar(x,v,NULL))) return p1;
142
143 p1 = gen_0; Q = p2 = gen_1; n = lg(x)-1;
144 x_k = monomial(gen_1, 1, v);
145 for (k=0; k<=n; k++)
146 {
147 GEN mk = stoi(-k);
148 gel(x_k,2) = mk;
149 p3 = det(gaddmat_i(mk, x));
150 p1 = gadd(gmul(p1, x_k), gmul(gmul(p2, p3), Q));
151 if (k == n) break;
152
153 Q = gmul(Q, x_k);
154 p2 = divis(mulsi(k-n,p2), k+1); /* (-1)^{k} binomial(n,k) */
155 }
156 return gerepileupto(av, gdiv(p1, mpfact(n)));
157 }
158
159 GEN
caradj0(GEN x,long v)160 caradj0(GEN x, long v)
161 {
162 return caradj(x,v,NULL);
163 }
164
165 /* assume x square matrice */
166 static GEN
mattrace(GEN x)167 mattrace(GEN x)
168 {
169 pari_sp av;
170 long i, lx = lg(x);
171 GEN t;
172 if (lx < 3) return lx == 1? gen_0: gcopy(gcoeff(x,1,1));
173 av = avma; t = gcoeff(x,1,1);
174 for (i = 2; i < lx; i++) t = gadd(t, gcoeff(x,i,i));
175 return gerepileupto(av, t);
176 }
177
178 /* Using traces: return the characteristic polynomial of x (in variable v).
179 * If py != NULL, the adjoint matrix is put there. */
180 GEN
caradj(GEN x,long v,GEN * py)181 caradj(GEN x, long v, GEN *py)
182 {
183 pari_sp av, av0;
184 long i, k, l;
185 GEN p, y, t;
186
187 if ((p = easychar(x, v, py))) return p;
188
189 l = lg(x); av0 = avma;
190 p = cgetg(l+2,t_POL); p[1] = evalsigne(1) | evalvarn(v);
191 gel(p,l+1) = gen_1;
192 if (l == 1) { if (py) *py = cgetg(1,t_MAT); return p; }
193 av = avma; t = gerepileupto(av, gneg(mattrace(x)));
194 gel(p,l) = t;
195 if (l == 2) { if (py) *py = matid(1); return p; }
196 if (l == 3) {
197 GEN a = gcoeff(x,1,1), b = gcoeff(x,1,2);
198 GEN c = gcoeff(x,2,1), d = gcoeff(x,2,2);
199 if (py) {
200 y = cgetg(3, t_MAT);
201 gel(y,1) = mkcol2(gcopy(d), gneg(c));
202 gel(y,2) = mkcol2(gneg(b), gcopy(a));
203 *py = y;
204 }
205 av = avma;
206 gel(p,2) = gerepileupto(av, gadd(gmul(a,d), gneg(gmul(b,c))));
207 return p;
208 }
209 av = avma; y = shallowcopy(x);
210 for (i = 1; i < l; i++) gcoeff(y,i,i) = gadd(gcoeff(y,i,i), t);
211 for (k = 2; k < l-1; k++)
212 {
213 GEN y0 = y;
214 y = gmul(y, x);
215 t = gdivgs(mattrace(y), -k);
216 for (i = 1; i < l; i++) gcoeff(y,i,i) = gadd(gcoeff(y,i,i), t);
217 y = gclone(y);
218 /* beware: since y is a clone and t computed from it some components
219 * may be out of stack (eg. INTMOD/POLMOD) */
220 gel(p,l-k+1) = gerepileupto(av, gcopy(t)); av = avma;
221 if (k > 2) gunclone(y0);
222 }
223 t = gen_0;
224 for (i=1; i<l; i++) t = gadd(t, gmul(gcoeff(x,1,i),gcoeff(y,i,1)));
225 gel(p,2) = gerepileupto(av, gneg(t));
226 i = gvar2(p);
227 if (i == v) pari_err(talker,"incorrect variable in caradj");
228 if (i < v) p = gerepileupto(av0, poleval(p, pol_x[v]));
229 if (py) *py = (l & 1)? gneg(y): gcopy(y);
230 gunclone(y); return p;
231 }
232
233 GEN
adj(GEN x)234 adj(GEN x)
235 {
236 GEN y;
237 (void)caradj(x,MAXVARN,&y); return y;
238 }
239
240 /*******************************************************************/
241 /* */
242 /* MINIMAL POLYNOMIAL */
243 /* */
244 /*******************************************************************/
245
246 static GEN
easymin(GEN x,long v)247 easymin(GEN x, long v)
248 {
249 pari_sp ltop=avma;
250 GEN G, R, dR;
251 if (typ(x)==t_POLMOD && !issquarefree(gel(x,1)))
252 return NULL;
253 R=easychar(x, v, NULL);
254 if (!R) return R;
255 dR=derivpol(R);
256 if (!lgpol(dR)) {avma=ltop; return NULL;}
257 G=srgcd(R,dR);
258 G=gdiv(G,leading_term(G));
259 G=gdeuc(R,G);
260 return gerepileupto(ltop,G);
261 }
262
263 GEN
minpoly(GEN x,long v)264 minpoly(GEN x, long v)
265 {
266 pari_sp ltop=avma;
267 GEN P;
268 if (v<0) v = 0;
269 P=easymin(x,v);
270 if (P) return P;
271 if (typ(x)==t_POLMOD)
272 {
273 P = gcopy(RgXQ_minpoly_naive(gel(x,2), gel(x,1)));
274 setvarn(P,v);
275 return gerepileupto(ltop,P);
276 }
277 if (typ(x)!=t_MAT) pari_err(typeer,"minpoly");
278 if (lg(x) == 1) return pol_1[v];
279 return gerepilecopy(ltop,gel(matfrobenius(x,1,v),1));
280 }
281
282 /*******************************************************************/
283 /* */
284 /* HESSENBERG FORM */
285 /* */
286 /*******************************************************************/
287 GEN
hess(GEN x)288 hess(GEN x)
289 {
290 pari_sp av = avma, av2 = avma, lim;
291 long lx = lg(x), m, i, j;
292 GEN p, p1, p2;
293
294 if (typ(x) != t_MAT) pari_err(mattype1,"hess");
295 if (lx == 1) return cgetg(1,t_MAT);
296 if (lg(x[1]) != lx) pari_err(mattype1,"hess");
297 x = shallowcopy(x); lim = stack_lim(av,1);
298
299 for (m=2; m<lx-1; m++)
300 for (i=m+1; i<lx; i++)
301 {
302 p = gcoeff(x,i,m-1);
303 if (gcmp0(p)) continue;
304
305 for (j=m-1; j<lx; j++) swap(gcoeff(x,i,j), gcoeff(x,m,j));
306 lswap(x[i], x[m]); p = ginv(p);
307 for (i=m+1; i<lx; i++)
308 {
309 p1 = gcoeff(x,i,m-1);
310 if (gcmp0(p1)) continue;
311
312 p1 = gmul(p1,p); p2 = gneg_i(p1);
313 gcoeff(x,i,m-1) = gen_0;
314 for (j=m; j<lx; j++)
315 gcoeff(x,i,j) = gadd(gcoeff(x,i,j), gmul(p2,gcoeff(x,m,j)));
316 for (j=1; j<lx; j++)
317 gcoeff(x,j,m) = gadd(gcoeff(x,j,m), gmul(p1,gcoeff(x,j,i)));
318 }
319 if (low_stack(lim, stack_lim(av,1)))
320 {
321 if (DEBUGMEM>1) pari_warn(warnmem,"hess, m = %ld", m);
322 x = gerepilecopy(av2, x);
323 }
324 break;
325 }
326 return gerepilecopy(av,x);
327 }
328
329 GEN
carhess(GEN x,long v)330 carhess(GEN x, long v)
331 {
332 pari_sp av;
333 long lx, r, i;
334 GEN y, H, X_h;
335
336 if ((H = easychar(x,v,NULL))) return H;
337
338 lx = lg(x); av = avma; y = cgetg(lx+1, t_VEC);
339 gel(y,1) = pol_1[v]; H = hess(x);
340 X_h = monomial(gen_1, 1, v);
341 for (r = 1; r < lx; r++)
342 {
343 GEN p3 = gen_1, p4 = gen_0;
344 for (i = r-1; i; i--)
345 {
346 p3 = gmul(p3, gcoeff(H,i+1,i));
347 p4 = gadd(p4, gmul(gmul(p3,gcoeff(H,i,r)), gel(y,i)));
348 }
349 gel(X_h,2) = gneg(gcoeff(H,r,r));
350 gel(y,r+1) = gsub(gmul(gel(y,r), X_h), p4);
351 }
352 return gerepileupto(av, gel(y,lx));
353 }
354
355 /*******************************************************************/
356 /* */
357 /* NORM */
358 /* */
359 /*******************************************************************/
360 GEN
cxnorm(GEN x)361 cxnorm(GEN x) { return gadd(gsqr(gel(x,1)), gsqr(gel(x,2))); }
362 GEN
quadnorm(GEN x)363 quadnorm(GEN x)
364 {
365 GEN u = gel(x,3), v = gel(x,2);
366 GEN X = gel(x,1), b = gel(X,3), c = gel(X,2);
367 GEN z = signe(b)? gmul(v, gadd(u,v)): gsqr(v);
368 return gadd(z, gmul(c, gsqr(u)));
369 }
370
371 GEN
RgXQ_norm(GEN x,GEN T)372 RgXQ_norm(GEN x, GEN T)
373 {
374 pari_sp av;
375 GEN L, y;
376 if (typ(x) != t_POL) return gpowgs(x, degpol(T));
377
378 av = avma; y = subres(T, x);
379 L = leading_term(T);
380 if (gcmp1(L) || gcmp0(x)) return y;
381 return gerepileupto(av, gdiv(y, gpowgs(L, degpol(x))));
382 }
383
384 GEN
gnorm(GEN x)385 gnorm(GEN x)
386 {
387 pari_sp av;
388 long lx, i, tx=typ(x);
389 GEN y;
390
391 switch(tx)
392 {
393 case t_INT: return sqri(x);
394 case t_REAL: return mulrr(x,x);
395 case t_FRAC: return gsqr(x);
396 case t_COMPLEX: av = avma; return gerepileupto(av, cxnorm(x));
397 case t_QUAD: av = avma; return gerepileupto(av, quadnorm(x));
398
399 case t_POL: case t_SER: case t_RFRAC: av = avma;
400 return gerepileupto(av, greal(gmul(gconj(x),x)));
401
402 case t_POLMOD: return RgXQ_norm(gel(x,2), gel(x,1));
403
404 case t_VEC: case t_COL: case t_MAT:
405 lx=lg(x); y=cgetg(lx,tx);
406 for (i=1; i<lx; i++) gel(y,i) = gnorm(gel(x,i));
407 return y;
408 }
409 pari_err(typeer,"gnorm");
410 return NULL; /* not reached */
411 }
412
413 GEN
gnorml2(GEN x)414 gnorml2(GEN x)
415 {
416 pari_sp av,lim;
417 GEN y;
418 long i,tx=typ(x),lx;
419
420 if (! is_matvec_t(tx)) return gnorm(x);
421 lx=lg(x); if (lx==1) return gen_0;
422
423 av=avma; lim = stack_lim(av,1); y = gnorml2(gel(x,1));
424 for (i=2; i<lx; i++)
425 {
426 y = gadd(y, gnorml2(gel(x,i)));
427 if (low_stack(lim, stack_lim(av,1)))
428 {
429 if (DEBUGMEM>1) pari_warn(warnmem,"gnorml2");
430 y = gerepileupto(av, y);
431 }
432 }
433 return gerepileupto(av,y);
434 }
435
436 GEN
QuickNormL2(GEN x,long prec)437 QuickNormL2(GEN x, long prec)
438 {
439 pari_sp av = avma;
440 GEN y = gmul(x, real_1(prec));
441 if (typ(x) == t_POL) *++y = evaltyp(t_VEC) | evallg(lg(x)-1);
442 return gerepileupto(av, gnorml2(y));
443 }
444
445 GEN
gnorml1(GEN x,long prec)446 gnorml1(GEN x,long prec)
447 {
448 pari_sp av = avma;
449 long lx,i;
450 GEN s;
451 switch(typ(x))
452 {
453 case t_INT: case t_REAL: case t_COMPLEX: case t_FRAC: case t_QUAD:
454 return gabs(x,prec);
455
456 case t_POL:
457 lx = lg(x); s = gen_0;
458 for (i=2; i<lx; i++) s = gadd(s, gabs(gel(x,i),prec));
459 break;
460
461 case t_VEC: case t_COL: case t_MAT:
462 lx = lg(x); s = gen_0;
463 for (i=1; i<lx; i++) s = gadd(s, gabs(gel(x,i),prec));
464 break;
465
466 default: pari_err(typeer,"gnorml1");
467 return NULL; /* not reached */
468 }
469 return gerepileupto(av, s);
470 }
471
472 GEN
QuickNormL1(GEN x,long prec)473 QuickNormL1(GEN x,long prec)
474 {
475 pari_sp av = avma;
476 long lx,i;
477 GEN p1,p2,s;
478 switch(typ(x))
479 {
480 case t_INT: case t_REAL: case t_FRAC:
481 return gabs(x,prec);
482
483 case t_INTMOD: case t_PADIC: case t_POLMOD: case t_SER: case t_RFRAC:
484 return gcopy(x);
485
486 case t_COMPLEX:
487 p1=gabs(gel(x,1),prec); p2=gabs(gel(x,2),prec);
488 return gerepileupto(av, gadd(p1,p2));
489
490 case t_QUAD:
491 p1=gabs(gel(x,2),prec); p2=gabs(gel(x,3),prec);
492 return gerepileupto(av, gadd(p1,p2));
493
494 case t_POL:
495 lx=lg(x); s=gen_0;
496 for (i=2; i<lx; i++) s=gadd(s,QuickNormL1(gel(x,i),prec));
497 return gerepileupto(av, s);
498
499 case t_VEC: case t_COL: case t_MAT:
500 lx=lg(x); s=gen_0;
501 for (i=1; i<lx; i++) s=gadd(s,QuickNormL1(gel(x,i),prec));
502 return gerepileupto(av, s);
503 }
504 pari_err(typeer,"QuickNormL1");
505 return NULL; /* not reached */
506 }
507
508 /*******************************************************************/
509 /* */
510 /* CONJUGATION */
511 /* */
512 /*******************************************************************/
513
514 /* lift( conj(Mod(x, y)) ), assuming degpol(y) = 2, degpol(x) < 2 */
515 GEN
quad_polmod_conj(GEN x,GEN y)516 quad_polmod_conj(GEN x, GEN y)
517 {
518 GEN z, u, v, a, b;
519 pari_sp av;
520 if (typ(x) != t_POL || varn(x) != varn(y) || degpol(x) <= 0)
521 return gcopy(x);
522 a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
523 b = gel(y,3); v = gel(x,2);
524 z = cgetg(4, t_POL); z[1] = x[1]; av = avma;
525 gel(z,2) = gerepileupto(av, gadd(v, gdiv(gmul(u,gneg(b)), a)));
526 gel(z,3) = gneg(u); return z;
527 }
528 GEN
quad_polmod_norm(GEN x,GEN y)529 quad_polmod_norm(GEN x, GEN y)
530 {
531 GEN z, u, v, a, b, c;
532 pari_sp av;
533 if (typ(x) != t_POL || varn(x) != varn(y) || degpol(x) <= 0)
534 return gsqr(x);
535 a = gel(y,4); u = gel(x,3); /*Mod(ux + v, ax^2 + bx + c)*/
536 b = gel(y,3); v = gel(x,2);
537 c = gel(y,2); av = avma;
538 z = gmul(u, gadd( gmul(c, u), gmul(gneg(b), v)));
539 if (!gcmp1(a)) z = gdiv(z, a);
540 return gerepileupto(av, gadd(z, gsqr(v)));
541 }
542
543 GEN
gconj(GEN x)544 gconj(GEN x)
545 {
546 long lx, i, tx = typ(x);
547 GEN z;
548
549 switch(tx)
550 {
551 case t_INT: case t_REAL: case t_INTMOD: case t_FRAC: case t_PADIC:
552 return gcopy(x);
553
554 case t_COMPLEX:
555 z = cgetg(3,t_COMPLEX);
556 gel(z,1) = gcopy(gel(x,1));
557 gel(z,2) = gneg(gel(x,2));
558 break;
559
560 case t_QUAD:
561 z = cgetg(4,t_QUAD);
562 gel(z,1) = gcopy(gel(x,1));
563 gel(z,2) = gcmp0(gmael(x,1,3))? gcopy(gel(x,2))
564 : gadd(gel(x,2), gel(x,3));
565 gel(z,3) = gneg(gel(x,3));
566 break;
567
568 case t_POL: case t_SER: case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
569 z = init_gen_op(x, tx, &lx, &i);
570 for (; i<lx; i++) gel(z,i) = gconj(gel(x,i));
571 break;
572
573 case t_POLMOD:
574 {
575 GEN z, X = gel(x,1);
576 long d = degpol(X);
577 if (d < 2) return gcopy(x);
578 if (d == 2) {
579 z = cgetg(3, t_POLMOD);
580 gel(z,1) = gcopy(X);
581 gel(z,2) = quad_polmod_conj(gel(x,2), X); return z;
582 }
583 }
584 default:
585 pari_err(typeer,"gconj");
586 return NULL; /* not reached */
587 }
588 return z;
589 }
590
591 GEN
conjvec(GEN x,long prec)592 conjvec(GEN x,long prec)
593 {
594 pari_sp av,tetpil;
595 long lx,s,i,tx=typ(x);
596 GEN z,y,p1,p2,p;
597
598 switch(tx)
599 {
600 case t_INT: case t_INTMOD: case t_FRAC:
601 return mkcolcopy(x);
602
603 case t_COMPLEX: case t_QUAD:
604 z=cgetg(3,t_COL); gel(z,1) = gcopy(x); gel(z,2) = gconj(x); break;
605
606 case t_VEC: case t_COL:
607 lx=lg(x); z=cgetg(lx,t_MAT);
608 for (i=1; i<lx; i++) gel(z,i) = conjvec(gel(x,i),prec);
609 if (lx == 1) break;
610 s = lg(z[1]);
611 for (i=2; i<lx; i++)
612 if (lg(z[i])!=s) pari_err(talker,"incompatible field degrees in conjvec");
613 break;
614
615 case t_POLMOD:
616 y=gel(x,1); lx=lg(y);
617 if (lx<=3) return cgetg(1,t_COL);
618 av=avma; p=NULL;
619 for (i=2; i<lx; i++)
620 {
621 tx=typ(y[i]);
622 if (tx==t_INTMOD) p = gmael(y,i,1);
623 else
624 if (!is_rational_t(tx))
625 pari_err(talker,"not a rational polynomial in conjvec");
626 }
627 if (!p)
628 {
629 p1=roots(y,prec); x = gel(x,2); tetpil=avma;
630 z=cgetg(lx-2,t_COL);
631 for (i=1; i<=lx-3; i++)
632 {
633 p2=gel(p1,i); if (gcmp0(gel(p2,2))) p2 = gel(p2,1);
634 gel(z,i) = poleval(x, p2);
635 }
636 return gerepile(av,tetpil,z);
637 }
638 z=cgetg(lx-2,t_COL); gel(z,1) = gcopy(x);
639 for (i=2; i<=lx-3; i++) gel(z,i) = gpow(gel(z,i-1),p,prec);
640 break;
641
642 default:
643 pari_err(typeer,"conjvec");
644 return NULL; /* not reached */
645 }
646 return z;
647 }
648
649 /*******************************************************************/
650 /* */
651 /* TRACES */
652 /* */
653 /*******************************************************************/
654
655 GEN
assmat(GEN x)656 assmat(GEN x)
657 {
658 long lx,i,j;
659 GEN y,p1,p2;
660
661 if (typ(x)!=t_POL) pari_err(notpoler,"assmat");
662 if (gcmp0(x)) pari_err(zeropoler,"assmat");
663
664 lx=lg(x)-2; y=cgetg(lx,t_MAT);
665 if (lx == 1) return y;
666 for (i=1; i<lx-1; i++)
667 {
668 p1=cgetg(lx,t_COL); gel(y,i) = p1;
669 for (j=1; j<lx; j++) gel(p1,j) = (j==(i+1))? gen_1: gen_0;
670 }
671 p1=cgetg(lx,t_COL); gel(y,i) = p1;
672 if (gcmp1(gel(x,lx+1)))
673 for (j=1; j<lx; j++)
674 gel(p1,j) = gneg(gel(x,j+1));
675 else
676 {
677 pari_sp av = avma;
678 p2 = gclone(gneg(gel(x,lx+1)));
679 avma = av;
680 for (j=1; j<lx; j++)
681 gel(p1,j) = gdiv(gel(x,j+1),p2);
682 gunclone(p2);
683 }
684 return y;
685 }
686
687 GEN
gtrace(GEN x)688 gtrace(GEN x)
689 {
690 pari_sp av;
691 long i, lx, tx = typ(x);
692 GEN y, z;
693
694 switch(tx)
695 {
696 case t_INT: case t_REAL: case t_FRAC:
697 return gmul2n(x,1);
698
699 case t_COMPLEX:
700 return gmul2n(gel(x,1),1);
701
702 case t_QUAD:
703 y = gel(x,1);
704 if (!gcmp0(gel(y,3)))
705 { /* assume quad. polynomial is either x^2 + d or x^2 - x + d */
706 av = avma;
707 return gerepileupto(av, gadd(gel(x,3), gmul2n(gel(x,2),1)));
708 }
709 return gmul2n(gel(x,2),1);
710
711 case t_POL: case t_SER:
712 lx = lg(x); y = cgetg(lx,tx); y[1] = x[1];
713 for (i=2; i<lx; i++) gel(y,i) = gtrace(gel(x,i));
714 return y;
715
716 case t_POLMOD:
717 y = gel(x,1); z = gel(x,2);
718 if (typ(z) != t_POL || varn(y) != varn(z)) return gmulsg(degpol(y), z);
719 av = avma;
720 return gerepileupto(av, quicktrace(z, polsym(y, degpol(y)-1)));
721
722 case t_RFRAC:
723 return gadd(x, gconj(x));
724
725 case t_VEC: case t_COL:
726 lx = lg(x); y = cgetg(lx,tx);
727 for (i=1; i<lx; i++) gel(y,i) = gtrace(gel(x,i));
728 return y;
729
730 case t_MAT:
731 lx = lg(x); if (lx == 1) return gen_0;
732 /*now lx >= 2*/
733 if (lx != lg(x[1])) pari_err(mattype1,"gtrace");
734 return mattrace(x);
735 }
736 pari_err(typeer,"gtrace");
737 return NULL; /* not reached */
738 }
739
740 /* Cholesky Decomposition for positive definite matrix a [matrix Q, Algo 2.7.6]
741 * If a is not positive definite:
742 * if flag is zero: raise an error
743 * else: return NULL. */
744 GEN
sqred1intern(GEN a)745 sqred1intern(GEN a)
746 {
747 pari_sp av = avma, lim=stack_lim(av,1);
748 GEN b,p;
749 long i,j,k, n = lg(a);
750
751 if (typ(a)!=t_MAT) pari_err(typeer,"sqred1");
752 if (n == 1) return cgetg(1, t_MAT);
753 if (lg(a[1])!=n) pari_err(mattype1,"sqred1");
754 b=cgetg(n,t_MAT);
755 for (j=1; j<n; j++)
756 {
757 GEN p1=cgetg(n,t_COL), p2=gel(a,j);
758
759 gel(b,j) = p1;
760 for (i=1; i<=j; i++) p1[i] = p2[i];
761 for ( ; i<n ; i++) gel(p1,i) = gen_0;
762 }
763 for (k=1; k<n; k++)
764 {
765 p = gcoeff(b,k,k);
766 if (gsigne(p)<=0) { avma = av; return NULL; } /* not positive definite */
767 p = ginv(p);
768 for (i=k+1; i<n; i++)
769 for (j=i; j<n; j++)
770 gcoeff(b,i,j) = gsub(gcoeff(b,i,j),
771 gmul(gmul(gcoeff(b,k,i),gcoeff(b,k,j)), p));
772 for (j=k+1; j<n; j++)
773 gcoeff(b,k,j) = gmul(gcoeff(b,k,j), p);
774 if (low_stack(lim, stack_lim(av,1)))
775 {
776 if (DEBUGMEM>1) pari_warn(warnmem,"sqred1");
777 b=gerepilecopy(av,b);
778 }
779 }
780 return gerepilecopy(av,b);
781 }
782
783 GEN
sqred1(GEN a)784 sqred1(GEN a)
785 {
786 GEN x = sqred1intern(a);
787 if (!x) pari_err(talker,"not a positive definite matrix in sqred1");
788 return x;
789 }
790
791 GEN
sqred3(GEN a)792 sqred3(GEN a)
793 {
794 pari_sp av = avma, lim = stack_lim(av,1);
795 long i,j,k,l, n = lg(a);
796 GEN p1,b;
797
798 if (typ(a)!=t_MAT) pari_err(typeer,"sqred3");
799 if (lg(a[1])!=n) pari_err(mattype1,"sqred3");
800 av=avma; b=cgetg(n,t_MAT);
801 for (j=1; j<n; j++)
802 {
803 p1=cgetg(n,t_COL); gel(b,j) = p1;
804 for (i=1; i<n; i++) gel(p1,i) = gen_0;
805 }
806 for (i=1; i<n; i++)
807 {
808 for (k=1; k<i; k++)
809 {
810 p1=gen_0;
811 for (l=1; l<k; l++)
812 p1=gadd(p1, gmul(gmul(gcoeff(b,l,l),gcoeff(b,k,l)), gcoeff(b,i,l)));
813 gcoeff(b,i,k) = gdiv(gsub(gcoeff(a,i,k),p1),gcoeff(b,k,k));
814 }
815 p1=gen_0;
816 for (l=1; l<i; l++)
817 p1=gadd(p1, gmul(gcoeff(b,l,l), gsqr(gcoeff(b,i,l))));
818 gcoeff(b,i,k) = gsub(gcoeff(a,i,i),p1);
819 if (low_stack(lim, stack_lim(av,1)))
820 {
821 if (DEBUGMEM>1) pari_warn(warnmem,"sqred3");
822 b=gerepilecopy(av,b);
823 }
824 }
825 return gerepilecopy(av,b);
826 }
827
828 /* Gauss reduction (arbitrary symetric matrix, only the part above the
829 * diagonal is considered). If signature is zero, return only the
830 * signature, in which case gsigne() should be defined for elements of a. */
831 static GEN
sqred2(GEN a,long signature)832 sqred2(GEN a, long signature)
833 {
834 GEN r, p;
835 pari_sp av, av1, lim;
836 long n, i, j, k, l, sp, sn, t;
837
838 if (typ(a) != t_MAT) pari_err(typeer,"sqred2");
839 n = lg(a); if (n > 1 && lg(a[1]) != n) pari_err(mattype1,"sqred2");
840 n--;
841
842 av = avma;
843 r = const_vecsmall(n, 1);
844 av1= avma; lim = stack_lim(av1,1);
845 a = shallowcopy(a);
846 t = n; sp = sn = 0;
847 while (t)
848 {
849 k=1; while (k<=n && (!r[k] || gcmp0(gcoeff(a,k,k)))) k++;
850 if (k<=n)
851 {
852 p = gcoeff(a,k,k);
853 if (signature) { /* don't check if signature = 0: gsigne may fail ! */
854 if (gsigne(p) > 0) sp++; else sn++;
855 }
856 r[k] = 0; t--;
857 for (j=1; j<=n; j++)
858 gcoeff(a,k,j) = r[j] ? gdiv(gcoeff(a,k,j),p) : gen_0;
859
860 for (i=1; i<=n; i++) if (r[i])
861 for (j=1; j<=n; j++)
862 gcoeff(a,i,j) = r[j] ? gsub(gcoeff(a,i,j),
863 gmul(gmul(gcoeff(a,k,i),gcoeff(a,k,j)),p))
864 : gen_0;
865 gcoeff(a,k,k) = p;
866 }
867 else
868 {
869 for (k=1; k<=n; k++) if (r[k])
870 {
871 l=k+1; while (l<=n && (!r[l] || gcmp0(gcoeff(a,k,l)))) l++;
872 if (l <= n)
873 {
874 p = gcoeff(a,k,l); r[k] = r[l] = 0;
875 sp++; sn++; t -= 2;
876 for (i=1; i<=n; i++) if (r[i])
877 for (j=1; j<=n; j++)
878 gcoeff(a,i,j) =
879 r[j]? gsub(gcoeff(a,i,j),
880 gdiv(gadd(gmul(gcoeff(a,k,i),gcoeff(a,l,j)),
881 gmul(gcoeff(a,k,j),gcoeff(a,l,i))),
882 p))
883 : gen_0;
884 for (i=1; i<=n; i++) if (r[i])
885 {
886 GEN u = gcoeff(a,k,i);
887 gcoeff(a,k,i) = gdiv(gadd(u, gcoeff(a,l,i)), p);
888 gcoeff(a,l,i) = gdiv(gsub(u, gcoeff(a,l,i)), p);
889 }
890 gcoeff(a,k,l) = gen_1;
891 gcoeff(a,l,k) = gen_m1;
892 gcoeff(a,k,k) = gmul2n(p,-1);
893 gcoeff(a,l,l) = gneg(gcoeff(a,k,k));
894 if (low_stack(lim, stack_lim(av1,1)))
895 {
896 if(DEBUGMEM>1) pari_warn(warnmem,"sqred2");
897 a = gerepilecopy(av1, a);
898 }
899 break;
900 }
901 }
902 if (k > n) break;
903 }
904 }
905 if (!signature) return gerepilecopy(av, a);
906 avma = av; return mkvec2s(sp, sn);
907 }
908
909 GEN
sqred(GEN a)910 sqred(GEN a) { return sqred2(a,0); }
911
912 GEN
signat(GEN a)913 signat(GEN a) { return sqred2(a,1); }
914
915 static void
rot(GEN x,GEN y,GEN s,GEN u)916 rot(GEN x, GEN y, GEN s, GEN u) {
917 GEN x1 = subrr(x,mulrr(s,addrr(y,mulrr(u,x))));
918 GEN y1 = addrr(y,mulrr(s,subrr(x,mulrr(u,y))));
919 affrr(x1,x); affrr(y1,y);
920 }
921
922 /* Diagonalization of a REAL symetric matrix. Return a vector [L, r]:
923 * L = vector of eigenvalues
924 * r = matrix of eigenvectors */
925 GEN
jacobi(GEN a,long prec)926 jacobi(GEN a, long prec)
927 {
928 pari_sp av1;
929 long de, e, e1, e2, i, j, p, q, l = lg(a);
930 GEN c, s, t, u, ja, L, r, unr, x, y;
931
932 if (typ(a) != t_MAT) pari_err(mattype1,"jacobi");
933 ja = cgetg(3,t_VEC);
934 L = cgetg(l,t_COL); gel(ja,1) = L;
935 r = cgetg(l,t_MAT); gel(ja,2) = r;
936 if (l == 1) return ja;
937 if (lg(a[1]) != l) pari_err(mattype1,"jacobi");
938
939 e1 = HIGHEXPOBIT-1;
940 for (j=1; j<l; j++)
941 {
942 gel(L,j) = gtofp(gcoeff(a,j,j), prec);
943 e = expo(L[j]); if (e < e1) e1 = e;
944 }
945 for (j=1; j<l; j++)
946 {
947 gel(r,j) = cgetg(l,t_COL);
948 for (i=1; i<l; i++) gcoeff(r,i,j) = stor(i==j, prec);
949 }
950 av1 = avma;
951
952 e2 = -(long)HIGHEXPOBIT; p = q = 1;
953 c = cgetg(l,t_MAT);
954 for (j=1; j<l; j++)
955 {
956 gel(c,j) = cgetg(j,t_COL);
957 for (i=1; i<j; i++)
958 {
959 gcoeff(c,i,j) = gtofp(gcoeff(a,i,j), prec);
960 e = expo(gcoeff(c,i,j)); if (e > e2) { e2 = e; p = i; q = j; }
961 }
962 }
963 a = c; unr = real_1(prec);
964 de = bit_accuracy(prec);
965
966 /* e1 = min expo(a[i,i])
967 * e2 = max expo(a[i,j]), i != j */
968 while (e1-e2 < de)
969 {
970 pari_sp av2 = avma;
971 /* compute associated rotation in the plane formed by basis vectors number
972 * p and q */
973 x = divrr(subrr(gel(L,q),gel(L,p)), shiftr(gcoeff(a,p,q),1));
974 y = sqrtr(addrr(unr, mulrr(x,x)));
975 t = divrr(unr, (signe(x)>0)? addrr(x,y): subrr(x,y));
976 c = sqrtr(addrr(unr,mulrr(t,t)));
977 s = divrr(t,c);
978 u = divrr(t,addrr(unr,c));
979
980 /* compute successive transforms of a and the matrix of accumulated
981 * rotations (r) */
982 for (i=1; i<p; i++) rot(gcoeff(a,i,p), gcoeff(a,i,q), s,u);
983 for (i=p+1; i<q; i++) rot(gcoeff(a,p,i), gcoeff(a,i,q), s,u);
984 for (i=q+1; i<l; i++) rot(gcoeff(a,p,i), gcoeff(a,q,i), s,u);
985 y = gcoeff(a,p,q);
986 t = mulrr(t, y); setexpo(y, expo(y)-de-1);
987 x = gel(L,p); subrrz(x,t, x);
988 y = gel(L,q); addrrz(y,t, y);
989 for (i=1; i<l; i++) rot(gcoeff(r,i,p), gcoeff(r,i,q), s,u);
990
991 e2 = expo(gcoeff(a,1,2)); p = 1; q = 2;
992 for (j=1; j<l; j++)
993 {
994 for (i=1; i<j; i++)
995 if ((e=expo(gcoeff(a,i,j))) > e2) { e2=e; p=i; q=j; }
996 for (i=j+1; i<l; i++)
997 if ((e=expo(gcoeff(a,j,i))) > e2) { e2=e; p=j; q=i; }
998 }
999 avma = av2;
1000 }
1001 avma = av1; return ja;
1002 }
1003
1004 /*************************************************************************/
1005 /** **/
1006 /** MATRICE RATIONNELLE --> ENTIERE **/
1007 /** **/
1008 /*************************************************************************/
1009
1010 GEN
matrixqz0(GEN x,GEN p)1011 matrixqz0(GEN x,GEN p)
1012 {
1013 if (typ(p)!=t_INT) pari_err(typeer,"matrixqz0");
1014 if (signe(p)>=0) return matrixqz(x,p);
1015 if (equaliu(p,1)) return matrixqz2(x);
1016 if (equaliu(p,2)) return matrixqz3(x);
1017 pari_err(flagerr,"matrixqz"); return NULL; /* not reached */
1018 }
1019
1020 static int
ZV_isin(GEN x)1021 ZV_isin(GEN x)
1022 {
1023 long i,l = lg(x);
1024 for (i=1; i<l; i++)
1025 if (typ(x[i]) != t_INT) return 0;
1026 return 1;
1027 }
1028
1029 GEN
matrixqz(GEN x,GEN p)1030 matrixqz(GEN x, GEN p)
1031 {
1032 pari_sp av = avma, av1, lim;
1033 long i,j,j1,m,n,nfact;
1034 GEN p1,p2,p3;
1035
1036 if (typ(x) != t_MAT) pari_err(typeer,"matrixqz");
1037 n = lg(x)-1; if (!n) return gcopy(x);
1038 m = lg(x[1])-1;
1039 if (n > m) pari_err(talker,"more rows than columns in matrixqz");
1040 if (n==m)
1041 {
1042 p1 = det(x);
1043 if (gcmp0(p1)) pari_err(talker,"matrix of non-maximal rank in matrixqz");
1044 avma = av; return matid(n);
1045 }
1046 /* m > n */
1047 p1 = x; x = cgetg(n+1,t_MAT);
1048 for (j=1; j<=n; j++)
1049 {
1050 gel(x,j) = primpart(gel(p1,j));
1051 if (!ZV_isin(gel(x,j))) pari_err(talker, "not a rational matrix in matrixqz");
1052 }
1053 /* x integral */
1054
1055 if (gcmp0(p))
1056 {
1057 p1 = gtrans(x); setlg(p1,n+1);
1058 p2 = det(p1); p1[n] = p1[n+1]; p2 = ggcd(p2,det(p1));
1059 if (!signe(p2))
1060 pari_err(impl,"matrixqz when the first 2 dets are zero");
1061 if (gcmp1(p2)) return gerepilecopy(av,x);
1062
1063 p1 = (GEN)factor(p2)[1];
1064 }
1065 else p1 = mkvec(p);
1066 nfact = lg(p1)-1;
1067 av1 = avma; lim = stack_lim(av1,1);
1068 for (i=1; i<=nfact; i++)
1069 {
1070 p = gel(p1,i);
1071 for(;;)
1072 {
1073 p2 = FpM_ker(x, p);
1074 if (lg(p2)==1) break;
1075
1076 p2 = centermod(p2,p); p3 = gdiv(gmul(x,p2), p);
1077 for (j=1; j<lg(p2); j++)
1078 {
1079 j1=n; while (gcmp0(gcoeff(p2,j1,j))) j1--;
1080 x[j1] = p3[j];
1081 }
1082 if (low_stack(lim, stack_lim(av1,1)))
1083 {
1084 if (DEBUGMEM>1) pari_warn(warnmem,"matrixqz");
1085 x = gerepilecopy(av1,x);
1086 }
1087 }
1088 }
1089 return gerepilecopy(av,x);
1090 }
1091
1092 static GEN
Z_V_mul(GEN u,GEN A)1093 Z_V_mul(GEN u, GEN A)
1094 {
1095 if (gcmp1(u)) return A;
1096 if (gcmp_1(u)) return gneg(A);
1097 if (gcmp0(u)) return zerocol(lg(A)-1);
1098 return gmul(u,A);
1099 }
1100
1101 static GEN
QV_lincomb(GEN u,GEN v,GEN A,GEN B)1102 QV_lincomb(GEN u, GEN v, GEN A, GEN B)
1103 {
1104 if (!signe(u)) return Z_V_mul(v,B);
1105 if (!signe(v)) return Z_V_mul(u,A);
1106 return gadd(Z_V_mul(u,A), Z_V_mul(v,B));
1107 }
1108
1109 /* cf ZV_elem */
1110 /* zero aj = Aij (!= 0) using ak = Aik (maybe 0), via linear combination of
1111 * A[j] and A[k] of determinant 1. */
1112 static void
QV_elem(GEN aj,GEN ak,GEN A,long j,long k)1113 QV_elem(GEN aj, GEN ak, GEN A, long j, long k)
1114 {
1115 GEN p1,u,v,d, D;
1116
1117 if (gcmp0(ak)) { lswap(A[j],A[k]); return; }
1118 D = lcmii(denom(aj), denom(ak));
1119 if (!is_pm1(D)) { aj = gmul(aj,D); ak = gmul(ak,D); }
1120 d = bezout(aj,ak,&u,&v);
1121 /* frequent special case (u,v) = (1,0) or (0,1) */
1122 if (!signe(u))
1123 { /* ak | aj */
1124 p1 = negi(diviiexact(aj,ak));
1125 gel(A,j) = QV_lincomb(gen_1, p1, gel(A,j), gel(A,k));
1126 return;
1127 }
1128 if (!signe(v))
1129 { /* aj | ak */
1130 p1 = negi(diviiexact(ak,aj));
1131 gel(A,k) = QV_lincomb(gen_1, p1, gel(A,k), gel(A,j));
1132 lswap(A[j], A[k]);
1133 return;
1134 }
1135
1136 if (!is_pm1(d)) { aj = diviiexact(aj,d); ak = diviiexact(ak,d); }
1137 p1 = gel(A,k); aj = negi(aj);
1138 gel(A,k) = QV_lincomb(u,v, gel(A,j),p1);
1139 gel(A,j) = QV_lincomb(aj,ak, p1,gel(A,j));
1140 }
1141
1142 static GEN
matrixqz_aux(GEN A)1143 matrixqz_aux(GEN A)
1144 {
1145 pari_sp av = avma, lim = stack_lim(av,1);
1146 long i,j,k,n,m;
1147 GEN a;
1148
1149 n = lg(A);
1150 if (n == 1) return cgetg(1,t_MAT);
1151 if (n == 2) return hnf(A); /* 1 col, maybe 0 */
1152 m = lg(A[1]);
1153 for (i=1; i<m; i++)
1154 {
1155 for (j = k = 1; j<n; j++)
1156 {
1157 GEN a = gcoeff(A,i,j);
1158 if (gcmp0(a)) continue;
1159
1160 k = j+1; if (k == n) k = 1;
1161 /* zero a = Aij using b = Aik */
1162 QV_elem(a, gcoeff(A,i,k), A, j,k);
1163 }
1164 a = gcoeff(A,i,k);
1165 if (!gcmp0(a))
1166 {
1167 a = denom(a);
1168 if (!is_pm1(a)) gel(A,k) = gmul(gel(A,k), a);
1169 }
1170 if (low_stack(lim, stack_lim(av,1)))
1171 {
1172 if(DEBUGMEM>1) pari_warn(warnmem,"matrixqz_aux");
1173 A = gerepilecopy(av,A);
1174 }
1175 }
1176 return m > 100? hnfall_i(A,NULL,1): hnf(A);
1177 }
1178
1179 GEN
matrixqz2(GEN x)1180 matrixqz2(GEN x)
1181 {
1182 pari_sp av = avma;
1183 if (typ(x)!=t_MAT) pari_err(typeer,"matrixqz2");
1184 x = shallowcopy(x);
1185 return gerepileupto(av, matrixqz_aux(x));
1186 }
1187
1188 GEN
matrixqz3(GEN x)1189 matrixqz3(GEN x)
1190 {
1191 pari_sp av = avma, av1, lim;
1192 long j,j1,k,m,n;
1193 GEN c;
1194
1195 if (typ(x)!=t_MAT) pari_err(typeer,"matrixqz3");
1196 n = lg(x); if (n==1) return gcopy(x);
1197 m = lg(x[1]); x = shallowcopy(x);
1198 c = cgetg(n,t_VECSMALL);
1199 for (j=1; j<n; j++) c[j] = 0;
1200 av1 = avma; lim = stack_lim(av1,1);
1201 for (k=1; k<m; k++)
1202 {
1203 j=1; while (j<n && (c[j] || gcmp0(gcoeff(x,k,j)))) j++;
1204 if (j==n) continue;
1205
1206 c[j]=k; gel(x,j) = gdiv(gel(x,j),gcoeff(x,k,j));
1207 for (j1=1; j1<n; j1++)
1208 if (j1!=j)
1209 {
1210 GEN t = gcoeff(x,k,j1);
1211 if (!gcmp0(t)) gel(x,j1) = gsub(gel(x,j1),gmul(t,gel(x,j)));
1212 }
1213 if (low_stack(lim, stack_lim(av1,1)))
1214 {
1215 if(DEBUGMEM>1) pari_warn(warnmem,"matrixqz3");
1216 x = gerepilecopy(av1,x);
1217 }
1218 }
1219 return gerepileupto(av, matrixqz_aux(x));
1220 }
1221
1222 GEN
intersect(GEN x,GEN y)1223 intersect(GEN x, GEN y)
1224 {
1225 pari_sp av,tetpil;
1226 long j, lx = lg(x);
1227 GEN z;
1228
1229 if (typ(x)!=t_MAT || typ(y)!=t_MAT) pari_err(typeer,"intersect");
1230 if (lx==1 || lg(y)==1) return cgetg(1,t_MAT);
1231
1232 av=avma; z=ker(shallowconcat(x,y));
1233 for (j=lg(z)-1; j; j--) setlg(z[j],lx);
1234 tetpil=avma; return gerepile(av,tetpil,gmul(x,z));
1235 }
1236
1237 /**************************************************************/
1238 /** **/
1239 /** HERMITE NORMAL FORM REDUCTION **/
1240 /** **/
1241 /**************************************************************/
1242 /* negate in place, except universal constants */
1243 static GEN
mynegi(GEN x)1244 mynegi(GEN x)
1245 {
1246 long s = signe(x);
1247
1248 if (!s) return gen_0;
1249 if (is_pm1(x)) return (s > 0)? gen_m1: gen_1;
1250 setsigne(x,-s); return x;
1251 }
1252 void
ZV_neg_ip(GEN M)1253 ZV_neg_ip(GEN M)
1254 {
1255 long i;
1256 for (i = lg(M)-1; i; i--)
1257 gel(M,i) = mynegi(gel(M,i));
1258 }
1259
1260 GEN
mathnf0(GEN x,long flag)1261 mathnf0(GEN x, long flag)
1262 {
1263 switch(flag)
1264 {
1265 case 0: return hnf(x);
1266 case 1: return hnfall(x);
1267 case 3: return hnfperm(x);
1268 case 4: return hnflll(x);
1269 default: pari_err(flagerr,"mathnf");
1270 }
1271 return NULL; /* not reached */
1272 }
1273
1274 static GEN
init_hnf(GEN x,GEN * denx,long * co,long * li,pari_sp * av)1275 init_hnf(GEN x, GEN *denx, long *co, long *li, pari_sp *av)
1276 {
1277 if (typ(x) != t_MAT) pari_err(typeer,"mathnf");
1278 *co=lg(x); if (*co==1) return NULL; /* empty matrix */
1279 *li=lg(x[1]); *denx=denom(x); *av=avma;
1280
1281 if (gcmp1(*denx)) /* no denominator */
1282 { *denx = NULL; return shallowcopy(x); }
1283 return Q_muli_to_int(x, *denx);
1284 }
1285
1286 /* check whether x is upper trinagular with positive diagonal coeffs */
1287 int
RgM_ishnf(GEN x)1288 RgM_ishnf(GEN x)
1289 {
1290 long i,j, lx = lg(x);
1291 for (i=2; i<lx; i++)
1292 {
1293 if (gsigne(gcoeff(x,i,i)) <= 0) return 0;
1294 for (j=1; j<i; j++)
1295 if (!gcmp0(gcoeff(x,i,j))) return 0;
1296 }
1297 return (gsigne(gcoeff(x,1,1)) > 0);
1298 }
1299 /* same x is known to be integral */
1300 int
ZM_ishnf(GEN x)1301 ZM_ishnf(GEN x)
1302 {
1303 long i,j, lx = lg(x);
1304 for (i=2; i<lx; i++)
1305 {
1306 if (signe(gcoeff(x,i,i)) <= 0) return 0;
1307 for (j=1; j<i; j++)
1308 if (signe(gcoeff(x,i,j))) return 0;
1309 }
1310 return (signe(gcoeff(x,1,1)) > 0);
1311 }
1312
1313 GEN
ZV_to_nv(GEN z)1314 ZV_to_nv(GEN z)
1315 {
1316 long i, l = lg(z);
1317 GEN x = cgetg(l, t_VECSMALL);
1318 for (i=1; i<l; i++) x[i] = itou(gel(z,i));
1319 return x;
1320 }
1321
1322 GEN
ZV_copy(GEN x)1323 ZV_copy(GEN x)
1324 {
1325 long i, lx = lg(x);
1326 GEN y = cgetg(lx, t_COL);
1327 for (i=1; i<lx; i++) gel(y,i) = signe(x[i])? icopy(gel(x,i)): gen_0;
1328 return y;
1329 }
1330
1331 GEN
ZM_copy(GEN x)1332 ZM_copy(GEN x)
1333 {
1334 long i, lx = lg(x);
1335 GEN y = cgetg(lx, t_MAT);
1336
1337 for (i=1; i<lx; i++)
1338 gel(y,i) = ZV_copy(gel(x,i));
1339 return y;
1340 }
1341
1342 /* return c * X. Not memory clean if c = 1 */
1343 GEN
ZV_Z_mul(GEN X,GEN c)1344 ZV_Z_mul(GEN X, GEN c)
1345 {
1346 long i,m = lg(X);
1347 GEN A = new_chunk(m);
1348 if (signe(c) && is_pm1(c))
1349 {
1350 if (signe(c) > 0)
1351 { for (i=1; i<m; i++) A[i] = X[i]; }
1352 else
1353 { for (i=1; i<m; i++) gel(A,i) = negi(gel(X,i)); }
1354 }
1355 else /* c = 0 should be exceedingly rare */
1356 { for (i=1; i<m; i++) gel(A,i) = mulii(c,gel(X,i)); }
1357 A[0] = X[0]; return A;
1358 }
1359
1360 /* X + v Y */
1361 static GEN
ZV_lincomb1(GEN v,GEN X,GEN Y)1362 ZV_lincomb1(GEN v, GEN X, GEN Y)
1363 {
1364 long i, lx = lg(X);
1365 GEN p1, p2, A = cgetg(lx,t_COL);
1366 if (is_bigint(v))
1367 {
1368 long m = lgefint(v);
1369 for (i=1; i<lx; i++)
1370 {
1371 p1=gel(X,i); p2=gel(Y,i);
1372 if (!signe(p1)) gel(A,i) = mulii(v,p2);
1373 else if (!signe(p2)) gel(A,i) = icopy(p1);
1374 else
1375 {
1376 pari_sp av = avma; (void)new_chunk(m+lgefint(p1)+lgefint(p2)); /*HACK*/
1377 p2 = mulii(v,p2);
1378 avma = av; gel(A,i) = addii(p1,p2);
1379 }
1380 }
1381 }
1382 else
1383 {
1384 long w = itos(v);
1385 for (i=1; i<lx; i++)
1386 {
1387 p1=gel(X,i); p2=gel(Y,i);
1388 if (!signe(p1)) gel(A,i) = mulsi(w,p2);
1389 else if (!signe(p2)) gel(A,i) = icopy(p1);
1390 else
1391 {
1392 pari_sp av = avma; (void)new_chunk(1+lgefint(p1)+lgefint(p2)); /*HACK*/
1393 p2 = mulsi(w,p2);
1394 avma = av; gel(A,i) = addii(p1,p2);
1395 }
1396 }
1397 }
1398 return A;
1399 }
1400 /* -X + vY */
1401 static GEN
ZV_lincomb_1(GEN v,GEN X,GEN Y)1402 ZV_lincomb_1(GEN v, GEN X, GEN Y)
1403 {
1404 long i, m, lx = lg(X);
1405 GEN p1, p2, A = cgetg(lx,t_COL);
1406 m = lgefint(v);
1407 for (i=1; i<lx; i++)
1408 {
1409 p1=gel(X,i); p2=gel(Y,i);
1410 if (!signe(p1)) gel(A,i) = mulii(v,p2);
1411 else if (!signe(p2)) gel(A,i) = negi(p1);
1412 else
1413 {
1414 pari_sp av = avma; (void)new_chunk(m+lgefint(p1)+lgefint(p2)); /* HACK */
1415 p2 = mulii(v,p2);
1416 avma = av; gel(A,i) = subii(p2,p1);
1417 }
1418 }
1419 return A;
1420 }
1421 GEN
ZV_add(GEN x,GEN y)1422 ZV_add(GEN x, GEN y)
1423 {
1424 long i, lx = lg(x);
1425 GEN A = cgetg(lx, t_COL);
1426 for (i=1; i<lx; i++) gel(A,i) = addii(gel(x,i), gel(y,i));
1427 return A;
1428 }
1429 GEN
ZV_sub(GEN x,GEN y)1430 ZV_sub(GEN x, GEN y)
1431 {
1432 long i, lx = lg(x);
1433 GEN A = cgetg(lx, t_COL);
1434 for (i=1; i<lx; i++) gel(A,i) = subii(gel(x,i), gel(y,i));
1435 return A;
1436 }
1437 /* X,Y columns; u,v scalars; everybody is integral. return A = u*X + v*Y
1438 * Not memory clean if (u,v) = (1,0) or (0,1) */
1439 GEN
ZV_lincomb(GEN u,GEN v,GEN X,GEN Y)1440 ZV_lincomb(GEN u, GEN v, GEN X, GEN Y)
1441 {
1442 pari_sp av;
1443 long i, lx, m, su, sv;
1444 GEN p1, p2, A;
1445
1446 su = signe(u); if (!su) return ZV_Z_mul(Y, v);
1447 sv = signe(v); if (!sv) return ZV_Z_mul(X, u);
1448 if (is_pm1(v))
1449 {
1450 if (is_pm1(u))
1451 {
1452 if (su != sv) A = ZV_sub(X, Y);
1453 else A = ZV_add(X, Y);
1454 if (su < 0) ZV_neg_ip(A);
1455 }
1456 else
1457 {
1458 if (sv > 0) A = ZV_lincomb1 (u, Y, X);
1459 else A = ZV_lincomb_1(u, Y, X);
1460 }
1461 }
1462 else if (is_pm1(u))
1463 {
1464 if (su > 0) A = ZV_lincomb1 (v, X, Y);
1465 else A = ZV_lincomb_1(v, X, Y);
1466 }
1467 else
1468 {
1469 lx = lg(X); A = cgetg(lx,t_COL); m = lgefint(u)+lgefint(v);
1470 for (i=1; i<lx; i++)
1471 {
1472 p1=gel(X,i); p2=gel(Y,i);
1473 if (!signe(p1)) gel(A,i) = mulii(v,p2);
1474 else if (!signe(p2)) gel(A,i) = mulii(u,p1);
1475 else
1476 {
1477 av = avma; (void)new_chunk(m+lgefint(p1)+lgefint(p2)); /* HACK */
1478 p1 = mulii(u,p1);
1479 p2 = mulii(v,p2);
1480 avma = av; gel(A,i) = addii(p1,p2);
1481 }
1482 }
1483 }
1484 return A;
1485 }
1486
1487 /* x = [A,U], nbcol(A) = nbcol(U), A integral. Return [AV, UV], with AV HNF */
1488 GEN
hnf_special(GEN x,long remove)1489 hnf_special(GEN x, long remove)
1490 {
1491 pari_sp av0,av,tetpil,lim;
1492 long s,li,co,i,j,k,def,ldef;
1493 GEN p1,u,v,d,denx,a,b, x2,res;
1494
1495 if (typ(x) != t_VEC || lg(x) !=3) pari_err(typeer,"hnf_special");
1496 res = cgetg(3,t_VEC);
1497 av0 = avma;
1498 x2 = gel(x,2);
1499 x = gel(x,1);
1500 x = init_hnf(x,&denx,&co,&li,&av);
1501 if (!x) return cgetg(1,t_MAT);
1502
1503 lim = stack_lim(av,1);
1504 def=co-1; ldef=(li>co)? li-co: 0;
1505 if (lg(x2) != co) pari_err(talker,"incompatible matrices in hnf_special");
1506 x2 = shallowcopy(x2);
1507 for (i=li-1; i>ldef; i--)
1508 {
1509 for (j=def-1; j; j--)
1510 {
1511 a = gcoeff(x,i,j);
1512 if (!signe(a)) continue;
1513
1514 k = (j==1)? def: j-1;
1515 b = gcoeff(x,i,k); d = bezout(a,b,&u,&v);
1516 if (!is_pm1(d)) { a = diviiexact(a,d); b = diviiexact(b,d); }
1517 p1 = gel(x,j); b = negi(b);
1518 gel(x,j) = ZV_lincomb(a,b, gel(x,k), p1);
1519 gel(x,k) = ZV_lincomb(u,v, p1, gel(x,k));
1520 p1 = gel(x2,j);
1521 gel(x2,j) = gadd(gmul(a, gel(x2,k)), gmul(b,p1));
1522 gel(x2,k) = gadd(gmul(u,p1), gmul(v, gel(x2,k)));
1523 if (low_stack(lim, stack_lim(av,1)))
1524 {
1525 GEN *gptr[2]; gptr[0]=&x; gptr[1]=&x2;
1526 if (DEBUGMEM>1) pari_warn(warnmem,"hnf_special[1]. i=%ld",i);
1527 gerepilemany(av,gptr,2);
1528 }
1529 }
1530 p1 = gcoeff(x,i,def); s = signe(p1);
1531 if (s)
1532 {
1533 if (s < 0)
1534 {
1535 gel(x,def) = gneg(gel(x,def)); p1 = gcoeff(x,i,def);
1536 gel(x2,def)= gneg(gel(x2,def));
1537 }
1538 for (j=def+1; j<co; j++)
1539 {
1540 b = negi(gdivent(gcoeff(x,i,j),p1));
1541 gel(x,j) = ZV_lincomb(gen_1,b, gel(x,j), gel(x,def));
1542 gel(x2,j) = gadd(gel(x2,j), gmul(b, gel(x2,def)));
1543 }
1544 def--;
1545 }
1546 else
1547 if (ldef && i==ldef+1) ldef--;
1548 if (low_stack(lim, stack_lim(av,1)))
1549 {
1550 GEN *gptr[2]; gptr[0]=&x; gptr[1]=&x2;
1551 if (DEBUGMEM>1) pari_warn(warnmem,"hnf_special[2]. i=%ld",i);
1552 gerepilemany(av,gptr,2);
1553 }
1554 }
1555 if (remove)
1556 { /* remove null columns */
1557 for (i=1,j=1; j<co; j++)
1558 if (!gcmp0(gel(x,j)))
1559 {
1560 x[i] = x[j];
1561 x2[i] = x2[j]; i++;
1562 }
1563 setlg(x,i);
1564 setlg(x2,i);
1565 }
1566 tetpil=avma;
1567 x = denx? gdiv(x,denx): ZM_copy(x);
1568 x2 = gcopy(x2);
1569 {
1570 GEN *gptr[2]; gptr[0]=&x; gptr[1]=&x2;
1571 gerepilemanysp(av0,tetpil,gptr,2);
1572 }
1573 gel(res,1) = x;
1574 gel(res,2) = x2;
1575 return res;
1576 }
1577
1578 /*******************************************************************/
1579 /* */
1580 /* SPECIAL HNF (FOR INTERNAL USE !!!) */
1581 /* */
1582 /*******************************************************************/
1583 static int
count(long ** mat,long row,long len,long * firstnonzero)1584 count(long **mat, long row, long len, long *firstnonzero)
1585 {
1586 long j, n = 0;
1587
1588 for (j=1; j<=len; j++)
1589 {
1590 long p = mat[j][row];
1591 if (p)
1592 {
1593 if (labs(p)!=1) return -1;
1594 n++; *firstnonzero=j;
1595 }
1596 }
1597 return n;
1598 }
1599
1600 static int
count2(long ** mat,long row,long len)1601 count2(long **mat, long row, long len)
1602 {
1603 long j;
1604 for (j=len; j; j--)
1605 if (labs(mat[j][row]) == 1) return j;
1606 return 0;
1607 }
1608
1609 static GEN
hnffinal(GEN matgen,GEN perm,GEN * ptdep,GEN * ptB,GEN * ptC)1610 hnffinal(GEN matgen,GEN perm,GEN* ptdep,GEN* ptB,GEN* ptC)
1611 {
1612 GEN p1,p2,U,H,Hnew,Bnew,Cnew,diagH1;
1613 GEN B = *ptB, C = *ptC, dep = *ptdep, depnew;
1614 pari_sp av, lim;
1615 long i,j,k,s,i1,j1,zc;
1616 long co = lg(C);
1617 long col = lg(matgen)-1;
1618 long lnz, nlze, lig;
1619
1620 if (col == 0) return matgen;
1621 lnz = lg(matgen[1])-1; /* was called lnz-1 - nr in hnfspec */
1622 nlze = lg(dep[1])-1;
1623 lig = nlze + lnz;
1624 if (DEBUGLEVEL>5)
1625 {
1626 fprintferr("Entering hnffinal:\n");
1627 if (DEBUGLEVEL>6)
1628 {
1629 if (nlze) fprintferr("dep = %Z\n",dep);
1630 fprintferr("mit = %Z\n",matgen);
1631 fprintferr("B = %Z\n",B);
1632 }
1633 }
1634 /* H: lnz x lnz [disregarding initial 0 cols], U: col x col */
1635 H = hnflll_i(matgen, &U, 0);
1636 H += (lg(H)-1 - lnz); H[0] = evaltyp(t_MAT) | evallg(lnz+1);
1637 /* Only keep the part above the H (above the 0s is 0 since the dep rows
1638 * are dependent from the ones in matgen) */
1639 zc = col - lnz; /* # of 0 columns, correspond to units */
1640 if (nlze) { dep = gmul(dep,U); dep += zc; }
1641
1642 diagH1 = new_chunk(lnz+1); /* diagH1[i] = 0 iff H[i,i] != 1 (set later) */
1643
1644 av = avma; lim = stack_lim(av,1);
1645 Cnew = cgetg(co, typ(C));
1646 setlg(C, col+1); p1 = gmul(C,U);
1647 for (j=1; j<=col; j++) Cnew[j] = p1[j];
1648 for ( ; j<co ; j++) Cnew[j] = C[j];
1649 if (DEBUGLEVEL>5) fprintferr(" hnflll done\n");
1650
1651 /* Clean up B using new H */
1652 for (s=0,i=lnz; i; i--)
1653 {
1654 GEN Di = gel(dep,i), Hi = gel(H,i);
1655 GEN h = gel(Hi,i); /* H[i,i] */
1656 if ( (diagH1[i] = gcmp1(h)) ) { h = NULL; s++; }
1657 for (j=col+1; j<co; j++)
1658 {
1659 GEN z = gel(B,j-col);
1660 p1 = gel(z,i+nlze); if (!signe(p1)) continue;
1661
1662 if (h) p1 = gdivent(p1,h);
1663 for (k=1; k<=nlze; k++) gel(z,k) = subii(gel(z,k), mulii(p1, gel(Di,k)));
1664 for ( ; k<=lig; k++) gel(z,k) = subii(gel(z,k), mulii(p1, gel(Hi,k-nlze)));
1665 gel(Cnew,j) = gsub(gel(Cnew,j), gmul(p1, gel(Cnew,i+zc)));
1666 }
1667 if (low_stack(lim, stack_lim(av,1)))
1668 {
1669 if(DEBUGMEM>1) pari_warn(warnmem,"hnffinal, i = %ld",i);
1670 gerepileall(av, 2, &Cnew, &B);
1671 }
1672 }
1673 p1 = cgetg(lnz+1,t_VEC); p2 = perm + nlze;
1674 for (i1=0, j1=lnz-s, i=1; i<=lnz; i++) /* push the 1 rows down */
1675 if (diagH1[i])
1676 p1[++j1] = p2[i];
1677 else
1678 p2[++i1] = p2[i];
1679 for (i=i1+1; i<=lnz; i++) p2[i] = p1[i];
1680 if (DEBUGLEVEL>5) fprintferr(" first pass in hnffinal done\n");
1681
1682 /* s = # extra redundant generators taken from H
1683 * zc col-s co zc = col - lnz
1684 * [ 0 |dep | ] i = lnze + lnz - s = lig - s
1685 * nlze [--------| B' ]
1686 * [ 0 | H' | ] H' = H minus the s rows with a 1 on diagonal
1687 * i [--------|-----] lig-s (= "1-rows")
1688 * [ 0 | Id ]
1689 * [ | ] li */
1690 lig -= s; col -= s; lnz -= s;
1691 Hnew = cgetg(lnz+1,t_MAT);
1692 depnew = cgetg(lnz+1,t_MAT); /* only used if nlze > 0 */
1693 Bnew = cgetg(co-col,t_MAT);
1694 C = shallowcopy(Cnew);
1695 for (j=1,i1=j1=0; j<=lnz+s; j++)
1696 {
1697 GEN z = gel(H,j);
1698 if (diagH1[j])
1699 { /* hit exactly s times */
1700 i1++; C[i1+col] = Cnew[j+zc];
1701 p1 = cgetg(lig+1,t_COL); gel(Bnew,i1) = p1;
1702 for (i=1; i<=nlze; i++) gel(p1,i) = gcoeff(dep,i,j);
1703 p1 += nlze;
1704 }
1705 else
1706 {
1707 j1++; C[j1+zc] = Cnew[j+zc];
1708 p1 = cgetg(lnz+1,t_COL); gel(Hnew,j1) = p1;
1709 depnew[j1] = dep[j];
1710 }
1711 for (i=k=1; k<=lnz; i++)
1712 if (!diagH1[i]) p1[k++] = z[i];
1713 }
1714 for (j=s+1; j<co-col; j++)
1715 {
1716 GEN z = gel(B,j-s);
1717 p1 = cgetg(lig+1,t_COL); gel(Bnew,j) = p1;
1718 for (i=1; i<=nlze; i++) p1[i] = z[i];
1719 z += nlze; p1 += nlze;
1720 for (i=k=1; k<=lnz; i++)
1721 if (!diagH1[i]) p1[k++] = z[i];
1722 }
1723 if (DEBUGLEVEL>5)
1724 {
1725 fprintferr("Leaving hnffinal\n");
1726 if (DEBUGLEVEL>6)
1727 {
1728 if (nlze) fprintferr("dep = %Z\n",depnew);
1729 fprintferr("mit = %Z\nB = %Z\nC = %Z\n", Hnew, Bnew, C);
1730 }
1731 }
1732 *ptdep = depnew;
1733 *ptC = C;
1734 *ptB = Bnew; return Hnew;
1735 }
1736
1737 /* for debugging */
1738 static void
p_mat(long ** mat,GEN perm,long k)1739 p_mat(long **mat, GEN perm, long k)
1740 {
1741 pari_sp av = avma;
1742 perm = vecslice(perm, k+1, lg(perm)-1);
1743 fprintferr("Permutation: %Z\n",perm);
1744 if (DEBUGLEVEL > 6)
1745 fprintferr("matgen = %Z\n", zm_to_ZM( rowpermute((GEN)mat, perm) ));
1746 avma = av;
1747 }
1748
1749 /* M_k <- M_k + q M_i (col operations) */
1750 static void
elt_col(GEN Mk,GEN Mi,GEN q)1751 elt_col(GEN Mk, GEN Mi, GEN q)
1752 {
1753 long j;
1754 if (is_pm1(q))
1755 {
1756 if (signe(q) > 0)
1757 {
1758 for (j = lg(Mk)-1; j; j--)
1759 if (signe(Mi[j])) gel(Mk,j) = addii(gel(Mk,j), gel(Mi,j));
1760 }
1761 else
1762 {
1763 for (j = lg(Mk)-1; j; j--)
1764 if (signe(Mi[j])) gel(Mk,j) = subii(gel(Mk,j), gel(Mi,j));
1765 }
1766 }
1767 else
1768 for (j = lg(Mk)-1; j; j--)
1769 if (signe(Mi[j])) gel(Mk,j) = addii(gel(Mk,j), mulii(q, gel(Mi,j)));
1770 }
1771
1772 static GEN
col_dup(long l,GEN col)1773 col_dup(long l, GEN col)
1774 {
1775 GEN c = new_chunk(l);
1776 memcpy(c,col,l * sizeof(long)); return c;
1777 }
1778
1779 /* HNF reduce a relation matrix (column operations + row permutation)
1780 ** Input:
1781 ** mat = (li-1) x (co-1) matrix of long
1782 ** C = r x (co-1) matrix of GEN
1783 ** perm= permutation vector (length li-1), indexing the rows of mat: easier
1784 ** to maintain perm than to copy rows. For columns we can do it directly
1785 ** using e.g. lswap(mat[i], mat[j])
1786 ** k0 = integer. The k0 first lines of mat are dense, the others are sparse.
1787 ** Output: cf ASCII art in the function body
1788 **
1789 ** row permutations applied to perm
1790 ** column operations applied to C. IN PLACE
1791 **/
1792 GEN
hnfspec_i(long ** mat0,GEN perm,GEN * ptdep,GEN * ptB,GEN * ptC,long k0)1793 hnfspec_i(long** mat0, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, long k0)
1794 {
1795 pari_sp av, lim;
1796 long co, n, s, nlze, lnz, nr, i, j, k, lk0, col, lig, *p, *matj, **mat;
1797 GEN p1, p2, matb, matbnew, vmax, matt, T, extramat, B, C, H, dep, permpro;
1798 const long li = lg(perm); /* = lg(mat0[1]) */
1799 const long CO = lg(mat0);
1800
1801 n = 0; /* -Wall */
1802
1803 C = *ptC; co = CO;
1804 if (co > 300 && co > 1.5 * li)
1805 { /* treat the rest at the end */
1806 co = (long)(1.2 * li);
1807 setlg(C, co);
1808 }
1809
1810 if (DEBUGLEVEL>5)
1811 {
1812 fprintferr("Entering hnfspec\n");
1813 p_mat(mat0,perm,0);
1814 }
1815 matt = cgetg(co, t_MAT); /* dense part of mat (top) */
1816 mat = (long**)cgetg(co, t_MAT);
1817 for (j = 1; j < co; j++)
1818 {
1819 mat[j] = col_dup(li, gel(mat0,j));
1820 p1 = cgetg(k0+1,t_COL); gel(matt,j) = p1; matj = mat[j];
1821 for (i=1; i<=k0; i++) gel(p1,i) = stoi(matj[perm[i]]);
1822 }
1823 vmax = cgetg(co,t_VECSMALL);
1824 av = avma; lim = stack_lim(av,1);
1825
1826 i = lig = li-1; col = co-1; lk0 = k0;
1827 T = (k0 || (lg(C) > 1 && lg(C[1]) > 1))? matid(col): NULL;
1828 /* Look for lines with a single non-0 entry, equal to 1 in absolute value */
1829 while (i > lk0)
1830 switch( count(mat,perm[i],col,&n) )
1831 {
1832 case 0: /* move zero lines between k0+1 and lk0 */
1833 lk0++; lswap(perm[i], perm[lk0]);
1834 i = lig; continue;
1835
1836 case 1: /* move trivial generator between lig+1 and li */
1837 lswap(perm[i], perm[lig]);
1838 if (T) lswap(T[n], T[col]);
1839 swap(mat[n], mat[col]); p = mat[col];
1840 if (p[perm[lig]] < 0) /* = -1 */
1841 { /* convert relation -g = 0 to g = 0 */
1842 for (i=lk0+1; i<lig; i++) p[perm[i]] = -p[perm[i]];
1843 if (T)
1844 {
1845 p1 = gel(T,col);
1846 for (i=1; ; i++) /* T is a permuted identity: single non-0 entry */
1847 if (signe(gel(p1,i))) { gel(p1,i) = negi(gel(p1,i)); break; }
1848 }
1849 }
1850 lig--; col--; i = lig; continue;
1851
1852 default: i--;
1853 }
1854 if (DEBUGLEVEL>5)
1855 {
1856 fprintferr(" after phase1:\n");
1857 p_mat(mat,perm,0);
1858 }
1859
1860 #define absmax(s,z) {long _z; _z = labs(z); if (_z > s) s = _z;}
1861 /* Get rid of all lines containing only 0 and +/- 1, keeping track of column
1862 * operations in T. Leave the rows 1..lk0 alone [up to k0, coefficient
1863 * explosion, between k0+1 and lk0, row is 0] */
1864 s = 0;
1865 while (lig > lk0 && s < (long)(HIGHBIT>>1))
1866 {
1867 for (i=lig; i>lk0; i--)
1868 if (count(mat,perm[i],col,&n) > 0) break;
1869 if (i == lk0) break;
1870
1871 /* only 0, +/- 1 entries, at least 2 of them non-zero */
1872 lswap(perm[i], perm[lig]);
1873 swap(mat[n], mat[col]); p = mat[col];
1874 if (T) lswap(T[n], T[col]);
1875 if (p[perm[lig]] < 0)
1876 {
1877 for (i=lk0+1; i<=lig; i++) p[perm[i]] = -p[perm[i]];
1878 if (T) ZV_neg_ip(gel(T,col));
1879 }
1880 for (j=1; j<col; j++)
1881 {
1882 long t;
1883 matj = mat[j];
1884 if (! (t = matj[perm[lig]]) ) continue;
1885 if (t == 1) {
1886 for (i=lk0+1; i<=lig; i++) absmax(s, matj[perm[i]] -= p[perm[i]]);
1887 }
1888 else { /* t = -1 */
1889 for (i=lk0+1; i<=lig; i++) absmax(s, matj[perm[i]] += p[perm[i]]);
1890 }
1891 if (T) elt_col(gel(T,j), gel(T,col), stoi(-t));
1892 }
1893 lig--; col--;
1894 if (low_stack(lim, stack_lim(av,1)))
1895 {
1896 if(DEBUGMEM>1) pari_warn(warnmem,"hnfspec[1]");
1897 if (T) T = gerepilecopy(av, T); else avma = av;
1898 }
1899 }
1900 /* As above with lines containing a +/- 1 (no other assumption).
1901 * Stop when single precision becomes dangerous */
1902 for (j=1; j<=col; j++)
1903 {
1904 matj = mat[j];
1905 for (s=0, i=lk0+1; i<=lig; i++) absmax(s, matj[i]);
1906 vmax[j] = s;
1907 }
1908 while (lig > lk0)
1909 {
1910 for (i=lig; i>lk0; i--)
1911 if ( (n = count2(mat,perm[i],col)) ) break;
1912 if (i == lk0) break;
1913
1914 lswap(vmax[n], vmax[col]);
1915 lswap(perm[i], perm[lig]);
1916 swap(mat[n], mat[col]); p = mat[col];
1917 if (T) lswap(T[n], T[col]);
1918 if (p[perm[lig]] < 0)
1919 {
1920 for (i=lk0+1; i<=lig; i++) p[perm[i]] = -p[perm[i]];
1921 if (T) ZV_neg_ip(gel(T,col));
1922 }
1923 for (j=1; j<col; j++)
1924 {
1925 long t;
1926 matj = mat[j];
1927 if (! (t = matj[perm[lig]]) ) continue;
1928 if (vmax[col] && (ulong)labs(t) >= (HIGHBIT-vmax[j]) / vmax[col])
1929 goto END2;
1930
1931 for (s=0, i=lk0+1; i<=lig; i++) absmax(s, matj[perm[i]] -= t*p[perm[i]]);
1932 vmax[j] = s;
1933 if (T) elt_col(gel(T,j), gel(T,col), stoi(-t));
1934 }
1935 lig--; col--;
1936 if (low_stack(lim, stack_lim(av,1)))
1937 {
1938 if(DEBUGMEM>1) pari_warn(warnmem,"hnfspec[2]");
1939 if (T) T = gerepilecopy(av,T); else avma = av;
1940 }
1941 }
1942
1943 END2: /* clean up mat: remove everything to the right of the 1s on diagonal */
1944 /* go multiprecision first */
1945 matb = cgetg(co,t_MAT); /* bottom part (complement of matt) */
1946 for (j=1; j<co; j++)
1947 {
1948 p1 = cgetg(li-k0,t_COL); gel(matb,j) = p1;
1949 p1 -= k0; matj = mat[j];
1950 for (i=k0+1; i<li; i++) gel(p1,i) = stoi(matj[perm[i]]);
1951 }
1952 if (DEBUGLEVEL>5)
1953 {
1954 fprintferr(" after phase2:\n");
1955 p_mat(mat,perm,lk0);
1956 }
1957 for (i=li-2; i>lig; i--)
1958 {
1959 long h, i0 = i - k0, k = i + co-li;
1960 GEN Bk = gel(matb,k);
1961 for (j=k+1; j<co; j++)
1962 {
1963 GEN Bj = gel(matb,j), v = gel(Bj,i0);
1964 s = signe(v); if (!s) continue;
1965
1966 gel(Bj,i0) = gen_0;
1967 if (is_pm1(v))
1968 {
1969 if (s > 0) /* v = 1 */
1970 { for (h=1; h<i0; h++) gel(Bj,h) = subii(gel(Bj,h), gel(Bk,h)); }
1971 else /* v = -1 */
1972 { for (h=1; h<i0; h++) gel(Bj,h) = addii(gel(Bj,h), gel(Bk,h)); }
1973 }
1974 else {
1975 for (h=1; h<i0; h++) gel(Bj,h) = subii(gel(Bj,h), mulii(v,gel(Bk,h)));
1976 }
1977 if (T) elt_col(gel(T,j), gel(T,k), negi(v));
1978 if (low_stack(lim, stack_lim(av,1)))
1979 {
1980 if(DEBUGMEM>1) pari_warn(warnmem,"hnfspec[3], (i,j) = %ld,%ld", i,j);
1981 for (h=1; h<co; h++) setlg(matb[h], i0+1); /* bottom can be forgotten */
1982 gerepileall(av, T? 2: 1, &matb, &T);
1983 Bk = gel(matb,k);
1984 }
1985 }
1986 }
1987 for (j=1; j<co; j++) setlg(matb[j], lig-k0+1); /* bottom can be forgotten */
1988 gerepileall(av, T? 2: 1, &matb, &T);
1989 if (DEBUGLEVEL>5) fprintferr(" matb cleaned up (using Id block)\n");
1990
1991 nlze = lk0 - k0; /* # of 0 rows */
1992 lnz = lig-nlze+1; /* 1 + # of non-0 rows (!= 0...0 1 0 ... 0) */
1993 if (T) matt = gmul(matt,T); /* update top rows */
1994 extramat = cgetg(col+1,t_MAT); /* = new C minus the 0 rows */
1995 for (j=1; j<=col; j++)
1996 {
1997 GEN z = gel(matt,j);
1998 GEN t = (gel(matb,j)) + nlze - k0;
1999 p2=cgetg(lnz,t_COL); gel(extramat,j) = p2;
2000 for (i=1; i<=k0; i++) p2[i] = z[i]; /* top k0 rows */
2001 for ( ; i<lnz; i++) p2[i] = t[i]; /* other non-0 rows */
2002 }
2003 permpro = imagecomplspec(extramat, &nr); /* lnz = lg(permpro) */
2004
2005 if (nlze)
2006 { /* put the nlze 0 rows (trivial generators) at the top */
2007 p1 = new_chunk(lk0+1);
2008 for (i=1; i<=nlze; i++) p1[i] = perm[i + k0];
2009 for ( ; i<=lk0; i++) p1[i] = perm[i - nlze];
2010 for (i=1; i<=lk0; i++) perm[i] = p1[i];
2011 }
2012 /* sort other rows according to permpro (nr redundant generators first) */
2013 p1 = new_chunk(lnz); p2 = perm + nlze;
2014 for (i=1; i<lnz; i++) p1[i] = p2[permpro[i]];
2015 for (i=1; i<lnz; i++) p2[i] = p1[i];
2016 /* perm indexes the rows of mat
2017 * |_0__|__redund__|__dense__|__too big__|_____done______|
2018 * 0 nlze lig li
2019 * \___nr___/ \___k0__/
2020 * \____________lnz ______________/
2021 *
2022 * col co
2023 * [dep | ]
2024 * i0 [--------| B ] (i0 = nlze + nr)
2025 * [matbnew | ] matbnew has maximal rank = lnz-1 - nr
2026 * mat = [--------|-----] lig
2027 * [ 0 | Id ]
2028 * [ | ] li */
2029
2030 matbnew = cgetg(col+1,t_MAT); /* dense+toobig, maximal rank. For hnffinal */
2031 dep = cgetg(col+1,t_MAT); /* rows dependent from the ones in matbnew */
2032 for (j=1; j<=col; j++)
2033 {
2034 GEN z = gel(extramat,j);
2035 p1 = cgetg(nlze+nr+1,t_COL); gel(dep,j) = p1;
2036 p2 = cgetg(lnz-nr,t_COL); gel(matbnew,j) = p2;
2037 for (i=1; i<=nlze; i++) gel(p1,i) = gen_0;
2038 p1 += nlze; for (i=1; i<=nr; i++) p1[i] = z[permpro[i]];
2039 p2 -= nr; for ( ; i<lnz; i++) p2[i] = z[permpro[i]];
2040 }
2041
2042 /* redundant generators in terms of the genuine generators
2043 * (x_i) = - (g_i) B */
2044 B = cgetg(co-col,t_MAT);
2045 for (j=col+1; j<co; j++)
2046 {
2047 GEN y = gel(matt,j);
2048 GEN z = gel(matb,j);
2049 p1=cgetg(lig+1,t_COL); gel(B,j-col) = p1;
2050 for (i=1; i<=nlze; i++) p1[i] = z[i];
2051 p1 += nlze; z += nlze-k0;
2052 for (k=1; k<lnz; k++)
2053 {
2054 i = permpro[k];
2055 p1[k] = (i <= k0)? y[i]: z[i];
2056 }
2057 }
2058 if (T) C = gmul(C,T);
2059 *ptdep = dep;
2060 *ptB = B;
2061 H = hnffinal(matbnew, perm, ptdep, ptB, &C);
2062 if (DEBUGLEVEL)
2063 msgtimer("hnfspec [%ld x %ld] --> [%ld x %ld]",li-1,co-1, lig-1,col-1);
2064 if (CO > co)
2065 { /* treat the rest, N cols at a time (hnflll slow otherwise) */
2066 const long N = 300;
2067 long a, L = CO - co, l = min(L, N); /* L columns to add */
2068 GEN CC = *ptC, m0 = (GEN)mat0;
2069 setlg(CC, CO); /* restore */
2070 CC += co-1;
2071 m0 += co-1;
2072 for (a = l;;)
2073 {
2074 GEN mat = cgetg(l + 1, t_MAT), emb = cgetg(l + 1, t_MAT);
2075 for (j = 1 ; j <= l; j++)
2076 {
2077 mat[j] = m0[j];
2078 emb[j] = CC[j];
2079 }
2080 H = hnfadd_i(H, perm, ptdep, ptB, &C, mat, emb);
2081 if (a == L) break;
2082 CC += l;
2083 m0 += l;
2084 a += l; if (a > L) { l = L - (a - l); a = L; }
2085 gerepileall(av, 4, &H,&C,ptB,ptdep);
2086 }
2087 }
2088 *ptC = C; return H;
2089 }
2090
2091 GEN
hnfspec(long ** mat,GEN perm,GEN * ptdep,GEN * ptB,GEN * ptC,long k0)2092 hnfspec(long** mat, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, long k0)
2093 {
2094 pari_sp av = avma;
2095 GEN H = hnfspec_i(mat, perm, ptdep, ptB, ptC, k0);
2096 gerepileall(av, 4, ptC, ptdep, ptB, &H); return H;
2097 }
2098
2099 /* HNF reduce x, apply same transforms to C */
2100 GEN
mathnfspec(GEN x,GEN * ptperm,GEN * ptdep,GEN * ptB,GEN * ptC)2101 mathnfspec(GEN x, GEN *ptperm, GEN *ptdep, GEN *ptB, GEN *ptC)
2102 {
2103 long i,j,k,ly,lx = lg(x);
2104 GEN p1,p2,z,perm;
2105 if (lx == 1) return gcopy(x);
2106 ly = lg(x[1]);
2107 z = cgetg(lx,t_MAT);
2108 perm = cgetg(ly,t_VECSMALL); *ptperm = perm;
2109 for (i=1; i<ly; i++) perm[i] = i;
2110 for (i=1; i<lx; i++)
2111 {
2112 p1 = cgetg(ly,t_COL); gel(z,i) = p1;
2113 p2 = gel(x,i);
2114 for (j=1; j<ly; j++)
2115 {
2116 if (is_bigint(p2[j])) goto TOOLARGE;
2117 p1[j] = itos(gel(p2,j));
2118 }
2119 }
2120 /* [ dep | ]
2121 * [-----| B ]
2122 * [ H | ]
2123 * [-----|-----]
2124 * [ 0 | Id ] */
2125 return hnfspec((long**)z,perm, ptdep, ptB, ptC, 0);
2126
2127 TOOLARGE:
2128 if (lg(*ptC) > 1 && lg((*ptC)[1]) > 1)
2129 pari_err(impl,"mathnfspec with large entries");
2130 x = hnf(x); lx = lg(x); j = ly; k = 0;
2131 for (i=1; i<ly; i++)
2132 {
2133 if (gcmp1(gcoeff(x,i,i + lx-ly)))
2134 perm[--j] = i;
2135 else
2136 perm[++k] = i;
2137 }
2138 setlg(perm,k+1);
2139 x = rowpermute(x, perm); /* upper part */
2140 setlg(perm,ly);
2141 *ptB = vecslice(x, j+lx-ly, lx-1);
2142 setlg(x, j);
2143 *ptdep = rowslice(x, 1, lx-ly);
2144 return rowslice(x, lx-ly+1, k); /* H */
2145 }
2146
2147 /* same as Flv_to_ZV, Flc_to_ZC, Flm_to_ZM but do not assume positivity */
2148 GEN
zx_to_ZX(GEN z)2149 zx_to_ZX(GEN z)
2150 {
2151 long i, l = lg(z);
2152 GEN x = cgetg(l,t_POL);
2153 for (i=2; i<l; i++) gel(x,i) = stoi(z[i]);
2154 x[1] = evalsigne(l-2!=0)| z[1]; return x;
2155 }
2156
2157 GEN
zm_to_ZM(GEN z)2158 zm_to_ZM(GEN z)
2159 {
2160 long i, l = lg(z);
2161 GEN x = cgetg(l,t_MAT);
2162 for (i=1; i<l; i++) gel(x,i) = zc_to_ZC(gel(z,i));
2163 return x;
2164 }
2165
2166 GEN
ZM_to_zm(GEN z)2167 ZM_to_zm(GEN z)
2168 {
2169 long i, l = lg(z);
2170 GEN x = cgetg(l,t_MAT);
2171 for (i=1; i<l; i++) gel(x,i) = ZV_to_zv(gel(z,i));
2172 return x;
2173 }
2174
2175 /* add new relations to a matrix treated by hnfspec (extramat / extraC) */
2176 GEN
hnfadd_i(GEN H,GEN perm,GEN * ptdep,GEN * ptB,GEN * ptC,GEN extramat,GEN extraC)2177 hnfadd_i(GEN H, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, /* cf hnfspec */
2178 GEN extramat,GEN extraC)
2179 {
2180 GEN matb, extratop, Cnew, permpro, B = *ptB, C = *ptC, dep = *ptdep;
2181 long i;
2182 long lH = lg(H)-1;
2183 long lB = lg(B)-1;
2184 long li = lg(perm)-1, lig = li - lB;
2185 long co = lg(C)-1, col = co - lB;
2186 long nlze = lH? lg(dep[1])-1: lg(B[1])-1;
2187
2188 /* col co
2189 * [ 0 |dep | ]
2190 * nlze [--------| B ]
2191 * [ 0 | H | ]
2192 * [--------|-----] lig
2193 * [ 0 | Id ]
2194 * [ | ] li */
2195 extratop = zm_to_ZM( rowslicepermute(extramat, perm, 1, lig) );
2196 if (li != lig)
2197 { /* zero out bottom part, using the Id block */
2198 GEN A = vecslice(C, col+1, co);
2199 GEN c = rowslicepermute(extramat, perm, lig+1, li);
2200 extraC = gsub(extraC, typ(A)==t_MAT? RgM_zm_mul(A, c): RgV_zm_mul(A,c));
2201 extratop = gsub(extratop, ZM_zm_mul(B, c));
2202 }
2203
2204 extramat = shallowconcat(extratop, vconcat(dep, H));
2205 Cnew = shallowconcat(extraC, vecslice(C, col-lH+1, co));
2206 if (DEBUGLEVEL>5) fprintferr(" 1st phase done\n");
2207 permpro = imagecomplspec(extramat, &nlze);
2208 extramat = rowpermute(extramat, permpro);
2209 *ptB = rowpermute(B, permpro);
2210 permpro = vecpermute(perm, permpro);
2211 for (i=1; i<=lig; i++) perm[i] = permpro[i]; /* perm o= permpro */
2212
2213 *ptdep = rowslice(extramat, 1, nlze);
2214 matb = rowslice(extramat, nlze+1, lig);
2215 if (DEBUGLEVEL>5) fprintferr(" 2nd phase done\n");
2216 H = hnffinal(matb,perm,ptdep,ptB,&Cnew);
2217 *ptC = shallowconcat(vecslice(C, 1, col-lH), Cnew);
2218 if (DEBUGLEVEL)
2219 {
2220 msgtimer("hnfadd (%ld + %ld)", lg(extratop)-1, lg(dep)-1);
2221 if (DEBUGLEVEL>7) fprintferr("H = %Z\nC = %Z\n",H,*ptC);
2222 }
2223 return H;
2224 }
2225
2226 GEN
hnfadd(GEN H,GEN perm,GEN * ptdep,GEN * ptB,GEN * ptC,GEN extramat,GEN extraC)2227 hnfadd(GEN H, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, /* cf hnfspec */
2228 GEN extramat,GEN extraC)
2229 {
2230 pari_sp av = avma;
2231 H = hnfadd_i(H, perm, ptdep, ptB, ptC, ZM_to_zm(extramat), extraC);
2232 gerepileall(av, 4, ptC, ptdep, ptB, &H); return H;
2233 }
2234
2235 static void
ZV_neg(GEN x)2236 ZV_neg(GEN x)
2237 {
2238 long i, lx = lg(x);
2239 for (i=1; i<lx ; i++) gel(x,i) = negi(gel(x,i));
2240 }
2241
2242 /* zero aj = Aij (!= 0) using ak = Aik (maybe 0), via linear combination of
2243 * A[j] and A[k] of determinant 1. If U != NULL, likewise update its columns */
2244 static void
ZV_elem(GEN aj,GEN ak,GEN A,GEN U,long j,long k)2245 ZV_elem(GEN aj, GEN ak, GEN A, GEN U, long j, long k)
2246 {
2247 GEN p1,u,v,d;
2248
2249 if (!signe(ak)) { lswap(A[j],A[k]); if (U) lswap(U[j],U[k]); return; }
2250 d = bezout(aj,ak,&u,&v);
2251 /* frequent special case (u,v) = (1,0) or (0,1) */
2252 if (!signe(u))
2253 { /* ak | aj */
2254 p1 = negi(diviiexact(aj,ak));
2255 gel(A,j) = ZV_lincomb(gen_1, p1, gel(A,j), gel(A,k));
2256 if (U)
2257 gel(U,j) = ZV_lincomb(gen_1, p1, gel(U,j), gel(U,k));
2258 return;
2259 }
2260 if (!signe(v))
2261 { /* aj | ak */
2262 p1 = negi(diviiexact(ak,aj));
2263 gel(A,k) = ZV_lincomb(gen_1, p1, gel(A,k), gel(A,j));
2264 lswap(A[j], A[k]);
2265 if (U) {
2266 gel(U,k) = ZV_lincomb(gen_1, p1, gel(U,k), gel(U,j));
2267 lswap(U[j], U[k]);
2268 }
2269 return;
2270 }
2271
2272 if (!is_pm1(d)) { aj = diviiexact(aj, d); ak = diviiexact(ak, d); }
2273 p1 = gel(A,k); aj = negi(aj);
2274 gel(A,k) = ZV_lincomb(u,v, gel(A,j),p1);
2275 gel(A,j) = ZV_lincomb(aj,ak, p1,gel(A,j));
2276 if (U)
2277 {
2278 p1 = gel(U,k);
2279 gel(U,k) = ZV_lincomb(u,v, gel(U,j),p1);
2280 gel(U,j) = ZV_lincomb(aj,ak, p1,gel(U,j));
2281 }
2282 }
2283
2284 /* reduce A[i,j] mod A[i,j0] for j=j0+1... via column operations */
2285 static void
ZM_reduce(GEN A,GEN U,long i,long j0)2286 ZM_reduce(GEN A, GEN U, long i, long j0)
2287 {
2288 long j, lA = lg(A);
2289 GEN d = gcoeff(A,i,j0);
2290 if (signe(d) < 0)
2291 {
2292 ZV_neg(gel(A,j0));
2293 if (U) ZV_neg(gel(U,j0));
2294 d = gcoeff(A,i,j0);
2295 }
2296 for (j=j0+1; j<lA; j++)
2297 {
2298 GEN q = truedivii(gcoeff(A,i,j), d);
2299 if (!signe(q)) continue;
2300
2301 q = negi(q);
2302 gel(A,j) = ZV_lincomb(gen_1,q, gel(A,j), gel(A,j0));
2303 if (U)
2304 gel(U,j) = ZV_lincomb(gen_1,q, gel(U,j), gel(U,j0));
2305 }
2306 }
2307
2308 /* A,B integral ideals in HNF. Return Au in Z^n (v in Z^n not computed), such
2309 * that Au + Bv = 1 */
2310 GEN
hnfmerge_get_1(GEN A,GEN B)2311 hnfmerge_get_1(GEN A, GEN B)
2312 {
2313 pari_sp av = avma;
2314 long j, k, c, l = lg(A), lb;
2315 GEN b, t, U = cgetg(l + 1, t_MAT), C = cgetg(l + 1, t_VEC);
2316
2317 t = NULL; /* -Wall */
2318 b = gcoeff(B,1,1); lb = lgefint(b);
2319 if (!signe(b)) {
2320 if (gcmp1(gcoeff(A,1,1))) return gscalcol_i(gen_1, l-1);
2321 l = 0; /* trigger error */
2322 }
2323 for (j = 1; j < l; j++)
2324 {
2325 c = j+1;
2326 gel(U,j) = col_ei(l-1, j);
2327 gel(U,c) = zerocol(l-1); /* dummy */
2328 gel(C,j) = vecslice(gel(A,j), 1,j);
2329 gel(C,c) = vecslice(gel(B,j), 1,j);
2330 for (k = j; k > 0; k--)
2331 {
2332 t = gcoeff(C,k,c);
2333 if (gcmp0(t)) continue;
2334 setlg(C[c], k+1);
2335 ZV_elem(t, gcoeff(C,k,k), C, U, c, k);
2336 if (lgefint(gcoeff(C,k,k)) > lb) gel(C,k) = FpC_red(gel(C,k), b);
2337 if (j > 4)
2338 {
2339 GEN u = gel(U,k);
2340 long h;
2341 for (h=1; h<l; h++)
2342 if (lgefint(u[h]) > lb) gel(u,h) = remii(gel(u,h), b);
2343 }
2344 }
2345 if (j == 1)
2346 t = gcoeff(C,1,1);
2347 else
2348 {
2349 GEN u,v;
2350 t = bezout(b, gcoeff(C,1,1), &u, &v); /* >= 0 */
2351 if (signe(v) && !gcmp1(v)) gel(U,1) = ZV_Z_mul(gel(U,1), v);
2352 gcoeff(C,1,1) = t;
2353 }
2354 if (signe(t) && is_pm1(t)) break;
2355 }
2356 if (j >= l) pari_err(talker, "non coprime ideals in hnfmerge");
2357 return gerepileupto(av, gmul(A,gel(U,1)));
2358
2359 }
2360
2361 /* remove: throw away lin.dep.columns */
2362 GEN
hnf0(GEN A,int remove)2363 hnf0(GEN A, int remove)
2364 {
2365 pari_sp av0 = avma, av, lim;
2366 long s,li,co,i,j,k,def,ldef;
2367 GEN denx, a;
2368
2369 if (typ(A) == t_VEC) return hnf_special(A,remove);
2370 A = init_hnf(A,&denx,&co,&li,&av);
2371 if (!A) return cgetg(1,t_MAT);
2372
2373 lim = stack_lim(av,1);
2374 def=co-1; ldef=(li>co)? li-co: 0;
2375 for (i=li-1; i>ldef; i--)
2376 {
2377 for (j=def-1; j; j--)
2378 {
2379 a = gcoeff(A,i,j);
2380 if (!signe(a)) continue;
2381
2382 /* zero a = Aij using b = Aik */
2383 k = (j==1)? def: j-1;
2384 ZV_elem(a,gcoeff(A,i,k), A,NULL, j,k);
2385
2386 if (low_stack(lim, stack_lim(av,1)))
2387 {
2388 if (DEBUGMEM>1) pari_warn(warnmem,"hnf[1]. i=%ld",i);
2389 A = gerepilecopy(av, A);
2390 }
2391 }
2392 s = signe(gcoeff(A,i,def));
2393 if (s)
2394 {
2395 if (s < 0) ZV_neg(gel(A,def));
2396 ZM_reduce(A, NULL, i,def);
2397 def--;
2398 }
2399 else
2400 if (ldef) ldef--;
2401 if (low_stack(lim, stack_lim(av,1)))
2402 {
2403 if (DEBUGMEM>1) pari_warn(warnmem,"hnf[2]. i=%ld",i);
2404 A = gerepilecopy(av, A);
2405 }
2406 }
2407 if (remove)
2408 { /* remove null columns */
2409 for (i=1,j=1; j<co; j++)
2410 if (!gcmp0(gel(A,j))) A[i++] = A[j];
2411 setlg(A,i);
2412 }
2413 A = denx? gdiv(A,denx): ZM_copy(A);
2414 return gerepileupto(av0, A);
2415 }
2416
2417 GEN
hnf(GEN x)2418 hnf(GEN x) { return hnf0(x,1); }
2419
2420 static GEN
ZC_Cei(long n,long i,GEN c)2421 ZC_Cei(long n, long i, GEN c) { GEN e = zerocol(n); gel(e,i) = c; return e; }
2422
2423 enum { hnf_MODID = 1, hnf_PART = 2 };
2424
2425 /* u*z[1..k] mod b, in place */
2426 static void
FpV_Fp_mul_part_ip(GEN z,GEN u,GEN p,long k)2427 FpV_Fp_mul_part_ip(GEN z, GEN u, GEN p, long k)
2428 {
2429 long i;
2430 if (is_pm1(u)) {
2431 if (signe(u) > 0) {
2432 for (i = 1; i <= k; i++)
2433 if (signe(z[i])) gel(z,i) = modii(gel(z,i), p);
2434 } else {
2435 for (i = 1; i <= k; i++)
2436 if (signe(z[i])) gel(z,i) = modii(negi(gel(z,i)), p);
2437 }
2438 }
2439 else {
2440 for (i = 1; i <= k; i++)
2441 if (signe(z[i])) gel(z,i) = modii(mulii(u,gel(z,i)), p);
2442 }
2443 }
2444 static void
FpV_red_part_ip(GEN z,GEN p,long k)2445 FpV_red_part_ip(GEN z, GEN p, long k)
2446 {
2447 long i;
2448 for (i = 1; i <= k; i++) gel(z,i) = modii(gel(z,i), p);
2449 }
2450 /* dm = multiple of diag element (usually detint(x))
2451 * flag & MODID: reduce mod dm * matid [ otherwise as above ].
2452 * flag & PART: don't reduce once diagonal is known; */
2453 static GEN
allhnfmod(GEN x,GEN dm,int flag)2454 allhnfmod(GEN x, GEN dm, int flag)
2455 {
2456 pari_sp av, lim;
2457 const int modid = (flag & hnf_MODID);
2458 long li, co, i, j, k, def, ldef, ldm;
2459 GEN a, b, p1, p2, u, v;
2460
2461 if (typ(dm)!=t_INT) pari_err(typeer,"allhnfmod");
2462 if (!signe(dm)) return hnf(x);
2463 if (typ(x)!=t_MAT) pari_err(typeer,"allhnfmod");
2464
2465 co = lg(x); if (co == 1) return cgetg(1,t_MAT);
2466 li = lg(x[1]);
2467 av = avma; lim = stack_lim(av,1);
2468 x = shallowcopy(x);
2469
2470 ldef = 0;
2471 if (li > co)
2472 {
2473 ldef = li - co;
2474 if (!modid) pari_err(talker,"nb lines > nb columns in hnfmod");
2475 }
2476 /* To prevent coeffs explosion, only reduce mod dm when lg() > ldm */
2477 ldm = lgefint(dm);
2478 for (def = co-1,i = li-1; i > ldef; i--,def--)
2479 {
2480 gcoeff(x,i,def) = remii(gcoeff(x,i,def), dm);
2481 for (j = def-1; j; j--)
2482 {
2483 gcoeff(x,i,j) = remii(gcoeff(x,i,j), dm);
2484 a = gcoeff(x,i,j);
2485 if (!signe(a)) continue;
2486
2487 k = (j==1)? def: j-1;
2488 gcoeff(x,i,k) = remii(gcoeff(x,i,k), dm);
2489 ZV_elem(a,gcoeff(x,i,k), x,NULL, j,k);
2490 p1 = gel(x,j);
2491 p2 = gel(x,k);
2492 for (k = 1; k < i; k++)
2493 {
2494 if (lgefint(p1[k]) > ldm) gel(p1,k) = remii(gel(p1,k), dm);
2495 if (lgefint(p2[k]) > ldm) gel(p2,k) = remii(gel(p2,k), dm);
2496 }
2497 if (low_stack(lim, stack_lim(av,1)))
2498 {
2499 if (DEBUGMEM>1) pari_warn(warnmem,"allhnfmod[1]. i=%ld",i);
2500 x = gerepilecopy(av, x);
2501 }
2502 }
2503 if (modid && !signe(gcoeff(x,i,def)))
2504 { /* missing pivot on line i, insert column */
2505 GEN a = cgetg(co + 1, t_MAT);
2506 for (k = 1; k <= def; k++) a[k] = x[k];
2507 gel(a,k++) = ZC_Cei(li-1, i, dm);
2508 for ( ; k <= co; k++) a[k] = x[k-1];
2509 ldef--; if (ldef < 0) ldef = 0;
2510 co++; def++; x = a;
2511 }
2512 }
2513 x += co - li;
2514 if (modid)
2515 { /* w[li] is an accumulator, discarded at the end */
2516 GEN w = cgetg(li+1, t_MAT);
2517 for (i = li-1; i > ldef; i--) w[i] = x[i];
2518 for ( ; i > 0; i--) gel(w,i) = ZC_Cei(li-1, i, dm);
2519 x = w;
2520 for (i = li-1; i > 0; i--)
2521 { /* check that dm*Id \subset L + add up missing dm*Id components */
2522 GEN d, c;
2523 d = bezout(gcoeff(x,i,i),dm, &u,&v);
2524 gcoeff(x,i,i) = d;
2525 if (is_pm1(d))
2526 {
2527 FpV_Fp_mul_part_ip(gel(x,i), u, dm, i-1);
2528 continue;
2529 }
2530 c = cgetg(li, t_COL);
2531 for (j = 1; j < i; j++) gel(c,j) = remii(gcoeff(x,j,i),d);
2532 for ( ; j <li; j++) gel(c,j) = gen_0;
2533 if (!equalii(dm, d)) c = gmul(c, diviiexact(dm, d));
2534 gel(x,li) = c;
2535 FpV_Fp_mul_part_ip(gel(x,i), u, dm, i-1);
2536 for (j = i - 1; j > ldef; j--)
2537 {
2538 GEN a = gcoeff(x, j, li);
2539 if (!signe(a)) continue;
2540 ZV_elem(a, gcoeff(x,j,j), x, NULL, li,j);
2541 FpV_red_part_ip(gel(x,li), dm, j-1);
2542 FpV_red_part_ip(gel(x,j), dm, j-1);
2543 if (low_stack(lim, stack_lim(av,1)))
2544 {
2545 if (DEBUGMEM>1) pari_warn(warnmem,"allhnfmod[2]. i=%ld", i);
2546 x = gerepilecopy(av, x);
2547 }
2548 }
2549 }
2550 }
2551 else
2552 {
2553 b = dm;
2554 for (i = li-1; i > 0; i--)
2555 {
2556 GEN d = bezout(gcoeff(x,i,i),b, &u,&v);
2557 gcoeff(x,i,i) = d;
2558 FpV_Fp_mul_part_ip(gel(x,i), u, b, i-1);
2559 if (i > 1) b = diviiexact(b,d);
2560 }
2561 }
2562 x[0] = evaltyp(t_MAT) | evallg(li); /* kill 0 columns / discard accumulator */
2563 if (flag & hnf_PART) return x;
2564
2565 if (modid) b = const_vec(li-1, dm);
2566 else
2567 { /* compute optimal value for dm */
2568 b = cgetg(li, t_VEC); gel(b,1) = gcoeff(x,1,1);
2569 for (i = 2; i < li; i++) gel(b,i) = mulii(gel(b,i-1), gcoeff(x,i,i));
2570 }
2571 dm = b;
2572
2573 for (i = li-2; i > 0; i--)
2574 {
2575 GEN diag = gcoeff(x,i,i);
2576 ldm = lgefint(dm[i]);
2577 for (j = i+1; j < li; j++)
2578 {
2579 b = negi(truedivii(gcoeff(x,i,j), diag));
2580 p1 = ZV_lincomb(gen_1,b, gel(x,j), gel(x,i));
2581 gel(x,j) = p1;
2582 for (k=1; k<i; k++)
2583 if (lgefint(p1[k]) > ldm) gel(p1,k) = remii(gel(p1,k), gel(dm,i));
2584 if (low_stack(lim, stack_lim(av,1)))
2585 {
2586 if (DEBUGMEM>1) pari_warn(warnmem,"allhnfmod[2]. i=%ld", i);
2587 gerepileall(av, 2, &x, &dm); diag = gcoeff(x,i,i);
2588 }
2589 }
2590 }
2591 return gerepilecopy(av, x);
2592 }
2593
2594 GEN
hnfmod(GEN x,GEN detmat)2595 hnfmod(GEN x, GEN detmat) { return allhnfmod(x,detmat, 0); }
2596 GEN
hnfmodid(GEN x,GEN p)2597 hnfmodid(GEN x, GEN p) { return allhnfmod(x, p, hnf_MODID); }
2598 GEN
hnfmodidpart(GEN x,GEN p)2599 hnfmodidpart(GEN x, GEN p) { return allhnfmod(x, p, hnf_MODID|hnf_PART); }
2600
2601 /***********************************************************************/
2602 /* */
2603 /* HNFLLL (Havas, Majewski, Mathews) */
2604 /* */
2605 /***********************************************************************/
2606
2607 static void
Minus(long j,GEN ** lambda)2608 Minus(long j, GEN **lambda)
2609 {
2610 long k, n = lg(lambda);
2611
2612 for (k=1 ; k<j; k++) lambda[j][k] = mynegi(lambda[j][k]);
2613 for (k=j+1; k<n; k++) lambda[k][j] = mynegi(lambda[k][j]);
2614 }
2615
2616 /* index of first non-zero entry */
2617 static long
findi(GEN M)2618 findi(GEN M)
2619 {
2620 long i, n = lg(M);
2621 for (i=1; i<n; i++)
2622 if (signe(M[i])) return i;
2623 return 0;
2624 }
2625
2626 static long
findi_normalize(GEN Aj,GEN B,long j,GEN ** lambda)2627 findi_normalize(GEN Aj, GEN B, long j, GEN **lambda)
2628 {
2629 long r = findi(Aj);
2630 if (r && signe(Aj[r]) < 0)
2631 {
2632 ZV_neg_ip(Aj); if (B) ZV_neg_ip(gel(B,j));
2633 Minus(j,lambda);
2634 }
2635 return r;
2636 }
2637
2638 static void
reduce2(GEN A,GEN B,long k,long j,long * row0,long * row1,GEN ** lambda,GEN * D)2639 reduce2(GEN A, GEN B, long k, long j, long *row0, long *row1, GEN **lambda, GEN *D)
2640 {
2641 GEN q;
2642 long i;
2643
2644 *row0 = findi_normalize(gel(A,j), B,j,lambda);
2645 *row1 = findi_normalize(gel(A,k), B,k,lambda);
2646 if (*row0)
2647 q = truedivii(gcoeff(A,*row0,k), gcoeff(A,*row0,j));
2648 else if (absi_cmp(shifti(lambda[k][j], 1), D[j]) > 0)
2649 q = diviiround(lambda[k][j], D[j]);
2650 else
2651 return;
2652
2653 if (signe(q))
2654 {
2655 GEN *Lk = lambda[k], *Lj = lambda[j];
2656 q = mynegi(q);
2657 if (*row0) elt_col(gel(A,k),gel(A,j),q);
2658 if (B) elt_col(gel(B,k),gel(B,j),q);
2659 Lk[j] = addii(Lk[j], mulii(q,D[j]));
2660 if (is_pm1(q))
2661 {
2662 if (signe(q) > 0)
2663 {
2664 for (i=1; i<j; i++)
2665 if (signe(Lj[i])) Lk[i] = addii(Lk[i], Lj[i]);
2666 }
2667 else
2668 {
2669 for (i=1; i<j; i++)
2670 if (signe(Lj[i])) Lk[i] = subii(Lk[i], Lj[i]);
2671 }
2672 }
2673 else
2674 {
2675 for (i=1; i<j; i++)
2676 if (signe(Lj[i])) Lk[i] = addii(Lk[i], mulii(q,Lj[i]));
2677 }
2678 }
2679 }
2680
2681 static void
hnfswap(GEN A,GEN B,long k,GEN ** lambda,GEN * D)2682 hnfswap(GEN A, GEN B, long k, GEN **lambda, GEN *D)
2683 {
2684 GEN t,p1,p2;
2685 long i,j,n = lg(A);
2686
2687 lswap(A[k],A[k-1]);
2688 if (B) lswap(B[k],B[k-1]);
2689 for (j=k-2; j; j--) swap(lambda[k-1][j], lambda[k][j]);
2690 for (i=k+1; i<n; i++)
2691 {
2692 p1 = mulii(lambda[i][k-1], D[k]);
2693 p2 = mulii(lambda[i][k], lambda[k][k-1]);
2694 t = subii(p1,p2);
2695
2696 p1 = mulii(lambda[i][k], D[k-2]);
2697 p2 = mulii(lambda[i][k-1], lambda[k][k-1]);
2698 lambda[i][k-1] = diviiexact(addii(p1,p2), D[k-1]);
2699
2700 lambda[i][k] = diviiexact(t, D[k-1]);
2701 }
2702 p1 = mulii(D[k-2],D[k]);
2703 p2 = sqri(lambda[k][k-1]);
2704 D[k-1] = diviiexact(addii(p1,p2), D[k-1]);
2705 }
2706
2707 /* reverse row order in matrix A */
2708 static GEN
fix_rows(GEN A)2709 fix_rows(GEN A)
2710 {
2711 long i,j, h,n = lg(A);
2712 GEN cB,cA,B = cgetg(n,t_MAT);
2713 if (n == 1) return B;
2714 h = lg(A[1]);
2715 for (j=1; j<n; j++)
2716 {
2717 cB = cgetg(h, t_COL);
2718 cA = gel(A,j); gel(B,j) = cB;
2719 for (i=h>>1; i; i--) { cB[h-i] = cA[i]; cB[i] = cA[h-i]; }
2720 }
2721 return B;
2722 }
2723
2724 GEN
hnflll_i(GEN A,GEN * ptB,int remove)2725 hnflll_i(GEN A, GEN *ptB, int remove)
2726 {
2727 pari_sp av = avma, lim = stack_lim(av,3);
2728 long m1 = 1, n1 = 1; /* alpha = m1/n1. Maybe 3/4 here ? */
2729 long do_swap,i,n,k;
2730 GEN z,B, **lambda, *D;
2731
2732 if (typ(A) != t_MAT) pari_err(typeer,"hnflll");
2733 n = lg(A);
2734 A = ZM_copy(fix_rows(A));
2735 B = ptB? matid(n-1): NULL;
2736 D = (GEN*)cgetg(n+1,t_VEC); lambda = (GEN**) cgetg(n,t_MAT);
2737 D++;
2738 for (i=0; i<n; i++) D[i] = gen_1;
2739 for (i=1; i<n; i++) lambda[i] = (GEN*)zerocol(n-1);
2740 k = 2;
2741 while (k < n)
2742 {
2743 long row0, row1;
2744 reduce2(A,B,k,k-1,&row0,&row1,lambda,D);
2745 if (row0)
2746 do_swap = (!row1 || row0 <= row1);
2747 else if (!row1)
2748 { /* row0 == row1 == 0 */
2749 pari_sp av1 = avma;
2750 z = addii(mulii(D[k-2],D[k]), sqri(lambda[k][k-1]));
2751 do_swap = (cmpii(mulsi(n1,z), mulsi(m1,sqri(D[k-1]))) < 0);
2752 avma = av1;
2753 }
2754 else
2755 do_swap = 0;
2756 if (do_swap)
2757 {
2758 hnfswap(A,B,k,lambda,D);
2759 if (k > 2) k--;
2760 }
2761 else
2762 {
2763 for (i=k-2; i; i--)
2764 {
2765 long row0, row1;
2766 reduce2(A,B,k,i,&row0,&row1,lambda,D);
2767 if (low_stack(lim, stack_lim(av,3)))
2768 {
2769 GEN b = (GEN)(D-1);
2770 if (DEBUGMEM) pari_warn(warnmem,"hnflll (reducing), i = %ld",i);
2771 gerepileall(av, B? 4: 3, &A, &lambda, &b, &B);
2772 D = (GEN*)(b+1);
2773 }
2774 }
2775 k++;
2776 }
2777 if (low_stack(lim, stack_lim(av,3)))
2778 {
2779 GEN b = (GEN)(D-1);
2780 if (DEBUGMEM) pari_warn(warnmem,"hnflll, k = %ld / %ld",k,n-1);
2781 gerepileall(av, B? 4: 3, &A, &lambda, &b, &B);
2782 D = (GEN*)(b+1);
2783 }
2784 }
2785 /* handle trivial case: return negative diag coefficient otherwise */
2786 if (n == 2) (void)findi_normalize(gel(A,1), B,1,lambda);
2787 A = fix_rows(A);
2788 if (remove)
2789 {
2790 for (i = 1; i < n; i++)
2791 if (findi(gel(A,i))) break;
2792 i--;
2793 A += i; A[0] = evaltyp(t_MAT) | evallg(n-i);
2794 }
2795 gerepileall(av, B? 2: 1, &A, &B);
2796 if (B) *ptB = B;
2797 return A;
2798 }
2799
2800 GEN
hnflll(GEN A)2801 hnflll(GEN A)
2802 {
2803 GEN B, z = cgetg(3, t_VEC);
2804 gel(z,1) = hnflll_i(A, &B, 0);
2805 gel(z,2) = B; return z;
2806 }
2807
2808 /* Variation on HNFLLL: Extended GCD */
2809
2810 static void
reduce1(GEN A,GEN B,long k,long j,GEN ** lambda,GEN * D)2811 reduce1(GEN A, GEN B, long k, long j, GEN **lambda, GEN *D)
2812 {
2813 GEN q;
2814 long i;
2815
2816 if (signe(A[j]))
2817 q = diviiround(gel(A,k),gel(A,j));
2818 else if (absi_cmp(shifti(lambda[k][j], 1), D[j]) > 0)
2819 q = diviiround(lambda[k][j], D[j]);
2820 else
2821 return;
2822
2823 if (signe(q))
2824 {
2825 GEN *Lk = lambda[k], *Lj = lambda[j];
2826 q = mynegi(q);
2827 gel(A,k) = addii(gel(A,k), mulii(q,gel(A,j)));
2828 elt_col(gel(B,k),gel(B,j),q);
2829 Lk[j] = addii(Lk[j], mulii(q,D[j]));
2830 for (i=1; i<j; i++)
2831 if (signe(Lj[i])) Lk[i] = addii(Lk[i], mulii(q,Lj[i]));
2832 }
2833 }
2834
2835 GEN
extendedgcd(GEN A)2836 extendedgcd(GEN A)
2837 {
2838 long m1 = 1, n1 = 1; /* alpha = m1/n1. Maybe 3/4 here ? */
2839 pari_sp av = avma;
2840 long do_swap,i,n,k;
2841 GEN z,B, **lambda, *D;
2842
2843 n = lg(A);
2844 for (i=1; i<n; i++)
2845 if (typ(A[i]) != t_INT) pari_err(typeer,"extendedgcd");
2846 A = shallowcopy(A);
2847 B = matid(n-1);
2848 D = (GEN*)new_chunk(n); lambda = (GEN**) cgetg(n,t_MAT);
2849 for (i=0; i<n; i++) D[i] = gen_1;
2850 for (i=1; i<n; i++) lambda[i] = (GEN*)zerocol(n-1);
2851 k = 2;
2852 while (k < n)
2853 {
2854 reduce1(A,B,k,k-1,lambda,D);
2855 if (signe(A[k-1])) do_swap = 1;
2856 else if (!signe(A[k]))
2857 {
2858 pari_sp av1 = avma;
2859 z = addii(mulii(D[k-2],D[k]), sqri(lambda[k][k-1]));
2860 do_swap = (cmpii(mulsi(n1,z), mulsi(m1,sqri(D[k-1]))) < 0);
2861 avma = av1;
2862 }
2863 else do_swap = 0;
2864
2865 if (do_swap)
2866 {
2867 hnfswap(A,B,k,lambda,D);
2868 if (k > 2) k--;
2869 }
2870 else
2871 {
2872 for (i=k-2; i; i--) reduce1(A,B,k,i,lambda,D);
2873 k++;
2874 }
2875 }
2876 if (signe(A[n-1]) < 0)
2877 {
2878 gel(A,n-1) = mynegi(gel(A,n-1));
2879 ZV_neg_ip(gel(B,n-1));
2880 }
2881 return gerepilecopy(av, mkvec2(gel(A,n-1), B));
2882 }
2883
2884 /* HNF with permutation. */
2885 GEN
hnfperm_i(GEN A,GEN * ptU,GEN * ptperm)2886 hnfperm_i(GEN A, GEN *ptU, GEN *ptperm)
2887 {
2888 GEN U, c, l, perm, d, p, q, b;
2889 pari_sp av = avma, av1, lim;
2890 long r, t, i, j, j1, k, m, n;
2891
2892 if (typ(A) != t_MAT) pari_err(typeer,"hnfperm");
2893 n = lg(A)-1;
2894 if (!n)
2895 {
2896 if (ptU) *ptU = cgetg(1,t_MAT);
2897 if (ptperm) *ptperm = cgetg(1,t_VEC);
2898 return cgetg(1, t_MAT);
2899 }
2900 m = lg(A[1])-1;
2901 c = const_vecsmall(m, 0);
2902 l = const_vecsmall(n, 0);
2903 perm = cgetg(m+1, t_VECSMALL);
2904 av1 = avma; lim = stack_lim(av1,1);
2905 A = shallowcopy(A);
2906 U = ptU? matid(n): NULL;
2907 /* U base change matrix : A0*U = A all along */
2908 for (r=0, k=1; k <= n; k++)
2909 {
2910 for (j=1; j<k; j++)
2911 {
2912 if (!l[j]) continue;
2913 t=l[j]; b=gcoeff(A,t,k);
2914 if (!signe(b)) continue;
2915
2916 ZV_elem(b,gcoeff(A,t,j), A,U,k,j);
2917 d = gcoeff(A,t,j);
2918 if (signe(d) < 0)
2919 {
2920 ZV_neg(gel(A,j));
2921 if (U) ZV_neg(gel(U,j));
2922 d = gcoeff(A,t,j);
2923 }
2924 for (j1=1; j1<j; j1++)
2925 {
2926 if (!l[j1]) continue;
2927 q = truedivii(gcoeff(A,t,j1),d);
2928 if (!signe(q)) continue;
2929
2930 q = negi(q);
2931 gel(A,j1) = ZV_lincomb(gen_1,q,gel(A,j1),gel(A,j));
2932 if (U) gel(U,j1) = ZV_lincomb(gen_1,q,gel(U,j1),gel(U,j));
2933 }
2934 }
2935 t = m; while (t && (c[t] || !signe(gcoeff(A,t,k)))) t--;
2936 if (t)
2937 {
2938 p = gcoeff(A,t,k);
2939 for (i=t-1; i; i--)
2940 {
2941 q = gcoeff(A,i,k);
2942 if (signe(q) && absi_cmp(p,q) > 0) { p = q; t = i; }
2943 }
2944 perm[++r] = l[k] = t; c[t] = k;
2945 if (signe(p) < 0)
2946 {
2947 ZV_neg(gel(A,k));
2948 if (U) ZV_neg(gel(U,k));
2949 p = gcoeff(A,t,k);
2950 }
2951 /* p > 0 */
2952 for (j=1; j<k; j++)
2953 {
2954 if (!l[j]) continue;
2955 q = truedivii(gcoeff(A,t,j),p);
2956 if (!signe(q)) continue;
2957
2958 q = negi(q);
2959 gel(A,j) = ZV_lincomb(gen_1,q,gel(A,j),gel(A,k));
2960 if (U) gel(U,j) = ZV_lincomb(gen_1,q,gel(U,j),gel(U,k));
2961 }
2962 }
2963 if (low_stack(lim, stack_lim(av1,1)))
2964 {
2965 if (DEBUGMEM>1) pari_warn(warnmem,"hnfperm");
2966 gerepileall(av1, U? 2: 1, &A, &U);
2967 }
2968 }
2969 if (r < m)
2970 {
2971 for (i=1,k=r; i<=m; i++)
2972 if (!c[i]) perm[++k] = i;
2973 }
2974
2975 /* We have A0*U=A, U in Gl(n,Z)
2976 * basis for Im(A): columns of A s.t l[j]>0 (r cols)
2977 * basis for Ker(A): columns of U s.t l[j]=0 (n-r cols) */
2978 p = cgetg(r+1,t_MAT);
2979 for (i=1; i<=m/2; i++) lswap(perm[i], perm[m+1-i]);
2980 if (U)
2981 {
2982 GEN u = cgetg(n+1,t_MAT);
2983 for (t=1,k=r,j=1; j<=n; j++)
2984 if (l[j])
2985 {
2986 u[k + n-r] = U[j];
2987 gel(p,k--) = vecpermute(gel(A,j), perm);
2988 }
2989 else
2990 u[t++] = U[j];
2991 *ptU = u;
2992 *ptperm = perm;
2993 gerepileall(av, 3, &p, ptU, ptperm);
2994 }
2995 else
2996 {
2997 for (k=r,j=1; j<=n; j++)
2998 if (l[j]) gel(p,k--) = vecpermute(gel(A,j), perm);
2999 if (ptperm) *ptperm = perm;
3000 gerepileall(av, ptperm? 2: 1, &p, ptperm);
3001 }
3002 return p;
3003 }
3004
3005 GEN
hnfperm(GEN A)3006 hnfperm(GEN A)
3007 {
3008 GEN U, perm, y = cgetg(4, t_VEC);
3009 gel(y,1) = hnfperm_i(A, &U, &perm);
3010 gel(y,2) = U;
3011 gel(y,3) = vecsmall_to_vec(perm); return y;
3012 }
3013
3014 /* Hermite Normal Form */
3015 GEN
hnfall_i(GEN A,GEN * ptB,long remove)3016 hnfall_i(GEN A, GEN *ptB, long remove)
3017 {
3018 GEN B,c,h,x,p,a;
3019 pari_sp av = avma, av1, lim;
3020 long m,n,r,i,j,k,li;
3021
3022 if (typ(A)!=t_MAT) pari_err(typeer,"hnfall");
3023 n=lg(A)-1;
3024 if (!n)
3025 {
3026 if (ptB) *ptB = cgetg(1,t_MAT);
3027 return cgetg(1,t_MAT);
3028 }
3029 m = lg(A[1])-1;
3030 c = cgetg(m+1,t_VECSMALL); for (i=1; i<=m; i++) c[i]=0;
3031 h = cgetg(n+1,t_VECSMALL); for (j=1; j<=n; j++) h[j]=m;
3032 av1 = avma; lim = stack_lim(av1,1);
3033 A = shallowcopy(A);
3034 B = ptB? matid(n): NULL;
3035 r = n+1;
3036 for (li=m; li; li--)
3037 {
3038 for (j=1; j<r; j++)
3039 {
3040 for (i=h[j]; i>li; i--)
3041 {
3042 a = gcoeff(A,i,j);
3043 if (!signe(a)) continue;
3044
3045 k = c[i]; /* zero a = Aij using Aik */
3046 ZV_elem(a,gcoeff(A,i,k), A,B,j,k);
3047 ZM_reduce(A,B, i,k);
3048 if (low_stack(lim, stack_lim(av1,1)))
3049 {
3050 if (DEBUGMEM>1) pari_warn(warnmem,"hnfall[1], li = %ld", li);
3051 gerepileall(av1, B? 2: 1, &A, &B);
3052 }
3053 }
3054 x = gcoeff(A,li,j); if (signe(x)) break;
3055 h[j] = li-1;
3056 }
3057 if (j == r) continue;
3058 r--;
3059 if (j < r) /* A[j] != 0 */
3060 {
3061 lswap(A[j], A[r]);
3062 if (B) lswap(B[j], B[r]);
3063 h[j]=h[r]; h[r]=li; c[li]=r;
3064 }
3065 p = gcoeff(A,li,r);
3066 if (signe(p) < 0)
3067 {
3068 ZV_neg(gel(A,r));
3069 if (B) ZV_neg(gel(B,r));
3070 p = gcoeff(A,li,r);
3071 }
3072 ZM_reduce(A,B, li,r);
3073 if (low_stack(lim, stack_lim(av1,1)))
3074 {
3075 if (DEBUGMEM>1) pari_warn(warnmem,"hnfall[2], li = %ld", li);
3076 gerepileall(av1, B? 2: 1, &A, &B);
3077 }
3078 }
3079 if (DEBUGLEVEL>5) fprintferr("\nhnfall, final phase: ");
3080 r--; /* first r cols are in the image the n-r (independent) last ones */
3081 for (j=1; j<=r; j++)
3082 for (i=h[j]; i; i--)
3083 {
3084 a = gcoeff(A,i,j);
3085 k = c[i];
3086 if (signe(a)) ZV_elem(a,gcoeff(A,i,k), A,B, j,k);
3087 ZM_reduce(A,B, i,k); /* ensure non-negative entries, even if a = 0 */
3088 if (low_stack(lim, stack_lim(av1,1)))
3089 {
3090 if (DEBUGMEM>1) pari_warn(warnmem,"hnfall[3], j = %ld", j);
3091 gerepileall(av1, B? 2: 1, &A, &B);
3092 }
3093 }
3094 if (DEBUGLEVEL>5) fprintferr("\n");
3095 /* remove the first r columns */
3096 if (remove) { A += r; A[0] = evaltyp(t_MAT) | evallg(n-r+1); }
3097 gerepileall(av, B? 2: 1, &A, &B);
3098 if (B) *ptB = B;
3099 return A;
3100 }
3101
3102 GEN
hnfall0(GEN A,long remove)3103 hnfall0(GEN A, long remove)
3104 {
3105 GEN B, z = cgetg(3, t_VEC);
3106 gel(z,1) = hnfall_i(A, &B, remove);
3107 gel(z,2) = B; return z;
3108 }
3109
3110 GEN
hnfall(GEN x)3111 hnfall(GEN x) {return hnfall0(x,1);}
3112
3113 /***************************************************************/
3114 /** **/
3115 /** SMITH NORMAL FORM REDUCTION **/
3116 /** **/
3117 /***************************************************************/
3118
3119 static GEN
col_mul(GEN x,GEN c)3120 col_mul(GEN x, GEN c)
3121 {
3122 if (typ(x) == t_INT)
3123 {
3124 long s = signe(x);
3125 GEN xc = NULL;
3126 if (s)
3127 {
3128 if (!is_pm1(x)) xc = gmul(x,c);
3129 else xc = (s>0)? c: gneg_i(c);
3130 }
3131 return xc;
3132 }
3133 else
3134 return gmul(x, c);
3135 }
3136
3137 static void
do_zero(GEN x)3138 do_zero(GEN x)
3139 {
3140 long i, lx = lg(x);
3141 for (i=1; i<lx; i++) gel(x,i) = gen_0;
3142 }
3143
3144 /* c1 <-- u.c1 + v.c2; c2 <-- a.c2 - b.c1 */
3145 static void
update(GEN u,GEN v,GEN a,GEN b,GEN * c1,GEN * c2)3146 update(GEN u, GEN v, GEN a, GEN b, GEN *c1, GEN *c2)
3147 {
3148 GEN p1,p2;
3149
3150 u = col_mul(u,*c1);
3151 v = col_mul(v,*c2);
3152 if (u) p1 = v? gadd(u,v): u;
3153 else p1 = v? v: (GEN)NULL;
3154
3155 a = col_mul(a,*c2);
3156 b = col_mul(gneg_i(b),*c1);
3157 if (a) p2 = b? gadd(a,b): a;
3158 else p2 = b? b: (GEN)NULL;
3159
3160 if (!p1) do_zero(*c1); else *c1 = p1;
3161 if (!p2) do_zero(*c2); else *c2 = p2;
3162 }
3163
3164 static GEN
trivsmith(long all)3165 trivsmith(long all)
3166 {
3167 GEN z;
3168 if (!all) return cgetg(1,t_VEC);
3169 z=cgetg(4,t_VEC);
3170 gel(z,1) = cgetg(1,t_MAT);
3171 gel(z,2) = cgetg(1,t_MAT);
3172 gel(z,3) = cgetg(1,t_MAT); return z;
3173 }
3174
3175 static void
snf_pile(pari_sp av,GEN * x,GEN * U,GEN * V)3176 snf_pile(pari_sp av, GEN *x, GEN *U, GEN *V)
3177 {
3178 GEN *gptr[3];
3179 int c = 1; gptr[0]=x;
3180 if (*U) gptr[c++] = U;
3181 if (*V) gptr[c++] = V;
3182 gerepilemany(av,gptr,c);
3183 }
3184
3185 static GEN
bezout_step(GEN * pa,GEN * pb,GEN * pu,GEN * pv)3186 bezout_step(GEN *pa, GEN *pb, GEN *pu, GEN *pv)
3187 {
3188 GEN a = *pa, b = *pb, d;
3189 if (absi_equal(a,b))
3190 {
3191 long sa = signe(a), sb = signe(b);
3192 *pv = gen_0;
3193 if (sb == sa) {
3194 *pa = *pb = gen_1;
3195 if (sa > 0) {*pu=gen_1; return a;} else {*pu=gen_m1; return absi(a);}
3196 }
3197 if (sa > 0) { *pa = *pu = gen_1; *pb = gen_m1; return a; }
3198 *pa = *pu = gen_m1; *pb = gen_1; return b;
3199 }
3200 d = bezout(a,b, pu,pv);
3201 *pa = diviiexact(a, d);
3202 *pb = diviiexact(b, d); return d;
3203 }
3204
3205 static int
negcmpii(GEN x,GEN y)3206 negcmpii(GEN x, GEN y) { return -cmpii(x,y); }
3207
3208 /* Return the SNF D of matrix X. If ptU/ptV non-NULL set them to U/V
3209 * to that D = UXV */
3210 GEN
smithall(GEN x,GEN * ptU,GEN * ptV)3211 smithall(GEN x, GEN *ptU, GEN *ptV)
3212 {
3213 pari_sp av0 = avma, av, lim = stack_lim(av0,1);
3214 long i, j, k, m0, m, n0, n;
3215 GEN p1, u, v, U, V, V0, mdet, ys, perm = NULL;
3216
3217 if (typ(x)!=t_MAT) pari_err(typeer,"smithall");
3218 n0 = n = lg(x)-1;
3219 if (!n) {
3220 if (ptU) *ptU = cgetg(1,t_MAT);
3221 if (ptV) *ptV = cgetg(1,t_MAT);
3222 return cgetg(1,(ptV||ptU)?t_MAT:t_VEC);
3223 }
3224 av = avma;
3225 m0 = m = lg(x[1])-1;
3226 for (j=1; j<=n; j++)
3227 for (i=1; i<=m; i++)
3228 if (typ(gcoeff(x,i,j)) != t_INT)
3229 pari_err(talker,"non integral matrix in smithall");
3230
3231 U = ptU? gen_1: NULL; /* TRANSPOSE of row transform matrix [act on columns]*/
3232 V = ptV? gen_1: NULL;
3233 V0 = NULL;
3234 x = shallowcopy(x);
3235 if (m == n && ZM_ishnf(x))
3236 {
3237 mdet = dethnf_i(x);
3238 if (V) *ptV = matid(n);
3239 }
3240 else
3241 {
3242 mdet = detint(x);
3243 if (signe(mdet))
3244 {
3245 if (!V)
3246 p1 = hnfmod(x,mdet);
3247 else
3248 {
3249 if (m == n)
3250 {
3251 p1 = hnfmod(x,mdet);
3252 *ptV = gauss(x,p1);
3253 }
3254 else
3255 p1 = hnfperm_i(x, ptV, ptU? &perm: NULL);
3256 }
3257 mdet = dethnf_i(p1);
3258 }
3259 else
3260 p1 = hnfperm_i(x, ptV, ptU? &perm: NULL);
3261 x = p1;
3262 }
3263 n = lg(x)-1;
3264 if (V)
3265 {
3266 V = *ptV;
3267 if (n != n0)
3268 {
3269 V0 = vecslice(V, 1, n0 - n); /* kernel */
3270 V = vecslice(V, n0-n+1, n0);
3271 av = avma;
3272 }
3273 }
3274 /* independent columns */
3275 if (!signe(mdet))
3276 {
3277 if (n)
3278 {
3279 x = smithall(shallowtrans(x), ptV, ptU); /* ptV, ptU swapped! */
3280 if (typ(x) == t_MAT && n != m) x = shallowtrans(x);
3281 if (V) V = gmul(V, shallowtrans(*ptV));
3282 if (U) U = *ptU; /* TRANSPOSE */
3283 }
3284 else /* 0 matrix */
3285 {
3286 x = cgetg(1,t_MAT);
3287 if (V) V = cgetg(1, t_MAT);
3288 if (U) U = matid(m);
3289 }
3290 goto THEEND;
3291 }
3292 if (U) U = matid(n);
3293
3294 /* square, maximal rank n */
3295 p1 = gen_sort(mattodiagonal_i(x), cmp_IND | cmp_C, &negcmpii);
3296 ys = cgetg(n+1,t_MAT);
3297 for (j=1; j<=n; j++) gel(ys,j) = vecpermute(gel(x, p1[j]), p1);
3298 x = ys;
3299 if (U) U = vecpermute(U, p1);
3300 if (V) V = vecpermute(V, p1);
3301
3302 p1 = hnfmod(x, mdet);
3303 if (V) V = gmul(V, gauss(x,p1));
3304 x = p1;
3305
3306 if (DEBUGLEVEL>7) fprintferr("starting SNF loop");
3307 for (i=n; i>1; i--)
3308 {
3309 if (DEBUGLEVEL>7) fprintferr("\ni = %ld: ",i);
3310 for(;;)
3311 {
3312 int c = 0;
3313 GEN a, b;
3314 for (j=i-1; j>=1; j--)
3315 {
3316 b = gcoeff(x,i,j); if (!signe(b)) continue;
3317 a = gcoeff(x,i,i);
3318 ZV_elem(b, a, x,V, j,i);
3319 if (low_stack(lim, stack_lim(av,1)))
3320 {
3321 if (DEBUGMEM>1) pari_warn(warnmem,"[1]: smithall i = %ld", i);
3322 snf_pile(av, &x,&U,&V);
3323 }
3324 }
3325 if (DEBUGLEVEL>7) fprintferr("; ");
3326 for (j=i-1; j>=1; j--)
3327 {
3328 GEN d;
3329 b = gcoeff(x,j,i); if (!signe(b)) continue;
3330 a = gcoeff(x,i,i);
3331 d = bezout_step(&a, &b, &u, &v);
3332 for (k = 1; k < i; k++)
3333 {
3334 GEN t = addii(mulii(u,gcoeff(x,i,k)),mulii(v,gcoeff(x,j,k)));
3335 gcoeff(x,j,k) = subii(mulii(a,gcoeff(x,j,k)),
3336 mulii(b,gcoeff(x,i,k)));
3337 gcoeff(x,i,k) = t;
3338 }
3339 gcoeff(x,j,i) = gen_0;
3340 gcoeff(x,i,i) = d;
3341 if (U) update(u,v,a,b,(GEN*)(U+i),(GEN*)(U+j));
3342 if (low_stack(lim, stack_lim(av,1)))
3343 {
3344 if (DEBUGMEM>1) pari_warn(warnmem,"[2]: smithall, i = %ld", i);
3345 snf_pile(av, &x,&U,&V);
3346 }
3347 c = 1;
3348 }
3349 if (!c)
3350 {
3351 b = gcoeff(x,i,i); if (!signe(b) || is_pm1(b)) break;
3352
3353 for (k=1; k<i; k++)
3354 {
3355 for (j=1; j<i; j++)
3356 if (signe(remii(gcoeff(x,k,j),b))) break;
3357 if (j != i) break;
3358 }
3359 if (k == i) break;
3360
3361 /* x[k,j] != 0 mod b */
3362 for (j=1; j<=i; j++)
3363 gcoeff(x,i,j) = addii(gcoeff(x,i,j),gcoeff(x,k,j));
3364 if (U) gel(U,i) = gadd(gel(U,i),gel(U,k));
3365 }
3366 if (low_stack(lim, stack_lim(av,1)))
3367 {
3368 if (DEBUGMEM>1) pari_warn(warnmem,"[3]: smithall");
3369 snf_pile(av, &x,&U,&V);
3370 }
3371 }
3372 }
3373 if (DEBUGLEVEL>7) fprintferr("\n");
3374 for (k=1; k<=n; k++)
3375 if (signe(gcoeff(x,k,k)) < 0)
3376 {
3377 if (V) gel(V,k) = gneg(gel(V,k));
3378 gcoeff(x,k,k) = negi(gcoeff(x,k,k));
3379 }
3380 THEEND:
3381 if (!U && !V)
3382 {
3383 if (typ(x) == t_MAT) x = mattodiagonal_i(x);
3384 m = lg(x)-1;
3385 if (m0 > m) x = shallowconcat(zerovec(m0-m), x);
3386 return gerepilecopy(av0, x);
3387 }
3388
3389 if (V0)
3390 {
3391 x = shallowconcat(zeromat(m,n0-n), x);
3392 if (V) V = shallowconcat(V0, V);
3393 }
3394 if (U)
3395 {
3396 U = shallowtrans(U);
3397 if (perm) U = vecpermute(U, perm_inv(perm));
3398 }
3399 snf_pile(av0, &x,&U,&V);
3400 if (ptU) *ptU = U;
3401 if (ptV) *ptV = V;
3402 return x;
3403 }
3404
3405 GEN
smith(GEN x)3406 smith(GEN x) { return smithall(x, NULL,NULL); }
3407
3408 GEN
smith2(GEN x)3409 smith2(GEN x)
3410 {
3411 GEN z = cgetg(4, t_VEC);
3412 gel(z,3) = smithall(x, (GEN*)(z+1),(GEN*)(z+2));
3413 return z;
3414 }
3415
3416 /* Assume z was computed by [g]smithall(). Remove the 1s on the diagonal */
3417 GEN
smithclean(GEN z)3418 smithclean(GEN z)
3419 {
3420 long i,j,l,c;
3421 GEN u,v,y,d,p1;
3422
3423 if (typ(z) != t_VEC) pari_err(typeer,"smithclean");
3424 l = lg(z); if (l == 1) return cgetg(1,t_VEC);
3425 u=gel(z,1);
3426 if (l != 4 || typ(u) != t_MAT)
3427 {
3428 for (c=1; c<l; c++)
3429 if (gcmp1(gel(z,c))) break;
3430 return gcopy_i(z, c);
3431 }
3432 v=gel(z,2); d=gel(z,3); l = lg(d);
3433 for (c=1; c<l; c++)
3434 if (gcmp1(gcoeff(d,c,c))) break;
3435
3436 y=cgetg(4,t_VEC);
3437 gel(y,1) = (p1 = cgetg(l,t_MAT));
3438 for (i=1; i<l; i++) gel(p1,i) = gcopy_i(gel(u,i), c);
3439 gel(y,2) = gcopy_i(v, c);
3440 gel(y,3) = (p1 = cgetg(c,t_MAT));
3441 for (i=1; i<c; i++)
3442 {
3443 GEN p2 = cgetg(c,t_COL); gel(p1,i) = p2;
3444 for (j=1; j<c; j++) gel(p2,j) = i==j? gcopy(gcoeff(d,i,i)): gen_0;
3445 }
3446 return y;
3447 }
3448
3449 static GEN
gbezout_step(GEN * pa,GEN * pb,GEN * pu,GEN * pv)3450 gbezout_step(GEN *pa, GEN *pb, GEN *pu, GEN *pv)
3451 {
3452 GEN a = *pa, b = *pb, d;
3453 if (gcmp0(a))
3454 {
3455 *pa = gen_0; *pu = gen_0;
3456 *pb = gen_1; *pv = gen_1; return b;
3457 }
3458 if (typ(a) == t_POL) a = RgX_renormalize(a);
3459 if (typ(b) == t_POL) b = RgX_renormalize(b);
3460 d = RgX_extgcd(a,b, pu,pv);
3461 if (typ(d) == t_POL)
3462 {
3463 if (degpol(d)) { a = RgX_div(a, d); b = RgX_div(b, d); }
3464 else if (typ(d[2]) == t_REAL && lg(d[2]) == 3)
3465 #if 1
3466 { /* possible accuracy problem */
3467 GEN D = RgX_gcd_simple(a,b);
3468 if (degpol(D)) {
3469 D = gdiv(D, leading_term(D));
3470 a = RgX_div(a, D); b = RgX_div(b, D);
3471 d = RgX_extgcd(a,b, pu,pv); /* retry now */
3472 d = gmul(d, D);
3473 }
3474 }
3475 #else
3476 { /* less stable */
3477 d = RgX_extgcd_simple(a,b, pu,pv);
3478 if (degpol(d)) { a = RgX_div(a, d); b = RgX_div(b, d); }
3479 }
3480 #endif
3481 }
3482 *pa = a;
3483 *pb = b; return d;
3484 }
3485
3486 static GEN
gsmithall(GEN x,long all)3487 gsmithall(GEN x,long all)
3488 {
3489 pari_sp av, lim;
3490 long i, j, k, n;
3491 GEN z, u, v, U, V;
3492
3493 if (typ(x)!=t_MAT) pari_err(typeer,"gsmithall");
3494 n = lg(x)-1;
3495 if (!n) return trivsmith(all);
3496 if (lg(x[1]) != n+1) pari_err(mattype1,"gsmithall");
3497 av = avma; lim = stack_lim(av,1);
3498 x = shallowcopy(x);
3499 if (all) { U = matid(n); V = matid(n); }
3500 for (i=n; i>=2; i--)
3501 {
3502 for(;;)
3503 {
3504 GEN a, b, d;
3505 int c = 0;
3506 for (j=i-1; j>=1; j--)
3507 {
3508 b = gcoeff(x,i,j); if (gcmp0(b)) continue;
3509 a = gcoeff(x,i,i);
3510 d = gbezout_step(&b, &a, &v, &u);
3511 for (k = 1; k < i; k++)
3512 {
3513 GEN t = gadd(gmul(u,gcoeff(x,k,i)),gmul(v,gcoeff(x,k,j)));
3514 gcoeff(x,k,j) = gsub(gmul(a,gcoeff(x,k,j)),gmul(b,gcoeff(x,k,i)));
3515 gcoeff(x,k,i) = t;
3516 }
3517 gcoeff(x,i,j) = gen_0;
3518 gcoeff(x,i,i) = d;
3519 if (all) update(u,v,a,b,(GEN*)(V+i),(GEN*)(V+j));
3520 }
3521 for (j=i-1; j>=1; j--)
3522 {
3523 b = gcoeff(x,j,i); if (gcmp0(b)) continue;
3524 a = gcoeff(x,i,i);
3525 d = gbezout_step(&b, &a, &v, &u);
3526 for (k = 1; k < i; k++)
3527 {
3528 GEN t = gadd(gmul(u,gcoeff(x,i,k)),gmul(v,gcoeff(x,j,k)));
3529 gcoeff(x,j,k) = gsub(gmul(a,gcoeff(x,j,k)),gmul(b,gcoeff(x,i,k)));
3530 gcoeff(x,i,k) = t;
3531 }
3532 gcoeff(x,j,i) = gen_0;
3533 gcoeff(x,i,i) = d;
3534 if (all) update(u,v,a,b,(GEN*)(U+i),(GEN*)(U+j));
3535 c = 1;
3536 }
3537 if (!c)
3538 {
3539 b = gcoeff(x,i,i); if (degpol(b) <= 0) break;
3540
3541 for (k=1; k<i; k++)
3542 {
3543 for (j=1; j<i; j++)
3544 {
3545 GEN r = gmod(gcoeff(x,k,j), b);
3546 if (signe(r) && (! isinexactreal(r) ||
3547 gexpo(r) > 16 + gexpo(b) - bit_accuracy(gprecision(r)))
3548 ) break;
3549 }
3550 if (j != i) break;
3551 }
3552 if (k == i) break;
3553
3554 for (j=1; j<=i; j++)
3555 gcoeff(x,i,j) = gadd(gcoeff(x,i,j),gcoeff(x,k,j));
3556 if (all) gel(U,i) = gadd(gel(U,i),gel(U,k));
3557 }
3558 if (low_stack(lim, stack_lim(av,1)))
3559 {
3560 if (DEBUGMEM>1) pari_warn(warnmem,"gsmithall");
3561 gerepileall(av, all? 3: 1, &x, &U, &V);
3562 }
3563 }
3564 }
3565 for (k=1; k<=n; k++)
3566 {
3567 GEN T = gcoeff(x,k,k);
3568 if (typ(T) == t_POL && signe(T))
3569 {
3570 GEN d = leading_term(T);
3571 while (gcmp0(d) || ( typ(d) == t_REAL && lg(d) == 3
3572 && gexpo(T) - expo(d) > (long)BITS_IN_LONG)) {
3573 T = normalizepol_i(T, lg(T)-1);
3574 if (!signe(T)) { gcoeff(x,k,k) = T; continue; }
3575 d = leading_term(T);
3576 }
3577 if (!gcmp1(d))
3578 {
3579 gcoeff(x,k,k) = gdiv(T,d);
3580 if (all) gel(V,k) = gdiv(gel(V,k), d);
3581 }
3582 }
3583 }
3584 z = all? mkvec3(shallowtrans(U), V, x): mattodiagonal_i(x);
3585 return gerepilecopy(av, z);
3586 }
3587
3588 GEN
matsnf0(GEN x,long flag)3589 matsnf0(GEN x,long flag)
3590 {
3591 pari_sp av = avma;
3592 if (flag > 7) pari_err(flagerr,"matsnf");
3593 if (typ(x) == t_VEC && flag & 4) return smithclean(x);
3594 if (flag & 2) x = flag&1 ? gsmith2(x): gsmith(x);
3595 else x = flag&1 ? smith2(x): smith(x);
3596 if (flag & 4) x = gerepileupto(av, smithclean(x));
3597 return x;
3598 }
3599
3600 GEN
gsmith(GEN x)3601 gsmith(GEN x) { return gsmithall(x,0); }
3602
3603 GEN
gsmith2(GEN x)3604 gsmith2(GEN x) { return gsmithall(x,1); }
3605
3606 /* H relation matrix among row of generators g in HNF. Let URV = D its SNF,
3607 * newU R newV = newD its clean SNF (no 1 in Dnew). Return the diagonal of
3608 * newD, newU and newUi such that 1/U = (newUi, ?).
3609 * Rationale: let (G,0) = g Ui be the new generators then
3610 * 0 = G U R --> G D = 0, g = G newU and G = g newUi */
3611 GEN
smithrel(GEN H,GEN * newU,GEN * newUi)3612 smithrel(GEN H, GEN *newU, GEN *newUi)
3613 {
3614 GEN U, V, Ui, D = smithall(H, &U, newUi? &V: NULL);
3615 long i, j, c, l = lg(D);
3616
3617 for (c=1; c<l; c++)
3618 {
3619 GEN t = gcoeff(D,c,c);
3620 if (is_pm1(t)) break;
3621 }
3622 setlg(D, c); D = mattodiagonal_i(D);
3623 if (newU) {
3624 U = rowslice(U, 1, c-1);
3625 for (i = 1; i < c; i++)
3626 {
3627 GEN d = gel(D,i), d2 = shifti(d, 1);
3628 for (j = 1; j < lg(U); j++)
3629 gcoeff(U,i,j) = centermodii(gcoeff(U,i,j), d, d2);
3630 }
3631 *newU = U;
3632 }
3633 if (newUi) { /* UHV = D --> U^-1 mod H = H(VD^-1 mod 1) mod H */
3634 if (c == 1) *newUi = cgetg(1, t_MAT);
3635 else
3636 { /* Ui = ZM_inv(U, gen_1); setlg(Ui, c); */
3637 setlg(V, c);
3638 V = FpM_red(V, gel(D,1));
3639 Ui = gmul(H, V);
3640 for (i = 1; i < c; i++)
3641 gel(Ui,i) = gdivexact(gel(Ui,i), gel(D,i));
3642 *newUi = reducemodHNF(Ui, H, NULL);
3643 }
3644 }
3645 return D;
3646 }
3647
3648 /***********************************************************************
3649 **** ****
3650 **** Frobenius form and Jordan form of a matrix ****
3651 **** ****
3652 ***********************************************************************/
3653
3654 GEN
Frobeniusform(GEN V,long n)3655 Frobeniusform(GEN V, long n)
3656 {
3657 long i,j,k;
3658 GEN M = cgetg(n+1, t_MAT);
3659 for(i=1; i<=n; i++)
3660 gel(M,i) = zerocol(n);
3661 for (k=1,i=1;i<lg(V);i++,k++)
3662 {
3663 GEN P = gel(V,i);
3664 long d = degpol(P);
3665 if (k+d-1 > n) pari_err(talker, "accuracy lost in matfrobenius");
3666 for (j=0; j<d-1; j++, k++)
3667 gcoeff(M,k+1,k) = gen_1;
3668 for (j=0; j<d; j++)
3669 gcoeff(M,k-j,k) = gneg(gel(P, 1+d-j));
3670 }
3671 return M;
3672 }
3673
3674 static GEN
build_frobeniusbc(GEN V,long n)3675 build_frobeniusbc(GEN V, long n)
3676 {
3677 long i,j,k,l;
3678 GEN z;
3679 long m=lg(V)-1;
3680 GEN M = cgetg(n+1, t_MAT);
3681 for(i=1; i<=n; i++)
3682 gel(M,i) = zerocol(n);
3683 z = monomial(gen_m1, 1, 0);
3684 for (k=1,l=1+m,i=1;i<=m;i++,k++)
3685 {
3686 GEN P = gel(V,i);
3687 long d = degpol(P);
3688 if (d <= 0) continue;
3689 if (l+d-2 > n) pari_err(talker, "accuracy lost in matfrobenius");
3690 gcoeff(M,k,i) = gen_1;
3691 for (j=1; j<d; j++,k++,l++)
3692 {
3693 gcoeff(M,k,l) = z;
3694 gcoeff(M,k+1,l) = gen_1;
3695 }
3696 }
3697 return M;
3698 }
3699
3700 static GEN
build_basischange(GEN N,GEN U)3701 build_basischange(GEN N, GEN U)
3702 {
3703 long n = lg(N);
3704 long i, j;
3705 GEN p2 = cgetg(n, t_MAT);
3706 for (j = 1; j < n; ++j)
3707 {
3708 pari_sp btop=avma;
3709 GEN p3 = gen_0;
3710 for (i = 1; i < n; ++i)
3711 p3 = gadd(p3, (GEN) gsubst(gcoeff(U, i, j), 0, N)[i]);
3712 gel(p2,j) = gerepileupto(btop, p3);
3713 }
3714 return p2;
3715 }
3716
3717 GEN
matfrobenius(GEN M,long flag,long v)3718 matfrobenius(GEN M, long flag, long v)
3719 {
3720 pari_sp ltop=avma;
3721 long n;
3722 GEN D, A, N, B, R, M_x;
3723 if (typ(M)!=t_MAT) pari_err(typeer,"matfrobenius");
3724 if (v<0) v=0;
3725 if (gvar(M)<=v)
3726 pari_err(talker,"variable must have higher priority in matfrobenius");
3727 n = lg(M)-1;
3728 if (n && lg(M[1])!=n+1) pari_err(mattype1,"matfrobenius");
3729 M_x = gaddmat(monomial(gen_m1, 1, v), M);
3730 if (flag<2)
3731 {
3732 D = matsnf0(M_x,6);
3733 if (flag != 1) D = Frobeniusform(D, n);
3734 return gerepileupto(ltop, D);
3735 }
3736 if (flag>2) pari_err(flagerr,"matfrobenius");
3737 A = matsnf0(M_x,3);
3738 D = smithclean(mattodiagonal_i(gel(A,3)));
3739 N = Frobeniusform(D, n);
3740 B = build_frobeniusbc(D, n);
3741 R = build_basischange(N, gmul(B,gel(A,1)));
3742 return gerepilecopy(ltop, mkvec2(N,R));
3743 }
3744