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