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 : static GEN
304 4540 : 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 3770 : Baker(baker_s *BS)
309 : {
310 : GEN tmp, B0, hb0, c9, Indc11;
311 : long i1, i2;
312 :
313 3770 : other_roots(BS->iroot, &i1,&i2);
314 : /* Compute a bound for the h_0 */
315 3770 : hb0 = gadd(gmul2n(BS->hal,2), gmul2n(gadd(BS->Hmu,mplog2(DEFAULTPREC)), 1));
316 3770 : tmp = gmul(BS->divro, gdiv(gel(BS->NE,i1), gel(BS->NE,i2)));
317 3770 : tmp = gmax_shallow(gen_1, abslog(tmp));
318 3770 : hb0 = gmax_shallow(hb0, gdiv(tmp, BS->bak));
319 3770 : c9 = gmul(BS->c91,hb0);
320 3770 : c9 = gprec_w(myround(c9, 1), DEFAULTPREC);
321 3770 : Indc11 = rtor(mulir(BS->Ind,BS->c11), DEFAULTPREC);
322 : /* Compute B0 according to Lemma 2.3.3 */
323 3770 : 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 3770 : B0 = gmax_shallow(B0, dbltor(2.71828183));
327 3770 : B0 = gmax_shallow(B0, mulrr(divir(BS->Ind, BS->c10),
328 : mplog(divrr(Indc11, BS->Pi2))));
329 3770 : if (DEBUGLEVEL>1) {
330 0 : err_printf(" B0 = %Ps\n",B0);
331 0 : err_printf(" Baker = %Ps\n",c9);
332 : }
333 3770 : return B0;
334 : }
335 :
336 : /* || x d ||, x t_REAL, d t_INT */
337 : static GEN
338 8519 : errnum(GEN x, GEN d)
339 : {
340 8519 : GEN dx = mulir(d, x), D = subri(dx, roundr(dx));
341 8519 : setabssign(D); return D;
342 : }
343 :
344 : /* Try to reduce the bound through continued fractions. *B0 t_REAL. */
345 : static int
346 4177 : CF_1stPass(GEN *B0, GEN kappa, baker_s *BS)
347 : {
348 4177 : GEN a, b, q, ql, qd, l0, denbound = mulri(*B0, kappa);
349 :
350 4177 : if (cmprr(mulrr(dbltor(0.1),sqrr(denbound)), BS->inverrdelta) > 0)
351 9 : return -1;
352 :
353 4168 : q = denom_i( bestappr(BS->delta, denbound) );
354 4168 : qd = errnum(BS->delta, q);
355 4168 : ql = errnum(BS->lambda,q);
356 :
357 4168 : l0 = subrr(ql, addrr(mulrr(qd, *B0), divri(dbltor(0.1),kappa)));
358 4168 : if (signe(l0) <= 0) return 0;
359 :
360 3488 : if (BS->r > 1) {
361 2732 : 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 3488 : *B0 = divrr(mplog(divrr(mulir(q,a), l0)), b);
368 3488 : if (DEBUGLEVEL>1) err_printf(" B0 -> %Ps\n",*B0);
369 3488 : return 1;
370 : }
371 :
372 : static void
373 10832 : get_B0Bx(baker_s *BS, GEN l0, GEN *B0, GEN *Bx)
374 : {
375 10832 : GEN t = divrr(mulir(BS->Ind, BS->c15), l0);
376 10832 : *B0 = divrr(mulir(BS->Ind, mplog(t)), BS->c13);
377 10832 : *Bx = sqrtnr(shiftr(t,1), BS->deg);
378 10832 : }
379 :
380 : static int
381 20178 : LLL_1stPass(GEN *pB0, GEN kappa, baker_s *BS, GEN *pBx)
382 : {
383 20178 : GEN B0 = *pB0, Bx = *pBx, lllmat, C, l0, l1, triv;
384 : long e;
385 :
386 20178 : C = grndtoi(mulir(mulii(BS->Ind, kappa),
387 : gpow(B0, dbltor(2.2), DEFAULTPREC)), NULL);
388 :
389 20178 : if (DEBUGLEVEL > 1) err_printf("C (bitsize) : %d\n", expi(C));
390 20178 : lllmat = matid(3);
391 20178 : if (cmpri(B0, BS->Ind) > 0)
392 : {
393 14324 : gcoeff(lllmat, 1, 1) = grndtoi(divri(B0, BS->Ind), NULL);
394 14324 : triv = shiftr(sqrr(B0), 1);
395 : }
396 : else
397 5854 : triv = addir(sqri(BS->Ind), sqrr(B0));
398 :
399 20178 : gcoeff(lllmat, 3, 1) = grndtoi(negr(mulir(C, BS->lambda)), &e);
400 20178 : if (e >= 0) return -1;
401 20169 : gcoeff(lllmat, 3, 2) = grndtoi(negr(mulir(C, BS->delta)), &e);
402 20169 : if (e >= 0) return -1;
403 20169 : gcoeff(lllmat, 3, 3) = C;
404 20169 : lllmat = ZM_lll(lllmat, 0.99, LLL_IM|LLL_INPLACE);
405 :
406 20169 : l0 = gnorml2(gel(lllmat,1));
407 20169 : l0 = subrr(divir(l0, dbltor(1.8262)), triv); /* delta = 0.99 */
408 20169 : if (signe(l0) <= 0) return 0;
409 :
410 13667 : l1 = shiftr(addri(shiftr(B0,1), BS->Ind), -1);
411 13667 : l0 = divri(subrr(sqrtr(l0), l1), C);
412 :
413 13667 : if (signe(l0) <= 0) return 0;
414 :
415 10649 : get_B0Bx(BS, l0, &B0, &Bx);
416 10649 : 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 10649 : *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 29476 : round_to_b(GEN v, long B, long b, GEN Delta2, long i1, GEN L)
467 : {
468 29476 : long i, l = lg(v);
469 29476 : if (!b)
470 : {
471 2212 : for (i = 1; i < l; i++)
472 : {
473 : long c;
474 1929 : if (i == i1)
475 1065 : c = 0;
476 : else
477 : {
478 864 : GEN d = gneg(gel(L,i));
479 : long e;
480 864 : d = grndtoi(d,&e);
481 864 : if (e > -EXPO1 || is_bigint(d)) return 0;
482 45 : c = itos(d); if (labs(c) > B) return 0;
483 : }
484 1110 : v[i] = c;
485 : }
486 : }
487 : else
488 : {
489 56727 : for (i = 1; i < l; i++)
490 : {
491 : long c;
492 36197 : if (i == i1)
493 28207 : c = b;
494 : else
495 : {
496 7990 : GEN d = gsub(gmulgs(gel(Delta2,i), b), gel(L,i));
497 : long e;
498 7990 : d = grndtoi(d,&e);
499 7990 : if (e > -EXPO1 || is_bigint(d)) return 0;
500 146 : c = itos(d); if (labs(c) > B) return 0;
501 : }
502 28353 : v[i] = c;
503 : }
504 : }
505 20813 : return 1;
506 : }
507 :
508 : /* mu \prod U[i]^b[i] */
509 : static ulong
510 62439 : Fl_factorback(ulong mu, GEN U, GEN b, ulong p)
511 : {
512 62439 : long i, l = lg(U);
513 62439 : ulong r = mu;
514 125451 : for (i = 1; i < l; i++)
515 : {
516 63012 : long c = b[i];
517 63012 : ulong u = U[i];
518 63012 : if (!c) continue;
519 61761 : if (c < 0) { u = Fl_inv(u,p); c = -c; }
520 61761 : r = Fl_mul(r, Fl_powu(u,c,p), p);
521 : }
522 62439 : 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 20813 : check_pr(GEN bi, GEN Lmu, GEN L)
529 : {
530 20813 : GEN A = gel(L,1), U = gel(L,2);
531 20813 : ulong a = A[1], b = A[2], c = A[3], p = A[4];
532 20813 : ulong r = Fl_mul(Fl_sub(c,b,p), Fl_factorback(Lmu[1],gel(U,1),bi, p), p);
533 20813 : ulong s = Fl_mul(Fl_sub(a,c,p), Fl_factorback(Lmu[2],gel(U,2),bi, p), p);
534 20813 : ulong t = Fl_mul(Fl_sub(b,a,p), Fl_factorback(Lmu[3],gel(U,3),bi, p), p);
535 20813 : return Fl_add(Fl_add(r,s,p),t,p) == 0;
536 : }
537 : static int
538 20813 : check_prinfo(GEN b, GEN Lmu, GEN prinfo)
539 : {
540 : long i;
541 20890 : for (i = 1; i < lg(prinfo); i++)
542 20813 : 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 1102 : 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 1102 : long bi1, i, B = itos(gceil(B0)), l = lg(Delta2);
552 1102 : GEN b = cgetg(l,t_VECSMALL), L = cgetg(l,t_VEC);
553 :
554 3075 : for (i = 1; i < l; i++)
555 : {
556 1973 : if (i == i1)
557 1102 : gel(L,i) = gen_0;
558 : else
559 : {
560 871 : GEN delta = gel(Delta2,i);
561 871 : gel(L, i) = gsub(gmul(delta,gel(Lambda,i1)), gel(Lambda,i));
562 : }
563 : }
564 30578 : for (bi1 = -B; bi1 <= B; bi1++)
565 : {
566 : GEN z1, z2;
567 29476 : if (!round_to_b(b, B, bi1, Delta2, i1, L)) continue;
568 20813 : 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 1102 : return 1;
579 : }
580 :
581 : /* b and c t_REAL; find q1,q2,q3 st q1 b + q2 c + q3 ~ 0 */
582 : static GEN
583 183 : GuessQi(GEN b, GEN c, GEN *eps)
584 : {
585 183 : const long shift = 65;
586 : GEN Q, L;
587 :
588 183 : L = matid(3);
589 183 : gcoeff(L,3,1) = roundr(shiftr(b, shift));
590 183 : gcoeff(L,3,2) = roundr(shiftr(c, shift));
591 183 : gcoeff(L,3,3) = int2n(shift);
592 :
593 183 : Q = gel(lllint(L),1);
594 183 : if (!signe(gel(Q,2))) return NULL; /* FAIL */
595 :
596 183 : *eps = addrr(addir(gel(Q,3), mulir(gel(Q,1),b)), mulir(gel(Q,2),c));
597 183 : setabssign(*eps); return Q;
598 : }
599 :
600 : /* x a t_REAL */
601 : static GEN
602 362 : 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 152 : nmax = 3 + (long)(dbllog2(bound) * 1.44042009041256);
620 152 : bound = myfloor(bound);
621 :
622 481 : for (k = 1; k <= s; k++)
623 : {
624 329 : GEN ro = real_i(gel(roo,k)), t = gboundcf(ro, nmax);
625 : GEN pm1, qm1, p0, q0;
626 :
627 329 : pm1 = gen_0; p0 = gen_1;
628 329 : qm1 = gen_1; q0 = gen_0;
629 :
630 13947 : for (j = 1; j < lg(t); j++)
631 : {
632 : GEN p, q, z, Q, R;
633 : pari_sp av;
634 13941 : p = addii(mulii(p0, gel(t,j)), pm1); pm1 = p0; p0 = p;
635 13941 : q = addii(mulii(q0, gel(t,j)), qm1); qm1 = q0; q0 = q;
636 13941 : if (cmpii(q, bound) > 0) break;
637 13618 : if (DEBUGLEVEL >= 2) err_printf("Checking (+/- %Ps, +/- %Ps)\n",p, q);
638 13618 : av = avma;
639 13618 : z = ZX_Z_eval(ZX_rescale(P,q), p); /* = P(p/q) q^dep(P) */
640 13618 : Q = dvmdii(rhs, z, &R);
641 13618 : 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 329 : if (j == lg(t))
650 : {
651 : long prec;
652 6 : if (j > nmax) pari_err_BUG("thue [short continued fraction]");
653 : /* the theoretical value is bit_prec = gexpo(ro)+1+log2(bound) */
654 6 : prec = precdbl(precision(ro));
655 6 : if (DEBUGLEVEL>1) pari_warn(warnprec,"thue",prec);
656 6 : roo = ZX_realroots_irred(P, prec);
657 6 : if (lg(roo)-1 != s) pari_err_BUG("thue [realroots]");
658 6 : k--;
659 : }
660 : }
661 152 : 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 : (void)gc_all(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 gc_GEN(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 gc_GEN(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 gc_GEN(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 3770 : 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 3770 : if (BS->r > 1)
876 : {
877 3454 : delta = gel(Delta2,i2);
878 3454 : lambda = gsub(gmul(delta,gel(Lambda,i1)), gel(Lambda,i2));
879 3454 : 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), 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 3770 : BS->delta = delta;
893 3770 : BS->lambda = lambda;
894 3770 : }
895 :
896 : static GEN
897 1111 : get_B0(long i1, GEN Delta2, GEN Lambda, GEN Deps5, long prec, baker_s *BS)
898 : {
899 1111 : GEN B0 = Baker(BS); /* t_REAL */
900 1111 : long step = 0, i2 = (i1 == 1)? 2: 1;
901 : for(;;) /* i2 from 1 to r unless r = 1 [then i2 = 2] */
902 : {
903 1111 : 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 2393 : {
907 3504 : GEN oldB0 = B0, kappa = utoipos(10);
908 : long cf;
909 :
910 4184 : for (cf = 0; cf < 10; cf++, kappa = muliu(kappa,10))
911 : {
912 4177 : int res = CF_1stPass(&B0, kappa, BS);
913 4177 : if (res < 0) return NULL; /* prec problem */
914 4168 : if (res) break;
915 680 : if (DEBUGLEVEL>1) err_printf("CF failed. Increasing kappa\n");
916 : }
917 3495 : 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 = addii(floorr(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 3495 : if (gcmp(oldB0, addrr(B0, dbltor(0.1))) <= 0)
933 1102 : return gmin_shallow(oldB0, B0);
934 2393 : 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 2659 : get_Bx_LLL(long i1, GEN Delta2, GEN Lambda, long prec, baker_s *BS)
945 : {
946 2659 : GEN B0 = Baker(BS), Bx = NULL;
947 2659 : long step = 0, i2 = (i1 == 1)? 2: 1;
948 : for(;;) /* i2 from 1 to r unless r = 1 [then i2 = 2] */
949 : {
950 2659 : init_get_B(i1,i2, Delta2,Lambda,NULL, BS, prec);
951 2659 : if (DEBUGLEVEL>1) err_printf(" Entering LLL...\n");
952 : /* Reduce B0 as long as we make progress: newB0 < oldB0 - 0.1 */
953 : for (;;)
954 8386 : {
955 11045 : GEN oldBx = Bx, kappa = utoipos(10);
956 11045 : const long cfMAX = 10;
957 : long cf;
958 :
959 20565 : for (cf = 0; cf < cfMAX; cf++, kappa = muliu(kappa,10))
960 : {
961 20178 : int res = LLL_1stPass(&B0, kappa, BS, &Bx);
962 20178 : if (res < 0) return NULL;
963 20169 : if (res) break;
964 9520 : if (DEBUGLEVEL>1) err_printf("LLL failed. Increasing kappa\n");
965 : }
966 :
967 : /* FIXME: TO BE COMPLETED */
968 11036 : if (!step && cf == cfMAX)
969 : { /* Semirational or totally rational case */
970 : GEN Q, Q1, Q2, ep, q, l0, denbound;
971 :
972 183 : if (! (Q = GuessQi(BS->delta, BS->lambda, &ep)) ) break;
973 :
974 : /* Q[2] != 0 */
975 183 : Q1 = absi_shallow(gel(Q,1));
976 183 : Q2 = absi_shallow(gel(Q,2));
977 183 : denbound = gadd(mulri(B0, Q1), mulii(BS->Ind, Q2));
978 183 : q = denom_i( bestappr(BS->delta, denbound) );
979 183 : l0 = divri(subrr(errnum(BS->delta, q), ep), Q2);
980 183 : if (signe(l0) <= 0) break;
981 :
982 183 : get_B0Bx(BS, l0, &B0, &Bx);
983 183 : if (DEBUGLEVEL>1)
984 0 : err_printf("Semirat. reduction: B0 -> %Ps x <= %Ps\n",B0, Bx);
985 : }
986 : /* if no progress, stop */
987 11036 : 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 200 : START:
1011 200 : if (S) /* restart from precision problems */
1012 : {
1013 18 : S = gc_GEN(av, S);
1014 18 : prec = precdbl(prec);
1015 18 : if (DEBUGLEVEL) pari_warn(warnprec,"thue",prec);
1016 18 : tnf = inithue(P, bnf, 0, prec);
1017 : }
1018 : else
1019 182 : S = cgetg(1, t_VEC);
1020 200 : bnf= gel(tnf,2);
1021 200 : nf = bnf_get_nf(bnf);
1022 200 : csts = gel(tnf,7);
1023 200 : nf_get_sign(nf, &s, &t);
1024 200 : BS.r = r = s+t-1; n = degpol(P);
1025 200 : ro = gel(tnf,3);
1026 200 : ALH = gel(tnf,4);
1027 200 : MatFU = gel(tnf,5);
1028 200 : A = gel(tnf,6);
1029 200 : c1 = mpmul(absi_shallow(rhs), gel(csts,1));
1030 200 : c2 = gel(csts,2);
1031 200 : BS.hal = gel(csts,3);
1032 200 : x0 = gel(csts,4);
1033 200 : eps5 = gel(csts,5);
1034 200 : Prec = itos(gel(csts,6));
1035 200 : BS.Ind = gel(csts,7);
1036 200 : BS.MatFU = MatFU;
1037 200 : BS.bak = muluu(n, (n-1)*(n-2)); /* safe */
1038 200 : BS.deg = n;
1039 200 : prinfo = gel(csts,8);
1040 :
1041 200 : if (t) x0 = gmul(x0, absisqrtn(rhs, n, Prec));
1042 200 : tmp = divrr(c1,c2);
1043 200 : c3 = mulrr(dbltor(1.39), tmp);
1044 200 : c4 = mulur(n-1, c3);
1045 200 : c14 = mulrr(c4, vecmax_shallow(RgM_sumcol(gabs(A,DEFAULTPREC))));
1046 :
1047 200 : x1 = gmax_shallow(x0, sqrtnr(shiftr(tmp,1),n));
1048 200 : x2 = gmax_shallow(x1, sqrtnr(mulur(10,c14), n));
1049 200 : x3 = gmax_shallow(x2, sqrtnr(shiftr(c14, EXPO1+1),n));
1050 400 : c90 = gmul(shiftr(mulur(18,mppi(DEFAULTPREC)), 5*(4+r)),
1051 200 : gmul(gmul(mpfact(r+3), powiu(muliu(BS.bak,r+2), r+3)),
1052 200 : glog(muliu(BS.bak,2*(r+2)),DEFAULTPREC)));
1053 :
1054 200 : dP = ZX_deriv(P);
1055 200 : Delta0 = RgM_sumcol(A);
1056 :
1057 571 : for (; iroot<=s; iroot++)
1058 : {
1059 389 : GEN Delta = Delta0, Delta2, D, Deps5, MatNE, Hmu, diffRo, c5, c7, ro0;
1060 : long i1, iroot1, iroot2, k;
1061 :
1062 389 : if (iroot <= r) Delta = RgC_add(Delta, RgC_Rg_mul(gel(A,iroot), stoi(-n)));
1063 389 : D = gabs(Delta,Prec); i1 = vecindexmax(D);
1064 389 : c5 = gel(D, i1);
1065 389 : Delta2 = RgC_Rg_div(Delta, gel(Delta, i1));
1066 389 : c5 = myround(gprec_w(c5,DEFAULTPREC), 1);
1067 389 : Deps5 = divrr(subrr(c5,eps5), eps5);
1068 389 : c7 = mulur(r,c5);
1069 389 : BS.c10 = divur(n,c7);
1070 389 : BS.c13 = divur(n,c5);
1071 389 : 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 389 : prec = Prec;
1078 : for (;;)
1079 : {
1080 389 : 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 389 : ro0 = gel(ro,iroot);
1086 389 : BS.ro = ro;
1087 389 : BS.iroot = iroot;
1088 389 : BS.Pi = mppi(prec);
1089 389 : BS.Pi2 = Pi2n(1,prec);
1090 389 : diffRo = cgetg(r+1, t_VEC);
1091 1159 : for (k=1; k<=r; k++)
1092 : {
1093 770 : GEN z = gel(ro,k);
1094 770 : z = (k == iroot)? gdiv(rhs, poleval(dP, z)): gsub(ro0, z);
1095 770 : gel(diffRo,k) = gabs(z, prec);
1096 : }
1097 389 : other_roots(iroot, &iroot1,&iroot2);
1098 389 : BS.divro = gdiv(gsub(ro0, gel(ro,iroot2)), gsub(ro0, gel(ro,iroot1)));
1099 : /* Compute h_1....h_r */
1100 389 : c91 = c90;
1101 1159 : for (k=1; k<=r; k++)
1102 : {
1103 770 : GEN z = gdiv(gcoeff(MatFU,iroot1,k), gcoeff(MatFU,iroot2,k));
1104 770 : z = gmax_shallow(gen_1, abslog(z));
1105 770 : c91 = gmul(c91, gmax_shallow(gel(ALH,k), gdiv(z, BS.bak)));
1106 : }
1107 389 : BS.c91 = c91;
1108 :
1109 4141 : for (; ine<lg(ne); ine++)
1110 : {
1111 3770 : pari_sp av2 = avma;
1112 3770 : long lS = lg(S);
1113 : GEN Lambda, B0, c6, c8;
1114 3770 : GEN NE = gel(MatNE,ine), v = cgetg(r+1,t_COL);
1115 :
1116 3770 : if (DEBUGLEVEL>1) err_printf(" - norm sol. no %ld/%ld\n",ine,lg(ne)-1);
1117 11652 : for (k=1; k<=r; k++)
1118 : {
1119 7882 : GEN z = gdiv(gel(diffRo,k), gabs(gel(NE,k), prec));
1120 7882 : gel(v,k) = glog(z, prec);
1121 : }
1122 3770 : Lambda = RgM_RgC_mul(A,v);
1123 :
1124 3770 : c6 = addrr(dbltor(0.1), vecmax_shallow(gabs(Lambda,DEFAULTPREC)));
1125 3770 : c6 = myround(c6, 1);
1126 3770 : c8 = addrr(dbltor(1.23), mulur(r,c6));
1127 3770 : BS.c11= mulrr(shiftr(c3,1) , mpexp(divrr(mulur(n,c8),c7)));
1128 3770 : BS.c15= mulrr(shiftr(c14,1), mpexp(divrr(mulur(n,c6),c5)));
1129 3770 : BS.NE = NE;
1130 3770 : BS.Hmu = gel(Hmu,ine);
1131 3770 : if (is_pm1(BS.Ind))
1132 : {
1133 1111 : GEN mu = gel(ne,ine), Lmu = cgetg(lg(prinfo),t_VEC);
1134 : long i, j;
1135 :
1136 2229 : for (i = 1; i < lg(prinfo); i++)
1137 : {
1138 1118 : GEN v = gel(prinfo,i), PR = gel(v,3), L = cgetg(4, t_VECSMALL);
1139 4472 : for (j = 1; j <= 3; j++) L[j] = itou(nf_to_Fq(nf, mu, gel(PR,j)));
1140 1118 : gel(Lmu, i) = L;
1141 : }
1142 2213 : if (! (B0 = get_B0(i1, Delta2, Lambda, Deps5, prec, &BS)) ||
1143 1102 : !TrySol(&S, B0, i1, Delta2, Lambda, ro, Lmu, NE,MatFU,prinfo,
1144 : P,rhs))
1145 18 : goto START;
1146 1102 : if (lg(S) == lS) set_avma(av2);
1147 : }
1148 : else
1149 : {
1150 2659 : if (! (Bx = get_Bx_LLL(i1, Delta2, Lambda, prec, &BS)) )
1151 9 : goto START;
1152 2650 : x3 = gc_upto(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, long flag);
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)), 0);
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 x)
1232 749 : { pari_APPLY_same(ne2_to_xy(gel(x,i))); }
1233 :
1234 : static GEN
1235 21 : sol_0(void) { retmkvec( mkvec2(gen_0,gen_0) ); }
1236 : static void
1237 140 : sols_from_R(GEN Rab, GEN *pS, GEN P, GEN POL, GEN rhs)
1238 : {
1239 140 : GEN ry = nfrootsQ(Rab);
1240 140 : long k, l = lg(ry);
1241 280 : for (k = 1; k < l; k++)
1242 140 : if (typ(gel(ry,k)) == t_INT) check_y(pS, P, POL, gel(ry,k), rhs);
1243 140 : }
1244 : static GEN
1245 70 : absZ_factor_if_easy(GEN rhs, GEN x3)
1246 : {
1247 : GEN F, U;
1248 70 : if (expi(rhs) < 150 || expo(x3) >= BITS_IN_LONG) return absZ_factor(rhs);
1249 0 : F = absZ_factor_limit_strict(rhs, 500000, &U);
1250 0 : return U? NULL: F;
1251 : }
1252 : /* Given a tnf structure as returned by thueinit, a RHS and
1253 : * optionally the solutions to the norm equation, returns the solutions to
1254 : * the Thue equation F(x,y)=a */
1255 : GEN
1256 392 : thue(GEN tnf, GEN rhs, GEN ne)
1257 : {
1258 392 : pari_sp av = avma;
1259 : GEN POL, C, L, x3, S;
1260 :
1261 392 : if (typ(tnf) == t_POL) tnf = thueinit(tnf, 0, DEFAULTPREC);
1262 392 : if (!checktnf(tnf)) pari_err_TYPE("thue [please apply thueinit()]", tnf);
1263 392 : if (typ(rhs) != t_INT) pari_err_TYPE("thue",rhs);
1264 392 : if (ne && typ(ne) != t_VEC) pari_err_TYPE("thue",ne);
1265 :
1266 : /* solve P(x,y) = rhs <=> POL(L x, y) = C rhs, with POL monic in Z[X] */
1267 392 : POL = gel(tnf,1);
1268 392 : C = gel(POL,2); rhs = gmul(C, rhs);
1269 392 : if (typ(rhs) != t_INT) retgc_const(av, cgetg(1, t_VEC));
1270 392 : if (!signe(rhs))
1271 : {
1272 28 : GEN v = gel(tnf,2);
1273 28 : set_avma(av);
1274 : /* at least 2 irreducible factors, one of which has degree 1 */
1275 28 : if (typ(v) == t_VECSMALL && v[1] ==1)
1276 7 : pari_err_DOMAIN("thue","#sols","=",strtoGENstr("oo"),rhs);
1277 21 : return sol_0();
1278 : }
1279 364 : L = gel(POL,3);
1280 364 : POL = gel(POL,1);
1281 364 : if (lg(tnf) == 8)
1282 : {
1283 182 : if (!ne)
1284 : {
1285 154 : GEN F = absZ_factor(rhs);
1286 154 : ne = get_ne(tnf_get_bnf(tnf), rhs, F, tnf_get_Ind(tnf));
1287 : }
1288 182 : if (lg(ne) == 1) retgc_const(av, cgetg(1, t_VEC));
1289 182 : S = LargeSols(POL, tnf, rhs, ne);
1290 : }
1291 182 : else if (typ(gel(tnf,3)) == t_REAL)
1292 : { /* Case s=0. All solutions are "small". */
1293 112 : GEN bnf = tnf_get_bnf(tnf);
1294 112 : GEN c0 = gel(tnf,3), F;
1295 112 : x3 = sqrtnr(mulir(absi_shallow(rhs),c0), degpol(POL));
1296 112 : x3 = addrr(x3, dbltor(0.1)); /* guard from round-off errors */
1297 112 : S = cgetg(1,t_VEC);
1298 112 : if (!ne && typ(bnf) == t_VEC && expo(x3) > 10)
1299 : {
1300 70 : F = absZ_factor_if_easy(rhs, x3);
1301 70 : if (F) ne = get_ne(bnf, rhs, F, gen_1);
1302 : }
1303 112 : if (ne)
1304 : {
1305 91 : if (lg(ne) == 1) retgc_const(av, cgetg(1, t_VEC));
1306 77 : if (degpol(POL) == 2) /* quadratic imaginary */
1307 : {
1308 70 : GEN u = NULL;
1309 70 : long w = 2;
1310 70 : if (typ(bnf) == t_VEC)
1311 : {
1312 56 : u = bnf_get_tuU(bnf);
1313 56 : w = bnf_get_tuN(bnf);
1314 : }
1315 : else
1316 : {
1317 14 : GEN D = coredisc(ZX_disc(POL));
1318 14 : if (cmpis(D, -4) >= 0)
1319 : {
1320 14 : GEN F, T = quadpoly_i(D);
1321 14 : w = equalis(D, -4)? 4: 6;
1322 14 : setvarn(T, fetch_var_higher());
1323 14 : F = gcoeff(nffactor(POL, T), 1, 1);
1324 14 : u = gneg(lift_shallow(gel(F,2))); delete_var();
1325 : }
1326 : }
1327 70 : if (w == 4) /* u = I */
1328 28 : ne = shallowconcat(ne, RgXQV_RgXQ_mul(ne,u,POL));
1329 42 : else if (w == 6) /* u = j */
1330 : {
1331 28 : GEN u2 = RgXQ_sqr(u,POL);
1332 28 : ne = shallowconcat1(mkvec3(ne, RgXQV_RgXQ_mul(ne,u,POL),
1333 : RgXQV_RgXQ_mul(ne,u2,POL)));
1334 : }
1335 70 : S = ne2V_to_xyV(ne);
1336 70 : S = filter_sol_Z(S);
1337 70 : S = shallowconcat(S, RgV_neg(S));
1338 : }
1339 : }
1340 98 : if (lg(S) == 1) S = SmallSols(S, x3, POL, rhs);
1341 : }
1342 70 : else if (typ(gel(tnf,3)) == t_INT) /* reducible case, pure power*/
1343 : {
1344 35 : GEN bnf, ne1 = NULL, ne2 = NULL;
1345 35 : long e = itos( gel(tnf,3) );
1346 35 : if (!Z_ispowerall(rhs, e, &rhs)) retgc_const(av, cgetg(1, t_VEC));
1347 28 : tnf = gel(tnf,2);
1348 28 : bnf = tnf_get_bnf(tnf);
1349 28 : ne = get_neabs(bnf, rhs, lg(tnf)==8?tnf_get_Ind(tnf): gen_1);
1350 28 : ne1= bnfisintnorm_i(bnf,rhs,1,ne,0);
1351 28 : S = thue(tnf, rhs, ne1);
1352 28 : if (!odd(e) && lg(tnf)==8) /* if s=0, norms are positive */
1353 : {
1354 7 : ne2 = bnfisintnorm_i(bnf,rhs,-1,ne,0);
1355 7 : S = shallowconcat(S, thue(tnf, negi(rhs), ne2));
1356 : }
1357 : }
1358 : else /* other reducible cases */
1359 : { /* solve f^e * g = rhs, f irreducible factor of smallest degree */
1360 35 : GEN P, D, v = gel(tnf, 2), R = gel(tnf, 3);
1361 35 : long i, l, e = v[2], va = v[3], vb = v[4];
1362 35 : P = cgetg(lg(POL), t_POL); P[1] = POL[1];
1363 35 : D = divisors(rhs); l = lg(D);
1364 35 : S = cgetg(1,t_VEC);
1365 217 : for (i = 1; i < l; i++)
1366 : {
1367 182 : GEN Rab, df = gel(D,i), dg = gel(D,l-i); /* df*dg=|rhs| */
1368 182 : if (e > 1 && !Z_ispowerall(df, e, &df)) continue;
1369 : /* Rab: univariate polynomial in Z[Y], whose roots are the possible y. */
1370 : /* Here and below, Rab != 0 */
1371 70 : if (signe(rhs) < 0) dg = negi(dg); /* df*dg=rhs */
1372 70 : Rab = gsubst(gsubst(R, va, df), vb, dg);
1373 70 : sols_from_R(Rab, &S,P,POL,rhs);
1374 70 : Rab = gsubst(gsubst(R, va, negi(df)), vb, odd(e)? negi(dg): dg);
1375 70 : sols_from_R(Rab, &S,P,POL,rhs);
1376 : }
1377 : }
1378 343 : S = filter_sol_x(S, L);
1379 343 : S = gen_sort_uniq(S, (void*)lexcmp, cmp_nodata);
1380 343 : return gc_upto(av, S);
1381 : }
1382 :
1383 : /********************************************************************/
1384 : /** **/
1385 : /** BNFISINTNORM (K. Belabas) **/
1386 : /** **/
1387 : /********************************************************************/
1388 : struct sol_abs
1389 : {
1390 : GEN rel; /* Primes PR[i] above a, expressed on generators of Cl(K) */
1391 : GEN partrel; /* list of vectors, partrel[i] = rel[1..i] * u[1..i] */
1392 : GEN cyc; /* orders of generators of Cl(K) given in bnf */
1393 :
1394 : long *f; /* f[i] = f(PR[i]/p), inertia degree */
1395 : long *n; /* a = prod p^{ n_p }. n[i]=n_p if PR[i] divides p */
1396 : long *next; /* index of first P above next p, 0 if p is last */
1397 : long *S; /* S[i] = n[i] - sum_{ 1<=k<=i } f[k]*u[k] */
1398 : long *u; /* We want principal ideals I = prod PR[i]^u[i] */
1399 : GEN normsol;/* lists of copies of the u[] which are solutions */
1400 :
1401 : long nPR; /* length(T->rel) = #PR */
1402 : long sindex, smax; /* current index in T->normsol; max. index */
1403 : };
1404 :
1405 : /* u[1..i] has been filled. Norm(u) is correct.
1406 : * Check relations in class group then save it. */
1407 : static void
1408 30639 : test_sol(struct sol_abs *T, long i)
1409 : {
1410 : long k, l;
1411 : GEN s;
1412 :
1413 30639 : if (T->partrel && !ZV_dvd(gel(T->partrel, i), T->cyc)) return;
1414 14546 : if (T->sindex == T->smax)
1415 : { /* no more room in solution list: enlarge */
1416 0 : long new_smax = T->smax << 1;
1417 0 : GEN new_normsol = new_chunk(new_smax+1);
1418 :
1419 0 : for (k=1; k<=T->smax; k++) gel(new_normsol,k) = gel(T->normsol,k);
1420 0 : T->normsol = new_normsol; T->smax = new_smax;
1421 : }
1422 14546 : gel(T->normsol, ++T->sindex) = s = cgetg_copy(T->u, &l);
1423 60389 : for (k=1; k <= i; k++) s[k] = T->u[k];
1424 26453 : for ( ; k < l; k++) s[k] = 0;
1425 14546 : if (DEBUGLEVEL>2)
1426 : {
1427 0 : err_printf("sol = %Ps\n",s);
1428 0 : if (T->partrel) err_printf("T->partrel = %Ps\n",T->partrel);
1429 : }
1430 : }
1431 : /* partrel[i] <-- partrel[i-1] + u[i] * rel[i] */
1432 : static void
1433 21665 : fix_partrel(struct sol_abs *T, long i)
1434 : {
1435 21665 : pari_sp av = avma;
1436 21665 : GEN part1 = gel(T->partrel,i);
1437 21665 : GEN part0 = gel(T->partrel,i-1);
1438 21665 : GEN rel = gel(T->rel, i);
1439 21665 : ulong u = T->u[i];
1440 21665 : long k, l = lg(part1);
1441 64729 : for (k=1; k < l; k++)
1442 43064 : affii(addii(gel(part0,k), muliu(gel(rel,k), u)), gel(part1,k));
1443 21665 : set_avma(av);
1444 21665 : }
1445 :
1446 : /* assume compatible dimensions */
1447 : static void
1448 63273 : ZV_affect(GEN x, GEN y)
1449 : {
1450 63273 : long i, l = lg(y);
1451 166264 : for (i = 1; i < l; i++) affii(gel(x,i), gel(y,i));
1452 63273 : }
1453 :
1454 : /* Recursive loop. Suppose u[1..i] has been filled
1455 : * Find possible solutions u such that, Norm(prod PR[i]^u[i]) = a, taking
1456 : * into account:
1457 : * 1) the relations in the class group if need be.
1458 : * 2) the factorization of a. */
1459 : static void
1460 80122 : isintnorm_loop(struct sol_abs *T, long i)
1461 : {
1462 80122 : if (T->S[i] == 0) /* sum u[i].f[i] = n[i], do another prime */
1463 : {
1464 34146 : long k, next = T->next[i];
1465 34146 : if (next == 0) { test_sol(T, i); return; } /* no primes left */
1466 :
1467 : /* some primes left */
1468 13986 : if (T->partrel) ZV_affect(gel(T->partrel,i), gel(T->partrel, next-1));
1469 21091 : for (k=i+1; k < next; k++) T->u[k] = 0;
1470 13986 : i = next-1;
1471 : }
1472 45976 : else if (i == T->next[i]-2 || i == T->nPR-1)
1473 : { /* only one Prime left above prime; change prime, fix u[i+1] */
1474 : long q;
1475 28455 : if (T->S[i] % T->f[i+1]) return;
1476 14420 : q = T->S[i] / T->f[i+1];
1477 14420 : i++; T->u[i] = q;
1478 14420 : if (T->partrel) fix_partrel(T,i);
1479 14420 : if (T->next[i] == 0) { test_sol(T,i); return; }
1480 : }
1481 :
1482 35448 : i++; T->u[i] = 0;
1483 35448 : if (T->partrel) ZV_affect(gel(T->partrel,i-1), gel(T->partrel,i));
1484 35448 : if (i == T->next[i-1])
1485 : { /* change prime */
1486 31304 : if (T->next[i] == i+1 || i == T->nPR) /* only one Prime above p */
1487 : {
1488 10598 : T->S[i] = 0;
1489 10598 : T->u[i] = T->n[i] / T->f[i]; /* we already know this is exact */
1490 10598 : if (T->partrel) fix_partrel(T, i);
1491 : }
1492 20706 : else T->S[i] = T->n[i];
1493 : }
1494 4144 : else T->S[i] = T->S[i-1]; /* same prime, different Prime */
1495 : for(;;)
1496 : {
1497 66304 : isintnorm_loop(T, i);
1498 66304 : T->S[i] -= T->f[i]; if (T->S[i] < 0) break;
1499 30856 : T->u[i]++;
1500 30856 : if (T->partrel) {
1501 22568 : pari_sp av = avma;
1502 22568 : ZV_affect(ZC_add(gel(T->partrel,i), gel(T->rel,i)), gel(T->partrel,i));
1503 22568 : set_avma(av);
1504 : }
1505 : }
1506 : }
1507 :
1508 : static int
1509 21903 : get_sol_abs(struct sol_abs *T, GEN bnf, GEN nf, GEN fact, GEN *ptPR)
1510 : {
1511 21903 : GEN P = gel(fact,1), E = gel(fact,2), PR;
1512 21903 : long N = nf_get_degree(nf), nP = lg(P)-1, Ngen, max, nPR, i, j;
1513 :
1514 21903 : max = nP*N; /* upper bound for T->nPR */
1515 21903 : T->f = new_chunk(max+1);
1516 21903 : T->n = new_chunk(max+1);
1517 21903 : T->next = new_chunk(max+1);
1518 21903 : *ptPR = PR = cgetg(max+1, t_COL); /* length to be fixed later */
1519 :
1520 21903 : nPR = 0;
1521 55776 : for (i = 1; i <= nP; i++)
1522 : {
1523 41958 : GEN L = idealprimedec(nf, gel(P,i));
1524 41958 : long lL = lg(L), gcd, k, v;
1525 41958 : ulong vn = itou(gel(E,i));
1526 :
1527 : /* check that gcd_{P | p} f_P divides n_p */
1528 41958 : gcd = pr_get_f(gel(L,1));
1529 41986 : for (j=2; gcd > 1 && j < lL; j++) gcd = ugcd(gcd, pr_get_f(gel(L,j)));
1530 41958 : if (gcd > 1 && vn % gcd)
1531 : {
1532 8085 : if (DEBUGLEVEL>2) err_printf("gcd f_P does not divide n_p\n");
1533 8085 : return 0;
1534 : }
1535 33873 : v = (i==nP)? 0: nPR + lL;
1536 93051 : for (k = 1; k < lL; k++)
1537 : {
1538 59178 : GEN pr = gel(L,k);
1539 59178 : gel(PR, ++nPR) = pr;
1540 59178 : T->f[nPR] = pr_get_f(pr) / gcd;
1541 59178 : T->n[nPR] = vn / gcd;
1542 59178 : T->next[nPR] = v;
1543 : }
1544 : }
1545 13818 : T->nPR = nPR;
1546 13818 : setlg(PR, nPR + 1);
1547 :
1548 13818 : T->u = cgetg(nPR+1, t_VECSMALL);
1549 13818 : T->S = new_chunk(nPR+1);
1550 13818 : if (bnf) { T->cyc = bnf_get_cyc(bnf); Ngen = lg(T->cyc)-1; }
1551 119 : else { T->cyc = NULL; Ngen = 0; }
1552 13818 : if (Ngen == 0)
1553 861 : T->rel = T->partrel = NULL; /* trivial Cl(K), no relations to check */
1554 : else
1555 : {
1556 12957 : int triv = 1;
1557 12957 : T->partrel = new_chunk(nPR+1);
1558 12957 : T->rel = new_chunk(nPR+1);
1559 58842 : for (i=1; i <= nPR; i++)
1560 : {
1561 45885 : GEN c = isprincipal(bnf, gel(PR,i));
1562 45885 : gel(T->rel,i) = c;
1563 45885 : if (triv && !ZV_equal0(c)) triv = 0; /* non trivial relations in Cl(K)*/
1564 : }
1565 : /* triv = 1: all ideals dividing a are principal */
1566 12957 : if (triv) T->rel = T->partrel = NULL;
1567 : }
1568 13818 : if (T->partrel)
1569 : {
1570 10640 : long B = ZV_max_lg(T->cyc) + 3;
1571 60130 : for (i = 0; i <= nPR; i++)
1572 : { /* T->partrel[0] also needs to be initialized */
1573 49490 : GEN c = cgetg(Ngen+1, t_COL); gel(T->partrel,i) = c;
1574 134694 : for (j=1; j<=Ngen; j++)
1575 : {
1576 85204 : GEN z = cgeti(B); gel(c,j) = z;
1577 85204 : z[1] = evalsigne(0)|evallgefint(B);
1578 : }
1579 : }
1580 : }
1581 13818 : T->smax = 511;
1582 13818 : T->normsol = new_chunk(T->smax+1);
1583 13818 : T->S[0] = T->n[1];
1584 13818 : T->next[0] = 1;
1585 13818 : T->sindex = 0;
1586 13818 : isintnorm_loop(T, 0); return 1;
1587 : }
1588 :
1589 : /* Return unit of norm -1 (NULL if it doesn't exit). */
1590 : static GEN
1591 3038 : get_unit_1(GEN bnf, long flag)
1592 : {
1593 : GEN v;
1594 : long i;
1595 :
1596 3038 : if (DEBUGLEVEL > 2) err_printf("looking for a fundamental unit of norm -1\n");
1597 3038 : if (odd(nf_get_degree(bnf_get_nf(bnf)))) return gen_m1;
1598 1169 : v = nfsign_fu(bnf, NULL);
1599 1190 : for (i = 1; i < lg(v); i++)
1600 1183 : if (Flv_sum( gel(v,i), 2))
1601 : {
1602 1162 : GEN fu = NULL;
1603 1162 : if (flag)
1604 : {
1605 7 : fu = bnf_build_cheapfu(bnf);
1606 7 : if (!fu) fu = bnf_compactfu(bnf);
1607 : }
1608 1162 : if (!fu) fu = bnf_get_fu(bnf);
1609 1162 : return gel(fu, i);
1610 : }
1611 7 : return NULL;
1612 : }
1613 :
1614 : GEN
1615 21896 : bnfisintnormabs0(GEN bnf, GEN a, long flag)
1616 : {
1617 : struct sol_abs T;
1618 : GEN nf, res, PR, F;
1619 21896 : long i, fl = nf_FORCE | nf_GEN_IF_PRINCIPAL | (flag ? nf_GENMAT: 0);
1620 :
1621 21896 : if ((F = check_arith_all(a,"bnfisintnormabs")))
1622 : {
1623 826 : a = typ(a) == t_VEC? gel(a,1): factorback(F);
1624 826 : if (signe(a) < 0) F = clean_Z_factor(F);
1625 : }
1626 21896 : nf = bnf_get_nf(bnf);
1627 21896 : if (!signe(a)) return mkvec(gen_0);
1628 21882 : if (is_pm1(a)) return mkvec(gen_1);
1629 21777 : if (!F) F = absZ_factor(a);
1630 21777 : if (!get_sol_abs(&T, bnf, nf, F, &PR)) return cgetg(1, t_VEC);
1631 : /* |a| > 1 => T.nPR > 0 */
1632 13699 : res = cgetg(T.sindex+1, t_VEC);
1633 27839 : for (i=1; i<=T.sindex; i++)
1634 : {
1635 14140 : GEN x = vecsmall_to_col( gel(T.normsol,i) );
1636 14140 : x = isprincipalfact(bnf, NULL, PR, x, fl);
1637 14140 : if (!flag) x = nf_to_scalar_or_alg(nf, x);
1638 14140 : gel(res,i) = x; /* solution, up to sign */
1639 : }
1640 13699 : return res;
1641 : }
1642 : GEN
1643 252 : bnfisintnormabs(GEN bnf, GEN a)
1644 252 : { return bnfisintnormabs0(bnf, a, 0); }
1645 :
1646 : /* true nf */
1647 : GEN
1648 161 : ideals_by_norm(GEN nf, GEN a)
1649 : {
1650 : struct sol_abs T;
1651 : GEN res, PR, F;
1652 : long i;
1653 :
1654 161 : if ((F = check_arith_pos(a,"ideals_by_norm")))
1655 0 : a = typ(a) == t_VEC? gel(a,1): factorback(F);
1656 161 : if (is_pm1(a)) return mkvec(trivial_fact());
1657 126 : if (!F) F = absZ_factor(a);
1658 126 : if (!get_sol_abs(&T, NULL, nf, F, &PR)) return cgetg(1, t_VEC);
1659 119 : res = cgetg(T.sindex+1, t_VEC);
1660 525 : for (i=1; i<=T.sindex; i++)
1661 : {
1662 406 : GEN x = vecsmall_to_col( gel(T.normsol,i) );
1663 406 : gel(res,i) = famat_remove_trivial(mkmat2(PR, x));
1664 : }
1665 119 : return res;
1666 : }
1667 :
1668 : /* largest prime used in factorbase */
1669 : static ulong
1670 28 : bnf_get_lastp(GEN bnf)
1671 : {
1672 28 : GEN vbase = gel(bnf,5);
1673 28 : long l = lg(vbase), i;
1674 28 : ulong P = 0;
1675 1393 : for (i = 1; i < l; i++)
1676 : {
1677 1365 : GEN pr = gel(vbase,i);
1678 1365 : ulong p = itou(pr_get_p(pr));
1679 1365 : if (p > P) P = p;
1680 : }
1681 28 : return P;
1682 : }
1683 :
1684 : /* true bnf; z = bnfisintnormabs0(bnf,a,flag), sa = 1 or -1,
1685 : * return bnfisintnorm0(bnf,sa*|a|,flag). If flag is set, allow returning
1686 : * elements in factored form */
1687 : static GEN
1688 21903 : bnfisintnorm_i(GEN bnf, GEN a, long sa, GEN z, long flag)
1689 : {
1690 21903 : GEN nf = bnf_get_nf(bnf), T = nf_get_pol(nf), f = nf_get_index(nf), unit = gen_0;
1691 21903 : GEN Tp, A = signe(a) == sa? a: negi(a);
1692 21903 : long sNx, i, j, N = degpol(T), l = lg(z);
1693 21903 : ulong p, Ap = 0; /* gcc -Wall */
1694 : forprime_t S;
1695 21903 : if (!signe(a)) return z;
1696 21889 : u_forprime_init(&S, flag? bnf_get_lastp(bnf): 3, ULONG_MAX);
1697 30492 : while((p = u_forprime_next(&S)))
1698 30492 : if (umodiu(f,p) && (Ap = umodiu(A,p))) break;
1699 21889 : Tp = ZX_to_Flx(T,p);
1700 : /* p > 2 doesn't divide A nor Q_denom(z in Z_K)*/
1701 :
1702 : /* update z in place to get correct signs: multiply by unit of norm -1 if
1703 : * it exists, otherwise delete solution with wrong sign */
1704 36141 : for (i = j = 1; i < l; i++)
1705 : {
1706 14252 : GEN x = gel(z,i);
1707 14252 : long tx = typ(x);
1708 :
1709 14252 : switch(tx)
1710 : {
1711 13307 : case t_POL:
1712 : {
1713 13307 : GEN dx, y = Q_remove_denom(x,&dx);
1714 13307 : ulong Np = Flx_resultant(Tp, ZX_to_Flx(y,p), p);
1715 13307 : ulong dA = dx? Fl_mul(Ap, Fl_powu(umodiu(dx,p), N, p), p): Ap;
1716 : /* Nx = Res(T,y) / dx^N = A or -A. Check mod p */
1717 13307 : sNx = dA == Np? sa: -sa; break;
1718 : }
1719 273 : case t_MAT:
1720 : {
1721 273 : GEN G = gel(x,1), E = gel(x,2);
1722 273 : long k, lG = lg(G);
1723 273 : GEN g = cgetg(lG, t_VECSMALL), e = cgetg(lG, t_VECSMALL);
1724 : ulong Np;
1725 30849 : for (k = 1; k < lG; k++)
1726 : {
1727 30576 : GEN NGk = nfnorm(nf, gel(G,k));
1728 30576 : (void)Z_lvalrem(NGk, p, &NGk);
1729 30576 : g[k] = umodiu(NGk, p);
1730 30576 : e[k] = umodiu(gel(E,k), p-1);
1731 : }
1732 : /* N.B. p can appear in Norm(G_k), where G = \prod G_k^{e_k},
1733 : * but total v_p(Norm(G)) = v_p(A) = 0 */
1734 273 : Np = Flv_factorback(g, e, p);
1735 273 : sNx = Np == Ap? sa: -sa; break;
1736 : }
1737 672 : default: sNx = gsigne(x) < 0 && odd(N) ? -1 : 1;
1738 : }
1739 14252 : if (sNx != sa)
1740 : {
1741 5054 : if (unit == gen_0) unit = get_unit_1(bnf, flag);
1742 5054 : if (!unit)
1743 : {
1744 77 : if (DEBUGLEVEL > 2) err_printf("%Ps eliminated because of sign\n",x);
1745 77 : continue;
1746 : }
1747 4977 : switch(tx)
1748 : {
1749 4781 : case t_POL: x = (unit == gen_m1)? RgX_neg(x): RgXQ_mul(unit,x,T); break;
1750 14 : case t_MAT: x = famat_mul(x, unit); break;
1751 182 : default: x = (unit == gen_m1)? gneg(x): RgX_Rg_mul(unit,x); break;
1752 : }
1753 : }
1754 14175 : gel(z,j++) = x;
1755 : }
1756 21889 : setlg(z, j); return z;
1757 : }
1758 : GEN
1759 602 : bnfisintnorm(GEN bnf, GEN a)
1760 602 : { return bnfisintnorm0(bnf, a, 0); }
1761 : GEN
1762 21644 : bnfisintnorm0(GEN bnf, GEN a, long flag)
1763 : {
1764 21644 : pari_sp av = avma;
1765 : GEN ne;
1766 21644 : bnf = checkbnf(bnf);
1767 21644 : if (flag < 0 || flag > 1) pari_err_FLAG("bnfisintnorm");
1768 21644 : ne = bnfisintnormabs0(bnf, a, flag);
1769 21644 : switch(typ(a))
1770 : {
1771 602 : case t_VEC: a = gel(a,1); break;
1772 0 : case t_MAT: a = factorback(a); break;
1773 : }
1774 21644 : return gc_GEN(av, bnfisintnorm_i(bnf, a, signe(a), ne, flag));
1775 : }
|
