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
check_quaddisc(GEN x,long * s,long * r,const char * f)23 check_quaddisc(GEN x, long *s, long *r, const char *f)
24 {
25 if (typ(x) != t_INT) pari_err_TYPE(f,x);
26 *s = signe(x);
27 if (Z_issquare(x)) pari_err_DOMAIN(f,"issquare(disc)","=", gen_1,x);
28 *r = mod4(x); if (*s < 0 && *r) *r = 4 - *r;
29 if (*r > 1) pari_err_DOMAIN(f,"disc % 4",">", gen_1,x);
30 }
31 void
check_quaddisc_real(GEN x,long * r,const char * f)32 check_quaddisc_real(GEN x, long *r, const char *f)
33 {
34 long sx; check_quaddisc(x, &sx, r, f);
35 if (sx < 0) pari_err_DOMAIN(f, "disc","<",gen_0,x);
36 }
37 void
check_quaddisc_imag(GEN x,long * r,const char * f)38 check_quaddisc_imag(GEN x, long *r, const char *f)
39 {
40 long sx; check_quaddisc(x, &sx, r, f);
41 if (sx > 0) pari_err_DOMAIN(f, "disc",">",gen_0,x);
42 }
43
44 /* X^2 + b X + c is the canonical quadratic t_POL of discriminant D.
45 * Dodd is nonzero iff D is odd */
46 static void
quadpoly_bc(GEN D,long Dodd,GEN * b,GEN * c)47 quadpoly_bc(GEN D, long Dodd, GEN *b, GEN *c)
48 {
49 if (Dodd)
50 {
51 pari_sp av = avma;
52 *b = gen_m1;
53 *c = gerepileuptoint(av, shifti(subui(1,D), -2));
54 }
55 else
56 {
57 *b = gen_0;
58 *c = shifti(D,-2); togglesign(*c);
59 }
60 }
61 /* X^2 - X - (D-1)/4 or X^2 - D/4 */
62 GEN
quadpoly(GEN D)63 quadpoly(GEN D)
64 {
65 long Dmod4, s;
66 GEN b, c, y = cgetg(5,t_POL);
67 check_quaddisc(D, &s, &Dmod4, "quadpoly");
68 y[1] = evalsigne(1) | evalvarn(0);
69 quadpoly_bc(D, Dmod4, &b,&c);
70 gel(y,2) = c;
71 gel(y,3) = b;
72 gel(y,4) = gen_1; return y;
73 }
74
75 GEN
quadpoly0(GEN x,long v)76 quadpoly0(GEN x, long v)
77 {
78 GEN T = quadpoly(x);
79 if (v > 0) setvarn(T, v);
80 return T;
81 }
82
83 GEN
quadgen(GEN x)84 quadgen(GEN x)
85 { retmkquad(quadpoly(x), gen_0, gen_1); }
86
87 GEN
quadgen0(GEN x,long v)88 quadgen0(GEN x, long v)
89 {
90 if (v==-1) v = fetch_user_var("w");
91 retmkquad(quadpoly0(x, v), gen_0, gen_1);
92 }
93
94 /***********************************************************************/
95 /** **/
96 /** BINARY QUADRATIC FORMS **/
97 /** **/
98 /***********************************************************************/
99 GEN
qfi(GEN x,GEN y,GEN z)100 qfi(GEN x, GEN y, GEN z)
101 {
102 if (signe(x) < 0) pari_err_IMPL("negative definite t_QFI");
103 retmkqfi(icopy(x),icopy(y),icopy(z));
104 }
105 GEN
qfr(GEN x,GEN y,GEN z,GEN d)106 qfr(GEN x, GEN y, GEN z, GEN d)
107 {
108 if (typ(d) != t_REAL) pari_err_TYPE("qfr",d);
109 retmkqfr(icopy(x),icopy(y),icopy(z),rcopy(d));
110 }
111
112 GEN
Qfb0(GEN x,GEN y,GEN z,GEN d,long prec)113 Qfb0(GEN x, GEN y, GEN z, GEN d, long prec)
114 {
115 pari_sp av = avma;
116 GEN D;
117 long s, r;
118 if (typ(x)!=t_INT) pari_err_TYPE("Qfb",x);
119 if (typ(y)!=t_INT) pari_err_TYPE("Qfb",y);
120 if (typ(z)!=t_INT) pari_err_TYPE("Qfb",z);
121 D = qfb_disc3(x,y,z);
122 check_quaddisc(D, &s, &r, "Qfb");
123 set_avma(av);
124 if (s < 0) return qfi(x, y, z);
125
126 d = d? gtofp(d,prec): real_0(prec);
127 return qfr(x,y,z,d);
128 }
129
130 /* composition */
131 static void
qfb_sqr(GEN z,GEN x)132 qfb_sqr(GEN z, GEN x)
133 {
134 GEN c, d1, x2, v1, v2, c3, m, p1, r;
135
136 d1 = bezout(gel(x,2),gel(x,1),&x2, NULL); /* usually 1 */
137 c = gel(x,3);
138 m = mulii(c,x2);
139 if (equali1(d1))
140 v1 = v2 = gel(x,1);
141 else
142 {
143 v1 = diviiexact(gel(x,1),d1);
144 v2 = mulii(v1, gcdii(d1,c)); /* = v1 iff x primitive */
145 c = mulii(c, d1);
146 }
147 togglesign(m);
148 r = modii(m,v2);
149 p1 = mulii(r, v1);
150 c3 = addii(c, mulii(r,addii(gel(x,2),p1)));
151 gel(z,1) = mulii(v1,v2);
152 gel(z,2) = addii(gel(x,2), shifti(p1,1));
153 gel(z,3) = diviiexact(c3,v2);
154 }
155 /* z <- x * y */
156 static void
qfb_comp(GEN z,GEN x,GEN y)157 qfb_comp(GEN z, GEN x, GEN y)
158 {
159 GEN n, c, d, y1, v1, v2, c3, m, p1, r;
160
161 if (x == y) { qfb_sqr(z,x); return; }
162 n = shifti(subii(gel(y,2),gel(x,2)), -1);
163 v1 = gel(x,1);
164 v2 = gel(y,1);
165 c = gel(y,3);
166 d = bezout(v2,v1,&y1,NULL);
167 if (equali1(d))
168 m = mulii(y1,n);
169 else
170 {
171 GEN s = subii(gel(y,2), n);
172 GEN x2, y2, d1 = bezout(s,d,&x2,&y2); /* x2 s + y2 (x1 v1 + y1 v2) = d1 */
173 if (!equali1(d1))
174 {
175 v1 = diviiexact(v1,d1);
176 v2 = diviiexact(v2,d1); /* gcd = 1 iff x or y primitive */
177 v1 = mulii(v1, gcdii(c,gcdii(gel(x,3),gcdii(d1,n))));
178 c = mulii(c, d1);
179 }
180 m = addii(mulii(mulii(y1,y2),n), mulii(gel(y,3),x2));
181 }
182 togglesign(m);
183 r = modii(m, v1);
184 p1 = mulii(r, v2);
185 c3 = addii(c, mulii(r,addii(gel(y,2),p1)));
186 gel(z,1) = mulii(v1,v2);
187 gel(z,2) = addii(gel(y,2), shifti(p1,1));
188 gel(z,3) = diviiexact(c3,v1);
189 }
190
191 static GEN redimag_av(pari_sp av, GEN q);
192 static GEN
qficomp0(GEN x,GEN y,int raw)193 qficomp0(GEN x, GEN y, int raw)
194 {
195 pari_sp av = avma;
196 GEN z = cgetg(4,t_QFI);
197 qfb_comp(z, x,y);
198 if (raw) return gerepilecopy(av,z);
199 return redimag_av(av, z);
200 }
201 static GEN
qfrcomp0(GEN x,GEN y,int raw)202 qfrcomp0(GEN x, GEN y, int raw)
203 {
204 pari_sp av = avma;
205 GEN z = cgetg(5,t_QFR);
206 qfb_comp(z, x,y); gel(z,4) = addrr(gel(x,4),gel(y,4));
207 if (raw) return gerepilecopy(av,z);
208 return gerepileupto(av, redreal(z));
209 }
210 GEN
qfrcomp(GEN x,GEN y)211 qfrcomp(GEN x, GEN y) { return qfrcomp0(x,y,0); }
212 GEN
qfrcompraw(GEN x,GEN y)213 qfrcompraw(GEN x, GEN y) { return qfrcomp0(x,y,1); }
214 GEN
qficomp(GEN x,GEN y)215 qficomp(GEN x, GEN y) { return qficomp0(x,y,0); }
216 GEN
qficompraw(GEN x,GEN y)217 qficompraw(GEN x, GEN y) { return qficomp0(x,y,1); }
218 GEN
qfbcompraw(GEN x,GEN y)219 qfbcompraw(GEN x, GEN y)
220 {
221 long tx = typ(x);
222 if (typ(y) != tx) pari_err_TYPE2("*",x,y);
223 switch(tx) {
224 case t_QFI: return qficompraw(x,y);
225 case t_QFR: return qfrcompraw(x,y);
226 }
227 pari_err_TYPE("composition",x);
228 return NULL; /* LCOV_EXCL_LINE */
229 }
230
231 static GEN
qfisqr0(GEN x,long raw)232 qfisqr0(GEN x, long raw)
233 {
234 pari_sp av = avma;
235 GEN z = cgetg(4,t_QFI);
236
237 if (typ(x)!=t_QFI) pari_err_TYPE("composition",x);
238 qfb_sqr(z,x);
239 if (raw) return gerepilecopy(av,z);
240 return redimag_av(av, z);
241 }
242 static GEN
qfrsqr0(GEN x,long raw)243 qfrsqr0(GEN x, long raw)
244 {
245 pari_sp av = avma;
246 GEN z = cgetg(5,t_QFR);
247
248 if (typ(x)!=t_QFR) pari_err_TYPE("composition",x);
249 qfb_sqr(z,x); gel(z,4) = shiftr(gel(x,4),1);
250 if (raw) return gerepilecopy(av,z);
251 return gerepileupto(av, redreal(z));
252 }
253 GEN
qfrsqr(GEN x)254 qfrsqr(GEN x) { return qfrsqr0(x,0); }
255 GEN
qfrsqrraw(GEN x)256 qfrsqrraw(GEN x) { return qfrsqr0(x,1); }
257 GEN
qfisqr(GEN x)258 qfisqr(GEN x) { return qfisqr0(x,0); }
259 GEN
qfisqrraw(GEN x)260 qfisqrraw(GEN x) { return qfisqr0(x,1); }
261
262 static GEN
qfr_1_by_disc(GEN D,long prec)263 qfr_1_by_disc(GEN D, long prec)
264 {
265 GEN y = cgetg(5,t_QFR), isqrtD;
266 pari_sp av = avma;
267 long r;
268
269 check_quaddisc_real(D, &r, "qfr_1_by_disc");
270 gel(y,1) = gen_1; isqrtD = sqrti(D);
271 if ((r & 1) != mod2(isqrtD)) /* we know isqrtD > 0 */
272 isqrtD = gerepileuptoint(av, subiu(isqrtD,1));
273 gel(y,2) = isqrtD; av = avma;
274 gel(y,3) = gerepileuptoint(av, shifti(subii(sqri(isqrtD), D),-2));
275 gel(y,4) = real_0(prec); return y;
276 }
277 GEN
qfr_1(GEN x)278 qfr_1(GEN x)
279 {
280 if (typ(x) != t_QFR) pari_err_TYPE("qfr_1",x);
281 return qfr_1_by_disc(qfb_disc(x), precision(gel(x,4)));
282 }
283
284 static void
qfr_1_fill(GEN y,struct qfr_data * S)285 qfr_1_fill(GEN y, struct qfr_data *S)
286 {
287 pari_sp av = avma;
288 GEN y2 = S->isqrtD;
289 gel(y,1) = gen_1;
290 if (mod2(S->D) != mod2(y2)) y2 = subiu(y,1);
291 gel(y,2) = y2; av = avma;
292 gel(y,3) = gerepileuptoint(av, shifti(subii(sqri(y2), S->D),-2));
293 }
294 static GEN
qfr5_1(struct qfr_data * S,long prec)295 qfr5_1(struct qfr_data *S, long prec)
296 {
297 GEN y = cgetg(6, t_VEC);
298 qfr_1_fill(y, S);
299 gel(y,4) = gen_0;
300 gel(y,5) = real_1(prec); return y;
301 }
302 static GEN
qfr3_1(struct qfr_data * S)303 qfr3_1(struct qfr_data *S)
304 {
305 GEN y = cgetg(4, t_VEC);
306 qfr_1_fill(y, S); return y;
307 }
308
309 /* Assume D < 0 is the discriminant of a t_QFI */
310 static GEN
qfi_1_by_disc(GEN D)311 qfi_1_by_disc(GEN D)
312 {
313 GEN b,c, y = cgetg(4,t_QFI);
314 quadpoly_bc(D, mod2(D), &b,&c);
315 if (b == gen_m1) b = gen_1;
316 gel(y,1) = gen_1;
317 gel(y,2) = b;
318 gel(y,3) = c; return y;
319 }
320 GEN
qfi_1(GEN x)321 qfi_1(GEN x)
322 {
323 if (typ(x) != t_QFI) pari_err_TYPE("qfi_1",x);
324 return qfi_1_by_disc(qfb_disc(x));
325 }
326
327 static GEN
invraw(GEN x)328 invraw(GEN x)
329 {
330 GEN y = gcopy(x);
331 if (typ(y) == t_QFR) togglesign(gel(y,4));
332 togglesign(gel(y,2)); return y;
333 }
334 GEN
qfrpowraw(GEN x,long n)335 qfrpowraw(GEN x, long n)
336 {
337 pari_sp av = avma;
338 long m;
339 GEN y;
340
341 if (typ(x) != t_QFR) pari_err_TYPE("qfrpowraw",x);
342 if (!n) return qfr_1(x);
343 if (n== 1) return gcopy(x);
344 if (n==-1) return invraw(x);
345
346 y = NULL; m = labs(n);
347 for (; m>1; m>>=1)
348 {
349 if (m&1) y = y? qfrcompraw(y,x): x;
350 x = qfrsqrraw(x);
351 }
352 y = y? qfrcompraw(y,x): x;
353 if (n < 0) y = invraw(y);
354 return gerepileupto(av,y);
355 }
356 GEN
qfipowraw(GEN x,long n)357 qfipowraw(GEN x, long n)
358 {
359 pari_sp av = avma;
360 long m;
361 GEN y;
362
363 if (typ(x) != t_QFI) pari_err_TYPE("qfipow",x);
364 if (!n) return qfi_1(x);
365 if (n== 1) return gcopy(x);
366 if (n==-1) return invraw(x);
367
368 y = NULL; m = labs(n);
369 for (; m>1; m>>=1)
370 {
371 if (m&1) y = y? qficompraw(y,x): x;
372 x = qfisqrraw(x);
373 }
374 y = y? qficompraw(y,x): x;
375 if (n < 0) y = invraw(y);
376 return gerepileupto(av,y);
377 }
378
379 GEN
qfbpowraw(GEN x,long n)380 qfbpowraw(GEN x, long n)
381 { return (typ(x)==t_QFI)? qfipowraw(x,n): qfrpowraw(x,n); }
382
383 static long
parteucl(GEN L,GEN * d,GEN * v3,GEN * v,GEN * v2)384 parteucl(GEN L, GEN *d, GEN *v3, GEN *v, GEN *v2)
385 {
386 long z;
387 *v = gen_0; *v2 = gen_1;
388 for (z=0; abscmpii(*v3,L) > 0; z++)
389 {
390 GEN t3, t2 = subii(*v, mulii(truedvmdii(*d,*v3,&t3),*v2));
391 *v = *v2; *d = *v3; *v2 = t2; *v3 = t3;
392 }
393 return z;
394 }
395
396 /* composition: Shanks' NUCOMP & NUDUPL */
397 /* L = floor((|d|/4)^(1/4)) */
398 GEN
nucomp(GEN x,GEN y,GEN L)399 nucomp(GEN x, GEN y, GEN L)
400 {
401 pari_sp av = avma;
402 long z;
403 GEN a, a1, a2, b2, b, d, d1, g, n, p1, q1, q2, s, u, u1, v, v2, v3, Q;
404
405 if (x==y) return nudupl(x,L);
406 if (typ(x) != t_QFI) pari_err_TYPE("nucomp",x);
407 if (typ(y) != t_QFI) pari_err_TYPE("nucomp",y);
408
409 if (abscmpii(gel(x,1),gel(y,1)) < 0) swap(x, y);
410 s = shifti(addii(gel(x,2),gel(y,2)), -1);
411 n = subii(gel(y,2), s);
412 a1 = gel(x,1);
413 a2 = gel(y,1); d = bezout(a2,a1,&u,&v);
414 if (equali1(d)) { a = negi(mulii(u,n)); d1 = d; }
415 else if (dvdii(s,d)) /* d | s */
416 {
417 a = negi(mulii(u,n)); d1 = d;
418 a1 = diviiexact(a1, d1);
419 a2 = diviiexact(a2, d1);
420 s = diviiexact(s, d1);
421 }
422 else
423 {
424 GEN p2, l;
425 d1 = bezout(s,d,&u1,NULL);
426 if (!equali1(d1))
427 {
428 a1 = diviiexact(a1,d1);
429 a2 = diviiexact(a2,d1);
430 s = diviiexact(s,d1);
431 d = diviiexact(d,d1);
432 }
433 p1 = remii(gel(x,3),d);
434 p2 = remii(gel(y,3),d);
435 l = modii(mulii(negi(u1), addii(mulii(u,p1),mulii(v,p2))), d);
436 a = subii(mulii(l,diviiexact(a1,d)), mulii(u,diviiexact(n,d)));
437 }
438 a = modii(a,a1); p1 = subii(a,a1); if (abscmpii(a,p1) > 0) a = p1;
439 d = a1; v3 = a; z = parteucl(L, &d,&v3, &v,&v2);
440 Q = cgetg(4,t_QFI);
441 if (!z) {
442 g = diviiexact(addii(mulii(v3,s),gel(y,3)), d);
443 b = a2;
444 b2 = gel(y,2);
445 v2 = d1;
446 gel(Q,1) = mulii(d,b);
447 } else {
448 GEN e, q3, q4;
449 if (z&1) { v3 = negi(v3); v2 = negi(v2); }
450 b = diviiexact(addii(mulii(a2,d), mulii(n,v)), a1);
451 e = diviiexact(addii(mulii(s,d),mulii(gel(y,3),v)), a1);
452 q3 = mulii(e,v2);
453 q4 = subii(q3,s);
454 b2 = addii(q3,q4);
455 g = diviiexact(q4,v);
456 if (!equali1(d1)) { v2 = mulii(d1,v2); v = mulii(d1,v); b2 = mulii(d1,b2); }
457 gel(Q,1) = addii(mulii(d,b), mulii(e,v));
458 }
459 q1 = mulii(b, v3);
460 q2 = addii(q1,n);
461 gel(Q,2) = addii(b2, z? addii(q1,q2): shifti(q1, 1));
462 gel(Q,3) = addii(mulii(v3,diviiexact(q2,d)), mulii(g,v2));
463 return redimag_av(av, Q);
464 }
465
466 GEN
nudupl(GEN x,GEN L)467 nudupl(GEN x, GEN L)
468 {
469 pari_sp av = avma;
470 long z;
471 GEN u, v, d, d1, p1, a, b, c, a2, b2, c2, Q, v2, v3, g;
472
473 if (typ(x) != t_QFI) pari_err_TYPE("nudupl",x);
474 a = gel(x,1);
475 b = gel(x,2);
476 d1 = bezout(b,a, &u,NULL);
477 if (!equali1(d1))
478 {
479 a = diviiexact(a, d1);
480 b = diviiexact(b, d1);
481 }
482 c = modii(negi(mulii(u,gel(x,3))), a);
483 p1 = subii(c,a); if (abscmpii(c,p1) > 0) c = p1;
484 d = a; v3 = c; z = parteucl(L, &d,&v3, &v,&v2);
485 a2 = sqri(d);
486 c2 = sqri(v3);
487 Q = cgetg(4,t_QFI);
488 if (!z) {
489 g = diviiexact(addii(mulii(v3,b),gel(x,3)), d);
490 b2 = gel(x,2);
491 v2 = d1;
492 gel(Q,1) = a2;
493 } else {
494 GEN e;
495 if (z&1) { v = negi(v); d = negi(d); }
496 e = diviiexact(addii(mulii(gel(x,3),v), mulii(b,d)), a);
497 g = diviiexact(subii(mulii(e,v2), b), v);
498 b2 = addii(mulii(e,v2), mulii(v,g));
499 if (!equali1(d1)) { b2 = mulii(d1,b2); v = mulii(d1,v); v2 = mulii(d1,v2); }
500 gel(Q,1) = addii(a2, mulii(e,v));
501 }
502 gel(Q,2) = addii(b2, subii(sqri(addii(d,v3)), addii(a2,c2)));
503 gel(Q,3) = addii(c2, mulii(g,v2));
504 return redimag_av(av, Q);
505 }
506
507 static GEN
mul_nucomp(void * l,GEN x,GEN y)508 mul_nucomp(void *l, GEN x, GEN y) { return nucomp(x, y, (GEN)l); }
509 static GEN
mul_nudupl(void * l,GEN x)510 mul_nudupl(void *l, GEN x) { return nudupl(x, (GEN)l); }
511 GEN
nupow(GEN x,GEN n,GEN L)512 nupow(GEN x, GEN n, GEN L)
513 {
514 pari_sp av;
515 GEN y, D;
516
517 if (typ(n) != t_INT) pari_err_TYPE("nupow",n);
518 if (typ(x) != t_QFI) pari_err_TYPE("nupow",x);
519 if (gequal1(n)) return gcopy(x);
520 av = avma;
521 D = qfb_disc(x);
522 y = qfi_1_by_disc(D);
523 if (!signe(n)) return y;
524 if (!L) L = sqrtnint(absi_shallow(D), 4);
525 y = gen_pow_i(x, n, (void*)L, &mul_nudupl, &mul_nucomp);
526 if (signe(n) < 0
527 && !absequalii(gel(y,1),gel(y,2))
528 && !absequalii(gel(y,1),gel(y,3))) togglesign(gel(y,2));
529 return gerepilecopy(av, y);
530 }
531
532 /* Reduction */
533
534 /* assume a > 0. Write b = q*2a + r, with -a < r <= a */
535 static GEN
dvmdii_round(GEN b,GEN a,GEN * r)536 dvmdii_round(GEN b, GEN a, GEN *r)
537 {
538 GEN a2 = shifti(a, 1), q = dvmdii(b, a2, r);
539 if (signe(b) >= 0) {
540 if (abscmpii(*r, a) > 0) { q = addiu(q, 1); *r = subii(*r, a2); }
541 } else { /* r <= 0 */
542 if (abscmpii(*r, a) >= 0){ q = subiu(q, 1); *r = addii(*r, a2); }
543 }
544 return q;
545 }
546 /* Assume 0 < a <= LONG_MAX. Ensure no overflow */
547 static long
dvmdsu_round(long b,ulong a,long * r)548 dvmdsu_round(long b, ulong a, long *r)
549 {
550 ulong a2 = a << 1, q, ub, ur;
551 if (b >= 0) {
552 ub = b;
553 q = ub / a2;
554 ur = ub % a2;
555 if (ur > a) { ur -= a; q++; *r = (long)ur; *r -= (long)a; }
556 else *r = (long)ur;
557 return (long)q;
558 } else { /* r <= 0 */
559 ub = (ulong)-b; /* |b| */
560 q = ub / a2;
561 ur = ub % a2;
562 if (ur >= a) { ur -= a; q++; *r = (long)ur; *r = (long)a - *r; }
563 else *r = -(long)ur;
564 return -(long)q;
565 }
566 }
567 /* reduce b mod 2*a. Update b,c */
568 static void
REDB(GEN a,GEN * b,GEN * c)569 REDB(GEN a, GEN *b, GEN *c)
570 {
571 GEN r, q = dvmdii_round(*b, a, &r);
572 if (!signe(q)) return;
573 *c = subii(*c, mulii(q, shifti(addii(*b, r),-1)));
574 *b = r;
575 }
576 /* Assume a > 0. Reduce b mod 2*a. Update b,c */
577 static void
sREDB(ulong a,long * b,ulong * c)578 sREDB(ulong a, long *b, ulong *c)
579 {
580 long r, q;
581 ulong uz;
582 if (a > LONG_MAX) return; /* b already reduced */
583 q = dvmdsu_round(*b, a, &r);
584 if (q == 0) return;
585 /* Final (a,r,c2) satisfies |r| <= |b| hence c2 <= c, c2 = c - q*z,
586 * where z = (b+r) / 2, representable as long, has the same sign as q. */
587 if (*b < 0)
588 { /* uz = -z >= 0, q < 0 */
589 if (r >= 0) /* different signs=>no overflow, exact division */
590 uz = (ulong)-((*b + r)>>1);
591 else
592 {
593 ulong ub = (ulong)-*b, ur = (ulong)-r;
594 uz = (ub + ur) >> 1;
595 }
596 *c -= (-q) * uz; /* c -= qz */
597 }
598 else
599 { /* uz = z >= 0, q > 0 */
600 if (r <= 0)
601 uz = (*b + r)>>1;
602 else
603 {
604 ulong ub = (ulong)*b, ur = (ulong)r;
605 uz = ((ub + ur) >> 1);
606 }
607 *c -= q * uz; /* c -= qz */
608 }
609 *b = r;
610 }
611 static void
REDBU(GEN a,GEN * b,GEN * c,GEN u1,GEN * u2)612 REDBU(GEN a, GEN *b, GEN *c, GEN u1, GEN *u2)
613 { /* REDB(a,b,c) */
614 GEN r, q = dvmdii_round(*b, a, &r);
615 *c = subii(*c, mulii(q, shifti(addii(*b, r),-1)));
616 *b = r;
617 *u2 = subii(*u2, mulii(q, u1));
618 }
619
620 /* q t_QFI, return reduced representative and set base change U in Sl2(Z) */
621 GEN
redimagsl2(GEN q,GEN * U)622 redimagsl2(GEN q, GEN *U)
623 {
624 GEN Q = cgetg(4, t_QFI);
625 pari_sp av = avma, av2;
626 GEN z, u1,u2,v1,v2, a = gel(q,1), b = gel(q,2), c = gel(q,3);
627 long cmp;
628 /* upper bound for size of final (a,b,c) */
629 (void)new_chunk(2*(lgefint(a) + lgefint(b) + lgefint(c) + 3));
630 av2 = avma;
631 u1 = gen_1; u2 = gen_0;
632 cmp = abscmpii(a, b);
633 if (cmp < 0)
634 REDBU(a,&b,&c, u1,&u2);
635 else if (cmp == 0 && signe(b) < 0)
636 { /* b = -a */
637 b = negi(b);
638 u2 = gen_1;
639 }
640 for(;;)
641 {
642 cmp = abscmpii(a, c); if (cmp <= 0) break;
643 swap(a,c); b = negi(b);
644 z = u1; u1 = u2; u2 = negi(z);
645 REDBU(a,&b,&c, u1,&u2);
646 if (gc_needed(av, 1)) {
647 if (DEBUGMEM>1) pari_warn(warnmem, "redimagsl2");
648 gerepileall(av2, 5, &a,&b,&c, &u1,&u2);
649 }
650 }
651 if (cmp == 0 && signe(b) < 0)
652 {
653 b = negi(b);
654 z = u1; u1 = u2; u2 = negi(z);
655 }
656 set_avma(av);
657 a = icopy(a); gel(Q,1) = a;
658 b = icopy(b); gel(Q,2) = b;
659 c = icopy(c); gel(Q,3) = c;
660 u1 = icopy(u1);
661 u2 = icopy(u2); av = avma;
662
663 /* Let q = (A,B,C). q o [u1,u2; v1,v2] = Q implies
664 * [v1,v2] = (1/C) [(b-B)/2 u1 - a u2, c u1 - (b+B)/2 u2] */
665 z = shifti(subii(b, gel(q,2)), -1);
666 v1 = subii(mulii(z, u1), mulii(a, u2)); v1 = diviiexact(v1, gel(q,3));
667 z = subii(z, b);
668 v2 = addii(mulii(z, u2), mulii(c, u1)); v2 = diviiexact(v2, gel(q,3));
669 set_avma(av);
670 v1 = icopy(v1);
671 v2 = icopy(v2);
672 *U = mkmat2(mkcol2(u1,v1), mkcol2(u2,v2)); return Q;
673 }
674
675 static GEN
setq_b0(ulong a,ulong c)676 setq_b0(ulong a, ulong c)
677 { retmkqfi( utoipos(a), gen_0, utoipos(c) ); }
678 /* assume |sb| = 1 */
679 static GEN
setq(ulong a,ulong b,ulong c,long sb)680 setq(ulong a, ulong b, ulong c, long sb)
681 { retmkqfi( utoipos(a), sb == 1? utoipos(b): utoineg(b), utoipos(c) ); }
682 /* 0 < a, c < 2^BIL, b = 0 */
683 static GEN
redimag_1_b0(ulong a,ulong c)684 redimag_1_b0(ulong a, ulong c)
685 { return (a <= c)? setq_b0(a, c): setq_b0(c, a); }
686
687 /* 0 < a, c < 2^BIL: single word affair */
688 static GEN
redimag_1(pari_sp av,GEN a,GEN b,GEN c)689 redimag_1(pari_sp av, GEN a, GEN b, GEN c)
690 {
691 ulong ua, ub, uc;
692 long sb;
693 for(;;)
694 { /* at most twice */
695 long lb = lgefint(b); /* <= 3 after first loop */
696 if (lb == 2) return redimag_1_b0(a[2],c[2]);
697 if (lb == 3 && uel(b,2) <= (ulong)LONG_MAX) break;
698 REDB(a,&b,&c);
699 if (uel(a,2) <= uel(c,2))
700 { /* lg(b) <= 3 but may be too large for itos */
701 long s = signe(b);
702 set_avma(av);
703 if (!s) return redimag_1_b0(a[2], c[2]);
704 if (a[2] == c[2]) s = 1;
705 return setq(a[2], b[2], c[2], s);
706 }
707 swap(a,c); b = negi(b);
708 }
709 /* b != 0 */
710 set_avma(av);
711 ua = a[2];
712 ub = sb = b[2]; if (signe(b) < 0) sb = -sb;
713 uc = c[2];
714 if (ua < ub)
715 sREDB(ua, &sb, &uc);
716 else if (ua == ub && sb < 0) sb = (long)ub;
717 while(ua > uc)
718 {
719 lswap(ua,uc); sb = -sb;
720 sREDB(ua, &sb, &uc);
721 }
722 if (!sb) return setq_b0(ua, uc);
723 else
724 {
725 long s = 1;
726 if (sb < 0)
727 {
728 sb = -sb;
729 if (ua != uc) s = -1;
730 }
731 return setq(ua, sb, uc, s);
732 }
733 }
734
735 static GEN
redimag_av(pari_sp av,GEN q)736 redimag_av(pari_sp av, GEN q)
737 {
738 GEN a = gel(q,1), b = gel(q,2), c = gel(q,3);
739 long cmp, lc = lgefint(c);
740
741 if (lgefint(a) == 3 && lc == 3) return redimag_1(av, a, b, c);
742 cmp = abscmpii(a, b);
743 if (cmp < 0)
744 REDB(a,&b,&c);
745 else if (cmp == 0 && signe(b) < 0)
746 b = negi(b);
747 for(;;)
748 {
749 cmp = abscmpii(a, c); if (cmp <= 0) break;
750 lc = lgefint(a); /* lg(future c): we swap a & c next */
751 if (lc == 3) return redimag_1(av, a, b, c);
752 swap(a,c); b = negi(b); /* apply rho */
753 REDB(a,&b,&c);
754 }
755 if (cmp == 0 && signe(b) < 0) b = negi(b);
756 /* size of reduced Qfb(a,b,c) <= 3 lg(c) + 4 <= 4 lg(c) */
757 (void)new_chunk(lc<<2);
758 a = icopy(a); b = icopy(b); c = icopy(c);
759 set_avma(av);
760 retmkqfi(icopy(a), icopy(b), icopy(c));
761 }
762 GEN
redimag(GEN q)763 redimag(GEN q) { return redimag_av(avma, q); }
764
765 static GEN
rhoimag(GEN x)766 rhoimag(GEN x)
767 {
768 GEN a = gel(x,1), b = gel(x,2), c = gel(x,3);
769 int fl = abscmpii(a, c);
770 if (fl <= 0) {
771 int fg = abscmpii(a, b);
772 if (fg >= 0) {
773 x = qfi(a,b,c);
774 if ((!fl || !fg) && signe(gel(x,2)) < 0) setsigne(gel(x,2), 1);
775 return x;
776 }
777 }
778 x = cgetg(4, t_QFI);
779 (void)new_chunk(lgefint(a) + lgefint(b) + lgefint(c) + 3);
780 swap(a,c); b = negi(b);
781 REDB(a, &b, &c); set_avma((pari_sp)x);
782 gel(x,1) = icopy(a);
783 gel(x,2) = icopy(b);
784 gel(x,3) = icopy(c); return x;
785 }
786
787 /* qfr3 / qfr5 */
788
789 /* t_QFR are unusable: D, sqrtD, isqrtD are recomputed all the time and the
790 * logarithmic Shanks's distance is costly and hard to control.
791 * qfr3 / qfr5 routines take a container of t_INTs (e.g a t_VEC) as argument,
792 * at least 3 (resp. 5) components [it is a feature that they do not check the
793 * precise type or length of the input]. They return a vector of length 3
794 * (resp. 5). A qfr3 [a,b,c] contains the form coeffs, in a qfr5 [a,b,c, e,d]
795 * the t_INT e is a binary exponent, d a t_REAL, coding the distance in
796 * multiplicative form: the true distance is obtained from qfr5_dist.
797 * All other qfr routines are obsolete (inefficient) wrappers */
798
799 /* static functions are not stack-clean. Unless mentionned otherwise, public
800 * functions are. */
801
802 #define EMAX 22
803 static void
fix_expo(GEN x)804 fix_expo(GEN x)
805 {
806 if (expo(gel(x,5)) >= (1L << EMAX)) {
807 gel(x,4) = addiu(gel(x,4), 1);
808 shiftr_inplace(gel(x,5), - (1L << EMAX));
809 }
810 }
811
812 /* (1/2) log (d * 2^{e * 2^EMAX}). Not stack clean if e != 0 */
813 GEN
qfr5_dist(GEN e,GEN d,long prec)814 qfr5_dist(GEN e, GEN d, long prec)
815 {
816 GEN t = logr_abs(d);
817 if (signe(e)) {
818 GEN u = mulir(e, mplog2(prec));
819 shiftr_inplace(u, EMAX); t = addrr(t, u);
820 }
821 shiftr_inplace(t, -1); return t;
822 }
823
824 static void
rho_get_BC(GEN * B,GEN * C,GEN b,GEN c,struct qfr_data * S)825 rho_get_BC(GEN *B, GEN *C, GEN b, GEN c, struct qfr_data *S)
826 {
827 GEN t, u;
828 u = shifti(c,1);
829 t = (abscmpii(S->isqrtD,c) >= 0)? S->isqrtD: c;
830 u = remii(addii_sign(t,1, b,signe(b)), u);
831 *B = addii_sign(t, 1, u, -signe(u)); /* |t| - (|t|+b) % |2c| */
832 if (*B == gen_0)
833 { u = shifti(S->D, -2); setsigne(u, -1); }
834 else
835 u = shifti(addii_sign(sqri(*B),1, S->D,-1), -2);
836 *C = diviiexact(u, c); /* = (B^2-D)/4c */
837 }
838 /* Not stack-clean */
839 GEN
qfr3_rho(GEN x,struct qfr_data * S)840 qfr3_rho(GEN x, struct qfr_data *S)
841 {
842 GEN B, C, b = gel(x,2), c = gel(x,3);
843 rho_get_BC(&B, &C, b, c, S);
844 return mkvec3(c,B,C);
845 }
846 /* Not stack-clean */
847 GEN
qfr5_rho(GEN x,struct qfr_data * S)848 qfr5_rho(GEN x, struct qfr_data *S)
849 {
850 GEN B, C, y, b = gel(x,2), c = gel(x,3);
851 long sb = signe(b);
852
853 rho_get_BC(&B, &C, b, c, S);
854 y = mkvec5(c,B,C, gel(x,4), gel(x,5));
855 if (sb) {
856 GEN t = subii(sqri(b), S->D);
857 if (sb < 0)
858 t = divir(t, sqrr(subir(b,S->sqrtD)));
859 else
860 t = divri(sqrr(addir(b,S->sqrtD)), t);
861 /* t = (b + sqrt(D)) / (b - sqrt(D)), evaluated stably */
862 gel(y,5) = mulrr(t, gel(y,5)); fix_expo(y);
863 }
864 return y;
865 }
866
867 /* Not stack-clean */
868 GEN
qfr_to_qfr5(GEN x,long prec)869 qfr_to_qfr5(GEN x, long prec)
870 { return mkvec5(gel(x,1),gel(x,2),gel(x,3),gen_0,real_1(prec)); }
871
872 /* d0 = initial distance, x = [a,b,c, expo(d), d], d = exp(2*distance) */
873 GEN
qfr5_to_qfr(GEN x,GEN d0)874 qfr5_to_qfr(GEN x, GEN d0)
875 {
876 GEN y;
877 if (lg(x) == 6)
878 {
879 GEN n = gel(x,4), d = absr(gel(x,5));
880 if (signe(n))
881 {
882 n = addis(shifti(n, EMAX), expo(d));
883 setexpo(d, 0); d = logr_abs(d);
884 if (signe(n)) d = addrr(d, mulir(n, mplog2(lg(d0))));
885 shiftr_inplace(d, -1);
886 d0 = addrr(d0, d);
887 }
888 else if (!gequal1(d)) /* avoid loss of precision */
889 {
890 d = logr_abs(d);
891 shiftr_inplace(d, -1);
892 d0 = addrr(d0, d);
893 }
894 }
895 y = cgetg(5, t_QFR);
896 gel(y,1) = gel(x,1);
897 gel(y,2) = gel(x,2);
898 gel(y,3) = gel(x,3);
899 gel(y,4) = d0; return y;
900 }
901
902 /* Not stack-clean */
903 GEN
qfr3_to_qfr(GEN x,GEN d)904 qfr3_to_qfr(GEN x, GEN d)
905 {
906 GEN z = cgetg(5, t_QFR);
907 gel(z,1) = gel(x,1);
908 gel(z,2) = gel(x,2);
909 gel(z,3) = gel(x,3);
910 gel(z,4) = d; return z;
911 }
912
913 static int
ab_isreduced(GEN a,GEN b,GEN isqrtD)914 ab_isreduced(GEN a, GEN b, GEN isqrtD)
915 {
916 if (signe(b) > 0 && abscmpii(b, isqrtD) <= 0)
917 {
918 GEN t = addii_sign(isqrtD,1, shifti(a,1),-1);
919 long l = abscmpii(b, t); /* compare |b| and |floor(sqrt(D)) - |2a|| */
920 if (l > 0 || (l == 0 && signe(t) < 0)) return 1;
921 }
922 return 0;
923 }
924
925 INLINE int
qfr_isreduced(GEN x,GEN isqrtD)926 qfr_isreduced(GEN x, GEN isqrtD)
927 {
928 return ab_isreduced(gel(x,1),gel(x,2),isqrtD);
929 }
930
931 /* Not stack-clean */
932 GEN
qfr5_red(GEN x,struct qfr_data * S)933 qfr5_red(GEN x, struct qfr_data *S)
934 {
935 pari_sp av = avma;
936 while (!qfr_isreduced(x, S->isqrtD))
937 {
938 x = qfr5_rho(x, S);
939 if (gc_needed(av,2))
940 {
941 if (DEBUGMEM>1) pari_warn(warnmem,"qfr5_red");
942 x = gerepilecopy(av, x);
943 }
944 }
945 return x;
946 }
947 /* Not stack-clean */
948 GEN
qfr3_red(GEN x,struct qfr_data * S)949 qfr3_red(GEN x, struct qfr_data *S)
950 {
951 pari_sp av = avma;
952 while (!qfr_isreduced(x, S->isqrtD))
953 {
954 x = qfr3_rho(x, S);
955 if (gc_needed(av,2))
956 {
957 if (DEBUGMEM>1) pari_warn(warnmem,"qfr3_red");
958 x = gerepilecopy(av, x);
959 }
960 }
961 return x;
962 }
963
964 static void
get_disc(GEN x,struct qfr_data * S)965 get_disc(GEN x, struct qfr_data *S)
966 {
967 if (!S->D) S->D = qfb_disc(x);
968 else if (typ(S->D) != t_INT) pari_err_TYPE("qfr_init",S->D);
969 if (!signe(S->D)) pari_err_DOMAIN("qfr_init", "disc", "=", gen_0,x);
970 }
971
972 void
qfr_data_init(GEN D,long prec,struct qfr_data * S)973 qfr_data_init(GEN D, long prec, struct qfr_data *S)
974 {
975 S->D = D;
976 S->sqrtD = sqrtr(itor(S->D,prec));
977 S->isqrtD = truncr(S->sqrtD);
978 }
979
980 static GEN
qfr5_init(GEN x,struct qfr_data * S)981 qfr5_init(GEN x, struct qfr_data *S)
982 {
983 GEN d = gel(x,4);
984 long prec = realprec(d), l = -expo(d);
985 if (l < BITS_IN_LONG) l = BITS_IN_LONG;
986 prec = maxss(prec, nbits2prec(l));
987 x = qfr_to_qfr5(x,prec);
988
989 get_disc(x, S);
990 if (!S->sqrtD) S->sqrtD = sqrtr(itor(S->D,prec));
991 else if (typ(S->sqrtD) != t_REAL) pari_err_TYPE("qfr_init",S->sqrtD);
992
993 if (!S->isqrtD)
994 {
995 pari_sp av=avma;
996 long e;
997 S->isqrtD = gcvtoi(S->sqrtD,&e);
998 if (e>-2) { set_avma(av); S->isqrtD = sqrti(S->D); }
999 }
1000 else if (typ(S->isqrtD) != t_INT) pari_err_TYPE("qfr_init",S->isqrtD);
1001 return x;
1002 }
1003 static GEN
qfr3_init(GEN x,struct qfr_data * S)1004 qfr3_init(GEN x, struct qfr_data *S)
1005 {
1006 get_disc(x, S);
1007 if (!S->isqrtD) S->isqrtD = sqrti(S->D);
1008 else if (typ(S->isqrtD) != t_INT) pari_err_TYPE("qfr_init",S->isqrtD);
1009 return x;
1010 }
1011
1012 #define qf_NOD 2
1013 #define qf_STEP 1
1014
1015 static GEN
redreal0(GEN x,long flag,GEN D,GEN isqrtD,GEN sqrtD)1016 redreal0(GEN x, long flag, GEN D, GEN isqrtD, GEN sqrtD)
1017 {
1018 pari_sp av = avma;
1019 struct qfr_data S;
1020 GEN d;
1021 if (typ(x) != t_QFR) pari_err_TYPE("redreal",x);
1022 d = gel(x,4);
1023 S.D = D;
1024 S.sqrtD = sqrtD;
1025 S.isqrtD = isqrtD;
1026 x = (flag & qf_NOD)? qfr3_init(x, &S): qfr5_init(x, &S);
1027 switch(flag) {
1028 case 0: x = qfr5_red(x,&S); break;
1029 case qf_NOD: x = qfr3_red(x,&S); break;
1030 case qf_STEP: x = qfr5_rho(x,&S); break;
1031 case qf_STEP|qf_NOD: x = qfr3_rho(x,&S); break;
1032 default: pari_err_FLAG("qfbred");
1033 }
1034 return gerepilecopy(av, qfr5_to_qfr(x,d));
1035 }
1036 GEN
redreal(GEN x)1037 redreal(GEN x)
1038 { return redreal0(x,0,NULL,NULL,NULL); }
1039 GEN
rhoreal(GEN x)1040 rhoreal(GEN x)
1041 { return redreal0(x,qf_STEP,NULL,NULL,NULL); }
1042 GEN
redrealnod(GEN x,GEN isqrtD)1043 redrealnod(GEN x, GEN isqrtD)
1044 { return redreal0(x,qf_NOD,NULL,isqrtD,NULL); }
1045 GEN
rhorealnod(GEN x,GEN isqrtD)1046 rhorealnod(GEN x, GEN isqrtD)
1047 { return redreal0(x,qf_STEP|qf_NOD,NULL,isqrtD,NULL); }
1048 GEN
qfbred0(GEN x,long flag,GEN D,GEN isqrtD,GEN sqrtD)1049 qfbred0(GEN x, long flag, GEN D, GEN isqrtD, GEN sqrtD)
1050 {
1051 if (typ(x) == t_QFI)
1052 return (flag & qf_STEP)? rhoimag(x): redimag(x);
1053 return redreal0(x,flag,D,isqrtD,sqrtD);
1054 }
1055
1056 GEN
qfr5_comp(GEN x,GEN y,struct qfr_data * S)1057 qfr5_comp(GEN x, GEN y, struct qfr_data *S)
1058 {
1059 pari_sp av = avma;
1060 GEN z = cgetg(6,t_VEC); qfb_comp(z,x,y);
1061 if (x == y)
1062 {
1063 gel(z,4) = shifti(gel(x,4),1);
1064 gel(z,5) = sqrr(gel(x,5));
1065 }
1066 else
1067 {
1068 gel(z,4) = addii(gel(x,4),gel(y,4));
1069 gel(z,5) = mulrr(gel(x,5),gel(y,5));
1070 }
1071 fix_expo(z); z = qfr5_red(z,S);
1072 return gerepilecopy(av,z);
1073 }
1074 /* Not stack-clean */
1075 GEN
qfr3_comp(GEN x,GEN y,struct qfr_data * S)1076 qfr3_comp(GEN x, GEN y, struct qfr_data *S)
1077 {
1078 GEN z = cgetg(4,t_VEC); qfb_comp(z,x,y);
1079 return qfr3_red(z, S);
1080 }
1081
1082 /* valid for t_QFR, qfr3, qfr5 */
1083 static GEN
qfr_inv(GEN x)1084 qfr_inv(GEN x) {
1085 GEN z = shallowcopy(x);
1086 gel(z,2) = negi(gel(z,2));
1087 return z;
1088 }
1089
1090 /* return x^n. Not stack-clean */
1091 GEN
qfr5_pow(GEN x,GEN n,struct qfr_data * S)1092 qfr5_pow(GEN x, GEN n, struct qfr_data *S)
1093 {
1094 GEN y = NULL;
1095 long i, m, s = signe(n);
1096 if (!s) return qfr5_1(S, lg(gel(x,5)));
1097 for (i=lgefint(n)-1; i>1; i--)
1098 {
1099 m = n[i];
1100 for (; m; m>>=1)
1101 {
1102 if (m&1) y = y? qfr5_comp(y,x,S): x;
1103 if (m == 1 && i == 2) break;
1104 x = qfr5_comp(x,x,S);
1105 }
1106 }
1107 return y;
1108 }
1109 /* return x^n. Not stack-clean */
1110 GEN
qfr3_pow(GEN x,GEN n,struct qfr_data * S)1111 qfr3_pow(GEN x, GEN n, struct qfr_data *S)
1112 {
1113 GEN y = NULL;
1114 long i, m, s = signe(n);
1115 if (!s) return qfr3_1(S);
1116 if (s < 0) x = qfr_inv(x);
1117 for (i=lgefint(n)-1; i>1; i--)
1118 {
1119 m = n[i];
1120 for (; m; m>>=1)
1121 {
1122 if (m&1) y = y? qfr3_comp(y,x,S): x;
1123 if (m == 1 && i == 2) break;
1124 x = qfr3_comp(x,x,S);
1125 }
1126 }
1127 return y;
1128 }
1129 GEN
qfrpow(GEN x,GEN n)1130 qfrpow(GEN x, GEN n)
1131 {
1132 struct qfr_data S = { NULL, NULL, NULL };
1133 long s = signe(n);
1134 pari_sp av = avma;
1135 GEN d0;
1136
1137 if (!s) return qfr_1(x);
1138 if (is_pm1(n)) return s > 0? redreal(x): ginv(x);
1139 if (s < 0) x = qfr_inv(x);
1140 d0 = gel(x,4);
1141 if (!signe(d0)) {
1142 x = qfr3_init(x, &S);
1143 x = qfr3_pow(x, n, &S);
1144 x = qfr3_to_qfr(x, d0);
1145 } else {
1146 x = qfr5_init(x, &S);
1147 x = qfr5_pow(qfr_to_qfr5(x, lg(S.sqrtD)), n, &S);
1148 x = qfr5_to_qfr(x, mulri(d0,n));
1149 }
1150 return gerepilecopy(av, x);
1151 }
1152
1153 /* Prime forms attached to prime ideals of degree 1 */
1154
1155 /* assume x != 0 a t_INT, p > 0
1156 * Return a t_QFI, but discriminant sign is not checked: can be used for
1157 * real forms as well */
1158 GEN
primeform_u(GEN x,ulong p)1159 primeform_u(GEN x, ulong p)
1160 {
1161 GEN c, y = cgetg(4, t_QFI);
1162 pari_sp av = avma;
1163 ulong b;
1164 long s;
1165
1166 s = mod8(x); if (signe(x) < 0 && s) s = 8-s;
1167 /* 2 or 3 mod 4 */
1168 if (s & 2) pari_err_DOMAIN("primeform", "disc % 4", ">",gen_1, x);
1169 if (p == 2) {
1170 switch(s) {
1171 case 0: b = 0; break;
1172 case 1: b = 1; break;
1173 case 4: b = 2; break;
1174 default: pari_err_SQRTN("primeform", mkintmod(x,utoi(p)) );
1175 b = 0; /* -Wall */
1176 }
1177 c = shifti(subsi(s,x), -3);
1178 } else {
1179 b = Fl_sqrt(umodiu(x,p), p);
1180 if (b == ~0UL) pari_err_SQRTN("primeform", mkintmod(x,utoi(p)) );
1181 /* mod(b) != mod2(x) ? */
1182 if ((b ^ s) & 1) b = p - b;
1183 c = diviuexact(shifti(subii(sqru(b), x), -2), p);
1184 }
1185 gel(y,3) = gerepileuptoint(av, c);
1186 gel(y,2) = utoi(b);
1187 gel(y,1) = utoipos(p); return y;
1188 }
1189
1190 /* special case: p = 1 return unit form */
1191 GEN
primeform(GEN x,GEN p,long prec)1192 primeform(GEN x, GEN p, long prec)
1193 {
1194 const char *f = "primeform";
1195 pari_sp av;
1196 long s, sx = signe(x), sp = signe(p);
1197 GEN y, b, absp;
1198
1199 if (typ(x) != t_INT) pari_err_TYPE(f,x);
1200 if (typ(p) != t_INT) pari_err_TYPE(f,p);
1201 if (!sp) pari_err_DOMAIN(f,"p","=",gen_0,p);
1202 if (!sx) pari_err_DOMAIN(f,"D","=",gen_0,x);
1203 if (lgefint(p) == 3)
1204 {
1205 ulong pp = p[2];
1206 if (pp == 1) {
1207 if (sx < 0) {
1208 long r;
1209 if (sp < 0) pari_err_IMPL("negative definite t_QFI");
1210 r = mod4(x);
1211 if (r && r != 3) pari_err_DOMAIN(f,"disc % 4",">", gen_1,x);
1212 return qfi_1_by_disc(x);
1213 }
1214 y = qfr_1_by_disc(x,prec);
1215 if (sp < 0) { gel(y,1) = gen_m1; togglesign(gel(y,3)); }
1216 return y;
1217 }
1218 y = primeform_u(x, pp);
1219 if (sx < 0) {
1220 if (sp < 0) pari_err_IMPL("negative definite t_QFI");
1221 return y;
1222 }
1223 if (sp < 0) { togglesign(gel(y,1)); togglesign(gel(y,3)); }
1224 return gcopy( qfr3_to_qfr(y, real_0(prec)) );
1225 }
1226 s = mod8(x);
1227 if (sx < 0)
1228 {
1229 if (sp < 0) pari_err_IMPL("negative definite t_QFI");
1230 if (s) s = 8-s;
1231 y = cgetg(4, t_QFI);
1232 }
1233 else
1234 {
1235 y = cgetg(5, t_QFR);
1236 gel(y,4) = real_0(prec);
1237 }
1238 /* 2 or 3 mod 4 */
1239 if (s & 2) pari_err_DOMAIN(f, "disc % 4", ">",gen_1, x);
1240 absp = absi_shallow(p); av = avma;
1241 b = Fp_sqrt(x, absp); if (!b) pari_err_SQRTN(f, mkintmod(x,absp));
1242 s &= 1; /* s = x mod 2 */
1243 /* mod(b) != mod2(x) ? [Warning: we may have b == 0] */
1244 if ((!signe(b) && s) || mod2(b) != s) b = gerepileuptoint(av, subii(absp,b));
1245
1246 av = avma;
1247 gel(y,3) = gerepileuptoint(av, diviiexact(shifti(subii(sqri(b), x), -2), p));
1248 gel(y,2) = b;
1249 gel(y,1) = icopy(p);
1250 return y;
1251 }
1252
1253 static GEN
normforms(GEN D,GEN fa)1254 normforms(GEN D, GEN fa)
1255 {
1256 long i, j, k, lB, aN, sa;
1257 GEN a, L, V, B, N, N2;
1258 int D_odd = mpodd(D);
1259 a = typ(fa) == t_INT ? fa: typ(fa) == t_VEC? gel(fa,1): factorback(fa);
1260 sa = signe(a);
1261 if (sa==0 || (signe(D)<0 && sa<0)) return NULL;
1262 V = D_odd? Zn_quad_roots(fa, gen_1, shifti(subsi(1, D), -2))
1263 : Zn_quad_roots(fa, gen_0, negi(shifti(D, -2)));
1264 if (!V) return NULL;
1265 N = gel(V,1); B = gel(V,2); lB = lg(B);
1266 N2 = shifti(N,1);
1267 aN = itou(diviiexact(a, N)); /* |a|/N */
1268 L = cgetg((lB-1)*aN+1, t_VEC);
1269 for (k = 1, i = 1; i < lB; i++)
1270 {
1271 GEN b = shifti(gel(B,i), 1), c, C;
1272 if (D_odd) b = addiu(b , 1);
1273 c = diviiexact(shifti(subii(sqri(b), D), -2), a);
1274 for (j = 0;; b = addii(b, N2))
1275 {
1276 gel(L, k++) = mkvec3(a,b,c);
1277 if (++j == aN) break;
1278 C = addii(b, N); if (aN > 1) C = diviuexact(C, aN);
1279 c = sa > 0? addii(c, C): subii(c, C);
1280 }
1281 }
1282 return L;
1283 }
1284
1285 /* Let M and N in SL2(Z), return (N*M^-1)[,1] */
1286 static GEN
SL2_div_mul_e1(GEN N,GEN M)1287 SL2_div_mul_e1(GEN N, GEN M)
1288 {
1289 GEN b = gcoeff(M,2,1), d = gcoeff(M,2,2);
1290 GEN p = subii(mulii(gcoeff(N,1,1), d), mulii(gcoeff(N,1,2), b));
1291 GEN q = subii(mulii(gcoeff(N,2,1), d), mulii(gcoeff(N,2,2), b));
1292 return mkvec2(p,q);
1293 }
1294
1295 /* Let M and N in SL2(Z), return (N*[1,0;0,-1]*M^-1)[,1] */
1296 static GEN
SL2_swap_div_mul_e1(GEN N,GEN M)1297 SL2_swap_div_mul_e1(GEN N, GEN M)
1298 {
1299 GEN b = gcoeff(M,2,1), d = gcoeff(M,2,2);
1300 GEN p = addii(mulii(gcoeff(N,1,1), d), mulii(gcoeff(N,1,2), b));
1301 GEN q = addii(mulii(gcoeff(N,2,1), d), mulii(gcoeff(N,2,2), b));
1302 return mkvec2(p,q);
1303 }
1304
1305 /* Test equality modulo GL2 of two reduced forms */
1306 static int
GL2_qfb_equal(GEN a,GEN b)1307 GL2_qfb_equal(GEN a, GEN b)
1308 {
1309 return equalii(gel(a,1),gel(b,1))
1310 && absequalii(gel(a,2),gel(b,2))
1311 && equalii(gel(a,3),gel(b,3));
1312 }
1313
1314 static GEN
qfisolve_normform(GEN Q,GEN P)1315 qfisolve_normform(GEN Q, GEN P)
1316 {
1317 GEN a = gel(Q,1), N = gel(Q,2);
1318 GEN M, b = redimagsl2(P, &M);
1319 if (!GL2_qfb_equal(a,b) || signe(gel(a,2))!=signe(gel(b,2)))
1320 return NULL;
1321 return SL2_div_mul_e1(N,M);
1322 }
1323
1324 static GEN
qfbsolve_cornacchia(GEN c,GEN p,int swap)1325 qfbsolve_cornacchia(GEN c, GEN p, int swap)
1326 {
1327 pari_sp av = avma;
1328 GEN M, N;
1329 if (kronecker(negi(c), p) < 0 || !cornacchia(c, p, &M,&N))
1330 { set_avma(av); return gen_0; }
1331 return gerepilecopy(av, swap? mkvec2(N,M): mkvec2(M,N));
1332 }
1333
1334 GEN
qfisolvep(GEN Q,GEN p)1335 qfisolvep(GEN Q, GEN p)
1336 {
1337 GEN M, N, x,y, a,b,c, d;
1338 pari_sp av = avma;
1339 if (!signe(gel(Q,2)))
1340 {
1341 a = gel(Q,1);
1342 c = gel(Q,3); /* if principal form, use faster cornacchia */
1343 if (equali1(a)) return qfbsolve_cornacchia(c, p, 0);
1344 if (equali1(c)) return qfbsolve_cornacchia(a, p, 1);
1345 }
1346 d = qfb_disc(Q); if (kronecker(d,p) < 0) return gen_0;
1347 a = redimagsl2(Q, &N);
1348 if (equali1(gel(a,1))) /* principal form */
1349 {
1350 long r;
1351 if (!signe(gel(a,2)))
1352 {
1353 a = qfbsolve_cornacchia(gel(a,3), p, 0);
1354 if (a == gen_0) { set_avma(av); return gen_0; }
1355 a = ZM_ZC_mul(N, a);
1356 a[0] = evaltyp(t_VEC) | _evallg(3); /* transpose */
1357 return gerepileupto(av, a);
1358 }
1359 /* x^2 + xy + ((1-d)/4)y^2 = p <==> (2x + y)^2 - d y^2 = 4p */
1360 if (!cornacchia2(negi(d), p, &x, &y)) { set_avma(av); return gen_0; }
1361 x = divis_rem(subii(x,y), 2, &r); if (r) { set_avma(av); return gen_0; }
1362 a = ZM_ZC_mul(N, mkvec2(x,y));
1363 a[0] = evaltyp(t_VEC) | _evallg(3); /* transpose */
1364 return gerepileupto(av, a);
1365 }
1366 b = redimagsl2(primeform(d, p, 0), &M);
1367 if (!GL2_qfb_equal(a,b)) { set_avma(av); return gen_0; }
1368 if (signe(gel(a,2))==signe(gel(b,2)))
1369 x = SL2_div_mul_e1(N,M);
1370 else
1371 x = SL2_swap_div_mul_e1(N,M);
1372 return gerepilecopy(av, x);
1373 }
1374
1375 GEN
redrealsl2step(GEN A,GEN d,GEN rd)1376 redrealsl2step(GEN A, GEN d, GEN rd)
1377 {
1378 pari_sp ltop = avma;
1379 GEN N, V = gel(A,1), M = gel(A,2);
1380 GEN a = gel(V,1), b = gel(V,2), c = gel(V,3);
1381 GEN C = mpabs_shallow(c);
1382 GEN t = addii(b, gmax_shallow(rd, C));
1383 GEN r, q = truedvmdii(t, shifti(C,1), &r);
1384 b = subii(t, addii(r,b));
1385 a = c;
1386 c = truedivii(subii(sqri(b), d), shifti(c,2));
1387 if (signe(a) < 0) togglesign(q);
1388 N = mkmat2(gel(M,2),
1389 mkcol2(subii(mulii(q, gcoeff(M, 1, 2)), gcoeff(M, 1, 1)),
1390 subii(mulii(q, gcoeff(M, 2, 2)), gcoeff(M, 2, 1))));
1391 return gerepilecopy(ltop, mkvec2(mkvec3(a,b,c),N));
1392 }
1393
1394 GEN
redrealsl2(GEN V,GEN d,GEN rd)1395 redrealsl2(GEN V, GEN d, GEN rd)
1396 {
1397 pari_sp ltop = avma;
1398 GEN M, u1, u2, v1, v2;
1399 GEN a = gel(V,1), b = gel(V,2), c = gel(V,3);
1400 u1 = v2 = gen_1; v1 = u2 = gen_0;
1401 while (!ab_isreduced(a,b,rd))
1402 {
1403 GEN C = mpabs_shallow(c);
1404 GEN t = addii(b, gmax_shallow(rd,C));
1405 GEN r, q = truedvmdii(t, shifti(C,1), &r);
1406 b = subii(t, addii(r,b));
1407 a = c;
1408 c = truedivii(subii(sqri(b), d), shifti(c,2));
1409 if (signe(a) < 0) togglesign(q);
1410 r = u1; u1 = v1; v1 = subii(mulii(q, v1), r);
1411 r = u2; u2 = v2; v2 = subii(mulii(q, v2), r);
1412 if (gc_needed(ltop, 1))
1413 {
1414 if (DEBUGMEM>1) pari_warn(warnmem,"redrealsl2");
1415 gerepileall(ltop, 7, &a,&b,&c,&u1,&u2,&v1,&v2);
1416 }
1417 }
1418 M = mkmat2(mkcol2(u1,u2), mkcol2(v1,v2));
1419 return gerepilecopy(ltop, mkvec2(mkvec3(a,b,c), M));
1420 }
1421
1422 GEN
qfbredsl2(GEN q,GEN S)1423 qfbredsl2(GEN q, GEN S)
1424 {
1425 GEN v, D, isD;
1426 pari_sp av;
1427 switch(typ(q))
1428 {
1429 case t_QFI:
1430 if (S) pari_err_TYPE("qfbredsl2",S);
1431 v = cgetg(3,t_VEC);
1432 gel(v,1) = redimagsl2(q, &gel(v,2));
1433 return v;
1434 case t_QFR:
1435 av = avma;
1436 if (S) {
1437 if (typ(S) != t_VEC || lg(S) != 3) pari_err_TYPE("qfbredsl2",S);
1438 D = gel(S,1);
1439 isD = gel(S,2);
1440 if (typ(D) != t_INT || signe(D) <= 0 || typ(isD) != t_INT)
1441 pari_err_TYPE("qfbredsl2",S);
1442 }
1443 else
1444 {
1445 D = qfb_disc(q);
1446 isD = sqrti(D);
1447 }
1448 v = redrealsl2(q,D,isD);
1449 gel(v,1) = qfr3_to_qfr(gel(v,1), real_0(precision(gel(q,4))));
1450 return gerepilecopy(av, v);
1451
1452 default:
1453 pari_err_TYPE("qfbredsl2",q);
1454 return NULL;/*LCOV_EXCL_LINE*/
1455 }
1456 }
1457
1458 static GEN
qfrsolve_normform(GEN N,GEN Ps,GEN d,GEN rd)1459 qfrsolve_normform(GEN N, GEN Ps, GEN d, GEN rd)
1460 {
1461 pari_sp av = avma, btop;
1462 GEN M = N, P = redrealsl2(Ps, d,rd), Q = P;
1463 btop = avma;
1464 for(;;)
1465 {
1466 if (ZV_equal(gel(M,1), gel(P,1)))
1467 return gerepilecopy(av, SL2_div_mul_e1(gel(M,2),gel(P,2)));
1468 if (ZV_equal(gel(N,1), gel(Q,1)))
1469 return gerepilecopy(av, SL2_div_mul_e1(gel(N,2),gel(Q,2)));
1470 M = redrealsl2step(M, d,rd);
1471 Q = redrealsl2step(Q, d,rd);
1472 if (ZV_equal(gel(M,1), gel(N,1))) return gc_NULL(av);
1473 if (ZV_equal(gel(P,1), gel(Q,1))) return gc_NULL(av);
1474 if (gc_needed(btop, 1))
1475 gerepileall(btop, 2, &M, &Q);
1476 }
1477 }
1478
1479 GEN
qfrsolvep(GEN Q,GEN p)1480 qfrsolvep(GEN Q, GEN p)
1481 {
1482 pari_sp ltop = avma;
1483 GEN P, N, x, rd, d = qfb_disc(Q);
1484 if (kronecker(d, p) < 0) { set_avma(ltop); return gen_0; }
1485 P = primeform(d, p, DEFAULTPREC);
1486 rd = sqrti(d);
1487 N = redrealsl2(Q, d,rd);
1488 x = qfrsolve_normform(N, P, d,rd);
1489 if (!x) { set_avma(ltop); return gen_0; }
1490 return gerepileupto(ltop, x);
1491 }
1492
1493 static GEN
redsl2(GEN Q,GEN d)1494 redsl2(GEN Q, GEN d)
1495 {
1496 GEN P, U;
1497 if( signe(d)>0) return redrealsl2(Q, d, sqrti(d));
1498 P = redimagsl2(Q, &U); return mkvec2(P,U);
1499 }
1500
1501 static GEN
qfsolve_normform(GEN Q,GEN f,GEN d,GEN rd)1502 qfsolve_normform(GEN Q, GEN f, GEN d, GEN rd)
1503 { return rd? qfrsolve_normform(Q, f, d, rd): qfisolve_normform(Q, f); }
1504 static GEN
qfbsolve_primitive_i(GEN Q,GEN d,GEN rd,GEN * Qr,GEN fa,long all)1505 qfbsolve_primitive_i(GEN Q, GEN d, GEN rd, GEN *Qr, GEN fa, long all)
1506 {
1507 GEN x, W, F = normforms(d, fa);
1508 long i, j, l;
1509 if (!F) return NULL;
1510 if (!*Qr) *Qr = redsl2(Q, d);
1511 l = lg(F); W = all? cgetg(l, t_VEC): NULL;
1512 for (j = i = 1; i < l; i++)
1513 if ((x = qfsolve_normform(*Qr, gel(F,i), d, rd)))
1514 {
1515 if (!all) return x;
1516 gel(W,j++) = x;
1517 }
1518 if (j == 1) return NULL;
1519 setlg(W,j); return W;
1520 }
1521
1522 static void
qfb_initrd(GEN Q,GEN * pd,GEN * prd)1523 qfb_initrd(GEN Q, GEN *pd, GEN *prd)
1524 {
1525 *pd = qfb_disc(Q);
1526 *prd = signe(*pd) > 0? sqrti(*pd): NULL;
1527 }
1528 static GEN
qfbsolve_primitive(GEN Q,GEN fa,long all)1529 qfbsolve_primitive(GEN Q, GEN fa, long all)
1530 {
1531 pari_sp av = avma;
1532 GEN Qr = NULL, dQ, rdQ, x;
1533 qfb_initrd(Q, &dQ, &rdQ);
1534 x = qfbsolve_primitive_i(Q, dQ, rdQ, &Qr, fa, all);
1535 if (!x) { set_avma(av); return cgetg(1, t_VEC); }
1536 return gerepilecopy(av, x);
1537 }
1538
1539 /* f / g^2 */
1540 static GEN
famat_divsqr(GEN f,GEN g)1541 famat_divsqr(GEN f, GEN g)
1542 { return famat_reduce(famat_div_shallow(f, famat_pows_shallow(g,2))); }
1543 static GEN
qfbsolve_all(GEN Q,GEN n,long all)1544 qfbsolve_all(GEN Q, GEN n, long all)
1545 {
1546 pari_sp av = avma;
1547 GEN W, dQ, rdQ, Qr = NULL, fa = factorint(n, 0);
1548 GEN D = divisors_factored(mkmat2(gel(fa,1), gshift(gel(fa,2),-1)));
1549 long i, j, l = lg(D);
1550 qfb_initrd(Q, &dQ, &rdQ);
1551 W = all? cgetg(l, t_VEC): NULL;
1552 for (i = j = 1; i < l; i++)
1553 {
1554 GEN w, d = gel(D,i), FA = i == 1? fa: famat_divsqr(fa, gel(d,2));
1555 if ((w = qfbsolve_primitive_i(Q, dQ, rdQ, &Qr, FA, all)))
1556 {
1557 if (i != 1) w = RgV_Rg_mul(w, gel(d,1));
1558 if (!all) return gerepilecopy(av, w);
1559 gel(W,j++) = w;
1560 }
1561 }
1562 if (j == 1) { set_avma(av); return cgetg(1, t_VEC); }
1563 setlg(W,j); return gerepilecopy(av, shallowconcat1(W));
1564 }
1565
1566 GEN
qfbsolve(GEN Q,GEN n,long fl)1567 qfbsolve(GEN Q, GEN n, long fl)
1568 {
1569 if (!is_qfb_t(typ(Q))) pari_err_TYPE("qfbsolve",Q);
1570 if (fl < 0 || fl > 3) pari_err_FLAG("qfbsolve");
1571 return (fl & 2)? qfbsolve_all(Q, n, fl & 1)
1572 : qfbsolve_primitive(Q, n, fl & 1);
1573 }
1574
1575 /* 1 if there exists x,y such that x^2 + dy^2 = p [prime], 0 otherwise */
1576 long
cornacchia(GEN d,GEN p,GEN * px,GEN * py)1577 cornacchia(GEN d, GEN p, GEN *px, GEN *py)
1578 {
1579 pari_sp av = avma;
1580 GEN a, b, c, r;
1581
1582 if (typ(d) != t_INT) pari_err_TYPE("cornacchia", d);
1583 if (typ(p) != t_INT) pari_err_TYPE("cornacchia", p);
1584 if (signe(d) <= 0) pari_err_DOMAIN("cornacchia", "d","<=",gen_0,d);
1585 *px = *py = gen_0;
1586 b = subii(p, d);
1587 if (signe(b) < 0) return 0;
1588 if (signe(b) == 0) { *py = gen_1; return gc_long(av,1); }
1589 b = Fp_sqrt(b, p); /* sqrt(-d) */
1590 if (!b) return gc_long(av,0);
1591 b = gmael(halfgcdii(p, b), 2, 2);
1592 a = subii(p, sqri(b));
1593 c = dvmdii(a, d, &r);
1594 if (r != gen_0 || !Z_issquareall(c, &c)) return gc_long(av,0);
1595 set_avma(av);
1596 *px = icopy(b);
1597 *py = icopy(c); return 1;
1598 }
1599
1600 static GEN
lastqi(GEN Q)1601 lastqi(GEN Q)
1602 {
1603 GEN s = gcoeff(Q,1,1), q = gcoeff(Q,1,2), p = absi_shallow(gcoeff(Q,2,2));
1604 if (signe(q)==0) return gen_0;
1605 if (signe(s)==0) return p;
1606 if (is_pm1(q)) return subiu(p,1);
1607 return divii(p, absi_shallow(q));
1608 }
1609
1610 static long
cornacchia2_helper(long av,GEN d,GEN p,GEN b,GEN px4,GEN * px,GEN * py)1611 cornacchia2_helper(long av, GEN d, GEN p, GEN b, GEN px4, GEN *px, GEN *py)
1612 {
1613 GEN M, Q, V, a, c, r;
1614 if (!signe(b)) { /* d = p,2p,3p,4p */
1615 set_avma(av);
1616 if (absequalii(d, px4)){ *py = gen_1; return 1; }
1617 if (absequalii(d, p)) { *py = gen_2; return 1; }
1618 return 0;
1619 }
1620 if (mod2(b) != mod2(d)) b = subii(p,b);
1621 M = halfgcdii(shifti(p,1), b); Q = gel(M,1); V = gel(M, 2);
1622 b = addii(mulii(gel(V,1), lastqi(Q)),gel(V,2));
1623 if (cmpii(sqri(b),px4) > 0) b = gel(V,1);
1624 if (cmpii(sqri(b),px4) > 0) b = gel(V,2);
1625 a = subii(px4, sqri(b));
1626 c = dvmdii(a, d, &r);
1627 if (r != gen_0 || !Z_issquareall(c, &c)) return gc_long(av,0);
1628 set_avma(av);
1629 *px = icopy(b);
1630 *py = icopy(c); return 1;
1631 }
1632
1633 /* 1 if there exists x,y such that x^2 + dy^2 = 4p [p prime], 0 otherwise */
1634 long
cornacchia2(GEN d,GEN p,GEN * px,GEN * py)1635 cornacchia2(GEN d, GEN p, GEN *px, GEN *py)
1636 {
1637 pari_sp av = avma;
1638 GEN b, px4;
1639 long k;
1640
1641 if (typ(d) != t_INT) pari_err_TYPE("cornacchia2", d);
1642 if (typ(p) != t_INT) pari_err_TYPE("cornacchia2", p);
1643 if (signe(d) <= 0) pari_err_DOMAIN("cornacchia2", "d","<=",gen_0,d);
1644 *px = *py = gen_0;
1645 k = mod4(d);
1646 if (k == 1 || k == 2) pari_err_DOMAIN("cornacchia2","-d mod 4", ">",gen_1,d);
1647 px4 = shifti(p,2);
1648 if (abscmpii(px4, d) < 0) return gc_long(av,0);
1649 if (absequaliu(p, 2))
1650 {
1651 set_avma(av);
1652 switch (itou_or_0(d)) {
1653 case 4: *px = gen_2; break;
1654 case 7: *px = gen_1; break;
1655 default: return 0;
1656 }
1657 *py = gen_1; return 1;
1658 }
1659 b = Fp_sqrt(negi(d), p);
1660 if (!b) return gc_long(av,0);
1661 return cornacchia2_helper(av, d, p, b, px4, px, py);
1662 }
1663
1664 /* 1 if there exists x,y such that x^2 + dy^2 = 4p [p prime], 0 otherwise */
1665 long
cornacchia2_sqrt(GEN d,GEN p,GEN b,GEN * px,GEN * py)1666 cornacchia2_sqrt(GEN d, GEN p, GEN b, GEN *px, GEN *py)
1667 {
1668 pari_sp av = avma;
1669 GEN px4;
1670 *px = *py = gen_0;
1671 px4 = shifti(p,2);
1672 if (abscmpii(px4, d) < 0) return gc_long(av,0);
1673 return cornacchia2_helper(av, d, p, b, px4, px, py);
1674 }
1675