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.18.0 lcov report (development 29806-4d001396c7) Lines: 800 938 85.3 %
Date: 2024-12-21 09:08:57 Functions: 122 139 87.8 %
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; 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             : static int
      19     1932418 : check_ZV(GEN x, long l)
      20             : {
      21             :   long i;
      22    13037678 :   for (i=1; i<l; i++)
      23    11105316 :     if (typ(gel(x,i)) != t_INT) return 0;
      24     1932362 :   return 1;
      25             : }
      26             : void
      27     1423591 : RgV_check_ZV(GEN A, const char *s)
      28             : {
      29     1423591 :   if (!RgV_is_ZV(A)) pari_err_TYPE(stack_strcat(s," [integer vector]"), A);
      30     1423584 : }
      31             : void
      32      651157 : RgM_check_ZM(GEN A, const char *s)
      33             : {
      34      651157 :   long n = lg(A);
      35      651157 :   if (n != 1)
      36             :   {
      37      651024 :     long j, m = lgcols(A);
      38     2583387 :     for (j=1; j<n; j++)
      39     1932418 :       if (!check_ZV(gel(A,j), m))
      40          56 :         pari_err_TYPE(stack_strcat(s," [integer matrix]"), A);
      41             :   }
      42      651102 : }
      43             : 
      44             : /* assume m > 1 */
      45             : static long
      46   114302827 : ZV_max_lg_i(GEN x, long m)
      47             : {
      48   114302827 :   long i, l = lgefint(gel(x,1));
      49   890708561 :   for (i = 2; i < m; i++) l = maxss(l, lgefint(gel(x,i)));
      50   114302136 :   return l;
      51             : }
      52             : long
      53       10640 : ZV_max_lg(GEN x)
      54             : {
      55       10640 :   long m = lg(x);
      56       10640 :   return m == 1? 2: ZV_max_lg_i(x, m);
      57             : }
      58             : 
      59             : /* assume n > 1 and m > 1 */
      60             : static long
      61    28537222 : ZM_max_lg_i(GEN x, long n, long m)
      62             : {
      63    28537222 :   long j, l = ZV_max_lg_i(gel(x,1), m);
      64   114292895 :   for (j = 2; j < n; j++) l = maxss(l, ZV_max_lg_i(gel(x,j), m));
      65    28537401 :   return l;
      66             : }
      67             : long
      68       22184 : ZM_max_lg(GEN x)
      69             : {
      70       22184 :   long n = lg(x), m;
      71       22184 :   if (n == 1) return 2;
      72       22184 :   m = lgcols(x); return m == 1? 2: ZM_max_lg_i(x, n, m);
      73             : }
      74             : 
      75             : /* assume m > 1 */
      76             : static long
      77           0 : ZV_max_expi_i(GEN x, long m)
      78             : {
      79           0 :   long i, prec = expi(gel(x,1));
      80           0 :   for (i = 2; i < m; i++) prec = maxss(prec, expi(gel(x,i)));
      81           0 :   return prec;
      82             : }
      83             : long
      84           0 : ZV_max_expi(GEN x)
      85             : {
      86           0 :   long m = lg(x);
      87           0 :   return m == 1? 2: ZV_max_expi_i(x, m);
      88             : }
      89             : 
      90             : /* assume n > 1 and m > 1 */
      91             : static long
      92           0 : ZM_max_expi_i(GEN x, long n, long m)
      93             : {
      94           0 :   long j, prec = ZV_max_expi_i(gel(x,1), m);
      95           0 :   for (j = 2; j < n; j++) prec = maxss(prec, ZV_max_expi_i(gel(x,j), m));
      96           0 :   return prec;
      97             : }
      98             : long
      99           0 : ZM_max_expi(GEN x)
     100             : {
     101           0 :   long n = lg(x), m;
     102           0 :   if (n == 1) return 2;
     103           0 :   m = lgcols(x); return m == 1? 2: ZM_max_expi_i(x, n, m);
     104             : }
     105             : 
     106             : GEN
     107        4162 : ZM_supnorm(GEN x)
     108             : {
     109        4162 :   long i, j, h, lx = lg(x);
     110        4162 :   GEN s = gen_0;
     111        4162 :   if (lx == 1) return gen_1;
     112        4162 :   h = lgcols(x);
     113       25549 :   for (j=1; j<lx; j++)
     114             :   {
     115       21387 :     GEN xj = gel(x,j);
     116      290468 :     for (i=1; i<h; i++)
     117             :     {
     118      269081 :       GEN c = gel(xj,i);
     119      269081 :       if (abscmpii(c, s) > 0) s = c;
     120             :     }
     121             :   }
     122        4162 :   return absi(s);
     123             : }
     124             : 
     125             : /********************************************************************/
     126             : /**                                                                **/
     127             : /**                           MULTIPLICATION                       **/
     128             : /**                                                                **/
     129             : /********************************************************************/
     130             : /* x nonempty ZM, y a compatible nc (dimension > 0). */
     131             : static GEN
     132           0 : ZM_nc_mul_i(GEN x, GEN y, long c, long l)
     133             : {
     134             :   long i, j;
     135             :   pari_sp av;
     136           0 :   GEN z = cgetg(l,t_COL), s;
     137             : 
     138           0 :   for (i=1; i<l; i++)
     139             :   {
     140           0 :     av = avma; s = muliu(gcoeff(x,i,1),y[1]);
     141           0 :     for (j=2; j<c; j++)
     142           0 :       if (y[j]) s = addii(s, muliu(gcoeff(x,i,j),y[j]));
     143           0 :     gel(z,i) = gerepileuptoint(av,s);
     144             :   }
     145           0 :   return z;
     146             : }
     147             : 
     148             : /* x ZV, y a compatible zc. */
     149             : GEN
     150     2229528 : ZV_zc_mul(GEN x, GEN y)
     151             : {
     152     2229528 :   long j, l = lg(x);
     153     2229528 :   pari_sp av = avma;
     154     2229528 :   GEN s = mulis(gel(x,1),y[1]);
     155    50292998 :   for (j=2; j<l; j++)
     156    48063470 :     if (y[j]) s = addii(s, mulis(gel(x,j),y[j]));
     157     2229528 :   return gerepileuptoint(av,s);
     158             : }
     159             : 
     160             : /* x nonempty ZM, y a compatible zc (dimension > 0). */
     161             : static GEN
     162    19278718 : ZM_zc_mul_i(GEN x, GEN y, long c, long l)
     163             : {
     164             :   long i, j;
     165    19278718 :   GEN z = cgetg(l,t_COL);
     166             : 
     167   121762344 :   for (i=1; i<l; i++)
     168             :   {
     169   102484362 :     pari_sp av = avma;
     170   102484362 :     GEN s = mulis(gcoeff(x,i,1),y[1]);
     171  1171029800 :     for (j=2; j<c; j++)
     172  1068537022 :       if (y[j]) s = addii(s, mulis(gcoeff(x,i,j),y[j]));
     173   102492778 :     gel(z,i) = gerepileuptoint(av,s);
     174             :   }
     175    19277982 :   return z;
     176             : }
     177             : GEN
     178    17214552 : ZM_zc_mul(GEN x, GEN y) {
     179    17214552 :   long l = lg(x);
     180    17214552 :   if (l == 1) return cgetg(1, t_COL);
     181    17214552 :   return ZM_zc_mul_i(x,y, l, lgcols(x));
     182             : }
     183             : 
     184             : /* y nonempty ZM, x a compatible zv (dimension > 0). */
     185             : GEN
     186        1736 : zv_ZM_mul(GEN x, GEN y) {
     187        1736 :   long i,j, lx = lg(x), ly = lg(y);
     188             :   GEN z;
     189        1736 :   if (lx == 1) return zerovec(ly-1);
     190        1736 :   z = cgetg(ly,t_VEC);
     191        4046 :   for (j=1; j<ly; j++)
     192             :   {
     193        2310 :     pari_sp av = avma;
     194        2310 :     GEN s = mulsi(x[1], gcoeff(y,1,j));
     195        3990 :     for (i=2; i<lx; i++)
     196        1680 :       if (x[i]) s = addii(s, mulsi(x[i], gcoeff(y,i,j)));
     197        2310 :     gel(z,j) = gerepileuptoint(av,s);
     198             :   }
     199        1736 :   return z;
     200             : }
     201             : 
     202             : /* x ZM, y a compatible zm (dimension > 0). */
     203             : GEN
     204      975485 : ZM_zm_mul(GEN x, GEN y)
     205             : {
     206      975485 :   long j, c, l = lg(x), ly = lg(y);
     207      975485 :   GEN z = cgetg(ly, t_MAT);
     208      975485 :   if (l == 1) return z;
     209      975478 :   c = lgcols(x);
     210     3039615 :   for (j = 1; j < ly; j++) gel(z,j) = ZM_zc_mul_i(x, gel(y,j), l,c);
     211      975481 :   return z;
     212             : }
     213             : /* x ZM, y a compatible zn (dimension > 0). */
     214             : GEN
     215           0 : ZM_nm_mul(GEN x, GEN y)
     216             : {
     217           0 :   long j, c, l = lg(x), ly = lg(y);
     218           0 :   GEN z = cgetg(ly, t_MAT);
     219           0 :   if (l == 1) return z;
     220           0 :   c = lgcols(x);
     221           0 :   for (j = 1; j < ly; j++) gel(z,j) = ZM_nc_mul_i(x, gel(y,j), l,c);
     222           0 :   return z;
     223             : }
     224             : 
     225             : /* Strassen-Winograd algorithm */
     226             : 
     227             : /* Return A[ma+1..ma+da, na+1..na+ea] - B[mb+1..mb+db, nb+1..nb+eb]
     228             :  * as an (m x n)-matrix, padding the input with zeroes as necessary. */
     229             : static GEN
     230      586464 : add_slices(long m, long n,
     231             :            GEN A, long ma, long da, long na, long ea,
     232             :            GEN B, long mb, long db, long nb, long eb)
     233             : {
     234      586464 :   long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
     235      586464 :   GEN M = cgetg(n + 1, t_MAT), C;
     236             : 
     237     4926874 :   for (j = 1; j <= min_e; j++) {
     238     4340410 :     gel(M, j) = C = cgetg(m + 1, t_COL);
     239    84801437 :     for (i = 1; i <= min_d; i++)
     240    80461027 :       gel(C, i) = addii(gcoeff(A, ma + i, na + j),
     241    80461027 :                         gcoeff(B, mb + i, nb + j));
     242     4406779 :     for (; i <= da; i++) gel(C, i) = gcoeff(A, ma + i, na + j);
     243     4340410 :     for (; i <= db; i++) gel(C, i) = gcoeff(B, mb + i, nb + j);
     244     4340410 :     for (; i <= m; i++)  gel(C, i) = gen_0;
     245             :   }
     246      645860 :   for (; j <= ea; j++) {
     247       59396 :     gel(M, j) = C = cgetg(m + 1, t_COL);
     248      208387 :     for (i = 1; i <= da; i++) gel(C, i) = gcoeff(A, ma + i, na + j);
     249       59396 :     for (; i <= m; i++) gel(C, i) = gen_0;
     250             :   }
     251      586464 :   for (; j <= eb; j++) {
     252           0 :     gel(M, j) = C = cgetg(m + 1, t_COL);
     253           0 :     for (i = 1; i <= db; i++) gel(C, i) = gcoeff(B, mb + i, nb + j);
     254           0 :     for (; i <= m; i++) gel(C, i) = gen_0;
     255             :   }
     256      586464 :   for (; j <= n; j++) gel(M, j) = zerocol(m);
     257      586464 :   return M;
     258             : }
     259             : 
     260             : /* Return A[ma+1..ma+da, na+1..na+ea] - B[mb+1..mb+db, nb+1..nb+eb]
     261             :  * as an (m x n)-matrix, padding the input with zeroes as necessary. */
     262             : static GEN
     263      513156 : subtract_slices(long m, long n,
     264             :                 GEN A, long ma, long da, long na, long ea,
     265             :                 GEN B, long mb, long db, long nb, long eb)
     266             : {
     267      513156 :   long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
     268      513156 :   GEN M = cgetg(n + 1, t_MAT), C;
     269             : 
     270     4327305 :   for (j = 1; j <= min_e; j++) {
     271     3814149 :     gel(M, j) = C = cgetg(m + 1, t_COL);
     272    75909498 :     for (i = 1; i <= min_d; i++)
     273    72095349 :       gel(C, i) = subii(gcoeff(A, ma + i, na + j),
     274    72095349 :                         gcoeff(B, mb + i, nb + j));
     275     3871250 :     for (; i <= da; i++) gel(C, i) = gcoeff(A, ma + i, na + j);
     276     3902202 :     for (; i <= db; i++) gel(C, i) = negi(gcoeff(B, mb + i, nb + j));
     277     3814149 :     for (; i <= m; i++) gel(C, i) = gen_0;
     278             :   }
     279      513156 :   for (; j <= ea; j++) {
     280           0 :     gel(M, j) = C = cgetg(m + 1, t_COL);
     281           0 :     for (i = 1; i <= da; i++) gel(C, i) = gcoeff(A, ma + i, na + j);
     282           0 :     for (; i <= m; i++) gel(C, i) = gen_0;
     283             :   }
     284      548145 :   for (; j <= eb; j++) {
     285       34989 :     gel(M, j) = C = cgetg(m + 1, t_COL);
     286      127034 :     for (i = 1; i <= db; i++) gel(C, i) = negi(gcoeff(B, mb + i, nb + j));
     287       34989 :     for (; i <= m; i++) gel(C, i) = gen_0;
     288             :   }
     289      548145 :   for (; j <= n; j++) gel(M, j) = zerocol(m);
     290      513156 :   return M;
     291             : }
     292             : 
     293             : static GEN ZM_mul_i(GEN x, GEN y, long l, long lx, long ly);
     294             : 
     295             : /* Strassen-Winograd matrix product A (m x n) * B (n x p) */
     296             : static GEN
     297       73308 : ZM_mul_sw(GEN A, GEN B, long m, long n, long p)
     298             : {
     299       73308 :   pari_sp av = avma;
     300       73308 :   long m1 = (m + 1)/2, m2 = m/2,
     301       73308 :     n1 = (n + 1)/2, n2 = n/2,
     302       73308 :     p1 = (p + 1)/2, p2 = p/2;
     303             :   GEN A11, A12, A22, B11, B21, B22,
     304             :     S1, S2, S3, S4, T1, T2, T3, T4,
     305             :     M1, M2, M3, M4, M5, M6, M7,
     306             :     V1, V2, V3, C11, C12, C21, C22, C;
     307             : 
     308       73308 :   T2 = subtract_slices(n1, p2, B, 0, n1, p1, p2, B, n1, n2, p1, p2);
     309       73308 :   S1 = subtract_slices(m2, n1, A, m1, m2, 0, n1, A, 0, m2, 0, n1);
     310       73308 :   M2 = ZM_mul_i(S1, T2, m2 + 1, n1 + 1, p2 + 1);
     311       73308 :   if (gc_needed(av, 1))
     312           0 :     gerepileall(av, 2, &T2, &M2);  /* destroy S1 */
     313       73308 :   T3 = subtract_slices(n1, p1, T2, 0, n1, 0, p2, B, 0, n1, 0, p1);
     314       73308 :   if (gc_needed(av, 1))
     315           0 :     gerepileall(av, 2, &M2, &T3);  /* destroy T2 */
     316       73308 :   S2 = add_slices(m2, n1, A, m1, m2, 0, n1, A, m1, m2, n1, n2);
     317       73308 :   T1 = subtract_slices(n1, p1, B, 0, n1, p1, p2, B, 0, n1, 0, p2);
     318       73308 :   M3 = ZM_mul_i(S2, T1, m2 + 1, n1 + 1, p2 + 1);
     319       73308 :   if (gc_needed(av, 1))
     320           0 :     gerepileall(av, 4, &M2, &T3, &S2, &M3);  /* destroy T1 */
     321       73308 :   S3 = subtract_slices(m1, n1, S2, 0, m2, 0, n1, A, 0, m1, 0, n1);
     322       73308 :   if (gc_needed(av, 1))
     323           0 :     gerepileall(av, 4, &M2, &T3, &M3, &S3);  /* destroy S2 */
     324       73308 :   A11 = matslice(A, 1, m1, 1, n1);
     325       73308 :   B11 = matslice(B, 1, n1, 1, p1);
     326       73308 :   M1 = ZM_mul_i(A11, B11, m1 + 1, n1 + 1, p1 + 1);
     327       73308 :   if (gc_needed(av, 1))
     328           0 :     gerepileall(av, 5, &M2, &T3, &M3, &S3, &M1);  /* destroy A11, B11 */
     329       73308 :   A12 = matslice(A, 1, m1, n1 + 1, n);
     330       73308 :   B21 = matslice(B, n1 + 1, n, 1, p1);
     331       73308 :   M4 = ZM_mul_i(A12, B21, m1 + 1, n2 + 1, p1 + 1);
     332       73308 :   if (gc_needed(av, 1))
     333           0 :     gerepileall(av, 6, &M2, &T3, &M3, &S3, &M1, &M4);  /* destroy A12, B21 */
     334       73308 :   C11 = add_slices(m1, p1, M1, 0, m1, 0, p1, M4, 0, m1, 0, p1);
     335       73308 :   if (gc_needed(av, 1))
     336           0 :     gerepileall(av, 6, &M2, &T3, &M3, &S3, &M1, &C11);  /* destroy M4 */
     337       73308 :   M5 = ZM_mul_i(S3, T3, m1 + 1, n1 + 1, p1 + 1);
     338       73308 :   S4 = subtract_slices(m1, n2, A, 0, m1, n1, n2, S3, 0, m1, 0, n2);
     339       73308 :   if (gc_needed(av, 1))
     340           5 :     gerepileall(av, 7, &M2, &T3, &M3, &M1, &C11, &M5, &S4);  /* destroy S3 */
     341       73308 :   T4 = add_slices(n2, p1, B, n1, n2, 0, p1, T3, 0, n2, 0, p1);
     342       73308 :   if (gc_needed(av, 1))
     343           0 :     gerepileall(av, 7, &M2, &M3, &M1, &C11, &M5, &S4, &T4);  /* destroy T3 */
     344       73308 :   V1 = subtract_slices(m1, p1, M1, 0, m1, 0, p1, M5, 0, m1, 0, p1);
     345       73308 :   if (gc_needed(av, 1))
     346           1 :     gerepileall(av, 6, &M2, &M3, &S4, &T4, &C11, &V1);  /* destroy M1, M5 */
     347       73308 :   B22 = matslice(B, n1 + 1, n, p1 + 1, p);
     348       73308 :   M6 = ZM_mul_i(S4, B22, m1 + 1, n2 + 1, p2 + 1);
     349       73308 :   if (gc_needed(av, 1))
     350           6 :     gerepileall(av, 6, &M2, &M3, &T4, &C11, &V1, &M6);  /* destroy S4, B22 */
     351       73308 :   A22 = matslice(A, m1 + 1, m, n1 + 1, n);
     352       73308 :   M7 = ZM_mul_i(A22, T4, m2 + 1, n2 + 1, p1 + 1);
     353       73308 :   if (gc_needed(av, 1))
     354           0 :     gerepileall(av, 6, &M2, &M3, &C11, &V1, &M6, &M7);  /* destroy A22, T4 */
     355       73308 :   V3 = add_slices(m1, p2, V1, 0, m1, 0, p2, M3, 0, m2, 0, p2);
     356       73308 :   C12 = add_slices(m1, p2, V3, 0, m1, 0, p2, M6, 0, m1, 0, p2);
     357       73308 :   if (gc_needed(av, 1))
     358           6 :     gerepileall(av, 6, &M2, &M3, &C11, &V1, &M7, &C12);  /* destroy V3, M6 */
     359       73308 :   V2 = add_slices(m2, p1, V1, 0, m2, 0, p1, M2, 0, m2, 0, p2);
     360       73308 :   if (gc_needed(av, 1))
     361           0 :     gerepileall(av, 5, &M3, &C11, &M7, &C12, &V2);  /* destroy V1, M2 */
     362       73308 :   C21 = add_slices(m2, p1, V2, 0, m2, 0, p1, M7, 0, m2, 0, p1);
     363       73308 :   if (gc_needed(av, 1))
     364           6 :     gerepileall(av, 5, &M3, &C11, &C12, &V2, &C21);  /* destroy M7 */
     365       73308 :   C22 = add_slices(m2, p2, V2, 0, m2, 0, p2, M3, 0, m2, 0, p2);
     366       73308 :   if (gc_needed(av, 1))
     367           0 :     gerepileall(av, 4, &C11, &C12, &C21, &C22);  /* destroy V2, M3 */
     368       73308 :   C = shallowconcat(vconcat(C11, C21), vconcat(C12, C22));
     369       73308 :   return gerepilecopy(av, C);
     370             : }
     371             : 
     372             : /* x[i,]*y. Assume lg(x) > 1 and 0 < i < lgcols(x) */
     373             : static GEN
     374   585524351 : ZMrow_ZC_mul_i(GEN x, GEN y, long i, long lx)
     375             : {
     376   585524351 :   pari_sp av = avma;
     377   585524351 :   GEN c = mulii(gcoeff(x,i,1), gel(y,1)), ZERO = gen_0;
     378             :   long k;
     379  6691624526 :   for (k = 2; k < lx; k++)
     380             :   {
     381  6106476449 :     GEN t = mulii(gcoeff(x,i,k), gel(y,k));
     382  6105059499 :     if (t != ZERO) c = addii(c, t);
     383             :   }
     384   585148077 :   return gerepileuptoint(av, c);
     385             : }
     386             : GEN
     387   135063066 : ZMrow_ZC_mul(GEN x, GEN y, long i)
     388   135063066 : { return ZMrow_ZC_mul_i(x, y, i, lg(x)); }
     389             : 
     390             : /* return x * y, 1 < lx = lg(x), l = lgcols(x) */
     391             : static GEN
     392    77655127 : ZM_ZC_mul_i(GEN x, GEN y, long lx, long l)
     393             : {
     394    77655127 :   GEN z = cgetg(l,t_COL);
     395             :   long i;
     396   528140456 :   for (i=1; i<l; i++) gel(z,i) = ZMrow_ZC_mul_i(x,y,i,lx);
     397    77614479 :   return z;
     398             : }
     399             : 
     400             : static GEN
     401    14195263 : ZM_mul_classical(GEN x, GEN y, long l, long lx, long ly)
     402             : {
     403             :   long j;
     404    14195263 :   GEN z = cgetg(ly, t_MAT);
     405    67755609 :   for (j = 1; j < ly; j++)
     406    53563267 :     gel(z, j) = ZM_ZC_mul_i(x, gel(y, j), lx, l);
     407    14192342 :   return z;
     408             : }
     409             : 
     410             : static GEN
     411        1120 : ZM_mul_slice(GEN A, GEN B, GEN P, GEN *mod)
     412             : {
     413        1120 :   pari_sp av = avma;
     414        1120 :   long i, n = lg(P)-1;
     415             :   GEN H, T;
     416        1120 :   if (n == 1)
     417             :   {
     418           0 :     ulong p = uel(P,1);
     419           0 :     GEN a = ZM_to_Flm(A, p);
     420           0 :     GEN b = ZM_to_Flm(B, p);
     421           0 :     GEN Hp = gerepileupto(av, Flm_to_ZM(Flm_mul(a, b, p)));
     422           0 :     *mod = utoi(p); return Hp;
     423             :   }
     424        1120 :   T = ZV_producttree(P);
     425        1120 :   A = ZM_nv_mod_tree(A, P, T);
     426        1120 :   B = ZM_nv_mod_tree(B, P, T);
     427        1120 :   H = cgetg(n+1, t_VEC);
     428        7938 :   for(i=1; i <= n; i++)
     429        6818 :     gel(H,i) = Flm_mul(gel(A,i),gel(B,i),P[i]);
     430        1120 :   H = nmV_chinese_center_tree_seq(H, P, T, ZV_chinesetree(P, T));
     431        1120 :   *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
     432             : }
     433             : 
     434             : GEN
     435        1120 : ZM_mul_worker(GEN P, GEN A, GEN B)
     436             : {
     437        1120 :   GEN V = cgetg(3, t_VEC);
     438        1120 :   gel(V,1) = ZM_mul_slice(A, B, P, &gel(V,2));
     439        1120 :   return V;
     440             : }
     441             : 
     442             : static GEN
     443         854 : ZM_mul_fast(GEN A, GEN B, long lA, long lB, long sA, long sB)
     444             : {
     445         854 :   pari_sp av = avma;
     446             :   forprime_t S;
     447             :   GEN H, worker;
     448             :   long h;
     449         854 :   if (sA == 2 || sB == 2) return zeromat(nbrows(A),lB-1);
     450         842 :   h = 1 + (sA + sB - 4) * BITS_IN_LONG + expu(lA-1);
     451         842 :   init_modular_big(&S);
     452         842 :   worker = snm_closure(is_entry("_ZM_mul_worker"), mkvec2(A,B));
     453         842 :   H = gen_crt("ZM_mul", worker, &S, NULL, h, 0, NULL,
     454             :               nmV_chinese_center, FpM_center);
     455         842 :   return gerepileupto(av, H);
     456             : }
     457             : 
     458             : /* s = min(log_BIL |x|, log_BIL |y|), use Strassen-Winograd when
     459             :  * min(dims) > B */
     460             : static long
     461    14268571 : sw_bound(long s)
     462    14268571 : { return s > 60 ? 2: s > 25 ? 4: s > 15 ? 8 : s > 8 ? 16 : 32; }
     463             : 
     464             : /* assume lx > 1 and ly > 1; x is (l-1) x (lx-1), y is (lx-1) x (ly-1).
     465             :  * Return x * y */
     466             : static GEN
     467    21573561 : ZM_mul_i(GEN x, GEN y, long l, long lx, long ly)
     468             : {
     469             :   long sx, sy, B;
     470             : #ifdef LONG_IS_64BIT /* From Flm_mul_i */
     471    18318957 :   long Flm_sw_bound = 70;
     472             : #else
     473     3254604 :   long Flm_sw_bound = 120;
     474             : #endif
     475    21573561 :   if (l == 1) return zeromat(0, ly-1);
     476    21571545 :   if (lx==2 && l==2 && ly==2)
     477      359209 :   { retmkmat(mkcol(mulii(gcoeff(x,1,1), gcoeff(y,1,1)))); }
     478    21212336 :   if (lx==3 && l==3 && ly==3) return ZM2_mul(x, y);
     479    14245238 :   sx = ZM_max_lg_i(x, lx, l);
     480    14246009 :   sy = ZM_max_lg_i(y, ly, lx);
     481             :   /* Use modular reconstruction if Flm_mul would use Strassen and the input
     482             :    * sizes look balanced */
     483    14246005 :   if (lx > Flm_sw_bound && ly > Flm_sw_bound && l > Flm_sw_bound
     484         866 :       && sx <= 10 * sy && sy <= 10 * sx) return ZM_mul_fast(x,y, lx,ly, sx,sy);
     485             : 
     486    14245151 :   B = sw_bound(minss(sx, sy));
     487    14245126 :   if (l <= B || lx <= B || ly <= B)
     488    14171848 :     return ZM_mul_classical(x, y, l, lx, ly);
     489             :   else
     490       73278 :     return ZM_mul_sw(x, y, l - 1, lx - 1, ly - 1);
     491             : }
     492             : 
     493             : GEN
     494    21197036 : ZM_mul(GEN x, GEN y)
     495             : {
     496    21197036 :   long lx = lg(x), ly = lg(y);
     497    21197036 :   if (ly == 1) return cgetg(1,t_MAT);
     498    21061947 :   if (lx == 1) return zeromat(0, ly-1);
     499    21060365 :   return ZM_mul_i(x, y, lgcols(x), lx, ly);
     500             : }
     501             : 
     502             : static GEN
     503           0 : ZM_sqr_slice(GEN A, GEN P, GEN *mod)
     504             : {
     505           0 :   pari_sp av = avma;
     506           0 :   long i, n = lg(P)-1;
     507             :   GEN H, T;
     508           0 :   if (n == 1)
     509             :   {
     510           0 :     ulong p = uel(P,1);
     511           0 :     GEN a = ZM_to_Flm(A, p);
     512           0 :     GEN Hp = gerepileupto(av, Flm_to_ZM(Flm_sqr(a, p)));
     513           0 :     *mod = utoi(p); return Hp;
     514             :   }
     515           0 :   T = ZV_producttree(P);
     516           0 :   A = ZM_nv_mod_tree(A, P, T);
     517           0 :   H = cgetg(n+1, t_VEC);
     518           0 :   for(i=1; i <= n; i++)
     519           0 :     gel(H,i) = Flm_sqr(gel(A,i), P[i]);
     520           0 :   H = nmV_chinese_center_tree_seq(H, P, T, ZV_chinesetree(P, T));
     521           0 :   *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
     522             : }
     523             : 
     524             : GEN
     525           0 : ZM_sqr_worker(GEN P, GEN A)
     526             : {
     527           0 :   GEN V = cgetg(3, t_VEC);
     528           0 :   gel(V,1) = ZM_sqr_slice(A, P, &gel(V,2));
     529           0 :   return V;
     530             : }
     531             : 
     532             : static GEN
     533           0 : ZM_sqr_fast(GEN A, long l, long s)
     534             : {
     535           0 :   pari_sp av = avma;
     536             :   forprime_t S;
     537             :   GEN H, worker;
     538             :   long h;
     539           0 :   if (s == 2) return zeromat(l-1,l-1);
     540           0 :   h = 1 + (2*s - 4) * BITS_IN_LONG + expu(l-1);
     541           0 :   init_modular_big(&S);
     542           0 :   worker = snm_closure(is_entry("_ZM_sqr_worker"), mkvec(A));
     543           0 :   H = gen_crt("ZM_sqr", worker, &S, NULL, h, 0, NULL,
     544             :               nmV_chinese_center, FpM_center);
     545           0 :   return gerepileupto(av, H);
     546             : }
     547             : 
     548             : GEN
     549      553869 : QM_mul(GEN x, GEN y)
     550             : {
     551      553869 :   GEN dx, nx = Q_primitive_part(x, &dx);
     552      553869 :   GEN dy, ny = Q_primitive_part(y, &dy);
     553      553869 :   GEN z = ZM_mul(nx, ny);
     554      553869 :   if (dx || dy)
     555             :   {
     556      469754 :     GEN d = dx ? dy ? gmul(dx, dy): dx : dy;
     557      469754 :     if (!gequal1(d)) z = ZM_Q_mul(z, d);
     558             :   }
     559      553869 :   return z;
     560             : }
     561             : 
     562             : GEN
     563         700 : QM_sqr(GEN x)
     564             : {
     565         700 :   GEN dx, nx = Q_primitive_part(x, &dx);
     566         700 :   GEN z = ZM_sqr(nx);
     567         700 :   if (dx)
     568         700 :     z = ZM_Q_mul(z, gsqr(dx));
     569         700 :   return z;
     570             : }
     571             : 
     572             : GEN
     573      143013 : QM_QC_mul(GEN x, GEN y)
     574             : {
     575      143013 :   GEN dx, nx = Q_primitive_part(x, &dx);
     576      143013 :   GEN dy, ny = Q_primitive_part(y, &dy);
     577      143013 :   GEN z = ZM_ZC_mul(nx, ny);
     578      143013 :   if (dx || dy)
     579             :   {
     580      142985 :     GEN d = dx ? dy ? gmul(dx, dy): dx : dy;
     581      142985 :     if (!gequal1(d)) z = ZC_Q_mul(z, d);
     582             :   }
     583      143013 :   return z;
     584             : }
     585             : 
     586             : /* assume result is symmetric */
     587             : GEN
     588           0 : ZM_multosym(GEN x, GEN y)
     589             : {
     590           0 :   long j, lx, ly = lg(y);
     591             :   GEN M;
     592           0 :   if (ly == 1) return cgetg(1,t_MAT);
     593           0 :   lx = lg(x); /* = lgcols(y) */
     594           0 :   if (lx == 1) return cgetg(1,t_MAT);
     595             :   /* ly = lgcols(x) */
     596           0 :   M = cgetg(ly, t_MAT);
     597           0 :   for (j=1; j<ly; j++)
     598             :   {
     599           0 :     GEN z = cgetg(ly,t_COL), yj = gel(y,j);
     600             :     long i;
     601           0 :     for (i=1; i<j; i++) gel(z,i) = gcoeff(M,j,i);
     602           0 :     for (i=j; i<ly; i++)gel(z,i) = ZMrow_ZC_mul_i(x,yj,i,lx);
     603           0 :     gel(M,j) = z;
     604             :   }
     605           0 :   return M;
     606             : }
     607             : 
     608             : /* compute m*diagonal(d), assume lg(d) = lg(m). Shallow */
     609             : GEN
     610          21 : ZM_mul_diag(GEN m, GEN d)
     611             : {
     612             :   long j, l;
     613          21 :   GEN y = cgetg_copy(m, &l);
     614          77 :   for (j=1; j<l; j++)
     615             :   {
     616          56 :     GEN c = gel(d,j);
     617          56 :     gel(y,j) = equali1(c)? gel(m,j): ZC_Z_mul(gel(m,j), c);
     618             :   }
     619          21 :   return y;
     620             : }
     621             : /* compute diagonal(d)*m, assume lg(d) = lg(m~). Shallow */
     622             : GEN
     623      533867 : ZM_diag_mul(GEN d, GEN m)
     624             : {
     625      533867 :   long i, j, l = lg(d), lm = lg(m);
     626      533867 :   GEN y = cgetg(lm, t_MAT);
     627     1521174 :   for (j=1; j<lm; j++) gel(y,j) = cgetg(l, t_COL);
     628     1647597 :   for (i=1; i<l; i++)
     629             :   {
     630     1114065 :     GEN c = gel(d,i);
     631     1114065 :     if (equali1(c))
     632      236605 :       for (j=1; j<lm; j++) gcoeff(y,i,j) = gcoeff(m,i,j);
     633             :     else
     634     3399719 :       for (j=1; j<lm; j++) gcoeff(y,i,j) = mulii(gcoeff(m,i,j), c);
     635             :   }
     636      533532 :   return y;
     637             : }
     638             : 
     639             : /* assume lx > 1 is lg(x) = lg(y) */
     640             : static GEN
     641    19174702 : ZV_dotproduct_i(GEN x, GEN y, long lx)
     642             : {
     643    19174702 :   pari_sp av = avma;
     644    19174702 :   GEN c = mulii(gel(x,1), gel(y,1));
     645             :   long i;
     646   144134148 :   for (i = 2; i < lx; i++)
     647             :   {
     648   124960842 :     GEN t = mulii(gel(x,i), gel(y,i));
     649   124960237 :     if (t != gen_0) c = addii(c, t);
     650             :   }
     651    19173306 :   return gerepileuptoint(av, c);
     652             : }
     653             : 
     654             : /* x~ * y, assuming result is symmetric */
     655             : GEN
     656      532221 : ZM_transmultosym(GEN x, GEN y)
     657             : {
     658      532221 :   long i, j, l, ly = lg(y);
     659             :   GEN M;
     660      532221 :   if (ly == 1) return cgetg(1,t_MAT);
     661             :   /* lg(x) = ly */
     662      532221 :   l = lgcols(y); /* = lgcols(x) */
     663      532221 :   M = cgetg(ly, t_MAT);
     664     2708730 :   for (i=1; i<ly; i++)
     665             :   {
     666     2176509 :     GEN xi = gel(x,i), c = cgetg(ly,t_COL);
     667     2176509 :     gel(M,i) = c;
     668     7109000 :     for (j=1; j<i; j++)
     669     4932491 :       gcoeff(M,i,j) = gel(c,j) = ZV_dotproduct_i(xi,gel(y,j),l);
     670     2176509 :     gel(c,i) = ZV_dotproduct_i(xi,gel(y,i),l);
     671             :   }
     672      532221 :   return M;
     673             : }
     674             : /* x~ * y */
     675             : GEN
     676        2289 : ZM_transmul(GEN x, GEN y)
     677             : {
     678        2289 :   long i, j, l, lx, ly = lg(y);
     679             :   GEN M;
     680        2289 :   if (ly == 1) return cgetg(1,t_MAT);
     681        2289 :   lx = lg(x);
     682        2289 :   l = lgcols(y);
     683        2289 :   if (lgcols(x) != l) pari_err_OP("operation 'ZM_transmul'", x,y);
     684        2289 :   M = cgetg(ly, t_MAT);
     685        6993 :   for (i=1; i<ly; i++)
     686             :   {
     687        4704 :     GEN yi = gel(y,i), c = cgetg(lx,t_COL);
     688        4704 :     gel(M,i) = c;
     689       12229 :     for (j=1; j<lx; j++) gel(c,j) = ZV_dotproduct_i(yi,gel(x,j),l);
     690             :   }
     691        2289 :   return M;
     692             : }
     693             : 
     694             : /* assume l > 1; x is (l-1) x (l-1), return x^2.
     695             :  * FIXME: we ultimately rely on Strassen-Winograd which uses 7M + 15A.
     696             :  * Should use Bodrato's variant of Winograd, using 3M + 4S + 11A */
     697             : static GEN
     698       24213 : ZM_sqr_i(GEN x, long l)
     699             : {
     700             :   long s;
     701       24213 :   if (l == 2) { retmkmat(mkcol(sqri(gcoeff(x,1,1)))); }
     702       24213 :   if (l == 3) return ZM2_sqr(x);
     703       23436 :   s = ZM_max_lg_i(x, l, l);
     704       23436 :   if (l > 70) return ZM_sqr_fast(x, l, s);
     705       23436 :   if (l <= sw_bound(s))
     706       23406 :     return ZM_mul_classical(x, x, l, l, l);
     707             :   else
     708          30 :     return ZM_mul_sw(x, x, l - 1, l - 1, l - 1);
     709             : }
     710             : 
     711             : GEN
     712       24213 : ZM_sqr(GEN x)
     713             : {
     714       24213 :   long lx=lg(x);
     715       24213 :   if (lx==1) return cgetg(1,t_MAT);
     716       24213 :   return ZM_sqr_i(x, lx);
     717             : }
     718             : GEN
     719    24181454 : ZM_ZC_mul(GEN x, GEN y)
     720             : {
     721    24181454 :   long lx = lg(x);
     722    24181454 :   return lx==1? cgetg(1,t_COL): ZM_ZC_mul_i(x, y, lx, lgcols(x));
     723             : }
     724             : 
     725             : GEN
     726     4527374 : ZC_Z_div(GEN x, GEN c)
     727    20035228 : { pari_APPLY_type(t_COL, Qdivii(gel(x,i), c)) }
     728             : 
     729             : GEN
     730       17553 : ZM_Z_div(GEN x, GEN c)
     731      202068 : { pari_APPLY_same(ZC_Z_div(gel(x, i), c)) }
     732             : 
     733             : GEN
     734     2723355 : ZC_Q_mul(GEN A, GEN z)
     735             : {
     736     2723355 :   pari_sp av = avma;
     737     2723355 :   long i, l = lg(A);
     738             :   GEN d, n, Ad, B, u;
     739     2723355 :   if (typ(z)==t_INT) return ZC_Z_mul(A,z);
     740     2717097 :   n = gel(z, 1); d = gel(z, 2);
     741     2717097 :   Ad = FpC_red(A, d);
     742     2717051 :   u = gcdii(d, FpV_factorback(Ad, NULL, d));
     743     2717027 :   B = cgetg(l, t_COL);
     744     2717026 :   if (equali1(u))
     745             :   {
     746      421961 :     for(i=1; i<l; i++)
     747      355443 :       gel(B, i) = mkfrac(mulii(n, gel(A,i)), d);
     748             :   } else
     749             :   {
     750    18462431 :     for(i=1; i<l; i++)
     751             :     {
     752    15812033 :       GEN di = gcdii(gel(Ad, i), u), ni = mulii(n, diviiexact(gel(A,i), di));
     753    15811924 :       if (equalii(d, di))
     754    10849952 :         gel(B, i) = ni;
     755             :       else
     756     4961953 :         gel(B, i) = mkfrac(ni, diviiexact(d, di));
     757             :     }
     758             :   }
     759     2716916 :   return gerepilecopy(av, B);
     760             : }
     761             : 
     762             : GEN
     763     1102011 : ZM_Q_mul(GEN x, GEN z)
     764             : {
     765     1102011 :   if (typ(z)==t_INT) return ZM_Z_mul(x,z);
     766     3245370 :   pari_APPLY_same(ZC_Q_mul(gel(x, i), z));
     767             : }
     768             : 
     769             : long
     770   196399119 : zv_dotproduct(GEN x, GEN y)
     771             : {
     772   196399119 :   long i, lx = lg(x);
     773             :   ulong c;
     774   196399119 :   if (lx == 1) return 0;
     775   196399119 :   c = uel(x,1)*uel(y,1);
     776  3047359154 :   for (i = 2; i < lx; i++)
     777  2850960035 :     c += uel(x,i)*uel(y,i);
     778   196399119 :   return c;
     779             : }
     780             : 
     781             : GEN
     782      230629 : ZV_ZM_mul(GEN x, GEN y)
     783             : {
     784      230629 :   long i, lx = lg(x), ly = lg(y);
     785             :   GEN z;
     786      230629 :   if (lx == 1) return zerovec(ly-1);
     787      230510 :   z = cgetg(ly, t_VEC);
     788      882063 :   for (i = 1; i < ly; i++) gel(z,i) = ZV_dotproduct_i(x, gel(y,i), lx);
     789      230510 :   return z;
     790             : }
     791             : 
     792             : GEN
     793           0 : ZC_ZV_mul(GEN x, GEN y)
     794             : {
     795           0 :   long i,j, lx=lg(x), ly=lg(y);
     796             :   GEN z;
     797           0 :   if (ly==1) return cgetg(1,t_MAT);
     798           0 :   z = cgetg(ly,t_MAT);
     799           0 :   for (j=1; j < ly; j++)
     800             :   {
     801           0 :     gel(z,j) = cgetg(lx,t_COL);
     802           0 :     for (i=1; i<lx; i++) gcoeff(z,i,j) = mulii(gel(x,i),gel(y,j));
     803             :   }
     804           0 :   return z;
     805             : }
     806             : 
     807             : GEN
     808     6659916 : ZV_dotsquare(GEN x)
     809             : {
     810             :   long i, lx;
     811             :   pari_sp av;
     812             :   GEN z;
     813     6659916 :   lx = lg(x);
     814     6659916 :   if (lx == 1) return gen_0;
     815     6659916 :   av = avma; z = sqri(gel(x,1));
     816    26354818 :   for (i=2; i<lx; i++) z = addii(z, sqri(gel(x,i)));
     817     6658234 :   return gerepileuptoint(av,z);
     818             : }
     819             : 
     820             : GEN
     821    16185758 : ZV_dotproduct(GEN x,GEN y)
     822             : {
     823             :   long lx;
     824    16185758 :   if (x == y) return ZV_dotsquare(x);
     825    11405704 :   lx = lg(x);
     826    11405704 :   if (lx == 1) return gen_0;
     827    11405704 :   return ZV_dotproduct_i(x, y, lx);
     828             : }
     829             : 
     830             : static GEN
     831         280 : _ZM_mul(void *data /*ignored*/, GEN x, GEN y)
     832         280 : { (void)data; return ZM_mul(x,y); }
     833             : static GEN
     834       23197 : _ZM_sqr(void *data /*ignored*/, GEN x)
     835       23197 : { (void)data; return ZM_sqr(x); }
     836             : /* FIXME: Using Bodrato's squaring, precomputations attached to fixed
     837             :  * multiplicand should be reused. And some postcomputations can be fused */
     838             : GEN
     839           0 : ZM_pow(GEN x, GEN n)
     840             : {
     841           0 :   pari_sp av = avma;
     842           0 :   if (!signe(n)) return matid(lg(x)-1);
     843           0 :   return gerepilecopy(av, gen_pow_i(x, n, NULL, &_ZM_sqr, &_ZM_mul));
     844             : }
     845             : GEN
     846       22672 : ZM_powu(GEN x, ulong n)
     847             : {
     848       22672 :   pari_sp av = avma;
     849       22672 :   if (!n) return matid(lg(x)-1);
     850       22672 :   return gerepilecopy(av, gen_powu_i(x, n, NULL, &_ZM_sqr, &_ZM_mul));
     851             : }
     852             : /********************************************************************/
     853             : /**                                                                **/
     854             : /**                           ADD, SUB                             **/
     855             : /**                                                                **/
     856             : /********************************************************************/
     857             : static GEN
     858    36280730 : ZC_add_i(GEN x, GEN y, long lx)
     859             : {
     860    36280730 :   GEN A = cgetg(lx, t_COL);
     861             :   long i;
     862   491322598 :   for (i=1; i<lx; i++) gel(A,i) = addii(gel(x,i), gel(y,i));
     863    36273500 :   return A;
     864             : }
     865             : GEN
     866    26483217 : ZC_add(GEN x, GEN y) { return ZC_add_i(x, y, lg(x)); }
     867             : GEN
     868      364883 : ZC_Z_add(GEN x, GEN y)
     869             : {
     870      364883 :   long k, lx = lg(x);
     871      364883 :   GEN z = cgetg(lx, t_COL);
     872      364880 :   if (lx == 1) pari_err_TYPE2("+",x,y);
     873      364880 :   gel(z,1) = addii(y,gel(x,1));
     874     2511925 :   for (k = 2; k < lx; k++) gel(z,k) = icopy(gel(x,k));
     875      364871 :   return z;
     876             : }
     877             : 
     878             : static GEN
     879     9054355 : ZC_sub_i(GEN x, GEN y, long lx)
     880             : {
     881             :   long i;
     882     9054355 :   GEN A = cgetg(lx, t_COL);
     883    64422581 :   for (i=1; i<lx; i++) gel(A,i) = subii(gel(x,i), gel(y,i));
     884     9053771 :   return A;
     885             : }
     886             : GEN
     887     8754741 : ZC_sub(GEN x, GEN y) { return ZC_sub_i(x, y, lg(x)); }
     888             : GEN
     889           0 : ZC_Z_sub(GEN x, GEN y)
     890             : {
     891           0 :   long k, lx = lg(x);
     892           0 :   GEN z = cgetg(lx, t_COL);
     893           0 :   if (lx == 1) pari_err_TYPE2("+",x,y);
     894           0 :   gel(z,1) = subii(gel(x,1), y);
     895           0 :   for (k = 2; k < lx; k++) gel(z,k) = icopy(gel(x,k));
     896           0 :   return z;
     897             : }
     898             : GEN
     899      545801 : Z_ZC_sub(GEN a, GEN x)
     900             : {
     901      545801 :   long k, lx = lg(x);
     902      545801 :   GEN z = cgetg(lx, t_COL);
     903      545802 :   if (lx == 1) pari_err_TYPE2("-",a,x);
     904      545802 :   gel(z,1) = subii(a, gel(x,1));
     905     1511918 :   for (k = 2; k < lx; k++) gel(z,k) = negi(gel(x,k));
     906      545784 :   return z;
     907             : }
     908             : 
     909             : GEN
     910      855046 : ZM_add(GEN x, GEN y)
     911             : {
     912      855046 :   long lx = lg(x), l, j;
     913             :   GEN z;
     914      855046 :   if (lx == 1) return cgetg(1, t_MAT);
     915      775527 :   z = cgetg(lx, t_MAT); l = lgcols(x);
     916    10572816 :   for (j = 1; j < lx; j++) gel(z,j) = ZC_add_i(gel(x,j), gel(y,j), l);
     917      775531 :   return z;
     918             : }
     919             : GEN
     920       91596 : ZM_sub(GEN x, GEN y)
     921             : {
     922       91596 :   long lx = lg(x), l, j;
     923             :   GEN z;
     924       91596 :   if (lx == 1) return cgetg(1, t_MAT);
     925       91596 :   z = cgetg(lx, t_MAT); l = lgcols(x);
     926      391155 :   for (j = 1; j < lx; j++) gel(z,j) = ZC_sub_i(gel(x,j), gel(y,j), l);
     927       91596 :   return z;
     928             : }
     929             : /********************************************************************/
     930             : /**                                                                **/
     931             : /**                         LINEAR COMBINATION                     **/
     932             : /**                                                                **/
     933             : /********************************************************************/
     934             : /* return X/c assuming division is exact */
     935             : GEN
     936     5134828 : ZC_Z_divexact(GEN x, GEN c)
     937    72257863 : { pari_APPLY_type(t_COL, diviiexact(gel(x,i), c)) }
     938             : GEN
     939        2261 : ZC_divexactu(GEN x, ulong c)
     940       11375 : { pari_APPLY_type(t_COL, diviuexact(gel(x,i), c)) }
     941             : 
     942             : GEN
     943      461514 : ZM_Z_divexact(GEN x, GEN c)
     944     3527721 : { pari_APPLY_same(ZC_Z_divexact(gel(x,i), c)) }
     945             : 
     946             : GEN
     947         441 : ZM_divexactu(GEN x, ulong c)
     948        2702 : { pari_APPLY_same(ZC_divexactu(gel(x,i), c)) }
     949             : 
     950             : GEN
     951    35528556 : ZC_Z_mul(GEN x, GEN c)
     952             : {
     953    35528556 :   if (!signe(c)) return zerocol(lg(x)-1);
     954    34160857 :   if (is_pm1(c)) return (signe(c) > 0)? ZC_copy(x): ZC_neg(x);
     955   255707667 :   pari_APPLY_type(t_COL, mulii(gel(x,i), c))
     956             : }
     957             : 
     958             : GEN
     959       61925 : ZC_z_mul(GEN x, long c)
     960             : {
     961       61925 :   if (!c) return zerocol(lg(x)-1);
     962       52766 :   if (c == 1) return ZC_copy(x);
     963       47912 :   if (c ==-1) return ZC_neg(x);
     964      478926 :   pari_APPLY_type(t_COL, mulsi(c, gel(x,i)))
     965             : }
     966             : 
     967             : GEN
     968       28498 : zv_z_mul(GEN x, long n)
     969      432192 : { pari_APPLY_long(x[i]*n) }
     970             : 
     971             : /* return a ZM */
     972             : GEN
     973           0 : nm_Z_mul(GEN X, GEN c)
     974             : {
     975           0 :   long i, j, h, l = lg(X), s = signe(c);
     976             :   GEN A;
     977           0 :   if (l == 1) return cgetg(1, t_MAT);
     978           0 :   h = lgcols(X);
     979           0 :   if (!s) return zeromat(h-1, l-1);
     980           0 :   if (is_pm1(c)) {
     981           0 :     if (s > 0) return Flm_to_ZM(X);
     982           0 :     X = Flm_to_ZM(X); ZM_togglesign(X); return X;
     983             :   }
     984           0 :   A = cgetg(l, t_MAT);
     985           0 :   for (j = 1; j < l; j++)
     986             :   {
     987           0 :     GEN a = cgetg(h, t_COL), x = gel(X, j);
     988           0 :     for (i = 1; i < h; i++) gel(a,i) = muliu(c, x[i]);
     989           0 :     gel(A,j) = a;
     990             :   }
     991           0 :   return A;
     992             : }
     993             : GEN
     994     3210912 : ZM_Z_mul(GEN X, GEN c)
     995             : {
     996     3210912 :   long i, j, h, l = lg(X);
     997             :   GEN A;
     998     3210912 :   if (l == 1) return cgetg(1, t_MAT);
     999     3210912 :   h = lgcols(X);
    1000     3210894 :   if (!signe(c)) return zeromat(h-1, l-1);
    1001     3165546 :   if (is_pm1(c)) return (signe(c) > 0)? ZM_copy(X): ZM_neg(X);
    1002     2786519 :   A = cgetg(l, t_MAT);
    1003    15321770 :   for (j = 1; j < l; j++)
    1004             :   {
    1005    12535583 :     GEN a = cgetg(h, t_COL), x = gel(X, j);
    1006   221283335 :     for (i = 1; i < h; i++) gel(a,i) = mulii(c, gel(x,i));
    1007    12535264 :     gel(A,j) = a;
    1008             :   }
    1009     2786187 :   return A;
    1010             : }
    1011             : void
    1012    73076029 : ZC_lincomb1_inplace_i(GEN X, GEN Y, GEN v, long n)
    1013             : {
    1014             :   long i;
    1015  1617372702 :   for (i = n; i; i--) gel(X,i) = addmulii_inplace(gel(X,i), gel(Y,i), v);
    1016    73016074 : }
    1017             : /* X <- X + v Y (elementary col operation) */
    1018             : void
    1019    62173260 : ZC_lincomb1_inplace(GEN X, GEN Y, GEN v)
    1020             : {
    1021    62173260 :   if (lgefint(v) != 2) return ZC_lincomb1_inplace_i(X, Y, v, lg(X)-1);
    1022             : }
    1023             : void
    1024    29488948 : Flc_lincomb1_inplace(GEN X, GEN Y, ulong v, ulong q)
    1025             : {
    1026             :   long i;
    1027    29488948 :   if (!v) return; /* v = 0 */
    1028   742146509 :   for (i = lg(X)-1; i; i--) X[i] = Fl_add(X[i], Fl_mul(Y[i], v, q), q);
    1029             : }
    1030             : 
    1031             : /* X + v Y, wasteful if (v = 0) */
    1032             : static GEN
    1033    15393259 : ZC_lincomb1(GEN v, GEN x, GEN y)
    1034   122030790 : { pari_APPLY_type(t_COL, addmulii(gel(x,i), gel(y,i), v)) }
    1035             : 
    1036             : /* -X + vY */
    1037             : static GEN
    1038      741409 : ZC_lincomb_1(GEN v, GEN x, GEN y)
    1039     4586351 : { pari_APPLY_type(t_COL, mulsubii(gel(y,i), v, gel(x,i))) }
    1040             : 
    1041             : /* X,Y compatible ZV; u,v in Z. Returns A = u*X + v*Y */
    1042             : GEN
    1043    33294974 : ZC_lincomb(GEN u, GEN v, GEN X, GEN Y)
    1044             : {
    1045             :   long su, sv;
    1046             :   GEN A;
    1047             : 
    1048    33294974 :   su = signe(u); if (!su) return ZC_Z_mul(Y, v);
    1049    33294022 :   sv = signe(v); if (!sv) return ZC_Z_mul(X, u);
    1050    32864258 :   if (is_pm1(v))
    1051             :   {
    1052    11085243 :     if (is_pm1(u))
    1053             :     {
    1054     9757271 :       if (su != sv) A = ZC_sub(X, Y);
    1055     2709087 :       else          A = ZC_add(X, Y);
    1056     9756843 :       if (su < 0) ZV_togglesign(A); /* in place but was created above */
    1057             :     }
    1058             :     else
    1059             :     {
    1060     1327894 :       if (sv > 0) A = ZC_lincomb1 (u, Y, X);
    1061      607419 :       else        A = ZC_lincomb_1(u, Y, X);
    1062             :     }
    1063             :   }
    1064    21780102 :   else if (is_pm1(u))
    1065             :   {
    1066    14806009 :     if (su > 0) A = ZC_lincomb1 (v, X, Y);
    1067      133211 :     else        A = ZC_lincomb_1(v, X, Y);
    1068             :   }
    1069             :   else
    1070             :   { /* not cgetg_copy: x may be a t_VEC */
    1071     6979018 :     long i, lx = lg(X);
    1072     6979018 :     A = cgetg(lx,t_COL);
    1073    44799072 :     for (i=1; i<lx; i++) gel(A,i) = lincombii(u,v,gel(X,i),gel(Y,i));
    1074             :   }
    1075    32857010 :   return A;
    1076             : }
    1077             : 
    1078             : /********************************************************************/
    1079             : /**                                                                **/
    1080             : /**                           CONVERSIONS                          **/
    1081             : /**                                                                **/
    1082             : /********************************************************************/
    1083             : GEN
    1084      466196 : ZV_to_nv(GEN x)
    1085      883211 : { pari_APPLY_ulong(itou(gel(x,i))) }
    1086             : 
    1087             : GEN
    1088      156492 : zm_to_ZM(GEN x)
    1089      912278 : { pari_APPLY_type(t_MAT, zc_to_ZC(gel(x,i))) }
    1090             : 
    1091             : GEN
    1092         126 : zmV_to_ZMV(GEN x)
    1093         791 : { pari_APPLY_type(t_VEC, zm_to_ZM(gel(x,i))) }
    1094             : 
    1095             : /* same as Flm_to_ZM but do not assume positivity */
    1096             : GEN
    1097        1022 : ZM_to_zm(GEN x)
    1098       17472 : { pari_APPLY_same(ZV_to_zv(gel(x,i))) }
    1099             : 
    1100             : GEN
    1101      366646 : zv_to_Flv(GEN x, ulong p)
    1102     5418812 : { pari_APPLY_ulong(umodsu(x[i], p)) }
    1103             : 
    1104             : GEN
    1105       22694 : zm_to_Flm(GEN x, ulong p)
    1106      351750 : { pari_APPLY_same(zv_to_Flv(gel(x,i),p)) }
    1107             : 
    1108             : GEN
    1109          49 : ZMV_to_zmV(GEN x)
    1110         399 : { pari_APPLY_type(t_VEC, ZM_to_zm(gel(x,i))) }
    1111             : 
    1112             : /********************************************************************/
    1113             : /**                                                                **/
    1114             : /**                         COPY, NEGATION                         **/
    1115             : /**                                                                **/
    1116             : /********************************************************************/
    1117             : GEN
    1118    14714492 : ZC_copy(GEN x)
    1119             : {
    1120    14714492 :   long i, lx = lg(x);
    1121    14714492 :   GEN y = cgetg(lx, t_COL);
    1122   122631566 :   for (i=1; i<lx; i++)
    1123             :   {
    1124   107916697 :     GEN c = gel(x,i);
    1125   107916697 :     gel(y,i) = lgefint(c) == 2? gen_0: icopy(c);
    1126             :   }
    1127    14714869 :   return y;
    1128             : }
    1129             : 
    1130             : GEN
    1131      633471 : ZM_copy(GEN x)
    1132     5809726 : { pari_APPLY_same(ZC_copy(gel(x,i))) }
    1133             : 
    1134             : void
    1135      345096 : ZV_neg_inplace(GEN M)
    1136             : {
    1137      345096 :   long l = lg(M);
    1138     1300157 :   while (--l > 0) gel(M,l) = negi(gel(M,l));
    1139      345147 : }
    1140             : GEN
    1141     6305851 : ZC_neg(GEN x)
    1142    36975264 : { pari_APPLY_type(t_COL, negi(gel(x,i))) }
    1143             : 
    1144             : GEN
    1145       51607 : zv_neg(GEN x)
    1146      662084 : { pari_APPLY_long(-x[i]) }
    1147             : GEN
    1148         126 : zv_neg_inplace(GEN M)
    1149             : {
    1150         126 :   long l = lg(M);
    1151         421 :   while (--l > 0) M[l] = -M[l];
    1152         126 :   return M;
    1153             : }
    1154             : GEN
    1155          77 : zv_abs(GEN x)
    1156        5446 : { pari_APPLY_ulong(labs(x[i])) }
    1157             : GEN
    1158     1657875 : ZM_neg(GEN x)
    1159     5055833 : { pari_APPLY_same(ZC_neg(gel(x,i))) }
    1160             : 
    1161             : void
    1162     5038644 : ZV_togglesign(GEN M)
    1163             : {
    1164     5038644 :   long l = lg(M);
    1165    75519620 :   while (--l > 0) togglesign_safe(&gel(M,l));
    1166     5038761 : }
    1167             : void
    1168           0 : ZM_togglesign(GEN M)
    1169             : {
    1170           0 :   long l = lg(M);
    1171           0 :   while (--l > 0) ZV_togglesign(gel(M,l));
    1172           0 : }
    1173             : 
    1174             : /********************************************************************/
    1175             : /**                                                                **/
    1176             : /**                        "DIVISION" mod HNF                      **/
    1177             : /**                                                                **/
    1178             : /********************************************************************/
    1179             : /* Reduce ZC x modulo ZM y in HNF, may return x itself (not a copy) */
    1180             : GEN
    1181    10737505 : ZC_hnfdivrem(GEN x, GEN y, GEN *Q)
    1182             : {
    1183    10737505 :   long i, l = lg(x);
    1184             :   GEN q;
    1185             : 
    1186    10737505 :   if (Q) *Q = cgetg(l,t_COL);
    1187    10737505 :   if (l == 1) return cgetg(1,t_COL);
    1188    61727898 :   for (i = l-1; i>0; i--)
    1189             :   {
    1190    50993616 :     q = diviiround(gel(x,i), gcoeff(y,i,i));
    1191    50992364 :     if (signe(q)) {
    1192    23013964 :       togglesign(q);
    1193    23014004 :       x = ZC_lincomb(gen_1, q, x, gel(y,i));
    1194             :     }
    1195    50990393 :     if (Q) gel(*Q, i) = q;
    1196             :   }
    1197    10734282 :   return x;
    1198             : }
    1199             : 
    1200             : /* x = y Q + R, may return some columns of x (not copies) */
    1201             : GEN
    1202      454381 : ZM_hnfdivrem(GEN x, GEN y, GEN *Q)
    1203             : {
    1204      454381 :   long lx = lg(x), i;
    1205      454381 :   GEN R = cgetg(lx, t_MAT);
    1206      454399 :   if (Q)
    1207             :   {
    1208      127577 :     GEN q = cgetg(lx, t_MAT); *Q = q;
    1209      188518 :     for (i=1; i<lx; i++) gel(R,i) = ZC_hnfdivrem(gel(x,i),y,(GEN*)(q+i));
    1210             :   }
    1211             :   else
    1212      824074 :     for (i=1; i<lx; i++)
    1213             :     {
    1214      497245 :       pari_sp av = avma;
    1215      497245 :       GEN z = ZC_hnfrem(gel(x,i),y);
    1216      497195 :       gel(R,i) = (avma == av)? ZC_copy(z): gerepileupto(av, z);
    1217             :     }
    1218      454405 :   return R;
    1219             : }
    1220             : 
    1221             : /********************************************************************/
    1222             : /**                                                                **/
    1223             : /**                               TESTS                            **/
    1224             : /**                                                                **/
    1225             : /********************************************************************/
    1226             : int
    1227    23103438 : zv_equal0(GEN V)
    1228             : {
    1229    23103438 :   long l = lg(V);
    1230    37579784 :   while (--l > 0)
    1231    30805579 :     if (V[l]) return 0;
    1232     6774205 :   return 1;
    1233             : }
    1234             : 
    1235             : int
    1236    14456130 : ZV_equal0(GEN V)
    1237             : {
    1238    14456130 :   long l = lg(V);
    1239    25520103 :   while (--l > 0)
    1240    25093902 :     if (signe(gel(V,l))) return 0;
    1241      426201 :   return 1;
    1242             : }
    1243             : int
    1244       16231 : ZMrow_equal0(GEN V, long i)
    1245             : {
    1246       16231 :   long l = lg(V);
    1247       25183 :   while (--l > 0)
    1248       21679 :     if (signe(gcoeff(V,i,l))) return 0;
    1249        3504 :   return 1;
    1250             : }
    1251             : 
    1252             : static int
    1253     6265314 : ZV_equal_lg(GEN V, GEN W, long l)
    1254             : {
    1255    25959025 :   while (--l > 0)
    1256    20175147 :     if (!equalii(gel(V,l), gel(W,l))) return 0;
    1257     5783878 :   return 1;
    1258             : }
    1259             : int
    1260      295030 : ZV_equal(GEN V, GEN W)
    1261             : {
    1262      295030 :   long l = lg(V);
    1263      295030 :   if (lg(W) != l) return 0;
    1264      295023 :   return ZV_equal_lg(V, W, l);
    1265             : }
    1266             : int
    1267     3354796 : ZM_equal(GEN A, GEN B)
    1268             : {
    1269     3354796 :   long i, m, l = lg(A);
    1270     3354796 :   if (lg(B) != l) return 0;
    1271     3354740 :   if (l == 1) return 1;
    1272     3354740 :   m = lgcols(A);
    1273     3354742 :   if (lgcols(B) != m) return 0;
    1274     9062990 :   for (i = 1; i < l; i++)
    1275     5970285 :     if (!ZV_equal_lg(gel(A,i), gel(B,i), m)) return 0;
    1276     3092705 :   return 1;
    1277             : }
    1278             : int
    1279       72442 : ZM_equal0(GEN A)
    1280             : {
    1281       72442 :   long i, j, m, l = lg(A);
    1282       72442 :   if (l == 1) return 1;
    1283       72442 :   m = lgcols(A);
    1284      123568 :   for (j = 1; j < l; j++)
    1285     2718915 :     for (i = 1; i < m; i++)
    1286     2667789 :       if (signe(gcoeff(A,i,j))) return 0;
    1287       15627 :   return 1;
    1288             : }
    1289             : int
    1290    30877388 : zv_equal(GEN V, GEN W)
    1291             : {
    1292    30877388 :   long l = lg(V);
    1293    30877388 :   if (lg(W) != l) return 0;
    1294   262512337 :   while (--l > 0)
    1295   232752057 :     if (V[l] != W[l]) return 0;
    1296    29760280 :   return 1;
    1297             : }
    1298             : 
    1299             : int
    1300        1638 : zvV_equal(GEN V, GEN W)
    1301             : {
    1302        1638 :   long l = lg(V);
    1303        1638 :   if (lg(W) != l) return 0;
    1304       80388 :   while (--l > 0)
    1305       79912 :     if (!zv_equal(gel(V,l),gel(W,l))) return 0;
    1306         476 :   return 1;
    1307             : }
    1308             : 
    1309             : int
    1310      760716 : ZM_ishnf(GEN x)
    1311             : {
    1312      760716 :   long i,j, lx = lg(x);
    1313     2286789 :   for (i=1; i<lx; i++)
    1314             :   {
    1315     1638890 :     GEN xii = gcoeff(x,i,i);
    1316     1638890 :     if (signe(xii) <= 0) return 0;
    1317     3181966 :     for (j=1; j<i; j++)
    1318     1580546 :       if (signe(gcoeff(x,i,j))) return 0;
    1319     3240498 :     for (j=i+1; j<lx; j++)
    1320             :     {
    1321     1714425 :       GEN xij = gcoeff(x,i,j);
    1322     1714425 :       if (signe(xij)<0 || cmpii(xij,xii)>=0) return 0;
    1323             :     }
    1324             :   }
    1325      647899 :   return 1;
    1326             : }
    1327             : int
    1328      658914 : ZM_isidentity(GEN x)
    1329             : {
    1330      658914 :   long i,j, lx = lg(x);
    1331             : 
    1332      658914 :   if (lx == 1) return 1;
    1333      658907 :   if (lx != lgcols(x)) return 0;
    1334     3220593 :   for (j=1; j<lx; j++)
    1335             :   {
    1336     2563744 :     GEN c = gel(x,j);
    1337     8111593 :     for (i=1; i<j; )
    1338     5547885 :       if (signe(gel(c,i++))) return 0;
    1339             :     /* i = j */
    1340     2563708 :     if (!equali1(gel(c,i++))) return 0;
    1341     8117663 :     for (   ; i<lx; )
    1342     5555963 :       if (signe(gel(c,i++))) return 0;
    1343             :   }
    1344      656849 :   return 1;
    1345             : }
    1346             : int
    1347      556561 : ZM_isdiagonal(GEN x)
    1348             : {
    1349      556561 :   long i,j, lx = lg(x);
    1350      556561 :   if (lx == 1) return 1;
    1351      556561 :   if (lx != lgcols(x)) return 0;
    1352             : 
    1353     1438768 :   for (j=1; j<lx; j++)
    1354             :   {
    1355     1209054 :     GEN c = gel(x,j);
    1356     1697558 :     for (i=1; i<j; i++)
    1357      815309 :       if (signe(gel(c,i))) return 0;
    1358     2021450 :     for (i++; i<lx; i++)
    1359     1139236 :       if (signe(gel(c,i))) return 0;
    1360             :   }
    1361      229714 :   return 1;
    1362             : }
    1363             : int
    1364      156428 : ZM_isscalar(GEN x, GEN s)
    1365             : {
    1366      156428 :   long i, j, lx = lg(x);
    1367             : 
    1368      156428 :   if (lx == 1) return 1;
    1369      156428 :   if (!s) s = gcoeff(x,1,1);
    1370      156428 :   if (equali1(s)) return ZM_isidentity(x);
    1371      154972 :   if (lx != lgcols(x)) return 0;
    1372      679687 :   for (j=1; j<lx; j++)
    1373             :   {
    1374      585697 :     GEN c = gel(x,j);
    1375     2641647 :     for (i=1; i<j; )
    1376     2114840 :       if (signe(gel(c,i++))) return 0;
    1377             :     /* i = j */
    1378      526807 :     if (!equalii(gel(c,i++), s)) return 0;
    1379     2668927 :     for (   ; i<lx; )
    1380     2144212 :       if (signe(gel(c,i++))) return 0;
    1381             :   }
    1382       93990 :   return 1;
    1383             : }
    1384             : 
    1385             : long
    1386      159983 : ZC_is_ei(GEN x)
    1387             : {
    1388      159983 :   long i, j = 0, l = lg(x);
    1389     1682600 :   for (i = 1; i < l; i++)
    1390             :   {
    1391     1522618 :     GEN c = gel(x,i);
    1392     1522618 :     long s = signe(c);
    1393     1522618 :     if (!s) continue;
    1394      159977 :     if (s < 0 || !is_pm1(c) || j) return 0;
    1395      159976 :     j = i;
    1396             :   }
    1397      159982 :   return j;
    1398             : }
    1399             : 
    1400             : /********************************************************************/
    1401             : /**                                                                **/
    1402             : /**                       MISCELLANEOUS                            **/
    1403             : /**                                                                **/
    1404             : /********************************************************************/
    1405             : /* assume lg(x) = lg(y), x,y in Z^n */
    1406             : int
    1407     3146987 : ZV_cmp(GEN x, GEN y)
    1408             : {
    1409     3146987 :   long fl,i, lx = lg(x);
    1410     6374612 :   for (i=1; i<lx; i++)
    1411     5079020 :     if (( fl = cmpii(gel(x,i), gel(y,i)) )) return fl;
    1412     1295592 :   return 0;
    1413             : }
    1414             : /* assume lg(x) = lg(y), x,y in Z^n */
    1415             : int
    1416       19740 : ZV_abscmp(GEN x, GEN y)
    1417             : {
    1418       19740 :   long fl,i, lx = lg(x);
    1419       54207 :   for (i=1; i<lx; i++)
    1420       54078 :     if (( fl = abscmpii(gel(x,i), gel(y,i)) )) return fl;
    1421         129 :   return 0;
    1422             : }
    1423             : 
    1424             : long
    1425    18456601 : zv_content(GEN x)
    1426             : {
    1427    18456601 :   long i, s, l = lg(x);
    1428    18456601 :   if (l == 1) return 0;
    1429    18456594 :   s = labs(x[1]);
    1430    42010502 :   for (i = 2; i < l && s != 1; i++) s = ugcd(s, labs(x[i]));
    1431    18456632 :   return s;
    1432             : }
    1433             : GEN
    1434      297309 : ZV_content(GEN x)
    1435             : {
    1436      297309 :   long i, l = lg(x);
    1437      297309 :   pari_sp av = avma;
    1438             :   GEN c;
    1439      297309 :   if (l == 1) return gen_0;
    1440      297309 :   if (l == 2) return absi(gel(x,1));
    1441      204783 :   c = gel(x,1);
    1442      557841 :   for (i = 2; i < l; i++)
    1443             :   {
    1444      403579 :     c = gcdii(c, gel(x,i));
    1445      403579 :     if (is_pm1(c)) { set_avma(av); return gen_1; }
    1446             :   }
    1447      154262 :   return gerepileuptoint(av, c);
    1448             : }
    1449             : 
    1450             : GEN
    1451     3952203 : ZM_det_triangular(GEN mat)
    1452             : {
    1453             :   pari_sp av;
    1454     3952203 :   long i,l = lg(mat);
    1455             :   GEN s;
    1456             : 
    1457     3952203 :   if (l<3) return l<2? gen_1: icopy(gcoeff(mat,1,1));
    1458     3522118 :   av = avma; s = gcoeff(mat,1,1);
    1459     9539607 :   for (i=2; i<l; i++) s = mulii(s,gcoeff(mat,i,i));
    1460     3521645 :   return gerepileuptoint(av,s);
    1461             : }
    1462             : 
    1463             : /* assumes no overflow */
    1464             : long
    1465      950591 : zv_prod(GEN v)
    1466             : {
    1467      950591 :   long n, i, l = lg(v);
    1468      950591 :   if (l == 1) return 1;
    1469      960626 :   n = v[1]; for (i = 2; i < l; i++) n *= v[i];
    1470      772140 :   return n;
    1471             : }
    1472             : 
    1473             : static GEN
    1474   318571124 : _mulii(void *E, GEN a, GEN b)
    1475   318571124 : { (void) E; return mulii(a, b); }
    1476             : 
    1477             : /* product of ulongs */
    1478             : GEN
    1479     1864414 : zv_prod_Z(GEN v)
    1480             : {
    1481             :   pari_sp av;
    1482     1864414 :   long k, m, n = lg(v)-1;
    1483     1864414 :   int stop = 0;
    1484             :   GEN V;
    1485     1864414 :   switch(n) {
    1486       21035 :     case 0: return gen_1;
    1487      130095 :     case 1: return utoi(v[1]);
    1488      977008 :     case 2: return muluu(v[1], v[2]);
    1489             :   }
    1490      736276 :   av = avma; m = n >> 1;
    1491      736276 :   V = cgetg(m + (odd(n)? 2: 1), t_VEC);
    1492   153068415 :   for (k = n; k; k--) /* start from the end: v is usually sorted */
    1493   152333193 :     if (v[k] & HIGHMASK) { stop = 1; break; }
    1494     2420156 :   while (!stop)
    1495             :   { /* HACK: handle V as a t_VECSMALL; gain a few iterations */
    1496    82327157 :     for (k = 1; k <= m; k++)
    1497             :     {
    1498    80034929 :       V[k] = uel(v,k<<1) * uel(v,(k<<1)-1);
    1499    80034929 :       if (V[k] & HIGHMASK) stop = 1; /* last "free" iteration */
    1500             :     }
    1501     2292228 :     if (odd(n))
    1502             :     {
    1503     1350919 :       if (n == 1) { set_avma(av); return utoi(v[1]); }
    1504      742571 :       V[++m] = v[n];
    1505             :     }
    1506     1683880 :     v = V; n = m; m = n >> 1;
    1507             :   }
    1508             :   /* n > 1; m > 0 */
    1509      127928 :   if (n == 2) { set_avma(av); return muluu(v[1], v[2]); }
    1510    46793132 :   for (k = 1; k <= m; k++) gel(V,k) = muluu(v[k<<1], v[(k<<1)-1]);
    1511       87790 :   if (odd(n)) gel(V, ++m) = utoipos(v[n]);
    1512       87794 :   setlg(V, m+1); /* HACK: now V is a bona fide t_VEC */
    1513       87793 :   return gerepileuptoint(av, gen_product(V, NULL, &_mulii));
    1514             : }
    1515             : GEN
    1516    14694393 : vecsmall_prod(GEN v)
    1517             : {
    1518    14694393 :   pari_sp av = avma;
    1519    14694393 :   long k, m, n = lg(v)-1;
    1520             :   GEN V;
    1521    14694393 :   switch (n) {
    1522           0 :     case 0: return gen_1;
    1523           0 :     case 1: return stoi(v[1]);
    1524          21 :     case 2: return mulss(v[1], v[2]);
    1525             :   }
    1526    14694372 :   m = n >> 1;
    1527    14694372 :   V = cgetg(m + (odd(n)? 2: 1), t_VEC);
    1528   161556906 :   for (k = 1; k <= m; k++) gel(V,k) = mulss(v[k<<1], v[(k<<1)-1]);
    1529    14694372 :   if (odd(n)) gel(V,k) = stoi(v[n]);
    1530    14694372 :   return gerepileuptoint(av, gen_product(V, NULL, &_mulii));
    1531             : }
    1532             : 
    1533             : GEN
    1534     8829206 : ZV_prod(GEN v)
    1535             : {
    1536     8829206 :   pari_sp av = avma;
    1537     8829206 :   long i, l = lg(v);
    1538             :   GEN n;
    1539     8829206 :   if (l == 1) return gen_1;
    1540     8644299 :   if (l > 7) return gerepileuptoint(av, gen_product(v, NULL, _mulii));
    1541     1301312 :   n = gel(v,1);
    1542     1301312 :   if (l == 2) return icopy(n);
    1543     2103940 :   for (i = 2; i < l; i++) n = mulii(n, gel(v,i));
    1544      846289 :   return gerepileuptoint(av, n);
    1545             : }
    1546             : /* assumes no overflow */
    1547             : long
    1548       16519 : zv_sum(GEN v)
    1549             : {
    1550       16519 :   long n, i, l = lg(v);
    1551       16519 :   if (l == 1) return 0;
    1552       95898 :   n = v[1]; for (i = 2; i < l; i++) n += v[i];
    1553       16498 :   return n;
    1554             : }
    1555             : /* assumes no overflow and 0 <= n <= #v */
    1556             : long
    1557           0 : zv_sumpart(GEN v, long n)
    1558             : {
    1559             :   long i, p;
    1560           0 :   if (!n) return 0;
    1561           0 :   p = v[1]; for (i = 2; i <= n; i++) p += v[i];
    1562           0 :   return p;
    1563             : }
    1564             : GEN
    1565          77 : ZV_sum(GEN v)
    1566             : {
    1567          77 :   pari_sp av = avma;
    1568          77 :   long i, l = lg(v);
    1569             :   GEN n;
    1570          77 :   if (l == 1) return gen_0;
    1571          77 :   n = gel(v,1);
    1572          77 :   if (l == 2) return icopy(n);
    1573         581 :   for (i = 2; i < l; i++) n = addii(n, gel(v,i));
    1574          77 :   return gerepileuptoint(av, n);
    1575             : }
    1576             : 
    1577             : /********************************************************************/
    1578             : /**                                                                **/
    1579             : /**         GRAM SCHMIDT REDUCTION (integer matrices)              **/
    1580             : /**                                                                **/
    1581             : /********************************************************************/
    1582             : 
    1583             : /* L[k,] += q * L[l,], l < k. Inefficient if q = 0 */
    1584             : static void
    1585      360113 : Zupdate_row(long k, long l, GEN q, GEN L, GEN B)
    1586             : {
    1587      360113 :   long i, qq = itos_or_0(q);
    1588      360113 :   if (!qq)
    1589             :   {
    1590       33369 :     for(i=1;i<l;i++)  gcoeff(L,k,i) = addii(gcoeff(L,k,i),mulii(q,gcoeff(L,l,i)));
    1591        7063 :     gcoeff(L,k,l) = addii(gcoeff(L,k,l), mulii(q,B));
    1592        7063 :     return;
    1593             :   }
    1594      353050 :   if (qq == 1) {
    1595      148699 :     for (i=1;i<l; i++) gcoeff(L,k,i) = addii(gcoeff(L,k,i),gcoeff(L,l,i));
    1596      102229 :     gcoeff(L,k,l) = addii(gcoeff(L,k,l), B);
    1597      250820 :   } else if (qq == -1) {
    1598      149888 :     for (i=1;i<l; i++) gcoeff(L,k,i) = subii(gcoeff(L,k,i),gcoeff(L,l,i));
    1599       87801 :     gcoeff(L,k,l) = addii(gcoeff(L,k,l), negi(B));
    1600             :   } else {
    1601      288498 :     for(i=1;i<l;i++) gcoeff(L,k,i) = addii(gcoeff(L,k,i),mulsi(qq,gcoeff(L,l,i)));
    1602      163040 :     gcoeff(L,k,l) = addii(gcoeff(L,k,l), mulsi(qq,B));
    1603             :   }
    1604             : }
    1605             : 
    1606             : /* update L[k,] */
    1607             : static void
    1608     1035856 : ZRED(long k, long l, GEN x, GEN L, GEN B)
    1609             : {
    1610     1035856 :   GEN q = truedivii(addii(B,shifti(gcoeff(L,k,l),1)), shifti(B,1));
    1611     1035970 :   if (!signe(q)) return;
    1612      360064 :   q = negi(q);
    1613      360093 :   Zupdate_row(k,l,q,L,B);
    1614      360074 :   gel(x,k) = ZC_lincomb(gen_1, q, gel(x,k), gel(x,l));
    1615             : }
    1616             : 
    1617             : /* Gram-Schmidt reduction, x a ZM */
    1618             : static void
    1619     1187833 : ZincrementalGS(GEN x, GEN L, GEN B, long k)
    1620             : {
    1621             :   long i, j;
    1622     3792105 :   for (j=1; j<=k; j++)
    1623             :   {
    1624     2605209 :     pari_sp av = avma;
    1625     2605209 :     GEN u = ZV_dotproduct(gel(x,k), gel(x,j));
    1626     5626082 :     for (i=1; i<j; i++)
    1627             :     {
    1628     3022597 :       u = subii(mulii(gel(B,i+1), u), mulii(gcoeff(L,k,i), gcoeff(L,j,i)));
    1629     3022055 :       u = diviiexact(u, gel(B,i));
    1630             :     }
    1631     2603485 :     gcoeff(L,k,j) = gerepileuptoint(av, u);
    1632             :   }
    1633     1186896 :   gel(B,k+1) = gcoeff(L,k,k); gcoeff(L,k,k) = gen_1;
    1634     1186896 : }
    1635             : 
    1636             : /* Variant reducemodinvertible(ZC v, ZM y), when y singular.
    1637             :  * Very inefficient if y is not LLL-reduced of maximal rank */
    1638             : static GEN
    1639      111732 : ZC_reducemodmatrix_i(GEN v, GEN y)
    1640             : {
    1641      111732 :   GEN B, L, x = shallowconcat(y, v);
    1642      111729 :   long k, lx = lg(x), nx = lx-1;
    1643             : 
    1644      111729 :   B = scalarcol_shallow(gen_1, lx);
    1645      111730 :   L = zeromatcopy(nx, nx);
    1646      459190 :   for (k=1; k <= nx; k++) ZincrementalGS(x, L, B, k);
    1647      347458 :   for (k = nx-1; k >= 1; k--) ZRED(nx,k, x,L,gel(B,k+1));
    1648      111716 :   return gel(x,nx);
    1649             : }
    1650             : GEN
    1651      111732 : ZC_reducemodmatrix(GEN v, GEN y) {
    1652      111732 :   pari_sp av = avma;
    1653      111732 :   return gerepilecopy(av, ZC_reducemodmatrix_i(v,y));
    1654             : }
    1655             : 
    1656             : /* Variant reducemodinvertible(ZM v, ZM y), when y singular.
    1657             :  * Very inefficient if y is not LLL-reduced of maximal rank */
    1658             : static GEN
    1659      227007 : ZM_reducemodmatrix_i(GEN v, GEN y)
    1660             : {
    1661             :   GEN B, L, V;
    1662      227007 :   long j, k, lv = lg(v), nx = lg(y), lx = nx+1;
    1663             : 
    1664      227007 :   V = cgetg(lv, t_MAT);
    1665      227022 :   B = scalarcol_shallow(gen_1, lx);
    1666      227027 :   L = zeromatcopy(nx, nx);
    1667      601474 :   for (k=1; k < nx; k++) ZincrementalGS(y, L, B, k);
    1668      692915 :   for (j = 1; j < lg(v); j++)
    1669             :   {
    1670      465957 :     GEN x = shallowconcat(y, gel(v,j));
    1671      465987 :     ZincrementalGS(x, L, B, nx); /* overwrite last */
    1672     1266121 :     for (k = nx-1; k >= 1; k--) ZRED(nx,k, x,L,gel(B,k+1));
    1673      465897 :     gel(V,j) = gel(x,nx);
    1674             :   }
    1675      226958 :   return V;
    1676             : }
    1677             : GEN
    1678      227006 : ZM_reducemodmatrix(GEN v, GEN y) {
    1679      227006 :   pari_sp av = avma;
    1680      227006 :   return gerepilecopy(av, ZM_reducemodmatrix_i(v,y));
    1681             : }
    1682             : 
    1683             : GEN
    1684       98584 : ZC_reducemodlll(GEN x,GEN y)
    1685             : {
    1686       98584 :   pari_sp av = avma;
    1687       98584 :   GEN z = ZC_reducemodmatrix(x, ZM_lll(y, 0.75, LLL_INPLACE));
    1688       98589 :   return gerepilecopy(av, z);
    1689             : }
    1690             : GEN
    1691           0 : ZM_reducemodlll(GEN x,GEN y)
    1692             : {
    1693           0 :   pari_sp av = avma;
    1694           0 :   GEN z = ZM_reducemodmatrix(x, ZM_lll(y, 0.75, LLL_INPLACE));
    1695           0 :   return gerepilecopy(av, z);
    1696             : }

Generated by: LCOV version 1.16