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.12.1 lcov report (development 25406-bf255ab81b) Lines: 444 480 92.5 %
Date: 2020-06-04 05:59:24 Functions: 31 33 93.9 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2008  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : #include "pari.h"
      15             : #include "paripriv.h"
      16             : 
      17             : /* default quality ratio for LLL */
      18             : static const double LLLDFT = 0.99;
      19             : 
      20             : /* assume flag & (LLL_KER|LLL_IM|LLL_ALL). LLL_INPLACE implies LLL_IM */
      21             : static GEN
      22       34276 : lll_trivial(GEN x, long flag)
      23             : {
      24             :   GEN y;
      25       34276 :   if (lg(x) == 1)
      26             :   { /* dim x = 0 */
      27        7799 :     if (! (flag & LLL_ALL)) return cgetg(1,t_MAT);
      28          28 :     y=cgetg(3,t_VEC);
      29          28 :     gel(y,1) = cgetg(1,t_MAT);
      30          28 :     gel(y,2) = cgetg(1,t_MAT); return y;
      31             :   }
      32             :   /* dim x = 1 */
      33       26477 :   if (gequal0(gel(x,1)))
      34             :   {
      35          84 :     if (flag & LLL_KER) return matid(1);
      36          84 :     if (flag & (LLL_IM|LLL_INPLACE)) return cgetg(1,t_MAT);
      37          28 :     y = cgetg(3,t_VEC);
      38          28 :     gel(y,1) = matid(1);
      39          28 :     gel(y,2) = cgetg(1,t_MAT); return y;
      40             :   }
      41       26393 :   if (flag & LLL_INPLACE) return gcopy(x);
      42        9719 :   if (flag & LLL_KER) return cgetg(1,t_MAT);
      43        9719 :   if (flag & LLL_IM)  return matid(1);
      44          28 :   y=cgetg(3,t_VEC);
      45          28 :   gel(y,1) = cgetg(1,t_MAT);
      46          28 :   gel(y,2) = (flag & LLL_GRAM)? gcopy(x): matid(1);
      47          28 :   return y;
      48             : }
      49             : 
      50             : /* vecslice(h,#h-k,#h) in place. Works for t_MAT, t_VEC/t_COL */
      51             : static GEN
      52     2559444 : lll_get_im(GEN h, long k)
      53             : {
      54     2559444 :   ulong mask = h[0] & ~LGBITS;
      55     2559444 :   long l = lg(h) - k;
      56     2559444 :   h += k; h[0] = mask | evallg(l);
      57     2559444 :   return h;
      58             : }
      59             : 
      60             : /* k = dim Kernel */
      61             : static GEN
      62     2559458 : lll_finish(GEN h, long k, long flag)
      63             : {
      64             :   GEN g;
      65     2559458 :   if (flag & LLL_KER) { setlg(h,k+1); return h; }
      66     2559430 :   if (flag & LLL_IM) return lll_get_im(h, k);
      67          70 :   g = vecslice(h,1,k);
      68          70 :   return mkvec2(g, lll_get_im(h, k));
      69             : }
      70             : 
      71             : INLINE GEN
      72     4114917 : mulshift(GEN y, GEN z, long e)
      73             : {
      74     4114917 :   long ly = lgefint(y), lz;
      75             :   pari_sp av;
      76             :   GEN t;
      77     4114917 :   if (ly == 2) return gen_0;
      78       97925 :   lz = lgefint(z);
      79       97925 :   av = avma; (void)new_chunk(ly+lz+nbits2lg(e)); /* HACK */
      80       97925 :   t = mulii(z, y);
      81       97925 :   set_avma(av); return shifti(t, e);
      82             : }
      83             : 
      84             : INLINE GEN
      85    26300512 : submulshift(GEN x, GEN y, GEN z, long e)
      86             : {
      87    26300512 :   long lx = lgefint(x), ly, lz;
      88             :   pari_sp av;
      89             :   GEN t;
      90    26300512 :   if (!e) return submulii(x, y, z);
      91    25550931 :   if (lx == 2) { t = mulshift(y, z, e); togglesign(t); return t; }
      92    21436014 :   ly = lgefint(y);
      93    21436014 :   if (ly == 2) return icopy(x);
      94    20281248 :   lz = lgefint(z);
      95    20281248 :   av = avma; (void)new_chunk(lx+ly+lz+nbits2lg(e)); /* HACK */
      96    20281248 :   t = shifti(mulii(z, y), e);
      97    20281248 :   set_avma(av); return subii(x, t);
      98             : }
      99             : 
     100             : /********************************************************************/
     101             : /**                                                                **/
     102             : /**                   FPLLL (adapted from D. Stehle's code)        **/
     103             : /**                                                                **/
     104             : /********************************************************************/
     105             : /* Babai() and fplll() are a conversion to libpari API and data types
     106             :    of the file proved.c in fplll-1.3 by Damien Stehle'.
     107             : 
     108             :   Copyright 2005, 2006 Damien Stehle'.
     109             : 
     110             :   This program is free software; you can redistribute it and/or modify it
     111             :   under the terms of the GNU General Public License as published by the
     112             :   Free Software Foundation; either version 2 of the License, or (at your
     113             :   option) any later version.
     114             : 
     115             :   This program implements ideas from the paper "Floating-point LLL Revisited",
     116             :   by Phong Nguyen and Damien Stehle', in the Proceedings of Eurocrypt'2005,
     117             :   Springer-Verlag; and was partly inspired by Shoup's NTL library:
     118             :   http://www.shoup.net/ntl/
     119             : */
     120             : 
     121             : /***********************************************/
     122             : /* Babai's Nearest Plane algorithm (iterative) */
     123             : /***********************************************/
     124             : /* Size-reduces b_kappa using mu_{i,j} and r_{i,j} for j<=i <kappa
     125             : Updates B (kappa); computes mu_{kappa,j}, r_{kappa,j} for j<=kappa, and s(kappa)
     126             : mu, r, s updated in place (affrr).
     127             : */
     128             : static long
     129    38236820 : Babai(pari_sp av, long kappa, GEN *pG, GEN *pB, GEN *pU, GEN mu, GEN r, GEN s,
     130             :       long a, long zeros, long maxG, long n, GEN eta, GEN halfplus1, long prec)
     131             : {
     132    38236820 :   pari_sp av0 = avma;
     133    38236820 :   GEN B = *pB, G = *pG, U = *pU, tmp, rtmp, ztmp;
     134    38236820 :   long k, aa = (a > zeros)? a : zeros+1;
     135    38236820 :   GEN maxmu = gen_0, max2mu = gen_0;
     136             :   /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
     137    38236820 :   const long d = U ? lg(U)-1: 0;
     138             : 
     139    38236820 :   if (gc_needed(av,2))
     140             :   {
     141        4360 :     if(DEBUGMEM>1) pari_warn(warnmem,"Babai[0], a=%ld", aa);
     142        4360 :     gerepileall(av,U?3:2,&B,&G,&U);
     143             :   }
     144    28131315 :   for (;;) {
     145    66368135 :     int go_on = 0;
     146             :     GEN max3mu;
     147             :     long i, j;
     148             : 
     149    66368135 :     if (gc_needed(av0,2))
     150             :     {
     151           7 :       if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
     152           7 :       gerepileall(av,U?5:4,&B,&G,&maxmu,&max2mu,&U);
     153             :     }
     154             :     /* Step2: compute the GSO for stage kappa */
     155    66368135 :     max3mu = max2mu;
     156    66368135 :     max2mu = maxmu;
     157    66368135 :     maxmu = real_0(prec);
     158   231157862 :     for (j=aa; j<kappa; j++)
     159             :     {
     160   164789751 :       pari_sp btop = avma;
     161   164789751 :       k = zeros+1;
     162   164789751 :       if (j > k)
     163             :       {
     164   109154312 :         tmp  = mulrr(gmael(mu,j,k), gmael(r,kappa,k));
     165   109154311 :         rtmp = subir(gmael(G,kappa,j), tmp);
     166   820417266 :         for (k++; k<j; k++)
     167             :         {
     168   711262955 :           tmp  = mulrr(gmael(mu,j,k), gmael(r,kappa,k));
     169   711262955 :           rtmp = subrr(rtmp,tmp);
     170             :         }
     171   109154311 :         affrr(rtmp, gmael(r,kappa,j));
     172             :       }
     173             :       else
     174    55635439 :         affir(gmael(G,kappa,j), gmael(r,kappa,j));
     175   164789722 :       affrr(divrr(gmael(r,kappa,j), gmael(r,j,j)), gmael(mu,kappa,j));
     176   164789594 :       if (abscmprr(maxmu, gmael(mu,kappa,j))<0)
     177    88845581 :         maxmu = gmael(mu,kappa,j);
     178   164789720 :       set_avma(btop);
     179             :     }
     180    66368111 :     maxmu = absr(maxmu);
     181    66368125 :     if (typ(max3mu)==t_REAL && abscmprr(max3mu, shiftr(max2mu, 5))<=0)
     182             :     {
     183           0 :       *pB = B; *pG = G; *pU = U;
     184           0 :       if (DEBUGLEVEL>5) err_printf("prec too low\n");
     185           0 :       return kappa;
     186             :     }
     187             : 
     188             :     /* Step3--5: compute the X_j's  */
     189   299413423 :     for (j=kappa-1; j>zeros; j--)
     190             :     { /* The code below seemingly handles U = NULL, but in this case d = 0 */
     191   233045288 :       tmp = gmael(mu,kappa,j);
     192   233045288 :       if (abscmprr(tmp, eta) <= 0) continue; /* (essentially) size-reduced */
     193             : 
     194    49296920 :       if (gc_needed(av0,2))
     195             :       {
     196          83 :         if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
     197          83 :         gerepileall(av,U?5:4,&B,&G,&maxmu,&max2mu,&U);
     198             :       }
     199    49296920 :       go_on = 1;
     200             :       /* we consider separately the case |X| = 1 */
     201    49296920 :       if (abscmprr(tmp, halfplus1) <= 0)
     202             :       {
     203    24705671 :         if (signe(tmp) > 0) { /* in this case, X = 1 */
     204    13022843 :           pari_sp btop = avma;
     205    83008787 :           for (k=zeros+1; k<j; k++)
     206    69985944 :             affrr(subrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
     207    13022843 :           set_avma(btop);
     208             : 
     209   234394970 :           for (i=1; i<=n; i++)
     210   221372127 :             gmael(B,kappa,i) = subii(gmael(B,kappa,i), gmael(B,j,i));
     211   132621005 :           for (i=1; i<=d; i++)
     212   119598162 :             gmael(U,kappa,i) = subii(gmael(U,kappa,i), gmael(U,j,i));
     213    13022843 :           btop = avma;
     214    13022843 :           ztmp = subii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
     215    13022843 :           ztmp = addii(gmael(G,kappa,kappa), ztmp);
     216    13022843 :           gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
     217    96367588 :           for (i=1; i<=j; i++)
     218    83344745 :             gmael(G,kappa,i) = subii(gmael(G,kappa,i), gmael(G,j,i));
     219    89318860 :           for (i=j+1; i<kappa; i++)
     220    76296017 :             gmael(G,kappa,i) = subii(gmael(G,kappa,i), gmael(G,i,j));
     221    69534055 :           for (i=kappa+1; i<=maxG; i++)
     222    56511212 :             gmael(G,i,kappa) = subii(gmael(G,i,kappa), gmael(G,i,j));
     223             :         } else { /* otherwise X = -1 */
     224    11682828 :           pari_sp btop = avma;
     225    81399161 :           for (k=zeros+1; k<j; k++)
     226    69716333 :             affrr(addrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
     227    11682828 :           set_avma(btop);
     228             : 
     229   229795771 :           for (i=1; i<=n; i++)
     230   218112942 :             gmael(B,kappa,i) = addii(gmael(B,kappa,i), gmael(B,j,i));
     231   127468204 :           for (i=1; i<=d; i++)
     232   115785375 :             gmael(U,kappa,i) = addii(gmael(U,kappa,i),gmael(U,j,i));
     233    11682829 :           btop = avma;
     234    11682829 :           ztmp = addii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
     235    11682826 :           ztmp = addii(gmael(G,kappa,kappa), ztmp);
     236    11682828 :           gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
     237    93304708 :           for (i=1; i<=j; i++)
     238    81621880 :             gmael(G,kappa,i) = addii(gmael(G,kappa,i), gmael(G,j,i));
     239    88035825 :           for (i=j+1; i<kappa; i++)
     240    76352997 :             gmael(G,kappa,i) = addii(gmael(G,kappa,i), gmael(G,i,j));
     241    68289885 :           for (i=kappa+1; i<=maxG; i++)
     242    56607057 :             gmael(G,i,kappa) = addii(gmael(G,i,kappa), gmael(G,i,j));
     243             :         }
     244    24705671 :         continue;
     245             :       }
     246             :       /* we have |X| >= 2 */
     247    24591256 :       ztmp = roundr_safe(tmp);
     248    24591261 :       if (lgefint(ztmp) == 3)
     249             :       {
     250    23901776 :         pari_sp btop = avma;
     251    23901776 :         ulong xx = ztmp[2]; /* X fits in an ulong */
     252    23901776 :         if (signe(ztmp) > 0) /* = xx */
     253             :         {
     254    21730012 :           for (k=zeros+1; k<j; k++)
     255             :           {
     256     9755877 :             rtmp = subrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k)));
     257     9755877 :             affrr(rtmp, gmael(mu,kappa,k));
     258             :           }
     259    11974135 :           set_avma(btop);
     260    75332981 :           for (i=1; i<=n; i++)
     261    63358861 :             gmael(B,kappa,i) = submuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
     262    60342092 :           for (i=1; i<=d; i++)
     263    48367972 :             gmael(U,kappa,i) = submuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
     264    11974120 :           btop = avma;
     265    11974120 :           ztmp = shifti(muliu(gmael(G,kappa,j), xx), 1);
     266    11974133 :           ztmp = subii(mulii(gmael(G,j,j), sqru(xx)), ztmp);
     267    11974136 :           ztmp = addii(gmael(G,kappa,kappa), ztmp);
     268    11974134 :           gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
     269    34155490 :           for (i=1; i<=j; i++)
     270    22181359 :             gmael(G,kappa,i) = submuliu_inplace(gmael(G,kappa,i), gmael(G,j,i), xx);
     271    32645686 :           for (i=j+1; i<kappa; i++)
     272    20671555 :             gmael(G,kappa,i) = submuliu_inplace(gmael(G,kappa,i), gmael(G,i,j), xx);
     273    17763010 :           for (i=kappa+1; i<=maxG; i++)
     274     5788879 :             gmael(G,i,kappa) = submuliu_inplace(gmael(G,i,kappa), gmael(G,i,j), xx);
     275             :         }
     276             :         else /* = -xx */
     277             :         {
     278    21523301 :           for (k=zeros+1; k<j; k++)
     279             :           {
     280     9595660 :             rtmp = addrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k)));
     281     9595660 :             affrr(rtmp, gmael(mu,kappa,k));
     282             :           }
     283    11927641 :           set_avma(btop);
     284    75237413 :           for (i=1; i<=n; i++)
     285    63309794 :             gmael(B,kappa,i) = addmuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
     286    59170595 :           for (i=1; i<=d; i++)
     287    47242976 :             gmael(U,kappa,i) = addmuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
     288    11927619 :           btop = avma;
     289    11927619 :           ztmp = shifti(muliu(gmael(G,kappa,j), xx), 1);
     290    11927638 :           ztmp = addii(mulii(gmael(G,j,j), sqru(xx)), ztmp);
     291    11927633 :           ztmp = addii(gmael(G,kappa,kappa), ztmp);
     292    11927638 :           gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
     293    33695234 :           for (i=1; i<=j; i++)
     294    21767595 :             gmael(G,kappa,i) = addmuliu_inplace(gmael(G,kappa,i), gmael(G,j,i), xx);
     295    32523260 :           for (i=j+1; i<kappa; i++)
     296    20595622 :             gmael(G,kappa,i) = addmuliu_inplace(gmael(G,kappa,i), gmael(G,i,j), xx);
     297    17375720 :           for (i=kappa+1; i<=maxG; i++)
     298     5448082 :             gmael(G,i,kappa) = addmuliu_inplace(gmael(G,i,kappa), gmael(G,i,j), xx);
     299             :         }
     300             :       }
     301             :       else
     302             :       {
     303             :         pari_sp btop;
     304      689485 :         GEN X, tmp2  = itor(ztmp,prec);
     305      689489 :         long e = expo(tmp2) - prec2nbits(prec);
     306             : 
     307      689489 :         X = trunc2nr(tmp2, -e); if (e < 0) { X = shifti(X,e); e = 0; }
     308      689489 :         btop = avma;
     309     5473129 :         for (k=zeros+1; k<j; k++)
     310             :         {
     311     4783640 :           rtmp = subrr(gmael(mu,kappa,k), mulir(ztmp, gmael(mu,j,k)));
     312     4783640 :           affrr(rtmp, gmael(mu,kappa,k));
     313             :         }
     314      689489 :         set_avma(btop);
     315    15497706 :         for (i=1; i<=n; i++)
     316    14808217 :           gmael(B,kappa,i) = submulshift(gmael(B,kappa,i), gmael(B,j,i), X, e);
     317     1690283 :         for (i=1; i<=d; i++)
     318     1000794 :           gmael(U,kappa,i) = submulshift(gmael(U,kappa,i), gmael(U,j,i), X, e);
     319      689489 :         btop = avma;
     320      689489 :         ztmp = shifti(mulii(gmael(G,kappa,j), X), e+1);
     321      689489 :         ztmp = subii(shifti(mulii(gmael(G,j,j), sqri(X)), 2*e), ztmp);
     322      689489 :         ztmp = addii(gmael(G,kappa,kappa), ztmp);
     323      689489 :         gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
     324     6162618 :         for (i=1; i<=j; i++)
     325     5473129 :           gmael(G,kappa,i) = submulshift(gmael(G,kappa,i), gmael(G,j,i), X, e);
     326     5616335 :         for (   ; i<kappa; i++)
     327     4926846 :           gmael(G,kappa,i) = submulshift(gmael(G,kappa,i), gmael(G,i,j), X, e);
     328      781015 :         for (i=kappa+1; i<=maxG; i++)
     329       91526 :           gmael(G,i,kappa) = submulshift(gmael(G,i,kappa), gmael(G,i,j), X, e);
     330             :       }
     331             :     }
     332    66368135 :     if (!go_on) break; /* Anything happened? */
     333    28131315 :     aa = zeros+1;
     334             :   }
     335             : 
     336    38236820 :   affir(gmael(G,kappa,kappa), gel(s,zeros+1));
     337             :   /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
     338    38236810 :   av = avma;
     339   149934022 :   for (k=zeros+1; k<=kappa-2; k++)
     340             :   {
     341   111697205 :     tmp = subrr(gel(s,k), mulrr(gmael(mu,kappa,k), gmael(r,kappa,k)));
     342   111697205 :     affrr(tmp, gel(s,k+1));
     343             :   }
     344    38236817 :   *pB = B; *pG = G; *pU = U; set_avma(av);
     345    38236821 :   return 0;
     346             : }
     347             : 
     348             : static void
     349    90276754 : rotate(GEN mu, long kappa2, long kappa, long d)
     350             : {
     351             :   long i, j;
     352    90276754 :   pari_sp av = avma;
     353    90276754 :   GEN mutmp = leafcopy(gel(mu,kappa2));
     354   216686534 :   for (i=kappa2; i>kappa; i--)
     355  1200076831 :     for (j=1;j<=d;j++) gmael(mu,i,j) = gmael(mu,i-1,j);
     356   528084852 :   for (j=1;j<=d;j++)   gmael(mu,kappa,j) = gel(mutmp,j);
     357    90276760 :   set_avma(av);
     358    90276765 : }
     359             : 
     360             : /* ****************** */
     361             : /* The LLL Algorithm  */
     362             : /* ****************** */
     363             : 
     364             : /* LLL-reduces the integer matrix(ces) (G,B,U)? "in place" */
     365             : static GEN
     366     2678020 : fplll(GEN *ptrB, GEN *ptrU, GEN *ptrr, double DELTA, double ETA, long flag, long prec)
     367             : {
     368     2678020 :   const long gram = flag & LLL_GRAM; /*Gram matrix*/
     369     2678020 :   const long keepfirst = flag & LLL_KEEP_FIRST; /*never swap with first vector*/
     370             :   pari_sp av, av2;
     371             :   long kappa, kappa2, d, n, i, j, zeros, kappamax, maxG, bab;
     372             :   GEN G, mu, r, s, tmp, SPtmp, alpha;
     373     2678020 :   GEN delta = dbltor(DELTA), eta = dbltor(ETA), halfplus1 = dbltor(1.5);
     374     2678020 :   const long triangular = 0;
     375             :   pari_timer T;
     376     2678020 :   GEN B = *ptrB, U;
     377     2678020 :   long cnt = 0;
     378             : 
     379     2678020 :   d = lg(B)-1;
     380     2678020 :   if (gram)
     381             :   {
     382       26706 :     G = B;
     383       26706 :     n = d;
     384       26706 :     B = cgetg(1, t_VECSMALL); /* dummy */
     385             :   }
     386             :   else
     387             :   {
     388     2651314 :     G = zeromatcopy(d,d);
     389     2651314 :     n = nbrows(B);
     390             :   }
     391     2678020 :   U = *ptrU; /* NULL if inplace */
     392             : 
     393     2678020 :   if(DEBUGLEVEL>=4)
     394             :   {
     395           0 :     timer_start(&T);
     396           0 :     err_printf("Entering L^2: LLL-parameters (%P.3f,%.3Pf), working precision %d words\n",delta,eta, prec);
     397             :   }
     398             : 
     399     2678020 :   mu = cgetg(d+1, t_MAT);
     400     2678020 :   r  = cgetg(d+1, t_MAT);
     401     2678020 :   s  = cgetg(d+1, t_VEC);
     402     9567163 :   for (j = 1; j <= d; j++)
     403             :   {
     404     6889143 :     GEN M = cgetg(d+1, t_COL), R = cgetg(d+1, t_COL);
     405     6889143 :     gel(mu,j)= M;
     406     6889143 :     gel(r,j) = R;
     407     6889143 :     gel(s,j) = cgetr(prec);
     408    31269404 :     for (i = 1; i <= d; i++) {
     409    24380261 :       gel(R,i) = cgetr(prec);
     410    24380261 :       gel(M,i) = cgetr(prec);
     411             :     }
     412             :   }
     413     2678020 :   SPtmp = zerovec(d+1);
     414     2678020 :   alpha = cgetg(d+1, t_VECSMALL);
     415     2678020 :   av = avma;
     416             : 
     417             :   /* Step2: Initializing the main loop */
     418     2678020 :   kappamax = 1;
     419     2678020 :   i = 1;
     420     2678020 :   maxG = d; /* later updated to kappamax if (!gram) */
     421             : 
     422             :   do {
     423     2683319 :     if (!gram) gmael(G,i,i) = ZV_dotsquare(gel(B,i));
     424     2683319 :     affir(gmael(G,i,i), gmael(r,i,i));
     425     2683319 :   } while (signe(gmael(G,i,i)) == 0 && (++i <=d));
     426     2678020 :   zeros = i-1; /* all vectors B[i] with i <= zeros are zero vectors */
     427     2678020 :   kappa = i;
     428     2678020 :   if (zeros < d) affir(gmael(G,zeros+1,zeros+1), gmael(r,zeros+1,zeros+1));
     429     9561857 :   for (i=zeros+1; i<=d; i++) alpha[i]=1;
     430             : 
     431    40914838 :   while (++kappa <= d)
     432             :   {
     433    38236818 :     if (kappa>kappamax)
     434             :     {
     435     4205824 :       if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
     436     4205824 :       kappamax = kappa;
     437     4205824 :       if (!gram) {
     438    16102188 :         for (i=zeros+1; i<=kappa; i++)
     439    12024006 :           gmael(G,kappa,i) = ZV_dotproduct(gel(B,kappa), gel(B,i));
     440     4078182 :         maxG = kappamax;
     441             :       }
     442             :     }
     443             :     /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
     444    74129213 :     bab = Babai(av, kappa, &G,&B,&U, mu,r,s, alpha[kappa], zeros, maxG,
     445    35892395 :       gram? 0 : ((triangular && kappamax <= n) ? kappamax: n),
     446             :       eta, halfplus1, prec);
     447    38236822 :     if (bab) {*ptrB=(gram?G:B); *ptrU=U; return NULL; }
     448             : 
     449    38236822 :     av2 = avma;
     450    76453740 :     if ((keepfirst && kappa == 2) ||
     451    38216928 :         cmprr(mulrr(gmael(r,kappa-1,kappa-1), delta), gel(s,kappa-1)) <= 0)
     452             :     { /* Step4: Success of Lovasz's condition */
     453    14763790 :       alpha[kappa] = kappa;
     454    14763790 :       tmp = mulrr(gmael(mu,kappa,kappa-1), gmael(r,kappa,kappa-1));
     455    14763793 :       affrr(subrr(gel(s,kappa-1), tmp), gmael(r,kappa,kappa));
     456    14763792 :       set_avma(av2);
     457             :     }
     458             :     else
     459             :     { /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
     460    23473022 :       if (DEBUGLEVEL>=4 && kappa==kappamax && signe(gel(s,kappa-1)))
     461           0 :         if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", expo(gel(s,1))); }
     462    23473021 :       kappa2 = kappa;
     463             :       do {
     464    34066189 :         kappa--;
     465    34066189 :         if (kappa<zeros+2 + (keepfirst ? 1: 0)) break;
     466    14801297 :         tmp = mulrr(gmael(r,kappa-1,kappa-1), delta);
     467    14801298 :       } while (cmprr(gel(s,kappa-1), tmp) <=0 );
     468    23473022 :       set_avma(av2);
     469             : 
     470    57539221 :       for (i=kappa; i<kappa2; i++)
     471    34066197 :         if (kappa <= alpha[i]) alpha[i] = kappa;
     472    57539221 :       for (i=kappa2; i>kappa; i--) alpha[i] = alpha[i-1];
     473    41510794 :       for (i=kappa2+1; i<=kappamax; i++)
     474    18037770 :         if (kappa < alpha[i]) alpha[i] = kappa;
     475    23473024 :       alpha[kappa] = kappa;
     476             : 
     477             :       /* Step6: Update the mu's and r's */
     478    23473024 :       rotate(mu,kappa2,kappa,d);
     479    23473030 :       rotate(r,kappa2,kappa,d);
     480    23473029 :       affrr(gel(s,kappa), gmael(r,kappa,kappa));
     481             : 
     482             :       /* Step7: Update B, G, U */
     483    23473028 :       if (!gram) rotate(B,kappa2,kappa,n);
     484    23473030 :       if (U) rotate(U,kappa2,kappa,d);
     485             : 
     486   108297163 :       for (i=1; i<=kappa2; i++) gel(SPtmp,i) = gmael(G,kappa2,i);
     487    42477881 :       for (i=kappa2+1; i<=maxG; i++) gel(SPtmp,i) = gmael(G,i,kappa2);
     488    57539232 :       for (i=kappa2; i>kappa; i--)
     489             :       {
     490   105487973 :         for (j=1; j<kappa; j++) gmael(G,i,j) = gmael(G,i-1,j);
     491    34066202 :         gmael(G,i,kappa) = gel(SPtmp,i-1);
     492   106060009 :         for (j=kappa+1; j<=i; j++) gmael(G,i,j) = gmael(G,i-1,j-1);
     493    80223845 :         for (j=kappa2+1; j<=maxG; j++) gmael(G,j,i) = gmael(G,j,i-1);
     494             :       }
     495    50757933 :       for (i=1; i<kappa; i++) gmael(G,kappa,i) = gel(SPtmp,i);
     496    23473030 :       gmael(G,kappa,kappa) = gel(SPtmp,kappa2);
     497    42477881 :       for (i=kappa2+1; i<=maxG; i++) gmael(G,i,kappa) = gel(SPtmp,i);
     498             : 
     499             :       /* Step8: Prepare the next loop iteration */
     500    23473030 :       if (kappa == zeros+1 && !signe(gmael(G,kappa,kappa)))
     501             :       {
     502       35203 :         zeros++; kappa++;
     503       35203 :         affir(gmael(G,kappa,kappa), gmael(r,kappa,kappa));
     504             :       }
     505             :     }
     506             :   }
     507             : 
     508     2678020 :   if (DEBUGLEVEL>=4) timer_printf(&T,"LLL");
     509     2678020 :   if (ptrr) *ptrr = RgM_diagonal_shallow(r);
     510     2678020 :   if (!U)
     511             :   {
     512      103030 :     if (zeros) {
     513          14 :       if (gram) {
     514           0 :         G = lll_get_im(G, zeros);
     515           0 :         d -= zeros;
     516           0 :         for (i = 1; i <= d; i++) gel(G,i) = lll_get_im(gel(G,i), zeros);
     517             :       }
     518             :       else
     519          14 :         B = lll_get_im(B, zeros);
     520             :     }
     521             :   }
     522     2574990 :   else if (flag & (LLL_IM|LLL_KER|LLL_ALL))
     523     2559395 :     U = lll_finish(U, zeros, flag);
     524     2678020 :   if (gram)
     525             :   {
     526       26706 :     if (U) return U;
     527           0 :     for (i = 1; i <= d; i++)
     528           0 :       for (j = i+1; j <= d; j++) gmael(G,i,j) = gmael(G,j,i);
     529           0 :     return G;
     530             :   }
     531     2651314 :   return U? U: B;
     532             : }
     533             : 
     534             : /* Assume x a ZM, if ptB != NULL, set it to Gram-Schmidt (squared) norms */
     535             : GEN
     536     2708964 : ZM_lll_norms(GEN x, double DELTA, long flag, GEN *B)
     537             : {
     538     2708964 :   pari_sp ltop = avma;
     539     2708964 :   const double ETA = 0.51;
     540     2708964 :   long p, n = lg(x)-1;
     541             :   GEN U;
     542     2708964 :   if (n <= 1) return lll_trivial(x, flag);
     543     2678020 :   x = RgM_shallowcopy(x);
     544     2678020 :   U = (flag & LLL_INPLACE)? NULL: matid(n);
     545     2678020 :   for (p = DEFAULTPREC;;)
     546           0 :   {
     547     2678020 :     GEN m = fplll(&x, &U, B, DELTA, ETA, flag, p);
     548     2678020 :     if (m) return m;
     549           0 :     p += DEFAULTPREC-2;
     550           0 :     gerepileall(ltop, U? 2: 1, &x, &U);
     551             :   }
     552             : }
     553             : 
     554             : /********************************************************************/
     555             : /**                                                                **/
     556             : /**                        LLL OVER K[X]                           **/
     557             : /**                                                                **/
     558             : /********************************************************************/
     559             : static int
     560         378 : pslg(GEN x)
     561             : {
     562             :   long tx;
     563         378 :   if (gequal0(x)) return 2;
     564         336 :   tx = typ(x); return is_scalar_t(tx)? 3: lg(x);
     565             : }
     566             : 
     567             : static int
     568         147 : REDgen(long k, long l, GEN h, GEN L, GEN B)
     569             : {
     570         147 :   GEN q, u = gcoeff(L,k,l);
     571             :   long i;
     572             : 
     573         147 :   if (pslg(u) < pslg(B)) return 0;
     574             : 
     575         105 :   q = gneg(gdeuc(u,B));
     576         105 :   gel(h,k) = gadd(gel(h,k), gmul(q,gel(h,l)));
     577         105 :   for (i=1; i<l; i++) gcoeff(L,k,i) = gadd(gcoeff(L,k,i), gmul(q,gcoeff(L,l,i)));
     578         105 :   gcoeff(L,k,l) = gadd(gcoeff(L,k,l), gmul(q,B)); return 1;
     579             : }
     580             : 
     581             : static int
     582         147 : do_SWAPgen(GEN h, GEN L, GEN B, long k, GEN fl, int *flc)
     583             : {
     584             :   GEN p1, la, la2, Bk;
     585             :   long ps1, ps2, i, j, lx;
     586             : 
     587         147 :   if (!fl[k-1]) return 0;
     588             : 
     589         105 :   la = gcoeff(L,k,k-1); la2 = gsqr(la);
     590         105 :   Bk = gel(B,k);
     591         105 :   if (fl[k])
     592             :   {
     593          42 :     GEN q = gadd(la2, gmul(gel(B,k-1),gel(B,k+1)));
     594          42 :     ps1 = pslg(gsqr(Bk));
     595          42 :     ps2 = pslg(q);
     596          42 :     if (ps1 <= ps2 && (ps1 < ps2 || !*flc)) return 0;
     597          21 :     *flc = (ps1 != ps2);
     598          21 :     gel(B,k) = gdiv(q, Bk);
     599             :   }
     600             : 
     601          84 :   swap(gel(h,k-1), gel(h,k)); lx = lg(L);
     602          84 :   for (j=1; j<k-1; j++) swap(gcoeff(L,k-1,j), gcoeff(L,k,j));
     603          84 :   if (fl[k])
     604             :   {
     605          21 :     for (i=k+1; i<lx; i++)
     606             :     {
     607           0 :       GEN t = gcoeff(L,i,k);
     608           0 :       p1 = gsub(gmul(gel(B,k+1),gcoeff(L,i,k-1)), gmul(la,t));
     609           0 :       gcoeff(L,i,k) = gdiv(p1, Bk);
     610           0 :       p1 = gadd(gmul(la,gcoeff(L,i,k-1)), gmul(gel(B,k-1),t));
     611           0 :       gcoeff(L,i,k-1) = gdiv(p1, Bk);
     612             :     }
     613             :   }
     614          63 :   else if (!gequal0(la))
     615             :   {
     616          21 :     p1 = gdiv(la2, Bk);
     617          21 :     gel(B,k+1) = gel(B,k) = p1;
     618          21 :     for (i=k+2; i<=lx; i++) gel(B,i) = gdiv(gmul(p1,gel(B,i)),Bk);
     619          21 :     for (i=k+1; i<lx; i++)
     620           0 :       gcoeff(L,i,k-1) = gdiv(gmul(la,gcoeff(L,i,k-1)), Bk);
     621          21 :     for (j=k+1; j<lx-1; j++)
     622           0 :       for (i=j+1; i<lx; i++)
     623           0 :         gcoeff(L,i,j) = gdiv(gmul(p1,gcoeff(L,i,j)), Bk);
     624             :   }
     625             :   else
     626             :   {
     627          42 :     gcoeff(L,k,k-1) = gen_0;
     628          42 :     for (i=k+1; i<lx; i++)
     629             :     {
     630           0 :       gcoeff(L,i,k) = gcoeff(L,i,k-1);
     631           0 :       gcoeff(L,i,k-1) = gen_0;
     632             :     }
     633          42 :     B[k] = B[k-1]; fl[k] = 1; fl[k-1] = 0;
     634             :   }
     635          84 :   return 1;
     636             : }
     637             : 
     638             : static void
     639         126 : incrementalGSgen(GEN x, GEN L, GEN B, long k, GEN fl)
     640             : {
     641         126 :   GEN u = NULL; /* gcc -Wall */
     642             :   long i, j, tu;
     643         315 :   for (j=1; j<=k; j++)
     644         189 :     if (j==k || fl[j])
     645             :     {
     646         189 :       u = gcoeff(x,k,j); tu = typ(u);
     647         189 :       if (! is_extscalar_t(tu)) pari_err_TYPE("incrementalGSgen",u);
     648         252 :       for (i=1; i<j; i++)
     649          63 :         if (fl[i])
     650             :         {
     651          63 :           u = gsub(gmul(gel(B,i+1),u), gmul(gcoeff(L,k,i),gcoeff(L,j,i)));
     652          63 :           u = gdiv(u, gel(B,i));
     653             :         }
     654         189 :       gcoeff(L,k,j) = u;
     655             :     }
     656         126 :   if (gequal0(u)) B[k+1] = B[k];
     657             :   else
     658             :   {
     659          84 :     gel(B,k+1) = gcoeff(L,k,k); gcoeff(L,k,k) = gen_1; fl[k] = 1;
     660             :   }
     661         126 : }
     662             : 
     663             : static GEN
     664         126 : lllgramallgen(GEN x, long flag)
     665             : {
     666         126 :   long lx = lg(x), i, j, k, l, n;
     667             :   pari_sp av;
     668             :   GEN B, L, h, fl;
     669             :   int flc;
     670             : 
     671         126 :   n = lx-1; if (n<=1) return lll_trivial(x,flag);
     672          63 :   if (lgcols(x) != lx) pari_err_DIM("lllgramallgen");
     673             : 
     674          63 :   fl = cgetg(lx, t_VECSMALL);
     675             : 
     676          63 :   av = avma;
     677          63 :   B = scalarcol_shallow(gen_1, lx);
     678          63 :   L = cgetg(lx,t_MAT);
     679         189 :   for (j=1; j<lx; j++) { gel(L,j) = zerocol(n); fl[j] = 0; }
     680             : 
     681          63 :   h = matid(n);
     682         189 :   for (i=1; i<lx; i++)
     683         126 :     incrementalGSgen(x, L, B, i, fl);
     684          63 :   flc = 0;
     685          63 :   for(k=2;;)
     686             :   {
     687         147 :     if (REDgen(k, k-1, h, L, gel(B,k))) flc = 1;
     688         147 :     if (do_SWAPgen(h, L, B, k, fl, &flc)) { if (k > 2) k--; }
     689             :     else
     690             :     {
     691          63 :       for (l=k-2; l>=1; l--)
     692           0 :         if (REDgen(k, l, h, L, gel(B,l+1))) flc = 1;
     693          63 :       if (++k > n) break;
     694             :     }
     695          84 :     if (gc_needed(av,1))
     696             :     {
     697           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"lllgramallgen");
     698           0 :       gerepileall(av,3,&B,&L,&h);
     699             :     }
     700             :   }
     701         105 :   k=1; while (k<lx && !fl[k]) k++;
     702          63 :   return lll_finish(h,k-1,flag);
     703             : }
     704             : 
     705             : static int
     706       21981 : RgM_square(GEN x) { long l = lg(x); return l == 1 || l == lgcols(x); }
     707             : static GEN
     708         126 : lllallgen(GEN x, long flag)
     709             : {
     710         126 :   pari_sp av = avma;
     711         126 :   if ((flag & LLL_GRAM) == 0) x = gram_matrix(x);
     712          42 :   else if (!RgM_square(x)) pari_err_DIM("qflllgram");
     713         126 :   return gerepilecopy(av, lllgramallgen(x, flag));
     714             : }
     715             : GEN
     716          42 : lllgen(GEN x) { return lllallgen(x, LLL_IM); }
     717             : GEN
     718          42 : lllkerimgen(GEN x) { return lllallgen(x, LLL_ALL); }
     719             : GEN
     720           0 : lllgramgen(GEN x)  { return lllallgen(x, LLL_IM|LLL_GRAM); }
     721             : GEN
     722          42 : lllgramkerimgen(GEN x)  { return lllallgen(x, LLL_ALL|LLL_GRAM); }
     723             : 
     724             : static GEN
     725       28297 : lllall(GEN x, long flag)
     726             : {
     727       28297 :   pari_sp av = avma;
     728       28297 :   if ((flag & LLL_GRAM) && !RgM_square(x)) pari_err_DIM("qflllgram");
     729       28297 :   return gerepilecopy(av, ZM_lll(x, LLLDFT, flag));
     730             : }
     731             : GEN
     732        6861 : lllint(GEN x) { return lllall(x, LLL_IM); }
     733             : GEN
     734          35 : lllkerim(GEN x) { return lllall(x, LLL_ALL); }
     735             : GEN
     736       21366 : lllgramint(GEN x) { return lllall(x, LLL_IM | LLL_GRAM); }
     737             : GEN
     738          35 : lllgramkerim(GEN x) { return lllall(x, LLL_ALL | LLL_GRAM); }
     739             : 
     740             : GEN
     741      229688 : lllfp(GEN x, double D, long flag)
     742             : {
     743      229688 :   long n = lg(x)-1;
     744      229688 :   pari_sp av = avma;
     745             :   GEN h;
     746      229688 :   if (n <= 1) return lll_trivial(x,flag);
     747      226419 :   if ((flag & LLL_GRAM) && !RgM_square(x)) pari_err_DIM("qflllgram");
     748      226405 :   h = ZM_lll(RgM_rescale_to_int(x), D, flag);
     749      226377 :   return gerepilecopy(av, h);
     750             : }
     751             : 
     752             : GEN
     753          63 : lllgram(GEN x) { return lllfp(x,LLLDFT,LLL_GRAM|LLL_IM); }
     754             : GEN
     755      214023 : lll(GEN x) { return lllfp(x,LLLDFT,LLL_IM); }
     756             : 
     757             : GEN
     758         294 : qflll0(GEN x, long flag)
     759             : {
     760         294 :   if (typ(x) != t_MAT) pari_err_TYPE("qflll",x);
     761         294 :   switch(flag)
     762             :   {
     763          49 :     case 0: return lll(x);
     764          63 :     case 1: RgM_check_ZM(x,"qflll"); return lllint(x);
     765          49 :     case 2: RgM_check_ZM(x,"qflll"); return lllintpartial(x);
     766          49 :     case 4: RgM_check_ZM(x,"qflll"); return lllkerim(x);
     767          42 :     case 5: return lllkerimgen(x);
     768          42 :     case 8: return lllgen(x);
     769           0 :     default: pari_err_FLAG("qflll");
     770             :   }
     771             :   return NULL; /* LCOV_EXCL_LINE */
     772             : }
     773             : 
     774             : GEN
     775         203 : qflllgram0(GEN x, long flag)
     776             : {
     777         203 :   if (typ(x) != t_MAT) pari_err_TYPE("qflllgram",x);
     778         203 :   switch(flag)
     779             :   {
     780          63 :     case 0: return lllgram(x);
     781          49 :     case 1: RgM_check_ZM(x,"qflllgram"); return lllgramint(x);
     782          49 :     case 4: RgM_check_ZM(x,"qflllgram"); return lllgramkerim(x);
     783          42 :     case 5: return lllgramkerimgen(x);
     784           0 :     case 8: return lllgramgen(x);
     785           0 :     default: pari_err_FLAG("qflllgram");
     786             :   }
     787             :   return NULL; /* LCOV_EXCL_LINE */
     788             : }
     789             : 
     790             : /********************************************************************/
     791             : /**                                                                **/
     792             : /**                   INTEGRAL KERNEL (LLL REDUCED)                **/
     793             : /**                                                                **/
     794             : /********************************************************************/
     795             : static GEN
     796          28 : kerint0(GEN M)
     797             : {
     798             :   /* return ZM_lll(M, LLLDFT, LLL_KER); */
     799          28 :   GEN U, H = ZM_hnflll(M,&U,1);
     800          28 :   long d = lg(M)-lg(H);
     801          28 :   if (!d) return cgetg(1, t_MAT);
     802          28 :   return ZM_lll(vecslice(U,1,d), LLLDFT, LLL_INPLACE);
     803             : }
     804             : GEN
     805           0 : kerint(GEN M)
     806             : {
     807           0 :   pari_sp av = avma;
     808           0 :   return gerepilecopy(av, kerint0(M));
     809             : }
     810             : /* OBSOLETE: use kerint */
     811             : GEN
     812          28 : matkerint0(GEN M, long flag)
     813             : {
     814          28 :   pari_sp av = avma;
     815          28 :   if (typ(M) != t_MAT) pari_err_TYPE("matkerint",M);
     816          28 :   M = Q_primpart(M);
     817          28 :   RgM_check_ZM(M, "kerint");
     818          28 :   switch(flag)
     819             :   {
     820          28 :     case 0:
     821          28 :     case 1: return gerepilecopy(av, kerint0(M));
     822           0 :     default: pari_err_FLAG("matkerint");
     823             :   }
     824             :   return NULL; /* LCOV_EXCL_LINE */
     825             : }

Generated by: LCOV version 1.13