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 : /** TRANSCENDENTAL FUNCTIONS **/
18 : /** (part 2) **/
19 : /** **/
20 : /********************************************************************/
21 : #include "pari.h"
22 : #include "paripriv.h"
23 :
24 : #define DEBUGLEVEL DEBUGLEVEL_trans
25 :
26 : GEN
27 251016 : trans_fix_arg(long *prec, GEN *s0, GEN *sig, GEN *tau, pari_sp *av, GEN *res)
28 : {
29 251016 : GEN p1, s = *s0 = cxtoreal(*s0);
30 : long l;
31 251016 : l = precision(s); if (!l) l = *prec;
32 251015 : if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
33 251015 : *res = cgetc(l); *av = avma;
34 251015 : if (typ(s) == t_COMPLEX)
35 : { /* s = sig + i t */
36 206787 : s = cxtofp(s, l+EXTRAPREC64);
37 206787 : *sig = gel(s,1);
38 206787 : *tau = gel(s,2);
39 : }
40 : else /* real number */
41 : {
42 44228 : *sig = s = gtofp(s, l+EXTRAPREC64);
43 44227 : *tau = gen_0;
44 44227 : p1 = trunc2nr(s, 0);
45 44224 : if (!signe(subri(s,p1))) *s0 = p1;
46 : }
47 251013 : *prec = l; return s;
48 : }
49 :
50 : /********************************************************************/
51 : /** **/
52 : /** ARCTANGENT **/
53 : /** **/
54 : /********************************************************************/
55 : /* atan(b/a), real a and b, suitable for gc_upto */
56 : static GEN
57 184 : atan2_agm(GEN a, GEN b, long prec)
58 184 : { return gel(logagmcx(mkcomplex(a, b), prec), 2); }
59 : GEN
60 4495733 : mpatan(GEN x)
61 : {
62 4495733 : long l, l1, l2, n, m, i, lp, e, s, sx = signe(x);
63 : pari_sp av0, av;
64 : double alpha, beta, delta;
65 : GEN y, p1, p2, p3, p4, p5, unr;
66 : int inv;
67 :
68 4495733 : if (!sx) return real_0_bit(expo(x));
69 4495698 : l = lp = realprec(x);
70 4495698 : if (absrnz_equal1(x)) { /* |x| = 1 */
71 18257 : y = Pi2n(-2, l+EXTRAPREC64); if (sx < 0) setsigne(y,-1);
72 18258 : return y;
73 : }
74 4477432 : if (l > AGM_ATAN_LIMIT)
75 172 : { av = avma; return gc_leaf(av, atan2_agm(gen_1, x, l)); }
76 :
77 4477260 : e = expo(x); inv = (e >= 0); /* = (|x| > 1 ) */
78 4477260 : if (e > 0) lp += nbits2extraprec(e);
79 :
80 4477260 : y = cgetr(lp); av0 = avma;
81 4477255 : p1 = rtor(x, l+EXTRAPREC64); setabssign(p1); /* p1 = |x| */
82 4477295 : if (inv) p1 = invr(p1);
83 4477288 : e = expo(p1);
84 4477288 : if (e < -100)
85 42965 : alpha = 1.65149612947 - e; /* log_2(Pi) - e */
86 : else
87 4434323 : alpha = log2(M_PI / atan(rtodbl(p1)));
88 4477287 : beta = (double)(prec2nbits(l)>>1);
89 4477307 : delta = 1 + beta - alpha/2;
90 4477307 : if (delta <= 0) { n = 1; m = 0; }
91 : else
92 : {
93 4442422 : double fi = alpha-2;
94 4442422 : if (delta >= fi*fi)
95 : {
96 4204419 : double t = 1 + sqrt(delta);
97 4204419 : n = (long)t;
98 4204419 : m = (long)(t - fi);
99 : }
100 : else
101 : {
102 238003 : n = (long)(1+beta/fi);
103 238003 : m = 0;
104 : }
105 : }
106 4477307 : l2 = l + nbits2extraprec(m);
107 4477304 : p2 = rtor(p1, l2); av = avma;
108 49594406 : for (i=1; i<=m; i++)
109 : {
110 45117143 : p5 = addsr(1, sqrr(p2)); setprec(p5,l2);
111 45115353 : p5 = addsr(1, sqrtr_abs(p5)); setprec(p5,l2);
112 45116891 : affrr(divrr(p2,p5), p2); set_avma(av);
113 : }
114 4477263 : p3 = sqrr(p2); l1 = minss(LOWDEFAULTPREC+EXTRAPREC64, l2); /* l1 increases to l2 */;
115 4477308 : unr = real_1(l2); setprec(unr,l1);
116 4477301 : p4 = cgetr(l2); setprec(p4,l1);
117 4477302 : affrr(divru(unr,2*n+1), p4);
118 4477281 : s = 0; e = expo(p3); av = avma;
119 55413561 : for (i = n; i > 1; i--) /* n >= 1. i = 1 done outside for efficiency */
120 : {
121 50936298 : setprec(p3,l1); p5 = mulrr(p4,p3);
122 50941394 : l1 += nbits2extraprec(dvmdsBIL(s - e, &s)<<TWOPOTBITS_IN_LONG);
123 50939970 : if (l1 > l2) l1 = l2;
124 50939970 : setprec(unr,l1); p5 = subrr(divru(unr,2*i-1), p5);
125 50936427 : setprec(p4,l1); affrr(p5,p4); set_avma(av);
126 : }
127 4477263 : setprec(p3, l2); p5 = mulrr(p4,p3); /* i = 1 */
128 4477307 : setprec(unr,l2); p4 = subrr(unr, p5);
129 :
130 4477270 : p4 = mulrr(p2,p4); shiftr_inplace(p4, m);
131 4477298 : if (inv) p4 = subrr(Pi2n(-1, lp), p4);
132 4477307 : if (sx < 0) togglesign(p4);
133 4477303 : affrr_fixlg(p4,y); set_avma(av0); return y;
134 : }
135 :
136 : GEN
137 59887 : gatan(GEN x, long prec)
138 : {
139 : pari_sp av;
140 : GEN a, y;
141 :
142 59887 : switch(typ(x))
143 : {
144 22640 : case t_REAL: return mpatan(x);
145 36141 : case t_COMPLEX: /* atan(x) = -i atanh(ix) */
146 36141 : if (ismpzero(gel(x,2))) return gatan(gel(x,1), prec);
147 33128 : av = avma; return gc_GEN(av, mulcxmI(gatanh(mulcxI(x),prec)));
148 1106 : default:
149 1106 : av = avma; if (!(y = toser_i(x))) break;
150 28 : if (valser(y) < 0) pari_err_DOMAIN("atan","valuation", "<", gen_0, x);
151 21 : if (lg(y)==2) return gc_GEN(av, y);
152 : /* lg(y) > 2 */
153 14 : a = integser(gdiv(derivser(y), gaddsg(1,gsqr(y))));
154 14 : if (!valser(y)) a = gadd(a, gatan(gel(y,2),prec));
155 14 : return gc_upto(av, a);
156 : }
157 1078 : return trans_eval("atan",gatan,x,prec);
158 : }
159 : /********************************************************************/
160 : /** **/
161 : /** ARCSINE **/
162 : /** **/
163 : /********************************************************************/
164 : /* |x| <= 1 */
165 : GEN
166 710 : mpasin(GEN x)
167 : {
168 710 : pari_sp av = avma;
169 710 : long sx = signe(x);
170 : GEN z, a;
171 :
172 710 : if (!sx) return rcopy(x);
173 710 : if (absrnz_equal1(x)) { /* |x| = 1 */
174 0 : if (sx > 0) return Pi2n(-1, realprec(x)); /* 1 */
175 0 : z = Pi2n(-1, realprec(x)); setsigne(z, -1); return z; /* -1 */
176 : }
177 710 : a = sqrtr(subsr(1, sqrr(x)));
178 710 : if (realprec(x) > AGM_ATAN_LIMIT)
179 4 : z = atan2_agm(a, x, realprec(x));
180 : else
181 706 : z = mpatan(divrr(x, a));
182 710 : return gc_leaf(av, z);
183 : }
184 :
185 : GEN
186 9005 : gasin(GEN x, long prec)
187 : {
188 : long sx;
189 : pari_sp av;
190 : GEN a, y, p1;
191 :
192 9005 : switch(typ(x))
193 : {
194 1095 : case t_REAL: sx = signe(x);
195 1095 : if (!sx) return rcopy(x);
196 1088 : if (absrnz_equal1(x)) { /* |x| = 1 */
197 28 : if (sx > 0) return Pi2n(-1, realprec(x)); /* 1 */
198 14 : y = Pi2n(-1, realprec(x)); setsigne(y, -1); return y; /* -1 */
199 : }
200 1060 : if (expo(x) < 0) return mpasin(x);
201 350 : y = cgetg(3,t_COMPLEX);
202 350 : gel(y,1) = Pi2n(-1, realprec(x));
203 350 : gel(y,2) = mpacosh(x);
204 350 : if (sx < 0) togglesign(gel(y,1)); else togglesign(gel(y,2));
205 350 : return y;
206 :
207 7406 : case t_COMPLEX: /* asin(z) = -i asinh(iz) */
208 7406 : if (ismpzero(gel(x,2))) return gasin(gel(x,1), prec);
209 7406 : av = avma;
210 7406 : return gc_GEN(av, mulcxmI(gasinh(mulcxI(x), prec)));
211 504 : default:
212 504 : av = avma; if (!(y = toser_i(x))) break;
213 42 : if (gequal0(y)) return gc_GEN(av, y);
214 : /* lg(y) > 2*/
215 35 : if (valser(y) < 0) pari_err_DOMAIN("asin","valuation", "<", gen_0, x);
216 28 : p1 = gsubsg(1,gsqr(y));
217 28 : if (gequal0(p1))
218 : {
219 21 : GEN t = Pi2n(-1,prec);
220 21 : if (gsigne(gel(y,2)) < 0) setsigne(t, -1);
221 21 : return gc_upto(av, scalarser(t, varn(y), valser(p1)>>1));
222 : }
223 7 : p1 = gdiv(derivser(y), gsqrt(p1,prec));
224 7 : a = integser(p1);
225 7 : if (!valser(y)) a = gadd(a, gasin(gel(y,2),prec));
226 7 : return gc_upto(av, a);
227 : }
228 462 : return trans_eval("asin",gasin,x,prec);
229 : }
230 : /********************************************************************/
231 : /** **/
232 : /** ARCCOSINE **/
233 : /** **/
234 : /********************************************************************/
235 : static GEN
236 14 : acos0(long e) { return Pi2n(-1, nbits2prec(e<0? -e: 1)); }
237 :
238 : /* |x| <= 1 */
239 : GEN
240 105 : mpacos(GEN x)
241 : {
242 105 : pari_sp av = avma;
243 105 : long sx = signe(x);
244 : GEN z, a;
245 :
246 105 : if (!sx) return acos0(expo(x));
247 105 : if (absrnz_equal1(x)) /* |x| = 1 */
248 0 : return sx > 0? real_0_bit( -(bit_prec(x)>>1) ) : mppi(realprec(x));
249 105 : a = sqrtr(subsr(1, sqrr(x)));
250 105 : if (realprec(x) > AGM_ATAN_LIMIT)
251 8 : z = atan2_agm(x, a, realprec(x));
252 : else
253 : {
254 97 : z = mpatan(divrr(a, x));
255 97 : if (signe(x) < 0) z = addrr(mppi(realprec(z)), z);
256 : }
257 105 : return gc_leaf(av, z);
258 : }
259 :
260 : GEN
261 7931 : gacos(GEN x, long prec)
262 : {
263 : long sx;
264 : pari_sp av;
265 : GEN a, y, p1;
266 :
267 7931 : switch(typ(x))
268 : {
269 252 : case t_REAL: sx = signe(x);
270 252 : if (!sx) return acos0(expo(x));
271 245 : if (absrnz_equal1(x)) /* |x| = 1 */
272 14 : return sx > 0? real_0_bit( -(bit_prec(x)>>1) ) : mppi(realprec(x));
273 231 : if (expo(x) < 0) return mpacos(x);
274 :
275 175 : y = cgetg(3,t_COMPLEX); p1 = mpacosh(x);
276 175 : if (sx < 0) { gel(y,1) = mppi(realprec(x)); togglesign(p1); }
277 91 : else gel(y,1) = gen_0;
278 175 : gel(y,2) = p1; return y;
279 :
280 7406 : case t_COMPLEX:
281 7406 : if (ismpzero(gel(x,2))) return gacos(gel(x,1), prec);
282 7406 : av = avma;
283 7406 : p1 = gadd(x, mulcxI(gsqrt(gsubsg(1,gsqr(x)), prec)));
284 7406 : y = glog(p1,prec); /* log(x + I*sqrt(1-x^2)) */
285 7406 : return gc_GEN(av, mulcxmI(y));
286 273 : default:
287 273 : av = avma; if (!(y = toser_i(x))) break;
288 35 : if (valser(y) < 0) pari_err_DOMAIN("acos","valuation", "<", gen_0, x);
289 28 : if (lg(y) > 2)
290 : {
291 21 : p1 = gsubsg(1,gsqr(y));
292 21 : if (gequal0(p1)) { set_avma(av); return zeroser(varn(y), valser(p1)>>1); }
293 7 : p1 = integser(gdiv(gneg(derivser(y)), gsqrt(p1,prec)));
294 : /*y(t) = 1+O(t)*/
295 7 : if (gequal1(gel(y,2)) && !valser(y)) return gc_upto(av, p1);
296 : }
297 7 : else p1 = y;
298 14 : a = (lg(y)==2 || valser(y))? Pi2n(-1, prec): gacos(gel(y,2),prec);
299 14 : return gc_upto(av, gadd(a,p1));
300 : }
301 238 : return trans_eval("acos",gacos,x,prec);
302 : }
303 : /********************************************************************/
304 : /** **/
305 : /** ARGUMENT **/
306 : /** **/
307 : /********************************************************************/
308 :
309 : /* we know that x and y are not both 0 */
310 : static GEN
311 4472328 : mparg(GEN x, GEN y)
312 : {
313 4472328 : long prec, sx = signe(x), sy = signe(y);
314 : GEN z;
315 :
316 4472328 : if (!sy)
317 : {
318 0 : if (sx > 0) return real_0_bit(expo(y) - expo(x));
319 0 : return mppi(realprec(x));
320 : }
321 4472328 : prec = realprec(y); if (prec < realprec(x)) prec = realprec(x);
322 4472328 : if (!sx)
323 : {
324 27 : z = Pi2n(-1, prec); if (sy < 0) setsigne(z,-1);
325 27 : return z;
326 : }
327 :
328 4472301 : if (expo(x)-expo(y) > -2)
329 : {
330 3525173 : z = mpatan(divrr(y,x)); if (sx > 0) return z;
331 1477087 : return addrr_sign(z, signe(z), mppi(prec), sy);
332 : }
333 947128 : z = mpatan(divrr(x,y));
334 947145 : return addrr_sign(z, -signe(z), Pi2n(-1, prec), sy);
335 : }
336 :
337 : static GEN
338 8944660 : rfix(GEN x,long prec)
339 : {
340 8944660 : switch(typ(x))
341 : {
342 38027 : case t_INT: return itor(x, prec);
343 641946 : case t_FRAC: return fractor(x, prec);
344 8264687 : case t_REAL: break;
345 0 : default: pari_err_TYPE("rfix (conversion to t_REAL)",x);
346 : }
347 8264687 : return x;
348 : }
349 :
350 : static GEN
351 4472338 : cxarg(GEN x, GEN y, long prec)
352 : {
353 4472338 : pari_sp av = avma;
354 4472338 : x = rfix(x,prec);
355 4472336 : y = rfix(y,prec); return gc_leaf(av, mparg(x,y));
356 : }
357 :
358 : GEN
359 4489868 : garg(GEN x, long prec)
360 : {
361 : long l;
362 4489868 : if (gequal0(x)) pari_err_DOMAIN("arg", "argument", "=", gen_0, x);
363 4489867 : switch(typ(x))
364 : {
365 17528 : case t_REAL: prec = realprec(x); /* fall through */
366 17528 : case t_INT: case t_FRAC: return (gsigne(x)>0)? real_0(prec): mppi(prec);
367 4472339 : case t_COMPLEX:
368 4472339 : l = precision(x); if (l) prec = l;
369 4472338 : return cxarg(gel(x,1),gel(x,2),prec);
370 : }
371 0 : return trans_eval("arg",garg,x,prec);
372 : }
373 :
374 : /********************************************************************/
375 : /** **/
376 : /** HYPERBOLIC COSINE **/
377 : /** **/
378 : /********************************************************************/
379 : /* 1 + x */
380 : static GEN
381 7 : mpcosh0(long e) { return e >= 0? real_0_bit(e): real_1_bit(-e); }
382 : GEN
383 3528 : mpcosh(GEN x)
384 : {
385 : pari_sp av;
386 : GEN z;
387 :
388 3528 : if (!signe(x)) return mpcosh0(expo(x));
389 3521 : av = avma;
390 3521 : z = mpexp(x); z = addrr(z, invr(z)); shiftr_inplace(z, -1);
391 3521 : return gc_leaf(av, z);
392 : }
393 :
394 : GEN
395 3619 : gcosh(GEN x, long prec)
396 : {
397 : pari_sp av;
398 : GEN y, p1;
399 : long v;
400 :
401 3619 : switch(typ(x))
402 : {
403 3528 : case t_REAL: return mpcosh(x);
404 21 : case t_COMPLEX:
405 21 : if (isintzero(gel(x,1))) return gcos(gel(x,2),prec);
406 : /* fall through */
407 : case t_PADIC:
408 21 : av = avma; p1 = gexp(x,prec); p1 = gadd(p1, ginv(p1));
409 21 : return gc_upto(av, gmul2n(p1,-1));
410 56 : default:
411 56 : av = avma; if (!(y = toser_i(x))) break;
412 35 : if (gequal0(y) && valser(y) == 0) return gc_GEN(av, y);
413 35 : v = valser(y);
414 35 : if (v > 0) y = sertoser(y, lg(y) - 2 + v);
415 35 : p1 = gexp(y,prec); p1 = gadd(p1, ginv(p1));
416 35 : return gc_upto(av, gmul2n(p1,-1));
417 : }
418 21 : return trans_eval("cosh",gcosh,x,prec);
419 : }
420 : /********************************************************************/
421 : /** **/
422 : /** HYPERBOLIC SINE **/
423 : /** **/
424 : /********************************************************************/
425 : static GEN
426 0 : mpsinh0(long e) { return real_0_bit(e); }
427 : GEN
428 6377 : mpsinh(GEN x)
429 : {
430 : pari_sp av;
431 : long lx;
432 : GEN z, res;
433 :
434 6377 : if (!signe(x)) return mpsinh0(expo(x));
435 6377 : lx = realprec(x); res = cgetr(lx); av = avma;
436 6377 : if (expo(x) + BITS_IN_LONG < 1)
437 : { /* y = e^x-1; e^x - e^(-x) = y(1 + 1/(y+1)) */
438 7 : GEN y = mpexpm1(x);
439 7 : lx += EXTRAPRECWORD;
440 7 : z = addrs(y, 1); if (realprec(z) > lx) z = rtor(z,lx); /* e^x */
441 7 : z = mulrr(y, addsr(1, invr(z)));
442 : }
443 : else
444 : {
445 6370 : z = mpexp(x);
446 6370 : z = subrr(z, invr(z));
447 : }
448 6377 : shiftr_inplace(z, -1);
449 6377 : affrr(z, res); set_avma(av); return res;
450 : }
451 :
452 : void
453 411002 : mpsinhcosh(GEN x, GEN *s, GEN *c)
454 : {
455 : pari_sp av;
456 : long lx, ex;
457 : GEN z, zi, S, C;
458 411002 : if (!signe(x))
459 : {
460 0 : ex = expo(x);
461 0 : *c = mpcosh0(ex);
462 0 : *s = mpsinh0(ex); return;
463 : }
464 411002 : lx = realprec(x);
465 411002 : *c = cgetr(lx);
466 411002 : *s = cgetr(lx); av = avma;
467 411002 : if (expo(x) + BITS_IN_LONG < 1)
468 : { /* y = e^x-1; e^x - e^(-x) = y(1 + 1/(y+1)) */
469 403 : GEN y = mpexpm1(x);
470 403 : lx += EXTRAPRECWORD;
471 403 : z = addrs(y,1); if (realprec(z) > lx) z = rtor(z, lx); /* e^x */
472 403 : zi = invr(z); /* z = exp(x), zi = exp(-x) */
473 403 : S = mulrr(y, addsr(1,zi));
474 : }
475 : else
476 : {
477 410599 : z = mpexp(x);
478 410599 : zi = invr(z);
479 410599 : S = subrr(z, zi);
480 : }
481 411002 : C = addrr(z, zi);
482 411002 : shiftr_inplace(S, -1); affrr(S, *s);
483 411002 : shiftr_inplace(C, -1); affrr(C, *c); set_avma(av);
484 : }
485 :
486 : GEN
487 12642 : gsinh(GEN x, long prec)
488 : {
489 : pari_sp av;
490 : GEN y, p1;
491 :
492 12642 : switch(typ(x))
493 : {
494 6027 : case t_REAL: return mpsinh(x);
495 21 : case t_COMPLEX:
496 21 : if (isintzero(gel(x,1))) retmkcomplex(gen_0, gsin(gel(x,2),prec));
497 : /* fall through */
498 : case t_PADIC:
499 14 : av = avma; p1 = gexp(x,prec); p1 = gsub(p1, ginv(p1));
500 14 : return gc_upto(av, gmul2n(p1,-1));
501 6587 : default:
502 6587 : av = avma; if (!(y = toser_i(x))) break;
503 6559 : if (gequal0(y) && valser(y) == 0) return gc_GEN(av, y);
504 6559 : p1 = gexp(y, prec); p1 = gsub(p1, ginv(p1));
505 6559 : return gc_upto(av, gmul2n(p1,-1));
506 : }
507 28 : return trans_eval("sinh",gsinh,x,prec);
508 : }
509 : /********************************************************************/
510 : /** **/
511 : /** HYPERBOLIC TANGENT **/
512 : /** **/
513 : /********************************************************************/
514 :
515 : GEN
516 77056 : mptanh(GEN x)
517 : {
518 77056 : long lx, s = signe(x);
519 : GEN y;
520 :
521 77056 : if (!s) return real_0_bit(expo(x));
522 77056 : lx = realprec(x);
523 77056 : if (abscmprr(x, utor(prec2nbits(lx), LOWDEFAULTPREC)) >= 0) {
524 24840 : y = real_1(lx);
525 : } else {
526 52216 : pari_sp av = avma;
527 52216 : long e = expo(x) + BITS_IN_LONG;
528 : GEN t;
529 52216 : if (e < 0) x = rtor(x, lx + nbits2extraprec(-e));
530 52216 : t = exp1r_abs(gmul2n(x,1)); /* exp(|2x|) - 1 */
531 52216 : y = gc_leaf(av, divrr(t, addsr(2,t)));
532 : }
533 77056 : if (s < 0) togglesign(y); /* tanh is odd */
534 77056 : return y;
535 : }
536 :
537 : GEN
538 77161 : gtanh(GEN x, long prec)
539 : {
540 : pari_sp av;
541 : GEN y, t;
542 :
543 77161 : switch(typ(x))
544 : {
545 77056 : case t_REAL: return mptanh(x);
546 35 : case t_COMPLEX:
547 35 : if (isintzero(gel(x,1))) retmkcomplex(gen_0, gtan(gel(x,2),prec));
548 : /* fall through */
549 : case t_PADIC:
550 28 : av = avma;
551 28 : t = gexp(gmul2n(x,1),prec);
552 28 : t = gdivsg(-2, gaddgs(t,1));
553 28 : return gc_upto(av, gaddsg(1,t));
554 63 : default:
555 63 : av = avma; if (!(y = toser_i(x))) break;
556 28 : if (gequal0(y)) return gc_GEN(av, y);
557 14 : t = gexp(gmul2n(y, 1),prec);
558 14 : t = gdivsg(-2, gaddgs(t,1));
559 14 : return gc_upto(av, gaddsg(1,t));
560 : }
561 35 : return trans_eval("tanh",gtanh,x,prec);
562 : }
563 :
564 : GEN
565 7 : mpcotanh(GEN x)
566 : {
567 7 : long lx, s = signe(x);
568 : GEN y;
569 :
570 7 : if (!s) pari_err_DOMAIN("cotan", "argument", "=", gen_0, x);
571 :
572 7 : lx = realprec(x);
573 7 : if (abscmprr(x, utor(prec2nbits(lx), LOWDEFAULTPREC)) >= 0) {
574 0 : y = real_1(lx);
575 : } else {
576 7 : pari_sp av = avma;
577 7 : long e = expo(x) + BITS_IN_LONG;
578 : GEN t;
579 7 : if (e < 0) x = rtor(x, lx + nbits2extraprec(-e));
580 7 : t = exp1r_abs(gmul2n(x,1)); /* exp(|2x|) - 1 */
581 7 : y = gc_leaf(av, divrr(addsr(2,t), t));
582 : }
583 7 : if (s < 0) togglesign(y); /* cotanh is odd */
584 7 : return y;
585 : }
586 :
587 : GEN
588 63 : gcotanh(GEN x, long prec)
589 : {
590 : pari_sp av;
591 : GEN y, t;
592 :
593 63 : switch(typ(x))
594 : {
595 7 : case t_REAL: return mpcotanh(x);
596 14 : case t_COMPLEX:
597 14 : if (isintzero(gel(x,1))) retmkcomplex(gen_0, gcotan(gel(x,2),prec));
598 : /* fall through */
599 : case t_PADIC:
600 14 : av = avma;
601 14 : t = gexpm1(gmul2n(x,1),prec);
602 14 : return gc_upto(av, gaddsg(1, gdivsg(2,t)));
603 35 : default:
604 35 : av = avma; if (!(y = toser_i(x))) break;
605 28 : if (gequal0(y)) return gc_GEN(av, y);
606 14 : t = gexpm1(gmul2n(y,1),prec);
607 14 : return gc_upto(av, gaddsg(1, gdivsg(2,t)));
608 : }
609 7 : return trans_eval("cotanh",gcotanh,x,prec);
610 : }
611 :
612 : /********************************************************************/
613 : /** **/
614 : /** AREA HYPERBOLIC SINE **/
615 : /** **/
616 : /********************************************************************/
617 :
618 : GEN
619 3346 : mpasinh(GEN x)
620 : {
621 3346 : long lx, e, s = signe(x);
622 : GEN z, res;
623 : pari_sp av;
624 :
625 3346 : if (!s) return rcopy(x);
626 3346 : lx = realprec(x); e = expo(x) + BITS_IN_LONG;
627 3346 : res = cgetr(lx);
628 3346 : av = avma;
629 3346 : if (e < 0) x = rtor(x, lx + nbits2extraprec(-e));
630 3346 : z = logr_abs( addrr_sign(x,1, sqrtr_abs( addrs(sqrr(x), 1) ), 1) );
631 3346 : if (signe(x) < 0) togglesign(z);
632 3346 : affrr(z, res); return gc_const(av, res);
633 : }
634 :
635 : GEN
636 39228 : gasinh(GEN x, long prec)
637 : {
638 : pari_sp av;
639 : GEN a, y, p1;
640 :
641 39228 : switch(typ(x))
642 : {
643 3031 : case t_REAL:
644 3031 : if (!signe(x)) return rcopy(x);
645 2996 : return mpasinh(x);
646 :
647 35560 : case t_COMPLEX: {
648 : GEN a, b, d;
649 35560 : if (ismpzero(gel(x,2))) return gasinh(gel(x,1), prec);
650 35560 : av = avma;
651 35560 : if (ismpzero(gel(x,1))) /* avoid cancellation */
652 843 : return gc_GEN(av, mulcxI(gasin(gel(x,2), prec)));
653 34717 : d = gsqrt(gaddsg(1,gsqr(x)), prec); /* Re(d) >= 0 */
654 34717 : a = gadd(d, x);
655 34717 : b = gsub(d, x);
656 : /* avoid cancellation as much as possible */
657 34717 : if (gprecision(a) < gprecision(b))
658 378 : y = gneg(glog(b,prec));
659 : else
660 34339 : y = glog(a,prec);
661 34717 : return gc_upto(av, y); /* log (x + sqrt(1+x^2)) */
662 : }
663 637 : default:
664 637 : av = avma; if (!(y = toser_i(x))) break;
665 168 : if (gequal0(y)) return gc_GEN(av, y);
666 161 : if (valser(y) < 0) pari_err_DOMAIN("asinh","valuation", "<", gen_0, x);
667 154 : p1 = gaddsg(1,gsqr(y));
668 154 : if (gequal0(p1))
669 : {
670 14 : GEN t = PiI2n(-1,prec);
671 14 : if ( gsigne(imag_i(gel(y,2))) < 0 ) setsigne(gel(t,2), -1);
672 14 : return gc_upto(av, scalarser(t, varn(y), valser(p1)>>1));
673 : }
674 140 : p1 = gdiv(derivser(y), gsqrt(p1,prec));
675 140 : a = integser(p1);
676 140 : if (!valser(y)) a = gadd(a, gasinh(gel(y,2),prec));
677 140 : return gc_upto(av, a);
678 : }
679 469 : return trans_eval("asinh",gasinh,x,prec);
680 : }
681 : /********************************************************************/
682 : /** **/
683 : /** AREA HYPERBOLIC COSINE **/
684 : /** **/
685 : /********************************************************************/
686 :
687 : /* |x| >= 1, return ach(|x|) */
688 : GEN
689 742 : mpacosh(GEN x)
690 : {
691 742 : long lx = realprec(x), e;
692 742 : GEN z, res = cgetr(lx);
693 742 : pari_sp av = avma;
694 742 : e = expo(signe(x) > 0? subrs(x,1): addrs(x,1));
695 742 : if (e == -(long)HIGHEXPOBIT)
696 0 : return gc_const((pari_sp)(res + lx), real_0_bit(- bit_prec(x) >> 1));
697 742 : if (e < -5) x = rtor(x, realprec(x) + nbits2extraprec(-e));
698 742 : z = logr_abs( addrr_sign(x, 1, sqrtr( subrs(sqrr(x), 1) ), 1) );
699 742 : affrr(z, res); return gc_const(av, res);
700 : }
701 :
702 : GEN
703 7994 : gacosh(GEN x, long prec)
704 : {
705 : pari_sp av;
706 : GEN y;
707 :
708 7994 : switch(typ(x))
709 : {
710 280 : case t_REAL: {
711 280 : long s = signe(x), e = expo(x);
712 : GEN a, b;
713 280 : if (s > 0 && e >= 0) return mpacosh(x);
714 : /* x < 1 */
715 147 : y = cgetg(3,t_COMPLEX); a = gen_0;
716 147 : if (s == 0) b = acos0(e);
717 140 : else if (e < 0) b = mpacos(x); /* -1 < x < 1 */
718 : else {
719 91 : if (!absrnz_equal1(x)) a = mpacosh(x);
720 91 : b = mppi(realprec(x));
721 : }
722 147 : gel(y,1) = a;
723 147 : gel(y,2) = b; return y;
724 : }
725 7413 : case t_COMPLEX: {
726 : GEN a, b, d;
727 7413 : if (ismpzero(gel(x,2))) return gacosh(gel(x,1), prec);
728 7413 : av = avma;
729 7413 : d = gsqrt(gaddsg(-1,gsqr(x)), prec); /* Re(d) >= 0 */
730 7413 : a = gadd(x, d);
731 7413 : b = gsub(x, d);
732 : /* avoid cancellation as much as possible */
733 7413 : if (gprecision(a) < gprecision(b))
734 7 : y = glog(b,prec);
735 : else
736 7406 : y = glog(a,prec);
737 : /* y = \pm log(x + sqrt(x^2-1)) */
738 7413 : if (gsigne(real_i(y)) < 0) y = gneg(y);
739 7413 : return gc_upto(av, y);
740 : }
741 301 : default: {
742 : GEN a, d;
743 : long v;
744 301 : av = avma; if (!(y = toser_i(x))) break;
745 49 : v = valser(y);
746 49 : if (v < 0) pari_err_DOMAIN("acosh","valuation", "<", gen_0, x);
747 42 : if (gequal0(y))
748 : {
749 7 : if (!v) return gc_GEN(av, y);
750 7 : return gc_upto(av, gadd(y, PiI2n(-1, prec)));
751 : }
752 35 : d = gsubgs(gsqr(y),1);
753 35 : if (gequal0(d)) { set_avma(av); return zeroser(varn(y), valser(d)>>1); }
754 21 : d = gdiv(derivser(y), gsqrt(d,prec));
755 21 : a = integser(d);
756 21 : if (v)
757 7 : d = PiI2n(-1, prec); /* I Pi/2 */
758 : else
759 : {
760 14 : d = gel(y,2); if (gequal1(d)) return gc_upto(av,a);
761 7 : d = gacosh(d, prec);
762 : }
763 14 : return gc_upto(av, gadd(d,a));
764 : }
765 : }
766 252 : return trans_eval("acosh",gacosh,x,prec);
767 : }
768 : /********************************************************************/
769 : /** **/
770 : /** AREA HYPERBOLIC TANGENT **/
771 : /** **/
772 : /********************************************************************/
773 :
774 : /* |x| < 1 */
775 : GEN
776 7665 : mpatanh(GEN x)
777 : {
778 7665 : pari_sp av = avma;
779 7665 : long e, s = signe(x);
780 : GEN z;
781 7665 : if (!s) return rcopy(x);
782 7665 : z = s > 0? subsr(1,x): addsr(1,x); e = expo(z);
783 7665 : if (e < -5)
784 : {
785 1218 : x = rtor(x, realprec(x) + nbits2extraprec(-e)-EXTRAPRECWORD);
786 1218 : z = s > 0? subsr(1,x): addsr(1,x); e = expo(z);
787 : }
788 7665 : z = invr(z); shiftr_inplace(z, 1); /* 2/(1-|x|) */
789 7665 : z = logr_abs( addrs(z,-1) ); if (s < 0) togglesign(z);
790 7665 : shiftr_inplace(z, -1); return gc_leaf(av, z);
791 : }
792 :
793 : static long
794 5505488 : get_nmax(double u, double v, long prec)
795 : {
796 5505488 : double d = 2 * log2(((double)v) / u); /* can be 0 due to rounding */
797 5505488 : long nmax = -1;
798 5505488 : if (d)
799 : {
800 5505522 : d = ceil(prec2nbits(prec) / d);
801 5505120 : if (dblexpo(d) < BITS_IN_LONG) nmax = (long)d;
802 : }
803 5504931 : return nmax;
804 : }
805 : /* atanh(u/v) using binary splitting, 0 < u < v */
806 : GEN
807 5516470 : atanhuu(ulong u, ulong v, long prec)
808 : {
809 5516470 : GEN u2 = sqru(u), v2 = sqru(v);
810 5505707 : long i, nmax = get_nmax((double)u, (double)v, prec);
811 : struct abpq_res R;
812 : struct abpq A;
813 5504776 : if (nmax < 0) pari_err_OVERFLOW("atanhuu");
814 5504769 : abpq_init(&A, nmax); /* nmax satisfies (2n+1) (v/u)^2n > 2^bitprec */
815 5515771 : A.a[0] = A.b[0] = gen_1;
816 5515771 : A.p[0] = utoipos(u);
817 5515831 : A.q[0] = utoipos(v);
818 106349501 : for (i = 1; i <= nmax; i++)
819 : {
820 100830416 : A.a[i] = gen_1;
821 100830416 : A.b[i] = utoipos((i<<1)+1);
822 100828540 : A.p[i] = u2;
823 100828540 : A.q[i] = v2;
824 : }
825 5519085 : abpq_sum(&R, 0, nmax, &A);
826 5516788 : return rdivii(R.T, mulii(R.B,R.Q),prec);
827 : }
828 : /* atanh(u/v) using binary splitting, 0 < u < v */
829 : GEN
830 14 : atanhui(ulong u, GEN v, long prec)
831 : {
832 14 : GEN u2 = sqru(u), v2 = sqri(v);
833 14 : long i, nmax = get_nmax((double)u, gtodouble(v), prec);
834 : struct abpq_res R;
835 : struct abpq A;
836 14 : if (nmax < 0) pari_err_OVERFLOW("atanhui");
837 14 : abpq_init(&A, nmax);
838 14 : A.a[0] = A.b[0] = gen_1;
839 14 : A.p[0] = utoipos(u);
840 14 : A.q[0] = v;
841 35 : for (i = 1; i <= nmax; i++)
842 : {
843 21 : A.a[i] = gen_1;
844 21 : A.b[i] = utoipos((i<<1)+1);
845 21 : A.p[i] = u2;
846 21 : A.q[i] = v2;
847 : }
848 14 : abpq_sum(&R, 0, nmax, &A);
849 14 : return rdivii(R.T, mulii(R.B,R.Q),prec);
850 : }
851 :
852 : static void
853 28 : err_atanh(GEN x, GEN bad) { pari_err_DOMAIN("atanh", "x", "=", bad, x); }
854 :
855 : GEN
856 68651 : gatanh(GEN x, long prec)
857 : {
858 : long sx;
859 : pari_sp av;
860 : GEN a, y, z;
861 :
862 68651 : switch(typ(x))
863 : {
864 126 : case t_INT:
865 126 : sx = signe(x);
866 126 : if (!sx) return real_0(prec);
867 119 : z = cgetg(3, t_COMPLEX); av = avma;
868 119 : if (lgefint(x) == 3)
869 : {
870 112 : if (x[2] == 1) err_atanh(x, sx == 1? gen_1: gen_m1);
871 84 : a = atanhuu(1, x[2], prec);
872 : }
873 : else
874 7 : a = atanhui(1, x, prec);
875 91 : gel(z,1) = gc_leaf(av, a);
876 91 : gel(z,2) = Pi2n(-1, prec);
877 91 : togglesign(sx > 0? gel(z,2): gel(z,1));
878 91 : return z;
879 350 : case t_FRAC:
880 : {
881 : long ly, lz, e;
882 :
883 350 : y = gel(x,1); ly = lgefint(y);
884 350 : z = gel(x,2); lz = lgefint(z); if (ly > 3 && lz > 3) break;
885 350 : if (abscmpii(y, z) > 0) /* |y| > z; lz = 3 */
886 : {
887 252 : ulong u = z[2];
888 252 : av = avma; e = expi((signe(y) < 0)? addii(y, z): subii(y, z));
889 252 : set_avma(av); if (e < - prec2nbits(prec)) break;
890 252 : z = cgetg(3, t_COMPLEX); av = avma;
891 252 : a = ly == 3? atanhuu(u, y[2], prec): atanhui(u, y, prec);
892 252 : gel(z,1) = gc_leaf(av, a);
893 252 : gel(z,2) = Pi2n(-1, prec);
894 252 : togglesign(signe(y) > 0? gel(z,2): gel(z,1));
895 : }
896 : else
897 : { /* |y| < z; ly = 3 */
898 98 : av = avma; e = expi((signe(y) < 0)? addii(y, z): subii(y, z));
899 98 : set_avma(av); if (e < - prec2nbits(prec)) break;
900 98 : a = lz == 3? atanhuu(y[2], z[2], prec): atanhui(y[2], z, prec);
901 91 : z = gc_leaf(av, a);
902 91 : if (signe(y) < 0) togglesign(z);
903 : }
904 343 : return z;
905 : }
906 21831 : case t_REAL:
907 21831 : sx = signe(x);
908 21831 : if (!sx) return rcopy(x);
909 21803 : if (expo(x) < 0) return mpatanh(x);
910 :
911 14138 : y = cgetg(3,t_COMPLEX);
912 14138 : av = avma;
913 14138 : z = subrs(x,1);
914 14138 : if (!signe(z)) err_atanh(x, gen_1);
915 14138 : z = invr(z); shiftr_inplace(z, 1); /* 2/(x-1)*/
916 14138 : z = addrs(z,1);
917 14138 : if (!signe(z)) err_atanh(x, gen_m1);
918 14138 : z = logr_abs(z);
919 14138 : shiftr_inplace(z, -1); /* (1/2)log((1+x)/(x-1)) */
920 14138 : gel(y,1) = gc_leaf(av, z);
921 14138 : gel(y,2) = Pi2n(-1, realprec(x));
922 14138 : if (sx > 0) togglesign(gel(y,2));
923 14138 : return y;
924 :
925 46309 : case t_COMPLEX: /* 2/(1-z) - 1 = (1+z) / (1-z) */
926 46309 : if (ismpzero(gel(x,2))) return gatanh(gel(x,1), prec);
927 37428 : av = avma; z = glog( gaddgs(gdivsg(2,gsubsg(1,x)),-1), prec );
928 37428 : return gc_upto(av, gmul2n(z,-1));
929 :
930 35 : default:
931 35 : av = avma; if (!(y = toser_i(x))) break;
932 28 : if (valser(y) < 0) pari_err_DOMAIN("atanh","valuation", "<", gen_0, x);
933 21 : z = gdiv(derivser(y), gsubsg(1,gsqr(y)));
934 14 : a = integser(z);
935 14 : if (!valser(y)) a = gadd(a, gatanh(gel(y,2),prec));
936 14 : return gc_upto(av, a);
937 : }
938 7 : return trans_eval("atanh",gatanh,x,prec);
939 : }
940 : /********************************************************************/
941 : /** **/
942 : /** EULER'S GAMMA **/
943 : /** **/
944 : /********************************************************************/
945 : /* 0 < a < b */
946 : static GEN
947 27481 : mulu_interval_step_i(ulong a, ulong b, ulong step)
948 : {
949 : ulong k, l, N, n;
950 : long lx;
951 : GEN x;
952 :
953 27481 : n = 1 + (b-a) / step;
954 27481 : b -= (b-a) % step;
955 : /* step | b-a */
956 27481 : lx = 1; x = cgetg(2 + n/2, t_VEC);
957 27481 : N = b + a;
958 27481 : for (k = a;; k += step)
959 : {
960 180421 : l = N - k; if (l <= k) break;
961 152940 : gel(x,lx++) = muluu(k,l);
962 : }
963 27481 : if (l == k) gel(x,lx++) = utoipos(k);
964 27481 : setlg(x, lx); return x;
965 : }
966 : static GEN
967 150033 : _mul(void *data, GEN x, GEN y)
968 : {
969 150033 : long prec = (long)data;
970 : /* switch to t_REAL ? */
971 150033 : if (typ(x) == t_INT && lg2prec(lgefint(x)) > prec) x = itor(x, prec);
972 150033 : if (typ(y) == t_INT && lg2prec(lgefint(y)) > prec) y = itor(y, prec);
973 150033 : return mpmul(x, y);
974 : }
975 : static GEN
976 27481 : mulu_interval_step_prec(long l, long m, long s, long prec)
977 : {
978 27481 : GEN v = mulu_interval_step_i(l, m, s);
979 27481 : return gen_product(v, (void*)prec, &_mul);
980 : }
981 :
982 : /* x * (i*(i+1)) */
983 : static GEN
984 7807354 : muliunextu(GEN x, ulong i)
985 : {
986 7807354 : if (i & HIGHMASK) /* i(i+1) >= 2^BITS_IN_LONG*/
987 0 : return mulii(x, muluu(i, i+1));
988 : else
989 7807354 : return muliu(x, i*(i+1));
990 : }
991 : /* arg(s + it); in principle s + it != 0 but currently arg(0. + i*0.) := 0 */
992 : double
993 233372 : dblcarg(double s, double t)
994 : {
995 : double x;
996 233372 : if (!t) return (s >= 0)? 0.: M_PI;
997 202062 : if (!s) return (t >= 0)? M_PI/2: -M_PI/2;
998 202055 : x = atan(t/s);
999 202055 : return (s >= 0)? x: (t >= 0)? x + M_PI : x - M_PI;
1000 : }
1001 :
1002 : /* Let z = s + it, set a = Re(log z), b = Im(log z) */
1003 : void
1004 233368 : dblclog(double s, double t, double *a, double *b)
1005 : {
1006 233368 : *a = log(dblcnorm(s, t)) / 2;
1007 233371 : *b = dblcarg(s, t);
1008 233372 : }
1009 :
1010 : double
1011 16716 : dblcabs(double s, double t) { return sqrt(dblcnorm(s, t)); }
1012 : double
1013 270730 : dblcnorm(double s, double t) { return s*s + t*t; }
1014 :
1015 : #if 0
1016 : /* x, z t_REAL. Compute unique x in ]-z,z] congruent to x mod 2z */
1017 : static GEN
1018 : red_mod_2z(GEN x, GEN z)
1019 : {
1020 : GEN Z = gmul2n(z, 1), d = subrr(z, x);
1021 : /* require little accuracy */
1022 : if (!signe(d)) return x;
1023 : setprec(d, nbits2prec(expo(d) - expo(Z)));
1024 : return addrr(mulir(floorr(divrr(d, Z)), Z), x);
1025 : }
1026 : #endif
1027 :
1028 : static GEN
1029 10661 : negeuler(long prec) { GEN g = mpeuler(prec); setsigne(g, -1); return g; }
1030 : /* lngamma(1+z) = -Euler*z + sum_{i > 1} zeta(i)/i (-z)^i
1031 : * at relative precision prec, |z| <= 1/2 is small */
1032 : static GEN
1033 15656 : lngamma1(GEN z, long prec)
1034 : { /* sum_{i > l} |z|^(i-1) = |z|^l / (1-|z|) < 2^-B
1035 : * for l > (B+1) / |log2(|z|)| */
1036 15656 : long i, l = ceil((prec2nbits(prec) + 1) / - dbllog2(z));
1037 : GEN s, vz;
1038 :
1039 15656 : if (l <= 1) return gmul(negeuler(prec), z);
1040 15481 : vz = constzeta(l, prec);
1041 1053254 : for (i = l, s = gen_0; i > 0; i--)
1042 : {
1043 1037773 : GEN c = divru(gel(vz,i), i);
1044 1039062 : if (odd(i)) setsigne(c, -1);
1045 1039054 : s = gadd(gmul(s,z), c);
1046 : }
1047 15481 : return gmul(z, s);
1048 : }
1049 : /* B_i / (i(i-1)), i even. Sometimes NOT reduced (but gadd/gmul won't care)!*/
1050 : static GEN
1051 7807537 : bern_unextu(long i)
1052 7807537 : { GEN B = bernfrac(i); return mkfrac(gel(B,1), muliunextu(gel(B,2), i-1)); }
1053 : /* B_i / i, i even. Sometimes NOT reduced (but gadd/gmul won't care)!*/
1054 : static GEN
1055 211323 : bern_u(long i)
1056 211323 : { GEN B = bernfrac(i); return mkfrac(gel(B,1), muliu(gel(B,2), i)); }
1057 : /* sum_{i > 0} B_{2i}/(2i(2i-1)) * a^(i-1) */
1058 : static GEN
1059 213302 : lngamma_sum(GEN a, long N)
1060 : {
1061 213302 : pari_sp av = avma;
1062 213302 : GEN S = bern_unextu(2*N);
1063 : long i;
1064 7807559 : for (i = 2*N-2; i > 0; i -= 2)
1065 : {
1066 7594287 : S = gadd(bern_unextu(i), gmul(a,S));
1067 7594234 : if (gc_needed(av,3))
1068 : {
1069 0 : if(DEBUGMEM>1) pari_warn(warnmem,"gamma: i = %ld", i);
1070 0 : S = gc_upto(av, S);
1071 : }
1072 : }
1073 213272 : return S;
1074 : }
1075 : /* sum_{i > 0} B_{2i}/(2i) * a^i */
1076 : static GEN
1077 4249 : psi_sum(GEN a, long N)
1078 : {
1079 4249 : pari_sp av = avma;
1080 4249 : GEN S = bern_u(2*N);
1081 : long i;
1082 211323 : for (i = 2*N-2; i > 0; i -= 2)
1083 : {
1084 207074 : S = gadd(bern_u(i), gmul(a,S));
1085 207074 : if (gc_needed(av,3))
1086 : {
1087 0 : if(DEBUGMEM>1) pari_warn(warnmem,"psi: i = %ld", i);
1088 0 : S = gc_upto(av, S);
1089 : }
1090 : }
1091 4249 : return gmul(a,S);
1092 : }
1093 : static void
1094 225360 : gamma_optim(double ssig, double st, long prec, long *plim, long *pN)
1095 : {
1096 : double la, l,l2,u,v, rlogs, ilogs;
1097 225360 : long N = 1, lim;
1098 225360 : dblclog(ssig,st, &rlogs,&ilogs);
1099 : /* Re (s - 1/2) log(s) */
1100 225364 : u = (ssig - 0.5)*rlogs - st * ilogs;
1101 : /* Im (s - 1/2) log(s) */
1102 225364 : v = (ssig - 0.5)*ilogs + st * rlogs;
1103 : /* l2 = | (s - 1/2) log(s) - s + log(2Pi)/2 |^2 ~ |lngamma(s))|^2 */
1104 225364 : u = u - ssig + log(2.*M_PI)/2;
1105 225364 : v = v - st;
1106 225364 : l2 = u*u + v*v;
1107 225364 : if (l2 < 0.000001) l2 = 0.000001;
1108 225364 : l = (prec2nbits_mul(prec, M_LN2) - log(l2)/2) / 2.;
1109 225364 : if (l < 0) l = 0.;
1110 :
1111 225364 : if (st > 1 && l > 0)
1112 67095 : {
1113 67095 : double t = st * M_PI / l;
1114 67095 : la = t * log(t);
1115 67095 : if (la < 4.) la = 4.;
1116 67095 : if (la > 150) la = t;
1117 : }
1118 : else
1119 158269 : la = 4.; /* heuristic */
1120 225364 : lim = (long)ceil(l / (1.+ log(la)));
1121 225364 : if (lim == 0) lim = 1;
1122 :
1123 225364 : u = (lim-0.5) * la / M_PI;
1124 225364 : l2 = u*u - st*st;
1125 225364 : if (l2 > 0)
1126 : {
1127 212525 : double t = ceil(sqrt(l2) - ssig);
1128 212525 : if (t > 1) N = (long)t;
1129 : }
1130 225364 : *plim = lim; *pN = N;
1131 225364 : }
1132 : /* do we use lngamma1 instead of Euler-Maclaurin ? */
1133 : static int
1134 228153 : gamma_use_1(double s, double t, long prec, long *plim, long *pN)
1135 : {
1136 228153 : double a = s-1, d = fabs(a) + fabs(t);
1137 : long k;
1138 228153 : if (d < 1e-16) return 1;
1139 225360 : gamma_optim(s, t, prec, plim, pN);
1140 225364 : if (d >= 0.5) return 0;
1141 16407 : k = prec2nbits(prec) / -log2(dblcnorm(a, t)); /* 2k = lngamma1 bound */
1142 16408 : return (t ? k: 1.5*k) < *plim + *pN;
1143 : }
1144 : static GEN
1145 228190 : cxgamma(GEN s0, int dolog, long prec)
1146 : {
1147 : GEN s, a, y, res, sig, tau, B, nnx, pi, pi2;
1148 228190 : long i, esig, et, lim, N = 1;
1149 : pari_sp av, av2;
1150 228190 : int funeq = 0;
1151 : pari_timer T;
1152 :
1153 228190 : if (DEBUGLEVEL>5) timer_start(&T);
1154 228190 : s = trans_fix_arg(&prec,&s0,&sig,&tau,&av,&res);
1155 :
1156 228186 : esig = expo(sig);
1157 228186 : et = signe(tau)? expo(tau): 0;
1158 228186 : if ((signe(sig) <= 0 || esig < -1) && et <= 16)
1159 : { /* s <--> 1-s */
1160 21756 : funeq = 1; s = gsubsg(1, s); sig = real_i(s);
1161 : }
1162 :
1163 : /* find "optimal" parameters [lim, N] */
1164 228186 : if (esig > 300 || et > 300)
1165 35 : { /* |s| is HUGE ! Play safe and avoid inf / NaN */
1166 : GEN S, iS, l2, la, u;
1167 : double logla, l;
1168 :
1169 35 : S = gprec_w(s,LOWDEFAULTPREC);
1170 : /* l2 ~ |lngamma(s))|^2 */
1171 35 : l2 = gnorm(gmul(S, glog(S, LOWDEFAULTPREC)));
1172 35 : l = (prec2nbits_mul(prec, M_LN2) - rtodbl(glog(l2,LOWDEFAULTPREC))/2) / 2.;
1173 35 : if (l < 0) l = 0.;
1174 :
1175 35 : iS = imag_i(S);
1176 35 : if (et > 0 && l > 0)
1177 21 : {
1178 21 : GEN t = gmul(iS, dbltor(M_PI / l)), logt = glog(t,LOWDEFAULTPREC);
1179 21 : la = gmul(t, logt);
1180 21 : if (gcmpgs(la, 3) < 0) { logla = log(3.); la = stoi(3); }
1181 14 : else if (gcmpgs(la, 150) > 0) { logla = rtodbl(logt); la = t; }
1182 7 : else logla = rtodbl(mplog(la));
1183 : }
1184 : else
1185 : {
1186 14 : logla = log(3.); la = stoi(3);
1187 : }
1188 35 : lim = (long)ceil(l / (1.+ logla));
1189 35 : if (lim == 0) lim = 1;
1190 :
1191 35 : u = gmul(la, dbltor((lim-0.5)/M_PI));
1192 35 : l2 = gsub(gsqr(u), gsqr(iS));
1193 35 : if (signe(l2) > 0)
1194 : {
1195 14 : l2 = gsub(gsqrt(l2,3), sig);
1196 14 : if (signe(l2) > 0) N = itos( gceil(l2) );
1197 : }
1198 : }
1199 : else
1200 : { /* |s| is moderate. Use floats */
1201 228151 : double ssig = rtodbl(sig);
1202 228152 : double st = typ(s) == t_REAL? 0.0: rtodbl(imag_i(s));
1203 :
1204 228153 : if (gamma_use_1(ssig, st, prec, &lim, &N))
1205 : { /* s ~ 1: loggamma(1+u) ~ - Euler * u, cancellation */
1206 14889 : if (funeq) /* s0 ~ 0: use lngamma(s0)+log(s0) = lngamma(s0+1) */
1207 119 : y = dolog? gsub(lngamma1(s0,prec), glog(s0,prec))
1208 119 : : gdiv(gexp(lngamma1(s0,prec), prec), s0);
1209 : else
1210 : {
1211 14770 : if (isint1(s0))
1212 : {
1213 1683 : set_avma(av);
1214 1683 : return dolog? real_0(prec): real_1(prec);
1215 : }
1216 13087 : y = lngamma1(gsubgs(s0,1),prec);
1217 13087 : if (!dolog) y = gexp(y,prec);
1218 : }
1219 13206 : set_avma(av); return affc_fixlg(y, res);
1220 : }
1221 : }
1222 213303 : if (DEBUGLEVEL>5) err_printf("lim, N: [%ld, %ld]\n",lim,N);
1223 213303 : incrprec(prec);
1224 :
1225 213303 : av2 = avma;
1226 213303 : y = s;
1227 213303 : if (typ(s0) == t_INT)
1228 : {
1229 2591 : ulong ss = itou_or_0(s0);
1230 2591 : if (signe(s0) <= 0)
1231 0 : pari_err_DOMAIN("gamma","argument", "=",
1232 : strtoGENstr("nonpositive integer"), s0);
1233 2591 : if (!ss || ss + (ulong)N < ss) {
1234 7 : for (i=1; i < N; i++)
1235 : {
1236 0 : y = mulri(y, addiu(s0, i));
1237 0 : if (gc_needed(av2,3))
1238 : {
1239 0 : if(DEBUGMEM>1) pari_warn(warnmem,"gamma");
1240 0 : y = gc_leaf(av2, y);
1241 : }
1242 : }
1243 : } else {
1244 33846 : for (i=1; i < N; i++)
1245 : {
1246 31262 : y = mulru(y, ss + i);
1247 31262 : if (gc_needed(av2,3))
1248 : {
1249 0 : if(DEBUGMEM>1) pari_warn(warnmem,"gamma");
1250 0 : y = gc_leaf(av2, y);
1251 : }
1252 : }
1253 : }
1254 : }
1255 : else
1256 : { /* Compute lngamma mod 2 I Pi */
1257 210712 : GEN sq = gsqr(s);
1258 210712 : pari_sp av3 = avma;
1259 4306170 : for (i = 1; i < N - 1; i += 2)
1260 : {
1261 4095469 : y = gmul(y, gaddsg(i*(i + 1), gadd(gmulsg(2*i + 1, s), sq)));
1262 4095463 : if (gc_needed(av2,3))
1263 : {
1264 0 : if(DEBUGMEM>1) pari_warn(warnmem,"gamma");
1265 0 : y = gc_upto(av3, y);
1266 : }
1267 : }
1268 210701 : if (!odd(N)) y = gmul(y, gaddsg(N - 1, s));
1269 : }
1270 213301 : if (DEBUGLEVEL>5) timer_printf(&T,"product from 0 to N-1");
1271 213301 : constbern(lim);
1272 213301 : nnx = gaddgs(s, N); a = ginv(nnx);
1273 213300 : B = gadd(gsub(gmul(gsub(nnx, ghalf), glog(nnx,prec)), nnx),
1274 : gmul(a, lngamma_sum(gsqr(a), lim)));
1275 213299 : if (DEBUGLEVEL>5) timer_printf(&T,"Bernoulli sum");
1276 :
1277 213299 : pi = mppi(prec); pi2 = shiftr(pi, 1);
1278 213303 : if (dolog)
1279 : {
1280 15750 : if (typ(s) == t_REAL)
1281 : {
1282 12411 : if (!funeq) y = logr_abs(divrr(sqrtr(pi2), y));
1283 : else
1284 : {
1285 7 : GEN T = shiftr(sqrtr(pi2),-1); /* sqrt(Pi/2) */
1286 : /* s0 < 0, step (*) simplifies: imag(lngamma(s0)) = - Pi * floor(s0) */
1287 7 : y = logr_abs(divrr(mulrr(y, T), mpsin(gmul(pi,s0))));
1288 7 : y = mkcomplex(y, mulri(pi, gfloor(s0)));
1289 7 : B = gneg(B);
1290 : }
1291 : }
1292 : else
1293 : { /* log(y), fixing imaginary part */
1294 3339 : long prec2 = LOWDEFAULTPREC;
1295 3339 : GEN k, s2 = gprec_w(s, prec2), y2 = garg(s2, prec2); /* ~ Im log(s) */
1296 10424 : for (i=1; i < N; i++) y2 = gadd(y2, garg(gaddgs(s2,i), prec2));
1297 3339 : y = glog(y, prec);
1298 3339 : k = ground( gdiv(gsub(y2, imag_i(y)), Pi2n(1,prec2)) );
1299 3339 : if (signe(k)) y = gadd(y, mulcxI(mulir(k, Pi2n(1, prec))));
1300 3339 : if (!funeq) y = gsub(shiftr(logr_abs(pi2),-1), y); /* y -> sqrt(2Pi)/y */
1301 : else
1302 : { /* recall that s = 1 - s0 */
1303 273 : GEN T = shiftr(sqrtr(pi2),-1); /* sqrt(Pi/2) */
1304 : /* (*) Compute log(sin(Pi s0)) so that it has branch cuts along
1305 : * (-oo, 0] and [1, oo). To do this in a numerically stable way
1306 : * we must compute the log first then mangle its imaginary part.
1307 : * The rounding operation below is stable because we're rounding
1308 : * a number which is already within 1/4 of an integer. */
1309 :
1310 : /* z = log(sin(Pi s0) / sqrt(Pi/2)) */
1311 273 : GEN z = glog(gdiv(gsin(gmul(pi,s0),prec), T), prec);
1312 273 : GEN b = shiftr(subrs(shiftr(sig, 1), 1), -2); /* (2 Re(s)-1) / 4 */
1313 273 : y = gsub(y, z);
1314 273 : if (gsigne(imag_i(s)) > 0) togglesign(b);
1315 273 : z = roundr(gsub(gdiv(imag_i(z), pi2), b)); /* round( Im(z)/2Pi - b ) */
1316 273 : if (signe(z)) { /* y += I*z, z a t_REAL */
1317 0 : z = mulir(z, pi2);
1318 0 : if (typ(y) == t_COMPLEX) gel(y,2) = gadd(gel(y,2), z);
1319 0 : else y = mkcomplex(y, z);
1320 : }
1321 273 : B = gneg(B);
1322 : }
1323 : }
1324 15750 : y = gadd(B, y);
1325 : }
1326 : else
1327 : {
1328 197553 : GEN sqrtpi2 = sqrtr(pi2);
1329 197551 : if (funeq)
1330 : { /* y --> y Pi/(sin(Pi s) * sqrt(2Pi)) = y sqrt(Pi/2)/sin(Pi s) */
1331 21357 : y = gdiv(gmul(shiftr(sqrtpi2,-1),y), gsin(gmul(pi,s0), prec));
1332 : /* don't use s above: sin(pi s0) = sin(pi s) and the former is
1333 : * more accurate, esp. if s0 ~ 0 */
1334 21357 : B = gneg(B);
1335 : }
1336 : else /* y --> sqrt(2Pi) / y */
1337 176194 : y = gdiv(sqrtpi2, y);
1338 197553 : y = gmul(gexp(B, prec), y);
1339 : }
1340 213303 : set_avma(av); return affc_fixlg(y, res);
1341 : }
1342 :
1343 : /* Theory says n > C * b^1.5 / log(b). Timings:
1344 : * b = 64*[1, 2, 3, 4, 5, 6, 7, 10, 20, 30, 50, 100, 200, 500];
1345 : * n = [1450, 1930, 2750, 3400, 4070, 5000, 6000, 8800, 26000, 50000, 130000,
1346 : * 380000, 1300000, 6000000]; */
1347 : static long
1348 37779 : gamma2_n(long prec)
1349 : {
1350 37779 : long b = prec2nbits(prec);
1351 37779 : if (b <= 64) return 1450;
1352 37156 : if (b <= 128) return 1930;
1353 31171 : if (b <= 192) return 2750;
1354 16963 : if (b <= 256) return 3400;
1355 7341 : if (b <= 320) return 4070;
1356 6925 : if (b <= 384) return 5000;
1357 4115 : if (b <= 448) return 6000;
1358 3933 : return 10.0 * b * sqrt(b) / log(b);
1359 : }
1360 :
1361 : /* m even, Gamma((m+1) / 2) */
1362 : static GEN
1363 37779 : gammahs(long m, long prec)
1364 : {
1365 37779 : GEN y = cgetr(prec), z;
1366 37779 : pari_sp av = avma;
1367 37779 : long ma = labs(m);
1368 :
1369 37779 : if (ma > gamma2_n(prec))
1370 : {
1371 0 : z = stor(m + 1, prec); shiftr_inplace(z, -1);
1372 0 : affrr(cxgamma(z,0,prec), y);
1373 0 : set_avma(av); return y;
1374 : }
1375 37779 : z = sqrtr( mppi(prec) );
1376 37779 : if (m)
1377 : {
1378 23030 : GEN t = mulu_interval_step_prec(1, ma-1, 2, prec + EXTRAPREC64);
1379 23030 : if (m >= 0) z = mpmul(z,t);
1380 : else
1381 : {
1382 217 : z = mpdiv(z,t);
1383 217 : if ((m&3) == 2) setsigne(z,-1);
1384 : }
1385 23030 : shiftr_inplace(z, -m/2);
1386 : }
1387 37779 : affrr(z, y); set_avma(av); return y;
1388 : }
1389 : GEN
1390 28 : ggammah(GEN x, long prec)
1391 : {
1392 28 : switch(typ(x))
1393 : {
1394 21 : case t_INT:
1395 : {
1396 21 : long k = itos_or_0(x);
1397 21 : if (!k && signe(x)) pari_err_OVERFLOW("gamma");
1398 21 : return gammahs(k * 2, prec);
1399 : }
1400 7 : case t_REAL: case t_COMPLEX: case t_PADIC: case t_SER: {
1401 7 : pari_sp av = avma;
1402 7 : return gc_upto(av, ggamma(gadd(x,ghalf), prec));
1403 : }
1404 : }
1405 0 : return trans_eval("gammah",ggammah,x,prec);
1406 : }
1407 :
1408 : /* find n such that n+v_p(n!)>=k p^2/(p-1)^2 */
1409 : static long
1410 215746 : nboft(long k, long p)
1411 : {
1412 215746 : pari_sp av = avma;
1413 : long s, n;
1414 :
1415 215746 : if (k <= 0) return 0;
1416 215746 : k = itou( gceil(gdiv(mului(k, sqru(p)), sqru(p-1))) );
1417 215746 : set_avma(av);
1418 2526574 : for (s=0, n=0; n+s < k; n++, s += u_lval(n, p));
1419 215746 : return n;
1420 : }
1421 :
1422 : /* Using Dwork's expansion, compute \Gamma(px+1)=-\Gamma(px) with x a unit.
1423 : * See p-Adic Gamma Functions and Dwork Cohomology, Maurizio Boyarsky
1424 : * Transactions of the AMS, Vol. 257, No. 2. (Feb., 1980), pp. 359-369.
1425 : * Inspired by a GP script by Fernando Rodriguez-Villegas */
1426 : static GEN
1427 215746 : gadw(GEN x, long p)
1428 : {
1429 215746 : pari_sp ltop = avma;
1430 215746 : GEN s, t, u = cgetg(p+1, t_VEC);
1431 215746 : long j, k, kp, n = nboft(precp(x)+valp(x)+1, p);
1432 :
1433 215746 : t = s = cvtop(gen_1, padic_p(x), n);
1434 215748 : gel(u, 1) = s;
1435 215748 : gel(u, 2) = s;
1436 891733 : for (j = 2; j < p; ++j)
1437 675987 : gel(u, j+1) = gdivgu(gel(u, j), j);
1438 2310828 : for (k = 1, kp = p; k < n; ++k, kp += p) /* kp = k*p */
1439 : {
1440 : GEN c;
1441 2095082 : gel(u, 1) = gdivgu(gadd(gel(u, 1), gel(u, p)), kp);
1442 10344742 : for (j = 1; j < p; ++j)
1443 8249660 : gel(u, j+1) = gdivgu(gadd(gel(u, j), gel(u, j+1)), kp + j);
1444 :
1445 2095082 : t = gmul(t, gaddgs(x, k-1));
1446 2095081 : c = leafcopy(gel(u,1)); setvalp(c, valp(c) + k); /* c = u[1] * p^k */
1447 2095082 : s = gadd(s, gmul(c, t));
1448 2095082 : if ((k&0xFL)==0) (void)gc_all(ltop, 3, &u,&s,&t);
1449 : }
1450 215746 : return gneg(s);
1451 : }
1452 :
1453 : /*Use Dwork expansion*/
1454 : /*This is a O(p*e*log(pe)) algorithm, should be used when p small
1455 : * If p==2 this is a O(pe) algorithm. */
1456 : static GEN
1457 215746 : Qp_gamma_Dwork(GEN x, long p)
1458 : {
1459 215746 : pari_sp ltop = avma;
1460 215746 : long k = padic_to_Fl(x, p);
1461 : GEN p1;
1462 : long j;
1463 215746 : long px = precp(x);
1464 215746 : if (p==2 && px)
1465 : {
1466 3010 : x = shallowcopy(x);
1467 3010 : setprecp(x, px+1);
1468 3010 : padic_pd(x) = shifti(padic_pd(x), 1);
1469 : }
1470 215746 : if (k)
1471 : {
1472 170036 : GEN x_k = gsubgs(x,k);
1473 170036 : x = gdivgu(x_k, p);
1474 170036 : p1 = gadw(x, p); if (!odd(k)) p1 = gneg(p1);
1475 447721 : for (j = 1; j < k; ++j) p1 = gmul(p1, gaddgs(x_k, j));
1476 : }
1477 : else
1478 45710 : p1 = gneg(gadw(gdivgu(x, p), p));
1479 215750 : return gc_upto(ltop, p1);
1480 : }
1481 :
1482 : /* Compute Qp_gamma using the definition. This is a O(x*M(log(pe))) algorithm.
1483 : * This should be used if x is very small. */
1484 : static GEN
1485 490 : Qp_gamma_Morita(long n, GEN p, long e)
1486 : {
1487 490 : pari_sp av = avma;
1488 490 : GEN p2 = cvtop((n&1)? gen_m1: gen_1, p, e);
1489 490 : long i, pp = is_bigint(p)? 0: itos(p);
1490 7749 : for (i = 2; i < n; i++)
1491 7259 : if (!pp || i%pp)
1492 : {
1493 5215 : p2 = gmulgu(p2, i);
1494 5215 : if ((i&0xFL) == 0xFL) p2 = gc_upto(av, p2);
1495 : }
1496 490 : return gc_upto(av, p2);
1497 : }
1498 :
1499 : /* x\in\N: Gamma(-x)=(-1)^(1+x+x\p)*Gamma(1+x) */
1500 : static GEN
1501 238 : Qp_gamma_neg_Morita(long n, GEN p, long e)
1502 : {
1503 238 : GEN g = ginv(Qp_gamma_Morita(n+1, p, e));
1504 238 : return ((n^sdivsi(n,p)) & 1)? g: gneg(g);
1505 : }
1506 :
1507 : /* p-adic Gamma function for x a p-adic integer */
1508 : /* If n < p*e : use Morita's definition.
1509 : * Else : use Dwork's expansion.
1510 : * If both n and p are big : itos(p) will fail.
1511 : * TODO: handle p=2 better (Qp_gamma_Dwork is slow for p=2). */
1512 : GEN
1513 216243 : Qp_gamma(GEN x)
1514 : {
1515 216243 : GEN n, m, N, p = padic_p(x);
1516 216243 : long s, e = valp(x) + precp(x);
1517 216243 : if (absequaliu(p, 2) && e == 2) e = 1;
1518 216243 : if (valp(x) < 0) pari_err_DOMAIN("gamma","v_p(x)", "<", gen_0, x);
1519 216236 : n = gtrunc(x);
1520 216236 : m = gtrunc(gneg(x));
1521 216236 : N = cmpii(n,m)<=0?n:m;
1522 216236 : s = itos_or_0(N);
1523 216236 : if (s && cmpsi(s, muliu(p,e)) < 0) /* s < p*e */
1524 490 : return (N == n) ? Qp_gamma_Morita(s,p,e): Qp_gamma_neg_Morita(s,p,e);
1525 215746 : return Qp_gamma_Dwork(x, itos(p));
1526 : }
1527 :
1528 : static GEN
1529 14 : Qp_lngamma(GEN x)
1530 : {
1531 : GEN s, y, Y;
1532 14 : long v = valp(x), e, k, K;
1533 14 : if (v >= 0) return Qp_log(Qp_gamma(x));
1534 7 : e = precp(x) + v; K = (2 + (e + 4) / (-v)) >> 1;
1535 7 : s = gen_0; Y = y = ginv(x); y = gsqr(y); constbern(K);
1536 63 : for (k = 1; k <= K; k++)
1537 : {
1538 56 : s = gadd(s, gmul(gdivgunextu(bernfrac(2*k), 2*k-1), Y));
1539 56 : if (k < K) Y = gmul(Y, y); /* x^(1-2k) */
1540 : }
1541 7 : return gadd(s, gsub(gmul(gsub(x, ghalf), Qp_log(x)), x));
1542 : }
1543 :
1544 : /* gamma(1+x) - 1, |x| < 1 is "small" */
1545 : GEN
1546 1211 : ggamma1m1(GEN x, long prec) { return gexpm1(lngamma1(x, prec), prec); }
1547 :
1548 : /* lngamma(y) with 0 constant term, using (lngamma y)' = y' psi(y) */
1549 : static GEN
1550 27755 : serlngamma0(GEN y, long prec)
1551 : {
1552 : GEN t;
1553 27755 : if (valser(y)) pari_err_DOMAIN("lngamma","valuation", "!=", gen_0, y);
1554 27748 : t = derivser(y);
1555 : /* don't compute psi if y'=0 */
1556 27748 : if (signe(t)) t = gmul(t, gpsi(y,prec));
1557 27748 : return integser(t);
1558 : }
1559 :
1560 : static GEN
1561 27720 : sergamma(GEN y, long prec)
1562 : {
1563 : GEN z, y0, Y;
1564 27720 : if (lg(y) == 2) pari_err_DOMAIN("gamma", "argument", "=", gen_0,y);
1565 : /* exp(lngamma) */
1566 27713 : if (valser(y) > 0) return gdiv(gexp(glngamma(gaddgs(y,1),prec),prec),y);
1567 27426 : y0 = simplify_shallow(gel(y,2));
1568 27426 : z = NULL; Y = y;
1569 27426 : if (isint(y0, &y0))
1570 : { /* fun eq. avoids log singularity of lngamma at negative ints */
1571 13545 : long s = signe(y0);
1572 : /* possible if y[2] is an inexact 0 */
1573 13545 : if (!s) return gdiv(gexp(glngamma(gaddgs(y,1),prec),prec),y);
1574 13538 : if (signe(y0) < 0) { Y = gsubsg(1, y); y0 = subsi(1, y0); }
1575 13538 : if (abscmpiu(y0, 50) < 0) z = mpfact(itos(y0)-1); /* more precise */
1576 : }
1577 27419 : if (!z) z = ggamma(y0,prec);
1578 27419 : z = gmul(z, gexp(serlngamma0(Y,prec),prec));
1579 27419 : if (Y != y)
1580 : {
1581 98 : GEN pi = mppi(prec);
1582 98 : z = gdiv(mpodd(y0)? pi: negr(pi),
1583 : gmul(z, gsin(gmul(pi,serchop0(y)), prec)));
1584 : }
1585 27419 : return z;
1586 : }
1587 :
1588 : static GEN
1589 9422 : sqrtu(ulong a, long prec) { return sqrtr_abs(utor(a, prec)); }
1590 : static GEN
1591 245 : cbrtu(ulong a, long prec) { return sqrtnr_abs(utor(a, prec), 3); }
1592 : /* N | 6 */
1593 : static GEN
1594 6020 : ellkprime(long N, GEN s2, GEN s3)
1595 : {
1596 : GEN z;
1597 6020 : switch(N)
1598 : {
1599 2072 : case 1: return shiftr(s2, -1);
1600 49 : case 2: return sqrtr_abs(shiftr(subrs(s2,1), 1));
1601 3794 : case 3: return shiftr(mulrr(s2, addrs(s3, 1)), -2);
1602 105 : default: /* 6 */
1603 105 : z = mulrr(subrr(s3,s2), subsr(2,s3));
1604 105 : return mulrr(addsr(2,s2), sqrtr_abs(z));
1605 : }
1606 : }
1607 :
1608 : static GEN
1609 6020 : ellKk(long N, GEN s2, GEN s3, long prec)
1610 6020 : { return gdiv(Pi2n(-1,prec), agm(ellkprime(N,s2,s3), gen_1, prec)); }
1611 :
1612 : /* Gamma(1/3) */
1613 : static GEN
1614 3689 : G3(GEN s2, GEN s3, long prec)
1615 : {
1616 3689 : GEN A = ellKk(3, s2,s3, prec), pi = mppi(prec);
1617 3689 : A = shiftr(divrs(powrs(mulrr(pi, A), 12), 27), 28);
1618 3689 : return sqrtnr_abs(A, 36);
1619 : }
1620 : /* Gamma(1/4) */
1621 : static GEN
1622 1918 : G4(GEN s2, long prec)
1623 : {
1624 1918 : GEN A = ellKk(1, s2,NULL, prec), pi = mppi(prec);
1625 1918 : return shiftr(sqrtr_abs(mulrr(sqrtr_abs(pi), A)), 1);
1626 : }
1627 :
1628 : /* Gamma(n / 24), n = 1,5,7,11 */
1629 : static GEN
1630 105 : Gn24(long n, GEN s2, GEN s3, long prec)
1631 : {
1632 105 : GEN A, B, C, t, t1, t2, t3, t4, pi = mppi(prec);
1633 105 : A = ellKk(1, s2,s3, prec);
1634 105 : B = ellKk(3, s2,s3, prec);
1635 105 : C = ellKk(6, s2,s3, prec);
1636 105 : t1 = sqrtr_abs(mulur(3, addsr(2, s3)));
1637 105 : t2 = sqrtr_abs(divrr(s3, pi));
1638 105 : t2 = mulrr(t2, shiftr(mulrr(addrr(s2,s3), A), 2));
1639 105 : t3 = mulrr(divur(3,pi), sqrr(B));
1640 105 : t3 = mulrr(addsr(2,s2), sqrtnr_abs(shiftr(powrs(t3, 3), 8), 9));
1641 105 : t4 = mulrr(mulrr(addsr(1, s2), subrr(s3, s2)), subsr(2, s3));
1642 105 : t4 = mulrr(mulrr(mulur(384, t4), pi), sqrr(C));
1643 105 : switch (n)
1644 : {
1645 63 : case 1: t = mulrr(mulrr(t1, t2), mulrr(t3, t4)); break;
1646 14 : case 5: t = divrr(mulrr(t2, t4), mulrr(t1, t3)); break;
1647 14 : case 7: t = divrr(mulrr(t3, t4), mulrr(t1, t2)); break;
1648 14 : default:t = divrr(mulrr(t1, t4), mulrr(t2, t3)); break;
1649 : }
1650 105 : return sqrtnr_abs(t, 4);
1651 : }
1652 : /* sin(x/2) = sqrt((1-c) / 2) > 0 given c = cos(x) */
1653 : static GEN
1654 28 : sinx2(GEN c)
1655 28 : { c = subsr(1, c); shiftr_inplace(c,-1); return sqrtr_abs(c); }
1656 : /* sin(Pi/12), given sqrt(3) */
1657 : static GEN
1658 21 : sin12(GEN s3)
1659 21 : { GEN t = subsr(2, s3); shiftr_inplace(t, -2); return sqrtr_abs(t); }
1660 : /* cos(Pi/12) = sin(5Pi/12), given sqrt(3) */
1661 : static GEN
1662 49 : cos12(GEN s3)
1663 49 : { GEN t = addsr(2, s3); shiftr_inplace(t, -2); return sqrtr_abs(t); }
1664 : /* 0 < n < d, (n, d) = 1, 2 < d | 24 */
1665 : static GEN
1666 5628 : gammafrac24_s(long n, long d, long prec)
1667 : {
1668 5628 : GEN A, B, s2, s3, pi = mppi(prec);
1669 5628 : s2 = sqrtu(2, prec);
1670 5628 : s3 = d % 3? NULL: sqrtu(3, prec);
1671 5628 : switch(d)
1672 : {
1673 3311 : case 3:
1674 3311 : A = G3(s2,s3,prec);
1675 3311 : if (n == 1) return A;
1676 2849 : return divrr(Pi2n(1, prec), mulrr(s3, A));
1677 1785 : case 4:
1678 1785 : A = G4(s2,prec);
1679 1785 : if (n == 1) return A;
1680 1183 : return divrr(mulrr(pi, s2), A);
1681 245 : case 6:
1682 245 : A = sqrr(G3(s2,s3,prec));
1683 245 : A = mulrr(A, sqrtr_abs(divsr(3, pi)));
1684 245 : A = divrr(A, cbrtu(2, prec));
1685 245 : if (n == 1) return A;
1686 140 : return divrr(Pi2n(1, prec), A);
1687 49 : case 8:
1688 49 : A = ellKk(1, s2,s3, prec);
1689 49 : B = ellKk(2, s2,s3, prec);
1690 49 : A = shiftr(sqrtr_abs(divrr(mulrr(addsr(1, s2), A), sqrtr_abs(pi))), 1);
1691 49 : B = shiftr(mulrr(sqrtr_abs(gmul(subrs(s2, 1), mulrr(s2, pi))), B), 3);
1692 49 : switch (n)
1693 : {
1694 : GEN t;
1695 28 : case 1: return sqrtr_abs(mulrr(A, B));
1696 7 : case 3: return sqrtr_abs(divrr(B, A));
1697 7 : case 5: A = sqrtr_abs(divrr(B, A));
1698 7 : t = sqrtr_abs(shiftr(addsr(1, shiftr(s2, -1)), -1)); /*sin(3Pi/8)*/
1699 7 : return divrr(pi, mulrr(t, A));
1700 7 : default: A = sqrtr_abs(mulrr(A, B));
1701 7 : t = sqrtr_abs(shiftr(subsr(1, shiftr(s2, -1)), -1)); /*sin(Pi/8)*/
1702 7 : return divrr(pi, mulrr(t, A));
1703 : }
1704 133 : case 12:
1705 133 : A = G3(s2,s3,prec);
1706 133 : B = G4(s2,prec);
1707 : switch (n)
1708 : {
1709 : GEN t2;
1710 77 : case 1: case 11:
1711 77 : t2 = shiftr(mulur(27, powrs(divrr(addsr(1,s3), pi), 4)), -2);
1712 77 : t2 = mulrr(sqrtnr_abs(t2, 8), mulrr(A, B));
1713 77 : if (n == 1) return t2;
1714 7 : return divrr(pi, mulrr(sin12(s3), t2));
1715 56 : case 5: case 7:
1716 56 : t2 = shiftr(divrs(powrs(mulrr(subrs(s3,1), pi), 4), 3), 2);
1717 56 : t2 = mulrr(sqrtnr_abs(t2, 8), divrr(B, A));
1718 56 : if (n == 5) return t2;
1719 35 : return divrr(pi, mulrr(cos12(s3), t2));
1720 : }
1721 : default: /* n = 24 */
1722 105 : if (n > 12)
1723 : {
1724 : GEN t;
1725 28 : n = 24 - n;
1726 28 : A = Gn24(n, s2,s3, prec);
1727 : switch(n)
1728 : { /* t = sin(n*Pi/24) */
1729 7 : case 1: t = cos12(s3); t = sinx2(t); break;
1730 7 : case 5: t = sin12(s3); t = sinx2(t); break;
1731 7 : case 7: t = sin12(s3); togglesign(t); t = sinx2(t); break;
1732 7 : default:t = cos12(s3); togglesign(t); t = sinx2(t); break; /* n=11 */
1733 : }
1734 28 : return divrr(pi, mulrr(A, t));
1735 : }
1736 77 : return Gn24(n, s2,s3, prec);
1737 : }
1738 : }
1739 :
1740 : /* (a,b) = 1. If 0 < x < b, m >= 0
1741 : gamma(x/b + m) = gamma(x/b) * mulu_interval_step(x, x+(m-1)*b, b) / b^m
1742 : gamma(x/b - m) = gamma(x/b) / mulu_interval_step(b-x, b*m-x, b) * (-b)^m */
1743 : static GEN
1744 46102 : gammafrac24(GEN a, GEN b, long prec)
1745 : {
1746 : pari_sp av;
1747 : long A, B, m, x, bit;
1748 : GEN z0, z, t;
1749 46102 : if (!(A = itos_or_0(a)) || !(B = itos_or_0(b)) || B > 24) return NULL;
1750 44737 : switch(B)
1751 : {
1752 37758 : case 2: return gammahs(A-1, prec);
1753 5642 : case 3: case 4: case 6: case 8: case 12: case 24:
1754 5642 : m = A / B;
1755 5642 : x = A % B; /* = A - m*B */
1756 5642 : if (x < 0) { x += B; m--; } /* now 0 < x < B, A/B = x/B + m */
1757 5642 : bit = prec2nbits(prec);
1758 : /* Depending on B and prec, we must experimentally replace the 0.5
1759 : * by 0.4 to 2.0 for optimal value. Play safe. */
1760 5642 : if (labs(m) > 0.5 * bit * sqrt(bit) / log(bit)) return NULL;
1761 5628 : z0 = cgetr(prec); av = avma;
1762 5628 : prec += EXTRAPREC64;
1763 5628 : z = gammafrac24_s(x, B, prec);
1764 5628 : if (m)
1765 : {
1766 3808 : if (m > 0)
1767 3759 : t = mpdiv(mulu_interval_step(x, (m-1)*B + x, B), rpowuu(B,m,prec));
1768 : else
1769 : {
1770 49 : m = -m;
1771 49 : t = mpdiv(rpowuu(B,m,prec), mulu_interval_step(B-x, m*B - x, B));
1772 49 : if (odd(m)) togglesign(t);
1773 : }
1774 3808 : z = mpmul(z,t);
1775 : }
1776 5628 : affrr(z, z0); set_avma(av); return z0;
1777 : }
1778 1337 : return NULL;
1779 : }
1780 : GEN
1781 634198 : ggamma(GEN x, long prec)
1782 : {
1783 : pari_sp av;
1784 : GEN y;
1785 :
1786 634198 : switch(typ(x))
1787 : {
1788 364337 : case t_INT:
1789 364337 : if (signe(x) <= 0)
1790 0 : pari_err_DOMAIN("gamma","argument", "=",
1791 : strtoGENstr("nonpositive integer"), x);
1792 364337 : return mpfactr(itos(x) - 1, prec);
1793 :
1794 204517 : case t_REAL: case t_COMPLEX:
1795 204517 : return cxgamma(x, 0, prec);
1796 :
1797 37541 : case t_FRAC:
1798 : {
1799 37541 : GEN a = gel(x,1), b = gel(x,2), c = gammafrac24(a, b, prec);
1800 37541 : if (c) return c;
1801 1652 : av = avma; c = subii(a,b);
1802 1652 : if (signe(a) < 0)
1803 : { /* gamma will use functional equation x -> z = 1-x = -c/b >= 1/2.
1804 : * Gamma(x) = Pi / (sin(Pi z) * Gamma(z)) */
1805 49 : GEN z = mkfrac(negi(c), b), q = ground(z), r = gsub(z,q);
1806 49 : GEN pi = mppi(prec); /* |r| <= 1/2 */
1807 49 : z = fractor(z, prec+EXTRAPREC64);
1808 49 : y = divrr(pi, mulrr(mpsin(gmul(pi, r)), cxgamma(z, 0, prec)));
1809 49 : if (mpodd(q)) togglesign(y);
1810 49 : return gc_upto(av, y);
1811 : }
1812 1603 : if (cmpii(shifti(a,1), b) < 0)
1813 : { /* 0 < x < 1/2 gamma would use funeq: adding 1 is cheaper. */
1814 343 : if (expi(a) - expi(b) < -3) /* close to 0 */
1815 : {
1816 14 : if (lg2prec(lgefint(b)) >= prec) x = fractor(x,prec);
1817 14 : y = mpexp(lngamma1(x, prec));
1818 : }
1819 : else
1820 329 : y = cxgamma(fractor(mkfrac(addii(a,b), b), prec), 0, prec);
1821 343 : return gc_upto(av, gdiv(y, x));
1822 : }
1823 1260 : if (expi(c) - expi(b) < -3)
1824 : { /* x = 1 + c/b is close to 1 */
1825 336 : x = mkfrac(c,b);
1826 336 : if (lg2prec(lgefint(b)) >= prec) x = fractor(x,prec);
1827 336 : y = mpexp(lngamma1(x, prec));
1828 : }
1829 : else
1830 924 : y = cxgamma(fractor(x, prec), 0, prec);
1831 1260 : return gc_upto(av, y);
1832 : }
1833 :
1834 84 : case t_PADIC: return Qp_gamma(x);
1835 27719 : default:
1836 27719 : av = avma; if (!(y = toser_i(x))) break;
1837 27720 : return gc_upto(av, sergamma(y, prec));
1838 : }
1839 0 : return trans_eval("gamma",ggamma,x,prec);
1840 : }
1841 :
1842 : static GEN
1843 524 : mpfactr_basecase(long n, long prec)
1844 : {
1845 524 : GEN v = cgetg(expu(n) + 2, t_VEC);
1846 524 : long k, prec2 = prec + EXTRAPREC64;
1847 : GEN a;
1848 524 : for (k = 1;; k++)
1849 4451 : {
1850 4975 : long m = n >> (k-1), l;
1851 4975 : if (m <= 2) break;
1852 4451 : l = (1 + (n >> k)) | 1;
1853 : /* product of odd numbers in ]n / 2^k, 2 / 2^(k-1)] */
1854 4451 : a = mulu_interval_step_prec(l, m, 2, prec2);
1855 4451 : gel(v,k) = k == 1? a: gpowgs(a, k);
1856 : }
1857 4451 : a = gel(v,--k); while (--k) a = mpmul(a, gel(v,k));
1858 524 : if (typ(a) == t_INT) a = itor(a, prec); else a = gprec_wtrunc(a, prec);
1859 524 : shiftr_inplace(a, factorial_lval(n, 2));
1860 524 : return a;
1861 : }
1862 : /* Theory says n > C * b^1.5 / log(b). Timings:
1863 : * b = [64, 128, 192, 256, 512, 1024, 2048, 4096, 8192, 16384]
1864 : * n = [1930, 2650, 3300, 4270, 9000, 23000, 75000, 210000, 750000, 2400000] */
1865 : static long
1866 643 : mpfactr_n(long prec)
1867 : {
1868 643 : long b = prec2nbits(prec);
1869 643 : if (b <= 64) return 1930;
1870 97 : if (b <= 128) return 2650;
1871 76 : if (b <= 192) return 3300;
1872 76 : return b * sqrt(b);
1873 : }
1874 : static GEN
1875 7889 : mpfactr_small(long n, long prec)
1876 : {
1877 7889 : GEN f = cgetr(prec);
1878 7889 : pari_sp av = avma;
1879 7889 : if (n < 410)
1880 7889 : affir(mpfact(n), f);
1881 : else
1882 0 : affrr(mpfactr_basecase(n, prec), f);
1883 7889 : set_avma(av); return f;
1884 : }
1885 : GEN
1886 407955 : mpfactr(long n, long prec)
1887 : {
1888 407955 : GEN f = cgetr(prec);
1889 407955 : pari_sp av = avma;
1890 :
1891 407955 : if (n < 410)
1892 407312 : affir(mpfact(n), f);
1893 : else
1894 : {
1895 643 : long N = mpfactr_n(prec);
1896 524 : GEN z = n <= N? mpfactr_basecase(n, prec)
1897 643 : : cxgamma(utor(n+1, prec), 0, prec);
1898 643 : affrr(z, f);
1899 : }
1900 407955 : set_avma(av); return f;
1901 : }
1902 :
1903 : /* First a little worse than mpfactr_n because of the extra logarithm.
1904 : * Asymptotically same. */
1905 : static ulong
1906 7889 : lngamma_n(long prec)
1907 : {
1908 7889 : long b = prec2nbits(prec);
1909 : double N;
1910 7889 : if (b <= 64) return 1450UL;
1911 7889 : if (b <= 128) return 2010UL;
1912 308 : if (b <= 192) return 2870UL;
1913 308 : N = b * sqrt(b);
1914 308 : if (b <= 256) return N/1.25;
1915 0 : if (b <= 512) return N/1.2;
1916 0 : if (b <= 2048) return N/1.1;
1917 0 : return N;
1918 : }
1919 :
1920 : GEN
1921 38885 : glngamma(GEN x, long prec)
1922 : {
1923 38885 : pari_sp av = avma;
1924 : GEN y, y0, t;
1925 :
1926 38885 : switch(typ(x))
1927 : {
1928 7896 : case t_INT:
1929 : {
1930 : ulong n;
1931 7896 : if (signe(x) <= 0)
1932 0 : pari_err_DOMAIN("lngamma","argument", "=",
1933 : strtoGENstr("nonpositive integer"), x);
1934 7896 : n = itou_or_0(x);
1935 7896 : if (!n || n > lngamma_n(prec)) return cxgamma(x, 1, prec);
1936 7889 : return gc_leaf(av, logr_abs( mpfactr_small(n-1, prec) ));
1937 : }
1938 8561 : case t_FRAC:
1939 : {
1940 8561 : GEN a = gel(x,1), b = gel(x,2), c = gammafrac24(a, b, prec);
1941 : long e;
1942 8561 : if (c) return glog(c, prec);
1943 1064 : c = subii(a,b); e = expi(b) - expi(c);
1944 1064 : if (signe(a) < 0)
1945 : { /* gamma will use functional equation x -> z = 1-x = -c/b >= 1/2.
1946 : * lngamma(x) = log |Pi / (sin(Pi z) * Gamma(z))| + I*Pi * floor(x) */
1947 7 : GEN z = mkfrac(negi(c), b), q = ground(z), r = gsub(z,q);
1948 7 : GEN pi = mppi(prec); /* |r| <= 1/2 */
1949 7 : z = fractor(z, prec+EXTRAPREC64);
1950 7 : y = subrr(logr_abs(divrr(pi, mpsin(gmul(pi,r)))), cxgamma(z, 1, prec));
1951 7 : y = gadd(y, mkcomplex(gen_0, mulri(pi, gfloor(x))));
1952 7 : return gc_upto(av, y);
1953 : }
1954 1057 : if (cmpii(shifti(a,1), b) < 0)
1955 : { /* 0 < x < 1/2 gamma would use funeq: adding 1 is cheaper. */
1956 14 : if (expi(a) - expi(b) < -3) /* close to 0 */
1957 : {
1958 14 : if (lg2prec(lgefint(b)) >= prec) x = fractor(x,prec);
1959 14 : y = lngamma1(x, prec);
1960 : }
1961 : else
1962 0 : y = cxgamma(fractor(mkfrac(addii(a,b), b), prec), 1, prec);
1963 14 : return gc_upto(av, gsub(y, glog(x, prec)));
1964 : }
1965 1043 : if (e > 3)
1966 : {
1967 875 : x = mkfrac(c,b);
1968 875 : if (lg2prec(lgefint(b)) >= prec)
1969 7 : x = fractor(x, prec + nbits2extraprec(e));
1970 875 : y = lngamma1(x, prec);
1971 : }
1972 : else
1973 : {
1974 168 : x = fractor(x, e > 1? prec+EXTRAPREC64: prec);
1975 168 : y = cxgamma(x, 1, prec);
1976 : }
1977 1043 : return gc_upto(av, y);
1978 : }
1979 :
1980 22071 : case t_REAL: case t_COMPLEX:
1981 22071 : return cxgamma(x, 1, prec);
1982 :
1983 343 : default:
1984 343 : if (!(y = toser_i(x))) break;
1985 343 : if (lg(y) == 2) pari_err_DOMAIN("lngamma", "argument", "=", gen_0,y);
1986 336 : t = serlngamma0(y,prec);
1987 322 : y0 = simplify_shallow(gel(y,2));
1988 : /* no constant term if y0 = 1 or 2 */
1989 322 : if (!isint(y0,&y0) || signe(y0) <= 0 || abscmpiu(y0,2) > 2)
1990 7 : t = gadd(t, glngamma(y0,prec));
1991 322 : return gc_upto(av, t);
1992 :
1993 14 : case t_PADIC: return gc_upto(av, Qp_lngamma(x));
1994 : }
1995 0 : return trans_eval("lngamma",glngamma,x,prec);
1996 : }
1997 : /********************************************************************/
1998 : /** **/
1999 : /** PSI(x) = GAMMA'(x)/GAMMA(x) **/
2000 : /** **/
2001 : /********************************************************************/
2002 : static void
2003 7 : err_psi(GEN s)
2004 : {
2005 7 : pari_err_DOMAIN("psi","argument", "=",
2006 : strtoGENstr("nonpositive integer"), s);
2007 0 : }
2008 : /* L ~ |log s|^2 */
2009 : static long
2010 4249 : psi_lim(double L, double la, long prec)
2011 : {
2012 4249 : double d = (prec2nbits_mul(prec, 2*M_LN2) - log(L)) / (4*(1+log(la)));
2013 4249 : return (d < 2)? 2: 2 + (long)ceil(d);
2014 : }
2015 : /* max(|log (s + it - Euler)|, 1e-6) */
2016 : static double
2017 4235 : dlogE(double s, double t)
2018 : {
2019 : double rlog, ilog;
2020 4235 : dblclog(s - 0.57721566, t, &rlog,&ilog);
2021 4235 : return maxdd(dblcnorm(rlog,ilog), 1e-6);
2022 : }
2023 : static GEN
2024 4452 : cxpsi(GEN s0, long der, long prec)
2025 : {
2026 : pari_sp av, av2;
2027 : GEN sum, z, a, res, sig, tau, s, unr, s2, sq;
2028 : long lim, nn, k;
2029 4452 : const long la = 3;
2030 4452 : int funeq = 0;
2031 : pari_timer T;
2032 :
2033 4452 : if (der)
2034 : {
2035 203 : av = avma;
2036 203 : res = zetahurwitz(stoi(der + 1), s0, 0, prec2nbits(prec));
2037 203 : if(!odd(der)) res = gneg(res);
2038 203 : return gc_upto(av, gmul(mpfact(der), res));
2039 : }
2040 4249 : if (DEBUGLEVEL>2) timer_start(&T);
2041 4249 : s = trans_fix_arg(&prec,&s0,&sig,&tau,&av,&res);
2042 4249 : if (signe(sig) <= 0) { funeq = 1; s = gsub(gen_1, s); sig = real_i(s); }
2043 4249 : if (typ(s0) == t_INT && signe(s0) <= 0) err_psi(s0);
2044 :
2045 4249 : if (expo(sig) > 300 || (typ(s) == t_COMPLEX && gexpo(gel(s,2)) > 300))
2046 14 : { /* |s| is HUGE. Play safe */
2047 14 : GEN L, S = gprec_w(s,LOWDEFAULTPREC), rS = real_i(S), iS = imag_i(S);
2048 : double l;
2049 14 : lim = psi_lim(rtodbl(gnorm(glog(S,LOWDEFAULTPREC))), la, prec);
2050 14 : l = (2*lim-1)*la / (2.*M_PI);
2051 14 : L = gsub(dbltor(l*l), gsqr(iS));
2052 14 : if (signe(L) < 0) L = gen_0;
2053 14 : L = gsub(gsqrt(L, LOWDEFAULTPREC), rS);
2054 14 : if (signe(L) > 0) nn = (long)ceil(rtodbl(L)); else nn = 1;
2055 : }
2056 : else
2057 : {
2058 4235 : double l, rS = rtodbl(sig), iS = typ(s) == t_REAL? 0.0: rtodbl(imag_i(s));
2059 4235 : lim = psi_lim(dlogE(rS, iS), la, prec);
2060 4235 : l = (2*lim-1)*la / (2.*M_PI);
2061 4235 : l = l*l - iS*iS;
2062 4235 : if (l < 0.) l = 0.;
2063 4235 : nn = (long)ceil( sqrt(l) - rS );
2064 4235 : if (nn < 1) nn = 1;
2065 : }
2066 4249 : if (DEBUGLEVEL>2) err_printf("lim, nn: [%ld, %ld]\n",lim,nn);
2067 4249 : incrprec(prec); unr = real_1(prec); /* one extra word of precision */
2068 4249 : s2 = gmul2n(s, 1); sq = gsqr(s);
2069 4249 : a = gdiv(unr, gaddgs(s, nn)); /* 1 / (s+n) */
2070 4249 : av2 = avma; sum = gmul2n(a, -1);
2071 99610 : for (k = 0; k < nn - 1; k += 2)
2072 : {
2073 95361 : GEN tmp = gaddsg(k*(k + 1), gadd(gmulsg(2*k + 1, s), sq));
2074 95361 : sum = gadd(sum, gdiv(gaddsg(2*k + 1, s2), tmp));
2075 95361 : if ((k & 1023) == 0) sum = gc_upto(av2, sum);
2076 : }
2077 4249 : if (odd(nn)) sum = gadd(sum, gdiv(unr, gaddsg(nn - 1, s)));
2078 4249 : z = gsub(glog(gaddgs(s, nn), prec), sum);
2079 4249 : if (DEBUGLEVEL>2) timer_printf(&T,"sum from 0 to N - 1");
2080 4249 : constbern(lim);
2081 4249 : z = gsub(z, psi_sum(gsqr(a), lim));
2082 4249 : if (DEBUGLEVEL>2) timer_printf(&T,"Bernoulli sum");
2083 4249 : if (funeq)
2084 : {
2085 4004 : GEN pi = mppi(prec);
2086 4004 : z = gadd(z, gmul(pi, gcotan(gmul(pi,s), prec)));
2087 : }
2088 4249 : set_avma(av); return affc_fixlg(z, res);
2089 : }
2090 :
2091 : /* n >= 0; return psi(1+x) + O(x^n), x = pol_x(v) */
2092 : GEN
2093 13685 : psi1series(long n, long v, long prec)
2094 : {
2095 13685 : long i, l = n+3;
2096 13685 : GEN s = cgetg(l, t_SER), z = constzeta(n + 1, prec);
2097 :
2098 13685 : s[1] = evalsigne(1)|evalvalser(0)|evalvarn(v);
2099 71568 : for (i = 1; i <= n+1; i++)
2100 : {
2101 57883 : GEN c = gel(z,i); /* zeta(i) */
2102 57883 : gel(s,i+1) = odd(i)? negr(c): c;
2103 : }
2104 13685 : return s;
2105 : }
2106 : /* T an RgX, return T(X + z0) + O(X^L) */
2107 : static GEN
2108 2374141 : tr(GEN T, GEN z0, long L)
2109 : {
2110 2374141 : GEN s = RgX_to_ser(RgX_Rg_translate(T, z0), L+3);
2111 2374141 : setvarn(s, 0); return s;
2112 : }
2113 : /* z0 a complex number with Re(z0) > 1/2; return psi(z0+x) + O(x^L)
2114 : * using Luke's rational approximation for psi(x) */
2115 : static GEN
2116 10486 : serpsiz0(GEN z0, long L, long v, long prec)
2117 : {
2118 : pari_sp av;
2119 : GEN A,A1,A2, B,B1,B2, Q;
2120 : long n;
2121 10486 : n = gprecision(z0); if (n) prec = n;
2122 10486 : z0 = gtofp(z0, prec + EXTRAPREC64);
2123 : /* Start from n = 3; in Luke's notation, A2 := A_{n-2}, A1 := A_{n-1},
2124 : * A := A_n. Same for B */
2125 10486 : av = avma;
2126 10486 : A2= gdivgu(mkpoln(2, gen_1, utoipos(6)), 2);
2127 10486 : B2 = scalarpol_shallow(utoipos(4), 0);
2128 10486 : A1= gdivgu(mkpoln(3, gen_1, utoipos(82), utoipos(96)), 6);
2129 10486 : B1 = mkpoln(2, utoipos(8), utoipos(28));
2130 10486 : A = gdivgu(mkpoln(4, gen_1, utoipos(387), utoipos(2906), utoipos(1920)), 12);
2131 10486 : B = mkpoln(3, utoipos(14), utoipos(204), utoipos(310));
2132 10486 : A2= tr(A2,z0, L);
2133 10486 : B2= tr(B2,z0, L);
2134 10486 : A1= tr(A1,z0, L);
2135 10486 : B1= tr(B1,z0, L);
2136 10486 : A = tr(A, z0, L);
2137 10486 : B = tr(B, z0, L); Q = gdiv(A, B);
2138 : /* work with z0+x as a variable */
2139 10486 : for (n = 4;; n++)
2140 756427 : {
2141 766913 : GEN Q0 = Q, a, b, r, c3,c2,c1,c0 = muluu(2*n-3, n+1);
2142 766913 : GEN u = subiu(muluu(n, 7*n-9), 6);
2143 766913 : GEN t = addiu(muluu(n, 7*n-19), 4);
2144 : /* c1=(2*n-1)*(3*(n-1)*z+7*n^2-9*n-6);
2145 : * c2=(2*n-3)*(z-n-1)*(-3*(n-1)*z+7*n^2-19*n+4);
2146 : * c3=(2*n-1)*(n-3)*(z-n)*(z-(n+1))*(z+(n-4)); */
2147 766913 : c1 = deg1pol_shallow(muluu(3*(n-1),2*n-1), muliu(u,2*n-1), 0);
2148 766913 : c2 = ZX_mul(deg1pol_shallow(utoipos(2*n-3), negi(muluu(2*n-3,n+1)), 0),
2149 766913 : deg1pol_shallow(utoineg(3*(n-1)), t, 0));
2150 766913 : r = mkvec3(utoipos(n), utoipos(n+1), stoi(4-n));
2151 766913 : c3 = ZX_Z_mul(roots_to_pol(r,0), muluu(2*n-1,n-3));
2152 766913 : c1 = tr(c1, z0, L+3);
2153 766913 : c2 = tr(c2, z0, L+3);
2154 766913 : c3 = tr(c3, z0, L+3);
2155 :
2156 : /* A_{n+1}, B_{n+1} */
2157 766913 : a = gdiv(gadd(gadd(gmul(c1,A),gmul(c2,A1)),gmul(c3,A2)), c0);
2158 766913 : b = gdiv(gadd(gadd(gmul(c1,B),gmul(c2,B1)),gmul(c3,B2)), c0);
2159 766913 : Q = gdiv(a,b);
2160 766913 : if (gexpo(gsub(Q,Q0)) < -prec2nbits(prec)) break;
2161 756427 : A2 = A1; A1 = A; A = a;
2162 756427 : B2 = B1; B1 = B; B = b;
2163 756427 : if (gc_needed(av,1))
2164 : {
2165 0 : if(DEBUGMEM>1) pari_warn(warnmem,"serpsiz0, n = %ld", n);
2166 0 : (void)gc_all(av, 7, &A,&A1,&A2, &B,&B1,&B2, &Q);
2167 : }
2168 : }
2169 10486 : Q = gmul(Q, gmul2n(gsubsg(1, ginv(tr(pol_x(v),z0, L))), 1));
2170 10486 : setvarn(Q, v);
2171 10486 : return gadd(negeuler(prec), Q);
2172 : }
2173 : /* sum (-1)^k*H(m,k)x^k + O(x^L); L > 0;
2174 : * H(m,k) = (-1)^{k * \delta_{m > 0}} sum_{1<=i<m} 1/i^(k+1) */
2175 : static GEN
2176 1575 : Hseries(long m, long L, long v, long prec)
2177 : {
2178 1575 : long i, k, bit, l = L+3, M = m < 0? 1-m: m;
2179 1575 : pari_sp av = avma;
2180 1575 : GEN H = cgetg(l, t_SER);
2181 1575 : H[1] = evalsigne(1)|evalvarn(v)|evalvalser(0);
2182 1575 : prec += EXTRAPREC64;
2183 1575 : bit = -prec2nbits(prec);
2184 8239 : for(k = 2; k < l; k++) gel(H,k) = gen_1; /* i=1 */
2185 1911 : for (i = 2; i < M; i++)
2186 : {
2187 336 : GEN ik = invr(utor(i, prec));
2188 2051 : for (k = 2; k < l; k++)
2189 : {
2190 1715 : if (k > 2) { ik = divru(ik, i); if (expo(ik) < bit) break; }
2191 1715 : gel(H,k) = gadd(gel(H,k), ik);
2192 : }
2193 336 : if (gc_needed(av,3))
2194 : {
2195 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Hseries, i = %ld/%ld", i,M);
2196 0 : H = gc_GEN(av, H);
2197 : }
2198 : }
2199 1575 : if (m > 0)
2200 4613 : for (k = 3; k < l; k+=2) togglesign_safe(&gel(H,k));
2201 1575 : return H;
2202 : }
2203 :
2204 : static GEN
2205 23877 : serpsi(GEN y, long prec)
2206 : {
2207 23877 : GEN Q = NULL, z0, Y = y, Y2;
2208 23877 : long L = lg(y)-2, v = varn(y), vy = valser(y);
2209 :
2210 23877 : if (!L) pari_err_DOMAIN("psi", "argument", "=", gen_0,y);
2211 23870 : if (vy < 0) pari_err_DOMAIN("psi", "series valuation", "<", gen_0,y);
2212 23870 : if (vy)
2213 14 : z0 = gen_0;
2214 : else
2215 : {
2216 23856 : z0 = simplify_shallow(gel(y,2));
2217 23856 : (void)isint(z0, &z0);
2218 : }
2219 23870 : if (typ(z0) == t_INT && !is_bigint(z0))
2220 : {
2221 13384 : long m = itos(z0);
2222 13384 : if (abscmpiu(muluu(prec2nbits(prec),L), labs(m)) > 0)
2223 : { /* psi(m+x) = psi(1+x) + sum_{1 <= i < m} 1/(i+x) for m > 0
2224 : psi(1+x) - sum_{0 <= i < -m} 1/(i+x) for m <= 0 */
2225 13384 : GEN H = NULL;
2226 13384 : if (m <= 0) L--; /* lose series accuracy due to 1/x term */
2227 13384 : if (L)
2228 : {
2229 13377 : Q = psi1series(L, v, prec);
2230 13377 : if (m && m != 1) { H = Hseries(m, L, v, prec); Q = gadd(Q, H); }
2231 13377 : if (m <= 0) Q = gsub(Q, ginv(pol_x(v)));
2232 : }
2233 : else
2234 : {
2235 7 : Q = scalarser(gen_m1, v, 1);
2236 7 : setvalser(Q,-1);
2237 : }
2238 : }
2239 : }
2240 23870 : if (!Q)
2241 : { /* use psi(1-y)=psi(y)+Pi*cotan(Pi*y) ? */
2242 10486 : if (gcmp(real_i(z0),ghalf) < 0) { z0 = gsubsg(1,z0); Y = gsubsg(1,y); }
2243 10486 : Q = serpsiz0(z0, L, v, prec);
2244 : }
2245 23870 : Y2 = serchop0(Y); if (signe(Y2)) Q = gsubst(Q, v, Y2);
2246 : /* psi(z0 + Y2) = psi(Y) */
2247 23870 : if (Y != y)
2248 : { /* psi(y) = psi(Y) + Pi cotan(Pi Y) */
2249 98 : GEN pi = mppi(prec);
2250 98 : if (typ(z0) == t_INT) Y = Y2; /* in this case cotan(Pi*Y2) = cotan(Pi*Y) */
2251 98 : Q = gadd(Q, gmul(pi, gcotan(gmul(pi,Y), prec)));
2252 : }
2253 23870 : return Q;
2254 : }
2255 :
2256 : static ulong
2257 21322 : psi_n(ulong b)
2258 : {
2259 21322 : if (b <= 64) return 50;
2260 21322 : if (b <= 128) return 85;
2261 21322 : if (b <= 192) return 122;
2262 21175 : if (b <= 256) return 150;
2263 12477 : if (b <= 512) return 320;
2264 7 : if (b <= 1024) return 715;
2265 0 : return 0.010709 * pow((double)b, 1.631); /* 1.631 ~ log_3(6) */
2266 : }
2267 :
2268 : GEN
2269 245 : Qp_psi(GEN x, long der)
2270 : {
2271 245 : pari_sp av = avma;
2272 245 : GEN p = padic_p(x), p1 = subis(p,1), z;
2273 245 : long e = valp(x) + precp(x);
2274 245 : if (valp(x) < 0) pari_err_DOMAIN("psi","v_p(x)", "<", gen_0, x);
2275 238 : if (der < 0) pari_err_DOMAIN("psi","der","<", gen_0, stoi(der));
2276 238 : x = cvtop(x, p, e + 1);
2277 238 : z = gmul(mpfact(der), Qp_zetahurwitz(cvtop(stoi(der + 1), p, e + sdivsi(e,p1)), x, -der));
2278 238 : if (!odd(der)) z = gneg(z);
2279 238 : if (!der) z = gadd(mkfrac(p1,p), z);
2280 238 : return gc_upto(av, z);
2281 : }
2282 :
2283 : GEN
2284 49567 : gpsi(GEN x, long prec)
2285 : {
2286 : pari_sp av;
2287 : ulong n;
2288 : GEN y;
2289 49567 : switch(typ(x))
2290 : {
2291 21315 : case t_INT:
2292 21315 : if (signe(x) <= 0) err_psi(x);
2293 21315 : if (lgefint(x) > 3 || (n = itou(x)) > psi_n(prec2nbits(prec))) break;
2294 21315 : av = avma; y = mpeuler(prec);
2295 21315 : return gc_leaf(av, n == 1? negr(y): gsub(harmonic(n-1), y));
2296 4200 : case t_REAL: case t_COMPLEX: return cxpsi(x,0,prec);
2297 0 : case t_PADIC: return Qp_psi(x, 0);
2298 24052 : default:
2299 24052 : av = avma; if (!(y = toser_i(x))) break;
2300 23828 : return gc_upto(av, serpsi(y,prec));
2301 : }
2302 224 : return trans_eval("psi",gpsi,x,prec);
2303 : }
2304 :
2305 : static GEN
2306 84 : _gpsi_der(void *E, GEN x, long prec)
2307 : {
2308 84 : return gpsi_der(x, (long) E, prec);
2309 : }
2310 :
2311 : GEN
2312 784 : gpsi_der(GEN x, long der, long prec)
2313 : {
2314 : pari_sp av;
2315 : ulong n;
2316 : GEN y;
2317 784 : if (der < 0) pari_err_DOMAIN("gpsi", "der", "<", gen_0, stoi(der));
2318 784 : switch(typ(x))
2319 : {
2320 84 : case t_INT:
2321 84 : if (signe(x) <= 0) err_psi(x);
2322 77 : if (lgefint(x) > 3 || (n = itou(x)) > psi_n(prec2nbits(prec))) break;
2323 7 : av = avma;
2324 7 : y = der ? szeta(der + 1, prec): mpeuler(prec);
2325 7 : if (n > 1)
2326 : {
2327 0 : y = gsub(y, harmonic0(n - 1, stoi(der + 1)));
2328 0 : if (!odd(der)) y = gneg(y);
2329 0 : y = gmul(mpfact(der), y);
2330 0 : return gc_leaf(av, y);
2331 : }
2332 252 : case t_REAL: case t_COMPLEX: return cxpsi(x, der, prec);
2333 245 : case t_PADIC: return Qp_psi(x, der);
2334 210 : default:
2335 210 : av = avma; if (!(y = toser_i(x))) break;
2336 196 : if (!der) y = serpsi(y,prec);
2337 : else
2338 : {
2339 147 : y = zetahurwitz(stoi(der + 1), x, 0, prec2nbits(prec));
2340 147 : if(!odd(der)) y = gneg(y);
2341 147 : y = gmul(mpfact(der), y);
2342 : }
2343 189 : return gc_upto(av, y);
2344 : }
2345 84 : return trans_evalgen("psi",(void*)der,_gpsi_der,x,prec);
2346 : }
|
