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.16.2 lcov report (development 29115-f22e516b23) Lines: 1301 1607 81.0 %
Date: 2024-04-22 08:08:03 Functions: 121 124 97.6 %
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 int
      21       50537 : RgM_is_square_mat(GEN x) { long l = lg(x); return l == 1 || l == lgcols(x); }
      22             : 
      23             : static long
      24     4223928 : ZM_is_upper(GEN R)
      25             : {
      26     4223928 :   long i,j, l = lg(R);
      27     4223928 :   if (l != lgcols(R)) return 0;
      28     8164744 :   for(i = 1; i < l; i++)
      29     8855405 :     for(j = 1; j < i; j++)
      30     4574803 :       if (signe(gcoeff(R,i,j))) return 0;
      31      251379 :   return 1;
      32             : }
      33             : 
      34             : static long
      35      610140 : ZM_is_knapsack(GEN R)
      36             : {
      37      610140 :   long i,j, l = lg(R);
      38      610140 :   if (l != lgcols(R)) return 0;
      39      857818 :   for(i = 2; i < l; i++)
      40     3044835 :     for(j = 1; j < l; j++)
      41     2797157 :       if ( i!=j && signe(gcoeff(R,i,j))) return 0;
      42       95146 :   return 1;
      43             : }
      44             : 
      45             : static long
      46     1191166 : ZM_is_lower(GEN R)
      47             : {
      48     1191166 :   long i,j, l = lg(R);
      49     1191166 :   if (l != lgcols(R)) return 0;
      50     2078498 :   for(i = 1; i < l; i++)
      51     2395699 :     for(j = 1; j < i; j++)
      52     1294627 :       if (signe(gcoeff(R,j,i))) return 0;
      53       34647 :   return 1;
      54             : }
      55             : 
      56             : static GEN
      57       34648 : RgM_flip(GEN R)
      58             : {
      59             :   GEN M;
      60             :   long i,j,l;
      61       34648 :   M = cgetg_copy(R, &l);
      62      179283 :   for(i = 1; i < l; i++)
      63             :   {
      64      144635 :     gel(M,i) = cgetg(l, t_COL);
      65      890768 :     for(j = 1; j < l; j++)
      66      746133 :       gmael(M,i,j) = gmael(R,l-i, l-j);
      67             :   }
      68       34648 :   return M;
      69             : }
      70             : 
      71             : static GEN
      72           0 : RgM_flop(GEN R)
      73             : {
      74             :   GEN M;
      75             :   long i,j,l;
      76           0 :   M = cgetg_copy(R, &l);
      77           0 :   for(i = 1; i < l; i++)
      78             :   {
      79           0 :     gel(M,i) = cgetg(l, t_COL);
      80           0 :     for(j = 1; j < l; j++)
      81           0 :       gmael(M,i,j) = gmael(R,i, l-j);
      82             :   }
      83           0 :   return M;
      84             : }
      85             : 
      86             : /* Assume x and y has same type! */
      87             : INLINE int
      88    23089076 : mpabscmp(GEN x, GEN y)
      89             : {
      90    23089076 :   return (typ(x)==t_INT) ? abscmpii(x,y) : abscmprr(x,y);
      91             : }
      92             : 
      93             : /****************************************************************************/
      94             : /***                             FLATTER                                  ***/
      95             : /****************************************************************************/
      96             : 
      97             : /* Implementation of "FLATTER" algorithm based on
      98             :    <https://eprint.iacr.org/2023/237>
      99             :    Fast Practical Lattice Reduction through Iterated Compression
     100             : 
     101             :    Keegan Ryan, University of California, San Diego
     102             :    Nadia Heninger, University of California, San Diego
     103             : BA20230925 */
     104             : 
     105             : static long
     106     5080754 : drop(GEN R)
     107             : {
     108     5080754 :   long i, n = lg(R)-1;
     109     5080754 :   long s = 0, m = mpexpo(gcoeff(R, 1, 1));
     110    20798056 :   for (i = 2; i <= n; ++i)
     111             :   {
     112    15717304 :     if (mpabscmp(gcoeff(R, i, i), gcoeff(R, i - 1, i - 1)) >= 0)
     113             :     {
     114     9597343 :       s += m - mpexpo(gcoeff(R, i - 1, i - 1));
     115     9597343 :       m = mpexpo(gcoeff(R, i, i));
     116             :     }
     117             :   }
     118     5080752 :   s += m - mpexpo(gcoeff(R, n, n));
     119     5080752 :   return s;
     120             : }
     121             : 
     122             : static long
     123     4183170 : gsisinv(GEN M)
     124             : {
     125     4183170 :   long i, l = lg(M);
     126    21602115 :   for (i = 1; i < l; ++i)
     127    17419349 :     if (signe(gmael(M, i, i))==0)
     128         404 :       return 0;
     129     4182766 :   return 1;
     130             : }
     131             : 
     132             : INLINE long
     133     6488465 : nbits2prec64(long n)
     134             : {
     135     6488465 :   return nbits2prec(((n+63)>>6)<<6);
     136             : }
     137             : 
     138             : static GEN
     139        4057 : RgM_Cholesky_dynprec(GEN M)
     140             : {
     141        4057 :   pari_sp ltop = avma;
     142             :   GEN L;
     143        4057 :   long minprec = 3*(lg(M)-1) + 30, bitprec = minprec, prec;
     144             :   while (1)
     145        4307 :   {
     146             :     long mbitprec;
     147        8364 :     prec = nbits2prec64(bitprec);
     148        8364 :     L = RgM_Cholesky(RgM_gtofp(M, prec), prec);
     149        8364 :     if (!L)
     150             :     {
     151        1807 :       bitprec *= 2;
     152        1807 :       set_avma(ltop);
     153        1807 :       continue;
     154             :     }
     155        6557 :     mbitprec = minprec + 2*drop(L);
     156        6557 :     if (bitprec >= mbitprec)
     157        4057 :       break;
     158        2500 :     bitprec = mbitprec;
     159        2500 :     set_avma(ltop);
     160             :   }
     161        4057 :   return gerepilecopy(ltop, L);
     162             : }
     163             : 
     164             : static GEN
     165        1455 : gramschmidt_upper(GEN M)
     166             : {
     167        1455 :   long bitprec = 3*(lg(M)-1) + 30 + 2*drop(M);
     168        1455 :   return RgM_gtofp(M, nbits2prec64(bitprec));
     169             : }
     170             : 
     171             : static GEN
     172     2709990 : gramschmidt_dynprec(GEN M)
     173             : {
     174     2709990 :   pari_sp ltop = avma;
     175     2709990 :   long minprec = 3*(lg(M)-1) + 30, bitprec = minprec;
     176     2709990 :   if (ZM_is_upper(M)) return gramschmidt_upper(M);
     177             :   while (1)
     178     2626613 :   {
     179             :     GEN B, Q, L;
     180     5335148 :     long prec = nbits2prec64(bitprec), mbitprec;
     181     5335146 :     if (!QR_init(RgM_gtofp(M, prec), &B, &Q, &L, prec) || !gsisinv(L))
     182             :     {
     183     1621445 :       bitprec *= 2;
     184     1621445 :       set_avma(ltop);
     185     1621455 :       continue;
     186             :     }
     187     3713691 :     mbitprec = minprec + 2*drop(L);
     188     3713691 :     if (bitprec >= mbitprec)
     189     2708533 :       return gerepilecopy(ltop, shallowtrans(L));
     190     1005158 :     bitprec = mbitprec;
     191     1005158 :     set_avma(ltop);
     192             :   }
     193             : }
     194             : /* return -T1 * round(T1^-1*(R1^-1*R2)*T3) */
     195             : static GEN
     196     1354996 : sizered(GEN T1, GEN T3, GEN R1, GEN R2)
     197             : {
     198     1354996 :   pari_sp ltop = avma;
     199             :   long e;
     200     1354996 :   return gerepileupto(ltop, ZM_mul(ZM_neg(T1), grndtoi(gmul(ZM_inv(T1,NULL),
     201             :          RgM_mul(RgM_mul(RgM_inv_upper(R1), R2), T3)), &e)));
     202             : }
     203             : 
     204             : static GEN
     205     1354996 : flat(GEN M, long flag, GEN *pt_T, long *pt_s)
     206             : {
     207     1354996 :   pari_sp ltop = avma;
     208             :   GEN R, R1, R2, R3, T1, T2, T3, T, S;
     209     1354996 :   long k = lg(M)-1, n = k>>1, n2 = k - n, m = n>>1;
     210     1354996 :   long keepfirst = flag & LLL_KEEP_FIRST, inplace = flag & LLL_INPLACE;
     211             :   /* for k = 3, we want n = 1; n2  = 2; m = 0 */
     212             :   /* for k = 5,         n = 2; n2 = 3; m = 1 */
     213     1354996 :   R = gramschmidt_dynprec(M);
     214     1354996 :   R1 = matslice(R, 1, n, 1, n);
     215     1354995 :   R2 = matslice(R, 1, n, n + 1, k);
     216     1354995 :   R3 = matslice(R, n + 1, k, n + 1, k);
     217     1354996 :   T1 = lllfp(R1, 0.99, LLL_IM| LLL_UPPER| LLL_NOCERTIFY| (keepfirst ? LLL_KEEP_FIRST: 0));
     218     1354995 :   T3 = lllfp(R3, 0.99, LLL_IM| LLL_UPPER| LLL_NOCERTIFY);
     219     1354996 :   T2 = sizered(T1, T3, R1, R2);
     220     1354996 :   T = shallowmatconcat(mkmat22(T1,T2,gen_0,T3));
     221     1354996 :   M = ZM_mul(M, T);
     222     1354994 :   R = gramschmidt_dynprec(M);
     223     1354996 :   R3 = matslice(R, m + 1, m + n2, m + 1, m + n2);
     224     1354996 :   T3 = lllfp(R3, 0.99, LLL_IM| LLL_UPPER| LLL_NOCERTIFY);
     225     2709991 :   S = shallowmatconcat(diagonal(
     226      575795 :        m == 0     ? mkvec2(T3, matid(k - m - n2))
     227           0 :      : m+n2 == k  ? mkvec2(matid(m), T3)
     228      779201 :                   : mkvec3(matid(m), T3, matid(k - m - n2))));
     229     1354996 :   M = ZM_mul(M, S);
     230     1354996 :   if (!inplace) *pt_T = ZM_mul(T, S);
     231     1354994 :   *pt_s = drop(R);
     232     1354994 :   return gc_all(ltop, inplace ? 1: 2, &M, pt_T);
     233             : }
     234             : 
     235             : static GEN
     236      628010 : ZM_flatter(GEN M, long flag)
     237             : {
     238      628010 :   pari_sp ltop = avma, btop;
     239      628010 :   long i, n = lg(M)-1;
     240      628010 :   long s = -1;
     241             :   GEN T;
     242             :   pari_timer ti;
     243      628010 :   long lti = 1;
     244      628010 :   long inplace = flag & LLL_INPLACE, cert = !(flag & LLL_NOCERTIFY);
     245      628010 :   T = matid(n);
     246      628010 :   if (DEBUGLEVEL>=3)
     247             :   {
     248           0 :     timer_start(&ti);
     249           0 :     if (cert)
     250           0 :       err_printf("flatter dim = %ld size = %ld\n", n, ZM_max_expi(M));
     251             :   }
     252      628010 :   btop = avma;
     253      628010 :   for (i = 1;;i++)
     254      726986 :   {
     255             :     long t;
     256     1354996 :     GEN U, M2 = flat(M, flag, &U, &t);
     257     1354995 :     if (t==0) { s = t; break; }
     258      771008 :     if (s >= 0)
     259             :     {
     260      440570 :       if (s==t)
     261       44023 :         break;
     262      396547 :       if (s<t && i > 20)
     263             :       {
     264           0 :         if (DEBUGLEVEL >= 3) err_printf("BACK:%ld:%ld:%g\n", n, i, s);
     265           0 :         break;
     266             :       }
     267             :     }
     268      726985 :     if (DEBUGLEVEL>=3 && (cert || timer_get(&ti) > 1000))
     269             :     {
     270           0 :       if (i==lti)
     271           0 :         timer_printf(&ti, "FLATTER, dim %ld, step %ld: \t slope=%0.10g", n, i, ((double)t)/n);
     272             :       else
     273           0 :         timer_printf(&ti, "FLATTER, dim %ld, steps %ld-%ld: \t slope=%0.10g", n, lti,i,((double)t)/n);
     274           0 :       lti = i+1;
     275             :     }
     276      726985 :     s = t;
     277      726985 :     M = M2;
     278      726985 :     if (!inplace) T = ZM_mul(T, U);
     279      726985 :     if (gc_needed(btop, 1))
     280           0 :       gerepileall(btop, 2, &M, &T);
     281             :   }
     282      628010 :   if (DEBUGLEVEL>=3 && (cert || timer_get(&ti) > 1000))
     283             :   {
     284           0 :     if (i==lti)
     285           0 :       timer_printf(&ti, "FLATTER, dim %ld, final \t slope=%0.10g", n, ((double)s)/n);
     286             :     else
     287           0 :       timer_printf(&ti, "FLATTER, dim %ld, steps %ld-final:\t slope=%0.10g", n, lti, ((double)s)/n);
     288             :   }
     289      628010 :   return  gerepilecopy(ltop, inplace ? M: T);
     290             : }
     291             : 
     292             : static GEN
     293      626029 : ZM_flatter_rank(GEN M, long rank, long flag)
     294             : {
     295             :   pari_timer ti;
     296      626029 :   pari_sp ltop = avma;
     297             :   GEN T;
     298      626029 :   long i, n = lg(M)-1;
     299      626029 :   if (rank == n)  return ZM_flatter(M, flag);
     300        3764 :   T = matid(n);
     301        3764 :   if (DEBUGLEVEL>=3) timer_start(&ti);
     302        3764 :   for (i = 1;; i++)
     303        1981 :   {
     304        5745 :     GEN S = ZM_flatter(vconcat(gshift(M,i),matid(n)), flag);
     305        5745 :     if (DEBUGLEVEL>=3)
     306           0 :       timer_printf(&ti,"FLATTERRANK step %ld: %ld",i,expi(gnorml2(S)));
     307        5745 :     if (ZM_isidentity(S)) break;
     308        1981 :     T = ZM_mul(T, S);
     309        1981 :     M = ZM_mul(M, S);
     310        1981 :     if (gc_needed(ltop, 1))
     311           0 :       gerepileall(ltop, 2, &M, &T);
     312             :   }
     313        3764 :   return gerepilecopy(ltop, T);
     314             : }
     315             : 
     316             : static GEN
     317        4057 : flattergram_i(GEN M, long flag, long *pt_s)
     318             : {
     319        4057 :   pari_sp ltop = avma;
     320             :   GEN R, T;
     321        4057 :   R = RgM_Cholesky_dynprec(M);
     322        4057 :   *pt_s = drop(R);
     323        4057 :   T =  lllfp(R, 0.99, LLL_IM| LLL_UPPER| LLL_NOCERTIFY | (flag&LLL_KEEP_FIRST));
     324        4057 :   return gerepilecopy(ltop, T);
     325             : }
     326             : 
     327             : static GEN
     328        1322 : ZM_flattergram(GEN M, long flag)
     329             : {
     330        1322 :   pari_sp ltop = avma, btop;
     331             :   GEN T;
     332        1322 :   long i, n = lg(M)-1;
     333             :   pari_timer ti;
     334        1322 :   long s = -1;
     335        1322 :   if (DEBUGLEVEL>=3)
     336             :   {
     337           0 :     timer_start(&ti);
     338           0 :     err_printf("FLATTERGRAM dim = %ld size = %ld\n", n, expi(gnorml2(M)));
     339             :   }
     340        1322 :   btop = avma;
     341        1322 :   T = matid(n);
     342        1322 :   for (i = 1;; i++)
     343        2735 :   {
     344             :     long t;
     345        4057 :     GEN S = flattergram_i(M, flag, &t);
     346        4057 :     t = expi(gnorml2(S));
     347        4057 :     if (t==0) { s = t;  break; }
     348        4057 :     if (s)
     349             :     {
     350        4057 :       double st = s-t;
     351        4057 :       if (st == 0)
     352        1322 :         break;
     353        2735 :       if (st < 0 && i > 20)
     354             :       {
     355           0 :         if (DEBUGLEVEL >= 3) err_printf("BACK:%ld:%ld:%0.10g\n", n, i, ((double)s)/n);
     356           0 :         break;
     357             :       }
     358             :     }
     359        2735 :     T = ZM_mul(T, S);
     360        2735 :     M = ZM_mul(shallowtrans(S), ZM_mul(M, S));
     361        2735 :     s = t;
     362        2735 :     if (DEBUGLEVEL >= 3)
     363           0 :       timer_printf(&ti, "FLATTERGRAM, dim %ld step %ld, slope=%0.10g", n, i, ((double)s)/n);
     364        2735 :     if (gc_needed(btop, 1))
     365           0 :       gerepileall(btop, 2, &M, &T);
     366             :   }
     367        1322 :   if (DEBUGLEVEL >= 3)
     368           0 :     timer_printf(&ti, "FLATTERGRAM, dim %ld, slope=%0.10g", n, ((double)s)/n);
     369        1322 :   return gerepilecopy(ltop, T);
     370             : }
     371             : 
     372             : static GEN
     373        1322 : ZM_flattergram_rank(GEN M, long rank, long flag)
     374             : {
     375             :   pari_timer ti;
     376        1322 :   pari_sp ltop = avma;
     377             :   GEN T;
     378        1322 :   long i, n = lg(M)-1;
     379        1322 :   if (rank == n)  return ZM_flattergram(M, flag);
     380           0 :   T = matid(n);
     381           0 :   if (DEBUGLEVEL>=3) timer_start(&ti);
     382           0 :   for (i = 1;; i++)
     383           0 :   {
     384           0 :     GEN S = ZM_flattergram(RgM_Rg_add(gshift(M, i), gen_1), flag);
     385           0 :     if (DEBUGLEVEL>=3)
     386           0 :       timer_printf(&ti,"FLATTERGRAMRANK step %ld: %ld",i,expi(gnorml2(S)));
     387           0 :     if (ZM_isidentity(S)) break;
     388           0 :     T = ZM_mul(T, S);
     389           0 :     M = ZM_mul(shallowtrans(S), ZM_mul(M, S));
     390           0 :     if (gc_needed(ltop, 1))
     391           0 :       gerepileall(ltop, 2, &M, &T);
     392             :   }
     393           0 :   return gerepilecopy(ltop, T);
     394             : }
     395             : 
     396             : static double
     397    10935199 : pari_rint(double a)
     398             : {
     399             : #ifdef HAS_RINT
     400    10935199 :   return rint(a);
     401             : #else
     402             :   const double two_to_52 = 4.5035996273704960e+15;
     403             :   double fa = fabs(a);
     404             :   double r = two_to_52 + fa;
     405             :   if (fa >= two_to_52) {
     406             :     r = a;
     407             :   } else {
     408             :     r = r - two_to_52;
     409             :     if (a < 0) r = -r;
     410             :   }
     411             :   return r;
     412             : #endif
     413             : }
     414             : 
     415             : /* default quality ratio for LLL */
     416             : static const double LLLDFT = 0.99;
     417             : 
     418             : /* assume flag & (LLL_KER|LLL_IM|LLL_ALL). LLL_INPLACE implies LLL_IM */
     419             : static GEN
     420      768771 : lll_trivial(GEN x, long flag)
     421             : {
     422      768771 :   if (lg(x) == 1)
     423             :   { /* dim x = 0 */
     424       15449 :     if (! (flag & LLL_ALL)) return cgetg(1,t_MAT);
     425          28 :     retmkvec2(cgetg(1,t_MAT), cgetg(1,t_MAT));
     426             :   }
     427             :   /* dim x = 1 */
     428      753322 :   if (gequal0(gel(x,1)))
     429             :   {
     430          91 :     if (flag & LLL_KER) return matid(1);
     431          91 :     if (flag & (LLL_IM|LLL_INPLACE)) return cgetg(1,t_MAT);
     432          28 :     retmkvec2(matid(1), cgetg(1,t_MAT));
     433             :   }
     434      753227 :   if (flag & LLL_INPLACE) return gcopy(x);
     435      649899 :   if (flag & LLL_KER) return cgetg(1,t_MAT);
     436      649899 :   if (flag & LLL_IM)  return matid(1);
     437          28 :   retmkvec2(cgetg(1,t_MAT), (flag & LLL_GRAM)? gcopy(x): matid(1));
     438             : }
     439             : 
     440             : /* vecslice(x,#x-k,#x) in place. Works for t_MAT, t_VEC/t_COL */
     441             : static GEN
     442     2067603 : vectail_inplace(GEN x, long k)
     443             : {
     444     2067603 :   if (!k) return x;
     445       57604 :   x[k] = ((ulong)x[0] & ~LGBITS) | _evallg(lg(x) - k);
     446       57604 :   return x + k;
     447             : }
     448             : 
     449             : /* k = dim Kernel */
     450             : static GEN
     451     2140689 : lll_finish(GEN h, long k, long flag)
     452             : {
     453             :   GEN g;
     454     2140689 :   if (!(flag & (LLL_IM|LLL_KER|LLL_ALL|LLL_INPLACE))) return h;
     455     2067632 :   if (flag & (LLL_IM|LLL_INPLACE)) return vectail_inplace(h, k);
     456          98 :   if (flag & LLL_KER) { setlg(h,k+1); return h; }
     457          70 :   g = vecslice(h,1,k); /* done first: vectail_inplace kills h */
     458          70 :   return mkvec2(g, vectail_inplace(h, k));
     459             : }
     460             : 
     461             : /* y * z * 2^e, e >= 0; y,z t_INT */
     462             : INLINE GEN
     463      440248 : mulshift(GEN y, GEN z, long e)
     464             : {
     465      440248 :   long ly = lgefint(y), lz;
     466             :   pari_sp av;
     467             :   GEN t;
     468      440248 :   if (ly == 2) return gen_0;
     469       55256 :   lz = lgefint(z);
     470       55256 :   av = avma; (void)new_chunk(ly+lz+nbits2lg(e)); /* HACK */
     471       55256 :   t = mulii(z, y);
     472       55256 :   set_avma(av); return shifti(t, e);
     473             : }
     474             : 
     475             : /* x - y * z * 2^e, e >= 0; x,y,z t_INT */
     476             : INLINE GEN
     477     2075086 : submulshift(GEN x, GEN y, GEN z, long e)
     478             : {
     479     2075086 :   long lx = lgefint(x), ly, lz;
     480             :   pari_sp av;
     481             :   GEN t;
     482     2075086 :   if (!e) return submulii(x, y, z);
     483     2055591 :   if (lx == 2) { t = mulshift(y, z, e); togglesign(t); return t; }
     484     1639513 :   ly = lgefint(y);
     485     1639513 :   if (ly == 2) return icopy(x);
     486     1184972 :   lz = lgefint(z);
     487     1184972 :   av = avma; (void)new_chunk(lx+ly+lz+nbits2lg(e)); /* HACK */
     488     1184972 :   t = shifti(mulii(z, y), e);
     489     1184972 :   set_avma(av); return subii(x, t);
     490             : }
     491             : static void
     492    18624585 : subzi(GEN *a, GEN b)
     493             : {
     494    18624585 :   pari_sp av = avma;
     495    18624585 :   b = subii(*a, b);
     496    18624609 :   if (lgefint(b)<=lg(*a) && isonstack(*a)) { affii(b,*a); set_avma(av); }
     497     2128138 :   else *a = b;
     498    18624663 : }
     499             : 
     500             : static void
     501    17909523 : addzi(GEN *a, GEN b)
     502             : {
     503    17909523 :   pari_sp av = avma;
     504    17909523 :   b = addii(*a, b);
     505    17909536 :   if (lgefint(b)<=lg(*a) && isonstack(*a)) { affii(b,*a); set_avma(av); }
     506     1908945 :   else *a = b;
     507    17909573 : }
     508             : 
     509             : /* x - u*y * 2^e */
     510             : INLINE GEN
     511     4220299 : submuliu2n(GEN x, GEN y, ulong u, long e)
     512             : {
     513             :   pari_sp av;
     514     4220299 :   long ly = lgefint(y);
     515     4220299 :   if (ly == 2) return x;
     516     2863568 :   av = avma;
     517     2863568 :   (void)new_chunk(3+ly+lgefint(x)+nbits2lg(e)); /* HACK */
     518     2865235 :   y = shifti(mului(u,y), e);
     519     2864179 :   set_avma(av); return subii(x, y);
     520             : }
     521             : /* *x -= u*y * 2^e */
     522             : INLINE void
     523      320022 : submulzu2n(GEN *x, GEN y, ulong u, long e)
     524             : {
     525             :   pari_sp av;
     526      320022 :   long ly = lgefint(y);
     527      320022 :   if (ly == 2) return;
     528      209921 :   av = avma;
     529      209921 :   (void)new_chunk(3+ly+lgefint(*x)+nbits2lg(e)); /* HACK */
     530      209921 :   y = shifti(mului(u,y), e);
     531      209921 :   set_avma(av); return subzi(x, y);
     532             : }
     533             : 
     534             : /* x + u*y * 2^e */
     535             : INLINE GEN
     536     4100390 : addmuliu2n(GEN x, GEN y, ulong u, long e)
     537             : {
     538             :   pari_sp av;
     539     4100390 :   long ly = lgefint(y);
     540     4100390 :   if (ly == 2) return x;
     541     2798185 :   av = avma;
     542     2798185 :   (void)new_chunk(3+ly+lgefint(x)+nbits2lg(e)); /* HACK */
     543     2799670 :   y = shifti(mului(u,y), e);
     544     2798680 :   set_avma(av); return addii(x, y);
     545             : }
     546             : 
     547             : /* *x += u*y * 2^e */
     548             : INLINE void
     549      326940 : addmulzu2n(GEN *x, GEN y, ulong u, long e)
     550             : {
     551             :   pari_sp av;
     552      326940 :   long ly = lgefint(y);
     553      326940 :   if (ly == 2) return;
     554      212692 :   av = avma;
     555      212692 :   (void)new_chunk(3+ly+lgefint(*x)+nbits2lg(e)); /* HACK */
     556      212692 :   y = shifti(mului(u,y), e);
     557      212692 :   set_avma(av); return addzi(x, y);
     558             : }
     559             : 
     560             : /* n < 10; gerepileall supporting &NULL arguments. Maybe rename and export ? */
     561             : INLINE void
     562        4640 : gc_lll(pari_sp av, int n, ...)
     563             : {
     564             :   int i, j;
     565             :   GEN *gptr[10];
     566             :   size_t s;
     567        4640 :   va_list a; va_start(a, n);
     568       13920 :   for (i=j=0; i<n; i++)
     569             :   {
     570        9280 :     GEN *x = va_arg(a,GEN*);
     571        9280 :     if (*x) { gptr[j++] = x; *x = (GEN)copy_bin(*x); }
     572             :   }
     573        4640 :   va_end(a); set_avma(av);
     574       11539 :   for (--j; j>=0; j--) *gptr[j] = bin_copy((GENbin*)*gptr[j]);
     575        4640 :   s = pari_mainstack->top - pari_mainstack->bot;
     576             :   /* size of saved objects ~ stacksize / 4 => overflow */
     577        4640 :   if (av - avma > (s >> 2))
     578             :   {
     579           0 :     size_t t = avma - pari_mainstack->bot;
     580           0 :     av = avma; new_chunk((s + t) / sizeof(long)); set_avma(av); /* double */
     581             :   }
     582        4640 : }
     583             : 
     584             : /********************************************************************/
     585             : /**                                                                **/
     586             : /**                   FPLLL (adapted from D. Stehle's code)        **/
     587             : /**                                                                **/
     588             : /********************************************************************/
     589             : /* Babai* and fplll* are a conversion to libpari API and data types
     590             :    of fplll-1.3 by Damien Stehle'.
     591             : 
     592             :   Copyright 2005, 2006 Damien Stehle'.
     593             : 
     594             :   This program is free software; you can redistribute it and/or modify it
     595             :   under the terms of the GNU General Public License as published by the
     596             :   Free Software Foundation; either version 2 of the License, or (at your
     597             :   option) any later version.
     598             : 
     599             :   This program implements ideas from the paper "Floating-point LLL Revisited",
     600             :   by Phong Nguyen and Damien Stehle', in the Proceedings of Eurocrypt'2005,
     601             :   Springer-Verlag; and was partly inspired by Shoup's NTL library:
     602             :   http://www.shoup.net/ntl/ */
     603             : 
     604             : /* x t_REAL, |x| >= 1/2. Test whether |x| <= 3/2 */
     605             : static int
     606      469696 : absrsmall2(GEN x)
     607             : {
     608      469696 :   long e = expo(x), l, i;
     609      469696 :   if (e < 0) return 1;
     610      224807 :   if (e > 0 || (ulong)x[2] > (3UL << (BITS_IN_LONG-2))) return 0;
     611             :   /* line above assumes l > 2. OK since x != 0 */
     612       77329 :   l = lg(x); for (i = 3; i < l; i++) if (x[i]) return 0;
     613       66467 :   return 1;
     614             : }
     615             : /* x t_REAL; test whether |x| <= 1/2 */
     616             : static int
     617      883924 : absrsmall(GEN x)
     618             : {
     619             :   long e, l, i;
     620      883924 :   if (!signe(x)) return 1;
     621      878279 :   e = expo(x); if (e < -1) return 1;
     622      476293 :   if (e > -1 || (ulong)x[2] > HIGHBIT) return 0;
     623        7528 :   l = lg(x); for (i = 3; i < l; i++) if (x[i]) return 0;
     624        6597 :   return 1;
     625             : }
     626             : 
     627             : static void
     628    31957320 : rotate(GEN A, long k2, long k)
     629             : {
     630             :   long i;
     631    31957320 :   GEN B = gel(A,k2);
     632   101032629 :   for (i = k2; i > k; i--) gel(A,i) = gel(A,i-1);
     633    31957320 :   gel(A,k) = B;
     634    31957320 : }
     635             : 
     636             : /************************* FAST version (double) ************************/
     637             : #define dmael(x,i,j) ((x)[i][j])
     638             : #define del(x,i) ((x)[i])
     639             : 
     640             : static double *
     641    34692068 : cget_dblvec(long d)
     642    34692068 : { return (double*) stack_malloc_align(d*sizeof(double), sizeof(double)); }
     643             : 
     644             : static double **
     645     8320580 : cget_dblmat(long d) { return (double **) cgetg(d, t_VECSMALL); }
     646             : 
     647             : static double
     648   159992823 : itodbl_exp(GEN x, long *e)
     649             : {
     650   159992823 :   pari_sp av = avma;
     651   159992823 :   GEN r = itor(x,DEFAULTPREC);
     652   159985172 :   *e = expo(r); setexpo(r,0);
     653   159981735 :   return gc_double(av, rtodbl(r));
     654             : }
     655             : 
     656             : static double
     657   116338381 : dbldotproduct(double *x, double *y, long n)
     658             : {
     659             :   long i;
     660   116338381 :   double sum = del(x,1) * del(y,1);
     661  1366417109 :   for (i=2; i<=n; i++) sum += del(x,i) * del(y,i);
     662   116338381 :   return sum;
     663             : }
     664             : 
     665             : static double
     666     2449260 : dbldotsquare(double *x, long n)
     667             : {
     668             :   long i;
     669     2449260 :   double sum = del(x,1) * del(x,1);
     670     8147097 :   for (i=2; i<=n; i++) sum += del(x,i) * del(x,i);
     671     2449260 :   return sum;
     672             : }
     673             : 
     674             : static long
     675    24738601 : set_line(double *appv, GEN v, long n)
     676             : {
     677    24738601 :   long i, maxexp = 0;
     678    24738601 :   pari_sp av = avma;
     679    24738601 :   GEN e = cgetg(n+1, t_VECSMALL);
     680   184719379 :   for (i = 1; i <= n; i++)
     681             :   {
     682   159990432 :     del(appv,i) = itodbl_exp(gel(v,i), e+i);
     683   159979473 :     if (e[i] > maxexp) maxexp = e[i];
     684             :   }
     685   184806299 :   for (i = 1; i <= n; i++) del(appv,i) = ldexp(del(appv,i), e[i]-maxexp);
     686    24728947 :   set_avma(av); return maxexp;
     687             : }
     688             : 
     689             : static void
     690    34605904 : dblrotate(double **A, long k2, long k)
     691             : {
     692             :   long i;
     693    34605904 :   double *B = del(A,k2);
     694   108507779 :   for (i = k2; i > k; i--) del(A,i) = del(A,i-1);
     695    34605904 :   del(A,k) = B;
     696    34605904 : }
     697             : /* update G[kappa][i] from appB */
     698             : static void
     699    22499969 : setG_fast(double **appB, long n, double **G, long kappa, long a, long b)
     700             : { long i;
     701   100648858 :   for (i = a; i <= b; i++)
     702    78149074 :     dmael(G,kappa,i) = dbldotproduct(del(appB,kappa), del(appB,i), n);
     703    22499784 : }
     704             : /* update G[i][kappa] from appB */
     705             : static void
     706    16883592 : setG2_fast(double **appB, long n, double **G, long kappa, long a, long b)
     707             : { long i;
     708    55077018 :   for (i = a; i <= b; i++)
     709    38193392 :     dmael(G,i,kappa) = dbldotproduct(del(appB,kappa), del(appB,i), n);
     710    16883626 : }
     711             : const long EX0 = -2; /* uninitialized; any value less than expo(0.51) = -1 */
     712             : 
     713             : #ifdef LONG_IS_64BIT
     714             : typedef long s64;
     715             : #define addmuliu64_inplace addmuliu_inplace
     716             : #define submuliu64_inplace submuliu_inplace
     717             : #define submuliu642n submuliu2n
     718             : #define addmuliu642n addmuliu2n
     719             : #else
     720             : typedef long long s64;
     721             : typedef unsigned long long u64;
     722             : 
     723             : INLINE GEN
     724    20062829 : u64toi(u64 x)
     725             : {
     726             :   GEN y;
     727             :   ulong h;
     728    20062829 :   if (!x) return gen_0;
     729    20062829 :   h = x>>32;
     730    20062829 :   if (!h) return utoipos(x);
     731     1149255 :   y = cgetipos(4);
     732     1149255 :   *int_LSW(y) = x&0xFFFFFFFF;
     733     1149255 :   *int_MSW(y) = x>>32;
     734     1149255 :   return y;
     735             : }
     736             : 
     737             : INLINE GEN
     738      665952 : u64toineg(u64 x)
     739             : {
     740             :   GEN y;
     741             :   ulong h;
     742      665952 :   if (!x) return gen_0;
     743      665952 :   h = x>>32;
     744      665952 :   if (!h) return utoineg(x);
     745      665952 :   y = cgetineg(4);
     746      665952 :   *int_LSW(y) = x&0xFFFFFFFF;
     747      665952 :   *int_MSW(y) = x>>32;
     748      665952 :   return y;
     749             : }
     750             : INLINE GEN
     751     9641715 : addmuliu64_inplace(GEN x, GEN y, u64 u) { return addmulii(x, y, u64toi(u)); }
     752             : 
     753             : INLINE GEN
     754     9732654 : submuliu64_inplace(GEN x, GEN y, u64 u) { return submulii(x, y, u64toi(u)); }
     755             : 
     756             : INLINE GEN
     757      665952 : addmuliu642n(GEN x, GEN y, u64 u, long e) { return submulshift(x, y, u64toineg(u), e); }
     758             : 
     759             : INLINE GEN
     760      688460 : submuliu642n(GEN x, GEN y, u64 u, long e) { return submulshift(x, y, u64toi(u), e); }
     761             : 
     762             : #endif
     763             : 
     764             : /* Babai's Nearest Plane algorithm (iterative); see Babai() */
     765             : static int
     766    29811660 : Babai_fast(pari_sp av, long kappa, GEN *pB, GEN *pU, double **mu, double **r,
     767             :            double *s, double **appB, GEN expoB, double **G,
     768             :            long a, long zeros, long maxG, double eta)
     769             : {
     770    29811660 :   GEN B = *pB, U = *pU;
     771    29811660 :   const long n = nbrows(B), d = U ? lg(U)-1: 0;
     772    29811505 :   long k, aa = (a > zeros)? a : zeros+1;
     773    29811505 :   long emaxmu = EX0, emax2mu = EX0;
     774             :   s64 xx;
     775    29811505 :   int did_something = 0;
     776             :   /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
     777             : 
     778    17108124 :   for (;;) {
     779    46919629 :     int go_on = 0;
     780    46919629 :     long i, j, emax3mu = emax2mu;
     781             : 
     782    46919629 :     if (gc_needed(av,2))
     783             :     {
     784          86 :       if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
     785          86 :       gc_lll(av,2,&B,&U);
     786             :     }
     787             :     /* Step2: compute the GSO for stage kappa */
     788    46918707 :     emax2mu = emaxmu; emaxmu = EX0;
     789   178530976 :     for (j=aa; j<kappa; j++)
     790             :     {
     791   131612954 :       double g = dmael(G,kappa,j);
     792   567154069 :       for (k = zeros+1; k < j; k++) g -= dmael(mu,j,k) * dmael(r,kappa,k);
     793   131612954 :       dmael(r,kappa,j) = g;
     794   131612954 :       dmael(mu,kappa,j) = dmael(r,kappa,j) / dmael(r,j,j);
     795   131612954 :       emaxmu = maxss(emaxmu, expoB[kappa]-expoB[j]);
     796             :     }
     797             :     /* maxmu doesn't decrease fast enough */
     798    46918022 :     if (emax3mu != EX0 && emax3mu <= emax2mu + 5) {*pB = B; *pU = U; return 1;}
     799             : 
     800   165676122 :     for (j=kappa-1; j>zeros; j--)
     801             :     {
     802   135870578 :       double tmp = fabs(ldexp (dmael(mu,kappa,j), expoB[kappa]-expoB[j]));
     803   135870578 :       if (tmp>eta) { go_on = 1; break; }
     804             :     }
     805             : 
     806             :     /* Step3--5: compute the X_j's  */
     807    46913768 :     if (go_on)
     808    76897584 :       for (j=kappa-1; j>zeros; j--)
     809             :       { /* The code below seemingly handles U = NULL, but in this case d = 0 */
     810    59789898 :         int e = expoB[j] - expoB[kappa];
     811    59789898 :         double tmp = ldexp(dmael(mu,kappa,j), -e), atmp = fabs(tmp);
     812             :         /* tmp = Inf is allowed */
     813    59789898 :         if (atmp <= .5) continue; /* size-reduced */
     814    33711566 :         if (gc_needed(av,2))
     815             :         {
     816         241 :           if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
     817         241 :           gc_lll(av,2,&B,&U);
     818             :         }
     819    33713797 :         did_something = 1;
     820             :         /* we consider separately the case |X| = 1 */
     821    33713797 :         if (atmp <= 1.5)
     822             :         {
     823    22401321 :           if (dmael(mu,kappa,j) > 0) { /* in this case, X = 1 */
     824    46165972 :             for (k=zeros+1; k<j; k++)
     825    34734502 :               dmael(mu,kappa,k) -= ldexp(dmael(mu,j,k), e);
     826   155486154 :             for (i=1; i<=n; i++)
     827   144056079 :               gmael(B,kappa,i) = subii(gmael(B,kappa,i), gmael(B,j,i));
     828   102427325 :             for (i=1; i<=d; i++)
     829    90997112 :               gmael(U,kappa,i) = subii(gmael(U,kappa,i), gmael(U,j,i));
     830             :           } else { /* otherwise X = -1 */
     831    45335419 :             for (k=zeros+1; k<j; k++)
     832    34365568 :               dmael(mu,kappa,k) += ldexp(dmael(mu,j,k), e);
     833   152764799 :             for (i=1; i<=n; i++)
     834   141796708 :               gmael(B,kappa,i) = addii(gmael(B,kappa,i), gmael(B,j,i));
     835    99818935 :             for (i=1; i<=d; i++)
     836    88850763 :               gmael(U,kappa,i) = addii(gmael(U,kappa,i), gmael(U,j,i));
     837             :           }
     838    22398385 :           continue;
     839             :         }
     840             :         /* we have |X| >= 2 */
     841    11312476 :         if (atmp < 9007199254740992.)
     842             :         {
     843    10460250 :           tmp = pari_rint(tmp);
     844    24716156 :           for (k=zeros+1; k<j; k++)
     845    14255894 :             dmael(mu,kappa,k) -= ldexp(tmp * dmael(mu,j,k), e);
     846    10460262 :           xx = (s64) tmp;
     847    10460262 :           if (xx > 0) /* = xx */
     848             :           {
     849    46144169 :             for (i=1; i<=n; i++)
     850    40891050 :               gmael(B,kappa,i) = submuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
     851    33840987 :             for (i=1; i<=d; i++)
     852    28587764 :               gmael(U,kappa,i) = submuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
     853             :           }
     854             :           else /* = -xx */
     855             :           {
     856    45892769 :             for (i=1; i<=n; i++)
     857    40685627 :               gmael(B,kappa,i) = addmuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), -xx);
     858    33427339 :             for (i=1; i<=d; i++)
     859    28220126 :               gmael(U,kappa,i) = addmuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), -xx);
     860             :           }
     861             :         }
     862             :         else
     863             :         {
     864             :           int E;
     865      852226 :           xx = (s64) ldexp(frexp(dmael(mu,kappa,j), &E), 53);
     866      852226 :           E -= e + 53;
     867      852226 :           if (E <= 0)
     868             :           {
     869           0 :             xx = xx << -E;
     870           0 :             for (k=zeros+1; k<j; k++)
     871           0 :               dmael(mu,kappa,k) -= ldexp(((double)xx) * dmael(mu,j,k), e);
     872           0 :             if (xx > 0) /* = xx */
     873             :             {
     874           0 :               for (i=1; i<=n; i++)
     875           0 :                 gmael(B,kappa,i) = submuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
     876           0 :               for (i=1; i<=d; i++)
     877           0 :                 gmael(U,kappa,i) = submuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
     878             :             }
     879             :             else /* = -xx */
     880             :             {
     881           0 :               for (i=1; i<=n; i++)
     882           0 :                 gmael(B,kappa,i) = addmuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), -xx);
     883           0 :               for (i=1; i<=d; i++)
     884           0 :                 gmael(U,kappa,i) = addmuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), -xx);
     885             :             }
     886             :           } else
     887             :           {
     888     2775638 :             for (k=zeros+1; k<j; k++)
     889     1923412 :               dmael(mu,kappa,k) -= ldexp(((double)xx) * dmael(mu,j,k), E + e);
     890      852226 :             if (xx > 0) /* = xx */
     891             :             {
     892     3938160 :               for (i=1; i<=n; i++)
     893     3508497 :                 gmael(B,kappa,i) = submuliu642n(gmael(B,kappa,i), gmael(B,j,i), xx, E);
     894     1570336 :               for (i=1; i<=d; i++)
     895     1140673 :                 gmael(U,kappa,i) = submuliu642n(gmael(U,kappa,i), gmael(U,j,i), xx, E);
     896             :             }
     897             :             else /* = -xx */
     898             :             {
     899     3870701 :               for (i=1; i<=n; i++)
     900     3448160 :                 gmael(B,kappa,i) = addmuliu642n(gmael(B,kappa,i), gmael(B,j,i), -xx, E);
     901     1525384 :               for (i=1; i<=d; i++)
     902     1102843 :                 gmael(U,kappa,i) = addmuliu642n(gmael(U,kappa,i), gmael(U,j,i), -xx, E);
     903             :             }
     904             :           }
     905             :         }
     906             :       }
     907    46913227 :     if (!go_on) break; /* Anything happened? */
     908    17106247 :     expoB[kappa] = set_line(del(appB,kappa), gel(B,kappa), n);
     909    17107724 :     setG_fast(appB, n, G, kappa, zeros+1, kappa-1);
     910    17108124 :     aa = zeros+1;
     911             :   }
     912    29806980 :   if (did_something) setG2_fast(appB, n, G, kappa, kappa, maxG);
     913             : 
     914    29807276 :   del(s,zeros+1) = dmael(G,kappa,kappa);
     915             :   /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
     916   107785227 :   for (k=zeros+1; k<=kappa-2; k++)
     917    77977951 :     del(s,k+1) = del(s,k) - dmael(mu,kappa,k)*dmael(r,kappa,k);
     918    29807276 :   *pB = B; *pU = U; return 0;
     919             : }
     920             : 
     921             : static void
     922    11982546 : update_alpha(GEN alpha, long kappa, long kappa2, long kappamax)
     923             : {
     924             :   long i;
     925    37735070 :   for (i = kappa; i < kappa2; i++)
     926    25752524 :     if (kappa <= alpha[i]) alpha[i] = kappa;
     927    37735141 :   for (i = kappa2; i > kappa; i--) alpha[i] = alpha[i-1];
     928    22979203 :   for (i = kappa2+1; i <= kappamax; i++)
     929    10996657 :     if (kappa < alpha[i]) alpha[i] = kappa;
     930    11982546 :   alpha[kappa] = kappa;
     931    11982546 : }
     932             : static void
     933      447130 : rotateG(GEN G, long kappa2, long kappa, long maxG, GEN Gtmp)
     934             : {
     935             :   long i, j;
     936     3569746 :   for (i=1; i<=kappa2; i++) gel(Gtmp,i) = gmael(G,kappa2,i);
     937     1874375 :   for (   ; i<=maxG; i++)   gel(Gtmp,i) = gmael(G,i,kappa2);
     938     1565462 :   for (i=kappa2; i>kappa; i--)
     939             :     {
     940     5648596 :       for (j=1; j<kappa; j++) gmael(G,i,j) = gmael(G,i-1,j);
     941     1118332 :       gmael(G,i,kappa) = gel(Gtmp,i-1);
     942     4082101 :       for (j=kappa+1; j<=i; j++) gmael(G,i,j) = gmael(G,i-1,j-1);
     943     4920603 :       for (j=kappa2+1; j<=maxG; j++) gmael(G,j,i) = gmael(G,j,i-1);
     944             :     }
     945     2004284 :   for (i=1; i<kappa; i++) gmael(G,kappa,i) = gel(Gtmp,i);
     946      447130 :   gmael(G,kappa,kappa) = gel(Gtmp,kappa2);
     947     1874375 :   for (i=kappa2+1; i<=maxG; i++) gmael(G,i,kappa) = gel(Gtmp,i);
     948      447130 : }
     949             : static void
     950    11535339 : rotateG_fast(double **G, long kappa2, long kappa, long maxG, double *Gtmp)
     951             : {
     952             :   long i, j;
     953    66654683 :   for (i=1; i<=kappa2; i++) del(Gtmp,i) = dmael(G,kappa2,i);
     954    21978617 :   for (   ; i<=maxG; i++) del(Gtmp,i) = dmael(G,i,kappa2);
     955    36169465 :   for (i=kappa2; i>kappa; i--)
     956             :   {
     957    69194179 :     for (j=1; j<kappa; j++) dmael(G,i,j) = dmael(G,i-1,j);
     958    24634126 :     dmael(G,i,kappa) = del(Gtmp,i-1);
     959    83346855 :     for (j=kappa+1; j<=i; j++) dmael(G,i,j) = dmael(G,i-1,j-1);
     960    46152047 :     for (j=kappa2+1; j<=maxG; j++) dmael(G,j,i) = dmael(G,j,i-1);
     961             :   }
     962    30485526 :   for (i=1; i<kappa; i++) dmael(G,kappa,i) = del(Gtmp,i);
     963    11535339 :   dmael(G,kappa,kappa) = del(Gtmp,kappa2);
     964    21978699 :   for (i=kappa2+1; i<=maxG; i++) dmael(G,i,kappa) = del(Gtmp,i);
     965    11535339 : }
     966             : 
     967             : /* LLL-reduces (B,U) in place [apply base change transforms to B and U].
     968             :  * Gram matrix, and GSO performed on matrices of 'double'.
     969             :  * If (keepfirst), never swap with first vector.
     970             :  * Return -1 on failure, else zeros = dim Kernel (>= 0) */
     971             : static long
     972     2080148 : fplll_fast(GEN *pB, GEN *pU, double delta, double eta, long keepfirst)
     973             : {
     974             :   pari_sp av;
     975             :   long kappa, kappa2, d, n, i, j, zeros, kappamax, maxG;
     976             :   double **mu, **r, *s, tmp, *Gtmp, **G, **appB;
     977     2080148 :   GEN alpha, expoB, B = *pB, U;
     978     2080148 :   long cnt = 0;
     979             : 
     980     2080148 :   d = lg(B)-1;
     981     2080148 :   n = nbrows(B);
     982     2080149 :   U = *pU; /* NULL if inplace */
     983             : 
     984     2080149 :   G = cget_dblmat(d+1);
     985     2080149 :   appB = cget_dblmat(d+1);
     986     2080149 :   mu = cget_dblmat(d+1);
     987     2080148 :   r  = cget_dblmat(d+1);
     988     2080148 :   s  = cget_dblvec(d+1);
     989     9713165 :   for (j = 1; j <= d; j++)
     990             :   {
     991     7633017 :     del(mu,j) = cget_dblvec(d+1);
     992     7633007 :     del(r,j) = cget_dblvec(d+1);
     993     7633002 :     del(appB,j) = cget_dblvec(n+1);
     994     7632998 :     del(G,j) = cget_dblvec(d+1);
     995    47685608 :     for (i=1; i<=d; i++) dmael(G,j,i) = 0.;
     996             :   }
     997     2080148 :   expoB = cgetg(d+1, t_VECSMALL);
     998     9712996 :   for (i=1; i<=d; i++) expoB[i] = set_line(del(appB,i), gel(B,i), n);
     999     2080078 :   Gtmp = cget_dblvec(d+1);
    1000     2080130 :   alpha = cgetg(d+1, t_VECSMALL);
    1001     2080129 :   av = avma;
    1002             : 
    1003             :   /* Step2: Initializing the main loop */
    1004     2080129 :   kappamax = 1;
    1005     2080129 :   i = 1;
    1006     2080129 :   maxG = d; /* later updated to kappamax */
    1007             : 
    1008             :   do {
    1009     2236613 :     dmael(G,i,i) = dbldotsquare(del(appB,i),n);
    1010     2236617 :   } while (dmael(G,i,i) <= 0 && (++i <=d));
    1011     2080133 :   zeros = i-1; /* all vectors B[i] with i <= zeros are zero vectors */
    1012     2080133 :   kappa = i;
    1013     2080133 :   if (zeros < d) dmael(r,zeros+1,zeros+1) = dmael(G,zeros+1,zeros+1);
    1014     9556642 :   for (i=zeros+1; i<=d; i++) alpha[i]=1;
    1015    31887632 :   while (++kappa <= d)
    1016             :   {
    1017    29811735 :     if (kappa > kappamax)
    1018             :     {
    1019     5392265 :       if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
    1020     5392265 :       maxG = kappamax = kappa;
    1021     5392265 :       setG_fast(appB, n, G, kappa, zeros+1, kappa);
    1022             :     }
    1023             :     /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
    1024    29811728 :     if (Babai_fast(av, kappa, &B,&U, mu,r,s, appB, expoB, G, alpha[kappa],
    1025        4254 :                    zeros, maxG, eta)) { *pB=B; *pU=U; return -1; }
    1026             : 
    1027    29807242 :     tmp = ldexp(r[kappa-1][kappa-1] * delta, 2*(expoB[kappa-1]-expoB[kappa]));
    1028    29807242 :     if ((keepfirst && kappa == 2) || tmp <= del(s,kappa-1))
    1029             :     { /* Step4: Success of Lovasz's condition */
    1030    18271828 :       alpha[kappa] = kappa;
    1031    18271828 :       tmp = dmael(mu,kappa,kappa-1) * dmael(r,kappa,kappa-1);
    1032    18271828 :       dmael(r,kappa,kappa) = del(s,kappa-1)- tmp;
    1033    18271828 :       continue;
    1034             :     }
    1035             :     /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
    1036    11535414 :     if (DEBUGLEVEL>=4 && kappa==kappamax && del(s,kappa-1)!=0)
    1037           0 :       if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", 2*expoB[1] + dblexpo(del(s,1))); }
    1038    11535409 :     kappa2 = kappa;
    1039             :     do {
    1040    24634184 :       kappa--;
    1041    24634184 :       if (kappa<zeros+2 + (keepfirst ? 1: 0)) break;
    1042    17955188 :       tmp = dmael(r,kappa-1,kappa-1) * delta;
    1043    17955188 :       tmp = ldexp(tmp, 2*(expoB[kappa-1]-expoB[kappa2]));
    1044    17955188 :     } while (del(s,kappa-1) <= tmp);
    1045    11535409 :     update_alpha(alpha, kappa, kappa2, kappamax);
    1046             : 
    1047             :     /* Step6: Update the mu's and r's */
    1048    11535433 :     dblrotate(mu,kappa2,kappa);
    1049    11535414 :     dblrotate(r,kappa2,kappa);
    1050    11535383 :     dmael(r,kappa,kappa) = del(s,kappa);
    1051             : 
    1052             :     /* Step7: Update B, appB, U, G */
    1053    11535383 :     rotate(B,kappa2,kappa);
    1054    11535377 :     dblrotate(appB,kappa2,kappa);
    1055    11535354 :     if (U) rotate(U,kappa2,kappa);
    1056    11535354 :     rotate(expoB,kappa2,kappa);
    1057    11535339 :     rotateG_fast(G,kappa2,kappa, maxG, Gtmp);
    1058             : 
    1059             :     /* Step8: Prepare the next loop iteration */
    1060    11535671 :     if (kappa == zeros+1 && dmael(G,kappa,kappa)<= 0)
    1061             :     {
    1062      212645 :       zeros++; kappa++;
    1063      212645 :       dmael(G,kappa,kappa) = dbldotsquare(del(appB,kappa),n);
    1064      212645 :       dmael(r,kappa,kappa) = dmael(G,kappa,kappa);
    1065             :     }
    1066             :   }
    1067     2075897 :   *pB = B; *pU = U; return zeros;
    1068             : }
    1069             : 
    1070             : /***************** HEURISTIC version (reduced precision) ****************/
    1071             : static GEN
    1072      186486 : realsqrdotproduct(GEN x)
    1073             : {
    1074      186486 :   long i, l = lg(x);
    1075      186486 :   GEN z = sqrr(gel(x,1));
    1076     1474993 :   for (i=2; i<l; i++) z = addrr(z, sqrr(gel(x,i)));
    1077      186486 :   return z;
    1078             : }
    1079             : /* x, y non-empty vector of t_REALs, same length */
    1080             : static GEN
    1081     1397760 : realdotproduct(GEN x, GEN y)
    1082             : {
    1083             :   long i, l;
    1084             :   GEN z;
    1085     1397760 :   if (x == y) return realsqrdotproduct(x);
    1086     1211274 :   l = lg(x); z = mulrr(gel(x,1),gel(y,1));
    1087    14242421 :   for (i=2; i<l; i++) z = addrr(z, mulrr(gel(x,i), gel(y,i)));
    1088     1211274 :   return z;
    1089             : }
    1090             : static void
    1091      198575 : setG_heuristic(GEN appB, GEN G, long kappa, long a, long b)
    1092      198575 : { pari_sp av = avma;
    1093             :   long i;
    1094     1121435 :   for (i = a; i <= b; i++)
    1095      922860 :     affrr(realdotproduct(gel(appB,kappa),gel(appB,i)), gmael(G,kappa,i));
    1096      198575 :   set_avma(av);
    1097      198575 : }
    1098             : static void
    1099      175753 : setG2_heuristic(GEN appB, GEN G, long kappa, long a, long b)
    1100      175753 : { pari_sp av = avma;
    1101             :   long i;
    1102      650653 :   for (i = a; i <= b; i++)
    1103      474900 :     affrr(realdotproduct(gel(appB,kappa),gel(appB,i)), gmael(G,i,kappa));
    1104      175753 :   set_avma(av);
    1105      175753 : }
    1106             : 
    1107             : /* approximate t_REAL x as m * 2^e, where |m| < 2^bit */
    1108             : static GEN
    1109       32828 : truncexpo(GEN x, long bit, long *e)
    1110             : {
    1111       32828 :   *e = expo(x) + 1 - bit;
    1112       32828 :   if (*e >= 0) return mantissa2nr(x, 0);
    1113        1156 :   *e = 0; return roundr_safe(x);
    1114             : }
    1115             : /* Babai's Nearest Plane algorithm (iterative); see Babai() */
    1116             : static int
    1117      285708 : Babai_heuristic(pari_sp av, long kappa, GEN *pB, GEN *pU, GEN mu, GEN r, GEN s,
    1118             :                 GEN appB, GEN G, long a, long zeros, long maxG,
    1119             :                 GEN eta, long prec)
    1120             : {
    1121      285708 :   GEN B = *pB, U = *pU;
    1122      285708 :   const long n = nbrows(B), d = U ? lg(U)-1: 0, bit = prec2nbits(prec);
    1123      285708 :   long k, aa = (a > zeros)? a : zeros+1;
    1124      285708 :   int did_something = 0;
    1125      285708 :   long emaxmu = EX0, emax2mu = EX0;
    1126             :   /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
    1127             : 
    1128      187842 :   for (;;) {
    1129      473550 :     int go_on = 0;
    1130      473550 :     long i, j, emax3mu = emax2mu;
    1131             : 
    1132      473550 :     if (gc_needed(av,2))
    1133             :     {
    1134          25 :       if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
    1135          25 :       gc_lll(av,2,&B,&U);
    1136             :     }
    1137             :     /* Step2: compute the GSO for stage kappa */
    1138      473550 :     emax2mu = emaxmu; emaxmu = EX0;
    1139     2100814 :     for (j=aa; j<kappa; j++)
    1140             :     {
    1141     1627264 :       pari_sp btop = avma;
    1142     1627264 :       GEN g = gmael(G,kappa,j);
    1143     7831794 :       for (k = zeros+1; k<j; k++)
    1144     6204530 :         g = subrr(g, mulrr(gmael(mu,j,k), gmael(r,kappa,k)));
    1145     1627264 :       affrr(g, gmael(r,kappa,j));
    1146     1627264 :       affrr(divrr(gmael(r,kappa,j), gmael(r,j,j)), gmael(mu,kappa,j));
    1147     1627264 :       emaxmu = maxss(emaxmu, expo(gmael(mu,kappa,j)));
    1148     1627264 :       set_avma(btop);
    1149             :     }
    1150      473550 :     if (emax3mu != EX0 && emax3mu <= emax2mu + 5)
    1151        1334 :     { *pB = B; *pU = U; return 1; }
    1152             : 
    1153     2078167 :     for (j=kappa-1; j>zeros; j--)
    1154     1793793 :       if (abscmprr(gmael(mu,kappa,j), eta) > 0) { go_on = 1; break; }
    1155             : 
    1156             :     /* Step3--5: compute the X_j's  */
    1157      472216 :     if (go_on)
    1158     1071766 :       for (j=kappa-1; j>zeros; j--)
    1159             :       { /* The code below seemingly handles U = NULL, but in this case d = 0 */
    1160             :         pari_sp btop;
    1161      883924 :         GEN tmp = gmael(mu,kappa,j);
    1162      883924 :         if (absrsmall(tmp)) continue; /* size-reduced */
    1163             : 
    1164      469696 :         if (gc_needed(av,2))
    1165             :         {
    1166          34 :           if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
    1167          34 :           gc_lll(av,2,&B,&U);
    1168             :         }
    1169      469696 :         btop = avma; did_something = 1;
    1170             :         /* we consider separately the case |X| = 1 */
    1171      469696 :         if (absrsmall2(tmp))
    1172             :         {
    1173      311356 :           if (signe(tmp) > 0) { /* in this case, X = 1 */
    1174      657803 :             for (k=zeros+1; k<j; k++)
    1175      502615 :               affrr(subrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
    1176      155188 :             set_avma(btop);
    1177     1988579 :             for (i=1; i<=n; i++)
    1178     1833391 :               gmael(B,kappa,i) = subii(gmael(B,kappa,i), gmael(B,j,i));
    1179     1570446 :             for (i=1; i<=d; i++)
    1180     1415258 :               gmael(U,kappa,i) = subii(gmael(U,kappa,i), gmael(U,j,i));
    1181             :           } else { /* otherwise X = -1 */
    1182      665547 :             for (k=zeros+1; k<j; k++)
    1183      509379 :               affrr(addrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
    1184      156168 :             set_avma(btop);
    1185     2005481 :             for (i=1; i<=n; i++)
    1186     1849313 :               gmael(B,kappa,i) = addii(gmael(B,kappa,i), gmael(B,j,i));
    1187     1574693 :             for (i=1; i<=d; i++)
    1188     1418525 :               gmael(U,kappa,i) = addii(gmael(U,kappa,i),gmael(U,j,i));
    1189             :           }
    1190      311356 :           continue;
    1191             :         }
    1192             :         /* we have |X| >= 2 */
    1193      158340 :         if (expo(tmp) < BITS_IN_LONG)
    1194             :         {
    1195      125512 :           ulong xx = roundr_safe(tmp)[2]; /* X fits in an ulong */
    1196      125512 :           if (signe(tmp) > 0) /* = xx */
    1197             :           {
    1198      152389 :             for (k=zeros+1; k<j; k++)
    1199       89272 :               affrr(subrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
    1200       89272 :                   gmael(mu,kappa,k));
    1201       63117 :             set_avma(btop);
    1202      463660 :             for (i=1; i<=n; i++)
    1203      400543 :               gmael(B,kappa,i) = submuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
    1204      359024 :             for (i=1; i<=d; i++)
    1205      295907 :               gmael(U,kappa,i) = submuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
    1206             :           }
    1207             :           else /* = -xx */
    1208             :           {
    1209      150016 :             for (k=zeros+1; k<j; k++)
    1210       87621 :               affrr(addrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
    1211       87621 :                   gmael(mu,kappa,k));
    1212       62395 :             set_avma(btop);
    1213      460046 :             for (i=1; i<=n; i++)
    1214      397651 :               gmael(B,kappa,i) = addmuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
    1215      349356 :             for (i=1; i<=d; i++)
    1216      286961 :               gmael(U,kappa,i) = addmuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
    1217             :           }
    1218             :         }
    1219             :         else
    1220             :         {
    1221             :           long e;
    1222       32828 :           GEN X = truncexpo(tmp, bit, &e); /* tmp ~ X * 2^e */
    1223       32828 :           btop = avma;
    1224      131412 :           for (k=zeros+1; k<j; k++)
    1225             :           {
    1226       98584 :             GEN x = mulir(X, gmael(mu,j,k));
    1227       98584 :             if (e) shiftr_inplace(x, e);
    1228       98584 :             affrr(subrr(gmael(mu,kappa,k), x), gmael(mu,kappa,k));
    1229             :           }
    1230       32828 :           set_avma(btop);
    1231      405295 :           for (i=1; i<=n; i++)
    1232      372467 :             gmael(B,kappa,i) = submulshift(gmael(B,kappa,i), gmael(B,j,i), X, e);
    1233      381035 :           for (i=1; i<=d; i++)
    1234      348207 :             gmael(U,kappa,i) = submulshift(gmael(U,kappa,i), gmael(U,j,i), X, e);
    1235             :         }
    1236             :       }
    1237      472216 :     if (!go_on) break; /* Anything happened? */
    1238     1681903 :     for (i=1 ; i<=n; i++) affir(gmael(B,kappa,i), gmael(appB,kappa,i));
    1239      187842 :     setG_heuristic(appB, G, kappa, zeros+1, kappa-1);
    1240      187842 :     aa = zeros+1;
    1241             :   }
    1242      284374 :   if (did_something) setG2_heuristic(appB, G, kappa, kappa, maxG);
    1243      284374 :   affrr(gmael(G,kappa,kappa), gel(s,zeros+1));
    1244             :   /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
    1245      284374 :   av = avma;
    1246     1458300 :   for (k=zeros+1; k<=kappa-2; k++)
    1247     1173926 :     affrr(subrr(gel(s,k), mulrr(gmael(mu,kappa,k), gmael(r,kappa,k))),
    1248     1173926 :           gel(s,k+1));
    1249      284374 :   *pB = B; *pU = U; return gc_bool(av, 0);
    1250             : }
    1251             : 
    1252             : static GEN
    1253       16931 : ZC_to_RC(GEN x, long prec)
    1254      120832 : { pari_APPLY_type(t_COL,itor(gel(x,i),prec)) }
    1255             : 
    1256             : static GEN
    1257        4254 : ZM_to_RM(GEN x, long prec)
    1258       21185 : { pari_APPLY_same(ZC_to_RC(gel(x,i),prec)) }
    1259             : 
    1260             : /* LLL-reduces (B,U) in place [apply base change transforms to B and U].
    1261             :  * Gram matrix made of t_REAL at precision prec2, performe GSO at prec.
    1262             :  * If (keepfirst), never swap with first vector.
    1263             :  * Return -1 on failure, else zeros = dim Kernel (>= 0) */
    1264             : static long
    1265        4254 : fplll_heuristic(GEN *pB, GEN *pU, double DELTA, double ETA, long keepfirst,
    1266             :                 long prec, long prec2)
    1267             : {
    1268             :   pari_sp av, av2;
    1269             :   long kappa, kappa2, d, i, j, zeros, kappamax, maxG;
    1270        4254 :   GEN mu, r, s, tmp, Gtmp, alpha, G, appB, B = *pB, U;
    1271        4254 :   GEN delta = dbltor(DELTA), eta = dbltor(ETA);
    1272        4254 :   long cnt = 0;
    1273             : 
    1274        4254 :   d = lg(B)-1;
    1275        4254 :   U = *pU; /* NULL if inplace */
    1276             : 
    1277        4254 :   G = cgetg(d+1, t_MAT);
    1278        4254 :   mu = cgetg(d+1, t_MAT);
    1279        4254 :   r  = cgetg(d+1, t_MAT);
    1280        4254 :   s  = cgetg(d+1, t_VEC);
    1281        4254 :   appB = ZM_to_RM(B, prec2);
    1282       21185 :   for (j = 1; j <= d; j++)
    1283             :   {
    1284       16931 :     GEN M = cgetg(d+1, t_COL), R = cgetg(d+1, t_COL), S = cgetg(d+1, t_COL);
    1285       16931 :     gel(mu,j)= M;
    1286       16931 :     gel(r,j) = R;
    1287       16931 :     gel(G,j) = S;
    1288       16931 :     gel(s,j) = cgetr(prec);
    1289      112508 :     for (i = 1; i <= d; i++)
    1290             :     {
    1291       95577 :       gel(R,i) = cgetr(prec);
    1292       95577 :       gel(M,i) = cgetr(prec);
    1293       95577 :       gel(S,i) = cgetr(prec2);
    1294             :     }
    1295             :   }
    1296        4254 :   Gtmp = cgetg(d+1, t_VEC);
    1297        4254 :   alpha = cgetg(d+1, t_VECSMALL);
    1298        4254 :   av = avma;
    1299             : 
    1300             :   /* Step2: Initializing the main loop */
    1301        4254 :   kappamax = 1;
    1302        4254 :   i = 1;
    1303        4254 :   maxG = d; /* later updated to kappamax */
    1304             : 
    1305             :   do {
    1306        4272 :     affrr(RgV_dotsquare(gel(appB,i)), gmael(G,i,i));
    1307        4272 :   } while (signe(gmael(G,i,i)) == 0 && (++i <=d));
    1308        4254 :   zeros = i-1; /* all vectors B[i] with i <= zeros are zero vectors */
    1309        4254 :   kappa = i;
    1310        4254 :   if (zeros < d) affrr(gmael(G,zeros+1,zeros+1), gmael(r,zeros+1,zeros+1));
    1311       21167 :   for (i=zeros+1; i<=d; i++) alpha[i]=1;
    1312             : 
    1313      288628 :   while (++kappa <= d)
    1314             :   {
    1315      285708 :     if (kappa > kappamax)
    1316             :     {
    1317       10733 :       if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
    1318       10733 :       maxG = kappamax = kappa;
    1319       10733 :       setG_heuristic(appB, G, kappa, zeros+1, kappa);
    1320             :     }
    1321             :     /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
    1322      285708 :     if (Babai_heuristic(av, kappa, &B,&U, mu,r,s, appB, G, alpha[kappa], zeros,
    1323        1334 :                         maxG, eta, prec)) { *pB = B; *pU = U; return -1; }
    1324      284374 :     av2 = avma;
    1325      568662 :     if ((keepfirst && kappa == 2) ||
    1326      284288 :         cmprr(mulrr(gmael(r,kappa-1,kappa-1), delta), gel(s,kappa-1)) <= 0)
    1327             :     { /* Step4: Success of Lovasz's condition */
    1328      170045 :       alpha[kappa] = kappa;
    1329      170045 :       tmp = mulrr(gmael(mu,kappa,kappa-1), gmael(r,kappa,kappa-1));
    1330      170045 :       affrr(subrr(gel(s,kappa-1), tmp), gmael(r,kappa,kappa));
    1331      170045 :       set_avma(av2); continue;
    1332             :     }
    1333             :     /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
    1334      114329 :     if (DEBUGLEVEL>=4 && kappa==kappamax && signe(gel(s,kappa-1)))
    1335           0 :       if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", expo(gel(s,1))); }
    1336      114329 :     kappa2 = kappa;
    1337             :     do {
    1338      274992 :       kappa--;
    1339      274992 :       if (kappa < zeros+2 + (keepfirst ? 1: 0)) break;
    1340      250815 :       tmp = mulrr(gmael(r,kappa-1,kappa-1), delta);
    1341      250815 :     } while (cmprr(gel(s,kappa-1), tmp) <= 0 );
    1342      114329 :     set_avma(av2);
    1343      114329 :     update_alpha(alpha, kappa, kappa2, kappamax);
    1344             : 
    1345             :     /* Step6: Update the mu's and r's */
    1346      114329 :     rotate(mu,kappa2,kappa);
    1347      114329 :     rotate(r,kappa2,kappa);
    1348      114329 :     affrr(gel(s,kappa), gmael(r,kappa,kappa));
    1349             : 
    1350             :     /* Step7: Update B, appB, U, G */
    1351      114329 :     rotate(B,kappa2,kappa);
    1352      114329 :     rotate(appB,kappa2,kappa);
    1353      114329 :     if (U) rotate(U,kappa2,kappa);
    1354      114329 :     rotateG(G,kappa2,kappa, maxG, Gtmp);
    1355             : 
    1356             :     /* Step8: Prepare the next loop iteration */
    1357      114329 :     if (kappa == zeros+1 && !signe(gmael(G,kappa,kappa)))
    1358             :     {
    1359           0 :       zeros++; kappa++;
    1360           0 :       affrr(RgV_dotsquare(gel(appB,kappa)), gmael(G,kappa,kappa));
    1361           0 :       affrr(gmael(G,kappa,kappa), gmael(r,kappa,kappa));
    1362             :     }
    1363             :   }
    1364        2920 :   *pB=B; *pU=U; return zeros;
    1365             : }
    1366             : 
    1367             : /************************* PROVED version (t_INT) ***********************/
    1368             : /* dpe inspired by dpe.h by Patrick Pelissier, Paul Zimmermann
    1369             :  * https://gforge.inria.fr/projects/dpe/
    1370             :  */
    1371             : 
    1372             : typedef struct
    1373             : {
    1374             :   double d;  /* significand */
    1375             :   long e; /* exponent */
    1376             : } dpe_t;
    1377             : 
    1378             : #define Dmael(x,i,j) (&((x)[i][j]))
    1379             : #define Del(x,i) (&((x)[i]))
    1380             : 
    1381             : static void
    1382      665602 : dperotate(dpe_t **A, long k2, long k)
    1383             : {
    1384             :   long i;
    1385      665602 :   dpe_t *B = A[k2];
    1386     2352282 :   for (i = k2; i > k; i--) A[i] = A[i-1];
    1387      665602 :   A[k] = B;
    1388      665602 : }
    1389             : 
    1390             : static void
    1391   108265804 : dpe_normalize0(dpe_t *x)
    1392             : {
    1393             :   int e;
    1394   108265804 :   x->d = frexp(x->d, &e);
    1395   108265804 :   x->e += e;
    1396   108265804 : }
    1397             : 
    1398             : static void
    1399    47949432 : dpe_normalize(dpe_t *x)
    1400             : {
    1401    47949432 :   if (x->d == 0.0)
    1402      500296 :     x->e = -LONG_MAX;
    1403             :   else
    1404    47449136 :     dpe_normalize0(x);
    1405    47949510 : }
    1406             : 
    1407             : static GEN
    1408       93410 : dpetor(dpe_t *x)
    1409             : {
    1410       93410 :   GEN r = dbltor(x->d);
    1411       93410 :   if (signe(r)==0) return r;
    1412       93410 :   setexpo(r, x->e-1);
    1413       93410 :   return r;
    1414             : }
    1415             : 
    1416             : static void
    1417    25749116 : affdpe(dpe_t *y, dpe_t *x)
    1418             : {
    1419    25749116 :   x->d = y->d;
    1420    25749116 :   x->e = y->e;
    1421    25749116 : }
    1422             : 
    1423             : static void
    1424    20660441 : affidpe(GEN y, dpe_t *x)
    1425             : {
    1426    20660441 :   pari_sp av = avma;
    1427    20660441 :   GEN r = itor(y, DEFAULTPREC);
    1428    20660038 :   x->e = expo(r)+1;
    1429    20660038 :   setexpo(r,-1);
    1430    20659949 :   x->d = rtodbl(r);
    1431    20659869 :   set_avma(av);
    1432    20659817 : }
    1433             : 
    1434             : static void
    1435     3150213 : affdbldpe(double y, dpe_t *x)
    1436             : {
    1437     3150213 :   x->d = (double)y;
    1438     3150213 :   x->e = 0;
    1439     3150213 :   dpe_normalize(x);
    1440     3150214 : }
    1441             : 
    1442             : static void
    1443    56860313 : dpe_mulz(dpe_t *x, dpe_t *y, dpe_t *z)
    1444             : {
    1445    56860313 :   z->d = x->d * y->d;
    1446    56860313 :   if (z->d == 0.0)
    1447     8175200 :     z->e = -LONG_MAX;
    1448             :   else
    1449             :   {
    1450    48685113 :     z->e = x->e + y->e;
    1451    48685113 :     dpe_normalize0(z);
    1452             :   }
    1453    56860688 : }
    1454             : 
    1455             : static void
    1456    14053027 : dpe_divz(dpe_t *x, dpe_t *y, dpe_t *z)
    1457             : {
    1458    14053027 :   z->d = x->d / y->d;
    1459    14053027 :   if (z->d == 0.0)
    1460     1921010 :     z->e = -LONG_MAX;
    1461             :   else
    1462             :   {
    1463    12132017 :     z->e = x->e - y->e;
    1464    12132017 :     dpe_normalize0(z);
    1465             :   }
    1466    14053105 : }
    1467             : 
    1468             : static void
    1469      242927 : dpe_negz(dpe_t *y, dpe_t *x)
    1470             : {
    1471      242927 :   x->d = - y->d;
    1472      242927 :   x->e = y->e;
    1473      242927 : }
    1474             : 
    1475             : static void
    1476     1952819 : dpe_addz(dpe_t *y, dpe_t *z, dpe_t *x)
    1477             : {
    1478     1952819 :   if (y->e > z->e + 53)
    1479      113382 :     affdpe(y, x);
    1480     1839437 :   else if (z->e > y->e + 53)
    1481       41842 :     affdpe(z, x);
    1482             :   else
    1483             :   {
    1484     1797595 :     long d = y->e - z->e;
    1485             : 
    1486     1797595 :     if (d >= 0)
    1487             :     {
    1488     1349018 :       x->d = y->d + ldexp(z->d, -d);
    1489     1349018 :       x->e  = y->e;
    1490             :     }
    1491             :     else
    1492             :     {
    1493      448577 :       x->d = z->d + ldexp(y->d, d);
    1494      448577 :       x->e = z->e;
    1495             :     }
    1496     1797595 :     dpe_normalize(x);
    1497             :   }
    1498     1952828 : }
    1499             : static void
    1500    53653084 : dpe_subz(dpe_t *y, dpe_t *z, dpe_t *x)
    1501             : {
    1502    53653084 :   if (y->e > z->e + 53)
    1503    11208217 :     affdpe(y, x);
    1504    42444867 :   else if (z->e > y->e + 53)
    1505      242927 :     dpe_negz(z, x);
    1506             :   else
    1507             :   {
    1508    42201940 :     long d = y->e - z->e;
    1509             : 
    1510    42201940 :     if (d >= 0)
    1511             :     {
    1512    39428288 :       x->d = y->d - ldexp(z->d, -d);
    1513    39428288 :       x->e = y->e;
    1514             :     }
    1515             :     else
    1516             :     {
    1517     2773652 :       x->d = ldexp(y->d, d) - z->d;
    1518     2773652 :       x->e = z->e;
    1519             :     }
    1520    42201940 :     dpe_normalize(x);
    1521             :   }
    1522    53653463 : }
    1523             : 
    1524             : static void
    1525      799435 : dpe_muluz(dpe_t *y, ulong t, dpe_t *x)
    1526             : {
    1527      799435 :   x->d = y->d * (double)t;
    1528      799435 :   x->e = y->e;
    1529      799435 :   dpe_normalize(x);
    1530      799435 : }
    1531             : 
    1532             : static void
    1533      343412 : dpe_addmuluz(dpe_t *y,  dpe_t *z, ulong t, dpe_t *x)
    1534             : {
    1535             :   dpe_t tmp;
    1536      343412 :   dpe_muluz(z, t, &tmp);
    1537      343412 :   dpe_addz(y, &tmp, x);
    1538      343412 : }
    1539             : 
    1540             : static void
    1541      408825 : dpe_submuluz(dpe_t *y,  dpe_t *z, ulong t, dpe_t *x)
    1542             : {
    1543             :   dpe_t tmp;
    1544      408825 :   dpe_muluz(z, t, &tmp);
    1545      408825 :   dpe_subz(y, &tmp, x);
    1546      408825 : }
    1547             : 
    1548             : static void
    1549    51616441 : dpe_submulz(dpe_t *y,  dpe_t *z, dpe_t *t, dpe_t *x)
    1550             : {
    1551             :   dpe_t tmp;
    1552    51616441 :   dpe_mulz(z, t, &tmp);
    1553    51616513 :   dpe_subz(y, &tmp, x);
    1554    51616613 : }
    1555             : 
    1556             : static int
    1557     5244020 : dpe_cmp(dpe_t *x, dpe_t *y)
    1558             : {
    1559     5244020 :   int sx = x->d < 0. ? -1: x->d > 0.;
    1560     5244020 :   int sy = y->d < 0. ? -1: y->d > 0.;
    1561     5244020 :   int d  = sx - sy;
    1562             : 
    1563     5244020 :   if (d != 0)
    1564      141657 :     return d;
    1565     5102363 :   else if (x->e > y->e)
    1566      488818 :     return (sx > 0) ? 1 : -1;
    1567     4613545 :   else if (y->e > x->e)
    1568     2522597 :     return (sx > 0) ? -1 : 1;
    1569             :   else
    1570     2090948 :     return (x->d < y->d) ? -1 : (x->d > y->d);
    1571             : }
    1572             : 
    1573             : static int
    1574    14489456 : dpe_abscmp(dpe_t *x, dpe_t *y)
    1575             : {
    1576    14489456 :   if (x->e > y->e)
    1577      277204 :     return 1;
    1578    14212252 :   else if (y->e > x->e)
    1579    13357712 :     return -1;
    1580             :   else
    1581      854540 :     return (fabs(x->d) < fabs(y->d)) ? -1 : (fabs(x->d) > fabs(y->d));
    1582             : }
    1583             : 
    1584             : static int
    1585     1406343 : dpe_abssmall(dpe_t *x)
    1586             : {
    1587     1406343 :   return (x->e <= 0) || (x->e == 1 && fabs(x->d) <= .75);
    1588             : }
    1589             : 
    1590             : static int
    1591     5244023 : dpe_cmpmul(dpe_t *x, dpe_t *y, dpe_t *z)
    1592             : {
    1593             :   dpe_t t;
    1594     5244023 :   dpe_mulz(x,y,&t);
    1595     5244020 :   return dpe_cmp(&t, z);
    1596             : }
    1597             : 
    1598             : static dpe_t *
    1599    13103747 : cget_dpevec(long d)
    1600    13103747 : { return (dpe_t*) stack_malloc_align(d*sizeof(dpe_t), sizeof(dpe_t)); }
    1601             : 
    1602             : static dpe_t **
    1603     3150208 : cget_dpemat(long d) { return (dpe_t **) cgetg(d, t_VECSMALL); }
    1604             : 
    1605             : static GEN
    1606       20028 : dpeM_diagonal_shallow(dpe_t **m, long d)
    1607             : {
    1608             :   long i;
    1609       20028 :   GEN y = cgetg(d+1,t_VEC);
    1610      113438 :   for (i=1; i<=d; i++) gel(y, i) = dpetor(Dmael(m,i,i));
    1611       20028 :   return y;
    1612             : }
    1613             : 
    1614             : static void
    1615     1406319 : affii_or_copy_gc(pari_sp av, GEN x, GEN *y)
    1616             : {
    1617     1406319 :   long l = lg(*y);
    1618     1406319 :   if (lgefint(x) <= l && isonstack(*y))
    1619             :   {
    1620     1405884 :     affii(x,*y);
    1621     1405883 :     set_avma(av);
    1622             :   }
    1623             :   else
    1624         433 :     *y = gerepileuptoint(av, x);
    1625     1406322 : }
    1626             : 
    1627             : /* *x -= u*y */
    1628             : INLINE void
    1629     5907963 : submulziu(GEN *x, GEN y, ulong u)
    1630             : {
    1631             :   pari_sp av;
    1632     5907963 :   long ly = lgefint(y);
    1633     5907963 :   if (ly == 2) return;
    1634     3253634 :   av = avma;
    1635     3253634 :   (void)new_chunk(3+ly+lgefint(*x)); /* HACK */
    1636     3253678 :   y = mului(u,y);
    1637     3253661 :   set_avma(av); subzi(x, y);
    1638             : }
    1639             : 
    1640             : /* *x += u*y */
    1641             : INLINE void
    1642     4599844 : addmulziu(GEN *x, GEN y, ulong u)
    1643             : {
    1644             :   pari_sp av;
    1645     4599844 :   long ly = lgefint(y);
    1646     4599844 :   if (ly == 2) return;
    1647     2776629 :   av = avma;
    1648     2776629 :   (void)new_chunk(3+ly+lgefint(*x)); /* HACK */
    1649     2776675 :   y = mului(u,y);
    1650     2776650 :   set_avma(av); addzi(x, y);
    1651             : }
    1652             : 
    1653             : /************************** PROVED version (dpe) *************************/
    1654             : 
    1655             : /* Babai's Nearest Plane algorithm (iterative).
    1656             :  * Size-reduces b_kappa using mu_{i,j} and r_{i,j} for j<=i <kappa
    1657             :  * Update B[,kappa]; compute mu_{kappa,j}, r_{kappa,j} for j<=kappa and s[kappa]
    1658             :  * mu, r, s updated in place (affrr). Return 1 on failure, else 0. */
    1659             : static int
    1660     4622960 : Babai_dpe(pari_sp av, long kappa, GEN *pG, GEN *pB, GEN *pU, dpe_t **mu, dpe_t **r, dpe_t *s,
    1661             :       long a, long zeros, long maxG, dpe_t *eta)
    1662             : {
    1663     4622960 :   GEN G = *pG, B = *pB, U = *pU, ztmp;
    1664     4622960 :   long k, d, n, aa = a > zeros? a: zeros+1;
    1665     4622960 :   long emaxmu = EX0, emax2mu = EX0;
    1666             :   /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
    1667     4622960 :   d = U? lg(U)-1: 0;
    1668     4622960 :   n = B? nbrows(B): 0;
    1669      532745 :   for (;;) {
    1670     5155728 :     int go_on = 0;
    1671     5155728 :     long i, j, emax3mu = emax2mu;
    1672             : 
    1673     5155728 :     if (gc_needed(av,2))
    1674             :     {
    1675           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
    1676           0 :       gc_lll(av,3,&G,&B,&U);
    1677             :     }
    1678             :     /* Step2: compute the GSO for stage kappa */
    1679     5155704 :     emax2mu = emaxmu; emaxmu = EX0;
    1680    19208801 :     for (j=aa; j<kappa; j++)
    1681             :     {
    1682             :       dpe_t g;
    1683    14053085 :       affidpe(gmael(G,kappa,j), &g);
    1684    52464351 :       for (k = zeros+1; k < j; k++)
    1685    38411333 :         dpe_submulz(&g, Dmael(mu,j,k), Dmael(r,kappa,k), &g);
    1686    14053018 :       affdpe(&g, Dmael(r,kappa,j));
    1687    14053019 :       dpe_divz(Dmael(r,kappa,j), Dmael(r,j,j), Dmael(mu,kappa,j));
    1688    14053040 :       emaxmu = maxss(emaxmu, Dmael(mu,kappa,j)->e);
    1689             :     }
    1690     5155716 :     if (emax3mu != EX0 && emax3mu <= emax2mu + 5) /* precision too low */
    1691           0 :     { *pG = G; *pB = B; *pU = U; return 1; }
    1692             : 
    1693    19112434 :     for (j=kappa-1; j>zeros; j--)
    1694    14489454 :       if (dpe_abscmp(Dmael(mu,kappa,j), eta) > 0) { go_on = 1; break; }
    1695             : 
    1696             :     /* Step3--5: compute the X_j's  */
    1697     5155717 :     if (go_on)
    1698     3079275 :       for (j=kappa-1; j>zeros; j--)
    1699             :       {
    1700             :         pari_sp btop;
    1701     2546527 :         dpe_t *tmp = Dmael(mu,kappa,j);
    1702     2546527 :         if (tmp->e < 0) continue; /* (essentially) size-reduced */
    1703             : 
    1704     1406341 :         if (gc_needed(av,2))
    1705             :         {
    1706           0 :           if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
    1707           0 :           gc_lll(av,3,&G,&B,&U);
    1708             :         }
    1709             :         /* we consider separately the case |X| = 1 */
    1710     1406341 :         if (dpe_abssmall(tmp))
    1711             :         {
    1712      931414 :           if (tmp->d > 0) { /* in this case, X = 1 */
    1713     2070998 :             for (k=zeros+1; k<j; k++)
    1714     1604384 :               dpe_subz(Dmael(mu,kappa,k), Dmael(mu,j,k), Dmael(mu,kappa,k));
    1715     3042196 :             for (i=1; i<=n; i++)
    1716     2575582 :               subzi(&gmael(B,kappa,i), gmael(B,j,i));
    1717     7001996 :             for (i=1; i<=d; i++)
    1718     6535391 :               subzi(&gmael(U,kappa,i), gmael(U,j,i));
    1719      466605 :             btop = avma;
    1720      466605 :             ztmp = subii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
    1721      466606 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    1722      466612 :             affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
    1723     2877989 :             for (i=1; i<=j; i++)
    1724     2411381 :               subzi(&gmael(G,kappa,i), gmael(G,j,i));
    1725     2202583 :             for (i=j+1; i<kappa; i++)
    1726     1735974 :               subzi(&gmael(G,kappa,i), gmael(G,i,j));
    1727     2370021 :             for (i=kappa+1; i<=maxG; i++)
    1728     1903411 :               subzi(&gmael(G,i,kappa), gmael(G,i,j));
    1729             :           } else { /* otherwise X = -1 */
    1730     2050436 :             for (k=zeros+1; k<j; k++)
    1731     1585630 :               dpe_addz(Dmael(mu,kappa,k), Dmael(mu,j,k), Dmael(mu,kappa,k));
    1732     3039210 :             for (i=1; i<=n; i++)
    1733     2574404 :               addzi(&gmael(B,kappa,i),gmael(B,j,i));
    1734     6881179 :             for (i=1; i<=d; i++)
    1735     6416386 :               addzi(&gmael(U,kappa,i),gmael(U,j,i));
    1736      464793 :             btop = avma;
    1737      464793 :             ztmp = addii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
    1738      464800 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    1739      464800 :             affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
    1740     2790836 :             for (i=1; i<=j; i++)
    1741     2326037 :               addzi(&gmael(G,kappa,i), gmael(G,j,i));
    1742     2209193 :             for (i=j+1; i<kappa; i++)
    1743     1744394 :               addzi(&gmael(G,kappa,i), gmael(G,i,j));
    1744     2324465 :             for (i=kappa+1; i<=maxG; i++)
    1745     1859663 :               addzi(&gmael(G,i,kappa), gmael(G,i,j));
    1746             :           }
    1747      931412 :           continue;
    1748             :         }
    1749             :         /* we have |X| >= 2 */
    1750      474931 :         if (tmp->e < BITS_IN_LONG-1)
    1751             :         {
    1752      450772 :           if (tmp->d > 0)
    1753             :           {
    1754      247570 :             ulong xx = (ulong) pari_rint(ldexp(tmp->d, tmp->e)); /* X fits in an ulong */
    1755      656396 :             for (k=zeros+1; k<j; k++)
    1756      408825 :               dpe_submuluz(Dmael(mu,kappa,k), Dmael(mu,j,k), xx, Dmael(mu,kappa,k));
    1757      734442 :             for (i=1; i<=n; i++)
    1758      486871 :               submulziu(&gmael(B,kappa,i), gmael(B,j,i), xx);
    1759     3096045 :             for (i=1; i<=d; i++)
    1760     2848476 :               submulziu(&gmael(U,kappa,i), gmael(U,j,i), xx);
    1761      247569 :             btop = avma;
    1762      247569 :             ztmp = submuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
    1763      247571 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    1764      247571 :             affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
    1765     1311519 :             for (i=1; i<=j; i++)
    1766     1063951 :               submulziu(&gmael(G,kappa,i), gmael(G,j,i), xx);
    1767      804784 :             for (i=j+1; i<kappa; i++)
    1768      557217 :               submulziu(&gmael(G,kappa,i), gmael(G,i,j), xx);
    1769     1199087 :             for (i=kappa+1; i<=maxG; i++)
    1770      951520 :               submulziu(&gmael(G,i,kappa), gmael(G,i,j), xx);
    1771             :           }
    1772             :           else
    1773             :           {
    1774      203202 :             ulong xx = (ulong) pari_rint(ldexp(-tmp->d, tmp->e)); /* X fits in an ulong */
    1775      546617 :             for (k=zeros+1; k<j; k++)
    1776      343412 :               dpe_addmuluz(Dmael(mu,kappa,k), Dmael(mu,j,k), xx, Dmael(mu,kappa,k));
    1777      703496 :             for (i=1; i<=n; i++)
    1778      500291 :               addmulziu(&gmael(B,kappa,i), gmael(B,j,i), xx);
    1779     2367409 :             for (i=1; i<=d; i++)
    1780     2164206 :               addmulziu(&gmael(U,kappa,i), gmael(U,j,i), xx);
    1781      203203 :             btop = avma;
    1782      203203 :             ztmp = addmuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
    1783      203204 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    1784      203202 :             affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
    1785      995321 :             for (i=1; i<=j; i++)
    1786      792117 :               addmulziu(&gmael(G,kappa,i), gmael(G,j,i), xx);
    1787      664541 :             for (i=j+1; i<kappa; i++)
    1788      461338 :               addmulziu(&gmael(G,kappa,i), gmael(G,i,j), xx);
    1789      885150 :             for (i=kappa+1; i<=maxG; i++)
    1790      681947 :               addmulziu(&gmael(G,i,kappa), gmael(G,i,j), xx);
    1791             :           }
    1792             :         }
    1793             :         else
    1794             :         {
    1795       24159 :           long e = tmp->e - BITS_IN_LONG + 1;
    1796       24159 :           if (tmp->d > 0)
    1797             :           {
    1798       12035 :             ulong xx = (ulong) pari_rint(ldexp(tmp->d, BITS_IN_LONG - 1));
    1799       35455 :             for (k=zeros+1; k<j; k++)
    1800             :             {
    1801             :               dpe_t x;
    1802       23420 :               dpe_muluz(Dmael(mu,j,k), xx, &x);
    1803       23420 :               x.e += e;
    1804       23420 :               dpe_subz(Dmael(mu,kappa,k), &x, Dmael(mu,kappa,k));
    1805             :             }
    1806      148590 :             for (i=1; i<=n; i++)
    1807      136555 :               submulzu2n(&gmael(B,kappa,i), gmael(B,j,i), xx, e);
    1808      111267 :             for (i=1; i<=d; i++)
    1809       99232 :               submulzu2n(&gmael(U,kappa,i), gmael(U,j,i), xx, e);
    1810       12035 :             btop = avma;
    1811       12035 :             ztmp = submuliu2n(mulshift(gmael(G,j,j), sqru(xx), 2*e),
    1812       12035 :                 gmael(G,kappa,j), xx, e+1);
    1813       12035 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    1814       12035 :             affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
    1815       48096 :             for (i=1; i<=j; i++)
    1816       36061 :               submulzu2n(&gmael(G,kappa,i), gmael(G,j,i), xx, e);
    1817       58095 :             for (   ; i<kappa; i++)
    1818       46060 :               submulzu2n(&gmael(G,kappa,i), gmael(G,i,j), xx, e);
    1819       14149 :             for (i=kappa+1; i<=maxG; i++)
    1820        2114 :               submulzu2n(&gmael(G,i,kappa), gmael(G,i,j), xx, e);
    1821             :           } else
    1822             :           {
    1823       12124 :             ulong xx = (ulong) pari_rint(ldexp(-tmp->d, BITS_IN_LONG - 1));
    1824       35914 :             for (k=zeros+1; k<j; k++)
    1825             :             {
    1826             :               dpe_t x;
    1827       23779 :               dpe_muluz(Dmael(mu,j,k), xx, &x);
    1828       23779 :               x.e += e;
    1829       23779 :               dpe_addz(Dmael(mu,kappa,k), &x, Dmael(mu,kappa,k));
    1830             :             }
    1831      152427 :             for (i=1; i<=n; i++)
    1832      140292 :               addmulzu2n(&gmael(B,kappa,i), gmael(B,j,i), xx, e);
    1833      113043 :             for (i=1; i<=d; i++)
    1834      100908 :               addmulzu2n(&gmael(U,kappa,i), gmael(U,j,i), xx, e);
    1835       12135 :             btop = avma;
    1836       12135 :             ztmp = addmuliu2n(mulshift(gmael(G,j,j), sqru(xx), 2*e),
    1837       12135 :                 gmael(G,kappa,j), xx, e+1);
    1838       12135 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    1839       12135 :             affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
    1840       48823 :             for (i=1; i<=j; i++)
    1841       36688 :               addmulzu2n(&gmael(G,kappa,i), gmael(G,j,i), xx, e);
    1842       59402 :             for (   ; i<kappa; i++)
    1843       47267 :               addmulzu2n(&gmael(G,kappa,i), gmael(G,i,j), xx, e);
    1844       13920 :             for (i=kappa+1; i<=maxG; i++)
    1845        1785 :               addmulzu2n(&gmael(G,i,kappa), gmael(G,i,j), xx, e);
    1846             :           }
    1847             :         }
    1848             :       }
    1849     5155728 :     if (!go_on) break; /* Anything happened? */
    1850      532745 :     aa = zeros+1;
    1851             :   }
    1852             : 
    1853     4622983 :   affidpe(gmael(G,kappa,kappa), Del(s,zeros+1));
    1854             :   /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
    1855    13538237 :   for (k=zeros+1; k<=kappa-2; k++)
    1856     8915265 :     dpe_submulz(Del(s,k), Dmael(mu,kappa,k), Dmael(r,kappa,k), Del(s,k+1));
    1857     4622972 :   *pG = G; *pB = B; *pU = U; return 0;
    1858             : }
    1859             : 
    1860             : /* G integral Gram matrix, LLL-reduces (G,B,U) in place [apply base change
    1861             :  * transforms to B and U]. If (keepfirst), never swap with first vector.
    1862             :  * If G = NULL, we compute the Gram matrix incrementally.
    1863             :  * Return -1 on failure, else zeros = dim Kernel (>= 0) */
    1864             : static long
    1865     1575107 : fplll_dpe(GEN *pG, GEN *pB, GEN *pU, GEN *pr, double DELTA, double ETA,
    1866             :       long keepfirst)
    1867             : {
    1868             :   pari_sp av;
    1869     1575107 :   GEN Gtmp, alpha, G = *pG, B = *pB, U = *pU;
    1870     1575107 :   long d, maxG, kappa, kappa2, i, j, zeros, kappamax, incgram = !G, cnt = 0;
    1871             :   dpe_t delta, eta, **mu, **r, *s;
    1872     1575107 :   affdbldpe(DELTA,&delta);
    1873     1575107 :   affdbldpe(ETA,&eta);
    1874             : 
    1875     1575108 :   if (incgram)
    1876             :   { /* incremental Gram matrix */
    1877     1514639 :     maxG = 2; d = lg(B)-1;
    1878     1514639 :     G = zeromatcopy(d, d);
    1879             :   }
    1880             :   else
    1881       60469 :     maxG = d = lg(G)-1;
    1882             : 
    1883     1575107 :   mu = cget_dpemat(d+1);
    1884     1575106 :   r  = cget_dpemat(d+1);
    1885     1575101 :   s  = cget_dpevec(d+1);
    1886     7339459 :   for (j = 1; j <= d; j++)
    1887             :   {
    1888     5764354 :     mu[j]= cget_dpevec(d+1);
    1889     5764356 :     r[j] = cget_dpevec(d+1);
    1890             :   }
    1891     1575105 :   Gtmp = cgetg(d+1, t_VEC);
    1892     1575100 :   alpha = cgetg(d+1, t_VECSMALL);
    1893     1575093 :   av = avma;
    1894             : 
    1895             :   /* Step2: Initializing the main loop */
    1896     1575093 :   kappamax = 1;
    1897     1575093 :   i = 1;
    1898             :   do {
    1899     1949504 :     if (incgram) gmael(G,i,i) = ZV_dotsquare(gel(B,i));
    1900     1949451 :     affidpe(gmael(G,i,i), Dmael(r,i,i));
    1901     1949460 :   } while (!signe(gmael(G,i,i)) && ++i <= d);
    1902     1575049 :   zeros = i-1; /* all basis vectors b_i with i <= zeros are zero vectors */
    1903     1575049 :   kappa = i;
    1904     6964864 :   for (i=zeros+1; i<=d; i++) alpha[i]=1;
    1905             : 
    1906     6198051 :   while (++kappa <= d)
    1907             :   {
    1908     4622958 :     if (kappa > kappamax)
    1909             :     {
    1910     3814779 :       if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
    1911     3814779 :       kappamax = kappa;
    1912     3814779 :       if (incgram)
    1913             :       {
    1914    15973316 :         for (i=zeros+1; i<=kappa; i++)
    1915    12360331 :           gmael(G,kappa,i) = ZV_dotproduct(gel(B,kappa), gel(B,i));
    1916     3612985 :         maxG = kappamax;
    1917             :       }
    1918             :     }
    1919             :     /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
    1920     4622704 :     if (Babai_dpe(av, kappa, &G,&B,&U, mu,r,s, alpha[kappa], zeros, maxG, &eta))
    1921           0 :     { *pG = incgram? NULL: G; *pB = B; *pU = U; return -1; }
    1922     9147945 :     if ((keepfirst && kappa == 2) ||
    1923     4524959 :         dpe_cmpmul(Dmael(r,kappa-1,kappa-1), &delta, Del(s,kappa-1)) <= 0)
    1924             :     { /* Step4: Success of Lovasz's condition */
    1925     4290185 :       alpha[kappa] = kappa;
    1926     4290185 :       dpe_submulz(Del(s,kappa-1), Dmael(mu,kappa,kappa-1), Dmael(r,kappa,kappa-1), Dmael(r,kappa,kappa));
    1927     4290195 :       continue;
    1928             :     }
    1929             :     /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
    1930      332801 :     if (DEBUGLEVEL>=4 && kappa==kappamax && Del(s,kappa-1)->d)
    1931           0 :       if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", Del(s,1)->e-1); }
    1932      332801 :     kappa2 = kappa;
    1933             :     do {
    1934      843340 :       kappa--;
    1935      843340 :       if (kappa < zeros+2 + (keepfirst ? 1: 0)) break;
    1936      719058 :     } while (dpe_cmpmul(Dmael(r,kappa-1,kappa-1), &delta, Del(s,kappa-1)) >= 0);
    1937      332801 :     update_alpha(alpha, kappa, kappa2, kappamax);
    1938             : 
    1939             :     /* Step6: Update the mu's and r's */
    1940      332801 :     dperotate(mu, kappa2, kappa);
    1941      332801 :     dperotate(r, kappa2, kappa);
    1942      332801 :     affdpe(Del(s,kappa), Dmael(r,kappa,kappa));
    1943             : 
    1944             :     /* Step7: Update G, B, U */
    1945      332801 :     if (U) rotate(U, kappa2, kappa);
    1946      332801 :     if (B) rotate(B, kappa2, kappa);
    1947      332801 :     rotateG(G,kappa2,kappa, maxG, Gtmp);
    1948             : 
    1949             :     /* Step8: Prepare the next loop iteration */
    1950      332801 :     if (kappa == zeros+1 && !signe(gmael(G,kappa,kappa)))
    1951             :     {
    1952       35161 :       zeros++; kappa++;
    1953       35161 :       affidpe(gmael(G,kappa,kappa), Dmael(r,kappa,kappa));
    1954             :     }
    1955             :   }
    1956     1575093 :   if (pr) *pr = dpeM_diagonal_shallow(r,d);
    1957     1575093 :   *pG = G; *pB = B; *pU = U; return zeros; /* success */
    1958             : }
    1959             : 
    1960             : 
    1961             : /************************** PROVED version (t_INT) *************************/
    1962             : 
    1963             : /* Babai's Nearest Plane algorithm (iterative).
    1964             :  * Size-reduces b_kappa using mu_{i,j} and r_{i,j} for j<=i <kappa
    1965             :  * Update B[,kappa]; compute mu_{kappa,j}, r_{kappa,j} for j<=kappa and s[kappa]
    1966             :  * mu, r, s updated in place (affrr). Return 1 on failure, else 0. */
    1967             : static int
    1968           0 : Babai(pari_sp av, long kappa, GEN *pG, GEN *pB, GEN *pU, GEN mu, GEN r, GEN s,
    1969             :       long a, long zeros, long maxG, GEN eta, long prec)
    1970             : {
    1971           0 :   GEN G = *pG, B = *pB, U = *pU, ztmp;
    1972           0 :   long k, aa = a > zeros? a: zeros+1;
    1973           0 :   const long n = B? nbrows(B): 0, d = U ? lg(U)-1: 0, bit = prec2nbits(prec);
    1974           0 :   long emaxmu = EX0, emax2mu = EX0;
    1975             :   /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
    1976             : 
    1977           0 :   for (;;) {
    1978           0 :     int go_on = 0;
    1979           0 :     long i, j, emax3mu = emax2mu;
    1980             : 
    1981           0 :     if (gc_needed(av,2))
    1982             :     {
    1983           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
    1984           0 :       gc_lll(av,3,&G,&B,&U);
    1985             :     }
    1986             :     /* Step2: compute the GSO for stage kappa */
    1987           0 :     emax2mu = emaxmu; emaxmu = EX0;
    1988           0 :     for (j=aa; j<kappa; j++)
    1989             :     {
    1990           0 :       pari_sp btop = avma;
    1991           0 :       GEN g = gmael(G,kappa,j);
    1992           0 :       for (k = zeros+1; k < j; k++)
    1993           0 :         g = mpsub(g, mulrr(gmael(mu,j,k), gmael(r,kappa,k)));
    1994           0 :       mpaff(g, gmael(r,kappa,j));
    1995           0 :       affrr(divrr(gmael(r,kappa,j), gmael(r,j,j)), gmael(mu,kappa,j));
    1996           0 :       emaxmu = maxss(emaxmu, expo(gmael(mu,kappa,j)));
    1997           0 :       set_avma(btop);
    1998             :     }
    1999           0 :     if (emax3mu != EX0 && emax3mu <= emax2mu + 5) /* precision too low */
    2000           0 :     { *pG = G; *pB = B; *pU = U; return 1; }
    2001             : 
    2002           0 :     for (j=kappa-1; j>zeros; j--)
    2003           0 :       if (abscmprr(gmael(mu,kappa,j), eta) > 0) { go_on = 1; break; }
    2004             : 
    2005             :     /* Step3--5: compute the X_j's  */
    2006           0 :     if (go_on)
    2007           0 :       for (j=kappa-1; j>zeros; j--)
    2008             :       {
    2009             :         pari_sp btop;
    2010           0 :         GEN tmp = gmael(mu,kappa,j);
    2011           0 :         if (absrsmall(tmp)) continue; /* size-reduced */
    2012             : 
    2013           0 :         if (gc_needed(av,2))
    2014             :         {
    2015           0 :           if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
    2016           0 :           gc_lll(av,3,&G,&B,&U);
    2017             :         }
    2018           0 :         btop = avma;
    2019             :         /* we consider separately the case |X| = 1 */
    2020           0 :         if (absrsmall2(tmp))
    2021             :         {
    2022           0 :           if (signe(tmp) > 0) { /* in this case, X = 1 */
    2023           0 :             for (k=zeros+1; k<j; k++)
    2024           0 :               affrr(subrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
    2025           0 :             set_avma(btop);
    2026           0 :             for (i=1; i<=n; i++)
    2027           0 :               gmael(B,kappa,i) = subii(gmael(B,kappa,i), gmael(B,j,i));
    2028           0 :             for (i=1; i<=d; i++)
    2029           0 :               gmael(U,kappa,i) = subii(gmael(U,kappa,i), gmael(U,j,i));
    2030           0 :             btop = avma;
    2031           0 :             ztmp = subii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
    2032           0 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    2033           0 :             gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
    2034           0 :             for (i=1; i<=j; i++)
    2035           0 :               gmael(G,kappa,i) = subii(gmael(G,kappa,i), gmael(G,j,i));
    2036           0 :             for (i=j+1; i<kappa; i++)
    2037           0 :               gmael(G,kappa,i) = subii(gmael(G,kappa,i), gmael(G,i,j));
    2038           0 :             for (i=kappa+1; i<=maxG; i++)
    2039           0 :               gmael(G,i,kappa) = subii(gmael(G,i,kappa), gmael(G,i,j));
    2040             :           } else { /* otherwise X = -1 */
    2041           0 :             for (k=zeros+1; k<j; k++)
    2042           0 :               affrr(addrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
    2043           0 :             set_avma(btop);
    2044           0 :             for (i=1; i<=n; i++)
    2045           0 :               gmael(B,kappa,i) = addii(gmael(B,kappa,i),gmael(B,j,i));
    2046           0 :             for (i=1; i<=d; i++)
    2047           0 :               gmael(U,kappa,i) = addii(gmael(U,kappa,i),gmael(U,j,i));
    2048           0 :             btop = avma;
    2049           0 :             ztmp = addii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
    2050           0 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    2051           0 :             gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
    2052           0 :             for (i=1; i<=j; i++)
    2053           0 :               gmael(G,kappa,i) = addii(gmael(G,kappa,i), gmael(G,j,i));
    2054           0 :             for (i=j+1; i<kappa; i++)
    2055           0 :               gmael(G,kappa,i) = addii(gmael(G,kappa,i), gmael(G,i,j));
    2056           0 :             for (i=kappa+1; i<=maxG; i++)
    2057           0 :               gmael(G,i,kappa) = addii(gmael(G,i,kappa), gmael(G,i,j));
    2058             :           }
    2059           0 :           continue;
    2060             :         }
    2061             :         /* we have |X| >= 2 */
    2062           0 :         if (expo(tmp) < BITS_IN_LONG)
    2063             :         {
    2064           0 :           ulong xx = roundr_safe(tmp)[2]; /* X fits in an ulong */
    2065           0 :           if (signe(tmp) > 0) /* = xx */
    2066             :           {
    2067           0 :             for (k=zeros+1; k<j; k++)
    2068           0 :               affrr(subrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
    2069           0 :                   gmael(mu,kappa,k));
    2070           0 :             set_avma(btop);
    2071           0 :             for (i=1; i<=n; i++)
    2072           0 :               gmael(B,kappa,i) = submuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
    2073           0 :             for (i=1; i<=d; i++)
    2074           0 :               gmael(U,kappa,i) = submuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
    2075           0 :             btop = avma;
    2076           0 :             ztmp = submuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
    2077           0 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    2078           0 :             gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
    2079           0 :             for (i=1; i<=j; i++)
    2080           0 :               gmael(G,kappa,i) = submuliu_inplace(gmael(G,kappa,i), gmael(G,j,i), xx);
    2081           0 :             for (i=j+1; i<kappa; i++)
    2082           0 :               gmael(G,kappa,i) = submuliu_inplace(gmael(G,kappa,i), gmael(G,i,j), xx);
    2083           0 :             for (i=kappa+1; i<=maxG; i++)
    2084           0 :               gmael(G,i,kappa) = submuliu_inplace(gmael(G,i,kappa), gmael(G,i,j), xx);
    2085             :           }
    2086             :           else /* = -xx */
    2087             :           {
    2088           0 :             for (k=zeros+1; k<j; k++)
    2089           0 :               affrr(addrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
    2090           0 :                   gmael(mu,kappa,k));
    2091           0 :             set_avma(btop);
    2092           0 :             for (i=1; i<=n; i++)
    2093           0 :               gmael(B,kappa,i) = addmuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
    2094           0 :             for (i=1; i<=d; i++)
    2095           0 :               gmael(U,kappa,i) = addmuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
    2096           0 :             btop = avma;
    2097           0 :             ztmp = addmuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
    2098           0 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    2099           0 :             gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
    2100           0 :             for (i=1; i<=j; i++)
    2101           0 :               gmael(G,kappa,i) = addmuliu_inplace(gmael(G,kappa,i), gmael(G,j,i), xx);
    2102           0 :             for (i=j+1; i<kappa; i++)
    2103           0 :               gmael(G,kappa,i) = addmuliu_inplace(gmael(G,kappa,i), gmael(G,i,j), xx);
    2104           0 :             for (i=kappa+1; i<=maxG; i++)
    2105           0 :               gmael(G,i,kappa) = addmuliu_inplace(gmael(G,i,kappa), gmael(G,i,j), xx);
    2106             :           }
    2107             :         }
    2108             :         else
    2109             :         {
    2110             :           long e;
    2111           0 :           GEN X = truncexpo(tmp, bit, &e); /* tmp ~ X * 2^e */
    2112           0 :           btop = avma;
    2113           0 :           for (k=zeros+1; k<j; k++)
    2114             :           {
    2115           0 :             GEN x = mulir(X, gmael(mu,j,k));
    2116           0 :             if (e) shiftr_inplace(x, e);
    2117           0 :             affrr(subrr(gmael(mu,kappa,k), x), gmael(mu,kappa,k));
    2118             :           }
    2119           0 :           set_avma(btop);
    2120           0 :           for (i=1; i<=n; i++)
    2121           0 :             gmael(B,kappa,i) = submulshift(gmael(B,kappa,i), gmael(B,j,i), X, e);
    2122           0 :           for (i=1; i<=d; i++)
    2123           0 :             gmael(U,kappa,i) = submulshift(gmael(U,kappa,i), gmael(U,j,i), X, e);
    2124           0 :           btop = avma;
    2125           0 :           ztmp = submulshift(mulshift(gmael(G,j,j), sqri(X), 2*e),
    2126           0 :               gmael(G,kappa,j), X, e+1);
    2127           0 :           ztmp = addii(gmael(G,kappa,kappa), ztmp);
    2128           0 :           gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
    2129           0 :           for (i=1; i<=j; i++)
    2130           0 :             gmael(G,kappa,i) = submulshift(gmael(G,kappa,i), gmael(G,j,i), X, e);
    2131           0 :           for (   ; i<kappa; i++)
    2132           0 :             gmael(G,kappa,i) = submulshift(gmael(G,kappa,i), gmael(G,i,j), X, e);
    2133           0 :           for (i=kappa+1; i<=maxG; i++)
    2134           0 :             gmael(G,i,kappa) = submulshift(gmael(G,i,kappa), gmael(G,i,j), X, e);
    2135             :         }
    2136             :       }
    2137           0 :     if (!go_on) break; /* Anything happened? */
    2138           0 :     aa = zeros+1;
    2139             :   }
    2140             : 
    2141           0 :   affir(gmael(G,kappa,kappa), gel(s,zeros+1));
    2142             :   /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
    2143           0 :   av = avma;
    2144           0 :   for (k=zeros+1; k<=kappa-2; k++)
    2145           0 :     affrr(subrr(gel(s,k), mulrr(gmael(mu,kappa,k), gmael(r,kappa,k))),
    2146           0 :           gel(s,k+1));
    2147           0 :   *pG = G; *pB = B; *pU = U; return gc_bool(av, 0);
    2148             : }
    2149             : 
    2150             : /* G integral Gram matrix, LLL-reduces (G,B,U) in place [apply base change
    2151             :  * transforms to B and U]. If (keepfirst), never swap with first vector.
    2152             :  * If G = NULL, we compute the Gram matrix incrementally.
    2153             :  * Return -1 on failure, else zeros = dim Kernel (>= 0) */
    2154             : static long
    2155           0 : fplll(GEN *pG, GEN *pB, GEN *pU, GEN *pr, double DELTA, double ETA,
    2156             :       long keepfirst, long prec)
    2157             : {
    2158             :   pari_sp av, av2;
    2159           0 :   GEN mu, r, s, tmp, Gtmp, alpha, G = *pG, B = *pB, U = *pU;
    2160           0 :   GEN delta = dbltor(DELTA), eta = dbltor(ETA);
    2161           0 :   long d, maxG, kappa, kappa2, i, j, zeros, kappamax, incgram = !G, cnt = 0;
    2162             : 
    2163           0 :   if (incgram)
    2164             :   { /* incremental Gram matrix */
    2165           0 :     maxG = 2; d = lg(B)-1;
    2166           0 :     G = zeromatcopy(d, d);
    2167             :   }
    2168             :   else
    2169           0 :     maxG = d = lg(G)-1;
    2170             : 
    2171           0 :   mu = cgetg(d+1, t_MAT);
    2172           0 :   r  = cgetg(d+1, t_MAT);
    2173           0 :   s  = cgetg(d+1, t_VEC);
    2174           0 :   for (j = 1; j <= d; j++)
    2175             :   {
    2176           0 :     GEN M = cgetg(d+1, t_COL), R = cgetg(d+1, t_COL);
    2177           0 :     gel(mu,j)= M;
    2178           0 :     gel(r,j) = R;
    2179           0 :     gel(s,j) = cgetr(prec);
    2180           0 :     for (i = 1; i <= d; i++)
    2181             :     {
    2182           0 :       gel(R,i) = cgetr(prec);
    2183           0 :       gel(M,i) = cgetr(prec);
    2184             :     }
    2185             :   }
    2186           0 :   Gtmp = cgetg(d+1, t_VEC);
    2187           0 :   alpha = cgetg(d+1, t_VECSMALL);
    2188           0 :   av = avma;
    2189             : 
    2190             :   /* Step2: Initializing the main loop */
    2191           0 :   kappamax = 1;
    2192           0 :   i = 1;
    2193             :   do {
    2194           0 :     if (incgram) gmael(G,i,i) = ZV_dotsquare(gel(B,i));
    2195           0 :     affir(gmael(G,i,i), gmael(r,i,i));
    2196           0 :   } while (!signe(gmael(G,i,i)) && ++i <= d);
    2197           0 :   zeros = i-1; /* all basis vectors b_i with i <= zeros are zero vectors */
    2198           0 :   kappa = i;
    2199           0 :   for (i=zeros+1; i<=d; i++) alpha[i]=1;
    2200             : 
    2201           0 :   while (++kappa <= d)
    2202             :   {
    2203           0 :     if (kappa > kappamax)
    2204             :     {
    2205           0 :       if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
    2206           0 :       kappamax = kappa;
    2207           0 :       if (incgram)
    2208             :       {
    2209           0 :         for (i=zeros+1; i<=kappa; i++)
    2210           0 :           gmael(G,kappa,i) = ZV_dotproduct(gel(B,kappa), gel(B,i));
    2211           0 :         maxG = kappamax;
    2212             :       }
    2213             :     }
    2214             :     /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
    2215           0 :     if (Babai(av, kappa, &G,&B,&U, mu,r,s, alpha[kappa], zeros, maxG, eta, prec))
    2216           0 :     { *pG = incgram? NULL: G; *pB = B; *pU = U; return -1; }
    2217           0 :     av2 = avma;
    2218           0 :     if ((keepfirst && kappa == 2) ||
    2219           0 :         cmprr(mulrr(gmael(r,kappa-1,kappa-1), delta), gel(s,kappa-1)) <= 0)
    2220             :     { /* Step4: Success of Lovasz's condition */
    2221           0 :       alpha[kappa] = kappa;
    2222           0 :       tmp = mulrr(gmael(mu,kappa,kappa-1), gmael(r,kappa,kappa-1));
    2223           0 :       affrr(subrr(gel(s,kappa-1), tmp), gmael(r,kappa,kappa));
    2224           0 :       set_avma(av2); continue;
    2225             :     }
    2226             :     /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
    2227           0 :     if (DEBUGLEVEL>=4 && kappa==kappamax && signe(gel(s,kappa-1)))
    2228           0 :       if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", expo(gel(s,1))); }
    2229           0 :     kappa2 = kappa;
    2230             :     do {
    2231           0 :       kappa--;
    2232           0 :       if (kappa < zeros+2 + (keepfirst ? 1: 0)) break;
    2233           0 :       tmp = mulrr(gmael(r,kappa-1,kappa-1), delta);
    2234           0 :     } while (cmprr(gel(s,kappa-1), tmp) <= 0);
    2235           0 :     set_avma(av2);
    2236           0 :     update_alpha(alpha, kappa, kappa2, kappamax);
    2237             : 
    2238             :     /* Step6: Update the mu's and r's */
    2239           0 :     rotate(mu, kappa2, kappa);
    2240           0 :     rotate(r, kappa2, kappa);
    2241           0 :     affrr(gel(s,kappa), gmael(r,kappa,kappa));
    2242             : 
    2243             :     /* Step7: Update G, B, U */
    2244           0 :     if (U) rotate(U, kappa2, kappa);
    2245           0 :     if (B) rotate(B, kappa2, kappa);
    2246           0 :     rotateG(G,kappa2,kappa, maxG, Gtmp);
    2247             : 
    2248             :     /* Step8: Prepare the next loop iteration */
    2249           0 :     if (kappa == zeros+1 && !signe(gmael(G,kappa,kappa)))
    2250             :     {
    2251           0 :       zeros++; kappa++;
    2252           0 :       affir(gmael(G,kappa,kappa), gmael(r,kappa,kappa));
    2253             :     }
    2254             :   }
    2255           0 :   if (pr) *pr = RgM_diagonal_shallow(r);
    2256           0 :   *pG = G; *pB = B; *pU = U; return zeros; /* success */
    2257             : }
    2258             : 
    2259             : /* do not support LLL_KER, LLL_ALL, LLL_KEEP_FIRST */
    2260             : static GEN
    2261     4685173 : ZM2_lll_norms(GEN x, long flag, GEN *pN)
    2262             : {
    2263             :   GEN a,b,c,d;
    2264             :   GEN G, U;
    2265     4685173 :   if (flag & LLL_GRAM)
    2266       15373 :     G = x;
    2267             :   else
    2268     4669800 :     G = gram_matrix(x);
    2269     4685161 :   a = gcoeff(G,1,1); b = shifti(gcoeff(G,1,2),1); c = gcoeff(G,2,2);
    2270     4685133 :   d = qfb_disc3(a,b,c);
    2271     4685120 :   if (signe(d)>=0) return NULL;
    2272     4684735 :   G = redimagsl2(mkqfb(a,b,c,d),&U);
    2273     4684790 :   if (pN) (void) RgM_gram_schmidt(G, pN);
    2274     4684790 :   if (flag & LLL_INPLACE) return ZM2_mul(x,U);
    2275     4684790 :   return U;
    2276             : }
    2277             : 
    2278             : static void
    2279      627351 : fplll_flatter(GEN *pG, GEN *pB, GEN *pU, long rank, long flag)
    2280             : {
    2281      627351 :   if (!*pG)
    2282             :   {
    2283      626029 :     GEN T = ZM_flatter_rank(*pB, rank, flag);
    2284      626029 :     if (*pU)
    2285             :     {
    2286      611731 :       *pU = ZM_mul(*pU, T);
    2287      611731 :       *pB = ZM_mul(*pB, T);
    2288       14298 :     } else *pB = T;
    2289             :   } else
    2290             :   {
    2291        1322 :     GEN T, G = *pG;
    2292        1322 :     long i, j, l = lg(G);
    2293        8716 :     for (i = 1; i < l; i++)
    2294       45817 :       for(j = 1; j < i; j++)
    2295       38423 :         gmael(G,j,i) = gmael(G,i,j);
    2296        1322 :     T = ZM_flattergram_rank(G, rank, flag);
    2297        1322 :     if (*pU) *pU = ZM_mul(*pU, T);
    2298        1322 :     *pG = ZM_mul(shallowtrans(T), ZM_mul(*pG,T));
    2299             :   }
    2300      627351 : }
    2301             : 
    2302             : static long
    2303     1127472 : spread(GEN R)
    2304             : {
    2305     1127472 :   long i, n = lg(R)-1;
    2306     1127472 :   GEN m = gcoeff(R, 1, 1), M = m;
    2307     4813391 :   for (i = 2; i <= n; ++i)
    2308             :   {
    2309     3685917 :     if (mpabscmp(gcoeff(R, i, i), m) < 0)
    2310     1291067 :       m = gcoeff(R, i, i);
    2311     3685918 :     if (mpabscmp(gcoeff(R, i, i), M) > 0)
    2312      857925 :       M = gcoeff(R, i, i);
    2313             :   }
    2314     1127474 :   return mpexpo(M)-mpexpo(m);
    2315             : }
    2316             : 
    2317             : static GEN
    2318     1093934 : get_gramschmidt(GEN M, long rank)
    2319             : {
    2320             :   GEN B, Q, L;
    2321     1093934 :   long l = lg(M), bitprec = 3*(l-1) + 30;
    2322     1093934 :   long prec = nbits2prec64(bitprec);
    2323     1093934 :   if (rank < l-1) M = vconcat(gshift(M,1), matid(l-1));
    2324     1093934 :   if (!QR_init(RgM_gtofp(M, prec), &B, &Q, &L, prec) || !gsisinv(L)) return NULL;
    2325      469084 :   return L;
    2326             : }
    2327             : 
    2328             : static GEN
    2329       49569 : get_gaussred(GEN M, long rank)
    2330             : {
    2331       49569 :   pari_sp ltop = avma;
    2332       49569 :   long lM = lg(M);
    2333       49569 :   long bitprec = 3*(lM-1) + 30, prec = nbits2prec64(bitprec);
    2334             :   GEN R;
    2335       49569 :   if (rank < lM-1) M = RgM_Rg_add(gshift(M, 1), gen_1);
    2336       49569 :   R = RgM_Cholesky(RgM_gtofp(M, prec), prec);
    2337       49570 :   if (!R) return NULL;
    2338       48248 :   return gerepilecopy(ltop, R);
    2339             : }
    2340             : 
    2341             : /* Assume x a ZM, if pN != NULL, set it to Gram-Schmidt (squared) norms
    2342             :  * The following modes are supported:
    2343             :  * - flag & LLL_INPLACE: x a lattice basis, return x*U
    2344             :  * - flag & LLL_GRAM: x a Gram matrix / else x a lattice basis; return
    2345             :  *     LLL base change matrix U [LLL_IM]
    2346             :  *     kernel basis [LLL_KER, nonreduced]
    2347             :  *     both [LLL_ALL] */
    2348             : GEN
    2349     6949323 : ZM_lll_norms(GEN x, double DELTA, long flag, GEN *pN)
    2350             : {
    2351     6949323 :   pari_sp av = avma;
    2352     6949323 :   const double ETA = 0.51;
    2353     6949323 :   const long keepfirst = flag & LLL_KEEP_FIRST;
    2354     6949323 :   long p, zeros = -1, n = lg(x)-1, is_upper, is_lower, useflatter = 0, rank;
    2355     6949323 :   GEN G, B, U, L = NULL;
    2356             :   pari_timer T;
    2357     6949323 :   long thre[]={31783,34393,20894,22525,13533,1928,672,671,
    2358             :                 422,506,315,313,222,205,167,154,139,138,
    2359             :                 110,120,98,94,81,75,74,64,74,74,
    2360             :                 79,96,112,111,105,104,96,86,84,78,75,70,66,62,62,57,56,47,45,52,50,44,48,42,36,35,35,34,40,33,34,32,36,31,
    2361             :                 38,38,40,38,38,37,35,31,34,36,34,32,34,32,28,27,25,31,25,27,28,26,25,21,21,25,25,22,21,24,24,22,21,23,22,22,22,22,21,24,21,22,19,20,19,20,19,19,19,18,19,18,18,20,19,20,18,19,18,21,18,20,18,18};
    2362     6949323 :   long thsn[]={23280,30486,50077,44136,78724,15690,1801,1611,
    2363             :                981,1359,978,1042,815,866,788,775,726,712,
    2364             :                626,613,548,564,474,481,504,447,453,508,
    2365             :                705,794,1008,946,767,898,886,763,842,757,
    2366             :                725,774,639,655,705,627,635,704,511,613,
    2367             :                583,595,568,640,541,640,567,540,577,584,
    2368             :                546,509,526,572,637,746,772,743,743,742,800,708,832,768,707,692,692,768,696,635,709,694,768,719,655,569,590,644,685,623,627,720,633,636,602,635,575,631,642,647,632,656,573,511,688,640,528,616,511,559,601,620,635,688,608,768,658,582,644,704,555,673,600,601,641,661,601,670};
    2369     6949323 :   if (n <= 1) return lll_trivial(x, flag);
    2370     6840451 :   if (nbrows(x)==0)
    2371             :   {
    2372       15190 :     if (flag & LLL_KER) return matid(n);
    2373       15190 :     if (flag & (LLL_INPLACE|LLL_IM)) return cgetg(1,t_MAT);
    2374           0 :     retmkvec2(matid(n), cgetg(1,t_MAT));
    2375             :   }
    2376     6825412 :   if (n==2 && nbrows(x)==2  && (flag&LLL_IM) && !keepfirst)
    2377             :   {
    2378     4685175 :     U = ZM2_lll_norms(x, flag, pN);
    2379     4685175 :     if (U) return U;
    2380             :   }
    2381     2140619 :   rank = ZM_rank(x);
    2382     2140610 :   if (flag & LLL_GRAM)
    2383       60461 :   { G = x; B = NULL; U = matid(n); is_upper = 0; is_lower = 0; }
    2384             :   else
    2385             :   {
    2386     2080149 :     G = NULL; B = x; U = (flag & LLL_INPLACE)? NULL: matid(n);
    2387     2080152 :     is_upper = (flag & LLL_UPPER) || ZM_is_upper(B);
    2388     2080153 :     is_lower = !B || is_upper || keepfirst ? 0: ZM_is_lower(B);
    2389     2080153 :     if (is_lower) L = RgM_flip(B);
    2390             :   }
    2391     2140621 :   if(DEBUGLEVEL>=4) timer_start(&T);
    2392     2140621 :   if (n > 2)
    2393             :   {
    2394     1753642 :     GEN R = B ? is_upper ? B : is_lower ? L : get_gramschmidt(B, rank) : get_gaussred(G, rank);
    2395     1753637 :     if (R)
    2396             :     {
    2397     1127472 :       long spr = spread(R), sz = mpexpo(gsupnorm(R, DEFAULTPREC)), thr;
    2398     1127475 :       if (DEBUGLEVEL>=5) err_printf("LLL: dim %ld, size %ld, spread %ld\n",n, sz, spr);
    2399     1127475 :       if ((is_upper && ZM_is_knapsack(B)) || (is_lower && ZM_is_knapsack(L)))
    2400       95146 :         thr = thsn[minss(n-3,numberof(thsn)-1)];
    2401             :       else
    2402             :       {
    2403     1032329 :         thr = thre[minss(n-3,numberof(thre)-1)];
    2404     1032329 :         if (n>=10) sz = spr;
    2405             :       }
    2406     1127475 :       useflatter = sz >= thr;
    2407             :     } else
    2408      626165 :       useflatter = 1;
    2409      386979 :   } else useflatter = 0;
    2410     2140619 :   if (useflatter)
    2411             :   {
    2412      627351 :     if (is_lower)
    2413             :     {
    2414           0 :       fplll_flatter(&G, &L, &U, rank, flag | LLL_UPPER);
    2415           0 :       B = RgM_flop(L);
    2416           0 :       if (U) U = RgM_flop(U);
    2417             :     }
    2418             :     else
    2419      627351 :       fplll_flatter(&G, &B, &U, rank, flag | (is_upper? LLL_UPPER:0));
    2420      627351 :     if (DEBUGLEVEL>=4  && !(flag & LLL_NOCERTIFY))
    2421           0 :       timer_printf(&T, "FLATTER");
    2422             :   }
    2423     2140619 :   if (!(flag & LLL_GRAM))
    2424             :   {
    2425             :     long t;
    2426     2080150 :     B = gcopy(B);
    2427     2080147 :     if(DEBUGLEVEL>=4)
    2428           0 :       err_printf("Entering L^2 (double): dim %ld, LLL-parameters (%.3f,%.3f)\n",
    2429             :                  n, DELTA,ETA);
    2430     2080147 :     zeros = fplll_fast(&B, &U, DELTA, ETA, keepfirst);
    2431     2080149 :     if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
    2432     2084403 :     for (p = DEFAULTPREC, t = 0; zeros < 0 && t < 1 ; p += EXTRAPREC64, t++)
    2433             :     {
    2434        4254 :       if (DEBUGLEVEL>=4)
    2435           0 :         err_printf("Entering L^2 (heuristic): LLL-parameters (%.3f,%.3f), prec = %d/%d\n", DELTA, ETA, p, p);
    2436        4254 :       zeros = fplll_heuristic(&B, &U, DELTA, ETA, keepfirst, p, p);
    2437        4254 :       gc_lll(av, 2, &B, &U);
    2438        4254 :       if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
    2439             :     }
    2440             :   } else
    2441       60469 :     G = gcopy(G);
    2442     2140618 :   if (zeros < 0 || !(flag & LLL_NOCERTIFY))
    2443             :   {
    2444     1575107 :     if(DEBUGLEVEL>=4)
    2445           0 :       err_printf("Entering L^2 (dpe): LLL-parameters (%.3f,%.3f)\n", DELTA,ETA);
    2446     1575107 :     zeros = fplll_dpe(&G, &B, &U, pN, DELTA, ETA, keepfirst);
    2447     1575092 :     if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
    2448     1575099 :     if (zeros < 0)
    2449           0 :       for (p = DEFAULTPREC;; p += EXTRAPREC64)
    2450             :       {
    2451           0 :         if (DEBUGLEVEL>=4)
    2452           0 :           err_printf("Entering L^2: LLL-parameters (%.3f,%.3f), prec = %d\n",
    2453             :               DELTA,ETA, p);
    2454           0 :         zeros = fplll(&G, &B, &U, pN, DELTA, ETA, keepfirst, p);
    2455           0 :         if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
    2456           0 :         if (zeros >= 0) break;
    2457           0 :         gc_lll(av, 3, &G, &B, &U);
    2458             :       }
    2459             :   }
    2460     2140610 :   return lll_finish(U? U: B, zeros, flag);
    2461             : }
    2462             : 
    2463             : /********************************************************************/
    2464             : /**                                                                **/
    2465             : /**                        LLL OVER K[X]                           **/
    2466             : /**                                                                **/
    2467             : /********************************************************************/
    2468             : static int
    2469         504 : pslg(GEN x)
    2470             : {
    2471             :   long tx;
    2472         504 :   if (gequal0(x)) return 2;
    2473         448 :   tx = typ(x); return is_scalar_t(tx)? 3: lg(x);
    2474             : }
    2475             : 
    2476             : static int
    2477         196 : REDgen(long k, long l, GEN h, GEN L, GEN B)
    2478             : {
    2479         196 :   GEN q, u = gcoeff(L,k,l);
    2480             :   long i;
    2481             : 
    2482         196 :   if (pslg(u) < pslg(B)) return 0;
    2483             : 
    2484         140 :   q = gneg(gdeuc(u,B));
    2485         140 :   gel(h,k) = gadd(gel(h,k), gmul(q,gel(h,l)));
    2486         140 :   for (i=1; i<l; i++) gcoeff(L,k,i) = gadd(gcoeff(L,k,i), gmul(q,gcoeff(L,l,i)));
    2487         140 :   gcoeff(L,k,l) = gadd(gcoeff(L,k,l), gmul(q,B)); return 1;
    2488             : }
    2489             : 
    2490             : static int
    2491         196 : do_SWAPgen(GEN h, GEN L, GEN B, long k, GEN fl, int *flc)
    2492             : {
    2493             :   GEN p1, la, la2, Bk;
    2494             :   long ps1, ps2, i, j, lx;
    2495             : 
    2496         196 :   if (!fl[k-1]) return 0;
    2497             : 
    2498         140 :   la = gcoeff(L,k,k-1); la2 = gsqr(la);
    2499         140 :   Bk = gel(B,k);
    2500         140 :   if (fl[k])
    2501             :   {
    2502          56 :     GEN q = gadd(la2, gmul(gel(B,k-1),gel(B,k+1)));
    2503          56 :     ps1 = pslg(gsqr(Bk));
    2504          56 :     ps2 = pslg(q);
    2505          56 :     if (ps1 <= ps2 && (ps1 < ps2 || !*flc)) return 0;
    2506          28 :     *flc = (ps1 != ps2);
    2507          28 :     gel(B,k) = gdiv(q, Bk);
    2508             :   }
    2509             : 
    2510         112 :   swap(gel(h,k-1), gel(h,k)); lx = lg(L);
    2511         112 :   for (j=1; j<k-1; j++) swap(gcoeff(L,k-1,j), gcoeff(L,k,j));
    2512         112 :   if (fl[k])
    2513             :   {
    2514          28 :     for (i=k+1; i<lx; i++)
    2515             :     {
    2516           0 :       GEN t = gcoeff(L,i,k);
    2517           0 :       p1 = gsub(gmul(gel(B,k+1),gcoeff(L,i,k-1)), gmul(la,t));
    2518           0 :       gcoeff(L,i,k) = gdiv(p1, Bk);
    2519           0 :       p1 = gadd(gmul(la,gcoeff(L,i,k-1)), gmul(gel(B,k-1),t));
    2520           0 :       gcoeff(L,i,k-1) = gdiv(p1, Bk);
    2521             :     }
    2522             :   }
    2523          84 :   else if (!gequal0(la))
    2524             :   {
    2525          28 :     p1 = gdiv(la2, Bk);
    2526          28 :     gel(B,k+1) = gel(B,k) = p1;
    2527          28 :     for (i=k+2; i<=lx; i++) gel(B,i) = gdiv(gmul(p1,gel(B,i)),Bk);
    2528          28 :     for (i=k+1; i<lx; i++)
    2529           0 :       gcoeff(L,i,k-1) = gdiv(gmul(la,gcoeff(L,i,k-1)), Bk);
    2530          28 :     for (j=k+1; j<lx-1; j++)
    2531           0 :       for (i=j+1; i<lx; i++)
    2532           0 :         gcoeff(L,i,j) = gdiv(gmul(p1,gcoeff(L,i,j)), Bk);
    2533             :   }
    2534             :   else
    2535             :   {
    2536          56 :     gcoeff(L,k,k-1) = gen_0;
    2537          56 :     for (i=k+1; i<lx; i++)
    2538             :     {
    2539           0 :       gcoeff(L,i,k) = gcoeff(L,i,k-1);
    2540           0 :       gcoeff(L,i,k-1) = gen_0;
    2541             :     }
    2542          56 :     gel(B,k) = gel(B,k-1); fl[k] = 1; fl[k-1] = 0;
    2543             :   }
    2544         112 :   return 1;
    2545             : }
    2546             : 
    2547             : static void
    2548         168 : incrementalGSgen(GEN x, GEN L, GEN B, long k, GEN fl)
    2549             : {
    2550         168 :   GEN u = NULL; /* gcc -Wall */
    2551             :   long i, j;
    2552         420 :   for (j = 1; j <= k; j++)
    2553         252 :     if (j==k || fl[j])
    2554             :     {
    2555         252 :       u = gcoeff(x,k,j);
    2556         252 :       if (!is_extscalar_t(typ(u))) pari_err_TYPE("incrementalGSgen",u);
    2557         336 :       for (i=1; i<j; i++)
    2558          84 :         if (fl[i])
    2559             :         {
    2560          84 :           u = gsub(gmul(gel(B,i+1),u), gmul(gcoeff(L,k,i),gcoeff(L,j,i)));
    2561          84 :           u = gdiv(u, gel(B,i));
    2562             :         }
    2563         252 :       gcoeff(L,k,j) = u;
    2564             :     }
    2565         168 :   if (gequal0(u)) gel(B,k+1) = gel(B,k);
    2566             :   else
    2567             :   {
    2568         112 :     gel(B,k+1) = gcoeff(L,k,k); gcoeff(L,k,k) = gen_1; fl[k] = 1;
    2569             :   }
    2570         168 : }
    2571             : 
    2572             : static GEN
    2573         168 : lllgramallgen(GEN x, long flag)
    2574             : {
    2575         168 :   long lx = lg(x), i, j, k, l, n;
    2576             :   pari_sp av;
    2577             :   GEN B, L, h, fl;
    2578             :   int flc;
    2579             : 
    2580         168 :   n = lx-1; if (n<=1) return lll_trivial(x,flag);
    2581          84 :   if (lgcols(x) != lx) pari_err_DIM("lllgramallgen");
    2582             : 
    2583          84 :   fl = cgetg(lx, t_VECSMALL);
    2584             : 
    2585          84 :   av = avma;
    2586          84 :   B = scalarcol_shallow(gen_1, lx);
    2587          84 :   L = cgetg(lx,t_MAT);
    2588         252 :   for (j=1; j<lx; j++) { gel(L,j) = zerocol(n); fl[j] = 0; }
    2589             : 
    2590          84 :   h = matid(n);
    2591         252 :   for (i=1; i<lx; i++)
    2592         168 :     incrementalGSgen(x, L, B, i, fl);
    2593          84 :   flc = 0;
    2594          84 :   for(k=2;;)
    2595             :   {
    2596         196 :     if (REDgen(k, k-1, h, L, gel(B,k))) flc = 1;
    2597         196 :     if (do_SWAPgen(h, L, B, k, fl, &flc)) { if (k > 2) k--; }
    2598             :     else
    2599             :     {
    2600          84 :       for (l=k-2; l>=1; l--)
    2601           0 :         if (REDgen(k, l, h, L, gel(B,l+1))) flc = 1;
    2602          84 :       if (++k > n) break;
    2603             :     }
    2604         112 :     if (gc_needed(av,1))
    2605             :     {
    2606           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"lllgramallgen");
    2607           0 :       gerepileall(av,3,&B,&L,&h);
    2608             :     }
    2609             :   }
    2610         140 :   k=1; while (k<lx && !fl[k]) k++;
    2611          84 :   return lll_finish(h,k-1,flag);
    2612             : }
    2613             : 
    2614             : static GEN
    2615         168 : lllallgen(GEN x, long flag)
    2616             : {
    2617         168 :   pari_sp av = avma;
    2618         168 :   if (!(flag & LLL_GRAM)) x = gram_matrix(x);
    2619          84 :   else if (!RgM_is_square_mat(x)) pari_err_DIM("qflllgram");
    2620         168 :   return gerepilecopy(av, lllgramallgen(x, flag));
    2621             : }
    2622             : GEN
    2623          42 : lllgen(GEN x) { return lllallgen(x, LLL_IM); }
    2624             : GEN
    2625          42 : lllkerimgen(GEN x) { return lllallgen(x, LLL_ALL); }
    2626             : GEN
    2627          42 : lllgramgen(GEN x)  { return lllallgen(x, LLL_IM|LLL_GRAM); }
    2628             : GEN
    2629          42 : lllgramkerimgen(GEN x)  { return lllallgen(x, LLL_ALL|LLL_GRAM); }
    2630             : 
    2631             : static GEN
    2632       56123 : lllall(GEN x, long flag)
    2633       56123 : { pari_sp av = avma; return gerepilecopy(av, ZM_lll(x, LLLDFT, flag)); }
    2634             : GEN
    2635       14693 : lllint(GEN x) { return lllall(x, LLL_IM); }
    2636             : GEN
    2637          35 : lllkerim(GEN x) { return lllall(x, LLL_ALL); }
    2638             : GEN
    2639       41353 : lllgramint(GEN x)
    2640       41353 : { if (!RgM_is_square_mat(x)) pari_err_DIM("qflllgram");
    2641       41353 :   return lllall(x, LLL_IM | LLL_GRAM); }
    2642             : GEN
    2643          35 : lllgramkerim(GEN x)
    2644          35 : { if (!RgM_is_square_mat(x)) pari_err_DIM("qflllgram");
    2645          35 :   return lllall(x, LLL_ALL | LLL_GRAM); }
    2646             : 
    2647             : GEN
    2648     5324835 : lllfp(GEN x, double D, long flag)
    2649             : {
    2650     5324835 :   long n = lg(x)-1;
    2651     5324835 :   pari_sp av = avma;
    2652             :   GEN h;
    2653     5324835 :   if (n <= 1) return lll_trivial(x,flag);
    2654     4665020 :   if ((flag & LLL_GRAM) && !RgM_is_square_mat(x)) pari_err_DIM("qflllgram");
    2655     4665006 :   h = ZM_lll(RgM_rescale_to_int(x), D, flag);
    2656     4664983 :   return gerepilecopy(av, h);
    2657             : }
    2658             : 
    2659             : GEN
    2660        8949 : lllgram(GEN x) { return lllfp(x,LLLDFT,LLL_GRAM|LLL_IM); }
    2661             : GEN
    2662     1173381 : lll(GEN x) { return lllfp(x,LLLDFT,LLL_IM); }
    2663             : 
    2664             : GEN
    2665         301 : qflll0(GEN x, long flag)
    2666             : {
    2667         301 :   if (typ(x) != t_MAT) pari_err_TYPE("qflll",x);
    2668         301 :   switch(flag)
    2669             :   {
    2670          49 :     case 0: return lll(x);
    2671          63 :     case 1: RgM_check_ZM(x,"qflll"); return lllint(x);
    2672          49 :     case 2: RgM_check_ZM(x,"qflll"); return lllintpartial(x);
    2673           7 :     case 3: RgM_check_ZM(x,"qflll"); return lllall(x, LLL_INPLACE);
    2674          49 :     case 4: RgM_check_ZM(x,"qflll"); return lllkerim(x);
    2675          42 :     case 5: return lllkerimgen(x);
    2676          42 :     case 8: return lllgen(x);
    2677           0 :     default: pari_err_FLAG("qflll");
    2678             :   }
    2679             :   return NULL; /* LCOV_EXCL_LINE */
    2680             : }
    2681             : 
    2682             : GEN
    2683         245 : qflllgram0(GEN x, long flag)
    2684             : {
    2685         245 :   if (typ(x) != t_MAT) pari_err_TYPE("qflllgram",x);
    2686         245 :   switch(flag)
    2687             :   {
    2688          63 :     case 0: return lllgram(x);
    2689          49 :     case 1: RgM_check_ZM(x,"qflllgram"); return lllgramint(x);
    2690          49 :     case 4: RgM_check_ZM(x,"qflllgram"); return lllgramkerim(x);
    2691          42 :     case 5: return lllgramkerimgen(x);
    2692          42 :     case 8: return lllgramgen(x);
    2693           0 :     default: pari_err_FLAG("qflllgram");
    2694             :   }
    2695             :   return NULL; /* LCOV_EXCL_LINE */
    2696             : }
    2697             : 
    2698             : /********************************************************************/
    2699             : /**                                                                **/
    2700             : /**                   INTEGRAL KERNEL (LLL REDUCED)                **/
    2701             : /**                                                                **/
    2702             : /********************************************************************/
    2703             : static GEN
    2704          70 : kerint0(GEN M)
    2705             : {
    2706             :   /* return ZM_lll(M, LLLDFT, LLL_KER); */
    2707          70 :   GEN U, H = ZM_hnflll(M,&U,1);
    2708          70 :   long d = lg(M)-lg(H);
    2709          70 :   if (!d) return cgetg(1, t_MAT);
    2710          70 :   return ZM_lll(vecslice(U,1,d), LLLDFT, LLL_INPLACE);
    2711             : }
    2712             : GEN
    2713          42 : kerint(GEN M)
    2714             : {
    2715          42 :   pari_sp av = avma;
    2716          42 :   return gerepilecopy(av, kerint0(M));
    2717             : }
    2718             : /* OBSOLETE: use kerint */
    2719             : GEN
    2720          28 : matkerint0(GEN M, long flag)
    2721             : {
    2722          28 :   pari_sp av = avma;
    2723          28 :   if (typ(M) != t_MAT) pari_err_TYPE("matkerint",M);
    2724          28 :   M = Q_primpart(M);
    2725          28 :   RgM_check_ZM(M, "kerint");
    2726          28 :   switch(flag)
    2727             :   {
    2728          28 :     case 0:
    2729          28 :     case 1: return gerepilecopy(av, kerint0(M));
    2730           0 :     default: pari_err_FLAG("matkerint");
    2731             :   }
    2732             :   return NULL; /* LCOV_EXCL_LINE */
    2733             : }

Generated by: LCOV version 1.14