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.13.0 lcov report (development 25825-edc109b529) Lines: 434 462 93.9 %
Date: 2020-09-21 06:08:33 Functions: 33 35 94.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       35403 : lll_trivial(GEN x, long flag)
      23             : {
      24       35403 :   if (lg(x) == 1)
      25             :   { /* dim x = 0 */
      26        7834 :     if (! (flag & LLL_ALL)) return cgetg(1,t_MAT);
      27          28 :     retmkvec2(cgetg(1,t_MAT), cgetg(1,t_MAT));
      28             :   }
      29             :   /* dim x = 1 */
      30       27569 :   if (gequal0(gel(x,1)))
      31             :   {
      32         742 :     if (flag & LLL_KER) return matid(1);
      33         742 :     if (flag & (LLL_IM|LLL_INPLACE)) return cgetg(1,t_MAT);
      34          28 :     retmkvec2(matid(1), cgetg(1,t_MAT));
      35             :   }
      36       26827 :   if (flag & LLL_INPLACE) return gcopy(x);
      37        9985 :   if (flag & LLL_KER) return cgetg(1,t_MAT);
      38        9985 :   if (flag & LLL_IM)  return matid(1);
      39          28 :   retmkvec2(cgetg(1,t_MAT), (flag & LLL_GRAM)? gcopy(x): matid(1));
      40             : }
      41             : 
      42             : /* vecslice(x,#x-k,#x) in place. Works for t_MAT, t_VEC/t_COL */
      43             : static GEN
      44     2711590 : vectail_inplace(GEN x, long k)
      45             : {
      46     2711590 :   if (!k) return x;
      47       16717 :   x[k] = ((ulong)x[0] & ~LGBITS) | evallg(lg(x) - k);
      48       16717 :   return x + k;
      49             : }
      50             : 
      51             : /* k = dim Kernel */
      52             : static GEN
      53     2727696 : lll_finish(GEN h, long k, long flag)
      54             : {
      55             :   GEN g;
      56     2727696 :   if (!(flag & (LLL_IM|LLL_KER|LLL_ALL|LLL_INPLACE))) return h;
      57     2711618 :   if (flag & (LLL_IM|LLL_INPLACE)) return vectail_inplace(h, k);
      58          98 :   if (flag & LLL_KER) { setlg(h,k+1); return h; }
      59          70 :   g = vecslice(h,1,k); /* done first: vectail_inplace kills h */
      60          70 :   return mkvec2(g, vectail_inplace(h, k));
      61             : }
      62             : 
      63             : INLINE GEN
      64     8043968 : mulshift(GEN y, GEN z, long e)
      65             : {
      66     8043968 :   long ly = lgefint(y), lz;
      67             :   pari_sp av;
      68             :   GEN t;
      69     8043968 :   if (ly == 2) return gen_0;
      70      120951 :   lz = lgefint(z);
      71      120951 :   av = avma; (void)new_chunk(ly+lz+nbits2lg(e)); /* HACK */
      72      120951 :   t = mulii(z, y);
      73      120951 :   set_avma(av); return shifti(t, e);
      74             : }
      75             : 
      76             : INLINE GEN
      77    50307402 : submulshift(GEN x, GEN y, GEN z, long e)
      78             : {
      79    50307402 :   long lx = lgefint(x), ly, lz;
      80             :   pari_sp av;
      81             :   GEN t;
      82    50307402 :   if (!e) return submulii(x, y, z);
      83    49292133 :   if (lx == 2) { t = mulshift(y, z, e); togglesign(t); return t; }
      84    41248165 :   ly = lgefint(y);
      85    41248165 :   if (ly == 2) return icopy(x);
      86    39283417 :   lz = lgefint(z);
      87    39283417 :   av = avma; (void)new_chunk(lx+ly+lz+nbits2lg(e)); /* HACK */
      88    39283417 :   t = shifti(mulii(z, y), e);
      89    39283417 :   set_avma(av); return subii(x, t);
      90             : }
      91             : 
      92             : /********************************************************************/
      93             : /**                                                                **/
      94             : /**                   FPLLL (adapted from D. Stehle's code)        **/
      95             : /**                                                                **/
      96             : /********************************************************************/
      97             : /* n < 10; gerepileall supporting &NULL arguments. Maybe rename and export ? */
      98             : INLINE void
      99        6228 : gc_lll(pari_sp av, int n, ...)
     100             : {
     101             :   int i, j;
     102             :   GEN *gptr[10];
     103        6228 :   va_list a; va_start(a, n);
     104       25364 :   for (i=j=0; i<n; i++)
     105             :   {
     106       19136 :     GEN *x = va_arg(a,GEN*);
     107       19136 :     if (*x) { gptr[j++] = x; *x = (GEN)copy_bin(*x); }
     108             :   }
     109        6228 :   set_avma(av);
     110       21688 :   for (--j; j>=0; j--) *gptr[j] = bin_copy((GENbin*)*gptr[j]);
     111        6228 :   va_end(a);
     112        6228 : }
     113             : /* Babai() and fplll() are a conversion to libpari API and data types
     114             :    of the file proved.c in fplll-1.3 by Damien Stehle'.
     115             : 
     116             :   Copyright 2005, 2006 Damien Stehle'.
     117             : 
     118             :   This program is free software; you can redistribute it and/or modify it
     119             :   under the terms of the GNU General Public License as published by the
     120             :   Free Software Foundation; either version 2 of the License, or (at your
     121             :   option) any later version.
     122             : 
     123             :   This program implements ideas from the paper "Floating-point LLL Revisited",
     124             :   by Phong Nguyen and Damien Stehle', in the Proceedings of Eurocrypt'2005,
     125             :   Springer-Verlag; and was partly inspired by Shoup's NTL library:
     126             :   http://www.shoup.net/ntl/
     127             : */
     128             : 
     129             : /* x t_REAL, |x| >= 1/2. Test whether |x| <= 1.5 */
     130             : static int
     131    62159230 : absrsmall(GEN x)
     132             : {
     133    62159230 :   long e = expo(x);
     134    62159230 :   return (e < 0) || (e == 0 && (ulong)x[2] <= (3UL << (BITS_IN_LONG-2)));
     135             : }
     136             : 
     137             : /***********************************************/
     138             : /* Babai's Nearest Plane algorithm (iterative) */
     139             : /***********************************************/
     140             : /* Size-reduces b_kappa using mu_{i,j} and r_{i,j} for j<=i <kappa
     141             : Updates B (kappa); computes mu_{kappa,j}, r_{kappa,j} for j<=kappa, and s(kappa)
     142             : mu, r, s updated in place (affrr).
     143             : */
     144             : static long
     145    44893299 : Babai(pari_sp av, long kappa, GEN *pG, GEN *pB, GEN *pU, GEN mu, GEN r, GEN s,
     146             :       long a, long zeros, long maxG, GEN eta, long prec)
     147             : {
     148    44893299 :   pari_sp av0 = avma;
     149    44893299 :   GEN G = *pG, B = *pB, U = *pU, tmp, rtmp, ztmp;
     150    44893299 :   long k, d, n, aa = a > zeros? a: zeros+1;
     151    44893299 :   GEN maxmu = gen_0, max2mu = gen_0;
     152             : 
     153             :   /* N.B: we set d = 0 (resp. n = 0) to avoid updating U (resp. B) */
     154    44893299 :   d = U? lg(U)-1: 0;
     155    44893299 :   n = B? nbrows(B): 0;
     156    44893298 :   if (gc_needed(av,2))
     157             :   {
     158        6002 :     if(DEBUGMEM>1) pari_warn(warnmem,"Babai[0], a=%ld", aa);
     159        6002 :     gc_lll(av,3,&G,&B,&U);
     160             :   }
     161    31348097 :   for (;;) {
     162    76241395 :     int go_on = 0;
     163             :     GEN max3mu;
     164             :     long i, j;
     165             : 
     166    76241395 :     if (gc_needed(av0,2))
     167             :     {
     168           1 :       if(DEBUGMEM>1) pari_warn(warnmem,"Babai[1], a=%ld", aa);
     169           1 :       gc_lll(av,5,&G,&B,&U,&maxmu,&max2mu);
     170             :     }
     171             :     /* Step2: compute the GSO for stage kappa */
     172    76241395 :     max3mu = max2mu;
     173    76241395 :     max2mu = maxmu;
     174    76241395 :     maxmu = real_0(prec);
     175   296488511 :     for (j=aa; j<kappa; j++)
     176             :     {
     177   220247148 :       pari_sp btop = avma;
     178   220247148 :       k = zeros+1;
     179   220247148 :       if (j > k)
     180             :       {
     181   159506942 :         tmp  = mulrr(gmael(mu,j,k), gmael(r,kappa,k));
     182   159506942 :         rtmp = subir(gmael(G,kappa,j), tmp);
     183  1202541029 :         for (k++; k<j; k++)
     184             :         {
     185  1043034092 :           tmp  = mulrr(gmael(mu,j,k), gmael(r,kappa,k));
     186  1043034092 :           rtmp = subrr(rtmp,tmp);
     187             :         }
     188   159506937 :         affrr(rtmp, gmael(r,kappa,j));
     189             :       }
     190             :       else
     191    60740206 :         affir(gmael(G,kappa,j), gmael(r,kappa,j));
     192   220247113 :       affrr(divrr(gmael(r,kappa,j), gmael(r,j,j)), gmael(mu,kappa,j));
     193   220246995 :       if (abscmprr(maxmu, gmael(mu,kappa,j))<0)
     194   109369018 :         maxmu = gmael(mu,kappa,j);
     195   220247110 :       set_avma(btop);
     196             :     }
     197    76241363 :     maxmu = absr(maxmu); /* copy needed: we 'affrr' to mu later */
     198    76241378 :     if (typ(max3mu)==t_REAL && abscmprr(max3mu, shiftr(max2mu, 5)) <= 0)
     199             :     { /* precision too low */
     200           0 :       *pG = G; *pB = B; *pU = U; return kappa;
     201             :     }
     202             : 
     203             :     /* Step3--5: compute the X_j's  */
     204   398547201 :     for (j=kappa-1; j>zeros; j--)
     205             :     {
     206   322305814 :       tmp = gmael(mu,kappa,j);
     207   322305814 :       if (abscmprr(tmp, eta) <= 0) continue; /* (essentially) size-reduced */
     208             : 
     209    62159216 :       if (gc_needed(av0,2))
     210             :       {
     211         225 :         if(DEBUGMEM>1) pari_warn(warnmem,"Babai[2], a=%ld, j=%ld", aa,j);
     212         225 :         gc_lll(av,5,&G,&B,&U,&maxmu,&max2mu);
     213             :       }
     214    62159216 :       go_on = 1;
     215             :       /* we consider separately the case |X| = 1 */
     216    62159216 :       if (absrsmall(tmp))
     217             :       {
     218    34744538 :         if (signe(tmp) > 0) { /* in this case, X = 1 */
     219    18060916 :           pari_sp btop = avma;
     220   118206745 :           for (k=zeros+1; k<j; k++)
     221   100145829 :             affrr(subrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
     222    18060916 :           set_avma(btop);
     223             : 
     224   329136970 :           for (i=1; i<=n; i++)
     225   311076055 :             gmael(B,kappa,i) = subii(gmael(B,kappa,i), gmael(B,j,i));
     226   191056941 :           for (i=1; i<=d; i++)
     227   172996026 :             gmael(U,kappa,i) = subii(gmael(U,kappa,i), gmael(U,j,i));
     228    18060915 :           btop = avma;
     229    18060915 :           ztmp = subii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
     230    18060914 :           ztmp = addii(gmael(G,kappa,kappa), ztmp);
     231    18060915 :           gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
     232   136796385 :           for (i=1; i<=j; i++)
     233   118735468 :             gmael(G,kappa,i) = subii(gmael(G,kappa,i), gmael(G,j,i));
     234   128031498 :           for (i=j+1; i<kappa; i++)
     235   109970581 :             gmael(G,kappa,i) = subii(gmael(G,kappa,i), gmael(G,i,j));
     236   100946428 :           for (i=kappa+1; i<=maxG; i++)
     237    82885511 :             gmael(G,i,kappa) = subii(gmael(G,i,kappa), gmael(G,i,j));
     238             :         } else { /* otherwise X = -1 */
     239    16683622 :           pari_sp btop = avma;
     240   116574787 :           for (k=zeros+1; k<j; k++)
     241    99891165 :             affrr(addrr(gmael(mu,kappa,k), gmael(mu,j,k)), gmael(mu,kappa,k));
     242    16683622 :           set_avma(btop);
     243             : 
     244   324114394 :           for (i=1; i<=n; i++)
     245   307430772 :             gmael(B,kappa,i) = addii(gmael(B,kappa,i),gmael(B,j,i));
     246   185562915 :           for (i=1; i<=d; i++)
     247   168879293 :             gmael(U,kappa,i) = addii(gmael(U,kappa,i),gmael(U,j,i));
     248    16683622 :           btop = avma;
     249    16683622 :           ztmp = addii(gmael(G,j,j), shifti(gmael(G,kappa,j), 1));
     250    16683622 :           ztmp = addii(gmael(G,kappa,kappa), ztmp);
     251    16683622 :           gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
     252   133695903 :           for (i=1; i<=j; i++)
     253   117012282 :             gmael(G,kappa,i) = addii(gmael(G,kappa,i), gmael(G,j,i));
     254   126656036 :           for (i=j+1; i<kappa; i++)
     255   109972415 :             gmael(G,kappa,i) = addii(gmael(G,kappa,i), gmael(G,i,j));
     256    99633723 :           for (i=kappa+1; i<=maxG; i++)
     257    82950102 :             gmael(G,i,kappa) = addii(gmael(G,i,kappa), gmael(G,i,j));
     258             :         }
     259    34744538 :         continue;
     260             :       }
     261             :       /* we have |X| >= 2 */
     262    27414692 :       ztmp = roundr_safe(tmp);
     263    27414687 :       if (lgefint(ztmp) == 3)
     264             :       {
     265    26273174 :         pari_sp btop = avma;
     266    26273174 :         ulong xx = ztmp[2]; /* X fits in an ulong */
     267    26273174 :         if (signe(ztmp) > 0) /* = xx */
     268             :         {
     269    25413925 :           for (k=zeros+1; k<j; k++)
     270             :           {
     271    12240965 :             rtmp = subrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k)));
     272    12240965 :             affrr(rtmp, gmael(mu,kappa,k));
     273             :           }
     274    13172960 :           set_avma(btop);
     275    91423624 :           for (i=1; i<=n; i++)
     276    78250687 :             gmael(B,kappa,i) = submuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
     277    73977775 :           for (i=1; i<=d; i++)
     278    60804838 :             gmael(U,kappa,i) = submuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
     279    13172937 :           btop = avma;
     280    13172937 :           ztmp = shifti(muliu(gmael(G,kappa,j), xx), 1);
     281    13172959 :           ztmp = subii(mulii(gmael(G,j,j), sqru(xx)), ztmp);
     282    13172955 :           ztmp = addii(gmael(G,kappa,kappa), ztmp);
     283    13172954 :           gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
     284    39941129 :           for (i=1; i<=j; i++)
     285    26768173 :             gmael(G,kappa,i) = submuliu_inplace(gmael(G,kappa,i), gmael(G,j,i), xx);
     286    39808113 :           for (i=j+1; i<kappa; i++)
     287    26635158 :             gmael(G,kappa,i) = submuliu_inplace(gmael(G,kappa,i), gmael(G,i,j), xx);
     288    20320891 :           for (i=kappa+1; i<=maxG; i++)
     289     7147936 :             gmael(G,i,kappa) = submuliu_inplace(gmael(G,i,kappa), gmael(G,i,j), xx);
     290             :         }
     291             :         else /* = -xx */
     292             :         {
     293    25103647 :           for (k=zeros+1; k<j; k++)
     294             :           {
     295    12003433 :             rtmp = addrr(gmael(mu,kappa,k), mulur(xx, gmael(mu,j,k)));
     296    12003433 :             affrr(rtmp, gmael(mu,kappa,k));
     297             :           }
     298    13100214 :           set_avma(btop);
     299    90958043 :           for (i=1; i<=n; i++)
     300    77857855 :             gmael(B,kappa,i) = addmuliu_inplace(gmael(B,kappa,i), gmael(B,j,i), xx);
     301    72629580 :           for (i=1; i<=d; i++)
     302    59529392 :             gmael(U,kappa,i) = addmuliu_inplace(gmael(U,kappa,i), gmael(U,j,i), xx);
     303    13100188 :           btop = avma;
     304    13100188 :           ztmp = shifti(muliu(gmael(G,kappa,j), xx), 1);
     305    13100213 :           ztmp = addii(mulii(gmael(G,j,j), sqru(xx)), ztmp);
     306    13100210 :           ztmp = addii(gmael(G,kappa,kappa), ztmp);
     307    13100210 :           gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
     308    39403104 :           for (i=1; i<=j; i++)
     309    26302892 :             gmael(G,kappa,i) = addmuliu_inplace(gmael(G,kappa,i), gmael(G,j,i), xx);
     310    39521390 :           for (i=j+1; i<kappa; i++)
     311    26421179 :             gmael(G,kappa,i) = addmuliu_inplace(gmael(G,kappa,i), gmael(G,i,j), xx);
     312    19836746 :           for (i=kappa+1; i<=maxG; i++)
     313     6736535 :             gmael(G,i,kappa) = addmuliu_inplace(gmael(G,i,kappa), gmael(G,i,j), xx);
     314             :         }
     315             :       }
     316             :       else
     317             :       {
     318             :         pari_sp btop;
     319     1141513 :         long e = expi(ztmp) - prec2nbits(prec);
     320     1141518 :         GEN X = ztmp;
     321             : 
     322     1141518 :         if (e <= 0) e = 0; else X = shifti(X, -e);
     323     1141518 :         btop = avma;
     324    10355782 :         for (k=zeros+1; k<j; k++)
     325             :         {
     326     9214264 :           rtmp = subrr(gmael(mu,kappa,k), mulir(ztmp, gmael(mu,j,k)));
     327     9214264 :           affrr(rtmp, gmael(mu,kappa,k));
     328             :         }
     329     1141518 :         set_avma(btop);
     330    29435968 :         for (i=1; i<=n; i++)
     331    28294450 :           gmael(B,kappa,i) = submulshift(gmael(B,kappa,i), gmael(B,j,i), X, e);
     332     3262994 :         for (i=1; i<=d; i++)
     333     2121476 :           gmael(U,kappa,i) = submulshift(gmael(U,kappa,i), gmael(U,j,i), X, e);
     334     1141518 :         btop = avma;
     335     1141518 :         ztmp = shifti(mulii(gmael(G,kappa,j), X), e+1);
     336     1141518 :         ztmp = subii(shifti(mulii(gmael(G,j,j), sqri(X)), 2*e), ztmp);
     337     1141518 :         ztmp = addii(gmael(G,kappa,kappa), ztmp);
     338     1141518 :         gmael(G,kappa,kappa) = gerepileuptoint(btop, ztmp);
     339    11514240 :         for (i=1; i<=j; i++)
     340    10372722 :           gmael(G,kappa,i) = submulshift(gmael(G,kappa,i), gmael(G,j,i), X, e);
     341    10535250 :         for (   ; i<kappa; i++)
     342     9393732 :           gmael(G,kappa,i) = submulshift(gmael(G,kappa,i), gmael(G,i,j), X, e);
     343     1266540 :         for (i=kappa+1; i<=maxG; i++)
     344      125022 :           gmael(G,i,kappa) = submulshift(gmael(G,i,kappa), gmael(G,i,j), X, e);
     345             :       }
     346             :     }
     347    76241387 :     if (!go_on) break; /* Anything happened? */
     348    31348097 :     aa = zeros+1;
     349             :   }
     350             : 
     351    44893290 :   affir(gmael(G,kappa,kappa), gel(s,zeros+1));
     352             :   /* the last s[kappa-1]=r[kappa][kappa] is computed only if kappa increases */
     353    44893286 :   av = avma;
     354   210857438 :   for (k=zeros+1; k<=kappa-2; k++)
     355             :   {
     356   165964159 :     tmp = subrr(gel(s,k), mulrr(gmael(mu,kappa,k), gmael(r,kappa,k)));
     357   165964156 :     affrr(tmp, gel(s,k+1));
     358             :   }
     359    44893279 :   *pG = G; *pB = B; *pU = U; return gc_long(av, 0);
     360             : }
     361             : 
     362             : static void
     363    97162975 : rotate(GEN A, long k2, long k)
     364             : {
     365             :   long i;
     366    97162975 :   GEN B = gel(A,k2);
     367   244178945 :   for (i = k2; i > k; i--) gel(A,i) = gel(A,i-1);
     368    97162975 :   gel(A,k) = B;
     369    97162975 : }
     370             : 
     371             : /* ****************** */
     372             : /* The LLL Algorithm  */
     373             : /* ****************** */
     374             : /* G integral Gram matrix, LLL-reduces (G,B,U) in place [apply base change
     375             :  * transforms to B and U]. If (keepfirst), never swap with first vector.
     376             :  * If G = NULL, we compute the Gram matrix incrementally.
     377             :  * Return -1 on failure, else zeros = dim Kernel (>= 0) */
     378             : static long
     379     2727633 : fplll(GEN *pG, GEN *pB, GEN *pU, GEN *pr, double DELTA, double ETA,
     380             :       long keepfirst, long prec)
     381             : {
     382             :   pari_sp av, av2;
     383     2727633 :   GEN mu, r, s, tmp, SPtmp, alpha, G = *pG, B = *pB, U = *pU;
     384     2727633 :   GEN delta = dbltor(DELTA), eta = dbltor(ETA);
     385     2727633 :   long d, maxG, kappa, kappa2, i, j, zeros, kappamax, incgram = !G, cnt = 0;
     386             : 
     387     2727633 :   if (incgram)
     388             :   { /* incremental Gram matrix */
     389     2700752 :     maxG = 2; d = lg(B)-1;
     390     2700752 :     G = zeromatcopy(d, d);
     391             :   }
     392             :   else
     393       26881 :     maxG = d = lg(G)-1;
     394             : 
     395     2727633 :   mu = cgetg(d+1, t_MAT);
     396     2727633 :   r  = cgetg(d+1, t_MAT);
     397     2727633 :   s  = cgetg(d+1, t_VEC);
     398    10148824 :   for (j = 1; j <= d; j++)
     399             :   {
     400     7421191 :     GEN M = cgetg(d+1, t_COL), R = cgetg(d+1, t_COL);
     401     7421191 :     gel(mu,j)= M;
     402     7421191 :     gel(r,j) = R;
     403     7421191 :     gel(s,j) = cgetr(prec);
     404    37861700 :     for (i = 1; i <= d; i++)
     405             :     {
     406    30440509 :       gel(R,i) = cgetr(prec);
     407    30440509 :       gel(M,i) = cgetr(prec);
     408             :     }
     409             :   }
     410     2727633 :   SPtmp = cgetg(d+1, t_VEC);
     411     2727633 :   alpha = cgetg(d+1, t_VECSMALL);
     412     2727633 :   av = avma;
     413             : 
     414             :   /* Step2: Initializing the main loop */
     415     2727633 :   kappamax = 1;
     416     2727633 :   i = 1;
     417             :   do {
     418     2769810 :     if (incgram) gmael(G,i,i) = ZV_dotsquare(gel(B,i));
     419     2769810 :     affir(gmael(G,i,i), gmael(r,i,i));
     420     2769810 :   } while (!signe(gmael(G,i,i)) && ++i <= d);
     421     2727633 :   zeros = i-1; /* all basis vectors b_i with i <= zeros are zero vectors */
     422     2727633 :   kappa = i;
     423     2727633 :   if (zeros < d) affir(gmael(G,zeros+1,zeros+1), gmael(r,zeros+1,zeros+1));
     424    10099534 :   for (i=zeros+1; i<=d; i++) alpha[i]=1;
     425             : 
     426    47620928 :   while (++kappa <= d)
     427             :   {
     428    44893295 :     if (kappa > kappamax)
     429             :     {
     430     4651381 :       if (DEBUGLEVEL>=4) err_printf("K%ld ",kappa);
     431     4651381 :       kappamax = kappa;
     432     4651381 :       if (incgram)
     433             :       {
     434    18873277 :         for (i=zeros+1; i<=kappa; i++)
     435    14350364 :           gmael(G,kappa,i) = ZV_dotproduct(gel(B,kappa), gel(B,i));
     436     4522913 :         maxG = kappamax;
     437             :       }
     438             :     }
     439             :     /* Step3: Call to the Babai algorithm, mu,r,s updated in place */
     440    44893295 :     if (Babai(av, kappa, &G,&B,&U, mu,r,s, alpha[kappa], zeros, maxG, eta,prec))
     441             :     {
     442           0 :       if (incgram) G = NULL;
     443           0 :       *pG = G; *pB = B; *pU = U; return -1;
     444             :     }
     445    44893289 :     av2 = avma;
     446    89766090 :     if ((keepfirst && kappa == 2) ||
     447    44872800 :         cmprr(mulrr(gmael(r,kappa-1,kappa-1), delta), gel(s,kappa-1)) <= 0)
     448             :     { /* Step4: Success of Lovasz's condition */
     449    19374464 :       alpha[kappa] = kappa;
     450    19374464 :       tmp = mulrr(gmael(mu,kappa,kappa-1), gmael(r,kappa,kappa-1));
     451    19374467 :       affrr(subrr(gel(s,kappa-1), tmp), gmael(r,kappa,kappa));
     452    19374465 :       set_avma(av2); continue;
     453             :     }
     454             :     /* Step5: Find the right insertion index kappa, kappa2 = initial kappa */
     455    25518826 :     if (DEBUGLEVEL>=4 && kappa==kappamax && signe(gel(s,kappa-1)))
     456           0 :       if (++cnt > 20) { cnt = 0; err_printf("(%ld) ", expo(gel(s,1))); }
     457    25518826 :     kappa2 = kappa;
     458             :     do {
     459    40305582 :       kappa--;
     460    40305582 :       if (kappa < zeros+2 + (keepfirst ? 1: 0)) break;
     461    20570389 :       tmp = mulrr(gmael(r,kappa-1,kappa-1), delta);
     462    20570395 :     } while (cmprr(gel(s,kappa-1), tmp) <= 0);
     463    25518832 :     set_avma(av2);
     464             : 
     465    65824417 :     for (i=kappa; i<kappa2; i++)
     466    40305590 :       if (kappa <= alpha[i]) alpha[i] = kappa;
     467    65824417 :     for (i=kappa2; i>kappa; i--) alpha[i] = alpha[i-1];
     468    52888633 :     for (i=kappa2+1; i<=kappamax; i++)
     469    27369806 :       if (kappa < alpha[i]) alpha[i] = kappa;
     470    25518827 :     alpha[kappa] = kappa;
     471             : 
     472             :     /* Step6: Update the mu's and r's */
     473    25518827 :     rotate(mu, kappa2, kappa);
     474    25518830 :     rotate(r, kappa2, kappa);
     475    25518828 :     affrr(gel(s,kappa), gmael(r,kappa,kappa));
     476             : 
     477             :     /* Step7: Update G, B, U */
     478    25518829 :     if (U) rotate(U, kappa2, kappa);
     479    25518829 :     if (B) rotate(B, kappa2, kappa);
     480   132739261 :     for (i=1; i<=kappa2; i++) gel(SPtmp,i) = gmael(G,kappa2,i);
     481    53859246 :     for (   ; i<=maxG;   i++) gel(SPtmp,i) = gmael(G,i,kappa2);
     482    65824419 :     for (i=kappa2; i>kappa; i--)
     483             :     {
     484   149159279 :       for (j=1; j<kappa; j++) gmael(G,i,j) = gmael(G,i-1,j);
     485    40305590 :       gmael(G,i,kappa) = gel(SPtmp,i-1);
     486   136682048 :       for (j=kappa+1; j<=i; j++) gmael(G,i,j) = gmael(G,i-1,j-1);
     487   108622551 :       for (j=kappa2+1; j<=maxG; j++) gmael(G,j,i) = gmael(G,j,i-1);
     488             :     }
     489    66914842 :     for (i=1; i<kappa; i++) gmael(G,kappa,i) = gel(SPtmp,i);
     490    25518829 :     gmael(G,kappa,kappa) = gel(SPtmp,kappa2);
     491    53859246 :     for (i=kappa2+1; i<=maxG; i++) gmael(G,i,kappa) = gel(SPtmp,i);
     492             : 
     493             :     /* Step8: Prepare the next loop iteration */
     494    25518829 :     if (kappa == zeros+1 && !signe(gmael(G,kappa,kappa)))
     495             :     {
     496       63676 :       zeros++; kappa++;
     497       63676 :       affir(gmael(G,kappa,kappa), gmael(r,kappa,kappa));
     498             :     }
     499             :   }
     500     2727633 :   if (pr) *pr = RgM_diagonal_shallow(r);
     501     2727633 :   *pG = G; *pB = B; *pU = U; return zeros; /* success */
     502             : }
     503             : 
     504             : /* Assume x a ZM, if pN != NULL, set it to Gram-Schmidt (squared) norms
     505             :  * The following modes are supported:
     506             :  * - flag & LLL_INPLACE: x a lattice basis, return x*U
     507             :  * - flag & LLL_GRAM: x a Gram matrix / else x a lattice basis; return
     508             :  *     LLL base change matrix U [LLL_IM]
     509             :  *     kernel basis [LLL_KER, non-reduced]
     510             :  *     both [LLL_ALL] */
     511             : GEN
     512     2748628 : ZM_lll_norms(GEN x, double DELTA, long flag, GEN *pN)
     513             : {
     514     2748628 :   pari_sp av = avma;
     515     2748628 :   const double ETA = 0.51;
     516     2748628 :   long p, zeros, n = lg(x)-1;
     517             :   GEN G, B, U;
     518             :   pari_timer T;
     519             : 
     520     2748628 :   if (n <= 1) return lll_trivial(x, flag);
     521     2727632 :   x = RgM_shallowcopy(x);
     522     2727633 :   if (flag & LLL_GRAM)
     523       26881 :   { G = x; B = NULL; U = matid(n); }
     524             :   else
     525     2700752 :   { G = NULL; B = x; U = (flag & LLL_INPLACE)? NULL: matid(n); }
     526     2727633 :   for (p = DEFAULTPREC;; p += EXTRAPREC64)
     527             :   {
     528     2727633 :     if (DEBUGLEVEL>=4)
     529             :     {
     530           0 :       err_printf("Entering L^2: LLL-parameters (%.3f,%.3f), prec = %d\n",
     531             :                  DELTA,ETA, p);
     532           0 :       timer_start(&T);
     533             :     }
     534     2727633 :     zeros = fplll(&G, &B, &U, pN, DELTA, ETA, flag & LLL_KEEP_FIRST, p);
     535     2727633 :     if (zeros >= 0) break;
     536           0 :     gc_lll(av, 3, &G, &B, &U);
     537             :   }
     538     2727633 :   if (DEBUGLEVEL>=4) timer_printf(&T,"LLL");
     539     2727633 :   return lll_finish(U? U: B, zeros, flag);
     540             : }
     541             : 
     542             : /********************************************************************/
     543             : /**                                                                **/
     544             : /**                        LLL OVER K[X]                           **/
     545             : /**                                                                **/
     546             : /********************************************************************/
     547             : static int
     548         378 : pslg(GEN x)
     549             : {
     550             :   long tx;
     551         378 :   if (gequal0(x)) return 2;
     552         336 :   tx = typ(x); return is_scalar_t(tx)? 3: lg(x);
     553             : }
     554             : 
     555             : static int
     556         147 : REDgen(long k, long l, GEN h, GEN L, GEN B)
     557             : {
     558         147 :   GEN q, u = gcoeff(L,k,l);
     559             :   long i;
     560             : 
     561         147 :   if (pslg(u) < pslg(B)) return 0;
     562             : 
     563         105 :   q = gneg(gdeuc(u,B));
     564         105 :   gel(h,k) = gadd(gel(h,k), gmul(q,gel(h,l)));
     565         105 :   for (i=1; i<l; i++) gcoeff(L,k,i) = gadd(gcoeff(L,k,i), gmul(q,gcoeff(L,l,i)));
     566         105 :   gcoeff(L,k,l) = gadd(gcoeff(L,k,l), gmul(q,B)); return 1;
     567             : }
     568             : 
     569             : static int
     570         147 : do_SWAPgen(GEN h, GEN L, GEN B, long k, GEN fl, int *flc)
     571             : {
     572             :   GEN p1, la, la2, Bk;
     573             :   long ps1, ps2, i, j, lx;
     574             : 
     575         147 :   if (!fl[k-1]) return 0;
     576             : 
     577         105 :   la = gcoeff(L,k,k-1); la2 = gsqr(la);
     578         105 :   Bk = gel(B,k);
     579         105 :   if (fl[k])
     580             :   {
     581          42 :     GEN q = gadd(la2, gmul(gel(B,k-1),gel(B,k+1)));
     582          42 :     ps1 = pslg(gsqr(Bk));
     583          42 :     ps2 = pslg(q);
     584          42 :     if (ps1 <= ps2 && (ps1 < ps2 || !*flc)) return 0;
     585          21 :     *flc = (ps1 != ps2);
     586          21 :     gel(B,k) = gdiv(q, Bk);
     587             :   }
     588             : 
     589          84 :   swap(gel(h,k-1), gel(h,k)); lx = lg(L);
     590          84 :   for (j=1; j<k-1; j++) swap(gcoeff(L,k-1,j), gcoeff(L,k,j));
     591          84 :   if (fl[k])
     592             :   {
     593          21 :     for (i=k+1; i<lx; i++)
     594             :     {
     595           0 :       GEN t = gcoeff(L,i,k);
     596           0 :       p1 = gsub(gmul(gel(B,k+1),gcoeff(L,i,k-1)), gmul(la,t));
     597           0 :       gcoeff(L,i,k) = gdiv(p1, Bk);
     598           0 :       p1 = gadd(gmul(la,gcoeff(L,i,k-1)), gmul(gel(B,k-1),t));
     599           0 :       gcoeff(L,i,k-1) = gdiv(p1, Bk);
     600             :     }
     601             :   }
     602          63 :   else if (!gequal0(la))
     603             :   {
     604          21 :     p1 = gdiv(la2, Bk);
     605          21 :     gel(B,k+1) = gel(B,k) = p1;
     606          21 :     for (i=k+2; i<=lx; i++) gel(B,i) = gdiv(gmul(p1,gel(B,i)),Bk);
     607          21 :     for (i=k+1; i<lx; i++)
     608           0 :       gcoeff(L,i,k-1) = gdiv(gmul(la,gcoeff(L,i,k-1)), Bk);
     609          21 :     for (j=k+1; j<lx-1; j++)
     610           0 :       for (i=j+1; i<lx; i++)
     611           0 :         gcoeff(L,i,j) = gdiv(gmul(p1,gcoeff(L,i,j)), Bk);
     612             :   }
     613             :   else
     614             :   {
     615          42 :     gcoeff(L,k,k-1) = gen_0;
     616          42 :     for (i=k+1; i<lx; i++)
     617             :     {
     618           0 :       gcoeff(L,i,k) = gcoeff(L,i,k-1);
     619           0 :       gcoeff(L,i,k-1) = gen_0;
     620             :     }
     621          42 :     gel(B,k) = gel(B,k-1); fl[k] = 1; fl[k-1] = 0;
     622             :   }
     623          84 :   return 1;
     624             : }
     625             : 
     626             : static void
     627         126 : incrementalGSgen(GEN x, GEN L, GEN B, long k, GEN fl)
     628             : {
     629         126 :   GEN u = NULL; /* gcc -Wall */
     630             :   long i, j;
     631         315 :   for (j = 1; j <= k; j++)
     632         189 :     if (j==k || fl[j])
     633             :     {
     634         189 :       u = gcoeff(x,k,j);
     635         189 :       if (!is_extscalar_t(typ(u))) pari_err_TYPE("incrementalGSgen",u);
     636         252 :       for (i=1; i<j; i++)
     637          63 :         if (fl[i])
     638             :         {
     639          63 :           u = gsub(gmul(gel(B,i+1),u), gmul(gcoeff(L,k,i),gcoeff(L,j,i)));
     640          63 :           u = gdiv(u, gel(B,i));
     641             :         }
     642         189 :       gcoeff(L,k,j) = u;
     643             :     }
     644         126 :   if (gequal0(u)) gel(B,k+1) = gel(B,k);
     645             :   else
     646             :   {
     647          84 :     gel(B,k+1) = gcoeff(L,k,k); gcoeff(L,k,k) = gen_1; fl[k] = 1;
     648             :   }
     649         126 : }
     650             : 
     651             : static GEN
     652         126 : lllgramallgen(GEN x, long flag)
     653             : {
     654         126 :   long lx = lg(x), i, j, k, l, n;
     655             :   pari_sp av;
     656             :   GEN B, L, h, fl;
     657             :   int flc;
     658             : 
     659         126 :   n = lx-1; if (n<=1) return lll_trivial(x,flag);
     660          63 :   if (lgcols(x) != lx) pari_err_DIM("lllgramallgen");
     661             : 
     662          63 :   fl = cgetg(lx, t_VECSMALL);
     663             : 
     664          63 :   av = avma;
     665          63 :   B = scalarcol_shallow(gen_1, lx);
     666          63 :   L = cgetg(lx,t_MAT);
     667         189 :   for (j=1; j<lx; j++) { gel(L,j) = zerocol(n); fl[j] = 0; }
     668             : 
     669          63 :   h = matid(n);
     670         189 :   for (i=1; i<lx; i++)
     671         126 :     incrementalGSgen(x, L, B, i, fl);
     672          63 :   flc = 0;
     673          63 :   for(k=2;;)
     674             :   {
     675         147 :     if (REDgen(k, k-1, h, L, gel(B,k))) flc = 1;
     676         147 :     if (do_SWAPgen(h, L, B, k, fl, &flc)) { if (k > 2) k--; }
     677             :     else
     678             :     {
     679          63 :       for (l=k-2; l>=1; l--)
     680           0 :         if (REDgen(k, l, h, L, gel(B,l+1))) flc = 1;
     681          63 :       if (++k > n) break;
     682             :     }
     683          84 :     if (gc_needed(av,1))
     684             :     {
     685           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"lllgramallgen");
     686           0 :       gerepileall(av,3,&B,&L,&h);
     687             :     }
     688             :   }
     689         105 :   k=1; while (k<lx && !fl[k]) k++;
     690          63 :   return lll_finish(h,k-1,flag);
     691             : }
     692             : 
     693             : static int
     694       22037 : RgM_square(GEN x) { long l = lg(x); return l == 1 || l == lgcols(x); }
     695             : static GEN
     696         126 : lllallgen(GEN x, long flag)
     697             : {
     698         126 :   pari_sp av = avma;
     699         126 :   if (!(flag & LLL_GRAM)) x = gram_matrix(x);
     700          42 :   else if (!RgM_square(x)) pari_err_DIM("qflllgram");
     701         126 :   return gerepilecopy(av, lllgramallgen(x, flag));
     702             : }
     703             : GEN
     704          42 : lllgen(GEN x) { return lllallgen(x, LLL_IM); }
     705             : GEN
     706          42 : lllkerimgen(GEN x) { return lllallgen(x, LLL_ALL); }
     707             : GEN
     708           0 : lllgramgen(GEN x)  { return lllallgen(x, LLL_IM|LLL_GRAM); }
     709             : GEN
     710          42 : lllgramkerimgen(GEN x)  { return lllallgen(x, LLL_ALL|LLL_GRAM); }
     711             : 
     712             : static GEN
     713       28353 : lllall(GEN x, long flag)
     714       28353 : { pari_sp av = avma; return gerepilecopy(av, ZM_lll(x, LLLDFT, flag)); }
     715             : GEN
     716        6861 : lllint(GEN x) { return lllall(x, LLL_IM); }
     717             : GEN
     718          35 : lllkerim(GEN x) { return lllall(x, LLL_ALL); }
     719             : GEN
     720       21415 : lllgramint(GEN x)
     721       21415 : { if (!RgM_square(x)) pari_err_DIM("qflllgram");
     722       21415 :   return lllall(x, LLL_IM | LLL_GRAM); }
     723             : GEN
     724          35 : lllgramkerim(GEN x)
     725          35 : { if (!RgM_square(x)) pari_err_DIM("qflllgram");
     726          35 :   return lllall(x, LLL_ALL | LLL_GRAM); }
     727             : 
     728             : GEN
     729      242406 : lllfp(GEN x, double D, long flag)
     730             : {
     731      242406 :   long n = lg(x)-1;
     732      242406 :   pari_sp av = avma;
     733             :   GEN h;
     734      242406 :   if (n <= 1) return lll_trivial(x,flag);
     735      228062 :   if ((flag & LLL_GRAM) && !RgM_square(x)) pari_err_DIM("qflllgram");
     736      228048 :   h = ZM_lll(RgM_rescale_to_int(x), D, flag);
     737      228020 :   return gerepilecopy(av, h);
     738             : }
     739             : 
     740             : GEN
     741          63 : lllgram(GEN x) { return lllfp(x,LLLDFT,LLL_GRAM|LLL_IM); }
     742             : GEN
     743      226258 : lll(GEN x) { return lllfp(x,LLLDFT,LLL_IM); }
     744             : 
     745             : GEN
     746         301 : qflll0(GEN x, long flag)
     747             : {
     748         301 :   if (typ(x) != t_MAT) pari_err_TYPE("qflll",x);
     749         301 :   switch(flag)
     750             :   {
     751          49 :     case 0: return lll(x);
     752          63 :     case 1: RgM_check_ZM(x,"qflll"); return lllint(x);
     753          49 :     case 2: RgM_check_ZM(x,"qflll"); return lllintpartial(x);
     754           7 :     case 3: RgM_check_ZM(x,"qflll"); return lllall(x, LLL_INPLACE);
     755          49 :     case 4: RgM_check_ZM(x,"qflll"); return lllkerim(x);
     756          42 :     case 5: return lllkerimgen(x);
     757          42 :     case 8: return lllgen(x);
     758           0 :     default: pari_err_FLAG("qflll");
     759             :   }
     760             :   return NULL; /* LCOV_EXCL_LINE */
     761             : }
     762             : 
     763             : GEN
     764         203 : qflllgram0(GEN x, long flag)
     765             : {
     766         203 :   if (typ(x) != t_MAT) pari_err_TYPE("qflllgram",x);
     767         203 :   switch(flag)
     768             :   {
     769          63 :     case 0: return lllgram(x);
     770          49 :     case 1: RgM_check_ZM(x,"qflllgram"); return lllgramint(x);
     771          49 :     case 4: RgM_check_ZM(x,"qflllgram"); return lllgramkerim(x);
     772          42 :     case 5: return lllgramkerimgen(x);
     773           0 :     case 8: return lllgramgen(x);
     774           0 :     default: pari_err_FLAG("qflllgram");
     775             :   }
     776             :   return NULL; /* LCOV_EXCL_LINE */
     777             : }
     778             : 
     779             : /********************************************************************/
     780             : /**                                                                **/
     781             : /**                   INTEGRAL KERNEL (LLL REDUCED)                **/
     782             : /**                                                                **/
     783             : /********************************************************************/
     784             : static GEN
     785          28 : kerint0(GEN M)
     786             : {
     787             :   /* return ZM_lll(M, LLLDFT, LLL_KER); */
     788          28 :   GEN U, H = ZM_hnflll(M,&U,1);
     789          28 :   long d = lg(M)-lg(H);
     790          28 :   if (!d) return cgetg(1, t_MAT);
     791          28 :   return ZM_lll(vecslice(U,1,d), LLLDFT, LLL_INPLACE);
     792             : }
     793             : GEN
     794           0 : kerint(GEN M)
     795             : {
     796           0 :   pari_sp av = avma;
     797           0 :   return gerepilecopy(av, kerint0(M));
     798             : }
     799             : /* OBSOLETE: use kerint */
     800             : GEN
     801          28 : matkerint0(GEN M, long flag)
     802             : {
     803          28 :   pari_sp av = avma;
     804          28 :   if (typ(M) != t_MAT) pari_err_TYPE("matkerint",M);
     805          28 :   M = Q_primpart(M);
     806          28 :   RgM_check_ZM(M, "kerint");
     807          28 :   switch(flag)
     808             :   {
     809          28 :     case 0:
     810          28 :     case 1: return gerepilecopy(av, kerint0(M));
     811           0 :     default: pari_err_FLAG("matkerint");
     812             :   }
     813             :   return NULL; /* LCOV_EXCL_LINE */
     814             : }

Generated by: LCOV version 1.13