Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is to exceed 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - basemath - ZV.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.12.1 lcov report (development 25819-e703fe1174) Lines: 801 905 88.5 %
Date: 2020-09-18 06:10:04 Functions: 120 131 91.6 %
Legend: Lines: hit not hit

          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. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : #include "pari.h"
      15             : #include "paripriv.h"
      16             : 
      17             : static int
      18     1225446 : check_ZV(GEN x, long l)
      19             : {
      20             :   long i;
      21     7939049 :   for (i=1; i<l; i++)
      22     6713694 :     if (typ(gel(x,i)) != t_INT) return 0;
      23     1225355 :   return 1;
      24             : }
      25             : void
      26      516876 : RgV_check_ZV(GEN A, const char *s)
      27             : {
      28      516876 :   if (!RgV_is_ZV(A)) pari_err_TYPE(stack_strcat(s," [integer vector]"), A);
      29      516869 : }
      30             : void
      31      348285 : RgM_check_ZM(GEN A, const char *s)
      32             : {
      33      348285 :   long n = lg(A);
      34      348285 :   if (n != 1)
      35             :   {
      36      348229 :     long j, m = lgcols(A);
      37     1573584 :     for (j=1; j<n; j++)
      38     1225446 :       if (!check_ZV(gel(A,j), m))
      39          91 :         pari_err_TYPE(stack_strcat(s," [integer matrix]"), A);
      40             :   }
      41      348194 : }
      42             : 
      43             : static long
      44    65408111 : ZV_max_lg_i(GEN x, long m)
      45             : {
      46    65408111 :   long i, prec = 2;
      47   483838836 :   for (i=1; i<m; i++) { long l = lgefint(gel(x,i)); if (l > prec) prec = l; }
      48    65408111 :   return prec;
      49             : }
      50             : 
      51             : long
      52       10612 : ZV_max_lg(GEN x)
      53       10612 : { return ZV_max_lg_i(x, lg(x)); }
      54             : 
      55             : static long
      56    22486712 : ZM_max_lg_i(GEN x, long n, long m)
      57             : {
      58    22486712 :   long prec = 2;
      59    22486712 :   if (n != 1)
      60             :   {
      61             :     long j;
      62    87884211 :     for (j=1; j<n; j++)
      63             :     {
      64    65397499 :       long l = ZV_max_lg_i(gel(x,j), m);
      65    65397499 :       if (l > prec) prec = l;
      66             :     }
      67             :   }
      68    22486712 :   return prec;
      69             : }
      70             : 
      71             : long
      72      424458 : ZM_max_lg(GEN x)
      73             : {
      74      424458 :   long n = lg(x);
      75      424458 :   if (n==1) return 2;
      76      424458 :   return ZM_max_lg_i(x, n, lgcols(x));
      77             : }
      78             : 
      79             : GEN
      80        3394 : ZM_supnorm(GEN x)
      81             : {
      82        3394 :   long i, j, h, lx = lg(x);
      83        3394 :   GEN s = gen_0;
      84        3394 :   if (lx == 1) return gen_1;
      85        3394 :   h = lgcols(x);
      86       20533 :   for (j=1; j<lx; j++)
      87             :   {
      88       17139 :     GEN xj = gel(x,j);
      89      240192 :     for (i=1; i<h; i++)
      90             :     {
      91      223053 :       GEN c = gel(xj,i);
      92      223053 :       if (abscmpii(c, s) > 0) s = c;
      93             :     }
      94             :   }
      95        3394 :   return absi(s);
      96             : }
      97             : 
      98             : /********************************************************************/
      99             : /**                                                                **/
     100             : /**                           MULTIPLICATION                       **/
     101             : /**                                                                **/
     102             : /********************************************************************/
     103             : /* x non-empty ZM, y a compatible nc (dimension > 0). */
     104             : static GEN
     105           0 : ZM_nc_mul_i(GEN x, GEN y, long c, long l)
     106             : {
     107             :   long i, j;
     108             :   pari_sp av;
     109           0 :   GEN z = cgetg(l,t_COL), s;
     110             : 
     111           0 :   for (i=1; i<l; i++)
     112             :   {
     113           0 :     av = avma; s = muliu(gcoeff(x,i,1),y[1]);
     114           0 :     for (j=2; j<c; j++)
     115           0 :       if (y[j]) s = addii(s, muliu(gcoeff(x,i,j),y[j]));
     116           0 :     gel(z,i) = gerepileuptoint(av,s);
     117             :   }
     118           0 :   return z;
     119             : }
     120             : 
     121             : /* x ZV, y a compatible zc. */
     122             : GEN
     123     2214660 : ZV_zc_mul(GEN x, GEN y)
     124             : {
     125     2214660 :   long j, l = lg(x);
     126     2214660 :   pari_sp av = avma;
     127     2214660 :   GEN s = mulis(gel(x,1),y[1]);
     128    50028650 :   for (j=2; j<l; j++)
     129    47813990 :     if (y[j]) s = addii(s, mulis(gel(x,j),y[j]));
     130     2214660 :   return gerepileuptoint(av,s);
     131             : }
     132             : 
     133             : /* x non-empty ZM, y a compatible zc (dimension > 0). */
     134             : static GEN
     135     6065742 : ZM_zc_mul_i(GEN x, GEN y, long c, long l)
     136             : {
     137             :   long i, j;
     138     6065742 :   GEN z = cgetg(l,t_COL);
     139             : 
     140    38033224 :   for (i=1; i<l; i++)
     141             :   {
     142    31967482 :     pari_sp av = avma;
     143    31967482 :     GEN s = mulis(gcoeff(x,i,1),y[1]);
     144   608883446 :     for (j=2; j<c; j++)
     145   576915964 :       if (y[j]) s = addii(s, mulis(gcoeff(x,i,j),y[j]));
     146    31967482 :     gel(z,i) = gerepileuptoint(av,s);
     147             :   }
     148     6065742 :   return z;
     149             : }
     150             : GEN
     151     4496546 : ZM_zc_mul(GEN x, GEN y) {
     152     4496546 :   long l = lg(x);
     153     4496546 :   if (l == 1) return cgetg(1, t_COL);
     154     4496546 :   return ZM_zc_mul_i(x,y, l, lgcols(x));
     155             : }
     156             : 
     157             : /* y non-empty ZM, x a compatible zv (dimension > 0). */
     158             : GEN
     159        1372 : zv_ZM_mul(GEN x, GEN y) {
     160        1372 :   long i,j, lx = lg(x), ly = lg(y);
     161             :   GEN z;
     162        1372 :   if (lx == 1) return zerovec(ly-1);
     163        1372 :   z = cgetg(ly,t_VEC);
     164        3304 :   for (j=1; j<ly; j++)
     165             :   {
     166        1932 :     pari_sp av = avma;
     167        1932 :     GEN s = mulsi(x[1], gcoeff(y,1,j));
     168        3584 :     for (i=2; i<lx; i++)
     169        1652 :       if (x[i]) s = addii(s, mulsi(x[i], gcoeff(y,i,j)));
     170        1932 :     gel(z,j) = gerepileuptoint(av,s);
     171             :   }
     172        1372 :   return z;
     173             : }
     174             : 
     175             : /* x ZM, y a compatible zm (dimension > 0). */
     176             : GEN
     177      759617 : ZM_zm_mul(GEN x, GEN y)
     178             : {
     179      759617 :   long j, c, l = lg(x), ly = lg(y);
     180      759617 :   GEN z = cgetg(ly, t_MAT);
     181      759617 :   if (l == 1) return z;
     182      759617 :   c = lgcols(x);
     183     2328813 :   for (j = 1; j < ly; j++) gel(z,j) = ZM_zc_mul_i(x, gel(y,j), l,c);
     184      759617 :   return z;
     185             : }
     186             : /* x ZM, y a compatible zn (dimension > 0). */
     187             : GEN
     188           0 : ZM_nm_mul(GEN x, GEN y)
     189             : {
     190           0 :   long j, c, l = lg(x), ly = lg(y);
     191           0 :   GEN z = cgetg(ly, t_MAT);
     192           0 :   if (l == 1) return z;
     193           0 :   c = lgcols(x);
     194           0 :   for (j = 1; j < ly; j++) gel(z,j) = ZM_nc_mul_i(x, gel(y,j), l,c);
     195           0 :   return z;
     196             : }
     197             : 
     198             : /* Strassen-Winograd algorithm */
     199             : 
     200             : /*
     201             :   Return A[ma+1..ma+da, na+1..na+ea] - B[mb+1..mb+db, nb+1..nb+eb]
     202             :   as an (m x n)-matrix, padding the input with zeroes as necessary.
     203             : */
     204             : static GEN
     205      407488 : add_slices(long m, long n,
     206             :            GEN A, long ma, long da, long na, long ea,
     207             :            GEN B, long mb, long db, long nb, long eb)
     208             : {
     209      407488 :   long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
     210      407488 :   GEN M = cgetg(n + 1, t_MAT), C;
     211             : 
     212     2379749 :   for (j = 1; j <= min_e; j++) {
     213     1972261 :     gel(M, j) = C = cgetg(m + 1, t_COL);
     214    28200578 :     for (i = 1; i <= min_d; i++)
     215    26228317 :       gel(C, i) = addii(gcoeff(A, ma + i, na + j),
     216    26228317 :                         gcoeff(B, mb + i, nb + j));
     217     2004129 :     for (; i <= da; i++)
     218       31868 :       gel(C, i) = gcoeff(A, ma + i, na + j);
     219     1972261 :     for (; i <= db; i++)
     220           0 :       gel(C, i) = gcoeff(B, mb + i, nb + j);
     221     1972261 :     for (; i <= m; i++)
     222           0 :       gel(C, i) = gen_0;
     223             :   }
     224      435047 :   for (; j <= ea; j++) {
     225       27559 :     gel(M, j) = C = cgetg(m + 1, t_COL);
     226       96612 :     for (i = 1; i <= da; i++)
     227       69053 :       gel(C, i) = gcoeff(A, ma + i, na + j);
     228       27559 :     for (; i <= m; i++)
     229           0 :       gel(C, i) = gen_0;
     230             :   }
     231      407488 :   for (; j <= eb; j++) {
     232           0 :     gel(M, j) = C = cgetg(m + 1, t_COL);
     233           0 :     for (i = 1; i <= db; i++)
     234           0 :       gel(C, i) = gcoeff(B, mb + i, nb + j);
     235           0 :     for (; i <= m; i++)
     236           0 :       gel(C, i) = gen_0;
     237             :   }
     238      407488 :   for (; j <= n; j++)
     239           0 :     gel(M, j) = zerocol(m);
     240      407488 :   return M;
     241             : }
     242             : 
     243             : /*
     244             :   Return A[ma+1..ma+da, na+1..na+ea] - B[mb+1..mb+db, nb+1..nb+eb]
     245             :   as an (m x n)-matrix, padding the input with zeroes as necessary.
     246             : */
     247             : static GEN
     248      356552 : subtract_slices(long m, long n,
     249             :                 GEN A, long ma, long da, long na, long ea,
     250             :                 GEN B, long mb, long db, long nb, long eb)
     251             : {
     252      356552 :   long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
     253      356552 :   GEN M = cgetg(n + 1, t_MAT), C;
     254             : 
     255     2077894 :   for (j = 1; j <= min_e; j++) {
     256     1721342 :     gel(M, j) = C = cgetg(m + 1, t_COL);
     257    25766481 :     for (i = 1; i <= min_d; i++)
     258    24045139 :       gel(C, i) = subii(gcoeff(A, ma + i, na + j),
     259    24045139 :                         gcoeff(B, mb + i, nb + j));
     260     1744315 :     for (; i <= da; i++)
     261       22973 :       gel(C, i) = gcoeff(A, ma + i, na + j);
     262     1758187 :     for (; i <= db; i++)
     263       36845 :       gel(C, i) = negi(gcoeff(B, mb + i, nb + j));
     264     1721342 :     for (; i <= m; i++)
     265           0 :       gel(C, i) = gen_0;
     266             :   }
     267      356552 :   for (; j <= ea; j++) {
     268           0 :     gel(M, j) = C = cgetg(m + 1, t_COL);
     269           0 :     for (i = 1; i <= da; i++)
     270           0 :       gel(C, i) = gcoeff(A, ma + i, na + j);
     271           0 :     for (; i <= m; i++)
     272           0 :       gel(C, i) = gen_0;
     273             :   }
     274      376792 :   for (; j <= eb; j++) {
     275       20240 :     gel(M, j) = C = cgetg(m + 1, t_COL);
     276       65873 :     for (i = 1; i <= db; i++)
     277       45633 :       gel(C, i) = negi(gcoeff(B, mb + i, nb + j));
     278       20240 :     for (; i <= m; i++)
     279           0 :       gel(C, i) = gen_0;
     280             :   }
     281      376792 :   for (; j <= n; j++)
     282       20240 :     gel(M, j) = zerocol(m);
     283      356552 :   return M;
     284             : }
     285             : 
     286             : static GEN ZM_mul_i(GEN x, GEN y, long l, long lx, long ly);
     287             : 
     288             : /* Strassen-Winograd matrix product A (m x n) * B (n x p) */
     289             : static GEN
     290       50936 : ZM_mul_sw(GEN A, GEN B, long m, long n, long p)
     291             : {
     292       50936 :   pari_sp av = avma;
     293       50936 :   long m1 = (m + 1)/2, m2 = m/2,
     294       50936 :     n1 = (n + 1)/2, n2 = n/2,
     295       50936 :     p1 = (p + 1)/2, p2 = p/2;
     296             :   GEN A11, A12, A22, B11, B21, B22,
     297             :     S1, S2, S3, S4, T1, T2, T3, T4,
     298             :     M1, M2, M3, M4, M5, M6, M7,
     299             :     V1, V2, V3, C11, C12, C21, C22, C;
     300             : 
     301       50936 :   T2 = subtract_slices(n1, p2, B, 0, n1, p1, p2, B, n1, n2, p1, p2);
     302       50936 :   S1 = subtract_slices(m2, n1, A, m1, m2, 0, n1, A, 0, m2, 0, n1);
     303       50936 :   M2 = ZM_mul_i(S1, T2, m2 + 1, n1 + 1, p2 + 1);
     304       50936 :   if (gc_needed(av, 1))
     305           0 :     gerepileall(av, 2, &T2, &M2);  /* destroy S1 */
     306       50936 :   T3 = subtract_slices(n1, p1, T2, 0, n1, 0, p2, B, 0, n1, 0, p1);
     307       50936 :   if (gc_needed(av, 1))
     308           0 :     gerepileall(av, 2, &M2, &T3);  /* destroy T2 */
     309       50936 :   S2 = add_slices(m2, n1, A, m1, m2, 0, n1, A, m1, m2, n1, n2);
     310       50936 :   T1 = subtract_slices(n1, p1, B, 0, n1, p1, p2, B, 0, n1, 0, p2);
     311       50936 :   M3 = ZM_mul_i(S2, T1, m2 + 1, n1 + 1, p2 + 1);
     312       50936 :   if (gc_needed(av, 1))
     313           0 :     gerepileall(av, 4, &M2, &T3, &S2, &M3);  /* destroy T1 */
     314       50936 :   S3 = subtract_slices(m1, n1, S2, 0, m2, 0, n1, A, 0, m1, 0, n1);
     315       50936 :   if (gc_needed(av, 1))
     316           0 :     gerepileall(av, 4, &M2, &T3, &M3, &S3);  /* destroy S2 */
     317       50936 :   A11 = matslice(A, 1, m1, 1, n1);
     318       50936 :   B11 = matslice(B, 1, n1, 1, p1);
     319       50936 :   M1 = ZM_mul_i(A11, B11, m1 + 1, n1 + 1, p1 + 1);
     320       50936 :   if (gc_needed(av, 1))
     321           0 :     gerepileall(av, 5, &M2, &T3, &M3, &S3, &M1);  /* destroy A11, B11 */
     322       50936 :   A12 = matslice(A, 1, m1, n1 + 1, n);
     323       50936 :   B21 = matslice(B, n1 + 1, n, 1, p1);
     324       50936 :   M4 = ZM_mul_i(A12, B21, m1 + 1, n2 + 1, p1 + 1);
     325       50936 :   if (gc_needed(av, 1))
     326           0 :     gerepileall(av, 6, &M2, &T3, &M3, &S3, &M1, &M4);  /* destroy A12, B21 */
     327       50936 :   C11 = add_slices(m1, p1, M1, 0, m1, 0, p1, M4, 0, m1, 0, p1);
     328       50936 :   if (gc_needed(av, 1))
     329           0 :     gerepileall(av, 6, &M2, &T3, &M3, &S3, &M1, &C11);  /* destroy M4 */
     330       50936 :   M5 = ZM_mul_i(S3, T3, m1 + 1, n1 + 1, p1 + 1);
     331       50936 :   S4 = subtract_slices(m1, n2, A, 0, m1, n1, n2, S3, 0, m1, 0, n2);
     332       50936 :   if (gc_needed(av, 1))
     333           5 :     gerepileall(av, 7, &M2, &T3, &M3, &M1, &C11, &M5, &S4);  /* destroy S3 */
     334       50936 :   T4 = add_slices(n2, p1, B, n1, n2, 0, p1, T3, 0, n2, 0, p1);
     335       50936 :   if (gc_needed(av, 1))
     336           0 :     gerepileall(av, 7, &M2, &M3, &M1, &C11, &M5, &S4, &T4);  /* destroy T3 */
     337       50936 :   V1 = subtract_slices(m1, p1, M1, 0, m1, 0, p1, M5, 0, m1, 0, p1);
     338       50936 :   if (gc_needed(av, 1))
     339           1 :     gerepileall(av, 6, &M2, &M3, &S4, &T4, &C11, &V1);  /* destroy M1, M5 */
     340       50936 :   B22 = matslice(B, n1 + 1, n, p1 + 1, p);
     341       50936 :   M6 = ZM_mul_i(S4, B22, m1 + 1, n2 + 1, p2 + 1);
     342       50936 :   if (gc_needed(av, 1))
     343           6 :     gerepileall(av, 6, &M2, &M3, &T4, &C11, &V1, &M6);  /* destroy S4, B22 */
     344       50936 :   A22 = matslice(A, m1 + 1, m, n1 + 1, n);
     345       50936 :   M7 = ZM_mul_i(A22, T4, m2 + 1, n2 + 1, p1 + 1);
     346       50936 :   if (gc_needed(av, 1))
     347           0 :     gerepileall(av, 6, &M2, &M3, &C11, &V1, &M6, &M7);  /* destroy A22, T4 */
     348       50936 :   V3 = add_slices(m1, p2, V1, 0, m1, 0, p2, M3, 0, m2, 0, p2);
     349       50936 :   C12 = add_slices(m1, p2, V3, 0, m1, 0, p2, M6, 0, m1, 0, p2);
     350       50936 :   if (gc_needed(av, 1))
     351           6 :     gerepileall(av, 6, &M2, &M3, &C11, &V1, &M7, &C12);  /* destroy V3, M6 */
     352       50936 :   V2 = add_slices(m2, p1, V1, 0, m2, 0, p1, M2, 0, m2, 0, p2);
     353       50936 :   if (gc_needed(av, 1))
     354           0 :     gerepileall(av, 5, &M3, &C11, &M7, &C12, &V2);  /* destroy V1, M2 */
     355       50936 :   C21 = add_slices(m2, p1, V2, 0, m2, 0, p1, M7, 0, m2, 0, p1);
     356       50936 :   if (gc_needed(av, 1))
     357           7 :     gerepileall(av, 5, &M3, &C11, &C12, &V2, &C21);  /* destroy M7 */
     358       50936 :   C22 = add_slices(m2, p2, V2, 0, m2, 0, p2, M3, 0, m2, 0, p2);
     359       50936 :   if (gc_needed(av, 1))
     360           0 :     gerepileall(av, 4, &C11, &C12, &C21, &C22);  /* destroy V2, M3 */
     361       50936 :   C = shallowconcat(vconcat(C11, C21), vconcat(C12, C22));
     362       50936 :   return gerepilecopy(av, C);
     363             : }
     364             : 
     365             : /* x[i,]*y. Assume lg(x) > 1 and 0 < i < lgcols(x) */
     366             : static GEN
     367   328645228 : ZMrow_ZC_mul_i(GEN x, GEN y, long i, long lx)
     368             : {
     369   328645228 :   pari_sp av = avma;
     370   328645228 :   GEN c = mulii(gcoeff(x,i,1), gel(y,1)), ZERO = gen_0;
     371             :   long k;
     372  3381809748 :   for (k = 2; k < lx; k++)
     373             :   {
     374  3053164906 :     GEN t = mulii(gcoeff(x,i,k), gel(y,k));
     375  3053164656 :     if (t != ZERO) c = addii(c, t);
     376             :   }
     377   328644842 :   return gerepileuptoint(av, c);
     378             : }
     379             : GEN
     380    82024655 : ZMrow_ZC_mul(GEN x, GEN y, long i)
     381    82024655 : { return ZMrow_ZC_mul_i(x, y, i, lg(x)); }
     382             : 
     383             : /* return x * y, 1 < lx = lg(x), l = lgcols(x) */
     384             : static GEN
     385    47196133 : ZM_ZC_mul_i(GEN x, GEN y, long lx, long l)
     386             : {
     387    47196133 :   GEN z = cgetg(l,t_COL);
     388             :   long i;
     389   293816676 :   for (i=1; i<l; i++) gel(z,i) = ZMrow_ZC_mul_i(x,y,i,lx);
     390    47196102 :   return z;
     391             : }
     392             : 
     393             : static GEN
     394    10982851 : ZM_mul_classical(GEN x, GEN y, long l, long lx, long ly)
     395             : {
     396             :   long j;
     397    10982851 :   GEN z = cgetg(ly, t_MAT);
     398    41481210 :   for (j = 1; j < ly; j++)
     399    30498359 :     gel(z, j) = ZM_ZC_mul_i(x, gel(y, j), lx, l);
     400    10982851 :   return z;
     401             : }
     402             : 
     403             : static GEN
     404        1372 : ZM_mul_slice(GEN A, GEN B, GEN P, GEN *mod)
     405             : {
     406        1372 :   pari_sp av = avma;
     407        1372 :   long i, n = lg(P)-1;
     408             :   GEN H, T;
     409        1372 :   if (n == 1)
     410             :   {
     411         742 :     ulong p = uel(P,1);
     412         742 :     GEN a = ZM_to_Flm(A, p);
     413         742 :     GEN b = ZM_to_Flm(B, p);
     414         742 :     GEN Hp = gerepileupto(av, Flm_to_ZM(Flm_mul(a, b, p)));
     415         742 :     *mod = utoi(p); return Hp;
     416             :   }
     417         630 :   T = ZV_producttree(P);
     418         630 :   A = ZM_nv_mod_tree(A, P, T);
     419         630 :   B = ZM_nv_mod_tree(B, P, T);
     420         630 :   H = cgetg(n+1, t_VEC);
     421       10426 :   for(i=1; i <= n; i++)
     422        9796 :     gel(H,i) = Flm_mul(gel(A,i),gel(B,i),P[i]);
     423         630 :   H = nmV_chinese_center_tree_seq(H, P, T, ZV_chinesetree(P, T));
     424         630 :   *mod = gmael(T, lg(T)-1, 1);
     425         630 :   gerepileall(av, 2, &H, mod);
     426         630 :   return H;
     427             : }
     428             : 
     429             : GEN
     430        1372 : ZM_mul_worker(GEN P, GEN A, GEN B)
     431             : {
     432        1372 :   GEN V = cgetg(3, t_VEC);
     433        1372 :   gel(V,1) = ZM_mul_slice(A, B, P, &gel(V,2));
     434        1372 :   return V;
     435             : }
     436             : 
     437             : static long
     438      167497 : ZC_expi(GEN x)
     439             : {
     440      167497 :   long i, l = lg(x), m=-1;
     441    17089684 :   for(i = 1; i < l; i++)
     442             :   {
     443    16922187 :     GEN xi = gel(x,i);
     444    16922187 :     long e = signe(xi) ? expi(xi): -1;
     445    16922187 :     if (e > m) m = e;
     446             :   }
     447      167497 :   return m;
     448             : }
     449             : 
     450             : static long
     451        1740 : ZM_expi(GEN x)
     452             : {
     453        1740 :   long i, l = lg(x), m=-1;
     454      169237 :   for(i = 1; i < l; i++)
     455             :   {
     456      167497 :     long e = ZC_expi(gel(x,i));
     457      167497 :     if (e > m) m = e;
     458             :   }
     459        1740 :   return m;
     460             : }
     461             : 
     462             : static GEN
     463         870 : ZM_mul_fast(GEN A, GEN B, long lA, long lB)
     464             : {
     465         870 :   pari_sp av = avma;
     466             :   forprime_t S;
     467             :   GEN H, worker;
     468         870 :   long h, mA = ZM_expi(A), mB = ZM_expi(B);
     469         870 :   if (mA==-1 || mB==-1) return zeromat(nbrows(A),lB-1);
     470         856 :   h = 3 + mA + mB + expu(lA-1);
     471         856 :   init_modular_big(&S);
     472         856 :   worker = snm_closure(is_entry("_ZM_mul_worker"), mkvec2(A,B));
     473         856 :   H = gen_crt("ZM_mul", worker, &S, NULL, h, 0, NULL,
     474             :               nmV_chinese_center, FpM_center);
     475         856 :   return gerepileupto(av, H);
     476             : }
     477             : 
     478             : /* Strassen-Winograd used for dim >= ZM_sw_bound */
     479             : static GEN
     480    11029337 : ZM_mul_i(GEN x, GEN y, long l, long lx, long ly)
     481             : {
     482    11029337 :   if (lx > 70 && ly > 70 && l > 70)
     483         870 :     return ZM_mul_fast(x, y, lx, ly);
     484             :   else
     485             :   {
     486    11028467 :     long s = maxss(ZM_max_lg_i(x,lx,l), ZM_max_lg_i(y,ly,lx));
     487    11028467 :     long ZM_sw_bound = s > 60 ? 2: s > 25 ? 4: s>15 ? 8 : s > 8 ? 16 : 32;
     488    11028467 :     if (l <= ZM_sw_bound || lx <= ZM_sw_bound || ly <= ZM_sw_bound)
     489    10977540 :       return ZM_mul_classical(x, y, l, lx, ly);
     490             :     else
     491       50927 :       return ZM_mul_sw(x, y, l - 1, lx - 1, ly - 1);
     492             :   }
     493             : }
     494             : 
     495             : GEN
     496    10694115 : ZM_mul(GEN x, GEN y)
     497             : {
     498    10694115 :   long lx=lg(x), ly=lg(y);
     499    10694115 :   if (ly==1) return cgetg(1,t_MAT);
     500    10673961 :   if (lx==1) return zeromat(0, ly-1);
     501    10672785 :   return ZM_mul_i(x, y, lgcols(x), lx, ly);
     502             : }
     503             : 
     504             : GEN
     505      139165 : QM_mul(GEN x, GEN y)
     506             : {
     507      139165 :   GEN dx, nx = Q_primitive_part(x, &dx);
     508      139165 :   GEN dy, ny = Q_primitive_part(y, &dy);
     509      139165 :   GEN z = ZM_mul(nx, ny);
     510      139165 :   if (dx || dy)
     511             :   {
     512      139165 :     GEN d = dx ? dy ? gmul(dx, dy): dx : dy;
     513      139165 :     if (!gequal1(d)) z = ZM_Q_mul(z, d);
     514             :   }
     515      139165 :   return z;
     516             : }
     517             : 
     518             : GEN
     519         700 : QM_sqr(GEN x)
     520             : {
     521         700 :   GEN dx, nx = Q_primitive_part(x, &dx);
     522         700 :   GEN z = ZM_sqr(nx);
     523         700 :   if (dx)
     524         700 :     z = ZM_Q_mul(z, gsqr(dx));
     525         700 :   return z;
     526             : }
     527             : 
     528             : GEN
     529      116226 : QM_QC_mul(GEN x, GEN y)
     530             : {
     531      116226 :   GEN dx, nx = Q_primitive_part(x, &dx);
     532      116226 :   GEN dy, ny = Q_primitive_part(y, &dy);
     533      116226 :   GEN z = ZM_ZC_mul(nx, ny);
     534      116226 :   if (dx || dy)
     535             :   {
     536      116198 :     GEN d = dx ? dy ? gmul(dx, dy): dx : dy;
     537      116198 :     if (!gequal1(d)) z = ZC_Q_mul(z, d);
     538             :   }
     539      116226 :   return z;
     540             : }
     541             : 
     542             : /* assume result is symmetric */
     543             : GEN
     544           0 : ZM_multosym(GEN x, GEN y)
     545             : {
     546           0 :   long j, lx, ly = lg(y);
     547             :   GEN M;
     548           0 :   if (ly == 1) return cgetg(1,t_MAT);
     549           0 :   lx = lg(x); /* = lgcols(y) */
     550           0 :   if (lx == 1) return cgetg(1,t_MAT);
     551             :   /* ly = lgcols(x) */
     552           0 :   M = cgetg(ly, t_MAT);
     553           0 :   for (j=1; j<ly; j++)
     554             :   {
     555           0 :     GEN z = cgetg(ly,t_COL), yj = gel(y,j);
     556             :     long i;
     557           0 :     for (i=1; i<j; i++) gel(z,i) = gcoeff(M,j,i);
     558           0 :     for (i=j; i<ly; i++)gel(z,i) = ZMrow_ZC_mul_i(x,yj,i,lx);
     559           0 :     gel(M,j) = z;
     560             :   }
     561           0 :   return M;
     562             : }
     563             : 
     564             : /* compute m*diagonal(d), assume lg(d) = lg(m). Shallow */
     565             : GEN
     566          14 : ZM_mul_diag(GEN m, GEN d)
     567             : {
     568             :   long j, l;
     569          14 :   GEN y = cgetg_copy(m, &l);
     570          56 :   for (j=1; j<l; j++)
     571             :   {
     572          42 :     GEN c = gel(d,j);
     573          42 :     gel(y,j) = equali1(c)? gel(m,j): ZC_Z_mul(gel(m,j), c);
     574             :   }
     575          14 :   return y;
     576             : }
     577             : /* compute diagonal(d)*m, assume lg(d) = lg(m~). Shallow */
     578             : GEN
     579      132902 : ZM_diag_mul(GEN d, GEN m)
     580             : {
     581      132902 :   long i, j, l = lg(d), lm = lg(m);
     582      132902 :   GEN y = cgetg(lm, t_MAT);
     583      402963 :   for (j=1; j<lm; j++) gel(y,j) = cgetg(l, t_COL);
     584      420883 :   for (i=1; i<l; i++)
     585             :   {
     586      287981 :     GEN c = gel(d,i);
     587      287981 :     if (equali1(c))
     588       33866 :       for (j=1; j<lm; j++) gcoeff(y,i,j) = gcoeff(m,i,j);
     589             :     else
     590      968902 :       for (j=1; j<lm; j++) gcoeff(y,i,j) = mulii(gcoeff(m,i,j), c);
     591             :   }
     592      132902 :   return y;
     593             : }
     594             : 
     595             : /* assume lx > 1 is lg(x) = lg(y) */
     596             : static GEN
     597    16634665 : ZV_dotproduct_i(GEN x, GEN y, long lx)
     598             : {
     599    16634665 :   pari_sp av = avma;
     600    16634665 :   GEN c = mulii(gel(x,1), gel(y,1));
     601             :   long i;
     602   131600389 :   for (i = 2; i < lx; i++)
     603             :   {
     604   114965722 :     GEN t = mulii(gel(x,i), gel(y,i));
     605   114965725 :     if (t != gen_0) c = addii(c, t);
     606             :   }
     607    16634667 :   return gerepileuptoint(av, c);
     608             : }
     609             : 
     610             : /* x~ * y, assuming result is symmetric */
     611             : GEN
     612      161775 : ZM_transmultosym(GEN x, GEN y)
     613             : {
     614      161775 :   long i, j, l, ly = lg(y);
     615             :   GEN M;
     616      161775 :   if (ly == 1) return cgetg(1,t_MAT);
     617             :   /* lg(x) = ly */
     618      161775 :   l = lgcols(y); /* = lgcols(x) */
     619      161775 :   M = cgetg(ly, t_MAT);
     620     1266443 :   for (i=1; i<ly; i++)
     621             :   {
     622     1104668 :     GEN xi = gel(x,i), c = cgetg(ly,t_COL);
     623     1104668 :     gel(M,i) = c;
     624     5067071 :     for (j=1; j<i; j++)
     625     3962403 :       gcoeff(M,i,j) = gel(c,j) = ZV_dotproduct_i(xi,gel(y,j),l);
     626     1104668 :     gel(c,i) = ZV_dotproduct_i(xi,gel(y,i),l);
     627             :   }
     628      161775 :   return M;
     629             : }
     630             : /* x~ * y */
     631             : GEN
     632         497 : ZM_transmul(GEN x, GEN y)
     633             : {
     634         497 :   long i, j, l, lx, ly = lg(y);
     635             :   GEN M;
     636         497 :   if (ly == 1) return cgetg(1,t_MAT);
     637         497 :   lx = lg(x);
     638         497 :   l = lgcols(y);
     639         497 :   if (lgcols(x) != l) pari_err_OP("operation 'ZM_transmul'", x,y);
     640         497 :   M = cgetg(ly, t_MAT);
     641        1617 :   for (i=1; i<ly; i++)
     642             :   {
     643        1120 :     GEN yi = gel(y,i), c = cgetg(lx,t_COL);
     644        1120 :     gel(M,i) = c;
     645        5061 :     for (j=1; j<lx; j++) gel(c,j) = ZV_dotproduct_i(yi,gel(x,j),l);
     646             :   }
     647         497 :   return M;
     648             : }
     649             : 
     650             : static GEN
     651        5320 : ZM_sqr_i(GEN x, long l)
     652             : {
     653        5320 :   if (l > 70)
     654           0 :     return ZM_mul_fast(x, x, l, l);
     655             :   else
     656             :   {
     657        5320 :     long s = ZM_max_lg_i(x,l,l);
     658        5320 :     long ZM_sw_bound = s > 60 ? 2: s > 25 ? 4: s>15 ? 8 : s > 8 ? 16 : 32;
     659        5320 :     if (l <= ZM_sw_bound)
     660        5311 :       return ZM_mul_classical(x, x, l, l, l);
     661             :     else
     662           9 :       return ZM_mul_sw(x, x, l - 1, l - 1, l - 1);
     663             :   }
     664             : }
     665             : 
     666             : GEN
     667        5320 : ZM_sqr(GEN x)
     668             : {
     669        5320 :   long lx=lg(x);
     670        5320 :   if (lx==1) return cgetg(1,t_MAT);
     671        5320 :   return ZM_sqr_i(x, lx);
     672             : }
     673             : GEN
     674    16704283 : ZM_ZC_mul(GEN x, GEN y)
     675             : {
     676    16704283 :   long lx = lg(x);
     677    16704283 :   return lx==1? cgetg(1,t_COL): ZM_ZC_mul_i(x, y, lx, lgcols(x));
     678             : }
     679             : 
     680             : GEN
     681     1292495 : ZC_Z_div(GEN x, GEN c)
     682     6302380 : { pari_APPLY_type(t_COL, Qdivii(gel(x,i), c)) }
     683             : 
     684             : GEN
     685       14560 : ZM_Z_div(GEN x, GEN c)
     686      164724 : { pari_APPLY_same(ZC_Z_div(gel(x, i), c)) }
     687             : 
     688             : GEN
     689     1316168 : ZC_Q_mul(GEN A, GEN z)
     690             : {
     691     1316168 :   pari_sp av = avma;
     692     1316168 :   long i, l = lg(A);
     693             :   GEN d, n, Ad, B, u;
     694     1316168 :   if (typ(z)==t_INT) return ZC_Z_mul(A,z);
     695     1313683 :   n = gel(z, 1); d = gel(z, 2);
     696     1313683 :   Ad = FpC_red(A, d);
     697     1313683 :   u = gcdii(d, FpV_factorback(Ad, NULL, d));
     698     1313683 :   B = cgetg(l, t_COL);
     699     1313683 :   if (equali1(u))
     700             :   {
     701      249069 :     for(i=1; i<l; i++)
     702      205831 :       gel(B, i) = mkfrac(mulii(n, gel(A,i)), d);
     703             :   } else
     704             :   {
     705    12252798 :     for(i=1; i<l; i++)
     706             :     {
     707    10982353 :       GEN di = gcdii(gel(Ad, i), u), ni = mulii(n, diviiexact(gel(A,i), di));
     708    10982353 :       if (equalii(d, di))
     709     7878181 :         gel(B, i) = ni;
     710             :       else
     711     3104172 :         gel(B, i) = mkfrac(ni, diviiexact(d, di));
     712             :     }
     713             :   }
     714     1313683 :   return gerepilecopy(av, B);
     715             : }
     716             : 
     717             : GEN
     718      286863 : ZM_Q_mul(GEN x, GEN z)
     719             : {
     720      286863 :   if (typ(z)==t_INT) return ZM_Z_mul(x,z);
     721     1396317 :   pari_APPLY_same(ZC_Q_mul(gel(x, i), z));
     722             : }
     723             : 
     724             : long
     725   182226205 : zv_dotproduct(GEN x, GEN y)
     726             : {
     727   182226205 :   long i, lx = lg(x);
     728             :   ulong c;
     729   182226205 :   if (lx == 1) return 0;
     730   182226205 :   c = uel(x,1)*uel(y,1);
     731  2825755660 :   for (i = 2; i < lx; i++)
     732  2643529455 :     c += uel(x,i)*uel(y,i);
     733   182226205 :   return c;
     734             : }
     735             : 
     736             : GEN
     737      166087 : ZV_ZM_mul(GEN x, GEN y)
     738             : {
     739      166087 :   long i, lx = lg(x), ly = lg(y);
     740             :   GEN z;
     741      166087 :   if (lx == 1) return zerovec(ly-1);
     742      166073 :   z = cgetg(ly, t_VEC);
     743      625017 :   for (i = 1; i < ly; i++) gel(z,i) = ZV_dotproduct_i(x, gel(y,i), lx);
     744      166073 :   return z;
     745             : }
     746             : 
     747             : GEN
     748           0 : ZC_ZV_mul(GEN x, GEN y)
     749             : {
     750           0 :   long i,j, lx=lg(x), ly=lg(y);
     751             :   GEN z;
     752           0 :   if (ly==1) return cgetg(1,t_MAT);
     753           0 :   z = cgetg(ly,t_MAT);
     754           0 :   for (j=1; j < ly; j++)
     755             :   {
     756           0 :     gel(z,j) = cgetg(lx,t_COL);
     757           0 :     for (i=1; i<lx; i++) gcoeff(z,i,j) = mulii(gel(x,i),gel(y,j));
     758             :   }
     759           0 :   return z;
     760             : }
     761             : 
     762             : GEN
     763     7462024 : ZV_dotsquare(GEN x)
     764             : {
     765             :   long i, lx;
     766             :   pari_sp av;
     767             :   GEN z;
     768     7462024 :   lx = lg(x);
     769     7462024 :   if (lx == 1) return gen_0;
     770     7432880 :   av = avma; z = sqri(gel(x,1));
     771    27699771 :   for (i=2; i<lx; i++) z = addii(z, sqri(gel(x,i)));
     772     7432880 :   return gerepileuptoint(av,z);
     773             : }
     774             : 
     775             : GEN
     776    15829061 : ZV_dotproduct(GEN x,GEN y)
     777             : {
     778             :   long lx;
     779    15829061 :   if (x == y) return ZV_dotsquare(x);
     780    11104709 :   lx = lg(x);
     781    11104709 :   if (lx == 1) return gen_0;
     782    11104709 :   return ZV_dotproduct_i(x, y, lx);
     783             : }
     784             : 
     785             : static GEN
     786         294 : _ZM_mul(void *data /*ignored*/, GEN x, GEN y)
     787         294 : { (void)data; return ZM_mul(x,y); }
     788             : static GEN
     789        3381 : _ZM_sqr(void *data /*ignored*/, GEN x)
     790        3381 : { (void)data; return ZM_sqr(x); }
     791             : GEN
     792           0 : ZM_pow(GEN x, GEN n)
     793             : {
     794           0 :   pari_sp av = avma;
     795           0 :   if (!signe(n)) return matid(lg(x)-1);
     796           0 :   return gerepilecopy(av, gen_pow_i(x, n, NULL, &_ZM_sqr, &_ZM_mul));
     797             : }
     798             : GEN
     799        2947 : ZM_powu(GEN x, ulong n)
     800             : {
     801        2947 :   pari_sp av = avma;
     802        2947 :   if (!n) return matid(lg(x)-1);
     803        2947 :   return gerepilecopy(av, gen_powu_i(x, n, NULL, &_ZM_sqr, &_ZM_mul));
     804             : }
     805             : /********************************************************************/
     806             : /**                                                                **/
     807             : /**                           ADD, SUB                             **/
     808             : /**                                                                **/
     809             : /********************************************************************/
     810             : static GEN
     811    21072928 : ZC_add_i(GEN x, GEN y, long lx)
     812             : {
     813    21072928 :   GEN A = cgetg(lx, t_COL);
     814             :   long i;
     815   417500960 :   for (i=1; i<lx; i++) gel(A,i) = addii(gel(x,i), gel(y,i));
     816    21072380 :   return A;
     817             : }
     818             : GEN
     819    11576145 : ZC_add(GEN x, GEN y) { return ZC_add_i(x, y, lg(x)); }
     820             : GEN
     821      145880 : ZC_Z_add(GEN x, GEN y)
     822             : {
     823      145880 :   long k, lx = lg(x);
     824      145880 :   GEN z = cgetg(lx, t_COL);
     825      145880 :   if (lx == 1) pari_err_TYPE2("+",x,y);
     826      145880 :   gel(z,1) = addii(y,gel(x,1));
     827      848596 :   for (k = 2; k < lx; k++) gel(z,k) = icopy(gel(x,k));
     828      145880 :   return z;
     829             : }
     830             : 
     831             : static GEN
     832     2568493 : ZC_sub_i(GEN x, GEN y, long lx)
     833             : {
     834             :   long i;
     835     2568493 :   GEN A = cgetg(lx, t_COL);
     836    25358497 :   for (i=1; i<lx; i++) gel(A,i) = subii(gel(x,i), gel(y,i));
     837     2568492 :   return A;
     838             : }
     839             : GEN
     840     2430560 : ZC_sub(GEN x, GEN y) { return ZC_sub_i(x, y, lg(x)); }
     841             : GEN
     842           0 : ZC_Z_sub(GEN x, GEN y)
     843             : {
     844           0 :   long k, lx = lg(x);
     845           0 :   GEN z = cgetg(lx, t_COL);
     846           0 :   if (lx == 1) pari_err_TYPE2("+",x,y);
     847           0 :   gel(z,1) = subii(gel(x,1), y);
     848           0 :   for (k = 2; k < lx; k++) gel(z,k) = icopy(gel(x,k));
     849           0 :   return z;
     850             : }
     851             : GEN
     852      117526 : Z_ZC_sub(GEN a, GEN x)
     853             : {
     854      117526 :   long k, lx = lg(x);
     855      117526 :   GEN z = cgetg(lx, t_COL);
     856      117535 :   if (lx == 1) pari_err_TYPE2("-",a,x);
     857      117535 :   gel(z,1) = subii(a, gel(x,1));
     858      325088 :   for (k = 2; k < lx; k++) gel(z,k) = negi(gel(x,k));
     859      117540 :   return z;
     860             : }
     861             : 
     862             : GEN
     863      669619 : ZM_add(GEN x, GEN y)
     864             : {
     865      669619 :   long lx = lg(x), l, j;
     866             :   GEN z;
     867      669619 :   if (lx == 1) return cgetg(1, t_MAT);
     868      661401 :   z = cgetg(lx, t_MAT); l = lgcols(x);
     869    10157783 :   for (j = 1; j < lx; j++) gel(z,j) = ZC_add_i(gel(x,j), gel(y,j), l);
     870      661401 :   return z;
     871             : }
     872             : GEN
     873       52680 : ZM_sub(GEN x, GEN y)
     874             : {
     875       52680 :   long lx = lg(x), l, j;
     876             :   GEN z;
     877       52680 :   if (lx == 1) return cgetg(1, t_MAT);
     878       52680 :   z = cgetg(lx, t_MAT); l = lgcols(x);
     879      190613 :   for (j = 1; j < lx; j++) gel(z,j) = ZC_sub_i(gel(x,j), gel(y,j), l);
     880       52680 :   return z;
     881             : }
     882             : /********************************************************************/
     883             : /**                                                                **/
     884             : /**                         LINEAR COMBINATION                     **/
     885             : /**                                                                **/
     886             : /********************************************************************/
     887             : /* return X/c assuming division is exact */
     888             : GEN
     889     3754139 : ZC_Z_divexact(GEN x, GEN c)
     890    53594560 : { pari_APPLY_type(t_COL, diviiexact(gel(x,i), c)) }
     891             : GEN
     892           0 : ZC_u_divexact(GEN x, ulong c)
     893           0 : { pari_APPLY_type(t_COL, diviuexact(gel(x,i), c)) }
     894             : 
     895             : GEN
     896      606532 : ZM_Z_divexact(GEN x, GEN c)
     897     3699228 : { pari_APPLY_same(ZC_Z_divexact(gel(x,i), c)) }
     898             : 
     899             : GEN
     900    18387489 : ZC_Z_mul(GEN x, GEN c)
     901             : {
     902    18387489 :   if (!signe(c)) return zerocol(lg(x)-1);
     903    17911477 :   if (is_pm1(c)) return (signe(c) > 0)? ZC_copy(x): ZC_neg(x);
     904   212116001 :   pari_APPLY_type(t_COL, mulii(gel(x,i), c))
     905             : }
     906             : 
     907             : GEN
     908       24253 : ZC_z_mul(GEN x, long c)
     909             : {
     910       24253 :   if (!c) return zerocol(lg(x)-1);
     911       16931 :   if (c == 1) return ZC_copy(x);
     912       12360 :   if (c ==-1) return ZC_neg(x);
     913      168273 :   pari_APPLY_type(t_COL, mulsi(c, gel(x,i)))
     914             : }
     915             : 
     916             : GEN
     917       18072 : zv_z_mul(GEN x, long n)
     918      414886 : { pari_APPLY_long(x[i]*n) }
     919             : 
     920             : /* return a ZM */
     921             : GEN
     922           0 : nm_Z_mul(GEN X, GEN c)
     923             : {
     924           0 :   long i, j, h, l = lg(X), s = signe(c);
     925             :   GEN A;
     926           0 :   if (l == 1) return cgetg(1, t_MAT);
     927           0 :   h = lgcols(X);
     928           0 :   if (!s) return zeromat(h-1, l-1);
     929           0 :   if (is_pm1(c)) {
     930           0 :     if (s > 0) return Flm_to_ZM(X);
     931           0 :     X = Flm_to_ZM(X); ZM_togglesign(X); return X;
     932             :   }
     933           0 :   A = cgetg(l, t_MAT);
     934           0 :   for (j = 1; j < l; j++)
     935             :   {
     936           0 :     GEN a = cgetg(h, t_COL), x = gel(X, j);
     937           0 :     for (i = 1; i < h; i++) gel(a,i) = muliu(c, x[i]);
     938           0 :     gel(A,j) = a;
     939             :   }
     940           0 :   return A;
     941             : }
     942             : GEN
     943     1508489 : ZM_Z_mul(GEN X, GEN c)
     944             : {
     945     1508489 :   long i, j, h, l = lg(X);
     946             :   GEN A;
     947     1508489 :   if (l == 1) return cgetg(1, t_MAT);
     948     1508489 :   h = lgcols(X);
     949     1508489 :   if (!signe(c)) return zeromat(h-1, l-1);
     950     1506956 :   if (is_pm1(c)) return (signe(c) > 0)? ZM_copy(X): ZM_neg(X);
     951     1269341 :   A = cgetg(l, t_MAT);
     952    12469689 :   for (j = 1; j < l; j++)
     953             :   {
     954    11200348 :     GEN a = cgetg(h, t_COL), x = gel(X, j);
     955   226044464 :     for (i = 1; i < h; i++) gel(a,i) = mulii(c, gel(x,i));
     956    11200348 :     gel(A,j) = a;
     957             :   }
     958     1269341 :   return A;
     959             : }
     960             : void
     961    59119047 : ZC_lincomb1_inplace_i(GEN X, GEN Y, GEN v, long n)
     962             : {
     963             :   long i;
     964  1304665144 :   for (i = n; i; i--) gel(X,i) = addmulii_inplace(gel(X,i), gel(Y,i), v);
     965    59111310 : }
     966             : /* X <- X + v Y (elementary col operation) */
     967             : void
     968    50851590 : ZC_lincomb1_inplace(GEN X, GEN Y, GEN v)
     969             : {
     970    50851590 :   if (lgefint(v) != 2) return ZC_lincomb1_inplace_i(X, Y, v, lg(X)-1);
     971             : }
     972             : void
     973    12775652 : Flc_lincomb1_inplace(GEN X, GEN Y, ulong v, ulong q)
     974             : {
     975             :   long i;
     976    12775652 :   if (!v) return; /* v = 0 */
     977   353947324 :   for (i = lg(X)-1; i; i--) X[i] = Fl_add(X[i], Fl_mul(Y[i], v, q), q);
     978             : }
     979             : 
     980             : /* X + v Y, wasteful if (v = 0) */
     981             : static GEN
     982     3385115 : ZC_lincomb1(GEN v, GEN x, GEN y)
     983    45654655 : { pari_APPLY_type(t_COL, addmulii(gel(x,i), gel(y,i), v)) }
     984             : 
     985             : /* -X + vY */
     986             : static GEN
     987      232335 : ZC_lincomb_1(GEN v, GEN x, GEN y)
     988     1534641 : { pari_APPLY_type(t_COL, mulsubii(gel(y,i), v, gel(x,i))) }
     989             : 
     990             : /* X,Y compatible ZV; u,v in Z. Returns A = u*X + v*Y */
     991             : GEN
     992    12704549 : ZC_lincomb(GEN u, GEN v, GEN X, GEN Y)
     993             : {
     994             :   long su, sv;
     995             :   GEN A;
     996             : 
     997    12704549 :   su = signe(u); if (!su) return ZC_Z_mul(Y, v);
     998    12704227 :   sv = signe(v); if (!sv) return ZC_Z_mul(X, u);
     999    12574363 :   if (is_pm1(v))
    1000             :   {
    1001     2933453 :     if (is_pm1(u))
    1002             :     {
    1003     2508365 :       if (su != sv) A = ZC_sub(X, Y);
    1004      677496 :       else          A = ZC_add(X, Y);
    1005     2508365 :       if (su < 0) ZV_togglesign(A); /* in place but was created above */
    1006             :     }
    1007             :     else
    1008             :     {
    1009      425088 :       if (sv > 0) A = ZC_lincomb1 (u, Y, X);
    1010      176823 :       else        A = ZC_lincomb_1(u, Y, X);
    1011             :     }
    1012             :   }
    1013     9640988 :   else if (is_pm1(u))
    1014             :   {
    1015     3192362 :     if (su > 0) A = ZC_lincomb1 (v, X, Y);
    1016       55512 :     else        A = ZC_lincomb_1(v, X, Y);
    1017             :   }
    1018             :   else
    1019             :   { /* not cgetg_copy: x may be a t_VEC */
    1020     6448657 :     long i, lx = lg(X);
    1021     6448657 :     A = cgetg(lx,t_COL);
    1022    30819921 :     for (i=1; i<lx; i++) gel(A,i) = lincombii(u,v,gel(X,i),gel(Y,i));
    1023             :   }
    1024    12573967 :   return A;
    1025             : }
    1026             : 
    1027             : /********************************************************************/
    1028             : /**                                                                **/
    1029             : /**                           CONVERSIONS                          **/
    1030             : /**                                                                **/
    1031             : /********************************************************************/
    1032             : GEN
    1033       77159 : ZV_to_nv(GEN x)
    1034      231355 : { pari_APPLY_ulong(itou(gel(x,i))) }
    1035             : 
    1036             : GEN
    1037       37574 : zm_to_ZM(GEN x)
    1038      412556 : { pari_APPLY_type(t_MAT, zc_to_ZC(gel(x,i))) }
    1039             : 
    1040             : GEN
    1041         112 : zmV_to_ZMV(GEN x)
    1042         812 : { pari_APPLY_type(t_VEC, zm_to_ZM(gel(x,i))) }
    1043             : 
    1044             : /* same as Flm_to_ZM but do not assume positivity */
    1045             : GEN
    1046         980 : ZM_to_zm(GEN x)
    1047       17346 : { pari_APPLY_same(ZV_to_zv(gel(x,i))) }
    1048             : 
    1049             : GEN
    1050      284480 : zv_to_Flv(GEN x, ulong p)
    1051     4486328 : { pari_APPLY_ulong(umodsu(x[i], p)) }
    1052             : 
    1053             : GEN
    1054       19894 : zm_to_Flm(GEN x, ulong p)
    1055      304374 : { pari_APPLY_same(zv_to_Flv(gel(x,i),p)) }
    1056             : 
    1057             : GEN
    1058          35 : ZMV_to_zmV(GEN x)
    1059         357 : { pari_APPLY_type(t_VEC, ZM_to_zm(gel(x,i))) }
    1060             : 
    1061             : /********************************************************************/
    1062             : /**                                                                **/
    1063             : /**                         COPY, NEGATION                         **/
    1064             : /**                                                                **/
    1065             : /********************************************************************/
    1066             : GEN
    1067     5208761 : ZC_copy(GEN x)
    1068             : {
    1069     5208761 :   long i, lx = lg(x);
    1070     5208761 :   GEN y = cgetg(lx, t_COL);
    1071    47695527 :   for (i=1; i<lx; i++)
    1072             :   {
    1073    42486758 :     GEN c = gel(x,i);
    1074    42486758 :     gel(y,i) = lgefint(c) == 2? gen_0: icopy(c);
    1075             :   }
    1076     5208769 :   return y;
    1077             : }
    1078             : 
    1079             : GEN
    1080      279306 : ZM_copy(GEN x)
    1081     2652864 : { pari_APPLY_same(ZC_copy(gel(x,i))) }
    1082             : 
    1083             : void
    1084       68923 : ZV_neg_inplace(GEN M)
    1085             : {
    1086       68923 :   long l = lg(M);
    1087      287742 :   while (--l > 0) gel(M,l) = negi(gel(M,l));
    1088       68923 : }
    1089             : GEN
    1090     1834266 : ZC_neg(GEN x)
    1091    16782579 : { pari_APPLY_type(t_COL, negi(gel(x,i))) }
    1092             : GEN
    1093       31699 : zv_neg(GEN x)
    1094      369896 : { pari_APPLY_long(-x[i]) }
    1095             : GEN
    1096         126 : zv_neg_inplace(GEN M)
    1097             : {
    1098         126 :   long l = lg(M);
    1099         534 :   while (--l > 0) M[l] = -M[l];
    1100         126 :   return M;
    1101             : }
    1102             : GEN
    1103      303865 : ZM_neg(GEN x)
    1104     1250300 : { pari_APPLY_same(ZC_neg(gel(x,i))) }
    1105             : 
    1106             : void
    1107     1528179 : ZV_togglesign(GEN M)
    1108             : {
    1109     1528179 :   long l = lg(M);
    1110    34144744 :   while (--l > 0) togglesign_safe(&gel(M,l));
    1111     1528179 : }
    1112             : void
    1113           0 : ZM_togglesign(GEN M)
    1114             : {
    1115           0 :   long l = lg(M);
    1116           0 :   while (--l > 0) ZV_togglesign(gel(M,l));
    1117           0 : }
    1118             : 
    1119             : /********************************************************************/
    1120             : /**                                                                **/
    1121             : /**                        "DIVISION" mod HNF                      **/
    1122             : /**                                                                **/
    1123             : /********************************************************************/
    1124             : /* Reduce ZC x modulo ZM y in HNF, may return x itself (not a copy) */
    1125             : GEN
    1126     1756967 : ZC_hnfremdiv(GEN x, GEN y, GEN *Q)
    1127             : {
    1128     1756967 :   long i, l = lg(x);
    1129             :   GEN q;
    1130             : 
    1131     1756967 :   if (Q) *Q = cgetg(l,t_COL);
    1132     1756967 :   if (l == 1) return cgetg(1,t_COL);
    1133    13906902 :   for (i = l-1; i>0; i--)
    1134             :   {
    1135    12149935 :     q = diviiround(gel(x,i), gcoeff(y,i,i));
    1136    12149935 :     if (signe(q)) {
    1137     5231505 :       togglesign(q);
    1138     5231505 :       x = ZC_lincomb(gen_1, q, x, gel(y,i));
    1139             :     }
    1140    12149935 :     if (Q) gel(*Q, i) = q;
    1141             :   }
    1142     1756967 :   return x;
    1143             : }
    1144             : 
    1145             : /* x = y Q + R, may return some columns of x (not copies) */
    1146             : GEN
    1147       66246 : ZM_hnfdivrem(GEN x, GEN y, GEN *Q)
    1148             : {
    1149       66246 :   long lx = lg(x), i;
    1150       66246 :   GEN R = cgetg(lx, t_MAT);
    1151       66246 :   if (Q)
    1152             :   {
    1153       24080 :     GEN q = cgetg(lx, t_MAT); *Q = q;
    1154       45614 :     for (i=1; i<lx; i++) gel(R,i) = ZC_hnfremdiv(gel(x,i),y,(GEN*)(q+i));
    1155             :   }
    1156             :   else
    1157      109538 :     for (i=1; i<lx; i++)
    1158             :     {
    1159       67372 :       pari_sp av = avma;
    1160       67372 :       GEN z = ZC_hnfrem(gel(x,i),y);
    1161       67372 :       gel(R,i) = (avma == av)? ZC_copy(z): gerepileupto(av, z);
    1162             :     }
    1163       66246 :   return R;
    1164             : }
    1165             : 
    1166             : /********************************************************************/
    1167             : /**                                                                **/
    1168             : /**                               TESTS                            **/
    1169             : /**                                                                **/
    1170             : /********************************************************************/
    1171             : int
    1172    21185244 : zv_equal0(GEN V)
    1173             : {
    1174    21185244 :   long l = lg(V);
    1175    34376027 :   while (--l > 0)
    1176    28050224 :     if (V[l]) return 0;
    1177     6325803 :   return 1;
    1178             : }
    1179             : 
    1180             : int
    1181     1028193 : ZV_equal0(GEN V)
    1182             : {
    1183     1028193 :   long l = lg(V);
    1184     4129431 :   while (--l > 0)
    1185     3930165 :     if (signe(gel(V,l))) return 0;
    1186      199266 :   return 1;
    1187             : }
    1188             : int
    1189       14170 : ZMrow_equal0(GEN V, long i)
    1190             : {
    1191       14170 :   long l = lg(V);
    1192       27855 :   while (--l > 0)
    1193       22769 :     if (signe(gcoeff(V,i,l))) return 0;
    1194        5086 :   return 1;
    1195             : }
    1196             : 
    1197             : static int
    1198    15955733 : ZV_equal_lg(GEN V, GEN W, long l)
    1199             : {
    1200    26150330 :   while (--l > 0)
    1201    24057108 :     if (!equalii(gel(V,l), gel(W,l))) return 0;
    1202     2093222 :   return 1;
    1203             : }
    1204             : int
    1205    14620393 : ZV_equal(GEN V, GEN W)
    1206             : {
    1207    14620393 :   long l = lg(V);
    1208    14620393 :   if (lg(W) != l) return 0;
    1209    14620393 :   return ZV_equal_lg(V, W, l);
    1210             : }
    1211             : int
    1212      467827 : ZM_equal(GEN A, GEN B)
    1213             : {
    1214      467827 :   long i, m, l = lg(A);
    1215      467827 :   if (lg(B) != l) return 0;
    1216      467673 :   if (l == 1) return 1;
    1217      467673 :   m = lgcols(A);
    1218      467673 :   if (lgcols(B) != m) return 0;
    1219     1731826 :   for (i = 1; i < l; i++)
    1220     1335340 :     if (!ZV_equal_lg(gel(A,i), gel(B,i), m)) return 0;
    1221      396486 :   return 1;
    1222             : }
    1223             : int
    1224       69304 : ZM_equal0(GEN A)
    1225             : {
    1226       69304 :   long i, j, m, l = lg(A);
    1227       69304 :   if (l == 1) return 1;
    1228       69304 :   m = lgcols(A);
    1229      118339 :   for (j = 1; j < l; j++)
    1230     2582332 :     for (i = 1; i < m; i++)
    1231     2533297 :       if (signe(gcoeff(A,i,j))) return 0;
    1232       14636 :   return 1;
    1233             : }
    1234             : int
    1235    30134249 : zv_equal(GEN V, GEN W)
    1236             : {
    1237    30134249 :   long l = lg(V);
    1238    30134249 :   if (lg(W) != l) return 0;
    1239   254147199 :   while (--l > 0)
    1240   224530614 :     if (V[l] != W[l]) return 0;
    1241    29616585 :   return 1;
    1242             : }
    1243             : 
    1244             : int
    1245        1337 : zvV_equal(GEN V, GEN W)
    1246             : {
    1247        1337 :   long l = lg(V);
    1248        1337 :   if (lg(W) != l) return 0;
    1249       77707 :   while (--l > 0)
    1250       77301 :     if (!zv_equal(gel(V,l),gel(W,l))) return 0;
    1251         406 :   return 1;
    1252             : }
    1253             : 
    1254             : int
    1255      156175 : ZM_ishnf(GEN x)
    1256             : {
    1257      156175 :   long i,j, lx = lg(x);
    1258      443875 :   for (i=1; i<lx; i++)
    1259             :   {
    1260      329028 :     GEN xii = gcoeff(x,i,i);
    1261      329028 :     if (signe(xii) <= 0) return 0;
    1262      631292 :     for (j=1; j<i; j++)
    1263      327429 :       if (signe(gcoeff(x,i,j))) return 0;
    1264      660594 :     for (j=i+1; j<lx; j++)
    1265             :     {
    1266      372894 :       GEN xij = gcoeff(x,i,j);
    1267      372894 :       if (signe(xij)<0 || cmpii(xij,xii)>=0) return 0;
    1268             :     }
    1269             :   }
    1270      114847 :   return 1;
    1271             : }
    1272             : int
    1273      123323 : ZM_isidentity(GEN x)
    1274             : {
    1275      123323 :   long i,j, lx = lg(x);
    1276             : 
    1277      123323 :   if (lx == 1) return 1;
    1278      123316 :   if (lx != lgcols(x)) return 0;
    1279      886047 :   for (j=1; j<lx; j++)
    1280             :   {
    1281      762780 :     GEN c = gel(x,j);
    1282     3543955 :     for (i=1; i<j; )
    1283     2781196 :       if (signe(gel(c,i++))) return 0;
    1284             :     /* i = j */
    1285      762759 :     if (!equali1(gel(c,i++))) return 0;
    1286     3543969 :     for (   ; i<lx; )
    1287     2781231 :       if (signe(gel(c,i++))) return 0;
    1288             :   }
    1289      123267 :   return 1;
    1290             : }
    1291             : int
    1292       75745 : ZM_isdiagonal(GEN x)
    1293             : {
    1294       75745 :   long i,j, lx = lg(x);
    1295       75745 :   if (lx == 1) return 1;
    1296       75745 :   if (lx != lgcols(x)) return 0;
    1297             : 
    1298      190698 :   for (j=1; j<lx; j++)
    1299             :   {
    1300      157147 :     GEN c = gel(x,j);
    1301      222394 :     for (i=1; i<j; i++)
    1302      107406 :       if (signe(gel(c,i))) return 0;
    1303      279160 :     for (i++; i<lx; i++)
    1304      164207 :       if (signe(gel(c,i))) return 0;
    1305             :   }
    1306       33551 :   return 1;
    1307             : }
    1308             : int
    1309       48105 : ZM_isscalar(GEN x, GEN s)
    1310             : {
    1311       48105 :   long i, j, lx = lg(x);
    1312             : 
    1313       48105 :   if (lx == 1) return 1;
    1314       48105 :   if (!s) s = gcoeff(x,1,1);
    1315       48105 :   if (equali1(s)) return ZM_isidentity(x);
    1316       48105 :   if (lx != lgcols(x)) return 0;
    1317      256069 :   for (j=1; j<lx; j++)
    1318             :   {
    1319      235109 :     GEN c = gel(x,j);
    1320     1674671 :     for (i=1; i<j; )
    1321     1465766 :       if (signe(gel(c,i++))) return 0;
    1322             :     /* i = j */
    1323      208905 :     if (!equalii(gel(c,i++), s)) return 0;
    1324     1696340 :     for (   ; i<lx; )
    1325     1488376 :       if (signe(gel(c,i++))) return 0;
    1326             :   }
    1327       20960 :   return 1;
    1328             : }
    1329             : 
    1330             : long
    1331      102870 : ZC_is_ei(GEN x)
    1332             : {
    1333      102870 :   long i, j = 0, l = lg(x);
    1334      969907 :   for (i = 1; i < l; i++)
    1335             :   {
    1336      867038 :     GEN c = gel(x,i);
    1337      867038 :     long s = signe(c);
    1338      867038 :     if (!s) continue;
    1339      102864 :     if (s < 0 || !is_pm1(c) || j) return 0;
    1340      102863 :     j = i;
    1341             :   }
    1342      102869 :   return j;
    1343             : }
    1344             : 
    1345             : /********************************************************************/
    1346             : /**                                                                **/
    1347             : /**                       MISCELLANEOUS                            **/
    1348             : /**                                                                **/
    1349             : /********************************************************************/
    1350             : /* assume lg(x) = lg(y), x,y in Z^n */
    1351             : int
    1352      871291 : ZV_cmp(GEN x, GEN y)
    1353             : {
    1354      871291 :   long fl,i, lx = lg(x);
    1355     1609741 :   for (i=1; i<lx; i++)
    1356     1317067 :     if (( fl = cmpii(gel(x,i), gel(y,i)) )) return fl;
    1357      292674 :   return 0;
    1358             : }
    1359             : /* assume lg(x) = lg(y), x,y in Z^n */
    1360             : int
    1361        4607 : ZV_abscmp(GEN x, GEN y)
    1362             : {
    1363        4607 :   long fl,i, lx = lg(x);
    1364       22108 :   for (i=1; i<lx; i++)
    1365       22028 :     if (( fl = abscmpii(gel(x,i), gel(y,i)) )) return fl;
    1366          80 :   return 0;
    1367             : }
    1368             : 
    1369             : long
    1370     4706688 : zv_content(GEN x)
    1371             : {
    1372     4706688 :   long i, s, l = lg(x);
    1373     4706688 :   if (l == 1) return 0;
    1374     4706681 :   s = labs(x[1]);
    1375     8747047 :   for (i = 2; i < l && s != 1; i++) s = ugcd(s, labs(x[i]));
    1376     4706681 :   return s;
    1377             : }
    1378             : GEN
    1379      280508 : ZV_content(GEN x)
    1380             : {
    1381      280508 :   long i, l = lg(x);
    1382      280508 :   pari_sp av = avma;
    1383             :   GEN c;
    1384      280508 :   if (l == 1) return gen_0;
    1385      280508 :   if (l == 2) return absi(gel(x,1));
    1386      212202 :   c = gel(x,1);
    1387      513870 :   for (i = 2; i < l; i++)
    1388             :   {
    1389      370120 :     c = gcdii(c, gel(x,i));
    1390      370120 :     if (is_pm1(c)) { set_avma(av); return gen_1; }
    1391             :   }
    1392      143750 :   return gerepileuptoint(av, c);
    1393             : }
    1394             : 
    1395             : GEN
    1396     1085290 : ZM_det_triangular(GEN mat)
    1397             : {
    1398             :   pari_sp av;
    1399     1085290 :   long i,l = lg(mat);
    1400             :   GEN s;
    1401             : 
    1402     1085290 :   if (l<3) return l<2? gen_1: icopy(gcoeff(mat,1,1));
    1403      976861 :   av = avma; s = gcoeff(mat,1,1);
    1404     3156337 :   for (i=2; i<l; i++) s = mulii(s,gcoeff(mat,i,i));
    1405      976861 :   return gerepileuptoint(av,s);
    1406             : }
    1407             : 
    1408             : /* assumes no overflow */
    1409             : long
    1410      397537 : zv_prod(GEN v)
    1411             : {
    1412      397537 :   long n, i, l = lg(v);
    1413      397537 :   if (l == 1) return 1;
    1414      436653 :   n = v[1]; for (i = 2; i < l; i++) n *= v[i];
    1415      297332 :   return n;
    1416             : }
    1417             : 
    1418             : static GEN
    1419   278326003 : _mulii(void *E, GEN a, GEN b)
    1420   278326003 : { (void) E; return mulii(a, b); }
    1421             : 
    1422             : /* product of ulongs */
    1423             : GEN
    1424      375000 : zv_prod_Z(GEN v)
    1425             : {
    1426      375000 :   pari_sp av = avma;
    1427      375000 :   long k, m, n = lg(v)-1;
    1428             :   GEN V;
    1429      375000 :   switch(n) {
    1430       19607 :     case 0: return gen_1;
    1431       59787 :     case 1: return utoi(v[1]);
    1432      130683 :     case 2: return muluu(v[1], v[2]);
    1433             :   }
    1434      164923 :   m = n >> 1;
    1435      164923 :   V = cgetg(m + (odd(n)? 2: 1), t_VEC);
    1436     6653327 :   for (k = 1; k <= m; k++) gel(V,k) = muluu(v[k<<1], v[(k<<1)-1]);
    1437      164923 :   if (odd(n)) gel(V,k) = utoipos(v[n]);
    1438      164923 :   return gerepileuptoint(av, gen_product(V, NULL, &_mulii));
    1439             : }
    1440             : GEN
    1441    14694393 : vecsmall_prod(GEN v)
    1442             : {
    1443    14694393 :   pari_sp av = avma;
    1444    14694393 :   long k, m, n = lg(v)-1;
    1445             :   GEN V;
    1446    14694393 :   switch (n) {
    1447           0 :     case 0: return gen_1;
    1448           0 :     case 1: return stoi(v[1]);
    1449          21 :     case 2: return mulss(v[1], v[2]);
    1450             :   }
    1451    14694372 :   m = n >> 1;
    1452    14694372 :   V = cgetg(m + (odd(n)? 2: 1), t_VEC);
    1453   161556906 :   for (k = 1; k <= m; k++) gel(V,k) = mulss(v[k<<1], v[(k<<1)-1]);
    1454    14694372 :   if (odd(n)) gel(V,k) = stoi(v[n]);
    1455    14694372 :   return gerepileuptoint(av, gen_product(V, NULL, &_mulii));
    1456             : }
    1457             : 
    1458             : GEN
    1459     7607339 : ZV_prod(GEN v)
    1460             : {
    1461     7607339 :   pari_sp av = avma;
    1462     7607339 :   long i, l = lg(v);
    1463             :   GEN n;
    1464     7607339 :   if (l == 1) return gen_1;
    1465     7584603 :   if (l > 7) return gerepileuptoint(av, gen_product(v, NULL, _mulii));
    1466      242891 :   n = gel(v,1);
    1467      242891 :   if (l == 2) return icopy(n);
    1468      365437 :   for (i = 2; i < l; i++) n = mulii(n, gel(v,i));
    1469      145156 :   return gerepileuptoint(av, n);
    1470             : }
    1471             : /* assumes no overflow */
    1472             : long
    1473       11249 : zv_sum(GEN v)
    1474             : {
    1475       11249 :   long n, i, l = lg(v);
    1476       11249 :   if (l == 1) return 0;
    1477       64414 :   n = v[1]; for (i = 2; i < l; i++) n += v[i];
    1478       11228 :   return n;
    1479             : }
    1480             : /* assumes no overflow and 0 <= n <= #v */
    1481             : long
    1482           0 : zv_sumpart(GEN v, long n)
    1483             : {
    1484             :   long i, p;
    1485           0 :   if (!n) return 0;
    1486           0 :   p = v[1]; for (i = 2; i <= n; i++) p += v[i];
    1487           0 :   return p;
    1488             : }
    1489             : GEN
    1490          63 : ZV_sum(GEN v)
    1491             : {
    1492          63 :   pari_sp av = avma;
    1493          63 :   long i, l = lg(v);
    1494             :   GEN n;
    1495          63 :   if (l == 1) return gen_0;
    1496          63 :   n = gel(v,1);
    1497          63 :   if (l == 2) return icopy(n);
    1498         532 :   for (i = 2; i < l; i++) n = addii(n, gel(v,i));
    1499          63 :   return gerepileuptoint(av, n);
    1500             : }
    1501             : 
    1502             : /********************************************************************/
    1503             : /**                                                                **/
    1504             : /**         GRAM SCHMIDT REDUCTION (integer matrices)              **/
    1505             : /**                                                                **/
    1506             : /********************************************************************/
    1507             : 
    1508             : /* L[k,] += q * L[l,], l < k. Inefficient if q = 0 */
    1509             : static void
    1510       70557 : Zupdate_row(long k, long l, GEN q, GEN L, GEN B)
    1511             : {
    1512       70557 :   long i, qq = itos_or_0(q);
    1513       70557 :   if (!qq)
    1514             :   {
    1515        6932 :     for(i=1;i<l;i++)  gcoeff(L,k,i) = addii(gcoeff(L,k,i),mulii(q,gcoeff(L,l,i)));
    1516        1349 :     gcoeff(L,k,l) = addii(gcoeff(L,k,l), mulii(q,B));
    1517        1349 :     return;
    1518             :   }
    1519       69208 :   if (qq == 1) {
    1520       31250 :     for (i=1;i<l; i++) gcoeff(L,k,i) = addii(gcoeff(L,k,i),gcoeff(L,l,i));
    1521       19758 :     gcoeff(L,k,l) = addii(gcoeff(L,k,l), B);
    1522       49450 :   } else if (qq == -1) {
    1523       35509 :     for (i=1;i<l; i++) gcoeff(L,k,i) = subii(gcoeff(L,k,i),gcoeff(L,l,i));
    1524       19024 :     gcoeff(L,k,l) = addii(gcoeff(L,k,l), negi(B));
    1525             :   } else {
    1526       65028 :     for(i=1;i<l;i++) gcoeff(L,k,i) = addii(gcoeff(L,k,i),mulsi(qq,gcoeff(L,l,i)));
    1527       30426 :     gcoeff(L,k,l) = addii(gcoeff(L,k,l), mulsi(qq,B));
    1528             :   }
    1529             : }
    1530             : 
    1531             : /* update L[k,] */
    1532             : static void
    1533      175980 : ZRED(long k, long l, GEN x, GEN L, GEN B)
    1534             : {
    1535      175980 :   GEN q = truedivii(addii(B,shifti(gcoeff(L,k,l),1)), shifti(B,1));
    1536      175980 :   if (!signe(q)) return;
    1537       70557 :   q = negi(q);
    1538       70557 :   Zupdate_row(k,l,q,L,B);
    1539       70557 :   gel(x,k) = ZC_lincomb(gen_1, q, gel(x,k), gel(x,l));
    1540             : }
    1541             : 
    1542             : /* Gram-Schmidt reduction, x a ZM */
    1543             : static void
    1544      201390 : ZincrementalGS(GEN x, GEN L, GEN B, long k)
    1545             : {
    1546             :   long i, j;
    1547      660814 :   for (j=1; j<=k; j++)
    1548             :   {
    1549      459424 :     pari_sp av = avma;
    1550      459424 :     GEN u = ZV_dotproduct(gel(x,k), gel(x,j));
    1551      962031 :     for (i=1; i<j; i++)
    1552             :     {
    1553      502607 :       u = subii(mulii(gel(B,i+1), u), mulii(gcoeff(L,k,i), gcoeff(L,j,i)));
    1554      502607 :       u = diviiexact(u, gel(B,i));
    1555             :     }
    1556      459424 :     gcoeff(L,k,j) = gerepileuptoint(av, u);
    1557             :   }
    1558      201390 :   gel(B,k+1) = gcoeff(L,k,k); gcoeff(L,k,k) = gen_1;
    1559      201390 : }
    1560             : 
    1561             : /* Variant reducemodinvertible(ZC v, ZM y), when y singular.
    1562             :  * Very inefficient if y is not LLL-reduced of maximal rank */
    1563             : static GEN
    1564       21378 : ZC_reducemodmatrix_i(GEN v, GEN y)
    1565             : {
    1566       21378 :   GEN B, L, x = shallowconcat(y, v);
    1567       21378 :   long k, lx = lg(x), nx = lx-1;
    1568             : 
    1569       21378 :   B = scalarcol_shallow(gen_1, lx);
    1570       21378 :   L = zeromatcopy(nx, nx);
    1571       91847 :   for (k=1; k <= nx; k++) ZincrementalGS(x, L, B, k);
    1572       70469 :   for (k = nx-1; k >= 1; k--) ZRED(nx,k, x,L,gel(B,k+1));
    1573       21378 :   return gel(x,nx);
    1574             : }
    1575             : GEN
    1576       21378 : ZC_reducemodmatrix(GEN v, GEN y) {
    1577       21378 :   pari_sp av = avma;
    1578       21378 :   return gerepilecopy(av, ZC_reducemodmatrix_i(v,y));
    1579             : }
    1580             : 
    1581             : /* Variant reducemodinvertible(ZM v, ZM y), when y singular.
    1582             :  * Very inefficient if y is not LLL-reduced of maximal rank */
    1583             : static GEN
    1584       37198 : ZM_reducemodmatrix_i(GEN v, GEN y)
    1585             : {
    1586             :   GEN B, L, V;
    1587       37198 :   long j, k, lv = lg(v), nx = lg(y), lx = nx+1;
    1588             : 
    1589       37198 :   V = cgetg(lv, t_MAT);
    1590       37198 :   B = scalarcol_shallow(gen_1, lx);
    1591       37198 :   L = zeromatcopy(nx, nx);
    1592       97461 :   for (k=1; k < nx; k++) ZincrementalGS(y, L, B, k);
    1593      107856 :   for (j = 1; j < lg(v); j++)
    1594             :   {
    1595       70658 :     GEN x = shallowconcat(y, gel(v,j));
    1596       70658 :     ZincrementalGS(x, L, B, nx); /* overwrite last */
    1597      197547 :     for (k = nx-1; k >= 1; k--) ZRED(nx,k, x,L,gel(B,k+1));
    1598       70658 :     gel(V,j) = gel(x,nx);
    1599             :   }
    1600       37198 :   return V;
    1601             : }
    1602             : GEN
    1603       37198 : ZM_reducemodmatrix(GEN v, GEN y) {
    1604       37198 :   pari_sp av = avma;
    1605       37198 :   return gerepilecopy(av, ZM_reducemodmatrix_i(v,y));
    1606             : }
    1607             : 
    1608             : GEN
    1609       11256 : ZC_reducemodlll(GEN x,GEN y)
    1610             : {
    1611       11256 :   pari_sp av = avma;
    1612       11256 :   GEN z = ZC_reducemodmatrix(x, ZM_lll(y, 0.75, LLL_INPLACE));
    1613       11256 :   return gerepilecopy(av, z);
    1614             : }
    1615             : GEN
    1616           0 : ZM_reducemodlll(GEN x,GEN y)
    1617             : {
    1618           0 :   pari_sp av = avma;
    1619           0 :   GEN z = ZM_reducemodmatrix(x, ZM_lll(y, 0.75, LLL_INPLACE));
    1620           0 :   return gerepilecopy(av, z);
    1621             : }

Generated by: LCOV version 1.13