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