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 : /** Batch p-adic logarithms **/
19 : /* a/b mod q */
20 : static ulong
21 1518 : divmodulo(ulong a, ulong b, ulong p, ulong q)
22 : {
23 1518 : long v = u_lvalrem(b, p, &b);
24 1518 : if (v) a /= upowuu(p,v);
25 : /* a/b is now a p-integer */
26 1518 : return Fl_div(a, b, q);
27 : }
28 : /* to compute log_p(a) mod q = p^n, p < 2^31 */
29 : static GEN
30 732 : initQplog(long p, ulong q, long n)
31 : {
32 : long i, nn, nt;
33 : ulong a, P;
34 : GEN C;
35 1128 : for(nn = n, nt = n + 1; nn >= p; nn /= p) nt++;
36 732 : if (p == 2)
37 : {
38 60 : P = q - 8;
39 324 : while (3 * (nt - 1) > u_lval(nt-1, p) + n) nt--;
40 : }
41 : else
42 : {
43 672 : P = q - p;
44 1500 : while (nt > u_lval(nt-1, p) + n) nt--;
45 : }
46 732 : C = cgetg(nt, t_VECSMALL);
47 : /* [ P^(k-1) / k, k = 1..nt ] / (p - 1) */
48 2250 : for(i = 1, a = 1; i < nt; i++, a = Fl_mul(a, P, q))
49 1518 : C[i] = divmodulo(a, i * (p-1), p, q);
50 732 : return C;
51 : }
52 : /* compute log_p(a) / p (resp. log_2(a) / 4) mod p^n, 'a' a p-unit,
53 : * C = initQplog(p,q,n), q = p^n */
54 : static ulong
55 1404492 : logp(GEN C, ulong a, ulong p, ulong q, ulong pn)
56 : {
57 1404492 : long i, b, z, n = lg(C)-1;
58 1404492 : a %= q;
59 : /* Euclidean quotient (a^2 = 1 mod 8, a^(p-1) = 1 mod p) */
60 1404492 : if (p == 2)
61 408 : b = Fl_sqr(a, q << 1) >> 3;
62 : else
63 1404084 : b = Fl_powu(a, p-1, q) / p;
64 1404492 : z = Fl_mul(b, C[n], pn);
65 4589820 : for (i = n-1; i > 0; i--) z = Fl_mul(z + C[i], b, pn);
66 1404492 : return z;
67 : }
68 :
69 : /* p-adic integration against d mu_E, mod p^n, m > 0. D fundamental,
70 : * W a msinit for k = 2, xpm an eigensymbol.
71 : * Assume |D p^n| < MAX_LONG. Much cheaper than oms at low precision.
72 : * Return 1/2 * correct value [half loop over a]. Assume xpm belongs to
73 : * (-1)^i (D/-1) - part */
74 : static GEN
75 1038 : polPn(GEN W, GEN xpm, long p, long D, long R, long n)
76 : {
77 1038 : pari_sp av = avma, av2;
78 1038 : long N = (p == 2)? n + 2: n + 1;
79 1038 : ulong aD = labs(D), pn = upowuu(p, n), q = upowuu(p, N), ic = 0;
80 1038 : GEN v, Dq = muluu(aD, q), nc = icopy(gen_1), c = mkfrac(nc, Dq);
81 1038 : GEN C = n? initQplog(p, q, N): NULL;
82 1038 : GEN tab = R? ZV_to_Flv(teichmullerinit(p, N), q): NULL;
83 1038 : long a, A = itou(shifti(Dq,-1));
84 :
85 1038 : if (n) ic = Fl_inv(logp(C, 1 + (p == 2? 4: p), p, q, pn), pn);
86 1038 : v = zero_zv(pn); av2 = avma;
87 3102702 : for (a = 1; a <= A; a++, set_avma(av2))
88 : {
89 : GEN X;
90 3101664 : long s, ap = a % p;
91 : ulong x, j;
92 3101664 : if (!ap || !(s = kross(D,a))) continue;
93 1895112 : nc[2] = (long)a;
94 1895112 : X = mseval2_ooQ(W, xpm, c); /* xpm(a / Dq) */
95 1895112 : x = umodiu(X, q); if (!x) continue;
96 : /* log(a) / log(c) */
97 1420890 : j = n? Fl_mul(logp(C, a, p, q, pn), ic, pn): 0;
98 1420890 : if (R) x = Fl_mul(x, tab[Fl_powu(ap, R, p)], q);
99 1420890 : if (s < 0) x = Fl_neg(x, q);
100 1420890 : v[j + 1] = Fl_add(v[j + 1], x, q);
101 : }
102 1038 : v = Flv_to_Flx(v, evalvarn(0));
103 1038 : v = zlx_translate1(v, p, N);
104 1038 : return gc_upto(av, Flx_to_ZX(v));
105 : }
106 : /* return lambda, assuming mu = 0 [else infinite loop] */
107 : static long
108 684 : lambda_ss(GEN W, GEN xpm, long v, long p, long D, long R, long n)
109 : {
110 324 : for (;; n += 2)
111 324 : {
112 684 : GEN P = polPn(W, xpm, p, D, R, n);
113 : long mu;
114 684 : if (!signe(P)) continue;
115 492 : mu = ZX_lvalrem(P, p, &P) + v;
116 492 : if (!mu)
117 : {
118 360 : long M = upowuu(p,n);
119 360 : if (odd(n)) M -= p; else M--;
120 360 : M /= p + 1;
121 360 : return Flx_val(ZX_to_Flx(P, p)) - M;
122 : }
123 : }
124 : }
125 : /* return lambda, assuming mu = 0 [else infinite loop] */
126 : static long
127 126 : lambda_ord(GEN W, GEN xpm, long v, long p, long D, long R, GEN ap)
128 : {
129 126 : GEN P1, P0 = polPn(W, xpm, p, D, R, 0);
130 : long n;
131 126 : for (n = 1;; n++, P0 = P1)
132 102 : {
133 228 : long mu, pn = upowuu(p,n);
134 228 : GEN P, xi, alpha, Q = utoipos(pn);
135 :
136 228 : P1 = polPn(W, xpm, p, D, R, n);
137 228 : alpha = mspadic_unit_eigenvalue(ap, 2, utoipos(p), n+1);
138 228 : alpha = padic_to_Q(ginv(alpha));
139 228 : xi = FpX_Fp_translate(polcyclo(pn, 0), gen_1, Q);
140 228 : P = ZX_sub(P1, ZX_Z_mul(FpX_mul(P0, xi, Q), alpha)); /* mod p^n */
141 :
142 228 : if (!signe(P) || n + v <= 0) continue;
143 126 : mu = ZX_lvalrem(P, p, &P) + v;
144 126 : if (!mu) return Flx_val(ZX_to_Flx(P, p));
145 : }
146 : }
147 : GEN
148 324 : ellpadiclambdamu(GEN E, long p, long D, long R)
149 : {
150 324 : pari_sp av = avma;
151 324 : long vC, s, muadd = 0;
152 : GEN W, xpm, C, ap;
153 :
154 324 : if (!sisfundamental(D))
155 6 : pari_err_DOMAIN("ellpadiclambdamu", "isfundamental(D)","=", gen_0, stoi(D));
156 318 : s = D < 0? -1: 1;
157 318 : if (odd(R)) s = -s;
158 :
159 318 : ap = ellap(E, utoi(p));
160 318 : if (ell_get_type(E) != t_ELL_Q)
161 6 : pari_err_TYPE("ellpadiclambdamu", E);
162 312 : if (!umodiu(ap, p))
163 : {
164 186 : if (Z_lval(ellQ_get_N(E), p) >= 2)
165 6 : pari_err_IMPL("additive reduction in ellpadiclambdamu");
166 180 : ap = NULL; /* supersingular */
167 : }
168 306 : if (ap)
169 : { /* ordinary */
170 126 : GEN v = ellisomat(E, p, 1), vE = gel(v,1), M = gel(v,2);
171 126 : if (lg(M) != 2) /* some p-isogeny */
172 : {
173 60 : long i, imax = 0, l = lg(vE);
174 60 : GEN wmax = NULL, vw = cgetg(l, t_VEC);
175 210 : for (i = 1; i < l; i++)
176 : {
177 150 : GEN w, e = ellinit(gel(vE,i), gen_1, 0), em = ellminimalmodel(e, NULL);
178 150 : gel(vE,i) = em; obj_free(e);
179 150 : w = ellQtwist_bsdperiod(gel(vE,i), s);
180 150 : if (s < 0) w = gel(w,2);
181 150 : gel(vw,i) = w;
182 : /* w is a positive real number, either Omega^+ or Omega^- / I */
183 150 : if (!imax || gcmp(w, wmax) > 0) { imax = i; wmax = w; }
184 : }
185 60 : if (imax != 1) muadd = Z_lval(ground(gdiv(gel(vw,imax), gel(vw,1))), p);
186 210 : for (i = 1; i < l; i++) obj_free(gel(vE,i));
187 60 : E = gel(vE, imax);
188 : }
189 : }
190 :
191 306 : W = msfromell(E, s);
192 306 : xpm = gel(W,2); W = gel(W,1);
193 306 : xpm = Q_primitive_part(xpm, &C);
194 306 : vC = C? Q_pval(C, utoipos(p)): 0;
195 306 : if (p == 2) vC++; /* due to half-loop in polPn */
196 306 : if (vC > 0) pari_err_BUG("ellpadiclambdamu [mu > 0]");
197 :
198 306 : if (ap)
199 : { /* ordinary */
200 126 : long L = lambda_ord(W, xpm, vC, p, D, R, ap);
201 126 : set_avma(av); return mkvec2s(L, muadd);
202 : }
203 : else
204 : {
205 180 : long Lp = lambda_ss(W, xpm, vC, p, D, R, 0);
206 180 : long Lm = lambda_ss(W, xpm, vC, p, D, R, 1);
207 180 : set_avma(av); retmkvec2(mkvec2s(Lp, Lm), zerovec(2));
208 : }
209 : }
|
