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 151934 : check_quaddisc(GEN x, long *s, long *pr, const char *f)
24 : {
25 : long r;
26 151934 : if (typ(x) != t_INT) pari_err_TYPE(f,x);
27 151920 : *s = signe(x);
28 151920 : if (Z_issquare(x)) pari_err_DOMAIN(f,"issquare(disc)","=", gen_1,x);
29 151920 : r = mod4(x); if (*s < 0 && r) r = 4 - r;
30 151920 : if (r > 1) pari_err_DOMAIN(f,"disc % 4",">", gen_1,x);
31 151906 : if (pr) *pr = r;
32 151906 : }
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 2177 : check_quaddisc_imag(GEN x, long *r, const char *f)
41 : {
42 2177 : long sx; check_quaddisc(x, &sx, r, f);
43 2170 : if (sx > 0) pari_err_DOMAIN(f, "disc",">",gen_0,x);
44 2170 : }
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 1018505 : quadpoly_bc(GEN D, long Dodd, GEN *b, GEN *c)
50 : {
51 1018505 : if (Dodd)
52 : {
53 920225 : pari_sp av = avma;
54 920225 : *b = gen_m1;
55 920225 : *c = gc_INT(av, shifti(subui(1,D), -2));
56 : }
57 : else
58 : {
59 98280 : *b = gen_0;
60 98280 : *c = shifti(D,-2); togglesign(*c);
61 : }
62 1018505 : }
63 : /* X^2 - X - (D-1)/4 or X^2 - D/4 */
64 : static GEN
65 243348 : quadpoly_ii(GEN D, long Dmod4)
66 : {
67 243348 : GEN b, c, y = cgetg(5,t_POL);
68 243348 : y[1] = evalsigne(1) | evalvarn(0);
69 243348 : quadpoly_bc(D, Dmod4, &b,&c);
70 243348 : gel(y,2) = c;
71 243348 : gel(y,3) = b;
72 243348 : gel(y,4) = gen_1; return y;
73 : }
74 : GEN
75 2051 : quadpoly(GEN D)
76 : {
77 : long s, Dmod4;
78 2051 : check_quaddisc(D, &s, &Dmod4, "quadpoly");
79 2044 : return quadpoly_ii(D, Dmod4);
80 : }
81 : GEN /* no checks */
82 241304 : 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 574 : quadgen0(GEN x, long v)
98 : {
99 574 : if (v==-1) v = fetch_user_var("w");
100 574 : 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 1945664 : check_qfbext(const char *fun, GEN x)
113 : {
114 1945664 : long t = typ(x);
115 1945664 : if (t == t_QFB) return x;
116 196 : 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 139081 : qfb3(GEN x, GEN y, GEN z)
127 139081 : { 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 1592908 : && equalii(gel(x,2),gel(y,2))
134 25376536 : && equalii(gel(x,3),gel(y,3));
135 : }
136 :
137 : /* valid for t_QFB, qfr3, qfr5; shallow */
138 : static GEN
139 932593 : qfb_inv(GEN x) {
140 932593 : GEN z = shallowcopy(x);
141 932593 : gel(z,2) = negi(gel(z,2));
142 932593 : return z;
143 : }
144 : /* valid for t_QFB, GC clean */
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 77231 : Qfb0(GEN a, GEN b, GEN c)
151 : {
152 : GEN q, D;
153 77231 : if (!b)
154 : {
155 49 : if (c) pari_err_TYPE("Qfb",c);
156 42 : if (typ(a) == t_VEC && lg(a) == 4)
157 21 : { b = gel(a,2); c = gel(a,3); a = gel(a,1); }
158 21 : else if (typ(a) == t_POL && degpol(a) == 2)
159 7 : { b = gel(a,3); c = gel(a,2); a = gel(a,4); }
160 14 : else if (typ(a) == t_MAT && lg(a)==3 && lgcols(a)==3)
161 : {
162 7 : b = gadd(gcoeff(a,2,1), gcoeff(a,1,2));
163 7 : c = gcoeff(a,2,2); a = gcoeff(a,1,1);
164 : }
165 : else
166 7 : pari_err_TYPE("Qfb",a);
167 : }
168 77182 : else if (!c)
169 7 : pari_err_TYPE("Qfb",b);
170 77210 : if (typ(a)!=t_INT) pari_err_TYPE("Qfb",a);
171 77203 : if (typ(b)!=t_INT) pari_err_TYPE("Qfb",b);
172 77203 : if (typ(c)!=t_INT) pari_err_TYPE("Qfb",c);
173 77203 : q = qfb3(a, b, c); D = qfb_disc(q);
174 77203 : if (signe(D) < 0)
175 42392 : { if (signe(a) < 0) pari_err_IMPL("negative definite t_QFB"); }
176 34811 : else if (Z_issquare(D)) pari_err_DOMAIN("Qfb","issquare(disc)","=", gen_1,q);
177 77196 : return q;
178 : }
179 :
180 : /***********************************************************************/
181 : /** **/
182 : /** Reduction **/
183 : /** **/
184 : /***********************************************************************/
185 :
186 : /* assume a > 0. Write b = q*2a + r, with -a < r <= a */
187 : static GEN
188 16908786 : dvmdii_round(GEN b, GEN a, GEN *r)
189 : {
190 16908786 : GEN a2 = shifti(a, 1), q = dvmdii(b, a2, r);
191 16908424 : if (signe(b) >= 0) {
192 9290093 : if (abscmpii(*r, a) > 0) { q = addiu(q, 1); *r = subii(*r, a2); }
193 : } else { /* r <= 0 */
194 7618331 : if (abscmpii(*r, a) >= 0){ q = subiu(q, 1); *r = addii(*r, a2); }
195 : }
196 16907997 : return q;
197 : }
198 : /* Assume 0 < a <= LONG_MAX. Ensure no overflow */
199 : static long
200 114624510 : dvmdsu_round(long b, ulong a, long *r)
201 : {
202 114624510 : ulong a2 = a << 1, q, ub, ur;
203 114624510 : if (b >= 0) {
204 73223884 : ub = b;
205 73223884 : q = ub / a2;
206 73223884 : ur = ub % a2;
207 73223884 : if (ur > a) { ur -= a; q++; *r = (long)ur; *r -= (long)a; }
208 26713124 : else *r = (long)ur;
209 73223884 : return (long)q;
210 : } else { /* r <= 0 */
211 41400626 : ub = (ulong)-b; /* |b| */
212 41400626 : q = ub / a2;
213 41400626 : ur = ub % a2;
214 41400626 : if (ur >= a) { ur -= a; q++; *r = (long)ur; *r = (long)a - *r; }
215 22036737 : else *r = -(long)ur;
216 41400626 : return -(long)q;
217 : }
218 : }
219 : /* reduce b mod 2*a. Update b,c */
220 : static void
221 2783988 : REDB(GEN a, GEN *b, GEN *c)
222 : {
223 2783988 : GEN r, q = dvmdii_round(*b, a, &r);
224 2783988 : if (!signe(q)) return;
225 2714529 : *c = subii(*c, mulii(q, shifti(addii(*b, r),-1)));
226 2714529 : *b = r;
227 : }
228 : /* Assume a > 0. Reduce b mod 2*a. Update b,c */
229 : static void
230 114617513 : sREDB(ulong a, long *b, ulong *c)
231 : {
232 : long r, q;
233 : ulong uz;
234 124468677 : if (a > LONG_MAX) return; /* b already reduced */
235 114617513 : q = dvmdsu_round(*b, a, &r);
236 114639006 : if (q == 0) return;
237 : /* Final (a,r,c2) satisfies |r| <= |b| hence c2 <= c, c2 = c - q*z,
238 : * where z = (b+r) / 2, representable as long, has the same sign as q. */
239 104787842 : if (*b < 0)
240 : { /* uz = -z >= 0, q < 0 */
241 36417792 : if (r >= 0) /* different signs=>no overflow, exact division */
242 19419241 : uz = (ulong)-((*b + r)>>1);
243 : else
244 : {
245 16998551 : ulong ub = (ulong)-*b, ur = (ulong)-r;
246 16998551 : uz = (ub + ur) >> 1;
247 : }
248 36417792 : *c -= (-q) * uz; /* c -= qz */
249 : }
250 : else
251 : { /* uz = z >= 0, q > 0 */
252 68370050 : if (r <= 0)
253 46570070 : uz = (*b + r)>>1;
254 : else
255 : {
256 21799980 : ulong ub = (ulong)*b, ur = (ulong)r;
257 21799980 : uz = ((ub + ur) >> 1);
258 : }
259 68370050 : *c -= q * uz; /* c -= qz */
260 : }
261 104787842 : *b = r;
262 : }
263 : static void
264 14124796 : REDBU(GEN a, GEN *b, GEN *c, GEN u1, GEN *u2)
265 : { /* REDB(a,b,c) */
266 14124796 : GEN r, q = dvmdii_round(*b, a, &r);
267 14123994 : *c = subii(*c, mulii(q, shifti(addii(*b, r),-1)));
268 14123927 : *b = r;
269 14123927 : *u2 = subii(*u2, mulii(q, u1));
270 14123940 : }
271 :
272 : /* q t_QFB, return reduced representative and set base change U in Sl2(Z) */
273 : static GEN
274 6773308 : qfi_redsl2_basecase(GEN q, GEN *U)
275 : {
276 6773308 : pari_sp av = avma;
277 : GEN z, u1,u2,v1,v2,Q;
278 6773308 : GEN a = gel(q,1), b = gel(q,2), c = gel(q,3);
279 : long cmp;
280 6773308 : u1 = gen_1; u2 = gen_0;
281 6773308 : cmp = abscmpii(a, b);
282 6773310 : if (cmp < 0)
283 2195082 : REDBU(a,&b,&c, u1,&u2);
284 4578228 : else if (cmp == 0 && signe(b) < 0)
285 : { /* b = -a */
286 11778 : b = negi(b);
287 11778 : u2 = gen_1;
288 : }
289 : for(;;)
290 : {
291 18702362 : cmp = abscmpii(a, c); if (cmp <= 0) break;
292 11928596 : swap(a,c); b = negi(b);
293 11929674 : z = u1; u1 = u2; u2 = negi(z);
294 11929701 : REDBU(a,&b,&c, u1,&u2);
295 11929058 : if (gc_needed(av, 1)) {
296 7 : if (DEBUGMEM>1) pari_warn(warnmem, "qfbredsl2");
297 7 : (void)gc_all(av, 5, &a,&b,&c, &u1,&u2);
298 : }
299 : }
300 6773270 : if (cmp == 0 && signe(b) < 0)
301 : {
302 17677 : b = negi(b);
303 17677 : z = u1; u1 = u2; u2 = negi(z);
304 : }
305 : /* Let q = (A,B,C). q o [u1,u2; v1,v2] = Q implies
306 : * [v1,v2] = (1/C) [(b-B)/2 u1 - a u2, c u1 - (b+B)/2 u2] */
307 6773270 : z = shifti(subii(b, gel(q,2)), -1);
308 6773250 : v1 = subii(mulii(z, u1), mulii(a, u2)); v1 = diviiexact(v1, gel(q,3));
309 6773271 : z = subii(z, b);
310 6773267 : v2 = addii(mulii(z, u2), mulii(c, u1)); v2 = diviiexact(v2, gel(q,3));
311 6773259 : *U = mkmat2(mkcol2(u1,v1), mkcol2(u2,v2));
312 6773327 : Q = mkqfb(a,b,c,gel(q,4));
313 6773330 : return gc_all(av, 2, &Q, U);
314 : }
315 :
316 : static GEN
317 1069563 : setq_b0(ulong a, ulong c, GEN D)
318 1069563 : { retmkqfb(utoipos(a), gen_0, utoipos(c), icopy(D)); }
319 : /* assume |sb| = 1 */
320 : static GEN
321 99091260 : setq(ulong a, ulong b, ulong c, long sb, GEN D)
322 99091260 : { retmkqfb(utoipos(a), sb==1? utoipos(b): utoineg(b), utoipos(c), icopy(D)); }
323 : /* 0 < a, c < 2^BIL, b = 0 */
324 : static GEN
325 958129 : qfi_red_1_b0(ulong a, ulong c, GEN D)
326 958129 : { return (a <= c)? setq_b0(a, c, D): setq_b0(c, a, D); }
327 :
328 : /* 0 < a, c < 2^BIL: single word affair */
329 : static GEN
330 100261353 : qfi_red_1(pari_sp av, GEN a, GEN b, GEN c, GEN D)
331 : {
332 : ulong ua, ub, uc;
333 : long sb;
334 : for(;;)
335 140925 : { /* at most twice */
336 100261353 : long lb = lgefint(b); /* <= 3 after first loop */
337 100261353 : if (lb == 2) return qfi_red_1_b0(a[2],c[2], D);
338 99303224 : if (lb == 3 && uel(b,2) <= (ulong)LONG_MAX) break;
339 139482 : REDB(a,&b,&c);
340 140504 : if (uel(a,2) <= uel(c,2))
341 : { /* lg(b) <= 3 but may be too large for itos */
342 0 : long s = signe(b);
343 0 : set_avma(av);
344 0 : if (!s) return qfi_red_1_b0(a[2], c[2], D);
345 0 : if (a[2] == c[2]) s = 1;
346 0 : return setq(a[2], b[2], c[2], s, D);
347 : }
348 140504 : swap(a,c); b = negi(b);
349 : }
350 : /* b != 0 */
351 99164583 : set_avma(av);
352 99164094 : ua = a[2];
353 99164094 : ub = sb = b[2]; if (signe(b) < 0) sb = -sb;
354 99164094 : uc = c[2];
355 99164094 : if (ua < ub)
356 37256666 : sREDB(ua, &sb, &uc);
357 61907428 : else if (ua == ub && sb < 0) sb = (long)ub;
358 176582281 : while(ua > uc)
359 : {
360 77382195 : lswap(ua,uc); sb = -sb;
361 77382195 : sREDB(ua, &sb, &uc);
362 : }
363 99200086 : if (!sb) return setq_b0(ua, uc, D);
364 : else
365 : {
366 99088652 : long s = 1;
367 99088652 : if (sb < 0)
368 : {
369 37182901 : sb = -sb;
370 37182901 : if (ua != uc) s = -1;
371 : }
372 99088652 : return setq(ua, sb, uc, s, D);
373 : }
374 : }
375 :
376 : static GEN
377 7 : qfi_rho(GEN x)
378 : {
379 7 : pari_sp av = avma;
380 7 : GEN a = gel(x,1), b = gel(x,2), c = gel(x,3);
381 7 : int fl = abscmpii(a, c);
382 7 : if (fl <= 0)
383 : {
384 7 : int fg = abscmpii(a, b);
385 7 : if (fg >= 0)
386 : {
387 7 : x = gcopy(x);
388 7 : if ((!fl || !fg) && signe(gel(x,2)) < 0) setsigne(gel(x,2), 1);
389 7 : return x;
390 : }
391 : }
392 0 : swap(a,c); b = negi(b);
393 0 : REDB(a, &b, &c);
394 0 : return gc_GEN(av, mkqfb(a,b,c, qfb_disc(x)));
395 : }
396 :
397 : /* qfr3 / qfr5 */
398 :
399 : /* t_QFB are unusable: D, sqrtD, isqrtD are recomputed all the time and the
400 : * logarithmic Shanks's distance is costly and hard to control.
401 : * qfr3 / qfr5 routines take a container of t_INTs (e.g a t_VEC) as argument,
402 : * at least 3 (resp. 5) components [it is a feature that they do not check the
403 : * precise type or length of the input]. They return a vector of length 3
404 : * (resp. 5). A qfr3 [a,b,c] contains the form coeffs, in a qfr5 [a,b,c, e,d]
405 : * the t_INT e is a binary exponent, d a t_REAL, coding the distance in
406 : * multiplicative form: the true distance is obtained from qfr5_dist.
407 : * All other qfr routines are obsolete (inefficient) wrappers */
408 :
409 : /* static functions are not stack-clean. Unless mentionned otherwise, public
410 : * functions are. */
411 :
412 : #define EMAX 22
413 : static void
414 10217928 : fix_expo(GEN x)
415 : {
416 10217928 : if (expo(gel(x,5)) >= (1L << EMAX)) {
417 0 : gel(x,4) = addiu(gel(x,4), 1);
418 0 : shiftr_inplace(gel(x,5), - (1L << EMAX));
419 : }
420 10217928 : }
421 :
422 : /* (1/2) log (|d| * 2^{e * 2^EMAX}). Not stack clean if e != 0 */
423 : GEN
424 184688 : qfr5_dist(GEN e, GEN d, long prec)
425 : {
426 184688 : GEN t = logr_abs(d);
427 184688 : if (signe(e)) {
428 0 : GEN u = mulir(e, mplog2(prec));
429 0 : shiftr_inplace(u, EMAX); t = addrr(t, u);
430 : }
431 184688 : shiftr_inplace(t, -1); return t;
432 : }
433 :
434 : static void
435 14152531 : rho_get_BC(GEN *B, GEN *C, GEN a, GEN b, GEN c, struct qfr_data *S)
436 : {
437 : GEN t, u, q;
438 14152531 : t = (abscmpii(S->isqrtD,c) >= 0)? S->isqrtD: absi_shallow(c);
439 14152531 : q = truedvmdii(addii(t, b), shifti(c,1), &u);
440 14152531 : *B = subii(t, u); /* t - ((t+b) % 2c) */
441 14152531 : *C = subii(a, mulii(q, subii(b, mulii(q,c))));
442 14152531 : }
443 : /* Not stack-clean */
444 : GEN
445 1139446 : qfr3_rho(GEN x, struct qfr_data *S)
446 : {
447 1139446 : GEN B, C, a = gel(x,1), b = gel(x,2), c = gel(x,3);
448 1139446 : rho_get_BC(&B, &C, a, b, c, S);
449 1139446 : return mkvec3(c, B, C);
450 : }
451 :
452 : /* Not stack-clean */
453 : GEN
454 8486681 : qfr5_rho(GEN x, struct qfr_data *S)
455 : {
456 8486681 : GEN B, C, a = gel(x,1), b = gel(x,2), c = gel(x,3), y;
457 8486681 : long sb = signe(b);
458 8486681 : rho_get_BC(&B, &C, a, b, c, S);
459 8486681 : y = mkvec5(c, B, C, gel(x,4), gel(x,5));
460 8486681 : if (sb) {
461 8482698 : GEN t = subii(sqri(b), S->D);
462 8482698 : if (sb < 0)
463 2509500 : t = divir(t, sqrr(subir(b,S->sqrtD)));
464 : else
465 5973198 : t = divri(sqrr(addir(b,S->sqrtD)), t);
466 : /* t = (b + sqrt(D)) / (b - sqrt(D)), evaluated stably */
467 8482698 : gel(y,5) = mulrr(t, gel(y,5)); fix_expo(y);
468 3983 : } else gel(y,5) = negr(gel(y,5));
469 8486681 : return y;
470 : }
471 :
472 : /* Not stack-clean */
473 : GEN
474 217728 : qfr_to_qfr5(GEN x, long prec)
475 217728 : { return mkvec5(gel(x,1),gel(x,2),gel(x,3),gen_0,real_1(prec)); }
476 :
477 : /* d0 = initial distance, x = [a,b,c, expo(d), d], d = exp(2*distance) */
478 : GEN
479 483 : qfr5_to_qfr(GEN x, GEN D, GEN d0)
480 : {
481 483 : if (d0)
482 : {
483 140 : GEN n = gel(x,4), d = absr(gel(x,5));
484 140 : if (signe(n))
485 : {
486 0 : n = addis(shifti(n, EMAX), expo(d));
487 0 : setexpo(d, 0); d = logr_abs(d);
488 0 : if (signe(n)) d = addrr(d, mulir(n, mplog2(lg(d0))));
489 0 : shiftr_inplace(d, -1);
490 0 : d0 = addrr(d0, d);
491 : }
492 140 : else if (!gequal1(d)) /* avoid loss of precision */
493 : {
494 91 : d = logr_abs(d);
495 91 : shiftr_inplace(d, -1);
496 91 : d0 = addrr(d0, d);
497 : }
498 : }
499 483 : x = qfr3_to_qfr(x, D);
500 483 : return d0 ? mkvec2(x,d0): x;
501 : }
502 :
503 : /* Not stack-clean */
504 : GEN
505 31920 : qfr3_to_qfr(GEN x, GEN d) { retmkqfb(gel(x,1), gel(x,2), gel(x,3), d); }
506 :
507 : static int
508 17712973 : ab_isreduced(GEN a, GEN b, GEN isqrtD)
509 : {
510 : GEN t;
511 17712973 : if (signe(b) <= 0 || abscmpii(b, isqrtD) > 0) return 0;
512 5271643 : t = addii_sign(isqrtD,1, shifti(a,1),-1); /* floor(sqrt(D)) - |2a| */
513 1099055 : return signe(t) < 0 ? abscmpii(b, t) >= 0
514 6370698 : : abscmpii(b, t) > 0;
515 : }
516 :
517 : /* Not stack-clean */
518 : GEN
519 1952895 : qfr5_red(GEN x, struct qfr_data *S)
520 : {
521 1952895 : pari_sp av = avma;
522 8465338 : while (!ab_isreduced(gel(x,1), gel(x,2), S->isqrtD))
523 : {
524 6512443 : x = qfr5_rho(x, S);
525 6512443 : if (gc_needed(av,2))
526 : {
527 0 : if (DEBUGMEM>1) pari_warn(warnmem,"qfr5_red");
528 0 : x = gc_GEN(av, x);
529 : }
530 : }
531 1952895 : return x;
532 : }
533 : /* Not stack-clean */
534 : GEN
535 1172847 : qfr3_red(GEN x, struct qfr_data *S)
536 : {
537 1172847 : pari_sp av = avma;
538 1172847 : GEN a = gel(x,1), b = gel(x,2), c = gel(x,3);
539 5699251 : while (!ab_isreduced(a, b, S->isqrtD))
540 : {
541 : GEN B, C;
542 4526404 : rho_get_BC(&B, &C, a, b, c, S);
543 4526404 : a = c; b = B; c = C;
544 4526404 : if (gc_needed(av,2))
545 : {
546 0 : if (DEBUGMEM>1) pari_warn(warnmem,"qfr3_red");
547 0 : (void)gc_all(av, 3, &a, &b, &c);
548 : }
549 : }
550 1172847 : return mkvec3(a, b, c);
551 : }
552 :
553 : void
554 2170 : qfr_data_init(GEN D, long prec, struct qfr_data *S)
555 : {
556 2170 : S->D = D;
557 2170 : S->sqrtD = sqrtr(itor(S->D,prec));
558 2170 : S->isqrtD = truncr(S->sqrtD);
559 2170 : }
560 :
561 : static GEN
562 140 : qfr5_init(GEN x, GEN d, struct qfr_data *S)
563 : {
564 140 : long prec = realprec(d), l = -expo(d);
565 140 : if (l < BITS_IN_LONG) l = BITS_IN_LONG;
566 140 : prec = maxss(prec, nbits2prec(l));
567 140 : S->D = qfb_disc(x);
568 140 : x = qfr_to_qfr5(x,prec);
569 140 : if (!S->sqrtD) S->sqrtD = sqrtr(itor(S->D,prec));
570 0 : else if (typ(S->sqrtD) != t_REAL) pari_err_TYPE("qfr_init",S->sqrtD);
571 :
572 140 : if (!S->isqrtD)
573 : {
574 126 : pari_sp av=avma;
575 : long e;
576 126 : S->isqrtD = gcvtoi(S->sqrtD,&e);
577 126 : if (e>-2) { set_avma(av); S->isqrtD = sqrti(S->D); }
578 : }
579 14 : else if (typ(S->isqrtD) != t_INT) pari_err_TYPE("qfr_init",S->isqrtD);
580 140 : return x;
581 : }
582 : static GEN
583 371 : qfr3_init(GEN x, struct qfr_data *S)
584 : {
585 371 : S->D = qfb_disc(x);
586 371 : if (!S->isqrtD) S->isqrtD = sqrti(S->D);
587 280 : else if (typ(S->isqrtD) != t_INT) pari_err_TYPE("qfr_init",S->isqrtD);
588 371 : return x;
589 : }
590 :
591 : #define qf_NOD 2
592 : #define qf_STEP 1
593 :
594 : static GEN
595 427 : qfr_red_basecase_i(GEN x, long flag, GEN isqrtD, GEN sqrtD)
596 : {
597 : struct qfr_data S;
598 427 : GEN d = NULL, y;
599 427 : if (typ(x)==t_VEC) { d = gel(x,2); x = gel(x,1); } else flag |= qf_NOD;
600 427 : S.sqrtD = sqrtD;
601 427 : S.isqrtD = isqrtD;
602 427 : y = (flag & qf_NOD)? qfr3_init(x, &S): qfr5_init(x, d, &S);
603 427 : switch(flag) {
604 63 : case 0: y = qfr5_red(y,&S); break;
605 336 : case qf_NOD: y = qfr3_red(y,&S); break;
606 21 : case qf_STEP: y = qfr5_rho(y,&S); break;
607 7 : case qf_STEP|qf_NOD: y = qfr3_rho(y,&S); break;
608 0 : default: pari_err_FLAG("qfbred");
609 : }
610 427 : return qfr5_to_qfr(y, qfb_disc(x), d);
611 : }
612 :
613 : static void
614 13379357 : qfr_rhosl2_i(GEN *pa, GEN *pb, GEN *pc, GEN *pu1, GEN *pu2, GEN *pv1,
615 : GEN *pv2, GEN rd)
616 : {
617 13379357 : GEN C = mpabs_shallow(*pc), t = addii(*pb, gmax_shallow(rd,C));
618 13379357 : GEN r, q = truedvmdii(t, shifti(C,1), &r);
619 13379357 : GEN a = *pa, b = *pb, c = *pc;
620 13379357 : if (signe(c) < 0) togglesign(q);
621 13379357 : *pa = *pc;
622 13379357 : *pb = subii(t, addii(r, *pb));
623 13379357 : *pc = subii(a, mulii(q, subii(b, mulii(q,c))));
624 13379357 : r = *pu1; *pu1 = *pv1; *pv1 = subii(mulii(q, *pv1), r);
625 13379357 : r = *pu2; *pu2 = *pv2; *pv2 = subii(mulii(q, *pv2), r);
626 13379357 : }
627 :
628 : static GEN
629 10810674 : qfr_rhosl2(GEN A, GEN rd)
630 : {
631 10810674 : GEN V = gel(A,1), M = gel(A,2);
632 10810674 : GEN a = gel(V,1), b = gel(V,2), c = gel(V,3), d = qfb_disc(V);
633 10810674 : GEN u1 = gcoeff(M,1,1), v1 = gcoeff(M,1,2);
634 10810674 : GEN u2 = gcoeff(M,2,1), v2 = gcoeff(M,2,2);
635 10810674 : qfr_rhosl2_i(&a,&b,&c, &u1,&u2,&v1,&v2, rd);
636 10810674 : return mkvec2(mkqfb(a,b,c,d), mkmat22(u1,v1,u2,v2));
637 : }
638 :
639 : static GEN
640 979701 : qfr_redsl2_basecase(GEN V, GEN rd)
641 : {
642 979701 : pari_sp av = avma;
643 979701 : GEN u1 = gen_1, u2 = gen_0, v1 = gen_0, v2 = gen_1;
644 979701 : GEN a = gel(V,1), b = gel(V,2), c = gel(V,3), d = qfb_disc(V);
645 3548384 : while (!ab_isreduced(a,b,rd))
646 : {
647 2568683 : qfr_rhosl2_i(&a,&b,&c, &u1,&u2,&v1,&v2, rd);
648 2568683 : if (gc_needed(av, 1))
649 : {
650 0 : if (DEBUGMEM>1) pari_warn(warnmem,"qfbredsl2");
651 0 : (void)gc_all(av, 7, &a,&b,&c,&u1,&u2,&v1,&v2);
652 : }
653 : }
654 979701 : return gc_GEN(av, mkvec2(mkqfb(a,b,c,d), mkmat22(u1,v1,u2,v2)));
655 : }
656 :
657 : /* fast reduction of qfb with positive coefficients, based on
658 : Arnold Schoenhage, Fast reduction and composition of binary quadratic forms,
659 : Proc. of Intern. Symp. on Symbolic and Algebraic Computation (Bonn) (S. M.
660 : Watt, ed.), ACM Press, 1991, pp. 128-133.
661 : <https://dl.acm.org/doi/pdf/10.1145/120694.120711>
662 : With thanks to Keegan Ryan
663 : BA20230927
664 : */
665 :
666 : /* pqfb: qf with positive coefficients */
667 :
668 : static int
669 4940497 : lti2n(GEN a, long m) { return signe(a) < 0 || expi(a) < m;}
670 :
671 : static GEN
672 1921185 : pqfbred_1(GEN Q, long m, GEN U)
673 : {
674 1921185 : GEN a = gel(Q,1), b = gel(Q,2), c = gel(Q,3), d = gel(Q,4);
675 1921185 : if (abscmpii(a, c) < 0)
676 : {
677 : GEN t, at, r;
678 960317 : GEN r2 = addii(shifti(a, m + 2), d);
679 960317 : long e2 = expi(r2);
680 960317 : r = int2n(signe(r2) < 0 || e2 < 2*m+2 ? m+1 : e2>>1);
681 960317 : t = truedivii(subii(b, r), shifti(a,1));
682 960317 : if (signe(t)==0) pari_err_BUG("pqfbred_1");
683 960317 : at = mulii(a,t);
684 960317 : c = addii(subii(c, mulii(b, t)), mulii(at, t));
685 960317 : b = subii(b, shifti(at,1));
686 960317 : gcoeff(U,1,2) = subii( gcoeff(U,1,2), mulii(gcoeff(U,1,1), t));
687 960317 : gcoeff(U,2,2) = subii( gcoeff(U,2,2), mulii(gcoeff(U,2,1), t));
688 : } else
689 : {
690 : GEN t, ct, r;
691 960868 : GEN r2 = addii(shifti(c, m + 2), d);
692 960868 : long e2 = expi(r2);
693 960868 : r = int2n(signe(r2) < 0 || e2 < 2*m+2 ? m+1 : e2>>1);
694 960868 : t = truedivii(subii(b, r), shifti(c,1));
695 960868 : if (signe(t)==0) pari_err_BUG("pqfbred_1");
696 960868 : ct = mulii(c, t);
697 960868 : a = addii(subii(a, mulii(b, t)), mulii(ct, t));
698 960868 : b = subii(b, shifti(ct, 1));
699 960868 : gcoeff(U,1,1) = subii(gcoeff(U,1,1), mulii(gcoeff(U,1,2), t));
700 960868 : gcoeff(U,2,1) = subii(gcoeff(U,2,1), mulii(gcoeff(U,2,2), t));
701 : }
702 1921185 : return mkqfb(a,b,c,d);
703 : }
704 :
705 : static int
706 2046352 : is_minimal(GEN Q, long m)
707 : {
708 2046352 : pari_sp av = avma;
709 2046352 : GEN a = gel(Q,1), b = gel(Q,2), c = gel(Q,3);
710 4940497 : return gc_bool(av, lti2n(addii(subii(a,b), c), m)
711 1927231 : || (lti2n(subii(b, shifti(a,1)), m+1)
712 966914 : && lti2n(subii(b, shifti(c,1)), m+1)));
713 : }
714 :
715 : static GEN
716 124020 : pqfbred_iter_1(GEN Q, ulong m, GEN U)
717 : {
718 124020 : pari_sp av = avma;
719 1923484 : while (!is_minimal(Q,m))
720 : {
721 1799464 : Q = pqfbred_1(Q, m, U);
722 1799464 : if (gc_needed(av, 1))
723 : {
724 0 : if (DEBUGMEM>1) pari_warn(warnmem,"pqfbred_iter_1, lc = %ld", expi(gel(Q,3)));
725 0 : (void)gc_all(av, 3, &Q, &gel(U,1), &gel(U,2));
726 : }
727 : }
728 124020 : return Q;
729 : }
730 :
731 : static GEN
732 61492 : pqfbred_basecase(GEN Q, ulong m, GEN *pt_U)
733 : {
734 61492 : pari_sp av = avma;
735 61492 : GEN U = matid(2);
736 61492 : Q = pqfbred_iter_1(Q, m, U);
737 61492 : *pt_U = U;
738 61492 : return gc_all(av, 2, &Q, pt_U);
739 : }
740 :
741 : static long
742 105804515 : qfb_maxexpi(GEN Q)
743 105804515 : { return 1+maxss(expi(gel(Q,1)), maxss(expi(gel(Q,2)), expi(gel(Q,3)))); }
744 :
745 : /* use asymptotically fast reduction ? */
746 : static int
747 105494496 : qfi_red_fast(GEN Q)
748 : {
749 105494496 : const long QFBRED_LIMIT = 9000;
750 105494496 : return 2*qfb_maxexpi(Q) - expi(gel(Q,4)) > QFBRED_LIMIT;
751 : }
752 :
753 : static long
754 125056 : qfb_minexpi(GEN Q)
755 : {
756 125056 : long m = minss(expi(gel(Q,1)), minss(expi(gel(Q,2)), expi(gel(Q,3))));
757 125056 : return m < 0 ? 0: m;
758 : }
759 :
760 : GEN
761 61913 : qfb3_SL2_apply(GEN q, GEN M)
762 : {
763 61913 : GEN a = gel(q,1), b = gel(q,2), c = gel(q,3);
764 61913 : GEN x = gcoeff(M,1,1), y = gcoeff(M,2,1);
765 61913 : GEN z = gcoeff(M,1,2), t = gcoeff(M,2,2);
766 61913 : GEN by = mulii(b,y), bt = mulii(b,t), bz = mulii(b,z);
767 61913 : GEN a2 = shifti(a,1), c2 = shifti(c,1);
768 :
769 61913 : GEN A1 = mulii(x, addii(mulii(a,x), by));
770 61913 : GEN A2 = mulii(c, sqri(y));
771 61913 : GEN B1 = mulii(x, addii(mulii(a2,z), bt));
772 61913 : GEN B2 = mulii(y, addii(mulii(c2,t), bz));
773 61913 : GEN C1 = mulii(z, addii(mulii(a,z), bt));
774 61913 : GEN C2 = mulii(c, sqri(t));
775 61913 : retmkvec3(addii(A1,A2), addii(B1,B2), addii(C1, C2));
776 : }
777 :
778 : static GEN
779 124020 : pqfbred_rec(GEN Q, long m, GEN *pt_U)
780 : {
781 124020 : pari_sp av = avma;
782 124020 : GEN U, Q0, Q1, QR, d = qfb_disc(Q);
783 124020 : long h, n = qfb_maxexpi(Q) - m;
784 124020 : int going_to_r8 = 0;
785 :
786 124020 : if (n < 170) return pqfbred_basecase(Q, m, pt_U);
787 62528 : if (qfb_minexpi(Q) <= m + 2) { U = matid(2); QR = Q; }
788 : else
789 : {
790 : long p, mm;
791 62528 : if (m <= n) { mm = m; p = 0; Q1 = Q; }
792 : else
793 : {
794 61878 : mm = n; p = m + 1 - n;
795 61878 : Q0 = mkvec3(remi2n(gel(Q,1),p), remi2n(gel(Q,2),p), remi2n(gel(Q,3),p));
796 61878 : Q1 = qfb3(shifti(gel(Q,1),-p), shifti(gel(Q,2),-p), shifti(gel(Q,3),-p));
797 : }
798 62528 : h = mm + (n>>1);
799 62528 : if (qfb_minexpi(Q1) <= h) { U = matid(2); QR = Q1; }
800 : else
801 62321 : QR = pqfbred_rec(Q1, h, &U);
802 184249 : while (qfb_maxexpi(QR) > h)
803 : {
804 122868 : if (is_minimal(QR, mm)) { going_to_r8 = 1; break; }
805 121721 : QR = pqfbred_1(QR, mm, U);
806 : }
807 62528 : if (!going_to_r8)
808 : {
809 : GEN V;
810 61381 : QR = pqfbred_rec(QR, mm, &V);
811 61381 : U = ZM2_mul(U,V);
812 : }
813 62528 : if (p > 0)
814 : {
815 61878 : GEN Q0U = qfb3_SL2_apply(Q0,U);
816 123756 : QR = mkqfb(addii(shifti(gel(QR,1), p), gel(Q0U,1)),
817 61878 : addii(shifti(gel(QR,2), p), gel(Q0U,2)),
818 61878 : addii(shifti(gel(QR,3), p), gel(Q0U,3)), d);
819 : }
820 : }
821 62528 : QR = pqfbred_iter_1(QR, m, U);
822 62528 : *pt_U = U; return gc_all(av, 2, &QR, pt_U);
823 : }
824 :
825 : static GEN
826 209575 : qfr_redsl2(GEN Q, GEN isqrtD)
827 : {
828 209575 : pari_sp av = avma;
829 209575 : if (!qfi_red_fast(Q))
830 209575 : return qfr_redsl2_basecase(Q, isqrtD);
831 : else
832 : {
833 0 : GEN a = gel(Q,1), b = gel(Q,2), c = gel(Q,3), d = gel(Q,4);
834 0 : GEN Qf, Qr, W, U, t = NULL;
835 0 : long sa = signe(a), sb;
836 0 : if (sa < 0) { a = negi(a); b = negi(b); c = negi(c); }
837 0 : if (signe(c) < 0)
838 : {
839 : GEN at;
840 0 : t = addiu(truedivii(subii(isqrtD,b),shifti(a,1)),1);
841 0 : at = mulii(a,t);
842 0 : c = addii(subii(c, mulii(b, t)), mulii(at, t));
843 0 : b = subii(b, shifti(at,1));
844 : }
845 0 : sb = signe(b);
846 0 : Qr = pqfbred_rec(mkqfb(a, sb < 0 ? negi(b): b, c, d), 0, &U);
847 0 : if (sa < 0)
848 0 : Qr = mkqfb(negi(gel(Qr,1)), negi(gel(Qr,2)), negi(gel(Qr,3)), gel(Qr,4));
849 0 : if (sb < 0)
850 : {
851 0 : gcoeff(U,2,1) = negi(gcoeff(U,2,1));
852 0 : gcoeff(U,2,2) = negi(gcoeff(U,2,2));
853 : }
854 0 : if (t)
855 : {
856 0 : gcoeff(U,1,1) = subii( gcoeff(U,1,1), mulii(gcoeff(U,2,1), t));
857 0 : gcoeff(U,1,2) = subii( gcoeff(U,1,2), mulii(gcoeff(U,2,2), t));
858 : }
859 0 : W = qfr_redsl2_basecase(Qr, isqrtD);
860 0 : Qf = gel(W,1);
861 0 : U = ZM2_mul(U,gel(W,2));
862 0 : return gc_GEN(av, mkvec2(Qf,U));
863 : }
864 : }
865 :
866 : static GEN
867 5182964 : qfi_redsl2(GEN Q)
868 : {
869 5182964 : pari_sp av = avma;
870 : GEN Qt, U;
871 5182964 : if (!qfi_red_fast(Q))
872 5182659 : Qt = qfi_redsl2_basecase(Q, &U);
873 : else
874 : {
875 311 : long sb = signe(gel(Q,2));
876 : GEN W;
877 311 : if (sb < 0) Q = mkqfb(gel(Q,1), negi(gel(Q,2)), gel(Q,3), gel(Q,4));
878 311 : Q = pqfbred_rec(Q, 0, &U);
879 311 : Qt = qfi_redsl2_basecase(Q, &W);
880 311 : U = ZM2_mul(U,W);
881 311 : if (sb < 0)
882 : {
883 173 : gcoeff(U,2,1) = negi(gcoeff(U,2,1));
884 173 : gcoeff(U,2,2) = negi(gcoeff(U,2,2));
885 : }
886 : }
887 5182993 : return gc_GEN(av, mkvec2(Qt,U));
888 : }
889 :
890 : GEN
891 4872647 : redimagsl2(GEN Q, GEN *U)
892 : {
893 4872647 : GEN q = qfi_redsl2(Q);
894 4872688 : *U = gel(q,2); return gel(q,1);
895 : }
896 :
897 : GEN
898 519893 : qfbredsl2(GEN q, GEN isD)
899 : {
900 : pari_sp av;
901 519893 : if (typ(q) != t_QFB) pari_err_TYPE("qfbredsl2",q);
902 519893 : if (qfb_is_qfi(q))
903 : {
904 310311 : if (isD) pari_err_TYPE("qfbredsl2", isD);
905 310311 : return qfi_redsl2(q);
906 : }
907 209582 : av = avma;
908 209582 : if (!isD) isD = sqrti(qfb_disc(q));
909 208068 : else if (typ(isD) != t_INT) pari_err_TYPE("qfbredsl2",isD);
910 209575 : return gc_upto(av, qfr_redsl2(q, isD));
911 : }
912 :
913 : /* not gc-clean */
914 : static GEN
915 427 : qfr_red_i(GEN Q, long flag, GEN isqrtD, GEN sqrtD)
916 : {
917 427 : if (typ(Q) != t_QFB || !qfi_red_fast(Q))
918 427 : return qfr_red_basecase_i(Q, flag, isqrtD, sqrtD);
919 : else
920 : {
921 0 : GEN a = gel(Q,1), b = gel(Q,2), c = gel(Q,3), d = gel(Q,4);
922 0 : GEN Qr, W, U, t = NULL;
923 0 : long sa = signe(a);
924 0 : if (sa < 0) { a = negi(a); b = negi(b); c = negi(c); }
925 0 : if (signe(c) < 0)
926 : {
927 : GEN at;
928 0 : t = addiu(truedivii(subii(isqrtD,b),shifti(a,1)),1);
929 0 : at = mulii(a,t);
930 0 : c = addii(subii(c, mulii(b, t)), mulii(at, t));
931 0 : b = subii(b, shifti(at,1));
932 : }
933 0 : Qr = pqfbred_rec(mkqfb(a, absi_shallow(b), c, d), 0, &U);
934 0 : if (sa < 0)
935 0 : Qr = mkqfb(negi(gel(Qr,1)), negi(gel(Qr,2)), negi(gel(Qr,3)), gel(Qr,4));
936 0 : W = qfr_red_basecase_i(Qr, flag, isqrtD, sqrtD);
937 0 : return gel(W,1);
938 : }
939 : }
940 :
941 : static GEN
942 63 : qfr_red_av(pari_sp av, GEN x)
943 63 : { return gc_GEN(av, qfr_red_i(x,0,NULL,NULL)); }
944 : GEN
945 0 : qfr_red(GEN x) { return qfr_red_av(avma, x); }
946 :
947 : static GEN
948 100294812 : qfi_red_basecase_av(pari_sp av, GEN q)
949 : {
950 100294812 : GEN a = gel(q,1), b = gel(q,2), c = gel(q,3), D = gel(q,4);
951 100294812 : long cmp, lc = lgefint(c);
952 :
953 100294812 : if (lgefint(a) == 3 && lc == 3) return qfi_red_1(av, a, b, c, D);
954 912889 : cmp = abscmpii(a, b);
955 919092 : if (cmp < 0)
956 436231 : REDB(a,&b,&c);
957 482861 : else if (cmp == 0 && signe(b) < 0)
958 27 : b = negi(b);
959 : for(;;)
960 : {
961 3126345 : cmp = abscmpii(a, c); if (cmp <= 0) break;
962 2932679 : lc = lgefint(a); /* lg(future c): we swap a & c next */
963 2932679 : if (lc == 3) return qfi_red_1(av, a, b, c, D);
964 2207253 : swap(a,c); b = negi(b); /* apply rho */
965 2207253 : REDB(a,&b,&c);
966 2207253 : if (gc_needed(av, 2))
967 : {
968 0 : if (DEBUGMEM>1) pari_warn(warnmem,"qfi_red, lc = %ld", lc);
969 0 : (void)gc_all(av, 3, &a,&b,&c);
970 : }
971 : }
972 193666 : if (cmp == 0 && signe(b) < 0) b = negi(b);
973 193666 : return gc_GEN(av, mkqfb(a, b, c, D));
974 : }
975 : static GEN
976 100103175 : qfi_red_av(pari_sp av, GEN Q)
977 : {
978 100103175 : if (qfi_red_fast(Q))
979 : {
980 : GEN U;
981 7 : if (signe(gel(Q,2)) < 0)
982 0 : Q = mkqfb(gel(Q,1), negi(gel(Q,2)), gel(Q,3), gel(Q,4));
983 7 : Q = pqfbred_rec(Q, 0, &U);
984 : }
985 100290787 : return qfi_red_basecase_av(av, Q);
986 : }
987 :
988 : GEN
989 23689044 : qfi_red(GEN q) { return qfi_red_av(avma, q); }
990 :
991 : GEN
992 55001 : qfbred0(GEN x, long flag, GEN isqrtD, GEN sqrtD)
993 : {
994 : pari_sp av;
995 55001 : GEN q = check_qfbext("qfbred",x);
996 55001 : if (qfb_is_qfi(q)) return (flag & qf_STEP)? qfi_rho(x): qfi_red(x);
997 364 : if (typ(x)==t_QFB) flag |= qf_NOD;
998 49 : else flag &= ~qf_NOD;
999 364 : av = avma;
1000 364 : return gc_GEN(av, qfr_red_i(x,flag,isqrtD,sqrtD));
1001 : }
1002 : /* t_QFB */
1003 : GEN
1004 15995408 : qfbred_i(GEN x) { return qfb_is_qfi(x)? qfi_red(x): qfr_red(x); }
1005 : GEN
1006 53188 : qfbred(GEN x) { return qfbred0(x, 0, NULL, NULL); }
1007 : /***********************************************************************/
1008 : /** **/
1009 : /** Composition **/
1010 : /** **/
1011 : /***********************************************************************/
1012 :
1013 : static void
1014 35124565 : qfb_sqr(GEN z, GEN x)
1015 : {
1016 : GEN c, d1, x2, v1, v2, c3, m, p1, r;
1017 :
1018 35124565 : d1 = bezout(gel(x,2),gel(x,1),&x2, NULL); /* usually 1 */
1019 35124448 : c = gel(x,3);
1020 35124448 : m = mulii(c,x2);
1021 35124139 : if (equali1(d1))
1022 26172741 : v1 = v2 = gel(x,1);
1023 : else
1024 : {
1025 8951449 : v1 = diviiexact(gel(x,1),d1);
1026 8951449 : v2 = mulii(v1, gcdii(d1,c)); /* = v1 iff x primitive */
1027 8951450 : c = mulii(c, d1);
1028 : }
1029 35124191 : togglesign(m);
1030 35124114 : r = modii(m,v2);
1031 35124019 : p1 = mulii(r, v1);
1032 35124148 : c3 = addii(c, mulii(r,addii(gel(x,2),p1)));
1033 35124134 : gel(z,1) = mulii(v1,v2);
1034 35124116 : gel(z,2) = addii(gel(x,2), shifti(p1,1));
1035 35124090 : gel(z,3) = diviiexact(c3,v2);
1036 35123981 : }
1037 : /* z <- x * y */
1038 : static void
1039 78161950 : qfb_comp(GEN z, GEN x, GEN y)
1040 : {
1041 : GEN n, c, d, y1, v1, v2, c3, m, p1, r;
1042 :
1043 78161950 : if (x == y) { qfb_sqr(z,x); return; }
1044 43872438 : n = shifti(subii(gel(y,2),gel(x,2)), -1);
1045 43678657 : v1 = gel(x,1);
1046 43678657 : v2 = gel(y,1);
1047 43678657 : c = gel(y,3);
1048 43678657 : d = bezout(v2,v1,&y1,NULL);
1049 43783675 : if (equali1(d))
1050 23138566 : m = mulii(y1,n);
1051 : else
1052 : {
1053 20682735 : GEN s = subii(gel(y,2), n);
1054 20680726 : GEN x2, y2, d1 = bezout(s,d,&x2,&y2); /* x2 s + y2 (x1 v1 + y1 v2) = d1 */
1055 20684556 : if (!equali1(d1))
1056 : {
1057 11430895 : v1 = diviiexact(v1,d1);
1058 11429832 : v2 = diviiexact(v2,d1); /* gcd = 1 iff x or y primitive */
1059 11429399 : v1 = mulii(v1, gcdii(c,gcdii(gel(x,3),gcdii(d1,n))));
1060 11429990 : c = mulii(c, d1);
1061 : }
1062 20683251 : m = addii(mulii(mulii(y1,y2),n), mulii(gel(y,3),x2));
1063 : }
1064 43758053 : togglesign(m);
1065 43686332 : r = modii(m, v1);
1066 43645463 : p1 = mulii(r, v2);
1067 43694619 : c3 = addii(c, mulii(r,addii(gel(y,2),p1)));
1068 43679205 : gel(z,1) = mulii(v1,v2);
1069 43689183 : gel(z,2) = addii(gel(y,2), shifti(p1,1));
1070 43674891 : gel(z,3) = diviiexact(c3,v1);
1071 : }
1072 :
1073 : /* not meant to be efficient */
1074 : static GEN
1075 84 : qfb_comp_gen(GEN x, GEN y)
1076 : {
1077 84 : GEN d1 = qfb_disc(x), d2 = qfb_disc(y);
1078 84 : GEN a1 = gel(x,1), b1 = gel(x,2), c1 = gel(x,3), n1;
1079 84 : GEN a2 = gel(y,1), b2 = gel(y,2), c2 = gel(y,3), n2;
1080 84 : GEN cx = content(x), cy = content(y), A, B, C, D, U, m, m2;
1081 :
1082 84 : if (!is_pm1(cx))
1083 : {
1084 14 : a1 = diviiexact(a1, cx); b1 = diviiexact(b1, cx);
1085 14 : c1 = diviiexact(c1, cx); d1 = diviiexact(d1, sqri(cx));
1086 : }
1087 84 : if (!is_pm1(cy))
1088 : {
1089 28 : a2 = diviiexact(a2, cy); c2 = diviiexact(c2, cy);
1090 28 : b2 = diviiexact(b2, cy); d2 = diviiexact(d2, sqri(cy));
1091 : }
1092 84 : D = gcdii(d1, d2); if (signe(d1) < 0) setsigne(D, -1);
1093 133 : if (!Z_issquareall(diviiexact(d1, D), &n1) ||
1094 84 : !Z_issquareall(diviiexact(d2, D), &n2)) return NULL;
1095 49 : A = mulii(a1, n2);
1096 49 : B = mulii(a2, n1);
1097 49 : C = shifti(addii(mulii(b1, n2), mulii(b2, n1)), -1);
1098 49 : U = ZV_extgcd(mkvec3(A, B, C));
1099 49 : m = gel(U,1); U = gmael(U,2,3);
1100 49 : A = mulii(diviiexact(mulii(a1,b2),m), gel(U,1));
1101 49 : B = mulii(diviiexact(mulii(a2,b1),m), gel(U,2));
1102 49 : C = addii(mulii(b1,b2), mulii(D, mulii(n1,n2)));
1103 49 : C = mulii(diviiexact(shifti(C,-1), m), gel(U,3));
1104 49 : B = addii(A, addii(B, C));
1105 49 : m2 = sqri(m);
1106 49 : A = diviiexact(mulii(a1, a2), m2);
1107 49 : C = diviiexact(shifti(subii(sqri(B),D), -2), A);
1108 49 : cx = mulii(cx, cy);
1109 49 : if (!is_pm1(cx))
1110 : {
1111 14 : A = mulii(A, cx); B = mulii(B, cx);
1112 14 : C = mulii(C, cx); D = mulii(D, sqri(cx));
1113 : }
1114 49 : return mkqfb(A, B, C, D);
1115 : }
1116 :
1117 : static GEN
1118 75421451 : qficomp0(GEN x, GEN y, int raw)
1119 : {
1120 75421451 : pari_sp av = avma;
1121 75421451 : GEN z = cgetg(5,t_QFB);
1122 75418931 : gel(z,4) = gel(x,4);
1123 75418931 : qfb_comp(z, x,y);
1124 75160883 : if (raw) return gc_GEN(av,z);
1125 75159105 : return qfi_red_av(av, z);
1126 : }
1127 : static GEN
1128 441 : qfrcomp0(GEN x, GEN y, int raw)
1129 : {
1130 441 : pari_sp av = avma;
1131 441 : GEN dx = NULL, dy = NULL;
1132 441 : GEN z = cgetg(5,t_QFB);
1133 441 : if (typ(x)==t_VEC) { dx = gel(x,2); x = gel(x,1); }
1134 441 : if (typ(y)==t_VEC) { dy = gel(y,2); y = gel(y,1); }
1135 441 : gel(z,4) = gel(x,4);
1136 441 : qfb_comp(z, x,y);
1137 441 : if (dx) z = mkvec2(z, dy? addrr(dx, dy): dx); else if (dy) z = mkvec2(z, dy);
1138 441 : if (raw) return gc_GEN(av, z);
1139 28 : return qfr_red_av(av, z);
1140 : }
1141 : /* same discriminant, no distance, no checks */
1142 : GEN
1143 27838982 : qfbcomp_i(GEN x, GEN y)
1144 27838982 : { return qfb_is_qfi(x)? qficomp0(x,y,0): qfrcomp0(x,y,0); }
1145 : GEN
1146 137554 : qfbcomp(GEN x, GEN y)
1147 : {
1148 137554 : GEN qx = check_qfbext("qfbcomp", x);
1149 137554 : GEN qy = check_qfbext("qfbcomp", y);
1150 137554 : if (!equalii(gel(qx,4),gel(qy,4)))
1151 : {
1152 63 : pari_sp av = avma;
1153 63 : GEN z = qfb_comp_gen(qx, qy);
1154 63 : if (typ(x) == t_VEC || typ(y) == t_VEC)
1155 7 : pari_err_IMPL("Shanks's distance in general composition");
1156 56 : if (!z) pari_err_OP("*",x,y);
1157 21 : return gc_upto(av, qfbred(z));
1158 : }
1159 137491 : return qfb_is_qfi(qx)? qficomp0(x,y,0): qfrcomp0(x,y,0);
1160 : }
1161 : /* same discriminant, no distance, no checks */
1162 : GEN
1163 0 : qfbcompraw_i(GEN x, GEN y)
1164 0 : { return qfb_is_qfi(x)? qficomp0(x,y,1): qfrcomp0(x,y,1); }
1165 : GEN
1166 2198 : qfbcompraw(GEN x, GEN y)
1167 : {
1168 2198 : GEN qx = check_qfbext("qfbcompraw", x);
1169 2198 : GEN qy = check_qfbext("qfbcompraw", y);
1170 2198 : if (!equalii(gel(qx,4),gel(qy,4)))
1171 : {
1172 21 : pari_sp av = avma;
1173 21 : GEN z = qfb_comp_gen(qx, qy);
1174 21 : if (typ(x) == t_VEC || typ(y) == t_VEC)
1175 0 : pari_err_IMPL("Shanks's distance in general composition");
1176 21 : if (!z) pari_err_OP("qfbcompraw",x,y);
1177 21 : return gc_GEN(av, z);
1178 : }
1179 2177 : if (!equalii(gel(qx,4),gel(qy,4))) pari_err_OP("qfbcompraw",x,y);
1180 2177 : return qfb_is_qfi(qx)? qficomp0(x,y,1): qfrcomp0(x,y,1);
1181 : }
1182 :
1183 : static GEN
1184 835011 : qfisqr0(GEN x, long raw)
1185 : {
1186 835011 : pari_sp av = avma;
1187 835011 : GEN z = cgetg(5,t_QFB);
1188 835011 : gel(z,4) = gel(x,4);
1189 835011 : qfb_sqr(z,x);
1190 835011 : if (raw) return gc_GEN(av,z);
1191 835011 : return qfi_red_av(av, z);
1192 : }
1193 : static GEN
1194 35 : qfrsqr0(GEN x, long raw)
1195 : {
1196 35 : pari_sp av = avma;
1197 35 : GEN dx = NULL, z = cgetg(5,t_QFB);
1198 35 : if (typ(x) == t_VEC) { dx = gel(x,2); x = gel(x,1); }
1199 35 : gel(z,4) = gel(x,4); qfb_sqr(z,x);
1200 35 : if (dx) z = mkvec2(z, shiftr(dx,1));
1201 35 : if (raw) return gc_GEN(av, z);
1202 35 : return qfr_red_av(av, z);
1203 : }
1204 : /* same discriminant, no distance, no checks */
1205 : GEN
1206 705056 : qfbsqr_i(GEN x)
1207 705056 : { return qfb_is_qfi(x)? qfisqr0(x,0): qfrsqr0(x,0); }
1208 : GEN
1209 129990 : qfbsqr(GEN x)
1210 : {
1211 129990 : GEN qx = check_qfbext("qfbsqr", x);
1212 129990 : return qfb_is_qfi(qx)? qfisqr0(x,0): qfrsqr0(x,0);
1213 : }
1214 :
1215 : static GEN
1216 6867 : qfr_1_by_disc(GEN D)
1217 : {
1218 : GEN y, r, s;
1219 6867 : check_quaddisc_real(D, NULL, "qfr_1_by_disc");
1220 6867 : y = cgetg(5,t_QFB);
1221 6867 : s = sqrtremi(D, &r); togglesign(r); /* s^2 - r = D */
1222 6867 : if (mpodd(r))
1223 : {
1224 3535 : s = subiu(s,1);
1225 3535 : r = subii(r, addiu(shifti(s, 1), 1));
1226 3535 : r = shifti(r, -2); set_avma((pari_sp)y); s = icopy(s);
1227 : }
1228 : else
1229 3332 : { r = shifti(r, -2); set_avma((pari_sp)s); }
1230 6867 : gel(y,1) = gen_1;
1231 6867 : gel(y,2) = s;
1232 6867 : gel(y,3) = icopy(r);
1233 6867 : gel(y,4) = icopy(D); return y;
1234 : }
1235 :
1236 : static GEN
1237 35 : qfr_disc(GEN x)
1238 35 : { return qfb_disc(typ(x)==t_VEC ? gel(x,1): x); }
1239 :
1240 : static GEN
1241 35 : qfr_1(GEN x)
1242 35 : { return qfr_1_by_disc(qfr_disc(x)); }
1243 :
1244 : static void
1245 0 : qfr_1_fill(GEN y, struct qfr_data *S)
1246 : {
1247 0 : pari_sp av = avma;
1248 0 : GEN y2 = S->isqrtD;
1249 0 : gel(y,1) = gen_1;
1250 0 : if (mod2(S->D) != mod2(y2)) y2 = subiu(y,1);
1251 0 : gel(y,2) = y2; av = avma;
1252 0 : gel(y,3) = gc_INT(av, shifti(subii(sqri(y2), S->D),-2));
1253 0 : }
1254 : static GEN
1255 0 : qfr5_1(struct qfr_data *S, long prec)
1256 : {
1257 0 : GEN y = cgetg(6, t_VEC);
1258 0 : qfr_1_fill(y, S);
1259 0 : gel(y,4) = gen_0;
1260 0 : gel(y,5) = real_1(prec); return y;
1261 : }
1262 : static GEN
1263 0 : qfr3_1(struct qfr_data *S)
1264 : {
1265 0 : GEN y = cgetg(4, t_VEC);
1266 0 : qfr_1_fill(y, S); return y;
1267 : }
1268 :
1269 : /* Assume D < 0 is the discriminant of a t_QFB */
1270 : static GEN
1271 775157 : qfi_1_by_disc(GEN D)
1272 : {
1273 775157 : GEN b,c, y = cgetg(5,t_QFB);
1274 775157 : quadpoly_bc(D, mod2(D), &b,&c);
1275 775157 : if (b == gen_m1) b = gen_1;
1276 775157 : gel(y,1) = gen_1;
1277 775157 : gel(y,2) = b;
1278 775157 : gel(y,3) = c;
1279 775157 : gel(y,4) = icopy(D); return y;
1280 : }
1281 : static GEN
1282 762876 : qfi_1(GEN x)
1283 : {
1284 762876 : if (typ(x) != t_QFB) pari_err_TYPE("qfi_1",x);
1285 762876 : return qfi_1_by_disc(qfb_disc(x));
1286 : }
1287 :
1288 : GEN
1289 0 : qfb_1(GEN x) { return qfb_is_qfi(x) ? qfi_1(x): qfr_1(x); }
1290 :
1291 : static GEN
1292 14284612 : _qfimul(void *E, GEN x, GEN y) { (void) E; return qficomp0(x,y,0); }
1293 : static GEN
1294 33161253 : _qfisqr(void *E, GEN x) { (void) E; return qficomp0(x,x,0); }
1295 : static GEN
1296 7 : _qfimulraw(void *E, GEN x, GEN y) { (void) E; return qficomp0(x,y,1); }
1297 : static GEN
1298 7 : _qfisqrraw(void *E, GEN x) { (void) E; return qficomp0(x,x,1); }
1299 :
1300 : static GEN
1301 7 : qfipowraw(GEN x, long n)
1302 : {
1303 7 : pari_sp av = avma;
1304 : GEN y;
1305 7 : if (!n) return qfi_1(x);
1306 7 : if (n== 1) return gcopy(x);
1307 7 : if (n==-1) { x = gcopy(x); togglesign(gel(x,2)); return x; }
1308 7 : if (n < 0) x = qfb_inv(x);
1309 7 : y = gen_powu(x, labs(n), NULL, &_qfisqrraw, &_qfimulraw);
1310 7 : return gc_GEN(av,y);
1311 : }
1312 :
1313 : static GEN
1314 16758283 : qfipow(GEN x, GEN n)
1315 : {
1316 16758283 : pari_sp av = avma;
1317 : GEN y;
1318 16758283 : long s = signe(n);
1319 16758283 : if (!s) return qfi_1(x);
1320 15995407 : if (s < 0) x = qfb_inv(x);
1321 15995407 : y = gen_pow(qfbred_i(x), n, NULL, &_qfisqr, &_qfimul);
1322 15995406 : return gc_GEN(av,y);
1323 : }
1324 :
1325 : static long
1326 412328 : parteucl(GEN L, GEN *d, GEN *v3, GEN *v, GEN *v2)
1327 : {
1328 : long z;
1329 412328 : *v = gen_0; *v2 = gen_1;
1330 4351417 : for (z=0; abscmpii(*v3,L) > 0; z++)
1331 : {
1332 3939089 : GEN t3, t2 = subii(*v, mulii(truedvmdii(*d,*v3,&t3),*v2));
1333 3939089 : *v = *v2; *d = *v3; *v2 = t2; *v3 = t3;
1334 : }
1335 412328 : return z;
1336 : }
1337 :
1338 : /* composition: Shanks' NUCOMP & NUDUPL */
1339 : /* L = floor((|d|/4)^(1/4)) */
1340 : GEN
1341 400722 : nucomp(GEN x, GEN y, GEN L)
1342 : {
1343 400722 : pari_sp av = avma;
1344 : long z;
1345 : GEN a, a1, a2, b2, b, d, d1, g, n, p1, q1, q2, s, u, u1, v, v2, v3, Q;
1346 :
1347 400722 : if (x==y) return nudupl(x,L);
1348 400680 : if (!is_qfi(x)) pari_err_TYPE("nucomp",x);
1349 400680 : if (!is_qfi(y)) pari_err_TYPE("nucomp",y);
1350 :
1351 400680 : if (abscmpii(gel(x,1),gel(y,1)) < 0) swap(x, y);
1352 400680 : s = shifti(addii(gel(x,2),gel(y,2)), -1);
1353 400680 : n = subii(gel(y,2), s);
1354 400680 : a1 = gel(x,1);
1355 400680 : a2 = gel(y,1); d = bezout(a2,a1,&u,&v);
1356 400680 : if (equali1(d)) { a = negi(mulii(u,n)); d1 = d; }
1357 163576 : else if (dvdii(s,d)) /* d | s */
1358 : {
1359 83503 : a = negi(mulii(u,n)); d1 = d;
1360 83503 : a1 = diviiexact(a1, d1);
1361 83503 : a2 = diviiexact(a2, d1);
1362 83503 : s = diviiexact(s, d1);
1363 : }
1364 : else
1365 : {
1366 : GEN p2, l;
1367 80073 : d1 = bezout(s,d,&u1,NULL);
1368 80073 : if (!equali1(d1))
1369 : {
1370 2044 : a1 = diviiexact(a1,d1);
1371 2044 : a2 = diviiexact(a2,d1);
1372 2044 : s = diviiexact(s,d1);
1373 2044 : d = diviiexact(d,d1);
1374 : }
1375 80073 : p1 = remii(gel(x,3),d);
1376 80073 : p2 = remii(gel(y,3),d);
1377 80073 : l = modii(mulii(negi(u1), addii(mulii(u,p1),mulii(v,p2))), d);
1378 80073 : a = subii(mulii(l,diviiexact(a1,d)), mulii(u,diviiexact(n,d)));
1379 : }
1380 400680 : a = modii(a,a1); p1 = subii(a,a1); if (abscmpii(a,p1) > 0) a = p1;
1381 400680 : d = a1; v3 = a; z = parteucl(L, &d,&v3, &v,&v2);
1382 400680 : Q = cgetg(5,t_QFB);
1383 400680 : if (!z) {
1384 37632 : g = diviiexact(addii(mulii(v3,s),gel(y,3)), d);
1385 37632 : b = a2;
1386 37632 : b2 = gel(y,2);
1387 37632 : v2 = d1;
1388 37632 : gel(Q,1) = mulii(d,b);
1389 : } else {
1390 : GEN e, q3, q4;
1391 363048 : if (z&1) { v3 = negi(v3); v2 = negi(v2); }
1392 363048 : b = diviiexact(addii(mulii(a2,d), mulii(n,v)), a1);
1393 363048 : e = diviiexact(addii(mulii(s,d),mulii(gel(y,3),v)), a1);
1394 363048 : q3 = mulii(e,v2);
1395 363048 : q4 = subii(q3,s);
1396 363048 : b2 = addii(q3,q4);
1397 363048 : g = diviiexact(q4,v);
1398 363048 : if (!equali1(d1)) { v2 = mulii(d1,v2); v = mulii(d1,v); b2 = mulii(d1,b2); }
1399 363048 : gel(Q,1) = addii(mulii(d,b), mulii(e,v));
1400 : }
1401 400680 : q1 = mulii(b, v3);
1402 400680 : q2 = addii(q1,n);
1403 400680 : gel(Q,2) = addii(b2, z? addii(q1,q2): shifti(q1, 1));
1404 400680 : gel(Q,3) = addii(mulii(v3,diviiexact(q2,d)), mulii(g,v2));
1405 400680 : gel(Q,4) = gel(x,4);
1406 400680 : return qfi_red_av(av, Q);
1407 : }
1408 :
1409 : GEN
1410 11648 : nudupl(GEN x, GEN L)
1411 : {
1412 11648 : pari_sp av = avma;
1413 : long z;
1414 : GEN u, v, d, d1, p1, a, b, c, a2, b2, c2, Q, v2, v3, g;
1415 :
1416 11648 : if (!is_qfi(x)) pari_err_TYPE("nudupl",x);
1417 11648 : a = gel(x,1);
1418 11648 : b = gel(x,2);
1419 11648 : d1 = bezout(b,a, &u,NULL);
1420 11648 : if (!equali1(d1))
1421 : {
1422 4620 : a = diviiexact(a, d1);
1423 4620 : b = diviiexact(b, d1);
1424 : }
1425 11648 : c = modii(negi(mulii(u,gel(x,3))), a);
1426 11648 : p1 = subii(c,a); if (abscmpii(c,p1) > 0) c = p1;
1427 11648 : d = a; v3 = c; z = parteucl(L, &d,&v3, &v,&v2);
1428 11648 : a2 = sqri(d);
1429 11648 : c2 = sqri(v3);
1430 11648 : Q = cgetg(5,t_QFB);
1431 11648 : if (!z) {
1432 1281 : g = diviiexact(addii(mulii(v3,b),gel(x,3)), d);
1433 1281 : b2 = gel(x,2);
1434 1281 : v2 = d1;
1435 1281 : gel(Q,1) = a2;
1436 : } else {
1437 : GEN e;
1438 10367 : if (z&1) { v = negi(v); d = negi(d); }
1439 10367 : e = diviiexact(addii(mulii(gel(x,3),v), mulii(b,d)), a);
1440 10367 : g = diviiexact(subii(mulii(e,v2), b), v);
1441 10367 : b2 = addii(mulii(e,v2), mulii(v,g));
1442 10367 : if (!equali1(d1)) { b2 = mulii(d1,b2); v = mulii(d1,v); v2 = mulii(d1,v2); }
1443 10367 : gel(Q,1) = addii(a2, mulii(e,v));
1444 : }
1445 11648 : gel(Q,2) = addii(b2, subii(sqri(addii(d,v3)), addii(a2,c2)));
1446 11648 : gel(Q,3) = addii(c2, mulii(g,v2));
1447 11648 : gel(Q,4) = gel(x,4);
1448 11648 : return qfi_red_av(av, Q);
1449 : }
1450 :
1451 : static GEN
1452 4739 : mul_nucomp(void *l, GEN x, GEN y) { return nucomp(x, y, (GEN)l); }
1453 : static GEN
1454 11606 : mul_nudupl(void *l, GEN x) { return nudupl(x, (GEN)l); }
1455 : GEN
1456 1008 : nupow(GEN x, GEN n, GEN L)
1457 : {
1458 : pari_sp av;
1459 : GEN y, D;
1460 :
1461 1008 : if (typ(n) != t_INT) pari_err_TYPE("nupow",n);
1462 1008 : if (!is_qfi(x)) pari_err_TYPE("nupow",x);
1463 1008 : if (gequal1(n)) return gcopy(x);
1464 1008 : av = avma;
1465 1008 : D = qfb_disc(x);
1466 1008 : y = qfi_1_by_disc(D);
1467 1008 : if (!signe(n)) return y;
1468 959 : if (!L) L = sqrtnint(absi_shallow(D), 4);
1469 959 : y = gen_pow_i(x, n, (void*)L, &mul_nudupl, &mul_nucomp);
1470 959 : if (signe(n) < 0
1471 35 : && !absequalii(gel(y,1),gel(y,2))
1472 35 : && !absequalii(gel(y,1),gel(y,3))) togglesign(gel(y,2));
1473 959 : return gc_GEN(av, y);
1474 : }
1475 :
1476 : /* Not stack-clean */
1477 : GEN
1478 1735230 : qfr5_compraw(GEN x, GEN y)
1479 : {
1480 1735230 : GEN z = cgetg(6,t_VEC); qfb_comp(z,x,y);
1481 1735230 : if (x == y)
1482 : {
1483 34552 : gel(z,4) = shifti(gel(x,4),1);
1484 34552 : gel(z,5) = sqrr(gel(x,5));
1485 : }
1486 : else
1487 : {
1488 1700678 : gel(z,4) = addii(gel(x,4),gel(y,4));
1489 1700678 : gel(z,5) = mulrr(gel(x,5),gel(y,5));
1490 : }
1491 1735230 : fix_expo(z); return z;
1492 : }
1493 : GEN
1494 1735216 : qfr5_comp(GEN x, GEN y, struct qfr_data *S)
1495 1735216 : { return qfr5_red(qfr5_compraw(x, y), S); }
1496 : /* Not stack-clean */
1497 : GEN
1498 1009397 : qfr3_compraw(GEN x, GEN y)
1499 : {
1500 1009397 : GEN z = cgetg(4,t_VEC); qfb_comp(z,x,y);
1501 1009397 : return z;
1502 : }
1503 : GEN
1504 1009397 : qfr3_comp(GEN x, GEN y, struct qfr_data *S)
1505 1009397 : { return qfr3_red(qfr3_compraw(x,y), S); }
1506 :
1507 : /* m > 0. Not stack-clean */
1508 : static GEN
1509 7 : qfr5_powraw(GEN x, long m)
1510 : {
1511 7 : GEN y = NULL;
1512 14 : for (; m; m >>= 1)
1513 : {
1514 14 : if (m&1) y = y? qfr5_compraw(y,x): x;
1515 14 : if (m == 1) break;
1516 7 : x = qfr5_compraw(x,x);
1517 : }
1518 7 : return y;
1519 : }
1520 :
1521 : /* return x^n. Not stack-clean */
1522 : GEN
1523 21 : qfr5_pow(GEN x, GEN n, struct qfr_data *S)
1524 : {
1525 21 : GEN y = NULL;
1526 21 : long i, m, s = signe(n);
1527 21 : if (!s) return qfr5_1(S, lg(gel(x,5)));
1528 21 : if (s < 0) x = qfb_inv(x);
1529 42 : for (i=lgefint(n)-1; i>1; i--)
1530 : {
1531 21 : m = n[i];
1532 56 : for (; m; m>>=1)
1533 : {
1534 56 : if (m&1) y = y? qfr5_comp(y,x,S): x;
1535 56 : if (m == 1 && i == 2) break;
1536 35 : x = qfr5_comp(x,x,S);
1537 : }
1538 : }
1539 21 : return y;
1540 : }
1541 : /* m > 0; return x^m. Not stack-clean */
1542 : static GEN
1543 0 : qfr3_powraw(GEN x, long m)
1544 : {
1545 0 : GEN y = NULL;
1546 0 : for (; m; m>>=1)
1547 : {
1548 0 : if (m&1) y = y? qfr3_compraw(y,x): x;
1549 0 : if (m == 1) break;
1550 0 : x = qfr3_compraw(x,x);
1551 : }
1552 0 : return y;
1553 : }
1554 : /* return x^n. Not stack-clean */
1555 : GEN
1556 4557 : qfr3_pow(GEN x, GEN n, struct qfr_data *S)
1557 : {
1558 4557 : GEN y = NULL;
1559 4557 : long i, m, s = signe(n);
1560 4557 : if (!s) return qfr3_1(S);
1561 4557 : if (s < 0) x = qfb_inv(x);
1562 9130 : for (i=lgefint(n)-1; i>1; i--)
1563 : {
1564 4573 : m = n[i];
1565 5312 : for (; m; m>>=1)
1566 : {
1567 5292 : if (m&1) y = y? qfr3_comp(y,x,S): x;
1568 5292 : if (m == 1 && i == 2) break;
1569 739 : x = qfr3_comp(x,x,S);
1570 : }
1571 : }
1572 4557 : return y;
1573 : }
1574 :
1575 : static GEN
1576 7 : qfrinvraw(GEN x)
1577 : {
1578 7 : if (typ(x) == t_VEC) retmkvec2(qfbinv(gel(x,1)), negr(gel(x,2)));
1579 7 : return qfbinv(x);
1580 : }
1581 : static GEN
1582 14 : qfrpowraw(GEN x, long n)
1583 : {
1584 14 : struct qfr_data S = { NULL, NULL, NULL };
1585 14 : pari_sp av = avma;
1586 14 : if (n==1) return gcopy(x);
1587 14 : if (n==-1) return qfrinvraw(x);
1588 7 : if (typ(x)==t_QFB)
1589 : {
1590 0 : GEN D = qfb_disc(x);
1591 0 : if (!n) return qfr_1(x);
1592 0 : if (n < 0) { x = qfb_inv(x); n = -n; }
1593 0 : x = qfr3_powraw(x, n);
1594 0 : x = qfr3_to_qfr(x, D);
1595 : }
1596 : else
1597 : {
1598 7 : GEN d0 = gel(x,2);
1599 7 : x = gel(x,1);
1600 7 : if (!n) retmkvec2(qfr_1(x), real_0(precision(d0)));
1601 7 : if (n < 0) { x = qfb_inv(x); n = -n; }
1602 7 : x = qfr5_init(x, d0, &S);
1603 7 : if (labs(n) != 1) x = qfr5_powraw(x, n);
1604 7 : x = qfr5_to_qfr(x, S.D, mulrs(d0,n));
1605 : }
1606 7 : return gc_GEN(av, x);
1607 : }
1608 : static GEN
1609 112 : qfrpow(GEN x, GEN n)
1610 : {
1611 112 : struct qfr_data S = { NULL, NULL, NULL };
1612 112 : long s = signe(n);
1613 112 : pari_sp av = avma;
1614 112 : if (typ(x)==t_QFB)
1615 : {
1616 42 : if (!s) return qfr_1(x);
1617 28 : if (s < 0) x = qfb_inv(x);
1618 28 : x = qfr3_init(x, &S);
1619 28 : x = is_pm1(n)? qfr3_red(x, &S): qfr3_pow(x, n, &S);
1620 28 : x = qfr3_to_qfr(x, S.D);
1621 : }
1622 : else
1623 : {
1624 70 : GEN d0 = gel(x,2);
1625 70 : x = gel(x,1);
1626 70 : if (!s) retmkvec2(qfr_1(x), real_0(precision(d0)));
1627 49 : if (s < 0) x = qfb_inv(x);
1628 49 : x = qfr5_init(x, d0, &S);
1629 49 : x = is_pm1(n)? qfr5_red(x, &S): qfr5_pow(x, n, &S);
1630 49 : x = qfr5_to_qfr(x, S.D, mulri(d0,n));
1631 : }
1632 77 : return gc_GEN(av, x);
1633 : }
1634 : GEN
1635 21 : qfbpowraw(GEN x, long n)
1636 : {
1637 21 : GEN q = check_qfbext("qfbpowraw",x);
1638 21 : return qfb_is_qfi(q)? qfipowraw(x,n): qfrpowraw(x,n);
1639 : }
1640 : /* same discriminant, no distance, no checks */
1641 : GEN
1642 15277246 : qfbpow_i(GEN x, GEN n) { return qfb_is_qfi(x)? qfipow(x,n): qfrpow(x,n); }
1643 : GEN
1644 1481149 : qfbpow(GEN x, GEN n)
1645 : {
1646 1481149 : GEN q = check_qfbext("qfbpow",x);
1647 1481149 : return qfb_is_qfi(q)? qfipow(x,n): qfrpow(x,n);
1648 : }
1649 : GEN
1650 1382173 : qfbpows(GEN x, long n)
1651 : {
1652 1382173 : long N[] = { evaltyp(t_INT) | _evallg(3), 0, 0};
1653 1382173 : affsi(n, N); return qfbpow(x, N);
1654 : }
1655 :
1656 : /* Prime forms attached to prime ideals of degree 1 */
1657 :
1658 : /* assume x != 0 a t_INT, p > 0
1659 : * Return a t_QFB, but discriminant sign is not checked: can be used for
1660 : * real forms as well */
1661 : GEN
1662 12682908 : primeform_u(GEN x, ulong p)
1663 : {
1664 12682908 : GEN c, y = cgetg(5, t_QFB);
1665 12682826 : pari_sp av = avma;
1666 : ulong b;
1667 : long s;
1668 :
1669 12682826 : s = mod8(x); if (signe(x) < 0 && s) s = 8-s;
1670 : /* 2 or 3 mod 4 */
1671 12682730 : if (s & 2) pari_err_DOMAIN("primeform", "disc % 4", ">",gen_1, x);
1672 12682960 : if (p == 2) {
1673 4170304 : switch(s) {
1674 589604 : case 0: b = 0; break;
1675 3229229 : case 1: b = 1; break;
1676 351472 : case 4: b = 2; break;
1677 0 : default: pari_err_SQRTN("primeform", mkintmod(x,utoi(p)) );
1678 0 : b = 0; /* -Wall */
1679 : }
1680 4170305 : c = shifti(subsi(s,x), -3);
1681 : } else {
1682 8512656 : b = Fl_sqrt(umodiu(x,p), p);
1683 8514070 : if (b == ~0UL) pari_err_SQRTN("primeform", mkintmod(x,utoi(p)) );
1684 : /* mod(b) != mod2(x) ? */
1685 8514926 : if ((b ^ s) & 1) b = p - b;
1686 8514926 : c = diviuexact(shifti(subii(sqru(b), x), -2), p);
1687 : }
1688 12670127 : gel(y,3) = gc_INT(av, c);
1689 12679112 : gel(y,4) = icopy(x);
1690 12682266 : gel(y,2) = utoi(b);
1691 12682266 : gel(y,1) = utoipos(p); return y;
1692 : }
1693 :
1694 : /* special case: p = 1 return unit form */
1695 : GEN
1696 135782 : primeform(GEN x, GEN p)
1697 : {
1698 135782 : const char *f = "primeform";
1699 : pari_sp av;
1700 135782 : long s, sx = signe(x), sp = signe(p);
1701 : GEN y, b, absp;
1702 :
1703 135782 : if (typ(x) != t_INT) pari_err_TYPE(f,x);
1704 135782 : if (typ(p) != t_INT) pari_err_TYPE(f,p);
1705 135782 : if (!sp) pari_err_DOMAIN(f,"p","=",gen_0,p);
1706 135782 : if (!sx) pari_err_DOMAIN(f,"D","=",gen_0,x);
1707 135782 : if (lgefint(p) == 3)
1708 : {
1709 135768 : ulong pp = p[2];
1710 135768 : if (pp == 1) {
1711 18105 : if (sx < 0) {
1712 : long r;
1713 11273 : if (sp < 0) pari_err_IMPL("negative definite t_QFB");
1714 11273 : r = mod4(x);
1715 11273 : if (r && r != 3) pari_err_DOMAIN(f,"disc % 4",">", gen_1,x);
1716 11273 : return qfi_1_by_disc(x);
1717 : }
1718 6832 : y = qfr_1_by_disc(x);
1719 6832 : if (sp < 0) { gel(y,1) = gen_m1; togglesign(gel(y,3)); }
1720 6832 : return y;
1721 : }
1722 117663 : y = primeform_u(x, pp);
1723 117656 : if (sx < 0) {
1724 89957 : if (sp < 0) pari_err_IMPL("negative definite t_QFB");
1725 89957 : return y;
1726 : }
1727 27699 : if (sp < 0) { togglesign(gel(y,1)); togglesign(gel(y,3)); }
1728 27699 : return gcopy( qfr3_to_qfr(y, x) );
1729 : }
1730 14 : s = mod8(x);
1731 14 : if (sx < 0)
1732 : {
1733 7 : if (sp < 0) pari_err_IMPL("negative definite t_QFB");
1734 7 : if (s) s = 8-s;
1735 : }
1736 14 : y = cgetg(5, t_QFB);
1737 : /* 2 or 3 mod 4 */
1738 14 : if (s & 2) pari_err_DOMAIN(f, "disc % 4", ">",gen_1, x);
1739 14 : absp = absi_shallow(p); av = avma;
1740 14 : b = Fp_sqrt(x, absp); if (!b) pari_err_SQRTN(f, mkintmod(x,absp));
1741 14 : s &= 1; /* s = x mod 2 */
1742 : /* mod(b) != mod2(x) ? [Warning: we may have b == 0] */
1743 14 : if ((!signe(b) && s) || mod2(b) != s) b = gc_INT(av, subii(absp,b));
1744 :
1745 14 : av = avma;
1746 14 : gel(y,3) = gc_INT(av, diviiexact(shifti(subii(sqri(b), x), -2), p));
1747 14 : gel(y,4) = icopy(x);
1748 14 : gel(y,2) = b;
1749 14 : gel(y,1) = icopy(p);
1750 14 : return y;
1751 : }
1752 :
1753 : static GEN
1754 2620771 : normforms(GEN D, GEN fa)
1755 : {
1756 : long i, j, k, lB, aN, sa;
1757 : GEN a, L, V, B, N, N2;
1758 2620771 : int D_odd = mpodd(D);
1759 2620771 : a = typ(fa) == t_INT ? fa: typ(fa) == t_VEC? gel(fa,1): factorback(fa);
1760 2620771 : sa = signe(a);
1761 2620771 : if (sa==0 || (signe(D)<0 && sa<0)) return NULL;
1762 1203972 : V = D_odd? Zn_quad_roots(fa, gen_1, shifti(subsi(1, D), -2))
1763 2551765 : : Zn_quad_roots(fa, gen_0, negi(shifti(D, -2)));
1764 2551766 : if (!V) return NULL;
1765 511966 : N = gel(V,1); B = gel(V,2); lB = lg(B);
1766 511966 : N2 = shifti(N,1);
1767 511966 : aN = itou(diviiexact(a, N)); /* |a|/N */
1768 511965 : L = cgetg((lB-1)*aN+1, t_VEC);
1769 2360556 : for (k = 1, i = 1; i < lB; i++)
1770 : {
1771 1848590 : GEN b = shifti(gel(B,i), 1), c, C;
1772 1848588 : if (D_odd) b = addiu(b , 1);
1773 1848588 : c = diviiexact(shifti(subii(sqri(b), D), -2), a);
1774 1848585 : for (j = 0;; b = addii(b, N2))
1775 : {
1776 2216659 : gel(L, k++) = mkqfb(a, b, c, D);
1777 2216668 : if (++j == aN) break;
1778 368074 : C = addii(b, N); if (aN > 1) C = diviuexact(C, aN);
1779 368074 : c = sa > 0? addii(c, C): subii(c, C);
1780 : }
1781 : }
1782 511966 : return L;
1783 : }
1784 :
1785 : /* Let M and N in SL2(Z), return (N*M^-1)[,1] */
1786 : static GEN
1787 344317 : SL2_div_mul_e1(GEN N, GEN M)
1788 : {
1789 344317 : GEN b = gcoeff(M,2,1), d = gcoeff(M,2,2);
1790 344317 : GEN A = mulii(gcoeff(N,1,1), d), B = mulii(gcoeff(N,1,2), b);
1791 344311 : GEN C = mulii(gcoeff(N,2,1), d), D = mulii(gcoeff(N,2,2), b);
1792 344310 : retmkvec2(subii(A,B), subii(C,D));
1793 : }
1794 : static GEN
1795 1445667 : qfisolve_normform(GEN Q, GEN P)
1796 : {
1797 1445667 : GEN a = gel(Q,1), N = gel(Q,2);
1798 1445667 : GEN M, b = qfi_redsl2_basecase(P, &M);
1799 1445682 : if (!qfb_equal(a,b)) return NULL;
1800 102124 : return SL2_div_mul_e1(N,M);
1801 : }
1802 :
1803 : /* Test equality modulo GL2 of two reduced forms */
1804 : static int
1805 61068 : GL2_qfb_equal(GEN a, GEN b)
1806 : {
1807 61068 : return equalii(gel(a,1),gel(b,1))
1808 11361 : && absequalii(gel(a,2),gel(b,2))
1809 72429 : && equalii(gel(a,3),gel(b,3));
1810 : }
1811 :
1812 : /* Q(u,v) = p; if s < 0 return that solution; else the set of all solutions */
1813 : static GEN
1814 48042 : allsols(GEN Q, long s, GEN u, GEN v)
1815 : {
1816 48042 : GEN w = mkvec2(u, v), b;
1817 48030 : if (signe(v) < 0) { u = negi(u); v = negi(v); } /* normalize for v >= 0 */
1818 48030 : w = mkvec2(u, v); if (s < 0) return w;
1819 41386 : if (!s) return mkvec(w);
1820 38957 : b = gel(Q,2); /* sum of the 2 solutions (if they exist) is -bv / a */
1821 38957 : if (signe(b))
1822 : { /* something to check */
1823 : GEN r, t;
1824 13433 : t = dvmdii(mulii(b, v), gel(Q,1), &r);
1825 13433 : if (signe(r)) return mkvec(w);
1826 1820 : u = addii(u, t);
1827 : }
1828 27344 : return mkvec2(w, mkvec2(negi(u), v));
1829 : }
1830 : static GEN
1831 223053 : qfisolvep_all(GEN Q, GEN p, long all)
1832 : {
1833 223053 : GEN R, U, V, M, N, x, q, D = qfb_disc(Q);
1834 223051 : long s = kronecker(D, p);
1835 :
1836 223056 : if (s < 0) return NULL;
1837 126981 : if (!all) s = -1; /* to indicate we want a single solution */
1838 : /* Solutions iff a class of maximal ideal above p is the class of Q;
1839 : * Two solutions iff (s > 0 and the class has order > 2), else one */
1840 126981 : if (!signe(gel(Q,2)))
1841 : { /* if principal form, use faster cornacchia */
1842 43660 : GEN a = gel(Q,1), c = gel(Q,3);
1843 43660 : if (equali1(a))
1844 : {
1845 38183 : if (!cornacchia(c, p, &M,&N)) return NULL;
1846 33725 : return allsols(Q, s, M, N);
1847 : }
1848 5474 : if (equali1(c))
1849 : {
1850 5194 : if (!cornacchia(a, p, &M,&N)) return NULL;
1851 721 : return allsols(Q, s, N, M);
1852 : }
1853 : }
1854 83601 : R = qfi_redsl2_basecase(Q, &U);
1855 83601 : if (equali1(gel(R,1)))
1856 : { /* principal form */
1857 22533 : if (!signe(gel(R,2)))
1858 : {
1859 4396 : if (!cornacchia(gel(R,3), p, &M,&N)) return NULL;
1860 812 : x = mkvec2(M,N);
1861 : }
1862 : else
1863 : { /* x^2 + xy + ((1-D)/4)y^2 = p <==> (2x + y)^2 - D y^2 = 4p */
1864 18137 : if (!cornacchia2(negi(D), p, &M, &N)) return NULL;
1865 2331 : x = subii(M,N); if (mpodd(x)) return NULL;
1866 2331 : x = mkvec2(shifti(x,-1), N);
1867 : }
1868 3143 : x = ZM_ZC_mul(U, x); x[0] = evaltyp(t_VEC) | _evallg(3); /* transpose */
1869 3143 : return allsols(Q, s, gel(x,1), gel(x,2));
1870 : }
1871 61068 : q = qfi_redsl2_basecase(primeform(D, p), &V);
1872 61068 : if (!GL2_qfb_equal(R,q)) return NULL;
1873 10451 : if (signe(gel(R,2)) != signe(gel(q,2))) gcoeff(V,2,1) = negi(gcoeff(V,2,1));
1874 10451 : x = SL2_div_mul_e1(U,V); return allsols(Q, s, gel(x,1), gel(x,2));
1875 : }
1876 : GEN
1877 0 : qfisolvep(GEN Q, GEN p)
1878 : {
1879 0 : pari_sp av = avma;
1880 0 : GEN x = qfisolvep_all(Q, p, 0);
1881 0 : return x? gc_GEN(av, x): gc_const(av, gen_0);
1882 : }
1883 :
1884 : static GEN
1885 770126 : qfrsolve_normform(GEN N, GEN Ps, GEN rd)
1886 : {
1887 770126 : pari_sp av = avma, btop;
1888 770126 : GEN M = N, P = qfr_redsl2_basecase(Ps, rd), Q = P;
1889 :
1890 770126 : btop = avma;
1891 : for(;;)
1892 : {
1893 5840681 : if (qfb_equal(gel(M,1), gel(P,1)))
1894 154084 : return gc_upto(av, SL2_div_mul_e1(gel(M,2),gel(P,2)));
1895 5686597 : if (qfb_equal(gel(N,1), gel(Q,1)))
1896 77658 : return gc_upto(av, SL2_div_mul_e1(gel(N,2),gel(Q,2)));
1897 5608939 : M = qfr_rhosl2(M, rd);
1898 5608939 : if (qfb_equal(gel(M,1), gel(N,1))) return gc_NULL(av);
1899 5201735 : Q = qfr_rhosl2(Q, rd);
1900 5201735 : if (qfb_equal(gel(P,1), gel(Q,1))) return gc_NULL(av);
1901 5070555 : if (gc_needed(btop, 1)) (void)gc_all(btop, 2, &M, &Q);
1902 : }
1903 : }
1904 :
1905 : GEN
1906 0 : qfrsolvep(GEN Q, GEN p)
1907 : {
1908 0 : pari_sp av = avma;
1909 0 : GEN N, x, rd, d = qfb_disc(Q);
1910 :
1911 0 : if (kronecker(d, p) < 0) return gc_const(av, gen_0);
1912 0 : rd = sqrti(d);
1913 0 : N = qfr_redsl2(Q, rd);
1914 0 : x = qfrsolve_normform(N, primeform(d, p), rd);
1915 0 : return x? gc_upto(av, x): gc_const(av, gen_0);
1916 : }
1917 :
1918 : static GEN
1919 1862957 : known_prime(GEN v)
1920 : {
1921 1862957 : GEN p, e, fa = check_arith_all(v, "qfbsolve");
1922 1862953 : if (!fa) return BPSW_psp(v)? v: NULL;
1923 42086 : if (lg(gel(fa,1)) != 2) return NULL;
1924 29361 : p = gcoeff(fa,1,1);
1925 29361 : e = gcoeff(fa,1,2);
1926 29361 : return (equali1(e) && !is_pm1(p) && signe(p) > 0)? p: NULL;
1927 : }
1928 : static GEN
1929 2215793 : qfsolve_normform(GEN Q, GEN f, GEN rd)
1930 2215793 : { return rd? qfrsolve_normform(Q, f, rd): qfisolve_normform(Q, f); }
1931 : static GEN
1932 2843834 : qfbsolve_primitive_i(GEN Q, GEN rd, GEN *Qr, GEN fa, long all)
1933 : {
1934 : GEN x, W, F, p;
1935 : long i, j, l;
1936 2843834 : if (!rd && (p = known_prime(fa))) return qfisolvep_all(Q, p, all);
1937 2620766 : F = normforms(qfb_disc(Q), fa);
1938 2620772 : if (!F) return NULL;
1939 511966 : if (!*Qr) *Qr = qfbredsl2(Q, rd);
1940 511966 : l = lg(F); W = all? cgetg(l, t_VEC): NULL;
1941 2727250 : for (j = i = 1; i < l; i++)
1942 2215793 : if ((x = qfsolve_normform(*Qr, gel(F,i), rd)))
1943 : {
1944 333859 : if (!all) return x;
1945 333348 : gel(W,j++) = x;
1946 : }
1947 511457 : if (j == 1) return NULL;
1948 127458 : setlg(W,j); return lexsort(W);
1949 : }
1950 :
1951 : static GEN
1952 2838530 : qfb_initrd(GEN Q) { GEN d = qfb_disc(Q); return signe(d) > 0? sqrti(d): NULL; }
1953 : static GEN
1954 2828303 : qfbsolve_primitive(GEN Q, GEN fa, long all)
1955 : {
1956 2828303 : GEN x, Qr = NULL, rdQ = qfb_initrd(Q);
1957 2828304 : x = qfbsolve_primitive_i(Q, rdQ, &Qr, fa, all);
1958 2828314 : if (!x) return cgetg(1, t_VEC);
1959 174775 : return x;
1960 : }
1961 :
1962 : /* f / g^2 */
1963 : static GEN
1964 5299 : famat_divsqr(GEN f, GEN g)
1965 5299 : { return famat_reduce(famat_div_shallow(f, famat_pows_shallow(g,2))); }
1966 : static GEN
1967 10227 : qfbsolve_all(GEN Q, GEN n, long all)
1968 : {
1969 10227 : GEN W, Qr = NULL, fa = factorint(n, 0), rdQ = qfb_initrd(Q);
1970 10227 : GEN D = divisors_factored(mkmat2(gel(fa,1), gshift(gel(fa,2),-1)));
1971 10227 : long i, j, l = lg(D);
1972 10227 : W = all? cgetg(l, t_VEC): NULL;
1973 25151 : for (i = j = 1; i < l; i++)
1974 : {
1975 15526 : GEN w, d = gel(D,i), FA = i == 1? fa: famat_divsqr(fa, gel(d,2));
1976 15526 : if ((w = qfbsolve_primitive_i(Q, rdQ, &Qr, FA, all)))
1977 : {
1978 1218 : if (i != 1) w = RgV_Rg_mul(w, gel(d,1));
1979 1218 : if (!all) return w;
1980 616 : gel(W,j++) = w;
1981 : }
1982 : }
1983 9625 : if (j == 1) return cgetg(1, t_VEC);
1984 525 : setlg(W,j); return lexsort(shallowconcat1(W));
1985 : }
1986 :
1987 : GEN
1988 2838533 : qfbsolve(GEN Q, GEN n, long fl)
1989 : {
1990 2838533 : pari_sp av = avma;
1991 2838533 : if (typ(Q) != t_QFB) pari_err_TYPE("qfbsolve",Q);
1992 2838533 : if (fl < 0 || fl > 3) pari_err_FLAG("qfbsolve");
1993 5666842 : return gc_GEN(av, (fl & 2)? qfbsolve_all(Q, n, fl & 1)
1994 2828301 : : qfbsolve_primitive(Q, n, fl & 1));
1995 : }
1996 :
1997 : /* 1 if there exists x,y such that x^2 + dy^2 = p, 0 otherwise;
1998 : * Assume d > 0 and p is prime */
1999 : long
2000 55258 : cornacchia(GEN d, GEN p, GEN *px, GEN *py)
2001 : {
2002 55258 : pari_sp av = avma;
2003 : GEN b, c, r;
2004 :
2005 55258 : *px = *py = gen_0;
2006 55258 : b = subii(p, d);
2007 55210 : if (signe(b) < 0) return gc_long(av,0);
2008 55000 : if (signe(b) == 0) { *py = gen_1; return gc_long(av,1); }
2009 54993 : b = Fp_sqrt(b, p); /* sqrt(-d) */
2010 55051 : if (!b) return gc_long(av,0);
2011 51320 : b = gmael(halfgcdii(p, b), 2, 2);
2012 51364 : c = dvmdii(subii(p, sqri(b)), d, &r);
2013 51249 : if (r != gen_0 || !Z_issquareall(c, &c)) return gc_long(av,0);
2014 35501 : set_avma(av);
2015 35497 : *px = icopy(b);
2016 35489 : *py = icopy(c); return 1;
2017 : }
2018 :
2019 : static GEN
2020 1423667 : lastqi(GEN Q)
2021 : {
2022 1423667 : GEN s = gcoeff(Q,1,1), q = gcoeff(Q,1,2), p = absi_shallow(gcoeff(Q,2,2));
2023 1423664 : if (!signe(q)) return gen_0;
2024 1423475 : if (!signe(s)) return p;
2025 1417091 : if (is_pm1(q)) return subiu(p,1);
2026 1417090 : return divii(p, absi_shallow(q));
2027 : }
2028 :
2029 : static long
2030 1423658 : cornacchia2_i(long av, GEN d, GEN p, GEN b, GEN px4, GEN *px, GEN *py)
2031 : {
2032 : GEN M, Q, V, c, r, b2;
2033 1423658 : if (!signe(b)) { /* d = p,2p,3p,4p */
2034 14 : set_avma(av);
2035 14 : if (absequalii(d, px4)){ *py = gen_1; return 1; }
2036 14 : if (absequalii(d, p)) { *py = gen_2; return 1; }
2037 0 : return 0;
2038 : }
2039 1423644 : if (mod2(b) != mod2(d)) b = subii(p,b);
2040 1423621 : M = halfgcdii(shifti(p,1), b); Q = gel(M,1); V = gel(M, 2);
2041 1423671 : b = addii(mulii(gel(V,1), lastqi(Q)), gel(V,2));
2042 1423565 : b2 = sqri(b);
2043 1423554 : if (cmpii(b2,px4) > 0)
2044 : {
2045 1415180 : b = gel(V,1); b2 = sqri(b);
2046 1415208 : if (cmpii(b2,px4) > 0) { b = gel(V,2); b2 = sqri(b); }
2047 : }
2048 1423564 : c = dvmdii(subii(px4, b2), d, &r);
2049 1423571 : if (r != gen_0 || !Z_issquareall(c, &c)) return gc_long(av,0);
2050 1378148 : set_avma(av);
2051 1378143 : *px = icopy(b);
2052 1378141 : *py = icopy(c); return 1;
2053 : }
2054 :
2055 : /* 1 if there exists x,y such that x^2 + dy^2 = 4p, 0 otherwise;
2056 : * Assume d > 0 is congruent to 0 or 3 mod 4 and p is prime */
2057 : long
2058 1381715 : cornacchia2(GEN d, GEN p, GEN *px, GEN *py)
2059 : {
2060 1381715 : pari_sp av = avma;
2061 1381715 : GEN b, p4 = shifti(p,2);
2062 :
2063 1381666 : *px = *py = gen_0;
2064 1381666 : if (abscmpii(p4, d) < 0) return gc_long(av,0);
2065 1380836 : if (absequaliu(p, 2))
2066 : {
2067 7 : set_avma(av);
2068 7 : switch (itou_or_0(d)) {
2069 0 : case 4: *px = gen_2; break;
2070 0 : case 7: *px = gen_1; break;
2071 7 : default: return 0;
2072 : }
2073 0 : *py = gen_1; return 1;
2074 : }
2075 1380833 : b = Fp_sqrt(negi(d), p);
2076 1380917 : if (!b) return gc_long(av,0);
2077 1380833 : return cornacchia2_i(av, d, p, b, p4, px, py);
2078 : }
2079 :
2080 : /* 1 if there exists x,y such that x^2 + dy^2 = 4p [p prime], 0 otherwise */
2081 : long
2082 42826 : cornacchia2_sqrt(GEN d, GEN p, GEN b, GEN *px, GEN *py)
2083 : {
2084 42826 : pari_sp av = avma;
2085 42826 : GEN p4 = shifti(p,2);
2086 42826 : *px = *py = gen_0;
2087 42826 : if (abscmpii(p4, d) < 0) return gc_long(av,0);
2088 42826 : return cornacchia2_i(av, d, p, b, p4, px, py);
2089 : }
2090 :
2091 : GEN
2092 7630 : qfbcornacchia(GEN d, GEN p)
2093 : {
2094 7630 : pari_sp av = avma;
2095 : GEN x, y;
2096 7630 : if (typ(d) != t_INT || signe(d) <= 0) pari_err_TYPE("qfbcornacchia", d);
2097 7630 : if (typ(p) != t_INT || cmpiu(p, 2) < 0) pari_err_TYPE("qfbcornacchia", p);
2098 7630 : if (mod4(p)? cornacchia(d, p, &x, &y): cornacchia2(d, shifti(p, -2), &x, &y))
2099 287 : return gc_GEN(av, mkvec2(x, y));
2100 7343 : retgc_const(av, cgetg(1, t_VEC));
2101 : }
|
