## "Fossies" - the Fresh Open Source Software Archive

### Member "pari-2.13.1/src/basemath/RgV.c" (14 Jan 2021, 23737 Bytes) of package /linux/misc/pari-2.13.1.tar.gz:

As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "RgV.c" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 2.13.0_vs_2.13.1.

```    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 RgM_is_ZM(GEN x)
20 hlOpen(20,1);{
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 hlClose(1, 29);}
30
31 int
32 RgM_is_QM(GEN x)
33 hlOpen(33,1);{
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 hlClose(2, 42);}
43
44 int
45 RgV_is_ZMV(GEN V)
46 hlOpen(46,1);{
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 hlClose(3, 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 RgMrow_zc_mul_i(GEN x, GEN y, long c, long i)
64 hlOpen(64,1);{
65   pari_sp av = avma;
66   GEN s = NULL;
67   long j;
68   for (j=1; j<c; j++)
69   hlOpen(69,2);{
70     long t = y[j];
71     if (!t) continue;
72     if (!s) hlOpen(72,3);{ s = gmulgs(gcoeff(x,i,j),t); continue; hlClose(6, 72);}
73     switch(t)
74     hlOpen(74,3);{
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     hlClose(7, 78);}
79   hlClose(5, 79);}
80   return s? gerepileupto(av, s): gc_const(av, gen_0);
81 hlClose(4, 81);}
82 GEN
83 RgMrow_zc_mul(GEN x, GEN y, long i) hlOpen(83,1);{ return RgMrow_zc_mul_i(x,y,lg(y),i); hlClose(8, 83);}
84 /* x nonempty t_MAT, y a compatible zc (dimension > 0). */
85 static GEN
86 RgM_zc_mul_i(GEN x, GEN y, long c, long l)
87 hlOpen(87,1);{
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 hlClose(9, 92);}
93 GEN
94 RgM_zc_mul(GEN x, GEN y) hlOpen(94,1);{ return RgM_zc_mul_i(x,y, lg(x), lgcols(x)); hlClose(10, 94);}
95 /* x t_MAT, y a compatible zm (dimension > 0). */
96 GEN
97 RgM_zm_mul(GEN x, GEN y)
98 hlOpen(98,1);{
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 hlClose(11, 105);}
106
107 /* x[i,]*y, l = lg(y) > 1 */
108 static GEN
109 RgMrow_ZC_mul_i(GEN x, GEN y, long i, long l)
110 hlOpen(110,1);{
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 hlClose(12, 117);}
118
119 /* compatible t_MAT * t_COL, lx = lg(x) = lg(y) > 1, l = lgcols(x) */
120 static GEN
121 RgM_ZC_mul_i(GEN x, GEN y, long lx, long l)
122 hlOpen(122,1);{
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 hlClose(13, 127);}
128
129 /* mostly useful when y is sparse */
130 GEN
131 RgM_ZM_mul(GEN x, GEN y)
132 hlOpen(132,1);{
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 hlClose(14, 139);}
140
141 static GEN
142 RgV_zc_mul_i(GEN x, GEN y, long l)
143 hlOpen(143,1);{
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 hlClose(15, 149);}
150 GEN
151 RgV_zc_mul(GEN x, GEN y) hlOpen(151,1);{ return RgV_zc_mul_i(x, y, lg(x)); hlClose(16, 151);}
152
153 GEN
154 RgV_zm_mul(GEN x, GEN y)
155 hlOpen(155,1);{
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 hlClose(17, 160);}
161
162 /* scalar product x.x */
163 GEN
164 RgV_dotsquare(GEN x)
165 hlOpen(165,1);{
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   hlOpen(172,2);{
173     z = gadd(z, gsqr(gel(x,i)));
174     if (gc_needed(av,3))
175     hlOpen(175,3);{
176       if(DEBUGMEM>1) pari_warn(warnmem,"RgV_dotsquare, i = %ld",i);
177       z = gerepileupto(av, z);
178     hlClose(20, 178);}
179   hlClose(19, 179);}
180   return gerepileupto(av,z);
181 hlClose(18, 181);}
182
183 /* scalar product x.y, lx = lg(x) = lg(y) */
184 static GEN
185 RgV_dotproduct_i(GEN x, GEN y, long lx)
186 hlOpen(186,1);{
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   hlOpen(193,2);{
194     z = gadd(z, gmul(gel(x,i), gel(y,i)));
195     if (gc_needed(av,3))
196     hlOpen(196,3);{
197       if(DEBUGMEM>1) pari_warn(warnmem,"RgV_dotproduct, i = %ld",i);
198       z = gerepileupto(av, z);
199     hlClose(23, 199);}
200   hlClose(22, 200);}
201   return gerepileupto(av,z);
202 hlClose(21, 202);}
203 GEN
204 RgV_dotproduct(GEN x,GEN y)
205 hlOpen(205,1);{
206   if (x == y) return RgV_dotsquare(x);
207   return RgV_dotproduct_i(x, y, lg(x));
208 hlClose(24, 208);}
209 /* v[1] + ... + v[lg(v)-1] */
210 GEN
211 RgV_sum(GEN v)
212 hlOpen(212,1);{
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 hlClose(25, 218);}
219 /* v[1] + ... + v[n]. Assume lg(v) > n. */
220 GEN
221 RgV_sumpart(GEN v, long n)
222 hlOpen(222,1);{
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 hlClose(26, 228);}
229 /* v[m] + ... + v[n]. Assume lg(v) > n, m > 0. */
230 GEN
231 RgV_sumpart2(GEN v, long m, long n)
232 hlOpen(232,1);{
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 hlClose(27, 238);}
239 GEN
240 RgM_sumcol(GEN A)
241 hlOpen(241,1);{
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   hlOpen(250,2);{
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   hlClose(29, 255);}
256   return v;
257 hlClose(28, 257);}
258
259 static GEN
260 _gmul(void *data, GEN x, GEN y)
261 hlOpen(261,1);{ (void)data; return gmul(x,y); hlClose(30, 261);}
262
263 GEN
264 RgV_prod(GEN x)
265 hlOpen(265,1);{
266   return gen_product(x, NULL, _gmul);
267 hlClose(31, 267);}
268
269 /*                    ADDITION SCALAR + MATRIX                     */
270 /* x square matrix, y scalar; create the square matrix x + y*Id */
271 GEN
272 RgM_Rg_add(GEN x, GEN y)
273 hlOpen(273,1);{
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   hlOpen(281,2);{
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   hlClose(33, 286);}
287   return z;
288 hlClose(32, 288);}
289 GEN
290 RgM_Rg_sub(GEN x, GEN y)
291 hlOpen(291,1);{
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   hlOpen(299,2);{
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   hlClose(35, 304);}
305   return z;
306 hlClose(34, 306);}
307 GEN
308 RgM_Rg_add_shallow(GEN x, GEN y)
309 hlOpen(309,1);{
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   hlOpen(316,2);{
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   hlClose(37, 321);}
322   return z;
323 hlClose(36, 323);}
324 GEN
325 RgM_Rg_sub_shallow(GEN x, GEN y)
326 hlOpen(326,1);{
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   hlOpen(333,2);{
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   hlClose(39, 338);}
339   return z;
340 hlClose(38, 340);}
341
342 GEN
343 RgC_Rg_add(GEN x, GEN y)
344 hlOpen(344,1);{
345   long k, lx = lg(x);
346   GEN z = cgetg(lx, t_COL);
347   if (lx == 1)
348   hlOpen(348,2);{
349     if (isintzero(y)) return z;
350     pari_err_TYPE2("+",x,y);
351   hlClose(41, 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 hlClose(40, 355);}
356 GEN
357 RgC_Rg_sub(GEN x, GEN y)
358 hlOpen(358,1);{
359   long k, lx = lg(x);
360   GEN z = cgetg(lx, t_COL);
361   if (lx == 1)
362   hlOpen(362,2);{
363     if (isintzero(y)) return z;
364     pari_err_TYPE2("-",x,y);
365   hlClose(43, 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 hlClose(42, 369);}
370 /* a - x */
371 GEN
372 Rg_RgC_sub(GEN a, GEN x)
373 hlOpen(373,1);{
374   long k, lx = lg(x);
375   GEN z = cgetg(lx,t_COL);
376   if (lx == 1)
377   hlOpen(377,2);{
378     if (isintzero(a)) return z;
379     pari_err_TYPE2("-",a,x);
380   hlClose(45, 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 hlClose(44, 384);}
385
386 static GEN
387 RgC_add_i(GEN x, GEN y, long lx)
388 hlOpen(388,1);{
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 hlClose(46, 393);}
394 GEN
395 RgC_add(GEN x, GEN y) hlOpen(395,1);{ return RgC_add_i(x, y, lg(x)); hlClose(47, 395);}
396 GEN
397 RgV_add(GEN x, GEN y)
398 hlOpen(398,1);{ pari_APPLY_type(t_VEC, gadd(gel(x,i), gel(y,i))) hlClose(48, 398);}
399
400 static GEN
401 RgC_sub_i(GEN x, GEN y, long lx)
402 hlOpen(402,1);{
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 hlClose(49, 407);}
408 GEN
409 RgC_sub(GEN x, GEN y) hlOpen(409,1);{ return RgC_sub_i(x, y, lg(x)); hlClose(50, 409);}
410 GEN
411 RgV_sub(GEN x, GEN y)
412 hlOpen(412,1);{ pari_APPLY_type(t_VEC, gsub(gel(x,i), gel(y,i))) hlClose(51, 412);}
413
414 GEN
415 RgM_add(GEN x, GEN y)
416 hlOpen(416,1);{
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 hlClose(52, 423);}
424 GEN
425 RgM_sub(GEN x, GEN y)
426 hlOpen(426,1);{
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 hlClose(53, 433);}
434
435 static GEN
436 RgC_neg_i(GEN x, long lx)
437 hlOpen(437,1);{
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 hlClose(54, 442);}
443 GEN
444 RgC_neg(GEN x) hlOpen(444,1);{ return RgC_neg_i(x, lg(x)); hlClose(55, 444);}
445 GEN
446 RgV_neg(GEN x)
447 hlOpen(447,1);{ pari_APPLY_type(t_VEC, gneg(gel(x,i))) hlClose(56, 447);}
448 GEN
449 RgM_neg(GEN x)
450 hlOpen(450,1);{
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 hlClose(57, 457);}
458
459 GEN
460 RgV_RgC_mul(GEN x, GEN y)
461 hlOpen(461,1);{
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 hlClose(58, 465);}
466 GEN
467 RgC_RgV_mul(GEN x, GEN y)
468 hlOpen(468,1);{
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 hlClose(59, 473);}
474 GEN
475 RgC_RgM_mul(GEN x, GEN y)
476 hlOpen(476,1);{
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 hlClose(60, 482);}
483 GEN
484 RgM_RgV_mul(GEN x, GEN y)
485 hlOpen(485,1);{
486   if (lg(x) != 2) pari_err_OP("operation 'RgM_RgV_mul'", x,y);
487   return RgC_RgV_mul(gel(x,1), y);
488 hlClose(61, 488);}
489
490 /* x[i,]*y, l = lg(y) > 1 */
491 static GEN
492 RgMrow_RgC_mul_i(GEN x, GEN y, long i, long l)
493 hlOpen(493,1);{
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   hlOpen(498,2);{
499     GEN c = gcoeff(x,i,j);
500     if (!isintzero(c)) t = gadd(t, gmul(c, gel(y,j)));
501   hlClose(63, 501);}
502   return gerepileupto(av,t);
503 hlClose(62, 503);}
504 GEN
505 RgMrow_RgC_mul(GEN x, GEN y, long i)
506 hlOpen(506,1);{ return RgMrow_RgC_mul_i(x, y, i, lg(x)); hlClose(64, 506);}
507
508 static GEN
509 RgM_RgC_mul_FpM(GEN x, GEN y, GEN p)
510 hlOpen(510,1);{
511   pari_sp av = avma;
512   GEN r;
513   if (lgefint(p) == 3)
514   hlOpen(514,2);{
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   hlClose(66, 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 hlClose(65, 522);}
523
524 static GEN
525 RgM_RgC_mul_FqM(GEN x, GEN y, GEN pol, GEN p)
526 hlOpen(526,1);{
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 hlClose(67, 532);}
533
534 #define code(t1,t2) ((t1 << 6) | t2)
535 static GEN
536 RgM_RgC_mul_fast(GEN x, GEN y)
537 hlOpen(537,1);{
538   GEN p, pol;
539   long pa;
540   long t = RgM_RgC_type(x,y, &p,&pol,&pa);
541   switch(t)
542   hlOpen(542,2);{
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   hlClose(69, 550);}
551 hlClose(68, 551);}
552 #undef code
553
554 /* compatible t_MAT * t_COL, lx = lg(x) = lg(y) > 1, l = lgcols(x) */
555 static GEN
556 RgM_RgC_mul_i(GEN x, GEN y, long lx, long l)
557 hlOpen(557,1);{
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 hlClose(70, 562);}
563
564 GEN
565 RgM_RgC_mul(GEN x, GEN y)
566 hlOpen(566,1);{
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 hlClose(71, 574);}
575
576 GEN
577 RgV_RgM_mul(GEN x, GEN y)
578 hlOpen(578,1);{
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 hlClose(72, 587);}
588
589 static GEN
590 RgM_mul_FpM(GEN x, GEN y, GEN p)
591 hlOpen(591,1);{
592   pari_sp av = avma;
593   GEN r;
594   if (lgefint(p) == 3)
595   hlOpen(595,2);{
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   hlClose(74, 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 hlClose(73, 603);}
604
605 static GEN
606 RgM_mul_FqM(GEN x, GEN y, GEN pol, GEN p)
607 hlOpen(607,1);{
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 hlClose(75, 613);}
614
615 static GEN
616 RgM_liftred(GEN x, GEN T)
617 hlOpen(617,1);{ return RgXQM_red(liftpol_shallow(x), T); hlClose(76, 617);}
618
619 static GEN
620 RgM_mul_ZXQM(GEN x, GEN y, GEN T)
621 hlOpen(621,1);{
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 hlClose(77, 625);}
626
627 static GEN
628 RgM_sqr_ZXQM(GEN x, GEN T)
629 hlOpen(629,1);{
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 hlClose(78, 633);}
634
635 static GEN
636 RgM_mul_QXQM(GEN x, GEN y, GEN T)
637 hlOpen(637,1);{
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 hlClose(79, 641);}
642
643 static GEN
644 RgM_sqr_QXQM(GEN x, GEN T)
645 hlOpen(645,1);{
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 hlClose(80, 649);}
650
651 INLINE int
652 RgX_is_monic_ZX(GEN pol)
653 hlOpen(653,1);{ return RgX_is_ZX(pol) && ZX_is_monic(pol); hlClose(81, 653);}
654
655 #define code(t1,t2) ((t1 << 6) | t2)
656 static GEN
657 RgM_mul_fast(GEN x, GEN y)
658 hlOpen(658,1);{
659   GEN p, pol;
660   long pa;
661   long t = RgM_type2(x,y, &p,&pol,&pa);
662   switch(t)
663   hlOpen(663,2);{
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   hlClose(83, 675);}
676 hlClose(82, 676);}
677
678 static GEN
679 RgM_sqr_fast(GEN x)
680 hlOpen(680,1);{
681   GEN p, pol;
682   long pa;
683   long t = RgM_type(x, &p,&pol,&pa);
684   switch(t)
685   hlOpen(685,2);{
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   hlClose(85, 697);}
698 hlClose(84, 698);}
699
700 #undef code
701
702 GEN
703 RgM_mul(GEN x, GEN y)
704 hlOpen(704,1);{
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 hlClose(86, 717);}
718
719 GEN
720 RgM_sqr(GEN x)
721 hlOpen(721,1);{
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 hlClose(87, 731);}
732
733 /* assume result is symmetric */
734 GEN
735 RgM_multosym(GEN x, GEN y)
736 hlOpen(736,1);{
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   hlOpen(746,2);{
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   hlClose(89, 752);}
753   return M;
754 hlClose(88, 754);}
755 /* x~ * y, assuming result is symmetric */
756 GEN
757 RgM_transmultosym(GEN x, GEN y)
758 hlOpen(758,1);{
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   hlOpen(767,2);{
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   hlClose(91, 773);}
774   return M;
775 hlClose(90, 775);}
776 /* x~ * y */
777 GEN
778 RgM_transmul(GEN x, GEN y)
779 hlOpen(779,1);{
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   hlOpen(788,2);{
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   hlClose(93, 792);}
793   return M;
794 hlClose(92, 794);}
795
796 GEN
797 gram_matrix(GEN x)
798 hlOpen(798,1);{
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   hlOpen(806,2);{
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   hlClose(95, 812);}
813   return M;
814 hlClose(94, 814);}
815
816 static GEN
817 _RgM_add(void *E, GEN x, GEN y) hlOpen(817,1);{ (void)E; return RgM_add(x, y); hlClose(96, 817);}
818
819 static GEN
820 _RgM_sub(void *E, GEN x, GEN y) hlOpen(820,1);{ (void)E; return RgM_sub(x, y); hlClose(97, 820);}
821
822 static GEN
823 _RgM_cmul(void *E, GEN P, long a, GEN x) hlOpen(823,1);{ (void)E; return RgM_Rg_mul(x,gel(P,a+2)); hlClose(98, 823);}
824
825 static GEN
826 _RgM_sqr(void *E, GEN x) hlOpen(826,1);{ (void) E; return RgM_sqr(x); hlClose(99, 826);}
827
828 static GEN
829 _RgM_mul(void *E, GEN x, GEN y) hlOpen(829,1);{ (void) E; return RgM_mul(x, y); hlClose(100, 829);}
830
831 static GEN
832 _RgM_one(void *E) hlOpen(832,1);{ long *n = (long*) E; return matid(*n); hlClose(101, 832);}
833
834 static GEN
835 _RgM_zero(void *E) hlOpen(835,1);{ long *n = (long*) E; return zeromat(*n,*n); hlClose(102, 835);}
836
837 static GEN
838 _RgM_red(void *E, GEN x) hlOpen(838,1);{ (void)E; return x; hlClose(103, 838);}
839
840 static struct bb_algebra RgM_algebra = hlOpen(840,1);{ _RgM_red, _RgM_add, _RgM_sub,
841        _RgM_mul, _RgM_sqr, _RgM_one, _RgM_zero hlClose(104, 841);};
842
843 /* generates the list of powers of x of degree 0,1,2,...,l*/
844 GEN
845 RgM_powers(GEN x, long l)
846 hlOpen(846,1);{
847   long n = lg(x)-1;
848   return gen_powers(x,l,1,(void *) &n, &_RgM_sqr, &_RgM_mul, &_RgM_one);
849 hlClose(105, 849);}
850
851 GEN
852 RgX_RgMV_eval(GEN Q, GEN x)
853 hlOpen(853,1);{
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 hlClose(106, 856);}
857
858 GEN
859 RgX_RgM_eval(GEN Q, GEN x)
860 hlOpen(860,1);{
861   long n = lg(x)-1;
862   return gen_bkeval(Q,degpol(Q),x,1,(void*)&n,&RgM_algebra,&_RgM_cmul);
863 hlClose(107, 863);}
864
865 GEN
866 RgC_Rg_div(GEN x, GEN y)
867 hlOpen(867,1);{ pari_APPLY_type(t_COL, gdiv(gel(x,i),y)) hlClose(108, 867);}
868
869 GEN
870 RgC_Rg_mul(GEN x, GEN y)
871 hlOpen(871,1);{ pari_APPLY_type(t_COL, gmul(gel(x,i),y)) hlClose(109, 871);}
872
873 GEN
874 RgV_Rg_mul(GEN x, GEN y)
875 hlOpen(875,1);{ pari_APPLY_type(t_VEC, gmul(gel(x,i),y)) hlClose(110, 875);}
876
877 GEN
878 RgM_Rg_div(GEN X, GEN c) hlOpen(878,1);{
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   hlOpen(884,2);{
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   hlClose(112, 888);}
889   return A;
890 hlClose(111, 890);}
891 GEN
892 RgM_Rg_mul(GEN X, GEN c) hlOpen(892,1);{
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   hlOpen(898,2);{
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   hlClose(114, 902);}
903   return A;
904 hlClose(113, 904);}
905
906 /********************************************************************/
907 /*                                                                  */
908 /*                    SCALAR TO MATRIX/VECTOR                       */
909 /*                                                                  */
910 /********************************************************************/
911 /* fill the square nxn matrix equal to t*Id */
912 static void
913 fill_scalmat(GEN y, GEN t, long n)
914 hlOpen(914,1);{
915   long i;
916   for (i = 1; i <= n; i++)
917   hlOpen(917,2);{
918     gel(y,i) = zerocol(n);
919     gcoeff(y,i,i) = t;
920   hlClose(116, 920);}
921 hlClose(115, 921);}
922
923 GEN
924 scalarmat(GEN x, long n) hlOpen(924,1);{
925   GEN y = cgetg(n+1, t_MAT);
926   if (!n) return y;
927   fill_scalmat(y, gcopy(x), n); return y;
928 hlClose(117, 928);}
929 GEN
930 scalarmat_shallow(GEN x, long n) hlOpen(930,1);{
931   GEN y = cgetg(n+1, t_MAT);
932   fill_scalmat(y, x, n); return y;
933 hlClose(118, 933);}
934 GEN
935 scalarmat_s(long x, long n) hlOpen(935,1);{
936   GEN y = cgetg(n+1, t_MAT);
937   if (!n) return y;
938   fill_scalmat(y, stoi(x), n); return y;
939 hlClose(119, 939);}
940 GEN
941 matid(long n) hlOpen(941,1);{
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 hlClose(120, 946);}
947
948 INLINE GEN
949 scalarcol_i(GEN x, long n, long c)
950 hlOpen(950,1);{
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 hlClose(121, 957);}
958
959 GEN
960 scalarcol(GEN x, long n) hlOpen(960,1);{ return scalarcol_i(x,n,1); hlClose(122, 960);}
961
962 GEN
963 scalarcol_shallow(GEN x, long n) hlOpen(963,1);{ return scalarcol_i(x,n,0); hlClose(123, 963);}
964
965 int
966 RgM_isscalar(GEN x, GEN s)
967 hlOpen(967,1);{
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   hlOpen(975,2);{
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   hlClose(125, 983);}
984   return 1;
985 hlClose(124, 985);}
986
987 int
988 RgM_isidentity(GEN x)
989 hlOpen(989,1);{
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   hlOpen(995,2);{
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   hlClose(127, 1003);}
1004   return 1;
1005 hlClose(126, 1005);}
1006
1007 long
1008 RgC_is_ei(GEN x)
1009 hlOpen(1009,1);{
1010   long i, j = 0, l = lg(x);
1011   for (i = 1; i < l; i++)
1012   hlOpen(1012,2);{
1013     GEN c = gel(x,i);
1014     if (gequal0(c)) continue;
1015     if (!gequal1(c) || j) return 0;
1016     j = i;
1017   hlClose(129, 1017);}
1018   return j;
1019 hlClose(128, 1019);}
1020
1021 int
1022 RgM_isdiagonal(GEN x)
1023 hlOpen(1023,1);{
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   hlOpen(1029,2);{
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   hlClose(131, 1035);}
1036   return 1;
1037 hlClose(130, 1037);}
1038 int
1039 isdiagonal(GEN x) hlOpen(1039,1);{ return (typ(x)==t_MAT) && RgM_isdiagonal(x); hlClose(132, 1039);}
1040
1041 GEN
1042 RgM_det_triangular(GEN mat)
1043 hlOpen(1043,1);{
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 hlClose(133, 1052);}
1053
1054 GEN
1055 RgV_kill0(GEN v)
1056 hlOpen(1056,1);{
1057   long i, l;
1058   GEN w = cgetg_copy(v, &l);
1059   for (i = 1; i < l; i++)
1060   hlOpen(1060,2);{
1061     GEN a = gel(v,i);
1062     gel(w,i) = gequal0(a) ? NULL: a;
1063   hlClose(135, 1063);}
1064   return w;
1065 hlClose(134, 1065);}
```