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 8486 : myround(GEN x, long dir)
54 : {
55 8486 : GEN eps = powis(stoi(dir > 0? 10: -10), -10);
56 8486 : return gmul(x, gadd(gen_1, eps));
57 : }
58 :
59 : /* v a t_VEC/t_VEC */
60 : static GEN
61 4095 : vecmax_shallow(GEN v) { return gel(v, vecindexmax(v)); }
62 :
63 : static GEN
64 240 : tnf_get_roots(GEN poly, long prec, long S, long T)
65 : {
66 240 : GEN R0 = QX_complex_roots(poly, prec), R = cgetg(lg(R0), t_COL);
67 : long k;
68 :
69 685 : for (k=1; k<=S; k++) gel(R,k) = gel(R0,k);
70 : /* swap roots to get the usual order */
71 423 : 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 240 : return R;
77 : }
78 :
79 : /* Computation of the logarithmic height of x (given by embeddings) */
80 : static GEN
81 4391 : LogHeight(GEN x, long prec)
82 : {
83 4391 : pari_sp av = avma;
84 4391 : long i, n = lg(x)-1;
85 4391 : GEN LH = gen_1;
86 18684 : for (i=1; i<=n; i++)
87 : {
88 14293 : GEN t = gabs(gel(x,i), prec);
89 14293 : if (gcmpgs(t,1) > 0) LH = gmul(LH, t);
90 : }
91 4391 : 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 4205 : get_emb(GEN x, GEN r)
101 : {
102 4205 : long l = lg(r), i;
103 : GEN y;
104 :
105 4205 : if (typ(x) == t_INT) return const_col(l-1, x);
106 4163 : y = cgetg(l, t_COL);
107 17613 : for (i=1; i<l; i++)
108 : {
109 13477 : GEN e = poleval(x, gel(r,i));
110 13477 : if (gequal0(e) || (typ(e) != t_INT && precision(e) <= LOWDEFAULTPREC ))
111 27 : return NULL;
112 13450 : gel(y,i) = e;
113 : }
114 4136 : return y;
115 : }
116 :
117 : /* Computation of the conjugates (sigma_i(v_j)), and log. heights, of elts of v */
118 : static GEN
119 628 : Conj_LH(GEN v, GEN *H, GEN r, long prec)
120 : {
121 628 : long j, l = lg(v);
122 628 : GEN e, M = cgetg(l,t_MAT);
123 :
124 628 : (*H) = cgetg(l,t_COL);
125 4806 : for (j = 1; j < l; j++)
126 : {
127 4205 : if (! (e = get_emb(gel(v,j), r)) ) return NULL; /* FAIL */
128 4178 : gel(M,j) = e;
129 4178 : gel(*H,j) = LogHeight(e, prec);
130 : }
131 601 : return M;
132 : }
133 :
134 : /* Computation of M, its inverse A and precision check (see paper) */
135 : static GEN
136 213 : T_A_Matrices(GEN MatFU, long r, GEN *eps5, long prec)
137 : {
138 : GEN A, p1, m1, IntM, nia, eps3, eps2;
139 213 : long e = prec2nbits(prec);
140 :
141 213 : m1 = matslice(MatFU, 1,r, 1,r); /* minor order r */
142 213 : m1 = glog(gabs(m1, prec), prec); /* HIGH accuracy */
143 213 : A = RgM_inv(m1); if (!A) pari_err_PREC("thue");
144 213 : IntM = RgM_Rg_add_shallow(RgM_mul(A,m1), gen_m1);
145 213 : IntM = gabs(IntM, 0);
146 :
147 213 : eps2 = gadd(vecmax(IntM), real2n(-e, LOWDEFAULTPREC)); /* t_REAL */
148 213 : nia = vecmax(gabs(A, 0));
149 :
150 : /* Check for the precision in matrix inversion. See paper, Lemma 2.4.2. */
151 213 : p1 = addrr(mulur(r, gmul2n(nia, e)), eps2); /* t_REAL */
152 213 : if (expo(p1) < -2*r) pari_err_PREC("thue");
153 :
154 213 : p1 = addrr(mulur(r, gmul2n(nia,-e)), eps2);
155 213 : eps3 = mulrr(mulur(2*r*r,nia), p1);
156 213 : if (!signe(eps3))
157 79 : eps3 = real2n(expo(eps3), LOWDEFAULTPREC);
158 : else
159 134 : eps3 = myround(eps3, 1);
160 :
161 213 : if (DEBUGLEVEL>1) err_printf("epsilon_3 -> %Ps\n",eps3);
162 213 : *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 213 : get_prime_info(GEN bnf)
170 : {
171 213 : long n = 1, e = gexpo(bnf_get_reg(bnf)), nbp = e < 20? 1: 2;
172 213 : GEN L = cgetg(nbp+1, t_VEC), nf = bnf_get_nf(bnf), fu = bnf_get_fu(bnf);
173 213 : GEN X = pol_x(nf_get_varn(nf));
174 : ulong p;
175 3417 : for(p = 2147483659UL; n <= nbp; p = unextprime(p+1))
176 : {
177 3204 : GEN PR, A, U, LP = idealprimedec_limit_f(bnf, utoipos(p), 1);
178 : long i;
179 3204 : if (lg(LP) < 4) continue;
180 227 : A = cgetg(5, t_VECSMALL);
181 227 : U = cgetg(4, t_VEC);
182 227 : PR = cgetg(4, t_VEC);
183 908 : for (i = 1; i <= 3; i++)
184 : {
185 681 : GEN modpr = zkmodprinit(nf, gel(LP,i));
186 681 : GEN a = nf_to_Fq(nf, X, modpr);
187 681 : GEN u = nfV_to_FqV(fu, nf, modpr);
188 681 : A[i] = itou(a);
189 681 : gel(U,i) = ZV_to_Flv(u,p);
190 681 : gel(PR,i) = modpr;
191 : }
192 227 : A[4] = p;
193 227 : gel(L,n++) = mkvec3(A,U,PR);
194 : }
195 213 : 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 213 : 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 213 : GEN Ind = gen_1;
212 213 : long s,t, k,j, prec_roots, n = degpol(P);
213 :
214 213 : if (!bnf)
215 : {
216 196 : bnf = Buchall(P, nf_FORCE, maxss(prec, DEFAULTPREC));
217 196 : if (flag) (void)bnfcertify(bnf);
218 : else
219 161 : Ind = floorr(mulru(bnf_get_reg(bnf), 5));
220 : }
221 :
222 213 : nf_get_sign(bnf_get_nf(bnf), &s, &t);
223 213 : fu = bnf_get_fu(bnf);
224 213 : prec_roots = prec + nbits2extraprec(gexpo(Q_primpart(fu)));
225 : for(;;)
226 : {
227 240 : ro = tnf_get_roots(P, prec_roots, s, t);
228 240 : MatFU = Conj_LH(fu, &ALH, ro, prec);
229 240 : if (MatFU) break;
230 27 : prec_roots = precdbl(prec_roots);
231 27 : if (DEBUGLEVEL>1) pari_warn(warnprec, "inithue", prec_roots);
232 : }
233 :
234 213 : dP = ZX_deriv(P);
235 213 : c1 = NULL; /* min |P'(r_i)|, i <= s */
236 619 : for (k=1; k<=s; k++)
237 : {
238 406 : tmp = gabs(poleval(dP,gel(ro,k)),prec);
239 406 : if (!c1 || gcmp(tmp,c1) < 0) c1 = tmp;
240 : }
241 213 : c1 = gdiv(int2n(n-1), c1);
242 213 : c1 = gprec_w(myround(c1, 1), DEFAULTPREC);
243 :
244 213 : c2 = NULL; /* max |r_i - r_j|, i!=j */
245 943 : for (k=1; k<=n; k++)
246 1747 : for (j=k+1; j<=n; j++)
247 : {
248 1017 : tmp = gabs(gsub(gel(ro,j),gel(ro,k)), prec);
249 1017 : if (!c2 || gcmp(c2,tmp) > 0) c2 = tmp;
250 : }
251 213 : c2 = gprec_w(myround(c2, -1), DEFAULTPREC);
252 :
253 213 : if (t==0)
254 86 : 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 213 : if (DEBUGLEVEL>1)
270 0 : err_printf("c1 = %Ps\nc2 = %Ps\nIndice <= %Ps\n", c1, c2, Ind);
271 :
272 213 : ALH = gmul2n(ALH, 1);
273 213 : tnf = cgetg(8,t_VEC); csts = cgetg(9,t_VEC);
274 213 : gel(tnf,1) = P;
275 213 : gel(tnf,2) = bnf;
276 213 : gel(tnf,3) = ro;
277 213 : gel(tnf,4) = ALH;
278 213 : gel(tnf,5) = MatFU;
279 213 : gel(tnf,6) = T_A_Matrices(MatFU, s+t-1, &eps5, prec);
280 213 : gel(tnf,7) = csts;
281 213 : gel(csts,1) = c1; gel(csts,2) = c2; gel(csts,3) = LogHeight(ro, prec);
282 213 : gel(csts,4) = x0; gel(csts,5) = eps5; gel(csts,6) = utoipos(prec);
283 213 : gel(csts,7) = Ind;
284 213 : gel(csts,8) = get_prime_info(bnf);
285 213 : 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 4157 : other_roots(long iroot, long *i1, long *i2)
296 : {
297 4157 : switch (iroot) {
298 1630 : case 1: *i1=2; *i2=3; break;
299 1190 : case 2: *i1=1; *i2=3; break;
300 1337 : default: *i1=1; *i2=2; break;
301 : }
302 4157 : }
303 : /* low precision */
304 4537 : 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 3769 : Baker(baker_s *BS)
309 : {
310 : GEN tmp, B0, hb0, c9, Indc11;
311 : long i1, i2;
312 :
313 3769 : other_roots(BS->iroot, &i1,&i2);
314 : /* Compute a bound for the h_0 */
315 3769 : hb0 = gadd(gmul2n(BS->hal,2), gmul2n(gadd(BS->Hmu,mplog2(DEFAULTPREC)), 1));
316 3769 : tmp = gmul(BS->divro, gdiv(gel(BS->NE,i1), gel(BS->NE,i2)));
317 3769 : tmp = gmax_shallow(gen_1, abslog(tmp));
318 3769 : hb0 = gmax_shallow(hb0, gdiv(tmp, BS->bak));
319 3769 : c9 = gmul(BS->c91,hb0);
320 3769 : c9 = gprec_w(myround(c9, 1), DEFAULTPREC);
321 3769 : Indc11 = rtor(mulir(BS->Ind,BS->c11), DEFAULTPREC);
322 : /* Compute B0 according to Lemma 2.3.3 */
323 3769 : 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 3769 : B0 = gmax_shallow(B0, dbltor(2.71828183));
327 3769 : B0 = gmax_shallow(B0, mulrr(divir(BS->Ind, BS->c10),
328 : mplog(divrr(Indc11, BS->Pi2))));
329 3769 : if (DEBUGLEVEL>1) {
330 0 : err_printf(" B0 = %Ps\n",B0);
331 0 : err_printf(" Baker = %Ps\n",c9);
332 : }
333 3769 : return B0;
334 : }
335 :
336 : /* || x d ||, x t_REAL, d t_INT */
337 : static GEN
338 8208 : errnum(GEN x, GEN d)
339 : {
340 8208 : GEN dx = mulir(d, x), D = subri(dx, roundr(dx));
341 8208 : setabssign(D); return D;
342 : }
343 :
344 : /* Try to reduce the bound through continued fractions; see paper. */
345 : static int
346 4020 : CF_1stPass(GEN *B0, GEN kappa, baker_s *BS)
347 : {
348 4020 : GEN a, b, q, ql, qd, l0, denbound = mulri(*B0, kappa);
349 :
350 4020 : if (cmprr(mulrr(dbltor(0.1),sqrr(denbound)), BS->inverrdelta) > 0)
351 9 : 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 11013 : get_B0Bx(baker_s *BS, GEN l0, GEN *B0, GEN *Bx)
374 : {
375 11013 : GEN t = divrr(mulir(BS->Ind, BS->c15), l0);
376 11013 : *B0 = divrr(mulir(BS->Ind, mplog(t)), BS->c13);
377 11013 : *Bx = sqrtnr(shiftr(t,1), BS->deg);
378 11013 : }
379 :
380 : static int
381 20454 : LLL_1stPass(GEN *pB0, GEN kappa, baker_s *BS, GEN *pBx)
382 : {
383 20454 : GEN B0 = *pB0, Bx = *pBx, lllmat, C, l0, l1, triv;
384 : long e;
385 :
386 20454 : C = grndtoi(mulir(mulii(BS->Ind, kappa),
387 : gpow(B0, dbltor(2.2), DEFAULTPREC)), NULL);
388 :
389 20454 : if (DEBUGLEVEL > 1) err_printf("C (bitsize) : %d\n", expi(C));
390 20454 : lllmat = matid(3);
391 20454 : if (cmpri(B0, BS->Ind) > 0)
392 : {
393 14570 : gcoeff(lllmat, 1, 1) = grndtoi(divri(B0, BS->Ind), NULL);
394 14570 : triv = shiftr(sqrr(B0), 1);
395 : }
396 : else
397 5884 : triv = addir(sqri(BS->Ind), sqrr(B0));
398 :
399 20454 : gcoeff(lllmat, 3, 1) = grndtoi(negr(mulir(C, BS->lambda)), &e);
400 20454 : if (e >= 0) return -1;
401 20446 : gcoeff(lllmat, 3, 2) = grndtoi(negr(mulir(C, BS->delta)), &e);
402 20446 : if (e >= 0) return -1;
403 20446 : gcoeff(lllmat, 3, 3) = C;
404 20446 : lllmat = ZM_lll(lllmat, 0.99, LLL_IM|LLL_INPLACE);
405 :
406 20446 : l0 = gnorml2(gel(lllmat,1));
407 20446 : l0 = subrr(divir(l0, dbltor(1.8262)), triv); /* delta = 0.99 */
408 20446 : if (signe(l0) <= 0) return 0;
409 :
410 13862 : l1 = shiftr(addri(shiftr(B0,1), BS->Ind), -1);
411 13862 : l0 = divri(subrr(sqrtr(l0), l1), C);
412 :
413 13862 : if (signe(l0) <= 0) return 0;
414 :
415 10827 : get_B0Bx(BS, l0, &B0, &Bx);
416 10827 : 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 10827 : *pB0 = B0; *pBx = Bx; return 1;
423 : }
424 :
425 : /* add solution (x,y) if not already known */
426 : static void
427 693 : 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 182 : add_pm(GEN *pS, GEN p, GEN q, GEN z, long d, GEN rhs)
432 : {
433 182 : if (signe(z) == signe(rhs))
434 : {
435 98 : add_sol(pS, p, q);
436 98 : if (!odd(d)) add_sol(pS, negi(p), negi(q));
437 : }
438 : else
439 84 : if (odd(d)) add_sol(pS, negi(p), negi(q));
440 182 : }
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 186 : GuessQi(GEN b, GEN c, GEN *eps)
584 : {
585 186 : const long shift = 65;
586 : GEN Q, Lat;
587 :
588 186 : Lat = matid(3);
589 186 : gcoeff(Lat,3,1) = ground(gmul2n(b, shift));
590 186 : gcoeff(Lat,3,2) = ground(gmul2n(c, shift));
591 186 : gcoeff(Lat,3,3) = int2n(shift);
592 :
593 186 : Q = gel(lllint(Lat),1);
594 186 : if (gequal0(gel(Q,2))) return NULL; /* FAIL */
595 :
596 186 : *eps = mpadd(mpadd(gel(Q,3), mpmul(gel(Q,1),b)), mpmul(gel(Q,2),c));
597 186 : *eps = mpabs_shallow(*eps); return Q;
598 : }
599 :
600 : /* x a t_REAL */
601 : static GEN
602 363 : 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 182 : 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 182 : if (expo(bound) < 0) return bound;
612 168 : d = degpol(P);
613 168 : bndcf = sqrtnr(shiftr(c1,1), d - 2);
614 168 : 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 153 : nmax = 3 + (long)(dbllog2(bound) * 1.44042009041256);
620 153 : bound = myfloor(bound);
621 :
622 486 : for (k = 1; k <= s; k++)
623 : {
624 333 : GEN ro = real_i(gel(roo,k)), t = gboundcf(ro, nmax);
625 : GEN pm1, qm1, p0, q0;
626 :
627 333 : pm1 = gen_0; p0 = gen_1;
628 333 : qm1 = gen_1; q0 = gen_0;
629 :
630 15395 : for (j = 1; j < lg(t); j++)
631 : {
632 : GEN p, q, z, Q, R;
633 : pari_sp av;
634 15388 : p = addii(mulii(p0, gel(t,j)), pm1); pm1 = p0; p0 = p;
635 15388 : q = addii(mulii(q0, gel(t,j)), qm1); qm1 = q0; q0 = q;
636 15388 : if (cmpii(q, bound) > 0) break;
637 15062 : if (DEBUGLEVEL >= 2) err_printf("Checking (+/- %Ps, +/- %Ps)\n",p, q);
638 15062 : av = avma;
639 15062 : z = ZX_Z_eval(ZX_rescale(P,q), p); /* = P(p/q) q^dep(P) */
640 15062 : Q = dvmdii(rhs, z, &R);
641 15062 : if (R != gen_0) { set_avma(av); continue; }
642 252 : setabssign(Q);
643 252 : if (Z_ispowerall(Q, d, &Q))
644 : {
645 112 : if (!is_pm1(Q)) { p = mulii(p, Q); q = mulii(q, Q); }
646 112 : add_pm(pS, p, q, z, d, rhs);
647 : }
648 : }
649 333 : 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 153 : return bndcf;
662 : }
663 :
664 : static void
665 672700 : check_y_root(GEN *pS, GEN P, GEN Y)
666 : {
667 672700 : GEN r = nfrootsQ(P);
668 : long j;
669 673134 : for (j = 1; j < lg(r); j++)
670 434 : if (typ(gel(r,j)) == t_INT) add_sol(pS, gel(r,j), Y);
671 672700 : }
672 :
673 : static void
674 336406 : check_y(GEN *pS, GEN P, GEN poly, GEN Y, GEN rhs)
675 : {
676 336406 : long j, l = lg(poly);
677 336406 : GEN Yn = Y;
678 336406 : gel(P, l-1) = gel(poly, l-1);
679 1345582 : for (j = l-2; j >= 2; j--)
680 : {
681 1009176 : gel(P,j) = mulii(Yn, gel(poly,j));
682 1009176 : if (j > 2) Yn = mulii(Yn, Y);
683 : }
684 336406 : gel(P,2) = subii(gel(P,2), rhs); /* P = poly(Y/y)*y^deg(poly) - rhs */
685 336406 : check_y_root(pS, P, Y);
686 336406 : }
687 :
688 : /* Check for solutions under a small bound (see paper) */
689 : static GEN
690 210 : SmallSols(GEN S, GEN x3, GEN poly, GEN rhs)
691 : {
692 210 : pari_sp av = avma;
693 : GEN X, P, rhs2;
694 210 : long j, l = lg(poly), n = degpol(poly);
695 : ulong y, By;
696 :
697 210 : x3 = myfloor(x3);
698 :
699 210 : if (DEBUGLEVEL>1) err_printf("* Checking for small solutions <= %Ps\n", x3);
700 210 : if (lgefint(x3) > 3)
701 0 : pari_err_OVERFLOW(stack_sprintf("thue (SmallSols): y <= %Ps", x3));
702 210 : By = itou(x3);
703 : /* y = 0 first: solve X^n = rhs */
704 210 : if (odd(n))
705 : {
706 161 : 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 210 : rhs2 = shifti(rhs,1);
715 : /* y != 0 */
716 210 : P = cgetg(l, t_POL); P[1] = poly[1];
717 336504 : for (y = 1; y <= By; y++)
718 : {
719 336294 : pari_sp av2 = avma;
720 336294 : long lS = lg(S);
721 336294 : GEN Y = utoipos(y);
722 : /* try y */
723 336294 : check_y(&S, P, poly, Y, rhs);
724 : /* try -y */
725 1008224 : for (j = l-2; j >= 2; j -= 2) togglesign( gel(P,j) );
726 336294 : if (j == 0) gel(P,2) = subii(gel(P,2), rhs2);
727 336294 : check_y_root(&S, P, utoineg(y));
728 336294 : if (lS == lg(S)) { set_avma(av2); continue; } /* no solution found */
729 :
730 154 : 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 210 : return S;
738 : }
739 :
740 : /* Computes [x]! */
741 : static double
742 196 : fact(double x)
743 : {
744 196 : double ft = 1.0;
745 910 : x = floor(x); while (x>1) { ft *= x; x--; }
746 196 : 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 413 : thueinit(GEN pol, long flag, long prec)
753 : {
754 413 : GEN POL, C, L, fa, tnf, bnf = NULL;
755 413 : pari_sp av = avma;
756 : long s, lfa, dpol;
757 :
758 413 : if (checktnf(pol)) { bnf = checkbnf(gel(pol,2)); pol = gel(pol,1); }
759 413 : if (typ(pol)!=t_POL) pari_err_TYPE("thueinit",pol);
760 413 : dpol = degpol(pol);
761 413 : if (dpol <= 0) pari_err_CONSTPOL("thueinit");
762 406 : RgX_check_ZX(pol, "thueinit");
763 406 : if (varn(pol)) { pol = leafcopy(pol); setvarn(pol, 0); }
764 :
765 406 : POL = Q_primitive_part(pol, &C);
766 406 : L = gen_1;
767 406 : fa = ZX_factor(POL); lfa = lgcols(fa);
768 406 : 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 315 : POL = ZX_primitive_to_monic(POL, &L);
797 315 : C = gdiv(powiu(L, dpol), gel(pol, dpol+2));
798 315 : pol = POL;
799 :
800 315 : s = ZX_sturm_irred(pol);
801 315 : if (s)
802 : {
803 196 : long PREC, n = degpol(pol);
804 196 : double d, dr, dn = (double)n;
805 :
806 196 : if (dpol <= 2) pari_err_DOMAIN("thueinit", "P","=",pol,pol);
807 196 : dr = (double)((s+n-2)>>1); /* s+t-1 */
808 196 : 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 196 : PREC = nbits2prec((long)((5.83 + (dr+4)*5 + log(fact(dr+3)) + (dr+3)*log(dr+2) +
813 196 : (dr+3)*log(d) + log(log(2*d*(dr+2))) + (dr+1))
814 196 : /10.)*32+32);
815 :
816 196 : if (flag == 0) PREC = (long)(2.2 * PREC); /* Lazy, to be improved */
817 196 : if (PREC < prec) PREC = prec;
818 196 : if (DEBUGLEVEL >=2) err_printf("prec = %d\n", PREC);
819 :
820 : for (;;)
821 : {
822 196 : 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 315 : gel(tnf,1) = mkvec3(gel(tnf,1), C, L);
848 315 : 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 359 : v = subrs(u,1); /* ]-1,0] */
860 359 : if (abscmprr(v,u) < 0) u = v;
861 : }
862 : else
863 : {
864 273 : v = addrs(u,1);/* ]0,1] */
865 273 : if (abscmprr(v,u) <= 0) u = v;
866 : }
867 632 : return u;
868 : }
869 : /* i1 != i2 */
870 : static void
871 3769 : 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 3769 : if (BS->r > 1)
876 : {
877 3453 : delta = gel(Delta2,i2);
878 3453 : lambda = gsub(gmul(delta,gel(Lambda,i1)), gel(Lambda,i2));
879 3453 : 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 3769 : BS->delta = delta;
893 3769 : BS->lambda = lambda;
894 3769 : }
895 :
896 : static GEN
897 1066 : get_B0(long i1, GEN Delta2, GEN Lambda, GEN Deps5, long prec, baker_s *BS)
898 : {
899 1066 : GEN B0 = Baker(BS);
900 1066 : long step = 0, i2 = (i1 == 1)? 2: 1;
901 : for(;;) /* i2 from 1 to r unless r = 1 [then i2 = 2] */
902 : {
903 1066 : 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 3369 : GEN oldB0 = B0, kappa = utoipos(10);
908 : long cf;
909 :
910 4027 : for (cf = 0; cf < 10; cf++, kappa = muliu(kappa,10))
911 : {
912 4020 : int res = CF_1stPass(&B0, kappa, BS);
913 4020 : 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(logr_abs(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 2703 : get_Bx_LLL(long i1, GEN Delta2, GEN Lambda, long prec, baker_s *BS)
945 : {
946 2703 : GEN B0 = Baker(BS), Bx = NULL;
947 2703 : long step = 0, i2 = (i1 == 1)? 2: 1;
948 : for(;;) /* i2 from 1 to r unless r = 1 [then i2 = 2] */
949 : {
950 2703 : init_get_B(i1,i2, Delta2,Lambda,NULL, BS, prec);
951 2703 : if (DEBUGLEVEL>1) err_printf(" Entering LLL...\n");
952 : /* Reduce B0 as long as we make progress: newB0 < oldB0 - 0.1 */
953 : for (;;)
954 8525 : {
955 11228 : GEN oldBx = Bx, kappa = utoipos(10);
956 11228 : const long cfMAX = 10;
957 : long cf;
958 :
959 20847 : for (cf = 0; cf < cfMAX; cf++, kappa = muliu(kappa,10))
960 : {
961 20454 : int res = LLL_1stPass(&B0, kappa, BS, &Bx);
962 20454 : if (res < 0) return NULL;
963 20446 : if (res) break;
964 9619 : if (DEBUGLEVEL>1) err_printf("LLL failed. Increasing kappa\n");
965 : }
966 :
967 : /* FIXME: TO BE COMPLETED */
968 11220 : if (!step && cf == cfMAX)
969 : { /* Semirational or totally rational case */
970 : GEN Q, Q1, Q2, ep, q, l0, denbound;
971 :
972 186 : if (! (Q = GuessQi(BS->delta, BS->lambda, &ep)) ) break;
973 :
974 : /* Q[2] != 0 */
975 186 : Q1 = absi_shallow(gel(Q,1));
976 186 : Q2 = absi_shallow(gel(Q,2));
977 186 : denbound = gadd(mulri(B0, Q1), mulii(BS->Ind, Q2));
978 186 : q = denom_i( bestappr(BS->delta, denbound) );
979 186 : l0 = divri(subrr(errnum(BS->delta, q), ep), Q2);
980 186 : if (signe(l0) <= 0) break;
981 :
982 186 : get_B0Bx(BS, l0, &B0, &Bx);
983 186 : if (DEBUGLEVEL>1)
984 0 : err_printf("Semirat. reduction: B0 -> %Ps x <= %Ps\n",B0, Bx);
985 : }
986 : /* if no progress, stop */
987 11220 : 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 182 : LargeSols(GEN P, GEN tnf, GEN rhs, GEN ne)
998 : {
999 182 : 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 182 : pari_sp av = avma;
1004 :
1005 182 : prec = 0; /*-Wall*/
1006 182 : bnf = NULL; /*-Wall*/
1007 182 : iroot = 1;
1008 182 : ine = 1;
1009 :
1010 199 : START:
1011 199 : if (S) /* restart from precision problems */
1012 : {
1013 17 : S = gerepilecopy(av, S);
1014 17 : prec = precdbl(prec);
1015 17 : if (DEBUGLEVEL) pari_warn(warnprec,"thue",prec);
1016 17 : tnf = inithue(P, bnf, 0, prec);
1017 : }
1018 : else
1019 182 : S = cgetg(1, t_VEC);
1020 199 : bnf= gel(tnf,2);
1021 199 : nf = bnf_get_nf(bnf);
1022 199 : csts = gel(tnf,7);
1023 199 : nf_get_sign(nf, &s, &t);
1024 199 : BS.r = r = s+t-1; n = degpol(P);
1025 199 : ro = gel(tnf,3);
1026 199 : ALH = gel(tnf,4);
1027 199 : MatFU = gel(tnf,5);
1028 199 : A = gel(tnf,6);
1029 199 : c1 = mpmul(absi_shallow(rhs), gel(csts,1));
1030 199 : c2 = gel(csts,2);
1031 199 : BS.hal = gel(csts,3);
1032 199 : x0 = gel(csts,4);
1033 199 : eps5 = gel(csts,5);
1034 199 : Prec = itos(gel(csts,6));
1035 199 : BS.Ind = gel(csts,7);
1036 199 : BS.MatFU = MatFU;
1037 199 : BS.bak = muluu(n, (n-1)*(n-2)); /* safe */
1038 199 : BS.deg = n;
1039 199 : prinfo = gel(csts,8);
1040 :
1041 199 : if (t) x0 = gmul(x0, absisqrtn(rhs, n, Prec));
1042 199 : tmp = divrr(c1,c2);
1043 199 : c3 = mulrr(dbltor(1.39), tmp);
1044 199 : c4 = mulur(n-1, c3);
1045 199 : c14 = mulrr(c4, vecmax_shallow(RgM_sumcol(gabs(A,DEFAULTPREC))));
1046 :
1047 199 : x1 = gmax_shallow(x0, sqrtnr(shiftr(tmp,1),n));
1048 199 : x2 = gmax_shallow(x1, sqrtnr(mulur(10,c14), n));
1049 199 : x3 = gmax_shallow(x2, sqrtnr(shiftr(c14, EXPO1+1),n));
1050 398 : c90 = gmul(shiftr(mulur(18,mppi(DEFAULTPREC)), 5*(4+r)),
1051 199 : gmul(gmul(mpfact(r+3), powiu(muliu(BS.bak,r+2), r+3)),
1052 199 : glog(muliu(BS.bak,2*(r+2)),DEFAULTPREC)));
1053 :
1054 199 : dP = ZX_deriv(P);
1055 199 : Delta0 = RgM_sumcol(A);
1056 :
1057 570 : for (; iroot<=s; iroot++)
1058 : {
1059 388 : GEN Delta = Delta0, Delta2, D, Deps5, MatNE, Hmu, diffRo, c5, c7, ro0;
1060 : long i1, iroot1, iroot2, k;
1061 :
1062 388 : if (iroot <= r) Delta = RgC_add(Delta, RgC_Rg_mul(gel(A,iroot), stoi(-n)));
1063 388 : D = gabs(Delta,Prec); i1 = vecindexmax(D);
1064 388 : c5 = gel(D, i1);
1065 388 : Delta2 = RgC_Rg_div(Delta, gel(Delta, i1));
1066 388 : c5 = myround(gprec_w(c5,DEFAULTPREC), 1);
1067 388 : Deps5 = divrr(subrr(c5,eps5), eps5);
1068 388 : c7 = mulur(r,c5);
1069 388 : BS.c10 = divur(n,c7);
1070 388 : BS.c13 = divur(n,c5);
1071 388 : 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 388 : prec = Prec;
1078 : for (;;)
1079 : {
1080 388 : 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 388 : ro0 = gel(ro,iroot);
1086 388 : BS.ro = ro;
1087 388 : BS.iroot = iroot;
1088 388 : BS.Pi = mppi(prec);
1089 388 : BS.Pi2 = Pi2n(1,prec);
1090 388 : diffRo = cgetg(r+1, t_VEC);
1091 1156 : for (k=1; k<=r; k++)
1092 : {
1093 768 : GEN z = gel(ro,k);
1094 768 : z = (k == iroot)? gdiv(rhs, poleval(dP, z)): gsub(ro0, z);
1095 768 : gel(diffRo,k) = gabs(z, prec);
1096 : }
1097 388 : other_roots(iroot, &iroot1,&iroot2);
1098 388 : BS.divro = gdiv(gsub(ro0, gel(ro,iroot2)), gsub(ro0, gel(ro,iroot1)));
1099 : /* Compute h_1....h_r */
1100 388 : c91 = c90;
1101 1156 : for (k=1; k<=r; k++)
1102 : {
1103 768 : GEN z = gdiv(gcoeff(MatFU,iroot1,k), gcoeff(MatFU,iroot2,k));
1104 768 : z = gmax_shallow(gen_1, abslog(z));
1105 768 : c91 = gmul(c91, gmax_shallow(gel(ALH,k), gdiv(z, BS.bak)));
1106 : }
1107 388 : BS.c91 = c91;
1108 :
1109 4140 : for (; ine<lg(ne); ine++)
1110 : {
1111 3769 : pari_sp av2 = avma;
1112 3769 : long lS = lg(S);
1113 : GEN Lambda, B0, c6, c8;
1114 3769 : GEN NE = gel(MatNE,ine), v = cgetg(r+1,t_COL);
1115 :
1116 3769 : if (DEBUGLEVEL>1) err_printf(" - norm sol. no %ld/%ld\n",ine,lg(ne)-1);
1117 11649 : for (k=1; k<=r; k++)
1118 : {
1119 7880 : GEN z = gdiv(gel(diffRo,k), gabs(gel(NE,k), prec));
1120 7880 : gel(v,k) = glog(z, prec);
1121 : }
1122 3769 : Lambda = RgM_RgC_mul(A,v);
1123 :
1124 3769 : c6 = addrr(dbltor(0.1), vecmax_shallow(gabs(Lambda,DEFAULTPREC)));
1125 3769 : c6 = myround(c6, 1);
1126 3769 : c8 = addrr(dbltor(1.23), mulur(r,c6));
1127 3769 : BS.c11= mulrr(shiftr(c3,1) , mpexp(divrr(mulur(n,c8),c7)));
1128 3769 : BS.c15= mulrr(shiftr(c14,1), mpexp(divrr(mulur(n,c6),c5)));
1129 3769 : BS.NE = NE;
1130 3769 : BS.Hmu = gel(Hmu,ine);
1131 3769 : if (is_pm1(BS.Ind))
1132 : {
1133 1066 : GEN mu = gel(ne,ine), Lmu = cgetg(lg(prinfo),t_VEC);
1134 : long i, j;
1135 :
1136 2139 : for (i = 1; i < lg(prinfo); i++)
1137 : {
1138 1073 : GEN v = gel(prinfo,i), PR = gel(v,3), L = cgetg(4, t_VECSMALL);
1139 4292 : for (j = 1; j <= 3; j++) L[j] = itou(nf_to_Fq(nf, mu, gel(PR,j)));
1140 1073 : gel(Lmu, i) = L;
1141 : }
1142 2123 : 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 17 : goto START;
1146 1057 : if (lg(S) == lS) set_avma(av2);
1147 : }
1148 : else
1149 : {
1150 2703 : if (! (Bx = get_Bx_LLL(i1, Delta2, Lambda, prec, &BS)) )
1151 8 : goto START;
1152 2695 : x3 = gerepileupto(av2, gmax_shallow(Bx, x3));
1153 : }
1154 : }
1155 371 : ine = 1;
1156 : }
1157 182 : x3 = gmax_shallow(x0, MiddleSols(&S, x3, ro, P, rhs, s, c1));
1158 182 : 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 343 : filter_sol_x(GEN S, GEN L)
1164 : {
1165 : long i, k, l;
1166 343 : 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 175 : tnf_get_Ind(GEN tnf) { return gmael(tnf,7,7); }
1191 : static GEN
1192 294 : 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 224 : get_ne(GEN bnf, GEN a, GEN fa, GEN Ind)
1203 : {
1204 224 : if (DEBUGLEVEL) maybe_warn(bnf,a,Ind);
1205 224 : 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 392 : thue(GEN tnf, GEN rhs, GEN ne)
1262 : {
1263 392 : pari_sp av = avma;
1264 : GEN POL, C, L, x3, S;
1265 :
1266 392 : if (typ(tnf) == t_POL) tnf = thueinit(tnf, 0, DEFAULTPREC);
1267 392 : if (!checktnf(tnf)) pari_err_TYPE("thue [please apply thueinit()]", tnf);
1268 392 : if (typ(rhs) != t_INT) pari_err_TYPE("thue",rhs);
1269 392 : 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 392 : POL = gel(tnf,1);
1273 392 : C = gel(POL,2); rhs = gmul(C, rhs);
1274 392 : if (typ(rhs) != t_INT) { set_avma(av); return cgetg(1, t_VEC); }
1275 392 : 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 364 : L = gel(POL,3);
1285 364 : POL = gel(POL,1);
1286 364 : if (lg(tnf) == 8)
1287 : {
1288 182 : if (!ne)
1289 : {
1290 154 : GEN F = absZ_factor(rhs);
1291 154 : ne = get_ne(tnf_get_bnf(tnf), rhs, F, tnf_get_Ind(tnf));
1292 : }
1293 182 : if (lg(ne) == 1) { set_avma(av); return cgetg(1, t_VEC); }
1294 182 : 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 343 : S = filter_sol_x(S, L);
1384 343 : S = gen_sort_uniq(S, (void*)lexcmp, cmp_nodata);
1385 343 : 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 29918 : test_sol(struct sol_abs *T, long i)
1414 : {
1415 : long k, l;
1416 : GEN s;
1417 :
1418 29918 : if (T->partrel && !ZV_dvd(gel(T->partrel, i), T->cyc)) return;
1419 13825 : 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 13825 : gel(T->normsol, ++T->sindex) = s = cgetg_copy(T->u, &l);
1428 57708 : for (k=1; k <= i; k++) s[k] = T->u[k];
1429 25081 : for ( ; k < l; k++) s[k] = 0;
1430 13825 : 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 78806 : isintnorm_loop(struct sol_abs *T, long i)
1458 : {
1459 78806 : 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 45255 : 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 28084 : if (T->S[i] % T->f[i+1]) return;
1473 14049 : q = T->S[i] / T->f[i+1];
1474 14049 : i++; T->u[i] = q;
1475 14049 : if (T->partrel) fix_partrel(T,i);
1476 14049 : 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 21770 : get_sol_abs(struct sol_abs *T, GEN bnf, GEN nf, GEN fact, GEN *ptPR)
1507 : {
1508 21770 : GEN P = gel(fact,1), E = gel(fact,2), PR;
1509 21770 : long N = nf_get_degree(nf), nP = lg(P)-1, Ngen, max, nPR, i, j;
1510 :
1511 21770 : max = nP*N; /* upper bound for T->nPR */
1512 21770 : T->f = new_chunk(max+1);
1513 21770 : T->n = new_chunk(max+1);
1514 21770 : T->next = new_chunk(max+1);
1515 21770 : *ptPR = PR = cgetg(max+1, t_COL); /* length to be fixed later */
1516 :
1517 21770 : nPR = 0;
1518 55482 : for (i = 1; i <= nP; i++)
1519 : {
1520 41797 : GEN L = idealprimedec(nf, gel(P,i));
1521 41797 : long lL = lg(L), gcd, k, v;
1522 41797 : ulong vn = itou(gel(E,i));
1523 :
1524 : /* check that gcd_{P | p} f_P divides n_p */
1525 41797 : gcd = pr_get_f(gel(L,1));
1526 41818 : for (j=2; gcd > 1 && j < lL; j++) gcd = ugcd(gcd, pr_get_f(gel(L,j)));
1527 41797 : if (gcd > 1 && vn % gcd)
1528 : {
1529 8085 : if (DEBUGLEVEL>2) err_printf("gcd f_P does not divide n_p\n");
1530 8085 : return 0;
1531 : }
1532 33712 : v = (i==nP)? 0: nPR + lL;
1533 92393 : for (k = 1; k < lL; k++)
1534 : {
1535 58681 : GEN pr = gel(L,k);
1536 58681 : gel(PR, ++nPR) = pr;
1537 58681 : T->f[nPR] = pr_get_f(pr) / gcd;
1538 58681 : T->n[nPR] = vn / gcd;
1539 58681 : T->next[nPR] = v;
1540 : }
1541 : }
1542 13685 : T->nPR = nPR;
1543 13685 : setlg(PR, nPR + 1);
1544 :
1545 13685 : T->u = cgetg(nPR+1, t_VECSMALL);
1546 13685 : T->S = new_chunk(nPR+1);
1547 13685 : if (bnf) { T->cyc = bnf_get_cyc(bnf); Ngen = lg(T->cyc)-1; }
1548 105 : else { T->cyc = NULL; Ngen = 0; }
1549 13685 : if (Ngen == 0)
1550 749 : 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 13685 : 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 13685 : T->smax = 511;
1579 13685 : T->normsol = new_chunk(T->smax+1);
1580 13685 : T->S[0] = T->n[1];
1581 13685 : T->next[0] = 1;
1582 13685 : T->sindex = 0;
1583 13685 : 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 3108 : get_unit_1(GEN bnf, GEN *unit)
1589 : {
1590 3108 : GEN v, nf = bnf_get_nf(bnf);
1591 3108 : long i, n = nf_get_degree(nf);
1592 :
1593 3108 : if (DEBUGLEVEL > 2) err_printf("looking for a fundamental unit of norm -1\n");
1594 3108 : 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 21742 : bnfisintnormabs(GEN bnf, GEN a)
1603 : {
1604 : struct sol_abs T;
1605 : GEN nf, res, PR, F;
1606 : long i;
1607 :
1608 21742 : if ((F = check_arith_all(a,"bnfisintnormabs")))
1609 : {
1610 700 : a = typ(a) == t_VEC? gel(a,1): factorback(F);
1611 700 : if (signe(a) < 0) F = clean_Z_factor(F);
1612 : }
1613 21742 : nf = bnf_get_nf(bnf);
1614 21742 : if (!signe(a)) return mkvec(gen_0);
1615 21728 : if (is_pm1(a)) return mkvec(gen_1);
1616 21658 : if (!F) F = absZ_factor(a);
1617 21658 : if (!get_sol_abs(&T, bnf, nf, F, &PR)) return cgetg(1, t_VEC);
1618 : /* |a| > 1 => T.nPR > 0 */
1619 13580 : res = cgetg(T.sindex+1, t_VEC);
1620 27083 : for (i=1; i<=T.sindex; i++)
1621 : {
1622 13503 : GEN x = vecsmall_to_col( gel(T.normsol,i) );
1623 13502 : x = isprincipalfact(bnf, NULL, PR, x, nf_FORCE | nf_GEN_IF_PRINCIPAL);
1624 13503 : gel(res,i) = nf_to_scalar_or_alg(nf, x); /* x solution, up to sign */
1625 : }
1626 13580 : return res;
1627 : }
1628 :
1629 : /* true nf */
1630 : GEN
1631 133 : ideals_by_norm(GEN nf, GEN a)
1632 : {
1633 : struct sol_abs T;
1634 : GEN res, PR, F;
1635 : long i;
1636 :
1637 133 : if ((F = check_arith_pos(a,"ideals_by_norm")))
1638 0 : a = typ(a) == t_VEC? gel(a,1): factorback(F);
1639 133 : if (is_pm1(a)) return mkvec(trivial_fact());
1640 112 : if (!F) F = absZ_factor(a);
1641 112 : 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 21749 : bnfisintnorm_i(GEN bnf, GEN a, long sa, GEN z)
1655 : {
1656 21749 : GEN nf = bnf_get_nf(bnf), T = nf_get_pol(nf), f = nf_get_index(nf), unit=NULL;
1657 21749 : GEN Tp, A = signe(a) == sa? a: negi(a);
1658 21749 : long sNx, i, j, N = degpol(T), l = lg(z);
1659 21749 : long norm_1 = 0; /* gcc -Wall */
1660 21749 : ulong p, Ap = 0; /* gcc -Wall */
1661 : forprime_t S;
1662 21749 : if (!signe(a)) return z;
1663 21735 : u_forprime_init(&S,3,ULONG_MAX);
1664 30338 : while((p = u_forprime_next(&S)))
1665 30338 : if (umodiu(f,p)) { Ap = umodiu(A,p); if (Ap) break; }
1666 21735 : 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 35315 : for (i = j = 1; i < l; i++)
1672 : {
1673 13580 : GEN x = gel(z,i);
1674 13580 : int xpol = (typ(x) == t_POL);
1675 :
1676 13580 : if (xpol)
1677 : {
1678 12943 : GEN dx, y = Q_remove_denom(x,&dx);
1679 12943 : ulong Np = Flx_resultant(Tp, ZX_to_Flx(y,p), p);
1680 12943 : 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 12943 : sNx = dA == Np? sa: -sa;
1683 : }
1684 : else
1685 637 : sNx = gsigne(x) < 0 && odd(N) ? -1 : 1;
1686 13580 : if (sNx != sa)
1687 : {
1688 5054 : if (! unit) norm_1 = get_unit_1(bnf, &unit);
1689 5054 : if (!norm_1)
1690 : {
1691 77 : if (DEBUGLEVEL > 2) err_printf("%Ps eliminated because of sign\n",x);
1692 77 : continue;
1693 : }
1694 4977 : 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 13503 : gel(z,j++) = x;
1698 : }
1699 21735 : setlg(z, j); return z;
1700 : }
1701 : GEN
1702 21490 : bnfisintnorm(GEN bnf, GEN a)
1703 : {
1704 21490 : pari_sp av = avma;
1705 : GEN ne;
1706 21490 : bnf = checkbnf(bnf);
1707 21490 : 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 : }
|
