"Fossies" - the Fresh Open Source Software Archive

Member "pari-2.13.1/src/basemath/bibli1.c" (12 Dec 2020, 55776 Bytes) of package /linux/misc/pari-2.13.1.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "bibli1.c" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 2.13.0_vs_2.13.1.

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