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.18.0 lcov report (development 29806-4d001396c7) Lines: 643 664 96.8 %
Date: 2024-12-21 09:08:57 Functions: 53 55 96.4 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : #include "pari.h"
      15             : #include "paripriv.h"
      16             : 
      17             : #define DEBUGLEVEL DEBUGLEVEL_mathnf
      18             : 
      19             : #define dbg_printf(lvl) if (DEBUGLEVEL >= (lvl) + 3) err_printf
      20             : 
      21             : /********************************************************************/
      22             : /**                                                                **/
      23             : /**          BLACK BOX HERMITE RINGS AND HOWELL NORMAL FORM        **/
      24             : /**                 contributed by Aurel Page (2017)               **/
      25             : /**                                                                **/
      26             : /********************************************************************/
      27             : /* bb_hermite R:
      28             :  *  - add(a,b): a+b
      29             :  *  - neg(a): -a
      30             :  *  - mul(a,b): a*b
      31             :  *  - extgcd(a,b,&small): [d,U] with d in R and U in GL_2(R)
      32             :  *      such that [0;d] = [a;b]*U. Set small==1 to assert that U is a 'small'
      33             :  *      operation (no red needed).
      34             :  *  - rann(a): b in R such that b*R = {x in R | a*x==0}
      35             :  *  - lquo(a,b,&r): q in R such that r=a-b*q is a canonical representative
      36             :  *      of the image of a in R/b*R. The canonical lift of 0 must be 0.
      37             :  *  - unit(a): u unit in R^* such that a*u is a canonical generator of
      38             :  *      the ideal a*R
      39             :  *  - equal0(a): a==0?
      40             :  *  - equal1(a): a==1?
      41             :  *  - s(n): image of the small integer n in R
      42             :  *  - red(a): unique representative of a as an element of R
      43             : 
      44             :  * op encoding of elementary operations:
      45             :  *  - t_VECSMALL: the corresponding permutation (vecpermute)
      46             :  *  - [Vecsmall([i,j])]: the transposition Ci <-> Cj
      47             :  *  - [Vecsmall([i]),u], u in R^*: Ci <- Ci*u
      48             :  *  - [Vecsmall([i,j]),a], a in R: Ci <- Ci + Cj*a
      49             :  *  - [Vecsmall([i,j,0]),U], U in GL_2(R): (Ci|Cj) <- (Ci|Cj)*U */
      50             : 
      51             : struct bb_hermite
      52             : {
      53             :   GEN (*add)(void*, GEN, GEN);
      54             :   GEN (*neg)(void*, GEN);
      55             :   GEN (*mul)(void*, GEN, GEN);
      56             :   GEN (*extgcd)(void*, GEN, GEN, int*);
      57             :   GEN (*rann)(void*, GEN);
      58             :   GEN (*lquo)(void*, GEN, GEN, GEN*);
      59             :   GEN (*unit)(void*, GEN);
      60             :   int (*equal0)(GEN);
      61             :   int (*equal1)(GEN);
      62             :   GEN (*s)(void*, long);
      63             :   GEN (*red)(void*, GEN);
      64             : };
      65             : 
      66             : static GEN
      67    39768273 : _Fp_add(void *data, GEN x, GEN y) { (void) data; return addii(x,y); }
      68             : 
      69             : static GEN
      70     4710391 : _Fp_neg(void *data, GEN x) { (void) data; return negi(x); }
      71             : 
      72             : static GEN
      73    47497050 : _Fp_mul(void *data, GEN x, GEN y) { (void) data; return mulii(x,y); }
      74             : 
      75             : static GEN
      76     2225820 : _Fp_rann(void *data, GEN x)
      77             : {
      78     2225820 :   GEN d, N = (GEN)data;
      79     2225820 :   if (!signe(x)) return gen_1;
      80     2058263 :   d = gcdii(x,N);
      81     2058082 :   return modii(diviiexact(N,d),N);
      82             : }
      83             : 
      84             : static GEN
      85     2388242 : _Fp_lquo(void *data, GEN x, GEN y, GEN* r) { (void) data; return truedvmdii(x,y,r); }
      86             : 
      87             : /* D=MN, p|M => !p|a, p|N => p|a, return M */
      88             : static GEN
      89          28 : Z_split(GEN D, GEN a)
      90             : {
      91             :   long i, n;
      92             :   GEN N;
      93          28 :   n = expi(D);
      94          28 :   n = n<2 ? 1 : expu(n)+1;
      95         196 :   for (i=1;i<=n;i++)
      96         168 :     a = Fp_sqr(a,D);
      97          28 :   N = gcdii(a,D);
      98          28 :   return diviiexact(D,N);
      99             : }
     100             : 
     101             : /* c s.t. gcd(a+cb,N) = gcd(a,b,N) without factoring */
     102             : static GEN
     103          28 : Z_stab(GEN a, GEN b, GEN N)
     104             : {
     105             :   GEN g, a2, N2;
     106          28 :   g = gcdii(a,b);
     107          28 :   g = gcdii(g,N);
     108          28 :   N2 = diviiexact(N,g);
     109          28 :   a2 = diviiexact(a,g);
     110          28 :   return Z_split(N2,a2);
     111             : }
     112             : 
     113             : static GEN
     114     6931354 : _Fp_unit(void *data, GEN x)
     115             : {
     116     6931354 :   GEN g,s,v,d,N=(GEN)data,N2;
     117             :   long i;
     118     6931354 :   if (!signe(x)) return NULL;
     119     6469174 :   g = bezout(x,N,&s,&v);
     120     6470047 :   if (equali1(g) || equali1(gcdii(s,N))) return mkvec2(g,s);
     121       44743 :   N2 = diviiexact(N,g);
     122       61411 :   for (i=0; i<5; i++)
     123             :   {
     124       61383 :     s = addii(s,N2);
     125       61383 :     if (equali1(gcdii(s,N))) return mkvec2(g,s);
     126             :   }
     127          28 :   d = Z_stab(s,N2,N);
     128          28 :   d = mulii(d,N2);
     129          28 :   v = Fp_add(s,d,N);
     130          28 :   if (equali1(v)) return NULL;
     131          28 :   return mkvec2(g,v);
     132             : }
     133             : 
     134             : static GEN
     135     3828529 : _Fp_extgcd(void *data, GEN x, GEN y, int* smallop)
     136             : {
     137             :   GEN d,u,v,m;
     138     3828529 :   if (equali1(y))
     139             :   {
     140     1176010 :     *smallop = 1;
     141     1176010 :     return mkvec2(y,mkmat2(
     142             :           mkcol2(gen_1,Fp_neg(x,(GEN)data)),
     143             :           mkcol2(gen_0,gen_1)));
     144             :   }
     145     2652519 :   *smallop = 0;
     146     2652519 :   d = bezout(x,y,&u,&v);
     147     2652517 :   if (!signe(d)) return mkvec2(d,matid(2));
     148     2652517 :   m = cgetg(3,t_MAT);
     149     2652515 :   m = mkmat2(
     150             :     mkcol2(diviiexact(y,d),negi(diviiexact(x,d))),
     151             :     mkcol2(u,v));
     152     2652521 :   return mkvec2(d,m);
     153             : }
     154             : 
     155             : static int
     156   102658170 : _Fp_equal0(GEN x) { return !signe(x); }
     157             : 
     158             : static int
     159    29258513 : _Fp_equal1(GEN x) { return equali1(x); }
     160             : 
     161             : static GEN
     162    15084703 : _Fp_s(void *data, long x)
     163             : {
     164    15084703 :   if (!x) return gen_0;
     165     1912720 :   if (x==1) return gen_1;
     166           0 :   return modsi(x,(GEN)data);
     167             : }
     168             : 
     169             : static GEN
     170    47331955 : _Fp_red(void *data, GEN x) { return Fp_red(x, (GEN)data); }
     171             : 
     172             : /* p not necessarily prime */
     173             : static const struct bb_hermite Fp_hermite=
     174             :   {_Fp_add,_Fp_neg,_Fp_mul,_Fp_extgcd,_Fp_rann,_Fp_lquo,_Fp_unit,_Fp_equal0,_Fp_equal1,_Fp_s,_Fp_red};
     175             : 
     176             : static const struct bb_hermite*
     177      849688 : get_Fp_hermite(void **data, GEN p)
     178             : {
     179      849688 :   *data = (void*)p; return &Fp_hermite;
     180             : }
     181             : 
     182             : static void
     183    17578468 : gen_redcol(GEN C, long lim, void* data, const struct bb_hermite *R)
     184             : {
     185             :   long i;
     186    65308891 :   for (i=1; i<=lim; i++)
     187    47733177 :     if (!R->equal0(gel(C,i)))
     188    24493081 :       gel(C,i) = R->red(data, gel(C,i));
     189    17575714 : }
     190             : 
     191             : static GEN
     192             : /* return NULL if a==0 */
     193             : /* assume C*a is zero after lim */
     194    20349215 : gen_rightmulcol(GEN C, GEN a, long lim, int fillzeros, void* data, const struct bb_hermite *R)
     195             : {
     196             :   GEN Ca,zero;
     197             :   long i;
     198    20349215 :   if (R->equal1(a)) return C;
     199    11731612 :   if (R->equal0(a)) return NULL;
     200     8114326 :   Ca = cgetg(lg(C),t_COL);
     201    34523666 :   for (i=1; i<=lim; i++)
     202    26409347 :     gel(Ca,i) = R->mul(data, gel(C,i), a);
     203     8114319 :   if (fillzeros && lim+1 < lg(C))
     204             :   {
     205     6265375 :     zero = R->s(data,0);
     206    29376907 :     for (i=lim+1; i<lg(C); i++)
     207    23111535 :       gel(Ca,i) = zero;
     208             :   }
     209     8114316 :   return Ca;
     210             : }
     211             : 
     212             : static void
     213             : /* C1 <- C1 + C2 */
     214             : /* assume C2[i]==0 for i>lim */
     215     5801643 : gen_addcol(GEN C1, GEN C2, long lim, void* data, const struct bb_hermite *R)
     216             : {
     217             :   long i;
     218    30081437 :   for (i=1; i<=lim; i++)
     219    24279795 :     gel(C1,i) = R->add(data, gel(C1,i), gel(C2,i));
     220     5801642 : }
     221             : 
     222             : static void
     223             : /* H[,i] <- H[,i] + C*a */
     224             : /* assume C is zero after lim */
     225     4138427 : gen_addrightmul(GEN H, GEN C, GEN a, long i, long lim, void* data, const struct bb_hermite *R)
     226             : {
     227             :   GEN Ca;
     228     4138427 :   if (R->equal0(a)) return;
     229     1473110 :   Ca = gen_rightmulcol(C, a, lim, 0, data, R);
     230     1473112 :   gen_addcol(gel(H,i), Ca, lim, data, R);
     231             : }
     232             : 
     233             : static GEN
     234     2952905 : gen_zerocol(long n, void* data, const struct bb_hermite *R)
     235             : {
     236     2952905 :   GEN C = cgetg(n+1,t_COL), zero = R->s(data, 0);
     237             :   long i;
     238    14227040 :   for (i=1; i<=n; i++) gel(C,i) = zero;
     239     2952874 :   return C;
     240             : }
     241             : 
     242             : static GEN
     243     1562089 : gen_zeromat(long m, long n, void* data, const struct bb_hermite *R)
     244             : {
     245     1562089 :   GEN M = cgetg(n+1,t_MAT);
     246             :   long i;
     247     4371079 :   for (i=1; i<=n; i++) gel(M,i) = gen_zerocol(m, data, R);
     248     1562034 :   return M;
     249             : }
     250             : 
     251             : static GEN
     252      371371 : gen_colei(long n, long i, void* data, const struct bb_hermite *R)
     253             : {
     254      371371 :   GEN C = cgetg(n+1,t_COL), zero = R->s(data, 0);
     255             :   long j;
     256     2881291 :   for (j=1; j<=n; j++)
     257     2509920 :     if (i!=j)   gel(C,j) = zero;
     258      371371 :     else        gel(C,j) = R->s(data,1);
     259      371371 :   return C;
     260             : }
     261             : 
     262             : static GEN
     263       56084 : gen_matid_hermite(long n, void* data, const struct bb_hermite *R)
     264             : {
     265       56084 :   GEN M = cgetg(n+1,t_MAT);
     266             :   long i;
     267      196196 :   for (i=1; i<=n; i++) gel(M,i) = gen_colei(n, i, data, R);
     268       56084 :   return M;
     269             : }
     270             : 
     271             : static GEN
     272     2633082 : gen_matmul_hermite(GEN A, GEN B, void* data, const struct bb_hermite *R)
     273             : {
     274     2633082 :   GEN M,sum,prod,zero = R->s(data,0);
     275             :   long a,b,c,c2,i,j,k;
     276     2633082 :   RgM_dimensions(A,&a,&c);
     277     2633083 :   RgM_dimensions(B,&c2,&b);
     278     2633085 :   if (c!=c2) pari_err_DIM("gen_matmul_hermite");
     279     2633085 :   M = cgetg(b+1,t_MAT);
     280     5686620 :   for (j=1; j<=b; j++)
     281             :   {
     282     3053568 :     gel(M,j) = cgetg(a+1,t_COL);
     283     8524155 :     for (i=1; i<=a; i++)
     284             :     {
     285     5470661 :       sum = zero;
     286    19840624 :       for (k=1; k<=c; k++)
     287             :       {
     288    14370035 :         prod = R->mul(data, gcoeff(A,i,k), gcoeff(B,k,j));
     289    14369875 :         sum = R->add(data, sum, prod);
     290             :       }
     291     5470589 :       gcoeff(M,i,j) = sum;
     292             :     }
     293     3053494 :     gen_redcol(gel(M,j), a, data, R);
     294             :   }
     295     2633052 :   return M;
     296             : }
     297             : 
     298             : static void
     299             : /* U = [u1,u2]~, C <- A*u1 + B*u2 */
     300             : /* assume both A, B and C are zero after lim */
     301     7834725 : gen_rightlincomb(GEN A, GEN B, GEN U, GEN *C, long lim, void* data, const struct bb_hermite *R)
     302             : {
     303             :   GEN Au1, Bu2;
     304     7834725 :   Au1 = gen_rightmulcol(A, gel(U,1), lim, 1, data, R);
     305     7834731 :   Bu2 = gen_rightmulcol(B, gel(U,2), lim, 1, data, R);
     306     7834728 :   if (!Au1 && !Bu2) { *C = gen_zerocol(lg(A)-1, data, R); return; }
     307     7834728 :   if (!Au1) { *C = Bu2; return; }
     308     4712658 :   if (!Bu2) { *C = Au1; return; }
     309     4328407 :   gen_addcol(Au1, Bu2, lim, data, R);
     310     4328408 :   *C = Au1;
     311             : }
     312             : 
     313             : static void
     314             : /* (H[,i] | H[,j]) <- (H[,i] | H[,j]) * U */
     315             : /* assume both columns are zero after lim */
     316     3917368 : gen_elem(GEN H, GEN U, long i, long j, long lim, void* data, const struct bb_hermite *R)
     317             : {
     318             :   GEN Hi, Hj;
     319     3917368 :   Hi = shallowcopy(gel(H,i));
     320     3917369 :   Hj = shallowcopy(gel(H,j));
     321     3917369 :   gen_rightlincomb(Hi, Hj, gel(U,1), &gel(H,i), lim, data, R);
     322     3917362 :   gen_rightlincomb(Hi, Hj, gel(U,2), &gel(H,j), lim, data, R);
     323     3917364 : }
     324             : 
     325             : static int
     326             : /* assume C is zero after lim */
     327     4305294 : gen_is_zerocol(GEN C, long lim, void* data, const struct bb_hermite *R)
     328             : {
     329             :   long i;
     330             :   (void) data;
     331     7909133 :   for (i=1; i<=lim; i++)
     332     6127499 :     if (!R->equal0(gel(C,i))) return 0;
     333     1781634 :   return 1;
     334             : }
     335             : 
     336             : /* The mkop* functions return NULL if the corresponding operation is the identity */
     337             : 
     338             : static GEN
     339             : /* Ci <- Ci + Cj*a */
     340     2890631 : mkoptransv(long i, long j, GEN a, void* data, const struct bb_hermite *R)
     341             : {
     342     2890631 :   a = R->red(data,a);
     343     2890614 :   if (R->equal0(a)) return NULL;
     344     1124804 :   return mkvec2(mkvecsmall2(i,j),a);
     345             : }
     346             : 
     347             : static GEN
     348             : /* (Ci|Cj) <- (Ci|Cj)*U */
     349      955430 : mkopU(long i, long j, GEN U, void* data, const struct bb_hermite *R)
     350             : {
     351      955430 :   if (R->equal1(gcoeff(U,1,1)) && R->equal0(gcoeff(U,1,2))
     352      432900 :       && R->equal1(gcoeff(U,2,2))) return mkoptransv(i,j,gcoeff(U,2,1),data,R);
     353      522530 :   return mkvec2(mkvecsmall3(i,j,0),U);
     354             : }
     355             : 
     356             : static GEN
     357             : /* Ci <- Ci*u */
     358     1882264 : mkopmul(long i, GEN u, const struct bb_hermite *R)
     359             : {
     360     1882264 :   if (R->equal1(u)) return NULL;
     361      191607 :   return mkvec2(mkvecsmall(i),u);
     362             : }
     363             : 
     364             : static GEN
     365             : /* Ci <-> Cj */
     366       45843 : mkopswap(long i, long j)
     367             : {
     368       45843 :   return mkvec(mkvecsmall2(i,j));
     369             : }
     370             : 
     371             : /* M: t_MAT. Apply the operation op to M by right multiplication. */
     372             : static void
     373      374465 : gen_rightapply(GEN M, GEN op, void* data, const struct bb_hermite *R)
     374             : {
     375             :   GEN M2, ind, X;
     376      374465 :   long i, j, m = lg(gel(M,1))-1;
     377      374465 :   switch (typ(op))
     378             :   {
     379       56014 :     case t_VECSMALL:
     380       56014 :       M2 = vecpermute(M,op);
     381      266056 :       for (i=1; i<lg(M); i++) gel(M,i) = gel(M2,i);
     382       56014 :       return;
     383      318451 :     case t_VEC:
     384      318451 :       ind = gel(op,1);
     385      318451 :       switch (lg(op))
     386             :       {
     387       16135 :         case 2:
     388       16135 :           swap(gel(M,ind[1]),gel(M,ind[2]));
     389       16135 :           return;
     390      302316 :         case 3:
     391      302316 :           X = gel(op,2);
     392      302316 :           i = ind[1];
     393      302316 :           switch (lg(ind))
     394             :           {
     395       37786 :             case 2:
     396       37786 :               gel(M,i) = gen_rightmulcol(gel(M,i), X, m, 0, data, R);
     397       37786 :               gen_redcol(gel(M,i), m, data, R);
     398       37786 :               return;
     399      175693 :             case 3:
     400      175693 :               gen_addrightmul(M, gel(M,ind[2]), X, i, m, data, R);
     401      175693 :               gen_redcol(gel(M,i), m, data, R);
     402      175693 :               return;
     403       88837 :             case 4:
     404       88837 :               j = ind[2];
     405       88837 :               gen_elem(M, X, i, j, m, data, R);
     406       88837 :               gen_redcol(gel(M,i), m, data, R);
     407       88837 :               gen_redcol(gel(M,j), m, data, R);
     408       88837 :               return;
     409             :           }
     410             :       }
     411             :   }
     412             : }
     413             : 
     414             : /* C: t_COL. Apply the operation op to C by left multiplication. */
     415             : static void
     416    15152832 : gen_leftapply(GEN C, GEN op, void* data, const struct bb_hermite *R)
     417             : {
     418             :   GEN C2, ind, X;
     419             :   long i, j;
     420    15152832 :   switch (typ(op))
     421             :   {
     422      591024 :     case t_VECSMALL:
     423      591024 :       C2 = vecpermute(C,perm_inv(op));
     424     5503743 :       for (i=1; i<lg(C); i++) gel(C,i) = gel(C2,i);
     425      591024 :       return;
     426    14561808 :     case t_VEC:
     427    14561808 :       ind = gel(op,1);
     428    14561808 :       switch (lg(op))
     429             :       {
     430      169134 :         case 2:
     431      169134 :           swap(gel(C,ind[1]),gel(C,ind[2]));
     432      169134 :           return;
     433    14392675 :         case 3:
     434    14392675 :           X = gel(op,2);
     435    14392675 :           i = ind[1];
     436    14392675 :           switch (lg(ind))
     437             :           {
     438      701424 :             case 2:
     439      701424 :               gel(C,i) = R->mul(data, X, gel(C,i));
     440      701419 :               gel(C,i) = R->red(data, gel(C,i));
     441      701421 :               return;
     442    11274152 :             case 3:
     443    11274152 :               j = ind[2];
     444    11274152 :               if (R->equal0(gel(C,i))) return;
     445     1118575 :               gel(C,j) = R->add(data, gel(C,j), R->mul(data, X, gel(C,i)));
     446     1118576 :               return;
     447     2417123 :             case 4:
     448     2417123 :               j = ind[2];
     449     2417123 :               C2 = gen_matmul_hermite(X, mkmat(mkcol2(gel(C,i),gel(C,j))), data, R);
     450     2417099 :               gel(C,i) = gcoeff(C2,1,1);
     451     2417099 :               gel(C,j) = gcoeff(C2,2,1);
     452     2417099 :               return;
     453             :           }
     454             :       }
     455             :   }
     456             : }
     457             : 
     458             : /* \prod_i det ops[i]. Only makes sense if R is commutative. */
     459             : static GEN
     460          42 : gen_detops(GEN ops, void* data, const struct bb_hermite *R)
     461             : {
     462          42 :   GEN d = R->s(data,1);
     463          42 :   long i, l = lg(ops);
     464         231 :   for (i = 1; i < l; i++)
     465             :   {
     466         189 :     GEN X, op = gel(ops,i);
     467         189 :     switch (typ(op))
     468             :     {
     469           0 :       case t_VECSMALL:
     470           0 :         if (perm_sign(op) < 0) d = R->neg(data,d);
     471           0 :         break;
     472         189 :       case t_VEC:
     473         189 :         switch (lg(op))
     474             :         {
     475           0 :           case 2:
     476           0 :             d = R->neg(data,d);
     477           0 :             break;
     478         189 :           case 3:
     479         189 :             X = gel(op,2);
     480         189 :             switch (lg(gel(op,1)))
     481             :             {
     482           0 :               case 2:
     483           0 :                  d = R->mul(data, d, X);
     484           0 :                  d = R->red(data, d);
     485           0 :                  break;
     486         105 :               case 4:
     487             :                {
     488         105 :                  GEN A = gcoeff(X,1,1), B = gcoeff(X,1,2);
     489         105 :                  GEN C = gcoeff(X,2,1), D = gcoeff(X,2,2);
     490         105 :                  GEN AD = R->mul(data,A,D);
     491         105 :                  GEN BC = R->mul(data,B,C);
     492         105 :                  d = R->mul(data, d, R->add(data, AD, R->neg(data,BC)));
     493         105 :                  d = R->red(data, d);
     494         105 :                  break;
     495             :                }
     496             :             }
     497         189 :             break;
     498             :         }
     499         189 :         break;
     500             :     }
     501             :   }
     502          42 :   return d;
     503             : }
     504             : 
     505             : static int
     506      220689 : gen_is_inv(GEN x, void* data, const struct bb_hermite *R)
     507             : {
     508      220689 :   GEN u = R->unit(data, x);
     509      220689 :   if (!u) return R->equal1(x);
     510       75698 :   return R->equal1(gel(u,1));
     511             : }
     512             : 
     513             : static long
     514      201691 : gen_last_inv_diago(GEN A, void* data, const struct bb_hermite *R)
     515             : {
     516             :   long i,m,n,j;
     517      201691 :   RgM_dimensions(A,&m,&n);
     518      235669 :   for (i=1,j=n-m+1; i<=m; i++,j++)
     519      220689 :     if (!gen_is_inv(gcoeff(A,i,j),data,R)) return i-1;
     520       14980 :   return m;
     521             : }
     522             : 
     523             : static GEN
     524             : /* remove_zerocols: 0 none, 1 until square, 2 all */
     525             : /* early abort: if not right-invertible, abort, return NULL, and set ops to the
     526             :  * noninvertible pivot */
     527      949637 : 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)
     528             : {
     529      949637 :   pari_sp av = avma, av1;
     530      949637 :   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);
     531      949634 :   long m,n,i,j,s,si,i2,si2,nbz,lim,extra,maxop=0,nbop=0,lastinv=0;
     532             :   int smallop;
     533             : 
     534      949634 :   av1 = avma;
     535             : 
     536      949634 :   RgM_dimensions(A,&m,&n);
     537      949640 :   if (early_abort && n<m)
     538             :   {
     539       14000 :     if (ops) *ops = zero;
     540       14000 :     return NULL;
     541             :   }
     542      935640 :   if (n<m+1)
     543             :   {
     544      833258 :     extra = m+1-n;
     545      833258 :     H = shallowmatconcat(mkvec2(gen_zeromat(m,extra,data,R),A));
     546             :   }
     547             :   else
     548             :   {
     549      102382 :     extra = 0;
     550      102382 :     H = RgM_shallowcopy(A);
     551             :   }
     552      935655 :   RgM_dimensions(H,&m,&n);
     553      935656 :   s = n-m; /* shift */
     554             : 
     555      935656 :   if(ops)
     556             :   {
     557      733964 :     maxop = m*n + (m*(m+1))/2 + 1;
     558      733964 :     *ops = zerovec(maxop); /* filled with placeholders so gerepile can handle it */
     559             :   }
     560             : 
     561             :   /* put in triangular form */
     562     3638164 :   for (i=m,si=s+m; i>0 && si>extra; i--,si--) /* si = s+i */
     563             :   {
     564     2707512 :     if (R->red) gcoeff(H,i,si) = R->red(data, gcoeff(H,i,si));
     565             :     /* bottom-right diagonal */
     566    12839781 :     for (j = extra+1; j < si; j++)
     567             :     {
     568    10132640 :       if (R->red) gcoeff(H,i,j) = R->red(data, gcoeff(H,i,j));
     569    10132582 :       if (R->equal0(gcoeff(H,i,j))) continue;
     570     3674725 :       U = R->extgcd(data, gcoeff(H,i,j), gcoeff(H,i,si), &smallop);
     571     3674796 :       d = gel(U,1);
     572     3674796 :       U = gel(U,2);
     573     3674796 :       if (n>10)
     574             :       {
     575             :         /* normalize diagonal coefficient -> faster reductions on this row */
     576     2448766 :         u = R->unit(data, d);
     577     2448766 :         if (u)
     578             :         {
     579     2448766 :           g = gel(u,1);
     580     2448766 :           u = gel(u,2);
     581     2448766 :           gcoeff(U,1,2) = R->mul(data, gcoeff(U,1,2), u);
     582     2448766 :           gcoeff(U,2,2) = R->mul(data, gcoeff(U,2,2), u);
     583     2448766 :           d = g;
     584             :         }
     585             :       }
     586     3674796 :       gen_elem(H, U, j, si, i-1, data, R);
     587     3674792 :       if (ops)
     588             :       {
     589      919387 :         op =  mkopU(j,si,U,data,R);
     590      919388 :         if (op) { nbop++; gel(*ops, nbop) = op; }
     591             :       }
     592     3674793 :       gcoeff(H,i,j) = zero;
     593     3674793 :       gcoeff(H,i,si) = d;
     594     3674793 :       if (R->red && !smallop)
     595             :       {
     596     2512286 :         gen_redcol(gel(H,si), i-1, data, R);
     597     2512286 :         gen_redcol(gel(H,j), i-1, data, R);
     598             :       }
     599             :     }
     600             : 
     601     2707141 :     if (early_abort)
     602             :     {
     603     1456350 :       d = gcoeff(H,i,si);
     604     1456350 :       u = R->unit(data, d);
     605     1456643 :       if (u) d = gel(u,1);
     606     1456643 :       if (!R->equal1(d))
     607             :       {
     608        4914 :         if (ops) *ops = d;
     609        4914 :         return NULL;
     610             :       }
     611             :     }
     612             : 
     613     2702507 :     if (gc_needed(av,1))
     614             :     {
     615          20 :       if (DEBUGMEM>1) pari_warn(warnmem,"gen_howell[1]. i=%ld",i);
     616          20 :       gerepileall(av1,ops?2:1,&H,ops);
     617             :     }
     618             :   }
     619             : 
     620      930652 :   if (!ops)
     621      201691 :     lastinv = gen_last_inv_diago(H, data, R);
     622             : 
     623             :   /* put in reduced Howell form */
     624      930652 :   if (!only_triangular)
     625             :   {
     626     3770853 :     for (i=m,si=s+m; i>0; i--,si--) /* si = s+i */
     627             :     {
     628             :       /* normalize diagonal coefficient */
     629     2840186 :       if (i<=lastinv) /* lastinv>0 => !ops */
     630       33978 :         gcoeff(H,i,si) = one;
     631             :       else
     632             :       {
     633     2806208 :         u = R->unit(data,gcoeff(H,i,si));
     634     2806455 :         if (u)
     635             :         {
     636     2489406 :           g = gel(u,1);
     637     2489406 :           u = gel(u,2);
     638     2489406 :           gel(H,si) = gen_rightmulcol(gel(H,si), u, i-1, 1, data, R);
     639     2489377 :           gcoeff(H,i,si) = g;
     640     2489377 :           if (R->red) gen_redcol(gel(H,si), i-1, data, R);
     641     2489357 :           if (ops)
     642             :           {
     643     1882268 :             op = mkopmul(si,u,R);
     644     1882266 :             if (op) { nbop++; gel(*ops,nbop) = op; }
     645             :           }
     646             :         }
     647      317049 :         else if (R->red) gcoeff(H,i,si) = R->red(data, gcoeff(H,i,si));
     648             :       }
     649     2840382 :       piv = gcoeff(H,i,si);
     650             : 
     651             :       /* reduce above diagonal */
     652     2840382 :       if (!R->equal0(piv))
     653             :       {
     654     2523320 :         C = gel(H,si);
     655     6535618 :         for (j=si+1; j<=n; j++)
     656             :         {
     657     4012302 :           if (i<=lastinv) /* lastinv>0 => !ops */
     658       49581 :             gcoeff(H,i,j) = zero;
     659             :           else
     660             :           {
     661     3962721 :             gcoeff(H,i,j) = R->red(data, gcoeff(H,i,j));
     662     3962716 :             if (R->equal1(piv)) { q = gcoeff(H,i,j); r = zero; }
     663     1640663 :             else                q = R->lquo(data, gcoeff(H,i,j), piv, &r);
     664     3962704 :             q = R->neg(data,q);
     665     3962728 :             gen_addrightmul(H, C, q, j, i-1, data, R);
     666     3962733 :             if (ops)
     667             :             {
     668     2309250 :               op = mkoptransv(j,si,q,data,R);
     669     2309234 :               if (op) { nbop++; gel(*ops,nbop) = op; }
     670             :             }
     671     3962717 :             gcoeff(H,i,j) = r;
     672             :           }
     673             :         }
     674             :       }
     675             : 
     676             :       /* ensure Howell property */
     677     2840355 :       if (i>1)
     678             :       {
     679     1909958 :         a = R->rann(data, piv);
     680     1909684 :         if (!R->equal0(a))
     681             :         {
     682      568381 :           gel(H,1) = gen_rightmulcol(gel(H,si), a, i-1, 1, data, R);
     683      568381 :           if (gel(H,1) == gel(H,si)) gel(H,1) = shallowcopy(gel(H,1)); /* in case rightmulcol cheated */
     684      568381 :           if (ops)
     685             :           {
     686      148484 :             op = mkoptransv(1,si,a,data,R);
     687      148484 :             if (op) { nbop++; gel(*ops,nbop) = op; }
     688             :           }
     689     2414923 :           for (i2=i-1,si2=s+i2; i2>0; i2--,si2--)
     690             :           {
     691     1846542 :             if (R->red) gcoeff(H,i2,1) = R->red(data, gcoeff(H,i2,1));
     692     1846542 :             if (R->equal0(gcoeff(H,i2,1))) continue;
     693      285694 :             if (R->red) gcoeff(H,i2,si2) = R->red(data, gcoeff(H,i2,si2));
     694      285694 :             if (R->equal0(gcoeff(H,i2,si2)))
     695             :             {
     696      131959 :               swap(gel(H,1), gel(H,si2));
     697      131959 :               if (ops) { nbop++; gel(*ops,nbop) = mkopswap(1,si2); }
     698      131959 :               continue;
     699             :             }
     700      153735 :             U = R->extgcd(data, gcoeff(H,i2,1), gcoeff(H,i2,si2), &smallop);
     701      153735 :             d = gel(U,1);
     702      153735 :             U = gel(U,2);
     703      153735 :             gen_elem(H, U, 1, si2, i2-1, data, R);
     704      153735 :             if (ops)
     705             :             {
     706       36043 :               op = mkopU(1,si2,U,data,R);
     707       36043 :               if (op) { nbop++; gel(*ops,nbop) = op; }
     708             :             }
     709      153735 :             gcoeff(H,i2,1) = zero;
     710      153735 :             gcoeff(H,i2,si2) = d;
     711      153735 :             if (R->red && !smallop)
     712             :             {
     713      140232 :               gen_redcol(gel(H,si2), i2, data, R);
     714      140232 :               gen_redcol(gel(H,1), i2-1, data, R);
     715             :             }
     716             :           }
     717             :         }
     718             :       }
     719             : 
     720     2840129 :       if (gc_needed(av,1))
     721             :       {
     722          12 :         if (DEBUGMEM>1) pari_warn(warnmem,"gen_howell[2]. i=%ld",i);
     723          12 :         gerepileall(av1,ops?3:2,&H,&piv,ops);
     724             :       }
     725             :     }
     726             :   }
     727             : 
     728      930638 :   if (R->red)
     729     5191633 :     for (j=1; j<=n; j++)
     730             :     {
     731     4261153 :       lim = maxss(0,m-n+j);
     732     4261129 :       gen_redcol(gel(H,j), lim, data, R);
     733    17292609 :       for (i=lim+1; i<=m; i++) gcoeff(H,i,j) = zero;
     734             :     }
     735             : 
     736             :   /* put zero columns first */
     737      930406 :   iszero = cgetg(n+1,t_VECSMALL);
     738             : 
     739      930656 :   nbz = 0;
     740     5191989 :   for (i=1; i<=n; i++)
     741             :   {
     742     4261271 :     iszero[i] = gen_is_zerocol(gel(H,i), maxss(0,m-n+i), data, R);
     743     4261333 :     if (iszero[i]) nbz++;
     744             :   }
     745             : 
     746      930718 :   j = 1;
     747      930718 :   if (permute_zerocols)
     748             :   {
     749      156093 :     perm = cgetg(n+1, t_VECSMALL);
     750      913563 :     for (i=1; i<=n; i++)
     751      757470 :       if (iszero[i])
     752             :       {
     753      320439 :         perm[j] = i;
     754      320439 :         j++;
     755             :       }
     756             :   }
     757      774625 :   else perm = cgetg(n-nbz+1, t_VECSMALL);
     758     5192217 :   for (i=1; i<=n; i++)
     759     4261500 :     if (!iszero[i])
     760             :     {
     761     2523686 :       perm[j] = i;
     762     2523686 :       j++;
     763             :     }
     764             : 
     765      930717 :   if (permute_zerocols || remove_zerocols==2) H = vecpermute(H, perm);
     766      930721 :   if (permute_zerocols && remove_zerocols==2) H = vecslice(H, nbz+1, n);
     767      930721 :   if (remove_zerocols==1) H = vecslice(H, s+1, n);
     768      930721 :   if (permute_zerocols && ops) { nbop++; gel(*ops,nbop) = perm; }
     769             : 
     770      930721 :   if (ops) { setlg(*ops, nbop+1); } /* should have nbop <= maxop */
     771             : 
     772      930719 :   return H;
     773             : }
     774             : 
     775             : static GEN
     776      101780 : 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)
     777             : {
     778      101780 :   pari_sp av = avma;
     779      101780 :   GEN H = gen_howell_i(A, remove_zerocols, permute_zerocols, early_abort, only_triangular, ops, data, R);
     780      101780 :   return gc_all(av, ops?2:1, &H, ops);
     781             : }
     782             : 
     783             : static GEN
     784      157794 : gen_matimage(GEN A, GEN* U, void *data, const struct bb_hermite *R)
     785             : {
     786             :   GEN ops, H;
     787      157794 :   if (U)
     788             :   {
     789       56014 :     pari_sp av = avma, av1;
     790             :     long m, n, i, r, n2, pergc;
     791       56014 :     RgM_dimensions(A,&m,&n);
     792       56014 :     H = gen_howell_i(A, 2, 1, 0, 0, &ops, data, R);
     793       56014 :     av1 = avma;
     794       56014 :     r = lg(H)-1;
     795       56014 :     *U = shallowmatconcat(mkvec2(gen_zeromat(n, maxss(0,m-n+1), data, R), gen_matid_hermite(n, data, R)));
     796       56014 :     n2 = lg(*U)-1;
     797       56014 :     pergc = maxss(m,n);
     798      430479 :     for (i=1; i<lg(ops); i++)
     799             :     {
     800      374465 :       gen_rightapply(*U, gel(ops,i), data, R);
     801      374465 :       if (!(i%pergc) && gc_needed(av,1))
     802             :       {
     803           0 :         if (DEBUGMEM>1) pari_warn(warnmem,"gen_matimage. i=%ld",i);
     804           0 :         gerepileall(av1,1,U);
     805             :       }
     806             :     }
     807       56014 :     if (r<n2) *U = vecslice(*U, n2-r+1, n2);
     808       56014 :     return gc_all(av, 2, &H, U);
     809             :   }
     810      101780 :   return gen_howell(A, 2, 0, 0, 0, NULL, data, R);
     811             : }
     812             : 
     813             : static GEN
     814             : /* H in true Howell form: no zero columns */
     815       99911 : gen_kernel_howell(GEN H, void *data, const struct bb_hermite *R)
     816             : {
     817             :   GEN K, piv, FK;
     818             :   long m, n, j, j2, i;
     819       99911 :   RgM_dimensions(H,&m,&n);
     820       99911 :   K = gen_zeromat(n, n, data, R);
     821      415772 :   for (j=n,i=m; j>0; j--)
     822             :   {
     823      535731 :     while (R->equal0(gcoeff(H,i,j))) i--;
     824      315861 :     piv = gcoeff(H,i,j);
     825      315861 :     if (R->equal0(piv)) continue;
     826      315861 :     gcoeff(K,j,j) = R->rann(data, piv);
     827      315861 :     if (j<n)
     828             :     {
     829      215950 :       FK = gen_matmul_hermite(matslice(H,i,i,j+1,n), matslice(K, j+1, n, j+1, n), data, R);
     830      852383 :       for (j2=j+1; j2<=n; j2++)
     831      636433 :         gcoeff(K,j,j2) = R->neg(data, R->lquo(data, gcoeff(FK,1,j2-j), piv, NULL));
     832             :         /* remainder has to be zero */
     833             :     }
     834             :   }
     835       99911 :   return K;
     836             : }
     837             : 
     838             : static GEN
     839             : /* (H,ops) Howell form of A, n = number of columns of A, return a kernel of A */
     840       99981 : gen_kernel_from_howell(GEN H, GEN ops, long n, void *data, const struct bb_hermite *R)
     841             : {
     842             :   pari_sp av;
     843             :   GEN K, KH, zC;
     844             :   long m, r, n2, nbz, i, o, extra, j;
     845       99981 :   RgM_dimensions(H,&m,&r);
     846       99981 :   if (!r) return gen_matid_hermite(n, data, R); /* zerology: what if 0==1 in R? */
     847       99911 :   n2 = maxss(n,m+1);
     848       99911 :   extra = n2-n;
     849       99911 :   nbz = n2-r;
     850             :   /* compute kernel of augmented matrix */
     851       99911 :   KH = gen_kernel_howell(H, data, R);
     852       99911 :   zC = gen_zerocol(nbz, data, R);
     853       99911 :   K = cgetg(nbz+r+1, t_MAT);
     854      331170 :   for (i=1; i<=nbz; i++)
     855      231259 :     gel(K,i) = gen_colei(nbz+r, i, data, R);
     856      415772 :   for (i=1; i<=r; i++)
     857      315861 :     gel(K,nbz+i) = shallowconcat(zC, gel(KH,i));
     858      647031 :   for (i=1; i<lg(K); i++)
     859             :   {
     860      547120 :     av = avma;
     861    13654860 :     for (o=lg(ops)-1; o>0; o--)
     862    13107740 :       gen_leftapply(gel(K,i), gel(ops,o), data, R);
     863      547120 :     gen_redcol(gel(K,i), nbz+r, data, R);
     864      547120 :     gerepileall(av, 1, &gel(K,i));
     865             :   }
     866             :   /* deduce kernel of original matrix */
     867       99911 :   K = rowpermute(K, cyclic_perm(n2,extra));
     868       99911 :   K = gen_howell_i(K, 2, 0, 0, 0, NULL, data, R);
     869      226702 :   for (j=lg(K)-1, i=n2; j>0; j--)
     870             :   {
     871      301098 :     while (R->equal0(gcoeff(K,i,j))) i--;
     872      198058 :     if (i<=n) return matslice(K, 1, n, 1, j);
     873             :   }
     874       28644 :   return cgetg(1,t_MAT);
     875             : }
     876             : 
     877             : /* not GC-clean */
     878             : static GEN
     879       56126 : gen_kernel(GEN A, GEN* im, void *data, const struct bb_hermite *R)
     880             : {
     881       56126 :   pari_sp av = avma;
     882       56126 :   long n = lg(A)-1;
     883             :   GEN H, ops, K;
     884       56126 :   H = gen_howell_i(A, 2, 1, 0, 0, &ops, data, R);
     885       56126 :   gerepileall(av,2,&H,&ops);
     886       56126 :   K = gen_kernel_from_howell(H, ops, n, data, R);
     887       56126 :   if (im) *im = H;
     888       56126 :   return K;
     889             : }
     890             : 
     891             : /* right inverse, not GC-clean */
     892             : static GEN
     893      591814 : gen_inv(GEN A, void* data, const struct bb_hermite *R)
     894             : {
     895             :   pari_sp av;
     896      591814 :   GEN ops, H, U, un=R->s(data,1);
     897             :   long m,n,j,o,n2;
     898      591811 :   RgM_dimensions(A,&m,&n);
     899      591811 :   av = avma;
     900      591811 :   H = gen_howell_i(A, 0, 0, 1, 0, &ops, data, R);
     901      591807 :   if (!H) pari_err_INV("gen_inv", ops);
     902      572893 :   n2 = lg(H)-1;
     903      572893 :   ops = gerepilecopy(av,ops); /* get rid of H */
     904      572922 :   U = gen_zeromat(n2, m, data, R);
     905     2018369 :   for (j=1; j<=m; j++)
     906     1445459 :     gcoeff(U,j+n2-m,j) = un;
     907     2018369 :   for (j=1; j<=m; j++)
     908             :   {
     909     1445444 :     av = avma;
     910     3092250 :     for (o=lg(ops)-1; o>0; o--)
     911     1646881 :       gen_leftapply(gel(U,j), gel(ops,o), data, R);
     912     1445369 :     gen_redcol(gel(U,j), n2, data, R);
     913     1444969 :     gerepileall(av, 1, &gel(U,j));
     914             :   }
     915      572925 :   if (n2>n) U = rowslice(U, n2-n+1, n2);
     916      572909 :   return U;
     917             : }
     918             : 
     919             : /* H true Howell form (no zero columns).
     920             :  * Compute Z = Y - HX canonical representative of R^m mod H(R^n) */
     921             : static GEN
     922       43953 : gen_reduce_mod_howell(GEN H, GEN Y, GEN *X, void* data, const struct bb_hermite *R)
     923             : {
     924       43953 :   pari_sp av = avma;
     925             :   long m,n,i,j;
     926             :   GEN Z, q, r, C;
     927       43953 :   RgM_dimensions(H,&m,&n);
     928       43953 :   if (X) *X = gen_zerocol(n,data,R);
     929       43953 :   Z = shallowcopy(Y);
     930       43953 :   i = m;
     931      155099 :   for (j=n; j>0; j--)
     932             :   {
     933      180348 :     while (R->equal0(gcoeff(H,i,j))) i--;
     934      111146 :     q = R->lquo(data, gel(Z,i), gcoeff(H,i,j), &r);
     935      111146 :     gel(Z,i) = r;
     936      111146 :     C = gen_rightmulcol(gel(H,j), R->neg(data,q), i, 0, data, R);
     937      111146 :     if (C) gen_addcol(Z, C, i-1, data, R);
     938      111146 :     if (X) gel(*X,j) = q;
     939             :   }
     940       43953 :   gen_redcol(Z, lg(Z)-1, data, R);
     941       43953 :   return gc_all(av, X?2:1, &Z, X);
     942             : }
     943             : 
     944             : /* H: Howell form of A */
     945             : /* (m,n): dimensions of original matrix A */
     946             : /* return NULL if no solution */
     947             : static GEN
     948       43953 : gen_solve_from_howell(GEN H, GEN ops, long m, long n, GEN Y, void* data, const struct bb_hermite *R)
     949             : {
     950       43953 :   pari_sp av = avma;
     951             :   GEN Z, X;
     952             :   long n2, mH, nH, i;
     953       43953 :   RgM_dimensions(H,&mH,&nH);
     954       43953 :   n2 = maxss(n,m+1);
     955       43953 :   Z = gen_reduce_mod_howell(H, Y, &X, data, R);
     956       43953 :   if (!gen_is_zerocol(Z,m,data,R)) { set_avma(av); return NULL; }
     957             : 
     958       43904 :   X = shallowconcat(zerocol(n2-nH),X);
     959      442113 :   for (i=lg(ops)-1; i>0; i--) gen_leftapply(X, gel(ops,i), data, R);
     960       43904 :   X = vecslice(X, n2-n+1, n2);
     961       43904 :   gen_redcol(X, n, data, R);
     962       43904 :   return gerepilecopy(av, X);
     963             : }
     964             : 
     965             : /* return NULL if no solution, not GC-clean */
     966             : static GEN
     967       43953 : gen_solve(GEN A, GEN Y, GEN* K, void* data, const struct bb_hermite *R)
     968             : {
     969             :   GEN H, ops, X;
     970             :   long m,n;
     971             : 
     972       43953 :   RgM_dimensions(A,&m,&n);
     973       43953 :   if (!n) m = lg(Y)-1;
     974       43953 :   H = gen_howell_i(A, 2, 1, 0, 0, &ops, data, R);
     975       43953 :   X = gen_solve_from_howell(H, ops, m, n, Y, data, R);
     976       43953 :   if (!X) return NULL;
     977       43904 :   if (K) *K = gen_kernel_from_howell(H, ops, n, data, R);
     978       43904 :   return X;
     979             : }
     980             : 
     981             : GEN
     982      157829 : matimagemod(GEN A, GEN d, GEN* U)
     983             : {
     984             :   void *data;
     985             :   const struct bb_hermite* R;
     986      157829 :   if (typ(A)!=t_MAT || !RgM_is_ZM(A)) pari_err_TYPE("matimagemod", A);
     987      157815 :   if (typ(d)!=t_INT) pari_err_TYPE("matimagemod", d);
     988      157808 :   if (signe(d)<=0) pari_err_DOMAIN("matimagemod", "d", "<=", gen_0, d);
     989      157801 :   if (equali1(d)) return cgetg(1,t_MAT);
     990      157794 :   R = get_Fp_hermite(&data, d);
     991      157794 :   return gen_matimage(A, U, data, R);
     992             : }
     993             : 
     994             : #if 0
     995             : /* for testing purpose */
     996             : GEN
     997             : ZM_hnfmodid2(GEN A, GEN d)
     998             : {
     999             :   pari_sp av = avma;
    1000             :   void *data;
    1001             :   long i;
    1002             :   const struct bb_hermite* R = get_Fp_hermite(&data, d);
    1003             :   GEN H;
    1004             :   if (typ(A)!=t_MAT || !RgM_is_ZM(A)) pari_err_TYPE("ZM_hnfmodid2", A);
    1005             :   if (typ(d)!=t_INT) pari_err_TYPE("ZM_hnfmodid2", d);
    1006             :   H = gen_howell_i(A, 1, 0, 0, 0, NULL, data, R);
    1007             :   for (i=1; i<lg(H); i++)
    1008             :     if (!signe(gcoeff(H,i,i))) gcoeff(H,i,i) = d;
    1009             :   return gerepilecopy(av,H);
    1010             : }
    1011             : #endif
    1012             : 
    1013             : GEN
    1014          84 : matdetmod(GEN A, GEN d)
    1015             : {
    1016          84 :   pari_sp av = avma;
    1017             :   void *data;
    1018             :   const struct bb_hermite* R;
    1019          84 :   long n = lg(A)-1, i;
    1020             :   GEN D, H, ops;
    1021          84 :   if (typ(A)!=t_MAT || !RgM_is_ZM(A)) pari_err_TYPE("matdetmod", A);
    1022          70 :   if (typ(d)!=t_INT) pari_err_TYPE("matdetmod", d);
    1023          70 :   if (signe(d)<=0) pari_err_DOMAIN("matdetmod", "d", "<=", gen_0, d);
    1024          63 :   if (!n) return equali1(d) ? gen_0 : gen_1;
    1025          56 :   if (n != nbrows(A)) pari_err_DIM("matdetmod");
    1026          49 :   if (equali1(d)) return gen_0;
    1027          42 :   R = get_Fp_hermite(&data, d);
    1028          42 :   H = gen_howell_i(A, 1, 0, 0, 1, &ops, data, R);
    1029          42 :   D = gen_detops(ops, data, R);
    1030          42 :   D = Fp_inv(D, d);
    1031         203 :   for (i = 1; i <= n; i++) D = Fp_mul(D, gcoeff(H,i,i), d);
    1032          42 :   return gerepileuptoint(av, D);
    1033             : }
    1034             : 
    1035             : GEN
    1036       56161 : matkermod(GEN A, GEN d, GEN* im)
    1037             : {
    1038       56161 :   pari_sp av = avma;
    1039             :   void *data;
    1040             :   const struct bb_hermite* R;
    1041             :   long m,n;
    1042             :   GEN K;
    1043       56161 :   if (typ(A)!=t_MAT || !RgM_is_ZM(A)) pari_err_TYPE("matkermod", A);
    1044       56147 :   if (typ(d)!=t_INT) pari_err_TYPE("matkermod", d);
    1045       56140 :   if (signe(d)<=0) pari_err_DOMAIN("makermod", "d", "<=", gen_0, d);
    1046       56133 :   if (equali1(d)) return cgetg(1,t_MAT);
    1047       56126 :   R = get_Fp_hermite(&data, d);
    1048       56126 :   RgM_dimensions(A,&m,&n);
    1049       56126 :   if (!im && m>2*n) /* TODO tune treshold */
    1050        4592 :     A = shallowtrans(matimagemod(shallowtrans(A),d,NULL));
    1051       56126 :   K = gen_kernel(A, im, data, R); return gc_all(av,im?2:1,&K,im);
    1052             : }
    1053             : 
    1054             : /* left inverse */
    1055             : GEN
    1056      678052 : matinvmod(GEN A, GEN d)
    1057             : {
    1058      678052 :   pari_sp av = avma;
    1059             :   void *data;
    1060             :   const struct bb_hermite* R;
    1061             :   GEN U;
    1062      678052 :   if (typ(A)!=t_MAT || !RgM_is_ZM(A)) pari_err_TYPE("matinvmod", A);
    1063      678061 :   if (typ(d)!=t_INT) pari_err_TYPE("matinvmod", d);
    1064      678054 :   if (signe(d)<=0) pari_err_DOMAIN("matinvmod", "d", "<=", gen_0, d);
    1065      678047 :   if (equali1(d)) {
    1066             :     long m,n;
    1067       86316 :     RgM_dimensions(A,&m,&n);
    1068       86317 :     if (m<n) pari_err_INV("matinvmod",A);
    1069       86310 :     return zeromatcopy(n,m);
    1070             :   }
    1071      591745 :   R = get_Fp_hermite(&data, d);
    1072      591772 :   U = gen_inv(shallowtrans(A), data, R);
    1073      572909 :   return gerepilecopy(av, shallowtrans(U));
    1074             : }
    1075             : 
    1076             : /* assume all D[i]>0, not GC-clean */
    1077             : static GEN
    1078       43953 : matsolvemod_finite(GEN M, GEN D, GEN Y, long flag)
    1079             : {
    1080             :   void *data;
    1081             :   const struct bb_hermite* R;
    1082             :   GEN X, K, d;
    1083             :   long m, n;
    1084             : 
    1085       43953 :   RgM_dimensions(M,&m,&n);
    1086       43953 :   if (typ(D)==t_COL)
    1087             :   { /* create extra variables for the D[i] */
    1088             :     GEN MD;
    1089          84 :     long i, i2, extra = 0;
    1090          84 :     d = gen_1;
    1091         231 :     for (i=1; i<lg(D); i++) d = lcmii(d,gel(D,i));
    1092         231 :     for (i=1; i<lg(D); i++) if (!equalii(gel(D,i),d)) extra++;
    1093          84 :     MD = cgetg(extra+1,t_MAT);
    1094          84 :     i2 = 1;
    1095         231 :     for (i=1; i<lg(D); i++)
    1096         147 :       if (!equalii(gel(D,i),d))
    1097             :       {
    1098          77 :         gel(MD,i2) = Rg_col_ei(gel(D,i),m,i);
    1099          77 :         i2++;
    1100             :       }
    1101          84 :     M = shallowconcat(M,MD);
    1102             :   }
    1103       43869 :   else d = D;
    1104             : 
    1105       43953 :   R = get_Fp_hermite(&data, d);
    1106       43953 :   X = gen_solve(M, Y, flag? &K: NULL, data, R);
    1107       43953 :   if (!X) return gen_0;
    1108       43904 :   X = vecslice(X,1,n);
    1109             : 
    1110       43904 :   if (flag)
    1111             :   {
    1112       43855 :     K = rowslice(K,1,n);
    1113       43855 :     K = hnfmodid(shallowconcat(zerocol(n),K),d);
    1114       43855 :     X = mkvec2(X,K);
    1115             :   }
    1116       43904 :   return X;
    1117             : }
    1118             : 
    1119             : /* Return a solution of congruence system sum M[i,j] X_j = Y_i mod D_i
    1120             :  * If pU != NULL, put in *pU a Z-basis of the homogeneous system.
    1121             :  * Works for all non negative D_i but inefficient compared to
    1122             :  * matsolvemod_finite; to be used only when one D_i is 0 */
    1123             : static GEN
    1124          70 : gaussmoduloall(GEN M, GEN D, GEN Y, GEN *pU)
    1125             : {
    1126          70 :   pari_sp av = avma;
    1127          70 :   long n, m, j, l, lM = lg(M);
    1128             :   GEN delta, H, U, u1, u2, x;
    1129             : 
    1130          70 :   if (lM == 1)
    1131             :   {
    1132          28 :     long lY = 0;
    1133          28 :     switch(typ(Y))
    1134             :     {
    1135           0 :       case t_INT: break;
    1136          28 :       case t_COL: lY = lg(Y); break;
    1137           0 :       default: pari_err_TYPE("gaussmodulo",Y);
    1138             :     }
    1139          28 :     switch(typ(D))
    1140             :     {
    1141          14 :       case t_INT: break;
    1142          14 :       case t_COL: if (lY && lY != lg(D)) pari_err_DIM("gaussmodulo");
    1143          14 :                   break;
    1144           0 :       default: pari_err_TYPE("gaussmodulo",D);
    1145             :     }
    1146          28 :     if (pU) *pU = cgetg(1, t_MAT);
    1147          28 :     return cgetg(1,t_COL);
    1148             :   }
    1149          42 :   n = nbrows(M);
    1150          42 :   switch(typ(D))
    1151             :   {
    1152          14 :     case t_COL:
    1153          14 :       if (lg(D)-1!=n) pari_err_DIM("gaussmodulo");
    1154          14 :       delta = diagonal_shallow(D); break;
    1155          28 :     case t_INT: delta = scalarmat_shallow(D,n); break;
    1156           0 :     default: pari_err_TYPE("gaussmodulo",D);
    1157             :       return NULL; /* LCOV_EXCL_LINE */
    1158             :   }
    1159          42 :   switch(typ(Y))
    1160             :   {
    1161           0 :     case t_INT: Y = const_col(n, Y); break;
    1162          42 :     case t_COL:
    1163          42 :       if (lg(Y)-1!=n) pari_err_DIM("gaussmodulo");
    1164          42 :       break;
    1165           0 :     default: pari_err_TYPE("gaussmodulo",Y);
    1166             :       return NULL; /* LCOV_EXCL_LINE */
    1167             :   }
    1168          42 :   H = ZM_hnfall_i(shallowconcat(M,delta), &U, 1);
    1169          42 :   Y = hnf_solve(H,Y); if (!Y) return gen_0;
    1170          35 :   l = lg(H); /* may be smaller than lM if some moduli are 0 */
    1171          35 :   n = l-1;
    1172          35 :   m = lg(U)-1 - n;
    1173          35 :   u1 = cgetg(m+1,t_MAT);
    1174          35 :   u2 = cgetg(n+1,t_MAT);
    1175          84 :   for (j=1; j<=m; j++) { GEN c = gel(U,j); setlg(c,lM); gel(u1,j) = c; }
    1176          35 :   U += m;
    1177          77 :   for (j=1; j<=n; j++) { GEN c = gel(U,j); setlg(c,lM); gel(u2,j) = c; }
    1178             :   /*       (u1 u2)
    1179             :    * (M D) (*  * ) = (0 H) */
    1180          35 :   u1 = ZM_lll(u1, 0.75, LLL_INPLACE);
    1181          35 :   Y = ZM_ZC_mul(u2,Y);
    1182          35 :   x = ZC_reducemodmatrix(Y, u1);
    1183          35 :   if (!pU) x = gerepileupto(av, x);
    1184             :   else
    1185             :   {
    1186          14 :     gerepileall(av, 2, &x, &u1);
    1187          14 :     *pU = u1;
    1188             :   }
    1189          35 :   return x;
    1190             : }
    1191             : /* to be used when one D_i is 0 */
    1192             : static GEN
    1193          70 : msolvemod0(GEN M, GEN D, GEN Y, long flag)
    1194             : {
    1195          70 :   pari_sp av = avma;
    1196             :   GEN y, x, U;
    1197          70 :   if (!flag) return gaussmoduloall(M,D,Y,NULL);
    1198          28 :   y = cgetg(3,t_VEC);
    1199          28 :   x = gaussmoduloall(M,D,Y,&U);
    1200          28 :   if (x == gen_0) { set_avma(av); return gen_0; }
    1201          21 :   gel(y,1) = x;
    1202          21 :   gel(y,2) = U; return y;
    1203             : 
    1204             : }
    1205             : GEN
    1206       44149 : matsolvemod(GEN M, GEN D, GEN Y, long flag)
    1207             : {
    1208       44149 :   pari_sp av = avma;
    1209       44149 :   long m, n, i, char0 = 0;
    1210       44149 :   if (typ(M)!=t_MAT || !RgM_is_ZM(M)) pari_err_TYPE("matsolvemod (M)",M);
    1211       44135 :   RgM_dimensions(M,&m,&n);
    1212       44135 :   if (typ(D)!=t_INT && (typ(D)!=t_COL || !RgV_is_ZV(D)))
    1213          28 :     pari_err_TYPE("matsolvemod (D)",D);
    1214       44107 :   if (n)
    1215       43960 :     { if (typ(D)==t_COL && lg(D)!=m+1) pari_err_DIM("matsolvemod [1]"); }
    1216             :   else
    1217         147 :     { if (typ(D)==t_COL) m = lg(D)-1; }
    1218       44100 :   if (typ(Y)==t_INT)
    1219       43890 :     Y = const_col(m,Y);
    1220         210 :   else if (typ(Y)!=t_COL || !RgV_is_ZV(Y)) pari_err_TYPE("matsolvemod (Y)",Y);
    1221       44072 :   if (typ(D)==t_COL && typ(Y)==t_COL && lg(D)!=lg(Y))
    1222          28 :     pari_err_DIM("matsolvemod [2]");
    1223       44044 :   if (!n && !m) m = lg(Y)-1;
    1224       43988 :   else if (m != lg(Y)-1) pari_err_DIM("matsolvemod [3]");
    1225       44037 :   if (typ(D)==t_INT)
    1226             :   {
    1227       43918 :     if (signe(D)<0) pari_err_DOMAIN("matsolvemod","D","<",gen_0,D);
    1228       43911 :     if (!signe(D)) char0 = 1;
    1229             :   }
    1230             :   else /*typ(D)==t_COL*/
    1231         301 :     for (i=1; i<=m; i++)
    1232             :     {
    1233         189 :       if (signe(gel(D,i))<0)
    1234           7 :         pari_err_DOMAIN("matsolvemod","D[i]","<",gen_0,gel(D,i));
    1235         182 :       if (!signe(gel(D,i))) char0 = 1;
    1236             :     }
    1237       87976 :   return gerepilecopy(av, char0? msolvemod0(M,D,Y,flag)
    1238       43953 :                                : matsolvemod_finite(M,D,Y,flag));
    1239             : }
    1240             : GEN
    1241           0 : gaussmodulo2(GEN M, GEN D, GEN Y) { return matsolvemod(M,D,Y, 1); }
    1242             : GEN
    1243           0 : gaussmodulo(GEN M, GEN D, GEN Y) { return matsolvemod(M,D,Y, 0); }

Generated by: LCOV version 1.16