Line data Source code
1 : /* Copyright (C) 2010 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation; either version 2 of the License, or (at your option) any later
8 : version. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 :
15 : /********************************************************************/
16 : /** **/
17 : /** L functions of elliptic curves **/
18 : /** **/
19 : /********************************************************************/
20 : #include "pari.h"
21 : #include "paripriv.h"
22 :
23 : #define DEBUGLEVEL DEBUGLEVEL_ellanal
24 :
25 : struct baby_giant
26 : {
27 : GEN baby, giant, sum;
28 : GEN bnd, rbnd;
29 : };
30 :
31 : /* Generic Buhler-Gross algorithm */
32 :
33 : struct bg_data
34 : {
35 : GEN E, N; /* ell, conductor */
36 : GEN bnd; /* t_INT; will need all an for n <= bnd */
37 : ulong rootbnd; /* sqrt(bnd) */
38 : GEN an; /* t_VECSMALL: cache of an, n <= rootbnd */
39 : GEN p; /* t_VECSMALL: primes <= rootbnd */
40 : };
41 :
42 : typedef void bg_fun(void*el, GEN n, GEN a);
43 :
44 : /* a = a_n, where p = bg->pp[i] divides n, and lasta = a_{n/p}.
45 : * Call fun(E, N, a_N), for all N, n | N, P^+(N) <= p, a_N != 0,
46 : * i.e. assumes that fun accumulates a_N * w(N) */
47 :
48 : static void
49 3431482 : gen_BG_add(void *E, bg_fun *fun, struct bg_data *bg, GEN n, long i, GEN a, GEN lasta)
50 : {
51 3431482 : pari_sp av = avma;
52 : long j;
53 3431482 : ulong nn = itou_or_0(n);
54 3431482 : if (nn && nn <= bg->rootbnd) bg->an[nn] = itos(a);
55 :
56 3431482 : if (signe(a))
57 : {
58 867121 : fun(E, n, a);
59 867121 : j = 1;
60 : }
61 : else
62 2564361 : j = i;
63 6857939 : for(; j <= i; j++)
64 : {
65 5484935 : ulong p = bg->p[j];
66 5484935 : GEN nexta, pn = mului(p, n);
67 5484935 : if (cmpii(pn, bg->bnd) > 0) return;
68 3426457 : nexta = mulis(a, bg->an[p]);
69 3426457 : if (i == j && umodiu(bg->N, p)) nexta = subii(nexta, mului(p, lasta));
70 3426457 : gen_BG_add(E, fun, bg, pn, j, nexta, a);
71 3426457 : set_avma(av);
72 : }
73 : }
74 :
75 : static void
76 42 : gen_BG_init(struct bg_data *bg, GEN E, GEN N, GEN bnd)
77 : {
78 42 : bg->E = E;
79 42 : bg->N = N;
80 42 : bg->bnd = bnd;
81 42 : bg->rootbnd = itou(sqrtint(bnd));
82 42 : bg->p = primes_upto_zv(bg->rootbnd);
83 42 : bg->an = ellanQ_zv(E, bg->rootbnd);
84 42 : }
85 :
86 : static void
87 42 : gen_BG_rec(void *E, bg_fun *fun, struct bg_data *bg)
88 : {
89 42 : long i, j, lp = lg(bg->p)-1;
90 42 : GEN bndov2 = shifti(bg->bnd, -1);
91 42 : pari_sp av = avma, av2;
92 : GEN p;
93 : forprime_t S;
94 42 : (void)forprime_init(&S, utoipos(bg->p[lp]+1), bg->bnd);
95 42 : av2 = avma;
96 42 : if (DEBUGLEVEL)
97 0 : err_printf("1st stage, using recursion for p <= %ld\n", bg->p[lp]);
98 5067 : for (i = 1; i <= lp; i++)
99 : {
100 5025 : ulong pp = bg->p[i];
101 5025 : long ap = bg->an[pp];
102 5025 : gen_BG_add(E, fun, bg, utoipos(pp), i, stoi(ap), gen_1);
103 5025 : set_avma(av2);
104 : }
105 42 : if (DEBUGLEVEL) err_printf("2nd stage, looping for p <= %Ps\n", bndov2);
106 2393586 : while ( (p = forprime_next(&S)) )
107 : {
108 : long jmax;
109 2393586 : GEN ap = ellap(bg->E, p);
110 2393586 : pari_sp av3 = avma;
111 2393586 : if (!signe(ap)) continue;
112 :
113 1196490 : jmax = itou( divii(bg->bnd, p) ); /* 2 <= jmax <= el->rootbound */
114 1196490 : fun(E, p, ap);
115 20342791 : for (j = 2; j <= jmax; j++)
116 : {
117 19146301 : long aj = bg->an[j];
118 : GEN a, n;
119 19146301 : if (!aj) continue;
120 2942771 : a = mulis(ap, aj);
121 2942771 : n = muliu(p, j);
122 2942771 : fun(E, n, a);
123 2942771 : set_avma(av3);
124 : }
125 1196490 : set_avma(av2);
126 1196490 : if (abscmpii(p, bndov2) >= 0) break;
127 : }
128 42 : if (DEBUGLEVEL) err_printf("3nd stage, looping for p <= %Ps\n", bg->bnd);
129 2164631 : while ( (p = forprime_next(&S)) )
130 : {
131 2164589 : GEN ap = ellap(bg->E, p);
132 2164589 : if (!signe(ap)) continue;
133 1082114 : fun(E, p, ap);
134 1082114 : set_avma(av2);
135 : }
136 42 : set_avma(av);
137 42 : }
138 :
139 : /******************************************************************
140 : *
141 : * L functions of elliptic curves
142 : * Pascal Molin (molin.maths@gmail.com) 2014
143 : *
144 : ******************************************************************/
145 :
146 : struct lcritical
147 : {
148 : GEN h; /* real */
149 : long cprec; /* computation prec */
150 : long L; /* number of points */
151 : GEN K; /* length of series */
152 : long real;
153 : };
154 :
155 : static double
156 245 : logboundG0(long e, double aY)
157 : {
158 : double cla, loggam;
159 245 : cla = 1 + 1/sqrt(aY);
160 245 : if (e) cla = ( cla + 1/(2*aY) ) / (2*sqrt(aY));
161 245 : loggam = (e) ? M_LN2-aY : -aY + log( log( 1+1/aY) );
162 245 : return log(cla) + loggam;
163 : }
164 :
165 : static void
166 245 : param_points(GEN N, double Y, double tmax, long bprec, long *cprec, long *L,
167 : GEN *K, double *h)
168 : {
169 : double D, a, aY, X, logM;
170 245 : long d = 2, w = 1;
171 245 : tmax *= d;
172 245 : D = bprec * M_LN2 + M_PI/4*tmax + 2;
173 245 : *cprec = nbits2prec(ceil(D / M_LN2) + 5);
174 245 : a = 2 * M_PI / sqrt(gtodouble(N));
175 245 : aY = a * cos(M_PI/2*Y);
176 245 : logM = 2*M_LN2 + logboundG0(w+1, aY) + tmax * Y * M_PI/2;
177 245 : *h = ( 2 * M_PI * M_PI / 2 * Y ) / ( D + logM );
178 245 : X = log( D / a);
179 245 : *L = ceil( X / *h);
180 245 : *K = ceil_safe(dbltor( D / a ));
181 245 : }
182 :
183 : static GEN
184 245 : vecF2_lk(GEN E, GEN K, GEN rbnd, GEN Q, GEN sleh, long prec)
185 : {
186 : pari_sp av;
187 245 : long l, L = lg(K)-1;
188 245 : GEN a = ellanQ_zv(E, itos(gel(K,1)));
189 245 : GEN S = cgetg(L+1, t_VEC);
190 :
191 13776 : for (l = 1; l <= L; l++) gel(S,l) = cgetr(prec);
192 245 : av = avma;
193 13776 : for (l = 1; l <= L; l++)
194 : {
195 : GEN e1, Sl, z, zB;
196 13531 : long aB, b, A, B, Kl = itou(gel(K,l));
197 : pari_sp av2;
198 : /* FIXME: could reduce prec here (useful for large prec) */
199 13531 : e1 = gel(Q, l);
200 13531 : Sl = real_0(prec);;
201 : /* baby-step giant step */
202 13531 : B = A = rbnd[l];
203 13531 : z = powersr(e1, B); zB = gel(z, B+1);
204 13531 : av2 = avma;
205 314174 : for (aB = A*B; aB >= 0; aB -= B)
206 : {
207 300643 : GEN s = real_0(prec); /* could change also prec here */
208 38032519 : for (b = B; b > 0; --b)
209 : {
210 37731876 : long k = aB+b;
211 37731876 : if (k <= Kl && a[k]) s = addrr(s, mulsr(a[k], gel(z, b+1)));
212 37731876 : if (gc_needed(av2, 1)) (void)gc_all(av2, 2, &s, &Sl);
213 : }
214 300643 : Sl = addrr(mulrr(Sl, zB), s);
215 : }
216 13531 : affrr(mulrr(Sl, gel(sleh,l)), gel(S, l)); /* to avoid copying all S */
217 13531 : set_avma(av);
218 : }
219 245 : return S;
220 : }
221 :
222 : /* Return C, C[i][j] = Q[j]^i, i = 1..nb */
223 : static void
224 0 : baby_init(struct baby_giant *bb, GEN Q, GEN bnd, GEN rbnd, long prec)
225 : {
226 0 : long i, j, l = lg(Q);
227 : GEN R, C;
228 0 : C = cgetg(l,t_VEC);
229 0 : for (i = 1; i < l; ++i) gel(C, i) = powersr(gel(Q, i), rbnd[i]);
230 0 : R = cgetg(l,t_VEC);
231 0 : for (i = 1; i < l; ++i)
232 : {
233 0 : gel(R, i) = cgetg(rbnd[i]+1, t_VEC);
234 0 : gmael(R, i, 1) = rtor(gmael(C, i, 2), prec);
235 0 : for (j = 2; j <= rbnd[i]; j++) gmael(R, i, j) = stor(0, prec);
236 : }
237 0 : bb->baby = C; bb->giant = R;
238 0 : bb->bnd = bnd; bb->rbnd = rbnd;
239 0 : }
240 :
241 : static long
242 245 : baby_size(GEN rbnd, long Ks, long prec)
243 : {
244 245 : long i, s, m, l = lg(rbnd);
245 13776 : for (s = 0, i = 1; i < l; ++i)
246 13531 : s += rbnd[i];
247 245 : m = 2*s*prec + 3*l + s;
248 245 : if (DEBUGLEVEL)
249 0 : err_printf("ellL1: BG_add: %ld words, ellan: %ld words\n", m, Ks);
250 245 : return m;
251 : }
252 :
253 : static void
254 0 : ellL1_add(void *E, GEN n, GEN a)
255 : {
256 0 : pari_sp av = avma;
257 0 : struct baby_giant *bb = (struct baby_giant*) E;
258 0 : long j, l = lg(bb->giant);
259 0 : for (j = 1; j < l; j++)
260 0 : if (cmpii(n, gel(bb->bnd,j)) <= 0)
261 : {
262 0 : ulong r, q = uabsdiviu_rem(n, bb->rbnd[j], &r);
263 0 : GEN giant = gel(bb->giant, j), baby = gel(bb->baby, j);
264 0 : affrr(addrr(gel(giant, q+1), mulri(gel(baby, r+1), a)), gel(giant, q+1));
265 0 : set_avma(av);
266 0 : } else break;
267 0 : }
268 :
269 : static GEN
270 0 : vecF2_lk_bsgs(GEN E, GEN bnd, GEN rbnd, GEN Q, GEN sleh, GEN N, long prec)
271 : {
272 : struct bg_data bg;
273 : struct baby_giant bb;
274 0 : long k, L = lg(bnd)-1;
275 : GEN S;
276 0 : baby_init(&bb, Q, bnd, rbnd, prec);
277 0 : gen_BG_init(&bg, E, N, gel(bnd,1));
278 0 : gen_BG_rec((void*) &bb, ellL1_add, &bg);
279 0 : S = cgetg(L+1, t_VEC);
280 0 : for (k = 1; k <= L; ++k)
281 : {
282 0 : pari_sp av = avma;
283 0 : long j, g = rbnd[k];
284 0 : GEN giant = gmael(bb.baby, k, g+1), Sl = gmael(bb.giant, k, g);
285 0 : for (j = g-1; j >=1; j--) Sl = addrr(mulrr(Sl, giant), gmael(bb.giant,k,j));
286 0 : gel(S, k) = gc_leaf(av, mulrr(gel(sleh,k), Sl));
287 : }
288 0 : return S;
289 : }
290 :
291 : static long
292 13531 : _sqrt(GEN x) { pari_sp av = avma; return gc_long(av, itou(sqrtint(x))); }
293 :
294 : static GEN
295 245 : vecF(struct lcritical *C, GEN E)
296 : {
297 245 : pari_sp av = avma;
298 245 : long prec = C->cprec, Ks = itos_or_0(C->K), L = C->L, l;
299 245 : GEN N = ellQ_get_N(E), PiN;
300 245 : GEN e = mpexp(C->h), elh = powersr(e, L-1), Q, bnd, rbnd, vec;
301 :
302 245 : PiN = divrr(Pi2n(1,prec), sqrtr_abs(itor(N, prec)));
303 245 : setsigne(PiN, -1); /* - 2Pi/sqrt(N) */
304 245 : bnd = gpowers0(invr(e), L-1, C->K); /* bnd[i] = K exp(-(i-1)h) */
305 245 : rbnd = cgetg(L+1, t_VECSMALL);
306 245 : Q = cgetg(L+1, t_VEC);
307 13776 : for (l = 1; l <= L; l++)
308 : {
309 13531 : gel(bnd,l) = ceil_safe(gel(bnd,l));
310 13531 : rbnd[l] = _sqrt(gel(bnd,l)) + 1;
311 13531 : gel(Q, l) = mpexp(mulrr(PiN, gel(elh, l)));
312 : }
313 245 : if (Ks && baby_size(rbnd, Ks, prec) > (Ks>>1))
314 245 : vec = vecF2_lk(E, bnd, rbnd, Q, elh, prec);
315 : else
316 0 : vec = vecF2_lk_bsgs(E, bnd, rbnd, Q, elh, N, prec);
317 245 : return gc_upto(av, vec);
318 : }
319 :
320 : /* Lambda function by Fourier inversion. vec is a grid, t a scalar or t_SER */
321 : static GEN
322 273 : glambda(GEN t, GEN vec, GEN h, long real, long prec)
323 : {
324 273 : GEN z, r, e = gexp(gmul(mkcomplex(gen_0,h), t), prec);
325 273 : long n = lg(vec)-1, i;
326 :
327 273 : r = real == 1? gmul2n(real_i(gel(vec, 1)), -1): gen_0;
328 273 : z = real == 1? e: gmul(powIs(3), e);
329 : /* FIXME: summing backward may be more stable */
330 16268 : for (i = 2; i <= n; i++)
331 : {
332 15995 : r = gadd(r, real_i(gmul(gel(vec,i), z)));
333 15995 : if (i < n) z = gmul(z, e);
334 : }
335 273 : return gmul(mulsr(4, h), r);
336 : }
337 :
338 : static GEN
339 245 : Lpoints(struct lcritical *C, GEN e, double tmax, long bprec)
340 : {
341 245 : double h = 0, Y = .97;
342 245 : GEN N = ellQ_get_N(e);
343 245 : param_points(N, Y, tmax, bprec, &C->cprec, &C->L, &C->K, &h);
344 245 : C->real = ellrootno_global(e);
345 245 : C->h = rtor(dbltor(h), C->cprec);
346 245 : return vecF(C, e);
347 : }
348 :
349 : static GEN
350 273 : Llambda(GEN vec, struct lcritical *C, GEN t, long prec)
351 : {
352 273 : GEN lambda = glambda(gprec_w(t, C->cprec), vec, C->h, C->real, C->cprec);
353 273 : return gprec_w(lambda, prec);
354 : }
355 :
356 : /* 2*(2*Pi)^(-s)*gamma(s)*N^(s/2); */
357 : static GEN
358 273 : ellgammafactor(GEN N, GEN s, long prec)
359 : {
360 273 : GEN c = gpow(divrr(gsqrt(N,prec), Pi2n(1,prec)), s, prec);
361 273 : return gmul(gmul2n(c,1), ggamma(s, prec));
362 : }
363 :
364 : static GEN
365 273 : ellL1_eval(GEN e, GEN vec, struct lcritical *C, GEN t, long prec)
366 : {
367 273 : GEN g = ellgammafactor(ellQ_get_N(e), gaddgs(gmul(gen_I(),t), 1), prec);
368 273 : return gdiv(Llambda(vec, C, t, prec), g);
369 : }
370 :
371 : static GEN
372 273 : ellL1_der(GEN e, GEN vec, struct lcritical *C, GEN t, long der, long prec)
373 : {
374 273 : GEN r = polcoef_i(ellL1_eval(e, vec, C, t, prec), der, 0);
375 273 : r = gmul(r,powIs(C->real == 1 ? -der: 1-der));
376 273 : return gmul(real_i(r), mpfact(der));
377 : }
378 :
379 : GEN
380 231 : ellL1(GEN E, long r, long bitprec)
381 : {
382 231 : pari_sp av = avma;
383 : struct lcritical C;
384 231 : long prec = nbits2prec(bitprec);
385 : GEN e, vec, t;
386 231 : if (r < 0)
387 7 : pari_err_DOMAIN("ellL1", "derivative order", "<", gen_0, stoi(r));
388 224 : e = ellanal_globalred(E, NULL);
389 224 : if (r == 0 && ellrootno_global(e) < 0) { set_avma(av); return gen_0; }
390 210 : vec = Lpoints(&C, e, 0., bitprec);
391 210 : t = r ? scalarser(gen_1, 0, r): zeroser(0, 0);
392 210 : setvalser(t, 1);
393 210 : return gc_upto(av, ellL1_der(e, vec, &C, t, r, prec));
394 : }
395 :
396 : GEN
397 35 : ellanalyticrank(GEN E, GEN eps, long bitprec)
398 : {
399 35 : pari_sp av = avma, av2;
400 35 : long prec = nbits2prec(bitprec);
401 : struct lcritical C;
402 : pari_timer ti;
403 : GEN e, vec;
404 : long rk;
405 35 : if (DEBUGLEVEL) timer_start(&ti);
406 35 : if (!eps)
407 35 : eps = real2n(-bitprec/2+1, DEFAULTPREC);
408 : else
409 0 : if (typ(eps) != t_REAL) {
410 0 : eps = gtofp(eps, DEFAULTPREC);
411 0 : if (typ(eps) != t_REAL) pari_err_TYPE("ellanalyticrank", eps);
412 : }
413 35 : e = ellanal_globalred(E, NULL);
414 35 : vec = Lpoints(&C, e, 0., bitprec);
415 35 : if (DEBUGLEVEL) timer_printf(&ti, "init L");
416 35 : av2 = avma;
417 35 : for (rk = C.real>0 ? 0: 1; ; rk += 2)
418 28 : {
419 : GEN Lrk;
420 63 : GEN t = rk ? scalarser(gen_1, 0, rk): zeroser(0, 0);
421 63 : setvalser(t, 1);
422 63 : Lrk = ellL1_der(e, vec, &C, t, rk, prec);
423 63 : if (DEBUGLEVEL) timer_printf(&ti, "L^(%ld)=%Ps", rk, Lrk);
424 63 : if (abscmprr(Lrk, eps) > 0)
425 35 : return gc_GEN(av, mkvec2(stoi(rk), Lrk));
426 28 : set_avma(av2);
427 : }
428 : }
429 :
430 : /* Heegner point computation
431 :
432 : This section is a C version by Bill Allombert of a GP script by
433 : Christophe Delaunay which was based on a GP script by John Cremona.
434 : Reference: Henri Cohen's book GTM 239.
435 : */
436 :
437 : static void
438 0 : heegner_L1_bg(void*E, GEN n, GEN a)
439 : {
440 0 : struct baby_giant *bb = (struct baby_giant*) E;
441 0 : long j, l = lg(bb->giant);
442 0 : for (j = 1; j < l; j++)
443 0 : if (cmpii(n, gel(bb->bnd,j)) <= 0)
444 : {
445 0 : ulong r, q = uabsdiviu_rem(n, bb->rbnd[j], &r);
446 0 : GEN giant = gel(bb->giant, j), baby = gel(bb->baby, j);
447 0 : affgc(gadd(gel(giant, q+1), gdiv(gmul(gel(baby, r+1), a), n)),
448 0 : gel(giant, q+1));
449 : }
450 0 : }
451 :
452 : static void
453 6088496 : heegner_L1(void*E, GEN n, GEN a)
454 : {
455 6088496 : struct baby_giant *bb = (struct baby_giant*) E;
456 6088496 : long j, l = lg(bb->giant);
457 36228875 : for (j = 1; j < l; j++)
458 30140379 : if (cmpii(n, gel(bb->bnd,j)) <= 0)
459 : {
460 25949953 : ulong r, q = uabsdiviu_rem(n, bb->rbnd[j], &r);
461 25949953 : GEN giant = gel(bb->giant, j), baby = gel(bb->baby, j);
462 25949953 : GEN ex = mulreal(gel(baby, r+1), gel(giant, q+1));
463 25949953 : affrr(addrr(gel(bb->sum, j), divri(mulri(ex, a), n)),
464 25949953 : gel(bb->sum, j));
465 : }
466 6088496 : }
467 : /* export ? */
468 : static GEN
469 0 : ctoc(GEN x, long prec)
470 0 : { GEN y = cgetc(prec); affgc(x, y); return y; }
471 :
472 : /* [powers(x[i], n[i]), i=1..#x] */
473 : static GEN
474 42 : RgV_zv_powers(GEN x, GEN n)
475 189 : { pari_APPLY_same(gpowers(gel(x,i), n[i])); }
476 :
477 : /* Return C, C[i][j] = Q[j]^i, i = 1..nb */
478 : static void
479 0 : baby_init2(struct baby_giant *bb, GEN Q, GEN bnd, GEN rbnd, long prec)
480 : {
481 0 : long i, j, l = lg(Q);
482 0 : GEN C = RgV_zv_powers(Q, rbnd), R = cgetg(l,t_VEC);
483 0 : for (i = 1; i < l; ++i)
484 : {
485 0 : gel(R, i) = cgetg(rbnd[i]+1, t_VEC);
486 0 : gmael(R, i, 1) = ctoc(gmael(C, i, 2), prec);
487 0 : for (j = 2; j <= rbnd[i]; j++) gmael(R, i, j) = ctoc(gen_0, prec);
488 : }
489 0 : bb->baby = C; bb->giant = R;
490 0 : bb->bnd = bnd; bb->rbnd = rbnd;
491 0 : }
492 :
493 : /* Return C, C[i][j] = Q[j]^i, i = 1..nb */
494 : static void
495 42 : baby_init3(struct baby_giant *bb, GEN Q, GEN bnd, GEN rbnd, long prec)
496 : {
497 42 : long i, l = lg(Q);
498 42 : GEN S, C = RgV_zv_powers(Q, rbnd), R = cgetg(l,t_VEC);
499 189 : for (i = 1; i < l; ++i)
500 147 : gel(R, i) = gpowers(gmael(C, i, 1+rbnd[i]), rbnd[i]);
501 42 : S = cgetg(l,t_VEC);
502 189 : for (i = 1; i < l; ++i) gel(S, i) = rtor(real_i(gmael(C, i, 2)), prec);
503 42 : bb->baby = C; bb->giant = R; bb->sum = S;
504 42 : bb->bnd = bnd; bb->rbnd = rbnd;
505 42 : }
506 :
507 : /* ymin a t_REAL */
508 : static GEN
509 42 : heegner_psi(GEN E, GEN N, GEN points, long bitprec)
510 : {
511 42 : pari_sp av = avma, av2;
512 : struct baby_giant bb;
513 : struct bg_data bg;
514 42 : long l, k, L = lg(points)-1, prec = nbits2prec(bitprec)+EXTRAPREC64;
515 42 : GEN Q, pi2 = Pi2n(1, prec), bnd, rbnd, bndmax;
516 42 : GEN B = divrr(mulur(bitprec,mplog2(DEFAULTPREC)), pi2);
517 :
518 42 : rbnd = cgetg(L+1, t_VECSMALL); av2 = avma;
519 42 : bnd = cgetg(L+1, t_VEC);
520 42 : Q = cgetg(L+1, t_VEC);
521 189 : for (l = 1; l <= L; ++l)
522 : {
523 147 : gel(bnd,l) = ceil_safe(divrr(B,imag_i(gel(points, l))));
524 147 : rbnd[l] = itou(sqrtint(gel(bnd,l)))+1;
525 147 : gel(Q, l) = expIxy(pi2, gel(points, l), prec);
526 : }
527 42 : (void)gc_all(av2, 2, &bnd, &Q);
528 42 : bndmax = gel(bnd,vecindexmax(bnd));
529 42 : gen_BG_init(&bg, E, N, bndmax);
530 42 : if (bitprec >= 1900)
531 : {
532 0 : GEN S = cgetg(L+1, t_VEC);
533 0 : baby_init2(&bb, Q, bnd, rbnd, prec);
534 0 : gen_BG_rec((void*)&bb, heegner_L1_bg, &bg);
535 0 : for (k = 1; k <= L; ++k)
536 : {
537 0 : pari_sp av2 = avma;
538 0 : long j, g = rbnd[k];
539 0 : GEN giant = gmael(bb.baby, k, g+1), Sl = real_0(prec);
540 0 : for (j = g; j >=1; j--) Sl = gadd(gmul(Sl, giant), gmael(bb.giant,k,j));
541 0 : gel(S, k) = gc_upto(av2, real_i(Sl));
542 : }
543 0 : return gc_upto(av, S);
544 : }
545 : else
546 : {
547 42 : baby_init3(&bb, Q, bnd, rbnd, prec);
548 42 : gen_BG_rec((void*)&bb, heegner_L1, &bg);
549 42 : return gc_GEN(av, bb.sum);
550 : }
551 : }
552 :
553 : /*Returns lambda_bad list for one prime p, nv = localred(E, p) */
554 : static GEN
555 91 : lambda1(GEN E, GEN nv, GEN p, long prec)
556 : {
557 : GEN res, lp;
558 91 : long kod = itos(gel(nv, 2));
559 91 : if (kod == 2 || kod == -2) return NULL;
560 91 : lp = glog(p, prec);
561 91 : if (kod > 4)
562 : { /* I_m */
563 14 : long n = Z_pval(ell_get_disc(E), p);
564 14 : long j, m = kod - 4, nl = 1 + (m >> 1L);
565 14 : res = cgetg(nl, t_VEC);
566 35 : for (j = 1; j < nl; j++) /* log(p) (j/n) (j - n) */
567 21 : gel(res, j) = divru(mulri(lp, mulss(j, j-n)), n);
568 : }
569 77 : else if (kod < -4) /* I^*_m */
570 14 : res = mkvec2(negr(lp), shiftr(mulrs(lp, kod), -2));
571 : else
572 : {
573 63 : const long lam[] = {8,9,0,6,0,0,0,3,4};
574 63 : long m = -lam[kod+4];
575 63 : res = mkvec(divru(mulrs(lp, m), 6));
576 : }
577 91 : return res;
578 : }
579 :
580 : /* lambda[1] = gen_0, other components are t_REAL */
581 : static GEN
582 42 : lambdalist(GEN E, long prec)
583 : {
584 42 : pari_sp ltop = avma;
585 42 : GEN glob = ellglobalred(E), P = gmael(glob,4,1), L = gel(glob,5);
586 42 : GEN res, v, D = ell_get_disc(E);
587 42 : long i, j, k, l, m, n, lv = lg(P), lr = 1;
588 42 : v = cgetg(lv, t_VEC);
589 147 : for (j = 1, i = 1 ; j < lv; j++)
590 : {
591 105 : GEN p = gel(P, j);
592 105 : if (dvdii(D, sqri(p)))
593 : {
594 91 : GEN la = lambda1(E, gel(L,j), p, prec);
595 91 : if (la) { gel(v, i++) = la; lr *= lg(la); }
596 : }
597 : }
598 42 : lv = i;
599 42 : res = cgetg(lr+1, t_VEC); gel(res, 1) = gen_0;
600 133 : for (j = n = 1; j < lv; j++)
601 : {
602 91 : GEN w = gel(v, j); /* vector of t_REAL */
603 91 : long lw = lg(w);
604 : /* k = 1 */
605 203 : for (l = 1, m = n; l < lw; l++, m+=n)
606 112 : gel(res, 1 + m) = gel(w, l);
607 217 : for (k = 2; k <= n; k++)
608 : {
609 126 : GEN t = gel(res, k);
610 252 : for (l = 1, m = n; l < lw; l++, m+=n)
611 126 : gel(res, k + m) = addrr(t, gel(w, l));
612 : }
613 91 : n = m;
614 : }
615 42 : return gc_GEN(ltop, res);
616 : }
617 :
618 : /* P a t_INT or t_FRAC, return its logarithmic height */
619 : static GEN
620 98 : heightQ(GEN P, long prec)
621 : {
622 : long s;
623 98 : if (typ(P) == t_FRAC)
624 : {
625 56 : GEN a = gel(P,1), b = gel(P,2);
626 56 : P = abscmpii(a,b) > 0 ? a: b;
627 : }
628 98 : s = signe(P);
629 98 : if (!s) return real_0(prec);
630 84 : if (s < 0) P = negi(P);
631 84 : return glog(P, prec);
632 : }
633 :
634 : /* t a t_INT or t_FRAC, returns max(1, log |t|), returns a t_REAL */
635 : static GEN
636 147 : logplusQ(GEN t, long prec)
637 : {
638 147 : if (typ(t) == t_INT)
639 : {
640 42 : if (!signe(t)) return real_1(prec);
641 28 : t = absi_shallow(t);
642 : }
643 : else
644 : {
645 105 : GEN a = gel(t,1), b = gel(t,2);
646 105 : if (abscmpii(a, b) < 0) return real_1(prec);
647 56 : if (signe(a) < 0) t = gneg(t);
648 : }
649 84 : return glog(t, prec);
650 : }
651 :
652 : /* See GTM239, p532, Th 8.1.18
653 : * Return M such that h_naive <= M */
654 : GEN
655 98 : hnaive_max(GEN ell, GEN ht)
656 : {
657 98 : const long prec = LOWDEFAULTPREC; /* minimal accuracy */
658 98 : GEN b2 = ell_get_b2(ell), j = ell_get_j(ell);
659 98 : GEN logd = glog(absi_shallow(ell_get_disc(ell)), prec);
660 98 : GEN logj = logplusQ(j, prec);
661 98 : GEN hj = heightQ(j, prec);
662 49 : GEN logb2p = signe(b2)? addrr(logplusQ(gdivgu(b2, 12),prec), mplog2(prec))
663 98 : : real_1(prec);
664 98 : GEN mu = addrr(divru(addrr(logd, logj),6), logb2p);
665 98 : return addrs(addrr(addrr(ht, divru(hj,12)), mu), 2);
666 : }
667 :
668 : static GEN
669 147 : qfb_root(GEN Q, GEN vDi)
670 : {
671 147 : GEN a2 = shifti(gel(Q, 1),1), b = gel(Q, 2);
672 147 : return mkcomplex(gdiv(negi(b),a2),divri(vDi,a2));
673 : }
674 :
675 : static GEN
676 24668 : qimag2(GEN Q)
677 : {
678 24668 : pari_sp av = avma;
679 24668 : GEN z = gdiv(negi(qfb_disc(Q)), shifti(sqri(gel(Q, 1)),2));
680 24668 : return gc_upto(av, z);
681 : }
682 :
683 : /***************************************************/
684 : /*Routines for increasing the imaginary parts using*/
685 : /*Atkin-Lehner operators */
686 : /***************************************************/
687 :
688 : static GEN
689 24668 : qfb_mult(GEN Q, GEN a, GEN b, GEN c, GEN d)
690 : {
691 24668 : GEN A = gel(Q, 1), B = gel(Q, 2), C = gel(Q, 3), D = qfb_disc(Q);
692 24668 : GEN a2 = sqri(a), b2 = sqri(b), c2 = sqri(c), d2 = sqri(d);
693 24668 : GEN ad = mulii(d, a), bc = mulii(b, c), e = subii(ad, bc);
694 24668 : GEN W1 = addii(addii(mulii(a2, A), mulii(mulii(c, a), B)), mulii(c2, C));
695 24668 : GEN W3 = addii(addii(mulii(b2, A), mulii(mulii(d, b), B)), mulii(d2, C));
696 24668 : GEN W2 = addii(addii(mulii(mulii(shifti(b,1), a), A),
697 : mulii(addii(ad, bc), B)),
698 : mulii(mulii(shifti(d,1), c), C));
699 24668 : if (!equali1(e)) {
700 22190 : W1 = diviiexact(W1,e);
701 22190 : W2 = diviiexact(W2,e);
702 22190 : W3 = diviiexact(W3,e);
703 : }
704 24668 : return mkqfb(W1, W2, W3, D);
705 : }
706 :
707 : #ifdef DEBUG
708 : static void
709 : best_point_old(GEN Q, GEN NQ, GEN f, GEN *u, GEN *v)
710 : {
711 : long n, k;
712 : GEN U, c, d, A = gel(f,1), B = gel(f,2), C = gel(f,3), D = qfb_disc(f);
713 : GEN q = mkqfb(mulii(NQ, C), negi(B), diviiexact(A, NQ), D);
714 : redimagsl2(q, &U);
715 : *u = c = gcoeff(U, 1, 1);
716 : *v = d = gcoeff(U, 2, 1);
717 : if (equali1(gcdii(mulii(*u, NQ), mulii(*v, Q)))) return;
718 : for (n = 1;; n++)
719 : {
720 : for (k = -n; k <= n; k++)
721 : {
722 : *u = addis(c, k); *v = addiu(d, n);
723 : if (equali1(gcdii(mulii(*u, NQ), mulii(*v, Q)))) return;
724 : *v = subiu(d, n);
725 : if (equali1(gcdii(mulii(*u, NQ), mulii(*v, Q)))) return;
726 : *u = addiu(c, n); *v = addis(d, k);
727 : if (equali1(gcdii(mulii(*u, NQ), mulii(*v, Q)))) return;
728 : *u = subiu(c, n);
729 : if (equali1(gcdii(mulii(*u, NQ), mulii(*v, Q)))) return;
730 : }
731 : }
732 : }
733 : /* q(x,y) = ax^2 + bxy + cy^2 */
734 : static GEN
735 : qfb_eval(GEN q, GEN x, GEN y)
736 : {
737 : GEN a = gel(q,1), b = gel(q,2), c = gel(q,3);
738 : GEN x2 = sqri(x), y2 = sqri(y), xy = mulii(x,y);
739 : return addii(addii(mulii(a, x2), mulii(b,xy)), mulii(c, y2));
740 : }
741 : #endif
742 :
743 : static long
744 6594 : nexti(long i) { return i>0 ? -i : 1-i; }
745 :
746 : /* q0 + i q1 + i^2 q2 */
747 : static GEN
748 12334 : qfmin_eval(GEN q0, GEN q1, GEN q2, long i)
749 12334 : { return addii(mulis(addii(mulis(q2, i), q1), i), q0); }
750 :
751 : /* assume a > 0, return gcd(a,b,c) */
752 : static ulong
753 16443 : gcduii(ulong a, GEN b, GEN c)
754 : {
755 16443 : a = ugcdiu(b, a);
756 16443 : return a == 1? 1: ugcdiu(c, a);
757 : }
758 :
759 : static void
760 24668 : best_point(GEN Q, GEN NQ, GEN f, GEN *pu, GEN *pv)
761 : {
762 24668 : GEN a = mulii(NQ, gel(f,3)), b = negi(gel(f,2)), c = diviiexact(gel(f,1), NQ);
763 24668 : GEN D = qfb_disc(f);
764 24668 : GEN U, qr = redimagsl2(mkqfb(a, b, c, D), &U);
765 24668 : GEN A = gel(qr,1), B = gel(qr,2), A2 = shifti(A,1), AA4 = sqri(A2);
766 : GEN V, best;
767 : long y;
768 :
769 24668 : D = absi_shallow(D);
770 : /* 4A qr(x,y) = (2A x + By)^2 + D y^2
771 : * Write x = x0(y) + i, where x0 is an integer minimum
772 : * (the smallest in case of tie) of x-> qr(x,y), for given y.
773 : * 4A qr(x,y) = ((2A x0 + By)^2 + Dy^2) + 4A i (2A x0 + By) + 4A^2 i^2
774 : * = q0(y) + q1(y) i + q2 i^2
775 : * Loop through (x,y), y>0 by (roughly) increasing values of qr(x,y) */
776 :
777 : /* We must test whether [X,Y]~ := U * [x,y]~ satisfy (X NQ, Y Q) = 1
778 : * This is equivalent to (X,Y) = 1 (note that (X,Y) = (x,y)), and
779 : * (X, Q) = (Y, NQ) = 1.
780 : * We have U * [x0+i, y]~ = U * [x0,y]~ + i U[,1] =: V0 + i U[,1] */
781 :
782 : /* try [1,0]~ = first minimum */
783 24668 : V = gel(U,1); /* U *[1,0]~ */
784 24668 : *pu = gel(V,1);
785 24668 : *pv = gel(V,2);
786 30947 : if (is_pm1(gcdii(*pu, Q)) && is_pm1(gcdii(*pv, NQ))) return;
787 :
788 : /* try [0,1]~ = second minimum */
789 11935 : V = gel(U,2); /* U *[0,1]~ */
790 11935 : *pu = gel(V,1);
791 11935 : *pv = gel(V,2);
792 11935 : if (is_pm1(gcdii(*pu, Q)) && is_pm1(gcdii(*pv, NQ))) return;
793 :
794 : /* (X,Y) = (1, \pm1) always works. Try to do better now */
795 5656 : best = subii(addii(a, c), absi_shallow(b));
796 5656 : *pu = gen_1;
797 5656 : *pv = signe(b) < 0? gen_1: gen_m1;
798 :
799 5656 : for (y = 1;; y++)
800 9135 : {
801 : GEN Dy2, r, By, x0, q0, q1, V0;
802 : long i;
803 14791 : if (y > 1)
804 : {
805 10703 : if (gcduii(y, gcoeff(U,1,1), Q) != 1) continue;
806 7308 : if (gcduii(y, gcoeff(U,2,1), NQ) != 1) continue;
807 : }
808 11403 : Dy2 = mulii(D, sqru(y));
809 11403 : if (cmpii(Dy2, best) >= 0) break; /* we won't improve. STOP */
810 5747 : By = muliu(B,y), x0 = truedvmdii(negi(By), A2, &r);
811 5747 : if (cmpii(r, A) >= 0) { x0 = subiu(x0,1); r = subii(r, A2); }
812 : /* (2A x + By)^2 + Dy^2, minimal at x = x0. Assume A2 > 0 */
813 : /* r = 2A x0 + By */
814 5747 : q0 = addii(sqri(r), Dy2); /* minimal value for this y, at x0 */
815 5747 : if (cmpii(q0, best) >= 0) continue; /* we won't improve for this y */
816 5740 : q1 = shifti(mulii(A2, r), 1);
817 :
818 5740 : V0 = ZM_ZC_mul(U, mkcol2(x0, utoipos(y)));
819 12334 : for (i = 0;; i = nexti(i))
820 6594 : {
821 12334 : pari_sp av2 = avma;
822 12334 : GEN x, N = qfmin_eval(q0, q1, AA4, i);
823 12334 : if (cmpii(N, best) >= 0) break;
824 12299 : x = addis(x0, i);
825 12299 : if (ugcdiu(x, y) == 1)
826 : {
827 : GEN u, v;
828 12257 : V = ZC_add(V0, ZC_z_mul(gel(U,1), i)); /* [X, Y] */
829 12257 : u = gel(V,1);
830 12257 : v = gel(V,2);
831 12257 : if (is_pm1(gcdii(u, Q)) && is_pm1(gcdii(v, NQ)))
832 : {
833 5705 : *pu = u;
834 5705 : *pv = v;
835 5705 : best = N; break;
836 : }
837 : }
838 6594 : set_avma(av2);
839 : }
840 : }
841 : #ifdef DEBUG
842 : {
843 : GEN oldu, oldv, F = mkqfb(a, b, c, qfb_disc(f));
844 : best_point_old(Q, NQ, f, &oldu, &oldv);
845 : if (!equalii(oldu, *pu) || !equalii(oldv, *pv))
846 : {
847 : if (!equali1(gcdii(mulii(*pu, NQ), mulii(*pv, Q))))
848 : pari_err_BUG("best_point (gcd)");
849 : if (cmpii(qfb_eval(F, *pu,*pv), qfb_eval(F, oldu, oldv)) > 0)
850 : {
851 : pari_warn(warner, "%Ps,%Ps,%Ps, %Ps > %Ps",
852 : Q,NQ,f, mkvec2(*pu,*pv), mkvec2(oldu,oldv));
853 : pari_err_BUG("best_point (too large)");
854 : }
855 : }
856 : }
857 : #endif
858 : }
859 :
860 : static GEN
861 24668 : best_lift(GEN Q, GEN NQ, GEN f)
862 : {
863 : GEN a, b, c, d, dQ, cNQ;
864 24668 : best_point(Q, NQ, f, &c, &d);
865 24668 : dQ = mulii(d, Q); cNQ = mulii(NQ, c);
866 24668 : (void)bezout(dQ, cNQ, &a, &b);
867 24668 : return qfb_mult(f, dQ, b, mulii(negi(Q),cNQ), mulii(a,Q));
868 : }
869 :
870 : static GEN
871 2478 : lift_points(GEN listQ, GEN f, GEN *pt, GEN *pQ)
872 : {
873 2478 : pari_sp av = avma;
874 2478 : GEN yf = gen_0, tf = NULL, Qf = NULL;
875 2478 : long k, l = lg(listQ);
876 27146 : for (k = 1; k < l; ++k)
877 : {
878 24668 : GEN c = gel(listQ, k), Q = gel(c,1), NQ = gel(c,2);
879 24668 : GEN t = best_lift(Q, NQ, f), y = qimag2(t);
880 24668 : if (gcmp(y, yf) > 0) { yf = y; Qf = Q; tf = t; }
881 : }
882 2478 : *pt = tf; *pQ = Qf; return gc_all(av, 3, &yf, pt, pQ);
883 : }
884 :
885 : /***************************/
886 : /* Twists */
887 : /***************************/
888 :
889 : static GEN
890 56 : ltwist1(GEN E, GEN d, long bitprec)
891 : {
892 56 : pari_sp av = avma;
893 56 : GEN Ed = elltwist(E, d), z = ellL1(Ed, 0, bitprec);
894 56 : obj_free(Ed); return gc_leaf(av, z);
895 : }
896 :
897 : /* omega(gcd(D, N)), given faN = factor(N) */
898 : static long
899 168 : omega_N_D(GEN faN, ulong D)
900 : {
901 168 : GEN P = gel(faN, 1);
902 168 : long i, l = lg(P), w = 0;
903 658 : for (i = 1; i < l; i++)
904 490 : if (dvdui(D, gel(P,i))) w++;
905 168 : return w;
906 : }
907 :
908 : static long
909 112 : get_w(long D)
910 : {
911 112 : switch(D)
912 : {
913 0 : case -3: return 3; break;
914 0 : case -4: return 2; break;
915 112 : default: return 1;
916 : }
917 : }
918 :
919 : static GEN
920 56 : heegner_indexmultD(GEN faN, GEN a, long D, GEN sqrtD)
921 : {
922 56 : pari_sp av = avma;
923 : GEN c;
924 : long w;
925 56 : switch(D)
926 : {
927 0 : case -3: w = 9; break;
928 0 : case -4: w = 4; break;
929 56 : default: w = 1;
930 : }
931 56 : c = shifti(stoi(w), omega_N_D(faN,-D)); /* (w(D)/2)^2 * 2^omega(gcd(D,N)) */
932 56 : return gc_upto(av, mulri(mulrr(a, sqrtD), c));
933 : }
934 :
935 : static GEN
936 399 : nf_to_basis(GEN nf, GEN x)
937 : {
938 399 : x = nf_to_scalar_or_basis(nf, x);
939 399 : if (typ(x)!=t_COL)
940 301 : x = scalarcol(x, nf_get_degree(nf));
941 399 : return x;
942 : }
943 :
944 : static GEN
945 196 : etnf_to_basis(GEN et, GEN x)
946 : {
947 196 : long i, l = lg(et);
948 196 : GEN V = cgetg(l, t_VEC);
949 595 : for (i = 1; i < l; i++)
950 399 : gel(V,i) = nf_to_basis(gel(et,i), x);
951 196 : return shallowconcat1(V);
952 : }
953 :
954 : static long
955 42 : etnf_get_varn(GEN et)
956 : {
957 42 : return nf_get_varn(gel(et,1));
958 : }
959 :
960 : static GEN
961 98 : heegner_descent_try_point(GEN nfA, GEN z, GEN den, long prec)
962 : {
963 98 : pari_sp av = avma;
964 98 : GEN etal = gel(nfA,1), A = gel(nfA,2), cb = gel(nfA,3);
965 98 : GEN u2 = gsqr(gel(cb,1)), r = gel(cb,2), zz = gdiv(gsub(z,r), u2);
966 98 : GEN al = gel(nfA,4), th = gel(nfA, 5), M = gel(nfA,6);
967 98 : GEN zk = gel(etal, 2), T = gel(etal,3), den2 = sqri(den);
968 98 : long i, j, n = lg(th)-1, l = lg(al);
969 98 : GEN be = cgetg(n+1, t_COL);
970 :
971 154 : for (j = 1; j < l; j++)
972 : {
973 98 : GEN aj = gel(al, j), Aj = gel(A,j);
974 371 : for (i = 1; i <= n; i++)
975 : {
976 273 : GEN b = gsqrt(gmul(gsub(zz, gel(th,i)), gel(aj,i)), prec);
977 273 : gel(be,i) = gmul(b, den);
978 : }
979 336 : for (i = 0; i < (1<<(n-1)); i++)
980 : { /* n <= 3, by altering be[1] and be[2] we try all sign choices */
981 280 : GEN V, S = RgM_solve_realimag(M, be);
982 : long eps;
983 280 : gel(be,1+odd(i)) = gneg(gel(be,1+odd(i)));
984 280 : S = grndtoi(S, &eps); if (eps > -7) continue;
985 42 : V = QXQ_mul(Aj, QXQ_sqr(RgV_RgC_mul(zk, S), T), T);
986 42 : if (typ(V) != t_POL || degpol(V) != 1) continue;
987 42 : if (gequal0(gadd(gel(V,3), den2)))
988 : {
989 42 : GEN x = gdiv(gel(V,2), den2);
990 42 : return gc_upto(av, gadd(gmul(x, u2), r));
991 : }
992 : }
993 : }
994 56 : return gc_NULL(av);
995 : }
996 :
997 : /* lambdas[1] = gen_0, the others are t_REAL */
998 : static GEN
999 1785 : heegner_try_point(GEN E, GEN nfA, GEN lambdas, GEN ht, GEN z, long prec)
1000 : {
1001 1785 : long l = lg(lambdas), i, eps;
1002 1785 : GEN P = real_i(pointell(E, z, prec)), x = gel(P,1);
1003 1785 : GEN rh = subrr(ht, shiftr(ellheightoo(E, P, prec),1));
1004 26572 : for (i = 1; i < l; ++i)
1005 : {
1006 24829 : GEN logd = shiftr(i == 1? rh: subrr(rh, gel(lambdas, i)), -1);
1007 24829 : GEN approxd = gexp(logd, prec), d = grndtoi(approxd, &eps);
1008 24829 : if (signe(d) > 0 && eps < -10)
1009 : {
1010 : GEN X, ylist;
1011 98 : if (DEBUGLEVEL > 2)
1012 0 : err_printf("\nTrying lambda number %ld, logd=%Ps, approxd=%Ps\n",
1013 : i, logd, approxd);
1014 98 : X = heegner_descent_try_point(nfA, x, d, prec);
1015 98 : if (X)
1016 : {
1017 42 : ylist = ellordinate(E, X, prec);
1018 42 : if (lg(ylist) > 1)
1019 : {
1020 42 : GEN P = mkvec2(X, gel(ylist, 1));
1021 42 : GEN hp = ellheight(E,P,prec);
1022 42 : if (signe(hp) && cmprr(hp, shiftr(ht,1)) < 0
1023 42 : && cmprr(hp, shiftr(ht,-1)) > 0)
1024 42 : return P;
1025 0 : if (DEBUGLEVEL)
1026 0 : err_printf("found non-Heegner point %Ps\n", P);
1027 : }
1028 : }
1029 : }
1030 : }
1031 1743 : return NULL;
1032 : }
1033 :
1034 : static GEN
1035 42 : heegner_find_point(GEN e, GEN nfA, GEN om, GEN ht, GEN z1, long k, long prec)
1036 : {
1037 42 : GEN lambdas = lambdalist(e, prec);
1038 42 : pari_sp av = avma;
1039 : long m;
1040 42 : GEN Ore = gel(om, 1), Oim = gel(om, 2);
1041 42 : if (DEBUGLEVEL)
1042 0 : err_printf("%ld*%ld multipliers to test: ",k,lg(lambdas)-1);
1043 966 : for (m = 0; m < k; m++)
1044 : {
1045 966 : GEN P, z2 = divru(addrr(z1, mulsr(m, Ore)), k);
1046 966 : if (DEBUGLEVEL > 2)
1047 0 : err_printf("%ld ",m);
1048 966 : P = heegner_try_point(e, nfA, lambdas, ht, z2, prec);
1049 966 : if (P) return P;
1050 931 : if (signe(ell_get_disc(e)) > 0)
1051 : {
1052 819 : z2 = gadd(z2, gmul2n(Oim, -1));
1053 819 : P = heegner_try_point(e, nfA, lambdas, ht, z2, prec);
1054 819 : if (P) return P;
1055 : }
1056 924 : set_avma(av);
1057 : }
1058 0 : pari_err_BUG("ellheegner, point not found");
1059 : return NULL; /* LCOV_EXCL_LINE */
1060 : }
1061 :
1062 : /* N > 1, fa = factor(N), return factor(4*N) */
1063 : static GEN
1064 42 : fa_shift2(GEN fa)
1065 : {
1066 42 : GEN P = gel(fa,1), E = gel(fa,2);
1067 42 : if (absequaliu(gcoeff(fa,1,1), 2))
1068 : {
1069 21 : E = shallowcopy(E);
1070 21 : gel(E,1) = addiu(gel(E,1), 2);
1071 : }
1072 : else
1073 : {
1074 21 : P = shallowconcat(gen_2, P);
1075 21 : E = shallowconcat(gen_2, E);
1076 : }
1077 42 : return mkmat2(P, E);
1078 : }
1079 :
1080 : /* P = prime divisors of N(E). Return the product of primes p in P, a_p != -1
1081 : * HACK: restrict to small primes since large ones won't divide our C-long
1082 : * discriminants */
1083 : static GEN
1084 70 : get_bad(GEN E, GEN F)
1085 : {
1086 70 : GEN P = gel(F,1), e = gel(F,2);
1087 70 : long k, l = lg(P), ibad = 1;
1088 70 : GEN B = cgetg(l, t_VECSMALL);
1089 224 : for (k = 1; k < l; k++)
1090 154 : if (is_pm1(gel(e,k)))
1091 : {
1092 63 : GEN p = gel(P,k);
1093 63 : long pp = itos_or_0(p);
1094 63 : if (!pp) break;
1095 63 : if (! equalim1(ellap(E,p))) B[ibad++] = pp;
1096 : }
1097 70 : setlg(B, ibad); return ibad == 1? NULL: zv_prod_Z(B);
1098 : }
1099 :
1100 : /* factorization into primary factors */
1101 : static GEN
1102 42 : to_primary(GEN fa)
1103 : {
1104 42 : GEN Q, P = gel(fa,1), E = gel(fa,2);
1105 42 : long i, l = lg(P);
1106 42 : Q = cgetg(l, t_COL);
1107 147 : for (i = 1; i < l; i++) gel(Q,i) = powii(gel(P,i), gel(E,i));
1108 42 : return mkmat2(Q, const_col(l-1, gen_1));
1109 : }
1110 :
1111 : /* list of pairs [Q,N/Q], where Q || N, sorted by increasing Q. */
1112 : static GEN
1113 42 : find_div(GEN faN)
1114 : {
1115 42 : GEN L, Q = divisors(to_primary(faN));
1116 42 : long k, l = lg(Q);
1117 42 : L = cgetg(l, t_VEC);
1118 322 : for (k = 1; k < l; k++) gel(L, k) = mkvec2(gel(Q,k), gel(Q,l-k));
1119 42 : return L;
1120 : }
1121 :
1122 : static long
1123 8736 : testDisc(GEN bad, long d) { return !bad || ugcdiu(bad, -d) == 1; }
1124 : /* bad = product of bad primes. Return the NDISC largest fundamental
1125 : * discriminants D < d such that (D,bad) = 1 and d is a square mod 4N */
1126 : static GEN
1127 42 : listDisc(GEN fa4N, GEN bad, long d, long ndisc)
1128 : {
1129 42 : GEN v = cgetg(ndisc+1, t_VECSMALL);
1130 42 : pari_sp av = avma;
1131 42 : long j = 1;
1132 : for(;;)
1133 : {
1134 8652 : d -= odd(d)? 1: 3;
1135 8652 : if (testDisc(bad,d) && unegisfundamental(-d) && Zn_issquare(stoi(d), fa4N))
1136 : {
1137 420 : v[j++] = d;
1138 420 : if (j > ndisc) break;
1139 : }
1140 8610 : set_avma(av);
1141 : }
1142 42 : set_avma(av); return v;
1143 : }
1144 :
1145 : /* as normforms, but only return solution so that b = beta [mod Q] */
1146 : static GEN
1147 112483 : normformsbeta(GEN D, GEN fa, GEN beta, GEN Q)
1148 : {
1149 112483 : pari_sp av = avma;
1150 : long i, j, k, lB, aN, sa;
1151 : GEN a, L, V, B, N, N2;
1152 112483 : int D_odd = mpodd(D);
1153 112483 : a = typ(fa) == t_INT ? fa: typ(fa) == t_VEC? gel(fa,1): factorback(fa);
1154 112483 : sa = signe(a);
1155 112483 : if (sa==0 || (signe(D)<0 && sa<0)) return NULL;
1156 5530 : V = D_odd? Zn_quad_roots(fa, gen_1, shifti(subsi(1, D), -2))
1157 112483 : : Zn_quad_roots(fa, gen_0, negi(shifti(D, -2)));
1158 112483 : if (!V) return NULL;
1159 32403 : N = gel(V,1); B = gel(V,2); lB = lg(B);
1160 32403 : N2 = shifti(N,1);
1161 32403 : aN = itou(diviiexact(a, N)); /* |a|/N */
1162 32403 : L = cgetg((lB-1)*aN+1, t_VEC);
1163 188657 : for (k = 1, i = 1; i < lB; i++)
1164 : {
1165 156254 : GEN b = shifti(gel(B,i), 1), c, C;
1166 156254 : if (D_odd) b = addiu(b, 1);
1167 156254 : c = diviiexact(shifti(subii(sqri(b), D), -2), a);
1168 156254 : for (j = 0;; b = addii(b, N2))
1169 : {
1170 172382 : if (dvdii(subii(b,beta),Q))
1171 105833 : gel(L, k++) = mkqfb(a, b, c, D);
1172 172382 : if (++j == aN) break;
1173 16128 : C = addii(b, N); if (aN > 1) C = diviuexact(C, aN);
1174 16128 : c = sa > 0? addii(c, C): subii(c, C);
1175 : }
1176 : }
1177 32403 : setlg(L, k); return gc_GEN(av, L);
1178 : }
1179 :
1180 : /* faN4 = factor(4*N) */
1181 : static GEN
1182 420 : listheegner(GEN N, GEN faN4, GEN listQ, GEN D)
1183 : {
1184 420 : pari_sp av = avma;
1185 : hashtable H;
1186 420 : ulong h = itou(quadclassno(D));
1187 420 : GEN ymin, beta = Zn_sqrt(D, faN4), N2 = shifti(N, 1), L, V;
1188 : long l, k, n;
1189 420 : hash_init_GEN(&H, h, gequal, 1);
1190 5740 : for (k = 1; H.nb < h ; k++)
1191 : {
1192 5320 : GEN LF = normformsbeta(D, mulis(N, k), beta, N2);
1193 5320 : if (LF)
1194 : {
1195 3990 : long i, l = lg(LF);
1196 8953 : for (i = 1; i < l && H.nb < h; i++)
1197 : {
1198 4963 : GEN F = gel(LF, i), f = qfi_red(F), k = hash_haskey_GEN(&H, f);
1199 4963 : if (!k)
1200 : {
1201 2478 : GEN a = gel(F,1), b = gel(F,2), c = gel(F,3);
1202 2478 : GEN F2 = mkqfb(diviiexact(a,N), negi(b), mulii(c,N), D);
1203 2478 : GEN f2 = qfi_red(F2);
1204 2478 : if (gequal(f, f2))
1205 322 : hash_insert(&H, f, mkvec(F));
1206 : else
1207 : {
1208 2156 : hash_insert(&H, f, mkvec2(F, F2));
1209 2156 : hash_insert(&H, f2, cgetg(1,t_VEC));
1210 : }
1211 : }
1212 : }
1213 : }
1214 : }
1215 420 : ymin = NULL;
1216 420 : V = hash_values_GEN(&H); l = lg(V);
1217 420 : L = cgetg(H.nb+1, t_VEC);
1218 5054 : for (k = n = 1; k < l; k++)
1219 : {
1220 4634 : GEN t, Q, Vk = gel(V,k), y;
1221 4634 : long nk = lg(Vk) - 1;
1222 4634 : if (nk)
1223 : {
1224 2478 : y = lift_points(listQ, gel(Vk,1), &t, &Q);
1225 2478 : gel(L, n++) = mkvec3(t, stoi(nk), Q);
1226 2478 : if (!ymin || gcmp(y, ymin) < 0) ymin = y;
1227 : }
1228 : }
1229 420 : setlg(L, n); /* H.nb/2 <= n-1 <= H.nb */
1230 420 : if (DEBUGLEVEL > 1)
1231 0 : err_printf("Disc %Ps : N*ymin = %Pg\n", D,
1232 : gmul(gsqrt(ymin, DEFAULTPREC),N));
1233 420 : return gc_GEN(av, mkvec3(ymin, L, D));
1234 : }
1235 :
1236 : /* Q | N, P = prime divisors of N, R[i] = local epsilon-factor at P[i].
1237 : * Return \prod_{p | Q} R[i] */
1238 : static long
1239 147 : rootno(GEN Q, GEN P, GEN R)
1240 : {
1241 147 : long s = 1, i, l = lg(P);
1242 581 : for (i = 1; i < l; i++)
1243 434 : if (dvdii(Q, gel(P,i))) s *= R[i];
1244 147 : return s;
1245 : }
1246 :
1247 : static void
1248 42 : heegner_find_disc(GEN *points, GEN *coefs, long *pind, GEN E,
1249 : GEN indmult, long ndisc, long prec)
1250 : {
1251 42 : long d = 0;
1252 : GEN faN4, bad, N, faN, listQ, listR;
1253 :
1254 42 : ellQ_get_Nfa(E, &N, &faN);
1255 42 : faN4 = fa_shift2(faN);
1256 42 : listQ = find_div(faN);
1257 42 : bad = get_bad(E, faN);
1258 42 : listR = gel(obj_check(E, Q_ROOTNO), 2);
1259 : for(;;)
1260 0 : {
1261 42 : pari_sp av = avma;
1262 42 : GEN list, listD = listDisc(faN4, bad, d, ndisc);
1263 42 : long k, l = lg(listD);
1264 42 : list = cgetg(l, t_VEC);
1265 462 : for (k = 1; k < l; ++k)
1266 420 : gel(list, k) = listheegner(N, faN4, listQ, stoi(listD[k]));
1267 42 : list = vecsort0(list, gen_1, 0);
1268 56 : for (k = l-1; k > 0; --k)
1269 : {
1270 56 : long bprec = 8;
1271 56 : GEN Lk = gel(list,k), D = gel(Lk,3);
1272 56 : GEN sqrtD = sqrtr_abs(itor(D, prec)); /* sqrt(|D|) */
1273 56 : GEN indmultD = heegner_indexmultD(faN, indmult, itos(D), sqrtD);
1274 : do
1275 : {
1276 : GEN mulf, indr;
1277 : pari_timer ti;
1278 56 : if (DEBUGLEVEL) timer_start(&ti);
1279 56 : mulf = ltwist1(E, D, bprec+expo(indmultD));
1280 56 : if (DEBUGLEVEL) timer_printf(&ti,"ellL1twist");
1281 56 : indr = mulrr(indmultD, mulf);
1282 56 : if (DEBUGLEVEL) err_printf("Disc = %Ps, Index^2 = %Ps\n", D, indr);
1283 56 : if (signe(indr)>0 && expo(indr) >= -1) /* indr >=.5 */
1284 : {
1285 : long e, i, l;
1286 42 : GEN pts, cfs, L, indi = grndtoi(sqrtr_abs(indr), &e);
1287 42 : if (e > expi(indi)-7)
1288 : {
1289 0 : bprec++;
1290 0 : pari_warn(warnprec, "ellL1",bprec);
1291 0 : continue;
1292 : }
1293 42 : *pind = itos(indi);
1294 42 : L = gel(Lk, 2); l = lg(L);
1295 42 : pts = cgetg(l, t_VEC);
1296 42 : cfs = cgetg(l, t_VECSMALL);
1297 189 : for (i = 1; i < l; ++i)
1298 : {
1299 147 : GEN P = gel(L,i), z = gel(P,2), Q = gel(P,3); /* [1 or 2, Q] */
1300 : long c;
1301 147 : gel(pts, i) = qfb_root(gel(P,1), sqrtD);
1302 147 : c = rootno(Q, gel(faN,1), listR);
1303 147 : if (!equali1(z)) c *= 2;
1304 147 : cfs[i] = c;
1305 : }
1306 42 : if (DEBUGLEVEL)
1307 0 : err_printf("N = %Ps, ymin*N = %Ps\n",N,
1308 0 : gmul(gsqrt(gel(Lk, 1), prec),N));
1309 42 : *coefs = cfs; *points = pts; return;
1310 : }
1311 : } while(0);
1312 : }
1313 0 : d = listD[l-1]; set_avma(av);
1314 : }
1315 : }
1316 :
1317 : GEN
1318 224 : ellanal_globalred_all(GEN e, GEN *cb, GEN *N, GEN *tam)
1319 : {
1320 224 : GEN E = ellanal_globalred(e, cb), red = obj_check(E, Q_GLOBALRED);
1321 224 : *N = gel(red, 1);
1322 224 : if (tam)
1323 : {
1324 154 : *tam = gel(red,2);
1325 154 : if (signe(ell_get_disc(E))>0) *tam = shifti(*tam,1);
1326 : }
1327 224 : return E;
1328 : }
1329 :
1330 : static GEN
1331 42 : vecelnfembed(GEN x, GEN M, GEN et)
1332 91 : { pari_APPLY_same(gmul(M, etnf_to_basis(et, gel(x,i)))) }
1333 :
1334 : static GEN
1335 42 : QXQV_inv(GEN x, GEN T)
1336 91 : { pari_APPLY_same(QXQ_inv(gel(x,i), T)) }
1337 :
1338 : static GEN
1339 42 : etnfnewprec(GEN x, long prec)
1340 126 : { pari_APPLY_same(nfnewprec(gel(x,i),prec)) }
1341 :
1342 : static GEN
1343 49 : vec_etnf_to_basis(GEN et, GEN x)
1344 154 : { pari_APPLY_same(etnf_to_basis(et,gel(x,i))) }
1345 :
1346 : static GEN
1347 42 : etnf_M(GEN et)
1348 : {
1349 42 : long i, l = lg(et);
1350 42 : GEN V = cgetg(l, t_VEC);
1351 126 : for (i = 1; i < l; i++) gel(V,i) = nf_get_M(gel(et,i));
1352 42 : return shallowmatconcat(diagonal_shallow(V));
1353 : }
1354 :
1355 : static GEN
1356 42 : makenfA(GEN sel, GEN A, GEN cb)
1357 : {
1358 42 : GEN etal = gel(sel,1), T = gel(etal,3);
1359 42 : GEN et = gel(etal,1), M = etnf_M(et);
1360 42 : GEN al = vecelnfembed(A, M, et);
1361 42 : GEN th = gmul(M, etnf_to_basis(et, pol_x(etnf_get_varn(et))));
1362 42 : return mkvec6(etal,QXQV_inv(A, T),cb,al,th,M);
1363 : }
1364 :
1365 : static GEN
1366 42 : ellheegner_z_i(GEN E, double ht, long *pindx,
1367 : GEN N, GEN tam, long wtor, long etor, long prec)
1368 : {
1369 42 : GEN om = ellR_omega(E,prec), w1 = gel(om,1), w2 = gel(om,2);
1370 : GEN pts, cfs, s, z, indmult;
1371 42 : long ndisc = maxss(10, (long)(ht? ht: prec*M_LN2) / 10), ind, lint, k, l;
1372 : pari_timer ti;
1373 :
1374 : /* O_re*c(E) / (4*O_vol*|E_t|^2) */
1375 42 : indmult = divir(tam, mulru(gel(w2,2), 4*wtor*wtor)); setsigne(indmult,1);
1376 42 : heegner_find_disc(&pts, &cfs, &ind, E, indmult, ndisc, prec);
1377 :
1378 42 : if (DEBUGLEVEL) timer_start(&ti);
1379 42 : s = heegner_psi(E, N, pts, prec);
1380 42 : if (DEBUGLEVEL) timer_printf(&ti,"heegner_psi");
1381 :
1382 42 : l = lg(pts); z = mulsr(cfs[1], gel(s, 1));
1383 147 : for (k = 2; k < l; k++) z = addrr(z, mulsr(cfs[k], gel(s, k)));
1384 42 : z = subrr(z, mulri(w1, roundr(divrr(z, w1))));
1385 42 : if (DEBUGLEVEL) err_printf("z=%.*Pg\n", nbits2ndec(prec), z);
1386 :
1387 42 : lint = etor > 1 ? ugcd(ind, etor): 1;
1388 42 : *pindx = 2*lint*ind; return gmulsg(2*lint, z);
1389 : }
1390 :
1391 : /* set cb,N,tam from globalred, w cardinal of E(Q)_tor and e its exponent */
1392 : static GEN
1393 56 : ellheegner_init(GEN E, GEN *pcb, GEN *pN, GEN *ptam, long *pw, long *pe)
1394 : {
1395 : GEN T;
1396 : long w;
1397 :
1398 56 : E = ellanal_globalred_all(E, pcb, pN, ptam);
1399 56 : if (ellrootno_global(E) == 1)
1400 7 : pari_err_DOMAIN("ellheegner", "(analytic rank)%2", "=", gen_0, E);
1401 49 : T = elltors(E); w = itou(abgrp_get_no(T));
1402 49 : *pe = lg(abgrp_get_cyc(T)) == 2? w: (w >> 1);
1403 49 : *pw = w; return E;
1404 : }
1405 :
1406 : GEN
1407 0 : ellheegner_z(GEN E, long prec)
1408 : {
1409 0 : pari_sp av = avma;
1410 : long indx, etor, wtor;
1411 : GEN N, tam, z;
1412 :
1413 0 : E = ellheegner_init(E, NULL, &N, &tam, &wtor, &etor);
1414 0 : z = ellheegner_z_i(E, 0., &indx, N, tam, wtor, etor, prec);
1415 0 : return gerepilecopy(av, mkvec2(z, utoi(indx)));
1416 : }
1417 :
1418 : GEN
1419 56 : ellheegner(GEN E)
1420 : {
1421 56 : pari_sp av = avma;
1422 : GEN N, tam, z, cb, P, ht, om, nfA, sel, etal, et, cbb, sbase, dAi, T, A, Ag;
1423 56 : long bitprec = 16, prec = nbits2prec(bitprec) + EXTRAPRECWORD;
1424 : long indx, wtor, etor, selrank;
1425 : pari_timer ti;
1426 :
1427 56 : E = ellheegner_init(E, &cb, &N, &tam, &wtor, &etor);
1428 49 : sel = ell2selmer_basis(E, &cbb, prec);
1429 49 : etal = gel(sel,1); sbase = gel(sel,2); et = gel(etal,1); T = gel(etal,3);
1430 49 : dAi = gsupnorm(vec_etnf_to_basis(et,sbase),prec);
1431 : while (1)
1432 42 : {
1433 : GEN hnaive, l1;
1434 : long bitneeded;
1435 91 : if (DEBUGLEVEL) timer_start(&ti);
1436 91 : l1 = ellL1(E, 1, bitprec);
1437 91 : if (DEBUGLEVEL) timer_printf(&ti,"ellL1");
1438 91 : if (expo(l1) < 1 - bitprec/2)
1439 7 : pari_err_DOMAIN("ellheegner", "analytic rank",">",gen_1,E);
1440 84 : om = ellR_omega(E,prec);
1441 84 : ht = divrr(mulru(l1, wtor * wtor), mulri(gel(om,1), tam));
1442 84 : if (DEBUGLEVEL) err_printf("Expected height=%Ps\n", ht);
1443 84 : hnaive = hnaive_max(E, ht);
1444 84 : if (DEBUGLEVEL) err_printf("Naive height <= %Ps\n", hnaive);
1445 84 : hnaive = gadd(shiftr(hnaive,-1), glog(dAi,prec));
1446 84 : bitneeded = itos(gceil(divrr(hnaive, mplog2(prec)))) + 32;
1447 84 : if (DEBUGLEVEL) err_printf("precision = %ld\n", bitneeded);
1448 84 : if (bitprec >= bitneeded) break;
1449 42 : bitprec = bitneeded;
1450 42 : prec = nbits2prec(bitprec) + EXTRAPRECWORD;
1451 : }
1452 42 : z = ellheegner_z_i(E, rtodbl(ht), &indx, N, tam, wtor, etor, prec);
1453 :
1454 42 : selrank = lg(sbase)-1;
1455 42 : Ag = selrank > lg(et)-1 ? pol_1(etnf_get_varn(et)): gel(sbase,selrank);
1456 42 : if (vals(indx) >= vals(etor))
1457 35 : A = mkvec(Ag);
1458 : else
1459 7 : A = mkvec2(Ag, QXQ_mul(Ag, gel(sbase,1), T));
1460 42 : gmael(sel,1,1) = etnfnewprec(et, prec);
1461 42 : nfA = makenfA(sel, A, cbb);
1462 42 : if (DEBUGLEVEL) timer_start(&ti);
1463 42 : P = heegner_find_point(E, nfA, om, ht, z, indx, prec);
1464 42 : if (DEBUGLEVEL) timer_printf(&ti,"heegner_find_point");
1465 42 : if (cb) P = ellchangepointinv(P, cb);
1466 42 : return gc_GEN(av, P);
1467 : }
1468 :
1469 : /* ellheegnertwist */
1470 :
1471 : static GEN
1472 105 : listheegnertwist(GEN Q, GEN D, GEN P, GEN Dt, long *kmax)
1473 : {
1474 105 : pari_sp av = avma;
1475 : hashtable H;
1476 105 : GEN beta = Zn_sqrt(D, shifti(Q, 2)), Q2 = shifti(Q, 1);
1477 105 : GEN P2 = sqri(P), DP2 = mulii(D, P2);
1478 105 : long h = itos(quadclassno(DP2));
1479 105 : long k, s = 0;
1480 105 : hash_init_GEN(&H, h, gequal, 1);
1481 107268 : for (k = 1; s < h ; k++)
1482 : {
1483 107163 : GEN LF = normformsbeta(D, mulis(Q, k), beta, Q2);
1484 107163 : if (LF)
1485 : {
1486 28413 : long i, l = lg(LF);
1487 129066 : for (i = 1; i < l && s < h; i++)
1488 : {
1489 100653 : GEN F = gel(LF, i), a = gel(F,1), b = gel(F,2), c = gel(F,3);
1490 100653 : GEN bi = mulii(P,b), ci = mulii(P2, c);
1491 100653 : if (is_pm1(gcdii(gcdii(a, bi), ci)))
1492 : {
1493 65730 : GEN f = qfi_red(mkqfb(a,bi,ci,DP2)), k = hash_haskey_GEN(&H, f);
1494 65730 : long e = kronecker(Dt,a);
1495 65730 : if (!k)
1496 : {
1497 2667 : hash_insert(&H, f, mkvec2(mkvec3(a,b,c),stoi(e)));
1498 2667 : s += !!e;
1499 63063 : } else if (e && signe(gel(k,2))==0)
1500 : {
1501 700 : gel(k,2) = stoi(e); s++;
1502 : }
1503 : }
1504 : }
1505 : }
1506 : }
1507 105 : *kmax = k-1;
1508 105 : return gc_GEN(av, hash_values_GEN(&H));
1509 : }
1510 :
1511 : static GEN
1512 35 : findtwist(GEN E, GEN *pt_D)
1513 : {
1514 35 : GEN d = ellminimaltwistcond(E), N, F, P, Ex, D = *pt_D;
1515 : long i, l;
1516 35 : D = coredisc(mulii(D,d));
1517 35 : E = elltwist(E, d);
1518 35 : ellanal_globalred_all(E, NULL, &N, NULL);
1519 35 : F = Z_factor(gcdii(absi(D),N));
1520 35 : P = gel(F,1); Ex = gel(F,2); l = lg(P);
1521 42 : for (i = 1; i < l; i++)
1522 : {
1523 7 : GEN p = gel(P,i), e = gel(Ex,i), q = powii(p, e);
1524 7 : if (is_pm1(gcdii(q,diviiexact(N,q)))) continue;
1525 0 : d = cmpiu(p,2)>0 ? mod4(p)== 1 ? p : negi(p) : shifti( Mod8(D)==4 ? gen_m1: Mod32(D)==8 ? gen_1:gen_m1,vali(D));
1526 0 : E = elltwist(E, d); D = coredisc(mulii(D,d));
1527 : }
1528 35 : *pt_D = D;
1529 35 : return ellminimalmodel(E,NULL);
1530 : }
1531 :
1532 : /* ellheegnertwist */
1533 :
1534 : static GEN
1535 770 : finddivis(GEN E, GEN T, GEN z, long n, long prec)
1536 : {
1537 770 : pari_sp av = avma;
1538 770 : long i, j, k, l = lg(T);
1539 770 : GEN om = ellR_omega(E, prec);
1540 770 : GEN om1 = gdivgs(gel(om,1),n), om2 = gmul2n(gel(om,2),-1);
1541 770 : z = gdivgs(z, n); T = gdivgs(T,n);
1542 46249 : for (i = 0; i < n; i++)
1543 : {
1544 135709 : for (j = 1; j < l; j++)
1545 : {
1546 90230 : GEN w = gadd(z,gel(T,j));
1547 270277 : for(k = 0; k < (odd(n)?1:2); k++)
1548 : {
1549 180082 : GEN x2, y2, P = pointell(E,w,prec);
1550 180082 : if (ell_is_inf(P)) return P;
1551 180082 : x2 = bestappr(real_i(gel(P,1)), int2n((prec>>1)-1));
1552 180082 : if (lg(x2)==1) continue;
1553 180082 : y2 = ellordinate(E, x2, prec);
1554 180082 : if (lg(y2)>1) return gc_GEN(av, mkvec2(x2, gel(y2,1)));
1555 180047 : w = gadd(z, om2);
1556 : }
1557 : }
1558 45479 : z = gadd(z, om1);
1559 : }
1560 735 : return gc_NULL(av);
1561 : }
1562 :
1563 : static GEN
1564 147 : listpointstwist(GEN D, GEN L, long prec)
1565 : {
1566 147 : GEN vDi = gsqrt(D, prec);
1567 : long k, l;
1568 147 : GEN V = cgetg_copy(L, &l);
1569 3360 : for (k = 1; k < l; k++)
1570 : {
1571 3213 : GEN Lk = gel(L,k), v = gel(Lk, 1);
1572 3213 : gel(V,k) = mkvec2(gdiv(gsub(vDi,gel(v,2)),shifti(gel(v,1),1)),gel(Lk,2));
1573 : }
1574 147 : return V;
1575 : }
1576 :
1577 : static GEN
1578 5836413 : _gen_cmul(void *E, GEN an, long i, GEN x)
1579 : {
1580 5836413 : (void) E; return i ? gdivgs(gmulgs(x, an[i]), i): gen_0;
1581 : }
1582 :
1583 : static GEN
1584 2505 : heegnersum(GEN an, long nb, GEN z, long r)
1585 : {
1586 2505 : GEN q = gen_bkeval(an, nb, z, 1, NULL, get_Rg_algebra(), _gen_cmul);
1587 2506 : return r>0 ? real_i(q) : imag_i(q);
1588 : }
1589 :
1590 : GEN
1591 2503 : heegnersum_worker(GEN l, GEN an, long r, long prec)
1592 : {
1593 2503 : pari_sp av = avma;
1594 2503 : long nb = ceil(prec*M_LN2/(2*M_PI*rtodbl(imag_i(gel(l,1)))));
1595 2502 : return gc_upto(av, gmul(gel(l,2), heegnersum(an, nb, expIPiC(gmul2n(gel(l,1),1), prec), r)));
1596 : }
1597 :
1598 : static GEN
1599 112 : sumheegner(GEN an, GEN L, long r, long prec)
1600 : {
1601 112 : pari_sp av = avma;
1602 112 : GEN worker = snm_closure(is_entry("_heegnersum_worker"),mkvec3(an,stoi(r),stoi(prec)));
1603 : GEN V;
1604 112 : if (DEBUGLEVEL>1) err_printf("ellheegnertwist: Computing sum with prec %ld: ",prec);
1605 112 : V = gen_parapply_percent(worker, L, DEBUGLEVEL>1);
1606 112 : if (DEBUGLEVEL>1) err_printf(" done.\n");
1607 112 : return gc_upto(av, vecsum(V));
1608 : }
1609 :
1610 : static GEN
1611 112 : torstoz(GEN E, GEN ET, long prec)
1612 : {
1613 112 : long n = itos(abgrp_get_no(ET));
1614 : GEN cyc, gen, g1, g2;
1615 : long o1, i;
1616 112 : GEN V = cgetg(n+1, t_VEC);
1617 112 : gel(V,1) = gen_0;
1618 112 : if (n == 1) return V;
1619 63 : cyc = abgrp_get_cyc(ET);
1620 63 : gen = abgrp_get_gen(ET);
1621 63 : g1 = zell(E, gel(gen,1), prec);
1622 63 : o1 = itos(gel(cyc,1));
1623 126 : for (i = 2; i <= o1; i++)
1624 63 : gel(V,i) = gmulgs(g1,i-1);
1625 63 : if (lg(cyc)==2) return V;
1626 7 : g2 = zell(E, gel(gen,2), prec);
1627 21 : for (i = o1+1; i <=n; i++)
1628 14 : gel(V,i) = gadd(g2, gel(V,i-o1));
1629 7 : return V;
1630 : }
1631 :
1632 : static GEN
1633 35 : findpoint(GEN E, GEN N, GEN D, GEN d, GEN f, long manin, GEN s, long prec)
1634 : {
1635 : long kmax, m;
1636 : GEN AL, LH, Ed, ET, tam, faN, R;
1637 : double mR;
1638 35 : AL = diviiexact(N,f);
1639 35 : if (!Zn_issquare(D,shifti(AL,2))) pari_err_BUG("findpoint");
1640 35 : LH = listheegnertwist(AL, D, f, d, &kmax);
1641 35 : Ed = elltwist(E, d); ET= elltors(Ed);
1642 35 : m = ellmanintable_heuristic(Ed);
1643 35 : tam = elltamagawa(Ed); faN = Z_factor(mulii(N,sqri(d)));
1644 35 : R = imag_i(row(listpointstwist(D, LH, MEDDEFAULTPREC),1));
1645 35 : mR = M_LN2/(2*M_PI*rtodbl(vecmin(R)));
1646 35 : if (DEBUGLEVEL > 2)
1647 0 : err_printf("Heegnertwist: min imag: %.9g, h = %ld\n",mR, lg(LH)-1);
1648 77 : for (;;prec*=2)
1649 77 : {
1650 : long i, lD;
1651 112 : GEN T = torstoz(Ed, ET, prec), dv, M, z;
1652 112 : GEN L = listpointstwist(D, LH, prec+EXTRAPREC64);
1653 112 : long nb = ceil(prec*mR);
1654 112 : if (DEBUGLEVEL > 2) err_printf("Heegnertwist: nb = %ld, prec = %ld\n",nb, prec);
1655 112 : z = sumheegner(ellanQ_zv(E, nb), L, signe(d), prec);
1656 112 : z = gmulgs(gdiv(z, gsqrt(absi(d), prec)), 2*manin);
1657 112 : M = mulii(shifti(tam, omega_N_D(faN,-itos(D))), mulsi(m,s));
1658 112 : M = mulis(M, get_w(itos_or_0(d)));
1659 112 : dv = divisors(M); lD = lg(dv);
1660 847 : for (i = 1; i < lD; i++)
1661 : {
1662 770 : ulong di = itou_or_0(gel(dv,i));
1663 770 : if (di)
1664 : {
1665 770 : GEN P = finddivis(Ed, T, z, di, prec);
1666 770 : if (!P) continue;
1667 35 : if (!signe(ellorder(Ed, P, NULL)))
1668 35 : return P;
1669 0 : else if (i==1) return NULL;
1670 : }
1671 : }
1672 : }
1673 : }
1674 :
1675 : static GEN
1676 35 : heegnertwistdisc(GEN *ptD, GEN N, GEN d, GEN f, GEN g, GEN F, GEN bad, long n)
1677 : {
1678 35 : GEN D = *ptD, p, f2 = sqri(f), AL = diviiexact(N,f);
1679 35 : GEN L = cgetg(n+1, t_VEC);
1680 35 : GEN h = cgetg(n+1, t_VEC);
1681 35 : GEN lh = cgetg(n+1, t_VEC);
1682 : long i, kmax;
1683 105 : for (i = 1; i <= n; i++)
1684 : {
1685 : for (;;)
1686 : {
1687 224 : D = subii(D, g);
1688 224 : if (Mod4(D)<=1 && Zn_issquare(D, F) && testDisc(bad, itos(D)))
1689 : {
1690 84 : GEN r, q = dvmdii(mulii(D, f2), d, &r);
1691 84 : if (signe(r)==0 && Mod4(q)<=1) break;
1692 : }
1693 : }
1694 70 : gel(L,i) = D;
1695 70 : gel(lh,i)= listheegnertwist(AL, D, f, d, &kmax);
1696 70 : gel(h,i) = gdiv(gdivgs(D,kmax),gsqr(quadclassno(D)));
1697 : }
1698 35 : *ptD = D;
1699 35 : p = indexsort(h);
1700 35 : return vecpermute(L, p);
1701 : }
1702 :
1703 : /* Algorithm by BA inspired by
1704 : Mock Heegner Points and Congruent Numbers
1705 : Paul Monsky
1706 : Mathematische Zeitschrift (1990) Volume: 204, Issue: 1, page 45-68
1707 : ISSN: 0025-5874; 1432-1823
1708 : https://gdz.sub.uni-goettingen.de/id/PPN266833020_0204
1709 : https://link.springer.com/article/10.1007/BF02570859
1710 :
1711 : Regulators of rank one quadratic twists
1712 : Christophe Delaunay, Xavier-Francois Roblot
1713 : Journal de theorie des nombres de Bordeaux, Tome 20 (2008) no. 3, pp. 601-624
1714 : https://jtnb.centre-mersenne.org/item/10.5802/jtnb.643.pdf
1715 : */
1716 :
1717 : GEN
1718 35 : ellheegnertwist(GEN E, GEN d, GEN s)
1719 : {
1720 35 : pari_sp av = avma;
1721 35 : GEN Et, Ed, iso, N, D1 = gen_0, P, f, f2, g, F, bad, F2;
1722 35 : long m, prec = DEFAULTPREC;
1723 35 : checkell(E);
1724 35 : if (!s) s = gen_1;
1725 7 : else if (typ(s)!=t_INT || signe(s)<=0)
1726 0 : pari_err_TYPE("ellheegnertwist",s);
1727 35 : if (d)
1728 : {
1729 21 : if (typ(d) != t_INT)
1730 0 : pari_err_TYPE("ellheegnertwist",d);
1731 21 : Et = elltwist(E,d);
1732 : }
1733 14 : else { d = gen_1; Et = E; }
1734 35 : if (ellrootno_global(Et) == 1)
1735 0 : pari_err_DOMAIN("ellheegnertwist", "(analytic rank)%2","=",gen_0,E);
1736 35 : E = findtwist(E, &d);
1737 35 : if (DEBUGLEVEL>1) err_printf("ellheegnertwist: using disc. %Ps\n",d);
1738 35 : m = ellmanintable_heuristic(E);
1739 35 : Ed = elltwist(E, d);
1740 35 : ellanal_globalred_all(E, NULL, &N, NULL);
1741 35 : iso = ellisisom(Et, Ed);
1742 35 : f = gcdii(N,d); f2 = sqri(f);
1743 35 : g = diviiexact(absi(d), gcdii(d, f2));
1744 35 : F = Z_factor(diviiexact(N, f));
1745 35 : F2 = famat_reduce(famat_mul(mkmat2(mkcol(gen_2), mkcol(gen_2)),F));
1746 35 : bad = equalii(N,f) ? NULL: get_bad(Ed, F);
1747 : for(;;)
1748 0 : {
1749 35 : const long n = 2;
1750 35 : GEN L = heegnertwistdisc(&D1, N, d, f, g, F2, bad, n);
1751 : long i;
1752 35 : for(i = 1; i <= n; i++)
1753 : {
1754 35 : GEN Di = gel(L,i);
1755 35 : if (DEBUGLEVEL > 1)
1756 0 : err_printf("ellheegnertwist: using secondary disc. %Ps\n", Di);
1757 35 : P = findpoint(E, N, Di, d, f, m, s, prec);
1758 35 : if (P) return gc_GEN(av, ellchangepointinv(P, iso));
1759 0 : if (DEBUGLEVEL)
1760 0 : err_printf("ellheegnertwist: point is torsion for D=%Ps\n", Di);
1761 : }
1762 : }
1763 : }
1764 :
1765 : /* Modular degree */
1766 :
1767 : static GEN
1768 70 : ellisobound(GEN e)
1769 : {
1770 70 : GEN M = gel(ellisomat(e,0,1),2);
1771 70 : return vecmax(gel(M,1));
1772 : }
1773 : /* 4Pi^2 / E.area */
1774 : static GEN
1775 140 : getA(GEN E, long prec) { return divrr(sqrr(Pi2n(1,prec)), ellR_area(E, prec)); }
1776 :
1777 : /* Modular degree of elliptic curve e over Q, assuming Manin constant = 1
1778 : * (otherwise multiply by square of Manin constant). */
1779 : GEN
1780 70 : ellmoddegree(GEN E)
1781 : {
1782 70 : pari_sp av = avma;
1783 : GEN N, tam, mc2, d;
1784 : long b;
1785 70 : E = ellanal_globalred_all(E, NULL, &N, &tam);
1786 70 : mc2 = sqri(ellisobound(E));
1787 70 : b = expi(mulii(N, mc2)) + maxss(0, expo(getA(E, LOWDEFAULTPREC))) + 16;
1788 : for(;;)
1789 0 : {
1790 70 : long prec = nbits2prec(b), e, s;
1791 70 : GEN deg = mulri(mulrr(lfunellmfpeters(E, b), getA(E, prec)), mc2);
1792 70 : d = grndtoi(deg, &e);
1793 70 : if (DEBUGLEVEL) err_printf("ellmoddegree: %Ps, bit=%ld, err=%ld\n",deg,b,e);
1794 70 : s = expo(deg);
1795 70 : if (e <= -8 && s <= b-8) return gc_upto(av, gdiv(d,mc2));
1796 0 : b = maxss(s, b+e) + 16;
1797 : }
1798 : }
|
