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 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.10.0 lcov report (development 21183-045fc66) Lines: 435 458 95.0 %
Date: 2017-10-19 06:23:01 Functions: 42 44 95.5 %
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             : struct bb_hermite
      26             : {
      27             :   GEN (*add)(void*, GEN, GEN);
      28             :   GEN (*neg)(void*, GEN);
      29             :   GEN (*mul)(void*, GEN, GEN);
      30             :   GEN (*extgcd)(void*, GEN, GEN, int*);
      31             :   GEN (*rann)(void*, GEN);
      32             :   GEN (*lquo)(void*, GEN, GEN, GEN*);
      33             :   GEN (*unit)(void*, GEN);
      34             :   int (*equal0)(GEN);
      35             :   int (*equal1)(GEN);
      36             :   GEN (*s)(void*, long);
      37             :   GEN (*red)(void*, GEN);
      38             : };
      39             : 
      40             : static GEN
      41   121169308 : _Fp_add(void *data, GEN x, GEN y) { (void) data; return addii(x,y); }
      42             : 
      43             : static GEN
      44    10314304 : _Fp_neg(void *data, GEN x) { (void) data; return negi(x); }
      45             : 
      46             : static GEN
      47   145301880 : _Fp_mul(void *data, GEN x, GEN y) { (void) data; return mulii(x,y); }
      48             : 
      49             : static GEN
      50     1358413 : _Fp_rann(void *data, GEN x)
      51             : {
      52     1358413 :   GEN d, N = (GEN)data;
      53     1358413 :   if (!signe(x)) return gen_1;
      54     1088402 :   d = gcdii(x,N);
      55     1088402 :   return modii(diviiexact(N,d),N);
      56             : }
      57             : 
      58             : static GEN
      59     4259178 : _Fp_lquo(void *data, GEN x, GEN y, GEN* r) { (void) data; return truedvmdii(x,y,r); }
      60             : 
      61             : /* D=MN, p|M => !p|a, p|N => p|a, return M */
      62             : static GEN
      63           0 : Z_split(GEN D, GEN a)
      64             : {
      65             :   long i, n;
      66             :   GEN N;
      67           0 :   n = expi(D);
      68           0 :   n = n<2 ? 1 : expu(n)+1;
      69           0 :   for (i=1;i<=n;i++)
      70           0 :     a = Fp_sqr(a,D);
      71           0 :   N = gcdii(a,D);
      72           0 :   return diviiexact(D,N);
      73             : }
      74             : 
      75             : /* c s.t. gcd(a+cb,N) = gcd(a,b,N) without factoring */
      76             : static GEN
      77           0 : Z_stab(GEN a, GEN b, GEN N)
      78             : {
      79             :   GEN g, a2, N2;
      80           0 :   g = gcdii(a,b);
      81           0 :   g = gcdii(g,N);
      82           0 :   N2 = diviiexact(N,g);
      83           0 :   a2 = diviiexact(a,g);
      84           0 :   return Z_split(N2,a2);
      85             : }
      86             : 
      87             : static GEN
      88     3161298 : _Fp_unit(void *data, GEN x)
      89             : {
      90     3161298 :   GEN g,s,v,d,N=(GEN)data,N2;
      91             :   long i;
      92     3161298 :   if (!signe(x)) return NULL;
      93     2687275 :   g = bezout(x,N,&s,&v);
      94     2687275 :   if (equali1(g) || equali1(gcdii(s,N))) return mkvec2(g,s);
      95       40628 :   N2 = diviiexact(N,g);
      96       56378 :   for (i=0; i<5; i++)
      97             :   {
      98       56378 :     s = addii(s,N2);
      99       56378 :     if (equali1(gcdii(s,N))) return mkvec2(g,s);
     100             :   }
     101           0 :   d = Z_stab(s,N2,N);
     102           0 :   d = mulii(d,N2);
     103           0 :   v = Fp_add(s,d,N);
     104           0 :   if (equali1(v)) return NULL;
     105           0 :   return mkvec2(g,v);
     106             : }
     107             : 
     108             : static GEN
     109     2573846 : _Fp_extgcd(void *data, GEN x, GEN y, int* smallop)
     110             : {
     111             :   GEN d,u,v,m;
     112     2573846 :   if (equali1(y))
     113             :   {
     114     1213028 :     *smallop = 1;
     115     1213028 :     return mkvec2(y,mkmat2(
     116             :           mkcol2(gen_1,Fp_neg(x,(GEN)data)),
     117             :           mkcol2(gen_0,gen_1)));
     118             :   }
     119     1360818 :   *smallop = 0;
     120     1360818 :   d = bezout(x,y,&u,&v);
     121     1360818 :   if (!signe(d)) return mkvec2(d,matid(2));
     122     1360818 :   m = cgetg(3,t_MAT);
     123     1360818 :   m = mkmat2(
     124             :     mkcol2(diviiexact(y,d),negi(diviiexact(x,d))),
     125             :     mkcol2(u,v));
     126     1360818 :   return mkvec2(d,m);
     127             : }
     128             : 
     129             : static int
     130   210257538 : _Fp_equal0(GEN x) { return !signe(x); }
     131             : 
     132             : static int
     133    25772858 : _Fp_equal1(GEN x) { return equali1(x); }
     134             : 
     135             : static GEN
     136    11646121 : _Fp_s(void *data, long x)
     137             : {
     138    11646121 :   if (!x) return gen_0;
     139      777602 :   if (x==1) return gen_1;
     140           0 :   return modsi(x,(GEN)data);
     141             : }
     142             : 
     143             : static GEN
     144    78002125 : _Fp_red(void *data, GEN x) { return Fp_red(x, (GEN)data); }
     145             : 
     146             : /* p not necessarily prime */
     147             : static const struct bb_hermite Fp_hermite=
     148             :   {_Fp_add,_Fp_neg,_Fp_mul,_Fp_extgcd,_Fp_rann,_Fp_lquo,_Fp_unit,_Fp_equal0,_Fp_equal1,_Fp_s,_Fp_red};
     149             : 
     150             : static const struct bb_hermite*
     151      236915 : get_Fp_hermite(void **data, GEN p)
     152             : {
     153      236915 :   *data = (void*)p; return &Fp_hermite;
     154             : }
     155             : 
     156             : static void
     157    11491880 : gen_redcol(GEN C, long lim, void* data, const struct bb_hermite *R)
     158             : {
     159             :   long i;
     160   114774184 :   for (i=1; i<=lim; i++)
     161   103282304 :     if (!R->equal0(gel(C,i)))
     162    45624227 :       gel(C,i) = R->red(data, gel(C,i));
     163    11491880 : }
     164             : 
     165             : static GEN
     166             : /* return NULL if a==0 */
     167             : /* assume C*a is zero after lim */
     168    14518461 : gen_rightmulcol(GEN C, GEN a, long lim, int fillzeros, void* data, const struct bb_hermite *R)
     169             : {
     170             :   GEN Ca,zero;
     171             :   long i;
     172    14518461 :   if (R->equal1(a)) return C;
     173     9651497 :   if (R->equal0(a)) return NULL;
     174     7404427 :   Ca = cgetg(lg(C),t_COL);
     175   121662910 :   for (i=1; i<=lim; i++)
     176   114258483 :     gel(Ca,i) = R->mul(data, gel(C,i), a);
     177     7404427 :   if (fillzeros && lim+1 < lg(C))
     178             :   {
     179     4817096 :     zero = R->s(data,0);
     180    63633969 :     for (i=lim+1; i<lg(C); i++)
     181    58816873 :       gel(Ca,i) = zero;
     182             :   }
     183     7404427 :   return Ca;
     184             : }
     185             : 
     186             : static void
     187             : /* C1 <- C1 + C2 */
     188             : /* assume C2[i]==0 for i>lim */
     189     5289902 : gen_addcol(GEN C1, GEN C2, long lim, void* data, const struct bb_hermite *R)
     190             : {
     191             :   long i;
     192    99092101 :   for (i=1; i<=lim; i++)
     193    93802199 :     gel(C1,i) = R->add(data, gel(C1,i), gel(C2,i));
     194     5289902 : }
     195             : 
     196             : static void
     197             : /* H[,i] <- H[,i] + C*a */
     198             : /* assume C is zero after lim */
     199     9722643 : gen_addrightmul(GEN H, GEN C, GEN a, long i, long lim, void* data, const struct bb_hermite *R)
     200             : {
     201             :   GEN Ca;
     202    19445286 :   if (R->equal0(a)) return;
     203     2211620 :   Ca = gen_rightmulcol(C, a, lim, 0, data, R);
     204     2211620 :   gen_addcol(gel(H,i), Ca, lim, data, R);
     205             : }
     206             : 
     207             : static GEN
     208      891191 : gen_zerocol(long n, void* data, const struct bb_hermite *R)
     209             : {
     210      891191 :   GEN C = cgetg(n+1,t_COL), zero = R->s(data, 0);
     211             :   long i;
     212      891191 :   for (i=1; i<=n; i++) gel(C,i) = zero;
     213      891191 :   return C;
     214             : }
     215             : 
     216             : static GEN
     217      326928 : gen_zeromat(long m, long n, void* data, const struct bb_hermite *R)
     218             : {
     219      326928 :   GEN M = cgetg(n+1,t_MAT);
     220             :   long i;
     221      326928 :   for (i=1; i<=n; i++) gel(M,i) = gen_zerocol(m, data, R);
     222      326928 :   return M;
     223             : }
     224             : 
     225             : static GEN
     226      451157 : gen_colei(long n, long i, void* data, const struct bb_hermite *R)
     227             : {
     228      451157 :   GEN C = cgetg(n+1,t_COL), zero = R->s(data, 0);
     229             :   long j;
     230    13386800 :   for (j=1; j<=n; j++)
     231    12935643 :     if (i!=j)   gel(C,j) = zero;
     232      451157 :     else        gel(C,j) = R->s(data,1);
     233      451157 :   return C;
     234             : }
     235             : 
     236             : static GEN
     237       56021 : gen_matid_hermite(long n, void* data, const struct bb_hermite *R)
     238             : {
     239       56021 :   GEN M = cgetg(n+1,t_MAT);
     240             :   long i;
     241       56021 :   for (i=1; i<=n; i++) gel(M,i) = gen_colei(n, i, data, R);
     242       56021 :   return M;
     243             : }
     244             : 
     245             : static GEN
     246     4417637 : gen_matmul_hermite(GEN A, GEN B, void* data, const struct bb_hermite *R)
     247             : {
     248     4417637 :   GEN M,sum,prod,zero = R->s(data,0);
     249             :   long a,b,c,c2,i,j,k;
     250     4417637 :   RgM_dimensions(A,&a,&c);
     251     4417637 :   RgM_dimensions(B,&c2,&b);
     252     4417637 :   if (c!=c2) pari_err_DIM("gen_matmul_hermite");
     253     4417637 :   M = cgetg(b+1,t_MAT);
     254     9432682 :   for (j=1; j<=b; j++)
     255             :   {
     256     5015045 :     gel(M,j) = cgetg(a+1,t_COL);
     257    14277788 :     for (i=1; i<=a; i++)
     258             :     {
     259     9262743 :       sum = zero;
     260    32739574 :       for (k=1; k<=c; k++)
     261             :       {
     262    23476831 :         prod = R->mul(data, gcoeff(A,i,k), gcoeff(B,k,j));
     263    23476831 :         sum = R->add(data, sum, prod);
     264             :       }
     265     9262743 :       gcoeff(M,i,j) = sum;
     266             :     }
     267     5015045 :     gen_redcol(gel(M,j), a, data, R);
     268             :   }
     269     4417637 :   return M;
     270             : }
     271             : 
     272             : static void
     273             : /* U = [u1,u2]~, C <- A*u1 + B*u2 */
     274             : /* assume both A, B and C are zero after lim */
     275     5325352 : gen_rightlincomb(GEN A, GEN B, GEN U, GEN *C, long lim, void* data, const struct bb_hermite *R)
     276             : {
     277             :   GEN Au1, Bu2;
     278     5325352 :   Au1 = gen_rightmulcol(A, gel(U,1), lim, 1, data, R);
     279     5325352 :   Bu2 = gen_rightmulcol(B, gel(U,2), lim, 1, data, R);
     280     5325352 :   if (!Au1 && !Bu2) { *C = gen_zerocol(lg(A)-1, data, R); return; }
     281     5325352 :   if (!Au1) { *C = Bu2; return; }
     282     3426758 :   if (!Bu2) { *C = Au1; return; }
     283     3078282 :   gen_addcol(Au1, Bu2, lim, data, R);
     284     3078282 :   *C = Au1;
     285             : }
     286             : 
     287             : static void
     288             : /* (H[,i] | H[,j]) <- (H[,i] | H[,j]) * U */
     289             : /* assume both columns are zero after lim */
     290     2662676 : gen_elem(GEN H, GEN U, long i, long j, long lim, void* data, const struct bb_hermite *R)
     291             : {
     292             :   GEN Hi, Hj;
     293     2662676 :   Hi = shallowcopy(gel(H,i));
     294     2662676 :   Hj = shallowcopy(gel(H,j));
     295     2662676 :   gen_rightlincomb(Hi, Hj, gel(U,1), &gel(H,i), lim, data, R);
     296     2662676 :   gen_rightlincomb(Hi, Hj, gel(U,2), &gel(H,j), lim, data, R);
     297     2662676 : }
     298             : 
     299             : static int
     300             : /* assume C is zero after lim */
     301     1790866 : gen_is_zerocol(GEN C, long lim, void* data, const struct bb_hermite *R)
     302             : {
     303             :   long i;
     304             :   (void) data;
     305    14796285 :   for (i=1; i<=lim; i++)
     306    14037800 :     if (!R->equal0(gel(C,i))) return 0;
     307      758485 :   return 1;
     308             : }
     309             : 
     310             : static GEN
     311             : /* Ci <- Ci + Cj*a */
     312     2083606 : mkoptransv(long i, long j, GEN a, void* data, const struct bb_hermite *R)
     313             : {
     314     2083606 :   a = R->red(data,a);
     315     2083606 :   if (R->equal0(a)) return NULL;
     316     1358483 :   return mkvec2(mkvecsmall2(i,j),a);
     317             : }
     318             : 
     319             : static GEN
     320             : /* (Ci|Cj) <- (Ci|Cj)*U */
     321      718375 : mkopU(long i, long j, GEN U, void* data, const struct bb_hermite *R)
     322             : {
     323      718375 :   if (R->equal1(gcoeff(U,1,1)) && R->equal0(gcoeff(U,1,2))
     324      400533 :       && R->equal1(gcoeff(U,2,2))) return mkoptransv(i,j,gcoeff(U,2,1),data,R);
     325      318066 :   return mkvec2(mkvecsmall3(i,j,0),U);
     326             : }
     327             : 
     328             : static GEN
     329             : /* Ci <- Ci*u */
     330      380632 : mkopmul(long i, GEN u, const struct bb_hermite *R)
     331             : {
     332      380632 :   if (R->equal1(u)) return NULL;
     333       80003 :   return mkvec2(mkvecsmall(i),u);
     334             : }
     335             : 
     336             : static GEN
     337             : /* Ci <-> Cj */
     338       59346 : mkopswap(long i, long j)
     339             : {
     340       59346 :   return mkvec(mkvecsmall2(i,j));
     341             : }
     342             : 
     343             : static void
     344      374444 : gen_rightapply(GEN M, GEN op, void* data, const struct bb_hermite *R)
     345             : {
     346             :   GEN M2, ind, X;
     347      374444 :   long i, j, m = lg(gel(M,1))-1;
     348      374444 :   switch (typ(op))
     349             :   {
     350             :     case t_VECSMALL:
     351       56007 :       M2 = vecpermute(M,op);
     352       56007 :       for (i=1; i<lg(M); i++) gel(M,i) = gel(M2,i);
     353       56007 :       return;
     354             :     case t_VEC:
     355      318437 :       ind = gel(op,1);
     356      318437 :       switch (lg(op))
     357             :       {
     358             :         case 2:
     359       16135 :           swap(gel(M,ind[1]),gel(M,ind[2]));
     360       16135 :           return;
     361             :         case 3:
     362      302302 :           X = gel(op,2);
     363      302302 :           i = ind[1];
     364      302302 :           switch (lg(ind))
     365             :           {
     366             :             case 2:
     367       37786 :               gel(M,i) = gen_rightmulcol(gel(M,i), X, m, 0, data, R);
     368       37786 :               gen_redcol(gel(M,i), m, data, R);
     369       37786 :               return;
     370             :             case 3:
     371      175686 :               gen_addrightmul(M, gel(M,ind[2]), X, i, m, data, R);
     372      175686 :               gen_redcol(gel(M,i), m, data, R);
     373      175686 :               return;
     374             :             case 4:
     375       88830 :               j = ind[2];
     376       88830 :               gen_elem(M, X, i, j, m, data, R);
     377       88830 :               gen_redcol(gel(M,i), m, data, R);
     378       88830 :               gen_redcol(gel(M,j), m, data, R);
     379       88830 :               return;
     380             :           }
     381             :       }
     382             :   }
     383             : }
     384             : 
     385             : static void
     386    55726349 : gen_leftapply(GEN C, GEN op, void* data, const struct bb_hermite *R)
     387             : {
     388             :   GEN C2, ind, X;
     389             :   long i, j;
     390    55726349 :   switch (typ(op))
     391             :   {
     392             :     case t_VECSMALL:
     393      535591 :       C2 = vecpermute(C,perm_inv(op));
     394      535591 :       for (i=1; i<lg(C); i++) gel(C,i) = gel(C2,i);
     395      535591 :       return;
     396             :     case t_VEC:
     397    55190758 :       ind = gel(op,1);
     398    55190758 :       switch (lg(op))
     399             :       {
     400             :         case 2:
     401     1921794 :           swap(gel(C,ind[1]),gel(C,ind[2]));
     402     1921794 :           return;
     403             :         case 3:
     404    53268964 :           X = gel(op,2);
     405    53268964 :           i = ind[1];
     406    53268964 :           switch (lg(ind))
     407             :           {
     408             :             case 2:
     409      522802 :               gel(C,i) = R->mul(data, X, gel(C,i));
     410      522802 :               gel(C,i) = R->red(data, gel(C,i));
     411      522802 :               return;
     412             :             case 3:
     413    48498464 :               j = ind[2];
     414    48498464 :               if (R->equal0(gel(C,i))) return;
     415     3890278 :               gel(C,j) = R->add(data, gel(C,j), R->mul(data, X, gel(C,i)));
     416     3890278 :               return;
     417             :             case 4:
     418     4247698 :               j = ind[2];
     419     4247698 :               C2 = gen_matmul_hermite(X, mkmat(mkcol2(gel(C,i),gel(C,j))), data, R);
     420     4247698 :               gel(C,i) = gcoeff(C2,1,1);
     421     4247698 :               gel(C,j) = gcoeff(C2,2,1);
     422     4247698 :               return;
     423             :           }
     424             :       }
     425             :   }
     426             : }
     427             : 
     428             : static int
     429      161483 : gen_is_inv(GEN x, void* data, const struct bb_hermite *R)
     430             : {
     431      161483 :   GEN u = R->unit(data, x);
     432      161483 :   if (!u) return R->equal1(x);
     433       61701 :   return R->equal1(gel(u,1));
     434             : }
     435             : 
     436             : static long
     437      145887 : gen_last_inv_diago(GEN A, void* data, const struct bb_hermite *R)
     438             : {
     439             :   long i,m,n,j;
     440      145887 :   RgM_dimensions(A,&m,&n);
     441      175707 :   for (i=1,j=n-m+1; i<=m; i++,j++)
     442      161483 :     if (!gen_is_inv(gcoeff(A,i,j),data,R)) return i-1;
     443       14224 :   return m;
     444             : }
     445             : 
     446             : static GEN
     447             : /* remove_zerocols: 0 none, 1 until square, 2 all */
     448             : /* early abort: if not right-invertible, abort, return NULL, and set ops to the
     449             :  * non-invertible pivot */
     450      291438 : gen_howell_i(GEN A, long remove_zerocols, long permute_zerocols, long early_abort, GEN* ops, void *data, const struct bb_hermite *R)
     451             : {
     452      291438 :   pari_sp av = avma;
     453      291438 :   GEN H,U,piv,u,q,a,perm,iszero,C,zero=R->s(data,0),d,g,r,op,one=R->s(data,1);
     454      291438 :   long m,n,i,j,s,si,i2,si2,nbz,lim,extra,maxop=0,nbop=0,lastinv=0;
     455             :   int smallop;
     456             : 
     457      291438 :   RgM_dimensions(A,&m,&n);
     458      291438 :   if (early_abort && n<m)
     459             :   {
     460       14000 :     if (ops) *ops = zero;
     461       14000 :     return NULL;
     462             :   }
     463      277438 :   if (n<m+1)
     464             :   {
     465      200305 :     extra = m+1-n;
     466      200305 :     H = shallowmatconcat(mkvec2(gen_zeromat(m,extra,data,R),A));
     467             :   }
     468             :   else
     469             :   {
     470       77133 :     extra = 0;
     471       77133 :     H = RgM_shallowcopy(A);
     472             :   }
     473      277438 :   RgM_dimensions(H,&m,&n);
     474      277438 :   s = n-m; /* shift */
     475             : 
     476      277438 :   if(ops)
     477             :   {
     478      131551 :     maxop = m*n + (m*(m+1))/2 + 1;
     479      131551 :     *ops = zerovec(maxop); /* filled with placeholders so gerepile can handle it */
     480             :   }
     481             : 
     482             :   /* put in triangular form */
     483     1383522 :   for (i=m,si=s+m; i>0 && si>extra; i--,si--) /* si = s+i */
     484             :   {
     485     1110998 :     if (R->red) gcoeff(H,i,si) = R->red(data, gcoeff(H,i,si));
     486             :     /* bottom-right diagonal */
     487    11124141 :     for (j = extra+1; j < si; j++)
     488             :     {
     489    10013143 :       if (R->red) gcoeff(H,i,j) = R->red(data, gcoeff(H,i,j));
     490    10013143 :       if (R->equal0(gcoeff(H,i,j))) continue;
     491     2364029 :       U = R->extgcd(data, gcoeff(H,i,j), gcoeff(H,i,si), &smallop);
     492     2364029 :       d = gel(U,1);
     493     2364029 :       U = gel(U,2);
     494     2364029 :       if (n>10)
     495             :       {
     496             :         /* normalize diagonal coefficient -> faster reductions on this row */
     497     1576743 :         u = R->unit(data, d);
     498     1576743 :         if (u)
     499             :         {
     500     1576743 :           g = gel(u,1);
     501     1576743 :           u = gel(u,2);
     502     1576743 :           gcoeff(U,1,2) = R->mul(data, gcoeff(U,1,2), u);
     503     1576743 :           gcoeff(U,2,2) = R->mul(data, gcoeff(U,2,2), u);
     504     1576743 :           d = g;
     505             :         }
     506             :       }
     507     2364029 :       gen_elem(H, U, j, si, i-1, data, R);
     508     2364029 :       if (ops)
     509             :       {
     510      609721 :         op =  mkopU(j,si,U,data,R);
     511      609721 :         if (op) { nbop++; gel(*ops, nbop) = op; }
     512             :       }
     513     2364029 :       gcoeff(H,i,j) = zero;
     514     2364029 :       gcoeff(H,i,si) = d;
     515     2364029 :       if (R->red && !smallop)
     516             :       {
     517     1159318 :         gen_redcol(gel(H,si), i-1, data, R);
     518     1159318 :         gen_redcol(gel(H,j), i-1, data, R);
     519             :       }
     520             :     }
     521             : 
     522             : 
     523     1110998 :     if (early_abort)
     524             :     {
     525       46417 :       d = gcoeff(H,i,si);
     526       46417 :       u = R->unit(data, d);
     527       46417 :       if (u) d = gel(u,1);
     528       46417 :       if (!R->equal1(d))
     529             :       {
     530        4914 :         if (ops) *ops = d;
     531        4914 :         return NULL;
     532             :       }
     533             :     }
     534             : 
     535     1106084 :     if (gc_needed(av,1))
     536             :     {
     537           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"gen_howell[1]. i=%ld",i);
     538           0 :       gerepileall(av,ops?2:1,&H,ops);
     539             :     }
     540             :   }
     541             : 
     542      272524 :   if (!ops)
     543      145887 :     lastinv = gen_last_inv_diago(H, data, R);
     544             : 
     545             :   /* put in reduced Howell form */
     546     1678999 :   for (i=m,si=s+m; i>0; i--,si--) /* si = s+i */
     547             :   {
     548             :     /* normalize diagonal coefficient */
     549     1406475 :     if (i<=lastinv) /* lastinv>0 => !ops */
     550       29820 :       gcoeff(H,i,si) = one;
     551             :     else
     552             :     {
     553     1376655 :       u = R->unit(data,gcoeff(H,i,si));
     554     1376655 :       if (u)
     555             :       {
     556     1002561 :         g = gel(u,1);
     557     1002561 :         u = gel(u,2);
     558     1002561 :         gel(H,si) = gen_rightmulcol(gel(H,si), u, i-1, 1, data, R);
     559     1002561 :         gcoeff(H,i,si) = g;
     560     1002561 :         if (R->red) gen_redcol(gel(H,si), i-1, data, R);
     561     1002561 :         if (ops)
     562             :         {
     563      380632 :           op = mkopmul(si,u,R);
     564      380632 :           if (op) { nbop++; gel(*ops,nbop) = op; }
     565             :         }
     566             :       }
     567      374094 :       else if (R->red) gcoeff(H,i,si) = R->red(data, gcoeff(H,i,si));
     568             :     }
     569     1406475 :     piv = gcoeff(H,i,si);
     570             : 
     571             :     /* reduce above diagonal */
     572     1406475 :     if (!R->equal0(piv))
     573             :     {
     574     1032381 :       C = gel(H,si);
     575    10604804 :       for (j=si+1; j<=n; j++)
     576             :       {
     577     9572423 :         if (i<=lastinv) /* lastinv>0 => !ops */
     578       25466 :           gcoeff(H,i,j) = zero;
     579             :         else
     580             :         {
     581     9546957 :           gcoeff(H,i,j) = R->red(data, gcoeff(H,i,j));
     582     9546957 :           if (R->equal1(piv)) { q = gcoeff(H,i,j); r = zero; }
     583     3491831 :           else                q = R->lquo(data, gcoeff(H,i,j), piv, &r);
     584     9546957 :           q = R->neg(data,q);
     585     9546957 :           gen_addrightmul(H, C, q, j, i-1, data, R);
     586     9546957 :           if (ops)
     587             :           {
     588     1374345 :             op = mkoptransv(j,si,q,data,R);
     589     1374345 :             if (op) { nbop++; gel(*ops,nbop) = op; }
     590             :           }
     591     9546957 :           gcoeff(H,i,j) = r;
     592             :         }
     593             :       }
     594             :     }
     595             : 
     596             :     /* ensure Howell property */
     597     1406475 :     if (i>1)
     598             :     {
     599     1133951 :       a = R->rann(data, piv);
     600     1133951 :       if (!R->equal0(a))
     601             :       {
     602      615790 :         gel(H,1) = gen_rightmulcol(gel(H,si), a, i-1, 1, data, R);
     603      615790 :         if (gel(H,1) == gel(H,si)) gel(H,1) = shallowcopy(gel(H,1)); /* in case rightmulcol cheated */
     604      615790 :         if (ops)
     605             :         {
     606      308952 :           op = mkoptransv(1,si,a,data,R);
     607      308952 :           if (op) { nbop++; gel(*ops,nbop) = op; }
     608             :         }
     609     9001230 :         for (i2=i-1,si2=s+i2; i2>0; i2--,si2--)
     610             :         {
     611     8385440 :           if (R->red) gcoeff(H,i2,1) = R->red(data, gcoeff(H,i2,1));
     612     8385440 :           if (R->equal0(gcoeff(H,i2,1))) continue;
     613      340858 :           if (R->red) gcoeff(H,i2,si2) = R->red(data, gcoeff(H,i2,si2));
     614      340858 :           if (R->equal0(gcoeff(H,i2,si2)))
     615             :           {
     616      131041 :             swap(gel(H,1), gel(H,si2));
     617      131041 :             if (ops) { nbop++; gel(*ops,nbop) = mkopswap(1,si2); }
     618      131041 :             continue;
     619             :           }
     620      209817 :           U = R->extgcd(data, gcoeff(H,i2,1), gcoeff(H,i2,si2), &smallop);
     621      209817 :           d = gel(U,1);
     622      209817 :           U = gel(U,2);
     623      209817 :           gen_elem(H, U, 1, si2, i2-1, data, R);
     624      209817 :           if (ops)
     625             :           {
     626      108654 :             op = mkopU(1,si2,U,data,R);
     627      108654 :             if (op) { nbop++; gel(*ops,nbop) = op; }
     628             :           }
     629      209817 :           gcoeff(H,i2,1) = zero;
     630      209817 :           gcoeff(H,i2,si2) = d;
     631      209817 :           if (R->red && !smallop)
     632             :           {
     633      201500 :             gen_redcol(gel(H,si2), i2, data, R);
     634      201500 :             gen_redcol(gel(H,1), i2-1, data, R);
     635             :           }
     636             :         }
     637             :       }
     638             :     }
     639             : 
     640     1406475 :     if (gc_needed(av,1))
     641             :     {
     642           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"gen_howell[2]. i=%ld",i);
     643           0 :       gerepileall(av,ops?3:2,&H,&piv,ops);
     644             :     }
     645             :   }
     646             : 
     647      272524 :   if (R->red)
     648     2063390 :     for (j=1; j<=n; j++)
     649             :     {
     650     1790866 :       lim = maxss(0,m-n+j);
     651     1790866 :       gen_redcol(gel(H,j), lim, data, R);
     652     1790866 :       for (i=lim+1; i<=m; i++) gcoeff(H,i,j) = zero;
     653             :     }
     654             : 
     655             :   /* put zero columns first */
     656      272524 :   iszero = cgetg(n+1,t_VECSMALL);
     657             : 
     658      272524 :   nbz = 0;
     659     2063390 :   for (i=1; i<=n; i++)
     660             :   {
     661     1790866 :     iszero[i] = gen_is_zerocol(gel(H,i), maxss(0,m-n+i), data, R);
     662     1790866 :     if (iszero[i]) nbz++;
     663             :   }
     664             : 
     665      272524 :   j = 1;
     666      272524 :   if (permute_zerocols)
     667             :   {
     668      110544 :     perm = cgetg(n+1, t_VECSMALL);
     669      856191 :     for (i=1; i<=n; i++)
     670      745647 :       if (iszero[i])
     671             :       {
     672      400064 :         perm[j] = i;
     673      400064 :         j++;
     674             :       }
     675             :   }
     676      161980 :   else perm = cgetg(n-nbz+1, t_VECSMALL);
     677     2063390 :   for (i=1; i<=n; i++)
     678     1790866 :     if (!iszero[i])
     679             :     {
     680     1032381 :       perm[j] = i;
     681     1032381 :       j++;
     682             :     }
     683             : 
     684      272524 :   if (permute_zerocols || remove_zerocols==2) H = vecpermute(H, perm);
     685      272524 :   if (permute_zerocols && remove_zerocols==2) H = vecslice(H, nbz+1, n);
     686      272524 :   if (remove_zerocols==1) H = vecslice(H, s+1, n);
     687      272524 :   if (permute_zerocols && ops) { nbop++; gel(*ops,nbop) = perm; }
     688             : 
     689      272524 :   if (ops) { setlg(*ops, nbop+1); } /* should have nbop <= maxop */
     690             : 
     691      272524 :   return H;
     692             : }
     693             : 
     694             : static GEN
     695       91364 : gen_howell(GEN A, long remove_zerocols, long permute_zerocols, long early_abort, GEN* ops, void *data, const struct bb_hermite *R)
     696             : {
     697       91364 :   pari_sp av = avma;
     698       91364 :   GEN H = gen_howell_i(A, remove_zerocols, permute_zerocols, early_abort, ops, data, R);
     699       91364 :   gerepileall(av, ops?2:1, &H, ops);
     700       91364 :   return H;
     701             : }
     702             : 
     703             : static GEN
     704      147371 : gen_matimage(GEN A, GEN* U, void *data, const struct bb_hermite *R)
     705             : {
     706             :   GEN ops, H;
     707      147371 :   if (U)
     708             :   {
     709       56007 :     pari_sp av = avma;
     710             :     long m, n, i, r, n2;
     711       56007 :     RgM_dimensions(A,&m,&n);
     712       56007 :     H = gen_howell_i(A, 2, 1, 0, &ops, data, R);
     713       56007 :     r = lg(H)-1;
     714       56007 :     *U = shallowmatconcat(mkvec2(gen_zeromat(n, maxss(0,m-n+1), data, R), gen_matid_hermite(n, data, R)));
     715       56007 :     n2 = lg(*U)-1;
     716      430451 :     for (i=1; i<lg(ops); i++)
     717      374444 :       gen_rightapply(*U, gel(ops,i), data, R);
     718       56007 :     if (r<n2) *U = vecslice(*U, n2-r+1, n2);
     719       56007 :     gerepileall(av, 2, &H, U);
     720       56007 :     return H;
     721             :   }
     722       91364 :   else return gen_howell(A, 2, 0, 0, NULL, data, R);
     723             : }
     724             : 
     725             : static GEN
     726             : /* H in true Howell form: no zero columns */
     727       54523 : gen_kernel_howell(GEN H, void *data, const struct bb_hermite *R)
     728             : {
     729             :   GEN K, piv, FK;
     730             :   long m, n, j, j2, i;
     731       54523 :   RgM_dimensions(H,&m,&n);
     732       54523 :   K = gen_zeromat(n, n, data, R);
     733      278985 :   for (j=n,i=m; j>0; j--)
     734             :   {
     735      224462 :     while (R->equal0(gcoeff(H,i,j))) i--;
     736      224462 :     piv = gcoeff(H,i,j);
     737      224462 :     if (R->equal0(piv)) continue;
     738      224462 :     gcoeff(K,j,j) = R->rann(data, piv);
     739      224462 :     if (j<n)
     740             :     {
     741      169939 :       FK = gen_matmul_hermite(matslice(H,i,i,j+1,n), matslice(K, j+1, n, j+1, n), data, R);
     742      937286 :       for (j2=j+1; j2<=n; j2++)
     743      767347 :         gcoeff(K,j,j2) = R->neg(data, R->lquo(data, gcoeff(FK,1,j2-j), piv, NULL));
     744             :         /* remainder has to be zero */
     745             :     }
     746             :   }
     747       54523 :   return K;
     748             : }
     749             : 
     750             : static GEN
     751             : /* (H,ops) Howell form of A, n = number of columns of A, return a kernel of A */
     752       54537 : gen_kernel_from_howell(GEN H, GEN ops, long n, void *data, const struct bb_hermite *R)
     753             : {
     754             :   pari_sp av;
     755             :   GEN K, KH, zC;
     756             :   long m, r, n2, nbz, i, o, extra, j;
     757       54537 :   RgM_dimensions(H,&m,&r);
     758       54537 :   if (!r) return gen_matid_hermite(n, data, R); /* zerology: what if 0==1 in R? */
     759       54523 :   n2 = maxss(n,m+1);
     760       54523 :   extra = n2-n;
     761       54523 :   nbz = n2-r;
     762             :   /* compute kernel of augmented matrix */
     763       54523 :   KH = gen_kernel_howell(H, data, R);
     764       54523 :   zC = gen_zerocol(nbz, data, R);
     765       54523 :   K = cgetg(nbz+r+1, t_MAT);
     766      365652 :   for (i=1; i<=nbz; i++)
     767      311129 :     gel(K,i) = gen_colei(nbz+r, i, data, R);
     768      278985 :   for (i=1; i<=r; i++)
     769      224462 :     gel(K,nbz+i) = shallowconcat(zC, gel(KH,i));
     770      590114 :   for (i=1; i<lg(K); i++)
     771             :   {
     772      535591 :     av = avma;
     773    55920746 :     for (o=lg(ops)-1; o>0; o--)
     774    55385155 :       gen_leftapply(gel(K,i), gel(ops,o), data, R);
     775      535591 :     gen_redcol(gel(K,i), nbz+r, data, R);
     776      535591 :     gerepileall(av, 1, &gel(K,i));
     777             :   }
     778             :   /* deduce kernel of original matrix */
     779       54523 :   K = rowpermute(K, cyclic_perm(n2,extra));
     780       54523 :   K = gen_howell_i(K, 2, 0, 0, NULL, data, R);
     781      317037 :   for (j=lg(K)-1, i=n2; j>0; j--)
     782             :   {
     783      302792 :     while (R->equal0(gcoeff(K,i,j))) i--;
     784      302792 :     if (i<=n) return matslice(K, 1, n, 1, j);
     785             :   }
     786       14245 :   return cgetg(1,t_MAT);
     787             : }
     788             : 
     789             : static GEN
     790       54537 : gen_kernel(GEN A, GEN* im, void *data, const struct bb_hermite *R)
     791             : {
     792       54537 :   pari_sp av = avma;
     793       54537 :   long n = lg(A)-1;
     794             :   GEN H, ops, K;
     795       54537 :   H = gen_howell_i(A, 2, 1, 0, &ops, data, R);
     796       54537 :   gerepileall(av,2,&H,&ops);
     797       54537 :   K = gen_kernel_from_howell(H, ops, n, data, R);
     798       54537 :   if (im) *im = H;
     799       54537 :   gerepileall(av,im?2:1,&K,im);
     800       54537 :   return K;
     801             : }
     802             : 
     803             : /* right inverse */
     804             : static GEN
     805       35007 : gen_inv(GEN A, void* data, const struct bb_hermite *R)
     806             : {
     807       35007 :   pari_sp av0 = avma, av;
     808       35007 :   GEN ops, H, U, un=R->s(data,1);
     809             :   long m,n,j,o,n2;
     810       35007 :   RgM_dimensions(A,&m,&n);
     811       35007 :   av = avma;
     812       35007 :   H = gen_howell_i(A, 0, 0, 1, &ops, data, R);
     813       35007 :   if (!H) pari_err_INV("gen_inv", ops);
     814       16093 :   n2 = lg(H)-1;
     815       16093 :   ops = gerepilecopy(av,ops); /* get rid of H */
     816       16093 :   U = gen_zeromat(n2, m, data, R);
     817       51142 :   for (j=1; j<=m; j++)
     818       35049 :     gcoeff(U,j+n2-m,j) = un;
     819       51142 :   for (j=1; j<=m; j++)
     820             :   {
     821       35049 :     av = avma;
     822      376243 :     for (o=lg(ops)-1; o>0; o--)
     823      341194 :       gen_leftapply(gel(U,j), gel(ops,o), data, R);
     824       35049 :     gen_redcol(gel(U,j), n2, data, R);
     825       35049 :     gerepileall(av, 1, &gel(U,j));
     826             :   }
     827       16093 :   if (n2>n) U = rowslice(U, n2-n+1, n2);
     828       16093 :   return gerepilecopy(av0,U);
     829             : }
     830             : 
     831             : GEN
     832      147392 : matimagemod(GEN A, GEN d, GEN* U)
     833             : {
     834             :   void* data;
     835             :   long i,j;
     836             :   const struct bb_hermite* R;
     837      147392 :   if (typ(A)!=t_MAT) pari_err_TYPE("matimagemod[1]", A);
     838      147385 :   if (typ(d)!=t_INT) pari_err_TYPE("matimagemod[2]", d);
     839      565005 :   for (j=1; j<lg(A); j++)
     840     1718206 :     for (i=1; i<lg(gel(A,j)); i++)
     841     1300579 :       if (typ(gcoeff(A,i,j))!=t_INT) pari_err_TYPE("matimagemod[3]", gcoeff(A,i,j));
     842      147371 :   R = get_Fp_hermite(&data, d);
     843      147371 :   return gen_matimage(A, U, data, R);
     844             : }
     845             : 
     846             : /* for testing purpose */
     847             : /*
     848             : GEN
     849             : ZM_hnfmodid2(GEN A, GEN d)
     850             : {
     851             :   pari_sp av = avma;
     852             :   void* data;
     853             :   long i, j;
     854             :   GEN H;
     855             :   if (typ(A)!=t_MAT) pari_err_TYPE("ZM_hnfmodid2", A);
     856             :   if (typ(d)!=t_INT) pari_err_TYPE("ZM_hnfmodid2", d);
     857             :   for (j=1; j<lg(A); j++)
     858             :     for (i=1; i<lg(gel(A,j)); i++)
     859             :       if (typ(gcoeff(A,i,j))!=t_INT) pari_err_TYPE("ZM_hnfmodid2", gcoeff(A,i,j));
     860             :   H = gen_howell_i(A, 1, 0, 0, NULL, data, get_Fp_hermite(&data, d));
     861             :   for (i=1; i<lg(H); i++)
     862             :     if (!signe(gcoeff(H,i,i))) gcoeff(H,i,i) = d;
     863             :   return gerepilecopy(av,H);
     864             : }
     865             : */
     866             : 
     867             : GEN
     868       54558 : matkermod(GEN A, GEN d, GEN* im)
     869             : {
     870             :   void* data;
     871             :   long i,j;
     872             :   const struct bb_hermite* R;
     873       54558 :   if (typ(A)!=t_MAT) pari_err_TYPE("matkermod[1]", A);
     874       54551 :   if (typ(d)!=t_INT) pari_err_TYPE("matkermod[2]", d);
     875      327642 :   for (j=1; j<lg(A); j++)
     876     3380398 :     for (i=1; i<lg(gel(A,j)); i++)
     877     3107300 :       if (typ(gcoeff(A,i,j))!=t_INT) pari_err_TYPE("matkermod[3]", gcoeff(A,i,j));
     878       54537 :   R = get_Fp_hermite(&data, d);
     879       54537 :   return gen_kernel(A, im, data, R);
     880             : }
     881             : 
     882             : /* left inverse */
     883             : GEN
     884       35028 : matinvmod(GEN A, GEN d)
     885             : {
     886       35028 :   pari_sp av = avma;
     887             :   void* data;
     888             :   long i,j;
     889             :   const struct bb_hermite* R;
     890             :   GEN B, U;
     891       35028 :   if (typ(A)!=t_MAT) pari_err_TYPE("matinvmod[1]", A);
     892       35021 :   if (typ(d)!=t_INT) pari_err_TYPE("matinvmod[2]", d);
     893      140042 :   for (j=1; j<lg(A); j++)
     894      420112 :     for (i=1; i<lg(gel(A,j)); i++)
     895      315084 :       if (typ(gcoeff(A,i,j))!=t_INT) pari_err_TYPE("matinvmod[3]", gcoeff(A,i,j));
     896       35007 :   B = shallowtrans(A);
     897       35007 :   R = get_Fp_hermite(&data, d);
     898       35007 :   U = gen_inv(B, data, R);
     899       16093 :   U = shallowtrans(U);
     900       16093 :   return gerepilecopy(av,U);
     901             : }
     902             : 

Generated by: LCOV version 1.11