1 /* $Id: trans3.c 10287 2008-06-09 22:03:21Z kb $
2
3 Copyright (C) 2000 The PARI group.
4
5 This file is part of the PARI/GP package.
6
7 PARI/GP is free software; you can redistribute it and/or modify it under the
8 terms of the GNU General Public License as published by the Free Software
9 Foundation. It is distributed in the hope that it will be useful, but WITHOUT
10 ANY WARRANTY WHATSOEVER.
11
12 Check the License for details. You should have received a copy of it, along
13 with the package; see the file 'COPYING'. If not, write to the Free Software
14 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
15
16 /********************************************************************/
17 /** **/
18 /** TRANSCENDENTAL FONCTIONS **/
19 /** (part 3) **/
20 /** **/
21 /********************************************************************/
22 #include "pari.h"
23 #include "paripriv.h"
24
25 /***********************************************************************/
26 /** **/
27 /** BESSEL FUNCTIONS **/
28 /** **/
29 /***********************************************************************/
30
31 /* n! sum_{k=0}^m Z^k / (k!*(k+n)!), with Z := (-1)^flag*z^2/4 */
32 static GEN
_jbessel(GEN n,GEN z,long flag,long m)33 _jbessel(GEN n, GEN z, long flag, long m)
34 {
35 long k;
36 pari_sp av, lim;
37 GEN Z,s;
38
39 Z = gmul2n(gsqr(z),-2); if (flag & 1) Z = gneg(Z);
40 if (typ(z) == t_SER)
41 {
42 long v = valp(z);
43 k = lg(Z)-2 - v;
44 if (v < 0) pari_err(negexper,"jbessel");
45 if (v == 0) pari_err(impl,"jbessel around a!=0");
46 if (k <= 0) return gadd(gen_1, zeroser(varn(z), 2*v));
47 Z = gprec(Z, k);
48 }
49 s = gen_1;
50 av = avma; lim = stack_lim(av,1);
51 for (k=m; k>=1; k--)
52 {
53 s = gaddsg(1, gdiv(gdivgs(gmul(Z,s),k),gaddsg(k,n)));
54 if (low_stack(lim, stack_lim(av,1)))
55 {
56 if (DEBUGMEM>1) pari_warn(warnmem,"jbessel");
57 s = gerepilecopy(av, s);
58 }
59 }
60 return s;
61 }
62
63 /* return L * approximate solution to x log x = B */
64 static long
bessel_get_lim(double B,double L)65 bessel_get_lim(double B, double L)
66 {
67 long lim;
68 double x = 1 + B; /* 3 iterations are enough except in pathological cases */
69 x = (x + B)/(log(x)+1);
70 x = (x + B)/(log(x)+1);
71 x = (x + B)/(log(x)+1); x = L*x;
72 lim = (long)x; if (lim < 2) lim = 2;
73 return lim;
74 }
75
76 static GEN
jbesselintern(GEN n,GEN z,long flag,long prec)77 jbesselintern(GEN n, GEN z, long flag, long prec)
78 {
79 long i, lz, lim, k, ki, precnew;
80 pari_sp av = avma;
81 double B, L;
82 GEN p1, p2, y;
83
84 switch(typ(z))
85 {
86 case t_INT: case t_FRAC: case t_QUAD:
87 case t_REAL: case t_COMPLEX:
88 i = precision(z); if (i) prec = i;
89 p2 = gdiv(gpow(gmul2n(z,-1),n,prec), ggamma(gaddgs(n,1),prec));
90 if (gcmp0(z)) return gerepilecopy(av, p2);
91
92 L = 1.3591409 * gtodouble(gabs(z,prec));
93 precnew = prec;
94 if (L >= 1.0) precnew += 1 + (long)(L/(1.3591409*LOG2*BITS_IN_LONG));
95 if (issmall(n,&ki))
96 {
97 k = labs(ki);
98 n = stoi(k);
99 } else {
100 i = precision(n);
101 if (i && i < precnew) n = gtofp(n,precnew);
102 }
103 z = gtofp(z,precnew);
104 B = bit_accuracy_mul(prec, LOG2/2) / L;
105 lim = bessel_get_lim(B, L);
106 p1 = gprec_wtrunc(_jbessel(n,z,flag,lim), prec);
107 return gerepileupto(av, gmul(p2,p1));
108
109 case t_VEC: case t_COL: case t_MAT:
110 lz=lg(z); y=cgetg(lz,typ(z));
111 for (i=1; i<lz; i++)
112 gel(y,i) = jbesselintern(n,gel(z,i),flag,prec);
113 return y;
114
115 case t_POLMOD:
116 y = cleanroots(gel(z,1), prec); lz = lg(y);
117 for (i=1; i<lz; i++) {
118 GEN t = poleval(gel(z,2), gel(y,i));
119 gel(y,i) = jbesselintern(n,t,flag,prec);
120 }
121 return gerepilecopy(av,y);
122
123 case t_PADIC: pari_err(impl,"p-adic jbessel function");
124 default:
125 if (!(y = toser_i(z))) break;
126 if (issmall(n,&ki)) n = stoi(labs(ki));
127 return gerepilecopy(av, _jbessel(n,y,flag,lg(y)-2));
128 }
129 pari_err(typeer,"jbessel");
130 return NULL; /* not reached */
131 }
132 GEN
jbessel(GEN n,GEN z,long prec)133 jbessel(GEN n, GEN z, long prec) { return jbesselintern(n,z,1,prec); }
134 GEN
ibessel(GEN n,GEN z,long prec)135 ibessel(GEN n, GEN z, long prec) { return jbesselintern(n,z,0,prec); }
136
137 static GEN
_jbesselh(long k,GEN z,long prec)138 _jbesselh(long k, GEN z, long prec)
139 {
140 GEN s,c,p0,p1,p2, zinv = ginv(z);
141 long i;
142
143 gsincos(z,&s,&c,prec);
144 p1 = gmul(zinv,s);
145 if (k)
146 {
147 p0 = p1; p1 = gmul(zinv,gsub(p0,c));
148 for (i=2; i<=k; i++)
149 {
150 p2 = gsub(gmul(gmulsg(2*i-1,zinv),p1), p0);
151 p0 = p1; p1 = p2;
152 }
153 }
154 return p1;
155 }
156
157 GEN
jbesselh(GEN n,GEN z,long prec)158 jbesselh(GEN n, GEN z, long prec)
159 {
160 long gz, k, l, linit, i, lz;
161 pari_sp av;
162 GEN res, y, p1;
163
164 if (typ(n)!=t_INT) pari_err(talker,"not an integer index in jbesselh");
165 k = itos(n);
166 if (k < 0) return jbessel(gadd(ghalf,n), z, prec);
167
168 switch(typ(z))
169 {
170 case t_INT: case t_FRAC: case t_QUAD:
171 case t_REAL: case t_COMPLEX:
172 if (gcmp0(z))
173 {
174 av = avma;
175 p1 = gmul(gsqrt(gdiv(z,mppi(prec)),prec),gpowgs(z,k));
176 p1 = gdiv(p1, seq_umul(k+1, 2*k+1)); /* x k! / (2k+1)! */
177 return gerepileupto(av, gmul2n(p1,2*k));
178 }
179 gz = gexpo(z);
180 linit = precision(z); if (!linit) linit = prec;
181 res = cgetc(linit);
182 av = avma;
183 if (gz>=0) l = linit;
184 else l = linit - 1 + ((-2*k*gz)>>TWOPOTBITS_IN_LONG);
185 if (l>prec) prec = l;
186 prec += (-gz)>>TWOPOTBITS_IN_LONG;
187 if (prec < 3) prec = 3;
188 z = gadd(z, real_0(prec));
189 if (typ(z) == t_COMPLEX) gel(z,2) = gadd(gel(z,2), real_0(prec));
190 p1 = gmul(_jbesselh(k,z,prec), gsqrt(gdiv(z,Pi2n(-1,prec)),prec));
191 avma = av;
192 if (typ(p1) == t_COMPLEX)
193 {
194 affr_fixlg(gel(p1,1), gel(res,1));
195 affr_fixlg(gel(p1,2), gel(res,2));
196 }
197 else
198 {
199 res = cgetr(linit);
200 affr_fixlg(p1, res);
201 }
202 return res;
203
204 case t_VEC: case t_COL: case t_MAT:
205 lz=lg(z); y=cgetg(lz,typ(z));
206 for (i=1; i<lz; i++) gel(y,i) = jbesselh(n,gel(z,i),prec);
207 return y;
208
209 case t_POLMOD:
210 av = avma;
211 y = cleanroots(gel(z,1), prec); lz = lg(y);
212 for (i=1; i<lz; i++) {
213 GEN t = poleval(gel(z,2), gel(y,i));
214 gel(y,i) = jbesselh(n,t,prec);
215 }
216 return gerepilecopy(av, y);
217
218 case t_PADIC: pari_err(impl,"p-adic jbesselh function");
219 default:
220 av = avma;
221 if (!(y = toser_i(z))) break;
222 if (gcmp0(y)) return gpowgs(y,k);
223 if (valp(y) < 0) pari_err(negexper,"jbesselh");
224 y = gprec(y, lg(y)-2 + (2*k+1)*valp(y));
225 p1 = gdiv(_jbesselh(k,y,prec),gpowgs(y,k));
226 for (i=1; i<=k; i++) p1 = gmulgs(p1,2*i+1);
227 return gerepilecopy(av,p1);
228 }
229 pari_err(typeer,"jbesselh");
230 return NULL; /* not reached */
231 }
232
233 GEN
kbessel2(GEN nu,GEN x,long prec)234 kbessel2(GEN nu, GEN x, long prec)
235 {
236 pari_sp av = avma;
237 GEN p1, x2, a;
238
239 if (typ(x)==t_REAL) prec = lg(x);
240 x2 = gshift(x,1);
241 a = gcmp0(imag_i(nu))? cgetr(prec): cgetc(prec);
242 gaddz(gen_1,gshift(nu,1), a);
243 p1 = hyperu(gshift(a,-1),a,x2,prec);
244 p1 = gmul(gmul(p1,gpow(x2,nu,prec)), sqrtr(mppi(prec)));
245 return gerepileupto(av, gdiv(p1,gexp(x,prec)));
246 }
247
248 GEN
kbessel(GEN nu,GEN gx,long prec)249 kbessel(GEN nu, GEN gx, long prec)
250 {
251 GEN x,y,yfin,p1,p2,zf,zz,s,t,q,r,u,v,e,f,c,d,ak,pitemp,nu2,w;
252 long l, lnew, k, k2, l1, n2, n, ex, rab=0;
253 pari_sp av, av1;
254
255 if (typ(nu)==t_COMPLEX) return kbessel2(nu,gx,prec);
256 l = (typ(gx)==t_REAL)? lg(gx): prec;
257 ex = gexpo(gx);
258 if (ex < 0)
259 {
260 rab = 1 + ((-ex)>>TWOPOTBITS_IN_LONG);
261 lnew = l + rab; prec += rab;
262 }
263 else lnew = l;
264 yfin=cgetr(l); l1=lnew+1;
265 av=avma; x = gtofp(gx, lnew); y=cgetr(lnew);
266 u=cgetr(l1); v=cgetr(l1); c=cgetr(l1); d=cgetr(l1);
267 e=cgetr(l1); f=cgetr(l1);
268 nu2=gmulgs(gsqr(nu),-4);
269 n = (long) (bit_accuracy_mul(l,LOG2) + PI*sqrt(gtodouble(gnorm(nu)))) / 2;
270 n2=(n<<1); pitemp=mppi(l1);
271 av1=avma;
272 if (cmprs(x, n) < 0)
273 {
274 zf = gsqrt(gdivgs(pitemp,n2),prec);
275 zz = ginv(stor(n2<<2, prec));
276 s=gen_1; t=gen_0;
277 for (k=n2,k2=2*n2-1; k > 0; k--,k2-=2)
278 {
279 if (k2 & ~(MAXHALFULONG>>1))
280 p1 = gadd(mulss(k2,k2),nu2);
281 else
282 p1 = gaddsg(k2*k2,nu2);
283 ak = gdivgs(gmul(p1,zz),-k);
284 s = gaddsg(1, gmul(ak,s));
285 t = gaddsg(k2,gmul(ak,t));
286 }
287 gmulz(s,zf,u); t=gmul2n(t,-1);
288 gdivgsz(gadd(gmul(t,zf),gmul(u,nu)),-n2,v); avma=av1;
289 q = stor(n2, l1);
290 r=gmul2n(x,1); av1=avma;
291 for(;;)
292 {
293 p1=divsr(5,q); if (expo(p1) >= -1) p1=ghalf;
294 p2=subsr(1,divrr(r,q));
295 if (gcmp(p1,p2)>0) p1=p2;
296 gnegz(p1,c); gaffsg(1,d); affrr(u,e); affrr(v,f);
297 for (k=1; ; k++)
298 {
299 w=gadd(gmul(gsubsg(k,ghalf),u), gmul(subrs(q,k),v));
300 w=gadd(w, gmul(nu,gsub(u,gmul2n(v,1))));
301 gdivgsz(gmul(q,v),k,u);
302 gdivgsz(w,k,v);
303 gmulz(d,c,d);
304 gaddz(e,gmul(d,u),e); p1=gmul(d,v);
305 gaddz(f,p1,f);
306 if (gcmp0(p1) || gexpo(p1) - gexpo(f) <= 1-bit_accuracy(precision(p1)))
307 break;
308 avma=av1;
309 }
310 p1=u; u=e; e=p1;
311 p1=v; v=f; f=p1;
312 gmulz(q,gaddsg(1,c),q);
313 if (expo(subrr(q,r)) - expo(r) <= 1-bit_accuracy(lnew)) break;
314 avma=av1;
315 }
316 gmulz(u,gpow(divrs(x,n),nu,prec),y);
317 }
318 else
319 {
320 p2=gmul2n(x,1);
321 zf=gsqrt(gdiv(pitemp,p2),prec);
322 zz=ginv(gmul2n(p2,2)); s=gen_1;
323 for (k=n2,k2=2*n2-1; k > 0; k--,k2-=2)
324 {
325 if (k2 & ~(MAXHALFULONG>>1))
326 p1=gadd(mulss(k2,k2),nu2);
327 else
328 p1=gaddsg(k2*k2,nu2);
329 ak = gdivgs(gmul(p1,zz),k);
330 s = gsub(gen_1,gmul(ak,s));
331 }
332 gmulz(s,zf,y);
333 }
334 gdivz(y,mpexp(x),yfin);
335 avma=av; return yfin;
336 }
337
338 /* sum_{k=0}^m Z^k (H(k)+H(k+n)) / (k! (k+n)!)
339 * + sum_{k=0}^{n-1} (-Z)^(k-n) (n-k-1)!/k! with Z := (-1)^flag*z^2/4.
340 * Warning: contrary to _jbessel, no n! in front.
341 * When flag > 1, compute exactly the H(k) and factorials (slow) */
342 static GEN
_kbessel(long n,GEN z,long flag,long m,long prec)343 _kbessel(long n, GEN z, long flag, long m, long prec)
344 {
345 long k, limit;
346 pari_sp av;
347 GEN Z, p1, p2, s, H;
348
349 Z = gmul2n(gsqr(z),-2); if (flag & 1) Z = gneg(Z);
350 if (typ(z) == t_SER)
351 {
352 long v = valp(z);
353 k = lg(Z)-2 - v;
354 if (v < 0) pari_err(negexper,"kbessel");
355 if (v == 0) pari_err(impl,"kbessel around a!=0");
356 if (k <= 0) return gadd(gen_1, zeroser(varn(z), 2*v));
357 Z = gprec(Z, k);
358 }
359 H = cgetg(m+n+2,t_VEC); gel(H,1) = gen_0;
360 if (flag <= 1)
361 {
362 gel(H,2) = s = real_1(prec);
363 for (k=2; k<=m+n; k++) gel(H,k+1) = s = divrs(addsr(1,mulsr(k,s)),k);
364 }
365 else
366 {
367 gel(H,2) = s = gen_1;
368 for (k=2; k<=m+n; k++) gel(H,k+1) = s = gdivgs(gaddsg(1,gmulsg(k,s)),k);
369 }
370 s = gadd(gel(H,m+1), gel(H,m+n+1));
371 av = avma; limit = stack_lim(av,1);
372 for (k=m; k>0; k--)
373 {
374 s = gadd(gadd(gel(H,k),gel(H,k+n)),gdiv(gmul(Z,s),mulss(k,k+n)));
375 if (low_stack(limit,stack_lim(av,1)))
376 {
377 if (DEBUGMEM>1) pari_warn(warnmem,"kbessel");
378 s = gerepilecopy(av, s);
379 }
380 }
381 p1 = (flag <= 1) ? mpfactr(n,prec) : mpfact(n);
382 s = gdiv(s,p1);
383 if (n)
384 {
385 Z = gneg(ginv(Z));
386 p2 = gmulsg(n, gdiv(Z,p1));
387 s = gadd(s,p2);
388 for (k=n-1; k>0; k--)
389 {
390 p2 = gmul(p2, gmul(mulss(k,n-k),Z));
391 s = gadd(s,p2);
392 }
393 }
394 return s;
395 }
396
397 static GEN
kbesselintern(GEN n,GEN z,long flag,long prec)398 kbesselintern(GEN n, GEN z, long flag, long prec)
399 {
400 long i, k, ki, lz, lim, precnew, fl, fl2, ex;
401 pari_sp av = avma;
402 GEN p1, p2, y, p3, pp, pm, s, c;
403 double B, L;
404
405 fl = (flag & 1) == 0;
406 switch(typ(z))
407 {
408 case t_INT: case t_FRAC: case t_QUAD:
409 case t_REAL: case t_COMPLEX:
410 if (gcmp0(z)) pari_err(talker,"zero argument in a k/n bessel function");
411 i = precision(z); if (i) prec = i;
412 i = precision(n); if (i && prec > i) prec = i;
413 ex = gexpo(z);
414 /* experimental */
415 if (!flag && ex > bit_accuracy(prec)/16 + gexpo(n))
416 return kbessel(n,z,prec);
417 L = 1.3591409 * gtodouble(gabs(z,prec));
418 precnew = prec;
419 if (L >= 1.3591409) {
420 long rab = (long)(L/(1.3591409*LOG2*BITS_IN_LONG));
421 if (fl) rab *= 2;
422 precnew += 1 + rab;
423 }
424 z = gtofp(z, precnew);
425 if (issmall(n,&ki))
426 {
427 GEN z2 = gmul2n(z, -1);
428 k = labs(ki);
429 B = bit_accuracy_mul(prec,LOG2/2) / L;
430 if (fl) B += 0.367879;
431 lim = bessel_get_lim(B, L);
432 p1 = gmul(gpowgs(z2,k), _kbessel(k,z,flag,lim,precnew));
433 p2 = gadd(mpeuler(precnew), glog(z2,precnew));
434 p3 = jbesselintern(stoi(k),z,flag,precnew);
435 p2 = gsub(gmul2n(p1,-1),gmul(p2,p3));
436 p2 = gprec_wtrunc(p2, prec);
437 if (fl) {
438 if (k & 1) p2 = gneg(p2);
439 }
440 else
441 {
442 p2 = gdiv(p2, Pi2n(-1,prec));
443 if (ki >= 0 || (k&1)==0) p2 = gneg(p2);
444 }
445 return gerepilecopy(av, p2);
446 }
447
448 n = gtofp(n, precnew);
449 gsincos(gmul(n,mppi(precnew)), &s,&c,precnew);
450 ex = gexpo(s);
451 if (ex < 0)
452 {
453 long rab = (-ex) >> TWOPOTBITS_IN_LONG;
454 if (fl) rab *= 2;
455 precnew += 1 + rab;
456 }
457 if (i && i < precnew) {
458 n = gtofp(n,precnew);
459 z = gtofp(z,precnew);
460 gsincos(gmul(n,mppi(precnew)), &s,&c,precnew);
461 }
462
463 pp = jbesselintern(n, z,flag,precnew);
464 pm = jbesselintern(gneg(n),z,flag,precnew);
465 if (fl)
466 p1 = gmul(gsub(pm,pp), Pi2n(-1,precnew));
467 else
468 p1 = gsub(gmul(c,pp),pm);
469 p1 = gdiv(p1, s);
470 return gerepilecopy(av, gprec_wtrunc(p1,prec));
471
472 case t_VEC: case t_COL: case t_MAT:
473 lz=lg(z); y=cgetg(lz,typ(z));
474 for (i=1; i<lz; i++) gel(y,i) = kbesselintern(n,gel(z,i),flag,prec);
475 return y;
476
477 case t_POLMOD:
478 y = cleanroots(gel(z,1), prec); lz = lg(y);
479 for (i=1; i<lz; i++) {
480 GEN t = poleval(gel(z,2), gel(y,i));
481 gel(y,i) = kbesselintern(n,t,flag,prec);
482 }
483 return gerepilecopy(av, y);
484
485 case t_PADIC: pari_err(impl,"p-adic kbessel function");
486 default:
487 if (!(y = toser_i(z))) break;
488 if (issmall(n,&ki))
489 {
490 k = labs(ki);
491 return gerepilecopy(av, _kbessel(k,y,flag+2,lg(y)-2,prec));
492 }
493 if (!issmall(gmul2n(n,1),&ki))
494 pari_err(talker,"cannot give a power series result in k/n bessel function");
495 k = labs(ki); n = gmul2n(stoi(k),-1);
496 fl2 = (k&3)==1;
497 pm = jbesselintern(gneg(n),y,flag,prec);
498 if (fl)
499 {
500 pp = jbesselintern(n,y,flag,prec);
501 p2 = gpowgs(y,-k); if (fl2 == 0) p2 = gneg(p2);
502 p3 = gmul2n(diviiexact(mpfact(k + 1),mpfact((k + 1) >> 1)),-(k + 1));
503 p3 = gdivgs(gmul2n(gsqr(p3),1),k);
504 p2 = gmul(p2,p3);
505 p1 = gsub(pp,gmul(p2,pm));
506 }
507 else p1 = pm;
508 return gerepileupto(av, fl2? gneg(p1): gcopy(p1));
509 }
510 pari_err(typeer,"kbessel");
511 return NULL; /* not reached */
512 }
513
514 GEN
kbesselnew(GEN n,GEN z,long prec)515 kbesselnew(GEN n, GEN z, long prec) { return kbesselintern(n,z,0,prec); }
516 GEN
nbessel(GEN n,GEN z,long prec)517 nbessel(GEN n, GEN z, long prec) { return kbesselintern(n,z,1,prec); }
518 /* J + iN */
519 GEN
hbessel1(GEN n,GEN z,long prec)520 hbessel1(GEN n, GEN z, long prec)
521 {
522 pari_sp av = avma;
523 return gerepileupto(av, gadd(jbessel(n,z,prec), mulcxI(nbessel(n,z,prec))));
524 }
525 /* J - iN */
526 GEN
hbessel2(GEN n,GEN z,long prec)527 hbessel2(GEN n, GEN z, long prec)
528 {
529 pari_sp av = avma;
530 return gerepileupto(av, gadd(jbessel(n,z,prec), mulcxmI(nbessel(n,z,prec))));
531 }
532
533 GEN
kbessel0(GEN nu,GEN gx,long flag,long prec)534 kbessel0(GEN nu, GEN gx, long flag, long prec)
535 {
536 switch(flag)
537 {
538 case 0: return kbesselnew(nu,gx,prec);
539 case 1: return kbessel(nu,gx,prec);
540 case 2: return kbessel2(nu,gx,prec);
541 default: pari_err(flagerr,"besselk");
542 }
543 return NULL; /* not reached */
544 }
545
546 /***********************************************************************/
547 /* **/
548 /** FONCTION U(a,b,z) GENERALE **/
549 /** ET CAS PARTICULIERS **/
550 /** **/
551 /***********************************************************************/
552 /* Assume gx > 0 and a,b complex */
553 /* This might one day be extended to handle complex gx */
554 /* see Temme, N. M. "The numerical computation of the confluent */
555 /* hypergeometric function U(a,b,z)" in Numer. Math. 41 (1983), */
556 /* no. 1, 63--82. */
557 GEN
hyperu(GEN a,GEN b,GEN gx,long prec)558 hyperu(GEN a, GEN b, GEN gx, long prec)
559 {
560 GEN x,y,p1,p2,p3,zf,zz,s,t,q,r,u,v,e,f,c,d,w,a1,gn;
561 long l, k, l1, n, ex;
562 pari_sp av, av1, av2;
563
564 if(gsigne(gx) <= 0) pari_err(talker,"hyperu's third argument must be positive");
565 ex = (iscomplex(a) || iscomplex(b));
566
567 l = (typ(gx)==t_REAL)? lg(gx): prec;
568 if (ex) y=cgetc(l); else y=cgetr(l);
569 l1=l+1; av=avma; x = gtofp(gx, l);
570 a1=gaddsg(1,gsub(a,b));
571 if (ex)
572 {
573 u=cgetc(l1); v=cgetc(l1); c=cgetc(l1);
574 d=cgetc(l1); e=cgetc(l1); f=cgetc(l1);
575 }
576 else
577 {
578 u=cgetr(l1); v=cgetr(l1); c=cgetr(l1);
579 d=cgetr(l1); e=cgetr(l1); f=cgetr(l1);
580 }
581 n=(long)(bit_accuracy_mul(l, LOG2) + PI*sqrt(gtodouble(gabs(gmul(a,a1),l1))));
582 av1=avma;
583 if (cmprs(x,n)<0)
584 {
585 gn=stoi(n); zf=gpow(gn,gneg_i(a),l1);
586 zz=gdivsg(-1,gn); s=gen_1; t=gen_0;
587 for (k=n-1; k>=0; k--)
588 {
589 p1=gdivgs(gmul(gmul(gaddgs(a,k),gaddgs(a1,k)),zz),k+1);
590 s=gaddsg(1,gmul(p1,s));
591 t=gadd(gaddsg(k,a),gmul(p1,t));
592 }
593 gmulz(s,zf,u); t=gmul(t,zz); gmulz(t,zf,v); avma=av1;
594 q = stor(n, l1); r=x; av1=avma;
595 for(;;)
596 {
597 p1=divsr(5,q); if (expo(p1)>= -1) p1=ghalf;
598 p2=subsr(1,divrr(r,q)); if (gcmp(p1,p2)>0) p1=p2;
599 gnegz(p1,c); avma=av1;
600 gaffsg(1,d); gaffect(u,e); gaffect(v,f);
601 p3=gsub(q,b); av2 = avma;
602 for(k=1;;k++)
603 {
604 w=gadd(gmul(gaddsg(k-1,a),u),gmul(gaddsg(1-k,p3),v));
605 gdivgsz(gmul(q,v),k,u);
606 gdivgsz(w,k,v);
607 gmulz(d,c,d);
608 gaddz(e,gmul(d,u),e); p1=gmul(d,v);
609 gaddz(f,p1,f);
610 if (gcmp0(p1) || gexpo(p1) - gexpo(f) <= 1-bit_accuracy(precision(p1))) break;
611 avma=av2;
612 }
613 p1=u; u=e; e=p1;
614 p1=v; v=f; f=p1;
615 gmulz(q,gadd(gen_1,c),q);
616 if (expo(subrr(q,r)) - expo(r) <= 1-bit_accuracy(l)) break;
617 avma=av1;
618 }
619 }
620 else
621 {
622 zf=gpow(x,gneg_i(a),l1);
623 zz=divsr(-1,x); s=gen_1;
624 for (k=n-1; k>=0; k--)
625 {
626 p1=gdivgs(gmul(gmul(gaddgs(a,k),gaddgs(a1,k)),zz),k+1);
627 s=gadd(gen_1,gmul(p1,s));
628 }
629 u = gmul(s,zf);
630 }
631 gaffect(u,y); avma=av; return y;
632 }
633
634 /* = incgam2(0, x, prec). typ(x) = t_REAL. Optimized for eint1 */
635 static GEN
incgam2_0(GEN x,GEN expx)636 incgam2_0(GEN x, GEN expx)
637 {
638 long l = lg(x), n, i;
639 GEN z;
640
641 if (expo(x) >= 4)
642 {
643 double m, mx = rtodbl(x);
644 m = (bit_accuracy_mul(l,LOG2) + mx)/4;
645 n = (long)(1+m*m/mx);
646 z = divsr(-n, addsr(n<<1,x));
647 for (i=n-1; i >= 1; i--)
648 z = divsr(-i, addrr(addsr(i<<1,x), mulsr(i,z))); /* -1 / (2 + z + x/i) */
649 return divrr(addrr(real_1(l),z), mulrr(expx, x));
650 }
651 else
652 {
653 GEN S, t, H, run = real_1(l+1);
654 n = -bit_accuracy(l)-1;
655 x = rtor(x, l+1);
656 S = z = t = H = run;
657 for (i = 2; expo(t) - expo(S) >= n; i++)
658 {
659 H = addrr(H, divrs(run,i)); /* H = sum_{i=1} 1/i */
660 z = divrs(mulrr(x,z), i); /* z = sum_{i=1} x^(i-1)/i */
661 t = mulrr(z, H); S = addrr(S, t);
662 }
663 return subrr(mulrr(x, divrr(S,expx)), addrr(mplog(x), mpeuler(l)));
664 }
665 }
666
667 /* assume x != 0 */
668 GEN
incgam2(GEN s,GEN x,long prec)669 incgam2(GEN s, GEN x, long prec)
670 {
671 GEN b, x_s, S, y;
672 long l, n, i;
673 pari_sp av = avma, av2, avlim;
674 double m,mx;
675
676 if (typ(x) != t_REAL) x = gtofp(x, prec);
677 if (gcmp0(s)) {
678 if (typ(x) == t_REAL && signe(x) > 0)
679 return gerepileuptoleaf(av, incgam2_0(x, mpexp(x)));
680 }
681 if (typ(x) == t_COMPLEX)
682 {
683 double a = rtodbl(gel(x,1));
684 double b = rtodbl(gel(x,2));
685 l = precision(x);
686 mx = sqrt(a*a + b*b);
687 }
688 else
689 {
690 l = lg(x);
691 mx = fabs(rtodbl(x));
692 }
693 m = (bit_accuracy_mul(l,LOG2) + mx)/4;
694 n = (long)(1+m*m/mx);
695 i = typ(s);
696 if (i == t_REAL) b = addsr(-1,s);
697 else
698 { /* keep b integral : final powering more efficient */
699 GEN z = gtofp(s, prec);
700 b = (i == t_INT)? addsi(-1,s): gaddsg(-1,z);
701 s = z;
702 }
703 y = gmul(gexp(gneg(x), prec), gpow(x,b,prec));
704 x_s = gsub(x, s);
705 av2 = avma; avlim = stack_lim(av2,3);
706 S = gdiv(gaddsg(-n,s), gaddgs(x_s,n<<1));
707 for (i=n-1; i>=1; i--)
708 {
709 S = gdiv(gaddsg(-i,s), gadd(gaddgs(x_s,i<<1),gmulsg(i,S)));
710 if (low_stack(avlim,stack_lim(av2,3)))
711 {
712 if(DEBUGMEM>1) pari_warn(warnmem,"incgam2");
713 S = gerepileupto(av2, S);
714 }
715 }
716 return gerepileupto(av, gmul(y, gaddsg(1,S)));
717 }
718
719 /* use exp(-x) * (x^s/s) * sum_{k >= 0} x^k / prod(i=1,k, s+i)
720 * ( = exp(-x) * x^s * sum_{k >= 0} x^k / k!(k+s) but the above is more
721 * efficient ) */
722 GEN
incgamc(GEN s,GEN x,long prec)723 incgamc(GEN s, GEN x, long prec)
724 {
725 GEN b, S, t, y;
726 long l, n, i;
727 pari_sp av = avma, av2, avlim;
728
729 if (typ(x) != t_REAL) x = gtofp(x, prec);
730 if (gcmp0(x)) return rcopy(x);
731
732 l = precision(x); n = -bit_accuracy(l)-1;
733 i = typ(s); b = s;
734 if (i != t_REAL)
735 {
736 s = gtofp(s, prec);
737 if (i != t_INT) b = s;
738 }
739 av2 = avma; avlim = stack_lim(av2,3);
740 S = t = real_1(l);
741 for (i=1; gexpo(S) >= n; i++)
742 {
743 S = gdiv(gmul(x,S), gaddsg(i,s)); /* x^i / ((s+1)...(s+i)) */
744 t = gadd(S,t);
745 if (low_stack(avlim,stack_lim(av2,3)))
746 {
747 if(DEBUGMEM>1) pari_warn(warnmem,"incgamc");
748 gerepileall(av2, 2, &S, &t);
749 }
750 }
751 y = gdiv(gmul(gexp(gneg(x),prec), gpow(x,b,prec)), s);
752 return gerepileupto(av, gmul(y,t));
753 }
754
755 /* If g != NULL, assume that g=gamma(s,prec). */
756 GEN
incgam0(GEN s,GEN x,GEN g,long prec)757 incgam0(GEN s, GEN x, GEN g, long prec)
758 {
759 pari_sp av = avma;
760 long es, e;
761 GEN z;
762
763 if (gcmp0(x)) { avma = av; return g? gcopy(g): ggamma(s,prec); }
764 es = gexpo(s); e = max(es, 0);
765 if (gsigne(real_i(s)) <= 0 || gexpo(x) > e)
766 z = incgam2(s,x,prec);
767 else
768 {
769 if (es < 0) {
770 long l = precision(s);
771 if (!l) l = prec;
772 prec = l + nbits2nlong(-es) + 1;
773 s = gtofp(s, prec);
774 x = gtofp(x, prec);
775 }
776 z = gadd(g? g: ggamma(s,prec), gneg(incgamc(s,x,prec)));
777 }
778 return gerepileupto(av, z);
779 }
780
781 GEN
incgam(GEN s,GEN x,long prec)782 incgam(GEN s, GEN x, long prec) { return incgam0(s, x, NULL, prec); }
783
784 GEN
eint1(GEN x,long prec)785 eint1(GEN x, long prec)
786 {
787 long l, n, i;
788 pari_sp av = avma;
789 GEN p1, t, S, y;
790
791 if (typ(x) != t_REAL) {
792 x = gtofp(x, prec);
793 if (typ(x) != t_REAL) pari_err(impl,"non-real argument in eint1");
794 }
795 if (signe(x) >= 0) return gerepileuptoleaf(av, incgam2_0(x, mpexp(x)));
796 /* rewritten from code contributed by Manfred Radimersky */
797 l = lg(x);
798 n = bit_accuracy(l);
799 y = negr(x);
800 if (cmprs(y, (3*n)/4) < 0) {
801 p1 = t = S = y;
802 for (i = 2; expo(t) - expo(S) >= -n; i++) {
803 p1 = mulrr(y, divrs(p1, i));
804 t = divrs(p1, i); S = addrr(S, t);
805 }
806 y = addrr(S, addrr(mplog(y), mpeuler(l)));
807 } else {
808 p1 = divsr(1, y);
809 t = S = real_1(l);
810 for (i = 1; expo(t) - expo(S) >= -n; i++) {
811 t = mulrr(p1, mulrs(t, i)); S = addrr(S, t);
812 }
813 y = mulrr(S, mulrr(p1, mpexp(y)));
814 }
815 return gerepileuptoleaf(av, negr(y));
816 }
817
818 GEN
veceint1(GEN C,GEN nmax,long prec)819 veceint1(GEN C, GEN nmax, long prec)
820 {
821 long i, n, nstop, nmin, G, chkpoint;
822 pari_sp av, av1;
823 GEN y, e1, e2, eC, F0, unr;
824
825 if (!nmax) return eint1(C,prec);
826 if (typ(nmax) != t_INT) pari_err(typeer,"veceint1");
827
828 if (signe(nmax)<=0) return cgetg(1,t_VEC);
829 if (DEBUGLEVEL>1) fprintferr("Entering veceint1:\n");
830 if (typ(C) != t_REAL || lg(C) > prec) {
831 C = gtofp(C, prec);
832 if (typ(C) != t_REAL) pari_err(typeer,"veceint1");
833 }
834 if (signe(C) <= 0) pari_err(talker,"negative or zero constant in veceint1");
835
836 n = itos(nmax); y = cgetg(n+1,t_VEC);
837 for(i=1; i<=n; i++) gel(y,i) = cgetr(prec);
838 av = avma; G = expo(C);
839 if (G >= 0) nstop = n;
840 else
841 {
842 nstop = itos(gceil(divsr(4,C))); /* >= 4 ~ 4 / C */
843 if (nstop > n) nstop = n;
844 }
845
846 eC = mpexp(C);
847 e1 = gpowgs(eC, -n);
848 e2 = gpowgs(eC, 10);
849 unr = real_1(prec);
850 av1 = avma;
851 if(DEBUGLEVEL>1) fprintferr("nstop = %ld\n",nstop);
852 if (nstop == n) goto END;
853
854 G = -bit_accuracy(prec);
855 F0 = gel(y,n); chkpoint = n;
856 affrr(eint1(mulsr(n,C),prec), F0);
857 nmin = n;
858 for(;;)
859 {
860 GEN minvn = divrs(unr,-n), My = subrr(minvn,C);
861 GEN mcn = divrs(C, -n), Mx = mcn;
862 GEN t = divrs(e1,-n), D = mkvec2( t, mulrr(My,t) );
863 long a, k, cD = 2; /* cD = #D */
864
865 /* D = [ e1/-n, (-1/n-C) * (e1/-n) ] */
866 nmin -= 10; if (nmin < nstop) nmin = nstop;
867 My = addrr(My, minvn);
868 if (DEBUGLEVEL>1 && n < chkpoint)
869 { fprintferr("%ld ",n) ; chkpoint -= nstop/20; }
870 for (a=1,n--; n>=nmin; n--,a++)
871 {
872 GEN F = F0, den = stor(-a, prec);
873 for (k=1;;)
874 {
875 GEN add;
876 if (k > cD)
877 {
878 GEN z = addrr(mulrr(My, gel(D,cD)), mulrr(Mx,gel(D,cD-1)));
879 Mx = addrr(Mx,mcn);
880 My = addrr(My,minvn);
881 D = shallowconcat(D, z); cD = k;
882 /* My = -C - k/n, Mx = -C k/n */
883 }
884 add = mulrr(den, gel(D,k));
885 if (expo(add) < G) { affrr(F,gel(y,n)); break; }
886 F = addrr(F,add); k++;
887 den = mulrs(divrs(den, k), -a);
888 /* den = prod(i=1,k, -a/i)*/
889 }
890 }
891 avma = av1; F0 = gel(y, ++n);
892 if (n <= nstop) break;
893 affrr(mulrr(e1,e2), e1);
894 }
895 END:
896 affrr(eC, e1);
897 for(i=1;; i++)
898 { /* e1 = exp(iC) */
899 affrr(incgam2_0(mulsr(i,C), e1), gel(y,i));
900 if (i == nstop) break;
901 affrr(mulrr(e1, eC), e1); avma = av1;
902 }
903 if (DEBUGLEVEL>1) fprintferr("\n");
904 avma = av; return y;
905 }
906
907 GEN
gerfc(GEN x,long prec)908 gerfc(GEN x, long prec)
909 {
910 pari_sp av;
911 GEN z, sqrtpi;
912
913 if (typ(x) != t_REAL) {
914 x = gtofp(x, prec);
915 if (typ(x) != t_REAL) pari_err(typeer,"erfc");
916 }
917 if (!signe(x)) return real_1(prec);
918 av = avma; sqrtpi = sqrtr(mppi(lg(x)));
919 z = incgam0(ghalf, gsqr(x), sqrtpi, prec);
920 z = divrr(z, sqrtpi);
921 if (signe(x) < 0) z = subsr(2,z);
922 return gerepileupto(av,z);
923 }
924
925 /***********************************************************************/
926 /** **/
927 /** FONCTION ZETA DE RIEMANN **/
928 /** **/
929 /***********************************************************************/
930 static const double log2PI = 1.83787706641;
931
932 static double
get_xinf(double beta)933 get_xinf(double beta)
934 {
935 const double maxbeta = 0.06415003; /* 3^(-2.5) */
936 double x0, y0, x1;
937
938 if (beta < maxbeta) return beta + pow(3*beta, 1.0/3.0);
939 x0 = beta + PI/2.0;
940 for(;;)
941 {
942 y0 = x0*x0;
943 x1 = (beta+atan(x0)) * (1+y0) / y0 - 1/x0;
944 if (0.99*x0 < x1) return x1;
945 x0 = x1;
946 }
947 }
948 /* optimize for zeta( s + it, prec ) */
949 static void
optim_zeta(GEN S,long prec,long * pp,long * pn)950 optim_zeta(GEN S, long prec, long *pp, long *pn)
951 {
952 double s, t, alpha, beta, n, B;
953 long p;
954 if (typ(S) == t_REAL) {
955 s = rtodbl(S);
956 t = 0.;
957 } else {
958 s = rtodbl(gel(S,1));
959 t = fabs( rtodbl(gel(S,2)) );
960 }
961
962 B = bit_accuracy_mul(prec, LOG2);
963 if (s <= 0) /* may occur if S ~ 0, and we don't use the func. eq. */
964 { /* TODO: the crude bounds below are generally valid. Optimize ? */
965 double l,l2, la = 1.; /* heuristic */
966 if (dnorm(s-1,t) < 0.1) /* |S - 1|^2 < 0.1 */
967 l2 = -(s - 0.5);
968 else
969 {
970 double rlog, ilog; dcxlog(s-1,t, &rlog,&ilog);
971 l2 = (s - 0.5)*rlog - t*ilog; /* = Re( (S - 1/2) log (S-1) ) */
972 }
973 l = (B - l2 + s*log2PI) / (2. * (1.+ log((double)la)));
974 l2 = dabs(s, t)/2;
975 if (l < l2) l = l2;
976 p = (long) ceil(l); if (p < 2) p = 2;
977 l2 = (p + s/2. - .25);
978 n = 1 + dabs(l2, t/2) * la / PI;
979 }
980 else if (t != 0)
981 {
982 double sn = dabs(s, t), L = log(sn/s);
983 alpha = B - 0.39 + L + s*(log2PI - log(sn));
984 beta = (alpha+s)/t - atan(s/t);
985 if (beta <= 0)
986 {
987 if (s >= 1.0)
988 {
989 p = 0;
990 n = exp((B - LOG2 + L) / s);
991 }
992 else
993 {
994 p = 1;
995 n = dabs(s + 1, t) / (2*PI);
996 }
997 }
998 else
999 {
1000 beta = 1.0 - s + t * get_xinf(beta);
1001 if (beta > 0)
1002 {
1003 p = (long)ceil(beta / 2.0);
1004 n = dabs(s + 2*p-1, t) / (2*PI);
1005 }
1006 else
1007 {
1008 p = 0;
1009 n = exp((B - LOG2 + L) / s);
1010 }
1011 }
1012 }
1013 else
1014 {
1015 double sn = fabs(s);
1016 beta = B + 0.61 + s*(log2PI - log(s));
1017 if (beta > 0)
1018 {
1019 p = (long)ceil(beta / 2.0);
1020 n = fabs(s + 2*p-1)/(2*PI);
1021 }
1022 else
1023 {
1024 p = 0;
1025 n = exp((B - LOG2 + log(sn/s)) / s);
1026 }
1027 }
1028 *pp = p;
1029 *pn = (long)ceil(n);
1030 if (DEBUGLEVEL) fprintferr("lim, nn: [%ld, %ld]\n", *pp, *pn);
1031 }
1032
1033 #if 0
1034 static GEN
1035 czeta(GEN s0, long prec)
1036 {
1037 int funeq = 0;
1038 long n, p, n1, i, i2;
1039 pari_sp av;
1040 GEN s, y, z, res, sig, ms, p1, p2, p3, p31, pitemp;
1041
1042 s = trans_fix_arg(&prec,&s0,&sig,&av,&res);
1043
1044 if (signe(sig) <= 0 || expo(sig) <= -2)
1045 {
1046 if (typ(s0) == t_INT)
1047 {
1048 gaffect(szeta(itos(s0), prec), res);
1049 avma = av; return res;
1050 }
1051 funeq = 1; s = gsub(gen_1,s);
1052 }
1053 optim_zeta(s, prec, &p, &n);
1054
1055 n1 = (n < 46340)? n*n: 0;
1056 y=gen_1; ms=gneg_i(s); p1=cgetr(prec+1); p2=gen_1;
1057 for (i=2; i<=n; i++)
1058 {
1059 affsr(i,p1); p2 = gexp(gmul(ms,mplog(p1)), prec+1);
1060 y = gadd(y,p2);
1061 }
1062 mpbern(p,prec); p31=cgetr(prec+1); z=gen_0;
1063 for (i=p; i>=1; i--)
1064 {
1065 i2 = i<<1;
1066 p1 = gmul(gaddsg(i2-1,s),gaddsg(i2,s));
1067 p1 = divgsns(p1, i2);
1068 p1 = n1? gdivgs(p1,n1): gdivgs(gdivgs(p1,n),n);
1069 p3 = bern(i);
1070 if (bernzone[2]>prec+1) { affrr(p3,p31); p3=p31; }
1071 z = gadd(divrs(p3,i),gmul(p1,z));
1072 }
1073 p1 = gsub(gdivsg(n,gsubgs(s,1)),ghalf);
1074 z = gmul(gadd(p1,gmul(s,gdivgs(z,n<<1))),p2);
1075 y = gadd(y,z);
1076 if (funeq)
1077 {
1078 pitemp=mppi(prec+1); setexpo(pitemp,2);
1079 y = gmul(gmul(y,ggamma(s,prec+1)), gpow(pitemp,ms,prec+1));
1080 setexpo(pitemp, 0);
1081 y = gmul2n(gmul(y, gcos(gmul(pitemp,s),prec+1)),1);
1082 }
1083 gaffect(y,res); avma = av; return res;
1084 }
1085 #endif
1086
1087 /* 1/zeta(n) using Euler product. Assume n > 0.
1088 * if (lba != 0) it is log(bit_accuracy) we _really_ require */
1089 GEN
inv_szeta_euler(long n,double lba,long prec)1090 inv_szeta_euler(long n, double lba, long prec)
1091 {
1092 GEN z, res = cgetr(prec);
1093 pari_sp av = avma, avlim = stack_lim(av, 1);
1094 byteptr d = diffptr + 2;
1095 double A = n / (LOG2*BITS_IN_LONG), D;
1096 ulong p, lim;
1097
1098 if (n > bit_accuracy(prec)) return real_1(prec);
1099 if (!lba) lba = bit_accuracy_mul(prec, LOG2);
1100 D = exp((lba - log(n-1)) / (n-1));
1101 lim = 1 + (ulong)ceil(D);
1102 maxprime_check(lim);
1103
1104 prec++;
1105 z = gsub(gen_1, real2n(-n, prec));
1106 for (p = 3; p <= lim;)
1107 {
1108 long l = prec + 1 - (long)floor(A * log(p));
1109 GEN h;
1110
1111 if (l < 3) l = 3;
1112 else if (l > prec) l = prec;
1113 h = divrr(z, rpowuu(p, (ulong)n, l));
1114 z = subrr(z, h);
1115 if (low_stack(avlim, stack_lim(av,1)))
1116 {
1117 if (DEBUGMEM>1) pari_warn(warnmem,"inv_szeta_euler, p = %lu/%lu", p,lim);
1118 affrr(z, res); avma = av;
1119 }
1120 NEXT_PRIME_VIADIFF(p,d);
1121 }
1122 affrr(z, res); avma = av; return res;
1123 }
1124
1125 /* assume n even > 0, if iz != NULL, assume iz = 1/zeta(n) */
1126 GEN
bernreal_using_zeta(long n,GEN iz,long prec)1127 bernreal_using_zeta(long n, GEN iz, long prec)
1128 {
1129 long l = prec + 1;
1130 GEN z;
1131
1132 if (!iz) iz = inv_szeta_euler(n, 0., l);
1133 z = divrr(mpfactr(n, l), mulrr(gpowgs(Pi2n(1, l), n), iz));
1134 setexpo(z, expo(z) + 1); /* 2 * n! * zeta(n) / (2Pi)^n */
1135 if ((n & 3) == 0) setsigne(z, -1);
1136 return z;
1137 }
1138
1139 /* assume n even > 0. Faster than standard bernfrac for n >= 6 */
1140 GEN
bernfrac_using_zeta(long n)1141 bernfrac_using_zeta(long n)
1142 {
1143 pari_sp av = avma;
1144 GEN iz, a, d, D = divisors(utoipos( n/2 ));
1145 long i, prec, l = lg(D);
1146 double t, u;
1147
1148 d = utoipos(6); /* 2 * 3 */
1149 for (i = 2; i < l; i++) /* skip 1 */
1150 { /* Clausen - von Staudt */
1151 ulong p = 2*itou(gel(D,i)) + 1;
1152 if (uisprime(p)) d = muliu(d, p);
1153 }
1154 /* 1.712086 = ??? */
1155 t = log( gtodouble(d) ) + (n + 0.5) * log(n) - n*(1+log2PI) + 1.712086;
1156 u = t / (LOG2*BITS_IN_LONG); prec = (long)ceil(u);
1157 prec += 3;
1158 iz = inv_szeta_euler(n, t, prec);
1159 a = roundr( mulir(d, bernreal_using_zeta(n, iz, prec)) );
1160 return gerepilecopy(av, mkfrac(a, d));
1161 }
1162
1163 /* y = binomial(n,k-2). Return binomial(n,k) */
1164 static GEN
next_bin(GEN y,long n,long k)1165 next_bin(GEN y, long n, long k)
1166 {
1167 y = divrs(mulrs(y, n-k+2), k-1);
1168 return divrs(mulrs(y, n-k+1), k);
1169 }
1170
1171 /* assume k > 1 odd */
1172 static GEN
szeta_odd(long k,long prec)1173 szeta_odd(long k, long prec)
1174 {
1175 long kk, n, li = -(1+bit_accuracy(prec));
1176 pari_sp av = avma, av2, limit;
1177 GEN y, p1, qn, z, q, pi2 = Pi2n(1, prec), binom= real_1(prec+1);
1178
1179 q = mpexp(pi2); kk = k+1; /* >= 4 */
1180 y = NULL; /* gcc -Wall */
1181 if ((k&3)==3)
1182 {
1183 for (n=0; n <= kk>>1; n+=2)
1184 {
1185 p1 = mulrr(bernreal(kk-n,prec),bernreal(n,prec));
1186 if (n) { binom = next_bin(binom,kk,n); setlg(binom,prec+1); }
1187 p1 = mulrr(binom,p1);
1188 if (n == kk>>1) setexpo(p1, expo(p1)-1);
1189 if ((n>>1)&1) setsigne(p1,-signe(p1));
1190 y = n? addrr(y,p1): p1;
1191 }
1192 y = mulrr(divrr(gpowgs(pi2,k),mpfactr(kk,prec)),y);
1193
1194 av2 = avma; limit = stack_lim(av2,1);
1195 qn = gsqr(q); z = ginv( addrs(q,-1) );
1196 for (n=2; ; n++)
1197 {
1198 p1 = ginv( mulir(powuu(n,k),addrs(qn,-1)) );
1199
1200 z = addrr(z,p1); if (expo(p1)< li) break;
1201 qn = mulrr(qn,q);
1202 if (low_stack(limit,stack_lim(av2,1)))
1203 {
1204 if (DEBUGMEM>1) pari_warn(warnmem,"szeta");
1205 gerepileall(av2,2, &z, &qn);
1206 }
1207 }
1208 setexpo(z, expo(z)+1);
1209 y = addrr(y,z); setsigne(y,-signe(y));
1210 }
1211 else
1212 {
1213 GEN p2 = divrs(pi2, k-1);
1214 for (n=0; n <= k>>1; n+=2)
1215 {
1216 p1 = mulrr(bernreal(kk-n,prec),bernreal(n,prec));
1217 if (n) binom = next_bin(binom,kk,n);
1218 p1 = mulrr(binom,p1);
1219 p1 = mulsr(kk-(n<<1),p1);
1220 if ((n>>1)&1) setsigne(p1,-signe(p1));
1221 y = n? addrr(y,p1): p1;
1222 }
1223 y = mulrr(divrr(gpowgs(pi2,k),mpfactr(kk,prec)),y);
1224 y = divrs(y,k-1);
1225 av2 = avma; limit = stack_lim(av2,1);
1226 qn = q; z=gen_0;
1227 for (n=1; ; n++)
1228 {
1229 p1=mulir(powuu(n,k),gsqr(addrs(qn,-1)));
1230 p1=divrr(addrs(mulrr(qn,addsr(1,mulsr(n<<1,p2))),-1),p1);
1231
1232 z=addrr(z,p1); if (expo(p1) < li) break;
1233 qn=mulrr(qn,q);
1234 if (low_stack(limit,stack_lim(av2,1)))
1235 {
1236 if (DEBUGMEM>1) pari_warn(warnmem,"szeta");
1237 gerepileall(av2,2, &z, &qn);
1238 }
1239 }
1240 setexpo(z, expo(z)+1);
1241 y = subrr(y,z);
1242 }
1243 return gerepileuptoleaf(av, y);
1244 }
1245
1246 /* assume k > 0 even. Return B_k */
1247 static GEN
single_bern(long k,long prec)1248 single_bern(long k, long prec)
1249 {
1250 pari_sp av;
1251 GEN B;
1252 if (OK_bern(k >> 1, prec)) B = bernreal(k, prec);
1253 else if (k * (log(k) - 2.83) > bit_accuracy_mul(prec, LOG2))
1254 B = bernreal_using_zeta(k, NULL, prec);
1255 else
1256 {
1257 B = cgetr(prec);
1258 av = avma; gaffect(bernfrac(k), B);
1259 avma = av;
1260 }
1261 return B;
1262 }
1263
1264 /* assume k != 1 */
1265 GEN
szeta(long k,long prec)1266 szeta(long k, long prec)
1267 {
1268 pari_sp av = avma;
1269 GEN y;
1270
1271 /* treat trivial cases */
1272 if (!k) { y = real2n(-1, prec); setsigne(y,-1); return y; }
1273 if (k < 0)
1274 {
1275 if ((k&1) == 0) return gen_0;
1276 /* the one value such that k < 0 and 1 - k < 0, due to overflow */
1277 if ((ulong)k == (HIGHBIT | 1))
1278 pari_err(talker, "too large negative arg %ld in zeta", k);
1279 k = 1-k;
1280 y = single_bern(k, prec);
1281 return gerepileuptoleaf(av, divrs(y, -k));
1282 }
1283 if (k > bit_accuracy(prec)+1) return real_1(prec);
1284 if ((k&1) == 0)
1285 {
1286 if (!OK_bern(k >> 1, prec)
1287 && (k * (log(k) - 2.83) > bit_accuracy_mul(prec, LOG2)))
1288 y = ginv( inv_szeta_euler(k, 0, prec) ); /* would use zeta above */
1289 else
1290 {
1291 y = mulrr(gpowgs(Pi2n(1, prec), k), single_bern(k, prec));
1292 y = divrr(y, mpfactr(k,prec));
1293 y[1] = evalsigne(1) | evalexpo(expo(y)-1);
1294 }
1295 return gerepileuptoleaf(av, y);
1296 }
1297 /* k > 1 odd */
1298 if (k * log(k) > bit_accuracy_mul(prec, LOG2)) /* heuristic */
1299 return gerepileuptoleaf(av, ginv( inv_szeta_euler(k, 0, prec) ));
1300 return szeta_odd(k, prec);
1301 }
1302
1303 /* return x^n, assume n > 0 */
1304 static long
pows(long x,long n)1305 pows(long x, long n)
1306 {
1307 long i, y = x;
1308 for (i=1; i<n; i++) y *= x;
1309 return y;
1310 }
1311
1312 /* return n^-s, n > 1 odd. tab[q] := q^-s, q prime power */
1313 static GEN
n_s(ulong n,GEN * tab)1314 n_s(ulong n, GEN *tab)
1315 {
1316 byteptr d = diffptr + 2;
1317 GEN x = NULL;
1318 long p, e;
1319
1320 for (p = 3; n > 1; )
1321 {
1322 e = u_lvalrem(n, p, &n);
1323 if (e)
1324 {
1325 GEN y = tab[pows(p,e)];
1326 if (!x) x = y; else x = gmul(x,y);
1327 }
1328 NEXT_PRIME_VIADIFF_CHECK(p,d);
1329 }
1330 return x;
1331 }
1332
1333 /* s0 a t_INT, t_REAL or t_COMPLEX.
1334 * If a t_INT, assume it's not a trivial case (i.e we have s0 > 1, odd)
1335 * */
1336 GEN
czeta(GEN s0,long prec)1337 czeta(GEN s0, long prec)
1338 {
1339 GEN s, u, a, y, res, tes, sig, invn2, unr;
1340 GEN sim, *tab, tabn;
1341 ulong p, sqn;
1342 long i, nn, lim, lim2, ct;
1343 pari_sp av, av2 = avma, avlim;
1344 int funeq = 0;
1345 byteptr d;
1346
1347 if (DEBUGLEVEL>2) (void)timer2();
1348 s = trans_fix_arg(&prec,&s0,&sig,&av,&res);
1349 if (gcmp0(s)) { y = gneg(ghalf); goto END; }
1350 if (gexpo(gsub(s, gen_1)) < -5 ||
1351 (gexpo(s) > -5 && (signe(sig) <= 0 || expo(sig) < -1)))
1352 { /* s <--> 1-s */
1353 if (typ(s0) == t_INT)
1354 {
1355 i = itos(s0); avma = av2;
1356 return szeta(i, prec);
1357 }
1358 funeq = 1; s = gsub(gen_1, s); sig = real_i(s);
1359 }
1360 if (gcmpgs(sig, bit_accuracy(prec) + 1) > 0) { y = gen_1; goto END; }
1361 optim_zeta(s, prec, &lim, &nn);
1362 maxprime_check((ulong)nn);
1363 prec++; unr = real_1(prec); /* one extra word of precision */
1364
1365 tab = (GEN*)cgetg(nn, t_VEC); /* table of q^(-s), q = p^e */
1366 d = diffptr + 1;
1367 if (typ(s0) == t_INT)
1368 { /* no explog for 1/p^s */
1369 ulong k = itou(s0);
1370 for (p=2; p < (ulong)nn;) {
1371 tab[p] = divrr(unr, rpowuu(p, k, prec));
1372 NEXT_PRIME_VIADIFF(p,d);
1373 }
1374 a = divrr(unr, rpowuu((ulong)nn, k, prec));
1375 }
1376 else
1377 { /* general case */
1378 GEN ms = gneg(s), rp = cgetr(prec);
1379 for (p=2; p < (ulong)nn;)
1380 {
1381 affur(p, rp);
1382 tab[p] = gexp(gmul(ms, mplog(rp)), prec);
1383 NEXT_PRIME_VIADIFF(p,d);
1384 }
1385 affsr(nn, rp);
1386 a = gexp(gmul(ms, mplog(rp)), prec);
1387 }
1388 sqn = (ulong)sqrt(nn-1.); maxprime_check(sqn);
1389 d = diffptr + 2; /* fill in odd prime powers */
1390 for (p=3; p <= sqn; )
1391 {
1392 ulong oldq = p, q = p*p;
1393 while (q<(ulong)nn) { tab[q] = gmul(tab[p], tab[oldq]); oldq = q; q *= p; }
1394 NEXT_PRIME_VIADIFF(p,d);
1395 }
1396 if (DEBUGLEVEL>2) msgtimer("tab[q^-s] from 1 to N-1");
1397
1398 tabn = cgetg(nn, t_VECSMALL); ct = 0;
1399 for (i = nn-1; i; i>>=1) tabn[++ct] = (i-1)>>1;
1400 sim = y = unr;
1401 for (i=ct; i > 1; i--)
1402 {
1403 long j;
1404 pari_sp av2 = avma;
1405 for (j=tabn[i]+1; j<=tabn[i-1]; j++)
1406 sim = gadd(sim, n_s(2*j+1, tab));
1407 sim = gerepileupto(av2, sim);
1408 y = gadd(sim, gmul(tab[2],y));
1409 }
1410 y = gadd(y, gmul2n(a,-1));
1411 if (DEBUGLEVEL>2) msgtimer("sum from 1 to N-1");
1412
1413 invn2 = divri(unr, mulss(nn,nn)); lim2 = lim<<1;
1414 tes = bernreal(lim2, prec);
1415 if (typ(s0) == t_INT)
1416 {
1417 av2 = avma; avlim = stack_lim(av2,3);
1418 for (i=lim2-2; i>=2; i-=2)
1419 { /* using single prec (when (s0 + i) < 2^31) not faster (even at \p28) */
1420 u = mulri(mulrr(tes,invn2), mulii(addsi(i,s0), addsi(i-1,s0)));
1421 tes = addrr(bernreal(i,prec), divrsns(u, i+1)); /* u / (i+1)(i+2) */
1422 if (low_stack(avlim,stack_lim(av2,3)))
1423 {
1424 if(DEBUGMEM>1) pari_warn(warnmem,"czeta");
1425 tes = gerepileuptoleaf(av2, tes);
1426 }
1427 }
1428 u = gmul(gmul(tes,invn2), gmul2n(mulii(s0, addsi(-1,s0)), -1));
1429 tes = gmulsg(nn, gaddsg(1, u));
1430 }
1431 else /* typ(s0) != t_INT */
1432 {
1433 GEN s1, s2, s3, s4, s5;
1434 s1 = gsub(gmul2n(s,1), unr);
1435 s2 = gmul(s, gsub(s,unr));
1436 s3 = gmul2n(invn2,3);
1437 av2 = avma; avlim = stack_lim(av2,3);
1438 s4 = gmul(invn2, gmul2n(gaddsg(4*lim-2,s1),1));
1439 s5 = gmul(invn2, gadd(s2, gmulsg(lim2, gaddgs(s1, lim2))));
1440 for (i = lim2-2; i>=2; i -= 2)
1441 {
1442 s5 = gsub(s5, s4);
1443 s4 = gsub(s4, s3);
1444 tes = gadd(bernreal(i,prec), divgsns(gmul(s5,tes), i+1));
1445 if (low_stack(avlim,stack_lim(av2,3)))
1446 {
1447 if(DEBUGMEM>1) pari_warn(warnmem,"czeta");
1448 gerepileall(av2,3, &tes,&s5,&s4);
1449 }
1450 }
1451 u = gmul(gmul(tes,invn2), gmul2n(s2, -1));
1452 tes = gmulsg(nn, gaddsg(1, u));
1453 }
1454 if (DEBUGLEVEL>2) msgtimer("Bernoulli sum");
1455 /* y += tes n^(-s) / (s-1) */
1456 y = gadd(y, gmul(tes, gdiv(a, gsub(s, unr))));
1457
1458 END:
1459 if (funeq)
1460 {
1461 y = gmul(gmul(y, ggamma(gprec_w(s,prec),prec)),
1462 gpow(Pi2n(1,prec), gneg(s), prec));
1463 y = gmul2n(gmul(y, gcos(gmul(Pi2n(-1,prec),s), prec)), 1);
1464 }
1465 gaffect(y,res); avma = av; return res;
1466 }
1467
1468 /* return P mod x^n where P is polynomial in x */
1469 static GEN
pol_mod_xn(GEN P,long n)1470 pol_mod_xn(GEN P, long n)
1471 {
1472 long j;
1473 GEN R = cgetg(n+2, t_POL);
1474 R[1] = evalvarn(0);
1475 for (j = 0; j < n; j++)
1476 R[j+2] = (long)polcoeff0(P, j, 0);
1477 return normalizepol_i(R, n+2);
1478 }
1479
1480 /* compute the values of the twisted partial
1481 zeta function Z_f(a, c, s) for a in va */
1482 GEN
twistpartialzeta(GEN p,GEN q,long f,long c,GEN va,GEN cff)1483 twistpartialzeta(GEN p, GEN q, long f, long c, GEN va, GEN cff)
1484 {
1485 long j, k, lva = lg(va)-1, N = lg(cff)-1;
1486 pari_sp av, av2, lim;
1487 GEN Ax, Cx, Bx, Dx, x = pol_x[0], y = pol_x[fetch_user_var("y")], eta;
1488 GEN cyc, psm, rep;
1489
1490 cyc = gdiv(gsubgs(gpowgs(y, c), 1), gsubgs(y, 1));
1491 psm = polsym(cyc, degpol(cyc) - 1);
1492 eta = gmodulo(y, cyc);
1493 av = avma;
1494 Ax = gsubgs(gpowgs(gaddgs(x, 1), f), 1);
1495 Ax = gdiv(gmul(Ax, gpowgs(eta, f)), gsubsg(1, gpowgs(eta, f)));
1496 Ax = gerepileupto(av, RgX_to_FqX(Ax, cyc, q));
1497 Cx = Ax;
1498 Bx = gen_1;
1499 av = avma; lim = stack_lim(av, 1);
1500 for (j = 2; j <= N; j++)
1501 {
1502 Bx = gadd(Bx, Cx);
1503 Bx = FpXQX_red(Bx, cyc, q);
1504 Cx = FpXQX_mul(Cx, Ax, cyc, q);
1505 Cx = pol_mod_xn(Cx, N);
1506 if (gcmp0(Cx)) break;
1507 if (low_stack(lim, stack_lim(av, 1)))
1508 {
1509 if(DEBUGMEM>1) pari_warn(warnmem, "twistpartialzeta (1), j = %ld/%ld", j, N);
1510 gerepileall(av, 2, &Cx, &Bx);
1511 }
1512 }
1513 Bx = lift(gmul(ginv(gsubsg(1, gpowgs(eta, f))), Bx));
1514 Bx = gerepileupto(av, RgX_to_FqX(Bx, cyc, q));
1515 Cx = lift(gmul(eta, gaddsg(1, x)));
1516 Dx = pol_1[varn(x)];
1517 av2 = avma; lim = stack_lim(av2, 1);
1518 for (j = lva; j > 1; j--)
1519 {
1520 GEN Ex;
1521 long e = va[j] - va[j-1];
1522 if (e == 1)
1523 Ex = Cx;
1524 else
1525 /* e is very small in general and actually very rarely different
1526 to 1, it is always 1 for zetap (so it should be OK not to store
1527 them or to compute them in a smart way) */
1528 Ex = gpowgs(Cx, e);
1529 Dx = gaddsg(1, FpXQX_mul(Dx, Ex, cyc, q));
1530 if (low_stack(lim, stack_lim(av2, 1)))
1531 {
1532 if(DEBUGMEM>1)
1533 pari_warn(warnmem, "twistpartialzeta (2), j = %ld/%ld", lva-j, lva);
1534 Dx = gerepileupto(av2, FpXQX_red(Dx, cyc, q));
1535 }
1536 }
1537 Dx = FpXQX_mul(Dx, Cx, cyc, q); /* va[1] = 1 */
1538 Bx = gerepileupto(av, FpXQX_mul(Dx, Bx, cyc, q));
1539 rep = gen_0;
1540 av2 = avma; lim = stack_lim(av2, 1);
1541 for (k = 1; k <= N; k++)
1542 {
1543 GEN p2, ak = polcoeff_i(Bx, k, 0);
1544 p2 = quicktrace(ak, psm);
1545 rep = modii(addii(rep, mulii(gel(cff, k), p2)), q);
1546 if (low_stack(lim, stack_lim(av2, 1)))
1547 {
1548 if(DEBUGMEM>1) pari_warn(warnmem, "twistpartialzeta (3), j = %ld/%ld", k, N);
1549 rep = gerepileupto(av2, rep);
1550 }
1551 }
1552 return rep;
1553 }
1554
1555 #if 0
1556 /* initialize the roots of unity for the computation
1557 of the Teichmuller character (also the values of f and c) */
1558 GEN
1559 init_teich(ulong p, GEN q, long prec)
1560 {
1561 GEN vz, gp = utoipos(p);
1562 pari_sp av = avma;
1563 long j;
1564
1565 if (p == 2UL)
1566 return NULL;
1567 else
1568 { /* primitive (p-1)-th root of 1 */
1569 GEN z, z0 = padicsqrtnlift(gen_1, utoipos(p-1), gener_Fp(gp), gp, prec);
1570 z = z0;
1571 vz = cgetg(p, t_VEC);
1572 for (j = 1; j < (long)p-2; j++)
1573 {
1574 gel(vz, umodiu(z, p)) = z; /* z = z0^i */
1575 z = modii(mulii(z, z0), q);
1576 }
1577 gel(vz, umodiu(z, p)) = z; /* z = z0^(p-2) */
1578 gel(vz,1) = gen_1; /* z0^(p-1) */
1579 }
1580 return gerepileupto(av, gcopy(vz));
1581 }
1582
1583 /* compute phi^(m)_s(x); s must be an integer */
1584 GEN
1585 phi_ms(ulong p, GEN q, long m, GEN s, long x, GEN vz)
1586 {
1587 long xp = x % p;
1588 GEN p1, p2;
1589
1590 if (!xp) return gen_0;
1591 if (vz)
1592 p1 =gel(vz,xp); /* vz[x] = Teichmuller(x) */
1593 else
1594 p1 = (x & 2)? gen_m1: gen_1;
1595 p1 = Fp_pow(p1, addis(s, m), q);
1596 p2 = Fp_pow(stoi(x), negi(s), q);
1597 return modii(mulii(p1,p2), q);
1598 }
1599
1600 /* compute the first N coefficients of the Mahler expansion
1601 of phi^m_s skipping the first one (which is zero) */
1602 GEN
1603 coeff_of_phi_ms(ulong p, GEN q, long m, GEN s, long N, GEN vz)
1604 {
1605 GEN qs2 = shifti(q, -1), cff = zerovec(N);
1606 pari_sp av, lim;
1607 long k, j;
1608
1609 av = avma; lim = stack_lim(av, 2);
1610 for (k = 1; k <= N; k++)
1611 {
1612 gel(cff, k) = phi_ms(p, q, m, s, k, vz);
1613 if (low_stack(lim, stack_lim(av, 2)))
1614 {
1615 if(DEBUGMEM>1)
1616 pari_warn(warnmem, "coeff_of_phi_ms (1), k = %ld/%ld", N-k, N);
1617 cff = gerepileupto(av, gcopy(cff));
1618 }
1619 }
1620 for (j = N; j > 1; j--)
1621 {
1622 GEN b = subii(gel(cff, j), gel(cff, j-1));
1623 gel(cff, j) = centermodii(b, q, qs2);
1624 if (low_stack(lim, stack_lim(av, 2)))
1625 {
1626 if(DEBUGMEM>1)
1627 pari_warn(warnmem, "coeff_of_phi_ms (2), j = %ld/%ld", N-j, N);
1628 cff = gerepileupto(av, gcopy(cff));
1629 }
1630 }
1631 for (k = 1; k < N; k++)
1632 for (j = N; j > k; j--)
1633 {
1634 GEN b = subii(gel(cff, j), gel(cff, j-1));
1635 gel(cff, j) = centermodii(b, q, qs2);
1636 if (low_stack(lim, stack_lim(av, 2)))
1637 {
1638 if(DEBUGMEM>1)
1639 pari_warn(warnmem, "coeff_of_phi_ms (3), (k,j) = (%ld,%ld)/%ld",
1640 k, N-j, N);
1641 cff = gerepileupto(av, gcopy(cff));
1642 }
1643 }
1644 k = N; while(gcmp0(gel(cff, k))) k--;
1645 setlg(cff, k+1);
1646 if (DEBUGLEVEL > 2)
1647 fprintferr(" coeff_of_phi_ms: %ld coefficients kept out of %ld\n",
1648 k, N);
1649 return gerepileupto(av, cff);
1650 }
1651
1652 static long
1653 valfact(long N, ulong p)
1654 {
1655 long f = 0;
1656 while (N > 1)
1657 {
1658 N /= p;
1659 f += N;
1660 }
1661 return f;
1662 }
1663
1664 static long
1665 number_of_terms(ulong p, long prec)
1666 {
1667 long N, f;
1668
1669 if (prec == 0) return p;
1670 N = (long)((p-1)*prec + (p>>1)*(log2(prec)/log2(p)));
1671 N = p*(N/p);
1672 f = valfact(N, p);
1673 while (f > prec)
1674 {
1675 N = p*(N/p) - 1;
1676 f -= u_lval(N+1, p);
1677 }
1678 while (f < prec)
1679 {
1680 N = p*(N/p+1);
1681 f += u_lval(N, p);
1682 }
1683 return N;
1684 }
1685
1686 static GEN
1687 zetap(GEN s)
1688 {
1689 ulong p;
1690 long N, f, c, prec = precp(s);
1691 pari_sp av = avma;
1692 GEN gp, q, vz, is, cff, val, va, cft;
1693
1694 if (valp(s) < 0)
1695 pari_err(talker, "argument must be a p-adic integer");
1696
1697 gp = gel(s,2); p = itou(gp);
1698 is = gtrunc(s); /* make s an integer */
1699
1700 N = number_of_terms(p, prec? prec:1);
1701 q = gpowgs(gp, prec? prec:1);
1702
1703 /* initialize the roots of unity for the computation
1704 of the Teichmuller character (also the values of f and c) */
1705 if (DEBUGLEVEL > 1) fprintferr("zetap: computing (p-1)th roots of 1\n");
1706 vz = init_teich(p, q, prec? prec:1);
1707 if (p == 2UL) { f = 4; c = 3; } else { f = (long)p; c = 2; }
1708
1709 /* compute the first N coefficients of the Mahler expansion
1710 of phi^(-1)_s skipping the first one (which is zero) */
1711 if (DEBUGLEVEL > 1)
1712 fprintferr("zetap: computing Mahler expansion of phi^(-1)_s\n");
1713 cff = coeff_of_phi_ms(p, q, -1, is, N, vz);
1714
1715 /* compute the coefficients of the power series corresponding
1716 to the twisted partial zeta function Z_f(a, c, s) for a in va */
1717 /* The line below looks a bit stupid but it is to keep the
1718 possibility of later adding p-adic Dirichlet L-functions */
1719 va = perm_identity(f - 1);
1720 if (DEBUGLEVEL > 1)
1721 fprintferr("zetap: computing values of twisted partial zeta functions\n");
1722 val = twistpartialzeta(gp, q, f, c, va, cff);
1723
1724 /* sum over all a's the coefficients of the twisted
1725 partial zeta functions and integrate */
1726 if (DEBUGLEVEL > 1)
1727 fprintferr("zetap: multiplying by correcting factor\n");
1728
1729 /* multiply by the corrective factor */
1730 cft = gsubgs(gmulsg(c, phi_ms(p, q, -1, is, c, vz)), 1);
1731 val = gdiv(val, cft);
1732
1733 /* adjust the precision and return */
1734 return gerepileupto(av, cvtop(val, gp, prec));
1735 }
1736 #else
1737 static GEN
hurwitz_p(GEN cache,GEN s,GEN x,GEN p,long prec)1738 hurwitz_p(GEN cache, GEN s, GEN x, GEN p, long prec)
1739 {
1740 GEN S, x2, x2j, s_1 = gsubgs(s,1);
1741 long j, J = lg(cache)-2;
1742 x = ginv(gadd(x, zeropadic(p, prec)));
1743 x2 = gsqr(x);
1744 S = gmul2n(gmul(s_1, x), -1);
1745 x2j = gen_1;
1746 for (j = 0;; j++)
1747 {
1748 S = gadd(S, gmul(gel(cache, j+1), x2j));
1749 if (j == J) break;
1750 x2j = gmul(x2, x2j);
1751 }
1752 return gmul(gdiv(S, s_1), gexp(gmul(s_1, glog(x, 0)), 0));
1753 }
1754
1755 static GEN
init_cache(long J,GEN s)1756 init_cache(long J, GEN s)
1757 {
1758 GEN C = gen_1, cache = bernvec(J);
1759 long j;
1760
1761 for (j = 1; j <= J; j++)
1762 { /* B_{2j} * binomial(1-s, 2j) */
1763 GEN t = gmul(gaddgs(s, 2*j-3), gaddgs(s, 2*j-2));
1764 C = gdiv(gmul(C, t), mulss(2*j, 2*j-1));
1765 gel(cache, j+1) = gmul(gel(cache, j+1), C);
1766 }
1767 return cache;
1768 }
1769
1770 static GEN
zetap(GEN s)1771 zetap(GEN s)
1772 {
1773 pari_sp av = avma;
1774 GEN cache, S, gp = gel(s,2);
1775 ulong a, p = itou(gp);
1776 long J, prec = valp(s) + precp(s);
1777
1778 if (prec <= 0) prec = 1;
1779 if (p == 2) {
1780 J = ((long)(1+ceil((prec+1.)/2))) >> 1;
1781 cache = init_cache(J, s);
1782 S = gmul2n(hurwitz_p(cache, s, gmul2n(gen_1, -2), gen_2, prec), -1);
1783 } else {
1784 J = (prec+2) >> 1;
1785 cache = init_cache(J, s);
1786 S = gen_0;
1787 for (a = 1; a <= (p-1)>>1; a++)
1788 S = gadd(S, hurwitz_p(cache, s, gdivsg(a, gp), gp, prec));
1789 S = gdiv(gmul2n(S, 1), gp);
1790 }
1791 return gerepileupto(av, S);
1792 }
1793 #endif
1794
1795 GEN
gzeta(GEN x,long prec)1796 gzeta(GEN x, long prec)
1797 {
1798 if (gcmp1(x)) pari_err(talker, "argument equal to one in zeta");
1799 switch(typ(x))
1800 {
1801 case t_INT:
1802 if (is_bigint(x))
1803 {
1804 if (signe(x) > 0) return real_1(prec);
1805 if (signe(x) < 0 && mod2(x) == 0) return real_0(prec);
1806 }
1807 return szeta(itos(x),prec);
1808
1809 case t_REAL: case t_COMPLEX:
1810 return czeta(x,prec);
1811
1812 case t_INTMOD: pari_err(typeer,"gzeta");
1813
1814 case t_PADIC:
1815 return zetap(x);
1816
1817 case t_SER: pari_err(impl,"zeta of power series");
1818 }
1819 return transc(gzeta,x,prec);
1820 }
1821
1822 /***********************************************************************/
1823 /** **/
1824 /** FONCTIONS POLYLOGARITHME **/
1825 /** **/
1826 /***********************************************************************/
1827
1828 /* validity domain contains .005 < |x| < 230
1829 * Li_m(x = e^z) = sum_n=0 zeta(m-n) z^n / n!
1830 * with zeta(1) := H_m - log(-z) */
1831 static GEN
cxpolylog(long m,GEN x,long prec)1832 cxpolylog(long m, GEN x, long prec)
1833 {
1834 long li, i, n, bern_upto;
1835 pari_sp av=avma;
1836 GEN p1,z,h,q,s;
1837 int real;
1838
1839 if (gcmp1(x)) return szeta(m,prec);
1840 real = typ(x) == t_REAL;
1841
1842 z = glog(x,prec); h = gen_1;
1843 for (i=2; i<m; i++) h = gadd(h, ginv(utoipos(i)));
1844 h = gadd(h, gneg_i( glog(gneg_i(z),prec) ));
1845
1846 bern_upto = m+50; mpbern(bern_upto,prec);
1847 q = gen_1; s = szeta(m,prec);
1848 for (n=1; n<=m+1; n++)
1849 {
1850 q = gdivgs(gmul(q,z),n);
1851 if (n == m-1) {
1852 p1 = gmul(h, q);
1853 if (real) p1 = real_i(p1);
1854 }
1855 else
1856 p1 = gmul(szeta(m-n,prec), real? real_i(q): q);
1857 s = gadd(s, p1);
1858 }
1859
1860 z = gsqr(z); li = -(bit_accuracy(prec)+1);
1861 for(n = m+3;; n += 2)
1862 {
1863 GEN zet = szeta(m-n,prec);
1864 q = divgsns(gmul(q,z), n-1);
1865 s = gadd(s, gmul(zet, real? real_i(q): q));
1866 if (gexpo(q) + expo(zet) < li) break;
1867 if (n >= bern_upto) { bern_upto += 50; mpbern(bern_upto,prec); }
1868 }
1869 return gerepileupto(av,s);
1870 }
1871
1872 GEN
polylog(long m,GEN x,long prec)1873 polylog(long m, GEN x, long prec)
1874 {
1875 long l, e, i, G, sx;
1876 pari_sp av, av1, limpile;
1877 GEN X, Xn, z, p1, p2, y;
1878
1879 if (m<0) pari_err(talker,"negative index in polylog");
1880 if (!m) return gneg(ghalf);
1881 if (gcmp0(x)) return gcopy(x);
1882 av = avma;
1883 if (m==1)
1884 return gerepileupto(av, gneg(glog(gsub(gen_1,x), prec)));
1885
1886 l = precision(x);
1887 if (!l) { l=prec; x=gmul(x, real_1(l)); }
1888 e = gexpo(gnorm(x)); if (!e || e== -1) return cxpolylog(m,x,prec);
1889 X = (e > 0)? ginv(x): x;
1890 G = -bit_accuracy(l);
1891 av1=avma; limpile=stack_lim(av1,1);
1892 y = Xn = X;
1893 for (i=2; ; i++)
1894 {
1895 Xn = gmul(X,Xn); p2 = gdiv(Xn,powuu(i,m));
1896 y = gadd(y,p2);
1897 if (gexpo(p2) <= G) break;
1898
1899 if (low_stack(limpile, stack_lim(av1,1)))
1900 {
1901 if(DEBUGMEM>1) pari_warn(warnmem,"polylog");
1902 gerepileall(av1,2, &y, &Xn);
1903 }
1904 }
1905 if (e < 0) return gerepileupto(av, y);
1906
1907 sx = gsigne(imag_i(x));
1908 if (!sx)
1909 {
1910 if (m&1) sx = gsigne(gsub(gen_1, real_i(x)));
1911 else sx = - gsigne(real_i(x));
1912 }
1913 z = pureimag( divri(mppi(l), mpfact(m-1)) );
1914 setsigne(z[2], sx);
1915
1916 if (m == 2)
1917 { /* same formula as below, written more efficiently */
1918 y = gneg_i(y);
1919 if (typ(x) == t_REAL && signe(x) < 0)
1920 p1 = logr_abs(x);
1921 else
1922 p1 = gsub(glog(x,l), z);
1923 p1 = gmul2n(gsqr(p1), -1); /* = (log(-x))^2 / 2 */
1924
1925 p1 = gadd(p1, divrs(gsqr(mppi(l)), 6));
1926 p1 = gneg_i(p1);
1927 }
1928 else
1929 {
1930 GEN logx = glog(x,l), logx2 = gsqr(logx); p1 = gneg_i(ghalf);
1931 for (i=m-2; i>=0; i-=2)
1932 p1 = gadd(szeta(m-i,l), gmul(p1,gdivgs(logx2,(i+1)*(i+2))));
1933 if (m&1) p1 = gmul(logx,p1); else y = gneg_i(y);
1934 p1 = gadd(gmul2n(p1,1), gmul(z,gpowgs(logx,m-1)));
1935 if (typ(x) == t_REAL && signe(x) < 0) p1 = real_i(p1);
1936 }
1937 return gerepileupto(av, gadd(y,p1));
1938 }
1939
1940 GEN
dilog(GEN x,long prec)1941 dilog(GEN x, long prec)
1942 {
1943 return gpolylog(2, x, prec);
1944 }
1945
1946 GEN
polylogd0(long m,GEN x,long flag,long prec)1947 polylogd0(long m, GEN x, long flag, long prec)
1948 {
1949 long k, l, fl, m2;
1950 pari_sp av;
1951 GEN p1,p2,p3,y;
1952
1953 m2=m&1; av=avma;
1954 if (gcmp0(x)) return gcopy(x);
1955 if (gcmp1(x) && m>=2) return m2?szeta(m,prec):gen_0;
1956 l = precision(x);
1957 if (!l) { l=prec; x=gmul(x,real_1(l)); }
1958 p1 = gabs(x,prec); fl=0;
1959 if (expo(p1) >= 0) { x=ginv(x); p1=gabs(x,prec); fl=!m2; }
1960
1961 p1=gneg_i(glog(p1,prec)); p2=gen_1;
1962 y=polylog(m,x,prec); y = m2? real_i(y): imag_i(y);
1963 for (k=1; k<m; k++)
1964 {
1965 p2=gdivgs(gmul(p2,p1),k);
1966 p3=m2? real_i(polylog(m-k,x,prec)): imag_i(polylog(m-k,x,prec));
1967 y=gadd(y,gmul(p2,p3));
1968 }
1969 if (m2)
1970 {
1971 if (flag)
1972 p2 = gdivgs(gmul(p2,p1),-2*m);
1973 else
1974 p2 = gdivgs(gmul(glog(gnorm(gsub(gen_1,x)),prec),p2),2*m);
1975 y=gadd(y,p2);
1976 }
1977 if (fl) y = gneg(y);
1978 return gerepileupto(av, y);
1979 }
1980
1981 GEN
polylogd(long m,GEN x,long prec)1982 polylogd(long m, GEN x, long prec)
1983 {
1984 return polylogd0(m,x,0,prec);
1985 }
1986
1987 GEN
polylogdold(long m,GEN x,long prec)1988 polylogdold(long m, GEN x, long prec)
1989 {
1990 return polylogd0(m,x,1,prec);
1991 }
1992
1993 GEN
polylogp(long m,GEN x,long prec)1994 polylogp(long m, GEN x, long prec)
1995 {
1996 long k, l, fl, m2;
1997 pari_sp av;
1998 GEN p1,y;
1999
2000 m2=m&1; av=avma;
2001 if (gcmp0(x)) return gcopy(x);
2002 if (gcmp1(x) && m>=2) return m2?szeta(m,prec):gen_0;
2003 l=precision(x);
2004 if (!l) { l=prec; x=gmul(x,real_1(l)); }
2005 p1=gabs(x,prec); fl=0;
2006 if (expo(p1) >= 0) { x=ginv(x); p1=gabs(x,prec); fl=!m2; }
2007
2008 p1=gmul2n(glog(p1,prec),1); mpbern(m>>1,prec);
2009 y=polylog(m,x,prec); y=m2?real_i(y):imag_i(y);
2010
2011 if (m==1)
2012 {
2013 p1=gmul2n(p1,-2); y=gadd(y,p1);
2014 }
2015 else
2016 {
2017 GEN p2=gen_1, p3, p4, p5, p51=cgetr(prec);
2018
2019 for (k=1; k<m; k++)
2020 {
2021 p2=gdivgs(gmul(p2,p1),k);
2022 if (!(k&1) || k==1)
2023 {
2024 if (k!=1)
2025 {
2026 p5=(GEN)bern(k>>1);
2027 if (bernzone[2]>prec) { affrr(p5,p51); p5=p51; }
2028 p4=gmul(p2,p5);
2029 }
2030 else p4=gneg_i(gmul2n(p2,-1));
2031 p3=polylog(m-k,x,prec); p3=m2?real_i(p3):imag_i(p3);
2032 y=gadd(y,gmul(p4,p3));
2033 }
2034 }
2035 }
2036 if (fl) y = gneg(y);
2037 return gerepileupto(av, y);
2038 }
2039
2040 GEN
gpolylog(long m,GEN x,long prec)2041 gpolylog(long m, GEN x, long prec)
2042 {
2043 long i, lx, n, v;
2044 pari_sp av = avma;
2045 GEN a, y, p1;
2046
2047 if (m <= 0)
2048 {
2049 GEN t = mkpoln(2, gen_m1, gen_1); /* 1 - X */
2050 p1 = pol_x[0];
2051 for (i=2; i <= -m; i++)
2052 p1 = gmul(pol_x[0], gadd(gmul(t,derivpol(p1)), gmulsg(i,p1)));
2053 p1 = gdiv(p1, gpowgs(t,1-m));
2054 return gerepileupto(av, poleval(p1,x));
2055 }
2056
2057 switch(typ(x))
2058 {
2059 case t_INT: case t_REAL: case t_FRAC: case t_COMPLEX: case t_QUAD:
2060 return polylog(m,x,prec);
2061
2062 case t_POLMOD:
2063 p1=cleanroots(gel(x,1),prec); lx=lg(p1);
2064 for (i=1; i<lx; i++) gel(p1,i) = poleval(gel(x,2),gel(p1,i));
2065 y=cgetg(lx,t_COL);
2066 for (i=1; i<lx; i++) gel(y,i) = polylog(m,gel(p1,i),prec);
2067 return gerepileupto(av, y);
2068
2069 case t_INTMOD: case t_PADIC: pari_err(impl, "padic polylogarithm");
2070 default:
2071 av = avma; if (!(y = toser_i(x))) break;
2072 if (!m) { avma = av; return gneg(ghalf); }
2073 if (m==1) return gerepileupto(av, gneg( glog(gsub(gen_1,y),prec) ));
2074 if (gcmp0(y)) return gcopy(y);
2075 v = valp(y);
2076 if (v <= 0) pari_err(impl,"polylog around a!=0");
2077 n = (lg(y)-3 + v) / v;
2078 a = zeroser(varn(y), lg(y)-2);
2079 for (i=n; i>=1; i--)
2080 a = gmul(y, gadd(a, gpowgs(utoipos(i),-m)));
2081 return gerepileupto(av, a);
2082
2083 case t_VEC: case t_COL: case t_MAT:
2084 lx=lg(x); y=cgetg(lx,typ(x));
2085 for (i=1; i<lx; i++)
2086 gel(y,i) = gpolylog(m,gel(x,i),prec);
2087 return y;
2088 }
2089 pari_err(typeer,"gpolylog");
2090 return NULL; /* not reached */
2091 }
2092
2093 void
gpolylogz(long m,GEN x,GEN y)2094 gpolylogz(long m, GEN x, GEN y)
2095 {
2096 long prec = precision(y);
2097 pari_sp av=avma;
2098
2099 if (!prec) pari_err(infprecer,"gpolylogz");
2100 gaffect(gpolylog(m,x,prec),y); avma=av;
2101 }
2102
2103 GEN
polylog0(long m,GEN x,long flag,long prec)2104 polylog0(long m, GEN x, long flag, long prec)
2105 {
2106 switch(flag)
2107 {
2108 case 0: return gpolylog(m,x,prec);
2109 case 1: return polylogd(m,x,prec);
2110 case 2: return polylogdold(m,x,prec);
2111 case 3: return polylogp(m,x,prec);
2112 default: pari_err(flagerr,"polylog");
2113 }
2114 return NULL; /* not reached */
2115 }
2116
2117 static GEN
upper_half(GEN x,long * prec)2118 upper_half(GEN x, long *prec)
2119 {
2120 long tx = typ(x), l;
2121 if (tx == t_QUAD) { x = quadtoc(x, *prec); tx = typ(x); }
2122 if (tx != t_COMPLEX || gsigne(gel(x,2)) <= 0)
2123 pari_err(talker,"argument must belong to upper half-plane");
2124 l = precision(x); if (l) *prec = l;
2125 return x;
2126 }
2127
2128 static GEN
qq(GEN x,long prec)2129 qq(GEN x, long prec)
2130 {
2131 long tx = typ(x);
2132
2133 if (is_scalar_t(tx))
2134 {
2135 if (tx == t_PADIC) return x;
2136 x = upper_half(x, &prec);
2137 return gexp(gmul(mulcxI(x), Pi2n(1,prec)), prec); /* e(x) */
2138 }
2139 if (! ( x = toser_i(x)) ) pari_err(talker,"bad argument for modular function");
2140 return x;
2141 }
2142
2143 static GEN
inteta(GEN q)2144 inteta(GEN q)
2145 {
2146 long tx=typ(q);
2147 GEN p1,ps,qn,y;
2148
2149 y=gen_1; qn=gen_1; ps=gen_1;
2150 if (tx==t_PADIC)
2151 {
2152 if (valp(q) <= 0) pari_err(talker,"non-positive valuation in eta");
2153 for(;;)
2154 {
2155 p1 = gneg_i(gmul(ps,gmul(q,gsqr(qn))));
2156 y = gadd(y,p1); qn = gmul(qn,q); ps = gmul(p1,qn);
2157 p1 = y;
2158 y = gadd(y,ps); if (gequal(p1,y)) return y;
2159 }
2160 }
2161 else
2162 {
2163 long l, v = 0; /* gcc -Wall */
2164 pari_sp av = avma, lim = stack_lim(av, 3);
2165
2166 if (is_scalar_t(tx)) l = -bit_accuracy(precision(q));
2167 else
2168 {
2169 v = gvar(q); l = lg(q)-2; tx = 0;
2170 if (valp(q) <= 0) pari_err(talker,"non-positive valuation in eta");
2171 }
2172 for(;;)
2173 {
2174 p1 = gneg_i(gmul(ps,gmul(q,gsqr(qn))));
2175 /* qn = q^n
2176 * ps = (-1)^n q^(n(3n+1)/2)
2177 * p1 = (-1)^(n+1) q^(n(3n+1)/2 + 2n+1) */
2178 y = gadd(y,p1); qn = gmul(qn,q); ps = gmul(p1,qn);
2179 y = gadd(y,ps);
2180 if (tx)
2181 { if (gexpo(ps)-gexpo(y) < l) return y; }
2182 else
2183 { if (gval(ps,v) >= l) return y; }
2184 if (low_stack(lim, stack_lim(av,3)))
2185 {
2186 if(DEBUGMEM>1) pari_warn(warnmem,"eta");
2187 gerepileall(av, 3, &y, &qn, &ps);
2188 }
2189 }
2190 }
2191 }
2192
2193 GEN
eta(GEN x,long prec)2194 eta(GEN x, long prec)
2195 {
2196 pari_sp av = avma;
2197 GEN q = qq(x,prec);
2198 return gerepileupto(av,inteta(q));
2199 }
2200
2201 /* sqrt(3)/2 */
2202 static GEN
sqrt32(long prec)2203 sqrt32(long prec) { GEN z = sqrtr(stor(3, prec)); setexpo(z, -1); return z; }
2204
2205 /* exp(i x), x = k pi/12, assume 0 <= k < 24 */
2206 static GEN
e12(ulong k,long prec)2207 e12(ulong k, long prec)
2208 {
2209 int s, sPi, sPiov2;
2210 GEN z, t;
2211 if (k >12) { s = 1; k = 24 - k; } else s = 0; /* x -> 2pi - x */
2212 if (k > 6) { sPi = 1; k = 12 - k; } else sPi = 0; /* x -> pi - x */
2213 if (k > 3) { sPiov2 = 1; k = 6 - k; } else sPiov2 = 0; /* x -> pi/2 - x */
2214 z = cgetg(3, t_COMPLEX);
2215 switch(k)
2216 {
2217 case 0: gel(z,1) = icopy(gen_1); gel(z,2) = gen_0; break;
2218 case 1: t = gmul2n(addrs(sqrt32(prec), 1), -1);
2219 gel(z,1) = sqrtr(t);
2220 gel(z,2) = gmul2n(ginv(gel(z,1)), -2); break;
2221
2222 case 2: gel(z,1) = sqrt32(prec);
2223 gel(z,2) = real2n(-1, prec); break;
2224
2225 case 3: gel(z,1) = ginv( gsqrt(gen_2, prec) );
2226 gel(z,2) = mpcopy(gel(z,1)); break;
2227 }
2228 if (sPiov2) lswap(z[1], z[2]);
2229 if (sPi) setsigne(z[1], -signe(z[1]));
2230 if (s) setsigne(z[2], -signe(z[2]));
2231 return z;
2232 }
2233
2234 /* returns the true value of eta(x) for Im(x) > 0, using reduction to
2235 * standard fundamental domain */
2236 GEN
trueeta(GEN x,long prec)2237 trueeta(GEN x, long prec)
2238 {
2239 long tx = typ(x);
2240 ulong Nmod24;
2241 pari_sp av = avma;
2242 GEN q24, N, n, m, run;
2243
2244 if (!is_scalar_t(tx)) pari_err(typeer,"trueeta");
2245 x = upper_half(x, &prec);
2246 run = dbltor(1 - 1e-8);
2247 m = gen_1;
2248 N = gen_0;
2249 for(;;)
2250 {
2251 n = ground( real_i(x) );
2252 if (signe(n)) { x = gsub(x,n); N = addii(N, n); }
2253 if (gcmp(cxnorm(x), run) > 0) break;
2254 x = gdivsg(-1,x);
2255 m = gmul(m, gsqrt(mulcxmI(x), prec));
2256 }
2257 Nmod24 = umodiu(N, 24);
2258 if (Nmod24) m = gmul(m, e12(Nmod24, prec));
2259 q24 = gexp(gdivgs(gmul(Pi2n(1,prec), mulcxI(x)), 24),prec); /* e(x/24) */
2260 m = gmul(q24, m);
2261 if (24 * gexpo(q24) >= -bit_accuracy(prec))
2262 m = gmul(m, inteta( gpowgs(q24,24) ));
2263 return gerepileupto(av, m);
2264 }
2265
2266 GEN
eta0(GEN x,long flag,long prec)2267 eta0(GEN x, long flag,long prec)
2268 {
2269 return flag? trueeta(x,prec): eta(x,prec);
2270 }
2271
2272 GEN
jell(GEN x,long prec)2273 jell(GEN x, long prec)
2274 {
2275 long tx = typ(x);
2276 pari_sp av = avma;
2277 GEN p1;
2278
2279 if (!is_scalar_t(tx) || tx == t_PADIC)
2280 {
2281 GEN p2, q = qq(x,prec);
2282 p1 = gdiv(inteta(gsqr(q)), inteta(q));
2283 p1 = gmul2n(gsqr(p1),1);
2284 p1 = gmul(q,gpowgs(p1,12));
2285 p2 = gaddsg(768,gadd(gsqr(p1),gdivsg(4096,p1)));
2286 p1 = gmulsg(48,p1);
2287 return gerepileupto(av, gadd(p2,p1));
2288 }
2289 p1 = gdiv(trueeta(gmul2n(x,1),prec), trueeta(x,prec));
2290 p1 = gsqr(gsqr(gsqr(p1)));
2291 p1 = gadd(gmul2n(gsqr(p1),8), ginv(p1));
2292 return gerepileupto(av, gpowgs(p1,3));
2293 }
2294
2295 GEN
weberf2(GEN x,long prec)2296 weberf2(GEN x, long prec)
2297 {
2298 pari_sp av=avma, tetpil;
2299 GEN p1,p2;
2300
2301 p1=gsqrt(gen_2,prec);
2302 p2=gdiv(trueeta(gmul2n(x,1),prec),trueeta(x,prec));
2303 tetpil=avma;
2304 return gerepile(av,tetpil,gmul(p1,p2));
2305 }
2306
2307 GEN
weberf1(GEN x,long prec)2308 weberf1(GEN x, long prec)
2309 {
2310 pari_sp av=avma, tetpil;
2311 GEN p1,p2;
2312
2313 p1=trueeta(gmul2n(x,-1),prec); p2=trueeta(x,prec);
2314 tetpil=avma;
2315 return gerepile(av,tetpil,gdiv(p1,p2));
2316 }
2317
2318 GEN
weberf(GEN x,long prec)2319 weberf(GEN x, long prec)
2320 {
2321 pari_sp av = avma, tetpil;
2322 GEN p1, p2;
2323
2324 p1 = gdiv(trueeta(gmul2n(gaddgs(x,1),-1),prec),trueeta(x,prec));
2325 p2 = exp_Ir(divrs(mppi(prec),-24));
2326 tetpil = avma;
2327 return gerepile(av,tetpil, gmul(p1,p2));
2328 }
2329
2330 GEN
weber0(GEN x,long flag,long prec)2331 weber0(GEN x, long flag,long prec)
2332 {
2333 switch(flag)
2334 {
2335 case 0: return weberf(x,prec);
2336 case 1: return weberf1(x,prec);
2337 case 2: return weberf2(x,prec);
2338 default: pari_err(flagerr,"weber");
2339 }
2340 return NULL; /* not reached */
2341 }
2342
2343 GEN
theta(GEN q,GEN z,long prec)2344 theta(GEN q, GEN z, long prec)
2345 {
2346 long l, n;
2347 pari_sp av = avma;
2348 GEN ps, qn, y, zy, ps2, k, zold;
2349
2350 l = precision(q);
2351 n = precision(z); if (n && n < l) l = n;
2352 if (l) prec = l;
2353 z = gtofp(z, prec);
2354 q = gtofp(q, prec); if (gexpo(q) >= 0) pari_err(talker,"q >= 1 in theta");
2355 zold = NULL; /* gcc -Wall */
2356 zy = imag_i(z);
2357 if (gcmp0(zy)) k = gen_0;
2358 else
2359 {
2360 GEN lq = glog(q,prec); k = roundr(divrr(zy, real_i(lq)));
2361 if (signe(k)) { zold = z; z = gadd(z, mulcxmI(gmul(lq,k))); }
2362 }
2363 qn = gen_1;
2364 ps2 = gsqr(q);
2365 ps = gneg_i(ps2);
2366 y = gsin(z,prec);
2367 for (n = 1;; n++)
2368 {
2369 GEN t;
2370 qn = gmul(qn,ps);
2371 ps = gmul(ps,ps2);
2372 t = gmul(qn, gsin(gmulsg(2*n+1,z),prec)); y = gadd(y, t);
2373 if (gexpo(t) < -bit_accuracy(prec)) break;
2374 }
2375 if (signe(k))
2376 {
2377 y = gmul(y, gmul(powgi(q,sqri(k)),
2378 gexp(gmul(mulcxI(zold),shifti(k,1)), prec)));
2379 if (mod2(k)) y = gneg_i(y);
2380 }
2381 return gerepileupto(av, gmul(y, gmul2n(gsqrt(gsqrt(q,prec),prec),1)));
2382 }
2383
2384 GEN
thetanullk(GEN q,long k,long prec)2385 thetanullk(GEN q, long k, long prec)
2386 {
2387 long l, n;
2388 pari_sp av = avma;
2389 GEN p1, ps, qn, y, ps2;
2390
2391 if (k < 0) pari_err(talker,"k < 0 in thetanullk");
2392 l = precision(q);
2393 if (l) prec = l;
2394 q = gtofp(q, prec); if (gexpo(q) >= 0) pari_err(talker,"q >= 1 in theta");
2395
2396 if (!(k&1)) { avma = av; return gen_0; }
2397 qn = gen_1;
2398 ps2 = gsqr(q);
2399 ps = gneg_i(ps2);
2400 y = gen_1;
2401 for (n = 1;; n++)
2402 {
2403 GEN t;
2404 qn = gmul(qn,ps);
2405 ps = gmul(ps,ps2);
2406 t = gmul(qn, powuu(2*n+1, k)); y = gadd(y, t);
2407 if (gexpo(t) < -bit_accuracy(prec)) break;
2408 }
2409 p1 = gmul2n(gsqrt(gsqrt(q,prec),prec),1);
2410 if (k&2) y = gneg_i(y);
2411 return gerepileupto(av, gmul(p1, y));
2412 }
2413
2414 /* [d^i theta/dz^i(q, 0), i = 1, 3, .., 2*k - 1] */
2415 GEN
vecthetanullk(GEN q,long k,long prec)2416 vecthetanullk(GEN q, long k, long prec)
2417 {
2418 long i, l, n;
2419 pari_sp av = avma;
2420 GEN p1, ps, qn, y, ps2;
2421
2422 l = precision(q); if (l) prec = l;
2423 q = gtofp(q, prec); if (gexpo(q) >= 0) pari_err(talker,"q >= 1 in theta");
2424
2425 qn = gen_1;
2426 ps2 = gsqr(q);
2427 ps = gneg_i(ps2);
2428 y = cgetg(k+1, t_VEC); for (i = 1; i <= k; i++) gel(y,i) = gen_1;
2429 for (n = 1;; n++)
2430 {
2431 ulong N = 2*n + 1;
2432 GEN t = NULL/*-Wall*/, P = utoipos(N), N2 = sqru(N);
2433 qn = gmul(qn,ps);
2434 ps = gmul(ps,ps2);
2435 for (i = 1; i <= k; i++)
2436 {
2437 t = gmul(qn, P); gel(y,i) = gadd(gel(y,i), t);
2438 P = mulii(P, N2);
2439 }
2440 if (gexpo(t) < -bit_accuracy(prec)) break;
2441 }
2442 p1 = gmul2n(gsqrt(gsqrt(q,prec),prec),1);
2443 for (i = 2; i <= k; i+=2) gel(y,i) = gneg_i(gel(y,i));
2444 return gerepileupto(av, gmul(p1, y));
2445 }
2446