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 : /* bb_hermite R:
28 : * - add(a,b): a+b
29 : * - neg(a): -a
30 : * - mul(a,b): a*b
31 : * - extgcd(a,b,&small): [d,U] with d in R and U in GL_2(R)
32 : * such that [0;d] = [a;b]*U. Set small==1 to assert that U is a 'small'
33 : * operation (no red needed).
34 : * - rann(a): b in R such that b*R = {x in R | a*x==0}
35 : * - lquo(a,b,&r): q in R such that r=a-b*q is a canonical representative
36 : * of the image of a in R/b*R. The canonical lift of 0 must be 0.
37 : * - unit(a): u unit in R^* such that a*u is a canonical generator of
38 : * 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 : struct bb_hermite
52 : {
53 : GEN (*add)(void*, GEN, GEN);
54 : GEN (*neg)(void*, GEN);
55 : GEN (*mul)(void*, GEN, GEN);
56 : GEN (*extgcd)(void*, GEN, GEN, int*);
57 : GEN (*rann)(void*, GEN);
58 : GEN (*lquo)(void*, GEN, GEN, GEN*);
59 : GEN (*unit)(void*, GEN);
60 : int (*equal0)(GEN);
61 : int (*equal1)(GEN);
62 : GEN (*s)(void*, long);
63 : GEN (*red)(void*, GEN);
64 : };
65 :
66 : static GEN
67 40659490 : _Fp_add(void *data, GEN x, GEN y) { (void) data; return addii(x,y); }
68 :
69 : static GEN
70 5040142 : _Fp_neg(void *data, GEN x) { (void) data; return negi(x); }
71 :
72 : static GEN
73 48438943 : _Fp_mul(void *data, GEN x, GEN y) { (void) data; return mulii(x,y); }
74 :
75 : static GEN
76 2405927 : _Fp_rann(void *data, GEN x)
77 : {
78 2405927 : GEN d, N = (GEN)data;
79 2405927 : if (!signe(x)) return gen_1;
80 2238552 : d = gcdii(x,N);
81 2238552 : return modii(diviiexact(N,d),N);
82 : }
83 :
84 : static GEN
85 2393394 : _Fp_lquo(void *data, GEN x, GEN y, GEN* r) { (void) data; return truedvmdii(x,y,r); }
86 :
87 : /* D=MN, p|M => !p|a, p|N => p|a, return M */
88 : static GEN
89 28 : Z_split(GEN D, GEN a)
90 : {
91 : long i, n;
92 : GEN N;
93 28 : n = expi(D);
94 28 : n = n<2 ? 1 : expu(n)+1;
95 196 : for (i=1;i<=n;i++)
96 168 : a = Fp_sqr(a,D);
97 28 : N = gcdii(a,D);
98 28 : return diviiexact(D,N);
99 : }
100 :
101 : /* c s.t. gcd(a+cb,N) = gcd(a,b,N) without factoring */
102 : static GEN
103 28 : Z_stab(GEN a, GEN b, GEN N)
104 : {
105 : GEN g, a2, N2;
106 28 : g = gcdii(a,b);
107 28 : g = gcdii(g,N);
108 28 : N2 = diviiexact(N,g);
109 28 : a2 = diviiexact(a,g);
110 28 : return Z_split(N2,a2);
111 : }
112 :
113 : static GEN
114 7576617 : _Fp_unit(void *data, GEN x)
115 : {
116 7576617 : GEN g,s,v,d,N=(GEN)data,N2;
117 : long i;
118 7576617 : if (!signe(x)) return NULL;
119 7114577 : g = bezout(x,N,&s,&v);
120 7114577 : if (equali1(g) || equali1(gcdii(s,N))) return mkvec2(g,s);
121 44729 : N2 = diviiexact(N,g);
122 61355 : for (i=0; i<5; i++)
123 : {
124 61327 : s = addii(s,N2);
125 61327 : if (equali1(gcdii(s,N))) return mkvec2(g,s);
126 : }
127 28 : d = Z_stab(s,N2,N);
128 28 : d = mulii(d,N2);
129 28 : v = Fp_add(s,d,N);
130 28 : if (equali1(v)) return NULL;
131 28 : return mkvec2(g,v);
132 : }
133 :
134 : static GEN
135 3894802 : _Fp_extgcd(void *data, GEN x, GEN y, int* smallop)
136 : {
137 : GEN d,u,v,m;
138 3894802 : if (equali1(y))
139 : {
140 1188737 : *smallop = 1;
141 1188737 : return mkvec2(y,mkmat2(
142 : mkcol2(gen_1,Fp_neg(x,(GEN)data)),
143 : mkcol2(gen_0,gen_1)));
144 : }
145 2706065 : *smallop = 0;
146 2706065 : d = bezout(x,y,&u,&v);
147 2706065 : if (!signe(d)) return mkvec2(d,matid(2));
148 2706065 : m = cgetg(3,t_MAT);
149 2706065 : m = mkmat2(
150 : mkcol2(diviiexact(y,d),negi(diviiexact(x,d))),
151 : mkcol2(u,v));
152 2706065 : return mkvec2(d,m);
153 : }
154 :
155 : static int
156 108211941 : _Fp_equal0(GEN x) { return !signe(x); }
157 :
158 : static int
159 30831288 : _Fp_equal1(GEN x) { return equali1(x); }
160 :
161 : static GEN
162 16039028 : _Fp_s(void *data, long x)
163 : {
164 16039028 : if (!x) return gen_0;
165 2154704 : if (x==1) return gen_1;
166 0 : return modsi(x,(GEN)data);
167 : }
168 :
169 : static GEN
170 49965381 : _Fp_red(void *data, GEN x) { return Fp_red(x, (GEN)data); }
171 :
172 : /* p not necessarily prime */
173 : static const struct bb_hermite Fp_hermite=
174 : {_Fp_add,_Fp_neg,_Fp_mul,_Fp_extgcd,_Fp_rann,_Fp_lquo,_Fp_unit,_Fp_equal0,_Fp_equal1,_Fp_s,_Fp_red};
175 :
176 : static const struct bb_hermite*
177 970602 : get_Fp_hermite(void **data, GEN p)
178 : {
179 970602 : *data = (void*)p; return &Fp_hermite;
180 : }
181 :
182 : static void
183 18793932 : gen_redcol(GEN C, long lim, void* data, const struct bb_hermite *R)
184 : {
185 : long i;
186 69566969 : for (i=1; i<=lim; i++)
187 50773037 : if (!R->equal0(gel(C,i)))
188 25815048 : gel(C,i) = R->red(data, gel(C,i));
189 18793932 : }
190 :
191 : static GEN
192 : /* return NULL if a==0 */
193 : /* assume C*a is zero after lim */
194 20964572 : gen_rightmulcol(GEN C, GEN a, long lim, int fillzeros, void* data, const struct bb_hermite *R)
195 : {
196 : GEN Ca,zero;
197 : long i;
198 20964572 : if (R->equal1(a)) return C;
199 11934653 : if (R->equal0(a)) return NULL;
200 8250460 : Ca = cgetg(lg(C),t_COL);
201 35128565 : for (i=1; i<=lim; i++)
202 26878105 : gel(Ca,i) = R->mul(data, gel(C,i), a);
203 8250460 : if (fillzeros && lim+1 < lg(C))
204 : {
205 6352775 : zero = R->s(data,0);
206 30021597 : for (i=lim+1; i<lg(C); i++)
207 23668822 : gel(Ca,i) = zero;
208 : }
209 8250460 : return Ca;
210 : }
211 :
212 : static void
213 : /* C1 <- C1 + C2 */
214 : /* assume C2[i]==0 for i>lim */
215 5916086 : gen_addcol(GEN C1, GEN C2, long lim, void* data, const struct bb_hermite *R)
216 : {
217 : long i;
218 30710881 : for (i=1; i<=lim; i++)
219 24794795 : gel(C1,i) = R->add(data, gel(C1,i), gel(C2,i));
220 5916086 : }
221 :
222 : static void
223 : /* H[,i] <- H[,i] + C*a */
224 : /* assume C is zero after lim */
225 4468053 : gen_addrightmul(GEN H, GEN C, GEN a, long i, long lim, void* data, const struct bb_hermite *R)
226 : {
227 : GEN Ca;
228 4468053 : if (R->equal0(a)) return;
229 1521855 : Ca = gen_rightmulcol(C, a, lim, 0, data, R);
230 1521855 : gen_addcol(gel(H,i), Ca, lim, data, R);
231 : }
232 :
233 : static GEN
234 3374600 : gen_zerocol(long n, void* data, const struct bb_hermite *R)
235 : {
236 3374600 : GEN C = cgetg(n+1,t_COL), zero = R->s(data, 0);
237 : long i;
238 16200152 : for (i=1; i<=n; i++) gel(C,i) = zero;
239 3374600 : return C;
240 : }
241 :
242 : static GEN
243 1803815 : gen_zeromat(long m, long n, void* data, const struct bb_hermite *R)
244 : {
245 1803815 : GEN M = cgetg(n+1,t_MAT);
246 : long i;
247 5034530 : for (i=1; i<=n; i++) gel(M,i) = gen_zerocol(m, data, R);
248 1803815 : return M;
249 : }
250 :
251 : static GEN
252 371469 : gen_colei(long n, long i, void* data, const struct bb_hermite *R)
253 : {
254 371469 : GEN C = cgetg(n+1,t_COL), zero = R->s(data, 0);
255 : long j;
256 2882033 : for (j=1; j<=n; j++)
257 2510564 : if (i!=j) gel(C,j) = zero;
258 371469 : else gel(C,j) = R->s(data,1);
259 371469 : return C;
260 : }
261 :
262 : static GEN
263 56091 : gen_matid_hermite(long n, void* data, const struct bb_hermite *R)
264 : {
265 56091 : GEN M = cgetg(n+1,t_MAT);
266 : long i;
267 196231 : for (i=1; i<=n; i++) gel(M,i) = gen_colei(n, i, data, R);
268 56091 : return M;
269 : }
270 :
271 : static GEN
272 2714946 : gen_matmul_hermite(GEN A, GEN B, void* data, const struct bb_hermite *R)
273 : {
274 2714946 : GEN M,sum,prod,zero = R->s(data,0);
275 : long a,b,c,c2,i,j,k;
276 2714946 : RgM_dimensions(A,&a,&c);
277 2714946 : RgM_dimensions(B,&c2,&b);
278 2714946 : if (c!=c2) pari_err_DIM("gen_matmul_hermite");
279 2714946 : M = cgetg(b+1,t_MAT);
280 5850424 : for (j=1; j<=b; j++)
281 : {
282 3135478 : gel(M,j) = cgetg(a+1,t_COL);
283 8769903 : for (i=1; i<=a; i++)
284 : {
285 5634425 : sum = zero;
286 20332204 : for (k=1; k<=c; k++)
287 : {
288 14697779 : prod = R->mul(data, gcoeff(A,i,k), gcoeff(B,k,j));
289 14697779 : sum = R->add(data, sum, prod);
290 : }
291 5634425 : gcoeff(M,i,j) = sum;
292 : }
293 3135478 : gen_redcol(gel(M,j), a, data, R);
294 : }
295 2714946 : return M;
296 : }
297 :
298 : static void
299 : /* U = [u1,u2]~, C <- A*u1 + B*u2 */
300 : /* assume both A, B and C are zero after lim */
301 7967278 : gen_rightlincomb(GEN A, GEN B, GEN U, GEN *C, long lim, void* data, const struct bb_hermite *R)
302 : {
303 : GEN Au1, Bu2;
304 7967278 : Au1 = gen_rightmulcol(A, gel(U,1), lim, 1, data, R);
305 7967278 : Bu2 = gen_rightmulcol(B, gel(U,2), lim, 1, data, R);
306 7967278 : if (!Au1 && !Bu2) { *C = gen_zerocol(lg(A)-1, data, R); return; }
307 7967278 : if (!Au1) { *C = Bu2; return; }
308 4790890 : if (!Bu2) { *C = Au1; return; }
309 4394105 : gen_addcol(Au1, Bu2, lim, data, R);
310 4394105 : *C = Au1;
311 : }
312 :
313 : static void
314 : /* (H[,i] | H[,j]) <- (H[,i] | H[,j]) * U */
315 : /* assume both columns are zero after lim */
316 3983639 : gen_elem(GEN H, GEN U, long i, long j, long lim, void* data, const struct bb_hermite *R)
317 : {
318 : GEN Hi, Hj;
319 3983639 : Hi = shallowcopy(gel(H,i));
320 3983639 : Hj = shallowcopy(gel(H,j));
321 3983639 : gen_rightlincomb(Hi, Hj, gel(U,1), &gel(H,i), lim, data, R);
322 3983639 : gen_rightlincomb(Hi, Hj, gel(U,2), &gel(H,j), lim, data, R);
323 3983639 : }
324 :
325 : static int
326 : /* assume C is zero after lim */
327 4727336 : gen_is_zerocol(GEN C, long lim, void* data, const struct bb_hermite *R)
328 : {
329 : long i;
330 : (void) data;
331 8657922 : for (i=1; i<=lim; i++)
332 6755405 : if (!R->equal0(gel(C,i))) return 0;
333 1902517 : return 1;
334 : }
335 :
336 : /* The mkop* functions return NULL if the corresponding operation is the identity */
337 :
338 : static GEN
339 : /* Ci <- Ci + Cj*a */
340 3217365 : mkoptransv(long i, long j, GEN a, void* data, const struct bb_hermite *R)
341 : {
342 3217365 : a = R->red(data,a);
343 3217365 : if (R->equal0(a)) return NULL;
344 1170261 : return mkvec2(mkvecsmall2(i,j),a);
345 : }
346 :
347 : static GEN
348 : /* (Ci|Cj) <- (Ci|Cj)*U */
349 977682 : mkopU(long i, long j, GEN U, void* data, const struct bb_hermite *R)
350 : {
351 977682 : if (R->equal1(gcoeff(U,1,1)) && R->equal0(gcoeff(U,1,2))
352 435163 : && R->equal1(gcoeff(U,2,2))) return mkoptransv(i,j,gcoeff(U,2,1),data,R);
353 542519 : return mkvec2(mkvecsmall3(i,j,0),U);
354 : }
355 :
356 : static GEN
357 : /* Ci <- Ci*u */
358 2183225 : mkopmul(long i, GEN u, const struct bb_hermite *R)
359 : {
360 2183225 : if (R->equal1(u)) return NULL;
361 194178 : return mkvec2(mkvecsmall(i),u);
362 : }
363 :
364 : static GEN
365 : /* Ci <-> Cj */
366 45843 : mkopswap(long i, long j)
367 : {
368 45843 : return mkvec(mkvecsmall2(i,j));
369 : }
370 :
371 : /* M: t_MAT. Apply the operation op to M by right multiplication. */
372 : static void
373 374465 : gen_rightapply(GEN M, GEN op, void* data, const struct bb_hermite *R)
374 : {
375 : GEN M2, ind, X;
376 374465 : long i, j, m = lg(gel(M,1))-1;
377 374465 : switch (typ(op))
378 : {
379 56014 : case t_VECSMALL:
380 56014 : M2 = vecpermute(M,op);
381 266056 : for (i=1; i<lg(M); i++) gel(M,i) = gel(M2,i);
382 56014 : return;
383 318451 : case t_VEC:
384 318451 : ind = gel(op,1);
385 318451 : switch (lg(op))
386 : {
387 16135 : case 2:
388 16135 : swap(gel(M,ind[1]),gel(M,ind[2]));
389 16135 : return;
390 302316 : case 3:
391 302316 : X = gel(op,2);
392 302316 : i = ind[1];
393 302316 : switch (lg(ind))
394 : {
395 37786 : case 2:
396 37786 : gel(M,i) = gen_rightmulcol(gel(M,i), X, m, 0, data, R);
397 37786 : gen_redcol(gel(M,i), m, data, R);
398 37786 : return;
399 175693 : case 3:
400 175693 : gen_addrightmul(M, gel(M,ind[2]), X, i, m, data, R);
401 175693 : gen_redcol(gel(M,i), m, data, R);
402 175693 : return;
403 88837 : case 4:
404 88837 : j = ind[2];
405 88837 : gen_elem(M, X, i, j, m, data, R);
406 88837 : gen_redcol(gel(M,i), m, data, R);
407 88837 : gen_redcol(gel(M,j), m, data, R);
408 88837 : return;
409 : }
410 : }
411 : }
412 : }
413 :
414 : /* C: t_COL. Apply the operation op to C by left multiplication. */
415 : static void
416 15435898 : gen_leftapply(GEN C, GEN op, void* data, const struct bb_hermite *R)
417 : {
418 : GEN C2, ind, X;
419 : long i, j;
420 15435898 : switch (typ(op))
421 : {
422 591164 : case t_VECSMALL:
423 591164 : C2 = vecpermute(C,perm_inv(op));
424 5504947 : for (i=1; i<lg(C); i++) gel(C,i) = gel(C2,i);
425 591164 : return;
426 14844734 : case t_VEC:
427 14844734 : ind = gel(op,1);
428 14844734 : switch (lg(op))
429 : {
430 169134 : case 2:
431 169134 : swap(gel(C,ind[1]),gel(C,ind[2]));
432 169134 : return;
433 14675600 : case 3:
434 14675600 : X = gel(op,2);
435 14675600 : i = ind[1];
436 14675600 : switch (lg(ind))
437 : {
438 714127 : case 2:
439 714127 : gel(C,i) = R->mul(data, X, gel(C,i));
440 714127 : gel(C,i) = R->red(data, gel(C,i));
441 714127 : return;
442 11462526 : case 3:
443 11462526 : j = ind[2];
444 11462526 : if (R->equal0(gel(C,i))) return;
445 1166811 : gel(C,j) = R->add(data, gel(C,j), R->mul(data, X, gel(C,i)));
446 1166811 : return;
447 2498947 : case 4:
448 2498947 : j = ind[2];
449 2498947 : C2 = gen_matmul_hermite(X, mkmat(mkcol2(gel(C,i),gel(C,j))), data, R);
450 2498947 : gel(C,i) = gcoeff(C2,1,1);
451 2498947 : gel(C,j) = gcoeff(C2,2,1);
452 2498947 : return;
453 : }
454 : }
455 : }
456 : }
457 :
458 : /* \prod_i det ops[i]. Only makes sense if R is commutative. */
459 : static GEN
460 42 : gen_detops(GEN ops, void* data, const struct bb_hermite *R)
461 : {
462 42 : GEN d = R->s(data,1);
463 42 : long i, l = lg(ops);
464 231 : for (i = 1; i < l; i++)
465 : {
466 189 : GEN X, op = gel(ops,i);
467 189 : switch (typ(op))
468 : {
469 0 : case t_VECSMALL:
470 0 : if (perm_sign(op) < 0) d = R->neg(data,d);
471 0 : break;
472 189 : case t_VEC:
473 189 : switch (lg(op))
474 : {
475 0 : case 2:
476 0 : d = R->neg(data,d);
477 0 : break;
478 189 : case 3:
479 189 : X = gel(op,2);
480 189 : switch (lg(gel(op,1)))
481 : {
482 0 : case 2:
483 0 : d = R->mul(data, d, X);
484 0 : d = R->red(data, d);
485 0 : break;
486 105 : case 4:
487 : {
488 105 : GEN A = gcoeff(X,1,1), B = gcoeff(X,1,2);
489 105 : GEN C = gcoeff(X,2,1), D = gcoeff(X,2,2);
490 105 : GEN AD = R->mul(data,A,D);
491 105 : GEN BC = R->mul(data,B,C);
492 105 : d = R->mul(data, d, R->add(data, AD, R->neg(data,BC)));
493 105 : d = R->red(data, d);
494 105 : break;
495 : }
496 : }
497 189 : break;
498 : }
499 189 : break;
500 : }
501 : }
502 42 : return d;
503 : }
504 :
505 : static int
506 220794 : gen_is_inv(GEN x, void* data, const struct bb_hermite *R)
507 : {
508 220794 : GEN u = R->unit(data, x);
509 220794 : if (!u) return R->equal1(x);
510 75782 : return R->equal1(gel(u,1));
511 : }
512 :
513 : static long
514 201712 : gen_last_inv_diago(GEN A, void* data, const struct bb_hermite *R)
515 : {
516 : long i,m,n,j;
517 201712 : RgM_dimensions(A,&m,&n);
518 235774 : for (i=1,j=n-m+1; i<=m; i++,j++)
519 220794 : if (!gen_is_inv(gcoeff(A,i,j),data,R)) return i-1;
520 14980 : return m;
521 : }
522 :
523 : static GEN
524 : /* remove_zerocols: 0 none, 1 until square, 2 all */
525 : /* early abort: if not right-invertible, abort, return NULL, and set ops to the
526 : * noninvertible pivot */
527 1070534 : 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)
528 : {
529 1070534 : pari_sp av = avma, av1;
530 1070534 : 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);
531 1070534 : long m,n,i,j,s,si,i2,si2,nbz,lim,extra,maxop=0,nbop=0,lastinv=0;
532 : int smallop;
533 :
534 1070534 : av1 = avma;
535 :
536 1070534 : RgM_dimensions(A,&m,&n);
537 1070534 : if (early_abort && n<m)
538 : {
539 14000 : if (ops) *ops = zero;
540 14000 : return NULL;
541 : }
542 1056534 : if (n<m+1)
543 : {
544 954124 : extra = m+1-n;
545 954124 : H = shallowmatconcat(mkvec2(gen_zeromat(m,extra,data,R),A));
546 : }
547 : else
548 : {
549 102410 : extra = 0;
550 102410 : H = RgM_shallowcopy(A);
551 : }
552 1056534 : RgM_dimensions(H,&m,&n);
553 1056534 : s = n-m; /* shift */
554 :
555 1056534 : if(ops)
556 : {
557 854822 : maxop = m*n + (m*(m+1))/2 + 1;
558 854822 : *ops = zerovec(maxop); /* filled with placeholders for GC */
559 : }
560 :
561 : /* put in triangular form */
562 4060117 : for (i=m,si=s+m; i>0 && si>extra; i--,si--) /* si = s+i */
563 : {
564 3008497 : if (R->red) gcoeff(H,i,si) = R->red(data, gcoeff(H,i,si));
565 : /* bottom-right diagonal */
566 13466430 : for (j = extra+1; j < si; j++)
567 : {
568 10457933 : if (R->red) gcoeff(H,i,j) = R->red(data, gcoeff(H,i,j));
569 10457933 : if (R->equal0(gcoeff(H,i,j))) continue;
570 3738421 : U = R->extgcd(data, gcoeff(H,i,j), gcoeff(H,i,si), &smallop);
571 3738421 : d = gel(U,1);
572 3738421 : U = gel(U,2);
573 3738421 : if (n>10)
574 : {
575 : /* normalize diagonal coefficient -> faster reductions on this row */
576 2490903 : u = R->unit(data, d);
577 2490903 : if (u)
578 : {
579 2490903 : g = gel(u,1);
580 2490903 : u = gel(u,2);
581 2490903 : gcoeff(U,1,2) = R->mul(data, gcoeff(U,1,2), u);
582 2490903 : gcoeff(U,2,2) = R->mul(data, gcoeff(U,2,2), u);
583 2490903 : d = g;
584 : }
585 : }
586 3738421 : gen_elem(H, U, j, si, i-1, data, R);
587 3738421 : if (ops)
588 : {
589 941639 : op = mkopU(j,si,U,data,R);
590 941639 : if (op) { nbop++; gel(*ops, nbop) = op; }
591 : }
592 3738421 : gcoeff(H,i,j) = zero;
593 3738421 : gcoeff(H,i,si) = d;
594 3738421 : if (R->red && !smallop)
595 : {
596 2563614 : gen_redcol(gel(H,si), i-1, data, R);
597 2563614 : gen_redcol(gel(H,j), i-1, data, R);
598 : }
599 : }
600 :
601 3008497 : if (early_abort)
602 : {
603 1757492 : d = gcoeff(H,i,si);
604 1757492 : u = R->unit(data, d);
605 1757492 : if (u) d = gel(u,1);
606 1757492 : if (!R->equal1(d))
607 : {
608 4914 : if (ops) *ops = d;
609 4914 : return NULL;
610 : }
611 : }
612 :
613 3003583 : if (gc_needed(av,1))
614 : {
615 20 : if (DEBUGMEM>1) pari_warn(warnmem,"gen_howell[1]. i=%ld",i);
616 20 : (void)gc_all(av1,ops?2:1,&H,ops);
617 : }
618 : }
619 :
620 1051620 : if (!ops)
621 201712 : lastinv = gen_last_inv_diago(H, data, R);
622 :
623 : /* put in reduced Howell form */
624 1051620 : if (!only_triangular)
625 : {
626 4193068 : for (i=m,si=s+m; i>0; i--,si--) /* si = s+i */
627 : {
628 : /* normalize diagonal coefficient */
629 3141490 : if (i<=lastinv) /* lastinv>0 => !ops */
630 34062 : gcoeff(H,i,si) = one;
631 : else
632 : {
633 3107428 : u = R->unit(data,gcoeff(H,i,si));
634 3107428 : if (u)
635 : {
636 2790547 : g = gel(u,1);
637 2790547 : u = gel(u,2);
638 2790547 : gel(H,si) = gen_rightmulcol(gel(H,si), u, i-1, 1, data, R);
639 2790547 : gcoeff(H,i,si) = g;
640 2790547 : if (R->red) gen_redcol(gel(H,si), i-1, data, R);
641 2790547 : if (ops)
642 : {
643 2183225 : op = mkopmul(si,u,R);
644 2183225 : if (op) { nbop++; gel(*ops,nbop) = op; }
645 : }
646 : }
647 316881 : else if (R->red) gcoeff(H,i,si) = R->red(data, gcoeff(H,i,si));
648 : }
649 3141490 : piv = gcoeff(H,i,si);
650 :
651 : /* reduce above diagonal */
652 3141490 : if (!R->equal0(piv))
653 : {
654 2824609 : C = gel(H,si);
655 7167572 : for (j=si+1; j<=n; j++)
656 : {
657 4342963 : if (i<=lastinv) /* lastinv>0 => !ops */
658 50603 : gcoeff(H,i,j) = zero;
659 : else
660 : {
661 4292360 : gcoeff(H,i,j) = R->red(data, gcoeff(H,i,j));
662 4292360 : if (R->equal1(piv)) { q = gcoeff(H,i,j); r = zero; }
663 1645717 : else q = R->lquo(data, gcoeff(H,i,j), piv, &r);
664 4292360 : q = R->neg(data,q);
665 4292360 : gen_addrightmul(H, C, q, j, i-1, data, R);
666 4292360 : if (ops)
667 : {
668 2633711 : op = mkoptransv(j,si,q,data,R);
669 2633711 : if (op) { nbop++; gel(*ops,nbop) = op; }
670 : }
671 4292360 : gcoeff(H,i,j) = r;
672 : }
673 : }
674 : }
675 :
676 : /* ensure Howell property */
677 3141490 : if (i>1)
678 : {
679 2089996 : a = R->rann(data, piv);
680 2089996 : if (!R->equal0(a))
681 : {
682 568682 : gel(H,1) = gen_rightmulcol(gel(H,si), a, i-1, 1, data, R);
683 568682 : if (gel(H,1) == gel(H,si)) gel(H,1) = shallowcopy(gel(H,1)); /* in case rightmulcol cheated */
684 568682 : if (ops)
685 : {
686 148491 : op = mkoptransv(1,si,a,data,R);
687 148491 : if (op) { nbop++; gel(*ops,nbop) = op; }
688 : }
689 2423036 : for (i2=i-1,si2=s+i2; i2>0; i2--,si2--)
690 : {
691 1854354 : if (R->red) gcoeff(H,i2,1) = R->red(data, gcoeff(H,i2,1));
692 1854354 : if (R->equal0(gcoeff(H,i2,1))) continue;
693 288711 : if (R->red) gcoeff(H,i2,si2) = R->red(data, gcoeff(H,i2,si2));
694 288711 : if (R->equal0(gcoeff(H,i2,si2)))
695 : {
696 132330 : swap(gel(H,1), gel(H,si2));
697 132330 : if (ops) { nbop++; gel(*ops,nbop) = mkopswap(1,si2); }
698 132330 : continue;
699 : }
700 156381 : U = R->extgcd(data, gcoeff(H,i2,1), gcoeff(H,i2,si2), &smallop);
701 156381 : d = gel(U,1);
702 156381 : U = gel(U,2);
703 156381 : gen_elem(H, U, 1, si2, i2-1, data, R);
704 156381 : if (ops)
705 : {
706 36043 : op = mkopU(1,si2,U,data,R);
707 36043 : if (op) { nbop++; gel(*ops,nbop) = op; }
708 : }
709 156381 : gcoeff(H,i2,1) = zero;
710 156381 : gcoeff(H,i2,si2) = d;
711 156381 : if (R->red && !smallop)
712 : {
713 142451 : gen_redcol(gel(H,si2), i2, data, R);
714 142451 : gen_redcol(gel(H,1), i2-1, data, R);
715 : }
716 : }
717 : }
718 : }
719 :
720 3141490 : if (gc_needed(av,1))
721 : {
722 12 : if (DEBUGMEM>1) pari_warn(warnmem,"gen_howell[2]. i=%ld",i);
723 12 : (void)gc_all(av1,ops?3:2,&H,&piv,ops);
724 : }
725 : }
726 : }
727 :
728 1051620 : if (R->red)
729 5735003 : for (j=1; j<=n; j++)
730 : {
731 4683383 : lim = maxss(0,m-n+j);
732 4683383 : gen_redcol(gel(H,j), lim, data, R);
733 18341466 : for (i=lim+1; i<=m; i++) gcoeff(H,i,j) = zero;
734 : }
735 :
736 : /* put zero columns first */
737 1051620 : iszero = cgetg(n+1,t_VECSMALL);
738 :
739 1051620 : nbz = 0;
740 5735003 : for (i=1; i<=n; i++)
741 : {
742 4683383 : iszero[i] = gen_is_zerocol(gel(H,i), maxss(0,m-n+i), data, R);
743 4683383 : if (iszero[i]) nbz++;
744 : }
745 :
746 1051620 : j = 1;
747 1051620 : if (permute_zerocols)
748 : {
749 156121 : perm = cgetg(n+1, t_VECSMALL);
750 913759 : for (i=1; i<=n; i++)
751 757638 : if (iszero[i])
752 : {
753 320537 : perm[j] = i;
754 320537 : j++;
755 : }
756 : }
757 895499 : else perm = cgetg(n-nbz+1, t_VECSMALL);
758 5735003 : for (i=1; i<=n; i++)
759 4683383 : if (!iszero[i])
760 : {
761 2824770 : perm[j] = i;
762 2824770 : j++;
763 : }
764 :
765 1051620 : if (permute_zerocols || remove_zerocols==2) H = vecpermute(H, perm);
766 1051620 : if (permute_zerocols && remove_zerocols==2) H = vecslice(H, nbz+1, n);
767 1051620 : if (remove_zerocols==1) H = vecslice(H, s+1, n);
768 1051620 : if (permute_zerocols && ops) { nbop++; gel(*ops,nbop) = perm; }
769 :
770 1051620 : if (ops) { setlg(*ops, nbop+1); } /* should have nbop <= maxop */
771 :
772 1051620 : return H;
773 : }
774 :
775 : static GEN
776 101780 : 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)
777 : {
778 101780 : pari_sp av = avma;
779 101780 : GEN H = gen_howell_i(A, remove_zerocols, permute_zerocols, early_abort, only_triangular, ops, data, R);
780 101780 : return gc_all(av, ops?2:1, &H, ops);
781 : }
782 :
783 : static GEN
784 157794 : gen_matimage(GEN A, GEN* U, void *data, const struct bb_hermite *R)
785 : {
786 : GEN ops, H;
787 157794 : if (U)
788 : {
789 56014 : pari_sp av = avma, av1;
790 : long m, n, i, r, n2, pergc;
791 56014 : RgM_dimensions(A,&m,&n);
792 56014 : H = gen_howell_i(A, 2, 1, 0, 0, &ops, data, R);
793 56014 : av1 = avma;
794 56014 : r = lg(H)-1;
795 56014 : *U = shallowmatconcat(mkvec2(gen_zeromat(n, maxss(0,m-n+1), data, R), gen_matid_hermite(n, data, R)));
796 56014 : n2 = lg(*U)-1;
797 56014 : pergc = maxss(m,n);
798 430479 : for (i=1; i<lg(ops); i++)
799 : {
800 374465 : gen_rightapply(*U, gel(ops,i), data, R);
801 374465 : if (!(i%pergc) && gc_needed(av,1))
802 : {
803 0 : if (DEBUGMEM>1) pari_warn(warnmem,"gen_matimage. i=%ld",i);
804 0 : (void)gc_all(av1,1,U);
805 : }
806 : }
807 56014 : if (r<n2) *U = vecslice(*U, n2-r+1, n2);
808 56014 : return gc_all(av, 2, &H, U);
809 : }
810 101780 : return gen_howell(A, 2, 0, 0, 0, NULL, data, R);
811 : }
812 :
813 : static GEN
814 : /* H in true Howell form: no zero columns */
815 99932 : gen_kernel_howell(GEN H, void *data, const struct bb_hermite *R)
816 : {
817 : GEN K, piv, FK;
818 : long m, n, j, j2, i;
819 99932 : RgM_dimensions(H,&m,&n);
820 99932 : K = gen_zeromat(n, n, data, R);
821 415863 : for (j=n,i=m; j>0; j--)
822 : {
823 535850 : while (R->equal0(gcoeff(H,i,j))) i--;
824 315931 : piv = gcoeff(H,i,j);
825 315931 : if (R->equal0(piv)) continue;
826 315931 : gcoeff(K,j,j) = R->rann(data, piv);
827 315931 : if (j<n)
828 : {
829 215999 : FK = gen_matmul_hermite(matslice(H,i,i,j+1,n), matslice(K, j+1, n, j+1, n), data, R);
830 852530 : for (j2=j+1; j2<=n; j2++)
831 636531 : gcoeff(K,j,j2) = R->neg(data, R->lquo(data, gcoeff(FK,1,j2-j), piv, NULL));
832 : /* remainder has to be zero */
833 : }
834 : }
835 99932 : return K;
836 : }
837 :
838 : static GEN
839 : /* (H,ops) Howell form of A, n = number of columns of A, return a kernel of A */
840 100009 : gen_kernel_from_howell(GEN H, GEN ops, long n, void *data, const struct bb_hermite *R)
841 : {
842 : pari_sp av;
843 : GEN K, KH, zC;
844 : long m, r, n2, nbz, i, o, extra, j;
845 100009 : RgM_dimensions(H,&m,&r);
846 100009 : if (!r) return gen_matid_hermite(n, data, R); /* zerology: what if 0==1 in R? */
847 99932 : n2 = maxss(n,m+1);
848 99932 : extra = n2-n;
849 99932 : nbz = n2-r;
850 : /* compute kernel of augmented matrix */
851 99932 : KH = gen_kernel_howell(H, data, R);
852 99932 : zC = gen_zerocol(nbz, data, R);
853 99932 : K = cgetg(nbz+r+1, t_MAT);
854 331261 : for (i=1; i<=nbz; i++)
855 231329 : gel(K,i) = gen_colei(nbz+r, i, data, R);
856 415863 : for (i=1; i<=r; i++)
857 315931 : gel(K,nbz+i) = shallowconcat(zC, gel(KH,i));
858 647192 : for (i=1; i<lg(K); i++)
859 : {
860 547260 : av = avma;
861 13698766 : for (o=lg(ops)-1; o>0; o--)
862 13151506 : gen_leftapply(gel(K,i), gel(ops,o), data, R);
863 547260 : gen_redcol(gel(K,i), nbz+r, data, R);
864 547260 : (void)gc_all(av, 1, &gel(K,i));
865 : }
866 : /* deduce kernel of original matrix */
867 99932 : K = rowpermute(K, cyclic_perm(n2,extra));
868 99932 : K = gen_howell_i(K, 2, 0, 0, 0, NULL, data, R);
869 226723 : for (j=lg(K)-1, i=n2; j>0; j--)
870 : {
871 301126 : while (R->equal0(gcoeff(K,i,j))) i--;
872 198079 : if (i<=n) return matslice(K, 1, n, 1, j);
873 : }
874 28644 : return cgetg(1,t_MAT);
875 : }
876 :
877 : /* not GC-clean */
878 : static GEN
879 56154 : gen_kernel(GEN A, GEN* im, void *data, const struct bb_hermite *R)
880 : {
881 56154 : pari_sp av = avma;
882 56154 : long n = lg(A)-1;
883 : GEN H, ops, K;
884 56154 : H = gen_howell_i(A, 2, 1, 0, 0, &ops, data, R);
885 56154 : (void)gc_all(av,2,&H,&ops);
886 56154 : K = gen_kernel_from_howell(H, ops, n, data, R);
887 56154 : if (im) *im = H;
888 56154 : return K;
889 : }
890 :
891 : /* right inverse, not GC-clean */
892 : static GEN
893 712659 : gen_inv(GEN A, void* data, const struct bb_hermite *R)
894 : {
895 : pari_sp av;
896 712659 : GEN ops, H, U, un=R->s(data,1);
897 : long m,n,j,o,n2;
898 712659 : RgM_dimensions(A,&m,&n);
899 712659 : av = avma;
900 712659 : H = gen_howell_i(A, 0, 0, 1, 0, &ops, data, R);
901 712659 : if (!H) pari_err_INV("gen_inv", ops);
902 693745 : n2 = lg(H)-1;
903 693745 : ops = gc_GEN(av,ops); /* get rid of H */
904 693745 : U = gen_zeromat(n2, m, data, R);
905 2439869 : for (j=1; j<=m; j++)
906 1746124 : gcoeff(U,j+n2-m,j) = un;
907 2439869 : for (j=1; j<=m; j++)
908 : {
909 1746124 : av = avma;
910 3632307 : for (o=lg(ops)-1; o>0; o--)
911 1886183 : gen_leftapply(gel(U,j), gel(ops,o), data, R);
912 1746124 : gen_redcol(gel(U,j), n2, data, R);
913 1746124 : (void)gc_all(av, 1, &gel(U,j));
914 : }
915 693745 : if (n2>n) U = rowslice(U, n2-n+1, n2);
916 693745 : return U;
917 : }
918 :
919 : /* H true Howell form (no zero columns).
920 : * Compute Z = Y - HX canonical representative of R^m mod H(R^n) */
921 : static GEN
922 43953 : gen_reduce_mod_howell(GEN H, GEN Y, GEN *X, void* data, const struct bb_hermite *R)
923 : {
924 43953 : pari_sp av = avma;
925 : long m,n,i,j;
926 : GEN Z, q, r, C;
927 43953 : RgM_dimensions(H,&m,&n);
928 43953 : if (X) *X = gen_zerocol(n,data,R);
929 43953 : Z = shallowcopy(Y);
930 43953 : i = m;
931 155099 : for (j=n; j>0; j--)
932 : {
933 180348 : while (R->equal0(gcoeff(H,i,j))) i--;
934 111146 : q = R->lquo(data, gel(Z,i), gcoeff(H,i,j), &r);
935 111146 : gel(Z,i) = r;
936 111146 : C = gen_rightmulcol(gel(H,j), R->neg(data,q), i, 0, data, R);
937 111146 : if (C) gen_addcol(Z, C, i-1, data, R);
938 111146 : if (X) gel(*X,j) = q;
939 : }
940 43953 : gen_redcol(Z, lg(Z)-1, data, R);
941 43953 : return gc_all(av, X?2:1, &Z, X);
942 : }
943 :
944 : /* H: Howell form of A */
945 : /* (m,n): dimensions of original matrix A */
946 : /* return NULL if no solution */
947 : static GEN
948 43953 : gen_solve_from_howell(GEN H, GEN ops, long m, long n, GEN Y, void* data, const struct bb_hermite *R)
949 : {
950 43953 : pari_sp av = avma;
951 : GEN Z, X;
952 : long n2, mH, nH, i;
953 43953 : RgM_dimensions(H,&mH,&nH);
954 43953 : n2 = maxss(n,m+1);
955 43953 : Z = gen_reduce_mod_howell(H, Y, &X, data, R);
956 43953 : if (!gen_is_zerocol(Z,m,data,R)) { set_avma(av); return NULL; }
957 :
958 43904 : X = shallowconcat(zerocol(n2-nH),X);
959 442113 : for (i=lg(ops)-1; i>0; i--) gen_leftapply(X, gel(ops,i), data, R);
960 43904 : X = vecslice(X, n2-n+1, n2);
961 43904 : gen_redcol(X, n, data, R);
962 43904 : return gc_GEN(av, X);
963 : }
964 :
965 : /* return NULL if no solution, not GC-clean */
966 : static GEN
967 43953 : gen_solve(GEN A, GEN Y, GEN* K, void* data, const struct bb_hermite *R)
968 : {
969 : GEN H, ops, X;
970 : long m,n;
971 :
972 43953 : RgM_dimensions(A,&m,&n);
973 43953 : if (!n) m = lg(Y)-1;
974 43953 : H = gen_howell_i(A, 2, 1, 0, 0, &ops, data, R);
975 43953 : X = gen_solve_from_howell(H, ops, m, n, Y, data, R);
976 43953 : if (!X) return NULL;
977 43904 : if (K) *K = gen_kernel_from_howell(H, ops, n, data, R);
978 43904 : return X;
979 : }
980 :
981 : GEN
982 157829 : matimagemod(GEN A, GEN d, GEN* U)
983 : {
984 : void *data;
985 : const struct bb_hermite* R;
986 157829 : if (typ(A)!=t_MAT || !RgM_is_ZM(A)) pari_err_TYPE("matimagemod", A);
987 157815 : if (typ(d)!=t_INT) pari_err_TYPE("matimagemod", d);
988 157808 : if (signe(d)<=0) pari_err_DOMAIN("matimagemod", "d", "<=", gen_0, d);
989 157801 : if (equali1(d)) return cgetg(1,t_MAT);
990 157794 : R = get_Fp_hermite(&data, d);
991 157794 : return gen_matimage(A, U, data, R);
992 : }
993 :
994 : #if 0
995 : /* for testing purpose */
996 : GEN
997 : ZM_hnfmodid2(GEN A, GEN d)
998 : {
999 : pari_sp av = avma;
1000 : void *data;
1001 : long i;
1002 : const struct bb_hermite* R = get_Fp_hermite(&data, d);
1003 : GEN H;
1004 : if (typ(A)!=t_MAT || !RgM_is_ZM(A)) pari_err_TYPE("ZM_hnfmodid2", A);
1005 : if (typ(d)!=t_INT) pari_err_TYPE("ZM_hnfmodid2", d);
1006 : H = gen_howell_i(A, 1, 0, 0, 0, NULL, data, R);
1007 : for (i=1; i<lg(H); i++)
1008 : if (!signe(gcoeff(H,i,i))) gcoeff(H,i,i) = d;
1009 : return gc_GEN(av,H);
1010 : }
1011 : #endif
1012 :
1013 : GEN
1014 84 : matdetmod(GEN A, GEN d)
1015 : {
1016 84 : pari_sp av = avma;
1017 : void *data;
1018 : const struct bb_hermite* R;
1019 84 : long n = lg(A)-1, i;
1020 : GEN D, H, ops;
1021 84 : if (typ(A)!=t_MAT || !RgM_is_ZM(A)) pari_err_TYPE("matdetmod", A);
1022 70 : if (typ(d)!=t_INT) pari_err_TYPE("matdetmod", d);
1023 70 : if (signe(d)<=0) pari_err_DOMAIN("matdetmod", "d", "<=", gen_0, d);
1024 63 : if (!n) return equali1(d) ? gen_0 : gen_1;
1025 56 : if (n != nbrows(A)) pari_err_DIM("matdetmod");
1026 49 : if (equali1(d)) return gen_0;
1027 42 : R = get_Fp_hermite(&data, d);
1028 42 : H = gen_howell_i(A, 1, 0, 0, 1, &ops, data, R);
1029 42 : D = gen_detops(ops, data, R);
1030 42 : D = Fp_inv(D, d);
1031 203 : for (i = 1; i <= n; i++) D = Fp_mul(D, gcoeff(H,i,i), d);
1032 42 : return gc_INT(av, D);
1033 : }
1034 :
1035 : GEN
1036 56189 : matkermod(GEN A, GEN d, GEN* im)
1037 : {
1038 56189 : pari_sp av = avma;
1039 : void *data;
1040 : const struct bb_hermite* R;
1041 : long m,n;
1042 : GEN K;
1043 56189 : if (typ(A)!=t_MAT || !RgM_is_ZM(A)) pari_err_TYPE("matkermod", A);
1044 56175 : if (typ(d)!=t_INT) pari_err_TYPE("matkermod", d);
1045 56168 : if (signe(d)<=0) pari_err_DOMAIN("makermod", "d", "<=", gen_0, d);
1046 56161 : if (equali1(d)) return cgetg(1,t_MAT);
1047 56154 : R = get_Fp_hermite(&data, d);
1048 56154 : RgM_dimensions(A,&m,&n);
1049 56154 : if (!im && m>2*n) /* TODO tune treshold */
1050 4592 : A = shallowtrans(matimagemod(shallowtrans(A),d,NULL));
1051 56154 : K = gen_kernel(A, im, data, R); return gc_all(av,im?2:1,&K,im);
1052 : }
1053 :
1054 : /* left inverse */
1055 : GEN
1056 799039 : matinvmod(GEN A, GEN d)
1057 : {
1058 799039 : pari_sp av = avma;
1059 : void *data;
1060 : const struct bb_hermite* R;
1061 : GEN U;
1062 799039 : if (typ(A)!=t_MAT || !RgM_is_ZM(A)) pari_err_TYPE("matinvmod", A);
1063 799025 : if (typ(d)!=t_INT) pari_err_TYPE("matinvmod", d);
1064 799018 : if (signe(d)<=0) pari_err_DOMAIN("matinvmod", "d", "<=", gen_0, d);
1065 799011 : if (equali1(d)) {
1066 : long m,n;
1067 86352 : RgM_dimensions(A,&m,&n);
1068 86352 : if (m<n) pari_err_INV("matinvmod",A);
1069 86345 : return zeromatcopy(n,m);
1070 : }
1071 712659 : R = get_Fp_hermite(&data, d);
1072 712659 : U = gen_inv(shallowtrans(A), data, R);
1073 693745 : return gc_GEN(av, shallowtrans(U));
1074 : }
1075 :
1076 : /* assume all D[i]>0, not GC-clean */
1077 : static GEN
1078 43953 : matsolvemod_finite(GEN M, GEN D, GEN Y, long flag)
1079 : {
1080 : void *data;
1081 : const struct bb_hermite* R;
1082 : GEN X, K, d;
1083 : long m, n;
1084 :
1085 43953 : RgM_dimensions(M,&m,&n);
1086 43953 : if (typ(D)==t_COL)
1087 : { /* create extra variables for the D[i] */
1088 : GEN MD;
1089 84 : long i, i2, extra = 0;
1090 84 : d = gen_1;
1091 231 : for (i=1; i<lg(D); i++) d = lcmii(d,gel(D,i));
1092 231 : for (i=1; i<lg(D); i++) if (!equalii(gel(D,i),d)) extra++;
1093 84 : MD = cgetg(extra+1,t_MAT);
1094 84 : i2 = 1;
1095 231 : for (i=1; i<lg(D); i++)
1096 147 : if (!equalii(gel(D,i),d))
1097 : {
1098 77 : gel(MD,i2) = Rg_col_ei(gel(D,i),m,i);
1099 77 : i2++;
1100 : }
1101 84 : M = shallowconcat(M,MD);
1102 : }
1103 43869 : else d = D;
1104 :
1105 43953 : R = get_Fp_hermite(&data, d);
1106 43953 : X = gen_solve(M, Y, flag? &K: NULL, data, R);
1107 43953 : if (!X) return gen_0;
1108 43904 : X = vecslice(X,1,n);
1109 :
1110 43904 : if (flag)
1111 : {
1112 43855 : K = rowslice(K,1,n);
1113 43855 : K = hnfmodid(shallowconcat(zerocol(n),K),d);
1114 43855 : X = mkvec2(X,K);
1115 : }
1116 43904 : return X;
1117 : }
1118 :
1119 : /* Return a solution of congruence system sum M[i,j] X_j = Y_i mod D_i
1120 : * If pU != NULL, put in *pU a Z-basis of the homogeneous system.
1121 : * Works for all non negative D_i but inefficient compared to
1122 : * matsolvemod_finite; to be used only when one D_i is 0 */
1123 : static GEN
1124 77 : gaussmoduloall(GEN M, GEN D, GEN Y, GEN *pU)
1125 : {
1126 77 : pari_sp av = avma;
1127 77 : long n, m, j, l, lM = lg(M);
1128 : GEN delta, H, U, u1, u2, x;
1129 :
1130 77 : if (lM == 1)
1131 : {
1132 28 : long lY = 0;
1133 28 : switch(typ(Y))
1134 : {
1135 0 : case t_INT: break;
1136 28 : case t_COL: lY = lg(Y); break;
1137 0 : default: pari_err_TYPE("gaussmodulo",Y);
1138 : }
1139 28 : switch(typ(D))
1140 : {
1141 14 : case t_INT: break;
1142 14 : case t_COL: if (lY && lY != lg(D)) pari_err_DIM("gaussmodulo");
1143 14 : break;
1144 0 : default: pari_err_TYPE("gaussmodulo",D);
1145 : }
1146 28 : if (pU) *pU = cgetg(1, t_MAT);
1147 28 : return cgetg(1,t_COL);
1148 : }
1149 49 : n = nbrows(M);
1150 49 : switch(typ(D))
1151 : {
1152 21 : case t_COL:
1153 21 : if (lg(D)-1!=n) pari_err_DIM("gaussmodulo");
1154 21 : delta = diagonal_shallow(D); break;
1155 28 : case t_INT: delta = scalarmat_shallow(D,n); break;
1156 0 : default: pari_err_TYPE("gaussmodulo",D);
1157 : return NULL; /* LCOV_EXCL_LINE */
1158 : }
1159 49 : switch(typ(Y))
1160 : {
1161 0 : case t_INT: Y = const_col(n, Y); break;
1162 49 : case t_COL:
1163 49 : if (lg(Y)-1!=n) pari_err_DIM("gaussmodulo");
1164 49 : break;
1165 0 : default: pari_err_TYPE("gaussmodulo",Y);
1166 : return NULL; /* LCOV_EXCL_LINE */
1167 : }
1168 49 : H = ZM_hnfall_i(shallowconcat(M,delta), &U, 1);
1169 49 : Y = hnf_solve(H,Y); if (!Y) return gen_0;
1170 42 : l = lg(H); /* may be smaller than lM if some moduli are 0 */
1171 42 : n = l-1;
1172 42 : m = lg(U)-1 - n;
1173 42 : u1 = cgetg(m+1,t_MAT);
1174 42 : u2 = cgetg(n+1,t_MAT);
1175 203 : for (j=1; j<=m; j++) { GEN c = gel(U,j); setlg(c,lM); gel(u1,j) = c; }
1176 42 : U += m;
1177 140 : for (j=1; j<=n; j++) { GEN c = gel(U,j); setlg(c,lM); gel(u2,j) = c; }
1178 : /* (u1 u2)
1179 : * (M D) (* * ) = (0 H) */
1180 42 : u1 = ZM_lll(u1, 0.75, LLL_INPLACE);
1181 42 : Y = ZM_ZC_mul(u2,Y);
1182 42 : x = ZC_reducemodmatrix(Y, u1);
1183 42 : if (!pU) x = gc_upto(av, x);
1184 : else
1185 : {
1186 14 : (void)gc_all(av, 2, &x, &u1);
1187 14 : *pU = u1;
1188 : }
1189 42 : return x;
1190 : }
1191 : /* to be used when one D_i is 0 */
1192 : static GEN
1193 77 : msolvemod0(GEN M, GEN D, GEN Y, long flag)
1194 : {
1195 77 : pari_sp av = avma;
1196 : GEN y, x, U;
1197 77 : if (!flag) return gaussmoduloall(M,D,Y,NULL);
1198 28 : y = cgetg(3,t_VEC);
1199 28 : x = gaussmoduloall(M,D,Y,&U);
1200 28 : if (x == gen_0) { set_avma(av); return gen_0; }
1201 21 : gel(y,1) = x;
1202 21 : gel(y,2) = U; return y;
1203 :
1204 : }
1205 : GEN
1206 44156 : matsolvemod(GEN M, GEN D, GEN Y, long flag)
1207 : {
1208 44156 : pari_sp av = avma;
1209 44156 : long m, n, i, char0 = 0;
1210 44156 : if (typ(M)!=t_MAT || !RgM_is_ZM(M)) pari_err_TYPE("matsolvemod (M)",M);
1211 44142 : RgM_dimensions(M,&m,&n);
1212 44142 : if (typ(D)!=t_INT && (typ(D)!=t_COL || !RgV_is_ZV(D)))
1213 28 : pari_err_TYPE("matsolvemod (D)",D);
1214 44114 : if (n)
1215 43967 : { if (typ(D)==t_COL && lg(D)!=m+1) pari_err_DIM("matsolvemod [1]"); }
1216 : else
1217 147 : { if (typ(D)==t_COL) m = lg(D)-1; }
1218 44107 : if (typ(Y)==t_INT)
1219 43890 : Y = const_col(m,Y);
1220 217 : else if (typ(Y)!=t_COL || !RgV_is_ZV(Y)) pari_err_TYPE("matsolvemod (Y)",Y);
1221 44079 : if (typ(D)==t_COL && typ(Y)==t_COL && lg(D)!=lg(Y))
1222 28 : pari_err_DIM("matsolvemod [2]");
1223 44051 : if (!n && !m) m = lg(Y)-1;
1224 43995 : else if (m != lg(Y)-1) pari_err_DIM("matsolvemod [3]");
1225 44044 : if (typ(D)==t_INT)
1226 : {
1227 43918 : if (signe(D)<0) pari_err_DOMAIN("matsolvemod","D","<",gen_0,D);
1228 43911 : if (!signe(D)) char0 = 1;
1229 : }
1230 : else /*typ(D)==t_COL*/
1231 364 : for (i=1; i<=m; i++)
1232 : {
1233 245 : if (signe(gel(D,i))<0)
1234 7 : pari_err_DOMAIN("matsolvemod","D[i]","<",gen_0,gel(D,i));
1235 238 : if (!signe(gel(D,i))) char0 = 1;
1236 : }
1237 87983 : return gc_GEN(av, char0? msolvemod0(M,D,Y,flag)
1238 43953 : : matsolvemod_finite(M,D,Y,flag));
1239 : }
1240 : GEN
1241 0 : gaussmodulo2(GEN M, GEN D, GEN Y) { return matsolvemod(M,D,Y, 1); }
1242 : GEN
1243 0 : gaussmodulo(GEN M, GEN D, GEN Y) { return matsolvemod(M,D,Y, 0); }
|
