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 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.10.0 lcov report (development 21682-493a494) Lines: 581 605 96.0 %
Date: 2018-01-16 06:18:33 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. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
      13             : 
      14             : #include "pari.h"
      15             : #include "paripriv.h"
      16             : 
      17             : static GEN
      18         266 : precp_fix(GEN h, long v)
      19         266 : { return (precp(h) > v)? cvtop(h,gel(h,2),v): h; }
      20             : 
      21             : /* TATE CURVE */
      22             : 
      23             : /* a1, b1 are t_PADICs, a1/b1 = 1 (mod p) if p odd, (mod 2^4) otherwise.
      24             :  * Let (A_n, B_n) be defined by A_1 = a1/p^v, B_1 = b1/p^v, v=v(a1)=v(a2);
      25             :  *   A_{n+1} = (A_n + B_n + 2 B_{n+1}) / 4
      26             :  *   B_{n+1} = B_n sqrt(A_n / B_n) = square root of A_n B_n congruent to B_n
      27             :  *   R_n = p^v( A_n - B_n ) = r_{n+1}
      28             :  * Return [An,Bn,Rn]. N.B. lim An = M2(a1,b1) = M(sqrt(a1),sqrt(b1))^2 */
      29             : GEN
      30         245 : Qp_agm2_sequence(GEN a1, GEN b1)
      31             : {
      32         245 :   GEN bp, pmod, p = gel(a1,2), q = gel(a1,3), An, Bn, Rn;
      33         245 :   long pp = precp(a1), v = valp(a1), i;
      34         245 :   int pis2 = absequaliu(p,2);
      35         245 :   a1 = gel(a1,4);
      36         245 :   b1 = gel(b1,4);
      37         245 :   if (pis2)
      38         154 :     pmod = utoipos(8);
      39             :   else
      40          91 :     pmod = p;
      41         245 :   bp = modii(b1, pmod);
      42         245 :   An = cgetg(pp+1, t_VEC); /* overestimate: rather log_2(pp) */
      43         245 :   Bn = cgetg(pp+1, t_VEC);
      44         245 :   Rn = cgetg(pp+1, t_VEC);
      45         854 :   for(i = 1;; i++)
      46             :   {
      47         854 :     GEN a = a1, b = b1, r;
      48             :     long vr;
      49         854 :     gel(An, i) = a;
      50         854 :     gel(Bn, i) = b;
      51         854 :     r = subii(a,b);
      52        1099 :     if (!signe(r)) break;
      53         637 :     vr = Z_pvalrem(r,p,&r);
      54         637 :     if (vr >= pp) break;
      55         609 :     r = cvtop(r, p, pp - vr); setvalp(r, vr+v);
      56         609 :     gel(Rn, i) = r;
      57             : 
      58         609 :     b1 = Zp_sqrt(Fp_mul(a,b,q), p, pp);
      59         609 :     if (!b1) pari_err_PREC("p-adic AGM");
      60         609 :     if (!equalii(modii(b1,pmod), bp)) b1 = Fp_neg(b1, q);
      61             :     /* a1 = (a+b+2sqrt(ab))/4 */
      62         609 :     if (pis2)
      63             :     {
      64         343 :       b1 = remi2n(b1, pp-1);
      65         343 :       a1 = shifti(addii(addii(a,b), shifti(b1,1)),-2);
      66         343 :       a1 = remi2n(a1, pp-2);
      67         343 :       pp -= 2;
      68             :     }
      69             :     else
      70         266 :       a1 = modii(Fp_halve(addii(Fp_halve(addii(a,b),q), b1), q), q);
      71         609 :   }
      72         245 :   setlg(An,i+1);
      73         245 :   setlg(Bn,i+1);
      74         245 :   setlg(Rn,i); return mkvec4(An, Bn, Rn, stoi(v));
      75             : }
      76             : void
      77         336 : Qp_descending_Landen(GEN AB, GEN *ptx, GEN *pty)
      78             : {
      79         336 :   GEN R = gel(AB,3);
      80         336 :   long i, n = lg(R)-1;
      81         336 :   GEN x = *ptx;
      82         336 :   if (isintzero(x))
      83             :   {
      84         245 :     i = 2;
      85         245 :     x = gmul2n(gel(R,1),-2);
      86         245 :     if (pty)
      87             :     {
      88           0 :       GEN A = gel(AB,1);
      89           0 :       if (n == 1)
      90           0 :         *pty = gmul(x, Qp_sqrt(gadd(x,gel(A,2))));
      91             :       else
      92           0 :         *pty = Qp_sqrt(gmul(gmul(x, gadd(x,gel(A,2))), gadd(x,gel(R,2))));
      93           0 :       if (!*pty) pari_err_PREC("Qp_descending_Landen");
      94             :     }
      95             :   }
      96             :   else
      97          91 :     i = 1;
      98         938 :   for (; i <= n; i++)
      99             :   {
     100         602 :     GEN r = gel(R,i), t;
     101         602 :     if (gequal0(x)) pari_err_PREC("Qp_descending_Landen");
     102         602 :     t = Qp_sqrt(gaddsg(1, gdiv(r,x))); /* = 1 (mod p) */
     103         602 :     if (!t) pari_err_PREC("Qp_descending_Landen");
     104         602 :     if (i == n)
     105             :     {
     106         287 :       GEN p = gel(r,2);
     107         287 :       long v, vx = valp(x), vr = valp(r);
     108         287 :       if (vx >= vr) pari_err_PREC("Qp_descending_Landen");
     109             :       /* last loop, take into account loss of accuracy from multiplication
     110             :        * by \prod_{j > n} sqrt(1+r_j/x_j); since vx < vr, j = n+1 is enough */
     111         287 :       v = 2*vr - vx;
     112             :       /* |r_{n+1}| <= |(r_n)^2 / 8| + 1 bit for sqrt loss */
     113         287 :       if (absequaliu(p,2)) v -= 4;
     114             :       /* tail is 1 + O(p^v) */
     115         287 :       if (v < precp(x)) x = cvtop(x,p,v);
     116             :     }
     117             :     /* x_{n+1} = x_n  ((1 + sqrt(1 + r_n/x_n)) / 2)^2 */
     118         602 :     x = gmul(x, gsqr(gmul2n(gaddsg(1,t),-1)));
     119             :     /* y_{n+1} = y_n / (1 - (r_n/4x_{n+1})^2) */
     120         602 :     if (pty) *pty = gdiv(*pty, gsubsg(1, gsqr(gdiv(r,gmul2n(x,2)))));
     121             :   }
     122         336 :   *ptx = x;
     123         336 : }
     124             : void
     125          56 : Qp_ascending_Landen(GEN AB, GEN *ptx, GEN *pty)
     126             : {
     127          56 :   GEN A = gel(AB,1), R = gel(AB,3), x = *ptx, p, r;
     128          56 :   long n = lg(R)-1, va = itos(gel(AB,4)), v, i;
     129             : 
     130          56 :   r = gel(R,n);
     131          56 :   v = 2*valp(r) + va;
     132          56 :   if (typ(x) == t_PADIC)
     133          35 :     v -= 2*valp(x);
     134             :   else
     135          21 :     v -= valp(gnorm(x)); /* v(x) = v(Nx) / (e*f), here ef = 2 */
     136          56 :   p = gel(r,2);
     137          56 :   if (absequaliu(p,2)) v -= 3; /* |r_{n+1}| <= |(r_n)^2 / 8| */
     138             :   /* v = v(A[n+1] R[n+1] / x_{n+1}^2) */
     139          56 :   if (v <= 0) pari_err_PREC("Qp_ascending_Landen");
     140             :   /* v > 0 => v = v(x_oo) = ... = v(x_{n+1}) */
     141          56 :   x = gsub(x, gmul2n(r,-1));
     142          56 :   if (padicprec_relative(x) > v) x = gcvtop(x, p, v);
     143             :   /* x = x_n */
     144         154 :   for (i = n; i > 1; i--)
     145             :   {
     146          98 :     GEN ar = gmul(gel(A,i),gel(R,i)), xp;
     147          98 :     setvalp(ar, valp(ar)+va); /* A_i = A[i] * p^va */
     148             :     /* x_{i-1} = x_i + a_i r_i / x_i - r_{i-1}/2 */
     149          98 :     xp = gsub(gadd(x, gdiv(ar, x)), gmul2n(gel(R,i-1),-1));
     150             :     /* y_{i-1} = y_i (1 - a_i r_i / x^2) */
     151          98 :     if (pty) *pty = gmul(*pty, gsubsg(1, gdiv(ar,gsqr(x))));
     152          98 :     x = xp;
     153             :   }
     154          56 :   *ptx = x;
     155          56 : }
     156             : 
     157             : /* Let T = 4x^3 + b2 x^2 + 2b4 x + b6, where T has a unique p-adic root 'a'.
     158             :  * Return a lift of a to padic accuracy prec. We have
     159             :  * 216 T = 864 X^3 - 18 c4X - c6, where X = x + b2/12 */
     160             : static GEN
     161         266 : doellQp_root(GEN E, long prec)
     162             : {
     163         266 :   GEN c4=ell_get_c4(E), c6=ell_get_c6(E), j=ell_get_j(E), p=ellQp_get_p(E);
     164             :   GEN c6p, T, a;
     165             :   long alpha;
     166         266 :   int pis2 = absequaliu(p, 2);
     167         266 :   if (Q_pval(j, p) >= 0) pari_err_DOMAIN(".root", "v_p(j)", ">=", gen_0, j);
     168             :   /* v(j) < 0 => v(c4^3) = v(c6^2) = 2 alpha */
     169         266 :   alpha = Q_pvalrem(ell_get_c4(E), p, &c4) >> 1;
     170         266 :   if (alpha) (void)Q_pvalrem(ell_get_c6(E), p, &c6);
     171             :   /* Renormalized so that v(c4) = v(c6) = 0; multiply by p^alpha at the end */
     172         266 :   if (prec < 4 && pis2) prec = 4;
     173         266 :   c6p = modii(c6,p);
     174         266 :   if (pis2)
     175             :   { /* Use 432T(X/4) = 27X^3 - 9c4 X - 2c6 to have integral root; a=0 mod 2 */
     176         140 :     T = mkpoln(4, utoipos(27), gen_0, mulis(c4,-9), mulis(c6, -2));
     177             :     /* v_2(root a) = 1, i.e. will lose one bit of accuracy: prec+1 */
     178         140 :     a = ZpX_liftroot(T, gen_0, p, prec+1);
     179         140 :     alpha -= 2;
     180             :   }
     181         126 :   else if (absequaliu(p, 3))
     182             :   { /* Use 216T(X/3) = 32X^3 - 6c4 X - c6 to have integral root; a=-c6 mod 3 */
     183          56 :     a = Fp_neg(c6p, p);
     184          56 :     T = mkpoln(4, utoipos(32), gen_0, mulis(c4, -6), negi(c6));
     185          56 :     a = ZX_Zp_root(T, a, p, prec);
     186          56 :     switch(lg(a)-1)
     187             :     {
     188             :       case 1: /* single root */
     189          28 :         a = gel(a,1); break;
     190             :       case 3: /* three roots, e.g. "15a1", choose the right one */
     191             :       {
     192          28 :         GEN a1 = gel(a,1), a2 = gel(a,2), a3 = gel(a,3);
     193          28 :         long v1 = Z_lval(subii(a2, a3), 3);
     194          28 :         long v2 = Z_lval(subii(a1, a3), 3);
     195          28 :         long v3 = Z_lval(subii(a1, a2), 3);
     196          28 :         if      (v1 == v2) a = a3;
     197           0 :         else if (v1 == v3) a = a2;
     198           0 :         else a = a1;
     199             :       }
     200          28 :       break;
     201             :     }
     202          56 :     alpha--;
     203             :   }
     204             :   else
     205             :   { /* p != 2,3: T = 4(x-a)(x-b)^2 = 4x^3 - 3a^2 x - a^3 when b = -a/2
     206             :      * (so that the trace coefficient vanishes) => a = c6/6c4 (mod p)*/
     207          70 :     GEN c4p = modii(c4,p);
     208          70 :     a = Fp_div(c6p, Fp_mulu(c4p, 6, p), p);
     209          70 :     T = mkpoln(4, utoipos(864), gen_0, mulis(c4, -18), negi(c6));
     210          70 :     a = ZpX_liftroot(T, a, p, prec);
     211             :   }
     212         266 :   a = cvtop(a, p, prec);
     213         266 :   if (alpha) setvalp(a, valp(a)+alpha);
     214         266 :   return gsub(a, gdivgs(ell_get_b2(E), 12));
     215             : }
     216             : GEN
     217         497 : ellQp_root(GEN E, long prec)
     218         497 : { return obj_checkbuild_padicprec(E, Qp_ROOT, &doellQp_root, prec); }
     219             : 
     220             : /* compute a,b such that E1: y^2 = x(x-a)(x+a-b) ~ E */
     221             : static void
     222         315 : doellQp_ab(GEN E, GEN *pta, GEN *ptb, long prec)
     223             : {
     224         315 :   GEN b2 = ell_get_b2(E), b4 = ell_get_b4(E), e1 = ellQp_root(E, prec);
     225         315 :   GEN w, u, t = gadd(gdivgs(b2,4), gmulsg(3,e1)), p = ellQp_get_p(E);
     226         315 :   w = Qp_sqrt(gmul2n(gadd(b4,gmul(e1,gadd(b2,gmulsg(6,e1)))),1));
     227         315 :   u = gadd(t,w);
     228             :   /* Decide between w and -w: we want v(a-b) > v(b) */
     229         315 :   if (absequaliu(p,2))
     230         203 :   { if (valp(u)-1 <= valp(w)) w = gneg_i(w); }
     231             :   else
     232         112 :   { if (valp(u) <= valp(w)) w = gneg_i(w); }
     233             : 
     234             :   /* w^2 = 2b4 + 2b2 e1 + 12 e1^2 = 4(e1-e2)(e1-e3) */
     235         315 :   *pta = gmul2n(gsub(w,t),-2);
     236         315 :   *ptb = gmul2n(w,-1);
     237         315 : }
     238             : 
     239             : static GEN
     240         112 : doellQp_Tate(GEN E, long prec0)
     241             : {
     242         112 :   GEN p = ellQp_get_p(E), j = ell_get_j(E);
     243             :   GEN L, u, u2, q, x1, a, b, d, s, t, AB, A, M2;
     244         112 :   long v, n, pp, prec = prec0+3;
     245         112 :   int split = -1; /* unknown */
     246         112 :   int pis2 = equaliu(p,2);
     247             : 
     248         112 :   if (Q_pval(j, p) >= 0) pari_err_DOMAIN(".tate", "v_p(j)", ">=", gen_0, j);
     249             : START:
     250         315 :   doellQp_ab(E, &a, &b, prec);
     251         315 :   d = gsub(a,b);
     252         315 :   v = prec0 - precp(d);
     253         315 :   if (v > 0) { prec += v; goto START; }
     254         245 :   AB = Qp_agm2_sequence(a,b);
     255         245 :   A = gel(AB,1);
     256         245 :   n = lg(A)-1; /* AGM iterations */
     257         245 :   pp = minss(precp(a),precp(b));
     258         245 :   M2 = cvtop(gel(A,n), p, pis2? pp-2*n: pp);
     259         245 :   setvalp(M2, valp(a));
     260         245 :   u2 = ginv(gmul2n(M2, 2));
     261         245 :   if (split < 0) split = issquare(u2);
     262         245 :   x1 = gen_0;
     263         245 :   Qp_descending_Landen(AB,&x1,NULL);
     264             : 
     265         245 :   t = gaddsg(1, ginv(gmul2n(gmul(u2,x1),1)));
     266         245 :   s = Qp_sqrt(gsubgs(gsqr(t), 1));
     267         245 :   q = gadd(t,s);
     268         245 :   if (gequal0(q)) q = gsub(t,s);
     269         245 :   v = prec0 - precp(q);
     270         245 :   if (split)
     271             :   { /* we want log q at precision prec0 */
     272          98 :     GEN q0 = leafcopy(q); setvalp(q0, 0);
     273          98 :     v +=  valp(gsubgs(q0,1));
     274             :   }
     275         245 :   if (v > 0) { prec += v; goto START; }
     276         112 :   if (valp(q) < 0) q = ginv(q);
     277         112 :   if (split)
     278             :   {
     279          56 :     u = Qp_sqrt(u2);
     280          56 :     L = gdivgs(Qp_log(q), valp(q));
     281             :   }
     282             :   else
     283             :   {
     284          56 :     GEN T = mkpoln(3, gen_1, gen_0, gneg(u2));
     285          56 :     u = mkpolmod(pol_x(0), T);
     286          56 :     L = gen_1;
     287             :   }
     288         112 :   return mkvecn(6, u2, u, q, mkvec2(a, b), L, AB);
     289             : }
     290             : static long
     291         609 : Tate_prec(GEN T) { return padicprec_relative(gel(T,3)); }
     292             : GEN
     293         714 : ellQp_Tate_uniformization(GEN E, long prec)
     294         714 : { return obj_checkbuild_prec(E,Qp_TATE,&doellQp_Tate, &Tate_prec,prec); }
     295             : GEN
     296         168 : ellQp_u(GEN E, long prec)
     297         168 : { GEN T = ellQp_Tate_uniformization(E, prec); return gel(T,2); }
     298             : GEN
     299          77 : ellQp_u2(GEN E, long prec)
     300          77 : { GEN T = ellQp_Tate_uniformization(E, prec); return gel(T,1); }
     301             : GEN
     302         126 : ellQp_q(GEN E, long prec)
     303         126 : { GEN T = ellQp_Tate_uniformization(E, prec); return gel(T,3); }
     304             : GEN
     305         105 : ellQp_ab(GEN E, long prec)
     306         105 : { GEN T = ellQp_Tate_uniformization(E, prec); return gel(T,4); }
     307             : GEN
     308           7 : ellQp_L(GEN E, long prec)
     309           7 : { GEN T = ellQp_Tate_uniformization(E, prec); return gel(T,5); }
     310             : GEN
     311         147 : ellQp_AGM(GEN E, long prec)
     312         147 : { GEN T = ellQp_Tate_uniformization(E, prec); return gel(T,6); }
     313             : 
     314             : /* FORMAL GROUP */
     315             : 
     316             : /* t to w := -1/y */
     317             : GEN
     318         497 : ellformalw(GEN e, long n, long v)
     319             : {
     320         497 :   pari_sp av = avma, av2;
     321             :   GEN a1,a2,a3,a4,a6, a63;
     322         497 :   GEN w = cgetg(3, t_SER), t, U, V, W, U2;
     323         497 :   ulong mask, nold = 1;
     324         497 :   if (v < 0) v = 0;
     325         497 :   if (n <= 0) pari_err_DOMAIN("ellformalw","precision","<=",gen_0,stoi(n));
     326         483 :   mask = quadratic_prec_mask(n);
     327         483 :   t = pol_x(v);
     328         483 :   checkell(e);
     329         483 :   a1 = ell_get_a1(e); a2 = ell_get_a2(e); a3 = ell_get_a3(e);
     330         483 :   a4 = ell_get_a4(e); a6 = ell_get_a6(e); a63 = gmulgs(a6,3);
     331         483 :   w[1] = evalsigne(1)|evalvarn(v)|evalvalp(3);
     332         483 :   gel(w,2) = gen_1; /* t^3 + O(t^4) */
     333             :   /* use Newton iteration, doubling accuracy at each step
     334             :    *
     335             :    *            w^3 a6 + w^2(a4 t + a3) + w (a2 t^2 + a1 t - 1) + t^3
     336             :    * w  <-  w - -----------------------------------------------------
     337             :    *              w^2 (3a6) + w (2a4 t + 2a3) + (a2 t^2 + a1 t - 1)
     338             :    *
     339             :    *              w^3 a6 + w^2 U + w V + W
     340             :    *      =: w -  -----------------------
     341             :    *                w^2 (3a6) + 2w U + V
     342             :    */
     343         483 :   U = gadd(gmul(a4,t), a3);
     344         483 :   U2 = gmul2n(U,1);
     345         483 :   V = gsubgs(gadd(gmul(a2,gsqr(t)), gmul(a1,t)), 1);
     346         483 :   W = gpowgs(t,3);
     347         483 :   av2 = avma;
     348        2744 :   while (mask > 1)
     349             :   { /* nold correct terms in w */
     350        1778 :     ulong i, nnew = nold << 1;
     351             :     GEN num, den, wnew, w2, w3;
     352        1778 :     if (mask & 1) nnew--;
     353        1778 :     mask >>= 1;
     354        1778 :     wnew = cgetg(nnew+2, t_SER);
     355        1778 :     wnew[1] = w[1];
     356        1778 :     for (i = 2; i < nold+2; i++) gel(wnew,i) = gel(w,i);
     357        1778 :     for (     ; i < nnew+2; i++) gel(wnew,i) = gen_0;
     358        1778 :     w = wnew;
     359        1778 :     w2 = gsqr(w); w3 = gmul(w2,w);
     360        1778 :     num = gadd(gmul(a6,w3), gadd(gmul(U,w2), gadd(gmul(V,w), W)));
     361        1778 :     den = gadd(gmul(a63,w2), gadd(gmul(w,U2), V));
     362             : 
     363        1778 :     w = gerepileupto(av2, gsub(w, gdiv(num, den)));
     364        1778 :     nold = nnew;
     365             :   }
     366         483 :   return gerepilecopy(av, w);
     367             : }
     368             : 
     369             : static GEN
     370         266 : ellformalpoint_i(GEN w, GEN wi)
     371         266 : { return mkvec2(gmul(pol_x(varn(w)),wi), gneg(wi)); }
     372             : 
     373             : /* t to [x,y] */
     374             : GEN
     375          70 : ellformalpoint(GEN e, long n, long v)
     376             : {
     377          70 :   pari_sp av = avma;
     378          70 :   GEN w = ellformalw(e, n, v), wi = ser_inv(w);
     379          70 :   return gerepilecopy(av, ellformalpoint_i(w, wi));
     380             : }
     381             : 
     382             : static GEN
     383         343 : ellformaldifferential_i(GEN e, GEN w, GEN wi, GEN *px)
     384             : {
     385             :   GEN x, w1;
     386         343 :   if (gequal0(ell_get_a1(e)) && gequal0(ell_get_a3(e)))
     387             :   { /* dx/2y = dx * -w/2, avoid division */
     388         147 :     x = gmul(pol_x(varn(w)), wi);
     389         147 :     w1 = gmul(derivser(x), gneg(gmul2n(w,-1)));
     390             :   }
     391             :   else
     392             :   {
     393         196 :     GEN P = ellformalpoint_i(w, wi);
     394         196 :     x = gel(P,1);
     395         196 :     w1 = gdiv(derivser(x), ec_dmFdy_evalQ(e, P));
     396             :   }
     397         343 :   *px = x; return w1;
     398             : }
     399             : /* t to [ dx / (2y + a1 x + a3), x * ... ]*/
     400             : GEN
     401          70 : ellformaldifferential(GEN e, long n, long v)
     402             : {
     403          70 :   pari_sp av = avma;
     404          70 :   GEN w = ellformalw(e, n, v), wi = ser_inv(w), x;
     405          70 :   GEN w1 = ellformaldifferential_i(e, w, wi, &x);
     406          70 :   return gerepilecopy(av, mkvec2(w1,gmul(x,w1)));
     407             : }
     408             : 
     409             : /* t to z, dz = w1 dt */
     410             : GEN
     411         147 : ellformallog(GEN e, long n, long v)
     412             : {
     413         147 :   pari_sp av = avma;
     414         147 :   GEN w = ellformalw(e, n, v), wi = ser_inv(w), x;
     415         147 :   GEN w1 = ellformaldifferential_i(e, w, wi, &x);
     416         147 :   return gerepileupto(av, integser(w1));
     417             : }
     418             : /* z to t */
     419             : GEN
     420          70 : ellformalexp(GEN e, long n, long v)
     421             : {
     422          70 :   pari_sp av = avma;
     423          70 :   return gerepileupto(av, serreverse(ellformallog(e,n,v)));
     424             : }
     425             : /* [log_p (sigma(t) / t), log_E t], as power series, d (log_E t) := w1 dt;
     426             :  * As a fonction of z: odd, = e.b2/12 * z + O(z^3).
     427             :  *   sigma(z) = ellsigma(e) exp(e.b2/24*z^2)
     428             :  * log_p(sigma(t)/t)=log(subst(sigma(z), x, ellformallog(e))/x) */
     429             : static GEN
     430         126 : ellformallogsigma_t(GEN e, long n)
     431             : {
     432         126 :   pari_sp av = avma;
     433         126 :   GEN w = ellformalw(e, n, 0), wi = ser_inv(w), t = pol_x(0);
     434         126 :   GEN x, s = ellformaldifferential_i(e, w, wi, &x);
     435         126 :   GEN f = gmul(s, gadd(integser(gmul(x,s)), gmul2n(ell_get_a1(e),-1)));
     436         126 :   return gerepilecopy(av, mkvec2(integser( gsub(ginv(gneg(t)), f) ),
     437             :                                  integser(s)));
     438             : }
     439             : 
     440             : /* P-ADIC HEIGHTS */
     441             : 
     442             : /* m >= 0, T = b6^2, g4 = b6^2 - b4 b8, return g_m(xP) mod N, in Mazur-Tate's
     443             :  * notation (Duke 1991)*/
     444             : static GEN
     445        4102 : rellg(hashtable *H, GEN m, GEN T, GEN g4, GEN b8, GEN N)
     446             : {
     447             :   hashentry *h;
     448             :   GEN n, z, np2, np1, nm2, nm1, fp2, fp1, fm2, fm1, f;
     449             :   ulong m4;
     450        4102 :   if (abscmpiu(m, 4) <= 0) switch(itou(m))
     451             :   {
     452         196 :     case 0: return gen_0;
     453         532 :     case 1: return gen_1;
     454         630 :     case 2: return subiu(N,1);
     455         728 :     case 3: return b8;
     456         756 :     case 4: return g4;
     457             :   }
     458        1260 :   if ((h = hash_search(H, (void*)m))) return (GEN)h->val;
     459         686 :   m4 = mod4(m);
     460         686 :   n = shifti(m, -1); f   = rellg(H,n,T,g4,b8,N);
     461         686 :   np2 = addiu(n, 2); fp2 = rellg(H,np2,T,g4,b8,N);
     462         686 :   np1 = addiu(n, 1); fp1 = rellg(H,np1,T,g4,b8,N);
     463         686 :   nm2 = subiu(n, 2); fm2 = rellg(H,nm2,T,g4,b8,N);
     464         686 :   nm1 = subiu(n, 1); fm1 = rellg(H,nm1,T,g4,b8,N);
     465         686 :   if (odd(m4))
     466             :   {
     467         406 :     GEN t1 = Fp_mul(fp2, Fp_powu(f,3,N), N);
     468         406 :     GEN t2 = Fp_mul(fm1, Fp_powu(fp1,3,N), N);
     469         406 :     if (mpodd(n)) t2 = Fp_mul(T,t2,N); else t1 = Fp_mul(T,t1,N);
     470         406 :     z = Fp_sub(t1, t2, N);
     471             :   }
     472             :   else
     473             :   {
     474         280 :     GEN t1 = Fp_mul(fm2, Fp_sqr(fp1,N), N);
     475         280 :     GEN t2 = Fp_mul(fp2, Fp_sqr(fm1,N), N);
     476         280 :     z = Fp_mul(f, Fp_sub(t1, t2, N), N);
     477             :   }
     478         686 :   hash_insert(H, (void*)m, (void*)z);
     479         686 :   return z;
     480             : }
     481             : 
     482             : static GEN
     483         672 : addii3(GEN x, GEN y, GEN z) { return addii(x,addii(y,z)); }
     484             : static GEN
     485         448 : addii4(GEN x, GEN y, GEN z, GEN t) { return addii(x,addii3(y,z,t)); }
     486             : static GEN
     487         224 : addii5(GEN x, GEN y, GEN z, GEN t, GEN u) { return addii(x,addii4(y,z,t,u)); }
     488             : 
     489             : /* xP = [n,d] (corr. to n/d, coprime), such that the reduction of the point
     490             :  * P = [xP,yP] is non singular at all places. Return x([m] P) mod N as
     491             :  * [num,den] (coprime) */
     492             : static GEN
     493         224 : xmP(GEN e, GEN xP, GEN m, GEN N)
     494             : {
     495         224 :   pari_sp av = avma;
     496         224 :   ulong k = expi(m);
     497         224 :   hashtable *H = hash_create((5+k)*k, (ulong(*)(void*))&hash_GEN,
     498             :                                       (int(*)(void*,void*))&gidentical, 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 gerepilecopy(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(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(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*LOG2/(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 - LOG2/(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) || !signe(remii(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 (signe(remii(t,p))) { avma = av; return 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 :   avma = av; return signe(remii(t,p))? 0: 1;
     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 non-singular; 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 :   checkellpt(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 gerepilecopy(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 :   checkellpt(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 = gel(P,1);
     695         126 :   x = mkvec2(numer(x), denom(x));
     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         112 :     }
     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 gerepilecopy(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 :   checkellpt(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 gerepileupto(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 = gel(q,3);
     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 = gerepileuptoint(av, s);
     782             :   }
     783          14 :   s = subui(1, muliu(s,24));
     784          14 :   s = gdivgs(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, utoi(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",utoi(p));
     798             :   }
     799         182 :   return hyperellpadicfrobenius(ec_bmodel(E), 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(p,n), 0);
     817           7 :     l = ellpadics2_tate(Ep, n);
     818           7 :     if (Ep != E) obj_free(Ep);
     819           7 :     return gerepilecopy(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(ap, q, shifti(q,-1));
     834             :   }
     835          21 :   l = mspadic_unit_eigenvalue(ap, 2, p, n);
     836          21 :   return gerepileupto(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 gerepileupto(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 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 :   W = ellpadicL_symbol(E, p, gen_0, D);
     920          28 :   ED = D? ellinit(elltwist(E,D), gen_1, 0): E;
     921          28 :   ED = ellanal_globalred_all(ED, NULL, &ND, &tam);
     922          28 :   vN = Z_pval(ND, p); /* additive reduction ? */
     923          28 :   if (vN >= 2) pari_err_DOMAIN("ellpadicbsd","v_p(N(E_D))",">",gen_1,stoi(vN));
     924          28 :   if (n < 5) n = 5;
     925           0 :   for(;; n <<= 1)
     926             :   {
     927          28 :     GEN v = ellpadicL_init(W, n);
     928          63 :     for (r = 0; r < MAXR; r++)
     929             :     {
     930          63 :       U = ellpadic_i(v, r);
     931          63 :       if (!gequal0(U)) break;
     932             :     }
     933          28 :     if (r < MAXR) break;
     934           0 :   }
     935          28 :   apD = ellap(ED, p);
     936          28 :   if (typ(U) == t_COL)
     937             :   { /* p | a_p(E_D), frobenius on E_D */
     938           0 :     GEN F = mkmat22(gen_0, negi(p), gen_1, apD);
     939           0 :     U = RgM_RgC_mul(gpowgs(gsubsg(1, gdiv(F,p)), -2), U);
     940           0 :     settyp(U, t_VEC);
     941             :   }
     942          28 :   else if (dvdii(ND,p))
     943             :   {
     944          21 :     if (equalim1(apD)) /* divide by 1-1/a_p */
     945          14 :       U = gdivgs(U, 2);
     946             :     else
     947             :     { /* ap = 1 */
     948           7 :       GEN EDp = ellinit(ED, zeropadic(p,n), 0);
     949           7 :       U = gdiv(U, ellQp_L(EDp,n));
     950           7 :       obj_free(EDp);
     951             :     }
     952             :   }
     953             :   else
     954             :   {
     955           7 :     GEN a = mspadic_unit_eigenvalue(apD, 2, p, n);
     956           7 :     U = gmul(U, gpowgs(gsubsg(1, ginv(a)), -2));
     957             :   }
     958          28 :   C = mulii(tam, mpfact(r));
     959          28 :   if (D) C = gmul(C, get_Euler(ED, D));
     960          28 :   C = gdiv(sqru(torsion_order(ED)), C);
     961          28 :   if (D) obj_free(ED);
     962          28 :   return gerepilecopy(av, mkvec2(utoi(r), gmul(U, C)));
     963             : }
     964             : 
     965             : GEN
     966           7 : ellpadicregulator(GEN E, GEN p, long n, GEN S)
     967             : {
     968           7 :   pari_sp av = avma;
     969           7 :   GEN FG = ellpadicheightmatrix(E,p,n,S); /* forbids additive reduction */
     970             :   /* [F,G]: height in basis [omega, eta] */
     971           7 :   GEN R, F = gel(FG,1), G = gel(FG,2), ap = ellap(E,p);
     972           7 :   if (dvdii(ap, p))
     973             :   { /* supersingular */
     974           0 :     GEN f = ellpadicfrobenius(E, itou(p), n);
     975           0 :     GEN x = gcoeff(f,1,1), y = gcoeff(f,2,1);
     976             :     /* [A,B]: regulator height in basis [omega, eta] */
     977           0 :     GEN A = det(F), B = det(gadd(F,G)), C;
     978           0 :     C = gdiv(gsub(B,A), y);
     979             :     /* R: regulator height in basis [omega, F.omega] */
     980           0 :     R = mkvec2(gsub(A, gmul(x,C)), C);
     981             :   }
     982             :   else
     983             :   {
     984             :     GEN s2;
     985           7 :     if (equali1(ap) && dvdii(ell_get_disc(E),p))
     986           7 :     { /* split multiplicative reduction */
     987           7 :       GEN Ep = ellinit(E, zeropadic(p,n), 0);
     988           7 :       GEN q = ellQp_q(Ep,n), u2 = ellQp_u2(Ep,n);
     989           7 :       s2 = ellpadics2_tate(Ep, n);
     990           7 :       s2 = gsub(s2, ginv(gmul(Qp_log(q), u2))); /*extended MW group contrib*/
     991           7 :       obj_free(Ep);
     992             :     }
     993             :     else
     994           0 :       s2 = ellpadics2(E,p,n);
     995           7 :     R = det( RgM_sub(F, RgM_Rg_mul(G,s2)) );
     996             :   }
     997           7 :   return gerepilecopy(av, R);
     998             : }

Generated by: LCOV version 1.11