Line data Source code
1 : /* Copyright (C) 2011 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., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
14 :
15 : #include "pari.h"
16 : #include "paripriv.h"
17 :
18 : static GEN
19 266 : precp_fix(GEN h, long v)
20 266 : { return (precp(h) > v)? cvtop(h,padic_p(h),v): h; }
21 :
22 : /* TATE CURVE */
23 :
24 : /* a1, b1 are t_PADICs, a1/b1 = 1 (mod p) if p odd, (mod 2^4) otherwise.
25 : * Let (A_n, B_n) be defined by A_1 = a1/p^v, B_1 = b1/p^v, v=v(a1)=v(a2);
26 : * A_{n+1} = (A_n + B_n + 2 B_{n+1}) / 4
27 : * B_{n+1} = B_n sqrt(A_n / B_n) = square root of A_n B_n congruent to B_n
28 : * R_n = p^v( A_n - B_n ) = r_{n+1}
29 : * Return [An,Bn,Rn]. N.B. lim An = M2(a1,b1) = M(sqrt(a1),sqrt(b1))^2 */
30 : GEN
31 259 : Qp_agm2_sequence(GEN a1, GEN b1)
32 : {
33 259 : GEN bp, pmod, p = padic_p(a1), q = padic_pd(a1), An, Bn, Rn;
34 259 : long pp = precp(a1), v = valp(a1), i;
35 259 : int pis2 = absequaliu(p,2);
36 259 : a1 = padic_u(a1);
37 259 : b1 = padic_u(b1);
38 259 : if (pis2)
39 154 : pmod = utoipos(8);
40 : else
41 105 : pmod = p;
42 259 : bp = modii(b1, pmod);
43 259 : An = cgetg(pp+1, t_VEC); /* overestimate: rather log_2(pp) */
44 259 : Bn = cgetg(pp+1, t_VEC);
45 259 : Rn = cgetg(pp+1, t_VEC);
46 259 : for(i = 1;; i++)
47 651 : {
48 910 : GEN a = a1, b = b1, r;
49 : long vr;
50 910 : gel(An, i) = a;
51 910 : gel(Bn, i) = b;
52 910 : r = subii(a,b);
53 910 : if (!signe(r)) break;
54 679 : vr = Z_pvalrem(r,p,&r);
55 679 : if (vr >= pp) break;
56 651 : r = cvtop(r, p, pp - vr); setvalp(r, vr+v);
57 651 : gel(Rn, i) = r;
58 :
59 651 : b1 = Zp_sqrt(Fp_mul(a,b,q), p, pp);
60 651 : if (!b1) pari_err_PREC("p-adic AGM");
61 651 : if (!equalii(modii(b1,pmod), bp)) b1 = Fp_neg(b1, q);
62 : /* a1 = (a+b+2sqrt(ab))/4 */
63 651 : if (pis2)
64 : {
65 343 : b1 = remi2n(b1, pp-1);
66 343 : a1 = shifti(addii(addii(a,b), shifti(b1,1)),-2);
67 343 : a1 = remi2n(a1, pp-2);
68 343 : pp -= 2;
69 : }
70 : else
71 308 : a1 = modii(Fp_halve(addii(Fp_halve(addii(a,b),q), b1), q), q);
72 : }
73 259 : setlg(An,i+1);
74 259 : setlg(Bn,i+1);
75 259 : setlg(Rn,i); return mkvec4(An, Bn, Rn, stoi(v));
76 : }
77 : void
78 357 : Qp_descending_Landen(GEN AB, GEN *ptx, GEN *pty)
79 : {
80 357 : GEN R = gel(AB,3);
81 357 : long i, n = lg(R)-1;
82 357 : GEN x = *ptx;
83 357 : if (isintzero(x))
84 : {
85 259 : i = 2;
86 259 : x = gmul2n(gel(R,1),-2);
87 259 : if (pty)
88 : {
89 0 : GEN A = gel(AB,1);
90 0 : if (n == 1)
91 0 : *pty = gmul(x, Qp_sqrt(gadd(x,gel(A,2))));
92 : else
93 0 : *pty = Qp_sqrt(gmul(gmul(x, gadd(x,gel(A,2))), gadd(x,gel(R,2))));
94 0 : if (!*pty) pari_err_PREC("Qp_descending_Landen");
95 : }
96 : }
97 : else
98 98 : i = 1;
99 1008 : for (; i <= n; i++)
100 : {
101 651 : GEN r = gel(R,i), t;
102 651 : if (gequal0(x)) pari_err_PREC("Qp_descending_Landen");
103 651 : t = Qp_sqrt(gaddsg(1, gdiv(r,x))); /* = 1 (mod p) */
104 651 : if (!t) pari_err_PREC("Qp_descending_Landen");
105 651 : if (i == n)
106 : {
107 308 : GEN p = padic_p(r);
108 308 : long v, vx = valp(x), vr = valp(r);
109 308 : if (vx >= vr) pari_err_PREC("Qp_descending_Landen");
110 : /* last loop, take into account loss of accuracy from multiplication
111 : * by \prod_{j > n} sqrt(1+r_j/x_j); since vx < vr, j = n+1 is enough */
112 308 : v = 2*vr - vx;
113 : /* |r_{n+1}| <= |(r_n)^2 / 8| + 1 bit for sqrt loss */
114 308 : if (absequaliu(p,2)) v -= 4;
115 : /* tail is 1 + O(p^v) */
116 308 : if (v < precp(x)) x = cvtop(x,p,v);
117 : }
118 : /* x_{n+1} = x_n ((1 + sqrt(1 + r_n/x_n)) / 2)^2 */
119 651 : x = gmul(x, gsqr(gmul2n(gaddsg(1,t),-1)));
120 : /* y_{n+1} = y_n / (1 - (r_n/4x_{n+1})^2) */
121 651 : if (pty) *pty = gdiv(*pty, gsubsg(1, gsqr(gdiv(r,gmul2n(x,2)))));
122 : }
123 357 : *ptx = x;
124 357 : }
125 : void
126 56 : Qp_ascending_Landen(GEN AB, GEN *ptx, GEN *pty)
127 : {
128 56 : GEN A = gel(AB,1), R = gel(AB,3), x = *ptx, p, r;
129 56 : long n = lg(R)-1, va = itos(gel(AB,4)), v, i;
130 :
131 56 : r = gel(R,n);
132 56 : v = 2*valp(r) + va;
133 56 : if (typ(x) == t_PADIC)
134 35 : v -= 2*valp(x);
135 : else
136 21 : v -= valp(gnorm(x)); /* v(x) = v(Nx) / (e*f), here ef = 2 */
137 56 : p = padic_p(r);
138 56 : if (absequaliu(p,2)) v -= 3; /* |r_{n+1}| <= |(r_n)^2 / 8| */
139 : /* v = v(A[n+1] R[n+1] / x_{n+1}^2) */
140 56 : if (v <= 0) pari_err_PREC("Qp_ascending_Landen");
141 : /* v > 0 => v = v(x_oo) = ... = v(x_{n+1}) */
142 56 : x = gsub(x, gmul2n(r,-1));
143 56 : if (padicprec_relative(x) > v) x = gcvtop(x, p, v);
144 : /* x = x_n */
145 154 : for (i = n; i > 1; i--)
146 : {
147 98 : GEN ar = gmul(gel(A,i),gel(R,i)), xp;
148 98 : setvalp(ar, valp(ar)+va); /* A_i = A[i] * p^va */
149 : /* x_{i-1} = x_i + a_i r_i / x_i - r_{i-1}/2 */
150 98 : xp = gsub(gadd(x, gdiv(ar, x)), gmul2n(gel(R,i-1),-1));
151 : /* y_{i-1} = y_i (1 - a_i r_i / x^2) */
152 98 : if (pty) *pty = gmul(*pty, gsubsg(1, gdiv(ar,gsqr(x))));
153 98 : x = xp;
154 : }
155 56 : *ptx = x;
156 56 : }
157 :
158 : /* Let T = 4x^3 + b2 x^2 + 2b4 x + b6, where T has a unique p-adic root 'a'.
159 : * Return a lift of a to padic accuracy prec. We have
160 : * 216 T = 864 X^3 - 18 c4X - c6, where X = x + b2/12 */
161 : static GEN
162 287 : doellQp_root(GEN E, long prec)
163 : {
164 287 : GEN c4=ell_get_c4(E), c6=ell_get_c6(E), j=ell_get_j(E), p=ellQp_get_p(E);
165 : GEN c6p, T, a;
166 : long alpha;
167 287 : int pis2 = absequaliu(p, 2);
168 287 : if (Q_pval(j, p) >= 0) pari_err_DOMAIN(".root", "v_p(j)", ">=", gen_0, j);
169 : /* v(j) < 0 => v(c4^3) = v(c6^2) = 2 alpha */
170 287 : alpha = Q_pvalrem(ell_get_c4(E), p, &c4) >> 1;
171 287 : if (alpha) (void)Q_pvalrem(ell_get_c6(E), p, &c6);
172 : /* Renormalized so that v(c4) = v(c6) = 0; multiply by p^alpha at the end */
173 287 : if (prec < 4 && pis2) prec = 4;
174 287 : c6p = modii(c6,p);
175 287 : if (pis2)
176 : { /* Use 432T(X/4) = 27X^3 - 9c4 X - 2c6 to have integral root; a=0 mod 2 */
177 140 : T = mkpoln(4, utoipos(27), gen_0, mulis(c4,-9), mulis(c6, -2));
178 : /* v_2(root a) = 1, i.e. will lose one bit of accuracy: prec+1 */
179 140 : a = ZpX_liftroot(T, gen_0, p, prec+1);
180 140 : alpha -= 2;
181 : }
182 147 : else if (absequaliu(p, 3))
183 : { /* Use 216T(X/3) = 32X^3 - 6c4 X - c6 to have integral root; a=-c6 mod 3 */
184 77 : a = Fp_neg(c6p, p);
185 77 : T = mkpoln(4, utoipos(32), gen_0, mulis(c4, -6), negi(c6));
186 77 : a = ZX_Zp_root(T, a, p, prec);
187 77 : switch(lg(a)-1)
188 : {
189 28 : case 1: /* single root */
190 28 : a = gel(a,1); break;
191 49 : case 3: /* three roots, e.g. "15a1", choose the right one */
192 : {
193 49 : GEN a1 = gel(a,1), a2 = gel(a,2), a3 = gel(a,3);
194 49 : long v1 = Z_lval(subii(a2, a3), 3);
195 49 : long v2 = Z_lval(subii(a1, a3), 3);
196 49 : long v3 = Z_lval(subii(a1, a2), 3);
197 49 : if (v1 == v2) a = a3;
198 21 : else if (v1 == v3) a = a2;
199 21 : else a = a1;
200 : }
201 49 : break;
202 : }
203 77 : alpha--;
204 : }
205 : else
206 : { /* p != 2,3: T = 4(x-a)(x-b)^2 = 4x^3 - 3a^2 x - a^3 when b = -a/2
207 : * (so that the trace coefficient vanishes) => a = c6/6c4 (mod p)*/
208 70 : GEN c4p = modii(c4,p);
209 70 : a = Fp_div(c6p, Fp_mulu(c4p, 6, p), p);
210 70 : T = mkpoln(4, utoipos(864), gen_0, mulis(c4, -18), negi(c6));
211 70 : a = ZpX_liftroot(T, a, p, prec);
212 : }
213 287 : a = cvtop(a, p, prec);
214 287 : if (alpha) setvalp(a, valp(a)+alpha);
215 287 : return gsub(a, gdivgu(ell_get_b2(E), 12));
216 : }
217 : GEN
218 539 : ellQp_root(GEN E, long prec)
219 539 : { return obj_checkbuild_padicprec(E, Qp_ROOT, &doellQp_root, prec); }
220 :
221 : /* compute a,b such that E1: y^2 = x(x-a)(x+a-b) ~ E */
222 : static void
223 336 : doellQp_ab(GEN E, GEN *pta, GEN *ptb, long prec)
224 : {
225 336 : GEN b2 = ell_get_b2(E), b4 = ell_get_b4(E), e1 = ellQp_root(E, prec);
226 336 : GEN w, u, t = gadd(gdivgu(b2,4), gmulsg(3,e1)), p = ellQp_get_p(E);
227 336 : w = Qp_sqrt(gmul2n(gadd(b4,gmul(e1,gadd(b2,gmulsg(6,e1)))),1));
228 336 : u = gadd(t,w);
229 : /* Decide between w and -w: we want v(a-b) > v(b) */
230 336 : if (absequaliu(p,2))
231 203 : { if (valp(u)-1 <= valp(w)) w = gneg_i(w); }
232 : else
233 133 : { if (valp(u) <= valp(w)) w = gneg_i(w); }
234 :
235 : /* w^2 = 2b4 + 2b2 e1 + 12 e1^2 = 4(e1-e2)(e1-e3) */
236 336 : *pta = gmul2n(gsub(w,t),-2);
237 336 : *ptb = gmul2n(w,-1);
238 336 : }
239 :
240 : static GEN
241 119 : doellQp_Tate(GEN E, long prec0)
242 : {
243 119 : GEN p = ellQp_get_p(E), j = ell_get_j(E);
244 : GEN L, u, u2, q, x1, a, b, d, s, t, AB, A, M2;
245 119 : long v, n, pp, prec = prec0+3;
246 119 : int split = -1; /* unknown */
247 119 : int pis2 = equaliu(p,2);
248 :
249 119 : if (Q_pval(j, p) >= 0) pari_err_DOMAIN(".tate", "v_p(j)", ">=", gen_0, j);
250 119 : START:
251 336 : doellQp_ab(E, &a, &b, prec);
252 336 : d = gsub(a,b);
253 336 : v = prec0 - precp(d);
254 336 : if (v > 0) { prec += v; goto START; }
255 259 : AB = Qp_agm2_sequence(a,b);
256 259 : A = gel(AB,1);
257 259 : n = lg(A)-1; /* AGM iterations */
258 259 : pp = minss(precp(a),precp(b));
259 259 : M2 = cvtop(gel(A,n), p, pis2? pp-2*n: pp);
260 259 : setvalp(M2, valp(a));
261 259 : u2 = ginv(gmul2n(M2, 2));
262 259 : if (split < 0) split = issquare(u2);
263 259 : x1 = gen_0;
264 259 : Qp_descending_Landen(AB,&x1,NULL);
265 :
266 259 : t = gaddsg(1, ginv(gmul2n(gmul(u2,x1),1)));
267 259 : s = Qp_sqrt(gsubgs(gsqr(t), 1));
268 259 : q = gadd(t,s);
269 259 : if (gequal0(q)) q = gsub(t,s);
270 259 : v = prec0 - precp(q);
271 259 : if (split)
272 : { /* we want log q at precision prec0 */
273 98 : GEN q0 = leafcopy(q); setvalp(q0, 0);
274 98 : v += valp(gsubgs(q0,1));
275 : }
276 259 : if (v > 0) { prec += v; goto START; }
277 119 : if (valp(q) < 0) q = ginv(q);
278 119 : if (split)
279 : {
280 56 : u = Qp_sqrt(u2);
281 56 : L = gdivgs(Qp_log(q), valp(q));
282 : }
283 : else
284 : {
285 63 : GEN T = mkpoln(3, gen_1, gen_0, gneg(u2));
286 63 : u = mkpolmod(pol_x(0), T);
287 63 : L = gen_1;
288 : }
289 119 : return mkvecn(6, u2, u, q, mkvec2(a, b), L, AB);
290 : }
291 : static long
292 686 : Tate_prec(GEN T) { return padicprec_relative(gel(T,3)); }
293 : GEN
294 798 : ellQp_Tate_uniformization(GEN E, long prec)
295 798 : { return obj_checkbuild_prec(E,Qp_TATE,&doellQp_Tate, &Tate_prec,prec); }
296 : GEN
297 203 : ellQp_u(GEN E, long prec)
298 203 : { GEN T = ellQp_Tate_uniformization(E, prec); return gel(T,2); }
299 : GEN
300 77 : ellQp_u2(GEN E, long prec)
301 77 : { GEN T = ellQp_Tate_uniformization(E, prec); return gel(T,1); }
302 : GEN
303 147 : ellQp_q(GEN E, long prec)
304 147 : { GEN T = ellQp_Tate_uniformization(E, prec); return gel(T,3); }
305 : GEN
306 126 : ellQp_ab(GEN E, long prec)
307 126 : { GEN T = ellQp_Tate_uniformization(E, prec); return gel(T,4); }
308 : GEN
309 7 : ellQp_L(GEN E, long prec)
310 7 : { GEN T = ellQp_Tate_uniformization(E, prec); return gel(T,5); }
311 : GEN
312 154 : ellQp_AGM(GEN E, long prec)
313 154 : { GEN T = ellQp_Tate_uniformization(E, prec); return gel(T,6); }
314 :
315 : /* FORMAL GROUP */
316 :
317 : /* t to w := -1/y */
318 : GEN
319 497 : ellformalw(GEN e, long n, long v)
320 : {
321 497 : pari_sp av = avma, av2;
322 : GEN a1,a2,a3,a4,a6, a63;
323 497 : GEN w = cgetg(3, t_SER), t, U, V, W, U2;
324 497 : ulong mask, nold = 1;
325 497 : if (v < 0) v = 0;
326 497 : if (n <= 0) pari_err_DOMAIN("ellformalw","precision","<=",gen_0,stoi(n));
327 483 : mask = quadratic_prec_mask(n);
328 483 : t = pol_x(v);
329 483 : checkell(e);
330 483 : a1 = ell_get_a1(e); a2 = ell_get_a2(e); a3 = ell_get_a3(e);
331 483 : a4 = ell_get_a4(e); a6 = ell_get_a6(e); a63 = gmulgu(a6,3);
332 483 : w[1] = evalsigne(1)|evalvarn(v)|evalvalser(3);
333 483 : gel(w,2) = gen_1; /* t^3 + O(t^4) */
334 : /* use Newton iteration, doubling accuracy at each step
335 : *
336 : * w^3 a6 + w^2(a4 t + a3) + w (a2 t^2 + a1 t - 1) + t^3
337 : * w <- w - -----------------------------------------------------
338 : * w^2 (3a6) + w (2a4 t + 2a3) + (a2 t^2 + a1 t - 1)
339 : *
340 : * w^3 a6 + w^2 U + w V + W
341 : * =: w - -----------------------
342 : * w^2 (3a6) + 2w U + V
343 : */
344 483 : U = gadd(gmul(a4,t), a3);
345 483 : U2 = gmul2n(U,1);
346 483 : V = gsubgs(gadd(gmul(a2,gsqr(t)), gmul(a1,t)), 1);
347 483 : W = gpowgs(t,3);
348 483 : av2 = avma;
349 2261 : while (mask > 1)
350 : { /* nold correct terms in w */
351 1778 : ulong i, nnew = nold << 1;
352 : GEN num, den, wnew, w2, w3;
353 1778 : if (mask & 1) nnew--;
354 1778 : mask >>= 1;
355 1778 : wnew = cgetg(nnew+2, t_SER);
356 1778 : wnew[1] = w[1];
357 7574 : for (i = 2; i < nold+2; i++) gel(wnew,i) = gel(w,i);
358 7210 : for ( ; i < nnew+2; i++) gel(wnew,i) = gen_0;
359 1778 : w = wnew;
360 1778 : w2 = gsqr(w); w3 = gmul(w2,w);
361 1778 : num = gadd(gmul(a6,w3), gadd(gmul(U,w2), gadd(gmul(V,w), W)));
362 1778 : den = gadd(gmul(a63,w2), gadd(gmul(w,U2), V));
363 :
364 1778 : w = gc_upto(av2, gsub(w, gdiv(num, den)));
365 1778 : nold = nnew;
366 : }
367 483 : return gc_GEN(av, w);
368 : }
369 :
370 : static GEN
371 266 : ellformalpoint_i(GEN w, GEN wi)
372 266 : { return mkvec2(gmul(pol_x(varn(w)),wi), gneg(wi)); }
373 :
374 : /* t to [x,y] */
375 : GEN
376 70 : ellformalpoint(GEN e, long n, long v)
377 : {
378 70 : pari_sp av = avma;
379 70 : GEN w = ellformalw(e, n, v), wi = ser_inv(w);
380 70 : return gc_GEN(av, ellformalpoint_i(w, wi));
381 : }
382 :
383 : static GEN
384 343 : ellformaldifferential_i(GEN e, GEN w, GEN wi, GEN *px)
385 : {
386 : GEN x, w1;
387 343 : if (gequal0(ell_get_a1(e)) && gequal0(ell_get_a3(e)))
388 : { /* dx/2y = dx * -w/2, avoid division */
389 147 : x = gmul(pol_x(varn(w)), wi);
390 147 : w1 = gmul(derivser(x), gneg(gmul2n(w,-1)));
391 : }
392 : else
393 : {
394 196 : GEN P = ellformalpoint_i(w, wi);
395 196 : x = gel(P,1);
396 196 : w1 = gdiv(derivser(x), ec_dmFdy_evalQ(e, P));
397 : }
398 343 : *px = x; return w1;
399 : }
400 : /* t to [ dx / (2y + a1 x + a3), x * ... ]*/
401 : GEN
402 70 : ellformaldifferential(GEN e, long n, long v)
403 : {
404 70 : pari_sp av = avma;
405 70 : GEN w = ellformalw(e, n, v), wi = ser_inv(w), x;
406 70 : GEN w1 = ellformaldifferential_i(e, w, wi, &x);
407 70 : return gc_GEN(av, mkvec2(w1,gmul(x,w1)));
408 : }
409 :
410 : /* t to z, dz = w1 dt */
411 : GEN
412 147 : ellformallog(GEN e, long n, long v)
413 : {
414 147 : pari_sp av = avma;
415 147 : GEN w = ellformalw(e, n, v), wi = ser_inv(w), x;
416 147 : GEN w1 = ellformaldifferential_i(e, w, wi, &x);
417 147 : return gc_upto(av, integser(w1));
418 : }
419 : /* z to t */
420 : GEN
421 70 : ellformalexp(GEN e, long n, long v)
422 : {
423 70 : pari_sp av = avma;
424 70 : return gc_upto(av, serreverse(ellformallog(e,n,v)));
425 : }
426 : /* [log_p (sigma(t) / t), log_E t], as power series, d (log_E t) := w1 dt;
427 : * As a fonction of z: odd, = e.b2/12 * z + O(z^3).
428 : * sigma(z) = ellsigma(e) exp(e.b2/24*z^2)
429 : * log_p(sigma(t)/t)=log(subst(sigma(z), x, ellformallog(e))/x) */
430 : static GEN
431 126 : ellformallogsigma_t(GEN e, long n)
432 : {
433 126 : pari_sp av = avma;
434 126 : GEN w = ellformalw(e, n, 0), wi = ser_inv(w), t = pol_x(0);
435 126 : GEN x, s = ellformaldifferential_i(e, w, wi, &x);
436 126 : GEN f = gmul(s, gadd(integser(gmul(x,s)), gmul2n(ell_get_a1(e),-1)));
437 126 : return gc_GEN(av, mkvec2(integser( gsub(ginv(gneg(t)), f) ),
438 : integser(s)));
439 : }
440 :
441 : /* P-ADIC HEIGHTS */
442 :
443 : /* m >= 0, T = b6^2, g4 = b6^2 - b4 b8, return g_m(xP) mod N, in Mazur-Tate's
444 : * notation (Duke 1991)*/
445 : static GEN
446 4102 : rellg(hashtable *H, GEN m, GEN T, GEN g4, GEN b8, GEN N)
447 : {
448 : hashentry *h;
449 : GEN n, z, np2, np1, nm2, nm1, fp2, fp1, fm2, fm1, f;
450 : ulong m4;
451 4102 : if (abscmpiu(m, 4) <= 0) switch(itou(m))
452 : {
453 196 : case 0: return gen_0;
454 532 : case 1: return gen_1;
455 630 : case 2: return subiu(N,1);
456 728 : case 3: return b8;
457 756 : case 4: return g4;
458 : }
459 1260 : if ((h = hash_search(H, (void*)m))) return (GEN)h->val;
460 686 : m4 = mod4(m);
461 686 : n = shifti(m, -1); f = rellg(H,n,T,g4,b8,N);
462 686 : np2 = addiu(n, 2); fp2 = rellg(H,np2,T,g4,b8,N);
463 686 : np1 = addiu(n, 1); fp1 = rellg(H,np1,T,g4,b8,N);
464 686 : nm2 = subiu(n, 2); fm2 = rellg(H,nm2,T,g4,b8,N);
465 686 : nm1 = subiu(n, 1); fm1 = rellg(H,nm1,T,g4,b8,N);
466 686 : if (odd(m4))
467 : {
468 406 : GEN t1 = Fp_mul(fp2, Fp_powu(f,3,N), N);
469 406 : GEN t2 = Fp_mul(fm1, Fp_powu(fp1,3,N), N);
470 406 : if (mpodd(n)) t2 = Fp_mul(T,t2,N); else t1 = Fp_mul(T,t1,N);
471 406 : z = Fp_sub(t1, t2, N);
472 : }
473 : else
474 : {
475 280 : GEN t1 = Fp_mul(fm2, Fp_sqr(fp1,N), N);
476 280 : GEN t2 = Fp_mul(fp2, Fp_sqr(fm1,N), N);
477 280 : z = Fp_mul(f, Fp_sub(t1, t2, N), N);
478 : }
479 686 : hash_insert(H, (void*)m, (void*)z);
480 686 : return z;
481 : }
482 :
483 : static GEN
484 672 : addii3(GEN x, GEN y, GEN z) { return addii(x,addii(y,z)); }
485 : static GEN
486 448 : addii4(GEN x, GEN y, GEN z, GEN t) { return addii(x,addii3(y,z,t)); }
487 : static GEN
488 224 : addii5(GEN x, GEN y, GEN z, GEN t, GEN u) { return addii(x,addii4(y,z,t,u)); }
489 :
490 : /* xP = [n,d] (corr. to n/d, coprime), such that the reduction of the point
491 : * P = [xP,yP] is non singular at all places. Return x([m] P) mod N as
492 : * [num,den] (coprime) */
493 : static GEN
494 224 : xmP(GEN e, GEN xP, GEN m, GEN N)
495 : {
496 224 : pari_sp av = avma;
497 224 : ulong k = expi(m);
498 224 : hashtable *H = hash_create_GEN((5+k)*k, 1);
499 224 : GEN b2 = ell_get_b2(e), b4 = ell_get_b4(e), n = gel(xP,1), d = gel(xP,2);
500 224 : GEN b6 = ell_get_b6(e), b8 = ell_get_b8(e);
501 : GEN B4, B6, B8, T, g4;
502 224 : GEN d2 = Fp_sqr(d,N), d3 = Fp_mul(d2,d,N), d4 = Fp_sqr(d2,N);
503 224 : GEN n2 = Fp_sqr(n,N), n3 = Fp_mul(n2,n,N), n4 = Fp_sqr(n2,N);
504 224 : GEN nd = Fp_mul(n,d,N), n2d2 = Fp_sqr(nd,N);
505 224 : GEN b2nd = Fp_mul(b2,nd, N), b2n2d = Fp_mul(b2nd,n,N);
506 224 : GEN b6d3 = Fp_mul(b6,d3,N), g,gp1,gm1, C,D;
507 224 : B8 = addii5(muliu(n4,3), mulii(b2n2d,n), mulii(muliu(b4,3), n2d2),
508 : mulii(muliu(b6d3,3), n), mulii(b8,d4));
509 224 : B6 = addii4(muliu(n3,4), mulii(b2nd,n),
510 : shifti(mulii(b4,Fp_mul(n,d2,N)), 1), b6d3);
511 224 : B4 = addii3(muliu(n2,6), b2nd, mulii(b4,d2));
512 224 : B4 = modii(B4,N);
513 224 : B6 = modii(B6,N);
514 224 : B8 = modii(B8,N);
515 224 : g4 = Fp_sub(sqri(B6), mulii(B4,B8), N);
516 224 : T = Fp_sqr(B6,N);
517 :
518 224 : g = rellg(H, m, T,g4,B8, N);
519 224 : gp1 = rellg(H, addiu(m,1), T,g4,B8, N);
520 224 : gm1 = rellg(H, subiu(m,1), T,g4,B8, N);
521 224 : C = Fp_sqr(g, N);
522 224 : D = Fp_mul(gp1,gm1, N);
523 224 : if (mpodd(m))
524 : {
525 112 : n = Fp_sub(mulii(C,n), mulii(D,B6), N);
526 112 : d = Fp_mul(C,d, N);
527 : }
528 : else
529 : {
530 112 : n = Fp_sub(Fp_mul(Fp_mul(B6,C,N), n, N), D, N);
531 112 : d = Fp_mul(Fp_mul(C,d,N), B6, N);
532 : }
533 224 : return gc_GEN(av, mkvec2(n,d));
534 : }
535 : /* given [n,d2], x = n/d2 (coprime, d2 = d^2), p | den,
536 : * return t = -x/y + O(p^v) */
537 : static GEN
538 126 : tfromx(GEN e, GEN x, GEN p, long v, GEN N, GEN *pd)
539 : {
540 126 : GEN n = gel(x,1), d2 = gel(x,2), d;
541 : GEN a1, a3, b2, b4, b6, B, C, d4, d6, Y;
542 126 : if (!signe(n)) { *pd = gen_1; return zeropadic_shallow(p, v); }
543 126 : a1 = ell_get_a1(e);
544 126 : b2 = ell_get_b2(e);
545 126 : a3 = ell_get_a3(e);
546 126 : b4 = ell_get_b4(e);
547 126 : b6 = ell_get_b6(e);
548 126 : d = Qp_sqrt(cvtop(d2, p, v - Z_pval(d2,p)));
549 126 : if (!d) pari_err_BUG("ellpadicheight");
550 : /* Solve Y^2 = 4n^3 + b2 n^2 d2+ 2b4 n d2^2 + b6 d2^3,
551 : * Y = 2y + a1 n d + a3 d^3 */
552 126 : d4 = Fp_sqr(d2, N);
553 126 : d6 = Fp_mul(d4, d2, N);
554 126 : B = gmul(d, Fp_add(mulii(a1,n), mulii(a3,d2), N));
555 126 : C = mkpoln(4, utoipos(4), Fp_mul(b2, d2, N),
556 : Fp_mul(shifti(b4,1), d4, N),
557 : Fp_mul(b6,d6,N));
558 126 : C = FpX_eval(C, n, N);
559 126 : if (!signe(C))
560 0 : Y = zeropadic_shallow(p, v >> 1);
561 : else
562 126 : Y = Qp_sqrt(cvtop(C, p, v - Z_pval(C,p)));
563 126 : if (!Y) pari_err_BUG("ellpadicheight");
564 126 : *pd = d;
565 126 : return gdiv(gmulgs(gmul(n,d), -2), gsub(Y,B));
566 : }
567 :
568 : /* return minimal i s.t. -v_p(j+1) - log_p(j-1) + (j+1)*t >= v for all j>=i */
569 : static long
570 126 : logsigma_prec(GEN p, long v, long t)
571 : {
572 126 : double log2p = dbllog2(p);
573 126 : long j, i = ceil((v - t) / (t - 2*M_LN2/(3*log2p)) + 0.01);
574 126 : if (absequaliu(p,2) && i < 5) i = 5;
575 : /* guaranteed to work, now optimize */
576 182 : for (j = i-1; j >= 2; j--)
577 : {
578 182 : if (- u_pval(j+1,p) - log2(j-1)/log2p + (j+1)*t + 0.01 < v) break;
579 56 : i = j;
580 : }
581 126 : if (j == 1)
582 : {
583 0 : if (- absequaliu(p,2) + 2*t + 0.01 >= v) i = 1;
584 : }
585 126 : return i;
586 : }
587 : /* return minimal i s.t. -v_p(j+1) + (j+1)*t >= v for all j>=i */
588 : static long
589 7 : log_prec(GEN p, long v, long t)
590 : {
591 7 : double log2p = dbllog2(p);
592 7 : long j, i = ceil(v / (t - M_LN2/(2*log2p)) + 0.01);
593 : /* guaranteed to work, now optimize */
594 28 : for (j = i-1; j >= 1; j--)
595 : {
596 28 : if (- u_pval(j+1,p) + (j+1)*t + 0.01 < v) break;
597 21 : i = j;
598 : }
599 7 : return i;
600 : }
601 :
602 : /* P = rational point of exact denominator d. Is Q singular on E(Fp) ? */
603 : static int
604 182 : FpE_issingular(GEN E, GEN P, GEN d, GEN p)
605 : {
606 182 : pari_sp av = avma;
607 : GEN t, x, y, a1, a2, a3, a4;
608 182 : if (ell_is_inf(E) || dvdii(d,p)) return 0; /* 0_E is smooth */
609 175 : P = Q_muli_to_int(P,d);
610 175 : x = gel(P,1);
611 175 : y = gel(P,2);
612 175 : a1 = ell_get_a1(E);
613 175 : a3 = ell_get_a3(E);
614 175 : t = addii(shifti(y,1), addii(mulii(a1,x), mulii(a3,d)));
615 175 : if (!dvdii(t,p)) return gc_bool(av, 0);
616 35 : a2 = ell_get_a2(E);
617 35 : a4 = ell_get_a4(E);
618 35 : d = Fp_inv(d, p);
619 35 : x = Fp_mul(x,d,p);
620 35 : y = Fp_mul(y,d,p);
621 35 : t = subii(mulii(a1,y), addii(a4, mulii(x, addii(gmul2n(a2,1), muliu(x,3)))));
622 35 : return gc_bool(av, dvdii(t,p));
623 : }
624 : /* E/Q, P on E(Q). Let g > 0 minimal such that the image of R = [g]P in a
625 : * minimal model is everywhere nonsingular; return [R,g] */
626 : GEN
627 140 : ellnonsingularmultiple(GEN e, GEN P)
628 : {
629 140 : pari_sp av = avma;
630 140 : GEN ch, E = ellanal_globalred(e, &ch), NP, L, S, d, g = gen_1;
631 : long i, l;
632 140 : if (!checkellpt_i(P)) pari_err_TYPE("ellnonsingularmultiple",P);
633 140 : if (ell_is_inf(P)) retmkvec2(gcopy(P), gen_1);
634 140 : if (E != e) P = ellchangepoint(P, ch);
635 140 : S = obj_check(E, Q_GLOBALRED);
636 140 : NP = gmael(S,3,1);
637 140 : L = gel(S,4);
638 140 : l = lg(NP);
639 140 : d = Q_denom(P);
640 315 : for (i = 1; i < l; i++)
641 : {
642 175 : GEN G = gel(L,i), p = gel(NP,i);/* prime of bad reduction */
643 : long kod;
644 175 : if (!FpE_issingular(E, P, d, p)) continue;
645 21 : kod = itos(gel(G,2)); /* Kodaira type */
646 21 : if (kod >= 5) /* I_nu */
647 : {
648 7 : long nu = kod - 4;
649 7 : long n = minss(Q_pval(ec_dmFdy_evalQ(E,P), p), nu/2);
650 7 : nu /= ugcd(nu, n);
651 7 : g = muliu(g, nu);
652 7 : P = ellmul(E, P, utoipos(nu)); d = Q_denom(P);
653 14 : } else if (kod <= -5) /* I^*_nu */
654 : { /* either 2 or 4 */
655 7 : long nu = - kod - 4;
656 7 : P = elladd(E, P,P); d = Q_denom(P);
657 7 : g = shifti(g,1);
658 7 : if (odd(nu) && FpE_issingular(E, P, d, p))
659 : { /* it's 4 */
660 7 : P = elladd(E, P,P); d = Q_denom(P);
661 7 : g = shifti(g,1);
662 : }
663 : } else {
664 7 : GEN c = gel(G, 4); /* Tamagawa number at p */
665 7 : if (absequaliu(c, 4)) c = gen_2;
666 7 : P = ellmul(E, P, c); d = Q_denom(P);
667 7 : g = mulii(g, c);
668 : }
669 : }
670 140 : if (E != e) P = ellchangepointinv(P, ch);
671 140 : return gc_GEN(av, mkvec2(P,g));
672 : }
673 :
674 : GEN
675 140 : ellpadicheight(GEN e, GEN p, long v0, GEN P)
676 : {
677 140 : pari_sp av = avma;
678 : GEN N, H, h, t, ch, g, E, x, D, ls, lt, S, a,b;
679 : long v, vd;
680 : int is2;
681 140 : if (!checkellpt_i(P)) pari_err_TYPE("ellpadicheight",P);
682 140 : if (v0<=0) pari_err_DOMAIN("ellpadicheight","precision","<=",gen_0,stoi(v0));
683 133 : checkell_Q(e);
684 133 : if (typ(p) != t_INT) pari_err_TYPE("ellpadicheight",p);
685 133 : if (cmpiu(p,2) < 0) pari_err_PRIME("ellpadicheight",p);
686 133 : if (ellorder_Q(e,P)) return mkvec2(gen_0,gen_0);
687 126 : E = ellanal_globalred(e, &ch);
688 126 : if (E != e) P = ellchangepoint(P, ch);
689 126 : S = ellnonsingularmultiple(E, P);
690 126 : P = gel(S,1);
691 126 : g = gel(S,2);
692 126 : v = v0 + 2*Z_pval(g, p);
693 126 : is2 = absequaliu(p,2); if (is2) v += 2;
694 126 : x = Q_remove_denom(gel(P,1), &b);
695 126 : x = mkvec2(x, b? b: gen_1);
696 126 : vd = Z_pval(gel(x,2), p);
697 126 : if (!vd)
698 : { /* P not in kernel of reduction mod p */
699 112 : GEN d, m, X, Pp, Ep = ellinit(E, mkintmod(gen_0,p), 0);
700 112 : long w = v+2;
701 112 : Pp = RgV_to_FpV(P, p);
702 112 : if (lg(Ep) != 1)
703 105 : m = ellorder(Ep, Pp, NULL);
704 : else
705 : { /* E has bad reduction at p */
706 7 : m = ellcard(E, p);
707 7 : if (equalii(m, p)) pari_err_TYPE("ellpadicheight: additive reduction", E);
708 : }
709 112 : g = mulii(g,m);
710 : for(;;)
711 : {
712 224 : N = powiu(p, w);
713 224 : X = xmP(E, x, m, N);
714 224 : d = gel(X,2);
715 224 : if (!signe(d)) w <<= 1;
716 : else
717 : {
718 224 : vd = Z_pval(d, p);
719 224 : if (w >= v+2*vd + is2) break;
720 112 : w = v+2*vd + is2;
721 : }
722 : }
723 112 : x = X;
724 : }
725 : /* we will want t mod p^(v+vd) because of t/D in H later, and
726 : * we lose p^vd in tfromx because of sqrt(d) (p^(vd+1) if p=2)*/
727 126 : v += 2*vd + is2;
728 126 : N = powiu(p,v);
729 126 : t = tfromx(E, x, p, v, N, &D); /* D^2=denom(x)=x[2] */
730 126 : S = ellformallogsigma_t(E, logsigma_prec(p, v-vd, valp(t)) + 1);
731 126 : ls = ser2rfrac_i(gel(S,1)); /* log_p (sigma(T)/T) */
732 126 : lt = ser2rfrac_i(gel(S,2)); /* log_E (T) */
733 : /* evaluate our formal power series at t */
734 126 : H = gadd(poleval(ls, t), glog(gdiv(t, D), 0));
735 126 : h = gsqr(poleval(lt, t));
736 126 : g = sqri(g);
737 126 : a = gdiv(gmulgs(H,-2), g);
738 126 : b = gdiv(gneg(h), g);
739 126 : if (E != e)
740 : {
741 21 : GEN u = gel(ch,1), r = gel(ch,2);
742 21 : a = gdiv(gadd(a, gmul(r,b)), u);
743 21 : b = gmul(u,b);
744 : }
745 126 : H = mkvec2(a,b);
746 126 : gel(H,1) = precp_fix(gel(H,1),v0);
747 126 : gel(H,2) = precp_fix(gel(H,2),v0);
748 126 : return gc_GEN(av, H);
749 : }
750 :
751 : GEN
752 14 : ellpadiclog(GEN E, GEN p, long n, GEN P)
753 : {
754 14 : pari_sp av = avma;
755 : long vt;
756 : GEN t, x, y, L;
757 14 : if (!checkellpt_i(P)) pari_err_TYPE("ellpadiclog",P);
758 14 : if (ell_is_inf(P)) return gen_0;
759 14 : x = gel(P,1);
760 14 : y = gel(P,2); t = gneg(gdiv(x,y));
761 14 : vt = gvaluation(t, p); /* can be a t_INT, t_FRAC or t_PADIC */
762 14 : if (vt <= 0)
763 7 : pari_err_DOMAIN("ellpadiclog","P","not in the kernel of reduction at",p,P);
764 7 : L = ser2rfrac_i(ellformallog(E, log_prec(p, n, vt) + 1, 0));
765 7 : return gc_upto(av, poleval(L, cvtop(t, p, n)));
766 : }
767 :
768 : /* E/Qp has multiplicative reduction, Tate curve */
769 : static GEN
770 14 : ellpadics2_tate(GEN Ep, long n)
771 : {
772 : pari_sp av;
773 14 : GEN u2 = ellQp_u2(Ep, n), q = ellQp_q(Ep, n), pn = padic_pd(q);
774 14 : GEN qm, s, b2 = ell_get_b2(Ep), v = vecfactoru_i(1, n);
775 : long m;
776 14 : qm = Fp_powers(padic_to_Fp(q, pn), n, pn);
777 14 : s = gel(qm, 2); av = avma;
778 70 : for (m = 2; m <= n; m++) /* sum sigma(m) q^m */
779 : {
780 56 : s = addii(s, mulii(gel(qm,m+1), usumdiv_fact(gel(v,m))));
781 56 : if ((m & 31) == 0) s = gc_INT(av, s);
782 : }
783 14 : s = subui(1, muliu(s,24));
784 14 : s = gdivgu(gsub(b2, gdiv(s,u2)), 12);
785 14 : return precp_fix(s,n);
786 : }
787 :
788 : GEN
789 189 : ellpadicfrobenius(GEN E, ulong p, long n)
790 : {
791 189 : checkell(E);
792 189 : if (p < 2) pari_err_DOMAIN("ellpadicfrobenius","p", "<", gen_2, utoipos(p));
793 189 : switch(ell_get_type(E))
794 : {
795 175 : case t_ELL_Q: break;
796 14 : case t_ELL_Qp: if (equaliu(ellQp_get_p(E), p)) break;
797 7 : default: pari_err_TYPE("ellpadicfrobenius", E);
798 : }
799 182 : return hyperellpadicfrobenius(ec_bmodel(E,0), p, n);
800 : }
801 :
802 : /* s2 = (b_2-E_2)/12 */
803 : GEN
804 35 : ellpadics2(GEN E, GEN p, long n)
805 : {
806 35 : pari_sp av = avma;
807 : GEN l, F, a,b,d, ap;
808 : ulong pp;
809 35 : if (typ(p) != t_INT) pari_err_TYPE("ellpadics2",p);
810 35 : if (cmpis(p,2) < 0) pari_err_PRIME("ellpadics2",p);
811 35 : checkell(E);
812 35 : if (Q_pval(ell_get_j(E), p) < 0)
813 : {
814 : GEN Ep;
815 7 : if (ell_get_type(E) == t_ELL_Qp) Ep = E;
816 0 : else Ep = ellinit(E, zeropadic_shallow(p,n), 0);
817 7 : l = ellpadics2_tate(Ep, n);
818 7 : if (Ep != E) obj_free(Ep);
819 7 : return gc_GEN(av, l);
820 : }
821 28 : pp = itou(p);
822 28 : F = ellpadicfrobenius(E, pp, n);
823 21 : a = gcoeff(F,1,1);
824 21 : b = gcoeff(F,1,2);
825 21 : d = gcoeff(F,2,2); ap = gadd(a,d);
826 21 : if (valp(ap) > 0) pari_err_DOMAIN("ellpadics2","E","is supersingular at",p,E);
827 21 : if (pp == 2 || (pp <= 13 && n == 1)) /* 2sqrt(p) > p/2: ambiguity */
828 0 : ap = ellap(E,p);
829 : else
830 : { /* either 2sqrt(p) < p/2 or n > 1 and 2sqrt(p) < p^2/2 (since p!=2) */
831 21 : GEN q = pp <= 13? utoipos(pp * pp): p;
832 21 : ap = padic_to_Fp(ap, q);
833 21 : ap = Fp_center_i(ap, q, shifti(q,-1));
834 : }
835 21 : l = mspadic_unit_eigenvalue(ap, 2, p, n);
836 21 : return gc_upto(av, gdiv(b, gsub(l, a))); /* slope of eigenvector */
837 : }
838 :
839 : /* symbol and modular symbol space attached to E to later compute
840 : * ellpadicL(E,p,, s,,D) */
841 : static GEN
842 91 : ellpadicL_symbol(GEN E, GEN p, GEN s, GEN D)
843 : {
844 : GEN s1, s2, ap;
845 : long sign;
846 91 : checkell(E);
847 91 : if (ell_get_type(E) != t_ELL_Q) pari_err_TYPE("ellpadicL",E);
848 91 : ap = ellap(E,p);
849 91 : if (D && typ(D) != t_INT) pari_err_TYPE("ellpadicL",D);
850 91 : if (D && !Z_isfundamental(D))
851 0 : pari_err_DOMAIN("ellpadicL", "isfundamental(D)", "=", gen_0, D);
852 91 : if (!D) D = gen_1;
853 91 : if (Z_pval(ellQ_get_N(E), p) >= 2)
854 0 : pari_err_IMPL("additive reduction in ellpadicL");
855 91 : mspadic_parse_chi(s, &s1,&s2);
856 91 : sign = signe(D); if (mpodd(s2)) sign = -sign;
857 91 : return shallowconcat(msfromell(E, sign), mkvec4(ap, p, s, D));
858 : }
859 : /* W an ellpadicL_symbol, initialize for ellpadicL(E,p,n,s,,D) */
860 : static GEN
861 91 : ellpadicL_init(GEN W, long n)
862 : {
863 91 : GEN Wp, den, M = gel(W,1), xpm = gel(W,2), ap = gel(W,3), s = gel(W,5);
864 91 : long p = itos(gel(W,4)), D = itos(gel(W,6));
865 :
866 91 : xpm = Q_remove_denom(xpm,&den); if (!den) den = gen_1;
867 91 : n += Z_lval(den,p);
868 91 : Wp = mspadicinit(M, p, n, Z_lval(ap,p));
869 91 : return mkvec3(mspadicmoments(Wp,xpm,D), den, s);
870 : }
871 : /* v from ellpadicL_init, compute ellpadicL(E,p,n,s,r,D) */
872 : static GEN
873 126 : ellpadic_i(GEN v, long r)
874 : {
875 126 : GEN oms = gel(v,1), den = gel(v,2), s = gel(v,3);
876 126 : return gdiv(mspadicL(oms,s,r), den);
877 : }
878 : GEN
879 63 : ellpadicL(GEN E, GEN p, long n, GEN s, long r, GEN D)
880 : {
881 63 : pari_sp av = avma;
882 : GEN W, v;
883 63 : if (r < 0) pari_err_DOMAIN("ellpadicL","r","<",gen_0,stoi(r));
884 63 : if (n <= 0) pari_err_DOMAIN("ellpadicL","precision","<=",gen_0,stoi(n));
885 63 : W = ellpadicL_symbol(E, p, s, D);
886 63 : v = ellpadicL_init(W, n);
887 63 : return gc_upto(av, ellpadic_i(v, r));
888 : }
889 :
890 : static long
891 28 : torsion_order(GEN E) { GEN T = elltors(E); return itos(gel(T,1)); }
892 :
893 : /* E given by a minimal model; D != 0. Compare Euler factor of L(E,(D/.),1)
894 : * with L(E^D,1). Return
895 : * \prod_{p|D} (p-a_p(E)+eps_{E}(p)) / p,
896 : * where eps(p) = 0 if p | N_E and 1 otherwise */
897 : static GEN
898 7 : get_Euler(GEN E, GEN D)
899 : {
900 7 : GEN a = gen_1, b = gen_1, P = gel(absZ_factor(D), 1);
901 7 : long i, l = lg(P);
902 14 : for (i = 1; i < l; i++)
903 : {
904 7 : GEN p = gel(P,i);
905 7 : a = mulii(a, ellcard(E, p));
906 7 : b = mulii(b, p);
907 : }
908 7 : return Qdivii(a, b);
909 : }
910 : GEN
911 28 : ellpadicbsd(GEN E, GEN p, long n, GEN D)
912 : {
913 28 : const long MAXR = 30;
914 28 : pari_sp av = avma;
915 28 : GEN D0, W, ED, tam, ND, C, apD, U = NULL;/*-Wall*/
916 : long r, vN;
917 28 : checkell(E);
918 28 : if (D && isint1(D)) D = NULL;
919 28 : D0 = ellminimaltwistcond(E); if (isint1(D0)) D0 = NULL;
920 28 : if (D0) { E = elltwist(E, D0); D = D? coredisc(mulii(D, D0)): D0; }
921 28 : W = ellpadicL_symbol(E, p, gen_0, D);
922 28 : ED = D? elltwist(E, D): E;
923 28 : ED = ellanal_globalred_all(ED, NULL, &ND, &tam);
924 28 : vN = Z_pval(ND, p); /* additive reduction ? */
925 28 : if (vN >= 2) pari_err_DOMAIN("ellpadicbsd","v_p(N(E_D))",">",gen_1,stoi(vN));
926 28 : if (n < 5) n = 5;
927 0 : for(;; n <<= 1)
928 0 : {
929 28 : pari_sp av2 = avma;
930 28 : GEN v = ellpadicL_init(W, n);
931 63 : for (r = 0; r < MAXR; r++)
932 : {
933 63 : U = ellpadic_i(v, r);
934 63 : if (!gequal0(U)) break;
935 : }
936 28 : if (r < MAXR) break;
937 0 : set_avma(av2);
938 : }
939 28 : apD = ellap(ED, p);
940 28 : if (typ(U) == t_COL)
941 : { /* p | a_p(E_D), frobenius on E_D */
942 0 : GEN F = mkmat22(gen_0, negi(p), gen_1, apD);
943 0 : U = RgM_RgC_mul(gpowgs(gsubsg(1, gdiv(F,p)), -2), U);
944 0 : settyp(U, t_VEC);
945 : }
946 28 : else if (dvdii(ND,p))
947 : {
948 21 : if (equalim1(apD)) /* divide by 1-1/a_p */
949 14 : U = gdivgu(U, 2);
950 : else
951 : { /* ap = 1 */
952 7 : GEN EDp = ellinit(ED, zeropadic_shallow(p,n), 0);
953 7 : U = gdiv(U, ellQp_L(EDp,n));
954 7 : obj_free(EDp);
955 : }
956 : }
957 : else
958 : {
959 7 : GEN a = mspadic_unit_eigenvalue(apD, 2, p, n);
960 7 : U = gmul(U, gpowgs(gsubsg(1, ginv(a)), -2));
961 : }
962 28 : C = mulii(tam, mpfact(r));
963 28 : if (D) C = gmul(C, get_Euler(ED, D));
964 28 : C = gdiv(sqru(torsion_order(ED)), C);
965 28 : if (D) obj_free(ED);
966 28 : return gc_GEN(av, mkvec2(utoi(r), gmul(U, C)));
967 : }
968 :
969 : GEN
970 7 : ellpadicregulator(GEN E, GEN p, long n, GEN S)
971 : {
972 7 : pari_sp av = avma;
973 7 : GEN FG = ellpadicheightmatrix(E,p,n,S); /* forbids additive reduction */
974 : /* [F,G]: height in basis [omega, eta] */
975 7 : GEN R, F = gel(FG,1), G = gel(FG,2), ap = ellap(E,p);
976 7 : if (dvdii(ap, p))
977 : { /* supersingular */
978 0 : GEN f = ellpadicfrobenius(E, itou(p), n);
979 0 : GEN x = gcoeff(f,1,1), y = gcoeff(f,2,1);
980 : /* [A,B]: regulator height in basis [omega, eta] */
981 0 : GEN A = det(F), B = det(gadd(F,G)), C;
982 0 : C = gdiv(gsub(B,A), y);
983 : /* R: regulator height in basis [omega, F.omega] */
984 0 : R = mkvec2(gsub(A, gmul(x,C)), C);
985 : }
986 : else
987 : {
988 : GEN s2;
989 7 : if (equali1(ap) && dvdii(ell_get_disc(E),p))
990 7 : { /* split multiplicative reduction */
991 7 : GEN Ep = ellinit(E, zeropadic_shallow(p,n), 0);
992 7 : GEN q = ellQp_q(Ep,n), u2 = ellQp_u2(Ep,n);
993 7 : s2 = ellpadics2_tate(Ep, n);
994 7 : s2 = gsub(s2, ginv(gmul(Qp_log(q), u2))); /*extended MW group contrib*/
995 7 : obj_free(Ep);
996 : }
997 : else
998 0 : s2 = ellpadics2(E,p,n);
999 7 : R = det( RgM_sub(F, RgM_Rg_mul(G,s2)) );
1000 : }
1001 7 : return gc_GEN(av, R);
1002 : }
|
