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