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 26101097 : _Fp_add(void *data, GEN x, GEN y) { (void) data; return addii(x,y); }
69 :
70 : static GEN
71 4305178 : _Fp_neg(void *data, GEN x) { (void) data; return negi(x); }
72 :
73 : static GEN
74 37382729 : _Fp_mul(void *data, GEN x, GEN y) { (void) data; return mulii(x,y); }
75 :
76 : static GEN
77 2168283 : _Fp_rann(void *data, GEN x)
78 : {
79 2168283 : GEN d, N = (GEN)data;
80 2168283 : if (!signe(x)) return gen_1;
81 2021577 : d = gcdii(x,N);
82 2021383 : return modii(diviiexact(N,d),N);
83 : }
84 :
85 : static GEN
86 2207975 : _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 28 : Z_split(GEN D, GEN a)
91 : {
92 : long i, n;
93 : GEN N;
94 28 : n = expi(D);
95 28 : n = n<2 ? 1 : expu(n)+1;
96 196 : for (i=1;i<=n;i++)
97 168 : a = Fp_sqr(a,D);
98 28 : N = gcdii(a,D);
99 28 : return diviiexact(D,N);
100 : }
101 :
102 : /* c s.t. gcd(a+cb,N) = gcd(a,b,N) without factoring */
103 : static GEN
104 28 : Z_stab(GEN a, GEN b, GEN N)
105 : {
106 : GEN g, a2, N2;
107 28 : g = gcdii(a,b);
108 28 : g = gcdii(g,N);
109 28 : N2 = diviiexact(N,g);
110 28 : a2 = diviiexact(a,g);
111 28 : return Z_split(N2,a2);
112 : }
113 :
114 : static GEN
115 6249258 : _Fp_unit(void *data, GEN x)
116 : {
117 6249258 : GEN g,s,v,d,N=(GEN)data,N2;
118 : long i;
119 6249258 : if (!signe(x)) return NULL;
120 5817792 : g = bezout(x,N,&s,&v);
121 5818345 : if (equali1(g) || equali1(gcdii(s,N))) return mkvec2(g,s);
122 44591 : N2 = diviiexact(N,g);
123 61196 : for (i=0; i<5; i++)
124 : {
125 61168 : s = addii(s,N2);
126 61168 : if (equali1(gcdii(s,N))) return mkvec2(g,s);
127 : }
128 28 : d = Z_stab(s,N2,N);
129 28 : d = mulii(d,N2);
130 28 : v = Fp_add(s,d,N);
131 28 : if (equali1(v)) return NULL;
132 28 : return mkvec2(g,v);
133 : }
134 :
135 : static GEN
136 3203028 : _Fp_extgcd(void *data, GEN x, GEN y, int* smallop)
137 : {
138 : GEN d,u,v,m;
139 3203028 : if (equali1(y))
140 : {
141 601826 : *smallop = 1;
142 601826 : return mkvec2(y,mkmat2(
143 : mkcol2(gen_1,Fp_neg(x,(GEN)data)),
144 : mkcol2(gen_0,gen_1)));
145 : }
146 2601202 : *smallop = 0;
147 2601202 : d = bezout(x,y,&u,&v);
148 2601204 : if (!signe(d)) return mkvec2(d,matid(2));
149 2601204 : m = cgetg(3,t_MAT);
150 2601203 : m = mkmat2(
151 : mkcol2(diviiexact(y,d),negi(diviiexact(x,d))),
152 : mkcol2(u,v));
153 2601205 : return mkvec2(d,m);
154 : }
155 :
156 : static int
157 89700674 : _Fp_equal0(GEN x) { return !signe(x); }
158 :
159 : static int
160 26232743 : _Fp_equal1(GEN x) { return equali1(x); }
161 :
162 : static GEN
163 14535021 : _Fp_s(void *data, long x)
164 : {
165 14535021 : if (!x) return gen_0;
166 1899302 : if (x==1) return gen_1;
167 0 : return modsi(x,(GEN)data);
168 : }
169 :
170 : static GEN
171 42231192 : _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 843106 : get_Fp_hermite(void **data, GEN p)
179 : {
180 843106 : *data = (void*)p; return &Fp_hermite;
181 : }
182 :
183 : static void
184 17073165 : gen_redcol(GEN C, long lim, void* data, const struct bb_hermite *R)
185 : {
186 : long i;
187 61384958 : for (i=1; i<=lim; i++)
188 44314178 : if (!R->equal0(gel(C,i)))
189 23667880 : gel(C,i) = R->red(data, gel(C,i));
190 17070780 : }
191 :
192 : static GEN
193 : /* return NULL if a==0 */
194 : /* assume C*a is zero after lim */
195 17714041 : 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 17714041 : if (R->equal1(a)) return C;
200 10576534 : if (R->equal0(a)) return NULL;
201 7586814 : Ca = cgetg(lg(C),t_COL);
202 27933252 : for (i=1; i<=lim; i++)
203 20346456 : gel(Ca,i) = R->mul(data, gel(C,i), a);
204 7586796 : if (fillzeros && lim+1 < lg(C))
205 : {
206 5817231 : zero = R->s(data,0);
207 24870761 : for (i=lim+1; i<lg(C); i++)
208 19053531 : gel(Ca,i) = zero;
209 : }
210 7586795 : return Ca;
211 : }
212 :
213 : static void
214 : /* C1 <- C1 + C2 */
215 : /* assume C2[i]==0 for i>lim */
216 5098837 : gen_addcol(GEN C1, GEN C2, long lim, void* data, const struct bb_hermite *R)
217 : {
218 : long i;
219 18527468 : for (i=1; i<=lim; i++)
220 13428630 : gel(C1,i) = R->add(data, gel(C1,i), gel(C2,i));
221 5098838 : }
222 :
223 : static void
224 : /* H[,i] <- H[,i] + C*a */
225 : /* assume C is zero after lim */
226 3825390 : gen_addrightmul(GEN H, GEN C, GEN a, long i, long lim, void* data, const struct bb_hermite *R)
227 : {
228 : GEN Ca;
229 3825390 : if (R->equal0(a)) return;
230 1393736 : Ca = gen_rightmulcol(C, a, lim, 0, data, R);
231 1393741 : gen_addcol(gel(H,i), Ca, lim, data, R);
232 : }
233 :
234 : static GEN
235 2940471 : gen_zerocol(long n, void* data, const struct bb_hermite *R)
236 : {
237 2940471 : GEN C = cgetg(n+1,t_COL), zero = R->s(data, 0);
238 : long i;
239 13973395 : for (i=1; i<=n; i++) gel(C,i) = zero;
240 2940444 : return C;
241 : }
242 :
243 : static GEN
244 1556438 : gen_zeromat(long m, long n, void* data, const struct bb_hermite *R)
245 : {
246 1556438 : GEN M = cgetg(n+1,t_MAT);
247 : long i;
248 4353353 : for (i=1; i<=n; i++) gel(M,i) = gen_zerocol(m, data, R);
249 1556369 : return M;
250 : }
251 :
252 : static GEN
253 365253 : gen_colei(long n, long i, void* data, const struct bb_hermite *R)
254 : {
255 365253 : GEN C = cgetg(n+1,t_COL), zero = R->s(data, 0);
256 : long j;
257 2494233 : for (j=1; j<=n; j++)
258 2128980 : if (i!=j) gel(C,j) = zero;
259 365253 : else gel(C,j) = R->s(data,1);
260 365253 : 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 2570460 : gen_matmul_hermite(GEN A, GEN B, void* data, const struct bb_hermite *R)
274 : {
275 2570460 : GEN M,sum,prod,zero = R->s(data,0);
276 : long a,b,c,c2,i,j,k;
277 2570459 : RgM_dimensions(A,&a,&c);
278 2570461 : RgM_dimensions(B,&c2,&b);
279 2570462 : if (c!=c2) pari_err_DIM("gen_matmul_hermite");
280 2570462 : M = cgetg(b+1,t_MAT);
281 5474920 : for (j=1; j<=b; j++)
282 : {
283 2904525 : gel(M,j) = cgetg(a+1,t_COL);
284 8169170 : for (i=1; i<=a; i++)
285 : {
286 5264704 : sum = zero;
287 16849874 : for (k=1; k<=c; k++)
288 : {
289 11585224 : prod = R->mul(data, gcoeff(A,i,k), gcoeff(B,k,j));
290 11585001 : sum = R->add(data, sum, prod);
291 : }
292 5264650 : gcoeff(M,i,j) = sum;
293 : }
294 2904466 : gen_redcol(gel(M,j), a, data, R);
295 : }
296 2570395 : 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 6583720 : 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 6583720 : Au1 = gen_rightmulcol(A, gel(U,1), lim, 1, data, R);
306 6583737 : Bu2 = gen_rightmulcol(B, gel(U,2), lim, 1, data, R);
307 6583727 : if (!Au1 && !Bu2) { *C = gen_zerocol(lg(A)-1, data, R); return; }
308 6583727 : if (!Au1) { *C = Bu2; return; }
309 4080023 : if (!Bu2) { *C = Au1; return; }
310 3704972 : gen_addcol(Au1, Bu2, lim, data, R);
311 3704973 : *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 3291868 : 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 3291868 : Hi = shallowcopy(gel(H,i));
321 3291868 : Hj = shallowcopy(gel(H,j));
322 3291868 : gen_rightlincomb(Hi, Hj, gel(U,1), &gel(H,i), lim, data, R);
323 3291858 : gen_rightlincomb(Hi, Hj, gel(U,2), &gel(H,j), lim, data, R);
324 3291868 : }
325 :
326 : static int
327 : /* assume C is zero after lim */
328 4094595 : gen_is_zerocol(GEN C, long lim, void* data, const struct bb_hermite *R)
329 : {
330 : long i;
331 : (void) data;
332 7357126 : for (i=1; i<=lim; i++)
333 5752990 : if (!R->equal0(gel(C,i))) return 0;
334 1604136 : 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 2767302 : mkoptransv(long i, long j, GEN a, void* data, const struct bb_hermite *R)
342 : {
343 2767302 : a = R->red(data,a);
344 2767287 : if (R->equal0(a)) return NULL;
345 1077801 : return mkvec2(mkvecsmall2(i,j),a);
346 : }
347 :
348 : static GEN
349 : /* (Ci|Cj) <- (Ci|Cj)*U */
350 924587 : mkopU(long i, long j, GEN U, void* data, const struct bb_hermite *R)
351 : {
352 924587 : if (R->equal1(gcoeff(U,1,1)) && R->equal0(gcoeff(U,1,2))
353 403263 : && R->equal1(gcoeff(U,2,2))) return mkoptransv(i,j,gcoeff(U,2,1),data,R);
354 521324 : return mkvec2(mkvecsmall3(i,j,0),U);
355 : }
356 :
357 : static GEN
358 : /* Ci <- Ci*u */
359 1875163 : mkopmul(long i, GEN u, const struct bb_hermite *R)
360 : {
361 1875163 : if (R->equal1(u)) return NULL;
362 191476 : 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 11421740 : gen_leftapply(GEN C, GEN op, void* data, const struct bb_hermite *R)
418 : {
419 : GEN C2, ind, X;
420 : long i, j;
421 11421740 : switch (typ(op))
422 : {
423 578788 : case t_VECSMALL:
424 578788 : C2 = vecpermute(C,perm_inv(op));
425 4729627 : for (i=1; i<lg(C); i++) gel(C,i) = gel(C2,i);
426 578788 : return;
427 10842951 : case t_VEC:
428 10842951 : ind = gel(op,1);
429 10842951 : switch (lg(op))
430 : {
431 169134 : case 2:
432 169134 : swap(gel(C,ind[1]),gel(C,ind[2]));
433 169134 : return;
434 10673817 : case 3:
435 10673817 : X = gel(op,2);
436 10673817 : i = ind[1];
437 10673817 : switch (lg(ind))
438 : {
439 700717 : case 2:
440 700717 : gel(C,i) = R->mul(data, X, gel(C,i));
441 700713 : gel(C,i) = R->red(data, gel(C,i));
442 700709 : return;
443 7612878 : case 3:
444 7612878 : j = ind[2];
445 7612878 : if (R->equal0(gel(C,i))) return;
446 1087427 : gel(C,j) = R->add(data, gel(C,j), R->mul(data, X, gel(C,i)));
447 1087428 : return;
448 2360258 : case 4:
449 2360258 : j = ind[2];
450 2360258 : C2 = gen_matmul_hermite(X, mkmat(mkcol2(gel(C,i),gel(C,j))), data, R);
451 2360199 : gel(C,i) = gcoeff(C2,1,1);
452 2360199 : gel(C,j) = gcoeff(C2,2,1);
453 2360199 : 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 : }
503 42 : return d;
504 : }
505 :
506 : static int
507 211589 : gen_is_inv(GEN x, void* data, const struct bb_hermite *R)
508 : {
509 211589 : GEN u = R->unit(data, x);
510 211589 : if (!u) return R->equal1(x);
511 72219 : return R->equal1(gel(u,1));
512 : }
513 :
514 : static long
515 195510 : gen_last_inv_diago(GEN A, void* data, const struct bb_hermite *R)
516 : {
517 : long i,m,n,j;
518 195510 : RgM_dimensions(A,&m,&n);
519 226009 : for (i=1,j=n-m+1; i<=m; i++,j++)
520 211589 : if (!gen_is_inv(gcoeff(A,i,j),data,R)) return i-1;
521 14420 : 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 942715 : 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 942715 : pari_sp av = avma, av1;
531 942715 : 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 942716 : long m,n,i,j,s,si,i2,si2,nbz,lim,extra,maxop=0,nbop=0,lastinv=0;
533 : int smallop;
534 :
535 942716 : av1 = avma;
536 :
537 942716 : RgM_dimensions(A,&m,&n);
538 942717 : if (early_abort && n<m)
539 : {
540 14000 : if (ops) *ops = zero;
541 14000 : return NULL;
542 : }
543 928717 : if (n<m+1)
544 : {
545 828358 : extra = m+1-n;
546 828358 : H = shallowmatconcat(mkvec2(gen_zeromat(m,extra,data,R),A));
547 : }
548 : else
549 : {
550 100359 : extra = 0;
551 100359 : H = RgM_shallowcopy(A);
552 : }
553 928720 : RgM_dimensions(H,&m,&n);
554 928719 : s = n-m; /* shift */
555 :
556 928719 : if(ops)
557 : {
558 733207 : maxop = m*n + (m*(m+1))/2 + 1;
559 733207 : *ops = zerovec(maxop); /* filled with placeholders so gerepile can handle it */
560 : }
561 :
562 : /* put in triangular form */
563 3572880 : for (i=m,si=s+m; i>0 && si>extra; i--,si--) /* si = s+i */
564 : {
565 2649184 : if (R->red) gcoeff(H,i,si) = R->red(data, gcoeff(H,i,si));
566 : /* bottom-right diagonal */
567 9249530 : for (j = extra+1; j < si; j++)
568 : {
569 6600721 : if (R->red) gcoeff(H,i,j) = R->red(data, gcoeff(H,i,j));
570 6600669 : if (R->equal0(gcoeff(H,i,j))) continue;
571 3055331 : U = R->extgcd(data, gcoeff(H,i,j), gcoeff(H,i,si), &smallop);
572 3055401 : d = gel(U,1);
573 3055401 : U = gel(U,2);
574 3055401 : if (n>10)
575 : {
576 : /* normalize diagonal coefficient -> faster reductions on this row */
577 1831409 : u = R->unit(data, d);
578 1831409 : if (u)
579 : {
580 1831409 : g = gel(u,1);
581 1831409 : u = gel(u,2);
582 1831409 : gcoeff(U,1,2) = R->mul(data, gcoeff(U,1,2), u);
583 1831409 : gcoeff(U,2,2) = R->mul(data, gcoeff(U,2,2), u);
584 1831409 : d = g;
585 : }
586 : }
587 3055401 : gen_elem(H, U, j, si, i-1, data, R);
588 3055401 : if (ops)
589 : {
590 888544 : op = mkopU(j,si,U,data,R);
591 888544 : if (op) { nbop++; gel(*ops, nbop) = op; }
592 : }
593 3055401 : gcoeff(H,i,j) = zero;
594 3055401 : gcoeff(H,i,si) = d;
595 3055401 : if (R->red && !smallop)
596 : {
597 2463739 : gen_redcol(gel(H,si), i-1, data, R);
598 2463737 : gen_redcol(gel(H,j), i-1, data, R);
599 : }
600 : }
601 :
602 2648809 : if (early_abort)
603 : {
604 1455376 : d = gcoeff(H,i,si);
605 1455376 : u = R->unit(data, d);
606 1455593 : if (u) d = gel(u,1);
607 1455593 : if (!R->equal1(d))
608 : {
609 4914 : if (ops) *ops = d;
610 4914 : return NULL;
611 : }
612 : }
613 :
614 2644112 : 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 923696 : if (!ops)
622 195510 : lastinv = gen_last_inv_diago(H, data, R);
623 :
624 : /* put in reduced Howell form */
625 923696 : if (!only_triangular)
626 : {
627 3705569 : for (i=m,si=s+m; i>0; i--,si--) /* si = s+i */
628 : {
629 : /* normalize diagonal coefficient */
630 2781865 : if (i<=lastinv) /* lastinv>0 => !ops */
631 30499 : gcoeff(H,i,si) = one;
632 : else
633 : {
634 2751366 : u = R->unit(data,gcoeff(H,i,si));
635 2751571 : if (u)
636 : {
637 2459618 : g = gel(u,1);
638 2459618 : u = gel(u,2);
639 2459618 : gel(H,si) = gen_rightmulcol(gel(H,si), u, i-1, 1, data, R);
640 2459627 : gcoeff(H,i,si) = g;
641 2459627 : if (R->red) gen_redcol(gel(H,si), i-1, data, R);
642 2459632 : if (ops)
643 : {
644 1875161 : op = mkopmul(si,u,R);
645 1875168 : if (op) { nbop++; gel(*ops,nbop) = op; }
646 : }
647 : }
648 291953 : else if (R->red) gcoeff(H,i,si) = R->red(data, gcoeff(H,i,si));
649 : }
650 2782091 : piv = gcoeff(H,i,si);
651 :
652 : /* reduce above diagonal */
653 2782091 : if (!R->equal0(piv))
654 : {
655 2490136 : C = gel(H,si);
656 6166799 : for (j=si+1; j<=n; j++)
657 : {
658 3676662 : if (i<=lastinv) /* lastinv>0 => !ops */
659 26978 : gcoeff(H,i,j) = zero;
660 : else
661 : {
662 3649684 : gcoeff(H,i,j) = R->red(data, gcoeff(H,i,j));
663 3649682 : if (R->equal1(piv)) { q = gcoeff(H,i,j); r = zero; }
664 1552572 : else q = R->lquo(data, gcoeff(H,i,j), piv, &r);
665 3649673 : q = R->neg(data,q);
666 3649689 : gen_addrightmul(H, C, q, j, i-1, data, R);
667 3649696 : if (ops)
668 : {
669 2215557 : op = mkoptransv(j,si,q,data,R);
670 2215546 : if (op) { nbop++; gel(*ops,nbop) = op; }
671 : }
672 3649685 : gcoeff(H,i,j) = r;
673 : }
674 : }
675 : }
676 :
677 : /* ensure Howell property */
678 2782085 : if (i>1)
679 : {
680 1858541 : a = R->rann(data, piv);
681 1858271 : if (!R->equal0(a))
682 : {
683 544390 : gel(H,1) = gen_rightmulcol(gel(H,si), a, i-1, 1, data, R);
684 544390 : if (gel(H,1) == gel(H,si)) gel(H,1) = shallowcopy(gel(H,1)); /* in case rightmulcol cheated */
685 544390 : if (ops)
686 : {
687 148484 : op = mkoptransv(1,si,a,data,R);
688 148484 : if (op) { nbop++; gel(*ops,nbop) = op; }
689 : }
690 2174746 : for (i2=i-1,si2=s+i2; i2>0; i2--,si2--)
691 : {
692 1630356 : if (R->red) gcoeff(H,i2,1) = R->red(data, gcoeff(H,i2,1));
693 1630356 : if (R->equal0(gcoeff(H,i2,1))) continue;
694 278054 : if (R->red) gcoeff(H,i2,si2) = R->red(data, gcoeff(H,i2,si2));
695 278054 : if (R->equal0(gcoeff(H,i2,si2)))
696 : {
697 130424 : swap(gel(H,1), gel(H,si2));
698 130424 : if (ops) { nbop++; gel(*ops,nbop) = mkopswap(1,si2); }
699 130424 : continue;
700 : }
701 147630 : U = R->extgcd(data, gcoeff(H,i2,1), gcoeff(H,i2,si2), &smallop);
702 147630 : d = gel(U,1);
703 147630 : U = gel(U,2);
704 147630 : gen_elem(H, U, 1, si2, i2-1, data, R);
705 147630 : if (ops)
706 : {
707 36043 : op = mkopU(1,si2,U,data,R);
708 36043 : if (op) { nbop++; gel(*ops,nbop) = op; }
709 : }
710 147630 : gcoeff(H,i2,1) = zero;
711 147630 : gcoeff(H,i2,si2) = d;
712 147630 : if (R->red && !smallop)
713 : {
714 137466 : gen_redcol(gel(H,si2), i2, data, R);
715 137466 : gen_redcol(gel(H,1), i2-1, data, R);
716 : }
717 : }
718 : }
719 : }
720 :
721 2781878 : 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 923671 : if (R->red)
730 4974010 : for (j=1; j<=n; j++)
731 : {
732 4050415 : lim = maxss(0,m-n+j);
733 4050396 : gen_redcol(gel(H,j), lim, data, R);
734 13510459 : for (i=lim+1; i<=m; i++) gcoeff(H,i,j) = zero;
735 : }
736 :
737 : /* put zero columns first */
738 923478 : iszero = cgetg(n+1,t_VECSMALL);
739 :
740 923745 : nbz = 0;
741 4974384 : for (i=1; i<=n; i++)
742 : {
743 4050580 : iszero[i] = gen_is_zerocol(gel(H,i), maxss(0,m-n+i), data, R);
744 4050639 : if (iszero[i]) nbz++;
745 : }
746 :
747 923804 : j = 1;
748 923804 : if (permute_zerocols)
749 : {
750 155729 : perm = cgetg(n+1, t_VECSMALL);
751 900963 : for (i=1; i<=n; i++)
752 745234 : if (iszero[i])
753 : {
754 314321 : perm[j] = i;
755 314321 : j++;
756 : }
757 : }
758 768075 : else perm = cgetg(n-nbz+1, t_VECSMALL);
759 4974556 : for (i=1; i<=n; i++)
760 4050761 : if (!iszero[i])
761 : {
762 2490456 : perm[j] = i;
763 2490456 : j++;
764 : }
765 :
766 923795 : if (permute_zerocols || remove_zerocols==2) H = vecpermute(H, perm);
767 923803 : if (permute_zerocols && remove_zerocols==2) H = vecslice(H, nbz+1, n);
768 923803 : if (remove_zerocols==1) H = vecslice(H, s+1, n);
769 923803 : if (permute_zerocols && ops) { nbop++; gel(*ops,nbop) = perm; }
770 :
771 923803 : if (ops) { setlg(*ops, nbop+1); } /* should have nbop <= maxop */
772 :
773 923797 : return H;
774 : }
775 :
776 : static GEN
777 95963 : 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 95963 : pari_sp av = avma;
780 95963 : GEN H = gen_howell_i(A, remove_zerocols, permute_zerocols, early_abort, only_triangular, ops, data, R);
781 95963 : return gc_all(av, ops?2:1, &H, ops);
782 : }
783 :
784 : static GEN
785 151977 : gen_matimage(GEN A, GEN* U, void *data, const struct bb_hermite *R)
786 : {
787 : GEN ops, H;
788 151977 : 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 95963 : 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 99547 : 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 99547 : RgM_dimensions(H,&m,&n);
821 99547 : K = gen_zeromat(n, n, data, R);
822 409290 : for (j=n,i=m; j>0; j--)
823 : {
824 523859 : while (R->equal0(gcoeff(H,i,j))) i--;
825 309743 : piv = gcoeff(H,i,j);
826 309743 : if (R->equal0(piv)) continue;
827 309743 : gcoeff(K,j,j) = R->rann(data, piv);
828 309743 : if (j<n)
829 : {
830 210196 : FK = gen_matmul_hermite(matslice(H,i,i,j+1,n), matslice(K, j+1, n, j+1, n), data, R);
831 754453 : for (j2=j+1; j2<=n; j2++)
832 544257 : 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 99547 : 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 99617 : 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 99617 : RgM_dimensions(H,&m,&r);
847 99617 : if (!r) return gen_matid_hermite(n, data, R); /* zerology: what if 0==1 in R? */
848 99547 : n2 = maxss(n,m+1);
849 99547 : extra = n2-n;
850 99547 : nbz = n2-r;
851 : /* compute kernel of augmented matrix */
852 99547 : KH = gen_kernel_howell(H, data, R);
853 99547 : zC = gen_zerocol(nbz, data, R);
854 99547 : K = cgetg(nbz+r+1, t_MAT);
855 324688 : for (i=1; i<=nbz; i++)
856 225141 : gel(K,i) = gen_colei(nbz+r, i, data, R);
857 409290 : for (i=1; i<=r; i++)
858 309743 : gel(K,nbz+i) = shallowconcat(zC, gel(KH,i));
859 634431 : for (i=1; i<lg(K); i++)
860 : {
861 534884 : av = avma;
862 9913666 : for (o=lg(ops)-1; o>0; o--)
863 9378782 : gen_leftapply(gel(K,i), gel(ops,o), data, R);
864 534884 : gen_redcol(gel(K,i), nbz+r, data, R);
865 534884 : gerepileall(av, 1, &gel(K,i));
866 : }
867 : /* deduce kernel of original matrix */
868 99547 : K = rowpermute(K, cyclic_perm(n2,extra));
869 99547 : K = gen_howell_i(K, 2, 0, 0, 0, NULL, data, R);
870 226338 : for (j=lg(K)-1, i=n2; j>0; j--)
871 : {
872 300734 : while (R->equal0(gcoeff(K,i,j))) i--;
873 197694 : if (i<=n) return matslice(K, 1, n, 1, j);
874 : }
875 28644 : return cgetg(1,t_MAT);
876 : }
877 :
878 : /* not GC-clean */
879 : static GEN
880 55762 : gen_kernel(GEN A, GEN* im, void *data, const struct bb_hermite *R)
881 : {
882 55762 : pari_sp av = avma;
883 55762 : long n = lg(A)-1;
884 : GEN H, ops, K;
885 55762 : H = gen_howell_i(A, 2, 1, 0, 0, &ops, data, R);
886 55762 : gerepileall(av,2,&H,&ops);
887 55762 : K = gen_kernel_from_howell(H, ops, n, data, R);
888 55762 : if (im) *im = H;
889 55762 : return K;
890 : }
891 :
892 : /* right inverse, not GC-clean */
893 : static GEN
894 591434 : gen_inv(GEN A, void* data, const struct bb_hermite *R)
895 : {
896 : pari_sp av;
897 591434 : GEN ops, H, U, un=R->s(data,1);
898 : long m,n,j,o,n2;
899 591434 : RgM_dimensions(A,&m,&n);
900 591434 : av = avma;
901 591434 : H = gen_howell_i(A, 0, 0, 1, 0, &ops, data, R);
902 591430 : if (!H) pari_err_INV("gen_inv", ops);
903 572516 : n2 = lg(H)-1;
904 572516 : ops = gerepilecopy(av,ops); /* get rid of H */
905 572536 : U = gen_zeromat(n2, m, data, R);
906 2016936 : for (j=1; j<=m; j++)
907 1444419 : gcoeff(U,j+n2-m,j) = un;
908 2016968 : for (j=1; j<=m; j++)
909 : {
910 1444427 : av = avma;
911 3089067 : for (o=lg(ops)-1; o>0; o--)
912 1644748 : gen_leftapply(gel(U,j), gel(ops,o), data, R);
913 1444319 : gen_redcol(gel(U,j), n2, data, R);
914 1444035 : gerepileall(av, 1, &gel(U,j));
915 : }
916 572541 : if (n2>n) U = rowslice(U, n2-n+1, n2);
917 572522 : 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 152012 : matimagemod(GEN A, GEN d, GEN* U)
986 : {
987 : void *data;
988 : const struct bb_hermite* R;
989 152012 : if (typ(A)!=t_MAT || !RgM_is_ZM(A)) pari_err_TYPE("matimagemod", A);
990 151998 : if (typ(d)!=t_INT) pari_err_TYPE("matimagemod", d);
991 151991 : if (signe(d)<=0) pari_err_DOMAIN("matimagemod", "d", "<=", gen_0, d);
992 151984 : if (equali1(d)) return cgetg(1,t_MAT);
993 151977 : R = get_Fp_hermite(&data, d);
994 151977 : 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 55797 : matkermod(GEN A, GEN d, GEN* im)
1040 : {
1041 55797 : pari_sp av = avma;
1042 : void *data;
1043 : const struct bb_hermite* R;
1044 : long m,n;
1045 : GEN K;
1046 55797 : if (typ(A)!=t_MAT || !RgM_is_ZM(A)) pari_err_TYPE("matkermod", A);
1047 55783 : if (typ(d)!=t_INT) pari_err_TYPE("matkermod", d);
1048 55776 : if (signe(d)<=0) pari_err_DOMAIN("makermod", "d", "<=", gen_0, d);
1049 55769 : if (equali1(d)) return cgetg(1,t_MAT);
1050 55762 : R = get_Fp_hermite(&data, d);
1051 55762 : RgM_dimensions(A,&m,&n);
1052 55762 : if (!im && m>2*n) /* TODO tune treshold */
1053 4592 : A = shallowtrans(matimagemod(shallowtrans(A),d,NULL));
1054 55762 : K = gen_kernel(A, im, data, R); return gc_all(av,im?2:1,&K,im);
1055 : }
1056 :
1057 : /* left inverse */
1058 : GEN
1059 677661 : matinvmod(GEN A, GEN d)
1060 : {
1061 677661 : pari_sp av = avma;
1062 : void *data;
1063 : const struct bb_hermite* R;
1064 : GEN U;
1065 677661 : if (typ(A)!=t_MAT || !RgM_is_ZM(A)) pari_err_TYPE("matinvmod", A);
1066 677668 : if (typ(d)!=t_INT) pari_err_TYPE("matinvmod", d);
1067 677661 : if (signe(d)<=0) pari_err_DOMAIN("matinvmod", "d", "<=", gen_0, d);
1068 677654 : if (equali1(d)) {
1069 : long m,n;
1070 86315 : RgM_dimensions(A,&m,&n);
1071 86315 : if (m<n) pari_err_INV("matinvmod",A);
1072 86308 : return zeromatcopy(n,m);
1073 : }
1074 591347 : R = get_Fp_hermite(&data, d);
1075 591373 : U = gen_inv(shallowtrans(A), data, R);
1076 572522 : 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 44149 : matsolvemod(GEN M, GEN D, GEN Y, long flag)
1210 : {
1211 44149 : pari_sp av = avma;
1212 44149 : long m, n, i, char0 = 0;
1213 44149 : if (typ(M)!=t_MAT || !RgM_is_ZM(M)) pari_err_TYPE("matsolvemod (M)",M);
1214 44135 : RgM_dimensions(M,&m,&n);
1215 44135 : if (typ(D)!=t_INT && (typ(D)!=t_COL || !RgV_is_ZV(D)))
1216 28 : pari_err_TYPE("matsolvemod (D)",D);
1217 44107 : if (n)
1218 43960 : { if (typ(D)==t_COL && lg(D)!=m+1) pari_err_DIM("matsolvemod [1]"); }
1219 : else
1220 147 : { if (typ(D)==t_COL) m = lg(D)-1; }
1221 44100 : if (typ(Y)==t_INT)
1222 43890 : Y = const_col(m,Y);
1223 210 : else if (typ(Y)!=t_COL || !RgV_is_ZV(Y)) pari_err_TYPE("matsolvemod (Y)",Y);
1224 44072 : if (typ(D)==t_COL && typ(Y)==t_COL && lg(D)!=lg(Y))
1225 28 : pari_err_DIM("matsolvemod [2]");
1226 44044 : if (!n && !m) m = lg(Y)-1;
1227 43988 : else if (m != lg(Y)-1) pari_err_DIM("matsolvemod [3]");
1228 44037 : if (typ(D)==t_INT)
1229 : {
1230 43918 : if (signe(D)<0) pari_err_DOMAIN("matsolvemod","D","<",gen_0,D);
1231 43911 : if (!signe(D)) char0 = 1;
1232 : }
1233 : else /*typ(D)==t_COL*/
1234 301 : for (i=1; i<=m; i++)
1235 : {
1236 189 : if (signe(gel(D,i))<0)
1237 7 : pari_err_DOMAIN("matsolvemod","D[i]","<",gen_0,gel(D,i));
1238 182 : if (!signe(gel(D,i))) char0 = 1;
1239 : }
1240 87976 : return gerepilecopy(av, char0? msolvemod0(M,D,Y,flag)
1241 43953 : : matsolvemod_finite(M,D,Y,flag));
1242 : }
1243 : GEN
1244 0 : gaussmodulo2(GEN M, GEN D, GEN Y) { return matsolvemod(M,D,Y, 1); }
1245 : GEN
1246 0 : gaussmodulo(GEN M, GEN D, GEN Y) { return matsolvemod(M,D,Y, 0); }
|
