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 16624-25b9976) Lines: 929 976 95.2 %
Date: 2014-06-24 Functions: 55 59 93.2 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 637 791 80.5 %

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

Generated by: LCOV version 1.9