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 1528191 : gen_BG_add(void *E, bg_fun *fun, struct bg_data *bg, GEN n, long i, GEN a, GEN lasta)
50 : {
51 1528191 : pari_sp av = avma;
52 : long j;
53 1528191 : ulong nn = itou_or_0(n);
54 1528191 : if (nn && nn <= bg->rootbnd) bg->an[nn] = itos(a);
55 :
56 1528191 : if (signe(a))
57 : {
58 367465 : fun(E, n, a);
59 367465 : j = 1;
60 : }
61 : else
62 1160726 : j = i;
63 3051076 : for(; j <= i; j++)
64 : {
65 2420509 : ulong p = bg->p[j];
66 2420509 : GEN nexta, pn = mului(p, n);
67 2420509 : if (cmpii(pn, bg->bnd) > 0) return;
68 1522885 : nexta = mulis(a, bg->an[p]);
69 1522885 : if (i == j && umodiu(bg->N, p)) nexta = subii(nexta, mului(p, lasta));
70 1522885 : gen_BG_add(E, fun, bg, pn, j, nexta, a);
71 1522885 : set_avma(av);
72 : }
73 : }
74 :
75 : static void
76 77 : gen_BG_init(struct bg_data *bg, GEN E, GEN N, GEN bnd)
77 : {
78 77 : bg->E = E;
79 77 : bg->N = N;
80 77 : bg->bnd = bnd;
81 77 : bg->rootbnd = itou(sqrtint(bnd));
82 77 : bg->p = primes_upto_zv(bg->rootbnd);
83 77 : bg->an = ellanQ_zv(E, bg->rootbnd);
84 77 : }
85 :
86 : static void
87 77 : gen_BG_rec(void *E, bg_fun *fun, struct bg_data *bg)
88 : {
89 77 : long i, j, lp = lg(bg->p)-1;
90 77 : GEN bndov2 = shifti(bg->bnd, -1);
91 77 : pari_sp av = avma, av2;
92 : GEN p;
93 : forprime_t S;
94 77 : (void)forprime_init(&S, utoipos(bg->p[lp]+1), bg->bnd);
95 77 : av2 = avma;
96 77 : if (DEBUGLEVEL)
97 0 : err_printf("1st stage, using recursion for p <= %ld\n", bg->p[lp]);
98 5383 : for (i = 1; i <= lp; i++)
99 : {
100 5306 : ulong pp = bg->p[i];
101 5306 : long ap = bg->an[pp];
102 5306 : gen_BG_add(E, fun, bg, utoipos(pp), i, stoi(ap), gen_1);
103 5306 : set_avma(av2);
104 : }
105 77 : if (DEBUGLEVEL) err_printf("2nd stage, looping for p <= %Ps\n", bndov2);
106 1113672 : while ( (p = forprime_next(&S)) )
107 : {
108 : long jmax;
109 1113672 : GEN ap = ellap(bg->E, p);
110 1113672 : pari_sp av3 = avma;
111 1113672 : if (!signe(ap)) continue;
112 :
113 556178 : jmax = itou( divii(bg->bnd, p) ); /* 2 <= jmax <= el->rootbound */
114 556178 : fun(E, p, ap);
115 8337322 : for (j = 2; j <= jmax; j++)
116 : {
117 7781144 : long aj = bg->an[j];
118 : GEN a, n;
119 7781144 : if (!aj) continue;
120 1136044 : a = mulis(ap, aj);
121 1136044 : n = muliu(p, j);
122 1136044 : fun(E, n, a);
123 1136044 : set_avma(av3);
124 : }
125 556178 : set_avma(av2);
126 556178 : if (abscmpii(p, bndov2) >= 0) break;
127 : }
128 77 : if (DEBUGLEVEL) err_printf("3nd stage, looping for p <= %Ps\n", bg->bnd);
129 997157 : while ( (p = forprime_next(&S)) )
130 : {
131 997080 : GEN ap = ellap(bg->E, p);
132 997080 : if (!signe(ap)) continue;
133 497896 : fun(E, p, ap);
134 497896 : set_avma(av2);
135 : }
136 77 : set_avma(av);
137 77 : }
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 210 : vecF2_lk(GEN E, GEN K, GEN rbnd, GEN Q, GEN sleh, long prec)
185 : {
186 : pari_sp av;
187 210 : long l, L = lg(K)-1;
188 210 : GEN a = ellanQ_zv(E, itos(gel(K,1)));
189 210 : GEN S = cgetg(L+1, t_VEC);
190 :
191 12131 : for (l = 1; l <= L; l++) gel(S,l) = cgetr(prec);
192 210 : av = avma;
193 12131 : for (l = 1; l <= L; l++)
194 : {
195 : GEN e1, Sl, z, zB;
196 11921 : 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 11921 : e1 = gel(Q, l);
200 11921 : Sl = real_0(prec);;
201 : /* baby-step giant step */
202 11921 : B = A = rbnd[l];
203 11921 : z = powersr(e1, B); zB = gel(z, B+1);
204 11921 : av2 = avma;
205 223293 : for (aB = A*B; aB >= 0; aB -= B)
206 : {
207 211372 : GEN s = real_0(prec); /* could change also prec here */
208 14834204 : for (b = B; b > 0; --b)
209 : {
210 14622832 : long k = aB+b;
211 14622832 : if (k <= Kl && a[k]) s = addrr(s, mulsr(a[k], gel(z, b+1)));
212 14622832 : if (gc_needed(av2, 1)) gerepileall(av2, 2, &s, &Sl);
213 : }
214 211372 : Sl = addrr(mulrr(Sl, zB), s);
215 : }
216 11921 : affrr(mulrr(Sl, gel(sleh,l)), gel(S, l)); /* to avoid copying all S */
217 11921 : set_avma(av);
218 : }
219 210 : return S;
220 : }
221 :
222 : /* Return C, C[i][j] = Q[j]^i, i = 1..nb */
223 : static void
224 35 : baby_init(struct baby_giant *bb, GEN Q, GEN bnd, GEN rbnd, long prec)
225 : {
226 35 : long i, j, l = lg(Q);
227 : GEN R, C, r0;
228 35 : C = cgetg(l,t_VEC);
229 1218 : for (i = 1; i < l; ++i)
230 1183 : gel(C, i) = powersr(gel(Q, i), rbnd[i]);
231 35 : R = cgetg(l,t_VEC);
232 35 : r0 = real_0(prec);
233 1218 : for (i = 1; i < l; ++i)
234 : {
235 1183 : gel(R, i) = cgetg(rbnd[i]+1, t_VEC);
236 1183 : gmael(R, i, 1) = cgetr(prec);
237 1183 : affrr(gmael(C, i, 2),gmael(R, i, 1));
238 80234 : for (j = 2; j <= rbnd[i]; j++)
239 : {
240 79051 : gmael(R, i, j) = cgetr(prec);
241 79051 : affrr(r0, gmael(R, i, j));
242 : }
243 : }
244 35 : bb->baby = C; bb->giant = R;
245 35 : bb->bnd = bnd; bb->rbnd = rbnd;
246 35 : }
247 :
248 : static long
249 245 : baby_size(GEN rbnd, long Ks, long prec)
250 : {
251 245 : long i, s, m, l = lg(rbnd);
252 13349 : for (s = 0, i = 1; i < l; ++i)
253 13104 : s += rbnd[i];
254 245 : m = 2*s*prec + 3*l + s;
255 245 : if (DEBUGLEVEL)
256 0 : err_printf("ellL1: BG_add: %ld words, ellan: %ld words\n", m, Ks);
257 245 : return m;
258 : }
259 :
260 : static void
261 454972 : ellL1_add(void *E, GEN n, GEN a)
262 : {
263 454972 : pari_sp av = avma;
264 454972 : struct baby_giant *bb = (struct baby_giant*) E;
265 454972 : long j, l = lg(bb->giant);
266 2292864 : for (j = 1; j < l; j++)
267 2292864 : if (cmpii(n, gel(bb->bnd,j)) <= 0)
268 : {
269 1837892 : ulong r, q = uabsdiviu_rem(n, bb->rbnd[j], &r);
270 1837892 : GEN giant = gel(bb->giant, j), baby = gel(bb->baby, j);
271 1837892 : affrr(addrr(gel(giant, q+1), mulri(gel(baby, r+1), a)), gel(giant, q+1));
272 1837892 : set_avma(av);
273 454972 : } else break;
274 454972 : }
275 :
276 : static GEN
277 35 : vecF2_lk_bsgs(GEN E, GEN bnd, GEN rbnd, GEN Q, GEN sleh, GEN N, long prec)
278 : {
279 : struct bg_data bg;
280 : struct baby_giant bb;
281 35 : long k, L = lg(bnd)-1;
282 : GEN S;
283 35 : baby_init(&bb, Q, bnd, rbnd, prec);
284 35 : gen_BG_init(&bg, E, N, gel(bnd,1));
285 35 : gen_BG_rec((void*) &bb, ellL1_add, &bg);
286 35 : S = cgetg(L+1, t_VEC);
287 1218 : for (k = 1; k <= L; ++k)
288 : {
289 1183 : pari_sp av = avma;
290 1183 : long j, g = rbnd[k];
291 1183 : GEN giant = gmael(bb.baby, k, g+1), Sl = gmael(bb.giant, k, g);
292 80234 : for (j = g-1; j >=1; j--) Sl = addrr(mulrr(Sl, giant), gmael(bb.giant,k,j));
293 1183 : gel(S, k) = gerepileuptoleaf(av, mulrr(gel(sleh,k), Sl));
294 : }
295 35 : return S;
296 : }
297 :
298 : static long
299 13104 : _sqrt(GEN x) { pari_sp av = avma; return gc_long(av, itou(sqrtint(x))); }
300 :
301 : static GEN
302 245 : vecF(struct lcritical *C, GEN E)
303 : {
304 245 : pari_sp av = avma;
305 245 : long prec = C->cprec, Ks = itos_or_0(C->K), L = C->L, l;
306 245 : GEN N = ellQ_get_N(E), PiN;
307 245 : GEN e = mpexp(C->h), elh = powersr(e, L-1), Q, bnd, rbnd, vec;
308 :
309 245 : PiN = divrr(Pi2n(1,prec), sqrtr_abs(itor(N, prec)));
310 245 : setsigne(PiN, -1); /* - 2Pi/sqrt(N) */
311 245 : bnd = gpowers0(invr(e), L-1, C->K); /* bnd[i] = K exp(-(i-1)h) */
312 245 : rbnd = cgetg(L+1, t_VECSMALL);
313 245 : Q = cgetg(L+1, t_VEC);
314 13349 : for (l = 1; l <= L; l++)
315 : {
316 13104 : gel(bnd,l) = ceil_safe(gel(bnd,l));
317 13104 : rbnd[l] = _sqrt(gel(bnd,l)) + 1;
318 13104 : gel(Q, l) = mpexp(mulrr(PiN, gel(elh, l)));
319 : }
320 245 : if (Ks && baby_size(rbnd, Ks, prec) > (Ks>>1))
321 210 : vec = vecF2_lk(E, bnd, rbnd, Q, elh, prec);
322 : else
323 35 : vec = vecF2_lk_bsgs(E, bnd, rbnd, Q, elh, N, prec);
324 245 : return gerepileupto(av, vec);
325 : }
326 :
327 : /* Lambda function by Fourier inversion. vec is a grid, t a scalar or t_SER */
328 : static GEN
329 273 : glambda(GEN t, GEN vec, GEN h, long real, long prec)
330 : {
331 273 : GEN z, r, e = gexp(gmul(mkcomplex(gen_0,h), t), prec);
332 273 : long n = lg(vec)-1, i;
333 :
334 273 : r = real == 1? gmul2n(real_i(gel(vec, 1)), -1): gen_0;
335 273 : z = real == 1? e: gmul(powIs(3), e);
336 : /* FIXME: summing backward may be more stable */
337 15841 : for (i = 2; i <= n; i++)
338 : {
339 15568 : r = gadd(r, real_i(gmul(gel(vec,i), z)));
340 15568 : if (i < n) z = gmul(z, e);
341 : }
342 273 : return gmul(mulsr(4, h), r);
343 : }
344 :
345 : static GEN
346 245 : Lpoints(struct lcritical *C, GEN e, double tmax, long bprec)
347 : {
348 245 : double h = 0, Y = .97;
349 245 : GEN N = ellQ_get_N(e);
350 245 : param_points(N, Y, tmax, bprec, &C->cprec, &C->L, &C->K, &h);
351 245 : C->real = ellrootno_global(e);
352 245 : C->h = rtor(dbltor(h), C->cprec);
353 245 : return vecF(C, e);
354 : }
355 :
356 : static GEN
357 273 : Llambda(GEN vec, struct lcritical *C, GEN t, long prec)
358 : {
359 273 : GEN lambda = glambda(gprec_w(t, C->cprec), vec, C->h, C->real, C->cprec);
360 273 : return gprec_w(lambda, prec);
361 : }
362 :
363 : /* 2*(2*Pi)^(-s)*gamma(s)*N^(s/2); */
364 : static GEN
365 273 : ellgammafactor(GEN N, GEN s, long prec)
366 : {
367 273 : GEN c = gpow(divrr(gsqrt(N,prec), Pi2n(1,prec)), s, prec);
368 273 : return gmul(gmul2n(c,1), ggamma(s, prec));
369 : }
370 :
371 : static GEN
372 273 : ellL1_eval(GEN e, GEN vec, struct lcritical *C, GEN t, long prec)
373 : {
374 273 : GEN g = ellgammafactor(ellQ_get_N(e), gaddgs(gmul(gen_I(),t), 1), prec);
375 273 : return gdiv(Llambda(vec, C, t, prec), g);
376 : }
377 :
378 : static GEN
379 273 : ellL1_der(GEN e, GEN vec, struct lcritical *C, GEN t, long der, long prec)
380 : {
381 273 : GEN r = polcoef_i(ellL1_eval(e, vec, C, t, prec), der, 0);
382 273 : r = gmul(r,powIs(C->real == 1 ? -der: 1-der));
383 273 : return gmul(real_i(r), mpfact(der));
384 : }
385 :
386 : GEN
387 231 : ellL1_bitprec(GEN E, long r, long bitprec)
388 : {
389 231 : pari_sp av = avma;
390 : struct lcritical C;
391 231 : long prec = nbits2prec(bitprec);
392 : GEN e, vec, t;
393 231 : if (r < 0)
394 7 : pari_err_DOMAIN("ellL1", "derivative order", "<", gen_0, stoi(r));
395 224 : e = ellanal_globalred(E, NULL);
396 224 : if (r == 0 && ellrootno_global(e) < 0) { set_avma(av); return gen_0; }
397 210 : vec = Lpoints(&C, e, 0., bitprec);
398 210 : t = r ? scalarser(gen_1, 0, r): zeroser(0, 0);
399 210 : setvalser(t, 1);
400 210 : return gerepileupto(av, ellL1_der(e, vec, &C, t, r, prec));
401 : }
402 :
403 : GEN
404 0 : ellL1(GEN E, long r, long prec) { return ellL1_bitprec(E, r, prec2nbits(prec)); }
405 :
406 : GEN
407 35 : ellanalyticrank_bitprec(GEN E, GEN eps, long bitprec)
408 : {
409 35 : pari_sp av = avma, av2;
410 35 : long prec = nbits2prec(bitprec);
411 : struct lcritical C;
412 : pari_timer ti;
413 : GEN e, vec;
414 : long rk;
415 35 : if (DEBUGLEVEL) timer_start(&ti);
416 35 : if (!eps)
417 35 : eps = real2n(-bitprec/2+1, DEFAULTPREC);
418 : else
419 0 : if (typ(eps) != t_REAL) {
420 0 : eps = gtofp(eps, DEFAULTPREC);
421 0 : if (typ(eps) != t_REAL) pari_err_TYPE("ellanalyticrank", eps);
422 : }
423 35 : e = ellanal_globalred(E, NULL);
424 35 : vec = Lpoints(&C, e, 0., bitprec);
425 35 : if (DEBUGLEVEL) timer_printf(&ti, "init L");
426 35 : av2 = avma;
427 35 : for (rk = C.real>0 ? 0: 1; ; rk += 2)
428 28 : {
429 : GEN Lrk;
430 63 : GEN t = rk ? scalarser(gen_1, 0, rk): zeroser(0, 0);
431 63 : setvalser(t, 1);
432 63 : Lrk = ellL1_der(e, vec, &C, t, rk, prec);
433 63 : if (DEBUGLEVEL) timer_printf(&ti, "L^(%ld)=%Ps", rk, Lrk);
434 63 : if (abscmprr(Lrk, eps) > 0)
435 35 : return gerepilecopy(av, mkvec2(stoi(rk), Lrk));
436 28 : set_avma(av2);
437 : }
438 : }
439 :
440 : GEN
441 0 : ellanalyticrank(GEN E, GEN eps, long prec)
442 : {
443 0 : return ellanalyticrank_bitprec(E, eps, prec2nbits(prec));
444 : }
445 :
446 : /* Heegner point computation
447 :
448 : This section is a C version by Bill Allombert of a GP script by
449 : Christophe Delaunay which was based on a GP script by John Cremona.
450 : Reference: Henri Cohen's book GTM 239.
451 : */
452 :
453 : static void
454 0 : heegner_L1_bg(void*E, GEN n, GEN a)
455 : {
456 0 : struct baby_giant *bb = (struct baby_giant*) E;
457 0 : long j, l = lg(bb->giant);
458 0 : for (j = 1; j < l; j++)
459 0 : if (cmpii(n, gel(bb->bnd,j)) <= 0)
460 : {
461 0 : ulong r, q = uabsdiviu_rem(n, bb->rbnd[j], &r);
462 0 : GEN giant = gel(bb->giant, j), baby = gel(bb->baby, j);
463 0 : gaffect(gadd(gel(giant, q+1), gdiv(gmul(gel(baby, r+1), a), n)), gel(giant, q+1));
464 : }
465 0 : }
466 :
467 : static void
468 2102611 : heegner_L1(void*E, GEN n, GEN a)
469 : {
470 2102611 : struct baby_giant *bb = (struct baby_giant*) E;
471 2102611 : long j, l = lg(bb->giant);
472 12033980 : for (j = 1; j < l; j++)
473 9931369 : if (cmpii(n, gel(bb->bnd,j)) <= 0)
474 : {
475 8349474 : ulong r, q = uabsdiviu_rem(n, bb->rbnd[j], &r);
476 8349474 : GEN giant = gel(bb->giant, j), baby = gel(bb->baby, j);
477 8349474 : GEN ex = mulreal(gel(baby, r+1), gel(giant, q+1));
478 8349474 : affrr(addrr(gel(bb->sum, j), divri(mulri(ex, a), n)), gel(bb->sum, j));
479 : }
480 2102611 : }
481 : /* Return C, C[i][j] = Q[j]^i, i = 1..nb */
482 : static void
483 0 : baby_init2(struct baby_giant *bb, GEN Q, GEN bnd, GEN rbnd, long prec)
484 : {
485 0 : long i, j, l = lg(Q);
486 : GEN R, C, r0;
487 0 : C = cgetg(l,t_VEC);
488 0 : for (i = 1; i < l; ++i)
489 0 : gel(C, i) = gpowers(gel(Q, i), rbnd[i]);
490 0 : R = cgetg(l,t_VEC);
491 0 : r0 = mkcomplex(real_0(prec),real_0(prec));
492 0 : for (i = 1; i < l; ++i)
493 : {
494 0 : gel(R, i) = cgetg(rbnd[i]+1, t_VEC);
495 0 : gmael(R, i, 1) = cgetc(prec);
496 0 : gaffect(gmael(C, i, 2),gmael(R, i, 1));
497 0 : for (j = 2; j <= rbnd[i]; j++)
498 : {
499 0 : gmael(R, i, j) = cgetc(prec);
500 0 : gaffect(r0, gmael(R, i, j));
501 : }
502 : }
503 0 : bb->baby = C; bb->giant = R;
504 0 : bb->bnd = bnd; bb->rbnd = rbnd;
505 0 : }
506 :
507 : /* Return C, C[i][j] = Q[j]^i, i = 1..nb */
508 : static void
509 42 : baby_init3(struct baby_giant *bb, GEN Q, GEN bnd, GEN rbnd, long prec)
510 : {
511 42 : long i, l = lg(Q);
512 : GEN R, C, S;
513 42 : C = cgetg(l,t_VEC);
514 189 : for (i = 1; i < l; ++i)
515 147 : gel(C, i) = gpowers(gel(Q, i), rbnd[i]);
516 42 : R = cgetg(l,t_VEC);
517 189 : for (i = 1; i < l; ++i)
518 147 : gel(R, i) = gpowers(gmael(C, i, 1+rbnd[i]), rbnd[i]);
519 42 : S = cgetg(l,t_VEC);
520 189 : for (i = 1; i < l; ++i)
521 : {
522 147 : gel(S, i) = cgetr(prec);
523 147 : affrr(real_i(gmael(C, i, 2)), gel(S, i));
524 : }
525 42 : bb->baby = C; bb->giant = R; bb->sum = S;
526 42 : bb->bnd = bnd; bb->rbnd = rbnd;
527 42 : }
528 :
529 : /* ymin a t_REAL */
530 : static GEN
531 42 : heegner_psi(GEN E, GEN N, GEN points, long bitprec)
532 : {
533 42 : pari_sp av = avma, av2;
534 : struct baby_giant bb;
535 : struct bg_data bg;
536 42 : long l, k, L = lg(points)-1, prec = nbits2prec(bitprec)+EXTRAPREC64;
537 42 : GEN Q, pi2 = Pi2n(1, prec), bnd, rbnd, bndmax;
538 42 : GEN B = divrr(mulur(bitprec,mplog2(DEFAULTPREC)), pi2);
539 :
540 42 : rbnd = cgetg(L+1, t_VECSMALL); av2 = avma;
541 42 : bnd = cgetg(L+1, t_VEC);
542 42 : Q = cgetg(L+1, t_VEC);
543 189 : for (l = 1; l <= L; ++l)
544 : {
545 147 : gel(bnd,l) = ceil_safe(divrr(B,imag_i(gel(points, l))));
546 147 : rbnd[l] = itou(sqrtint(gel(bnd,l)))+1;
547 147 : gel(Q, l) = expIxy(pi2, gel(points, l), prec);
548 : }
549 42 : gerepileall(av2, 2, &bnd, &Q);
550 42 : bndmax = gel(bnd,vecindexmax(bnd));
551 42 : gen_BG_init(&bg, E, N, bndmax);
552 42 : if (bitprec >= 1900)
553 : {
554 0 : GEN S = cgetg(L+1, t_VEC);
555 0 : baby_init2(&bb, Q, bnd, rbnd, prec);
556 0 : gen_BG_rec((void*)&bb, heegner_L1_bg, &bg);
557 0 : for (k = 1; k <= L; ++k)
558 : {
559 0 : pari_sp av2 = avma;
560 0 : long j, g = rbnd[k];
561 0 : GEN giant = gmael(bb.baby, k, g+1), Sl = real_0(prec);
562 0 : for (j = g; j >=1; j--) Sl = gadd(gmul(Sl, giant), gmael(bb.giant,k,j));
563 0 : gel(S, k) = gerepileupto(av2, real_i(Sl));
564 : }
565 0 : return gerepileupto(av, S);
566 : }
567 : else
568 : {
569 42 : baby_init3(&bb, Q, bnd, rbnd, prec);
570 42 : gen_BG_rec((void*)&bb, heegner_L1, &bg);
571 42 : return gerepilecopy(av, bb.sum);
572 : }
573 : }
574 :
575 : /*Returns lambda_bad list for one prime p, nv = localred(E, p) */
576 : static GEN
577 91 : lambda1(GEN E, GEN nv, GEN p, long prec)
578 : {
579 : pari_sp av;
580 : GEN res, lp;
581 91 : long kod = itos(gel(nv, 2));
582 91 : if (kod==2 || kod ==-2) return cgetg(1,t_VEC);
583 91 : av = avma; lp = glog(p, prec);
584 91 : if (kod > 4)
585 : {
586 14 : long n = Z_pval(ell_get_disc(E), p);
587 14 : long j, m = kod - 4, nl = 1 + (m >> 1L);
588 14 : res = cgetg(nl, t_VEC);
589 35 : for (j = 1; j < nl; j++)
590 21 : gel(res, j) = gmul(lp, gsubgs(gdivgu(sqru(j), n), j)); /* j^2/n - j */
591 : }
592 77 : else if (kod < -4)
593 14 : res = mkvec2(negr(lp), shiftr(mulrs(lp, kod), -2));
594 : else
595 : {
596 63 : const long lam[] = {8,9,0,6,0,0,0,3,4};
597 63 : long m = -lam[kod+4];
598 63 : res = mkvec(divru(mulrs(lp, m), 6));
599 : }
600 91 : return gerepilecopy(av, res);
601 : }
602 :
603 : static GEN
604 42 : lambdalist(GEN E, long prec)
605 : {
606 42 : pari_sp ltop = avma;
607 42 : GEN glob = ellglobalred(E), plist = gmael(glob,4,1), L = gel(glob,5);
608 42 : GEN res, v, D = ell_get_disc(E);
609 42 : long i, j, k, l, m, n, np = lg(plist), lr = 1;
610 42 : v = cgetg(np, t_VEC);
611 147 : for (j = 1, i = 1 ; j < np; ++j)
612 : {
613 105 : GEN p = gel(plist, j);
614 105 : if (dvdii(D, sqri(p)))
615 : {
616 91 : GEN la = lambda1(E, gel(L,j), p, prec);
617 91 : gel(v, i++) = la;
618 91 : lr *= lg(la);
619 : }
620 : }
621 42 : np = i;
622 42 : res = cgetg(lr+1, t_VEC);
623 42 : gel(res, 1) = gen_0; n = 1; m = 1;
624 133 : for (j = 1; j < np; ++j)
625 : {
626 91 : GEN w = gel(v, j);
627 91 : long lw = lg(w);
628 308 : for (k = 1; k <= n; k++)
629 : {
630 217 : GEN t = gel(res, k);
631 455 : for (l = 1, m = n; l < lw; l++, m+=n)
632 238 : gel(res, k + m) = mpadd(t, gel(w, l));
633 : }
634 91 : n = m;
635 : }
636 42 : return gerepilecopy(ltop, res);
637 : }
638 :
639 : /* P a t_INT or t_FRAC, return its logarithmic height */
640 : static GEN
641 98 : heightQ(GEN P, long prec)
642 : {
643 : long s;
644 98 : if (typ(P) == t_FRAC)
645 : {
646 56 : GEN a = gel(P,1), b = gel(P,2);
647 56 : P = abscmpii(a,b) > 0 ? a: b;
648 : }
649 98 : s = signe(P);
650 98 : if (!s) return real_0(prec);
651 84 : if (s < 0) P = negi(P);
652 84 : return glog(P, prec);
653 : }
654 :
655 : /* t a t_INT or t_FRAC, returns max(1, log |t|), returns a t_REAL */
656 : static GEN
657 147 : logplusQ(GEN t, long prec)
658 : {
659 147 : if (typ(t) == t_INT)
660 : {
661 42 : if (!signe(t)) return real_1(prec);
662 28 : if (signe(t) < 0) t = negi(t);
663 : }
664 : else
665 : {
666 105 : GEN a = gel(t,1), b = gel(t,2);
667 105 : if (abscmpii(a, b) < 0) return real_1(prec);
668 56 : if (signe(a) < 0) t = gneg(t);
669 : }
670 84 : return glog(t, prec);
671 : }
672 :
673 : /* See GTM239, p532, Th 8.1.18
674 : * Return M such that h_naive <= M */
675 : GEN
676 98 : hnaive_max(GEN ell, GEN ht)
677 : {
678 98 : const long prec = LOWDEFAULTPREC; /* minimal accuracy */
679 98 : GEN b2 = ell_get_b2(ell), j = ell_get_j(ell);
680 98 : GEN logd = glog(absi_shallow(ell_get_disc(ell)), prec);
681 98 : GEN logj = logplusQ(j, prec);
682 98 : GEN hj = heightQ(j, prec);
683 49 : GEN logb2p = signe(b2)? addrr(logplusQ(gdivgu(b2, 12),prec), mplog2(prec))
684 98 : : real_1(prec);
685 98 : GEN mu = addrr(divru(addrr(logd, logj),6), logb2p);
686 98 : return addrs(addrr(addrr(ht, divru(hj,12)), mu), 2);
687 : }
688 :
689 : static GEN
690 147 : qfb_root(GEN Q, GEN vDi)
691 : {
692 147 : GEN a2 = shifti(gel(Q, 1),1), b = gel(Q, 2);
693 147 : return mkcomplex(gdiv(negi(b),a2),divri(vDi,a2));
694 : }
695 :
696 : static GEN
697 24668 : qimag2(GEN Q)
698 : {
699 24668 : pari_sp av = avma;
700 24668 : GEN z = gdiv(negi(qfb_disc(Q)), shifti(sqri(gel(Q, 1)),2));
701 24668 : return gerepileupto(av, z);
702 : }
703 :
704 : /***************************************************/
705 : /*Routines for increasing the imaginary parts using*/
706 : /*Atkin-Lehner operators */
707 : /***************************************************/
708 :
709 : static GEN
710 24668 : qfb_mult(GEN Q, GEN a, GEN b, GEN c, GEN d)
711 : {
712 24668 : GEN A = gel(Q, 1) , B = gel(Q, 2), C = gel(Q, 3), D = qfb_disc(Q);
713 24668 : GEN a2 = sqri(a), b2 = sqri(b), c2 = sqri(c), d2 = sqri(d);
714 24668 : GEN ad = mulii(d, a), bc = mulii(b, c), e = subii(ad, bc);
715 24668 : GEN W1 = addii(addii(mulii(a2, A), mulii(mulii(c, a), B)), mulii(c2, C));
716 24668 : GEN W3 = addii(addii(mulii(b2, A), mulii(mulii(d, b), B)), mulii(d2, C));
717 24668 : GEN W2 = addii(addii(mulii(mulii(shifti(b,1), a), A),
718 : mulii(addii(ad, bc), B)),
719 : mulii(mulii(shifti(d,1), c), C));
720 24668 : if (!equali1(e)) {
721 22190 : W1 = diviiexact(W1,e);
722 22190 : W2 = diviiexact(W2,e);
723 22190 : W3 = diviiexact(W3,e);
724 : }
725 24668 : return mkqfb(W1, W2, W3, D);
726 : }
727 :
728 : #ifdef DEBUG
729 : static void
730 : best_point_old(GEN Q, GEN NQ, GEN f, GEN *u, GEN *v)
731 : {
732 : long n, k;
733 : GEN U, c, d, A = gel(f,1), B = gel(f,2), C = gel(f,3), D = qfb_disc(f);
734 : GEN q = mkqfb(mulii(NQ, C), negi(B), diviiexact(A, NQ), D);
735 : redimagsl2(q, &U);
736 : *u = c = gcoeff(U, 1, 1);
737 : *v = d = gcoeff(U, 2, 1);
738 : if (equali1(gcdii(mulii(*u, NQ), mulii(*v, Q)))) return;
739 : for (n = 1;; n++)
740 : {
741 : for (k = -n; k <= n; k++)
742 : {
743 : *u = addis(c, k); *v = addiu(d, n);
744 : if (equali1(gcdii(mulii(*u, NQ), mulii(*v, Q)))) return;
745 : *v = subiu(d, n);
746 : if (equali1(gcdii(mulii(*u, NQ), mulii(*v, Q)))) return;
747 : *u = addiu(c, n); *v = addis(d, k);
748 : if (equali1(gcdii(mulii(*u, NQ), mulii(*v, Q)))) return;
749 : *u = subiu(c, n);
750 : if (equali1(gcdii(mulii(*u, NQ), mulii(*v, Q)))) return;
751 : }
752 : }
753 : }
754 : /* q(x,y) = ax^2 + bxy + cy^2 */
755 : static GEN
756 : qfb_eval(GEN q, GEN x, GEN y)
757 : {
758 : GEN a = gel(q,1), b = gel(q,2), c = gel(q,3);
759 : GEN x2 = sqri(x), y2 = sqri(y), xy = mulii(x,y);
760 : return addii(addii(mulii(a, x2), mulii(b,xy)), mulii(c, y2));
761 : }
762 : #endif
763 :
764 : static long
765 6580 : nexti(long i) { return i>0 ? -i : 1-i; }
766 :
767 : /* q0 + i q1 + i^2 q2 */
768 : static GEN
769 12313 : qfmin_eval(GEN q0, GEN q1, GEN q2, long i)
770 12313 : { return addii(mulis(addii(mulis(q2, i), q1), i), q0); }
771 :
772 : /* assume a > 0, return gcd(a,b,c) */
773 : static ulong
774 16436 : gcduii(ulong a, GEN b, GEN c)
775 : {
776 16436 : a = ugcdiu(b, a);
777 16436 : return a == 1? 1: ugcdiu(c, a);
778 : }
779 :
780 : static void
781 24668 : best_point(GEN Q, GEN NQ, GEN f, GEN *pu, GEN *pv)
782 : {
783 24668 : GEN a = mulii(NQ, gel(f,3)), b = negi(gel(f,2)), c = diviiexact(gel(f,1), NQ);
784 24668 : GEN D = qfb_disc(f);
785 24668 : GEN U, qr = redimagsl2(mkqfb(a, b, c, D), &U);
786 24668 : GEN A = gel(qr,1), B = gel(qr,2), A2 = shifti(A,1), AA4 = sqri(A2);
787 : GEN V, best;
788 : long y;
789 :
790 24668 : D = absi_shallow(D);
791 : /* 4A qr(x,y) = (2A x + By)^2 + D y^2
792 : * Write x = x0(y) + i, where x0 is an integer minimum
793 : * (the smallest in case of tie) of x-> qr(x,y), for given y.
794 : * 4A qr(x,y) = ((2A x0 + By)^2 + Dy^2) + 4A i (2A x0 + By) + 4A^2 i^2
795 : * = q0(y) + q1(y) i + q2 i^2
796 : * Loop through (x,y), y>0 by (roughly) increasing values of qr(x,y) */
797 :
798 : /* We must test whether [X,Y]~ := U * [x,y]~ satisfy (X NQ, Y Q) = 1
799 : * This is equivalent to (X,Y) = 1 (note that (X,Y) = (x,y)), and
800 : * (X, Q) = (Y, NQ) = 1.
801 : * We have U * [x0+i, y]~ = U * [x0,y]~ + i U[,1] =: V0 + i U[,1] */
802 :
803 : /* try [1,0]~ = first minimum */
804 24668 : V = gel(U,1); /* U *[1,0]~ */
805 24668 : *pu = gel(V,1);
806 24668 : *pv = gel(V,2);
807 30947 : if (is_pm1(gcdii(*pu, Q)) && is_pm1(gcdii(*pv, NQ))) return;
808 :
809 : /* try [0,1]~ = second minimum */
810 11935 : V = gel(U,2); /* U *[0,1]~ */
811 11935 : *pu = gel(V,1);
812 11935 : *pv = gel(V,2);
813 11935 : if (is_pm1(gcdii(*pu, Q)) && is_pm1(gcdii(*pv, NQ))) return;
814 :
815 : /* (X,Y) = (1, \pm1) always works. Try to do better now */
816 5656 : best = subii(addii(a, c), absi_shallow(b));
817 5656 : *pu = gen_1;
818 5656 : *pv = signe(b) < 0? gen_1: gen_m1;
819 :
820 5656 : for (y = 1;; y++)
821 9128 : {
822 : GEN Dy2, r, By, x0, q0, q1, V0;
823 : long i;
824 14784 : if (y > 1)
825 : {
826 10703 : if (gcduii(y, gcoeff(U,1,1), Q) != 1) continue;
827 7308 : if (gcduii(y, gcoeff(U,2,1), NQ) != 1) continue;
828 : }
829 11396 : Dy2 = mulii(D, sqru(y));
830 11396 : if (cmpii(Dy2, best) >= 0) break; /* we won't improve. STOP */
831 5740 : By = muliu(B,y), x0 = truedvmdii(negi(By), A2, &r);
832 5740 : if (cmpii(r, A) >= 0) { x0 = subiu(x0,1); r = subii(r, A2); }
833 : /* (2A x + By)^2 + Dy^2, minimal at x = x0. Assume A2 > 0 */
834 : /* r = 2A x0 + By */
835 5740 : q0 = addii(sqri(r), Dy2); /* minimal value for this y, at x0 */
836 5740 : if (cmpii(q0, best) >= 0) continue; /* we won't improve for this y */
837 5733 : q1 = shifti(mulii(A2, r), 1);
838 :
839 5733 : V0 = ZM_ZC_mul(U, mkcol2(x0, utoipos(y)));
840 12313 : for (i = 0;; i = nexti(i))
841 6580 : {
842 12313 : pari_sp av2 = avma;
843 12313 : GEN x, N = qfmin_eval(q0, q1, AA4, i);
844 12313 : if (cmpii(N , best) >= 0) break;
845 12271 : x = addis(x0, i);
846 12271 : if (ugcdiu(x, y) == 1)
847 : {
848 : GEN u, v;
849 12229 : V = ZC_add(V0, ZC_z_mul(gel(U,1), i)); /* [X, Y] */
850 12229 : u = gel(V,1);
851 12229 : v = gel(V,2);
852 12229 : if (is_pm1(gcdii(u, Q)) && is_pm1(gcdii(v, NQ)))
853 : {
854 5691 : *pu = u;
855 5691 : *pv = v;
856 5691 : best = N; break;
857 : }
858 : }
859 6580 : set_avma(av2);
860 : }
861 : }
862 : #ifdef DEBUG
863 : {
864 : GEN oldu, oldv, F = mkqfb(a, b, c, qfb_disc(f));
865 : best_point_old(Q, NQ, f, &oldu, &oldv);
866 : if (!equalii(oldu, *pu) || !equalii(oldv, *pv))
867 : {
868 : if (!equali1(gcdii(mulii(*pu, NQ), mulii(*pv, Q))))
869 : pari_err_BUG("best_point (gcd)");
870 : if (cmpii(qfb_eval(F, *pu,*pv), qfb_eval(F, oldu, oldv)) > 0)
871 : {
872 : pari_warn(warner, "%Ps,%Ps,%Ps, %Ps > %Ps",
873 : Q,NQ,f, mkvec2(*pu,*pv), mkvec2(oldu,oldv));
874 : pari_err_BUG("best_point (too large)");
875 : }
876 : }
877 : }
878 : #endif
879 : }
880 :
881 : static GEN
882 24668 : best_lift(GEN Q, GEN NQ, GEN f)
883 : {
884 : GEN a, b, c, d, dQ, cNQ;
885 24668 : best_point(Q, NQ, f, &c, &d);
886 24668 : dQ = mulii(d, Q); cNQ = mulii(NQ, c);
887 24668 : (void)bezout(dQ, cNQ, &a, &b);
888 24668 : return qfb_mult(f, dQ, b, mulii(negi(Q),cNQ), mulii(a,Q));
889 : }
890 :
891 : static GEN
892 2478 : lift_points(GEN listQ, GEN f, GEN *pt, GEN *pQ)
893 : {
894 2478 : pari_sp av = avma;
895 2478 : GEN yf = gen_0, tf = NULL, Qf = NULL;
896 2478 : long k, l = lg(listQ);
897 27146 : for (k = 1; k < l; ++k)
898 : {
899 24668 : GEN c = gel(listQ, k), Q = gel(c,1), NQ = gel(c,2);
900 24668 : GEN t = best_lift(Q, NQ, f), y = qimag2(t);
901 24668 : if (gcmp(y, yf) > 0) { yf = y; Qf = Q; tf = t; }
902 : }
903 2478 : *pt = tf; *pQ = Qf; return gc_all(av, 3, &yf, pt, pQ);
904 : }
905 :
906 : /***************************/
907 : /* Twists */
908 : /***************************/
909 :
910 : static GEN
911 56 : ltwist1(GEN E, GEN d, long bitprec)
912 : {
913 56 : pari_sp av = avma;
914 56 : GEN Ed = elltwist(E, d), z = ellL1_bitprec(Ed, 0, bitprec);
915 56 : obj_free(Ed); return gerepileuptoleaf(av, z);
916 : }
917 :
918 : /* Return O_re*c(E)/(4*O_vol*|E_t|^2) */
919 :
920 : static GEN
921 42 : heegner_indexmult(GEN om, long t, GEN tam, long prec)
922 : {
923 42 : pari_sp av = avma;
924 42 : GEN Ovr = gabs(imag_i(gel(om, 2)), prec); /* O_vol/O_re, t_REAL */
925 42 : return gerepileupto(av, divru(divir(tam, Ovr), 4*t*t));
926 : }
927 :
928 : /* omega(gcd(D, N)), given faN = factor(N) */
929 : static long
930 56 : omega_N_D(GEN faN, ulong D)
931 : {
932 56 : GEN P = gel(faN, 1);
933 56 : long i, l = lg(P), w = 0;
934 196 : for (i = 1; i < l; i++)
935 140 : if (dvdui(D, gel(P,i))) w++;
936 56 : return w;
937 : }
938 :
939 : static GEN
940 56 : heegner_indexmultD(GEN faN, GEN a, long D, GEN sqrtD)
941 : {
942 56 : pari_sp av = avma;
943 : GEN c;
944 : long w;
945 56 : switch(D)
946 : {
947 0 : case -3: w = 9; break;
948 0 : case -4: w = 4; break;
949 56 : default: w = 1;
950 : }
951 56 : c = shifti(stoi(w), omega_N_D(faN,-D)); /* (w(D)/2)^2 * 2^omega(gcd(D,N)) */
952 56 : return gerepileupto(av, mulri(mulrr(a, sqrtD), c));
953 : }
954 :
955 : static GEN
956 399 : nf_to_basis(GEN nf, GEN x)
957 : {
958 399 : x = nf_to_scalar_or_basis(nf, x);
959 399 : if (typ(x)!=t_COL)
960 287 : x = scalarcol(x, nf_get_degree(nf));
961 399 : return x;
962 : }
963 :
964 : static GEN
965 196 : etnf_to_basis(GEN et, GEN x)
966 : {
967 196 : long i, l = lg(et);
968 196 : GEN V = cgetg(l, t_VEC);
969 595 : for (i = 1; i < l; i++)
970 399 : gel(V,i) = nf_to_basis(gel(et,i), x);
971 196 : return shallowconcat1(V);
972 : }
973 :
974 : static GEN
975 140 : etnf_get_M(GEN et)
976 : {
977 140 : long i, l = lg(et);
978 140 : GEN V = cgetg(l, t_VEC);
979 448 : for (i = 1; i < l; i++)
980 308 : gel(V,i)=nf_get_M(gel(et,i));
981 140 : return shallowmatconcat(diagonal(V));
982 : }
983 :
984 : static long
985 49 : etnf_get_varn(GEN et)
986 : {
987 49 : return nf_get_varn(gel(et,1));
988 : }
989 :
990 : static GEN
991 98 : heegner_descent_try_point(GEN nfA, GEN z, GEN den, long prec)
992 : {
993 98 : pari_sp av = avma;
994 98 : GEN etal = gel(nfA,1), A = gel(nfA,2), cb = gel(nfA,3);
995 98 : GEN al = gel(nfA,4), th = gel(nfA, 5);
996 98 : GEN et = gel(etal,1), zk = gel(etal, 2), T = gel(etal,3);
997 98 : GEN M = etnf_get_M(et);
998 98 : long i, j, n = lg(th)-1, l = lg(al);
999 98 : GEN u2 = gsqr(gel(cb,1)), r = gel(cb,2);
1000 98 : GEN zz = gdiv(gsub(z,r), u2);
1001 98 : GEN be = cgetg(n+1, t_COL);
1002 161 : for (j = 1; j < l; j++)
1003 : {
1004 105 : GEN aj = gel(al, j), Aj = gel(A,j);
1005 392 : for (i = 1; i <= n; i++)
1006 287 : gel(be,i) = gsqrt(gmul(gsub(zz, gel(th,i)), gel(aj,i)), prec);
1007 357 : for (i = 0; i <= (1<<(n-1))-1; i++)
1008 : {
1009 : long eps;
1010 294 : GEN s = gmul(den, RgM_solve_realimag(M, be));
1011 294 : GEN S = grndtoi(s, &eps), V, S2;
1012 294 : gel(be,1+odd(i)) = gneg(gel(be,1+odd(i)));
1013 294 : if (eps > -7) continue;
1014 42 : S2 = QXQ_sqr(RgV_RgC_mul(zk, S), T);
1015 42 : V = gdiv(QXQ_mul(S2, Aj, T), sqri(den));
1016 42 : if (typ(V) != t_POL || degpol(V) != 1) continue;
1017 42 : if (gequalm1(gel(V,3)))
1018 42 : return gerepileupto(av,gadd(gmul(gel(V,2),u2),r));
1019 : }
1020 : }
1021 56 : return gc_NULL(av);
1022 : }
1023 :
1024 : static GEN
1025 1785 : heegner_try_point(GEN E, GEN nfA, GEN lambdas, GEN ht, GEN z, long prec)
1026 : {
1027 1785 : long l = lg(lambdas);
1028 : long i, eps;
1029 1785 : GEN P = real_i(pointell(E, z, prec)), x = gel(P,1);
1030 1785 : GEN rh = subrr(ht, shiftr(ellheightoo(E, P, prec),1));
1031 26572 : for (i = 1; i < l; ++i)
1032 : {
1033 24829 : GEN logd = shiftr(gsub(rh, gel(lambdas, i)), -1);
1034 24829 : GEN d, approxd = gexp(logd, prec);
1035 24829 : d = grndtoi(approxd, &eps);
1036 24829 : if (signe(d) > 0 && eps<-10)
1037 : {
1038 : GEN X, ylist;
1039 98 : if (DEBUGLEVEL > 2)
1040 0 : err_printf("\nTrying lambda number %ld, logd=%Ps, approxd=%Ps\n", i, logd, approxd);
1041 98 : X = heegner_descent_try_point(nfA, x, d, prec);
1042 98 : if (X)
1043 : {
1044 42 : ylist = ellordinate(E, X, prec);
1045 42 : if (lg(ylist) > 1)
1046 : {
1047 42 : GEN P = mkvec2(X, gel(ylist, 1));
1048 42 : GEN hp = ellheight(E,P,prec);
1049 42 : if (signe(hp) && cmprr(hp, shiftr(ht,1)) < 0 && cmprr(hp, shiftr(ht,-1)) > 0)
1050 42 : return P;
1051 0 : if (DEBUGLEVEL)
1052 0 : err_printf("found non-Heegner point %Ps\n", P);
1053 : }
1054 : }
1055 : }
1056 : }
1057 1743 : return NULL;
1058 : }
1059 :
1060 : static GEN
1061 42 : heegner_find_point(GEN e, GEN nfA, GEN om, GEN ht, GEN z1, long k, long prec)
1062 : {
1063 42 : GEN lambdas = lambdalist(e, prec);
1064 42 : pari_sp av = avma;
1065 : long m;
1066 42 : GEN Ore = gel(om, 1), Oim = gel(om, 2);
1067 42 : if (DEBUGLEVEL)
1068 0 : err_printf("%ld*%ld multipliers to test: ",k,lg(lambdas)-1);
1069 966 : for (m = 0; m < k; m++)
1070 : {
1071 966 : GEN P, z2 = divru(addrr(z1, mulsr(m, Ore)), k);
1072 966 : if (DEBUGLEVEL > 2)
1073 0 : err_printf("%ld ",m);
1074 966 : P = heegner_try_point(e, nfA, lambdas, ht, z2, prec);
1075 966 : if (P) return P;
1076 931 : if (signe(ell_get_disc(e)) > 0)
1077 : {
1078 819 : z2 = gadd(z2, gmul2n(Oim, -1));
1079 819 : P = heegner_try_point(e, nfA, lambdas, ht, z2, prec);
1080 819 : if (P) return P;
1081 : }
1082 924 : set_avma(av);
1083 : }
1084 0 : pari_err_BUG("ellheegner, point not found");
1085 : return NULL; /* LCOV_EXCL_LINE */
1086 : }
1087 :
1088 : /* N > 1, fa = factor(N), return factor(4*N) */
1089 : static GEN
1090 42 : fa_shift2(GEN fa)
1091 : {
1092 42 : GEN P = gel(fa,1), E = gel(fa,2);
1093 42 : if (absequaliu(gcoeff(fa,1,1), 2))
1094 : {
1095 21 : E = shallowcopy(E);
1096 21 : gel(E,1) = addiu(gel(E,1), 2);
1097 : }
1098 : else
1099 : {
1100 21 : P = shallowconcat(gen_2, P);
1101 21 : E = shallowconcat(gen_2, E);
1102 : }
1103 42 : return mkmat2(P, E);
1104 : }
1105 :
1106 : /* P = prime divisors of N(E). Return the product of primes p in P, a_p != -1
1107 : * HACK: restrict to small primes since large ones won't divide our C-long
1108 : * discriminants */
1109 : static GEN
1110 42 : get_bad(GEN E, GEN P)
1111 : {
1112 42 : long k, l = lg(P), ibad = 1;
1113 42 : GEN B = cgetg(l, t_VECSMALL);
1114 147 : for (k = 1; k < l; k++)
1115 : {
1116 105 : GEN p = gel(P,k);
1117 105 : long pp = itos_or_0(p);
1118 105 : if (!pp) break;
1119 105 : if (! equalim1(ellap(E,p))) B[ibad++] = pp;
1120 : }
1121 42 : setlg(B, ibad); return ibad == 1? NULL: zv_prod_Z(B);
1122 : }
1123 :
1124 : /* list of pairs [Q,N/Q], where Q | N and gcd(Q,N/Q) = 1 */
1125 : static GEN
1126 42 : find_div(GEN N, GEN faN)
1127 : {
1128 42 : GEN listQ = divisors(faN);
1129 42 : long j, k, l = lg(listQ);
1130 :
1131 42 : gel(listQ, 1) = mkvec2(gen_1, N);
1132 1624 : for (j = k = 2; k < l; ++k)
1133 : {
1134 1582 : GEN Q = gel(listQ, k), NQ = diviiexact(N, Q);
1135 1582 : if (is_pm1(gcdii(Q,NQ))) gel(listQ, j++) = mkvec2(Q,NQ);
1136 : }
1137 42 : setlg(listQ, j); return listQ;
1138 : }
1139 :
1140 : static long
1141 8652 : testDisc(GEN bad, long d) { return !bad || ugcdiu(bad, -d) == 1; }
1142 : /* bad = product of bad primes. Return the NDISC largest fundamental
1143 : * discriminants D < d such that (D,bad) = 1 and d is a square mod 4N */
1144 : static GEN
1145 42 : listDisc(GEN fa4N, GEN bad, long d, long ndisc)
1146 : {
1147 42 : GEN v = cgetg(ndisc+1, t_VECSMALL);
1148 42 : pari_sp av = avma;
1149 42 : long j = 1;
1150 : for(;;)
1151 : {
1152 8652 : d -= odd(d)? 1: 3;
1153 8652 : if (testDisc(bad,d) && unegisfundamental(-d) && Zn_issquare(stoi(d), fa4N))
1154 : {
1155 420 : v[j++] = d;
1156 420 : if (j > ndisc) break;
1157 : }
1158 8610 : set_avma(av);
1159 : }
1160 42 : set_avma(av); return v;
1161 : }
1162 : /* L = vector of [q1,q2] or [q1,q2,q2']
1163 : * cd = (b^2 - D)/(4N) */
1164 : static void
1165 166880 : listfill(GEN N, GEN b, GEN c, GEN d, GEN D, GEN L, long *s)
1166 : {
1167 166880 : long k, l = lg(L);
1168 166880 : GEN add, frm2, a = mulii(d, N), V = mkqfb(a,b,c,D), frm = qfbred_i(V);
1169 600089 : for (k = 1; k < l; ++k)
1170 : { /* Lk = [v,frm] or [v,frm,frm2] */
1171 597611 : GEN Lk = gel(L,k);
1172 : long i;
1173 1509515 : for (i = 2; i < lg(Lk); i++) /* 1 or 2 elements */
1174 1076306 : if (gequal(frm, gel(Lk,i)))
1175 : {
1176 164402 : GEN v = gel(Lk, 1);
1177 164402 : if (cmpii(a, gel(v,1)) < 0) gel(Lk,1) = V;
1178 164402 : return;
1179 : }
1180 : }
1181 2478 : frm2 = qfbred_i(mkqfb(d, negi(b), mulii(c,N), D));
1182 2478 : add = gequal(frm, frm2)? mkvec2(V,frm): mkvec3(V,frm,frm2);
1183 2478 : vectrunc_append(L, add);
1184 2478 : *s += lg(add) - 2;
1185 : }
1186 : /* faN4 = factor(4*N) */
1187 : static GEN
1188 420 : listheegner(GEN N, GEN faN4, GEN listQ, GEN D)
1189 : {
1190 420 : pari_sp av = avma;
1191 420 : const long kmin = 30;
1192 420 : long h = itos(quadclassno(D));
1193 420 : GEN ymin, b = Zn_sqrt(D, faN4), L = vectrunc_init(h+1);
1194 420 : long l, k, s = 0;
1195 13020 : for (k = 0; k < kmin || s < h; k++)
1196 : {
1197 12600 : GEN bk = addii(b, mulsi(2*k, N));
1198 12600 : GEN C = diviiexact(shifti(subii(sqri(bk), D), -2), N);
1199 12600 : GEN div = divisors(C);
1200 12600 : long i, l = lg(div);
1201 179480 : for (i = 1; i < l; i++)
1202 : {
1203 166880 : GEN d = gel(div, i), c = gel(div, l-i); /* cd = C */
1204 166880 : listfill(N, bk, c, d, D, L, &s);
1205 : }
1206 : }
1207 420 : l = lg(L); ymin = NULL;
1208 2898 : for (k = 1; k < l; k++)
1209 : {
1210 2478 : GEN t, Q, Lk = gel(L,k), f = gel(Lk,1);
1211 2478 : GEN y = lift_points(listQ, f, &t, &Q);
1212 2478 : gel(L, k) = mkvec3(t, stoi(lg(Lk) - 2), Q);
1213 2478 : if (!ymin || gcmp(y, ymin) < 0) ymin = y;
1214 : }
1215 420 : if (DEBUGLEVEL > 1)
1216 0 : err_printf("Disc %Ps : N*ymin = %Pg\n", D,
1217 : gmul(gsqrt(ymin, DEFAULTPREC),N));
1218 420 : return gerepilecopy(av, mkvec3(ymin, L, D));
1219 : }
1220 :
1221 : /* Q | N, P = prime divisors of N, R[i] = local epsilon-factor at P[i].
1222 : * Return \prod_{p | Q} R[i] */
1223 : static long
1224 147 : rootno(GEN Q, GEN P, GEN R)
1225 : {
1226 147 : long s = 1, i, l = lg(P);
1227 581 : for (i = 1; i < l; i++)
1228 434 : if (dvdii(Q, gel(P,i))) s *= R[i];
1229 147 : return s;
1230 : }
1231 :
1232 : static void
1233 42 : heegner_find_disc(GEN *points, GEN *coefs, long *pind, GEN E,
1234 : GEN indmult, long ndisc, long prec)
1235 : {
1236 42 : long d = 0;
1237 : GEN faN4, bad, N, faN, listQ, listR;
1238 :
1239 42 : ellQ_get_Nfa(E, &N, &faN);
1240 42 : faN4 = fa_shift2(faN);
1241 42 : listQ = find_div(N, faN);
1242 42 : bad = get_bad(E, gel(faN, 1));
1243 42 : listR = gel(obj_check(E, Q_ROOTNO), 2);
1244 : for(;;)
1245 0 : {
1246 42 : pari_sp av = avma;
1247 42 : GEN list, listD = listDisc(faN4, bad, d, ndisc);
1248 42 : long k, l = lg(listD);
1249 42 : list = cgetg(l, t_VEC);
1250 462 : for (k = 1; k < l; ++k)
1251 420 : gel(list, k) = listheegner(N, faN4, listQ, stoi(listD[k]));
1252 42 : list = vecsort0(list, gen_1, 0);
1253 56 : for (k = l-1; k > 0; --k)
1254 : {
1255 56 : long bprec = 8;
1256 56 : GEN Lk = gel(list,k), D = gel(Lk,3);
1257 56 : GEN sqrtD = sqrtr_abs(itor(D, prec)); /* sqrt(|D|) */
1258 56 : GEN indmultD = heegner_indexmultD(faN, indmult, itos(D), sqrtD);
1259 : do
1260 : {
1261 : GEN mulf, indr;
1262 : pari_timer ti;
1263 56 : if (DEBUGLEVEL) timer_start(&ti);
1264 56 : mulf = ltwist1(E, D, bprec+expo(indmultD));
1265 56 : if (DEBUGLEVEL) timer_printf(&ti,"ellL1twist");
1266 56 : indr = mulrr(indmultD, mulf);
1267 56 : if (DEBUGLEVEL) err_printf("Disc = %Ps, Index^2 = %Ps\n", D, indr);
1268 56 : if (signe(indr)>0 && expo(indr) >= -1) /* indr >=.5 */
1269 : {
1270 : long e, i, l;
1271 42 : GEN pts, cfs, L, indi = grndtoi(sqrtr_abs(indr), &e);
1272 42 : if (e > expi(indi)-7)
1273 : {
1274 0 : bprec++;
1275 0 : pari_warn(warnprec, "ellL1",bprec);
1276 0 : continue;
1277 : }
1278 42 : *pind = itos(indi);
1279 42 : L = gel(Lk, 2); l = lg(L);
1280 42 : pts = cgetg(l, t_VEC);
1281 42 : cfs = cgetg(l, t_VECSMALL);
1282 189 : for (i = 1; i < l; ++i)
1283 : {
1284 147 : GEN P = gel(L,i), z = gel(P,2), Q = gel(P,3); /* [1 or 2, Q] */
1285 : long c;
1286 147 : gel(pts, i) = qfb_root(gel(P,1), sqrtD);
1287 147 : c = rootno(Q, gel(faN,1), listR);
1288 147 : if (!equali1(z)) c *= 2;
1289 147 : cfs[i] = c;
1290 : }
1291 42 : if (DEBUGLEVEL)
1292 0 : err_printf("N = %Ps, ymin*N = %Ps\n",N,
1293 0 : gmul(gsqrt(gel(Lk, 1), prec),N));
1294 42 : *coefs = cfs; *points = pts; return;
1295 : }
1296 : } while(0);
1297 : }
1298 0 : d = listD[l-1]; set_avma(av);
1299 : }
1300 : }
1301 :
1302 : GEN
1303 154 : ellanal_globalred_all(GEN e, GEN *cb, GEN *N, GEN *tam)
1304 : {
1305 154 : GEN E = ellanal_globalred(e, cb), red = obj_check(E, Q_GLOBALRED);
1306 154 : *N = gel(red, 1);
1307 154 : *tam = gel(red,2);
1308 154 : if (signe(ell_get_disc(E))>0) *tam = shifti(*tam,1);
1309 154 : return E;
1310 : }
1311 :
1312 : static GEN
1313 42 : vecelnfembed(GEN x, GEN M, GEN et)
1314 91 : { pari_APPLY_same(gmul(M, etnf_to_basis(et, gel(x,i)))) }
1315 :
1316 : static GEN
1317 42 : QXQV_inv(GEN x, GEN T)
1318 91 : { pari_APPLY_same(QXQ_inv(gel(x,i), T)) }
1319 :
1320 : static GEN
1321 42 : etnfnewprec(GEN x, long prec)
1322 126 : { pari_APPLY_same(nfnewprec(gel(x,i),prec)) }
1323 :
1324 : static GEN
1325 49 : vec_etnf_to_basis(GEN et, GEN x)
1326 154 : { pari_APPLY_same(etnf_to_basis(et,gel(x,i))) }
1327 :
1328 : static GEN
1329 42 : makenfA(GEN sel, GEN A, GEN cb)
1330 : {
1331 42 : GEN etal = gel(sel,1), T = gel(etal,3);
1332 42 : GEN et = gel(etal,1), M = etnf_get_M(et);
1333 42 : long v = etnf_get_varn(et);
1334 42 : GEN al = vecelnfembed(A, M, et);
1335 42 : GEN th = gmul(M, etnf_to_basis(et, pol_x(v)));
1336 42 : return mkvec5(etal,QXQV_inv(A, T),cb,al,th);
1337 : }
1338 :
1339 : GEN
1340 56 : ellheegner(GEN E)
1341 : {
1342 56 : pari_sp av = avma;
1343 : GEN z, P, ht, points, coefs, s, om, indmult;
1344 : GEN sel, etal, et, cbb, A, dAi, T, Ag, At;
1345 : long ind, indx, lint, k, l, wtor, etor, ndisc, ltors2, selrank;
1346 56 : long bitprec = 16, prec = nbits2prec(bitprec)+1;
1347 : pari_timer ti;
1348 : GEN N, cb, tam, torsion, nfA;
1349 56 : E = ellanal_globalred_all(E, &cb, &N, &tam);
1350 56 : if (ellrootno_global(E) == 1)
1351 7 : pari_err_DOMAIN("ellheegner", "(analytic rank)%2","=",gen_0,E);
1352 49 : torsion = elltors(E);
1353 49 : wtor = itos( gel(torsion,1) ); /* #E(Q)_tor */
1354 49 : etor = wtor > 1? itou(gmael(torsion, 2, 1)): 1; /* exponent of E(Q)_tor */
1355 49 : sel = ell2selmer_basis(E, &cbb, prec);
1356 49 : etal = gel(sel,1); A = gel(sel,2); et = gel(etal,1); T = gel(etal,3);
1357 49 : ltors2 = lg(et)-2; selrank = lg(A)-1;
1358 49 : Ag = selrank > ltors2+1 ? pol_1(etnf_get_varn(et)): gel(A,selrank);
1359 49 : At = vecslice(A,1,ltors2);
1360 49 : dAi = gsupnorm(vec_etnf_to_basis(et,A),prec);
1361 : while (1)
1362 42 : {
1363 : GEN hnaive, l1;
1364 : long bitneeded;
1365 91 : if (DEBUGLEVEL) timer_start(&ti);
1366 91 : l1 = ellL1_bitprec(E, 1, bitprec);
1367 91 : if (DEBUGLEVEL) timer_printf(&ti,"ellL1");
1368 91 : if (expo(l1) < 1 - bitprec/2)
1369 7 : pari_err_DOMAIN("ellheegner", "analytic rank",">",gen_1,E);
1370 84 : om = ellR_omega(E,prec);
1371 84 : ht = divrr(mulru(l1, wtor * wtor), mulri(gel(om,1), tam));
1372 84 : if (DEBUGLEVEL) err_printf("Expected height=%Ps\n", ht);
1373 84 : hnaive = hnaive_max(E, ht);
1374 84 : if (DEBUGLEVEL) err_printf("Naive height <= %Ps\n", hnaive);
1375 84 : hnaive = gadd(shiftr(hnaive,-1),glog(dAi,prec));
1376 84 : bitneeded = itos(gceil(divrr(hnaive, mplog2(prec)))) + 12;
1377 84 : if (DEBUGLEVEL) err_printf("precision = %ld\n", bitneeded);
1378 84 : if (bitprec>=bitneeded) break;
1379 42 : bitprec = bitneeded;
1380 42 : prec = nbits2prec(bitprec)+1;
1381 : }
1382 42 : indmult = heegner_indexmult(om, wtor, tam, prec);
1383 42 : ndisc = maxss(10, (long) rtodbl(ht)/10);
1384 42 : heegner_find_disc(&points, &coefs, &ind, E, indmult, ndisc, prec);
1385 42 : if (DEBUGLEVEL) timer_start(&ti);
1386 42 : s = heegner_psi(E, N, points, bitprec);
1387 42 : if (DEBUGLEVEL) timer_printf(&ti,"heegner_psi");
1388 42 : l = lg(points);
1389 42 : z = mulsr(coefs[1], gel(s, 1));
1390 147 : for (k = 2; k < l; ++k) z = addrr(z, mulsr(coefs[k], gel(s, k)));
1391 42 : z = gsub(z, gmul(gel(om,1), ground(gdiv(z, gel(om,1)))));
1392 42 : if (DEBUGLEVEL) err_printf("z=%.*Pg\n",nbits2ndec(bitprec), z);
1393 42 : lint = wtor > 1 ? ugcd(ind, etor): 1;
1394 42 : indx = lint*2*ind;
1395 42 : if (vals(indx) >= vals(etor))
1396 35 : A = mkvec(Ag);
1397 : else
1398 7 : A = mkvec2(Ag, QXQ_mul(Ag, gel(At,1), T));
1399 42 : gmael(sel,1,1) = etnfnewprec(et, prec);
1400 42 : nfA = makenfA(sel, A, cbb);
1401 42 : P = heegner_find_point(E, nfA, om, ht, gmulsg(2*lint, z), indx, prec);
1402 42 : if (DEBUGLEVEL) timer_printf(&ti,"heegner_find_point");
1403 42 : if (cb) P = ellchangepointinv(P, cb);
1404 42 : return gerepilecopy(av, P);
1405 : }
1406 :
1407 : /* Modular degree */
1408 :
1409 : static GEN
1410 70 : ellisobound(GEN e)
1411 : {
1412 70 : GEN M = gel(ellisomat(e,0,1),2);
1413 70 : return vecmax(gel(M,1));
1414 : }
1415 : /* 4Pi^2 / E.area */
1416 : static GEN
1417 140 : getA(GEN E, long prec) { return mpdiv(sqrr(Pi2n(1,prec)), ellR_area(E, prec)); }
1418 :
1419 : /* Modular degree of elliptic curve e over Q, assuming Manin constant = 1
1420 : * (otherwise multiply by square of Manin constant). */
1421 : GEN
1422 70 : ellmoddegree(GEN E)
1423 : {
1424 70 : pari_sp av = avma;
1425 : GEN N, tam, mc2, d;
1426 : long b;
1427 70 : E = ellanal_globalred_all(E, NULL, &N, &tam);
1428 70 : mc2 = sqri(ellisobound(E));
1429 70 : b = expi(mulii(N, mc2)) + maxss(0, expo(getA(E, LOWDEFAULTPREC))) + 16;
1430 : for(;;)
1431 0 : {
1432 70 : long prec = nbits2prec(b), e, s;
1433 70 : GEN deg = mulri(mulrr(lfunellmfpeters(E, b), getA(E, prec)), mc2);
1434 70 : d = grndtoi(deg, &e);
1435 70 : if (DEBUGLEVEL) err_printf("ellmoddegree: %Ps, bit=%ld, err=%ld\n",deg,b,e);
1436 70 : s = expo(deg);
1437 70 : if (e <= -8 && s <= b-8) return gerepileupto(av, gdiv(d,mc2));
1438 0 : b = maxss(s, b+e) + 16;
1439 : }
1440 : }
|
