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