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.2 lcov report (development 29115-f22e516b23) Lines: 1085 1137 95.4 %
Date: 2024-03-29 08:06:26 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    21096819 : no_prec_pb(GEN x)
      30             : {
      31    21007594 :   return (typ(x) != t_REAL || realprec(x) > DEFAULTPREC
      32    42104413 :                            || expo(x) < DEFAULTPREC>>1);
      33             : }
      34             : /* Find a Householder transformation which, applied to x[k..#x], zeroes
      35             :  * x[k+1..#x]; fill L = (mu_{i,j}). Return 0 if precision problem [obtained
      36             :  * a 0 vector], 1 otherwise */
      37             : static int
      38    21106115 : FindApplyQ(GEN x, GEN L, GEN B, long k, GEN Q, long prec)
      39             : {
      40    21106115 :   long i, lx = lg(x)-1;
      41    21106115 :   GEN x2, x1, xd = x + (k-1);
      42             : 
      43    21106115 :   x1 = gel(xd,1);
      44    21106115 :   x2 = mpsqr(x1);
      45    21105038 :   if (k < lx)
      46             :   {
      47    16887310 :     long lv = lx - (k-1) + 1;
      48    16887310 :     GEN beta, Nx, v = cgetg(lv, t_VEC);
      49    69664824 :     for (i=2; i<lv; i++) {
      50    52777899 :       x2 = mpadd(x2, mpsqr(gel(xd,i)));
      51    52777092 :       gel(v,i) = gel(xd,i);
      52             :     }
      53    16886925 :     if (!signe(x2)) return 0;
      54    16878339 :     Nx = gsqrt(x2, prec); if (signe(x1) < 0) setsigne(Nx, -1);
      55    16879145 :     gel(v,1) = mpadd(x1, Nx);
      56             : 
      57    16878498 :     if (!signe(x1))
      58      647799 :       beta = gtofp(x2, prec); /* make sure typ(beta) != t_INT */
      59             :     else
      60    16230699 :       beta = mpadd(x2, mpmul(Nx,x1));
      61    16878828 :     gel(Q,k) = mkvec2(invr(beta), v);
      62             : 
      63    16879085 :     togglesign(Nx);
      64    16878854 :     gcoeff(L,k,k) = Nx;
      65             :   }
      66             :   else
      67     4217728 :     gcoeff(L,k,k) = gel(x,k);
      68    21096582 :   gel(B,k) = x2;
      69    63005059 :   for (i=1; i<k; i++) gcoeff(L,k,i) = gel(x,i);
      70    21096582 :   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    41918628 : ApplyQ(GEN Q, GEN r)
      77             : {
      78    41918628 :   GEN s, rd, beta = gel(Q,1), v = gel(Q,2);
      79    41918628 :   long i, l = lg(v), lr = lg(r);
      80             : 
      81    41918628 :   rd = r + (lr - l);
      82    41918628 :   s = mpmul(gel(v,1), gel(rd,1));
      83   493093271 :   for (i=2; i<l; i++) s = mpadd(s, mpmul(gel(v,i) ,gel(rd,i)));
      84    41914081 :   s = mpmul(beta, s);
      85   535500167 :   for (i=1; i<l; i++)
      86   493588813 :     if (signe(gel(v,i))) gel(rd,i) = mpsub(gel(rd,i), mpmul(s, gel(v,i)));
      87    41911354 : }
      88             : /* apply Q[1], ..., Q[j-1] to r */
      89             : static GEN
      90    14421263 : ApplyAllQ(GEN Q, GEN r, long j)
      91             : {
      92    14421263 :   pari_sp av = avma;
      93             :   long i;
      94    14421263 :   r = leafcopy(r);
      95    56337233 :   for (i=1; i<j; i++) ApplyQ(gel(Q,i), r);
      96    14418965 :   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     6684554 : QR_init(GEN x, GEN *pB, GEN *pQ, GEN *pL, long prec)
     188             : {
     189     6684554 :   long j, k = lg(x)-1;
     190     6684554 :   GEN B = cgetg(k+1, t_VEC), Q = cgetg(k, t_VEC), L = zeromatcopy(k,k);
     191    25545650 :   for (j=1; j<=k; j++)
     192             :   {
     193    21105757 :     GEN r = gel(x,j);
     194    21105757 :     if (j > 1) r = ApplyAllQ(Q, r, j);
     195    21106079 :     if (!FindApplyQ(r, L, B, j, Q, prec)) return 0;
     196             :   }
     197     4439893 :   *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      244148 : gaussred_from_QR(GEN x, long prec)
     203             : {
     204      244148 :   long j, k = lg(x)-1;
     205             :   GEN B, Q, L;
     206      244148 :   if (!QR_init(x, &B,&Q,&L, prec)) return NULL;
     207      976285 :   for (j=1; j<k; j++)
     208             :   {
     209      732134 :     GEN m = gel(L,j), invNx = invr(gel(m,j));
     210             :     long i;
     211      732117 :     gel(m,j) = gel(B,j);
     212     2922258 :     for (i=j+1; i<=k; i++) gel(m,i) = mpmul(invNx, gel(m,i));
     213             :   }
     214      244151 :   gcoeff(L,j,j) = gel(B,j);
     215      244151 :   return shallowtrans(L);
     216             : }
     217             : GEN
     218       14056 : R_from_QR(GEN x, long prec)
     219             : {
     220             :   GEN B, Q, L;
     221       14056 :   if (!QR_init(x, &B,&Q,&L, prec)) return NULL;
     222       14042 :   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       47498 : RgM_gram_schmidt(GEN e, GEN *ptB)
     233             : {
     234       47498 :   long i,j,lx = lg(e);
     235       47498 :   GEN f = RgM_shallowcopy(e), B, iB;
     236             : 
     237       47498 :   B = cgetg(lx, t_VEC);
     238       47498 :   iB= cgetg(lx, t_VEC);
     239             : 
     240      101606 :   for (i=1;i<lx;i++)
     241             :   {
     242       54108 :     GEN p1 = NULL;
     243       54108 :     pari_sp av = avma;
     244      116089 :     for (j=1; j<i; j++)
     245             :     {
     246       61981 :       GEN mu = gmul(RgV_dotproduct(gel(e,i),gel(f,j)), gel(iB,j));
     247       61981 :       GEN p2 = gmul(mu, gel(f,j));
     248       61981 :       p1 = p1? gadd(p1,p2): p2;
     249             :     }
     250       54108 :     p1 = p1? gerepileupto(av, gsub(gel(e,i), p1)): gel(e,i);
     251       54108 :     gel(f,i) = p1;
     252       54108 :     gel(B,i) = RgV_dotsquare(gel(f,i));
     253       54108 :     gel(iB,i) = ginv(gel(B,i));
     254             :   }
     255       47498 :   *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       47498 : RgM_Babai(GEN B, GEN t)
     262             : {
     263       47498 :   GEN C, N, G = RgM_gram_schmidt(B, &N), b = t;
     264       47498 :   long j, n = lg(B)-1;
     265             : 
     266       47498 :   C = cgetg(n+1,t_COL);
     267      101606 :   for (j = n; j > 0; j--)
     268             :   {
     269       54108 :     GEN c = gdiv( RgV_dotproduct(b, gel(G,j)), gel(N,j) );
     270             :     long e;
     271       54108 :     c = grndtoi(c,&e);
     272       54108 :     if (e >= 0) return NULL;
     273       54108 :     if (signe(c)) b = RgC_sub(b, RgC_Rg_mul(gel(B,j), c));
     274       54108 :     gel(C,j) = c;
     275             :   }
     276       47498 :   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        8669 : lindepfull_bit(GEN x, long bit)
     706             : {
     707        8669 :   long lx = lg(x), ly, i, j;
     708             :   GEN re, im, M;
     709             : 
     710        8669 :   if (! is_vec_t(typ(x))) pari_err_TYPE("lindep2",x);
     711        8669 :   if (lx <= 2)
     712             :   {
     713          21 :     if (lx == 2 && gequal0(x)) return matid(1);
     714          14 :     return NULL;
     715             :   }
     716        8648 :   re = real_i(x);
     717        8648 :   im = imag_i(x);
     718             :   /* independent over R ? */
     719        8648 :   if (lx == 3 && real_indep(re,im,bit)) return NULL;
     720        8634 :   if (gequal0(im)) im = NULL;
     721        8634 :   ly = im? lx+2: lx+1;
     722        8634 :   M = cgetg(lx,t_MAT);
     723       40630 :   for (i=1; i<lx; i++)
     724             :   {
     725       31996 :     GEN c = cgetg(ly,t_COL); gel(M,i) = c;
     726      168480 :     for (j=1; j<lx; j++) gel(c,j) = gen_0;
     727       31996 :     gel(c,i) = gen_1;
     728       31996 :     gel(c,lx)           = gtrunc2n(gel(re,i), bit);
     729       31996 :     if (im) gel(c,lx+1) = gtrunc2n(gel(im,i), bit);
     730             :   }
     731        8634 :   return ZM_lll(M, 0.99, LLL_INPLACE);
     732             : }
     733             : GEN
     734        3181 : lindep_bit(GEN x, long bit)
     735             : {
     736        3181 :   pari_sp av = avma;
     737        3181 :   GEN v, M = lindepfull_bit(x,bit);
     738        3181 :   if (!M) { set_avma(av); return cgetg(1, t_COL); }
     739        3153 :   v = gel(M,1); setlg(v, lg(M));
     740        3153 :   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          14 :   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      107428 : minim_alloc(long n, double ***q, GEN *x, double **y,  double **z, double **v)
    1083             : {
    1084             :   long i, s;
    1085             : 
    1086      107428 :   *x = cgetg(n, t_VECSMALL);
    1087      107428 :   *q = (double**) new_chunk(n);
    1088      107429 :   s = n * sizeof(double);
    1089      107429 :   *y = (double*) stack_malloc_align(s, sizeof(double));
    1090      107429 :   *z = (double*) stack_malloc_align(s, sizeof(double));
    1091      107429 :   *v = (double*) stack_malloc_align(s, sizeof(double));
    1092      505288 :   for (i=1; i<n; i++) (*q)[i] = (double*) stack_malloc_align(s, sizeof(double));
    1093      107429 : }
    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_ZM_apply(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)-1, 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 == 0) return;
    1178          28 :   minim_alloc(n+1, &q, &x, &y, &z, &v);
    1179          98 :   for (j=1; j<=n; j++)
    1180             :   {
    1181          70 :     v[j] = rtodbl(gcoeff(r,j,j));
    1182         133 :     for (i=1; i<j; i++) q[i][j] = rtodbl(gcoeff(r,i,j));
    1183             :   }
    1184             : 
    1185          28 :   if (gequal0(BORNE))
    1186             :   {
    1187             :     double c;
    1188          14 :     p = rtodbl(gcoeff(a,1,1));
    1189          42 :     for (i=2; i<=n; i++) { c = rtodbl(gcoeff(a,i,i)); if (c < p) p = c; }
    1190          14 :     BORNE = roundr(dbltor(p));
    1191             :   }
    1192             :   else
    1193          14 :     p = gtodouble(BORNE);
    1194          28 :   BOUND = p * (1 + eps);
    1195          28 :   if (BOUND > (double)ULONG_MAX || (ulong)BOUND != (ulong)p)
    1196           7 :     pari_err_PREC("forqfvec");
    1197             : 
    1198          21 :   k = n; y[n] = z[n] = 0;
    1199          21 :   x[n] = (long)sqrt(BOUND/v[n]);
    1200          56 :   for(;;x[1]--)
    1201             :   {
    1202             :     do
    1203             :     {
    1204         140 :       if (k>1)
    1205             :       {
    1206          84 :         long l = k-1;
    1207          84 :         z[l] = 0;
    1208         245 :         for (j=k; j<=n; j++) z[l] += q[l][j]*x[j];
    1209          84 :         p = (double)x[k] + z[k];
    1210          84 :         y[l] = y[k] + p*p*v[k];
    1211          84 :         x[l] = (long)floor(sqrt((BOUND-y[l])/v[l])-z[l]);
    1212          84 :         k = l;
    1213             :       }
    1214             :       for(;;)
    1215             :       {
    1216         189 :         p = (double)x[k] + z[k];
    1217         189 :         if (y[k] + p*p*v[k] <= BOUND) break;
    1218          49 :         k++; x[k]--;
    1219             :       }
    1220         140 :     } while (k > 1);
    1221          77 :     if (! x[1] && y[1]<=eps) break;
    1222             : 
    1223          56 :     p = (double)x[1] + z[1]; p = y[1] + p*p*v[1]; /* norm(x) */
    1224          56 :     if (fun(E, u, x, p)) break;
    1225             :   }
    1226             : }
    1227             : 
    1228             : void
    1229           0 : forqfvec(void *E, long (*fun)(void *, GEN, GEN, double), GEN a, GEN BORNE)
    1230             : {
    1231           0 :   pari_sp av = avma;
    1232             :   struct qfvec qv;
    1233           0 :   forqfvec_init(&qv, a);
    1234           0 :   forqfvec_i(E, fun, &qv, BORNE);
    1235           0 :   set_avma(av);
    1236           0 : }
    1237             : 
    1238             : struct qfvecwrap
    1239             : {
    1240             :   void *E;
    1241             :   long (*fun)(void *, GEN);
    1242             : };
    1243             : 
    1244             : static long
    1245          56 : forqfvec_wrap(void *E, GEN u, GEN x, double d)
    1246             : {
    1247          56 :   pari_sp av = avma;
    1248          56 :   struct qfvecwrap *W = (struct qfvecwrap *) E;
    1249             :   (void) d;
    1250          56 :   return gc_long(av, W->fun(W->E, ZM_zc_mul_canon(u, x)));
    1251             : }
    1252             : 
    1253             : void
    1254          42 : forqfvec1(void *E, long (*fun)(void *, GEN), GEN a, GEN BORNE)
    1255             : {
    1256          42 :   pari_sp av = avma;
    1257             :   struct qfvecwrap wr;
    1258             :   struct qfvec qv;
    1259          42 :   wr.E = E; wr.fun = fun;
    1260          42 :   forqfvec_init(&qv, a);
    1261          42 :   forqfvec_i((void*) &wr, forqfvec_wrap, &qv, BORNE);
    1262          35 :   set_avma(av);
    1263          35 : }
    1264             : 
    1265             : void
    1266          42 : forqfvec0(GEN a, GEN BORNE, GEN code)
    1267          42 : { EXPRVOID_WRAP(code, forqfvec1(EXPR_ARGVOID, a,  BORNE)) }
    1268             : 
    1269             : enum { min_ALL = 0, min_FIRST, min_VECSMALL, min_VECSMALL2 };
    1270             : 
    1271             : /* Minimal vectors for the integral definite quadratic form: a.
    1272             :  * Result u:
    1273             :  *   u[1]= Number of vectors of square norm <= BORNE
    1274             :  *   u[2]= maximum norm found
    1275             :  *   u[3]= list of vectors found (at most STOCKMAX, unless NULL)
    1276             :  *
    1277             :  *  If BORNE = NULL: Minimal nonzero vectors.
    1278             :  *  flag = min_ALL,   as above
    1279             :  *  flag = min_FIRST, exits when first suitable vector is found.
    1280             :  *  flag = min_VECSMALL, return a t_VECSMALL of (half) the number of vectors
    1281             :  *  for each norm
    1282             :  *  flag = min_VECSMALL2, same but count only vectors with even norm, and shift
    1283             :  *  the answer */
    1284             : static GEN
    1285         847 : minim0_dolll(GEN a, GEN BORNE, GEN STOCKMAX, long flag, long dolll)
    1286             : {
    1287             :   GEN x, u, r, L, gnorme;
    1288         847 :   long n = lg(a)-1, i, j, k, s, maxrank, sBORNE;
    1289         847 :   pari_sp av = avma, av1;
    1290             :   double p,maxnorm,BOUND,*v,*y,*z,**q;
    1291         847 :   const double eps = 1e-10;
    1292         847 :   int stockall = 0;
    1293             :   struct qfvec qv;
    1294             : 
    1295         847 :   if (!BORNE)
    1296          56 :     sBORNE = 0;
    1297             :   else
    1298             :   {
    1299         791 :     BORNE = gfloor(BORNE);
    1300         791 :     if (typ(BORNE) != t_INT) pari_err_TYPE("minim0",BORNE);
    1301         791 :     if (is_bigint(BORNE)) pari_err_PREC( "qfminim");
    1302         790 :     sBORNE = itos(BORNE); set_avma(av);
    1303         790 :     if (sBORNE < 0) sBORNE = 0;
    1304             :   }
    1305         846 :   if (!STOCKMAX)
    1306             :   {
    1307         335 :     stockall = 1;
    1308         335 :     maxrank = 200;
    1309             :   }
    1310             :   else
    1311             :   {
    1312         511 :     STOCKMAX = gfloor(STOCKMAX);
    1313         511 :     if (typ(STOCKMAX) != t_INT) pari_err_TYPE("minim0",STOCKMAX);
    1314         511 :     maxrank = itos(STOCKMAX);
    1315         511 :     if (maxrank < 0)
    1316           0 :       pari_err_TYPE("minim0 [negative number of vectors]",STOCKMAX);
    1317             :   }
    1318             : 
    1319         846 :   switch(flag)
    1320             :   {
    1321         462 :     case min_VECSMALL:
    1322             :     case min_VECSMALL2:
    1323         462 :       if (sBORNE <= 0) return cgetg(1, t_VECSMALL);
    1324         434 :       L = zero_zv(sBORNE);
    1325         434 :       if (flag == min_VECSMALL2) sBORNE <<= 1;
    1326         434 :       if (n == 0) return L;
    1327         434 :       break;
    1328          35 :     case min_FIRST:
    1329          35 :       if (n == 0 || (!sBORNE && BORNE)) return cgetg(1,t_VEC);
    1330          21 :       L = NULL; /* gcc -Wall */
    1331          21 :       break;
    1332         349 :     case min_ALL:
    1333         349 :       if (n == 0 || (!sBORNE && BORNE))
    1334          14 :         retmkvec3(gen_0, gen_0, cgetg(1, t_MAT));
    1335         335 :       L = new_chunk(1+maxrank);
    1336         335 :       break;
    1337           0 :     default:
    1338           0 :       return NULL;
    1339             :   }
    1340         790 :   minim_alloc(n+1, &q, &x, &y, &z, &v);
    1341             : 
    1342         790 :   forqfvec_init_dolll(&qv, &a, dolll);
    1343         790 :   av1 = avma;
    1344         790 :   r = qv.r;
    1345         790 :   u = qv.u;
    1346        5912 :   for (j=1; j<=n; j++)
    1347             :   {
    1348        5122 :     v[j] = rtodbl(gcoeff(r,j,j));
    1349       29579 :     for (i=1; i<j; i++) q[i][j] = rtodbl(gcoeff(r,i,j)); /* |.| <= 1/2 */
    1350             :   }
    1351             : 
    1352         790 :   if (sBORNE) maxnorm = 0.;
    1353             :   else
    1354             :   {
    1355          56 :     GEN B = gcoeff(a,1,1);
    1356          56 :     long t = 1;
    1357         616 :     for (i=2; i<=n; i++)
    1358             :     {
    1359         560 :       GEN c = gcoeff(a,i,i);
    1360         560 :       if (cmpii(c, B) < 0) { B = c; t = i; }
    1361             :     }
    1362          56 :     if (flag == min_FIRST) return gerepilecopy(av, mkvec2(B, gel(u,t)));
    1363          49 :     maxnorm = -1.; /* don't update maxnorm */
    1364          49 :     if (is_bigint(B)) return NULL;
    1365          48 :     sBORNE = itos(B);
    1366             :   }
    1367         782 :   BOUND = sBORNE * (1 + eps);
    1368         782 :   if ((long)BOUND != sBORNE) return NULL;
    1369             : 
    1370         770 :   s = 0;
    1371         770 :   k = n; y[n] = z[n] = 0;
    1372         770 :   x[n] = (long)sqrt(BOUND/v[n]);
    1373     1223264 :   for(;;x[1]--)
    1374             :   {
    1375             :     do
    1376             :     {
    1377     2245614 :       if (k>1)
    1378             :       {
    1379     1022259 :         long l = k-1;
    1380     1022259 :         z[l] = 0;
    1381    11756080 :         for (j=k; j<=n; j++) z[l] += q[l][j]*x[j];
    1382     1022259 :         p = (double)x[k] + z[k];
    1383     1022259 :         y[l] = y[k] + p*p*v[k];
    1384     1022259 :         x[l] = (long)floor(sqrt((BOUND-y[l])/v[l])-z[l]);
    1385     1022259 :         k = l;
    1386             :       }
    1387             :       for(;;)
    1388             :       {
    1389     3263729 :         p = (double)x[k] + z[k];
    1390     3263729 :         if (y[k] + p*p*v[k] <= BOUND) break;
    1391     1018115 :         k++; x[k]--;
    1392             :       }
    1393             :     }
    1394     2245614 :     while (k > 1);
    1395     1224034 :     if (! x[1] && y[1]<=eps) break;
    1396             : 
    1397     1223271 :     p = (double)x[1] + z[1]; p = y[1] + p*p*v[1]; /* norm(x) */
    1398     1223271 :     if (maxnorm >= 0)
    1399             :     {
    1400     1220723 :       if (p > maxnorm) maxnorm = p;
    1401             :     }
    1402             :     else
    1403             :     { /* maxnorm < 0 : only look for minimal vectors */
    1404        2548 :       pari_sp av2 = avma;
    1405        2548 :       gnorme = roundr(dbltor(p));
    1406        2548 :       if (cmpis(gnorme, sBORNE) >= 0) set_avma(av2);
    1407             :       else
    1408             :       {
    1409          14 :         sBORNE = itos(gnorme); set_avma(av1);
    1410          14 :         BOUND = sBORNE * (1+eps);
    1411          14 :         L = new_chunk(maxrank+1);
    1412          14 :         s = 0;
    1413             :       }
    1414             :     }
    1415     1223271 :     s++;
    1416             : 
    1417     1223271 :     switch(flag)
    1418             :     {
    1419           7 :       case min_FIRST:
    1420           7 :         if (dolll) x = ZM_zc_mul_canon(u,x);
    1421           7 :         return gerepilecopy(av, mkvec2(roundr(dbltor(p)), x));
    1422             : 
    1423      248241 :       case min_ALL:
    1424      248241 :         if (s > maxrank && stockall) /* overflow */
    1425             :         {
    1426         490 :           long maxranknew = maxrank << 1;
    1427         490 :           GEN Lnew = new_chunk(1 + maxranknew);
    1428      344890 :           for (i=1; i<=maxrank; i++) Lnew[i] = L[i];
    1429         490 :           L = Lnew; maxrank = maxranknew;
    1430             :         }
    1431      248241 :         if (s<=maxrank) gel(L,s) = leafcopy(x);
    1432      248241 :         break;
    1433             : 
    1434       39200 :       case min_VECSMALL:
    1435       39200 :         { ulong norm = (ulong)(p + 0.5); L[norm]++; }
    1436       39200 :         break;
    1437             : 
    1438      935823 :       case min_VECSMALL2:
    1439      935823 :         { ulong norm = (ulong)(p + 0.5); if (!odd(norm)) L[norm>>1]++; }
    1440      935823 :         break;
    1441             : 
    1442             :     }
    1443     1223264 :   }
    1444         763 :   switch(flag)
    1445             :   {
    1446           7 :     case min_FIRST:
    1447           7 :       set_avma(av); return cgetg(1,t_VEC);
    1448         434 :     case min_VECSMALL:
    1449             :     case min_VECSMALL2:
    1450         434 :       set_avma((pari_sp)L); return L;
    1451             :   }
    1452         322 :   r = (maxnorm >= 0) ? roundr(dbltor(maxnorm)): stoi(sBORNE);
    1453         322 :   k = minss(s,maxrank);
    1454         322 :   L[0] = evaltyp(t_MAT) | evallg(k + 1);
    1455         322 :   if (dolll)
    1456      246092 :     for (j=1; j<=k; j++)
    1457      245805 :       gel(L,j) = dolll==1 ? ZM_zc_mul_canon(u, gel(L,j))
    1458      245805 :                           : ZM_zc_mul_canon_zm(u, gel(L,j));
    1459         322 :   return gerepilecopy(av, mkvec3(stoi(s<<1), r, L));
    1460             : }
    1461             : 
    1462             : static GEN
    1463         553 : minim0(GEN a, GEN BORNE, GEN STOCKMAX, long flag)
    1464             : {
    1465         553 :   GEN v = minim0_dolll(a, BORNE, STOCKMAX, flag, 1);
    1466         552 :   if (!v) pari_err_PREC("qfminim");
    1467         546 :   return v;
    1468             : }
    1469             : 
    1470             : static GEN
    1471         252 : minim0_zm(GEN a, GEN BORNE, GEN STOCKMAX, long flag)
    1472             : {
    1473         252 :   GEN v = minim0_dolll(a, BORNE, STOCKMAX, flag, 2);
    1474         252 :   if (!v) pari_err_PREC("qfminim");
    1475         252 :   return v;
    1476             : }
    1477             : 
    1478             : GEN
    1479         462 : qfrep0(GEN a, GEN borne, long flag)
    1480         462 : { return minim0(a, borne, gen_0, (flag & 1)? min_VECSMALL2: min_VECSMALL); }
    1481             : 
    1482             : GEN
    1483         133 : qfminim0(GEN a, GEN borne, GEN stockmax, long flag, long prec)
    1484             : {
    1485         133 :   switch(flag)
    1486             :   {
    1487          49 :     case 0: return minim0(a,borne,stockmax,min_ALL);
    1488          35 :     case 1: return minim0(a,borne,gen_0   ,min_FIRST);
    1489          49 :     case 2:
    1490             :     {
    1491          49 :       long maxnum = -1;
    1492          49 :       if (typ(a) != t_MAT) pari_err_TYPE("qfminim",a);
    1493          49 :       if (stockmax) {
    1494          14 :         if (typ(stockmax) != t_INT) pari_err_TYPE("qfminim",stockmax);
    1495          14 :         maxnum = itos(stockmax);
    1496             :       }
    1497          49 :       a = fincke_pohst(a,borne,maxnum,prec,NULL);
    1498          42 :       if (!a) pari_err_PREC("qfminim");
    1499          42 :       return a;
    1500             :     }
    1501           0 :     default: pari_err_FLAG("qfminim");
    1502             :   }
    1503             :   return NULL; /* LCOV_EXCL_LINE */
    1504             : }
    1505             : 
    1506             : GEN
    1507           7 : minim(GEN a, GEN borne, GEN stockmax)
    1508           7 : { return minim0(a,borne,stockmax,min_ALL); }
    1509             : 
    1510             : GEN
    1511         252 : minim_zm(GEN a, GEN borne, GEN stockmax)
    1512         252 : { return minim0_zm(a,borne,stockmax,min_ALL); }
    1513             : 
    1514             : GEN
    1515          42 : minim_raw(GEN a, GEN BORNE, GEN STOCKMAX)
    1516          42 : { return minim0_dolll(a, BORNE, STOCKMAX, min_ALL, 0); }
    1517             : 
    1518             : GEN
    1519           0 : minim2(GEN a, GEN borne, GEN stockmax)
    1520           0 : { return minim0(a,borne,stockmax,min_FIRST); }
    1521             : 
    1522             : /* If V depends linearly from the columns of the matrix, return 0.
    1523             :  * Otherwise, update INVP and L and return 1. No GC. */
    1524             : static int
    1525        1652 : addcolumntomatrix(GEN V, GEN invp, GEN L)
    1526             : {
    1527        1652 :   long i,j,k, n = lg(invp);
    1528        1652 :   GEN a = cgetg(n, t_COL), ak = NULL, mak;
    1529             : 
    1530       84231 :   for (k = 1; k < n; k++)
    1531       83706 :     if (!L[k])
    1532             :     {
    1533       27902 :       ak = RgMrow_zc_mul(invp, V, k);
    1534       27902 :       if (!gequal0(ak)) break;
    1535             :     }
    1536        1652 :   if (k == n) return 0;
    1537        1127 :   L[k] = 1;
    1538        1127 :   mak = gneg_i(ak);
    1539       43253 :   for (i=k+1; i<n; i++)
    1540       42126 :     gel(a,i) = gdiv(RgMrow_zc_mul(invp, V, i), mak);
    1541       43883 :   for (j=1; j<=k; j++)
    1542             :   {
    1543       42756 :     GEN c = gel(invp,j), ck = gel(c,k);
    1544       42756 :     if (gequal0(ck)) continue;
    1545        8757 :     gel(c,k) = gdiv(ck, ak);
    1546        8757 :     if (j==k)
    1547       43253 :       for (i=k+1; i<n; i++)
    1548       42126 :         gel(c,i) = gmul(gel(a,i), ck);
    1549             :     else
    1550      184814 :       for (i=k+1; i<n; i++)
    1551      177184 :         gel(c,i) = gadd(gel(c,i), gmul(gel(a,i), ck));
    1552             :   }
    1553        1127 :   return 1;
    1554             : }
    1555             : 
    1556             : GEN
    1557          42 : qfperfection(GEN a)
    1558             : {
    1559          42 :   pari_sp av = avma;
    1560             :   GEN u, L;
    1561          42 :   long r, s, k, l, n = lg(a)-1;
    1562             : 
    1563          42 :   if (!n) return gen_0;
    1564          42 :   if (typ(a) != t_MAT || !RgM_is_ZM(a)) pari_err_TYPE("qfperfection",a);
    1565          42 :   a = minim_lll(a, &u);
    1566          42 :   L = minim_raw(a,NULL,NULL);
    1567          42 :   r = (n*(n+1)) >> 1;
    1568          42 :   if (L)
    1569             :   {
    1570             :     GEN D, V, invp;
    1571          35 :     L = gel(L, 3); l = lg(L);
    1572          35 :     if (l == 2) { set_avma(av); return gen_1; }
    1573             :     /* |L[i]|^2 fits  into a long for all i */
    1574          21 :     D = zero_zv(r);
    1575          21 :     V = cgetg(r+1, t_VECSMALL);
    1576          21 :     invp = matid(r);
    1577          21 :     s = 0;
    1578        1659 :     for (k = 1; k < l; k++)
    1579             :     {
    1580        1652 :       pari_sp av2 = avma;
    1581        1652 :       GEN x = gel(L,k);
    1582             :       long i, j, I;
    1583       21098 :       for (i = I = 1; i<=n; i++)
    1584      145278 :         for (j=i; j<=n; j++,I++) V[I] = x[i]*x[j];
    1585        1652 :       if (!addcolumntomatrix(V,invp,D)) set_avma(av2);
    1586        1127 :       else if (++s == r) break;
    1587             :     }
    1588             :   }
    1589             :   else
    1590             :   {
    1591             :     GEN M;
    1592           7 :     L = fincke_pohst(a,NULL,-1, DEFAULTPREC, NULL);
    1593           7 :     if (!L) pari_err_PREC("qfminim");
    1594           7 :     L = gel(L, 3); l = lg(L);
    1595           7 :     if (l == 2) { set_avma(av); return gen_1; }
    1596           7 :     M = cgetg(l, t_MAT);
    1597         959 :     for (k = 1; k < l; k++)
    1598             :     {
    1599         952 :       GEN x = gel(L,k), c = cgetg(r+1, t_COL);
    1600             :       long i, I, j;
    1601       16184 :       for (i = I = 1; i<=n; i++)
    1602      144704 :         for (j=i; j<=n; j++,I++) gel(c,I) = mulii(gel(x,i), gel(x,j));
    1603         952 :       gel(M,k) = c;
    1604             :     }
    1605           7 :     s = ZM_rank(M);
    1606             :   }
    1607          28 :   return gc_utoipos(av, s);
    1608             : }
    1609             : 
    1610             : static GEN
    1611         159 : clonefill(GEN S, long s, long t)
    1612             : { /* initialize to dummy values */
    1613         159 :   GEN T = S, dummy = cgetg(1, t_STR);
    1614             :   long i;
    1615      324485 :   for (i = s+1; i <= t; i++) gel(S,i) = dummy;
    1616         159 :   S = gclone(S); if (isclone(T)) gunclone(T);
    1617         159 :   return S;
    1618             : }
    1619             : 
    1620             : /* increment ZV x, by incrementing cell of index k. Initial value x0[k] was
    1621             :  * chosen to minimize qf(x) for given x0[1..k-1] and x0[k+1,..] = 0
    1622             :  * The last nonzero entry must be positive and goes through x0[k]+1,2,3,...
    1623             :  * Others entries go through: x0[k]+1,-1,2,-2,...*/
    1624             : INLINE void
    1625     4271178 : step(GEN x, GEN y, GEN inc, long k)
    1626             : {
    1627     4271178 :   if (!signe(gel(y,k))) /* x[k+1..] = 0 */
    1628      156649 :     gel(x,k) = addiu(gel(x,k), 1); /* leading coeff > 0 */
    1629             :   else
    1630             :   {
    1631     4114529 :     long i = inc[k];
    1632     4114529 :     gel(x,k) = addis(gel(x,k), i),
    1633     4114548 :     inc[k] = (i > 0)? -1-i: 1-i;
    1634             :   }
    1635     4271197 : }
    1636             : 
    1637             : /* 1 if we are "sure" that x < y, up to few rounding errors, i.e.
    1638             :  * x < y - epsilon. More precisely :
    1639             :  * - sign(x - y) < 0
    1640             :  * - lgprec(x-y) > 3 || expo(x - y) - expo(x) > -24 */
    1641             : static int
    1642     1848704 : mplessthan(GEN x, GEN y)
    1643             : {
    1644     1848704 :   pari_sp av = avma;
    1645     1848704 :   GEN z = mpsub(x, y);
    1646     1848703 :   set_avma(av);
    1647     1848703 :   if (typ(z) == t_INT) return (signe(z) < 0);
    1648     1848703 :   if (signe(z) >= 0) return 0;
    1649       20326 :   if (realprec(z) > LOWDEFAULTPREC) return 1;
    1650       20326 :   return ( expo(z) - mpexpo(x) > -24 );
    1651             : }
    1652             : 
    1653             : /* 1 if we are "sure" that x > y, up to few rounding errors, i.e.
    1654             :  * x > y + epsilon */
    1655             : static int
    1656     6568050 : mpgreaterthan(GEN x, GEN y)
    1657             : {
    1658     6568050 :   pari_sp av = avma;
    1659     6568050 :   GEN z = mpsub(x, y);
    1660     6568043 :   set_avma(av);
    1661     6568094 :   if (typ(z) == t_INT) return (signe(z) > 0);
    1662     6568094 :   if (signe(z) <= 0) return 0;
    1663     4005980 :   if (realprec(z) > LOWDEFAULTPREC) return 1;
    1664      468926 :   return ( expo(z) - mpexpo(x) > -24 );
    1665             : }
    1666             : 
    1667             : /* x a t_INT, y  t_INT or t_REAL */
    1668             : INLINE GEN
    1669     1860565 : mulimp(GEN x, GEN y)
    1670             : {
    1671     1860565 :   if (typ(y) == t_INT) return mulii(x,y);
    1672     1860565 :   return signe(x) ? mulir(x,y): gen_0;
    1673             : }
    1674             : /* x + y*z, x,z two mp's, y a t_INT */
    1675             : INLINE GEN
    1676    22470482 : addmulimp(GEN x, GEN y, GEN z)
    1677             : {
    1678    22470482 :   if (!signe(y)) return x;
    1679     9691331 :   if (typ(z) == t_INT) return mpadd(x, mulii(y, z));
    1680     9691331 :   return mpadd(x, mulir(y, z));
    1681             : }
    1682             : 
    1683             : /* yk + vk * (xk + zk)^2 */
    1684             : static GEN
    1685     8359682 : norm_aux(GEN xk, GEN yk, GEN zk, GEN vk)
    1686             : {
    1687     8359682 :   GEN t = mpadd(xk, zk);
    1688     8359701 :   if (typ(t) == t_INT) { /* probably gen_0, avoid loss of accuracy */
    1689      298217 :     yk = addmulimp(yk, sqri(t), vk);
    1690             :   } else {
    1691     8061484 :     yk = mpadd(yk, mpmul(sqrr(t), vk));
    1692             :   }
    1693     8359632 :   return yk;
    1694             : }
    1695             : /* yk + vk * (xk + zk)^2 < B + epsilon */
    1696             : static int
    1697     6117975 : check_bound(GEN B, GEN xk, GEN yk, GEN zk, GEN vk)
    1698             : {
    1699     6117975 :   pari_sp av = avma;
    1700     6117975 :   int f = mpgreaterthan(norm_aux(xk,yk,zk,vk), B);
    1701     6117959 :   return gc_bool(av, !f);
    1702             : }
    1703             : 
    1704             : /* q(k-th canonical basis vector), where q is given in Cholesky form
    1705             :  * q(x) = sum_{i = 1}^n q[i,i] (x[i] + sum_{j > i} q[i,j] x[j])^2.
    1706             :  * Namely q(e_k) = q[k,k] + sum_{i < k} q[i,i] q[i,k]^2
    1707             :  * Assume 1 <= k <= n. */
    1708             : static GEN
    1709         182 : cholesky_norm_ek(GEN q, long k)
    1710             : {
    1711         182 :   GEN t = gcoeff(q,k,k);
    1712             :   long i;
    1713        1484 :   for (i = 1; i < k; i++) t = norm_aux(gen_0, t, gcoeff(q,i,k), gcoeff(q,i,i));
    1714         182 :   return t;
    1715             : }
    1716             : 
    1717             : /* q is the Cholesky decomposition of a quadratic form
    1718             :  * Enumerate vectors whose norm is less than BORNE (Algo 2.5.7),
    1719             :  * minimal vectors if BORNE = NULL (implies check = NULL).
    1720             :  * If (check != NULL) consider only vectors passing the check, and assumes
    1721             :  *   we only want the smallest possible vectors */
    1722             : static GEN
    1723       14438 : smallvectors(GEN q, GEN BORNE, long maxnum, FP_chk_fun *CHECK)
    1724             : {
    1725       14438 :   long N = lg(q), n = N-1, i, j, k, s, stockmax, checkcnt = 1;
    1726             :   pari_sp av, av1;
    1727             :   GEN inc, S, x, y, z, v, p1, alpha, norms;
    1728             :   GEN norme1, normax1, borne1, borne2;
    1729       14438 :   GEN (*check)(void *,GEN) = CHECK? CHECK->f: NULL;
    1730       14438 :   void *data = CHECK? CHECK->data: NULL;
    1731       14438 :   const long skipfirst = CHECK? CHECK->skipfirst: 0;
    1732       14438 :   const int stockall = (maxnum == -1);
    1733             : 
    1734       14438 :   alpha = dbltor(0.95);
    1735       14438 :   normax1 = gen_0;
    1736             : 
    1737       14438 :   v = cgetg(N,t_VEC);
    1738       14438 :   inc = const_vecsmall(n, 1);
    1739             : 
    1740       14438 :   av = avma;
    1741       14438 :   stockmax = stockall? 2000: maxnum;
    1742       14438 :   norms = cgetg(check?(stockmax+1): 1,t_VEC); /* unused if (!check) */
    1743       14438 :   S = cgetg(stockmax+1,t_VEC);
    1744       14438 :   x = cgetg(N,t_COL);
    1745       14438 :   y = cgetg(N,t_COL);
    1746       14438 :   z = cgetg(N,t_COL);
    1747       95685 :   for (i=1; i<N; i++) {
    1748       81247 :     gel(v,i) = gcoeff(q,i,i);
    1749       81247 :     gel(x,i) = gel(y,i) = gel(z,i) = gen_0;
    1750             :   }
    1751       14438 :   if (BORNE)
    1752             :   {
    1753       14417 :     borne1 = BORNE;
    1754       14417 :     if (gsigne(borne1) <= 0) retmkvec3(gen_0, gen_0, cgetg(1,t_MAT));
    1755       14403 :     if (typ(borne1) != t_REAL)
    1756             :     {
    1757             :       long prec;
    1758         368 :       prec = nbits2prec(gexpo(borne1) + 10);
    1759         368 :       borne1 = gtofp(borne1, maxss(prec, DEFAULTPREC));
    1760             :     }
    1761             :   }
    1762             :   else
    1763             :   {
    1764          21 :     borne1 = gcoeff(q,1,1);
    1765         203 :     for (i=2; i<N; i++)
    1766             :     {
    1767         182 :       GEN b = cholesky_norm_ek(q, i);
    1768         182 :       if (gcmp(b, borne1) < 0) borne1 = b;
    1769             :     }
    1770             :     /* borne1 = norm of smallest basis vector */
    1771             :   }
    1772       14424 :   borne2 = mulrr(borne1,alpha);
    1773       14424 :   if (DEBUGLEVEL>2)
    1774           0 :     err_printf("smallvectors looking for norm < %P.4G\n",borne1);
    1775       14424 :   s = 0; k = n;
    1776      380542 :   for(;; step(x,y,inc,k)) /* main */
    1777             :   { /* x (supposedly) small vector, ZV.
    1778             :      * For all t >= k, we have
    1779             :      *   z[t] = sum_{j > t} q[t,j] * x[j]
    1780             :      *   y[t] = sum_{i > t} q[i,i] * (x[i] + z[i])^2
    1781             :      *        = 0 <=> x[i]=0 for all i>t */
    1782             :     do
    1783             :     {
    1784     2241099 :       int skip = 0;
    1785     2241099 :       if (k > 1)
    1786             :       {
    1787     1860566 :         long l = k-1;
    1788     1860566 :         av1 = avma;
    1789     1860566 :         p1 = mulimp(gel(x,k), gcoeff(q,l,k));
    1790    24032868 :         for (j=k+1; j<N; j++) p1 = addmulimp(p1, gel(x,j), gcoeff(q,l,j));
    1791     1860578 :         gel(z,l) = gerepileuptoleaf(av1,p1);
    1792             : 
    1793     1860571 :         av1 = avma;
    1794     1860571 :         p1 = norm_aux(gel(x,k), gel(y,k), gel(z,k), gel(v,k));
    1795     1860568 :         gel(y,l) = gerepileuptoleaf(av1, p1);
    1796             :         /* skip the [x_1,...,x_skipfirst,0,...,0] */
    1797     1860569 :         if ((l <= skipfirst && !signe(gel(y,skipfirst)))
    1798     1860569 :          || mplessthan(borne1, gel(y,l))) skip = 1;
    1799             :         else /* initial value, minimizing (x[l] + z[l])^2, hence qf(x) for
    1800             :                 the given x[1..l-1] */
    1801     1846813 :           gel(x,l) = mpround( mpneg(gel(z,l)) );
    1802     1860570 :         k = l;
    1803             :       }
    1804     1860563 :       for(;; step(x,y,inc,k))
    1805             :       { /* at most 2n loops */
    1806     4101667 :         if (!skip)
    1807             :         {
    1808     4087906 :           if (check_bound(borne1, gel(x,k),gel(y,k),gel(z,k),gel(v,k))) break;
    1809     2030086 :           step(x,y,inc,k);
    1810     2030104 :           if (check_bound(borne1, gel(x,k),gel(y,k),gel(z,k),gel(v,k))) break;
    1811             :         }
    1812     1874987 :         skip = 0; inc[k] = 1;
    1813     1874987 :         if (++k > n) goto END;
    1814             :       }
    1815             : 
    1816     2226686 :       if (gc_needed(av,2))
    1817             :       {
    1818          33 :         if(DEBUGMEM>1) pari_warn(warnmem,"smallvectors");
    1819          33 :         if (stockmax) S = clonefill(S, s, stockmax);
    1820          33 :         if (check) {
    1821          33 :           GEN dummy = cgetg(1, t_STR);
    1822       23191 :           for (i=s+1; i<=stockmax; i++) gel(norms,i) = dummy;
    1823             :         }
    1824          33 :         gerepileall(av,7,&x,&y,&z,&normax1,&borne1,&borne2,&norms);
    1825             :       }
    1826             :     }
    1827     2226686 :     while (k > 1);
    1828      380543 :     if (!signe(gel(x,1)) && !signe(gel(y,1))) continue; /* exclude 0 */
    1829             : 
    1830      379875 :     av1 = avma;
    1831      379875 :     norme1 = norm_aux(gel(x,1),gel(y,1),gel(z,1),gel(v,1));
    1832      379872 :     if (mpgreaterthan(norme1,borne1)) { set_avma(av1); continue; /* main */ }
    1833             : 
    1834      379874 :     norme1 = gerepileuptoleaf(av1,norme1);
    1835      379874 :     if (check)
    1836             :     {
    1837      311288 :       if (checkcnt < 5 && mpcmp(norme1, borne2) < 0)
    1838             :       {
    1839        4222 :         if (!check(data,x)) { checkcnt++ ; continue; /* main */}
    1840         412 :         if (DEBUGLEVEL>4) err_printf("New bound: %Ps", norme1);
    1841         412 :         borne1 = norme1;
    1842         412 :         borne2 = mulrr(borne1, alpha);
    1843         412 :         s = 0; checkcnt = 0;
    1844             :       }
    1845             :     }
    1846             :     else
    1847             :     {
    1848       68586 :       if (!BORNE) /* find minimal vectors */
    1849             :       {
    1850        1890 :         if (mplessthan(norme1, borne1))
    1851             :         { /* strictly smaller vector than previously known */
    1852           0 :           borne1 = norme1; /* + epsilon */
    1853           0 :           s = 0;
    1854             :         }
    1855             :       }
    1856             :       else
    1857       66696 :         if (mpcmp(norme1,normax1) > 0) normax1 = norme1;
    1858             :     }
    1859      376064 :     if (++s > stockmax) continue; /* too many vectors: no longer remember */
    1860      375133 :     if (check) gel(norms,s) = norme1;
    1861      375133 :     gel(S,s) = leafcopy(x);
    1862      375133 :     if (s != stockmax) continue; /* still room, get next vector */
    1863             : 
    1864         126 :     if (check)
    1865             :     { /* overflow, eliminate vectors failing "check" */
    1866         105 :       pari_sp av2 = avma;
    1867             :       long imin, imax;
    1868         105 :       GEN per = indexsort(norms), S2 = cgetg(stockmax+1, t_VEC);
    1869         105 :       if (DEBUGLEVEL>2) err_printf("sorting... [%ld elts]\n",s);
    1870             :       /* let N be the minimal norm so far for x satisfying 'check'. Keep
    1871             :        * all elements of norm N */
    1872       26593 :       for (i = 1; i <= s; i++)
    1873             :       {
    1874       26586 :         long k = per[i];
    1875       26586 :         if (check(data,gel(S,k))) { borne1 = gel(norms,k); break; }
    1876             :       }
    1877         105 :       imin = i;
    1878       20937 :       for (; i <= s; i++)
    1879       20916 :         if (mpgreaterthan(gel(norms,per[i]), borne1)) break;
    1880         105 :       imax = i;
    1881       20937 :       for (i=imin, s=0; i < imax; i++) gel(S2,++s) = gel(S,per[i]);
    1882       20937 :       for (i = 1; i <= s; i++) gel(S,i) = gel(S2,i);
    1883         105 :       set_avma(av2);
    1884         105 :       if (s) { borne2 = mulrr(borne1, alpha); checkcnt = 0; }
    1885         105 :       if (!stockall) continue;
    1886         105 :       if (s > stockmax/2) stockmax <<= 1;
    1887         105 :       norms = cgetg(stockmax+1, t_VEC);
    1888       20937 :       for (i = 1; i <= s; i++) gel(norms,i) = borne1;
    1889             :     }
    1890             :     else
    1891             :     {
    1892          21 :       if (!stockall && BORNE) goto END;
    1893          21 :       if (!stockall) continue;
    1894          21 :       stockmax <<= 1;
    1895             :     }
    1896             : 
    1897             :     {
    1898         126 :       GEN Snew = clonefill(vec_lengthen(S,stockmax), s, stockmax);
    1899         126 :       if (isclone(S)) gunclone(S);
    1900         126 :       S = Snew;
    1901             :     }
    1902             :   }
    1903       14424 : END:
    1904       14424 :   if (s < stockmax) stockmax = s;
    1905       14424 :   if (check)
    1906             :   {
    1907             :     GEN per, alph, pols, p;
    1908       14396 :     if (DEBUGLEVEL>2) err_printf("final sort & check...\n");
    1909       14396 :     setlg(norms,stockmax+1); per = indexsort(norms);
    1910       14396 :     alph = cgetg(stockmax+1,t_VEC);
    1911       14396 :     pols = cgetg(stockmax+1,t_VEC);
    1912       81955 :     for (j=0,i=1; i<=stockmax; i++)
    1913             :     {
    1914       67798 :       long t = per[i];
    1915       67798 :       GEN N = gel(norms,t);
    1916       67798 :       if (j && mpgreaterthan(N, borne1)) break;
    1917       67559 :       if ((p = check(data,gel(S,t))))
    1918             :       {
    1919       55024 :         if (!j) borne1 = N;
    1920       55024 :         j++;
    1921       55024 :         gel(pols,j) = p;
    1922       55024 :         gel(alph,j) = gel(S,t);
    1923             :       }
    1924             :     }
    1925       14396 :     setlg(pols,j+1);
    1926       14396 :     setlg(alph,j+1);
    1927       14396 :     if (stockmax && isclone(S)) { alph = gcopy(alph); gunclone(S); }
    1928       14396 :     return mkvec2(pols, alph);
    1929             :   }
    1930          28 :   if (stockmax)
    1931             :   {
    1932          21 :     setlg(S,stockmax+1);
    1933          21 :     settyp(S,t_MAT);
    1934          21 :     if (isclone(S)) { p1 = S; S = gcopy(S); gunclone(p1); }
    1935             :   }
    1936             :   else
    1937           7 :     S = cgetg(1,t_MAT);
    1938          28 :   return mkvec3(utoi(s<<1), borne1, S);
    1939             : }
    1940             : 
    1941             : /* solve q(x) = x~.a.x <= bound, a > 0.
    1942             :  * If check is non-NULL keep x only if check(x).
    1943             :  * If a is a vector, assume a[1] is the LLL-reduced Cholesky form of q */
    1944             : GEN
    1945       14459 : fincke_pohst(GEN a, GEN B0, long stockmax, long PREC, FP_chk_fun *CHECK)
    1946             : {
    1947       14459 :   pari_sp av = avma;
    1948             :   VOLATILE long i,j,l;
    1949       14459 :   VOLATILE GEN r,rinv,rinvtrans,u,v,res,z,vnorm,rperm,perm,uperm, bound = B0;
    1950             : 
    1951       14459 :   if (typ(a) == t_VEC)
    1952             :   {
    1953       14042 :     r = gel(a,1);
    1954       14042 :     u = NULL;
    1955             :   }
    1956             :   else
    1957             :   {
    1958         417 :     long prec = PREC;
    1959         417 :     l = lg(a);
    1960         417 :     if (l == 1)
    1961             :     {
    1962           7 :       if (CHECK) pari_err_TYPE("fincke_pohst [dimension 0]", a);
    1963           7 :       retmkvec3(gen_0, gen_0, cgetg(1,t_MAT));
    1964             :     }
    1965         410 :     u = lllfp(a, 0.75, LLL_GRAM | LLL_IM);
    1966         403 :     if (lg(u) != lg(a)) return gc_NULL(av);
    1967         403 :     r = qf_RgM_apply(a,u);
    1968         403 :     i = gprecision(r);
    1969         403 :     if (i)
    1970         361 :       prec = i;
    1971             :     else {
    1972          42 :       prec = DEFAULTPREC + nbits2extraprec(gexpo(r));
    1973          42 :       if (prec < PREC) prec = PREC;
    1974             :     }
    1975         403 :     if (DEBUGLEVEL>2) err_printf("first LLL: prec = %ld\n", prec);
    1976         403 :     r = qfgaussred_positive(r);
    1977         403 :     if (!r) return gc_NULL(av);
    1978        1780 :     for (i=1; i<l; i++)
    1979             :     {
    1980        1377 :       GEN s = gsqrt(gcoeff(r,i,i), prec);
    1981        1377 :       gcoeff(r,i,i) = s;
    1982        3930 :       for (j=i+1; j<l; j++) gcoeff(r,i,j) = gmul(s, gcoeff(r,i,j));
    1983             :     }
    1984             :   }
    1985             :   /* now r~ * r = a in LLL basis */
    1986       14445 :   rinv = RgM_inv_upper(r);
    1987       14445 :   if (!rinv) return gc_NULL(av);
    1988       14445 :   rinvtrans = shallowtrans(rinv);
    1989       14445 :   if (DEBUGLEVEL>2)
    1990           0 :     err_printf("Fincke-Pohst, final LLL: prec = %ld\n", gprecision(rinvtrans));
    1991       14445 :   v = lll(rinvtrans);
    1992       14445 :   if (lg(v) != lg(rinvtrans)) return gc_NULL(av);
    1993             : 
    1994       14445 :   rinvtrans = RgM_mul(rinvtrans, v);
    1995       14445 :   v = ZM_inv(shallowtrans(v),NULL);
    1996       14445 :   r = RgM_mul(r,v);
    1997       14445 :   u = u? ZM_mul(u,v): v;
    1998             : 
    1999       14445 :   l = lg(r);
    2000       14445 :   vnorm = cgetg(l,t_VEC);
    2001       95719 :   for (j=1; j<l; j++) gel(vnorm,j) = gnorml2(gel(rinvtrans,j));
    2002       14444 :   rperm = cgetg(l,t_MAT);
    2003       14444 :   uperm = cgetg(l,t_MAT); perm = indexsort(vnorm);
    2004       95715 :   for (i=1; i<l; i++) { uperm[l-i] = u[perm[i]]; rperm[l-i] = r[perm[i]]; }
    2005       14444 :   u = uperm;
    2006       14444 :   r = rperm; res = NULL;
    2007       14444 :   pari_CATCH(e_PREC) { }
    2008             :   pari_TRY {
    2009             :     GEN q;
    2010       14445 :     if (CHECK && CHECK->f_init) bound = CHECK->f_init(CHECK, r, u);
    2011       14438 :     q = gaussred_from_QR(r, gprecision(vnorm));
    2012       14438 :     if (q) res = smallvectors(q, bound, stockmax, CHECK);
    2013       14438 :   } pari_ENDCATCH;
    2014       14445 :   if (!res) return gc_NULL(av);
    2015       14438 :   if (CHECK)
    2016             :   {
    2017       14396 :     if (CHECK->f_post) res = CHECK->f_post(CHECK, res, u);
    2018       14396 :     return res;
    2019             :   }
    2020             : 
    2021          42 :   z = cgetg(4,t_VEC);
    2022          42 :   gel(z,1) = gcopy(gel(res,1));
    2023          42 :   gel(z,2) = gcopy(gel(res,2));
    2024          42 :   gel(z,3) = ZM_mul(u, gel(res,3)); return gerepileupto(av,z);
    2025             : }

Generated by: LCOV version 1.14