Line data Source code
1 : /* Copyright (C) 2000-2005 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 : #include "pari.h"
15 : #include "paripriv.h"
16 : /*******************************************************************/
17 : /* */
18 : /* QUADRATIC POLYNOMIAL ASSOCIATED TO A DISCRIMINANT */
19 : /* */
20 : /*******************************************************************/
21 :
22 : void
23 151511 : check_quaddisc(GEN x, long *s, long *pr, const char *f)
24 : {
25 : long r;
26 151511 : if (typ(x) != t_INT) pari_err_TYPE(f,x);
27 151497 : *s = signe(x);
28 151497 : if (Z_issquare(x)) pari_err_DOMAIN(f,"issquare(disc)","=", gen_1,x);
29 151500 : r = mod4(x); if (*s < 0 && r) r = 4 - r;
30 151500 : if (r > 1) pari_err_DOMAIN(f,"disc % 4",">", gen_1,x);
31 151486 : if (pr) *pr = r;
32 151486 : }
33 : void
34 6916 : check_quaddisc_real(GEN x, long *r, const char *f)
35 : {
36 6916 : long sx; check_quaddisc(x, &sx, r, f);
37 6916 : if (sx < 0) pari_err_DOMAIN(f, "disc","<",gen_0,x);
38 6916 : }
39 : void
40 2170 : check_quaddisc_imag(GEN x, long *r, const char *f)
41 : {
42 2170 : long sx; check_quaddisc(x, &sx, r, f);
43 2163 : if (sx > 0) pari_err_DOMAIN(f, "disc",">",gen_0,x);
44 2163 : }
45 :
46 : /* X^2 + b X + c is the canonical quadratic t_POL of discriminant D.
47 : * Dodd is nonzero iff D is odd */
48 : static void
49 963477 : quadpoly_bc(GEN D, long Dodd, GEN *b, GEN *c)
50 : {
51 963477 : if (Dodd)
52 : {
53 865043 : pari_sp av = avma;
54 865043 : *b = gen_m1;
55 865043 : *c = gerepileuptoint(av, shifti(subui(1,D), -2));
56 : }
57 : else
58 : {
59 98434 : *b = gen_0;
60 98434 : *c = shifti(D,-2); togglesign(*c);
61 : }
62 963477 : }
63 : /* X^2 - X - (D-1)/4 or X^2 - D/4 */
64 : static GEN
65 243362 : quadpoly_ii(GEN D, long Dmod4)
66 : {
67 243362 : GEN b, c, y = cgetg(5,t_POL);
68 243362 : y[1] = evalsigne(1) | evalvarn(0);
69 243362 : quadpoly_bc(D, Dmod4, &b,&c);
70 243362 : gel(y,2) = c;
71 243362 : gel(y,3) = b;
72 243362 : gel(y,4) = gen_1; return y;
73 : }
74 : GEN
75 2044 : quadpoly(GEN D)
76 : {
77 : long s, Dmod4;
78 2044 : check_quaddisc(D, &s, &Dmod4, "quadpoly");
79 2037 : return quadpoly_ii(D, Dmod4);
80 : }
81 : GEN /* no checks */
82 241325 : quadpoly_i(GEN D) { return quadpoly_ii(D, Mod4(D)); }
83 :
84 : GEN
85 1036 : quadpoly0(GEN x, long v)
86 : {
87 1036 : GEN T = quadpoly(x);
88 1029 : if (v > 0) setvarn(T, v);
89 1029 : return T;
90 : }
91 :
92 : GEN
93 0 : quadgen(GEN x)
94 0 : { retmkquad(quadpoly(x), gen_0, gen_1); }
95 :
96 : GEN
97 560 : quadgen0(GEN x, long v)
98 : {
99 560 : if (v==-1) v = fetch_user_var("w");
100 560 : retmkquad(quadpoly0(x, v), gen_0, gen_1);
101 : }
102 :
103 : /***********************************************************************/
104 : /** **/
105 : /** BINARY QUADRATIC FORMS **/
106 : /** **/
107 : /***********************************************************************/
108 : static int
109 814212 : is_qfi(GEN q) { return typ(q)==t_QFB && qfb_is_qfi(q); }
110 :
111 : static GEN
112 1806971 : check_qfbext(const char *fun, GEN x)
113 : {
114 1806971 : long t = typ(x);
115 1806971 : if (t == t_QFB) return x;
116 195 : if (t == t_VEC && lg(x)==3)
117 : {
118 196 : GEN q = gel(x,1);
119 196 : if (!is_qfi(q) && typ(gel(x,2))==t_REAL) return q;
120 : }
121 0 : pari_err_TYPE(fun, x);
122 : return NULL;/* LCOV_EXCL_LINE */
123 : }
124 :
125 : static GEN
126 77161 : qfb3(GEN x, GEN y, GEN z)
127 77161 : { retmkqfb(icopy(x), icopy(y), icopy(z), qfb_disc3(x,y,z)); }
128 :
129 : static int
130 23783633 : qfb_equal(GEN x, GEN y)
131 : {
132 23783633 : return equalii(gel(x,1),gel(y,1))
133 1592912 : && equalii(gel(x,2),gel(y,2))
134 25376544 : && equalii(gel(x,3),gel(y,3));
135 : }
136 :
137 : /* valid for t_QFB, qfr3, qfr5; shallow */
138 : static GEN
139 878230 : qfb_inv(GEN x) {
140 878230 : GEN z = shallowcopy(x);
141 878232 : gel(z,2) = negi(gel(z,2));
142 878232 : return z;
143 : }
144 : /* valid for t_QFB, gerepile-safe */
145 : static GEN
146 7 : qfbinv(GEN x)
147 7 : { retmkqfb(icopy(gel(x,1)),negi(gel(x,2)),icopy(gel(x,3)), icopy(gel(x,4))); }
148 :
149 : GEN
150 77189 : Qfb0(GEN a, GEN b, GEN c)
151 : {
152 : GEN q, D;
153 77189 : if (!b)
154 : {
155 28 : if (c) pari_err_TYPE("Qfb",c);
156 21 : if (typ(a) == t_VEC && lg(a) == 4)
157 14 : { b = gel(a,2); c = gel(a,3); a = gel(a,1); }
158 7 : else if (typ(a) == t_POL && degpol(a) == 2)
159 0 : { b = gel(a,3); c = gel(a,2); a = gel(a,4); }
160 7 : else pari_err_TYPE("Qfb",a);
161 : }
162 77161 : else if (!c)
163 7 : pari_err_TYPE("Qfb",b);
164 77168 : if (typ(a)!=t_INT) pari_err_TYPE("Qfb",a);
165 77161 : if (typ(b)!=t_INT) pari_err_TYPE("Qfb",b);
166 77161 : if (typ(c)!=t_INT) pari_err_TYPE("Qfb",c);
167 77161 : q = qfb3(a, b, c); D = qfb_disc(q);
168 77161 : if (signe(D) < 0)
169 42357 : { if (signe(a) < 0) pari_err_IMPL("negative definite t_QFB"); }
170 34804 : else if (Z_issquare(D)) pari_err_DOMAIN("Qfb","issquare(disc)","=", gen_1,q);
171 77154 : return q;
172 : }
173 :
174 : /* composition */
175 : static void
176 27276632 : qfb_sqr(GEN z, GEN x)
177 : {
178 : GEN c, d1, x2, v1, v2, c3, m, p1, r;
179 :
180 27276632 : d1 = bezout(gel(x,2),gel(x,1),&x2, NULL); /* usually 1 */
181 27276383 : c = gel(x,3);
182 27276383 : m = mulii(c,x2);
183 27275979 : if (equali1(d1))
184 20653550 : v1 = v2 = gel(x,1);
185 : else
186 : {
187 6622555 : v1 = diviiexact(gel(x,1),d1);
188 6622554 : v2 = mulii(v1, gcdii(d1,c)); /* = v1 iff x primitive */
189 6622549 : c = mulii(c, d1);
190 : }
191 27276099 : togglesign(m);
192 27276216 : r = modii(m,v2);
193 27276045 : p1 = mulii(r, v1);
194 27276011 : c3 = addii(c, mulii(r,addii(gel(x,2),p1)));
195 27275947 : gel(z,1) = mulii(v1,v2);
196 27275905 : gel(z,2) = addii(gel(x,2), shifti(p1,1));
197 27275973 : gel(z,3) = diviiexact(c3,v2);
198 27276024 : }
199 : /* z <- x * y */
200 : static void
201 65103359 : qfb_comp(GEN z, GEN x, GEN y)
202 : {
203 : GEN n, c, d, y1, v1, v2, c3, m, p1, r;
204 :
205 65103359 : if (x == y) { qfb_sqr(z,x); return; }
206 38620402 : n = shifti(subii(gel(y,2),gel(x,2)), -1);
207 38413425 : v1 = gel(x,1);
208 38413425 : v2 = gel(y,1);
209 38413425 : c = gel(y,3);
210 38413425 : d = bezout(v2,v1,&y1,NULL);
211 38439868 : if (equali1(d))
212 22330103 : m = mulii(y1,n);
213 : else
214 : {
215 16174599 : GEN s = subii(gel(y,2), n);
216 16173669 : GEN x2, y2, d1 = bezout(s,d,&x2,&y2); /* x2 s + y2 (x1 v1 + y1 v2) = d1 */
217 16177730 : if (!equali1(d1))
218 : {
219 7903799 : v1 = diviiexact(v1,d1);
220 7903297 : v2 = diviiexact(v2,d1); /* gcd = 1 iff x or y primitive */
221 7903148 : v1 = mulii(v1, gcdii(c,gcdii(gel(x,3),gcdii(d1,n))));
222 7902577 : c = mulii(c, d1);
223 : }
224 16175964 : m = addii(mulii(mulii(y1,y2),n), mulii(gel(y,3),x2));
225 : }
226 38445555 : togglesign(m);
227 38369751 : r = modii(m, v1);
228 38368458 : p1 = mulii(r, v2);
229 38354714 : c3 = addii(c, mulii(r,addii(gel(y,2),p1)));
230 38345685 : gel(z,1) = mulii(v1,v2);
231 38347067 : gel(z,2) = addii(gel(y,2), shifti(p1,1));
232 38351562 : gel(z,3) = diviiexact(c3,v1);
233 : }
234 :
235 : /* not meant to be efficient */
236 : static GEN
237 84 : qfb_comp_gen(GEN x, GEN y)
238 : {
239 84 : GEN d1 = qfb_disc(x), d2 = qfb_disc(y);
240 84 : GEN a1 = gel(x,1), b1 = gel(x,2), c1 = gel(x,3), n1;
241 84 : GEN a2 = gel(y,1), b2 = gel(y,2), c2 = gel(y,3), n2;
242 84 : GEN cx = content(x), cy = content(y), A, B, C, D, U, m, m2;
243 :
244 84 : if (!is_pm1(cx))
245 : {
246 14 : a1 = diviiexact(a1, cx); b1 = diviiexact(b1, cx);
247 14 : c1 = diviiexact(c1, cx); d1 = diviiexact(d1, sqri(cx));
248 : }
249 84 : if (!is_pm1(cy))
250 : {
251 28 : a2 = diviiexact(a2, cy); c2 = diviiexact(c2, cy);
252 28 : b2 = diviiexact(b2, cy); d2 = diviiexact(d2, sqri(cy));
253 : }
254 84 : D = gcdii(d1, d2); if (signe(d1) < 0) setsigne(D, -1);
255 133 : if (!Z_issquareall(diviiexact(d1, D), &n1) ||
256 84 : !Z_issquareall(diviiexact(d2, D), &n2)) return NULL;
257 49 : A = mulii(a1, n2);
258 49 : B = mulii(a2, n1);
259 49 : C = shifti(addii(mulii(b1, n2), mulii(b2, n1)), -1);
260 49 : U = ZV_extgcd(mkvec3(A, B, C));
261 49 : m = gel(U,1); U = gmael(U,2,3);
262 49 : A = mulii(diviiexact(mulii(a1,b2),m), gel(U,1));
263 49 : B = mulii(diviiexact(mulii(a2,b1),m), gel(U,2));
264 49 : C = addii(mulii(b1,b2), mulii(D, mulii(n1,n2)));
265 49 : C = mulii(diviiexact(shifti(C,-1), m), gel(U,3));
266 49 : B = addii(A, addii(B, C));
267 49 : m2 = sqri(m);
268 49 : A = diviiexact(mulii(a1, a2), m2);
269 49 : C = diviiexact(shifti(subii(sqri(B),D), -2), A);
270 49 : cx = mulii(cx, cy);
271 49 : if (!is_pm1(cx))
272 : {
273 14 : A = mulii(A, cx); B = mulii(B, cx);
274 14 : C = mulii(C, cx); D = mulii(D, sqri(cx));
275 : }
276 49 : return mkqfb(A, B, C, D);
277 : }
278 :
279 : static GEN redimag_av(pari_sp av, GEN q);
280 : static GEN
281 62378261 : qficomp0(GEN x, GEN y, int raw)
282 : {
283 62378261 : pari_sp av = avma;
284 62378261 : GEN z = cgetg(5,t_QFB);
285 62373300 : gel(z,4) = gel(x,4);
286 62373300 : qfb_comp(z, x,y);
287 62094197 : if (raw) return gerepilecopy(av,z);
288 62092426 : return redimag_av(av, z);
289 : }
290 : static GEN redreal(GEN x);
291 : static GEN
292 441 : qfrcomp0(GEN x, GEN y, int raw)
293 : {
294 441 : pari_sp av = avma;
295 441 : GEN dx = NULL, dy = NULL;
296 441 : GEN z = cgetg(5,t_QFB);
297 441 : if (typ(x)==t_VEC) { dx = gel(x,2); x = gel(x,1); }
298 441 : if (typ(y)==t_VEC) { dy = gel(y,2); y = gel(y,1); }
299 441 : gel(z,4) = gel(x,4);
300 441 : qfb_comp(z, x,y);
301 441 : if (dx) z = mkvec2(z, dy? addrr(dx, dy): dx); else if (dy) z = mkvec2(z, dy);
302 441 : if (!raw) z = redreal(z);
303 441 : return gerepilecopy(av, z);
304 : }
305 : /* same discriminant, no distance, no checks */
306 : GEN
307 27252157 : qfbcomp_i(GEN x, GEN y)
308 27252157 : { return qfb_is_qfi(x)? qficomp0(x,y,0): qfrcomp0(x,y,0); }
309 : GEN
310 129398 : qfbcomp(GEN x, GEN y)
311 : {
312 129398 : GEN qx = check_qfbext("qfbcomp", x);
313 129398 : GEN qy = check_qfbext("qfbcomp", y);
314 129398 : if (!equalii(gel(qx,4),gel(qy,4)))
315 : {
316 63 : pari_sp av = avma;
317 63 : GEN z = qfb_comp_gen(qx, qy);
318 63 : if (typ(x) == t_VEC || typ(y) == t_VEC)
319 7 : pari_err_IMPL("Shanks's distance in general composition");
320 56 : if (!z) pari_err_OP("*",x,y);
321 21 : return gerepileupto(av, qfbred(z));
322 : }
323 129335 : return qfb_is_qfi(qx)? qficomp0(x,y,0): qfrcomp0(x,y,0);
324 : }
325 : /* same discriminant, no distance, no checks */
326 : GEN
327 0 : qfbcompraw_i(GEN x, GEN y)
328 0 : { return qfb_is_qfi(x)? qficomp0(x,y,1): qfrcomp0(x,y,1); }
329 : GEN
330 2191 : qfbcompraw(GEN x, GEN y)
331 : {
332 2191 : GEN qx = check_qfbext("qfbcompraw", x);
333 2191 : GEN qy = check_qfbext("qfbcompraw", y);
334 2191 : if (!equalii(gel(qx,4),gel(qy,4)))
335 : {
336 21 : pari_sp av = avma;
337 21 : GEN z = qfb_comp_gen(qx, qy);
338 21 : if (typ(x) == t_VEC || typ(y) == t_VEC)
339 0 : pari_err_IMPL("Shanks's distance in general composition");
340 21 : if (!z) pari_err_OP("qfbcompraw",x,y);
341 21 : return gerepilecopy(av, z);
342 : }
343 2170 : if (!equalii(gel(qx,4),gel(qy,4))) pari_err_OP("qfbcompraw",x,y);
344 2170 : return qfb_is_qfi(qx)? qficomp0(x,y,1): qfrcomp0(x,y,1);
345 : }
346 :
347 : static GEN
348 793666 : qfisqr0(GEN x, long raw)
349 : {
350 793666 : pari_sp av = avma;
351 793666 : GEN z = cgetg(5,t_QFB);
352 793666 : gel(z,4) = gel(x,4);
353 793666 : qfb_sqr(z,x);
354 793666 : if (raw) return gerepilecopy(av,z);
355 793666 : return redimag_av(av, z);
356 : }
357 : static GEN
358 14 : qfrsqr0(GEN x, long raw)
359 : {
360 14 : pari_sp av = avma;
361 14 : GEN dx = NULL, z = cgetg(5,t_QFB);
362 14 : if (typ(x) == t_VEC) { dx = gel(x,2); x = gel(x,1); }
363 14 : gel(z,4) = gel(x,4); qfb_sqr(z,x);
364 14 : if (dx) z = mkvec2(z, shiftr(dx,1));
365 14 : if (!raw) z = redreal(z);
366 14 : return gerepilecopy(av, z);
367 : }
368 : /* same discriminant, no distance, no checks */
369 : GEN
370 698063 : qfbsqr_i(GEN x)
371 698063 : { return qfb_is_qfi(x)? qfisqr0(x,0): qfrsqr0(x,0); }
372 : GEN
373 95617 : qfbsqr(GEN x)
374 : {
375 95617 : GEN qx = check_qfbext("qfbsqr", x);
376 95617 : return qfb_is_qfi(qx)? qfisqr0(x,0): qfrsqr0(x,0);
377 : }
378 :
379 : static GEN
380 6867 : qfr_1_by_disc(GEN D)
381 : {
382 : GEN y, r, s;
383 6867 : check_quaddisc_real(D, NULL, "qfr_1_by_disc");
384 6867 : y = cgetg(5,t_QFB);
385 6867 : s = sqrtremi(D, &r); togglesign(r); /* s^2 - r = D */
386 6867 : if (mpodd(r))
387 : {
388 3535 : s = subiu(s,1);
389 3535 : r = subii(r, addiu(shifti(s, 1), 1));
390 3535 : r = shifti(r, -2); set_avma((pari_sp)y); s = icopy(s);
391 : }
392 : else
393 3332 : { r = shifti(r, -2); set_avma((pari_sp)s); }
394 6867 : gel(y,1) = gen_1;
395 6867 : gel(y,2) = s;
396 6867 : gel(y,3) = icopy(r);
397 6867 : gel(y,4) = icopy(D); return y;
398 : }
399 :
400 : static GEN
401 35 : qfr_disc(GEN x)
402 35 : { return qfb_disc(typ(x)==t_VEC ? gel(x,1): x); }
403 :
404 : static GEN
405 35 : qfr_1(GEN x)
406 35 : { return qfr_1_by_disc(qfr_disc(x)); }
407 :
408 : static void
409 0 : qfr_1_fill(GEN y, struct qfr_data *S)
410 : {
411 0 : pari_sp av = avma;
412 0 : GEN y2 = S->isqrtD;
413 0 : gel(y,1) = gen_1;
414 0 : if (mod2(S->D) != mod2(y2)) y2 = subiu(y,1);
415 0 : gel(y,2) = y2; av = avma;
416 0 : gel(y,3) = gerepileuptoint(av, shifti(subii(sqri(y2), S->D),-2));
417 0 : }
418 : static GEN
419 0 : qfr5_1(struct qfr_data *S, long prec)
420 : {
421 0 : GEN y = cgetg(6, t_VEC);
422 0 : qfr_1_fill(y, S);
423 0 : gel(y,4) = gen_0;
424 0 : gel(y,5) = real_1(prec); return y;
425 : }
426 : static GEN
427 0 : qfr3_1(struct qfr_data *S)
428 : {
429 0 : GEN y = cgetg(4, t_VEC);
430 0 : qfr_1_fill(y, S); return y;
431 : }
432 :
433 : /* Assume D < 0 is the discriminant of a t_QFB */
434 : static GEN
435 720115 : qfi_1_by_disc(GEN D)
436 : {
437 720115 : GEN b,c, y = cgetg(5,t_QFB);
438 720115 : quadpoly_bc(D, mod2(D), &b,&c);
439 720115 : if (b == gen_m1) b = gen_1;
440 720115 : gel(y,1) = gen_1;
441 720115 : gel(y,2) = b;
442 720115 : gel(y,3) = c;
443 720115 : gel(y,4) = icopy(D); return y;
444 : }
445 : static GEN
446 708157 : qfi_1(GEN x)
447 : {
448 708157 : if (typ(x) != t_QFB) pari_err_TYPE("qfi_1",x);
449 708157 : return qfi_1_by_disc(qfb_disc(x));
450 : }
451 :
452 : GEN
453 0 : qfb_1(GEN x) { return qfb_is_qfi(x) ? qfi_1(x): qfr_1(x); }
454 :
455 : static GEN
456 9584153 : _qfimul(void *E, GEN x, GEN y) { (void) E; return qficomp0(x,y,0); }
457 : static GEN
458 25411373 : _qfisqr(void *E, GEN x) { (void) E; return qficomp0(x,x,0); }
459 : static GEN
460 7 : _qfimulraw(void *E, GEN x, GEN y) { (void) E; return qficomp0(x,y,1); }
461 : static GEN
462 7 : _qfisqrraw(void *E, GEN x) { (void) E; return qficomp0(x,x,1); }
463 :
464 : static GEN
465 7 : qfipowraw(GEN x, long n)
466 : {
467 7 : pari_sp av = avma;
468 : GEN y;
469 7 : if (!n) return qfi_1(x);
470 7 : if (n== 1) return gcopy(x);
471 7 : if (n==-1) { x = gcopy(x); togglesign(gel(x,2)); return x; }
472 7 : if (n < 0) x = qfb_inv(x);
473 7 : y = gen_powu(x, labs(n), NULL, &_qfisqrraw, &_qfimulraw);
474 7 : return gerepilecopy(av,y);
475 : }
476 :
477 : static GEN
478 11475824 : qfipow(GEN x, GEN n)
479 : {
480 11475824 : pari_sp av = avma;
481 : GEN y;
482 11475824 : long s = signe(n);
483 11475824 : if (!s) return qfi_1(x);
484 10767667 : if (s < 0) x = qfb_inv(x);
485 10767669 : y = gen_pow(qfbred_i(x), n, NULL, &_qfisqr, &_qfimul);
486 10767673 : return gerepilecopy(av,y);
487 : }
488 :
489 : static long
490 412328 : parteucl(GEN L, GEN *d, GEN *v3, GEN *v, GEN *v2)
491 : {
492 : long z;
493 412328 : *v = gen_0; *v2 = gen_1;
494 4351417 : for (z=0; abscmpii(*v3,L) > 0; z++)
495 : {
496 3939089 : GEN t3, t2 = subii(*v, mulii(truedvmdii(*d,*v3,&t3),*v2));
497 3939089 : *v = *v2; *d = *v3; *v2 = t2; *v3 = t3;
498 : }
499 412328 : return z;
500 : }
501 :
502 : /* composition: Shanks' NUCOMP & NUDUPL */
503 : /* L = floor((|d|/4)^(1/4)) */
504 : GEN
505 400722 : nucomp(GEN x, GEN y, GEN L)
506 : {
507 400722 : pari_sp av = avma;
508 : long z;
509 : GEN a, a1, a2, b2, b, d, d1, g, n, p1, q1, q2, s, u, u1, v, v2, v3, Q;
510 :
511 400722 : if (x==y) return nudupl(x,L);
512 400680 : if (!is_qfi(x)) pari_err_TYPE("nucomp",x);
513 400680 : if (!is_qfi(y)) pari_err_TYPE("nucomp",y);
514 :
515 400680 : if (abscmpii(gel(x,1),gel(y,1)) < 0) swap(x, y);
516 400680 : s = shifti(addii(gel(x,2),gel(y,2)), -1);
517 400680 : n = subii(gel(y,2), s);
518 400680 : a1 = gel(x,1);
519 400680 : a2 = gel(y,1); d = bezout(a2,a1,&u,&v);
520 400680 : if (equali1(d)) { a = negi(mulii(u,n)); d1 = d; }
521 163576 : else if (dvdii(s,d)) /* d | s */
522 : {
523 83503 : a = negi(mulii(u,n)); d1 = d;
524 83503 : a1 = diviiexact(a1, d1);
525 83503 : a2 = diviiexact(a2, d1);
526 83503 : s = diviiexact(s, d1);
527 : }
528 : else
529 : {
530 : GEN p2, l;
531 80073 : d1 = bezout(s,d,&u1,NULL);
532 80073 : if (!equali1(d1))
533 : {
534 2044 : a1 = diviiexact(a1,d1);
535 2044 : a2 = diviiexact(a2,d1);
536 2044 : s = diviiexact(s,d1);
537 2044 : d = diviiexact(d,d1);
538 : }
539 80073 : p1 = remii(gel(x,3),d);
540 80073 : p2 = remii(gel(y,3),d);
541 80073 : l = modii(mulii(negi(u1), addii(mulii(u,p1),mulii(v,p2))), d);
542 80073 : a = subii(mulii(l,diviiexact(a1,d)), mulii(u,diviiexact(n,d)));
543 : }
544 400680 : a = modii(a,a1); p1 = subii(a,a1); if (abscmpii(a,p1) > 0) a = p1;
545 400680 : d = a1; v3 = a; z = parteucl(L, &d,&v3, &v,&v2);
546 400680 : Q = cgetg(5,t_QFB);
547 400680 : if (!z) {
548 37632 : g = diviiexact(addii(mulii(v3,s),gel(y,3)), d);
549 37632 : b = a2;
550 37632 : b2 = gel(y,2);
551 37632 : v2 = d1;
552 37632 : gel(Q,1) = mulii(d,b);
553 : } else {
554 : GEN e, q3, q4;
555 363048 : if (z&1) { v3 = negi(v3); v2 = negi(v2); }
556 363048 : b = diviiexact(addii(mulii(a2,d), mulii(n,v)), a1);
557 363048 : e = diviiexact(addii(mulii(s,d),mulii(gel(y,3),v)), a1);
558 363048 : q3 = mulii(e,v2);
559 363048 : q4 = subii(q3,s);
560 363048 : b2 = addii(q3,q4);
561 363048 : g = diviiexact(q4,v);
562 363048 : if (!equali1(d1)) { v2 = mulii(d1,v2); v = mulii(d1,v); b2 = mulii(d1,b2); }
563 363048 : gel(Q,1) = addii(mulii(d,b), mulii(e,v));
564 : }
565 400680 : q1 = mulii(b, v3);
566 400680 : q2 = addii(q1,n);
567 400680 : gel(Q,2) = addii(b2, z? addii(q1,q2): shifti(q1, 1));
568 400680 : gel(Q,3) = addii(mulii(v3,diviiexact(q2,d)), mulii(g,v2));
569 400680 : gel(Q,4) = gel(x,4);
570 400680 : return redimag_av(av, Q);
571 : }
572 :
573 : GEN
574 11648 : nudupl(GEN x, GEN L)
575 : {
576 11648 : pari_sp av = avma;
577 : long z;
578 : GEN u, v, d, d1, p1, a, b, c, a2, b2, c2, Q, v2, v3, g;
579 :
580 11648 : if (!is_qfi(x)) pari_err_TYPE("nudupl",x);
581 11648 : a = gel(x,1);
582 11648 : b = gel(x,2);
583 11648 : d1 = bezout(b,a, &u,NULL);
584 11648 : if (!equali1(d1))
585 : {
586 4620 : a = diviiexact(a, d1);
587 4620 : b = diviiexact(b, d1);
588 : }
589 11648 : c = modii(negi(mulii(u,gel(x,3))), a);
590 11648 : p1 = subii(c,a); if (abscmpii(c,p1) > 0) c = p1;
591 11648 : d = a; v3 = c; z = parteucl(L, &d,&v3, &v,&v2);
592 11648 : a2 = sqri(d);
593 11648 : c2 = sqri(v3);
594 11648 : Q = cgetg(5,t_QFB);
595 11648 : if (!z) {
596 1281 : g = diviiexact(addii(mulii(v3,b),gel(x,3)), d);
597 1281 : b2 = gel(x,2);
598 1281 : v2 = d1;
599 1281 : gel(Q,1) = a2;
600 : } else {
601 : GEN e;
602 10367 : if (z&1) { v = negi(v); d = negi(d); }
603 10367 : e = diviiexact(addii(mulii(gel(x,3),v), mulii(b,d)), a);
604 10367 : g = diviiexact(subii(mulii(e,v2), b), v);
605 10367 : b2 = addii(mulii(e,v2), mulii(v,g));
606 10367 : if (!equali1(d1)) { b2 = mulii(d1,b2); v = mulii(d1,v); v2 = mulii(d1,v2); }
607 10367 : gel(Q,1) = addii(a2, mulii(e,v));
608 : }
609 11648 : gel(Q,2) = addii(b2, subii(sqri(addii(d,v3)), addii(a2,c2)));
610 11648 : gel(Q,3) = addii(c2, mulii(g,v2));
611 11648 : gel(Q,4) = gel(x,4);
612 11648 : return redimag_av(av, Q);
613 : }
614 :
615 : static GEN
616 4739 : mul_nucomp(void *l, GEN x, GEN y) { return nucomp(x, y, (GEN)l); }
617 : static GEN
618 11606 : mul_nudupl(void *l, GEN x) { return nudupl(x, (GEN)l); }
619 : GEN
620 1008 : nupow(GEN x, GEN n, GEN L)
621 : {
622 : pari_sp av;
623 : GEN y, D;
624 :
625 1008 : if (typ(n) != t_INT) pari_err_TYPE("nupow",n);
626 1008 : if (!is_qfi(x)) pari_err_TYPE("nupow",x);
627 1008 : if (gequal1(n)) return gcopy(x);
628 1008 : av = avma;
629 1008 : D = qfb_disc(x);
630 1008 : y = qfi_1_by_disc(D);
631 1008 : if (!signe(n)) return y;
632 959 : if (!L) L = sqrtnint(absi_shallow(D), 4);
633 959 : y = gen_pow_i(x, n, (void*)L, &mul_nudupl, &mul_nucomp);
634 959 : if (signe(n) < 0
635 35 : && !absequalii(gel(y,1),gel(y,2))
636 35 : && !absequalii(gel(y,1),gel(y,3))) togglesign(gel(y,2));
637 959 : return gerepilecopy(av, y);
638 : }
639 :
640 : /* Reduction */
641 :
642 : /* assume a > 0. Write b = q*2a + r, with -a < r <= a */
643 : static GEN
644 8066727 : dvmdii_round(GEN b, GEN a, GEN *r)
645 : {
646 8066727 : GEN a2 = shifti(a, 1), q = dvmdii(b, a2, r);
647 8066721 : if (signe(b) >= 0) {
648 4594073 : if (abscmpii(*r, a) > 0) { q = addiu(q, 1); *r = subii(*r, a2); }
649 : } else { /* r <= 0 */
650 3472648 : if (abscmpii(*r, a) >= 0){ q = subiu(q, 1); *r = addii(*r, a2); }
651 : }
652 8066709 : return q;
653 : }
654 : /* Assume 0 < a <= LONG_MAX. Ensure no overflow */
655 : static long
656 97108332 : dvmdsu_round(long b, ulong a, long *r)
657 : {
658 97108332 : ulong a2 = a << 1, q, ub, ur;
659 97108332 : if (b >= 0) {
660 61930000 : ub = b;
661 61930000 : q = ub / a2;
662 61930000 : ur = ub % a2;
663 61930000 : if (ur > a) { ur -= a; q++; *r = (long)ur; *r -= (long)a; }
664 21808475 : else *r = (long)ur;
665 61930000 : return (long)q;
666 : } else { /* r <= 0 */
667 35178332 : ub = (ulong)-b; /* |b| */
668 35178332 : q = ub / a2;
669 35178332 : ur = ub % a2;
670 35178332 : if (ur >= a) { ur -= a; q++; *r = (long)ur; *r = (long)a - *r; }
671 19236011 : else *r = -(long)ur;
672 35178332 : return -(long)q;
673 : }
674 : }
675 : /* reduce b mod 2*a. Update b,c */
676 : static void
677 2768824 : REDB(GEN a, GEN *b, GEN *c)
678 : {
679 2768824 : GEN r, q = dvmdii_round(*b, a, &r);
680 2768824 : if (!signe(q)) return;
681 2699896 : *c = subii(*c, mulii(q, shifti(addii(*b, r),-1)));
682 2699896 : *b = r;
683 : }
684 : /* Assume a > 0. Reduce b mod 2*a. Update b,c */
685 : static void
686 97109394 : sREDB(ulong a, long *b, ulong *c)
687 : {
688 : long r, q;
689 : ulong uz;
690 105324759 : if (a > LONG_MAX) return; /* b already reduced */
691 97109394 : q = dvmdsu_round(*b, a, &r);
692 97133317 : if (q == 0) return;
693 : /* Final (a,r,c2) satisfies |r| <= |b| hence c2 <= c, c2 = c - q*z,
694 : * where z = (b+r) / 2, representable as long, has the same sign as q. */
695 88917952 : if (*b < 0)
696 : { /* uz = -z >= 0, q < 0 */
697 31078589 : if (r >= 0) /* different signs=>no overflow, exact division */
698 15997152 : uz = (ulong)-((*b + r)>>1);
699 : else
700 : {
701 15081437 : ulong ub = (ulong)-*b, ur = (ulong)-r;
702 15081437 : uz = (ub + ur) >> 1;
703 : }
704 31078589 : *c -= (-q) * uz; /* c -= qz */
705 : }
706 : else
707 : { /* uz = z >= 0, q > 0 */
708 57839363 : if (r <= 0)
709 40183437 : uz = (*b + r)>>1;
710 : else
711 : {
712 17655926 : ulong ub = (ulong)*b, ur = (ulong)r;
713 17655926 : uz = ((ub + ur) >> 1);
714 : }
715 57839363 : *c -= q * uz; /* c -= qz */
716 : }
717 88917952 : *b = r;
718 : }
719 : static void
720 5297903 : REDBU(GEN a, GEN *b, GEN *c, GEN u1, GEN *u2)
721 : { /* REDB(a,b,c) */
722 5297903 : GEN r, q = dvmdii_round(*b, a, &r);
723 5297882 : *c = subii(*c, mulii(q, shifti(addii(*b, r),-1)));
724 5297873 : *b = r;
725 5297873 : *u2 = subii(*u2, mulii(q, u1));
726 5297885 : }
727 :
728 : /* q t_QFB, return reduced representative and set base change U in Sl2(Z) */
729 : GEN
730 1933542 : redimagsl2(GEN q, GEN *U)
731 : {
732 1933542 : pari_sp av = avma;
733 : GEN z, u1,u2,v1,v2,Q;
734 1933542 : GEN a = gel(q,1), b = gel(q,2), c = gel(q,3);
735 : long cmp;
736 1933542 : u1 = gen_1; u2 = gen_0;
737 1933542 : cmp = abscmpii(a, b);
738 1933544 : if (cmp < 0)
739 879483 : REDBU(a,&b,&c, u1,&u2);
740 1054061 : else if (cmp == 0 && signe(b) < 0)
741 : { /* b = -a */
742 13 : b = negi(b);
743 13 : u2 = gen_1;
744 : }
745 : for(;;)
746 : {
747 6351958 : cmp = abscmpii(a, c); if (cmp <= 0) break;
748 4418400 : swap(a,c); b = negi(b);
749 4418424 : z = u1; u1 = u2; u2 = negi(z);
750 4418420 : REDBU(a,&b,&c, u1,&u2);
751 4418415 : if (gc_needed(av, 1)) {
752 0 : if (DEBUGMEM>1) pari_warn(warnmem, "redimagsl2");
753 0 : gerepileall(av, 5, &a,&b,&c, &u1,&u2);
754 : }
755 : }
756 1933542 : if (cmp == 0 && signe(b) < 0)
757 : {
758 15603 : b = negi(b);
759 15603 : z = u1; u1 = u2; u2 = negi(z);
760 : }
761 : /* Let q = (A,B,C). q o [u1,u2; v1,v2] = Q implies
762 : * [v1,v2] = (1/C) [(b-B)/2 u1 - a u2, c u1 - (b+B)/2 u2] */
763 1933542 : z = shifti(subii(b, gel(q,2)), -1);
764 1933541 : v1 = subii(mulii(z, u1), mulii(a, u2)); v1 = diviiexact(v1, gel(q,3));
765 1933543 : z = subii(z, b);
766 1933538 : v2 = addii(mulii(z, u2), mulii(c, u1)); v2 = diviiexact(v2, gel(q,3));
767 1933536 : *U = mkmat2(mkcol2(u1,v1), mkcol2(u2,v2));
768 1933547 : Q = lg(q)==5 ? mkqfb(a,b,c,gel(q,4)): mkvec3(a,b,c);
769 1933547 : return gc_all(av, 2, &Q, U);
770 : }
771 :
772 : static GEN
773 1068755 : setq_b0(ulong a, ulong c, GEN D)
774 1068755 : { retmkqfb(utoipos(a), gen_0, utoipos(c), icopy(D)); }
775 : /* assume |sb| = 1 */
776 : static GEN
777 80262301 : setq(ulong a, ulong b, ulong c, long sb, GEN D)
778 80262301 : { retmkqfb(utoipos(a), sb==1? utoipos(b): utoineg(b), utoipos(c), icopy(D)); }
779 : /* 0 < a, c < 2^BIL, b = 0 */
780 : static GEN
781 959190 : redimag_1_b0(ulong a, ulong c, GEN D)
782 959190 : { return (a <= c)? setq_b0(a, c, D): setq_b0(c, a, D); }
783 :
784 : /* 0 < a, c < 2^BIL: single word affair */
785 : static GEN
786 81428558 : redimag_1(pari_sp av, GEN a, GEN b, GEN c, GEN D)
787 : {
788 : ulong ua, ub, uc;
789 : long sb;
790 : for(;;)
791 140729 : { /* at most twice */
792 81428558 : long lb = lgefint(b); /* <= 3 after first loop */
793 81428558 : if (lb == 2) return redimag_1_b0(a[2],c[2], D);
794 80469368 : if (lb == 3 && uel(b,2) <= (ulong)LONG_MAX) break;
795 141121 : REDB(a,&b,&c);
796 140142 : if (uel(a,2) <= uel(c,2))
797 : { /* lg(b) <= 3 but may be too large for itos */
798 0 : long s = signe(b);
799 0 : set_avma(av);
800 0 : if (!s) return redimag_1_b0(a[2], c[2], D);
801 0 : if (a[2] == c[2]) s = 1;
802 0 : return setq(a[2], b[2], c[2], s, D);
803 : }
804 140142 : swap(a,c); b = negi(b);
805 : }
806 : /* b != 0 */
807 80328247 : set_avma(av);
808 80334989 : ua = a[2];
809 80334989 : ub = sb = b[2]; if (signe(b) < 0) sb = -sb;
810 80334989 : uc = c[2];
811 80334989 : if (ua < ub)
812 29131334 : sREDB(ua, &sb, &uc);
813 51203655 : else if (ua == ub && sb < 0) sb = (long)ub;
814 148393407 : while(ua > uc)
815 : {
816 68015153 : lswap(ua,uc); sb = -sb;
817 68015153 : sREDB(ua, &sb, &uc);
818 : }
819 80378254 : if (!sb) return setq_b0(ua, uc, D);
820 : else
821 : {
822 80268687 : long s = 1;
823 80268687 : if (sb < 0)
824 : {
825 30107368 : sb = -sb;
826 30107368 : if (ua != uc) s = -1;
827 : }
828 80268687 : return setq(ua, sb, uc, s, D);
829 : }
830 : }
831 :
832 : static GEN
833 81234665 : redimag_av(pari_sp av, GEN q)
834 : {
835 81234665 : GEN a = gel(q,1), b = gel(q,2), c = gel(q,3), D = gel(q,4);
836 81234665 : long cmp, lc = lgefint(c);
837 :
838 81234665 : if (lgefint(a) == 3 && lc == 3) return redimag_1(av, a, b, c, D);
839 906428 : cmp = abscmpii(a, b);
840 914075 : if (cmp < 0)
841 434190 : REDB(a,&b,&c);
842 479885 : else if (cmp == 0 && signe(b) < 0)
843 30 : b = negi(b);
844 : for(;;)
845 : {
846 3108567 : cmp = abscmpii(a, c); if (cmp <= 0) break;
847 2918664 : lc = lgefint(a); /* lg(future c): we swap a & c next */
848 2918664 : if (lc == 3) return redimag_1(av, a, b, c, D);
849 2194492 : swap(a,c); b = negi(b); /* apply rho */
850 2194492 : REDB(a,&b,&c);
851 : }
852 189903 : if (cmp == 0 && signe(b) < 0) b = negi(b);
853 189903 : return gerepilecopy(av, mkqfb(a, b, c, D));
854 : }
855 : static GEN
856 17936044 : redimag(GEN q) { return redimag_av(avma, q); }
857 :
858 : static GEN
859 7 : rhoimag(GEN x)
860 : {
861 7 : pari_sp av = avma;
862 7 : GEN a = gel(x,1), b = gel(x,2), c = gel(x,3);
863 7 : int fl = abscmpii(a, c);
864 7 : if (fl <= 0)
865 : {
866 7 : int fg = abscmpii(a, b);
867 7 : if (fg >= 0)
868 : {
869 7 : x = gcopy(x);
870 7 : if ((!fl || !fg) && signe(gel(x,2)) < 0) setsigne(gel(x,2), 1);
871 7 : return x;
872 : }
873 : }
874 0 : swap(a,c); b = negi(b);
875 0 : REDB(a, &b, &c);
876 0 : return gerepilecopy(av, mkqfb(a,b,c, qfb_disc(x)));
877 : }
878 :
879 : /* qfr3 / qfr5 */
880 :
881 : /* t_QFB are unusable: D, sqrtD, isqrtD are recomputed all the time and the
882 : * logarithmic Shanks's distance is costly and hard to control.
883 : * qfr3 / qfr5 routines take a container of t_INTs (e.g a t_VEC) as argument,
884 : * at least 3 (resp. 5) components [it is a feature that they do not check the
885 : * precise type or length of the input]. They return a vector of length 3
886 : * (resp. 5). A qfr3 [a,b,c] contains the form coeffs, in a qfr5 [a,b,c, e,d]
887 : * the t_INT e is a binary exponent, d a t_REAL, coding the distance in
888 : * multiplicative form: the true distance is obtained from qfr5_dist.
889 : * All other qfr routines are obsolete (inefficient) wrappers */
890 :
891 : /* static functions are not stack-clean. Unless mentionned otherwise, public
892 : * functions are. */
893 :
894 : #define EMAX 22
895 : static void
896 12022738 : fix_expo(GEN x)
897 : {
898 12022738 : if (expo(gel(x,5)) >= (1L << EMAX)) {
899 0 : gel(x,4) = addiu(gel(x,4), 1);
900 0 : shiftr_inplace(gel(x,5), - (1L << EMAX));
901 : }
902 12022738 : }
903 :
904 : /* (1/2) log (d * 2^{e * 2^EMAX}). Not stack clean if e != 0 */
905 : GEN
906 184254 : qfr5_dist(GEN e, GEN d, long prec)
907 : {
908 184254 : GEN t = logr_abs(d);
909 184254 : if (signe(e)) {
910 0 : GEN u = mulir(e, mplog2(prec));
911 0 : shiftr_inplace(u, EMAX); t = addrr(t, u);
912 : }
913 184254 : shiftr_inplace(t, -1); return t;
914 : }
915 :
916 : static void
917 17290519 : rho_get_BC(GEN *B, GEN *C, GEN b, GEN c, struct qfr_data *S)
918 : {
919 : GEN t, u;
920 17290519 : u = shifti(c,1);
921 17290519 : t = (abscmpii(S->isqrtD,c) >= 0)? S->isqrtD: c;
922 17290519 : u = remii(addii_sign(t,1, b,signe(b)), u);
923 17290519 : *B = addii_sign(t, 1, u, -signe(u)); /* |t| - (|t|+b) % |2c| */
924 17290519 : if (*B == gen_0)
925 3413 : { u = shifti(S->D, -2); setsigne(u, -1); }
926 : else
927 17287106 : u = shifti(addii_sign(sqri(*B),1, S->D,-1), -2);
928 17290519 : *C = diviiexact(u, c); /* = (B^2-D)/4c */
929 17290519 : }
930 : /* Not stack-clean */
931 : GEN
932 6994478 : qfr3_rho(GEN x, struct qfr_data *S)
933 : {
934 6994478 : GEN B, C, b = gel(x,2), c = gel(x,3);
935 6994478 : rho_get_BC(&B, &C, b, c, S);
936 6994478 : return mkvec3(c,B,C);
937 : }
938 : /* Not stack-clean */
939 : GEN
940 10296041 : qfr5_rho(GEN x, struct qfr_data *S)
941 : {
942 10296041 : GEN B, C, y, b = gel(x,2), c = gel(x,3);
943 10296041 : long sb = signe(b);
944 :
945 10296041 : rho_get_BC(&B, &C, b, c, S);
946 10296041 : y = mkvec5(c,B,C, gel(x,4), gel(x,5));
947 10296041 : if (sb) {
948 10292324 : GEN t = subii(sqri(b), S->D);
949 10292324 : if (sb < 0)
950 2299122 : t = divir(t, sqrr(subir(b,S->sqrtD)));
951 : else
952 7993202 : t = divri(sqrr(addir(b,S->sqrtD)), t);
953 : /* t = (b + sqrt(D)) / (b - sqrt(D)), evaluated stably */
954 10292324 : gel(y,5) = mulrr(t, gel(y,5)); fix_expo(y);
955 : }
956 10296041 : return y;
957 : }
958 :
959 : /* Not stack-clean */
960 : GEN
961 217084 : qfr_to_qfr5(GEN x, long prec)
962 217084 : { return mkvec5(gel(x,1),gel(x,2),gel(x,3),gen_0,real_1(prec)); }
963 :
964 : /* d0 = initial distance, x = [a,b,c, expo(d), d], d = exp(2*distance) */
965 : GEN
966 462 : qfr5_to_qfr(GEN x, GEN D, GEN d0)
967 : {
968 462 : if (d0)
969 : {
970 140 : GEN n = gel(x,4), d = absr(gel(x,5));
971 140 : if (signe(n))
972 : {
973 0 : n = addis(shifti(n, EMAX), expo(d));
974 0 : setexpo(d, 0); d = logr_abs(d);
975 0 : if (signe(n)) d = addrr(d, mulir(n, mplog2(lg(d0))));
976 0 : shiftr_inplace(d, -1);
977 0 : d0 = addrr(d0, d);
978 : }
979 140 : else if (!gequal1(d)) /* avoid loss of precision */
980 : {
981 91 : d = logr_abs(d);
982 91 : shiftr_inplace(d, -1);
983 91 : d0 = addrr(d0, d);
984 : }
985 : }
986 462 : x = qfr3_to_qfr(x, D);
987 462 : return d0 ? mkvec2(x,d0): x;
988 : }
989 :
990 : /* Not stack-clean */
991 : GEN
992 31913 : qfr3_to_qfr(GEN x, GEN d) { retmkqfb(gel(x,1), gel(x,2), gel(x,3), d); }
993 :
994 : static int
995 20824050 : ab_isreduced(GEN a, GEN b, GEN isqrtD)
996 : {
997 20824050 : if (signe(b) > 0 && abscmpii(b, isqrtD) <= 0)
998 : {
999 5505704 : GEN t = addii_sign(isqrtD,1, shifti(a,1),-1);
1000 5505704 : long l = abscmpii(b, t); /* compare |b| and |floor(sqrt(D)) - |2a|| */
1001 5505704 : if (l > 0 || (l == 0 && signe(t) < 0)) return 1;
1002 : }
1003 16727053 : return 0;
1004 : }
1005 :
1006 : INLINE int
1007 17275795 : qfr_isreduced(GEN x, GEN isqrtD)
1008 : {
1009 17275795 : return ab_isreduced(gel(x,1),gel(x,2),isqrtD);
1010 : }
1011 :
1012 : /* Not stack-clean */
1013 : GEN
1014 1947428 : qfr5_red(GEN x, struct qfr_data *S)
1015 : {
1016 1947428 : pari_sp av = avma;
1017 10247874 : while (!qfr_isreduced(x, S->isqrtD))
1018 : {
1019 8300446 : x = qfr5_rho(x, S);
1020 8300446 : if (gc_needed(av,2))
1021 : {
1022 0 : if (DEBUGMEM>1) pari_warn(warnmem,"qfr5_red");
1023 0 : x = gerepilecopy(av, x);
1024 : }
1025 : }
1026 1947428 : return x;
1027 : }
1028 : /* Not stack-clean */
1029 : GEN
1030 1169914 : qfr3_red(GEN x, struct qfr_data *S)
1031 : {
1032 1169914 : pari_sp av = avma;
1033 7027921 : while (!qfr_isreduced(x, S->isqrtD))
1034 : {
1035 5858007 : x = qfr3_rho(x, S);
1036 5858007 : if (gc_needed(av,2))
1037 : {
1038 0 : if (DEBUGMEM>1) pari_warn(warnmem,"qfr3_red");
1039 0 : x = gerepilecopy(av, x);
1040 : }
1041 : }
1042 1169914 : return x;
1043 : }
1044 :
1045 : void
1046 2170 : qfr_data_init(GEN D, long prec, struct qfr_data *S)
1047 : {
1048 2170 : S->D = D;
1049 2170 : S->sqrtD = sqrtr(itor(S->D,prec));
1050 2170 : S->isqrtD = truncr(S->sqrtD);
1051 2170 : }
1052 :
1053 : static GEN
1054 140 : qfr5_init(GEN x, GEN d, struct qfr_data *S)
1055 : {
1056 140 : long prec = realprec(d), l = -expo(d);
1057 140 : if (l < BITS_IN_LONG) l = BITS_IN_LONG;
1058 140 : prec = maxss(prec, nbits2prec(l));
1059 140 : S->D = qfb_disc(x);
1060 140 : x = qfr_to_qfr5(x,prec);
1061 140 : if (!S->sqrtD) S->sqrtD = sqrtr(itor(S->D,prec));
1062 0 : else if (typ(S->sqrtD) != t_REAL) pari_err_TYPE("qfr_init",S->sqrtD);
1063 :
1064 140 : if (!S->isqrtD)
1065 : {
1066 126 : pari_sp av=avma;
1067 : long e;
1068 126 : S->isqrtD = gcvtoi(S->sqrtD,&e);
1069 126 : if (e>-2) { set_avma(av); S->isqrtD = sqrti(S->D); }
1070 : }
1071 14 : else if (typ(S->isqrtD) != t_INT) pari_err_TYPE("qfr_init",S->isqrtD);
1072 140 : return x;
1073 : }
1074 : static GEN
1075 364 : qfr3_init(GEN x, struct qfr_data *S)
1076 : {
1077 364 : S->D = qfb_disc(x);
1078 364 : if (!S->isqrtD) S->isqrtD = sqrti(S->D);
1079 280 : else if (typ(S->isqrtD) != t_INT) pari_err_TYPE("qfr_init",S->isqrtD);
1080 364 : return x;
1081 : }
1082 :
1083 : #define qf_NOD 2
1084 : #define qf_STEP 1
1085 :
1086 : static GEN
1087 399 : redreal_i(GEN x, long flag, GEN isqrtD, GEN sqrtD)
1088 : {
1089 : struct qfr_data S;
1090 399 : GEN d = NULL, y;
1091 399 : if (typ(x)==t_VEC) { d = gel(x,2); x = gel(x,1); } else flag |= qf_NOD;
1092 399 : S.sqrtD = sqrtD;
1093 399 : S.isqrtD = isqrtD;
1094 399 : y = (flag & qf_NOD)? qfr3_init(x, &S): qfr5_init(x, d, &S);
1095 399 : switch(flag) {
1096 56 : case 0: y = qfr5_red(y,&S); break;
1097 315 : case qf_NOD: y = qfr3_red(y,&S); break;
1098 21 : case qf_STEP: y = qfr5_rho(y,&S); break;
1099 7 : case qf_STEP|qf_NOD: y = qfr3_rho(y,&S); break;
1100 0 : default: pari_err_FLAG("qfbred");
1101 : }
1102 399 : return qfr5_to_qfr(y, qfb_disc(x), d);
1103 : }
1104 : static GEN
1105 42 : redreal(GEN x) { return redreal_i(x,0,NULL,NULL); }
1106 :
1107 : GEN
1108 54774 : qfbred0(GEN x, long flag, GEN isqrtD, GEN sqrtD)
1109 : {
1110 : pari_sp av;
1111 54774 : GEN q = check_qfbext("qfbred",x);
1112 54774 : if (qfb_is_qfi(q)) return (flag & qf_STEP)? rhoimag(x): redimag(x);
1113 357 : if (typ(x)==t_QFB) flag |= qf_NOD;
1114 49 : else flag &= ~qf_NOD;
1115 357 : av = avma;
1116 357 : return gerepilecopy(av, redreal_i(x,flag,isqrtD,sqrtD));
1117 : }
1118 : /* t_QFB */
1119 : GEN
1120 17881636 : qfbred_i(GEN x) { return qfb_is_qfi(x)? redimag(x): redreal(x); }
1121 : GEN
1122 52982 : qfbred(GEN x) { return qfbred0(x, 0, NULL, NULL); }
1123 :
1124 : /* Not stack-clean */
1125 : GEN
1126 1730414 : qfr5_compraw(GEN x, GEN y)
1127 : {
1128 1730414 : GEN z = cgetg(6,t_VEC); qfb_comp(z,x,y);
1129 1730414 : if (x == y)
1130 : {
1131 34104 : gel(z,4) = shifti(gel(x,4),1);
1132 34104 : gel(z,5) = sqrr(gel(x,5));
1133 : }
1134 : else
1135 : {
1136 1696310 : gel(z,4) = addii(gel(x,4),gel(y,4));
1137 1696310 : gel(z,5) = mulrr(gel(x,5),gel(y,5));
1138 : }
1139 1730414 : fix_expo(z); return z;
1140 : }
1141 : GEN
1142 1730400 : qfr5_comp(GEN x, GEN y, struct qfr_data *S)
1143 1730400 : { return qfr5_red(qfr5_compraw(x, y), S); }
1144 : /* Not stack-clean */
1145 : GEN
1146 1006919 : qfr3_compraw(GEN x, GEN y)
1147 : {
1148 1006919 : GEN z = cgetg(4,t_VEC); qfb_comp(z,x,y);
1149 1006919 : return z;
1150 : }
1151 : GEN
1152 1006919 : qfr3_comp(GEN x, GEN y, struct qfr_data *S)
1153 1006919 : { return qfr3_red(qfr3_compraw(x,y), S); }
1154 :
1155 : /* m > 0. Not stack-clean */
1156 : static GEN
1157 7 : qfr5_powraw(GEN x, long m)
1158 : {
1159 7 : GEN y = NULL;
1160 14 : for (; m; m >>= 1)
1161 : {
1162 14 : if (m&1) y = y? qfr5_compraw(y,x): x;
1163 14 : if (m == 1) break;
1164 7 : x = qfr5_compraw(x,x);
1165 : }
1166 7 : return y;
1167 : }
1168 :
1169 : /* return x^n. Not stack-clean */
1170 : GEN
1171 28 : qfr5_pow(GEN x, GEN n, struct qfr_data *S)
1172 : {
1173 28 : GEN y = NULL;
1174 28 : long i, m, s = signe(n);
1175 28 : if (!s) return qfr5_1(S, lg(gel(x,5)));
1176 28 : if (s < 0) x = qfb_inv(x);
1177 56 : for (i=lgefint(n)-1; i>1; i--)
1178 : {
1179 28 : m = n[i];
1180 70 : for (; m; m>>=1)
1181 : {
1182 70 : if (m&1) y = y? qfr5_comp(y,x,S): x;
1183 70 : if (m == 1 && i == 2) break;
1184 42 : x = qfr5_comp(x,x,S);
1185 : }
1186 : }
1187 28 : return y;
1188 : }
1189 : /* m > 0; return x^m. Not stack-clean */
1190 : static GEN
1191 0 : qfr3_powraw(GEN x, long m)
1192 : {
1193 0 : GEN y = NULL;
1194 0 : for (; m; m>>=1)
1195 : {
1196 0 : if (m&1) y = y? qfr3_compraw(y,x): x;
1197 0 : if (m == 1) break;
1198 0 : x = qfr3_compraw(x,x);
1199 : }
1200 0 : return y;
1201 : }
1202 : /* return x^n. Not stack-clean */
1203 : GEN
1204 4529 : qfr3_pow(GEN x, GEN n, struct qfr_data *S)
1205 : {
1206 4529 : GEN y = NULL;
1207 4529 : long i, m, s = signe(n);
1208 4529 : if (!s) return qfr3_1(S);
1209 4529 : if (s < 0) x = qfb_inv(x);
1210 9074 : for (i=lgefint(n)-1; i>1; i--)
1211 : {
1212 4545 : m = n[i];
1213 5137 : for (; m; m>>=1)
1214 : {
1215 5117 : if (m&1) y = y? qfr3_comp(y,x,S): x;
1216 5117 : if (m == 1 && i == 2) break;
1217 592 : x = qfr3_comp(x,x,S);
1218 : }
1219 : }
1220 4529 : return y;
1221 : }
1222 :
1223 : static GEN
1224 7 : qfrinvraw(GEN x)
1225 : {
1226 7 : if (typ(x) == t_VEC) retmkvec2(qfbinv(gel(x,1)), negr(gel(x,2)));
1227 7 : return qfbinv(x);
1228 : }
1229 : static GEN
1230 14 : qfrpowraw(GEN x, long n)
1231 : {
1232 14 : struct qfr_data S = { NULL, NULL, NULL };
1233 14 : pari_sp av = avma;
1234 14 : if (n==1) return gcopy(x);
1235 14 : if (n==-1) return qfrinvraw(x);
1236 7 : if (typ(x)==t_QFB)
1237 : {
1238 0 : GEN D = qfb_disc(x);
1239 0 : if (!n) return qfr_1(x);
1240 0 : if (n < 0) { x = qfb_inv(x); n = -n; }
1241 0 : x = qfr3_powraw(x, n);
1242 0 : x = qfr3_to_qfr(x, D);
1243 : }
1244 : else
1245 : {
1246 7 : GEN d0 = gel(x,2);
1247 7 : x = gel(x,1);
1248 7 : if (!n) retmkvec2(qfr_1(x), real_0(precision(d0)));
1249 7 : if (n < 0) { x = qfb_inv(x); n = -n; }
1250 7 : x = qfr5_init(x, d0, &S);
1251 7 : if (labs(n) != 1) x = qfr5_powraw(x, n);
1252 7 : x = qfr5_to_qfr(x, S.D, mulrs(d0,n));
1253 : }
1254 7 : return gerepilecopy(av, x);
1255 : }
1256 : static GEN
1257 133 : qfrpow(GEN x, GEN n)
1258 : {
1259 133 : struct qfr_data S = { NULL, NULL, NULL };
1260 133 : long s = signe(n);
1261 133 : pari_sp av = avma;
1262 133 : if (typ(x)==t_QFB)
1263 : {
1264 56 : if (!s) return qfr_1(x);
1265 42 : if (s < 0) x = qfb_inv(x);
1266 42 : x = qfr3_init(x, &S);
1267 42 : x = is_pm1(n)? qfr3_red(x, &S): qfr3_pow(x, n, &S);
1268 42 : x = qfr3_to_qfr(x, S.D);
1269 : }
1270 : else
1271 : {
1272 77 : GEN d0 = gel(x,2);
1273 77 : x = gel(x,1);
1274 77 : if (!s) retmkvec2(qfr_1(x), real_0(precision(d0)));
1275 56 : if (s < 0) x = qfb_inv(x);
1276 56 : x = qfr5_init(x, d0, &S);
1277 56 : x = is_pm1(n)? qfr5_red(x, &S): qfr5_pow(x, n, &S);
1278 56 : x = qfr5_to_qfr(x, S.D, mulri(d0,n));
1279 : }
1280 98 : return gerepilecopy(av, x);
1281 : }
1282 : GEN
1283 21 : qfbpowraw(GEN x, long n)
1284 : {
1285 21 : GEN q = check_qfbext("qfbpowraw",x);
1286 21 : return qfb_is_qfi(q)? qfipowraw(x,n): qfrpowraw(x,n);
1287 : }
1288 : /* same discriminant, no distance, no checks */
1289 : GEN
1290 10082572 : qfbpow_i(GEN x, GEN n) { return qfb_is_qfi(x)? qfipow(x,n): qfrpow(x,n); }
1291 : GEN
1292 1393386 : qfbpow(GEN x, GEN n)
1293 : {
1294 1393386 : GEN q = check_qfbext("qfbpow",x);
1295 1393386 : return qfb_is_qfi(q)? qfipow(x,n): qfrpow(x,n);
1296 : }
1297 : GEN
1298 1293837 : qfbpows(GEN x, long n)
1299 : {
1300 1293837 : long N[] = { evaltyp(t_INT) | _evallg(3), 0, 0};
1301 1293837 : affsi(n, N); return qfbpow(x, N);
1302 : }
1303 :
1304 : /* Prime forms attached to prime ideals of degree 1 */
1305 :
1306 : /* assume x != 0 a t_INT, p > 0
1307 : * Return a t_QFB, but discriminant sign is not checked: can be used for
1308 : * real forms as well */
1309 : GEN
1310 12208691 : primeform_u(GEN x, ulong p)
1311 : {
1312 12208691 : GEN c, y = cgetg(5, t_QFB);
1313 12208646 : pari_sp av = avma;
1314 : ulong b;
1315 : long s;
1316 :
1317 12208646 : s = mod8(x); if (signe(x) < 0 && s) s = 8-s;
1318 : /* 2 or 3 mod 4 */
1319 12208447 : if (s & 2) pari_err_DOMAIN("primeform", "disc % 4", ">",gen_1, x);
1320 12208645 : if (p == 2) {
1321 3880181 : switch(s) {
1322 589232 : case 0: b = 0; break;
1323 2940071 : case 1: b = 1; break;
1324 350878 : case 4: b = 2; break;
1325 0 : default: pari_err_SQRTN("primeform", mkintmod(x,utoi(p)) );
1326 0 : b = 0; /* -Wall */
1327 : }
1328 3880181 : c = shifti(subsi(s,x), -3);
1329 : } else {
1330 8328464 : b = Fl_sqrt(umodiu(x,p), p);
1331 8330160 : if (b == ~0UL) pari_err_SQRTN("primeform", mkintmod(x,utoi(p)) );
1332 : /* mod(b) != mod2(x) ? */
1333 8330614 : if ((b ^ s) & 1) b = p - b;
1334 8330614 : c = diviuexact(shifti(subii(sqru(b), x), -2), p);
1335 : }
1336 12195041 : gel(y,3) = gerepileuptoint(av, c);
1337 12208092 : gel(y,4) = icopy(x);
1338 12207952 : gel(y,2) = utoi(b);
1339 12208022 : gel(y,1) = utoipos(p); return y;
1340 : }
1341 :
1342 : /* special case: p = 1 return unit form */
1343 : GEN
1344 135459 : primeform(GEN x, GEN p)
1345 : {
1346 135459 : const char *f = "primeform";
1347 : pari_sp av;
1348 135459 : long s, sx = signe(x), sp = signe(p);
1349 : GEN y, b, absp;
1350 :
1351 135459 : if (typ(x) != t_INT) pari_err_TYPE(f,x);
1352 135459 : if (typ(p) != t_INT) pari_err_TYPE(f,p);
1353 135459 : if (!sp) pari_err_DOMAIN(f,"p","=",gen_0,p);
1354 135459 : if (!sx) pari_err_DOMAIN(f,"D","=",gen_0,x);
1355 135459 : if (lgefint(p) == 3)
1356 : {
1357 135445 : ulong pp = p[2];
1358 135445 : if (pp == 1) {
1359 17782 : if (sx < 0) {
1360 : long r;
1361 10950 : if (sp < 0) pari_err_IMPL("negative definite t_QFB");
1362 10950 : r = mod4(x);
1363 10950 : if (r && r != 3) pari_err_DOMAIN(f,"disc % 4",">", gen_1,x);
1364 10950 : return qfi_1_by_disc(x);
1365 : }
1366 6832 : y = qfr_1_by_disc(x);
1367 6832 : if (sp < 0) { gel(y,1) = gen_m1; togglesign(gel(y,3)); }
1368 6832 : return y;
1369 : }
1370 117663 : y = primeform_u(x, pp);
1371 117656 : if (sx < 0) {
1372 89957 : if (sp < 0) pari_err_IMPL("negative definite t_QFB");
1373 89957 : return y;
1374 : }
1375 27699 : if (sp < 0) { togglesign(gel(y,1)); togglesign(gel(y,3)); }
1376 27699 : return gcopy( qfr3_to_qfr(y, x) );
1377 : }
1378 14 : s = mod8(x);
1379 14 : if (sx < 0)
1380 : {
1381 7 : if (sp < 0) pari_err_IMPL("negative definite t_QFB");
1382 7 : if (s) s = 8-s;
1383 : }
1384 14 : y = cgetg(5, t_QFB);
1385 : /* 2 or 3 mod 4 */
1386 14 : if (s & 2) pari_err_DOMAIN(f, "disc % 4", ">",gen_1, x);
1387 14 : absp = absi_shallow(p); av = avma;
1388 14 : b = Fp_sqrt(x, absp); if (!b) pari_err_SQRTN(f, mkintmod(x,absp));
1389 14 : s &= 1; /* s = x mod 2 */
1390 : /* mod(b) != mod2(x) ? [Warning: we may have b == 0] */
1391 14 : if ((!signe(b) && s) || mod2(b) != s) b = gerepileuptoint(av, subii(absp,b));
1392 :
1393 14 : av = avma;
1394 14 : gel(y,3) = gerepileuptoint(av, diviiexact(shifti(subii(sqri(b), x), -2), p));
1395 14 : gel(y,4) = icopy(x);
1396 14 : gel(y,2) = b;
1397 14 : gel(y,1) = icopy(p);
1398 14 : return y;
1399 : }
1400 :
1401 : static GEN
1402 2620771 : normforms(GEN D, GEN fa)
1403 : {
1404 : long i, j, k, lB, aN, sa;
1405 : GEN a, L, V, B, N, N2;
1406 2620771 : int D_odd = mpodd(D);
1407 2620771 : a = typ(fa) == t_INT ? fa: typ(fa) == t_VEC? gel(fa,1): factorback(fa);
1408 2620771 : sa = signe(a);
1409 2620771 : if (sa==0 || (signe(D)<0 && sa<0)) return NULL;
1410 1203972 : V = D_odd? Zn_quad_roots(fa, gen_1, shifti(subsi(1, D), -2))
1411 2551765 : : Zn_quad_roots(fa, gen_0, negi(shifti(D, -2)));
1412 2551766 : if (!V) return NULL;
1413 511966 : N = gel(V,1); B = gel(V,2); lB = lg(B);
1414 511966 : N2 = shifti(N,1);
1415 511966 : aN = itou(diviiexact(a, N)); /* |a|/N */
1416 511965 : L = cgetg((lB-1)*aN+1, t_VEC);
1417 2360563 : for (k = 1, i = 1; i < lB; i++)
1418 : {
1419 1848598 : GEN b = shifti(gel(B,i), 1), c, C;
1420 1848599 : if (D_odd) b = addiu(b , 1);
1421 1848599 : c = diviiexact(shifti(subii(sqri(b), D), -2), a);
1422 1848595 : for (j = 0;; b = addii(b, N2))
1423 : {
1424 2216669 : gel(L, k++) = mkqfb(a, b, c, D);
1425 2216672 : if (++j == aN) break;
1426 368074 : C = addii(b, N); if (aN > 1) C = diviuexact(C, aN);
1427 368074 : c = sa > 0? addii(c, C): subii(c, C);
1428 : }
1429 : }
1430 511965 : return L;
1431 : }
1432 :
1433 : /* Let M and N in SL2(Z), return (N*M^-1)[,1] */
1434 : static GEN
1435 344321 : SL2_div_mul_e1(GEN N, GEN M)
1436 : {
1437 344321 : GEN b = gcoeff(M,2,1), d = gcoeff(M,2,2);
1438 344321 : GEN A = mulii(gcoeff(N,1,1), d), B = mulii(gcoeff(N,1,2), b);
1439 344315 : GEN C = mulii(gcoeff(N,2,1), d), D = mulii(gcoeff(N,2,2), b);
1440 344318 : retmkvec2(subii(A,B), subii(C,D));
1441 : }
1442 : static GEN
1443 1445675 : qfisolve_normform(GEN Q, GEN P)
1444 : {
1445 1445675 : GEN a = gel(Q,1), N = gel(Q,2);
1446 1445675 : GEN M, b = redimagsl2(P, &M);
1447 1445681 : if (!qfb_equal(a,b)) return NULL;
1448 102128 : return SL2_div_mul_e1(N,M);
1449 : }
1450 :
1451 : /* Test equality modulo GL2 of two reduced forms */
1452 : static int
1453 61068 : GL2_qfb_equal(GEN a, GEN b)
1454 : {
1455 61068 : return equalii(gel(a,1),gel(b,1))
1456 11361 : && absequalii(gel(a,2),gel(b,2))
1457 72429 : && equalii(gel(a,3),gel(b,3));
1458 : }
1459 :
1460 : /* Q(u,v) = p; if s < 0 return that solution; else the set of all solutions */
1461 : static GEN
1462 48018 : allsols(GEN Q, long s, GEN u, GEN v)
1463 : {
1464 48018 : GEN w = mkvec2(u, v), b;
1465 48010 : if (signe(v) < 0) { u = negi(u); v = negi(v); } /* normalize for v >= 0 */
1466 48010 : w = mkvec2(u, v); if (s < 0) return w;
1467 41370 : if (!s) return mkvec(w);
1468 38955 : b = gel(Q,2); /* sum of the 2 solutions (if they exist) is -bv / a */
1469 38955 : if (signe(b))
1470 : { /* something to check */
1471 : GEN r, t;
1472 13433 : t = dvmdii(mulii(b, v), gel(Q,1), &r);
1473 13433 : if (signe(r)) return mkvec(w);
1474 1820 : u = addii(u, t);
1475 : }
1476 27342 : return mkvec2(w, mkvec2(negi(u), v));
1477 : }
1478 : static GEN
1479 223051 : qfisolvep_all(GEN Q, GEN p, long all)
1480 : {
1481 223051 : GEN R, U, V, M, N, x, q, D = qfb_disc(Q);
1482 223050 : long s = kronecker(D, p);
1483 :
1484 223046 : if (s < 0) return NULL;
1485 126971 : if (!all) s = -1; /* to indicate we want a single solution */
1486 : /* Solutions iff a class of maximal ideal above p is the class of Q;
1487 : * Two solutions iff (s > 0 and the class has order > 2), else one */
1488 126971 : if (!signe(gel(Q,2)))
1489 : { /* if principal form, use faster cornacchia */
1490 43650 : GEN a = gel(Q,1), c = gel(Q,3);
1491 43650 : if (equali1(a))
1492 : {
1493 38174 : if (!cornacchia(c, p, &M,&N)) return NULL;
1494 33701 : return allsols(Q, s, M, N);
1495 : }
1496 5474 : if (equali1(c))
1497 : {
1498 5194 : if (!cornacchia(a, p, &M,&N)) return NULL;
1499 721 : return allsols(Q, s, N, M);
1500 : }
1501 : }
1502 83601 : R = redimagsl2(Q, &U);
1503 83601 : if (equali1(gel(R,1)))
1504 : { /* principal form */
1505 22533 : if (!signe(gel(R,2)))
1506 : {
1507 4396 : if (!cornacchia(gel(R,3), p, &M,&N)) return NULL;
1508 812 : x = mkvec2(M,N);
1509 : }
1510 : else
1511 : { /* x^2 + xy + ((1-D)/4)y^2 = p <==> (2x + y)^2 - D y^2 = 4p */
1512 18137 : if (!cornacchia2(negi(D), p, &M, &N)) return NULL;
1513 2331 : x = subii(M,N); if (mpodd(x)) return NULL;
1514 2331 : x = mkvec2(shifti(x,-1), N);
1515 : }
1516 3143 : x = ZM_ZC_mul(U, x); x[0] = evaltyp(t_VEC) | _evallg(3); /* transpose */
1517 3143 : return allsols(Q, s, gel(x,1), gel(x,2));
1518 : }
1519 61068 : q = redimagsl2(primeform(D, p), &V);
1520 61068 : if (!GL2_qfb_equal(R,q)) return NULL;
1521 10451 : if (signe(gel(R,2)) != signe(gel(q,2))) gcoeff(V,2,1) = negi(gcoeff(V,2,1));
1522 10451 : x = SL2_div_mul_e1(U,V); return allsols(Q, s, gel(x,1), gel(x,2));
1523 : }
1524 : GEN
1525 0 : qfisolvep(GEN Q, GEN p)
1526 : {
1527 0 : pari_sp av = avma;
1528 0 : GEN x = qfisolvep_all(Q, p, 0);
1529 0 : return x? gerepilecopy(av, x): gc_const(av, gen_0);
1530 : }
1531 :
1532 : static void
1533 13379274 : _rhorealsl2(GEN *pa, GEN *pb, GEN *pc, GEN *pu1, GEN *pu2, GEN *pv1,
1534 : GEN *pv2, GEN d, GEN rd)
1535 : {
1536 13379274 : GEN C = mpabs_shallow(*pc), t = addii(*pb, gmax_shallow(rd,C));
1537 13379274 : GEN r, q = truedvmdii(t, shifti(C,1), &r);
1538 13379274 : *pb = subii(t, addii(r, *pb));
1539 13379274 : *pa = *pc;
1540 13379274 : *pc = diviiexact(subii(sqri(*pb), d), shifti(*pa, 2));
1541 13379274 : if (signe(*pa) < 0) togglesign(q);
1542 13379274 : r = *pu1; *pu1 = *pv1; *pv1 = subii(mulii(q, *pv1), r);
1543 13379274 : r = *pu2; *pu2 = *pv2; *pv2 = subii(mulii(q, *pv2), r);
1544 13379274 : }
1545 :
1546 : static GEN
1547 10810674 : rhorealsl2(GEN A, GEN rd)
1548 : {
1549 10810674 : GEN V = gel(A,1), M = gel(A,2);
1550 10810674 : GEN a = gel(V,1), b = gel(V,2), c = gel(V,3), d = qfb_disc(V);
1551 10810674 : GEN u1 = gcoeff(M,1,1), v1 = gcoeff(M,1,2);
1552 10810674 : GEN u2 = gcoeff(M,2,1), v2 = gcoeff(M,2,2);
1553 10810674 : _rhorealsl2(&a,&b,&c, &u1,&u2,&v1,&v2, d, rd);
1554 10810674 : return mkvec2(mkqfb(a,b,c,d), mkmat22(u1,v1,u2,v2));
1555 : }
1556 :
1557 : static GEN
1558 979655 : redrealsl2(GEN V, GEN rd)
1559 : {
1560 979655 : pari_sp av = avma;
1561 979655 : GEN u1 = gen_1, u2 = gen_0, v1 = gen_0, v2 = gen_1;
1562 979655 : GEN a = gel(V,1), b = gel(V,2), c = gel(V,3), d = qfb_disc(V);
1563 3548255 : while (!ab_isreduced(a,b,rd))
1564 : {
1565 2568600 : _rhorealsl2(&a,&b,&c, &u1,&u2,&v1,&v2, d, rd);
1566 2568600 : if (gc_needed(av, 1))
1567 : {
1568 0 : if (DEBUGMEM>1) pari_warn(warnmem,"redrealsl2");
1569 0 : gerepileall(av, 7, &a,&b,&c,&u1,&u2,&v1,&v2);
1570 : }
1571 : }
1572 979655 : return gerepilecopy(av, mkvec2(mkqfb(a,b,c,d), mkmat22(u1,v1,u2,v2)));
1573 : }
1574 :
1575 : GEN
1576 519693 : qfbredsl2(GEN q, GEN isD)
1577 : {
1578 : pari_sp av;
1579 519693 : if (typ(q) != t_QFB) pari_err_TYPE("qfbredsl2",q);
1580 519693 : if (qfb_is_qfi(q))
1581 : {
1582 310157 : GEN v = cgetg(3,t_VEC);
1583 310157 : if (isD) pari_err_TYPE("qfbredsl2", isD);
1584 310157 : gel(v,1) = redimagsl2(q, &gel(v,2)); return v;
1585 : }
1586 209536 : av = avma;
1587 209536 : if (!isD) isD = sqrti(qfb_disc(q));
1588 208068 : else if (typ(isD) != t_INT) pari_err_TYPE("qfbredsl2",isD);
1589 209529 : return gerepileupto(av, redrealsl2(q, isD));
1590 : }
1591 :
1592 : static GEN
1593 770126 : qfrsolve_normform(GEN N, GEN Ps, GEN rd)
1594 : {
1595 770126 : pari_sp av = avma, btop;
1596 770126 : GEN M = N, P = redrealsl2(Ps, rd), Q = P;
1597 :
1598 770126 : btop = avma;
1599 : for(;;)
1600 : {
1601 5840681 : if (qfb_equal(gel(M,1), gel(P,1)))
1602 154084 : return gerepileupto(av, SL2_div_mul_e1(gel(M,2),gel(P,2)));
1603 5686597 : if (qfb_equal(gel(N,1), gel(Q,1)))
1604 77658 : return gerepileupto(av, SL2_div_mul_e1(gel(N,2),gel(Q,2)));
1605 5608939 : M = rhorealsl2(M, rd);
1606 5608939 : if (qfb_equal(gel(M,1), gel(N,1))) return gc_NULL(av);
1607 5201735 : Q = rhorealsl2(Q, rd);
1608 5201735 : if (qfb_equal(gel(P,1), gel(Q,1))) return gc_NULL(av);
1609 5070555 : if (gc_needed(btop, 1)) gerepileall(btop, 2, &M, &Q);
1610 : }
1611 : }
1612 :
1613 : GEN
1614 0 : qfrsolvep(GEN Q, GEN p)
1615 : {
1616 0 : pari_sp av = avma;
1617 0 : GEN N, x, rd, d = qfb_disc(Q);
1618 :
1619 0 : if (kronecker(d, p) < 0) return gc_const(av, gen_0);
1620 0 : rd = sqrti(d);
1621 0 : N = redrealsl2(Q, rd);
1622 0 : x = qfrsolve_normform(N, primeform(d, p), rd);
1623 0 : return x? gerepileupto(av, x): gc_const(av, gen_0);
1624 : }
1625 :
1626 : static GEN
1627 1862943 : known_prime(GEN v)
1628 : {
1629 1862943 : GEN p, e, fa = check_arith_all(v, "qfbsolve");
1630 1862940 : if (!fa) return BPSW_psp(v)? v: NULL;
1631 42082 : if (lg(gel(fa,1)) != 2) return NULL;
1632 29357 : p = gcoeff(fa,1,1);
1633 29357 : e = gcoeff(fa,1,2);
1634 29357 : return (equali1(e) && !is_pm1(p) && signe(p) > 0)? p: NULL;
1635 : }
1636 : static GEN
1637 2215801 : qfsolve_normform(GEN Q, GEN f, GEN rd)
1638 2215801 : { return rd? qfrsolve_normform(Q, f, rd): qfisolve_normform(Q, f); }
1639 : static GEN
1640 2843814 : qfbsolve_primitive_i(GEN Q, GEN rd, GEN *Qr, GEN fa, long all)
1641 : {
1642 : GEN x, W, F, p;
1643 : long i, j, l;
1644 2843814 : if (!rd && (p = known_prime(fa))) return qfisolvep_all(Q, p, all);
1645 2620760 : F = normforms(qfb_disc(Q), fa);
1646 2620771 : if (!F) return NULL;
1647 511965 : if (!*Qr) *Qr = qfbredsl2(Q, rd);
1648 511966 : l = lg(F); W = all? cgetg(l, t_VEC): NULL;
1649 2727257 : for (j = i = 1; i < l; i++)
1650 2215801 : if ((x = qfsolve_normform(*Qr, gel(F,i), rd)))
1651 : {
1652 333867 : if (!all) return x;
1653 333356 : gel(W,j++) = x;
1654 : }
1655 511456 : if (j == 1) return NULL;
1656 127457 : setlg(W,j); return W;
1657 : }
1658 :
1659 : static GEN
1660 2838508 : qfb_initrd(GEN Q) { GEN d = qfb_disc(Q); return signe(d) > 0? sqrti(d): NULL; }
1661 : static GEN
1662 2828281 : qfbsolve_primitive(GEN Q, GEN fa, long all)
1663 : {
1664 2828281 : GEN x, Qr = NULL, rdQ = qfb_initrd(Q);
1665 2828281 : x = qfbsolve_primitive_i(Q, rdQ, &Qr, fa, all);
1666 2828295 : if (!x) return cgetg(1, t_VEC);
1667 174756 : return x;
1668 : }
1669 :
1670 : /* f / g^2 */
1671 : static GEN
1672 5299 : famat_divsqr(GEN f, GEN g)
1673 5299 : { return famat_reduce(famat_div_shallow(f, famat_pows_shallow(g,2))); }
1674 : static GEN
1675 10227 : qfbsolve_all(GEN Q, GEN n, long all)
1676 : {
1677 10227 : GEN W, Qr = NULL, fa = factorint(n, 0), rdQ = qfb_initrd(Q);
1678 10227 : GEN D = divisors_factored(mkmat2(gel(fa,1), gshift(gel(fa,2),-1)));
1679 10227 : long i, j, l = lg(D);
1680 10227 : W = all? cgetg(l, t_VEC): NULL;
1681 25151 : for (i = j = 1; i < l; i++)
1682 : {
1683 15526 : GEN w, d = gel(D,i), FA = i == 1? fa: famat_divsqr(fa, gel(d,2));
1684 15526 : if ((w = qfbsolve_primitive_i(Q, rdQ, &Qr, FA, all)))
1685 : {
1686 1218 : if (i != 1) w = RgV_Rg_mul(w, gel(d,1));
1687 1218 : if (!all) return w;
1688 616 : gel(W,j++) = w;
1689 : }
1690 : }
1691 9625 : if (j == 1) return cgetg(1, t_VEC);
1692 525 : setlg(W,j); return shallowconcat1(W);
1693 : }
1694 :
1695 : GEN
1696 2838516 : qfbsolve(GEN Q, GEN n, long fl)
1697 : {
1698 2838516 : pari_sp av = avma;
1699 2838516 : if (typ(Q) != t_QFB) pari_err_TYPE("qfbsolve",Q);
1700 2838516 : if (fl < 0 || fl > 3) pari_err_FLAG("qfbsolve");
1701 5666801 : return gerepilecopy(av, (fl & 2)? qfbsolve_all(Q, n, fl & 1)
1702 2828281 : : qfbsolve_primitive(Q, n, fl & 1));
1703 : }
1704 :
1705 : /* 1 if there exists x,y such that x^2 + dy^2 = p, 0 otherwise;
1706 : * Assume d > 0 and p is prime */
1707 : long
1708 55244 : cornacchia(GEN d, GEN p, GEN *px, GEN *py)
1709 : {
1710 55244 : pari_sp av = avma;
1711 : GEN b, c, r;
1712 :
1713 55244 : *px = *py = gen_0;
1714 55244 : b = subii(p, d);
1715 55176 : if (signe(b) < 0) return gc_long(av,0);
1716 54966 : if (signe(b) == 0) { *py = gen_1; return gc_long(av,1); }
1717 54959 : b = Fp_sqrt(b, p); /* sqrt(-d) */
1718 55026 : if (!b) return gc_long(av,0);
1719 51295 : b = gmael(halfgcdii(p, b), 2, 2);
1720 51343 : c = dvmdii(subii(p, sqri(b)), d, &r);
1721 51216 : if (r != gen_0 || !Z_issquareall(c, &c)) return gc_long(av,0);
1722 35488 : set_avma(av);
1723 35482 : *px = icopy(b);
1724 35474 : *py = icopy(c); return 1;
1725 : }
1726 :
1727 : static GEN
1728 1128545 : lastqi(GEN Q)
1729 : {
1730 1128545 : GEN s = gcoeff(Q,1,1), q = gcoeff(Q,1,2), p = absi_shallow(gcoeff(Q,2,2));
1731 1128541 : if (!signe(q)) return gen_0;
1732 1128352 : if (!signe(s)) return p;
1733 1122528 : if (is_pm1(q)) return subiu(p,1);
1734 1122527 : return divii(p, absi_shallow(q));
1735 : }
1736 :
1737 : static long
1738 1128555 : cornacchia2_i(long av, GEN d, GEN p, GEN b, GEN px4, GEN *px, GEN *py)
1739 : {
1740 : GEN M, Q, V, c, r, b2;
1741 1128555 : if (!signe(b)) { /* d = p,2p,3p,4p */
1742 14 : set_avma(av);
1743 14 : if (absequalii(d, px4)){ *py = gen_1; return 1; }
1744 14 : if (absequalii(d, p)) { *py = gen_2; return 1; }
1745 0 : return 0;
1746 : }
1747 1128541 : if (mod2(b) != mod2(d)) b = subii(p,b);
1748 1128503 : M = halfgcdii(shifti(p,1), b); Q = gel(M,1); V = gel(M, 2);
1749 1128548 : b = addii(mulii(gel(V,1), lastqi(Q)), gel(V,2));
1750 1128431 : b2 = sqri(b);
1751 1128430 : if (cmpii(b2,px4) > 0)
1752 : {
1753 1120770 : b = gel(V,1); b2 = sqri(b);
1754 1120749 : if (cmpii(b2,px4) > 0) { b = gel(V,2); b2 = sqri(b); }
1755 : }
1756 1128476 : c = dvmdii(subii(px4, b2), d, &r);
1757 1128457 : if (r != gen_0 || !Z_issquareall(c, &c)) return gc_long(av,0);
1758 1083735 : set_avma(av);
1759 1083736 : *px = icopy(b);
1760 1083738 : *py = icopy(c); return 1;
1761 : }
1762 :
1763 : /* 1 if there exists x,y such that x^2 + dy^2 = 4p, 0 otherwise;
1764 : * Assume d > 0 is congruent to 0 or 3 mod 4 and p is prime */
1765 : long
1766 1088770 : cornacchia2(GEN d, GEN p, GEN *px, GEN *py)
1767 : {
1768 1088770 : pari_sp av = avma;
1769 1088770 : GEN b, p4 = shifti(p,2);
1770 :
1771 1088709 : *px = *py = gen_0;
1772 1088709 : if (abscmpii(p4, d) < 0) return gc_long(av,0);
1773 1087935 : if (absequaliu(p, 2))
1774 : {
1775 7 : set_avma(av);
1776 7 : switch (itou_or_0(d)) {
1777 0 : case 4: *px = gen_2; break;
1778 0 : case 7: *px = gen_1; break;
1779 7 : default: return 0;
1780 : }
1781 0 : *py = gen_1; return 1;
1782 : }
1783 1087926 : b = Fp_sqrt(negi(d), p);
1784 1087961 : if (!b) return gc_long(av,0);
1785 1087877 : return cornacchia2_i(av, d, p, b, p4, px, py);
1786 : }
1787 :
1788 : /* 1 if there exists x,y such that x^2 + dy^2 = 4p [p prime], 0 otherwise */
1789 : long
1790 40677 : cornacchia2_sqrt(GEN d, GEN p, GEN b, GEN *px, GEN *py)
1791 : {
1792 40677 : pari_sp av = avma;
1793 40677 : GEN p4 = shifti(p,2);
1794 40677 : *px = *py = gen_0;
1795 40677 : if (abscmpii(p4, d) < 0) return gc_long(av,0);
1796 40677 : return cornacchia2_i(av, d, p, b, p4, px, py);
1797 : }
1798 :
1799 : GEN
1800 7630 : qfbcornacchia(GEN d, GEN p)
1801 : {
1802 7630 : pari_sp av = avma;
1803 : GEN x, y;
1804 7630 : if (typ(d) != t_INT || signe(d) <= 0) pari_err_TYPE("qfbcornacchia", d);
1805 7630 : if (typ(p) != t_INT || cmpiu(p, 2) < 0) pari_err_TYPE("qfbcornacchia", p);
1806 7630 : if (mod4(p)? cornacchia(d, p, &x, &y): cornacchia2(d, shifti(p, -2), &x, &y))
1807 287 : return gerepilecopy(av, mkvec2(x, y));
1808 7343 : set_avma(av); return cgetg(1, t_VEC);
1809 : }
|
