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 29115-f22e516b23) Lines: 779 888 87.7 %
Date: 2024-03-18 08:03:28 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     1713904 : check_ZV(GEN x, long l)
      20             : {
      21             :   long i;
      22    10117633 :   for (i=1; i<l; i++)
      23     8403820 :     if (typ(gel(x,i)) != t_INT) return 0;
      24     1713813 :   return 1;
      25             : }
      26             : void
      27     1421286 : RgV_check_ZV(GEN A, const char *s)
      28             : {
      29     1421286 :   if (!RgV_is_ZV(A)) pari_err_TYPE(stack_strcat(s," [integer vector]"), A);
      30     1421277 : }
      31             : void
      32      599298 : RgM_check_ZM(GEN A, const char *s)
      33             : {
      34      599298 :   long n = lg(A);
      35      599298 :   if (n != 1)
      36             :   {
      37      599151 :     long j, m = lgcols(A);
      38     2312964 :     for (j=1; j<n; j++)
      39     1713904 :       if (!check_ZV(gel(A,j), m))
      40          91 :         pari_err_TYPE(stack_strcat(s," [integer matrix]"), A);
      41             :   }
      42      599207 : }
      43             : 
      44             : static long
      45   110283277 : ZV_max_lg_i(GEN x, long m)
      46             : {
      47   110283277 :   long i, l = 2;
      48   873075749 :   for (i = 1; i < m; i++) l = maxss(l, lgefint(gel(x,i)));
      49   110284537 :   return l;
      50             : }
      51             : long
      52       10619 : ZV_max_lg(GEN x) { return ZV_max_lg_i(x, lg(x)); }
      53             : 
      54             : static long
      55    28873375 : ZM_max_lg_i(GEN x, long n, long m)
      56             : {
      57    28873375 :   long j, l = 2;
      58   139145697 :   for (j = 1; j < n; j++) l = maxss(l, ZV_max_lg_i(gel(x,j), m));
      59    28874044 :   return l;
      60             : }
      61             : long
      62       21598 : ZM_max_lg(GEN x)
      63             : {
      64       21598 :   long n = lg(x);
      65       21598 :   if (n == 1) return 2;
      66       21598 :   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        3791 : ZM_supnorm(GEN x)
      96             : {
      97        3791 :   long i, j, h, lx = lg(x);
      98        3791 :   GEN s = gen_0;
      99        3791 :   if (lx == 1) return gen_1;
     100        3791 :   h = lgcols(x);
     101       23246 :   for (j=1; j<lx; j++)
     102             :   {
     103       19455 :     GEN xj = gel(x,j);
     104      272254 :     for (i=1; i<h; i++)
     105             :     {
     106      252799 :       GEN c = gel(xj,i);
     107      252799 :       if (abscmpii(c, s) > 0) s = c;
     108             :     }
     109             :   }
     110        3791 :   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     2214660 : ZV_zc_mul(GEN x, GEN y)
     139             : {
     140     2214660 :   long j, l = lg(x);
     141     2214660 :   pari_sp av = avma;
     142     2214660 :   GEN s = mulis(gel(x,1),y[1]);
     143    50028650 :   for (j=2; j<l; j++)
     144    47813990 :     if (y[j]) s = addii(s, mulis(gel(x,j),y[j]));
     145     2214660 :   return gerepileuptoint(av,s);
     146             : }
     147             : 
     148             : /* x nonempty ZM, y a compatible zc (dimension > 0). */
     149             : static GEN
     150    20650902 : ZM_zc_mul_i(GEN x, GEN y, long c, long l)
     151             : {
     152             :   long i, j;
     153    20650902 :   GEN z = cgetg(l,t_COL);
     154             : 
     155   125857473 :   for (i=1; i<l; i++)
     156             :   {
     157   105212273 :     pari_sp av = avma;
     158   105212273 :     GEN s = mulis(gcoeff(x,i,1),y[1]);
     159  1171095412 :     for (j=2; j<c; j++)
     160  1065885057 :       if (y[j]) s = addii(s, mulis(gcoeff(x,i,j),y[j]));
     161   105210355 :     gel(z,i) = gerepileuptoint(av,s);
     162             :   }
     163    20645200 :   return z;
     164             : }
     165             : GEN
     166    18930079 : ZM_zc_mul(GEN x, GEN y) {
     167    18930079 :   long l = lg(x);
     168    18930079 :   if (l == 1) return cgetg(1, t_COL);
     169    18930079 :   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      810503 : ZM_zm_mul(GEN x, GEN y)
     193             : {
     194      810503 :   long j, c, l = lg(x), ly = lg(y);
     195      810503 :   GEN z = cgetg(ly, t_MAT);
     196      810504 :   if (l == 1) return z;
     197      810497 :   c = lgcols(x);
     198     2531359 :   for (j = 1; j < ly; j++) gel(z,j) = ZM_zc_mul_i(x, gel(y,j), l,c);
     199      810496 :   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      348552 : 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      348552 :   long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
     223      348552 :   GEN M = cgetg(n + 1, t_MAT), C;
     224             : 
     225     2762472 :   for (j = 1; j <= min_e; j++) {
     226     2413920 :     gel(M, j) = C = cgetg(m + 1, t_COL);
     227    40765346 :     for (i = 1; i <= min_d; i++)
     228    38351426 :       gel(C, i) = addii(gcoeff(A, ma + i, na + j),
     229    38351426 :                         gcoeff(B, mb + i, nb + j));
     230     2451124 :     for (; i <= da; i++) gel(C, i) = gcoeff(A, ma + i, na + j);
     231     2413920 :     for (; i <= db; i++) gel(C, i) = gcoeff(B, mb + i, nb + j);
     232     2413920 :     for (; i <= m; i++)  gel(C, i) = gen_0;
     233             :   }
     234      380664 :   for (; j <= ea; j++) {
     235       32112 :     gel(M, j) = C = cgetg(m + 1, t_COL);
     236      122111 :     for (i = 1; i <= da; i++) gel(C, i) = gcoeff(A, ma + i, na + j);
     237       32112 :     for (; i <= m; i++) gel(C, i) = gen_0;
     238             :   }
     239      348552 :   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      348552 :   for (; j <= n; j++) gel(M, j) = zerocol(m);
     245      348552 :   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      304983 : 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      304983 :   long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
     256      304983 :   GEN M = cgetg(n + 1, t_MAT), C;
     257             : 
     258     2408597 :   for (j = 1; j <= min_e; j++) {
     259     2103614 :     gel(M, j) = C = cgetg(m + 1, t_COL);
     260    36227909 :     for (i = 1; i <= min_d; i++)
     261    34124295 :       gel(C, i) = subii(gcoeff(A, ma + i, na + j),
     262    34124295 :                         gcoeff(B, mb + i, nb + j));
     263     2131363 :     for (; i <= da; i++) gel(C, i) = gcoeff(A, ma + i, na + j);
     264     2148823 :     for (; i <= db; i++) gel(C, i) = negi(gcoeff(B, mb + i, nb + j));
     265     2103614 :     for (; i <= m; i++) gel(C, i) = gen_0;
     266             :   }
     267      304983 :   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      326423 :   for (; j <= eb; j++) {
     273       21440 :     gel(M, j) = C = cgetg(m + 1, t_COL);
     274       69677 :     for (i = 1; i <= db; i++) gel(C, i) = negi(gcoeff(B, mb + i, nb + j));
     275       21440 :     for (; i <= m; i++) gel(C, i) = gen_0;
     276             :   }
     277      326423 :   for (; j <= n; j++) gel(M, j) = zerocol(m);
     278      304983 :   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       43569 : ZM_mul_sw(GEN A, GEN B, long m, long n, long p)
     286             : {
     287       43569 :   pari_sp av = avma;
     288       43569 :   long m1 = (m + 1)/2, m2 = m/2,
     289       43569 :     n1 = (n + 1)/2, n2 = n/2,
     290       43569 :     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       43569 :   T2 = subtract_slices(n1, p2, B, 0, n1, p1, p2, B, n1, n2, p1, p2);
     297       43569 :   S1 = subtract_slices(m2, n1, A, m1, m2, 0, n1, A, 0, m2, 0, n1);
     298       43569 :   M2 = ZM_mul_i(S1, T2, m2 + 1, n1 + 1, p2 + 1);
     299       43569 :   if (gc_needed(av, 1))
     300           0 :     gerepileall(av, 2, &T2, &M2);  /* destroy S1 */
     301       43569 :   T3 = subtract_slices(n1, p1, T2, 0, n1, 0, p2, B, 0, n1, 0, p1);
     302       43569 :   if (gc_needed(av, 1))
     303           0 :     gerepileall(av, 2, &M2, &T3);  /* destroy T2 */
     304       43569 :   S2 = add_slices(m2, n1, A, m1, m2, 0, n1, A, m1, m2, n1, n2);
     305       43569 :   T1 = subtract_slices(n1, p1, B, 0, n1, p1, p2, B, 0, n1, 0, p2);
     306       43569 :   M3 = ZM_mul_i(S2, T1, m2 + 1, n1 + 1, p2 + 1);
     307       43569 :   if (gc_needed(av, 1))
     308           0 :     gerepileall(av, 4, &M2, &T3, &S2, &M3);  /* destroy T1 */
     309       43569 :   S3 = subtract_slices(m1, n1, S2, 0, m2, 0, n1, A, 0, m1, 0, n1);
     310       43569 :   if (gc_needed(av, 1))
     311           0 :     gerepileall(av, 4, &M2, &T3, &M3, &S3);  /* destroy S2 */
     312       43569 :   A11 = matslice(A, 1, m1, 1, n1);
     313       43569 :   B11 = matslice(B, 1, n1, 1, p1);
     314       43569 :   M1 = ZM_mul_i(A11, B11, m1 + 1, n1 + 1, p1 + 1);
     315       43569 :   if (gc_needed(av, 1))
     316           0 :     gerepileall(av, 5, &M2, &T3, &M3, &S3, &M1);  /* destroy A11, B11 */
     317       43569 :   A12 = matslice(A, 1, m1, n1 + 1, n);
     318       43569 :   B21 = matslice(B, n1 + 1, n, 1, p1);
     319       43569 :   M4 = ZM_mul_i(A12, B21, m1 + 1, n2 + 1, p1 + 1);
     320       43569 :   if (gc_needed(av, 1))
     321           0 :     gerepileall(av, 6, &M2, &T3, &M3, &S3, &M1, &M4);  /* destroy A12, B21 */
     322       43569 :   C11 = add_slices(m1, p1, M1, 0, m1, 0, p1, M4, 0, m1, 0, p1);
     323       43569 :   if (gc_needed(av, 1))
     324           0 :     gerepileall(av, 6, &M2, &T3, &M3, &S3, &M1, &C11);  /* destroy M4 */
     325       43569 :   M5 = ZM_mul_i(S3, T3, m1 + 1, n1 + 1, p1 + 1);
     326       43569 :   S4 = subtract_slices(m1, n2, A, 0, m1, n1, n2, S3, 0, m1, 0, n2);
     327       43569 :   if (gc_needed(av, 1))
     328           5 :     gerepileall(av, 7, &M2, &T3, &M3, &M1, &C11, &M5, &S4);  /* destroy S3 */
     329       43569 :   T4 = add_slices(n2, p1, B, n1, n2, 0, p1, T3, 0, n2, 0, p1);
     330       43569 :   if (gc_needed(av, 1))
     331           0 :     gerepileall(av, 7, &M2, &M3, &M1, &C11, &M5, &S4, &T4);  /* destroy T3 */
     332       43569 :   V1 = subtract_slices(m1, p1, M1, 0, m1, 0, p1, M5, 0, m1, 0, p1);
     333       43569 :   if (gc_needed(av, 1))
     334           1 :     gerepileall(av, 6, &M2, &M3, &S4, &T4, &C11, &V1);  /* destroy M1, M5 */
     335       43569 :   B22 = matslice(B, n1 + 1, n, p1 + 1, p);
     336       43569 :   M6 = ZM_mul_i(S4, B22, m1 + 1, n2 + 1, p2 + 1);
     337       43569 :   if (gc_needed(av, 1))
     338           6 :     gerepileall(av, 6, &M2, &M3, &T4, &C11, &V1, &M6);  /* destroy S4, B22 */
     339       43569 :   A22 = matslice(A, m1 + 1, m, n1 + 1, n);
     340       43569 :   M7 = ZM_mul_i(A22, T4, m2 + 1, n2 + 1, p1 + 1);
     341       43569 :   if (gc_needed(av, 1))
     342           0 :     gerepileall(av, 6, &M2, &M3, &C11, &V1, &M6, &M7);  /* destroy A22, T4 */
     343       43569 :   V3 = add_slices(m1, p2, V1, 0, m1, 0, p2, M3, 0, m2, 0, p2);
     344       43569 :   C12 = add_slices(m1, p2, V3, 0, m1, 0, p2, M6, 0, m1, 0, p2);
     345       43569 :   if (gc_needed(av, 1))
     346          12 :     gerepileall(av, 6, &M2, &M3, &C11, &V1, &M7, &C12);  /* destroy V3, M6 */
     347       43569 :   V2 = add_slices(m2, p1, V1, 0, m2, 0, p1, M2, 0, m2, 0, p2);
     348       43569 :   if (gc_needed(av, 1))
     349           0 :     gerepileall(av, 5, &M3, &C11, &M7, &C12, &V2);  /* destroy V1, M2 */
     350       43569 :   C21 = add_slices(m2, p1, V2, 0, m2, 0, p1, M7, 0, m2, 0, p1);
     351       43569 :   if (gc_needed(av, 1))
     352           6 :     gerepileall(av, 5, &M3, &C11, &C12, &V2, &C21);  /* destroy M7 */
     353       43569 :   C22 = add_slices(m2, p2, V2, 0, m2, 0, p2, M3, 0, m2, 0, p2);
     354       43569 :   if (gc_needed(av, 1))
     355           0 :     gerepileall(av, 4, &C11, &C12, &C21, &C22);  /* destroy V2, M3 */
     356       43569 :   C = shallowconcat(vconcat(C11, C21), vconcat(C12, C22));
     357       43569 :   return gerepilecopy(av, C);
     358             : }
     359             : 
     360             : /* x[i,]*y. Assume lg(x) > 1 and 0 < i < lgcols(x) */
     361             : static GEN
     362   550839021 : ZMrow_ZC_mul_i(GEN x, GEN y, long i, long lx)
     363             : {
     364   550839021 :   pari_sp av = avma;
     365   550839021 :   GEN c = mulii(gcoeff(x,i,1), gel(y,1)), ZERO = gen_0;
     366             :   long k;
     367  5446173229 :   for (k = 2; k < lx; k++)
     368             :   {
     369  4895737787 :     GEN t = mulii(gcoeff(x,i,k), gel(y,k));
     370  4893977361 :     if (t != ZERO) c = addii(c, t);
     371             :   }
     372   550435442 :   return gerepileuptoint(av, c);
     373             : }
     374             : GEN
     375   154654179 : ZMrow_ZC_mul(GEN x, GEN y, long i)
     376   154654179 : { 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    75431002 : ZM_ZC_mul_i(GEN x, GEN y, long lx, long l)
     381             : {
     382    75431002 :   GEN z = cgetg(l,t_COL);
     383             :   long i;
     384   471664214 :   for (i=1; i<l; i++) gel(z,i) = ZMrow_ZC_mul_i(x,y,i,lx);
     385    75386436 :   return z;
     386             : }
     387             : 
     388             : static GEN
     389    14394408 : ZM_mul_classical(GEN x, GEN y, long l, long lx, long ly)
     390             : {
     391             :   long j;
     392    14394408 :   GEN z = cgetg(ly, t_MAT);
     393    66766912 :   for (j = 1; j < ly; j++)
     394    52376321 :     gel(z, j) = ZM_ZC_mul_i(x, gel(y, j), lx, l);
     395    14390591 :   return z;
     396             : }
     397             : 
     398             : static GEN
     399        1308 : ZM_mul_slice(GEN A, GEN B, GEN P, GEN *mod)
     400             : {
     401        1308 :   pari_sp av = avma;
     402        1308 :   long i, n = lg(P)-1;
     403             :   GEN H, T;
     404        1308 :   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        1308 :   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        1308 : ZM_mul_worker(GEN P, GEN A, GEN B)
     424             : {
     425        1308 :   GEN V = cgetg(3, t_VEC);
     426        1308 :   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    14437940 : sw_bound(long s)
     450    14437940 : { return s > 60 ? 2: s > 25 ? 4: s > 15 ? 8 : s > 8 ? 16 : 32; }
     451             : 
     452             : static GEN
     453    21210328 : ZM_mul_i(GEN x, GEN y, long l, long lx, long ly)
     454             : {
     455             :   long sx, sy, B;
     456    21210328 :   if (lx==3 && l==3 && ly==3) return ZM2_mul(x,y);
     457    14412847 :   sx = ZM_max_lg_i(x, lx, l);
     458    14413732 :   sy = ZM_max_lg_i(y, ly, lx);
     459    14413756 :   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    14412782 :   B = sw_bound(minss(sx, sy));
     463    14412782 :   if (l <= B || lx <= B || ly <= B)
     464    14369243 :     return ZM_mul_classical(x, y, l, lx, ly);
     465             :   else
     466       43539 :     return ZM_mul_sw(x, y, l - 1, lx - 1, ly - 1);
     467             : }
     468             : 
     469             : GEN
     470    21041555 : ZM_mul(GEN x, GEN y)
     471             : {
     472    21041555 :   long lx = lg(x), ly = lg(y);
     473    21041555 :   if (ly == 1) return cgetg(1,t_MAT);
     474    20906857 :   if (lx == 1) return zeromat(0, ly-1);
     475    20905275 :   return ZM_mul_i(x, y, lgcols(x), lx, ly);
     476             : }
     477             : 
     478             : GEN
     479      536873 : QM_mul(GEN x, GEN y)
     480             : {
     481      536873 :   GEN dx, nx = Q_primitive_part(x, &dx);
     482      536873 :   GEN dy, ny = Q_primitive_part(y, &dy);
     483      536873 :   GEN z = ZM_mul(nx, ny);
     484      536873 :   if (dx || dy)
     485             :   {
     486      453274 :     GEN d = dx ? dy ? gmul(dx, dy): dx : dy;
     487      453274 :     if (!gequal1(d)) z = ZM_Q_mul(z, d);
     488             :   }
     489      536873 :   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      132989 : QM_QC_mul(GEN x, GEN y)
     504             : {
     505      132989 :   GEN dx, nx = Q_primitive_part(x, &dx);
     506      132989 :   GEN dy, ny = Q_primitive_part(y, &dy);
     507      132989 :   GEN z = ZM_ZC_mul(nx, ny);
     508      132989 :   if (dx || dy)
     509             :   {
     510      132961 :     GEN d = dx ? dy ? gmul(dx, dy): dx : dy;
     511      132961 :     if (!gequal1(d)) z = ZC_Q_mul(z, d);
     512             :   }
     513      132989 :   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      530928 : ZM_diag_mul(GEN d, GEN m)
     554             : {
     555      530928 :   long i, j, l = lg(d), lm = lg(m);
     556      530928 :   GEN y = cgetg(lm, t_MAT);
     557     1514787 :   for (j=1; j<lm; j++) gel(y,j) = cgetg(l, t_COL);
     558     1641009 :   for (i=1; i<l; i++)
     559             :   {
     560     1110504 :     GEN c = gel(d,i);
     561     1110504 :     if (equali1(c))
     562      236584 :       for (j=1; j<lm; j++) gcoeff(y,i,j) = gcoeff(m,i,j);
     563             :     else
     564     3390136 :       for (j=1; j<lm; j++) gcoeff(y,i,j) = mulii(gcoeff(m,i,j), c);
     565             :   }
     566      530505 :   return y;
     567             : }
     568             : 
     569             : /* assume lx > 1 is lg(x) = lg(y) */
     570             : static GEN
     571    19062729 : ZV_dotproduct_i(GEN x, GEN y, long lx)
     572             : {
     573    19062729 :   pari_sp av = avma;
     574    19062729 :   GEN c = mulii(gel(x,1), gel(y,1));
     575             :   long i;
     576   142982441 :   for (i = 2; i < lx; i++)
     577             :   {
     578   123921006 :     GEN t = mulii(gel(x,i), gel(y,i));
     579   123920459 :     if (t != gen_0) c = addii(c, t);
     580             :   }
     581    19061435 :   return gerepileuptoint(av, c);
     582             : }
     583             : 
     584             : /* x~ * y, assuming result is symmetric */
     585             : GEN
     586      525460 : ZM_transmultosym(GEN x, GEN y)
     587             : {
     588      525460 :   long i, j, l, ly = lg(y);
     589             :   GEN M;
     590      525460 :   if (ly == 1) return cgetg(1,t_MAT);
     591             :   /* lg(x) = ly */
     592      525460 :   l = lgcols(y); /* = lgcols(x) */
     593      525460 :   M = cgetg(ly, t_MAT);
     594     2674439 :   for (i=1; i<ly; i++)
     595             :   {
     596     2148979 :     GEN xi = gel(x,i), c = cgetg(ly,t_COL);
     597     2148979 :     gel(M,i) = c;
     598     7011023 :     for (j=1; j<i; j++)
     599     4862044 :       gcoeff(M,i,j) = gel(c,j) = ZV_dotproduct_i(xi,gel(y,j),l);
     600     2148979 :     gel(c,i) = ZV_dotproduct_i(xi,gel(y,i),l);
     601             :   }
     602      525460 :   return M;
     603             : }
     604             : /* x~ * y */
     605             : GEN
     606        2485 : ZM_transmul(GEN x, GEN y)
     607             : {
     608        2485 :   long i, j, l, lx, ly = lg(y);
     609             :   GEN M;
     610        2485 :   if (ly == 1) return cgetg(1,t_MAT);
     611        2485 :   lx = lg(x);
     612        2485 :   l = lgcols(y);
     613        2485 :   if (lgcols(x) != l) pari_err_OP("operation 'ZM_transmul'", x,y);
     614        2485 :   M = cgetg(ly, t_MAT);
     615        7581 :   for (i=1; i<ly; i++)
     616             :   {
     617        5096 :     GEN yi = gel(y,i), c = cgetg(lx,t_COL);
     618        5096 :     gel(M,i) = c;
     619       13013 :     for (j=1; j<lx; j++) gel(c,j) = ZV_dotproduct_i(yi,gel(x,j),l);
     620             :   }
     621        2485 :   return M;
     622             : }
     623             : 
     624             : static GEN
     625       25165 : ZM_sqr_i(GEN x, long l)
     626             : {
     627       25165 :   long s = ZM_max_lg_i(x, l, l);
     628       25165 :   if (l > 70) return ZM_mul_fast(x, x, l, l, s, s);
     629       25165 :   if (l <= sw_bound(s))
     630       25135 :     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       25165 : ZM_sqr(GEN x)
     637             : {
     638       25165 :   long lx=lg(x);
     639       25165 :   if (lx==1) return cgetg(1,t_MAT);
     640       25165 :   return ZM_sqr_i(x, lx);
     641             : }
     642             : GEN
     643    23144214 : ZM_ZC_mul(GEN x, GEN y)
     644             : {
     645    23144214 :   long lx = lg(x);
     646    23144214 :   return lx==1? cgetg(1,t_COL): ZM_ZC_mul_i(x, y, lx, lgcols(x));
     647             : }
     648             : 
     649             : GEN
     650     4470784 : ZC_Z_div(GEN x, GEN c)
     651    19214401 : { pari_APPLY_type(t_COL, Qdivii(gel(x,i), c)) }
     652             : 
     653             : GEN
     654       14935 : ZM_Z_div(GEN x, GEN c)
     655      168293 : { pari_APPLY_same(ZC_Z_div(gel(x, i), c)) }
     656             : 
     657             : GEN
     658     2625975 : ZC_Q_mul(GEN A, GEN z)
     659             : {
     660     2625975 :   pari_sp av = avma;
     661     2625975 :   long i, l = lg(A);
     662             :   GEN d, n, Ad, B, u;
     663     2625975 :   if (typ(z)==t_INT) return ZC_Z_mul(A,z);
     664     2623462 :   n = gel(z, 1); d = gel(z, 2);
     665     2623462 :   Ad = FpC_red(A, d);
     666     2623407 :   u = gcdii(d, FpV_factorback(Ad, NULL, d));
     667     2623425 :   B = cgetg(l, t_COL);
     668     2623424 :   if (equali1(u))
     669             :   {
     670      319265 :     for(i=1; i<l; i++)
     671      257584 :       gel(B, i) = mkfrac(mulii(n, gel(A,i)), d);
     672             :   } else
     673             :   {
     674    17558166 :     for(i=1; i<l; i++)
     675             :     {
     676    14996500 :       GEN di = gcdii(gel(Ad, i), u), ni = mulii(n, diviiexact(gel(A,i), di));
     677    14996426 :       if (equalii(d, di))
     678    10334775 :         gel(B, i) = ni;
     679             :       else
     680     4661630 :         gel(B, i) = mkfrac(ni, diviiexact(d, di));
     681             :     }
     682             :   }
     683     2623347 :   return gerepilecopy(av, B);
     684             : }
     685             : 
     686             : GEN
     687     1098426 : ZM_Q_mul(GEN x, GEN z)
     688             : {
     689     1098426 :   if (typ(z)==t_INT) return ZM_Z_mul(x,z);
     690     3158618 :   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      212436 : ZV_ZM_mul(GEN x, GEN y)
     707             : {
     708      212436 :   long i, lx = lg(x), ly = lg(y);
     709             :   GEN z;
     710      212436 :   if (lx == 1) return zerovec(ly-1);
     711      212317 :   z = cgetg(ly, t_VEC);
     712      845215 :   for (i = 1; i < ly; i++) gel(z,i) = ZV_dotproduct_i(x, gel(y,i), lx);
     713      212317 :   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     6680516 : ZV_dotsquare(GEN x)
     733             : {
     734             :   long i, lx;
     735             :   pari_sp av;
     736             :   GEN z;
     737     6680516 :   lx = lg(x);
     738     6680516 :   if (lx == 1) return gen_0;
     739     6680516 :   av = avma; z = sqri(gel(x,1));
     740    26425435 :   for (i=2; i<lx; i++) z = addii(z, sqri(gel(x,i)));
     741     6678903 :   return gerepileuptoint(av,z);
     742             : }
     743             : 
     744             : GEN
     745    16206749 : ZV_dotproduct(GEN x,GEN y)
     746             : {
     747             :   long lx;
     748    16206749 :   if (x == y) return ZV_dotsquare(x);
     749    11409966 :   lx = lg(x);
     750    11409966 :   if (lx == 1) return gen_0;
     751    11409966 :   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    34406125 : ZC_add_i(GEN x, GEN y, long lx)
     781             : {
     782    34406125 :   GEN A = cgetg(lx, t_COL);
     783             :   long i;
     784   440408592 :   for (i=1; i<lx; i++) gel(A,i) = addii(gel(x,i), gel(y,i));
     785    34399418 :   return A;
     786             : }
     787             : GEN
     788    26407821 : ZC_add(GEN x, GEN y) { return ZC_add_i(x, y, lg(x)); }
     789             : GEN
     790      363006 : ZC_Z_add(GEN x, GEN y)
     791             : {
     792      363006 :   long k, lx = lg(x);
     793      363006 :   GEN z = cgetg(lx, t_COL);
     794      363007 :   if (lx == 1) pari_err_TYPE2("+",x,y);
     795      363007 :   gel(z,1) = addii(y,gel(x,1));
     796     2499930 :   for (k = 2; k < lx; k++) gel(z,k) = icopy(gel(x,k));
     797      362995 :   return z;
     798             : }
     799             : 
     800             : static GEN
     801     9550159 : ZC_sub_i(GEN x, GEN y, long lx)
     802             : {
     803             :   long i;
     804     9550159 :   GEN A = cgetg(lx, t_COL);
     805    66772975 :   for (i=1; i<lx; i++) gel(A,i) = subii(gel(x,i), gel(y,i));
     806     9549713 :   return A;
     807             : }
     808             : GEN
     809     9292319 : 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      545207 : Z_ZC_sub(GEN a, GEN x)
     822             : {
     823      545207 :   long k, lx = lg(x);
     824      545207 :   GEN z = cgetg(lx, t_COL);
     825      545216 :   if (lx == 1) pari_err_TYPE2("-",a,x);
     826      545216 :   gel(z,1) = subii(a, gel(x,1));
     827     1510793 :   for (k = 2; k < lx; k++) gel(z,k) = negi(gel(x,k));
     828      545215 :   return z;
     829             : }
     830             : 
     831             : GEN
     832      741338 : ZM_add(GEN x, GEN y)
     833             : {
     834      741338 :   long lx = lg(x), l, j;
     835             :   GEN z;
     836      741338 :   if (lx == 1) return cgetg(1, t_MAT);
     837      662058 :   z = cgetg(lx, t_MAT); l = lgcols(x);
     838     8660190 :   for (j = 1; j < lx; j++) gel(z,j) = ZC_add_i(gel(x,j), gel(y,j), l);
     839      662059 :   return z;
     840             : }
     841             : GEN
     842       77338 : ZM_sub(GEN x, GEN y)
     843             : {
     844       77338 :   long lx = lg(x), l, j;
     845             :   GEN z;
     846       77338 :   if (lx == 1) return cgetg(1, t_MAT);
     847       77338 :   z = cgetg(lx, t_MAT); l = lgcols(x);
     848      335065 :   for (j = 1; j < lx; j++) gel(z,j) = ZC_sub_i(gel(x,j), gel(y,j), l);
     849       77338 :   return z;
     850             : }
     851             : /********************************************************************/
     852             : /**                                                                **/
     853             : /**                         LINEAR COMBINATION                     **/
     854             : /**                                                                **/
     855             : /********************************************************************/
     856             : /* return X/c assuming division is exact */
     857             : GEN
     858     4409815 : ZC_Z_divexact(GEN x, GEN c)
     859    45587391 : { 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      427350 : ZM_Z_divexact(GEN x, GEN c)
     866     2778562 : { 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    35182255 : ZC_Z_mul(GEN x, GEN c)
     874             : {
     875    35182255 :   if (!signe(c)) return zerocol(lg(x)-1);
     876    33827279 :   if (is_pm1(c)) return (signe(c) > 0)? ZC_copy(x): ZC_neg(x);
     877   254459758 :   pari_APPLY_type(t_COL, mulii(gel(x,i), c))
     878             : }
     879             : 
     880             : GEN
     881       61069 : ZC_z_mul(GEN x, long c)
     882             : {
     883       61069 :   if (!c) return zerocol(lg(x)-1);
     884       52582 :   if (c == 1) return ZC_copy(x);
     885       47728 :   if (c ==-1) return ZC_neg(x);
     886      476030 :   pari_APPLY_type(t_COL, mulsi(c, gel(x,i)))
     887             : }
     888             : 
     889             : GEN
     890       28287 : zv_z_mul(GEN x, long n)
     891      430480 : { 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     3006401 : ZM_Z_mul(GEN X, GEN c)
     917             : {
     918     3006401 :   long i, j, h, l = lg(X);
     919             :   GEN A;
     920     3006401 :   if (l == 1) return cgetg(1, t_MAT);
     921     3006401 :   h = lgcols(X);
     922     3006393 :   if (!signe(c)) return zeromat(h-1, l-1);
     923     2961124 :   if (is_pm1(c)) return (signe(c) > 0)? ZM_copy(X): ZM_neg(X);
     924     2675432 :   A = cgetg(l, t_MAT);
     925    13459641 :   for (j = 1; j < l; j++)
     926             :   {
     927    10784628 :     GEN a = cgetg(h, t_COL), x = gel(X, j);
     928   168729340 :     for (i = 1; i < h; i++) gel(a,i) = mulii(c, gel(x,i));
     929    10784206 :     gel(A,j) = a;
     930             :   }
     931     2675013 :   return A;
     932             : }
     933             : void
     934    93604828 : ZC_lincomb1_inplace_i(GEN X, GEN Y, GEN v, long n)
     935             : {
     936             :   long i;
     937  2033512365 :   for (i = n; i; i--) gel(X,i) = addmulii_inplace(gel(X,i), gel(Y,i), v);
     938    93180442 : }
     939             : /* X <- X + v Y (elementary col operation) */
     940             : void
     941    82804819 : ZC_lincomb1_inplace(GEN X, GEN Y, GEN v)
     942             : {
     943    82804819 :   if (lgefint(v) != 2) return ZC_lincomb1_inplace_i(X, Y, v, lg(X)-1);
     944             : }
     945             : void
     946    29504029 : Flc_lincomb1_inplace(GEN X, GEN Y, ulong v, ulong q)
     947             : {
     948             :   long i;
     949    29504029 :   if (!v) return; /* v = 0 */
     950   742595938 :   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    17675424 : ZC_lincomb1(GEN v, GEN x, GEN y)
     956   137205625 : { pari_APPLY_type(t_COL, addmulii(gel(x,i), gel(y,i), v)) }
     957             : 
     958             : /* -X + vY */
     959             : static GEN
     960      733082 : ZC_lincomb_1(GEN v, GEN x, GEN y)
     961     4531410 : { 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    35944506 : ZC_lincomb(GEN u, GEN v, GEN X, GEN Y)
     966             : {
     967             :   long su, sv;
     968             :   GEN A;
     969             : 
     970    35944506 :   su = signe(u); if (!su) return ZC_Z_mul(Y, v);
     971    35943547 :   sv = signe(v); if (!sv) return ZC_Z_mul(X, u);
     972    35514457 :   if (is_pm1(v))
     973             :   {
     974    11647624 :     if (is_pm1(u))
     975             :     {
     976    10332527 :       if (su != sv) A = ZC_sub(X, Y);
     977     2741977 :       else          A = ZC_add(X, Y);
     978    10332176 :       if (su < 0) ZV_togglesign(A); /* in place but was created above */
     979             :     }
     980             :     else
     981             :     {
     982     1315005 :       if (sv > 0) A = ZC_lincomb1 (u, Y, X);
     983      603323 :       else        A = ZC_lincomb_1(u, Y, X);
     984             :     }
     985             :   }
     986    23867659 :   else if (is_pm1(u))
     987             :   {
     988    17092769 :     if (su > 0) A = ZC_lincomb1 (v, X, Y);
     989      128966 :     else        A = ZC_lincomb_1(v, X, Y);
     990             :   }
     991             :   else
     992             :   { /* not cgetg_copy: x may be a t_VEC */
     993     6779810 :     long i, lx = lg(X);
     994     6779810 :     A = cgetg(lx,t_COL);
     995    43978360 :     for (i=1; i<lx; i++) gel(A,i) = lincombii(u,v,gel(X,i),gel(Y,i));
     996             :   }
     997    35507280 :   return A;
     998             : }
     999             : 
    1000             : /********************************************************************/
    1001             : /**                                                                **/
    1002             : /**                           CONVERSIONS                          **/
    1003             : /**                                                                **/
    1004             : /********************************************************************/
    1005             : GEN
    1006      393946 : ZV_to_nv(GEN x)
    1007      739953 : { pari_APPLY_ulong(itou(gel(x,i))) }
    1008             : 
    1009             : GEN
    1010      137173 : zm_to_ZM(GEN x)
    1011      759352 : { 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    13301862 : ZC_copy(GEN x)
    1041             : {
    1042    13301862 :   long i, lx = lg(x);
    1043    13301862 :   GEN y = cgetg(lx, t_COL);
    1044    83537866 :   for (i=1; i<lx; i++)
    1045             :   {
    1046    70235441 :     GEN c = gel(x,i);
    1047    70235441 :     gel(y,i) = lgefint(c) == 2? gen_0: icopy(c);
    1048             :   }
    1049    13302425 :   return y;
    1050             : }
    1051             : 
    1052             : GEN
    1053      518734 : ZM_copy(GEN x)
    1054     4310704 : { pari_APPLY_same(ZC_copy(gel(x,i))) }
    1055             : 
    1056             : void
    1057      344209 : ZV_neg_inplace(GEN M)
    1058             : {
    1059      344209 :   long l = lg(M);
    1060     1265752 :   while (--l > 0) gel(M,l) = negi(gel(M,l));
    1061      344299 : }
    1062             : GEN
    1063     6043483 : ZC_neg(GEN x)
    1064    35516430 : { pari_APPLY_type(t_COL, negi(gel(x,i))) }
    1065             : 
    1066             : GEN
    1067       51084 : zv_neg(GEN x)
    1068      655372 : { 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     1674762 : ZM_neg(GEN x)
    1081     5081955 : { pari_APPLY_same(ZC_neg(gel(x,i))) }
    1082             : 
    1083             : void
    1084     4987214 : ZV_togglesign(GEN M)
    1085             : {
    1086     4987214 :   long l = lg(M);
    1087    74103017 :   while (--l > 0) togglesign_safe(&gel(M,l));
    1088     4987376 : }
    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    11314066 : ZC_hnfremdiv(GEN x, GEN y, GEN *Q)
    1104             : {
    1105    11314066 :   long i, l = lg(x);
    1106             :   GEN q;
    1107             : 
    1108    11314066 :   if (Q) *Q = cgetg(l,t_COL);
    1109    11314066 :   if (l == 1) return cgetg(1,t_COL);
    1110    64632920 :   for (i = l-1; i>0; i--)
    1111             :   {
    1112    53322094 :     q = diviiround(gel(x,i), gcoeff(y,i,i));
    1113    53319900 :     if (signe(q)) {
    1114    25896705 :       togglesign(q);
    1115    25896901 :       x = ZC_lincomb(gen_1, q, x, gel(y,i));
    1116             :     }
    1117    53318854 :     if (Q) gel(*Q, i) = q;
    1118             :   }
    1119    11310826 :   return x;
    1120             : }
    1121             : 
    1122             : /* x = y Q + R, may return some columns of x (not copies) */
    1123             : GEN
    1124      453663 : ZM_hnfdivrem(GEN x, GEN y, GEN *Q)
    1125             : {
    1126      453663 :   long lx = lg(x), i;
    1127      453663 :   GEN R = cgetg(lx, t_MAT);
    1128      453697 :   if (Q)
    1129             :   {
    1130      127242 :     GEN q = cgetg(lx, t_MAT); *Q = q;
    1131      188038 :     for (i=1; i<lx; i++) gel(R,i) = ZC_hnfremdiv(gel(x,i),y,(GEN*)(q+i));
    1132             :   }
    1133             :   else
    1134      822990 :     for (i=1; i<lx; i++)
    1135             :     {
    1136      496529 :       pari_sp av = avma;
    1137      496529 :       GEN z = ZC_hnfrem(gel(x,i),y);
    1138      496501 :       gel(R,i) = (avma == av)? ZC_copy(z): gerepileupto(av, z);
    1139             :     }
    1140      453703 :   return R;
    1141             : }
    1142             : 
    1143             : /********************************************************************/
    1144             : /**                                                                **/
    1145             : /**                               TESTS                            **/
    1146             : /**                                                                **/
    1147             : /********************************************************************/
    1148             : int
    1149    23099592 : zv_equal0(GEN V)
    1150             : {
    1151    23099592 :   long l = lg(V);
    1152    37561367 :   while (--l > 0)
    1153    30788055 :     if (V[l]) return 0;
    1154     6773312 :   return 1;
    1155             : }
    1156             : 
    1157             : int
    1158    14503196 : ZV_equal0(GEN V)
    1159             : {
    1160    14503196 :   long l = lg(V);
    1161    25632319 :   while (--l > 0)
    1162    25205979 :     if (signe(gel(V,l))) return 0;
    1163      426340 :   return 1;
    1164             : }
    1165             : int
    1166       16185 : ZMrow_equal0(GEN V, long i)
    1167             : {
    1168       16185 :   long l = lg(V);
    1169       25283 :   while (--l > 0)
    1170       21717 :     if (signe(gcoeff(V,i,l))) return 0;
    1171        3566 :   return 1;
    1172             : }
    1173             : 
    1174             : static int
    1175     6315120 : ZV_equal_lg(GEN V, GEN W, long l)
    1176             : {
    1177    26657068 :   while (--l > 0)
    1178    20822496 :     if (!equalii(gel(V,l), gel(W,l))) return 0;
    1179     5834572 :   return 1;
    1180             : }
    1181             : int
    1182      291963 : ZV_equal(GEN V, GEN W)
    1183             : {
    1184      291963 :   long l = lg(V);
    1185      291963 :   if (lg(W) != l) return 0;
    1186      291956 :   return ZV_equal_lg(V, W, l);
    1187             : }
    1188             : int
    1189     3347679 : ZM_equal(GEN A, GEN B)
    1190             : {
    1191     3347679 :   long i, m, l = lg(A);
    1192     3347679 :   if (lg(B) != l) return 0;
    1193     3347544 :   if (l == 1) return 1;
    1194     3347544 :   m = lgcols(A);
    1195     3347539 :   if (lgcols(B) != m) return 0;
    1196     9110974 :   for (i = 1; i < l; i++)
    1197     6023161 :     if (!ZV_equal_lg(gel(A,i), gel(B,i), m)) return 0;
    1198     3087813 :   return 1;
    1199             : }
    1200             : int
    1201       70223 : ZM_equal0(GEN A)
    1202             : {
    1203       70223 :   long i, j, m, l = lg(A);
    1204       70223 :   if (l == 1) return 1;
    1205       70223 :   m = lgcols(A);
    1206      121034 :   for (j = 1; j < l; j++)
    1207     2616071 :     for (i = 1; i < m; i++)
    1208     2565260 :       if (signe(gcoeff(A,i,j))) return 0;
    1209       15550 :   return 1;
    1210             : }
    1211             : int
    1212    30833181 : zv_equal(GEN V, GEN W)
    1213             : {
    1214    30833181 :   long l = lg(V);
    1215    30833181 :   if (lg(W) != l) return 0;
    1216   262055904 :   while (--l > 0)
    1217   232340629 :     if (V[l] != W[l]) return 0;
    1218    29715275 :   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      749103 : ZM_ishnf(GEN x)
    1233             : {
    1234      749103 :   long i,j, lx = lg(x);
    1235     2267320 :   for (i=1; i<lx; i++)
    1236             :   {
    1237     1621771 :     GEN xii = gcoeff(x,i,i);
    1238     1621771 :     if (signe(xii) <= 0) return 0;
    1239     3163846 :     for (j=1; j<i; j++)
    1240     1572343 :       if (signe(gcoeff(x,i,j))) return 0;
    1241     3221225 :     for (j=i+1; j<lx; j++)
    1242             :     {
    1243     1703008 :       GEN xij = gcoeff(x,i,j);
    1244     1703008 :       if (signe(xij)<0 || cmpii(xij,xii)>=0) return 0;
    1245             :     }
    1246             :   }
    1247      645549 :   return 1;
    1248             : }
    1249             : int
    1250      638790 : ZM_isidentity(GEN x)
    1251             : {
    1252      638790 :   long i,j, lx = lg(x);
    1253             : 
    1254      638790 :   if (lx == 1) return 1;
    1255      638783 :   if (lx != lgcols(x)) return 0;
    1256     3120719 :   for (j=1; j<lx; j++)
    1257             :   {
    1258     2483974 :     GEN c = gel(x,j);
    1259     7700736 :     for (i=1; i<j; )
    1260     5216798 :       if (signe(gel(c,i++))) return 0;
    1261             :     /* i = j */
    1262     2483938 :     if (!equali1(gel(c,i++))) return 0;
    1263     7706539 :     for (   ; i<lx; )
    1264     5224589 :       if (signe(gel(c,i++))) return 0;
    1265             :   }
    1266      636745 :   return 1;
    1267             : }
    1268             : int
    1269      555778 : ZM_isdiagonal(GEN x)
    1270             : {
    1271      555778 :   long i,j, lx = lg(x);
    1272      555778 :   if (lx == 1) return 1;
    1273      555778 :   if (lx != lgcols(x)) return 0;
    1274             : 
    1275     1436648 :   for (j=1; j<lx; j++)
    1276             :   {
    1277     1207355 :     GEN c = gel(x,j);
    1278     1694690 :     for (i=1; i<j; i++)
    1279      813778 :       if (signe(gel(c,i))) return 0;
    1280     2017830 :     for (i++; i<lx; i++)
    1281     1136953 :       if (signe(gel(c,i))) return 0;
    1282             :   }
    1283      229293 :   return 1;
    1284             : }
    1285             : int
    1286      113620 : ZM_isscalar(GEN x, GEN s)
    1287             : {
    1288      113620 :   long i, j, lx = lg(x);
    1289             : 
    1290      113620 :   if (lx == 1) return 1;
    1291      113620 :   if (!s) s = gcoeff(x,1,1);
    1292      113620 :   if (equali1(s)) return ZM_isidentity(x);
    1293      112409 :   if (lx != lgcols(x)) return 0;
    1294      565455 :   for (j=1; j<lx; j++)
    1295             :   {
    1296      476382 :     GEN c = gel(x,j);
    1297     2265579 :     for (i=1; i<j; )
    1298     1810082 :       if (signe(gel(c,i++))) return 0;
    1299             :     /* i = j */
    1300      455497 :     if (!equalii(gel(c,i++), s)) return 0;
    1301     2315962 :     for (   ; i<lx; )
    1302     1862916 :       if (signe(gel(c,i++))) return 0;
    1303             :   }
    1304       89073 :   return 1;
    1305             : }
    1306             : 
    1307             : long
    1308      105829 : ZC_is_ei(GEN x)
    1309             : {
    1310      105829 :   long i, j = 0, l = lg(x);
    1311      979320 :   for (i = 1; i < l; i++)
    1312             :   {
    1313      873492 :     GEN c = gel(x,i);
    1314      873492 :     long s = signe(c);
    1315      873492 :     if (!s) continue;
    1316      105823 :     if (s < 0 || !is_pm1(c) || j) return 0;
    1317      105822 :     j = i;
    1318             :   }
    1319      105828 :   return j;
    1320             : }
    1321             : 
    1322             : /********************************************************************/
    1323             : /**                                                                **/
    1324             : /**                       MISCELLANEOUS                            **/
    1325             : /**                                                                **/
    1326             : /********************************************************************/
    1327             : /* assume lg(x) = lg(y), x,y in Z^n */
    1328             : int
    1329     3133870 : ZV_cmp(GEN x, GEN y)
    1330             : {
    1331     3133870 :   long fl,i, lx = lg(x);
    1332     6347338 :   for (i=1; i<lx; i++)
    1333     5059522 :     if (( fl = cmpii(gel(x,i), gel(y,i)) )) return fl;
    1334     1287816 :   return 0;
    1335             : }
    1336             : /* assume lg(x) = lg(y), x,y in Z^n */
    1337             : int
    1338       19677 : ZV_abscmp(GEN x, GEN y)
    1339             : {
    1340       19677 :   long fl,i, lx = lg(x);
    1341       53895 :   for (i=1; i<lx; i++)
    1342       53768 :     if (( fl = abscmpii(gel(x,i), gel(y,i)) )) return fl;
    1343         127 :   return 0;
    1344             : }
    1345             : 
    1346             : long
    1347    20568949 : zv_content(GEN x)
    1348             : {
    1349    20568949 :   long i, s, l = lg(x);
    1350    20568949 :   if (l == 1) return 0;
    1351    20568942 :   s = labs(x[1]);
    1352    46236234 :   for (i = 2; i < l && s != 1; i++) s = ugcd(s, labs(x[i]));
    1353    20570246 :   return s;
    1354             : }
    1355             : GEN
    1356      405308 : ZV_content(GEN x)
    1357             : {
    1358      405308 :   long i, l = lg(x);
    1359      405308 :   pari_sp av = avma;
    1360             :   GEN c;
    1361      405308 :   if (l == 1) return gen_0;
    1362      405308 :   if (l == 2) return absi(gel(x,1));
    1363      330667 :   c = gel(x,1);
    1364      886117 :   for (i = 2; i < l; i++)
    1365             :   {
    1366      657630 :     c = gcdii(c, gel(x,i));
    1367      657630 :     if (is_pm1(c)) { set_avma(av); return gen_1; }
    1368             :   }
    1369      228487 :   return gerepileuptoint(av, c);
    1370             : }
    1371             : 
    1372             : GEN
    1373     3868454 : ZM_det_triangular(GEN mat)
    1374             : {
    1375             :   pari_sp av;
    1376     3868454 :   long i,l = lg(mat);
    1377             :   GEN s;
    1378             : 
    1379     3868454 :   if (l<3) return l<2? gen_1: icopy(gcoeff(mat,1,1));
    1380     3465798 :   av = avma; s = gcoeff(mat,1,1);
    1381     9423248 :   for (i=2; i<l; i++) s = mulii(s,gcoeff(mat,i,i));
    1382     3465347 :   return gerepileuptoint(av,s);
    1383             : }
    1384             : 
    1385             : /* assumes no overflow */
    1386             : long
    1387      939941 : zv_prod(GEN v)
    1388             : {
    1389      939941 :   long n, i, l = lg(v);
    1390      939941 :   if (l == 1) return 1;
    1391      952265 :   n = v[1]; for (i = 2; i < l; i++) n *= v[i];
    1392      763905 :   return n;
    1393             : }
    1394             : 
    1395             : static GEN
    1396   356789841 : _mulii(void *E, GEN a, GEN b)
    1397   356789841 : { (void) E; return mulii(a, b); }
    1398             : 
    1399             : /* product of ulongs */
    1400             : GEN
    1401     1862856 : zv_prod_Z(GEN v)
    1402             : {
    1403     1862856 :   pari_sp av = avma;
    1404     1862856 :   long k, m, n = lg(v)-1;
    1405             :   GEN V;
    1406     1862856 :   switch(n) {
    1407       20902 :     case 0: return gen_1;
    1408      129878 :     case 1: return utoi(v[1]);
    1409      975949 :     case 2: return muluu(v[1], v[2]);
    1410             :   }
    1411      736127 :   m = n >> 1;
    1412      736127 :   V = cgetg(m + (odd(n)? 2: 1), t_VEC);
    1413    85829300 :   for (k = 1; k <= m; k++) gel(V,k) = muluu(v[k<<1], v[(k<<1)-1]);
    1414      736100 :   if (odd(n)) gel(V,k) = utoipos(v[n]);
    1415      736123 :   return gerepileuptoint(av, gen_product(V, NULL, &_mulii));
    1416             : }
    1417             : GEN
    1418    14694393 : vecsmall_prod(GEN v)
    1419             : {
    1420    14694393 :   pari_sp av = avma;
    1421    14694393 :   long k, m, n = lg(v)-1;
    1422             :   GEN V;
    1423    14694393 :   switch (n) {
    1424           0 :     case 0: return gen_1;
    1425           0 :     case 1: return stoi(v[1]);
    1426          21 :     case 2: return mulss(v[1], v[2]);
    1427             :   }
    1428    14694372 :   m = n >> 1;
    1429    14694372 :   V = cgetg(m + (odd(n)? 2: 1), t_VEC);
    1430   161556906 :   for (k = 1; k <= m; k++) gel(V,k) = mulss(v[k<<1], v[(k<<1)-1]);
    1431    14694372 :   if (odd(n)) gel(V,k) = stoi(v[n]);
    1432    14694372 :   return gerepileuptoint(av, gen_product(V, NULL, &_mulii));
    1433             : }
    1434             : 
    1435             : GEN
    1436     8814271 : ZV_prod(GEN v)
    1437             : {
    1438     8814271 :   pari_sp av = avma;
    1439     8814271 :   long i, l = lg(v);
    1440             :   GEN n;
    1441     8814271 :   if (l == 1) return gen_1;
    1442     8632653 :   if (l > 7) return gerepileuptoint(av, gen_product(v, NULL, _mulii));
    1443     1289826 :   n = gel(v,1);
    1444     1289826 :   if (l == 2) return icopy(n);
    1445     2090609 :   for (i = 2; i < l; i++) n = mulii(n, gel(v,i));
    1446      840136 :   return gerepileuptoint(av, n);
    1447             : }
    1448             : /* assumes no overflow */
    1449             : long
    1450       16428 : zv_sum(GEN v)
    1451             : {
    1452       16428 :   long n, i, l = lg(v);
    1453       16428 :   if (l == 1) return 0;
    1454       95624 :   n = v[1]; for (i = 2; i < l; i++) n += v[i];
    1455       16407 :   return n;
    1456             : }
    1457             : /* assumes no overflow and 0 <= n <= #v */
    1458             : long
    1459           0 : zv_sumpart(GEN v, long n)
    1460             : {
    1461             :   long i, p;
    1462           0 :   if (!n) return 0;
    1463           0 :   p = v[1]; for (i = 2; i <= n; i++) p += v[i];
    1464           0 :   return p;
    1465             : }
    1466             : GEN
    1467          77 : ZV_sum(GEN v)
    1468             : {
    1469          77 :   pari_sp av = avma;
    1470          77 :   long i, l = lg(v);
    1471             :   GEN n;
    1472          77 :   if (l == 1) return gen_0;
    1473          77 :   n = gel(v,1);
    1474          77 :   if (l == 2) return icopy(n);
    1475         581 :   for (i = 2; i < l; i++) n = addii(n, gel(v,i));
    1476          77 :   return gerepileuptoint(av, n);
    1477             : }
    1478             : 
    1479             : /********************************************************************/
    1480             : /**                                                                **/
    1481             : /**         GRAM SCHMIDT REDUCTION (integer matrices)              **/
    1482             : /**                                                                **/
    1483             : /********************************************************************/
    1484             : 
    1485             : /* L[k,] += q * L[l,], l < k. Inefficient if q = 0 */
    1486             : static void
    1487      358766 : Zupdate_row(long k, long l, GEN q, GEN L, GEN B)
    1488             : {
    1489      358766 :   long i, qq = itos_or_0(q);
    1490      358765 :   if (!qq)
    1491             :   {
    1492       33347 :     for(i=1;i<l;i++)  gcoeff(L,k,i) = addii(gcoeff(L,k,i),mulii(q,gcoeff(L,l,i)));
    1493        7069 :     gcoeff(L,k,l) = addii(gcoeff(L,k,l), mulii(q,B));
    1494        7069 :     return;
    1495             :   }
    1496      351696 :   if (qq == 1) {
    1497      150148 :     for (i=1;i<l; i++) gcoeff(L,k,i) = addii(gcoeff(L,k,i),gcoeff(L,l,i));
    1498      102031 :     gcoeff(L,k,l) = addii(gcoeff(L,k,l), B);
    1499      249665 :   } else if (qq == -1) {
    1500      150839 :     for (i=1;i<l; i++) gcoeff(L,k,i) = subii(gcoeff(L,k,i),gcoeff(L,l,i));
    1501       87138 :     gcoeff(L,k,l) = addii(gcoeff(L,k,l), negi(B));
    1502             :   } else {
    1503      282606 :     for(i=1;i<l;i++) gcoeff(L,k,i) = addii(gcoeff(L,k,i),mulsi(qq,gcoeff(L,l,i)));
    1504      162549 :     gcoeff(L,k,l) = addii(gcoeff(L,k,l), mulsi(qq,B));
    1505             :   }
    1506             : }
    1507             : 
    1508             : /* update L[k,] */
    1509             : static void
    1510     1032961 : ZRED(long k, long l, GEN x, GEN L, GEN B)
    1511             : {
    1512     1032961 :   GEN q = truedivii(addii(B,shifti(gcoeff(L,k,l),1)), shifti(B,1));
    1513     1033067 :   if (!signe(q)) return;
    1514      358728 :   q = negi(q);
    1515      358758 :   Zupdate_row(k,l,q,L,B);
    1516      358748 :   gel(x,k) = ZC_lincomb(gen_1, q, gel(x,k), gel(x,l));
    1517             : }
    1518             : 
    1519             : /* Gram-Schmidt reduction, x a ZM */
    1520             : static void
    1521     1183476 : ZincrementalGS(GEN x, GEN L, GEN B, long k)
    1522             : {
    1523             :   long i, j;
    1524     3778673 :   for (j=1; j<=k; j++)
    1525             :   {
    1526     2596206 :     pari_sp av = avma;
    1527     2596206 :     GEN u = ZV_dotproduct(gel(x,k), gel(x,j));
    1528     5610354 :     for (i=1; i<j; i++)
    1529             :     {
    1530     3015996 :       u = subii(mulii(gel(B,i+1), u), mulii(gcoeff(L,k,i), gcoeff(L,j,i)));
    1531     3015378 :       u = diviiexact(u, gel(B,i));
    1532             :     }
    1533     2594358 :     gcoeff(L,k,j) = gerepileuptoint(av, u);
    1534             :   }
    1535     1182467 :   gel(B,k+1) = gcoeff(L,k,k); gcoeff(L,k,k) = gen_1;
    1536     1182467 : }
    1537             : 
    1538             : /* Variant reducemodinvertible(ZC v, ZM y), when y singular.
    1539             :  * Very inefficient if y is not LLL-reduced of maximal rank */
    1540             : static GEN
    1541      110410 : ZC_reducemodmatrix_i(GEN v, GEN y)
    1542             : {
    1543      110410 :   GEN B, L, x = shallowconcat(y, v);
    1544      110408 :   long k, lx = lg(x), nx = lx-1;
    1545             : 
    1546      110408 :   B = scalarcol_shallow(gen_1, lx);
    1547      110408 :   L = zeromatcopy(nx, nx);
    1548      454176 :   for (k=1; k <= nx; k++) ZincrementalGS(x, L, B, k);
    1549      343774 :   for (k = nx-1; k >= 1; k--) ZRED(nx,k, x,L,gel(B,k+1));
    1550      110403 :   return gel(x,nx);
    1551             : }
    1552             : GEN
    1553      110410 : ZC_reducemodmatrix(GEN v, GEN y) {
    1554      110410 :   pari_sp av = avma;
    1555      110410 :   return gerepilecopy(av, ZC_reducemodmatrix_i(v,y));
    1556             : }
    1557             : 
    1558             : /* Variant reducemodinvertible(ZM v, ZM y), when y singular.
    1559             :  * Very inefficient if y is not LLL-reduced of maximal rank */
    1560             : static GEN
    1561      226860 : ZM_reducemodmatrix_i(GEN v, GEN y)
    1562             : {
    1563             :   GEN B, L, V;
    1564      226860 :   long j, k, lv = lg(v), nx = lg(y), lx = nx+1;
    1565             : 
    1566      226860 :   V = cgetg(lv, t_MAT);
    1567      226880 :   B = scalarcol_shallow(gen_1, lx);
    1568      226881 :   L = zeromatcopy(nx, nx);
    1569      600933 :   for (k=1; k < nx; k++) ZincrementalGS(y, L, B, k);
    1570      692532 :   for (j = 1; j < lg(v); j++)
    1571             :   {
    1572      465704 :     GEN x = shallowconcat(y, gel(v,j));
    1573      465723 :     ZincrementalGS(x, L, B, nx); /* overwrite last */
    1574     1265350 :     for (k = nx-1; k >= 1; k--) ZRED(nx,k, x,L,gel(B,k+1));
    1575      465670 :     gel(V,j) = gel(x,nx);
    1576             :   }
    1577      226828 :   return V;
    1578             : }
    1579             : GEN
    1580      226860 : ZM_reducemodmatrix(GEN v, GEN y) {
    1581      226860 :   pari_sp av = avma;
    1582      226860 :   return gerepilecopy(av, ZM_reducemodmatrix_i(v,y));
    1583             : }
    1584             : 
    1585             : GEN
    1586       98485 : ZC_reducemodlll(GEN x,GEN y)
    1587             : {
    1588       98485 :   pari_sp av = avma;
    1589       98485 :   GEN z = ZC_reducemodmatrix(x, ZM_lll(y, 0.75, LLL_INPLACE));
    1590       98486 :   return gerepilecopy(av, z);
    1591             : }
    1592             : GEN
    1593           0 : ZM_reducemodlll(GEN x,GEN y)
    1594             : {
    1595           0 :   pari_sp av = avma;
    1596           0 :   GEN z = ZM_reducemodmatrix(x, ZM_lll(y, 0.75, LLL_INPLACE));
    1597           0 :   return gerepilecopy(av, z);
    1598             : }

Generated by: LCOV version 1.14