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.18.1 lcov report (development 30005-fc14bb602a) Lines: 1335 1641 81.4 %
Date: 2025-02-18 09:22:46 Functions: 125 130 96.2 %
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       45828 : RgM_is_square_mat(GEN x) { long l = lg(x); return l == 1 || l == lgcols(x); }
      22             : 
      23             : static long
      24     4180212 : ZM_is_upper(GEN R)
      25             : {
      26     4180212 :   long i,j, l = lg(R);
      27     4180212 :   if (l != lgcols(R)) return 0;
      28     8075694 :   for(i = 1; i < l; i++)
      29     8720464 :     for(j = 1; j < i; j++)
      30     4484088 :       if (signe(gcoeff(R,i,j))) return 0;
      31      252077 :   return 1;
      32             : }
      33             : 
      34             : static long
      35      606566 : ZM_is_knapsack(GEN R)
      36             : {
      37      606566 :   long i,j, l = lg(R);
      38      606566 :   if (l != lgcols(R)) return 0;
      39      843805 :   for(i = 2; i < l; i++)
      40     2902170 :     for(j = 1; j < l; j++)
      41     2664931 :       if ( i!=j && signe(gcoeff(R,i,j))) return 0;
      42       92364 :   return 1;
      43             : }
      44             : 
      45             : static long
      46     1182717 : ZM_is_lower(GEN R)
      47             : {
      48     1182717 :   long i,j, l = lg(R);
      49     1182717 :   if (l != lgcols(R)) return 0;
      50     2059270 :   for(i = 1; i < l; i++)
      51     2384311 :     for(j = 1; j < i; j++)
      52     1292083 :       if (signe(gcoeff(R,j,i))) return 0;
      53       34903 :   return 1;
      54             : }
      55             : 
      56             : static GEN
      57       34904 : RgM_flip(GEN R)
      58             : {
      59             :   GEN M;
      60             :   long i,j,l;
      61       34904 :   M = cgetg_copy(R, &l);
      62      181322 :   for(i = 1; i < l; i++)
      63             :   {
      64      146418 :     gel(M,i) = cgetg(l, t_COL);
      65      910402 :     for(j = 1; j < l; j++)
      66      763984 :       gmael(M,i,j) = gmael(R,l-i, l-j);
      67             :   }
      68       34904 :   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     4066913 : mpabscmp(GEN x, GEN y)
      89             : {
      90     4066913 :   return (typ(x)==t_INT) ? abscmpii(x,y) : abscmprr(x,y);
      91             : }
      92             : 
      93             : /****************************************************************************/
      94             : /***                             FLATTER                                  ***/
      95             : /****************************************************************************/
      96             : /* Implementation of "FLATTER" algorithm based on
      97             :  * <https://eprint.iacr.org/2023/237>
      98             :  * Fast Practical Lattice Reduction through Iterated Compression
      99             :  *
     100             :  * Keegan Ryan, University of California, San Diego
     101             :  * Nadia Heninger, University of California, San Diego. BA20230925 */
     102             : static long
     103     1336404 : drop(GEN R)
     104             : {
     105     1336404 :   long i, n = lg(R)-1;
     106     1336404 :   long s = 0, m = mpexpo(gcoeff(R, 1, 1));
     107     5403317 :   for (i = 2; i <= n; ++i)
     108             :   {
     109     4066913 :     if (mpabscmp(gcoeff(R, i, i), gcoeff(R, i - 1, i - 1)) >= 0)
     110             :     {
     111     2754727 :       s += m - mpexpo(gcoeff(R, i - 1, i - 1));
     112     2754727 :       m = mpexpo(gcoeff(R, i, i));
     113             :     }
     114             :   }
     115     1336404 :   s += m - mpexpo(gcoeff(R, n, n));
     116     1336406 :   return s;
     117             : }
     118             : 
     119             : static long
     120     1336406 : potential(GEN R)
     121             : {
     122     1336406 :   long i, n = lg(R)-1;
     123     1336406 :   long s = 0, mul = n-1;;
     124     6739727 :   for (i = 1; i <= n; i++, mul-=2) s += mul * mpexpo(gcoeff(R,i,i));
     125     1336405 :   return s;
     126             : }
     127             : 
     128             : /* U upper-triangular invertible:
     129             :  * Bound on the exponent of the condition number of U.
     130             :  * Algo 8.13 in Higham, Accuracy and stability of numercal algorithms. */
     131             : static long
     132     4691663 : condition_bound(GEN U, int lower)
     133             : {
     134     4691663 :   long n = lg(U)-1, e, i, j;
     135             :   GEN y;
     136     4691663 :   pari_sp av = avma;
     137     4691663 :   y = cgetg(n+1, t_VECSMALL);
     138     4691663 :   e = y[n] = -gexpo(gcoeff(U,n,n));
     139    18652305 :   for (i=n-1; i>0; i--)
     140             :   {
     141    13960641 :     long s = 0;
     142    49645769 :     for (j=i+1; j<=n; j++)
     143    35685128 :       s = maxss(s, (lower? gexpo(gcoeff(U,j,i)): gexpo(gcoeff(U,i,j))) + y[j]);
     144    13960641 :     y[i] = s - gexpo(gcoeff(U,i,i));
     145    13960641 :     e = maxss(e, y[i]);
     146             :   }
     147     4691664 :   return gc_long(av, gexpo(U) + e);
     148             : }
     149             : 
     150             : static long
     151     5151838 : gsisinv(GEN M)
     152             : {
     153     5151838 :   long i, l = lg(M);
     154    25920515 :   for (i = 1; i < l; ++i)
     155    20769085 :     if (! signe(gmael(M, i, i))) return 0;
     156     5151430 :   return 1;
     157             : }
     158             : 
     159             : INLINE long
     160     7418708 : nbits2prec64(long n)
     161             : {
     162     7418708 :   return nbits2prec(((n+63)>>6)<<6);
     163             : }
     164             : 
     165             : static long
     166     5809081 : spread(GEN R)
     167             : {
     168     5809081 :   long i, n = lg(R)-1, m = mpexpo(gcoeff(R, 1, 1)), M = m;
     169    23377406 :   for (i = 2; i <= n; ++i)
     170             :   {
     171    17568323 :     long e = mpexpo(gcoeff(R, i, i));
     172    17568322 :     if (e < m) m = e;
     173    17568322 :     if (e > M) M = e;
     174             :   }
     175     5809083 :   return M - m;
     176             : }
     177             : 
     178             : static long
     179     4691663 : GS_extraprec(GEN L, int lower)
     180             : {
     181     4691663 :   long C = condition_bound(L, lower), S = spread(L), n = lg(L)-1;
     182     4691665 :   return maxss(2*S+2*n, C-S-2*n); /* = 2*S + 2*n + maxss(0, C-3*S-4*n) */
     183             : }
     184             : 
     185             : static GEN
     186        2967 : RgM_Cholesky_dynprec(GEN M)
     187             : {
     188        2967 :   pari_sp ltop = avma;
     189             :   GEN L;
     190        2967 :   long minprec = lg(M) + 30, bitprec = minprec, prec;
     191             :   while (1)
     192        4877 :   {
     193             :     long mbitprec;
     194        7844 :     prec = nbits2prec64(bitprec);
     195        7844 :     L = RgM_Cholesky(RgM_gtofp(M, prec), prec); /* upper-triangular */
     196        7844 :     if (!L)
     197             :     {
     198        1458 :       bitprec *= 2;
     199        1458 :       set_avma(ltop);
     200        1458 :       continue;
     201             :     }
     202        6386 :     mbitprec = minprec + GS_extraprec(L, 0);
     203        6386 :     if (bitprec >= mbitprec)
     204        2967 :       break;
     205        3419 :     bitprec = maxss((4*bitprec)/3, mbitprec);
     206        3419 :     set_avma(ltop);
     207             :   }
     208        2967 :   return gerepilecopy(ltop, L);
     209             : }
     210             : 
     211             : static GEN
     212        1336 : gramschmidt_upper(GEN M)
     213             : {
     214        1336 :   long bitprec = lg(M)-1 + 31 + GS_extraprec(M, 0);
     215        1336 :   return RgM_gtofp(M, nbits2prec64(bitprec));
     216             : }
     217             : 
     218             : static GEN
     219     2672807 : gramschmidt_dynprec(GEN M)
     220             : {
     221     2672807 :   pari_sp ltop = avma;
     222     2672807 :   long minprec = lg(M) + 30, bitprec = minprec;
     223     2672807 :   if (ZM_is_upper(M)) return gramschmidt_upper(M);
     224             :   while (1)
     225     3610739 :   {
     226             :     GEN B, Q, L;
     227     6282210 :     long prec = nbits2prec64(bitprec), mbitprec;
     228     6282207 :     if (!QR_init(RgM_gtofp(M, prec), &B, &Q, &L, prec) || !gsisinv(L))
     229             :     {
     230     1598263 :       bitprec *= 2;
     231     1598263 :       set_avma(ltop);
     232     1598269 :       continue;
     233             :     }
     234     4683939 :     mbitprec = minprec + GS_extraprec(L, 1);
     235     4683943 :     if (bitprec >= mbitprec)
     236     2671473 :       return gerepilecopy(ltop, shallowtrans(L));
     237     2012470 :     bitprec = maxss((4*bitprec)/3, mbitprec);
     238     2012470 :     set_avma(ltop);
     239             :   }
     240             : }
     241             : /* return -T1 * round(T1^-1*(R1^-1*R2)*T3) */
     242             : static GEN
     243     1336406 : sizered(GEN T1, GEN T3, GEN R1, GEN R2)
     244             : {
     245     1336406 :   pari_sp ltop = avma;
     246             :   long e;
     247     1336406 :   return gerepileupto(ltop, ZM_mul(ZM_neg(T1), grndtoi(gmul(ZM_inv(T1,NULL),
     248             :          RgM_mul(RgM_mul(RgM_inv_upper(R1), R2), T3)), &e)));
     249             : }
     250             : 
     251             : static GEN
     252     1336403 : flat(GEN M, long flag, GEN *pt_T, long *pt_s, long *pt_pot)
     253             : {
     254     1336403 :   pari_sp ltop = avma;
     255             :   GEN R, R1, R2, R3, T1, T2, T3, T, S;
     256     1336403 :   long k = lg(M)-1, n = k>>1, n2 = k - n, m = n>>1;
     257     1336403 :   long keepfirst = flag & LLL_KEEP_FIRST, inplace = flag & LLL_INPLACE;
     258             :   /* for k = 3, we want n = 1; n2  = 2; m = 0 */
     259             :   /* for k = 5,         n = 2; n2 = 3; m = 1 */
     260     1336403 :   R = gramschmidt_dynprec(M);
     261     1336405 :   R1 = matslice(R, 1, n, 1, n);
     262     1336403 :   R2 = matslice(R, 1, n, n + 1, k);
     263     1336404 :   R3 = matslice(R, n + 1, k, n + 1, k);
     264     1336405 :   T1 = lllfp(R1, 0.99, LLL_IM| LLL_UPPER| LLL_NOCERTIFY| (keepfirst ? LLL_KEEP_FIRST: 0));
     265     1336405 :   T3 = lllfp(R3, 0.99, LLL_IM| LLL_UPPER| LLL_NOCERTIFY);
     266     1336406 :   T2 = sizered(T1, T3, R1, R2);
     267     1336405 :   T = shallowmatconcat(mkmat22(T1,T2,gen_0,T3));
     268     1336405 :   M = ZM_mul(M, T);
     269     1336404 :   R = gramschmidt_dynprec(M);
     270     1336406 :   R3 = matslice(R, m + 1, m + n2, m + 1, m + n2);
     271     1336405 :   T3 = lllfp(R3, 0.99, LLL_IM| LLL_UPPER| LLL_NOCERTIFY);
     272     2672809 :   S = shallowmatconcat(diagonal(
     273      576108 :        m == 0     ? mkvec2(T3, matid(k - m - n2))
     274           0 :      : m+n2 == k  ? mkvec2(matid(m), T3)
     275      760296 :                   : mkvec3(matid(m), T3, matid(k - m - n2))));
     276     1336406 :   M = ZM_mul(M, S);
     277     1336404 :   if (!inplace) *pt_T = ZM_mul(T, S);
     278     1336404 :   *pt_s = drop(R);
     279     1336406 :   *pt_pot = potential(R);
     280     1336405 :   return gc_all(ltop, inplace ? 1: 2, &M, pt_T);
     281             : }
     282             : 
     283             : static void
     284           0 : dbg_flatter(pari_timer *ti, long n, long i, long lti, double t, double pot2)
     285             : {
     286           0 :   double s = t / n, p = pot2 / (n*(n+1));
     287             :   const char *str;
     288           0 :   if (i == -1)
     289           0 :     str = (i == lti)? "final"
     290           0 :                     : stack_sprintf("steps %ld-final", lti);
     291             :   else
     292           0 :     str = (i == lti)? stack_sprintf("step %ld", i)
     293           0 :                     : stack_sprintf("steps %ld-%ld", lti, i);
     294           0 :   timer_printf(ti, "FLATTER, dim %ld, %s: \t slope=%0.10g \t pot=%0.10g",
     295             :                n, str, s, p);
     296           0 : }
     297             : 
     298             : static GEN
     299      618594 : ZM_flatter(GEN M, long flag)
     300             : {
     301      618594 :   pari_sp av = avma;
     302      618594 :   long i, n = lg(M)-1, s = -1, lti = 1, pot = LONG_MAX;
     303      618594 :   GEN T = NULL;
     304             :   pari_timer ti;
     305      618594 :   long inplace = flag & LLL_INPLACE, cert = !(flag & LLL_NOCERTIFY);
     306             : 
     307      618594 :   if (DEBUGLEVEL>=3)
     308             :   {
     309           0 :     timer_start(&ti);
     310           0 :     if (cert) err_printf("FLATTER dim = %ld size = %ld\n", n, ZM_max_expi(M));
     311             :   }
     312      618594 :   for (i = 1;;i++)
     313      717809 :   {
     314             :     long t, pot2;
     315     1336403 :     GEN U, M2 = flat(M, flag, &U, &t, &pot2);
     316     1336406 :     if (t == 0) { s = t; break; }
     317      761144 :     if (s >= 0)
     318             :     {
     319      435870 :       if (s == t && pot>=pot2) break;
     320      392536 :       if (s < t && i > 20)
     321             :       {
     322           0 :         if (DEBUGLEVEL >= 3) err_printf("BACK:%ld:%ld:%g\n", n, i, s);
     323           0 :         break;
     324             :       }
     325             :     }
     326      717810 :     if (DEBUGLEVEL>=3 && (cert || timer_get(&ti) > 1000))
     327           0 :       dbg_flatter(&ti, n, i, lti, t, pot2);
     328      717810 :     s = t;
     329      717810 :     pot = pot2;
     330      717810 :     M = M2;
     331      717810 :     if (!inplace)
     332             :     {
     333      690178 :       T = T? ZM_mul(T, U): U;
     334      690177 :       if (gc_needed(av, 1)) gerepileall(av, 2, &M, &T);
     335             :     }
     336             :     else
     337       27632 :       if (gc_needed(av, 1)) M = gerepilecopy(av, M);
     338             :   }
     339      618596 :   if (DEBUGLEVEL>=3 && (cert || timer_get(&ti) > 1000))
     340           0 :     dbg_flatter(&ti, n, -1, i == lti? -1: lti, s, pot);
     341      618596 :   if (!inplace)
     342             :   {
     343      604588 :     if (!T) return gc_NULL(av);
     344      311434 :     return gerepilecopy(av, T);
     345             :   }
     346       14008 :   return  gerepilecopy(av, M);
     347             : }
     348             : 
     349             : static GEN
     350      616580 : ZM_flatter_rank(GEN M, long rank, long flag)
     351             : {
     352             :   pari_timer ti;
     353      616580 :   pari_sp av = avma;
     354      616580 :   GEN T = NULL;
     355      616580 :   long i, n = lg(M)-1, sm = LONG_MAX;
     356      616580 :   long inplace = flag & LLL_INPLACE;
     357             : 
     358      616580 :   if (rank == n) return ZM_flatter(M, flag);
     359        3785 :   if (DEBUGLEVEL>=3) timer_start(&ti);
     360        3785 :   for (i = 1;; i++)
     361        2014 :   {
     362        5799 :     GEN S = ZM_flatter(vconcat(gshift(M,i),matid(n)), flag);
     363             :     long s;
     364        5799 :     if (!S || (s = expi(gnorml2(S))) >= sm) break;
     365        2014 :     sm = s;
     366        2014 :     if (DEBUGLEVEL>=3) timer_printf(&ti,"FLATTERRANK step %ld: %ld",i,sm);
     367        2014 :     T = T? ZM_mul(T, S): S;
     368        2014 :     M = ZM_mul(M, S);
     369        2014 :     if (gc_needed(av, 1)) gerepileall(av, 2, &M, &T);
     370             :   }
     371        3785 :   if (!inplace)
     372             :   {
     373        3778 :     if (!T) { set_avma(av); return matid(n); }
     374        1951 :     return gerepilecopy(av, T);
     375             :   }
     376           7 :   return  gerepilecopy(av, M);
     377             : }
     378             : 
     379             : static GEN
     380        2967 : flattergram_i(GEN M, long flag)
     381             : {
     382        2967 :   pari_sp av = avma;
     383        2967 :   GEN T, R = RgM_Cholesky_dynprec(M);
     384        2967 :   T = lllfp(R, 0.99, LLL_IM|LLL_UPPER|LLL_NOCERTIFY | (flag&LLL_KEEP_FIRST));
     385        2967 :   return gerepileupto(av, T);
     386             : }
     387             : 
     388             : static void
     389           0 : dbg_flattergram(pari_timer *t, long n, long i, long s)
     390           0 : { timer_printf(t, "FLATTERGRAM, dim %ld step %ld, slope=%0.10g", n, i,
     391           0 :                ((double)s)/n); }
     392             : /* return base change, NULL if identity */
     393             : static GEN
     394         961 : ZM_flattergram(GEN M, long flag)
     395             : {
     396         961 :   pari_sp av = avma;
     397         961 :   GEN T = NULL;
     398         961 :   long i, n = lg(M)-1, s = -1;
     399             : 
     400             :   pari_timer ti;
     401         961 :   if (DEBUGLEVEL>=3)
     402             :   {
     403           0 :     timer_start(&ti);
     404           0 :     err_printf("FLATTERGRAM dim = %ld size = %ld\n", n, ZM_max_expi(M));
     405             :   }
     406         961 :   for (i = 1;; i++)
     407        2006 :   {
     408        2967 :     GEN S = flattergram_i(M, flag);
     409        2967 :     long t = expi(gnorml2(S));
     410        2967 :     if (t == 0) { s = t;  break; }
     411        2967 :     if (s)
     412             :     {
     413        2967 :       double st = s - t;
     414        2967 :       if (st == 0) break;
     415        2006 :       if (st < 0 && i > 20)
     416             :       {
     417           0 :         if (DEBUGLEVEL >= 3)
     418           0 :           err_printf("BACK:%ld:%ld:%0.10g\n", n, i, ((double)s)/n);
     419           0 :         break;
     420             :       }
     421             :     }
     422        2006 :     T = T? ZM_mul(T, S): S;
     423        2006 :     M = qf_ZM_apply(M, S);
     424        2006 :     s = t;
     425        2006 :     if (DEBUGLEVEL >= 3) dbg_flattergram(&ti, n, i, s);
     426        2006 :     if (gc_needed(av, 1)) gerepileall(av, 2, &M, &T);
     427             :   }
     428         961 :   if (DEBUGLEVEL >= 3) dbg_flattergram(&ti, n, i, s);
     429         961 :   if (!T && ZM_isidentity(T)) return gc_NULL(av);
     430         961 :   return gerepilecopy(av, T);
     431             : }
     432             : 
     433             : /* return base change, NULL if identity */
     434             : static GEN
     435         961 : ZM_flattergram_rank(GEN M, long rank, long flag)
     436             : {
     437             :   pari_timer ti;
     438         961 :   pari_sp av = avma;
     439         961 :   GEN T = NULL;
     440         961 :   long i, n = lg(M)-1;
     441         961 :   if (rank == n) return ZM_flattergram(M, flag);
     442           0 :   if (DEBUGLEVEL>=3) timer_start(&ti);
     443           0 :   for (i = 1;; i++)
     444           0 :   {
     445           0 :     GEN S = ZM_flattergram(RgM_Rg_add(gshift(M, i), gen_1), flag);
     446           0 :     if (DEBUGLEVEL>=3)
     447           0 :       timer_printf(&ti,"FLATTERGRAMRANK step %ld: %ld",i,expi(gnorml2(S)));
     448           0 :     if (!S) break;
     449           0 :     T = T? ZM_mul(T, S): S;
     450           0 :     M = qf_ZM_apply(M, S);
     451           0 :     if (gc_needed(av, 1)) gerepileall(av, 2, &M, &T);
     452             :   }
     453           0 :   if (!T || ZM_isidentity(T)) return gc_NULL(av);
     454           0 :   return gerepilecopy(av, T);
     455             : }
     456             : 
     457             : /* round to closest integer (as a double). If |a| >= 2^52, return it */
     458             : static double
     459    10721532 : pari_rint(double a)
     460             : {
     461             : #ifdef HAS_RINT
     462    10721532 :   return rint(a);
     463             : #else
     464             :   const double pow2 = 4.5035996273704960e+15; /* 2^52 */
     465             :   double r, fa = fabs(a);
     466             :   if (fa >= pow2) return a;
     467             :   r = (pow2 + fa) - pow2;
     468             :   if (a < 0) r = -r;
     469             :   return r;
     470             : #endif
     471             : }
     472             : 
     473             : /* default quality ratio for LLL */
     474             : static const double LLLDFT = 0.99;
     475             : 
     476             : /* assume flag & (LLL_KER|LLL_IM|LLL_ALL). LLL_INPLACE implies LLL_IM */
     477             : static GEN
     478      770184 : lll_trivial(GEN x, long flag)
     479             : {
     480      770184 :   if (lg(x) == 1)
     481             :   { /* dim x = 0 */
     482       15512 :     if (! (flag & LLL_ALL)) return cgetg(1,t_MAT);
     483          28 :     retmkvec2(cgetg(1,t_MAT), cgetg(1,t_MAT));
     484             :   }
     485             :   /* dim x = 1 */
     486      754672 :   if (gequal0(gel(x,1)))
     487             :   {
     488         146 :     if (flag & LLL_KER) return matid(1);
     489         146 :     if (flag & (LLL_IM|LLL_INPLACE)) return cgetg(1,t_MAT);
     490          28 :     retmkvec2(matid(1), cgetg(1,t_MAT));
     491             :   }
     492      754524 :   if (flag & LLL_INPLACE) return gcopy(x);
     493      651123 :   if (flag & LLL_KER) return cgetg(1,t_MAT);
     494      651123 :   if (flag & LLL_IM)  return matid(1);
     495          28 :   retmkvec2(cgetg(1,t_MAT), (flag & LLL_GRAM)? gcopy(x): matid(1));
     496             : }
     497             : 
     498             : /* vecslice(x,#x-k,#x) in place. Works for t_MAT, t_VEC/t_COL */
     499             : static GEN
     500     2056032 : vectail_inplace(GEN x, long k)
     501             : {
     502     2056032 :   if (!k) return x;
     503       57766 :   x[k] = ((ulong)x[0] & ~LGBITS) | _evallg(lg(x) - k);
     504       57766 :   return x + k;
     505             : }
     506             : 
     507             : /* k = dim Kernel */
     508             : static GEN
     509     2130215 : lll_finish(GEN h, long k, long flag)
     510             : {
     511             :   GEN g;
     512     2130215 :   if (!(flag & (LLL_IM|LLL_KER|LLL_ALL|LLL_INPLACE))) return h;
     513     2056057 :   if (flag & (LLL_IM|LLL_INPLACE)) return vectail_inplace(h, k);
     514          96 :   if (flag & LLL_KER) { setlg(h,k+1); return h; }
     515          68 :   g = vecslice(h,1,k); /* done first: vectail_inplace kills h */
     516          70 :   return mkvec2(g, vectail_inplace(h, k));
     517             : }
     518             : 
     519             : /* y * z * 2^e, e >= 0; y,z t_INT */
     520             : INLINE GEN
     521      317482 : mulshift(GEN y, GEN z, long e)
     522             : {
     523      317482 :   long ly = lgefint(y), lz;
     524             :   pari_sp av;
     525             :   GEN t;
     526      317482 :   if (ly == 2) return gen_0;
     527       46556 :   lz = lgefint(z);
     528       46556 :   av = avma; (void)new_chunk(ly+lz+nbits2lg(e)); /* HACK */
     529       46556 :   t = mulii(z, y);
     530       46556 :   set_avma(av); return shifti(t, e);
     531             : }
     532             : 
     533             : /* x - y * z * 2^e, e >= 0; x,y,z t_INT */
     534             : INLINE GEN
     535     1502834 : submulshift(GEN x, GEN y, GEN z, long e)
     536             : {
     537     1502834 :   long lx = lgefint(x), ly, lz;
     538             :   pari_sp av;
     539             :   GEN t;
     540     1502834 :   if (!e) return submulii(x, y, z);
     541     1492887 :   if (lx == 2) { t = mulshift(y, z, e); togglesign(t); return t; }
     542     1194416 :   ly = lgefint(y);
     543     1194416 :   if (ly == 2) return icopy(x);
     544      956829 :   lz = lgefint(z);
     545      956829 :   av = avma; (void)new_chunk(lx+ly+lz+nbits2lg(e)); /* HACK */
     546      956829 :   t = shifti(mulii(z, y), e);
     547      956829 :   set_avma(av); return subii(x, t);
     548             : }
     549             : static void
     550    18517758 : subzi(GEN *a, GEN b)
     551             : {
     552    18517758 :   pari_sp av = avma;
     553    18517758 :   b = subii(*a, b);
     554    18517747 :   if (lgefint(b)<=lg(*a) && isonstack(*a)) { affii(b,*a); set_avma(av); }
     555     2104191 :   else *a = b;
     556    18517811 : }
     557             : 
     558             : static void
     559    17776086 : addzi(GEN *a, GEN b)
     560             : {
     561    17776086 :   pari_sp av = avma;
     562    17776086 :   b = addii(*a, b);
     563    17776080 :   if (lgefint(b)<=lg(*a) && isonstack(*a)) { affii(b,*a); set_avma(av); }
     564     1876697 :   else *a = b;
     565    17776143 : }
     566             : 
     567             : /* x - u*y * 2^e */
     568             : INLINE GEN
     569     4178991 : submuliu2n(GEN x, GEN y, ulong u, long e)
     570             : {
     571             :   pari_sp av;
     572     4178991 :   long ly = lgefint(y);
     573     4178991 :   if (ly == 2) return x;
     574     2866785 :   av = avma;
     575     2866785 :   (void)new_chunk(3+ly+lgefint(x)+nbits2lg(e)); /* HACK */
     576     2868032 :   y = shifti(mului(u,y), e);
     577     2867227 :   set_avma(av); return subii(x, y);
     578             : }
     579             : /* *x -= u*y * 2^e */
     580             : INLINE void
     581      262549 : submulzu2n(GEN *x, GEN y, ulong u, long e)
     582             : {
     583             :   pari_sp av;
     584      262549 :   long ly = lgefint(y);
     585      262549 :   if (ly == 2) return;
     586      172674 :   av = avma;
     587      172674 :   (void)new_chunk(3+ly+lgefint(*x)+nbits2lg(e)); /* HACK */
     588      172674 :   y = shifti(mului(u,y), e);
     589      172674 :   set_avma(av); return subzi(x, y);
     590             : }
     591             : 
     592             : /* x + u*y * 2^e */
     593             : INLINE GEN
     594     4093435 : addmuliu2n(GEN x, GEN y, ulong u, long e)
     595             : {
     596             :   pari_sp av;
     597     4093435 :   long ly = lgefint(y);
     598     4093435 :   if (ly == 2) return x;
     599     2812592 :   av = avma;
     600     2812592 :   (void)new_chunk(3+ly+lgefint(x)+nbits2lg(e)); /* HACK */
     601     2813779 :   y = shifti(mului(u,y), e);
     602     2813041 :   set_avma(av); return addii(x, y);
     603             : }
     604             : 
     605             : /* *x += u*y * 2^e */
     606             : INLINE void
     607      271563 : addmulzu2n(GEN *x, GEN y, ulong u, long e)
     608             : {
     609             :   pari_sp av;
     610      271563 :   long ly = lgefint(y);
     611      271563 :   if (ly == 2) return;
     612      178632 :   av = avma;
     613      178632 :   (void)new_chunk(3+ly+lgefint(*x)+nbits2lg(e)); /* HACK */
     614      178632 :   y = shifti(mului(u,y), e);
     615      178632 :   set_avma(av); return addzi(x, y);
     616             : }
     617             : 
     618             : /* n < 10; gerepileall supporting &NULL arguments. Maybe rename and export ? */
     619             : INLINE void
     620        4722 : gc_lll(pari_sp av, int n, ...)
     621             : {
     622             :   int i, j;
     623             :   GEN *gptr[10];
     624             :   size_t s;
     625        4722 :   va_list a; va_start(a, n);
     626       14166 :   for (i=j=0; i<n; i++)
     627             :   {
     628        9444 :     GEN *x = va_arg(a,GEN*);
     629        9444 :     if (*x) { gptr[j++] = x; *x = (GEN)copy_bin(*x); }
     630             :   }
     631        4722 :   va_end(a); set_avma(av);
     632       11764 :   for (--j; j>=0; j--) *gptr[j] = bin_copy((GENbin*)*gptr[j]);
     633        4722 :   s = pari_mainstack->top - pari_mainstack->bot;
     634             :   /* size of saved objects ~ stacksize / 4 => overflow */
     635        4722 :   if (av - avma > (s >> 2))
     636             :   {
     637           0 :     size_t t = avma - pari_mainstack->bot;
     638           0 :     av = avma; new_chunk((s + t) / sizeof(long)); set_avma(av); /* double */
     639             :   }
     640        4722 : }
     641             : 
     642             : /********************************************************************/
     643             : /**                                                                **/
     644             : /**                   FPLLL (adapted from D. Stehle's code)        **/
     645             : /**                                                                **/
     646             : /********************************************************************/
     647             : /* Babai* and fplll* are a conversion to libpari API and data types
     648             :    of fplll-1.3 by Damien Stehle'.
     649             : 
     650             :   Copyright 2005, 2006 Damien Stehle'.
     651             : 
     652             :   This program is free software; you can redistribute it and/or modify it
     653             :   under the terms of the GNU General Public License as published by the
     654             :   Free Software Foundation; either version 2 of the License, or (at your
     655             :   option) any later version.
     656             : 
     657             :   This program implements ideas from the paper "Floating-point LLL Revisited",
     658             :   by Phong Nguyen and Damien Stehle', in the Proceedings of Eurocrypt'2005,
     659             :   Springer-Verlag; and was partly inspired by Shoup's NTL library:
     660             :   http://www.shoup.net/ntl/ */
     661             : 
     662             : /* x t_REAL, |x| >= 1/2. Test whether |x| <= 3/2 */
     663             : static int
     664      394782 : absrsmall2(GEN x)
     665             : {
     666      394782 :   long e = expo(x), l, i;
     667      394782 :   if (e < 0) return 1;
     668      198103 :   if (e > 0 || (ulong)x[2] > (3UL << (BITS_IN_LONG-2))) return 0;
     669             :   /* line above assumes l > 2. OK since x != 0 */
     670       73398 :   l = lg(x); for (i = 3; i < l; i++) if (x[i]) return 0;
     671       62863 :   return 1;
     672             : }
     673             : /* x t_REAL; test whether |x| <= 1/2 */
     674             : static int
     675      687537 : absrsmall(GEN x)
     676             : {
     677             :   long e, l, i;
     678      687537 :   if (!signe(x)) return 1;
     679      683370 :   e = expo(x); if (e < -1) return 1;
     680      400192 :   if (e > -1 || (ulong)x[2] > HIGHBIT) return 0;
     681        6123 :   l = lg(x); for (i = 3; i < l; i++) if (x[i]) return 0;
     682        5410 :   return 1;
     683             : }
     684             : 
     685             : static void
     686    31771143 : rotate(GEN A, long k2, long k)
     687             : {
     688             :   long i;
     689    31771143 :   GEN B = gel(A,k2);
     690   101538828 :   for (i = k2; i > k; i--) gel(A,i) = gel(A,i-1);
     691    31771143 :   gel(A,k) = B;
     692    31771143 : }
     693             : 
     694             : /************************* FAST version (double) ************************/
     695             : #define dmael(x,i,j) ((x)[i][j])
     696             : #define del(x,i) ((x)[i])
     697             : 
     698             : static double *
     699    34500158 : cget_dblvec(long d)
     700    34500158 : { return (double*) stack_malloc_align(d*sizeof(double), sizeof(double)); }
     701             : 
     702             : static double **
     703     8278850 : cget_dblmat(long d) { return (double **) cgetg(d, t_VECSMALL); }
     704             : 
     705             : static double
     706   161934464 : itodbl_exp(GEN x, long *e)
     707             : {
     708   161934464 :   pari_sp av = avma;
     709   161934464 :   GEN r = itor(x,DEFAULTPREC);
     710   161924931 :   *e = expo(r); setexpo(r,0);
     711   161921379 :   return gc_double(av, rtodbl(r));
     712             : }
     713             : 
     714             : static double
     715   117788307 : dbldotproduct(double *x, double *y, long n)
     716             : {
     717             :   long i;
     718   117788307 :   double sum = del(x,1) * del(y,1);
     719  1383700306 :   for (i=2; i<=n; i++) sum += del(x,i) * del(y,i);
     720   117788307 :   return sum;
     721             : }
     722             : 
     723             : static double
     724     2441027 : dbldotsquare(double *x, long n)
     725             : {
     726             :   long i;
     727     2441027 :   double sum = del(x,1) * del(x,1);
     728     8092494 :   for (i=2; i<=n; i++) sum += del(x,i) * del(x,i);
     729     2441027 :   return sum;
     730             : }
     731             : 
     732             : static long
     733    24721788 : set_line(double *appv, GEN v, long n)
     734             : {
     735    24721788 :   long i, maxexp = 0;
     736    24721788 :   pari_sp av = avma;
     737    24721788 :   GEN e = cgetg(n+1, t_VECSMALL);
     738   186644477 :   for (i = 1; i <= n; i++)
     739             :   {
     740   161932038 :     del(appv,i) = itodbl_exp(gel(v,i), e+i);
     741   161921334 :     if (e[i] > maxexp) maxexp = e[i];
     742             :   }
     743   186702172 :   for (i = 1; i <= n; i++) del(appv,i) = ldexp(del(appv,i), e[i]-maxexp);
     744    24712439 :   set_avma(av); return maxexp;
     745             : }
     746             : 
     747             : static void
     748    34460854 : dblrotate(double **A, long k2, long k)
     749             : {
     750             :   long i;
     751    34460854 :   double *B = del(A,k2);
     752   109179516 :   for (i = k2; i > k; i--) del(A,i) = del(A,i-1);
     753    34460854 :   del(A,k) = B;
     754    34460854 : }
     755             : /* update G[kappa][i] from appB */
     756             : static void
     757    22484886 : setG_fast(double **appB, long n, double **G, long kappa, long a, long b)
     758             : { long i;
     759   101398199 :   for (i = a; i <= b; i++)
     760    78913581 :     dmael(G,kappa,i) = dbldotproduct(del(appB,kappa), del(appB,i), n);
     761    22484618 : }
     762             : /* update G[i][kappa] from appB */
     763             : static void
     764    16926760 : setG2_fast(double **appB, long n, double **G, long kappa, long a, long b)
     765             : { long i;
     766    55803739 :   for (i = a; i <= b; i++)
     767    38877007 :     dmael(G,i,kappa) = dbldotproduct(del(appB,kappa), del(appB,i), n);
     768    16926732 : }
     769             : const long EX0 = -2; /* uninitialized; any value less than expo(0.51) = -1 */
     770             : 
     771             : #ifdef LONG_IS_64BIT
     772             : typedef long s64;
     773             : #define addmuliu64_inplace addmuliu_inplace
     774             : #define submuliu64_inplace submuliu_inplace
     775             : #define submuliu642n submuliu2n
     776             : #define addmuliu642n addmuliu2n
     777             : #else
     778             : typedef long long s64;
     779             : typedef unsigned long long u64;
     780             : 
     781             : INLINE GEN
     782    19746472 : u64toi(u64 x)
     783             : {
     784             :   GEN y;
     785             :   ulong h;
     786    19746472 :   if (!x) return gen_0;
     787    19746472 :   h = x>>32;
     788    19746472 :   if (!h) return utoipos(x);
     789     1135979 :   y = cgetipos(4);
     790     1135979 :   *int_LSW(y) = x&0xFFFFFFFF;
     791     1135979 :   *int_MSW(y) = x>>32;
     792     1135979 :   return y;
     793             : }
     794             : 
     795             : INLINE GEN
     796      668825 : u64toineg(u64 x)
     797             : {
     798             :   GEN y;
     799             :   ulong h;
     800      668825 :   if (!x) return gen_0;
     801      668825 :   h = x>>32;
     802      668825 :   if (!h) return utoineg(x);
     803      668825 :   y = cgetineg(4);
     804      668825 :   *int_LSW(y) = x&0xFFFFFFFF;
     805      668825 :   *int_MSW(y) = x>>32;
     806      668825 :   return y;
     807             : }
     808             : INLINE GEN
     809     9502895 : addmuliu64_inplace(GEN x, GEN y, u64 u) { return addmulii(x, y, u64toi(u)); }
     810             : 
     811             : INLINE GEN
     812     9557744 : submuliu64_inplace(GEN x, GEN y, u64 u) { return submulii(x, y, u64toi(u)); }
     813             : 
     814             : INLINE GEN
     815      668825 : addmuliu642n(GEN x, GEN y, u64 u, long e) { return submulshift(x, y, u64toineg(u), e); }
     816             : 
     817             : INLINE GEN
     818      685833 : submuliu642n(GEN x, GEN y, u64 u, long e) { return submulshift(x, y, u64toi(u), e); }
     819             : 
     820             : #endif
     821             : 
     822             : /* Babai's Nearest Plane algorithm (iterative); see Babai() */
     823             : static int
     824    30050064 : Babai_fast(pari_sp av, long kappa, GEN *pB, GEN *pU, double **mu, double **r,
     825             :            double *s, double **appB, GEN expoB, double **G,
     826             :            long a, long zeros, long maxG, double eta)
     827             : {
     828    30050064 :   GEN B = *pB, U = *pU;
     829    30050064 :   const long n = nbrows(B), d = U ? lg(U)-1: 0;
     830    30049964 :   long k, aa = (a > zeros)? a : zeros+1;
     831    30049964 :   long emaxmu = EX0, emax2mu = EX0;
     832             :   s64 xx;
     833    30049964 :   int did_something = 0;
     834             :   /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
     835             : 
     836    17133791 :   for (;;) {
     837    47183755 :     int go_on = 0;
     838    47183755 :     long i, j, emax3mu = emax2mu;
     839             : 
     840    47183755 :     if (gc_needed(av,2))
     841             :     {
     842         197 :       if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
     843         197 :       gc_lll(av,2,&B,&U);
     844             :     }
     845             :     /* Step2: compute the GSO for stage kappa */
     846    47183079 :     emax2mu = emaxmu; emaxmu = EX0;
     847   180886617 :     for (j=aa; j<kappa; j++)
     848             :     {
     849   133704063 :       double g = dmael(G,kappa,j);
     850   574228840 :       for (k = zeros+1; k < j; k++) g -= dmael(mu,j,k) * dmael(r,kappa,k);
     851   133704063 :       dmael(r,kappa,j) = g;
     852   133704063 :       dmael(mu,kappa,j) = dmael(r,kappa,j) / dmael(r,j,j);
     853   133704063 :       emaxmu = maxss(emaxmu, expoB[kappa]-expoB[j]);
     854             :     }
     855             :     /* maxmu doesn't decrease fast enough */
     856    47182554 :     if (emax3mu != EX0 && emax3mu <= emax2mu + 5) {*pB = B; *pU = U; return 1;}
     857             : 
     858   167943622 :     for (j=kappa-1; j>zeros; j--)
     859             :     {
     860   137899154 :       double tmp = fabs(ldexp (dmael(mu,kappa,j), expoB[kappa]-expoB[j]));
     861   137899154 :       if (tmp>eta) { go_on = 1; break; }
     862             :     }
     863             : 
     864             :     /* Step3--5: compute the X_j's  */
     865    47178421 :     if (go_on)
     866    77820611 :       for (j=kappa-1; j>zeros; j--)
     867             :       { /* The code below seemingly handles U = NULL, but in this case d = 0 */
     868    60688998 :         int e = expoB[j] - expoB[kappa];
     869    60688998 :         double tmp = ldexp(dmael(mu,kappa,j), -e), atmp = fabs(tmp);
     870             :         /* tmp = Inf is allowed */
     871    60688998 :         if (atmp <= .5) continue; /* size-reduced */
     872    34006022 :         if (gc_needed(av,2))
     873             :         {
     874         346 :           if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
     875         346 :           gc_lll(av,2,&B,&U);
     876             :         }
     877    34007768 :         did_something = 1;
     878             :         /* we consider separately the case |X| = 1 */
     879    34007768 :         if (atmp <= 1.5)
     880             :         {
     881    22914752 :           if (dmael(mu,kappa,j) > 0) { /* in this case, X = 1 */
     882    46910076 :             for (k=zeros+1; k<j; k++)
     883    35214548 :               dmael(mu,kappa,k) -= ldexp(dmael(mu,j,k), e);
     884   158294051 :             for (i=1; i<=n; i++)
     885   146600327 :               gmael(B,kappa,i) = subii(gmael(B,kappa,i), gmael(B,j,i));
     886   104532374 :             for (i=1; i<=d; i++)
     887    92838602 :               gmael(U,kappa,i) = subii(gmael(U,kappa,i), gmael(U,j,i));
     888             :           } else { /* otherwise X = -1 */
     889    46075810 :             for (k=zeros+1; k<j; k++)
     890    34856586 :               dmael(mu,kappa,k) += ldexp(dmael(mu,j,k), e);
     891   155602454 :             for (i=1; i<=n; i++)
     892   144385518 :               gmael(B,kappa,i) = addii(gmael(B,kappa,i), gmael(B,j,i));
     893   101814114 :             for (i=1; i<=d; i++)
     894    90597060 :               gmael(U,kappa,i) = addii(gmael(U,kappa,i), gmael(U,j,i));
     895             :           }
     896    22910826 :           continue;
     897             :         }
     898             :         /* we have |X| >= 2 */
     899    11093016 :         if (atmp < 9007199254740992.)
     900             :         {
     901    10254293 :           tmp = pari_rint(tmp);
     902    24479394 :           for (k=zeros+1; k<j; k++)
     903    14225086 :             dmael(mu,kappa,k) -= ldexp(tmp * dmael(mu,j,k), e);
     904    10254308 :           xx = (s64) tmp;
     905    10254308 :           if (xx > 0) /* = xx */
     906             :           {
     907    46060708 :             for (i=1; i<=n; i++)
     908    40903121 :               gmael(B,kappa,i) = submuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
     909    32767204 :             for (i=1; i<=d; i++)
     910    27609585 :               gmael(U,kappa,i) = submuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
     911             :           }
     912             :           else /* = -xx */
     913             :           {
     914    45728452 :             for (i=1; i<=n; i++)
     915    40631960 :               gmael(B,kappa,i) = addmuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), -xx);
     916    32397545 :             for (i=1; i<=d; i++)
     917    27300991 :               gmael(U,kappa,i) = addmuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), -xx);
     918             :           }
     919             :         }
     920             :         else
     921             :         {
     922             :           int E;
     923      838723 :           xx = (s64) ldexp(frexp(dmael(mu,kappa,j), &E), 53);
     924      838723 :           E -= e + 53;
     925      838723 :           if (E <= 0)
     926             :           {
     927           0 :             xx = xx << -E;
     928           0 :             for (k=zeros+1; k<j; k++)
     929           0 :               dmael(mu,kappa,k) -= ldexp(((double)xx) * dmael(mu,j,k), e);
     930           0 :             if (xx > 0) /* = xx */
     931             :             {
     932           0 :               for (i=1; i<=n; i++)
     933           0 :                 gmael(B,kappa,i) = submuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
     934           0 :               for (i=1; i<=d; i++)
     935           0 :                 gmael(U,kappa,i) = submuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
     936             :             }
     937             :             else /* = -xx */
     938             :             {
     939           0 :               for (i=1; i<=n; i++)
     940           0 :                 gmael(B,kappa,i) = addmuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), -xx);
     941           0 :               for (i=1; i<=d; i++)
     942           0 :                 gmael(U,kappa,i) = addmuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), -xx);
     943             :             }
     944             :           } else
     945             :           {
     946     2792351 :             for (k=zeros+1; k<j; k++)
     947     1953628 :               dmael(mu,kappa,k) -= ldexp(((double)xx) * dmael(mu,j,k), E + e);
     948      838723 :             if (xx > 0) /* = xx */
     949             :             {
     950     3939334 :               for (i=1; i<=n; i++)
     951     3517304 :                 gmael(B,kappa,i) = submuliu642n(gmael(B,kappa,i), gmael(B,j,i), xx, E);
     952     1512926 :               for (i=1; i<=d; i++)
     953     1090896 :                 gmael(U,kappa,i) = submuliu642n(gmael(U,kappa,i), gmael(U,j,i), xx, E);
     954             :             }
     955             :             else /* = -xx */
     956             :             {
     957     3892499 :               for (i=1; i<=n; i++)
     958     3475820 :                 gmael(B,kappa,i) = addmuliu642n(gmael(B,kappa,i), gmael(B,j,i), -xx, E);
     959     1492466 :               for (i=1; i<=d; i++)
     960     1075787 :                 gmael(U,kappa,i) = addmuliu642n(gmael(U,kappa,i), gmael(U,j,i), -xx, E);
     961             :             }
     962             :           }
     963             :         }
     964             :       }
     965    47176107 :     if (!go_on) break; /* Anything happened? */
     966    17130413 :     expoB[kappa] = set_line(del(appB,kappa), gel(B,kappa), n);
     967    17133648 :     setG_fast(appB, n, G, kappa, zeros+1, kappa-1);
     968    17133791 :     aa = zeros+1;
     969             :   }
     970    30045694 :   if (did_something) setG2_fast(appB, n, G, kappa, kappa, maxG);
     971             : 
     972    30045919 :   del(s,zeros+1) = dmael(G,kappa,kappa);
     973             :   /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
     974   109478058 :   for (k=zeros+1; k<=kappa-2; k++)
     975    79432139 :     del(s,k+1) = del(s,k) - dmael(mu,kappa,k)*dmael(r,kappa,k);
     976    30045919 :   *pB = B; *pU = U; return 0;
     977             : }
     978             : 
     979             : static void
     980    11927542 : update_alpha(GEN alpha, long kappa, long kappa2, long kappamax)
     981             : {
     982             :   long i;
     983    37937120 :   for (i = kappa; i < kappa2; i++)
     984    26009578 :     if (kappa <= alpha[i]) alpha[i] = kappa;
     985    37937115 :   for (i = kappa2; i > kappa; i--) alpha[i] = alpha[i-1];
     986    23152724 :   for (i = kappa2+1; i <= kappamax; i++)
     987    11225182 :     if (kappa < alpha[i]) alpha[i] = kappa;
     988    11927542 :   alpha[kappa] = kappa;
     989    11927542 : }
     990             : static void
     991      440490 : rotateG(GEN G, long kappa2, long kappa, long maxG, GEN Gtmp)
     992             : {
     993             :   long i, j;
     994     3414087 :   for (i=1; i<=kappa2; i++) gel(Gtmp,i) = gmael(G,kappa2,i);
     995     1816012 :   for (   ; i<=maxG; i++)   gel(Gtmp,i) = gmael(G,i,kappa2);
     996     1543597 :   for (i=kappa2; i>kappa; i--)
     997             :     {
     998     5247270 :       for (j=1; j<kappa; j++) gmael(G,i,j) = gmael(G,i-1,j);
     999     1103107 :       gmael(G,i,kappa) = gel(Gtmp,i-1);
    1000     3977758 :       for (j=kappa+1; j<=i; j++) gmael(G,i,j) = gmael(G,i-1,j-1);
    1001     4787566 :       for (j=kappa2+1; j<=maxG; j++) gmael(G,j,i) = gmael(G,j,i-1);
    1002             :     }
    1003     1870490 :   for (i=1; i<kappa; i++) gmael(G,kappa,i) = gel(Gtmp,i);
    1004      440490 :   gmael(G,kappa,kappa) = gel(Gtmp,kappa2);
    1005     1816012 :   for (i=kappa2+1; i<=maxG; i++) gmael(G,i,kappa) = gel(Gtmp,i);
    1006      440490 : }
    1007             : static void
    1008    11486936 : rotateG_fast(double **G, long kappa2, long kappa, long maxG, double *Gtmp)
    1009             : {
    1010             :   long i, j;
    1011    66697825 :   for (i=1; i<=kappa2; i++) del(Gtmp,i) = dmael(G,kappa2,i);
    1012    22198966 :   for (   ; i<=maxG; i++) del(Gtmp,i) = dmael(G,i,kappa2);
    1013    36393209 :   for (i=kappa2; i>kappa; i--)
    1014             :   {
    1015    69091549 :     for (j=1; j<kappa; j++) dmael(G,i,j) = dmael(G,i-1,j);
    1016    24906273 :     dmael(G,i,kappa) = del(Gtmp,i-1);
    1017    84887536 :     for (j=kappa+1; j<=i; j++) dmael(G,i,j) = dmael(G,i-1,j-1);
    1018    46959245 :     for (j=kappa2+1; j<=maxG; j++) dmael(G,j,i) = dmael(G,j,i-1);
    1019             :   }
    1020    30304933 :   for (i=1; i<kappa; i++) dmael(G,kappa,i) = del(Gtmp,i);
    1021    11486936 :   dmael(G,kappa,kappa) = del(Gtmp,kappa2);
    1022    22198999 :   for (i=kappa2+1; i<=maxG; i++) dmael(G,i,kappa) = del(Gtmp,i);
    1023    11486936 : }
    1024             : 
    1025             : /* LLL-reduces (B,U) in place [apply base change transforms to B and U].
    1026             :  * Gram matrix, and GSO performed on matrices of 'double'.
    1027             :  * If (keepfirst), never swap with first vector.
    1028             :  * Return -1 on failure, else zeros = dim Kernel (>= 0) */
    1029             : static long
    1030     2069717 : fplll_fast(GEN *pB, GEN *pU, double delta, double eta, long keepfirst)
    1031             : {
    1032             :   pari_sp av;
    1033             :   long kappa, kappa2, d, n, i, j, zeros, kappamax, maxG;
    1034             :   double **mu, **r, *s, tmp, *Gtmp, **G, **appB;
    1035     2069717 :   GEN alpha, expoB, B = *pB, U;
    1036     2069717 :   long cnt = 0;
    1037             : 
    1038     2069717 :   d = lg(B)-1;
    1039     2069717 :   n = nbrows(B);
    1040     2069717 :   U = *pU; /* NULL if inplace */
    1041             : 
    1042     2069717 :   G = cget_dblmat(d+1);
    1043     2069714 :   appB = cget_dblmat(d+1);
    1044     2069712 :   mu = cget_dblmat(d+1);
    1045     2069713 :   r  = cget_dblmat(d+1);
    1046     2069713 :   s  = cget_dblvec(d+1);
    1047     9659966 :   for (j = 1; j <= d; j++)
    1048             :   {
    1049     7590244 :     del(mu,j) = cget_dblvec(d+1);
    1050     7590246 :     del(r,j) = cget_dblvec(d+1);
    1051     7590247 :     del(appB,j) = cget_dblvec(n+1);
    1052     7590245 :     del(G,j) = cget_dblvec(d+1);
    1053    46676239 :     for (i=1; i<=d; i++) dmael(G,j,i) = 0.;
    1054             :   }
    1055     2069722 :   expoB = cgetg(d+1, t_VECSMALL);
    1056     9659875 :   for (i=1; i<=d; i++) expoB[i] = set_line(del(appB,i), gel(B,i), n);
    1057     2069678 :   Gtmp = cget_dblvec(d+1);
    1058     2069716 :   alpha = cgetg(d+1, t_VECSMALL);
    1059     2069713 :   av = avma;
    1060             : 
    1061             :   /* Step2: Initializing the main loop */
    1062     2069713 :   kappamax = 1;
    1063     2069713 :   i = 1;
    1064     2069713 :   maxG = d; /* later updated to kappamax */
    1065             : 
    1066             :   do {
    1067     2235455 :     dmael(G,i,i) = dbldotsquare(del(appB,i),n);
    1068     2235457 :   } while (dmael(G,i,i) <= 0 && (++i <=d));
    1069     2069715 :   zeros = i-1; /* all vectors B[i] with i <= zeros are zero vectors */
    1070     2069715 :   kappa = i;
    1071     2069715 :   if (zeros < d) dmael(r,zeros+1,zeros+1) = dmael(G,zeros+1,zeros+1);
    1072     9494246 :   for (i=zeros+1; i<=d; i++) alpha[i]=1;
    1073    32115729 :   while (++kappa <= d)
    1074             :   {
    1075    30050139 :     if (kappa > kappamax)
    1076             :     {
    1077     5351263 :       if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
    1078     5351263 :       maxG = kappamax = kappa;
    1079     5351263 :       setG_fast(appB, n, G, kappa, zeros+1, kappa);
    1080             :     }
    1081             :     /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
    1082    30050131 :     if (Babai_fast(av, kappa, &B,&U, mu,r,s, appB, expoB, G, alpha[kappa],
    1083        4133 :                    zeros, maxG, eta)) { *pB=B; *pU=U; return -1; }
    1084             : 
    1085    30045892 :     tmp = ldexp(r[kappa-1][kappa-1] * delta, 2*(expoB[kappa-1]-expoB[kappa]));
    1086    30045892 :     if ((keepfirst && kappa == 2) || tmp <= del(s,kappa-1))
    1087             :     { /* Step4: Success of Lovasz's condition */
    1088    18558838 :       alpha[kappa] = kappa;
    1089    18558838 :       tmp = dmael(mu,kappa,kappa-1) * dmael(r,kappa,kappa-1);
    1090    18558838 :       dmael(r,kappa,kappa) = del(s,kappa-1)- tmp;
    1091    18558838 :       continue;
    1092             :     }
    1093             :     /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
    1094    11487054 :     if (DEBUGLEVEL>=4 && kappa==kappamax && del(s,kappa-1)!=0)
    1095           0 :       if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", 2*expoB[1] + dblexpo(del(s,1))); }
    1096    11487049 :     kappa2 = kappa;
    1097             :     do {
    1098    24906464 :       kappa--;
    1099    24906464 :       if (kappa<zeros+2 + (keepfirst ? 1: 0)) break;
    1100    18362856 :       tmp = dmael(r,kappa-1,kappa-1) * delta;
    1101    18362856 :       tmp = ldexp(tmp, 2*(expoB[kappa-1]-expoB[kappa2]));
    1102    18362856 :     } while (del(s,kappa-1) <= tmp);
    1103    11487049 :     update_alpha(alpha, kappa, kappa2, kappamax);
    1104             : 
    1105             :     /* Step6: Update the mu's and r's */
    1106    11487049 :     dblrotate(mu,kappa2,kappa);
    1107    11487027 :     dblrotate(r,kappa2,kappa);
    1108    11486998 :     dmael(r,kappa,kappa) = del(s,kappa);
    1109             : 
    1110             :     /* Step7: Update B, appB, U, G */
    1111    11486998 :     rotate(B,kappa2,kappa);
    1112    11486984 :     dblrotate(appB,kappa2,kappa);
    1113    11486963 :     if (U) rotate(U,kappa2,kappa);
    1114    11486960 :     rotate(expoB,kappa2,kappa);
    1115    11486932 :     rotateG_fast(G,kappa2,kappa, maxG, Gtmp);
    1116             : 
    1117             :     /* Step8: Prepare the next loop iteration */
    1118    11487177 :     if (kappa == zeros+1 && dmael(G,kappa,kappa)<= 0)
    1119             :     {
    1120      205570 :       zeros++; kappa++;
    1121      205570 :       dmael(G,kappa,kappa) = dbldotsquare(del(appB,kappa),n);
    1122      205569 :       dmael(r,kappa,kappa) = dmael(G,kappa,kappa);
    1123             :     }
    1124             :   }
    1125     2065590 :   *pB = B; *pU = U; return zeros;
    1126             : }
    1127             : 
    1128             : /***************** HEURISTIC version (reduced precision) ****************/
    1129             : static GEN
    1130      193281 : realsqrdotproduct(GEN x)
    1131             : {
    1132      193281 :   long i, l = lg(x);
    1133      193281 :   GEN z = sqrr(gel(x,1));
    1134     1254740 :   for (i=2; i<l; i++) z = addrr(z, sqrr(gel(x,i)));
    1135      193281 :   return z;
    1136             : }
    1137             : /* x, y non-empty vector of t_REALs, same length */
    1138             : static GEN
    1139     1167169 : realdotproduct(GEN x, GEN y)
    1140             : {
    1141             :   long i, l;
    1142             :   GEN z;
    1143     1167169 :   if (x == y) return realsqrdotproduct(x);
    1144      973888 :   l = lg(x); z = mulrr(gel(x,1),gel(y,1));
    1145     8478293 :   for (i=2; i<l; i++) z = addrr(z, mulrr(gel(x,i), gel(y,i)));
    1146      973888 :   return z;
    1147             : }
    1148             : static void
    1149      201454 : setG_heuristic(GEN appB, GEN G, long kappa, long a, long b)
    1150      201454 : { pari_sp av = avma;
    1151             :   long i;
    1152      919693 :   for (i = a; i <= b; i++)
    1153      718239 :     affrr(realdotproduct(gel(appB,kappa),gel(appB,i)), gmael(G,kappa,i));
    1154      201454 :   set_avma(av);
    1155      201454 : }
    1156             : static void
    1157      183605 : setG2_heuristic(GEN appB, GEN G, long kappa, long a, long b)
    1158      183605 : { pari_sp av = avma;
    1159             :   long i;
    1160      632535 :   for (i = a; i <= b; i++)
    1161      448930 :     affrr(realdotproduct(gel(appB,kappa),gel(appB,i)), gmael(G,i,kappa));
    1162      183605 :   set_avma(av);
    1163      183605 : }
    1164             : 
    1165             : /* approximate t_REAL x as m * 2^e, where |m| < 2^bit */
    1166             : static GEN
    1167       12132 : truncexpo(GEN x, long bit, long *e)
    1168             : {
    1169       12132 :   *e = expo(x) + 1 - bit;
    1170       12132 :   if (*e >= 0) return mantissa2nr(x, 0);
    1171         938 :   *e = 0; return roundr_safe(x);
    1172             : }
    1173             : /* Babai's Nearest Plane algorithm (iterative); see Babai() */
    1174             : static int
    1175      283467 : Babai_heuristic(pari_sp av, long kappa, GEN *pB, GEN *pU, GEN mu, GEN r, GEN s,
    1176             :                 GEN appB, GEN G, long a, long zeros, long maxG,
    1177             :                 GEN eta, long prec)
    1178             : {
    1179      283467 :   GEN B = *pB, U = *pU;
    1180      283467 :   const long n = nbrows(B), d = U ? lg(U)-1: 0, bit = prec2nbits(prec);
    1181      283467 :   long k, aa = (a > zeros)? a : zeros+1;
    1182      283467 :   int did_something = 0;
    1183      283467 :   long emaxmu = EX0, emax2mu = EX0;
    1184             :   /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
    1185             : 
    1186      191778 :   for (;;) {
    1187      475245 :     int go_on = 0;
    1188      475245 :     long i, j, emax3mu = emax2mu;
    1189             : 
    1190      475245 :     if (gc_needed(av,2))
    1191             :     {
    1192          27 :       if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
    1193          27 :       gc_lll(av,2,&B,&U);
    1194             :     }
    1195             :     /* Step2: compute the GSO for stage kappa */
    1196      475245 :     emax2mu = emaxmu; emaxmu = EX0;
    1197     1834675 :     for (j=aa; j<kappa; j++)
    1198             :     {
    1199     1359430 :       pari_sp btop = avma;
    1200     1359430 :       GEN g = gmael(G,kappa,j);
    1201     4429043 :       for (k = zeros+1; k<j; k++)
    1202     3069613 :         g = subrr(g, mulrr(gmael(mu,j,k), gmael(r,kappa,k)));
    1203     1359430 :       affrr(g, gmael(r,kappa,j));
    1204     1359430 :       affrr(divrr(gmael(r,kappa,j), gmael(r,j,j)), gmael(mu,kappa,j));
    1205     1359430 :       emaxmu = maxss(emaxmu, expo(gmael(mu,kappa,j)));
    1206     1359430 :       set_avma(btop);
    1207             :     }
    1208      475245 :     if (emax3mu != EX0 && emax3mu <= emax2mu + 5)
    1209        1328 :     { *pB = B; *pU = U; return 1; }
    1210             : 
    1211     1615560 :     for (j=kappa-1; j>zeros; j--)
    1212     1333421 :       if (abscmprr(gmael(mu,kappa,j), eta) > 0) { go_on = 1; break; }
    1213             : 
    1214             :     /* Step3--5: compute the X_j's  */
    1215      473917 :     if (go_on)
    1216      879315 :       for (j=kappa-1; j>zeros; j--)
    1217             :       { /* The code below seemingly handles U = NULL, but in this case d = 0 */
    1218             :         pari_sp btop;
    1219      687537 :         GEN tmp = gmael(mu,kappa,j);
    1220      687537 :         if (absrsmall(tmp)) continue; /* size-reduced */
    1221             : 
    1222      394782 :         if (gc_needed(av,2))
    1223             :         {
    1224          19 :           if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
    1225          19 :           gc_lll(av,2,&B,&U);
    1226             :         }
    1227      394782 :         btop = avma; did_something = 1;
    1228             :         /* we consider separately the case |X| = 1 */
    1229      394782 :         if (absrsmall2(tmp))
    1230             :         {
    1231      259542 :           if (signe(tmp) > 0) { /* in this case, X = 1 */
    1232      382670 :             for (k=zeros+1; k<j; k++)
    1233      253422 :               affrr(subrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
    1234      129248 :             set_avma(btop);
    1235     1229037 :             for (i=1; i<=n; i++)
    1236     1099789 :               gmael(B,kappa,i) = subii(gmael(B,kappa,i), gmael(B,j,i));
    1237      801586 :             for (i=1; i<=d; i++)
    1238      672338 :               gmael(U,kappa,i) = subii(gmael(U,kappa,i), gmael(U,j,i));
    1239             :           } else { /* otherwise X = -1 */
    1240      387593 :             for (k=zeros+1; k<j; k++)
    1241      257299 :               affrr(addrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
    1242      130294 :             set_avma(btop);
    1243     1243223 :             for (i=1; i<=n; i++)
    1244     1112929 :               gmael(B,kappa,i) = addii(gmael(B,kappa,i), gmael(B,j,i));
    1245      803235 :             for (i=1; i<=d; i++)
    1246      672941 :               gmael(U,kappa,i) = addii(gmael(U,kappa,i),gmael(U,j,i));
    1247             :           }
    1248      259542 :           continue;
    1249             :         }
    1250             :         /* we have |X| >= 2 */
    1251      135240 :         if (expo(tmp) < BITS_IN_LONG)
    1252             :         {
    1253      123108 :           ulong xx = roundr_safe(tmp)[2]; /* X fits in an ulong */
    1254      123108 :           if (signe(tmp) > 0) /* = xx */
    1255             :           {
    1256      135547 :             for (k=zeros+1; k<j; k++)
    1257       73418 :               affrr(subrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
    1258       73418 :                   gmael(mu,kappa,k));
    1259       62129 :             set_avma(btop);
    1260      413125 :             for (i=1; i<=n; i++)
    1261      350996 :               gmael(B,kappa,i) = submuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
    1262      304501 :             for (i=1; i<=d; i++)
    1263      242372 :               gmael(U,kappa,i) = submuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
    1264             :           }
    1265             :           else /* = -xx */
    1266             :           {
    1267      132116 :             for (k=zeros+1; k<j; k++)
    1268       71137 :               affrr(addrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
    1269       71137 :                   gmael(mu,kappa,k));
    1270       60979 :             set_avma(btop);
    1271      405589 :             for (i=1; i<=n; i++)
    1272      344610 :               gmael(B,kappa,i) = addmuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
    1273      290707 :             for (i=1; i<=d; i++)
    1274      229728 :               gmael(U,kappa,i) = addmuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
    1275             :           }
    1276             :         }
    1277             :         else
    1278             :         {
    1279             :           long e;
    1280       12132 :           GEN X = truncexpo(tmp, bit, &e); /* tmp ~ X * 2^e */
    1281       12132 :           btop = avma;
    1282       29542 :           for (k=zeros+1; k<j; k++)
    1283             :           {
    1284       17410 :             GEN x = mulir(X, gmael(mu,j,k));
    1285       17410 :             if (e) shiftr_inplace(x, e);
    1286       17410 :             affrr(subrr(gmael(mu,kappa,k), x), gmael(mu,kappa,k));
    1287             :           }
    1288       12132 :           set_avma(btop);
    1289       99229 :           for (i=1; i<=n; i++)
    1290       87097 :             gmael(B,kappa,i) = submulshift(gmael(B,kappa,i), gmael(B,j,i), X, e);
    1291       73211 :           for (i=1; i<=d; i++)
    1292       61079 :             gmael(U,kappa,i) = submulshift(gmael(U,kappa,i), gmael(U,j,i), X, e);
    1293             :         }
    1294             :       }
    1295      473917 :     if (!go_on) break; /* Anything happened? */
    1296     1440134 :     for (i=1 ; i<=n; i++) affir(gmael(B,kappa,i), gmael(appB,kappa,i));
    1297      191778 :     setG_heuristic(appB, G, kappa, zeros+1, kappa-1);
    1298      191778 :     aa = zeros+1;
    1299             :   }
    1300      282139 :   if (did_something) setG2_heuristic(appB, G, kappa, kappa, maxG);
    1301      282139 :   affrr(gmael(G,kappa,kappa), gel(s,zeros+1));
    1302             :   /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
    1303      282139 :   av = avma;
    1304     1013731 :   for (k=zeros+1; k<=kappa-2; k++)
    1305      731592 :     affrr(subrr(gel(s,k), mulrr(gmael(mu,kappa,k), gmael(r,kappa,k))),
    1306      731592 :           gel(s,k+1));
    1307      282139 :   *pB = B; *pU = U; return gc_bool(av, 0);
    1308             : }
    1309             : 
    1310             : static GEN
    1311       15707 : ZC_to_RC(GEN x, long prec)
    1312      103203 : { pari_APPLY_type(t_COL,itor(gel(x,i),prec)) }
    1313             : 
    1314             : static GEN
    1315        4133 : ZM_to_RM(GEN x, long prec)
    1316       19840 : { pari_APPLY_same(ZC_to_RC(gel(x,i),prec)) }
    1317             : 
    1318             : /* LLL-reduces (B,U) in place [apply base change transforms to B and U].
    1319             :  * Gram matrix made of t_REAL at precision prec2, performe GSO at prec.
    1320             :  * If (keepfirst), never swap with first vector.
    1321             :  * Return -1 on failure, else zeros = dim Kernel (>= 0) */
    1322             : static long
    1323        4133 : fplll_heuristic(GEN *pB, GEN *pU, double DELTA, double ETA, long keepfirst,
    1324             :                 long prec, long prec2)
    1325             : {
    1326             :   pari_sp av, av2;
    1327             :   long kappa, kappa2, d, i, j, zeros, kappamax, maxG;
    1328        4133 :   GEN mu, r, s, tmp, Gtmp, alpha, G, appB, B = *pB, U;
    1329        4133 :   GEN delta = dbltor(DELTA), eta = dbltor(ETA);
    1330        4133 :   long cnt = 0;
    1331             : 
    1332        4133 :   d = lg(B)-1;
    1333        4133 :   U = *pU; /* NULL if inplace */
    1334             : 
    1335        4133 :   G = cgetg(d+1, t_MAT);
    1336        4133 :   mu = cgetg(d+1, t_MAT);
    1337        4133 :   r  = cgetg(d+1, t_MAT);
    1338        4133 :   s  = cgetg(d+1, t_VEC);
    1339        4133 :   appB = ZM_to_RM(B, prec2);
    1340       19840 :   for (j = 1; j <= d; j++)
    1341             :   {
    1342       15707 :     GEN M = cgetg(d+1, t_COL), R = cgetg(d+1, t_COL), S = cgetg(d+1, t_COL);
    1343       15707 :     gel(mu,j)= M;
    1344       15707 :     gel(r,j) = R;
    1345       15707 :     gel(G,j) = S;
    1346       15707 :     gel(s,j) = cgetr(prec);
    1347       94954 :     for (i = 1; i <= d; i++)
    1348             :     {
    1349       79247 :       gel(R,i) = cgetr(prec);
    1350       79247 :       gel(M,i) = cgetr(prec);
    1351       79247 :       gel(S,i) = cgetr(prec2);
    1352             :     }
    1353             :   }
    1354        4133 :   Gtmp = cgetg(d+1, t_VEC);
    1355        4133 :   alpha = cgetg(d+1, t_VECSMALL);
    1356        4133 :   av = avma;
    1357             : 
    1358             :   /* Step2: Initializing the main loop */
    1359        4133 :   kappamax = 1;
    1360        4133 :   i = 1;
    1361        4133 :   maxG = d; /* later updated to kappamax */
    1362             : 
    1363             :   do {
    1364        4136 :     affrr(RgV_dotsquare(gel(appB,i)), gmael(G,i,i));
    1365        4136 :   } while (signe(gmael(G,i,i)) == 0 && (++i <=d));
    1366        4133 :   zeros = i-1; /* all vectors B[i] with i <= zeros are zero vectors */
    1367        4133 :   kappa = i;
    1368        4133 :   if (zeros < d) affrr(gmael(G,zeros+1,zeros+1), gmael(r,zeros+1,zeros+1));
    1369       19837 :   for (i=zeros+1; i<=d; i++) alpha[i]=1;
    1370             : 
    1371      286272 :   while (++kappa <= d)
    1372             :   {
    1373      283467 :     if (kappa > kappamax)
    1374             :     {
    1375        9676 :       if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
    1376        9676 :       maxG = kappamax = kappa;
    1377        9676 :       setG_heuristic(appB, G, kappa, zeros+1, kappa);
    1378             :     }
    1379             :     /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
    1380      283467 :     if (Babai_heuristic(av, kappa, &B,&U, mu,r,s, appB, G, alpha[kappa], zeros,
    1381        1328 :                         maxG, eta, prec)) { *pB = B; *pU = U; return -1; }
    1382      282139 :     av2 = avma;
    1383      564171 :     if ((keepfirst && kappa == 2) ||
    1384      282032 :         cmprr(mulrr(gmael(r,kappa-1,kappa-1), delta), gel(s,kappa-1)) <= 0)
    1385             :     { /* Step4: Success of Lovasz's condition */
    1386      167391 :       alpha[kappa] = kappa;
    1387      167391 :       tmp = mulrr(gmael(mu,kappa,kappa-1), gmael(r,kappa,kappa-1));
    1388      167391 :       affrr(subrr(gel(s,kappa-1), tmp), gmael(r,kappa,kappa));
    1389      167391 :       set_avma(av2); continue;
    1390             :     }
    1391             :     /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
    1392      114748 :     if (DEBUGLEVEL>=4 && kappa==kappamax && signe(gel(s,kappa-1)))
    1393           0 :       if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", expo(gel(s,1))); }
    1394      114748 :     kappa2 = kappa;
    1395             :     do {
    1396      273819 :       kappa--;
    1397      273819 :       if (kappa < zeros+2 + (keepfirst ? 1: 0)) break;
    1398      244891 :       tmp = mulrr(gmael(r,kappa-1,kappa-1), delta);
    1399      244891 :     } while (cmprr(gel(s,kappa-1), tmp) <= 0 );
    1400      114748 :     set_avma(av2);
    1401      114748 :     update_alpha(alpha, kappa, kappa2, kappamax);
    1402             : 
    1403             :     /* Step6: Update the mu's and r's */
    1404      114748 :     rotate(mu,kappa2,kappa);
    1405      114748 :     rotate(r,kappa2,kappa);
    1406      114748 :     affrr(gel(s,kappa), gmael(r,kappa,kappa));
    1407             : 
    1408             :     /* Step7: Update B, appB, U, G */
    1409      114748 :     rotate(B,kappa2,kappa);
    1410      114748 :     rotate(appB,kappa2,kappa);
    1411      114748 :     if (U) rotate(U,kappa2,kappa);
    1412      114748 :     rotateG(G,kappa2,kappa, maxG, Gtmp);
    1413             : 
    1414             :     /* Step8: Prepare the next loop iteration */
    1415      114748 :     if (kappa == zeros+1 && !signe(gmael(G,kappa,kappa)))
    1416             :     {
    1417           0 :       zeros++; kappa++;
    1418           0 :       affrr(RgV_dotsquare(gel(appB,kappa)), gmael(G,kappa,kappa));
    1419           0 :       affrr(gmael(G,kappa,kappa), gmael(r,kappa,kappa));
    1420             :     }
    1421             :   }
    1422        2805 :   *pB=B; *pU=U; return zeros;
    1423             : }
    1424             : 
    1425             : /************************* PROVED version (t_INT) ***********************/
    1426             : /* dpe inspired by dpe.h by Patrick Pelissier, Paul Zimmermann
    1427             :  * https://gforge.inria.fr/projects/dpe/
    1428             :  */
    1429             : 
    1430             : typedef struct
    1431             : {
    1432             :   double d;  /* significand */
    1433             :   long e; /* exponent */
    1434             : } dpe_t;
    1435             : 
    1436             : #define Dmael(x,i,j) (&((x)[i][j]))
    1437             : #define Del(x,i) (&((x)[i]))
    1438             : 
    1439             : static void
    1440      651484 : dperotate(dpe_t **A, long k2, long k)
    1441             : {
    1442             :   long i;
    1443      651484 :   dpe_t *B = A[k2];
    1444     2310060 :   for (i = k2; i > k; i--) A[i] = A[i-1];
    1445      651484 :   A[k] = B;
    1446      651484 : }
    1447             : 
    1448             : static void
    1449   108048399 : dpe_normalize0(dpe_t *x)
    1450             : {
    1451             :   int e;
    1452   108048399 :   x->d = frexp(x->d, &e);
    1453   108048399 :   x->e += e;
    1454   108048399 : }
    1455             : 
    1456             : static void
    1457    47910070 : dpe_normalize(dpe_t *x)
    1458             : {
    1459    47910070 :   if (x->d == 0.0)
    1460      496437 :     x->e = -LONG_MAX;
    1461             :   else
    1462    47413633 :     dpe_normalize0(x);
    1463    47910120 : }
    1464             : 
    1465             : static GEN
    1466       94015 : dpetor(dpe_t *x)
    1467             : {
    1468       94015 :   GEN r = dbltor(x->d);
    1469       94015 :   if (signe(r)==0) return r;
    1470       94015 :   setexpo(r, x->e-1);
    1471       94015 :   return r;
    1472             : }
    1473             : 
    1474             : static void
    1475    25570064 : affdpe(dpe_t *y, dpe_t *x)
    1476             : {
    1477    25570064 :   x->d = y->d;
    1478    25570064 :   x->e = y->e;
    1479    25570064 : }
    1480             : 
    1481             : static void
    1482    20527375 : affidpe(GEN y, dpe_t *x)
    1483             : {
    1484    20527375 :   pari_sp av = avma;
    1485    20527375 :   GEN r = itor(y, DEFAULTPREC);
    1486    20527089 :   x->e = expo(r)+1;
    1487    20527089 :   setexpo(r,-1);
    1488    20527027 :   x->d = rtodbl(r);
    1489    20526967 :   set_avma(av);
    1490    20526919 : }
    1491             : 
    1492             : static void
    1493     3137026 : affdbldpe(double y, dpe_t *x)
    1494             : {
    1495     3137026 :   x->d = (double)y;
    1496     3137026 :   x->e = 0;
    1497     3137026 :   dpe_normalize(x);
    1498     3137026 : }
    1499             : 
    1500             : static void
    1501    56725292 : dpe_mulz(dpe_t *x, dpe_t *y, dpe_t *z)
    1502             : {
    1503    56725292 :   z->d = x->d * y->d;
    1504    56725292 :   if (z->d == 0.0)
    1505     8148502 :     z->e = -LONG_MAX;
    1506             :   else
    1507             :   {
    1508    48576790 :     z->e = x->e + y->e;
    1509    48576790 :     dpe_normalize0(z);
    1510             :   }
    1511    56725546 : }
    1512             : 
    1513             : static void
    1514    13959826 : dpe_divz(dpe_t *x, dpe_t *y, dpe_t *z)
    1515             : {
    1516    13959826 :   z->d = x->d / y->d;
    1517    13959826 :   if (z->d == 0.0)
    1518     1901571 :     z->e = -LONG_MAX;
    1519             :   else
    1520             :   {
    1521    12058255 :     z->e = x->e - y->e;
    1522    12058255 :     dpe_normalize0(z);
    1523             :   }
    1524    13959880 : }
    1525             : 
    1526             : static void
    1527      244090 : dpe_negz(dpe_t *y, dpe_t *x)
    1528             : {
    1529      244090 :   x->d = - y->d;
    1530      244090 :   x->e = y->e;
    1531      244090 : }
    1532             : 
    1533             : static void
    1534     1942243 : dpe_addz(dpe_t *y, dpe_t *z, dpe_t *x)
    1535             : {
    1536     1942243 :   if (y->e > z->e + 53)
    1537      112062 :     affdpe(y, x);
    1538     1830181 :   else if (z->e > y->e + 53)
    1539       41675 :     affdpe(z, x);
    1540             :   else
    1541             :   {
    1542     1788506 :     long d = y->e - z->e;
    1543             : 
    1544     1788506 :     if (d >= 0)
    1545             :     {
    1546     1344379 :       x->d = y->d + ldexp(z->d, -d);
    1547     1344379 :       x->e  = y->e;
    1548             :     }
    1549             :     else
    1550             :     {
    1551      444127 :       x->d = z->d + ldexp(y->d, d);
    1552      444127 :       x->e = z->e;
    1553             :     }
    1554     1788506 :     dpe_normalize(x);
    1555             :   }
    1556     1942245 : }
    1557             : static void
    1558    53560169 : dpe_subz(dpe_t *y, dpe_t *z, dpe_t *x)
    1559             : {
    1560    53560169 :   if (y->e > z->e + 53)
    1561    11130894 :     affdpe(y, x);
    1562    42429275 :   else if (z->e > y->e + 53)
    1563      244090 :     dpe_negz(z, x);
    1564             :   else
    1565             :   {
    1566    42185185 :     long d = y->e - z->e;
    1567             : 
    1568    42185185 :     if (d >= 0)
    1569             :     {
    1570    39424565 :       x->d = y->d - ldexp(z->d, -d);
    1571    39424565 :       x->e = y->e;
    1572             :     }
    1573             :     else
    1574             :     {
    1575     2760620 :       x->d = ldexp(y->d, d) - z->d;
    1576     2760620 :       x->e = z->e;
    1577             :     }
    1578    42185185 :     dpe_normalize(x);
    1579             :   }
    1580    53560394 : }
    1581             : 
    1582             : static void
    1583      799191 : dpe_muluz(dpe_t *y, ulong t, dpe_t *x)
    1584             : {
    1585      799191 :   x->d = y->d * (double)t;
    1586      799191 :   x->e = y->e;
    1587      799191 :   dpe_normalize(x);
    1588      799191 : }
    1589             : 
    1590             : static void
    1591      342689 : dpe_addmuluz(dpe_t *y,  dpe_t *z, ulong t, dpe_t *x)
    1592             : {
    1593             :   dpe_t tmp;
    1594      342689 :   dpe_muluz(z, t, &tmp);
    1595      342689 :   dpe_addz(y, &tmp, x);
    1596      342689 : }
    1597             : 
    1598             : static void
    1599      411952 : dpe_submuluz(dpe_t *y,  dpe_t *z, ulong t, dpe_t *x)
    1600             : {
    1601             :   dpe_t tmp;
    1602      411952 :   dpe_muluz(z, t, &tmp);
    1603      411952 :   dpe_subz(y, &tmp, x);
    1604      411952 : }
    1605             : 
    1606             : static void
    1607    51529917 : dpe_submulz(dpe_t *y,  dpe_t *z, dpe_t *t, dpe_t *x)
    1608             : {
    1609             :   dpe_t tmp;
    1610    51529917 :   dpe_mulz(z, t, &tmp);
    1611    51529958 :   dpe_subz(y, &tmp, x);
    1612    51530005 : }
    1613             : 
    1614             : static int
    1615     5195500 : dpe_cmp(dpe_t *x, dpe_t *y)
    1616             : {
    1617     5195500 :   int sx = x->d < 0. ? -1: x->d > 0.;
    1618     5195500 :   int sy = y->d < 0. ? -1: y->d > 0.;
    1619     5195500 :   int d  = sx - sy;
    1620             : 
    1621     5195500 :   if (d != 0)
    1622      141626 :     return d;
    1623     5053874 :   else if (x->e > y->e)
    1624      481189 :     return (sx > 0) ? 1 : -1;
    1625     4572685 :   else if (y->e > x->e)
    1626     2502201 :     return (sx > 0) ? -1 : 1;
    1627             :   else
    1628     2070484 :     return (x->d < y->d) ? -1 : (x->d > y->d);
    1629             : }
    1630             : 
    1631             : static int
    1632    14400508 : dpe_abscmp(dpe_t *x, dpe_t *y)
    1633             : {
    1634    14400508 :   if (x->e > y->e)
    1635      271246 :     return 1;
    1636    14129262 :   else if (y->e > x->e)
    1637    13278345 :     return -1;
    1638             :   else
    1639      850917 :     return (fabs(x->d) < fabs(y->d)) ? -1 : (fabs(x->d) > fabs(y->d));
    1640             : }
    1641             : 
    1642             : static int
    1643     1388062 : dpe_abssmall(dpe_t *x)
    1644             : {
    1645     1388062 :   return (x->e <= 0) || (x->e == 1 && fabs(x->d) <= .75);
    1646             : }
    1647             : 
    1648             : static int
    1649     5195493 : dpe_cmpmul(dpe_t *x, dpe_t *y, dpe_t *z)
    1650             : {
    1651             :   dpe_t t;
    1652     5195493 :   dpe_mulz(x,y,&t);
    1653     5195496 :   return dpe_cmp(&t, z);
    1654             : }
    1655             : 
    1656             : static dpe_t *
    1657    13045610 : cget_dpevec(long d)
    1658    13045610 : { return (dpe_t*) stack_malloc_align(d*sizeof(dpe_t), sizeof(dpe_t)); }
    1659             : 
    1660             : static dpe_t **
    1661     3137030 : cget_dpemat(long d) { return (dpe_t **) cgetg(d, t_VECSMALL); }
    1662             : 
    1663             : static GEN
    1664       20292 : dpeM_diagonal_shallow(dpe_t **m, long d)
    1665             : {
    1666             :   long i;
    1667       20292 :   GEN y = cgetg(d+1,t_VEC);
    1668      114307 :   for (i=1; i<=d; i++) gel(y, i) = dpetor(Dmael(m,i,i));
    1669       20292 :   return y;
    1670             : }
    1671             : 
    1672             : static void
    1673     1388053 : affii_or_copy_gc(pari_sp av, GEN x, GEN *y)
    1674             : {
    1675     1388053 :   long l = lg(*y);
    1676     1388053 :   if (lgefint(x) <= l && isonstack(*y))
    1677             :   {
    1678     1388039 :     affii(x,*y);
    1679     1388039 :     set_avma(av);
    1680             :   }
    1681             :   else
    1682          12 :     *y = gerepileuptoint(av, x);
    1683     1388057 : }
    1684             : 
    1685             : /* *x -= u*y */
    1686             : INLINE void
    1687     5920105 : submulziu(GEN *x, GEN y, ulong u)
    1688             : {
    1689             :   pari_sp av;
    1690     5920105 :   long ly = lgefint(y);
    1691     5920105 :   if (ly == 2) return;
    1692     3251980 :   av = avma;
    1693     3251980 :   (void)new_chunk(3+ly+lgefint(*x)); /* HACK */
    1694     3252008 :   y = mului(u,y);
    1695     3251996 :   set_avma(av); subzi(x, y);
    1696             : }
    1697             : 
    1698             : /* *x += u*y */
    1699             : INLINE void
    1700     4577417 : addmulziu(GEN *x, GEN y, ulong u)
    1701             : {
    1702             :   pari_sp av;
    1703     4577417 :   long ly = lgefint(y);
    1704     4577417 :   if (ly == 2) return;
    1705     2755761 :   av = avma;
    1706     2755761 :   (void)new_chunk(3+ly+lgefint(*x)); /* HACK */
    1707     2755789 :   y = mului(u,y);
    1708     2755784 :   set_avma(av); addzi(x, y);
    1709             : }
    1710             : 
    1711             : /************************** PROVED version (dpe) *************************/
    1712             : 
    1713             : /* Babai's Nearest Plane algorithm (iterative).
    1714             :  * Size-reduces b_kappa using mu_{i,j} and r_{i,j} for j<=i <kappa
    1715             :  * Update B[,kappa]; compute mu_{kappa,j}, r_{kappa,j} for j<=kappa and s[kappa]
    1716             :  * mu, r, s updated in place (affrr). Return 1 on failure, else 0. */
    1717             : static int
    1718     4587532 : Babai_dpe(pari_sp av, long kappa, GEN *pG, GEN *pB, GEN *pU, dpe_t **mu, dpe_t **r, dpe_t *s,
    1719             :       long a, long zeros, long maxG, dpe_t *eta)
    1720             : {
    1721     4587532 :   GEN G = *pG, B = *pB, U = *pU, ztmp;
    1722     4587532 :   long k, d, n, aa = a > zeros? a: zeros+1;
    1723     4587532 :   long emaxmu = EX0, emax2mu = EX0;
    1724             :   /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
    1725     4587532 :   d = U? lg(U)-1: 0;
    1726     4587532 :   n = B? nbrows(B): 0;
    1727      522328 :   for (;;) {
    1728     5109878 :     int go_on = 0;
    1729     5109878 :     long i, j, emax3mu = emax2mu;
    1730             : 
    1731     5109878 :     if (gc_needed(av,2))
    1732             :     {
    1733           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
    1734           0 :       gc_lll(av,3,&G,&B,&U);
    1735             :     }
    1736             :     /* Step2: compute the GSO for stage kappa */
    1737     5109877 :     emax2mu = emaxmu; emaxmu = EX0;
    1738    19069756 :     for (j=aa; j<kappa; j++)
    1739             :     {
    1740             :       dpe_t g;
    1741    13959886 :       affidpe(gmael(G,kappa,j), &g);
    1742    52342336 :       for (k = zeros+1; k < j; k++)
    1743    38382552 :         dpe_submulz(&g, Dmael(mu,j,k), Dmael(r,kappa,k), &g);
    1744    13959784 :       affdpe(&g, Dmael(r,kappa,j));
    1745    13959825 :       dpe_divz(Dmael(r,kappa,j), Dmael(r,j,j), Dmael(mu,kappa,j));
    1746    13959832 :       emaxmu = maxss(emaxmu, Dmael(mu,kappa,j)->e);
    1747             :     }
    1748     5109870 :     if (emax3mu != EX0 && emax3mu <= emax2mu + 5) /* precision too low */
    1749           0 :     { *pG = G; *pB = B; *pU = U; return 1; }
    1750             : 
    1751    18988033 :     for (j=kappa-1; j>zeros; j--)
    1752    14400508 :       if (dpe_abscmp(Dmael(mu,kappa,j), eta) > 0) { go_on = 1; break; }
    1753             : 
    1754             :     /* Step3--5: compute the X_j's  */
    1755     5109858 :     if (go_on)
    1756     3031509 :       for (j=kappa-1; j>zeros; j--)
    1757             :       {
    1758             :         pari_sp btop;
    1759     2509169 :         dpe_t *tmp = Dmael(mu,kappa,j);
    1760     2509169 :         if (tmp->e < 0) continue; /* (essentially) size-reduced */
    1761             : 
    1762     1388062 :         if (gc_needed(av,2))
    1763             :         {
    1764           0 :           if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
    1765           0 :           gc_lll(av,3,&G,&B,&U);
    1766             :         }
    1767             :         /* we consider separately the case |X| = 1 */
    1768     1388062 :         if (dpe_abssmall(tmp))
    1769             :         {
    1770      920842 :           if (tmp->d > 0) { /* in this case, X = 1 */
    1771     2057991 :             for (k=zeros+1; k<j; k++)
    1772     1596339 :               dpe_subz(Dmael(mu,kappa,k), Dmael(mu,j,k), Dmael(mu,kappa,k));
    1773     3017484 :             for (i=1; i<=n; i++)
    1774     2555832 :               subzi(&gmael(B,kappa,i), gmael(B,j,i));
    1775     6971493 :             for (i=1; i<=d; i++)
    1776     6509848 :               subzi(&gmael(U,kappa,i), gmael(U,j,i));
    1777      461645 :             btop = avma;
    1778      461645 :             ztmp = subii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
    1779      461652 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    1780      461652 :             affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
    1781     2859843 :             for (i=1; i<=j; i++)
    1782     2398192 :               subzi(&gmael(G,kappa,i), gmael(G,j,i));
    1783     2191780 :             for (i=j+1; i<kappa; i++)
    1784     1730129 :               subzi(&gmael(G,kappa,i), gmael(G,i,j));
    1785     2361183 :             for (i=kappa+1; i<=maxG; i++)
    1786     1899531 :               subzi(&gmael(G,i,kappa), gmael(G,i,j));
    1787             :           } else { /* otherwise X = -1 */
    1788     2036131 :             for (k=zeros+1; k<j; k++)
    1789     1576941 :               dpe_addz(Dmael(mu,kappa,k), Dmael(mu,j,k), Dmael(mu,kappa,k));
    1790     3012875 :             for (i=1; i<=n; i++)
    1791     2553685 :               addzi(&gmael(B,kappa,i),gmael(B,j,i));
    1792     6844739 :             for (i=1; i<=d; i++)
    1793     6385545 :               addzi(&gmael(U,kappa,i),gmael(U,j,i));
    1794      459194 :             btop = avma;
    1795      459194 :             ztmp = addii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
    1796      459188 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    1797      459189 :             affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
    1798     2770700 :             for (i=1; i<=j; i++)
    1799     2311511 :               addzi(&gmael(G,kappa,i), gmael(G,j,i));
    1800     2195763 :             for (i=j+1; i<kappa; i++)
    1801     1736573 :               addzi(&gmael(G,kappa,i), gmael(G,i,j));
    1802     2313907 :             for (i=kappa+1; i<=maxG; i++)
    1803     1854718 :               addzi(&gmael(G,i,kappa), gmael(G,i,j));
    1804             :           }
    1805      920841 :           continue;
    1806             :         }
    1807             :         /* we have |X| >= 2 */
    1808      467221 :         if (tmp->e < BITS_IN_LONG-1)
    1809             :         {
    1810      448218 :           if (tmp->d > 0)
    1811             :           {
    1812      247218 :             ulong xx = (ulong) pari_rint(ldexp(tmp->d, tmp->e)); /* X fits in an ulong */
    1813      659170 :             for (k=zeros+1; k<j; k++)
    1814      411952 :               dpe_submuluz(Dmael(mu,kappa,k), Dmael(mu,j,k), xx, Dmael(mu,kappa,k));
    1815      722192 :             for (i=1; i<=n; i++)
    1816      474974 :               submulziu(&gmael(B,kappa,i), gmael(B,j,i), xx);
    1817     3107448 :             for (i=1; i<=d; i++)
    1818     2860229 :               submulziu(&gmael(U,kappa,i), gmael(U,j,i), xx);
    1819      247219 :             btop = avma;
    1820      247219 :             ztmp = submuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
    1821      247217 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    1822      247218 :             affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
    1823     1313885 :             for (i=1; i<=j; i++)
    1824     1066668 :               submulziu(&gmael(G,kappa,i), gmael(G,j,i), xx);
    1825      807824 :             for (i=j+1; i<kappa; i++)
    1826      560607 :               submulziu(&gmael(G,kappa,i), gmael(G,i,j), xx);
    1827     1204885 :             for (i=kappa+1; i<=maxG; i++)
    1828      957668 :               submulziu(&gmael(G,i,kappa), gmael(G,i,j), xx);
    1829             :           }
    1830             :           else
    1831             :           {
    1832      201000 :             ulong xx = (ulong) pari_rint(ldexp(-tmp->d, tmp->e)); /* X fits in an ulong */
    1833      543689 :             for (k=zeros+1; k<j; k++)
    1834      342689 :               dpe_addmuluz(Dmael(mu,kappa,k), Dmael(mu,j,k), xx, Dmael(mu,kappa,k));
    1835      687646 :             for (i=1; i<=n; i++)
    1836      486646 :               addmulziu(&gmael(B,kappa,i), gmael(B,j,i), xx);
    1837     2359434 :             for (i=1; i<=d; i++)
    1838     2158434 :               addmulziu(&gmael(U,kappa,i), gmael(U,j,i), xx);
    1839      201000 :             btop = avma;
    1840      201000 :             ztmp = addmuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
    1841      201000 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    1842      201000 :             affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
    1843      990085 :             for (i=1; i<=j; i++)
    1844      789085 :               addmulziu(&gmael(G,kappa,i), gmael(G,j,i), xx);
    1845      662229 :             for (i=j+1; i<kappa; i++)
    1846      461230 :               addmulziu(&gmael(G,kappa,i), gmael(G,i,j), xx);
    1847      883057 :             for (i=kappa+1; i<=maxG; i++)
    1848      682057 :               addmulziu(&gmael(G,i,kappa), gmael(G,i,j), xx);
    1849             :           }
    1850             :         }
    1851             :         else
    1852             :         {
    1853       19003 :           long e = tmp->e - BITS_IN_LONG + 1;
    1854       19003 :           if (tmp->d > 0)
    1855             :           {
    1856        9381 :             ulong xx = (ulong) pari_rint(ldexp(tmp->d, BITS_IN_LONG - 1));
    1857       31318 :             for (k=zeros+1; k<j; k++)
    1858             :             {
    1859             :               dpe_t x;
    1860       21937 :               dpe_muluz(Dmael(mu,j,k), xx, &x);
    1861       21937 :               x.e += e;
    1862       21937 :               dpe_subz(Dmael(mu,kappa,k), &x, Dmael(mu,kappa,k));
    1863             :             }
    1864      124141 :             for (i=1; i<=n; i++)
    1865      114760 :               submulzu2n(&gmael(B,kappa,i), gmael(B,j,i), xx, e);
    1866       86854 :             for (i=1; i<=d; i++)
    1867       77473 :               submulzu2n(&gmael(U,kappa,i), gmael(U,j,i), xx, e);
    1868        9381 :             btop = avma;
    1869        9381 :             ztmp = submuliu2n(mulshift(gmael(G,j,j), sqru(xx), 2*e),
    1870        9381 :                 gmael(G,kappa,j), xx, e+1);
    1871        9381 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    1872        9381 :             affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
    1873       40897 :             for (i=1; i<=j; i++)
    1874       31516 :               submulzu2n(&gmael(G,kappa,i), gmael(G,j,i), xx, e);
    1875       47307 :             for (   ; i<kappa; i++)
    1876       37926 :               submulzu2n(&gmael(G,kappa,i), gmael(G,i,j), xx, e);
    1877       10255 :             for (i=kappa+1; i<=maxG; i++)
    1878         874 :               submulzu2n(&gmael(G,i,kappa), gmael(G,i,j), xx, e);
    1879             :           } else
    1880             :           {
    1881        9622 :             ulong xx = (ulong) pari_rint(ldexp(-tmp->d, BITS_IN_LONG - 1));
    1882       32244 :             for (k=zeros+1; k<j; k++)
    1883             :             {
    1884             :               dpe_t x;
    1885       22614 :               dpe_muluz(Dmael(mu,j,k), xx, &x);
    1886       22614 :               x.e += e;
    1887       22614 :               dpe_addz(Dmael(mu,kappa,k), &x, Dmael(mu,kappa,k));
    1888             :             }
    1889      128739 :             for (i=1; i<=n; i++)
    1890      119109 :               addmulzu2n(&gmael(B,kappa,i), gmael(B,j,i), xx, e);
    1891       89526 :             for (i=1; i<=d; i++)
    1892       79896 :               addmulzu2n(&gmael(U,kappa,i), gmael(U,j,i), xx, e);
    1893        9630 :             btop = avma;
    1894        9630 :             ztmp = addmuliu2n(mulshift(gmael(G,j,j), sqru(xx), 2*e),
    1895        9630 :                 gmael(G,kappa,j), xx, e+1);
    1896        9630 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    1897        9630 :             affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
    1898       42072 :             for (i=1; i<=j; i++)
    1899       32442 :               addmulzu2n(&gmael(G,kappa,i), gmael(G,j,i), xx, e);
    1900       48988 :             for (   ; i<kappa; i++)
    1901       39358 :               addmulzu2n(&gmael(G,kappa,i), gmael(G,i,j), xx, e);
    1902       10388 :             for (i=kappa+1; i<=maxG; i++)
    1903         758 :               addmulzu2n(&gmael(G,i,kappa), gmael(G,i,j), xx, e);
    1904             :           }
    1905             :         }
    1906             :       }
    1907     5109865 :     if (!go_on) break; /* Anything happened? */
    1908      522328 :     aa = zeros+1;
    1909             :   }
    1910             : 
    1911     4587537 :   affidpe(gmael(G,kappa,kappa), Del(s,zeros+1));
    1912             :   /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
    1913    13473319 :   for (k=zeros+1; k<=kappa-2; k++)
    1914     8885798 :     dpe_submulz(Del(s,k), Dmael(mu,kappa,k), Dmael(r,kappa,k), Del(s,k+1));
    1915     4587521 :   *pG = G; *pB = B; *pU = U; return 0;
    1916             : }
    1917             : 
    1918             : /* G integral Gram matrix, LLL-reduces (G,B,U) in place [apply base change
    1919             :  * transforms to B and U]. If (keepfirst), never swap with first vector.
    1920             :  * If G = NULL, we compute the Gram matrix incrementally.
    1921             :  * Return -1 on failure, else zeros = dim Kernel (>= 0) */
    1922             : static long
    1923     1568515 : fplll_dpe(GEN *pG, GEN *pB, GEN *pU, GEN *pr, double DELTA, double ETA,
    1924             :       long keepfirst)
    1925             : {
    1926             :   pari_sp av;
    1927     1568515 :   GEN Gtmp, alpha, G = *pG, B = *pB, U = *pU;
    1928     1568515 :   long d, maxG, kappa, kappa2, i, j, zeros, kappamax, incgram = !G, cnt = 0;
    1929             :   dpe_t delta, eta, **mu, **r, *s;
    1930     1568515 :   affdbldpe(DELTA,&delta);
    1931     1568513 :   affdbldpe(ETA,&eta);
    1932             : 
    1933     1568514 :   if (incgram)
    1934             :   { /* incremental Gram matrix */
    1935     1508084 :     maxG = 2; d = lg(B)-1;
    1936     1508084 :     G = zeromatcopy(d, d);
    1937             :   }
    1938             :   else
    1939       60430 :     maxG = d = lg(G)-1;
    1940             : 
    1941     1568518 :   mu = cget_dpemat(d+1);
    1942     1568515 :   r  = cget_dpemat(d+1);
    1943     1568511 :   s  = cget_dpevec(d+1);
    1944     7307087 :   for (j = 1; j <= d; j++)
    1945             :   {
    1946     5738575 :     mu[j]= cget_dpevec(d+1);
    1947     5738571 :     r[j] = cget_dpevec(d+1);
    1948             :   }
    1949     1568512 :   Gtmp = cgetg(d+1, t_VEC);
    1950     1568514 :   alpha = cgetg(d+1, t_VECSMALL);
    1951     1568509 :   av = avma;
    1952             : 
    1953             :   /* Step2: Initializing the main loop */
    1954     1568509 :   kappamax = 1;
    1955     1568509 :   i = 1;
    1956             :   do {
    1957     1945102 :     if (incgram) gmael(G,i,i) = ZV_dotsquare(gel(B,i));
    1958     1945034 :     affidpe(gmael(G,i,i), Dmael(r,i,i));
    1959     1945061 :   } while (!signe(gmael(G,i,i)) && ++i <= d);
    1960     1568468 :   zeros = i-1; /* all basis vectors b_i with i <= zeros are zero vectors */
    1961     1568468 :   kappa = i;
    1962     6930325 :   for (i=zeros+1; i<=d; i++) alpha[i]=1;
    1963             : 
    1964     6156031 :   while (++kappa <= d)
    1965             :   {
    1966     4587534 :     if (kappa > kappamax)
    1967             :     {
    1968     3793408 :       if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
    1969     3793408 :       kappamax = kappa;
    1970     3793408 :       if (incgram)
    1971             :       {
    1972    15917695 :         for (i=zeros+1; i<=kappa; i++)
    1973    12324089 :           gmael(G,kappa,i) = ZV_dotproduct(gel(B,kappa), gel(B,i));
    1974     3593606 :         maxG = kappamax;
    1975             :       }
    1976             :     }
    1977             :     /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
    1978     4587279 :     if (Babai_dpe(av, kappa, &G,&B,&U, mu,r,s, alpha[kappa], zeros, maxG, &eta))
    1979           0 :     { *pG = incgram? NULL: G; *pB = B; *pU = U; return -1; }
    1980     9075437 :     if ((keepfirst && kappa == 2) ||
    1981     4487887 :         dpe_cmpmul(Dmael(r,kappa-1,kappa-1), &delta, Del(s,kappa-1)) <= 0)
    1982             :     { /* Step4: Success of Lovasz's condition */
    1983     4261808 :       alpha[kappa] = kappa;
    1984     4261808 :       dpe_submulz(Del(s,kappa-1), Dmael(mu,kappa,kappa-1), Dmael(r,kappa,kappa-1), Dmael(r,kappa,kappa));
    1985     4261812 :       continue;
    1986             :     }
    1987             :     /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
    1988      325742 :     if (DEBUGLEVEL>=4 && kappa==kappamax && Del(s,kappa-1)->d)
    1989           0 :       if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", Del(s,1)->e-1); }
    1990      325742 :     kappa2 = kappa;
    1991             :     do {
    1992      829288 :       kappa--;
    1993      829288 :       if (kappa < zeros+2 + (keepfirst ? 1: 0)) break;
    1994      707592 :     } while (dpe_cmpmul(Dmael(r,kappa-1,kappa-1), &delta, Del(s,kappa-1)) >= 0);
    1995      325742 :     update_alpha(alpha, kappa, kappa2, kappamax);
    1996             : 
    1997             :     /* Step6: Update the mu's and r's */
    1998      325742 :     dperotate(mu, kappa2, kappa);
    1999      325742 :     dperotate(r, kappa2, kappa);
    2000      325742 :     affdpe(Del(s,kappa), Dmael(r,kappa,kappa));
    2001             : 
    2002             :     /* Step7: Update G, B, U */
    2003      325742 :     if (U) rotate(U, kappa2, kappa);
    2004      325742 :     if (B) rotate(B, kappa2, kappa);
    2005      325742 :     rotateG(G,kappa2,kappa, maxG, Gtmp);
    2006             : 
    2007             :     /* Step8: Prepare the next loop iteration */
    2008      325742 :     if (kappa == zeros+1 && !signe(gmael(G,kappa,kappa)))
    2009             :     {
    2010       35161 :       zeros++; kappa++;
    2011       35161 :       affidpe(gmael(G,kappa,kappa), Dmael(r,kappa,kappa));
    2012             :     }
    2013             :   }
    2014     1568497 :   if (pr) *pr = dpeM_diagonal_shallow(r,d);
    2015     1568497 :   *pG = G; *pB = B; *pU = U; return zeros; /* success */
    2016             : }
    2017             : 
    2018             : 
    2019             : /************************** PROVED version (t_INT) *************************/
    2020             : 
    2021             : /* Babai's Nearest Plane algorithm (iterative).
    2022             :  * Size-reduces b_kappa using mu_{i,j} and r_{i,j} for j<=i <kappa
    2023             :  * Update B[,kappa]; compute mu_{kappa,j}, r_{kappa,j} for j<=kappa and s[kappa]
    2024             :  * mu, r, s updated in place (affrr). Return 1 on failure, else 0. */
    2025             : static int
    2026           0 : Babai(pari_sp av, long kappa, GEN *pG, GEN *pB, GEN *pU, GEN mu, GEN r, GEN s,
    2027             :       long a, long zeros, long maxG, GEN eta, long prec)
    2028             : {
    2029           0 :   GEN G = *pG, B = *pB, U = *pU, ztmp;
    2030           0 :   long k, aa = a > zeros? a: zeros+1;
    2031           0 :   const long n = B? nbrows(B): 0, d = U ? lg(U)-1: 0, bit = prec2nbits(prec);
    2032           0 :   long emaxmu = EX0, emax2mu = EX0;
    2033             :   /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
    2034             : 
    2035           0 :   for (;;) {
    2036           0 :     int go_on = 0;
    2037           0 :     long i, j, emax3mu = emax2mu;
    2038             : 
    2039           0 :     if (gc_needed(av,2))
    2040             :     {
    2041           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
    2042           0 :       gc_lll(av,3,&G,&B,&U);
    2043             :     }
    2044             :     /* Step2: compute the GSO for stage kappa */
    2045           0 :     emax2mu = emaxmu; emaxmu = EX0;
    2046           0 :     for (j=aa; j<kappa; j++)
    2047             :     {
    2048           0 :       pari_sp btop = avma;
    2049           0 :       GEN g = gmael(G,kappa,j);
    2050           0 :       for (k = zeros+1; k < j; k++)
    2051           0 :         g = mpsub(g, mulrr(gmael(mu,j,k), gmael(r,kappa,k)));
    2052           0 :       mpaff(g, gmael(r,kappa,j));
    2053           0 :       affrr(divrr(gmael(r,kappa,j), gmael(r,j,j)), gmael(mu,kappa,j));
    2054           0 :       emaxmu = maxss(emaxmu, expo(gmael(mu,kappa,j)));
    2055           0 :       set_avma(btop);
    2056             :     }
    2057           0 :     if (emax3mu != EX0 && emax3mu <= emax2mu + 5) /* precision too low */
    2058           0 :     { *pG = G; *pB = B; *pU = U; return 1; }
    2059             : 
    2060           0 :     for (j=kappa-1; j>zeros; j--)
    2061           0 :       if (abscmprr(gmael(mu,kappa,j), eta) > 0) { go_on = 1; break; }
    2062             : 
    2063             :     /* Step3--5: compute the X_j's  */
    2064           0 :     if (go_on)
    2065           0 :       for (j=kappa-1; j>zeros; j--)
    2066             :       {
    2067             :         pari_sp btop;
    2068           0 :         GEN tmp = gmael(mu,kappa,j);
    2069           0 :         if (absrsmall(tmp)) continue; /* size-reduced */
    2070             : 
    2071           0 :         if (gc_needed(av,2))
    2072             :         {
    2073           0 :           if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
    2074           0 :           gc_lll(av,3,&G,&B,&U);
    2075             :         }
    2076           0 :         btop = avma;
    2077             :         /* we consider separately the case |X| = 1 */
    2078           0 :         if (absrsmall2(tmp))
    2079             :         {
    2080           0 :           if (signe(tmp) > 0) { /* in this case, X = 1 */
    2081           0 :             for (k=zeros+1; k<j; k++)
    2082           0 :               affrr(subrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
    2083           0 :             set_avma(btop);
    2084           0 :             for (i=1; i<=n; i++)
    2085           0 :               gmael(B,kappa,i) = subii(gmael(B,kappa,i), gmael(B,j,i));
    2086           0 :             for (i=1; i<=d; i++)
    2087           0 :               gmael(U,kappa,i) = subii(gmael(U,kappa,i), gmael(U,j,i));
    2088           0 :             btop = avma;
    2089           0 :             ztmp = subii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
    2090           0 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    2091           0 :             gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
    2092           0 :             for (i=1; i<=j; i++)
    2093           0 :               gmael(G,kappa,i) = subii(gmael(G,kappa,i), gmael(G,j,i));
    2094           0 :             for (i=j+1; i<kappa; i++)
    2095           0 :               gmael(G,kappa,i) = subii(gmael(G,kappa,i), gmael(G,i,j));
    2096           0 :             for (i=kappa+1; i<=maxG; i++)
    2097           0 :               gmael(G,i,kappa) = subii(gmael(G,i,kappa), gmael(G,i,j));
    2098             :           } else { /* otherwise X = -1 */
    2099           0 :             for (k=zeros+1; k<j; k++)
    2100           0 :               affrr(addrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
    2101           0 :             set_avma(btop);
    2102           0 :             for (i=1; i<=n; i++)
    2103           0 :               gmael(B,kappa,i) = addii(gmael(B,kappa,i),gmael(B,j,i));
    2104           0 :             for (i=1; i<=d; i++)
    2105           0 :               gmael(U,kappa,i) = addii(gmael(U,kappa,i),gmael(U,j,i));
    2106           0 :             btop = avma;
    2107           0 :             ztmp = addii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
    2108           0 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    2109           0 :             gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
    2110           0 :             for (i=1; i<=j; i++)
    2111           0 :               gmael(G,kappa,i) = addii(gmael(G,kappa,i), gmael(G,j,i));
    2112           0 :             for (i=j+1; i<kappa; i++)
    2113           0 :               gmael(G,kappa,i) = addii(gmael(G,kappa,i), gmael(G,i,j));
    2114           0 :             for (i=kappa+1; i<=maxG; i++)
    2115           0 :               gmael(G,i,kappa) = addii(gmael(G,i,kappa), gmael(G,i,j));
    2116             :           }
    2117           0 :           continue;
    2118             :         }
    2119             :         /* we have |X| >= 2 */
    2120           0 :         if (expo(tmp) < BITS_IN_LONG)
    2121             :         {
    2122           0 :           ulong xx = roundr_safe(tmp)[2]; /* X fits in an ulong */
    2123           0 :           if (signe(tmp) > 0) /* = xx */
    2124             :           {
    2125           0 :             for (k=zeros+1; k<j; k++)
    2126           0 :               affrr(subrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
    2127           0 :                   gmael(mu,kappa,k));
    2128           0 :             set_avma(btop);
    2129           0 :             for (i=1; i<=n; i++)
    2130           0 :               gmael(B,kappa,i) = submuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
    2131           0 :             for (i=1; i<=d; i++)
    2132           0 :               gmael(U,kappa,i) = submuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
    2133           0 :             btop = avma;
    2134           0 :             ztmp = submuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
    2135           0 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    2136           0 :             gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
    2137           0 :             for (i=1; i<=j; i++)
    2138           0 :               gmael(G,kappa,i) = submuliu_inplace(gmael(G,kappa,i), gmael(G,j,i), xx);
    2139           0 :             for (i=j+1; i<kappa; i++)
    2140           0 :               gmael(G,kappa,i) = submuliu_inplace(gmael(G,kappa,i), gmael(G,i,j), xx);
    2141           0 :             for (i=kappa+1; i<=maxG; i++)
    2142           0 :               gmael(G,i,kappa) = submuliu_inplace(gmael(G,i,kappa), gmael(G,i,j), xx);
    2143             :           }
    2144             :           else /* = -xx */
    2145             :           {
    2146           0 :             for (k=zeros+1; k<j; k++)
    2147           0 :               affrr(addrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
    2148           0 :                   gmael(mu,kappa,k));
    2149           0 :             set_avma(btop);
    2150           0 :             for (i=1; i<=n; i++)
    2151           0 :               gmael(B,kappa,i) = addmuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
    2152           0 :             for (i=1; i<=d; i++)
    2153           0 :               gmael(U,kappa,i) = addmuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
    2154           0 :             btop = avma;
    2155           0 :             ztmp = addmuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
    2156           0 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    2157           0 :             gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
    2158           0 :             for (i=1; i<=j; i++)
    2159           0 :               gmael(G,kappa,i) = addmuliu_inplace(gmael(G,kappa,i), gmael(G,j,i), xx);
    2160           0 :             for (i=j+1; i<kappa; i++)
    2161           0 :               gmael(G,kappa,i) = addmuliu_inplace(gmael(G,kappa,i), gmael(G,i,j), xx);
    2162           0 :             for (i=kappa+1; i<=maxG; i++)
    2163           0 :               gmael(G,i,kappa) = addmuliu_inplace(gmael(G,i,kappa), gmael(G,i,j), xx);
    2164             :           }
    2165             :         }
    2166             :         else
    2167             :         {
    2168             :           long e;
    2169           0 :           GEN X = truncexpo(tmp, bit, &e); /* tmp ~ X * 2^e */
    2170           0 :           btop = avma;
    2171           0 :           for (k=zeros+1; k<j; k++)
    2172             :           {
    2173           0 :             GEN x = mulir(X, gmael(mu,j,k));
    2174           0 :             if (e) shiftr_inplace(x, e);
    2175           0 :             affrr(subrr(gmael(mu,kappa,k), x), gmael(mu,kappa,k));
    2176             :           }
    2177           0 :           set_avma(btop);
    2178           0 :           for (i=1; i<=n; i++)
    2179           0 :             gmael(B,kappa,i) = submulshift(gmael(B,kappa,i), gmael(B,j,i), X, e);
    2180           0 :           for (i=1; i<=d; i++)
    2181           0 :             gmael(U,kappa,i) = submulshift(gmael(U,kappa,i), gmael(U,j,i), X, e);
    2182           0 :           btop = avma;
    2183           0 :           ztmp = submulshift(mulshift(gmael(G,j,j), sqri(X), 2*e),
    2184           0 :               gmael(G,kappa,j), X, e+1);
    2185           0 :           ztmp = addii(gmael(G,kappa,kappa), ztmp);
    2186           0 :           gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
    2187           0 :           for (i=1; i<=j; i++)
    2188           0 :             gmael(G,kappa,i) = submulshift(gmael(G,kappa,i), gmael(G,j,i), X, e);
    2189           0 :           for (   ; i<kappa; i++)
    2190           0 :             gmael(G,kappa,i) = submulshift(gmael(G,kappa,i), gmael(G,i,j), X, e);
    2191           0 :           for (i=kappa+1; i<=maxG; i++)
    2192           0 :             gmael(G,i,kappa) = submulshift(gmael(G,i,kappa), gmael(G,i,j), X, e);
    2193             :         }
    2194             :       }
    2195           0 :     if (!go_on) break; /* Anything happened? */
    2196           0 :     aa = zeros+1;
    2197             :   }
    2198             : 
    2199           0 :   affir(gmael(G,kappa,kappa), gel(s,zeros+1));
    2200             :   /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
    2201           0 :   av = avma;
    2202           0 :   for (k=zeros+1; k<=kappa-2; k++)
    2203           0 :     affrr(subrr(gel(s,k), mulrr(gmael(mu,kappa,k), gmael(r,kappa,k))),
    2204           0 :           gel(s,k+1));
    2205           0 :   *pG = G; *pB = B; *pU = U; return gc_bool(av, 0);
    2206             : }
    2207             : 
    2208             : /* G integral Gram matrix, LLL-reduces (G,B,U) in place [apply base change
    2209             :  * transforms to B and U]. If (keepfirst), never swap with first vector.
    2210             :  * If G = NULL, we compute the Gram matrix incrementally.
    2211             :  * Return -1 on failure, else zeros = dim Kernel (>= 0) */
    2212             : static long
    2213           0 : fplll(GEN *pG, GEN *pB, GEN *pU, GEN *pr, double DELTA, double ETA,
    2214             :       long keepfirst, long prec)
    2215             : {
    2216             :   pari_sp av, av2;
    2217           0 :   GEN mu, r, s, tmp, Gtmp, alpha, G = *pG, B = *pB, U = *pU;
    2218           0 :   GEN delta = dbltor(DELTA), eta = dbltor(ETA);
    2219           0 :   long d, maxG, kappa, kappa2, i, j, zeros, kappamax, incgram = !G, cnt = 0;
    2220             : 
    2221           0 :   if (incgram)
    2222             :   { /* incremental Gram matrix */
    2223           0 :     maxG = 2; d = lg(B)-1;
    2224           0 :     G = zeromatcopy(d, d);
    2225             :   }
    2226             :   else
    2227           0 :     maxG = d = lg(G)-1;
    2228             : 
    2229           0 :   mu = cgetg(d+1, t_MAT);
    2230           0 :   r  = cgetg(d+1, t_MAT);
    2231           0 :   s  = cgetg(d+1, t_VEC);
    2232           0 :   for (j = 1; j <= d; j++)
    2233             :   {
    2234           0 :     GEN M = cgetg(d+1, t_COL), R = cgetg(d+1, t_COL);
    2235           0 :     gel(mu,j)= M;
    2236           0 :     gel(r,j) = R;
    2237           0 :     gel(s,j) = cgetr(prec);
    2238           0 :     for (i = 1; i <= d; i++)
    2239             :     {
    2240           0 :       gel(R,i) = cgetr(prec);
    2241           0 :       gel(M,i) = cgetr(prec);
    2242             :     }
    2243             :   }
    2244           0 :   Gtmp = cgetg(d+1, t_VEC);
    2245           0 :   alpha = cgetg(d+1, t_VECSMALL);
    2246           0 :   av = avma;
    2247             : 
    2248             :   /* Step2: Initializing the main loop */
    2249           0 :   kappamax = 1;
    2250           0 :   i = 1;
    2251             :   do {
    2252           0 :     if (incgram) gmael(G,i,i) = ZV_dotsquare(gel(B,i));
    2253           0 :     affir(gmael(G,i,i), gmael(r,i,i));
    2254           0 :   } while (!signe(gmael(G,i,i)) && ++i <= d);
    2255           0 :   zeros = i-1; /* all basis vectors b_i with i <= zeros are zero vectors */
    2256           0 :   kappa = i;
    2257           0 :   for (i=zeros+1; i<=d; i++) alpha[i]=1;
    2258             : 
    2259           0 :   while (++kappa <= d)
    2260             :   {
    2261           0 :     if (kappa > kappamax)
    2262             :     {
    2263           0 :       if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
    2264           0 :       kappamax = kappa;
    2265           0 :       if (incgram)
    2266             :       {
    2267           0 :         for (i=zeros+1; i<=kappa; i++)
    2268           0 :           gmael(G,kappa,i) = ZV_dotproduct(gel(B,kappa), gel(B,i));
    2269           0 :         maxG = kappamax;
    2270             :       }
    2271             :     }
    2272             :     /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
    2273           0 :     if (Babai(av, kappa, &G,&B,&U, mu,r,s, alpha[kappa], zeros, maxG, eta, prec))
    2274           0 :     { *pG = incgram? NULL: G; *pB = B; *pU = U; return -1; }
    2275           0 :     av2 = avma;
    2276           0 :     if ((keepfirst && kappa == 2) ||
    2277           0 :         cmprr(mulrr(gmael(r,kappa-1,kappa-1), delta), gel(s,kappa-1)) <= 0)
    2278             :     { /* Step4: Success of Lovasz's condition */
    2279           0 :       alpha[kappa] = kappa;
    2280           0 :       tmp = mulrr(gmael(mu,kappa,kappa-1), gmael(r,kappa,kappa-1));
    2281           0 :       affrr(subrr(gel(s,kappa-1), tmp), gmael(r,kappa,kappa));
    2282           0 :       set_avma(av2); continue;
    2283             :     }
    2284             :     /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
    2285           0 :     if (DEBUGLEVEL>=4 && kappa==kappamax && signe(gel(s,kappa-1)))
    2286           0 :       if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", expo(gel(s,1))); }
    2287           0 :     kappa2 = kappa;
    2288             :     do {
    2289           0 :       kappa--;
    2290           0 :       if (kappa < zeros+2 + (keepfirst ? 1: 0)) break;
    2291           0 :       tmp = mulrr(gmael(r,kappa-1,kappa-1), delta);
    2292           0 :     } while (cmprr(gel(s,kappa-1), tmp) <= 0);
    2293           0 :     set_avma(av2);
    2294           0 :     update_alpha(alpha, kappa, kappa2, kappamax);
    2295             : 
    2296             :     /* Step6: Update the mu's and r's */
    2297           0 :     rotate(mu, kappa2, kappa);
    2298           0 :     rotate(r, kappa2, kappa);
    2299           0 :     affrr(gel(s,kappa), gmael(r,kappa,kappa));
    2300             : 
    2301             :     /* Step7: Update G, B, U */
    2302           0 :     if (U) rotate(U, kappa2, kappa);
    2303           0 :     if (B) rotate(B, kappa2, kappa);
    2304           0 :     rotateG(G,kappa2,kappa, maxG, Gtmp);
    2305             : 
    2306             :     /* Step8: Prepare the next loop iteration */
    2307           0 :     if (kappa == zeros+1 && !signe(gmael(G,kappa,kappa)))
    2308             :     {
    2309           0 :       zeros++; kappa++;
    2310           0 :       affir(gmael(G,kappa,kappa), gmael(r,kappa,kappa));
    2311             :     }
    2312             :   }
    2313           0 :   if (pr) *pr = RgM_diagonal_shallow(r);
    2314           0 :   *pG = G; *pB = B; *pU = U; return zeros; /* success */
    2315             : }
    2316             : 
    2317             : /* do not support LLL_KER, LLL_ALL, LLL_KEEP_FIRST */
    2318             : static GEN
    2319     4702971 : ZM2_lll_norms(GEN x, long flag, GEN *pN)
    2320             : {
    2321             :   GEN a,b,c,d;
    2322             :   GEN G, U;
    2323     4702971 :   if (flag & LLL_GRAM)
    2324        7355 :     G = x;
    2325             :   else
    2326     4695616 :     G = gram_matrix(x);
    2327     4702957 :   a = gcoeff(G,1,1); b = shifti(gcoeff(G,1,2),1); c = gcoeff(G,2,2);
    2328     4702917 :   d = qfb_disc3(a,b,c);
    2329     4702910 :   if (signe(d)>=0) return NULL;
    2330     4702525 :   G = redimagsl2(mkqfb(a,b,c,d),&U);
    2331     4702588 :   if (pN) (void) RgM_gram_schmidt(G, pN);
    2332     4702588 :   if (flag & LLL_INPLACE) return ZM2_mul(x,U);
    2333     4702588 :   return U;
    2334             : }
    2335             : 
    2336             : static void
    2337      617541 : fplll_flatter(GEN *pG, GEN *pB, GEN *pU, long rank, long flag)
    2338             : {
    2339      617541 :   if (!*pG)
    2340             :   {
    2341      616580 :     GEN T = ZM_flatter_rank(*pB, rank, flag);
    2342      616582 :     if (T)
    2343             :     {
    2344      327200 :       if (*pU)
    2345             :       {
    2346      313199 :         *pU = ZM_mul(*pU, T);
    2347      313199 :         *pB = ZM_mul(*pB, T);
    2348             :       }
    2349       14001 :       else *pB = T;
    2350             :     }
    2351             :   }
    2352             :   else
    2353             :   {
    2354         961 :     GEN T, G = *pG;
    2355         961 :     long i, j, l = lg(G);
    2356        7207 :     for (i = 1; i < l; i++)
    2357       43383 :       for(j = 1; j < i; j++) gmael(G,j,i) = gmael(G,i,j);
    2358         961 :     T = ZM_flattergram_rank(G, rank, flag);
    2359         961 :     if (T)
    2360             :     {
    2361         961 :       if (*pU) *pU = ZM_mul(*pU, T);
    2362         961 :       *pG = qf_ZM_apply(*pG, T);
    2363             :     }
    2364             :   }
    2365      617543 : }
    2366             : 
    2367             : static GEN
    2368     1083008 : get_gramschmidt(GEN M, long rank)
    2369             : {
    2370             :   GEN B, Q, L;
    2371     1083008 :   long r = lg(M)-1, bitprec = 3*r + 30;
    2372     1083008 :   long prec = nbits2prec64(bitprec);
    2373     1083009 :   if (rank < r) M = vconcat(gshift(M,1), matid(r));
    2374     1083008 :   if (!QR_init(RgM_gtofp(M, prec), &B, &Q, &L, prec) || !gsisinv(L)) return NULL;
    2375      467498 :   return L;
    2376             : }
    2377             : 
    2378             : static GEN
    2379       44316 : get_gaussred(GEN M, long rank)
    2380             : {
    2381       44316 :   pari_sp ltop = avma;
    2382       44316 :   long r = lg(M)-1, bitprec = 3*r + 30, prec = nbits2prec64(bitprec);
    2383             :   GEN R;
    2384       44316 :   if (rank < r) M = RgM_Rg_add(gshift(M, 1), gen_1);
    2385       44316 :   R = RgM_Cholesky(RgM_gtofp(M, prec), prec);
    2386       44316 :   if (!R) return NULL;
    2387       43355 :   return gerepilecopy(ltop, R);
    2388             : }
    2389             : 
    2390             : /* Assume x a ZM, if pN != NULL, set it to Gram-Schmidt (squared) norms
    2391             :  * The following modes are supported:
    2392             :  * - flag & LLL_INPLACE: x a lattice basis, return x*U
    2393             :  * - flag & LLL_GRAM: x a Gram matrix / else x a lattice basis; return
    2394             :  *     LLL base change matrix U [LLL_IM]
    2395             :  *     kernel basis [LLL_KER, nonreduced]
    2396             :  *     both [LLL_ALL] */
    2397             : GEN
    2398     6957379 : ZM_lll_norms(GEN x, double DELTA, long flag, GEN *pN)
    2399             : {
    2400     6957379 :   pari_sp av = avma;
    2401     6957379 :   const double ETA = 0.51;
    2402     6957379 :   const long keepfirst = flag & LLL_KEEP_FIRST;
    2403     6957379 :   long p, zeros = -1, n = lg(x)-1, is_upper, is_lower, useflatter = 0, rank;
    2404     6957379 :   GEN G, B, U, L = NULL;
    2405             :   pari_timer T;
    2406     6957379 :   long thre[]={31783,34393,20894,22525,13533,1928,672,671,
    2407             :                 422,506,315,313,222,205,167,154,139,138,
    2408             :                 110,120,98,94,81,75,74,64,74,74,
    2409             :                 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,
    2410             :                 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};
    2411     6957379 :   long thsn[]={23280,30486,50077,44136,78724,15690,1801,1611,
    2412             :                981,1359,978,1042,815,866,788,775,726,712,
    2413             :                626,613,548,564,474,481,504,447,453,508,
    2414             :                705,794,1008,946,767,898,886,763,842,757,
    2415             :                725,774,639,655,705,627,635,704,511,613,
    2416             :                583,595,568,640,541,640,567,540,577,584,
    2417             :                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};
    2418     6957379 :   if (n <= 1) return lll_trivial(x, flag);
    2419     6847729 :   if (nbrows(x)==0)
    2420             :   {
    2421       15184 :     if (flag & LLL_KER) return matid(n);
    2422       15184 :     if (flag & (LLL_INPLACE|LLL_IM)) return cgetg(1,t_MAT);
    2423           0 :     retmkvec2(matid(n), cgetg(1,t_MAT));
    2424             :   }
    2425     6832744 :   if (n==2 && nbrows(x)==2  && (flag&LLL_IM) && !keepfirst)
    2426             :   {
    2427     4702972 :     U = ZM2_lll_norms(x, flag, pN);
    2428     4702973 :     if (U) return U;
    2429             :   }
    2430     2130154 :   if (flag & LLL_GRAM)
    2431       60430 :   { G = x; B = NULL; U = matid(n); is_upper = 0; is_lower = 0; }
    2432             :   else
    2433             :   {
    2434     2069724 :     G = NULL; B = x; U = (flag & LLL_INPLACE)? NULL: matid(n);
    2435     2069728 :     is_upper = (flag & LLL_UPPER) || ZM_is_upper(B);
    2436     2069727 :     is_lower = !B || is_upper || keepfirst ? 0: ZM_is_lower(B);
    2437     2069726 :     if (is_lower) L = RgM_flip(B);
    2438             :   }
    2439     2130156 :   rank = (flag&LLL_NOFLATTER) ? 0: ZM_rank(x);
    2440     2130140 :   if (n > 2 && !(flag&LLL_NOFLATTER))
    2441     1733889 :   {
    2442     1689572 :     GEN R = B ? (is_upper ? B : (is_lower ? L : get_gramschmidt(B, rank)))
    2443     3423449 :               : get_gaussred(G, rank);
    2444     1733881 :     if (R)
    2445             :     {
    2446     1117418 :       long spr = spread(R), sz = mpexpo(gsupnorm(R, DEFAULTPREC)), thr;
    2447     1117426 :       if (DEBUGLEVEL>=5) err_printf("LLL: dim %ld, size %ld, spread %ld\n",n, sz, spr);
    2448     1117426 :       if ((is_upper && ZM_is_knapsack(B)) || (is_lower && ZM_is_knapsack(L)))
    2449       92364 :         thr = thsn[minss(n-3,numberof(thsn)-1)];
    2450             :       else
    2451             :       {
    2452     1025062 :         thr = thre[minss(n-3,numberof(thre)-1)];
    2453     1025062 :         if (n>=10) sz = spr;
    2454             :       }
    2455     1117426 :       useflatter = sz >= thr;
    2456             :     } else
    2457      616463 :       useflatter = 1;
    2458      396253 :   } else useflatter = 0;
    2459     2130142 :   if(DEBUGLEVEL>=4) timer_start(&T);
    2460     2130142 :   if (useflatter)
    2461             :   {
    2462      617541 :     if (is_lower)
    2463             :     {
    2464           0 :       fplll_flatter(&G, &L, &U, rank, flag | LLL_UPPER);
    2465           0 :       B = RgM_flop(L);
    2466           0 :       if (U) U = RgM_flop(U);
    2467             :     }
    2468             :     else
    2469      617541 :       fplll_flatter(&G, &B, &U, rank, flag | (is_upper? LLL_UPPER:0));
    2470      617543 :     if (DEBUGLEVEL>=4  && !(flag & LLL_NOCERTIFY))
    2471           0 :       timer_printf(&T, "FLATTER");
    2472             :   }
    2473     2130144 :   if (!(flag & LLL_GRAM))
    2474             :   {
    2475             :     long t;
    2476     2069714 :     B = gcopy(B);
    2477     2069716 :     if(DEBUGLEVEL>=4)
    2478           0 :       err_printf("Entering L^2 (double): dim %ld, LLL-parameters (%.3f,%.3f)\n",
    2479             :                  n, DELTA,ETA);
    2480     2069716 :     zeros = fplll_fast(&B, &U, DELTA, ETA, keepfirst);
    2481     2069724 :     if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
    2482     2073857 :     for (p = DEFAULTPREC, t = 0; zeros < 0 && t < 1 ; p += EXTRAPREC64, t++)
    2483             :     {
    2484        4133 :       if (DEBUGLEVEL>=4)
    2485           0 :         err_printf("Entering L^2 (heuristic): LLL-parameters (%.3f,%.3f), prec = %d/%d\n", DELTA, ETA, p, p);
    2486        4133 :       zeros = fplll_heuristic(&B, &U, DELTA, ETA, keepfirst, p, p);
    2487        4133 :       gc_lll(av, 2, &B, &U);
    2488        4133 :       if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
    2489             :     }
    2490             :   } else
    2491       60430 :     G = gcopy(G);
    2492     2130154 :   if (zeros < 0 || !(flag & LLL_NOCERTIFY))
    2493             :   {
    2494     1568515 :     if(DEBUGLEVEL>=4)
    2495           0 :       err_printf("Entering L^2 (dpe): LLL-parameters (%.3f,%.3f)\n", DELTA,ETA);
    2496     1568515 :     zeros = fplll_dpe(&G, &B, &U, pN, DELTA, ETA, keepfirst);
    2497     1568496 :     if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
    2498     1568497 :     if (zeros < 0)
    2499           0 :       for (p = DEFAULTPREC;; p += EXTRAPREC64)
    2500             :       {
    2501           0 :         if (DEBUGLEVEL>=4)
    2502           0 :           err_printf("Entering L^2: LLL-parameters (%.3f,%.3f), prec = %d\n",
    2503             :               DELTA,ETA, p);
    2504           0 :         zeros = fplll(&G, &B, &U, pN, DELTA, ETA, keepfirst, p);
    2505           0 :         if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
    2506           0 :         if (zeros >= 0) break;
    2507           0 :         gc_lll(av, 3, &G, &B, &U);
    2508             :       }
    2509             :   }
    2510     2130136 :   return lll_finish(U? U: B, zeros, flag);
    2511             : }
    2512             : 
    2513             : /********************************************************************/
    2514             : /**                                                                **/
    2515             : /**                        LLL OVER K[X]                           **/
    2516             : /**                                                                **/
    2517             : /********************************************************************/
    2518             : static int
    2519         504 : pslg(GEN x)
    2520             : {
    2521             :   long tx;
    2522         504 :   if (gequal0(x)) return 2;
    2523         448 :   tx = typ(x); return is_scalar_t(tx)? 3: lg(x);
    2524             : }
    2525             : 
    2526             : static int
    2527         196 : REDgen(long k, long l, GEN h, GEN L, GEN B)
    2528             : {
    2529         196 :   GEN q, u = gcoeff(L,k,l);
    2530             :   long i;
    2531             : 
    2532         196 :   if (pslg(u) < pslg(B)) return 0;
    2533             : 
    2534         140 :   q = gneg(gdeuc(u,B));
    2535         140 :   gel(h,k) = gadd(gel(h,k), gmul(q,gel(h,l)));
    2536         140 :   for (i=1; i<l; i++) gcoeff(L,k,i) = gadd(gcoeff(L,k,i), gmul(q,gcoeff(L,l,i)));
    2537         140 :   gcoeff(L,k,l) = gadd(gcoeff(L,k,l), gmul(q,B)); return 1;
    2538             : }
    2539             : 
    2540             : static int
    2541         196 : do_SWAPgen(GEN h, GEN L, GEN B, long k, GEN fl, int *flc)
    2542             : {
    2543             :   GEN p1, la, la2, Bk;
    2544             :   long ps1, ps2, i, j, lx;
    2545             : 
    2546         196 :   if (!fl[k-1]) return 0;
    2547             : 
    2548         140 :   la = gcoeff(L,k,k-1); la2 = gsqr(la);
    2549         140 :   Bk = gel(B,k);
    2550         140 :   if (fl[k])
    2551             :   {
    2552          56 :     GEN q = gadd(la2, gmul(gel(B,k-1),gel(B,k+1)));
    2553          56 :     ps1 = pslg(gsqr(Bk));
    2554          56 :     ps2 = pslg(q);
    2555          56 :     if (ps1 <= ps2 && (ps1 < ps2 || !*flc)) return 0;
    2556          28 :     *flc = (ps1 != ps2);
    2557          28 :     gel(B,k) = gdiv(q, Bk);
    2558             :   }
    2559             : 
    2560         112 :   swap(gel(h,k-1), gel(h,k)); lx = lg(L);
    2561         112 :   for (j=1; j<k-1; j++) swap(gcoeff(L,k-1,j), gcoeff(L,k,j));
    2562         112 :   if (fl[k])
    2563             :   {
    2564          28 :     for (i=k+1; i<lx; i++)
    2565             :     {
    2566           0 :       GEN t = gcoeff(L,i,k);
    2567           0 :       p1 = gsub(gmul(gel(B,k+1),gcoeff(L,i,k-1)), gmul(la,t));
    2568           0 :       gcoeff(L,i,k) = gdiv(p1, Bk);
    2569           0 :       p1 = gadd(gmul(la,gcoeff(L,i,k-1)), gmul(gel(B,k-1),t));
    2570           0 :       gcoeff(L,i,k-1) = gdiv(p1, Bk);
    2571             :     }
    2572             :   }
    2573          84 :   else if (!gequal0(la))
    2574             :   {
    2575          28 :     p1 = gdiv(la2, Bk);
    2576          28 :     gel(B,k+1) = gel(B,k) = p1;
    2577          28 :     for (i=k+2; i<=lx; i++) gel(B,i) = gdiv(gmul(p1,gel(B,i)),Bk);
    2578          28 :     for (i=k+1; i<lx; i++)
    2579           0 :       gcoeff(L,i,k-1) = gdiv(gmul(la,gcoeff(L,i,k-1)), Bk);
    2580          28 :     for (j=k+1; j<lx-1; j++)
    2581           0 :       for (i=j+1; i<lx; i++)
    2582           0 :         gcoeff(L,i,j) = gdiv(gmul(p1,gcoeff(L,i,j)), Bk);
    2583             :   }
    2584             :   else
    2585             :   {
    2586          56 :     gcoeff(L,k,k-1) = gen_0;
    2587          56 :     for (i=k+1; i<lx; i++)
    2588             :     {
    2589           0 :       gcoeff(L,i,k) = gcoeff(L,i,k-1);
    2590           0 :       gcoeff(L,i,k-1) = gen_0;
    2591             :     }
    2592          56 :     gel(B,k) = gel(B,k-1); fl[k] = 1; fl[k-1] = 0;
    2593             :   }
    2594         112 :   return 1;
    2595             : }
    2596             : 
    2597             : static void
    2598         168 : incrementalGSgen(GEN x, GEN L, GEN B, long k, GEN fl)
    2599             : {
    2600         168 :   GEN u = NULL; /* gcc -Wall */
    2601             :   long i, j;
    2602         420 :   for (j = 1; j <= k; j++)
    2603         252 :     if (j==k || fl[j])
    2604             :     {
    2605         252 :       u = gcoeff(x,k,j);
    2606         252 :       if (!is_extscalar_t(typ(u))) pari_err_TYPE("incrementalGSgen",u);
    2607         336 :       for (i=1; i<j; i++)
    2608          84 :         if (fl[i])
    2609             :         {
    2610          84 :           u = gsub(gmul(gel(B,i+1),u), gmul(gcoeff(L,k,i),gcoeff(L,j,i)));
    2611          84 :           u = gdiv(u, gel(B,i));
    2612             :         }
    2613         252 :       gcoeff(L,k,j) = u;
    2614             :     }
    2615         168 :   if (gequal0(u)) gel(B,k+1) = gel(B,k);
    2616             :   else
    2617             :   {
    2618         112 :     gel(B,k+1) = gcoeff(L,k,k); gcoeff(L,k,k) = gen_1; fl[k] = 1;
    2619             :   }
    2620         168 : }
    2621             : 
    2622             : static GEN
    2623         168 : lllgramallgen(GEN x, long flag)
    2624             : {
    2625         168 :   long lx = lg(x), i, j, k, l, n;
    2626             :   pari_sp av;
    2627             :   GEN B, L, h, fl;
    2628             :   int flc;
    2629             : 
    2630         168 :   n = lx-1; if (n<=1) return lll_trivial(x,flag);
    2631          84 :   if (lgcols(x) != lx) pari_err_DIM("lllgramallgen");
    2632             : 
    2633          84 :   fl = cgetg(lx, t_VECSMALL);
    2634             : 
    2635          84 :   av = avma;
    2636          84 :   B = scalarcol_shallow(gen_1, lx);
    2637          84 :   L = cgetg(lx,t_MAT);
    2638         252 :   for (j=1; j<lx; j++) { gel(L,j) = zerocol(n); fl[j] = 0; }
    2639             : 
    2640          84 :   h = matid(n);
    2641         252 :   for (i=1; i<lx; i++)
    2642         168 :     incrementalGSgen(x, L, B, i, fl);
    2643          84 :   flc = 0;
    2644          84 :   for(k=2;;)
    2645             :   {
    2646         196 :     if (REDgen(k, k-1, h, L, gel(B,k))) flc = 1;
    2647         196 :     if (do_SWAPgen(h, L, B, k, fl, &flc)) { if (k > 2) k--; }
    2648             :     else
    2649             :     {
    2650          84 :       for (l=k-2; l>=1; l--)
    2651           0 :         if (REDgen(k, l, h, L, gel(B,l+1))) flc = 1;
    2652          84 :       if (++k > n) break;
    2653             :     }
    2654         112 :     if (gc_needed(av,1))
    2655             :     {
    2656           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"lllgramallgen");
    2657           0 :       gerepileall(av,3,&B,&L,&h);
    2658             :     }
    2659             :   }
    2660         140 :   k=1; while (k<lx && !fl[k]) k++;
    2661          84 :   return lll_finish(h,k-1,flag);
    2662             : }
    2663             : 
    2664             : static GEN
    2665         168 : lllallgen(GEN x, long flag)
    2666             : {
    2667         168 :   pari_sp av = avma;
    2668         168 :   if (!(flag & LLL_GRAM)) x = gram_matrix(x);
    2669          84 :   else if (!RgM_is_square_mat(x)) pari_err_DIM("qflllgram");
    2670         168 :   return gerepilecopy(av, lllgramallgen(x, flag));
    2671             : }
    2672             : GEN
    2673          42 : lllgen(GEN x) { return lllallgen(x, LLL_IM); }
    2674             : GEN
    2675          42 : lllkerimgen(GEN x) { return lllallgen(x, LLL_ALL); }
    2676             : GEN
    2677          42 : lllgramgen(GEN x)  { return lllallgen(x, LLL_IM|LLL_GRAM); }
    2678             : GEN
    2679          42 : lllgramkerimgen(GEN x)  { return lllallgen(x, LLL_ALL|LLL_GRAM); }
    2680             : 
    2681             : static GEN
    2682       36699 : lllall(GEN x, long flag)
    2683       36699 : { pari_sp av = avma; return gerepilecopy(av, ZM_lll(x, LLLDFT, flag)); }
    2684             : GEN
    2685         183 : lllint(GEN x) { return lllall(x, LLL_IM); }
    2686             : GEN
    2687          35 : lllkerim(GEN x) { return lllall(x, LLL_ALL); }
    2688             : GEN
    2689       36439 : lllgramint(GEN x)
    2690       36439 : { if (!RgM_is_square_mat(x)) pari_err_DIM("qflllgram");
    2691       36439 :   return lllall(x, LLL_IM | LLL_GRAM); }
    2692             : GEN
    2693          35 : lllgramkerim(GEN x)
    2694          35 : { if (!RgM_is_square_mat(x)) pari_err_DIM("qflllgram");
    2695          35 :   return lllall(x, LLL_ALL | LLL_GRAM); }
    2696             : 
    2697             : GEN
    2698     5271033 : lllfp(GEN x, double D, long flag)
    2699             : {
    2700     5271033 :   long n = lg(x)-1;
    2701     5271033 :   pari_sp av = avma;
    2702             :   GEN h;
    2703     5271033 :   if (n <= 1) return lll_trivial(x,flag);
    2704     4610582 :   if (flag & LLL_GRAM)
    2705             :   {
    2706        9270 :     if (!RgM_is_square_mat(x)) pari_err_DIM("qflllgram");
    2707        9256 :     if (isinexact(x))
    2708             :     {
    2709        9165 :       x = RgM_Cholesky(x, gprecision(x));
    2710        9165 :       if (!x) return NULL;
    2711        9165 :       flag &= ~LLL_GRAM;
    2712             :     }
    2713             :   }
    2714     4610568 :   h = ZM_lll(RgM_rescale_to_int(x), D, flag);
    2715     4610520 :   return gerepilecopy(av, h);
    2716             : }
    2717             : 
    2718             : GEN
    2719        9089 : lllgram(GEN x) { return lllfp(x,LLLDFT,LLL_GRAM|LLL_IM); }
    2720             : GEN
    2721     1175039 : lll(GEN x) { return lllfp(x,LLLDFT,LLL_IM); }
    2722             : 
    2723             : static GEN
    2724          63 : qflllgram(GEN x)
    2725             : {
    2726          63 :   GEN T = lllgram(x);
    2727          42 :   if (!T) pari_err_PREC("qflllgram");
    2728          42 :   return T;
    2729             : }
    2730             : 
    2731             : GEN
    2732         301 : qflll0(GEN x, long flag)
    2733             : {
    2734         301 :   if (typ(x) != t_MAT) pari_err_TYPE("qflll",x);
    2735         301 :   switch(flag)
    2736             :   {
    2737          49 :     case 0: return lll(x);
    2738          63 :     case 1: return lllfp(x, LLLDFT, LLL_IM | LLL_NOFLATTER);
    2739          49 :     case 2: RgM_check_ZM(x,"qflll"); return lllintpartial(x);
    2740           7 :     case 3: RgM_check_ZM(x,"qflll"); return lllall(x, LLL_INPLACE);
    2741          49 :     case 4: RgM_check_ZM(x,"qflll"); return lllkerim(x);
    2742          42 :     case 5: return lllkerimgen(x);
    2743          42 :     case 8: return lllgen(x);
    2744           0 :     default: pari_err_FLAG("qflll");
    2745             :   }
    2746             :   return NULL; /* LCOV_EXCL_LINE */
    2747             : }
    2748             : 
    2749             : GEN
    2750         245 : qflllgram0(GEN x, long flag)
    2751             : {
    2752         245 :   if (typ(x) != t_MAT) pari_err_TYPE("qflllgram",x);
    2753         245 :   switch(flag)
    2754             :   {
    2755          63 :     case 0: return qflllgram(x);
    2756          49 :     case 1: return lllfp(x, LLLDFT, LLL_IM | LLL_GRAM | LLL_NOFLATTER);
    2757          49 :     case 4: RgM_check_ZM(x,"qflllgram"); return lllgramkerim(x);
    2758          42 :     case 5: return lllgramkerimgen(x);
    2759          42 :     case 8: return lllgramgen(x);
    2760           0 :     default: pari_err_FLAG("qflllgram");
    2761             :   }
    2762             :   return NULL; /* LCOV_EXCL_LINE */
    2763             : }
    2764             : 
    2765             : /********************************************************************/
    2766             : /**                                                                **/
    2767             : /**                   INTEGRAL KERNEL (LLL REDUCED)                **/
    2768             : /**                                                                **/
    2769             : /********************************************************************/
    2770             : static GEN
    2771          70 : kerint0(GEN M)
    2772             : {
    2773             :   /* return ZM_lll(M, LLLDFT, LLL_KER); */
    2774          70 :   GEN U, H = ZM_hnflll(M,&U,1);
    2775          70 :   long d = lg(M)-lg(H);
    2776          70 :   if (!d) return cgetg(1, t_MAT);
    2777          70 :   return ZM_lll(vecslice(U,1,d), LLLDFT, LLL_INPLACE);
    2778             : }
    2779             : GEN
    2780          42 : kerint(GEN M)
    2781             : {
    2782          42 :   pari_sp av = avma;
    2783          42 :   return gerepilecopy(av, kerint0(M));
    2784             : }
    2785             : /* OBSOLETE: use kerint */
    2786             : GEN
    2787          28 : matkerint0(GEN M, long flag)
    2788             : {
    2789          28 :   pari_sp av = avma;
    2790          28 :   if (typ(M) != t_MAT) pari_err_TYPE("matkerint",M);
    2791          28 :   M = Q_primpart(M);
    2792          28 :   RgM_check_ZM(M, "kerint");
    2793          28 :   switch(flag)
    2794             :   {
    2795          28 :     case 0:
    2796          28 :     case 1: return gerepilecopy(av, kerint0(M));
    2797           0 :     default: pari_err_FLAG("matkerint");
    2798             :   }
    2799             :   return NULL; /* LCOV_EXCL_LINE */
    2800             : }

Generated by: LCOV version 1.16