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 - bibli1.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.18.1 lcov report (development 30365-beea1ff998) Lines: 1175 1243 94.5 %
Date: 2025-07-01 09:21:48 Functions: 75 81 92.6 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : /********************************************************************/
      16             : /**                                                                **/
      17             : /**                 LLL Algorithm and close friends                **/
      18             : /**                                                                **/
      19             : /********************************************************************/
      20             : #include "pari.h"
      21             : #include "paripriv.h"
      22             : 
      23             : #define DEBUGLEVEL DEBUGLEVEL_qf
      24             : 
      25             : /********************************************************************/
      26             : /**             QR Factorization via Householder matrices          **/
      27             : /********************************************************************/
      28             : static int
      29    24692806 : no_prec_pb(GEN x)
      30             : {
      31    24617690 :   return (typ(x) != t_REAL || realprec(x) > DEFAULTPREC
      32    49310496 :                            || expo(x) < DEFAULTPREC>>1);
      33             : }
      34             : /* Find a Householder transformation which, applied to x[k..#x], zeroes
      35             :  * x[k+1..#x]; fill L = (mu_{i,j}). Return 0 if precision problem [obtained
      36             :  * a 0 vector], 1 otherwise */
      37             : static int
      38    24701673 : FindApplyQ(GEN x, GEN L, GEN B, long k, GEN Q, long prec)
      39             : {
      40    24701673 :   long i, lx = lg(x)-1;
      41    24701673 :   GEN x2, x1, xd = x + (k-1);
      42             : 
      43    24701673 :   x1 = gel(xd,1);
      44    24701673 :   x2 = mpsqr(x1);
      45    24700657 :   if (k < lx)
      46             :   {
      47    19468592 :     long lv = lx - (k-1) + 1;
      48    19468592 :     GEN beta, Nx, v = cgetg(lv, t_VEC);
      49    76903779 :     for (i=2; i<lv; i++) {
      50    57435574 :       x2 = mpadd(x2, mpsqr(gel(xd,i)));
      51    57434772 :       gel(v,i) = gel(xd,i);
      52             :     }
      53    19468205 :     if (!signe(x2)) return 0;
      54    19460040 :     Nx = gsqrt(x2, prec); if (signe(x1) < 0) setsigne(Nx, -1);
      55    19461072 :     gel(v,1) = mpadd(x1, Nx);
      56             : 
      57    19460269 :     if (!signe(x1))
      58      732544 :       beta = gtofp(x2, prec); /* make sure typ(beta) != t_INT */
      59             :     else
      60    18727725 :       beta = mpadd(x2, mpmul(Nx,x1));
      61    19460517 :     gel(Q,k) = mkvec2(invr(beta), v);
      62             : 
      63    19460830 :     togglesign(Nx);
      64    19460559 :     gcoeff(L,k,k) = Nx;
      65             :   }
      66             :   else
      67     5232065 :     gcoeff(L,k,k) = gel(x,k);
      68    24692624 :   gel(B,k) = x2;
      69    70729061 :   for (i=1; i<k; i++) gcoeff(L,k,i) = gel(x,i);
      70    24692624 :   return no_prec_pb(x2);
      71             : }
      72             : 
      73             : /* apply Householder transformation Q = [beta,v] to r with t_INT/t_REAL
      74             :  * coefficients, in place: r -= ((0|v).r * beta) v */
      75             : static void
      76    46046070 : ApplyQ(GEN Q, GEN r)
      77             : {
      78    46046070 :   GEN s, rd, beta = gel(Q,1), v = gel(Q,2);
      79    46046070 :   long i, l = lg(v), lr = lg(r);
      80             : 
      81    46046070 :   rd = r + (lr - l);
      82    46046070 :   s = mpmul(gel(v,1), gel(rd,1));
      83   479552287 :   for (i=2; i<l; i++) s = mpadd(s, mpmul(gel(v,i) ,gel(rd,i)));
      84    46041919 :   s = mpmul(beta, s);
      85   525786011 :   for (i=1; i<l; i++)
      86   479738111 :     if (signe(gel(v,i))) gel(rd,i) = mpsub(gel(rd,i), mpmul(s, gel(v,i)));
      87    46047900 : }
      88             : /* apply Q[1], ..., Q[j-1] to r */
      89             : static GEN
      90    16958613 : ApplyAllQ(GEN Q, GEN r, long j)
      91             : {
      92    16958613 :   pari_sp av = avma;
      93             :   long i;
      94    16958613 :   r = leafcopy(r);
      95    63002278 :   for (i=1; i<j; i++) ApplyQ(gel(Q,i), r);
      96    16956427 :   return gc_GEN(av, r);
      97             : }
      98             : 
      99             : /* same, arbitrary coefficients [20% slower for t_REAL at DEFAULTPREC] */
     100             : static void
     101       22113 : RgC_ApplyQ(GEN Q, GEN r)
     102             : {
     103       22113 :   GEN s, rd, beta = gel(Q,1), v = gel(Q,2);
     104       22113 :   long i, l = lg(v), lr = lg(r);
     105             : 
     106       22113 :   rd = r + (lr - l);
     107       22113 :   s = gmul(gel(v,1), gel(rd,1));
     108      464373 :   for (i=2; i<l; i++) s = gadd(s, gmul(gel(v,i) ,gel(rd,i)));
     109       22113 :   s = gmul(beta, s);
     110      486486 :   for (i=1; i<l; i++)
     111      464373 :     if (signe(gel(v,i))) gel(rd,i) = gsub(gel(rd,i), gmul(s, gel(v,i)));
     112       22113 : }
     113             : static GEN
     114         567 : RgC_ApplyAllQ(GEN Q, GEN r, long j)
     115             : {
     116         567 :   pari_sp av = avma;
     117             :   long i;
     118         567 :   r = leafcopy(r);
     119       22680 :   for (i=1; i<j; i++) RgC_ApplyQ(gel(Q,i), r);
     120         567 :   return gc_GEN(av, r);
     121             : }
     122             : 
     123             : int
     124          21 : RgM_QR_init(GEN x, GEN *pB, GEN *pQ, GEN *pL, long prec)
     125             : {
     126          21 :   x = RgM_gtomp(x, prec);
     127          21 :   return QR_init(x, pB, pQ, pL, prec);
     128             : }
     129             : 
     130             : static void
     131          35 : check_householder(GEN Q)
     132             : {
     133          35 :   long i, l = lg(Q);
     134          35 :   if (typ(Q) != t_VEC) pari_err_TYPE("mathouseholder", Q);
     135         854 :   for (i = 1; i < l; i++)
     136             :   {
     137         826 :     GEN q = gel(Q,i), v;
     138         826 :     if (typ(q) != t_VEC || lg(q) != 3) pari_err_TYPE("mathouseholder", Q);
     139         826 :     v = gel(q,2);
     140         826 :     if (typ(v) != t_VEC || lg(v)+i-2 != l) pari_err_TYPE("mathouseholder", Q);
     141             :   }
     142          28 : }
     143             : 
     144             : GEN
     145          35 : mathouseholder(GEN Q, GEN x)
     146             : {
     147          35 :   long l = lg(Q);
     148          35 :   check_householder(Q);
     149          28 :   switch(typ(x))
     150             :   {
     151          14 :     case t_MAT:
     152          14 :       if (lg(x) == 1) return cgetg(1, t_MAT);
     153          14 :       if (lgcols(x) != l+1) pari_err_TYPE("mathouseholder", x);
     154         574 :       pari_APPLY_same(RgC_ApplyAllQ(Q, gel(x,i), l));
     155           7 :     case t_COL:
     156           7 :       if (lg(x) == l+1) return RgC_ApplyAllQ(Q, x, l);
     157             :   }
     158           7 :   pari_err_TYPE("mathouseholder", x);
     159             :   return NULL; /* LCOV_EXCL_LINE */
     160             : }
     161             : 
     162             : GEN
     163          35 : matqr(GEN x, long flag, long prec)
     164             : {
     165          35 :   pari_sp av = avma;
     166             :   GEN B, Q, L;
     167          35 :   long n = lg(x)-1;
     168          35 :   if (typ(x) != t_MAT) pari_err_TYPE("matqr",x);
     169          35 :   if (!n)
     170             :   {
     171          14 :     if (!flag) retmkvec2(cgetg(1,t_MAT),cgetg(1,t_MAT));
     172           7 :     retmkvec2(cgetg(1,t_VEC),cgetg(1,t_MAT));
     173             :   }
     174          21 :   if (n != nbrows(x)) pari_err_DIM("matqr");
     175          21 :   if (!RgM_QR_init(x, &B,&Q,&L, prec)) pari_err_PREC("matqr");
     176          21 :   if (!flag) Q = shallowtrans(mathouseholder(Q, matid(n)));
     177          21 :   return gc_GEN(av, mkvec2(Q, shallowtrans(L)));
     178             : }
     179             : 
     180             : /* compute B = | x[k] |^2, Q = Householder transforms and L = mu_{i,j} */
     181             : int
     182     7742749 : QR_init(GEN x, GEN *pB, GEN *pQ, GEN *pL, long prec)
     183             : {
     184     7742749 :   long j, k = lg(x)-1;
     185     7742749 :   GEN B = cgetg(k+1, t_VEC), Q = cgetg(k, t_VEC), L = zeromatcopy(k,k);
     186    30205580 :   for (j=1; j<=k; j++)
     187             :   {
     188    24701314 :     GEN r = gel(x,j);
     189    24701314 :     if (j > 1) r = ApplyAllQ(Q, r, j);
     190    24701644 :     if (!FindApplyQ(r, L, B, j, Q, prec)) return 0;
     191             :   }
     192     5504266 :   *pB = B; *pQ = Q; *pL = L; return 1;
     193             : }
     194             : /* x a square t_MAT with t_INT / t_REAL entries and maximal rank. Return
     195             :  * qfgaussred(x~*x) */
     196             : GEN
     197      299166 : gaussred_from_QR(GEN x, long prec)
     198             : {
     199      299166 :   long j, k = lg(x)-1;
     200             :   GEN B, Q, L;
     201      299166 :   if (!QR_init(x, &B,&Q,&L, prec)) return NULL;
     202     1068647 :   for (j=1; j<k; j++)
     203             :   {
     204      769480 :     GEN m = gel(L,j), invNx = invr(gel(m,j));
     205             :     long i;
     206      769463 :     gel(m,j) = gel(B,j);
     207     2974111 :     for (i=j+1; i<=k; i++) gel(m,i) = mpmul(invNx, gel(m,i));
     208             :   }
     209      299167 :   gcoeff(L,j,j) = gel(B,j);
     210      299167 :   return shallowtrans(L);
     211             : }
     212             : GEN
     213       14280 : R_from_QR(GEN x, long prec)
     214             : {
     215             :   GEN B, Q, L;
     216       14280 :   if (!QR_init(x, &B,&Q,&L, prec)) return NULL;
     217       14266 :   return shallowtrans(L);
     218             : }
     219             : 
     220             : /********************************************************************/
     221             : /**             QR Factorization via Gram-Schmidt                  **/
     222             : /********************************************************************/
     223             : 
     224             : /* return Gram-Schmidt orthogonal basis (f) attached to (e), B is the
     225             :  * vector of the (f_i . f_i)*/
     226             : GEN
     227       56723 : RgM_gram_schmidt(GEN e, GEN *ptB)
     228             : {
     229       56723 :   long i,j,lx = lg(e);
     230       56723 :   GEN f = RgM_shallowcopy(e), B, iB;
     231             : 
     232       56723 :   B = cgetg(lx, t_VEC);
     233       56723 :   iB= cgetg(lx, t_VEC);
     234             : 
     235      120533 :   for (i=1;i<lx;i++)
     236             :   {
     237       63809 :     GEN p1 = NULL;
     238       63809 :     pari_sp av = avma;
     239      126546 :     for (j=1; j<i; j++)
     240             :     {
     241       62737 :       GEN mu = gmul(RgV_dotproduct(gel(e,i),gel(f,j)), gel(iB,j));
     242       62737 :       GEN p2 = gmul(mu, gel(f,j));
     243       62737 :       p1 = p1? gadd(p1,p2): p2;
     244             :     }
     245       63809 :     p1 = p1? gc_upto(av, gsub(gel(e,i), p1)): gel(e,i);
     246       63809 :     gel(f,i) = p1;
     247       63809 :     gel(B,i) = RgV_dotsquare(gel(f,i));
     248       63809 :     gel(iB,i) = ginv(gel(B,i));
     249             :   }
     250       56724 :   *ptB = B; return f;
     251             : }
     252             : 
     253             : /* B a Z-basis (which the caller should LLL-reduce for efficiency), t a vector.
     254             :  * Apply Babai's nearest plane algorithm to (B,t) */
     255             : GEN
     256       56723 : RgM_Babai(GEN B, GEN t)
     257             : {
     258       56723 :   GEN C, N, G = RgM_gram_schmidt(B, &N), b = t;
     259       56724 :   long j, n = lg(B)-1;
     260             : 
     261       56724 :   C = cgetg(n+1,t_COL);
     262      120533 :   for (j = n; j > 0; j--)
     263             :   {
     264       63809 :     GEN c = gdiv( RgV_dotproduct(b, gel(G,j)), gel(N,j) );
     265             :     long e;
     266       63809 :     c = grndtoi(c,&e);
     267       63809 :     if (e >= 0) return NULL;
     268       63809 :     if (signe(c)) b = RgC_sub(b, RgC_Rg_mul(gel(B,j), c));
     269       63809 :     gel(C,j) = c;
     270             :   }
     271       56724 :   return C;
     272             : }
     273             : 
     274             : /********************************************************************/
     275             : /**                                                                **/
     276             : /**                          LLL ALGORITHM                         **/
     277             : /**                                                                **/
     278             : /********************************************************************/
     279             : /* Def: a matrix M is said to be -partially reduced- if | m1 +- m2 | >= |m1|
     280             :  * for any two columns m1 != m2, in M.
     281             :  *
     282             :  * Input: an integer matrix mat whose columns are linearly independent. Find
     283             :  * another matrix T such that mat * T is partially reduced.
     284             :  *
     285             :  * Output: mat * T if flag = 0;  T if flag != 0,
     286             :  *
     287             :  * This routine is designed to quickly reduce lattices in which one row
     288             :  * is huge compared to the other rows.  For example, when searching for a
     289             :  * polynomial of degree 3 with root a mod N, the four input vectors might
     290             :  * be the coefficients of
     291             :  *     X^3 - (a^3 mod N), X^2 - (a^2 mod N), X - (a mod N), N.
     292             :  * All four constant coefficients are O(p) and the rest are O(1). By the
     293             :  * pigeon-hole principle, the coefficients of the smallest vector in the
     294             :  * lattice are O(p^(1/4)), hence significant reduction of vector lengths
     295             :  * can be anticipated.
     296             :  *
     297             :  * An improved algorithm would look only at the leading digits of dot*.  It
     298             :  * would use single-precision calculations as much as possible.
     299             :  *
     300             :  * Original code: Peter Montgomery (1994) */
     301             : static GEN
     302          35 : lllintpartialall(GEN m, long flag)
     303             : {
     304          35 :   const long ncol = lg(m)-1;
     305          35 :   const pari_sp av = avma;
     306             :   GEN tm1, tm2, mid;
     307             : 
     308          35 :   if (ncol <= 1) return flag? matid(ncol): gcopy(m);
     309             : 
     310          14 :   tm1 = flag? matid(ncol): NULL;
     311             :   {
     312          14 :     const pari_sp av2 = avma;
     313          14 :     GEN dot11 = ZV_dotsquare(gel(m,1));
     314          14 :     GEN dot22 = ZV_dotsquare(gel(m,2));
     315          14 :     GEN dot12 = ZV_dotproduct(gel(m,1), gel(m,2));
     316          14 :     GEN tm  = matid(2); /* For first two columns only */
     317             : 
     318          14 :     int progress = 0;
     319          14 :     long npass2 = 0;
     320             : 
     321             : /* Row reduce the first two columns of m. Our best result so far is
     322             :  * (first two columns of m)*tm.
     323             :  *
     324             :  * Initially tm = 2 x 2 identity matrix.
     325             :  * Inner products of the reduced matrix are in dot11, dot12, dot22. */
     326          49 :     while (npass2 < 2 || progress)
     327             :     {
     328          42 :       GEN dot12new, q = diviiround(dot12, dot22);
     329             : 
     330          35 :       npass2++; progress = signe(q);
     331          35 :       if (progress)
     332             :       {/* Conceptually replace (v1, v2) by (v1 - q*v2, v2), where v1 and v2
     333             :         * represent the reduced basis for the first two columns of the matrix.
     334             :         * We do this by updating tm and the inner products. */
     335          21 :         togglesign(q);
     336          21 :         dot12new = addii(dot12, mulii(q, dot22));
     337          21 :         dot11 = addii(dot11, mulii(q, addii(dot12, dot12new)));
     338          21 :         dot12 = dot12new;
     339          21 :         ZC_lincomb1_inplace(gel(tm,1), gel(tm,2), q);
     340             :       }
     341             : 
     342             :       /* Interchange the output vectors v1 and v2.  */
     343          35 :       swap(dot11,dot22);
     344          35 :       swap(gel(tm,1), gel(tm,2));
     345             : 
     346             :       /* Occasionally (including final pass) do garbage collection.  */
     347          35 :       if ((npass2 & 0xff) == 0 || !progress)
     348          14 :         (void)gc_all(av2, 4, &dot11,&dot12,&dot22,&tm);
     349             :     } /* while npass2 < 2 || progress */
     350             : 
     351             :     {
     352             :       long i;
     353           7 :       GEN det12 = subii(mulii(dot11, dot22), sqri(dot12));
     354             : 
     355           7 :       mid = cgetg(ncol+1, t_MAT);
     356          21 :       for (i = 1; i <= 2; i++)
     357             :       {
     358          14 :         GEN tmi = gel(tm,i);
     359          14 :         if (tm1)
     360             :         {
     361          14 :           GEN tm1i = gel(tm1,i);
     362          14 :           gel(tm1i,1) = gel(tmi,1);
     363          14 :           gel(tm1i,2) = gel(tmi,2);
     364             :         }
     365          14 :         gel(mid,i) = ZC_lincomb(gel(tmi,1),gel(tmi,2), gel(m,1),gel(m,2));
     366             :       }
     367          42 :       for (i = 3; i <= ncol; i++)
     368             :       {
     369          35 :         GEN c = gel(m,i);
     370          35 :         GEN dot1i = ZV_dotproduct(gel(mid,1), c);
     371          35 :         GEN dot2i = ZV_dotproduct(gel(mid,2), c);
     372             :        /* ( dot11  dot12 ) (q1)   ( dot1i )
     373             :         * ( dot12  dot22 ) (q2) = ( dot2i )
     374             :         *
     375             :         * Round -q1 and -q2 to nearest integer. Then compute
     376             :         *   c - q1*mid[1] - q2*mid[2].
     377             :         * This will be approximately orthogonal to the first two vectors in
     378             :         * the new basis. */
     379          35 :         GEN q1neg = subii(mulii(dot12, dot2i), mulii(dot22, dot1i));
     380          35 :         GEN q2neg = subii(mulii(dot12, dot1i), mulii(dot11, dot2i));
     381             : 
     382          35 :         q1neg = diviiround(q1neg, det12);
     383          35 :         q2neg = diviiround(q2neg, det12);
     384          35 :         if (tm1)
     385             :         {
     386          35 :           gcoeff(tm1,1,i) = addii(mulii(q1neg, gcoeff(tm,1,1)),
     387          35 :                                   mulii(q2neg, gcoeff(tm,1,2)));
     388          35 :           gcoeff(tm1,2,i) = addii(mulii(q1neg, gcoeff(tm,2,1)),
     389          35 :                                   mulii(q2neg, gcoeff(tm,2,2)));
     390             :         }
     391          35 :         gel(mid,i) = ZC_add(c, ZC_lincomb(q1neg,q2neg, gel(mid,1),gel(mid,2)));
     392             :       } /* for i */
     393             :     } /* local block */
     394             :   }
     395           7 :   if (DEBUGLEVEL>6)
     396             :   {
     397           0 :     if (tm1) err_printf("tm1 = %Ps",tm1);
     398           0 :     err_printf("mid = %Ps",mid); /* = m * tm1 */
     399             :   }
     400           7 :   (void)gc_all(av, tm1? 2: 1, &mid, &tm1);
     401             :   {
     402             :    /* For each pair of column vectors v and w in mid * tm2,
     403             :     * try to replace (v, w) by (v, v - q*w) for some q.
     404             :     * We compute all inner products and check them repeatedly. */
     405           7 :     const pari_sp av3 = avma;
     406           7 :     long i, j, npass = 0, e = LONG_MAX;
     407           7 :     GEN dot = cgetg(ncol+1, t_MAT); /* scalar products */
     408             : 
     409           7 :     tm2 = matid(ncol);
     410          56 :     for (i=1; i <= ncol; i++)
     411             :     {
     412          49 :       gel(dot,i) = cgetg(ncol+1,t_COL);
     413         245 :       for (j=1; j <= i; j++)
     414         196 :         gcoeff(dot,j,i) = gcoeff(dot,i,j) = ZV_dotproduct(gel(mid,i),gel(mid,j));
     415             :     }
     416             :     for(;;)
     417          35 :     {
     418          42 :       long reductions = 0, olde = e;
     419         336 :       for (i=1; i <= ncol; i++)
     420             :       {
     421             :         long ijdif;
     422        2058 :         for (ijdif=1; ijdif < ncol; ijdif++)
     423             :         {
     424             :           long d, k1, k2;
     425             :           GEN codi, q;
     426             : 
     427        1764 :           j = i + ijdif; if (j > ncol) j -= ncol;
     428             :           /* let k1, resp. k2,  index of larger, resp. smaller, column */
     429        1764 :           if (cmpii(gcoeff(dot,i,i), gcoeff(dot,j,j)) > 0) { k1 = i; k2 = j; }
     430        1022 :           else                                             { k1 = j; k2 = i; }
     431        1764 :           codi = gcoeff(dot,k2,k2);
     432        1764 :           q = signe(codi)? diviiround(gcoeff(dot,k1,k2), codi): gen_0;
     433        1764 :           if (!signe(q)) continue;
     434             : 
     435             :           /* Try to subtract a multiple of column k2 from column k1.  */
     436         700 :           reductions++; togglesign_safe(&q);
     437         700 :           ZC_lincomb1_inplace(gel(tm2,k1), gel(tm2,k2), q);
     438         700 :           ZC_lincomb1_inplace(gel(dot,k1), gel(dot,k2), q);
     439         700 :           gcoeff(dot,k1,k1) = addii(gcoeff(dot,k1,k1),
     440         700 :                                     mulii(q, gcoeff(dot,k2,k1)));
     441        5600 :           for (d = 1; d <= ncol; d++) gcoeff(dot,k1,d) = gcoeff(dot,d,k1);
     442             :         } /* for ijdif */
     443         294 :         if (gc_needed(av3,2))
     444             :         {
     445           0 :           if(DEBUGMEM>1) pari_warn(warnmem,"lllintpartialall");
     446           0 :           (void)gc_all(av3, 2, &dot,&tm2);
     447             :         }
     448             :       } /* for i */
     449          42 :       if (!reductions) break;
     450          35 :       e = 0;
     451         280 :       for (i = 1; i <= ncol; i++) e += expi( gcoeff(dot,i,i) );
     452          35 :       if (e == olde) break;
     453          35 :       if (DEBUGLEVEL>6)
     454             :       {
     455           0 :         npass++;
     456           0 :         err_printf("npass = %ld, red. last time = %ld, log_2(det) ~ %ld\n\n",
     457             :                     npass, reductions, e);
     458             :       }
     459             :     } /* for(;;)*/
     460             : 
     461             :    /* Sort columns so smallest comes first in m * tm1 * tm2.
     462             :     * Use insertion sort. */
     463          49 :     for (i = 1; i < ncol; i++)
     464             :     {
     465          42 :       long j, s = i;
     466             : 
     467         189 :       for (j = i+1; j <= ncol; j++)
     468         147 :         if (cmpii(gcoeff(dot,s,s),gcoeff(dot,j,j)) > 0) s = j;
     469          42 :       if (i != s)
     470             :       { /* Exchange with proper column; only the diagonal of dot is updated */
     471          28 :         swap(gel(tm2,i), gel(tm2,s));
     472          28 :         swap(gcoeff(dot,i,i), gcoeff(dot,s,s));
     473             :       }
     474             :     }
     475             :   } /* local block */
     476           7 :   return gc_upto(av, ZM_mul(tm1? tm1: mid, tm2));
     477             : }
     478             : 
     479             : GEN
     480          35 : lllintpartial(GEN mat) { return lllintpartialall(mat,1); }
     481             : 
     482             : GEN
     483           0 : lllintpartial_inplace(GEN mat) { return lllintpartialall(mat,0); }
     484             : 
     485             : /********************************************************************/
     486             : /**                                                                **/
     487             : /**                    COPPERSMITH ALGORITHM                       **/
     488             : /**           Finding small roots of univariate equations.         **/
     489             : /**                                                                **/
     490             : /********************************************************************/
     491             : 
     492             : static int
     493         882 : check(double b, double x, double rho, long d, long dim, long delta, long t)
     494             : {
     495         882 :   double cond = delta * (d * (delta+1) - 2*b*dim + rho * (delta-1 + 2*t))
     496         882 :                 + x*dim*(dim - 1);
     497         882 :   if (DEBUGLEVEL >= 4)
     498           0 :     err_printf("delta = %d, t = %d (%.1lf)\n", delta, t, cond);
     499         882 :   return (cond <= 0);
     500             : }
     501             : 
     502             : static void
     503          21 : choose_params(GEN P, GEN N, GEN X, GEN B, long *pdelta, long *pt)
     504             : {
     505          21 :   long d = degpol(P), dim;
     506          21 :   GEN P0 = leading_coeff(P);
     507          21 :   double logN = gtodouble(glog(N, DEFAULTPREC)), x, b, rho;
     508          21 :   x = gtodouble(glog(X, DEFAULTPREC)) / logN;
     509          21 :   b = B? gtodouble(glog(B, DEFAULTPREC)) / logN: 1.;
     510          21 :   if (x * d >= b * b) pari_err_OVERFLOW("zncoppersmith [bound too large]");
     511             :   /* TODO : remove P0 completely */
     512          14 :   rho = is_pm1(P0)? 0: gtodouble(glog(P0, DEFAULTPREC)) / logN;
     513             : 
     514             :   /* Enumerate (delta,t) by increasing lattice dimension */
     515          14 :   for(dim = d + 1;; dim++)
     516         161 :   {
     517             :     long delta, t; /* dim = d*delta + t in the loop */
     518        1043 :     for (delta = 0, t = dim; t >= 0; delta++, t -= d)
     519         882 :       if (check(b,x,rho,d,dim,delta,t)) { *pdelta = delta; *pt = t; return; }
     520             :   }
     521             : }
     522             : 
     523             : static int
     524       14021 : sol_OK(GEN x, GEN N, GEN B)
     525       14021 : { return B? (cmpii(gcdii(x,N),B) >= 0): dvdii(x,N); }
     526             : /* deg(P) > 0, x >= 0. Find all j such that gcd(P(j), N) >= B, |j| <= x */
     527             : static GEN
     528           7 : do_exhaustive(GEN P, GEN N, long x, GEN B)
     529             : {
     530           7 :   GEN Pe, Po, sol = vecsmalltrunc_init(2*x + 2);
     531             :   pari_sp av;
     532             :   long j;
     533           7 :   RgX_even_odd(P, &Pe,&Po); av = avma;
     534           7 :   if (sol_OK(gel(P,2), N,B)) vecsmalltrunc_append(sol, 0);
     535        7007 :   for (j = 1; j <= x; j++, set_avma(av))
     536             :   {
     537        7000 :     GEN j2 = sqru(j), E = FpX_eval(Pe,j2,N), O = FpX_eval(Po,j2,N);
     538        7000 :     if (sol_OK(addmuliu(E,O,j), N,B)) vecsmalltrunc_append(sol, j);
     539        7000 :     if (sol_OK(submuliu(E,O,j), N,B)) vecsmalltrunc_append(sol,-j);
     540             :   }
     541           7 :   vecsmall_sort(sol); return zv_to_ZV(sol);
     542             : }
     543             : 
     544             : /* General Coppersmith, look for a root x0 <= p, p >= B, p | N, |x0| <= X.
     545             :  * B = N coded as NULL */
     546             : GEN
     547          35 : zncoppersmith(GEN P, GEN N, GEN X, GEN B)
     548             : {
     549             :   GEN Q, R, N0, M, sh, short_pol, *Xpowers, sol, nsp, cP, Z;
     550          35 :   long delta, i, j, row, d, l, t, dim, bnd = 10;
     551          35 :   const ulong X_SMALL = 1000;
     552          35 :   pari_sp av = avma;
     553             : 
     554          35 :   if (typ(P) != t_POL || !RgX_is_ZX(P)) pari_err_TYPE("zncoppersmith",P);
     555          28 :   if (typ(N) != t_INT) pari_err_TYPE("zncoppersmith",N);
     556          28 :   if (typ(X) != t_INT) {
     557           7 :     X = gfloor(X);
     558           7 :     if (typ(X) != t_INT) pari_err_TYPE("zncoppersmith",X);
     559             :   }
     560          28 :   if (signe(X) < 0) pari_err_DOMAIN("zncoppersmith", "X", "<", gen_0, X);
     561          28 :   P = FpX_red(P, N); d = degpol(P);
     562          28 :   if (d == 0) retgc_const(av, cgetg(1, t_VEC));
     563          28 :   if (d < 0) pari_err_ROOTS0("zncoppersmith");
     564          28 :   if (B && typ(B) != t_INT) B = gceil(B);
     565          28 :   if (abscmpiu(X, X_SMALL) <= 0)
     566           7 :     return gc_upto(av, do_exhaustive(P, N, itos(X), B));
     567             : 
     568          21 :   if (B && equalii(B,N)) B = NULL;
     569          21 :   if (B) bnd = 1; /* bnd-hack is only for the case B = N */
     570          21 :   cP = gel(P,d+2);
     571          21 :   if (!gequal1(cP))
     572             :   {
     573             :     GEN r, z;
     574          14 :     gel(P,d+2) = cP = bezout(cP, N, &z, &r);
     575          35 :     for (j = 0; j < d; j++) gel(P,j+2) = Fp_mul(gel(P,j+2), z, N);
     576          14 :     if (!is_pm1(cP))
     577             :     {
     578           7 :       P = Q_primitive_part(P, &cP);
     579           7 :       if (cP) { N = diviiexact(N,cP); B = gceil(gdiv(B, cP)); }
     580             :     }
     581             :   }
     582          21 :   if (DEBUGLEVEL >= 2) err_printf("Modified P: %Ps\n", P);
     583             : 
     584          21 :   choose_params(P, N, X, B, &delta, &t);
     585          14 :   if (DEBUGLEVEL >= 2)
     586           0 :     err_printf("Init: trying delta = %d, t = %d\n", delta, t);
     587             :   for(;;)
     588             :   {
     589          14 :     dim = d * delta + t;
     590             :     /* TODO: In case of failure do not recompute the full vector */
     591          14 :     Xpowers = (GEN*)new_chunk(dim + 1);
     592          14 :     Xpowers[0] = gen_1;
     593         217 :     for (j = 1; j <= dim; j++) Xpowers[j] = mulii(Xpowers[j-1], X);
     594             : 
     595             :     /* TODO: in case of failure, use the part of the matrix already computed */
     596          14 :     M = zeromatcopy(dim,dim);
     597             : 
     598             :     /* Rows of M correspond to the polynomials
     599             :      * N^delta, N^delta Xi, ... N^delta (Xi)^d-1,
     600             :      * N^(delta-1)P(Xi), N^(delta-1)XiP(Xi), ... N^(delta-1)P(Xi)(Xi)^d-1,
     601             :      * ...
     602             :      * P(Xi)^delta, XiP(Xi)^delta, ..., P(Xi)^delta(Xi)^t-1 */
     603          42 :     for (j = 1; j <= d;   j++) gcoeff(M, j, j) = gel(Xpowers,j-1);
     604             : 
     605             :     /* P-part */
     606          14 :     if (delta) row = d + 1; else row = 0;
     607             : 
     608          14 :     Q = P;
     609          70 :     for (i = 1; i < delta; i++)
     610             :     {
     611         182 :       for (j = 0; j < d; j++,row++)
     612        1239 :         for (l = j + 1; l <= row; l++)
     613        1113 :           gcoeff(M, l, row) = mulii(Xpowers[l-1], gel(Q,l-j+1));
     614          56 :       Q = ZX_mul(Q, P);
     615             :     }
     616          63 :     for (j = 0; j < t; row++, j++)
     617         490 :       for (l = j + 1; l <= row; l++)
     618         441 :         gcoeff(M, l, row) = mulii(Xpowers[l-1], gel(Q,l-j+1));
     619             : 
     620             :     /* N-part */
     621          14 :     row = dim - t; N0 = N;
     622          84 :     while (row >= 1)
     623             :     {
     624         224 :       for (j = 0; j < d; j++,row--)
     625        1421 :         for (l = 1; l <= row; l++)
     626        1267 :           gcoeff(M, l, row) = mulii(gmael(M, row, l), N0);
     627          70 :       if (row >= 1) N0 = mulii(N0, N);
     628             :     }
     629             :     /* Z is the upper bound for the L^1 norm of the polynomial,
     630             :        ie. N^delta if B = N, B^delta otherwise */
     631          14 :     if (B) Z = powiu(B, delta); else Z = N0;
     632             : 
     633          14 :     if (DEBUGLEVEL >= 2)
     634             :     {
     635           0 :       if (DEBUGLEVEL >= 6) err_printf("Matrix to be reduced:\n%Ps\n", M);
     636           0 :       err_printf("Entering LLL\nbitsize bound: %ld\n", expi(Z));
     637           0 :       err_printf("expected shvector bitsize: %ld\n", expi(ZM_det_triangular(M))/dim);
     638             :     }
     639             : 
     640          14 :     sh = ZM_lll(M, 0.75, LLL_INPLACE);
     641             :     /* Take the first vector if it is non constant */
     642          14 :     short_pol = gel(sh,1);
     643          14 :     if (ZV_isscalar(short_pol)) short_pol = gel(sh, 2);
     644             : 
     645          14 :     nsp = gen_0;
     646         217 :     for (j = 1; j <= dim; j++) nsp = addii(nsp, absi_shallow(gel(short_pol,j)));
     647             : 
     648          14 :     if (DEBUGLEVEL >= 2)
     649             :     {
     650           0 :       err_printf("Candidate: %Ps\n", short_pol);
     651           0 :       err_printf("bitsize Norm: %ld\n", expi(nsp));
     652           0 :       err_printf("bitsize bound: %ld\n", expi(mului(bnd, Z)));
     653             :     }
     654          14 :     if (cmpii(nsp, mului(bnd, Z)) < 0) break; /* SUCCESS */
     655             : 
     656             :     /* Failed with the precomputed or supplied value */
     657           0 :     if (++t == d) { delta++; t = 1; }
     658           0 :     if (DEBUGLEVEL >= 2)
     659           0 :       err_printf("Increasing dim, delta = %d t = %d\n", delta, t);
     660             :   }
     661          14 :   bnd = itos(divii(nsp, Z)) + 1;
     662             : 
     663          14 :   while (!signe(gel(short_pol,dim))) dim--;
     664             : 
     665          14 :   R = cgetg(dim + 2, t_POL); R[1] = P[1];
     666         217 :   for (j = 1; j <= dim; j++)
     667         203 :     gel(R,j+1) = diviiexact(gel(short_pol,j), Xpowers[j-1]);
     668          14 :   gel(R,2) = subii(gel(R,2), mului(bnd - 1, N0));
     669             : 
     670          14 :   sol = cgetg(1, t_VEC);
     671          84 :   for (i = -bnd + 1; i < bnd; i++)
     672             :   {
     673          70 :     GEN r = nfrootsQ(R);
     674          70 :     if (DEBUGLEVEL >= 2) err_printf("Roots: %Ps\n", r);
     675          91 :     for (j = 1; j < lg(r); j++)
     676             :     {
     677          21 :       GEN z = gel(r,j);
     678          21 :       if (typ(z) == t_INT && sol_OK(FpX_eval(P,z,N), N,B))
     679          14 :         sol = shallowconcat(sol, z);
     680             :     }
     681          70 :     if (i < bnd) gel(R,2) = addii(gel(R,2), Z);
     682             :   }
     683          14 :   return gc_upto(av, ZV_sort_uniq(sol));
     684             : }
     685             : 
     686             : /********************************************************************/
     687             : /**                                                                **/
     688             : /**                   LINEAR & ALGEBRAIC DEPENDENCE                **/
     689             : /**                                                                **/
     690             : /********************************************************************/
     691             : 
     692             : static int
     693        8123 : real_indep(GEN re, GEN im, long bit)
     694             : {
     695        8123 :   GEN d = gsub(gmul(gel(re,1),gel(im,2)), gmul(gel(re,2),gel(im,1)));
     696        8123 :   return (!gequal0(d) && gexpo(d) > - bit);
     697             : }
     698             : 
     699             : GEN
     700       15302 : lindepfull_bit(GEN x, long bit)
     701             : {
     702       15302 :   long lx = lg(x), ly, i, j;
     703             :   GEN re, im, M;
     704             : 
     705       15302 :   if (! is_vec_t(typ(x))) pari_err_TYPE("lindep2",x);
     706       15302 :   if (lx <= 2)
     707             :   {
     708          21 :     if (lx == 2 && gequal0(x)) return matid(1);
     709          14 :     return NULL;
     710             :   }
     711       15281 :   re = real_i(x);
     712       15281 :   im = imag_i(x);
     713             :   /* independent over R ? */
     714       15281 :   if (lx == 3 && real_indep(re,im,bit)) return NULL;
     715       15267 :   if (gequal0(im)) im = NULL;
     716       15267 :   ly = im? lx+2: lx+1;
     717       15267 :   M = cgetg(lx,t_MAT);
     718       60701 :   for (i=1; i<lx; i++)
     719             :   {
     720       45434 :     GEN c = cgetg(ly,t_COL); gel(M,i) = c;
     721      209394 :     for (j=1; j<lx; j++) gel(c,j) = gen_0;
     722       45434 :     gel(c,i) = gen_1;
     723       45434 :     gel(c,lx)           = gtrunc2n(gel(re,i), bit);
     724       45434 :     if (im) gel(c,lx+1) = gtrunc2n(gel(im,i), bit);
     725             :   }
     726       15267 :   return ZM_lll(M, 0.99, LLL_INPLACE);
     727             : }
     728             : GEN
     729        3311 : lindep_bit(GEN x, long bit)
     730             : {
     731        3311 :   pari_sp av = avma;
     732        3311 :   GEN v, M = lindepfull_bit(x,bit);
     733        3311 :   if (!M) retgc_const(av, cgetg(1, t_COL));
     734        3283 :   v = gel(M,1); setlg(v, lg(M));
     735        3283 :   return gc_GEN(av, v);
     736             : }
     737             : /* deprecated */
     738             : GEN
     739         112 : lindep2(GEN x, long dig)
     740             : {
     741             :   long bit;
     742         112 :   if (dig < 0) pari_err_DOMAIN("lindep2", "accuracy", "<", gen_0, stoi(dig));
     743         112 :   if (dig) bit = (long) (dig/LOG10_2);
     744             :   else
     745             :   {
     746          98 :     bit = gprecision(x);
     747          98 :     if (!bit)
     748             :     {
     749          35 :       x = Q_primpart(x); /* left on stack */
     750          35 :       bit = 32 + gexpo(x);
     751             :     }
     752             :     else
     753          63 :       bit = (long)prec2nbits_mul(bit, 0.8);
     754             :   }
     755         112 :   return lindep_bit(x, bit);
     756             : }
     757             : 
     758             : /* x is a vector of elts of a p-adic field */
     759             : GEN
     760          28 : lindep_padic(GEN x)
     761             : {
     762          28 :   long i, j, prec = LONG_MAX, nx = lg(x)-1, v;
     763          28 :   pari_sp av = avma;
     764          28 :   GEN p = NULL, pn, m, a;
     765             : 
     766          28 :   if (nx < 2) return cgetg(1,t_COL);
     767         147 :   for (i=1; i<=nx; i++)
     768             :   {
     769         119 :     GEN c = gel(x,i), q;
     770         119 :     if (typ(c) != t_PADIC) continue;
     771             : 
     772          91 :     j = precp(c); if (j < prec) prec = j;
     773          91 :     q = padic_p(c);
     774          91 :     if (!p) p = q; else if (!equalii(p, q)) pari_err_MODULUS("lindep_padic", p, q);
     775             :   }
     776          28 :   if (!p) pari_err_TYPE("lindep_padic [not a p-adic vector]",x);
     777          28 :   v = gvaluation(x,p); pn = powiu(p,prec);
     778          28 :   if (v) x = gmul(x, powis(p, -v));
     779          28 :   x = RgV_to_FpV(x, pn);
     780             : 
     781          28 :   a = negi(gel(x,1));
     782          28 :   m = cgetg(nx,t_MAT);
     783         119 :   for (i=1; i<nx; i++)
     784             :   {
     785          91 :     GEN c = zerocol(nx);
     786          91 :     gel(c,1+i) = a;
     787          91 :     gel(c,1) = gel(x,i+1);
     788          91 :     gel(m,i) = c;
     789             :   }
     790          28 :   m = ZM_lll(ZM_hnfmodid(m, pn), 0.99, LLL_INPLACE);
     791          28 :   return gc_GEN(av, gel(m,1));
     792             : }
     793             : /* x is a vector of t_POL/t_SER */
     794             : GEN
     795          77 : lindep_Xadic(GEN x)
     796             : {
     797          77 :   long i, prec = LONG_MAX, deg = 0, lx = lg(x), vx, v;
     798          77 :   pari_sp av = avma;
     799             :   GEN m;
     800             : 
     801          77 :   if (lx == 1) return cgetg(1,t_COL);
     802          77 :   vx = gvar(x);
     803          77 :   if (gequal0(x)) return col_ei(lx-1,1);
     804          70 :   v = gvaluation(x, pol_x(vx));
     805          70 :   if (!v)         x = shallowcopy(x);
     806           0 :   else if (v > 0) x = gdiv(x, pol_xn(v, vx));
     807           0 :   else            x = gmul(x, pol_xn(-v, vx));
     808             :   /* all t_SER have valuation >= 0 */
     809         735 :   for (i=1; i<lx; i++)
     810             :   {
     811         665 :     GEN c = gel(x,i);
     812         665 :     if (gvar(c) != vx) { gel(x,i) = scalarpol_shallow(c, vx); continue; }
     813         658 :     switch(typ(c))
     814             :     {
     815         231 :       case t_POL: deg = maxss(deg, degpol(c)); break;
     816           0 :       case t_RFRAC: pari_err_TYPE("lindep_Xadic", c);
     817         427 :       case t_SER:
     818         427 :         prec = minss(prec, valser(c)+lg(c)-2);
     819         427 :         gel(x,i) = ser2rfrac_i(c);
     820             :     }
     821             :   }
     822          70 :   if (prec == LONG_MAX) prec = deg+1;
     823          70 :   m = RgXV_to_RgM(x, prec);
     824          70 :   return gc_upto(av, deplin(m));
     825             : }
     826             : static GEN
     827          35 : vec_lindep(GEN x)
     828             : {
     829          35 :   pari_sp av = avma;
     830          35 :   long i, l = lg(x); /* > 1 */
     831          35 :   long t = typ(gel(x,1)), h = lg(gel(x,1));
     832          35 :   GEN m = cgetg(l, t_MAT);
     833         126 :   for (i = 1; i < l; i++)
     834             :   {
     835          98 :     GEN y = gel(x,i);
     836          98 :     if (lg(y) != h || typ(y) != t) pari_err_TYPE("lindep",x);
     837          91 :     if (t != t_COL) y = shallowtrans(y); /* Sigh */
     838          91 :     gel(m,i) = y;
     839             :   }
     840          28 :   return gc_upto(av, deplin(m));
     841             : }
     842             : 
     843             : GEN
     844           0 : lindep(GEN x) { return lindep2(x, 0); }
     845             : 
     846             : GEN
     847         434 : lindep0(GEN x,long bit)
     848             : {
     849         434 :   long i, tx = typ(x);
     850         434 :   if (tx == t_MAT) return deplin(x);
     851         147 :   if (! is_vec_t(tx)) pari_err_TYPE("lindep",x);
     852         441 :   for (i = 1; i < lg(x); i++)
     853         357 :     switch(typ(gel(x,i)))
     854             :     {
     855           7 :       case t_PADIC: return lindep_padic(x);
     856          21 :       case t_POL:
     857             :       case t_RFRAC:
     858          21 :       case t_SER: return lindep_Xadic(x);
     859          35 :       case t_VEC:
     860          35 :       case t_COL: return vec_lindep(x);
     861             :     }
     862          84 :   return lindep2(x, bit);
     863             : }
     864             : 
     865             : GEN
     866          77 : algdep0(GEN x, long n, long bit)
     867             : {
     868          77 :   long tx = typ(x), i;
     869             :   pari_sp av;
     870             :   GEN y;
     871             : 
     872          77 :   if (! is_scalar_t(tx)) pari_err_TYPE("algdep0",x);
     873          77 :   if (tx == t_POLMOD)
     874             :   {
     875          14 :     av = avma; y = minpoly(x, 0);
     876          14 :     return (degpol(y) > n)? gc_const(av, gen_1): y;
     877             :   }
     878          63 :   if (gequal0(x)) return pol_x(0);
     879          63 :   if (n <= 0)
     880             :   {
     881          14 :     if (!n) return gen_1;
     882           7 :     pari_err_DOMAIN("algdep", "degree", "<", gen_0, stoi(n));
     883             :   }
     884             : 
     885          49 :   av = avma; y = cgetg(n+2,t_COL);
     886          49 :   gel(y,1) = gen_1;
     887          49 :   gel(y,2) = x; /* n >= 1 */
     888         210 :   for (i=3; i<=n+1; i++) gel(y,i) = gmul(gel(y,i-1),x);
     889          49 :   if (typ(x) == t_PADIC)
     890          21 :     y = lindep_padic(y);
     891             :   else
     892          28 :     y = lindep2(y, bit);
     893          49 :   if (lg(y) == 1) pari_err(e_DOMAIN,"algdep", "degree(x)",">", stoi(n), x);
     894          49 :   y = RgV_to_RgX(y, 0);
     895          49 :   if (signe(leading_coeff(y)) > 0) return gc_GEN(av, y);
     896          14 :   return gc_upto(av, ZX_neg(y));
     897             : }
     898             : 
     899             : GEN
     900           0 : algdep(GEN x, long n)
     901             : {
     902           0 :   return algdep0(x,n,0);
     903             : }
     904             : 
     905             : static GEN
     906          56 : sertomat(GEN S, long p, long r, long vy)
     907             : {
     908             :   long n, m;
     909          56 :   GEN v = cgetg(r*p+1, t_VEC); /* v[r*n+m+1] = s^n * y^m */
     910             :   /* n = 0 */
     911         245 :   for (m = 0; m < r; m++) gel(v, m+1) = pol_xn(m, vy);
     912         175 :   for(n=1; n < p; n++)
     913         546 :     for (m = 0; m < r; m++)
     914             :     {
     915         427 :       GEN c = gel(S,n);
     916         427 :       if (m)
     917             :       {
     918         308 :         c = shallowcopy(c);
     919         308 :         setvalser(c, valser(c) + m);
     920             :       }
     921         427 :       gel(v, r*n + m + 1) = c;
     922             :     }
     923          56 :   return v;
     924             : }
     925             : 
     926             : GEN
     927          42 : seralgdep(GEN s, long p, long r)
     928             : {
     929          42 :   pari_sp av = avma;
     930             :   long vy, i, n, prec;
     931             :   GEN S, v, D;
     932             : 
     933          42 :   if (typ(s) != t_SER) pari_err_TYPE("seralgdep",s);
     934          42 :   if (p <= 0) pari_err_DOMAIN("seralgdep", "p", "<=", gen_0, stoi(p));
     935          42 :   if (r < 0) pari_err_DOMAIN("seralgdep", "r", "<", gen_0, stoi(r));
     936          42 :   if (is_bigint(addiu(muluu(p, r), 1))) pari_err_OVERFLOW("seralgdep");
     937          42 :   vy = varn(s);
     938          42 :   if (!vy) pari_err_PRIORITY("seralgdep", s, ">", 0);
     939          42 :   r++; p++;
     940          42 :   prec = valser(s) + lg(s)-2;
     941          42 :   if (r > prec) r = prec;
     942          42 :   S = cgetg(p+1, t_VEC); gel(S, 1) = s;
     943         119 :   for (i = 2; i <= p; i++) gel(S,i) = gmul(gel(S,i-1), s);
     944          42 :   v = sertomat(S, p, r, vy);
     945          42 :   D = lindep_Xadic(v);
     946          42 :   if (lg(D) == 1) { set_avma(av); return gen_0; }
     947          35 :   v = cgetg(p+1, t_VEC);
     948         133 :   for (n = 0; n < p; n++)
     949          98 :     gel(v, n+1) = RgV_to_RgX(vecslice(D, r*n+1, r*n+r), vy);
     950          35 :   return gc_GEN(av, RgV_to_RgX(v, 0));
     951             : }
     952             : 
     953             : GEN
     954          14 : serdiffdep(GEN s, long p, long r)
     955             : {
     956          14 :   pari_sp av = avma;
     957             :   long vy, i, n, prec;
     958             :   GEN P, S, v, D;
     959             : 
     960          14 :   if (typ(s) != t_SER) pari_err_TYPE("serdiffdep",s);
     961          14 :   if (p <= 0) pari_err_DOMAIN("serdiffdep", "p", "<=", gen_0, stoi(p));
     962          14 :   if (r < 0) pari_err_DOMAIN("serdiffdep", "r", "<", gen_0, stoi(r));
     963          14 :   if (is_bigint(addiu(muluu(p, r), 1))) pari_err_OVERFLOW("serdiffdep");
     964          14 :   vy = varn(s);
     965          14 :   if (!vy) pari_err_PRIORITY("serdiffdep", s, ">", 0);
     966          14 :   r++; p++;
     967          14 :   prec = valser(s) + lg(s)-2;
     968          14 :   if (r > prec) r = prec;
     969          14 :   S = cgetg(p+1, t_VEC); gel(S, 1) = s;
     970          56 :   for (i = 2; i <= p; i++) gel(S,i) = derivser(gel(S,i-1));
     971          14 :   v = sertomat(S, p, r, vy);
     972          14 :   D = lindep_Xadic(v);
     973          14 :   if (lg(D) == 1) { set_avma(av); return gen_0; }
     974          14 :   P = RgV_to_RgX(vecslice(D, 1, r), vy);
     975          14 :   v = cgetg(p, t_VEC);
     976          56 :   for (n = 1; n < p; n++)
     977          42 :     gel(v, n) = RgV_to_RgX(vecslice(D, r*n+1, r*n+r), vy);
     978          14 :   return gc_GEN(av, mkvec2(RgV_to_RgX(v, 0), gneg(P)));
     979             : }
     980             : 
     981             : /* FIXME: could precompute ZM_lll attached to V[2..] */
     982             : static GEN
     983       11991 : lindepcx(GEN V, long bit)
     984             : {
     985       11991 :   GEN Vr = real_i(V), Vi = imag_i(V);
     986       11991 :   long d = gexpo(Vr) - gexpo(Vi);
     987       11991 :   if (d < -bit) V = Vi;
     988       11991 :   else if (d > bit) V = Vr;
     989       11991 :   return lindepfull_bit(V, bit);
     990             : }
     991             : /* c floating point t_REAL or t_COMPLEX, T ZX, recognize in Q[x]/(T).
     992             :  * V helper vector (containing complex roots of T), MODIFIED */
     993             : static GEN
     994       11991 : cx_bestapprnf(GEN c, GEN T, GEN V, long bit)
     995             : {
     996       11991 :   GEN M, a, v = NULL;
     997             :   long i, l;
     998       11991 :   gel(V,1) = gneg(c); M = lindepcx(V, bit);
     999       11991 :   if (!M) pari_err(e_MISC, "cannot rationalize coeff in bestapprnf");
    1000       11991 :   l = lg(M); a = NULL;
    1001       11991 :   for (i = 1; i < l; i ++) { v = gel(M,i); a = gel(v,1); if (signe(a)) break; }
    1002       11991 :   v = RgC_Rg_div(vecslice(v, 2, lg(M)-1), a);
    1003       11991 :   if (!T) return gel(v,1);
    1004        4830 :   v = RgV_to_RgX(v, varn(T)); l = lg(v);
    1005        4830 :   if (l == 2) return gen_0;
    1006        4165 :   if (l == 3) return gel(v,2);
    1007        3668 :   return mkpolmod(v, T);
    1008             : }
    1009             : static GEN
    1010       15624 : bestapprnf_i(GEN x, GEN T, GEN V, long bit)
    1011             : {
    1012       15624 :   long i, l, tx = typ(x);
    1013             :   GEN z;
    1014       15624 :   switch (tx)
    1015             :   {
    1016        1505 :     case t_INT: case t_FRAC: return x;
    1017       11991 :     case t_REAL: case t_COMPLEX: return cx_bestapprnf(x, T, V, bit);
    1018           0 :     case t_POLMOD: if (RgX_equal(gel(x,1),T)) return x;
    1019           0 :                    break;
    1020        2128 :     case t_POL: case t_SER: case t_VEC: case t_COL: case t_MAT:
    1021        2128 :       l = lg(x); z = cgetg(l, tx);
    1022        3654 :       for (i = 1; i < lontyp[tx]; i++) z[i] = x[i];
    1023       15589 :       for (; i < l; i++) gel(z,i) = bestapprnf_i(gel(x,i), T, V, bit);
    1024        2128 :       return z;
    1025             :   }
    1026           0 :   pari_err_TYPE("mfcxtoQ", x);
    1027             :   return NULL;/*LCOV_EXCL_LINE*/
    1028             : }
    1029             : 
    1030             : GEN
    1031        2163 : bestapprnf(GEN x, GEN T, GEN roT, long prec)
    1032             : {
    1033        2163 :   pari_sp av = avma;
    1034        2163 :   long tx = typ(x), dT = 1, bit;
    1035             :   GEN V;
    1036             : 
    1037        2163 :   if (T)
    1038             :   {
    1039        1610 :     if (typ(T) != t_POL) T = nf_get_pol(checknf(T));
    1040        1610 :     else if (!RgX_is_ZX(T)) pari_err_TYPE("bestapprnf", T);
    1041        1610 :     dT = degpol(T);
    1042             :   }
    1043        2163 :   if (is_rational_t(tx)) return gcopy(x);
    1044        2163 :   if (tx == t_POLMOD)
    1045             :   {
    1046           0 :     if (!T || !RgX_equal(T, gel(x,1))) pari_err_TYPE("bestapprnf",x);
    1047           0 :     return gcopy(x);
    1048             :   }
    1049             : 
    1050        2163 :   if (roT)
    1051             :   {
    1052         644 :     long l = gprecision(roT);
    1053         644 :     switch(typ(roT))
    1054             :     {
    1055         644 :       case t_INT: case t_FRAC: case t_REAL: case t_COMPLEX: break;
    1056           0 :       default: pari_err_TYPE("bestapprnf", roT);
    1057             :     }
    1058         644 :     if (prec < l) prec = l;
    1059             :   }
    1060        1519 :   else if (!T)
    1061         553 :     roT = gen_1;
    1062             :   else
    1063             :   {
    1064         966 :     long n = poliscyclo(T); /* cyclotomic is an important special case */
    1065         966 :     roT = n? rootsof1u_cx(n,prec): gel(QX_complex_roots(T,prec), 1);
    1066             :   }
    1067        2163 :   V = vec_prepend(gpowers(roT, dT-1), NULL);
    1068        2163 :   bit = prec2nbits_mul(prec, 0.8);
    1069        2163 :   return gc_GEN(av, bestapprnf_i(x, T, V, bit));
    1070             : }
    1071             : 
    1072             : /********************************************************************/
    1073             : /**                                                                **/
    1074             : /**                              MINIM                             **/
    1075             : /**                                                                **/
    1076             : /********************************************************************/
    1077             : void
    1078      123292 : minim_alloc(long n, double ***q, GEN *x, double **y,  double **z, double **v)
    1079             : {
    1080      123292 :   long i, s = n * sizeof(double);
    1081             : 
    1082      123292 :   *x = cgetg(n, t_VECSMALL);
    1083      123292 :   *q = (double**) new_chunk(n);
    1084      123292 :   *y = (double*) stack_malloc_align(s, sizeof(double));
    1085      123293 :   *z = (double*) stack_malloc_align(s, sizeof(double));
    1086      123293 :   *v = (double*) stack_malloc_align(s, sizeof(double));
    1087      532950 :   for (i=1; i<n; i++) (*q)[i] = (double*) stack_malloc_align(s, sizeof(double));
    1088      123291 : }
    1089             : 
    1090             : static void
    1091          70 : cvp_alloc(long n, double **t, double **tpre)
    1092             : {
    1093          70 :   long s = n * sizeof(double);
    1094          70 :   *t = (double*) stack_malloc_align(s, sizeof(double));
    1095          70 :   *tpre = (double*) stack_malloc_align(s, sizeof(double));
    1096          70 : }
    1097             : 
    1098             : static GEN
    1099        5502 : ZC_canon(GEN V)
    1100             : {
    1101        5502 :   long l = lg(V), j, s;
    1102       11242 :   for (j = 1; j < l; j++)
    1103       11242 :     if ((s = signe(gel(V,j)))) return s < 0? ZC_neg(V): V;
    1104           0 :   return V;
    1105             : }
    1106             : static GEN
    1107        5502 : ZM_zc_mul_canon(GEN u, GEN x) { return ZC_canon(ZM_zc_mul(u,x)); }
    1108             : static GEN
    1109      240366 : ZM_zc_mul_canon_zm(GEN u, GEN x)
    1110             : {
    1111      240366 :   pari_sp av = avma;
    1112      240366 :   GEN y = ZV_to_zv(ZM_zc_mul(u,x));
    1113      240366 :   zv_canon_inplace(y); return gc_upto(av, y);
    1114             : }
    1115             : 
    1116             : struct qfvec
    1117             : {
    1118             :   GEN a, r, u;
    1119             : };
    1120             : 
    1121             : static void
    1122           0 : err_minim(GEN a)
    1123             : {
    1124           0 :   pari_err_DOMAIN("minim0","form","is not",
    1125             :                   strtoGENstr("positive definite"),a);
    1126           0 : }
    1127             : 
    1128             : static GEN
    1129         902 : minim_lll(GEN a, GEN *u)
    1130             : {
    1131         902 :   *u = lllgramint(a);
    1132         902 :   if (lg(*u) != lg(a)) err_minim(a);
    1133         902 :   return qf_ZM_apply(a,*u);
    1134             : }
    1135             : 
    1136             : static void
    1137         902 : forqfvec_init_dolll(struct qfvec *qv, GEN *pa, long dolll)
    1138             : {
    1139         902 :   GEN r, u, a = *pa;
    1140         902 :   if (!dolll) u = NULL;
    1141             :   else
    1142             :   {
    1143         860 :     if (typ(a) != t_MAT || !RgM_is_ZM(a)) pari_err_TYPE("qfminim",a);
    1144         860 :     a = *pa = minim_lll(a, &u);
    1145             :   }
    1146         902 :   qv->a = RgM_gtofp(a, DEFAULTPREC);
    1147         902 :   r = qfgaussred_positive(qv->a);
    1148         902 :   if (!r)
    1149             :   {
    1150           0 :     r = qfgaussred_positive(a); /* exact computation */
    1151           0 :     if (!r) err_minim(a);
    1152           0 :     r = RgM_gtofp(r, DEFAULTPREC);
    1153             :   }
    1154         902 :   qv->r = r;
    1155         902 :   qv->u = u;
    1156         902 : }
    1157             : 
    1158             : static void
    1159          42 : forqfvec_init(struct qfvec *qv, GEN a)
    1160          42 : { forqfvec_init_dolll(qv, &a, 1); }
    1161             : 
    1162             : static void
    1163          42 : forqfvec_i(void *E, long (*fun)(void *, GEN, GEN, double), struct qfvec *qv, GEN BORNE)
    1164             : {
    1165          42 :   GEN x, a = qv->a, r = qv->r, u = qv->u;
    1166          42 :   long n = lg(a)-1, i, j, k;
    1167             :   double p,BOUND,*v,*y,*z,**q;
    1168          42 :   const double eps = 1e-10;
    1169          42 :   if (!BORNE) BORNE = gen_0;
    1170             :   else
    1171             :   {
    1172          28 :     BORNE = gfloor(BORNE);
    1173          28 :     if (typ(BORNE) != t_INT) pari_err_TYPE("minim0",BORNE);
    1174          35 :     if (signe(BORNE) <= 0) return;
    1175             :   }
    1176          35 :   if (n == 0) return;
    1177          28 :   minim_alloc(n+1, &q, &x, &y, &z, &v);
    1178          98 :   for (j=1; j<=n; j++)
    1179             :   {
    1180          70 :     v[j] = rtodbl(gcoeff(r,j,j));
    1181         133 :     for (i=1; i<j; i++) q[i][j] = rtodbl(gcoeff(r,i,j));
    1182             :   }
    1183             : 
    1184          28 :   if (gequal0(BORNE))
    1185             :   {
    1186             :     double c;
    1187          14 :     p = rtodbl(gcoeff(a,1,1));
    1188          42 :     for (i=2; i<=n; i++) { c = rtodbl(gcoeff(a,i,i)); if (c < p) p = c; }
    1189          14 :     BORNE = roundr(dbltor(p));
    1190             :   }
    1191             :   else
    1192          14 :     p = gtodouble(BORNE);
    1193          28 :   BOUND = p * (1 + eps);
    1194          28 :   if (BOUND > (double)ULONG_MAX || (ulong)BOUND != (ulong)p)
    1195           7 :     pari_err_PREC("forqfvec");
    1196             : 
    1197          21 :   k = n; y[n] = z[n] = 0;
    1198          21 :   x[n] = (long)sqrt(BOUND/v[n]);
    1199          56 :   for(;;x[1]--)
    1200             :   {
    1201             :     do
    1202             :     {
    1203         140 :       if (k>1)
    1204             :       {
    1205          84 :         long l = k-1;
    1206          84 :         z[l] = 0;
    1207         245 :         for (j=k; j<=n; j++) z[l] += q[l][j]*x[j];
    1208          84 :         p = (double)x[k] + z[k];
    1209          84 :         y[l] = y[k] + p*p*v[k];
    1210          84 :         x[l] = (long)floor(sqrt((BOUND-y[l])/v[l])-z[l]);
    1211          84 :         k = l;
    1212             :       }
    1213             :       for(;;)
    1214             :       {
    1215         189 :         p = (double)x[k] + z[k];
    1216         189 :         if (y[k] + p*p*v[k] <= BOUND) break;
    1217          49 :         k++; x[k]--;
    1218             :       }
    1219         140 :     } while (k > 1);
    1220          77 :     if (! x[1] && y[1]<=eps) break;
    1221             : 
    1222          56 :     p = (double)x[1] + z[1]; p = y[1] + p*p*v[1]; /* norm(x) */
    1223          56 :     if (fun(E, u, x, p)) break;
    1224             :   }
    1225             : }
    1226             : 
    1227             : void
    1228           0 : forqfvec(void *E, long (*fun)(void *, GEN, GEN, double), GEN a, GEN BORNE)
    1229             : {
    1230           0 :   pari_sp av = avma;
    1231             :   struct qfvec qv;
    1232           0 :   forqfvec_init(&qv, a);
    1233           0 :   forqfvec_i(E, fun, &qv, BORNE);
    1234           0 :   set_avma(av);
    1235           0 : }
    1236             : 
    1237             : struct qfvecwrap
    1238             : {
    1239             :   void *E;
    1240             :   long (*fun)(void *, GEN);
    1241             : };
    1242             : 
    1243             : static long
    1244          56 : forqfvec_wrap(void *E, GEN u, GEN x, double d)
    1245             : {
    1246          56 :   pari_sp av = avma;
    1247          56 :   struct qfvecwrap *W = (struct qfvecwrap *) E;
    1248             :   (void) d;
    1249          56 :   return gc_long(av, W->fun(W->E, ZM_zc_mul_canon(u, x)));
    1250             : }
    1251             : 
    1252             : void
    1253          42 : forqfvec1(void *E, long (*fun)(void *, GEN), GEN a, GEN BORNE)
    1254             : {
    1255          42 :   pari_sp av = avma;
    1256             :   struct qfvecwrap wr;
    1257             :   struct qfvec qv;
    1258          42 :   wr.E = E; wr.fun = fun;
    1259          42 :   forqfvec_init(&qv, a);
    1260          42 :   forqfvec_i((void*) &wr, forqfvec_wrap, &qv, BORNE);
    1261          35 :   set_avma(av);
    1262          35 : }
    1263             : 
    1264             : void
    1265          42 : forqfvec0(GEN a, GEN BORNE, GEN code)
    1266          42 : { EXPRVOID_WRAP(code, forqfvec1(EXPR_ARGVOID, a,  BORNE)) }
    1267             : 
    1268             : enum { min_ALL = 0, min_FIRST, min_VECSMALL, min_VECSMALL2 };
    1269             : 
    1270             : static int
    1271         923 : stockmax_init(const char *fun, GEN STOCKMAX, long *maxrank)
    1272             : {
    1273         923 :   long r = 200;
    1274         923 :   if (!STOCKMAX) { *maxrank = 200; return 1; }
    1275         511 :   STOCKMAX = gfloor(STOCKMAX);
    1276         511 :   if (typ(STOCKMAX) != t_INT) pari_err_TYPE(fun, STOCKMAX);
    1277         511 :   r = itos(STOCKMAX);
    1278         511 :   if (r < 0)
    1279             :   {
    1280           0 :     char *e = stack_strcat(fun, "[negative number of vectors]");
    1281           0 :     pari_err_TYPE(e, STOCKMAX);
    1282             :   }
    1283         511 :   *maxrank = r; return 0;
    1284             : }
    1285             : 
    1286             : /* Minimal vectors for the integral definite quadratic form: a.
    1287             :  * Result u:
    1288             :  *   u[1]= Number of vectors of square norm <= BORNE
    1289             :  *   u[2]= maximum norm found
    1290             :  *   u[3]= list of vectors found (at most STOCKMAX, unless NULL)
    1291             :  *
    1292             :  *  If BORNE = NULL: Minimal nonzero vectors.
    1293             :  *  flag = min_ALL,   as above
    1294             :  *  flag = min_FIRST, exits when first suitable vector is found.
    1295             :  *  flag = min_VECSMALL, return a t_VECSMALL of (half) the number of vectors
    1296             :  *  for each norm
    1297             :  *  flag = min_VECSMALL2, same but count only vectors with even norm, and shift
    1298             :  *  the answer */
    1299             : static GEN
    1300         847 : minim0_dolll(GEN a, GEN BORNE, GEN STOCKMAX, long flag, long dolll)
    1301             : {
    1302             :   GEN x, u, r, L, gnorme;
    1303         847 :   long n = lg(a)-1, i, j, k, s, maxrank, sBORNE;
    1304         847 :   pari_sp av = avma, av1;
    1305             :   double p,maxnorm,BOUND,*v,*y,*z,**q;
    1306         847 :   const double eps = 1e-10;
    1307             :   int stockall;
    1308             :   struct qfvec qv;
    1309             : 
    1310         847 :   if (!BORNE)
    1311          56 :     sBORNE = 0;
    1312             :   else
    1313             :   {
    1314         791 :     BORNE = gfloor(BORNE);
    1315         791 :     if (typ(BORNE) != t_INT) pari_err_TYPE("minim0",BORNE);
    1316         791 :     if (is_bigint(BORNE)) pari_err_PREC( "qfminim");
    1317         790 :     sBORNE = itos(BORNE); set_avma(av);
    1318         790 :     if (sBORNE < 0) sBORNE = 0;
    1319             :   }
    1320         846 :   stockall = stockmax_init("minim0", STOCKMAX, &maxrank);
    1321             : 
    1322         846 :   switch(flag)
    1323             :   {
    1324         462 :     case min_VECSMALL:
    1325             :     case min_VECSMALL2:
    1326         462 :       if (sBORNE <= 0) return cgetg(1, t_VECSMALL);
    1327         434 :       L = zero_zv(sBORNE);
    1328         434 :       if (flag == min_VECSMALL2) sBORNE <<= 1;
    1329         434 :       if (n == 0) return L;
    1330         434 :       break;
    1331          35 :     case min_FIRST:
    1332          35 :       if (n == 0 || (!sBORNE && BORNE)) return cgetg(1,t_VEC);
    1333          21 :       L = NULL; /* gcc -Wall */
    1334          21 :       break;
    1335         349 :     case min_ALL:
    1336         349 :       if (n == 0 || (!sBORNE && BORNE))
    1337          14 :         retmkvec3(gen_0, gen_0, cgetg(1, t_MAT));
    1338         335 :       L = new_chunk(1+maxrank);
    1339         335 :       break;
    1340           0 :     default:
    1341           0 :       return NULL;
    1342             :   }
    1343         790 :   minim_alloc(n+1, &q, &x, &y, &z, &v);
    1344             : 
    1345         790 :   forqfvec_init_dolll(&qv, &a, dolll);
    1346         790 :   av1 = avma;
    1347         790 :   r = qv.r;
    1348         790 :   u = qv.u;
    1349        5912 :   for (j=1; j<=n; j++)
    1350             :   {
    1351        5122 :     v[j] = rtodbl(gcoeff(r,j,j));
    1352       29579 :     for (i=1; i<j; i++) q[i][j] = rtodbl(gcoeff(r,i,j)); /* |.| <= 1/2 */
    1353             :   }
    1354             : 
    1355         790 :   if (sBORNE) maxnorm = 0.;
    1356             :   else
    1357             :   {
    1358          56 :     GEN B = gcoeff(a,1,1);
    1359          56 :     long t = 1;
    1360         616 :     for (i=2; i<=n; i++)
    1361             :     {
    1362         560 :       GEN c = gcoeff(a,i,i);
    1363         560 :       if (cmpii(c, B) < 0) { B = c; t = i; }
    1364             :     }
    1365          56 :     if (flag == min_FIRST) return gc_GEN(av, mkvec2(B, gel(u,t)));
    1366          49 :     maxnorm = -1.; /* don't update maxnorm */
    1367          49 :     if (is_bigint(B)) return NULL;
    1368          48 :     sBORNE = itos(B);
    1369             :   }
    1370         782 :   BOUND = sBORNE * (1 + eps);
    1371         782 :   if ((long)BOUND != sBORNE) return NULL;
    1372             : 
    1373         770 :   s = 0;
    1374         770 :   k = n; y[n] = z[n] = 0;
    1375         770 :   x[n] = (long)sqrt(BOUND/v[n]);
    1376     1223264 :   for(;;x[1]--)
    1377             :   {
    1378             :     do
    1379             :     {
    1380     2245614 :       if (k>1)
    1381             :       {
    1382     1022259 :         long l = k-1;
    1383     1022259 :         z[l] = 0;
    1384    11756080 :         for (j=k; j<=n; j++) z[l] += q[l][j]*x[j];
    1385     1022259 :         p = (double)x[k] + z[k];
    1386     1022259 :         y[l] = y[k] + p*p*v[k];
    1387     1022259 :         x[l] = (long)floor(sqrt((BOUND-y[l])/v[l])-z[l]);
    1388     1022259 :         k = l;
    1389             :       }
    1390             :       for(;;)
    1391             :       {
    1392     3263729 :         p = (double)x[k] + z[k];
    1393     3263729 :         if (y[k] + p*p*v[k] <= BOUND) break;
    1394     1018115 :         k++; x[k]--;
    1395             :       }
    1396             :     }
    1397     2245614 :     while (k > 1);
    1398     1224034 :     if (! x[1] && y[1]<=eps) break;
    1399             : 
    1400     1223271 :     p = (double)x[1] + z[1];
    1401     1223271 :     p = y[1] + p*p*v[1]; /* norm(x) */
    1402     1223271 :     if (maxnorm >= 0)
    1403             :     {
    1404     1220723 :       if (p > maxnorm) maxnorm = p;
    1405             :     }
    1406             :     else
    1407             :     { /* maxnorm < 0 : only look for minimal vectors */
    1408        2548 :       pari_sp av2 = avma;
    1409        2548 :       gnorme = roundr(dbltor(p));
    1410        2548 :       if (cmpis(gnorme, sBORNE) >= 0) set_avma(av2);
    1411             :       else
    1412             :       {
    1413          14 :         sBORNE = itos(gnorme); set_avma(av1);
    1414          14 :         BOUND = sBORNE * (1+eps);
    1415          14 :         L = new_chunk(maxrank+1);
    1416          14 :         s = 0;
    1417             :       }
    1418             :     }
    1419     1223271 :     s++;
    1420             : 
    1421     1223271 :     switch(flag)
    1422             :     {
    1423           7 :       case min_FIRST:
    1424           7 :         if (dolll) x = ZM_zc_mul_canon(u,x);
    1425           7 :         return gc_GEN(av, mkvec2(roundr(dbltor(p)), x));
    1426             : 
    1427      248241 :       case min_ALL:
    1428      248241 :         if (s > maxrank && stockall) /* overflow */
    1429             :         {
    1430         490 :           long maxranknew = maxrank << 1;
    1431         490 :           GEN Lnew = new_chunk(1 + maxranknew);
    1432      344890 :           for (i=1; i<=maxrank; i++) Lnew[i] = L[i];
    1433         490 :           L = Lnew; maxrank = maxranknew;
    1434             :         }
    1435      248241 :         if (s<=maxrank) gel(L,s) = leafcopy(x);
    1436      248241 :         break;
    1437             : 
    1438       39200 :       case min_VECSMALL:
    1439       39200 :         { ulong norm = (ulong)(p + 0.5); L[norm]++; }
    1440       39200 :         break;
    1441             : 
    1442      935823 :       case min_VECSMALL2:
    1443      935823 :         { ulong norm = (ulong)(p + 0.5); if (!odd(norm)) L[norm>>1]++; }
    1444      935823 :         break;
    1445             : 
    1446             :     }
    1447             :   }
    1448         763 :   switch(flag)
    1449             :   {
    1450           7 :     case min_FIRST:
    1451           7 :       retgc_const(av, cgetg(1, t_VEC));
    1452         434 :     case min_VECSMALL:
    1453             :     case min_VECSMALL2:
    1454         434 :       set_avma((pari_sp)L); return L;
    1455             :   }
    1456         322 :   r = (maxnorm >= 0) ? roundr(dbltor(maxnorm)): stoi(sBORNE);
    1457         322 :   k = minss(s,maxrank);
    1458         322 :   L[0] = evaltyp(t_MAT) | evallg(k + 1);
    1459         322 :   if (dolll)
    1460      246092 :     for (j=1; j<=k; j++)
    1461      245805 :       gel(L,j) = dolll==1 ? ZM_zc_mul_canon(u, gel(L,j))
    1462      245805 :                           : ZM_zc_mul_canon_zm(u, gel(L,j));
    1463         322 :   return gc_GEN(av, mkvec3(stoi(s<<1), r, L));
    1464             : }
    1465             : 
    1466             : /* Closest vectors for the integral definite quadratic form: a.
    1467             :  * Code bases on minim0_dolll
    1468             :  * Result u:
    1469             :  *   u[1]= Number of closest vectors of square distance <= BORNE
    1470             :  *   u[2]= maximum squared distance found
    1471             :  *   u[3]= list of vectors found (at most STOCKMAX, unless NULL)
    1472             :  *
    1473             :  *  If BORNE = NULL or <= 0.: returns closest vectors.
    1474             :  *  flag = min_ALL,   as above
    1475             :  *  flag = min_FIRST, exits when first suitable vector is found.
    1476             : */
    1477             : static GEN
    1478          91 : cvp0_dolll(GEN a, GEN target, GEN BORNE, GEN STOCKMAX, long flag, long dolll)
    1479             : {
    1480             :   GEN x, u, r, L;
    1481          91 :   long n = lg(a)-1, i, j, k, s, maxrank;
    1482          91 :   pari_sp av = avma, av1;
    1483             :   double p,maxnorm,BOUND,*v,*y,*z,*tt,**q, *tpre, sBORNE;
    1484          91 :   const double eps = 1e-10;
    1485             :   int stockall;
    1486             :   struct qfvec qv;
    1487          91 :   int done = 0;
    1488             : 
    1489          91 :   if (!is_vec_t(typ(target))) pari_err_TYPE("cvp0",target);
    1490          91 :   if (n != lg(target)-1) pari_err_TYPE("cvp0 [different dimensions]",target);
    1491          77 :   if (!BORNE)
    1492           0 :     sBORNE = 0.;
    1493             :   else
    1494             :   {
    1495          77 :     if (!is_real_t(typ(BORNE))) pari_err_TYPE("cvp0",BORNE);
    1496          77 :     sBORNE = gtodouble(BORNE);
    1497          77 :     if (sBORNE < 0.) sBORNE = 0.;
    1498             :   }
    1499          77 :   stockall = stockmax_init("cvp0", STOCKMAX, &maxrank);
    1500             : 
    1501          77 :   L = (flag==min_ALL) ? new_chunk(1+maxrank) : NULL;
    1502          77 :   if (n == 0)
    1503             :   {
    1504           7 :     if (flag==min_ALL) retmkvec3(gen_0, gen_0, cgetg(1, t_MAT));
    1505           0 :     return cgetg(1,t_VEC);
    1506             :   }
    1507             : 
    1508          70 :   minim_alloc(n+1, &q, &x, &y, &z, &v);
    1509          70 :   cvp_alloc(n+1, &tt, &tpre);
    1510             : 
    1511          70 :   forqfvec_init_dolll(&qv, &a, dolll);
    1512          70 :   av1 = avma;
    1513          70 :   r = qv.r;
    1514          70 :   u = qv.u;
    1515         392 :   for (j=1; j<=n; j++)
    1516             :   {
    1517         322 :     v[j] = rtodbl(gcoeff(r,j,j));
    1518        1729 :     for (i=1; i<j; i++) q[i][j] = rtodbl(gcoeff(r,i,j)); /* |.| <= 1/2 */
    1519             :   }
    1520             : 
    1521          70 :   if (dolll)
    1522             :   {
    1523          70 :     GEN tv = RgM_RgC_mul(ZM_inv(u, NULL), target);
    1524         392 :     for (j=1; j<=n; j++) tt[j] = gtodouble(gel(tv, j));
    1525             :   } else
    1526           0 :     for (j=1; j<=n; j++) tt[j] = gtodouble(gel(target, j));
    1527             :   /* precompute contribution of tt to z[l] */
    1528         392 :   for(k=1; k <= n; k++)
    1529             :   {
    1530         322 :     tpre[k] = -tt[k];
    1531        1729 :     for(j=k+1; j<=n; j++) tpre[k] -= q[k][j] * tt[j];
    1532             :   }
    1533             : 
    1534          70 :   if (sBORNE) maxnorm = 0.;
    1535             :   else
    1536             :   {
    1537          28 :     GEN B = gcoeff(a,1,1);
    1538         112 :     for (i = 2; i <= n; i++) B = addii(B, gcoeff(a,i,i));
    1539          28 :     maxnorm = -1.; /* don't update maxnorm */
    1540          28 :     if (is_bigint(B)) return NULL;
    1541          28 :     sBORNE = 0.;
    1542         140 :     for(i=1; i<=n; i++) sBORNE += v[i];
    1543             :   }
    1544          70 :   BOUND = sBORNE * (1 + eps);
    1545             : 
    1546          70 :   s = 0;
    1547          70 :   k = n; y[n] = 0;
    1548          70 :   z[n] = tpre[n];
    1549          70 :   x[n] = (long)floor(sqrt(BOUND/v[n])-z[n]);
    1550         889 :   for(;;x[1]--)
    1551             :   {
    1552             :     do
    1553             :     {
    1554        8582 :       if (k>1)
    1555             :       {
    1556        7665 :         long l = k-1;
    1557        7665 :         z[l] = tpre[l];
    1558       61488 :         for (j=k; j<=n; j++) z[l] += q[l][j]*x[j];
    1559        7665 :         p = (double)x[k] + z[k];
    1560        7665 :         y[l] = y[k] + p*p*v[k];
    1561        7665 :         x[l] = (long)floor(sqrt((BOUND-y[l])/v[l])-z[l]);
    1562        7665 :         k = l;
    1563             :       }
    1564             :       for(;;)
    1565             :       {
    1566       16247 :         p = (double)x[k] + z[k];
    1567       16247 :         if (y[k] + p*p*v[k] <= BOUND) break;
    1568        7735 :         if (k >= n) { done = 1; break; }
    1569        7665 :         k++; x[k]--;
    1570             :       }
    1571             :     }
    1572        8582 :     while (k > 1 && !done);
    1573         959 :     if (done) break;
    1574             : 
    1575         889 :     p = (double)x[1] + z[1];
    1576         889 :     p = y[1] + p*p*v[1]; /* norm(x-target) */
    1577         889 :     if (maxnorm >= 0)
    1578             :     {
    1579         175 :       if (p > maxnorm) maxnorm = p;
    1580             :     }
    1581             :     else
    1582             :     { /* maxnorm < 0 : only look for closest vectors */
    1583         714 :       if (p * (1+10*eps) < sBORNE) {
    1584         231 :         sBORNE = p; set_avma(av1);
    1585         231 :         BOUND = sBORNE * (1+eps);
    1586         231 :         L = new_chunk(maxrank+1);
    1587         231 :         s = 0;
    1588             :       }
    1589             :     }
    1590         889 :     s++;
    1591             : 
    1592         889 :     switch(flag)
    1593             :     {
    1594           0 :       case min_FIRST:
    1595           0 :         if (dolll) x = ZM_zc_mul(u,x);
    1596           0 :         return gc_GEN(av, mkvec2(dbltor(p), x));
    1597             : 
    1598         889 :       case min_ALL:
    1599         889 :         if (s > maxrank && stockall) /* overflow */
    1600             :         {
    1601           0 :           long maxranknew = maxrank << 1;
    1602           0 :           GEN Lnew = new_chunk(1 + maxranknew);
    1603           0 :           for (i=1; i<=maxrank; i++) Lnew[i] = L[i];
    1604           0 :           L = Lnew; maxrank = maxranknew;
    1605             :         }
    1606         889 :         if (s<=maxrank) gel(L,s) = leafcopy(x);
    1607         889 :         break;
    1608             :     }
    1609             :   }
    1610          70 :   switch(flag)
    1611             :   {
    1612           0 :     case min_FIRST:
    1613           0 :       retgc_const(av, cgetg(1, t_VEC));
    1614             :   }
    1615          70 :   r = (maxnorm >= 0) ? dbltor(maxnorm): dbltor(sBORNE);
    1616          70 :   k = minss(s,maxrank);
    1617          70 :   L[0] = evaltyp(t_MAT) | evallg(k + 1);
    1618         322 :   for (j=1; j<=k; j++)
    1619         252 :     gel(L,j) = dolll==1 ? ZM_zc_mul(u, gel(L,j))
    1620         252 :                         : zc_to_ZC(gel(L,j));
    1621          70 :   return gc_GEN(av, mkvec3(stoi(s), r, L));
    1622             : }
    1623             : 
    1624             : static GEN
    1625         553 : minim0(GEN a, GEN BORNE, GEN STOCKMAX, long flag)
    1626             : {
    1627         553 :   GEN v = minim0_dolll(a, BORNE, STOCKMAX, flag, 1);
    1628         552 :   if (!v) pari_err_PREC("qfminim");
    1629         546 :   return v;
    1630             : }
    1631             : 
    1632             : static GEN
    1633          91 : cvp0(GEN a, GEN target, GEN BORNE, GEN STOCKMAX, long flag)
    1634             : {
    1635          91 :   GEN v = cvp0_dolll(a, target, BORNE, STOCKMAX, flag, 1);
    1636          77 :   if (!v) pari_err_PREC("qfcvp");
    1637          77 :   return v;
    1638             : }
    1639             : 
    1640             : static GEN
    1641         252 : minim0_zm(GEN a, GEN BORNE, GEN STOCKMAX, long flag)
    1642             : {
    1643         252 :   GEN v = minim0_dolll(a, BORNE, STOCKMAX, flag, 2);
    1644         252 :   if (!v) pari_err_PREC("qfminim");
    1645         252 :   return v;
    1646             : }
    1647             : 
    1648             : GEN
    1649         462 : qfrep0(GEN a, GEN borne, long flag)
    1650         462 : { return minim0(a, borne, gen_0, (flag & 1)? min_VECSMALL2: min_VECSMALL); }
    1651             : 
    1652             : GEN
    1653         133 : qfminim0(GEN a, GEN borne, GEN stockmax, long flag, long prec)
    1654             : {
    1655         133 :   switch(flag)
    1656             :   {
    1657          49 :     case 0: return minim0(a,borne,stockmax,min_ALL);
    1658          35 :     case 1: return minim0(a,borne,gen_0   ,min_FIRST);
    1659          49 :     case 2:
    1660             :     {
    1661          49 :       long maxnum = -1;
    1662          49 :       if (typ(a) != t_MAT) pari_err_TYPE("qfminim",a);
    1663          49 :       if (stockmax) {
    1664          14 :         if (typ(stockmax) != t_INT) pari_err_TYPE("qfminim",stockmax);
    1665          14 :         maxnum = itos(stockmax);
    1666             :       }
    1667          49 :       a = fincke_pohst(a,borne,maxnum,prec,NULL);
    1668          42 :       if (!a) pari_err_PREC("qfminim");
    1669          42 :       return a;
    1670             :     }
    1671           0 :     default: pari_err_FLAG("qfminim");
    1672             :   }
    1673             :   return NULL; /* LCOV_EXCL_LINE */
    1674             : }
    1675             : 
    1676             : 
    1677             : GEN
    1678          91 : qfcvp0(GEN a, GEN target, GEN borne, GEN stockmax, long flag)
    1679             : {
    1680          91 :   switch(flag)
    1681             :   {
    1682          91 :     case 0: return cvp0(a,target,borne,stockmax,min_ALL);
    1683           0 :     case 1: return cvp0(a,target,borne,gen_0   ,min_FIRST);
    1684             :     /* case 2:
    1685             :        TODO: more robust finke_pohst enumeration */
    1686           0 :     default: pari_err_FLAG("qfcvp");
    1687             :   }
    1688             :   return NULL; /* LCOV_EXCL_LINE */
    1689             : }
    1690             : 
    1691             : GEN
    1692           7 : minim(GEN a, GEN borne, GEN stockmax)
    1693           7 : { return minim0(a,borne,stockmax,min_ALL); }
    1694             : 
    1695             : GEN
    1696         252 : minim_zm(GEN a, GEN borne, GEN stockmax)
    1697         252 : { return minim0_zm(a,borne,stockmax,min_ALL); }
    1698             : 
    1699             : GEN
    1700          42 : minim_raw(GEN a, GEN BORNE, GEN STOCKMAX)
    1701          42 : { return minim0_dolll(a, BORNE, STOCKMAX, min_ALL, 0); }
    1702             : 
    1703             : GEN
    1704           0 : minim2(GEN a, GEN borne, GEN stockmax)
    1705           0 : { return minim0(a,borne,stockmax,min_FIRST); }
    1706             : 
    1707             : /* If V depends linearly from the columns of the matrix, return 0.
    1708             :  * Otherwise, update INVP and L and return 1. No GC. */
    1709             : static int
    1710        1652 : addcolumntomatrix(GEN V, GEN invp, GEN L)
    1711             : {
    1712        1652 :   long i,j,k, n = lg(invp);
    1713        1652 :   GEN a = cgetg(n, t_COL), ak = NULL, mak;
    1714             : 
    1715       84231 :   for (k = 1; k < n; k++)
    1716       83706 :     if (!L[k])
    1717             :     {
    1718       27902 :       ak = RgMrow_zc_mul(invp, V, k);
    1719       27902 :       if (!gequal0(ak)) break;
    1720             :     }
    1721        1652 :   if (k == n) return 0;
    1722        1127 :   L[k] = 1;
    1723        1127 :   mak = gneg_i(ak);
    1724       43253 :   for (i=k+1; i<n; i++)
    1725       42126 :     gel(a,i) = gdiv(RgMrow_zc_mul(invp, V, i), mak);
    1726       43883 :   for (j=1; j<=k; j++)
    1727             :   {
    1728       42756 :     GEN c = gel(invp,j), ck = gel(c,k);
    1729       42756 :     if (gequal0(ck)) continue;
    1730        8757 :     gel(c,k) = gdiv(ck, ak);
    1731        8757 :     if (j==k)
    1732       43253 :       for (i=k+1; i<n; i++)
    1733       42126 :         gel(c,i) = gmul(gel(a,i), ck);
    1734             :     else
    1735      184814 :       for (i=k+1; i<n; i++)
    1736      177184 :         gel(c,i) = gadd(gel(c,i), gmul(gel(a,i), ck));
    1737             :   }
    1738        1127 :   return 1;
    1739             : }
    1740             : 
    1741             : GEN
    1742          42 : qfperfection(GEN a)
    1743             : {
    1744          42 :   pari_sp av = avma;
    1745             :   GEN u, L;
    1746          42 :   long r, s, k, l, n = lg(a)-1;
    1747             : 
    1748          42 :   if (!n) return gen_0;
    1749          42 :   if (typ(a) != t_MAT || !RgM_is_ZM(a)) pari_err_TYPE("qfperfection",a);
    1750          42 :   a = minim_lll(a, &u);
    1751          42 :   L = minim_raw(a,NULL,NULL);
    1752          42 :   r = (n*(n+1)) >> 1;
    1753          42 :   if (L)
    1754             :   {
    1755             :     GEN D, V, invp;
    1756          35 :     L = gel(L, 3); l = lg(L);
    1757          35 :     if (l == 2) { set_avma(av); return gen_1; }
    1758             :     /* |L[i]|^2 fits  into a long for all i */
    1759          21 :     D = zero_zv(r);
    1760          21 :     V = cgetg(r+1, t_VECSMALL);
    1761          21 :     invp = matid(r);
    1762          21 :     s = 0;
    1763        1659 :     for (k = 1; k < l; k++)
    1764             :     {
    1765        1652 :       pari_sp av2 = avma;
    1766        1652 :       GEN x = gel(L,k);
    1767             :       long i, j, I;
    1768       21098 :       for (i = I = 1; i<=n; i++)
    1769      145278 :         for (j=i; j<=n; j++,I++) V[I] = x[i]*x[j];
    1770        1652 :       if (!addcolumntomatrix(V,invp,D)) set_avma(av2);
    1771        1127 :       else if (++s == r) break;
    1772             :     }
    1773             :   }
    1774             :   else
    1775             :   {
    1776             :     GEN M;
    1777           7 :     L = fincke_pohst(a,NULL,-1, DEFAULTPREC, NULL);
    1778           7 :     if (!L) pari_err_PREC("qfminim");
    1779           7 :     L = gel(L, 3); l = lg(L);
    1780           7 :     if (l == 2) { set_avma(av); return gen_1; }
    1781           7 :     M = cgetg(l, t_MAT);
    1782         959 :     for (k = 1; k < l; k++)
    1783             :     {
    1784         952 :       GEN x = gel(L,k), c = cgetg(r+1, t_COL);
    1785             :       long i, I, j;
    1786       16184 :       for (i = I = 1; i<=n; i++)
    1787      144704 :         for (j=i; j<=n; j++,I++) gel(c,I) = mulii(gel(x,i), gel(x,j));
    1788         952 :       gel(M,k) = c;
    1789             :     }
    1790           7 :     s = ZM_rank(M);
    1791             :   }
    1792          28 :   return gc_utoipos(av, s);
    1793             : }
    1794             : 
    1795             : static GEN
    1796         141 : clonefill(GEN S, long s, long t)
    1797             : { /* initialize to dummy values */
    1798         141 :   GEN T = S, dummy = cgetg(1, t_STR);
    1799             :   long i;
    1800      310917 :   for (i = s+1; i <= t; i++) gel(S,i) = dummy;
    1801         141 :   S = gclone(S); if (isclone(T)) gunclone(T);
    1802         141 :   return S;
    1803             : }
    1804             : 
    1805             : /* increment ZV x, by incrementing cell of index k. Initial value x0[k] was
    1806             :  * chosen to minimize qf(x) for given x0[1..k-1] and x0[k+1,..] = 0
    1807             :  * The last nonzero entry must be positive and goes through x0[k]+1,2,3,...
    1808             :  * Others entries go through: x0[k]+1,-1,2,-2,...*/
    1809             : INLINE void
    1810     2952936 : step(GEN x, GEN y, GEN inc, long k)
    1811             : {
    1812     2952936 :   if (!signe(gel(y,k))) /* x[k+1..] = 0 */
    1813      160816 :     gel(x,k) = addiu(gel(x,k), 1); /* leading coeff > 0 */
    1814             :   else
    1815             :   {
    1816     2792120 :     long i = inc[k];
    1817     2792120 :     gel(x,k) = addis(gel(x,k), i),
    1818     2792124 :     inc[k] = (i > 0)? -1-i: 1-i;
    1819             :   }
    1820     2952938 : }
    1821             : 
    1822             : /* 1 if we are "sure" that x < y, up to few rounding errors, i.e.
    1823             :  * x < y - epsilon. More precisely :
    1824             :  * - sign(x - y) < 0
    1825             :  * - lgprec(x-y) > 3 || expo(x - y) - expo(x) > -24 */
    1826             : static int
    1827     1216483 : mplessthan(GEN x, GEN y)
    1828             : {
    1829     1216483 :   pari_sp av = avma;
    1830     1216483 :   GEN z = mpsub(x, y);
    1831     1216479 :   set_avma(av);
    1832     1216479 :   if (typ(z) == t_INT) return (signe(z) < 0);
    1833     1216479 :   if (signe(z) >= 0) return 0;
    1834       22157 :   if (realprec(z) > LOWDEFAULTPREC) return 1;
    1835       22157 :   return ( expo(z) - mpexpo(x) > -24 );
    1836             : }
    1837             : 
    1838             : /* 1 if we are "sure" that x > y, up to few rounding errors, i.e.
    1839             :  * x > y + epsilon */
    1840             : static int
    1841     4621588 : mpgreaterthan(GEN x, GEN y)
    1842             : {
    1843     4621588 :   pari_sp av = avma;
    1844     4621588 :   GEN z = mpsub(x, y);
    1845     4621607 :   set_avma(av);
    1846     4621680 :   if (typ(z) == t_INT) return (signe(z) > 0);
    1847     4621680 :   if (signe(z) <= 0) return 0;
    1848     2689986 :   if (realprec(z) > LOWDEFAULTPREC) return 1;
    1849      476056 :   return ( expo(z) - mpexpo(x) > -24 );
    1850             : }
    1851             : 
    1852             : /* x a t_INT, y  t_INT or t_REAL */
    1853             : INLINE GEN
    1854     1228574 : mulimp(GEN x, GEN y)
    1855             : {
    1856     1228574 :   if (typ(y) == t_INT) return mulii(x,y);
    1857     1228574 :   return signe(x) ? mulir(x,y): gen_0;
    1858             : }
    1859             : /* x + y*z, x,z two mp's, y a t_INT */
    1860             : INLINE GEN
    1861    13538822 : addmulimp(GEN x, GEN y, GEN z)
    1862             : {
    1863    13538822 :   if (!signe(y)) return x;
    1864     5831055 :   if (typ(z) == t_INT) return mpadd(x, mulii(y, z));
    1865     5831055 :   return mpadd(x, mulir(y, z));
    1866             : }
    1867             : 
    1868             : /* yk + vk * (xk + zk)^2 */
    1869             : static GEN
    1870     5780555 : norm_aux(GEN xk, GEN yk, GEN zk, GEN vk)
    1871             : {
    1872     5780555 :   GEN t = mpadd(xk, zk);
    1873     5780547 :   if (typ(t) == t_INT) { /* probably gen_0, avoid loss of accuracy */
    1874      306219 :     yk = addmulimp(yk, sqri(t), vk);
    1875             :   } else {
    1876     5474328 :     yk = mpadd(yk, mpmul(sqrr(t), vk));
    1877             :   }
    1878     5780484 :   return yk;
    1879             : }
    1880             : /* yk + vk * (xk + zk)^2 < B + epsilon */
    1881             : static int
    1882     4167488 : check_bound(GEN B, GEN xk, GEN yk, GEN zk, GEN vk)
    1883             : {
    1884     4167488 :   pari_sp av = avma;
    1885     4167488 :   int f = mpgreaterthan(norm_aux(xk,yk,zk,vk), B);
    1886     4167475 :   return gc_bool(av, !f);
    1887             : }
    1888             : 
    1889             : /* q(k-th canonical basis vector), where q is given in Cholesky form
    1890             :  * q(x) = sum_{i = 1}^n q[i,i] (x[i] + sum_{j > i} q[i,j] x[j])^2.
    1891             :  * Namely q(e_k) = q[k,k] + sum_{i < k} q[i,i] q[i,k]^2
    1892             :  * Assume 1 <= k <= n. */
    1893             : static GEN
    1894         182 : cholesky_norm_ek(GEN q, long k)
    1895             : {
    1896         182 :   GEN t = gcoeff(q,k,k);
    1897             :   long i;
    1898        1484 :   for (i = 1; i < k; i++) t = norm_aux(gen_0, t, gcoeff(q,i,k), gcoeff(q,i,i));
    1899         182 :   return t;
    1900             : }
    1901             : 
    1902             : /* q is the Cholesky decomposition of a quadratic form
    1903             :  * Enumerate vectors whose norm is less than BORNE (Algo 2.5.7),
    1904             :  * minimal vectors if BORNE = NULL (implies check = NULL).
    1905             :  * If (check != NULL) consider only vectors passing the check, and assumes
    1906             :  *   we only want the smallest possible vectors */
    1907             : static GEN
    1908       14713 : smallvectors(GEN q, GEN BORNE, long maxnum, FP_chk_fun *CHECK)
    1909             : {
    1910       14713 :   long N = lg(q), n = N-1, i, j, k, s, stockmax, checkcnt = 1;
    1911             :   pari_sp av, av1;
    1912             :   GEN inc, S, x, y, z, v, p1, alpha, norms;
    1913             :   GEN norme1, normax1, borne1, borne2;
    1914       14713 :   GEN (*check)(void *,GEN) = CHECK? CHECK->f: NULL;
    1915       14713 :   void *data = CHECK? CHECK->data: NULL;
    1916       14713 :   const long skipfirst = CHECK? CHECK->skipfirst: 0;
    1917       14713 :   const int stockall = (maxnum == -1);
    1918             : 
    1919       14713 :   alpha = dbltor(0.95);
    1920       14713 :   normax1 = gen_0;
    1921             : 
    1922       14713 :   v = cgetg(N,t_VEC);
    1923       14713 :   inc = const_vecsmall(n, 1);
    1924             : 
    1925       14713 :   av = avma;
    1926       14713 :   stockmax = stockall? 2000: maxnum;
    1927       14713 :   norms = cgetg(check?(stockmax+1): 1,t_VEC); /* unused if (!check) */
    1928       14713 :   S = cgetg(stockmax+1,t_VEC);
    1929       14713 :   x = cgetg(N,t_COL);
    1930       14713 :   y = cgetg(N,t_COL);
    1931       14713 :   z = cgetg(N,t_COL);
    1932       97807 :   for (i=1; i<N; i++) {
    1933       83094 :     gel(v,i) = gcoeff(q,i,i);
    1934       83094 :     gel(x,i) = gel(y,i) = gel(z,i) = gen_0;
    1935             :   }
    1936       14713 :   if (BORNE)
    1937             :   {
    1938       14692 :     borne1 = BORNE;
    1939       14692 :     if (gsigne(borne1) <= 0) retmkvec3(gen_0, gen_0, cgetg(1,t_MAT));
    1940       14678 :     if (typ(borne1) != t_REAL)
    1941             :     {
    1942             :       long prec;
    1943         419 :       prec = nbits2prec(gexpo(borne1) + 10);
    1944         419 :       borne1 = gtofp(borne1, maxss(prec, DEFAULTPREC));
    1945             :     }
    1946             :   }
    1947             :   else
    1948             :   {
    1949          21 :     borne1 = gcoeff(q,1,1);
    1950         203 :     for (i=2; i<N; i++)
    1951             :     {
    1952         182 :       GEN b = cholesky_norm_ek(q, i);
    1953         182 :       if (gcmp(b, borne1) < 0) borne1 = b;
    1954             :     }
    1955             :     /* borne1 = norm of smallest basis vector */
    1956             :   }
    1957       14699 :   borne2 = mulrr(borne1,alpha);
    1958       14699 :   if (DEBUGLEVEL>2)
    1959           0 :     err_printf("smallvectors looking for norm < %P.4G\n",borne1);
    1960       14699 :   s = 0; k = n;
    1961      383970 :   for(;; step(x,y,inc,k)) /* main */
    1962             :   { /* x (supposedly) small vector, ZV.
    1963             :      * For all t >= k, we have
    1964             :      *   z[t] = sum_{j > t} q[t,j] * x[j]
    1965             :      *   y[t] = sum_{i > t} q[i,i] * (x[i] + z[i])^2
    1966             :      *        = 0 <=> x[i]=0 for all i>t */
    1967             :     do
    1968             :     {
    1969     1612537 :       int skip = 0;
    1970     1612537 :       if (k > 1)
    1971             :       {
    1972     1228575 :         long l = k-1;
    1973     1228575 :         av1 = avma;
    1974     1228575 :         p1 = mulimp(gel(x,k), gcoeff(q,l,k));
    1975    14461211 :         for (j=k+1; j<N; j++) p1 = addmulimp(p1, gel(x,j), gcoeff(q,l,j));
    1976     1228584 :         gel(z,l) = gc_leaf(av1,p1);
    1977             : 
    1978     1228577 :         av1 = avma;
    1979     1228577 :         p1 = norm_aux(gel(x,k), gel(y,k), gel(z,k), gel(v,k));
    1980     1228574 :         gel(y,l) = gc_leaf(av1, p1);
    1981             :         /* skip the [x_1,...,x_skipfirst,0,...,0] */
    1982     1228576 :         if ((l <= skipfirst && !signe(gel(y,skipfirst)))
    1983     1228576 :          || mplessthan(borne1, gel(y,l))) skip = 1;
    1984             :         else /* initial value, minimizing (x[l] + z[l])^2, hence qf(x) for
    1985             :                 the given x[1..l-1] */
    1986     1214589 :           gel(x,l) = mpround( mpneg(gel(z,l)) );
    1987     1228577 :         k = l;
    1988             :       }
    1989     1228571 :       for(;; step(x,y,inc,k))
    1990             :       { /* at most 2n loops */
    1991     2841103 :         if (!skip)
    1992             :         {
    1993     2827117 :           if (check_bound(borne1, gel(x,k),gel(y,k),gel(z,k),gel(v,k))) break;
    1994     1340408 :           step(x,y,inc,k);
    1995     1340424 :           if (check_bound(borne1, gel(x,k),gel(y,k),gel(z,k),gel(v,k))) break;
    1996             :         }
    1997     1243270 :         skip = 0; inc[k] = 1;
    1998     1243270 :         if (++k > n) goto END;
    1999             :       }
    2000             : 
    2001     1597848 :       if (gc_needed(av,2))
    2002             :       {
    2003          15 :         if(DEBUGMEM>1) pari_warn(warnmem,"smallvectors");
    2004          15 :         if (stockmax) S = clonefill(S, s, stockmax);
    2005          15 :         if (check) {
    2006          15 :           GEN dummy = cgetg(1, t_STR);
    2007        9629 :           for (i=s+1; i<=stockmax; i++) gel(norms,i) = dummy;
    2008             :         }
    2009          15 :         (void)gc_all(av,7,&x,&y,&z,&normax1,&borne1,&borne2,&norms);
    2010             :       }
    2011             :     }
    2012     1597848 :     while (k > 1);
    2013      383969 :     if (!signe(gel(x,1)) && !signe(gel(y,1))) continue; /* exclude 0 */
    2014             : 
    2015      383254 :     av1 = avma;
    2016      383254 :     norme1 = norm_aux(gel(x,1),gel(y,1),gel(z,1),gel(v,1));
    2017      383255 :     if (mpgreaterthan(norme1,borne1)) { set_avma(av1); continue; /* main */ }
    2018             : 
    2019      383255 :     norme1 = gc_leaf(av1,norme1);
    2020      383255 :     if (check)
    2021             :     {
    2022      314669 :       if (checkcnt < 5 && mpcmp(norme1, borne2) < 0)
    2023             :       {
    2024        4418 :         if (!check(data,x)) { checkcnt++ ; continue; /* main */}
    2025         476 :         if (DEBUGLEVEL>4) err_printf("New bound: %Ps", norme1);
    2026         476 :         borne1 = norme1;
    2027         476 :         borne2 = mulrr(borne1, alpha);
    2028         476 :         s = 0; checkcnt = 0;
    2029             :       }
    2030             :     }
    2031             :     else
    2032             :     {
    2033       68586 :       if (!BORNE) /* find minimal vectors */
    2034             :       {
    2035        1890 :         if (mplessthan(norme1, borne1))
    2036             :         { /* strictly smaller vector than previously known */
    2037           0 :           borne1 = norme1; /* + epsilon */
    2038           0 :           s = 0;
    2039             :         }
    2040             :       }
    2041             :       else
    2042       66696 :         if (mpcmp(norme1,normax1) > 0) normax1 = norme1;
    2043             :     }
    2044      379313 :     if (++s > stockmax) continue; /* too many vectors: no longer remember */
    2045      378382 :     if (check) gel(norms,s) = norme1;
    2046      378382 :     gel(S,s) = leafcopy(x);
    2047      378382 :     if (s != stockmax) continue; /* still room, get next vector */
    2048             : 
    2049         126 :     if (check)
    2050             :     { /* overflow, eliminate vectors failing "check" */
    2051         105 :       pari_sp av2 = avma;
    2052             :       long imin, imax;
    2053         105 :       GEN per = indexsort(norms), S2 = cgetg(stockmax+1, t_VEC);
    2054         105 :       if (DEBUGLEVEL>2) err_printf("sorting... [%ld elts]\n",s);
    2055             :       /* let N be the minimal norm so far for x satisfying 'check'. Keep
    2056             :        * all elements of norm N */
    2057       26593 :       for (i = 1; i <= s; i++)
    2058             :       {
    2059       26586 :         long k = per[i];
    2060       26586 :         if (check(data,gel(S,k))) { borne1 = gel(norms,k); break; }
    2061             :       }
    2062         105 :       imin = i;
    2063       20943 :       for (; i <= s; i++)
    2064       20922 :         if (mpgreaterthan(gel(norms,per[i]), borne1)) break;
    2065         105 :       imax = i;
    2066       20943 :       for (i=imin, s=0; i < imax; i++) gel(S2,++s) = gel(S,per[i]);
    2067       20943 :       for (i = 1; i <= s; i++) gel(S,i) = gel(S2,i);
    2068         105 :       set_avma(av2);
    2069         105 :       if (s) { borne2 = mulrr(borne1, alpha); checkcnt = 0; }
    2070         105 :       if (!stockall) continue;
    2071         105 :       if (s > stockmax/2) stockmax <<= 1;
    2072         105 :       norms = cgetg(stockmax+1, t_VEC);
    2073       20943 :       for (i = 1; i <= s; i++) gel(norms,i) = borne1;
    2074             :     }
    2075             :     else
    2076             :     {
    2077          21 :       if (!stockall && BORNE) goto END;
    2078          21 :       if (!stockall) continue;
    2079          21 :       stockmax <<= 1;
    2080             :     }
    2081             : 
    2082             :     {
    2083         126 :       GEN Snew = clonefill(vec_lengthen(S,stockmax), s, stockmax);
    2084         126 :       if (isclone(S)) gunclone(S);
    2085         126 :       S = Snew;
    2086             :     }
    2087             :   }
    2088       14699 : END:
    2089       14699 :   if (s < stockmax) stockmax = s;
    2090       14699 :   if (check)
    2091             :   {
    2092             :     GEN per, alph, pols, p;
    2093       14671 :     if (DEBUGLEVEL>2) err_printf("final sort & check...\n");
    2094       14671 :     setlg(norms,stockmax+1); per = indexsort(norms);
    2095       14671 :     alph = cgetg(stockmax+1,t_VEC);
    2096       14671 :     pols = cgetg(stockmax+1,t_VEC);
    2097       84520 :     for (j=0,i=1; i<=stockmax; i++)
    2098             :     {
    2099       70116 :       long t = per[i];
    2100       70116 :       GEN N = gel(norms,t);
    2101       70116 :       if (j && mpgreaterthan(N, borne1)) break;
    2102       69849 :       if ((p = check(data,gel(S,t))))
    2103             :       {
    2104       55908 :         if (!j) borne1 = N;
    2105       55908 :         j++;
    2106       55908 :         gel(pols,j) = p;
    2107       55908 :         gel(alph,j) = gel(S,t);
    2108             :       }
    2109             :     }
    2110       14671 :     setlg(pols,j+1);
    2111       14671 :     setlg(alph,j+1);
    2112       14671 :     if (stockmax && isclone(S)) { alph = gcopy(alph); gunclone(S); }
    2113       14671 :     return mkvec2(pols, alph);
    2114             :   }
    2115          28 :   if (stockmax)
    2116             :   {
    2117          21 :     setlg(S,stockmax+1);
    2118          21 :     settyp(S,t_MAT);
    2119          21 :     if (isclone(S)) { p1 = S; S = gcopy(S); gunclone(p1); }
    2120             :   }
    2121             :   else
    2122           7 :     S = cgetg(1,t_MAT);
    2123          28 :   return mkvec3(utoi(s<<1), borne1, S);
    2124             : }
    2125             : 
    2126             : /* solve q(x) = x~.a.x <= bound, a > 0.
    2127             :  * If check is non-NULL keep x only if check(x).
    2128             :  * If a is a vector, assume a[1] is the LLL-reduced Cholesky form of q */
    2129             : GEN
    2130       14734 : fincke_pohst(GEN a, GEN B0, long stockmax, long PREC, FP_chk_fun *CHECK)
    2131             : {
    2132       14734 :   pari_sp av = avma;
    2133             :   VOLATILE long i,j,l;
    2134       14734 :   VOLATILE GEN r,rinv,rinvtrans,u,v,res,z,vnorm,rperm,perm,uperm, bound = B0;
    2135             : 
    2136       14734 :   if (typ(a) == t_VEC)
    2137             :   {
    2138       14266 :     r = gel(a,1);
    2139       14266 :     u = NULL;
    2140             :   }
    2141             :   else
    2142             :   {
    2143         468 :     long prec = PREC;
    2144         468 :     l = lg(a);
    2145         468 :     if (l == 1)
    2146             :     {
    2147           7 :       if (CHECK) pari_err_TYPE("fincke_pohst [dimension 0]", a);
    2148           7 :       retmkvec3(gen_0, gen_0, cgetg(1,t_MAT));
    2149             :     }
    2150         461 :     u = lllfp(a, 0.75, LLL_GRAM | LLL_IM);
    2151         454 :     if (!u || lg(u) != lg(a)) return gc_NULL(av);
    2152         454 :     r = qf_RgM_apply(a,u);
    2153         454 :     i = gprecision(r);
    2154         454 :     if (i)
    2155         412 :       prec = i;
    2156             :     else {
    2157          42 :       prec = DEFAULTPREC + nbits2extraprec(gexpo(r));
    2158          42 :       if (prec < PREC) prec = PREC;
    2159             :     }
    2160         454 :     if (DEBUGLEVEL>2) err_printf("first LLL: prec = %ld\n", prec);
    2161         454 :     r = qfgaussred_positive(r);
    2162         454 :     if (!r) return gc_NULL(av);
    2163        1984 :     for (i=1; i<l; i++)
    2164             :     {
    2165        1530 :       GEN s = gsqrt(gcoeff(r,i,i), prec);
    2166        1530 :       gcoeff(r,i,i) = s;
    2167        4236 :       for (j=i+1; j<l; j++) gcoeff(r,i,j) = gmul(s, gcoeff(r,i,j));
    2168             :     }
    2169             :   }
    2170             :   /* now r~ * r = a in LLL basis */
    2171       14720 :   rinv = RgM_inv_upper(r);
    2172       14720 :   if (!rinv) return gc_NULL(av);
    2173       14720 :   rinvtrans = shallowtrans(rinv);
    2174       14720 :   if (DEBUGLEVEL>2)
    2175           0 :     err_printf("Fincke-Pohst, final LLL: prec = %ld\n", gprecision(rinvtrans));
    2176       14720 :   v = lll(rinvtrans);
    2177       14720 :   if (lg(v) != lg(rinvtrans)) return gc_NULL(av);
    2178             : 
    2179       14720 :   rinvtrans = RgM_mul(rinvtrans, v);
    2180       14720 :   v = ZM_inv(shallowtrans(v),NULL);
    2181       14720 :   r = RgM_mul(r,v);
    2182       14720 :   u = u? ZM_mul(u,v): v;
    2183             : 
    2184       14720 :   l = lg(r);
    2185       14720 :   vnorm = cgetg(l,t_VEC);
    2186       97841 :   for (j=1; j<l; j++) gel(vnorm,j) = gnorml2(gel(rinvtrans,j));
    2187       14720 :   rperm = cgetg(l,t_MAT);
    2188       14720 :   uperm = cgetg(l,t_MAT); perm = indexsort(vnorm);
    2189       97842 :   for (i=1; i<l; i++) { uperm[l-i] = u[perm[i]]; rperm[l-i] = r[perm[i]]; }
    2190       14720 :   u = uperm;
    2191       14720 :   r = rperm; res = NULL;
    2192       14720 :   pari_CATCH(e_PREC) { }
    2193             :   pari_TRY {
    2194             :     GEN q;
    2195       14720 :     if (CHECK && CHECK->f_init) bound = CHECK->f_init(CHECK, r, u);
    2196       14713 :     q = gaussred_from_QR(r, gprecision(vnorm));
    2197       14713 :     if (q) res = smallvectors(q, bound, stockmax, CHECK);
    2198       14713 :   } pari_ENDCATCH;
    2199       14720 :   if (!res) return gc_NULL(av);
    2200       14713 :   if (CHECK)
    2201             :   {
    2202       14671 :     if (CHECK->f_post) res = CHECK->f_post(CHECK, res, u);
    2203       14671 :     return res;
    2204             :   }
    2205             : 
    2206          42 :   z = cgetg(4,t_VEC);
    2207          42 :   gel(z,1) = gcopy(gel(res,1));
    2208          42 :   gel(z,2) = gcopy(gel(res,2));
    2209          42 :   gel(z,3) = ZM_mul(u, gel(res,3)); return gc_upto(av,z);
    2210             : }

Generated by: LCOV version 1.16