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 - elltrans.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.18.1 lcov report (development 30074-db9b3c0519) Lines: 1274 1374 92.7 %
Date: 2025-03-13 08:38:27 Functions: 107 111 96.4 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : /********************************************************************/
      16             : /**                                                                **/
      17             : /**          ELLIPTIC and MODULAR FUNCTIONS                        **/
      18             : /**         (as complex or p-adic functions)                       **/
      19             : /**                                                                **/
      20             : /********************************************************************/
      21             : #include "pari.h"
      22             : #include "paripriv.h"
      23             : 
      24             : #define DEBUGLEVEL DEBUGLEVEL_ell
      25             : 
      26             : /********************************************************************/
      27             : /**                       Periods                                  **/
      28             : /********************************************************************/
      29             : /* The complex AGM, periods of elliptic curves over C and complex elliptic
      30             :  * logarithms; John E. Cremona, Thotsaphon Thongjunthug, arXiv:1011.0914 */
      31             : 
      32             : static GEN
      33       44881 : ellomega_agm(GEN a, GEN b, GEN c, long prec)
      34             : {
      35       44881 :   GEN pi = mppi(prec), mIpi = mkcomplex(gen_0, negr(pi));
      36       44881 :   GEN Mac = agm(a,c,prec), Mbc = agm(b,c,prec);
      37       44881 :   retmkvec2(gdiv(pi, Mac), gdiv(mIpi, Mbc));
      38             : }
      39             : 
      40             : static GEN
      41       36651 : ellomega_cx(GEN E, long prec)
      42             : {
      43       36651 :   pari_sp av = avma;
      44       36651 :   GEN roots = ellR_roots(E, prec + EXTRAPREC64);
      45       36651 :   GEN d1=gel(roots,4), d2=gel(roots,5), d3=gel(roots,6);
      46       36651 :   GEN a = gsqrt(d3,prec), b = gsqrt(d1,prec), c = gsqrt(d2,prec);
      47       36651 :   return gerepileupto(av, ellomega_agm(a,b,c,prec));
      48             : }
      49             : 
      50             : /* return [w1,w2] for E / R; w1 > 0 is real.
      51             :  * If e.disc > 0, w2 = -I r; else w2 = w1/2 - I r, for some real r > 0.
      52             :  * => tau = w1/w2 is in upper half plane */
      53             : static GEN
      54       44881 : doellR_omega(GEN E, long prec)
      55             : {
      56       44881 :   pari_sp av = avma;
      57             :   GEN roots, d2, z, a, b, c;
      58       44881 :   if (ellR_get_sign(E) >= 0) return ellomega_cx(E,prec);
      59        8230 :   roots = ellR_roots(E,prec + EXTRAPREC64);
      60        8230 :   d2 = gel(roots,5);
      61        8230 :   z = gsqrt(d2,prec); /* imag(e1-e3) > 0, so that b > 0*/
      62        8230 :   a = gel(z,1); /* >= 0 */
      63        8230 :   b = gel(z,2);
      64        8230 :   c = gabs(z, prec);
      65        8230 :   z = ellomega_agm(a,b,c,prec);
      66        8230 :   return gerepilecopy(av, mkvec2(gel(z,1),gmul2n(gadd(gel(z,1),gel(z,2)),-1)));
      67             : }
      68             : static GEN
      69          60 : doellR_eta(GEN E, long prec)
      70          60 : { GEN w = ellR_omega(E, prec + EXTRAPREC64); return elleta(w, prec); }
      71             : 
      72             : GEN
      73       79416 : ellR_omega(GEN E, long prec)
      74       79416 : { return obj_checkbuild_realprec(E, R_PERIODS, &doellR_omega, prec); }
      75             : GEN
      76          72 : ellR_eta(GEN E, long prec)
      77          72 : { return obj_checkbuild_realprec(E, R_ETA, &doellR_eta, prec); }
      78             : 
      79             : /* P = [x,0] is 2-torsion on y^2 = g(x). Return w1/2, (w1+w2)/2, or w2/2
      80             :  * depending on whether x is closest to e1,e2, or e3, the 3 complex root of g */
      81             : static GEN
      82          12 : zell_closest_0(GEN om, GEN x, GEN ro)
      83             : {
      84          12 :   GEN e1 = gel(ro,1), e2 = gel(ro,2), e3 = gel(ro,3);
      85          12 :   GEN d1 = gnorm(gsub(x,e1));
      86          12 :   GEN d2 = gnorm(gsub(x,e2));
      87          12 :   GEN d3 = gnorm(gsub(x,e3));
      88          12 :   GEN z = gel(om,2);
      89          12 :   if (gcmp(d1, d2) <= 0)
      90           0 :   { if (gcmp(d1, d3) <= 0) z = gel(om,1); }
      91             :   else
      92          12 :   { if (gcmp(d2, d3)<=0) z = gadd(gel(om,1),gel(om,2)); }
      93          12 :   return gmul2n(z, -1);
      94             : }
      95             : 
      96             : static GEN
      97       24630 : zellcx(GEN E, GEN P, long prec)
      98             : {
      99       24630 :   GEN R = ellR_roots(E, prec+EXTRAPREC64);
     100       24630 :   GEN x0 = gel(P,1), y0 = ec_dmFdy_evalQ(E,P);
     101       24630 :   if (gequal0(y0))
     102           0 :     return zell_closest_0(ellomega_cx(E,prec),x0,R);
     103             :   else
     104             :   {
     105       24630 :     GEN e2 = gel(R,2), e3 = gel(R,3), d2 = gel(R,5), d3 = gel(R,6);
     106       24630 :     GEN a = gsqrt(d2,prec), b = gsqrt(d3,prec);
     107       24630 :     GEN r = gsqrt(gdiv(gsub(x0,e3), gsub(x0,e2)),prec);
     108       24630 :     GEN t = gdiv(gneg(y0), gmul2n(gmul(r,gsub(x0,e2)),1));
     109       24630 :     GEN ar = real_i(a), br = real_i(b), ai = imag_i(a), bi = imag_i(b);
     110             :     /* |a+b| < |a-b| */
     111       24630 :     if (gcmp(gmul(ar,br), gneg(gmul(ai,bi))) < 0) b = gneg(b);
     112       24630 :     return zellagmcx(a,b,r,t,prec);
     113             :   }
     114             : }
     115             : 
     116             : /* Assume E/R, disc E < 0, and P \in E(R) ==> z \in R */
     117             : static GEN
     118           0 : zellrealneg(GEN E, GEN P, long prec)
     119             : {
     120           0 :   GEN x0 = gel(P,1), y0 = ec_dmFdy_evalQ(E,P);
     121           0 :   if (gequal0(y0)) return gmul2n(gel(ellR_omega(E,prec),1),-1);
     122             :   else
     123             :   {
     124           0 :     GEN R = ellR_roots(E, prec+EXTRAPREC64);
     125           0 :     GEN d2 = gel(R,5), e3 = gel(R,3);
     126           0 :     GEN a = gsqrt(d2,prec);
     127           0 :     GEN z = gsqrt(gsub(x0,e3), prec);
     128           0 :     GEN ar = real_i(a), zr = real_i(z), ai = imag_i(a), zi = imag_i(z);
     129           0 :     GEN t = gdiv(gneg(y0), gmul2n(gnorm(z),1));
     130           0 :     GEN r2 = ginv(gsqrt(gaddsg(1,gdiv(gmul(ai,zi),gmul(ar,zr))),prec));
     131           0 :     return zellagmcx(ar,gabs(a,prec),r2,gmul(t,r2),prec);
     132             :   }
     133             : }
     134             : 
     135             : /* Assume E/R, disc E > 0, and P \in E(R) */
     136             : static GEN
     137          24 : zellrealpos(GEN E, GEN P, long prec)
     138             : {
     139          24 :   GEN R = ellR_roots(E, prec+EXTRAPREC64);
     140          24 :   GEN d2,d3,e1,e2,e3, a,b, x0 = gel(P,1), y0 = ec_dmFdy_evalQ(E,P);
     141          24 :   if (gequal0(y0)) return zell_closest_0(ellR_omega(E,prec), x0,R);
     142          12 :   e1 = gel(R,1);
     143          12 :   e2 = gel(R,2);
     144          12 :   e3 = gel(R,3);
     145          12 :   d2 = gel(R,5);
     146          12 :   d3 = gel(R,6);
     147          12 :   a = gsqrt(d2,prec);
     148          12 :   b = gsqrt(d3,prec);
     149          12 :   if (gcmp(x0,e1)>0) {
     150           6 :     GEN r = gsqrt(gdiv(gsub(x0,e3), gsub(x0,e2)),prec);
     151           6 :     GEN t = gdiv(gneg(y0), gmul2n(gmul(r,gsub(x0,e2)),1));
     152           6 :     return zellagmcx(a,b,r,t,prec);
     153             :   } else {
     154           6 :     GEN om = ellR_omega(E,prec);
     155           6 :     GEN r = gdiv(a,gsqrt(gsub(e1,x0),prec));
     156           6 :     GEN t = gdiv(gmul(r,y0),gmul2n(gsub(x0,e3),1));
     157           6 :     return gsub(zellagmcx(a,b,r,t,prec),gmul2n(gel(om,2),-1));
     158             :   }
     159             : }
     160             : 
     161             : static void
     162          18 : ellQp_P2t_err(GEN E, GEN z)
     163             : {
     164          18 :   if (typ(ellQp_u(E,1)) == t_POLMOD)
     165          18 :     pari_err_IMPL("ellpointtoz when u not in Qp");
     166           0 :   pari_err_DOMAIN("ellpointtoz", "point", "not on", strtoGENstr("E"),z);
     167           0 : }
     168             : static GEN
     169         156 : get_r0(GEN E, long prec)
     170             : {
     171         156 :   GEN b2 = ell_get_b2(E), e1 = ellQp_root(E, prec);
     172         156 :   return gadd(e1,gmul2n(b2,-2));
     173             : }
     174             : static GEN
     175         114 : ellQp_P2t(GEN E, GEN P, long prec)
     176             : {
     177         114 :   pari_sp av = avma;
     178             :   GEN a, b, ab, c0, r0, ar, r, x, delta, x1, y1, t, u, q;
     179             :   long vq, vt, Q, R;
     180         114 :   if (ell_is_inf(P)) return gen_1;
     181         108 :   ab = ellQp_ab(E, prec); a = gel(ab,1); b = gel(ab,2);
     182         108 :   u = ellQp_u(E, prec);
     183         108 :   q = ellQp_q(E, prec);
     184         108 :   x = gel(P,1);
     185         108 :   r0 = get_r0(E, prec);
     186         108 :   c0 = gadd(x, gmul2n(r0,-1));
     187         108 :   if (typ(c0) != t_PADIC || !is_scalar_t(typ(gel(P,2))))
     188           6 :     pari_err_TYPE("ellpointtoz",P);
     189         102 :   r = gsub(a,b);
     190         102 :   ar = gmul(a, r);
     191         102 :   if (gequal0(c0))
     192             :   {
     193           6 :     x1 = Qp_sqrt(gneg(ar));
     194           6 :     if (!x1) ellQp_P2t_err(E,P);
     195             :   }
     196             :   else
     197             :   {
     198          96 :     delta = gdiv(ar, gsqr(c0));
     199          96 :     t = Qp_sqrt(gsubsg(1,gmul2n(delta,2)));
     200          96 :     if (!t) ellQp_P2t_err(E,P);
     201          90 :     x1 = gmul(gmul2n(c0,-1), gaddsg(1,t));
     202             :   }
     203          96 :   y1 = gsubsg(1, gdiv(ar, gsqr(x1)));
     204          96 :   if (gequal0(y1))
     205             :   {
     206          12 :     y1 = Qp_sqrt(gmul(x1, gmul(gadd(x1, a), gadd(x1, r))));
     207          12 :     if (!y1) ellQp_P2t_err(E,P);
     208             :   }
     209             :   else
     210          84 :     y1 = gdiv(gmul2n(ec_dmFdy_evalQ(E,P), -1), y1);
     211          84 :   Qp_descending_Landen(ellQp_AGM(E,prec), &x1,&y1);
     212             : 
     213          84 :   t = gmul(u, gmul2n(y1,1)); /* 2u y_oo */
     214          84 :   t = gdiv(gsub(t, x1), gadd(t, x1));
     215             :   /* Reduce mod q^Z: we want 0 <= v(t) < v(q) */
     216          84 :   if (typ(t) == t_PADIC)
     217          48 :     vt = valp(t);
     218             :   else
     219          36 :     vt = valp(gnorm(t)) / 2; /* v(t) = v(Nt) / (e*f) */
     220          84 :   vq = valp(q); /* > 0 */
     221          84 :   Q = vt / vq; R = vt % vq; if (R < 0) Q--;
     222          84 :   if (Q) t = gdiv(t, gpowgs(q,Q));
     223          84 :   if (padicprec_relative(t) > prec) t = gprec(t, prec);
     224          84 :   return gerepileupto(av, t);
     225             : }
     226             : 
     227             : static GEN
     228          48 : ellQp_t2P(GEN E, GEN t, long prec)
     229             : {
     230          48 :   pari_sp av = avma;
     231             :   GEN AB, A, R, x0,x1, y0,y1, u, u2, r0, s0, ar;
     232             :   long v;
     233          48 :   if (gequal1(t)) return ellinf();
     234             : 
     235          48 :   AB = ellQp_AGM(E,prec); A = gel(AB,1); R = gel(AB,3); v = itos(gel(AB,4));
     236          48 :   u = ellQp_u(E,prec);
     237          48 :   u2= ellQp_u2(E,prec);
     238          48 :   x1 = gdiv(t, gmul(u2, gsqr(gsubsg(1,t))));
     239          48 :   y1 = gdiv(gmul(x1,gaddsg(1,t)), gmul(gmul2n(u,1),gsubsg(1,t)));
     240          48 :   Qp_ascending_Landen(AB, &x1,&y1);
     241          48 :   r0 = get_r0(E, prec);
     242             : 
     243          48 :   ar = gmul(gel(A,1), gel(R,1)); setvalp(ar, valp(ar)+v);
     244          48 :   x0 = gsub(gadd(x1, gdiv(ar, x1)), gmul2n(r0,-1));
     245          48 :   s0 = gmul2n(ec_h_evalx(E, x0), -1);
     246          48 :   y0 = gsub(gmul(y1, gsubsg(1, gdiv(ar,gsqr(x1)))), s0);
     247          48 :   return gerepilecopy(av, mkvec2(x0,y0));
     248             : }
     249             : 
     250             : static GEN
     251       24654 : zell_i(GEN e, GEN z, long prec)
     252             : {
     253             :   GEN t;
     254             :   long s;
     255       24654 :   (void)ellR_omega(e, prec); /* type checking */
     256       24654 :   if (ell_is_inf(z)) return gen_0;
     257       24654 :   s = ellR_get_sign(e);
     258       24654 :   if (s && typ(gel(z,1))!=t_COMPLEX && typ(gel(z,2))!=t_COMPLEX)
     259          24 :     t = (s < 0)? zellrealneg(e,z,prec): zellrealpos(e,z,prec);
     260             :   else
     261       24630 :     t = zellcx(e,z,prec);
     262       24654 :   return t;
     263             : }
     264             : 
     265             : GEN
     266       24774 : zell(GEN E, GEN P, long prec)
     267             : {
     268       24774 :   pari_sp av = avma;
     269       24774 :   checkell(E);
     270       24774 :   if (!checkellpt_i(P)) pari_err_TYPE("ellpointtoz", P);
     271       24762 :   switch(ell_get_type(E))
     272             :   {
     273         114 :     case t_ELL_Qp:
     274         114 :       prec = minss(ellQp_get_prec(E), padicprec_relative(P));
     275         114 :       return ellQp_P2t(E, P, prec);
     276           6 :     case t_ELL_NF:
     277             :     {
     278           6 :       GEN Ee = ellnfembed(E, prec), Pe = ellpointnfembed(E, P, prec);
     279           6 :       long i, l = lg(Pe);
     280          18 :       for (i = 1; i < l; i++) gel(Pe,i) = zell_i(gel(Ee,i), gel(Pe,i), prec);
     281           6 :       ellnfembed_free(Ee); return gerepilecopy(av, Pe);
     282             :     }
     283          12 :     case t_ELL_Q: break;
     284       24630 :     case t_ELL_Rg: break;
     285           0 :     default: pari_err_TYPE("ellpointtoz", E);
     286             :   }
     287       24642 :   return gerepileupto(av, zell_i(E, P, prec));
     288             : }
     289             : 
     290             : /********************************************************************/
     291             : /**                COMPLEX ELLIPTIC FUNCTIONS                      **/
     292             : /********************************************************************/
     293             : 
     294             : enum period_type { t_PER_W, t_PER_WETA, t_PER_ELL };
     295             : /* normalization / argument reduction for elliptic functions */
     296             : typedef struct {
     297             :   enum period_type type;
     298             :   GEN in; /* original input */
     299             :   GEN w1,w2,tau; /* original basis for L = <w1,w2> = w2 <1,tau> */
     300             :   GEN W1,W2,Tau; /* new basis for L = <W1,W2> = W2 <1,tau> */
     301             :   GEN a,b,c,d; /* t_INT; tau in F = h/Sl2, tau = g.t, g=[a,b;c,d] in SL(2,Z) */
     302             :   GEN z,Z; /* z/w2 defined mod <1,tau>, Z = z/w2 + x*tau+y reduced mod <1,tau>*/
     303             :   GEN x,y; /* t_INT */
     304             :   int swap; /* 1 if we swapped w1 and w2 */
     305             :   int some_q_is_real; /* exp(2iPi g.tau) for some g \in SL(2,Z) */
     306             :   int some_z_is_real; /* z + xw1 + yw2 is real for some x,y \in Z */
     307             :   int some_z_is_pure_imag; /* z + xw1 + yw2 in i*R */
     308             :   int q_is_real; /* exp(2iPi tau) \in R */
     309             :   int abs_u_is_1; /* |exp(2iPi Z)| = 1 */
     310             :   long prec; /* precision(Z) */
     311             :   long prec0; /* required precision for result */
     312             : } ellred_t;
     313             : 
     314             : /* compute g in SL_2(Z), g.t is in the usual
     315             :    fundamental domain. Internal function no check, no garbage. */
     316             : static void
     317       95022 : set_gamma(GEN *pt, GEN *pa, GEN *pb, GEN *pc, GEN *pd)
     318             : {
     319       95022 :   GEN a, b, c, d, t, t0 = *pt, run = dbltor(1. - 1e-8);
     320       95022 :   long e = gexpo(gel(t0,2));
     321       95022 :   if (e < 0) t0 = gprec_wensure(t0, precision(t0)+nbits2extraprec(-e));
     322       95022 :   t = t0;
     323       95022 :   a = d = gen_1;
     324       95022 :   b = c = gen_0;
     325             :   for(;;)
     326       32058 :   {
     327      127080 :     GEN m, n = ground(gel(t,1));
     328      127080 :     if (signe(n))
     329             :     { /* apply T^n */
     330       40834 :       t = gsub(t,n);
     331       40834 :       a = subii(a, mulii(n,c));
     332       40834 :       b = subii(b, mulii(n,d));
     333             :     }
     334      127080 :     m = cxnorm(t); if (gcmp(m,run) > 0) break;
     335       32058 :     t = gneg_i(gdiv(conj_i(t), m)); /* apply S */
     336       32058 :     togglesign_safe(&c); swap(a,c);
     337       32058 :     togglesign_safe(&d); swap(b,d);
     338             :   }
     339       95022 :   if (e < 0 && (signe(b) || signe(c))) *pt = t0;
     340       95022 :   *pa = a; *pb = b; *pc = c; *pd = d;
     341       95022 : }
     342             : /* Im z > 0. Return U.z in PSl2(Z)'s standard fundamental domain.
     343             :  * Set *pU to U. */
     344             : GEN
     345       32508 : cxredsl2_i(GEN z, GEN *pU, GEN *czd)
     346             : {
     347             :   GEN a,b,c,d;
     348       32508 :   set_gamma(&z, &a, &b, &c, &d);
     349       32508 :   *pU = mkmat2(mkcol2(a,c), mkcol2(b,d));
     350       32508 :   *czd = gadd(gmul(c,z), d);
     351       32508 :   return gdiv(gadd(gmul(a,z), b), *czd);
     352             : }
     353             : GEN
     354       32472 : cxredsl2(GEN t, GEN *pU)
     355             : {
     356       32472 :   pari_sp av = avma;
     357             :   GEN czd;
     358       32472 :   t = cxredsl2_i(t, pU, &czd);
     359       32472 :   return gc_all(av, 2, &t, pU);
     360             : }
     361             : 
     362             : /* swap w1, w2 so that Im(t := w1/w2) > 0. Set tau = representative of t in
     363             :  * the standard fundamental domain, and g in Sl_2, such that tau = g.t */
     364             : static void
     365       62514 : red_modSL2(ellred_t *T, long prec)
     366             : {
     367             :   long s, p;
     368       62514 :   T->tau = gdiv(T->w1,T->w2);
     369       62514 :   if (isintzero(real_i(T->tau))) T->some_q_is_real = 1;
     370       62514 :   s = gsigne(imag_i(T->tau));
     371       62514 :   if (!s) pari_err_DOMAIN("elliptic function", "det(w1,w2)", "=", gen_0,
     372             :                           mkvec2(T->w1,T->w2));
     373       62514 :   T->swap = (s < 0);
     374       62514 :   if (T->swap) { swap(T->w1, T->w2); T->tau = ginv(T->tau); }
     375       62514 :   p = precision(T->tau); T->prec0 = p? p: prec;
     376       62514 :   set_gamma(&T->tau, &T->a, &T->b, &T->c, &T->d);
     377             :   /* update lattice */
     378       62514 :   p = precision(T->tau);
     379       62514 :   if (p)
     380             :   {
     381       62190 :     T->w1 = gprec_wensure(T->w1, p);
     382       62190 :     T->w2 = gprec_wensure(T->w2, p);
     383             :   }
     384       62514 :   T->W1 = gadd(gmul(T->a,T->w1), gmul(T->b,T->w2));
     385       62514 :   T->W2 = gadd(gmul(T->c,T->w1), gmul(T->d,T->w2));
     386       62514 :   T->Tau = gdiv(T->W1, T->W2);
     387       62514 :   if (isintzero(real_i(T->Tau))) T->some_q_is_real = T->q_is_real = 1;
     388       62514 :   p = precision(T->Tau); T->prec = p? p: prec;
     389       62514 : }
     390             : /* is z real or pure imaginary ? */
     391             : static void
     392       67728 : check_complex(GEN z, int *real, int *imag)
     393             : {
     394       67728 :   if (typ(z) != t_COMPLEX)      { *real = 1; *imag = 0; }
     395       55878 :   else if (isintzero(gel(z,1))) { *real = 0; *imag = 1; }
     396       50820 :   else *real = *imag = 0;
     397       67728 : }
     398             : static void
     399       33906 : reduce_z(GEN z, ellred_t *T)
     400             : {
     401             :   GEN x, Z;
     402             :   long p, e;
     403       33906 :   switch(typ(z))
     404             :   {
     405       33906 :     case t_INT: case t_REAL: case t_FRAC: case t_COMPLEX: break;
     406           0 :     case t_QUAD:
     407           0 :       z = isexactzero(gel(z,2))? gel(z,1): quadtofp(z, T->prec);
     408           0 :       break;
     409           0 :     default: pari_err_TYPE("reduction mod 2-dim lattice (reduce_z)", z);
     410             :   }
     411       33906 :   Z = gdiv(z, T->W2);
     412       33906 :   T->z = z;
     413       33906 :   x = gdiv(imag_i(Z), imag_i(T->Tau));
     414       33906 :   T->x = grndtoi(x, &e); /* |Im(Z - x*Tau)| <= Im(Tau)/2 */
     415             :   /* Avoid Im(Z) << 0; take 0 <= Im(Z - x*Tau) < Im(Tau) instead.
     416             :    * Leave round when Im(Z - x*Tau) ~ 0 to allow detecting Z in <1,Tau>
     417             :    * at the end */
     418       33906 :   if (e > -10) T->x = gfloor(x);
     419       33906 :   if (signe(T->x)) Z = gsub(Z, gmul(T->x,T->Tau));
     420       33906 :   T->y = ground(real_i(Z));/* |Re(Z - y)| <= 1/2 */
     421       33906 :   if (signe(T->y)) Z = gsub(Z, T->y);
     422       33906 :   T->abs_u_is_1 = (typ(Z) != t_COMPLEX);
     423             :   /* Z = - y - x tau + z/W2, x,y integers */
     424       33906 :   check_complex(z, &(T->some_z_is_real), &(T->some_z_is_pure_imag));
     425       33906 :   if (!T->some_z_is_real && !T->some_z_is_pure_imag)
     426             :   {
     427             :     int W2real, W2imag;
     428       25404 :     check_complex(T->W2,&W2real,&W2imag);
     429       25404 :     if (W2real)
     430        3402 :       check_complex(Z, &(T->some_z_is_real), &(T->some_z_is_pure_imag));
     431       22002 :     else if (W2imag)
     432        4956 :       check_complex(Z, &(T->some_z_is_pure_imag), &(T->some_z_is_real));
     433             :   }
     434       33906 :   p = precision(Z);
     435       33906 :   if (gequal0(Z) || (p && gexpo(Z) < 5 - p)) Z = NULL; /*z in L*/
     436       33906 :   if (p && p < T->prec) T->prec = p;
     437       33906 :   T->Z = Z;
     438       33906 : }
     439             : /* return x.eta1 + y.eta2 */
     440             : static GEN
     441       32244 : eta_period(ellred_t *T, GEN eta)
     442             : {
     443       32244 :   GEN y1 = NULL, y2 = NULL;
     444       32244 :   if (signe(T->x)) y1 = gmul(T->x, gel(eta,1));
     445       32244 :   if (signe(T->y)) y2 = gmul(T->y, gel(eta,2));
     446       32244 :   if (!y1) return y2? y2: gen_0;
     447       12086 :   return y2? gadd(y1, y2): y1;
     448             : }
     449             : /* e is either
     450             :  * - [w1,w2]
     451             :  * - [[w1,w2],[eta1,eta2]]
     452             :  * - an ellinit structure */
     453             : static void
     454       62514 : compute_periods(ellred_t *T, GEN z, long prec)
     455             : {
     456             :   GEN w, e;
     457       62514 :   T->q_is_real = 0;
     458       62514 :   T->some_q_is_real = 0;
     459       62514 :   switch(T->type)
     460             :   {
     461       26304 :     case t_PER_ELL:
     462             :     {
     463       26304 :       long pr, p = prec;
     464       26304 :       if (z && (pr = precision(z))) p = pr;
     465       26304 :       e = T->in;
     466       26304 :       w = ellR_omega(e, p);
     467       26304 :       T->some_q_is_real = T->q_is_real = 1;
     468       26304 :       break;
     469             :     }
     470       11454 :     case t_PER_W:
     471       11454 :       w = T->in; break;
     472       24756 :     default: /*t_PER_WETA*/
     473       24756 :       w = gel(T->in,1); break;
     474             :   }
     475       62514 :   T->w1 = gel(w,1);
     476       62514 :   T->w2 = gel(w,2);
     477       62514 :   red_modSL2(T, prec);
     478       62514 :   if (z) reduce_z(z, T);
     479       62514 : }
     480             : static int
     481       62520 : check_periods(GEN e, ellred_t *T)
     482             : {
     483             :   GEN w1;
     484       62520 :   if (typ(e) != t_VEC) return 0;
     485       62520 :   T->in = e;
     486       62520 :   switch(lg(e))
     487             :   {
     488       26310 :     case 17:
     489       26310 :       T->type = t_PER_ELL;
     490       26310 :       break;
     491       36210 :     case 3:
     492       36210 :       w1 = gel(e,1);
     493       36210 :       if (typ(w1) != t_VEC)
     494       11454 :         T->type = t_PER_W;
     495             :       else
     496             :       {
     497       24756 :         if (lg(w1) != 3) return 0;
     498       24756 :         T->type = t_PER_WETA;
     499             :       }
     500       36210 :       break;
     501           0 :     default: return 0;
     502             :   }
     503       62520 :   return 1;
     504             : }
     505             : static int
     506       62448 : get_periods(GEN e, GEN z, ellred_t *T, long prec)
     507             : {
     508       62448 :   if (!check_periods(e, T)) return 0;
     509       62448 :   compute_periods(T, z, prec); return 1;
     510             : }
     511             : 
     512             : /* 2iPi/x, more efficient when x pure imaginary */
     513             : static GEN
     514      119412 : PiI2div(GEN x, long prec) { return gdiv(Pi2n(1, prec), mulcxmI(x)); }
     515             : /* (2iPi/W2)^k E_k(W1/W2) */
     516             : static GEN
     517       60852 : _elleisnum(ellred_t *T, long k)
     518             : {
     519       60852 :   GEN z = gmul(cxEk(T->Tau, k, T->prec), gpowgs(PiI2div(T->W2, T->prec), k));
     520       60852 :   return cxtoreal(z);
     521             : }
     522             : 
     523             : /* Return (2iPi)^k E_k(L) = (2iPi/w2)^k E_k(tau), with L = <w1,w2>, k > 0 even
     524             :  * E_k(tau) = 1 + 2/zeta(1-k) * sum(n>=1, n^(k-1) q^n/(1-q^n))
     525             :  * If flag is != 0 and k=4 or 6, compute g2 = E4/12 or g3 = -E6/216 resp. */
     526             : GEN
     527        3822 : elleisnum(GEN om, long k, long prec)
     528             : {
     529        3822 :   pari_sp av = avma;
     530             :   GEN y;
     531             :   ellred_t T;
     532             : 
     533        3822 :   if (k<=0) pari_err_DOMAIN("elleisnum", "k", "<=", gen_0, stoi(k));
     534        3822 :   if (k&1) pari_err_DOMAIN("elleisnum", "k % 2", "!=", gen_0, stoi(k));
     535        3822 :   if (!get_periods(om, NULL, &T, prec)) pari_err_TYPE("elleisnum",om);
     536        3822 :   y = _elleisnum(&T, k);
     537        3822 :   if (k==2 && signe(T.c))
     538             :   {
     539        3450 :     GEN a = gmul(Pi2n(1,T.prec), mului(12, T.c));
     540        3450 :     y = gsub(y, mulcxI(gdiv(a, gmul(T.w2, T.W2))));
     541             :   }
     542        3822 :   return gerepilecopy(av, gprec_wtrunc(y, T.prec0));
     543             : }
     544             : 
     545             : /* return quasi-periods attached to [T->W1,T->W2] */
     546             : static GEN
     547       56874 : _elleta(ellred_t *T)
     548             : {
     549       56874 :   GEN y1, y2, e2 = gdivgs(_elleisnum(T,2), -12);
     550       56874 :   y2 = gmul(T->W2, e2);
     551       56874 :   y1 = gsub(gmul(T->W1,e2), PiI2div(T->W2, T->prec));
     552       56874 :   retmkvec2(y1, y2);
     553             : }
     554             : 
     555             : /* compute eta1, eta2 */
     556             : GEN
     557          72 : elleta(GEN om, long prec)
     558             : {
     559          72 :   pari_sp av = avma;
     560             :   GEN y1, y2, E2, pi;
     561             :   ellred_t T;
     562             : 
     563          72 :   if (!check_periods(om, &T))
     564             :   {
     565           0 :     pari_err_TYPE("elleta",om);
     566             :     return NULL;/*LCOV_EXCL_LINE*/
     567             :   }
     568          72 :   if (T.type == t_PER_ELL) return ellR_eta(om, prec);
     569             : 
     570          66 :   compute_periods(&T, NULL, prec);
     571          66 :   prec = T.prec;
     572          66 :   pi = mppi(prec);
     573          66 :   E2 = cxEk(T.Tau, 2, prec); /* E_2(Tau) */
     574          66 :   if (signe(T.c))
     575             :   {
     576          18 :     GEN u = gdiv(T.w2, T.W2);
     577             :     /* E2 := u^2 E2 + 6iuc/pi = E_2(tau) */
     578          18 :     E2 = gadd(gmul(gsqr(u), E2), mulcxI(gdiv(gmul(mului(6,T.c), u), pi)));
     579             :   }
     580          66 :   y2 = gdiv(gmul(E2, sqrr(pi)), gmulsg(3, T.w2));
     581          66 :   if (T.swap)
     582             :   {
     583           6 :     y1 = y2;
     584           6 :     y2 = gadd(gmul(T.tau,y1), PiI2div(T.w2, prec));
     585             :   }
     586             :   else
     587          60 :     y1 = gsub(gmul(T.tau,y2), PiI2div(T.w2, prec));
     588          66 :   switch(typ(T.w1))
     589             :   {
     590          42 :     case t_INT: case t_FRAC: case t_REAL:
     591          42 :       y1 = real_i(y1);
     592             :   }
     593          66 :   return gerepilecopy(av, mkvec2(y1,y2));
     594             : }
     595             : GEN
     596       24642 : ellperiods(GEN w, long flag, long prec)
     597             : {
     598       24642 :   pari_sp av = avma;
     599             :   ellred_t T;
     600       24642 :   if (!get_periods(w, NULL, &T, prec)) pari_err_TYPE("ellperiods",w);
     601       24642 :   switch(flag)
     602             :   {
     603          12 :     case 0: return gerepilecopy(av, mkvec2(T.W1, T.W2));
     604       24630 :     case 1: return gerepilecopy(av, mkvec2(mkvec2(T.W1, T.W2), _elleta(&T)));
     605           0 :     default: pari_err_FLAG("ellperiods");
     606             :              return NULL;/*LCOV_EXCL_LINE*/
     607             :   }
     608             : }
     609             : 
     610             : /* 2Pi Im(z)/log(2) */
     611             : static double
     612       33798 : get_toadd(GEN z) { return (2*M_PI/M_LN2)*gtodouble(imag_i(z)); }
     613             : 
     614             : /* computes the numerical value of wp(z | L), L = om1 Z + om2 Z
     615             :  * return NULL if z in L.  If flall=1, compute also wp' */
     616             : static GEN
     617        1638 : ellwpnum_all(GEN e, GEN z, long flall, long prec)
     618             : {
     619             :   long toadd;
     620        1638 :   pari_sp av = avma, av1;
     621             :   GEN q, u, y, yp, u1, u2, qn;
     622             :   ellred_t T;
     623             :   int simple_case;
     624             : 
     625        1638 :   if (!get_periods(e, z, &T, prec)) pari_err_TYPE("ellwp",e);
     626        1638 :   if (!T.Z) return NULL;
     627        1620 :   prec = T.prec;
     628             : 
     629             :   /* Now L,Z normalized to <1,tau>. Z in fund. domain of <1, tau> */
     630        1620 :   q = expIPiC(gmul2n(T.Tau,1), prec);
     631        1620 :   u = expIPiC(gmul2n(T.Z,1), prec);
     632        1620 :   u1 = gsubsg(1,u);
     633        1620 :   u2 = gsqr(u1); /* (1-u)^2 = -4u sin^2(Pi Z) */
     634        1620 :   if (gequal0(gnorm(u2))) return NULL; /* possible if loss of accuracy */
     635        1620 :   y = gdiv(u,u2); /* -1/4(sin^2(Pi Z)) */
     636        1620 :   if (T.abs_u_is_1) y = real_i(y);
     637        1620 :   simple_case = T.abs_u_is_1 && T.q_is_real;
     638        1620 :   y = gadd(mkfrac(gen_1, utoipos(12)), y);
     639        1620 :   yp = flall? gen_0: NULL;
     640        1620 :   toadd = (long)ceil(get_toadd(T.Z));
     641             : 
     642        1620 :   av1 = avma; qn = q;
     643             :   for(;;)
     644       20819 :   { /* y += u q^n [ 1/(1-q^n u)^2 + 1/(q^n-u)^2 ] - 2q^n /(1-q^n)^2 */
     645             :     /* analogous formula for yp */
     646       22439 :     GEN yadd, ypadd = NULL;
     647       22439 :     GEN qnu = gmul(qn,u); /* q^n u */
     648       22439 :     GEN a = gsubsg(1,qnu);/* 1 - q^n u */
     649       22439 :     GEN a2 = gsqr(a);     /* (1 - q^n u)^2 */
     650       22439 :     if (yp) ypadd = gdiv(gaddsg(1,qnu),gmul(a,a2));
     651       22439 :     if (simple_case) /* conj(u) = 1/u: formula simplifies */
     652         355 :       yadd = gmul2n(real_i(gdiv(u,a2)), 1);
     653             :     else
     654             :     {
     655       22084 :       GEN b = gsub(qn,u);/* q^n - u */
     656       22084 :       GEN b2 = gsqr(b);  /* (q^n - u)^2 */
     657       22084 :       yadd = gmul(u, gadd(ginv(a2),ginv(b2)));
     658       22084 :       if (yp) ypadd = gadd(ypadd, gdiv(gadd(qn,u),gmul(b,b2)));
     659             :     }
     660       22439 :     yadd = gsub(yadd, gmul2n(ginv(gsqr(gsubsg(1,qn))), 1));
     661       22439 :     y = gadd(y, gmul(qn,yadd));
     662       22439 :     if (yp) yp = gadd(yp, gmul(qn,ypadd));
     663             : 
     664       22439 :     qn = gmul(q,qn);
     665       22439 :     if (gexpo(qn) <= - prec - 5 - toadd) break;
     666       20819 :     if (gc_needed(av1,1))
     667             :     {
     668           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"ellwp");
     669           0 :       gerepileall(av1, flall? 3: 2, &y, &qn, &yp);
     670             :     }
     671             :   }
     672        1620 :   if (yp)
     673             :   {
     674        1566 :     if (simple_case) yp = gsub(yp, conj_i(gmul(yp,gsqr(u))));
     675        1566 :     yp = gadd(yp, gdiv(gaddsg(1,u), gmul(u1,u2)));
     676             :   }
     677             : 
     678        1620 :   u1 = PiI2div(T.W2, prec);
     679        1620 :   u2 = gsqr(u1);
     680        1620 :   y = gmul(u2,y); /* y *= (2i pi / w2)^2 */
     681        1620 :   if (T.some_q_is_real && (T.some_z_is_real || T.some_z_is_pure_imag))
     682         882 :     y = real_i(y);
     683        1620 :   if (yp)
     684             :   {
     685        1566 :     yp = gmul(u, gmul(gmul(u1,u2),yp));/* yp *= u (2i pi / w2)^3 */
     686        1566 :     if (T.some_q_is_real)
     687             :     {
     688        1566 :       if (T.some_z_is_real) yp = real_i(yp);
     689         726 :       else if (T.some_z_is_pure_imag) yp = mkcomplex(gen_0, imag_i(yp));
     690             :     }
     691        1566 :     y = mkvec2(y, yp);
     692             :   }
     693        1620 :   return gerepilecopy(av, gprec_wtrunc(y, T.prec0));
     694             : }
     695             : static GEN
     696         258 : ellwpseries_aux(GEN c4, GEN c6, long v, long PRECDL)
     697             : {
     698             :   long i, k, l;
     699             :   pari_sp av;
     700         258 :   GEN _1, t, res = cgetg(PRECDL+2,t_SER), *P = (GEN*)(res + 2);
     701             : 
     702         258 :   res[1] = evalsigne(1) | _evalvalser(-2) | evalvarn(v);
     703         258 :   if (!PRECDL) { setsigne(res,0); return res; }
     704             : 
     705        2160 :   for (i=1; i<PRECDL; i+=2) P[i]= gen_0;
     706         258 :   _1 = Rg_get_1(c4);
     707         258 :   switch(PRECDL)
     708             :   {
     709         258 :     default:P[6] = gdivgu(c6,6048);
     710         258 :     case 6:
     711         258 :     case 5: P[4] = gdivgu(c4, 240);
     712         258 :     case 4:
     713         258 :     case 3: P[2] = gmul(_1,gen_0);
     714         258 :     case 2:
     715         258 :     case 1: P[0] = _1;
     716             :   }
     717         258 :   if (PRECDL <= 8) return res;
     718         258 :   av = avma;
     719         258 :   P[8] = gerepileupto(av, gdivgu(gsqr(P[4]), 3));
     720         930 :   for (k=5; (k<<1) < PRECDL; k++)
     721             :   {
     722         672 :     av = avma;
     723         672 :     t = gmul(P[4], P[(k-2)<<1]);
     724        1062 :     for (l=3; (l<<1) < k; l++) t = gadd(t, gmul(P[l<<1], P[(k-l)<<1]));
     725         672 :     t = gmul2n(t, 1);
     726         672 :     if ((k & 1) == 0) t = gadd(gsqr(P[k]), t);
     727         672 :     if (k % 3 == 2)
     728         234 :       t = gdivgu(gmulsg(3, t), (k-3)*(2*k+1));
     729             :     else /* same value, more efficient */
     730         438 :       t = gdivgu(t, ((k-3)*(2*k+1)) / 3);
     731         672 :     P[k<<1] = gerepileupto(av, t);
     732             :   }
     733         258 :   return res;
     734             : }
     735             : 
     736             : static int
     737         252 : get_c4c6(GEN w, GEN *c4, GEN *c6, long prec)
     738             : {
     739         252 :   if (typ(w) == t_VEC) switch(lg(w))
     740             :   {
     741         174 :     case 17:
     742         174 :       *c4 = ell_get_c4(w);
     743         174 :       *c6 = ell_get_c6(w);
     744         174 :       return 1;
     745          78 :     case 3:
     746             :     {
     747             :       ellred_t T;
     748          78 :       if (!get_periods(w,NULL,&T, prec)) break;
     749          78 :       *c4 = _elleisnum(&T, 4);
     750          78 :       *c6 = gneg(_elleisnum(&T, 6));
     751          78 :       return 1;
     752             :     }
     753             :   }
     754           0 :   *c4 = *c6 = NULL;
     755           0 :   return 0;
     756             : }
     757             : 
     758             : GEN
     759          12 : ellwpseries(GEN e, long v, long PRECDL)
     760             : {
     761             :   GEN c4, c6;
     762          12 :   checkell(e);
     763          12 :   c4 = ell_get_c4(e);
     764          12 :   c6 = ell_get_c6(e); return ellwpseries_aux(c4,c6,v,PRECDL);
     765             : }
     766             : 
     767             : GEN
     768           0 : ellwp(GEN w, GEN z, long prec)
     769           0 : { return ellwp0(w,z,0,prec); }
     770             : 
     771             : GEN
     772         156 : ellwp0(GEN w, GEN z, long flag, long prec)
     773             : {
     774         156 :   pari_sp av = avma;
     775             :   GEN y;
     776             : 
     777         156 :   if (flag && flag != 1) pari_err_FLAG("ellwp");
     778         156 :   if (!z) z = pol_x(0);
     779         156 :   y = toser_i(z);
     780         156 :   if (y)
     781             :   {
     782          90 :     long vy = varn(y), v = valser(y);
     783             :     GEN P, Q, c4,c6;
     784          90 :     if (!get_c4c6(w,&c4,&c6,prec)) pari_err_TYPE("ellwp",w);
     785          90 :     if (v <= 0) pari_err(e_IMPL,"ellwp(t_SER) away from 0");
     786          90 :     if (gequal0(y)) {
     787           0 :       set_avma(av);
     788           0 :       if (!flag) return zeroser(vy, -2*v);
     789           0 :       retmkvec2(zeroser(vy, -2*v), zeroser(vy, -3*v));
     790             :     }
     791          90 :     P = ellwpseries_aux(c4,c6, vy, lg(y)-2);
     792          90 :     Q = gsubst(P, varn(P), y);
     793          90 :     if (!flag)
     794          90 :       return gerepileupto(av, Q);
     795             :     else
     796             :     {
     797           0 :       GEN R = mkvec2(Q, gdiv(derivser(Q), derivser(y)));
     798           0 :       return gerepilecopy(av, R);
     799             :     }
     800             :   }
     801          66 :   y = ellwpnum_all(w,z,flag,prec);
     802          66 :   if (!y) pari_err_DOMAIN("ellwp", "argument","=", gen_0,z);
     803          60 :   return gerepileupto(av, y);
     804             : }
     805             : 
     806             : GEN
     807         138 : ellzeta(GEN w, GEN z, long prec0)
     808             : {
     809             :   long prec;
     810         138 :   pari_sp av = avma;
     811         138 :   GEN pi2, q, y, et = NULL;
     812             :   ellred_t T;
     813             : 
     814         138 :   if (!z) z = pol_x(0);
     815         138 :   y = toser_i(z);
     816         138 :   if (y)
     817             :   {
     818          78 :     long vy = varn(y), v = valser(y);
     819             :     GEN P, Q, c4,c6;
     820          78 :     if (!get_c4c6(w,&c4,&c6,prec0)) pari_err_TYPE("ellzeta",w);
     821          78 :     if (v <= 0) pari_err(e_IMPL,"ellzeta(t_SER) away from 0");
     822          78 :     if (gequal0(y)) { set_avma(av); return zeroser(vy, -v); }
     823          78 :     P = ellwpseries_aux(c4,c6, vy, lg(y)-2);
     824          78 :     P = integser(gneg(P)); /* \zeta' = - \wp*/
     825          78 :     Q = gsubst(P, varn(P), y);
     826          78 :     return gerepileupto(av, Q);
     827             :   }
     828          60 :   if (!get_periods(w, z, &T, prec0)) pari_err_TYPE("ellzeta", w);
     829          60 :   if (!T.Z) pari_err_DOMAIN("ellzeta", "z", "=", gen_0, z);
     830          60 :   prec = T.prec;
     831          60 :   if (signe(T.x) || signe(T.y)) et = eta_period(&T, _elleta(&T));
     832             : 
     833          60 :   pi2 = Pi2n(1, prec);
     834          60 :   q = expIPiC(gmul2n(T.Tau,1), prec);
     835          60 :   y = mulcxI(gmul(cxEk(T.Tau,2,prec), gmul(T.Z,divrs(pi2,-12))));
     836          60 :   if (!T.abs_u_is_1 || (!gequal(T.Z,ghalf) && !gequal(T.Z,gneg(ghalf))))
     837             :   { /* else u = -1 and this vanishes */
     838          60 :     long toadd = (long)ceil(get_toadd(T.Z));
     839          60 :     GEN qn, u, v, S = gen_0;
     840             :     pari_sp av1;
     841          60 :     u = expIPiC(gmul2n(T.Z,1), prec);
     842          60 :     v = gadd(ghalf, ginv(gsubgs(u, 1)));
     843          60 :     if (T.abs_u_is_1) gel(v,1) = gen_0; /*v = (u+1)/2(u-1), pure imaginary*/
     844          60 :     y = gadd(y, v);
     845             :     /* add sum_n q^n ( u/(u*q^n - 1) + 1/(u - q^n) )
     846             :      *     = (u^2 - 1) sum_n q^n / (uq^n - 1)(u - q^n) */
     847          60 :     av1 = avma;
     848          60 :     for (qn = q;;)
     849             :     {
     850         738 :       S = gadd(S, gdiv(qn, gmul(gsubgs(gmul(qn,u),1), gsub(u,qn))));
     851         738 :       qn = gmul(q,qn);
     852         738 :       if (gexpo(qn) <= - prec - 5 - toadd) break;
     853         678 :       if (gc_needed(av1,1))
     854             :       {
     855           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"ellzeta");
     856           0 :         gerepileall(av1,2, &S,&qn);
     857             :       }
     858             :     }
     859          60 :     y = gadd(y, gmul(gsubgs(gsqr(u),1), S));
     860             :   }
     861          60 :   y = mulcxI(gmul(gdiv(pi2,T.W2), y));
     862          60 :   if (T.some_q_is_real)
     863             :   {
     864          60 :     if (T.some_z_is_real)
     865             :     {
     866          24 :       if (!et || typ(et) != t_COMPLEX) y = real_i(y);
     867             :     }
     868          36 :     else if (T.some_z_is_pure_imag)
     869             :     {
     870          18 :       if (!et || (typ(et) == t_COMPLEX && isintzero(gel(et,1))))
     871          18 :         gel(y,1) = gen_0;
     872             :     }
     873             :   }
     874          60 :   if (et) y = gadd(y, et);
     875          60 :   return gerepilecopy(av, gprec_wtrunc(y, T.prec0));
     876             : }
     877             : 
     878             : /* if flag=0, return ellsigma, otherwise return log(ellsigma) */
     879             : GEN
     880       32292 : ellsigma(GEN w, GEN z, long flag, long prec0)
     881             : {
     882             :   long toadd, prec, n;
     883       32292 :   pari_sp av = avma, av1;
     884             :   GEN u, urn, urninv, z0, pi, pi2, q, q8, qn2, qn, y, y1, uinv, et, etnew;
     885             :   ellred_t T;
     886             : 
     887       32292 :   if (flag < 0 || flag > 1) pari_err_FLAG("ellsigma");
     888             : 
     889       32292 :   if (!z) z = pol_x(0);
     890       32292 :   y = toser_i(z);
     891       32292 :   if (y)
     892             :   {
     893          84 :     long vy = varn(y), v = valser(y);
     894             :     GEN P, Q, c4,c6;
     895          84 :     if (!get_c4c6(w,&c4,&c6,prec0)) pari_err_TYPE("ellsigma",w);
     896          84 :     if (v <= 0) pari_err_IMPL("ellsigma(t_SER) away from 0");
     897          84 :     if (flag) pari_err_TYPE("log(ellsigma)",y);
     898          78 :     if (gequal0(y)) { set_avma(av); return zeroser(vy, -v); }
     899          78 :     P = ellwpseries_aux(c4,c6, vy, lg(y)-2);
     900          78 :     P = integser(gneg(P)); /* \zeta' = - \wp*/
     901             :     /* (log \sigma)' = \zeta; remove log-singularity first */
     902          78 :     P = integser(serchop0(P));
     903          78 :     P = gexp(P, prec0);
     904          78 :     setvalser(P, valser(P)+1);
     905          78 :     Q = gsubst(P, varn(P), y);
     906          78 :     return gerepileupto(av, Q);
     907             :   }
     908       32208 :   if (!get_periods(w, z, &T, prec0)) pari_err_TYPE("ellsigma",w);
     909       32208 :   if (!T.Z)
     910             :   {
     911           6 :     if (!flag) return gen_0;
     912           6 :     pari_err_DOMAIN("log(ellsigma)", "argument","=",gen_0,z);
     913             :   }
     914       32202 :   prec = T.prec;
     915       32202 :   pi2 = Pi2n(1,prec);
     916       32202 :   pi  = mppi(prec);
     917             : 
     918       32202 :   urninv = uinv = NULL;
     919       32202 :   if (typ(T.Z) == t_FRAC && equaliu(gel(T.Z,2), 2) && equalim1(gel(T.Z,1)))
     920             :   {
     921          84 :     toadd = 0;
     922          84 :     urn = mkcomplex(gen_0, gen_m1); /* Z = -1/2 => urn = -I */
     923          84 :     u = gen_1;
     924             :   }
     925             :   else
     926             :   {
     927       32118 :     toadd = (long)ceil(fabs( get_toadd(T.Z) ));
     928       32118 :     urn = expIPiC(T.Z, prec); /* exp(i Pi Z) */
     929       32118 :     u = gneg_i(gsqr(urn));
     930       32118 :     if (!T.abs_u_is_1) { urninv = ginv(urn); uinv = gneg_i(gsqr(urninv)); }
     931             :   }
     932       32202 :   q8 = expIPiC(gmul2n(T.Tau, -2), prec);
     933       32202 :   q = gpowgs(q8,8); av1 = avma;
     934       32202 :   y = gen_0; qn = q; qn2 = gen_1;
     935      204462 :   for(n=0;;n++)
     936             :   { /* qn = q^(n+1), qn2 = q^(n(n+1)/2), urn = u^((n+1)/2)
     937             :      * if |u| = 1, will multiply by 2*I at the end ! */
     938      204462 :     y = gadd(y, gmul(qn2, uinv? gsub(urn,urninv): imag_i(urn)));
     939      204462 :     qn2 = gmul(qn,qn2);
     940      204462 :     if (gexpo(qn2) + n*toadd <= - prec - 5) break;
     941      172260 :     qn  = gmul(q,qn);
     942      172260 :     urn = gmul(urn,u);
     943      172260 :     if (uinv) urninv = gmul(urninv,uinv);
     944      172260 :     if (gc_needed(av1,1))
     945             :     {
     946           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"ellsigma");
     947           0 :       gerepileall(av1,urninv? 5: 4, &y,&qn,&qn2,&urn,&urninv);
     948             :     }
     949             :   }
     950       32202 :   y = gmul(y, gdiv(q8, gmul(pi2, gpowgs(trueeta(T.Tau,prec),3))));
     951       32202 :   y = gmul(y, T.abs_u_is_1? gmul2n(T.W2,1): mulcxmI(T.W2));
     952             : 
     953       32202 :   et = _elleta(&T);
     954       32202 :   z0 = gmul(T.Z,T.W2);
     955       32202 :   y1 = gadd(z0, gmul2n(gadd(gmul(T.x,T.W1), gmul(T.y,T.W2)),-1));
     956       32202 :   etnew = gmul(eta_period(&T, et), y1);
     957       32202 :   y1 = gadd(etnew, gmul2n(gmul(gmul(T.Z,z0),gel(et,2)),-1));
     958       32202 :   if (flag)
     959             :   {
     960       32142 :     y = gadd(y1, glog(y,prec));
     961       32142 :     if (mpodd(T.x) || mpodd(T.y)) y = gadd(y, mulcxI(pi));
     962             :     /* log(real number): im(y) = 0 or Pi */
     963       32142 :     if (T.some_q_is_real && isintzero(imag_i(z)) && gexpo(imag_i(y)) < 1)
     964           6 :       y = real_i(y);
     965             :   }
     966             :   else
     967             :   {
     968          60 :     y = gmul(y, gexp(y1,prec));
     969          60 :     if (mpodd(T.x) || mpodd(T.y)) y = gneg_i(y);
     970          60 :     if (T.some_q_is_real)
     971             :     {
     972             :       int re, cx;
     973          60 :       check_complex(z,&re,&cx);
     974          60 :       if (re) y = real_i(y);
     975          42 :       else if (cx && typ(y) == t_COMPLEX) gel(y,1) = gen_0;
     976             :     }
     977             :   }
     978       32202 :   return gerepilecopy(av, gprec_wtrunc(y, T.prec0));
     979             : }
     980             : 
     981             : GEN
     982        1620 : pointell(GEN e, GEN z, long prec)
     983             : {
     984        1620 :   pari_sp av = avma;
     985             :   GEN v;
     986             : 
     987        1620 :   checkell(e);
     988        1620 :   if (ell_get_type(e) == t_ELL_Qp)
     989             :   {
     990          48 :     prec = minss(ellQp_get_prec(e), padicprec_relative(z));
     991          48 :     return ellQp_t2P(e, z, prec);
     992             :   }
     993        1572 :   v = ellwpnum_all(e,z,1,prec);
     994        1572 :   if (!v) { set_avma(av); return ellinf(); }
     995        1560 :   gel(v,1) = gsub(gel(v,1), gdivgu(ell_get_b2(e),12));
     996        1560 :   gel(v,2) = gmul2n(gsub(gel(v,2), ec_h_evalx(e,gel(v,1))),-1);
     997        1560 :   return gerepilecopy(av, v);
     998             : }
     999             : 
    1000             : /********************************************************************/
    1001             : /**        exp(I*Pi*x) with attention for rational arguments        **/
    1002             : /********************************************************************/
    1003             : 
    1004             : /* sqrt(3)/2 */
    1005             : static GEN
    1006        1931 : sqrt32(long prec) { GEN z = sqrtr_abs(utor(3,prec)); setexpo(z, -1); return z; }
    1007             : /* exp(i k pi/12)  */
    1008             : static GEN
    1009        3197 : e12(ulong k, long prec)
    1010             : {
    1011             :   int s, sPi, sPiov2;
    1012             :   GEN z, t;
    1013        3197 :   k %= 24;
    1014        3197 :   if (!k) return gen_1;
    1015        3191 :   if (k == 12) return gen_m1;
    1016        3191 :   if (k >12) { s = 1; k = 24 - k; } else s = 0; /* x -> 2pi - x */
    1017        3191 :   if (k > 6) { sPi = 1; k = 12 - k; } else sPi = 0; /* x -> pi  - x */
    1018        3191 :   if (k > 3) { sPiov2 = 1; k = 6 - k; } else sPiov2 = 0; /* x -> pi/2 - x */
    1019        3191 :   z = cgetg(3, t_COMPLEX);
    1020        3191 :   switch(k)
    1021             :   {
    1022         553 :     case 0: gel(z,1) = icopy(gen_1); gel(z,2) = gen_0; break;
    1023         701 :     case 1: t = gmul2n(addrs(sqrt32(prec), 1), -1);
    1024         701 :       gel(z,1) = sqrtr(t);
    1025         701 :       gel(z,2) = gmul2n(invr(gel(z,1)), -2); break;
    1026             : 
    1027        1230 :     case 2: gel(z,1) = sqrt32(prec);
    1028        1230 :             gel(z,2) = real2n(-1, prec); break;
    1029             : 
    1030         707 :     case 3: gel(z,1) = sqrtr_abs(real2n(-1,prec));
    1031         707 :             gel(z,2) = rcopy(gel(z,1)); break;
    1032             :   }
    1033        3191 :   if (sPiov2) swap(gel(z,1), gel(z,2));
    1034        3191 :   if (sPi) togglesign(gel(z,1));
    1035        3191 :   if (s)   togglesign(gel(z,2));
    1036        3191 :   return z;
    1037             : }
    1038             : /* z a t_FRAC */
    1039             : static GEN
    1040       12946 : expIPifrac(GEN z, long prec)
    1041             : {
    1042       12946 :   GEN n = gel(z,1), d = gel(z,2);
    1043       12946 :   ulong r, q = uabsdivui_rem(12, d, &r);
    1044       12946 :   if (!r) return e12(q * umodiu(n, 24), prec); /* d | 12 */
    1045        9845 :   n = centermodii(n, shifti(d,1), d);
    1046        9845 :   return expIr(divri(mulri(mppi(prec), n), d));
    1047             : }
    1048             : /* exp(i Pi z), z a t_INT or t_FRAC */
    1049             : GEN
    1050       34062 : expIPiQ(GEN z, long prec)
    1051             : {
    1052       34062 :   if (typ(z) == t_INT) return mpodd(z)? gen_m1: gen_1;
    1053        1727 :   return expIPifrac(z, prec);
    1054             : }
    1055             : 
    1056             : /* convert power of 2 t_REAL to rational */
    1057             : static GEN
    1058        8093 : real2nQ(GEN x)
    1059             : {
    1060        8093 :   long e = expo(x);
    1061             :   GEN z;
    1062        8093 :   if (e < 0)
    1063        2477 :     z = mkfrac(signe(x) < 0? gen_m1: gen_1, int2n(-e));
    1064             :   else
    1065             :   {
    1066        5616 :     z = int2n(e);
    1067        5616 :     if (signe(x) < 0) togglesign_safe(&z);
    1068             :   }
    1069        8093 :   return z;
    1070             : }
    1071             : /* x a real number */
    1072             : GEN
    1073      158155 : expIPiR(GEN x, long prec)
    1074             : {
    1075      158155 :   if (typ(x) == t_REAL && absrnz_equal2n(x)) x = real2nQ(x);
    1076      158155 :   switch(typ(x))
    1077             :   {
    1078        2770 :     case t_INT:  return mpodd(x)? gen_m1: gen_1;
    1079        1691 :     case t_FRAC: return expIPifrac(x, prec);
    1080             :   }
    1081      153694 :   return expIr(mulrr(mppi(prec), x));
    1082             : }
    1083             : /* z a t_COMPLEX */
    1084             : GEN
    1085      323763 : expIPiC(GEN z, long prec)
    1086             : {
    1087             :   GEN pi, r, x, y;
    1088      323763 :   if (typ(z) != t_COMPLEX) return expIPiR(z, prec);
    1089      166503 :   x = gel(z,1);
    1090      166503 :   y = gel(z,2); if (gequal0(y)) return expIPiR(x, prec);
    1091      165608 :   pi = mppi(prec);
    1092      165608 :   r = gmul(pi, y); togglesign(r); r = mpexp(r); /* exp(-pi y) */
    1093      165608 :   if (typ(x) == t_REAL && absrnz_equal2n(x)) x = real2nQ(x);
    1094      165608 :   switch(typ(x))
    1095             :   {
    1096       54404 :     case t_INT: if (mpodd(x)) togglesign(r);
    1097       54404 :                 return r;
    1098        9528 :     case t_FRAC: return gmul(r, expIPifrac(x, prec));
    1099             :   }
    1100      101676 :   return gmul(r, expIr(mulrr(pi, x)));
    1101             : }
    1102             : /* exp(I x y), more efficient for x in R, y pure imaginary */
    1103             : GEN
    1104      510458 : expIxy(GEN x, GEN y, long prec) { return gexp(gmul(x, mulcxI(y)), prec); }
    1105             : 
    1106             : /********************************************************************/
    1107             : /**                 Eta function(s) and j-invariant                **/
    1108             : /********************************************************************/
    1109             : GEN
    1110       93303 : upper_to_cx(GEN x, long *prec)
    1111             : {
    1112       93303 :   long tx = typ(x), l;
    1113       93303 :   if (tx == t_QUAD) { x = quadtofp(x, *prec); tx = typ(x); }
    1114       93303 :   switch(tx)
    1115             :   {
    1116       93285 :     case t_COMPLEX:
    1117       93285 :       if (gsigne(gel(x,2)) > 0) break; /*fall through*/
    1118             :     case t_REAL: case t_INT: case t_FRAC:
    1119          12 :       pari_err_DOMAIN("modular function", "Im(argument)", "<=", gen_0, x);
    1120           6 :     default:
    1121           6 :       pari_err_TYPE("modular function", x);
    1122             :   }
    1123       93285 :   l = precision(x); if (l) *prec = l;
    1124       93285 :   return x;
    1125             : }
    1126             : 
    1127             : static GEN
    1128       60831 : qq(GEN x, long prec)
    1129             : {
    1130       60831 :   long tx = typ(x);
    1131             :   GEN y;
    1132             : 
    1133       60831 :   if (is_scalar_t(tx))
    1134             :   {
    1135       60795 :     if (tx == t_PADIC) return x;
    1136       60783 :     x = upper_to_cx(x, &prec);
    1137       60771 :     return cxtoreal(expIPiC(gmul2n(x,1), prec)); /* e(x) */
    1138             :   }
    1139          36 :   if (! ( y = toser_i(x)) ) pari_err_TYPE("modular function", x);
    1140          36 :   return y;
    1141             : }
    1142             : 
    1143             : /* return (y * X^d) + x. Assume d > 0, x != 0, valser(x) = 0 */
    1144             : static GEN
    1145          18 : ser_addmulXn(GEN y, GEN x, long d)
    1146             : {
    1147          18 :   long i, lx, ly, l = valser(y) + d; /* > 0 */
    1148             :   GEN z;
    1149             : 
    1150          18 :   lx = lg(x);
    1151          18 :   ly = lg(y) + l; if (lx < ly) ly = lx;
    1152          18 :   if (l > lx-2) return gcopy(x);
    1153          18 :   z = cgetg(ly,t_SER);
    1154          66 :   for (i=2; i<=l+1; i++) gel(z,i) = gel(x,i);
    1155          60 :   for (   ; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i-l));
    1156          18 :   z[1] = x[1]; return z;
    1157             : }
    1158             : 
    1159             : /* q a t_POL s.t. q(0) != 0, v > 0, Q = x^v*q; return \prod_i (1-Q^i) */
    1160             : static GEN
    1161          24 : RgXn_eta(GEN q, long v, long lim)
    1162             : {
    1163          24 :   pari_sp av = avma;
    1164             :   GEN qn, ps, y;
    1165             :   ulong vps, vqn, n;
    1166             : 
    1167          24 :   if (!degpol(q) && isint1(gel(q,2))) return eta_ZXn(v, lim+v);
    1168           6 :   y = qn = ps = pol_1(0);
    1169           6 :   vps = vqn = 0;
    1170           6 :   for(n = 0;; n++)
    1171           6 :   { /* qn = q^n,  ps = (-1)^n q^(n(3n+1)/2),
    1172             :      * vps, vqn valuation of ps, qn HERE */
    1173          12 :     pari_sp av2 = avma;
    1174          12 :     ulong vt = vps + 2*vqn + v; /* valuation of t at END of loop body */
    1175             :     long k1, k2;
    1176             :     GEN t;
    1177          12 :     vqn += v; vps = vt + vqn; /* valuation of qn, ps at END of body */
    1178          12 :     k1 = lim + v - vt + 1;
    1179          12 :     k2 = k1 - vqn; /* = lim + v - vps + 1 */
    1180          12 :     if (k1 <= 0) break;
    1181          12 :     t = RgX_mul(q, RgX_sqr(qn));
    1182          12 :     t = RgXn_red_shallow(t, k1);
    1183          12 :     t = RgX_mul(ps,t);
    1184          12 :     t = RgXn_red_shallow(t, k1);
    1185          12 :     t = RgX_neg(t); /* t = (-1)^(n+1) q^(n(3n+1)/2 + 2n+1) */
    1186          12 :     t = gerepileupto(av2, t);
    1187          12 :     y = RgX_addmulXn_shallow(t, y, vt);
    1188          12 :     if (k2 <= 0) break;
    1189             : 
    1190           6 :     qn = RgX_mul(qn,q);
    1191           6 :     ps = RgX_mul(t,qn);
    1192           6 :     ps = RgXn_red_shallow(ps, k2);
    1193           6 :     y = RgX_addmulXn_shallow(ps, y, vps);
    1194             : 
    1195           6 :     if (gc_needed(av,1))
    1196             :     {
    1197           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"eta, n = %ld", n);
    1198           0 :       gerepileall(av, 3, &y, &qn, &ps);
    1199             :     }
    1200             :   }
    1201           6 :   return y;
    1202             : }
    1203             : 
    1204             : static GEN
    1205       38725 : inteta(GEN q)
    1206             : {
    1207       38725 :   long tx = typ(q);
    1208             :   GEN ps, qn, y;
    1209             : 
    1210       38725 :   y = gen_1; qn = gen_1; ps = gen_1;
    1211       38725 :   if (tx==t_PADIC)
    1212             :   {
    1213          24 :     if (valp(q) <= 0) pari_err_DOMAIN("eta", "v_p(q)", "<=",gen_0,q);
    1214             :     for(;;)
    1215          48 :     {
    1216          66 :       GEN t = gneg_i(gmul(ps,gmul(q,gsqr(qn))));
    1217          66 :       y = gadd(y,t); qn = gmul(qn,q); ps = gmul(t,qn);
    1218          66 :       t = y;
    1219          66 :       y = gadd(y,ps); if (gequal(t,y)) return y;
    1220             :     }
    1221             :   }
    1222             : 
    1223       38701 :   if (tx == t_SER)
    1224             :   {
    1225             :     ulong vps, vqn;
    1226          36 :     long l = lg(q), v, n;
    1227             :     pari_sp av;
    1228             : 
    1229          36 :     v = valser(q); /* handle valuation separately to avoid overflow */
    1230          36 :     if (v <= 0) pari_err_DOMAIN("eta", "v_p(q)", "<=",gen_0,q);
    1231          30 :     y = ser2pol_i(q, l); /* t_SER inefficient when input has low degree */
    1232          30 :     n = degpol(y);
    1233          30 :     if (n <= (l>>2))
    1234             :     {
    1235          24 :       GEN z = RgXn_eta(y, v, l-2);
    1236          24 :       setvarn(z, varn(y)); return RgX_to_ser(z, l+v);
    1237             :     }
    1238           6 :     q = leafcopy(q); av = avma;
    1239           6 :     setvalser(q, 0);
    1240           6 :     y = scalarser(gen_1, varn(q), l+v);
    1241           6 :     vps = vqn = 0;
    1242           6 :     for(n = 0;; n++)
    1243           6 :     { /* qn = q^n,  ps = (-1)^n q^(n(3n+1)/2) */
    1244          12 :       ulong vt = vps + 2*vqn + v;
    1245             :       long k;
    1246             :       GEN t;
    1247          12 :       t = gneg_i(gmul(ps,gmul(q,gsqr(qn))));
    1248             :       /* t = (-1)^(n+1) q^(n(3n+1)/2 + 2n+1) */
    1249          12 :       y = ser_addmulXn(t, y, vt);
    1250          12 :       vqn += v; vps = vt + vqn;
    1251          12 :       k = l+v - vps; if (k <= 2) return y;
    1252             : 
    1253           6 :       qn = gmul(qn,q); ps = gmul(t,qn);
    1254           6 :       y = ser_addmulXn(ps, y, vps);
    1255           6 :       setlg(q, k);
    1256           6 :       setlg(qn, k);
    1257           6 :       setlg(ps, k);
    1258           6 :       if (gc_needed(av,3))
    1259             :       {
    1260           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"eta");
    1261           0 :         gerepileall(av, 3, &y, &qn, &ps);
    1262             :       }
    1263             :     }
    1264             :   }
    1265             :   {
    1266       38665 :     long l = -prec2nbits(precision(q));
    1267       38665 :     pari_sp av = avma;
    1268             : 
    1269             :     for(;;)
    1270      109992 :     {
    1271      148657 :       GEN t = gneg_i(gmul(ps,gmul(q,gsqr(qn))));
    1272             :       /* qn = q^n
    1273             :        * ps = (-1)^n q^(n(3n+1)/2)
    1274             :        * t = (-1)^(n+1) q^(n(3n+1)/2 + 2n+1) */
    1275      148657 :       y = gadd(y,t); qn = gmul(qn,q); ps = gmul(t,qn);
    1276      148657 :       y = gadd(y,ps);
    1277      148657 :       if (gexpo(ps)-gexpo(y) < l) return y;
    1278      109992 :       if (gc_needed(av,3))
    1279             :       {
    1280           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"eta");
    1281           0 :         gerepileall(av, 3, &y, &qn, &ps);
    1282             :       }
    1283             :     }
    1284             :   }
    1285             : }
    1286             : 
    1287             : GEN
    1288          66 : eta(GEN x, long prec)
    1289             : {
    1290          66 :   pari_sp av = avma;
    1291          66 :   GEN z = inteta( qq(x,prec) );
    1292          42 :   if (typ(z) == t_SER) return gerepilecopy(av, z);
    1293          12 :   return gerepileupto(av, z);
    1294             : }
    1295             : 
    1296             : /* s(h,k) = sum(n = 1, k-1, (n/k)*(frac(h*n/k) - 1/2))
    1297             :  * Knuth's algorithm. h integer, k integer > 0, (h,k) = 1 */
    1298             : GEN
    1299        5400 : sumdedekind_coprime(GEN h, GEN k)
    1300             : {
    1301        5400 :   pari_sp av = avma;
    1302             :   GEN s2, s1, p, pp;
    1303             :   long s;
    1304        5400 :   if (lgefint(k) == 3 && uel(k,2) <= (2*(ulong)LONG_MAX) / 3)
    1305             :   {
    1306        5394 :     ulong kk = k[2], hh = umodiu(h, kk);
    1307             :     long s1, s2;
    1308             :     GEN v;
    1309        5394 :     if (signe(k) < 0) { k = negi(k); hh = Fl_neg(hh, kk); }
    1310        5394 :     v = u_sumdedekind_coprime(hh, kk);
    1311        5394 :     s1 = v[1]; s2 = v[2];
    1312        5394 :     return gerepileupto(av, gdiv(addis(mulis(k,s1), s2), muluu(12, kk)));
    1313             :   }
    1314           6 :   s = 1;
    1315           6 :   s1 = gen_0; p = gen_1; pp = gen_0;
    1316           6 :   s2 = h = modii(h, k);
    1317          30 :   while (signe(h)) {
    1318          24 :     GEN r, nexth, a = dvmdii(k, h, &nexth);
    1319          24 :     if (is_pm1(h)) s2 = s == 1? addii(s2, p): subii(s2, p);
    1320          24 :     s1 = s == 1? addii(s1, a): subii(s1, a);
    1321          24 :     s = -s;
    1322          24 :     k = h; h = nexth;
    1323          24 :     r = addii(mulii(a,p), pp); pp = p; p = r;
    1324             :   }
    1325             :   /* at this point p = original k */
    1326           6 :   if (s == -1) s1 = subiu(s1, 3);
    1327           6 :   return gerepileupto(av, gdiv(addii(mulii(p,s1), s2), muliu(p,12)));
    1328             : }
    1329             : /* as above, for ulong arguments.
    1330             :  * k integer > 0, 0 <= h < k, (h,k) = 1. Returns [s1,s2] such that
    1331             :  * s(h,k) = (s2 + k s1) / (12k). Requires max(h + k/2, k) < LONG_MAX
    1332             :  * to avoid overflow, in particular k <= LONG_MAX * 2/3 is fine */
    1333             : GEN
    1334        5394 : u_sumdedekind_coprime(long h, long k)
    1335             : {
    1336        5394 :   long s = 1, s1 = 0, s2 = h, p = 1, pp = 0;
    1337        9828 :   while (h) {
    1338        4434 :     long r, nexth = k % h, a = k / h; /* a >= 1, a >= 2 if h = 1 */
    1339        4434 :     if (h == 1) s2 += p * s; /* occurs exactly once, last step */
    1340        4434 :     s1 += a * s;
    1341        4434 :     s = -s;
    1342        4434 :     k = h; h = nexth;
    1343        4434 :     r = a*p + pp; pp = p; p = r; /* p >= pp >= 0 */
    1344             :   }
    1345             :   /* in the above computation, p increases from 1 to original k,
    1346             :    * -k/2 <= s2 <= h + k/2, and |s1| <= k */
    1347        5394 :   if (s < 0) s1 -= 3; /* |s1| <= k+3 ? */
    1348             :   /* But in fact, |s2 + p s1| <= k^2 + 1/2 - 3k; if (s < 0), we have
    1349             :    * |s2| <= k/2 and it follows that |s1| < k here as well */
    1350             :   /* p = k; s(h,k) = (s2 + p s1)/12p. */
    1351        5394 :   return mkvecsmall2(s1, s2);
    1352             : }
    1353             : GEN
    1354          24 : sumdedekind(GEN h, GEN k)
    1355             : {
    1356          24 :   pari_sp av = avma;
    1357             :   GEN d;
    1358          24 :   if (typ(h) != t_INT) pari_err_TYPE("sumdedekind",h);
    1359          24 :   if (typ(k) != t_INT) pari_err_TYPE("sumdedekind",k);
    1360          24 :   d = gcdii(h,k);
    1361          24 :   if (!is_pm1(d))
    1362             :   {
    1363           6 :     h = diviiexact(h, d);
    1364           6 :     k = diviiexact(k, d);
    1365             :   }
    1366          24 :   return gerepileupto(av, sumdedekind_coprime(h,k));
    1367             : }
    1368             : 
    1369             : /* eta(x); assume Im x >> 0 (e.g. x in SL2's standard fundamental domain) */
    1370             : static GEN
    1371       39468 : eta_reduced(GEN x, long prec)
    1372             : {
    1373       39468 :   GEN z = expIPiC(gdivgu(x, 12), prec); /* e(x/24) */
    1374       39468 :   if (24 * gexpo(z) >= -prec2nbits(prec))
    1375       38647 :     z = gmul(z, inteta( gpowgs(z,24) ));
    1376       39468 :   return z;
    1377             : }
    1378             : 
    1379             : /* x = U.z (flag = 1), or x = U^(-1).z (flag = 0)
    1380             :  * Return [s,t] such that eta(z) = eta(x) * sqrt(s) * exp(I Pi t) */
    1381             : static GEN
    1382       39480 : eta_correction(GEN x, GEN U, long flag)
    1383             : {
    1384             :   GEN a,b,c,d, s,t;
    1385             :   long sc;
    1386       39480 :   a = gcoeff(U,1,1);
    1387       39480 :   b = gcoeff(U,1,2);
    1388       39480 :   c = gcoeff(U,2,1);
    1389       39480 :   d = gcoeff(U,2,2);
    1390             :   /* replace U by U^(-1) */
    1391       39480 :   if (flag) {
    1392       32304 :     swap(a,d);
    1393       32304 :     togglesign_safe(&b);
    1394       32304 :     togglesign_safe(&c);
    1395             :   }
    1396       39480 :   sc = signe(c);
    1397       39480 :   if (!sc) {
    1398       34104 :     if (signe(d) < 0) togglesign_safe(&b);
    1399       34104 :     s = gen_1;
    1400       34104 :     t = uutoQ(umodiu(b, 24), 12);
    1401             :   } else {
    1402        5376 :     if (sc < 0) {
    1403        1512 :       togglesign_safe(&a);
    1404        1512 :       togglesign_safe(&b);
    1405        1512 :       togglesign_safe(&c);
    1406        1512 :       togglesign_safe(&d);
    1407             :     } /* now c > 0 */
    1408        5376 :     s = mulcxmI(gadd(gmul(c,x), d));
    1409        5376 :     t = gadd(gdiv(addii(a,d),muliu(c,12)), sumdedekind_coprime(negi(d),c));
    1410             :     /* correction : exp(I Pi (((a+d)/12c) + s(-d,c)) ) sqrt(-i(cx+d))  */
    1411             :   }
    1412       39480 :   return mkvec2(s, t);
    1413             : }
    1414             : 
    1415             : /* returns the true value of eta(x) for Im(x) > 0, using reduction to
    1416             :  * standard fundamental domain */
    1417             : GEN
    1418       32232 : trueeta(GEN x, long prec)
    1419             : {
    1420       32232 :   pari_sp av = avma;
    1421             :   GEN U, st, s, t;
    1422             : 
    1423       32232 :   if (!is_scalar_t(typ(x))) pari_err_TYPE("trueeta",x);
    1424       32232 :   x = upper_to_cx(x, &prec);
    1425       32232 :   x = cxredsl2(x, &U);
    1426       32232 :   st = eta_correction(x, U, 1);
    1427       32232 :   x = eta_reduced(x, prec);
    1428       32232 :   s = gel(st, 1);
    1429       32232 :   t = gel(st, 2);
    1430       32232 :   x = gmul(x, expIPiQ(t, prec));
    1431       32232 :   if (s != gen_1) x = gmul(x, gsqrt(s, prec));
    1432       32232 :   return gerepileupto(av, x);
    1433             : }
    1434             : 
    1435             : GEN
    1436          96 : eta0(GEN x, long flag,long prec)
    1437          96 : { return flag? trueeta(x,prec): eta(x,prec); }
    1438             : 
    1439             : /* eta(q) = 1 + \sum_{n>0} (-1)^n * (q^(n(3n-1)/2) + q^(n(3n+1)/2)) */
    1440             : static GEN
    1441           6 : ser_eta(long prec)
    1442             : {
    1443           6 :   GEN e = cgetg(prec+2, t_SER), ed = e+2;
    1444             :   long n, j;
    1445           6 :   e[1] = evalsigne(1)|_evalvalser(0)|evalvarn(0);
    1446           6 :   gel(ed,0) = gen_1;
    1447         414 :   for (n = 1; n < prec; n++) gel(ed,n) = gen_0;
    1448          42 :   for (n = 1, j = 0; n < prec; n++)
    1449             :   {
    1450             :     GEN s;
    1451          42 :     j += 3*n-2; /* = n*(3*n-1) / 2 */;
    1452          42 :     if (j >= prec) break;
    1453          36 :     s = odd(n)? gen_m1: gen_1;
    1454          36 :     gel(ed, j) = s;
    1455          36 :     if (j+n >= prec) break;
    1456          36 :     gel(ed, j+n) = s;
    1457             :   }
    1458           6 :   return e;
    1459             : }
    1460             : 
    1461             : static GEN
    1462         408 : coeffEu(GEN fa)
    1463             : {
    1464         408 :   pari_sp av = avma;
    1465         408 :   return gerepileuptoint(av, mului(65520, usumdivk_fact(fa,11)));
    1466             : }
    1467             : /* E12 = 1 + q*E/691 */
    1468             : static GEN
    1469           6 : ser_E(long prec)
    1470             : {
    1471           6 :   GEN e = cgetg(prec+2, t_SER), ed = e+2;
    1472           6 :   GEN F = vecfactoru_i(2, prec); /* F[n] = factoru(n+1) */
    1473             :   long n;
    1474           6 :   e[1] = evalsigne(1)|_evalvalser(0)|evalvarn(0);
    1475           6 :   gel(ed,0) = utoipos(65520);
    1476         414 :   for (n = 1; n < prec; n++) gel(ed,n) = coeffEu(gel(F,n));
    1477           6 :   return e;
    1478             : }
    1479             : /* j = E12/Delta + 432000/691, E12 = 1 + q*E/691 */
    1480             : static GEN
    1481           6 : ser_j2(long prec, long v)
    1482             : {
    1483           6 :   pari_sp av = avma;
    1484           6 :   GEN iD = gpowgs(ginv(ser_eta(prec)), 24); /* q/Delta */
    1485           6 :   GEN J = gmul(ser_E(prec), iD);
    1486           6 :   setvalser(iD,-1); /* now 1/Delta */
    1487           6 :   J = gadd(gdivgu(J, 691), iD);
    1488           6 :   J = gerepileupto(av, J);
    1489           6 :   if (prec > 1) gel(J,3) = utoipos(744);
    1490           6 :   setvarn(J,v); return J;
    1491             : }
    1492             : 
    1493             : /* j(q) = \sum_{n >= -1} c(n)q^n,
    1494             :  * \sum_{n = -1}^{N-1} c(n) (-10n \sigma_3(N-n) + 21 \sigma_5(N-n))
    1495             :  * = c(N) (N+1)/24 */
    1496             : static GEN
    1497          12 : ser_j(long prec, long v)
    1498             : {
    1499             :   GEN j, J, S3, S5, F;
    1500             :   long i, n;
    1501          12 :   if (prec > 64) return ser_j2(prec, v);
    1502           6 :   S3 = cgetg(prec+1, t_VEC);
    1503           6 :   S5 = cgetg(prec+1,t_VEC);
    1504           6 :   F = vecfactoru_i(1, prec);
    1505          30 :   for (n = 1; n <= prec; n++)
    1506             :   {
    1507          24 :     GEN fa = gel(F,n);
    1508          24 :     gel(S3,n) = mului(10, usumdivk_fact(fa,3));
    1509          24 :     gel(S5,n) = mului(21, usumdivk_fact(fa,5));
    1510             :   }
    1511           6 :   J = cgetg(prec+2, t_SER),
    1512           6 :   J[1] = evalvarn(v)|evalsigne(1)|evalvalser(-1);
    1513           6 :   j = J+3;
    1514           6 :   gel(j,-1) = gen_1;
    1515           6 :   gel(j,0) = utoipos(744);
    1516           6 :   gel(j,1) = utoipos(196884);
    1517          18 :   for (n = 2; n < prec; n++)
    1518             :   {
    1519          12 :     pari_sp av = avma;
    1520          12 :     GEN c, s3 = gel(S3,n+1), s5 = gel(S5,n+1);
    1521          12 :     c = addii(s3, s5);
    1522          42 :     for (i = 0; i < n; i++)
    1523             :     {
    1524          30 :       s3 = gel(S3,n-i); s5 = gel(S5,n-i);
    1525          30 :       c = addii(c, mulii(gel(j,i), subii(s5, mului(i,s3))));
    1526             :     }
    1527          12 :     gel(j,n) = gerepileuptoint(av, diviuexact(muliu(c,24), n+1));
    1528             :   }
    1529           6 :   return J;
    1530             : }
    1531             : 
    1532             : GEN
    1533          36 : jell(GEN x, long prec)
    1534             : {
    1535          36 :   long tx = typ(x);
    1536          36 :   pari_sp av = avma;
    1537             :   GEN q, h, U;
    1538             : 
    1539          36 :   if (!is_scalar_t(tx))
    1540             :   {
    1541             :     long v;
    1542          18 :     if (gequalX(x)) return ser_j(precdl, varn(x));
    1543          18 :     q = toser_i(x); if (!q) pari_err_TYPE("ellj",x);
    1544          12 :     v = fetch_var_higher();
    1545          12 :     h = ser_j(lg(q)-2, v);
    1546          12 :     h = gsubst(h, v, q);
    1547          12 :     delete_var(); return gerepileupto(av, h);
    1548             :   }
    1549          18 :   if (tx == t_PADIC)
    1550             :   {
    1551           6 :     GEN p2, p1 = gdiv(inteta(gsqr(x)), inteta(x));
    1552           6 :     p1 = gmul2n(gsqr(p1),1);
    1553           6 :     p1 = gmul(x,gpowgs(p1,12));
    1554           6 :     p2 = gaddsg(768,gadd(gsqr(p1),gdivsg(4096,p1)));
    1555           6 :     p1 = gmulsg(48,p1);
    1556           6 :     return gerepileupto(av, gadd(p2,p1));
    1557             :   }
    1558             :   /* Let h = Delta(2x) / Delta(x), then j(x) = (1 + 256h)^3 / h */
    1559          12 :   x = upper_to_cx(x, &prec);
    1560           6 :   x = cxredsl2(x, &U); /* forget about Ua : j has weight 0 */
    1561             :   { /* cf eta_reduced, raised to power 24
    1562             :      * Compute
    1563             :      *   t = (inteta(q(2x)) / inteta(q(x))) ^ 24;
    1564             :      * then
    1565             :      *   h = t * (q(2x) / q(x) = t * q(x);
    1566             :      * but inteta(q) costly and useless if expo(q) << 1  => inteta(q) = 1.
    1567             :      * log_2 ( exp(-2Pi Im tau) ) < -prec2nbits(prec)
    1568             :      * <=> Im tau > prec2nbits(prec) * log(2) / 2Pi */
    1569           6 :     long C = (long)prec2nbits_mul(prec, M_LN2/(2*M_PI));
    1570           6 :     q = expIPiC(gmul2n(x,1), prec); /* e(x) */
    1571           6 :     if (gcmpgs(gel(x,2), C) > 0) /* eta(q(x)) = 1 : no need to compute q(2x) */
    1572           0 :       h = q;
    1573             :     else
    1574             :     {
    1575           6 :       GEN t = gdiv(inteta(gsqr(q)), inteta(q));
    1576           6 :       h = gmul(q, gpowgs(t, 24));
    1577             :     }
    1578             :   }
    1579             :   /* real_1 important ! gaddgs(, 1) could increase the accuracy ! */
    1580           6 :   return gerepileupto(av, gdiv(gpowgs(gadd(gmul2n(h,8), real_1(prec)), 3), h));
    1581             : }
    1582             : 
    1583             : static GEN
    1584        7176 : to_form(GEN a, GEN w, GEN C, GEN D)
    1585        7176 : { return mkqfb(a, w, diviiexact(C, a), D); }
    1586             : static GEN
    1587        7176 : form_to_quad(GEN f, GEN sqrtD)
    1588             : {
    1589        7176 :   long a = itos(gel(f,1)), a2 = a << 1;
    1590        7176 :   GEN b = gel(f,2);
    1591        7176 :   return mkcomplex(gdivgs(b, -a2), gdivgs(sqrtD, a2));
    1592             : }
    1593             : static GEN
    1594        7176 : eta_form(GEN f, GEN sqrtD, GEN *s_t, long prec)
    1595             : {
    1596        7176 :   GEN U, t = form_to_quad(redimagsl2(f, &U), sqrtD);
    1597        7176 :   *s_t = eta_correction(t, U, 0);
    1598        7176 :   return eta_reduced(t, prec);
    1599             : }
    1600             : 
    1601             : /* eta(t/p)eta(t/q) / (eta(t)eta(t/pq)), t = (-w + sqrt(D)) / 2a */
    1602             : GEN
    1603        1794 : double_eta_quotient(GEN a, GEN w, GEN D, long p, long q, GEN pq, GEN sqrtD)
    1604             : {
    1605        1794 :   GEN C = shifti(subii(sqri(w), D), -2);
    1606             :   GEN d, t, z, zp, zq, zpq, s_t, s_tp, s_tpq, s, sp, spq;
    1607        1794 :   long prec = realprec(sqrtD);
    1608             : 
    1609        1794 :   z = eta_form(to_form(a, w, C, D), sqrtD, &s_t, prec);
    1610        1794 :   s = gel(s_t, 1);
    1611        1794 :   zp = eta_form(to_form(mului(p, a), w, C, D), sqrtD, &s_tp, prec);
    1612        1794 :   sp = gel(s_tp, 1);
    1613        1794 :   zpq = eta_form(to_form(mulii(pq, a), w, C, D), sqrtD, &s_tpq, prec);
    1614        1794 :   spq = gel(s_tpq, 1);
    1615        1794 :   if (p == q) {
    1616           0 :     z = gdiv(gsqr(zp), gmul(z, zpq));
    1617           0 :     t = gsub(gmul2n(gel(s_tp,2), 1),
    1618           0 :              gadd(gel(s_t,2), gel(s_tpq,2)));
    1619           0 :     if (sp != gen_1) z = gmul(z, sp);
    1620             :   } else {
    1621             :     GEN s_tq, sq;
    1622        1794 :     zq = eta_form(to_form(mului(q, a), w, C, D), sqrtD, &s_tq, prec);
    1623        1794 :     sq = gel(s_tq, 1);
    1624        1794 :     z = gdiv(gmul(zp, zq), gmul(z, zpq));
    1625        1794 :     t = gsub(gadd(gel(s_tp,2), gel(s_tq,2)),
    1626        1794 :              gadd(gel(s_t,2), gel(s_tpq,2)));
    1627        1794 :     if (sp != gen_1) z = gmul(z, gsqrt(sp, prec));
    1628        1794 :     if (sq != gen_1) z = gmul(z, gsqrt(sq, prec));
    1629             :   }
    1630        1794 :   d = NULL;
    1631        1794 :   if (s != gen_1) d = gsqrt(s, prec);
    1632        1794 :   if (spq != gen_1) {
    1633        1770 :     GEN x = gsqrt(spq, prec);
    1634        1770 :     d = d? gmul(d, x): x;
    1635             :   }
    1636        1794 :   if (d) z = gdiv(z, d);
    1637        1794 :   return gmul(z, expIPiQ(t, prec));
    1638             : }
    1639             : 
    1640             : typedef struct { GEN u; long v, t; } cxanalyze_t;
    1641             : 
    1642             : /* Check whether a t_COMPLEX, t_REAL or t_INT z != 0 can be written as
    1643             :  * z = u * 2^(v/2) * exp(I Pi/4 t), u > 0, v = 0,1 and -3 <= t <= 4.
    1644             :  * Allow z t_INT/t_REAL to simplify handling of eta_correction() output */
    1645             : static int
    1646          72 : cxanalyze(cxanalyze_t *T, GEN z)
    1647             : {
    1648             :   GEN a, b;
    1649             :   long ta, tb;
    1650             : 
    1651          72 :   T->u = z;
    1652          72 :   T->v = 0;
    1653          72 :   if (is_intreal_t(typ(z)))
    1654             :   {
    1655          60 :     T->u = mpabs_shallow(z);
    1656          60 :     T->t = signe(z) < 0? 4: 0;
    1657          60 :     return 1;
    1658             :   }
    1659          12 :   a = gel(z,1); ta = typ(a);
    1660          12 :   b = gel(z,2); tb = typ(b);
    1661             : 
    1662          12 :   T->t = 0;
    1663          12 :   if (ta == t_INT && !signe(a))
    1664             :   {
    1665           0 :     T->u = R_abs_shallow(b);
    1666           0 :     T->t = gsigne(b) < 0? -2: 2;
    1667           0 :     return 1;
    1668             :   }
    1669          12 :   if (tb == t_INT && !signe(b))
    1670             :   {
    1671           0 :     T->u = R_abs_shallow(a);
    1672           0 :     T->t = gsigne(a) < 0? 4: 0;
    1673           0 :     return 1;
    1674             :   }
    1675          12 :   if (ta != tb || ta == t_REAL) return 0;
    1676             :   /* a,b both non zero, both t_INT or t_FRAC */
    1677          12 :   if (ta == t_INT)
    1678             :   {
    1679           6 :     if (!absequalii(a, b)) return 0;
    1680           6 :     T->u = absi_shallow(a);
    1681           6 :     T->v = 1;
    1682           6 :     if (signe(a) == signe(b))
    1683           0 :     { T->t = signe(a) < 0? -3: 1; }
    1684             :     else
    1685           6 :     { T->t = signe(a) < 0? 3: -1; }
    1686             :   }
    1687             :   else
    1688             :   {
    1689           6 :     if (!absequalii(gel(a,2), gel(b,2)) || !absequalii(gel(a,1),gel(b,1)))
    1690           6 :       return 0;
    1691           0 :     T->u = absfrac_shallow(a);
    1692           0 :     T->v = 1;
    1693           0 :     a = gel(a,1);
    1694           0 :     b = gel(b,1);
    1695           0 :     if (signe(a) == signe(b))
    1696           0 :     { T->t = signe(a) < 0? -3: 1; }
    1697             :     else
    1698           0 :     { T->t = signe(a) < 0? 3: -1; }
    1699             :   }
    1700           6 :   return 1;
    1701             : }
    1702             : 
    1703             : /* z * sqrt(st_b) / sqrt(st_a) exp(I Pi (t + t0)). Assume that
    1704             :  * sqrt2 = gsqrt(gen_2, prec) or NULL */
    1705             : static GEN
    1706          36 : apply_eta_correction(GEN z, GEN st_a, GEN st_b, GEN t0, GEN sqrt2, long prec)
    1707             : {
    1708          36 :   GEN t, s_a = gel(st_a, 1), s_b = gel(st_b, 1);
    1709             :   cxanalyze_t Ta, Tb;
    1710             :   int ca, cb;
    1711             : 
    1712          36 :   t = gsub(gel(st_b,2), gel(st_a,2));
    1713          36 :   if (t0 != gen_0) t = gadd(t, t0);
    1714          36 :   ca = cxanalyze(&Ta, s_a);
    1715          36 :   cb = cxanalyze(&Tb, s_b);
    1716          36 :   if (ca || cb)
    1717          36 :   { /* compute sqrt(s_b) / sqrt(s_a) in a more efficient way:
    1718             :      * sb = ub sqrt(2)^vb exp(i Pi/4 tb) */
    1719          36 :     GEN u = gdiv(Tb.u,Ta.u);
    1720          36 :     switch(Tb.v - Ta.v)
    1721             :     {
    1722           0 :       case -1: u = gmul2n(u,-1); /* fall through: write 1/sqrt2 = sqrt2/2 */
    1723           6 :       case 1: u = gmul(u, sqrt2? sqrt2: sqrtr_abs(real2n(1, prec)));
    1724             :     }
    1725          36 :     if (!isint1(u)) z = gmul(z, gsqrt(u, prec));
    1726          36 :     t = gadd(t, gmul2n(stoi(Tb.t - Ta.t), -3));
    1727             :   }
    1728             :   else
    1729             :   {
    1730           0 :     z = gmul(z, gsqrt(s_b, prec));
    1731           0 :     z = gdiv(z, gsqrt(s_a, prec));
    1732             :   }
    1733          36 :   return gmul(z, expIPiQ(t, prec));
    1734             : }
    1735             : 
    1736             : /* sqrt(2) eta(2x) / eta(x) */
    1737             : GEN
    1738          12 : weberf2(GEN x, long prec)
    1739             : {
    1740          12 :   pari_sp av = avma;
    1741             :   GEN z, sqrt2, a,b, Ua,Ub, st_a,st_b;
    1742             : 
    1743          12 :   x = upper_to_cx(x, &prec);
    1744          12 :   a = cxredsl2(x, &Ua);
    1745          12 :   b = cxredsl2(gmul2n(x,1), &Ub);
    1746          12 :   if (gequal(a,b)) /* not infrequent */
    1747           0 :     z = gen_1;
    1748             :   else
    1749          12 :     z = gdiv(eta_reduced(b,prec), eta_reduced(a,prec));
    1750          12 :   st_a = eta_correction(a, Ua, 1);
    1751          12 :   st_b = eta_correction(b, Ub, 1);
    1752          12 :   sqrt2 = sqrtr_abs(real2n(1, prec));
    1753          12 :   z = apply_eta_correction(z, st_a, st_b, gen_0, sqrt2, prec);
    1754          12 :   return gerepileupto(av, gmul(z, sqrt2));
    1755             : }
    1756             : 
    1757             : /* eta(x/2) / eta(x) */
    1758             : GEN
    1759          12 : weberf1(GEN x, long prec)
    1760             : {
    1761          12 :   pari_sp av = avma;
    1762             :   GEN z, a,b, Ua,Ub, st_a,st_b;
    1763             : 
    1764          12 :   x = upper_to_cx(x, &prec);
    1765          12 :   a = cxredsl2(x, &Ua);
    1766          12 :   b = cxredsl2(gmul2n(x,-1), &Ub);
    1767          12 :   if (gequal(a,b)) /* not infrequent */
    1768           0 :     z = gen_1;
    1769             :   else
    1770          12 :     z = gdiv(eta_reduced(b,prec), eta_reduced(a,prec));
    1771          12 :   st_a = eta_correction(a, Ua, 1);
    1772          12 :   st_b = eta_correction(b, Ub, 1);
    1773          12 :   z = apply_eta_correction(z, st_a, st_b, gen_0, NULL, prec);
    1774          12 :   return gerepileupto(av, z);
    1775             : }
    1776             : /* exp(-I*Pi/24) * eta((x+1)/2) / eta(x) */
    1777             : GEN
    1778          12 : weberf(GEN x, long prec)
    1779             : {
    1780          12 :   pari_sp av = avma;
    1781             :   GEN z, t0, a,b, Ua,Ub, st_a,st_b;
    1782          12 :   x = upper_to_cx(x, &prec);
    1783          12 :   a = cxredsl2(x, &Ua);
    1784          12 :   b = cxredsl2(gmul2n(gaddgs(x,1),-1), &Ub);
    1785          12 :   if (gequal(a,b)) /* not infrequent */
    1786           6 :     z = gen_1;
    1787             :   else
    1788           6 :     z = gdiv(eta_reduced(b,prec), eta_reduced(a,prec));
    1789          12 :   st_a = eta_correction(a, Ua, 1);
    1790          12 :   st_b = eta_correction(b, Ub, 1);
    1791          12 :   t0 = mkfrac(gen_m1, utoipos(24));
    1792          12 :   z = apply_eta_correction(z, st_a, st_b, t0, NULL, prec);
    1793          12 :   if (typ(z) == t_COMPLEX && isexactzero(real_i(x)))
    1794           0 :     z = gerepilecopy(av, gel(z,1));
    1795             :   else
    1796          12 :     z = gerepileupto(av, z);
    1797          12 :   return z;
    1798             : }
    1799             : GEN
    1800          36 : weber0(GEN x, long flag,long prec)
    1801             : {
    1802          36 :   switch(flag)
    1803             :   {
    1804          12 :     case 0: return weberf(x,prec);
    1805          12 :     case 1: return weberf1(x,prec);
    1806          12 :     case 2: return weberf2(x,prec);
    1807           0 :     default: pari_err_FLAG("weber");
    1808             :   }
    1809             :   return NULL; /* LCOV_EXCL_LINE */
    1810             : }
    1811             : 
    1812             : /********************************************************************/
    1813             : /**                     Jacobi sine theta                          **/
    1814             : /********************************************************************/
    1815             : 
    1816             : /* check |q| < 1 */
    1817             : static GEN
    1818          18 : check_unit_disc(const char *fun, GEN q, long prec)
    1819             : {
    1820          18 :   GEN Q = gtofp(q, prec), Qlow;
    1821          18 :   Qlow = (prec > LOWDEFAULTPREC)? gtofp(Q,LOWDEFAULTPREC): Q;
    1822          18 :   if (gcmp(gnorm(Qlow), gen_1) >= 0)
    1823           0 :     pari_err_DOMAIN(fun, "abs(q)", ">=", gen_1, q);
    1824          18 :   return Q;
    1825             : }
    1826             : 
    1827             : GEN
    1828           6 : thetanullk(GEN q, long k, long prec)
    1829             : {
    1830             :   long l, n;
    1831           6 :   pari_sp av = avma;
    1832             :   GEN p1, ps, qn, y, ps2;
    1833             : 
    1834           6 :   if (k < 0)
    1835           0 :     pari_err_DOMAIN("thetanullk", "k", "<", gen_0, stoi(k));
    1836           6 :   l = precision(q);
    1837           6 :   if (l) prec = l;
    1838           6 :   q = check_unit_disc("thetanullk", q, prec);
    1839             : 
    1840           6 :   if (!odd(k)) { set_avma(av); return gen_0; }
    1841           6 :   qn = gen_1;
    1842           6 :   ps2 = gsqr(q);
    1843           6 :   ps = gneg_i(ps2);
    1844           6 :   y = gen_1;
    1845           6 :   for (n = 3;; n += 2)
    1846         240 :   {
    1847             :     GEN t;
    1848         246 :     qn = gmul(qn,ps);
    1849         246 :     ps = gmul(ps,ps2);
    1850         246 :     t = gmul(qn, powuu(n, k)); y = gadd(y, t);
    1851         246 :     if (gexpo(t) < -prec2nbits(prec)) break;
    1852             :   }
    1853           6 :   p1 = gmul2n(gsqrt(gsqrt(q,prec),prec),1);
    1854           6 :   if (k&2) y = gneg_i(y);
    1855           6 :   return gerepileupto(av, gmul(p1, y));
    1856             : }
    1857             : 
    1858             : /* q2 = q^2 */
    1859             : static GEN
    1860       60765 : vecthetanullk_loop(GEN q2, long k, long prec)
    1861             : {
    1862       60765 :   GEN ps, qn = gen_1, y = const_vec(k, gen_1);
    1863       60765 :   pari_sp av = avma;
    1864       60765 :   const long bit = prec2nbits(prec);
    1865             :   long i, n;
    1866             : 
    1867       60765 :   if (gexpo(q2) < -2*bit) return y;
    1868       60765 :   ps = gneg_i(q2);
    1869       60765 :   for (n = 3;; n += 2)
    1870      307618 :   {
    1871      368383 :     GEN t = NULL/*-Wall*/, P = utoipos(n), N2 = sqru(n);
    1872      368383 :     qn = gmul(qn,ps);
    1873      368383 :     ps = gmul(ps,q2);
    1874     1105149 :     for (i = 1; i <= k; i++)
    1875             :     {
    1876      736766 :       t = gmul(qn, P); gel(y,i) = gadd(gel(y,i), t);
    1877      736766 :       P = mulii(P, N2);
    1878             :     }
    1879      368383 :     if (gexpo(t) < -bit) return y;
    1880      307618 :     if (gc_needed(av,2))
    1881             :     {
    1882           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"vecthetanullk_loop, n = %ld",n);
    1883           0 :       gerepileall(av, 3, &qn, &ps, &y);
    1884             :     }
    1885             :   }
    1886             : }
    1887             : /* [d^i theta/dz^i(q, 0), i = 1, 3, .., 2*k - 1] */
    1888             : GEN
    1889           0 : vecthetanullk(GEN q, long k, long prec)
    1890             : {
    1891           0 :   long i, l = precision(q);
    1892           0 :   pari_sp av = avma;
    1893             :   GEN p1, y;
    1894             : 
    1895           0 :   if (l) prec = l;
    1896           0 :   q = check_unit_disc("vecthetanullk", q, prec);
    1897           0 :   y = vecthetanullk_loop(gsqr(q), k, prec);
    1898           0 :   p1 = gmul2n(gsqrt(gsqrt(q,prec),prec),1);
    1899           0 :   for (i = 2; i <= k; i += 2) gel(y,i) = gneg_i(gel(y,i));
    1900           0 :   return gerepileupto(av, gmul(p1, y));
    1901             : }
    1902             : 
    1903             : /* [d^i theta/dz^i(q, 0), i = 1, 3, .., 2*k - 1], q = exp(2iPi tau) */
    1904             : GEN
    1905           0 : vecthetanullk_tau(GEN tau, long k, long prec)
    1906             : {
    1907           0 :   long i, l = precision(tau);
    1908           0 :   pari_sp av = avma;
    1909             :   GEN q4, y;
    1910             : 
    1911           0 :   if (l) prec = l;
    1912           0 :   if (typ(tau) != t_COMPLEX || gsigne(gel(tau,2)) <= 0)
    1913           0 :     pari_err_DOMAIN("vecthetanullk_tau", "imag(tau)", "<=", gen_0, tau);
    1914           0 :   q4 = expIPiC(gmul2n(tau,-1), prec); /* q^(1/4) */
    1915           0 :   y = vecthetanullk_loop(gpowgs(q4,8), k, prec);
    1916           0 :   for (i = 2; i <= k; i += 2) gel(y,i) = gneg_i(gel(y,i));
    1917           0 :   return gerepileupto(av, gmul(gmul2n(q4,1), y));
    1918             : }
    1919             : 
    1920             : /* Return E_k(tau). Slow if tau is not in standard fundamental domain */
    1921             : GEN
    1922       60996 : cxEk(GEN tau, long k, long prec)
    1923             : {
    1924             :   pari_sp av;
    1925             :   GEN q, y, qn;
    1926       60996 :   long n, b, l = precision(tau);
    1927             : 
    1928       60996 :   if (l) prec = l;
    1929       60996 :   b = prec2nbits(prec);
    1930             :   /* sum n^(k-1) x^n <= x(1 + (k!-1)x) / (1-x)^k (cf Eulerian polynomials)
    1931             :    * S = \sum_{n > 0} n^(k-1) |q^n/(1-q^n)| <= x(1+(k!-1)x) / (1-x)^(k+1),
    1932             :    * where x = |q| = exp(-2Pi Im(tau)) < 1. Neglegt 2/zeta(1-k) * S if
    1933             :    * (2Pi)^k/(k-1)! x < 2^(-b-1) and k! x < 1. Use log2((2Pi)^k/(k-1)!) < 10 */
    1934       60996 :   if (gcmpgs(imag_i(tau), (M_LN2 / (2*M_PI)) * (b+1+10)) > 0)
    1935          15 :     return real_1(prec);
    1936       60981 :   if (k == 2)
    1937             :   { /* -theta^(3)(tau/2) / theta^(1)(tau/2). Assume that Im tau > 0 */
    1938       60765 :     y = vecthetanullk_loop(qq(tau,prec), 2, prec);
    1939       60765 :     return gdiv(gel(y,2), gel(y,1));
    1940             :   }
    1941         216 :   q = cxtoreal(expIPiC(gneg(gmul2n(tau,1)), prec));
    1942         216 :   av = avma; y = gen_0; qn = q;
    1943         216 :   for(n = 1;; n++)
    1944        2246 :   { /* compute y := sum_{n>0} n^(k-1) / (q^n-1) */
    1945        2462 :     GEN t = gdiv(powuu(n,k-1), gsubgs(qn,1));
    1946        2462 :     if (gequal0(t) || gexpo(t) <= - prec2nbits(prec) - 5) break;
    1947        2246 :     y = gadd(y, t);
    1948        2246 :     qn = gmul(q,qn);
    1949        2246 :     if (gc_needed(av,2))
    1950             :     {
    1951           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"elleisnum");
    1952           0 :       gerepileall(av, 2, &y,&qn);
    1953             :     }
    1954             :   }
    1955         216 :   return gadd(gen_1, gmul(y, gdiv(gen_2, szeta(1-k, prec))));
    1956             : }
    1957             : 
    1958             : /********************************************************************/
    1959             : /*      New code for 1-variable theta, does not use AGM             */
    1960             : /********************************************************************/
    1961             : 
    1962             : /* This code provides four new functions: theta (same name as before)
    1963             :    flag = NULL for backward compatibility; otherwise flags
    1964             :    1,2,3,4 old Jacobi notation, flags [0,0], [0,1], [1,0], [1,1],
    1965             :    Riemann notation corresponding to 3, 4, 2, and -1, and flag 0
    1966             :    all four [i,j] with 0 <= i,j <= 1. For instance
    1967             :    vector(4,n,theta(z,tau,[(n-1)\2,(n-1)%2])) and theta(z,tau,0)
    1968             :    should both be identical to riemann_theta([z]~,Mat(tau)) from Jean Kieffer.
    1969             : 
    1970             :    thetanull: theta Nullwerte, same flags, value at z=0 except for flag
    1971             :    1 and [1,1], first derivative at z=0 of [1,1] since value identically zero.
    1972             : 
    1973             :    elljacobi: Jacobi elliptic functions [sn,cn,dn]. Included since trivially
    1974             :    obtained from the four theta functions.
    1975             : 
    1976             :    ellweierstrass: implementation of all Weierstrass data from theta functions.
    1977             :    More precisely, to be written. */
    1978             : 
    1979             : static long
    1980          96 : equali01(GEN x)
    1981             : {
    1982          96 :   if (!signe(x)) return 0;
    1983          72 :   if (!equali1(x)) pari_err_FLAG("theta");
    1984          72 :   return 1;
    1985             : }
    1986             : 
    1987             : static long
    1988         180 : thetaflag(GEN v)
    1989             : {
    1990             :   long v1, v2;
    1991         180 :   if (!v) return 0;
    1992         180 :   switch(typ(v))
    1993             :   {
    1994         132 :     case t_INT:
    1995         132 :       if (signe(v) < 0 || cmpis(v, 4) > 0) pari_err_FLAG("theta");
    1996         132 :       return itou(v);
    1997          48 :     case t_VEC:
    1998          48 :       if (RgV_is_ZV(v) && lg(v) == 3) break;
    1999           0 :     default: pari_err_FLAG("theta");
    2000             :   }
    2001          48 :   v1 = equali01(gel(v,1));
    2002          48 :   v2 = equali01(gel(v,2)); return v1? (v2? -1: 2): (v2? 4: 3);
    2003             : }
    2004             : 
    2005             : /* Automorphy factor for bringing tau towards standard fundamental domain
    2006             :  * (we stop when im(tau) >= 1/2, no need to go all the way to sqrt(3)/2).
    2007             :  * At z = 0 if NULL */
    2008             : static GEN
    2009         240 : autojtau(GEN *pz, GEN *ptau, long *psumr, long *pct, long prec)
    2010             : {
    2011         240 :   GEN S = gen_1, z = *pz, tau = *ptau;
    2012         240 :   long ct = 0, sumr = 0;
    2013         240 :   if (z && gequal0(z)) z = NULL;
    2014         444 :   while (gexpo(imag_i(tau)) < -1)
    2015             :   {
    2016         204 :     GEN r = ground(real_i(tau)), taup;
    2017         204 :     tau = gsub(tau, r); taup = gneg(ginv(tau));
    2018         204 :     S = gdiv(S, gsqrt(mulcxmI(tau), prec));
    2019         204 :     if (z)
    2020             :     {
    2021         144 :       S = gmul(S, expIPiC(gmul(taup, gsqr(z)), prec));
    2022         144 :       z = gneg(gmul(z, taup));
    2023             :     }
    2024         204 :     ct++; tau = taup; sumr = (sumr + Mod8(r)) & 7;
    2025             :   }
    2026         240 :   if (pct) *pct = ct;
    2027         240 :   *psumr = sumr; *pz = z; *ptau = tau; return S;
    2028             : }
    2029             : 
    2030             : /* At 0 if z = NULL. Real(tau) = n is an integer; 4 | n if fl = 1 or 2 */
    2031             : static void
    2032         372 : clearim(GEN *v, GEN z, long fl)
    2033             : {
    2034         372 :   if (!z || gequal0(imag_i(z)) || (fl != 1 && gequal0(real_i(z))))
    2035         282 :     *v = real_i(*v);
    2036         372 : }
    2037             : 
    2038             : static GEN
    2039         234 : clearimall(GEN z, GEN tau, GEN VS)
    2040             : {
    2041             :   GEN n;
    2042         234 :   if (isint(real_i(tau), &n))
    2043             :   {
    2044         138 :     long nmod4 = Mod4(n);
    2045         138 :     clearim(&gel(VS,1), z, 3);
    2046         138 :     clearim(&gel(VS,2), z, 4);
    2047         138 :     if (!nmod4)
    2048             :     {
    2049          48 :       clearim(&gel(VS,3), z, 2);
    2050          48 :       clearim(&gel(VS,4), z, 1);
    2051             :     }
    2052             :   }
    2053         234 :   return VS;
    2054             : }
    2055             : 
    2056             : static GEN
    2057          18 : gadd3(GEN a, GEN b, GEN c) { return gadd(gadd(a, b), c); }
    2058             : static GEN
    2059          36 : gmul3(GEN a, GEN b, GEN c) { return gmul(gmul(a, b), c); }
    2060             : static GEN
    2061           6 : gmul4(GEN a, GEN b, GEN c, GEN d) { return gmul(gmul(a, b), gmul(c,d)); }
    2062             : 
    2063             : /* Implementation of all 4 theta functions */
    2064             : 
    2065             : /* If z = NULL, we are at 0 */
    2066             : static long
    2067         246 : thetaprec(GEN z, GEN tau, long prec)
    2068             : {
    2069         246 :   long l = precision(tau);
    2070         246 :   if (z)
    2071             :   {
    2072         192 :     long n = precision(z);
    2073         192 :     if (n && n < l) l = n;
    2074             :   }
    2075         246 :   return l? l: prec;
    2076             : }
    2077             : 
    2078             : static GEN
    2079          18 : redmodZ(GEN z)
    2080             : {
    2081          18 :   if (typ(z) != t_COMPLEX) return gfrac(z);
    2082           6 :   return mkcomplex(gfrac(gel(z,1)), gel(z,2));
    2083             : }
    2084             : static GEN
    2085         180 : redmod2Z(GEN z)
    2086             : {
    2087         180 :   GEN k = ground(gmul2n(real_i(z), -1));
    2088         180 :   if (signe(k)) z = gsub(z, shifti(k, 1));
    2089         180 :   return z;
    2090             : }
    2091             : 
    2092             : /* If z = NULL assume z = 0 and compute theta[1,1]' instead of theta[1,1] = 0 */
    2093             : static GEN
    2094         234 : thetaall_i(GEN z, GEN tau, long prec)
    2095             : {
    2096             :   GEN zold, tauold, k, u, un, q, q2, qd, qn;
    2097             :   GEN S, S00, S01, S10, S11, u2, ui2, uin;
    2098         234 :   long n, ct, eS, B, sumr, precold = prec;
    2099         234 :   int theta1p = !z;
    2100             : 
    2101         234 :   if (z) z = redmod2Z(z);
    2102         234 :   tau = upper_to_cx(tau, &prec);
    2103         234 :   prec = thetaprec(z, tau, prec);
    2104         234 :   z = zold = z? gtofp(z, prec): NULL;
    2105         234 :   tau = tauold = gtofp(tau, prec);
    2106         234 :   S = autojtau(&z, &tau, &sumr, &ct, prec);
    2107         234 :   k = gen_0; S00 = S01 = gen_1; S10 = S11 = gen_0;
    2108         234 :   if (z)
    2109             :   {
    2110         156 :     GEN y = imag_i(z);
    2111         156 :     if (!gequal0(y)) k = roundr(divrr(y, gneg(imag_i(tau))));
    2112         156 :     if (signe(k))
    2113             :     {
    2114          12 :       GEN Sz = expIPiC(gadd(gmul(sqri(k), tau), gmul(shifti(k,1), z)), prec);
    2115          12 :       S = gmul(S, Sz);
    2116          12 :       z = gadd(z, gmul(tau, k));
    2117             :     }
    2118             :   }
    2119         234 :   if ((eS = gexpo(S)) > 0)
    2120             :   {
    2121          42 :     prec = nbits2prec(eS + prec2nbits(prec));
    2122          42 :     if (z) z = gprec_w(z, prec);
    2123          42 :     tau = gprec_w(tau, prec);
    2124             :   }
    2125         234 :   q = expIPiC(gmul2n(tau,-2), prec); q2 = gsqr(q); qn = gen_1;
    2126         234 :   if (!z) u = u2 = ui2 = un = uin = NULL; /* constant, equal to 1 */
    2127             :   else
    2128             :   {
    2129         156 :     u = expIPiC(z, prec); u2 = gsqr(u); ui2 = ginv(u2);
    2130         156 :     un = uin = gen_1;
    2131             :   }
    2132         234 :   qd = q; B = prec2nbits(prec);
    2133         234 :   for (n = 1;; n++)
    2134        1230 :   { /* qd = q^(4n-3), qn = q^(4(n-1)^2), un = u^(2n-2), uin = 1/un */
    2135        1464 :     long e = 0, eqn, prec2;
    2136             :     GEN tmp;
    2137        1464 :     if (u) uin = gmul(uin, ui2);
    2138        1464 :     qn = gmul(qn, qd); /* q^((2n-1)^2) */
    2139        1464 :     tmp = u? gmul(qn, gadd(un, uin)): gmul2n(qn, 1);
    2140        1464 :     S10 = gadd(S10, tmp);
    2141        1464 :     if (z)
    2142             :     {
    2143         942 :       tmp = gmul(qn, gsub(un, uin));
    2144         942 :       S11 = odd(n)? gsub(S11, tmp): gadd(S11, tmp);
    2145         942 :       e = maxss(0, gexpo(un)); un = gmul(un, u2);
    2146         942 :       e = maxss(e, gexpo(un));
    2147             :     }
    2148         522 :     else if (theta1p) /* theta'[1,1] at 0 */
    2149             :     {
    2150         360 :       tmp = gmulsg(2*n-1, tmp);
    2151         360 :       S11 = odd(n)? gsub(S11, tmp): gadd(S11, tmp);
    2152             :     }
    2153        1464 :     qd = gmul(qd, q2); qn = gmul(qn, qd); /* q^(4n^2) */
    2154        1464 :     tmp = u? gmul(qn, gadd(un, uin)): gmul2n(qn, 1);
    2155        1464 :     S00 = gadd(S00, tmp);
    2156        1464 :     S01 = odd(n)? gsub(S01, tmp): gadd(S01, tmp);
    2157        1464 :     eqn = gexpo(qn) + e; if (eqn < -B) break;
    2158        1230 :     qd = gmul(qd, q2);
    2159        1230 :     prec2 = minss(prec, nbits2prec(eqn + B + 64));
    2160        1230 :     qn = gprec_w(qn, prec2); qd = gprec_w(qd, prec2);
    2161        1230 :     if (u) { un = gprec_w(un, prec2); uin = gprec_w(uin, prec2); }
    2162             :   }
    2163         234 :   if (u)
    2164             :   {
    2165         156 :     S10 = gmul(u, S10);
    2166         156 :     S11 = gmul(u, S11);
    2167             :   }
    2168             :   /* automorphic factor
    2169             :    *   theta[1,1]: I^ct
    2170             :    *   theta[1,0]: exp(-I*Pi/4*sumr)
    2171             :    *   theta[0,1]: (-1)^k
    2172             :    *   theta[1,1]: (-1)^k exp(-I*Pi/4*sumr) */
    2173         234 :   S11 = z? mulcxpowIs(S11, ct + 3): gmul(mppi(prec), S11);
    2174         234 :   if (ct&1L) swap(S10, S01);
    2175         234 :   if (sumr & 7)
    2176             :   {
    2177          90 :     GEN zet = e12(sumr * 3, prec); /* exp(I Pi sumr / 4) */
    2178          90 :     if (odd(sumr)) swap(S01, S00);
    2179          90 :     S10 = gmul(S10, zet);
    2180          90 :     S11 = gmul(S11, zet);
    2181             :   }
    2182         234 :   if (theta1p) S11 = gmul(gsqr(S), S11);
    2183         180 :   else if (mpodd(k)) { S01 = gneg(S01); S11 = gneg(S11); }
    2184         234 :   if (precold < prec)
    2185             :   {
    2186          54 :     prec = precold;
    2187          54 :     S00 = gprec_w(S00, prec); S01 = gprec_w(S01, prec);
    2188          54 :     S10 = gprec_w(S10, prec); S11 = gprec_w(S11, prec);
    2189             :   }
    2190         234 :   return clearimall(zold, tauold, gmul(S, mkvec4(S00, S01, S10, S11)));
    2191             : }
    2192             : 
    2193             : static GEN
    2194          54 : thetanull_i(GEN tau, long prec) { return thetaall_i(NULL, tau, prec); }
    2195             : 
    2196             : GEN
    2197         150 : theta(GEN z, GEN tau, GEN flag, long prec)
    2198             : {
    2199         150 :   pari_sp av = avma;
    2200             :   GEN T;
    2201         150 :   if (!flag)
    2202             :   { /* backward compatibility: sine theta */
    2203          12 :     GEN pi = mppi(prec), q = z; z = tau; /* input (q = exp(i pi tau), Pi*z) */
    2204          12 :     prec = thetaprec(z, tau, prec);
    2205          12 :     q = check_unit_disc("theta", q, prec);
    2206          12 :     z = gdiv(gtofp(z, prec), pi);
    2207          12 :     tau = gdiv(mulcxmI(glog(q, prec)), pi);
    2208          12 :     flag = gen_1;
    2209             :   }
    2210         150 :   T = thetaall_i(z, tau, prec);
    2211         150 :   switch (thetaflag(flag))
    2212             :   {
    2213          24 :     case -1: T = gel(T,4); break;
    2214          18 :     case 0: break;
    2215          72 :     case 1: T = gneg(gel(T,4)); break;
    2216          12 :     case 2: T = gel(T,3); break;
    2217          12 :     case 3: T = gel(T,1); break;
    2218          12 :     case 4: T = gel(T,2); break;
    2219           0 :     default: pari_err_FLAG("theta");
    2220             :   }
    2221         150 :   return gerepilecopy(av, T);
    2222             : }
    2223             : 
    2224             : /* Same as 2*Pi*eta(tau,1)^3 = - thetaall_i(NULL, tau)[4], faster than both. */
    2225             : static GEN
    2226           6 : thetanull11(GEN tau, long prec)
    2227             : {
    2228           6 :   GEN z = NULL, tauold, q, q8, qd, qn, S, S11;
    2229           6 :   long n, eS, B, sumr, precold = prec;
    2230             : 
    2231           6 :   tau = upper_to_cx(tau, &prec);
    2232           6 :   tau = tauold = gtofp(tau, prec);
    2233           6 :   S = autojtau(&z, &tau, &sumr, NULL, prec);
    2234           6 :   S11 = gen_1; ;
    2235           6 :   if ((eS = gexpo(S)) > 0)
    2236             :   {
    2237           0 :     prec += nbits2extraprec(eS);
    2238           0 :     tau = gprec_w(tau, prec);
    2239             :   }
    2240           6 :   q8 = expIPiC(gmul2n(tau,-2), prec); q = gpowgs(q8, 8);
    2241           6 :   qn = gen_1; qd = q; B = prec2nbits(prec);
    2242           6 :   for (n = 1;; n++)
    2243          36 :   { /* qd = q^n, qn = q^((n^2-n)/2) */
    2244             :     long eqn, prec2;
    2245             :     GEN tmp;
    2246          42 :     qn = gmul(qn, qd); tmp = gmulsg(2*n+1, qn); eqn = gexpo(tmp);
    2247          42 :     S11 = odd(n)? gsub(S11, tmp): gadd(S11, tmp);
    2248          42 :     if (eqn < -B) break;
    2249          36 :     qd = gmul(qd, q);
    2250          36 :     prec2 = minss(prec, nbits2prec(eqn + B + 32));
    2251          36 :     qn = gprec_w(qn, prec2); qd = gprec_w(qd, prec2);
    2252             :   }
    2253           6 :   if (precold < prec) prec = precold;
    2254           6 :   S11 = gmul3(S11, q8, e12(3*sumr, prec));
    2255           6 :   S11 = gmul3(Pi2n(1, prec), gpowgs(S, 3), S11);
    2256           6 :   if (isint(real_i(tauold), &q) && !Mod4(q)) clearim(&S11, z, 1);
    2257           6 :   return S11;
    2258             : }
    2259             : 
    2260             : GEN
    2261          30 : thetanull(GEN tau, GEN flag, long prec)
    2262             : {
    2263          30 :   pari_sp av = avma;
    2264          30 :   long fl = thetaflag(flag);
    2265             :   GEN T0;
    2266          30 :   if (fl == 1) T0 = thetanull11(tau, prec);
    2267          30 :   else if (fl == -1) T0 = gneg(thetanull11(tau, prec));
    2268             :   else
    2269             :   {
    2270          24 :     T0 = thetanull_i(tau, prec);
    2271          24 :     switch (fl)
    2272             :     {
    2273           6 :       case 0: break;
    2274           6 :       case 2: T0 = gel(T0,3); break;
    2275           6 :       case 3: T0 = gel(T0,1); break;
    2276           6 :       case 4: T0 = gel(T0,2); break;
    2277           0 :       default: pari_err_FLAG("thetanull");
    2278             :     }
    2279             :   }
    2280          30 :   return gerepilecopy(av, T0);
    2281             : }
    2282             : 
    2283             : /********************************************************************/
    2284             : /**                     Jacobi sn, cn, dn                          **/
    2285             : /********************************************************************/
    2286             : 
    2287             : /* Basic Jacobi elliptic functions in terms of thetas */
    2288             : GEN
    2289          12 : elljacobi(GEN z, GEN k, long prec)
    2290             : {
    2291          12 :   pari_sp av = avma;
    2292          12 :   GEN K = ellK(k, prec), Kp = ellK(gsqrt(gsubsg(1, gsqr(k)), prec), prec);
    2293          12 :   GEN zet = gdiv(gmul2n(z, -1), K), tau = mulcxI(gdiv(Kp, K));
    2294          12 :   GEN T = thetaall_i(zet, tau, prec), T0 = thetanull_i(tau, prec);
    2295          12 :   GEN t1 = gneg(gel(T,4)), t2 = gel(T,3), t3 = gel(T,1), t4 = gel(T,2);
    2296          12 :   GEN z2 = gel(T0,3), z3 = gel(T0,1), z4 = gel(T0,2), z2t4 = gmul(z2, t4);
    2297             :   GEN SN, CN, DN;
    2298          12 :   SN = gdiv(gmul(z3, t1), z2t4);
    2299          12 :   CN = gdiv(gmul(z4, t2), z2t4);
    2300          12 :   DN = gdiv(gmul(z4, t3), gmul(z3, t4));
    2301          12 :   return gerepilecopy(av, mkvec3(SN, CN, DN));
    2302             : }
    2303             : 
    2304             : /********************************************************************/
    2305             : /**           Weierstrass elliptic data in terms of thetas         **/
    2306             : /********************************************************************/
    2307             : 
    2308             : static GEN
    2309          18 : mfE2eval(GEN tau, long prec)
    2310             : {
    2311          18 :   GEN M, c, taured = cxredsl2(tau, &M), R = cxEk(taured, 2, prec);
    2312          18 :   c = gcoeff(M,2,1);
    2313          18 :   if (signe(c))
    2314             :   {
    2315          18 :     GEN d = gcoeff(M,2,2), J = gadd(d, gmul(c, tau));
    2316          18 :     R = gadd(R, gdiv(mulcxI(gmul(mului(6,c), J)), mppi(prec)));
    2317          18 :     R = gdiv(R, gsqr(J));
    2318             :   }
    2319          18 :   return R;
    2320             : }
    2321             : 
    2322             : static GEN
    2323          18 : weta1eta2(GEN tau, GEN e, long prec)
    2324          18 : { return mkvec2(gneg(e), gsub(gmul(tau, e), PiI2n(1, prec))); }
    2325             : 
    2326             : static GEN
    2327          18 : wg2g3(GEN T0, GEN a, GEN *pe)
    2328             : {
    2329          18 :   GEN z2 = gel(T0,3), z3 = gel(T0,1), z4 = gel(T0,2), g2, g3, e1, e2, e3;
    2330          18 :   z2 = gmul(a, gpowgs(z2, 4));
    2331          18 :   z3 = gmul(a, gpowgs(z3, 4));
    2332          18 :   z4 = gmul(a, gpowgs(z4, 4));
    2333          18 :   e1 = gadd(z3, z4); e2 = gsub(z2, z4); e3 = gneg(gadd(z2, z3));
    2334          18 :   g2 = gmulgs(gadd3(gsqr(z2), gsqr(z3), gsqr(z4)), 6);
    2335          18 :   g3 = gmul2n(gmul3(e1, e2, e3), 2);
    2336          18 :   *pe = mkvec3(e1, e2, e3); return mkvec2(g2, g3);
    2337             : }
    2338             : 
    2339             : static GEN
    2340           6 : wsigma(GEN z, GEN T, GEN T0, GEN e, long prec)
    2341             : {
    2342           6 :   GEN E = gexp(gmul(gsqr(z), gmul2n(e,-1)), prec);
    2343           6 :   return gdiv(gmul(E, gel(T,4)), gel(T0,4));
    2344             : }
    2345             : 
    2346             : static GEN
    2347          12 : _thetaz(void *E, GEN z, long prec)
    2348          12 : { GEN T = thetaall_i(z, (GEN)E, prec); return gel(T,4); }
    2349             : static GEN
    2350           6 : wzeta(GEN z, GEN tau, GEN T, GEN e, long prec)
    2351             : {
    2352           6 :   GEN TP = derivnum((void*)tau, _thetaz, z, prec);
    2353           6 :   return gadd(gmul(z, e), gdiv(TP, gel(T,4)));
    2354             : }
    2355             : 
    2356             : static GEN
    2357           6 : wp(GEN T, GEN T0, GEN a)
    2358             : {
    2359           6 :   GEN t4 = gel(T,2), z2 = gel(T0,3), z3 = gel(T0,1);
    2360           6 :   return gmul(a, gsub(gmulgs(gsqr(gdiv(gmul3(z2, z3, t4), gel(T,4))), 3),
    2361             :                       gadd(gpowgs(z2,4), gpowgs(z3,4))));
    2362             : }
    2363             : 
    2364             : static GEN
    2365           6 : wpprime(GEN T, GEN T0, long prec)
    2366             : {
    2367           6 :   GEN t2 = gel(T,3), t3 = gel(T,1), t4 = gel(T,2);
    2368           6 :   GEN c = gmul(Pi2n(1, prec), gsqr(gel(T0,4)));
    2369           6 :   return gdiv(gmul4(c, t2, t3, t4), gpowgs(gel(T,4), 3));
    2370             : }
    2371             : 
    2372             : GEN
    2373          18 : ellweierstrass(GEN z, GEN tau, long prec)
    2374             : {
    2375          18 :   pari_sp av = avma;
    2376          18 :   long prec2 = prec + EXTRAPREC64;
    2377          18 :   GEN T0 = thetanull_i(tau, prec), a = divru(sqrr(mppi(prec2)), 3);
    2378          18 :   GEN e = gmul(a, mfE2eval(tau, prec2));
    2379          18 :   GEN e1e2e3, g2g3 = wg2g3(T0, a, &e1e2e3);
    2380          18 :   GEN R = mkvec4(mkvec2(gen_m1, tau), g2g3, e1e2e3, weta1eta2(tau, e, prec));
    2381          18 :   z = redmodZ(z);
    2382          18 :   if (!gequal0(z))
    2383             :   {
    2384           6 :     GEN T = thetaall_i(z, tau, prec);
    2385           6 :     R = shallowconcat(R, mkvec4(wsigma(z, T, T0, e, prec),
    2386             :                          wzeta(z, tau, T, e, prec),
    2387             :                          wp(T, T0, a), wpprime(T, T0, prec)));
    2388             :   }
    2389          18 :   return gerepilecopy(av, R);
    2390             : }

Generated by: LCOV version 1.16