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 - lll.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.14.0 lcov report (development 27775-aca467eab2) Lines: 1203 1254 95.9 %
Date: 2022-07-03 07:33:15 Functions: 92 93 98.9 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2008  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : #include "pari.h"
      16             : #include "paripriv.h"
      17             : 
      18             : #define DEBUGLEVEL DEBUGLEVEL_qflll
      19             : 
      20             : static double
      21    18781733 : pari_rint(double a)
      22             : {
      23             : #ifdef HAS_RINT
      24    18781733 :   return rint(a);
      25             : #else
      26             :   const double two_to_52 = 4.5035996273704960e+15;
      27             :   double fa = fabs(a);
      28             :   double r = two_to_52 + fa;
      29             :   if (fa >= two_to_52) {
      30             :     r = a;
      31             :   } else {
      32             :     r = r - two_to_52;
      33             :     if (a < 0) r = -r;
      34             :   }
      35             :   return r;
      36             : #endif
      37             : }
      38             : 
      39             : /* default quality ratio for LLL */
      40             : static const double LLLDFT = 0.99;
      41             : 
      42             : /* assume flag & (LLL_KER|LLL_IM|LLL_ALL). LLL_INPLACE implies LLL_IM */
      43             : static GEN
      44      192883 : lll_trivial(GEN x, long flag)
      45             : {
      46      192883 :   if (lg(x) == 1)
      47             :   { /* dim x = 0 */
      48       15401 :     if (! (flag & LLL_ALL)) return cgetg(1,t_MAT);
      49          28 :     retmkvec2(cgetg(1,t_MAT), cgetg(1,t_MAT));
      50             :   }
      51             :   /* dim x = 1 */
      52      177482 :   if (gequal0(gel(x,1)))
      53             :   {
      54        1386 :     if (flag & LLL_KER) return matid(1);
      55        1386 :     if (flag & (LLL_IM|LLL_INPLACE)) return cgetg(1,t_MAT);
      56          28 :     retmkvec2(matid(1), cgetg(1,t_MAT));
      57             :   }
      58      176094 :   if (flag & LLL_INPLACE) return gcopy(x);
      59       73005 :   if (flag & LLL_KER) return cgetg(1,t_MAT);
      60       73005 :   if (flag & LLL_IM)  return matid(1);
      61          27 :   retmkvec2(cgetg(1,t_MAT), (flag & LLL_GRAM)? gcopy(x): matid(1));
      62             : }
      63             : 
      64             : /* vecslice(x,#x-k,#x) in place. Works for t_MAT, t_VEC/t_COL */
      65             : static GEN
      66     2335986 : vectail_inplace(GEN x, long k)
      67             : {
      68     2335986 :   if (!k) return x;
      69       57258 :   x[k] = ((ulong)x[0] & ~LGBITS) | evallg(lg(x) - k);
      70       57262 :   return x + k;
      71             : }
      72             : 
      73             : /* k = dim Kernel */
      74             : static GEN
      75     2409079 : lll_finish(GEN h, long k, long flag)
      76             : {
      77             :   GEN g;
      78     2409079 :   if (!(flag & (LLL_IM|LLL_KER|LLL_ALL|LLL_INPLACE))) return h;
      79     2336009 :   if (flag & (LLL_IM|LLL_INPLACE)) return vectail_inplace(h, k);
      80          90 :   if (flag & LLL_KER) { setlg(h,k+1); return h; }
      81          62 :   g = vecslice(h,1,k); /* done first: vectail_inplace kills h */
      82          70 :   return mkvec2(g, vectail_inplace(h, k));
      83             : }
      84             : 
      85             : /* y * z * 2^e, e >= 0; y,z t_INT */
      86             : INLINE GEN
      87     1946899 : mulshift(GEN y, GEN z, long e)
      88             : {
      89     1946899 :   long ly = lgefint(y), lz;
      90             :   pari_sp av;
      91             :   GEN t;
      92     1946899 :   if (ly == 2) return gen_0;
      93      206968 :   lz = lgefint(z);
      94      206968 :   av = avma; (void)new_chunk(ly+lz+nbits2lg(e)); /* HACK */
      95      206968 :   t = mulii(z, y);
      96      206968 :   set_avma(av); return shifti(t, e);
      97             : }
      98             : 
      99             : /* x - y * z * 2^e, e >= 0; x,y,z t_INT */
     100             : INLINE GEN
     101     6914675 : submulshift(GEN x, GEN y, GEN z, long e)
     102             : {
     103     6914675 :   long lx = lgefint(x), ly, lz;
     104             :   pari_sp av;
     105             :   GEN t;
     106     6914675 :   if (!e) return submulii(x, y, z);
     107     6879000 :   if (lx == 2) { t = mulshift(y, z, e); togglesign(t); return t; }
     108     5107627 :   ly = lgefint(y);
     109     5107627 :   if (ly == 2) return icopy(x);
     110     4629875 :   lz = lgefint(z);
     111     4629875 :   av = avma; (void)new_chunk(lx+ly+lz+nbits2lg(e)); /* HACK */
     112     4629875 :   t = shifti(mulii(z, y), e);
     113     4629875 :   set_avma(av); return subii(x, t);
     114             : }
     115             : 
     116             : /* x - u*y * 2^e */
     117             : INLINE GEN
     118    23535325 : submuliu2n(GEN x, GEN y, ulong u, long e)
     119             : {
     120             :   pari_sp av;
     121    23535325 :   long ly = lgefint(y);
     122    23535325 :   if (ly == 2) return x;
     123    16550041 :   av = avma;
     124    16550041 :   (void)new_chunk(3+ly+lgefint(x)+nbits2lg(e)); /* HACK */
     125    16551733 :   y = shifti(mului(u,y), e);
     126    16550806 :   set_avma(av); return subii(x, y);
     127             : }
     128             : /* x + u*y * 2^e */
     129             : INLINE GEN
     130    23802883 : addmuliu2n(GEN x, GEN y, ulong u, long e)
     131             : {
     132             :   pari_sp av;
     133    23802883 :   long ly = lgefint(y);
     134    23802883 :   if (ly == 2) return x;
     135    16666249 :   av = avma;
     136    16666249 :   (void)new_chunk(3+ly+lgefint(x)+nbits2lg(e)); /* HACK */
     137    16667836 :   y = shifti(mului(u,y), e);
     138    16667119 :   set_avma(av); return addii(x, y);
     139             : }
     140             : 
     141             : /* n < 10; gerepileall supporting &NULL arguments. Maybe rename and export ? */
     142             : INLINE void
     143        7526 : gc_lll(pari_sp av, int n, ...)
     144             : {
     145             :   int i, j;
     146             :   GEN *gptr[10];
     147             :   size_t s;
     148        7526 :   va_list a; va_start(a, n);
     149       24200 :   for (i=j=0; i<n; i++)
     150             :   {
     151       16674 :     GEN *x = va_arg(a,GEN*);
     152       16674 :     if (*x) { gptr[j++] = x; *x = (GEN)copy_bin(*x); }
     153             :   }
     154        7526 :   va_end(a); set_avma(av);
     155       18894 :   for (--j; j>=0; j--) *gptr[j] = bin_copy((GENbin*)*gptr[j]);
     156        7526 :   s = pari_mainstack->top - pari_mainstack->bot;
     157             :   /* size of saved objects ~ stacksize / 4 => overflow */
     158        7526 :   if (av - avma > (s >> 2))
     159             :   {
     160           0 :     size_t t = avma - pari_mainstack->bot;
     161           0 :     av = avma; new_chunk((s + t) / sizeof(long)); set_avma(av); /* double */
     162             :   }
     163        7526 : }
     164             : 
     165             : /********************************************************************/
     166             : /**                                                                **/
     167             : /**                   FPLLL (adapted from D. Stehle's code)        **/
     168             : /**                                                                **/
     169             : /********************************************************************/
     170             : /* Babai* and fplll* are a conversion to libpari API and data types
     171             :    of fplll-1.3 by Damien Stehle'.
     172             : 
     173             :   Copyright 2005, 2006 Damien Stehle'.
     174             : 
     175             :   This program is free software; you can redistribute it and/or modify it
     176             :   under the terms of the GNU General Public License as published by the
     177             :   Free Software Foundation; either version 2 of the License, or (at your
     178             :   option) any later version.
     179             : 
     180             :   This program implements ideas from the paper "Floating-point LLL Revisited",
     181             :   by Phong Nguyen and Damien Stehle', in the Proceedings of Eurocrypt'2005,
     182             :   Springer-Verlag; and was partly inspired by Shoup's NTL library:
     183             :   http://www.shoup.net/ntl/ */
     184             : 
     185             : /* x t_REAL, |x| >= 1/2. Test whether |x| <= 3/2 */
     186             : static int
     187     1095308 : absrsmall2(GEN x)
     188             : {
     189     1095308 :   long e = expo(x), l, i;
     190     1095308 :   if (e < 0) return 1;
     191      344999 :   if (e > 0 || (ulong)x[2] > (3UL << (BITS_IN_LONG-2))) return 0;
     192             :   /* line above assumes l > 2. OK since x != 0 */
     193      178320 :   l = lg(x); for (i = 3; i < l; i++) if (x[i]) return 0;
     194      151059 :   return 1;
     195             : }
     196             : /* x t_REAL; test whether |x| <= 1/2 */
     197             : static int
     198     2403848 : absrsmall(GEN x)
     199             : {
     200             :   long e, l, i;
     201     2403848 :   if (!signe(x)) return 1;
     202     2393033 :   e = expo(x); if (e < -1) return 1;
     203     1101635 :   if (e > -1 || (ulong)x[2] > HIGHBIT) return 0;
     204        7571 :   l = lg(x); for (i = 3; i < l; i++) if (x[i]) return 0;
     205        6327 :   return 1;
     206             : }
     207             : 
     208             : static void
     209    43794359 : rotate(GEN A, long k2, long k)
     210             : {
     211             :   long i;
     212    43794359 :   GEN B = gel(A,k2);
     213   140759355 :   for (i = k2; i > k; i--) gel(A,i) = gel(A,i-1);
     214    43794359 :   gel(A,k) = B;
     215    43794359 : }
     216             : 
     217             : /************************* FAST version (double) ************************/
     218             : #define dmael(x,i,j) ((x)[i][j])
     219             : #define del(x,i) ((x)[i])
     220             : 
     221             : static double *
     222    37037420 : cget_dblvec(long d)
     223    37037420 : { return (double*) stack_malloc_align(d*sizeof(double), sizeof(double)); }
     224             : 
     225             : static double **
     226     9307924 : cget_dblmat(long d) { return (double **) cgetg(d, t_VECSMALL); }
     227             : 
     228             : static double
     229   274753563 : itodbl_exp(GEN x, long *e)
     230             : {
     231   274753563 :   pari_sp av = avma;
     232   274753563 :   GEN r = itor(x,DEFAULTPREC);
     233   274685636 :   *e = expo(r); setexpo(r,0);
     234   274687389 :   return gc_double(av, rtodbl(r));
     235             : }
     236             : 
     237             : static double
     238   186155931 : dbldotproduct(double *x, double *y, long n)
     239             : {
     240             :   long i;
     241   186155931 :   double sum = del(x,1) * del(y,1);
     242  2969147867 :   for (i=2; i<=n; i++) sum += del(x,i) * del(y,i);
     243   186155931 :   return sum;
     244             : }
     245             : 
     246             : static double
     247     2674820 : dbldotsquare(double *x, long n)
     248             : {
     249             :   long i;
     250     2674820 :   double sum = del(x,1) * del(x,1);
     251     8606492 :   for (i=2; i<=n; i++) sum += del(x,i) * del(x,i);
     252     2674820 :   return sum;
     253             : }
     254             : 
     255             : static long
     256    32611705 : set_line(double *appv, GEN v, long n)
     257             : {
     258    32611705 :   long i, maxexp = 0;
     259    32611705 :   pari_sp av = avma;
     260    32611705 :   GEN e = cgetg(n+1, t_VECSMALL);
     261   307338846 :   for (i = 1; i <= n; i++)
     262             :   {
     263   274752809 :     del(appv,i) = itodbl_exp(gel(v,i), e+i);
     264   274725985 :     if (e[i] > maxexp) maxexp = e[i];
     265             :   }
     266   307468948 :   for (i = 1; i <= n; i++) del(appv,i) = ldexp(del(appv,i), e[i]-maxexp);
     267    32586037 :   set_avma(av); return maxexp;
     268             : }
     269             : 
     270             : static void
     271    47891975 : dblrotate(double **A, long k2, long k)
     272             : {
     273             :   long i;
     274    47891975 :   double *B = del(A,k2);
     275   153958996 :   for (i = k2; i > k; i--) del(A,i) = del(A,i-1);
     276    47891975 :   del(A,k) = B;
     277    47891975 : }
     278             : /* update G[kappa][i] from appB */
     279             : static void
     280    30127160 : setG_fast(double **appB, long n, double **G, long kappa, long a, long b)
     281             : { long i;
     282   152203063 :   for (i = a; i <= b; i++)
     283   122076433 :     dmael(G,kappa,i) = dbldotproduct(del(appB,kappa), del(appB,i), n);
     284    30126630 : }
     285             : /* update G[i][kappa] from appB */
     286             : static void
     287    24199319 : setG2_fast(double **appB, long n, double **G, long kappa, long a, long b)
     288             : { long i;
     289    88284639 :   for (i = a; i <= b; i++)
     290    64085570 :     dmael(G,i,kappa) = dbldotproduct(del(appB,kappa), del(appB,i), n);
     291    24199069 : }
     292             : const long EX0 = -2; /* uninitialized; any value less than expo(0.51) = -1 */
     293             : 
     294             : #ifdef LONG_IS_64BIT
     295             : typedef long s64;
     296             : #define addmuliu64_inplace addmuliu_inplace
     297             : #define submuliu64_inplace submuliu_inplace
     298             : #define submuliu642n submuliu2n
     299             : #define addmuliu642n addmuliu2n
     300             : #else
     301             : typedef long long s64;
     302             : typedef unsigned long long u64;
     303             : 
     304             : INLINE GEN
     305    42949151 : u64toi(u64 x)
     306             : {
     307             :   GEN y;
     308             :   ulong h;
     309    42949151 :   if (!x) return gen_0;
     310    42949151 :   h = x>>32;
     311    42949151 :   if (!h) return utoipos(x);
     312     4038329 :   y = cgetipos(4);
     313     4038329 :   *int_LSW(y) = x&0xFFFFFFFF;
     314     4038329 :   *int_MSW(y) = x>>32;
     315     4038329 :   return y;
     316             : }
     317             : 
     318             : INLINE GEN
     319     3309490 : u64toineg(u64 x)
     320             : {
     321             :   GEN y;
     322             :   ulong h;
     323     3309490 :   if (!x) return gen_0;
     324     3309490 :   h = x>>32;
     325     3309490 :   if (!h) return utoineg(x);
     326     3309490 :   y = cgetineg(4);
     327     3309490 :   *int_LSW(y) = x&0xFFFFFFFF;
     328     3309490 :   *int_MSW(y) = x>>32;
     329     3309490 :   return y;
     330             : }
     331             : INLINE GEN
     332    19770843 : addmuliu64_inplace(GEN x, GEN y, u64 u) { return addmulii(x, y, u64toi(u)); }
     333             : 
     334             : INLINE GEN
     335    19926645 : submuliu64_inplace(GEN x, GEN y, u64 u) { return submulii(x, y, u64toi(u)); }
     336             : 
     337             : INLINE GEN
     338     3309490 : addmuliu642n(GEN x, GEN y, u64 u, long e) { return submulshift(x, y, u64toineg(u), e); }
     339             : 
     340             : INLINE GEN
     341     3251663 : submuliu642n(GEN x, GEN y, u64 u, long e) { return submulshift(x, y, u64toi(u), e); }
     342             : 
     343             : #endif
     344             : 
     345             : /* Babai's Nearest Plane algorithm (iterative); see Babai() */
     346             : static int
     347    40768813 : Babai_fast(pari_sp av, long kappa, GEN *pB, GEN *pU, double **mu, double **r,
     348             :            double *s, double **appB, GEN expoB, double **G,
     349             :            long a, long zeros, long maxG, double eta)
     350             : {
     351    40768813 :   GEN B = *pB, U = *pU;
     352    40768813 :   const long n = nbrows(B), d = U ? lg(U)-1: 0;
     353    40768516 :   long k, aa = (a > zeros)? a : zeros+1;
     354    40768516 :   long emaxmu = EX0, emax2mu = EX0;
     355             :   s64 xx;
     356    40768516 :   int did_something = 0;
     357             :   /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
     358             : 
     359    24518052 :   for (;;) {
     360    65286568 :     int go_on = 0;
     361    65286568 :     long i, j, emax3mu = emax2mu;
     362             : 
     363    65286568 :     if (gc_needed(av,2))
     364             :     {
     365         479 :       if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
     366         479 :       gc_lll(av,2,&B,&U);
     367             :     }
     368             :     /* Step2: compute the GSO for stage kappa */
     369    65285140 :     emax2mu = emaxmu; emaxmu = EX0;
     370   276327412 :     for (j=aa; j<kappa; j++)
     371             :     {
     372   211043208 :       double g = dmael(G,kappa,j);
     373  1119687488 :       for (k = zeros+1; k < j; k++) g -= dmael(mu,j,k) * dmael(r,kappa,k);
     374   211043208 :       dmael(r,kappa,j) = g;
     375   211043208 :       dmael(mu,kappa,j) = dmael(r,kappa,j) / dmael(r,j,j);
     376   211043208 :       emaxmu = maxss(emaxmu, expoB[kappa]-expoB[j]);
     377             :     }
     378             :     /* maxmu doesn't decrease fast enough */
     379    65284204 :     if (emax3mu != EX0 && emax3mu <= emax2mu + 5) {*pB = B; *pU = U; return 1;}
     380             : 
     381   266441535 :     for (j=kappa-1; j>zeros; j--)
     382             :     {
     383   225679101 :       double tmp = fabs(ldexp (dmael(mu,kappa,j), expoB[kappa]-expoB[j]));
     384   225679101 :       if (tmp>eta) { go_on = 1; break; }
     385             :     }
     386             : 
     387             :     /* Step3--5: compute the X_j's  */
     388    65280757 :     if (go_on)
     389   126705921 :       for (j=kappa-1; j>zeros; j--)
     390             :       { /* The code below seemingly handles U = NULL, but in this case d = 0 */
     391   102207925 :         int e = expoB[j] - expoB[kappa];
     392   102207925 :         double tmp = ldexp(dmael(mu,kappa,j), -e), atmp = fabs(tmp);
     393             :         /* tmp = Inf is allowed */
     394   102207925 :         if (atmp <= .5) continue; /* size-reduced */
     395    56909007 :         if (gc_needed(av,2))
     396             :         {
     397        1880 :           if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
     398        1880 :           gc_lll(av,2,&B,&U);
     399             :         }
     400    56911505 :         did_something = 1;
     401             :         /* we consider separately the case |X| = 1 */
     402    56911505 :         if (atmp <= 1.5)
     403             :         {
     404    37176639 :           if (dmael(mu,kappa,j) > 0) { /* in this case, X = 1 */
     405    95564236 :             for (k=zeros+1; k<j; k++)
     406    76445474 :               dmael(mu,kappa,k) -= ldexp(dmael(mu,j,k), e);
     407   347538831 :             for (i=1; i<=n; i++)
     408   328431338 :               gmael(B,kappa,i) = subii(gmael(B,kappa,i), gmael(B,j,i));
     409   120374972 :             for (i=1; i<=d; i++)
     410   101267246 :               gmael(U,kappa,i) = subii(gmael(U,kappa,i), gmael(U,j,i));
     411             :           } else { /* otherwise X = -1 */
     412    93967246 :             for (k=zeros+1; k<j; k++)
     413    75909369 :               dmael(mu,kappa,k) += ldexp(dmael(mu,j,k), e);
     414   341533678 :             for (i=1; i<=n; i++)
     415   323486036 :               gmael(B,kappa,i) = addii(gmael(B,kappa,i), gmael(B,j,i));
     416   114799660 :             for (i=1; i<=d; i++)
     417    96751880 :               gmael(U,kappa,i) = addii(gmael(U,kappa,i), gmael(U,j,i));
     418             :           }
     419    37155506 :           continue;
     420             :         }
     421             :         /* we have |X| >= 2 */
     422    19734866 :         if (atmp < 9007199254740992.)
     423             :         {
     424    17343643 :           tmp = pari_rint(tmp);
     425    44504243 :           for (k=zeros+1; k<j; k++)
     426    27160601 :             dmael(mu,kappa,k) -= ldexp(tmp * dmael(mu,j,k), e);
     427    17343642 :           xx = (s64) tmp;
     428    17343642 :           if (xx > 0) /* = xx */
     429             :           {
     430    94706934 :             for (i=1; i<=n; i++)
     431    85989449 :               gmael(B,kappa,i) = submuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
     432    62242490 :             for (i=1; i<=d; i++)
     433    53524873 :               gmael(U,kappa,i) = submuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
     434             :           }
     435             :           else /* = -xx */
     436             :           {
     437    93911514 :             for (i=1; i<=n; i++)
     438    85286889 :               gmael(B,kappa,i) = addmuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), -xx);
     439    61497995 :             for (i=1; i<=d; i++)
     440    52873217 :               gmael(U,kappa,i) = addmuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), -xx);
     441             :           }
     442             :         }
     443             :         else
     444             :         {
     445             :           int E;
     446     2391223 :           xx = (s64) ldexp(frexp(dmael(mu,kappa,j), &E), 53);
     447     2391223 :           E -= e + 53;
     448     2391223 :           if (E <= 0)
     449             :           {
     450           0 :             xx = xx << -E;
     451           0 :             for (k=zeros+1; k<j; k++)
     452           0 :               dmael(mu,kappa,k) -= ldexp(((double)xx) * dmael(mu,j,k), e);
     453           0 :             if (xx > 0) /* = xx */
     454             :             {
     455           0 :               for (i=1; i<=n; i++)
     456           0 :                 gmael(B,kappa,i) = submuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
     457           0 :               for (i=1; i<=d; i++)
     458           0 :                 gmael(U,kappa,i) = submuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
     459             :             }
     460             :             else /* = -xx */
     461             :             {
     462           0 :               for (i=1; i<=n; i++)
     463           0 :                 gmael(B,kappa,i) = addmuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), -xx);
     464           0 :               for (i=1; i<=d; i++)
     465           0 :                 gmael(U,kappa,i) = addmuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), -xx);
     466             :             }
     467             :           } else
     468             :           {
     469    15961730 :             for (k=zeros+1; k<j; k++)
     470    13570507 :               dmael(mu,kappa,k) -= ldexp(((double)xx) * dmael(mu,j,k), E + e);
     471     2391223 :             if (xx > 0) /* = xx */
     472             :             {
     473    22707755 :               for (i=1; i<=n; i++)
     474    21518924 :                 gmael(B,kappa,i) = submuliu642n(gmael(B,kappa,i), gmael(B,j,i), xx, E);
     475     2265022 :               for (i=1; i<=d; i++)
     476     1076191 :                 gmael(U,kappa,i) = submuliu642n(gmael(U,kappa,i), gmael(U,j,i), xx, E);
     477             :             }
     478             :             else /* = -xx */
     479             :             {
     480    23320863 :               for (i=1; i<=n; i++)
     481    22118830 :                 gmael(B,kappa,i) = addmuliu642n(gmael(B,kappa,i), gmael(B,j,i), -xx, E);
     482     2258739 :               for (i=1; i<=d; i++)
     483     1056706 :                 gmael(U,kappa,i) = addmuliu642n(gmael(U,kappa,i), gmael(U,j,i), -xx, E);
     484             :             }
     485             :           }
     486             :         }
     487             :       }
     488    65260515 :     if (!go_on) break; /* Anything happened? */
     489    24495513 :     expoB[kappa] = set_line(del(appB,kappa), gel(B,kappa), n);
     490    24517433 :     setG_fast(appB, n, G, kappa, zeros+1, kappa-1);
     491    24518052 :     aa = zeros+1;
     492             :   }
     493    40765002 :   if (did_something) setG2_fast(appB, n, G, kappa, kappa, maxG);
     494             : 
     495    40765471 :   del(s,zeros+1) = dmael(G,kappa,kappa);
     496             :   /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
     497   183334809 :   for (k=zeros+1; k<=kappa-2; k++)
     498   142569338 :     del(s,k+1) = del(s,k) - dmael(mu,kappa,k)*dmael(r,kappa,k);
     499    40765471 :   *pB = B; *pU = U; return 0;
     500             : }
     501             : 
     502             : static void
     503    17701252 : update_alpha(GEN alpha, long kappa, long kappa2, long kappamax)
     504             : {
     505             :   long i;
     506    58481975 :   for (i = kappa; i < kappa2; i++)
     507    40780723 :     if (kappa <= alpha[i]) alpha[i] = kappa;
     508    58482028 :   for (i = kappa2; i > kappa; i--) alpha[i] = alpha[i-1];
     509    48023003 :   for (i = kappa2+1; i <= kappamax; i++)
     510    30321751 :     if (kappa < alpha[i]) alpha[i] = kappa;
     511    17701252 :   alpha[kappa] = kappa;
     512    17701252 : }
     513             : static void
     514     1737043 : rotateG(GEN G, long kappa2, long kappa, long maxG, GEN Gtmp)
     515             : {
     516             :   long i, j;
     517    20963591 :   for (i=1; i<=kappa2; i++) gel(Gtmp,i) = gmael(G,kappa2,i);
     518    10456122 :   for (   ; i<=maxG; i++)   gel(Gtmp,i) = gmael(G,i,kappa2);
     519     7161434 :   for (i=kappa2; i>kappa; i--)
     520             :     {
     521    38112995 :       for (j=1; j<kappa; j++) gmael(G,i,j) = gmael(G,i-1,j);
     522     5424391 :       gmael(G,i,kappa) = gel(Gtmp,i-1);
     523    26054763 :       for (j=kappa+1; j<=i; j++) gmael(G,i,j) = gmael(G,i-1,j-1);
     524    26994555 :       for (j=kappa2+1; j<=maxG; j++) gmael(G,j,i) = gmael(G,j,i-1);
     525             :     }
     526    13802157 :   for (i=1; i<kappa; i++) gmael(G,kappa,i) = gel(Gtmp,i);
     527     1737043 :   gmael(G,kappa,kappa) = gel(Gtmp,kappa2);
     528    10456122 :   for (i=kappa2+1; i<=maxG; i++) gmael(G,i,kappa) = gel(Gtmp,i);
     529     1737043 : }
     530             : static void
     531    15963995 : rotateG_fast(double **G, long kappa2, long kappa, long maxG, double *Gtmp)
     532             : {
     533             :   long i, j;
     534   101898432 :   for (i=1; i<=kappa2; i++) del(Gtmp,i) = dmael(G,kappa2,i);
     535    38551532 :   for (   ; i<=maxG; i++) del(Gtmp,i) = dmael(G,i,kappa2);
     536    51319961 :   for (i=kappa2; i>kappa; i--)
     537             :   {
     538   121033075 :     for (j=1; j<kappa; j++) dmael(G,i,j) = dmael(G,i-1,j);
     539    35355966 :     dmael(G,i,kappa) = del(Gtmp,i-1);
     540   128703744 :     for (j=kappa+1; j<=i; j++) dmael(G,i,j) = dmael(G,i-1,j-1);
     541    86313214 :     for (j=kappa2+1; j<=maxG; j++) dmael(G,j,i) = dmael(G,j,i-1);
     542             :   }
     543    50579098 :   for (i=1; i<kappa; i++) dmael(G,kappa,i) = del(Gtmp,i);
     544    15963995 :   dmael(G,kappa,kappa) = del(Gtmp,kappa2);
     545    38551602 :   for (i=kappa2+1; i<=maxG; i++) dmael(G,i,kappa) = del(Gtmp,i);
     546    15963995 : }
     547             : 
     548             : /* LLL-reduces (B,U) in place [apply base change transforms to B and U].
     549             :  * Gram matrix, and GSO performed on matrices of 'double'.
     550             :  * If (keepfirst), never swap with first vector.
     551             :  * Return -1 on failure, else zeros = dim Kernel (>= 0) */
     552             : static long
     553     2327010 : fplll_fast(GEN *pB, GEN *pU, double delta, double eta, long keepfirst)
     554             : {
     555             :   pari_sp av;
     556             :   long kappa, kappa2, d, n, i, j, zeros, kappamax, maxG;
     557             :   double **mu, **r, *s, tmp, *Gtmp, **G, **appB;
     558     2327010 :   GEN alpha, expoB, B = *pB, U;
     559     2327010 :   long cnt = 0;
     560             : 
     561     2327010 :   d = lg(B)-1;
     562     2327010 :   n = nbrows(B);
     563     2327010 :   U = *pU; /* NULL if inplace */
     564             : 
     565     2327010 :   G = cget_dblmat(d+1);
     566     2327000 :   appB = cget_dblmat(d+1);
     567     2327000 :   mu = cget_dblmat(d+1);
     568     2326995 :   r  = cget_dblmat(d+1);
     569     2326999 :   s  = cget_dblvec(d+1);
     570    10423047 :   for (j = 1; j <= d; j++)
     571             :   {
     572     8096043 :     del(mu,j) = cget_dblvec(d+1);
     573     8096045 :     del(r,j) = cget_dblvec(d+1);
     574     8096068 :     del(appB,j) = cget_dblvec(n+1);
     575     8096048 :     del(G,j) = cget_dblvec(d+1);
     576             :   }
     577     2327004 :   expoB = cgetg(d+1, t_VECSMALL);
     578    10422920 :   for (i=1; i<=d; i++) expoB[i] = set_line(del(appB,i), gel(B,i), n);
     579     2326938 :   Gtmp = cget_dblvec(d+1);
     580     2326993 :   alpha = cgetg(d+1, t_VECSMALL);
     581     2326983 :   av = avma;
     582             : 
     583             :   /* Step2: Initializing the main loop */
     584     2326983 :   kappamax = 1;
     585     2326983 :   i = 1;
     586     2326983 :   maxG = d; /* later updated to kappamax */
     587             : 
     588             :   do {
     589     2480046 :     dmael(G,i,i) = dbldotsquare(del(appB,i),n);
     590     2480053 :   } while (dmael(G,i,i) <= 0 && (++i <=d));
     591     2326990 :   zeros = i-1; /* all vectors B[i] with i <= zeros are zero vectors */
     592     2326990 :   kappa = i;
     593     2326990 :   if (zeros < d) dmael(r,zeros+1,zeros+1) = dmael(G,zeros+1,zeros+1);
     594    10269983 :   for (i=zeros+1; i<=d; i++) alpha[i]=1;
     595    43092596 :   while (++kappa <= d)
     596             :   {
     597    40769016 :     if (kappa > kappamax)
     598             :     {
     599     5609776 :       if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
     600     5609776 :       maxG = kappamax = kappa;
     601     5609776 :       setG_fast(appB, n, G, kappa, zeros+1, kappa);
     602             :     }
     603             :     /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
     604    40769009 :     if (Babai_fast(av, kappa, &B,&U, mu,r,s, appB, expoB, G, alpha[kappa],
     605        3447 :                    zeros, maxG, eta)) { *pB=B; *pU=U; return -1; }
     606             : 
     607    40765362 :     tmp = ldexp(r[kappa-1][kappa-1] * delta, 2*(expoB[kappa-1]-expoB[kappa]));
     608    40765362 :     if ((keepfirst && kappa == 2) || tmp <= del(s,kappa-1))
     609             :     { /* Step4: Success of Lovasz's condition */
     610    24801166 :       alpha[kappa] = kappa;
     611    24801166 :       tmp = dmael(mu,kappa,kappa-1) * dmael(r,kappa,kappa-1);
     612    24801166 :       dmael(r,kappa,kappa) = del(s,kappa-1)- tmp;
     613    24801166 :       continue;
     614             :     }
     615             :     /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
     616    15964196 :     if (DEBUGLEVEL>=4 && kappa==kappamax && del(s,kappa-1)!=0)
     617           0 :       if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", 2*expoB[1] + dblexpo(del(s,1))); }
     618    15964191 :     kappa2 = kappa;
     619             :     do {
     620    35356262 :       kappa--;
     621    35356262 :       if (kappa<zeros+2 + (keepfirst ? 1: 0)) break;
     622    26541908 :       tmp = dmael(r,kappa-1,kappa-1) * delta;
     623    26541908 :       tmp = ldexp(tmp, 2*(expoB[kappa-1]-expoB[kappa2]));
     624    26541908 :     } while (del(s,kappa-1) <= tmp);
     625    15964191 :     update_alpha(alpha, kappa, kappa2, kappamax);
     626             : 
     627             :     /* Step6: Update the mu's and r's */
     628    15964241 :     dblrotate(mu,kappa2,kappa);
     629    15964199 :     dblrotate(r,kappa2,kappa);
     630    15964152 :     dmael(r,kappa,kappa) = del(s,kappa);
     631             : 
     632             :     /* Step7: Update B, appB, U, G */
     633    15964152 :     rotate(B,kappa2,kappa);
     634    15964096 :     dblrotate(appB,kappa2,kappa);
     635    15964057 :     if (U) rotate(U,kappa2,kappa);
     636    15964050 :     rotate(expoB,kappa2,kappa);
     637    15964005 :     rotateG_fast(G,kappa2,kappa, maxG, Gtmp);
     638             : 
     639             :     /* Step8: Prepare the next loop iteration */
     640    15964440 :     if (kappa == zeros+1 && dmael(G,kappa,kappa)<= 0)
     641             :     {
     642      194770 :       zeros++; kappa++;
     643      194770 :       dmael(G,kappa,kappa) = dbldotsquare(del(appB,kappa),n);
     644      194770 :       dmael(r,kappa,kappa) = dmael(G,kappa,kappa);
     645             :     }
     646             :   }
     647     2323580 :   *pB = B; *pU = U; return zeros;
     648             : }
     649             : 
     650             : /***************** HEURISTIC version (reduced precision) ****************/
     651             : static GEN
     652      340849 : realsqrdotproduct(GEN x)
     653             : {
     654      340849 :   long i, l = lg(x);
     655      340849 :   GEN z = sqrr(gel(x,1));
     656     4236909 :   for (i=2; i<l; i++) z = addrr(z, sqrr(gel(x,i)));
     657      340849 :   return z;
     658             : }
     659             : /* x, y non-empty vector of t_REALs, same length */
     660             : static GEN
     661     3855818 : realdotproduct(GEN x, GEN y)
     662             : {
     663             :   long i, l;
     664             :   GEN z;
     665     3855818 :   if (x == y) return realsqrdotproduct(x);
     666     3514969 :   l = lg(x); z = mulrr(gel(x,1),gel(y,1));
     667    57087233 :   for (i=2; i<l; i++) z = addrr(z, mulrr(gel(x,i), gel(y,i)));
     668     3514969 :   return z;
     669             : }
     670             : static void
     671      348202 : setG_heuristic(GEN appB, GEN G, long kappa, long a, long b)
     672      348202 : { pari_sp av = avma;
     673             :   long i;
     674     2785548 :   for (i = a; i <= b; i++)
     675     2437346 :     affrr(realdotproduct(gel(appB,kappa),gel(appB,i)), gmael(G,kappa,i));
     676      348202 :   set_avma(av);
     677      348202 : }
     678             : static void
     679      329719 : setG2_heuristic(GEN appB, GEN G, long kappa, long a, long b)
     680      329719 : { pari_sp av = avma;
     681             :   long i;
     682     1748191 :   for (i = a; i <= b; i++)
     683     1418472 :     affrr(realdotproduct(gel(appB,kappa),gel(appB,i)), gmael(G,i,kappa));
     684      329719 :   set_avma(av);
     685      329719 : }
     686             : 
     687             : /* approximate t_REAL x as m * 2^e, where |m| < 2^bit */
     688             : static GEN
     689       21060 : truncexpo(GEN x, long bit, long *e)
     690             : {
     691       21060 :   *e = expo(x) + 1 - bit;
     692       21060 :   if (*e >= 0) return mantissa2nr(x, 0);
     693        1870 :   *e = 0; return roundr_safe(x);
     694             : }
     695             : /* Babai's Nearest Plane algorithm (iterative); see Babai() */
     696             : static int
     697      677423 : Babai_heuristic(pari_sp av, long kappa, GEN *pB, GEN *pU, GEN mu, GEN r, GEN s,
     698             :                 GEN appB, GEN G, long a, long zeros, long maxG,
     699             :                 GEN eta, long prec)
     700             : {
     701      677423 :   GEN B = *pB, U = *pU;
     702      677423 :   const long n = nbrows(B), d = U ? lg(U)-1: 0, bit = prec2nbits(prec);
     703      677423 :   long k, aa = (a > zeros)? a : zeros+1;
     704      677423 :   int did_something = 0;
     705      677423 :   long emaxmu = EX0, emax2mu = EX0;
     706             :   /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
     707             : 
     708      337072 :   for (;;) {
     709     1014495 :     int go_on = 0;
     710     1014495 :     long i, j, emax3mu = emax2mu;
     711             : 
     712     1014495 :     if (gc_needed(av,2))
     713             :     {
     714          14 :       if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
     715          14 :       gc_lll(av,2,&B,&U);
     716             :     }
     717             :     /* Step2: compute the GSO for stage kappa */
     718     1014495 :     emax2mu = emaxmu; emaxmu = EX0;
     719     5723022 :     for (j=aa; j<kappa; j++)
     720             :     {
     721     4708527 :       pari_sp btop = avma;
     722     4708527 :       GEN g = gmael(G,kappa,j);
     723    30245965 :       for (k = zeros+1; k<j; k++)
     724    25537438 :         g = subrr(g, mulrr(gmael(mu,j,k), gmael(r,kappa,k)));
     725     4708527 :       affrr(g, gmael(r,kappa,j));
     726     4708527 :       affrr(divrr(gmael(r,kappa,j), gmael(r,j,j)), gmael(mu,kappa,j));
     727     4708527 :       emaxmu = maxss(emaxmu, expo(gmael(mu,kappa,j)));
     728     4708527 :       set_avma(btop);
     729             :     }
     730     1014495 :     if (emax3mu != EX0 && emax3mu <= emax2mu + 5)
     731        1146 :     { *pB = B; *pU = U; return 1; }
     732             : 
     733     6250095 :     for (j=kappa-1; j>zeros; j--)
     734     5573818 :       if (abscmprr(gmael(mu,kappa,j), eta) > 0) { go_on = 1; break; }
     735             : 
     736             :     /* Step3--5: compute the X_j's  */
     737     1013349 :     if (go_on)
     738     2726189 :       for (j=kappa-1; j>zeros; j--)
     739             :       { /* The code below seemingly handles U = NULL, but in this case d = 0 */
     740             :         pari_sp btop;
     741     2389117 :         GEN tmp = gmael(mu,kappa,j);
     742     2389117 :         if (absrsmall(tmp)) continue; /* size-reduced */
     743             : 
     744     1086567 :         if (gc_needed(av,2))
     745             :         {
     746          84 :           if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
     747          84 :           gc_lll(av,2,&B,&U);
     748             :         }
     749     1086567 :         btop = avma; did_something = 1;
     750             :         /* we consider separately the case |X| = 1 */
     751     1086567 :         if (absrsmall2(tmp))
     752             :         {
     753      897647 :           if (signe(tmp) > 0) { /* in this case, X = 1 */
     754     2582780 :             for (k=zeros+1; k<j; k++)
     755     2134607 :               affrr(subrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
     756      448173 :             set_avma(btop);
     757     7894422 :             for (i=1; i<=n; i++)
     758     7446249 :               gmael(B,kappa,i) = subii(gmael(B,kappa,i), gmael(B,j,i));
     759     6676484 :             for (i=1; i<=d; i++)
     760     6228311 :               gmael(U,kappa,i) = subii(gmael(U,kappa,i), gmael(U,j,i));
     761             :           } else { /* otherwise X = -1 */
     762     2586266 :             for (k=zeros+1; k<j; k++)
     763     2136792 :               affrr(addrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
     764      449474 :             set_avma(btop);
     765     7921890 :             for (i=1; i<=n; i++)
     766     7472416 :               gmael(B,kappa,i) = addii(gmael(B,kappa,i), gmael(B,j,i));
     767     6700869 :             for (i=1; i<=d; i++)
     768     6251395 :               gmael(U,kappa,i) = addii(gmael(U,kappa,i),gmael(U,j,i));
     769             :           }
     770      897647 :           continue;
     771             :         }
     772             :         /* we have |X| >= 2 */
     773      188920 :         if (expo(tmp) < BITS_IN_LONG)
     774             :         {
     775      168380 :           ulong xx = roundr_safe(tmp)[2]; /* X fits in an ulong */
     776      168380 :           if (signe(tmp) > 0) /* = xx */
     777             :           {
     778      307368 :             for (k=zeros+1; k<j; k++)
     779      223569 :               affrr(subrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
     780      223569 :                   gmael(mu,kappa,k));
     781       83799 :             set_avma(btop);
     782     1110001 :             for (i=1; i<=n; i++)
     783     1026202 :               gmael(B,kappa,i) = submuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
     784      919889 :             for (i=1; i<=d; i++)
     785      836090 :               gmael(U,kappa,i) = submuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
     786             :           }
     787             :           else /* = -xx */
     788             :           {
     789      309539 :             for (k=zeros+1; k<j; k++)
     790      224958 :               affrr(addrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
     791      224958 :                   gmael(mu,kappa,k));
     792       84581 :             set_avma(btop);
     793     1113565 :             for (i=1; i<=n; i++)
     794     1028984 :               gmael(B,kappa,i) = addmuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
     795      915355 :             for (i=1; i<=d; i++)
     796      830774 :               gmael(U,kappa,i) = addmuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
     797             :           }
     798             :         }
     799             :         else
     800             :         {
     801             :           long e;
     802       20540 :           GEN X = truncexpo(tmp, bit, &e); /* tmp ~ X * 2^e */
     803       20540 :           btop = avma;
     804       88630 :           for (k=zeros+1; k<j; k++)
     805             :           {
     806       68090 :             GEN x = mulir(X, gmael(mu,j,k));
     807       68090 :             if (e) shiftr_inplace(x, e);
     808       68090 :             affrr(subrr(gmael(mu,kappa,k), x), gmael(mu,kappa,k));
     809             :           }
     810       20540 :           set_avma(btop);
     811      265641 :           for (i=1; i<=n; i++)
     812      245101 :             gmael(B,kappa,i) = submulshift(gmael(B,kappa,i), gmael(B,j,i), X, e);
     813       98271 :           for (i=1; i<=d; i++)
     814       77731 :             gmael(U,kappa,i) = submulshift(gmael(U,kappa,i), gmael(U,j,i), X, e);
     815             :         }
     816             :       }
     817     1013349 :     if (!go_on) break; /* Anything happened? */
     818     4546534 :     for (i=1 ; i<=n; i++) affir(gmael(B,kappa,i), gmael(appB,kappa,i));
     819      337072 :     setG_heuristic(appB, G, kappa, zeros+1, kappa-1);
     820      337072 :     aa = zeros+1;
     821             :   }
     822      676277 :   if (did_something) setG2_heuristic(appB, G, kappa, kappa, maxG);
     823      676277 :   affrr(gmael(G,kappa,kappa), gel(s,zeros+1));
     824             :   /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
     825      676277 :   av = avma;
     826     4790819 :   for (k=zeros+1; k<=kappa-2; k++)
     827     4114542 :     affrr(subrr(gel(s,k), mulrr(gmael(mu,kappa,k), gmael(r,kappa,k))),
     828     4114542 :           gel(s,k+1));
     829      676277 :   *pB = B; *pU = U; return gc_bool(av, 0);
     830             : }
     831             : 
     832             : static GEN
     833       17816 : ZC_to_RC(GEN x, long prec)
     834      177243 : { pari_APPLY_type(t_COL,itor(gel(x,i),prec)) }
     835             : 
     836             : static GEN
     837        3447 : ZM_to_RM(GEN x, long prec)
     838       21263 : { pari_APPLY_same(ZC_to_RC(gel(x,i),prec)) }
     839             : 
     840             : /* LLL-reduces (B,U) in place [apply base change transforms to B and U].
     841             :  * Gram matrix made of t_REAL at precision prec2, performe GSO at prec.
     842             :  * If (keepfirst), never swap with first vector.
     843             :  * Return -1 on failure, else zeros = dim Kernel (>= 0) */
     844             : static long
     845        3447 : fplll_heuristic(GEN *pB, GEN *pU, double DELTA, double ETA, long keepfirst,
     846             :                 long prec, long prec2)
     847             : {
     848             :   pari_sp av, av2;
     849             :   long kappa, kappa2, d, i, j, zeros, kappamax, maxG;
     850        3447 :   GEN mu, r, s, tmp, Gtmp, alpha, G, appB, B = *pB, U;
     851        3447 :   GEN delta = dbltor(DELTA), eta = dbltor(ETA);
     852        3447 :   long cnt = 0;
     853             : 
     854        3447 :   d = lg(B)-1;
     855        3447 :   U = *pU; /* NULL if inplace */
     856             : 
     857        3447 :   G = cgetg(d+1, t_MAT);
     858        3447 :   mu = cgetg(d+1, t_MAT);
     859        3447 :   r  = cgetg(d+1, t_MAT);
     860        3447 :   s  = cgetg(d+1, t_VEC);
     861        3447 :   appB = ZM_to_RM(B, prec2);
     862       21263 :   for (j = 1; j <= d; j++)
     863             :   {
     864       17816 :     GEN M = cgetg(d+1, t_COL), R = cgetg(d+1, t_COL), S = cgetg(d+1, t_COL);
     865       17816 :     gel(mu,j)= M;
     866       17816 :     gel(r,j) = R;
     867       17816 :     gel(G,j) = S;
     868       17816 :     gel(s,j) = cgetr(prec);
     869      168482 :     for (i = 1; i <= d; i++)
     870             :     {
     871      150666 :       gel(R,i) = cgetr(prec);
     872      150666 :       gel(M,i) = cgetr(prec);
     873      150666 :       gel(S,i) = cgetr(prec2);
     874             :     }
     875             :   }
     876        3447 :   Gtmp = cgetg(d+1, t_VEC);
     877        3447 :   alpha = cgetg(d+1, t_VECSMALL);
     878        3447 :   av = avma;
     879             : 
     880             :   /* Step2: Initializing the main loop */
     881        3447 :   kappamax = 1;
     882        3447 :   i = 1;
     883        3447 :   maxG = d; /* later updated to kappamax */
     884             : 
     885             :   do {
     886        3447 :     affrr(RgV_dotsquare(gel(appB,i)), gmael(G,i,i));
     887        3447 :   } while (signe(gmael(G,i,i)) == 0 && (++i <=d));
     888        3447 :   zeros = i-1; /* all vectors B[i] with i <= zeros are zero vectors */
     889        3447 :   kappa = i;
     890        3447 :   if (zeros < d) affrr(gmael(G,zeros+1,zeros+1), gmael(r,zeros+1,zeros+1));
     891       21263 :   for (i=zeros+1; i<=d; i++) alpha[i]=1;
     892             : 
     893      679724 :   while (++kappa <= d)
     894             :   {
     895      677423 :     if (kappa > kappamax)
     896             :     {
     897       11130 :       if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
     898       11130 :       maxG = kappamax = kappa;
     899       11130 :       setG_heuristic(appB, G, kappa, zeros+1, kappa);
     900             :     }
     901             :     /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
     902      677423 :     if (Babai_heuristic(av, kappa, &B,&U, mu,r,s, appB, G, alpha[kappa], zeros,
     903        1146 :                         maxG, eta, prec)) { *pB = B; *pU = U; return -1; }
     904      676277 :     av2 = avma;
     905     1352462 :     if ((keepfirst && kappa == 2) ||
     906      676185 :         cmprr(mulrr(gmael(r,kappa-1,kappa-1), delta), gel(s,kappa-1)) <= 0)
     907             :     { /* Step4: Success of Lovasz's condition */
     908      456943 :       alpha[kappa] = kappa;
     909      456943 :       tmp = mulrr(gmael(mu,kappa,kappa-1), gmael(r,kappa,kappa-1));
     910      456943 :       affrr(subrr(gel(s,kappa-1), tmp), gmael(r,kappa,kappa));
     911      456943 :       set_avma(av2); continue;
     912             :     }
     913             :     /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
     914      219334 :     if (DEBUGLEVEL>=4 && kappa==kappamax && signe(gel(s,kappa-1)))
     915           0 :       if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", expo(gel(s,1))); }
     916      219334 :     kappa2 = kappa;
     917             :     do {
     918      666320 :       kappa--;
     919      666320 :       if (kappa < zeros+2 + (keepfirst ? 1: 0)) break;
     920      636368 :       tmp = mulrr(gmael(r,kappa-1,kappa-1), delta);
     921      636368 :     } while (cmprr(gel(s,kappa-1), tmp) <= 0 );
     922      219334 :     set_avma(av2);
     923      219334 :     update_alpha(alpha, kappa, kappa2, kappamax);
     924             : 
     925             :     /* Step6: Update the mu's and r's */
     926      219334 :     rotate(mu,kappa2,kappa);
     927      219334 :     rotate(r,kappa2,kappa);
     928      219334 :     affrr(gel(s,kappa), gmael(r,kappa,kappa));
     929             : 
     930             :     /* Step7: Update B, appB, U, G */
     931      219334 :     rotate(B,kappa2,kappa);
     932      219334 :     rotate(appB,kappa2,kappa);
     933      219334 :     if (U) rotate(U,kappa2,kappa);
     934      219334 :     rotateG(G,kappa2,kappa, maxG, Gtmp);
     935             : 
     936             :     /* Step8: Prepare the next loop iteration */
     937      219334 :     if (kappa == zeros+1 && !signe(gmael(G,kappa,kappa)))
     938             :     {
     939           6 :       zeros++; kappa++;
     940           6 :       affrr(RgV_dotsquare(gel(appB,kappa)), gmael(G,kappa,kappa));
     941           6 :       affrr(gmael(G,kappa,kappa), gmael(r,kappa,kappa));
     942             :     }
     943             :   }
     944        2301 :   *pB=B; *pU=U; return zeros;
     945             : }
     946             : 
     947             : /************************* PROVED version (t_INT) ***********************/
     948             : /* dpe inspired by dpe.h by Patrick Pelissier, Paul Zimmermann
     949             :  * https://gforge.inria.fr/projects/dpe/
     950             :  */
     951             : 
     952             : typedef struct
     953             : {
     954             :   double d;  /* significand */
     955             :   long e; /* exponent */
     956             : } dpe_t;
     957             : 
     958             : #define Dmael(x,i,j) (&((x)[i][j]))
     959             : #define Del(x,i) (&((x)[i]))
     960             : 
     961             : static void
     962     3033930 : dperotate(dpe_t **A, long k2, long k)
     963             : {
     964             :   long i;
     965     3033930 :   dpe_t *B = A[k2];
     966    12545734 :   for (i = k2; i > k; i--) A[i] = A[i-1];
     967     3033930 :   A[k] = B;
     968     3033930 : }
     969             : 
     970             : static void
     971   859035435 : dpe_normalize0(dpe_t *x)
     972             : {
     973             :   int e;
     974   859035435 :   x->d = frexp(x->d, &e);
     975   859035435 :   x->e += e;
     976   859035435 : }
     977             : 
     978             : static void
     979   423805195 : dpe_normalize(dpe_t *x)
     980             : {
     981   423805195 :   if (x->d == 0.0)
     982     2301693 :     x->e = -LONG_MAX;
     983             :   else
     984   421503502 :     dpe_normalize0(x);
     985   423805335 : }
     986             : 
     987             : static GEN
     988       83121 : dpetor(dpe_t *x)
     989             : {
     990       83121 :   GEN r = dbltor(x->d);
     991       83121 :   if (signe(r)==0) return r;
     992       83121 :   setexpo(r, x->e-1);
     993       83121 :   return r;
     994             : }
     995             : 
     996             : static void
     997    96288774 : affdpe(dpe_t *y, dpe_t *x)
     998             : {
     999    96288774 :   x->d = y->d;
    1000    96288774 :   x->e = y->e;
    1001    96288774 : }
    1002             : 
    1003             : static void
    1004    67843691 : affidpe(GEN y, dpe_t *x)
    1005             : {
    1006    67843691 :   pari_sp av = avma;
    1007    67843691 :   GEN r = itor(y, DEFAULTPREC);
    1008    67842170 :   x->e = expo(r)+1;
    1009    67842170 :   setexpo(r,-1);
    1010    67842052 :   x->d = rtodbl(r);
    1011    67842831 :   set_avma(av);
    1012    67842996 : }
    1013             : 
    1014             : static void
    1015     4818148 : affdbldpe(double y, dpe_t *x)
    1016             : {
    1017     4818148 :   x->d = (double)y;
    1018     4818148 :   x->e = 0;
    1019     4818148 :   dpe_normalize(x);
    1020     4818150 : }
    1021             : 
    1022             : static void
    1023   408028197 : dpe_mulz(dpe_t *x, dpe_t *y, dpe_t *z)
    1024             : {
    1025   408028197 :   z->d = x->d * y->d;
    1026   408028197 :   if (z->d == 0.0)
    1027    22777805 :     z->e = -LONG_MAX;
    1028             :   else
    1029             :   {
    1030   385250392 :     z->e = x->e + y->e;
    1031   385250392 :     dpe_normalize0(z);
    1032             :   }
    1033   408028567 : }
    1034             : 
    1035             : static void
    1036    54674415 : dpe_divz(dpe_t *x, dpe_t *y, dpe_t *z)
    1037             : {
    1038    54674415 :   z->d = x->d / y->d;
    1039    54674415 :   if (z->d == 0.0)
    1040     2392122 :     z->e = -LONG_MAX;
    1041             :   else
    1042             :   {
    1043    52282293 :     z->e = x->e - y->e;
    1044    52282293 :     dpe_normalize0(z);
    1045             :   }
    1046    54674506 : }
    1047             : 
    1048             : static void
    1049      574136 : dpe_negz(dpe_t *y, dpe_t *x)
    1050             : {
    1051      574136 :   x->d = - y->d;
    1052      574136 :   x->e = y->e;
    1053      574136 : }
    1054             : 
    1055             : static void
    1056    30111819 : dpe_addz(dpe_t *y, dpe_t *z, dpe_t *x)
    1057             : {
    1058    30111819 :   if (y->e > z->e + 53)
    1059     1062345 :     affdpe(y, x);
    1060    29049474 :   else if (z->e > y->e + 53)
    1061      274741 :     affdpe(z, x);
    1062             :   else
    1063             :   {
    1064    28774733 :     long d = y->e - z->e;
    1065             : 
    1066    28774733 :     if (d >= 0)
    1067             :     {
    1068    21361901 :       x->d = y->d + ldexp(z->d, -d);
    1069    21361901 :       x->e  = y->e;
    1070             :     }
    1071             :     else
    1072             :     {
    1073     7412832 :       x->d = z->d + ldexp(y->d, d);
    1074     7412832 :       x->e = z->e;
    1075             :     }
    1076    28774733 :     dpe_normalize(x);
    1077             :   }
    1078    30111826 : }
    1079             : static void
    1080   423634128 : dpe_subz(dpe_t *y, dpe_t *z, dpe_t *x)
    1081             : {
    1082   423634128 :   if (y->e > z->e + 53)
    1083    38760506 :     affdpe(y, x);
    1084   384873622 :   else if (z->e > y->e + 53)
    1085      574136 :     dpe_negz(z, x);
    1086             :   else
    1087             :   {
    1088   384299486 :     long d = y->e - z->e;
    1089             : 
    1090   384299486 :     if (d >= 0)
    1091             :     {
    1092   351260835 :       x->d = y->d - ldexp(z->d, -d);
    1093   351260835 :       x->e = y->e;
    1094             :     }
    1095             :     else
    1096             :     {
    1097    33038651 :       x->d = ldexp(y->d, d) - z->d;
    1098    33038651 :       x->e = z->e;
    1099             :     }
    1100   384299486 :     dpe_normalize(x);
    1101             :   }
    1102   423634549 : }
    1103             : 
    1104             : static void
    1105     5912638 : dpe_muluz(dpe_t *y, ulong t, dpe_t *x)
    1106             : {
    1107     5912638 :   x->d = y->d * (double)t;
    1108     5912638 :   x->e = y->e;
    1109     5912638 :   dpe_normalize(x);
    1110     5912638 : }
    1111             : 
    1112             : static void
    1113     2301191 : dpe_addmuluz(dpe_t *y,  dpe_t *z, ulong t, dpe_t *x)
    1114             : {
    1115             :   dpe_t tmp;
    1116     2301191 :   dpe_muluz(z, t, &tmp);
    1117     2301191 :   dpe_addz(y, &tmp, x);
    1118     2301191 : }
    1119             : 
    1120             : static void
    1121     2430025 : dpe_submuluz(dpe_t *y,  dpe_t *z, ulong t, dpe_t *x)
    1122             : {
    1123             :   dpe_t tmp;
    1124     2430025 :   dpe_muluz(z, t, &tmp);
    1125     2430025 :   dpe_subz(y, &tmp, x);
    1126     2430025 : }
    1127             : 
    1128             : static void
    1129   393255644 : dpe_submulz(dpe_t *y,  dpe_t *z, dpe_t *t, dpe_t *x)
    1130             : {
    1131             :   dpe_t tmp;
    1132   393255644 :   dpe_mulz(z, t, &tmp);
    1133   393255656 :   dpe_subz(y, &tmp, x);
    1134   393255688 : }
    1135             : 
    1136             : static int
    1137    14772724 : dpe_cmp(dpe_t *x, dpe_t *y)
    1138             : {
    1139    14772724 :   int sx = x->d < 0. ? -1: x->d > 0.;
    1140    14772724 :   int sy = y->d < 0. ? -1: y->d > 0.;
    1141    14772724 :   int d  = sx - sy;
    1142             : 
    1143    14772724 :   if (d != 0)
    1144      142833 :     return d;
    1145    14629891 :   else if (x->e > y->e)
    1146     2335819 :     return (sx > 0) ? 1 : -1;
    1147    12294072 :   else if (y->e > x->e)
    1148     4319417 :     return (sx > 0) ? -1 : 1;
    1149             :   else
    1150     7974655 :     return (x->d < y->d) ? -1 : (x->d > y->d);
    1151             : }
    1152             : 
    1153             : static int
    1154    66774454 : dpe_abscmp(dpe_t *x, dpe_t *y)
    1155             : {
    1156    66774454 :   if (x->e > y->e)
    1157      448895 :     return 1;
    1158    66325559 :   else if (y->e > x->e)
    1159    63279381 :     return -1;
    1160             :   else
    1161     3046178 :     return (fabs(x->d) < fabs(y->d)) ? -1 : (fabs(x->d) > fabs(y->d));
    1162             : }
    1163             : 
    1164             : static int
    1165     9664545 : dpe_abssmall(dpe_t *x)
    1166             : {
    1167     9664545 :   return (x->e <= 0) || (x->e == 1 && fabs(x->d) <= .75);
    1168             : }
    1169             : 
    1170             : static int
    1171    14772717 : dpe_cmpmul(dpe_t *x, dpe_t *y, dpe_t *z)
    1172             : {
    1173             :   dpe_t t;
    1174    14772717 :   dpe_mulz(x,y,&t);
    1175    14772718 :   return dpe_cmp(&t, z);
    1176             : }
    1177             : 
    1178             : static dpe_t *
    1179    19236844 : cget_dpevec(long d)
    1180    19236844 : { return (dpe_t*) stack_malloc_align(d*sizeof(dpe_t), sizeof(dpe_t)); }
    1181             : 
    1182             : static dpe_t **
    1183     4818146 : cget_dpemat(long d) { return (dpe_t **) cgetg(d, t_VECSMALL); }
    1184             : 
    1185             : static GEN
    1186       17373 : dpeM_diagonal_shallow(dpe_t **m, long d)
    1187             : {
    1188             :   long i;
    1189       17373 :   GEN y = cgetg(d+1,t_VEC);
    1190      100494 :   for (i=1; i<=d; i++) gel(y, i) = dpetor(Dmael(m,i,i));
    1191       17373 :   return y;
    1192             : }
    1193             : 
    1194             : /************************** PROVED version (dpe) *************************/
    1195             : 
    1196             : /* Babai's Nearest Plane algorithm (iterative).
    1197             :  * Size-reduces b_kappa using mu_{i,j} and r_{i,j} for j<=i <kappa
    1198             :  * Update B[,kappa]; compute mu_{kappa,j}, r_{kappa,j} for j<=kappa and s[kappa]
    1199             :  * mu, r, s updated in place (affrr). Return 1 on failure, else 0. */
    1200             : static int
    1201    10372356 : Babai_dpe(pari_sp av, long kappa, GEN *pG, GEN *pB, GEN *pU, dpe_t **mu, dpe_t **r, dpe_t *s,
    1202             :       long a, long zeros, long maxG, dpe_t *eta)
    1203             : {
    1204    10372356 :   GEN G = *pG, B = *pB, U = *pU, ztmp;
    1205    10372356 :   long k, d, n, aa = a > zeros? a: zeros+1;
    1206    10372356 :   long emaxmu = EX0, emax2mu = EX0;
    1207             :   /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
    1208    10372356 :   d = U? lg(U)-1: 0;
    1209    10372356 :   n = B? nbrows(B): 0;
    1210     2165723 :   for (;;) {
    1211    12538124 :     int go_on = 0;
    1212    12538124 :     long i, j, emax3mu = emax2mu;
    1213             : 
    1214    12538124 :     if (gc_needed(av,2))
    1215             :     {
    1216         305 :       if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
    1217         305 :       gc_lll(av,3,&G,&B,&U);
    1218             :     }
    1219             :     /* Step2: compute the GSO for stage kappa */
    1220    12538087 :     emax2mu = emaxmu; emaxmu = EX0;
    1221    67212539 :     for (j=aa; j<kappa; j++)
    1222             :     {
    1223             :       dpe_t g;
    1224    54674455 :       affidpe(gmael(G,kappa,j), &g);
    1225   388233397 :       for (k = zeros+1; k < j; k++)
    1226   333559119 :         dpe_submulz(&g, Dmael(mu,j,k), Dmael(r,kappa,k), &g);
    1227    54674278 :       affdpe(&g, Dmael(r,kappa,j));
    1228    54674404 :       dpe_divz(Dmael(r,kappa,j), Dmael(r,j,j), Dmael(mu,kappa,j));
    1229    54674415 :       emaxmu = maxss(emaxmu, Dmael(mu,kappa,j)->e);
    1230             :     }
    1231    12538084 :     if (emax3mu != EX0 && emax3mu <= emax2mu + 5) /* precision too low */
    1232           2 :     { *pG = G; *pB = B; *pU = U; return 1; }
    1233             : 
    1234    77146818 :     for (j=kappa-1; j>zeros; j--)
    1235    66774456 :       if (dpe_abscmp(Dmael(mu,kappa,j), eta) > 0) { go_on = 1; break; }
    1236             : 
    1237             :     /* Step3--5: compute the X_j's  */
    1238    12538080 :     if (go_on)
    1239    23430840 :       for (j=kappa-1; j>zeros; j--)
    1240             :       {
    1241             :         pari_sp btop;
    1242    21265113 :         dpe_t *tmp = Dmael(mu,kappa,j);
    1243    21265113 :         if (tmp->e < 0) continue; /* (essentially) size-reduced */
    1244             : 
    1245     9664546 :         if (gc_needed(av,2))
    1246             :         {
    1247        1317 :           if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
    1248        1317 :           gc_lll(av,3,&G,&B,&U);
    1249             :         }
    1250             :         /* we consider separately the case |X| = 1 */
    1251     9664546 :         if (dpe_abssmall(tmp))
    1252             :         {
    1253     8226480 :           if (tmp->d > 0) { /* in this case, X = 1 */
    1254    31459342 :             for (k=zeros+1; k<j; k++)
    1255    27338692 :               dpe_subz(Dmael(mu,kappa,k), Dmael(mu,j,k), Dmael(mu,kappa,k));
    1256    42578884 :             for (i=1; i<=n; i++)
    1257    38458234 :               gmael(B,kappa,i) = subii(gmael(B,kappa,i), gmael(B,j,i));
    1258    68769151 :             for (i=1; i<=d; i++)
    1259    64648515 :               gmael(U,kappa,i) = subii(gmael(U,kappa,i), gmael(U,j,i));
    1260     4120636 :             btop = avma;
    1261     4120636 :             ztmp = subii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
    1262     4120650 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    1263     4120646 :             gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
    1264    35921505 :             for (i=1; i<=j; i++)
    1265    31800856 :               gmael(G,kappa,i) = subii(gmael(G,kappa,i), gmael(G,j,i));
    1266    34521402 :             for (i=j+1; i<kappa; i++)
    1267    30400753 :               gmael(G,kappa,i) = subii(gmael(G,kappa,i), gmael(G,i,j));
    1268    27765192 :             for (i=kappa+1; i<=maxG; i++)
    1269    23644546 :               gmael(G,i,kappa) = subii(gmael(G,i,kappa), gmael(G,i,j));
    1270             :           } else { /* otherwise X = -1 */
    1271    31344944 :             for (k=zeros+1; k<j; k++)
    1272    27239112 :               dpe_addz(Dmael(mu,kappa,k), Dmael(mu,j,k), Dmael(mu,kappa,k));
    1273    42535001 :             for (i=1; i<=n; i++)
    1274    38429169 :               gmael(B,kappa,i) = addii(gmael(B,kappa,i),gmael(B,j,i));
    1275    68415337 :             for (i=1; i<=d; i++)
    1276    64309524 :               gmael(U,kappa,i) = addii(gmael(U,kappa,i),gmael(U,j,i));
    1277     4105813 :             btop = avma;
    1278     4105813 :             ztmp = addii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
    1279     4105829 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    1280     4105831 :             gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
    1281    35727592 :             for (i=1; i<=j; i++)
    1282    31621763 :               gmael(G,kappa,i) = addii(gmael(G,kappa,i), gmael(G,j,i));
    1283    34404424 :             for (i=j+1; i<kappa; i++)
    1284    30298600 :               gmael(G,kappa,i) = addii(gmael(G,kappa,i), gmael(G,i,j));
    1285    27688372 :             for (i=kappa+1; i<=maxG; i++)
    1286    23582549 :               gmael(G,i,kappa) = addii(gmael(G,i,kappa), gmael(G,i,j));
    1287             :           }
    1288     8226469 :           continue;
    1289             :         }
    1290             :         /* we have |X| >= 2 */
    1291     1438064 :         if (tmp->e < BITS_IN_LONG-1)
    1292             :         {
    1293     1263083 :           long xx = (long) pari_rint(ldexp(tmp->d, tmp->e)); /* X fits in an long */
    1294     1263084 :           if (xx > 0) /* = xx */
    1295             :           {
    1296     3090006 :             for (k=zeros+1; k<j; k++)
    1297     2430025 :               dpe_submuluz(Dmael(mu,kappa,k), Dmael(mu,j,k), xx, Dmael(mu,kappa,k));
    1298     4299041 :             for (i=1; i<=n; i++)
    1299     3639060 :               gmael(B,kappa,i) = submuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
    1300    11189367 :             for (i=1; i<=d; i++)
    1301    10529386 :               gmael(U,kappa,i) = submuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
    1302      659981 :             btop = avma;
    1303      659981 :             ztmp = submuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
    1304      659983 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    1305      659983 :             gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
    1306     4158343 :             for (i=1; i<=j; i++)
    1307     3498363 :               gmael(G,kappa,i) = submuliu_inplace(gmael(G,kappa,i), gmael(G,j,i), xx);
    1308     6575620 :             for (i=j+1; i<kappa; i++)
    1309     5915638 :               gmael(G,kappa,i) = submuliu_inplace(gmael(G,kappa,i), gmael(G,i,j), xx);
    1310     3401484 :             for (i=kappa+1; i<=maxG; i++)
    1311     2741501 :               gmael(G,i,kappa) = submuliu_inplace(gmael(G,i,kappa), gmael(G,i,j), xx);
    1312             :           }
    1313             :           else /* = -xx */
    1314             :           {
    1315      603103 :             xx = -xx;
    1316     2904294 :             for (k=zeros+1; k<j; k++)
    1317     2301191 :               dpe_addmuluz(Dmael(mu,kappa,k), Dmael(mu,j,k), xx, Dmael(mu,kappa,k));
    1318     4233219 :             for (i=1; i<=n; i++)
    1319     3630116 :               gmael(B,kappa,i) = addmuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
    1320    10255241 :             for (i=1; i<=d; i++)
    1321     9652142 :               gmael(U,kappa,i) = addmuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
    1322      603099 :             btop = avma;
    1323      603099 :             ztmp = addmuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
    1324      603102 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    1325      603102 :             gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
    1326     3753591 :             for (i=1; i<=j; i++)
    1327     3150491 :               gmael(G,kappa,i) = addmuliu_inplace(gmael(G,kappa,i), gmael(G,j,i), xx);
    1328     6367387 :             for (i=j+1; i<kappa; i++)
    1329     5764289 :               gmael(G,kappa,i) = addmuliu_inplace(gmael(G,kappa,i), gmael(G,i,j), xx);
    1330     2989996 :             for (i=kappa+1; i<=maxG; i++)
    1331     2386899 :               gmael(G,i,kappa) = addmuliu_inplace(gmael(G,i,kappa), gmael(G,i,j), xx);
    1332             :           }
    1333             :         }
    1334             :         else
    1335             :         {
    1336      174981 :           long e = tmp->e - BITS_IN_LONG + 1;
    1337      174981 :           long xx = (long) pari_rint(ldexp(tmp->d, BITS_IN_LONG - 1));
    1338      175006 :           if (xx > 0)
    1339             :           {
    1340      700388 :             for (k=zeros+1; k<j; k++)
    1341             :             {
    1342             :               dpe_t x;
    1343      609904 :               dpe_muluz(Dmael(mu,j,k), xx, &x);
    1344      609904 :               x.e += e;
    1345      609904 :               dpe_subz(Dmael(mu,kappa,k), &x, Dmael(mu,kappa,k));
    1346             :             }
    1347     1732528 :             for (i=1; i<=n; i++)
    1348     1642044 :               gmael(B,kappa,i) = submuliu2n(gmael(B,kappa,i), gmael(B,j,i), xx, e);
    1349      463822 :             for (i=1; i<=d; i++)
    1350      373338 :               gmael(U,kappa,i) = submuliu2n(gmael(U,kappa,i), gmael(U,j,i), xx, e);
    1351       90484 :             btop = avma;
    1352       90484 :             ztmp = submuliu2n(mulshift(gmael(G,j,j), sqrs(xx), 2*e),
    1353       90484 :                 gmael(G,kappa,j), xx, e+1);
    1354       90484 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    1355       90484 :             gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
    1356      791428 :             for (i=1; i<=j; i++)
    1357      700944 :               gmael(G,kappa,i) = submuliu2n(gmael(G,kappa,i), gmael(G,j,i), xx, e);
    1358      751887 :             for (   ; i<kappa; i++)
    1359      661403 :               gmael(G,kappa,i) = submuliu2n(gmael(G,kappa,i), gmael(G,i,j), xx, e);
    1360      151612 :             for (i=kappa+1; i<=maxG; i++)
    1361       61128 :               gmael(G,i,kappa) = submuliu2n(gmael(G,i,kappa), gmael(G,i,j), xx, e);
    1362             :           } else
    1363             :           {
    1364       84522 :             xx = -xx;
    1365      656040 :             for (k=zeros+1; k<j; k++)
    1366             :             {
    1367             :               dpe_t x;
    1368      571518 :               dpe_muluz(Dmael(mu,j,k), xx, &x);
    1369      571518 :               x.e += e;
    1370      571518 :               dpe_addz(Dmael(mu,kappa,k), &x, Dmael(mu,kappa,k));
    1371             :             }
    1372     1704854 :             for (i=1; i<=n; i++)
    1373     1620332 :               gmael(B,kappa,i) = addmuliu2n(gmael(B,kappa,i), gmael(B,j,i), xx, e);
    1374      378247 :             for (i=1; i<=d; i++)
    1375      293725 :               gmael(U,kappa,i) = addmuliu2n(gmael(U,kappa,i), gmael(U,j,i), xx, e);
    1376       84522 :             btop = avma;
    1377       84522 :             ztmp = addmuliu2n(mulshift(gmael(G,j,j), sqrs(xx), 2*e),
    1378       84522 :                 gmael(G,kappa,j), xx, e+1);
    1379       84522 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    1380       84522 :             gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
    1381      741027 :             for (i=1; i<=j; i++)
    1382      656505 :               gmael(G,kappa,i) = addmuliu2n(gmael(G,kappa,i), gmael(G,j,i), xx, e);
    1383      719172 :             for (   ; i<kappa; i++)
    1384      634650 :               gmael(G,kappa,i) = addmuliu2n(gmael(G,kappa,i), gmael(G,i,j), xx, e);
    1385      125875 :             for (i=kappa+1; i<=maxG; i++)
    1386       41353 :               gmael(G,i,kappa) = addmuliu2n(gmael(G,i,kappa), gmael(G,i,j), xx, e);
    1387             :           }
    1388             :         }
    1389             :       }
    1390    12538089 :     if (!go_on) break; /* Anything happened? */
    1391     2165723 :     aa = zeros+1;
    1392             :   }
    1393             : 
    1394    10372366 :   affidpe(gmael(G,kappa,kappa), Del(s,zeros+1));
    1395             :   /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
    1396    61213763 :   for (k=zeros+1; k<=kappa-2; k++)
    1397    50841424 :     dpe_submulz(Del(s,k), Dmael(mu,kappa,k), Dmael(r,kappa,k), Del(s,k+1));
    1398    10372339 :   *pG = G; *pB = B; *pU = U; return 0;
    1399             : }
    1400             : 
    1401             : /* G integral Gram matrix, LLL-reduces (G,B,U) in place [apply base change
    1402             :  * transforms to B and U]. If (keepfirst), never swap with first vector.
    1403             :  * If G = NULL, we compute the Gram matrix incrementally.
    1404             :  * Return -1 on failure, else zeros = dim Kernel (>= 0) */
    1405             : static long
    1406     2409082 : fplll_dpe(GEN *pG, GEN *pB, GEN *pU, GEN *pr, double DELTA, double ETA,
    1407             :       long keepfirst)
    1408             : {
    1409             :   pari_sp av;
    1410     2409082 :   GEN Gtmp, alpha, G = *pG, B = *pB, U = *pU;
    1411     2409082 :   long d, maxG, kappa, kappa2, i, j, zeros, kappamax, incgram = !G, cnt = 0;
    1412             :   dpe_t delta, eta, **mu, **r, *s;
    1413     2409082 :   affdbldpe(DELTA,&delta);
    1414     2409078 :   affdbldpe(ETA,&eta);
    1415             : 
    1416     2409080 :   if (incgram)
    1417             :   { /* incremental Gram matrix */
    1418     2327020 :     maxG = 2; d = lg(B)-1;
    1419     2327020 :     G = zeromatcopy(d, d);
    1420             :   }
    1421             :   else
    1422       82060 :     maxG = d = lg(G)-1;
    1423             : 
    1424     2409084 :   mu = cget_dpemat(d+1);
    1425     2409078 :   r  = cget_dpemat(d+1);
    1426     2409083 :   s  = cget_dpevec(d+1);
    1427    10823045 :   for (j = 1; j <= d; j++)
    1428             :   {
    1429     8413960 :     mu[j]= cget_dpevec(d+1);
    1430     8413961 :     r[j] = cget_dpevec(d+1);
    1431             :   }
    1432     2409085 :   Gtmp = cgetg(d+1, t_VEC);
    1433     2409079 :   alpha = cgetg(d+1, t_VECSMALL);
    1434     2409068 :   av = avma;
    1435             : 
    1436             :   /* Step2: Initializing the main loop */
    1437     2409068 :   kappamax = 1;
    1438     2409068 :   i = 1;
    1439             :   do {
    1440     2762191 :     if (incgram) gmael(G,i,i) = ZV_dotsquare(gel(B,i));
    1441     2762097 :     affidpe(gmael(G,i,i), Dmael(r,i,i));
    1442     2762112 :   } while (!signe(gmael(G,i,i)) && ++i <= d);
    1443     2408989 :   zeros = i-1; /* all basis vectors b_i with i <= zeros are zero vectors */
    1444     2408989 :   kappa = i;
    1445    10469611 :   for (i=zeros+1; i<=d; i++) alpha[i]=1;
    1446             : 
    1447    12781383 :   while (++kappa <= d)
    1448             :   {
    1449    10372364 :     if (kappa > kappamax)
    1450             :     {
    1451     5651634 :       if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
    1452     5651634 :       kappamax = kappa;
    1453     5651634 :       if (incgram)
    1454             :       {
    1455    24749675 :         for (i=zeros+1; i<=kappa; i++)
    1456    19328856 :           gmael(G,kappa,i) = ZV_dotproduct(gel(B,kappa), gel(B,i));
    1457     5420819 :         maxG = kappamax;
    1458             :       }
    1459             :     }
    1460             :     /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
    1461    10372014 :     if (Babai_dpe(av, kappa, &G,&B,&U, mu,r,s, alpha[kappa], zeros, maxG, &eta))
    1462           2 :     { *pG = incgram? NULL: G; *pB = B; *pU = U; return -1; }
    1463    20646754 :     if ((keepfirst && kappa == 2) ||
    1464    10274377 :         dpe_cmpmul(Dmael(r,kappa-1,kappa-1), &delta, Del(s,kappa-1)) <= 0)
    1465             :     { /* Step4: Success of Lovasz's condition */
    1466     8855412 :       alpha[kappa] = kappa;
    1467     8855412 :       dpe_submulz(Del(s,kappa-1), Dmael(mu,kappa,kappa-1), Dmael(r,kappa,kappa-1), Dmael(r,kappa,kappa));
    1468     8855411 :       continue;
    1469             :     }
    1470             :     /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
    1471     1516965 :     if (DEBUGLEVEL>=4 && kappa==kappamax && Del(s,kappa-1)->d)
    1472           0 :       if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", Del(s,1)->e-1); }
    1473     1516965 :     kappa2 = kappa;
    1474             :     do {
    1475     4755902 :       kappa--;
    1476     4755902 :       if (kappa < zeros+2 + (keepfirst ? 1: 0)) break;
    1477     4498324 :     } while (dpe_cmpmul(Dmael(r,kappa-1,kappa-1), &delta, Del(s,kappa-1)) >= 0);
    1478     1516965 :     update_alpha(alpha, kappa, kappa2, kappamax);
    1479             : 
    1480             :     /* Step6: Update the mu's and r's */
    1481     1516965 :     dperotate(mu, kappa2, kappa);
    1482     1516965 :     dperotate(r, kappa2, kappa);
    1483     1516965 :     affdpe(Del(s,kappa), Dmael(r,kappa,kappa));
    1484             : 
    1485             :     /* Step7: Update G, B, U */
    1486     1516965 :     if (U) rotate(U, kappa2, kappa);
    1487     1516965 :     if (B) rotate(B, kappa2, kappa);
    1488     1516965 :     rotateG(G,kappa2,kappa, maxG, Gtmp);
    1489             : 
    1490             :     /* Step8: Prepare the next loop iteration */
    1491     1516965 :     if (kappa == zeros+1 && !signe(gmael(G,kappa,kappa)))
    1492             :     {
    1493       35173 :       zeros++; kappa++;
    1494       35173 :       affidpe(gmael(G,kappa,kappa), Dmael(r,kappa,kappa));
    1495             :     }
    1496             :   }
    1497     2409019 :   if (pr) *pr = dpeM_diagonal_shallow(r,d);
    1498     2409019 :   *pG = G; *pB = B; *pU = U; return zeros; /* success */
    1499             : }
    1500             : 
    1501             : 
    1502             : /************************** PROVED version (t_INT) *************************/
    1503             : 
    1504             : /* Babai's Nearest Plane algorithm (iterative).
    1505             :  * Size-reduces b_kappa using mu_{i,j} and r_{i,j} for j<=i <kappa
    1506             :  * Update B[,kappa]; compute mu_{kappa,j}, r_{kappa,j} for j<=kappa and s[kappa]
    1507             :  * mu, r, s updated in place (affrr). Return 1 on failure, else 0. */
    1508             : static int
    1509        2208 : Babai(pari_sp av, long kappa, GEN *pG, GEN *pB, GEN *pU, GEN mu, GEN r, GEN s,
    1510             :       long a, long zeros, long maxG, GEN eta, long prec)
    1511             : {
    1512        2208 :   GEN G = *pG, B = *pB, U = *pU, ztmp;
    1513        2208 :   long k, aa = a > zeros? a: zeros+1;
    1514        2208 :   const long n = B? nbrows(B): 0, d = U ? lg(U)-1: 0, bit = prec2nbits(prec);
    1515        2208 :   long emaxmu = EX0, emax2mu = EX0;
    1516             :   /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
    1517             : 
    1518        1262 :   for (;;) {
    1519        3470 :     int go_on = 0;
    1520        3470 :     long i, j, emax3mu = emax2mu;
    1521             : 
    1522        3470 :     if (gc_needed(av,2))
    1523             :     {
    1524           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
    1525           0 :       gc_lll(av,3,&G,&B,&U);
    1526             :     }
    1527             :     /* Step2: compute the GSO for stage kappa */
    1528        3470 :     emax2mu = emaxmu; emaxmu = EX0;
    1529       24881 :     for (j=aa; j<kappa; j++)
    1530             :     {
    1531       21411 :       pari_sp btop = avma;
    1532       21411 :       GEN g = gmael(G,kappa,j);
    1533      209529 :       for (k = zeros+1; k < j; k++)
    1534      188118 :         g = mpsub(g, mulrr(gmael(mu,j,k), gmael(r,kappa,k)));
    1535       21411 :       mpaff(g, gmael(r,kappa,j));
    1536       21411 :       affrr(divrr(gmael(r,kappa,j), gmael(r,j,j)), gmael(mu,kappa,j));
    1537       21411 :       emaxmu = maxss(emaxmu, expo(gmael(mu,kappa,j)));
    1538       21411 :       set_avma(btop);
    1539             :     }
    1540        3470 :     if (emax3mu != EX0 && emax3mu <= emax2mu + 5) /* precision too low */
    1541           0 :     { *pG = G; *pB = B; *pU = U; return 1; }
    1542             : 
    1543       25500 :     for (j=kappa-1; j>zeros; j--)
    1544       23292 :       if (abscmprr(gmael(mu,kappa,j), eta) > 0) { go_on = 1; break; }
    1545             : 
    1546             :     /* Step3--5: compute the X_j's  */
    1547        3470 :     if (go_on)
    1548       15993 :       for (j=kappa-1; j>zeros; j--)
    1549             :       {
    1550             :         pari_sp btop;
    1551       14731 :         GEN tmp = gmael(mu,kappa,j);
    1552       14731 :         if (absrsmall(tmp)) continue; /* size-reduced */
    1553             : 
    1554        8741 :         if (gc_needed(av,2))
    1555             :         {
    1556           0 :           if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
    1557           0 :           gc_lll(av,3,&G,&B,&U);
    1558             :         }
    1559        8741 :         btop = avma;
    1560             :         /* we consider separately the case |X| = 1 */
    1561        8741 :         if (absrsmall2(tmp))
    1562             :         {
    1563        3721 :           if (signe(tmp) > 0) { /* in this case, X = 1 */
    1564       17772 :             for (k=zeros+1; k<j; k++)
    1565       15910 :               affrr(subrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
    1566        1862 :             set_avma(btop);
    1567       45875 :             for (i=1; i<=n; i++)
    1568       44013 :               gmael(B,kappa,i) = subii(gmael(B,kappa,i), gmael(B,j,i));
    1569        7955 :             for (i=1; i<=d; i++)
    1570        6093 :               gmael(U,kappa,i) = subii(gmael(U,kappa,i), gmael(U,j,i));
    1571        1862 :             btop = avma;
    1572        1862 :             ztmp = subii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
    1573        1862 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    1574        1862 :             gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
    1575       19634 :             for (i=1; i<=j; i++)
    1576       17772 :               gmael(G,kappa,i) = subii(gmael(G,kappa,i), gmael(G,j,i));
    1577       17730 :             for (i=j+1; i<kappa; i++)
    1578       15868 :               gmael(G,kappa,i) = subii(gmael(G,kappa,i), gmael(G,i,j));
    1579        7280 :             for (i=kappa+1; i<=maxG; i++)
    1580        5418 :               gmael(G,i,kappa) = subii(gmael(G,i,kappa), gmael(G,i,j));
    1581             :           } else { /* otherwise X = -1 */
    1582       18073 :             for (k=zeros+1; k<j; k++)
    1583       16214 :               affrr(addrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
    1584        1859 :             set_avma(btop);
    1585       45822 :             for (i=1; i<=n; i++)
    1586       43963 :               gmael(B,kappa,i) = addii(gmael(B,kappa,i),gmael(B,j,i));
    1587        7934 :             for (i=1; i<=d; i++)
    1588        6075 :               gmael(U,kappa,i) = addii(gmael(U,kappa,i),gmael(U,j,i));
    1589        1859 :             btop = avma;
    1590        1859 :             ztmp = addii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
    1591        1859 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    1592        1859 :             gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
    1593       19932 :             for (i=1; i<=j; i++)
    1594       18073 :               gmael(G,kappa,i) = addii(gmael(G,kappa,i), gmael(G,j,i));
    1595       17334 :             for (i=j+1; i<kappa; i++)
    1596       15475 :               gmael(G,kappa,i) = addii(gmael(G,kappa,i), gmael(G,i,j));
    1597        7296 :             for (i=kappa+1; i<=maxG; i++)
    1598        5437 :               gmael(G,i,kappa) = addii(gmael(G,i,kappa), gmael(G,i,j));
    1599             :           }
    1600        3721 :           continue;
    1601             :         }
    1602             :         /* we have |X| >= 2 */
    1603        5020 :         if (expo(tmp) < BITS_IN_LONG)
    1604             :         {
    1605        4500 :           ulong xx = roundr_safe(tmp)[2]; /* X fits in an ulong */
    1606        4500 :           if (signe(tmp) > 0) /* = xx */
    1607             :           {
    1608       25486 :             for (k=zeros+1; k<j; k++)
    1609       23274 :               affrr(subrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
    1610       23274 :                   gmael(mu,kappa,k));
    1611        2212 :             set_avma(btop);
    1612       68948 :             for (i=1; i<=n; i++)
    1613       66736 :               gmael(B,kappa,i) = submuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
    1614        3796 :             for (i=1; i<=d; i++)
    1615        1584 :               gmael(U,kappa,i) = submuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
    1616        2212 :             btop = avma;
    1617        2212 :             ztmp = submuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
    1618        2212 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    1619        2212 :             gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
    1620       27698 :             for (i=1; i<=j; i++)
    1621       25486 :               gmael(G,kappa,i) = submuliu_inplace(gmael(G,kappa,i), gmael(G,j,i), xx);
    1622       30761 :             for (i=j+1; i<kappa; i++)
    1623       28549 :               gmael(G,kappa,i) = submuliu_inplace(gmael(G,kappa,i), gmael(G,i,j), xx);
    1624        6073 :             for (i=kappa+1; i<=maxG; i++)
    1625        3861 :               gmael(G,i,kappa) = submuliu_inplace(gmael(G,i,kappa), gmael(G,i,j), xx);
    1626             :           }
    1627             :           else /* = -xx */
    1628             :           {
    1629       26514 :             for (k=zeros+1; k<j; k++)
    1630       24226 :               affrr(addrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
    1631       24226 :                   gmael(mu,kappa,k));
    1632        2288 :             set_avma(btop);
    1633       71318 :             for (i=1; i<=n; i++)
    1634       69030 :               gmael(B,kappa,i) = addmuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
    1635        3926 :             for (i=1; i<=d; i++)
    1636        1638 :               gmael(U,kappa,i) = addmuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
    1637        2288 :             btop = avma;
    1638        2288 :             ztmp = addmuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
    1639        2288 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    1640        2288 :             gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
    1641       28802 :             for (i=1; i<=j; i++)
    1642       26514 :               gmael(G,kappa,i) = addmuliu_inplace(gmael(G,kappa,i), gmael(G,j,i), xx);
    1643       31397 :             for (i=j+1; i<kappa; i++)
    1644       29109 :               gmael(G,kappa,i) = addmuliu_inplace(gmael(G,kappa,i), gmael(G,i,j), xx);
    1645        6391 :             for (i=kappa+1; i<=maxG; i++)
    1646        4103 :               gmael(G,i,kappa) = addmuliu_inplace(gmael(G,i,kappa), gmael(G,i,j), xx);
    1647             :           }
    1648             :         }
    1649             :         else
    1650             :         {
    1651             :           long e;
    1652         520 :           GEN X = truncexpo(tmp, bit, &e); /* tmp ~ X * 2^e */
    1653         520 :           btop = avma;
    1654        6849 :           for (k=zeros+1; k<j; k++)
    1655             :           {
    1656        6329 :             GEN x = mulir(X, gmael(mu,j,k));
    1657        6329 :             if (e) shiftr_inplace(x, e);
    1658        6329 :             affrr(subrr(gmael(mu,kappa,k), x), gmael(mu,kappa,k));
    1659             :           }
    1660         520 :           set_avma(btop);
    1661       16677 :           for (i=1; i<=n; i++)
    1662       16157 :             gmael(B,kappa,i) = submulshift(gmael(B,kappa,i), gmael(B,j,i), X, e);
    1663         709 :           for (i=1; i<=d; i++)
    1664         189 :             gmael(U,kappa,i) = submulshift(gmael(U,kappa,i), gmael(U,j,i), X, e);
    1665         520 :           btop = avma;
    1666         520 :           ztmp = submulshift(mulshift(gmael(G,j,j), sqri(X), 2*e),
    1667         520 :               gmael(G,kappa,j), X, e+1);
    1668         520 :           ztmp = addii(gmael(G,kappa,kappa), ztmp);
    1669         520 :           gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
    1670        7369 :           for (i=1; i<=j; i++)
    1671        6849 :             gmael(G,kappa,i) = submulshift(gmael(G,kappa,i), gmael(G,j,i), X, e);
    1672        7435 :           for (   ; i<kappa; i++)
    1673        6915 :             gmael(G,kappa,i) = submulshift(gmael(G,kappa,i), gmael(G,i,j), X, e);
    1674         580 :           for (i=kappa+1; i<=maxG; i++)
    1675          60 :             gmael(G,i,kappa) = submulshift(gmael(G,i,kappa), gmael(G,i,j), X, e);
    1676             :         }
    1677             :       }
    1678        3470 :     if (!go_on) break; /* Anything happened? */
    1679        1262 :     aa = zeros+1;
    1680             :   }
    1681             : 
    1682        2208 :   affir(gmael(G,kappa,kappa), gel(s,zeros+1));
    1683             :   /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
    1684        2208 :   av = avma;
    1685       20786 :   for (k=zeros+1; k<=kappa-2; k++)
    1686       18578 :     affrr(subrr(gel(s,k), mulrr(gmael(mu,kappa,k), gmael(r,kappa,k))),
    1687       18578 :           gel(s,k+1));
    1688        2208 :   *pG = G; *pB = B; *pU = U; return gc_bool(av, 0);
    1689             : }
    1690             : 
    1691             : /* G integral Gram matrix, LLL-reduces (G,B,U) in place [apply base change
    1692             :  * transforms to B and U]. If (keepfirst), never swap with first vector.
    1693             :  * If G = NULL, we compute the Gram matrix incrementally.
    1694             :  * Return -1 on failure, else zeros = dim Kernel (>= 0) */
    1695             : static long
    1696           2 : fplll(GEN *pG, GEN *pB, GEN *pU, GEN *pr, double DELTA, double ETA,
    1697             :       long keepfirst, long prec)
    1698             : {
    1699             :   pari_sp av, av2;
    1700           2 :   GEN mu, r, s, tmp, Gtmp, alpha, G = *pG, B = *pB, U = *pU;
    1701           2 :   GEN delta = dbltor(DELTA), eta = dbltor(ETA);
    1702           2 :   long d, maxG, kappa, kappa2, i, j, zeros, kappamax, incgram = !G, cnt = 0;
    1703             : 
    1704           2 :   if (incgram)
    1705             :   { /* incremental Gram matrix */
    1706           2 :     maxG = 2; d = lg(B)-1;
    1707           2 :     G = zeromatcopy(d, d);
    1708             :   }
    1709             :   else
    1710           0 :     maxG = d = lg(G)-1;
    1711             : 
    1712           2 :   mu = cgetg(d+1, t_MAT);
    1713           2 :   r  = cgetg(d+1, t_MAT);
    1714           2 :   s  = cgetg(d+1, t_VEC);
    1715          43 :   for (j = 1; j <= d; j++)
    1716             :   {
    1717          41 :     GEN M = cgetg(d+1, t_COL), R = cgetg(d+1, t_COL);
    1718          41 :     gel(mu,j)= M;
    1719          41 :     gel(r,j) = R;
    1720          41 :     gel(s,j) = cgetr(prec);
    1721        1146 :     for (i = 1; i <= d; i++)
    1722             :     {
    1723        1105 :       gel(R,i) = cgetr(prec);
    1724        1105 :       gel(M,i) = cgetr(prec);
    1725             :     }
    1726             :   }
    1727           2 :   Gtmp = cgetg(d+1, t_VEC);
    1728           2 :   alpha = cgetg(d+1, t_VECSMALL);
    1729           2 :   av = avma;
    1730             : 
    1731             :   /* Step2: Initializing the main loop */
    1732           2 :   kappamax = 1;
    1733           2 :   i = 1;
    1734             :   do {
    1735           2 :     if (incgram) gmael(G,i,i) = ZV_dotsquare(gel(B,i));
    1736           2 :     affir(gmael(G,i,i), gmael(r,i,i));
    1737           2 :   } while (!signe(gmael(G,i,i)) && ++i <= d);
    1738           2 :   zeros = i-1; /* all basis vectors b_i with i <= zeros are zero vectors */
    1739           2 :   kappa = i;
    1740          43 :   for (i=zeros+1; i<=d; i++) alpha[i]=1;
    1741             : 
    1742        2210 :   while (++kappa <= d)
    1743             :   {
    1744        2208 :     if (kappa > kappamax)
    1745             :     {
    1746          39 :       if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
    1747          39 :       kappamax = kappa;
    1748          39 :       if (incgram)
    1749             :       {
    1750         610 :         for (i=zeros+1; i<=kappa; i++)
    1751         571 :           gmael(G,kappa,i) = ZV_dotproduct(gel(B,kappa), gel(B,i));
    1752          39 :         maxG = kappamax;
    1753             :       }
    1754             :     }
    1755             :     /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
    1756        2208 :     if (Babai(av, kappa, &G,&B,&U, mu,r,s, alpha[kappa], zeros, maxG, eta, prec))
    1757           0 :     { *pG = incgram? NULL: G; *pB = B; *pU = U; return -1; }
    1758        2208 :     av2 = avma;
    1759        4416 :     if ((keepfirst && kappa == 2) ||
    1760        2208 :         cmprr(mulrr(gmael(r,kappa-1,kappa-1), delta), gel(s,kappa-1)) <= 0)
    1761             :     { /* Step4: Success of Lovasz's condition */
    1762        1464 :       alpha[kappa] = kappa;
    1763        1464 :       tmp = mulrr(gmael(mu,kappa,kappa-1), gmael(r,kappa,kappa-1));
    1764        1464 :       affrr(subrr(gel(s,kappa-1), tmp), gmael(r,kappa,kappa));
    1765        1464 :       set_avma(av2); continue;
    1766             :     }
    1767             :     /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
    1768         744 :     if (DEBUGLEVEL>=4 && kappa==kappamax && signe(gel(s,kappa-1)))
    1769           0 :       if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", expo(gel(s,1))); }
    1770         744 :     kappa2 = kappa;
    1771             :     do {
    1772        2169 :       kappa--;
    1773        2169 :       if (kappa < zeros+2 + (keepfirst ? 1: 0)) break;
    1774        2033 :       tmp = mulrr(gmael(r,kappa-1,kappa-1), delta);
    1775        2033 :     } while (cmprr(gel(s,kappa-1), tmp) <= 0);
    1776         744 :     set_avma(av2);
    1777         744 :     update_alpha(alpha, kappa, kappa2, kappamax);
    1778             : 
    1779             :     /* Step6: Update the mu's and r's */
    1780         744 :     rotate(mu, kappa2, kappa);
    1781         744 :     rotate(r, kappa2, kappa);
    1782         744 :     affrr(gel(s,kappa), gmael(r,kappa,kappa));
    1783             : 
    1784             :     /* Step7: Update G, B, U */
    1785         744 :     if (U) rotate(U, kappa2, kappa);
    1786         744 :     if (B) rotate(B, kappa2, kappa);
    1787         744 :     rotateG(G,kappa2,kappa, maxG, Gtmp);
    1788             : 
    1789             :     /* Step8: Prepare the next loop iteration */
    1790         744 :     if (kappa == zeros+1 && !signe(gmael(G,kappa,kappa)))
    1791             :     {
    1792           0 :       zeros++; kappa++;
    1793           0 :       affir(gmael(G,kappa,kappa), gmael(r,kappa,kappa));
    1794             :     }
    1795             :   }
    1796           2 :   if (pr) *pr = RgM_diagonal_shallow(r);
    1797           2 :   *pG = G; *pB = B; *pU = U; return zeros; /* success */
    1798             : }
    1799             : 
    1800             : /* Assume x a ZM, if pN != NULL, set it to Gram-Schmidt (squared) norms
    1801             :  * The following modes are supported:
    1802             :  * - flag & LLL_INPLACE: x a lattice basis, return x*U
    1803             :  * - flag & LLL_GRAM: x a Gram matrix / else x a lattice basis; return
    1804             :  *     LLL base change matrix U [LLL_IM]
    1805             :  *     kernel basis [LLL_KER, nonreduced]
    1806             :  *     both [LLL_ALL] */
    1807             : GEN
    1808     2531887 : ZM_lll_norms(GEN x, double DELTA, long flag, GEN *pN)
    1809             : {
    1810     2531887 :   pari_sp av = avma;
    1811     2531887 :   const double ETA = 0.51;
    1812     2531887 :   const long keepfirst = flag & LLL_KEEP_FIRST;
    1813     2531887 :   long p, zeros = -1, n = lg(x)-1;
    1814             :   GEN G, B, U;
    1815             :   pari_timer T;
    1816     2531887 :   if (n <= 1) return lll_trivial(x, flag);
    1817     2422859 :   if (nbrows(x)==0)
    1818             :   {
    1819       13861 :     if (flag & LLL_KER) return matid(n);
    1820       13861 :     if (flag & (LLL_INPLACE|LLL_IM)) return cgetg(1,t_MAT);
    1821           0 :     retmkvec2(matid(n), cgetg(1,t_MAT));
    1822             :   }
    1823     2409045 :   x = RgM_shallowcopy(x);
    1824     2409056 :   if (flag & LLL_GRAM)
    1825       82040 :   { G = x; B = NULL; U = matid(n); }
    1826             :   else
    1827     2327016 :   { G = NULL; B = x; U = (flag & LLL_INPLACE)? NULL: matid(n); }
    1828     2409070 :   if(DEBUGLEVEL>=4) timer_start(&T);
    1829     2409070 :   if (!(flag & LLL_GRAM))
    1830             :   {
    1831             :     long t;
    1832     2327011 :     if(DEBUGLEVEL>=4)
    1833           0 :       err_printf("Entering L^2 (double): LLL-parameters (%.3f,%.3f)\n",
    1834             :                  DELTA,ETA);
    1835     2327011 :     zeros = fplll_fast(&B, &U, DELTA, ETA, keepfirst);
    1836     2327027 :     if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
    1837     2330472 :     for (p = DEFAULTPREC, t = 0; zeros < 0 && t < 1 ; p += EXTRAPREC64, t++)
    1838             :     {
    1839        3447 :       if (DEBUGLEVEL>=4)
    1840           0 :         err_printf("Entering L^2 (heuristic): LLL-parameters (%.3f,%.3f), prec = %d/%d\n", DELTA, ETA, p, p);
    1841        3447 :       zeros = fplll_heuristic(&B, &U, DELTA, ETA, keepfirst, p, p);
    1842        3447 :       gc_lll(av, 2, &B, &U);
    1843        3447 :       if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
    1844             :     }
    1845             :   }
    1846     2409084 :   if(DEBUGLEVEL>=4)
    1847           0 :     err_printf("Entering L^2 (dpe): LLL-parameters (%.3f,%.3f)\n", DELTA,ETA);
    1848     2409084 :   zeros = fplll_dpe(&G, &B, &U, pN, DELTA, ETA, keepfirst);
    1849     2409022 :   if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
    1850     2409048 :   if (zeros < 0)
    1851           2 :     for (p = DEFAULTPREC;; p += EXTRAPREC64)
    1852             :     {
    1853           2 :       if (DEBUGLEVEL>=4)
    1854           0 :         err_printf("Entering L^2: LLL-parameters (%.3f,%.3f), prec = %d\n",
    1855             :                     DELTA,ETA, p);
    1856           2 :       zeros = fplll(&G, &B, &U, pN, DELTA, ETA, keepfirst, p);
    1857           2 :       if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
    1858           2 :       if (zeros >= 0) break;
    1859           0 :       gc_lll(av, 3, &G, &B, &U);
    1860             :     }
    1861     2409048 :   return lll_finish(U? U: B, zeros, flag);
    1862             : }
    1863             : 
    1864             : /********************************************************************/
    1865             : /**                                                                **/
    1866             : /**                        LLL OVER K[X]                           **/
    1867             : /**                                                                **/
    1868             : /********************************************************************/
    1869             : static int
    1870         378 : pslg(GEN x)
    1871             : {
    1872             :   long tx;
    1873         378 :   if (gequal0(x)) return 2;
    1874         336 :   tx = typ(x); return is_scalar_t(tx)? 3: lg(x);
    1875             : }
    1876             : 
    1877             : static int
    1878         147 : REDgen(long k, long l, GEN h, GEN L, GEN B)
    1879             : {
    1880         147 :   GEN q, u = gcoeff(L,k,l);
    1881             :   long i;
    1882             : 
    1883         147 :   if (pslg(u) < pslg(B)) return 0;
    1884             : 
    1885         105 :   q = gneg(gdeuc(u,B));
    1886         105 :   gel(h,k) = gadd(gel(h,k), gmul(q,gel(h,l)));
    1887         105 :   for (i=1; i<l; i++) gcoeff(L,k,i) = gadd(gcoeff(L,k,i), gmul(q,gcoeff(L,l,i)));
    1888         105 :   gcoeff(L,k,l) = gadd(gcoeff(L,k,l), gmul(q,B)); return 1;
    1889             : }
    1890             : 
    1891             : static int
    1892         147 : do_SWAPgen(GEN h, GEN L, GEN B, long k, GEN fl, int *flc)
    1893             : {
    1894             :   GEN p1, la, la2, Bk;
    1895             :   long ps1, ps2, i, j, lx;
    1896             : 
    1897         147 :   if (!fl[k-1]) return 0;
    1898             : 
    1899         105 :   la = gcoeff(L,k,k-1); la2 = gsqr(la);
    1900         105 :   Bk = gel(B,k);
    1901         105 :   if (fl[k])
    1902             :   {
    1903          42 :     GEN q = gadd(la2, gmul(gel(B,k-1),gel(B,k+1)));
    1904          42 :     ps1 = pslg(gsqr(Bk));
    1905          42 :     ps2 = pslg(q);
    1906          42 :     if (ps1 <= ps2 && (ps1 < ps2 || !*flc)) return 0;
    1907          21 :     *flc = (ps1 != ps2);
    1908          21 :     gel(B,k) = gdiv(q, Bk);
    1909             :   }
    1910             : 
    1911          84 :   swap(gel(h,k-1), gel(h,k)); lx = lg(L);
    1912          84 :   for (j=1; j<k-1; j++) swap(gcoeff(L,k-1,j), gcoeff(L,k,j));
    1913          84 :   if (fl[k])
    1914             :   {
    1915          21 :     for (i=k+1; i<lx; i++)
    1916             :     {
    1917           0 :       GEN t = gcoeff(L,i,k);
    1918           0 :       p1 = gsub(gmul(gel(B,k+1),gcoeff(L,i,k-1)), gmul(la,t));
    1919           0 :       gcoeff(L,i,k) = gdiv(p1, Bk);
    1920           0 :       p1 = gadd(gmul(la,gcoeff(L,i,k-1)), gmul(gel(B,k-1),t));
    1921           0 :       gcoeff(L,i,k-1) = gdiv(p1, Bk);
    1922             :     }
    1923             :   }
    1924          63 :   else if (!gequal0(la))
    1925             :   {
    1926          21 :     p1 = gdiv(la2, Bk);
    1927          21 :     gel(B,k+1) = gel(B,k) = p1;
    1928          21 :     for (i=k+2; i<=lx; i++) gel(B,i) = gdiv(gmul(p1,gel(B,i)),Bk);
    1929          21 :     for (i=k+1; i<lx; i++)
    1930           0 :       gcoeff(L,i,k-1) = gdiv(gmul(la,gcoeff(L,i,k-1)), Bk);
    1931          21 :     for (j=k+1; j<lx-1; j++)
    1932           0 :       for (i=j+1; i<lx; i++)
    1933           0 :         gcoeff(L,i,j) = gdiv(gmul(p1,gcoeff(L,i,j)), Bk);
    1934             :   }
    1935             :   else
    1936             :   {
    1937          42 :     gcoeff(L,k,k-1) = gen_0;
    1938          42 :     for (i=k+1; i<lx; i++)
    1939             :     {
    1940           0 :       gcoeff(L,i,k) = gcoeff(L,i,k-1);
    1941           0 :       gcoeff(L,i,k-1) = gen_0;
    1942             :     }
    1943          42 :     gel(B,k) = gel(B,k-1); fl[k] = 1; fl[k-1] = 0;
    1944             :   }
    1945          84 :   return 1;
    1946             : }
    1947             : 
    1948             : static void
    1949         126 : incrementalGSgen(GEN x, GEN L, GEN B, long k, GEN fl)
    1950             : {
    1951         126 :   GEN u = NULL; /* gcc -Wall */
    1952             :   long i, j;
    1953         315 :   for (j = 1; j <= k; j++)
    1954         189 :     if (j==k || fl[j])
    1955             :     {
    1956         189 :       u = gcoeff(x,k,j);
    1957         189 :       if (!is_extscalar_t(typ(u))) pari_err_TYPE("incrementalGSgen",u);
    1958         252 :       for (i=1; i<j; i++)
    1959          63 :         if (fl[i])
    1960             :         {
    1961          63 :           u = gsub(gmul(gel(B,i+1),u), gmul(gcoeff(L,k,i),gcoeff(L,j,i)));
    1962          63 :           u = gdiv(u, gel(B,i));
    1963             :         }
    1964         189 :       gcoeff(L,k,j) = u;
    1965             :     }
    1966         126 :   if (gequal0(u)) gel(B,k+1) = gel(B,k);
    1967             :   else
    1968             :   {
    1969          84 :     gel(B,k+1) = gcoeff(L,k,k); gcoeff(L,k,k) = gen_1; fl[k] = 1;
    1970             :   }
    1971         126 : }
    1972             : 
    1973             : static GEN
    1974         126 : lllgramallgen(GEN x, long flag)
    1975             : {
    1976         126 :   long lx = lg(x), i, j, k, l, n;
    1977             :   pari_sp av;
    1978             :   GEN B, L, h, fl;
    1979             :   int flc;
    1980             : 
    1981         126 :   n = lx-1; if (n<=1) return lll_trivial(x,flag);
    1982          63 :   if (lgcols(x) != lx) pari_err_DIM("lllgramallgen");
    1983             : 
    1984          63 :   fl = cgetg(lx, t_VECSMALL);
    1985             : 
    1986          63 :   av = avma;
    1987          63 :   B = scalarcol_shallow(gen_1, lx);
    1988          63 :   L = cgetg(lx,t_MAT);
    1989         189 :   for (j=1; j<lx; j++) { gel(L,j) = zerocol(n); fl[j] = 0; }
    1990             : 
    1991          63 :   h = matid(n);
    1992         189 :   for (i=1; i<lx; i++)
    1993         126 :     incrementalGSgen(x, L, B, i, fl);
    1994          63 :   flc = 0;
    1995          63 :   for(k=2;;)
    1996             :   {
    1997         147 :     if (REDgen(k, k-1, h, L, gel(B,k))) flc = 1;
    1998         147 :     if (do_SWAPgen(h, L, B, k, fl, &flc)) { if (k > 2) k--; }
    1999             :     else
    2000             :     {
    2001          63 :       for (l=k-2; l>=1; l--)
    2002           0 :         if (REDgen(k, l, h, L, gel(B,l+1))) flc = 1;
    2003          63 :       if (++k > n) break;
    2004             :     }
    2005          84 :     if (gc_needed(av,1))
    2006             :     {
    2007           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"lllgramallgen");
    2008           0 :       gerepileall(av,3,&B,&L,&h);
    2009             :     }
    2010             :   }
    2011         105 :   k=1; while (k<lx && !fl[k]) k++;
    2012          63 :   return lll_finish(h,k-1,flag);
    2013             : }
    2014             : 
    2015             : static int
    2016       57163 : RgM_square(GEN x) { long l = lg(x); return l == 1 || l == lgcols(x); }
    2017             : static GEN
    2018         126 : lllallgen(GEN x, long flag)
    2019             : {
    2020         126 :   pari_sp av = avma;
    2021         126 :   if (!(flag & LLL_GRAM)) x = gram_matrix(x);
    2022          42 :   else if (!RgM_square(x)) pari_err_DIM("qflllgram");
    2023         126 :   return gerepilecopy(av, lllgramallgen(x, flag));
    2024             : }
    2025             : GEN
    2026          42 : lllgen(GEN x) { return lllallgen(x, LLL_IM); }
    2027             : GEN
    2028          42 : lllkerimgen(GEN x) { return lllallgen(x, LLL_ALL); }
    2029             : GEN
    2030           0 : lllgramgen(GEN x)  { return lllallgen(x, LLL_IM|LLL_GRAM); }
    2031             : GEN
    2032          42 : lllgramkerimgen(GEN x)  { return lllallgen(x, LLL_ALL|LLL_GRAM); }
    2033             : 
    2034             : static GEN
    2035       55922 : lllall(GEN x, long flag)
    2036       55922 : { pari_sp av = avma; return gerepilecopy(av, ZM_lll(x, LLLDFT, flag)); }
    2037             : GEN
    2038       14662 : lllint(GEN x) { return lllall(x, LLL_IM); }
    2039             : GEN
    2040          35 : lllkerim(GEN x) { return lllall(x, LLL_ALL); }
    2041             : GEN
    2042       41183 : lllgramint(GEN x)
    2043       41183 : { if (!RgM_square(x)) pari_err_DIM("qflllgram");
    2044       41183 :   return lllall(x, LLL_IM | LLL_GRAM); }
    2045             : GEN
    2046          35 : lllgramkerim(GEN x)
    2047          35 : { if (!RgM_square(x)) pari_err_DIM("qflllgram");
    2048          35 :   return lllall(x, LLL_ALL | LLL_GRAM); }
    2049             : 
    2050             : GEN
    2051      645602 : lllfp(GEN x, double D, long flag)
    2052             : {
    2053      645602 :   long n = lg(x)-1;
    2054      645602 :   pari_sp av = avma;
    2055             :   GEN h;
    2056      645602 :   if (n <= 1) return lll_trivial(x,flag);
    2057      561811 :   if ((flag & LLL_GRAM) && !RgM_square(x)) pari_err_DIM("qflllgram");
    2058      561797 :   h = ZM_lll(RgM_rescale_to_int(x), D, flag);
    2059      561768 :   return gerepilecopy(av, h);
    2060             : }
    2061             : 
    2062             : GEN
    2063       15673 : lllgram(GEN x) { return lllfp(x,LLLDFT,LLL_GRAM|LLL_IM); }
    2064             : GEN
    2065      556349 : lll(GEN x) { return lllfp(x,LLLDFT,LLL_IM); }
    2066             : 
    2067             : GEN
    2068         301 : qflll0(GEN x, long flag)
    2069             : {
    2070         301 :   if (typ(x) != t_MAT) pari_err_TYPE("qflll",x);
    2071         301 :   switch(flag)
    2072             :   {
    2073          49 :     case 0: return lll(x);
    2074          63 :     case 1: RgM_check_ZM(x,"qflll"); return lllint(x);
    2075          49 :     case 2: RgM_check_ZM(x,"qflll"); return lllintpartial(x);
    2076           7 :     case 3: RgM_check_ZM(x,"qflll"); return lllall(x, LLL_INPLACE);
    2077          49 :     case 4: RgM_check_ZM(x,"qflll"); return lllkerim(x);
    2078          42 :     case 5: return lllkerimgen(x);
    2079          42 :     case 8: return lllgen(x);
    2080           0 :     default: pari_err_FLAG("qflll");
    2081             :   }
    2082             :   return NULL; /* LCOV_EXCL_LINE */
    2083             : }
    2084             : 
    2085             : GEN
    2086         203 : qflllgram0(GEN x, long flag)
    2087             : {
    2088         203 :   if (typ(x) != t_MAT) pari_err_TYPE("qflllgram",x);
    2089         203 :   switch(flag)
    2090             :   {
    2091          63 :     case 0: return lllgram(x);
    2092          49 :     case 1: RgM_check_ZM(x,"qflllgram"); return lllgramint(x);
    2093          49 :     case 4: RgM_check_ZM(x,"qflllgram"); return lllgramkerim(x);
    2094          42 :     case 5: return lllgramkerimgen(x);
    2095           0 :     case 8: return lllgramgen(x);
    2096           0 :     default: pari_err_FLAG("qflllgram");
    2097             :   }
    2098             :   return NULL; /* LCOV_EXCL_LINE */
    2099             : }
    2100             : 
    2101             : /********************************************************************/
    2102             : /**                                                                **/
    2103             : /**                   INTEGRAL KERNEL (LLL REDUCED)                **/
    2104             : /**                                                                **/
    2105             : /********************************************************************/
    2106             : static GEN
    2107          70 : kerint0(GEN M)
    2108             : {
    2109             :   /* return ZM_lll(M, LLLDFT, LLL_KER); */
    2110          70 :   GEN U, H = ZM_hnflll(M,&U,1);
    2111          70 :   long d = lg(M)-lg(H);
    2112          70 :   if (!d) return cgetg(1, t_MAT);
    2113          70 :   return ZM_lll(vecslice(U,1,d), LLLDFT, LLL_INPLACE);
    2114             : }
    2115             : GEN
    2116          42 : kerint(GEN M)
    2117             : {
    2118          42 :   pari_sp av = avma;
    2119          42 :   return gerepilecopy(av, kerint0(M));
    2120             : }
    2121             : /* OBSOLETE: use kerint */
    2122             : GEN
    2123          28 : matkerint0(GEN M, long flag)
    2124             : {
    2125          28 :   pari_sp av = avma;
    2126          28 :   if (typ(M) != t_MAT) pari_err_TYPE("matkerint",M);
    2127          28 :   M = Q_primpart(M);
    2128          28 :   RgM_check_ZM(M, "kerint");
    2129          28 :   switch(flag)
    2130             :   {
    2131          28 :     case 0:
    2132          28 :     case 1: return gerepilecopy(av, kerint0(M));
    2133           0 :     default: pari_err_FLAG("matkerint");
    2134             :   }
    2135             :   return NULL; /* LCOV_EXCL_LINE */
    2136             : }

Generated by: LCOV version 1.13