Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is to exceed 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - basemath - lll.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.16.1 lcov report (development 28695-49bb1ac00f) Lines: 1046 1302 80.3 %
Date: 2023-09-24 07:47:42 Functions: 98 100 98.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2008  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : #include "pari.h"
      16             : #include "paripriv.h"
      17             : 
      18             : #define DEBUGLEVEL DEBUGLEVEL_qflll
      19             : 
      20             : static double
      21    17184670 : pari_rint(double a)
      22             : {
      23             : #ifdef HAS_RINT
      24    17184670 :   return rint(a);
      25             : #else
      26             :   const double two_to_52 = 4.5035996273704960e+15;
      27             :   double fa = fabs(a);
      28             :   double r = two_to_52 + fa;
      29             :   if (fa >= two_to_52) {
      30             :     r = a;
      31             :   } else {
      32             :     r = r - two_to_52;
      33             :     if (a < 0) r = -r;
      34             :   }
      35             :   return r;
      36             : #endif
      37             : }
      38             : 
      39             : /* default quality ratio for LLL */
      40             : static const double LLLDFT = 0.99;
      41             : 
      42             : /* assume flag & (LLL_KER|LLL_IM|LLL_ALL). LLL_INPLACE implies LLL_IM */
      43             : static GEN
      44      191958 : lll_trivial(GEN x, long flag)
      45             : {
      46      191958 :   if (lg(x) == 1)
      47             :   { /* dim x = 0 */
      48       15443 :     if (! (flag & LLL_ALL)) return cgetg(1,t_MAT);
      49          28 :     retmkvec2(cgetg(1,t_MAT), cgetg(1,t_MAT));
      50             :   }
      51             :   /* dim x = 1 */
      52      176515 :   if (gequal0(gel(x,1)))
      53             :   {
      54          91 :     if (flag & LLL_KER) return matid(1);
      55          91 :     if (flag & (LLL_IM|LLL_INPLACE)) return cgetg(1,t_MAT);
      56          28 :     retmkvec2(matid(1), cgetg(1,t_MAT));
      57             :   }
      58      176419 :   if (flag & LLL_INPLACE) return gcopy(x);
      59       73191 :   if (flag & LLL_KER) return cgetg(1,t_MAT);
      60       73191 :   if (flag & LLL_IM)  return matid(1);
      61          28 :   retmkvec2(cgetg(1,t_MAT), (flag & LLL_GRAM)? gcopy(x): matid(1));
      62             : }
      63             : 
      64             : /* vecslice(x,#x-k,#x) in place. Works for t_MAT, t_VEC/t_COL */
      65             : static GEN
      66     3256173 : vectail_inplace(GEN x, long k)
      67             : {
      68     3256173 :   if (!k) return x;
      69       57552 :   x[k] = ((ulong)x[0] & ~LGBITS) | evallg(lg(x) - k);
      70       57553 :   return x + k;
      71             : }
      72             : 
      73             : /* k = dim Kernel */
      74             : static GEN
      75     3329196 : lll_finish(GEN h, long k, long flag)
      76             : {
      77             :   GEN g;
      78     3329196 :   if (!(flag & (LLL_IM|LLL_KER|LLL_ALL|LLL_INPLACE))) return h;
      79     3256199 :   if (flag & (LLL_IM|LLL_INPLACE)) return vectail_inplace(h, k);
      80          96 :   if (flag & LLL_KER) { setlg(h,k+1); return h; }
      81          68 :   g = vecslice(h,1,k); /* done first: vectail_inplace kills h */
      82          70 :   return mkvec2(g, vectail_inplace(h, k));
      83             : }
      84             : 
      85             : /* y * z * 2^e, e >= 0; y,z t_INT */
      86             : INLINE GEN
      87     1926715 : mulshift(GEN y, GEN z, long e)
      88             : {
      89     1926715 :   long ly = lgefint(y), lz;
      90             :   pari_sp av;
      91             :   GEN t;
      92     1926715 :   if (ly == 2) return gen_0;
      93      204589 :   lz = lgefint(z);
      94      204589 :   av = avma; (void)new_chunk(ly+lz+nbits2lg(e)); /* HACK */
      95      204589 :   t = mulii(z, y);
      96      204589 :   set_avma(av); return shifti(t, e);
      97             : }
      98             : 
      99             : /* x - y * z * 2^e, e >= 0; x,y,z t_INT */
     100             : INLINE GEN
     101     6817078 : submulshift(GEN x, GEN y, GEN z, long e)
     102             : {
     103     6817078 :   long lx = lgefint(x), ly, lz;
     104             :   pari_sp av;
     105             :   GEN t;
     106     6817078 :   if (!e) return submulii(x, y, z);
     107     6786430 :   if (lx == 2) { t = mulshift(y, z, e); togglesign(t); return t; }
     108     5034992 :   ly = lgefint(y);
     109     5034992 :   if (ly == 2) return icopy(x);
     110     4556509 :   lz = lgefint(z);
     111     4556509 :   av = avma; (void)new_chunk(lx+ly+lz+nbits2lg(e)); /* HACK */
     112     4556509 :   t = shifti(mulii(z, y), e);
     113     4556509 :   set_avma(av); return subii(x, t);
     114             : }
     115             : static void
     116   211509099 : subzi(GEN *a, GEN b)
     117             : {
     118   211509099 :   pari_sp av = avma;
     119   211509099 :   b = subii(*a, b);
     120   211508951 :   if (lgefint(b)<=lg(*a) && isonstack(*a)) { affii(b,*a); set_avma(av); }
     121    10964882 :   else *a = b;
     122   211509060 : }
     123             : 
     124             : static void
     125   209919765 : addzi(GEN *a, GEN b)
     126             : {
     127   209919765 :   pari_sp av = avma;
     128   209919765 :   b = addii(*a, b);
     129   209919682 :   if (lgefint(b)<=lg(*a) && isonstack(*a)) { affii(b,*a); set_avma(av); }
     130    10743504 :   else *a = b;
     131   209919798 : }
     132             : 
     133             : /* x - u*y * 2^e */
     134             : INLINE GEN
     135    20099554 : submuliu2n(GEN x, GEN y, ulong u, long e)
     136             : {
     137             :   pari_sp av;
     138    20099554 :   long ly = lgefint(y);
     139    20099554 :   if (ly == 2) return x;
     140    13881631 :   av = avma;
     141    13881631 :   (void)new_chunk(3+ly+lgefint(x)+nbits2lg(e)); /* HACK */
     142    13882681 :   y = shifti(mului(u,y), e);
     143    13881940 :   set_avma(av); return subii(x, y);
     144             : }
     145             : /* *x -= u*y * 2^e */
     146             : INLINE void
     147     3451235 : submulzu2n(GEN *x, GEN y, ulong u, long e)
     148             : {
     149             :   pari_sp av;
     150     3451235 :   long ly = lgefint(y);
     151     3451235 :   if (ly == 2) return;
     152     2655728 :   av = avma;
     153     2655728 :   (void)new_chunk(3+ly+lgefint(*x)+nbits2lg(e)); /* HACK */
     154     2655728 :   y = shifti(mului(u,y), e);
     155     2655728 :   set_avma(av); return subzi(x, y);
     156             : }
     157             : 
     158             : /* x + u*y * 2^e */
     159             : INLINE GEN
     160    20585373 : addmuliu2n(GEN x, GEN y, ulong u, long e)
     161             : {
     162             :   pari_sp av;
     163    20585373 :   long ly = lgefint(y);
     164    20585373 :   if (ly == 2) return x;
     165    14143170 :   av = avma;
     166    14143170 :   (void)new_chunk(3+ly+lgefint(x)+nbits2lg(e)); /* HACK */
     167    14144185 :   y = shifti(mului(u,y), e);
     168    14143451 :   set_avma(av); return addii(x, y);
     169             : }
     170             : 
     171             : /* *x += u*y * 2^e */
     172             : INLINE void
     173     3258874 : addmulzu2n(GEN *x, GEN y, ulong u, long e)
     174             : {
     175             :   pari_sp av;
     176     3258874 :   long ly = lgefint(y);
     177     3258874 :   if (ly == 2) return;
     178     2530508 :   av = avma;
     179     2530508 :   (void)new_chunk(3+ly+lgefint(*x)+nbits2lg(e)); /* HACK */
     180     2530508 :   y = shifti(mului(u,y), e);
     181     2530508 :   set_avma(av); return addzi(x, y);
     182             : }
     183             : 
     184             : /* n < 10; gerepileall supporting &NULL arguments. Maybe rename and export ? */
     185             : INLINE void
     186        9425 : gc_lll(pari_sp av, int n, ...)
     187             : {
     188             :   int i, j;
     189             :   GEN *gptr[10];
     190             :   size_t s;
     191        9425 :   va_list a; va_start(a, n);
     192       28275 :   for (i=j=0; i<n; i++)
     193             :   {
     194       18850 :     GEN *x = va_arg(a,GEN*);
     195       18850 :     if (*x) { gptr[j++] = x; *x = (GEN)copy_bin(*x); }
     196             :   }
     197        9425 :   va_end(a); set_avma(av);
     198       21602 :   for (--j; j>=0; j--) *gptr[j] = bin_copy((GENbin*)*gptr[j]);
     199        9425 :   s = pari_mainstack->top - pari_mainstack->bot;
     200             :   /* size of saved objects ~ stacksize / 4 => overflow */
     201        9425 :   if (av - avma > (s >> 2))
     202             :   {
     203           0 :     size_t t = avma - pari_mainstack->bot;
     204           0 :     av = avma; new_chunk((s + t) / sizeof(long)); set_avma(av); /* double */
     205             :   }
     206        9425 : }
     207             : 
     208             : /********************************************************************/
     209             : /**                                                                **/
     210             : /**                   FPLLL (adapted from D. Stehle's code)        **/
     211             : /**                                                                **/
     212             : /********************************************************************/
     213             : /* Babai* and fplll* are a conversion to libpari API and data types
     214             :    of fplll-1.3 by Damien Stehle'.
     215             : 
     216             :   Copyright 2005, 2006 Damien Stehle'.
     217             : 
     218             :   This program is free software; you can redistribute it and/or modify it
     219             :   under the terms of the GNU General Public License as published by the
     220             :   Free Software Foundation; either version 2 of the License, or (at your
     221             :   option) any later version.
     222             : 
     223             :   This program implements ideas from the paper "Floating-point LLL Revisited",
     224             :   by Phong Nguyen and Damien Stehle', in the Proceedings of Eurocrypt'2005,
     225             :   Springer-Verlag; and was partly inspired by Shoup's NTL library:
     226             :   http://www.shoup.net/ntl/ */
     227             : 
     228             : /* x t_REAL, |x| >= 1/2. Test whether |x| <= 3/2 */
     229             : static int
     230     1024349 : absrsmall2(GEN x)
     231             : {
     232     1024349 :   long e = expo(x), l, i;
     233     1024349 :   if (e < 0) return 1;
     234      325585 :   if (e > 0 || (ulong)x[2] > (3UL << (BITS_IN_LONG-2))) return 0;
     235             :   /* line above assumes l > 2. OK since x != 0 */
     236      167104 :   l = lg(x); for (i = 3; i < l; i++) if (x[i]) return 0;
     237      152436 :   return 1;
     238             : }
     239             : /* x t_REAL; test whether |x| <= 1/2 */
     240             : static int
     241     2249244 : absrsmall(GEN x)
     242             : {
     243             :   long e, l, i;
     244     2249244 :   if (!signe(x)) return 1;
     245     2237422 :   e = expo(x); if (e < -1) return 1;
     246     1030721 :   if (e > -1 || (ulong)x[2] > HIGHBIT) return 0;
     247        7730 :   l = lg(x); for (i = 3; i < l; i++) if (x[i]) return 0;
     248        6372 :   return 1;
     249             : }
     250             : 
     251             : static void
     252    46427626 : rotate(GEN A, long k2, long k)
     253             : {
     254             :   long i;
     255    46427626 :   GEN B = gel(A,k2);
     256   143446722 :   for (i = k2; i > k; i--) gel(A,i) = gel(A,i-1);
     257    46427626 :   gel(A,k) = B;
     258    46427626 : }
     259             : 
     260             : /************************* FAST version (double) ************************/
     261             : #define dmael(x,i,j) ((x)[i][j])
     262             : #define del(x,i) ((x)[i])
     263             : 
     264             : static double *
     265    42238881 : cget_dblvec(long d)
     266    42238881 : { return (double*) stack_malloc_align(d*sizeof(double), sizeof(double)); }
     267             : 
     268             : static double **
     269    12986223 : cget_dblmat(long d) { return (double **) cgetg(d, t_VECSMALL); }
     270             : 
     271             : static double
     272   254262802 : itodbl_exp(GEN x, long *e)
     273             : {
     274   254262802 :   pari_sp av = avma;
     275   254262802 :   GEN r = itor(x,DEFAULTPREC);
     276   254231042 :   *e = expo(r); setexpo(r,0);
     277   254227712 :   return gc_double(av, rtodbl(r));
     278             : }
     279             : 
     280             : static double
     281   178606900 : dbldotproduct(double *x, double *y, long n)
     282             : {
     283             :   long i;
     284   178606900 :   double sum = del(x,1) * del(y,1);
     285  2877191432 :   for (i=2; i<=n; i++) sum += del(x,i) * del(y,i);
     286   178606900 :   return sum;
     287             : }
     288             : 
     289             : static double
     290     3608776 : dbldotsquare(double *x, long n)
     291             : {
     292             :   long i;
     293     3608776 :   double sum = del(x,1) * del(x,1);
     294     9433299 :   for (i=2; i<=n; i++) sum += del(x,i) * del(x,i);
     295     3608776 :   return sum;
     296             : }
     297             : 
     298             : static long
     299    33522812 : set_line(double *appv, GEN v, long n)
     300             : {
     301    33522812 :   long i, maxexp = 0;
     302    33522812 :   pari_sp av = avma;
     303    33522812 :   GEN e = cgetg(n+1, t_VECSMALL);
     304   287772459 :   for (i = 1; i <= n; i++)
     305             :   {
     306   254261535 :     del(appv,i) = itodbl_exp(gel(v,i), e+i);
     307   254248326 :     if (e[i] > maxexp) maxexp = e[i];
     308             :   }
     309   287850545 :   for (i = 1; i <= n; i++) del(appv,i) = ldexp(del(appv,i), e[i]-maxexp);
     310    33510924 :   set_avma(av); return maxexp;
     311             : }
     312             : 
     313             : static void
     314    51163542 : dblrotate(double **A, long k2, long k)
     315             : {
     316             :   long i;
     317    51163542 :   double *B = del(A,k2);
     318   158447476 :   for (i = k2; i > k; i--) del(A,i) = del(A,i-1);
     319    51163542 :   del(A,k) = B;
     320    51163542 : }
     321             : /* update G[kappa][i] from appB */
     322             : static void
     323    30117569 : setG_fast(double **appB, long n, double **G, long kappa, long a, long b)
     324             : { long i;
     325   144851854 :   for (i = a; i <= b; i++)
     326   114734555 :     dmael(G,kappa,i) = dbldotproduct(del(appB,kappa), del(appB,i), n);
     327    30117299 : }
     328             : /* update G[i][kappa] from appB */
     329             : static void
     330    24251753 : setG2_fast(double **appB, long n, double **G, long kappa, long a, long b)
     331             : { long i;
     332    88128071 :   for (i = a; i <= b; i++)
     333    63876403 :     dmael(G,i,kappa) = dbldotproduct(del(appB,kappa), del(appB,i), n);
     334    24251668 : }
     335             : const long EX0 = -2; /* uninitialized; any value less than expo(0.51) = -1 */
     336             : 
     337             : #ifdef LONG_IS_64BIT
     338             : typedef long s64;
     339             : #define addmuliu64_inplace addmuliu_inplace
     340             : #define submuliu64_inplace submuliu_inplace
     341             : #define submuliu642n submuliu2n
     342             : #define addmuliu642n addmuliu2n
     343             : #else
     344             : typedef long long s64;
     345             : typedef unsigned long long u64;
     346             : 
     347             : INLINE GEN
     348    33873301 : u64toi(u64 x)
     349             : {
     350             :   GEN y;
     351             :   ulong h;
     352    33873301 :   if (!x) return gen_0;
     353    33873301 :   h = x>>32;
     354    33873301 :   if (!h) return utoipos(x);
     355     3987690 :   y = cgetipos(4);
     356     3987690 :   *int_LSW(y) = x&0xFFFFFFFF;
     357     3987690 :   *int_MSW(y) = x>>32;
     358     3987690 :   return y;
     359             : }
     360             : 
     361             : INLINE GEN
     362     3195525 : u64toineg(u64 x)
     363             : {
     364             :   GEN y;
     365             :   ulong h;
     366     3195525 :   if (!x) return gen_0;
     367     3195525 :   h = x>>32;
     368     3195525 :   if (!h) return utoineg(x);
     369     3195525 :   y = cgetineg(4);
     370     3195525 :   *int_LSW(y) = x&0xFFFFFFFF;
     371     3195525 :   *int_MSW(y) = x>>32;
     372     3195525 :   return y;
     373             : }
     374             : INLINE GEN
     375    15253786 : addmuliu64_inplace(GEN x, GEN y, u64 u) { return addmulii(x, y, u64toi(u)); }
     376             : 
     377             : INLINE GEN
     378    15316291 : submuliu64_inplace(GEN x, GEN y, u64 u) { return submulii(x, y, u64toi(u)); }
     379             : 
     380             : INLINE GEN
     381     3195525 : addmuliu642n(GEN x, GEN y, u64 u, long e) { return submulshift(x, y, u64toineg(u), e); }
     382             : 
     383             : INLINE GEN
     384     3303224 : submuliu642n(GEN x, GEN y, u64 u, long e) { return submulshift(x, y, u64toi(u), e); }
     385             : 
     386             : #endif
     387             : 
     388             : /* Babai's Nearest Plane algorithm (iterative); see Babai() */
     389             : static int
     390    41081326 : Babai_fast(pari_sp av, long kappa, GEN *pB, GEN *pU, double **mu, double **r,
     391             :            double *s, double **appB, GEN expoB, double **G,
     392             :            long a, long zeros, long maxG, double eta)
     393             : {
     394    41081326 :   GEN B = *pB, U = *pU;
     395    41081326 :   const long n = nbrows(B), d = U ? lg(U)-1: 0;
     396    41081127 :   long k, aa = (a > zeros)? a : zeros+1;
     397    41081127 :   long emaxmu = EX0, emax2mu = EX0;
     398             :   s64 xx;
     399    41081127 :   int did_something = 0;
     400             :   /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
     401             : 
     402    24588596 :   for (;;) {
     403    65669723 :     int go_on = 0;
     404    65669723 :     long i, j, emax3mu = emax2mu;
     405             : 
     406    65669723 :     if (gc_needed(av,2))
     407             :     {
     408        1521 :       if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
     409        1521 :       gc_lll(av,2,&B,&U);
     410             :     }
     411             :     /* Step2: compute the GSO for stage kappa */
     412    65668888 :     emax2mu = emaxmu; emaxmu = EX0;
     413   269159926 :     for (j=aa; j<kappa; j++)
     414             :     {
     415   203491892 :       double g = dmael(G,kappa,j);
     416  1100248365 :       for (k = zeros+1; k < j; k++) g -= dmael(mu,j,k) * dmael(r,kappa,k);
     417   203491892 :       dmael(r,kappa,j) = g;
     418   203491892 :       dmael(mu,kappa,j) = dmael(r,kappa,j) / dmael(r,j,j);
     419   203491892 :       emaxmu = maxss(emaxmu, expoB[kappa]-expoB[j]);
     420             :     }
     421             :     /* maxmu doesn't decrease fast enough */
     422    65668034 :     if (emax3mu != EX0 && emax3mu <= emax2mu + 5) {*pB = B; *pU = U; return 1;}
     423             : 
     424   262072203 :     for (j=kappa-1; j>zeros; j--)
     425             :     {
     426   220996385 :       double tmp = fabs(ldexp (dmael(mu,kappa,j), expoB[kappa]-expoB[j]));
     427   220996385 :       if (tmp>eta) { go_on = 1; break; }
     428             :     }
     429             : 
     430             :     /* Step3--5: compute the X_j's  */
     431    65664503 :     if (go_on)
     432   122952281 :       for (j=kappa-1; j>zeros; j--)
     433             :       { /* The code below seemingly handles U = NULL, but in this case d = 0 */
     434    98372377 :         int e = expoB[j] - expoB[kappa];
     435    98372377 :         double tmp = ldexp(dmael(mu,kappa,j), -e), atmp = fabs(tmp);
     436             :         /* tmp = Inf is allowed */
     437    98372377 :         if (atmp <= .5) continue; /* size-reduced */
     438    54314714 :         if (gc_needed(av,2))
     439             :         {
     440        4279 :           if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
     441        4279 :           gc_lll(av,2,&B,&U);
     442             :         }
     443    54316683 :         did_something = 1;
     444             :         /* we consider separately the case |X| = 1 */
     445    54316683 :         if (atmp <= 1.5)
     446             :         {
     447    36134798 :           if (dmael(mu,kappa,j) > 0) { /* in this case, X = 1 */
     448    94912123 :             for (k=zeros+1; k<j; k++)
     449    76354026 :               dmael(mu,kappa,k) -= ldexp(dmael(mu,j,k), e);
     450   342736652 :             for (i=1; i<=n; i++)
     451   324183390 :               gmael(B,kappa,i) = subii(gmael(B,kappa,i), gmael(B,j,i));
     452   113185711 :             for (i=1; i<=d; i++)
     453    94632283 :               gmael(U,kappa,i) = subii(gmael(U,kappa,i), gmael(U,j,i));
     454             :           } else { /* otherwise X = -1 */
     455    93272326 :             for (k=zeros+1; k<j; k++)
     456    75695625 :               dmael(mu,kappa,k) += ldexp(dmael(mu,j,k), e);
     457   337773779 :             for (i=1; i<=n; i++)
     458   320202268 :               gmael(B,kappa,i) = addii(gmael(B,kappa,i), gmael(B,j,i));
     459   108659978 :             for (i=1; i<=d; i++)
     460    91088315 :               gmael(U,kappa,i) = addii(gmael(U,kappa,i), gmael(U,j,i));
     461             :           }
     462    36125091 :           continue;
     463             :         }
     464             :         /* we have |X| >= 2 */
     465    18181885 :         if (atmp < 9007199254740992.)
     466             :         {
     467    15741167 :           tmp = pari_rint(tmp);
     468    38347965 :           for (k=zeros+1; k<j; k++)
     469    22606803 :             dmael(mu,kappa,k) -= ldexp(tmp * dmael(mu,j,k), e);
     470    15741162 :           xx = (s64) tmp;
     471    15741162 :           if (xx > 0) /* = xx */
     472             :           {
     473    80178977 :             for (i=1; i<=n; i++)
     474    72277238 :               gmael(B,kappa,i) = submuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
     475    44941686 :             for (i=1; i<=d; i++)
     476    37039785 :               gmael(U,kappa,i) = submuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
     477             :           }
     478             :           else /* = -xx */
     479             :           {
     480    79660851 :             for (i=1; i<=n; i++)
     481    71822364 :               gmael(B,kappa,i) = addmuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), -xx);
     482    44281641 :             for (i=1; i<=d; i++)
     483    36443075 :               gmael(U,kappa,i) = addmuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), -xx);
     484             :           }
     485             :         }
     486             :         else
     487             :         {
     488             :           int E;
     489     2440718 :           xx = (s64) ldexp(frexp(dmael(mu,kappa,j), &E), 53);
     490     2440718 :           E -= e + 53;
     491     2440718 :           if (E <= 0)
     492             :           {
     493           0 :             xx = xx << -E;
     494           0 :             for (k=zeros+1; k<j; k++)
     495           0 :               dmael(mu,kappa,k) -= ldexp(((double)xx) * dmael(mu,j,k), e);
     496           0 :             if (xx > 0) /* = xx */
     497             :             {
     498           0 :               for (i=1; i<=n; i++)
     499           0 :                 gmael(B,kappa,i) = submuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
     500           0 :               for (i=1; i<=d; i++)
     501           0 :                 gmael(U,kappa,i) = submuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
     502             :             }
     503             :             else /* = -xx */
     504             :             {
     505           0 :               for (i=1; i<=n; i++)
     506           0 :                 gmael(B,kappa,i) = addmuliu64_inplace(gmael(B,kappa,i), gmael(B,j,i), -xx);
     507           0 :               for (i=1; i<=d; i++)
     508           0 :                 gmael(U,kappa,i) = addmuliu64_inplace(gmael(U,kappa,i), gmael(U,j,i), -xx);
     509             :             }
     510             :           } else
     511             :           {
     512    16082302 :             for (k=zeros+1; k<j; k++)
     513    13641584 :               dmael(mu,kappa,k) -= ldexp(((double)xx) * dmael(mu,j,k), E + e);
     514     2440718 :             if (xx > 0) /* = xx */
     515             :             {
     516    22970848 :               for (i=1; i<=n; i++)
     517    21753881 :                 gmael(B,kappa,i) = submuliu642n(gmael(B,kappa,i), gmael(B,j,i), xx, E);
     518     2112723 :               for (i=1; i<=d; i++)
     519      895756 :                 gmael(U,kappa,i) = submuliu642n(gmael(U,kappa,i), gmael(U,j,i), xx, E);
     520             :             }
     521             :             else /* = -xx */
     522             :             {
     523    23414147 :               for (i=1; i<=n; i++)
     524    22190756 :                 gmael(B,kappa,i) = addmuliu642n(gmael(B,kappa,i), gmael(B,j,i), -xx, E);
     525     2122694 :               for (i=1; i<=d; i++)
     526      899303 :                 gmael(U,kappa,i) = addmuliu642n(gmael(U,kappa,i), gmael(U,j,i), -xx, E);
     527             :             }
     528             :           }
     529             :         }
     530             :       }
     531    65655705 :     if (!go_on) break; /* Anything happened? */
     532    24578340 :     expoB[kappa] = set_line(del(appB,kappa), gel(B,kappa), n);
     533    24588286 :     setG_fast(appB, n, G, kappa, zeros+1, kappa-1);
     534    24588596 :     aa = zeros+1;
     535             :   }
     536    41077365 :   if (did_something) setG2_fast(appB, n, G, kappa, kappa, maxG);
     537             : 
     538    41077702 :   del(s,zeros+1) = dmael(G,kappa,kappa);
     539             :   /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
     540   179805316 :   for (k=zeros+1; k<=kappa-2; k++)
     541   138727614 :     del(s,k+1) = del(s,k) - dmael(mu,kappa,k)*dmael(r,kappa,k);
     542    41077702 :   *pB = B; *pU = U; return 0;
     543             : }
     544             : 
     545             : static void
     546    18788861 : update_alpha(GEN alpha, long kappa, long kappa2, long kappamax)
     547             : {
     548             :   long i;
     549    59957608 :   for (i = kappa; i < kappa2; i++)
     550    41168747 :     if (kappa <= alpha[i]) alpha[i] = kappa;
     551    59957668 :   for (i = kappa2; i > kappa; i--) alpha[i] = alpha[i-1];
     552    49195826 :   for (i = kappa2+1; i <= kappamax; i++)
     553    30406965 :     if (kappa < alpha[i]) alpha[i] = kappa;
     554    18788861 :   alpha[kappa] = kappa;
     555    18788861 : }
     556             : static void
     557     1734202 : rotateG(GEN G, long kappa2, long kappa, long maxG, GEN Gtmp)
     558             : {
     559             :   long i, j;
     560    20879724 :   for (i=1; i<=kappa2; i++) gel(Gtmp,i) = gmael(G,kappa2,i);
     561    10446183 :   for (   ; i<=maxG; i++)   gel(Gtmp,i) = gmael(G,i,kappa2);
     562     7141248 :   for (i=kappa2; i>kappa; i--)
     563             :     {
     564    37921853 :       for (j=1; j<kappa; j++) gmael(G,i,j) = gmael(G,i-1,j);
     565     5407046 :       gmael(G,i,kappa) = gel(Gtmp,i-1);
     566    25948538 :       for (j=kappa+1; j<=i; j++) gmael(G,i,j) = gmael(G,i-1,j-1);
     567    26940949 :       for (j=kappa2+1; j<=maxG; j++) gmael(G,j,i) = gmael(G,j,i-1);
     568             :     }
     569    13738476 :   for (i=1; i<kappa; i++) gmael(G,kappa,i) = gel(Gtmp,i);
     570     1734202 :   gmael(G,kappa,kappa) = gel(Gtmp,kappa2);
     571    10446183 :   for (i=kappa2+1; i<=maxG; i++) gmael(G,i,kappa) = gel(Gtmp,i);
     572     1734202 : }
     573             : static void
     574    17054580 : rotateG_fast(double **G, long kappa2, long kappa, long maxG, double *Gtmp)
     575             : {
     576             :   long i, j;
     577   103939611 :   for (i=1; i<=kappa2; i++) del(Gtmp,i) = dmael(G,kappa2,i);
     578    39734749 :   for (   ; i<=maxG; i++) del(Gtmp,i) = dmael(G,i,kappa2);
     579    52816149 :   for (i=kappa2; i>kappa; i--)
     580             :   {
     581   119435737 :     for (j=1; j<kappa; j++) dmael(G,i,j) = dmael(G,i-1,j);
     582    35761569 :     dmael(G,i,kappa) = del(Gtmp,i-1);
     583   128786832 :     for (j=kappa+1; j<=i; j++) dmael(G,i,j) = dmael(G,i-1,j-1);
     584    87068271 :     for (j=kappa2+1; j<=maxG; j++) dmael(G,j,i) = dmael(G,j,i-1);
     585             :   }
     586    51123721 :   for (i=1; i<kappa; i++) dmael(G,kappa,i) = del(Gtmp,i);
     587    17054580 :   dmael(G,kappa,kappa) = del(Gtmp,kappa2);
     588    39734788 :   for (i=kappa2+1; i<=maxG; i++) dmael(G,i,kappa) = del(Gtmp,i);
     589    17054580 : }
     590             : 
     591             : /* LLL-reduces (B,U) in place [apply base change transforms to B and U].
     592             :  * Gram matrix, and GSO performed on matrices of 'double'.
     593             :  * If (keepfirst), never swap with first vector.
     594             :  * Return -1 on failure, else zeros = dim Kernel (>= 0) */
     595             : static long
     596     3246572 : fplll_fast(GEN *pB, GEN *pU, double delta, double eta, long keepfirst)
     597             : {
     598             :   pari_sp av;
     599             :   long kappa, kappa2, d, n, i, j, zeros, kappamax, maxG;
     600             :   double **mu, **r, *s, tmp, *Gtmp, **G, **appB;
     601     3246572 :   GEN alpha, expoB, B = *pB, U;
     602     3246572 :   long cnt = 0;
     603             : 
     604     3246572 :   d = lg(B)-1;
     605     3246572 :   n = nbrows(B);
     606     3246569 :   U = *pU; /* NULL if inplace */
     607             : 
     608     3246569 :   G = cget_dblmat(d+1);
     609     3246570 :   appB = cget_dblmat(d+1);
     610     3246568 :   mu = cget_dblmat(d+1);
     611     3246565 :   r  = cget_dblmat(d+1);
     612     3246568 :   s  = cget_dblvec(d+1);
     613    12183178 :   for (j = 1; j <= d; j++)
     614             :   {
     615     8936593 :     del(mu,j) = cget_dblvec(d+1);
     616     8936604 :     del(r,j) = cget_dblvec(d+1);
     617     8936610 :     del(appB,j) = cget_dblvec(n+1);
     618     8936607 :     del(G,j) = cget_dblvec(d+1);
     619    45451473 :     for (i=1; i<=d; i++) dmael(G,j,i) = 0.;
     620             :   }
     621     3246585 :   expoB = cgetg(d+1, t_VECSMALL);
     622    12183071 :   for (i=1; i<=d; i++) expoB[i] = set_line(del(appB,i), gel(B,i), n);
     623     3246502 :   Gtmp = cget_dblvec(d+1);
     624     3246567 :   alpha = cgetg(d+1, t_VECSMALL);
     625     3246559 :   av = avma;
     626             : 
     627             :   /* Step2: Initializing the main loop */
     628     3246559 :   kappamax = 1;
     629     3246559 :   i = 1;
     630     3246559 :   maxG = d; /* later updated to kappamax */
     631             : 
     632             :   do {
     633     3401101 :     dmael(G,i,i) = dbldotsquare(del(appB,i),n);
     634     3401103 :   } while (dmael(G,i,i) <= 0 && (++i <=d));
     635     3246561 :   zeros = i-1; /* all vectors B[i] with i <= zeros are zero vectors */
     636     3246561 :   kappa = i;
     637     3246561 :   if (zeros < d) dmael(r,zeros+1,zeros+1) = dmael(G,zeros+1,zeros+1);
     638    12028646 :   for (i=zeros+1; i<=d; i++) alpha[i]=1;
     639    44324437 :   while (++kappa <= d)
     640             :   {
     641    41081392 :     if (kappa > kappamax)
     642             :     {
     643     5529353 :       if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
     644     5529353 :       maxG = kappamax = kappa;
     645     5529353 :       setG_fast(appB, n, G, kappa, zeros+1, kappa);
     646             :     }
     647             :     /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
     648    41081401 :     if (Babai_fast(av, kappa, &B,&U, mu,r,s, appB, expoB, G, alpha[kappa],
     649        3531 :                    zeros, maxG, eta)) { *pB=B; *pU=U; return -1; }
     650             : 
     651    41077642 :     tmp = ldexp(r[kappa-1][kappa-1] * delta, 2*(expoB[kappa-1]-expoB[kappa]));
     652    41077642 :     if ((keepfirst && kappa == 2) || tmp <= del(s,kappa-1))
     653             :     { /* Step4: Success of Lovasz's condition */
     654    24022980 :       alpha[kappa] = kappa;
     655    24022980 :       tmp = dmael(mu,kappa,kappa-1) * dmael(r,kappa,kappa-1);
     656    24022980 :       dmael(r,kappa,kappa) = del(s,kappa-1)- tmp;
     657    24022980 :       continue;
     658             :     }
     659             :     /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
     660    17054662 :     if (DEBUGLEVEL>=4 && kappa==kappamax && del(s,kappa-1)!=0)
     661           0 :       if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", 2*expoB[1] + dblexpo(del(s,1))); }
     662    17054645 :     kappa2 = kappa;
     663             :     do {
     664    35761678 :       kappa--;
     665    35761678 :       if (kappa<zeros+2 + (keepfirst ? 1: 0)) break;
     666    25447016 :       tmp = dmael(r,kappa-1,kappa-1) * delta;
     667    25447016 :       tmp = ldexp(tmp, 2*(expoB[kappa-1]-expoB[kappa2]));
     668    25447016 :     } while (del(s,kappa-1) <= tmp);
     669    17054645 :     update_alpha(alpha, kappa, kappa2, kappamax);
     670             : 
     671             :     /* Step6: Update the mu's and r's */
     672    17054687 :     dblrotate(mu,kappa2,kappa);
     673    17054642 :     dblrotate(r,kappa2,kappa);
     674    17054605 :     dmael(r,kappa,kappa) = del(s,kappa);
     675             : 
     676             :     /* Step7: Update B, appB, U, G */
     677    17054605 :     rotate(B,kappa2,kappa);
     678    17054596 :     dblrotate(appB,kappa2,kappa);
     679    17054566 :     if (U) rotate(U,kappa2,kappa);
     680    17054564 :     rotate(expoB,kappa2,kappa);
     681    17054575 :     rotateG_fast(G,kappa2,kappa, maxG, Gtmp);
     682             : 
     683             :     /* Step8: Prepare the next loop iteration */
     684    17054896 :     if (kappa == zeros+1 && dmael(G,kappa,kappa)<= 0)
     685             :     {
     686      207677 :       zeros++; kappa++;
     687      207677 :       dmael(G,kappa,kappa) = dbldotsquare(del(appB,kappa),n);
     688      207677 :       dmael(r,kappa,kappa) = dmael(G,kappa,kappa);
     689             :     }
     690             :   }
     691     3243045 :   *pB = B; *pU = U; return zeros;
     692             : }
     693             : 
     694             : /***************** HEURISTIC version (reduced precision) ****************/
     695             : static GEN
     696      329886 : realsqrdotproduct(GEN x)
     697             : {
     698      329886 :   long i, l = lg(x);
     699      329886 :   GEN z = sqrr(gel(x,1));
     700     4024724 :   for (i=2; i<l; i++) z = addrr(z, sqrr(gel(x,i)));
     701      329886 :   return z;
     702             : }
     703             : /* x, y non-empty vector of t_REALs, same length */
     704             : static GEN
     705     3656008 : realdotproduct(GEN x, GEN y)
     706             : {
     707             :   long i, l;
     708             :   GEN z;
     709     3656008 :   if (x == y) return realsqrdotproduct(x);
     710     3326122 :   l = lg(x); z = mulrr(gel(x,1),gel(y,1));
     711    52976165 :   for (i=2; i<l; i++) z = addrr(z, mulrr(gel(x,i), gel(y,i)));
     712     3326122 :   return z;
     713             : }
     714             : static void
     715      337385 : setG_heuristic(GEN appB, GEN G, long kappa, long a, long b)
     716      337385 : { pari_sp av = avma;
     717             :   long i;
     718     2634186 :   for (i = a; i <= b; i++)
     719     2296801 :     affrr(realdotproduct(gel(appB,kappa),gel(appB,i)), gmael(G,kappa,i));
     720      337385 :   set_avma(av);
     721      337385 : }
     722             : static void
     723      318663 : setG2_heuristic(GEN appB, GEN G, long kappa, long a, long b)
     724      318663 : { pari_sp av = avma;
     725             :   long i;
     726     1677870 :   for (i = a; i <= b; i++)
     727     1359207 :     affrr(realdotproduct(gel(appB,kappa),gel(appB,i)), gmael(G,i,kappa));
     728      318663 :   set_avma(av);
     729      318663 : }
     730             : 
     731             : /* approximate t_REAL x as m * 2^e, where |m| < 2^bit */
     732             : static GEN
     733       20439 : truncexpo(GEN x, long bit, long *e)
     734             : {
     735       20439 :   *e = expo(x) + 1 - bit;
     736       20439 :   if (*e >= 0) return mantissa2nr(x, 0);
     737        1787 :   *e = 0; return roundr_safe(x);
     738             : }
     739             : /* Babai's Nearest Plane algorithm (iterative); see Babai() */
     740             : static int
     741      655382 : Babai_heuristic(pari_sp av, long kappa, GEN *pB, GEN *pU, GEN mu, GEN r, GEN s,
     742             :                 GEN appB, GEN G, long a, long zeros, long maxG,
     743             :                 GEN eta, long prec)
     744             : {
     745      655382 :   GEN B = *pB, U = *pU;
     746      655382 :   const long n = nbrows(B), d = U ? lg(U)-1: 0, bit = prec2nbits(prec);
     747      655382 :   long k, aa = (a > zeros)? a : zeros+1;
     748      655382 :   int did_something = 0;
     749      655382 :   long emaxmu = EX0, emax2mu = EX0;
     750             :   /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
     751             : 
     752      326162 :   for (;;) {
     753      981544 :     int go_on = 0;
     754      981544 :     long i, j, emax3mu = emax2mu;
     755             : 
     756      981544 :     if (gc_needed(av,2))
     757             :     {
     758           6 :       if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
     759           6 :       gc_lll(av,2,&B,&U);
     760             :     }
     761             :     /* Step2: compute the GSO for stage kappa */
     762      981544 :     emax2mu = emaxmu; emaxmu = EX0;
     763     5446427 :     for (j=aa; j<kappa; j++)
     764             :     {
     765     4464883 :       pari_sp btop = avma;
     766     4464883 :       GEN g = gmael(G,kappa,j);
     767    27941574 :       for (k = zeros+1; k<j; k++)
     768    23476691 :         g = subrr(g, mulrr(gmael(mu,j,k), gmael(r,kappa,k)));
     769     4464883 :       affrr(g, gmael(r,kappa,j));
     770     4464883 :       affrr(divrr(gmael(r,kappa,j), gmael(r,j,j)), gmael(mu,kappa,j));
     771     4464883 :       emaxmu = maxss(emaxmu, expo(gmael(mu,kappa,j)));
     772     4464883 :       set_avma(btop);
     773             :     }
     774      981544 :     if (emax3mu != EX0 && emax3mu <= emax2mu + 5)
     775        1177 :     { *pB = B; *pU = U; return 1; }
     776             : 
     777     5932720 :     for (j=kappa-1; j>zeros; j--)
     778     5278515 :       if (abscmprr(gmael(mu,kappa,j), eta) > 0) { go_on = 1; break; }
     779             : 
     780             :     /* Step3--5: compute the X_j's  */
     781      980367 :     if (go_on)
     782     2575406 :       for (j=kappa-1; j>zeros; j--)
     783             :       { /* The code below seemingly handles U = NULL, but in this case d = 0 */
     784             :         pari_sp btop;
     785     2249244 :         GEN tmp = gmael(mu,kappa,j);
     786     2249244 :         if (absrsmall(tmp)) continue; /* size-reduced */
     787             : 
     788     1024349 :         if (gc_needed(av,2))
     789             :         {
     790          88 :           if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
     791          88 :           gc_lll(av,2,&B,&U);
     792             :         }
     793     1024349 :         btop = avma; did_something = 1;
     794             :         /* we consider separately the case |X| = 1 */
     795     1024349 :         if (absrsmall2(tmp))
     796             :         {
     797      851200 :           if (signe(tmp) > 0) { /* in this case, X = 1 */
     798     2381987 :             for (k=zeros+1; k<j; k++)
     799     1956009 :               affrr(subrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
     800      425978 :             set_avma(btop);
     801     7408988 :             for (i=1; i<=n; i++)
     802     6983010 :               gmael(B,kappa,i) = subii(gmael(B,kappa,i), gmael(B,j,i));
     803     6194282 :             for (i=1; i<=d; i++)
     804     5768304 :               gmael(U,kappa,i) = subii(gmael(U,kappa,i), gmael(U,j,i));
     805             :           } else { /* otherwise X = -1 */
     806     2378459 :             for (k=zeros+1; k<j; k++)
     807     1953237 :               affrr(addrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
     808      425222 :             set_avma(btop);
     809     7388388 :             for (i=1; i<=n; i++)
     810     6963166 :               gmael(B,kappa,i) = addii(gmael(B,kappa,i), gmael(B,j,i));
     811     6169988 :             for (i=1; i<=d; i++)
     812     5744766 :               gmael(U,kappa,i) = addii(gmael(U,kappa,i),gmael(U,j,i));
     813             :           }
     814      851200 :           continue;
     815             :         }
     816             :         /* we have |X| >= 2 */
     817      173149 :         if (expo(tmp) < BITS_IN_LONG)
     818             :         {
     819      152710 :           ulong xx = roundr_safe(tmp)[2]; /* X fits in an ulong */
     820      152710 :           if (signe(tmp) > 0) /* = xx */
     821             :           {
     822      251334 :             for (k=zeros+1; k<j; k++)
     823      175832 :               affrr(subrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
     824      175832 :                   gmael(mu,kappa,k));
     825       75502 :             set_avma(btop);
     826      912585 :             for (i=1; i<=n; i++)
     827      837083 :               gmael(B,kappa,i) = submuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
     828      722073 :             for (i=1; i<=d; i++)
     829      646571 :               gmael(U,kappa,i) = submuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
     830             :           }
     831             :           else /* = -xx */
     832             :           {
     833      259324 :             for (k=zeros+1; k<j; k++)
     834      182116 :               affrr(addrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
     835      182116 :                   gmael(mu,kappa,k));
     836       77208 :             set_avma(btop);
     837      941899 :             for (i=1; i<=n; i++)
     838      864691 :               gmael(B,kappa,i) = addmuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
     839      743053 :             for (i=1; i<=d; i++)
     840      665845 :               gmael(U,kappa,i) = addmuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
     841             :           }
     842             :         }
     843             :         else
     844             :         {
     845             :           long e;
     846       20439 :           GEN X = truncexpo(tmp, bit, &e); /* tmp ~ X * 2^e */
     847       20439 :           btop = avma;
     848       87782 :           for (k=zeros+1; k<j; k++)
     849             :           {
     850       67343 :             GEN x = mulir(X, gmael(mu,j,k));
     851       67343 :             if (e) shiftr_inplace(x, e);
     852       67343 :             affrr(subrr(gmael(mu,kappa,k), x), gmael(mu,kappa,k));
     853             :           }
     854       20439 :           set_avma(btop);
     855      263178 :           for (i=1; i<=n; i++)
     856      242739 :             gmael(B,kappa,i) = submulshift(gmael(B,kappa,i), gmael(B,j,i), X, e);
     857       96029 :           for (i=1; i<=d; i++)
     858       75590 :             gmael(U,kappa,i) = submulshift(gmael(U,kappa,i), gmael(U,j,i), X, e);
     859             :         }
     860             :       }
     861      980367 :     if (!go_on) break; /* Anything happened? */
     862     4326279 :     for (i=1 ; i<=n; i++) affir(gmael(B,kappa,i), gmael(appB,kappa,i));
     863      326162 :     setG_heuristic(appB, G, kappa, zeros+1, kappa-1);
     864      326162 :     aa = zeros+1;
     865             :   }
     866      654205 :   if (did_something) setG2_heuristic(appB, G, kappa, kappa, maxG);
     867      654205 :   affrr(gmael(G,kappa,kappa), gel(s,zeros+1));
     868             :   /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
     869      654205 :   av = avma;
     870     4521945 :   for (k=zeros+1; k<=kappa-2; k++)
     871     3867740 :     affrr(subrr(gel(s,k), mulrr(gmael(mu,kappa,k), gmael(r,kappa,k))),
     872     3867740 :           gel(s,k+1));
     873      654205 :   *pB = B; *pU = U; return gc_bool(av, 0);
     874             : }
     875             : 
     876             : static GEN
     877       17991 : ZC_to_RC(GEN x, long prec)
     878      175874 : { pari_APPLY_type(t_COL,itor(gel(x,i),prec)) }
     879             : 
     880             : static GEN
     881        3531 : ZM_to_RM(GEN x, long prec)
     882       21522 : { pari_APPLY_same(ZC_to_RC(gel(x,i),prec)) }
     883             : 
     884             : /* LLL-reduces (B,U) in place [apply base change transforms to B and U].
     885             :  * Gram matrix made of t_REAL at precision prec2, performe GSO at prec.
     886             :  * If (keepfirst), never swap with first vector.
     887             :  * Return -1 on failure, else zeros = dim Kernel (>= 0) */
     888             : static long
     889        3531 : fplll_heuristic(GEN *pB, GEN *pU, double DELTA, double ETA, long keepfirst,
     890             :                 long prec, long prec2)
     891             : {
     892             :   pari_sp av, av2;
     893             :   long kappa, kappa2, d, i, j, zeros, kappamax, maxG;
     894        3531 :   GEN mu, r, s, tmp, Gtmp, alpha, G, appB, B = *pB, U;
     895        3531 :   GEN delta = dbltor(DELTA), eta = dbltor(ETA);
     896        3531 :   long cnt = 0;
     897             : 
     898        3531 :   d = lg(B)-1;
     899        3531 :   U = *pU; /* NULL if inplace */
     900             : 
     901        3531 :   G = cgetg(d+1, t_MAT);
     902        3531 :   mu = cgetg(d+1, t_MAT);
     903        3531 :   r  = cgetg(d+1, t_MAT);
     904        3531 :   s  = cgetg(d+1, t_VEC);
     905        3531 :   appB = ZM_to_RM(B, prec2);
     906       21522 :   for (j = 1; j <= d; j++)
     907             :   {
     908       17991 :     GEN M = cgetg(d+1, t_COL), R = cgetg(d+1, t_COL), S = cgetg(d+1, t_COL);
     909       17991 :     gel(mu,j)= M;
     910       17991 :     gel(r,j) = R;
     911       17991 :     gel(G,j) = S;
     912       17991 :     gel(s,j) = cgetr(prec);
     913      167182 :     for (i = 1; i <= d; i++)
     914             :     {
     915      149191 :       gel(R,i) = cgetr(prec);
     916      149191 :       gel(M,i) = cgetr(prec);
     917      149191 :       gel(S,i) = cgetr(prec2);
     918             :     }
     919             :   }
     920        3531 :   Gtmp = cgetg(d+1, t_VEC);
     921        3531 :   alpha = cgetg(d+1, t_VECSMALL);
     922        3531 :   av = avma;
     923             : 
     924             :   /* Step2: Initializing the main loop */
     925        3531 :   kappamax = 1;
     926        3531 :   i = 1;
     927        3531 :   maxG = d; /* later updated to kappamax */
     928             : 
     929             :   do {
     930        3532 :     affrr(RgV_dotsquare(gel(appB,i)), gmael(G,i,i));
     931        3532 :   } while (signe(gmael(G,i,i)) == 0 && (++i <=d));
     932        3531 :   zeros = i-1; /* all vectors B[i] with i <= zeros are zero vectors */
     933        3531 :   kappa = i;
     934        3531 :   if (zeros < d) affrr(gmael(G,zeros+1,zeros+1), gmael(r,zeros+1,zeros+1));
     935       21521 :   for (i=zeros+1; i<=d; i++) alpha[i]=1;
     936             : 
     937      657736 :   while (++kappa <= d)
     938             :   {
     939      655382 :     if (kappa > kappamax)
     940             :     {
     941       11223 :       if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
     942       11223 :       maxG = kappamax = kappa;
     943       11223 :       setG_heuristic(appB, G, kappa, zeros+1, kappa);
     944             :     }
     945             :     /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
     946      655382 :     if (Babai_heuristic(av, kappa, &B,&U, mu,r,s, appB, G, alpha[kappa], zeros,
     947        1177 :                         maxG, eta, prec)) { *pB = B; *pU = U; return -1; }
     948      654205 :     av2 = avma;
     949     1308302 :     if ((keepfirst && kappa == 2) ||
     950      654097 :         cmprr(mulrr(gmael(r,kappa-1,kappa-1), delta), gel(s,kappa-1)) <= 0)
     951             :     { /* Step4: Success of Lovasz's condition */
     952      440182 :       alpha[kappa] = kappa;
     953      440182 :       tmp = mulrr(gmael(mu,kappa,kappa-1), gmael(r,kappa,kappa-1));
     954      440182 :       affrr(subrr(gel(s,kappa-1), tmp), gmael(r,kappa,kappa));
     955      440182 :       set_avma(av2); continue;
     956             :     }
     957             :     /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
     958      214023 :     if (DEBUGLEVEL>=4 && kappa==kappamax && signe(gel(s,kappa-1)))
     959           0 :       if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", expo(gel(s,1))); }
     960      214023 :     kappa2 = kappa;
     961             :     do {
     962      644194 :       kappa--;
     963      644194 :       if (kappa < zeros+2 + (keepfirst ? 1: 0)) break;
     964      614548 :       tmp = mulrr(gmael(r,kappa-1,kappa-1), delta);
     965      614548 :     } while (cmprr(gel(s,kappa-1), tmp) <= 0 );
     966      214023 :     set_avma(av2);
     967      214023 :     update_alpha(alpha, kappa, kappa2, kappamax);
     968             : 
     969             :     /* Step6: Update the mu's and r's */
     970      214023 :     rotate(mu,kappa2,kappa);
     971      214023 :     rotate(r,kappa2,kappa);
     972      214023 :     affrr(gel(s,kappa), gmael(r,kappa,kappa));
     973             : 
     974             :     /* Step7: Update B, appB, U, G */
     975      214023 :     rotate(B,kappa2,kappa);
     976      214023 :     rotate(appB,kappa2,kappa);
     977      214023 :     if (U) rotate(U,kappa2,kappa);
     978      214023 :     rotateG(G,kappa2,kappa, maxG, Gtmp);
     979             : 
     980             :     /* Step8: Prepare the next loop iteration */
     981      214023 :     if (kappa == zeros+1 && !signe(gmael(G,kappa,kappa)))
     982             :     {
     983           2 :       zeros++; kappa++;
     984           2 :       affrr(RgV_dotsquare(gel(appB,kappa)), gmael(G,kappa,kappa));
     985           2 :       affrr(gmael(G,kappa,kappa), gmael(r,kappa,kappa));
     986             :     }
     987             :   }
     988        2354 :   *pB=B; *pU=U; return zeros;
     989             : }
     990             : 
     991             : /************************* PROVED version (t_INT) ***********************/
     992             : /* dpe inspired by dpe.h by Patrick Pelissier, Paul Zimmermann
     993             :  * https://gforge.inria.fr/projects/dpe/
     994             :  */
     995             : 
     996             : typedef struct
     997             : {
     998             :   double d;  /* significand */
     999             :   long e; /* exponent */
    1000             : } dpe_t;
    1001             : 
    1002             : #define Dmael(x,i,j) (&((x)[i][j]))
    1003             : #define Del(x,i) (&((x)[i]))
    1004             : 
    1005             : static void
    1006     3040358 : dperotate(dpe_t **A, long k2, long k)
    1007             : {
    1008             :   long i;
    1009     3040358 :   dpe_t *B = A[k2];
    1010    12566062 :   for (i = k2; i > k; i--) A[i] = A[i-1];
    1011     3040358 :   A[k] = B;
    1012     3040358 : }
    1013             : 
    1014             : static void
    1015   833001992 : dpe_normalize0(dpe_t *x)
    1016             : {
    1017             :   int e;
    1018   833001992 :   x->d = frexp(x->d, &e);
    1019   833001992 :   x->e += e;
    1020   833001992 : }
    1021             : 
    1022             : static void
    1023   413237074 : dpe_normalize(dpe_t *x)
    1024             : {
    1025   413237074 :   if (x->d == 0.0)
    1026     2260279 :     x->e = -LONG_MAX;
    1027             :   else
    1028   410976795 :     dpe_normalize0(x);
    1029   413237142 : }
    1030             : 
    1031             : static GEN
    1032       93248 : dpetor(dpe_t *x)
    1033             : {
    1034       93248 :   GEN r = dbltor(x->d);
    1035       93248 :   if (signe(r)==0) return r;
    1036       93248 :   setexpo(r, x->e-1);
    1037       93248 :   return r;
    1038             : }
    1039             : 
    1040             : static void
    1041    93260932 : affdpe(dpe_t *y, dpe_t *x)
    1042             : {
    1043    93260932 :   x->d = y->d;
    1044    93260932 :   x->e = y->e;
    1045    93260932 : }
    1046             : 
    1047             : static void
    1048    65328654 : affidpe(GEN y, dpe_t *x)
    1049             : {
    1050    65328654 :   pari_sp av = avma;
    1051    65328654 :   GEN r = itor(y, DEFAULTPREC);
    1052    65327878 :   x->e = expo(r)+1;
    1053    65327878 :   setexpo(r,-1);
    1054    65327788 :   x->d = rtodbl(r);
    1055    65328110 :   set_avma(av);
    1056    65328178 : }
    1057             : 
    1058             : static void
    1059     6658256 : affdbldpe(double y, dpe_t *x)
    1060             : {
    1061     6658256 :   x->d = (double)y;
    1062     6658256 :   x->e = 0;
    1063     6658256 :   dpe_normalize(x);
    1064     6658258 : }
    1065             : 
    1066             : static void
    1067   395660729 : dpe_mulz(dpe_t *x, dpe_t *y, dpe_t *z)
    1068             : {
    1069   395660729 :   z->d = x->d * y->d;
    1070   395660729 :   if (z->d == 0.0)
    1071    22584710 :     z->e = -LONG_MAX;
    1072             :   else
    1073             :   {
    1074   373076019 :     z->e = x->e + y->e;
    1075   373076019 :     dpe_normalize0(z);
    1076             :   }
    1077   395661009 : }
    1078             : 
    1079             : static void
    1080    51310668 : dpe_divz(dpe_t *x, dpe_t *y, dpe_t *z)
    1081             : {
    1082    51310668 :   z->d = x->d / y->d;
    1083    51310668 :   if (z->d == 0.0)
    1084     2361067 :     z->e = -LONG_MAX;
    1085             :   else
    1086             :   {
    1087    48949601 :     z->e = x->e - y->e;
    1088    48949601 :     dpe_normalize0(z);
    1089             :   }
    1090    51310734 : }
    1091             : 
    1092             : static void
    1093      539807 : dpe_negz(dpe_t *y, dpe_t *x)
    1094             : {
    1095      539807 :   x->d = - y->d;
    1096      539807 :   x->e = y->e;
    1097      539807 : }
    1098             : 
    1099             : static void
    1100    30177587 : dpe_addz(dpe_t *y, dpe_t *z, dpe_t *x)
    1101             : {
    1102    30177587 :   if (y->e > z->e + 53)
    1103     1060602 :     affdpe(y, x);
    1104    29116985 :   else if (z->e > y->e + 53)
    1105      274619 :     affdpe(z, x);
    1106             :   else
    1107             :   {
    1108    28842366 :     long d = y->e - z->e;
    1109             : 
    1110    28842366 :     if (d >= 0)
    1111             :     {
    1112    21414458 :       x->d = y->d + ldexp(z->d, -d);
    1113    21414458 :       x->e  = y->e;
    1114             :     }
    1115             :     else
    1116             :     {
    1117     7427908 :       x->d = z->d + ldexp(y->d, d);
    1118     7427908 :       x->e = z->e;
    1119             :     }
    1120    28842366 :     dpe_normalize(x);
    1121             :   }
    1122    30177594 : }
    1123             : static void
    1124   411416015 : dpe_subz(dpe_t *y, dpe_t *z, dpe_t *x)
    1125             : {
    1126   411416015 :   if (y->e > z->e + 53)
    1127    39094972 :     affdpe(y, x);
    1128   372321043 :   else if (z->e > y->e + 53)
    1129      539807 :     dpe_negz(z, x);
    1130             :   else
    1131             :   {
    1132   371781236 :     long d = y->e - z->e;
    1133             : 
    1134   371781236 :     if (d >= 0)
    1135             :     {
    1136   339380038 :       x->d = y->d - ldexp(z->d, -d);
    1137   339380038 :       x->e = y->e;
    1138             :     }
    1139             :     else
    1140             :     {
    1141    32401198 :       x->d = ldexp(y->d, d) - z->d;
    1142    32401198 :       x->e = z->e;
    1143             :     }
    1144   371781236 :     dpe_normalize(x);
    1145             :   }
    1146   411416364 : }
    1147             : 
    1148             : static void
    1149     5954957 : dpe_muluz(dpe_t *y, ulong t, dpe_t *x)
    1150             : {
    1151     5954957 :   x->d = y->d * (double)t;
    1152     5954957 :   x->e = y->e;
    1153     5954957 :   dpe_normalize(x);
    1154     5954956 : }
    1155             : 
    1156             : static void
    1157     2317688 : dpe_addmuluz(dpe_t *y,  dpe_t *z, ulong t, dpe_t *x)
    1158             : {
    1159             :   dpe_t tmp;
    1160     2317688 :   dpe_muluz(z, t, &tmp);
    1161     2317688 :   dpe_addz(y, &tmp, x);
    1162     2317688 : }
    1163             : 
    1164             : static void
    1165     2446799 : dpe_submuluz(dpe_t *y,  dpe_t *z, ulong t, dpe_t *x)
    1166             : {
    1167             :   dpe_t tmp;
    1168     2446799 :   dpe_muluz(z, t, &tmp);
    1169     2446799 :   dpe_subz(y, &tmp, x);
    1170     2446799 : }
    1171             : 
    1172             : static void
    1173   380968643 : dpe_submulz(dpe_t *y,  dpe_t *z, dpe_t *t, dpe_t *x)
    1174             : {
    1175             :   dpe_t tmp;
    1176   380968643 :   dpe_mulz(z, t, &tmp);
    1177   380968664 :   dpe_subz(y, &tmp, x);
    1178   380968768 : }
    1179             : 
    1180             : static int
    1181    14692202 : dpe_cmp(dpe_t *x, dpe_t *y)
    1182             : {
    1183    14692202 :   int sx = x->d < 0. ? -1: x->d > 0.;
    1184    14692202 :   int sy = y->d < 0. ? -1: y->d > 0.;
    1185    14692202 :   int d  = sx - sy;
    1186             : 
    1187    14692202 :   if (d != 0)
    1188      142858 :     return d;
    1189    14549344 :   else if (x->e > y->e)
    1190     2340772 :     return (sx > 0) ? 1 : -1;
    1191    12208572 :   else if (y->e > x->e)
    1192     4814384 :     return (sx > 0) ? -1 : 1;
    1193             :   else
    1194     7394188 :     return (x->d < y->d) ? -1 : (x->d > y->d);
    1195             : }
    1196             : 
    1197             : static int
    1198    63395780 : dpe_abscmp(dpe_t *x, dpe_t *y)
    1199             : {
    1200    63395780 :   if (x->e > y->e)
    1201      450541 :     return 1;
    1202    62945239 :   else if (y->e > x->e)
    1203    59961275 :     return -1;
    1204             :   else
    1205     2983964 :     return (fabs(x->d) < fabs(y->d)) ? -1 : (fabs(x->d) > fabs(y->d));
    1206             : }
    1207             : 
    1208             : static int
    1209     9684531 : dpe_abssmall(dpe_t *x)
    1210             : {
    1211     9684531 :   return (x->e <= 0) || (x->e == 1 && fabs(x->d) <= .75);
    1212             : }
    1213             : 
    1214             : static int
    1215    14692197 : dpe_cmpmul(dpe_t *x, dpe_t *y, dpe_t *z)
    1216             : {
    1217             :   dpe_t t;
    1218    14692197 :   dpe_mulz(x,y,&t);
    1219    14692200 :   return dpe_cmp(&t, z);
    1220             : }
    1221             : 
    1222             : static dpe_t *
    1223    21840130 : cget_dpevec(long d)
    1224    21840130 : { return (dpe_t*) stack_malloc_align(d*sizeof(dpe_t), sizeof(dpe_t)); }
    1225             : 
    1226             : static dpe_t **
    1227     6658271 : cget_dpemat(long d) { return (dpe_t **) cgetg(d, t_VECSMALL); }
    1228             : 
    1229             : static GEN
    1230       19966 : dpeM_diagonal_shallow(dpe_t **m, long d)
    1231             : {
    1232             :   long i;
    1233       19966 :   GEN y = cgetg(d+1,t_VEC);
    1234      113214 :   for (i=1; i<=d; i++) gel(y, i) = dpetor(Dmael(m,i,i));
    1235       19966 :   return y;
    1236             : }
    1237             : 
    1238             : static void
    1239     9684509 : affii_or_copy_gc(pari_sp av, GEN x, GEN *y)
    1240             : {
    1241     9684509 :   long l = lg(*y);
    1242     9684509 :   if (lgefint(x) <= l && isonstack(*y))
    1243             :   {
    1244     9684081 :     affii(x,*y);
    1245     9684086 :     set_avma(av);
    1246             :   }
    1247             :   else
    1248         433 :     *y = gerepileuptoint(av, x);
    1249     9684520 : }
    1250             : 
    1251             : /* *x -= u*y */
    1252             : INLINE void
    1253    26431384 : submulziu(GEN *x, GEN y, ulong u)
    1254             : {
    1255             :   pari_sp av;
    1256    26431384 :   long ly = lgefint(y);
    1257    26431384 :   if (ly == 2) return;
    1258    19558641 :   av = avma;
    1259    19558641 :   (void)new_chunk(3+ly+lgefint(*x)); /* HACK */
    1260    19558718 :   y = mului(u,y);
    1261    19558674 :   set_avma(av); subzi(x, y);
    1262             : }
    1263             : 
    1264             : /* *x += u*y */
    1265             : INLINE void
    1266    24696100 : addmulziu(GEN *x, GEN y, ulong u)
    1267             : {
    1268             :   pari_sp av;
    1269    24696100 :   long ly = lgefint(y);
    1270    24696100 :   if (ly == 2) return;
    1271    18810341 :   av = avma;
    1272    18810341 :   (void)new_chunk(3+ly+lgefint(*x)); /* HACK */
    1273    18810386 :   y = mului(u,y);
    1274    18810351 :   set_avma(av); addzi(x, y);
    1275             : }
    1276             : 
    1277             : /************************** PROVED version (dpe) *************************/
    1278             : 
    1279             : /* Babai's Nearest Plane algorithm (iterative).
    1280             :  * Size-reduces b_kappa using mu_{i,j} and r_{i,j} for j<=i <kappa
    1281             :  * Update B[,kappa]; compute mu_{kappa,j}, r_{kappa,j} for j<=kappa and s[kappa]
    1282             :  * mu, r, s updated in place (affrr). Return 1 on failure, else 0. */
    1283             : static int
    1284    10286482 : Babai_dpe(pari_sp av, long kappa, GEN *pG, GEN *pB, GEN *pU, dpe_t **mu, dpe_t **r, dpe_t *s,
    1285             :       long a, long zeros, long maxG, dpe_t *eta)
    1286             : {
    1287    10286482 :   GEN G = *pG, B = *pB, U = *pU, ztmp;
    1288    10286482 :   long k, d, n, aa = a > zeros? a: zeros+1;
    1289    10286482 :   long emaxmu = EX0, emax2mu = EX0;
    1290             :   /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
    1291    10286482 :   d = U? lg(U)-1: 0;
    1292    10286482 :   n = B? nbrows(B): 0;
    1293     2170836 :   for (;;) {
    1294    12457354 :     int go_on = 0;
    1295    12457354 :     long i, j, emax3mu = emax2mu;
    1296             : 
    1297    12457354 :     if (gc_needed(av,2))
    1298             :     {
    1299           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
    1300           0 :       gc_lll(av,3,&G,&B,&U);
    1301             :     }
    1302             :     /* Step2: compute the GSO for stage kappa */
    1303    12457323 :     emax2mu = emaxmu; emaxmu = EX0;
    1304    63768035 :     for (j=aa; j<kappa; j++)
    1305             :     {
    1306             :       dpe_t g;
    1307    51310702 :       affidpe(gmael(G,kappa,j), &g);
    1308   375974441 :       for (k = zeros+1; k < j; k++)
    1309   324663793 :         dpe_submulz(&g, Dmael(mu,j,k), Dmael(r,kappa,k), &g);
    1310    51310648 :       affdpe(&g, Dmael(r,kappa,j));
    1311    51310667 :       dpe_divz(Dmael(r,kappa,j), Dmael(r,j,j), Dmael(mu,kappa,j));
    1312    51310671 :       emaxmu = maxss(emaxmu, Dmael(mu,kappa,j)->e);
    1313             :     }
    1314    12457333 :     if (emax3mu != EX0 && emax3mu <= emax2mu + 5) /* precision too low */
    1315           0 :     { *pG = G; *pB = B; *pU = U; return 1; }
    1316             : 
    1317    73682302 :     for (j=kappa-1; j>zeros; j--)
    1318    63395786 :       if (dpe_abscmp(Dmael(mu,kappa,j), eta) > 0) { go_on = 1; break; }
    1319             : 
    1320             :     /* Step3--5: compute the X_j's  */
    1321    12457333 :     if (go_on)
    1322    23470252 :       for (j=kappa-1; j>zeros; j--)
    1323             :       {
    1324             :         pari_sp btop;
    1325    21299420 :         dpe_t *tmp = Dmael(mu,kappa,j);
    1326    21299420 :         if (tmp->e < 0) continue; /* (essentially) size-reduced */
    1327             : 
    1328     9684530 :         if (gc_needed(av,2))
    1329             :         {
    1330           0 :           if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
    1331           0 :           gc_lll(av,3,&G,&B,&U);
    1332             :         }
    1333             :         /* we consider separately the case |X| = 1 */
    1334     9684530 :         if (dpe_abssmall(tmp))
    1335             :         {
    1336     8241037 :           if (tmp->d > 0) { /* in this case, X = 1 */
    1337    31513736 :             for (k=zeros+1; k<j; k++)
    1338    27386178 :               dpe_subz(Dmael(mu,kappa,k), Dmael(mu,j,k), Dmael(mu,kappa,k));
    1339    42750098 :             for (i=1; i<=n; i++)
    1340    38622540 :               subzi(&gmael(B,kappa,i), gmael(B,j,i));
    1341    68834844 :             for (i=1; i<=d; i++)
    1342    64707308 :               subzi(&gmael(U,kappa,i), gmael(U,j,i));
    1343     4127536 :             btop = avma;
    1344     4127536 :             ztmp = subii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
    1345     4127551 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    1346     4127553 :             affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
    1347    35983082 :             for (i=1; i<=j; i++)
    1348    31855526 :               subzi(&gmael(G,kappa,i), gmael(G,j,i));
    1349    34572353 :             for (i=j+1; i<kappa; i++)
    1350    30444799 :               subzi(&gmael(G,kappa,i), gmael(G,i,j));
    1351    27792725 :             for (i=kappa+1; i<=maxG; i++)
    1352    23665172 :               subzi(&gmael(G,i,kappa), gmael(G,i,j));
    1353             :           } else { /* otherwise X = -1 */
    1354    31397351 :             for (k=zeros+1; k<j; k++)
    1355    27283867 :               dpe_addz(Dmael(mu,kappa,k), Dmael(mu,j,k), Dmael(mu,kappa,k));
    1356    42696668 :             for (i=1; i<=n; i++)
    1357    38583184 :               addzi(&gmael(B,kappa,i),gmael(B,j,i));
    1358    68485756 :             for (i=1; i<=d; i++)
    1359    64372291 :               addzi(&gmael(U,kappa,i),gmael(U,j,i));
    1360     4113465 :             btop = avma;
    1361     4113465 :             ztmp = addii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
    1362     4113480 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    1363     4113479 :             affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
    1364    35787895 :             for (i=1; i<=j; i++)
    1365    31674413 :               addzi(&gmael(G,kappa,i), gmael(G,j,i));
    1366    34458589 :             for (i=j+1; i<kappa; i++)
    1367    30345106 :               addzi(&gmael(G,kappa,i), gmael(G,i,j));
    1368    27718097 :             for (i=kappa+1; i<=maxG; i++)
    1369    23604615 :               addzi(&gmael(G,i,kappa), gmael(G,i,j));
    1370             :           }
    1371     8241035 :           continue;
    1372             :         }
    1373             :         /* we have |X| >= 2 */
    1374     1443498 :         if (tmp->e < BITS_IN_LONG-1)
    1375             :         {
    1376     1268229 :           if (tmp->d > 0)
    1377             :           {
    1378      662302 :             ulong xx = (ulong) pari_rint(ldexp(tmp->d, tmp->e)); /* X fits in an ulong */
    1379     3109101 :             for (k=zeros+1; k<j; k++)
    1380     2446799 :               dpe_submuluz(Dmael(mu,kappa,k), Dmael(mu,j,k), xx, Dmael(mu,kappa,k));
    1381     4358711 :             for (i=1; i<=n; i++)
    1382     3696409 :               submulziu(&gmael(B,kappa,i), gmael(B,j,i), xx);
    1383    11195362 :             for (i=1; i<=d; i++)
    1384    10533060 :               submulziu(&gmael(U,kappa,i), gmael(U,j,i), xx);
    1385      662302 :             btop = avma;
    1386      662302 :             ztmp = submuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
    1387      662301 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    1388      662300 :             affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
    1389     4179472 :             for (i=1; i<=j; i++)
    1390     3517170 :               submulziu(&gmael(G,kappa,i), gmael(G,j,i), xx);
    1391     6602070 :             for (i=j+1; i<kappa; i++)
    1392     5939767 :               submulziu(&gmael(G,kappa,i), gmael(G,i,j), xx);
    1393     3407373 :             for (i=kappa+1; i<=maxG; i++)
    1394     2745069 :               submulziu(&gmael(G,i,kappa), gmael(G,i,j), xx);
    1395             :           }
    1396             :           else
    1397             :           {
    1398      605927 :             ulong xx = (ulong) pari_rint(ldexp(-tmp->d, tmp->e)); /* X fits in an ulong */
    1399     2923616 :             for (k=zeros+1; k<j; k++)
    1400     2317688 :               dpe_addmuluz(Dmael(mu,kappa,k), Dmael(mu,j,k), xx, Dmael(mu,kappa,k));
    1401     4297868 :             for (i=1; i<=n; i++)
    1402     3691940 :               addmulziu(&gmael(B,kappa,i), gmael(B,j,i), xx);
    1403    10260959 :             for (i=1; i<=d; i++)
    1404     9655030 :               addmulziu(&gmael(U,kappa,i), gmael(U,j,i), xx);
    1405      605929 :             btop = avma;
    1406      605929 :             ztmp = addmuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
    1407      605928 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    1408      605927 :             affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
    1409     3775519 :             for (i=1; i<=j; i++)
    1410     3169591 :               addmulziu(&gmael(G,kappa,i), gmael(G,j,i), xx);
    1411     6394449 :             for (i=j+1; i<kappa; i++)
    1412     5788520 :               addmulziu(&gmael(G,kappa,i), gmael(G,i,j), xx);
    1413     2997003 :             for (i=kappa+1; i<=maxG; i++)
    1414     2391074 :               addmulziu(&gmael(G,i,kappa), gmael(G,i,j), xx);
    1415             :           }
    1416             :         }
    1417             :         else
    1418             :         {
    1419      175269 :           long e = tmp->e - BITS_IN_LONG + 1;
    1420      175269 :           if (tmp->d > 0)
    1421             :           {
    1422       90612 :             ulong xx = (ulong) pari_rint(ldexp(tmp->d, BITS_IN_LONG - 1));
    1423      705052 :             for (k=zeros+1; k<j; k++)
    1424             :             {
    1425             :               dpe_t x;
    1426      614440 :               dpe_muluz(Dmael(mu,j,k), xx, &x);
    1427      614440 :               x.e += e;
    1428      614440 :               dpe_subz(Dmael(mu,kappa,k), &x, Dmael(mu,kappa,k));
    1429             :             }
    1430     1740992 :             for (i=1; i<=n; i++)
    1431     1650380 :               submulzu2n(&gmael(B,kappa,i), gmael(B,j,i), xx, e);
    1432      459717 :             for (i=1; i<=d; i++)
    1433      369105 :               submulzu2n(&gmael(U,kappa,i), gmael(U,j,i), xx, e);
    1434       90612 :             btop = avma;
    1435       90612 :             ztmp = submuliu2n(mulshift(gmael(G,j,j), sqru(xx), 2*e),
    1436       90612 :                 gmael(G,kappa,j), xx, e+1);
    1437       90612 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    1438       90612 :             affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
    1439      795857 :             for (i=1; i<=j; i++)
    1440      705245 :               submulzu2n(&gmael(G,kappa,i), gmael(G,j,i), xx, e);
    1441      757140 :             for (   ; i<kappa; i++)
    1442      666528 :               submulzu2n(&gmael(G,kappa,i), gmael(G,i,j), xx, e);
    1443      150589 :             for (i=kappa+1; i<=maxG; i++)
    1444       59977 :               submulzu2n(&gmael(G,i,kappa), gmael(G,i,j), xx, e);
    1445             :           } else
    1446             :           {
    1447       84657 :             ulong xx = (ulong) pari_rint(ldexp(-tmp->d, BITS_IN_LONG - 1));
    1448      660698 :             for (k=zeros+1; k<j; k++)
    1449             :             {
    1450             :               dpe_t x;
    1451      576033 :               dpe_muluz(Dmael(mu,j,k), xx, &x);
    1452      576033 :               x.e += e;
    1453      576033 :               dpe_addz(Dmael(mu,kappa,k), &x, Dmael(mu,kappa,k));
    1454             :             }
    1455     1714144 :             for (i=1; i<=n; i++)
    1456     1629479 :               addmulzu2n(&gmael(B,kappa,i), gmael(B,j,i), xx, e);
    1457      373223 :             for (i=1; i<=d; i++)
    1458      288558 :               addmulzu2n(&gmael(U,kappa,i), gmael(U,j,i), xx, e);
    1459       84665 :             btop = avma;
    1460       84665 :             ztmp = addmuliu2n(mulshift(gmael(G,j,j), sqru(xx), 2*e),
    1461       84665 :                 gmael(G,kappa,j), xx, e+1);
    1462       84665 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    1463       84665 :             affii_or_copy_gc(btop, ztmp, &gmael(G,kappa,kappa));
    1464      745524 :             for (i=1; i<=j; i++)
    1465      660859 :               addmulzu2n(&gmael(G,kappa,i), gmael(G,j,i), xx, e);
    1466      724460 :             for (   ; i<kappa; i++)
    1467      639795 :               addmulzu2n(&gmael(G,kappa,i), gmael(G,i,j), xx, e);
    1468      124848 :             for (i=kappa+1; i<=maxG; i++)
    1469       40183 :               addmulzu2n(&gmael(G,i,kappa), gmael(G,i,j), xx, e);
    1470             :           }
    1471             :         }
    1472             :       }
    1473    12457348 :     if (!go_on) break; /* Anything happened? */
    1474     2170836 :     aa = zeros+1;
    1475             :   }
    1476             : 
    1477    10286512 :   affidpe(gmael(G,kappa,kappa), Del(s,zeros+1));
    1478             :   /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
    1479    57825199 :   for (k=zeros+1; k<=kappa-2; k++)
    1480    47538701 :     dpe_submulz(Del(s,k), Dmael(mu,kappa,k), Dmael(r,kappa,k), Del(s,k+1));
    1481    10286498 :   *pG = G; *pB = B; *pU = U; return 0;
    1482             : }
    1483             : 
    1484             : /* G integral Gram matrix, LLL-reduces (G,B,U) in place [apply base change
    1485             :  * transforms to B and U]. If (keepfirst), never swap with first vector.
    1486             :  * If G = NULL, we compute the Gram matrix incrementally.
    1487             :  * Return -1 on failure, else zeros = dim Kernel (>= 0) */
    1488             : static long
    1489     3329129 : fplll_dpe(GEN *pG, GEN *pB, GEN *pU, GEN *pr, double DELTA, double ETA,
    1490             :       long keepfirst)
    1491             : {
    1492             :   pari_sp av;
    1493     3329129 :   GEN Gtmp, alpha, G = *pG, B = *pB, U = *pU;
    1494     3329129 :   long d, maxG, kappa, kappa2, i, j, zeros, kappamax, incgram = !G, cnt = 0;
    1495             :   dpe_t delta, eta, **mu, **r, *s;
    1496     3329129 :   affdbldpe(DELTA,&delta);
    1497     3329131 :   affdbldpe(ETA,&eta);
    1498             : 
    1499     3329150 :   if (incgram)
    1500             :   { /* incremental Gram matrix */
    1501     3246592 :     maxG = 2; d = lg(B)-1;
    1502     3246592 :     G = zeromatcopy(d, d);
    1503             :   }
    1504             :   else
    1505       82558 :     maxG = d = lg(G)-1;
    1506             : 
    1507     3329149 :   mu = cget_dpemat(d+1);
    1508     3329133 :   r  = cget_dpemat(d+1);
    1509     3329137 :   s  = cget_dpevec(d+1);
    1510    12584702 :   for (j = 1; j <= d; j++)
    1511             :   {
    1512     9255553 :     mu[j]= cget_dpevec(d+1);
    1513     9255553 :     r[j] = cget_dpevec(d+1);
    1514             :   }
    1515     3329149 :   Gtmp = cgetg(d+1, t_VEC);
    1516     3329142 :   alpha = cgetg(d+1, t_VECSMALL);
    1517     3329141 :   av = avma;
    1518             : 
    1519             :   /* Step2: Initializing the main loop */
    1520     3329141 :   kappamax = 1;
    1521     3329141 :   i = 1;
    1522             :   do {
    1523     3696637 :     if (incgram) gmael(G,i,i) = ZV_dotsquare(gel(B,i));
    1524     3696552 :     affidpe(gmael(G,i,i), Dmael(r,i,i));
    1525     3696571 :   } while (!signe(gmael(G,i,i)) && ++i <= d);
    1526     3329075 :   zeros = i-1; /* all basis vectors b_i with i <= zeros are zero vectors */
    1527     3329075 :   kappa = i;
    1528    12216960 :   for (i=zeros+1; i<=d; i++) alpha[i]=1;
    1529             : 
    1530    13615619 :   while (++kappa <= d)
    1531             :   {
    1532    10286501 :     if (kappa > kappamax)
    1533             :     {
    1534     5558823 :       if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
    1535     5558823 :       kappamax = kappa;
    1536     5558823 :       if (incgram)
    1537             :       {
    1538    21147530 :         for (i=zeros+1; i<=kappa; i++)
    1539    15819985 :           gmael(G,kappa,i) = ZV_dotproduct(gel(B,kappa), gel(B,i));
    1540     5327545 :         maxG = kappamax;
    1541             :       }
    1542             :     }
    1543             :     /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
    1544    10286240 :     if (Babai_dpe(av, kappa, &G,&B,&U, mu,r,s, alpha[kappa], zeros, maxG, &eta))
    1545           0 :     { *pG = incgram? NULL: G; *pB = B; *pU = U; return -1; }
    1546    20475109 :     if ((keepfirst && kappa == 2) ||
    1547    10188592 :         dpe_cmpmul(Dmael(r,kappa-1,kappa-1), &delta, Del(s,kappa-1)) <= 0)
    1548             :     { /* Step4: Success of Lovasz's condition */
    1549     8766339 :       alpha[kappa] = kappa;
    1550     8766339 :       dpe_submulz(Del(s,kappa-1), Dmael(mu,kappa,kappa-1), Dmael(r,kappa,kappa-1), Dmael(r,kappa,kappa));
    1551     8766362 :       continue;
    1552             :     }
    1553             :     /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
    1554     1520178 :     if (DEBUGLEVEL>=4 && kappa==kappamax && Del(s,kappa-1)->d)
    1555           0 :       if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", Del(s,1)->e-1); }
    1556     1520178 :     kappa2 = kappa;
    1557             :     do {
    1558     4762851 :       kappa--;
    1559     4762851 :       if (kappa < zeros+2 + (keepfirst ? 1: 0)) break;
    1560     4503597 :     } while (dpe_cmpmul(Dmael(r,kappa-1,kappa-1), &delta, Del(s,kappa-1)) >= 0);
    1561     1520178 :     update_alpha(alpha, kappa, kappa2, kappamax);
    1562             : 
    1563             :     /* Step6: Update the mu's and r's */
    1564     1520179 :     dperotate(mu, kappa2, kappa);
    1565     1520179 :     dperotate(r, kappa2, kappa);
    1566     1520179 :     affdpe(Del(s,kappa), Dmael(r,kappa,kappa));
    1567             : 
    1568             :     /* Step7: Update G, B, U */
    1569     1520179 :     if (U) rotate(U, kappa2, kappa);
    1570     1520179 :     if (B) rotate(B, kappa2, kappa);
    1571     1520179 :     rotateG(G,kappa2,kappa, maxG, Gtmp);
    1572             : 
    1573             :     /* Step8: Prepare the next loop iteration */
    1574     1520179 :     if (kappa == zeros+1 && !signe(gmael(G,kappa,kappa)))
    1575             :     {
    1576       35176 :       zeros++; kappa++;
    1577       35176 :       affidpe(gmael(G,kappa,kappa), Dmael(r,kappa,kappa));
    1578             :     }
    1579             :   }
    1580     3329118 :   if (pr) *pr = dpeM_diagonal_shallow(r,d);
    1581     3329118 :   *pG = G; *pB = B; *pU = U; return zeros; /* success */
    1582             : }
    1583             : 
    1584             : 
    1585             : /************************** PROVED version (t_INT) *************************/
    1586             : 
    1587             : /* Babai's Nearest Plane algorithm (iterative).
    1588             :  * Size-reduces b_kappa using mu_{i,j} and r_{i,j} for j<=i <kappa
    1589             :  * Update B[,kappa]; compute mu_{kappa,j}, r_{kappa,j} for j<=kappa and s[kappa]
    1590             :  * mu, r, s updated in place (affrr). Return 1 on failure, else 0. */
    1591             : static int
    1592           0 : Babai(pari_sp av, long kappa, GEN *pG, GEN *pB, GEN *pU, GEN mu, GEN r, GEN s,
    1593             :       long a, long zeros, long maxG, GEN eta, long prec)
    1594             : {
    1595           0 :   GEN G = *pG, B = *pB, U = *pU, ztmp;
    1596           0 :   long k, aa = a > zeros? a: zeros+1;
    1597           0 :   const long n = B? nbrows(B): 0, d = U ? lg(U)-1: 0, bit = prec2nbits(prec);
    1598           0 :   long emaxmu = EX0, emax2mu = EX0;
    1599             :   /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
    1600             : 
    1601           0 :   for (;;) {
    1602           0 :     int go_on = 0;
    1603           0 :     long i, j, emax3mu = emax2mu;
    1604             : 
    1605           0 :     if (gc_needed(av,2))
    1606             :     {
    1607           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
    1608           0 :       gc_lll(av,3,&G,&B,&U);
    1609             :     }
    1610             :     /* Step2: compute the GSO for stage kappa */
    1611           0 :     emax2mu = emaxmu; emaxmu = EX0;
    1612           0 :     for (j=aa; j<kappa; j++)
    1613             :     {
    1614           0 :       pari_sp btop = avma;
    1615           0 :       GEN g = gmael(G,kappa,j);
    1616           0 :       for (k = zeros+1; k < j; k++)
    1617           0 :         g = mpsub(g, mulrr(gmael(mu,j,k), gmael(r,kappa,k)));
    1618           0 :       mpaff(g, gmael(r,kappa,j));
    1619           0 :       affrr(divrr(gmael(r,kappa,j), gmael(r,j,j)), gmael(mu,kappa,j));
    1620           0 :       emaxmu = maxss(emaxmu, expo(gmael(mu,kappa,j)));
    1621           0 :       set_avma(btop);
    1622             :     }
    1623           0 :     if (emax3mu != EX0 && emax3mu <= emax2mu + 5) /* precision too low */
    1624           0 :     { *pG = G; *pB = B; *pU = U; return 1; }
    1625             : 
    1626           0 :     for (j=kappa-1; j>zeros; j--)
    1627           0 :       if (abscmprr(gmael(mu,kappa,j), eta) > 0) { go_on = 1; break; }
    1628             : 
    1629             :     /* Step3--5: compute the X_j's  */
    1630           0 :     if (go_on)
    1631           0 :       for (j=kappa-1; j>zeros; j--)
    1632             :       {
    1633             :         pari_sp btop;
    1634           0 :         GEN tmp = gmael(mu,kappa,j);
    1635           0 :         if (absrsmall(tmp)) continue; /* size-reduced */
    1636             : 
    1637           0 :         if (gc_needed(av,2))
    1638             :         {
    1639           0 :           if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
    1640           0 :           gc_lll(av,3,&G,&B,&U);
    1641             :         }
    1642           0 :         btop = avma;
    1643             :         /* we consider separately the case |X| = 1 */
    1644           0 :         if (absrsmall2(tmp))
    1645             :         {
    1646           0 :           if (signe(tmp) > 0) { /* in this case, X = 1 */
    1647           0 :             for (k=zeros+1; k<j; k++)
    1648           0 :               affrr(subrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
    1649           0 :             set_avma(btop);
    1650           0 :             for (i=1; i<=n; i++)
    1651           0 :               gmael(B,kappa,i) = subii(gmael(B,kappa,i), gmael(B,j,i));
    1652           0 :             for (i=1; i<=d; i++)
    1653           0 :               gmael(U,kappa,i) = subii(gmael(U,kappa,i), gmael(U,j,i));
    1654           0 :             btop = avma;
    1655           0 :             ztmp = subii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
    1656           0 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    1657           0 :             gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
    1658           0 :             for (i=1; i<=j; i++)
    1659           0 :               gmael(G,kappa,i) = subii(gmael(G,kappa,i), gmael(G,j,i));
    1660           0 :             for (i=j+1; i<kappa; i++)
    1661           0 :               gmael(G,kappa,i) = subii(gmael(G,kappa,i), gmael(G,i,j));
    1662           0 :             for (i=kappa+1; i<=maxG; i++)
    1663           0 :               gmael(G,i,kappa) = subii(gmael(G,i,kappa), gmael(G,i,j));
    1664             :           } else { /* otherwise X = -1 */
    1665           0 :             for (k=zeros+1; k<j; k++)
    1666           0 :               affrr(addrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
    1667           0 :             set_avma(btop);
    1668           0 :             for (i=1; i<=n; i++)
    1669           0 :               gmael(B,kappa,i) = addii(gmael(B,kappa,i),gmael(B,j,i));
    1670           0 :             for (i=1; i<=d; i++)
    1671           0 :               gmael(U,kappa,i) = addii(gmael(U,kappa,i),gmael(U,j,i));
    1672           0 :             btop = avma;
    1673           0 :             ztmp = addii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
    1674           0 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    1675           0 :             gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
    1676           0 :             for (i=1; i<=j; i++)
    1677           0 :               gmael(G,kappa,i) = addii(gmael(G,kappa,i), gmael(G,j,i));
    1678           0 :             for (i=j+1; i<kappa; i++)
    1679           0 :               gmael(G,kappa,i) = addii(gmael(G,kappa,i), gmael(G,i,j));
    1680           0 :             for (i=kappa+1; i<=maxG; i++)
    1681           0 :               gmael(G,i,kappa) = addii(gmael(G,i,kappa), gmael(G,i,j));
    1682             :           }
    1683           0 :           continue;
    1684             :         }
    1685             :         /* we have |X| >= 2 */
    1686           0 :         if (expo(tmp) < BITS_IN_LONG)
    1687             :         {
    1688           0 :           ulong xx = roundr_safe(tmp)[2]; /* X fits in an ulong */
    1689           0 :           if (signe(tmp) > 0) /* = xx */
    1690             :           {
    1691           0 :             for (k=zeros+1; k<j; k++)
    1692           0 :               affrr(subrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
    1693           0 :                   gmael(mu,kappa,k));
    1694           0 :             set_avma(btop);
    1695           0 :             for (i=1; i<=n; i++)
    1696           0 :               gmael(B,kappa,i) = submuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
    1697           0 :             for (i=1; i<=d; i++)
    1698           0 :               gmael(U,kappa,i) = submuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
    1699           0 :             btop = avma;
    1700           0 :             ztmp = submuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
    1701           0 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    1702           0 :             gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
    1703           0 :             for (i=1; i<=j; i++)
    1704           0 :               gmael(G,kappa,i) = submuliu_inplace(gmael(G,kappa,i), gmael(G,j,i), xx);
    1705           0 :             for (i=j+1; i<kappa; i++)
    1706           0 :               gmael(G,kappa,i) = submuliu_inplace(gmael(G,kappa,i), gmael(G,i,j), xx);
    1707           0 :             for (i=kappa+1; i<=maxG; i++)
    1708           0 :               gmael(G,i,kappa) = submuliu_inplace(gmael(G,i,kappa), gmael(G,i,j), xx);
    1709             :           }
    1710             :           else /* = -xx */
    1711             :           {
    1712           0 :             for (k=zeros+1; k<j; k++)
    1713           0 :               affrr(addrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k))),
    1714           0 :                   gmael(mu,kappa,k));
    1715           0 :             set_avma(btop);
    1716           0 :             for (i=1; i<=n; i++)
    1717           0 :               gmael(B,kappa,i) = addmuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
    1718           0 :             for (i=1; i<=d; i++)
    1719           0 :               gmael(U,kappa,i) = addmuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
    1720           0 :             btop = avma;
    1721           0 :             ztmp = addmuliu2n(mulii(gmael(G,j,j), sqru(xx)), gmael(G,kappa,j), xx, 1);
    1722           0 :             ztmp = addii(gmael(G,kappa,kappa), ztmp);
    1723           0 :             gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
    1724           0 :             for (i=1; i<=j; i++)
    1725           0 :               gmael(G,kappa,i) = addmuliu_inplace(gmael(G,kappa,i), gmael(G,j,i), xx);
    1726           0 :             for (i=j+1; i<kappa; i++)
    1727           0 :               gmael(G,kappa,i) = addmuliu_inplace(gmael(G,kappa,i), gmael(G,i,j), xx);
    1728           0 :             for (i=kappa+1; i<=maxG; i++)
    1729           0 :               gmael(G,i,kappa) = addmuliu_inplace(gmael(G,i,kappa), gmael(G,i,j), xx);
    1730             :           }
    1731             :         }
    1732             :         else
    1733             :         {
    1734             :           long e;
    1735           0 :           GEN X = truncexpo(tmp, bit, &e); /* tmp ~ X * 2^e */
    1736           0 :           btop = avma;
    1737           0 :           for (k=zeros+1; k<j; k++)
    1738             :           {
    1739           0 :             GEN x = mulir(X, gmael(mu,j,k));
    1740           0 :             if (e) shiftr_inplace(x, e);
    1741           0 :             affrr(subrr(gmael(mu,kappa,k), x), gmael(mu,kappa,k));
    1742             :           }
    1743           0 :           set_avma(btop);
    1744           0 :           for (i=1; i<=n; i++)
    1745           0 :             gmael(B,kappa,i) = submulshift(gmael(B,kappa,i), gmael(B,j,i), X, e);
    1746           0 :           for (i=1; i<=d; i++)
    1747           0 :             gmael(U,kappa,i) = submulshift(gmael(U,kappa,i), gmael(U,j,i), X, e);
    1748           0 :           btop = avma;
    1749           0 :           ztmp = submulshift(mulshift(gmael(G,j,j), sqri(X), 2*e),
    1750           0 :               gmael(G,kappa,j), X, e+1);
    1751           0 :           ztmp = addii(gmael(G,kappa,kappa), ztmp);
    1752           0 :           gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
    1753           0 :           for (i=1; i<=j; i++)
    1754           0 :             gmael(G,kappa,i) = submulshift(gmael(G,kappa,i), gmael(G,j,i), X, e);
    1755           0 :           for (   ; i<kappa; i++)
    1756           0 :             gmael(G,kappa,i) = submulshift(gmael(G,kappa,i), gmael(G,i,j), X, e);
    1757           0 :           for (i=kappa+1; i<=maxG; i++)
    1758           0 :             gmael(G,i,kappa) = submulshift(gmael(G,i,kappa), gmael(G,i,j), X, e);
    1759             :         }
    1760             :       }
    1761           0 :     if (!go_on) break; /* Anything happened? */
    1762           0 :     aa = zeros+1;
    1763             :   }
    1764             : 
    1765           0 :   affir(gmael(G,kappa,kappa), gel(s,zeros+1));
    1766             :   /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
    1767           0 :   av = avma;
    1768           0 :   for (k=zeros+1; k<=kappa-2; k++)
    1769           0 :     affrr(subrr(gel(s,k), mulrr(gmael(mu,kappa,k), gmael(r,kappa,k))),
    1770           0 :           gel(s,k+1));
    1771           0 :   *pG = G; *pB = B; *pU = U; return gc_bool(av, 0);
    1772             : }
    1773             : 
    1774             : /* G integral Gram matrix, LLL-reduces (G,B,U) in place [apply base change
    1775             :  * transforms to B and U]. If (keepfirst), never swap with first vector.
    1776             :  * If G = NULL, we compute the Gram matrix incrementally.
    1777             :  * Return -1 on failure, else zeros = dim Kernel (>= 0) */
    1778             : static long
    1779           0 : fplll(GEN *pG, GEN *pB, GEN *pU, GEN *pr, double DELTA, double ETA,
    1780             :       long keepfirst, long prec)
    1781             : {
    1782             :   pari_sp av, av2;
    1783           0 :   GEN mu, r, s, tmp, Gtmp, alpha, G = *pG, B = *pB, U = *pU;
    1784           0 :   GEN delta = dbltor(DELTA), eta = dbltor(ETA);
    1785           0 :   long d, maxG, kappa, kappa2, i, j, zeros, kappamax, incgram = !G, cnt = 0;
    1786             : 
    1787           0 :   if (incgram)
    1788             :   { /* incremental Gram matrix */
    1789           0 :     maxG = 2; d = lg(B)-1;
    1790           0 :     G = zeromatcopy(d, d);
    1791             :   }
    1792             :   else
    1793           0 :     maxG = d = lg(G)-1;
    1794             : 
    1795           0 :   mu = cgetg(d+1, t_MAT);
    1796           0 :   r  = cgetg(d+1, t_MAT);
    1797           0 :   s  = cgetg(d+1, t_VEC);
    1798           0 :   for (j = 1; j <= d; j++)
    1799             :   {
    1800           0 :     GEN M = cgetg(d+1, t_COL), R = cgetg(d+1, t_COL);
    1801           0 :     gel(mu,j)= M;
    1802           0 :     gel(r,j) = R;
    1803           0 :     gel(s,j) = cgetr(prec);
    1804           0 :     for (i = 1; i <= d; i++)
    1805             :     {
    1806           0 :       gel(R,i) = cgetr(prec);
    1807           0 :       gel(M,i) = cgetr(prec);
    1808             :     }
    1809             :   }
    1810           0 :   Gtmp = cgetg(d+1, t_VEC);
    1811           0 :   alpha = cgetg(d+1, t_VECSMALL);
    1812           0 :   av = avma;
    1813             : 
    1814             :   /* Step2: Initializing the main loop */
    1815           0 :   kappamax = 1;
    1816           0 :   i = 1;
    1817             :   do {
    1818           0 :     if (incgram) gmael(G,i,i) = ZV_dotsquare(gel(B,i));
    1819           0 :     affir(gmael(G,i,i), gmael(r,i,i));
    1820           0 :   } while (!signe(gmael(G,i,i)) && ++i <= d);
    1821           0 :   zeros = i-1; /* all basis vectors b_i with i <= zeros are zero vectors */
    1822           0 :   kappa = i;
    1823           0 :   for (i=zeros+1; i<=d; i++) alpha[i]=1;
    1824             : 
    1825           0 :   while (++kappa <= d)
    1826             :   {
    1827           0 :     if (kappa > kappamax)
    1828             :     {
    1829           0 :       if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
    1830           0 :       kappamax = kappa;
    1831           0 :       if (incgram)
    1832             :       {
    1833           0 :         for (i=zeros+1; i<=kappa; i++)
    1834           0 :           gmael(G,kappa,i) = ZV_dotproduct(gel(B,kappa), gel(B,i));
    1835           0 :         maxG = kappamax;
    1836             :       }
    1837             :     }
    1838             :     /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
    1839           0 :     if (Babai(av, kappa, &G,&B,&U, mu,r,s, alpha[kappa], zeros, maxG, eta, prec))
    1840           0 :     { *pG = incgram? NULL: G; *pB = B; *pU = U; return -1; }
    1841           0 :     av2 = avma;
    1842           0 :     if ((keepfirst && kappa == 2) ||
    1843           0 :         cmprr(mulrr(gmael(r,kappa-1,kappa-1), delta), gel(s,kappa-1)) <= 0)
    1844             :     { /* Step4: Success of Lovasz's condition */
    1845           0 :       alpha[kappa] = kappa;
    1846           0 :       tmp = mulrr(gmael(mu,kappa,kappa-1), gmael(r,kappa,kappa-1));
    1847           0 :       affrr(subrr(gel(s,kappa-1), tmp), gmael(r,kappa,kappa));
    1848           0 :       set_avma(av2); continue;
    1849             :     }
    1850             :     /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
    1851           0 :     if (DEBUGLEVEL>=4 && kappa==kappamax && signe(gel(s,kappa-1)))
    1852           0 :       if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", expo(gel(s,1))); }
    1853           0 :     kappa2 = kappa;
    1854             :     do {
    1855           0 :       kappa--;
    1856           0 :       if (kappa < zeros+2 + (keepfirst ? 1: 0)) break;
    1857           0 :       tmp = mulrr(gmael(r,kappa-1,kappa-1), delta);
    1858           0 :     } while (cmprr(gel(s,kappa-1), tmp) <= 0);
    1859           0 :     set_avma(av2);
    1860           0 :     update_alpha(alpha, kappa, kappa2, kappamax);
    1861             : 
    1862             :     /* Step6: Update the mu's and r's */
    1863           0 :     rotate(mu, kappa2, kappa);
    1864           0 :     rotate(r, kappa2, kappa);
    1865           0 :     affrr(gel(s,kappa), gmael(r,kappa,kappa));
    1866             : 
    1867             :     /* Step7: Update G, B, U */
    1868           0 :     if (U) rotate(U, kappa2, kappa);
    1869           0 :     if (B) rotate(B, kappa2, kappa);
    1870           0 :     rotateG(G,kappa2,kappa, maxG, Gtmp);
    1871             : 
    1872             :     /* Step8: Prepare the next loop iteration */
    1873           0 :     if (kappa == zeros+1 && !signe(gmael(G,kappa,kappa)))
    1874             :     {
    1875           0 :       zeros++; kappa++;
    1876           0 :       affir(gmael(G,kappa,kappa), gmael(r,kappa,kappa));
    1877             :     }
    1878             :   }
    1879           0 :   if (pr) *pr = RgM_diagonal_shallow(r);
    1880           0 :   *pG = G; *pB = B; *pU = U; return zeros; /* success */
    1881             : }
    1882             : 
    1883             : /* Assume x a ZM, if pN != NULL, set it to Gram-Schmidt (squared) norms
    1884             :  * The following modes are supported:
    1885             :  * - flag & LLL_INPLACE: x a lattice basis, return x*U
    1886             :  * - flag & LLL_GRAM: x a Gram matrix / else x a lattice basis; return
    1887             :  *     LLL base change matrix U [LLL_IM]
    1888             :  *     kernel basis [LLL_KER, nonreduced]
    1889             :  *     both [LLL_ALL] */
    1890             : GEN
    1891     3452156 : ZM_lll_norms(GEN x, double DELTA, long flag, GEN *pN)
    1892             : {
    1893     3452156 :   pari_sp av = avma;
    1894     3452156 :   const double ETA = 0.51;
    1895     3452156 :   const long keepfirst = flag & LLL_KEEP_FIRST;
    1896     3452156 :   long p, zeros = -1, n = lg(x)-1;
    1897             :   GEN G, B, U;
    1898             :   pari_timer T;
    1899     3452156 :   if (n <= 1) return lll_trivial(x, flag);
    1900     3344251 :   if (nbrows(x)==0)
    1901             :   {
    1902       15190 :     if (flag & LLL_KER) return matid(n);
    1903       15190 :     if (flag & (LLL_INPLACE|LLL_IM)) return cgetg(1,t_MAT);
    1904           0 :     retmkvec2(matid(n), cgetg(1,t_MAT));
    1905             :   }
    1906     3329104 :   x = gcopy(x);
    1907     3329133 :   if (flag & LLL_GRAM)
    1908       82548 :   { G = x; B = NULL; U = matid(n); }
    1909             :   else
    1910     3246585 :   { G = NULL; B = x; U = (flag & LLL_INPLACE)? NULL: matid(n); }
    1911     3329137 :   if(DEBUGLEVEL>=4) timer_start(&T);
    1912     3329137 :   if (!(flag & LLL_GRAM))
    1913             :   {
    1914             :     long t;
    1915     3246576 :     if(DEBUGLEVEL>=4)
    1916           0 :       err_printf("Entering L^2 (double): LLL-parameters (%.3f,%.3f)\n",
    1917             :                  DELTA,ETA);
    1918     3246576 :     zeros = fplll_fast(&B, &U, DELTA, ETA, keepfirst);
    1919     3246576 :     if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
    1920     3250106 :     for (p = DEFAULTPREC, t = 0; zeros < 0 && t < 1 ; p += EXTRAPREC64, t++)
    1921             :     {
    1922        3531 :       if (DEBUGLEVEL>=4)
    1923           0 :         err_printf("Entering L^2 (heuristic): LLL-parameters (%.3f,%.3f), prec = %d/%d\n", DELTA, ETA, p, p);
    1924        3531 :       zeros = fplll_heuristic(&B, &U, DELTA, ETA, keepfirst, p, p);
    1925        3531 :       gc_lll(av, 2, &B, &U);
    1926        3531 :       if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
    1927             :     }
    1928             :   }
    1929     3329136 :   if(DEBUGLEVEL>=4)
    1930           0 :     err_printf("Entering L^2 (dpe): LLL-parameters (%.3f,%.3f)\n", DELTA,ETA);
    1931     3329136 :   zeros = fplll_dpe(&G, &B, &U, pN, DELTA, ETA, keepfirst);
    1932     3329115 :   if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
    1933     3329128 :   if (zeros < 0)
    1934           0 :     for (p = DEFAULTPREC;; p += EXTRAPREC64)
    1935             :     {
    1936           0 :       if (DEBUGLEVEL>=4)
    1937           0 :         err_printf("Entering L^2: LLL-parameters (%.3f,%.3f), prec = %d\n",
    1938             :                     DELTA,ETA, p);
    1939           0 :       zeros = fplll(&G, &B, &U, pN, DELTA, ETA, keepfirst, p);
    1940           0 :       if (DEBUGLEVEL>=4) timer_printf(&T, zeros < 0? "LLL (failed)": "LLL");
    1941           0 :       if (zeros >= 0) break;
    1942           0 :       gc_lll(av, 3, &G, &B, &U);
    1943             :     }
    1944     3329128 :   return lll_finish(U? U: B, zeros, flag);
    1945             : }
    1946             : 
    1947             : /********************************************************************/
    1948             : /**                                                                **/
    1949             : /**                        LLL OVER K[X]                           **/
    1950             : /**                                                                **/
    1951             : /********************************************************************/
    1952             : static int
    1953         504 : pslg(GEN x)
    1954             : {
    1955             :   long tx;
    1956         504 :   if (gequal0(x)) return 2;
    1957         448 :   tx = typ(x); return is_scalar_t(tx)? 3: lg(x);
    1958             : }
    1959             : 
    1960             : static int
    1961         196 : REDgen(long k, long l, GEN h, GEN L, GEN B)
    1962             : {
    1963         196 :   GEN q, u = gcoeff(L,k,l);
    1964             :   long i;
    1965             : 
    1966         196 :   if (pslg(u) < pslg(B)) return 0;
    1967             : 
    1968         140 :   q = gneg(gdeuc(u,B));
    1969         140 :   gel(h,k) = gadd(gel(h,k), gmul(q,gel(h,l)));
    1970         140 :   for (i=1; i<l; i++) gcoeff(L,k,i) = gadd(gcoeff(L,k,i), gmul(q,gcoeff(L,l,i)));
    1971         140 :   gcoeff(L,k,l) = gadd(gcoeff(L,k,l), gmul(q,B)); return 1;
    1972             : }
    1973             : 
    1974             : static int
    1975         196 : do_SWAPgen(GEN h, GEN L, GEN B, long k, GEN fl, int *flc)
    1976             : {
    1977             :   GEN p1, la, la2, Bk;
    1978             :   long ps1, ps2, i, j, lx;
    1979             : 
    1980         196 :   if (!fl[k-1]) return 0;
    1981             : 
    1982         140 :   la = gcoeff(L,k,k-1); la2 = gsqr(la);
    1983         140 :   Bk = gel(B,k);
    1984         140 :   if (fl[k])
    1985             :   {
    1986          56 :     GEN q = gadd(la2, gmul(gel(B,k-1),gel(B,k+1)));
    1987          56 :     ps1 = pslg(gsqr(Bk));
    1988          56 :     ps2 = pslg(q);
    1989          56 :     if (ps1 <= ps2 && (ps1 < ps2 || !*flc)) return 0;
    1990          28 :     *flc = (ps1 != ps2);
    1991          28 :     gel(B,k) = gdiv(q, Bk);
    1992             :   }
    1993             : 
    1994         112 :   swap(gel(h,k-1), gel(h,k)); lx = lg(L);
    1995         112 :   for (j=1; j<k-1; j++) swap(gcoeff(L,k-1,j), gcoeff(L,k,j));
    1996         112 :   if (fl[k])
    1997             :   {
    1998          28 :     for (i=k+1; i<lx; i++)
    1999             :     {
    2000           0 :       GEN t = gcoeff(L,i,k);
    2001           0 :       p1 = gsub(gmul(gel(B,k+1),gcoeff(L,i,k-1)), gmul(la,t));
    2002           0 :       gcoeff(L,i,k) = gdiv(p1, Bk);
    2003           0 :       p1 = gadd(gmul(la,gcoeff(L,i,k-1)), gmul(gel(B,k-1),t));
    2004           0 :       gcoeff(L,i,k-1) = gdiv(p1, Bk);
    2005             :     }
    2006             :   }
    2007          84 :   else if (!gequal0(la))
    2008             :   {
    2009          28 :     p1 = gdiv(la2, Bk);
    2010          28 :     gel(B,k+1) = gel(B,k) = p1;
    2011          28 :     for (i=k+2; i<=lx; i++) gel(B,i) = gdiv(gmul(p1,gel(B,i)),Bk);
    2012          28 :     for (i=k+1; i<lx; i++)
    2013           0 :       gcoeff(L,i,k-1) = gdiv(gmul(la,gcoeff(L,i,k-1)), Bk);
    2014          28 :     for (j=k+1; j<lx-1; j++)
    2015           0 :       for (i=j+1; i<lx; i++)
    2016           0 :         gcoeff(L,i,j) = gdiv(gmul(p1,gcoeff(L,i,j)), Bk);
    2017             :   }
    2018             :   else
    2019             :   {
    2020          56 :     gcoeff(L,k,k-1) = gen_0;
    2021          56 :     for (i=k+1; i<lx; i++)
    2022             :     {
    2023           0 :       gcoeff(L,i,k) = gcoeff(L,i,k-1);
    2024           0 :       gcoeff(L,i,k-1) = gen_0;
    2025             :     }
    2026          56 :     gel(B,k) = gel(B,k-1); fl[k] = 1; fl[k-1] = 0;
    2027             :   }
    2028         112 :   return 1;
    2029             : }
    2030             : 
    2031             : static void
    2032         168 : incrementalGSgen(GEN x, GEN L, GEN B, long k, GEN fl)
    2033             : {
    2034         168 :   GEN u = NULL; /* gcc -Wall */
    2035             :   long i, j;
    2036         420 :   for (j = 1; j <= k; j++)
    2037         252 :     if (j==k || fl[j])
    2038             :     {
    2039         252 :       u = gcoeff(x,k,j);
    2040         252 :       if (!is_extscalar_t(typ(u))) pari_err_TYPE("incrementalGSgen",u);
    2041         336 :       for (i=1; i<j; i++)
    2042          84 :         if (fl[i])
    2043             :         {
    2044          84 :           u = gsub(gmul(gel(B,i+1),u), gmul(gcoeff(L,k,i),gcoeff(L,j,i)));
    2045          84 :           u = gdiv(u, gel(B,i));
    2046             :         }
    2047         252 :       gcoeff(L,k,j) = u;
    2048             :     }
    2049         168 :   if (gequal0(u)) gel(B,k+1) = gel(B,k);
    2050             :   else
    2051             :   {
    2052         112 :     gel(B,k+1) = gcoeff(L,k,k); gcoeff(L,k,k) = gen_1; fl[k] = 1;
    2053             :   }
    2054         168 : }
    2055             : 
    2056             : static GEN
    2057         168 : lllgramallgen(GEN x, long flag)
    2058             : {
    2059         168 :   long lx = lg(x), i, j, k, l, n;
    2060             :   pari_sp av;
    2061             :   GEN B, L, h, fl;
    2062             :   int flc;
    2063             : 
    2064         168 :   n = lx-1; if (n<=1) return lll_trivial(x,flag);
    2065          84 :   if (lgcols(x) != lx) pari_err_DIM("lllgramallgen");
    2066             : 
    2067          84 :   fl = cgetg(lx, t_VECSMALL);
    2068             : 
    2069          84 :   av = avma;
    2070          84 :   B = scalarcol_shallow(gen_1, lx);
    2071          84 :   L = cgetg(lx,t_MAT);
    2072         252 :   for (j=1; j<lx; j++) { gel(L,j) = zerocol(n); fl[j] = 0; }
    2073             : 
    2074          84 :   h = matid(n);
    2075         252 :   for (i=1; i<lx; i++)
    2076         168 :     incrementalGSgen(x, L, B, i, fl);
    2077          84 :   flc = 0;
    2078          84 :   for(k=2;;)
    2079             :   {
    2080         196 :     if (REDgen(k, k-1, h, L, gel(B,k))) flc = 1;
    2081         196 :     if (do_SWAPgen(h, L, B, k, fl, &flc)) { if (k > 2) k--; }
    2082             :     else
    2083             :     {
    2084          84 :       for (l=k-2; l>=1; l--)
    2085           0 :         if (REDgen(k, l, h, L, gel(B,l+1))) flc = 1;
    2086          84 :       if (++k > n) break;
    2087             :     }
    2088         112 :     if (gc_needed(av,1))
    2089             :     {
    2090           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"lllgramallgen");
    2091           0 :       gerepileall(av,3,&B,&L,&h);
    2092             :     }
    2093             :   }
    2094         140 :   k=1; while (k<lx && !fl[k]) k++;
    2095          84 :   return lll_finish(h,k-1,flag);
    2096             : }
    2097             : 
    2098             : static int
    2099       57688 : RgM_square(GEN x) { long l = lg(x); return l == 1 || l == lgcols(x); }
    2100             : static GEN
    2101         168 : lllallgen(GEN x, long flag)
    2102             : {
    2103         168 :   pari_sp av = avma;
    2104         168 :   if (!(flag & LLL_GRAM)) x = gram_matrix(x);
    2105          84 :   else if (!RgM_square(x)) pari_err_DIM("qflllgram");
    2106         168 :   return gerepilecopy(av, lllgramallgen(x, flag));
    2107             : }
    2108             : GEN
    2109          42 : lllgen(GEN x) { return lllallgen(x, LLL_IM); }
    2110             : GEN
    2111          42 : lllkerimgen(GEN x) { return lllallgen(x, LLL_ALL); }
    2112             : GEN
    2113          42 : lllgramgen(GEN x)  { return lllallgen(x, LLL_IM|LLL_GRAM); }
    2114             : GEN
    2115          42 : lllgramkerimgen(GEN x)  { return lllallgen(x, LLL_ALL|LLL_GRAM); }
    2116             : 
    2117             : static GEN
    2118       56139 : lllall(GEN x, long flag)
    2119       56139 : { pari_sp av = avma; return gerepilecopy(av, ZM_lll(x, LLLDFT, flag)); }
    2120             : GEN
    2121       14683 : lllint(GEN x) { return lllall(x, LLL_IM); }
    2122             : GEN
    2123          35 : lllkerim(GEN x) { return lllall(x, LLL_ALL); }
    2124             : GEN
    2125       41379 : lllgramint(GEN x)
    2126       41379 : { if (!RgM_square(x)) pari_err_DIM("qflllgram");
    2127       41379 :   return lllall(x, LLL_IM | LLL_GRAM); }
    2128             : GEN
    2129          35 : lllgramkerim(GEN x)
    2130          35 : { if (!RgM_square(x)) pari_err_DIM("qflllgram");
    2131          35 :   return lllall(x, LLL_ALL | LLL_GRAM); }
    2132             : 
    2133             : GEN
    2134     1255019 : lllfp(GEN x, double D, long flag)
    2135             : {
    2136     1255019 :   long n = lg(x)-1;
    2137     1255019 :   pari_sp av = avma;
    2138             :   GEN h;
    2139     1255019 :   if (n <= 1) return lll_trivial(x,flag);
    2140     1171050 :   if ((flag & LLL_GRAM) && !RgM_square(x)) pari_err_DIM("qflllgram");
    2141     1171036 :   h = ZM_lll(RgM_rescale_to_int(x), D, flag);
    2142     1171008 :   return gerepilecopy(av, h);
    2143             : }
    2144             : 
    2145             : GEN
    2146       15967 : lllgram(GEN x) { return lllfp(x,LLLDFT,LLL_GRAM|LLL_IM); }
    2147             : GEN
    2148     1165545 : lll(GEN x) { return lllfp(x,LLLDFT,LLL_IM); }
    2149             : 
    2150             : GEN
    2151         301 : qflll0(GEN x, long flag)
    2152             : {
    2153         301 :   if (typ(x) != t_MAT) pari_err_TYPE("qflll",x);
    2154         301 :   switch(flag)
    2155             :   {
    2156          49 :     case 0: return lll(x);
    2157          63 :     case 1: RgM_check_ZM(x,"qflll"); return lllint(x);
    2158          49 :     case 2: RgM_check_ZM(x,"qflll"); return lllintpartial(x);
    2159           7 :     case 3: RgM_check_ZM(x,"qflll"); return lllall(x, LLL_INPLACE);
    2160          49 :     case 4: RgM_check_ZM(x,"qflll"); return lllkerim(x);
    2161          42 :     case 5: return lllkerimgen(x);
    2162          42 :     case 8: return lllgen(x);
    2163           0 :     default: pari_err_FLAG("qflll");
    2164             :   }
    2165             :   return NULL; /* LCOV_EXCL_LINE */
    2166             : }
    2167             : 
    2168             : GEN
    2169         245 : qflllgram0(GEN x, long flag)
    2170             : {
    2171         245 :   if (typ(x) != t_MAT) pari_err_TYPE("qflllgram",x);
    2172         245 :   switch(flag)
    2173             :   {
    2174          63 :     case 0: return lllgram(x);
    2175          49 :     case 1: RgM_check_ZM(x,"qflllgram"); return lllgramint(x);
    2176          49 :     case 4: RgM_check_ZM(x,"qflllgram"); return lllgramkerim(x);
    2177          42 :     case 5: return lllgramkerimgen(x);
    2178          42 :     case 8: return lllgramgen(x);
    2179           0 :     default: pari_err_FLAG("qflllgram");
    2180             :   }
    2181             :   return NULL; /* LCOV_EXCL_LINE */
    2182             : }
    2183             : 
    2184             : /********************************************************************/
    2185             : /**                                                                **/
    2186             : /**                   INTEGRAL KERNEL (LLL REDUCED)                **/
    2187             : /**                                                                **/
    2188             : /********************************************************************/
    2189             : static GEN
    2190          70 : kerint0(GEN M)
    2191             : {
    2192             :   /* return ZM_lll(M, LLLDFT, LLL_KER); */
    2193          70 :   GEN U, H = ZM_hnflll(M,&U,1);
    2194          70 :   long d = lg(M)-lg(H);
    2195          70 :   if (!d) return cgetg(1, t_MAT);
    2196          70 :   return ZM_lll(vecslice(U,1,d), LLLDFT, LLL_INPLACE);
    2197             : }
    2198             : GEN
    2199          42 : kerint(GEN M)
    2200             : {
    2201          42 :   pari_sp av = avma;
    2202          42 :   return gerepilecopy(av, kerint0(M));
    2203             : }
    2204             : /* OBSOLETE: use kerint */
    2205             : GEN
    2206          28 : matkerint0(GEN M, long flag)
    2207             : {
    2208          28 :   pari_sp av = avma;
    2209          28 :   if (typ(M) != t_MAT) pari_err_TYPE("matkerint",M);
    2210          28 :   M = Q_primpart(M);
    2211          28 :   RgM_check_ZM(M, "kerint");
    2212          28 :   switch(flag)
    2213             :   {
    2214          28 :     case 0:
    2215          28 :     case 1: return gerepilecopy(av, kerint0(M));
    2216           0 :     default: pari_err_FLAG("matkerint");
    2217             :   }
    2218             :   return NULL; /* LCOV_EXCL_LINE */
    2219             : }

Generated by: LCOV version 1.14