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 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.10.0 lcov report (development 21059-cbe0d6a) Lines: 422 462 91.3 %
Date: 2017-09-22 06:24:58 Functions: 28 30 93.3 %
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       22041 : lll_trivial(GEN x, long flag)
      23             : {
      24             :   GEN y;
      25       22041 :   if (lg(x) == 1)
      26             :   { /* dim x = 0 */
      27        5495 :     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       16546 :   if (gequal0(gel(x,1)))
      34             :   {
      35          63 :     if (flag & LLL_KER) return matid(1);
      36          63 :     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       16483 :   if (flag & LLL_INPLACE) return gcopy(x);
      42        6858 :   if (flag & LLL_KER) return cgetg(1,t_MAT);
      43        6858 :   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     1056465 : lll_get_im(GEN h, long k)
      53             : {
      54     1056465 :   ulong mask = h[0] & ~LGBITS;
      55     1056465 :   long l = lg(h) - k;
      56     1056465 :   h += k; h[0] = mask | evallg(l);
      57     1056465 :   return h;
      58             : }
      59             : 
      60             : /* k = dim Kernel */
      61             : static GEN
      62     1056493 : lll_finish(GEN h, long k, long flag)
      63             : {
      64             :   GEN g;
      65     1056493 :   if (flag & LLL_KER) { setlg(h,k+1); return h; }
      66     1056465 :   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             : /********************************************************************/
      72             : /**                                                                **/
      73             : /**                   FPLLL (adapted from D. Stehle's code)        **/
      74             : /**                                                                **/
      75             : /********************************************************************/
      76             : /* Babai() and fplll() are a conversion to libpari API and data types
      77             :    of the file proved.c in fplll-1.3 by Damien Stehle'.
      78             : 
      79             :   Copyright 2005, 2006 Damien Stehle'.
      80             : 
      81             :   This program is free software; you can redistribute it and/or modify it
      82             :   under the terms of the GNU General Public License as published by the
      83             :   Free Software Foundation; either version 2 of the License, or (at your
      84             :   option) any later version.
      85             : 
      86             :   This program implements ideas from the paper "Floating-point LLL Revisited",
      87             :   by Phong Nguyen and Damien Stehle', in the Proceedings of Eurocrypt'2005,
      88             :   Springer-Verlag; and was partly inspired by Shoup's NTL library:
      89             :   http://www.shoup.net/ntl/
      90             : */
      91             : 
      92             : /***********************************************/
      93             : /* Babai's Nearest Plane algorithm (iterative) */
      94             : /***********************************************/
      95             : /* Size-reduces b_kappa using mu_{i,j} and r_{i,j} for j<=i <kappa
      96             : Updates B (kappa); computes mu_{kappa,j}, r_{kappa,j} for j<=kappa, and s(kappa)
      97             : mu, r, s updated in place (affrr).
      98             : */
      99             : static long
     100    15547446 : Babai(pari_sp av, long kappa, GEN *pG, GEN *pB, GEN *pU, GEN mu, GEN r, GEN s,
     101             :       long a, long zeros, long maxG, long n, GEN eta, GEN halfplus1, long prec)
     102             : {
     103    15547446 :   pari_sp av0 = avma;
     104    15547446 :   GEN B = *pB, G = *pG, U = *pU, tmp, rtmp, ztmp;
     105    15547446 :   long k, aa = (a > zeros)? a : zeros+1;
     106    15547446 :   GEN maxmu = gen_0, max2mu = gen_0;
     107             :   /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
     108    15547446 :   const long d = U ? lg(U)-1: 0;
     109             : 
     110    15547446 :   if (gc_needed(av,2))
     111             :   {
     112        3987 :     if(DEBUGMEM>1) pari_warn(warnmem,"Babai[0], a=%ld", aa);
     113        3987 :     gerepileall(av,U?3:2,&B,&G,&U);
     114             :   }
     115             :   for (;;) {
     116    24323592 :     int go_on = 0;
     117             :     GEN max3mu;
     118             :     long i, j;
     119             : 
     120    24323592 :     if (gc_needed(av0,2))
     121             :     {
     122          24 :       if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
     123          24 :       gerepileall(av,U?5:4,&B,&G,&maxmu,&max2mu,&U);
     124             :     }
     125             :     /* Step2: compute the GSO for stage kappa */
     126    24323592 :     max3mu = max2mu;
     127    24323592 :     max2mu = maxmu;
     128    24323592 :     maxmu = real_0(prec);
     129   116568169 :     for (j=aa; j<kappa; j++)
     130             :     {
     131    92244577 :       pari_sp btop = avma;
     132    92244577 :       k = zeros+1;
     133    92244577 :       if (j > k)
     134             :       {
     135    75450030 :         tmp  = mulrr(gmael(mu,j,k), gmael(r,kappa,k));
     136    75450030 :         rtmp = subir(gmael(G,kappa,j), tmp);
     137   588192299 :         for (k++; k<j; k++)
     138             :         {
     139   512742269 :           tmp  = mulrr(gmael(mu,j,k), gmael(r,kappa,k));
     140   512742269 :           rtmp = subrr(rtmp,tmp);
     141             :         }
     142    75450030 :         affrr(rtmp, gmael(r,kappa,j));
     143             :       }
     144             :       else
     145    16794547 :         affir(gmael(G,kappa,j), gmael(r,kappa,j));
     146    92244577 :       affrr(divrr(gmael(r,kappa,j), gmael(r,j,j)), gmael(mu,kappa,j));
     147    92244577 :       if (abscmprr(maxmu, gmael(mu,kappa,j))<0)
     148    39585832 :         maxmu = gmael(mu,kappa,j);
     149    92244577 :       avma = btop;
     150             :     }
     151    24323592 :     maxmu = absr(maxmu);
     152    24323592 :     if (typ(max3mu)==t_REAL && abscmprr(max3mu, shiftr(max2mu, 5))<=0)
     153             :     {
     154           0 :       *pB = B; *pG = G; *pU = U;
     155           0 :       if (DEBUGLEVEL>5) err_printf("prec too low\n");
     156           0 :       return kappa;
     157             :     }
     158             : 
     159             :     /* Step3--5: compute the X_j's  */
     160   164251465 :     for (j=kappa-1; j>zeros; j--)
     161             :     {
     162   139927873 :       tmp = gmael(mu,kappa,j);
     163   139927873 :       if (abscmprr(tmp, eta) <= 0) continue; /* (essentially) size-reduced */
     164             : 
     165    23740957 :       if (gc_needed(av0,2))
     166             :       {
     167         205 :         if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
     168         205 :         gerepileall(av,U?5:4,&B,&G,&maxmu,&max2mu,&U);
     169             :       }
     170    23740957 :       go_on = 1;
     171             :       /* we consider separately the case |X| = 1 */
     172    23740957 :       if (abscmprr(tmp, halfplus1) <= 0)
     173             :       {
     174    16883921 :         if (signe(tmp) > 0) { /* in this case, X = 1 */
     175     8740398 :           pari_sp btop = avma;
     176    60063282 :           for (k=zeros+1; k<j; k++)
     177    51322884 :             affrr(subrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
     178     8740398 :           avma = btop;
     179             : 
     180   156462768 :           for (i=1; i<=n; i++)
     181   147722370 :             gmael(B,kappa,i) = subii(gmael(B,kappa,i), gmael(B,j,i));
     182   107734883 :           for (i=1; i<=d; i++)
     183    98994485 :             gmael(U,kappa,i) = subii(gmael(U,kappa,i), gmael(U,j,i));
     184     8740398 :           btop = avma;
     185     8740398 :           ztmp = subii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
     186     8740398 :           ztmp = addii(gmael(G,kappa,kappa), ztmp);
     187     8740398 :           gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
     188    69050479 :           for (i=1; i<=j; i++)
     189    60310081 :             gmael(G,kappa,i) = subii(gmael(G,kappa,i), gmael(G,j,i));
     190    64125977 :           for (i=j+1; i<kappa; i++)
     191    55385579 :             gmael(G,kappa,i) = subii(gmael(G,kappa,i), gmael(G,i,j));
     192    50455809 :           for (i=kappa+1; i<=maxG; i++)
     193    41715411 :             gmael(G,i,kappa) = subii(gmael(G,i,kappa), gmael(G,i,j));
     194             :         } else { /* otherwise X = -1 */
     195     8143523 :           pari_sp btop = avma;
     196    59054704 :           for (k=zeros+1; k<j; k++)
     197    50911181 :             affrr(addrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
     198     8143523 :           avma = btop;
     199             : 
     200   153660805 :           for (i=1; i<=n; i++)
     201   145517282 :             gmael(B,kappa,i) = addii(gmael(B,kappa,i), gmael(B,j,i));
     202   104717959 :           for (i=1; i<=d; i++)
     203    96574436 :             gmael(U,kappa,i) = addii(gmael(U,kappa,i),gmael(U,j,i));
     204     8143523 :           btop = avma;
     205     8143523 :           ztmp = addii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
     206     8143523 :           ztmp = addii(gmael(G,kappa,kappa), ztmp);
     207     8143523 :           gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
     208    67344940 :           for (i=1; i<=j; i++)
     209    59201417 :             gmael(G,kappa,i) = addii(gmael(G,kappa,i), gmael(G,j,i));
     210    63228867 :           for (i=j+1; i<kappa; i++)
     211    55085344 :             gmael(G,kappa,i) = addii(gmael(G,kappa,i), gmael(G,i,j));
     212    49855240 :           for (i=kappa+1; i<=maxG; i++)
     213    41711717 :             gmael(G,i,kappa) = addii(gmael(G,i,kappa), gmael(G,i,j));
     214             :         }
     215    16883921 :         continue;
     216             :       }
     217             :       /* we have |X| >= 2 */
     218     6857036 :       ztmp = roundr_safe(tmp);
     219     6857036 :       if (lgefint(ztmp) == 3)
     220             :       {
     221     6459549 :         pari_sp btop = avma;
     222     6459549 :         ulong xx = ztmp[2]; /* X fits in an ulong */
     223     6459549 :         if (signe(ztmp) > 0) /* = xx */
     224             :         {
     225    13179141 :           for (k=zeros+1; k<j; k++)
     226             :           {
     227     9932384 :             rtmp = subrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k)));
     228     9932384 :             affrr(rtmp, gmael(mu,kappa,k));
     229             :           }
     230     3246757 :           avma = btop;
     231    49781220 :           for (i=1; i<=n; i++)
     232    46534463 :             gmael(B,kappa,i) = submuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
     233    28300026 :           for (i=1; i<=d; i++)
     234    25053269 :             gmael(U,kappa,i) = submuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
     235     3246757 :           btop = avma;
     236     3246757 :           ztmp = shifti(muliu(gmael(G,kappa,j), xx), 1);
     237     3246757 :           ztmp = subii(mulii(gmael(G,j,j), sqru(xx)), ztmp);
     238     3246757 :           ztmp = addii(gmael(G,kappa,kappa), ztmp);
     239     3246757 :           gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
     240    16751650 :           for (i=1; i<=j; i++)
     241    13504893 :             gmael(G,kappa,i) = submuliu_inplace(gmael(G,kappa,i), gmael(G,j,i), xx);
     242    24434235 :           for (i=j+1; i<kappa; i++)
     243    21187478 :             gmael(G,kappa,i) = submuliu_inplace(gmael(G,kappa,i), gmael(G,i,j), xx);
     244     9708924 :           for (i=kappa+1; i<=maxG; i++)
     245     6462167 :             gmael(G,i,kappa) = submuliu_inplace(gmael(G,i,kappa), gmael(G,i,j), xx);
     246             :         }
     247             :         else /* = -xx */
     248             :         {
     249    13052405 :           for (k=zeros+1; k<j; k++)
     250             :           {
     251     9839613 :             rtmp = addrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k)));
     252     9839613 :             affrr(rtmp, gmael(mu,kappa,k));
     253             :           }
     254     3212792 :           avma = btop;
     255    49750055 :           for (i=1; i<=n; i++)
     256    46537263 :             gmael(B,kappa,i) = addmuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
     257    27478770 :           for (i=1; i<=d; i++)
     258    24265978 :             gmael(U,kappa,i) = addmuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
     259     3212792 :           btop = avma;
     260     3212792 :           ztmp = shifti(muliu(gmael(G,kappa,j), xx), 1);
     261     3212792 :           ztmp = addii(mulii(gmael(G,j,j), sqru(xx)), ztmp);
     262     3212792 :           ztmp = addii(gmael(G,kappa,kappa), ztmp);
     263     3212792 :           gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
     264    16405673 :           for (i=1; i<=j; i++)
     265    13192881 :             gmael(G,kappa,i) = addmuliu_inplace(gmael(G,kappa,i), gmael(G,j,i), xx);
     266    24389495 :           for (i=j+1; i<kappa; i++)
     267    21176703 :             gmael(G,kappa,i) = addmuliu_inplace(gmael(G,kappa,i), gmael(G,i,j), xx);
     268     9389091 :           for (i=kappa+1; i<=maxG; i++)
     269     6176299 :             gmael(G,i,kappa) = addmuliu_inplace(gmael(G,i,kappa), gmael(G,i,j), xx);
     270             :         }
     271             :       }
     272             :       else
     273             :       {
     274      397487 :         GEN tmp2  = itor(ztmp,prec);
     275      397487 :         long e = expo(tmp2)-prec2nbits(prec);
     276      397487 :         GEN X = shifti(trunc2nr(tmp2, -e), e);
     277      397487 :         pari_sp btop = avma;
     278             : 
     279     2731316 :         for (k=zeros+1; k<j; k++)
     280             :         {
     281     2333829 :           rtmp = subrr(gmael(mu,kappa,k), mulir(ztmp, gmael(mu,j,k)));
     282     2333829 :           affrr(rtmp, gmael(mu,kappa,k));
     283             :         }
     284      397487 :         avma = btop;
     285     7346260 :         for (i=1; i<=n; i++)
     286     6948773 :           gmael(B,kappa,i) = submulii(gmael(B,kappa,i), gmael(B,j,i), X);
     287     1290209 :         for (i=1; i<=d; i++)
     288      892722 :           gmael(U,kappa,i) = submulii(gmael(U,kappa,i), gmael(U,j,i), X);
     289      397487 :         btop = avma;
     290      397487 :         ztmp = shifti(mulii(gmael(G,kappa,j), X), 1);
     291      397487 :         ztmp = subii(mulii(gmael(G,j,j), sqri(X)), ztmp);
     292      397487 :         ztmp = addii(gmael(G,kappa,kappa), ztmp);
     293      397487 :         gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
     294     3128803 :         for (i=1; i<=j; i++)
     295     2731316 :           gmael(G,kappa,i) = submulii(gmael(G,kappa,i), gmael(G,j,i), X);
     296     2936507 :         for (   ; i<kappa; i++)
     297     2539020 :           gmael(G,kappa,i) = submulii(gmael(G,kappa,i), gmael(G,i,j), X);
     298      494118 :         for (i=kappa+1; i<=maxG; i++)
     299       96631 :           gmael(G,i,kappa) = submulii(gmael(G,i,kappa), gmael(G,i,j), X);
     300             :       }
     301             :     }
     302    24323592 :     if (!go_on) break; /* Anything happened? */
     303     8776146 :     aa = zeros+1;
     304     8776146 :   }
     305             : 
     306    15547446 :   affir(gmael(G,kappa,kappa), gel(s,zeros+1));
     307             :   /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
     308    15547446 :   av = avma;
     309    93118211 :   for (k=zeros+1; k<=kappa-2; k++)
     310             :   {
     311    77570765 :     tmp = subrr(gel(s,k), mulrr(gmael(mu,kappa,k), gmael(r,kappa,k)));
     312    77570765 :     affrr(tmp, gel(s,k+1));
     313             :   }
     314    15547446 :   *pB = B; *pG = G; *pU = U; avma = av;
     315    15547446 :   return 0;
     316             : }
     317             : 
     318             : static void
     319    20581767 : rotate(GEN mu, long kappa2, long kappa, long d)
     320             : {
     321             :   long i, j;
     322    20581767 :   pari_sp av = avma;
     323    20581767 :   GEN mutmp = leafcopy(gel(mu,kappa2));
     324    67304718 :   for (i=kappa2; i>kappa; i--)
     325    46722951 :     for (j=1;j<=d;j++) gmael(mu,i,j) = gmael(mu,i-1,j);
     326    20581767 :   for (j=1;j<=d;j++)   gmael(mu,kappa,j) = gel(mutmp,j);
     327    20581767 :   avma = av;
     328    20581767 : }
     329             : 
     330             : /* ****************** */
     331             : /* The LLL Algorithm  */
     332             : /* ****************** */
     333             : 
     334             : /* LLL-reduces the integer matrix(ces) (G,B,U)? "in place" */
     335             : static GEN
     336     1147391 : fplll(GEN *ptrB, GEN *ptrU, GEN *ptrr, double DELTA, double ETA, long flag, long prec)
     337             : {
     338     1147391 :   const long gram = flag & LLL_GRAM; /*Gram matrix*/
     339     1147391 :   const long keepfirst = flag & LLL_KEEP_FIRST; /*never swap with first vector*/
     340             :   pari_sp av, av2;
     341             :   long kappa, kappa2, d, n, i, j, zeros, kappamax, maxG, bab;
     342             :   GEN G, mu, r, s, tmp, SPtmp, alpha;
     343     1147391 :   GEN delta = dbltor(DELTA), eta = dbltor(ETA), halfplus1 = dbltor(1.5);
     344     1147391 :   const long triangular = 0;
     345             :   pari_timer T;
     346     1147391 :   GEN B = *ptrB, U;
     347     1147391 :   long cnt = 0;
     348             : 
     349     1147391 :   d = lg(B)-1;
     350     1147391 :   if (gram)
     351             :   {
     352       23756 :     G = B;
     353       23756 :     n = d;
     354       23756 :     B = cgetg(1, t_VECSMALL); /* dummy */
     355             :   }
     356             :   else
     357             :   {
     358     1123635 :     G = zeromatcopy(d,d);
     359     1123635 :     n = nbrows(B);
     360             :   }
     361     1147391 :   U = *ptrU; /* NULL if inplace */
     362             : 
     363     1147391 :   if(DEBUGLEVEL>=4)
     364             :   {
     365           0 :     timer_start(&T);
     366           0 :     err_printf("Entering L^2: LLL-parameters (%P.3f,%.3Pf), working precision %d words\n",delta,eta, prec);
     367             :   }
     368             : 
     369     1147391 :   mu = cgetg(d+1, t_MAT);
     370     1147391 :   r  = cgetg(d+1, t_MAT);
     371     1147391 :   s  = cgetg(d+1, t_VEC);
     372     4648642 :   for (j = 1; j <= d; j++)
     373             :   {
     374     3501251 :     GEN M = cgetg(d+1, t_COL), R = cgetg(d+1, t_COL);
     375     3501251 :     gel(mu,j)= M;
     376     3501251 :     gel(r,j) = R;
     377     3501251 :     gel(s,j) = cgetr(prec);
     378    18570358 :     for (i = 1; i <= d; i++) {
     379    15069107 :       gel(R,i) = cgetr(prec);
     380    15069107 :       gel(M,i) = cgetr(prec);
     381             :     }
     382             :   }
     383     1147391 :   SPtmp = zerovec(d+1);
     384     1147391 :   alpha = cgetg(d+1, t_VECSMALL);
     385     1147391 :   av = avma;
     386             : 
     387             :   /* Step2: Initializing the main loop */
     388     1147391 :   kappamax = 1;
     389     1147391 :   i = 1;
     390     1147391 :   maxG = d; /* later updated to kappamax if (!gram) */
     391             : 
     392             :   do {
     393     1151577 :     if (!gram) gmael(G,i,i) = ZV_dotsquare(gel(B,i));
     394     1151577 :     affir(gmael(G,i,i), gmael(r,i,i));
     395     1151577 :   } while (signe(gmael(G,i,i)) == 0 && (++i <=d));
     396     1147391 :   zeros = i-1; /* all vectors B[i] with i <= zeros are zero vectors */
     397     1147391 :   kappa = i;
     398     1147391 :   if (zeros < d) affir(gmael(G,zeros+1,zeros+1), gmael(r,zeros+1,zeros+1));
     399     1147391 :   for (i=zeros+1; i<=d; i++) alpha[i]=1;
     400             : 
     401    17842228 :   while (++kappa <= d)
     402             :   {
     403    15547446 :     if (kappa>kappamax)
     404             :     {
     405     2349674 :       if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
     406     2349674 :       kappamax = kappa;
     407     2349674 :       if (!gram) {
     408     9609173 :         for (i=zeros+1; i<=kappa; i++)
     409     7373083 :           gmael(G,kappa,i) = ZV_dotproduct(gel(B,kappa), gel(B,i));
     410     2236090 :         maxG = kappamax;
     411             :       }
     412             :     }
     413             :     /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
     414    29037453 :     bab = Babai(av, kappa, &G,&B,&U, mu,r,s, alpha[kappa], zeros, maxG,
     415    13490007 :       gram? 0 : ((triangular && kappamax <= n) ? kappamax: n),
     416             :       eta, halfplus1, prec);
     417    15547446 :     if (bab) {*ptrB=(gram?G:B); *ptrU=U; return NULL; }
     418             : 
     419    15547446 :     av2 = avma;
     420    31082145 :     if ((keepfirst && kappa == 2) ||
     421    15534699 :         cmprr(mulrr(gmael(r,kappa-1,kappa-1), delta), gel(s,kappa-1)) <= 0)
     422             :     { /* Step4: Success of Lovasz's condition */
     423     9803060 :       alpha[kappa] = kappa;
     424     9803060 :       tmp = mulrr(gmael(mu,kappa,kappa-1), gmael(r,kappa,kappa-1));
     425     9803060 :       affrr(subrr(gel(s,kappa-1), tmp), gmael(r,kappa,kappa));
     426     9803060 :       avma = av2;
     427             :     }
     428             :     else
     429             :     { /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
     430     5744386 :       if (DEBUGLEVEL>=4 && kappa==kappamax && signe(gel(s,kappa-1)))
     431           0 :         if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", expo(gel(s,1))); }
     432     5744386 :       kappa2 = kappa;
     433             :       do {
     434    13225121 :         kappa--;
     435    13225121 :         if (kappa<zeros+2 + (keepfirst ? 1: 0)) break;
     436    10474199 :         tmp = mulrr(gmael(r,kappa-1,kappa-1), delta);
     437    10474199 :       } while (cmprr(gel(s,kappa-1), tmp) <=0 );
     438     5744386 :       avma = av2;
     439             : 
     440    18969507 :       for (i=kappa; i<kappa2; i++)
     441    13225121 :         if (kappa <= alpha[i]) alpha[i] = kappa;
     442     5744386 :       for (i=kappa2; i>kappa; i--) alpha[i] = alpha[i-1];
     443    17998296 :       for (i=kappa2+1; i<=kappamax; i++)
     444    12253910 :         if (kappa < alpha[i]) alpha[i] = kappa;
     445     5744386 :       alpha[kappa] = kappa;
     446             : 
     447             :       /* Step6: Update the mu's and r's */
     448     5744386 :       rotate(mu,kappa2,kappa,d);
     449     5744386 :       rotate(r,kappa2,kappa,d);
     450     5744386 :       affrr(gel(s,kappa), gmael(r,kappa,kappa));
     451             : 
     452             :       /* Step7: Update B, G, U */
     453     5744386 :       if (!gram) rotate(B,kappa2,kappa,n);
     454     5744386 :       if (U) rotate(U,kappa2,kappa,d);
     455             : 
     456     5744386 :       for (i=1; i<=kappa2; i++) gel(SPtmp,i) = gmael(G,kappa2,i);
     457     5744386 :       for (i=kappa2+1; i<=maxG; i++) gel(SPtmp,i) = gmael(G,i,kappa2);
     458    18969507 :       for (i=kappa2; i>kappa; i--)
     459             :       {
     460    13225121 :         for (j=1; j<kappa; j++) gmael(G,i,j) = gmael(G,i-1,j);
     461    13225121 :         gmael(G,i,kappa) = gel(SPtmp,i-1);
     462    13225121 :         for (j=kappa+1; j<=i; j++) gmael(G,i,j) = gmael(G,i-1,j-1);
     463    13225121 :         for (j=kappa2+1; j<=maxG; j++) gmael(G,j,i) = gmael(G,j,i-1);
     464             :       }
     465     5744386 :       for (i=1; i<kappa; i++) gmael(G,kappa,i) = gel(SPtmp,i);
     466     5744386 :       gmael(G,kappa,kappa) = gel(SPtmp,kappa2);
     467     5744386 :       for (i=kappa2+1; i<=maxG; i++) gmael(G,i,kappa) = gel(SPtmp,i);
     468             : 
     469             :       /* Step8: Prepare the next loop iteration */
     470     5744386 :       if (kappa == zeros+1 && !signe(gmael(G,kappa,kappa)))
     471             :       {
     472       27349 :         zeros++; kappa++;
     473       27349 :         affir(gmael(G,kappa,kappa), gmael(r,kappa,kappa));
     474             :       }
     475             :     }
     476             :   }
     477             : 
     478     1147391 :   if (DEBUGLEVEL>=4) timer_printf(&T,"LLL");
     479     1147391 :   if (ptrr) *ptrr = RgM_diagonal_shallow(r);
     480     1147391 :   if (!U)
     481             :   {
     482       80903 :     if (zeros) {
     483           0 :       if (gram) {
     484           0 :         G = lll_get_im(G, zeros);
     485           0 :         d -= zeros;
     486           0 :         for (i = 1; i <= d; i++) gel(G,i) = lll_get_im(gel(G,i), zeros);
     487             :       }
     488             :       else
     489           0 :         B = lll_get_im(B, zeros);
     490             :     }
     491             :   }
     492     1066488 :   else if (flag & (LLL_IM|LLL_KER|LLL_ALL))
     493     1056430 :     U = lll_finish(U, zeros, flag);
     494     1147391 :   if (gram)
     495             :   {
     496       23756 :     if (U) return U;
     497           0 :     for (i = 1; i <= d; i++)
     498           0 :       for (j = i+1; j <= d; j++) gmael(G,i,j) = gmael(G,j,i);
     499           0 :     return G;
     500             :   }
     501     1123635 :   return U? U: B;
     502             : }
     503             : 
     504             : /* Assume x a ZM, if ptB != NULL, set it to Gram-Schmidt (squared) norms */
     505             : GEN
     506     1167158 : ZM_lll_norms(GEN x, double DELTA, long flag, GEN *B)
     507             : {
     508     1167158 :   pari_sp ltop = avma;
     509     1167158 :   const long compat = flag & LLL_COMPATIBLE;
     510     1167158 :   const double ETA = 0.51;
     511     1167158 :   long p, n = lg(x)-1;
     512             :   GEN U;
     513     1167158 :   if (n <= 1) return lll_trivial(x, flag);
     514     1147391 :   x = RgM_shallowcopy(x);
     515     1147391 :   U = (flag & LLL_INPLACE)? NULL: matid(n);
     516     1147391 :   for (p = compat? DEFAULTPREC: LOWDEFAULTPREC;;)
     517             :   {
     518     1147391 :     GEN m = fplll(&x, &U, B, DELTA, ETA, flag, p);
     519     1147391 :     if (m) return m;
     520           0 :     if (compat)
     521           0 :       p += DEFAULTPREC-2;
     522             :     else
     523           0 :       incrprec(p);
     524           0 :     gerepileall(ltop, U? 2: 1, &x, &U);
     525           0 :   }
     526             : }
     527             : 
     528             : /********************************************************************/
     529             : /**                                                                **/
     530             : /**                        LLL OVER K[X]                           **/
     531             : /**                                                                **/
     532             : /********************************************************************/
     533             : static int
     534         378 : pslg(GEN x)
     535             : {
     536             :   long tx;
     537         378 :   if (gequal0(x)) return 2;
     538         336 :   tx = typ(x); return is_scalar_t(tx)? 3: lg(x);
     539             : }
     540             : 
     541             : static int
     542         147 : REDgen(long k, long l, GEN h, GEN L, GEN B)
     543             : {
     544         147 :   GEN q, u = gcoeff(L,k,l);
     545             :   long i;
     546             : 
     547         147 :   if (pslg(u) < pslg(B)) return 0;
     548             : 
     549         105 :   q = gneg(gdeuc(u,B));
     550         105 :   gel(h,k) = gadd(gel(h,k), gmul(q,gel(h,l)));
     551         105 :   for (i=1; i<l; i++) gcoeff(L,k,i) = gadd(gcoeff(L,k,i), gmul(q,gcoeff(L,l,i)));
     552         105 :   gcoeff(L,k,l) = gadd(gcoeff(L,k,l), gmul(q,B)); return 1;
     553             : }
     554             : 
     555             : static int
     556         147 : do_SWAPgen(GEN h, GEN L, GEN B, long k, GEN fl, int *flc)
     557             : {
     558             :   GEN p1, la, la2, Bk;
     559             :   long ps1, ps2, i, j, lx;
     560             : 
     561         147 :   if (!fl[k-1]) return 0;
     562             : 
     563         105 :   la = gcoeff(L,k,k-1); la2 = gsqr(la);
     564         105 :   Bk = gel(B,k);
     565         105 :   if (fl[k])
     566             :   {
     567          42 :     GEN q = gadd(la2, gmul(gel(B,k-1),gel(B,k+1)));
     568          42 :     ps1 = pslg(gsqr(Bk));
     569          42 :     ps2 = pslg(q);
     570          42 :     if (ps1 <= ps2 && (ps1 < ps2 || !*flc)) return 0;
     571          21 :     *flc = (ps1 != ps2);
     572          21 :     gel(B,k) = gdiv(q, Bk);
     573             :   }
     574             : 
     575          84 :   swap(gel(h,k-1), gel(h,k)); lx = lg(L);
     576          84 :   for (j=1; j<k-1; j++) swap(gcoeff(L,k-1,j), gcoeff(L,k,j));
     577          84 :   if (fl[k])
     578             :   {
     579          21 :     for (i=k+1; i<lx; i++)
     580             :     {
     581           0 :       GEN t = gcoeff(L,i,k);
     582           0 :       p1 = gsub(gmul(gel(B,k+1),gcoeff(L,i,k-1)), gmul(la,t));
     583           0 :       gcoeff(L,i,k) = gdiv(p1, Bk);
     584           0 :       p1 = gadd(gmul(la,gcoeff(L,i,k-1)), gmul(gel(B,k-1),t));
     585           0 :       gcoeff(L,i,k-1) = gdiv(p1, Bk);
     586             :     }
     587             :   }
     588          63 :   else if (!gequal0(la))
     589             :   {
     590          21 :     p1 = gdiv(la2, Bk);
     591          21 :     gel(B,k+1) = gel(B,k) = p1;
     592          21 :     for (i=k+2; i<=lx; i++) gel(B,i) = gdiv(gmul(p1,gel(B,i)),Bk);
     593          21 :     for (i=k+1; i<lx; i++)
     594           0 :       gcoeff(L,i,k-1) = gdiv(gmul(la,gcoeff(L,i,k-1)), Bk);
     595          21 :     for (j=k+1; j<lx-1; j++)
     596           0 :       for (i=j+1; i<lx; i++)
     597           0 :         gcoeff(L,i,j) = gdiv(gmul(p1,gcoeff(L,i,j)), Bk);
     598             :   }
     599             :   else
     600             :   {
     601          42 :     gcoeff(L,k,k-1) = gen_0;
     602          42 :     for (i=k+1; i<lx; i++)
     603             :     {
     604           0 :       gcoeff(L,i,k) = gcoeff(L,i,k-1);
     605           0 :       gcoeff(L,i,k-1) = gen_0;
     606             :     }
     607          42 :     B[k] = B[k-1]; fl[k] = 1; fl[k-1] = 0;
     608             :   }
     609          84 :   return 1;
     610             : }
     611             : 
     612             : static void
     613         126 : incrementalGSgen(GEN x, GEN L, GEN B, long k, GEN fl)
     614             : {
     615         126 :   GEN u = NULL; /* gcc -Wall */
     616             :   long i, j, tu;
     617         315 :   for (j=1; j<=k; j++)
     618         189 :     if (j==k || fl[j])
     619             :     {
     620         189 :       u = gcoeff(x,k,j); tu = typ(u);
     621         189 :       if (! is_extscalar_t(tu)) pari_err_TYPE("incrementalGSgen",u);
     622         252 :       for (i=1; i<j; i++)
     623          63 :         if (fl[i])
     624             :         {
     625          63 :           u = gsub(gmul(gel(B,i+1),u), gmul(gcoeff(L,k,i),gcoeff(L,j,i)));
     626          63 :           u = gdiv(u, gel(B,i));
     627             :         }
     628         189 :       gcoeff(L,k,j) = u;
     629             :     }
     630         126 :   if (gequal0(u)) B[k+1] = B[k];
     631             :   else
     632             :   {
     633          84 :     gel(B,k+1) = gcoeff(L,k,k); gcoeff(L,k,k) = gen_1; fl[k] = 1;
     634             :   }
     635         126 : }
     636             : 
     637             : static GEN
     638         126 : lllgramallgen(GEN x, long flag)
     639             : {
     640         126 :   long lx = lg(x), i, j, k, l, n;
     641             :   pari_sp av;
     642             :   GEN B, L, h, fl;
     643             :   int flc;
     644             : 
     645         126 :   n = lx-1; if (n<=1) return lll_trivial(x,flag);
     646          63 :   if (lgcols(x) != lx) pari_err_DIM("lllgramallgen");
     647             : 
     648          63 :   fl = cgetg(lx, t_VECSMALL);
     649             : 
     650          63 :   av = avma;
     651          63 :   B = scalarcol_shallow(gen_1, lx);
     652          63 :   L = cgetg(lx,t_MAT);
     653          63 :   for (j=1; j<lx; j++) { gel(L,j) = zerocol(n); fl[j] = 0; }
     654             : 
     655          63 :   h = matid(n);
     656         189 :   for (i=1; i<lx; i++)
     657         126 :     incrementalGSgen(x, L, B, i, fl);
     658          63 :   flc = 0;
     659          63 :   for(k=2;;)
     660             :   {
     661         147 :     if (REDgen(k, k-1, h, L, gel(B,k))) flc = 1;
     662         147 :     if (do_SWAPgen(h, L, B, k, fl, &flc)) { if (k > 2) k--; }
     663             :     else
     664             :     {
     665          63 :       for (l=k-2; l>=1; l--)
     666           0 :         if (REDgen(k, l, h, L, gel(B,l+1))) flc = 1;
     667          63 :       if (++k > n) break;
     668             :     }
     669          84 :     if (gc_needed(av,1))
     670             :     {
     671           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"lllgramallgen");
     672           0 :       gerepileall(av,3,&B,&L,&h);
     673             :     }
     674          84 :   }
     675          63 :   k=1; while (k<lx && !fl[k]) k++;
     676          63 :   return lll_finish(h,k-1,flag);
     677             : }
     678             : 
     679             : static GEN
     680         126 : lllallgen(GEN x, long flag)
     681             : {
     682         126 :   pari_sp av = avma;
     683         126 :   if ((flag & LLL_GRAM) == 0) x = gram_matrix(x);
     684         126 :   return gerepilecopy(av, lllgramallgen(x, flag));
     685             : }
     686             : GEN
     687          42 : lllgen(GEN x) { return lllallgen(x, LLL_IM); }
     688             : GEN
     689          42 : lllkerimgen(GEN x) { return lllallgen(x, LLL_ALL); }
     690             : GEN
     691           0 : lllgramgen(GEN x)  { return lllallgen(x, LLL_IM|LLL_GRAM); }
     692             : GEN
     693          42 : lllgramkerimgen(GEN x)  { return lllallgen(x, LLL_ALL|LLL_GRAM); }
     694             : 
     695             : static GEN
     696       26914 : lllall(GEN x, long flag)
     697             : {
     698       26914 :   pari_sp av = avma;
     699       26914 :   return gerepilecopy(av, ZM_lll(x, LLLDFT, flag));
     700             : }
     701             : GEN
     702        6860 : lllint(GEN x) { return lllall(x, LLL_IM); }
     703             : GEN
     704          35 : lllkerim(GEN x) { return lllall(x, LLL_ALL); }
     705             : GEN
     706       19984 : lllgramint(GEN x) { return lllall(x, LLL_IM | LLL_GRAM); }
     707             : GEN
     708          35 : lllgramkerim(GEN x) { return lllall(x, LLL_ALL | LLL_GRAM); }
     709             : 
     710             : GEN
     711      224991 : lllfp(GEN x, double D, long flag)
     712             : {
     713      224991 :   long n = lg(x)-1;
     714      224991 :   pari_sp av = avma;
     715             :   GEN h;
     716      224991 :   if (n <= 1) return lll_trivial(x,flag);
     717      222780 :   h = ZM_lll(RgM_rescale_to_int(x), D, flag);
     718      222752 :   return gerepilecopy(av, h);
     719             : }
     720             : 
     721             : GEN
     722          56 : lllgram(GEN x) { return lllfp(x,LLLDFT,LLL_GRAM|LLL_IM); }
     723             : GEN
     724      214877 : lll(GEN x) { return lllfp(x,LLLDFT,LLL_IM); }
     725             : 
     726             : GEN
     727         294 : qflll0(GEN x, long flag)
     728             : {
     729         294 :   if (typ(x) != t_MAT) pari_err_TYPE("qflll",x);
     730         294 :   switch(flag)
     731             :   {
     732          49 :     case 0: return lll(x);
     733          63 :     case 1: RgM_check_ZM(x,"qflll"); return lllint(x);
     734          49 :     case 2: RgM_check_ZM(x,"qflll"); return lllintpartial(x);
     735          49 :     case 4: RgM_check_ZM(x,"qflll"); return lllkerim(x);
     736          42 :     case 5: return lllkerimgen(x);
     737          42 :     case 8: return lllgen(x);
     738           0 :     default: pari_err_FLAG("qflll");
     739             :   }
     740             :   return NULL; /* LCOV_EXCL_LINE */
     741             : }
     742             : 
     743             : GEN
     744         196 : qflllgram0(GEN x, long flag)
     745             : {
     746         196 :   if (typ(x) != t_MAT) pari_err_TYPE("qflllgram",x);
     747         196 :   switch(flag)
     748             :   {
     749          56 :     case 0: return lllgram(x);
     750          49 :     case 1: RgM_check_ZM(x,"qflllgram"); return lllgramint(x);
     751          49 :     case 4: RgM_check_ZM(x,"qflllgram"); return lllgramkerim(x);
     752          42 :     case 5: return lllgramkerimgen(x);
     753           0 :     case 8: return lllgramgen(x);
     754           0 :     default: pari_err_FLAG("qflllgram");
     755             :   }
     756             :   return NULL; /* LCOV_EXCL_LINE */
     757             : }
     758             : 
     759             : /********************************************************************/
     760             : /**                                                                **/
     761             : /**                   INTEGRAL KERNEL (LLL REDUCED)                **/
     762             : /**                                                                **/
     763             : /********************************************************************/
     764             : static GEN
     765          28 : kerint0(GEN M)
     766             : {
     767          28 :   GEN U, H = ZM_hnfall(M,&U,1);
     768          28 :   long d = lg(M)-lg(H);
     769          28 :   if (!d) return cgetg(1, t_MAT);
     770          28 :   return ZM_lll(vecslice(U,1,d), LLLDFT, LLL_INPLACE);
     771             : }
     772             : GEN
     773           0 : kerint(GEN M)
     774             : {
     775           0 :   pari_sp av = avma;
     776           0 :   return gerepilecopy(av, kerint0(M));
     777             : }
     778             : /* OBSOLETE: use kerint */
     779             : GEN
     780          28 : matkerint0(GEN M, long flag)
     781             : {
     782          28 :   pari_sp av = avma;
     783          28 :   if (typ(M) != t_MAT) pari_err_TYPE("matkerint",M);
     784          28 :   M = Q_primpart(M);
     785          28 :   RgM_check_ZM(M, "kerint");
     786          28 :   switch(flag)
     787             :   {
     788             :     case 0:
     789          28 :     case 1: return gerepilecopy(av, kerint0(M));
     790           0 :     default: pari_err_FLAG("matkerint");
     791             :   }
     792             :   return NULL; /* LCOV_EXCL_LINE */
     793             : }

Generated by: LCOV version 1.11