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 1233258 : RgM_is_ZM(GEN x)
20 : {
21 1233258 : long i, j, h, l = lg(x);
22 1233258 : if (l == 1) return 1;
23 1232593 : h = lgcols(x);
24 1232579 : if (h == 1) return 1;
25 4670320 : for (j = l-1; j > 0; j--)
26 22322262 : for (i = h-1; i > 0; i--)
27 18884360 : if (typ(gcoeff(x,i,j)) != t_INT) return 0;
28 1074823 : 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 8351427 : RgMrow_zc_mul_i(GEN x, GEN y, long c, long i)
64 : {
65 8351427 : pari_sp av = avma;
66 8351427 : GEN s = NULL;
67 : long j;
68 835297806 : for (j=1; j<c; j++)
69 : {
70 826949523 : long t = y[j];
71 826949523 : if (!t) continue;
72 29016860 : if (!s) { s = gmulgs(gcoeff(x,i,j),t); continue; }
73 20880814 : switch(t)
74 : {
75 15092642 : case 1: s = gadd(s, gcoeff(x,i,j)); break;
76 239922 : case -1: s = gsub(s, gcoeff(x,i,j)); break;
77 5548250 : default: s = gadd(s, gmulgs(gcoeff(x,i,j), t)); break;
78 : }
79 : }
80 8348283 : 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 239837 : RgM_zc_mul_i(GEN x, GEN y, long c, long l)
87 : {
88 239837 : GEN z = cgetg(l,t_COL);
89 : long i;
90 8521216 : for (i = 1; i < l; i++) gel(z,i) = RgMrow_zc_mul_i(x,y,c,i);
91 239833 : 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 57113 : RgM_zm_mul(GEN x, GEN y)
98 : {
99 57113 : long j, c, l = lg(x), ly = lg(y);
100 57113 : GEN z = cgetg(ly, t_MAT);
101 57113 : if (l == 1) return z;
102 57113 : c = lgcols(x);
103 224968 : for (j = 1; j < ly; j++) gel(z,j) = RgM_zc_mul_i(x, gel(y,j), l,c);
104 57112 : return z;
105 : }
106 :
107 : /* x[i,]*y, l = lg(y) > 1 */
108 : static GEN
109 32305586 : RgMrow_ZC_mul_i(GEN x, GEN y, long i, long l)
110 : {
111 32305586 : pari_sp av = avma;
112 32305586 : GEN t = gmul(gcoeff(x,i,1), gel(y,1)); /* l > 1 ! */
113 : long j;
114 5021588583 : for (j=2; j<l; j++)
115 4989603110 : if (signe(gel(y,j))) t = gadd(t, gmul(gcoeff(x,i,j), gel(y,j)));
116 31985473 : 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 1928449 : RgM_ZC_mul_i(GEN x, GEN y, long lx, long l)
122 : {
123 1928449 : GEN z = cgetg(l,t_COL);
124 : long i;
125 34235001 : for (i=1; i<l; i++) gel(z,i) = RgMrow_ZC_mul_i(x,y,i,lx);
126 1925818 : return z;
127 : }
128 :
129 : GEN
130 1928483 : RgM_ZM_mul_worker(GEN y, GEN x)
131 1928483 : { return RgM_ZC_mul_i(x, y, lg(x), lgcols(x)); }
132 :
133 : /* mostly useful when y is sparse */
134 : GEN
135 354882 : RgM_ZM_mul(GEN x, GEN y)
136 : {
137 354882 : pari_sp av = avma;
138 : GEN worker;
139 354882 : if (lg(x) == 1) return cgetg(lg(y), t_MAT);
140 354882 : worker = snm_closure(is_entry("_RgM_ZM_mul_worker"),mkvec(x));
141 354888 : return gerepileupto(av, gen_parapply(worker, y));
142 : }
143 :
144 : static GEN
145 210025 : RgV_zc_mul_i(GEN x, GEN y, long l)
146 : {
147 : long i;
148 210025 : GEN z = gen_0;
149 210025 : pari_sp av = avma;
150 8213714 : for (i = 1; i < l; i++) z = gadd(z, gmulgs(gel(x,i), y[i]));
151 210024 : 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 34216 : RgV_zm_mul(GEN x, GEN y)
158 : {
159 34216 : long j, l = lg(x), ly = lg(y);
160 34216 : GEN z = cgetg(ly, t_VEC);
161 160856 : for (j = 1; j < ly; j++) gel(z,j) = RgV_zc_mul_i(x, gel(y,j), l);
162 34215 : return z;
163 : }
164 :
165 : /* scalar product x.x */
166 : GEN
167 9449959 : RgV_dotsquare(GEN x)
168 : {
169 9449959 : long i, lx = lg(x);
170 9449959 : pari_sp av = avma;
171 : GEN z;
172 9449959 : if (lx == 1) return gen_0;
173 9449959 : z = gsqr(gel(x,1));
174 21171017 : for (i=2; i<lx; i++)
175 : {
176 11721281 : z = gadd(z, gsqr(gel(x,i)));
177 11721252 : 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 9449736 : return gerepileupto(av,z);
184 : }
185 :
186 : /* scalar product x.y, lx = lg(x) = lg(y) */
187 : static GEN
188 13620969 : RgV_dotproduct_i(GEN x, GEN y, long lx)
189 : {
190 13620969 : pari_sp av = avma;
191 : long i;
192 : GEN z;
193 13620969 : if (lx == 1) return gen_0;
194 13620570 : z = gmul(gel(x,1),gel(y,1));
195 219050376 : for (i=2; i<lx; i++)
196 : {
197 205431354 : z = gadd(z, gmul(gel(x,i), gel(y,i)));
198 205431834 : 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 13619022 : return gerepileupto(av,z);
205 : }
206 : GEN
207 1013399 : RgV_dotproduct(GEN x,GEN y)
208 : {
209 1013399 : if (x == y) return RgV_dotsquare(x);
210 1013399 : return RgV_dotproduct_i(x, y, lg(x));
211 : }
212 : /* v[1] + ... + v[lg(v)-1] */
213 : GEN
214 3032830 : RgV_sum(GEN v)
215 : {
216 : GEN p;
217 3032830 : long i, l = lg(v);
218 3032830 : if (l == 1) return gen_0;
219 7381328 : p = gel(v,1); for (i=2; i<l; i++) p = gadd(p, gel(v,i));
220 3032787 : 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 3336121 : _gmul(void *data, GEN x, GEN y)
264 3336121 : { (void)data; return gmul(x,y); }
265 :
266 : GEN
267 453210 : RgV_prod(GEN x)
268 : {
269 453210 : 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 7623 : RgM_Rg_add(GEN x, GEN y)
276 : {
277 7623 : long l = lg(x), i, j;
278 7623 : GEN z = cgetg(l,t_MAT);
279 :
280 7623 : if (l==1) return z;
281 7623 : if (l != lgcols(x)) pari_err_OP("+", x, y);
282 57358 : for (i=1; i<l; i++)
283 : {
284 49735 : GEN zi = cgetg(l,t_COL), xi = gel(x,i);
285 49735 : gel(z,i) = zi;
286 2057170 : for (j=1; j<l; j++)
287 2007435 : gel(zi,j) = i==j? gadd(y,gel(xi,j)): gcopy(gel(xi,j));
288 : }
289 7623 : 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 314460 : RgM_Rg_sub_shallow(GEN x, GEN y)
327 : {
328 314460 : long l = lg(x), i, j;
329 314460 : GEN z = cgetg(l,t_MAT);
330 :
331 314458 : if (l==1) return z;
332 314458 : if (l != lgcols(x)) pari_err_OP( "-", x, y);
333 2029643 : for (i=1; i<l; i++)
334 : {
335 1715214 : GEN zi = cgetg(l,t_COL), xi = gel(x,i);
336 1715222 : gel(z,i) = zi;
337 27273821 : for (j=1; j<l; j++) gel(zi,j) = gel(xi,j);
338 1715222 : gel(zi,i) = gsub(gel(zi,i), y);
339 : }
340 314429 : return z;
341 : }
342 :
343 : GEN
344 2390850 : RgC_Rg_add(GEN x, GEN y)
345 : {
346 2390850 : long k, lx = lg(x);
347 2390850 : GEN z = cgetg(lx, t_COL);
348 2390851 : if (lx == 1)
349 : {
350 0 : if (isintzero(y)) return z;
351 0 : pari_err_TYPE2("+",x,y);
352 : }
353 2390851 : gel(z,1) = gadd(y,gel(x,1));
354 6746918 : for (k = 2; k < lx; k++) gel(z,k) = gcopy(gel(x,k));
355 2390842 : return z;
356 : }
357 : GEN
358 24029 : RgC_Rg_sub(GEN x, GEN y)
359 : {
360 24029 : long k, lx = lg(x);
361 24029 : GEN z = cgetg(lx, t_COL);
362 24029 : if (lx == 1)
363 : {
364 0 : if (isintzero(y)) return z;
365 0 : pari_err_TYPE2("-",x,y);
366 : }
367 24029 : gel(z,1) = gsub(gel(x,1), y);
368 69113 : for (k = 2; k < lx; k++) gel(z,k) = gcopy(gel(x,k));
369 24029 : return z;
370 : }
371 : /* a - x */
372 : GEN
373 32759 : Rg_RgC_sub(GEN a, GEN x)
374 : {
375 32759 : long k, lx = lg(x);
376 32759 : GEN z = cgetg(lx,t_COL);
377 32759 : if (lx == 1)
378 : {
379 0 : if (isintzero(a)) return z;
380 0 : pari_err_TYPE2("-",a,x);
381 : }
382 32759 : gel(z,1) = gsub(a, gel(x,1));
383 74745 : for (k = 2; k < lx; k++) gel(z,k) = gneg(gel(x,k));
384 32759 : return z;
385 : }
386 :
387 : static GEN
388 18683207 : RgC_add_i(GEN x, GEN y, long lx)
389 : {
390 18683207 : GEN A = cgetg(lx, t_COL);
391 : long i;
392 142201076 : for (i=1; i<lx; i++) gel(A,i) = gadd(gel(x,i), gel(y,i));
393 18682885 : return A;
394 : }
395 : GEN
396 15419974 : RgC_add(GEN x, GEN y) { return RgC_add_i(x, y, lg(x)); }
397 : GEN
398 2546889 : RgV_add(GEN x, GEN y)
399 10435388 : { pari_APPLY_type(t_VEC, gadd(gel(x,i), gel(y,i))) }
400 :
401 : static GEN
402 5754281 : RgC_sub_i(GEN x, GEN y, long lx)
403 : {
404 : long i;
405 5754281 : GEN A = cgetg(lx, t_COL);
406 634341351 : for (i=1; i<lx; i++) gel(A,i) = gsub(gel(x,i), gel(y,i));
407 5753036 : return A;
408 : }
409 : GEN
410 4846684 : RgC_sub(GEN x, GEN y) { return RgC_sub_i(x, y, lg(x)); }
411 : GEN
412 356526 : RgV_sub(GEN x, GEN y)
413 3032044 : { pari_APPLY_type(t_VEC, gsub(gel(x,i), gel(y,i))) }
414 :
415 : GEN
416 690767 : RgM_add(GEN x, GEN y)
417 : {
418 690767 : long lx = lg(x), l, j;
419 : GEN z;
420 690767 : if (lx == 1) return cgetg(1, t_MAT);
421 690767 : z = cgetg(lx, t_MAT); l = lgcols(x);
422 3953999 : for (j = 1; j < lx; j++) gel(z,j) = RgC_add_i(gel(x,j), gel(y,j), l);
423 690767 : return z;
424 : }
425 : GEN
426 173205 : RgM_sub(GEN x, GEN y)
427 : {
428 173205 : long lx = lg(x), l, j;
429 : GEN z;
430 173205 : if (lx == 1) return cgetg(1, t_MAT);
431 173205 : z = cgetg(lx, t_MAT); l = lgcols(x);
432 1080839 : for (j = 1; j < lx; j++) gel(z,j) = RgC_sub_i(gel(x,j), gel(y,j), l);
433 173200 : return z;
434 : }
435 :
436 : static GEN
437 3920438 : RgC_neg_i(GEN x, long lx)
438 : {
439 : long i;
440 3920438 : GEN y = cgetg(lx, t_COL);
441 29928715 : for (i=1; i<lx; i++) gel(y,i) = gneg(gel(x,i));
442 3920439 : return y;
443 : }
444 : GEN
445 932309 : 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 542310 : RgM_neg(GEN x)
451 : {
452 542310 : long i, hx, lx = lg(x);
453 542310 : GEN y = cgetg(lx, t_MAT);
454 542310 : if (lx == 1) return y;
455 542303 : hx = lgcols(x);
456 3530432 : for (i=1; i<lx; i++) gel(y,i) = RgC_neg_i(gel(x,i), hx);
457 542303 : return y;
458 : }
459 :
460 : GEN
461 4735512 : RgV_RgC_mul(GEN x, GEN y)
462 : {
463 4735512 : long lx = lg(x);
464 4735512 : if (lx != lg(y)) pari_err_OP("operation 'RgV_RgC_mul'", x, y);
465 4735428 : 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 115959228 : RgMrow_RgC_mul_i(GEN x, GEN y, long i, long l)
494 : {
495 115959228 : pari_sp av = avma;
496 115959228 : GEN t = gmul(gcoeff(x,i,1), gel(y,1)); /* l > 1 ! */
497 : long j;
498 987631715 : for (j=2; j<l; j++)
499 : {
500 871681623 : GEN c = gcoeff(x,i,j);
501 871681623 : if (!isintzero(c)) t = gadd(t, gmul(c, gel(y,j)));
502 : }
503 115950092 : 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 : static GEN
536 31184662 : RgM_RgC_mul_fast(GEN x, GEN y)
537 : {
538 : GEN p, pol;
539 : long pa;
540 31184662 : long t = RgM_RgC_type(x,y, &p,&pol,&pa);
541 31184902 : switch(t)
542 : {
543 8835435 : case t_INT: return ZM_ZC_mul(x,y);
544 126811 : case t_FRAC: return QM_QC_mul(x,y);
545 91 : case t_FFELT: return FFM_FFC_mul(x, y, pol);
546 28 : case t_INTMOD: return RgM_RgC_mul_FpM(x, y, p);
547 27 : case RgX_type_code(t_POLMOD, t_INTMOD):
548 27 : return RgM_RgC_mul_FqM(x, y, pol, p);
549 22222510 : default: return NULL;
550 : }
551 : }
552 :
553 : /* compatible t_MAT * t_COL, lx = lg(x) = lg(y) > 1, l = lgcols(x) */
554 : static GEN
555 34021973 : RgM_RgC_mul_i(GEN x, GEN y, long lx, long l)
556 : {
557 34021973 : GEN z = cgetg(l,t_COL);
558 : long i;
559 149971363 : for (i=1; i<l; i++) gel(z,i) = RgMrow_RgC_mul_i(x,y,i,lx);
560 34017269 : return z;
561 : }
562 :
563 : GEN
564 31184662 : RgM_RgC_mul(GEN x, GEN y)
565 : {
566 31184662 : long lx = lg(x);
567 : GEN z;
568 31184662 : if (lx != lg(y)) pari_err_OP("operation 'RgM_RgC_mul'", x,y);
569 31184668 : if (lx == 1) return cgetg(1,t_COL);
570 31184668 : z = RgM_RgC_mul_fast(x, y);
571 31184898 : if (z) return z;
572 22222505 : return RgM_RgC_mul_i(x, y, lx, lgcols(x));
573 : }
574 :
575 : GEN
576 296033 : RgV_RgM_mul(GEN x, GEN y)
577 : {
578 296033 : long i, lx, ly = lg(y);
579 : GEN z;
580 296033 : if (ly == 1) return cgetg(1,t_VEC);
581 296026 : lx = lg(x);
582 296026 : if (lx != lgcols(y)) pari_err_OP("operation 'RgV_RgM_mul'", x,y);
583 296019 : z = cgetg(ly, t_VEC);
584 3424920 : for (i=1; i<ly; i++) gel(z,i) = RgV_dotproduct_i(x, gel(y,i), lx);
585 295929 : return z;
586 : }
587 :
588 : static GEN
589 57148 : RgM_mul_FpM_i(GEN x, GEN y, GEN p)
590 : {
591 57148 : if (lgefint(p) == 3)
592 : {
593 57076 : ulong pp = uel(p, 2);
594 57076 : if (pp==2) return F2m_to_mod(F2m_mul(RgM_to_F2m(x), RgM_to_F2m(y)));
595 57027 : if (pp==3) return F3m_to_mod(F3m_mul(RgM_to_F3m(x), RgM_to_F3m(y)));
596 56873 : return Flm_to_mod(Flm_mul(RgM_to_Flm(x, pp), RgM_to_Flm(y, pp), pp), pp);
597 : }
598 72 : return FpM_to_mod(FpM_mul(RgM_to_FpM(x,p), RgM_to_FpM(y,p), p), p);
599 : }
600 : static GEN
601 57148 : RgM_mul_FpM(GEN x, GEN y, GEN p)
602 57148 : { pari_sp av = avma; return gerepileupto(av, RgM_mul_FpM_i(x, y, p)); }
603 :
604 : static GEN
605 26005 : RgM_mul_FqM(GEN x, GEN y, GEN pol, GEN p)
606 : {
607 26005 : pari_sp av = avma;
608 26005 : GEN b, T = RgX_to_FpX(pol, p);
609 26005 : if (signe(T) == 0) pari_err_OP("*", x, y);
610 26005 : b = FqM_mul(RgM_to_FqM(x, T, p), RgM_to_FqM(y, T, p), T, p);
611 26005 : return gerepileupto(av, FqM_to_mod(b, T, p));
612 : }
613 :
614 : static GEN
615 12516 : RgM_liftred(GEN x, GEN T)
616 12516 : { return RgXQM_red(liftpol_shallow(x), T); }
617 :
618 : static GEN
619 1463 : RgM_mul_ZXQM(GEN x, GEN y, GEN T)
620 : {
621 1463 : pari_sp av = avma;
622 1463 : GEN b = ZXQM_mul(RgM_liftred(x,T), RgM_liftred(y, T), T);
623 1463 : return gerepilecopy(av, QXQM_to_mod_shallow(b,T));
624 : }
625 :
626 : static GEN
627 133 : RgM_sqr_ZXQM(GEN x, GEN T)
628 : {
629 133 : pari_sp av = avma;
630 133 : GEN b = ZXQM_sqr(RgM_liftred(x, T), T);
631 133 : return gerepilecopy(av, QXQM_to_mod_shallow(b,T));
632 : }
633 :
634 : static GEN
635 4725 : RgM_mul_QXQM(GEN x, GEN y, GEN T)
636 : {
637 4725 : pari_sp av = avma;
638 4725 : GEN b = QXQM_mul(RgM_liftred(x, T), RgM_liftred(y, T), T);
639 4725 : return gerepilecopy(av, QXQM_to_mod_shallow(b,T));
640 : }
641 :
642 : static GEN
643 7 : RgM_sqr_QXQM(GEN x, GEN T)
644 : {
645 7 : pari_sp av = avma;
646 7 : GEN b = QXQM_sqr(RgM_liftred(x, T), T);
647 7 : return gerepilecopy(av, QXQM_to_mod_shallow(b,T));
648 : }
649 :
650 : INLINE int
651 4774 : RgX_is_monic_ZX(GEN pol)
652 4774 : { return RgX_is_ZX(pol) && ZX_is_monic(pol); }
653 :
654 : static GEN
655 7575452 : RgM_mul_fast(GEN x, GEN y)
656 : {
657 : GEN p, pol;
658 : long pa;
659 7575452 : long t = RgM_type2(x,y, &p,&pol,&pa);
660 7575508 : switch(t)
661 : {
662 2689608 : case t_INT: return ZM_mul(x,y);
663 117138 : case t_FRAC: return QM_mul(x,y);
664 4389 : case t_FFELT: return FFM_mul(x, y, pol);
665 57085 : case t_INTMOD: return RgM_mul_FpM(x, y, p);
666 1470 : case RgX_type_code(t_POLMOD, t_INT):
667 1470 : return ZX_is_monic(pol)? RgM_mul_ZXQM(x, y, pol): NULL;
668 4746 : case RgX_type_code(t_POLMOD, t_FRAC):
669 4746 : return RgX_is_monic_ZX(pol)? RgM_mul_QXQM(x, y, pol): NULL;
670 26003 : case RgX_type_code(t_POLMOD, t_INTMOD):
671 26003 : return RgM_mul_FqM(x, y, pol, p);
672 4675069 : default: return NULL;
673 : }
674 : }
675 :
676 : static GEN
677 1365 : RgM_sqr_fast(GEN x)
678 : {
679 : GEN p, pol;
680 : long pa;
681 1365 : long t = RgM_type(x, &p,&pol,&pa);
682 1365 : switch(t)
683 : {
684 231 : case t_INT: return ZM_sqr(x);
685 700 : case t_FRAC: return QM_sqr(x);
686 105 : case t_FFELT: return FFM_mul(x, x, pol);
687 63 : case t_INTMOD: return RgM_mul_FpM(x, x, p);
688 140 : case RgX_type_code(t_POLMOD, t_INT):
689 140 : return ZX_is_monic(pol)? RgM_sqr_ZXQM(x, pol): NULL;
690 28 : case RgX_type_code(t_POLMOD, t_FRAC):
691 28 : return RgX_is_monic_ZX(pol)? RgM_sqr_QXQM(x, pol): NULL;
692 0 : case RgX_type_code(t_POLMOD, t_INTMOD):
693 0 : return RgM_mul_FqM(x, x, pol, p);
694 98 : default: return NULL;
695 : }
696 : }
697 :
698 : /* lx, ly > 1 */
699 : static GEN
700 4675095 : RgM_mul_i(GEN x, GEN y, long lx, long ly)
701 : {
702 4675095 : GEN z = cgetg(ly, t_MAT);
703 4675089 : long j, l = lgcols(x);
704 16474374 : for (j = 1; j < ly; j++) gel(z,j) = RgM_RgC_mul_i(x, gel(y,j), lx, l);
705 4675041 : return z;
706 : }
707 : GEN
708 7591229 : RgM_mul(GEN x, GEN y)
709 : {
710 7591229 : long lx, ly = lg(y);
711 : GEN z;
712 7591229 : if (ly == 1) return cgetg(1,t_MAT);
713 7575479 : lx = lg(x);
714 7575479 : if (lx != lgcols(y)) pari_err_OP("operation 'RgM_mul'", x,y);
715 7575478 : if (lx == 1) return zeromat(0,ly-1);
716 7575450 : z = RgM_mul_fast(x, y);
717 7575506 : if (z) return z;
718 4675094 : return RgM_mul_i(x, y, lx, ly);
719 : }
720 :
721 : GEN
722 1400 : RgM_sqr(GEN x)
723 : {
724 1400 : long j, lx = lg(x);
725 : GEN z;
726 1400 : if (lx == 1) return cgetg(1, t_MAT);
727 1365 : if (lx != lgcols(x)) pari_err_OP("operation 'RgM_mul'", x,x);
728 1365 : z = RgM_sqr_fast(x);
729 1365 : if (z) return z;
730 126 : z = cgetg(lx, t_MAT);
731 497 : for (j=1; j<lx; j++) gel(z,j) = RgM_RgC_mul_i(x, gel(x,j), lx, lx);
732 126 : return z;
733 : }
734 :
735 : /* assume result is symmetric */
736 : GEN
737 0 : RgM_multosym(GEN x, GEN y)
738 : {
739 0 : long j, lx, ly = lg(y);
740 : GEN M;
741 0 : if (ly == 1) return cgetg(1,t_MAT);
742 0 : lx = lg(x);
743 0 : if (lx != lgcols(y)) pari_err_OP("operation 'RgM_multosym'", x,y);
744 0 : if (lx == 1) return cgetg(1,t_MAT);
745 0 : if (ly != lgcols(x)) pari_err_OP("operation 'RgM_multosym'", x,y);
746 0 : M = cgetg(ly, t_MAT);
747 0 : for (j=1; j<ly; j++)
748 : {
749 0 : GEN z = cgetg(ly,t_COL), yj = gel(y,j);
750 : long i;
751 0 : for (i=1; i<j; i++) gel(z,i) = gcoeff(M,j,i);
752 0 : for (i=j; i<ly; i++)gel(z,i) = RgMrow_RgC_mul_i(x,yj,i,lx);
753 0 : gel(M,j) = z;
754 : }
755 0 : return M;
756 : }
757 : /* x~ * y, assuming result is symmetric */
758 : GEN
759 7663 : RgM_transmultosym(GEN x, GEN y)
760 : {
761 7663 : long i, j, l, ly = lg(y);
762 : GEN M;
763 7663 : if (ly == 1) return cgetg(1,t_MAT);
764 7663 : if (lg(x) != ly) pari_err_OP("operation 'RgM_transmultosym'", x,y);
765 7663 : l = lgcols(y);
766 7663 : if (lgcols(x) != l) pari_err_OP("operation 'RgM_transmultosym'", x,y);
767 7663 : M = cgetg(ly, t_MAT);
768 30820 : for (i=1; i<ly; i++)
769 : {
770 23157 : GEN xi = gel(x,i), c = cgetg(ly,t_COL);
771 23157 : gel(M,i) = c;
772 47490 : for (j=1; j<i; j++)
773 24333 : gcoeff(M,i,j) = gel(c,j) = RgV_dotproduct_i(xi,gel(y,j),l);
774 23157 : gel(c,i) = RgV_dotproduct_i(xi,gel(y,i),l);
775 : }
776 7663 : return M;
777 : }
778 : /* x~ * y */
779 : GEN
780 0 : RgM_transmul(GEN x, GEN y)
781 : {
782 0 : long i, j, l, lx, ly = lg(y);
783 : GEN M;
784 0 : if (ly == 1) return cgetg(1,t_MAT);
785 0 : lx = lg(x);
786 0 : l = lgcols(y);
787 0 : if (lgcols(x) != l) pari_err_OP("operation 'RgM_transmul'", x,y);
788 0 : M = cgetg(ly, t_MAT);
789 0 : for (i=1; i<ly; i++)
790 : {
791 0 : GEN yi = gel(y,i), c = cgetg(lx,t_COL);
792 0 : gel(M,i) = c;
793 0 : for (j=1; j<lx; j++) gel(c,j) = RgV_dotproduct_i(yi,gel(x,j),l);
794 : }
795 0 : return M;
796 : }
797 :
798 : GEN
799 4695734 : gram_matrix(GEN x)
800 : {
801 4695734 : long i,j, l, lx = lg(x);
802 : GEN M;
803 4695734 : if (!is_matvec_t(typ(x))) pari_err_TYPE("gram",x);
804 4695730 : if (lx == 1) return cgetg(1,t_MAT);
805 4695716 : l = lgcols(x);
806 4695714 : M = cgetg(lx,t_MAT);
807 14087089 : for (i=1; i<lx; i++)
808 : {
809 9391382 : GEN xi = gel(x,i), c = cgetg(lx,t_COL);
810 9391379 : gel(M,i) = c;
811 14087060 : for (j=1; j<i; j++)
812 4695688 : gcoeff(M,i,j) = gel(c,j) = RgV_dotproduct_i(xi,gel(x,j),l);
813 9391372 : gel(c,i) = RgV_dotsquare(xi);
814 : }
815 4695707 : return M;
816 : }
817 :
818 : static GEN
819 3717 : _RgM_add(void *E, GEN x, GEN y) { (void)E; return RgM_add(x, y); }
820 :
821 : static GEN
822 0 : _RgM_sub(void *E, GEN x, GEN y) { (void)E; return RgM_sub(x, y); }
823 :
824 : static GEN
825 5600 : _RgM_cmul(void *E, GEN P, long a, GEN x) { (void)E; return RgM_Rg_mul(x,gel(P,a+2)); }
826 :
827 : static GEN
828 231 : _RgM_sqr(void *E, GEN x) { (void) E; return RgM_sqr(x); }
829 :
830 : static GEN
831 735 : _RgM_mul(void *E, GEN x, GEN y) { (void) E; return RgM_mul(x, y); }
832 :
833 : static GEN
834 4039 : _RgM_one(void *E) { long *n = (long*) E; return matid(*n); }
835 :
836 : static GEN
837 0 : _RgM_zero(void *E) { long *n = (long*) E; return zeromat(*n,*n); }
838 :
839 : static GEN
840 2821 : _RgM_red(void *E, GEN x) { (void)E; return x; }
841 :
842 : static struct bb_algebra RgM_algebra = { _RgM_red, _RgM_add, _RgM_sub,
843 : _RgM_mul, _RgM_sqr, _RgM_one, _RgM_zero };
844 :
845 : /* generates the list of powers of x of degree 0,1,2,...,l*/
846 : GEN
847 168 : RgM_powers(GEN x, long l)
848 : {
849 168 : long n = lg(x)-1;
850 168 : return gen_powers(x,l,1,(void *) &n, &_RgM_sqr, &_RgM_mul, &_RgM_one);
851 : }
852 :
853 : GEN
854 490 : RgX_RgMV_eval(GEN Q, GEN x)
855 : {
856 490 : long n = lg(x)>1 ? lg(gel(x,1))-1:0;
857 490 : return gen_bkeval_powers(Q,degpol(Q),x,(void*)&n,&RgM_algebra,&_RgM_cmul);
858 : }
859 :
860 : GEN
861 1393 : RgX_RgM_eval(GEN Q, GEN x)
862 : {
863 1393 : long n = lg(x)-1;
864 1393 : return gen_bkeval(Q,degpol(Q),x,1,(void*)&n,&RgM_algebra,&_RgM_cmul);
865 : }
866 :
867 : GEN
868 4413016 : RgC_Rg_div(GEN x, GEN y)
869 19050316 : { pari_APPLY_type(t_COL, gdiv(gel(x,i),y)) }
870 :
871 : GEN
872 11844615 : RgC_Rg_mul(GEN x, GEN y)
873 75470419 : { pari_APPLY_type(t_COL, gmul(gel(x,i),y)) }
874 :
875 : GEN
876 31157 : RgV_Rg_mul(GEN x, GEN y)
877 1200766 : { pari_APPLY_type(t_VEC, gmul(gel(x,i),y)) }
878 :
879 : GEN
880 604183 : RgM_Rg_div(GEN X, GEN c) {
881 604183 : long i, j, h, l = lg(X);
882 604183 : GEN A = cgetg(l, t_MAT);
883 604183 : if (l == 1) return A;
884 604183 : h = lgcols(X);
885 2784120 : for (j=1; j<l; j++)
886 : {
887 2179958 : GEN a = cgetg(h, t_COL), x = gel(X, j);
888 17358388 : for (i = 1; i < h; i++) gel(a,i) = gdiv(gel(x,i), c);
889 2179937 : gel(A,j) = a;
890 : }
891 604162 : return A;
892 : }
893 : GEN
894 1248841 : RgM_Rg_mul(GEN X, GEN c) {
895 1248841 : long i, j, h, l = lg(X);
896 1248841 : GEN A = cgetg(l, t_MAT);
897 1248851 : if (l == 1) return A;
898 1248851 : h = lgcols(X);
899 4084272 : for (j=1; j<l; j++)
900 : {
901 2835500 : GEN a = cgetg(h, t_COL), x = gel(X, j);
902 12930312 : for (i = 1; i < h; i++) gel(a,i) = gmul(gel(x,i), c);
903 2835418 : gel(A,j) = a;
904 : }
905 1248772 : return A;
906 : }
907 :
908 : /********************************************************************/
909 : /* */
910 : /* SCALAR TO MATRIX/VECTOR */
911 : /* */
912 : /********************************************************************/
913 : /* fill the square nxn matrix equal to t*Id */
914 : static void
915 12724438 : fill_scalmat(GEN y, GEN t, long n)
916 : {
917 : long i;
918 47966190 : for (i = 1; i <= n; i++)
919 : {
920 35242117 : gel(y,i) = zerocol(n);
921 35241752 : gcoeff(y,i,i) = t;
922 : }
923 12724073 : }
924 :
925 : GEN
926 869702 : scalarmat(GEN x, long n) {
927 869702 : GEN y = cgetg(n+1, t_MAT);
928 869679 : if (!n) return y;
929 869679 : fill_scalmat(y, gcopy(x), n); return y;
930 : }
931 : GEN
932 3154048 : scalarmat_shallow(GEN x, long n) {
933 3154048 : GEN y = cgetg(n+1, t_MAT);
934 3154120 : fill_scalmat(y, x, n); return y;
935 : }
936 : GEN
937 217 : scalarmat_s(long x, long n) {
938 217 : GEN y = cgetg(n+1, t_MAT);
939 217 : if (!n) return y;
940 217 : fill_scalmat(y, stoi(x), n); return y;
941 : }
942 : GEN
943 8701064 : matid(long n) {
944 : GEN y;
945 8701064 : if (n < 0) pari_err_DOMAIN("matid", "size", "<", gen_0, stoi(n));
946 8701057 : y = cgetg(n+1, t_MAT);
947 8701095 : fill_scalmat(y, gen_1, n); return y;
948 : }
949 :
950 : INLINE GEN
951 2265812 : scalarcol_i(GEN x, long n, long c)
952 : {
953 : long i;
954 2265812 : GEN y = cgetg(n+1,t_COL);
955 2265807 : if (!n) return y;
956 2265807 : gel(y,1) = c? gcopy(x): x;
957 6758249 : for (i=2; i<=n; i++) gel(y,i) = gen_0;
958 2265805 : return y;
959 : }
960 :
961 : GEN
962 785693 : scalarcol(GEN x, long n) { return scalarcol_i(x,n,1); }
963 :
964 : GEN
965 1480122 : scalarcol_shallow(GEN x, long n) { return scalarcol_i(x,n,0); }
966 :
967 : int
968 17597 : RgM_isscalar(GEN x, GEN s)
969 : {
970 17597 : long i, j, lx = lg(x);
971 :
972 17597 : if (lx == 1) return 1;
973 17597 : if (lx != lgcols(x)) return 0;
974 17597 : if (!s) s = gcoeff(x,1,1);
975 :
976 49797 : for (j=1; j<lx; j++)
977 : {
978 40505 : GEN c = gel(x,j);
979 154099 : for (i=1; i<j; )
980 120492 : if (!gequal0(gel(c,i++))) return 0;
981 : /* i = j */
982 33607 : if (!gequal(gel(c,i++),s)) return 0;
983 176003 : for ( ; i<lx; )
984 143803 : if (!gequal0(gel(c,i++))) return 0;
985 : }
986 9292 : return 1;
987 : }
988 :
989 : int
990 3664 : RgM_isidentity(GEN x)
991 : {
992 3664 : long i,j, lx = lg(x);
993 :
994 3664 : if (lx == 1) return 1;
995 3664 : if (lx != lgcols(x)) return 0;
996 8798 : for (j=1; j<lx; j++)
997 : {
998 8518 : GEN c = gel(x,j);
999 40052 : for (i=1; i<j; )
1000 33373 : if (!gequal0(gel(c,i++))) return 0;
1001 : /* i = j */
1002 6679 : if (!gequal1(gel(c,i++))) return 0;
1003 49750 : for ( ; i<lx; )
1004 44616 : if (!gequal0(gel(c,i++))) return 0;
1005 : }
1006 280 : return 1;
1007 : }
1008 :
1009 : long
1010 644 : RgC_is_ei(GEN x)
1011 : {
1012 644 : long i, j = 0, l = lg(x);
1013 3472 : for (i = 1; i < l; i++)
1014 : {
1015 2828 : GEN c = gel(x,i);
1016 2828 : if (gequal0(c)) continue;
1017 644 : if (!gequal1(c) || j) return 0;
1018 644 : j = i;
1019 : }
1020 644 : return j;
1021 : }
1022 :
1023 : int
1024 336 : RgM_isdiagonal(GEN x)
1025 : {
1026 336 : long i,j, lx = lg(x);
1027 336 : if (lx == 1) return 1;
1028 336 : if (lx != lgcols(x)) return 0;
1029 :
1030 3220 : for (j=1; j<lx; j++)
1031 : {
1032 2891 : GEN c = gel(x,j);
1033 18452 : for (i=1; i<j; i++)
1034 15561 : if (!gequal0(gel(c,i))) return 0;
1035 18452 : for (i++; i<lx; i++)
1036 15568 : if (!gequal0(gel(c,i))) return 0;
1037 : }
1038 329 : return 1;
1039 : }
1040 : int
1041 315 : isdiagonal(GEN x) { return (typ(x)==t_MAT) && RgM_isdiagonal(x); }
1042 :
1043 : GEN
1044 22358 : RgM_det_triangular(GEN mat)
1045 : {
1046 22358 : long i,l = lg(mat);
1047 : pari_sp av;
1048 : GEN s;
1049 :
1050 22358 : if (l<3) return l<2? gen_1: gcopy(gcoeff(mat,1,1));
1051 20671 : av = avma; s = gcoeff(mat,1,1);
1052 121695 : for (i=2; i<l; i++) s = gmul(s,gcoeff(mat,i,i));
1053 20671 : return av==avma? gcopy(s): gerepileupto(av,s);
1054 : }
1055 :
1056 : GEN
1057 478032 : RgV_kill0(GEN v)
1058 : {
1059 : long i, l;
1060 478032 : GEN w = cgetg_copy(v, &l);
1061 130113800 : for (i = 1; i < l; i++)
1062 : {
1063 129636000 : GEN a = gel(v,i);
1064 129636000 : gel(w,i) = gequal0(a) ? NULL: a;
1065 : }
1066 477800 : return w;
1067 : }
|
