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 1771 : divmodulo(ulong a, ulong b, ulong p, ulong q)
22 : {
23 1771 : long v = u_lvalrem(b, p, &b);
24 1771 : if (v) a /= upowuu(p,v);
25 : /* a/b is now a p-integer */
26 1771 : return Fl_div(a, b, q);
27 : }
28 : /* to compute log_p(a) mod q = p^n, p < 2^31 */
29 : static GEN
30 854 : initQplog(long p, ulong q, long n)
31 : {
32 : long i, nn, nt;
33 : ulong a, P;
34 : GEN C;
35 1316 : for(nn = n, nt = n + 1; nn >= p; nn /= p) nt++;
36 854 : if (p == 2)
37 : {
38 70 : P = q - 8;
39 378 : while (3 * (nt - 1) > u_lval(nt-1, p) + n) nt--;
40 : }
41 : else
42 : {
43 784 : P = q - p;
44 1750 : while (nt > u_lval(nt-1, p) + n) nt--;
45 : }
46 854 : C = cgetg(nt, t_VECSMALL);
47 : /* [ P^(k-1) / k, k = 1..nt ] / (p - 1) */
48 2625 : for(i = 1, a = 1; i < nt; i++, a = Fl_mul(a, P, q))
49 1771 : C[i] = divmodulo(a, i * (p-1), p, q);
50 854 : 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 1638574 : logp(GEN C, ulong a, ulong p, ulong q, ulong pn)
56 : {
57 1638574 : long i, b, z, n = lg(C)-1;
58 1638574 : a %= q;
59 : /* Euclidean quotient (a^2 = 1 mod 8, a^(p-1) = 1 mod p) */
60 1638574 : if (p == 2)
61 476 : b = Fl_sqr(a, q << 1) >> 3;
62 : else
63 1638098 : b = Fl_powu(a, p-1, q) / p;
64 1638574 : z = Fl_mul(b, C[n], pn);
65 5354790 : for (i = n-1; i > 0; i--) z = Fl_mul(z + C[i], b, pn);
66 1638574 : 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 1211 : polPn(GEN W, GEN xpm, long p, long D, long R, long n)
76 : {
77 1211 : pari_sp av = avma, av2;
78 1211 : long N = (p == 2)? n + 2: n + 1;
79 1211 : ulong aD = labs(D), pn = upowuu(p, n), q = upowuu(p, N), ic = 0;
80 1211 : GEN v, Dq = muluu(aD, q), nc = icopy(gen_1), c = mkfrac(nc, Dq);
81 1211 : GEN C = n? initQplog(p, q, N): NULL;
82 1211 : GEN tab = R? ZV_to_Flv(teichmullerinit(p, N), q): NULL;
83 1211 : long a, A = itou(shifti(Dq,-1));
84 :
85 1211 : if (n) ic = Fl_inv(logp(C, 1 + (p == 2? 4: p), p, q, pn), pn);
86 1211 : v = zero_zv(pn); av2 = avma;
87 3619819 : for (a = 1; a <= A; a++, set_avma(av2))
88 : {
89 : GEN X;
90 3618608 : long s, ap = a % p;
91 : ulong x, j;
92 3618608 : if (!ap || !(s = kross(D,a))) continue;
93 2210964 : nc[2] = (long)a;
94 2210964 : X = mseval2_ooQ(W, xpm, c); /* xpm(a / Dq) */
95 2210964 : x = umodiu(X, q); if (!x) continue;
96 : /* log(a) / log(c) */
97 1657705 : j = n? Fl_mul(logp(C, a, p, q, pn), ic, pn): 0;
98 1657705 : if (R) x = Fl_mul(x, tab[Fl_powu(ap, R, p)], q);
99 1657705 : if (s < 0) x = Fl_neg(x, q);
100 1657705 : v[j + 1] = Fl_add(v[j + 1], x, q);
101 : }
102 1211 : v = Flv_to_Flx(v, evalvarn(0));
103 1211 : v = zlx_translate1(v, p, N);
104 1211 : return gc_upto(av, Flx_to_ZX(v));
105 : }
106 : /* return lambda, assuming mu = 0 [else infinite loop] */
107 : static long
108 798 : lambda_ss(GEN W, GEN xpm, long v, long p, long D, long R, long n)
109 : {
110 378 : for (;; n += 2)
111 378 : {
112 798 : GEN P = polPn(W, xpm, p, D, R, n);
113 : long mu;
114 798 : if (!signe(P)) continue;
115 574 : mu = ZX_lvalrem(P, p, &P) + v;
116 574 : if (!mu)
117 : {
118 420 : long M = upowuu(p,n);
119 420 : if (odd(n)) M -= p; else M--;
120 420 : M /= p + 1;
121 420 : 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 147 : lambda_ord(GEN W, GEN xpm, long v, long p, long D, long R, GEN ap)
128 : {
129 147 : GEN P1, P0 = polPn(W, xpm, p, D, R, 0);
130 : long n;
131 147 : for (n = 1;; n++, P0 = P1)
132 119 : {
133 266 : long mu, pn = upowuu(p,n);
134 266 : GEN P, xi, alpha, Q = utoipos(pn);
135 :
136 266 : P1 = polPn(W, xpm, p, D, R, n);
137 266 : alpha = mspadic_unit_eigenvalue(ap, 2, utoipos(p), n+1);
138 266 : alpha = padic_to_Q(ginv(alpha));
139 266 : xi = FpX_Fp_translate(polcyclo(pn, 0), gen_1, Q);
140 266 : P = ZX_sub(P1, ZX_Z_mul(FpX_mul(P0, xi, Q), alpha)); /* mod p^n */
141 :
142 266 : if (!signe(P) || n + v <= 0) continue;
143 147 : mu = ZX_lvalrem(P, p, &P) + v;
144 147 : if (!mu) return Flx_val(ZX_to_Flx(P, p));
145 : }
146 : }
147 : GEN
148 378 : ellpadiclambdamu(GEN E, long p, long D, long R)
149 : {
150 378 : pari_sp av = avma;
151 378 : long vC, s, muadd = 0;
152 : GEN W, xpm, C, ap;
153 :
154 378 : if (!sisfundamental(D))
155 7 : pari_err_DOMAIN("ellpadiclambdamu", "isfundamental(D)","=", gen_0, stoi(D));
156 371 : s = D < 0? -1: 1;
157 371 : if (odd(R)) s = -s;
158 :
159 371 : ap = ellap(E, utoi(p));
160 371 : if (ell_get_type(E) != t_ELL_Q)
161 7 : pari_err_TYPE("ellpadiclambdamu", E);
162 364 : if (!umodiu(ap, p))
163 : {
164 217 : if (Z_lval(ellQ_get_N(E), p) >= 2)
165 7 : pari_err_IMPL("additive reduction in ellpadiclambdamu");
166 210 : ap = NULL; /* supersingular */
167 : }
168 357 : if (ap)
169 : { /* ordinary */
170 147 : GEN v = ellisomat(E, p, 1), vE = gel(v,1), M = gel(v,2);
171 147 : if (lg(M) != 2) /* some p-isogeny */
172 : {
173 70 : long i, imax = 0, l = lg(vE);
174 70 : GEN wmax = NULL, vw = cgetg(l, t_VEC);
175 245 : for (i = 1; i < l; i++)
176 : {
177 175 : GEN w, e = ellinit(gel(vE,i), gen_1, 0), em = ellminimalmodel(e, NULL);
178 175 : gel(vE,i) = em; obj_free(e);
179 175 : w = ellQtwist_bsdperiod(gel(vE,i), s);
180 175 : if (s < 0) w = gel(w,2);
181 175 : gel(vw,i) = w;
182 : /* w is a positive real number, either Omega^+ or Omega^- / I */
183 175 : if (!imax || gcmp(w, wmax) > 0) { imax = i; wmax = w; }
184 : }
185 70 : if (imax != 1) muadd = Z_lval(ground(gdiv(gel(vw,imax), gel(vw,1))), p);
186 245 : for (i = 1; i < l; i++) obj_free(gel(vE,i));
187 70 : E = gel(vE, imax);
188 : }
189 : }
190 :
191 357 : W = msfromell(E, s);
192 357 : xpm = gel(W,2); W = gel(W,1);
193 357 : xpm = Q_primitive_part(xpm, &C);
194 357 : vC = C? Q_pval(C, utoipos(p)): 0;
195 357 : if (p == 2) vC++; /* due to half-loop in polPn */
196 357 : if (vC > 0) pari_err_BUG("ellpadiclambdamu [mu > 0]");
197 :
198 357 : if (ap)
199 : { /* ordinary */
200 147 : long L = lambda_ord(W, xpm, vC, p, D, R, ap);
201 147 : set_avma(av); return mkvec2s(L, muadd);
202 : }
203 : else
204 : {
205 210 : long Lp = lambda_ss(W, xpm, vC, p, D, R, 0);
206 210 : long Lm = lambda_ss(W, xpm, vC, p, D, R, 1);
207 210 : set_avma(av); retmkvec2(mkvec2s(Lp, Lm), zerovec(2));
208 : }
209 : }
|
