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

Generated by: LCOV version 1.13