Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is to exceed 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - basemath - bibli1.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.16.2 lcov report (development 29401-b333f5c7d5) Lines: 1190 1262 94.3 %
Date: 2024-06-19 09:03:24 Functions: 74 80 92.5 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.16