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 - bb_hnf.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.12.0 lcov report (development 23328-a3379c31c) Lines: 600 639 93.9 %
Date: 2018-12-09 05:41:42 Functions: 51 55 92.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. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : #include "pari.h"
      14             : #include "paripriv.h"
      15             : 
      16             : #define dbg_printf(lvl) if (DEBUGLEVEL >= (lvl) + 3) err_printf
      17             : 
      18             : /********************************************************************/
      19             : /**                                                                **/
      20             : /**          BLACK BOX HERMITE RINGS AND HOWELL NORMAL FORM        **/
      21             : /**                 contributed by Aurel Page (2017)               **/
      22             : /**                                                                **/
      23             : /********************************************************************/
      24             : 
      25             : /*
      26             :   bb_hermite R:
      27             :     - add(a,b): a+b
      28             :     - neg(a): -a
      29             :     - mul(a,b): a*b
      30             :     - extgcd(a,b,&small): [d,U] with d in R and U in GL_2(R) such that [0;d] = [a;b]*U.
      31             :       set small==1 to assert that U is a 'small' operation (no red needed).
      32             :     - rann(a): b in R such that b*R = {x in R | a*x==0}
      33             :     - lquo(a,b,&r): q in R such that r=a-b*q is a canonical representative
      34             :       of the image of a in R/b*R. The canonical lift of 0 must be 0.
      35             :     - unit(a): u unit in R^* such that a*u is a canonical generator of the ideal a*R
      36             :     - equal0(a): a==0?
      37             :     - equal1(a): a==1?
      38             :     - s(n): image of the small integer n in R
      39             :     - red(a): unique representative of a as an element of R
      40             : 
      41             :   op encoding of elementary operations:
      42             :     - t_VECSMALL: the corresponding permutation (vecpermute)
      43             :     - [Vecsmall([i,j])]: the transposition Ci <-> Cj
      44             :     - [Vecsmall([i]),u], u in R^*: Ci <- Ci*u
      45             :     - [Vecsmall([i,j]),a], a in R: Ci <- Ci + Cj*a
      46             :     - [Vecsmall([i,j,0]),U], U in GL_2(R): (Ci|Cj) <- (Ci|Cj)*U
      47             : */
      48             : 
      49             : struct bb_hermite
      50             : {
      51             :   GEN (*add)(void*, GEN, GEN);
      52             :   GEN (*neg)(void*, GEN);
      53             :   GEN (*mul)(void*, GEN, GEN);
      54             :   GEN (*extgcd)(void*, GEN, GEN, int*);
      55             :   GEN (*rann)(void*, GEN);
      56             :   GEN (*lquo)(void*, GEN, GEN, GEN*);
      57             :   GEN (*unit)(void*, GEN);
      58             :   int (*equal0)(GEN);
      59             :   int (*equal1)(GEN);
      60             :   GEN (*s)(void*, long);
      61             :   GEN (*red)(void*, GEN);
      62             : };
      63             : 
      64             : static GEN
      65    22489516 : _Fp_add(void *data, GEN x, GEN y) { (void) data; return addii(x,y); }
      66             : 
      67             : static GEN
      68     2757713 : _Fp_neg(void *data, GEN x) { (void) data; return negi(x); }
      69             : 
      70             : static GEN
      71    32836587 : _Fp_mul(void *data, GEN x, GEN y) { (void) data; return mulii(x,y); }
      72             : 
      73             : static GEN
      74     1309658 : _Fp_rann(void *data, GEN x)
      75             : {
      76     1309658 :   GEN d, N = (GEN)data;
      77     1309658 :   if (!signe(x)) return gen_1;
      78     1163715 :   d = gcdii(x,N);
      79     1163715 :   return modii(diviiexact(N,d),N);
      80             : }
      81             : 
      82             : static GEN
      83     2197762 : _Fp_lquo(void *data, GEN x, GEN y, GEN* r) { (void) data; return truedvmdii(x,y,r); }
      84             : 
      85             : /* D=MN, p|M => !p|a, p|N => p|a, return M */
      86             : static GEN
      87           0 : Z_split(GEN D, GEN a)
      88             : {
      89             :   long i, n;
      90             :   GEN N;
      91           0 :   n = expi(D);
      92           0 :   n = n<2 ? 1 : expu(n)+1;
      93           0 :   for (i=1;i<=n;i++)
      94           0 :     a = Fp_sqr(a,D);
      95           0 :   N = gcdii(a,D);
      96           0 :   return diviiexact(D,N);
      97             : }
      98             : 
      99             : /* c s.t. gcd(a+cb,N) = gcd(a,b,N) without factoring */
     100             : static GEN
     101           0 : Z_stab(GEN a, GEN b, GEN N)
     102             : {
     103             :   GEN g, a2, N2;
     104           0 :   g = gcdii(a,b);
     105           0 :   g = gcdii(g,N);
     106           0 :   N2 = diviiexact(N,g);
     107           0 :   a2 = diviiexact(a,g);
     108           0 :   return Z_split(N2,a2);
     109             : }
     110             : 
     111             : static GEN
     112     3353035 : _Fp_unit(void *data, GEN x)
     113             : {
     114     3353035 :   GEN g,s,v,d,N=(GEN)data,N2;
     115             :   long i;
     116     3353035 :   if (!signe(x)) return NULL;
     117     2926343 :   g = bezout(x,N,&s,&v);
     118     2926343 :   if (equali1(g) || equali1(gcdii(s,N))) return mkvec2(g,s);
     119       44632 :   N2 = diviiexact(N,g);
     120       61110 :   for (i=0; i<5; i++)
     121             :   {
     122       61110 :     s = addii(s,N2);
     123       61110 :     if (equali1(gcdii(s,N))) return mkvec2(g,s);
     124             :   }
     125           0 :   d = Z_stab(s,N2,N);
     126           0 :   d = mulii(d,N2);
     127           0 :   v = Fp_add(s,d,N);
     128           0 :   if (equali1(v)) return NULL;
     129           0 :   return mkvec2(g,v);
     130             : }
     131             : 
     132             : static GEN
     133     2961371 : _Fp_extgcd(void *data, GEN x, GEN y, int* smallop)
     134             : {
     135             :   GEN d,u,v,m;
     136     2961371 :   if (equali1(y))
     137             :   {
     138      539658 :     *smallop = 1;
     139      539658 :     return mkvec2(y,mkmat2(
     140             :           mkcol2(gen_1,Fp_neg(x,(GEN)data)),
     141             :           mkcol2(gen_0,gen_1)));
     142             :   }
     143     2421713 :   *smallop = 0;
     144     2421713 :   d = bezout(x,y,&u,&v);
     145     2421713 :   if (!signe(d)) return mkvec2(d,matid(2));
     146     2421713 :   m = cgetg(3,t_MAT);
     147     2421713 :   m = mkmat2(
     148             :     mkcol2(diviiexact(y,d),negi(diviiexact(x,d))),
     149             :     mkcol2(u,v));
     150     2421713 :   return mkvec2(d,m);
     151             : }
     152             : 
     153             : static int
     154    65470874 : _Fp_equal0(GEN x) { return !signe(x); }
     155             : 
     156             : static int
     157    19201812 : _Fp_equal1(GEN x) { return equali1(x); }
     158             : 
     159             : static GEN
     160     9938530 : _Fp_s(void *data, long x)
     161             : {
     162     9938530 :   if (!x) return gen_0;
     163      782481 :   if (x==1) return gen_1;
     164           0 :   return modsi(x,(GEN)data);
     165             : }
     166             : 
     167             : static GEN
     168    31555636 : _Fp_red(void *data, GEN x) { return Fp_red(x, (GEN)data); }
     169             : 
     170             : /* p not necessarily prime */
     171             : static const struct bb_hermite Fp_hermite=
     172             :   {_Fp_add,_Fp_neg,_Fp_mul,_Fp_extgcd,_Fp_rann,_Fp_lquo,_Fp_unit,_Fp_equal0,_Fp_equal1,_Fp_s,_Fp_red};
     173             : 
     174             : static const struct bb_hermite*
     175      284578 : get_Fp_hermite(void **data, GEN p)
     176             : {
     177      284578 :   *data = (void*)p; return &Fp_hermite;
     178             : }
     179             : 
     180             : static void
     181    11323564 : gen_redcol(GEN C, long lim, void* data, const struct bb_hermite *R)
     182             : {
     183             :   long i;
     184    43000769 :   for (i=1; i<=lim; i++)
     185    31677205 :     if (!R->equal0(gel(C,i)))
     186    19465523 :       gel(C,i) = R->red(data, gel(C,i));
     187    11323564 : }
     188             : 
     189             : static GEN
     190             : /* return NULL if a==0 */
     191             : /* assume C*a is zero after lim */
     192    15241730 : gen_rightmulcol(GEN C, GEN a, long lim, int fillzeros, void* data, const struct bb_hermite *R)
     193             : {
     194             :   GEN Ca,zero;
     195             :   long i;
     196    15241730 :   if (R->equal1(a)) return C;
     197     9759701 :   if (R->equal0(a)) return NULL;
     198     7112476 :   Ca = cgetg(lg(C),t_COL);
     199    26116013 :   for (i=1; i<=lim; i++)
     200    19003537 :     gel(Ca,i) = R->mul(data, gel(C,i), a);
     201     7112476 :   if (fillzeros && lim+1 < lg(C))
     202             :   {
     203     5436550 :     zero = R->s(data,0);
     204    23496893 :     for (i=lim+1; i<lg(C); i++)
     205    18060343 :       gel(Ca,i) = zero;
     206             :   }
     207     7112476 :   return Ca;
     208             : }
     209             : 
     210             : static void
     211             : /* C1 <- C1 + C2 */
     212             : /* assume C2[i]==0 for i>lim */
     213     4864433 : gen_addcol(GEN C1, GEN C2, long lim, void* data, const struct bb_hermite *R)
     214             : {
     215             :   long i;
     216    17365285 :   for (i=1; i<=lim; i++)
     217    12500852 :     gel(C1,i) = R->add(data, gel(C1,i), gel(C2,i));
     218     4864433 : }
     219             : 
     220             : static void
     221             : /* H[,i] <- H[,i] + C*a */
     222             : /* assume C is zero after lim */
     223     2285535 : gen_addrightmul(GEN H, GEN C, GEN a, long i, long lim, void* data, const struct bb_hermite *R)
     224             : {
     225             :   GEN Ca;
     226     2285535 :   if (R->equal0(a)) return;
     227     1300096 :   Ca = gen_rightmulcol(C, a, lim, 0, data, R);
     228     1300096 :   gen_addcol(gel(H,i), Ca, lim, data, R);
     229             : }
     230             : 
     231             : static GEN
     232      970165 : gen_zerocol(long n, void* data, const struct bb_hermite *R)
     233             : {
     234      970165 :   GEN C = cgetg(n+1,t_COL), zero = R->s(data, 0);
     235             :   long i;
     236      970165 :   for (i=1; i<=n; i++) gel(C,i) = zero;
     237      970165 :   return C;
     238             : }
     239             : 
     240             : static GEN
     241      441231 : gen_zeromat(long m, long n, void* data, const struct bb_hermite *R)
     242             : {
     243      441231 :   GEN M = cgetg(n+1,t_MAT);
     244             :   long i;
     245      441231 :   for (i=1; i<=n; i++) gel(M,i) = gen_zerocol(m, data, R);
     246      441231 :   return M;
     247             : }
     248             : 
     249             : static GEN
     250      364420 : gen_colei(long n, long i, void* data, const struct bb_hermite *R)
     251             : {
     252      364420 :   GEN C = cgetg(n+1,t_COL), zero = R->s(data, 0);
     253             :   long j;
     254     2501548 :   for (j=1; j<=n; j++)
     255     2137128 :     if (i!=j)   gel(C,j) = zero;
     256      364420 :     else        gel(C,j) = R->s(data,1);
     257      364420 :   return C;
     258             : }
     259             : 
     260             : static GEN
     261       56056 : gen_matid_hermite(long n, void* data, const struct bb_hermite *R)
     262             : {
     263       56056 :   GEN M = cgetg(n+1,t_MAT);
     264             :   long i;
     265       56056 :   for (i=1; i<=n; i++) gel(M,i) = gen_colei(n, i, data, R);
     266       56056 :   return M;
     267             : }
     268             : 
     269             : static GEN
     270     2001902 : gen_matmul_hermite(GEN A, GEN B, void* data, const struct bb_hermite *R)
     271             : {
     272     2001902 :   GEN M,sum,prod,zero = R->s(data,0);
     273             :   long a,b,c,c2,i,j,k;
     274     2001902 :   RgM_dimensions(A,&a,&c);
     275     2001902 :   RgM_dimensions(B,&c2,&b);
     276     2001902 :   if (c!=c2) pari_err_DIM("gen_matmul_hermite");
     277     2001902 :   M = cgetg(b+1,t_MAT);
     278     4331250 :   for (j=1; j<=b; j++)
     279             :   {
     280     2329348 :     gel(M,j) = cgetg(a+1,t_COL);
     281     6451424 :     for (i=1; i<=a; i++)
     282             :     {
     283     4122076 :       sum = zero;
     284    13312558 :       for (k=1; k<=c; k++)
     285             :       {
     286     9190482 :         prod = R->mul(data, gcoeff(A,i,k), gcoeff(B,k,j));
     287     9190482 :         sum = R->add(data, sum, prod);
     288             :       }
     289     4122076 :       gcoeff(M,i,j) = sum;
     290             :     }
     291     2329348 :     gen_redcol(gel(M,j), a, data, R);
     292             :   }
     293     2001902 :   return M;
     294             : }
     295             : 
     296             : static void
     297             : /* U = [u1,u2]~, C <- A*u1 + B*u2 */
     298             : /* assume both A, B and C are zero after lim */
     299     6100416 : gen_rightlincomb(GEN A, GEN B, GEN U, GEN *C, long lim, void* data, const struct bb_hermite *R)
     300             : {
     301             :   GEN Au1, Bu2;
     302     6100416 :   Au1 = gen_rightmulcol(A, gel(U,1), lim, 1, data, R);
     303     6100416 :   Bu2 = gen_rightmulcol(B, gel(U,2), lim, 1, data, R);
     304     6100416 :   if (!Au1 && !Bu2) { *C = gen_zerocol(lg(A)-1, data, R); return; }
     305     6100416 :   if (!Au1) { *C = Bu2; return; }
     306     3799831 :   if (!Bu2) { *C = Au1; return; }
     307     3564211 :   gen_addcol(Au1, Bu2, lim, data, R);
     308     3564211 :   *C = Au1;
     309             : }
     310             : 
     311             : static void
     312             : /* (H[,i] | H[,j]) <- (H[,i] | H[,j]) * U */
     313             : /* assume both columns are zero after lim */
     314     3050208 : gen_elem(GEN H, GEN U, long i, long j, long lim, void* data, const struct bb_hermite *R)
     315             : {
     316             :   GEN Hi, Hj;
     317     3050208 :   Hi = shallowcopy(gel(H,i));
     318     3050208 :   Hj = shallowcopy(gel(H,j));
     319     3050208 :   gen_rightlincomb(Hi, Hj, gel(U,1), &gel(H,i), lim, data, R);
     320     3050208 :   gen_rightlincomb(Hi, Hj, gel(U,2), &gel(H,j), lim, data, R);
     321     3050208 : }
     322             : 
     323             : static int
     324             : /* assume C is zero after lim */
     325     2101967 : gen_is_zerocol(GEN C, long lim, void* data, const struct bb_hermite *R)
     326             : {
     327             :   long i;
     328             :   (void) data;
     329     3852324 :   for (i=1; i<=lim; i++)
     330     2827468 :     if (!R->equal0(gel(C,i))) return 0;
     331     1024856 :   return 1;
     332             : }
     333             : 
     334             : /* The mkop* functions return NULL if the corresponding operation is the identity */
     335             : 
     336             : static GEN
     337             : /* Ci <- Ci + Cj*a */
     338     1203685 : mkoptransv(long i, long j, GEN a, void* data, const struct bb_hermite *R)
     339             : {
     340     1203685 :   a = R->red(data,a);
     341     1203685 :   if (R->equal0(a)) return NULL;
     342      973077 :   return mkvec2(mkvecsmall2(i,j),a);
     343             : }
     344             : 
     345             : static GEN
     346             : /* (Ci|Cj) <- (Ci|Cj)*U */
     347      748573 : mkopU(long i, long j, GEN U, void* data, const struct bb_hermite *R)
     348             : {
     349      748573 :   if (R->equal1(gcoeff(U,1,1)) && R->equal0(gcoeff(U,1,2))
     350      382165 :       && R->equal1(gcoeff(U,2,2))) return mkoptransv(i,j,gcoeff(U,2,1),data,R);
     351      366422 :   return mkvec2(mkvecsmall3(i,j,0),U);
     352             : }
     353             : 
     354             : static GEN
     355             : /* Ci <- Ci*u */
     356      463827 : mkopmul(long i, GEN u, const struct bb_hermite *R)
     357             : {
     358      463827 :   if (R->equal1(u)) return NULL;
     359      103117 :   return mkvec2(mkvecsmall(i),u);
     360             : }
     361             : 
     362             : static GEN
     363             : /* Ci <-> Cj */
     364       45724 : mkopswap(long i, long j)
     365             : {
     366       45724 :   return mkvec(mkvecsmall2(i,j));
     367             : }
     368             : 
     369             : /* M: t_MAT. Apply the operation op to M by right multiplication. */
     370             : static void
     371      374465 : gen_rightapply(GEN M, GEN op, void* data, const struct bb_hermite *R)
     372             : {
     373             :   GEN M2, ind, X;
     374      374465 :   long i, j, m = lg(gel(M,1))-1;
     375      374465 :   switch (typ(op))
     376             :   {
     377             :     case t_VECSMALL:
     378       56014 :       M2 = vecpermute(M,op);
     379       56014 :       for (i=1; i<lg(M); i++) gel(M,i) = gel(M2,i);
     380       56014 :       return;
     381             :     case t_VEC:
     382      318451 :       ind = gel(op,1);
     383      318451 :       switch (lg(op))
     384             :       {
     385             :         case 2:
     386       16135 :           swap(gel(M,ind[1]),gel(M,ind[2]));
     387       16135 :           return;
     388             :         case 3:
     389      302316 :           X = gel(op,2);
     390      302316 :           i = ind[1];
     391      302316 :           switch (lg(ind))
     392             :           {
     393             :             case 2:
     394       37786 :               gel(M,i) = gen_rightmulcol(gel(M,i), X, m, 0, data, R);
     395       37786 :               gen_redcol(gel(M,i), m, data, R);
     396       37786 :               return;
     397             :             case 3:
     398      175693 :               gen_addrightmul(M, gel(M,ind[2]), X, i, m, data, R);
     399      175693 :               gen_redcol(gel(M,i), m, data, R);
     400      175693 :               return;
     401             :             case 4:
     402       88837 :               j = ind[2];
     403       88837 :               gen_elem(M, X, i, j, m, data, R);
     404       88837 :               gen_redcol(gel(M,i), m, data, R);
     405       88837 :               gen_redcol(gel(M,j), m, data, R);
     406       88837 :               return;
     407             :           }
     408             :       }
     409             :   }
     410             : }
     411             : 
     412             : /* C: t_COL. Apply the operation op to C by left multiplication. */
     413             : static void
     414     9582496 : gen_leftapply(GEN C, GEN op, void* data, const struct bb_hermite *R)
     415             : {
     416             :   GEN C2, ind, X;
     417             :   long i, j;
     418     9582496 :   switch (typ(op))
     419             :   {
     420             :     case t_VECSMALL:
     421      575841 :       C2 = vecpermute(C,perm_inv(op));
     422      575841 :       for (i=1; i<lg(C); i++) gel(C,i) = gel(C2,i);
     423      575841 :       return;
     424             :     case t_VEC:
     425     9006655 :       ind = gel(op,1);
     426     9006655 :       switch (lg(op))
     427             :       {
     428             :         case 2:
     429      166390 :           swap(gel(C,ind[1]),gel(C,ind[2]));
     430      166390 :           return;
     431             :         case 3:
     432     8840265 :           X = gel(op,2);
     433     8840265 :           i = ind[1];
     434     8840265 :           switch (lg(ind))
     435             :           {
     436             :             case 2:
     437      321258 :               gel(C,i) = R->mul(data, X, gel(C,i));
     438      321258 :               gel(C,i) = R->red(data, gel(C,i));
     439      321258 :               return;
     440             :             case 3:
     441     6726279 :               j = ind[2];
     442     6726279 :               if (R->equal0(gel(C,i))) return;
     443      798077 :               gel(C,j) = R->add(data, gel(C,j), R->mul(data, X, gel(C,i)));
     444      798077 :               return;
     445             :             case 4:
     446     1792728 :               j = ind[2];
     447     1792728 :               C2 = gen_matmul_hermite(X, mkmat(mkcol2(gel(C,i),gel(C,j))), data, R);
     448     1792728 :               gel(C,i) = gcoeff(C2,1,1);
     449     1792728 :               gel(C,j) = gcoeff(C2,2,1);
     450     1792728 :               return;
     451             :           }
     452             :       }
     453             :   }
     454             : }
     455             : 
     456             : /* \prod_i det ops[i]. Only makes sense if R is commutative. */
     457             : static GEN
     458          42 : gen_detops(GEN ops, void* data, const struct bb_hermite *R)
     459             : {
     460          42 :   GEN d = R->s(data,1);
     461          42 :   long i, l = lg(ops);
     462         231 :   for (i = 1; i < l; i++)
     463             :   {
     464         189 :     GEN X, op = gel(ops,i);
     465         189 :     switch (typ(op))
     466             :     {
     467             :       case t_VECSMALL:
     468           0 :         if (perm_sign(op) < 0) d = R->neg(data,d);
     469           0 :         break;
     470             :       case t_VEC:
     471         189 :         switch (lg(op))
     472             :         {
     473             :           case 2:
     474           0 :             d = R->neg(data,d);
     475           0 :             break;
     476             :           case 3:
     477         189 :             X = gel(op,2);
     478         189 :             switch (lg(gel(op,1)))
     479             :             {
     480             :               case 2:
     481           0 :                  d = R->mul(data, d, X);
     482           0 :                  d = R->red(data, d);
     483           0 :                  break;
     484             :               case 4:
     485             :                {
     486         105 :                  GEN A = gcoeff(X,1,1), B = gcoeff(X,1,2);
     487         105 :                  GEN C = gcoeff(X,2,1), D = gcoeff(X,2,2);
     488         105 :                  GEN AD = R->mul(data,A,D);
     489         105 :                  GEN BC = R->mul(data,B,C);
     490         105 :                  d = R->mul(data, d, R->add(data, AD, R->neg(data,BC)));
     491         105 :                  d = R->red(data, d);
     492         105 :                  break;
     493             :                }
     494             :             }
     495         189 :             break;
     496             :         }
     497         189 :         break;
     498             :     }
     499             :   }
     500          42 :   return d;
     501             : }
     502             : 
     503             : static int
     504      209258 : gen_is_inv(GEN x, void* data, const struct bb_hermite *R)
     505             : {
     506      209258 :   GEN u = R->unit(data, x);
     507      209258 :   if (!u) return R->equal1(x);
     508       71834 :   return R->equal1(gel(u,1));
     509             : }
     510             : 
     511             : static long
     512      193375 : gen_last_inv_diago(GEN A, void* data, const struct bb_hermite *R)
     513             : {
     514             :   long i,m,n,j;
     515      193375 :   RgM_dimensions(A,&m,&n);
     516      223496 :   for (i=1,j=n-m+1; i<=m; i++,j++)
     517      209258 :     if (!gen_is_inv(gcoeff(A,i,j),data,R)) return i-1;
     518       14238 :   return m;
     519             : }
     520             : 
     521             : static GEN
     522             : /* remove_zerocols: 0 none, 1 until square, 2 all */
     523             : /* early abort: if not right-invertible, abort, return NULL, and set ops to the
     524             :  * non-invertible pivot */
     525      383012 : gen_howell_i(GEN A, long remove_zerocols, long permute_zerocols, long early_abort, long only_triangular, GEN* ops, void *data, const struct bb_hermite *R)
     526             : {
     527      383012 :   pari_sp av = avma, av1;
     528      383012 :   GEN H,U,piv=gen_0,u,q,a,perm,iszero,C,zero=R->s(data,0),d,g,r,op,one=R->s(data,1);
     529      383012 :   long m,n,i,j,s,si,i2,si2,nbz,lim,extra,maxop=0,nbop=0,lastinv=0;
     530             :   int smallop;
     531             : 
     532      383012 :   av1 = avma;
     533             : 
     534      383012 :   RgM_dimensions(A,&m,&n);
     535      383012 :   if (early_abort && n<m)
     536             :   {
     537       14000 :     if (ops) *ops = zero;
     538       14000 :     return NULL;
     539             :   }
     540      369012 :   if (n<m+1)
     541             :   {
     542      270690 :     extra = m+1-n;
     543      270690 :     H = shallowmatconcat(mkvec2(gen_zeromat(m,extra,data,R),A));
     544             :   }
     545             :   else
     546             :   {
     547       98322 :     extra = 0;
     548       98322 :     H = RgM_shallowcopy(A);
     549             :   }
     550      369012 :   RgM_dimensions(H,&m,&n);
     551      369012 :   s = n-m; /* shift */
     552             : 
     553      369012 :   if(ops)
     554             :   {
     555      175637 :     maxop = m*n + (m*(m+1))/2 + 1;
     556      175637 :     *ops = zerovec(maxop); /* filled with placeholders so gerepile can handle it */
     557             :   }
     558             : 
     559             :   /* put in triangular form */
     560     1597204 :   for (i=m,si=s+m; i>0 && si>extra; i--,si--) /* si = s+i */
     561             :   {
     562     1233106 :     if (R->red) gcoeff(H,i,si) = R->red(data, gcoeff(H,i,si));
     563             :     /* bottom-right diagonal */
     564     6253653 :     for (j = extra+1; j < si; j++)
     565             :     {
     566     5020547 :       if (R->red) gcoeff(H,i,j) = R->red(data, gcoeff(H,i,j));
     567     5020547 :       if (R->equal0(gcoeff(H,i,j))) continue;
     568     2818718 :       U = R->extgcd(data, gcoeff(H,i,j), gcoeff(H,i,si), &smallop);
     569     2818718 :       d = gel(U,1);
     570     2818718 :       U = gel(U,2);
     571     2818718 :       if (n>10)
     572             :       {
     573             :         /* normalize diagonal coefficient -> faster reductions on this row */
     574     1761459 :         u = R->unit(data, d);
     575     1761459 :         if (u)
     576             :         {
     577     1761459 :           g = gel(u,1);
     578     1761459 :           u = gel(u,2);
     579     1761459 :           gcoeff(U,1,2) = R->mul(data, gcoeff(U,1,2), u);
     580     1761459 :           gcoeff(U,2,2) = R->mul(data, gcoeff(U,2,2), u);
     581     1761459 :           d = g;
     582             :         }
     583             :       }
     584     2818718 :       gen_elem(H, U, j, si, i-1, data, R);
     585     2818718 :       if (ops)
     586             :       {
     587      715099 :         op =  mkopU(j,si,U,data,R);
     588      715099 :         if (op) { nbop++; gel(*ops, nbop) = op; }
     589             :       }
     590     2818718 :       gcoeff(H,i,j) = zero;
     591     2818718 :       gcoeff(H,i,si) = d;
     592     2818718 :       if (R->red && !smallop)
     593             :       {
     594     2288629 :         gen_redcol(gel(H,si), i-1, data, R);
     595     2288629 :         gen_redcol(gel(H,j), i-1, data, R);
     596             :       }
     597             :     }
     598             : 
     599             : 
     600     1233106 :     if (early_abort)
     601             :     {
     602       46417 :       d = gcoeff(H,i,si);
     603       46417 :       u = R->unit(data, d);
     604       46417 :       if (u) d = gel(u,1);
     605       46417 :       if (!R->equal1(d))
     606             :       {
     607        4914 :         if (ops) *ops = d;
     608        4914 :         return NULL;
     609             :       }
     610             :     }
     611             : 
     612     1228192 :     if (gc_needed(av,1))
     613             :     {
     614           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"gen_howell[1]. i=%ld",i);
     615           0 :       gerepileall(av1,ops?2:1,&H,ops);
     616             :     }
     617             :   }
     618             : 
     619      364098 :   if (!ops)
     620      193375 :     lastinv = gen_last_inv_diago(H, data, R);
     621             : 
     622             :   /* put in reduced Howell form */
     623      364098 :   if (!only_triangular)
     624             :   {
     625     1730078 :     for (i=m,si=s+m; i>0; i--,si--) /* si = s+i */
     626             :     {
     627             :       /* normalize diagonal coefficient */
     628     1366022 :       if (i<=lastinv) /* lastinv>0 => !ops */
     629       30121 :         gcoeff(H,i,si) = one;
     630             :       else
     631             :       {
     632     1335901 :         u = R->unit(data,gcoeff(H,i,si));
     633     1335901 :         if (u)
     634             :         {
     635     1046780 :           g = gel(u,1);
     636     1046780 :           u = gel(u,2);
     637     1046780 :           gel(H,si) = gen_rightmulcol(gel(H,si), u, i-1, 1, data, R);
     638     1046780 :           gcoeff(H,i,si) = g;
     639     1046780 :           if (R->red) gen_redcol(gel(H,si), i-1, data, R);
     640     1046780 :           if (ops)
     641             :           {
     642      463827 :             op = mkopmul(si,u,R);
     643      463827 :             if (op) { nbop++; gel(*ops,nbop) = op; }
     644             :           }
     645             :         }
     646      289121 :         else if (R->red) gcoeff(H,i,si) = R->red(data, gcoeff(H,i,si));
     647             :       }
     648     1366022 :       piv = gcoeff(H,i,si);
     649             : 
     650             :       /* reduce above diagonal */
     651     1366022 :       if (!R->equal0(piv))
     652             :       {
     653     1076901 :         C = gel(H,si);
     654     3215975 :         for (j=si+1; j<=n; j++)
     655             :         {
     656     2139074 :           if (i<=lastinv) /* lastinv>0 => !ops */
     657       29232 :             gcoeff(H,i,j) = zero;
     658             :           else
     659             :           {
     660     2109842 :             gcoeff(H,i,j) = R->red(data, gcoeff(H,i,j));
     661     2109842 :             if (R->equal1(piv)) { q = gcoeff(H,i,j); r = zero; }
     662     1549996 :             else                q = R->lquo(data, gcoeff(H,i,j), piv, &r);
     663     2109842 :             q = R->neg(data,q);
     664     2109842 :             gen_addrightmul(H, C, q, j, i-1, data, R);
     665     2109842 :             if (ops)
     666             :             {
     667      672595 :               op = mkoptransv(j,si,q,data,R);
     668      672595 :               if (op) { nbop++; gel(*ops,nbop) = op; }
     669             :             }
     670     2109842 :             gcoeff(H,i,j) = r;
     671             :           }
     672             :         }
     673             :       }
     674             : 
     675             :       /* ensure Howell property */
     676     1366022 :       if (i>1)
     677             :       {
     678     1002050 :         a = R->rann(data, piv);
     679     1002050 :         if (!R->equal0(a))
     680             :         {
     681      545090 :           gel(H,1) = gen_rightmulcol(gel(H,si), a, i-1, 1, data, R);
     682      545090 :           if (gel(H,1) == gel(H,si)) gel(H,1) = shallowcopy(gel(H,1)); /* in case rightmulcol cheated */
     683      545090 :           if (ops)
     684             :           {
     685      148939 :             op = mkoptransv(1,si,a,data,R);
     686      148939 :             if (op) { nbop++; gel(*ops,nbop) = op; }
     687             :           }
     688     2184336 :           for (i2=i-1,si2=s+i2; i2>0; i2--,si2--)
     689             :           {
     690     1639246 :             if (R->red) gcoeff(H,i2,1) = R->red(data, gcoeff(H,i2,1));
     691     1639246 :             if (R->equal0(gcoeff(H,i2,1))) continue;
     692      273203 :             if (R->red) gcoeff(H,i2,si2) = R->red(data, gcoeff(H,i2,si2));
     693      273203 :             if (R->equal0(gcoeff(H,i2,si2)))
     694             :             {
     695      130550 :               swap(gel(H,1), gel(H,si2));
     696      130550 :               if (ops) { nbop++; gel(*ops,nbop) = mkopswap(1,si2); }
     697      130550 :               continue;
     698             :             }
     699      142653 :             U = R->extgcd(data, gcoeff(H,i2,1), gcoeff(H,i2,si2), &smallop);
     700      142653 :             d = gel(U,1);
     701      142653 :             U = gel(U,2);
     702      142653 :             gen_elem(H, U, 1, si2, i2-1, data, R);
     703      142653 :             if (ops)
     704             :             {
     705       33474 :               op = mkopU(1,si2,U,data,R);
     706       33474 :               if (op) { nbop++; gel(*ops,nbop) = op; }
     707             :             }
     708      142653 :             gcoeff(H,i2,1) = zero;
     709      142653 :             gcoeff(H,i2,si2) = d;
     710      142653 :             if (R->red && !smallop)
     711             :             {
     712      133084 :               gen_redcol(gel(H,si2), i2, data, R);
     713      133084 :               gen_redcol(gel(H,1), i2-1, data, R);
     714             :             }
     715             :           }
     716             :         }
     717             :       }
     718             : 
     719     1366022 :       if (gc_needed(av,1))
     720             :       {
     721           0 :         if (DEBUGMEM>1) pari_warn(warnmem,"gen_howell[2]. i=%ld",i);
     722           0 :         gerepileall(av1,ops?3:2,&H,&piv,ops);
     723             :       }
     724             :     }
     725             :   }
     726             : 
     727      364098 :   if (R->red)
     728     2422140 :     for (j=1; j<=n; j++)
     729             :     {
     730     2058042 :       lim = maxss(0,m-n+j);
     731     2058042 :       gen_redcol(gel(H,j), lim, data, R);
     732     2058042 :       for (i=lim+1; i<=m; i++) gcoeff(H,i,j) = zero;
     733             :     }
     734             : 
     735             :   /* put zero columns first */
     736      364098 :   iszero = cgetg(n+1,t_VECSMALL);
     737             : 
     738      364098 :   nbz = 0;
     739     2422140 :   for (i=1; i<=n; i++)
     740             :   {
     741     2058042 :     iszero[i] = gen_is_zerocol(gel(H,i), maxss(0,m-n+i), data, R);
     742     2058042 :     if (iszero[i]) nbz++;
     743             :   }
     744             : 
     745      364098 :   j = 1;
     746      364098 :   if (permute_zerocols)
     747             :   {
     748      154588 :     perm = cgetg(n+1, t_VECSMALL);
     749      896826 :     for (i=1; i<=n; i++)
     750      742238 :       if (iszero[i])
     751             :       {
     752      313460 :         perm[j] = i;
     753      313460 :         j++;
     754             :       }
     755             :   }
     756      209510 :   else perm = cgetg(n-nbz+1, t_VECSMALL);
     757     2422140 :   for (i=1; i<=n; i++)
     758     2058042 :     if (!iszero[i])
     759             :     {
     760     1077062 :       perm[j] = i;
     761     1077062 :       j++;
     762             :     }
     763             : 
     764      364098 :   if (permute_zerocols || remove_zerocols==2) H = vecpermute(H, perm);
     765      364098 :   if (permute_zerocols && remove_zerocols==2) H = vecslice(H, nbz+1, n);
     766      364098 :   if (remove_zerocols==1) H = vecslice(H, s+1, n);
     767      364098 :   if (permute_zerocols && ops) { nbop++; gel(*ops,nbop) = perm; }
     768             : 
     769      364098 :   if (ops) { setlg(*ops, nbop+1); } /* should have nbop <= maxop */
     770             : 
     771      364098 :   return H;
     772             : }
     773             : 
     774             : static GEN
     775       94941 : gen_howell(GEN A, long remove_zerocols, long permute_zerocols, long early_abort, long only_triangular, GEN* ops, void *data, const struct bb_hermite *R)
     776             : {
     777       94941 :   pari_sp av = avma;
     778       94941 :   GEN H = gen_howell_i(A, remove_zerocols, permute_zerocols, early_abort, only_triangular, ops, data, R);
     779       94941 :   gerepileall(av, ops?2:1, &H, ops);
     780       94941 :   return H;
     781             : }
     782             : 
     783             : static GEN
     784      150955 : gen_matimage(GEN A, GEN* U, void *data, const struct bb_hermite *R)
     785             : {
     786             :   GEN ops, H;
     787      150955 :   if (U)
     788             :   {
     789       56014 :     pari_sp av = avma;
     790             :     long m, n, i, r, n2;
     791       56014 :     RgM_dimensions(A,&m,&n);
     792       56014 :     H = gen_howell_i(A, 2, 1, 0, 0, &ops, data, R);
     793       56014 :     r = lg(H)-1;
     794       56014 :     *U = shallowmatconcat(mkvec2(gen_zeromat(n, maxss(0,m-n+1), data, R), gen_matid_hermite(n, data, R)));
     795       56014 :     n2 = lg(*U)-1;
     796      430479 :     for (i=1; i<lg(ops); i++)
     797      374465 :       gen_rightapply(*U, gel(ops,i), data, R);
     798       56014 :     if (r<n2) *U = vecslice(*U, n2-r+1, n2);
     799       56014 :     gerepileall(av, 2, &H, U);
     800       56014 :     return H;
     801             :   }
     802       94941 :   else return gen_howell(A, 2, 0, 0, 0, NULL, data, R);
     803             : }
     804             : 
     805             : static GEN
     806             : /* H in true Howell form: no zero columns */
     807       98434 : gen_kernel_howell(GEN H, void *data, const struct bb_hermite *R)
     808             : {
     809             :   GEN K, piv, FK;
     810             :   long m, n, j, j2, i;
     811       98434 :   RgM_dimensions(H,&m,&n);
     812       98434 :   K = gen_zeromat(n, n, data, R);
     813      406042 :   for (j=n,i=m; j>0; j--)
     814             :   {
     815      307608 :     while (R->equal0(gcoeff(H,i,j))) i--;
     816      307608 :     piv = gcoeff(H,i,j);
     817      307608 :     if (R->equal0(piv)) continue;
     818      307608 :     gcoeff(K,j,j) = R->rann(data, piv);
     819      307608 :     if (j<n)
     820             :     {
     821      209174 :       FK = gen_matmul_hermite(matslice(H,i,i,j+1,n), matslice(K, j+1, n, j+1, n), data, R);
     822      745794 :       for (j2=j+1; j2<=n; j2++)
     823      536620 :         gcoeff(K,j,j2) = R->neg(data, R->lquo(data, gcoeff(FK,1,j2-j), piv, NULL));
     824             :         /* remainder has to be zero */
     825             :     }
     826             :   }
     827       98434 :   return K;
     828             : }
     829             : 
     830             : static GEN
     831             : /* (H,ops) Howell form of A, n = number of columns of A, return a kernel of A */
     832       98476 : gen_kernel_from_howell(GEN H, GEN ops, long n, void *data, const struct bb_hermite *R)
     833             : {
     834             :   pari_sp av;
     835             :   GEN K, KH, zC;
     836             :   long m, r, n2, nbz, i, o, extra, j;
     837       98476 :   RgM_dimensions(H,&m,&r);
     838       98476 :   if (!r) return gen_matid_hermite(n, data, R); /* zerology: what if 0==1 in R? */
     839       98434 :   n2 = maxss(n,m+1);
     840       98434 :   extra = n2-n;
     841       98434 :   nbz = n2-r;
     842             :   /* compute kernel of augmented matrix */
     843       98434 :   KH = gen_kernel_howell(H, data, R);
     844       98434 :   zC = gen_zerocol(nbz, data, R);
     845       98434 :   K = cgetg(nbz+r+1, t_MAT);
     846      322791 :   for (i=1; i<=nbz; i++)
     847      224357 :     gel(K,i) = gen_colei(nbz+r, i, data, R);
     848      406042 :   for (i=1; i<=r; i++)
     849      307608 :     gel(K,nbz+i) = shallowconcat(zC, gel(KH,i));
     850      630399 :   for (i=1; i<lg(K); i++)
     851             :   {
     852      531965 :     av = avma;
     853     9375107 :     for (o=lg(ops)-1; o>0; o--)
     854     8843142 :       gen_leftapply(gel(K,i), gel(ops,o), data, R);
     855      531965 :     gen_redcol(gel(K,i), nbz+r, data, R);
     856      531965 :     gerepileall(av, 1, &gel(K,i));
     857             :   }
     858             :   /* deduce kernel of original matrix */
     859       98434 :   K = rowpermute(K, cyclic_perm(n2,extra));
     860       98434 :   K = gen_howell_i(K, 2, 0, 0, 0, NULL, data, R);
     861      224840 :   for (j=lg(K)-1, i=n2; j>0; j--)
     862             :   {
     863      196392 :     while (R->equal0(gcoeff(K,i,j))) i--;
     864      196392 :     if (i<=n) return matslice(K, 1, n, 1, j);
     865             :   }
     866       28448 :   return cgetg(1,t_MAT);
     867             : }
     868             : 
     869             : /* not GC-clean */
     870             : static GEN
     871       54649 : gen_kernel(GEN A, GEN* im, void *data, const struct bb_hermite *R)
     872             : {
     873       54649 :   pari_sp av = avma;
     874       54649 :   long n = lg(A)-1;
     875             :   GEN H, ops, K;
     876       54649 :   H = gen_howell_i(A, 2, 1, 0, 0, &ops, data, R);
     877       54649 :   gerepileall(av,2,&H,&ops);
     878       54649 :   K = gen_kernel_from_howell(H, ops, n, data, R);
     879       54649 :   if (im) *im = H;
     880       54649 :   return K;
     881             : }
     882             : 
     883             : /* right inverse, not GC-clean */
     884             : static GEN
     885       35007 : gen_inv(GEN A, void* data, const struct bb_hermite *R)
     886             : {
     887             :   pari_sp av;
     888       35007 :   GEN ops, H, U, un=R->s(data,1);
     889             :   long m,n,j,o,n2;
     890       35007 :   RgM_dimensions(A,&m,&n);
     891       35007 :   av = avma;
     892       35007 :   H = gen_howell_i(A, 0, 0, 1, 0, &ops, data, R);
     893       35007 :   if (!H) pari_err_INV("gen_inv", ops);
     894       16093 :   n2 = lg(H)-1;
     895       16093 :   ops = gerepilecopy(av,ops); /* get rid of H */
     896       16093 :   U = gen_zeromat(n2, m, data, R);
     897       51142 :   for (j=1; j<=m; j++)
     898       35049 :     gcoeff(U,j+n2-m,j) = un;
     899       51142 :   for (j=1; j<=m; j++)
     900             :   {
     901       35049 :     av = avma;
     902      376243 :     for (o=lg(ops)-1; o>0; o--)
     903      341194 :       gen_leftapply(gel(U,j), gel(ops,o), data, R);
     904       35049 :     gen_redcol(gel(U,j), n2, data, R);
     905       35049 :     gerepileall(av, 1, &gel(U,j));
     906             :   }
     907       16093 :   if (n2>n) U = rowslice(U, n2-n+1, n2);
     908       16093 :   return U;
     909             : }
     910             : 
     911             : /*
     912             :   H true Howell form (no zero columns).
     913             :   Compute Z = Y - HX canonical representative of R^m mod H(R^n)
     914             : */
     915             : static GEN
     916       43925 : gen_reduce_mod_howell(GEN H, GEN Y, GEN *X, void* data, const struct bb_hermite *R)
     917             : {
     918       43925 :   pari_sp av = avma;
     919             :   long m,n,i,j;
     920             :   GEN Z, q, r, C;
     921       43925 :   RgM_dimensions(H,&m,&n);
     922       43925 :   if (X) *X = gen_zerocol(n,data,R);
     923       43925 :   Z = shallowcopy(Y);
     924       43925 :   i = m;
     925      155071 :   for (j=n; j>0; j--)
     926             :   {
     927      111146 :     while (R->equal0(gcoeff(H,i,j))) i--;
     928      111146 :     q = R->lquo(data, gel(Z,i), gcoeff(H,i,j), &r);
     929      111146 :     gel(Z,i) = r;
     930      111146 :     C = gen_rightmulcol(gel(H,j), R->neg(data,q), i, 0, data, R);
     931      111146 :     if (C) gen_addcol(Z, C, i-1, data, R);
     932      111146 :     if (X) gel(*X,j) = q;
     933             :   }
     934       43925 :   gen_redcol(Z, lg(Z)-1, data, R);
     935       43925 :   if (X)
     936             :   {
     937       43925 :     gerepileall(av, 2, &Z, X);
     938       43925 :     return Z;
     939             :   }
     940           0 :   return gerepilecopy(av, Z);
     941             : }
     942             : 
     943             : /* H: Howell form of A */
     944             : /* (m,n): dimensions of original matrix A */
     945             : /* return NULL if no solution */
     946             : static GEN
     947       43925 : gen_solve_from_howell(GEN H, GEN ops, long m, long n, GEN Y, void* data, const struct bb_hermite *R)
     948             : {
     949       43925 :   pari_sp av = avma;
     950             :   GEN Z, X;
     951             :   long n2, mH, nH, i;
     952       43925 :   RgM_dimensions(H,&mH,&nH);
     953       43925 :   n2 = maxss(n,m+1);
     954       43925 :   Z = gen_reduce_mod_howell(H, Y, &X, data, R);
     955       43925 :   if (!gen_is_zerocol(Z,m,data,R)) { set_avma(av); return NULL; }
     956             : 
     957       43876 :   X = shallowconcat(zerocol(n2-nH),X);
     958       43876 :   for (i=lg(ops)-1; i>0; i--) gen_leftapply(X, gel(ops,i), data, R);
     959       43876 :   X = vecslice(X, n2-n+1, n2);
     960       43876 :   gen_redcol(X, n, data, R);
     961       43876 :   return gerepilecopy(av, X);
     962             : }
     963             : 
     964             : /* return NULL if no solution, not GC-clean */
     965             : static GEN
     966       43925 : gen_solve(GEN A, GEN Y, GEN* K, void* data, const struct bb_hermite *R)
     967             : {
     968             :   GEN H, ops, X;
     969             :   long m,n;
     970             : 
     971       43925 :   RgM_dimensions(A,&m,&n);
     972       43925 :   if (!n) m = lg(Y)-1;
     973       43925 :   H = gen_howell_i(A, 2, 1, 0, 0, &ops, data, R);
     974       43925 :   X = gen_solve_from_howell(H, ops, m, n, Y, data, R);
     975       43925 :   if (!X) return NULL;
     976       43876 :   if (K) *K = gen_kernel_from_howell(H, ops, n, data, R);
     977       43876 :   return X;
     978             : }
     979             : 
     980             : GEN
     981      150990 : matimagemod(GEN A, GEN d, GEN* U)
     982             : {
     983             :   void *data;
     984             :   const struct bb_hermite* R;
     985      150990 :   if (typ(A)!=t_MAT || !RgM_is_ZM(A)) pari_err_TYPE("matimagemod", A);
     986      150976 :   if (typ(d)!=t_INT) pari_err_TYPE("matimagemod", d);
     987      150969 :   if (signe(d)<=0) pari_err_DOMAIN("matimagemod", "d", "<=", gen_0, d);
     988      150962 :   if (equali1(d)) return cgetg(1,t_MAT);
     989      150955 :   R = get_Fp_hermite(&data, d);
     990      150955 :   return gen_matimage(A, U, data, R);
     991             : }
     992             : 
     993             : /* for testing purpose */
     994             : /*
     995             : GEN
     996             : ZM_hnfmodid2(GEN A, GEN d)
     997             : {
     998             :   pari_sp av = avma;
     999             :   void *data;
    1000             :   long i;
    1001             :   const struct bb_hermite* R = get_Fp_hermite(&data, d);
    1002             :   GEN H;
    1003             :   if (typ(A)!=t_MAT || !RgM_is_ZM(A)) pari_err_TYPE("ZM_hnfmodid2", A);
    1004             :   if (typ(d)!=t_INT) pari_err_TYPE("ZM_hnfmodid2", d);
    1005             :   H = gen_howell_i(A, 1, 0, 0, 0, NULL, data, R);
    1006             :   for (i=1; i<lg(H); i++)
    1007             :     if (!signe(gcoeff(H,i,i))) gcoeff(H,i,i) = d;
    1008             :   return gerepilecopy(av,H);
    1009             : }
    1010             : */
    1011             : 
    1012             : GEN
    1013          84 : matdetmod(GEN A, GEN d)
    1014             : {
    1015          84 :   pari_sp av = avma;
    1016             :   void *data;
    1017             :   const struct bb_hermite* R;
    1018          84 :   long n = lg(A)-1, i;
    1019             :   GEN D, H, ops;
    1020          84 :   if (typ(A)!=t_MAT || !RgM_is_ZM(A)) pari_err_TYPE("matdetmod", A);
    1021          70 :   if (typ(d)!=t_INT) pari_err_TYPE("matdetmod", d);
    1022          70 :   if (signe(d)<=0) pari_err_DOMAIN("matdetmod", "d", "<=", gen_0, d);
    1023          63 :   if (!n) return equali1(d) ? gen_0 : gen_1;
    1024          56 :   if (n != nbrows(A)) pari_err_DIM("matdetmod");
    1025          49 :   if (equali1(d)) return gen_0;
    1026          42 :   R = get_Fp_hermite(&data, d);
    1027          42 :   H = gen_howell_i(A, 1, 0, 0, 1, &ops, data, R);
    1028          42 :   D = gen_detops(ops, data, R);
    1029          42 :   D = Fp_inv(D, d);
    1030          42 :   for (i = 1; i <= n; i++) D = Fp_mul(D, gcoeff(H,i,i), d);
    1031          42 :   return gerepileuptoint(av, D);
    1032             : }
    1033             : 
    1034             : GEN
    1035       54684 : matkermod(GEN A, GEN d, GEN* im)
    1036             : {
    1037       54684 :   pari_sp av = avma;
    1038             :   void *data;
    1039             :   const struct bb_hermite* R;
    1040             :   long m,n;
    1041             :   GEN K;
    1042       54684 :   if (typ(A)!=t_MAT || !RgM_is_ZM(A)) pari_err_TYPE("matkermod", A);
    1043       54670 :   if (typ(d)!=t_INT) pari_err_TYPE("matkermod", d);
    1044       54663 :   if (signe(d)<=0) pari_err_DOMAIN("makermod", "d", "<=", gen_0, d);
    1045       54656 :   if (equali1(d)) return cgetg(1,t_MAT);
    1046       54649 :   R = get_Fp_hermite(&data, d);
    1047       54649 :   RgM_dimensions(A,&m,&n);
    1048       54649 :   if (!im && m>2*n) /* TODO tune treshold */
    1049        3570 :     A = shallowtrans(matimagemod(shallowtrans(A),d,NULL));
    1050       54649 :   K = gen_kernel(A, im, data, R);
    1051       54649 :   gerepileall(av,im?2:1,&K,im);
    1052       54649 :   return K;
    1053             : }
    1054             : 
    1055             : /* left inverse */
    1056             : GEN
    1057       35056 : matinvmod(GEN A, GEN d)
    1058             : {
    1059       35056 :   pari_sp av = avma;
    1060             :   void *data;
    1061             :   const struct bb_hermite* R;
    1062             :   GEN U;
    1063       35056 :   if (typ(A)!=t_MAT || !RgM_is_ZM(A)) pari_err_TYPE("matinvmod", A);
    1064       35042 :   if (typ(d)!=t_INT) pari_err_TYPE("matinvmod", d);
    1065       35035 :   if (signe(d)<=0) pari_err_DOMAIN("matinvmod", "d", "<=", gen_0, d);
    1066       35028 :   if (equali1(d)) {
    1067             :     long m,n;
    1068          21 :     RgM_dimensions(A,&m,&n);
    1069          21 :     if (m<n) pari_err_INV("matinvmod",A);
    1070          14 :     return zeromatcopy(n,m);
    1071             :   }
    1072       35007 :   R = get_Fp_hermite(&data, d);
    1073       35007 :   U = gen_inv(shallowtrans(A), data, R);
    1074       16093 :   return gerepilecopy(av, shallowtrans(U));
    1075             : }
    1076             : 
    1077             : /* assume all D[i]>0, not GC-clean */
    1078             : static GEN
    1079       43925 : matsolvemod_finite(GEN M, GEN D, GEN Y, long flag)
    1080             : {
    1081             :   void *data;
    1082             :   const struct bb_hermite* R;
    1083             :   GEN X, K, d;
    1084             :   long m, n;
    1085             : 
    1086       43925 :   RgM_dimensions(M,&m,&n);
    1087       43925 :   if (typ(D)==t_COL)
    1088             :   { /* create extra variables for the D[i] */
    1089             :     GEN MD;
    1090          84 :     long i, i2, extra = 0;
    1091          84 :     d = gen_1;
    1092          84 :     for (i=1; i<lg(D); i++) d = lcmii(d,gel(D,i));
    1093          84 :     for (i=1; i<lg(D); i++) if (!equalii(gel(D,i),d)) extra++;
    1094          84 :     MD = cgetg(extra+1,t_MAT);
    1095          84 :     i2 = 1;
    1096         231 :     for (i=1; i<lg(D); i++)
    1097         147 :       if (!equalii(gel(D,i),d))
    1098             :       {
    1099          77 :         gel(MD,i2) = Rg_col_ei(gel(D,i),m,i);
    1100          77 :         i2++;
    1101             :       }
    1102          84 :     M = shallowconcat(M,MD);
    1103             :   }
    1104       43841 :   else d = D;
    1105             : 
    1106       43925 :   R = get_Fp_hermite(&data, d);
    1107       43925 :   X = gen_solve(M, Y, flag? &K: NULL, data, R);
    1108       43925 :   if (!X) return gen_0;
    1109       43876 :   X = vecslice(X,1,n);
    1110             : 
    1111       43876 :   if (flag)
    1112             :   {
    1113       43827 :     K = rowslice(K,1,n);
    1114       43827 :     K = hnfmodid(shallowconcat(zerocol(n),K),d);
    1115       43827 :     X = mkvec2(X,K);
    1116             :   }
    1117       43876 :   return X;
    1118             : }
    1119             : 
    1120             : /* Return a solution of congruence system sum M[i,j] X_j = Y_i mod D_i
    1121             :  * If pU != NULL, put in *pU a Z-basis of the homogeneous system.
    1122             :  * Works for all non negative D_i but inefficient compared to
    1123             :  * matsolvemod_finite; to be used only when one D_i is 0 */
    1124             : static GEN
    1125          70 : gaussmoduloall(GEN M, GEN D, GEN Y, GEN *pU)
    1126             : {
    1127          70 :   pari_sp av = avma;
    1128          70 :   long n, m, j, l, lM = lg(M);
    1129             :   GEN delta, H, U, u1, u2, x;
    1130             : 
    1131          70 :   if (lM == 1)
    1132             :   {
    1133          28 :     long lY = 0;
    1134          28 :     switch(typ(Y))
    1135             :     {
    1136           0 :       case t_INT: break;
    1137          28 :       case t_COL: lY = lg(Y); break;
    1138           0 :       default: pari_err_TYPE("gaussmodulo",Y);
    1139             :     }
    1140          28 :     switch(typ(D))
    1141             :     {
    1142          14 :       case t_INT: break;
    1143          14 :       case t_COL: if (lY && lY != lg(D)) pari_err_DIM("gaussmodulo");
    1144          14 :                   break;
    1145           0 :       default: pari_err_TYPE("gaussmodulo",D);
    1146             :     }
    1147          28 :     if (pU) *pU = cgetg(1, t_MAT);
    1148          28 :     return cgetg(1,t_COL);
    1149             :   }
    1150          42 :   n = nbrows(M);
    1151          42 :   switch(typ(D))
    1152             :   {
    1153             :     case t_COL:
    1154          14 :       if (lg(D)-1!=n) pari_err_DIM("gaussmodulo");
    1155          14 :       delta = diagonal_shallow(D); break;
    1156          28 :     case t_INT: delta = scalarmat_shallow(D,n); break;
    1157           0 :     default: pari_err_TYPE("gaussmodulo",D);
    1158             :       return NULL; /* LCOV_EXCL_LINE */
    1159             :   }
    1160          42 :   switch(typ(Y))
    1161             :   {
    1162           0 :     case t_INT: Y = const_col(n, Y); break;
    1163             :     case t_COL:
    1164          42 :       if (lg(Y)-1!=n) pari_err_DIM("gaussmodulo");
    1165          42 :       break;
    1166           0 :     default: pari_err_TYPE("gaussmodulo",Y);
    1167             :       return NULL; /* LCOV_EXCL_LINE */
    1168             :   }
    1169          42 :   H = ZM_hnfall_i(shallowconcat(M,delta), &U, 1);
    1170          42 :   Y = hnf_solve(H,Y); if (!Y) return gen_0;
    1171          35 :   l = lg(H); /* may be smaller than lM if some moduli are 0 */
    1172          35 :   n = l-1;
    1173          35 :   m = lg(U)-1 - n;
    1174          35 :   u1 = cgetg(m+1,t_MAT);
    1175          35 :   u2 = cgetg(n+1,t_MAT);
    1176          35 :   for (j=1; j<=m; j++) { GEN c = gel(U,j); setlg(c,lM); gel(u1,j) = c; }
    1177          35 :   U += m;
    1178          35 :   for (j=1; j<=n; j++) { GEN c = gel(U,j); setlg(c,lM); gel(u2,j) = c; }
    1179             :   /*       (u1 u2)
    1180             :    * (M D) (*  * ) = (0 H) */
    1181          35 :   u1 = ZM_lll(u1, 0.75, LLL_INPLACE);
    1182          35 :   Y = ZM_ZC_mul(u2,Y);
    1183          35 :   x = ZC_reducemodmatrix(Y, u1);
    1184          35 :   if (!pU) x = gerepileupto(av, x);
    1185             :   else
    1186             :   {
    1187          14 :     gerepileall(av, 2, &x, &u1);
    1188          14 :     *pU = u1;
    1189             :   }
    1190          35 :   return x;
    1191             : }
    1192             : /* to be used when one D_i is 0 */
    1193             : static GEN
    1194          70 : msolvemod0(GEN M, GEN D, GEN Y, long flag)
    1195             : {
    1196          70 :   pari_sp av = avma;
    1197             :   GEN y, x, U;
    1198          70 :   if (!flag) return gaussmoduloall(M,D,Y,NULL);
    1199          28 :   y = cgetg(3,t_VEC);
    1200          28 :   x = gaussmoduloall(M,D,Y,&U);
    1201          28 :   if (x == gen_0) { set_avma(av); return gen_0; }
    1202          21 :   gel(y,1) = x;
    1203          21 :   gel(y,2) = U; return y;
    1204             : 
    1205             : }
    1206             : GEN
    1207       44072 : matsolvemod(GEN M, GEN D, GEN Y, long flag)
    1208             : {
    1209       44072 :   pari_sp av = avma;
    1210       44072 :   long m, n, i, char0 = 0;
    1211       44072 :   if (typ(M)!=t_MAT) pari_err_TYPE("matsolvemod (M)",M);
    1212       44065 :   RgM_dimensions(M,&m,&n);
    1213       44065 :   if (typ(D)!=t_COL && typ(D)!=t_INT) pari_err_TYPE("matsolvemod (D)",D);
    1214       44051 :   if (n)
    1215       43918 :     { if (typ(D)==t_COL && lg(D)!=m+1) pari_err_DIM("matsolvemod [1]"); }
    1216             :   else
    1217         133 :     { if (typ(D)==t_COL) m = lg(D)-1; }
    1218       44044 :   if (typ(Y)==t_INT) Y = const_col(m,Y);
    1219       44044 :   if (typ(Y)!=t_COL) pari_err_TYPE("matsolvemod (Y)",Y);
    1220       44030 :   if (!n && !m) m = lg(Y)-1;
    1221       43974 :   else if (m != lg(Y)-1) pari_err_DIM("matsolvemod [2]");
    1222       44009 :   if (typ(D)==t_INT)
    1223             :   {
    1224       43890 :     if (signe(D)<0) pari_err_DOMAIN("matsolvemod","D","<",gen_0,D);
    1225       43883 :     if (!signe(D)) char0 = 1;
    1226             :   }
    1227             :   else /*typ(D)==t_COL*/
    1228         301 :     for (i=1; i<=m; i++)
    1229             :     {
    1230         189 :       if (signe(gel(D,i))<0)
    1231           7 :         pari_err_DOMAIN("matsolvemod","D[i]","<",gen_0,gel(D,i));
    1232         182 :       if (!signe(gel(D,i))) char0 = 1;
    1233             :     }
    1234       43995 :   return gerepilecopy(av, char0? msolvemod0(M,D,Y,flag)
    1235             :                                : matsolvemod_finite(M,D,Y,flag));
    1236             : }
    1237             : GEN
    1238           0 : gaussmodulo2(GEN M, GEN D, GEN Y) { return matsolvemod(M,D,Y, 1); }
    1239             : GEN
    1240           0 : gaussmodulo(GEN M, GEN D, GEN Y) { return matsolvemod(M,D,Y, 0); }

Generated by: LCOV version 1.13