Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - modules - thue.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.14.0 lcov report (development 27783-affec94c65) Lines: 907 955 95.0 %
Date: 2022-07-07 07:34:25 Functions: 60 61 98.4 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : #include "pari.h"
      16             : #include "paripriv.h"
      17             : 
      18             : #define DEBUGLEVEL DEBUGLEVEL_thue
      19             : 
      20             : /********************************************************************/
      21             : /**                                                                **/
      22             : /**             THUE EQUATION SOLVER (G. Hanrot)                   **/
      23             : /**                                                                **/
      24             : /********************************************************************/
      25             : /* In all the forthcoming remarks, "paper" designs the paper "Thue Equations of
      26             :  * High Degree", by Yu. Bilu and G. Hanrot, J. Number Theory (1996). The
      27             :  * numbering of the constants corresponds to Hanrot's thesis rather than to the
      28             :  * paper. See also
      29             :  * "Solving Thue equations without the full unit group", Math. Comp. (2000) */
      30             : 
      31             : /* Check whether tnf is a valid structure */
      32             : static int
      33         791 : checktnf(GEN tnf)
      34             : {
      35         791 :   long l = lg(tnf);
      36             :   GEN v;
      37         791 :   if (typ(tnf)!=t_VEC || (l!=8 && l!=4)) return 0;
      38         385 :   v = gel(tnf,1);
      39         385 :   if (typ(v) != t_VEC || lg(v) != 4) return 0;
      40         385 :   if (l == 4) return 1; /* s=0 */
      41             : 
      42         182 :   (void)checkbnf(gel(tnf,2));
      43         182 :   return (typ(gel(tnf,3)) == t_COL
      44         182 :        && typ(gel(tnf,4)) == t_COL
      45         182 :        && typ(gel(tnf,5)) == t_MAT
      46         182 :        && typ(gel(tnf,6)) == t_MAT
      47         364 :        && typ(gel(tnf,7)) == t_VEC);
      48             : }
      49             : 
      50             : /* Compensates rounding errors for computation/display of the constants.
      51             :  * Round up if dir > 0, down otherwise */
      52             : static GEN
      53        8398 : myround(GEN x, long dir)
      54             : {
      55        8398 :   GEN eps = powis(stoi(dir > 0? 10: -10), -10);
      56        8398 :   return gmul(x, gadd(gen_1, eps));
      57             : }
      58             : 
      59             : /* v a t_VEC/t_VEC */
      60             : static GEN
      61        4065 : vecmax_shallow(GEN v) { return gel(v, vecindexmax(v)); }
      62             : 
      63             : static GEN
      64         233 : tnf_get_roots(GEN poly, long prec, long S, long T)
      65             : {
      66         233 :   GEN R0 = QX_complex_roots(poly, prec), R = cgetg(lg(R0), t_COL);
      67             :   long k;
      68             : 
      69         657 :   for (k=1; k<=S; k++) gel(R,k) = gel(R0,k);
      70             :   /* swap roots to get the usual order */
      71         416 :   for (k=1; k<=T; k++)
      72             :   {
      73         183 :     gel(R,k+S)  = gel(R0,2*k+S-1);
      74         183 :     gel(R,k+S+T)= gel(R0,2*k+S);
      75             :   }
      76         233 :   return R;
      77             : }
      78             : 
      79             : /* Computation of the logarithmic height of x (given by embeddings) */
      80             : static GEN
      81        4311 : LogHeight(GEN x, long prec)
      82             : {
      83        4311 :   pari_sp av = avma;
      84        4311 :   long i, n = lg(x)-1;
      85        4311 :   GEN LH = gen_1;
      86       18364 :   for (i=1; i<=n; i++)
      87             :   {
      88       14053 :     GEN t = gabs(gel(x,i), prec);
      89       14053 :     if (gcmpgs(t,1) > 0) LH = gmul(LH, t);
      90             :   }
      91        4311 :   return gerepileupto(av, gdivgu(glog(LH,prec), n));
      92             : }
      93             : 
      94             : /* |x|^(1/n), x t_INT */
      95             : static GEN
      96         113 : absisqrtn(GEN x, long n, long prec)
      97         113 : { GEN r = itor(x,prec); setabssign(r); return sqrtnr(r, n); }
      98             : 
      99             : static GEN
     100        4134 : get_emb(GEN x, GEN r)
     101             : {
     102        4134 :   long l = lg(r), i;
     103             :   GEN y;
     104             : 
     105        4134 :   if (typ(x) == t_INT) return const_col(l-1, x);
     106        4092 :   y = cgetg(l, t_COL);
     107       17327 :   for (i=1; i<l; i++)
     108             :   {
     109       13263 :     GEN e = poleval(x, gel(r,i));
     110       13263 :     if (gequal0(e) || (typ(e) != t_INT && precision(e) <= LOWDEFAULTPREC ))
     111          28 :       return NULL;
     112       13235 :     gel(y,i) = e;
     113             :   }
     114        4064 :   return y;
     115             : }
     116             : 
     117             : /* Computation of the conjugates (sigma_i(v_j)), and log. heights, of elts of v */
     118             : static GEN
     119         599 : Conj_LH(GEN v, GEN *H, GEN r, long prec)
     120             : {
     121         599 :   long j, l = lg(v);
     122         599 :   GEN e, M = cgetg(l,t_MAT);
     123             : 
     124         599 :   (*H) = cgetg(l,t_COL);
     125        4705 :   for (j = 1; j < l; j++)
     126             :   {
     127        4134 :     if (! (e = get_emb(gel(v,j), r)) ) return NULL; /* FAIL */
     128        4106 :     gel(M,j) = e;
     129        4106 :     gel(*H,j) = LogHeight(e, prec);
     130             :   }
     131         571 :   return M;
     132             : }
     133             : 
     134             : /* Computation of M, its inverse A and precision check (see paper) */
     135             : static GEN
     136         205 : T_A_Matrices(GEN MatFU, long r, GEN *eps5, long prec)
     137             : {
     138             :   GEN A, p1, m1, IntM, nia, eps3, eps2;
     139         205 :   long e = prec2nbits(prec);
     140             : 
     141         205 :   m1 = matslice(MatFU, 1,r, 1,r); /* minor order r */
     142         205 :   m1 = glog(gabs(m1, prec), prec); /* HIGH accuracy */
     143         205 :   A = RgM_inv(m1); if (!A) pari_err_PREC("thue");
     144         205 :   IntM = RgM_Rg_add_shallow(RgM_mul(A,m1), gen_m1);
     145         205 :   IntM = gabs(IntM, 0);
     146             : 
     147         205 :   eps2 = gadd(vecmax(IntM), real2n(-e, LOWDEFAULTPREC)); /* t_REAL */
     148         205 :   nia = vecmax(gabs(A, 0));
     149             : 
     150             :   /* Check for the precision in matrix inversion. See paper, Lemma 2.4.2. */
     151         205 :   p1 = addrr(mulur(r, gmul2n(nia, e)), eps2); /* t_REAL */
     152         205 :   if (expo(p1) < -2*r) pari_err_PREC("thue");
     153             : 
     154         205 :   p1 = addrr(mulur(r, gmul2n(nia,-e)), eps2);
     155         205 :   eps3 = mulrr(mulur(2*r*r,nia), p1);
     156         205 :   if (!signe(eps3))
     157          77 :     eps3 = real2n(expo(eps3), LOWDEFAULTPREC);
     158             :   else
     159         128 :     eps3 = myround(eps3, 1);
     160             : 
     161         205 :   if (DEBUGLEVEL>1) err_printf("epsilon_3 -> %Ps\n",eps3);
     162         205 :   *eps5 = mulur(r, eps3); return A;
     163             : }
     164             : 
     165             : /* find a few large primes such that p Z_K = P1 P2 P3 Q, with f(Pi/p) = 1
     166             :  * From x - \alpha y = \prod u_i^b_i we will deduce 3 equations in F_p
     167             :  * in check_prinfo. Eliminating x,y we get a stringent condition on (b_i). */
     168             : static GEN
     169         205 : get_prime_info(GEN bnf)
     170             : {
     171         205 :   long n = 1, e = gexpo(bnf_get_reg(bnf)), nbp = e < 20? 1: 2;
     172         205 :   GEN L = cgetg(nbp+1, t_VEC), nf = bnf_get_nf(bnf), fu = bnf_get_fu(bnf);
     173         205 :   GEN X = pol_x(nf_get_varn(nf));
     174             :   ulong p;
     175        3295 :   for(p = 2147483659UL; n <= nbp; p = unextprime(p+1))
     176             :   {
     177        3090 :     GEN PR, A, U, LP = idealprimedec_limit_f(bnf, utoipos(p), 1);
     178             :     long i;
     179        3090 :     if (lg(LP) < 4) continue;
     180         219 :     A = cgetg(5, t_VECSMALL);
     181         219 :     U = cgetg(4, t_VEC);
     182         219 :     PR = cgetg(4, t_VEC);
     183         876 :     for (i = 1; i <= 3; i++)
     184             :     {
     185         657 :       GEN modpr = zkmodprinit(nf, gel(LP,i));
     186         657 :       GEN a = nf_to_Fq(nf, X, modpr);
     187         657 :       GEN u = nfV_to_FqV(fu, nf, modpr);
     188         657 :       A[i] = itou(a);
     189         657 :       gel(U,i) = ZV_to_Flv(u,p);
     190         657 :       gel(PR,i) = modpr;
     191             :     }
     192         219 :     A[4] = p;
     193         219 :     gel(L,n++) = mkvec3(A,U,PR);
     194             :   }
     195         205 :   return L;
     196             : }
     197             : 
     198             : /* Performs basic computations concerning the equation.
     199             :  * Returns a "tnf" structure containing
     200             :  *  1) the polynomial
     201             :  *  2) the bnf (used to solve the norm equation)
     202             :  *  3) roots, with presumably enough precision
     203             :  *  4) The logarithmic heights of units
     204             :  *  5) The matrix of conjugates of units
     205             :  *  6) its inverse
     206             :  *  7) a few technical constants */
     207             : static GEN
     208         205 : inithue(GEN P, GEN bnf, long flag, long prec)
     209             : {
     210             :   GEN fu, MatFU, x0, tnf, tmp, gpmin, dP, csts, ALH, eps5, ro, c1, c2;
     211         205 :   GEN Ind = gen_1;
     212         205 :   long s,t, k,j, prec_roots, n = degpol(P);
     213             : 
     214         205 :   if (!bnf)
     215             :   {
     216         189 :     bnf = Buchall(P, nf_FORCE, maxss(prec, DEFAULTPREC));
     217         189 :     if (flag) (void)bnfcertify(bnf);
     218             :     else
     219         154 :       Ind = floorr(mulru(bnf_get_reg(bnf), 5));
     220             :   }
     221             : 
     222         205 :   nf_get_sign(bnf_get_nf(bnf), &s, &t);
     223         205 :   fu = bnf_get_fu(bnf);
     224         205 :   prec_roots = prec + nbits2extraprec(gexpo(Q_primpart(fu)));
     225             :   for(;;)
     226             :   {
     227         233 :     ro = tnf_get_roots(P, prec_roots, s, t);
     228         233 :     MatFU = Conj_LH(fu, &ALH, ro, prec);
     229         233 :     if (MatFU) break;
     230          28 :     prec_roots = precdbl(prec_roots);
     231          28 :     if (DEBUGLEVEL>1) pari_warn(warnprec, "inithue", prec_roots);
     232             :   }
     233             : 
     234         205 :   dP = ZX_deriv(P);
     235         205 :   c1 = NULL; /* min |P'(r_i)|, i <= s */
     236         587 :   for (k=1; k<=s; k++)
     237             :   {
     238         382 :     tmp = gabs(poleval(dP,gel(ro,k)),prec);
     239         382 :     if (!c1 || gcmp(tmp,c1) < 0) c1 = tmp;
     240             :   }
     241         205 :   c1 = gdiv(int2n(n-1), c1);
     242         205 :   c1 = gprec_w(myround(c1, 1), DEFAULTPREC);
     243             : 
     244         205 :   c2 = NULL; /* max |r_i - r_j|, i!=j */
     245         911 :   for (k=1; k<=n; k++)
     246        1699 :     for (j=k+1; j<=n; j++)
     247             :     {
     248         993 :       tmp = gabs(gsub(gel(ro,j),gel(ro,k)), prec);
     249         993 :       if (!c2 || gcmp(c2,tmp) > 0) c2 = tmp;
     250             :     }
     251         205 :   c2 = gprec_w(myround(c2, -1), DEFAULTPREC);
     252             : 
     253         205 :   if (t==0)
     254          78 :     x0 = real_1(DEFAULTPREC);
     255             :   else
     256             :   {
     257         127 :     gpmin = NULL; /* min |P'(r_i)|, i > s */
     258         289 :     for (k=1; k<=t; k++)
     259             :     {
     260         162 :       tmp = gabs(poleval(dP,gel(ro,s+k)), prec);
     261         162 :       if (!gpmin || gcmp(tmp,gpmin) < 0) gpmin = tmp;
     262             :     }
     263         127 :     gpmin = gprec_w(gpmin, DEFAULTPREC);
     264             : 
     265             :     /* Compute x0. See paper, Prop. 2.2.1 */
     266         127 :     x0 = gmul(gpmin, vecmax_shallow(gabs(imag_i(ro), prec)));
     267         127 :     x0 = sqrtnr(gdiv(int2n(n-1), x0), n);
     268             :   }
     269         205 :   if (DEBUGLEVEL>1)
     270           0 :     err_printf("c1 = %Ps\nc2 = %Ps\nIndice <= %Ps\n", c1, c2, Ind);
     271             : 
     272         205 :   ALH = gmul2n(ALH, 1);
     273         205 :   tnf = cgetg(8,t_VEC); csts = cgetg(9,t_VEC);
     274         205 :   gel(tnf,1) = P;
     275         205 :   gel(tnf,2) = bnf;
     276         205 :   gel(tnf,3) = ro;
     277         205 :   gel(tnf,4) = ALH;
     278         205 :   gel(tnf,5) = MatFU;
     279         205 :   gel(tnf,6) = T_A_Matrices(MatFU, s+t-1, &eps5, prec);
     280         205 :   gel(tnf,7) = csts;
     281         205 :   gel(csts,1) = c1; gel(csts,2) = c2;   gel(csts,3) = LogHeight(ro, prec);
     282         205 :   gel(csts,4) = x0; gel(csts,5) = eps5; gel(csts,6) = utoipos(prec);
     283         205 :   gel(csts,7) = Ind;
     284         205 :   gel(csts,8) = get_prime_info(bnf);
     285         205 :   return tnf;
     286             : }
     287             : 
     288             : typedef struct {
     289             :   GEN c10, c11, c13, c15, c91, bak, NE, Ind, hal, MatFU, divro, Hmu;
     290             :   GEN delta, lambda, inverrdelta, ro, Pi, Pi2;
     291             :   long r, iroot, deg;
     292             : } baker_s;
     293             : 
     294             : static void
     295        4113 : other_roots(long iroot, long *i1, long *i2)
     296             : {
     297        4113 :   switch (iroot) {
     298        1614 :     case 1: *i1=2; *i2=3; break;
     299        1176 :     case 2: *i1=1; *i2=3; break;
     300        1323 :    default: *i1=1; *i2=2; break;
     301             :   }
     302        4113 : }
     303             : /* low precision */
     304        4471 : static GEN abslog(GEN x) { return gabs(glog(gtofp(x,DEFAULTPREC),0), 0); }
     305             : 
     306             : /* Compute Baker's bound c9 and B_0, the bound for the b_i's. See Thm 2.3.1 */
     307             : static GEN
     308        3747 : Baker(baker_s *BS)
     309             : {
     310             :   GEN tmp, B0, hb0, c9, Indc11;
     311             :   long i1, i2;
     312             : 
     313        3747 :   other_roots(BS->iroot, &i1,&i2);
     314             :   /* Compute a bound for the h_0 */
     315        3747 :   hb0 = gadd(gmul2n(BS->hal,2), gmul2n(gadd(BS->Hmu,mplog2(DEFAULTPREC)), 1));
     316        3747 :   tmp = gmul(BS->divro, gdiv(gel(BS->NE,i1), gel(BS->NE,i2)));
     317        3747 :   tmp = gmax_shallow(gen_1, abslog(tmp));
     318        3747 :   hb0 = gmax_shallow(hb0, gdiv(tmp, BS->bak));
     319        3747 :   c9 = gmul(BS->c91,hb0);
     320        3747 :   c9 = gprec_w(myround(c9, 1), DEFAULTPREC);
     321        3747 :   Indc11 = rtor(mulir(BS->Ind,BS->c11), DEFAULTPREC);
     322             :   /* Compute B0 according to Lemma 2.3.3 */
     323        3747 :   B0 = mulir(shifti(BS->Ind,1),
     324             :              divrr(addrr(mulrr(c9,mplog(divrr(mulir(BS->Ind, c9),BS->c10))),
     325             :                          mplog(Indc11)), BS->c10));
     326        3747 :   B0 = gmax_shallow(B0, dbltor(2.71828183));
     327        3747 :   B0 = gmax_shallow(B0, mulrr(divir(BS->Ind, BS->c10),
     328             :                               mplog(divrr(Indc11, BS->Pi2))));
     329        3747 :   if (DEBUGLEVEL>1) {
     330           0 :     err_printf("  B0  = %Ps\n",B0);
     331           0 :     err_printf("  Baker = %Ps\n",c9);
     332             :   }
     333        3747 :   return B0;
     334             : }
     335             : 
     336             : /* || x d ||, x t_REAL, d t_INT */
     337             : static GEN
     338        8187 : errnum(GEN x, GEN d)
     339             : {
     340        8187 :   GEN dx = mulir(d, x), D = subri(dx, roundr(dx));
     341        8187 :   setabssign(D); return D;
     342             : }
     343             : 
     344             : /* Try to reduce the bound through continued fractions; see paper. */
     345             : static int
     346        4019 : CF_1stPass(GEN *B0, GEN kappa, baker_s *BS)
     347             : {
     348        4019 :   GEN a, b, q, ql, qd, l0, denbound = mulri(*B0, kappa);
     349             : 
     350        4019 :   if (cmprr(mulrr(dbltor(0.1),sqrr(denbound)), BS->inverrdelta) > 0)
     351           8 :     return -1;
     352             : 
     353        4011 :   q = denom_i( bestappr(BS->delta, denbound) );
     354        4011 :   qd = errnum(BS->delta, q);
     355        4011 :   ql = errnum(BS->lambda,q);
     356             : 
     357        4011 :   l0 = subrr(ql, addrr(mulrr(qd, *B0), divri(dbltor(0.1),kappa)));
     358        4011 :   if (signe(l0) <= 0) return 0;
     359             : 
     360        3353 :   if (BS->r > 1) {
     361        2597 :     a = BS->c15; b = BS->c13;
     362             :   }
     363             :   else {
     364         756 :     a = BS->c11; b = BS->c10;
     365         756 :     l0 = mulrr(l0, Pi2n(1, DEFAULTPREC));
     366             :   }
     367        3353 :   *B0 = divrr(mplog(divrr(mulir(q,a), l0)), b);
     368        3353 :   if (DEBUGLEVEL>1) err_printf("    B0 -> %Ps\n",*B0);
     369        3353 :   return 1;
     370             : }
     371             : 
     372             : static void
     373       10964 : get_B0Bx(baker_s *BS, GEN l0, GEN *B0, GEN *Bx)
     374             : {
     375       10964 :   GEN t = divrr(mulir(BS->Ind, BS->c15), l0);
     376       10964 :   *B0 = divrr(mulir(BS->Ind, mplog(t)), BS->c13);
     377       10964 :   *Bx = sqrtnr(shiftr(t,1), BS->deg);
     378       10964 : }
     379             : 
     380             : static int
     381       20020 : LLL_1stPass(GEN *pB0, GEN kappa, baker_s *BS, GEN *pBx)
     382             : {
     383       20020 :   GEN B0 = *pB0, Bx = *pBx, lllmat, C, l0, l1, triv;
     384             :   long e;
     385             : 
     386       20020 :   C = grndtoi(mulir(mulii(BS->Ind, kappa),
     387             :                     gpow(B0, dbltor(2.2), DEFAULTPREC)), NULL);
     388             : 
     389       20020 :   if (DEBUGLEVEL > 1) err_printf("C (bitsize) : %d\n", expi(C));
     390       20020 :   lllmat = matid(3);
     391       20020 :   if (cmpri(B0, BS->Ind) > 0)
     392             :   {
     393       14136 :     gcoeff(lllmat, 1, 1) = grndtoi(divri(B0, BS->Ind), NULL);
     394       14136 :     triv = shiftr(sqrr(B0), 1);
     395             :   }
     396             :   else
     397        5884 :     triv = addir(sqri(BS->Ind), sqrr(B0));
     398             : 
     399       20020 :   gcoeff(lllmat, 3, 1) = grndtoi(negr(mulir(C, BS->lambda)), &e);
     400       20020 :   if (e >= 0) return -1;
     401       20012 :   gcoeff(lllmat, 3, 2) = grndtoi(negr(mulir(C, BS->delta)), &e);
     402       20012 :   if (e >= 0) return -1;
     403       20012 :   gcoeff(lllmat, 3, 3) = C;
     404       20012 :   lllmat = ZM_lll(lllmat, 0.99, LLL_IM|LLL_INPLACE);
     405             : 
     406       20012 :   l0 = gnorml2(gel(lllmat,1));
     407       20012 :   l0 = subrr(divir(l0, dbltor(1.8262)), triv); /* delta = 0.99 */
     408       20012 :   if (signe(l0) <= 0) return 0;
     409             : 
     410       13848 :   l1 = shiftr(addri(shiftr(B0,1), BS->Ind), -1);
     411       13848 :   l0 = divri(subrr(sqrtr(l0), l1), C);
     412             : 
     413       13848 :   if (signe(l0) <= 0) return 0;
     414             : 
     415       10799 :   get_B0Bx(BS, l0, &B0, &Bx);
     416       10799 :   if (DEBUGLEVEL>=2)
     417             :   {
     418           0 :     err_printf("LLL_First_Pass successful\n");
     419           0 :     err_printf("B0 -> %Ps\n", B0);
     420           0 :     err_printf("x <= %Ps\n", Bx);
     421             :   }
     422       10799 :   *pB0 = B0; *pBx = Bx; return 1;
     423             : }
     424             : 
     425             : /* add solution (x,y) if not already known */
     426             : static void
     427         651 : add_sol(GEN *pS, GEN x, GEN y) { *pS = vec_append(*pS, mkvec2(x,y)); }
     428             : /* z = P(p,q), d = deg P, |z| = |rhs|. Check signs and (possibly)
     429             :  * add solutions (p,q), (-p,-q) */
     430             : static void
     431         161 : add_pm(GEN *pS, GEN p, GEN q, GEN z, long d, GEN rhs)
     432             : {
     433         161 :   if (signe(z) == signe(rhs))
     434             :   {
     435          91 :     add_sol(pS, p, q);
     436          91 :     if (!odd(d)) add_sol(pS, negi(p), negi(q));
     437             :   }
     438             :   else
     439          70 :     if (odd(d))  add_sol(pS, negi(p), negi(q));
     440         161 : }
     441             : 
     442             : /* Check whether a potential solution is a true solution. Return 0 if
     443             :  * truncation error (increase precision) */
     444             : static int
     445          77 : CheckSol(GEN *pS, GEN z1, GEN z2, GEN P, GEN rhs, GEN ro)
     446             : {
     447          77 :   GEN x, y, ro1 = gel(ro,1), ro2 = gel(ro,2);
     448             :   long e;
     449             : 
     450          77 :   y = grndtoi(real_i(gdiv(gsub(z2,z1), gsub(ro1,ro2))), &e);
     451          77 :   if (e > 0) return 0;
     452          77 :   if (!signe(y)) return 1; /* y = 0 taken care of in SmallSols */
     453          70 :   x = gadd(z1, gmul(ro1, y));
     454          70 :   x = grndtoi(real_i(x), &e);
     455          70 :   if (e > 0) return 0;
     456          70 :   if (e <= -13)
     457             :   { /* y != 0 and rhs != 0; check whether P(x,y) = rhs or P(-x,-y) = rhs */
     458          70 :     GEN z = ZX_Z_eval(ZX_rescale(P,y),x);
     459          70 :     if (absequalii(z, rhs)) add_pm(pS, x,y, z, degpol(P), rhs);
     460             :   }
     461          70 :   return 1;
     462             : }
     463             : 
     464             : static const long EXPO1 = 7;
     465             : static int
     466       29207 : round_to_b(GEN v, long B, long b, GEN Delta2, long i1, GEN L)
     467             : {
     468       29207 :   long i, l = lg(v);
     469       29207 :   if (!b)
     470             :   {
     471        2140 :     for (i = 1; i < l; i++)
     472             :     {
     473             :       long c;
     474        1860 :       if (i == i1)
     475        1041 :         c = 0;
     476             :       else
     477             :       {
     478         819 :         GEN d = gneg(gel(L,i));
     479             :         long e;
     480         819 :         d = grndtoi(d,&e);
     481         819 :         if (e > -EXPO1 || is_bigint(d)) return 0;
     482          42 :         c = itos(d); if (labs(c) > B) return 0;
     483             :       }
     484        1083 :       v[i] = c;
     485             :     }
     486             :   }
     487             :   else
     488             :   {
     489       56357 :     for (i = 1; i < l; i++)
     490             :     {
     491             :       long c;
     492       35840 :       if (i == i1)
     493       28074 :         c = b;
     494             :       else
     495             :       {
     496        7766 :         GEN d = gsub(gmulgs(gel(Delta2,i), b), gel(L,i));
     497             :         long e;
     498        7766 :         d = grndtoi(d,&e);
     499        7766 :         if (e > -EXPO1 || is_bigint(d)) return 0;
     500         133 :         c = itos(d); if (labs(c) > B) return 0;
     501             :       }
     502       28207 :       v[i] = c;
     503             :     }
     504             :   }
     505       20797 :   return 1;
     506             : }
     507             : 
     508             : /* mu \prod U[i]^b[i] */
     509             : static ulong
     510       62391 : Fl_factorback(ulong mu, GEN U, GEN b, ulong p)
     511             : {
     512       62391 :   long i, l = lg(U);
     513       62391 :   ulong r = mu;
     514      125307 :   for (i = 1; i < l; i++)
     515             :   {
     516       62916 :     long c = b[i];
     517       62916 :     ulong u = U[i];
     518       62916 :     if (!c) continue;
     519       61698 :     if (c < 0) { u = Fl_inv(u,p); c = -c; }
     520       61698 :     r = Fl_mul(r, Fl_powu(u,c,p), p);
     521             :   }
     522       62391 :   return r;
     523             : }
     524             : 
     525             : /* x - alpha y = \pm mu \prod \mu_i^{b_i}. Reduce mod 3 distinct primes of
     526             :  * degree 1 above the same p, and eliminate x,y => drastic conditions on b_i */
     527             : static int
     528       20797 : check_pr(GEN bi, GEN Lmu, GEN L)
     529             : {
     530       20797 :   GEN A = gel(L,1), U = gel(L,2);
     531       20797 :   ulong a = A[1], b = A[2], c = A[3], p = A[4];
     532       20797 :   ulong r = Fl_mul(Fl_sub(c,b,p), Fl_factorback(Lmu[1],gel(U,1),bi, p), p);
     533       20797 :   ulong s = Fl_mul(Fl_sub(a,c,p), Fl_factorback(Lmu[2],gel(U,2),bi, p), p);
     534       20797 :   ulong t = Fl_mul(Fl_sub(b,a,p), Fl_factorback(Lmu[3],gel(U,3),bi, p), p);
     535       20797 :   return Fl_add(Fl_add(r,s,p),t,p) == 0;
     536             : }
     537             : static int
     538       20797 : check_prinfo(GEN b, GEN Lmu, GEN prinfo)
     539             : {
     540             :   long i;
     541       20874 :   for (i = 1; i < lg(prinfo); i++)
     542       20797 :     if (!check_pr(b, gel(Lmu,i), gel(prinfo,i))) return 0;
     543          77 :   return 1;
     544             : }
     545             : /* For each possible value of b_i1, compute the b_i's
     546             : * and 2 conjugates of z = x - alpha y. Then check. */
     547             : static int
     548        1057 : TrySol(GEN *pS, GEN B0, long i1, GEN Delta2, GEN Lambda, GEN ro,
     549             :        GEN Lmu, GEN NE, GEN MatFU, GEN prinfo, GEN P, GEN rhs)
     550             : {
     551        1057 :   long bi1, i, B = itos(gceil(B0)), l = lg(Delta2);
     552        1057 :   GEN b = cgetg(l,t_VECSMALL), L = cgetg(l,t_VEC);
     553             : 
     554        2940 :   for (i = 1; i < l; i++)
     555             :   {
     556        1883 :     if (i == i1)
     557        1057 :       gel(L,i) = gen_0;
     558             :     else
     559             :     {
     560         826 :       GEN delta = gel(Delta2,i);
     561         826 :       gel(L, i) = gsub(gmul(delta,gel(Lambda,i1)), gel(Lambda,i));
     562             :     }
     563             :   }
     564       30264 :   for (bi1 = -B; bi1 <= B; bi1++)
     565             :   {
     566             :     GEN z1, z2;
     567       29207 :     if (!round_to_b(b, B, bi1, Delta2, i1, L)) continue;
     568       20797 :     if (!check_prinfo(b, Lmu, prinfo)) continue;
     569          77 :     z1 = gel(NE,1);
     570          77 :     z2 = gel(NE,2);
     571         182 :     for (i = 1; i < l; i++)
     572             :     {
     573         105 :       z1 = gmul(z1, gpowgs(gcoeff(MatFU,1,i), b[i]));
     574         105 :       z2 = gmul(z2, gpowgs(gcoeff(MatFU,2,i), b[i]));
     575             :     }
     576          77 :     if (!CheckSol(pS, z1,z2,P,rhs,ro)) return 0;
     577             :   }
     578        1057 :   return 1;
     579             : }
     580             : 
     581             : /* find q1,q2,q3 st q1 b + q2 c + q3 ~ 0 */
     582             : static GEN
     583         165 : GuessQi(GEN b, GEN c, GEN *eps)
     584             : {
     585         165 :   const long shift = 65;
     586             :   GEN Q, Lat;
     587             : 
     588         165 :   Lat = matid(3);
     589         165 :   gcoeff(Lat,3,1) = ground(gmul2n(b, shift));
     590         165 :   gcoeff(Lat,3,2) = ground(gmul2n(c, shift));
     591         165 :   gcoeff(Lat,3,3) = int2n(shift);
     592             : 
     593         165 :   Q = gel(lllint(Lat),1);
     594         165 :   if (gequal0(gel(Q,2))) return NULL; /* FAIL */
     595             : 
     596         165 :   *eps = mpadd(mpadd(gel(Q,3), mpmul(gel(Q,1),b)), mpmul(gel(Q,2),c));
     597         165 :   *eps = mpabs_shallow(*eps); return Q;
     598             : }
     599             : 
     600             : /* x a t_REAL */
     601             : static GEN
     602         349 : myfloor(GEN x) { return expo(x) > 30 ? ceil_safe(x): floorr(x); }
     603             : 
     604             : /* Check for not-so-small solutions. Return a t_REAL or NULL */
     605             : static GEN
     606         175 : MiddleSols(GEN *pS, GEN bound, GEN roo, GEN P, GEN rhs, long s, GEN c1)
     607             : {
     608             :   long j, k, nmax, d;
     609             :   GEN bndcf;
     610             : 
     611         175 :   if (expo(bound) < 0) return bound;
     612         161 :   d = degpol(P);
     613         161 :   bndcf = sqrtnr(shiftr(c1,1), d - 2);
     614         161 :   if (cmprr(bound, bndcf) < 0) return bound;
     615             :   /* divide by log2((1+sqrt(5))/2)
     616             :    * 1 + ==> ceil
     617             :    * 2 + ==> continued fraction is normalized if last entry is 1
     618             :    * 3 + ==> start at a0, not a1 */
     619         146 :   nmax = 3 + (long)(dbllog2(bound) * 1.44042009041256);
     620         146 :   bound = myfloor(bound);
     621             : 
     622         458 :   for (k = 1; k <= s; k++)
     623             :   {
     624         312 :     GEN ro = real_i(gel(roo,k)), t = gboundcf(ro, nmax);
     625             :     GEN pm1, qm1, p0, q0;
     626             : 
     627         312 :     pm1 = gen_0; p0 = gen_1;
     628         312 :     qm1 = gen_1; q0 = gen_0;
     629             : 
     630       14961 :     for (j = 1; j < lg(t); j++)
     631             :     {
     632             :       GEN p, q, z, Q, R;
     633             :       pari_sp av;
     634       14954 :       p = addii(mulii(p0, gel(t,j)), pm1); pm1 = p0; p0 = p;
     635       14954 :       q = addii(mulii(q0, gel(t,j)), qm1); qm1 = q0; q0 = q;
     636       14954 :       if (cmpii(q, bound) > 0) break;
     637       14649 :       if (DEBUGLEVEL >= 2) err_printf("Checking (+/- %Ps, +/- %Ps)\n",p, q);
     638       14649 :       av = avma;
     639       14649 :       z = ZX_Z_eval(ZX_rescale(P,q), p); /* = P(p/q) q^dep(P) */
     640       14649 :       Q = dvmdii(rhs, z, &R);
     641       14649 :       if (R != gen_0) { set_avma(av); continue; }
     642         196 :       setabssign(Q);
     643         196 :       if (Z_ispowerall(Q, d, &Q))
     644             :       {
     645          91 :         if (!is_pm1(Q)) { p = mulii(p, Q); q = mulii(q, Q); }
     646          91 :         add_pm(pS, p, q, z, d, rhs);
     647             :       }
     648             :     }
     649         312 :     if (j == lg(t))
     650             :     {
     651             :       long prec;
     652           7 :       if (j > nmax) pari_err_BUG("thue [short continued fraction]");
     653             :       /* the theoretical value is bit_prec = gexpo(ro)+1+log2(bound) */
     654           7 :       prec = precdbl(precision(ro));
     655           7 :       if (DEBUGLEVEL>1) pari_warn(warnprec,"thue",prec);
     656           7 :       roo = ZX_realroots_irred(P, prec);
     657           7 :       if (lg(roo)-1 != s) pari_err_BUG("thue [realroots]");
     658           7 :       k--;
     659             :     }
     660             :   }
     661         146 :   return bndcf;
     662             : }
     663             : 
     664             : static void
     665      579866 : check_y_root(GEN *pS, GEN P, GEN Y)
     666             : {
     667      579866 :   GEN r = nfrootsQ(P);
     668             :   long j;
     669      580279 :   for (j = 1; j < lg(r); j++)
     670         413 :     if (typ(gel(r,j)) == t_INT) add_sol(pS, gel(r,j), Y);
     671      579866 : }
     672             : 
     673             : static void
     674      289989 : check_y(GEN *pS, GEN P, GEN poly, GEN Y, GEN rhs)
     675             : {
     676      289989 :   long j, l = lg(poly);
     677      289989 :   GEN Yn = Y;
     678      289989 :   gel(P, l-1) = gel(poly, l-1);
     679     1159914 :   for (j = l-2; j >= 2; j--)
     680             :   {
     681      869925 :     gel(P,j) = mulii(Yn, gel(poly,j));
     682      869925 :     if (j > 2) Yn = mulii(Yn, Y);
     683             :   }
     684      289989 :   gel(P,2) = subii(gel(P,2), rhs); /* P = poly(Y/y)*y^deg(poly) - rhs */
     685      289989 :   check_y_root(pS, P, Y);
     686      289989 : }
     687             : 
     688             : /* Check for solutions under a small bound (see paper) */
     689             : static GEN
     690         203 : SmallSols(GEN S, GEN x3, GEN poly, GEN rhs)
     691             : {
     692         203 :   pari_sp av = avma;
     693             :   GEN X, P, rhs2;
     694         203 :   long j, l = lg(poly), n = degpol(poly);
     695             :   ulong y, By;
     696             : 
     697         203 :   x3 = myfloor(x3);
     698             : 
     699         203 :   if (DEBUGLEVEL>1) err_printf("* Checking for small solutions <= %Ps\n", x3);
     700         203 :   if (lgefint(x3) > 3)
     701           0 :     pari_err_OVERFLOW(stack_sprintf("thue (SmallSols): y <= %Ps", x3));
     702         203 :   By = itou(x3);
     703             :   /* y = 0 first: solve X^n = rhs */
     704         203 :   if (odd(n))
     705             :   {
     706         154 :     if (Z_ispowerall(absi_shallow(rhs), n, &X))
     707          28 :       add_sol(&S, signe(rhs) > 0? X: negi(X), gen_0);
     708             :   }
     709          49 :   else if (signe(rhs) > 0 && Z_ispowerall(rhs, n, &X))
     710             :   {
     711          21 :     add_sol(&S, X, gen_0);
     712          21 :     add_sol(&S, negi(X), gen_0);
     713             :   }
     714         203 :   rhs2 = shifti(rhs,1);
     715             :   /* y != 0 */
     716         203 :   P = cgetg(l, t_POL); P[1] = poly[1];
     717      290080 :   for (y = 1; y <= By; y++)
     718             :   {
     719      289877 :     pari_sp av2 = avma;
     720      289877 :     long lS = lg(S);
     721      289877 :     GEN Y = utoipos(y);
     722             :     /* try y */
     723      289877 :     check_y(&S, P, poly, Y, rhs);
     724             :     /* try -y */
     725      868973 :     for (j = l-2; j >= 2; j -= 2) togglesign( gel(P,j) );
     726      289877 :     if (j == 0) gel(P,2) = subii(gel(P,2), rhs2);
     727      289877 :     check_y_root(&S, P, utoineg(y));
     728      289877 :     if (lS == lg(S)) { set_avma(av2); continue; } /* no solution found */
     729             : 
     730         140 :     if (gc_needed(av,1))
     731             :     {
     732           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"SmallSols");
     733           0 :       gerepileall(av, 2, &S, &rhs2);
     734           0 :       P = cgetg(l, t_POL); P[1] = poly[1];
     735             :     }
     736             :   }
     737         203 :   return S;
     738             : }
     739             : 
     740             : /* Computes [x]! */
     741             : static double
     742         189 : fact(double x)
     743             : {
     744         189 :   double ft = 1.0;
     745         875 :   x = floor(x); while (x>1) { ft *= x; x--; }
     746         189 :   return ft ;
     747             : }
     748             : 
     749             : /* Compute all relevant constants needed to solve the equation P(x,y)=a given
     750             :  * the solutions of N_{K/Q}(x)=a (see inithue). */
     751             : GEN
     752         406 : thueinit(GEN pol, long flag, long prec)
     753             : {
     754         406 :   GEN POL, C, L, fa, tnf, bnf = NULL;
     755         406 :   pari_sp av = avma;
     756             :   long s, lfa, dpol;
     757             : 
     758         406 :   if (checktnf(pol)) { bnf = checkbnf(gel(pol,2)); pol = gel(pol,1); }
     759         406 :   if (typ(pol)!=t_POL) pari_err_TYPE("thueinit",pol);
     760         406 :   dpol = degpol(pol);
     761         406 :   if (dpol <= 0) pari_err_CONSTPOL("thueinit");
     762         399 :   RgX_check_ZX(pol, "thueinit");
     763         399 :   if (varn(pol)) { pol = leafcopy(pol); setvarn(pol, 0); }
     764             : 
     765         399 :   POL = Q_primitive_part(pol, &C);
     766         399 :   L = gen_1;
     767         399 :   fa = ZX_factor(POL); lfa = lgcols(fa);
     768         399 :   if (lfa > 2 || itos(gcoeff(fa,1,2)) > 1)
     769             :   { /* reducible polynomial, no need to reduce to the monic case */
     770          91 :     GEN P, Q, R, g, f = gcoeff(fa,1,1), E = gcoeff(fa,1,2);
     771          91 :     long e = itos(E);
     772          91 :     long vy = fetch_var();
     773          91 :     long va = fetch_var();
     774          91 :     long vb = fetch_var();
     775          91 :     C = C? ginv(C): gen_1;
     776          91 :     if (e != 1)
     777             :     {
     778          63 :       if (lfa == 2) {
     779          42 :         tnf = mkvec3(mkvec3(POL,C,L), thueinit(f, flag, prec), E);
     780          42 :         delete_var(); delete_var(); delete_var();
     781          42 :         return gerepilecopy(av, tnf);
     782             :       }
     783          21 :       P = gpowgs(f,e);
     784             :     }
     785             :     else
     786          28 :       P = f;
     787          49 :     g = RgX_div(POL, P);
     788          49 :     P = RgX_Rg_sub(RgX_homogenize(f, vy), pol_x(va));
     789          49 :     Q = RgX_Rg_sub(RgX_homogenize(g, vy), pol_x(vb));
     790          49 :     R = polresultant0(P, Q, -1, 0);
     791          49 :     tnf = mkvec3(mkvec3(POL,C,L), mkvecsmall4(degpol(f), e, va,vb),  R);
     792          49 :     delete_var(); delete_var(); delete_var();
     793          49 :     return gerepilecopy(av, tnf);
     794             :   }
     795             :   /* POL monic irreducible: POL(x) = C pol(x/L), L integer */
     796         308 :   POL = ZX_primitive_to_monic(POL, &L);
     797         308 :   C = gdiv(powiu(L, dpol), gel(pol, dpol+2));
     798         308 :   pol = POL;
     799             : 
     800         308 :   s = ZX_sturm_irred(pol);
     801         308 :   if (s)
     802             :   {
     803         189 :     long PREC, n = degpol(pol);
     804         189 :     double d, dr, dn = (double)n;
     805             : 
     806         189 :     if (dpol <= 2) pari_err_DOMAIN("thueinit", "P","=",pol,pol);
     807         189 :     dr = (double)((s+n-2)>>1); /* s+t-1 */
     808         189 :     d = dn*(dn-1)*(dn-2);
     809             :     /* Guess precision by approximating Baker's bound. The guess is most of
     810             :      * the time not sharp, ie 10 to 30 decimal digits above what is _really_
     811             :      * necessary. Note that the limiting step is the reduction. See paper. */
     812         189 :     PREC = nbits2prec((long)((5.83 + (dr+4)*5 + log(fact(dr+3)) + (dr+3)*log(dr+2) +
     813         378 :                      (dr+3)*log(d) + log(log(2*d*(dr+2))) + (dr+1))
     814         189 :                      /10.)*32+32);
     815             : 
     816         189 :     if (flag == 0) PREC = (long)(2.2 * PREC); /* Lazy, to be improved */
     817         189 :     if (PREC < prec) PREC = prec;
     818         189 :     if (DEBUGLEVEL >=2) err_printf("prec = %d\n", PREC);
     819             : 
     820             :     for (;;)
     821             :     {
     822         189 :       if (( tnf = inithue(pol, bnf, flag, PREC) )) break;
     823           0 :       PREC = precdbl(PREC);
     824           0 :       if (DEBUGLEVEL>1) pari_warn(warnprec,"thueinit",PREC);
     825           0 :       bnf = NULL; set_avma(av);
     826             :     }
     827             :   }
     828             :   else
     829             :   {
     830             :     GEN ro, c0;
     831             :     long k,l;
     832         119 :     if (!bnf)
     833             :     {
     834         119 :       bnf = gen_0;
     835         119 :       if (expi(ZX_disc(pol)) <= 50)
     836             :       {
     837          91 :         bnf = Buchall(pol, nf_FORCE, DEFAULTPREC);
     838          91 :         if (flag) (void)bnfcertify(bnf);
     839             :       }
     840             :     }
     841          91 :     ro = typ(bnf)==t_VEC? nf_get_roots(bnf_get_nf(bnf))
     842         119 :                         : QX_complex_roots(pol, DEFAULTPREC);
     843         119 :     l = lg(ro); c0 = imag_i(gel(ro,1));
     844         875 :     for (k = 2; k < l; k++) c0 = mulrr(c0, imag_i(gel(ro,k)));
     845         119 :     setsigne(c0,1); c0 = invr(c0); tnf = mkvec3(pol, bnf, c0);
     846             :   }
     847         308 :   gel(tnf,1) = mkvec3(gel(tnf,1), C, L);
     848         308 :   return gerepilecopy(av,tnf);
     849             : }
     850             : 
     851             : /* arg(t^2) / 2Pi; arg(t^2) = arg(t/conj(t)) */
     852             : static GEN
     853         632 : argsqr(GEN t, GEN Pi)
     854             : {
     855         632 :   GEN v, u = divrr(garg(t,0), Pi); /* in -1 < u <= 1 */
     856             :   /* reduce mod Z to -1/2 < u <= 1/2 */
     857         632 :   if (signe(u) > 0)
     858             :   {
     859         323 :     v = subrs(u,1); /* ]-1,0] */
     860         323 :     if (abscmprr(v,u) < 0) u = v;
     861             :   }
     862             :   else
     863             :   {
     864         309 :     v = addrs(u,1);/* ]0,1] */
     865         309 :     if (abscmprr(v,u) <= 0) u = v;
     866             :   }
     867         632 :   return u;
     868             : }
     869             : /* i1 != i2 */
     870             : static void
     871        3747 : init_get_B(long i1, long i2, GEN Delta2, GEN Lambda, GEN Deps5, baker_s *BS,
     872             :            long prec)
     873             : {
     874             :   GEN delta, lambda;
     875        3747 :   if (BS->r > 1)
     876             :   {
     877        3431 :     delta = gel(Delta2,i2);
     878        3431 :     lambda = gsub(gmul(delta,gel(Lambda,i1)), gel(Lambda,i2));
     879        3431 :     if (Deps5) BS->inverrdelta = divrr(Deps5, addsr(1,delta));
     880             :   }
     881             :   else
     882             :   { /* r == 1: i1 = s = t = 1; i2 = 2 */
     883         316 :     GEN fu = gel(BS->MatFU,1), ro = BS->ro, t;
     884             : 
     885         316 :     t = gel(fu,2);
     886         316 :     delta = argsqr(t, BS->Pi);
     887         316 :     if (Deps5) BS->inverrdelta = shiftr(gabs(t,prec), prec2nbits(prec)-1);
     888             : 
     889         316 :     t = gmul(gsub(gel(ro,1), gel(ro,2)), gel(BS->NE,3));
     890         316 :     lambda = argsqr(t, BS->Pi);
     891             :   }
     892        3747 :   BS->delta = delta;
     893        3747 :   BS->lambda = lambda;
     894        3747 : }
     895             : 
     896             : static GEN
     897        1065 : get_B0(long i1, GEN Delta2, GEN Lambda, GEN Deps5, long prec, baker_s *BS)
     898             : {
     899        1065 :   GEN B0 = Baker(BS);
     900        1065 :   long step = 0, i2 = (i1 == 1)? 2: 1;
     901             :   for(;;) /* i2 from 1 to r unless r = 1 [then i2 = 2] */
     902             :   {
     903        1065 :     init_get_B(i1,i2, Delta2,Lambda,Deps5, BS, prec);
     904             :     /* Reduce B0 as long as we make progress: newB0 < oldB0 - 0.1 */
     905             :     for (;;)
     906        2303 :     {
     907        3368 :       GEN oldB0 = B0, kappa = utoipos(10);
     908             :       long cf;
     909             : 
     910        4026 :       for (cf = 0; cf < 10; cf++, kappa = muliu(kappa,10))
     911             :       {
     912        4019 :         int res = CF_1stPass(&B0, kappa, BS);
     913        4019 :         if (res < 0) return NULL; /* prec problem */
     914        4011 :         if (res) break;
     915         658 :         if (DEBUGLEVEL>1) err_printf("CF failed. Increasing kappa\n");
     916             :       }
     917        3360 :       if (!step && cf == 10)
     918             :       { /* Semirational or totally rational case */
     919             :         GEN Q, ep, q, l0, denbound;
     920             : 
     921           0 :         if (! (Q = GuessQi(BS->delta, BS->lambda, &ep)) ) break;
     922             : 
     923           0 :         denbound = mpadd(B0, absi_shallow(gel(Q,1)));
     924           0 :         q = denom_i( bestappr(BS->delta, denbound) );
     925           0 :         l0 = subrr(errnum(BS->delta, q), ep);
     926           0 :         if (signe(l0) <= 0) break;
     927             : 
     928           0 :         B0 = divrr(mplog(divrr(mulir(gel(Q,2), BS->c15), l0)),  BS->c13);
     929           0 :         if (DEBUGLEVEL>1) err_printf("Semirat. reduction: B0 -> %Ps\n",B0);
     930             :       }
     931             :       /* if no progress, stop */
     932        3360 :       if (gcmp(oldB0, gadd(B0,dbltor(0.1))) <= 0)
     933        1057 :         return gmin_shallow(oldB0, B0);
     934        2303 :       else step++;
     935             :     }
     936           0 :     i2++; if (i2 == i1) i2++;
     937           0 :     if (i2 > BS->r) break;
     938             :   }
     939           0 :   pari_err_BUG("thue (totally rational case)");
     940             :   return NULL; /* LCOV_EXCL_LINE */
     941             : }
     942             : 
     943             : static GEN
     944        2682 : get_Bx_LLL(long i1, GEN Delta2, GEN Lambda, long prec, baker_s *BS)
     945             : {
     946        2682 :   GEN B0 = Baker(BS), Bx = NULL;
     947        2682 :   long step = 0, i2 = (i1 == 1)? 2: 1;
     948             :   for(;;) /* i2 from 1 to r unless r = 1 [then i2 = 2] */
     949             :   {
     950        2682 :     init_get_B(i1,i2, Delta2,Lambda,NULL, BS, prec);
     951        2682 :     if (DEBUGLEVEL>1) err_printf("  Entering LLL...\n");
     952             :     /* Reduce B0 as long as we make progress: newB0 < oldB0 - 0.1 */
     953             :     for (;;)
     954        8476 :     {
     955       11158 :       GEN oldBx = Bx, kappa = utoipos(10);
     956       11158 :       const long cfMAX = 10;
     957             :       long cf;
     958             : 
     959       20371 :       for (cf = 0; cf < cfMAX; cf++, kappa = muliu(kappa,10))
     960             :       {
     961       20020 :         int res = LLL_1stPass(&B0, kappa, BS, &Bx);
     962       20020 :         if (res < 0) return NULL;
     963       20012 :         if (res) break;
     964        9213 :         if (DEBUGLEVEL>1) err_printf("LLL failed. Increasing kappa\n");
     965             :       }
     966             : 
     967             :       /* FIXME: TO BE COMPLETED */
     968       11150 :       if (!step && cf == cfMAX)
     969             :       { /* Semirational or totally rational case */
     970             :         GEN Q, Q1, Q2, ep, q, l0, denbound;
     971             : 
     972         165 :         if (! (Q = GuessQi(BS->delta, BS->lambda, &ep)) ) break;
     973             : 
     974             :         /* Q[2] != 0 */
     975         165 :         Q1 = absi_shallow(gel(Q,1));
     976         165 :         Q2 = absi_shallow(gel(Q,2));
     977         165 :         denbound = gadd(mulri(B0, Q1), mulii(BS->Ind, Q2));
     978         165 :         q = denom_i( bestappr(BS->delta, denbound) );
     979         165 :         l0 = divri(subrr(errnum(BS->delta, q), ep), Q2);
     980         165 :         if (signe(l0) <= 0) break;
     981             : 
     982         165 :         get_B0Bx(BS, l0, &B0, &Bx);
     983         165 :         if (DEBUGLEVEL>1)
     984           0 :           err_printf("Semirat. reduction: B0 -> %Ps x <= %Ps\n",B0, Bx);
     985             :       }
     986             :       /* if no progress, stop */
     987       11150 :       if (oldBx && gcmp(oldBx, Bx) <= 0) return oldBx; else step++;
     988             :     }
     989           0 :     i2++; if (i2 == i1) i2++;
     990           0 :     if (i2 > BS->r) break;
     991             :   }
     992           0 :   pari_err_BUG("thue (totally rational case)");
     993             :   return NULL; /* LCOV_EXCL_LINE */
     994             : }
     995             : 
     996             : static GEN
     997         175 : LargeSols(GEN P, GEN tnf, GEN rhs, GEN ne)
     998             : {
     999         175 :   GEN S = NULL, Delta0, ro, ALH, bnf, nf, MatFU, A, csts, dP, Bx;
    1000             :   GEN c1,c2,c3,c4,c90,c91,c14, x0, x1, x2, x3, tmp, eps5, prinfo;
    1001             :   long iroot, ine, n, r, Prec, prec, s,t;
    1002             :   baker_s BS;
    1003         175 :   pari_sp av = avma;
    1004             : 
    1005         175 :   prec = 0; /*-Wall*/
    1006         175 :   bnf = NULL; /*-Wall*/
    1007         175 :   iroot = 1;
    1008         175 :   ine = 1;
    1009             : 
    1010         191 : START:
    1011         191 :   if (S) /* restart from precision problems */
    1012             :   {
    1013          16 :     S = gerepilecopy(av, S);
    1014          16 :     prec = precdbl(prec);
    1015          16 :     if (DEBUGLEVEL) pari_warn(warnprec,"thue",prec);
    1016          16 :     tnf = inithue(P, bnf, 0, prec);
    1017             :   }
    1018             :   else
    1019         175 :     S = cgetg(1, t_VEC);
    1020         191 :   bnf= gel(tnf,2);
    1021         191 :   nf = bnf_get_nf(bnf);
    1022         191 :   csts = gel(tnf,7);
    1023         191 :   nf_get_sign(nf, &s, &t);
    1024         191 :   BS.r = r = s+t-1; n = degpol(P);
    1025         191 :   ro     = gel(tnf,3);
    1026         191 :   ALH    = gel(tnf,4);
    1027         191 :   MatFU  = gel(tnf,5);
    1028         191 :   A      = gel(tnf,6);
    1029         191 :   c1     = mpmul(absi_shallow(rhs), gel(csts,1));
    1030         191 :   c2     = gel(csts,2);
    1031         191 :   BS.hal = gel(csts,3);
    1032         191 :   x0     = gel(csts,4);
    1033         191 :   eps5   = gel(csts,5);
    1034         191 :   Prec = itos(gel(csts,6));
    1035         191 :   BS.Ind = gel(csts,7);
    1036         191 :   BS.MatFU = MatFU;
    1037         191 :   BS.bak = muluu(n, (n-1)*(n-2)); /* safe */
    1038         191 :   BS.deg = n;
    1039         191 :   prinfo = gel(csts,8);
    1040             : 
    1041         191 :   if (t) x0 = gmul(x0, absisqrtn(rhs, n, Prec));
    1042         191 :   tmp = divrr(c1,c2);
    1043         191 :   c3 = mulrr(dbltor(1.39), tmp);
    1044         191 :   c4 = mulur(n-1, c3);
    1045         191 :   c14 = mulrr(c4, vecmax_shallow(RgM_sumcol(gabs(A,DEFAULTPREC))));
    1046             : 
    1047         191 :   x1 = gmax_shallow(x0, sqrtnr(shiftr(tmp,1),n));
    1048         191 :   x2 = gmax_shallow(x1, sqrtnr(mulur(10,c14), n));
    1049         191 :   x3 = gmax_shallow(x2, sqrtnr(shiftr(c14, EXPO1+1),n));
    1050         382 :   c90 = gmul(shiftr(mulur(18,mppi(DEFAULTPREC)), 5*(4+r)),
    1051         191 :                     gmul(gmul(mpfact(r+3), powiu(muliu(BS.bak,r+2), r+3)),
    1052         191 :                          glog(muliu(BS.bak,2*(r+2)),DEFAULTPREC)));
    1053             : 
    1054         191 :   dP = ZX_deriv(P);
    1055         191 :   Delta0 = RgM_sumcol(A);
    1056             : 
    1057         541 :   for (; iroot<=s; iroot++)
    1058             :   {
    1059         366 :     GEN Delta = Delta0, Delta2, D, Deps5, MatNE, Hmu, diffRo, c5, c7, ro0;
    1060             :     long i1, iroot1, iroot2, k;
    1061             : 
    1062         366 :     if (iroot <= r) Delta = RgC_add(Delta, RgC_Rg_mul(gel(A,iroot), stoi(-n)));
    1063         366 :     D = gabs(Delta,Prec); i1 = vecindexmax(D);
    1064         366 :     c5 = gel(D, i1);
    1065         366 :     Delta2 = RgC_Rg_div(Delta, gel(Delta, i1));
    1066         366 :     c5  = myround(gprec_w(c5,DEFAULTPREC), 1);
    1067         366 :     Deps5 = divrr(subrr(c5,eps5), eps5);
    1068         366 :     c7  = mulur(r,c5);
    1069         366 :     BS.c10 = divur(n,c7);
    1070         366 :     BS.c13 = divur(n,c5);
    1071         366 :     if (DEBUGLEVEL>1) {
    1072           0 :       err_printf("* real root no %ld/%ld\n", iroot,s);
    1073           0 :       err_printf("  c10 = %Ps\n",BS.c10);
    1074           0 :       err_printf("  c13 = %Ps\n",BS.c13);
    1075             :     }
    1076             : 
    1077         366 :     prec = Prec;
    1078             :     for (;;)
    1079             :     {
    1080         366 :       if (( MatNE = Conj_LH(ne, &Hmu, ro, prec) )) break;
    1081           0 :       prec = precdbl(prec);
    1082           0 :       if (DEBUGLEVEL>1) pari_warn(warnprec,"thue",prec);
    1083           0 :       ro = tnf_get_roots(P, prec, s, t);
    1084             :     }
    1085         366 :     ro0 = gel(ro,iroot);
    1086         366 :     BS.ro    = ro;
    1087         366 :     BS.iroot = iroot;
    1088         366 :     BS.Pi  = mppi(prec);
    1089         366 :     BS.Pi2 = Pi2n(1,prec);
    1090         366 :     diffRo = cgetg(r+1, t_VEC);
    1091        1090 :     for (k=1; k<=r; k++)
    1092             :     {
    1093         724 :       GEN z = gel(ro,k);
    1094         724 :       z = (k == iroot)? gdiv(rhs, poleval(dP, z)): gsub(ro0, z);
    1095         724 :       gel(diffRo,k) = gabs(z, prec);
    1096             :     }
    1097         366 :     other_roots(iroot, &iroot1,&iroot2);
    1098         366 :     BS.divro = gdiv(gsub(ro0, gel(ro,iroot2)), gsub(ro0, gel(ro,iroot1)));
    1099             :     /* Compute h_1....h_r */
    1100         366 :     c91 = c90;
    1101        1090 :     for (k=1; k<=r; k++)
    1102             :     {
    1103         724 :       GEN z = gdiv(gcoeff(MatFU,iroot1,k), gcoeff(MatFU,iroot2,k));
    1104         724 :       z = gmax_shallow(gen_1, abslog(z));
    1105         724 :       c91 = gmul(c91, gmax_shallow(gel(ALH,k), gdiv(z, BS.bak)));
    1106             :     }
    1107         366 :     BS.c91 = c91;
    1108             : 
    1109        4097 :     for (; ine<lg(ne); ine++)
    1110             :     {
    1111        3747 :       pari_sp av2 = avma;
    1112        3747 :       long lS = lg(S);
    1113             :       GEN Lambda, B0, c6, c8;
    1114        3747 :       GEN NE = gel(MatNE,ine), v = cgetg(r+1,t_COL);
    1115             : 
    1116        3747 :       if (DEBUGLEVEL>1) err_printf("  - norm sol. no %ld/%ld\n",ine,lg(ne)-1);
    1117       11583 :       for (k=1; k<=r; k++)
    1118             :       {
    1119        7836 :         GEN z = gdiv(gel(diffRo,k), gabs(gel(NE,k), prec));
    1120        7836 :         gel(v,k) = glog(z, prec);
    1121             :       }
    1122        3747 :       Lambda = RgM_RgC_mul(A,v);
    1123             : 
    1124        3747 :       c6 = addrr(dbltor(0.1), vecmax_shallow(gabs(Lambda,DEFAULTPREC)));
    1125        3747 :       c6 = myround(c6, 1);
    1126        3747 :       c8 = addrr(dbltor(1.23), mulur(r,c6));
    1127        3747 :       BS.c11= mulrr(shiftr(c3,1) , mpexp(divrr(mulur(n,c8),c7)));
    1128        3747 :       BS.c15= mulrr(shiftr(c14,1), mpexp(divrr(mulur(n,c6),c5)));
    1129        3747 :       BS.NE = NE;
    1130        3747 :       BS.Hmu = gel(Hmu,ine);
    1131        3747 :       if (is_pm1(BS.Ind))
    1132             :       {
    1133        1065 :         GEN mu = gel(ne,ine), Lmu = cgetg(lg(prinfo),t_VEC);
    1134             :         long i, j;
    1135             : 
    1136        2137 :         for (i = 1; i < lg(prinfo); i++)
    1137             :         {
    1138        1072 :           GEN v = gel(prinfo,i), PR = gel(v,3), L = cgetg(4, t_VECSMALL);
    1139        4288 :           for (j = 1; j <= 3; j++) L[j] = itou(nf_to_Fq(nf, mu, gel(PR,j)));
    1140        1072 :           gel(Lmu, i) = L;
    1141             :         }
    1142        2122 :         if (! (B0 = get_B0(i1, Delta2, Lambda, Deps5, prec, &BS)) ||
    1143        1057 :             !TrySol(&S, B0, i1, Delta2, Lambda, ro, Lmu, NE,MatFU,prinfo,
    1144             :                     P,rhs))
    1145          16 :           goto START;
    1146        1057 :         if (lg(S) == lS) set_avma(av2);
    1147             :       }
    1148             :       else
    1149             :       {
    1150        2682 :         if (! (Bx = get_Bx_LLL(i1, Delta2, Lambda, prec, &BS)) )
    1151           8 :            goto START;
    1152        2674 :         x3 = gerepileupto(av2, gmax_shallow(Bx, x3));
    1153             :       }
    1154             :     }
    1155         350 :     ine = 1;
    1156             :   }
    1157         175 :   x3 = gmax_shallow(x0, MiddleSols(&S, x3, ro, P, rhs, s, c1));
    1158         175 :   return SmallSols(S, x3, P, rhs);
    1159             : }
    1160             : 
    1161             : /* restrict to solutions (x,y) with L | x, replacing each by (x/L, y) */
    1162             : static GEN
    1163         336 : filter_sol_x(GEN S, GEN L)
    1164             : {
    1165             :   long i, k, l;
    1166         336 :   if (is_pm1(L)) return S;
    1167          84 :   l = lg(S); k = 1;
    1168         532 :   for (i = 1; i < l; i++)
    1169             :   {
    1170         448 :     GEN s = gel(S,i), r;
    1171         448 :     gel(s,1) = dvmdii(gel(s,1), L, &r);
    1172         448 :     if (r == gen_0) gel(S, k++) = s;
    1173             :   }
    1174          84 :   setlg(S, k); return S;
    1175             : }
    1176             : static GEN
    1177          70 : filter_sol_Z(GEN S)
    1178             : {
    1179          70 :   long i, k = 1, l = lg(S);
    1180         749 :   for (i = 1; i < l; i++)
    1181             :   {
    1182         679 :     GEN s = gel(S,i);
    1183         679 :     if (RgV_is_ZV(s)) gel(S, k++) = s;
    1184             :   }
    1185          70 :   setlg(S, k); return S;
    1186             : }
    1187             : 
    1188             : static GEN bnfisintnorm_i(GEN bnf, GEN a, long s, GEN z);
    1189             : static GEN
    1190         168 : tnf_get_Ind(GEN tnf) { return gmael(tnf,7,7); }
    1191             : static GEN
    1192         287 : tnf_get_bnf(GEN tnf) { return gel(tnf,2); }
    1193             : 
    1194             : static void
    1195           0 : maybe_warn(GEN bnf, GEN a, GEN Ind)
    1196             : {
    1197           0 :   if (!is_pm1(Ind) && !is_pm1(bnf_get_no(bnf)) && !is_pm1(a))
    1198           0 :     pari_warn(warner, "The result returned by 'thue' is conditional on the GRH");
    1199           0 : }
    1200             : /* true bnf; return solutions of Norm(x) = a mod U(K)^+ */
    1201             : static GEN
    1202         217 : get_ne(GEN bnf, GEN a, GEN fa, GEN Ind)
    1203             : {
    1204         217 :   if (DEBUGLEVEL) maybe_warn(bnf,a,Ind);
    1205         217 :   return bnfisintnorm_i(bnf, a, signe(a), bnfisintnormabs(bnf, mkvec2(a, fa)));
    1206             : }
    1207             : /* return solutions of |Norm(x)| = |a| mod U(K) */
    1208             : static GEN
    1209          28 : get_neabs(GEN bnf, GEN a, GEN Ind)
    1210             : {
    1211          28 :   if (DEBUGLEVEL) maybe_warn(bnf,a,Ind);
    1212          28 :   return bnfisintnormabs(bnf, a);
    1213             : }
    1214             : 
    1215             : /* Let P(z)=z^2+Bz+C, convert t=u+v*z (mod P) solution of norm equation N(t)=A
    1216             :  * to [x,y] = [u,-v] form: y^2P(x/y) = A */
    1217             : static GEN
    1218         679 : ne2_to_xy(GEN t)
    1219             : {
    1220             :   GEN u,v;
    1221         679 :   if (typ(t) != t_POL) { u = t; v = gen_0; }
    1222         651 :   else switch(degpol(t))
    1223             :   {
    1224           0 :     case -1: u = v = gen_0; break;
    1225          14 :     case 0: u = gel(t,2); v = gen_0; break;
    1226         637 :     default: u = gel(t,2); v = gneg(gel(t,3));
    1227             :   }
    1228         679 :   return mkvec2(u, v);
    1229             : }
    1230             : static GEN
    1231          70 : ne2V_to_xyV(GEN v)
    1232             : {
    1233             :   long i, l;
    1234          70 :   GEN w = cgetg_copy(v,&l);
    1235         749 :   for (i=1; i<l; i++) gel(w,i) = ne2_to_xy(gel(v,i));
    1236          70 :   return w;
    1237             : }
    1238             : 
    1239             : static GEN
    1240          21 : sol_0(void) { retmkvec( mkvec2(gen_0,gen_0) ); }
    1241             : static void
    1242         140 : sols_from_R(GEN Rab, GEN *pS, GEN P, GEN POL, GEN rhs)
    1243             : {
    1244         140 :   GEN ry = nfrootsQ(Rab);
    1245         140 :   long k, l = lg(ry);
    1246         280 :   for (k = 1; k < l; k++)
    1247         140 :     if (typ(gel(ry,k)) == t_INT) check_y(pS, P, POL, gel(ry,k), rhs);
    1248         140 : }
    1249             : static GEN
    1250          70 : absZ_factor_if_easy(GEN rhs, GEN x3)
    1251             : {
    1252             :   GEN F, U;
    1253          70 :   if (expi(rhs) < 150 || expo(x3) >= BITS_IN_LONG) return absZ_factor(rhs);
    1254           0 :   F = absZ_factor_limit_strict(rhs, 500000, &U);
    1255           0 :   return U? NULL: F;
    1256             : }
    1257             : /* Given a tnf structure as returned by thueinit, a RHS and
    1258             :  * optionally the solutions to the norm equation, returns the solutions to
    1259             :  * the Thue equation F(x,y)=a */
    1260             : GEN
    1261         385 : thue(GEN tnf, GEN rhs, GEN ne)
    1262             : {
    1263         385 :   pari_sp av = avma;
    1264             :   GEN POL, C, L, x3, S;
    1265             : 
    1266         385 :   if (typ(tnf) == t_POL) tnf = thueinit(tnf, 0, DEFAULTPREC);
    1267         385 :   if (!checktnf(tnf)) pari_err_TYPE("thue [please apply thueinit()]", tnf);
    1268         385 :   if (typ(rhs) != t_INT) pari_err_TYPE("thue",rhs);
    1269         385 :   if (ne && typ(ne) != t_VEC) pari_err_TYPE("thue",ne);
    1270             : 
    1271             :   /* solve P(x,y) = rhs <=> POL(L x, y) = C rhs, with POL monic in Z[X] */
    1272         385 :   POL = gel(tnf,1);
    1273         385 :   C = gel(POL,2); rhs = gmul(C, rhs);
    1274         385 :   if (typ(rhs) != t_INT) { set_avma(av); return cgetg(1, t_VEC); }
    1275         385 :   if (!signe(rhs))
    1276             :   {
    1277          28 :     GEN v = gel(tnf,2);
    1278          28 :     set_avma(av);
    1279             :     /* at least 2 irreducible factors, one of which has degree 1 */
    1280          28 :     if (typ(v) == t_VECSMALL && v[1] ==1)
    1281           7 :       pari_err_DOMAIN("thue","#sols","=",strtoGENstr("oo"),rhs);
    1282          21 :     return sol_0();
    1283             :   }
    1284         357 :   L = gel(POL,3);
    1285         357 :   POL = gel(POL,1);
    1286         357 :   if (lg(tnf) == 8)
    1287             :   {
    1288         175 :     if (!ne)
    1289             :     {
    1290         147 :       GEN F = absZ_factor(rhs);
    1291         147 :       ne = get_ne(tnf_get_bnf(tnf), rhs, F, tnf_get_Ind(tnf));
    1292             :     }
    1293         175 :     if (lg(ne) == 1) { set_avma(av); return cgetg(1, t_VEC); }
    1294         175 :     S = LargeSols(POL, tnf, rhs, ne);
    1295             :   }
    1296         182 :   else if (typ(gel(tnf,3)) == t_REAL)
    1297             :   { /* Case s=0. All solutions are "small". */
    1298         112 :     GEN bnf = tnf_get_bnf(tnf);
    1299         112 :     GEN c0 = gel(tnf,3), F;
    1300         112 :     x3 = sqrtnr(mulir(absi_shallow(rhs),c0), degpol(POL));
    1301         112 :     x3 = addrr(x3, dbltor(0.1)); /* guard from round-off errors */
    1302         112 :     S = cgetg(1,t_VEC);
    1303         112 :     if (!ne && typ(bnf) == t_VEC && expo(x3) > 10)
    1304             :     {
    1305          70 :       F = absZ_factor_if_easy(rhs, x3);
    1306          70 :       if (F) ne = get_ne(bnf, rhs, F, gen_1);
    1307             :     }
    1308         112 :     if (ne)
    1309             :     {
    1310          91 :       if (lg(ne) == 1) { set_avma(av); return cgetg(1,t_VEC); }
    1311          77 :       if (degpol(POL) == 2) /* quadratic imaginary */
    1312             :       {
    1313          70 :         GEN u = NULL;
    1314          70 :         long w = 2;
    1315          70 :         if (typ(bnf) == t_VEC)
    1316             :         {
    1317          56 :           u = bnf_get_tuU(bnf);
    1318          56 :           w = bnf_get_tuN(bnf);
    1319             :         }
    1320             :         else
    1321             :         {
    1322          14 :           GEN D = coredisc(ZX_disc(POL));
    1323          14 :           if (cmpis(D, -4) >= 0)
    1324             :           {
    1325          14 :             GEN F, T = quadpoly_i(D);
    1326          14 :             w = equalis(D, -4)? 4: 6;
    1327          14 :             setvarn(T, fetch_var_higher());
    1328          14 :             F = gcoeff(nffactor(POL, T), 1, 1);
    1329          14 :             u = gneg(lift_shallow(gel(F,2))); delete_var();
    1330             :           }
    1331             :         }
    1332          70 :         if (w == 4) /* u = I */
    1333          28 :           ne = shallowconcat(ne, RgXQV_RgXQ_mul(ne,u,POL));
    1334          42 :         else if (w == 6) /* u = j */
    1335             :         {
    1336          28 :           GEN u2 = RgXQ_sqr(u,POL);
    1337          28 :           ne = shallowconcat1(mkvec3(ne, RgXQV_RgXQ_mul(ne,u,POL),
    1338             :                                          RgXQV_RgXQ_mul(ne,u2,POL)));
    1339             :         }
    1340          70 :         S = ne2V_to_xyV(ne);
    1341          70 :         S = filter_sol_Z(S);
    1342          70 :         S = shallowconcat(S, RgV_neg(S));
    1343             :       }
    1344             :     }
    1345          98 :     if (lg(S) == 1) S = SmallSols(S, x3, POL, rhs);
    1346             :   }
    1347          70 :   else if (typ(gel(tnf,3)) == t_INT) /* reducible case, pure power*/
    1348             :   {
    1349          35 :     GEN bnf, ne1 = NULL, ne2 = NULL;
    1350          35 :     long e = itos( gel(tnf,3) );
    1351          35 :     if (!Z_ispowerall(rhs, e, &rhs)) { set_avma(av); return cgetg(1, t_VEC); }
    1352          28 :     tnf = gel(tnf,2);
    1353          28 :     bnf = tnf_get_bnf(tnf);
    1354          28 :     ne = get_neabs(bnf, rhs, lg(tnf)==8?tnf_get_Ind(tnf): gen_1);
    1355          28 :     ne1= bnfisintnorm_i(bnf,rhs,1,ne);
    1356          28 :     S = thue(tnf, rhs, ne1);
    1357          28 :     if (!odd(e) && lg(tnf)==8) /* if s=0, norms are positive */
    1358             :     {
    1359           7 :       ne2 = bnfisintnorm_i(bnf,rhs,-1,ne);
    1360           7 :       S = shallowconcat(S, thue(tnf, negi(rhs), ne2));
    1361             :     }
    1362             :   }
    1363             :   else /* other reducible cases */
    1364             :   { /* solve f^e * g = rhs, f irreducible factor of smallest degree */
    1365          35 :     GEN P, D, v = gel(tnf, 2), R = gel(tnf, 3);
    1366          35 :     long i, l, e = v[2], va = v[3], vb = v[4];
    1367          35 :     P = cgetg(lg(POL), t_POL); P[1] = POL[1];
    1368          35 :     D = divisors(rhs); l = lg(D);
    1369          35 :     S = cgetg(1,t_VEC);
    1370         217 :     for (i = 1; i < l; i++)
    1371             :     {
    1372         182 :       GEN Rab, df = gel(D,i), dg = gel(D,l-i); /* df*dg=|rhs| */
    1373         182 :       if (e > 1 && !Z_ispowerall(df, e, &df)) continue;
    1374             :       /* Rab: univariate polynomial in Z[Y], whose roots are the possible y. */
    1375             :       /* Here and below, Rab != 0 */
    1376          70 :       if (signe(rhs) < 0) dg = negi(dg); /* df*dg=rhs */
    1377          70 :       Rab = gsubst(gsubst(R, va, df), vb, dg);
    1378          70 :       sols_from_R(Rab, &S,P,POL,rhs);
    1379          70 :       Rab = gsubst(gsubst(R, va, negi(df)), vb, odd(e)? negi(dg): dg);
    1380          70 :       sols_from_R(Rab, &S,P,POL,rhs);
    1381             :     }
    1382             :   }
    1383         336 :   S = filter_sol_x(S, L);
    1384         336 :   S = gen_sort_uniq(S, (void*)lexcmp, cmp_nodata);
    1385         336 :   return gerepileupto(av, S);
    1386             : }
    1387             : 
    1388             : /********************************************************************/
    1389             : /**                                                                **/
    1390             : /**                      BNFISINTNORM (K. Belabas)                 **/
    1391             : /**                                                                **/
    1392             : /********************************************************************/
    1393             : struct sol_abs
    1394             : {
    1395             :   GEN rel; /* Primes PR[i] above a, expressed on generators of Cl(K) */
    1396             :   GEN partrel; /* list of vectors, partrel[i] = rel[1..i] * u[1..i] */
    1397             :   GEN cyc;     /* orders of generators of Cl(K) given in bnf */
    1398             : 
    1399             :   long *f;     /* f[i] = f(PR[i]/p), inertia degree */
    1400             :   long *n;     /* a = prod p^{ n_p }. n[i]=n_p if PR[i] divides p */
    1401             :   long *next;  /* index of first P above next p, 0 if p is last */
    1402             :   long *S;     /* S[i] = n[i] - sum_{ 1<=k<=i } f[k]*u[k] */
    1403             :   long *u;     /* We want principal ideals I = prod PR[i]^u[i] */
    1404             :   GEN  normsol;/* lists of copies of the u[] which are solutions */
    1405             : 
    1406             :   long nPR;    /* length(T->rel) = #PR */
    1407             :   long sindex, smax; /* current index in T->normsol; max. index */
    1408             : };
    1409             : 
    1410             : /* u[1..i] has been filled. Norm(u) is correct.
    1411             :  * Check relations in class group then save it. */
    1412             : static void
    1413       29911 : test_sol(struct sol_abs *T, long i)
    1414             : {
    1415             :   long k, l;
    1416             :   GEN s;
    1417             : 
    1418       29911 :   if (T->partrel && !ZV_dvd(gel(T->partrel, i),  T->cyc)) return;
    1419       13818 :   if (T->sindex == T->smax)
    1420             :   { /* no more room in solution list: enlarge */
    1421           0 :     long new_smax = T->smax << 1;
    1422           0 :     GEN  new_normsol = new_chunk(new_smax+1);
    1423             : 
    1424           0 :     for (k=1; k<=T->smax; k++) gel(new_normsol,k) = gel(T->normsol,k);
    1425           0 :     T->normsol = new_normsol; T->smax = new_smax;
    1426             :   }
    1427       13818 :   gel(T->normsol, ++T->sindex) = s = cgetg_copy(T->u, &l);
    1428       57694 :   for (k=1; k <= i; k++) s[k] = T->u[k];
    1429       25074 :   for (   ; k < l;  k++) s[k] = 0;
    1430       13818 :   if (DEBUGLEVEL>2)
    1431             :   {
    1432           0 :     err_printf("sol = %Ps\n",s);
    1433           0 :     if (T->partrel) err_printf("T->partrel = %Ps\n",T->partrel);
    1434             :   }
    1435             : }
    1436             : /* partrel[i] <-- partrel[i-1] + u[i] * rel[i] */
    1437             : static void
    1438       21210 : fix_partrel(struct sol_abs *T, long i)
    1439             : {
    1440       21210 :   pari_sp av = avma;
    1441       21210 :   GEN part1 = gel(T->partrel,i);
    1442       21210 :   GEN part0 = gel(T->partrel,i-1);
    1443       21210 :   GEN rel = gel(T->rel, i);
    1444       21210 :   ulong u = T->u[i];
    1445       21210 :   long k, l = lg(part1);
    1446       62538 :   for (k=1; k < l; k++)
    1447       41328 :     affii(addii(gel(part0,k), muliu(gel(rel,k), u)), gel(part1,k));
    1448       21210 :   set_avma(av);
    1449       21210 : }
    1450             : 
    1451             : /* Recursive loop. Suppose u[1..i] has been filled
    1452             :  * Find possible solutions u such that, Norm(prod PR[i]^u[i]) = a, taking
    1453             :  * into account:
    1454             :  *  1) the relations in the class group if need be.
    1455             :  *  2) the factorization of a. */
    1456             : static void
    1457       78799 : isintnorm_loop(struct sol_abs *T, long i)
    1458             : {
    1459       78799 :   if (T->S[i] == 0) /* sum u[i].f[i] = n[i], do another prime */
    1460             :   {
    1461       33551 :     long k, next = T->next[i];
    1462       33551 :     if (next == 0) { test_sol(T, i); return; } /* no primes left */
    1463             : 
    1464             :     /* some primes left */
    1465       13958 :     if (T->partrel) gaffect(gel(T->partrel,i), gel(T->partrel, next-1));
    1466       21049 :     for (k=i+1; k < next; k++) T->u[k] = 0;
    1467       13958 :     i = next-1;
    1468             :   }
    1469       45248 :   else if (i == T->next[i]-2 || i == T->nPR-1)
    1470             :   { /* only one Prime left above prime; change prime, fix u[i+1] */
    1471             :     long q;
    1472       28077 :     if (T->S[i] % T->f[i+1]) return;
    1473       14042 :     q = T->S[i] / T->f[i+1];
    1474       14042 :     i++; T->u[i] = q;
    1475       14042 :     if (T->partrel) fix_partrel(T,i);
    1476       14042 :     if (T->next[i] == 0) { test_sol(T,i); return; }
    1477             :   }
    1478             : 
    1479       34853 :   i++; T->u[i] = 0;
    1480       34853 :   if (T->partrel) gaffect(gel(T->partrel,i-1), gel(T->partrel,i));
    1481       34853 :   if (i == T->next[i-1])
    1482             :   { /* change prime */
    1483       30926 :     if (T->next[i] == i+1 || i == T->nPR) /* only one Prime above p */
    1484             :     {
    1485       10367 :       T->S[i] = 0;
    1486       10367 :       T->u[i] = T->n[i] / T->f[i]; /* we already know this is exact */
    1487       10367 :       if (T->partrel) fix_partrel(T, i);
    1488             :     }
    1489       20559 :     else T->S[i] = T->n[i];
    1490             :   }
    1491        3927 :   else T->S[i] = T->S[i-1]; /* same prime, different Prime */
    1492             :   for(;;)
    1493             :   {
    1494       65121 :     isintnorm_loop(T, i);
    1495       65121 :     T->S[i] -= T->f[i]; if (T->S[i] < 0) break;
    1496       30268 :     T->u[i]++;
    1497       30268 :     if (T->partrel) {
    1498       22344 :       pari_sp av = avma;
    1499       22344 :       gaffect(ZC_add(gel(T->partrel,i), gel(T->rel,i)), gel(T->partrel,i));
    1500       22344 :       set_avma(av);
    1501             :     }
    1502             :   }
    1503             : }
    1504             : 
    1505             : static int
    1506       21756 : get_sol_abs(struct sol_abs *T, GEN bnf, GEN nf, GEN fact, GEN *ptPR)
    1507             : {
    1508       21756 :   GEN P = gel(fact,1), E = gel(fact,2), PR;
    1509       21756 :   long N = nf_get_degree(nf), nP = lg(P)-1, Ngen, max, nPR, i, j;
    1510             : 
    1511       21756 :   max = nP*N; /* upper bound for T->nPR */
    1512       21756 :   T->f = new_chunk(max+1);
    1513       21756 :   T->n = new_chunk(max+1);
    1514       21756 :   T->next = new_chunk(max+1);
    1515       21756 :   *ptPR = PR = cgetg(max+1, t_COL); /* length to be fixed later */
    1516             : 
    1517       21756 :   nPR = 0;
    1518       55461 :   for (i = 1; i <= nP; i++)
    1519             :   {
    1520       41783 :     GEN L = idealprimedec(nf, gel(P,i));
    1521       41783 :     long lL = lg(L), gcd, k, v;
    1522       41783 :     ulong vn = itou(gel(E,i));
    1523             : 
    1524             :     /* check that gcd_{P | p} f_P  divides  n_p */
    1525       41783 :     gcd = pr_get_f(gel(L,1));
    1526       41804 :     for (j=2; gcd > 1 && j < lL; j++) gcd = ugcd(gcd, pr_get_f(gel(L,j)));
    1527       41783 :     if (gcd > 1 && vn % gcd)
    1528             :     {
    1529        8078 :       if (DEBUGLEVEL>2) err_printf("gcd f_P  does not divide n_p\n");
    1530        8078 :       return 0;
    1531             :     }
    1532       33705 :     v = (i==nP)? 0: nPR + lL;
    1533       92379 :     for (k = 1; k < lL; k++)
    1534             :     {
    1535       58674 :       GEN pr = gel(L,k);
    1536       58674 :       gel(PR, ++nPR) = pr;
    1537       58674 :       T->f[nPR] = pr_get_f(pr) / gcd;
    1538       58674 :       T->n[nPR] = vn / gcd;
    1539       58674 :       T->next[nPR] = v;
    1540             :     }
    1541             :   }
    1542       13678 :   T->nPR = nPR;
    1543       13678 :   setlg(PR, nPR + 1);
    1544             : 
    1545       13678 :   T->u = cgetg(nPR+1, t_VECSMALL);
    1546       13678 :   T->S = new_chunk(nPR+1);
    1547       13678 :   if (bnf) { T->cyc = bnf_get_cyc(bnf); Ngen = lg(T->cyc)-1; }
    1548         105 :   else { T->cyc = NULL; Ngen = 0; }
    1549       13678 :   if (Ngen == 0)
    1550         742 :     T->rel = T->partrel = NULL; /* trivial Cl(K), no relations to check */
    1551             :   else
    1552             :   {
    1553       12936 :     int triv = 1;
    1554       12936 :     T->partrel = new_chunk(nPR+1);
    1555       12936 :     T->rel = new_chunk(nPR+1);
    1556       58758 :     for (i=1; i <= nPR; i++)
    1557             :     {
    1558       45822 :       GEN c = isprincipal(bnf, gel(PR,i));
    1559       45822 :       gel(T->rel,i) = c;
    1560       45822 :       if (triv && !ZV_equal0(c)) triv = 0; /* non trivial relations in Cl(K)*/
    1561             :     }
    1562             :     /* triv = 1: all ideals dividing a are principal */
    1563       12936 :     if (triv) T->rel = T->partrel = NULL;
    1564             :   }
    1565       13678 :   if (T->partrel)
    1566             :   {
    1567       10619 :     long B = ZV_max_lg(T->cyc) + 3;
    1568       60025 :     for (i = 0; i <= nPR; i++)
    1569             :     { /* T->partrel[0] also needs to be initialized */
    1570       49406 :       GEN c = cgetg(Ngen+1, t_COL); gel(T->partrel,i) = c;
    1571      134442 :       for (j=1; j<=Ngen; j++)
    1572             :       {
    1573       85036 :         GEN z = cgeti(B); gel(c,j) = z;
    1574       85036 :         z[1] = evalsigne(0)|evallgefint(B);
    1575             :       }
    1576             :     }
    1577             :   }
    1578       13678 :   T->smax = 511;
    1579       13678 :   T->normsol = new_chunk(T->smax+1);
    1580       13678 :   T->S[0] = T->n[1];
    1581       13678 :   T->next[0] = 1;
    1582       13678 :   T->sindex = 0;
    1583       13678 :   isintnorm_loop(T, 0); return 1;
    1584             : }
    1585             : 
    1586             : /* Look for unit of norm -1. Return 1 if it exists and set *unit, 0 otherwise */
    1587             : static long
    1588        3178 : get_unit_1(GEN bnf, GEN *unit)
    1589             : {
    1590        3178 :   GEN v, nf = bnf_get_nf(bnf);
    1591        3178 :   long i, n = nf_get_degree(nf);
    1592             : 
    1593        3178 :   if (DEBUGLEVEL > 2) err_printf("looking for a fundamental unit of norm -1\n");
    1594        3178 :   if (odd(n)) { *unit = gen_m1; return 1; }
    1595        1232 :   v = nfsign_fu(bnf, NULL);
    1596        1323 :   for (i = 1; i < lg(v); i++)
    1597        1246 :     if ( Flv_sum( gel(v,i), 2) ) { *unit = gel(bnf_get_fu(bnf), i); return 1; }
    1598          77 :   return 0;
    1599             : }
    1600             : 
    1601             : GEN
    1602       21734 : bnfisintnormabs(GEN bnf, GEN a)
    1603             : {
    1604             :   struct sol_abs T;
    1605             :   GEN nf, res, PR, F;
    1606             :   long i;
    1607             : 
    1608       21734 :   if ((F = check_arith_all(a,"bnfisintnormabs")))
    1609             :   {
    1610         693 :     a = typ(a) == t_VEC? gel(a,1): factorback(F);
    1611         693 :     if (signe(a) < 0) F = clean_Z_factor(F);
    1612             :   }
    1613       21735 :   nf = bnf_get_nf(bnf);
    1614       21735 :   if (!signe(a)) return mkvec(gen_0);
    1615       21721 :   if (is_pm1(a)) return mkvec(gen_1);
    1616       21651 :   if (!F) F = absZ_factor(a);
    1617       21651 :   if (!get_sol_abs(&T, bnf, nf, F, &PR)) return cgetg(1, t_VEC);
    1618             :   /* |a| > 1 => T.nPR > 0 */
    1619       13573 :   res = cgetg(T.sindex+1, t_VEC);
    1620       27069 :   for (i=1; i<=T.sindex; i++)
    1621             :   {
    1622       13496 :     GEN x = vecsmall_to_col( gel(T.normsol,i) );
    1623       13496 :     x = isprincipalfact(bnf, NULL, PR, x, nf_FORCE | nf_GEN_IF_PRINCIPAL);
    1624       13496 :     gel(res,i) = nf_to_scalar_or_alg(nf, x); /* x solution, up to sign */
    1625             :   }
    1626       13573 :   return res;
    1627             : }
    1628             : 
    1629             : /* true nf */
    1630             : GEN
    1631         126 : ideals_by_norm(GEN nf, GEN a)
    1632             : {
    1633             :   struct sol_abs T;
    1634             :   GEN res, PR, F;
    1635             :   long i;
    1636             : 
    1637         126 :   if ((F = check_arith_pos(a,"ideals_by_norm")))
    1638           0 :     a = typ(a) == t_VEC? gel(a,1): factorback(F);
    1639         126 :   if (is_pm1(a)) return mkvec(trivial_fact());
    1640         105 :   if (!F) F = absZ_factor(a);
    1641         105 :   if (!get_sol_abs(&T, NULL, nf, F, &PR)) return cgetg(1, t_VEC);
    1642         105 :   res = cgetg(T.sindex+1, t_VEC);
    1643         427 :   for (i=1; i<=T.sindex; i++)
    1644             :   {
    1645         322 :     GEN x = vecsmall_to_col( gel(T.normsol,i) );
    1646         322 :     gel(res,i) = famat_remove_trivial(mkmat2(PR, x));
    1647             :   }
    1648         105 :   return res;
    1649             : }
    1650             : 
    1651             : /* true bnf; z = bnfisintnormabs(bnf,a), sa = 1 or -1,
    1652             :  * return bnfisintnorm(bnf,sa*|a|) */
    1653             : static GEN
    1654       21742 : bnfisintnorm_i(GEN bnf, GEN a, long sa, GEN z)
    1655             : {
    1656       21742 :   GEN nf = bnf_get_nf(bnf), T = nf_get_pol(nf), f = nf_get_index(nf), unit=NULL;
    1657       21742 :   GEN Tp, A = signe(a) == sa? a: negi(a);
    1658       21742 :   long sNx, i, j, N = degpol(T), l = lg(z);
    1659       21742 :   long norm_1 = 0; /* gcc -Wall */
    1660       21742 :   ulong p, Ap = 0; /* gcc -Wall */
    1661             :   forprime_t S;
    1662       21742 :   if (!signe(a)) return z;
    1663       21728 :   u_forprime_init(&S,3,ULONG_MAX);
    1664       30324 :   while((p = u_forprime_next(&S)))
    1665       30324 :     if (umodiu(f,p)) { Ap = umodiu(A,p); if (Ap) break; }
    1666       21728 :   Tp = ZX_to_Flx(T,p);
    1667             :   /* p > 2 doesn't divide A nor Q_denom(z in Z_K)*/
    1668             : 
    1669             :   /* update z in place to get correct signs: multiply by unit of norm -1 if
    1670             :    * it exists, otherwise delete solution with wrong sign */
    1671       35301 :   for (i = j = 1; i < l; i++)
    1672             :   {
    1673       13573 :     GEN x = gel(z,i);
    1674       13573 :     int xpol = (typ(x) == t_POL);
    1675             : 
    1676       13573 :     if (xpol)
    1677             :     {
    1678       12936 :       GEN dx, y = Q_remove_denom(x,&dx);
    1679       12936 :       ulong Np = Flx_resultant(Tp, ZX_to_Flx(y,p), p);
    1680       12936 :       ulong dA = dx? Fl_mul(Ap, Fl_powu(umodiu(dx,p), N, p), p): Ap;
    1681             :       /* Nx = Res(T,y) / dx^N = A or -A. Check mod p */
    1682       12936 :       sNx = dA == Np? sa: -sa;
    1683             :     }
    1684             :     else
    1685         637 :       sNx = gsigne(x) < 0 && odd(N) ? -1 : 1;
    1686       13573 :     if (sNx != sa)
    1687             :     {
    1688        5173 :       if (! unit) norm_1 = get_unit_1(bnf, &unit);
    1689        5173 :       if (!norm_1)
    1690             :       {
    1691          77 :         if (DEBUGLEVEL > 2) err_printf("%Ps eliminated because of sign\n",x);
    1692          77 :         continue;
    1693             :       }
    1694        5096 :       if (xpol) x = (unit == gen_m1)? RgX_neg(x): RgXQ_mul(unit,x,T);
    1695         182 :       else      x = (unit == gen_m1)? gneg(x): RgX_Rg_mul(unit,x);
    1696             :     }
    1697       13496 :     gel(z,j++) = x;
    1698             :   }
    1699       21728 :   setlg(z, j); return z;
    1700             : }
    1701             : GEN
    1702       21489 : bnfisintnorm(GEN bnf, GEN a)
    1703             : {
    1704       21489 :   pari_sp av = avma;
    1705             :   GEN ne;
    1706       21489 :   bnf = checkbnf(bnf);
    1707       21489 :   ne = bnfisintnormabs(bnf,a);
    1708       21490 :   switch(typ(a))
    1709             :   {
    1710         476 :     case t_VEC: a = gel(a,1); break;
    1711           0 :     case t_MAT: a = factorback(a); break;
    1712             :   }
    1713       21490 :   return gerepilecopy(av, bnfisintnorm_i(bnf, a, signe(a), ne));
    1714             : }

Generated by: LCOV version 1.13