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-bordeaux1.fr machine (x86_64 architecture), and agregate them in the final report:

The target is 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.8.0 lcov report (development 16912-212c0f0) Lines: 940 988 95.1 %
Date: 2014-10-20 Functions: 54 60 90.0 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 642 793 81.0 %

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

Generated by: LCOV version 1.9