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