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

Generated by: LCOV version 1.14