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.1 lcov report (development 28695-49bb1ac00f) Lines: 779 875 89.0 %
Date: 2023-09-24 07:47:42 Functions: 122 132 92.4 %
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     1723195 : check_ZV(GEN x, long l)
      20             : {
      21             :   long i;
      22    10247414 :   for (i=1; i<l; i++)
      23     8524310 :     if (typ(gel(x,i)) != t_INT) return 0;
      24     1723104 :   return 1;
      25             : }
      26             : void
      27     1415239 : RgV_check_ZV(GEN A, const char *s)
      28             : {
      29     1415239 :   if (!RgV_is_ZV(A)) pari_err_TYPE(stack_strcat(s," [integer vector]"), A);
      30     1415233 : }
      31             : void
      32      599031 : RgM_check_ZM(GEN A, const char *s)
      33             : {
      34      599031 :   long n = lg(A);
      35      599031 :   if (n != 1)
      36             :   {
      37      598884 :     long j, m = lgcols(A);
      38     2321988 :     for (j=1; j<n; j++)
      39     1723194 :       if (!check_ZV(gel(A,j), m))
      40          91 :         pari_err_TYPE(stack_strcat(s," [integer matrix]"), A);
      41             :   }
      42      598941 : }
      43             : 
      44             : static long
      45    68386521 : ZV_max_lg_i(GEN x, long m)
      46             : {
      47    68386521 :   long i, prec = 2;
      48   653955392 :   for (i = 1; i < m; i++) prec = maxss(prec, lgefint(gel(x,i)));
      49    68386898 :   return prec;
      50             : }
      51             : long
      52       10619 : ZV_max_lg(GEN x) { return ZV_max_lg_i(x, lg(x)); }
      53             : 
      54             : static long
      55    17824259 : ZM_max_lg_i(GEN x, long n, long m)
      56             : {
      57    17824259 :   long j, prec = 2;
      58    86199805 :   for (j = 1; j < n; j++) prec = maxss(prec, ZV_max_lg_i(gel(x,j), m));
      59    17824580 :   return prec;
      60             : }
      61             : long
      62     1594571 : ZM_max_lg(GEN x)
      63             : {
      64     1594571 :   long n = lg(x);
      65     1594571 :   if (n == 1) return 2;
      66     1594571 :   return ZM_max_lg_i(x, n, lgcols(x));
      67             : }
      68             : 
      69             : GEN
      70        3784 : ZM_supnorm(GEN x)
      71             : {
      72        3784 :   long i, j, h, lx = lg(x);
      73        3784 :   GEN s = gen_0;
      74        3784 :   if (lx == 1) return gen_1;
      75        3784 :   h = lgcols(x);
      76       23211 :   for (j=1; j<lx; j++)
      77             :   {
      78       19427 :     GEN xj = gel(x,j);
      79      272114 :     for (i=1; i<h; i++)
      80             :     {
      81      252687 :       GEN c = gel(xj,i);
      82      252687 :       if (abscmpii(c, s) > 0) s = c;
      83             :     }
      84             :   }
      85        3784 :   return absi(s);
      86             : }
      87             : 
      88             : /********************************************************************/
      89             : /**                                                                **/
      90             : /**                           MULTIPLICATION                       **/
      91             : /**                                                                **/
      92             : /********************************************************************/
      93             : /* x nonempty ZM, y a compatible nc (dimension > 0). */
      94             : static GEN
      95           0 : ZM_nc_mul_i(GEN x, GEN y, long c, long l)
      96             : {
      97             :   long i, j;
      98             :   pari_sp av;
      99           0 :   GEN z = cgetg(l,t_COL), s;
     100             : 
     101           0 :   for (i=1; i<l; i++)
     102             :   {
     103           0 :     av = avma; s = muliu(gcoeff(x,i,1),y[1]);
     104           0 :     for (j=2; j<c; j++)
     105           0 :       if (y[j]) s = addii(s, muliu(gcoeff(x,i,j),y[j]));
     106           0 :     gel(z,i) = gerepileuptoint(av,s);
     107             :   }
     108           0 :   return z;
     109             : }
     110             : 
     111             : /* x ZV, y a compatible zc. */
     112             : GEN
     113     2214660 : ZV_zc_mul(GEN x, GEN y)
     114             : {
     115     2214660 :   long j, l = lg(x);
     116     2214660 :   pari_sp av = avma;
     117     2214660 :   GEN s = mulis(gel(x,1),y[1]);
     118    50028650 :   for (j=2; j<l; j++)
     119    47813990 :     if (y[j]) s = addii(s, mulis(gel(x,j),y[j]));
     120     2214660 :   return gerepileuptoint(av,s);
     121             : }
     122             : 
     123             : /* x nonempty ZM, y a compatible zc (dimension > 0). */
     124             : static GEN
     125    23054836 : ZM_zc_mul_i(GEN x, GEN y, long c, long l)
     126             : {
     127             :   long i, j;
     128    23054836 :   GEN z = cgetg(l,t_COL);
     129             : 
     130   133470396 :   for (i=1; i<l; i++)
     131             :   {
     132   110418424 :     pari_sp av = avma;
     133   110418424 :     GEN s = mulis(gcoeff(x,i,1),y[1]);
     134  1195963634 :     for (j=2; j<c; j++)
     135  1085551083 :       if (y[j]) s = addii(s, mulis(gcoeff(x,i,j),y[j]));
     136   110412551 :     gel(z,i) = gerepileuptoint(av,s);
     137             :   }
     138    23051972 :   return z;
     139             : }
     140             : GEN
     141    21344912 : ZM_zc_mul(GEN x, GEN y) {
     142    21344912 :   long l = lg(x);
     143    21344912 :   if (l == 1) return cgetg(1, t_COL);
     144    21344912 :   return ZM_zc_mul_i(x,y, l, lgcols(x));
     145             : }
     146             : 
     147             : /* y nonempty ZM, x a compatible zv (dimension > 0). */
     148             : GEN
     149        1624 : zv_ZM_mul(GEN x, GEN y) {
     150        1624 :   long i,j, lx = lg(x), ly = lg(y);
     151             :   GEN z;
     152        1624 :   if (lx == 1) return zerovec(ly-1);
     153        1624 :   z = cgetg(ly,t_VEC);
     154        3822 :   for (j=1; j<ly; j++)
     155             :   {
     156        2198 :     pari_sp av = avma;
     157        2198 :     GEN s = mulsi(x[1], gcoeff(y,1,j));
     158        3878 :     for (i=2; i<lx; i++)
     159        1680 :       if (x[i]) s = addii(s, mulsi(x[i], gcoeff(y,i,j)));
     160        2198 :     gel(z,j) = gerepileuptoint(av,s);
     161             :   }
     162        1624 :   return z;
     163             : }
     164             : 
     165             : /* x ZM, y a compatible zm (dimension > 0). */
     166             : GEN
     167      806810 : ZM_zm_mul(GEN x, GEN y)
     168             : {
     169      806810 :   long j, c, l = lg(x), ly = lg(y);
     170      806810 :   GEN z = cgetg(ly, t_MAT);
     171      806810 :   if (l == 1) return z;
     172      806803 :   c = lgcols(x);
     173     2516946 :   for (j = 1; j < ly; j++) gel(z,j) = ZM_zc_mul_i(x, gel(y,j), l,c);
     174      806803 :   return z;
     175             : }
     176             : /* x ZM, y a compatible zn (dimension > 0). */
     177             : GEN
     178           0 : ZM_nm_mul(GEN x, GEN y)
     179             : {
     180           0 :   long j, c, l = lg(x), ly = lg(y);
     181           0 :   GEN z = cgetg(ly, t_MAT);
     182           0 :   if (l == 1) return z;
     183           0 :   c = lgcols(x);
     184           0 :   for (j = 1; j < ly; j++) gel(z,j) = ZM_nc_mul_i(x, gel(y,j), l,c);
     185           0 :   return z;
     186             : }
     187             : 
     188             : /* Strassen-Winograd algorithm */
     189             : 
     190             : /* Return A[ma+1..ma+da, na+1..na+ea] - B[mb+1..mb+db, nb+1..nb+eb]
     191             :  * as an (m x n)-matrix, padding the input with zeroes as necessary. */
     192             : static GEN
     193      274352 : add_slices(long m, long n,
     194             :            GEN A, long ma, long da, long na, long ea,
     195             :            GEN B, long mb, long db, long nb, long eb)
     196             : {
     197      274352 :   long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
     198      274352 :   GEN M = cgetg(n + 1, t_MAT), C;
     199             : 
     200     2607907 :   for (j = 1; j <= min_e; j++) {
     201     2333555 :     gel(M, j) = C = cgetg(m + 1, t_COL);
     202    41208305 :     for (i = 1; i <= min_d; i++)
     203    38874750 :       gel(C, i) = addii(gcoeff(A, ma + i, na + j),
     204    38874750 :                         gcoeff(B, mb + i, nb + j));
     205     2362960 :     for (; i <= da; i++) gel(C, i) = gcoeff(A, ma + i, na + j);
     206     2333555 :     for (; i <= db; i++) gel(C, i) = gcoeff(B, mb + i, nb + j);
     207     2333555 :     for (; i <= m; i++)  gel(C, i) = gen_0;
     208             :   }
     209      293470 :   for (; j <= ea; j++) {
     210       19118 :     gel(M, j) = C = cgetg(m + 1, t_COL);
     211       83048 :     for (i = 1; i <= da; i++) gel(C, i) = gcoeff(A, ma + i, na + j);
     212       19118 :     for (; i <= m; i++) gel(C, i) = gen_0;
     213             :   }
     214      274352 :   for (; j <= eb; j++) {
     215           0 :     gel(M, j) = C = cgetg(m + 1, t_COL);
     216           0 :     for (i = 1; i <= db; i++) gel(C, i) = gcoeff(B, mb + i, nb + j);
     217           0 :     for (; i <= m; i++) gel(C, i) = gen_0;
     218             :   }
     219      274352 :   for (; j <= n; j++) gel(M, j) = zerocol(m);
     220      274352 :   return M;
     221             : }
     222             : 
     223             : /* Return A[ma+1..ma+da, na+1..na+ea] - B[mb+1..mb+db, nb+1..nb+eb]
     224             :  * as an (m x n)-matrix, padding the input with zeroes as necessary. */
     225             : static GEN
     226      240058 : subtract_slices(long m, long n,
     227             :                 GEN A, long ma, long da, long na, long ea,
     228             :                 GEN B, long mb, long db, long nb, long eb)
     229             : {
     230      240058 :   long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
     231      240058 :   GEN M = cgetg(n + 1, t_MAT), C;
     232             : 
     233     2271095 :   for (j = 1; j <= min_e; j++) {
     234     2031037 :     gel(M, j) = C = cgetg(m + 1, t_COL);
     235    36710754 :     for (i = 1; i <= min_d; i++)
     236    34679717 :       gel(C, i) = subii(gcoeff(A, ma + i, na + j),
     237    34679717 :                         gcoeff(B, mb + i, nb + j));
     238     2050251 :     for (; i <= da; i++) gel(C, i) = gcoeff(A, ma + i, na + j);
     239     2062896 :     for (; i <= db; i++) gel(C, i) = negi(gcoeff(B, mb + i, nb + j));
     240     2031037 :     for (; i <= m; i++) gel(C, i) = gen_0;
     241             :   }
     242      240058 :   for (; j <= ea; j++) {
     243           0 :     gel(M, j) = C = cgetg(m + 1, t_COL);
     244           0 :     for (i = 1; i <= da; i++) gel(C, i) = gcoeff(A, ma + i, na + j);
     245           0 :     for (; i <= m; i++) gel(C, i) = gen_0;
     246             :   }
     247      254964 :   for (; j <= eb; j++) {
     248       14906 :     gel(M, j) = C = cgetg(m + 1, t_COL);
     249       46788 :     for (i = 1; i <= db; i++) gel(C, i) = negi(gcoeff(B, mb + i, nb + j));
     250       14906 :     for (; i <= m; i++) gel(C, i) = gen_0;
     251             :   }
     252      254964 :   for (; j <= n; j++) gel(M, j) = zerocol(m);
     253      240058 :   return M;
     254             : }
     255             : 
     256             : static GEN ZM_mul_i(GEN x, GEN y, long l, long lx, long ly);
     257             : 
     258             : /* Strassen-Winograd matrix product A (m x n) * B (n x p) */
     259             : static GEN
     260       34294 : ZM_mul_sw(GEN A, GEN B, long m, long n, long p)
     261             : {
     262       34294 :   pari_sp av = avma;
     263       34294 :   long m1 = (m + 1)/2, m2 = m/2,
     264       34294 :     n1 = (n + 1)/2, n2 = n/2,
     265       34294 :     p1 = (p + 1)/2, p2 = p/2;
     266             :   GEN A11, A12, A22, B11, B21, B22,
     267             :     S1, S2, S3, S4, T1, T2, T3, T4,
     268             :     M1, M2, M3, M4, M5, M6, M7,
     269             :     V1, V2, V3, C11, C12, C21, C22, C;
     270             : 
     271       34294 :   T2 = subtract_slices(n1, p2, B, 0, n1, p1, p2, B, n1, n2, p1, p2);
     272       34294 :   S1 = subtract_slices(m2, n1, A, m1, m2, 0, n1, A, 0, m2, 0, n1);
     273       34294 :   M2 = ZM_mul_i(S1, T2, m2 + 1, n1 + 1, p2 + 1);
     274       34294 :   if (gc_needed(av, 1))
     275           0 :     gerepileall(av, 2, &T2, &M2);  /* destroy S1 */
     276       34294 :   T3 = subtract_slices(n1, p1, T2, 0, n1, 0, p2, B, 0, n1, 0, p1);
     277       34294 :   if (gc_needed(av, 1))
     278           0 :     gerepileall(av, 2, &M2, &T3);  /* destroy T2 */
     279       34294 :   S2 = add_slices(m2, n1, A, m1, m2, 0, n1, A, m1, m2, n1, n2);
     280       34294 :   T1 = subtract_slices(n1, p1, B, 0, n1, p1, p2, B, 0, n1, 0, p2);
     281       34294 :   M3 = ZM_mul_i(S2, T1, m2 + 1, n1 + 1, p2 + 1);
     282       34294 :   if (gc_needed(av, 1))
     283           0 :     gerepileall(av, 4, &M2, &T3, &S2, &M3);  /* destroy T1 */
     284       34294 :   S3 = subtract_slices(m1, n1, S2, 0, m2, 0, n1, A, 0, m1, 0, n1);
     285       34294 :   if (gc_needed(av, 1))
     286           0 :     gerepileall(av, 4, &M2, &T3, &M3, &S3);  /* destroy S2 */
     287       34294 :   A11 = matslice(A, 1, m1, 1, n1);
     288       34294 :   B11 = matslice(B, 1, n1, 1, p1);
     289       34294 :   M1 = ZM_mul_i(A11, B11, m1 + 1, n1 + 1, p1 + 1);
     290       34294 :   if (gc_needed(av, 1))
     291           0 :     gerepileall(av, 5, &M2, &T3, &M3, &S3, &M1);  /* destroy A11, B11 */
     292       34294 :   A12 = matslice(A, 1, m1, n1 + 1, n);
     293       34294 :   B21 = matslice(B, n1 + 1, n, 1, p1);
     294       34294 :   M4 = ZM_mul_i(A12, B21, m1 + 1, n2 + 1, p1 + 1);
     295       34294 :   if (gc_needed(av, 1))
     296           0 :     gerepileall(av, 6, &M2, &T3, &M3, &S3, &M1, &M4);  /* destroy A12, B21 */
     297       34294 :   C11 = add_slices(m1, p1, M1, 0, m1, 0, p1, M4, 0, m1, 0, p1);
     298       34294 :   if (gc_needed(av, 1))
     299           0 :     gerepileall(av, 6, &M2, &T3, &M3, &S3, &M1, &C11);  /* destroy M4 */
     300       34294 :   M5 = ZM_mul_i(S3, T3, m1 + 1, n1 + 1, p1 + 1);
     301       34294 :   S4 = subtract_slices(m1, n2, A, 0, m1, n1, n2, S3, 0, m1, 0, n2);
     302       34294 :   if (gc_needed(av, 1))
     303           5 :     gerepileall(av, 7, &M2, &T3, &M3, &M1, &C11, &M5, &S4);  /* destroy S3 */
     304       34294 :   T4 = add_slices(n2, p1, B, n1, n2, 0, p1, T3, 0, n2, 0, p1);
     305       34294 :   if (gc_needed(av, 1))
     306           0 :     gerepileall(av, 7, &M2, &M3, &M1, &C11, &M5, &S4, &T4);  /* destroy T3 */
     307       34294 :   V1 = subtract_slices(m1, p1, M1, 0, m1, 0, p1, M5, 0, m1, 0, p1);
     308       34294 :   if (gc_needed(av, 1))
     309           1 :     gerepileall(av, 6, &M2, &M3, &S4, &T4, &C11, &V1);  /* destroy M1, M5 */
     310       34294 :   B22 = matslice(B, n1 + 1, n, p1 + 1, p);
     311       34294 :   M6 = ZM_mul_i(S4, B22, m1 + 1, n2 + 1, p2 + 1);
     312       34294 :   if (gc_needed(av, 1))
     313           6 :     gerepileall(av, 6, &M2, &M3, &T4, &C11, &V1, &M6);  /* destroy S4, B22 */
     314       34294 :   A22 = matslice(A, m1 + 1, m, n1 + 1, n);
     315       34294 :   M7 = ZM_mul_i(A22, T4, m2 + 1, n2 + 1, p1 + 1);
     316       34294 :   if (gc_needed(av, 1))
     317           0 :     gerepileall(av, 6, &M2, &M3, &C11, &V1, &M6, &M7);  /* destroy A22, T4 */
     318       34294 :   V3 = add_slices(m1, p2, V1, 0, m1, 0, p2, M3, 0, m2, 0, p2);
     319       34294 :   C12 = add_slices(m1, p2, V3, 0, m1, 0, p2, M6, 0, m1, 0, p2);
     320       34294 :   if (gc_needed(av, 1))
     321          12 :     gerepileall(av, 6, &M2, &M3, &C11, &V1, &M7, &C12);  /* destroy V3, M6 */
     322       34294 :   V2 = add_slices(m2, p1, V1, 0, m2, 0, p1, M2, 0, m2, 0, p2);
     323       34294 :   if (gc_needed(av, 1))
     324           0 :     gerepileall(av, 5, &M3, &C11, &M7, &C12, &V2);  /* destroy V1, M2 */
     325       34294 :   C21 = add_slices(m2, p1, V2, 0, m2, 0, p1, M7, 0, m2, 0, p1);
     326       34294 :   if (gc_needed(av, 1))
     327           6 :     gerepileall(av, 5, &M3, &C11, &C12, &V2, &C21);  /* destroy M7 */
     328       34294 :   C22 = add_slices(m2, p2, V2, 0, m2, 0, p2, M3, 0, m2, 0, p2);
     329       34294 :   if (gc_needed(av, 1))
     330           0 :     gerepileall(av, 4, &C11, &C12, &C21, &C22);  /* destroy V2, M3 */
     331       34294 :   C = shallowconcat(vconcat(C11, C21), vconcat(C12, C22));
     332       34294 :   return gerepilecopy(av, C);
     333             : }
     334             : 
     335             : /* x[i,]*y. Assume lg(x) > 1 and 0 < i < lgcols(x) */
     336             : static GEN
     337   444632540 : ZMrow_ZC_mul_i(GEN x, GEN y, long i, long lx)
     338             : {
     339   444632540 :   pari_sp av = avma;
     340   444632540 :   GEN c = mulii(gcoeff(x,i,1), gel(y,1)), ZERO = gen_0;
     341             :   long k;
     342  4880370806 :   for (k = 2; k < lx; k++)
     343             :   {
     344  4436147661 :     GEN t = mulii(gcoeff(x,i,k), gel(y,k));
     345  4434047926 :     if (t != ZERO) c = addii(c, t);
     346             :   }
     347   444223145 :   return gerepileuptoint(av, c);
     348             : }
     349             : GEN
     350   150688514 : ZMrow_ZC_mul(GEN x, GEN y, long i)
     351   150688514 : { return ZMrow_ZC_mul_i(x, y, i, lg(x)); }
     352             : 
     353             : /* return x * y, 1 < lx = lg(x), l = lgcols(x) */
     354             : static GEN
     355    51528431 : ZM_ZC_mul_i(GEN x, GEN y, long lx, long l)
     356             : {
     357    51528431 :   GEN z = cgetg(l,t_COL);
     358             :   long i;
     359   345489362 :   for (i=1; i<l; i++) gel(z,i) = ZMrow_ZC_mul_i(x,y,i,lx);
     360    51500439 :   return z;
     361             : }
     362             : 
     363             : static GEN
     364     8092590 : ZM_mul_classical(GEN x, GEN y, long l, long lx, long ly)
     365             : {
     366             :   long j;
     367     8092590 :   GEN z = cgetg(ly, t_MAT);
     368    36692520 :   for (j = 1; j < ly; j++)
     369    28601355 :     gel(z, j) = ZM_ZC_mul_i(x, gel(y, j), lx, l);
     370     8091165 :   return z;
     371             : }
     372             : 
     373             : static GEN
     374        1290 : ZM_mul_slice(GEN A, GEN B, GEN P, GEN *mod)
     375             : {
     376        1290 :   pari_sp av = avma;
     377        1290 :   long i, n = lg(P)-1;
     378             :   GEN H, T;
     379        1290 :   if (n == 1)
     380             :   {
     381           0 :     ulong p = uel(P,1);
     382           0 :     GEN a = ZM_to_Flm(A, p);
     383           0 :     GEN b = ZM_to_Flm(B, p);
     384           0 :     GEN Hp = gerepileupto(av, Flm_to_ZM(Flm_mul(a, b, p)));
     385           0 :     *mod = utoi(p); return Hp;
     386             :   }
     387        1290 :   T = ZV_producttree(P);
     388        1290 :   A = ZM_nv_mod_tree(A, P, T);
     389        1290 :   B = ZM_nv_mod_tree(B, P, T);
     390        1290 :   H = cgetg(n+1, t_VEC);
     391        9882 :   for(i=1; i <= n; i++)
     392        8592 :     gel(H,i) = Flm_mul(gel(A,i),gel(B,i),P[i]);
     393        1290 :   H = nmV_chinese_center_tree_seq(H, P, T, ZV_chinesetree(P, T));
     394        1290 :   *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
     395             : }
     396             : 
     397             : GEN
     398        1290 : ZM_mul_worker(GEN P, GEN A, GEN B)
     399             : {
     400        1290 :   GEN V = cgetg(3, t_VEC);
     401        1290 :   gel(V,1) = ZM_mul_slice(A, B, P, &gel(V,2));
     402        1290 :   return V;
     403             : }
     404             : 
     405             : static GEN
     406         919 : ZM_mul_fast(GEN A, GEN B, long lA, long lB, long sA, long sB)
     407             : {
     408         919 :   pari_sp av = avma;
     409             :   forprime_t S;
     410             :   GEN H, worker;
     411             :   long h;
     412         919 :   if (sA == 2 || sB == 2) return zeromat(nbrows(A),lB-1);
     413         905 :   h = 1 + (sA + sB - 4) * BITS_IN_LONG + expu(lA-1);
     414         905 :   init_modular_big(&S);
     415         905 :   worker = snm_closure(is_entry("_ZM_mul_worker"), mkvec2(A,B));
     416         905 :   H = gen_crt("ZM_mul", worker, &S, NULL, h, 0, NULL,
     417             :               nmV_chinese_center, FpM_center);
     418         905 :   return gerepileupto(av, H);
     419             : }
     420             : 
     421             : /* s = min(log_BIL |x|, log_BIL |y|), use Strassen-Winograd when
     422             :  * min(dims) > B */
     423             : static long
     424     8126854 : sw_bound(long s)
     425     8126854 : { return s > 60 ? 2: s > 25 ? 4: s > 15 ? 8 : s > 8 ? 16 : 32; }
     426             : 
     427             : static GEN
     428    15307810 : ZM_mul_i(GEN x, GEN y, long l, long lx, long ly)
     429             : {
     430             :   long sx, sy, B;
     431    15307810 :   if (lx==3 && l==3 && ly==3) return ZM2_mul(x,y);
     432     8101932 :   sx = ZM_max_lg_i(x, lx, l);
     433     8102549 :   sy = ZM_max_lg_i(y, ly, lx);
     434     8102630 :   if ((lx > 70 && ly > 70 && l > 70) && (sx <= 10 * sy && sy <= 10 * sx))
     435         919 :     return ZM_mul_fast(x, y, lx, ly, sx, sy);
     436             : 
     437     8101711 :   B = sw_bound(minss(sx, sy));
     438     8101688 :   if (l <= B || lx <= B || ly <= B)
     439     8067424 :     return ZM_mul_classical(x, y, l, lx, ly);
     440             :   else
     441       34264 :     return ZM_mul_sw(x, y, l - 1, lx - 1, ly - 1);
     442             : }
     443             : 
     444             : GEN
     445    15203665 : ZM_mul(GEN x, GEN y)
     446             : {
     447    15203665 :   long lx = lg(x), ly = lg(y);
     448    15203665 :   if (ly == 1) return cgetg(1,t_MAT);
     449    15068985 :   if (lx == 1) return zeromat(0, ly-1);
     450    15067405 :   return ZM_mul_i(x, y, lgcols(x), lx, ly);
     451             : }
     452             : 
     453             : GEN
     454      538453 : QM_mul(GEN x, GEN y)
     455             : {
     456      538453 :   GEN dx, nx = Q_primitive_part(x, &dx);
     457      538453 :   GEN dy, ny = Q_primitive_part(y, &dy);
     458      538453 :   GEN z = ZM_mul(nx, ny);
     459      538453 :   if (dx || dy)
     460             :   {
     461      452820 :     GEN d = dx ? dy ? gmul(dx, dy): dx : dy;
     462      452820 :     if (!gequal1(d)) z = ZM_Q_mul(z, d);
     463             :   }
     464      538453 :   return z;
     465             : }
     466             : 
     467             : GEN
     468         700 : QM_sqr(GEN x)
     469             : {
     470         700 :   GEN dx, nx = Q_primitive_part(x, &dx);
     471         700 :   GEN z = ZM_sqr(nx);
     472         700 :   if (dx)
     473         700 :     z = ZM_Q_mul(z, gsqr(dx));
     474         700 :   return z;
     475             : }
     476             : 
     477             : GEN
     478      132926 : QM_QC_mul(GEN x, GEN y)
     479             : {
     480      132926 :   GEN dx, nx = Q_primitive_part(x, &dx);
     481      132926 :   GEN dy, ny = Q_primitive_part(y, &dy);
     482      132926 :   GEN z = ZM_ZC_mul(nx, ny);
     483      132926 :   if (dx || dy)
     484             :   {
     485      132898 :     GEN d = dx ? dy ? gmul(dx, dy): dx : dy;
     486      132898 :     if (!gequal1(d)) z = ZC_Q_mul(z, d);
     487             :   }
     488      132926 :   return z;
     489             : }
     490             : 
     491             : /* assume result is symmetric */
     492             : GEN
     493           0 : ZM_multosym(GEN x, GEN y)
     494             : {
     495           0 :   long j, lx, ly = lg(y);
     496             :   GEN M;
     497           0 :   if (ly == 1) return cgetg(1,t_MAT);
     498           0 :   lx = lg(x); /* = lgcols(y) */
     499           0 :   if (lx == 1) return cgetg(1,t_MAT);
     500             :   /* ly = lgcols(x) */
     501           0 :   M = cgetg(ly, t_MAT);
     502           0 :   for (j=1; j<ly; j++)
     503             :   {
     504           0 :     GEN z = cgetg(ly,t_COL), yj = gel(y,j);
     505             :     long i;
     506           0 :     for (i=1; i<j; i++) gel(z,i) = gcoeff(M,j,i);
     507           0 :     for (i=j; i<ly; i++)gel(z,i) = ZMrow_ZC_mul_i(x,yj,i,lx);
     508           0 :     gel(M,j) = z;
     509             :   }
     510           0 :   return M;
     511             : }
     512             : 
     513             : /* compute m*diagonal(d), assume lg(d) = lg(m). Shallow */
     514             : GEN
     515          21 : ZM_mul_diag(GEN m, GEN d)
     516             : {
     517             :   long j, l;
     518          21 :   GEN y = cgetg_copy(m, &l);
     519          77 :   for (j=1; j<l; j++)
     520             :   {
     521          56 :     GEN c = gel(d,j);
     522          56 :     gel(y,j) = equali1(c)? gel(m,j): ZC_Z_mul(gel(m,j), c);
     523             :   }
     524          21 :   return y;
     525             : }
     526             : /* compute diagonal(d)*m, assume lg(d) = lg(m~). Shallow */
     527             : GEN
     528      530461 : ZM_diag_mul(GEN d, GEN m)
     529             : {
     530      530461 :   long i, j, l = lg(d), lm = lg(m);
     531      530461 :   GEN y = cgetg(lm, t_MAT);
     532     1513696 :   for (j=1; j<lm; j++) gel(y,j) = cgetg(l, t_COL);
     533     1639876 :   for (i=1; i<l; i++)
     534             :   {
     535     1109744 :     GEN c = gel(d,i);
     536     1109744 :     if (equali1(c))
     537      236155 :       for (j=1; j<lm; j++) gcoeff(y,i,j) = gcoeff(m,i,j);
     538             :     else
     539     3388730 :       for (j=1; j<lm; j++) gcoeff(y,i,j) = mulii(gcoeff(m,i,j), c);
     540             :   }
     541      530132 :   return y;
     542             : }
     543             : 
     544             : /* assume lx > 1 is lg(x) = lg(y) */
     545             : static GEN
     546    20888467 : ZV_dotproduct_i(GEN x, GEN y, long lx)
     547             : {
     548    20888467 :   pari_sp av = avma;
     549    20888467 :   GEN c = mulii(gel(x,1), gel(y,1));
     550             :   long i;
     551   147501445 :   for (i = 2; i < lx; i++)
     552             :   {
     553   126614779 :     GEN t = mulii(gel(x,i), gel(y,i));
     554   126614615 :     if (t != gen_0) c = addii(c, t);
     555             :   }
     556    20886666 :   return gerepileuptoint(av, c);
     557             : }
     558             : 
     559             : /* x~ * y, assuming result is symmetric */
     560             : GEN
     561      526314 : ZM_transmultosym(GEN x, GEN y)
     562             : {
     563      526314 :   long i, j, l, ly = lg(y);
     564             :   GEN M;
     565      526314 :   if (ly == 1) return cgetg(1,t_MAT);
     566             :   /* lg(x) = ly */
     567      526314 :   l = lgcols(y); /* = lgcols(x) */
     568      526314 :   M = cgetg(ly, t_MAT);
     569     2678336 :   for (i=1; i<ly; i++)
     570             :   {
     571     2152022 :     GEN xi = gel(x,i), c = cgetg(ly,t_COL);
     572     2152022 :     gel(M,i) = c;
     573     7018586 :     for (j=1; j<i; j++)
     574     4866564 :       gcoeff(M,i,j) = gel(c,j) = ZV_dotproduct_i(xi,gel(y,j),l);
     575     2152022 :     gel(c,i) = ZV_dotproduct_i(xi,gel(y,i),l);
     576             :   }
     577      526314 :   return M;
     578             : }
     579             : /* x~ * y */
     580             : GEN
     581        2374 : ZM_transmul(GEN x, GEN y)
     582             : {
     583        2374 :   long i, j, l, lx, ly = lg(y);
     584             :   GEN M;
     585        2374 :   if (ly == 1) return cgetg(1,t_MAT);
     586        2374 :   lx = lg(x);
     587        2374 :   l = lgcols(y);
     588        2374 :   if (lgcols(x) != l) pari_err_OP("operation 'ZM_transmul'", x,y);
     589        2374 :   M = cgetg(ly, t_MAT);
     590        7248 :   for (i=1; i<ly; i++)
     591             :   {
     592        4874 :     GEN yi = gel(y,i), c = cgetg(lx,t_COL);
     593        4874 :     gel(M,i) = c;
     594       12569 :     for (j=1; j<lx; j++) gel(c,j) = ZV_dotproduct_i(yi,gel(x,j),l);
     595             :   }
     596        2374 :   return M;
     597             : }
     598             : 
     599             : static GEN
     600       25186 : ZM_sqr_i(GEN x, long l)
     601             : {
     602       25186 :   long s = ZM_max_lg_i(x, l, l);
     603       25186 :   if (l > 70) return ZM_mul_fast(x, x, l, l, s, s);
     604       25186 :   if (l <= sw_bound(s))
     605       25156 :     return ZM_mul_classical(x, x, l, l, l);
     606             :   else
     607          30 :     return ZM_mul_sw(x, x, l - 1, l - 1, l - 1);
     608             : }
     609             : 
     610             : GEN
     611       25186 : ZM_sqr(GEN x)
     612             : {
     613       25186 :   long lx=lg(x);
     614       25186 :   if (lx==1) return cgetg(1,t_MAT);
     615       25186 :   return ZM_sqr_i(x, lx);
     616             : }
     617             : GEN
     618    23014861 : ZM_ZC_mul(GEN x, GEN y)
     619             : {
     620    23014861 :   long lx = lg(x);
     621    23014861 :   return lx==1? cgetg(1,t_COL): ZM_ZC_mul_i(x, y, lx, lgcols(x));
     622             : }
     623             : 
     624             : GEN
     625     1446932 : ZC_Z_div(GEN x, GEN c)
     626     7120621 : { pari_APPLY_type(t_COL, Qdivii(gel(x,i), c)) }
     627             : 
     628             : GEN
     629       14837 : ZM_Z_div(GEN x, GEN c)
     630      168167 : { pari_APPLY_same(ZC_Z_div(gel(x, i), c)) }
     631             : 
     632             : GEN
     633     2629652 : ZC_Q_mul(GEN A, GEN z)
     634             : {
     635     2629652 :   pari_sp av = avma;
     636     2629652 :   long i, l = lg(A);
     637             :   GEN d, n, Ad, B, u;
     638     2629652 :   if (typ(z)==t_INT) return ZC_Z_mul(A,z);
     639     2627139 :   n = gel(z, 1); d = gel(z, 2);
     640     2627139 :   Ad = FpC_red(A, d);
     641     2627082 :   u = gcdii(d, FpV_factorback(Ad, NULL, d));
     642     2627068 :   B = cgetg(l, t_COL);
     643     2627085 :   if (equali1(u))
     644             :   {
     645      295317 :     for(i=1; i<l; i++)
     646      236557 :       gel(B, i) = mkfrac(mulii(n, gel(A,i)), d);
     647             :   } else
     648             :   {
     649    17631077 :     for(i=1; i<l; i++)
     650             :     {
     651    15062844 :       GEN di = gcdii(gel(Ad, i), u), ni = mulii(n, diviiexact(gel(A,i), di));
     652    15062723 :       if (equalii(d, di))
     653    10420598 :         gel(B, i) = ni;
     654             :       else
     655     4642147 :         gel(B, i) = mkfrac(ni, diviiexact(d, di));
     656             :     }
     657             :   }
     658     2626993 :   return gerepilecopy(av, B);
     659             : }
     660             : 
     661             : GEN
     662     1095530 : ZM_Q_mul(GEN x, GEN z)
     663             : {
     664     1095530 :   if (typ(z)==t_INT) return ZM_Z_mul(x,z);
     665     3162372 :   pari_APPLY_same(ZC_Q_mul(gel(x, i), z));
     666             : }
     667             : 
     668             : long
     669   196399119 : zv_dotproduct(GEN x, GEN y)
     670             : {
     671   196399119 :   long i, lx = lg(x);
     672             :   ulong c;
     673   196399119 :   if (lx == 1) return 0;
     674   196399119 :   c = uel(x,1)*uel(y,1);
     675  3047359154 :   for (i = 2; i < lx; i++)
     676  2850960035 :     c += uel(x,i)*uel(y,i);
     677   196399119 :   return c;
     678             : }
     679             : 
     680             : GEN
     681      212254 : ZV_ZM_mul(GEN x, GEN y)
     682             : {
     683      212254 :   long i, lx = lg(x), ly = lg(y);
     684             :   GEN z;
     685      212254 :   if (lx == 1) return zerovec(ly-1);
     686      212135 :   z = cgetg(ly, t_VEC);
     687      843143 :   for (i = 1; i < ly; i++) gel(z,i) = ZV_dotproduct_i(x, gel(y,i), lx);
     688      212135 :   return z;
     689             : }
     690             : 
     691             : GEN
     692           0 : ZC_ZV_mul(GEN x, GEN y)
     693             : {
     694           0 :   long i,j, lx=lg(x), ly=lg(y);
     695             :   GEN z;
     696           0 :   if (ly==1) return cgetg(1,t_MAT);
     697           0 :   z = cgetg(ly,t_MAT);
     698           0 :   for (j=1; j < ly; j++)
     699             :   {
     700           0 :     gel(z,j) = cgetg(lx,t_COL);
     701           0 :     for (i=1; i<lx; i++) gcoeff(z,i,j) = mulii(gel(x,i),gel(y,j));
     702             :   }
     703           0 :   return z;
     704             : }
     705             : 
     706             : GEN
     707    10119216 : ZV_dotsquare(GEN x)
     708             : {
     709             :   long i, lx;
     710             :   pari_sp av;
     711             :   GEN z;
     712    10119216 :   lx = lg(x);
     713    10119216 :   if (lx == 1) return gen_0;
     714    10119216 :   av = avma; z = sqri(gel(x,1));
     715    33334084 :   for (i=2; i<lx; i++) z = addii(z, sqri(gel(x,i)));
     716    10116283 :   return gerepileuptoint(av,z);
     717             : }
     718             : 
     719             : GEN
     720    19740843 : ZV_dotproduct(GEN x,GEN y)
     721             : {
     722             :   long lx;
     723    19740843 :   if (x == y) return ZV_dotsquare(x);
     724    13230314 :   lx = lg(x);
     725    13230314 :   if (lx == 1) return gen_0;
     726    13230314 :   return ZV_dotproduct_i(x, y, lx);
     727             : }
     728             : 
     729             : static GEN
     730         294 : _ZM_mul(void *data /*ignored*/, GEN x, GEN y)
     731         294 : { (void)data; return ZM_mul(x,y); }
     732             : static GEN
     733       23205 : _ZM_sqr(void *data /*ignored*/, GEN x)
     734       23205 : { (void)data; return ZM_sqr(x); }
     735             : GEN
     736           0 : ZM_pow(GEN x, GEN n)
     737             : {
     738           0 :   pari_sp av = avma;
     739           0 :   if (!signe(n)) return matid(lg(x)-1);
     740           0 :   return gerepilecopy(av, gen_pow_i(x, n, NULL, &_ZM_sqr, &_ZM_mul));
     741             : }
     742             : GEN
     743       22659 : ZM_powu(GEN x, ulong n)
     744             : {
     745       22659 :   pari_sp av = avma;
     746       22659 :   if (!n) return matid(lg(x)-1);
     747       22659 :   return gerepilecopy(av, gen_powu_i(x, n, NULL, &_ZM_sqr, &_ZM_mul));
     748             : }
     749             : /********************************************************************/
     750             : /**                                                                **/
     751             : /**                           ADD, SUB                             **/
     752             : /**                                                                **/
     753             : /********************************************************************/
     754             : static GEN
     755    34237944 : ZC_add_i(GEN x, GEN y, long lx)
     756             : {
     757    34237944 :   GEN A = cgetg(lx, t_COL);
     758             :   long i;
     759   472043560 :   for (i=1; i<lx; i++) gel(A,i) = addii(gel(x,i), gel(y,i));
     760    34224582 :   return A;
     761             : }
     762             : GEN
     763    24926027 : ZC_add(GEN x, GEN y) { return ZC_add_i(x, y, lg(x)); }
     764             : GEN
     765      361337 : ZC_Z_add(GEN x, GEN y)
     766             : {
     767      361337 :   long k, lx = lg(x);
     768      361337 :   GEN z = cgetg(lx, t_COL);
     769      361336 :   if (lx == 1) pari_err_TYPE2("+",x,y);
     770      361336 :   gel(z,1) = addii(y,gel(x,1));
     771     2479872 :   for (k = 2; k < lx; k++) gel(z,k) = icopy(gel(x,k));
     772      361355 :   return z;
     773             : }
     774             : 
     775             : static GEN
     776     8933827 : ZC_sub_i(GEN x, GEN y, long lx)
     777             : {
     778             :   long i;
     779     8933827 :   GEN A = cgetg(lx, t_COL);
     780    55977652 :   for (i=1; i<lx; i++) gel(A,i) = subii(gel(x,i), gel(y,i));
     781     8932629 :   return A;
     782             : }
     783             : GEN
     784     8642823 : ZC_sub(GEN x, GEN y) { return ZC_sub_i(x, y, lg(x)); }
     785             : GEN
     786           0 : ZC_Z_sub(GEN x, GEN y)
     787             : {
     788           0 :   long k, lx = lg(x);
     789           0 :   GEN z = cgetg(lx, t_COL);
     790           0 :   if (lx == 1) pari_err_TYPE2("+",x,y);
     791           0 :   gel(z,1) = subii(gel(x,1), y);
     792           0 :   for (k = 2; k < lx; k++) gel(z,k) = icopy(gel(x,k));
     793           0 :   return z;
     794             : }
     795             : GEN
     796      444912 : Z_ZC_sub(GEN a, GEN x)
     797             : {
     798      444912 :   long k, lx = lg(x);
     799      444912 :   GEN z = cgetg(lx, t_COL);
     800      444929 :   if (lx == 1) pari_err_TYPE2("-",a,x);
     801      444929 :   gel(z,1) = subii(a, gel(x,1));
     802     1212737 :   for (k = 2; k < lx; k++) gel(z,k) = negi(gel(x,k));
     803      444952 :   return z;
     804             : }
     805             : 
     806             : GEN
     807      818087 : ZM_add(GEN x, GEN y)
     808             : {
     809      818087 :   long lx = lg(x), l, j;
     810             :   GEN z;
     811      818087 :   if (lx == 1) return cgetg(1, t_MAT);
     812      738819 :   z = cgetg(lx, t_MAT); l = lgcols(x);
     813    10049880 :   for (j = 1; j < lx; j++) gel(z,j) = ZC_add_i(gel(x,j), gel(y,j), l);
     814      738817 :   return z;
     815             : }
     816             : GEN
     817      113046 : ZM_sub(GEN x, GEN y)
     818             : {
     819      113046 :   long lx = lg(x), l, j;
     820             :   GEN z;
     821      113046 :   if (lx == 1) return cgetg(1, t_MAT);
     822      113046 :   z = cgetg(lx, t_MAT); l = lgcols(x);
     823      404020 :   for (j = 1; j < lx; j++) gel(z,j) = ZC_sub_i(gel(x,j), gel(y,j), l);
     824      113045 :   return z;
     825             : }
     826             : /********************************************************************/
     827             : /**                                                                **/
     828             : /**                         LINEAR COMBINATION                     **/
     829             : /**                                                                **/
     830             : /********************************************************************/
     831             : /* return X/c assuming division is exact */
     832             : GEN
     833     9212682 : ZC_Z_divexact(GEN x, GEN c)
     834    77907399 : { pari_APPLY_type(t_COL, diviiexact(gel(x,i), c)) }
     835             : GEN
     836        2240 : ZC_divexactu(GEN x, ulong c)
     837       11354 : { pari_APPLY_type(t_COL, diviuexact(gel(x,i), c)) }
     838             : 
     839             : GEN
     840     2067930 : ZM_Z_divexact(GEN x, GEN c)
     841     9223537 : { pari_APPLY_same(ZC_Z_divexact(gel(x,i), c)) }
     842             : 
     843             : GEN
     844         434 : ZM_divexactu(GEN x, ulong c)
     845        2674 : { pari_APPLY_same(ZC_divexactu(gel(x,i), c)) }
     846             : 
     847             : GEN
     848    33154284 : ZC_Z_mul(GEN x, GEN c)
     849             : {
     850    33154284 :   if (!signe(c)) return zerocol(lg(x)-1);
     851    31909875 :   if (is_pm1(c)) return (signe(c) > 0)? ZC_copy(x): ZC_neg(x);
     852   249913564 :   pari_APPLY_type(t_COL, mulii(gel(x,i), c))
     853             : }
     854             : 
     855             : GEN
     856       47563 : ZC_z_mul(GEN x, long c)
     857             : {
     858       47563 :   if (!c) return zerocol(lg(x)-1);
     859       40013 :   if (c == 1) return ZC_copy(x);
     860       35453 :   if (c ==-1) return ZC_neg(x);
     861      441132 :   pari_APPLY_type(t_COL, mulsi(c, gel(x,i)))
     862             : }
     863             : 
     864             : GEN
     865       28529 : zv_z_mul(GEN x, long n)
     866      431552 : { pari_APPLY_long(x[i]*n) }
     867             : 
     868             : /* return a ZM */
     869             : GEN
     870           0 : nm_Z_mul(GEN X, GEN c)
     871             : {
     872           0 :   long i, j, h, l = lg(X), s = signe(c);
     873             :   GEN A;
     874           0 :   if (l == 1) return cgetg(1, t_MAT);
     875           0 :   h = lgcols(X);
     876           0 :   if (!s) return zeromat(h-1, l-1);
     877           0 :   if (is_pm1(c)) {
     878           0 :     if (s > 0) return Flm_to_ZM(X);
     879           0 :     X = Flm_to_ZM(X); ZM_togglesign(X); return X;
     880             :   }
     881           0 :   A = cgetg(l, t_MAT);
     882           0 :   for (j = 1; j < l; j++)
     883             :   {
     884           0 :     GEN a = cgetg(h, t_COL), x = gel(X, j);
     885           0 :     for (i = 1; i < h; i++) gel(a,i) = muliu(c, x[i]);
     886           0 :     gel(A,j) = a;
     887             :   }
     888           0 :   return A;
     889             : }
     890             : GEN
     891     3248172 : ZM_Z_mul(GEN X, GEN c)
     892             : {
     893     3248172 :   long i, j, h, l = lg(X);
     894             :   GEN A;
     895     3248172 :   if (l == 1) return cgetg(1, t_MAT);
     896     3248172 :   h = lgcols(X);
     897     3248168 :   if (!signe(c)) return zeromat(h-1, l-1);
     898     3203103 :   if (is_pm1(c)) return (signe(c) > 0)? ZM_copy(X): ZM_neg(X);
     899     2919517 :   A = cgetg(l, t_MAT);
     900    18084035 :   for (j = 1; j < l; j++)
     901             :   {
     902    15165754 :     GEN a = cgetg(h, t_COL), x = gel(X, j);
     903   229517801 :     for (i = 1; i < h; i++) gel(a,i) = mulii(c, gel(x,i));
     904    15164512 :     gel(A,j) = a;
     905             :   }
     906     2918281 :   return A;
     907             : }
     908             : void
     909    90017350 : ZC_lincomb1_inplace_i(GEN X, GEN Y, GEN v, long n)
     910             : {
     911             :   long i;
     912  2038929935 :   for (i = n; i; i--) gel(X,i) = addmulii_inplace(gel(X,i), gel(Y,i), v);
     913    89220946 : }
     914             : /* X <- X + v Y (elementary col operation) */
     915             : void
     916    79542626 : ZC_lincomb1_inplace(GEN X, GEN Y, GEN v)
     917             : {
     918    79542626 :   if (lgefint(v) != 2) return ZC_lincomb1_inplace_i(X, Y, v, lg(X)-1);
     919             : }
     920             : void
     921    29271253 : Flc_lincomb1_inplace(GEN X, GEN Y, ulong v, ulong q)
     922             : {
     923             :   long i;
     924    29271253 :   if (!v) return; /* v = 0 */
     925   745078356 :   for (i = lg(X)-1; i; i--) X[i] = Fl_add(X[i], Fl_mul(Y[i], v, q), q);
     926             : }
     927             : 
     928             : /* X + v Y, wasteful if (v = 0) */
     929             : static GEN
     930    17203972 : ZC_lincomb1(GEN v, GEN x, GEN y)
     931   127234601 : { pari_APPLY_type(t_COL, addmulii(gel(x,i), gel(y,i), v)) }
     932             : 
     933             : /* -X + vY */
     934             : static GEN
     935      603214 : ZC_lincomb_1(GEN v, GEN x, GEN y)
     936     4079253 : { pari_APPLY_type(t_COL, mulsubii(gel(y,i), v, gel(x,i))) }
     937             : 
     938             : /* X,Y compatible ZV; u,v in Z. Returns A = u*X + v*Y */
     939             : GEN
     940    33904231 : ZC_lincomb(GEN u, GEN v, GEN X, GEN Y)
     941             : {
     942             :   long su, sv;
     943             :   GEN A;
     944             : 
     945    33904231 :   su = signe(u); if (!su) return ZC_Z_mul(Y, v);
     946    33903273 :   sv = signe(v); if (!sv) return ZC_Z_mul(X, u);
     947    33476201 :   if (is_pm1(v))
     948             :   {
     949    10547810 :     if (is_pm1(u))
     950             :     {
     951     9379381 :       if (su != sv) A = ZC_sub(X, Y);
     952     2479334 :       else          A = ZC_add(X, Y);
     953     9378794 :       if (su < 0) ZV_togglesign(A); /* in place but was created above */
     954             :     }
     955             :     else
     956             :     {
     957     1168400 :       if (sv > 0) A = ZC_lincomb1 (u, Y, X);
     958      471698 :       else        A = ZC_lincomb_1(u, Y, X);
     959             :     }
     960             :   }
     961    22929702 :   else if (is_pm1(u))
     962             :   {
     963    16638338 :     if (su > 0) A = ZC_lincomb1 (v, X, Y);
     964      130993 :     else        A = ZC_lincomb_1(v, X, Y);
     965             :   }
     966             :   else
     967             :   { /* not cgetg_copy: x may be a t_VEC */
     968     6294943 :     long i, lx = lg(X);
     969     6294943 :     A = cgetg(lx,t_COL);
     970    40420522 :     for (i=1; i<lx; i++) gel(A,i) = lincombii(u,v,gel(X,i),gel(Y,i));
     971             :   }
     972    33469010 :   return A;
     973             : }
     974             : 
     975             : /********************************************************************/
     976             : /**                                                                **/
     977             : /**                           CONVERSIONS                          **/
     978             : /**                                                                **/
     979             : /********************************************************************/
     980             : GEN
     981      393789 : ZV_to_nv(GEN x)
     982      739692 : { pari_APPLY_ulong(itou(gel(x,i))) }
     983             : 
     984             : GEN
     985      133485 : zm_to_ZM(GEN x)
     986      744964 : { pari_APPLY_type(t_MAT, zc_to_ZC(gel(x,i))) }
     987             : 
     988             : GEN
     989         126 : zmV_to_ZMV(GEN x)
     990         791 : { pari_APPLY_type(t_VEC, zm_to_ZM(gel(x,i))) }
     991             : 
     992             : /* same as Flm_to_ZM but do not assume positivity */
     993             : GEN
     994        1022 : ZM_to_zm(GEN x)
     995       17472 : { pari_APPLY_same(ZV_to_zv(gel(x,i))) }
     996             : 
     997             : GEN
     998      366646 : zv_to_Flv(GEN x, ulong p)
     999     5418812 : { pari_APPLY_ulong(umodsu(x[i], p)) }
    1000             : 
    1001             : GEN
    1002       22694 : zm_to_Flm(GEN x, ulong p)
    1003      351750 : { pari_APPLY_same(zv_to_Flv(gel(x,i),p)) }
    1004             : 
    1005             : GEN
    1006          49 : ZMV_to_zmV(GEN x)
    1007         399 : { pari_APPLY_type(t_VEC, ZM_to_zm(gel(x,i))) }
    1008             : 
    1009             : /********************************************************************/
    1010             : /**                                                                **/
    1011             : /**                         COPY, NEGATION                         **/
    1012             : /**                                                                **/
    1013             : /********************************************************************/
    1014             : GEN
    1015    12615181 : ZC_copy(GEN x)
    1016             : {
    1017    12615181 :   long i, lx = lg(x);
    1018    12615181 :   GEN y = cgetg(lx, t_COL);
    1019    81807744 :   for (i=1; i<lx; i++)
    1020             :   {
    1021    69192802 :     GEN c = gel(x,i);
    1022    69192802 :     gel(y,i) = lgefint(c) == 2? gen_0: icopy(c);
    1023             :   }
    1024    12614942 :   return y;
    1025             : }
    1026             : 
    1027             : GEN
    1028      512716 : ZM_copy(GEN x)
    1029     4308892 : { pari_APPLY_same(ZC_copy(gel(x,i))) }
    1030             : 
    1031             : void
    1032      340178 : ZV_neg_inplace(GEN M)
    1033             : {
    1034      340178 :   long l = lg(M);
    1035     1252026 :   while (--l > 0) gel(M,l) = negi(gel(M,l));
    1036      340200 : }
    1037             : GEN
    1038     3523854 : ZC_neg(GEN x)
    1039    27107321 : { pari_APPLY_type(t_COL, negi(gel(x,i))) }
    1040             : 
    1041             : GEN
    1042       51336 : zv_neg(GEN x)
    1043      658452 : { pari_APPLY_long(-x[i]) }
    1044             : GEN
    1045         126 : zv_neg_inplace(GEN M)
    1046             : {
    1047         126 :   long l = lg(M);
    1048         534 :   while (--l > 0) M[l] = -M[l];
    1049         126 :   return M;
    1050             : }
    1051             : GEN
    1052          77 : zv_abs(GEN x)
    1053        5446 : { pari_APPLY_ulong(labs(x[i])) }
    1054             : GEN
    1055      319744 : ZM_neg(GEN x)
    1056     1310562 : { pari_APPLY_same(ZC_neg(gel(x,i))) }
    1057             : 
    1058             : void
    1059     4988902 : ZV_togglesign(GEN M)
    1060             : {
    1061     4988902 :   long l = lg(M);
    1062    74560597 :   while (--l > 0) togglesign_safe(&gel(M,l));
    1063     4989025 : }
    1064             : void
    1065           0 : ZM_togglesign(GEN M)
    1066             : {
    1067           0 :   long l = lg(M);
    1068           0 :   while (--l > 0) ZV_togglesign(gel(M,l));
    1069           0 : }
    1070             : 
    1071             : /********************************************************************/
    1072             : /**                                                                **/
    1073             : /**                        "DIVISION" mod HNF                      **/
    1074             : /**                                                                **/
    1075             : /********************************************************************/
    1076             : /* Reduce ZC x modulo ZM y in HNF, may return x itself (not a copy) */
    1077             : GEN
    1078    11020721 : ZC_hnfremdiv(GEN x, GEN y, GEN *Q)
    1079             : {
    1080    11020721 :   long i, l = lg(x);
    1081             :   GEN q;
    1082             : 
    1083    11020721 :   if (Q) *Q = cgetg(l,t_COL);
    1084    11020721 :   if (l == 1) return cgetg(1,t_COL);
    1085    58465391 :   for (i = l-1; i>0; i--)
    1086             :   {
    1087    47447588 :     q = diviiround(gel(x,i), gcoeff(y,i,i));
    1088    47446117 :     if (signe(q)) {
    1089    24488096 :       togglesign(q);
    1090    24488329 :       x = ZC_lincomb(gen_1, q, x, gel(y,i));
    1091             :     }
    1092    47444670 :     if (Q) gel(*Q, i) = q;
    1093             :   }
    1094    11017803 :   return x;
    1095             : }
    1096             : 
    1097             : /* x = y Q + R, may return some columns of x (not copies) */
    1098             : GEN
    1099      453542 : ZM_hnfdivrem(GEN x, GEN y, GEN *Q)
    1100             : {
    1101      453542 :   long lx = lg(x), i;
    1102      453542 :   GEN R = cgetg(lx, t_MAT);
    1103      453550 :   if (Q)
    1104             :   {
    1105      127200 :     GEN q = cgetg(lx, t_MAT); *Q = q;
    1106      187983 :     for (i=1; i<lx; i++) gel(R,i) = ZC_hnfremdiv(gel(x,i),y,(GEN*)(q+i));
    1107             :   }
    1108             :   else
    1109      822652 :     for (i=1; i<lx; i++)
    1110             :     {
    1111      496299 :       pari_sp av = avma;
    1112      496299 :       GEN z = ZC_hnfrem(gel(x,i),y);
    1113      496283 :       gel(R,i) = (avma == av)? ZC_copy(z): gerepileupto(av, z);
    1114             :     }
    1115      453555 :   return R;
    1116             : }
    1117             : 
    1118             : /********************************************************************/
    1119             : /**                                                                **/
    1120             : /**                               TESTS                            **/
    1121             : /**                                                                **/
    1122             : /********************************************************************/
    1123             : int
    1124    23099474 : zv_equal0(GEN V)
    1125             : {
    1126    23099474 :   long l = lg(V);
    1127    37561040 :   while (--l > 0)
    1128    30787797 :     if (V[l]) return 0;
    1129     6773243 :   return 1;
    1130             : }
    1131             : 
    1132             : int
    1133     1693208 : ZV_equal0(GEN V)
    1134             : {
    1135     1693208 :   long l = lg(V);
    1136     4597394 :   while (--l > 0)
    1137     4394224 :     if (signe(gel(V,l))) return 0;
    1138      203170 :   return 1;
    1139             : }
    1140             : int
    1141       15974 : ZMrow_equal0(GEN V, long i)
    1142             : {
    1143       15974 :   long l = lg(V);
    1144       24891 :   while (--l > 0)
    1145       21429 :     if (signe(gcoeff(V,i,l))) return 0;
    1146        3462 :   return 1;
    1147             : }
    1148             : 
    1149             : static int
    1150     4462466 : ZV_equal_lg(GEN V, GEN W, long l)
    1151             : {
    1152    19239418 :   while (--l > 0)
    1153    15179813 :     if (!equalii(gel(V,l), gel(W,l))) return 0;
    1154     4059605 :   return 1;
    1155             : }
    1156             : int
    1157      238358 : ZV_equal(GEN V, GEN W)
    1158             : {
    1159      238358 :   long l = lg(V);
    1160      238358 :   if (lg(W) != l) return 0;
    1161      238358 :   return ZV_equal_lg(V, W, l);
    1162             : }
    1163             : int
    1164     1767985 : ZM_equal(GEN A, GEN B)
    1165             : {
    1166     1767985 :   long i, m, l = lg(A);
    1167     1767985 :   if (lg(B) != l) return 0;
    1168     1767627 :   if (l == 1) return 1;
    1169     1767627 :   m = lgcols(A);
    1170     1767626 :   if (lgcols(B) != m) return 0;
    1171     5766108 :   for (i = 1; i < l; i++)
    1172     4224109 :     if (!ZV_equal_lg(gel(A,i), gel(B,i), m)) return 0;
    1173     1541999 :   return 1;
    1174             : }
    1175             : int
    1176       70217 : ZM_equal0(GEN A)
    1177             : {
    1178       70217 :   long i, j, m, l = lg(A);
    1179       70217 :   if (l == 1) return 1;
    1180       70217 :   m = lgcols(A);
    1181      121019 :   for (j = 1; j < l; j++)
    1182     2616058 :     for (i = 1; i < m; i++)
    1183     2565256 :       if (signe(gcoeff(A,i,j))) return 0;
    1184       15541 :   return 1;
    1185             : }
    1186             : int
    1187    30764946 : zv_equal(GEN V, GEN W)
    1188             : {
    1189    30764946 :   long l = lg(V);
    1190    30764946 :   if (lg(W) != l) return 0;
    1191   261905846 :   while (--l > 0)
    1192   232214731 :     if (V[l] != W[l]) return 0;
    1193    29691115 :   return 1;
    1194             : }
    1195             : 
    1196             : int
    1197        1638 : zvV_equal(GEN V, GEN W)
    1198             : {
    1199        1638 :   long l = lg(V);
    1200        1638 :   if (lg(W) != l) return 0;
    1201       80388 :   while (--l > 0)
    1202       79912 :     if (!zv_equal(gel(V,l),gel(W,l))) return 0;
    1203         476 :   return 1;
    1204             : }
    1205             : 
    1206             : int
    1207      748494 : ZM_ishnf(GEN x)
    1208             : {
    1209      748494 :   long i,j, lx = lg(x);
    1210     2266804 :   for (i=1; i<lx; i++)
    1211             :   {
    1212     1621538 :     GEN xii = gcoeff(x,i,i);
    1213     1621538 :     if (signe(xii) <= 0) return 0;
    1214     3165460 :     for (j=1; j<i; j++)
    1215     1574162 :       if (signe(gcoeff(x,i,j))) return 0;
    1216     3221145 :     for (j=i+1; j<lx; j++)
    1217             :     {
    1218     1702835 :       GEN xij = gcoeff(x,i,j);
    1219     1702835 :       if (signe(xij)<0 || cmpii(xij,xii)>=0) return 0;
    1220             :     }
    1221             :   }
    1222      645266 :   return 1;
    1223             : }
    1224             : int
    1225      477833 : ZM_isidentity(GEN x)
    1226             : {
    1227      477833 :   long i,j, lx = lg(x);
    1228             : 
    1229      477833 :   if (lx == 1) return 1;
    1230      477826 :   if (lx != lgcols(x)) return 0;
    1231     2348533 :   for (j=1; j<lx; j++)
    1232             :   {
    1233     1870763 :     GEN c = gel(x,j);
    1234     5924740 :     for (i=1; i<j; )
    1235     4053998 :       if (signe(gel(c,i++))) return 0;
    1236             :     /* i = j */
    1237     1870742 :     if (!equali1(gel(c,i++))) return 0;
    1238     5924754 :     for (   ; i<lx; )
    1239     4054033 :       if (signe(gel(c,i++))) return 0;
    1240             :   }
    1241      477770 :   return 1;
    1242             : }
    1243             : int
    1244      555397 : ZM_isdiagonal(GEN x)
    1245             : {
    1246      555397 :   long i,j, lx = lg(x);
    1247      555397 :   if (lx == 1) return 1;
    1248      555397 :   if (lx != lgcols(x)) return 0;
    1249             : 
    1250     1434074 :   for (j=1; j<lx; j++)
    1251             :   {
    1252     1205032 :     GEN c = gel(x,j);
    1253     1687700 :     for (i=1; i<j; i++)
    1254      809006 :       if (signe(gel(c,i))) return 0;
    1255     2007392 :     for (i++; i<lx; i++)
    1256     1128733 :       if (signe(gel(c,i))) return 0;
    1257             :   }
    1258      229042 :   return 1;
    1259             : }
    1260             : int
    1261      106678 : ZM_isscalar(GEN x, GEN s)
    1262             : {
    1263      106678 :   long i, j, lx = lg(x);
    1264             : 
    1265      106678 :   if (lx == 1) return 1;
    1266      106678 :   if (!s) s = gcoeff(x,1,1);
    1267      106678 :   if (equali1(s)) return ZM_isidentity(x);
    1268      105425 :   if (lx != lgcols(x)) return 0;
    1269      555866 :   for (j=1; j<lx; j++)
    1270             :   {
    1271      467706 :     GEN c = gel(x,j);
    1272     2384584 :     for (i=1; i<j; )
    1273     1931924 :       if (signe(gel(c,i++))) return 0;
    1274             :     /* i = j */
    1275      452660 :     if (!equalii(gel(c,i++), s)) return 0;
    1276     2428013 :     for (   ; i<lx; )
    1277     1977572 :       if (signe(gel(c,i++))) return 0;
    1278             :   }
    1279       88160 :   return 1;
    1280             : }
    1281             : 
    1282             : long
    1283      106677 : ZC_is_ei(GEN x)
    1284             : {
    1285      106677 :   long i, j = 0, l = lg(x);
    1286      995426 :   for (i = 1; i < l; i++)
    1287             :   {
    1288      888750 :     GEN c = gel(x,i);
    1289      888750 :     long s = signe(c);
    1290      888750 :     if (!s) continue;
    1291      106671 :     if (s < 0 || !is_pm1(c) || j) return 0;
    1292      106670 :     j = i;
    1293             :   }
    1294      106676 :   return j;
    1295             : }
    1296             : 
    1297             : /********************************************************************/
    1298             : /**                                                                **/
    1299             : /**                       MISCELLANEOUS                            **/
    1300             : /**                                                                **/
    1301             : /********************************************************************/
    1302             : /* assume lg(x) = lg(y), x,y in Z^n */
    1303             : int
    1304     2963245 : ZV_cmp(GEN x, GEN y)
    1305             : {
    1306     2963245 :   long fl,i, lx = lg(x);
    1307     5694906 :   for (i=1; i<lx; i++)
    1308     4575009 :     if (( fl = cmpii(gel(x,i), gel(y,i)) )) return fl;
    1309     1119897 :   return 0;
    1310             : }
    1311             : /* assume lg(x) = lg(y), x,y in Z^n */
    1312             : int
    1313       21743 : ZV_abscmp(GEN x, GEN y)
    1314             : {
    1315       21743 :   long fl,i, lx = lg(x);
    1316       68161 :   for (i=1; i<lx; i++)
    1317       68039 :     if (( fl = abscmpii(gel(x,i), gel(y,i)) )) return fl;
    1318         122 :   return 0;
    1319             : }
    1320             : 
    1321             : long
    1322    23092930 : zv_content(GEN x)
    1323             : {
    1324    23092930 :   long i, s, l = lg(x);
    1325    23092930 :   if (l == 1) return 0;
    1326    23092923 :   s = labs(x[1]);
    1327    51268321 :   for (i = 2; i < l && s != 1; i++) s = ugcd(s, labs(x[i]));
    1328    23093814 :   return s;
    1329             : }
    1330             : GEN
    1331      406299 : ZV_content(GEN x)
    1332             : {
    1333      406299 :   long i, l = lg(x);
    1334      406299 :   pari_sp av = avma;
    1335             :   GEN c;
    1336      406299 :   if (l == 1) return gen_0;
    1337      406299 :   if (l == 2) return absi(gel(x,1));
    1338      332057 :   c = gel(x,1);
    1339      895281 :   for (i = 2; i < l; i++)
    1340             :   {
    1341      665957 :     c = gcdii(c, gel(x,i));
    1342      665952 :     if (is_pm1(c)) { set_avma(av); return gen_1; }
    1343             :   }
    1344      229324 :   return gerepileuptoint(av, c);
    1345             : }
    1346             : 
    1347             : GEN
    1348     3861950 : ZM_det_triangular(GEN mat)
    1349             : {
    1350             :   pari_sp av;
    1351     3861950 :   long i,l = lg(mat);
    1352             :   GEN s;
    1353             : 
    1354     3861950 :   if (l<3) return l<2? gen_1: icopy(gcoeff(mat,1,1));
    1355     3462076 :   av = avma; s = gcoeff(mat,1,1);
    1356     9389123 :   for (i=2; i<l; i++) s = mulii(s,gcoeff(mat,i,i));
    1357     3461614 :   return gerepileuptoint(av,s);
    1358             : }
    1359             : 
    1360             : /* assumes no overflow */
    1361             : long
    1362      939622 : zv_prod(GEN v)
    1363             : {
    1364      939622 :   long n, i, l = lg(v);
    1365      939622 :   if (l == 1) return 1;
    1366      951796 :   n = v[1]; for (i = 2; i < l; i++) n *= v[i];
    1367      763593 :   return n;
    1368             : }
    1369             : 
    1370             : static GEN
    1371   355912400 : _mulii(void *E, GEN a, GEN b)
    1372   355912400 : { (void) E; return mulii(a, b); }
    1373             : 
    1374             : /* product of ulongs */
    1375             : GEN
    1376     1861626 : zv_prod_Z(GEN v)
    1377             : {
    1378     1861626 :   pari_sp av = avma;
    1379     1861626 :   long k, m, n = lg(v)-1;
    1380             :   GEN V;
    1381     1861626 :   switch(n) {
    1382       20902 :     case 0: return gen_1;
    1383      121058 :     case 1: return utoi(v[1]);
    1384      978419 :     case 2: return muluu(v[1], v[2]);
    1385             :   }
    1386      741247 :   m = n >> 1;
    1387      741247 :   V = cgetg(m + (odd(n)? 2: 1), t_VEC);
    1388    84985339 :   for (k = 1; k <= m; k++) gel(V,k) = muluu(v[k<<1], v[(k<<1)-1]);
    1389      741180 :   if (odd(n)) gel(V,k) = utoipos(v[n]);
    1390      741241 :   return gerepileuptoint(av, gen_product(V, NULL, &_mulii));
    1391             : }
    1392             : GEN
    1393    14694393 : vecsmall_prod(GEN v)
    1394             : {
    1395    14694393 :   pari_sp av = avma;
    1396    14694393 :   long k, m, n = lg(v)-1;
    1397             :   GEN V;
    1398    14694393 :   switch (n) {
    1399           0 :     case 0: return gen_1;
    1400           0 :     case 1: return stoi(v[1]);
    1401          21 :     case 2: return mulss(v[1], v[2]);
    1402             :   }
    1403    14694372 :   m = n >> 1;
    1404    14694372 :   V = cgetg(m + (odd(n)? 2: 1), t_VEC);
    1405   161556906 :   for (k = 1; k <= m; k++) gel(V,k) = mulss(v[k<<1], v[(k<<1)-1]);
    1406    14694372 :   if (odd(n)) gel(V,k) = stoi(v[n]);
    1407    14694372 :   return gerepileuptoint(av, gen_product(V, NULL, &_mulii));
    1408             : }
    1409             : 
    1410             : GEN
    1411     8811185 : ZV_prod(GEN v)
    1412             : {
    1413     8811185 :   pari_sp av = avma;
    1414     8811185 :   long i, l = lg(v);
    1415             :   GEN n;
    1416     8811185 :   if (l == 1) return gen_1;
    1417     8629832 :   if (l > 7) return gerepileuptoint(av, gen_product(v, NULL, _mulii));
    1418     1287728 :   n = gel(v,1);
    1419     1287728 :   if (l == 2) return icopy(n);
    1420     2088155 :   for (i = 2; i < l; i++) n = mulii(n, gel(v,i));
    1421      839214 :   return gerepileuptoint(av, n);
    1422             : }
    1423             : /* assumes no overflow */
    1424             : long
    1425       16464 : zv_sum(GEN v)
    1426             : {
    1427       16464 :   long n, i, l = lg(v);
    1428       16464 :   if (l == 1) return 0;
    1429       95924 :   n = v[1]; for (i = 2; i < l; i++) n += v[i];
    1430       16443 :   return n;
    1431             : }
    1432             : /* assumes no overflow and 0 <= n <= #v */
    1433             : long
    1434           0 : zv_sumpart(GEN v, long n)
    1435             : {
    1436             :   long i, p;
    1437           0 :   if (!n) return 0;
    1438           0 :   p = v[1]; for (i = 2; i <= n; i++) p += v[i];
    1439           0 :   return p;
    1440             : }
    1441             : GEN
    1442          63 : ZV_sum(GEN v)
    1443             : {
    1444          63 :   pari_sp av = avma;
    1445          63 :   long i, l = lg(v);
    1446             :   GEN n;
    1447          63 :   if (l == 1) return gen_0;
    1448          63 :   n = gel(v,1);
    1449          63 :   if (l == 2) return icopy(n);
    1450         532 :   for (i = 2; i < l; i++) n = addii(n, gel(v,i));
    1451          63 :   return gerepileuptoint(av, n);
    1452             : }
    1453             : 
    1454             : /********************************************************************/
    1455             : /**                                                                **/
    1456             : /**         GRAM SCHMIDT REDUCTION (integer matrices)              **/
    1457             : /**                                                                **/
    1458             : /********************************************************************/
    1459             : 
    1460             : /* L[k,] += q * L[l,], l < k. Inefficient if q = 0 */
    1461             : static void
    1462      360456 : Zupdate_row(long k, long l, GEN q, GEN L, GEN B)
    1463             : {
    1464      360456 :   long i, qq = itos_or_0(q);
    1465      360457 :   if (!qq)
    1466             :   {
    1467       27955 :     for(i=1;i<l;i++)  gcoeff(L,k,i) = addii(gcoeff(L,k,i),mulii(q,gcoeff(L,l,i)));
    1468        6453 :     gcoeff(L,k,l) = addii(gcoeff(L,k,l), mulii(q,B));
    1469        6453 :     return;
    1470             :   }
    1471      354004 :   if (qq == 1) {
    1472      149995 :     for (i=1;i<l; i++) gcoeff(L,k,i) = addii(gcoeff(L,k,i),gcoeff(L,l,i));
    1473      101431 :     gcoeff(L,k,l) = addii(gcoeff(L,k,l), B);
    1474      252573 :   } else if (qq == -1) {
    1475      154396 :     for (i=1;i<l; i++) gcoeff(L,k,i) = subii(gcoeff(L,k,i),gcoeff(L,l,i));
    1476       89806 :     gcoeff(L,k,l) = addii(gcoeff(L,k,l), negi(B));
    1477             :   } else {
    1478      273721 :     for(i=1;i<l;i++) gcoeff(L,k,i) = addii(gcoeff(L,k,i),mulsi(qq,gcoeff(L,l,i)));
    1479      162784 :     gcoeff(L,k,l) = addii(gcoeff(L,k,l), mulsi(qq,B));
    1480             :   }
    1481             : }
    1482             : 
    1483             : /* update L[k,] */
    1484             : static void
    1485     1032821 : ZRED(long k, long l, GEN x, GEN L, GEN B)
    1486             : {
    1487     1032821 :   GEN q = truedivii(addii(B,shifti(gcoeff(L,k,l),1)), shifti(B,1));
    1488     1032860 :   if (!signe(q)) return;
    1489      360429 :   q = negi(q);
    1490      360452 :   Zupdate_row(k,l,q,L,B);
    1491      360443 :   gel(x,k) = ZC_lincomb(gen_1, q, gel(x,k), gel(x,l));
    1492             : }
    1493             : 
    1494             : /* Gram-Schmidt reduction, x a ZM */
    1495             : static void
    1496     1182910 : ZincrementalGS(GEN x, GEN L, GEN B, long k)
    1497             : {
    1498             :   long i, j;
    1499     3777913 :   for (j=1; j<=k; j++)
    1500             :   {
    1501     2595610 :     pari_sp av = avma;
    1502     2595610 :     GEN u = ZV_dotproduct(gel(x,k), gel(x,j));
    1503     5603748 :     for (i=1; i<j; i++)
    1504             :     {
    1505     3009416 :       u = subii(mulii(gel(B,i+1), u), mulii(gcoeff(L,k,i), gcoeff(L,j,i)));
    1506     3008853 :       u = diviiexact(u, gel(B,i));
    1507             :     }
    1508     2594332 :     gcoeff(L,k,j) = gerepileuptoint(av, u);
    1509             :   }
    1510     1182303 :   gel(B,k+1) = gcoeff(L,k,k); gcoeff(L,k,k) = gen_1;
    1511     1182303 : }
    1512             : 
    1513             : /* Variant reducemodinvertible(ZC v, ZM y), when y singular.
    1514             :  * Very inefficient if y is not LLL-reduced of maximal rank */
    1515             : static GEN
    1516      110114 : ZC_reducemodmatrix_i(GEN v, GEN y)
    1517             : {
    1518      110114 :   GEN B, L, x = shallowconcat(y, v);
    1519      110114 :   long k, lx = lg(x), nx = lx-1;
    1520             : 
    1521      110114 :   B = scalarcol_shallow(gen_1, lx);
    1522      110115 :   L = zeromatcopy(nx, nx);
    1523      453928 :   for (k=1; k <= nx; k++) ZincrementalGS(x, L, B, k);
    1524      343812 :   for (k = nx-1; k >= 1; k--) ZRED(nx,k, x,L,gel(B,k+1));
    1525      110101 :   return gel(x,nx);
    1526             : }
    1527             : GEN
    1528      110114 : ZC_reducemodmatrix(GEN v, GEN y) {
    1529      110114 :   pari_sp av = avma;
    1530      110114 :   return gerepilecopy(av, ZC_reducemodmatrix_i(v,y));
    1531             : }
    1532             : 
    1533             : /* Variant reducemodinvertible(ZM v, ZM y), when y singular.
    1534             :  * Very inefficient if y is not LLL-reduced of maximal rank */
    1535             : static GEN
    1536      226674 : ZM_reducemodmatrix_i(GEN v, GEN y)
    1537             : {
    1538             :   GEN B, L, V;
    1539      226674 :   long j, k, lv = lg(v), nx = lg(y), lx = nx+1;
    1540             : 
    1541      226674 :   V = cgetg(lv, t_MAT);
    1542      226676 :   B = scalarcol_shallow(gen_1, lx);
    1543      226681 :   L = zeromatcopy(nx, nx);
    1544      600467 :   for (k=1; k < nx; k++) ZincrementalGS(y, L, B, k);
    1545      692001 :   for (j = 1; j < lg(v); j++)
    1546             :   {
    1547      465347 :     GEN x = shallowconcat(y, gel(v,j));
    1548      465370 :     ZincrementalGS(x, L, B, nx); /* overwrite last */
    1549     1264494 :     for (k = nx-1; k >= 1; k--) ZRED(nx,k, x,L,gel(B,k+1));
    1550      465332 :     gel(V,j) = gel(x,nx);
    1551             :   }
    1552      226654 :   return V;
    1553             : }
    1554             : GEN
    1555      226671 : ZM_reducemodmatrix(GEN v, GEN y) {
    1556      226671 :   pari_sp av = avma;
    1557      226671 :   return gerepilecopy(av, ZM_reducemodmatrix_i(v,y));
    1558             : }
    1559             : 
    1560             : GEN
    1561       98861 : ZC_reducemodlll(GEN x,GEN y)
    1562             : {
    1563       98861 :   pari_sp av = avma;
    1564       98861 :   GEN z = ZC_reducemodmatrix(x, ZM_lll(y, 0.75, LLL_INPLACE));
    1565       98862 :   return gerepilecopy(av, z);
    1566             : }
    1567             : GEN
    1568           0 : ZM_reducemodlll(GEN x,GEN y)
    1569             : {
    1570           0 :   pari_sp av = avma;
    1571           0 :   GEN z = ZM_reducemodmatrix(x, ZM_lll(y, 0.75, LLL_INPLACE));
    1572           0 :   return gerepilecopy(av, z);
    1573             : }

Generated by: LCOV version 1.14