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.16.2 lcov report (development 29464-65cd88daf0) Lines: 800 911 87.8 %
Date: 2024-07-23 09:03:50 Functions: 122 136 89.7 %
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     1807292 : check_ZV(GEN x, long l)
      20             : {
      21             :   long i;
      22    10566528 :   for (i=1; i<l; i++)
      23     8759327 :     if (typ(gel(x,i)) != t_INT) return 0;
      24     1807201 :   return 1;
      25             : }
      26             : void
      27     1422199 : RgV_check_ZV(GEN A, const char *s)
      28             : {
      29     1422199 :   if (!RgV_is_ZV(A)) pari_err_TYPE(stack_strcat(s," [integer vector]"), A);
      30     1422191 : }
      31             : void
      32      640275 : RgM_check_ZM(GEN A, const char *s)
      33             : {
      34      640275 :   long n = lg(A);
      35      640275 :   if (n != 1)
      36             :   {
      37      640128 :     long j, m = lgcols(A);
      38     2447331 :     for (j=1; j<n; j++)
      39     1807292 :       if (!check_ZV(gel(A,j), m))
      40          91 :         pari_err_TYPE(stack_strcat(s," [integer matrix]"), A);
      41             :   }
      42      640186 : }
      43             : 
      44             : /* assume m > 1 */
      45             : static long
      46   109534117 : ZV_max_lg_i(GEN x, long m)
      47             : {
      48   109534117 :   long i, l = lgefint(gel(x,1));
      49   768950952 :   for (i = 2; i < m; i++) l = maxss(l, lgefint(gel(x,i)));
      50   109534375 :   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    28229832 : ZM_max_lg_i(GEN x, long n, long m)
      62             : {
      63    28229832 :   long j, l = ZV_max_lg_i(gel(x,1), m);
      64   109524215 :   for (j = 2; j < n; j++) l = maxss(l, ZV_max_lg_i(gel(x,j), m));
      65    28230145 :   return l;
      66             : }
      67             : long
      68       21647 : ZM_max_lg(GEN x)
      69             : {
      70       21647 :   long n = lg(x), m;
      71       21647 :   if (n == 1) return 2;
      72       21647 :   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        4001 : ZM_supnorm(GEN x)
     108             : {
     109        4001 :   long i, j, h, lx = lg(x);
     110        4001 :   GEN s = gen_0;
     111        4001 :   if (lx == 1) return gen_1;
     112        4001 :   h = lgcols(x);
     113       24408 :   for (j=1; j<lx; j++)
     114             :   {
     115       20407 :     GEN xj = gel(x,j);
     116      281690 :     for (i=1; i<h; i++)
     117             :     {
     118      261283 :       GEN c = gel(xj,i);
     119      261283 :       if (abscmpii(c, s) > 0) s = c;
     120             :     }
     121             :   }
     122        4001 :   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    20983965 : ZM_zc_mul_i(GEN x, GEN y, long c, long l)
     163             : {
     164             :   long i, j;
     165    20983965 :   GEN z = cgetg(l,t_COL);
     166             : 
     167   126938628 :   for (i=1; i<l; i++)
     168             :   {
     169   105960602 :     pari_sp av = avma;
     170   105960602 :     GEN s = mulis(gcoeff(x,i,1),y[1]);
     171  1178874598 :     for (j=2; j<c; j++)
     172  1072917603 :       if (y[j]) s = addii(s, mulis(gcoeff(x,i,j),y[j]));
     173   105956995 :     gel(z,i) = gerepileuptoint(av,s);
     174             :   }
     175    20978026 :   return z;
     176             : }
     177             : GEN
     178    18989793 : ZM_zc_mul(GEN x, GEN y) {
     179    18989793 :   long l = lg(x);
     180    18989793 :   if (l == 1) return cgetg(1, t_COL);
     181    18989793 :   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      949191 : ZM_zm_mul(GEN x, GEN y)
     205             : {
     206      949191 :   long j, c, l = lg(x), ly = lg(y);
     207      949191 :   GEN z = cgetg(ly, t_MAT);
     208      949191 :   if (l == 1) return z;
     209      949184 :   c = lgcols(x);
     210     2943270 :   for (j = 1; j < ly; j++) gel(z,j) = ZM_zc_mul_i(x, gel(y,j), l,c);
     211      949183 :   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      523928 : 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      523928 :   long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
     235      523928 :   GEN M = cgetg(n + 1, t_MAT), C;
     236             : 
     237     3310990 :   for (j = 1; j <= min_e; j++) {
     238     2787062 :     gel(M, j) = C = cgetg(m + 1, t_COL);
     239    44850181 :     for (i = 1; i <= min_d; i++)
     240    42063119 :       gel(C, i) = addii(gcoeff(A, ma + i, na + j),
     241    42063119 :                         gcoeff(B, mb + i, nb + j));
     242     2853431 :     for (; i <= da; i++) gel(C, i) = gcoeff(A, ma + i, na + j);
     243     2787062 :     for (; i <= db; i++) gel(C, i) = gcoeff(B, mb + i, nb + j);
     244     2787062 :     for (; i <= m; i++)  gel(C, i) = gen_0;
     245             :   }
     246      583324 :   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      523928 :   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      523928 :   for (; j <= n; j++) gel(M, j) = zerocol(m);
     257      523928 :   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      458437 : 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      458437 :   long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
     268      458437 :   GEN M = cgetg(n + 1, t_MAT), C;
     269             : 
     270     2911290 :   for (j = 1; j <= min_e; j++) {
     271     2452853 :     gel(M, j) = C = cgetg(m + 1, t_COL);
     272    40792162 :     for (i = 1; i <= min_d; i++)
     273    38339309 :       gel(C, i) = subii(gcoeff(A, ma + i, na + j),
     274    38339309 :                         gcoeff(B, mb + i, nb + j));
     275     2509954 :     for (; i <= da; i++) gel(C, i) = gcoeff(A, ma + i, na + j);
     276     2540906 :     for (; i <= db; i++) gel(C, i) = negi(gcoeff(B, mb + i, nb + j));
     277     2452853 :     for (; i <= m; i++) gel(C, i) = gen_0;
     278             :   }
     279      458437 :   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      493426 :   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      493426 :   for (; j <= n; j++) gel(M, j) = zerocol(m);
     290      458437 :   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       65491 : ZM_mul_sw(GEN A, GEN B, long m, long n, long p)
     298             : {
     299       65491 :   pari_sp av = avma;
     300       65491 :   long m1 = (m + 1)/2, m2 = m/2,
     301       65491 :     n1 = (n + 1)/2, n2 = n/2,
     302       65491 :     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       65491 :   T2 = subtract_slices(n1, p2, B, 0, n1, p1, p2, B, n1, n2, p1, p2);
     309       65491 :   S1 = subtract_slices(m2, n1, A, m1, m2, 0, n1, A, 0, m2, 0, n1);
     310       65491 :   M2 = ZM_mul_i(S1, T2, m2 + 1, n1 + 1, p2 + 1);
     311       65491 :   if (gc_needed(av, 1))
     312           0 :     gerepileall(av, 2, &T2, &M2);  /* destroy S1 */
     313       65491 :   T3 = subtract_slices(n1, p1, T2, 0, n1, 0, p2, B, 0, n1, 0, p1);
     314       65491 :   if (gc_needed(av, 1))
     315           0 :     gerepileall(av, 2, &M2, &T3);  /* destroy T2 */
     316       65491 :   S2 = add_slices(m2, n1, A, m1, m2, 0, n1, A, m1, m2, n1, n2);
     317       65491 :   T1 = subtract_slices(n1, p1, B, 0, n1, p1, p2, B, 0, n1, 0, p2);
     318       65491 :   M3 = ZM_mul_i(S2, T1, m2 + 1, n1 + 1, p2 + 1);
     319       65491 :   if (gc_needed(av, 1))
     320           0 :     gerepileall(av, 4, &M2, &T3, &S2, &M3);  /* destroy T1 */
     321       65491 :   S3 = subtract_slices(m1, n1, S2, 0, m2, 0, n1, A, 0, m1, 0, n1);
     322       65491 :   if (gc_needed(av, 1))
     323           0 :     gerepileall(av, 4, &M2, &T3, &M3, &S3);  /* destroy S2 */
     324       65491 :   A11 = matslice(A, 1, m1, 1, n1);
     325       65491 :   B11 = matslice(B, 1, n1, 1, p1);
     326       65491 :   M1 = ZM_mul_i(A11, B11, m1 + 1, n1 + 1, p1 + 1);
     327       65491 :   if (gc_needed(av, 1))
     328           0 :     gerepileall(av, 5, &M2, &T3, &M3, &S3, &M1);  /* destroy A11, B11 */
     329       65491 :   A12 = matslice(A, 1, m1, n1 + 1, n);
     330       65491 :   B21 = matslice(B, n1 + 1, n, 1, p1);
     331       65491 :   M4 = ZM_mul_i(A12, B21, m1 + 1, n2 + 1, p1 + 1);
     332       65491 :   if (gc_needed(av, 1))
     333           0 :     gerepileall(av, 6, &M2, &T3, &M3, &S3, &M1, &M4);  /* destroy A12, B21 */
     334       65491 :   C11 = add_slices(m1, p1, M1, 0, m1, 0, p1, M4, 0, m1, 0, p1);
     335       65491 :   if (gc_needed(av, 1))
     336           0 :     gerepileall(av, 6, &M2, &T3, &M3, &S3, &M1, &C11);  /* destroy M4 */
     337       65491 :   M5 = ZM_mul_i(S3, T3, m1 + 1, n1 + 1, p1 + 1);
     338       65491 :   S4 = subtract_slices(m1, n2, A, 0, m1, n1, n2, S3, 0, m1, 0, n2);
     339       65491 :   if (gc_needed(av, 1))
     340           5 :     gerepileall(av, 7, &M2, &T3, &M3, &M1, &C11, &M5, &S4);  /* destroy S3 */
     341       65491 :   T4 = add_slices(n2, p1, B, n1, n2, 0, p1, T3, 0, n2, 0, p1);
     342       65491 :   if (gc_needed(av, 1))
     343           0 :     gerepileall(av, 7, &M2, &M3, &M1, &C11, &M5, &S4, &T4);  /* destroy T3 */
     344       65491 :   V1 = subtract_slices(m1, p1, M1, 0, m1, 0, p1, M5, 0, m1, 0, p1);
     345       65491 :   if (gc_needed(av, 1))
     346           1 :     gerepileall(av, 6, &M2, &M3, &S4, &T4, &C11, &V1);  /* destroy M1, M5 */
     347       65491 :   B22 = matslice(B, n1 + 1, n, p1 + 1, p);
     348       65491 :   M6 = ZM_mul_i(S4, B22, m1 + 1, n2 + 1, p2 + 1);
     349       65491 :   if (gc_needed(av, 1))
     350           6 :     gerepileall(av, 6, &M2, &M3, &T4, &C11, &V1, &M6);  /* destroy S4, B22 */
     351       65491 :   A22 = matslice(A, m1 + 1, m, n1 + 1, n);
     352       65491 :   M7 = ZM_mul_i(A22, T4, m2 + 1, n2 + 1, p1 + 1);
     353       65491 :   if (gc_needed(av, 1))
     354           0 :     gerepileall(av, 6, &M2, &M3, &C11, &V1, &M6, &M7);  /* destroy A22, T4 */
     355       65491 :   V3 = add_slices(m1, p2, V1, 0, m1, 0, p2, M3, 0, m2, 0, p2);
     356       65491 :   C12 = add_slices(m1, p2, V3, 0, m1, 0, p2, M6, 0, m1, 0, p2);
     357       65491 :   if (gc_needed(av, 1))
     358           6 :     gerepileall(av, 6, &M2, &M3, &C11, &V1, &M7, &C12);  /* destroy V3, M6 */
     359       65491 :   V2 = add_slices(m2, p1, V1, 0, m2, 0, p1, M2, 0, m2, 0, p2);
     360       65491 :   if (gc_needed(av, 1))
     361           0 :     gerepileall(av, 5, &M3, &C11, &M7, &C12, &V2);  /* destroy V1, M2 */
     362       65491 :   C21 = add_slices(m2, p1, V2, 0, m2, 0, p1, M7, 0, m2, 0, p1);
     363       65491 :   if (gc_needed(av, 1))
     364           6 :     gerepileall(av, 5, &M3, &C11, &C12, &V2, &C21);  /* destroy M7 */
     365       65491 :   C22 = add_slices(m2, p2, V2, 0, m2, 0, p2, M3, 0, m2, 0, p2);
     366       65491 :   if (gc_needed(av, 1))
     367           0 :     gerepileall(av, 4, &C11, &C12, &C21, &C22);  /* destroy V2, M3 */
     368       65491 :   C = shallowconcat(vconcat(C11, C21), vconcat(C12, C22));
     369       65491 :   return gerepilecopy(av, C);
     370             : }
     371             : 
     372             : /* x[i,]*y. Assume lg(x) > 1 and 0 < i < lgcols(x) */
     373             : static GEN
     374   529796552 : ZMrow_ZC_mul_i(GEN x, GEN y, long i, long lx)
     375             : {
     376   529796552 :   pari_sp av = avma;
     377   529796552 :   GEN c = mulii(gcoeff(x,i,1), gel(y,1)), ZERO = gen_0;
     378             :   long k;
     379  5296723776 :   for (k = 2; k < lx; k++)
     380             :   {
     381  4767257943 :     GEN t = mulii(gcoeff(x,i,k), gel(y,k));
     382  4766050130 :     if (t != ZERO) c = addii(c, t);
     383             :   }
     384   529465833 :   return gerepileuptoint(av, c);
     385             : }
     386             : GEN
     387   134539740 : ZMrow_ZC_mul(GEN x, GEN y, long i)
     388   134539740 : { 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    74698769 : ZM_ZC_mul_i(GEN x, GEN y, long lx, long l)
     393             : {
     394    74698769 :   GEN z = cgetg(l,t_COL);
     395             :   long i;
     396   469974602 :   for (i=1; i<l; i++) gel(z,i) = ZMrow_ZC_mul_i(x,y,i,lx);
     397    74660350 :   return z;
     398             : }
     399             : 
     400             : static GEN
     401    14049751 : ZM_mul_classical(GEN x, GEN y, long l, long lx, long ly)
     402             : {
     403             :   long j;
     404    14049751 :   GEN z = cgetg(ly, t_MAT);
     405    65763272 :   for (j = 1; j < ly; j++)
     406    51716515 :     gel(z, j) = ZM_ZC_mul_i(x, gel(y, j), lx, l);
     407    14046757 :   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    14115246 : sw_bound(long s)
     462    14115246 : { 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    21285675 : 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    18071825 :   long Flm_sw_bound = 70;
     472             : #else
     473     3213850 :   long Flm_sw_bound = 120;
     474             : #endif
     475    21285675 :   if (l == 1) return zeromat(0, ly-1);
     476    21283659 :   if (lx==2 && l==2 && ly==2)
     477      357335 :   { retmkmat(mkcol(mulii(gcoeff(x,1,1), gcoeff(y,1,1)))); }
     478    20926324 :   if (lx==3 && l==3 && ly==3) return ZM2_mul(x, y);
     479    14091926 :   sx = ZM_max_lg_i(x, lx, l);
     480    14092620 :   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    14092697 :   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    14091843 :   B = sw_bound(minss(sx, sy));
     487    14091826 :   if (l <= B || lx <= B || ly <= B)
     488    14026365 :     return ZM_mul_classical(x, y, l, lx, ly);
     489             :   else
     490       65461 :     return ZM_mul_sw(x, y, l - 1, lx - 1, ly - 1);
     491             : }
     492             : 
     493             : GEN
     494    20963566 : ZM_mul(GEN x, GEN y)
     495             : {
     496    20963566 :   long lx = lg(x), ly = lg(y);
     497    20963566 :   if (ly == 1) return cgetg(1,t_MAT);
     498    20828774 :   if (lx == 1) return zeromat(0, ly-1);
     499    20827192 :   return ZM_mul_i(x, y, lgcols(x), lx, ly);
     500             : }
     501             : 
     502             : GEN
     503      548927 : QM_mul(GEN x, GEN y)
     504             : {
     505      548927 :   GEN dx, nx = Q_primitive_part(x, &dx);
     506      548927 :   GEN dy, ny = Q_primitive_part(y, &dy);
     507      548927 :   GEN z = ZM_mul(nx, ny);
     508      548927 :   if (dx || dy)
     509             :   {
     510      465316 :     GEN d = dx ? dy ? gmul(dx, dy): dx : dy;
     511      465316 :     if (!gequal1(d)) z = ZM_Q_mul(z, d);
     512             :   }
     513      548927 :   return z;
     514             : }
     515             : 
     516             : GEN
     517         700 : QM_sqr(GEN x)
     518             : {
     519         700 :   GEN dx, nx = Q_primitive_part(x, &dx);
     520         700 :   GEN z = ZM_sqr(nx);
     521         700 :   if (dx)
     522         700 :     z = ZM_Q_mul(z, gsqr(dx));
     523         700 :   return z;
     524             : }
     525             : 
     526             : GEN
     527      142082 : QM_QC_mul(GEN x, GEN y)
     528             : {
     529      142082 :   GEN dx, nx = Q_primitive_part(x, &dx);
     530      142082 :   GEN dy, ny = Q_primitive_part(y, &dy);
     531      142082 :   GEN z = ZM_ZC_mul(nx, ny);
     532      142082 :   if (dx || dy)
     533             :   {
     534      142054 :     GEN d = dx ? dy ? gmul(dx, dy): dx : dy;
     535      142054 :     if (!gequal1(d)) z = ZC_Q_mul(z, d);
     536             :   }
     537      142082 :   return z;
     538             : }
     539             : 
     540             : /* assume result is symmetric */
     541             : GEN
     542           0 : ZM_multosym(GEN x, GEN y)
     543             : {
     544           0 :   long j, lx, ly = lg(y);
     545             :   GEN M;
     546           0 :   if (ly == 1) return cgetg(1,t_MAT);
     547           0 :   lx = lg(x); /* = lgcols(y) */
     548           0 :   if (lx == 1) return cgetg(1,t_MAT);
     549             :   /* ly = lgcols(x) */
     550           0 :   M = cgetg(ly, t_MAT);
     551           0 :   for (j=1; j<ly; j++)
     552             :   {
     553           0 :     GEN z = cgetg(ly,t_COL), yj = gel(y,j);
     554             :     long i;
     555           0 :     for (i=1; i<j; i++) gel(z,i) = gcoeff(M,j,i);
     556           0 :     for (i=j; i<ly; i++)gel(z,i) = ZMrow_ZC_mul_i(x,yj,i,lx);
     557           0 :     gel(M,j) = z;
     558             :   }
     559           0 :   return M;
     560             : }
     561             : 
     562             : /* compute m*diagonal(d), assume lg(d) = lg(m). Shallow */
     563             : GEN
     564          21 : ZM_mul_diag(GEN m, GEN d)
     565             : {
     566             :   long j, l;
     567          21 :   GEN y = cgetg_copy(m, &l);
     568          77 :   for (j=1; j<l; j++)
     569             :   {
     570          56 :     GEN c = gel(d,j);
     571          56 :     gel(y,j) = equali1(c)? gel(m,j): ZC_Z_mul(gel(m,j), c);
     572             :   }
     573          21 :   return y;
     574             : }
     575             : /* compute diagonal(d)*m, assume lg(d) = lg(m~). Shallow */
     576             : GEN
     577      533515 : ZM_diag_mul(GEN d, GEN m)
     578             : {
     579      533515 :   long i, j, l = lg(d), lm = lg(m);
     580      533515 :   GEN y = cgetg(lm, t_MAT);
     581     1520170 :   for (j=1; j<lm; j++) gel(y,j) = cgetg(l, t_COL);
     582     1646690 :   for (i=1; i<l; i++)
     583             :   {
     584     1113431 :     GEN c = gel(d,i);
     585     1113431 :     if (equali1(c))
     586      236738 :       for (j=1; j<lm; j++) gcoeff(y,i,j) = gcoeff(m,i,j);
     587             :     else
     588     3396996 :       for (j=1; j<lm; j++) gcoeff(y,i,j) = mulii(gcoeff(m,i,j), c);
     589             :   }
     590      533259 :   return y;
     591             : }
     592             : 
     593             : /* assume lx > 1 is lg(x) = lg(y) */
     594             : static GEN
     595    19199113 : ZV_dotproduct_i(GEN x, GEN y, long lx)
     596             : {
     597    19199113 :   pari_sp av = avma;
     598    19199113 :   GEN c = mulii(gel(x,1), gel(y,1));
     599             :   long i;
     600   146117608 :   for (i = 2; i < lx; i++)
     601             :   {
     602   126919628 :     GEN t = mulii(gel(x,i), gel(y,i));
     603   126919155 :     if (t != gen_0) c = addii(c, t);
     604             :   }
     605    19197980 :   return gerepileuptoint(av, c);
     606             : }
     607             : 
     608             : /* x~ * y, assuming result is symmetric */
     609             : GEN
     610      527937 : ZM_transmultosym(GEN x, GEN y)
     611             : {
     612      527937 :   long i, j, l, ly = lg(y);
     613             :   GEN M;
     614      527937 :   if (ly == 1) return cgetg(1,t_MAT);
     615             :   /* lg(x) = ly */
     616      527937 :   l = lgcols(y); /* = lgcols(x) */
     617      527937 :   M = cgetg(ly, t_MAT);
     618     2696018 :   for (i=1; i<ly; i++)
     619             :   {
     620     2168081 :     GEN xi = gel(x,i), c = cgetg(ly,t_COL);
     621     2168081 :     gel(M,i) = c;
     622     7154255 :     for (j=1; j<i; j++)
     623     4986174 :       gcoeff(M,i,j) = gel(c,j) = ZV_dotproduct_i(xi,gel(y,j),l);
     624     2168081 :     gel(c,i) = ZV_dotproduct_i(xi,gel(y,i),l);
     625             :   }
     626      527937 :   return M;
     627             : }
     628             : /* x~ * y */
     629             : GEN
     630        2289 : ZM_transmul(GEN x, GEN y)
     631             : {
     632        2289 :   long i, j, l, lx, ly = lg(y);
     633             :   GEN M;
     634        2289 :   if (ly == 1) return cgetg(1,t_MAT);
     635        2289 :   lx = lg(x);
     636        2289 :   l = lgcols(y);
     637        2289 :   if (lgcols(x) != l) pari_err_OP("operation 'ZM_transmul'", x,y);
     638        2289 :   M = cgetg(ly, t_MAT);
     639        6993 :   for (i=1; i<ly; i++)
     640             :   {
     641        4704 :     GEN yi = gel(y,i), c = cgetg(lx,t_COL);
     642        4704 :     gel(M,i) = c;
     643       12229 :     for (j=1; j<lx; j++) gel(c,j) = ZV_dotproduct_i(yi,gel(x,j),l);
     644             :   }
     645        2289 :   return M;
     646             : }
     647             : 
     648             : /* assume l > 1; x is (l-1) x (l-1), return x^2.
     649             :  * FIXME: we ultimately rely on Strassen-Winograd which uses 7M + 15A.
     650             :  * Should use Bodrato's variant of Winograd, using 3M + 4S + 11A */
     651             : static GEN
     652       24192 : ZM_sqr_i(GEN x, long l)
     653             : {
     654             :   long s;
     655       24192 :   if (l == 2) { retmkmat(mkcol(sqri(gcoeff(x,1,1)))); }
     656       24192 :   if (l == 3) return ZM2_sqr(x);
     657       23422 :   s = ZM_max_lg_i(x, l, l);
     658       23422 :   if (l > 70) return ZM_mul_fast(x, x, l, l, s, s);
     659       23422 :   if (l <= sw_bound(s))
     660       23392 :     return ZM_mul_classical(x, x, l, l, l);
     661             :   else
     662          30 :     return ZM_mul_sw(x, x, l - 1, l - 1, l - 1);
     663             : }
     664             : 
     665             : GEN
     666       24192 : ZM_sqr(GEN x)
     667             : {
     668       24192 :   long lx=lg(x);
     669       24192 :   if (lx==1) return cgetg(1,t_MAT);
     670       24192 :   return ZM_sqr_i(x, lx);
     671             : }
     672             : GEN
     673    23070928 : ZM_ZC_mul(GEN x, GEN y)
     674             : {
     675    23070928 :   long lx = lg(x);
     676    23070928 :   return lx==1? cgetg(1,t_COL): ZM_ZC_mul_i(x, y, lx, lgcols(x));
     677             : }
     678             : 
     679             : GEN
     680     4493168 : ZC_Z_div(GEN x, GEN c)
     681    19443725 : { pari_APPLY_type(t_COL, Qdivii(gel(x,i), c)) }
     682             : 
     683             : GEN
     684       15782 : ZM_Z_div(GEN x, GEN c)
     685      177995 : { pari_APPLY_same(ZC_Z_div(gel(x, i), c)) }
     686             : 
     687             : GEN
     688     2697408 : ZC_Q_mul(GEN A, GEN z)
     689             : {
     690     2697408 :   pari_sp av = avma;
     691     2697408 :   long i, l = lg(A);
     692             :   GEN d, n, Ad, B, u;
     693     2697408 :   if (typ(z)==t_INT) return ZC_Z_mul(A,z);
     694     2691150 :   n = gel(z, 1); d = gel(z, 2);
     695     2691150 :   Ad = FpC_red(A, d);
     696     2691095 :   u = gcdii(d, FpV_factorback(Ad, NULL, d));
     697     2691085 :   B = cgetg(l, t_COL);
     698     2691081 :   if (equali1(u))
     699             :   {
     700      420633 :     for(i=1; i<l; i++)
     701      354570 :       gel(B, i) = mkfrac(mulii(n, gel(A,i)), d);
     702             :   } else
     703             :   {
     704    18180874 :     for(i=1; i<l; i++)
     705             :     {
     706    15555931 :       GEN di = gcdii(gel(Ad, i), u), ni = mulii(n, diviiexact(gel(A,i), di));
     707    15555884 :       if (equalii(d, di))
     708    10626567 :         gel(B, i) = ni;
     709             :       else
     710     4929281 :         gel(B, i) = mkfrac(ni, diviiexact(d, di));
     711             :     }
     712             :   }
     713     2691006 :   return gerepilecopy(av, B);
     714             : }
     715             : 
     716             : GEN
     717     1092204 : ZM_Q_mul(GEN x, GEN z)
     718             : {
     719     1092204 :   if (typ(z)==t_INT) return ZM_Z_mul(x,z);
     720     3214980 :   pari_APPLY_same(ZC_Q_mul(gel(x, i), z));
     721             : }
     722             : 
     723             : long
     724   196399119 : zv_dotproduct(GEN x, GEN y)
     725             : {
     726   196399119 :   long i, lx = lg(x);
     727             :   ulong c;
     728   196399119 :   if (lx == 1) return 0;
     729   196399119 :   c = uel(x,1)*uel(y,1);
     730  3047359154 :   for (i = 2; i < lx; i++)
     731  2850960035 :     c += uel(x,i)*uel(y,i);
     732   196399119 :   return c;
     733             : }
     734             : 
     735             : GEN
     736      230629 : ZV_ZM_mul(GEN x, GEN y)
     737             : {
     738      230629 :   long i, lx = lg(x), ly = lg(y);
     739             :   GEN z;
     740      230629 :   if (lx == 1) return zerovec(ly-1);
     741      230510 :   z = cgetg(ly, t_VEC);
     742      882063 :   for (i = 1; i < ly; i++) gel(z,i) = ZV_dotproduct_i(x, gel(y,i), lx);
     743      230510 :   return z;
     744             : }
     745             : 
     746             : GEN
     747           0 : ZC_ZV_mul(GEN x, GEN y)
     748             : {
     749           0 :   long i,j, lx=lg(x), ly=lg(y);
     750             :   GEN z;
     751           0 :   if (ly==1) return cgetg(1,t_MAT);
     752           0 :   z = cgetg(ly,t_MAT);
     753           0 :   for (j=1; j < ly; j++)
     754             :   {
     755           0 :     gel(z,j) = cgetg(lx,t_COL);
     756           0 :     for (i=1; i<lx; i++) gcoeff(z,i,j) = mulii(gel(x,i),gel(y,j));
     757             :   }
     758           0 :   return z;
     759             : }
     760             : 
     761             : GEN
     762     6636321 : ZV_dotsquare(GEN x)
     763             : {
     764             :   long i, lx;
     765             :   pari_sp av;
     766             :   GEN z;
     767     6636321 :   lx = lg(x);
     768     6636321 :   if (lx == 1) return gen_0;
     769     6636321 :   av = avma; z = sqri(gel(x,1));
     770    26297024 :   for (i=2; i<lx; i++) z = addii(z, sqri(gel(x,i)));
     771     6634961 :   return gerepileuptoint(av,z);
     772             : }
     773             : 
     774             : GEN
     775    16154047 : ZV_dotproduct(GEN x,GEN y)
     776             : {
     777             :   long lx;
     778    16154047 :   if (x == y) return ZV_dotsquare(x);
     779    11384978 :   lx = lg(x);
     780    11384978 :   if (lx == 1) return gen_0;
     781    11384978 :   return ZV_dotproduct_i(x, y, lx);
     782             : }
     783             : 
     784             : static GEN
     785         280 : _ZM_mul(void *data /*ignored*/, GEN x, GEN y)
     786         280 : { (void)data; return ZM_mul(x,y); }
     787             : static GEN
     788       23184 : _ZM_sqr(void *data /*ignored*/, GEN x)
     789       23184 : { (void)data; return ZM_sqr(x); }
     790             : /* FIXME: Using Bodrato's squaring, precomputations attached to fixed
     791             :  * multiplicand should be reused. And some postcomputations can be fused */
     792             : GEN
     793           0 : ZM_pow(GEN x, GEN n)
     794             : {
     795           0 :   pari_sp av = avma;
     796           0 :   if (!signe(n)) return matid(lg(x)-1);
     797           0 :   return gerepilecopy(av, gen_pow_i(x, n, NULL, &_ZM_sqr, &_ZM_mul));
     798             : }
     799             : GEN
     800       22638 : ZM_powu(GEN x, ulong n)
     801             : {
     802       22638 :   pari_sp av = avma;
     803       22638 :   if (!n) return matid(lg(x)-1);
     804       22638 :   return gerepilecopy(av, gen_powu_i(x, n, NULL, &_ZM_sqr, &_ZM_mul));
     805             : }
     806             : /********************************************************************/
     807             : /**                                                                **/
     808             : /**                           ADD, SUB                             **/
     809             : /**                                                                **/
     810             : /********************************************************************/
     811             : static GEN
     812    34468123 : ZC_add_i(GEN x, GEN y, long lx)
     813             : {
     814    34468123 :   GEN A = cgetg(lx, t_COL);
     815             :   long i;
     816   442108959 :   for (i=1; i<lx; i++) gel(A,i) = addii(gel(x,i), gel(y,i));
     817    34461100 :   return A;
     818             : }
     819             : GEN
     820    26386894 : ZC_add(GEN x, GEN y) { return ZC_add_i(x, y, lg(x)); }
     821             : GEN
     822      363165 : ZC_Z_add(GEN x, GEN y)
     823             : {
     824      363165 :   long k, lx = lg(x);
     825      363165 :   GEN z = cgetg(lx, t_COL);
     826      363164 :   if (lx == 1) pari_err_TYPE2("+",x,y);
     827      363164 :   gel(z,1) = addii(y,gel(x,1));
     828     2500489 :   for (k = 2; k < lx; k++) gel(z,k) = icopy(gel(x,k));
     829      363174 :   return z;
     830             : }
     831             : 
     832             : static GEN
     833     8952072 : ZC_sub_i(GEN x, GEN y, long lx)
     834             : {
     835             :   long i;
     836     8952072 :   GEN A = cgetg(lx, t_COL);
     837    64022684 :   for (i=1; i<lx; i++) gel(A,i) = subii(gel(x,i), gel(y,i));
     838     8951611 :   return A;
     839             : }
     840             : GEN
     841     8722465 : ZC_sub(GEN x, GEN y) { return ZC_sub_i(x, y, lg(x)); }
     842             : GEN
     843           0 : ZC_Z_sub(GEN x, GEN y)
     844             : {
     845           0 :   long k, lx = lg(x);
     846           0 :   GEN z = cgetg(lx, t_COL);
     847           0 :   if (lx == 1) pari_err_TYPE2("+",x,y);
     848           0 :   gel(z,1) = subii(gel(x,1), y);
     849           0 :   for (k = 2; k < lx; k++) gel(z,k) = icopy(gel(x,k));
     850           0 :   return z;
     851             : }
     852             : GEN
     853      544244 : Z_ZC_sub(GEN a, GEN x)
     854             : {
     855      544244 :   long k, lx = lg(x);
     856      544244 :   GEN z = cgetg(lx, t_COL);
     857      544240 :   if (lx == 1) pari_err_TYPE2("-",a,x);
     858      544240 :   gel(z,1) = subii(a, gel(x,1));
     859     1508449 :   for (k = 2; k < lx; k++) gel(z,k) = negi(gel(x,k));
     860      544239 :   return z;
     861             : }
     862             : 
     863             : GEN
     864      747481 : ZM_add(GEN x, GEN y)
     865             : {
     866      747481 :   long lx = lg(x), l, j;
     867             :   GEN z;
     868      747481 :   if (lx == 1) return cgetg(1, t_MAT);
     869      668143 :   z = cgetg(lx, t_MAT); l = lgcols(x);
     870     8749275 :   for (j = 1; j < lx; j++) gel(z,j) = ZC_add_i(gel(x,j), gel(y,j), l);
     871      668145 :   return z;
     872             : }
     873             : GEN
     874       65301 : ZM_sub(GEN x, GEN y)
     875             : {
     876       65301 :   long lx = lg(x), l, j;
     877             :   GEN z;
     878       65301 :   if (lx == 1) return cgetg(1, t_MAT);
     879       65301 :   z = cgetg(lx, t_MAT); l = lgcols(x);
     880      294812 :   for (j = 1; j < lx; j++) gel(z,j) = ZC_sub_i(gel(x,j), gel(y,j), l);
     881       65300 :   return z;
     882             : }
     883             : /********************************************************************/
     884             : /**                                                                **/
     885             : /**                         LINEAR COMBINATION                     **/
     886             : /**                                                                **/
     887             : /********************************************************************/
     888             : /* return X/c assuming division is exact */
     889             : GEN
     890     4432056 : ZC_Z_divexact(GEN x, GEN c)
     891    45853773 : { pari_APPLY_type(t_COL, diviiexact(gel(x,i), c)) }
     892             : GEN
     893        2261 : ZC_divexactu(GEN x, ulong c)
     894       11375 : { pari_APPLY_type(t_COL, diviuexact(gel(x,i), c)) }
     895             : 
     896             : GEN
     897      429812 : ZM_Z_divexact(GEN x, GEN c)
     898     2798525 : { pari_APPLY_same(ZC_Z_divexact(gel(x,i), c)) }
     899             : 
     900             : GEN
     901         441 : ZM_divexactu(GEN x, ulong c)
     902        2702 : { pari_APPLY_same(ZC_divexactu(gel(x,i), c)) }
     903             : 
     904             : GEN
     905    35108428 : ZC_Z_mul(GEN x, GEN c)
     906             : {
     907    35108428 :   if (!signe(c)) return zerocol(lg(x)-1);
     908    33752230 :   if (is_pm1(c)) return (signe(c) > 0)? ZC_copy(x): ZC_neg(x);
     909   254329475 :   pari_APPLY_type(t_COL, mulii(gel(x,i), c))
     910             : }
     911             : 
     912             : GEN
     913       61732 : ZC_z_mul(GEN x, long c)
     914             : {
     915       61732 :   if (!c) return zerocol(lg(x)-1);
     916       53245 :   if (c == 1) return ZC_copy(x);
     917       48391 :   if (c ==-1) return ZC_neg(x);
     918      479041 :   pari_APPLY_type(t_COL, mulsi(c, gel(x,i)))
     919             : }
     920             : 
     921             : GEN
     922       28346 : zv_z_mul(GEN x, long n)
     923      430551 : { pari_APPLY_long(x[i]*n) }
     924             : 
     925             : /* return a ZM */
     926             : GEN
     927           0 : nm_Z_mul(GEN X, GEN c)
     928             : {
     929           0 :   long i, j, h, l = lg(X), s = signe(c);
     930             :   GEN A;
     931           0 :   if (l == 1) return cgetg(1, t_MAT);
     932           0 :   h = lgcols(X);
     933           0 :   if (!s) return zeromat(h-1, l-1);
     934           0 :   if (is_pm1(c)) {
     935           0 :     if (s > 0) return Flm_to_ZM(X);
     936           0 :     X = Flm_to_ZM(X); ZM_togglesign(X); return X;
     937             :   }
     938           0 :   A = cgetg(l, t_MAT);
     939           0 :   for (j = 1; j < l; j++)
     940             :   {
     941           0 :     GEN a = cgetg(h, t_COL), x = gel(X, j);
     942           0 :     for (i = 1; i < h; i++) gel(a,i) = muliu(c, x[i]);
     943           0 :     gel(A,j) = a;
     944             :   }
     945           0 :   return A;
     946             : }
     947             : GEN
     948     2984575 : ZM_Z_mul(GEN X, GEN c)
     949             : {
     950     2984575 :   long i, j, h, l = lg(X);
     951             :   GEN A;
     952     2984575 :   if (l == 1) return cgetg(1, t_MAT);
     953     2984575 :   h = lgcols(X);
     954     2984569 :   if (!signe(c)) return zeromat(h-1, l-1);
     955     2939302 :   if (is_pm1(c)) return (signe(c) > 0)? ZM_copy(X): ZM_neg(X);
     956     2653938 :   A = cgetg(l, t_MAT);
     957    13337971 :   for (j = 1; j < l; j++)
     958             :   {
     959    10684276 :     GEN a = cgetg(h, t_COL), x = gel(X, j);
     960   168376746 :     for (i = 1; i < h; i++) gel(a,i) = mulii(c, gel(x,i));
     961    10684032 :     gel(A,j) = a;
     962             :   }
     963     2653695 :   return A;
     964             : }
     965             : void
     966    71944471 : ZC_lincomb1_inplace_i(GEN X, GEN Y, GEN v, long n)
     967             : {
     968             :   long i;
     969  1564288710 :   for (i = n; i; i--) gel(X,i) = addmulii_inplace(gel(X,i), gel(Y,i), v);
     970    71897677 : }
     971             : /* X <- X + v Y (elementary col operation) */
     972             : void
     973    61210079 : ZC_lincomb1_inplace(GEN X, GEN Y, GEN v)
     974             : {
     975    61210079 :   if (lgefint(v) != 2) return ZC_lincomb1_inplace_i(X, Y, v, lg(X)-1);
     976             : }
     977             : void
     978    29461923 : Flc_lincomb1_inplace(GEN X, GEN Y, ulong v, ulong q)
     979             : {
     980             :   long i;
     981    29461923 :   if (!v) return; /* v = 0 */
     982   742152953 :   for (i = lg(X)-1; i; i--) X[i] = Fl_add(X[i], Fl_mul(Y[i], v, q), q);
     983             : }
     984             : 
     985             : /* X + v Y, wasteful if (v = 0) */
     986             : static GEN
     987    15322779 : ZC_lincomb1(GEN v, GEN x, GEN y)
     988   121653058 : { pari_APPLY_type(t_COL, addmulii(gel(x,i), gel(y,i), v)) }
     989             : 
     990             : /* -X + vY */
     991             : static GEN
     992      729333 : ZC_lincomb_1(GEN v, GEN x, GEN y)
     993     4529180 : { pari_APPLY_type(t_COL, mulsubii(gel(y,i), v, gel(x,i))) }
     994             : 
     995             : /* X,Y compatible ZV; u,v in Z. Returns A = u*X + v*Y */
     996             : GEN
     997    32857096 : ZC_lincomb(GEN u, GEN v, GEN X, GEN Y)
     998             : {
     999             :   long su, sv;
    1000             :   GEN A;
    1001             : 
    1002    32857096 :   su = signe(u); if (!su) return ZC_Z_mul(Y, v);
    1003    32856137 :   sv = signe(v); if (!sv) return ZC_Z_mul(X, u);
    1004    32426928 :   if (is_pm1(v))
    1005             :   {
    1006    11023086 :     if (is_pm1(u))
    1007             :     {
    1008     9717633 :       if (su != sv) A = ZC_sub(X, Y);
    1009     2700557 :       else          A = ZC_add(X, Y);
    1010     9717238 :       if (su < 0) ZV_togglesign(A); /* in place but was created above */
    1011             :     }
    1012             :     else
    1013             :     {
    1014     1305347 :       if (sv > 0) A = ZC_lincomb1 (u, Y, X);
    1015      600164 :       else        A = ZC_lincomb_1(u, Y, X);
    1016             :     }
    1017             :   }
    1018    21404498 :   else if (is_pm1(u))
    1019             :   {
    1020    14746093 :     if (su > 0) A = ZC_lincomb1 (v, X, Y);
    1021      128478 :     else        A = ZC_lincomb_1(v, X, Y);
    1022             :   }
    1023             :   else
    1024             :   { /* not cgetg_copy: x may be a t_VEC */
    1025     6662592 :     long i, lx = lg(X);
    1026     6662592 :     A = cgetg(lx,t_COL);
    1027    43585859 :     for (i=1; i<lx; i++) gel(A,i) = lincombii(u,v,gel(X,i),gel(Y,i));
    1028             :   }
    1029    32420587 :   return A;
    1030             : }
    1031             : 
    1032             : /********************************************************************/
    1033             : /**                                                                **/
    1034             : /**                           CONVERSIONS                          **/
    1035             : /**                                                                **/
    1036             : /********************************************************************/
    1037             : GEN
    1038      400857 : ZV_to_nv(GEN x)
    1039      749174 : { pari_APPLY_ulong(itou(gel(x,i))) }
    1040             : 
    1041             : GEN
    1042      130197 : zm_to_ZM(GEN x)
    1043      815935 : { pari_APPLY_type(t_MAT, zc_to_ZC(gel(x,i))) }
    1044             : 
    1045             : GEN
    1046         126 : zmV_to_ZMV(GEN x)
    1047         791 : { pari_APPLY_type(t_VEC, zm_to_ZM(gel(x,i))) }
    1048             : 
    1049             : /* same as Flm_to_ZM but do not assume positivity */
    1050             : GEN
    1051        1022 : ZM_to_zm(GEN x)
    1052       17472 : { pari_APPLY_same(ZV_to_zv(gel(x,i))) }
    1053             : 
    1054             : GEN
    1055      366646 : zv_to_Flv(GEN x, ulong p)
    1056     5418812 : { pari_APPLY_ulong(umodsu(x[i], p)) }
    1057             : 
    1058             : GEN
    1059       22694 : zm_to_Flm(GEN x, ulong p)
    1060      351750 : { pari_APPLY_same(zv_to_Flv(gel(x,i),p)) }
    1061             : 
    1062             : GEN
    1063          49 : ZMV_to_zmV(GEN x)
    1064         399 : { pari_APPLY_type(t_VEC, ZM_to_zm(gel(x,i))) }
    1065             : 
    1066             : /********************************************************************/
    1067             : /**                                                                **/
    1068             : /**                         COPY, NEGATION                         **/
    1069             : /**                                                                **/
    1070             : /********************************************************************/
    1071             : GEN
    1072    13269746 : ZC_copy(GEN x)
    1073             : {
    1074    13269746 :   long i, lx = lg(x);
    1075    13269746 :   GEN y = cgetg(lx, t_COL);
    1076    83870768 :   for (i=1; i<lx; i++)
    1077             :   {
    1078    70601043 :     GEN c = gel(x,i);
    1079    70601043 :     gel(y,i) = lgefint(c) == 2? gen_0: icopy(c);
    1080             :   }
    1081    13269725 :   return y;
    1082             : }
    1083             : 
    1084             : GEN
    1085      509600 : ZM_copy(GEN x)
    1086     4285920 : { pari_APPLY_same(ZC_copy(gel(x,i))) }
    1087             : 
    1088             : void
    1089      342115 : ZV_neg_inplace(GEN M)
    1090             : {
    1091      342115 :   long l = lg(M);
    1092     1287816 :   while (--l > 0) gel(M,l) = negi(gel(M,l));
    1093      342115 : }
    1094             : GEN
    1095     6264348 : ZC_neg(GEN x)
    1096    36496641 : { pari_APPLY_type(t_COL, negi(gel(x,i))) }
    1097             : 
    1098             : GEN
    1099       51062 : zv_neg(GEN x)
    1100      655284 : { pari_APPLY_long(-x[i]) }
    1101             : GEN
    1102         126 : zv_neg_inplace(GEN M)
    1103             : {
    1104         126 :   long l = lg(M);
    1105         420 :   while (--l > 0) M[l] = -M[l];
    1106         126 :   return M;
    1107             : }
    1108             : GEN
    1109          77 : zv_abs(GEN x)
    1110        5446 : { pari_APPLY_ulong(labs(x[i])) }
    1111             : GEN
    1112     1655047 : ZM_neg(GEN x)
    1113     5039053 : { pari_APPLY_same(ZC_neg(gel(x,i))) }
    1114             : 
    1115             : void
    1116     4921453 : ZV_togglesign(GEN M)
    1117             : {
    1118     4921453 :   long l = lg(M);
    1119    73304754 :   while (--l > 0) togglesign_safe(&gel(M,l));
    1120     4921639 : }
    1121             : void
    1122           0 : ZM_togglesign(GEN M)
    1123             : {
    1124           0 :   long l = lg(M);
    1125           0 :   while (--l > 0) ZV_togglesign(gel(M,l));
    1126           0 : }
    1127             : 
    1128             : /********************************************************************/
    1129             : /**                                                                **/
    1130             : /**                        "DIVISION" mod HNF                      **/
    1131             : /**                                                                **/
    1132             : /********************************************************************/
    1133             : /* Reduce ZC x modulo ZM y in HNF, may return x itself (not a copy) */
    1134             : GEN
    1135    10690478 : ZC_hnfremdiv(GEN x, GEN y, GEN *Q)
    1136             : {
    1137    10690478 :   long i, l = lg(x);
    1138             :   GEN q;
    1139             : 
    1140    10690478 :   if (Q) *Q = cgetg(l,t_COL);
    1141    10690478 :   if (l == 1) return cgetg(1,t_COL);
    1142    61544457 :   for (i = l-1; i>0; i--)
    1143             :   {
    1144    50856895 :     q = diviiround(gel(x,i), gcoeff(y,i,i));
    1145    50855619 :     if (signe(q)) {
    1146    22935874 :       togglesign(q);
    1147    22935924 :       x = ZC_lincomb(gen_1, q, x, gel(y,i));
    1148             :     }
    1149    50853979 :     if (Q) gel(*Q, i) = q;
    1150             :   }
    1151    10687562 :   return x;
    1152             : }
    1153             : 
    1154             : /* x = y Q + R, may return some columns of x (not copies) */
    1155             : GEN
    1156      453829 : ZM_hnfdivrem(GEN x, GEN y, GEN *Q)
    1157             : {
    1158      453829 :   long lx = lg(x), i;
    1159      453829 :   GEN R = cgetg(lx, t_MAT);
    1160      453853 :   if (Q)
    1161             :   {
    1162      127357 :     GEN q = cgetg(lx, t_MAT); *Q = q;
    1163      188193 :     for (i=1; i<lx; i++) gel(R,i) = ZC_hnfremdiv(gel(x,i),y,(GEN*)(q+i));
    1164             :   }
    1165             :   else
    1166      823137 :     for (i=1; i<lx; i++)
    1167             :     {
    1168      496632 :       pari_sp av = avma;
    1169      496632 :       GEN z = ZC_hnfrem(gel(x,i),y);
    1170      496601 :       gel(R,i) = (avma == av)? ZC_copy(z): gerepileupto(av, z);
    1171             :     }
    1172      453862 :   return R;
    1173             : }
    1174             : 
    1175             : /********************************************************************/
    1176             : /**                                                                **/
    1177             : /**                               TESTS                            **/
    1178             : /**                                                                **/
    1179             : /********************************************************************/
    1180             : int
    1181    23103514 : zv_equal0(GEN V)
    1182             : {
    1183    23103514 :   long l = lg(V);
    1184    37579909 :   while (--l > 0)
    1185    30805743 :     if (V[l]) return 0;
    1186     6774166 :   return 1;
    1187             : }
    1188             : 
    1189             : int
    1190    14368004 : ZV_equal0(GEN V)
    1191             : {
    1192    14368004 :   long l = lg(V);
    1193    25399279 :   while (--l > 0)
    1194    24973085 :     if (signe(gel(V,l))) return 0;
    1195      426194 :   return 1;
    1196             : }
    1197             : int
    1198       16409 : ZMrow_equal0(GEN V, long i)
    1199             : {
    1200       16409 :   long l = lg(V);
    1201       25507 :   while (--l > 0)
    1202       21857 :     if (signe(gcoeff(V,i,l))) return 0;
    1203        3650 :   return 1;
    1204             : }
    1205             : 
    1206             : static int
    1207     6258022 : ZV_equal_lg(GEN V, GEN W, long l)
    1208             : {
    1209    25927251 :   while (--l > 0)
    1210    20151121 :     if (!equalii(gel(V,l), gel(W,l))) return 0;
    1211     5776130 :   return 1;
    1212             : }
    1213             : int
    1214      296305 : ZV_equal(GEN V, GEN W)
    1215             : {
    1216      296305 :   long l = lg(V);
    1217      296305 :   if (lg(W) != l) return 0;
    1218      296298 :   return ZV_equal_lg(V, W, l);
    1219             : }
    1220             : int
    1221     3347160 : ZM_equal(GEN A, GEN B)
    1222             : {
    1223     3347160 :   long i, m, l = lg(A);
    1224     3347160 :   if (lg(B) != l) return 0;
    1225     3347106 :   if (l == 1) return 1;
    1226     3347106 :   m = lgcols(A);
    1227     3347106 :   if (lgcols(B) != m) return 0;
    1228     9047771 :   for (i = 1; i < l; i++)
    1229     5961713 :     if (!ZV_equal_lg(gel(A,i), gel(B,i), m)) return 0;
    1230     3086058 :   return 1;
    1231             : }
    1232             : int
    1233       72400 : ZM_equal0(GEN A)
    1234             : {
    1235       72400 :   long i, j, m, l = lg(A);
    1236       72400 :   if (l == 1) return 1;
    1237       72400 :   m = lgcols(A);
    1238      123463 :   for (j = 1; j < l; j++)
    1239     2715086 :     for (i = 1; i < m; i++)
    1240     2664023 :       if (signe(gcoeff(A,i,j))) return 0;
    1241       15585 :   return 1;
    1242             : }
    1243             : int
    1244    30837903 : zv_equal(GEN V, GEN W)
    1245             : {
    1246    30837903 :   long l = lg(V);
    1247    30837903 :   if (lg(W) != l) return 0;
    1248   262401726 :   while (--l > 0)
    1249   232681666 :     if (V[l] != W[l]) return 0;
    1250    29720060 :   return 1;
    1251             : }
    1252             : 
    1253             : int
    1254        1638 : zvV_equal(GEN V, GEN W)
    1255             : {
    1256        1638 :   long l = lg(V);
    1257        1638 :   if (lg(W) != l) return 0;
    1258       80388 :   while (--l > 0)
    1259       79912 :     if (!zv_equal(gel(V,l),gel(W,l))) return 0;
    1260         476 :   return 1;
    1261             : }
    1262             : 
    1263             : int
    1264      759274 : ZM_ishnf(GEN x)
    1265             : {
    1266      759274 :   long i,j, lx = lg(x);
    1267     2282269 :   for (i=1; i<lx; i++)
    1268             :   {
    1269     1635733 :     GEN xii = gcoeff(x,i,i);
    1270     1635733 :     if (signe(xii) <= 0) return 0;
    1271     3175997 :     for (j=1; j<i; j++)
    1272     1577665 :       if (signe(gcoeff(x,i,j))) return 0;
    1273     3234507 :     for (j=i+1; j<lx; j++)
    1274             :     {
    1275     1711512 :       GEN xij = gcoeff(x,i,j);
    1276     1711512 :       if (signe(xij)<0 || cmpii(xij,xii)>=0) return 0;
    1277             :     }
    1278             :   }
    1279      646536 :   return 1;
    1280             : }
    1281             : int
    1282      637983 : ZM_isidentity(GEN x)
    1283             : {
    1284      637983 :   long i,j, lx = lg(x);
    1285             : 
    1286      637983 :   if (lx == 1) return 1;
    1287      637976 :   if (lx != lgcols(x)) return 0;
    1288     3116697 :   for (j=1; j<lx; j++)
    1289             :   {
    1290     2480779 :     GEN c = gel(x,j);
    1291     7694193 :     for (i=1; i<j; )
    1292     5213450 :       if (signe(gel(c,i++))) return 0;
    1293             :     /* i = j */
    1294     2480743 :     if (!equali1(gel(c,i++))) return 0;
    1295     7700263 :     for (   ; i<lx; )
    1296     5221528 :       if (signe(gel(c,i++))) return 0;
    1297             :   }
    1298      635918 :   return 1;
    1299             : }
    1300             : int
    1301      556029 : ZM_isdiagonal(GEN x)
    1302             : {
    1303      556029 :   long i,j, lx = lg(x);
    1304      556029 :   if (lx == 1) return 1;
    1305      556029 :   if (lx != lgcols(x)) return 0;
    1306             : 
    1307     1437423 :   for (j=1; j<lx; j++)
    1308             :   {
    1309     1207916 :     GEN c = gel(x,j);
    1310     1695722 :     for (i=1; i<j; i++)
    1311      814290 :       if (signe(gel(c,i))) return 0;
    1312     2018971 :     for (i++; i<lx; i++)
    1313     1137574 :       if (signe(gel(c,i))) return 0;
    1314             :   }
    1315      229507 :   return 1;
    1316             : }
    1317             : int
    1318      103974 : ZM_isscalar(GEN x, GEN s)
    1319             : {
    1320      103974 :   long i, j, lx = lg(x);
    1321             : 
    1322      103974 :   if (lx == 1) return 1;
    1323      103974 :   if (!s) s = gcoeff(x,1,1);
    1324      103974 :   if (equali1(s)) return ZM_isidentity(x);
    1325      102763 :   if (lx != lgcols(x)) return 0;
    1326      547692 :   for (j=1; j<lx; j++)
    1327             :   {
    1328      456924 :     GEN c = gel(x,j);
    1329     2251292 :     for (i=1; i<j; )
    1330     1804373 :       if (signe(gel(c,i++))) return 0;
    1331             :     /* i = j */
    1332      446919 :     if (!equalii(gel(c,i++), s)) return 0;
    1333     2278882 :     for (   ; i<lx; )
    1334     1833953 :       if (signe(gel(c,i++))) return 0;
    1335             :   }
    1336       90768 :   return 1;
    1337             : }
    1338             : 
    1339             : long
    1340      106942 : ZC_is_ei(GEN x)
    1341             : {
    1342      106942 :   long i, j = 0, l = lg(x);
    1343      990877 :   for (i = 1; i < l; i++)
    1344             :   {
    1345      883936 :     GEN c = gel(x,i);
    1346      883936 :     long s = signe(c);
    1347      883936 :     if (!s) continue;
    1348      106936 :     if (s < 0 || !is_pm1(c) || j) return 0;
    1349      106935 :     j = i;
    1350             :   }
    1351      106941 :   return j;
    1352             : }
    1353             : 
    1354             : /********************************************************************/
    1355             : /**                                                                **/
    1356             : /**                       MISCELLANEOUS                            **/
    1357             : /**                                                                **/
    1358             : /********************************************************************/
    1359             : /* assume lg(x) = lg(y), x,y in Z^n */
    1360             : int
    1361     3133085 : ZV_cmp(GEN x, GEN y)
    1362             : {
    1363     3133085 :   long fl,i, lx = lg(x);
    1364     6344748 :   for (i=1; i<lx; i++)
    1365     5057580 :     if (( fl = cmpii(gel(x,i), gel(y,i)) )) return fl;
    1366     1287168 :   return 0;
    1367             : }
    1368             : /* assume lg(x) = lg(y), x,y in Z^n */
    1369             : int
    1370       19740 : ZV_abscmp(GEN x, GEN y)
    1371             : {
    1372       19740 :   long fl,i, lx = lg(x);
    1373       54049 :   for (i=1; i<lx; i++)
    1374       53922 :     if (( fl = abscmpii(gel(x,i), gel(y,i)) )) return fl;
    1375         127 :   return 0;
    1376             : }
    1377             : 
    1378             : long
    1379    20720058 : zv_content(GEN x)
    1380             : {
    1381    20720058 :   long i, s, l = lg(x);
    1382    20720058 :   if (l == 1) return 0;
    1383    20720051 :   s = labs(x[1]);
    1384    46519862 :   for (i = 2; i < l && s != 1; i++) s = ugcd(s, labs(x[i]));
    1385    20721407 :   return s;
    1386             : }
    1387             : GEN
    1388      297322 : ZV_content(GEN x)
    1389             : {
    1390      297322 :   long i, l = lg(x);
    1391      297322 :   pari_sp av = avma;
    1392             :   GEN c;
    1393      297322 :   if (l == 1) return gen_0;
    1394      297322 :   if (l == 2) return absi(gel(x,1));
    1395      204796 :   c = gel(x,1);
    1396      557875 :   for (i = 2; i < l; i++)
    1397             :   {
    1398      403575 :     c = gcdii(c, gel(x,i));
    1399      403575 :     if (is_pm1(c)) { set_avma(av); return gen_1; }
    1400             :   }
    1401      154300 :   return gerepileuptoint(av, c);
    1402             : }
    1403             : 
    1404             : GEN
    1405     3868429 : ZM_det_triangular(GEN mat)
    1406             : {
    1407             :   pari_sp av;
    1408     3868429 :   long i,l = lg(mat);
    1409             :   GEN s;
    1410             : 
    1411     3868429 :   if (l<3) return l<2? gen_1: icopy(gcoeff(mat,1,1));
    1412     3465350 :   av = avma; s = gcoeff(mat,1,1);
    1413     9398958 :   for (i=2; i<l; i++) s = mulii(s,gcoeff(mat,i,i));
    1414     3464996 :   return gerepileuptoint(av,s);
    1415             : }
    1416             : 
    1417             : /* assumes no overflow */
    1418             : long
    1419      949867 : zv_prod(GEN v)
    1420             : {
    1421      949867 :   long n, i, l = lg(v);
    1422      949867 :   if (l == 1) return 1;
    1423      959951 :   n = v[1]; for (i = 2; i < l; i++) n *= v[i];
    1424      771535 :   return n;
    1425             : }
    1426             : 
    1427             : static GEN
    1428   318095621 : _mulii(void *E, GEN a, GEN b)
    1429   318095621 : { (void) E; return mulii(a, b); }
    1430             : 
    1431             : /* product of ulongs */
    1432             : GEN
    1433     1864321 : zv_prod_Z(GEN v)
    1434             : {
    1435             :   pari_sp av;
    1436     1864321 :   long k, m, n = lg(v)-1;
    1437     1864321 :   int stop = 0;
    1438             :   GEN V;
    1439     1864321 :   switch(n) {
    1440       20902 :     case 0: return gen_1;
    1441      129878 :     case 1: return utoi(v[1]);
    1442      976693 :     case 2: return muluu(v[1], v[2]);
    1443             :   }
    1444      736848 :   av = avma; m = n >> 1;
    1445      736848 :   V = cgetg(m + (odd(n)? 2: 1), t_VEC);
    1446   151571184 :   for (k = n; k; k--) /* start from the end: v is usually sorted */
    1447   150835379 :     if (v[k] & HIGHMASK) { stop = 1; break; }
    1448     2422557 :   while (!stop)
    1449             :   { /* HACK: handle V as a t_VECSMALL; gain a few iterations */
    1450    81551547 :     for (k = 1; k <= m; k++)
    1451             :     {
    1452    79256705 :       V[k] = uel(v,k<<1) * uel(v,(k<<1)-1);
    1453    79256705 :       if (V[k] & HIGHMASK) stop = 1; /* last "free" iteration */
    1454             :     }
    1455     2294842 :     if (odd(n))
    1456             :     {
    1457     1352286 :       if (n == 1) { set_avma(av); return utoi(v[1]); }
    1458      743152 :       V[++m] = v[n];
    1459             :     }
    1460     1685710 :     v = V; n = m; m = n >> 1;
    1461             :   }
    1462             :   /* n > 1; m > 0 */
    1463      127715 :   if (n == 2) { set_avma(av); return muluu(v[1], v[2]); }
    1464    46319407 :   for (k = 1; k <= m; k++) gel(V,k) = muluu(v[k<<1], v[(k<<1)-1]);
    1465       87482 :   if (odd(n)) gel(V, ++m) = utoipos(v[n]);
    1466       87480 :   setlg(V, m+1); /* HACK: now V is a bona fide t_VEC */
    1467       87480 :   return gerepileuptoint(av, gen_product(V, NULL, &_mulii));
    1468             : }
    1469             : GEN
    1470    14694393 : vecsmall_prod(GEN v)
    1471             : {
    1472    14694393 :   pari_sp av = avma;
    1473    14694393 :   long k, m, n = lg(v)-1;
    1474             :   GEN V;
    1475    14694393 :   switch (n) {
    1476           0 :     case 0: return gen_1;
    1477           0 :     case 1: return stoi(v[1]);
    1478          21 :     case 2: return mulss(v[1], v[2]);
    1479             :   }
    1480    14694372 :   m = n >> 1;
    1481    14694372 :   V = cgetg(m + (odd(n)? 2: 1), t_VEC);
    1482   161556906 :   for (k = 1; k <= m; k++) gel(V,k) = mulss(v[k<<1], v[(k<<1)-1]);
    1483    14694372 :   if (odd(n)) gel(V,k) = stoi(v[n]);
    1484    14694372 :   return gerepileuptoint(av, gen_product(V, NULL, &_mulii));
    1485             : }
    1486             : 
    1487             : GEN
    1488     8819005 : ZV_prod(GEN v)
    1489             : {
    1490     8819005 :   pari_sp av = avma;
    1491     8819005 :   long i, l = lg(v);
    1492             :   GEN n;
    1493     8819005 :   if (l == 1) return gen_1;
    1494     8635057 :   if (l > 7) return gerepileuptoint(av, gen_product(v, NULL, _mulii));
    1495     1292281 :   n = gel(v,1);
    1496     1292281 :   if (l == 2) return icopy(n);
    1497     2090993 :   for (i = 2; i < l; i++) n = mulii(n, gel(v,i));
    1498      840320 :   return gerepileuptoint(av, n);
    1499             : }
    1500             : /* assumes no overflow */
    1501             : long
    1502       16445 : zv_sum(GEN v)
    1503             : {
    1504       16445 :   long n, i, l = lg(v);
    1505       16445 :   if (l == 1) return 0;
    1506       95742 :   n = v[1]; for (i = 2; i < l; i++) n += v[i];
    1507       16424 :   return n;
    1508             : }
    1509             : /* assumes no overflow and 0 <= n <= #v */
    1510             : long
    1511           0 : zv_sumpart(GEN v, long n)
    1512             : {
    1513             :   long i, p;
    1514           0 :   if (!n) return 0;
    1515           0 :   p = v[1]; for (i = 2; i <= n; i++) p += v[i];
    1516           0 :   return p;
    1517             : }
    1518             : GEN
    1519          77 : ZV_sum(GEN v)
    1520             : {
    1521          77 :   pari_sp av = avma;
    1522          77 :   long i, l = lg(v);
    1523             :   GEN n;
    1524          77 :   if (l == 1) return gen_0;
    1525          77 :   n = gel(v,1);
    1526          77 :   if (l == 2) return icopy(n);
    1527         581 :   for (i = 2; i < l; i++) n = addii(n, gel(v,i));
    1528          77 :   return gerepileuptoint(av, n);
    1529             : }
    1530             : 
    1531             : /********************************************************************/
    1532             : /**                                                                **/
    1533             : /**         GRAM SCHMIDT REDUCTION (integer matrices)              **/
    1534             : /**                                                                **/
    1535             : /********************************************************************/
    1536             : 
    1537             : /* L[k,] += q * L[l,], l < k. Inefficient if q = 0 */
    1538             : static void
    1539      362302 : Zupdate_row(long k, long l, GEN q, GEN L, GEN B)
    1540             : {
    1541      362302 :   long i, qq = itos_or_0(q);
    1542      362302 :   if (!qq)
    1543             :   {
    1544       33377 :     for(i=1;i<l;i++)  gcoeff(L,k,i) = addii(gcoeff(L,k,i),mulii(q,gcoeff(L,l,i)));
    1545        7069 :     gcoeff(L,k,l) = addii(gcoeff(L,k,l), mulii(q,B));
    1546        7069 :     return;
    1547             :   }
    1548      355233 :   if (qq == 1) {
    1549      148802 :     for (i=1;i<l; i++) gcoeff(L,k,i) = addii(gcoeff(L,k,i),gcoeff(L,l,i));
    1550      102094 :     gcoeff(L,k,l) = addii(gcoeff(L,k,l), B);
    1551      253139 :   } else if (qq == -1) {
    1552      151505 :     for (i=1;i<l; i++) gcoeff(L,k,i) = subii(gcoeff(L,k,i),gcoeff(L,l,i));
    1553       89220 :     gcoeff(L,k,l) = addii(gcoeff(L,k,l), negi(B));
    1554             :   } else {
    1555      289622 :     for(i=1;i<l;i++) gcoeff(L,k,i) = addii(gcoeff(L,k,i),mulsi(qq,gcoeff(L,l,i)));
    1556      163930 :     gcoeff(L,k,l) = addii(gcoeff(L,k,l), mulsi(qq,B));
    1557             :   }
    1558             : }
    1559             : 
    1560             : /* update L[k,] */
    1561             : static void
    1562     1033405 : ZRED(long k, long l, GEN x, GEN L, GEN B)
    1563             : {
    1564     1033405 :   GEN q = truedivii(addii(B,shifti(gcoeff(L,k,l),1)), shifti(B,1));
    1565     1033522 :   if (!signe(q)) return;
    1566      362282 :   q = negi(q);
    1567      362297 :   Zupdate_row(k,l,q,L,B);
    1568      362282 :   gel(x,k) = ZC_lincomb(gen_1, q, gel(x,k), gel(x,l));
    1569             : }
    1570             : 
    1571             : /* Gram-Schmidt reduction, x a ZM */
    1572             : static void
    1573     1183973 : ZincrementalGS(GEN x, GEN L, GEN B, long k)
    1574             : {
    1575             :   long i, j;
    1576     3780708 :   for (j=1; j<=k; j++)
    1577             :   {
    1578     2597505 :     pari_sp av = avma;
    1579     2597505 :     GEN u = ZV_dotproduct(gel(x,k), gel(x,j));
    1580     5613344 :     for (i=1; i<j; i++)
    1581             :     {
    1582     3017230 :       u = subii(mulii(gel(B,i+1), u), mulii(gcoeff(L,k,i), gcoeff(L,j,i)));
    1583     3016601 :       u = diviiexact(u, gel(B,i));
    1584             :     }
    1585     2596114 :     gcoeff(L,k,j) = gerepileuptoint(av, u);
    1586             :   }
    1587     1183203 :   gel(B,k+1) = gcoeff(L,k,k); gcoeff(L,k,k) = gen_1;
    1588     1183203 : }
    1589             : 
    1590             : /* Variant reducemodinvertible(ZC v, ZM y), when y singular.
    1591             :  * Very inefficient if y is not LLL-reduced of maximal rank */
    1592             : static GEN
    1593      110485 : ZC_reducemodmatrix_i(GEN v, GEN y)
    1594             : {
    1595      110485 :   GEN B, L, x = shallowconcat(y, v);
    1596      110485 :   long k, lx = lg(x), nx = lx-1;
    1597             : 
    1598      110485 :   B = scalarcol_shallow(gen_1, lx);
    1599      110484 :   L = zeromatcopy(nx, nx);
    1600      454520 :   for (k=1; k <= nx; k++) ZincrementalGS(x, L, B, k);
    1601      344039 :   for (k = nx-1; k >= 1; k--) ZRED(nx,k, x,L,gel(B,k+1));
    1602      110483 :   return gel(x,nx);
    1603             : }
    1604             : GEN
    1605      110485 : ZC_reducemodmatrix(GEN v, GEN y) {
    1606      110485 :   pari_sp av = avma;
    1607      110485 :   return gerepilecopy(av, ZC_reducemodmatrix_i(v,y));
    1608             : }
    1609             : 
    1610             : /* Variant reducemodinvertible(ZM v, ZM y), when y singular.
    1611             :  * Very inefficient if y is not LLL-reduced of maximal rank */
    1612             : static GEN
    1613      226915 : ZM_reducemodmatrix_i(GEN v, GEN y)
    1614             : {
    1615             :   GEN B, L, V;
    1616      226915 :   long j, k, lv = lg(v), nx = lg(y), lx = nx+1;
    1617             : 
    1618      226915 :   V = cgetg(lv, t_MAT);
    1619      226926 :   B = scalarcol_shallow(gen_1, lx);
    1620      226926 :   L = zeromatcopy(nx, nx);
    1621      601046 :   for (k=1; k < nx; k++) ZincrementalGS(y, L, B, k);
    1622      692745 :   for (j = 1; j < lg(v); j++)
    1623             :   {
    1624      465853 :     GEN x = shallowconcat(y, gel(v,j));
    1625      465864 :     ZincrementalGS(x, L, B, nx); /* overwrite last */
    1626     1265773 :     for (k = nx-1; k >= 1; k--) ZRED(nx,k, x,L,gel(B,k+1));
    1627      465842 :     gel(V,j) = gel(x,nx);
    1628             :   }
    1629      226892 :   return V;
    1630             : }
    1631             : GEN
    1632      226914 : ZM_reducemodmatrix(GEN v, GEN y) {
    1633      226914 :   pari_sp av = avma;
    1634      226914 :   return gerepilecopy(av, ZM_reducemodmatrix_i(v,y));
    1635             : }
    1636             : 
    1637             : GEN
    1638       98543 : ZC_reducemodlll(GEN x,GEN y)
    1639             : {
    1640       98543 :   pari_sp av = avma;
    1641       98543 :   GEN z = ZC_reducemodmatrix(x, ZM_lll(y, 0.75, LLL_INPLACE));
    1642       98545 :   return gerepilecopy(av, z);
    1643             : }
    1644             : GEN
    1645           0 : ZM_reducemodlll(GEN x,GEN y)
    1646             : {
    1647           0 :   pari_sp av = avma;
    1648           0 :   GEN z = ZM_reducemodmatrix(x, ZM_lll(y, 0.75, LLL_INPLACE));
    1649           0 :   return gerepilecopy(av, z);
    1650             : }

Generated by: LCOV version 1.16