1 /* Copyright (C) 2000 The PARI group.
2
3 This file is part of the PARI/GP package.
4
5 PARI/GP is free software; you can redistribute it and/or modify it under the
6 terms of the GNU General Public License as published by the Free Software
7 Foundation; either version 2 of the License, or (at your option) any later
8 version. It is distributed in the hope that it will be useful, but WITHOUT
9 ANY WARRANTY WHATSOEVER.
10
11 Check the License for details. You should have received a copy of it, along
12 with the package; see the file 'COPYING'. If not, write to the Free Software
13 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14
15 #include "pari.h"
16 #include "paripriv.h"
17
18 int
RgM_is_ZM(GEN x)19 RgM_is_ZM(GEN x)
20 {
21 long i, j, h, l = lg(x);
22 if (l == 1) return 1;
23 h = lgcols(x);
24 if (h == 1) return 1;
25 for (j = l-1; j > 0; j--)
26 for (i = h-1; i > 0; i--)
27 if (typ(gcoeff(x,i,j)) != t_INT) return 0;
28 return 1;
29 }
30
31 int
RgM_is_QM(GEN x)32 RgM_is_QM(GEN x)
33 {
34 long i, j, h, l = lg(x);
35 if (l == 1) return 1;
36 h = lgcols(x);
37 if (h == 1) return 1;
38 for (j = l-1; j > 0; j--)
39 for (i = h-1; i > 0; i--)
40 if (!is_rational_t(typ(gcoeff(x,i,j)))) return 0;
41 return 1;
42 }
43
44 int
RgV_is_ZMV(GEN V)45 RgV_is_ZMV(GEN V)
46 {
47 long i, l = lg(V);
48 for (i=1; i<l; i++)
49 if (typ(gel(V,i))!=t_MAT || !RgM_is_ZM(gel(V,i)))
50 return 0;
51 return 1;
52 }
53
54 /********************************************************************/
55 /** **/
56 /** GENERIC LINEAR ALGEBRA **/
57 /** **/
58 /********************************************************************/
59 /* GENERIC MULTIPLICATION involving zc/zm */
60
61 /* x[i,] * y */
62 static GEN
RgMrow_zc_mul_i(GEN x,GEN y,long c,long i)63 RgMrow_zc_mul_i(GEN x, GEN y, long c, long i)
64 {
65 pari_sp av = avma;
66 GEN s = NULL;
67 long j;
68 for (j=1; j<c; j++)
69 {
70 long t = y[j];
71 if (!t) continue;
72 if (!s) { s = gmulgs(gcoeff(x,i,j),t); continue; }
73 switch(t)
74 {
75 case 1: s = gadd(s, gcoeff(x,i,j)); break;
76 case -1: s = gsub(s, gcoeff(x,i,j)); break;
77 default: s = gadd(s, gmulgs(gcoeff(x,i,j), t)); break;
78 }
79 }
80 return s? gerepileupto(av, s): gc_const(av, gen_0);
81 }
82 GEN
RgMrow_zc_mul(GEN x,GEN y,long i)83 RgMrow_zc_mul(GEN x, GEN y, long i) { return RgMrow_zc_mul_i(x,y,lg(y),i); }
84 /* x nonempty t_MAT, y a compatible zc (dimension > 0). */
85 static GEN
RgM_zc_mul_i(GEN x,GEN y,long c,long l)86 RgM_zc_mul_i(GEN x, GEN y, long c, long l)
87 {
88 GEN z = cgetg(l,t_COL);
89 long i;
90 for (i = 1; i < l; i++) gel(z,i) = RgMrow_zc_mul_i(x,y,c,i);
91 return z;
92 }
93 GEN
RgM_zc_mul(GEN x,GEN y)94 RgM_zc_mul(GEN x, GEN y) { return RgM_zc_mul_i(x,y, lg(x), lgcols(x)); }
95 /* x t_MAT, y a compatible zm (dimension > 0). */
96 GEN
RgM_zm_mul(GEN x,GEN y)97 RgM_zm_mul(GEN x, GEN y)
98 {
99 long j, c, l = lg(x), ly = lg(y);
100 GEN z = cgetg(ly, t_MAT);
101 if (l == 1) return z;
102 c = lgcols(x);
103 for (j = 1; j < ly; j++) gel(z,j) = RgM_zc_mul_i(x, gel(y,j), l,c);
104 return z;
105 }
106
107 /* x[i,]*y, l = lg(y) > 1 */
108 static GEN
RgMrow_ZC_mul_i(GEN x,GEN y,long i,long l)109 RgMrow_ZC_mul_i(GEN x, GEN y, long i, long l)
110 {
111 pari_sp av = avma;
112 GEN t = gmul(gcoeff(x,i,1), gel(y,1)); /* l > 1 ! */
113 long j;
114 for (j=2; j<l; j++)
115 if (signe(gel(y,j))) t = gadd(t, gmul(gcoeff(x,i,j), gel(y,j)));
116 return gerepileupto(av,t);
117 }
118
119 /* compatible t_MAT * t_COL, lx = lg(x) = lg(y) > 1, l = lgcols(x) */
120 static GEN
RgM_ZC_mul_i(GEN x,GEN y,long lx,long l)121 RgM_ZC_mul_i(GEN x, GEN y, long lx, long l)
122 {
123 GEN z = cgetg(l,t_COL);
124 long i;
125 for (i=1; i<l; i++) gel(z,i) = RgMrow_ZC_mul_i(x,y,i,lx);
126 return z;
127 }
128
129 /* mostly useful when y is sparse */
130 GEN
RgM_ZM_mul(GEN x,GEN y)131 RgM_ZM_mul(GEN x, GEN y)
132 {
133 long j, c, l = lg(x), ly = lg(y);
134 GEN z = cgetg(ly, t_MAT);
135 if (l == 1) return z;
136 c = lgcols(x);
137 for (j = 1; j < ly; j++) gel(z,j) = RgM_ZC_mul_i(x, gel(y,j), l,c);
138 return z;
139 }
140
141 static GEN
RgV_zc_mul_i(GEN x,GEN y,long l)142 RgV_zc_mul_i(GEN x, GEN y, long l)
143 {
144 long i;
145 GEN z = gen_0;
146 pari_sp av = avma;
147 for (i = 1; i < l; i++) z = gadd(z, gmulgs(gel(x,i), y[i]));
148 return gerepileupto(av, z);
149 }
150 GEN
RgV_zc_mul(GEN x,GEN y)151 RgV_zc_mul(GEN x, GEN y) { return RgV_zc_mul_i(x, y, lg(x)); }
152
153 GEN
RgV_zm_mul(GEN x,GEN y)154 RgV_zm_mul(GEN x, GEN y)
155 {
156 long j, l = lg(x), ly = lg(y);
157 GEN z = cgetg(ly, t_VEC);
158 for (j = 1; j < ly; j++) gel(z,j) = RgV_zc_mul_i(x, gel(y,j), l);
159 return z;
160 }
161
162 /* scalar product x.x */
163 GEN
RgV_dotsquare(GEN x)164 RgV_dotsquare(GEN x)
165 {
166 long i, lx = lg(x);
167 pari_sp av = avma;
168 GEN z;
169 if (lx == 1) return gen_0;
170 z = gsqr(gel(x,1));
171 for (i=2; i<lx; i++)
172 {
173 z = gadd(z, gsqr(gel(x,i)));
174 if (gc_needed(av,3))
175 {
176 if(DEBUGMEM>1) pari_warn(warnmem,"RgV_dotsquare, i = %ld",i);
177 z = gerepileupto(av, z);
178 }
179 }
180 return gerepileupto(av,z);
181 }
182
183 /* scalar product x.y, lx = lg(x) = lg(y) */
184 static GEN
RgV_dotproduct_i(GEN x,GEN y,long lx)185 RgV_dotproduct_i(GEN x, GEN y, long lx)
186 {
187 pari_sp av = avma;
188 long i;
189 GEN z;
190 if (lx == 1) return gen_0;
191 z = gmul(gel(x,1),gel(y,1));
192 for (i=2; i<lx; i++)
193 {
194 z = gadd(z, gmul(gel(x,i), gel(y,i)));
195 if (gc_needed(av,3))
196 {
197 if(DEBUGMEM>1) pari_warn(warnmem,"RgV_dotproduct, i = %ld",i);
198 z = gerepileupto(av, z);
199 }
200 }
201 return gerepileupto(av,z);
202 }
203 GEN
RgV_dotproduct(GEN x,GEN y)204 RgV_dotproduct(GEN x,GEN y)
205 {
206 if (x == y) return RgV_dotsquare(x);
207 return RgV_dotproduct_i(x, y, lg(x));
208 }
209 /* v[1] + ... + v[lg(v)-1] */
210 GEN
RgV_sum(GEN v)211 RgV_sum(GEN v)
212 {
213 GEN p;
214 long i, l = lg(v);
215 if (l == 1) return gen_0;
216 p = gel(v,1); for (i=2; i<l; i++) p = gadd(p, gel(v,i));
217 return p;
218 }
219 /* v[1] + ... + v[n]. Assume lg(v) > n. */
220 GEN
RgV_sumpart(GEN v,long n)221 RgV_sumpart(GEN v, long n)
222 {
223 GEN p;
224 long i;
225 if (!n) return gen_0;
226 p = gel(v,1); for (i=2; i<=n; i++) p = gadd(p, gel(v,i));
227 return p;
228 }
229 /* v[m] + ... + v[n]. Assume lg(v) > n, m > 0. */
230 GEN
RgV_sumpart2(GEN v,long m,long n)231 RgV_sumpart2(GEN v, long m, long n)
232 {
233 GEN p;
234 long i;
235 if (n < m) return gen_0;
236 p = gel(v,m); for (i=m+1; i<=n; i++) p = gadd(p, gel(v,i));
237 return p;
238 }
239 GEN
RgM_sumcol(GEN A)240 RgM_sumcol(GEN A)
241 {
242 long i,j,m,l = lg(A);
243 GEN v;
244
245 if (l == 1) return cgetg(1,t_MAT);
246 if (l == 2) return gcopy(gel(A,1));
247 m = lgcols(A);
248 v = cgetg(m, t_COL);
249 for (i = 1; i < m; i++)
250 {
251 pari_sp av = avma;
252 GEN s = gcoeff(A,i,1);
253 for (j = 2; j < l; j++) s = gadd(s, gcoeff(A,i,j));
254 gel(v, i) = gerepileupto(av, s);
255 }
256 return v;
257 }
258
259 static GEN
_gmul(void * data,GEN x,GEN y)260 _gmul(void *data, GEN x, GEN y)
261 { (void)data; return gmul(x,y); }
262
263 GEN
RgV_prod(GEN x)264 RgV_prod(GEN x)
265 {
266 return gen_product(x, NULL, _gmul);
267 }
268
269 /* ADDITION SCALAR + MATRIX */
270 /* x square matrix, y scalar; create the square matrix x + y*Id */
271 GEN
RgM_Rg_add(GEN x,GEN y)272 RgM_Rg_add(GEN x, GEN y)
273 {
274 long l = lg(x), i, j;
275 GEN z = cgetg(l,t_MAT);
276
277 if (l==1) return z;
278 if (l != lgcols(x)) pari_err_OP( "+", x, y);
279 z = cgetg(l,t_MAT);
280 for (i=1; i<l; i++)
281 {
282 GEN zi = cgetg(l,t_COL), xi = gel(x,i);
283 gel(z,i) = zi;
284 for (j=1; j<l; j++)
285 gel(zi,j) = i==j? gadd(y,gel(xi,j)): gcopy(gel(xi,j));
286 }
287 return z;
288 }
289 GEN
RgM_Rg_sub(GEN x,GEN y)290 RgM_Rg_sub(GEN x, GEN y)
291 {
292 long l = lg(x), i, j;
293 GEN z = cgetg(l,t_MAT);
294
295 if (l==1) return z;
296 if (l != lgcols(x)) pari_err_OP( "-", x, y);
297 z = cgetg(l,t_MAT);
298 for (i=1; i<l; i++)
299 {
300 GEN zi = cgetg(l,t_COL), xi = gel(x,i);
301 gel(z,i) = zi;
302 for (j=1; j<l; j++)
303 gel(zi,j) = i==j? gsub(gel(xi,j), y): gcopy(gel(xi,j));
304 }
305 return z;
306 }
307 GEN
RgM_Rg_add_shallow(GEN x,GEN y)308 RgM_Rg_add_shallow(GEN x, GEN y)
309 {
310 long l = lg(x), i, j;
311 GEN z = cgetg(l,t_MAT);
312
313 if (l==1) return z;
314 if (l != lgcols(x)) pari_err_OP( "+", x, y);
315 for (i=1; i<l; i++)
316 {
317 GEN zi = cgetg(l,t_COL), xi = gel(x,i);
318 gel(z,i) = zi;
319 for (j=1; j<l; j++) gel(zi,j) = gel(xi,j);
320 gel(zi,i) = gadd(gel(zi,i), y);
321 }
322 return z;
323 }
324 GEN
RgM_Rg_sub_shallow(GEN x,GEN y)325 RgM_Rg_sub_shallow(GEN x, GEN y)
326 {
327 long l = lg(x), i, j;
328 GEN z = cgetg(l,t_MAT);
329
330 if (l==1) return z;
331 if (l != lgcols(x)) pari_err_OP( "-", x, y);
332 for (i=1; i<l; i++)
333 {
334 GEN zi = cgetg(l,t_COL), xi = gel(x,i);
335 gel(z,i) = zi;
336 for (j=1; j<l; j++) gel(zi,j) = gel(xi,j);
337 gel(zi,i) = gsub(gel(zi,i), y);
338 }
339 return z;
340 }
341
342 GEN
RgC_Rg_add(GEN x,GEN y)343 RgC_Rg_add(GEN x, GEN y)
344 {
345 long k, lx = lg(x);
346 GEN z = cgetg(lx, t_COL);
347 if (lx == 1)
348 {
349 if (isintzero(y)) return z;
350 pari_err_TYPE2("+",x,y);
351 }
352 gel(z,1) = gadd(y,gel(x,1));
353 for (k = 2; k < lx; k++) gel(z,k) = gcopy(gel(x,k));
354 return z;
355 }
356 GEN
RgC_Rg_sub(GEN x,GEN y)357 RgC_Rg_sub(GEN x, GEN y)
358 {
359 long k, lx = lg(x);
360 GEN z = cgetg(lx, t_COL);
361 if (lx == 1)
362 {
363 if (isintzero(y)) return z;
364 pari_err_TYPE2("-",x,y);
365 }
366 gel(z,1) = gsub(gel(x,1), y);
367 for (k = 2; k < lx; k++) gel(z,k) = gcopy(gel(x,k));
368 return z;
369 }
370 /* a - x */
371 GEN
Rg_RgC_sub(GEN a,GEN x)372 Rg_RgC_sub(GEN a, GEN x)
373 {
374 long k, lx = lg(x);
375 GEN z = cgetg(lx,t_COL);
376 if (lx == 1)
377 {
378 if (isintzero(a)) return z;
379 pari_err_TYPE2("-",a,x);
380 }
381 gel(z,1) = gsub(a, gel(x,1));
382 for (k = 2; k < lx; k++) gel(z,k) = gneg(gel(x,k));
383 return z;
384 }
385
386 static GEN
RgC_add_i(GEN x,GEN y,long lx)387 RgC_add_i(GEN x, GEN y, long lx)
388 {
389 GEN A = cgetg(lx, t_COL);
390 long i;
391 for (i=1; i<lx; i++) gel(A,i) = gadd(gel(x,i), gel(y,i));
392 return A;
393 }
394 GEN
RgC_add(GEN x,GEN y)395 RgC_add(GEN x, GEN y) { return RgC_add_i(x, y, lg(x)); }
396 GEN
RgV_add(GEN x,GEN y)397 RgV_add(GEN x, GEN y)
398 { pari_APPLY_type(t_VEC, gadd(gel(x,i), gel(y,i))) }
399
400 static GEN
RgC_sub_i(GEN x,GEN y,long lx)401 RgC_sub_i(GEN x, GEN y, long lx)
402 {
403 long i;
404 GEN A = cgetg(lx, t_COL);
405 for (i=1; i<lx; i++) gel(A,i) = gsub(gel(x,i), gel(y,i));
406 return A;
407 }
408 GEN
RgC_sub(GEN x,GEN y)409 RgC_sub(GEN x, GEN y) { return RgC_sub_i(x, y, lg(x)); }
410 GEN
RgV_sub(GEN x,GEN y)411 RgV_sub(GEN x, GEN y)
412 { pari_APPLY_type(t_VEC, gsub(gel(x,i), gel(y,i))) }
413
414 GEN
RgM_add(GEN x,GEN y)415 RgM_add(GEN x, GEN y)
416 {
417 long lx = lg(x), l, j;
418 GEN z;
419 if (lx == 1) return cgetg(1, t_MAT);
420 z = cgetg(lx, t_MAT); l = lgcols(x);
421 for (j = 1; j < lx; j++) gel(z,j) = RgC_add_i(gel(x,j), gel(y,j), l);
422 return z;
423 }
424 GEN
RgM_sub(GEN x,GEN y)425 RgM_sub(GEN x, GEN y)
426 {
427 long lx = lg(x), l, j;
428 GEN z;
429 if (lx == 1) return cgetg(1, t_MAT);
430 z = cgetg(lx, t_MAT); l = lgcols(x);
431 for (j = 1; j < lx; j++) gel(z,j) = RgC_sub_i(gel(x,j), gel(y,j), l);
432 return z;
433 }
434
435 static GEN
RgC_neg_i(GEN x,long lx)436 RgC_neg_i(GEN x, long lx)
437 {
438 long i;
439 GEN y = cgetg(lx, t_COL);
440 for (i=1; i<lx; i++) gel(y,i) = gneg(gel(x,i));
441 return y;
442 }
443 GEN
RgC_neg(GEN x)444 RgC_neg(GEN x) { return RgC_neg_i(x, lg(x)); }
445 GEN
RgV_neg(GEN x)446 RgV_neg(GEN x)
447 { pari_APPLY_type(t_VEC, gneg(gel(x,i))) }
448 GEN
RgM_neg(GEN x)449 RgM_neg(GEN x)
450 {
451 long i, hx, lx = lg(x);
452 GEN y = cgetg(lx, t_MAT);
453 if (lx == 1) return y;
454 hx = lgcols(x);
455 for (i=1; i<lx; i++) gel(y,i) = RgC_neg_i(gel(x,i), hx);
456 return y;
457 }
458
459 GEN
RgV_RgC_mul(GEN x,GEN y)460 RgV_RgC_mul(GEN x, GEN y)
461 {
462 long lx = lg(x);
463 if (lx != lg(y)) pari_err_OP("operation 'RgV_RgC_mul'", x, y);
464 return RgV_dotproduct_i(x, y, lx);
465 }
466 GEN
RgC_RgV_mul(GEN x,GEN y)467 RgC_RgV_mul(GEN x, GEN y)
468 {
469 long i, ly = lg(y);
470 GEN z = cgetg(ly,t_MAT);
471 for (i=1; i<ly; i++) gel(z,i) = RgC_Rg_mul(x, gel(y,i));
472 return z;
473 }
474 GEN
RgC_RgM_mul(GEN x,GEN y)475 RgC_RgM_mul(GEN x, GEN y)
476 {
477 long i, ly = lg(y);
478 GEN z = cgetg(ly,t_MAT);
479 if (ly != 1 && lgcols(y) != 2) pari_err_OP("operation 'RgC_RgM_mul'",x,y);
480 for (i=1; i<ly; i++) gel(z,i) = RgC_Rg_mul(x, gcoeff(y,1,i));
481 return z;
482 }
483 GEN
RgM_RgV_mul(GEN x,GEN y)484 RgM_RgV_mul(GEN x, GEN y)
485 {
486 if (lg(x) != 2) pari_err_OP("operation 'RgM_RgV_mul'", x,y);
487 return RgC_RgV_mul(gel(x,1), y);
488 }
489
490 /* x[i,]*y, l = lg(y) > 1 */
491 static GEN
RgMrow_RgC_mul_i(GEN x,GEN y,long i,long l)492 RgMrow_RgC_mul_i(GEN x, GEN y, long i, long l)
493 {
494 pari_sp av = avma;
495 GEN t = gmul(gcoeff(x,i,1), gel(y,1)); /* l > 1 ! */
496 long j;
497 for (j=2; j<l; j++)
498 {
499 GEN c = gcoeff(x,i,j);
500 if (!isintzero(c)) t = gadd(t, gmul(c, gel(y,j)));
501 }
502 return gerepileupto(av,t);
503 }
504 GEN
RgMrow_RgC_mul(GEN x,GEN y,long i)505 RgMrow_RgC_mul(GEN x, GEN y, long i)
506 { return RgMrow_RgC_mul_i(x, y, i, lg(x)); }
507
508 static GEN
RgM_RgC_mul_FpM(GEN x,GEN y,GEN p)509 RgM_RgC_mul_FpM(GEN x, GEN y, GEN p)
510 {
511 pari_sp av = avma;
512 GEN r;
513 if (lgefint(p) == 3)
514 {
515 ulong pp = uel(p, 2);
516 r = Flc_to_ZC_inplace(Flm_Flc_mul(RgM_to_Flm(x, pp),
517 RgV_to_Flv(y, pp), pp));
518 }
519 else
520 r = FpM_FpC_mul(RgM_to_FpM(x, p), RgC_to_FpC(y, p), p);
521 return gerepileupto(av, FpC_to_mod(r, p));
522 }
523
524 static GEN
RgM_RgC_mul_FqM(GEN x,GEN y,GEN pol,GEN p)525 RgM_RgC_mul_FqM(GEN x, GEN y, GEN pol, GEN p)
526 {
527 pari_sp av = avma;
528 GEN b, T = RgX_to_FpX(pol, p);
529 if (signe(T) == 0) pari_err_OP("*", x, y);
530 b = FqM_FqC_mul(RgM_to_FqM(x, T, p), RgC_to_FqC(y, T, p), T, p);
531 return gerepileupto(av, FqC_to_mod(b, T, p));
532 }
533
534 #define code(t1,t2) ((t1 << 6) | t2)
535 static GEN
RgM_RgC_mul_fast(GEN x,GEN y)536 RgM_RgC_mul_fast(GEN x, GEN y)
537 {
538 GEN p, pol;
539 long pa;
540 long t = RgM_RgC_type(x,y, &p,&pol,&pa);
541 switch(t)
542 {
543 case t_INT: return ZM_ZC_mul(x,y);
544 case t_FRAC: return QM_QC_mul(x,y);
545 case t_FFELT: return FFM_FFC_mul(x, y, pol);
546 case t_INTMOD: return RgM_RgC_mul_FpM(x, y, p);
547 case code(t_POLMOD, t_INTMOD):
548 return RgM_RgC_mul_FqM(x, y, pol, p);
549 default: return NULL;
550 }
551 }
552 #undef code
553
554 /* compatible t_MAT * t_COL, lx = lg(x) = lg(y) > 1, l = lgcols(x) */
555 static GEN
RgM_RgC_mul_i(GEN x,GEN y,long lx,long l)556 RgM_RgC_mul_i(GEN x, GEN y, long lx, long l)
557 {
558 GEN z = cgetg(l,t_COL);
559 long i;
560 for (i=1; i<l; i++) gel(z,i) = RgMrow_RgC_mul_i(x,y,i,lx);
561 return z;
562 }
563
564 GEN
RgM_RgC_mul(GEN x,GEN y)565 RgM_RgC_mul(GEN x, GEN y)
566 {
567 long lx = lg(x);
568 GEN z;
569 if (lx != lg(y)) pari_err_OP("operation 'RgM_RgC_mul'", x,y);
570 if (lx == 1) return cgetg(1,t_COL);
571 z = RgM_RgC_mul_fast(x, y);
572 if (z) return z;
573 return RgM_RgC_mul_i(x, y, lx, lgcols(x));
574 }
575
576 GEN
RgV_RgM_mul(GEN x,GEN y)577 RgV_RgM_mul(GEN x, GEN y)
578 {
579 long i, lx, ly = lg(y);
580 GEN z;
581 if (ly == 1) return cgetg(1,t_VEC);
582 lx = lg(x);
583 if (lx != lgcols(y)) pari_err_OP("operation 'RgV_RgM_mul'", x,y);
584 z = cgetg(ly, t_VEC);
585 for (i=1; i<ly; i++) gel(z,i) = RgV_dotproduct_i(x, gel(y,i), lx);
586 return z;
587 }
588
589 static GEN
RgM_mul_FpM(GEN x,GEN y,GEN p)590 RgM_mul_FpM(GEN x, GEN y, GEN p)
591 {
592 pari_sp av = avma;
593 GEN r;
594 if (lgefint(p) == 3)
595 {
596 ulong pp = uel(p, 2);
597 r = Flm_to_ZM_inplace(Flm_mul(RgM_to_Flm(x, pp),
598 RgM_to_Flm(y, pp), pp));
599 }
600 else
601 r = FpM_mul(RgM_to_FpM(x, p), RgM_to_FpM(y, p), p);
602 return gerepileupto(av, FpM_to_mod(r, p));
603 }
604
605 static GEN
RgM_mul_FqM(GEN x,GEN y,GEN pol,GEN p)606 RgM_mul_FqM(GEN x, GEN y, GEN pol, GEN p)
607 {
608 pari_sp av = avma;
609 GEN b, T = RgX_to_FpX(pol, p);
610 if (signe(T) == 0) pari_err_OP("*", x, y);
611 b = FqM_mul(RgM_to_FqM(x, T, p), RgM_to_FqM(y, T, p), T, p);
612 return gerepileupto(av, FqM_to_mod(b, T, p));
613 }
614
615 static GEN
RgM_liftred(GEN x,GEN T)616 RgM_liftred(GEN x, GEN T)
617 { return RgXQM_red(liftpol_shallow(x), T); }
618
619 static GEN
RgM_mul_ZXQM(GEN x,GEN y,GEN T)620 RgM_mul_ZXQM(GEN x, GEN y, GEN T)
621 {
622 pari_sp av = avma;
623 GEN b = ZXQM_mul(RgM_liftred(x,T), RgM_liftred(y, T), T);
624 return gerepilecopy(av, QXQM_to_mod_shallow(b,T));
625 }
626
627 static GEN
RgM_sqr_ZXQM(GEN x,GEN T)628 RgM_sqr_ZXQM(GEN x, GEN T)
629 {
630 pari_sp av = avma;
631 GEN b = ZXQM_sqr(RgM_liftred(x, T), T);
632 return gerepilecopy(av, QXQM_to_mod_shallow(b,T));
633 }
634
635 static GEN
RgM_mul_QXQM(GEN x,GEN y,GEN T)636 RgM_mul_QXQM(GEN x, GEN y, GEN T)
637 {
638 pari_sp av = avma;
639 GEN b = QXQM_mul(RgM_liftred(x, T), RgM_liftred(y, T), T);
640 return gerepilecopy(av, QXQM_to_mod_shallow(b,T));
641 }
642
643 static GEN
RgM_sqr_QXQM(GEN x,GEN T)644 RgM_sqr_QXQM(GEN x, GEN T)
645 {
646 pari_sp av = avma;
647 GEN b = QXQM_sqr(RgM_liftred(x, T), T);
648 return gerepilecopy(av, QXQM_to_mod_shallow(b,T));
649 }
650
651 INLINE int
RgX_is_monic_ZX(GEN pol)652 RgX_is_monic_ZX(GEN pol)
653 { return RgX_is_ZX(pol) && ZX_is_monic(pol); }
654
655 #define code(t1,t2) ((t1 << 6) | t2)
656 static GEN
RgM_mul_fast(GEN x,GEN y)657 RgM_mul_fast(GEN x, GEN y)
658 {
659 GEN p, pol;
660 long pa;
661 long t = RgM_type2(x,y, &p,&pol,&pa);
662 switch(t)
663 {
664 case t_INT: return ZM_mul(x,y);
665 case t_FRAC: return QM_mul(x,y);
666 case t_FFELT: return FFM_mul(x, y, pol);
667 case t_INTMOD: return RgM_mul_FpM(x, y, p);
668 case code(t_POLMOD, t_INT):
669 return ZX_is_monic(pol)? RgM_mul_ZXQM(x, y, pol): NULL;
670 case code(t_POLMOD, t_FRAC):
671 return RgX_is_monic_ZX(pol)? RgM_mul_QXQM(x, y, pol): NULL;
672 case code(t_POLMOD, t_INTMOD):
673 return RgM_mul_FqM(x, y, pol, p);
674 default: return NULL;
675 }
676 }
677
678 static GEN
RgM_sqr_fast(GEN x)679 RgM_sqr_fast(GEN x)
680 {
681 GEN p, pol;
682 long pa;
683 long t = RgM_type(x, &p,&pol,&pa);
684 switch(t)
685 {
686 case t_INT: return ZM_sqr(x);
687 case t_FRAC: return QM_sqr(x);
688 case t_FFELT: return FFM_mul(x, x, pol);
689 case t_INTMOD: return RgM_mul_FpM(x, x, p);
690 case code(t_POLMOD, t_INT):
691 return ZX_is_monic(pol)? RgM_sqr_ZXQM(x, pol): NULL;
692 case code(t_POLMOD, t_FRAC):
693 return RgX_is_monic_ZX(pol)? RgM_sqr_QXQM(x, pol): NULL;
694 case code(t_POLMOD, t_INTMOD):
695 return RgM_mul_FqM(x, x, pol, p);
696 default: return NULL;
697 }
698 }
699
700 #undef code
701
702 GEN
RgM_mul(GEN x,GEN y)703 RgM_mul(GEN x, GEN y)
704 {
705 long j, l, lx, ly = lg(y);
706 GEN z;
707 if (ly == 1) return cgetg(1,t_MAT);
708 lx = lg(x);
709 if (lx != lgcols(y)) pari_err_OP("operation 'RgM_mul'", x,y);
710 if (lx == 1) return zeromat(0,ly-1);
711 z = RgM_mul_fast(x, y);
712 if (z) return z;
713 z = cgetg(ly, t_MAT);
714 l = lgcols(x);
715 for (j=1; j<ly; j++) gel(z,j) = RgM_RgC_mul_i(x, gel(y,j), lx, l);
716 return z;
717 }
718
719 GEN
RgM_sqr(GEN x)720 RgM_sqr(GEN x)
721 {
722 long j, lx = lg(x);
723 GEN z;
724 if (lx == 1) return cgetg(1, t_MAT);
725 if (lx != lgcols(x)) pari_err_OP("operation 'RgM_mul'", x,x);
726 z = RgM_sqr_fast(x);
727 if (z) return z;
728 z = cgetg(lx, t_MAT);
729 for (j=1; j<lx; j++) gel(z,j) = RgM_RgC_mul_i(x, gel(x,j), lx, lx);
730 return z;
731 }
732
733 /* assume result is symmetric */
734 GEN
RgM_multosym(GEN x,GEN y)735 RgM_multosym(GEN x, GEN y)
736 {
737 long j, lx, ly = lg(y);
738 GEN M;
739 if (ly == 1) return cgetg(1,t_MAT);
740 lx = lg(x);
741 if (lx != lgcols(y)) pari_err_OP("operation 'RgM_multosym'", x,y);
742 if (lx == 1) return cgetg(1,t_MAT);
743 if (ly != lgcols(x)) pari_err_OP("operation 'RgM_multosym'", x,y);
744 M = cgetg(ly, t_MAT);
745 for (j=1; j<ly; j++)
746 {
747 GEN z = cgetg(ly,t_COL), yj = gel(y,j);
748 long i;
749 for (i=1; i<j; i++) gel(z,i) = gcoeff(M,j,i);
750 for (i=j; i<ly; i++)gel(z,i) = RgMrow_RgC_mul_i(x,yj,i,lx);
751 gel(M,j) = z;
752 }
753 return M;
754 }
755 /* x~ * y, assuming result is symmetric */
756 GEN
RgM_transmultosym(GEN x,GEN y)757 RgM_transmultosym(GEN x, GEN y)
758 {
759 long i, j, l, ly = lg(y);
760 GEN M;
761 if (ly == 1) return cgetg(1,t_MAT);
762 if (lg(x) != ly) pari_err_OP("operation 'RgM_transmultosym'", x,y);
763 l = lgcols(y);
764 if (lgcols(x) != l) pari_err_OP("operation 'RgM_transmultosym'", x,y);
765 M = cgetg(ly, t_MAT);
766 for (i=1; i<ly; i++)
767 {
768 GEN xi = gel(x,i), c = cgetg(ly,t_COL);
769 gel(M,i) = c;
770 for (j=1; j<i; j++)
771 gcoeff(M,i,j) = gel(c,j) = RgV_dotproduct_i(xi,gel(y,j),l);
772 gel(c,i) = RgV_dotproduct_i(xi,gel(y,i),l);
773 }
774 return M;
775 }
776 /* x~ * y */
777 GEN
RgM_transmul(GEN x,GEN y)778 RgM_transmul(GEN x, GEN y)
779 {
780 long i, j, l, lx, ly = lg(y);
781 GEN M;
782 if (ly == 1) return cgetg(1,t_MAT);
783 lx = lg(x);
784 l = lgcols(y);
785 if (lgcols(x) != l) pari_err_OP("operation 'RgM_transmul'", x,y);
786 M = cgetg(ly, t_MAT);
787 for (i=1; i<ly; i++)
788 {
789 GEN yi = gel(y,i), c = cgetg(lx,t_COL);
790 gel(M,i) = c;
791 for (j=1; j<lx; j++) gel(c,j) = RgV_dotproduct_i(yi,gel(x,j),l);
792 }
793 return M;
794 }
795
796 GEN
gram_matrix(GEN x)797 gram_matrix(GEN x)
798 {
799 long i,j, l, lx = lg(x);
800 GEN M;
801 if (!is_matvec_t(typ(x))) pari_err_TYPE("gram",x);
802 if (lx == 1) return cgetg(1,t_MAT);
803 l = lgcols(x);
804 M = cgetg(lx,t_MAT);
805 for (i=1; i<lx; i++)
806 {
807 GEN xi = gel(x,i), c = cgetg(lx,t_COL);
808 gel(M,i) = c;
809 for (j=1; j<i; j++)
810 gcoeff(M,i,j) = gel(c,j) = RgV_dotproduct_i(xi,gel(x,j),l);
811 gel(c,i) = RgV_dotsquare(xi);
812 }
813 return M;
814 }
815
816 static GEN
_RgM_add(void * E,GEN x,GEN y)817 _RgM_add(void *E, GEN x, GEN y) { (void)E; return RgM_add(x, y); }
818
819 static GEN
_RgM_sub(void * E,GEN x,GEN y)820 _RgM_sub(void *E, GEN x, GEN y) { (void)E; return RgM_sub(x, y); }
821
822 static GEN
_RgM_cmul(void * E,GEN P,long a,GEN x)823 _RgM_cmul(void *E, GEN P, long a, GEN x) { (void)E; return RgM_Rg_mul(x,gel(P,a+2)); }
824
825 static GEN
_RgM_sqr(void * E,GEN x)826 _RgM_sqr(void *E, GEN x) { (void) E; return RgM_sqr(x); }
827
828 static GEN
_RgM_mul(void * E,GEN x,GEN y)829 _RgM_mul(void *E, GEN x, GEN y) { (void) E; return RgM_mul(x, y); }
830
831 static GEN
_RgM_one(void * E)832 _RgM_one(void *E) { long *n = (long*) E; return matid(*n); }
833
834 static GEN
_RgM_zero(void * E)835 _RgM_zero(void *E) { long *n = (long*) E; return zeromat(*n,*n); }
836
837 static GEN
_RgM_red(void * E,GEN x)838 _RgM_red(void *E, GEN x) { (void)E; return x; }
839
840 static struct bb_algebra RgM_algebra = { _RgM_red, _RgM_add, _RgM_sub,
841 _RgM_mul, _RgM_sqr, _RgM_one, _RgM_zero };
842
843 /* generates the list of powers of x of degree 0,1,2,...,l*/
844 GEN
RgM_powers(GEN x,long l)845 RgM_powers(GEN x, long l)
846 {
847 long n = lg(x)-1;
848 return gen_powers(x,l,1,(void *) &n, &_RgM_sqr, &_RgM_mul, &_RgM_one);
849 }
850
851 GEN
RgX_RgMV_eval(GEN Q,GEN x)852 RgX_RgMV_eval(GEN Q, GEN x)
853 {
854 long n = lg(x)>1 ? lg(gel(x,1))-1:0;
855 return gen_bkeval_powers(Q,degpol(Q),x,(void*)&n,&RgM_algebra,&_RgM_cmul);
856 }
857
858 GEN
RgX_RgM_eval(GEN Q,GEN x)859 RgX_RgM_eval(GEN Q, GEN x)
860 {
861 long n = lg(x)-1;
862 return gen_bkeval(Q,degpol(Q),x,1,(void*)&n,&RgM_algebra,&_RgM_cmul);
863 }
864
865 GEN
RgC_Rg_div(GEN x,GEN y)866 RgC_Rg_div(GEN x, GEN y)
867 { pari_APPLY_type(t_COL, gdiv(gel(x,i),y)) }
868
869 GEN
RgC_Rg_mul(GEN x,GEN y)870 RgC_Rg_mul(GEN x, GEN y)
871 { pari_APPLY_type(t_COL, gmul(gel(x,i),y)) }
872
873 GEN
RgV_Rg_mul(GEN x,GEN y)874 RgV_Rg_mul(GEN x, GEN y)
875 { pari_APPLY_type(t_VEC, gmul(gel(x,i),y)) }
876
877 GEN
RgM_Rg_div(GEN X,GEN c)878 RgM_Rg_div(GEN X, GEN c) {
879 long i, j, h, l = lg(X);
880 GEN A = cgetg(l, t_MAT);
881 if (l == 1) return A;
882 h = lgcols(X);
883 for (j=1; j<l; j++)
884 {
885 GEN a = cgetg(h, t_COL), x = gel(X, j);
886 for (i = 1; i < h; i++) gel(a,i) = gdiv(gel(x,i), c);
887 gel(A,j) = a;
888 }
889 return A;
890 }
891 GEN
RgM_Rg_mul(GEN X,GEN c)892 RgM_Rg_mul(GEN X, GEN c) {
893 long i, j, h, l = lg(X);
894 GEN A = cgetg(l, t_MAT);
895 if (l == 1) return A;
896 h = lgcols(X);
897 for (j=1; j<l; j++)
898 {
899 GEN a = cgetg(h, t_COL), x = gel(X, j);
900 for (i = 1; i < h; i++) gel(a,i) = gmul(gel(x,i), c);
901 gel(A,j) = a;
902 }
903 return A;
904 }
905
906 /********************************************************************/
907 /* */
908 /* SCALAR TO MATRIX/VECTOR */
909 /* */
910 /********************************************************************/
911 /* fill the square nxn matrix equal to t*Id */
912 static void
fill_scalmat(GEN y,GEN t,long n)913 fill_scalmat(GEN y, GEN t, long n)
914 {
915 long i;
916 for (i = 1; i <= n; i++)
917 {
918 gel(y,i) = zerocol(n);
919 gcoeff(y,i,i) = t;
920 }
921 }
922
923 GEN
scalarmat(GEN x,long n)924 scalarmat(GEN x, long n) {
925 GEN y = cgetg(n+1, t_MAT);
926 if (!n) return y;
927 fill_scalmat(y, gcopy(x), n); return y;
928 }
929 GEN
scalarmat_shallow(GEN x,long n)930 scalarmat_shallow(GEN x, long n) {
931 GEN y = cgetg(n+1, t_MAT);
932 fill_scalmat(y, x, n); return y;
933 }
934 GEN
scalarmat_s(long x,long n)935 scalarmat_s(long x, long n) {
936 GEN y = cgetg(n+1, t_MAT);
937 if (!n) return y;
938 fill_scalmat(y, stoi(x), n); return y;
939 }
940 GEN
matid(long n)941 matid(long n) {
942 GEN y;
943 if (n < 0) pari_err_DOMAIN("matid", "size", "<", gen_0, stoi(n));
944 y = cgetg(n+1, t_MAT);
945 fill_scalmat(y, gen_1, n); return y;
946 }
947
948 INLINE GEN
scalarcol_i(GEN x,long n,long c)949 scalarcol_i(GEN x, long n, long c)
950 {
951 long i;
952 GEN y = cgetg(n+1,t_COL);
953 if (!n) return y;
954 gel(y,1) = c? gcopy(x): x;
955 for (i=2; i<=n; i++) gel(y,i) = gen_0;
956 return y;
957 }
958
959 GEN
scalarcol(GEN x,long n)960 scalarcol(GEN x, long n) { return scalarcol_i(x,n,1); }
961
962 GEN
scalarcol_shallow(GEN x,long n)963 scalarcol_shallow(GEN x, long n) { return scalarcol_i(x,n,0); }
964
965 int
RgM_isscalar(GEN x,GEN s)966 RgM_isscalar(GEN x, GEN s)
967 {
968 long i, j, lx = lg(x);
969
970 if (lx == 1) return 1;
971 if (lx != lgcols(x)) return 0;
972 if (!s) s = gcoeff(x,1,1);
973
974 for (j=1; j<lx; j++)
975 {
976 GEN c = gel(x,j);
977 for (i=1; i<j; )
978 if (!gequal0(gel(c,i++))) return 0;
979 /* i = j */
980 if (!gequal(gel(c,i++),s)) return 0;
981 for ( ; i<lx; )
982 if (!gequal0(gel(c,i++))) return 0;
983 }
984 return 1;
985 }
986
987 int
RgM_isidentity(GEN x)988 RgM_isidentity(GEN x)
989 {
990 long i,j, lx = lg(x);
991
992 if (lx == 1) return 1;
993 if (lx != lgcols(x)) return 0;
994 for (j=1; j<lx; j++)
995 {
996 GEN c = gel(x,j);
997 for (i=1; i<j; )
998 if (!gequal0(gel(c,i++))) return 0;
999 /* i = j */
1000 if (!gequal1(gel(c,i++))) return 0;
1001 for ( ; i<lx; )
1002 if (!gequal0(gel(c,i++))) return 0;
1003 }
1004 return 1;
1005 }
1006
1007 long
RgC_is_ei(GEN x)1008 RgC_is_ei(GEN x)
1009 {
1010 long i, j = 0, l = lg(x);
1011 for (i = 1; i < l; i++)
1012 {
1013 GEN c = gel(x,i);
1014 if (gequal0(c)) continue;
1015 if (!gequal1(c) || j) return 0;
1016 j = i;
1017 }
1018 return j;
1019 }
1020
1021 int
RgM_isdiagonal(GEN x)1022 RgM_isdiagonal(GEN x)
1023 {
1024 long i,j, lx = lg(x);
1025 if (lx == 1) return 1;
1026 if (lx != lgcols(x)) return 0;
1027
1028 for (j=1; j<lx; j++)
1029 {
1030 GEN c = gel(x,j);
1031 for (i=1; i<j; i++)
1032 if (!gequal0(gel(c,i))) return 0;
1033 for (i++; i<lx; i++)
1034 if (!gequal0(gel(c,i))) return 0;
1035 }
1036 return 1;
1037 }
1038 int
isdiagonal(GEN x)1039 isdiagonal(GEN x) { return (typ(x)==t_MAT) && RgM_isdiagonal(x); }
1040
1041 GEN
RgM_det_triangular(GEN mat)1042 RgM_det_triangular(GEN mat)
1043 {
1044 long i,l = lg(mat);
1045 pari_sp av;
1046 GEN s;
1047
1048 if (l<3) return l<2? gen_1: gcopy(gcoeff(mat,1,1));
1049 av = avma; s = gcoeff(mat,1,1);
1050 for (i=2; i<l; i++) s = gmul(s,gcoeff(mat,i,i));
1051 return av==avma? gcopy(s): gerepileupto(av,s);
1052 }
1053
1054 GEN
RgV_kill0(GEN v)1055 RgV_kill0(GEN v)
1056 {
1057 long i, l;
1058 GEN w = cgetg_copy(v, &l);
1059 for (i = 1; i < l; i++)
1060 {
1061 GEN a = gel(v,i);
1062 gel(w,i) = gequal0(a) ? NULL: a;
1063 }
1064 return w;
1065 }
1066