Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is to exceed 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - basemath - ellpadic.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.14.0 lcov report (development 27775-aca467eab2) Lines: 588 611 96.2 %
Date: 2022-07-03 07:33:15 Functions: 47 47 100.0 %
Legend: Lines: hit not hit

          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,gel(h,2),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 = gel(a1,2), q = gel(a1,3), An, Bn, Rn;
      34         259 :   long pp = precp(a1), v = valp(a1), i;
      35         259 :   int pis2 = absequaliu(p,2);
      36         259 :   a1 = gel(a1,4);
      37         259 :   b1 = gel(b1,4);
      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 = gel(r,2);
     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 = gel(r,2);
     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         546 : ellQp_root(GEN E, long prec)
     219         546 : { 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         707 : Tate_prec(GEN T) { return padicprec_relative(gel(T,3)); }
     293             : GEN
     294         819 : ellQp_Tate_uniformization(GEN E, long prec)
     295         819 : { return obj_checkbuild_prec(E,Qp_TATE,&doellQp_Tate, &Tate_prec,prec); }
     296             : GEN
     297         210 : ellQp_u(GEN E, long prec)
     298         210 : { 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         154 : ellQp_q(GEN E, long prec)
     304         154 : { GEN T = ellQp_Tate_uniformization(E, prec); return gel(T,3); }
     305             : GEN
     306         133 : ellQp_ab(GEN E, long prec)
     307         133 : { 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)|evalvalp(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 = gerepileupto(av2, gsub(w, gdiv(num, den)));
     365        1778 :     nold = nnew;
     366             :   }
     367         483 :   return gerepilecopy(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 gerepilecopy(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 gerepilecopy(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 gerepileupto(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 gerepileupto(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 gerepilecopy(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((5+k)*k, (ulong(*)(void*))&hash_GEN,
     499             :                                       (int(*)(void*,void*))&gidentical, 1);
     500         224 :   GEN b2 = ell_get_b2(e), b4 = ell_get_b4(e), n = gel(xP,1), d = gel(xP,2);
     501         224 :   GEN b6 = ell_get_b6(e), b8 = ell_get_b8(e);
     502             :   GEN B4, B6, B8, T, g4;
     503         224 :   GEN d2 = Fp_sqr(d,N), d3 = Fp_mul(d2,d,N), d4 = Fp_sqr(d2,N);
     504         224 :   GEN n2 = Fp_sqr(n,N), n3 = Fp_mul(n2,n,N), n4 = Fp_sqr(n2,N);
     505         224 :   GEN nd = Fp_mul(n,d,N), n2d2 = Fp_sqr(nd,N);
     506         224 :   GEN b2nd = Fp_mul(b2,nd, N), b2n2d = Fp_mul(b2nd,n,N);
     507         224 :   GEN b6d3 = Fp_mul(b6,d3,N), g,gp1,gm1, C,D;
     508         224 :   B8 = addii5(muliu(n4,3), mulii(b2n2d,n), mulii(muliu(b4,3), n2d2),
     509             :               mulii(muliu(b6d3,3), n), mulii(b8,d4));
     510         224 :   B6 = addii4(muliu(n3,4), mulii(b2nd,n),
     511             :               shifti(mulii(b4,Fp_mul(n,d2,N)), 1), b6d3);
     512         224 :   B4 = addii3(muliu(n2,6), b2nd,  mulii(b4,d2));
     513         224 :   B4 = modii(B4,N);
     514         224 :   B6 = modii(B6,N);
     515         224 :   B8 = modii(B8,N);
     516         224 :   g4 = Fp_sub(sqri(B6), mulii(B4,B8), N);
     517         224 :   T = Fp_sqr(B6,N);
     518             : 
     519         224 :   g = rellg(H, m, T,g4,B8, N);
     520         224 :   gp1 = rellg(H, addiu(m,1), T,g4,B8, N);
     521         224 :   gm1 = rellg(H, subiu(m,1), T,g4,B8, N);
     522         224 :   C = Fp_sqr(g, N);
     523         224 :   D = Fp_mul(gp1,gm1, N);
     524         224 :   if (mpodd(m))
     525             :   {
     526         112 :     n = Fp_sub(mulii(C,n), mulii(D,B6), N);
     527         112 :     d = Fp_mul(C,d, N);
     528             :   }
     529             :   else
     530             :   {
     531         112 :     n = Fp_sub(Fp_mul(Fp_mul(B6,C,N), n, N), D, N);
     532         112 :     d = Fp_mul(Fp_mul(C,d,N), B6, N);
     533             :   }
     534         224 :   return gerepilecopy(av, mkvec2(n,d));
     535             : }
     536             : /* given [n,d2], x = n/d2 (coprime, d2 = d^2), p | den,
     537             :  * return t = -x/y + O(p^v) */
     538             : static GEN
     539         126 : tfromx(GEN e, GEN x, GEN p, long v, GEN N, GEN *pd)
     540             : {
     541         126 :   GEN n = gel(x,1), d2 = gel(x,2), d;
     542             :   GEN a1, a3, b2, b4, b6, B, C, d4, d6, Y;
     543         126 :   if (!signe(n)) { *pd = gen_1; return zeropadic(p, v); }
     544         126 :   a1 = ell_get_a1(e);
     545         126 :   b2 = ell_get_b2(e);
     546         126 :   a3 = ell_get_a3(e);
     547         126 :   b4 = ell_get_b4(e);
     548         126 :   b6 = ell_get_b6(e);
     549         126 :   d = Qp_sqrt(cvtop(d2, p, v - Z_pval(d2,p)));
     550         126 :   if (!d) pari_err_BUG("ellpadicheight");
     551             :   /* Solve Y^2 = 4n^3 + b2 n^2 d2+ 2b4 n d2^2 + b6 d2^3,
     552             :    * Y = 2y + a1 n d + a3 d^3 */
     553         126 :   d4 = Fp_sqr(d2, N);
     554         126 :   d6 = Fp_mul(d4, d2, N);
     555         126 :   B = gmul(d, Fp_add(mulii(a1,n), mulii(a3,d2), N));
     556         126 :   C = mkpoln(4, utoipos(4), Fp_mul(b2, d2, N),
     557             :                 Fp_mul(shifti(b4,1), d4, N),
     558             :                 Fp_mul(b6,d6,N));
     559         126 :   C = FpX_eval(C, n, N);
     560         126 :   if (!signe(C))
     561           0 :     Y = zeropadic(p, v >> 1);
     562             :   else
     563         126 :     Y = Qp_sqrt(cvtop(C, p, v - Z_pval(C,p)));
     564         126 :   if (!Y) pari_err_BUG("ellpadicheight");
     565         126 :   *pd = d;
     566         126 :   return gdiv(gmulgs(gmul(n,d), -2), gsub(Y,B));
     567             : }
     568             : 
     569             : /* return minimal i s.t. -v_p(j+1) - log_p(j-1) + (j+1)*t >= v for all j>=i */
     570             : static long
     571         126 : logsigma_prec(GEN p, long v, long t)
     572             : {
     573         126 :   double log2p = dbllog2(p);
     574         126 :   long j, i = ceil((v - t) / (t - 2*M_LN2/(3*log2p)) + 0.01);
     575         126 :   if (absequaliu(p,2) && i < 5) i = 5;
     576             :   /* guaranteed to work, now optimize */
     577         182 :   for (j = i-1; j >= 2; j--)
     578             :   {
     579         182 :     if (- u_pval(j+1,p) - log2(j-1)/log2p + (j+1)*t + 0.01 < v) break;
     580          56 :     i = j;
     581             :   }
     582         126 :   if (j == 1)
     583             :   {
     584           0 :     if (- absequaliu(p,2) + 2*t + 0.01 >= v) i = 1;
     585             :   }
     586         126 :   return i;
     587             : }
     588             : /* return minimal i s.t. -v_p(j+1) + (j+1)*t >= v for all j>=i */
     589             : static long
     590           7 : log_prec(GEN p, long v, long t)
     591             : {
     592           7 :   double log2p = dbllog2(p);
     593           7 :   long j, i = ceil(v / (t - M_LN2/(2*log2p)) + 0.01);
     594             :   /* guaranteed to work, now optimize */
     595          28 :   for (j = i-1; j >= 1; j--)
     596             :   {
     597          28 :     if (- u_pval(j+1,p) + (j+1)*t + 0.01 < v) break;
     598          21 :     i = j;
     599             :   }
     600           7 :   return i;
     601             : }
     602             : 
     603             : /* P = rational point of exact denominator d. Is Q singular on E(Fp) ? */
     604             : static int
     605         182 : FpE_issingular(GEN E, GEN P, GEN d, GEN p)
     606             : {
     607         182 :   pari_sp av = avma;
     608             :   GEN t, x, y, a1, a2, a3, a4;
     609         182 :   if (ell_is_inf(E) || dvdii(d,p)) return 0; /* 0_E is smooth */
     610         175 :   P = Q_muli_to_int(P,d);
     611         175 :   x = gel(P,1);
     612         175 :   y = gel(P,2);
     613         175 :   a1 = ell_get_a1(E);
     614         175 :   a3 = ell_get_a3(E);
     615         175 :   t = addii(shifti(y,1), addii(mulii(a1,x), mulii(a3,d)));
     616         175 :   if (!dvdii(t,p)) return gc_bool(av, 0);
     617          35 :   a2 = ell_get_a2(E);
     618          35 :   a4 = ell_get_a4(E);
     619          35 :   d = Fp_inv(d, p);
     620          35 :   x = Fp_mul(x,d,p);
     621          35 :   y = Fp_mul(y,d,p);
     622          35 :   t = subii(mulii(a1,y), addii(a4, mulii(x, addii(gmul2n(a2,1), muliu(x,3)))));
     623          35 :   return gc_bool(av, dvdii(t,p));
     624             : }
     625             : /* E/Q, P on E(Q). Let g > 0 minimal such that the image of R = [g]P in a
     626             :  * minimal model is everywhere nonsingular; return [R,g] */
     627             : GEN
     628         140 : ellnonsingularmultiple(GEN e, GEN P)
     629             : {
     630         140 :   pari_sp av = avma;
     631         140 :   GEN ch, E = ellanal_globalred(e, &ch), NP, L, S, d, g = gen_1;
     632             :   long i, l;
     633         140 :   checkellpt(P);
     634         140 :   if (ell_is_inf(P)) retmkvec2(gcopy(P), gen_1);
     635         140 :   if (E != e) P = ellchangepoint(P, ch);
     636         140 :   S = obj_check(E, Q_GLOBALRED);
     637         140 :   NP = gmael(S,3,1);
     638         140 :   L = gel(S,4);
     639         140 :   l = lg(NP);
     640         140 :   d = Q_denom(P);
     641         315 :   for (i = 1; i < l; i++)
     642             :   {
     643         175 :     GEN G = gel(L,i), p = gel(NP,i);/* prime of bad reduction */
     644             :     long kod;
     645         175 :     if (!FpE_issingular(E, P, d, p)) continue;
     646          21 :     kod = itos(gel(G,2)); /* Kodaira type */
     647          21 :     if (kod >= 5) /* I_nu */
     648             :     {
     649           7 :       long nu = kod - 4;
     650           7 :       long n = minss(Q_pval(ec_dmFdy_evalQ(E,P), p), nu/2);
     651           7 :       nu /= ugcd(nu, n);
     652           7 :       g = muliu(g, nu);
     653           7 :       P = ellmul(E, P, utoipos(nu)); d = Q_denom(P);
     654          14 :     } else if (kod <= -5) /* I^*_nu */
     655             :     { /* either 2 or 4 */
     656           7 :       long nu = - kod - 4;
     657           7 :       P = elladd(E, P,P); d = Q_denom(P);
     658           7 :       g = shifti(g,1);
     659           7 :       if (odd(nu) && FpE_issingular(E, P, d, p))
     660             :       { /* it's 4 */
     661           7 :         P = elladd(E, P,P); d = Q_denom(P);
     662           7 :         g = shifti(g,1);
     663             :       }
     664             :     } else {
     665           7 :       GEN c = gel(G, 4); /* Tamagawa number at p */
     666           7 :       if (absequaliu(c, 4)) c = gen_2;
     667           7 :       P = ellmul(E, P, c); d = Q_denom(P);
     668           7 :       g = mulii(g, c);
     669             :     }
     670             :   }
     671         140 :   if (E != e) P = ellchangepointinv(P, ch);
     672         140 :   return gerepilecopy(av, mkvec2(P,g));
     673             : }
     674             : 
     675             : GEN
     676         140 : ellpadicheight(GEN e, GEN p, long v0, GEN P)
     677             : {
     678         140 :   pari_sp av = avma;
     679             :   GEN N, H, h, t, ch, g, E, x, D, ls, lt, S, a,b;
     680             :   long v, vd;
     681             :   int is2;
     682         140 :   checkellpt(P);
     683         140 :   if (v0<=0) pari_err_DOMAIN("ellpadicheight","precision","<=",gen_0,stoi(v0));
     684         133 :   checkell_Q(e);
     685         133 :   if (typ(p) != t_INT) pari_err_TYPE("ellpadicheight",p);
     686         133 :   if (cmpiu(p,2) < 0) pari_err_PRIME("ellpadicheight",p);
     687         133 :   if (ellorder_Q(e,P)) return mkvec2(gen_0,gen_0);
     688         126 :   E = ellanal_globalred(e, &ch);
     689         126 :   if (E != e) P = ellchangepoint(P, ch);
     690         126 :   S = ellnonsingularmultiple(E, P);
     691         126 :   P = gel(S,1);
     692         126 :   g = gel(S,2);
     693         126 :   v = v0 + 2*Z_pval(g, p);
     694         126 :   is2 = absequaliu(p,2); if (is2) v += 2;
     695         126 :   x = Q_remove_denom(gel(P,1), &b);
     696         126 :   x = mkvec2(x, b? b: gen_1);
     697         126 :   vd = Z_pval(gel(x,2), p);
     698         126 :   if (!vd)
     699             :   { /* P not in kernel of reduction mod p */
     700         112 :     GEN d, m, X, Pp, Ep = ellinit(E, mkintmod(gen_0,p), 0);
     701         112 :     long w = v+2;
     702         112 :     Pp = RgV_to_FpV(P, p);
     703         112 :     if (lg(Ep) != 1)
     704         105 :       m = ellorder(Ep, Pp, NULL);
     705             :     else
     706             :     { /* E has bad reduction at p */
     707           7 :       m = ellcard(E, p);
     708           7 :       if (equalii(m, p)) pari_err_TYPE("ellpadicheight: additive reduction", E);
     709             :     }
     710         112 :     g = mulii(g,m);
     711             :     for(;;)
     712             :     {
     713         224 :       N = powiu(p, w);
     714         224 :       X = xmP(E, x, m, N);
     715         224 :       d = gel(X,2);
     716         224 :       if (!signe(d)) w <<= 1;
     717             :       else
     718             :       {
     719         224 :         vd = Z_pval(d, p);
     720         224 :         if (w >= v+2*vd + is2) break;
     721         112 :         w = v+2*vd + is2;
     722             :       }
     723             :     }
     724         112 :     x = X;
     725             :   }
     726             :   /* we will want t mod p^(v+vd) because of t/D in H later, and
     727             :    * we lose p^vd in tfromx because of sqrt(d) (p^(vd+1) if p=2)*/
     728         126 :   v += 2*vd + is2;
     729         126 :   N = powiu(p,v);
     730         126 :   t = tfromx(E, x, p, v, N, &D); /* D^2=denom(x)=x[2] */
     731         126 :   S = ellformallogsigma_t(E, logsigma_prec(p, v-vd, valp(t)) + 1);
     732         126 :   ls = ser2rfrac_i(gel(S,1)); /* log_p (sigma(T)/T) */
     733         126 :   lt = ser2rfrac_i(gel(S,2)); /* log_E (T) */
     734             :   /* evaluate our formal power series at t */
     735         126 :   H = gadd(poleval(ls, t), glog(gdiv(t, D), 0));
     736         126 :   h = gsqr(poleval(lt, t));
     737         126 :   g = sqri(g);
     738         126 :   a = gdiv(gmulgs(H,-2), g);
     739         126 :   b = gdiv(gneg(h), g);
     740         126 :   if (E != e)
     741             :   {
     742          21 :     GEN u = gel(ch,1), r = gel(ch,2);
     743          21 :     a = gdiv(gadd(a, gmul(r,b)), u);
     744          21 :     b = gmul(u,b);
     745             :   }
     746         126 :   H = mkvec2(a,b);
     747         126 :   gel(H,1) = precp_fix(gel(H,1),v0);
     748         126 :   gel(H,2) = precp_fix(gel(H,2),v0);
     749         126 :   return gerepilecopy(av, H);
     750             : }
     751             : 
     752             : GEN
     753          14 : ellpadiclog(GEN E, GEN p, long n, GEN P)
     754             : {
     755          14 :   pari_sp av = avma;
     756             :   long vt;
     757             :   GEN t, x, y, L;
     758          14 :   checkellpt(P);
     759          14 :   if (ell_is_inf(P)) return gen_0;
     760          14 :   x = gel(P,1);
     761          14 :   y = gel(P,2); t = gneg(gdiv(x,y));
     762          14 :   vt = gvaluation(t, p); /* can be a t_INT, t_FRAC or t_PADIC */
     763          14 :   if (vt <= 0)
     764           7 :     pari_err_DOMAIN("ellpadiclog","P","not in the kernel of reduction at",p,P);
     765           7 :   L = ser2rfrac_i(ellformallog(E, log_prec(p, n, vt) + 1, 0));
     766           7 :   return gerepileupto(av, poleval(L, cvtop(t, p, n)));
     767             : }
     768             : 
     769             : /* E/Qp has multiplicative reduction, Tate curve */
     770             : static GEN
     771          14 : ellpadics2_tate(GEN Ep, long n)
     772             : {
     773             :   pari_sp av;
     774          14 :   GEN u2 = ellQp_u2(Ep, n), q = ellQp_q(Ep, n), pn = gel(q,3);
     775          14 :   GEN qm, s, b2 = ell_get_b2(Ep), v = vecfactoru_i(1, n);
     776             :   long m;
     777          14 :   qm = Fp_powers(padic_to_Fp(q, pn), n, pn);
     778          14 :   s = gel(qm, 2); av = avma;
     779          70 :   for (m = 2; m <= n; m++) /* sum sigma(m) q^m */
     780             :   {
     781          56 :     s = addii(s, mulii(gel(qm,m+1), usumdiv_fact(gel(v,m))));
     782          56 :     if ((m & 31) == 0) s = gerepileuptoint(av, s);
     783             :   }
     784          14 :   s = subui(1, muliu(s,24));
     785          14 :   s = gdivgu(gsub(b2, gdiv(s,u2)), 12);
     786          14 :   return precp_fix(s,n);
     787             : }
     788             : 
     789             : GEN
     790         189 : ellpadicfrobenius(GEN E, ulong p, long n)
     791             : {
     792         189 :   checkell(E);
     793         189 :   if (p < 2) pari_err_DOMAIN("ellpadicfrobenius","p", "<", gen_2, utoipos(p));
     794         189 :   switch(ell_get_type(E))
     795             :   {
     796         175 :     case t_ELL_Q: break;
     797          14 :     case t_ELL_Qp: if (equaliu(ellQp_get_p(E), p)) break;
     798           7 :     default: pari_err_TYPE("ellpadicfrobenius", E);
     799             :   }
     800         182 :   return hyperellpadicfrobenius(ec_bmodel(E), p, n);
     801             : }
     802             : 
     803             : /* s2 = (b_2-E_2)/12 */
     804             : GEN
     805          35 : ellpadics2(GEN E, GEN p, long n)
     806             : {
     807          35 :   pari_sp av = avma;
     808             :   GEN l, F, a,b,d, ap;
     809             :   ulong pp;
     810          35 :   if (typ(p) != t_INT) pari_err_TYPE("ellpadics2",p);
     811          35 :   if (cmpis(p,2) < 0) pari_err_PRIME("ellpadics2",p);
     812          35 :   checkell(E);
     813          35 :   if (Q_pval(ell_get_j(E), p) < 0)
     814             :   {
     815             :     GEN Ep;
     816           7 :     if (ell_get_type(E) == t_ELL_Qp) Ep = E;
     817           0 :     else Ep = ellinit(E, zeropadic(p,n), 0);
     818           7 :     l = ellpadics2_tate(Ep, n);
     819           7 :     if (Ep != E) obj_free(Ep);
     820           7 :     return gerepilecopy(av, l);
     821             :   }
     822          28 :   pp = itou(p);
     823          28 :   F = ellpadicfrobenius(E, pp, n);
     824          21 :   a = gcoeff(F,1,1);
     825          21 :   b = gcoeff(F,1,2);
     826          21 :   d = gcoeff(F,2,2); ap = gadd(a,d);
     827          21 :   if (valp(ap) > 0) pari_err_DOMAIN("ellpadics2","E","is supersingular at",p,E);
     828          21 :   if (pp == 2 || (pp <= 13 && n == 1)) /* 2sqrt(p) > p/2: ambiguity */
     829           0 :     ap = ellap(E,p);
     830             :   else
     831             :   { /* either 2sqrt(p) < p/2 or n > 1 and 2sqrt(p) < p^2/2 (since p!=2) */
     832          21 :     GEN q = pp <= 13? utoipos(pp * pp): p;
     833          21 :     ap = padic_to_Fp(ap, q);
     834          21 :     ap = Fp_center_i(ap, q, shifti(q,-1));
     835             :   }
     836          21 :   l = mspadic_unit_eigenvalue(ap, 2, p, n);
     837          21 :   return gerepileupto(av, gdiv(b, gsub(l, a))); /* slope of eigenvector */
     838             : }
     839             : 
     840             : /* symbol and modular symbol space attached to E to later compute
     841             :  * ellpadicL(E,p,, s,,D) */
     842             : static GEN
     843          91 : ellpadicL_symbol(GEN E, GEN p, GEN s, GEN D)
     844             : {
     845             :   GEN s1, s2, ap;
     846             :   long sign;
     847          91 :   checkell(E);
     848          91 :   if (ell_get_type(E) != t_ELL_Q) pari_err_TYPE("ellpadicL",E);
     849          91 :   ap = ellap(E,p);
     850          91 :   if (D && typ(D) != t_INT) pari_err_TYPE("ellpadicL",D);
     851          91 :   if (D && !Z_isfundamental(D))
     852           0 :     pari_err_DOMAIN("ellpadicL", "isfundamental(D)", "=", gen_0, D);
     853          91 :   if (!D) D = gen_1;
     854          91 :   if (Z_pval(ellQ_get_N(E), p) >= 2)
     855           0 :     pari_err_IMPL("additive reduction in ellpadicL");
     856          91 :   mspadic_parse_chi(s, &s1,&s2);
     857          91 :   sign = signe(D); if (mpodd(s2)) sign = -sign;
     858          91 :   return shallowconcat(msfromell(E, sign), mkvec4(ap, p, s, D));
     859             : }
     860             : /* W an ellpadicL_symbol, initialize for ellpadicL(E,p,n,s,,D) */
     861             : static GEN
     862          91 : ellpadicL_init(GEN W, long n)
     863             : {
     864          91 :   GEN Wp, den, M = gel(W,1), xpm = gel(W,2), ap = gel(W,3), s = gel(W,5);
     865          91 :   long p = itos(gel(W,4)), D = itos(gel(W,6));
     866             : 
     867          91 :   xpm = Q_remove_denom(xpm,&den); if (!den) den = gen_1;
     868          91 :   n += Z_lval(den,p);
     869          91 :   Wp = mspadicinit(M, p, n, Z_lval(ap,p));
     870          91 :   return mkvec3(mspadicmoments(Wp,xpm,D), den, s);
     871             : }
     872             : /* v from ellpadicL_init, compute ellpadicL(E,p,n,s,r,D) */
     873             : static GEN
     874         126 : ellpadic_i(GEN v, long r)
     875             : {
     876         126 :   GEN oms = gel(v,1), den = gel(v,2), s = gel(v,3);
     877         126 :   return gdiv(mspadicL(oms,s,r), den);
     878             : }
     879             : GEN
     880          63 : ellpadicL(GEN E, GEN p, long n, GEN s, long r, GEN D)
     881             : {
     882          63 :   pari_sp av = avma;
     883             :   GEN W, v;
     884          63 :   if (r < 0) pari_err_DOMAIN("ellpadicL","r","<",gen_0,stoi(r));
     885          63 :   if (n <= 0) pari_err_DOMAIN("ellpadicL","precision","<=",gen_0,stoi(n));
     886          63 :   W = ellpadicL_symbol(E, p, s, D);
     887          63 :   v = ellpadicL_init(W, n);
     888          63 :   return gerepileupto(av, ellpadic_i(v, r));
     889             : }
     890             : 
     891             : static long
     892          28 : torsion_order(GEN E) { GEN T = elltors(E); return itos(gel(T,1)); }
     893             : 
     894             : /* E given by a minimal model; D != 0. Compare Euler factor of L(E,(D/.),1)
     895             :  * with L(E^D,1). Return
     896             :  *   \prod_{p|D} (p-a_p(E)+eps_{E}(p)) / p,
     897             :  * where eps(p) = 0 if p | N_E and 1 otherwise */
     898             : static GEN
     899           7 : get_Euler(GEN E, GEN D)
     900             : {
     901           7 :   GEN a = gen_1, b = gen_1, P = gel(absZ_factor(D), 1);
     902           7 :   long i, l = lg(P);
     903          14 :   for (i = 1; i < l; i++)
     904             :   {
     905           7 :     GEN p = gel(P,i);
     906           7 :     a = mulii(a, ellcard(E, p));
     907           7 :     b = mulii(b, p);
     908             :   }
     909           7 :   return Qdivii(a, b);
     910             : }
     911             : GEN
     912          28 : ellpadicbsd(GEN E, GEN p, long n, GEN D)
     913             : {
     914          28 :   const long MAXR = 30;
     915          28 :   pari_sp av = avma;
     916          28 :   GEN D0, W, ED, tam, ND, C, apD, U = NULL;/*-Wall*/
     917             :   long r, vN;
     918          28 :   checkell(E);
     919          28 :   if (D && isint1(D)) D = NULL;
     920          28 :   D0 = ellminimaltwistcond(E); if (isint1(D0)) D0 = NULL;
     921          28 :   if (D0) { E = elltwist(E, D0); D = D? coredisc(mulii(D, D0)): D0; }
     922          28 :   W = ellpadicL_symbol(E, p, gen_0, D);
     923          28 :   ED = D? elltwist(E, D): E;
     924          28 :   ED = ellanal_globalred_all(ED, NULL, &ND, &tam);
     925          28 :   vN = Z_pval(ND, p); /* additive reduction ? */
     926          28 :   if (vN >= 2) pari_err_DOMAIN("ellpadicbsd","v_p(N(E_D))",">",gen_1,stoi(vN));
     927          28 :   if (n < 5) n = 5;
     928           0 :   for(;; n <<= 1)
     929           0 :   {
     930          28 :     pari_sp av2 = avma;
     931          28 :     GEN v = ellpadicL_init(W, n);
     932          63 :     for (r = 0; r < MAXR; r++)
     933             :     {
     934          63 :       U = ellpadic_i(v, r);
     935          63 :       if (!gequal0(U)) break;
     936             :     }
     937          28 :     if (r < MAXR) break;
     938           0 :     set_avma(av2);
     939             :   }
     940          28 :   apD = ellap(ED, p);
     941          28 :   if (typ(U) == t_COL)
     942             :   { /* p | a_p(E_D), frobenius on E_D */
     943           0 :     GEN F = mkmat22(gen_0, negi(p), gen_1, apD);
     944           0 :     U = RgM_RgC_mul(gpowgs(gsubsg(1, gdiv(F,p)), -2), U);
     945           0 :     settyp(U, t_VEC);
     946             :   }
     947          28 :   else if (dvdii(ND,p))
     948             :   {
     949          21 :     if (equalim1(apD)) /* divide by 1-1/a_p */
     950          14 :       U = gdivgu(U, 2);
     951             :     else
     952             :     { /* ap = 1 */
     953           7 :       GEN EDp = ellinit(ED, zeropadic(p,n), 0);
     954           7 :       U = gdiv(U, ellQp_L(EDp,n));
     955           7 :       obj_free(EDp);
     956             :     }
     957             :   }
     958             :   else
     959             :   {
     960           7 :     GEN a = mspadic_unit_eigenvalue(apD, 2, p, n);
     961           7 :     U = gmul(U, gpowgs(gsubsg(1, ginv(a)), -2));
     962             :   }
     963          28 :   C = mulii(tam, mpfact(r));
     964          28 :   if (D) C = gmul(C, get_Euler(ED, D));
     965          28 :   C = gdiv(sqru(torsion_order(ED)), C);
     966          28 :   if (D) obj_free(ED);
     967          28 :   return gerepilecopy(av, mkvec2(utoi(r), gmul(U, C)));
     968             : }
     969             : 
     970             : GEN
     971           7 : ellpadicregulator(GEN E, GEN p, long n, GEN S)
     972             : {
     973           7 :   pari_sp av = avma;
     974           7 :   GEN FG = ellpadicheightmatrix(E,p,n,S); /* forbids additive reduction */
     975             :   /* [F,G]: height in basis [omega, eta] */
     976           7 :   GEN R, F = gel(FG,1), G = gel(FG,2), ap = ellap(E,p);
     977           7 :   if (dvdii(ap, p))
     978             :   { /* supersingular */
     979           0 :     GEN f = ellpadicfrobenius(E, itou(p), n);
     980           0 :     GEN x = gcoeff(f,1,1), y = gcoeff(f,2,1);
     981             :     /* [A,B]: regulator height in basis [omega, eta] */
     982           0 :     GEN A = det(F), B = det(gadd(F,G)), C;
     983           0 :     C = gdiv(gsub(B,A), y);
     984             :     /* R: regulator height in basis [omega, F.omega] */
     985           0 :     R = mkvec2(gsub(A, gmul(x,C)), C);
     986             :   }
     987             :   else
     988             :   {
     989             :     GEN s2;
     990           7 :     if (equali1(ap) && dvdii(ell_get_disc(E),p))
     991           7 :     { /* split multiplicative reduction */
     992           7 :       GEN Ep = ellinit(E, zeropadic(p,n), 0);
     993           7 :       GEN q = ellQp_q(Ep,n), u2 = ellQp_u2(Ep,n);
     994           7 :       s2 = ellpadics2_tate(Ep, n);
     995           7 :       s2 = gsub(s2, ginv(gmul(Qp_log(q), u2))); /*extended MW group contrib*/
     996           7 :       obj_free(Ep);
     997             :     }
     998             :     else
     999           0 :       s2 = ellpadics2(E,p,n);
    1000           7 :     R = det( RgM_sub(F, RgM_Rg_mul(G,s2)) );
    1001             :   }
    1002           7 :   return gerepilecopy(av, R);
    1003             : }

Generated by: LCOV version 1.13