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 203877 : check_quaddisc(GEN x, long *s, long *pr, const char *f)
24 : {
25 : long r;
26 203877 : if (typ(x) != t_INT) pari_err_TYPE(f,x);
27 203863 : *s = signe(x);
28 203863 : if (Z_issquare(x)) pari_err_DOMAIN(f,"issquare(disc)","=", gen_1,x);
29 203863 : r = mod4(x); if (*s < 0 && r) r = 4 - r;
30 203863 : if (r > 1) pari_err_DOMAIN(f,"disc % 4",">", gen_1,x);
31 203849 : if (pr) *pr = r;
32 203849 : }
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 2184 : check_quaddisc_imag(GEN x, long *r, const char *f)
41 : {
42 2184 : long sx; check_quaddisc(x, &sx, r, f);
43 2177 : if (sx > 0) pari_err_DOMAIN(f, "disc",">",gen_0,x);
44 2177 : }
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 1021000 : quadpoly_bc(GEN D, long Dodd, GEN *b, GEN *c)
50 : {
51 1021000 : if (Dodd)
52 : {
53 921964 : pari_sp av = avma;
54 921964 : *b = gen_m1;
55 921964 : *c = gc_INT(av, shifti(subui(1,D), -2));
56 : }
57 : else
58 : {
59 99036 : *b = gen_0;
60 99036 : *c = shifti(D,-2); togglesign(*c);
61 : }
62 1021000 : }
63 : /* X^2 - X - (D-1)/4 or X^2 - D/4 */
64 : static GEN
65 245665 : quadpoly_ii(GEN D, long Dmod4)
66 : {
67 245665 : GEN b, c, y = cgetg(5,t_POL);
68 245665 : y[1] = evalsigne(1) | evalvarn(0);
69 245665 : quadpoly_bc(D, Dmod4, &b,&c);
70 245665 : gel(y,2) = c;
71 245665 : gel(y,3) = b;
72 245665 : gel(y,4) = gen_1; return y;
73 : }
74 : GEN
75 2205 : quadpoly(GEN D)
76 : {
77 : long s, Dmod4;
78 2205 : check_quaddisc(D, &s, &Dmod4, "quadpoly");
79 2198 : return quadpoly_ii(D, Dmod4);
80 : }
81 : GEN /* no checks */
82 243467 : quadpoly_i(GEN D) { return quadpoly_ii(D, Mod4(D)); }
83 :
84 : GEN
85 1071 : quadpoly0(GEN x, long v)
86 : {
87 1071 : GEN T = quadpoly(x);
88 1064 : if (v > 0) setvarn(T, v);
89 1064 : return T;
90 : }
91 :
92 : GEN
93 0 : quadgen(GEN x)
94 0 : { retmkquad(quadpoly(x), gen_0, gen_1); }
95 :
96 : GEN
97 623 : quadgen0(GEN x, long v)
98 : {
99 623 : if (v==-1) v = fetch_user_var("w");
100 623 : 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 2131138 : check_qfbext(const char *fun, GEN x)
113 : {
114 2131138 : long t = typ(x);
115 2131138 : 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 142504 : qfb3(GEN x, GEN y, GEN z)
127 142504 : { retmkqfb(icopy(x), icopy(y), icopy(z), qfb_disc3(x,y,z)); }
128 :
129 : static int
130 23783634 : qfb_equal(GEN x, GEN y)
131 : {
132 23783634 : return equalii(gel(x,1),gel(y,1))
133 1592913 : && equalii(gel(x,2),gel(y,2))
134 25376547 : && equalii(gel(x,3),gel(y,3));
135 : }
136 :
137 : /* valid for t_QFB, qfr3, qfr5; shallow */
138 : static GEN
139 1042687 : qfb_inv(GEN x) {
140 1042687 : GEN z = shallowcopy(x);
141 1042687 : gel(z,2) = negi(gel(z,2));
142 1042687 : 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 77259 : Qfb0(GEN a, GEN b, GEN c)
151 : {
152 : GEN q, D;
153 77259 : 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 77210 : else if (!c)
169 7 : pari_err_TYPE("Qfb",b);
170 77238 : if (typ(a)!=t_INT) pari_err_TYPE("Qfb",a);
171 77231 : if (typ(b)!=t_INT) pari_err_TYPE("Qfb",b);
172 77231 : if (typ(c)!=t_INT) pari_err_TYPE("Qfb",c);
173 77231 : q = qfb3(a, b, c); D = qfb_disc(q);
174 77231 : if (signe(D) < 0)
175 42392 : { if (signe(a) < 0) pari_err_IMPL("negative definite t_QFB"); }
176 34839 : else if (Z_issquare(D)) pari_err_DOMAIN("Qfb","issquare(disc)","=", gen_1,q);
177 77224 : 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 16933943 : dvmdii_round(GEN b, GEN a, GEN *r)
189 : {
190 16933943 : GEN a2 = shifti(a, 1), q = dvmdii(b, a2, r);
191 16933943 : if (signe(b) >= 0) {
192 9305229 : if (abscmpii(*r, a) > 0) { q = addiu(q, 1); *r = subii(*r, a2); }
193 : } else { /* r <= 0 */
194 7628714 : if (abscmpii(*r, a) >= 0){ q = subiu(q, 1); *r = addii(*r, a2); }
195 : }
196 16933943 : return q;
197 : }
198 : /* Assume 0 < a <= LONG_MAX. Ensure no overflow */
199 : static long
200 120856302 : dvmdsu_round(long b, ulong a, long *r)
201 : {
202 120856302 : ulong a2 = a << 1, q, ub, ur;
203 120856302 : if (b >= 0) {
204 76862111 : ub = b;
205 76862111 : q = ub / a2;
206 76862111 : ur = ub % a2;
207 76862111 : if (ur > a) { ur -= a; q++; *r = (long)ur; *r -= (long)a; }
208 26973712 : else *r = (long)ur;
209 76862111 : return (long)q;
210 : } else { /* r <= 0 */
211 43994191 : ub = (ulong)-b; /* |b| */
212 43994191 : q = ub / a2;
213 43994191 : ur = ub % a2;
214 43994191 : if (ur >= a) { ur -= a; q++; *r = (long)ur; *r = (long)a - *r; }
215 24199455 : else *r = -(long)ur;
216 43994191 : return -(long)q;
217 : }
218 : }
219 : /* reduce b mod 2*a. Update b,c */
220 : static void
221 2776818 : REDB(GEN a, GEN *b, GEN *c)
222 : {
223 2776818 : GEN r, q = dvmdii_round(*b, a, &r);
224 2776818 : if (!signe(q)) return;
225 2707359 : *c = subii(*c, mulii(q, shifti(addii(*b, r),-1)));
226 2707359 : *b = r;
227 : }
228 : /* Assume a > 0. Reduce b mod 2*a. Update b,c */
229 : static void
230 120856302 : sREDB(ulong a, long *b, ulong *c)
231 : {
232 : long r, q;
233 : ulong uz;
234 131244105 : if (a > LONG_MAX) return; /* b already reduced */
235 120856302 : q = dvmdsu_round(*b, a, &r);
236 120856302 : 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 110468499 : if (*b < 0)
240 : { /* uz = -z >= 0, q < 0 */
241 38740464 : if (r >= 0) /* different signs=>no overflow, exact division */
242 19864556 : uz = (ulong)-((*b + r)>>1);
243 : else
244 : {
245 18875908 : ulong ub = (ulong)-*b, ur = (ulong)-r;
246 18875908 : uz = (ub + ur) >> 1;
247 : }
248 38740464 : *c -= (-q) * uz; /* c -= qz */
249 : }
250 : else
251 : { /* uz = z >= 0, q > 0 */
252 71728035 : if (r <= 0)
253 49973388 : uz = (*b + r)>>1;
254 : else
255 : {
256 21754647 : ulong ub = (ulong)*b, ur = (ulong)r;
257 21754647 : uz = ((ub + ur) >> 1);
258 : }
259 71728035 : *c -= q * uz; /* c -= qz */
260 : }
261 110468499 : *b = r;
262 : }
263 : static void
264 14157125 : REDBU(GEN a, GEN *b, GEN *c, GEN u1, GEN *u2)
265 : { /* REDB(a,b,c) */
266 14157125 : GEN r, q = dvmdii_round(*b, a, &r);
267 14157125 : *c = subii(*c, mulii(q, shifti(addii(*b, r),-1)));
268 14157125 : *b = r;
269 14157125 : *u2 = subii(*u2, mulii(q, u1));
270 14157125 : }
271 :
272 : /* q t_QFB, return reduced representative and set base change U in Sl2(Z) */
273 : static GEN
274 6784631 : qfi_redsl2_basecase(GEN q, GEN *U)
275 : {
276 6784631 : pari_sp av = avma;
277 : GEN z, u1,u2,v1,v2,Q;
278 6784631 : GEN a = gel(q,1), b = gel(q,2), c = gel(q,3);
279 : long cmp;
280 6784631 : u1 = gen_1; u2 = gen_0;
281 6784631 : cmp = abscmpii(a, b);
282 6784631 : if (cmp < 0)
283 2198892 : REDBU(a,&b,&c, u1,&u2);
284 4585739 : else if (cmp == 0 && signe(b) < 0)
285 : { /* b = -a */
286 11964 : b = negi(b);
287 11964 : u2 = gen_1;
288 : }
289 : for(;;)
290 : {
291 18742864 : cmp = abscmpii(a, c); if (cmp <= 0) break;
292 11958233 : swap(a,c); b = negi(b);
293 11958233 : z = u1; u1 = u2; u2 = negi(z);
294 11958233 : REDBU(a,&b,&c, u1,&u2);
295 11958233 : 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 6784631 : if (cmp == 0 && signe(b) < 0)
301 : {
302 17739 : b = negi(b);
303 17739 : 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 6784631 : z = shifti(subii(b, gel(q,2)), -1);
308 6784631 : v1 = subii(mulii(z, u1), mulii(a, u2)); v1 = diviiexact(v1, gel(q,3));
309 6784631 : z = subii(z, b);
310 6784631 : v2 = addii(mulii(z, u2), mulii(c, u1)); v2 = diviiexact(v2, gel(q,3));
311 6784631 : *U = mkmat2(mkcol2(u1,v1), mkcol2(u2,v2));
312 6784631 : Q = mkqfb(a,b,c,gel(q,4));
313 6784631 : return gc_all(av, 2, &Q, U);
314 : }
315 :
316 : static GEN
317 1137401 : setq_b0(ulong a, ulong c, GEN D)
318 1137401 : { retmkqfb(utoipos(a), gen_0, utoipos(c), icopy(D)); }
319 : /* assume |sb| = 1 */
320 : static GEN
321 92687076 : setq(ulong a, ulong b, ulong c, long sb, GEN D)
322 92687076 : { 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 982592 : qfi_red_1_b0(ulong a, ulong c, GEN D)
326 982592 : { 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 93964981 : 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 140504 : { /* at most twice */
336 93964981 : long lb = lgefint(b); /* <= 3 after first loop */
337 93964981 : if (lb == 2) return qfi_red_1_b0(a[2],c[2], D);
338 92982389 : if (lb == 3 && uel(b,2) <= (ulong)LONG_MAX) break;
339 140504 : 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 92841885 : set_avma(av);
352 92841885 : ua = a[2];
353 92841885 : ub = sb = b[2]; if (signe(b) < 0) sb = -sb;
354 92841885 : uc = c[2];
355 92841885 : if (ua < ub)
356 35429936 : sREDB(ua, &sb, &uc);
357 57411949 : else if (ua == ub && sb < 0) sb = (long)ub;
358 178268251 : while(ua > uc)
359 : {
360 85426366 : lswap(ua,uc); sb = -sb;
361 85426366 : sREDB(ua, &sb, &uc);
362 : }
363 92841885 : if (!sb) return setq_b0(ua, uc, D);
364 : else
365 : {
366 92687076 : long s = 1;
367 92687076 : if (sb < 0)
368 : {
369 36763503 : sb = -sb;
370 36763503 : if (ua != uc) s = -1;
371 : }
372 92687076 : 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 : /* cf rho_get_BC, with contfrac normalizations: applies rho^(-1) and
435 : * make sure a > 0 */
436 : static GEN
437 175 : rhoi_cf(GEN *A, GEN *B, GEN *C, GEN t)
438 : {
439 175 : GEN u, q, a = *A, b = *B, c = *C;
440 175 : q = truedvmdii(addii(t, b), shifti(a,1), &u);
441 175 : *A = subii(c, mulii(q, subii(b, mulii(q,a))));
442 175 : *B = subii(u, t);
443 175 : *C = a;
444 175 : if (signe(*A) < 0) { *A = negi(*A); *B = negi(*B); *C = negi(*C); }
445 175 : return q;
446 : }
447 : static void
448 14152587 : rho_get_BC(GEN *B, GEN *C, GEN a, GEN b, GEN c, struct qfr_data *S)
449 : {
450 : GEN t, u, q;
451 14152587 : t = (abscmpii(S->isqrtD,c) >= 0)? S->isqrtD: absi_shallow(c);
452 14152587 : q = truedvmdii(addii(t, b), shifti(c,1), &u);
453 14152587 : *B = subii(t, u); /* t - ((t+b) % 2c) */
454 14152587 : *C = subii(a, mulii(q, subii(b, mulii(q,c))));
455 14152587 : }
456 : /* Not stack-clean */
457 : GEN
458 1139460 : qfr3_rho(GEN x, struct qfr_data *S)
459 : {
460 1139460 : GEN B, C, a = gel(x,1), b = gel(x,2), c = gel(x,3);
461 1139460 : rho_get_BC(&B, &C, a, b, c, S);
462 1139460 : return mkvec3(c, B, C);
463 : }
464 :
465 : /* Not stack-clean */
466 : GEN
467 8486681 : qfr5_rho(GEN x, struct qfr_data *S)
468 : {
469 8486681 : GEN B, C, a = gel(x,1), b = gel(x,2), c = gel(x,3), y;
470 8486681 : long sb = signe(b);
471 8486681 : rho_get_BC(&B, &C, a, b, c, S);
472 8486681 : y = mkvec5(c, B, C, gel(x,4), gel(x,5));
473 8486681 : if (sb) {
474 8482698 : GEN t = subii(sqri(b), S->D);
475 8482698 : if (sb < 0)
476 2509500 : t = divir(t, sqrr(subir(b,S->sqrtD)));
477 : else
478 5973198 : t = divri(sqrr(addir(b,S->sqrtD)), t);
479 : /* t = (b + sqrt(D)) / (b - sqrt(D)), evaluated stably */
480 8482698 : gel(y,5) = mulrr(t, gel(y,5)); fix_expo(y);
481 3983 : } else gel(y,5) = negr(gel(y,5));
482 8486681 : return y;
483 : }
484 :
485 : /* Not stack-clean */
486 : GEN
487 217728 : qfr_to_qfr5(GEN x, long prec)
488 217728 : { return mkvec5(gel(x,1),gel(x,2),gel(x,3),gen_0,real_1(prec)); }
489 :
490 : /* d0 = initial distance, x = [a,b,c, expo(d), d], d = exp(2*distance) */
491 : GEN
492 532 : qfr5_to_qfr(GEN x, GEN D, GEN d0)
493 : {
494 532 : if (d0)
495 : {
496 140 : GEN n = gel(x,4), d = absr(gel(x,5));
497 140 : if (signe(n))
498 : {
499 0 : n = addis(shifti(n, EMAX), expo(d));
500 0 : setexpo(d, 0); d = logr_abs(d);
501 0 : if (signe(n)) d = addrr(d, mulir(n, mplog2(lg(d0))));
502 0 : shiftr_inplace(d, -1);
503 0 : d0 = addrr(d0, d);
504 : }
505 140 : else if (!gequal1(d)) /* avoid loss of precision */
506 : {
507 91 : d = logr_abs(d);
508 91 : shiftr_inplace(d, -1);
509 91 : d0 = addrr(d0, d);
510 : }
511 : }
512 532 : x = qfr3_to_qfr(x, D);
513 532 : return d0 ? mkvec2(x,d0): x;
514 : }
515 :
516 : /* Not stack-clean */
517 : GEN
518 31969 : qfr3_to_qfr(GEN x, GEN d) { retmkqfb(gel(x,1), gel(x,2), gel(x,3), d); }
519 :
520 : static int
521 17713057 : ab_isreduced(GEN a, GEN b, GEN isqrtD)
522 : {
523 : GEN t;
524 17713057 : if (signe(b) <= 0 || abscmpii(b, isqrtD) > 0) return 0;
525 5271685 : t = addii_sign(isqrtD,1, shifti(a,1),-1); /* floor(sqrt(D)) - |2a| */
526 1099055 : return signe(t) < 0 ? abscmpii(b, t) >= 0
527 6370740 : : abscmpii(b, t) > 0;
528 : }
529 :
530 : /* Not stack-clean */
531 : GEN
532 1952895 : qfr5_red(GEN x, struct qfr_data *S)
533 : {
534 1952895 : pari_sp av = avma;
535 8465338 : while (!ab_isreduced(gel(x,1), gel(x,2), S->isqrtD))
536 : {
537 6512443 : x = qfr5_rho(x, S);
538 6512443 : if (gc_needed(av,2))
539 : {
540 0 : if (DEBUGMEM>1) pari_warn(warnmem,"qfr5_red");
541 0 : x = gc_GEN(av, x);
542 : }
543 : }
544 1952895 : return x;
545 : }
546 : /* Not stack-clean */
547 : GEN
548 1172882 : qfr3_red(GEN x, struct qfr_data *S)
549 : {
550 1172882 : pari_sp av = avma;
551 1172882 : GEN a = gel(x,1), b = gel(x,2), c = gel(x,3);
552 5699328 : while (!ab_isreduced(a, b, S->isqrtD))
553 : {
554 : GEN B, C;
555 4526446 : rho_get_BC(&B, &C, a, b, c, S);
556 4526446 : a = c; b = B; c = C;
557 4526446 : if (gc_needed(av,2))
558 : {
559 0 : if (DEBUGMEM>1) pari_warn(warnmem,"qfr3_red");
560 0 : (void)gc_all(av, 3, &a, &b, &c);
561 : }
562 : }
563 1172882 : return mkvec3(a, b, c);
564 : }
565 :
566 : void
567 2170 : qfr_data_init(GEN D, long prec, struct qfr_data *S)
568 : {
569 2170 : S->D = D;
570 2170 : S->sqrtD = sqrtr(itor(S->D,prec));
571 2170 : S->isqrtD = truncr(S->sqrtD);
572 2170 : }
573 :
574 : static GEN
575 140 : qfr5_init(GEN x, GEN d, struct qfr_data *S)
576 : {
577 140 : long prec = realprec(d), l = -expo(d);
578 140 : if (l < BITS_IN_LONG) l = BITS_IN_LONG;
579 140 : prec = maxss(prec, nbits2prec(l));
580 140 : S->D = qfb_disc(x);
581 140 : x = qfr_to_qfr5(x,prec);
582 140 : if (!S->sqrtD) S->sqrtD = sqrtr(itor(S->D,prec));
583 0 : else if (typ(S->sqrtD) != t_REAL) pari_err_TYPE("qfr_init",S->sqrtD);
584 :
585 140 : if (!S->isqrtD)
586 : {
587 126 : pari_sp av=avma;
588 : long e;
589 126 : S->isqrtD = gcvtoi(S->sqrtD,&e);
590 126 : if (e>-2) { set_avma(av); S->isqrtD = sqrti(S->D); }
591 : }
592 14 : else if (typ(S->isqrtD) != t_INT) pari_err_TYPE("qfr_init",S->isqrtD);
593 140 : return x;
594 : }
595 : static GEN
596 420 : qfr3_init(GEN x, struct qfr_data *S)
597 : {
598 420 : S->D = qfb_disc(x);
599 420 : if (!S->isqrtD) S->isqrtD = sqrti(S->D);
600 294 : else if (typ(S->isqrtD) != t_INT) pari_err_TYPE("qfr_init",S->isqrtD);
601 420 : return x;
602 : }
603 :
604 : #define qf_NOD 2
605 : #define qf_STEP 1
606 :
607 : static GEN
608 476 : qfr_red_basecase_i(GEN x, long flag, GEN isqrtD, GEN sqrtD)
609 : {
610 : struct qfr_data S;
611 476 : GEN d = NULL, y;
612 476 : if (typ(x)==t_VEC) { d = gel(x,2); x = gel(x,1); } else flag |= qf_NOD;
613 476 : S.sqrtD = sqrtD;
614 476 : S.isqrtD = isqrtD;
615 476 : y = (flag & qf_NOD)? qfr3_init(x, &S): qfr5_init(x, d, &S);
616 476 : switch(flag) {
617 63 : case 0: y = qfr5_red(y,&S); break;
618 371 : case qf_NOD: y = qfr3_red(y,&S); break;
619 21 : case qf_STEP: y = qfr5_rho(y,&S); break;
620 21 : case qf_STEP|qf_NOD: y = qfr3_rho(y,&S); break;
621 0 : default: pari_err_FLAG("qfbred");
622 : }
623 476 : return qfr5_to_qfr(y, qfb_disc(x), d);
624 : }
625 :
626 : static void
627 13379357 : qfr_rhosl2_i(GEN *pa, GEN *pb, GEN *pc, GEN *pu1, GEN *pu2, GEN *pv1,
628 : GEN *pv2, GEN rd)
629 : {
630 13379357 : GEN C = mpabs_shallow(*pc), t = addii(*pb, gmax_shallow(rd,C));
631 13379357 : GEN r, q = truedvmdii(t, shifti(C,1), &r);
632 13379357 : GEN a = *pa, b = *pb, c = *pc;
633 13379357 : if (signe(c) < 0) togglesign(q);
634 13379357 : *pa = *pc;
635 13379357 : *pb = subii(t, addii(r, *pb));
636 13379357 : *pc = subii(a, mulii(q, subii(b, mulii(q,c))));
637 13379357 : r = *pu1; *pu1 = *pv1; *pv1 = subii(mulii(q, *pv1), r);
638 13379357 : r = *pu2; *pu2 = *pv2; *pv2 = subii(mulii(q, *pv2), r);
639 13379357 : }
640 :
641 : static GEN
642 10810674 : qfr_rhosl2(GEN A, GEN rd)
643 : {
644 10810674 : GEN V = gel(A,1), M = gel(A,2);
645 10810674 : GEN a = gel(V,1), b = gel(V,2), c = gel(V,3), d = qfb_disc(V);
646 10810674 : GEN u1 = gcoeff(M,1,1), v1 = gcoeff(M,1,2);
647 10810674 : GEN u2 = gcoeff(M,2,1), v2 = gcoeff(M,2,2);
648 10810674 : qfr_rhosl2_i(&a,&b,&c, &u1,&u2,&v1,&v2, rd);
649 10810674 : return mkvec2(mkqfb(a,b,c,d), mkmat22(u1,v1,u2,v2));
650 : }
651 :
652 : static GEN
653 979701 : qfr_redsl2_basecase(GEN V, GEN rd)
654 : {
655 979701 : pari_sp av = avma;
656 979701 : GEN u1 = gen_1, u2 = gen_0, v1 = gen_0, v2 = gen_1;
657 979701 : GEN a = gel(V,1), b = gel(V,2), c = gel(V,3), d = qfb_disc(V);
658 3548384 : while (!ab_isreduced(a,b,rd))
659 : {
660 2568683 : qfr_rhosl2_i(&a,&b,&c, &u1,&u2,&v1,&v2, rd);
661 2568683 : if (gc_needed(av, 1))
662 : {
663 0 : if (DEBUGMEM>1) pari_warn(warnmem,"qfbredsl2");
664 0 : (void)gc_all(av, 7, &a,&b,&c,&u1,&u2,&v1,&v2);
665 : }
666 : }
667 979701 : return gc_GEN(av, mkvec2(mkqfb(a,b,c,d), mkmat22(u1,v1,u2,v2)));
668 : }
669 :
670 : /* fast reduction of qfb with positive coefficients, based on
671 : Arnold Schoenhage, Fast reduction and composition of binary quadratic forms,
672 : Proc. of Intern. Symp. on Symbolic and Algebraic Computation (Bonn) (S. M.
673 : Watt, ed.), ACM Press, 1991, pp. 128-133.
674 : <https://dl.acm.org/doi/pdf/10.1145/120694.120711>
675 : With thanks to Keegan Ryan
676 : BA20230927
677 : */
678 :
679 : /* pqfb: qf with positive coefficients */
680 :
681 : static int
682 5357816 : lti2n(GEN a, long m) { return signe(a) < 0 || expi(a) < m;}
683 :
684 : static GEN
685 2085300 : pqfbred_1(GEN Q, long m, GEN U)
686 : {
687 2085300 : GEN a = gel(Q,1), b = gel(Q,2), c = gel(Q,3), d = gel(Q,4);
688 2085300 : if (abscmpii(a, c) < 0)
689 : {
690 : GEN t, at, r;
691 1042455 : GEN r2 = addii(shifti(a, m + 2), d);
692 1042455 : long e2 = expi(r2);
693 1042455 : r = int2n(signe(r2) < 0 || e2 < 2*m+2 ? m+1 : e2>>1);
694 1042455 : t = truedivii(subii(b, r), shifti(a,1));
695 1042455 : if (signe(t)==0) pari_err_BUG("pqfbred_1");
696 1042455 : at = mulii(a,t);
697 1042455 : c = addii(subii(c, mulii(b, t)), mulii(at, t));
698 1042455 : b = subii(b, shifti(at,1));
699 1042455 : gcoeff(U,1,2) = subii( gcoeff(U,1,2), mulii(gcoeff(U,1,1), t));
700 1042455 : gcoeff(U,2,2) = subii( gcoeff(U,2,2), mulii(gcoeff(U,2,1), t));
701 : } else
702 : {
703 : GEN t, ct, r;
704 1042845 : GEN r2 = addii(shifti(c, m + 2), d);
705 1042845 : long e2 = expi(r2);
706 1042845 : r = int2n(signe(r2) < 0 || e2 < 2*m+2 ? m+1 : e2>>1);
707 1042845 : t = truedivii(subii(b, r), shifti(c,1));
708 1042845 : if (signe(t)==0) pari_err_BUG("pqfbred_1");
709 1042845 : ct = mulii(c, t);
710 1042845 : a = addii(subii(a, mulii(b, t)), mulii(ct, t));
711 1042845 : b = subii(b, shifti(ct, 1));
712 1042845 : gcoeff(U,1,1) = subii(gcoeff(U,1,1), mulii(gcoeff(U,1,2), t));
713 1042845 : gcoeff(U,2,1) = subii(gcoeff(U,2,1), mulii(gcoeff(U,2,2), t));
714 : }
715 2085300 : return mkqfb(a,b,c,d);
716 : }
717 :
718 : static int
719 2217733 : is_minimal(GEN Q, long m)
720 : {
721 2217733 : pari_sp av = avma;
722 2217733 : GEN a = gel(Q,1), b = gel(Q,2), c = gel(Q,3);
723 5357816 : return gc_bool(av, lti2n(addii(subii(a,b), c), m)
724 2091269 : || (lti2n(subii(b, shifti(a,1)), m+1)
725 1048814 : && lti2n(subii(b, shifti(c,1)), m+1)));
726 : }
727 :
728 : static GEN
729 131244 : pqfbred_iter_1(GEN Q, ulong m, GEN U)
730 : {
731 131244 : pari_sp av = avma;
732 2087424 : while (!is_minimal(Q,m))
733 : {
734 1956180 : Q = pqfbred_1(Q, m, U);
735 1956180 : if (gc_needed(av, 1))
736 : {
737 0 : if (DEBUGMEM>1) pari_warn(warnmem,"pqfbred_iter_1, lc = %ld", expi(gel(Q,3)));
738 0 : (void)gc_all(av, 3, &Q, &gel(U,1), &gel(U,2));
739 : }
740 : }
741 131244 : return Q;
742 : }
743 :
744 : static GEN
745 65097 : pqfbred_basecase(GEN Q, ulong m, GEN *pt_U)
746 : {
747 65097 : pari_sp av = avma;
748 65097 : GEN U = matid(2);
749 65097 : Q = pqfbred_iter_1(Q, m, U);
750 65097 : *pt_U = U;
751 65097 : return gc_all(av, 2, &Q, pt_U);
752 : }
753 :
754 : static long
755 99746689 : qfb_maxexpi(GEN Q)
756 99746689 : { return 1+maxss(expi(gel(Q,1)), maxss(expi(gel(Q,2)), expi(gel(Q,3)))); }
757 :
758 : /* use asymptotically fast reduction ? */
759 : static int
760 99420178 : qfi_red_fast(GEN Q)
761 : {
762 99420178 : const long QFBRED_LIMIT = 9000;
763 99420178 : return 2*qfb_maxexpi(Q) - expi(gel(Q,4)) > QFBRED_LIMIT;
764 : }
765 :
766 : static long
767 132294 : qfb_minexpi(GEN Q)
768 : {
769 132294 : long m = minss(expi(gel(Q,1)), minss(expi(gel(Q,2)), expi(gel(Q,3))));
770 132294 : return m < 0 ? 0: m;
771 : }
772 :
773 : GEN
774 65308 : qfb3_SL2_apply(GEN q, GEN M)
775 : {
776 65308 : GEN a = gel(q,1), b = gel(q,2), c = gel(q,3);
777 65308 : GEN x = gcoeff(M,1,1), y = gcoeff(M,2,1);
778 65308 : GEN z = gcoeff(M,1,2), t = gcoeff(M,2,2);
779 65308 : GEN by = mulii(b,y), bt = mulii(b,t), bz = mulii(b,z);
780 65308 : GEN a2 = shifti(a,1), c2 = shifti(c,1);
781 :
782 65308 : GEN A1 = mulii(x, addii(mulii(a,x), by));
783 65308 : GEN A2 = mulii(c, sqri(y));
784 65308 : GEN B1 = mulii(x, addii(mulii(a2,z), bt));
785 65308 : GEN B2 = mulii(y, addii(mulii(c2,t), bz));
786 65308 : GEN C1 = mulii(z, addii(mulii(a,z), bt));
787 65308 : GEN C2 = mulii(c, sqri(t));
788 65308 : retmkvec3(addii(A1,A2), addii(B1,B2), addii(C1, C2));
789 : }
790 :
791 : static GEN
792 131244 : pqfbred_rec(GEN Q, long m, GEN *pt_U)
793 : {
794 131244 : pari_sp av = avma;
795 131244 : GEN U, Q0, Q1, QR, d = qfb_disc(Q);
796 131244 : long h, n = qfb_maxexpi(Q) - m;
797 131244 : int going_to_r8 = 0;
798 :
799 131244 : if (n < 170) return pqfbred_basecase(Q, m, pt_U);
800 66147 : if (qfb_minexpi(Q) <= m + 2) { U = matid(2); QR = Q; }
801 : else
802 : {
803 : long p, mm;
804 66147 : if (m <= n) { mm = m; p = 0; Q1 = Q; }
805 : else
806 : {
807 65273 : mm = n; p = m + 1 - n;
808 65273 : Q0 = mkvec3(remi2n(gel(Q,1),p), remi2n(gel(Q,2),p), remi2n(gel(Q,3),p));
809 65273 : Q1 = qfb3(shifti(gel(Q,1),-p), shifti(gel(Q,2),-p), shifti(gel(Q,3),-p));
810 : }
811 66147 : h = mm + (n>>1);
812 66147 : if (qfb_minexpi(Q1) <= h) { U = matid(2); QR = Q1; }
813 : else
814 65940 : QR = pqfbred_rec(Q1, h, &U);
815 195267 : while (qfb_maxexpi(QR) > h)
816 : {
817 130309 : if (is_minimal(QR, mm)) { going_to_r8 = 1; break; }
818 129120 : QR = pqfbred_1(QR, mm, U);
819 : }
820 66147 : if (!going_to_r8)
821 : {
822 : GEN V;
823 64958 : QR = pqfbred_rec(QR, mm, &V);
824 64958 : U = ZM2_mul(U,V);
825 : }
826 66147 : if (p > 0)
827 : {
828 65273 : GEN Q0U = qfb3_SL2_apply(Q0,U);
829 130546 : QR = mkqfb(addii(shifti(gel(QR,1), p), gel(Q0U,1)),
830 65273 : addii(shifti(gel(QR,2), p), gel(Q0U,2)),
831 65273 : addii(shifti(gel(QR,3), p), gel(Q0U,3)), d);
832 : }
833 : }
834 66147 : QR = pqfbred_iter_1(QR, m, U);
835 66147 : *pt_U = U; return gc_all(av, 2, &QR, pt_U);
836 : }
837 :
838 : static GEN
839 209575 : qfr_redsl2(GEN Q, GEN isqrtD)
840 : {
841 209575 : pari_sp av = avma;
842 209575 : if (!qfi_red_fast(Q))
843 209575 : return qfr_redsl2_basecase(Q, isqrtD);
844 : else
845 : {
846 0 : GEN a = gel(Q,1), b = gel(Q,2), c = gel(Q,3), d = gel(Q,4);
847 0 : GEN Qf, Qr, W, U, t = NULL;
848 0 : long sa = signe(a), sb;
849 0 : if (sa < 0) { a = negi(a); b = negi(b); c = negi(c); }
850 0 : if (signe(c) < 0)
851 : {
852 : GEN at;
853 0 : t = addiu(truedivii(subii(isqrtD,b),shifti(a,1)),1);
854 0 : at = mulii(a,t);
855 0 : c = addii(subii(c, mulii(b, t)), mulii(at, t));
856 0 : b = subii(b, shifti(at,1));
857 : }
858 0 : sb = signe(b);
859 0 : Qr = pqfbred_rec(mkqfb(a, sb < 0 ? negi(b): b, c, d), 0, &U);
860 0 : if (sa < 0)
861 0 : Qr = mkqfb(negi(gel(Qr,1)), negi(gel(Qr,2)), negi(gel(Qr,3)), gel(Qr,4));
862 0 : if (sb < 0)
863 : {
864 0 : gcoeff(U,2,1) = negi(gcoeff(U,2,1));
865 0 : gcoeff(U,2,2) = negi(gcoeff(U,2,2));
866 : }
867 0 : if (t)
868 : {
869 0 : gcoeff(U,1,1) = subii( gcoeff(U,1,1), mulii(gcoeff(U,2,1), t));
870 0 : gcoeff(U,1,2) = subii( gcoeff(U,1,2), mulii(gcoeff(U,2,2), t));
871 : }
872 0 : W = qfr_redsl2_basecase(Qr, isqrtD);
873 0 : Qf = gel(W,1);
874 0 : U = ZM2_mul(U,gel(W,2));
875 0 : return gc_GEN(av, mkvec2(Qf,U));
876 : }
877 : }
878 :
879 : static GEN
880 5194280 : qfi_redsl2(GEN Q)
881 : {
882 5194280 : pari_sp av = avma;
883 : GEN Qt, U;
884 5194280 : if (!qfi_red_fast(Q))
885 5193969 : Qt = qfi_redsl2_basecase(Q, &U);
886 : else
887 : {
888 311 : long sb = signe(gel(Q,2));
889 : GEN W;
890 311 : if (sb < 0) Q = mkqfb(gel(Q,1), negi(gel(Q,2)), gel(Q,3), gel(Q,4));
891 311 : Q = pqfbred_rec(Q, 0, &U);
892 311 : Qt = qfi_redsl2_basecase(Q, &W);
893 311 : U = ZM2_mul(U,W);
894 311 : if (sb < 0)
895 : {
896 173 : gcoeff(U,2,1) = negi(gcoeff(U,2,1));
897 173 : gcoeff(U,2,2) = negi(gcoeff(U,2,2));
898 : }
899 : }
900 5194280 : return gc_GEN(av, mkvec2(Qt,U));
901 : }
902 :
903 : GEN
904 4883969 : redimagsl2(GEN Q, GEN *U)
905 : {
906 4883969 : GEN q = qfi_redsl2(Q);
907 4883969 : *U = gel(q,2); return gel(q,1);
908 : }
909 :
910 : GEN
911 519893 : qfbredsl2(GEN q, GEN isD)
912 : {
913 : pari_sp av;
914 519893 : if (typ(q) != t_QFB) pari_err_TYPE("qfbredsl2",q);
915 519893 : if (qfb_is_qfi(q))
916 : {
917 310311 : if (isD) pari_err_TYPE("qfbredsl2", isD);
918 310311 : return qfi_redsl2(q);
919 : }
920 209582 : av = avma;
921 209582 : if (!isD) isD = sqrti(qfb_disc(q));
922 208068 : else if (typ(isD) != t_INT) pari_err_TYPE("qfbredsl2",isD);
923 209575 : return gc_upto(av, qfr_redsl2(q, isD));
924 : }
925 :
926 : /* not gc-clean */
927 : static GEN
928 476 : qfr_red_i(GEN Q, long flag, GEN isqrtD, GEN sqrtD)
929 : {
930 476 : if (typ(Q) == t_QFB && !(flag & qf_STEP) && qfi_red_fast(Q))
931 : {
932 28 : GEN U, a = gel(Q,1), b = gel(Q,2), c = gel(Q,3), d = gel(Q,4);
933 28 : long sa = signe(a);
934 28 : if (sa < 0) { a = negi(a); b = negi(b); c = negi(c); }
935 28 : if (signe(c) < 0)
936 : {
937 : GEN at, t;
938 14 : if (!isqrtD) isqrtD = sqrti(d);
939 14 : t = addiu(truedivii(subii(isqrtD,b),shifti(a,1)),1);
940 14 : at = mulii(a,t);
941 14 : c = addii(subii(c, mulii(b, t)), mulii(at, t));
942 14 : b = subii(b, shifti(at,1));
943 : }
944 28 : Q = pqfbred_rec(mkqfb(a, absi_shallow(b), c, d), 0, &U);
945 28 : if (sa < 0)
946 0 : Q = mkqfb(negi(gel(Q,1)), negi(gel(Q,2)), negi(gel(Q,3)), gel(Q,4));
947 : }
948 476 : return qfr_red_basecase_i(Q, flag, isqrtD, sqrtD);
949 : }
950 :
951 : GEN
952 7 : qfr_boundcf(GEN x, long n)
953 : {
954 7 : pari_sp av = avma;
955 7 : GEN a = gel(x,1), b = gel(x,2), c = gel(x,3), d = gel(x,4), V;
956 7 : GEN t = sqrtint(d);
957 : long i;
958 7 : if (signe(a) < 0) { a = negi(a); c = negi(c); } else b = negi(b);
959 7 : V = cgetg(n+1, t_VEC);
960 147 : for (i = 1; i <= n; i++) gel(V,i) = rhoi_cf(&a, &b, &c, t);
961 7 : return gc_GEN(av, V);
962 : }
963 :
964 : GEN
965 7 : qfr_cf(GEN x)
966 : {
967 7 : pari_sp av = avma;
968 7 : GEN a = gel(x,1), b = gel(x,2), c = gel(x,3), d = gel(x,4);
969 7 : GEN t = sqrtint(d);
970 7 : GEN a0 = NULL, b0 = NULL, V, W;
971 7 : long i, l = 16;
972 7 : if (signe(a) < 0) { a = negi(a); c = negi(c); }
973 7 : else b = negi(b);
974 7 : V = cgetg(l+1, t_VEC);
975 7 : for (i = 1;;i++)
976 : {
977 7 : if (ab_isreduced(a, b, t)) break;
978 0 : gel(V,i) = rhoi_cf(&a, &b, &c, t);
979 0 : if (i==l) { l *= 2; V = vec_lengthen(V, l); }
980 : }
981 7 : setlg(V, i); l = 16;
982 7 : a0 = a; b0 = b;
983 7 : W = cgetg(l+1, t_VEC);
984 7 : for (i = 1;; i++)
985 : {
986 35 : gel(W,i) = rhoi_cf(&a, &b, &c, t);
987 35 : if (equalii(a,a0) && equalii(b,b0)) break;
988 28 : if (i==l) { l *= 2; W = vec_lengthen(W, l); }
989 : }
990 7 : setlg(W, i+1); return gc_GEN(av, mkvec2(V,W));
991 : }
992 :
993 : static GEN
994 63 : qfr_red_av(pari_sp av, GEN x)
995 63 : { return gc_GEN(av, qfr_red_i(x,0,NULL,NULL)); }
996 : GEN
997 0 : qfr_red(GEN x) { return qfr_red_av(avma, x); }
998 :
999 : static GEN
1000 94015952 : qfi_red_basecase_av(pari_sp av, GEN q)
1001 : {
1002 94015952 : GEN a = gel(q,1), b = gel(q,2), c = gel(q,3), D = gel(q,4);
1003 94015952 : long cmp, lc = lgefint(c);
1004 :
1005 94015952 : if (lgefint(a) == 3 && lc == 3) return qfi_red_1(av, a, b, c, D);
1006 911922 : cmp = abscmpii(a, b);
1007 911922 : if (cmp < 0)
1008 436234 : REDB(a,&b,&c);
1009 475688 : else if (cmp == 0 && signe(b) < 0)
1010 27 : b = negi(b);
1011 : for(;;)
1012 : {
1013 3112002 : cmp = abscmpii(a, c); if (cmp <= 0) break;
1014 2920527 : lc = lgefint(a); /* lg(future c): we swap a & c next */
1015 2920527 : if (lc == 3) return qfi_red_1(av, a, b, c, D);
1016 2200080 : swap(a,c); b = negi(b); /* apply rho */
1017 2200080 : REDB(a,&b,&c);
1018 2200080 : if (gc_needed(av, 2))
1019 : {
1020 0 : if (DEBUGMEM>1) pari_warn(warnmem,"qfi_red, lc = %ld", lc);
1021 0 : (void)gc_all(av, 3, &a,&b,&c);
1022 : }
1023 : }
1024 191475 : if (cmp == 0 && signe(b) < 0) b = negi(b);
1025 191475 : return gc_GEN(av, mkqfb(a, b, c, D));
1026 : }
1027 : static GEN
1028 94015952 : qfi_red_av(pari_sp av, GEN Q)
1029 : {
1030 94015952 : if (qfi_red_fast(Q))
1031 : {
1032 : GEN U;
1033 7 : if (signe(gel(Q,2)) < 0)
1034 0 : Q = mkqfb(gel(Q,1), negi(gel(Q,2)), gel(Q,3), gel(Q,4));
1035 7 : Q = pqfbred_rec(Q, 0, &U);
1036 : }
1037 94015952 : return qfi_red_basecase_av(av, Q);
1038 : }
1039 :
1040 : GEN
1041 19339781 : qfi_red(GEN q) { return qfi_red_av(avma, q); }
1042 :
1043 : GEN
1044 94076 : qfbred0(GEN x, long flag, GEN isqrtD, GEN sqrtD)
1045 : {
1046 : pari_sp av;
1047 94076 : GEN q = check_qfbext("qfbred",x);
1048 94076 : if (qfb_is_qfi(q)) return (flag & qf_STEP)? qfi_rho(x): qfi_red(x);
1049 413 : if (typ(x)==t_QFB) flag |= qf_NOD;
1050 49 : else flag &= ~qf_NOD;
1051 413 : av = avma;
1052 413 : return gc_GEN(av, qfr_red_i(x,flag,isqrtD,sqrtD));
1053 : }
1054 : /* t_QFB */
1055 : GEN
1056 11761409 : qfbred_i(GEN x) { return qfb_is_qfi(x)? qfi_red(x): qfr_red(x); }
1057 : GEN
1058 92214 : qfbred(GEN x) { return qfbred0(x, 0, NULL, NULL); }
1059 : /***********************************************************************/
1060 : /** **/
1061 : /** Composition **/
1062 : /** **/
1063 : /***********************************************************************/
1064 :
1065 : static void
1066 27236445 : qfb_sqr(GEN z, GEN x)
1067 : {
1068 : GEN c, d1, x2, v1, v2, c3, m, p1, r;
1069 :
1070 27236445 : d1 = bezout(gel(x,2),gel(x,1),&x2, NULL); /* usually 1 */
1071 27236445 : c = gel(x,3);
1072 27236445 : m = mulii(c,x2);
1073 27236445 : if (equali1(d1))
1074 20481639 : v1 = v2 = gel(x,1);
1075 : else
1076 : {
1077 6754806 : v1 = diviiexact(gel(x,1),d1);
1078 6754806 : v2 = mulii(v1, gcdii(d1,c)); /* = v1 iff x primitive */
1079 6754806 : c = mulii(c, d1);
1080 : }
1081 27236445 : togglesign(m);
1082 27236445 : r = modii(m,v2);
1083 27236445 : p1 = mulii(r, v1);
1084 27236445 : c3 = addii(c, mulii(r,addii(gel(x,2),p1)));
1085 27236445 : gel(z,1) = mulii(v1,v2);
1086 27236445 : gel(z,2) = addii(gel(x,2), shifti(p1,1));
1087 27236445 : gel(z,3) = diviiexact(c3,v2);
1088 27236445 : }
1089 : /* z <- x * y */
1090 : static void
1091 76187759 : qfb_comp(GEN z, GEN x, GEN y)
1092 : {
1093 : GEN n, c, d, y1, v1, v2, c3, m, p1, r;
1094 :
1095 76187759 : if (x == y) { qfb_sqr(z,x); return; }
1096 49774279 : n = shifti(subii(gel(y,2),gel(x,2)), -1);
1097 49774279 : v1 = gel(x,1);
1098 49774279 : v2 = gel(y,1);
1099 49774279 : c = gel(y,3);
1100 49774279 : d = bezout(v2,v1,&y1,NULL);
1101 49774279 : if (equali1(d))
1102 30861456 : m = mulii(y1,n);
1103 : else
1104 : {
1105 18912823 : GEN s = subii(gel(y,2), n);
1106 18912823 : GEN x2, y2, d1 = bezout(s,d,&x2,&y2); /* x2 s + y2 (x1 v1 + y1 v2) = d1 */
1107 18912823 : if (!equali1(d1))
1108 : {
1109 9049670 : v1 = diviiexact(v1,d1);
1110 9049670 : v2 = diviiexact(v2,d1); /* gcd = 1 iff x or y primitive */
1111 9049670 : v1 = mulii(v1, gcdii(c,gcdii(gel(x,3),gcdii(d1,n))));
1112 9049670 : c = mulii(c, d1);
1113 : }
1114 18912823 : m = addii(mulii(mulii(y1,y2),n), mulii(gel(y,3),x2));
1115 : }
1116 49774279 : togglesign(m);
1117 49774279 : r = modii(m, v1);
1118 49774279 : p1 = mulii(r, v2);
1119 49774279 : c3 = addii(c, mulii(r,addii(gel(y,2),p1)));
1120 49774279 : gel(z,1) = mulii(v1,v2);
1121 49774279 : gel(z,2) = addii(gel(y,2), shifti(p1,1));
1122 49774279 : gel(z,3) = diviiexact(c3,v1);
1123 : }
1124 :
1125 : /* not meant to be efficient */
1126 : static GEN
1127 84 : qfb_comp_gen(GEN x, GEN y)
1128 : {
1129 84 : GEN d1 = qfb_disc(x), d2 = qfb_disc(y);
1130 84 : GEN a1 = gel(x,1), b1 = gel(x,2), c1 = gel(x,3), n1;
1131 84 : GEN a2 = gel(y,1), b2 = gel(y,2), c2 = gel(y,3), n2;
1132 84 : GEN cx = content(x), cy = content(y), A, B, C, D, U, m, m2;
1133 :
1134 84 : if (!is_pm1(cx))
1135 : {
1136 14 : a1 = diviiexact(a1, cx); b1 = diviiexact(b1, cx);
1137 14 : c1 = diviiexact(c1, cx); d1 = diviiexact(d1, sqri(cx));
1138 : }
1139 84 : if (!is_pm1(cy))
1140 : {
1141 28 : a2 = diviiexact(a2, cy); c2 = diviiexact(c2, cy);
1142 28 : b2 = diviiexact(b2, cy); d2 = diviiexact(d2, sqri(cy));
1143 : }
1144 84 : D = gcdii(d1, d2); if (signe(d1) < 0) setsigne(D, -1);
1145 133 : if (!Z_issquareall(diviiexact(d1, D), &n1) ||
1146 84 : !Z_issquareall(diviiexact(d2, D), &n2)) return NULL;
1147 49 : A = mulii(a1, n2);
1148 49 : B = mulii(a2, n1);
1149 49 : C = shifti(addii(mulii(b1, n2), mulii(b2, n1)), -1);
1150 49 : U = ZV_extgcd(mkvec3(A, B, C));
1151 49 : m = gel(U,1); U = gmael(U,2,3);
1152 49 : A = mulii(diviiexact(mulii(a1,b2),m), gel(U,1));
1153 49 : B = mulii(diviiexact(mulii(a2,b1),m), gel(U,2));
1154 49 : C = addii(mulii(b1,b2), mulii(D, mulii(n1,n2)));
1155 49 : C = mulii(diviiexact(shifti(C,-1), m), gel(U,3));
1156 49 : B = addii(A, addii(B, C));
1157 49 : m2 = sqri(m);
1158 49 : A = diviiexact(mulii(a1, a2), m2);
1159 49 : C = diviiexact(shifti(subii(sqri(B),D), -2), A);
1160 49 : cx = mulii(cx, cy);
1161 49 : if (!is_pm1(cx))
1162 : {
1163 14 : A = mulii(A, cx); B = mulii(B, cx);
1164 14 : C = mulii(C, cx); D = mulii(D, sqri(cx));
1165 : }
1166 49 : return mkqfb(A, B, C, D);
1167 : }
1168 :
1169 : static GEN
1170 73442691 : qficomp0(GEN x, GEN y, int raw)
1171 : {
1172 73442691 : pari_sp av = avma;
1173 73442691 : GEN z = cgetg(5,t_QFB);
1174 73442691 : gel(z,4) = gel(x,4);
1175 73442691 : qfb_comp(z, x,y);
1176 73442691 : if (raw) return gc_GEN(av,z);
1177 73440913 : return qfi_red_av(av, z);
1178 : }
1179 : static GEN
1180 441 : qfrcomp0(GEN x, GEN y, int raw)
1181 : {
1182 441 : pari_sp av = avma;
1183 441 : GEN dx = NULL, dy = NULL;
1184 441 : GEN z = cgetg(5,t_QFB);
1185 441 : if (typ(x)==t_VEC) { dx = gel(x,2); x = gel(x,1); }
1186 441 : if (typ(y)==t_VEC) { dy = gel(y,2); y = gel(y,1); }
1187 441 : gel(z,4) = gel(x,4);
1188 441 : qfb_comp(z, x,y);
1189 441 : if (dx) z = mkvec2(z, dy? addrr(dx, dy): dx); else if (dy) z = mkvec2(z, dy);
1190 441 : if (raw) return gc_GEN(av, z);
1191 28 : return qfr_red_av(av, z);
1192 : }
1193 : /* same discriminant, no distance, no checks */
1194 : GEN
1195 38663194 : qfbcomp_i(GEN x, GEN y)
1196 38663194 : { return qfb_is_qfi(x)? qficomp0(x,y,0): qfrcomp0(x,y,0); }
1197 : GEN
1198 138644 : qfbcomp(GEN x, GEN y)
1199 : {
1200 138644 : GEN qx = check_qfbext("qfbcomp", x);
1201 138644 : GEN qy = check_qfbext("qfbcomp", y);
1202 138644 : if (!equalii(gel(qx,4),gel(qy,4)))
1203 : {
1204 63 : pari_sp av = avma;
1205 63 : GEN z = qfb_comp_gen(qx, qy);
1206 63 : if (typ(x) == t_VEC || typ(y) == t_VEC)
1207 7 : pari_err_IMPL("Shanks's distance in general composition");
1208 56 : if (!z) pari_err_OP("*",x,y);
1209 21 : return gc_upto(av, qfbred(z));
1210 : }
1211 138581 : return qfb_is_qfi(qx)? qficomp0(x,y,0): qfrcomp0(x,y,0);
1212 : }
1213 : /* same discriminant, no distance, no checks */
1214 : GEN
1215 0 : qfbcompraw_i(GEN x, GEN y)
1216 0 : { return qfb_is_qfi(x)? qficomp0(x,y,1): qfrcomp0(x,y,1); }
1217 : GEN
1218 2198 : qfbcompraw(GEN x, GEN y)
1219 : {
1220 2198 : GEN qx = check_qfbext("qfbcompraw", x);
1221 2198 : GEN qy = check_qfbext("qfbcompraw", y);
1222 2198 : if (!equalii(gel(qx,4),gel(qy,4)))
1223 : {
1224 21 : pari_sp av = avma;
1225 21 : GEN z = qfb_comp_gen(qx, qy);
1226 21 : if (typ(x) == t_VEC || typ(y) == t_VEC)
1227 0 : pari_err_IMPL("Shanks's distance in general composition");
1228 21 : if (!z) pari_err_OP("qfbcompraw",x,y);
1229 21 : return gc_GEN(av, z);
1230 : }
1231 2177 : if (!equalii(gel(qx,4),gel(qy,4))) pari_err_OP("qfbcompraw",x,y);
1232 2177 : return qfb_is_qfi(qx)? qficomp0(x,y,1): qfrcomp0(x,y,1);
1233 : }
1234 :
1235 : static GEN
1236 822930 : qfisqr0(GEN x, long raw)
1237 : {
1238 822930 : pari_sp av = avma;
1239 822930 : GEN z = cgetg(5,t_QFB);
1240 822930 : gel(z,4) = gel(x,4);
1241 822930 : qfb_sqr(z,x);
1242 822930 : if (raw) return gc_GEN(av,z);
1243 822930 : return qfi_red_av(av, z);
1244 : }
1245 : static GEN
1246 35 : qfrsqr0(GEN x, long raw)
1247 : {
1248 35 : pari_sp av = avma;
1249 35 : GEN dx = NULL, z = cgetg(5,t_QFB);
1250 35 : if (typ(x) == t_VEC) { dx = gel(x,2); x = gel(x,1); }
1251 35 : gel(z,4) = gel(x,4); qfb_sqr(z,x);
1252 35 : if (dx) z = mkvec2(z, shiftr(dx,1));
1253 35 : if (raw) return gc_GEN(av, z);
1254 35 : return qfr_red_av(av, z);
1255 : }
1256 : /* same discriminant, no distance, no checks */
1257 : GEN
1258 694003 : qfbsqr_i(GEN x)
1259 694003 : { return qfb_is_qfi(x)? qfisqr0(x,0): qfrsqr0(x,0); }
1260 : GEN
1261 128962 : qfbsqr(GEN x)
1262 : {
1263 128962 : GEN qx = check_qfbext("qfbsqr", x);
1264 128962 : return qfb_is_qfi(qx)? qfisqr0(x,0): qfrsqr0(x,0);
1265 : }
1266 :
1267 : static GEN
1268 6867 : qfr_1_by_disc(GEN D)
1269 : {
1270 : GEN y, r, s;
1271 6867 : check_quaddisc_real(D, NULL, "qfr_1_by_disc");
1272 6867 : y = cgetg(5,t_QFB);
1273 6867 : s = sqrtremi(D, &r); togglesign(r); /* s^2 - r = D */
1274 6867 : if (mpodd(r))
1275 : {
1276 3535 : s = subiu(s,1);
1277 3535 : r = subii(r, addiu(shifti(s, 1), 1));
1278 3535 : r = shifti(r, -2); set_avma((pari_sp)y); s = icopy(s);
1279 : }
1280 : else
1281 3332 : { r = shifti(r, -2); set_avma((pari_sp)s); }
1282 6867 : gel(y,1) = gen_1;
1283 6867 : gel(y,2) = s;
1284 6867 : gel(y,3) = icopy(r);
1285 6867 : gel(y,4) = icopy(D); return y;
1286 : }
1287 :
1288 : static GEN
1289 35 : qfr_disc(GEN x)
1290 35 : { return qfb_disc(typ(x)==t_VEC ? gel(x,1): x); }
1291 :
1292 : static GEN
1293 35 : qfr_1(GEN x)
1294 35 : { return qfr_1_by_disc(qfr_disc(x)); }
1295 :
1296 : static void
1297 0 : qfr_1_fill(GEN y, struct qfr_data *S)
1298 : {
1299 0 : pari_sp av = avma;
1300 0 : GEN y2 = S->isqrtD;
1301 0 : gel(y,1) = gen_1;
1302 0 : if (mod2(S->D) != mod2(y2)) y2 = subiu(y,1);
1303 0 : gel(y,2) = y2; av = avma;
1304 0 : gel(y,3) = gc_INT(av, shifti(subii(sqri(y2), S->D),-2));
1305 0 : }
1306 : static GEN
1307 0 : qfr5_1(struct qfr_data *S, long prec)
1308 : {
1309 0 : GEN y = cgetg(6, t_VEC);
1310 0 : qfr_1_fill(y, S);
1311 0 : gel(y,4) = gen_0;
1312 0 : gel(y,5) = real_1(prec); return y;
1313 : }
1314 : static GEN
1315 0 : qfr3_1(struct qfr_data *S)
1316 : {
1317 0 : GEN y = cgetg(4, t_VEC);
1318 0 : qfr_1_fill(y, S); return y;
1319 : }
1320 :
1321 : /* Assume D < 0 is the discriminant of a t_QFB */
1322 : static GEN
1323 775335 : qfi_1_by_disc(GEN D)
1324 : {
1325 775335 : GEN b,c, y = cgetg(5,t_QFB);
1326 775335 : quadpoly_bc(D, mod2(D), &b,&c);
1327 775335 : if (b == gen_m1) b = gen_1;
1328 775335 : gel(y,1) = gen_1;
1329 775335 : gel(y,2) = b;
1330 775335 : gel(y,3) = c;
1331 775335 : gel(y,4) = icopy(D); return y;
1332 : }
1333 : static GEN
1334 763241 : qfi_1(GEN x)
1335 : {
1336 763241 : if (typ(x) != t_QFB) pari_err_TYPE("qfi_1",x);
1337 763241 : return qfi_1_by_disc(qfb_disc(x));
1338 : }
1339 :
1340 : GEN
1341 0 : qfb_1(GEN x) { return qfb_is_qfi(x) ? qfi_1(x): qfr_1(x); }
1342 :
1343 : static GEN
1344 9607916 : _qfimul(void *E, GEN x, GEN y) { (void) E; return qficomp0(x,y,0); }
1345 : static GEN
1346 25031250 : _qfisqr(void *E, GEN x) { (void) E; return qficomp0(x,x,0); }
1347 : static GEN
1348 7 : _qfimulraw(void *E, GEN x, GEN y) { (void) E; return qficomp0(x,y,1); }
1349 : static GEN
1350 7 : _qfisqrraw(void *E, GEN x) { (void) E; return qficomp0(x,x,1); }
1351 :
1352 : static GEN
1353 7 : qfipowraw(GEN x, long n)
1354 : {
1355 7 : pari_sp av = avma;
1356 : GEN y;
1357 7 : if (!n) return qfi_1(x);
1358 7 : if (n== 1) return gcopy(x);
1359 7 : if (n==-1) { x = gcopy(x); togglesign(gel(x,2)); return x; }
1360 7 : if (n < 0) x = qfb_inv(x);
1361 7 : y = gen_powu(x, labs(n), NULL, &_qfisqrraw, &_qfimulraw);
1362 7 : return gc_GEN(av,y);
1363 : }
1364 :
1365 : static GEN
1366 12524650 : qfipow(GEN x, GEN n)
1367 : {
1368 12524650 : pari_sp av = avma;
1369 : GEN y;
1370 12524650 : long s = signe(n);
1371 12524650 : if (!s) return qfi_1(x);
1372 11761409 : if (s < 0) x = qfb_inv(x);
1373 11761409 : y = gen_pow(qfbred_i(x), n, NULL, &_qfisqr, &_qfimul);
1374 11761409 : return gc_GEN(av,y);
1375 : }
1376 :
1377 : static long
1378 412328 : parteucl(GEN L, GEN *d, GEN *v3, GEN *v, GEN *v2)
1379 : {
1380 : long z;
1381 412328 : *v = gen_0; *v2 = gen_1;
1382 4351417 : for (z=0; abscmpii(*v3,L) > 0; z++)
1383 : {
1384 3939089 : GEN t3, t2 = subii(*v, mulii(truedvmdii(*d,*v3,&t3),*v2));
1385 3939089 : *v = *v2; *d = *v3; *v2 = t2; *v3 = t3;
1386 : }
1387 412328 : return z;
1388 : }
1389 :
1390 : /* composition: Shanks' NUCOMP & NUDUPL */
1391 : /* L = floor((|d|/4)^(1/4)) */
1392 : GEN
1393 400722 : nucomp(GEN x, GEN y, GEN L)
1394 : {
1395 400722 : pari_sp av = avma;
1396 : long z;
1397 : GEN a, a1, a2, b2, b, d, d1, g, n, p1, q1, q2, s, u, u1, v, v2, v3, Q;
1398 :
1399 400722 : if (x==y) return nudupl(x,L);
1400 400680 : if (!is_qfi(x)) pari_err_TYPE("nucomp",x);
1401 400680 : if (!is_qfi(y)) pari_err_TYPE("nucomp",y);
1402 :
1403 400680 : if (abscmpii(gel(x,1),gel(y,1)) < 0) swap(x, y);
1404 400680 : s = shifti(addii(gel(x,2),gel(y,2)), -1);
1405 400680 : n = subii(gel(y,2), s);
1406 400680 : a1 = gel(x,1);
1407 400680 : a2 = gel(y,1); d = bezout(a2,a1,&u,&v);
1408 400680 : if (equali1(d)) { a = negi(mulii(u,n)); d1 = d; }
1409 163576 : else if (dvdii(s,d)) /* d | s */
1410 : {
1411 83503 : a = negi(mulii(u,n)); d1 = d;
1412 83503 : a1 = diviiexact(a1, d1);
1413 83503 : a2 = diviiexact(a2, d1);
1414 83503 : s = diviiexact(s, d1);
1415 : }
1416 : else
1417 : {
1418 : GEN p2, l;
1419 80073 : d1 = bezout(s,d,&u1,NULL);
1420 80073 : if (!equali1(d1))
1421 : {
1422 2044 : a1 = diviiexact(a1,d1);
1423 2044 : a2 = diviiexact(a2,d1);
1424 2044 : s = diviiexact(s,d1);
1425 2044 : d = diviiexact(d,d1);
1426 : }
1427 80073 : p1 = remii(gel(x,3),d);
1428 80073 : p2 = remii(gel(y,3),d);
1429 80073 : l = modii(mulii(negi(u1), addii(mulii(u,p1),mulii(v,p2))), d);
1430 80073 : a = subii(mulii(l,diviiexact(a1,d)), mulii(u,diviiexact(n,d)));
1431 : }
1432 400680 : a = modii(a,a1); p1 = subii(a,a1); if (abscmpii(a,p1) > 0) a = p1;
1433 400680 : d = a1; v3 = a; z = parteucl(L, &d,&v3, &v,&v2);
1434 400680 : Q = cgetg(5,t_QFB);
1435 400680 : if (!z) {
1436 37632 : g = diviiexact(addii(mulii(v3,s),gel(y,3)), d);
1437 37632 : b = a2;
1438 37632 : b2 = gel(y,2);
1439 37632 : v2 = d1;
1440 37632 : gel(Q,1) = mulii(d,b);
1441 : } else {
1442 : GEN e, q3, q4;
1443 363048 : if (z&1) { v3 = negi(v3); v2 = negi(v2); }
1444 363048 : b = diviiexact(addii(mulii(a2,d), mulii(n,v)), a1);
1445 363048 : e = diviiexact(addii(mulii(s,d),mulii(gel(y,3),v)), a1);
1446 363048 : q3 = mulii(e,v2);
1447 363048 : q4 = subii(q3,s);
1448 363048 : b2 = addii(q3,q4);
1449 363048 : g = diviiexact(q4,v);
1450 363048 : if (!equali1(d1)) { v2 = mulii(d1,v2); v = mulii(d1,v); b2 = mulii(d1,b2); }
1451 363048 : gel(Q,1) = addii(mulii(d,b), mulii(e,v));
1452 : }
1453 400680 : q1 = mulii(b, v3);
1454 400680 : q2 = addii(q1,n);
1455 400680 : gel(Q,2) = addii(b2, z? addii(q1,q2): shifti(q1, 1));
1456 400680 : gel(Q,3) = addii(mulii(v3,diviiexact(q2,d)), mulii(g,v2));
1457 400680 : gel(Q,4) = gel(x,4);
1458 400680 : return qfi_red_av(av, Q);
1459 : }
1460 :
1461 : GEN
1462 11648 : nudupl(GEN x, GEN L)
1463 : {
1464 11648 : pari_sp av = avma;
1465 : long z;
1466 : GEN u, v, d, d1, p1, a, b, c, a2, b2, c2, Q, v2, v3, g;
1467 :
1468 11648 : if (!is_qfi(x)) pari_err_TYPE("nudupl",x);
1469 11648 : a = gel(x,1);
1470 11648 : b = gel(x,2);
1471 11648 : d1 = bezout(b,a, &u,NULL);
1472 11648 : if (!equali1(d1))
1473 : {
1474 4620 : a = diviiexact(a, d1);
1475 4620 : b = diviiexact(b, d1);
1476 : }
1477 11648 : c = modii(negi(mulii(u,gel(x,3))), a);
1478 11648 : p1 = subii(c,a); if (abscmpii(c,p1) > 0) c = p1;
1479 11648 : d = a; v3 = c; z = parteucl(L, &d,&v3, &v,&v2);
1480 11648 : a2 = sqri(d);
1481 11648 : c2 = sqri(v3);
1482 11648 : Q = cgetg(5,t_QFB);
1483 11648 : if (!z) {
1484 1281 : g = diviiexact(addii(mulii(v3,b),gel(x,3)), d);
1485 1281 : b2 = gel(x,2);
1486 1281 : v2 = d1;
1487 1281 : gel(Q,1) = a2;
1488 : } else {
1489 : GEN e;
1490 10367 : if (z&1) { v = negi(v); d = negi(d); }
1491 10367 : e = diviiexact(addii(mulii(gel(x,3),v), mulii(b,d)), a);
1492 10367 : g = diviiexact(subii(mulii(e,v2), b), v);
1493 10367 : b2 = addii(mulii(e,v2), mulii(v,g));
1494 10367 : if (!equali1(d1)) { b2 = mulii(d1,b2); v = mulii(d1,v); v2 = mulii(d1,v2); }
1495 10367 : gel(Q,1) = addii(a2, mulii(e,v));
1496 : }
1497 11648 : gel(Q,2) = addii(b2, subii(sqri(addii(d,v3)), addii(a2,c2)));
1498 11648 : gel(Q,3) = addii(c2, mulii(g,v2));
1499 11648 : gel(Q,4) = gel(x,4);
1500 11648 : return qfi_red_av(av, Q);
1501 : }
1502 :
1503 : static GEN
1504 4739 : mul_nucomp(void *l, GEN x, GEN y) { return nucomp(x, y, (GEN)l); }
1505 : static GEN
1506 11606 : mul_nudupl(void *l, GEN x) { return nudupl(x, (GEN)l); }
1507 : GEN
1508 1008 : nupow(GEN x, GEN n, GEN L)
1509 : {
1510 : pari_sp av;
1511 : GEN y, D;
1512 :
1513 1008 : if (typ(n) != t_INT) pari_err_TYPE("nupow",n);
1514 1008 : if (!is_qfi(x)) pari_err_TYPE("nupow",x);
1515 1008 : if (gequal1(n)) return gcopy(x);
1516 1008 : av = avma;
1517 1008 : D = qfb_disc(x);
1518 1008 : y = qfi_1_by_disc(D);
1519 1008 : if (!signe(n)) return y;
1520 959 : if (!L) L = sqrtnint(absi_shallow(D), 4);
1521 959 : y = gen_pow_i(x, n, (void*)L, &mul_nudupl, &mul_nucomp);
1522 959 : if (signe(n) < 0
1523 35 : && !absequalii(gel(y,1),gel(y,2))
1524 35 : && !absequalii(gel(y,1),gel(y,3))) togglesign(gel(y,2));
1525 959 : return gc_GEN(av, y);
1526 : }
1527 :
1528 : /* Not stack-clean */
1529 : GEN
1530 1735230 : qfr5_compraw(GEN x, GEN y)
1531 : {
1532 1735230 : GEN z = cgetg(6,t_VEC); qfb_comp(z,x,y);
1533 1735230 : if (x == y)
1534 : {
1535 34552 : gel(z,4) = shifti(gel(x,4),1);
1536 34552 : gel(z,5) = sqrr(gel(x,5));
1537 : }
1538 : else
1539 : {
1540 1700678 : gel(z,4) = addii(gel(x,4),gel(y,4));
1541 1700678 : gel(z,5) = mulrr(gel(x,5),gel(y,5));
1542 : }
1543 1735230 : fix_expo(z); return z;
1544 : }
1545 : GEN
1546 1735216 : qfr5_comp(GEN x, GEN y, struct qfr_data *S)
1547 1735216 : { return qfr5_red(qfr5_compraw(x, y), S); }
1548 : /* Not stack-clean */
1549 : GEN
1550 1009397 : qfr3_compraw(GEN x, GEN y)
1551 : {
1552 1009397 : GEN z = cgetg(4,t_VEC); qfb_comp(z,x,y);
1553 1009397 : return z;
1554 : }
1555 : GEN
1556 1009397 : qfr3_comp(GEN x, GEN y, struct qfr_data *S)
1557 1009397 : { return qfr3_red(qfr3_compraw(x,y), S); }
1558 :
1559 : /* m > 0. Not stack-clean */
1560 : static GEN
1561 7 : qfr5_powraw(GEN x, long m)
1562 : {
1563 7 : GEN y = NULL;
1564 14 : for (; m; m >>= 1)
1565 : {
1566 14 : if (m&1) y = y? qfr5_compraw(y,x): x;
1567 14 : if (m == 1) break;
1568 7 : x = qfr5_compraw(x,x);
1569 : }
1570 7 : return y;
1571 : }
1572 :
1573 : /* return x^n. Not stack-clean */
1574 : GEN
1575 21 : qfr5_pow(GEN x, GEN n, struct qfr_data *S)
1576 : {
1577 21 : GEN y = NULL;
1578 21 : long i, m, s = signe(n);
1579 21 : if (!s) return qfr5_1(S, lg(gel(x,5)));
1580 21 : if (s < 0) x = qfb_inv(x);
1581 42 : for (i=lgefint(n)-1; i>1; i--)
1582 : {
1583 21 : m = n[i];
1584 56 : for (; m; m>>=1)
1585 : {
1586 56 : if (m&1) y = y? qfr5_comp(y,x,S): x;
1587 56 : if (m == 1 && i == 2) break;
1588 35 : x = qfr5_comp(x,x,S);
1589 : }
1590 : }
1591 21 : return y;
1592 : }
1593 : /* m > 0; return x^m. Not stack-clean */
1594 : static GEN
1595 0 : qfr3_powraw(GEN x, long m)
1596 : {
1597 0 : GEN y = NULL;
1598 0 : for (; m; m>>=1)
1599 : {
1600 0 : if (m&1) y = y? qfr3_compraw(y,x): x;
1601 0 : if (m == 1) break;
1602 0 : x = qfr3_compraw(x,x);
1603 : }
1604 0 : return y;
1605 : }
1606 : /* return x^n. Not stack-clean */
1607 : GEN
1608 4557 : qfr3_pow(GEN x, GEN n, struct qfr_data *S)
1609 : {
1610 4557 : GEN y = NULL;
1611 4557 : long i, m, s = signe(n);
1612 4557 : if (!s) return qfr3_1(S);
1613 4557 : if (s < 0) x = qfb_inv(x);
1614 9130 : for (i=lgefint(n)-1; i>1; i--)
1615 : {
1616 4573 : m = n[i];
1617 5312 : for (; m; m>>=1)
1618 : {
1619 5292 : if (m&1) y = y? qfr3_comp(y,x,S): x;
1620 5292 : if (m == 1 && i == 2) break;
1621 739 : x = qfr3_comp(x,x,S);
1622 : }
1623 : }
1624 4557 : return y;
1625 : }
1626 :
1627 : static GEN
1628 7 : qfrinvraw(GEN x)
1629 : {
1630 7 : if (typ(x) == t_VEC) retmkvec2(qfbinv(gel(x,1)), negr(gel(x,2)));
1631 7 : return qfbinv(x);
1632 : }
1633 : static GEN
1634 14 : qfrpowraw(GEN x, long n)
1635 : {
1636 14 : struct qfr_data S = { NULL, NULL, NULL };
1637 14 : pari_sp av = avma;
1638 14 : if (n==1) return gcopy(x);
1639 14 : if (n==-1) return qfrinvraw(x);
1640 7 : if (typ(x)==t_QFB)
1641 : {
1642 0 : GEN D = qfb_disc(x);
1643 0 : if (!n) return qfr_1(x);
1644 0 : if (n < 0) { x = qfb_inv(x); n = -n; }
1645 0 : x = qfr3_powraw(x, n);
1646 0 : x = qfr3_to_qfr(x, D);
1647 : }
1648 : else
1649 : {
1650 7 : GEN d0 = gel(x,2);
1651 7 : x = gel(x,1);
1652 7 : if (!n) retmkvec2(qfr_1(x), real_0(precision(d0)));
1653 7 : if (n < 0) { x = qfb_inv(x); n = -n; }
1654 7 : x = qfr5_init(x, d0, &S);
1655 7 : if (labs(n) != 1) x = qfr5_powraw(x, n);
1656 7 : x = qfr5_to_qfr(x, S.D, mulrs(d0,n));
1657 : }
1658 7 : return gc_GEN(av, x);
1659 : }
1660 : static GEN
1661 112 : qfrpow(GEN x, GEN n)
1662 : {
1663 112 : struct qfr_data S = { NULL, NULL, NULL };
1664 112 : long s = signe(n);
1665 112 : pari_sp av = avma;
1666 112 : if (typ(x)==t_QFB)
1667 : {
1668 42 : if (!s) return qfr_1(x);
1669 28 : if (s < 0) x = qfb_inv(x);
1670 28 : x = qfr3_init(x, &S);
1671 28 : x = is_pm1(n)? qfr3_red(x, &S): qfr3_pow(x, n, &S);
1672 28 : x = qfr3_to_qfr(x, S.D);
1673 : }
1674 : else
1675 : {
1676 70 : GEN d0 = gel(x,2);
1677 70 : x = gel(x,1);
1678 70 : if (!s) retmkvec2(qfr_1(x), real_0(precision(d0)));
1679 49 : if (s < 0) x = qfb_inv(x);
1680 49 : x = qfr5_init(x, d0, &S);
1681 49 : x = is_pm1(n)? qfr5_red(x, &S): qfr5_pow(x, n, &S);
1682 49 : x = qfr5_to_qfr(x, S.D, mulri(d0,n));
1683 : }
1684 77 : return gc_GEN(av, x);
1685 : }
1686 : GEN
1687 21 : qfbpowraw(GEN x, long n)
1688 : {
1689 21 : GEN q = check_qfbext("qfbpowraw",x);
1690 21 : return qfb_is_qfi(q)? qfipowraw(x,n): qfrpowraw(x,n);
1691 : }
1692 : /* same discriminant, no distance, no checks */
1693 : GEN
1694 10898367 : qfbpow_i(GEN x, GEN n) { return qfb_is_qfi(x)? qfipow(x,n): qfrpow(x,n); }
1695 : GEN
1696 1626395 : qfbpow(GEN x, GEN n)
1697 : {
1698 1626395 : GEN q = check_qfbext("qfbpow",x);
1699 1626395 : return qfb_is_qfi(q)? qfipow(x,n): qfrpow(x,n);
1700 : }
1701 : GEN
1702 1472324 : qfbpows(GEN x, long n)
1703 : {
1704 1472324 : long N[] = { evaltyp(t_INT) | _evallg(3), 0, 0};
1705 1472324 : affsi(n, N); return qfbpow(x, N);
1706 : }
1707 :
1708 : /* Prime forms attached to prime ideals of degree 1 */
1709 :
1710 : /* assume x != 0 a t_INT, p > 0
1711 : * Return a t_QFB, but discriminant sign is not checked: can be used for
1712 : * real forms as well */
1713 : GEN
1714 15083067 : primeform_u(GEN x, ulong p)
1715 : {
1716 15083067 : GEN c, y = cgetg(5, t_QFB);
1717 15083067 : pari_sp av = avma;
1718 : ulong b;
1719 : long s;
1720 :
1721 15083067 : s = mod8(x); if (signe(x) < 0 && s) s = 8-s;
1722 : /* 2 or 3 mod 4 */
1723 15083067 : if (s & 2) pari_err_DOMAIN("primeform", "disc % 4", ">",gen_1, x);
1724 15083060 : if (p == 2) {
1725 4296174 : switch(s) {
1726 642754 : case 0: b = 0; break;
1727 3301766 : case 1: b = 1; break;
1728 351654 : case 4: b = 2; break;
1729 0 : default: pari_err_SQRTN("primeform", mkintmod(x,utoi(p)) );
1730 0 : b = 0; /* -Wall */
1731 : }
1732 4296174 : c = shifti(subsi(s,x), -3);
1733 : } else {
1734 10786886 : b = Fl_sqrt(umodiu(x,p), p);
1735 10786886 : if (b == ~0UL) pari_err_SQRTN("primeform", mkintmod(x,utoi(p)) );
1736 : /* mod(b) != mod2(x) ? */
1737 10786886 : if ((b ^ s) & 1) b = p - b;
1738 10786886 : c = diviuexact(shifti(subii(sqru(b), x), -2), p);
1739 : }
1740 15083060 : gel(y,3) = gc_INT(av, c);
1741 15083060 : gel(y,4) = icopy(x);
1742 15083060 : gel(y,2) = utoi(b);
1743 15083060 : gel(y,1) = utoipos(p); return y;
1744 : }
1745 :
1746 : /* special case: p = 1 return unit form */
1747 : GEN
1748 135595 : primeform(GEN x, GEN p)
1749 : {
1750 135595 : const char *f = "primeform";
1751 : pari_sp av;
1752 135595 : long s, sx = signe(x), sp = signe(p);
1753 : GEN y, b, absp;
1754 :
1755 135595 : if (typ(x) != t_INT) pari_err_TYPE(f,x);
1756 135595 : if (typ(p) != t_INT) pari_err_TYPE(f,p);
1757 135595 : if (!sp) pari_err_DOMAIN(f,"p","=",gen_0,p);
1758 135595 : if (!sx) pari_err_DOMAIN(f,"D","=",gen_0,x);
1759 135595 : if (lgefint(p) == 3)
1760 : {
1761 135581 : ulong pp = p[2];
1762 135581 : if (pp == 1) {
1763 17918 : if (sx < 0) {
1764 : long r;
1765 11086 : if (sp < 0) pari_err_IMPL("negative definite t_QFB");
1766 11086 : r = mod4(x);
1767 11086 : if (r && r != 3) pari_err_DOMAIN(f,"disc % 4",">", gen_1,x);
1768 11086 : return qfi_1_by_disc(x);
1769 : }
1770 6832 : y = qfr_1_by_disc(x);
1771 6832 : if (sp < 0) { gel(y,1) = gen_m1; togglesign(gel(y,3)); }
1772 6832 : return y;
1773 : }
1774 117663 : y = primeform_u(x, pp);
1775 117656 : if (sx < 0) {
1776 89957 : if (sp < 0) pari_err_IMPL("negative definite t_QFB");
1777 89957 : return y;
1778 : }
1779 27699 : if (sp < 0) { togglesign(gel(y,1)); togglesign(gel(y,3)); }
1780 27699 : return gcopy( qfr3_to_qfr(y, x) );
1781 : }
1782 14 : s = mod8(x);
1783 14 : if (sx < 0)
1784 : {
1785 7 : if (sp < 0) pari_err_IMPL("negative definite t_QFB");
1786 7 : if (s) s = 8-s;
1787 : }
1788 14 : y = cgetg(5, t_QFB);
1789 : /* 2 or 3 mod 4 */
1790 14 : if (s & 2) pari_err_DOMAIN(f, "disc % 4", ">",gen_1, x);
1791 14 : absp = absi_shallow(p); av = avma;
1792 14 : b = Fp_sqrt(x, absp); if (!b) pari_err_SQRTN(f, mkintmod(x,absp));
1793 14 : s &= 1; /* s = x mod 2 */
1794 : /* mod(b) != mod2(x) ? [Warning: we may have b == 0] */
1795 14 : if ((!signe(b) && s) || mod2(b) != s) b = gc_INT(av, subii(absp,b));
1796 :
1797 14 : av = avma;
1798 14 : gel(y,3) = gc_INT(av, diviiexact(shifti(subii(sqri(b), x), -2), p));
1799 14 : gel(y,4) = icopy(x);
1800 14 : gel(y,2) = b;
1801 14 : gel(y,1) = icopy(p);
1802 14 : return y;
1803 : }
1804 :
1805 : static GEN
1806 2620772 : normforms(GEN D, GEN fa)
1807 : {
1808 : long i, j, k, lB, aN, sa;
1809 : GEN a, L, V, B, N, N2;
1810 2620772 : int D_odd = mpodd(D);
1811 2620772 : a = typ(fa) == t_INT ? fa: typ(fa) == t_VEC? gel(fa,1): factorback(fa);
1812 2620772 : sa = signe(a);
1813 2620772 : if (sa==0 || (signe(D)<0 && sa<0)) return NULL;
1814 1203972 : V = D_odd? Zn_quad_roots(fa, gen_1, shifti(subsi(1, D), -2))
1815 2551766 : : Zn_quad_roots(fa, gen_0, negi(shifti(D, -2)));
1816 2551766 : if (!V) return NULL;
1817 511966 : N = gel(V,1); B = gel(V,2); lB = lg(B);
1818 511966 : N2 = shifti(N,1);
1819 511966 : aN = itou(diviiexact(a, N)); /* |a|/N */
1820 511966 : L = cgetg((lB-1)*aN+1, t_VEC);
1821 2360568 : for (k = 1, i = 1; i < lB; i++)
1822 : {
1823 1848602 : GEN b = shifti(gel(B,i), 1), c, C;
1824 1848602 : if (D_odd) b = addiu(b, 1);
1825 1848602 : c = diviiexact(shifti(subii(sqri(b), D), -2), a);
1826 1848602 : for (j = 0;; b = addii(b, N2))
1827 : {
1828 2216676 : gel(L, k++) = mkqfb(a, b, c, D);
1829 2216676 : if (++j == aN) break;
1830 368074 : C = addii(b, N); if (aN > 1) C = diviuexact(C, aN);
1831 368074 : c = sa > 0? addii(c, C): subii(c, C);
1832 : }
1833 : }
1834 511966 : return L;
1835 : }
1836 :
1837 : /* Let M and N in SL2(Z), return (N*M^-1)[,1] */
1838 : static GEN
1839 344323 : SL2_div_mul_e1(GEN N, GEN M)
1840 : {
1841 344323 : GEN b = gcoeff(M,2,1), d = gcoeff(M,2,2);
1842 344323 : GEN A = mulii(gcoeff(N,1,1), d), B = mulii(gcoeff(N,1,2), b);
1843 344323 : GEN C = mulii(gcoeff(N,2,1), d), D = mulii(gcoeff(N,2,2), b);
1844 344323 : retmkvec2(subii(A,B), subii(C,D));
1845 : }
1846 : static GEN
1847 1445682 : qfisolve_normform(GEN Q, GEN P)
1848 : {
1849 1445682 : GEN a = gel(Q,1), N = gel(Q,2);
1850 1445682 : GEN M, b = qfi_redsl2_basecase(P, &M);
1851 1445682 : if (!qfb_equal(a,b)) return NULL;
1852 102130 : return SL2_div_mul_e1(N,M);
1853 : }
1854 :
1855 : /* Test equality modulo GL2 of two reduced forms */
1856 : static int
1857 61068 : GL2_qfb_equal(GEN a, GEN b)
1858 : {
1859 61068 : return equalii(gel(a,1),gel(b,1))
1860 11361 : && absequalii(gel(a,2),gel(b,2))
1861 72429 : && equalii(gel(a,3),gel(b,3));
1862 : }
1863 :
1864 : /* Q(u,v) = p; if s < 0 return that solution; else the set of all solutions */
1865 : static GEN
1866 48083 : allsols(GEN Q, long s, GEN u, GEN v)
1867 : {
1868 48083 : GEN w = mkvec2(u, v), b;
1869 48083 : if (signe(v) < 0) { u = negi(u); v = negi(v); } /* normalize for v >= 0 */
1870 48083 : w = mkvec2(u, v); if (s < 0) return w;
1871 41447 : if (!s) return mkvec(w);
1872 39018 : b = gel(Q,2); /* sum of the 2 solutions (if they exist) is -bv / a */
1873 39018 : if (signe(b))
1874 : { /* something to check */
1875 : GEN r, t;
1876 13433 : t = dvmdii(mulii(b, v), gel(Q,1), &r);
1877 13433 : if (signe(r)) return mkvec(w);
1878 1820 : u = addii(u, t);
1879 : }
1880 27405 : return mkvec2(w, mkvec2(negi(u), v));
1881 : }
1882 : static GEN
1883 223125 : qfisolvep_all(GEN Q, GEN p, long all)
1884 : {
1885 223125 : GEN R, U, V, M, N, x, q, D = qfb_disc(Q);
1886 223125 : long s = kronecker(D, p);
1887 :
1888 223125 : if (s < 0) return NULL;
1889 127050 : if (!all) s = -1; /* to indicate we want a single solution */
1890 : /* Solutions iff a class of maximal ideal above p is the class of Q;
1891 : * Two solutions iff (s > 0 and the class has order > 2), else one */
1892 127050 : if (!signe(gel(Q,2)))
1893 : { /* if principal form, use faster cornacchia */
1894 43729 : GEN a = gel(Q,1), c = gel(Q,3);
1895 43729 : if (equali1(a))
1896 : {
1897 38255 : if (!cornacchia(c, p, &M,&N)) return NULL;
1898 33768 : return allsols(Q, s, M, N);
1899 : }
1900 5474 : if (equali1(c))
1901 : {
1902 5194 : if (!cornacchia(a, p, &M,&N)) return NULL;
1903 721 : return allsols(Q, s, N, M);
1904 : }
1905 : }
1906 83601 : R = qfi_redsl2_basecase(Q, &U);
1907 83601 : if (equali1(gel(R,1)))
1908 : { /* principal form */
1909 22533 : if (!signe(gel(R,2)))
1910 : {
1911 4396 : if (!cornacchia(gel(R,3), p, &M,&N)) return NULL;
1912 812 : x = mkvec2(M,N);
1913 : }
1914 : else
1915 : { /* x^2 + xy + ((1-D)/4)y^2 = p <==> (2x + y)^2 - D y^2 = 4p */
1916 18137 : if (!cornacchia2(negi(D), p, &M, &N)) return NULL;
1917 2331 : x = subii(M,N); if (mpodd(x)) return NULL;
1918 2331 : x = mkvec2(shifti(x,-1), N);
1919 : }
1920 3143 : x = ZM_ZC_mul(U, x); x[0] = evaltyp(t_VEC) | _evallg(3); /* transpose */
1921 3143 : return allsols(Q, s, gel(x,1), gel(x,2));
1922 : }
1923 61068 : q = qfi_redsl2_basecase(primeform(D, p), &V);
1924 61068 : if (!GL2_qfb_equal(R,q)) return NULL;
1925 10451 : if (signe(gel(R,2)) != signe(gel(q,2))) gcoeff(V,2,1) = negi(gcoeff(V,2,1));
1926 10451 : x = SL2_div_mul_e1(U,V); return allsols(Q, s, gel(x,1), gel(x,2));
1927 : }
1928 : GEN
1929 0 : qfisolvep(GEN Q, GEN p)
1930 : {
1931 0 : pari_sp av = avma;
1932 0 : GEN x = qfisolvep_all(Q, p, 0);
1933 0 : return x? gc_GEN(av, x): gc_const(av, gen_0);
1934 : }
1935 :
1936 : static GEN
1937 770126 : qfrsolve_normform(GEN N, GEN Ps, GEN rd)
1938 : {
1939 770126 : pari_sp av = avma, btop;
1940 770126 : GEN M = N, P = qfr_redsl2_basecase(Ps, rd), Q = P;
1941 :
1942 770126 : btop = avma;
1943 : for(;;)
1944 : {
1945 5840681 : if (qfb_equal(gel(M,1), gel(P,1)))
1946 154084 : return gc_upto(av, SL2_div_mul_e1(gel(M,2),gel(P,2)));
1947 5686597 : if (qfb_equal(gel(N,1), gel(Q,1)))
1948 77658 : return gc_upto(av, SL2_div_mul_e1(gel(N,2),gel(Q,2)));
1949 5608939 : M = qfr_rhosl2(M, rd);
1950 5608939 : if (qfb_equal(gel(M,1), gel(N,1))) return gc_NULL(av);
1951 5201735 : Q = qfr_rhosl2(Q, rd);
1952 5201735 : if (qfb_equal(gel(P,1), gel(Q,1))) return gc_NULL(av);
1953 5070555 : if (gc_needed(btop, 1)) (void)gc_all(btop, 2, &M, &Q);
1954 : }
1955 : }
1956 :
1957 : GEN
1958 0 : qfrsolvep(GEN Q, GEN p)
1959 : {
1960 0 : pari_sp av = avma;
1961 0 : GEN N, x, rd, d = qfb_disc(Q);
1962 :
1963 0 : if (kronecker(d, p) < 0) return gc_const(av, gen_0);
1964 0 : rd = sqrti(d);
1965 0 : N = qfr_redsl2(Q, rd);
1966 0 : x = qfrsolve_normform(N, primeform(d, p), rd);
1967 0 : return x? gc_upto(av, x): gc_const(av, gen_0);
1968 : }
1969 :
1970 : static GEN
1971 1863022 : known_prime(GEN v)
1972 : {
1973 1863022 : GEN p, e, fa = check_arith_all(v, "qfbsolve");
1974 1863022 : if (!fa) return BPSW_psp(v)? v: NULL;
1975 42154 : if (lg(gel(fa,1)) != 2) return NULL;
1976 29428 : p = gcoeff(fa,1,1);
1977 29428 : e = gcoeff(fa,1,2);
1978 29428 : return (equali1(e) && !is_pm1(p) && signe(p) > 0)? p: NULL;
1979 : }
1980 : static GEN
1981 2215808 : qfsolve_normform(GEN Q, GEN f, GEN rd)
1982 2215808 : { return rd? qfrsolve_normform(Q, f, rd): qfisolve_normform(Q, f); }
1983 : static GEN
1984 2843897 : qfbsolve_primitive_i(GEN Q, GEN rd, GEN *Qr, GEN fa, long all)
1985 : {
1986 : GEN x, W, F, p;
1987 : long i, j, l;
1988 2843897 : if (!rd && (p = known_prime(fa))) return qfisolvep_all(Q, p, all);
1989 2620772 : F = normforms(qfb_disc(Q), fa);
1990 2620772 : if (!F) return NULL;
1991 511966 : if (!*Qr) *Qr = qfbredsl2(Q, rd);
1992 511966 : l = lg(F); W = all? cgetg(l, t_VEC): NULL;
1993 2727263 : for (j = i = 1; i < l; i++)
1994 2215808 : if ((x = qfsolve_normform(*Qr, gel(F,i), rd)))
1995 : {
1996 333872 : if (!all) return x;
1997 333361 : gel(W,j++) = x;
1998 : }
1999 511455 : if (j == 1) return NULL;
2000 127456 : setlg(W,j); return lexsort(W);
2001 : }
2002 :
2003 : static GEN
2004 2838598 : qfb_initrd(GEN Q) { GEN d = qfb_disc(Q); return signe(d) > 0? sqrti(d): NULL; }
2005 : static GEN
2006 2828371 : qfbsolve_primitive(GEN Q, GEN fa, long all)
2007 : {
2008 2828371 : GEN x, Qr = NULL, rdQ = qfb_initrd(Q);
2009 2828371 : x = qfbsolve_primitive_i(Q, rdQ, &Qr, fa, all);
2010 2828371 : if (!x) return cgetg(1, t_VEC);
2011 174832 : return x;
2012 : }
2013 :
2014 : /* f / g^2 */
2015 : static GEN
2016 5299 : famat_divsqr(GEN f, GEN g)
2017 5299 : { return famat_reduce(famat_div_shallow(f, famat_pows_shallow(g,2))); }
2018 : static GEN
2019 10227 : qfbsolve_all(GEN Q, GEN n, long all)
2020 : {
2021 10227 : GEN W, Qr = NULL, fa = factorint(n, 0), rdQ = qfb_initrd(Q);
2022 10227 : GEN D = divisors_factored(mkmat2(gel(fa,1), gshift(gel(fa,2),-1)));
2023 10227 : long i, j, l = lg(D);
2024 10227 : W = all? cgetg(l, t_VEC): NULL;
2025 25151 : for (i = j = 1; i < l; i++)
2026 : {
2027 15526 : GEN w, d = gel(D,i), FA = i == 1? fa: famat_divsqr(fa, gel(d,2));
2028 15526 : if ((w = qfbsolve_primitive_i(Q, rdQ, &Qr, FA, all)))
2029 : {
2030 1218 : if (i != 1) w = RgV_Rg_mul(w, gel(d,1));
2031 1218 : if (!all) return w;
2032 616 : gel(W,j++) = w;
2033 : }
2034 : }
2035 9625 : if (j == 1) return cgetg(1, t_VEC);
2036 525 : setlg(W,j); return lexsort(shallowconcat1(W));
2037 : }
2038 :
2039 : GEN
2040 2838605 : qfbsolve(GEN Q, GEN n, long fl)
2041 : {
2042 2838605 : pari_sp av = avma;
2043 2838605 : if (typ(Q) != t_QFB) pari_err_TYPE("qfbsolve",Q);
2044 2838605 : if (fl < 0 || fl > 3) pari_err_FLAG("qfbsolve");
2045 5666969 : return gc_GEN(av, (fl & 2)? qfbsolve_all(Q, n, fl & 1)
2046 2828371 : : qfbsolve_primitive(Q, n, fl & 1));
2047 : }
2048 :
2049 : /* 1 if there exists x,y such that x^2 + dy^2 = p, 0 otherwise;
2050 : * Assume d > 0 and p is prime */
2051 : long
2052 55328 : cornacchia(GEN d, GEN p, GEN *px, GEN *py)
2053 : {
2054 55328 : pari_sp av = avma;
2055 : GEN b, c, r;
2056 :
2057 55328 : *px = *py = gen_0;
2058 55328 : b = subii(p, d);
2059 55328 : if (signe(b) < 0) return gc_long(av,0);
2060 55118 : if (signe(b) == 0) { *py = gen_1; return gc_long(av,1); }
2061 55111 : b = Fp_sqrt(b, p); /* sqrt(-d) */
2062 55111 : if (!b) return gc_long(av,0);
2063 51380 : b = gmael(halfgcdii(p, b), 2, 2);
2064 51380 : c = dvmdii(subii(p, sqri(b)), d, &r);
2065 51380 : if (r != gen_0 || !Z_issquareall(c, &c)) return gc_long(av,0);
2066 35532 : set_avma(av);
2067 35532 : *px = icopy(b);
2068 35532 : *py = icopy(c); return 1;
2069 : }
2070 :
2071 : static GEN
2072 2595837 : lastqi(GEN Q)
2073 : {
2074 2595837 : GEN s = gcoeff(Q,1,1), q = gcoeff(Q,1,2), p = absi_shallow(gcoeff(Q,2,2));
2075 2595837 : if (!signe(q)) return gen_0;
2076 2595648 : if (!signe(s)) return p;
2077 2588048 : if (is_pm1(q)) return subiu(p,1);
2078 2588048 : return divii(p, absi_shallow(q));
2079 : }
2080 :
2081 : static long
2082 2595851 : cornacchia2_i(long av, GEN d, GEN p, GEN b, GEN px4, GEN *px, GEN *py)
2083 : {
2084 : GEN M, Q, V, c, r, b2;
2085 2595851 : if (!signe(b)) { /* d = p,2p,3p,4p */
2086 14 : set_avma(av);
2087 14 : if (absequalii(d, px4)){ *py = gen_1; return 1; }
2088 14 : if (absequalii(d, p)) { *py = gen_2; return 1; }
2089 0 : return 0;
2090 : }
2091 2595837 : if (mod2(b) != mod2(d)) b = subii(p,b);
2092 2595837 : M = halfgcdii(shifti(p,1), b); Q = gel(M,1); V = gel(M, 2);
2093 2595837 : b = addii(mulii(gel(V,1), lastqi(Q)), gel(V,2));
2094 2595837 : b2 = sqri(b);
2095 2595837 : if (cmpii(b2,px4) > 0)
2096 : {
2097 2585851 : b = gel(V,1); b2 = sqri(b);
2098 2585851 : if (cmpii(b2,px4) > 0) { b = gel(V,2); b2 = sqri(b); }
2099 : }
2100 2595837 : c = dvmdii(subii(px4, b2), d, &r);
2101 2595837 : if (r != gen_0 || !Z_issquareall(c, &c)) return gc_long(av,0);
2102 2555027 : set_avma(av);
2103 2555027 : *px = icopy(b);
2104 2555027 : *py = icopy(c); return 1;
2105 : }
2106 :
2107 : /* 1 if there exists x,y such that x^2 + dy^2 = 4p, 0 otherwise;
2108 : * Assume d > 0 is congruent to 0 or 3 mod 4 and p is prime */
2109 : long
2110 2561278 : cornacchia2(GEN d, GEN p, GEN *px, GEN *py)
2111 : {
2112 2561278 : pari_sp av = avma;
2113 2561278 : GEN b, p4 = shifti(p,2);
2114 :
2115 2561278 : *px = *py = gen_0;
2116 2561278 : if (abscmpii(p4, d) < 0) return gc_long(av,0);
2117 2560459 : if (absequaliu(p, 2))
2118 : {
2119 7 : set_avma(av);
2120 7 : switch (itou_or_0(d)) {
2121 0 : case 4: *px = gen_2; break;
2122 0 : case 7: *px = gen_1; break;
2123 7 : default: return 0;
2124 : }
2125 0 : *py = gen_1; return 1;
2126 : }
2127 2560452 : b = Fp_sqrt(negi(d), p);
2128 2560452 : if (!b) return gc_long(av,0);
2129 2560368 : return cornacchia2_i(av, d, p, b, p4, px, py);
2130 : }
2131 :
2132 : /* 1 if there exists x,y such that x^2 + dy^2 = 4p [p prime], 0 otherwise */
2133 : long
2134 35483 : cornacchia2_sqrt(GEN d, GEN p, GEN b, GEN *px, GEN *py)
2135 : {
2136 35483 : pari_sp av = avma;
2137 35483 : GEN p4 = shifti(p,2);
2138 35483 : *px = *py = gen_0;
2139 35483 : if (abscmpii(p4, d) < 0) return gc_long(av,0);
2140 35483 : return cornacchia2_i(av, d, p, b, p4, px, py);
2141 : }
2142 :
2143 : GEN
2144 7630 : qfbcornacchia(GEN d, GEN p)
2145 : {
2146 7630 : pari_sp av = avma;
2147 : GEN x, y;
2148 7630 : if (typ(d) != t_INT || signe(d) <= 0) pari_err_TYPE("qfbcornacchia", d);
2149 7630 : if (typ(p) != t_INT || cmpiu(p, 2) < 0) pari_err_TYPE("qfbcornacchia", p);
2150 7630 : if (mod4(p)? cornacchia(d, p, &x, &y): cornacchia2(d, shifti(p, -2), &x, &y))
2151 287 : return gc_GEN(av, mkvec2(x, y));
2152 7343 : retgc_const(av, cgetg(1, t_VEC));
2153 : }
|
