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 :
15 : #include "pari.h"
16 : #include "paripriv.h"
17 :
18 : int
19 1244630 : RgM_is_ZM(GEN x)
20 : {
21 1244630 : long i, j, h, l = lg(x);
22 1244630 : if (l == 1) return 1;
23 1243965 : h = lgcols(x);
24 1243960 : if (h == 1) return 1;
25 4700343 : for (j = l-1; j > 0; j--)
26 22384320 : for (i = h-1; i > 0; i--)
27 18927776 : if (typ(gcoeff(x,i,j)) != t_INT) return 0;
28 1083458 : return 1;
29 : }
30 :
31 : int
32 203 : RgM_is_QM(GEN x)
33 : {
34 203 : long i, j, h, l = lg(x);
35 203 : if (l == 1) return 1;
36 189 : h = lgcols(x);
37 189 : if (h == 1) return 1;
38 1218 : for (j = l-1; j > 0; j--)
39 15057 : for (i = h-1; i > 0; i--)
40 14028 : if (!is_rational_t(typ(gcoeff(x,i,j)))) return 0;
41 175 : return 1;
42 : }
43 :
44 : int
45 49 : RgV_is_ZMV(GEN V)
46 : {
47 49 : long i, l = lg(V);
48 399 : for (i=1; i<l; i++)
49 350 : if (typ(gel(V,i))!=t_MAT || !RgM_is_ZM(gel(V,i)))
50 0 : return 0;
51 49 : 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
63 8293887 : RgMrow_zc_mul_i(GEN x, GEN y, long c, long i)
64 : {
65 8293887 : pari_sp av = avma;
66 8293887 : GEN s = NULL;
67 : long j;
68 829019422 : for (j=1; j<c; j++)
69 : {
70 820728975 : long t = y[j];
71 820728975 : if (!t) continue;
72 28456189 : if (!s) { s = gmulgs(gcoeff(x,i,j),t); continue; }
73 20375774 : switch(t)
74 : {
75 14797989 : case 1: s = gadd(s, gcoeff(x,i,j)); break;
76 239940 : case -1: s = gsub(s, gcoeff(x,i,j)); break;
77 5337845 : default: s = gadd(s, gmulgs(gcoeff(x,i,j), t)); break;
78 : }
79 : }
80 8290447 : return s? gerepileupto(av, s): gc_const(av, gen_0);
81 : }
82 : GEN
83 70028 : 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
86 239230 : RgM_zc_mul_i(GEN x, GEN y, long c, long l)
87 : {
88 239230 : GEN z = cgetg(l,t_COL);
89 : long i;
90 8463059 : for (i = 1; i < l; i++) gel(z,i) = RgMrow_zc_mul_i(x,y,c,i);
91 239236 : return z;
92 : }
93 : GEN
94 71981 : 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
97 56858 : RgM_zm_mul(GEN x, GEN y)
98 : {
99 56858 : long j, c, l = lg(x), ly = lg(y);
100 56858 : GEN z = cgetg(ly, t_MAT);
101 56858 : if (l == 1) return z;
102 56858 : c = lgcols(x);
103 224107 : for (j = 1; j < ly; j++) gel(z,j) = RgM_zc_mul_i(x, gel(y,j), l,c);
104 56858 : return z;
105 : }
106 :
107 : /* x[i,]*y, l = lg(y) > 1 */
108 : static GEN
109 32272472 : RgMrow_ZC_mul_i(GEN x, GEN y, long i, long l)
110 : {
111 32272472 : pari_sp av = avma;
112 32272472 : GEN t = gmul(gcoeff(x,i,1), gel(y,1)); /* l > 1 ! */
113 : long j;
114 5021538642 : for (j=2; j<l; j++)
115 4989835742 : if (signe(gel(y,j))) t = gadd(t, gmul(gcoeff(x,i,j), gel(y,j)));
116 31702900 : return gerepileupto(av,t);
117 : }
118 :
119 : /* compatible t_MAT * t_COL, lx = lg(x) = lg(y) > 1, l = lgcols(x) */
120 : static GEN
121 1927815 : RgM_ZC_mul_i(GEN x, GEN y, long lx, long l)
122 : {
123 1927815 : GEN z = cgetg(l,t_COL);
124 : long i;
125 34225966 : for (i=1; i<l; i++) gel(z,i) = RgMrow_ZC_mul_i(x,y,i,lx);
126 1924815 : return z;
127 : }
128 :
129 : GEN
130 1927847 : RgM_ZM_mul_worker(GEN y, GEN x)
131 1927847 : { return RgM_ZC_mul_i(x, y, lg(x), lgcols(x)); }
132 :
133 : /* mostly useful when y is sparse */
134 : GEN
135 354736 : RgM_ZM_mul(GEN x, GEN y)
136 : {
137 354736 : pari_sp av = avma;
138 : GEN worker;
139 354736 : if (lg(x) == 1) return cgetg(lg(y), t_MAT);
140 354736 : worker = snm_closure(is_entry("_RgM_ZM_mul_worker"),mkvec(x));
141 354739 : return gerepileupto(av, gen_parapply(worker, y));
142 : }
143 :
144 : static GEN
145 210269 : RgV_zc_mul_i(GEN x, GEN y, long l)
146 : {
147 : long i;
148 210269 : GEN z = gen_0;
149 210269 : pari_sp av = avma;
150 8217895 : for (i = 1; i < l; i++) z = gadd(z, gmulgs(gel(x,i), y[i]));
151 210269 : return gerepileupto(av, z);
152 : }
153 : GEN
154 83384 : RgV_zc_mul(GEN x, GEN y) { return RgV_zc_mul_i(x, y, lg(x)); }
155 :
156 : GEN
157 34268 : RgV_zm_mul(GEN x, GEN y)
158 : {
159 34268 : long j, l = lg(x), ly = lg(y);
160 34268 : GEN z = cgetg(ly, t_VEC);
161 161153 : for (j = 1; j < ly; j++) gel(z,j) = RgV_zc_mul_i(x, gel(y,j), l);
162 34268 : return z;
163 : }
164 :
165 : /* scalar product x.x */
166 : GEN
167 9444712 : RgV_dotsquare(GEN x)
168 : {
169 9444712 : long i, lx = lg(x);
170 9444712 : pari_sp av = avma;
171 : GEN z;
172 9444712 : if (lx == 1) return gen_0;
173 9444712 : z = gsqr(gel(x,1));
174 21159062 : for (i=2; i<lx; i++)
175 : {
176 11714602 : z = gadd(z, gsqr(gel(x,i)));
177 11714502 : if (gc_needed(av,3))
178 : {
179 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgV_dotsquare, i = %ld",i);
180 0 : z = gerepileupto(av, z);
181 : }
182 : }
183 9444460 : return gerepileupto(av,z);
184 : }
185 :
186 : /* scalar product x.y, lx = lg(x) = lg(y) */
187 : static GEN
188 13603334 : RgV_dotproduct_i(GEN x, GEN y, long lx)
189 : {
190 13603334 : pari_sp av = avma;
191 : long i;
192 : GEN z;
193 13603334 : if (lx == 1) return gen_0;
194 13602935 : z = gmul(gel(x,1),gel(y,1));
195 218972070 : for (i=2; i<lx; i++)
196 : {
197 205371333 : z = gadd(z, gmul(gel(x,i), gel(y,i)));
198 205369663 : if (gc_needed(av,3))
199 : {
200 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgV_dotproduct, i = %ld",i);
201 0 : z = gerepileupto(av, z);
202 : }
203 : }
204 13600737 : return gerepileupto(av,z);
205 : }
206 : GEN
207 1012094 : RgV_dotproduct(GEN x,GEN y)
208 : {
209 1012094 : if (x == y) return RgV_dotsquare(x);
210 1012094 : return RgV_dotproduct_i(x, y, lg(x));
211 : }
212 : /* v[1] + ... + v[lg(v)-1] */
213 : GEN
214 3032588 : RgV_sum(GEN v)
215 : {
216 : GEN p;
217 3032588 : long i, l = lg(v);
218 3032588 : if (l == 1) return gen_0;
219 7380829 : p = gel(v,1); for (i=2; i<l; i++) p = gadd(p, gel(v,i));
220 3032560 : return p;
221 : }
222 : /* v[1] + ... + v[n]. Assume lg(v) > n. */
223 : GEN
224 0 : RgV_sumpart(GEN v, long n)
225 : {
226 : GEN p;
227 : long i;
228 0 : if (!n) return gen_0;
229 0 : p = gel(v,1); for (i=2; i<=n; i++) p = gadd(p, gel(v,i));
230 0 : return p;
231 : }
232 : /* v[m] + ... + v[n]. Assume lg(v) > n, m > 0. */
233 : GEN
234 0 : RgV_sumpart2(GEN v, long m, long n)
235 : {
236 : GEN p;
237 : long i;
238 0 : if (n < m) return gen_0;
239 0 : p = gel(v,m); for (i=m+1; i<=n; i++) p = gadd(p, gel(v,i));
240 0 : return p;
241 : }
242 : GEN
243 400 : RgM_sumcol(GEN A)
244 : {
245 400 : long i,j,m,l = lg(A);
246 : GEN v;
247 :
248 400 : if (l == 1) return cgetg(1,t_MAT);
249 400 : if (l == 2) return gcopy(gel(A,1));
250 216 : m = lgcols(A);
251 216 : v = cgetg(m, t_COL);
252 718 : for (i = 1; i < m; i++)
253 : {
254 502 : pari_sp av = avma;
255 502 : GEN s = gcoeff(A,i,1);
256 1214 : for (j = 2; j < l; j++) s = gadd(s, gcoeff(A,i,j));
257 502 : gel(v, i) = gerepileupto(av, s);
258 : }
259 216 : return v;
260 : }
261 :
262 : static GEN
263 3335357 : _gmul(void *data, GEN x, GEN y)
264 3335357 : { (void)data; return gmul(x,y); }
265 :
266 : GEN
267 452293 : RgV_prod(GEN x)
268 : {
269 452293 : return gen_product(x, NULL, _gmul);
270 : }
271 :
272 : /* ADDITION SCALAR + MATRIX */
273 : /* x square matrix, y scalar; create the square matrix x + y*Id */
274 : GEN
275 7881 : RgM_Rg_add(GEN x, GEN y)
276 : {
277 7881 : long l = lg(x), i, j;
278 7881 : GEN z = cgetg(l,t_MAT);
279 :
280 7881 : if (l==1) return z;
281 7881 : if (l != lgcols(x)) pari_err_OP("+", x, y);
282 58146 : for (i=1; i<l; i++)
283 : {
284 50265 : GEN zi = cgetg(l,t_COL), xi = gel(x,i);
285 50265 : gel(z,i) = zi;
286 2058774 : for (j=1; j<l; j++)
287 2008509 : gel(zi,j) = i==j? gadd(y,gel(xi,j)): gcopy(gel(xi,j));
288 : }
289 7881 : return z;
290 : }
291 : GEN
292 0 : RgM_Rg_sub(GEN x, GEN y)
293 : {
294 0 : long l = lg(x), i, j;
295 0 : GEN z = cgetg(l,t_MAT);
296 :
297 0 : if (l==1) return z;
298 0 : if (l != lgcols(x)) pari_err_OP("-", x, y);
299 0 : for (i=1; i<l; i++)
300 : {
301 0 : GEN zi = cgetg(l,t_COL), xi = gel(x,i);
302 0 : gel(z,i) = zi;
303 0 : for (j=1; j<l; j++)
304 0 : gel(zi,j) = i==j? gsub(gel(xi,j), y): gcopy(gel(xi,j));
305 : }
306 0 : return z;
307 : }
308 : GEN
309 739 : RgM_Rg_add_shallow(GEN x, GEN y)
310 : {
311 739 : long l = lg(x), i, j;
312 739 : GEN z = cgetg(l,t_MAT);
313 :
314 739 : if (l==1) return z;
315 739 : if (l != lgcols(x)) pari_err_OP( "+", x, y);
316 2384 : for (i=1; i<l; i++)
317 : {
318 1645 : GEN zi = cgetg(l,t_COL), xi = gel(x,i);
319 1645 : gel(z,i) = zi;
320 9456 : for (j=1; j<l; j++) gel(zi,j) = gel(xi,j);
321 1645 : gel(zi,i) = gadd(gel(zi,i), y);
322 : }
323 739 : return z;
324 : }
325 : GEN
326 312637 : RgM_Rg_sub_shallow(GEN x, GEN y)
327 : {
328 312637 : long l = lg(x), i, j;
329 312637 : GEN z = cgetg(l,t_MAT);
330 :
331 312636 : if (l==1) return z;
332 312636 : if (l != lgcols(x)) pari_err_OP( "-", x, y);
333 2015408 : for (i=1; i<l; i++)
334 : {
335 1702791 : GEN zi = cgetg(l,t_COL), xi = gel(x,i);
336 1702796 : gel(z,i) = zi;
337 27145605 : for (j=1; j<l; j++) gel(zi,j) = gel(xi,j);
338 1702796 : gel(zi,i) = gsub(gel(zi,i), y);
339 : }
340 312617 : return z;
341 : }
342 :
343 : GEN
344 5184209 : RgC_Rg_add(GEN x, GEN y)
345 : {
346 5184209 : long k, lx = lg(x);
347 5184209 : GEN z = cgetg(lx, t_COL);
348 5184209 : if (lx == 1)
349 : {
350 0 : if (isintzero(y)) return z;
351 0 : pari_err_TYPE2("+",x,y);
352 : }
353 5184209 : gel(z,1) = gadd(y,gel(x,1));
354 13031842 : for (k = 2; k < lx; k++) gel(z,k) = gcopy(gel(x,k));
355 5184195 : return z;
356 : }
357 : GEN
358 42099 : RgC_Rg_sub(GEN x, GEN y)
359 : {
360 42099 : long k, lx = lg(x);
361 42099 : GEN z = cgetg(lx, t_COL);
362 42099 : if (lx == 1)
363 : {
364 0 : if (isintzero(y)) return z;
365 0 : pari_err_TYPE2("-",x,y);
366 : }
367 42099 : gel(z,1) = gsub(gel(x,1), y);
368 108533 : for (k = 2; k < lx; k++) gel(z,k) = gcopy(gel(x,k));
369 42099 : return z;
370 : }
371 : /* a - x */
372 : GEN
373 135484 : Rg_RgC_sub(GEN a, GEN x)
374 : {
375 135484 : long k, lx = lg(x);
376 135484 : GEN z = cgetg(lx,t_COL);
377 135484 : if (lx == 1)
378 : {
379 0 : if (isintzero(a)) return z;
380 0 : pari_err_TYPE2("-",a,x);
381 : }
382 135484 : gel(z,1) = gsub(a, gel(x,1));
383 306053 : for (k = 2; k < lx; k++) gel(z,k) = gneg(gel(x,k));
384 135484 : return z;
385 : }
386 :
387 : static GEN
388 19474682 : RgC_add_i(GEN x, GEN y, long lx)
389 : {
390 19474682 : GEN A = cgetg(lx, t_COL);
391 : long i;
392 144958265 : for (i=1; i<lx; i++) gel(A,i) = gadd(gel(x,i), gel(y,i));
393 19474354 : return A;
394 : }
395 : GEN
396 16208452 : RgC_add(GEN x, GEN y) { return RgC_add_i(x, y, lg(x)); }
397 : GEN
398 2546803 : RgV_add(GEN x, GEN y)
399 10435407 : { pari_APPLY_type(t_VEC, gadd(gel(x,i), gel(y,i))) }
400 :
401 : static GEN
402 5999201 : RgC_sub_i(GEN x, GEN y, long lx)
403 : {
404 : long i;
405 5999201 : GEN A = cgetg(lx, t_COL);
406 635063396 : for (i=1; i<lx; i++) gel(A,i) = gsub(gel(x,i), gel(y,i));
407 5997838 : return A;
408 : }
409 : GEN
410 5097415 : RgC_sub(GEN x, GEN y) { return RgC_sub_i(x, y, lg(x)); }
411 : GEN
412 356509 : RgV_sub(GEN x, GEN y)
413 3032118 : { pari_APPLY_type(t_VEC, gsub(gel(x,i), gel(y,i))) }
414 :
415 : GEN
416 690774 : RgM_add(GEN x, GEN y)
417 : {
418 690774 : long lx = lg(x), l, j;
419 : GEN z;
420 690774 : if (lx == 1) return cgetg(1, t_MAT);
421 690774 : z = cgetg(lx, t_MAT); l = lgcols(x);
422 3957002 : for (j = 1; j < lx; j++) gel(z,j) = RgC_add_i(gel(x,j), gel(y,j), l);
423 690774 : return z;
424 : }
425 : GEN
426 172495 : RgM_sub(GEN x, GEN y)
427 : {
428 172495 : long lx = lg(x), l, j;
429 : GEN z;
430 172495 : if (lx == 1) return cgetg(1, t_MAT);
431 172495 : z = cgetg(lx, t_MAT); l = lgcols(x);
432 1074329 : for (j = 1; j < lx; j++) gel(z,j) = RgC_sub_i(gel(x,j), gel(y,j), l);
433 172495 : return z;
434 : }
435 :
436 : static GEN
437 3974856 : RgC_neg_i(GEN x, long lx)
438 : {
439 : long i;
440 3974856 : GEN y = cgetg(lx, t_COL);
441 30154666 : for (i=1; i<lx; i++) gel(y,i) = gneg(gel(x,i));
442 3974863 : return y;
443 : }
444 : GEN
445 990583 : RgC_neg(GEN x) { return RgC_neg_i(x, lg(x)); }
446 : GEN
447 4515 : RgV_neg(GEN x)
448 60473 : { pari_APPLY_type(t_VEC, gneg(gel(x,i))) }
449 : GEN
450 541274 : RgM_neg(GEN x)
451 : {
452 541274 : long i, hx, lx = lg(x);
453 541274 : GEN y = cgetg(lx, t_MAT);
454 541274 : if (lx == 1) return y;
455 541267 : hx = lgcols(x);
456 3525504 : for (i=1; i<lx; i++) gel(y,i) = RgC_neg_i(gel(x,i), hx);
457 541267 : return y;
458 : }
459 :
460 : GEN
461 4723090 : RgV_RgC_mul(GEN x, GEN y)
462 : {
463 4723090 : long lx = lg(x);
464 4723090 : if (lx != lg(y)) pari_err_OP("operation 'RgV_RgC_mul'", x, y);
465 4723006 : return RgV_dotproduct_i(x, y, lx);
466 : }
467 : GEN
468 1687 : RgC_RgV_mul(GEN x, GEN y)
469 : {
470 1687 : long i, ly = lg(y);
471 1687 : GEN z = cgetg(ly,t_MAT);
472 5061 : for (i=1; i<ly; i++) gel(z,i) = RgC_Rg_mul(x, gel(y,i));
473 1687 : return z;
474 : }
475 : GEN
476 0 : RgC_RgM_mul(GEN x, GEN y)
477 : {
478 0 : long i, ly = lg(y);
479 0 : GEN z = cgetg(ly,t_MAT);
480 0 : if (ly != 1 && lgcols(y) != 2) pari_err_OP("operation 'RgC_RgM_mul'",x,y);
481 0 : for (i=1; i<ly; i++) gel(z,i) = RgC_Rg_mul(x, gcoeff(y,1,i));
482 0 : return z;
483 : }
484 : GEN
485 0 : RgM_RgV_mul(GEN x, GEN y)
486 : {
487 0 : if (lg(x) != 2) pari_err_OP("operation 'RgM_RgV_mul'", x,y);
488 0 : return RgC_RgV_mul(gel(x,1), y);
489 : }
490 :
491 : /* x[i,]*y, l = lg(y) > 1 */
492 : static GEN
493 115784674 : RgMrow_RgC_mul_i(GEN x, GEN y, long i, long l)
494 : {
495 115784674 : pari_sp av = avma;
496 115784674 : GEN t = gmul(gcoeff(x,i,1), gel(y,1)); /* l > 1 ! */
497 : long j;
498 985260808 : for (j=2; j<l; j++)
499 : {
500 869495613 : GEN c = gcoeff(x,i,j);
501 869495613 : if (!isintzero(c)) t = gadd(t, gmul(c, gel(y,j)));
502 : }
503 115765195 : return gerepileupto(av,t);
504 : }
505 : GEN
506 4655 : RgMrow_RgC_mul(GEN x, GEN y, long i)
507 4655 : { return RgMrow_RgC_mul_i(x, y, i, lg(x)); }
508 :
509 : static GEN
510 28 : RgM_RgC_mul_FpM(GEN x, GEN y, GEN p)
511 : {
512 28 : pari_sp av = avma;
513 : GEN r;
514 28 : if (lgefint(p) == 3)
515 : {
516 14 : ulong pp = uel(p, 2);
517 14 : r = Flm_Flc_mul(RgM_to_Flm(x, pp), RgV_to_Flv(y, pp), pp);
518 14 : r = Flc_to_ZC_inplace(r);
519 : }
520 : else
521 14 : r = FpM_FpC_mul(RgM_to_FpM(x, p), RgC_to_FpC(y, p), p);
522 28 : return gerepileupto(av, FpC_to_mod(r, p));
523 : }
524 :
525 : static GEN
526 28 : RgM_RgC_mul_FqM(GEN x, GEN y, GEN pol, GEN p)
527 : {
528 28 : pari_sp av = avma;
529 28 : GEN b, T = RgX_to_FpX(pol, p);
530 28 : if (signe(T) == 0) pari_err_OP("*", x, y);
531 28 : b = FqM_FqC_mul(RgM_to_FqM(x, T, p), RgC_to_FqC(y, T, p), T, p);
532 28 : return gerepileupto(av, FqC_to_mod(b, T, p));
533 : }
534 :
535 : #define code(t1,t2) ((t1 << 6) | t2)
536 : static GEN
537 31196188 : RgM_RgC_mul_fast(GEN x, GEN y)
538 : {
539 : GEN p, pol;
540 : long pa;
541 31196188 : long t = RgM_RgC_type(x,y, &p,&pol,&pa);
542 31196469 : switch(t)
543 : {
544 8838829 : case t_INT: return ZM_ZC_mul(x,y);
545 142985 : case t_FRAC: return QM_QC_mul(x,y);
546 91 : case t_FFELT: return FFM_FFC_mul(x, y, pol);
547 28 : case t_INTMOD: return RgM_RgC_mul_FpM(x, y, p);
548 27 : case code(t_POLMOD, t_INTMOD):
549 27 : return RgM_RgC_mul_FqM(x, y, pol, p);
550 22214509 : default: return NULL;
551 : }
552 : }
553 : #undef code
554 :
555 : /* compatible t_MAT * t_COL, lx = lg(x) = lg(y) > 1, l = lgcols(x) */
556 : static GEN
557 33998460 : RgM_RgC_mul_i(GEN x, GEN y, long lx, long l)
558 : {
559 33998460 : GEN z = cgetg(l,t_COL);
560 : long i;
561 149772550 : for (i=1; i<l; i++) gel(z,i) = RgMrow_RgC_mul_i(x,y,i,lx);
562 33993174 : return z;
563 : }
564 :
565 : GEN
566 31196183 : RgM_RgC_mul(GEN x, GEN y)
567 : {
568 31196183 : long lx = lg(x);
569 : GEN z;
570 31196183 : if (lx != lg(y)) pari_err_OP("operation 'RgM_RgC_mul'", x,y);
571 31196188 : if (lx == 1) return cgetg(1,t_COL);
572 31196188 : z = RgM_RgC_mul_fast(x, y);
573 31196456 : if (z) return z;
574 22214497 : return RgM_RgC_mul_i(x, y, lx, lgcols(x));
575 : }
576 :
577 : GEN
578 295533 : RgV_RgM_mul(GEN x, GEN y)
579 : {
580 295533 : long i, lx, ly = lg(y);
581 : GEN z;
582 295533 : if (ly == 1) return cgetg(1,t_VEC);
583 295526 : lx = lg(x);
584 295526 : if (lx != lgcols(y)) pari_err_OP("operation 'RgV_RgM_mul'", x,y);
585 295518 : z = cgetg(ly, t_VEC);
586 3423052 : for (i=1; i<ly; i++) gel(z,i) = RgV_dotproduct_i(x, gel(y,i), lx);
587 295367 : return z;
588 : }
589 :
590 : static GEN
591 57148 : RgM_mul_FpM_i(GEN x, GEN y, GEN p)
592 : {
593 57148 : if (lgefint(p) == 3)
594 : {
595 57076 : ulong pp = uel(p, 2);
596 57076 : if (pp==2) return F2m_to_mod(F2m_mul(RgM_to_F2m(x), RgM_to_F2m(y)));
597 57027 : if (pp==3) return F3m_to_mod(F3m_mul(RgM_to_F3m(x), RgM_to_F3m(y)));
598 56873 : return Flm_to_mod(Flm_mul(RgM_to_Flm(x, pp), RgM_to_Flm(y, pp), pp), pp);
599 : }
600 72 : return FpM_to_mod(FpM_mul(RgM_to_FpM(x,p), RgM_to_FpM(y,p), p), p);
601 : }
602 : static GEN
603 57148 : RgM_mul_FpM(GEN x, GEN y, GEN p)
604 57148 : { pari_sp av = avma; return gerepileupto(av, RgM_mul_FpM_i(x, y, p)); }
605 :
606 : static GEN
607 26005 : RgM_mul_FqM(GEN x, GEN y, GEN pol, GEN p)
608 : {
609 26005 : pari_sp av = avma;
610 26005 : GEN b, T = RgX_to_FpX(pol, p);
611 26005 : if (signe(T) == 0) pari_err_OP("*", x, y);
612 26005 : b = FqM_mul(RgM_to_FqM(x, T, p), RgM_to_FqM(y, T, p), T, p);
613 26005 : return gerepileupto(av, FqM_to_mod(b, T, p));
614 : }
615 :
616 : static GEN
617 12418 : RgM_liftred(GEN x, GEN T)
618 12418 : { return RgXQM_red(liftpol_shallow(x), T); }
619 :
620 : static GEN
621 1449 : RgM_mul_ZXQM(GEN x, GEN y, GEN T)
622 : {
623 1449 : pari_sp av = avma;
624 1449 : GEN b = ZXQM_mul(RgM_liftred(x,T), RgM_liftred(y, T), T);
625 1449 : return gerepilecopy(av, QXQM_to_mod_shallow(b,T));
626 : }
627 :
628 : static GEN
629 133 : RgM_sqr_ZXQM(GEN x, GEN T)
630 : {
631 133 : pari_sp av = avma;
632 133 : GEN b = ZXQM_sqr(RgM_liftred(x, T), T);
633 133 : return gerepilecopy(av, QXQM_to_mod_shallow(b,T));
634 : }
635 :
636 : static GEN
637 4690 : RgM_mul_QXQM(GEN x, GEN y, GEN T)
638 : {
639 4690 : pari_sp av = avma;
640 4690 : GEN b = QXQM_mul(RgM_liftred(x, T), RgM_liftred(y, T), T);
641 4690 : return gerepilecopy(av, QXQM_to_mod_shallow(b,T));
642 : }
643 :
644 : static GEN
645 7 : RgM_sqr_QXQM(GEN x, GEN T)
646 : {
647 7 : pari_sp av = avma;
648 7 : GEN b = QXQM_sqr(RgM_liftred(x, T), T);
649 7 : return gerepilecopy(av, QXQM_to_mod_shallow(b,T));
650 : }
651 :
652 : INLINE int
653 4739 : RgX_is_monic_ZX(GEN pol)
654 4739 : { return RgX_is_ZX(pol) && ZX_is_monic(pol); }
655 :
656 : #define code(t1,t2) ((t1 << 6) | t2)
657 : static GEN
658 7570913 : RgM_mul_fast(GEN x, GEN y)
659 : {
660 : GEN p, pol;
661 : long pa;
662 7570913 : long t = RgM_type2(x,y, &p,&pol,&pa);
663 7570968 : switch(t)
664 : {
665 2689556 : case t_INT: return ZM_mul(x,y);
666 117138 : case t_FRAC: return QM_mul(x,y);
667 4403 : case t_FFELT: return FFM_mul(x, y, pol);
668 57085 : case t_INTMOD: return RgM_mul_FpM(x, y, p);
669 1456 : case code(t_POLMOD, t_INT):
670 1456 : return ZX_is_monic(pol)? RgM_mul_ZXQM(x, y, pol): NULL;
671 4711 : case code(t_POLMOD, t_FRAC):
672 4711 : return RgX_is_monic_ZX(pol)? RgM_mul_QXQM(x, y, pol): NULL;
673 26002 : case code(t_POLMOD, t_INTMOD):
674 26002 : return RgM_mul_FqM(x, y, pol, p);
675 4670617 : default: return NULL;
676 : }
677 : }
678 :
679 : static GEN
680 1316 : RgM_sqr_fast(GEN x)
681 : {
682 : GEN p, pol;
683 : long pa;
684 1316 : long t = RgM_type(x, &p,&pol,&pa);
685 1316 : switch(t)
686 : {
687 175 : case t_INT: return ZM_sqr(x);
688 700 : case t_FRAC: return QM_sqr(x);
689 112 : case t_FFELT: return FFM_mul(x, x, pol);
690 63 : case t_INTMOD: return RgM_mul_FpM(x, x, p);
691 140 : case code(t_POLMOD, t_INT):
692 140 : return ZX_is_monic(pol)? RgM_sqr_ZXQM(x, pol): NULL;
693 28 : case code(t_POLMOD, t_FRAC):
694 28 : return RgX_is_monic_ZX(pol)? RgM_sqr_QXQM(x, pol): NULL;
695 0 : case code(t_POLMOD, t_INTMOD):
696 0 : return RgM_mul_FqM(x, x, pol, p);
697 98 : default: return NULL;
698 : }
699 : }
700 :
701 : #undef code
702 :
703 : /* lx, ly > 1 */
704 : static GEN
705 4670645 : RgM_mul_i(GEN x, GEN y, long lx, long ly)
706 : {
707 4670645 : GEN z = cgetg(ly, t_MAT);
708 4670642 : long j, l = lgcols(x);
709 16454443 : for (j = 1; j < ly; j++) gel(z,j) = RgM_RgC_mul_i(x, gel(y,j), lx, l);
710 4670599 : return z;
711 : }
712 : GEN
713 7586695 : RgM_mul(GEN x, GEN y)
714 : {
715 7586695 : long lx, ly = lg(y);
716 : GEN z;
717 7586695 : if (ly == 1) return cgetg(1,t_MAT);
718 7570945 : lx = lg(x);
719 7570945 : if (lx != lgcols(y)) pari_err_OP("operation 'RgM_mul'", x,y);
720 7570941 : if (lx == 1) return zeromat(0,ly-1);
721 7570913 : z = RgM_mul_fast(x, y);
722 7570966 : if (z) return z;
723 4670640 : return RgM_mul_i(x, y, lx, ly);
724 : }
725 :
726 : GEN
727 1351 : RgM_sqr(GEN x)
728 : {
729 1351 : long j, lx = lg(x);
730 : GEN z;
731 1351 : if (lx == 1) return cgetg(1, t_MAT);
732 1316 : if (lx != lgcols(x)) pari_err_OP("operation 'RgM_mul'", x,x);
733 1316 : z = RgM_sqr_fast(x);
734 1316 : if (z) return z;
735 126 : z = cgetg(lx, t_MAT);
736 497 : for (j=1; j<lx; j++) gel(z,j) = RgM_RgC_mul_i(x, gel(x,j), lx, lx);
737 126 : return z;
738 : }
739 :
740 : /* assume result is symmetric */
741 : GEN
742 0 : RgM_multosym(GEN x, GEN y)
743 : {
744 0 : long j, lx, ly = lg(y);
745 : GEN M;
746 0 : if (ly == 1) return cgetg(1,t_MAT);
747 0 : lx = lg(x);
748 0 : if (lx != lgcols(y)) pari_err_OP("operation 'RgM_multosym'", x,y);
749 0 : if (lx == 1) return cgetg(1,t_MAT);
750 0 : if (ly != lgcols(x)) pari_err_OP("operation 'RgM_multosym'", x,y);
751 0 : M = cgetg(ly, t_MAT);
752 0 : for (j=1; j<ly; j++)
753 : {
754 0 : GEN z = cgetg(ly,t_COL), yj = gel(y,j);
755 : long i;
756 0 : for (i=1; i<j; i++) gel(z,i) = gcoeff(M,j,i);
757 0 : for (i=j; i<ly; i++)gel(z,i) = RgMrow_RgC_mul_i(x,yj,i,lx);
758 0 : gel(M,j) = z;
759 : }
760 0 : return M;
761 : }
762 : /* x~ * y, assuming result is symmetric */
763 : GEN
764 7663 : RgM_transmultosym(GEN x, GEN y)
765 : {
766 7663 : long i, j, l, ly = lg(y);
767 : GEN M;
768 7663 : if (ly == 1) return cgetg(1,t_MAT);
769 7663 : if (lg(x) != ly) pari_err_OP("operation 'RgM_transmultosym'", x,y);
770 7663 : l = lgcols(y);
771 7663 : if (lgcols(x) != l) pari_err_OP("operation 'RgM_transmultosym'", x,y);
772 7663 : M = cgetg(ly, t_MAT);
773 30820 : for (i=1; i<ly; i++)
774 : {
775 23157 : GEN xi = gel(x,i), c = cgetg(ly,t_COL);
776 23157 : gel(M,i) = c;
777 47490 : for (j=1; j<i; j++)
778 24333 : gcoeff(M,i,j) = gel(c,j) = RgV_dotproduct_i(xi,gel(y,j),l);
779 23157 : gel(c,i) = RgV_dotproduct_i(xi,gel(y,i),l);
780 : }
781 7663 : return M;
782 : }
783 : /* x~ * y */
784 : GEN
785 0 : RgM_transmul(GEN x, GEN y)
786 : {
787 0 : long i, j, l, lx, ly = lg(y);
788 : GEN M;
789 0 : if (ly == 1) return cgetg(1,t_MAT);
790 0 : lx = lg(x);
791 0 : l = lgcols(y);
792 0 : if (lgcols(x) != l) pari_err_OP("operation 'RgM_transmul'", x,y);
793 0 : M = cgetg(ly, t_MAT);
794 0 : for (i=1; i<ly; i++)
795 : {
796 0 : GEN yi = gel(y,i), c = cgetg(lx,t_COL);
797 0 : gel(M,i) = c;
798 0 : for (j=1; j<lx; j++) gel(c,j) = RgV_dotproduct_i(yi,gel(x,j),l);
799 : }
800 0 : return M;
801 : }
802 :
803 : GEN
804 4693117 : gram_matrix(GEN x)
805 : {
806 4693117 : long i,j, l, lx = lg(x);
807 : GEN M;
808 4693117 : if (!is_matvec_t(typ(x))) pari_err_TYPE("gram",x);
809 4693115 : if (lx == 1) return cgetg(1,t_MAT);
810 4693101 : l = lgcols(x);
811 4693101 : M = cgetg(lx,t_MAT);
812 14079237 : for (i=1; i<lx; i++)
813 : {
814 9386157 : GEN xi = gel(x,i), c = cgetg(lx,t_COL);
815 9386153 : gel(M,i) = c;
816 14079222 : for (j=1; j<i; j++)
817 4693074 : gcoeff(M,i,j) = gel(c,j) = RgV_dotproduct_i(xi,gel(x,j),l);
818 9386148 : gel(c,i) = RgV_dotsquare(xi);
819 : }
820 4693080 : return M;
821 : }
822 :
823 : static GEN
824 3801 : _RgM_add(void *E, GEN x, GEN y) { (void)E; return RgM_add(x, y); }
825 :
826 : static GEN
827 0 : _RgM_sub(void *E, GEN x, GEN y) { (void)E; return RgM_sub(x, y); }
828 :
829 : static GEN
830 5684 : _RgM_cmul(void *E, GEN P, long a, GEN x) { (void)E; return RgM_Rg_mul(x,gel(P,a+2)); }
831 :
832 : static GEN
833 238 : _RgM_sqr(void *E, GEN x) { (void) E; return RgM_sqr(x); }
834 :
835 : static GEN
836 749 : _RgM_mul(void *E, GEN x, GEN y) { (void) E; return RgM_mul(x, y); }
837 :
838 : static GEN
839 4046 : _RgM_one(void *E) { long *n = (long*) E; return matid(*n); }
840 :
841 : static GEN
842 0 : _RgM_zero(void *E) { long *n = (long*) E; return zeromat(*n,*n); }
843 :
844 : static GEN
845 2828 : _RgM_red(void *E, GEN x) { (void)E; return x; }
846 :
847 : static struct bb_algebra RgM_algebra = { _RgM_red, _RgM_add, _RgM_sub,
848 : _RgM_mul, _RgM_sqr, _RgM_one, _RgM_zero };
849 :
850 : /* generates the list of powers of x of degree 0,1,2,...,l*/
851 : GEN
852 168 : RgM_powers(GEN x, long l)
853 : {
854 168 : long n = lg(x)-1;
855 168 : return gen_powers(x,l,1,(void *) &n, &_RgM_sqr, &_RgM_mul, &_RgM_one);
856 : }
857 :
858 : GEN
859 490 : RgX_RgMV_eval(GEN Q, GEN x)
860 : {
861 490 : long n = lg(x)>1 ? lg(gel(x,1))-1:0;
862 490 : return gen_bkeval_powers(Q,degpol(Q),x,(void*)&n,&RgM_algebra,&_RgM_cmul);
863 : }
864 :
865 : GEN
866 1393 : RgX_RgM_eval(GEN Q, GEN x)
867 : {
868 1393 : long n = lg(x)-1;
869 1393 : return gen_bkeval(Q,degpol(Q),x,1,(void*)&n,&RgM_algebra,&_RgM_cmul);
870 : }
871 :
872 : GEN
873 4458105 : RgC_Rg_div(GEN x, GEN y)
874 19197933 : { pari_APPLY_type(t_COL, gdiv(gel(x,i),y)) }
875 :
876 : GEN
877 13201820 : RgC_Rg_mul(GEN x, GEN y)
878 79787901 : { pari_APPLY_type(t_COL, gmul(gel(x,i),y)) }
879 :
880 : GEN
881 31157 : RgV_Rg_mul(GEN x, GEN y)
882 1200766 : { pari_APPLY_type(t_VEC, gmul(gel(x,i),y)) }
883 :
884 : GEN
885 609428 : RgM_Rg_div(GEN X, GEN c) {
886 609428 : long i, j, h, l = lg(X);
887 609428 : GEN A = cgetg(l, t_MAT);
888 609429 : if (l == 1) return A;
889 609429 : h = lgcols(X);
890 2790634 : for (j=1; j<l; j++)
891 : {
892 2181259 : GEN a = cgetg(h, t_COL), x = gel(X, j);
893 17265416 : for (i = 1; i < h; i++) gel(a,i) = gdiv(gel(x,i), c);
894 2181205 : gel(A,j) = a;
895 : }
896 609375 : return A;
897 : }
898 : GEN
899 248329 : RgM_Rg_mul(GEN X, GEN c) {
900 248329 : long i, j, h, l = lg(X);
901 248329 : GEN A = cgetg(l, t_MAT);
902 248330 : if (l == 1) return A;
903 248330 : h = lgcols(X);
904 1090579 : for (j=1; j<l; j++)
905 : {
906 842267 : GEN a = cgetg(h, t_COL), x = gel(X, j);
907 7072112 : for (i = 1; i < h; i++) gel(a,i) = gmul(gel(x,i), c);
908 842249 : gel(A,j) = a;
909 : }
910 248312 : return A;
911 : }
912 :
913 : /********************************************************************/
914 : /* */
915 : /* SCALAR TO MATRIX/VECTOR */
916 : /* */
917 : /********************************************************************/
918 : /* fill the square nxn matrix equal to t*Id */
919 : static void
920 13409945 : fill_scalmat(GEN y, GEN t, long n)
921 : {
922 : long i;
923 51077107 : for (i = 1; i <= n; i++)
924 : {
925 37667522 : gel(y,i) = zerocol(n);
926 37667162 : gcoeff(y,i,i) = t;
927 : }
928 13409585 : }
929 :
930 : GEN
931 930163 : scalarmat(GEN x, long n) {
932 930163 : GEN y = cgetg(n+1, t_MAT);
933 930133 : if (!n) return y;
934 930133 : fill_scalmat(y, gcopy(x), n); return y;
935 : }
936 : GEN
937 3155162 : scalarmat_shallow(GEN x, long n) {
938 3155162 : GEN y = cgetg(n+1, t_MAT);
939 3155238 : fill_scalmat(y, x, n); return y;
940 : }
941 : GEN
942 217 : scalarmat_s(long x, long n) {
943 217 : GEN y = cgetg(n+1, t_MAT);
944 217 : if (!n) return y;
945 217 : fill_scalmat(y, stoi(x), n); return y;
946 : }
947 : GEN
948 9325140 : matid(long n) {
949 : GEN y;
950 9325140 : if (n < 0) pari_err_DOMAIN("matid", "size", "<", gen_0, stoi(n));
951 9325133 : y = cgetg(n+1, t_MAT);
952 9325144 : fill_scalmat(y, gen_1, n); return y;
953 : }
954 :
955 : INLINE GEN
956 2260068 : scalarcol_i(GEN x, long n, long c)
957 : {
958 : long i;
959 2260068 : GEN y = cgetg(n+1,t_COL);
960 2260065 : if (!n) return y;
961 2260065 : gel(y,1) = c? gcopy(x): x;
962 6742170 : for (i=2; i<=n; i++) gel(y,i) = gen_0;
963 2260067 : return y;
964 : }
965 :
966 : GEN
967 777170 : scalarcol(GEN x, long n) { return scalarcol_i(x,n,1); }
968 :
969 : GEN
970 1482899 : scalarcol_shallow(GEN x, long n) { return scalarcol_i(x,n,0); }
971 :
972 : int
973 33511 : RgM_isscalar(GEN x, GEN s)
974 : {
975 33511 : long i, j, lx = lg(x);
976 :
977 33511 : if (lx == 1) return 1;
978 33511 : if (lx != lgcols(x)) return 0;
979 33511 : if (!s) s = gcoeff(x,1,1);
980 :
981 92671 : for (j=1; j<lx; j++)
982 : {
983 75244 : GEN c = gel(x,j);
984 202502 : for (i=1; i<j; )
985 140710 : if (!gequal0(gel(c,i++))) return 0;
986 : /* i = j */
987 61792 : if (!gequal(gel(c,i++),s)) return 0;
988 225389 : for ( ; i<lx; )
989 166229 : if (!gequal0(gel(c,i++))) return 0;
990 : }
991 17427 : return 1;
992 : }
993 :
994 : int
995 15441 : RgM_isidentity(GEN x)
996 : {
997 15441 : long i,j, lx = lg(x);
998 :
999 15441 : if (lx == 1) return 1;
1000 15441 : if (lx != lgcols(x)) return 0;
1001 25473 : for (j=1; j<lx; j++)
1002 : {
1003 25242 : GEN c = gel(x,j);
1004 58750 : for (i=1; i<j; )
1005 38477 : if (!gequal0(gel(c,i++))) return 0;
1006 : /* i = j */
1007 20273 : if (!gequal1(gel(c,i++))) return 0;
1008 61441 : for ( ; i<lx; )
1009 51409 : if (!gequal0(gel(c,i++))) return 0;
1010 : }
1011 231 : return 1;
1012 : }
1013 :
1014 : long
1015 539 : RgC_is_ei(GEN x)
1016 : {
1017 539 : long i, j = 0, l = lg(x);
1018 2947 : for (i = 1; i < l; i++)
1019 : {
1020 2408 : GEN c = gel(x,i);
1021 2408 : if (gequal0(c)) continue;
1022 539 : if (!gequal1(c) || j) return 0;
1023 539 : j = i;
1024 : }
1025 539 : return j;
1026 : }
1027 :
1028 : int
1029 336 : RgM_isdiagonal(GEN x)
1030 : {
1031 336 : long i,j, lx = lg(x);
1032 336 : if (lx == 1) return 1;
1033 336 : if (lx != lgcols(x)) return 0;
1034 :
1035 3220 : for (j=1; j<lx; j++)
1036 : {
1037 2891 : GEN c = gel(x,j);
1038 18452 : for (i=1; i<j; i++)
1039 15561 : if (!gequal0(gel(c,i))) return 0;
1040 18452 : for (i++; i<lx; i++)
1041 15568 : if (!gequal0(gel(c,i))) return 0;
1042 : }
1043 329 : return 1;
1044 : }
1045 : int
1046 315 : isdiagonal(GEN x) { return (typ(x)==t_MAT) && RgM_isdiagonal(x); }
1047 :
1048 : GEN
1049 21610 : RgM_det_triangular(GEN mat)
1050 : {
1051 21610 : long i,l = lg(mat);
1052 : pari_sp av;
1053 : GEN s;
1054 :
1055 21610 : if (l<3) return l<2? gen_1: gcopy(gcoeff(mat,1,1));
1056 19979 : av = avma; s = gcoeff(mat,1,1);
1057 120356 : for (i=2; i<l; i++) s = gmul(s,gcoeff(mat,i,i));
1058 19979 : return av==avma? gcopy(s): gerepileupto(av,s);
1059 : }
1060 :
1061 : GEN
1062 478034 : RgV_kill0(GEN v)
1063 : {
1064 : long i, l;
1065 478034 : GEN w = cgetg_copy(v, &l);
1066 130025385 : for (i = 1; i < l; i++)
1067 : {
1068 129547438 : GEN a = gel(v,i);
1069 129547438 : gel(w,i) = gequal0(a) ? NULL: a;
1070 : }
1071 477947 : return w;
1072 : }
|
