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