Line data Source code
1 : /* Copyright (C) 2000 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation; either version 2 of the License, or (at your option) any later
8 : version. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 : #include "pari.h"
15 : #include "paripriv.h"
16 :
17 : #define DEBUGLEVEL DEBUGLEVEL_mathnf
18 :
19 : #define dbg_printf(lvl) if (DEBUGLEVEL >= (lvl) + 3) err_printf
20 :
21 : /********************************************************************/
22 : /** **/
23 : /** BLACK BOX HERMITE RINGS AND HOWELL NORMAL FORM **/
24 : /** contributed by Aurel Page (2017) **/
25 : /** **/
26 : /********************************************************************/
27 :
28 : /*
29 : bb_hermite R:
30 : - add(a,b): a+b
31 : - neg(a): -a
32 : - mul(a,b): a*b
33 : - extgcd(a,b,&small): [d,U] with d in R and U in GL_2(R) such that [0;d] = [a;b]*U.
34 : set small==1 to assert that U is a 'small' operation (no red needed).
35 : - rann(a): b in R such that b*R = {x in R | a*x==0}
36 : - lquo(a,b,&r): q in R such that r=a-b*q is a canonical representative
37 : of the image of a in R/b*R. The canonical lift of 0 must be 0.
38 : - unit(a): u unit in R^* such that a*u is a canonical generator of the ideal a*R
39 : - equal0(a): a==0?
40 : - equal1(a): a==1?
41 : - s(n): image of the small integer n in R
42 : - red(a): unique representative of a as an element of R
43 :
44 : op encoding of elementary operations:
45 : - t_VECSMALL: the corresponding permutation (vecpermute)
46 : - [Vecsmall([i,j])]: the transposition Ci <-> Cj
47 : - [Vecsmall([i]),u], u in R^*: Ci <- Ci*u
48 : - [Vecsmall([i,j]),a], a in R: Ci <- Ci + Cj*a
49 : - [Vecsmall([i,j,0]),U], U in GL_2(R): (Ci|Cj) <- (Ci|Cj)*U
50 : */
51 :
52 : struct bb_hermite
53 : {
54 : GEN (*add)(void*, GEN, GEN);
55 : GEN (*neg)(void*, GEN);
56 : GEN (*mul)(void*, GEN, GEN);
57 : GEN (*extgcd)(void*, GEN, GEN, int*);
58 : GEN (*rann)(void*, GEN);
59 : GEN (*lquo)(void*, GEN, GEN, GEN*);
60 : GEN (*unit)(void*, GEN);
61 : int (*equal0)(GEN);
62 : int (*equal1)(GEN);
63 : GEN (*s)(void*, long);
64 : GEN (*red)(void*, GEN);
65 : };
66 :
67 : static GEN
68 25670382 : _Fp_add(void *data, GEN x, GEN y) { (void) data; return addii(x,y); }
69 :
70 : static GEN
71 4301241 : _Fp_neg(void *data, GEN x) { (void) data; return negi(x); }
72 :
73 : static GEN
74 36842599 : _Fp_mul(void *data, GEN x, GEN y) { (void) data; return mulii(x,y); }
75 :
76 : static GEN
77 2161949 : _Fp_rann(void *data, GEN x)
78 : {
79 2161949 : GEN d, N = (GEN)data;
80 2161949 : if (!signe(x)) return gen_1;
81 2016265 : d = gcdii(x,N);
82 2016058 : return modii(diviiexact(N,d),N);
83 : }
84 :
85 : static GEN
86 2212378 : _Fp_lquo(void *data, GEN x, GEN y, GEN* r) { (void) data; return truedvmdii(x,y,r); }
87 :
88 : /* D=MN, p|M => !p|a, p|N => p|a, return M */
89 : static GEN
90 0 : Z_split(GEN D, GEN a)
91 : {
92 : long i, n;
93 : GEN N;
94 0 : n = expi(D);
95 0 : n = n<2 ? 1 : expu(n)+1;
96 0 : for (i=1;i<=n;i++)
97 0 : a = Fp_sqr(a,D);
98 0 : N = gcdii(a,D);
99 0 : return diviiexact(D,N);
100 : }
101 :
102 : /* c s.t. gcd(a+cb,N) = gcd(a,b,N) without factoring */
103 : static GEN
104 0 : Z_stab(GEN a, GEN b, GEN N)
105 : {
106 : GEN g, a2, N2;
107 0 : g = gcdii(a,b);
108 0 : g = gcdii(g,N);
109 0 : N2 = diviiexact(N,g);
110 0 : a2 = diviiexact(a,g);
111 0 : return Z_split(N2,a2);
112 : }
113 :
114 : static GEN
115 6182870 : _Fp_unit(void *data, GEN x)
116 : {
117 6182870 : GEN g,s,v,d,N=(GEN)data,N2;
118 : long i;
119 6182870 : if (!signe(x)) return NULL;
120 5756304 : g = bezout(x,N,&s,&v);
121 5757072 : if (equali1(g) || equali1(gcdii(s,N))) return mkvec2(g,s);
122 44672 : N2 = diviiexact(N,g);
123 61164 : for (i=0; i<5; i++)
124 : {
125 61164 : s = addii(s,N2);
126 61164 : if (equali1(gcdii(s,N))) return mkvec2(g,s);
127 : }
128 0 : d = Z_stab(s,N2,N);
129 0 : d = mulii(d,N2);
130 0 : v = Fp_add(s,d,N);
131 0 : if (equali1(v)) return NULL;
132 0 : return mkvec2(g,v);
133 : }
134 :
135 : static GEN
136 3134463 : _Fp_extgcd(void *data, GEN x, GEN y, int* smallop)
137 : {
138 : GEN d,u,v,m;
139 3134463 : if (equali1(y))
140 : {
141 548821 : *smallop = 1;
142 548821 : return mkvec2(y,mkmat2(
143 : mkcol2(gen_1,Fp_neg(x,(GEN)data)),
144 : mkcol2(gen_0,gen_1)));
145 : }
146 2585642 : *smallop = 0;
147 2585642 : d = bezout(x,y,&u,&v);
148 2585640 : if (!signe(d)) return mkvec2(d,matid(2));
149 2585640 : m = cgetg(3,t_MAT);
150 2585643 : m = mkmat2(
151 : mkcol2(diviiexact(y,d),negi(diviiexact(x,d))),
152 : mkcol2(u,v));
153 2585648 : return mkvec2(d,m);
154 : }
155 :
156 : static int
157 89146023 : _Fp_equal0(GEN x) { return !signe(x); }
158 :
159 : static int
160 25914292 : _Fp_equal1(GEN x) { return equali1(x); }
161 :
162 : static GEN
163 14458316 : _Fp_s(void *data, long x)
164 : {
165 14458316 : if (!x) return gen_0;
166 1893433 : if (x==1) return gen_1;
167 0 : return modsi(x,(GEN)data);
168 : }
169 :
170 : static GEN
171 42004373 : _Fp_red(void *data, GEN x) { return Fp_red(x, (GEN)data); }
172 :
173 : /* p not necessarily prime */
174 : static const struct bb_hermite Fp_hermite=
175 : {_Fp_add,_Fp_neg,_Fp_mul,_Fp_extgcd,_Fp_rann,_Fp_lquo,_Fp_unit,_Fp_equal0,_Fp_equal1,_Fp_s,_Fp_red};
176 :
177 : static const struct bb_hermite*
178 840237 : get_Fp_hermite(void **data, GEN p)
179 : {
180 840237 : *data = (void*)p; return &Fp_hermite;
181 : }
182 :
183 : static void
184 17002635 : gen_redcol(GEN C, long lim, void* data, const struct bb_hermite *R)
185 : {
186 : long i;
187 61097602 : for (i=1; i<=lim; i++)
188 44100736 : if (!R->equal0(gel(C,i)))
189 23523366 : gel(C,i) = R->red(data, gel(C,i));
190 16996866 : }
191 :
192 : static GEN
193 : /* return NULL if a==0 */
194 : /* assume C*a is zero after lim */
195 17428477 : gen_rightmulcol(GEN C, GEN a, long lim, int fillzeros, void* data, const struct bb_hermite *R)
196 : {
197 : GEN Ca,zero;
198 : long i;
199 17428477 : if (R->equal1(a)) return C;
200 10440390 : if (R->equal0(a)) return NULL;
201 7519668 : Ca = cgetg(lg(C),t_COL);
202 27490988 : for (i=1; i<=lim; i++)
203 19971359 : gel(Ca,i) = R->mul(data, gel(C,i), a);
204 7519629 : if (fillzeros && lim+1 < lg(C))
205 : {
206 5757257 : zero = R->s(data,0);
207 24711737 : for (i=lim+1; i<lg(C); i++)
208 18954474 : gel(Ca,i) = zero;
209 : }
210 7519635 : return Ca;
211 : }
212 :
213 : static void
214 : /* C1 <- C1 + C2 */
215 : /* assume C2[i]==0 for i>lim */
216 5023536 : gen_addcol(GEN C1, GEN C2, long lim, void* data, const struct bb_hermite *R)
217 : {
218 : long i;
219 18068564 : for (i=1; i<=lim; i++)
220 13045030 : gel(C1,i) = R->add(data, gel(C1,i), gel(C2,i));
221 5023534 : }
222 :
223 : static void
224 : /* H[,i] <- H[,i] + C*a */
225 : /* assume C is zero after lim */
226 3823947 : gen_addrightmul(GEN H, GEN C, GEN a, long i, long lim, void* data, const struct bb_hermite *R)
227 : {
228 : GEN Ca;
229 3823947 : if (R->equal0(a)) return;
230 1386552 : Ca = gen_rightmulcol(C, a, lim, 0, data, R);
231 1386558 : gen_addcol(gel(H,i), Ca, lim, data, R);
232 : }
233 :
234 : static GEN
235 2933511 : gen_zerocol(long n, void* data, const struct bb_hermite *R)
236 : {
237 2933511 : GEN C = cgetg(n+1,t_COL), zero = R->s(data, 0);
238 : long i;
239 13942027 : for (i=1; i<=n; i++) gel(C,i) = zero;
240 2933499 : return C;
241 : }
242 :
243 : static GEN
244 1552478 : gen_zeromat(long m, long n, void* data, const struct bb_hermite *R)
245 : {
246 1552478 : GEN M = cgetg(n+1,t_MAT);
247 : long i;
248 4343509 : for (i=1; i<=n; i++) gel(M,i) = gen_zerocol(m, data, R);
249 1552440 : return M;
250 : }
251 :
252 : static GEN
253 364161 : gen_colei(long n, long i, void* data, const struct bb_hermite *R)
254 : {
255 364161 : GEN C = cgetg(n+1,t_COL), zero = R->s(data, 0);
256 : long j;
257 2490110 : for (j=1; j<=n; j++)
258 2125949 : if (i!=j) gel(C,j) = zero;
259 364161 : else gel(C,j) = R->s(data,1);
260 364161 : return C;
261 : }
262 :
263 : static GEN
264 56084 : gen_matid_hermite(long n, void* data, const struct bb_hermite *R)
265 : {
266 56084 : GEN M = cgetg(n+1,t_MAT);
267 : long i;
268 196196 : for (i=1; i<=n; i++) gel(M,i) = gen_colei(n, i, data, R);
269 56084 : return M;
270 : }
271 :
272 : static GEN
273 2571494 : gen_matmul_hermite(GEN A, GEN B, void* data, const struct bb_hermite *R)
274 : {
275 2571494 : GEN M,sum,prod,zero = R->s(data,0);
276 : long a,b,c,c2,i,j,k;
277 2571493 : RgM_dimensions(A,&a,&c);
278 2571495 : RgM_dimensions(B,&c2,&b);
279 2571494 : if (c!=c2) pari_err_DIM("gen_matmul_hermite");
280 2571494 : M = cgetg(b+1,t_MAT);
281 5475080 : for (j=1; j<=b; j++)
282 : {
283 2903660 : gel(M,j) = cgetg(a+1,t_COL);
284 8169024 : for (i=1; i<=a; i++)
285 : {
286 5265474 : sum = zero;
287 16835880 : for (k=1; k<=c; k++)
288 : {
289 11570514 : prod = R->mul(data, gcoeff(A,i,k), gcoeff(B,k,j));
290 11570252 : sum = R->add(data, sum, prod);
291 : }
292 5265366 : gcoeff(M,i,j) = sum;
293 : }
294 2903550 : gen_redcol(gel(M,j), a, data, R);
295 : }
296 2571420 : return M;
297 : }
298 :
299 : static void
300 : /* U = [u1,u2]~, C <- A*u1 + B*u2 */
301 : /* assume both A, B and C are zero after lim */
302 6446609 : gen_rightlincomb(GEN A, GEN B, GEN U, GEN *C, long lim, void* data, const struct bb_hermite *R)
303 : {
304 : GEN Au1, Bu2;
305 6446609 : Au1 = gen_rightmulcol(A, gel(U,1), lim, 1, data, R);
306 6446609 : Bu2 = gen_rightmulcol(B, gel(U,2), lim, 1, data, R);
307 6446604 : if (!Au1 && !Bu2) { *C = gen_zerocol(lg(A)-1, data, R); return; }
308 6446604 : if (!Au1) { *C = Bu2; return; }
309 4010798 : if (!Bu2) { *C = Au1; return; }
310 3636849 : gen_addcol(Au1, Bu2, lim, data, R);
311 3636852 : *C = Au1;
312 : }
313 :
314 : static void
315 : /* (H[,i] | H[,j]) <- (H[,i] | H[,j]) * U */
316 : /* assume both columns are zero after lim */
317 3223308 : gen_elem(GEN H, GEN U, long i, long j, long lim, void* data, const struct bb_hermite *R)
318 : {
319 : GEN Hi, Hj;
320 3223308 : Hi = shallowcopy(gel(H,i));
321 3223308 : Hj = shallowcopy(gel(H,j));
322 3223307 : gen_rightlincomb(Hi, Hj, gel(U,1), &gel(H,i), lim, data, R);
323 3223301 : gen_rightlincomb(Hi, Hj, gel(U,2), &gel(H,j), lim, data, R);
324 3223304 : }
325 :
326 : static int
327 : /* assume C is zero after lim */
328 4065401 : gen_is_zerocol(GEN C, long lim, void* data, const struct bb_hermite *R)
329 : {
330 : long i;
331 : (void) data;
332 7338664 : for (i=1; i<=lim; i++)
333 5758067 : if (!R->equal0(gel(C,i))) return 0;
334 1580597 : return 1;
335 : }
336 :
337 : /* The mkop* functions return NULL if the corresponding operation is the identity */
338 :
339 : static GEN
340 : /* Ci <- Ci + Cj*a */
341 2750194 : mkoptransv(long i, long j, GEN a, void* data, const struct bb_hermite *R)
342 : {
343 2750194 : a = R->red(data,a);
344 2750167 : if (R->equal0(a)) return NULL;
345 1065422 : return mkvec2(mkvecsmall2(i,j),a);
346 : }
347 :
348 : static GEN
349 : /* (Ci|Cj) <- (Ci|Cj)*U */
350 912877 : mkopU(long i, long j, GEN U, void* data, const struct bb_hermite *R)
351 : {
352 912877 : if (R->equal1(gcoeff(U,1,1)) && R->equal0(gcoeff(U,1,2))
353 391929 : && R->equal1(gcoeff(U,2,2))) return mkoptransv(i,j,gcoeff(U,2,1),data,R);
354 520948 : return mkvec2(mkvecsmall3(i,j,0),U);
355 : }
356 :
357 : static GEN
358 : /* Ci <- Ci*u */
359 1871358 : mkopmul(long i, GEN u, const struct bb_hermite *R)
360 : {
361 1871358 : if (R->equal1(u)) return NULL;
362 189905 : return mkvec2(mkvecsmall(i),u);
363 : }
364 :
365 : static GEN
366 : /* Ci <-> Cj */
367 45843 : mkopswap(long i, long j)
368 : {
369 45843 : return mkvec(mkvecsmall2(i,j));
370 : }
371 :
372 : /* M: t_MAT. Apply the operation op to M by right multiplication. */
373 : static void
374 374465 : gen_rightapply(GEN M, GEN op, void* data, const struct bb_hermite *R)
375 : {
376 : GEN M2, ind, X;
377 374465 : long i, j, m = lg(gel(M,1))-1;
378 374465 : switch (typ(op))
379 : {
380 56014 : case t_VECSMALL:
381 56014 : M2 = vecpermute(M,op);
382 266056 : for (i=1; i<lg(M); i++) gel(M,i) = gel(M2,i);
383 56014 : return;
384 318451 : case t_VEC:
385 318451 : ind = gel(op,1);
386 318451 : switch (lg(op))
387 : {
388 16135 : case 2:
389 16135 : swap(gel(M,ind[1]),gel(M,ind[2]));
390 16135 : return;
391 302316 : case 3:
392 302316 : X = gel(op,2);
393 302316 : i = ind[1];
394 302316 : switch (lg(ind))
395 : {
396 37786 : case 2:
397 37786 : gel(M,i) = gen_rightmulcol(gel(M,i), X, m, 0, data, R);
398 37786 : gen_redcol(gel(M,i), m, data, R);
399 37786 : return;
400 175693 : case 3:
401 175693 : gen_addrightmul(M, gel(M,ind[2]), X, i, m, data, R);
402 175693 : gen_redcol(gel(M,i), m, data, R);
403 175693 : return;
404 88837 : case 4:
405 88837 : j = ind[2];
406 88837 : gen_elem(M, X, i, j, m, data, R);
407 88837 : gen_redcol(gel(M,i), m, data, R);
408 88837 : gen_redcol(gel(M,j), m, data, R);
409 88837 : return;
410 : }
411 : }
412 : }
413 : }
414 :
415 : /* C: t_COL. Apply the operation op to C by left multiplication. */
416 : static void
417 11284101 : gen_leftapply(GEN C, GEN op, void* data, const struct bb_hermite *R)
418 : {
419 : GEN C2, ind, X;
420 : long i, j;
421 11284101 : switch (typ(op))
422 : {
423 576065 : case t_VECSMALL:
424 576065 : C2 = vecpermute(C,perm_inv(op));
425 4715305 : for (i=1; i<lg(C); i++) gel(C,i) = gel(C2,i);
426 576065 : return;
427 10708040 : case t_VEC:
428 10708040 : ind = gel(op,1);
429 10708040 : switch (lg(op))
430 : {
431 169134 : case 2:
432 169134 : swap(gel(C,ind[1]),gel(C,ind[2]));
433 169134 : return;
434 10538908 : case 3:
435 10538908 : X = gel(op,2);
436 10538908 : i = ind[1];
437 10538908 : switch (lg(ind))
438 : {
439 689904 : case 2:
440 689904 : gel(C,i) = R->mul(data, X, gel(C,i));
441 689904 : gel(C,i) = R->red(data, gel(C,i));
442 689895 : return;
443 7487183 : case 3:
444 7487183 : j = ind[2];
445 7487183 : if (R->equal0(gel(C,i))) return;
446 1055068 : gel(C,j) = R->add(data, gel(C,j), R->mul(data, X, gel(C,i)));
447 1055065 : return;
448 2361865 : case 4:
449 2361865 : j = ind[2];
450 2361865 : C2 = gen_matmul_hermite(X, mkmat(mkcol2(gel(C,i),gel(C,j))), data, R);
451 2361797 : gel(C,i) = gcoeff(C2,1,1);
452 2361797 : gel(C,j) = gcoeff(C2,2,1);
453 2361797 : return;
454 : }
455 : }
456 : }
457 : }
458 :
459 : /* \prod_i det ops[i]. Only makes sense if R is commutative. */
460 : static GEN
461 42 : gen_detops(GEN ops, void* data, const struct bb_hermite *R)
462 : {
463 42 : GEN d = R->s(data,1);
464 42 : long i, l = lg(ops);
465 231 : for (i = 1; i < l; i++)
466 : {
467 189 : GEN X, op = gel(ops,i);
468 189 : switch (typ(op))
469 : {
470 0 : case t_VECSMALL:
471 0 : if (perm_sign(op) < 0) d = R->neg(data,d);
472 0 : break;
473 189 : case t_VEC:
474 189 : switch (lg(op))
475 : {
476 0 : case 2:
477 0 : d = R->neg(data,d);
478 0 : break;
479 189 : case 3:
480 189 : X = gel(op,2);
481 189 : switch (lg(gel(op,1)))
482 : {
483 0 : case 2:
484 0 : d = R->mul(data, d, X);
485 0 : d = R->red(data, d);
486 0 : break;
487 105 : case 4:
488 : {
489 105 : GEN A = gcoeff(X,1,1), B = gcoeff(X,1,2);
490 105 : GEN C = gcoeff(X,2,1), D = gcoeff(X,2,2);
491 105 : GEN AD = R->mul(data,A,D);
492 105 : GEN BC = R->mul(data,B,C);
493 105 : d = R->mul(data, d, R->add(data, AD, R->neg(data,BC)));
494 105 : d = R->red(data, d);
495 105 : break;
496 : }
497 : }
498 189 : break;
499 : }
500 189 : break;
501 : }
502 189 : }
503 42 : return d;
504 : }
505 :
506 : static int
507 209202 : gen_is_inv(GEN x, void* data, const struct bb_hermite *R)
508 : {
509 209202 : GEN u = R->unit(data, x);
510 209202 : if (!u) return R->equal1(x);
511 71743 : return R->equal1(gel(u,1));
512 : }
513 :
514 : static long
515 193431 : gen_last_inv_diago(GEN A, void* data, const struct bb_hermite *R)
516 : {
517 : long i,m,n,j;
518 193431 : RgM_dimensions(A,&m,&n);
519 223461 : for (i=1,j=n-m+1; i<=m; i++,j++)
520 209202 : if (!gen_is_inv(gcoeff(A,i,j),data,R)) return i-1;
521 14259 : return m;
522 : }
523 :
524 : static GEN
525 : /* remove_zerocols: 0 none, 1 until square, 2 all */
526 : /* early abort: if not right-invertible, abort, return NULL, and set ops to the
527 : * noninvertible pivot */
528 938770 : gen_howell_i(GEN A, long remove_zerocols, long permute_zerocols, long early_abort, long only_triangular, GEN* ops, void *data, const struct bb_hermite *R)
529 : {
530 938770 : pari_sp av = avma, av1;
531 938770 : GEN H,U,piv=gen_0,u,q,a,perm,iszero,C,zero=R->s(data,0),d,g,r,op,one=R->s(data,1);
532 938770 : long m,n,i,j,s,si,i2,si2,nbz,lim,extra,maxop=0,nbop=0,lastinv=0;
533 : int smallop;
534 :
535 938770 : av1 = avma;
536 :
537 938770 : RgM_dimensions(A,&m,&n);
538 938776 : if (early_abort && n<m)
539 : {
540 14000 : if (ops) *ops = zero;
541 14000 : return NULL;
542 : }
543 924776 : if (n<m+1)
544 : {
545 826258 : extra = m+1-n;
546 826258 : H = shallowmatconcat(mkvec2(gen_zeromat(m,extra,data,R),A));
547 : }
548 : else
549 : {
550 98518 : extra = 0;
551 98518 : H = RgM_shallowcopy(A);
552 : }
553 924787 : RgM_dimensions(H,&m,&n);
554 924783 : s = n-m; /* shift */
555 :
556 924783 : if(ops)
557 : {
558 731352 : maxop = m*n + (m*(m+1))/2 + 1;
559 731352 : *ops = zerovec(maxop); /* filled with placeholders so gerepile can handle it */
560 : }
561 :
562 : /* put in triangular form */
563 3560343 : for (i=m,si=s+m; i>0 && si>extra; i--,si--) /* si = s+i */
564 : {
565 2640551 : if (R->red) gcoeff(H,i,si) = R->red(data, gcoeff(H,i,si));
566 : /* bottom-right diagonal */
567 9192906 : for (j = extra+1; j < si; j++)
568 : {
569 6552768 : if (R->red) gcoeff(H,i,j) = R->red(data, gcoeff(H,i,j));
570 6552800 : if (R->equal0(gcoeff(H,i,j))) continue;
571 2987720 : U = R->extgcd(data, gcoeff(H,i,j), gcoeff(H,i,si), &smallop);
572 2987814 : d = gel(U,1);
573 2987814 : U = gel(U,2);
574 2987814 : if (n>10)
575 : {
576 : /* normalize diagonal coefficient -> faster reductions on this row */
577 1777813 : u = R->unit(data, d);
578 1777813 : if (u)
579 : {
580 1777813 : g = gel(u,1);
581 1777813 : u = gel(u,2);
582 1777813 : gcoeff(U,1,2) = R->mul(data, gcoeff(U,1,2), u);
583 1777813 : gcoeff(U,2,2) = R->mul(data, gcoeff(U,2,2), u);
584 1777813 : d = g;
585 : }
586 : }
587 2987814 : gen_elem(H, U, j, si, i-1, data, R);
588 2987809 : if (ops)
589 : {
590 877071 : op = mkopU(j,si,U,data,R);
591 877074 : if (op) { nbop++; gel(*ops, nbop) = op; }
592 : }
593 2987812 : gcoeff(H,i,j) = zero;
594 2987812 : gcoeff(H,i,si) = d;
595 2987812 : if (R->red && !smallop)
596 : {
597 2448546 : gen_redcol(gel(H,si), i-1, data, R);
598 2448545 : gen_redcol(gel(H,j), i-1, data, R);
599 : }
600 : }
601 :
602 2640138 : if (early_abort)
603 : {
604 1453085 : d = gcoeff(H,i,si);
605 1453085 : u = R->unit(data, d);
606 1453439 : if (u) d = gel(u,1);
607 1453439 : if (!R->equal1(d))
608 : {
609 4914 : if (ops) *ops = d;
610 4914 : return NULL;
611 : }
612 : }
613 :
614 2635566 : if (gc_needed(av,1))
615 : {
616 14 : if (DEBUGMEM>1) pari_warn(warnmem,"gen_howell[1]. i=%ld",i);
617 14 : gerepileall(av1,ops?2:1,&H,ops);
618 : }
619 : }
620 :
621 919792 : if (!ops)
622 193431 : lastinv = gen_last_inv_diago(H, data, R);
623 :
624 : /* put in reduced Howell form */
625 919792 : if (!only_triangular)
626 : {
627 3693123 : for (i=m,si=s+m; i>0; i--,si--) /* si = s+i */
628 : {
629 : /* normalize diagonal coefficient */
630 2773329 : if (i<=lastinv) /* lastinv>0 => !ops */
631 30030 : gcoeff(H,i,si) = one;
632 : else
633 : {
634 2743299 : u = R->unit(data,gcoeff(H,i,si));
635 2743450 : if (u)
636 : {
637 2454490 : g = gel(u,1);
638 2454490 : u = gel(u,2);
639 2454490 : gel(H,si) = gen_rightmulcol(gel(H,si), u, i-1, 1, data, R);
640 2454484 : gcoeff(H,i,si) = g;
641 2454484 : if (R->red) gen_redcol(gel(H,si), i-1, data, R);
642 2454483 : if (ops)
643 : {
644 1871358 : op = mkopmul(si,u,R);
645 1871356 : if (op) { nbop++; gel(*ops,nbop) = op; }
646 : }
647 : }
648 288960 : else if (R->red) gcoeff(H,i,si) = R->red(data, gcoeff(H,i,si));
649 : }
650 2773471 : piv = gcoeff(H,i,si);
651 :
652 : /* reduce above diagonal */
653 2773471 : if (!R->equal0(piv))
654 : {
655 2484503 : C = gel(H,si);
656 6158909 : for (j=si+1; j<=n; j++)
657 : {
658 3674424 : if (i<=lastinv) /* lastinv>0 => !ops */
659 26187 : gcoeff(H,i,j) = zero;
660 : else
661 : {
662 3648237 : gcoeff(H,i,j) = R->red(data, gcoeff(H,i,j));
663 3648208 : if (R->equal1(piv)) { q = gcoeff(H,i,j); r = zero; }
664 1559439 : else q = R->lquo(data, gcoeff(H,i,j), piv, &r);
665 3648200 : q = R->neg(data,q);
666 3648246 : gen_addrightmul(H, C, q, j, i-1, data, R);
667 3648243 : if (ops)
668 : {
669 2209368 : op = mkoptransv(j,si,q,data,R);
670 2209344 : if (op) { nbop++; gel(*ops,nbop) = op; }
671 : }
672 3648219 : gcoeff(H,i,j) = r;
673 : }
674 : }
675 : }
676 :
677 : /* ensure Howell property */
678 2773444 : if (i>1)
679 : {
680 1853839 : a = R->rann(data, piv);
681 1853671 : if (!R->equal0(a))
682 : {
683 545328 : gel(H,1) = gen_rightmulcol(gel(H,si), a, i-1, 1, data, R);
684 545328 : if (gel(H,1) == gel(H,si)) gel(H,1) = shallowcopy(gel(H,1)); /* in case rightmulcol cheated */
685 545328 : if (ops)
686 : {
687 148897 : op = mkoptransv(1,si,a,data,R);
688 148897 : if (op) { nbop++; gel(*ops,nbop) = op; }
689 : }
690 2185421 : for (i2=i-1,si2=s+i2; i2>0; i2--,si2--)
691 : {
692 1640093 : if (R->red) gcoeff(H,i2,1) = R->red(data, gcoeff(H,i2,1));
693 1640093 : if (R->equal0(gcoeff(H,i2,1))) continue;
694 277354 : if (R->red) gcoeff(H,i2,si2) = R->red(data, gcoeff(H,i2,si2));
695 277354 : if (R->equal0(gcoeff(H,i2,si2)))
696 : {
697 130697 : swap(gel(H,1), gel(H,si2));
698 130697 : if (ops) { nbop++; gel(*ops,nbop) = mkopswap(1,si2); }
699 130697 : continue;
700 : }
701 146657 : U = R->extgcd(data, gcoeff(H,i2,1), gcoeff(H,i2,si2), &smallop);
702 146657 : d = gel(U,1);
703 146657 : U = gel(U,2);
704 146657 : gen_elem(H, U, 1, si2, i2-1, data, R);
705 146657 : if (ops)
706 : {
707 35805 : op = mkopU(1,si2,U,data,R);
708 35805 : if (op) { nbop++; gel(*ops,nbop) = op; }
709 : }
710 146657 : gcoeff(H,i2,1) = zero;
711 146657 : gcoeff(H,i2,si2) = d;
712 146657 : if (R->red && !smallop)
713 : {
714 137102 : gen_redcol(gel(H,si2), i2, data, R);
715 137102 : gen_redcol(gel(H,1), i2-1, data, R);
716 : }
717 : }
718 : }
719 : }
720 :
721 2773296 : if (gc_needed(av,1))
722 : {
723 12 : if (DEBUGMEM>1) pari_warn(warnmem,"gen_howell[2]. i=%ld",i);
724 12 : gerepileall(av1,ops?3:2,&H,&piv,ops);
725 : }
726 : }
727 : }
728 :
729 919772 : if (R->red)
730 4940680 : for (j=1; j<=n; j++)
731 : {
732 4021157 : lim = maxss(0,m-n+j);
733 4021163 : gen_redcol(gel(H,j), lim, data, R);
734 13427341 : for (i=lim+1; i<=m; i++) gcoeff(H,i,j) = zero;
735 : }
736 :
737 : /* put zero columns first */
738 919443 : iszero = cgetg(n+1,t_VECSMALL);
739 :
740 919805 : nbz = 0;
741 4941264 : for (i=1; i<=n; i++)
742 : {
743 4021416 : iszero[i] = gen_is_zerocol(gel(H,i), maxss(0,m-n+i), data, R);
744 4021459 : if (iszero[i]) nbz++;
745 : }
746 :
747 919848 : j = 1;
748 919848 : if (permute_zerocols)
749 : {
750 154665 : perm = cgetg(n+1, t_VECSMALL);
751 897176 : for (i=1; i<=n; i++)
752 742511 : if (iszero[i])
753 : {
754 313229 : perm[j] = i;
755 313229 : j++;
756 : }
757 : }
758 765183 : else perm = cgetg(n-nbz+1, t_VECSMALL);
759 4941473 : for (i=1; i<=n; i++)
760 4021608 : if (!iszero[i])
761 : {
762 2484825 : perm[j] = i;
763 2484825 : j++;
764 : }
765 :
766 919865 : if (permute_zerocols || remove_zerocols==2) H = vecpermute(H, perm);
767 919867 : if (permute_zerocols && remove_zerocols==2) H = vecslice(H, nbz+1, n);
768 919867 : if (remove_zerocols==1) H = vecslice(H, s+1, n);
769 919867 : if (permute_zerocols && ops) { nbop++; gel(*ops,nbop) = perm; }
770 :
771 919867 : if (ops) { setlg(*ops, nbop+1); } /* should have nbop <= maxop */
772 :
773 919863 : return H;
774 : }
775 :
776 : static GEN
777 94948 : gen_howell(GEN A, long remove_zerocols, long permute_zerocols, long early_abort, long only_triangular, GEN* ops, void *data, const struct bb_hermite *R)
778 : {
779 94948 : pari_sp av = avma;
780 94948 : GEN H = gen_howell_i(A, remove_zerocols, permute_zerocols, early_abort, only_triangular, ops, data, R);
781 94948 : return gc_all(av, ops?2:1, &H, ops);
782 : }
783 :
784 : static GEN
785 150962 : gen_matimage(GEN A, GEN* U, void *data, const struct bb_hermite *R)
786 : {
787 : GEN ops, H;
788 150962 : if (U)
789 : {
790 56014 : pari_sp av = avma, av1;
791 : long m, n, i, r, n2, pergc;
792 56014 : RgM_dimensions(A,&m,&n);
793 56014 : H = gen_howell_i(A, 2, 1, 0, 0, &ops, data, R);
794 56014 : av1 = avma;
795 56014 : r = lg(H)-1;
796 56014 : *U = shallowmatconcat(mkvec2(gen_zeromat(n, maxss(0,m-n+1), data, R), gen_matid_hermite(n, data, R)));
797 56014 : n2 = lg(*U)-1;
798 56014 : pergc = maxss(m,n);
799 430479 : for (i=1; i<lg(ops); i++)
800 : {
801 374465 : gen_rightapply(*U, gel(ops,i), data, R);
802 374465 : if (!(i%pergc) && gc_needed(av,1))
803 : {
804 0 : if (DEBUGMEM>1) pari_warn(warnmem,"gen_matimage. i=%ld",i);
805 0 : gerepileall(av1,1,U);
806 : }
807 : }
808 56014 : if (r<n2) *U = vecslice(*U, n2-r+1, n2);
809 56014 : return gc_all(av, 2, &H, U);
810 : }
811 94948 : return gen_howell(A, 2, 0, 0, 0, NULL, data, R);
812 : }
813 :
814 : static GEN
815 : /* H in true Howell form: no zero columns */
816 98483 : gen_kernel_howell(GEN H, void *data, const struct bb_hermite *R)
817 : {
818 : GEN K, piv, FK;
819 : long m, n, j, j2, i;
820 98483 : RgM_dimensions(H,&m,&n);
821 98483 : K = gen_zeromat(n, n, data, R);
822 406595 : for (j=n,i=m; j>0; j--)
823 : {
824 521619 : while (R->equal0(gcoeff(H,i,j))) i--;
825 308112 : piv = gcoeff(H,i,j);
826 308112 : if (R->equal0(piv)) continue;
827 308112 : gcoeff(K,j,j) = R->rann(data, piv);
828 308112 : if (j<n)
829 : {
830 209629 : FK = gen_matmul_hermite(matslice(H,i,i,j+1,n), matslice(K, j+1, n, j+1, n), data, R);
831 751422 : for (j2=j+1; j2<=n; j2++)
832 541793 : gcoeff(K,j,j2) = R->neg(data, R->lquo(data, gcoeff(FK,1,j2-j), piv, NULL));
833 : /* remainder has to be zero */
834 : }
835 : }
836 98483 : return K;
837 : }
838 :
839 : static GEN
840 : /* (H,ops) Howell form of A, n = number of columns of A, return a kernel of A */
841 98553 : gen_kernel_from_howell(GEN H, GEN ops, long n, void *data, const struct bb_hermite *R)
842 : {
843 : pari_sp av;
844 : GEN K, KH, zC;
845 : long m, r, n2, nbz, i, o, extra, j;
846 98553 : RgM_dimensions(H,&m,&r);
847 98553 : if (!r) return gen_matid_hermite(n, data, R); /* zerology: what if 0==1 in R? */
848 98483 : n2 = maxss(n,m+1);
849 98483 : extra = n2-n;
850 98483 : nbz = n2-r;
851 : /* compute kernel of augmented matrix */
852 98483 : KH = gen_kernel_howell(H, data, R);
853 98483 : zC = gen_zerocol(nbz, data, R);
854 98483 : K = cgetg(nbz+r+1, t_MAT);
855 322532 : for (i=1; i<=nbz; i++)
856 224049 : gel(K,i) = gen_colei(nbz+r, i, data, R);
857 406595 : for (i=1; i<=r; i++)
858 308112 : gel(K,nbz+i) = shallowconcat(zC, gel(KH,i));
859 630644 : for (i=1; i<lg(K); i++)
860 : {
861 532161 : av = avma;
862 9775220 : for (o=lg(ops)-1; o>0; o--)
863 9243059 : gen_leftapply(gel(K,i), gel(ops,o), data, R);
864 532161 : gen_redcol(gel(K,i), nbz+r, data, R);
865 532161 : gerepileall(av, 1, &gel(K,i));
866 : }
867 : /* deduce kernel of original matrix */
868 98483 : K = rowpermute(K, cyclic_perm(n2,extra));
869 98483 : K = gen_howell_i(K, 2, 0, 0, 0, NULL, data, R);
870 225029 : for (j=lg(K)-1, i=n2; j>0; j--)
871 : {
872 299453 : while (R->equal0(gcoeff(K,i,j))) i--;
873 196546 : if (i<=n) return matslice(K, 1, n, 1, j);
874 : }
875 28483 : return cgetg(1,t_MAT);
876 : }
877 :
878 : /* not GC-clean */
879 : static GEN
880 54698 : gen_kernel(GEN A, GEN* im, void *data, const struct bb_hermite *R)
881 : {
882 54698 : pari_sp av = avma;
883 54698 : long n = lg(A)-1;
884 : GEN H, ops, K;
885 54698 : H = gen_howell_i(A, 2, 1, 0, 0, &ops, data, R);
886 54698 : gerepileall(av,2,&H,&ops);
887 54698 : K = gen_kernel_from_howell(H, ops, n, data, R);
888 54698 : if (im) *im = H;
889 54698 : return K;
890 : }
891 :
892 : /* right inverse, not GC-clean */
893 : static GEN
894 590620 : gen_inv(GEN A, void* data, const struct bb_hermite *R)
895 : {
896 : pari_sp av;
897 590620 : GEN ops, H, U, un=R->s(data,1);
898 : long m,n,j,o,n2;
899 590619 : RgM_dimensions(A,&m,&n);
900 590622 : av = avma;
901 590622 : H = gen_howell_i(A, 0, 0, 1, 0, &ops, data, R);
902 590640 : if (!H) pari_err_INV("gen_inv", ops);
903 571726 : n2 = lg(H)-1;
904 571726 : ops = gerepilecopy(av,ops); /* get rid of H */
905 571746 : U = gen_zeromat(n2, m, data, R);
906 2014003 : for (j=1; j<=m; j++)
907 1442261 : gcoeff(U,j+n2-m,j) = un;
908 2013992 : for (j=1; j<=m; j++)
909 : {
910 1442246 : av = avma;
911 3084961 : for (o=lg(ops)-1; o>0; o--)
912 1642820 : gen_leftapply(gel(U,j), gel(ops,o), data, R);
913 1442141 : gen_redcol(gel(U,j), n2, data, R);
914 1441679 : gerepileall(av, 1, &gel(U,j));
915 : }
916 571746 : if (n2>n) U = rowslice(U, n2-n+1, n2);
917 571725 : return U;
918 : }
919 :
920 : /*
921 : H true Howell form (no zero columns).
922 : Compute Z = Y - HX canonical representative of R^m mod H(R^n)
923 : */
924 : static GEN
925 43953 : gen_reduce_mod_howell(GEN H, GEN Y, GEN *X, void* data, const struct bb_hermite *R)
926 : {
927 43953 : pari_sp av = avma;
928 : long m,n,i,j;
929 : GEN Z, q, r, C;
930 43953 : RgM_dimensions(H,&m,&n);
931 43953 : if (X) *X = gen_zerocol(n,data,R);
932 43953 : Z = shallowcopy(Y);
933 43953 : i = m;
934 155099 : for (j=n; j>0; j--)
935 : {
936 180348 : while (R->equal0(gcoeff(H,i,j))) i--;
937 111146 : q = R->lquo(data, gel(Z,i), gcoeff(H,i,j), &r);
938 111146 : gel(Z,i) = r;
939 111146 : C = gen_rightmulcol(gel(H,j), R->neg(data,q), i, 0, data, R);
940 111146 : if (C) gen_addcol(Z, C, i-1, data, R);
941 111146 : if (X) gel(*X,j) = q;
942 : }
943 43953 : gen_redcol(Z, lg(Z)-1, data, R);
944 43953 : return gc_all(av, X?2:1, &Z, X);
945 : }
946 :
947 : /* H: Howell form of A */
948 : /* (m,n): dimensions of original matrix A */
949 : /* return NULL if no solution */
950 : static GEN
951 43953 : gen_solve_from_howell(GEN H, GEN ops, long m, long n, GEN Y, void* data, const struct bb_hermite *R)
952 : {
953 43953 : pari_sp av = avma;
954 : GEN Z, X;
955 : long n2, mH, nH, i;
956 43953 : RgM_dimensions(H,&mH,&nH);
957 43953 : n2 = maxss(n,m+1);
958 43953 : Z = gen_reduce_mod_howell(H, Y, &X, data, R);
959 43953 : if (!gen_is_zerocol(Z,m,data,R)) { set_avma(av); return NULL; }
960 :
961 43904 : X = shallowconcat(zerocol(n2-nH),X);
962 442113 : for (i=lg(ops)-1; i>0; i--) gen_leftapply(X, gel(ops,i), data, R);
963 43904 : X = vecslice(X, n2-n+1, n2);
964 43904 : gen_redcol(X, n, data, R);
965 43904 : return gerepilecopy(av, X);
966 : }
967 :
968 : /* return NULL if no solution, not GC-clean */
969 : static GEN
970 43953 : gen_solve(GEN A, GEN Y, GEN* K, void* data, const struct bb_hermite *R)
971 : {
972 : GEN H, ops, X;
973 : long m,n;
974 :
975 43953 : RgM_dimensions(A,&m,&n);
976 43953 : if (!n) m = lg(Y)-1;
977 43953 : H = gen_howell_i(A, 2, 1, 0, 0, &ops, data, R);
978 43953 : X = gen_solve_from_howell(H, ops, m, n, Y, data, R);
979 43953 : if (!X) return NULL;
980 43904 : if (K) *K = gen_kernel_from_howell(H, ops, n, data, R);
981 43904 : return X;
982 : }
983 :
984 : GEN
985 150997 : matimagemod(GEN A, GEN d, GEN* U)
986 : {
987 : void *data;
988 : const struct bb_hermite* R;
989 150997 : if (typ(A)!=t_MAT || !RgM_is_ZM(A)) pari_err_TYPE("matimagemod", A);
990 150983 : if (typ(d)!=t_INT) pari_err_TYPE("matimagemod", d);
991 150976 : if (signe(d)<=0) pari_err_DOMAIN("matimagemod", "d", "<=", gen_0, d);
992 150969 : if (equali1(d)) return cgetg(1,t_MAT);
993 150962 : R = get_Fp_hermite(&data, d);
994 150962 : return gen_matimage(A, U, data, R);
995 : }
996 :
997 : /* for testing purpose */
998 : /*
999 : GEN
1000 : ZM_hnfmodid2(GEN A, GEN d)
1001 : {
1002 : pari_sp av = avma;
1003 : void *data;
1004 : long i;
1005 : const struct bb_hermite* R = get_Fp_hermite(&data, d);
1006 : GEN H;
1007 : if (typ(A)!=t_MAT || !RgM_is_ZM(A)) pari_err_TYPE("ZM_hnfmodid2", A);
1008 : if (typ(d)!=t_INT) pari_err_TYPE("ZM_hnfmodid2", d);
1009 : H = gen_howell_i(A, 1, 0, 0, 0, NULL, data, R);
1010 : for (i=1; i<lg(H); i++)
1011 : if (!signe(gcoeff(H,i,i))) gcoeff(H,i,i) = d;
1012 : return gerepilecopy(av,H);
1013 : }
1014 : */
1015 :
1016 : GEN
1017 84 : matdetmod(GEN A, GEN d)
1018 : {
1019 84 : pari_sp av = avma;
1020 : void *data;
1021 : const struct bb_hermite* R;
1022 84 : long n = lg(A)-1, i;
1023 : GEN D, H, ops;
1024 84 : if (typ(A)!=t_MAT || !RgM_is_ZM(A)) pari_err_TYPE("matdetmod", A);
1025 70 : if (typ(d)!=t_INT) pari_err_TYPE("matdetmod", d);
1026 70 : if (signe(d)<=0) pari_err_DOMAIN("matdetmod", "d", "<=", gen_0, d);
1027 63 : if (!n) return equali1(d) ? gen_0 : gen_1;
1028 56 : if (n != nbrows(A)) pari_err_DIM("matdetmod");
1029 49 : if (equali1(d)) return gen_0;
1030 42 : R = get_Fp_hermite(&data, d);
1031 42 : H = gen_howell_i(A, 1, 0, 0, 1, &ops, data, R);
1032 42 : D = gen_detops(ops, data, R);
1033 42 : D = Fp_inv(D, d);
1034 203 : for (i = 1; i <= n; i++) D = Fp_mul(D, gcoeff(H,i,i), d);
1035 42 : return gerepileuptoint(av, D);
1036 : }
1037 :
1038 : GEN
1039 54733 : matkermod(GEN A, GEN d, GEN* im)
1040 : {
1041 54733 : pari_sp av = avma;
1042 : void *data;
1043 : const struct bb_hermite* R;
1044 : long m,n;
1045 : GEN K;
1046 54733 : if (typ(A)!=t_MAT || !RgM_is_ZM(A)) pari_err_TYPE("matkermod", A);
1047 54719 : if (typ(d)!=t_INT) pari_err_TYPE("matkermod", d);
1048 54712 : if (signe(d)<=0) pari_err_DOMAIN("makermod", "d", "<=", gen_0, d);
1049 54705 : if (equali1(d)) return cgetg(1,t_MAT);
1050 54698 : R = get_Fp_hermite(&data, d);
1051 54698 : RgM_dimensions(A,&m,&n);
1052 54698 : if (!im && m>2*n) /* TODO tune treshold */
1053 3577 : A = shallowtrans(matimagemod(shallowtrans(A),d,NULL));
1054 54698 : K = gen_kernel(A, im, data, R); return gc_all(av,im?2:1,&K,im);
1055 : }
1056 :
1057 : /* left inverse */
1058 : GEN
1059 676687 : matinvmod(GEN A, GEN d)
1060 : {
1061 676687 : pari_sp av = avma;
1062 : void *data;
1063 : const struct bb_hermite* R;
1064 : GEN U;
1065 676687 : if (typ(A)!=t_MAT || !RgM_is_ZM(A)) pari_err_TYPE("matinvmod", A);
1066 676726 : if (typ(d)!=t_INT) pari_err_TYPE("matinvmod", d);
1067 676719 : if (signe(d)<=0) pari_err_DOMAIN("matinvmod", "d", "<=", gen_0, d);
1068 676712 : if (equali1(d)) {
1069 : long m,n;
1070 86156 : RgM_dimensions(A,&m,&n);
1071 86156 : if (m<n) pari_err_INV("matinvmod",A);
1072 86149 : return zeromatcopy(n,m);
1073 : }
1074 590558 : R = get_Fp_hermite(&data, d);
1075 590582 : U = gen_inv(shallowtrans(A), data, R);
1076 571726 : return gerepilecopy(av, shallowtrans(U));
1077 : }
1078 :
1079 : /* assume all D[i]>0, not GC-clean */
1080 : static GEN
1081 43953 : matsolvemod_finite(GEN M, GEN D, GEN Y, long flag)
1082 : {
1083 : void *data;
1084 : const struct bb_hermite* R;
1085 : GEN X, K, d;
1086 : long m, n;
1087 :
1088 43953 : RgM_dimensions(M,&m,&n);
1089 43953 : if (typ(D)==t_COL)
1090 : { /* create extra variables for the D[i] */
1091 : GEN MD;
1092 84 : long i, i2, extra = 0;
1093 84 : d = gen_1;
1094 231 : for (i=1; i<lg(D); i++) d = lcmii(d,gel(D,i));
1095 231 : for (i=1; i<lg(D); i++) if (!equalii(gel(D,i),d)) extra++;
1096 84 : MD = cgetg(extra+1,t_MAT);
1097 84 : i2 = 1;
1098 231 : for (i=1; i<lg(D); i++)
1099 147 : if (!equalii(gel(D,i),d))
1100 : {
1101 77 : gel(MD,i2) = Rg_col_ei(gel(D,i),m,i);
1102 77 : i2++;
1103 : }
1104 84 : M = shallowconcat(M,MD);
1105 : }
1106 43869 : else d = D;
1107 :
1108 43953 : R = get_Fp_hermite(&data, d);
1109 43953 : X = gen_solve(M, Y, flag? &K: NULL, data, R);
1110 43953 : if (!X) return gen_0;
1111 43904 : X = vecslice(X,1,n);
1112 :
1113 43904 : if (flag)
1114 : {
1115 43855 : K = rowslice(K,1,n);
1116 43855 : K = hnfmodid(shallowconcat(zerocol(n),K),d);
1117 43855 : X = mkvec2(X,K);
1118 : }
1119 43904 : return X;
1120 : }
1121 :
1122 : /* Return a solution of congruence system sum M[i,j] X_j = Y_i mod D_i
1123 : * If pU != NULL, put in *pU a Z-basis of the homogeneous system.
1124 : * Works for all non negative D_i but inefficient compared to
1125 : * matsolvemod_finite; to be used only when one D_i is 0 */
1126 : static GEN
1127 70 : gaussmoduloall(GEN M, GEN D, GEN Y, GEN *pU)
1128 : {
1129 70 : pari_sp av = avma;
1130 70 : long n, m, j, l, lM = lg(M);
1131 : GEN delta, H, U, u1, u2, x;
1132 :
1133 70 : if (lM == 1)
1134 : {
1135 28 : long lY = 0;
1136 28 : switch(typ(Y))
1137 : {
1138 0 : case t_INT: break;
1139 28 : case t_COL: lY = lg(Y); break;
1140 0 : default: pari_err_TYPE("gaussmodulo",Y);
1141 : }
1142 28 : switch(typ(D))
1143 : {
1144 14 : case t_INT: break;
1145 14 : case t_COL: if (lY && lY != lg(D)) pari_err_DIM("gaussmodulo");
1146 14 : break;
1147 0 : default: pari_err_TYPE("gaussmodulo",D);
1148 : }
1149 28 : if (pU) *pU = cgetg(1, t_MAT);
1150 28 : return cgetg(1,t_COL);
1151 : }
1152 42 : n = nbrows(M);
1153 42 : switch(typ(D))
1154 : {
1155 14 : case t_COL:
1156 14 : if (lg(D)-1!=n) pari_err_DIM("gaussmodulo");
1157 14 : delta = diagonal_shallow(D); break;
1158 28 : case t_INT: delta = scalarmat_shallow(D,n); break;
1159 0 : default: pari_err_TYPE("gaussmodulo",D);
1160 : return NULL; /* LCOV_EXCL_LINE */
1161 : }
1162 42 : switch(typ(Y))
1163 : {
1164 0 : case t_INT: Y = const_col(n, Y); break;
1165 42 : case t_COL:
1166 42 : if (lg(Y)-1!=n) pari_err_DIM("gaussmodulo");
1167 42 : break;
1168 0 : default: pari_err_TYPE("gaussmodulo",Y);
1169 : return NULL; /* LCOV_EXCL_LINE */
1170 : }
1171 42 : H = ZM_hnfall_i(shallowconcat(M,delta), &U, 1);
1172 42 : Y = hnf_solve(H,Y); if (!Y) return gen_0;
1173 35 : l = lg(H); /* may be smaller than lM if some moduli are 0 */
1174 35 : n = l-1;
1175 35 : m = lg(U)-1 - n;
1176 35 : u1 = cgetg(m+1,t_MAT);
1177 35 : u2 = cgetg(n+1,t_MAT);
1178 84 : for (j=1; j<=m; j++) { GEN c = gel(U,j); setlg(c,lM); gel(u1,j) = c; }
1179 35 : U += m;
1180 77 : for (j=1; j<=n; j++) { GEN c = gel(U,j); setlg(c,lM); gel(u2,j) = c; }
1181 : /* (u1 u2)
1182 : * (M D) (* * ) = (0 H) */
1183 35 : u1 = ZM_lll(u1, 0.75, LLL_INPLACE);
1184 35 : Y = ZM_ZC_mul(u2,Y);
1185 35 : x = ZC_reducemodmatrix(Y, u1);
1186 35 : if (!pU) x = gerepileupto(av, x);
1187 : else
1188 : {
1189 14 : gerepileall(av, 2, &x, &u1);
1190 14 : *pU = u1;
1191 : }
1192 35 : return x;
1193 : }
1194 : /* to be used when one D_i is 0 */
1195 : static GEN
1196 70 : msolvemod0(GEN M, GEN D, GEN Y, long flag)
1197 : {
1198 70 : pari_sp av = avma;
1199 : GEN y, x, U;
1200 70 : if (!flag) return gaussmoduloall(M,D,Y,NULL);
1201 28 : y = cgetg(3,t_VEC);
1202 28 : x = gaussmoduloall(M,D,Y,&U);
1203 28 : if (x == gen_0) { set_avma(av); return gen_0; }
1204 21 : gel(y,1) = x;
1205 21 : gel(y,2) = U; return y;
1206 :
1207 : }
1208 : GEN
1209 44135 : matsolvemod(GEN M, GEN D, GEN Y, long flag)
1210 : {
1211 44135 : pari_sp av = avma;
1212 44135 : long m, n, i, char0 = 0;
1213 44135 : if (typ(M)!=t_MAT || !RgM_is_ZM(M)) pari_err_TYPE("matsolvemod (M)",M);
1214 44121 : RgM_dimensions(M,&m,&n);
1215 44121 : if (typ(D)!=t_INT && (typ(D)!=t_COL || !RgV_is_ZV(D)))
1216 28 : pari_err_TYPE("matsolvemod (D)",D);
1217 44093 : if (n)
1218 43960 : { if (typ(D)==t_COL && lg(D)!=m+1) pari_err_DIM("matsolvemod [1]"); }
1219 : else
1220 133 : { if (typ(D)==t_COL) m = lg(D)-1; }
1221 44086 : if (typ(Y)==t_INT)
1222 43890 : Y = const_col(m,Y);
1223 196 : else if (typ(Y)!=t_COL || !RgV_is_ZV(Y)) pari_err_TYPE("matsolvemod (Y)",Y);
1224 44058 : if (!n && !m) m = lg(Y)-1;
1225 44002 : else if (m != lg(Y)-1) pari_err_DIM("matsolvemod [2]");
1226 44037 : if (typ(D)==t_INT)
1227 : {
1228 43918 : if (signe(D)<0) pari_err_DOMAIN("matsolvemod","D","<",gen_0,D);
1229 43911 : if (!signe(D)) char0 = 1;
1230 : }
1231 : else /*typ(D)==t_COL*/
1232 301 : for (i=1; i<=m; i++)
1233 : {
1234 189 : if (signe(gel(D,i))<0)
1235 7 : pari_err_DOMAIN("matsolvemod","D[i]","<",gen_0,gel(D,i));
1236 182 : if (!signe(gel(D,i))) char0 = 1;
1237 : }
1238 87976 : return gerepilecopy(av, char0? msolvemod0(M,D,Y,flag)
1239 43953 : : matsolvemod_finite(M,D,Y,flag));
1240 : }
1241 : GEN
1242 0 : gaussmodulo2(GEN M, GEN D, GEN Y) { return matsolvemod(M,D,Y, 1); }
1243 : GEN
1244 0 : gaussmodulo(GEN M, GEN D, GEN Y) { return matsolvemod(M,D,Y, 0); }
|
