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 29395-ef22f77854) Lines: 791 900 87.9 %
Date: 2024-06-14 09:03:06 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     1807287 : check_ZV(GEN x, long l)
      20             : {
      21             :   long i;
      22    10566509 :   for (i=1; i<l; i++)
      23     8759313 :     if (typ(gel(x,i)) != t_INT) return 0;
      24     1807196 :   return 1;
      25             : }
      26             : void
      27     1422197 : RgV_check_ZV(GEN A, const char *s)
      28             : {
      29     1422197 :   if (!RgV_is_ZV(A)) pari_err_TYPE(stack_strcat(s," [integer vector]"), A);
      30     1422190 : }
      31             : void
      32      640272 : RgM_check_ZM(GEN A, const char *s)
      33             : {
      34      640272 :   long n = lg(A);
      35      640272 :   if (n != 1)
      36             :   {
      37      640125 :     long j, m = lgcols(A);
      38     2447321 :     for (j=1; j<n; j++)
      39     1807287 :       if (!check_ZV(gel(A,j), m))
      40          91 :         pari_err_TYPE(stack_strcat(s," [integer matrix]"), A);
      41             :   }
      42      640181 : }
      43             : 
      44             : static long
      45   109496734 : ZV_max_lg_i(GEN x, long m)
      46             : {
      47   109496734 :   long i, l = 2;
      48   868791156 :   for (i = 1; i < m; i++) l = maxss(l, lgefint(gel(x,i)));
      49   109497087 :   return l;
      50             : }
      51             : long
      52       10640 : ZV_max_lg(GEN x) { return ZV_max_lg_i(x, lg(x)); }
      53             : 
      54             : static long
      55    28650840 : ZM_max_lg_i(GEN x, long n, long m)
      56             : {
      57    28650840 :   long j, l = 2;
      58   138136591 :   for (j = 1; j < n; j++) l = maxss(l, ZV_max_lg_i(gel(x,j), m));
      59    28651287 :   return l;
      60             : }
      61             : long
      62       21647 : ZM_max_lg(GEN x)
      63             : {
      64       21647 :   long n = lg(x);
      65       21647 :   if (n == 1) return 2;
      66       21647 :   return ZM_max_lg_i(x, n, lgcols(x));
      67             : }
      68             : 
      69             : static long
      70           0 : ZV_max_expi_i(GEN x, long m)
      71             : {
      72           0 :   long i, prec = 2;
      73           0 :   for (i = 1; i < m; i++) prec = maxss(prec, expi(gel(x,i)));
      74           0 :   return prec;
      75             : }
      76             : long
      77           0 : ZV_max_expi(GEN x) { return ZV_max_expi_i(x, lg(x)); }
      78             : 
      79             : static long
      80           0 : ZM_max_expi_i(GEN x, long n, long m)
      81             : {
      82           0 :   long j, prec = 2;
      83           0 :   for (j = 1; j < n; j++) prec = maxss(prec, ZV_max_expi_i(gel(x,j), m));
      84           0 :   return prec;
      85             : }
      86             : long
      87           0 : ZM_max_expi(GEN x)
      88             : {
      89           0 :   long n = lg(x);
      90           0 :   if (n == 1) return 2;
      91           0 :   return ZM_max_expi_i(x, n, lgcols(x));
      92             : }
      93             : 
      94             : GEN
      95        4001 : ZM_supnorm(GEN x)
      96             : {
      97        4001 :   long i, j, h, lx = lg(x);
      98        4001 :   GEN s = gen_0;
      99        4001 :   if (lx == 1) return gen_1;
     100        4001 :   h = lgcols(x);
     101       24408 :   for (j=1; j<lx; j++)
     102             :   {
     103       20407 :     GEN xj = gel(x,j);
     104      281690 :     for (i=1; i<h; i++)
     105             :     {
     106      261283 :       GEN c = gel(xj,i);
     107      261283 :       if (abscmpii(c, s) > 0) s = c;
     108             :     }
     109             :   }
     110        4001 :   return absi(s);
     111             : }
     112             : 
     113             : /********************************************************************/
     114             : /**                                                                **/
     115             : /**                           MULTIPLICATION                       **/
     116             : /**                                                                **/
     117             : /********************************************************************/
     118             : /* x nonempty ZM, y a compatible nc (dimension > 0). */
     119             : static GEN
     120           0 : ZM_nc_mul_i(GEN x, GEN y, long c, long l)
     121             : {
     122             :   long i, j;
     123             :   pari_sp av;
     124           0 :   GEN z = cgetg(l,t_COL), s;
     125             : 
     126           0 :   for (i=1; i<l; i++)
     127             :   {
     128           0 :     av = avma; s = muliu(gcoeff(x,i,1),y[1]);
     129           0 :     for (j=2; j<c; j++)
     130           0 :       if (y[j]) s = addii(s, muliu(gcoeff(x,i,j),y[j]));
     131           0 :     gel(z,i) = gerepileuptoint(av,s);
     132             :   }
     133           0 :   return z;
     134             : }
     135             : 
     136             : /* x ZV, y a compatible zc. */
     137             : GEN
     138     2229528 : ZV_zc_mul(GEN x, GEN y)
     139             : {
     140     2229528 :   long j, l = lg(x);
     141     2229528 :   pari_sp av = avma;
     142     2229528 :   GEN s = mulis(gel(x,1),y[1]);
     143    50292998 :   for (j=2; j<l; j++)
     144    48063470 :     if (y[j]) s = addii(s, mulis(gel(x,j),y[j]));
     145     2229528 :   return gerepileuptoint(av,s);
     146             : }
     147             : 
     148             : /* x nonempty ZM, y a compatible zc (dimension > 0). */
     149             : static GEN
     150    20982903 : ZM_zc_mul_i(GEN x, GEN y, long c, long l)
     151             : {
     152             :   long i, j;
     153    20982903 :   GEN z = cgetg(l,t_COL);
     154             : 
     155   126932968 :   for (i=1; i<l; i++)
     156             :   {
     157   105956351 :     pari_sp av = avma;
     158   105956351 :     GEN s = mulis(gcoeff(x,i,1),y[1]);
     159  1178804364 :     for (j=2; j<c; j++)
     160  1072852465 :       if (y[j]) s = addii(s, mulis(gcoeff(x,i,j),y[j]));
     161   105951899 :     gel(z,i) = gerepileuptoint(av,s);
     162             :   }
     163    20976617 :   return z;
     164             : }
     165             : GEN
     166    18989381 : ZM_zc_mul(GEN x, GEN y) {
     167    18989381 :   long l = lg(x);
     168    18989381 :   if (l == 1) return cgetg(1, t_COL);
     169    18989381 :   return ZM_zc_mul_i(x,y, l, lgcols(x));
     170             : }
     171             : 
     172             : /* y nonempty ZM, x a compatible zv (dimension > 0). */
     173             : GEN
     174        1736 : zv_ZM_mul(GEN x, GEN y) {
     175        1736 :   long i,j, lx = lg(x), ly = lg(y);
     176             :   GEN z;
     177        1736 :   if (lx == 1) return zerovec(ly-1);
     178        1736 :   z = cgetg(ly,t_VEC);
     179        4046 :   for (j=1; j<ly; j++)
     180             :   {
     181        2310 :     pari_sp av = avma;
     182        2310 :     GEN s = mulsi(x[1], gcoeff(y,1,j));
     183        3990 :     for (i=2; i<lx; i++)
     184        1680 :       if (x[i]) s = addii(s, mulsi(x[i], gcoeff(y,i,j)));
     185        2310 :     gel(z,j) = gerepileuptoint(av,s);
     186             :   }
     187        1736 :   return z;
     188             : }
     189             : 
     190             : /* x ZM, y a compatible zm (dimension > 0). */
     191             : GEN
     192      949053 : ZM_zm_mul(GEN x, GEN y)
     193             : {
     194      949053 :   long j, c, l = lg(x), ly = lg(y);
     195      949053 :   GEN z = cgetg(ly, t_MAT);
     196      949053 :   if (l == 1) return z;
     197      949046 :   c = lgcols(x);
     198     2942510 :   for (j = 1; j < ly; j++) gel(z,j) = ZM_zc_mul_i(x, gel(y,j), l,c);
     199      949046 :   return z;
     200             : }
     201             : /* x ZM, y a compatible zn (dimension > 0). */
     202             : GEN
     203           0 : ZM_nm_mul(GEN x, GEN y)
     204             : {
     205           0 :   long j, c, l = lg(x), ly = lg(y);
     206           0 :   GEN z = cgetg(ly, t_MAT);
     207           0 :   if (l == 1) return z;
     208           0 :   c = lgcols(x);
     209           0 :   for (j = 1; j < ly; j++) gel(z,j) = ZM_nc_mul_i(x, gel(y,j), l,c);
     210           0 :   return z;
     211             : }
     212             : 
     213             : /* Strassen-Winograd algorithm */
     214             : 
     215             : /* Return A[ma+1..ma+da, na+1..na+ea] - B[mb+1..mb+db, nb+1..nb+eb]
     216             :  * as an (m x n)-matrix, padding the input with zeroes as necessary. */
     217             : static GEN
     218      338136 : add_slices(long m, long n,
     219             :            GEN A, long ma, long da, long na, long ea,
     220             :            GEN B, long mb, long db, long nb, long eb)
     221             : {
     222      338136 :   long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
     223      338136 :   GEN M = cgetg(n + 1, t_MAT), C;
     224             : 
     225     2689689 :   for (j = 1; j <= min_e; j++) {
     226     2351553 :     gel(M, j) = C = cgetg(m + 1, t_COL);
     227    39340584 :     for (i = 1; i <= min_d; i++)
     228    36989031 :       gel(C, i) = addii(gcoeff(A, ma + i, na + j),
     229    36989031 :                         gcoeff(B, mb + i, nb + j));
     230     2386658 :     for (; i <= da; i++) gel(C, i) = gcoeff(A, ma + i, na + j);
     231     2351553 :     for (; i <= db; i++) gel(C, i) = gcoeff(B, mb + i, nb + j);
     232     2351553 :     for (; i <= m; i++)  gel(C, i) = gen_0;
     233             :   }
     234      368757 :   for (; j <= ea; j++) {
     235       30621 :     gel(M, j) = C = cgetg(m + 1, t_COL);
     236      116158 :     for (i = 1; i <= da; i++) gel(C, i) = gcoeff(A, ma + i, na + j);
     237       30621 :     for (; i <= m; i++) gel(C, i) = gen_0;
     238             :   }
     239      338136 :   for (; j <= eb; j++) {
     240           0 :     gel(M, j) = C = cgetg(m + 1, t_COL);
     241           0 :     for (i = 1; i <= db; i++) gel(C, i) = gcoeff(B, mb + i, nb + j);
     242           0 :     for (; i <= m; i++) gel(C, i) = gen_0;
     243             :   }
     244      338136 :   for (; j <= n; j++) gel(M, j) = zerocol(m);
     245      338136 :   return M;
     246             : }
     247             : 
     248             : /* Return A[ma+1..ma+da, na+1..na+ea] - B[mb+1..mb+db, nb+1..nb+eb]
     249             :  * as an (m x n)-matrix, padding the input with zeroes as necessary. */
     250             : static GEN
     251      295869 : subtract_slices(long m, long n,
     252             :                 GEN A, long ma, long da, long na, long ea,
     253             :                 GEN B, long mb, long db, long nb, long eb)
     254             : {
     255      295869 :   long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
     256      295869 :   GEN M = cgetg(n + 1, t_MAT), C;
     257             : 
     258     2344677 :   for (j = 1; j <= min_e; j++) {
     259     2048808 :     gel(M, j) = C = cgetg(m + 1, t_COL);
     260    35116260 :     for (i = 1; i <= min_d; i++)
     261    33067452 :       gel(C, i) = subii(gcoeff(A, ma + i, na + j),
     262    33067452 :                         gcoeff(B, mb + i, nb + j));
     263     2074925 :     for (; i <= da; i++) gel(C, i) = gcoeff(A, ma + i, na + j);
     264     2091134 :     for (; i <= db; i++) gel(C, i) = negi(gcoeff(B, mb + i, nb + j));
     265     2048808 :     for (; i <= m; i++) gel(C, i) = gen_0;
     266             :   }
     267      295869 :   for (; j <= ea; j++) {
     268           0 :     gel(M, j) = C = cgetg(m + 1, t_COL);
     269           0 :     for (i = 1; i <= da; i++) gel(C, i) = gcoeff(A, ma + i, na + j);
     270           0 :     for (; i <= m; i++) gel(C, i) = gen_0;
     271             :   }
     272      316550 :   for (; j <= eb; j++) {
     273       20681 :     gel(M, j) = C = cgetg(m + 1, t_COL);
     274       66473 :     for (i = 1; i <= db; i++) gel(C, i) = negi(gcoeff(B, mb + i, nb + j));
     275       20681 :     for (; i <= m; i++) gel(C, i) = gen_0;
     276             :   }
     277      316550 :   for (; j <= n; j++) gel(M, j) = zerocol(m);
     278      295869 :   return M;
     279             : }
     280             : 
     281             : static GEN ZM_mul_i(GEN x, GEN y, long l, long lx, long ly);
     282             : 
     283             : /* Strassen-Winograd matrix product A (m x n) * B (n x p) */
     284             : static GEN
     285       42267 : ZM_mul_sw(GEN A, GEN B, long m, long n, long p)
     286             : {
     287       42267 :   pari_sp av = avma;
     288       42267 :   long m1 = (m + 1)/2, m2 = m/2,
     289       42267 :     n1 = (n + 1)/2, n2 = n/2,
     290       42267 :     p1 = (p + 1)/2, p2 = p/2;
     291             :   GEN A11, A12, A22, B11, B21, B22,
     292             :     S1, S2, S3, S4, T1, T2, T3, T4,
     293             :     M1, M2, M3, M4, M5, M6, M7,
     294             :     V1, V2, V3, C11, C12, C21, C22, C;
     295             : 
     296       42267 :   T2 = subtract_slices(n1, p2, B, 0, n1, p1, p2, B, n1, n2, p1, p2);
     297       42267 :   S1 = subtract_slices(m2, n1, A, m1, m2, 0, n1, A, 0, m2, 0, n1);
     298       42267 :   M2 = ZM_mul_i(S1, T2, m2 + 1, n1 + 1, p2 + 1);
     299       42267 :   if (gc_needed(av, 1))
     300           0 :     gerepileall(av, 2, &T2, &M2);  /* destroy S1 */
     301       42267 :   T3 = subtract_slices(n1, p1, T2, 0, n1, 0, p2, B, 0, n1, 0, p1);
     302       42267 :   if (gc_needed(av, 1))
     303           0 :     gerepileall(av, 2, &M2, &T3);  /* destroy T2 */
     304       42267 :   S2 = add_slices(m2, n1, A, m1, m2, 0, n1, A, m1, m2, n1, n2);
     305       42267 :   T1 = subtract_slices(n1, p1, B, 0, n1, p1, p2, B, 0, n1, 0, p2);
     306       42267 :   M3 = ZM_mul_i(S2, T1, m2 + 1, n1 + 1, p2 + 1);
     307       42267 :   if (gc_needed(av, 1))
     308           0 :     gerepileall(av, 4, &M2, &T3, &S2, &M3);  /* destroy T1 */
     309       42267 :   S3 = subtract_slices(m1, n1, S2, 0, m2, 0, n1, A, 0, m1, 0, n1);
     310       42267 :   if (gc_needed(av, 1))
     311           0 :     gerepileall(av, 4, &M2, &T3, &M3, &S3);  /* destroy S2 */
     312       42267 :   A11 = matslice(A, 1, m1, 1, n1);
     313       42267 :   B11 = matslice(B, 1, n1, 1, p1);
     314       42267 :   M1 = ZM_mul_i(A11, B11, m1 + 1, n1 + 1, p1 + 1);
     315       42267 :   if (gc_needed(av, 1))
     316           0 :     gerepileall(av, 5, &M2, &T3, &M3, &S3, &M1);  /* destroy A11, B11 */
     317       42267 :   A12 = matslice(A, 1, m1, n1 + 1, n);
     318       42267 :   B21 = matslice(B, n1 + 1, n, 1, p1);
     319       42267 :   M4 = ZM_mul_i(A12, B21, m1 + 1, n2 + 1, p1 + 1);
     320       42267 :   if (gc_needed(av, 1))
     321           0 :     gerepileall(av, 6, &M2, &T3, &M3, &S3, &M1, &M4);  /* destroy A12, B21 */
     322       42267 :   C11 = add_slices(m1, p1, M1, 0, m1, 0, p1, M4, 0, m1, 0, p1);
     323       42267 :   if (gc_needed(av, 1))
     324           0 :     gerepileall(av, 6, &M2, &T3, &M3, &S3, &M1, &C11);  /* destroy M4 */
     325       42267 :   M5 = ZM_mul_i(S3, T3, m1 + 1, n1 + 1, p1 + 1);
     326       42267 :   S4 = subtract_slices(m1, n2, A, 0, m1, n1, n2, S3, 0, m1, 0, n2);
     327       42267 :   if (gc_needed(av, 1))
     328           5 :     gerepileall(av, 7, &M2, &T3, &M3, &M1, &C11, &M5, &S4);  /* destroy S3 */
     329       42267 :   T4 = add_slices(n2, p1, B, n1, n2, 0, p1, T3, 0, n2, 0, p1);
     330       42267 :   if (gc_needed(av, 1))
     331           0 :     gerepileall(av, 7, &M2, &M3, &M1, &C11, &M5, &S4, &T4);  /* destroy T3 */
     332       42267 :   V1 = subtract_slices(m1, p1, M1, 0, m1, 0, p1, M5, 0, m1, 0, p1);
     333       42267 :   if (gc_needed(av, 1))
     334           1 :     gerepileall(av, 6, &M2, &M3, &S4, &T4, &C11, &V1);  /* destroy M1, M5 */
     335       42267 :   B22 = matslice(B, n1 + 1, n, p1 + 1, p);
     336       42267 :   M6 = ZM_mul_i(S4, B22, m1 + 1, n2 + 1, p2 + 1);
     337       42267 :   if (gc_needed(av, 1))
     338           6 :     gerepileall(av, 6, &M2, &M3, &T4, &C11, &V1, &M6);  /* destroy S4, B22 */
     339       42267 :   A22 = matslice(A, m1 + 1, m, n1 + 1, n);
     340       42267 :   M7 = ZM_mul_i(A22, T4, m2 + 1, n2 + 1, p1 + 1);
     341       42267 :   if (gc_needed(av, 1))
     342           0 :     gerepileall(av, 6, &M2, &M3, &C11, &V1, &M6, &M7);  /* destroy A22, T4 */
     343       42267 :   V3 = add_slices(m1, p2, V1, 0, m1, 0, p2, M3, 0, m2, 0, p2);
     344       42267 :   C12 = add_slices(m1, p2, V3, 0, m1, 0, p2, M6, 0, m1, 0, p2);
     345       42267 :   if (gc_needed(av, 1))
     346          12 :     gerepileall(av, 6, &M2, &M3, &C11, &V1, &M7, &C12);  /* destroy V3, M6 */
     347       42267 :   V2 = add_slices(m2, p1, V1, 0, m2, 0, p1, M2, 0, m2, 0, p2);
     348       42267 :   if (gc_needed(av, 1))
     349           0 :     gerepileall(av, 5, &M3, &C11, &M7, &C12, &V2);  /* destroy V1, M2 */
     350       42267 :   C21 = add_slices(m2, p1, V2, 0, m2, 0, p1, M7, 0, m2, 0, p1);
     351       42267 :   if (gc_needed(av, 1))
     352           6 :     gerepileall(av, 5, &M3, &C11, &C12, &V2, &C21);  /* destroy M7 */
     353       42267 :   C22 = add_slices(m2, p2, V2, 0, m2, 0, p2, M3, 0, m2, 0, p2);
     354       42267 :   if (gc_needed(av, 1))
     355           0 :     gerepileall(av, 4, &C11, &C12, &C21, &C22);  /* destroy V2, M3 */
     356       42267 :   C = shallowconcat(vconcat(C11, C21), vconcat(C12, C22));
     357       42267 :   return gerepilecopy(av, C);
     358             : }
     359             : 
     360             : /* x[i,]*y. Assume lg(x) > 1 and 0 < i < lgcols(x) */
     361             : static GEN
     362   527801312 : ZMrow_ZC_mul_i(GEN x, GEN y, long i, long lx)
     363             : {
     364   527801312 :   pari_sp av = avma;
     365   527801312 :   GEN c = mulii(gcoeff(x,i,1), gel(y,1)), ZERO = gen_0;
     366             :   long k;
     367  5232608615 :   for (k = 2; k < lx; k++)
     368             :   {
     369  4705191759 :     GEN t = mulii(gcoeff(x,i,k), gel(y,k));
     370  4703685500 :     if (t != ZERO) c = addii(c, t);
     371             :   }
     372   527416856 :   return gerepileuptoint(av, c);
     373             : }
     374             : GEN
     375   134535947 : ZMrow_ZC_mul(GEN x, GEN y, long i)
     376   134535947 : { return ZMrow_ZC_mul_i(x, y, i, lg(x)); }
     377             : 
     378             : /* return x * y, 1 < lx = lg(x), l = lgcols(x) */
     379             : static GEN
     380    74811413 : ZM_ZC_mul_i(GEN x, GEN y, long lx, long l)
     381             : {
     382    74811413 :   GEN z = cgetg(l,t_COL);
     383             :   long i;
     384   468094636 :   for (i=1; i<l; i++) gel(z,i) = ZMrow_ZC_mul_i(x,y,i,lx);
     385    74769904 :   return z;
     386             : }
     387             : 
     388             : static GEN
     389    14284621 : ZM_mul_classical(GEN x, GEN y, long l, long lx, long ly)
     390             : {
     391             :   long j;
     392    14284621 :   GEN z = cgetg(ly, t_MAT);
     393    66111534 :   for (j = 1; j < ly; j++)
     394    51830346 :     gel(z, j) = ZM_ZC_mul_i(x, gel(y, j), lx, l);
     395    14281188 :   return z;
     396             : }
     397             : 
     398             : static GEN
     399        1307 : ZM_mul_slice(GEN A, GEN B, GEN P, GEN *mod)
     400             : {
     401        1307 :   pari_sp av = avma;
     402        1307 :   long i, n = lg(P)-1;
     403             :   GEN H, T;
     404        1307 :   if (n == 1)
     405             :   {
     406           0 :     ulong p = uel(P,1);
     407           0 :     GEN a = ZM_to_Flm(A, p);
     408           0 :     GEN b = ZM_to_Flm(B, p);
     409           0 :     GEN Hp = gerepileupto(av, Flm_to_ZM(Flm_mul(a, b, p)));
     410           0 :     *mod = utoi(p); return Hp;
     411             :   }
     412        1307 :   T = ZV_producttree(P);
     413        1308 :   A = ZM_nv_mod_tree(A, P, T);
     414        1308 :   B = ZM_nv_mod_tree(B, P, T);
     415        1308 :   H = cgetg(n+1, t_VEC);
     416       10008 :   for(i=1; i <= n; i++)
     417        8700 :     gel(H,i) = Flm_mul(gel(A,i),gel(B,i),P[i]);
     418        1308 :   H = nmV_chinese_center_tree_seq(H, P, T, ZV_chinesetree(P, T));
     419        1308 :   *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
     420             : }
     421             : 
     422             : GEN
     423        1307 : ZM_mul_worker(GEN P, GEN A, GEN B)
     424             : {
     425        1307 :   GEN V = cgetg(3, t_VEC);
     426        1307 :   gel(V,1) = ZM_mul_slice(A, B, P, &gel(V,2));
     427        1308 :   return V;
     428             : }
     429             : 
     430             : static GEN
     431         974 : ZM_mul_fast(GEN A, GEN B, long lA, long lB, long sA, long sB)
     432             : {
     433         974 :   pari_sp av = avma;
     434             :   forprime_t S;
     435             :   GEN H, worker;
     436             :   long h;
     437         974 :   if (sA == 2 || sB == 2) return zeromat(nbrows(A),lB-1);
     438         960 :   h = 1 + (sA + sB - 4) * BITS_IN_LONG + expu(lA-1);
     439         960 :   init_modular_big(&S);
     440         960 :   worker = snm_closure(is_entry("_ZM_mul_worker"), mkvec2(A,B));
     441         960 :   H = gen_crt("ZM_mul", worker, &S, NULL, h, 0, NULL,
     442             :               nmV_chinese_center, FpM_center);
     443         960 :   return gerepileupto(av, H);
     444             : }
     445             : 
     446             : /* s = min(log_BIL |x|, log_BIL |y|), use Strassen-Winograd when
     447             :  * min(dims) > B */
     448             : static long
     449    14326891 : sw_bound(long s)
     450    14326891 : { return s > 60 ? 2: s > 25 ? 4: s > 15 ? 8 : s > 8 ? 16 : 32; }
     451             : 
     452             : static GEN
     453    21124596 : ZM_mul_i(GEN x, GEN y, long l, long lx, long ly)
     454             : {
     455             :   long sx, sy, B;
     456    21124596 :   if (lx==3 && l==3 && ly==3) return ZM2_mul(x,y);
     457    14300644 :   sx = ZM_max_lg_i(x, lx, l);
     458    14302050 :   sy = ZM_max_lg_i(y, ly, lx);
     459    14302120 :   if ((lx > 70 && ly > 70 && l > 70) && (sx <= 10 * sy && sy <= 10 * sx))
     460         974 :     return ZM_mul_fast(x, y, lx, ly, sx, sy);
     461             : 
     462    14301146 :   B = sw_bound(minss(sx, sy));
     463    14301139 :   if (l <= B || lx <= B || ly <= B)
     464    14258902 :     return ZM_mul_classical(x, y, l, lx, ly);
     465             :   else
     466       42237 :     return ZM_mul_sw(x, y, l - 1, lx - 1, ly - 1);
     467             : }
     468             : 
     469             : GEN
     470    20965700 : ZM_mul(GEN x, GEN y)
     471             : {
     472    20965700 :   long lx = lg(x), ly = lg(y);
     473    20965700 :   if (ly == 1) return cgetg(1,t_MAT);
     474    20830908 :   if (lx == 1) return zeromat(0, ly-1);
     475    20829326 :   return ZM_mul_i(x, y, lgcols(x), lx, ly);
     476             : }
     477             : 
     478             : GEN
     479      548927 : QM_mul(GEN x, GEN y)
     480             : {
     481      548927 :   GEN dx, nx = Q_primitive_part(x, &dx);
     482      548927 :   GEN dy, ny = Q_primitive_part(y, &dy);
     483      548927 :   GEN z = ZM_mul(nx, ny);
     484      548927 :   if (dx || dy)
     485             :   {
     486      465316 :     GEN d = dx ? dy ? gmul(dx, dy): dx : dy;
     487      465316 :     if (!gequal1(d)) z = ZM_Q_mul(z, d);
     488             :   }
     489      548927 :   return z;
     490             : }
     491             : 
     492             : GEN
     493         700 : QM_sqr(GEN x)
     494             : {
     495         700 :   GEN dx, nx = Q_primitive_part(x, &dx);
     496         700 :   GEN z = ZM_sqr(nx);
     497         700 :   if (dx)
     498         700 :     z = ZM_Q_mul(z, gsqr(dx));
     499         700 :   return z;
     500             : }
     501             : 
     502             : GEN
     503      142082 : QM_QC_mul(GEN x, GEN y)
     504             : {
     505      142082 :   GEN dx, nx = Q_primitive_part(x, &dx);
     506      142082 :   GEN dy, ny = Q_primitive_part(y, &dy);
     507      142082 :   GEN z = ZM_ZC_mul(nx, ny);
     508      142082 :   if (dx || dy)
     509             :   {
     510      142054 :     GEN d = dx ? dy ? gmul(dx, dy): dx : dy;
     511      142054 :     if (!gequal1(d)) z = ZC_Q_mul(z, d);
     512             :   }
     513      142082 :   return z;
     514             : }
     515             : 
     516             : /* assume result is symmetric */
     517             : GEN
     518           0 : ZM_multosym(GEN x, GEN y)
     519             : {
     520           0 :   long j, lx, ly = lg(y);
     521             :   GEN M;
     522           0 :   if (ly == 1) return cgetg(1,t_MAT);
     523           0 :   lx = lg(x); /* = lgcols(y) */
     524           0 :   if (lx == 1) return cgetg(1,t_MAT);
     525             :   /* ly = lgcols(x) */
     526           0 :   M = cgetg(ly, t_MAT);
     527           0 :   for (j=1; j<ly; j++)
     528             :   {
     529           0 :     GEN z = cgetg(ly,t_COL), yj = gel(y,j);
     530             :     long i;
     531           0 :     for (i=1; i<j; i++) gel(z,i) = gcoeff(M,j,i);
     532           0 :     for (i=j; i<ly; i++)gel(z,i) = ZMrow_ZC_mul_i(x,yj,i,lx);
     533           0 :     gel(M,j) = z;
     534             :   }
     535           0 :   return M;
     536             : }
     537             : 
     538             : /* compute m*diagonal(d), assume lg(d) = lg(m). Shallow */
     539             : GEN
     540          21 : ZM_mul_diag(GEN m, GEN d)
     541             : {
     542             :   long j, l;
     543          21 :   GEN y = cgetg_copy(m, &l);
     544          77 :   for (j=1; j<l; j++)
     545             :   {
     546          56 :     GEN c = gel(d,j);
     547          56 :     gel(y,j) = equali1(c)? gel(m,j): ZC_Z_mul(gel(m,j), c);
     548             :   }
     549          21 :   return y;
     550             : }
     551             : /* compute diagonal(d)*m, assume lg(d) = lg(m~). Shallow */
     552             : GEN
     553      533463 : ZM_diag_mul(GEN d, GEN m)
     554             : {
     555      533463 :   long i, j, l = lg(d), lm = lg(m);
     556      533463 :   GEN y = cgetg(lm, t_MAT);
     557     1520044 :   for (j=1; j<lm; j++) gel(y,j) = cgetg(l, t_COL);
     558     1646454 :   for (i=1; i<l; i++)
     559             :   {
     560     1113330 :     GEN c = gel(d,i);
     561     1113330 :     if (equali1(c))
     562      236732 :       for (j=1; j<lm; j++) gcoeff(y,i,j) = gcoeff(m,i,j);
     563             :     else
     564     3396533 :       for (j=1; j<lm; j++) gcoeff(y,i,j) = mulii(gcoeff(m,i,j), c);
     565             :   }
     566      533124 :   return y;
     567             : }
     568             : 
     569             : /* assume lx > 1 is lg(x) = lg(y) */
     570             : static GEN
     571    19056310 : ZV_dotproduct_i(GEN x, GEN y, long lx)
     572             : {
     573    19056310 :   pari_sp av = avma;
     574    19056310 :   GEN c = mulii(gel(x,1), gel(y,1));
     575             :   long i;
     576   143195163 :   for (i = 2; i < lx; i++)
     577             :   {
     578   124140304 :     GEN t = mulii(gel(x,i), gel(y,i));
     579   124139750 :     if (t != gen_0) c = addii(c, t);
     580             :   }
     581    19054859 :   return gerepileuptoint(av, c);
     582             : }
     583             : 
     584             : /* x~ * y, assuming result is symmetric */
     585             : GEN
     586      524928 : ZM_transmultosym(GEN x, GEN y)
     587             : {
     588      524928 :   long i, j, l, ly = lg(y);
     589             :   GEN M;
     590      524928 :   if (ly == 1) return cgetg(1,t_MAT);
     591             :   /* lg(x) = ly */
     592      524928 :   l = lgcols(y); /* = lgcols(x) */
     593      524928 :   M = cgetg(ly, t_MAT);
     594     2672773 :   for (i=1; i<ly; i++)
     595             :   {
     596     2147845 :     GEN xi = gel(x,i), c = cgetg(ly,t_COL);
     597     2147845 :     gel(M,i) = c;
     598     7010463 :     for (j=1; j<i; j++)
     599     4862618 :       gcoeff(M,i,j) = gel(c,j) = ZV_dotproduct_i(xi,gel(y,j),l);
     600     2147845 :     gel(c,i) = ZV_dotproduct_i(xi,gel(y,i),l);
     601             :   }
     602      524928 :   return M;
     603             : }
     604             : /* x~ * y */
     605             : GEN
     606        2289 : ZM_transmul(GEN x, GEN y)
     607             : {
     608        2289 :   long i, j, l, lx, ly = lg(y);
     609             :   GEN M;
     610        2289 :   if (ly == 1) return cgetg(1,t_MAT);
     611        2289 :   lx = lg(x);
     612        2289 :   l = lgcols(y);
     613        2289 :   if (lgcols(x) != l) pari_err_OP("operation 'ZM_transmul'", x,y);
     614        2289 :   M = cgetg(ly, t_MAT);
     615        6993 :   for (i=1; i<ly; i++)
     616             :   {
     617        4704 :     GEN yi = gel(y,i), c = cgetg(lx,t_COL);
     618        4704 :     gel(M,i) = c;
     619       12229 :     for (j=1; j<lx; j++) gel(c,j) = ZV_dotproduct_i(yi,gel(x,j),l);
     620             :   }
     621        2289 :   return M;
     622             : }
     623             : 
     624             : static GEN
     625       25753 : ZM_sqr_i(GEN x, long l)
     626             : {
     627       25753 :   long s = ZM_max_lg_i(x, l, l);
     628       25753 :   if (l > 70) return ZM_mul_fast(x, x, l, l, s, s);
     629       25753 :   if (l <= sw_bound(s))
     630       25723 :     return ZM_mul_classical(x, x, l, l, l);
     631             :   else
     632          30 :     return ZM_mul_sw(x, x, l - 1, l - 1, l - 1);
     633             : }
     634             : 
     635             : GEN
     636       25753 : ZM_sqr(GEN x)
     637             : {
     638       25753 :   long lx=lg(x);
     639       25753 :   if (lx==1) return cgetg(1,t_MAT);
     640       25753 :   return ZM_sqr_i(x, lx);
     641             : }
     642             : GEN
     643    23070151 : ZM_ZC_mul(GEN x, GEN y)
     644             : {
     645    23070151 :   long lx = lg(x);
     646    23070151 :   return lx==1? cgetg(1,t_COL): ZM_ZC_mul_i(x, y, lx, lgcols(x));
     647             : }
     648             : 
     649             : GEN
     650     4493171 : ZC_Z_div(GEN x, GEN c)
     651    19443748 : { pari_APPLY_type(t_COL, Qdivii(gel(x,i), c)) }
     652             : 
     653             : GEN
     654       15782 : ZM_Z_div(GEN x, GEN c)
     655      177995 : { pari_APPLY_same(ZC_Z_div(gel(x, i), c)) }
     656             : 
     657             : GEN
     658     2697398 : ZC_Q_mul(GEN A, GEN z)
     659             : {
     660     2697398 :   pari_sp av = avma;
     661     2697398 :   long i, l = lg(A);
     662             :   GEN d, n, Ad, B, u;
     663     2697398 :   if (typ(z)==t_INT) return ZC_Z_mul(A,z);
     664     2691140 :   n = gel(z, 1); d = gel(z, 2);
     665     2691140 :   Ad = FpC_red(A, d);
     666     2691070 :   u = gcdii(d, FpV_factorback(Ad, NULL, d));
     667     2691084 :   B = cgetg(l, t_COL);
     668     2691071 :   if (equali1(u))
     669             :   {
     670      420633 :     for(i=1; i<l; i++)
     671      354570 :       gel(B, i) = mkfrac(mulii(n, gel(A,i)), d);
     672             :   } else
     673             :   {
     674    18180819 :     for(i=1; i<l; i++)
     675             :     {
     676    15555882 :       GEN di = gcdii(gel(Ad, i), u), ni = mulii(n, diviiexact(gel(A,i), di));
     677    15555827 :       if (equalii(d, di))
     678    10626513 :         gel(B, i) = ni;
     679             :       else
     680     4929275 :         gel(B, i) = mkfrac(ni, diviiexact(d, di));
     681             :     }
     682             :   }
     683     2691000 :   return gerepilecopy(av, B);
     684             : }
     685             : 
     686             : GEN
     687     1092269 : ZM_Q_mul(GEN x, GEN z)
     688             : {
     689     1092269 :   if (typ(z)==t_INT) return ZM_Z_mul(x,z);
     690     3215017 :   pari_APPLY_same(ZC_Q_mul(gel(x, i), z));
     691             : }
     692             : 
     693             : long
     694   196399119 : zv_dotproduct(GEN x, GEN y)
     695             : {
     696   196399119 :   long i, lx = lg(x);
     697             :   ulong c;
     698   196399119 :   if (lx == 1) return 0;
     699   196399119 :   c = uel(x,1)*uel(y,1);
     700  3047359154 :   for (i = 2; i < lx; i++)
     701  2850960035 :     c += uel(x,i)*uel(y,i);
     702   196399119 :   return c;
     703             : }
     704             : 
     705             : GEN
     706      230629 : ZV_ZM_mul(GEN x, GEN y)
     707             : {
     708      230629 :   long i, lx = lg(x), ly = lg(y);
     709             :   GEN z;
     710      230629 :   if (lx == 1) return zerovec(ly-1);
     711      230510 :   z = cgetg(ly, t_VEC);
     712      882063 :   for (i = 1; i < ly; i++) gel(z,i) = ZV_dotproduct_i(x, gel(y,i), lx);
     713      230510 :   return z;
     714             : }
     715             : 
     716             : GEN
     717           0 : ZC_ZV_mul(GEN x, GEN y)
     718             : {
     719           0 :   long i,j, lx=lg(x), ly=lg(y);
     720             :   GEN z;
     721           0 :   if (ly==1) return cgetg(1,t_MAT);
     722           0 :   z = cgetg(ly,t_MAT);
     723           0 :   for (j=1; j < ly; j++)
     724             :   {
     725           0 :     gel(z,j) = cgetg(lx,t_COL);
     726           0 :     for (i=1; i<lx; i++) gcoeff(z,i,j) = mulii(gel(x,i),gel(y,j));
     727             :   }
     728           0 :   return z;
     729             : }
     730             : 
     731             : GEN
     732     6636074 : ZV_dotsquare(GEN x)
     733             : {
     734             :   long i, lx;
     735             :   pari_sp av;
     736             :   GEN z;
     737     6636074 :   lx = lg(x);
     738     6636074 :   if (lx == 1) return gen_0;
     739     6636074 :   av = avma; z = sqri(gel(x,1));
     740    26295754 :   for (i=2; i<lx; i++) z = addii(z, sqri(gel(x,i)));
     741     6634390 :   return gerepileuptoint(av,z);
     742             : }
     743             : 
     744             : GEN
     745    16154868 : ZV_dotproduct(GEN x,GEN y)
     746             : {
     747             :   long lx;
     748    16154868 :   if (x == y) return ZV_dotsquare(x);
     749    11385937 :   lx = lg(x);
     750    11385937 :   if (lx == 1) return gen_0;
     751    11385937 :   return ZV_dotproduct_i(x, y, lx);
     752             : }
     753             : 
     754             : static GEN
     755         280 : _ZM_mul(void *data /*ignored*/, GEN x, GEN y)
     756         280 : { (void)data; return ZM_mul(x,y); }
     757             : static GEN
     758       23184 : _ZM_sqr(void *data /*ignored*/, GEN x)
     759       23184 : { (void)data; return ZM_sqr(x); }
     760             : GEN
     761           0 : ZM_pow(GEN x, GEN n)
     762             : {
     763           0 :   pari_sp av = avma;
     764           0 :   if (!signe(n)) return matid(lg(x)-1);
     765           0 :   return gerepilecopy(av, gen_pow_i(x, n, NULL, &_ZM_sqr, &_ZM_mul));
     766             : }
     767             : GEN
     768       22638 : ZM_powu(GEN x, ulong n)
     769             : {
     770       22638 :   pari_sp av = avma;
     771       22638 :   if (!n) return matid(lg(x)-1);
     772       22638 :   return gerepilecopy(av, gen_powu_i(x, n, NULL, &_ZM_sqr, &_ZM_mul));
     773             : }
     774             : /********************************************************************/
     775             : /**                                                                **/
     776             : /**                           ADD, SUB                             **/
     777             : /**                                                                **/
     778             : /********************************************************************/
     779             : static GEN
     780    34468643 : ZC_add_i(GEN x, GEN y, long lx)
     781             : {
     782    34468643 :   GEN A = cgetg(lx, t_COL);
     783             :   long i;
     784   442114811 :   for (i=1; i<lx; i++) gel(A,i) = addii(gel(x,i), gel(y,i));
     785    34459959 :   return A;
     786             : }
     787             : GEN
     788    26387280 : ZC_add(GEN x, GEN y) { return ZC_add_i(x, y, lg(x)); }
     789             : GEN
     790      363161 : ZC_Z_add(GEN x, GEN y)
     791             : {
     792      363161 :   long k, lx = lg(x);
     793      363161 :   GEN z = cgetg(lx, t_COL);
     794      363159 :   if (lx == 1) pari_err_TYPE2("+",x,y);
     795      363159 :   gel(z,1) = addii(y,gel(x,1));
     796     2500433 :   for (k = 2; k < lx; k++) gel(z,k) = icopy(gel(x,k));
     797      363146 :   return z;
     798             : }
     799             : 
     800             : static GEN
     801     8951324 : ZC_sub_i(GEN x, GEN y, long lx)
     802             : {
     803             :   long i;
     804     8951324 :   GEN A = cgetg(lx, t_COL);
     805    64017984 :   for (i=1; i<lx; i++) gel(A,i) = subii(gel(x,i), gel(y,i));
     806     8950696 :   return A;
     807             : }
     808             : GEN
     809     8722338 : ZC_sub(GEN x, GEN y) { return ZC_sub_i(x, y, lg(x)); }
     810             : GEN
     811           0 : ZC_Z_sub(GEN x, GEN y)
     812             : {
     813           0 :   long k, lx = lg(x);
     814           0 :   GEN z = cgetg(lx, t_COL);
     815           0 :   if (lx == 1) pari_err_TYPE2("+",x,y);
     816           0 :   gel(z,1) = subii(gel(x,1), y);
     817           0 :   for (k = 2; k < lx; k++) gel(z,k) = icopy(gel(x,k));
     818           0 :   return z;
     819             : }
     820             : GEN
     821      544224 : Z_ZC_sub(GEN a, GEN x)
     822             : {
     823      544224 :   long k, lx = lg(x);
     824      544224 :   GEN z = cgetg(lx, t_COL);
     825      544230 :   if (lx == 1) pari_err_TYPE2("-",a,x);
     826      544230 :   gel(z,1) = subii(a, gel(x,1));
     827     1508396 :   for (k = 2; k < lx; k++) gel(z,k) = negi(gel(x,k));
     828      544229 :   return z;
     829             : }
     830             : 
     831             : GEN
     832      747536 : ZM_add(GEN x, GEN y)
     833             : {
     834      747536 :   long lx = lg(x), l, j;
     835             :   GEN z;
     836      747536 :   if (lx == 1) return cgetg(1, t_MAT);
     837      668198 :   z = cgetg(lx, t_MAT); l = lgcols(x);
     838     8749501 :   for (j = 1; j < lx; j++) gel(z,j) = ZC_add_i(gel(x,j), gel(y,j), l);
     839      668196 :   return z;
     840             : }
     841             : GEN
     842       65162 : ZM_sub(GEN x, GEN y)
     843             : {
     844       65162 :   long lx = lg(x), l, j;
     845             :   GEN z;
     846       65162 :   if (lx == 1) return cgetg(1, t_MAT);
     847       65162 :   z = cgetg(lx, t_MAT); l = lgcols(x);
     848      294047 :   for (j = 1; j < lx; j++) gel(z,j) = ZC_sub_i(gel(x,j), gel(y,j), l);
     849       65163 :   return z;
     850             : }
     851             : /********************************************************************/
     852             : /**                                                                **/
     853             : /**                         LINEAR COMBINATION                     **/
     854             : /**                                                                **/
     855             : /********************************************************************/
     856             : /* return X/c assuming division is exact */
     857             : GEN
     858     4431850 : ZC_Z_divexact(GEN x, GEN c)
     859    45852827 : { pari_APPLY_type(t_COL, diviiexact(gel(x,i), c)) }
     860             : GEN
     861        2261 : ZC_divexactu(GEN x, ulong c)
     862       11375 : { pari_APPLY_type(t_COL, diviuexact(gel(x,i), c)) }
     863             : 
     864             : GEN
     865      429794 : ZM_Z_divexact(GEN x, GEN c)
     866     2798483 : { pari_APPLY_same(ZC_Z_divexact(gel(x,i), c)) }
     867             : 
     868             : GEN
     869         441 : ZM_divexactu(GEN x, ulong c)
     870        2702 : { pari_APPLY_same(ZC_divexactu(gel(x,i), c)) }
     871             : 
     872             : GEN
     873    35106444 : ZC_Z_mul(GEN x, GEN c)
     874             : {
     875    35106444 :   if (!signe(c)) return zerocol(lg(x)-1);
     876    33750220 :   if (is_pm1(c)) return (signe(c) > 0)? ZC_copy(x): ZC_neg(x);
     877   254331359 :   pari_APPLY_type(t_COL, mulii(gel(x,i), c))
     878             : }
     879             : 
     880             : GEN
     881       61732 : ZC_z_mul(GEN x, long c)
     882             : {
     883       61732 :   if (!c) return zerocol(lg(x)-1);
     884       53245 :   if (c == 1) return ZC_copy(x);
     885       48391 :   if (c ==-1) return ZC_neg(x);
     886      479015 :   pari_APPLY_type(t_COL, mulsi(c, gel(x,i)))
     887             : }
     888             : 
     889             : GEN
     890       28365 : zv_z_mul(GEN x, long n)
     891      430611 : { pari_APPLY_long(x[i]*n) }
     892             : 
     893             : /* return a ZM */
     894             : GEN
     895           0 : nm_Z_mul(GEN X, GEN c)
     896             : {
     897           0 :   long i, j, h, l = lg(X), s = signe(c);
     898             :   GEN A;
     899           0 :   if (l == 1) return cgetg(1, t_MAT);
     900           0 :   h = lgcols(X);
     901           0 :   if (!s) return zeromat(h-1, l-1);
     902           0 :   if (is_pm1(c)) {
     903           0 :     if (s > 0) return Flm_to_ZM(X);
     904           0 :     X = Flm_to_ZM(X); ZM_togglesign(X); return X;
     905             :   }
     906           0 :   A = cgetg(l, t_MAT);
     907           0 :   for (j = 1; j < l; j++)
     908             :   {
     909           0 :     GEN a = cgetg(h, t_COL), x = gel(X, j);
     910           0 :     for (i = 1; i < h; i++) gel(a,i) = muliu(c, x[i]);
     911           0 :     gel(A,j) = a;
     912             :   }
     913           0 :   return A;
     914             : }
     915             : GEN
     916     2984553 : ZM_Z_mul(GEN X, GEN c)
     917             : {
     918     2984553 :   long i, j, h, l = lg(X);
     919             :   GEN A;
     920     2984553 :   if (l == 1) return cgetg(1, t_MAT);
     921     2984553 :   h = lgcols(X);
     922     2984541 :   if (!signe(c)) return zeromat(h-1, l-1);
     923     2939274 :   if (is_pm1(c)) return (signe(c) > 0)? ZM_copy(X): ZM_neg(X);
     924     2653906 :   A = cgetg(l, t_MAT);
     925    13338277 :   for (j = 1; j < l; j++)
     926             :   {
     927    10684551 :     GEN a = cgetg(h, t_COL), x = gel(X, j);
     928   168378311 :     for (i = 1; i < h; i++) gel(a,i) = mulii(c, gel(x,i));
     929    10684377 :     gel(A,j) = a;
     930             :   }
     931     2653726 :   return A;
     932             : }
     933             : void
     934    71937713 : ZC_lincomb1_inplace_i(GEN X, GEN Y, GEN v, long n)
     935             : {
     936             :   long i;
     937  1564200412 :   for (i = n; i; i--) gel(X,i) = addmulii_inplace(gel(X,i), gel(Y,i), v);
     938    71892200 : }
     939             : /* X <- X + v Y (elementary col operation) */
     940             : void
     941    61204309 : ZC_lincomb1_inplace(GEN X, GEN Y, GEN v)
     942             : {
     943    61204309 :   if (lgefint(v) != 2) return ZC_lincomb1_inplace_i(X, Y, v, lg(X)-1);
     944             : }
     945             : void
     946    29465613 : Flc_lincomb1_inplace(GEN X, GEN Y, ulong v, ulong q)
     947             : {
     948             :   long i;
     949    29465613 :   if (!v) return; /* v = 0 */
     950   742209942 :   for (i = lg(X)-1; i; i--) X[i] = Fl_add(X[i], Fl_mul(Y[i], v, q), q);
     951             : }
     952             : 
     953             : /* X + v Y, wasteful if (v = 0) */
     954             : static GEN
     955    15324472 : ZC_lincomb1(GEN v, GEN x, GEN y)
     956   125190991 : { pari_APPLY_type(t_COL, addmulii(gel(x,i), gel(y,i), v)) }
     957             : 
     958             : /* -X + vY */
     959             : static GEN
     960      729417 : ZC_lincomb_1(GEN v, GEN x, GEN y)
     961     4529915 : { pari_APPLY_type(t_COL, mulsubii(gel(y,i), v, gel(x,i))) }
     962             : 
     963             : /* X,Y compatible ZV; u,v in Z. Returns A = u*X + v*Y */
     964             : GEN
     965    32857439 : ZC_lincomb(GEN u, GEN v, GEN X, GEN Y)
     966             : {
     967             :   long su, sv;
     968             :   GEN A;
     969             : 
     970    32857439 :   su = signe(u); if (!su) return ZC_Z_mul(Y, v);
     971    32856480 :   sv = signe(v); if (!sv) return ZC_Z_mul(X, u);
     972    32427297 :   if (is_pm1(v))
     973             :   {
     974    11023186 :     if (is_pm1(u))
     975             :     {
     976     9717665 :       if (su != sv) A = ZC_sub(X, Y);
     977     2700701 :       else          A = ZC_add(X, Y);
     978     9717118 :       if (su < 0) ZV_togglesign(A); /* in place but was created above */
     979             :     }
     980             :     else
     981             :     {
     982     1305436 :       if (sv > 0) A = ZC_lincomb1 (u, Y, X);
     983      600192 :       else        A = ZC_lincomb_1(u, Y, X);
     984             :     }
     985             :   }
     986    21404874 :   else if (is_pm1(u))
     987             :   {
     988    14747724 :     if (su > 0) A = ZC_lincomb1 (v, X, Y);
     989      128470 :     else        A = ZC_lincomb_1(v, X, Y);
     990             :   }
     991             :   else
     992             :   { /* not cgetg_copy: x may be a t_VEC */
     993     6661940 :     long i, lx = lg(X);
     994     6661940 :     A = cgetg(lx,t_COL);
     995    43587769 :     for (i=1; i<lx; i++) gel(A,i) = lincombii(u,v,gel(X,i),gel(Y,i));
     996             :   }
     997    32419971 :   return A;
     998             : }
     999             : 
    1000             : /********************************************************************/
    1001             : /**                                                                **/
    1002             : /**                           CONVERSIONS                          **/
    1003             : /**                                                                **/
    1004             : /********************************************************************/
    1005             : GEN
    1006      400844 : ZV_to_nv(GEN x)
    1007      749154 : { pari_APPLY_ulong(itou(gel(x,i))) }
    1008             : 
    1009             : GEN
    1010      130058 : zm_to_ZM(GEN x)
    1011      815171 : { pari_APPLY_type(t_MAT, zc_to_ZC(gel(x,i))) }
    1012             : 
    1013             : GEN
    1014         126 : zmV_to_ZMV(GEN x)
    1015         791 : { pari_APPLY_type(t_VEC, zm_to_ZM(gel(x,i))) }
    1016             : 
    1017             : /* same as Flm_to_ZM but do not assume positivity */
    1018             : GEN
    1019        1022 : ZM_to_zm(GEN x)
    1020       17472 : { pari_APPLY_same(ZV_to_zv(gel(x,i))) }
    1021             : 
    1022             : GEN
    1023      366646 : zv_to_Flv(GEN x, ulong p)
    1024     5418812 : { pari_APPLY_ulong(umodsu(x[i], p)) }
    1025             : 
    1026             : GEN
    1027       22694 : zm_to_Flm(GEN x, ulong p)
    1028      351750 : { pari_APPLY_same(zv_to_Flv(gel(x,i),p)) }
    1029             : 
    1030             : GEN
    1031          49 : ZMV_to_zmV(GEN x)
    1032         399 : { pari_APPLY_type(t_VEC, ZM_to_zm(gel(x,i))) }
    1033             : 
    1034             : /********************************************************************/
    1035             : /**                                                                **/
    1036             : /**                         COPY, NEGATION                         **/
    1037             : /**                                                                **/
    1038             : /********************************************************************/
    1039             : GEN
    1040    13268034 : ZC_copy(GEN x)
    1041             : {
    1042    13268034 :   long i, lx = lg(x);
    1043    13268034 :   GEN y = cgetg(lx, t_COL);
    1044    83864208 :   for (i=1; i<lx; i++)
    1045             :   {
    1046    70595187 :     GEN c = gel(x,i);
    1047    70595187 :     gel(y,i) = lgefint(c) == 2? gen_0: icopy(c);
    1048             :   }
    1049    13269021 :   return y;
    1050             : }
    1051             : 
    1052             : GEN
    1053      509468 : ZM_copy(GEN x)
    1054     4284495 : { pari_APPLY_same(ZC_copy(gel(x,i))) }
    1055             : 
    1056             : void
    1057      342088 : ZV_neg_inplace(GEN M)
    1058             : {
    1059      342088 :   long l = lg(M);
    1060     1287695 :   while (--l > 0) gel(M,l) = negi(gel(M,l));
    1061      342173 : }
    1062             : GEN
    1063     6264151 : ZC_neg(GEN x)
    1064    36496275 : { pari_APPLY_type(t_COL, negi(gel(x,i))) }
    1065             : 
    1066             : GEN
    1067       51082 : zv_neg(GEN x)
    1068      655368 : { pari_APPLY_long(-x[i]) }
    1069             : GEN
    1070         126 : zv_neg_inplace(GEN M)
    1071             : {
    1072         126 :   long l = lg(M);
    1073         432 :   while (--l > 0) M[l] = -M[l];
    1074         126 :   return M;
    1075             : }
    1076             : GEN
    1077          77 : zv_abs(GEN x)
    1078        5446 : { pari_APPLY_ulong(labs(x[i])) }
    1079             : GEN
    1080     1654943 : ZM_neg(GEN x)
    1081     5038814 : { pari_APPLY_same(ZC_neg(gel(x,i))) }
    1082             : 
    1083             : void
    1084     4920051 : ZV_togglesign(GEN M)
    1085             : {
    1086     4920051 :   long l = lg(M);
    1087    73290897 :   while (--l > 0) togglesign_safe(&gel(M,l));
    1088     4920093 : }
    1089             : void
    1090           0 : ZM_togglesign(GEN M)
    1091             : {
    1092           0 :   long l = lg(M);
    1093           0 :   while (--l > 0) ZV_togglesign(gel(M,l));
    1094           0 : }
    1095             : 
    1096             : /********************************************************************/
    1097             : /**                                                                **/
    1098             : /**                        "DIVISION" mod HNF                      **/
    1099             : /**                                                                **/
    1100             : /********************************************************************/
    1101             : /* Reduce ZC x modulo ZM y in HNF, may return x itself (not a copy) */
    1102             : GEN
    1103    10690465 : ZC_hnfremdiv(GEN x, GEN y, GEN *Q)
    1104             : {
    1105    10690465 :   long i, l = lg(x);
    1106             :   GEN q;
    1107             : 
    1108    10690465 :   if (Q) *Q = cgetg(l,t_COL);
    1109    10690467 :   if (l == 1) return cgetg(1,t_COL);
    1110    61543249 :   for (i = l-1; i>0; i--)
    1111             :   {
    1112    50856088 :     q = diviiround(gel(x,i), gcoeff(y,i,i));
    1113    50855000 :     if (signe(q)) {
    1114    22935652 :       togglesign(q);
    1115    22935714 :       x = ZC_lincomb(gen_1, q, x, gel(y,i));
    1116             :     }
    1117    50852782 :     if (Q) gel(*Q, i) = q;
    1118             :   }
    1119    10687161 :   return x;
    1120             : }
    1121             : 
    1122             : /* x = y Q + R, may return some columns of x (not copies) */
    1123             : GEN
    1124      453843 : ZM_hnfdivrem(GEN x, GEN y, GEN *Q)
    1125             : {
    1126      453843 :   long lx = lg(x), i;
    1127      453843 :   GEN R = cgetg(lx, t_MAT);
    1128      453872 :   if (Q)
    1129             :   {
    1130      127355 :     GEN q = cgetg(lx, t_MAT); *Q = q;
    1131      188189 :     for (i=1; i<lx; i++) gel(R,i) = ZC_hnfremdiv(gel(x,i),y,(GEN*)(q+i));
    1132             :   }
    1133             :   else
    1134      823179 :     for (i=1; i<lx; i++)
    1135             :     {
    1136      496662 :       pari_sp av = avma;
    1137      496662 :       GEN z = ZC_hnfrem(gel(x,i),y);
    1138      496613 :       gel(R,i) = (avma == av)? ZC_copy(z): gerepileupto(av, z);
    1139             :     }
    1140      453872 :   return R;
    1141             : }
    1142             : 
    1143             : /********************************************************************/
    1144             : /**                                                                **/
    1145             : /**                               TESTS                            **/
    1146             : /**                                                                **/
    1147             : /********************************************************************/
    1148             : int
    1149    23103510 : zv_equal0(GEN V)
    1150             : {
    1151    23103510 :   long l = lg(V);
    1152    37579930 :   while (--l > 0)
    1153    30805766 :     if (V[l]) return 0;
    1154     6774164 :   return 1;
    1155             : }
    1156             : 
    1157             : int
    1158    14366450 : ZV_equal0(GEN V)
    1159             : {
    1160    14366450 :   long l = lg(V);
    1161    25396410 :   while (--l > 0)
    1162    24970259 :     if (signe(gel(V,l))) return 0;
    1163      426151 :   return 1;
    1164             : }
    1165             : int
    1166       16409 : ZMrow_equal0(GEN V, long i)
    1167             : {
    1168       16409 :   long l = lg(V);
    1169       25507 :   while (--l > 0)
    1170       21857 :     if (signe(gcoeff(V,i,l))) return 0;
    1171        3650 :   return 1;
    1172             : }
    1173             : 
    1174             : static int
    1175     6257921 : ZV_equal_lg(GEN V, GEN W, long l)
    1176             : {
    1177    25926149 :   while (--l > 0)
    1178    20150130 :     if (!equalii(gel(V,l), gel(W,l))) return 0;
    1179     5776019 :   return 1;
    1180             : }
    1181             : int
    1182      296309 : ZV_equal(GEN V, GEN W)
    1183             : {
    1184      296309 :   long l = lg(V);
    1185      296309 :   if (lg(W) != l) return 0;
    1186      296302 :   return ZV_equal_lg(V, W, l);
    1187             : }
    1188             : int
    1189     3347001 : ZM_equal(GEN A, GEN B)
    1190             : {
    1191     3347001 :   long i, m, l = lg(A);
    1192     3347001 :   if (lg(B) != l) return 0;
    1193     3346947 :   if (l == 1) return 1;
    1194     3346947 :   m = lgcols(A);
    1195     3346954 :   if (lgcols(B) != m) return 0;
    1196     9047510 :   for (i = 1; i < l; i++)
    1197     5961590 :     if (!ZV_equal_lg(gel(A,i), gel(B,i), m)) return 0;
    1198     3085920 :   return 1;
    1199             : }
    1200             : int
    1201       72400 : ZM_equal0(GEN A)
    1202             : {
    1203       72400 :   long i, j, m, l = lg(A);
    1204       72400 :   if (l == 1) return 1;
    1205       72400 :   m = lgcols(A);
    1206      123463 :   for (j = 1; j < l; j++)
    1207     2715086 :     for (i = 1; i < m; i++)
    1208     2664023 :       if (signe(gcoeff(A,i,j))) return 0;
    1209       15585 :   return 1;
    1210             : }
    1211             : int
    1212    30831559 : zv_equal(GEN V, GEN W)
    1213             : {
    1214    30831559 :   long l = lg(V);
    1215    30831559 :   if (lg(W) != l) return 0;
    1216   262051519 :   while (--l > 0)
    1217   232337181 :     if (V[l] != W[l]) return 0;
    1218    29714338 :   return 1;
    1219             : }
    1220             : 
    1221             : int
    1222        1638 : zvV_equal(GEN V, GEN W)
    1223             : {
    1224        1638 :   long l = lg(V);
    1225        1638 :   if (lg(W) != l) return 0;
    1226       80388 :   while (--l > 0)
    1227       79912 :     if (!zv_equal(gel(V,l),gel(W,l))) return 0;
    1228         476 :   return 1;
    1229             : }
    1230             : 
    1231             : int
    1232      759277 : ZM_ishnf(GEN x)
    1233             : {
    1234      759277 :   long i,j, lx = lg(x);
    1235     2282309 :   for (i=1; i<lx; i++)
    1236             :   {
    1237     1635769 :     GEN xii = gcoeff(x,i,i);
    1238     1635769 :     if (signe(xii) <= 0) return 0;
    1239     3176077 :     for (j=1; j<i; j++)
    1240     1577709 :       if (signe(gcoeff(x,i,j))) return 0;
    1241     3234584 :     for (j=i+1; j<lx; j++)
    1242             :     {
    1243     1711552 :       GEN xij = gcoeff(x,i,j);
    1244     1711552 :       if (signe(xij)<0 || cmpii(xij,xii)>=0) return 0;
    1245             :     }
    1246             :   }
    1247      646540 :   return 1;
    1248             : }
    1249             : int
    1250      638019 : ZM_isidentity(GEN x)
    1251             : {
    1252      638019 :   long i,j, lx = lg(x);
    1253             : 
    1254      638019 :   if (lx == 1) return 1;
    1255      638012 :   if (lx != lgcols(x)) return 0;
    1256     3116872 :   for (j=1; j<lx; j++)
    1257             :   {
    1258     2480918 :     GEN c = gel(x,j);
    1259     7694533 :     for (i=1; i<j; )
    1260     5213651 :       if (signe(gel(c,i++))) return 0;
    1261             :     /* i = j */
    1262     2480882 :     if (!equali1(gel(c,i++))) return 0;
    1263     7700603 :     for (   ; i<lx; )
    1264     5221729 :       if (signe(gel(c,i++))) return 0;
    1265             :   }
    1266      635954 :   return 1;
    1267             : }
    1268             : int
    1269      556002 : ZM_isdiagonal(GEN x)
    1270             : {
    1271      556002 :   long i,j, lx = lg(x);
    1272      556002 :   if (lx == 1) return 1;
    1273      556002 :   if (lx != lgcols(x)) return 0;
    1274             : 
    1275     1437320 :   for (j=1; j<lx; j++)
    1276             :   {
    1277     1207865 :     GEN c = gel(x,j);
    1278     1695654 :     for (i=1; i<j; i++)
    1279      814296 :       if (signe(gel(c,i))) return 0;
    1280     2018893 :     for (i++; i<lx; i++)
    1281     1137570 :       if (signe(gel(c,i))) return 0;
    1282             :   }
    1283      229455 :   return 1;
    1284             : }
    1285             : int
    1286      103945 : ZM_isscalar(GEN x, GEN s)
    1287             : {
    1288      103945 :   long i, j, lx = lg(x);
    1289             : 
    1290      103945 :   if (lx == 1) return 1;
    1291      103945 :   if (!s) s = gcoeff(x,1,1);
    1292      103945 :   if (equali1(s)) return ZM_isidentity(x);
    1293      102734 :   if (lx != lgcols(x)) return 0;
    1294      547634 :   for (j=1; j<lx; j++)
    1295             :   {
    1296      456865 :     GEN c = gel(x,j);
    1297     2251233 :     for (i=1; i<j; )
    1298     1804343 :       if (signe(gel(c,i++))) return 0;
    1299             :     /* i = j */
    1300      446890 :     if (!equalii(gel(c,i++), s)) return 0;
    1301     2278760 :     for (   ; i<lx; )
    1302     1833860 :       if (signe(gel(c,i++))) return 0;
    1303             :   }
    1304       90769 :   return 1;
    1305             : }
    1306             : 
    1307             : long
    1308      106942 : ZC_is_ei(GEN x)
    1309             : {
    1310      106942 :   long i, j = 0, l = lg(x);
    1311      990877 :   for (i = 1; i < l; i++)
    1312             :   {
    1313      883936 :     GEN c = gel(x,i);
    1314      883936 :     long s = signe(c);
    1315      883936 :     if (!s) continue;
    1316      106936 :     if (s < 0 || !is_pm1(c) || j) return 0;
    1317      106935 :     j = i;
    1318             :   }
    1319      106941 :   return j;
    1320             : }
    1321             : 
    1322             : /********************************************************************/
    1323             : /**                                                                **/
    1324             : /**                       MISCELLANEOUS                            **/
    1325             : /**                                                                **/
    1326             : /********************************************************************/
    1327             : /* assume lg(x) = lg(y), x,y in Z^n */
    1328             : int
    1329     3133006 : ZV_cmp(GEN x, GEN y)
    1330             : {
    1331     3133006 :   long fl,i, lx = lg(x);
    1332     6344641 :   for (i=1; i<lx; i++)
    1333     5057521 :     if (( fl = cmpii(gel(x,i), gel(y,i)) )) return fl;
    1334     1287120 :   return 0;
    1335             : }
    1336             : /* assume lg(x) = lg(y), x,y in Z^n */
    1337             : int
    1338       19740 : ZV_abscmp(GEN x, GEN y)
    1339             : {
    1340       19740 :   long fl,i, lx = lg(x);
    1341       54049 :   for (i=1; i<lx; i++)
    1342       53922 :     if (( fl = abscmpii(gel(x,i), gel(y,i)) )) return fl;
    1343         127 :   return 0;
    1344             : }
    1345             : 
    1346             : long
    1347    20719073 : zv_content(GEN x)
    1348             : {
    1349    20719073 :   long i, s, l = lg(x);
    1350    20719073 :   if (l == 1) return 0;
    1351    20719066 :   s = labs(x[1]);
    1352    46518398 :   for (i = 2; i < l && s != 1; i++) s = ugcd(s, labs(x[i]));
    1353    20720422 :   return s;
    1354             : }
    1355             : GEN
    1356      297319 : ZV_content(GEN x)
    1357             : {
    1358      297319 :   long i, l = lg(x);
    1359      297319 :   pari_sp av = avma;
    1360             :   GEN c;
    1361      297319 :   if (l == 1) return gen_0;
    1362      297319 :   if (l == 2) return absi(gel(x,1));
    1363      204793 :   c = gel(x,1);
    1364      557863 :   for (i = 2; i < l; i++)
    1365             :   {
    1366      403566 :     c = gcdii(c, gel(x,i));
    1367      403566 :     if (is_pm1(c)) { set_avma(av); return gen_1; }
    1368             :   }
    1369      154297 :   return gerepileuptoint(av, c);
    1370             : }
    1371             : 
    1372             : GEN
    1373     3868248 : ZM_det_triangular(GEN mat)
    1374             : {
    1375             :   pari_sp av;
    1376     3868248 :   long i,l = lg(mat);
    1377             :   GEN s;
    1378             : 
    1379     3868248 :   if (l<3) return l<2? gen_1: icopy(gcoeff(mat,1,1));
    1380     3465226 :   av = avma; s = gcoeff(mat,1,1);
    1381     9398273 :   for (i=2; i<l; i++) s = mulii(s,gcoeff(mat,i,i));
    1382     3464821 :   return gerepileuptoint(av,s);
    1383             : }
    1384             : 
    1385             : /* assumes no overflow */
    1386             : long
    1387      944593 : zv_prod(GEN v)
    1388             : {
    1389      944593 :   long n, i, l = lg(v);
    1390      944593 :   if (l == 1) return 1;
    1391      954677 :   n = v[1]; for (i = 2; i < l; i++) n *= v[i];
    1392      766261 :   return n;
    1393             : }
    1394             : 
    1395             : static GEN
    1396   318095601 : _mulii(void *E, GEN a, GEN b)
    1397   318095601 : { (void) E; return mulii(a, b); }
    1398             : 
    1399             : /* product of ulongs */
    1400             : GEN
    1401     1864328 : zv_prod_Z(GEN v)
    1402             : {
    1403             :   pari_sp av;
    1404     1864328 :   long k, m, n = lg(v)-1;
    1405     1864328 :   int stop = 0;
    1406             :   GEN V;
    1407     1864328 :   switch(n) {
    1408       20902 :     case 0: return gen_1;
    1409      129878 :     case 1: return utoi(v[1]);
    1410      976693 :     case 2: return muluu(v[1], v[2]);
    1411             :   }
    1412      736855 :   av = avma; m = n >> 1;
    1413      736855 :   V = cgetg(m + (odd(n)? 2: 1), t_VEC);
    1414   151571292 :   for (k = n; k; k--) /* start from the end: v is usually sorted */
    1415   150835479 :     if (v[k] & HIGHMASK) { stop = 1; break; }
    1416     2422592 :   while (!stop)
    1417             :   { /* HACK: handle V as a t_VECSMALL; gain a few iterations */
    1418    81551657 :     for (k = 1; k <= m; k++)
    1419             :     {
    1420    79256784 :       V[k] = uel(v,k<<1) * uel(v,(k<<1)-1);
    1421    79256784 :       if (V[k] & HIGHMASK) stop = 1; /* last "free" iteration */
    1422             :     }
    1423     2294873 :     if (odd(n))
    1424             :     {
    1425     1352302 :       if (n == 1) { set_avma(av); return utoi(v[1]); }
    1426      743166 :       V[++m] = v[n];
    1427             :     }
    1428     1685737 :     v = V; n = m; m = n >> 1;
    1429             :   }
    1430             :   /* n > 1; m > 0 */
    1431      127719 :   if (n == 2) { set_avma(av); return muluu(v[1], v[2]); }
    1432    46319405 :   for (k = 1; k <= m; k++) gel(V,k) = muluu(v[k<<1], v[(k<<1)-1]);
    1433       87482 :   if (odd(n)) gel(V, ++m) = utoipos(v[n]);
    1434       87481 :   setlg(V, m+1); /* HACK: now V is a bona fide t_VEC */
    1435       87481 :   return gerepileuptoint(av, gen_product(V, NULL, &_mulii));
    1436             : }
    1437             : GEN
    1438    14694393 : vecsmall_prod(GEN v)
    1439             : {
    1440    14694393 :   pari_sp av = avma;
    1441    14694393 :   long k, m, n = lg(v)-1;
    1442             :   GEN V;
    1443    14694393 :   switch (n) {
    1444           0 :     case 0: return gen_1;
    1445           0 :     case 1: return stoi(v[1]);
    1446          21 :     case 2: return mulss(v[1], v[2]);
    1447             :   }
    1448    14694372 :   m = n >> 1;
    1449    14694372 :   V = cgetg(m + (odd(n)? 2: 1), t_VEC);
    1450   161556906 :   for (k = 1; k <= m; k++) gel(V,k) = mulss(v[k<<1], v[(k<<1)-1]);
    1451    14694372 :   if (odd(n)) gel(V,k) = stoi(v[n]);
    1452    14694372 :   return gerepileuptoint(av, gen_product(V, NULL, &_mulii));
    1453             : }
    1454             : 
    1455             : GEN
    1456     8819016 : ZV_prod(GEN v)
    1457             : {
    1458     8819016 :   pari_sp av = avma;
    1459     8819016 :   long i, l = lg(v);
    1460             :   GEN n;
    1461     8819016 :   if (l == 1) return gen_1;
    1462     8635072 :   if (l > 7) return gerepileuptoint(av, gen_product(v, NULL, _mulii));
    1463     1292294 :   n = gel(v,1);
    1464     1292294 :   if (l == 2) return icopy(n);
    1465     2090985 :   for (i = 2; i < l; i++) n = mulii(n, gel(v,i));
    1466      840297 :   return gerepileuptoint(av, n);
    1467             : }
    1468             : /* assumes no overflow */
    1469             : long
    1470       16430 : zv_sum(GEN v)
    1471             : {
    1472       16430 :   long n, i, l = lg(v);
    1473       16430 :   if (l == 1) return 0;
    1474       95678 :   n = v[1]; for (i = 2; i < l; i++) n += v[i];
    1475       16409 :   return n;
    1476             : }
    1477             : /* assumes no overflow and 0 <= n <= #v */
    1478             : long
    1479           0 : zv_sumpart(GEN v, long n)
    1480             : {
    1481             :   long i, p;
    1482           0 :   if (!n) return 0;
    1483           0 :   p = v[1]; for (i = 2; i <= n; i++) p += v[i];
    1484           0 :   return p;
    1485             : }
    1486             : GEN
    1487          77 : ZV_sum(GEN v)
    1488             : {
    1489          77 :   pari_sp av = avma;
    1490          77 :   long i, l = lg(v);
    1491             :   GEN n;
    1492          77 :   if (l == 1) return gen_0;
    1493          77 :   n = gel(v,1);
    1494          77 :   if (l == 2) return icopy(n);
    1495         581 :   for (i = 2; i < l; i++) n = addii(n, gel(v,i));
    1496          77 :   return gerepileuptoint(av, n);
    1497             : }
    1498             : 
    1499             : /********************************************************************/
    1500             : /**                                                                **/
    1501             : /**         GRAM SCHMIDT REDUCTION (integer matrices)              **/
    1502             : /**                                                                **/
    1503             : /********************************************************************/
    1504             : 
    1505             : /* L[k,] += q * L[l,], l < k. Inefficient if q = 0 */
    1506             : static void
    1507      362302 : Zupdate_row(long k, long l, GEN q, GEN L, GEN B)
    1508             : {
    1509      362302 :   long i, qq = itos_or_0(q);
    1510      362303 :   if (!qq)
    1511             :   {
    1512       33377 :     for(i=1;i<l;i++)  gcoeff(L,k,i) = addii(gcoeff(L,k,i),mulii(q,gcoeff(L,l,i)));
    1513        7069 :     gcoeff(L,k,l) = addii(gcoeff(L,k,l), mulii(q,B));
    1514        7069 :     return;
    1515             :   }
    1516      355234 :   if (qq == 1) {
    1517      148810 :     for (i=1;i<l; i++) gcoeff(L,k,i) = addii(gcoeff(L,k,i),gcoeff(L,l,i));
    1518      102096 :     gcoeff(L,k,l) = addii(gcoeff(L,k,l), B);
    1519      253138 :   } else if (qq == -1) {
    1520      151504 :     for (i=1;i<l; i++) gcoeff(L,k,i) = subii(gcoeff(L,k,i),gcoeff(L,l,i));
    1521       89217 :     gcoeff(L,k,l) = addii(gcoeff(L,k,l), negi(B));
    1522             :   } else {
    1523      289618 :     for(i=1;i<l;i++) gcoeff(L,k,i) = addii(gcoeff(L,k,i),mulsi(qq,gcoeff(L,l,i)));
    1524      163931 :     gcoeff(L,k,l) = addii(gcoeff(L,k,l), mulsi(qq,B));
    1525             :   }
    1526             : }
    1527             : 
    1528             : /* update L[k,] */
    1529             : static void
    1530     1033367 : ZRED(long k, long l, GEN x, GEN L, GEN B)
    1531             : {
    1532     1033367 :   GEN q = truedivii(addii(B,shifti(gcoeff(L,k,l),1)), shifti(B,1));
    1533     1033477 :   if (!signe(q)) return;
    1534      362278 :   q = negi(q);
    1535      362297 :   Zupdate_row(k,l,q,L,B);
    1536      362282 :   gel(x,k) = ZC_lincomb(gen_1, q, gel(x,k), gel(x,l));
    1537             : }
    1538             : 
    1539             : /* Gram-Schmidt reduction, x a ZM */
    1540             : static void
    1541     1183949 : ZincrementalGS(GEN x, GEN L, GEN B, long k)
    1542             : {
    1543             :   long i, j;
    1544     3780514 :   for (j=1; j<=k; j++)
    1545             :   {
    1546     2597425 :     pari_sp av = avma;
    1547     2597425 :     GEN u = ZV_dotproduct(gel(x,k), gel(x,j));
    1548     5612765 :     for (i=1; i<j; i++)
    1549             :     {
    1550     3016958 :       u = subii(mulii(gel(B,i+1), u), mulii(gcoeff(L,k,i), gcoeff(L,j,i)));
    1551     3016357 :       u = diviiexact(u, gel(B,i));
    1552             :     }
    1553     2595807 :     gcoeff(L,k,j) = gerepileuptoint(av, u);
    1554             :   }
    1555     1183089 :   gel(B,k+1) = gcoeff(L,k,k); gcoeff(L,k,k) = gen_1;
    1556     1183089 : }
    1557             : 
    1558             : /* Variant reducemodinvertible(ZC v, ZM y), when y singular.
    1559             :  * Very inefficient if y is not LLL-reduced of maximal rank */
    1560             : static GEN
    1561      110484 : ZC_reducemodmatrix_i(GEN v, GEN y)
    1562             : {
    1563      110484 :   GEN B, L, x = shallowconcat(y, v);
    1564      110484 :   long k, lx = lg(x), nx = lx-1;
    1565             : 
    1566      110484 :   B = scalarcol_shallow(gen_1, lx);
    1567      110482 :   L = zeromatcopy(nx, nx);
    1568      454504 :   for (k=1; k <= nx; k++) ZincrementalGS(x, L, B, k);
    1569      344019 :   for (k = nx-1; k >= 1; k--) ZRED(nx,k, x,L,gel(B,k+1));
    1570      110466 :   return gel(x,nx);
    1571             : }
    1572             : GEN
    1573      110484 : ZC_reducemodmatrix(GEN v, GEN y) {
    1574      110484 :   pari_sp av = avma;
    1575      110484 :   return gerepilecopy(av, ZC_reducemodmatrix_i(v,y));
    1576             : }
    1577             : 
    1578             : /* Variant reducemodinvertible(ZM v, ZM y), when y singular.
    1579             :  * Very inefficient if y is not LLL-reduced of maximal rank */
    1580             : static GEN
    1581      226902 : ZM_reducemodmatrix_i(GEN v, GEN y)
    1582             : {
    1583             :   GEN B, L, V;
    1584      226902 :   long j, k, lv = lg(v), nx = lg(y), lx = nx+1;
    1585             : 
    1586      226902 :   V = cgetg(lv, t_MAT);
    1587      226922 :   B = scalarcol_shallow(gen_1, lx);
    1588      226925 :   L = zeromatcopy(nx, nx);
    1589      601057 :   for (k=1; k < nx; k++) ZincrementalGS(y, L, B, k);
    1590      692738 :   for (j = 1; j < lg(v); j++)
    1591             :   {
    1592      465853 :     GEN x = shallowconcat(y, gel(v,j));
    1593      465870 :     ZincrementalGS(x, L, B, nx); /* overwrite last */
    1594     1265700 :     for (k = nx-1; k >= 1; k--) ZRED(nx,k, x,L,gel(B,k+1));
    1595      465825 :     gel(V,j) = gel(x,nx);
    1596             :   }
    1597      226885 :   return V;
    1598             : }
    1599             : GEN
    1600      226901 : ZM_reducemodmatrix(GEN v, GEN y) {
    1601      226901 :   pari_sp av = avma;
    1602      226901 :   return gerepilecopy(av, ZM_reducemodmatrix_i(v,y));
    1603             : }
    1604             : 
    1605             : GEN
    1606       98541 : ZC_reducemodlll(GEN x,GEN y)
    1607             : {
    1608       98541 :   pari_sp av = avma;
    1609       98541 :   GEN z = ZC_reducemodmatrix(x, ZM_lll(y, 0.75, LLL_INPLACE));
    1610       98543 :   return gerepilecopy(av, z);
    1611             : }
    1612             : GEN
    1613           0 : ZM_reducemodlll(GEN x,GEN y)
    1614             : {
    1615           0 :   pari_sp av = avma;
    1616           0 :   GEN z = ZM_reducemodmatrix(x, ZM_lll(y, 0.75, LLL_INPLACE));
    1617           0 :   return gerepilecopy(av, z);
    1618             : }

Generated by: LCOV version 1.16