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 - modules - thue.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 17110-9967e23) Lines: 802 856 93.7 %
Date: 2014-11-26 Functions: 54 55 98.2 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 481 591 81.4 %

           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                 :            : #include "pari.h"
      15                 :            : #include "paripriv.h"
      16                 :            : 
      17                 :            : /********************************************************************/
      18                 :            : /**                                                                **/
      19                 :            : /**             THUE EQUATION SOLVER (G. Hanrot)                   **/
      20                 :            : /**                                                                **/
      21                 :            : /********************************************************************/
      22                 :            : /* In all the forthcoming remarks, "paper" designs the paper "Thue Equations of
      23                 :            :  * High Degree", by Yu. Bilu and G. Hanrot, J. Number Theory (1996). The
      24                 :            :  * numbering of the constants corresponds to Hanrot's thesis rather than to the
      25                 :            :  * paper. See also
      26                 :            :  * "Solving Thue equations without the full unit group", Math. Comp. (2000) */
      27                 :            : 
      28                 :            : /* Check whether tnf is a valid structure */
      29                 :            : static int
      30                 :        470 : checktnf(GEN tnf)
      31                 :            : {
      32                 :        470 :   long l = lg(tnf);
      33                 :            :   GEN v;
      34 [ +  + ][ +  + ]:        470 :   if (typ(tnf)!=t_VEC || (l!=8 && l!=4)) return 0;
                 [ -  + ]
      35                 :        230 :   v = gel(tnf,1);
      36 [ +  - ][ -  + ]:        230 :   if (typ(v) != t_VEC || lg(v) != 4) return 0;
      37         [ +  + ]:        230 :   if (l == 4) return 1; /* s=0 */
      38                 :            : 
      39                 :        110 :   (void)checkbnf(gel(tnf,2));
      40                 :        470 :   return (typ(gel(tnf,3)) == t_COL
      41         [ +  - ]:        110 :        && typ(gel(tnf,4)) == t_COL
      42         [ +  - ]:        110 :        && typ(gel(tnf,5)) == t_MAT
      43         [ +  - ]:        110 :        && typ(gel(tnf,6)) == t_MAT
      44 [ +  - ][ +  - ]:        220 :        && typ(gel(tnf,7)) == t_VEC);
      45                 :            : }
      46                 :            : 
      47                 :            : /* Compensates rounding errors for computation/display of the constants.
      48                 :            :  * Round up if dir > 0, down otherwise */
      49                 :            : static GEN
      50                 :       8430 : myround(GEN x, long dir)
      51                 :            : {
      52         [ +  + ]:       8430 :   GEN eps = powis(stoi(dir > 0? 10: -10), -10);
      53                 :       8430 :   return gmul(x, gadd(gen_1, eps));
      54                 :            : }
      55                 :            : 
      56                 :            : /* v a t_VEC/t_VEC */
      57                 :            : static GEN
      58                 :       4102 : vecmax_shallow(GEN v) { return gel(v, vecindexmax(v)); }
      59                 :            : 
      60                 :            : static GEN
      61                 :        174 : tnf_get_roots(GEN poly, long prec, long S, long T)
      62                 :            : {
      63                 :        174 :   GEN R0 = QX_complex_roots(poly, prec), R = cgetg(lg(R0), t_COL);
      64                 :            :   long k;
      65                 :            : 
      66         [ +  + ]:        561 :   for (k=1; k<=S; k++) gel(R,k) = gel(R0,k);
      67                 :            :   /* swap roots to get the usual order */
      68         [ +  + ]:        274 :   for (k=1; k<=T; k++)
      69                 :            :   {
      70                 :        100 :     gel(R,k+S)  = gel(R0,2*k+S-1);
      71                 :        100 :     gel(R,k+S+T)= gel(R0,2*k+S);
      72                 :            :   }
      73                 :        174 :   return R;
      74                 :            : }
      75                 :            : 
      76                 :            : /* Computation of the logarithmic height of x (given by embeddings) */
      77                 :            : static GEN
      78                 :       5126 : LogHeight(GEN x, long prec)
      79                 :            : {
      80                 :       5126 :   pari_sp av = avma;
      81                 :       5126 :   long i, n = lg(x)-1;
      82                 :       5126 :   GEN LH = gen_1;
      83         [ +  + ]:      21244 :   for (i=1; i<=n; i++)
      84                 :            :   {
      85                 :      16118 :     GEN t = gabs(gel(x,i), prec);
      86         [ +  + ]:      16118 :     if (gcmpgs(t,1) > 0) LH = gmul(LH, t);
      87                 :            :   }
      88                 :       5126 :   return gerepileupto(av, gdivgs(glog(LH,prec), n));
      89                 :            : }
      90                 :            : 
      91                 :            : /* |x|^(1/n), x t_INT */
      92                 :            : static GEN
      93                 :         60 : absisqrtn(GEN x, long n, long prec)
      94                 :         60 : { GEN r = itor(x,prec); setabssign(r); return sqrtnr(r, n); }
      95                 :            : 
      96                 :            : static GEN
      97                 :       4986 : get_emb(GEN x, GEN r)
      98                 :            : {
      99                 :       4986 :   long l = lg(r), i;
     100                 :            :   GEN y;
     101                 :            : 
     102         [ +  + ]:       4986 :   if (typ(x) == t_INT) return const_col(l-1, x);
     103                 :       4971 :   y = cgetg(l, t_COL);
     104         [ +  + ]:      20548 :   for (i=1; i<l; i++)
     105                 :            :   {
     106                 :      15594 :     GEN e = poleval(x, gel(r,i));
     107 [ +  + ][ +  + ]:      15594 :     if (gequal0(e) || (typ(e) != t_INT && precision(e) <= LOWDEFAULTPREC ))
                 [ +  + ]
     108                 :         17 :       return NULL;
     109                 :      15577 :     gel(y,i) = e;
     110                 :            :   }
     111                 :       4986 :   return y;
     112                 :            : }
     113                 :            : 
     114                 :            : /* Computation of the conjugates (sigma_i(v_j)), and log. heights, of elts of v */
     115                 :            : static GEN
     116                 :        454 : Conj_LH(GEN v, GEN *H, GEN r, long prec)
     117                 :            : {
     118                 :        454 :   long j, l = lg(v);
     119                 :        454 :   GEN e, M = cgetg(l,t_MAT);
     120                 :            : 
     121                 :        454 :   (*H) = cgetg(l,t_COL);
     122         [ +  + ]:       5423 :   for (j = 1; j < l; j++)
     123                 :            :   {
     124         [ +  + ]:       4986 :     if (! (e = get_emb(gel(v,j), r)) ) return NULL; /* FAIL */
     125                 :       4969 :     gel(M,j) = e;
     126                 :       4969 :     gel(*H,j) = LogHeight(e, prec);
     127                 :            :   }
     128                 :        454 :   return M;
     129                 :            : }
     130                 :            : 
     131                 :      11970 : static GEN abslog(GEN x, long prec) { return gabs(glog(x,prec), prec); }
     132                 :        157 : static GEN logabs(GEN x, long prec) { return glog(gabs(x,prec), prec); }
     133                 :            : 
     134                 :            : /* Computation of M, its inverse A and precision check (see paper) */
     135                 :            : static GEN
     136                 :        157 : T_A_Matrices(GEN MatFU, long r, GEN *eps5, long prec)
     137                 :            : {
     138                 :            :   GEN A, p1, m1, IntM, nia, eps3, eps2;
     139                 :        157 :   long e = prec2nbits(prec);
     140                 :            : 
     141                 :        157 :   m1 = rowslice(vecslice(MatFU, 1,r), 1,r); /* minor order r */
     142                 :        157 :   m1 = logabs(m1, 3);
     143                 :            : 
     144         [ -  + ]:        157 :   A = RgM_inv(m1); if (!A) pari_err_PREC("thue");
     145                 :        157 :   IntM = RgM_Rg_add(RgM_mul(A,m1), gen_m1);
     146                 :            : 
     147                 :        157 :   eps2 = gadd(vecmax(gabs(IntM, 3)), real2n(-e, LOWDEFAULTPREC)); /* t_REAL */
     148                 :        157 :   nia = vecmax(gabs(A, 3));
     149         [ -  + ]:        157 :   if (typ(nia) != t_REAL) nia = gtofp(nia, LOWDEFAULTPREC);
     150                 :            : 
     151                 :            :   /* Check for the precision in matrix inversion. See paper, Lemma 2.4.2. */
     152                 :        157 :   p1 = addrr(mulsr(r, gmul2n(nia, e)), eps2); /* t_REAL */
     153         [ -  + ]:        157 :   if (expo(p1) < -2*r) pari_err_PREC("thue");
     154                 :            : 
     155                 :        157 :   p1 = addrr(mulsr(r, gmul2n(nia,-e)), eps2);
     156                 :        157 :   eps3 = mulrr(mulsr(2*r*r,nia), p1);
     157         [ +  + ]:        157 :   if (!signe(eps3))
     158                 :         91 :     eps3 = real2n(expo(eps3), LOWDEFAULTPREC);
     159                 :            :   else
     160                 :         66 :     eps3 = myround(eps3, 1);
     161                 :            : 
     162         [ -  + ]:        157 :   if (DEBUGLEVEL>1) err_printf("epsilon_3 -> %Ps\n",eps3);
     163                 :        157 :   *eps5 = mulur(r, eps3); return A;
     164                 :            : }
     165                 :            : 
     166                 :            : /* Performs basic computations concerning the equation.
     167                 :            :  * Returns a "tnf" structure containing
     168                 :            :  *  1) the polynomial
     169                 :            :  *  2) the bnf (used to solve the norm equation)
     170                 :            :  *  3) roots, with presumably enough precision
     171                 :            :  *  4) The logarithmic heights of units
     172                 :            :  *  5) The matrix of conjugates of units
     173                 :            :  *  6) its inverse
     174                 :            :  *  7) a few technical constants */
     175                 :            : static GEN
     176                 :        157 : inithue(GEN P, GEN bnf, long flag, long prec)
     177                 :            : {
     178                 :        157 :   GEN MatFU, x0, tnf, tmp, gpmin, dP, csts, ALH, eps5, ro, c1, c2, Ind = gen_1;
     179                 :        157 :   long k,j, n = degpol(P);
     180                 :            :   long s,t, prec_roots;
     181                 :            : 
     182         [ +  + ]:        157 :   if (!bnf)
     183                 :            :   {
     184                 :        115 :     bnf = Buchall(P, nf_FORCE, DEFAULTPREC);
     185         [ +  + ]:        115 :     if (flag) (void)bnfcertify(bnf);
     186                 :            :     else
     187                 :        105 :       Ind = floorr(mulru(bnf_get_reg(bnf), 5));
     188                 :            :   }
     189                 :            : 
     190                 :        157 :   nf_get_sign(bnf_get_nf(bnf), &s, &t);
     191                 :        157 :   prec_roots = prec;
     192                 :            :   for(;;)
     193                 :            :   {
     194                 :        174 :     ro = tnf_get_roots(P, prec_roots, s, t);
     195                 :        174 :     MatFU = Conj_LH(bnf_get_fu(bnf), &ALH, ro, prec);
     196         [ +  + ]:        174 :     if (MatFU) break;
     197                 :         17 :     prec_roots = precdbl(prec_roots);
     198         [ -  + ]:         17 :     if (DEBUGLEVEL>1) pari_warn(warnprec, "inithue", prec_roots);
     199                 :         17 :   }
     200                 :            : 
     201                 :        157 :   dP = ZX_deriv(P);
     202                 :        157 :   c1 = NULL; /* min |P'(r_i)|, i <= s */
     203         [ +  + ]:        503 :   for (k=1; k<=s; k++)
     204                 :            :   {
     205                 :        346 :     tmp = gabs(poleval(dP,gel(ro,k)),prec);
     206 [ +  + ][ +  + ]:        346 :     if (!c1 || gcmp(tmp,c1) < 0) c1 = tmp;
     207                 :            :   }
     208                 :        157 :   c1 = gdiv(int2n(n-1), c1);
     209                 :        157 :   c1 = gprec_w(myround(c1, 1), DEFAULTPREC);
     210                 :            : 
     211                 :        157 :   c2 = NULL; /* max |r_i - r_j|, i!=j */
     212         [ +  + ]:        673 :   for (k=1; k<=n; k++)
     213         [ +  + ]:       1167 :     for (j=k+1; j<=n; j++)
     214                 :            :     {
     215                 :        651 :       tmp = gabs(gsub(gel(ro,j),gel(ro,k)), prec);
     216 [ +  + ][ +  + ]:        651 :       if (!c2 || gcmp(c2,tmp) > 0) c2 = tmp;
     217                 :            :     }
     218                 :        157 :   c2 = gprec_w(myround(c2, -1), DEFAULTPREC);
     219                 :            : 
     220         [ +  + ]:        157 :   if (t==0)
     221                 :         87 :     x0 = real_1(DEFAULTPREC);
     222                 :            :   else
     223                 :            :   {
     224                 :         70 :     gpmin = NULL; /* min |P'(r_i)|, i > s */
     225         [ +  + ]:        155 :     for (k=1; k<=t; k++)
     226                 :            :     {
     227                 :         85 :       tmp = gabs(poleval(dP,gel(ro,s+k)), prec);
     228 [ +  + ][ +  + ]:         85 :       if (!gpmin || gcmp(tmp,gpmin) < 0) gpmin = tmp;
     229                 :            :     }
     230                 :         70 :     gpmin = gprec_w(gpmin, DEFAULTPREC);
     231                 :            : 
     232                 :            :     /* Compute x0. See paper, Prop. 2.2.1 */
     233                 :         70 :     x0 = gmul(gpmin, vecmax_shallow(gabs(imag_i(ro), prec)));
     234                 :         70 :     x0 = sqrtnr(gdiv(int2n(n-1), x0), n);
     235                 :            :   }
     236         [ -  + ]:        157 :   if (DEBUGLEVEL>1)
     237                 :          0 :     err_printf("c1 = %Ps\nc2 = %Ps\nIndice <= %Ps\n", c1, c2, Ind);
     238                 :            : 
     239                 :        157 :   ALH = gmul2n(ALH, 1);
     240                 :        157 :   tnf = cgetg(8,t_VEC); csts = cgetg(8,t_VEC);
     241                 :        157 :   gel(tnf,1) = P;
     242                 :        157 :   gel(tnf,2) = bnf;
     243                 :        157 :   gel(tnf,3) = ro;
     244                 :        157 :   gel(tnf,4) = ALH;
     245                 :        157 :   gel(tnf,5) = MatFU;
     246                 :        157 :   gel(tnf,6) = T_A_Matrices(MatFU, s+t-1, &eps5, prec);
     247                 :        157 :   gel(tnf,7) = csts;
     248                 :        157 :   gel(csts,1) = c1; gel(csts,2) = c2;   gel(csts,3) = LogHeight(ro, prec);
     249                 :        157 :   gel(csts,4) = x0; gel(csts,5) = eps5; gel(csts,6) = utoipos(prec);
     250                 :        157 :   gel(csts,7) = Ind;
     251                 :        157 :   return tnf;
     252                 :            : }
     253                 :            : 
     254                 :            : typedef struct {
     255                 :            :   GEN c10, c11, c13, c15, bak, NE, ALH, Ind, hal, MatFU, ro, Hmu;
     256                 :            :   GEN delta, lambda, inverrdelta, Pi, Pi2;
     257                 :            :   long r, iroot, deg;
     258                 :            : } baker_s;
     259                 :            : 
     260                 :            : /* Compute Baker's bound c9 and B_0, the bound for the b_i's. See Thm 2.3.1 */
     261                 :            : static GEN
     262                 :       3885 : Baker(baker_s *BS)
     263                 :            : {
     264                 :       3885 :   const long prec = DEFAULTPREC;
     265                 :       3885 :   GEN tmp, B0, hb0, c9 = gen_1, ro = BS->ro, ro0 = gel(ro,BS->iroot);
     266                 :       3885 :   long k, i1, i2, r = BS->r;
     267                 :            : 
     268      [ +  +  + ]:       3885 :   switch (BS->iroot) {
     269                 :       1883 :     case 1: i1=2; i2=3; break;
     270                 :       1142 :     case 2: i1=1; i2=3; break;
     271                 :        860 :    default: i1=1; i2=2; break;
     272                 :            :   }
     273                 :            : 
     274                 :            :   /* Compute h_1....h_r */
     275         [ +  + ]:      11970 :   for (k=1; k<=r; k++)
     276                 :            :   {
     277                 :       8085 :     tmp = gdiv(gcoeff(BS->MatFU,i1,k), gcoeff(BS->MatFU,i2,k));
     278                 :       8085 :     tmp = gmax(gen_1, abslog(tmp,prec));
     279                 :       8085 :     c9 = gmul(c9, gmax(gel(BS->ALH,k), gdiv(tmp, BS->bak)));
     280                 :            :   }
     281                 :            : 
     282                 :            :   /* Compute a bound for the h_0 */
     283                 :       3885 :   hb0 = gadd(gmul2n(BS->hal,2), gmul2n(gadd(BS->Hmu,mplog2(prec)), 1));
     284                 :       3885 :   tmp = gdiv(gmul(gsub(ro0, gel(ro,i2)), gel(BS->NE,i1)),
     285                 :       7770 :              gmul(gsub(ro0, gel(ro,i1)), gel(BS->NE,i2)));
     286                 :       3885 :   tmp = gmax(gen_1, abslog(tmp, prec));
     287                 :       3885 :   hb0 = gmax(hb0, gdiv(tmp, BS->bak));
     288                 :       3885 :   c9 = gmul(c9,hb0);
     289                 :            :   /* Multiply c9 by the "constant" factor */
     290                 :       3885 :   c9 = gmul(c9, gmul(shiftr(mulur(9,BS->Pi2), 5*(4+r)),
     291                 :       7770 :                      gmul(gmul(mpfact(r+3), powiu(muliu(BS->bak,r+2), r+3)),
     292                 :       3885 :                           glog(muliu(BS->bak,2*(r+2)),prec))));
     293                 :       3885 :   c9 = gprec_w(myround(c9, 1), DEFAULTPREC);
     294                 :            :   /* Compute B0 according to Lemma 2.3.3 */
     295                 :       3885 :   B0 = mulir(shifti(BS->Ind,1),
     296                 :            :              divrr(addrr(mulrr(c9,mplog(divrr(mulir(BS->Ind, c9),BS->c10))),
     297                 :            :                          mplog(mulir(BS->Ind, BS->c11))),
     298                 :            :                    BS->c10));
     299                 :       3885 :   B0 = gmax(B0, dbltor(2.71828183));
     300                 :       3885 :   B0 = gmax(B0, mulrr(divir(BS->Ind, BS->c10),
     301                 :            :                       mplog(divrr(mulir(BS->Ind, BS->c11), BS->Pi2))));
     302                 :            : 
     303         [ -  + ]:       3885 :   if (DEBUGLEVEL>1) {
     304                 :          0 :     err_printf("  B0  = %Ps\n",B0);
     305                 :          0 :     err_printf("  Baker = %Ps\n",c9);
     306                 :            :   }
     307                 :       3885 :   return B0;
     308                 :            : }
     309                 :            : 
     310                 :            : /* || x d ||, x t_REAL, d t_INT */
     311                 :            : static GEN
     312                 :      15958 : errnum(GEN x, GEN d)
     313                 :            : {
     314                 :      15958 :   GEN dx = mulir(d, x), D = subri(dx, roundr(dx));
     315                 :      15958 :   setabssign(D); return D;
     316                 :            : }
     317                 :            : 
     318                 :            : /* Try to reduce the bound through continued fractions; see paper. */
     319                 :            : static int
     320                 :       7939 : CF_1stPass(GEN *B0, GEN kappa, baker_s *BS)
     321                 :            : {
     322                 :       7939 :   GEN a, b, q, ql, qd, l0, denbound = mulri(*B0, kappa);
     323                 :            : 
     324         [ +  + ]:       7939 :   if (cmprr(mulrr(dbltor(0.1),sqrr(denbound)), BS->inverrdelta) > 0)
     325                 :          5 :     return -1;
     326                 :            : 
     327                 :       7934 :   q = denom( bestappr(BS->delta, denbound) );
     328                 :       7934 :   qd = errnum(BS->delta, q);
     329                 :       7934 :   ql = errnum(BS->lambda,q);
     330                 :            : 
     331                 :       7934 :   l0 = subrr(ql, addrr(mulrr(qd, *B0), divri(dbltor(0.1),kappa)));
     332         [ +  + ]:       7934 :   if (signe(l0) <= 0) return 0;
     333                 :            : 
     334         [ +  + ]:       6772 :   if (BS->r > 1) {
     335                 :       6372 :     a = BS->c15; b = BS->c13;
     336                 :            :   }
     337                 :            :   else {
     338                 :        400 :     a = BS->c11; b = BS->c10;
     339                 :        400 :     l0 = mulrr(l0, Pi2n(1, DEFAULTPREC));
     340                 :            :   }
     341                 :       6772 :   *B0 = divrr(mplog(divrr(mulir(q,a), l0)), b);
     342         [ -  + ]:       6772 :   if (DEBUGLEVEL>1) err_printf("    B0 -> %Ps\n",*B0);
     343                 :       7939 :   return 1;
     344                 :            : }
     345                 :            : 
     346                 :            : static void
     347                 :       7144 : get_B0Bx(baker_s *BS, GEN l0, GEN *B0, GEN *Bx)
     348                 :            : {
     349                 :       7144 :   GEN t = divrr(mulir(BS->Ind, BS->c15), l0);
     350                 :       7144 :   *B0 = divrr(mulir(BS->Ind, mplog(t)), BS->c13);
     351                 :       7144 :   *Bx = sqrtnr(shiftr(t,1), BS->deg);
     352                 :       7144 : }
     353                 :            : 
     354                 :            : static int
     355                 :      12689 : LLL_1stPass(GEN *pB0, GEN kappa, baker_s *BS, GEN *pBx)
     356                 :            : {
     357                 :      12689 :   GEN B0 = *pB0, Bx = *pBx, lllmat, C, l0, l1, triv;
     358                 :            :   long e;
     359                 :            : 
     360                 :      12689 :   C = grndtoi(mulir(mulii(BS->Ind, kappa),
     361                 :            :                     gpow(B0, dbltor(2.2), DEFAULTPREC)), &e);
     362                 :            : 
     363         [ -  + ]:      12689 :   if (DEBUGLEVEL > 1) err_printf("C (bitsize) : %d\n", expi(C));
     364                 :      12689 :   lllmat = matid(3);
     365         [ +  + ]:      12689 :   if (cmpri(B0, BS->Ind) > 0)
     366                 :            :   {
     367                 :       8864 :     gcoeff(lllmat, 1, 1) = grndtoi(divri(B0, BS->Ind), &e);
     368                 :       8864 :     triv = shiftr(sqrr(B0), 1);
     369                 :            :   }
     370                 :            :   else
     371                 :       3825 :     triv = addir(sqri(BS->Ind), sqrr(B0));
     372                 :            : 
     373                 :      12689 :   gcoeff(lllmat, 3, 1) = grndtoi(negr(mulir(C, BS->lambda)), &e);
     374         [ +  + ]:      12689 :   if (e >= 0) return -1;
     375                 :      12680 :   gcoeff(lllmat, 3, 2) = grndtoi(negr(mulir(C, BS->delta)), &e);
     376         [ -  + ]:      12680 :   if (e >= 0) return -1;
     377                 :      12680 :   gcoeff(lllmat, 3, 3) = C;
     378                 :      12680 :   lllmat = ZM_lll(lllmat, 0.99, LLL_IM|LLL_INPLACE);
     379                 :            : 
     380                 :      12680 :   l0 = gnorml2(gel(lllmat,1));
     381                 :      12680 :   l0 = subrr(divir(l0, dbltor(1.8262)), triv); /* delta = 0.99 */
     382         [ +  + ]:      12680 :   if (signe(l0) <= 0) return 0;
     383                 :            : 
     384                 :       9000 :   l1 = shiftr(addri(shiftr(B0,1), BS->Ind), -1);
     385                 :       9000 :   l0 = divri(subrr(sqrtr(l0), l1), C);
     386                 :            : 
     387         [ +  + ]:       9000 :   if (signe(l0) <= 0) return 0;
     388                 :            : 
     389                 :       7054 :   get_B0Bx(BS, l0, &B0, &Bx);
     390         [ -  + ]:       7054 :   if (DEBUGLEVEL>=2)
     391                 :            :   {
     392                 :          0 :     err_printf("LLL_First_Pass successful\n");
     393                 :          0 :     err_printf("B0 -> %Ps\n", B0);
     394                 :          0 :     err_printf("x <= %Ps\n", Bx);
     395                 :            :   }
     396                 :      12689 :   *pB0 = B0; *pBx = Bx; return 1;
     397                 :            : }
     398                 :            : 
     399                 :            : /* add solution (x,y) if not already known */
     400                 :            : static void
     401                 :        480 : add_sol(GEN *pS, GEN x, GEN y)
     402                 :        480 : { *pS = shallowconcat(*pS, mkvec(mkvec2(x,y))); }
     403                 :            : /* z = P(p,q), d = deg P, |z| = |rhs|. Check signs and (possibly)
     404                 :            :  * add solutions (p,q), (-p,-q) */
     405                 :            : static void
     406                 :        140 : add_pm(GEN *pS, GEN p, GEN q, GEN z, long d, GEN rhs)
     407                 :            : {
     408         [ +  + ]:        140 :   if (signe(z) == signe(rhs))
     409                 :            :   {
     410                 :        115 :     add_sol(pS, p, q);
     411         [ +  + ]:        115 :     if (!odd(d)) add_sol(pS, negi(p), negi(q));
     412                 :            :   }
     413                 :            :   else
     414         [ +  - ]:         25 :     if (odd(d))  add_sol(pS, negi(p), negi(q));
     415                 :        140 : }
     416                 :            : 
     417                 :            : /* Check whether a potential solution is a true solution. Return 0 if
     418                 :            :  * truncation error (increase precision) */
     419                 :            : static int
     420                 :      15443 : CheckSol(GEN *pS, GEN z1, GEN z2, GEN P, GEN rhs, GEN ro)
     421                 :            : {
     422                 :      15443 :   GEN x, y, ro1 = gel(ro,1), ro2 = gel(ro,2);
     423                 :            :   long e;
     424                 :            : 
     425                 :      15443 :   y = grndtoi(real_i(gdiv(gsub(z2,z1), gsub(ro1,ro2))), &e);
     426         [ +  + ]:      15443 :   if (e > 0) return 0;
     427         [ +  + ]:      15415 :   if (!signe(y)) return 1; /* y = 0 taken care of in SmallSols */
     428                 :      15332 :   x = gadd(z1, gmul(ro1, y));
     429                 :      15332 :   x = grndtoi(real_i(x), &e);
     430         [ -  + ]:      15332 :   if (e > 0) return 0;
     431         [ +  + ]:      15332 :   if (e <= -13)
     432                 :            :   { /* y != 0 and rhs != 0; check whether P(x,y) = rhs or P(-x,-y) = rhs */
     433                 :        100 :     GEN z = poleval(ZX_rescale(P,y),x);
     434         [ +  + ]:        100 :     if (absi_equal(z, rhs)) add_pm(pS, x,y, z, degpol(P), rhs);
     435                 :            :   }
     436                 :      15443 :   return 1;
     437                 :            : }
     438                 :            : 
     439                 :            : static const long EXPO1 = 6;
     440                 :            : static int
     441                 :      35611 : round_to_b(GEN v, long b, GEN Delta2, long i1, GEN L)
     442                 :            : {
     443                 :      35611 :   long i, l = lg(v);
     444         [ +  + ]:      71767 :   for (i = 1; i < l; i++)
     445                 :            :   {
     446                 :            :     GEN c;
     447         [ +  + ]:      56324 :     if (i == i1)
     448                 :      35363 :       c = stoi(b);
     449                 :            :     else
     450                 :            :     {
     451                 :            :       long e;
     452                 :      20961 :       c = gsub(gmulgs(gel(Delta2,i), b), gel(L,i));
     453                 :      20961 :       c = grndtoi(c,&e);
     454         [ +  + ]:      20961 :       if (e > -EXPO1) return 0;
     455                 :            :     }
     456                 :      36156 :     gel(v,i) = c;
     457                 :            :   }
     458                 :      35611 :   return 1;
     459                 :            : }
     460                 :            : /* For each possible value of b_i1, compute the b_i's
     461                 :            : * and 2 conjugates of z = x - alpha y. Then check. */
     462                 :            : static int
     463                 :       2126 : TrySol(GEN *pS, GEN B0, long i1, GEN Delta2, GEN Lambda, GEN ro,
     464                 :            :        GEN NE, GEN MatFU, GEN P, GEN rhs)
     465                 :            : {
     466                 :       2126 :   long bi1, i, B = gtolong(gceil(B0)), l = lg(Delta2);
     467                 :       2126 :   GEN b = cgetg(l,t_VEC), L = cgetg(l,t_VEC);
     468                 :            : 
     469         [ +  + ]:       6278 :   for (i = 1; i < l; i++)
     470                 :            :   {
     471         [ +  + ]:       4152 :     if (i == i1)
     472                 :       2126 :       gel(L,i) = gen_0;
     473                 :            :     else
     474                 :            :     {
     475                 :       2026 :       GEN delta = gel(Delta2,i);
     476                 :       2026 :       gel(L, i) = gsub(gmul(delta,gel(Lambda,i1)), gel(Lambda,i));
     477                 :            :     }
     478                 :            :   }
     479         [ +  + ]:      37709 :   for (bi1 = -B; bi1 <= B; bi1++)
     480                 :            :   {
     481                 :            :     GEN z1, z2;
     482         [ +  + ]:      35611 :     if (!round_to_b(b, bi1, Delta2, i1, L)) continue;
     483                 :      15443 :     z1 = gel(NE,1);
     484                 :      15443 :     z2 = gel(NE,2);
     485         [ +  + ]:      31679 :     for (i = 1; i < l; i++)
     486                 :            :     {
     487                 :      16236 :       GEN c = gel(b,i);
     488                 :      16236 :       z1 = gmul(z1, powgi(gcoeff(MatFU,1,i), c));
     489                 :      16236 :       z2 = gmul(z2, powgi(gcoeff(MatFU,2,i), c));
     490                 :            :     }
     491         [ +  + ]:      15443 :     if (!CheckSol(pS, z1,z2,P,rhs,ro)) return 0;
     492                 :            :   }
     493                 :       2126 :   return 1;
     494                 :            : }
     495                 :            : 
     496                 :            : /* find q1,q2,q3 st q1 + b q2 + c q3 ~ 0 */
     497                 :            : static GEN
     498                 :         90 : GuessQi(GEN b, GEN c, GEN *eps)
     499                 :            : {
     500                 :         90 :   const long shift = 33;
     501                 :         90 :   GEN Q, Lat, C = int2n(shift);
     502                 :            : 
     503                 :         90 :   Lat = matid(3);
     504                 :         90 :   gcoeff(Lat,3,1) = ground(gmul2n(b, shift));
     505                 :         90 :   gcoeff(Lat,3,2) = ground(gmul2n(c, shift));
     506                 :         90 :   gcoeff(Lat,3,3) = C;
     507                 :            : 
     508                 :         90 :   Q = gel(lllint(Lat),1);
     509         [ -  + ]:         90 :   if (gequal0(gel(Q,2))) return NULL; /* FAIL */
     510                 :            : 
     511                 :         90 :   *eps = gadd(gadd(gel(Q,3), gmul(gel(Q,1),b)), gmul(gel(Q,2),c));
     512                 :         90 :   *eps = mpabs(*eps); return Q;
     513                 :            : }
     514                 :            : 
     515                 :            : /* x a t_REAL */
     516                 :            : static GEN
     517         [ +  + ]:        201 : myfloor(GEN x) { return expo(x) > 30 ? ceil_safe(x): floorr(x); }
     518                 :            : 
     519                 :            : /* Check for not-so-small solutions. Return a t_REAL or NULL */
     520                 :            : static GEN
     521                 :        105 : MiddleSols(GEN *pS, GEN bound, GEN roo, GEN P, GEN rhs, long s, GEN c1)
     522                 :            : {
     523                 :            :   long j, k, nmax, d;
     524                 :            :   GEN bndcf;
     525                 :            : 
     526         [ -  + ]:        105 :   if (expo(bound) < 0) return bound;
     527                 :        105 :   d = degpol(P);
     528                 :        105 :   bndcf = sqrtnr(shiftr(c1,1), d - 2);
     529         [ +  + ]:        105 :   if (cmprr(bound, bndcf) < 0) return bound;
     530                 :            :   /* divide by log2((1+sqrt(5))/2)
     531                 :            :    * 1 + ==> ceil
     532                 :            :    * 2 + ==> continued fraction is normalized if last entry is 1
     533                 :            :    * 3 + ==> start at a0, not a1 */
     534                 :         81 :   nmax = 3 + (long)(dbllog2(bound) * 1.44042009041256);
     535                 :         81 :   bound = myfloor(bound);
     536                 :            : 
     537         [ +  + ]:        243 :   for (k = 1; k <= s; k++)
     538                 :            :   {
     539                 :        162 :     GEN ro = real_i(gel(roo,k)), t = gboundcf(ro, nmax);
     540                 :            :     GEN pm1, qm1, p0, q0;
     541                 :            : 
     542                 :        162 :     pm1 = gen_0; p0 = gen_1;
     543                 :        162 :     qm1 = gen_1; q0 = gen_0;
     544                 :            : 
     545         [ +  + ]:       7052 :     for (j = 1; j < lg(t); j++)
     546                 :            :     {
     547                 :            :       GEN p, q, z, Q, R;
     548                 :       7048 :       p = addii(mulii(p0, gel(t,j)), pm1); pm1 = p0; p0 = p;
     549                 :       7048 :       q = addii(mulii(q0, gel(t,j)), qm1); qm1 = q0; q0 = q;
     550         [ +  + ]:       7048 :       if (cmpii(q, bound) > 0) break;
     551         [ -  + ]:       6890 :       if (DEBUGLEVEL >= 2) err_printf("Checking (+/- %Ps, +/- %Ps)\n",p, q);
     552                 :            : 
     553                 :       6890 :       z = poleval(ZX_rescale(P,q), p); /* = P(p/q) q^dep(P) */
     554                 :       6890 :       Q = dvmdii(rhs, z, &R);
     555         [ +  + ]:       6890 :       if (R != gen_0) continue;
     556                 :        130 :       setabssign(Q);
     557         [ +  + ]:        130 :       if (Z_ispowerall(Q, d, &Q))
     558                 :            :       {
     559         [ +  + ]:         55 :         if (!is_pm1(Q)) { p = mulii(p, Q); q = mulii(q, Q); }
     560                 :        288 :         add_pm(pS, p, q, z, d, rhs);
     561                 :            :       }
     562                 :            :     }
     563         [ +  + ]:        162 :     if (j == lg(t))
     564                 :            :     {
     565                 :            :       long prec;
     566         [ -  + ]:          4 :       if (j > nmax) pari_err_BUG("thue [short continued fraction]");
     567                 :            :       /* the theoretical value is bit_prec = gexpo(ro)+1+log2(bound) */
     568                 :          4 :       prec = precdbl(precision(ro));
     569         [ -  + ]:          4 :       if (DEBUGLEVEL>1) pari_warn(warnprec,"thue",prec);
     570                 :          4 :       roo = realroots(P, NULL, prec);
     571         [ -  + ]:          4 :       if (lg(roo)-1 != s) pari_err_BUG("thue [realroots]");
     572                 :          4 :       k--;
     573                 :            :     }
     574                 :            :   }
     575                 :        105 :   return bndcf;
     576                 :            : }
     577                 :            : 
     578                 :            : static void
     579                 :     411150 : check_y_root(GEN *pS, GEN P, GEN Y)
     580                 :            : {
     581                 :     411150 :   GEN r = nfrootsQ(P);
     582                 :            :   long j;
     583         [ +  + ]:     411445 :   for (j = 1; j < lg(r); j++)
     584         [ +  - ]:        295 :     if (typ(gel(r,j)) == t_INT) add_sol(pS, gel(r,j), Y);
     585                 :     411150 : }
     586                 :            : 
     587                 :            : static void
     588                 :     205615 : check_y(GEN *pS, GEN P, GEN poly, GEN Y, GEN rhs)
     589                 :            : {
     590                 :     205615 :   long j, l = lg(poly);
     591                 :     205615 :   GEN Yn = Y;
     592                 :     205615 :   gel(P, l-1) = gel(poly, l-1);
     593         [ +  + ]:     822390 :   for (j = l-2; j >= 2; j--)
     594                 :            :   {
     595                 :     616775 :     gel(P,j) = mulii(Yn, gel(poly,j));
     596         [ +  + ]:     616775 :     if (j > 2) Yn = mulii(Yn, Y);
     597                 :            :   }
     598                 :     205615 :   gel(P,2) = subii(gel(P,2), rhs); /* P = poly(Y/y)*y^deg(poly) - rhs */
     599                 :     205615 :   check_y_root(pS, P, Y);
     600                 :     205615 : }
     601                 :            : 
     602                 :            : /* Check for solutions under a small bound (see paper) */
     603                 :            : static GEN
     604                 :        120 : SmallSols(GEN S, GEN x3, GEN poly, GEN rhs)
     605                 :            : {
     606                 :        120 :   pari_sp av = avma;
     607                 :            :   GEN X, P, rhs2;
     608                 :        120 :   long j, l = lg(poly), n = degpol(poly);
     609                 :            :   ulong y, By;
     610                 :            : 
     611                 :        120 :   x3 = myfloor(x3);
     612                 :            : 
     613         [ -  + ]:        120 :   if (DEBUGLEVEL>1) err_printf("* Checking for small solutions <= %Ps\n", x3);
     614         [ -  + ]:        120 :   if (lgefint(x3) > 3)
     615                 :          0 :     pari_err_OVERFLOW(stack_sprintf("thue (SmallSols): y <= %Ps", x3));
     616                 :        120 :   By = itou(x3);
     617                 :            :   /* y = 0 first: solve X^n = rhs */
     618         [ +  + ]:        120 :   if (odd(n))
     619                 :            :   {
     620         [ +  + ]:         90 :     if (Z_ispowerall(absi(rhs), n, &X))
     621         [ +  + ]:         20 :       add_sol(&S, signe(rhs) > 0? X: negi(X), gen_0);
     622                 :            :   }
     623 [ +  + ][ +  + ]:         30 :   else if (signe(rhs) > 0 && Z_ispowerall(rhs, n, &X))
     624                 :            :   {
     625                 :         10 :     add_sol(&S, X, gen_0);
     626                 :         10 :     add_sol(&S, negi(X), gen_0);
     627                 :            :   }
     628                 :        120 :   rhs2 = shifti(rhs,1);
     629                 :            :   /* y != 0 */
     630                 :        120 :   P = cgetg(l, t_POL); P[1] = poly[1];
     631         [ +  + ]:     205655 :   for (y = 1; y <= By; y++)
     632                 :            :   {
     633                 :     205535 :     pari_sp av2 = avma;
     634                 :     205535 :     long lS = lg(S);
     635                 :     205535 :     GEN Y = utoipos(y);
     636                 :            :     /* try y */
     637                 :     205535 :     check_y(&S, P, poly, Y, rhs);
     638                 :            :     /* try -y */
     639         [ +  + ]:     616095 :     for (j = l-2; j >= 2; j -= 2) togglesign( gel(P,j) );
     640         [ +  + ]:     205535 :     if (j == 0) gel(P,2) = subii(gel(P,2), rhs2);
     641                 :     205535 :     check_y_root(&S, P, utoineg(y));
     642         [ +  + ]:     205535 :     if (lS == lg(S)) { avma = av2; continue; } /* no solution found */
     643                 :            : 
     644         [ -  + ]:        100 :     if (gc_needed(av,1))
     645                 :            :     {
     646         [ #  # ]:          0 :       if(DEBUGMEM>1) pari_warn(warnmem,"SmallSols");
     647                 :          0 :       gerepileall(av, 2, &S, &rhs2);
     648                 :          0 :       P = cgetg(l, t_POL); P[1] = poly[1];
     649                 :            :     }
     650                 :            :   }
     651                 :        120 :   return S;
     652                 :            : }
     653                 :            : 
     654                 :            : /* Computes [x]! */
     655                 :            : static double
     656                 :        115 : fact(double x)
     657                 :            : {
     658                 :        115 :   double ft = 1.0;
     659         [ +  + ]:        535 :   x = floor(x); while (x>1) { ft *= x; x--; }
     660                 :        115 :   return ft ;
     661                 :            : }
     662                 :            : 
     663                 :            : static GEN
     664                 :         70 : RgX_homogenize(GEN P, long v)
     665                 :            : {
     666                 :         70 :   GEN Q = leafcopy(P);
     667                 :         70 :   long i, l = lg(P), d = degpol(P);
     668         [ +  + ]:        300 :   for (i = 2; i < l; i++) gel(Q,i) = monomial(gel(Q,i), d--, v);
     669                 :         70 :   return Q;
     670                 :            : }
     671                 :            : 
     672                 :            : /* Compute all relevant constants needed to solve the equation P(x,y)=a given
     673                 :            :  * the solutions of N_{K/Q}(x)=a (see inithue). */
     674                 :            : GEN
     675                 :        240 : thueinit(GEN pol, long flag, long prec)
     676                 :            : {
     677                 :        240 :   GEN POL, C, L, fa, tnf, bnf = NULL;
     678                 :        240 :   pari_sp av = avma;
     679                 :            :   long s, lfa, dpol;
     680                 :            : 
     681         [ -  + ]:        240 :   if (checktnf(pol)) { bnf = checkbnf(gel(pol,2)); pol = gel(pol,1); }
     682         [ -  + ]:        240 :   if (typ(pol)!=t_POL) pari_err_TYPE("thueinit",pol);
     683                 :        240 :   dpol = degpol(pol);
     684         [ +  + ]:        240 :   if (dpol <= 0) pari_err_CONSTPOL("thueinit");
     685                 :        235 :   RgX_check_ZX(pol, "thueinit");
     686         [ -  + ]:        235 :   if (varn(pol)) { pol = leafcopy(pol); setvarn(pol, 0); }
     687                 :            : 
     688                 :        235 :   POL = Q_primitive_part(pol, &C);
     689                 :        235 :   L = gen_1;
     690                 :        235 :   fa = ZX_factor(POL); lfa = lgcols(fa);
     691 [ +  + ][ +  + ]:        235 :   if (lfa > 2 || itos(gcoeff(fa,1,2)) > 1)
     692                 :            :   { /* reducible polynomial, no need to reduce to the monic case */
     693                 :         65 :     GEN P, Q, R, g, f = gcoeff(fa,1,1), E = gcoeff(fa,1,2);
     694                 :         65 :     long e = itos(E);
     695                 :         65 :     long vy = fetch_var();
     696                 :         65 :     long va = fetch_var();
     697                 :         65 :     long vb = fetch_var();
     698         [ +  + ]:         65 :     C = C? ginv(C): gen_1;
     699         [ +  + ]:         65 :     if (e != 1)
     700                 :            :     {
     701         [ +  + ]:         45 :       if (lfa == 2) {
     702                 :         30 :         tnf = mkvec3(mkvec3(POL,C,L), thueinit(f, flag, prec), E);
     703                 :         30 :         delete_var(); delete_var(); delete_var();
     704                 :         30 :         return gerepilecopy(av, tnf);
     705                 :            :       }
     706                 :         15 :       P = gpowgs(f,e);
     707                 :            :     }
     708                 :            :     else
     709                 :         20 :       P = f;
     710                 :         35 :     g = RgX_div(POL, P);
     711                 :         35 :     P = RgX_Rg_sub(RgX_homogenize(f, vy), pol_x(va));
     712                 :         35 :     Q = RgX_Rg_sub(RgX_homogenize(g, vy), pol_x(vb));
     713                 :         35 :     R = polresultant0(P, Q, -1, 0);
     714                 :         35 :     tnf = mkvec3(mkvec3(POL,C,L), mkvecsmall4(degpol(f), e, va,vb),  R);
     715                 :         35 :     delete_var(); delete_var(); delete_var();
     716                 :         35 :     return gerepilecopy(av, tnf);
     717                 :            :   }
     718                 :            :   /* POL monic: POL(x) = C pol(x/L), L integer */
     719                 :        170 :   POL = ZX_primitive_to_monic(POL, &L);
     720                 :        170 :   C = gdiv(powiu(L, dpol), gel(pol, dpol+2));
     721                 :        170 :   pol = POL;
     722                 :            : 
     723                 :        170 :   s = ZX_sturm(pol);
     724         [ +  + ]:        170 :   if (s)
     725                 :            :   {
     726                 :        115 :     long PREC, n = degpol(pol);
     727                 :        115 :     double d, dr, dn = (double)n;
     728                 :            : 
     729         [ -  + ]:        115 :     if (dpol <= 2) pari_err_DOMAIN("thueinit", "P","=",pol,pol);
     730                 :        115 :     dr = (double)((s+n-2)>>1); /* s+t-1 */
     731                 :        115 :     d = dn*(dn-1)*(dn-2);
     732                 :            :     /* Guess precision by approximating Baker's bound. The guess is most of
     733                 :            :      * the time not sharp, ie 10 to 30 decimal digits above what is _really_
     734                 :            :      * necessary. Note that the limiting step is the reduction. See paper. */
     735                 :        115 :     PREC = nbits2prec((long)((5.83 + (dr+4)*5 + log(fact(dr+3)) + (dr+3)*log(dr+2) +
     736                 :        230 :                      (dr+3)*log(d) + log(log(2*d*(dr+2))) + (dr+1))
     737                 :        115 :                      /10.)*32+32);
     738                 :            : 
     739         [ +  + ]:        115 :     if (flag == 0) PREC = (long)(2.2 * PREC); /* Lazy, to be improved */
     740         [ -  + ]:        115 :     if (PREC < prec) PREC = prec;
     741         [ -  + ]:        115 :     if (DEBUGLEVEL >=2) err_printf("prec = %d\n", PREC);
     742                 :            : 
     743                 :            :     for (;;)
     744                 :            :     {
     745         [ +  - ]:        115 :       if (( tnf = inithue(pol, bnf, flag, PREC) )) break;
     746                 :          0 :       PREC = precdbl(PREC);
     747         [ #  # ]:          0 :       if (DEBUGLEVEL>1) pari_warn(warnprec,"thueinit",PREC);
     748                 :          0 :       bnf = NULL; avma = av;
     749                 :        115 :     }
     750                 :            :   }
     751                 :            :   else
     752                 :            :   {
     753                 :            :     GEN ro, c0;
     754                 :            :     long k,l;
     755         [ +  - ]:         55 :     if (!bnf)
     756                 :            :     {
     757                 :         55 :       bnf = Buchall(pol, nf_FORCE, DEFAULTPREC);
     758         [ -  + ]:         55 :       if (flag) (void)bnfcertify(bnf);
     759                 :            :     }
     760                 :         55 :     ro = nf_get_roots(bnf_get_nf(bnf));
     761                 :         55 :     l = lg(ro);
     762                 :         55 :     c0 = imag_i(gel(ro,1));
     763         [ +  + ]:         75 :     for (k = 2; k < l; k++) c0 = mulrr(c0, imag_i(gel(ro,k)));
     764                 :         55 :     c0 = invr(sqrr(c0)); tnf = mkvec3(pol,bnf,c0);
     765                 :            :   }
     766                 :        170 :   gel(tnf,1) = mkvec3(gel(tnf,1), C, L);
     767                 :        235 :   return gerepilecopy(av,tnf);
     768                 :            : }
     769                 :            : 
     770                 :            : /* arg(t^2) / 2Pi; arg(t^2) = arg(t/conj(t)) */
     771                 :            : static GEN
     772                 :        300 : argsqr(GEN t, GEN Pi)
     773                 :            : {
     774                 :        300 :   GEN v, u = divrr(garg(t,0), Pi); /* in -1 < u <= 1 */
     775                 :            :   /* reduce mod Z to -1/2 < u <= 1/2 */
     776         [ +  + ]:        300 :   if (signe(u) > 0)
     777                 :            :   {
     778                 :        145 :     v = subrs(u,1); /* ]-1,0] */
     779         [ +  + ]:        145 :     if (absr_cmp(v,u) < 0) u = v;
     780                 :            :   }
     781                 :            :   else
     782                 :            :   {
     783                 :        155 :     v = addrs(u,1);/* ]0,1] */
     784         [ +  + ]:        155 :     if (absr_cmp(v,u) <= 0) u = v;
     785                 :            :   }
     786                 :        300 :   return u;
     787                 :            : }
     788                 :            : /* i1 != i2 */
     789                 :            : static void
     790                 :       3885 : init_get_B(long i1, long i2, GEN Delta2, GEN Lambda, GEN Deps5, baker_s *BS,
     791                 :            :            long prec)
     792                 :            : {
     793                 :            :   GEN delta, lambda;
     794         [ +  + ]:       3885 :   if (BS->r > 1)
     795                 :            :   {
     796                 :       3735 :     delta = gel(Delta2,i2);
     797                 :       3735 :     lambda = gsub(gmul(delta,gel(Lambda,i1)), gel(Lambda,i2));
     798         [ +  + ]:       3735 :     if (Deps5) BS->inverrdelta = divrr(Deps5, addsr(1,delta));
     799                 :            :   }
     800                 :            :   else
     801                 :            :   { /* r == 1: i1 = s = t = 1; i2 = 2 */
     802                 :        150 :     GEN fu = gel(BS->MatFU,1), ro = BS->ro, t;
     803                 :            : 
     804                 :        150 :     t = gel(fu,2);
     805                 :        150 :     delta = argsqr(t, BS->Pi);
     806         [ +  + ]:        150 :     if (Deps5) BS->inverrdelta = shiftr(gabs(t,prec), prec2nbits(prec)-1);
     807                 :            : 
     808                 :        150 :     t = gmul(gsub(gel(ro,1), gel(ro,2)), gel(BS->NE,3));
     809                 :        150 :     lambda = argsqr(t, BS->Pi);
     810                 :            :   }
     811                 :       3885 :   BS->delta = delta;
     812                 :       3885 :   BS->lambda = lambda;
     813                 :       3885 : }
     814                 :            : 
     815                 :            : static GEN
     816                 :       2131 : get_B0(long i1, GEN Delta2, GEN Lambda, GEN Deps5, long prec, baker_s *BS)
     817                 :            : {
     818                 :       2131 :   GEN B0 = Baker(BS);
     819         [ +  + ]:       2131 :   long step = 0, i2 = (i1 == 1)? 2: 1;
     820                 :            :   for(;;) /* i2 from 1 to r unless r = 1 [then i2 = 2] */
     821                 :            :   {
     822                 :       2131 :     init_get_B(i1,i2, Delta2,Lambda,Deps5, BS, prec);
     823         [ -  + ]:       2131 :     if (DEBUGLEVEL>1) err_printf("  inverrdelta = %Ps\n",BS->inverrdelta);
     824         [ -  + ]:       2131 :     if (DEBUGLEVEL>1) err_printf("  Entering CF...\n");
     825                 :            :     /* Reduce B0 as long as we make progress: newB0 < oldB0 - 0.1 */
     826                 :            :     for (;;)
     827                 :            :     {
     828                 :       6777 :       GEN oldB0 = B0, kappa = utoipos(10);
     829                 :            :       long cf;
     830                 :            : 
     831         [ +  - ]:       7939 :       for (cf = 0; cf < 10; cf++, kappa = muliu(kappa,10))
     832                 :            :       {
     833                 :       7939 :         int res = CF_1stPass(&B0, kappa, BS);
     834         [ +  + ]:       7939 :         if (res < 0) return NULL; /* prec problem */
     835         [ +  + ]:       7934 :         if (res) break;
     836         [ -  + ]:       1162 :         if (DEBUGLEVEL>1) err_printf("CF failed. Increasing kappa\n");
     837                 :            :       }
     838 [ +  + ][ -  + ]:       6772 :       if (!step && cf == 10)
     839                 :            :       { /* Semirational or totally rational case */
     840                 :            :         GEN Q, ep, q, l0, denbound;
     841                 :            : 
     842         [ #  # ]:          0 :         if (! (Q = GuessQi(BS->delta, BS->lambda, &ep)) ) break;
     843                 :            : 
     844                 :          0 :         denbound = gadd(B0, absi(gel(Q,1)));
     845                 :          0 :         q = denom( bestappr(BS->delta, denbound) );
     846                 :          0 :         l0 = subrr(errnum(BS->delta, q), ep);
     847         [ #  # ]:          0 :         if (signe(l0) <= 0) break;
     848                 :            : 
     849                 :          0 :         B0 = divrr(mplog(divrr(mulir(gel(Q,2), BS->c15), l0)),  BS->c13);
     850         [ #  # ]:          0 :         if (DEBUGLEVEL>1) err_printf("Semirat. reduction: B0 -> %Ps\n",B0);
     851                 :            :       }
     852                 :            :       /* if no progress, stop */
     853         [ +  + ]:       6772 :       if (gcmp(oldB0, gadd(B0,dbltor(0.1))) <= 0) return gmin(oldB0, B0);
     854                 :       4646 :       else step++;
     855                 :       4646 :     }
     856         [ #  # ]:          0 :     i2++; if (i2 == i1) i2++;
     857         [ #  # ]:          0 :     if (i2 > BS->r) break;
     858                 :          0 :   }
     859                 :          0 :   pari_err_BUG("thue (totally rational case)");
     860                 :       2131 :   return NULL; /* not reached */
     861                 :            : }
     862                 :            : 
     863                 :            : static GEN
     864                 :       1754 : get_Bx_LLL(long i1, GEN Delta2, GEN Lambda, long prec, baker_s *BS)
     865                 :            : {
     866                 :       1754 :   GEN B0 = Baker(BS), Bx = NULL;
     867         [ +  + ]:       1754 :   long step = 0, i2 = (i1 == 1)? 2: 1;
     868                 :            :   for(;;) /* i2 from 1 to r unless r = 1 [then i2 = 2] */
     869                 :            :   {
     870                 :       1754 :     init_get_B(i1,i2, Delta2,Lambda,NULL, BS, prec);
     871         [ -  + ]:       1754 :     if (DEBUGLEVEL>1) err_printf("  Entering LLL...\n");
     872                 :            :     /* Reduce B0 as long as we make progress: newB0 < oldB0 - 0.1 */
     873                 :            :     for (;;)
     874                 :            :     {
     875                 :       7258 :       GEN oldBx = Bx, kappa = utoipos(10);
     876                 :       7258 :       const long cfMAX = 10;
     877                 :            :       long cf;
     878                 :            : 
     879         [ +  + ]:      12884 :       for (cf = 0; cf < cfMAX; cf++, kappa = muliu(kappa,10))
     880                 :            :       {
     881                 :      12689 :         int res = LLL_1stPass(&B0, kappa, BS, &Bx);
     882         [ +  + ]:      12689 :         if (res < 0) return NULL;
     883         [ +  + ]:      12680 :         if (res) break;
     884         [ -  + ]:       5626 :         if (DEBUGLEVEL>1) err_printf("LLL failed. Increasing kappa\n");
     885                 :            :       }
     886                 :            : 
     887                 :            :       /* FIXME: TO BE COMPLETED */
     888 [ +  + ][ +  + ]:       7249 :       if (!step && cf == cfMAX)
     889                 :            :       { /* Semirational or totally rational case */
     890                 :            :         GEN Q, ep, q, l0, denbound;
     891                 :            : 
     892         [ +  - ]:         90 :         if (! (Q = GuessQi(BS->delta, BS->lambda, &ep)) ) break;
     893                 :            : 
     894                 :            :         /* Q[2] != 0 */
     895                 :         90 :         denbound = gadd(mulri(B0, absi(gel(Q,1))),
     896                 :         90 :                         mulii(BS->Ind, absi(gel(Q,2))));
     897                 :         90 :         q = denom( bestappr(BS->delta, denbound) );
     898                 :         90 :         l0 = divri(subrr(errnum(BS->delta, q), ep), absi(gel(Q,2)));
     899         [ +  - ]:         90 :         if (signe(l0) <= 0) break;
     900                 :            : 
     901                 :         90 :         get_B0Bx(BS, l0, &B0, &Bx);
     902         [ -  + ]:         90 :         if (DEBUGLEVEL>1)
     903                 :         90 :           err_printf("Semirat. reduction: B0 -> %Ps x <= %Ps\n",B0, Bx);
     904                 :            :       }
     905                 :            :       /* if no progress, stop */
     906 [ +  + ][ +  + ]:       7249 :       if (oldBx && gcmp(oldBx, Bx) <= 0) return oldBx; else step++;
     907                 :       5504 :     }
     908         [ #  # ]:          0 :     i2++; if (i2 == i1) i2++;
     909         [ #  # ]:          0 :     if (i2 > BS->r) break;
     910                 :          0 :   }
     911                 :          0 :   pari_err_BUG("thue (totally rational case)");
     912                 :       1754 :   return NULL; /* not reached */
     913                 :            : }
     914                 :            : 
     915                 :            : static GEN
     916                 :        105 : LargeSols(GEN P, GEN tnf, GEN rhs, GEN ne)
     917                 :            : {
     918                 :        105 :   GEN S = NULL, Delta0, ro, bnf, MatFU, A, csts, dP, Bx;
     919                 :            :   GEN c1,c2,c3,c4,c14, x0, x1, x2, x3, tmp, eps5;
     920                 :            :   long iroot, ine, n, r, Prec, prec, s,t;
     921                 :            :   baker_s BS;
     922                 :        105 :   pari_sp av = avma;
     923                 :            : 
     924                 :            : START:
     925         [ +  + ]:        147 :   if (S) /* restart from precision problems */
     926                 :            :   {
     927                 :         42 :     avma = av; prec += nbits2extraprec(320);
     928         [ -  + ]:         42 :     if (DEBUGLEVEL>1) pari_warn(warnprec,"thue",prec);
     929                 :         42 :     tnf = inithue(P, bnf, 0, prec);
     930                 :            :   }
     931                 :        147 :   S = cgetg(1, t_VEC);
     932                 :        147 :   bnf  = gel(tnf,2);
     933                 :        147 :   csts = gel(tnf,7);
     934                 :        147 :   nf_get_sign(bnf_get_nf(bnf), &s, &t);
     935                 :        147 :   BS.r = r = s+t-1; n = degpol(P);
     936                 :        147 :   ro     = gel(tnf,3);
     937                 :        147 :   BS.ALH = gel(tnf,4);
     938                 :        147 :   MatFU  = gel(tnf,5);
     939                 :        147 :   A      = gel(tnf,6);
     940                 :        147 :   c1     = gel(csts,1); c1 = gmul(absi(rhs), c1);
     941                 :        147 :   c2     = gel(csts,2);
     942                 :        147 :   BS.hal = gel(csts,3);
     943                 :        147 :   x0     = gel(csts,4);
     944                 :        147 :   eps5   = gel(csts,5);
     945                 :        147 :   Prec = gtolong(gel(csts,6));
     946                 :        147 :   BS.Ind = gel(csts,7);
     947                 :        147 :   BS.MatFU = MatFU;
     948                 :        147 :   BS.bak = muluu(n, (n-1)*(n-2)); /* safe */
     949                 :        147 :   BS.deg = n;
     950                 :            : 
     951         [ +  + ]:        147 :   if (t) x0 = gmul(x0, absisqrtn(rhs, n, Prec));
     952                 :        147 :   tmp = divrr(c1,c2);
     953                 :        147 :   c3 = mulrr(dbltor(1.39), tmp);
     954                 :        147 :   c4 = mulur(n-1, c3);
     955                 :        147 :   c14 = mulrr(c4, vecmax_shallow(RgM_sumcol(gabs(A,DEFAULTPREC))));
     956                 :            : 
     957                 :        147 :   x1 = gmax(x0, sqrtnr(shiftr(tmp,1),n));
     958                 :        147 :   x2 = gmax(x1, sqrtnr(mulur(10,c14), n));
     959                 :        147 :   x3 = gmax(x2, sqrtnr(shiftr(c14, EXPO1+1),n));
     960         [ -  + ]:        147 :   if (DEBUGLEVEL>1) {
     961                 :          0 :     err_printf("x1 -> %Ps\n",x1);
     962                 :          0 :     err_printf("x2 -> %Ps\n",x2);
     963                 :          0 :     err_printf("c14 = %Ps\n",c14);
     964                 :            :   }
     965                 :            : 
     966                 :        147 :   dP = ZX_deriv(P);
     967                 :        147 :   Delta0 = RgM_sumcol(A);
     968                 :            : 
     969         [ +  + ]:        427 :   for (iroot=1; iroot<=s; iroot++)
     970                 :            :   {
     971                 :        280 :     GEN Delta = Delta0, Delta2, D, Deps5, MatNE, Hmu, diffRo, c5, c7;
     972                 :            :     long i1, k;
     973                 :            : 
     974         [ +  + ]:        280 :     if (iroot <= r) Delta = RgC_add(Delta, RgC_Rg_mul(gel(A,iroot), stoi(-n)));
     975                 :        280 :     D = gabs(Delta,Prec); i1 = vecindexmax(D);
     976                 :        280 :     c5 = gel(D, i1);
     977                 :        280 :     Delta2 = RgC_Rg_div(Delta, gel(Delta, i1));
     978                 :        280 :     c5  = myround(gprec_w(c5,DEFAULTPREC), 1);
     979                 :        280 :     Deps5 = divrr(subrr(c5,eps5), eps5);
     980                 :        280 :     c7  = mulur(r,c5);
     981                 :        280 :     BS.c10 = divur(n,c7);
     982                 :        280 :     BS.c13 = divur(n,c5);
     983         [ -  + ]:        280 :     if (DEBUGLEVEL>1) {
     984                 :          0 :       err_printf("* real root no %ld/%ld\n", iroot,s);
     985                 :          0 :       err_printf("  c10 = %Ps\n",BS.c10);
     986                 :          0 :       err_printf("  c13 = %Ps\n",BS.c13);
     987                 :            :     }
     988                 :            : 
     989                 :        280 :     prec = Prec;
     990                 :            :     for (;;)
     991                 :            :     {
     992         [ +  - ]:        280 :       if (( MatNE = Conj_LH(ne, &Hmu, ro, prec) )) break;
     993                 :          0 :       prec = precdbl(prec);
     994         [ #  # ]:          0 :       if (DEBUGLEVEL>1) pari_warn(warnprec,"thue",prec);
     995                 :          0 :       ro = tnf_get_roots(P, prec, s, t);
     996                 :          0 :     }
     997                 :        280 :     BS.ro    = ro;
     998                 :        280 :     BS.iroot = iroot;
     999                 :        280 :     BS.Pi  = mppi(prec);
    1000                 :        280 :     BS.Pi2 = Pi2n(1,prec);
    1001                 :        280 :     diffRo = cgetg(r+1, t_VEC);
    1002         [ +  + ]:        845 :     for (k=1; k<=r; k++)
    1003                 :            :     {
    1004                 :        565 :       GEN z = gel(ro,k);
    1005                 :        235 :       z = (k == iroot)? gdiv(rhs, poleval(dP, z))
    1006         [ +  + ]:        800 :                       : gsub(gel(ro,iroot), z);
    1007                 :        565 :       gel(diffRo,k) = gabs(z, prec);
    1008                 :            :     }
    1009         [ +  + ]:       4123 :     for (ine=1; ine<lg(ne); ine++)
    1010                 :            :     {
    1011                 :       3885 :       pari_sp av2 = avma;
    1012                 :            :       GEN Lambda, B0, c6, c8;
    1013                 :       3885 :       GEN NE = gel(MatNE,ine), v = cgetg(r+1,t_COL);
    1014                 :            : 
    1015         [ -  + ]:       3885 :       if (DEBUGLEVEL>1) err_printf("  - norm sol. no %ld/%ld\n",ine,lg(ne)-1);
    1016         [ +  + ]:      11970 :       for (k=1; k<=r; k++)
    1017                 :            :       {
    1018                 :       8085 :         GEN z = gdiv(gel(diffRo,k), gabs(gel(NE,k), prec));
    1019                 :       8085 :         gel(v,k) = glog(z, prec);
    1020                 :            :       }
    1021                 :       3885 :       Lambda = RgM_RgC_mul(A,v);
    1022                 :            : 
    1023                 :       3885 :       c6 = addrr(dbltor(0.1), vecmax_shallow(gabs(Lambda,DEFAULTPREC)));
    1024                 :       3885 :       c6 = myround(c6, 1);
    1025                 :       3885 :       c8 = addrr(dbltor(1.23), mulur(r,c6));
    1026                 :       3885 :       BS.c11= mulrr(shiftr(c3,1) , mpexp(divrr(mulur(n,c8),c7)));
    1027                 :       3885 :       BS.c15= mulrr(shiftr(c14,1), mpexp(divrr(mulur(n,c6),c5)));
    1028                 :       3885 :       BS.NE = NE;
    1029                 :       3885 :       BS.Hmu = gel(Hmu,ine);
    1030         [ +  + ]:       3885 :       if (is_pm1(BS.Ind))
    1031                 :            :       {
    1032   [ +  +  +  + ]:       4257 :         if (! (B0 = get_B0(i1, Delta2, Lambda, Deps5, prec, &BS)) ||
    1033                 :       2126 :             !TrySol(&S, B0, i1, Delta2, Lambda, ro, NE,MatFU, P,rhs))
    1034                 :            :           goto START;
    1035                 :            :       }
    1036                 :            :       else
    1037                 :            :       {
    1038         [ +  + ]:       1754 :         if (! (Bx = get_Bx_LLL(i1, Delta2, Lambda, prec, &BS)) )
    1039                 :            :            goto START;
    1040                 :       1745 :         x3 = gerepileupto(av2, gmax(Bx, x3));
    1041                 :            :       }
    1042                 :            :     }
    1043                 :            :   }
    1044                 :        105 :   x3 = gmax(x0, MiddleSols(&S, x3, ro, P, rhs, s, c1));
    1045                 :        105 :   return SmallSols(S, x3, P, rhs);
    1046                 :            : }
    1047                 :            : 
    1048                 :            : /* restrict to solutions (x,y) with L | x, replacing each by (x/L, y) */
    1049                 :            : static GEN
    1050                 :        195 : filter_sol_x(GEN S, GEN L)
    1051                 :            : {
    1052                 :            :   long i, k, l;
    1053         [ +  + ]:        195 :   if (is_pm1(L)) return S;
    1054                 :         40 :   l = lg(S); k = 1;
    1055         [ +  + ]:        350 :   for (i = 1; i < l; i++)
    1056                 :            :   {
    1057                 :        310 :     GEN s = gel(S,i), r;
    1058                 :        310 :     gel(s,1) = dvmdii(gel(s,1), L, &r);
    1059         [ +  - ]:        310 :     if (r == gen_0) gel(S, k++) = s;
    1060                 :            :   }
    1061                 :        195 :   setlg(S, k); return S;
    1062                 :            : }
    1063                 :            : 
    1064                 :            : static GEN bnfisintnorm_i(GEN bnf, long s, GEN z);
    1065                 :            : static GEN
    1066                 :        100 : tnf_get_Ind(GEN tnf) { return gmael(tnf,7,7); }
    1067                 :            : static GEN
    1068                 :        160 : tnf_get_bnf(GEN tnf) { return gel(tnf,2); }
    1069                 :            : 
    1070                 :            : static void
    1071                 :          0 : maybe_warn(GEN bnf, GEN a, GEN Ind)
    1072                 :            : {
    1073 [ #  # ][ #  # ]:          0 :   if (!is_pm1(Ind) && !is_pm1(bnf_get_no(bnf)) && !is_pm1(a))
                 [ #  # ]
    1074                 :          0 :     pari_warn(warner, "The result returned by 'thue' is conditional on the GRH");
    1075                 :          0 : }
    1076                 :            : /* return solutions of Norm(x) = a mod U(K)^+ */
    1077                 :            : static GEN
    1078                 :        125 : get_ne(GEN bnf, GEN a, GEN Ind)
    1079                 :            : {
    1080         [ -  + ]:        125 :   if (DEBUGLEVEL) maybe_warn(bnf,a,Ind);
    1081                 :        125 :   return bnfisintnorm(bnf, a);
    1082                 :            : }
    1083                 :            : /* return solutions of |Norm(x)| = |a| mod U(K) */
    1084                 :            : static GEN
    1085                 :         20 : get_neabs(GEN bnf, GEN a, GEN Ind)
    1086                 :            : {
    1087         [ -  + ]:         20 :   if (DEBUGLEVEL) maybe_warn(bnf,a,Ind);
    1088                 :         20 :   return bnfisintnormabs(bnf, a);
    1089                 :            : }
    1090                 :            : 
    1091                 :            : /* Let P(z)=z^2+Bz+C, convert t=u+v*z (mod P) solution of norm equation N(t)=A
    1092                 :            :  * to [x,y] = [u,-v] form: y^2P(x/y) = A */
    1093                 :            : static GEN
    1094                 :        350 : ne2_to_xy(GEN t)
    1095                 :            : {
    1096                 :            :   GEN u,v;
    1097      [ -  +  + ]:        350 :   switch(degpol(t))
    1098                 :            :   {
    1099                 :          0 :     case -1: u = v = gen_0; break;
    1100                 :         20 :     case 0: u = gel(t,2); v = gen_0; break;
    1101                 :        330 :     default: u = gel(t,2); v = gel(t,3);
    1102                 :            :   }
    1103                 :        350 :   return mkvec2(u, gneg(v));
    1104                 :            : }
    1105                 :            : static GEN
    1106                 :         30 : ne2V_to_xyV(GEN v)
    1107                 :            : {
    1108                 :            :   long i, l;
    1109                 :         30 :   GEN w = cgetg_copy(v,&l);
    1110         [ +  + ]:        380 :   for (i=1; i<l; i++) gel(w,i) = ne2_to_xy(gel(v,i));
    1111                 :         30 :   return w;
    1112                 :            : }
    1113                 :            : 
    1114                 :            : static GEN
    1115                 :         15 : sol_0(void) { retmkvec( mkvec2(gen_0,gen_0) ); }
    1116                 :            : static void
    1117                 :        100 : sols_from_R(GEN Rab, GEN *pS, GEN P, GEN POL, GEN rhs)
    1118                 :            : {
    1119                 :        100 :   GEN ry = nfrootsQ(Rab);
    1120                 :        100 :   long k, l = lg(ry);
    1121         [ +  + ]:        200 :   for (k = 1; k < l; k++)
    1122         [ +  + ]:        100 :     if (typ(gel(ry,k)) == t_INT) check_y(pS, P, POL, gel(ry,k), rhs);
    1123                 :        100 : }
    1124                 :            : /* Given a tnf structure as returned by thueinit, a RHS and
    1125                 :            :  * optionally the solutions to the norm equation, returns the solutions to
    1126                 :            :  * the Thue equation F(x,y)=a */
    1127                 :            : GEN
    1128                 :        230 : thue(GEN tnf, GEN rhs, GEN ne)
    1129                 :            : {
    1130                 :        230 :   pari_sp av = avma;
    1131                 :            :   GEN POL, C, L, x3, S;
    1132                 :            : 
    1133         [ +  + ]:        230 :   if (typ(tnf) == t_POL) tnf = thueinit(tnf, 0, DEFAULTPREC);
    1134         [ -  + ]:        230 :   if (!checktnf(tnf)) pari_err_TYPE("thue [please apply thueinit()]", tnf);
    1135         [ -  + ]:        230 :   if (typ(rhs) != t_INT) pari_err_TYPE("thue",rhs);
    1136 [ +  + ][ -  + ]:        230 :   if (ne && typ(ne) != t_VEC) pari_err_TYPE("thue",ne);
    1137                 :            : 
    1138                 :            :   /* solve P(x,y) = rhs <=> POL(L x, y) = C rhs, with POL monic in Z[X] */
    1139                 :        230 :   POL = gel(tnf,1);
    1140                 :        230 :   C = gel(POL,2); rhs = gmul(C, rhs);
    1141         [ -  + ]:        230 :   if (typ(rhs) != t_INT) { avma = av; return cgetg(1, t_VEC); }
    1142         [ +  + ]:        230 :   if (!signe(rhs))
    1143                 :            :   {
    1144                 :         20 :     GEN v = gel(tnf,2);
    1145                 :         20 :     avma = av;
    1146                 :            :     /* at least 2 irreducible factors, one of which has degree 1 */
    1147 [ +  + ][ +  + ]:         20 :     if (typ(v) == t_VECSMALL && v[1] ==1)
    1148                 :          5 :       pari_err_DOMAIN("thue","#sols","=",strtoGENstr("oo"),rhs);
    1149                 :         15 :     return sol_0();
    1150                 :            :   }
    1151                 :        210 :   L = gel(POL,3);
    1152                 :        210 :   POL = gel(POL,1);
    1153         [ +  + ]:        210 :   if (lg(tnf) == 8)
    1154                 :            :   {
    1155         [ +  + ]:        105 :     if (!ne) ne = get_ne(tnf_get_bnf(tnf), rhs, tnf_get_Ind(tnf));
    1156         [ -  + ]:        105 :     if (lg(ne) == 1) { avma = av; return cgetg(1, t_VEC); }
    1157                 :        105 :     S = LargeSols(POL, tnf, rhs, ne);
    1158                 :            :   }
    1159         [ +  + ]:        105 :   else if (typ(gel(tnf,3)) == t_REAL)
    1160                 :            :   { /* Case s=0. All solutions are "small". */
    1161                 :         55 :     GEN bnf = tnf_get_bnf(tnf);
    1162                 :         55 :     GEN c0 = gel(tnf,3), F;
    1163                 :         55 :     x3 = sqrtnr(mulir(absi(rhs),c0), degpol(POL));
    1164                 :         55 :     x3 = addrr(x3, dbltor(0.1)); /* guard from round-off errors */
    1165                 :         55 :     S = cgetg(1,t_VEC);
    1166 [ +  + ][ +  + ]:         55 :     if (!ne && expo(x3) > 10)
    1167                 :            :     {
    1168                 :            :       long l;
    1169                 :            :       GEN P;
    1170         [ +  - ]:         40 :       if (expi(rhs) < 150) F = Z_factor(rhs);
    1171                 :            :       else
    1172                 :            :       {
    1173                 :          0 :         F = Z_factor_limit(rhs, 500000);
    1174                 :          0 :         P = gel(F,1); l = lg(P);
    1175 [ #  # ][ #  # ]:          0 :         if (!is_pm1(rhs) && !BPSW_psp(gel(P,l-1))) F = NULL;
    1176                 :            :       }
    1177         [ +  - ]:         40 :       if (F) ne = get_ne(bnf, F, gen_1);
    1178                 :            :     }
    1179         [ +  + ]:         55 :     if (ne)
    1180                 :            :     {
    1181         [ +  + ]:         45 :       if (lg(ne) == 1) { avma = av; return cgetg(1,t_VEC); }
    1182         [ +  + ]:         35 :       if (degpol(POL) == 2) /* quadratic imaginary */
    1183                 :            :       {
    1184                 :         30 :         GEN u = bnf_get_tuU(bnf);
    1185                 :         30 :         long w =  bnf_get_tuN(bnf);
    1186         [ +  + ]:         30 :         if (w == 4) /* u = I */
    1187                 :         10 :           ne = shallowconcat(ne, RgXQV_RgXQ_mul(ne,u,POL));
    1188         [ +  + ]:         20 :         else if (w == 6) /* u = j */
    1189                 :            :         {
    1190                 :         10 :           GEN u2 = RgXQ_sqr(u,POL);
    1191                 :         10 :           ne = shallowconcat1(mkvec3(ne, RgXQV_RgXQ_mul(ne,u,POL),
    1192                 :            :                                          RgXQV_RgXQ_mul(ne,u2,POL)));
    1193                 :            :         }
    1194                 :         30 :         ne = ne2V_to_xyV(ne);
    1195                 :         30 :         S = shallowconcat(ne, RgV_neg(ne));
    1196                 :            :       }
    1197                 :            :     }
    1198         [ +  + ]:         45 :     if (lg(S) == 1) S = SmallSols(S, x3, POL, rhs);
    1199                 :            :   }
    1200         [ +  + ]:         50 :   else if (typ(gel(tnf,3)) == t_INT) /* reducible case, pure power*/
    1201                 :            :   {
    1202                 :         25 :     GEN bnf, ne1 = NULL, ne2 = NULL;
    1203                 :         25 :     long e = itos( gel(tnf,3) );
    1204         [ +  + ]:         25 :     if (!Z_ispowerall(rhs, e, &rhs)) { avma = av; return cgetg(1, t_VEC); }
    1205                 :         20 :     tnf = gel(tnf,2);
    1206                 :         20 :     bnf = tnf_get_bnf(tnf);
    1207         [ +  + ]:         20 :     ne = get_neabs(bnf, rhs, lg(tnf)==8?tnf_get_Ind(tnf): gen_1);
    1208                 :         20 :     ne1= bnfisintnorm_i(bnf,1,ne);
    1209                 :         20 :     S = thue(tnf, rhs, ne1);
    1210 [ +  + ][ +  + ]:         20 :     if (!odd(e) && lg(tnf)==8) /* if s=0, norms are positive */
    1211                 :            :     {
    1212                 :          5 :       ne2 = bnfisintnorm_i(bnf,-1,ne);
    1213                 :          5 :       S = shallowconcat(S, thue(tnf, negi(rhs), ne2));
    1214                 :            :     }
    1215                 :            :   }
    1216                 :            :   else /* other reducible cases */
    1217                 :            :   { /* solve f^e * g = rhs, f irreducible factor of smallest degree */
    1218                 :         25 :     GEN P, D, v = gel(tnf, 2), R = gel(tnf, 3);
    1219                 :         25 :     long i, l, e = v[2], va = v[3], vb = v[4];
    1220                 :         25 :     P = cgetg(lg(POL), t_POL); P[1] = POL[1];
    1221                 :         25 :     D = divisors(rhs); l = lg(D);
    1222                 :         25 :     S = cgetg(1,t_VEC);
    1223         [ +  + ]:        155 :     for (i = 1; i < l; i++)
    1224                 :            :     {
    1225                 :        130 :       GEN Rab, df = gel(D,i), dg = gel(D,l-i); /* df*dg=|rhs| */
    1226 [ +  + ][ +  + ]:        130 :       if (e > 1 && !Z_ispowerall(df, e, &df)) continue;
    1227                 :            :       /* Rab: univariate polynomial in Z[Y], whose roots are the possible y. */
    1228                 :            :       /* Here and below, Rab != 0 */
    1229         [ +  + ]:         50 :       if (signe(rhs) < 0) dg = negi(dg); /* df*dg=rhs */
    1230                 :         50 :       Rab = gsubst(gsubst(R, va, df), vb, dg);
    1231                 :         50 :       sols_from_R(Rab, &S,P,POL,rhs);
    1232         [ +  + ]:         50 :       Rab = gsubst(gsubst(R, va, negi(df)), vb, odd(e)? negi(dg): dg);
    1233                 :         50 :       sols_from_R(Rab, &S,P,POL,rhs);
    1234                 :            :     }
    1235                 :            :   }
    1236                 :        195 :   S = filter_sol_x(S, L);
    1237                 :        195 :   S = gen_sort_uniq(S, (void*)lexcmp, cmp_nodata);
    1238                 :        225 :   return gerepileupto(av, S);
    1239                 :            : }
    1240                 :            : 
    1241                 :            : /********************************************************************/
    1242                 :            : /**                                                                **/
    1243                 :            : /**                      BNFISINTNORM (K. Belabas)                 **/
    1244                 :            : /**                                                                **/
    1245                 :            : /********************************************************************/
    1246                 :            : struct sol_abs
    1247                 :            : {
    1248                 :            :   GEN rel; /* Primes PR[i] above a, expressed on generators of Cl(K) */
    1249                 :            :   GEN partrel; /* list of vectors, partrel[i] = rel[1..i] * u[1..i] */
    1250                 :            :   GEN cyc;     /* orders of generators of Cl(K) given in bnf */
    1251                 :            : 
    1252                 :            :   long *f;     /* f[i] = f(PR[i]/p), inertia degree */
    1253                 :            :   long *n;     /* a = prod p^{ n_p }. n[i]=n_p if PR[i] divides p */
    1254                 :            :   long *next;  /* index of first P above next p, 0 if p is last */
    1255                 :            :   long *S;     /* S[i] = n[i] - sum_{ 1<=k<=i } f[k]*u[k] */
    1256                 :            :   long *u;     /* We want principal ideals I = prod PR[i]^u[i] */
    1257                 :            :   GEN  normsol;/* lists of copies of the u[] which are solutions */
    1258                 :            : 
    1259                 :            :   long nPR;    /* length(T->rel) = #PR */
    1260                 :            :   long sindex, smax; /* current index in T->normsol; max. index */
    1261                 :            : };
    1262                 :            : 
    1263                 :            : /* u[1..i] has been filled. Norm(u) is correct.
    1264                 :            :  * Check relations in class group then save it. */
    1265                 :            : static void
    1266                 :      19695 : test_sol(struct sol_abs *T, long i)
    1267                 :            : {
    1268                 :            :   long k, l;
    1269                 :            :   GEN s;
    1270                 :            : 
    1271 [ +  + ][ +  + ]:      39390 :   if (T->partrel && !ZV_dvd(gel(T->partrel, i),  T->cyc)) return;
    1272         [ -  + ]:       8220 :   if (T->sindex == T->smax)
    1273                 :            :   { /* no more room in solution list: enlarge */
    1274                 :          0 :     long new_smax = T->smax << 1;
    1275                 :          0 :     GEN  new_normsol = new_chunk(new_smax+1);
    1276                 :            : 
    1277         [ #  # ]:          0 :     for (k=1; k<=T->smax; k++) gel(new_normsol,k) = gel(T->normsol,k);
    1278                 :          0 :     T->normsol = new_normsol; T->smax = new_smax;
    1279                 :            :   }
    1280                 :       8220 :   gel(T->normsol, ++T->sindex) = s = cgetg_copy(T->u, &l);
    1281         [ +  + ]:      35145 :   for (k=1; k <= i; k++) s[k] = T->u[k];
    1282         [ +  + ]:      14090 :   for (   ; k < l;  k++) s[k] = 0;
    1283         [ -  + ]:       8220 :   if (DEBUGLEVEL>2)
    1284                 :            :   {
    1285                 :          0 :     err_printf("sol = %Ps\n",s);
    1286         [ #  # ]:          0 :     if (T->partrel) err_printf("T->partrel = %Ps\n",T->partrel);
    1287                 :      19695 :     err_flush();
    1288                 :            :   }
    1289                 :            : }
    1290                 :            : /* partrel[i] <-- partrel[i-1] + u[i] * rel[i] */
    1291                 :            : static void
    1292                 :      15135 : fix_partrel(struct sol_abs *T, long i)
    1293                 :            : {
    1294                 :      15135 :   pari_sp av = avma;
    1295                 :      15135 :   GEN part1 = gel(T->partrel,i);
    1296                 :      15135 :   GEN part0 = gel(T->partrel,i-1);
    1297                 :      15135 :   GEN rel = gel(T->rel, i);
    1298                 :      15135 :   ulong u = T->u[i];
    1299                 :      15135 :   long k, l = lg(part1);
    1300         [ +  + ]:      44640 :   for (k=1; k < l; k++)
    1301                 :      29505 :     affii(addii(gel(part0,k), muliu(gel(rel,k), u)), gel(part1,k));
    1302                 :      15135 :   avma = av;
    1303                 :      15135 : }
    1304                 :            : 
    1305                 :            : /* Recursive loop. Suppose u[1..i] has been filled
    1306                 :            :  * Find possible solutions u such that, Norm(prod PR[i]^u[i]) = a, taking
    1307                 :            :  * into account:
    1308                 :            :  *  1) the relations in the class group if need be.
    1309                 :            :  *  2) the factorization of a. */
    1310                 :            : static void
    1311                 :      53310 : isintnorm_loop(struct sol_abs *T, long i)
    1312                 :            : {
    1313         [ +  + ]:      53310 :   if (T->S[i] == 0) /* sum u[i].f[i] = n[i], do another prime */
    1314                 :            :   {
    1315                 :      22730 :     long k, next = T->next[i];
    1316         [ +  + ]:      22730 :     if (next == 0) { test_sol(T, i); return; } /* no primes left */
    1317                 :            : 
    1318                 :            :     /* some primes left */
    1319         [ +  + ]:       9905 :     if (T->partrel) gaffect(gel(T->partrel,i), gel(T->partrel, next-1));
    1320         [ +  + ]:      14935 :     for (k=i+1; k < next; k++) T->u[k] = 0;
    1321                 :       9905 :     i = next-1;
    1322                 :            :   }
    1323 [ +  + ][ +  + ]:      30580 :   else if (i == T->next[i]-2 || i == T->nPR-1)
    1324                 :            :   { /* only one Prime left above prime; change prime, fix u[i+1] */
    1325                 :            :     long q;
    1326         [ +  + ]:      19445 :     if (T->S[i] % T->f[i+1]) return;
    1327                 :       9490 :     q = T->S[i] / T->f[i+1];
    1328                 :       9490 :     i++; T->u[i] = q;
    1329         [ +  + ]:       9490 :     if (T->partrel) fix_partrel(T,i);
    1330         [ +  + ]:       9490 :     if (T->next[i] == 0) { test_sol(T,i); return; }
    1331                 :            :   }
    1332                 :            : 
    1333                 :      23660 :   i++; T->u[i] = 0;
    1334         [ +  + ]:      23660 :   if (T->partrel) gaffect(gel(T->partrel,i-1), gel(T->partrel,i));
    1335         [ +  + ]:      23660 :   if (i == T->next[i-1])
    1336                 :            :   { /* change prime */
    1337 [ +  + ][ +  + ]:      21565 :     if (T->next[i] == i+1 || i == T->nPR) /* only one Prime above p */
    1338                 :            :     {
    1339                 :       7375 :       T->S[i] = 0;
    1340                 :       7375 :       T->u[i] = T->n[i] / T->f[i]; /* we already know this is exact */
    1341         [ +  + ]:       7375 :       if (T->partrel) fix_partrel(T, i);
    1342                 :            :     }
    1343                 :      21565 :     else T->S[i] = T->n[i];
    1344                 :            :   }
    1345                 :       2095 :   else T->S[i] = T->S[i-1]; /* same prime, different Prime */
    1346                 :            :   for(;;)
    1347                 :            :   {
    1348                 :      43965 :     isintnorm_loop(T, i);
    1349         [ +  + ]:      43965 :     T->S[i] -= T->f[i]; if (T->S[i] < 0) break;
    1350                 :      20305 :     T->u[i]++;
    1351         [ +  + ]:      20305 :     if (T->partrel) {
    1352                 :      15935 :       pari_sp av = avma;
    1353                 :      15935 :       gaffect(ZC_add(gel(T->partrel,i), gel(T->rel,i)), gel(T->partrel,i));
    1354                 :      15935 :       avma = av;
    1355                 :            :     }
    1356                 :      73615 :   }
    1357                 :            : }
    1358                 :            : 
    1359                 :            : static int
    1360                 :      15115 : get_sol_abs(struct sol_abs *T, GEN bnf, GEN fact, GEN *ptPR)
    1361                 :            : {
    1362                 :      15115 :   GEN nf = bnf_get_nf(bnf);
    1363                 :      15115 :   GEN P = gel(fact,1), E = gel(fact,2), PR;
    1364                 :      15115 :   long N = nf_get_degree(nf), nP = lg(P)-1, Ngen, max, nPR, i, j;
    1365                 :            : 
    1366                 :      15115 :   max = nP*N; /* upper bound for T->nPR */
    1367                 :      15115 :   T->f = new_chunk(max+1);
    1368                 :      15115 :   T->n = new_chunk(max+1);
    1369                 :      15115 :   T->next = new_chunk(max+1);
    1370                 :      15115 :   *ptPR = PR = cgetg(max+1, t_VEC); /* length to be fixed later */
    1371                 :            : 
    1372                 :      15115 :   nPR = 0;
    1373         [ +  + ]:      38715 :   for (i = 1; i <= nP; i++)
    1374                 :            :   {
    1375                 :      29370 :     GEN L = idealprimedec(nf, gel(P,i));
    1376                 :      29370 :     long lL = lg(L), gcd, k, v;
    1377                 :      29370 :     ulong vn = itou(gel(E,i));
    1378                 :            : 
    1379                 :            :     /* check that gcd_{P | p} f_P  divides  n_p */
    1380                 :      29370 :     gcd = pr_get_f(gel(L,1));
    1381 [ +  + ][ +  + ]:      29385 :     for (j=2; gcd > 1 && j < lL; j++) gcd = ugcd(gcd, pr_get_f(gel(L,j)));
    1382 [ +  + ][ +  + ]:      29370 :     if (gcd > 1 && vn % gcd)
    1383                 :            :     {
    1384         [ -  + ]:       5770 :       if (DEBUGLEVEL>2)
    1385                 :          0 :       { err_printf("gcd f_P  does not divide n_p\n"); err_flush(); }
    1386                 :       5770 :       return 0;
    1387                 :            :     }
    1388         [ +  + ]:      23600 :     v = (i==nP)? 0: nPR + lL;
    1389         [ +  + ]:      63920 :     for (k = 1; k < lL; k++)
    1390                 :            :     {
    1391                 :      40320 :       GEN pr = gel(L,k);
    1392                 :      40320 :       gel(PR, ++nPR) = pr;
    1393                 :      40320 :       T->f[nPR] = pr_get_f(pr) / gcd;
    1394                 :      40320 :       T->n[nPR] = vn / gcd;
    1395                 :      40320 :       T->next[nPR] = v;
    1396                 :            :     }
    1397                 :            :   }
    1398                 :       9345 :   T->nPR = nPR;
    1399                 :       9345 :   setlg(PR, nPR + 1);
    1400                 :            : 
    1401                 :       9345 :   T->u = cgetg(nPR+1, t_VECSMALL);
    1402                 :       9345 :   T->S = new_chunk(nPR+1);
    1403                 :       9345 :   T->cyc = bnf_get_cyc(bnf);
    1404                 :       9345 :   Ngen = lg(T->cyc)-1;
    1405         [ +  + ]:       9345 :   if (Ngen == 0)
    1406                 :        115 :     T->rel = T->partrel = NULL; /* trivial Cl(K), no relations to check */
    1407                 :            :   else
    1408                 :            :   {
    1409                 :       9230 :     int triv = 1;
    1410                 :       9230 :     T->partrel = new_chunk(nPR+1);
    1411                 :       9230 :     T->rel = new_chunk(nPR+1);
    1412         [ +  + ]:      41915 :     for (i=1; i <= nPR; i++)
    1413                 :            :     {
    1414                 :      32685 :       GEN c = isprincipal(bnf, gel(PR,i));
    1415                 :      32685 :       gel(T->rel,i) = c;
    1416 [ +  + ][ +  + ]:      32685 :       if (triv && !ZV_equal0(c)) triv = 0; /* non trivial relations in Cl(K)*/
    1417                 :            :     }
    1418                 :            :     /* triv = 1: all ideals dividing a are principal */
    1419         [ +  + ]:       9230 :     if (triv) T->rel = T->partrel = NULL;
    1420                 :            :   }
    1421         [ +  + ]:       9345 :   if (T->partrel)
    1422                 :            :   {
    1423                 :       7580 :     long B = ZV_max_lg(T->cyc) + 3;
    1424         [ +  + ]:      42850 :     for (i = 0; i <= nPR; i++)
    1425                 :            :     { /* T->partrel[0] also needs to be initialized */
    1426                 :      35270 :       GEN c = cgetg(Ngen+1, t_COL); gel(T->partrel,i) = c;
    1427         [ +  + ]:      95990 :       for (j=1; j<=Ngen; j++)
    1428                 :            :       {
    1429                 :      60720 :         GEN z = cgeti(B); gel(c,j) = z;
    1430                 :      60720 :         z[1] = evalsigne(0)|evallgefint(B);
    1431                 :            :       }
    1432                 :            :     }
    1433                 :            :   }
    1434                 :       9345 :   T->smax = 511;
    1435                 :       9345 :   T->normsol = new_chunk(T->smax+1);
    1436                 :       9345 :   T->S[0] = T->n[1];
    1437                 :       9345 :   T->next[0] = 1;
    1438                 :       9345 :   T->sindex = 0;
    1439                 :      15115 :   isintnorm_loop(T, 0); return 1;
    1440                 :            : }
    1441                 :            : 
    1442                 :            : /* Look for unit of norm -1. Return 1 if it exists and set *unit, 0 otherwise */
    1443                 :            : static long
    1444                 :       2175 : get_unit_1(GEN bnf, GEN *unit)
    1445                 :            : {
    1446                 :       2175 :   GEN v, nf = bnf_get_nf(bnf);
    1447                 :       2175 :   long i, n = nf_get_degree(nf);
    1448                 :            : 
    1449         [ -  + ]:       2175 :   if (DEBUGLEVEL > 2) err_printf("looking for a fundamental unit of norm -1\n");
    1450         [ +  + ]:       2175 :   if (odd(n)) { *unit = gen_m1; return 1; }
    1451                 :        825 :   v = nfsign_units(bnf, NULL, 0);
    1452         [ +  + ]:        890 :   for (i = 1; i < lg(v); i++)
    1453         [ +  + ]:        835 :     if ( Flv_sum( gel(v,i), 2) ) { *unit = gel(bnf_get_fu(bnf), i); return 1; }
    1454                 :       2175 :   return 0;
    1455                 :            : }
    1456                 :            : 
    1457                 :            : GEN
    1458                 :      15155 : bnfisintnormabs(GEN bnf, GEN a)
    1459                 :            : {
    1460                 :            :   struct sol_abs T;
    1461                 :            :   GEN nf, res, PR, F;
    1462                 :            :   long i;
    1463                 :            : 
    1464         [ +  + ]:      15155 :   if ((F = check_arith_all(a,"bnfisintnormabs")))
    1465                 :            :   {
    1466         [ -  + ]:         40 :     a = typ(a) == t_VEC? gel(a,1): factorback(F);
    1467         [ +  + ]:         40 :     if (signe(a) < 0) F = clean_Z_factor(F);
    1468                 :            :   }
    1469                 :      15155 :   bnf = checkbnf(bnf); nf = bnf_get_nf(bnf);
    1470         [ +  + ]:      15155 :   if (!signe(a)) return mkvec(gen_0);
    1471         [ +  + ]:      15145 :   if (is_pm1(a)) return mkvec(gen_1);
    1472         [ +  + ]:      15115 :   if (!F) F = absi_factor(a);
    1473         [ +  + ]:      15115 :   if (!get_sol_abs(&T, bnf, F, &PR)) return cgetg(1, t_VEC);
    1474                 :            :   /* |a| > 1 => T.nPR > 0 */
    1475                 :       9345 :   res = cgetg(T.sindex+1, t_VEC);
    1476         [ +  + ]:      17565 :   for (i=1; i<=T.sindex; i++)
    1477                 :            :   {
    1478                 :       8220 :     GEN x = vecsmall_to_col( gel(T.normsol,i) );
    1479                 :       8220 :     x = isprincipalfact(bnf, NULL, PR, x, nf_FORCE | nf_GEN_IF_PRINCIPAL);
    1480                 :       8220 :     gel(res,i) = coltoliftalg(nf, x); /* x solution, up to sign */
    1481                 :            :   }
    1482                 :      15155 :   return res;
    1483                 :            : }
    1484                 :            : 
    1485                 :            : /* z = bnfisintnormabs(bnf,a), sa = 1 or -1, return bnfisintnorm(bnf,sa*|a|) */
    1486                 :            : static GEN
    1487                 :      15160 : bnfisintnorm_i(GEN bnf, long sa, GEN z)
    1488                 :            : {
    1489                 :      15160 :   GEN nf = checknf(bnf), T = nf_get_pol(nf), unit = NULL;
    1490                 :      15160 :   long sNx, i, j, N = degpol(T), l = lg(z);
    1491                 :      15160 :   long norm_1 = 0; /* gcc -Wall */
    1492                 :            : 
    1493                 :            :   /* update z in place to get correct signs: multiply by unit of norm -1 if
    1494                 :            :    * it exists, otherwise delete solution with wrong sign */
    1495         [ +  + ]:      23425 :   for (i = j = 1; i < l; i++)
    1496                 :            :   {
    1497                 :       8265 :     GEN x = gel(z,i);
    1498                 :       8265 :     int xpol = (typ(x) == t_POL);
    1499                 :            : 
    1500         [ +  + ]:       8265 :     if (xpol) sNx = signe(ZX_resultant(T, Q_primpart(x)));
    1501 [ -  + ][ #  # ]:         40 :     else      sNx = gsigne(x) < 0 && odd(N) ? -1 : 1;
    1502         [ +  + ]:       8265 :     if (sNx != sa)
    1503                 :            :     {
    1504         [ +  + ]:       3665 :       if (! unit) norm_1 = get_unit_1(bnf, &unit);
    1505         [ +  + ]:       3665 :       if (!norm_1)
    1506                 :            :       {
    1507         [ -  + ]:         55 :         if (DEBUGLEVEL > 2) err_printf("%Ps eliminated because of sign\n",x);
    1508                 :         55 :         continue;
    1509                 :            :       }
    1510 [ +  + ][ +  + ]:       3610 :       if (xpol) x = (unit == gen_m1)? RgX_neg(x): RgXQ_mul(unit,x,T);
    1511         [ +  + ]:         20 :       else      x = (unit == gen_m1)? gneg(x): RgX_Rg_mul(unit,x);
    1512                 :            :     }
    1513                 :       8210 :     gel(z,j++) = x;
    1514                 :            :   }
    1515                 :      15160 :   setlg(z, j); return z;
    1516                 :            : }
    1517                 :            : GEN
    1518                 :      15135 : bnfisintnorm(GEN bnf, GEN a)
    1519                 :            : {
    1520                 :      15135 :   pari_sp av = avma;
    1521                 :      15135 :   GEN ne = bnfisintnormabs(bnf,a);
    1522      [ -  +  + ]:      15135 :   switch(typ(a))
    1523                 :            :   {
    1524                 :          0 :     case t_VEC: a = gel(a,1); break;
    1525                 :         40 :     case t_MAT: a = factorback(a); break;
    1526                 :            :   }
    1527                 :      15135 :   return gerepilecopy(av, bnfisintnorm_i(bnf,signe(a), ne));
    1528                 :            : }

Generated by: LCOV version 1.9