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 251832 : trans_fix_arg(long *prec, GEN *s0, GEN *sig, GEN *tau, pari_sp *av, GEN *res)
28 : {
29 251832 : GEN p1, s = *s0 = cxtoreal(*s0);
30 : long l;
31 251832 : l = precision(s); if (!l) l = *prec;
32 251832 : if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
33 251832 : *res = cgetc(l); *av = avma;
34 251832 : if (typ(s) == t_COMPLEX)
35 : { /* s = sig + i t */
36 207305 : s = cxtofp(s, l+EXTRAPREC64);
37 207305 : *sig = gel(s,1);
38 207305 : *tau = gel(s,2);
39 : }
40 : else /* real number */
41 : {
42 44527 : *sig = s = gtofp(s, l+EXTRAPREC64);
43 44527 : *tau = gen_0;
44 44527 : p1 = trunc2nr(s, 0);
45 44527 : if (!signe(subri(s,p1))) *s0 = p1;
46 : }
47 251832 : *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 4500879 : mpatan(GEN x)
61 : {
62 4500879 : 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 4500879 : if (!sx) return real_0_bit(expo(x));
69 4500844 : l = lp = realprec(x);
70 4500844 : if (absrnz_equal1(x)) { /* |x| = 1 */
71 18254 : y = Pi2n(-2, l+EXTRAPREC64); if (sx < 0) setsigne(y,-1);
72 18254 : return y;
73 : }
74 4482590 : if (l > AGM_ATAN_LIMIT)
75 172 : { av = avma; return gc_leaf(av, atan2_agm(gen_1, x, l)); }
76 :
77 4482418 : e = expo(x); inv = (e >= 0); /* = (|x| > 1 ) */
78 4482418 : if (e > 0) lp += nbits2extraprec(e);
79 :
80 4482418 : y = cgetr(lp); av0 = avma;
81 4482418 : p1 = rtor(x, l+EXTRAPREC64); setabssign(p1); /* p1 = |x| */
82 4482418 : if (inv) p1 = invr(p1);
83 4482418 : e = expo(p1);
84 4482418 : if (e < -100)
85 42965 : alpha = 1.65149612947 - e; /* log_2(Pi) - e */
86 : else
87 4439453 : alpha = log2(M_PI / atan(rtodbl(p1)));
88 4482418 : beta = (double)(prec2nbits(l)>>1);
89 4482418 : delta = 1 + beta - alpha/2;
90 4482418 : if (delta <= 0) { n = 1; m = 0; }
91 : else
92 : {
93 4447533 : double fi = alpha-2;
94 4447533 : if (delta >= fi*fi)
95 : {
96 4209487 : double t = 1 + sqrt(delta);
97 4209487 : n = (long)t;
98 4209487 : m = (long)(t - fi);
99 : }
100 : else
101 : {
102 238046 : n = (long)(1+beta/fi);
103 238046 : m = 0;
104 : }
105 : }
106 4482418 : l2 = l + nbits2extraprec(m);
107 4482418 : p2 = rtor(p1, l2); av = avma;
108 49635609 : for (i=1; i<=m; i++)
109 : {
110 45153191 : p5 = addsr(1, sqrr(p2)); setprec(p5,l2);
111 45153191 : p5 = addsr(1, sqrtr_abs(p5)); setprec(p5,l2);
112 45153191 : affrr(divrr(p2,p5), p2); set_avma(av);
113 : }
114 4482418 : p3 = sqrr(p2); l1 = minss(LOWDEFAULTPREC+EXTRAPREC64, l2); /* l1 increases to l2 */;
115 4482418 : unr = real_1(l2); setprec(unr,l1);
116 4482418 : p4 = cgetr(l2); setprec(p4,l1);
117 4482418 : affrr(divru(unr,2*n+1), p4);
118 4482418 : s = 0; e = expo(p3); av = avma;
119 55459733 : for (i = n; i > 1; i--) /* n >= 1. i = 1 done outside for efficiency */
120 : {
121 50977315 : setprec(p3,l1); p5 = mulrr(p4,p3);
122 50977315 : l1 += nbits2extraprec(dvmdsBIL(s - e, &s)<<TWOPOTBITS_IN_LONG);
123 50977315 : if (l1 > l2) l1 = l2;
124 50977315 : setprec(unr,l1); p5 = subrr(divru(unr,2*i-1), p5);
125 50977315 : setprec(p4,l1); affrr(p5,p4); set_avma(av);
126 : }
127 4482418 : setprec(p3, l2); p5 = mulrr(p4,p3); /* i = 1 */
128 4482418 : setprec(unr,l2); p4 = subrr(unr, p5);
129 :
130 4482418 : p4 = mulrr(p2,p4); shiftr_inplace(p4, m);
131 4482418 : if (inv) p4 = subrr(Pi2n(-1, lp), p4);
132 4482418 : if (sx < 0) togglesign(p4);
133 4482418 : affrr_fixlg(p4,y); set_avma(av0); return y;
134 : }
135 :
136 : GEN
137 59915 : gatan(GEN x, long prec)
138 : {
139 : pari_sp av;
140 : GEN a, y;
141 :
142 59915 : switch(typ(x))
143 : {
144 22668 : 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 738 : mpasin(GEN x)
167 : {
168 738 : pari_sp av = avma;
169 738 : long sx = signe(x);
170 : GEN z, a;
171 :
172 738 : if (!sx) return rcopy(x);
173 738 : 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 738 : a = sqrtr(subsr(1, sqrr(x)));
178 738 : if (realprec(x) > AGM_ATAN_LIMIT)
179 4 : z = atan2_agm(a, x, realprec(x));
180 : else
181 734 : z = mpatan(divrr(x, a));
182 738 : return gc_leaf(av, z);
183 : }
184 :
185 : GEN
186 9033 : gasin(GEN x, long prec)
187 : {
188 : long sx;
189 : pari_sp av;
190 : GEN a, y, p1;
191 :
192 9033 : switch(typ(x))
193 : {
194 1123 : case t_REAL: sx = signe(x);
195 1123 : if (!sx) return rcopy(x);
196 1116 : 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 1088 : 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 4477407 : mparg(GEN x, GEN y)
312 : {
313 4477407 : long prec, sx = signe(x), sy = signe(y);
314 : GEN z;
315 :
316 4477407 : if (!sy)
317 : {
318 0 : if (sx > 0) return real_0_bit(expo(y) - expo(x));
319 0 : return mppi(realprec(x));
320 : }
321 4477407 : prec = realprec(y); if (prec < realprec(x)) prec = realprec(x);
322 4477407 : if (!sx)
323 : {
324 27 : z = Pi2n(-1, prec); if (sy < 0) setsigne(z,-1);
325 27 : return z;
326 : }
327 :
328 4477380 : if (expo(x)-expo(y) > -2)
329 : {
330 3529624 : z = mpatan(divrr(y,x)); if (sx > 0) return z;
331 1478045 : return addrr_sign(z, signe(z), mppi(prec), sy);
332 : }
333 947756 : z = mpatan(divrr(x,y));
334 947756 : return addrr_sign(z, -signe(z), Pi2n(-1, prec), sy);
335 : }
336 :
337 : static GEN
338 8954814 : rfix(GEN x,long prec)
339 : {
340 8954814 : switch(typ(x))
341 : {
342 38029 : case t_INT: return itor(x, prec);
343 642849 : case t_FRAC: return fractor(x, prec);
344 8273936 : case t_REAL: break;
345 0 : default: pari_err_TYPE("rfix (conversion to t_REAL)",x);
346 : }
347 8273936 : return x;
348 : }
349 :
350 : static GEN
351 4477407 : cxarg(GEN x, GEN y, long prec)
352 : {
353 4477407 : pari_sp av = avma;
354 4477407 : x = rfix(x,prec);
355 4477407 : y = rfix(y,prec); return gc_leaf(av, mparg(x,y));
356 : }
357 :
358 : GEN
359 4494935 : garg(GEN x, long prec)
360 : {
361 : long l;
362 4494935 : if (gequal0(x)) pari_err_DOMAIN("arg", "argument", "=", gen_0, x);
363 4494935 : 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 4477407 : case t_COMPLEX:
368 4477407 : l = precision(x); if (l) prec = l;
369 4477407 : 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 6391 : mpsinh(GEN x)
429 : {
430 : pari_sp av;
431 : long lx;
432 : GEN z, res;
433 :
434 6391 : if (!signe(x)) return mpsinh0(expo(x));
435 6391 : lx = realprec(x); res = cgetr(lx); av = avma;
436 6391 : 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 6384 : z = mpexp(x);
446 6384 : z = subrr(z, invr(z));
447 : }
448 6391 : shiftr_inplace(z, -1);
449 6391 : 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 12894 : gsinh(GEN x, long prec)
488 : {
489 : pari_sp av;
490 : GEN y, p1;
491 :
492 12894 : 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 6839 : default:
502 6839 : av = avma; if (!(y = toser_i(x))) break;
503 6811 : if (gequal0(y) && valser(y) == 0) return gc_GEN(av, y);
504 6811 : p1 = gexp(y, prec); p1 = gsub(p1, ginv(p1));
505 6811 : 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 3360 : mpasinh(GEN x)
620 : {
621 3360 : long lx, e, s = signe(x);
622 : GEN z, res;
623 : pari_sp av;
624 :
625 3360 : if (!s) return rcopy(x);
626 3360 : lx = realprec(x); e = expo(x) + BITS_IN_LONG;
627 3360 : res = cgetr(lx);
628 3360 : av = avma;
629 3360 : if (e < 0) x = rtor(x, lx + nbits2extraprec(-e));
630 3360 : z = logr_abs( addrr_sign(x,1, sqrtr_abs( addrs(sqrr(x), 1) ), 1) );
631 3360 : if (signe(x) < 0) togglesign(z);
632 3360 : affrr(z, res); return gc_const(av, res);
633 : }
634 :
635 : GEN
636 40698 : gasinh(GEN x, long prec)
637 : {
638 : pari_sp av;
639 : GEN a, y, p1;
640 :
641 40698 : switch(typ(x))
642 : {
643 3031 : case t_REAL:
644 3031 : if (!signe(x)) return rcopy(x);
645 2996 : return mpasinh(x);
646 :
647 37030 : case t_COMPLEX: {
648 : GEN a, b, d;
649 37030 : if (ismpzero(gel(x,2))) return gasinh(gel(x,1), prec);
650 37030 : av = avma;
651 37030 : if (ismpzero(gel(x,1))) /* avoid cancellation */
652 871 : return gc_GEN(av, mulcxI(gasin(gel(x,2), prec)));
653 36159 : d = gsqrt(gaddsg(1,gsqr(x)), prec); /* Re(d) >= 0 */
654 36159 : a = gadd(d, x);
655 36159 : b = gsub(d, x);
656 : /* avoid cancellation as much as possible */
657 36159 : if (gprecision(a) < gprecision(b))
658 392 : y = gneg(glog(b,prec));
659 : else
660 35767 : y = glog(a,prec);
661 36159 : 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 : GEN x1 = signe(x) > 0? subrs(x,1): addrs(x,1);
695 742 : e = expo(x1);
696 742 : if (e == -(long)HIGHEXPOBIT)
697 0 : return gc_const((pari_sp)(res + lx), real_0_bit(- bit_prec(x) >> 1));
698 742 : if (e < -5) x = rtor(x, realprec(x) + nbits2extraprec(-e));
699 742 : z = logr_abs( addrr_sign(x, 1, sqrtr( subrs(sqrr(x), 1) ), 1) );
700 742 : affrr(z, res); return gc_const(av, res);
701 : }
702 :
703 : GEN
704 7994 : gacosh(GEN x, long prec)
705 : {
706 : pari_sp av;
707 : GEN y;
708 :
709 7994 : switch(typ(x))
710 : {
711 280 : case t_REAL: {
712 280 : long s = signe(x), e = expo(x);
713 : GEN a, b;
714 280 : if (s > 0 && e >= 0) return mpacosh(x);
715 : /* x < 1 */
716 147 : y = cgetg(3,t_COMPLEX); a = gen_0;
717 147 : if (s == 0) b = acos0(e);
718 140 : else if (e < 0) b = mpacos(x); /* -1 < x < 1 */
719 : else {
720 91 : if (!absrnz_equal1(x)) a = mpacosh(x);
721 91 : b = mppi(realprec(x));
722 : }
723 147 : gel(y,1) = a;
724 147 : gel(y,2) = b; return y;
725 : }
726 7413 : case t_COMPLEX: {
727 : GEN a, b, d;
728 7413 : if (ismpzero(gel(x,2))) return gacosh(gel(x,1), prec);
729 7413 : av = avma;
730 7413 : d = gsqrt(gaddsg(-1,gsqr(x)), prec); /* Re(d) >= 0 */
731 7413 : a = gadd(x, d);
732 7413 : b = gsub(x, d);
733 : /* avoid cancellation as much as possible */
734 7413 : if (gprecision(a) < gprecision(b))
735 7 : y = glog(b,prec);
736 : else
737 7406 : y = glog(a,prec);
738 : /* y = \pm log(x + sqrt(x^2-1)) */
739 7413 : if (gsigne(real_i(y)) < 0) y = gneg(y);
740 7413 : return gc_upto(av, y);
741 : }
742 301 : default: {
743 : GEN a, d;
744 : long v;
745 301 : av = avma; if (!(y = toser_i(x))) break;
746 49 : v = valser(y);
747 49 : if (v < 0) pari_err_DOMAIN("acosh","valuation", "<", gen_0, x);
748 42 : if (gequal0(y))
749 : {
750 7 : if (!v) return gc_GEN(av, y);
751 7 : return gc_upto(av, gadd(y, PiI2n(-1, prec)));
752 : }
753 35 : d = gsubgs(gsqr(y),1);
754 35 : if (gequal0(d)) { set_avma(av); return zeroser(varn(y), valser(d)>>1); }
755 21 : d = gdiv(derivser(y), gsqrt(d,prec));
756 21 : a = integser(d);
757 21 : if (v)
758 7 : d = PiI2n(-1, prec); /* I Pi/2 */
759 : else
760 : {
761 14 : d = gel(y,2); if (gequal1(d)) return gc_upto(av,a);
762 7 : d = gacosh(d, prec);
763 : }
764 14 : return gc_upto(av, gadd(d,a));
765 : }
766 : }
767 252 : return trans_eval("acosh",gacosh,x,prec);
768 : }
769 : /********************************************************************/
770 : /** **/
771 : /** AREA HYPERBOLIC TANGENT **/
772 : /** **/
773 : /********************************************************************/
774 :
775 : /* |x| < 1 */
776 : GEN
777 7665 : mpatanh(GEN x)
778 : {
779 7665 : pari_sp av = avma;
780 7665 : long e, s = signe(x);
781 : GEN z;
782 7665 : if (!s) return rcopy(x);
783 7665 : z = s > 0? subsr(1,x): addsr(1,x); e = expo(z);
784 7665 : if (e < -5)
785 : {
786 1218 : x = rtor(x, realprec(x) + nbits2extraprec(-e)-EXTRAPRECWORD);
787 1218 : z = s > 0? subsr(1,x): addsr(1,x); e = expo(z);
788 : }
789 7665 : z = invr(z); shiftr_inplace(z, 1); /* 2/(1-|x|) */
790 7665 : z = logr_abs( addrs(z,-1) ); if (s < 0) togglesign(z);
791 7665 : shiftr_inplace(z, -1); return gc_leaf(av, z);
792 : }
793 :
794 : static long
795 5524155 : get_nmax(double u, double v, long prec)
796 : {
797 5524155 : double d = 2 * log2(((double)v) / u); /* can be 0 due to rounding */
798 5524155 : long nmax = -1;
799 5524155 : if (d)
800 : {
801 5524149 : d = ceil(prec2nbits(prec) / d);
802 5524149 : if (dblexpo(d) < BITS_IN_LONG) nmax = (long)d;
803 : }
804 5524155 : return nmax;
805 : }
806 : /* atanh(u/v) using binary splitting, 0 < u < v */
807 : GEN
808 5524141 : atanhuu(ulong u, ulong v, long prec)
809 : {
810 5524141 : GEN u2 = sqru(u), v2 = sqru(v);
811 5524141 : long i, nmax = get_nmax((double)u, (double)v, prec);
812 : struct abpq_res R;
813 : struct abpq A;
814 5524141 : if (nmax < 0) pari_err_OVERFLOW("atanhuu");
815 5524134 : abpq_init(&A, nmax); /* nmax satisfies (2n+1) (v/u)^2n > 2^bitprec */
816 5524134 : A.a[0] = A.b[0] = gen_1;
817 5524134 : A.p[0] = utoipos(u);
818 5524134 : A.q[0] = utoipos(v);
819 106551216 : for (i = 1; i <= nmax; i++)
820 : {
821 101027082 : A.a[i] = gen_1;
822 101027082 : A.b[i] = utoipos((i<<1)+1);
823 101027082 : A.p[i] = u2;
824 101027082 : A.q[i] = v2;
825 : }
826 5524134 : abpq_sum(&R, 0, nmax, &A);
827 5524134 : return rdivii(R.T, mulii(R.B,R.Q),prec);
828 : }
829 : /* atanh(u/v) using binary splitting, 0 < u < v */
830 : GEN
831 14 : atanhui(ulong u, GEN v, long prec)
832 : {
833 14 : GEN u2 = sqru(u), v2 = sqri(v);
834 14 : long i, nmax = get_nmax((double)u, gtodouble(v), prec);
835 : struct abpq_res R;
836 : struct abpq A;
837 14 : if (nmax < 0) pari_err_OVERFLOW("atanhui");
838 14 : abpq_init(&A, nmax);
839 14 : A.a[0] = A.b[0] = gen_1;
840 14 : A.p[0] = utoipos(u);
841 14 : A.q[0] = v;
842 35 : for (i = 1; i <= nmax; i++)
843 : {
844 21 : A.a[i] = gen_1;
845 21 : A.b[i] = utoipos((i<<1)+1);
846 21 : A.p[i] = u2;
847 21 : A.q[i] = v2;
848 : }
849 14 : abpq_sum(&R, 0, nmax, &A);
850 14 : return rdivii(R.T, mulii(R.B,R.Q),prec);
851 : }
852 :
853 : static void
854 28 : err_atanh(GEN x, GEN bad) { pari_err_DOMAIN("atanh", "x", "=", bad, x); }
855 :
856 : GEN
857 70121 : gatanh(GEN x, long prec)
858 : {
859 : long sx;
860 : pari_sp av;
861 : GEN a, y, z;
862 :
863 70121 : switch(typ(x))
864 : {
865 126 : case t_INT:
866 126 : sx = signe(x);
867 126 : if (!sx) return real_0(prec);
868 119 : z = cgetg(3, t_COMPLEX); av = avma;
869 119 : if (lgefint(x) == 3)
870 : {
871 112 : if (x[2] == 1) err_atanh(x, sx == 1? gen_1: gen_m1);
872 84 : a = atanhuu(1, x[2], prec);
873 : }
874 : else
875 7 : a = atanhui(1, x, prec);
876 91 : gel(z,1) = gc_leaf(av, a);
877 91 : gel(z,2) = Pi2n(-1, prec);
878 91 : togglesign(sx > 0? gel(z,2): gel(z,1));
879 91 : return z;
880 350 : case t_FRAC:
881 : {
882 : long ly, lz, e;
883 :
884 350 : y = gel(x,1); ly = lgefint(y);
885 350 : z = gel(x,2); lz = lgefint(z); if (ly > 3 && lz > 3) break;
886 350 : if (abscmpii(y, z) > 0) /* |y| > z; lz = 3 */
887 : {
888 252 : ulong u = z[2];
889 252 : av = avma; e = expi((signe(y) < 0)? addii(y, z): subii(y, z));
890 252 : set_avma(av); if (e < - prec2nbits(prec)) break;
891 252 : z = cgetg(3, t_COMPLEX); av = avma;
892 252 : a = ly == 3? atanhuu(u, y[2], prec): atanhui(u, y, prec);
893 252 : gel(z,1) = gc_leaf(av, a);
894 252 : gel(z,2) = Pi2n(-1, prec);
895 252 : togglesign(signe(y) > 0? gel(z,2): gel(z,1));
896 : }
897 : else
898 : { /* |y| < z; ly = 3 */
899 98 : av = avma; e = expi((signe(y) < 0)? addii(y, z): subii(y, z));
900 98 : set_avma(av); if (e < - prec2nbits(prec)) break;
901 98 : a = lz == 3? atanhuu(y[2], z[2], prec): atanhui(y[2], z, prec);
902 91 : z = gc_leaf(av, a);
903 91 : if (signe(y) < 0) togglesign(z);
904 : }
905 343 : return z;
906 : }
907 21831 : case t_REAL:
908 21831 : sx = signe(x);
909 21831 : if (!sx) return rcopy(x);
910 21803 : if (expo(x) < 0) return mpatanh(x);
911 :
912 14138 : y = cgetg(3,t_COMPLEX);
913 14138 : av = avma;
914 14138 : z = subrs(x,1);
915 14138 : if (!signe(z)) err_atanh(x, gen_1);
916 14138 : z = invr(z); shiftr_inplace(z, 1); /* 2/(x-1)*/
917 14138 : z = addrs(z,1);
918 14138 : if (!signe(z)) err_atanh(x, gen_m1);
919 14138 : z = logr_abs(z);
920 14138 : shiftr_inplace(z, -1); /* (1/2)log((1+x)/(x-1)) */
921 14138 : gel(y,1) = gc_leaf(av, z);
922 14138 : gel(y,2) = Pi2n(-1, realprec(x));
923 14138 : if (sx > 0) togglesign(gel(y,2));
924 14138 : return y;
925 :
926 47779 : case t_COMPLEX: /* 2/(1-z) - 1 = (1+z) / (1-z) */
927 47779 : if (ismpzero(gel(x,2))) return gatanh(gel(x,1), prec);
928 38898 : av = avma; z = glog( gaddgs(gdivsg(2,gsubsg(1,x)),-1), prec );
929 38898 : return gc_upto(av, gmul2n(z,-1));
930 :
931 35 : default:
932 35 : av = avma; if (!(y = toser_i(x))) break;
933 28 : if (valser(y) < 0) pari_err_DOMAIN("atanh","valuation", "<", gen_0, x);
934 21 : z = gdiv(derivser(y), gsubsg(1,gsqr(y)));
935 14 : a = integser(z);
936 14 : if (!valser(y)) a = gadd(a, gatanh(gel(y,2),prec));
937 14 : return gc_upto(av, a);
938 : }
939 7 : return trans_eval("atanh",gatanh,x,prec);
940 : }
941 : /********************************************************************/
942 : /** **/
943 : /** EULER'S GAMMA **/
944 : /** **/
945 : /********************************************************************/
946 : /* 0 < a < b */
947 : static GEN
948 27453 : mulu_interval_step_i(ulong a, ulong b, ulong step)
949 : {
950 : ulong k, l, N, n;
951 : long lx;
952 : GEN x;
953 :
954 27453 : n = 1 + (b-a) / step;
955 27453 : b -= (b-a) % step;
956 : /* step | b-a */
957 27453 : lx = 1; x = cgetg(2 + n/2, t_VEC);
958 27453 : N = b + a;
959 27453 : for (k = a;; k += step)
960 : {
961 180393 : l = N - k; if (l <= k) break;
962 152940 : gel(x,lx++) = muluu(k,l);
963 : }
964 27453 : if (l == k) gel(x,lx++) = utoipos(k);
965 27453 : setlg(x, lx); return x;
966 : }
967 : static GEN
968 150033 : _mul(void *data, GEN x, GEN y)
969 : {
970 150033 : long prec = (long)data;
971 : /* switch to t_REAL ? */
972 150033 : if (typ(x) == t_INT && lg2prec(lgefint(x)) > prec) x = itor(x, prec);
973 150033 : if (typ(y) == t_INT && lg2prec(lgefint(y)) > prec) y = itor(y, prec);
974 150033 : return mpmul(x, y);
975 : }
976 : static GEN
977 27453 : mulu_interval_step_prec(long l, long m, long s, long prec)
978 : {
979 27453 : GEN v = mulu_interval_step_i(l, m, s);
980 27453 : return gen_product(v, (void*)prec, &_mul);
981 : }
982 :
983 : /* x * (i*(i+1)) */
984 : static GEN
985 7829479 : muliunextu(GEN x, ulong i)
986 : {
987 7829479 : if (i & HIGHMASK) /* i(i+1) >= 2^BITS_IN_LONG*/
988 0 : return mulii(x, muluu(i, i+1));
989 : else
990 7829479 : return muliu(x, i*(i+1));
991 : }
992 : /* arg(s + it); in principle s + it != 0 but currently arg(0. + i*0.) := 0 */
993 : double
994 234185 : dblcarg(double s, double t)
995 : {
996 : double x;
997 234185 : if (!t) return (s >= 0)? 0.: M_PI;
998 202580 : if (!s) return (t >= 0)? M_PI/2: -M_PI/2;
999 202573 : x = atan(t/s);
1000 202573 : return (s >= 0)? x: (t >= 0)? x + M_PI : x - M_PI;
1001 : }
1002 :
1003 : /* Let z = s + it, set a = Re(log z), b = Im(log z) */
1004 : void
1005 234185 : dblclog(double s, double t, double *a, double *b)
1006 : {
1007 234185 : *a = log(dblcnorm(s, t)) / 2;
1008 234185 : *b = dblcarg(s, t);
1009 234185 : }
1010 :
1011 : double
1012 16716 : dblcabs(double s, double t) { return sqrt(dblcnorm(s, t)); }
1013 : double
1014 271670 : dblcnorm(double s, double t) { return s*s + t*t; }
1015 :
1016 : #if 0
1017 : /* x, z t_REAL. Compute unique x in ]-z,z] congruent to x mod 2z */
1018 : static GEN
1019 : red_mod_2z(GEN x, GEN z)
1020 : {
1021 : GEN Z = gmul2n(z, 1), d = subrr(z, x);
1022 : /* require little accuracy */
1023 : if (!signe(d)) return x;
1024 : setprec(d, nbits2prec(expo(d) - expo(Z)));
1025 : return addrr(mulir(floorr(divrr(d, Z)), Z), x);
1026 : }
1027 : #endif
1028 :
1029 : static GEN
1030 11081 : negeuler(long prec) { GEN g = mpeuler(prec); setsigne(g, -1); return g; }
1031 : /* lngamma(1+z) = -Euler*z + sum_{i > 1} zeta(i)/i (-z)^i
1032 : * at relative precision prec, |z| <= 1/2 is small */
1033 : static GEN
1034 15686 : lngamma1(GEN z, long prec)
1035 : { /* sum_{i > l} |z|^(i-1) = |z|^l / (1-|z|) < 2^-B
1036 : * for l > (B+1) / |log2(|z|)| */
1037 15686 : long i, l = ceil((prec2nbits(prec) + 1) / - dbllog2(z));
1038 : GEN s, vz;
1039 :
1040 15686 : if (l <= 1) return gmul(negeuler(prec), z);
1041 15511 : vz = constzeta(l, prec);
1042 1057384 : for (i = l, s = gen_0; i > 0; i--)
1043 : {
1044 1041873 : GEN c = divru(gel(vz,i), i);
1045 1041873 : if (odd(i)) setsigne(c, -1);
1046 1041873 : s = gadd(gmul(s,z), c);
1047 : }
1048 15511 : return gmul(z, s);
1049 : }
1050 : /* B_i / (i(i-1)), i even. Sometimes NOT reduced (but gadd/gmul won't care)!*/
1051 : static GEN
1052 7829479 : bern_unextu(long i)
1053 7829479 : { GEN B = bernfrac(i); return mkfrac(gel(B,1), muliunextu(gel(B,2), i-1)); }
1054 : /* B_i / i, i even. Sometimes NOT reduced (but gadd/gmul won't care)!*/
1055 : static GEN
1056 211323 : bern_u(long i)
1057 211323 : { GEN B = bernfrac(i); return mkfrac(gel(B,1), muliu(gel(B,2), i)); }
1058 : /* sum_{i > 0} B_{2i}/(2i(2i-1)) * a^(i-1) */
1059 : static GEN
1060 214086 : lngamma_sum(GEN a, long N)
1061 : {
1062 214086 : pari_sp av = avma;
1063 214086 : GEN S = bern_unextu(2*N);
1064 : long i;
1065 7829479 : for (i = 2*N-2; i > 0; i -= 2)
1066 : {
1067 7615393 : S = gadd(bern_unextu(i), gmul(a,S));
1068 7615393 : if (gc_needed(av,3))
1069 : {
1070 0 : if(DEBUGMEM>1) pari_warn(warnmem,"gamma: i = %ld", i);
1071 0 : S = gc_upto(av, S);
1072 : }
1073 : }
1074 214086 : return S;
1075 : }
1076 : /* sum_{i > 0} B_{2i}/(2i) * a^i */
1077 : static GEN
1078 4249 : psi_sum(GEN a, long N)
1079 : {
1080 4249 : pari_sp av = avma;
1081 4249 : GEN S = bern_u(2*N);
1082 : long i;
1083 211323 : for (i = 2*N-2; i > 0; i -= 2)
1084 : {
1085 207074 : S = gadd(bern_u(i), gmul(a,S));
1086 207074 : if (gc_needed(av,3))
1087 : {
1088 0 : if(DEBUGMEM>1) pari_warn(warnmem,"psi: i = %ld", i);
1089 0 : S = gc_upto(av, S);
1090 : }
1091 : }
1092 4249 : return gmul(a,S);
1093 : }
1094 : static void
1095 226177 : gamma_optim(double ssig, double st, long prec, long *plim, long *pN)
1096 : {
1097 : double la, l,l2,u,v, rlogs, ilogs;
1098 226177 : long N = 1, lim;
1099 226177 : dblclog(ssig,st, &rlogs,&ilogs);
1100 : /* Re (s - 1/2) log(s) */
1101 226177 : u = (ssig - 0.5)*rlogs - st * ilogs;
1102 : /* Im (s - 1/2) log(s) */
1103 226177 : v = (ssig - 0.5)*ilogs + st * rlogs;
1104 : /* l2 = | (s - 1/2) log(s) - s + log(2Pi)/2 |^2 ~ |lngamma(s))|^2 */
1105 226177 : u = u - ssig + log(2.*M_PI)/2;
1106 226177 : v = v - st;
1107 226177 : l2 = u*u + v*v;
1108 226177 : if (l2 < 0.000001) l2 = 0.000001;
1109 226177 : l = (prec2nbits_mul(prec, M_LN2) - log(l2)/2) / 2.;
1110 226177 : if (l < 0) l = 0.;
1111 :
1112 226177 : if (st > 1 && l > 0)
1113 67347 : {
1114 67347 : double t = st * M_PI / l;
1115 67347 : la = t * log(t);
1116 67347 : if (la < 4.) la = 4.;
1117 67347 : if (la > 150) la = t;
1118 : }
1119 : else
1120 158830 : la = 4.; /* heuristic */
1121 226177 : lim = (long)ceil(l / (1.+ log(la)));
1122 226177 : if (lim == 0) lim = 1;
1123 :
1124 226177 : u = (lim-0.5) * la / M_PI;
1125 226177 : l2 = u*u - st*st;
1126 226177 : if (l2 > 0)
1127 : {
1128 213338 : double t = ceil(sqrt(l2) - ssig);
1129 213338 : if (t > 1) N = (long)t;
1130 : }
1131 226177 : *plim = lim; *pN = N;
1132 226177 : }
1133 : /* do we use lngamma1 instead of Euler-Maclaurin ? */
1134 : static int
1135 228970 : gamma_use_1(double s, double t, long prec, long *plim, long *pN)
1136 : {
1137 228970 : double a = s-1, d = fabs(a) + fabs(t);
1138 : long k;
1139 228970 : if (d < 1e-16) return 1;
1140 226177 : gamma_optim(s, t, prec, plim, pN);
1141 226177 : if (d >= 0.5) return 0;
1142 16534 : k = prec2nbits(prec) / -log2(dblcnorm(a, t)); /* 2k = lngamma1 bound */
1143 16534 : return (t ? k: 1.5*k) < *plim + *pN;
1144 : }
1145 : static GEN
1146 229005 : cxgamma(GEN s0, int dolog, long prec)
1147 : {
1148 : GEN s, a, y, res, sig, tau, B, nnx, pi, pi2;
1149 229005 : long i, esig, et, lim, N = 1;
1150 : pari_sp av, av2;
1151 229005 : int funeq = 0;
1152 : pari_timer T;
1153 :
1154 229005 : if (DEBUGLEVEL>5) timer_start(&T);
1155 229005 : s = trans_fix_arg(&prec,&s0,&sig,&tau,&av,&res);
1156 :
1157 229005 : esig = expo(sig);
1158 229005 : et = signe(tau)? expo(tau): 0;
1159 229005 : if ((signe(sig) <= 0 || esig < -1) && et <= 16)
1160 : { /* s <--> 1-s */
1161 21770 : funeq = 1; s = gsubsg(1, s); sig = real_i(s);
1162 : }
1163 :
1164 : /* find "optimal" parameters [lim, N] */
1165 229005 : if (esig > 300 || et > 300)
1166 35 : { /* |s| is HUGE ! Play safe and avoid inf / NaN */
1167 : GEN S, iS, l2, la, u;
1168 : double logla, l;
1169 :
1170 35 : S = gprec_w(s,LOWDEFAULTPREC);
1171 : /* l2 ~ |lngamma(s))|^2 */
1172 35 : l2 = gnorm(gmul(S, glog(S, LOWDEFAULTPREC)));
1173 35 : l = (prec2nbits_mul(prec, M_LN2) - rtodbl(glog(l2,LOWDEFAULTPREC))/2) / 2.;
1174 35 : if (l < 0) l = 0.;
1175 :
1176 35 : iS = imag_i(S);
1177 35 : if (et > 0 && l > 0)
1178 21 : {
1179 21 : GEN t = gmul(iS, dbltor(M_PI / l)), logt = glog(t,LOWDEFAULTPREC);
1180 21 : la = gmul(t, logt);
1181 21 : if (gcmpgs(la, 3) < 0) { logla = log(3.); la = stoi(3); }
1182 14 : else if (gcmpgs(la, 150) > 0) { logla = rtodbl(logt); la = t; }
1183 7 : else logla = rtodbl(mplog(la));
1184 : }
1185 : else
1186 : {
1187 14 : logla = log(3.); la = stoi(3);
1188 : }
1189 35 : lim = (long)ceil(l / (1.+ logla));
1190 35 : if (lim == 0) lim = 1;
1191 :
1192 35 : u = gmul(la, dbltor((lim-0.5)/M_PI));
1193 35 : l2 = gsub(gsqr(u), gsqr(iS));
1194 35 : if (signe(l2) > 0)
1195 : {
1196 14 : l2 = gsub(gsqrt(l2,3), sig);
1197 14 : if (signe(l2) > 0) N = itos( gceil(l2) );
1198 : }
1199 : }
1200 : else
1201 : { /* |s| is moderate. Use floats */
1202 228970 : double ssig = rtodbl(sig);
1203 228970 : double st = typ(s) == t_REAL? 0.0: rtodbl(imag_i(s));
1204 :
1205 228970 : if (gamma_use_1(ssig, st, prec, &lim, &N))
1206 : { /* s ~ 1: loggamma(1+u) ~ - Euler * u, cancellation */
1207 14919 : if (funeq) /* s0 ~ 0: use lngamma(s0)+log(s0) = lngamma(s0+1) */
1208 119 : y = dolog? gsub(lngamma1(s0,prec), glog(s0,prec))
1209 119 : : gdiv(gexp(lngamma1(s0,prec), prec), s0);
1210 : else
1211 : {
1212 14800 : if (isint1(s0))
1213 : {
1214 1683 : set_avma(av);
1215 1683 : return dolog? real_0(prec): real_1(prec);
1216 : }
1217 13117 : y = lngamma1(gsubgs(s0,1),prec);
1218 13117 : if (!dolog) y = gexp(y,prec);
1219 : }
1220 13236 : set_avma(av); return affc_fixlg(y, res);
1221 : }
1222 : }
1223 214086 : if (DEBUGLEVEL>5) err_printf("lim, N: [%ld, %ld]\n",lim,N);
1224 214086 : incrprec(prec);
1225 :
1226 214086 : av2 = avma;
1227 214086 : y = s;
1228 214086 : if (typ(s0) == t_INT)
1229 : {
1230 2605 : ulong ss = itou_or_0(s0);
1231 2605 : if (signe(s0) <= 0)
1232 0 : pari_err_DOMAIN("gamma","argument", "=",
1233 : strtoGENstr("nonpositive integer"), s0);
1234 2605 : if (!ss || ss + (ulong)N < ss) {
1235 7 : for (i=1; i < N; i++)
1236 : {
1237 0 : y = mulri(y, addiu(s0, i));
1238 0 : if (gc_needed(av2,3))
1239 : {
1240 0 : if(DEBUGMEM>1) pari_warn(warnmem,"gamma");
1241 0 : y = gc_leaf(av2, y);
1242 : }
1243 : }
1244 : } else {
1245 34000 : for (i=1; i < N; i++)
1246 : {
1247 31402 : y = mulru(y, ss + i);
1248 31402 : if (gc_needed(av2,3))
1249 : {
1250 0 : if(DEBUGMEM>1) pari_warn(warnmem,"gamma");
1251 0 : y = gc_leaf(av2, y);
1252 : }
1253 : }
1254 : }
1255 : }
1256 : else
1257 : { /* Compute lngamma mod 2 I Pi */
1258 211481 : GEN sq = gsqr(s);
1259 211481 : pari_sp av3 = avma;
1260 4319628 : for (i = 1; i < N - 1; i += 2)
1261 : {
1262 4108147 : y = gmul(y, gaddsg(i*(i + 1), gadd(gmulsg(2*i + 1, s), sq)));
1263 4108147 : if (gc_needed(av2,3))
1264 : {
1265 0 : if(DEBUGMEM>1) pari_warn(warnmem,"gamma");
1266 0 : y = gc_upto(av3, y);
1267 : }
1268 : }
1269 211481 : if (!odd(N)) y = gmul(y, gaddsg(N - 1, s));
1270 : }
1271 214086 : if (DEBUGLEVEL>5) timer_printf(&T,"product from 0 to N-1");
1272 214086 : constbern(lim);
1273 214086 : nnx = gaddgs(s, N); a = ginv(nnx);
1274 214086 : B = gadd(gsub(gmul(gsub(nnx, ghalf), glog(nnx,prec)), nnx),
1275 : gmul(a, lngamma_sum(gsqr(a), lim)));
1276 214086 : if (DEBUGLEVEL>5) timer_printf(&T,"Bernoulli sum");
1277 :
1278 214086 : pi = mppi(prec); pi2 = shiftr(pi, 1);
1279 214086 : if (dolog)
1280 : {
1281 15764 : if (typ(s) == t_REAL)
1282 : {
1283 12411 : if (!funeq) y = logr_abs(divrr(sqrtr(pi2), y));
1284 : else
1285 : {
1286 7 : GEN T = shiftr(sqrtr(pi2),-1); /* sqrt(Pi/2) */
1287 : /* s0 < 0, step (*) simplifies: imag(lngamma(s0)) = - Pi * floor(s0) */
1288 7 : y = logr_abs(divrr(mulrr(y, T), mpsin(gmul(pi,s0))));
1289 7 : y = mkcomplex(y, mulri(pi, gfloor(s0)));
1290 7 : B = gneg(B);
1291 : }
1292 : }
1293 : else
1294 : { /* log(y), fixing imaginary part */
1295 3353 : long prec2 = LOWDEFAULTPREC;
1296 3353 : GEN k, s2 = gprec_w(s, prec2), y2 = garg(s2, prec2); /* ~ Im log(s) */
1297 10438 : for (i=1; i < N; i++) y2 = gadd(y2, garg(gaddgs(s2,i), prec2));
1298 3353 : y = glog(y, prec);
1299 3353 : k = ground( gdiv(gsub(y2, imag_i(y)), Pi2n(1,prec2)) );
1300 3353 : if (signe(k)) y = gadd(y, mulcxI(mulir(k, Pi2n(1, prec))));
1301 3353 : if (!funeq) y = gsub(shiftr(logr_abs(pi2),-1), y); /* y -> sqrt(2Pi)/y */
1302 : else
1303 : { /* recall that s = 1 - s0 */
1304 287 : GEN T = shiftr(sqrtr(pi2),-1); /* sqrt(Pi/2) */
1305 : /* (*) Compute log(sin(Pi s0)) so that it has branch cuts along
1306 : * (-oo, 0] and [1, oo). To do this in a numerically stable way
1307 : * we must compute the log first then mangle its imaginary part.
1308 : * The rounding operation below is stable because we're rounding
1309 : * a number which is already within 1/4 of an integer. */
1310 :
1311 : /* z = log(sin(Pi s0) / sqrt(Pi/2)) */
1312 287 : GEN z = glog(gdiv(gsin(gmul(pi,s0),prec), T), prec);
1313 287 : GEN b = shiftr(subrs(shiftr(sig, 1), 1), -2); /* (2 Re(s)-1) / 4 */
1314 287 : y = gsub(y, z);
1315 287 : if (gsigne(imag_i(s)) > 0) togglesign(b);
1316 287 : z = roundr(gsub(gdiv(imag_i(z), pi2), b)); /* round( Im(z)/2Pi - b ) */
1317 287 : if (signe(z)) { /* y += I*z, z a t_REAL */
1318 0 : z = mulir(z, pi2);
1319 0 : if (typ(y) == t_COMPLEX) gel(y,2) = gadd(gel(y,2), z);
1320 0 : else y = mkcomplex(y, z);
1321 : }
1322 287 : B = gneg(B);
1323 : }
1324 : }
1325 15764 : y = gadd(B, y);
1326 : }
1327 : else
1328 : {
1329 198322 : GEN sqrtpi2 = sqrtr(pi2);
1330 198322 : if (funeq)
1331 : { /* y --> y Pi/(sin(Pi s) * sqrt(2Pi)) = y sqrt(Pi/2)/sin(Pi s) */
1332 21357 : y = gdiv(gmul(shiftr(sqrtpi2,-1),y), gsin(gmul(pi,s0), prec));
1333 : /* don't use s above: sin(pi s0) = sin(pi s) and the former is
1334 : * more accurate, esp. if s0 ~ 0 */
1335 21357 : B = gneg(B);
1336 : }
1337 : else /* y --> sqrt(2Pi) / y */
1338 176965 : y = gdiv(sqrtpi2, y);
1339 198322 : y = gmul(gexp(B, prec), y);
1340 : }
1341 214086 : set_avma(av); return affc_fixlg(y, res);
1342 : }
1343 :
1344 : /* Theory says n > C * b^1.5 / log(b). Timings:
1345 : * b = 64*[1, 2, 3, 4, 5, 6, 7, 10, 20, 30, 50, 100, 200, 500];
1346 : * n = [1450, 1930, 2750, 3400, 4070, 5000, 6000, 8800, 26000, 50000, 130000,
1347 : * 380000, 1300000, 6000000]; */
1348 : static long
1349 37751 : gamma2_n(long prec)
1350 : {
1351 37751 : long b = prec2nbits(prec);
1352 37751 : if (b <= 64) return 1450;
1353 37128 : if (b <= 128) return 1930;
1354 31143 : if (b <= 192) return 2750;
1355 16935 : if (b <= 256) return 3400;
1356 7341 : if (b <= 320) return 4070;
1357 6925 : if (b <= 384) return 5000;
1358 4115 : if (b <= 448) return 6000;
1359 3933 : return 10.0 * b * sqrt(b) / log(b);
1360 : }
1361 :
1362 : /* m even, Gamma((m+1) / 2) */
1363 : static GEN
1364 37751 : gammahs(long m, long prec)
1365 : {
1366 37751 : GEN y = cgetr(prec), z;
1367 37751 : pari_sp av = avma;
1368 37751 : long ma = labs(m);
1369 :
1370 37751 : if (ma > gamma2_n(prec))
1371 : {
1372 0 : z = stor(m + 1, prec); shiftr_inplace(z, -1);
1373 0 : affrr(cxgamma(z,0,prec), y);
1374 0 : set_avma(av); return y;
1375 : }
1376 37751 : z = sqrtr( mppi(prec) );
1377 37751 : if (m)
1378 : {
1379 23002 : GEN t = mulu_interval_step_prec(1, ma-1, 2, prec + EXTRAPREC64);
1380 23002 : if (typ(t) == t_INT)
1381 22995 : z = m >= 0? mulri(z, t): divri(z, t);
1382 : else
1383 7 : z = m >= 0? mulrr(z, t): divrr(z, t);
1384 23002 : if (m < 0 && (m&3) == 2) setsigne(z,-1);
1385 23002 : shiftr_inplace(z, -m/2);
1386 : }
1387 37751 : 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 215746 : gel(u, 1) = s;
1435 215746 : gel(u, 2) = s;
1436 891798 : for (j = 2; j < p; ++j)
1437 676052 : 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 10345546 : for (j = 1; j < p; ++j)
1443 8250464 : 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 2095082 : 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 447739 : 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 215746 : 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 28427 : serlngamma0(GEN y, long prec)
1551 : {
1552 : GEN t;
1553 28427 : if (valser(y)) pari_err_DOMAIN("lngamma","valuation", "!=", gen_0, y);
1554 28420 : t = derivser(y);
1555 : /* don't compute psi if y'=0 */
1556 28420 : if (signe(t)) t = gmul(t, gpsi(y,prec));
1557 28420 : return integser(t);
1558 : }
1559 :
1560 : static GEN
1561 28392 : sergamma(GEN y, long prec)
1562 : {
1563 : GEN z, y0, Y;
1564 28392 : if (lg(y) == 2) pari_err_DOMAIN("gamma", "argument", "=", gen_0,y);
1565 : /* exp(lngamma) */
1566 28385 : if (valser(y) > 0) return gdiv(gexp(glngamma(gaddgs(y,1),prec),prec),y);
1567 28098 : y0 = simplify_shallow(gel(y,2));
1568 28098 : z = NULL; Y = y;
1569 28098 : if (isint(y0, &y0))
1570 : { /* fun eq. avoids log singularity of lngamma at negative ints */
1571 13797 : long s = signe(y0);
1572 : /* possible if y[2] is an inexact 0 */
1573 13797 : if (!s) return gdiv(gexp(glngamma(gaddgs(y,1),prec),prec),y);
1574 13790 : if (signe(y0) < 0) { Y = gsubsg(1, y); y0 = subsi(1, y0); }
1575 13790 : if (abscmpiu(y0, 50) < 0) z = mpfact(itos(y0)-1); /* more precise */
1576 : }
1577 28091 : if (!z) z = ggamma(y0,prec);
1578 28091 : z = gmul(z, gexp(serlngamma0(Y,prec),prec));
1579 28091 : 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 28091 : 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; return a t_REAL */
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 133 : 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 28 : 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 46074 : gammafrac24(GEN a, GEN b, long prec)
1745 : {
1746 : pari_sp av;
1747 : long A, B, m, am, x, bit;
1748 : GEN z0, z, t;
1749 46074 : if (!(A = itos_or_0(a)) || !(B = itos_or_0(b)) || B > 24) return NULL;
1750 44709 : switch(B)
1751 : {
1752 37730 : 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 : am = labs(m); 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 (am > 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 : GEN Bm = rpowuu(B, am, prec); /* B^|m|, t_REAL */
1767 3808 : if (m > 0)
1768 : {
1769 3759 : t = mulu_interval_step(x, (m-1)*B + x, B);
1770 3759 : t = divmpr(t, Bm);
1771 : }
1772 : else
1773 : {
1774 49 : t = mulu_interval_step(B-x, -m*B - x, B);
1775 49 : t = divrmp(Bm, t);
1776 : }
1777 3808 : if (m < 0 && odd(m)) togglesign(t);
1778 3808 : z = mulrr(z,t);
1779 : }
1780 5628 : affrr(z, z0); set_avma(av); return z0;
1781 : }
1782 1337 : return NULL;
1783 : }
1784 : GEN
1785 635643 : ggamma(GEN x, long prec)
1786 : {
1787 : pari_sp av;
1788 : GEN y;
1789 :
1790 635643 : switch(typ(x))
1791 : {
1792 364337 : case t_INT:
1793 364337 : if (signe(x) <= 0)
1794 0 : pari_err_DOMAIN("gamma","argument", "=",
1795 : strtoGENstr("nonpositive integer"), x);
1796 364337 : return mpfactr(itos(x) - 1, prec);
1797 :
1798 205317 : case t_REAL: case t_COMPLEX:
1799 205317 : return cxgamma(x, 0, prec);
1800 :
1801 37513 : case t_FRAC:
1802 : {
1803 37513 : GEN a = gel(x,1), b = gel(x,2), c = gammafrac24(a, b, prec);
1804 37513 : if (c) return c;
1805 1652 : av = avma; c = subii(a,b);
1806 1652 : if (signe(a) < 0)
1807 : { /* gamma will use functional equation x -> z = 1-x = -c/b >= 1/2.
1808 : * Gamma(x) = Pi / (sin(Pi z) * Gamma(z)) */
1809 49 : GEN z = mkfrac(negi(c), b), q = ground(z), r = gsub(z,q);
1810 49 : GEN pi = mppi(prec); /* |r| <= 1/2 */
1811 49 : z = fractor(z, prec+EXTRAPREC64);
1812 49 : y = divrr(pi, mulrr(mpsin(gmul(pi, r)), cxgamma(z, 0, prec)));
1813 49 : if (mpodd(q)) togglesign(y);
1814 49 : return gc_upto(av, y);
1815 : }
1816 1603 : if (cmpii(shifti(a,1), b) < 0)
1817 : { /* 0 < x < 1/2 gamma would use funeq: adding 1 is cheaper. */
1818 343 : if (expi(a) - expi(b) < -3) /* close to 0 */
1819 : {
1820 14 : if (lg2prec(lgefint(b)) >= prec) x = fractor(x,prec);
1821 14 : y = mpexp(lngamma1(x, prec));
1822 : }
1823 : else
1824 329 : y = cxgamma(fractor(mkfrac(addii(a,b), b), prec), 0, prec);
1825 343 : return gc_upto(av, gdiv(y, x));
1826 : }
1827 1260 : if (expi(c) - expi(b) < -3)
1828 : { /* x = 1 + c/b is close to 1 */
1829 336 : x = mkfrac(c,b);
1830 336 : if (lg2prec(lgefint(b)) >= prec) x = fractor(x,prec);
1831 336 : y = mpexp(lngamma1(x, prec));
1832 : }
1833 : else
1834 924 : y = cxgamma(fractor(x, prec), 0, prec);
1835 1260 : return gc_upto(av, y);
1836 : }
1837 :
1838 84 : case t_PADIC: return Qp_gamma(x);
1839 28392 : default:
1840 28392 : av = avma; if (!(y = toser_i(x))) break;
1841 28392 : return gc_upto(av, sergamma(y, prec));
1842 : }
1843 0 : return trans_eval("gamma",ggamma,x,prec);
1844 : }
1845 :
1846 : static GEN
1847 524 : mpfactr_basecase(long n, long prec)
1848 : {
1849 524 : GEN v = cgetg(expu(n) + 2, t_VEC);
1850 524 : long k, prec2 = prec + EXTRAPREC64;
1851 : GEN a;
1852 524 : for (k = 1;; k++)
1853 4451 : {
1854 4975 : long m = n >> (k-1), l;
1855 4975 : if (m <= 2) break;
1856 4451 : l = (1 + (n >> k)) | 1;
1857 : /* product of odd numbers in ]n / 2^k, 2 / 2^(k-1)] */
1858 4451 : a = mulu_interval_step_prec(l, m, 2, prec2);
1859 4451 : gel(v,k) = k == 1? a: gpowgs(a, k);
1860 : }
1861 4451 : a = gel(v,--k); while (--k) a = mpmul(a, gel(v,k));
1862 524 : if (typ(a) == t_INT) a = itor(a, prec); else a = gprec_wtrunc(a, prec);
1863 524 : shiftr_inplace(a, factorial_lval(n, 2));
1864 524 : return a;
1865 : }
1866 : /* Theory says n > C * b^1.5 / log(b). Timings:
1867 : * b = [64, 128, 192, 256, 512, 1024, 2048, 4096, 8192, 16384]
1868 : * n = [1930, 2650, 3300, 4270, 9000, 23000, 75000, 210000, 750000, 2400000] */
1869 : static long
1870 643 : mpfactr_n(long prec)
1871 : {
1872 643 : long b = prec2nbits(prec);
1873 643 : if (b <= 64) return 1930;
1874 97 : if (b <= 128) return 2650;
1875 76 : if (b <= 192) return 3300;
1876 76 : return b * sqrt(b);
1877 : }
1878 : static GEN
1879 7889 : mpfactr_small(long n, long prec)
1880 : {
1881 7889 : GEN f = cgetr(prec);
1882 7889 : pari_sp av = avma;
1883 7889 : if (n < 410)
1884 7889 : affir(mpfact(n), f);
1885 : else
1886 0 : affrr(mpfactr_basecase(n, prec), f);
1887 7889 : set_avma(av); return f;
1888 : }
1889 : GEN
1890 407955 : mpfactr(long n, long prec)
1891 : {
1892 407955 : GEN f = cgetr(prec);
1893 407955 : pari_sp av = avma;
1894 :
1895 407955 : if (n < 410)
1896 407312 : affir(mpfact(n), f);
1897 : else
1898 : {
1899 643 : long N = mpfactr_n(prec);
1900 524 : GEN z = n <= N? mpfactr_basecase(n, prec)
1901 643 : : cxgamma(utor(n+1, prec), 0, prec);
1902 643 : affrr(z, f);
1903 : }
1904 407955 : set_avma(av); return f;
1905 : }
1906 :
1907 : /* First a little worse than mpfactr_n because of the extra logarithm.
1908 : * Asymptotically same. */
1909 : static ulong
1910 7889 : lngamma_n(long prec)
1911 : {
1912 7889 : long b = prec2nbits(prec);
1913 : double N;
1914 7889 : if (b <= 64) return 1450UL;
1915 7889 : if (b <= 128) return 2010UL;
1916 308 : if (b <= 192) return 2870UL;
1917 308 : N = b * sqrt(b);
1918 308 : if (b <= 256) return N/1.25;
1919 0 : if (b <= 512) return N/1.2;
1920 0 : if (b <= 2048) return N/1.1;
1921 0 : return N;
1922 : }
1923 :
1924 : GEN
1925 38899 : glngamma(GEN x, long prec)
1926 : {
1927 38899 : pari_sp av = avma;
1928 : GEN y, y0, t;
1929 :
1930 38899 : switch(typ(x))
1931 : {
1932 7896 : case t_INT:
1933 : {
1934 : ulong n;
1935 7896 : if (signe(x) <= 0)
1936 0 : pari_err_DOMAIN("lngamma","argument", "=",
1937 : strtoGENstr("nonpositive integer"), x);
1938 7896 : n = itou_or_0(x);
1939 7896 : if (!n || n > lngamma_n(prec)) return cxgamma(x, 1, prec);
1940 7889 : return gc_leaf(av, logr_abs( mpfactr_small(n-1, prec) ));
1941 : }
1942 8561 : case t_FRAC:
1943 : {
1944 8561 : GEN a = gel(x,1), b = gel(x,2), c = gammafrac24(a, b, prec);
1945 : long e;
1946 8561 : if (c) return glog(c, prec);
1947 1064 : c = subii(a,b); e = expi(b) - expi(c);
1948 1064 : if (signe(a) < 0)
1949 : { /* gamma will use functional equation x -> z = 1-x = -c/b >= 1/2.
1950 : * lngamma(x) = log |Pi / (sin(Pi z) * Gamma(z))| + I*Pi * floor(x) */
1951 7 : GEN z = mkfrac(negi(c), b), q = ground(z), r = gsub(z,q);
1952 7 : GEN pi = mppi(prec); /* |r| <= 1/2 */
1953 7 : z = fractor(z, prec+EXTRAPREC64);
1954 7 : y = subrr(logr_abs(divrr(pi, mpsin(gmul(pi,r)))), cxgamma(z, 1, prec));
1955 7 : y = gadd(y, mkcomplex(gen_0, mulri(pi, gfloor(x))));
1956 7 : return gc_upto(av, y);
1957 : }
1958 1057 : if (cmpii(shifti(a,1), b) < 0)
1959 : { /* 0 < x < 1/2 gamma would use funeq: adding 1 is cheaper. */
1960 14 : if (expi(a) - expi(b) < -3) /* close to 0 */
1961 : {
1962 14 : if (lg2prec(lgefint(b)) >= prec) x = fractor(x,prec);
1963 14 : y = lngamma1(x, prec);
1964 : }
1965 : else
1966 0 : y = cxgamma(fractor(mkfrac(addii(a,b), b), prec), 1, prec);
1967 14 : return gc_upto(av, gsub(y, glog(x, prec)));
1968 : }
1969 1043 : if (e > 3)
1970 : {
1971 875 : x = mkfrac(c,b);
1972 875 : if (lg2prec(lgefint(b)) >= prec)
1973 7 : x = fractor(x, prec + nbits2extraprec(e));
1974 875 : y = lngamma1(x, prec);
1975 : }
1976 : else
1977 : {
1978 168 : x = fractor(x, e > 1? prec+EXTRAPREC64: prec);
1979 168 : y = cxgamma(x, 1, prec);
1980 : }
1981 1043 : return gc_upto(av, y);
1982 : }
1983 :
1984 22085 : case t_REAL: case t_COMPLEX:
1985 22085 : return cxgamma(x, 1, prec);
1986 :
1987 343 : default:
1988 343 : if (!(y = toser_i(x))) break;
1989 343 : if (lg(y) == 2) pari_err_DOMAIN("lngamma", "argument", "=", gen_0,y);
1990 336 : t = serlngamma0(y,prec);
1991 322 : y0 = simplify_shallow(gel(y,2));
1992 : /* no constant term if y0 = 1 or 2 */
1993 322 : if (!isint(y0,&y0) || signe(y0) <= 0 || abscmpiu(y0,2) > 2)
1994 7 : t = gadd(t, glngamma(y0,prec));
1995 322 : return gc_upto(av, t);
1996 :
1997 14 : case t_PADIC: return gc_upto(av, Qp_lngamma(x));
1998 : }
1999 0 : return trans_eval("lngamma",glngamma,x,prec);
2000 : }
2001 : /********************************************************************/
2002 : /** **/
2003 : /** PSI(x) = GAMMA'(x)/GAMMA(x) **/
2004 : /** **/
2005 : /********************************************************************/
2006 : static void
2007 7 : err_psi(GEN s)
2008 : {
2009 7 : pari_err_DOMAIN("psi","argument", "=",
2010 : strtoGENstr("nonpositive integer"), s);
2011 0 : }
2012 : /* L ~ |log s|^2 */
2013 : static long
2014 4249 : psi_lim(double L, double la, long prec)
2015 : {
2016 4249 : double d = (prec2nbits_mul(prec, 2*M_LN2) - log(L)) / (4*(1+log(la)));
2017 4249 : return (d < 2)? 2: 2 + (long)ceil(d);
2018 : }
2019 : /* max(|log (s + it - Euler)|, 1e-6) */
2020 : static double
2021 4235 : dlogE(double s, double t)
2022 : {
2023 : double rlog, ilog;
2024 4235 : dblclog(s - 0.57721566, t, &rlog,&ilog);
2025 4235 : return maxdd(dblcnorm(rlog,ilog), 1e-6);
2026 : }
2027 : static GEN
2028 4452 : cxpsi(GEN s0, long der, long prec)
2029 : {
2030 : pari_sp av, av2;
2031 : GEN sum, z, a, res, sig, tau, s, unr, s2, sq;
2032 : long lim, nn, k;
2033 4452 : const long la = 3;
2034 4452 : int funeq = 0;
2035 : pari_timer T;
2036 :
2037 4452 : if (der)
2038 : {
2039 203 : av = avma;
2040 203 : res = zetahurwitz(stoi(der + 1), s0, 0, prec2nbits(prec));
2041 203 : if(!odd(der)) res = gneg(res);
2042 203 : return gc_upto(av, gmul(mpfact(der), res));
2043 : }
2044 4249 : if (DEBUGLEVEL>2) timer_start(&T);
2045 4249 : s = trans_fix_arg(&prec,&s0,&sig,&tau,&av,&res);
2046 4249 : if (signe(sig) <= 0) { funeq = 1; s = gsub(gen_1, s); sig = real_i(s); }
2047 4249 : if (typ(s0) == t_INT && signe(s0) <= 0) err_psi(s0);
2048 :
2049 4249 : if (expo(sig) > 300 || (typ(s) == t_COMPLEX && gexpo(gel(s,2)) > 300))
2050 14 : { /* |s| is HUGE. Play safe */
2051 14 : GEN L, S = gprec_w(s,LOWDEFAULTPREC), rS = real_i(S), iS = imag_i(S);
2052 : double l;
2053 14 : lim = psi_lim(rtodbl(gnorm(glog(S,LOWDEFAULTPREC))), la, prec);
2054 14 : l = (2*lim-1)*la / (2.*M_PI);
2055 14 : L = gsub(dbltor(l*l), gsqr(iS));
2056 14 : if (signe(L) < 0) L = gen_0;
2057 14 : L = gsub(gsqrt(L, LOWDEFAULTPREC), rS);
2058 14 : if (signe(L) > 0) nn = (long)ceil(rtodbl(L)); else nn = 1;
2059 : }
2060 : else
2061 : {
2062 4235 : double l, rS = rtodbl(sig), iS = typ(s) == t_REAL? 0.0: rtodbl(imag_i(s));
2063 4235 : lim = psi_lim(dlogE(rS, iS), la, prec);
2064 4235 : l = (2*lim-1)*la / (2.*M_PI);
2065 4235 : l = l*l - iS*iS;
2066 4235 : if (l < 0.) l = 0.;
2067 4235 : nn = (long)ceil( sqrt(l) - rS );
2068 4235 : if (nn < 1) nn = 1;
2069 : }
2070 4249 : if (DEBUGLEVEL>2) err_printf("lim, nn: [%ld, %ld]\n",lim,nn);
2071 4249 : incrprec(prec); unr = real_1(prec); /* one extra word of precision */
2072 4249 : s2 = gmul2n(s, 1); sq = gsqr(s);
2073 4249 : a = gdiv(unr, gaddgs(s, nn)); /* 1 / (s+n) */
2074 4249 : av2 = avma; sum = gmul2n(a, -1);
2075 99610 : for (k = 0; k < nn - 1; k += 2)
2076 : {
2077 95361 : GEN tmp = gaddsg(k*(k + 1), gadd(gmulsg(2*k + 1, s), sq));
2078 95361 : sum = gadd(sum, gdiv(gaddsg(2*k + 1, s2), tmp));
2079 95361 : if ((k & 1023) == 0) sum = gc_upto(av2, sum);
2080 : }
2081 4249 : if (odd(nn)) sum = gadd(sum, gdiv(unr, gaddsg(nn - 1, s)));
2082 4249 : z = gsub(glog(gaddgs(s, nn), prec), sum);
2083 4249 : if (DEBUGLEVEL>2) timer_printf(&T,"sum from 0 to N - 1");
2084 4249 : constbern(lim);
2085 4249 : z = gsub(z, psi_sum(gsqr(a), lim));
2086 4249 : if (DEBUGLEVEL>2) timer_printf(&T,"Bernoulli sum");
2087 4249 : if (funeq)
2088 : {
2089 4004 : GEN pi = mppi(prec);
2090 4004 : z = gadd(z, gmul(pi, gcotan(gmul(pi,s), prec)));
2091 : }
2092 4249 : set_avma(av); return affc_fixlg(z, res);
2093 : }
2094 :
2095 : /* n >= 0; return psi(1+x) + O(x^n), x = pol_x(v) */
2096 : GEN
2097 13958 : psi1series(long n, long v, long prec)
2098 : {
2099 13958 : long i, l = n+3;
2100 13958 : GEN s = cgetg(l, t_SER), z = constzeta(n + 1, prec);
2101 :
2102 13958 : s[1] = evalsigne(1)|evalvalser(0)|evalvarn(v);
2103 72716 : for (i = 1; i <= n+1; i++)
2104 : {
2105 58758 : GEN c = gel(z,i); /* zeta(i) */
2106 58758 : gel(s,i+1) = odd(i)? negr(c): c;
2107 : }
2108 13958 : return s;
2109 : }
2110 : /* T an RgX, return T(X + z0) + O(X^L) */
2111 : static GEN
2112 2474896 : tr(GEN T, GEN z0, long L)
2113 : {
2114 2474896 : GEN s = RgX_to_ser(RgX_Rg_translate(T, z0), L+3);
2115 2474896 : setvarn(s, 0); return s;
2116 : }
2117 : /* z0 a complex number with Re(z0) > 1/2; return psi(z0+x) + O(x^L)
2118 : * using Luke's rational approximation for psi(x) */
2119 : static GEN
2120 10906 : serpsiz0(GEN z0, long L, long v, long prec)
2121 : {
2122 : pari_sp av;
2123 : GEN A,A1,A2, B,B1,B2, Q;
2124 : long n;
2125 10906 : n = gprecision(z0); if (n) prec = n;
2126 10906 : z0 = gtofp(z0, prec + EXTRAPREC64);
2127 : /* Start from n = 3; in Luke's notation, A2 := A_{n-2}, A1 := A_{n-1},
2128 : * A := A_n. Same for B */
2129 10906 : av = avma;
2130 10906 : A2= gdivgu(mkpoln(2, gen_1, utoipos(6)), 2);
2131 10906 : B2 = scalarpol_shallow(utoipos(4), 0);
2132 10906 : A1= gdivgu(mkpoln(3, gen_1, utoipos(82), utoipos(96)), 6);
2133 10906 : B1 = mkpoln(2, utoipos(8), utoipos(28));
2134 10906 : A = gdivgu(mkpoln(4, gen_1, utoipos(387), utoipos(2906), utoipos(1920)), 12);
2135 10906 : B = mkpoln(3, utoipos(14), utoipos(204), utoipos(310));
2136 10906 : A2= tr(A2,z0, L);
2137 10906 : B2= tr(B2,z0, L);
2138 10906 : A1= tr(A1,z0, L);
2139 10906 : B1= tr(B1,z0, L);
2140 10906 : A = tr(A, z0, L);
2141 10906 : B = tr(B, z0, L); Q = gdiv(A, B);
2142 : /* work with z0+x as a variable */
2143 10906 : for (n = 4;; n++)
2144 788612 : {
2145 799518 : GEN Q0 = Q, a, b, r, c3,c2,c1,c0 = muluu(2*n-3, n+1);
2146 799518 : GEN u = subiu(muluu(n, 7*n-9), 6);
2147 799518 : GEN t = addiu(muluu(n, 7*n-19), 4);
2148 : /* c1=(2*n-1)*(3*(n-1)*z+7*n^2-9*n-6);
2149 : * c2=(2*n-3)*(z-n-1)*(-3*(n-1)*z+7*n^2-19*n+4);
2150 : * c3=(2*n-1)*(n-3)*(z-n)*(z-(n+1))*(z+(n-4)); */
2151 799518 : c1 = deg1pol_shallow(muluu(3*(n-1),2*n-1), muliu(u,2*n-1), 0);
2152 799518 : c2 = ZX_mul(deg1pol_shallow(utoipos(2*n-3), negi(muluu(2*n-3,n+1)), 0),
2153 799518 : deg1pol_shallow(utoineg(3*(n-1)), t, 0));
2154 799518 : r = mkvec3(utoipos(n), utoipos(n+1), stoi(4-n));
2155 799518 : c3 = ZX_Z_mul(roots_to_pol(r,0), muluu(2*n-1,n-3));
2156 799518 : c1 = tr(c1, z0, L+3);
2157 799518 : c2 = tr(c2, z0, L+3);
2158 799518 : c3 = tr(c3, z0, L+3);
2159 :
2160 : /* A_{n+1}, B_{n+1} */
2161 799518 : a = gdiv(gadd(gadd(gmul(c1,A),gmul(c2,A1)),gmul(c3,A2)), c0);
2162 799518 : b = gdiv(gadd(gadd(gmul(c1,B),gmul(c2,B1)),gmul(c3,B2)), c0);
2163 799518 : Q = gdiv(a,b);
2164 799518 : if (gexpo(gsub(Q,Q0)) < -prec2nbits(prec)) break;
2165 788612 : A2 = A1; A1 = A; A = a;
2166 788612 : B2 = B1; B1 = B; B = b;
2167 788612 : if (gc_needed(av,1))
2168 : {
2169 0 : if(DEBUGMEM>1) pari_warn(warnmem,"serpsiz0, n = %ld", n);
2170 0 : (void)gc_all(av, 7, &A,&A1,&A2, &B,&B1,&B2, &Q);
2171 : }
2172 : }
2173 10906 : Q = gmul(Q, gmul2n(gsubsg(1, ginv(tr(pol_x(v),z0, L))), 1));
2174 10906 : setvarn(Q, v);
2175 10906 : return gadd(negeuler(prec), Q);
2176 : }
2177 : /* sum (-1)^k*H(m,k)x^k + O(x^L); L > 0;
2178 : * H(m,k) = (-1)^{k * \delta_{m > 0}} sum_{1<=i<m} 1/i^(k+1) */
2179 : static GEN
2180 1575 : Hseries(long m, long L, long v, long prec)
2181 : {
2182 1575 : long i, k, bit, l = L+3, M = m < 0? 1-m: m;
2183 1575 : pari_sp av = avma;
2184 1575 : GEN H = cgetg(l, t_SER);
2185 1575 : H[1] = evalsigne(1)|evalvarn(v)|evalvalser(0);
2186 1575 : prec += EXTRAPREC64;
2187 1575 : bit = -prec2nbits(prec);
2188 8239 : for(k = 2; k < l; k++) gel(H,k) = gen_1; /* i=1 */
2189 1911 : for (i = 2; i < M; i++)
2190 : {
2191 336 : GEN ik = invr(utor(i, prec));
2192 2051 : for (k = 2; k < l; k++)
2193 : {
2194 1715 : if (k > 2) { ik = divru(ik, i); if (expo(ik) < bit) break; }
2195 1715 : gel(H,k) = gadd(gel(H,k), ik);
2196 : }
2197 336 : if (gc_needed(av,3))
2198 : {
2199 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Hseries, i = %ld/%ld", i,M);
2200 0 : H = gc_GEN(av, H);
2201 : }
2202 : }
2203 1575 : if (m > 0)
2204 4613 : for (k = 3; k < l; k+=2) togglesign_safe(&gel(H,k));
2205 1575 : return H;
2206 : }
2207 :
2208 : static GEN
2209 24570 : serpsi(GEN y, long prec)
2210 : {
2211 24570 : GEN Q = NULL, z0, Y = y, Y2;
2212 24570 : long L = lg(y)-2, v = varn(y), vy = valser(y);
2213 :
2214 24570 : if (!L) pari_err_DOMAIN("psi", "argument", "=", gen_0,y);
2215 24563 : if (vy < 0) pari_err_DOMAIN("psi", "series valuation", "<", gen_0,y);
2216 24563 : if (vy)
2217 21 : z0 = gen_0;
2218 : else
2219 : {
2220 24542 : z0 = simplify_shallow(gel(y,2));
2221 24542 : (void)isint(z0, &z0);
2222 : }
2223 24563 : if (typ(z0) == t_INT && !is_bigint(z0))
2224 : {
2225 13657 : long m = itos(z0);
2226 13657 : if (abscmpiu(muluu(prec2nbits(prec),L), labs(m)) > 0)
2227 : { /* psi(m+x) = psi(1+x) + sum_{1 <= i < m} 1/(i+x) for m > 0
2228 : psi(1+x) - sum_{0 <= i < -m} 1/(i+x) for m <= 0 */
2229 13657 : GEN H = NULL;
2230 13657 : if (m <= 0) L--; /* lose series accuracy due to 1/x term */
2231 13657 : if (L)
2232 : {
2233 13650 : Q = psi1series(L, v, prec);
2234 13650 : if (m && m != 1) { H = Hseries(m, L, v, prec); Q = gadd(Q, H); }
2235 13650 : if (m <= 0) Q = gsub(Q, ginv(pol_x(v)));
2236 : }
2237 : else
2238 : {
2239 7 : Q = scalarser(gen_m1, v, 1);
2240 7 : setvalser(Q,-1);
2241 : }
2242 : }
2243 : }
2244 24563 : if (!Q)
2245 : { /* use psi(1-y)=psi(y)+Pi*cotan(Pi*y) ? */
2246 10906 : if (gcmp(real_i(z0),ghalf) < 0) { z0 = gsubsg(1,z0); Y = gsubsg(1,y); }
2247 10906 : Q = serpsiz0(z0, L, v, prec);
2248 : }
2249 24563 : Y2 = valser(Y)==0 ? serchop0(Y): Y; Q = gsubst(Q, v, Y2);
2250 : /* psi(z0 + Y2) = psi(Y) */
2251 24563 : if (Y != y)
2252 : { /* psi(y) = psi(Y) + Pi cotan(Pi Y) */
2253 98 : GEN pi = mppi(prec);
2254 98 : if (typ(z0) == t_INT) Y = Y2; /* in this case cotan(Pi*Y2) = cotan(Pi*Y) */
2255 98 : Q = gadd(Q, gmul(pi, gcotan(gmul(pi,Y), prec)));
2256 : }
2257 24563 : return Q;
2258 : }
2259 :
2260 : static ulong
2261 21322 : psi_n(ulong b)
2262 : {
2263 21322 : if (b <= 64) return 50;
2264 21322 : if (b <= 128) return 85;
2265 21322 : if (b <= 192) return 122;
2266 21175 : if (b <= 256) return 150;
2267 12477 : if (b <= 512) return 320;
2268 7 : if (b <= 1024) return 715;
2269 0 : return 0.010709 * pow((double)b, 1.631); /* 1.631 ~ log_3(6) */
2270 : }
2271 :
2272 : GEN
2273 245 : Qp_psi(GEN x, long der)
2274 : {
2275 245 : pari_sp av = avma;
2276 245 : GEN p = padic_p(x), p1 = subis(p,1), z;
2277 245 : long e = valp(x) + precp(x);
2278 245 : if (valp(x) < 0) pari_err_DOMAIN("psi","v_p(x)", "<", gen_0, x);
2279 238 : if (der < 0) pari_err_DOMAIN("psi","der","<", gen_0, stoi(der));
2280 238 : x = cvtop(x, p, e + 1);
2281 238 : z = gmul(mpfact(der), Qp_zetahurwitz(cvtop(stoi(der + 1), p, e + sdivsi(e,p1)), x, -der));
2282 238 : if (!odd(der)) z = gneg(z);
2283 238 : if (!der) z = gadd(mkfrac(p1,p), z);
2284 238 : return gc_upto(av, z);
2285 : }
2286 :
2287 : GEN
2288 50239 : gpsi(GEN x, long prec)
2289 : {
2290 : pari_sp av;
2291 : ulong n;
2292 : GEN y;
2293 50239 : switch(typ(x))
2294 : {
2295 21315 : case t_INT:
2296 21315 : if (signe(x) <= 0) err_psi(x);
2297 21315 : if (lgefint(x) > 3 || (n = itou(x)) > psi_n(prec2nbits(prec))) break;
2298 21315 : av = avma; y = mpeuler(prec);
2299 21315 : return gc_leaf(av, n == 1? negr(y): gsub(harmonic(n-1), y));
2300 4200 : case t_REAL: case t_COMPLEX: return cxpsi(x,0,prec);
2301 0 : case t_PADIC: return Qp_psi(x, 0);
2302 24724 : default:
2303 24724 : av = avma; if (!(y = toser_i(x))) break;
2304 24500 : return gc_upto(av, serpsi(y,prec));
2305 : }
2306 224 : return trans_eval("psi",gpsi,x,prec);
2307 : }
2308 :
2309 : static GEN
2310 84 : _gpsi_der(void *E, GEN x, long prec)
2311 : {
2312 84 : return gpsi_der(x, (long) E, prec);
2313 : }
2314 :
2315 : GEN
2316 805 : gpsi_der(GEN x, long der, long prec)
2317 : {
2318 : pari_sp av;
2319 : ulong n;
2320 : GEN y;
2321 805 : if (der < 0) pari_err_DOMAIN("gpsi", "der", "<", gen_0, stoi(der));
2322 805 : switch(typ(x))
2323 : {
2324 84 : case t_INT:
2325 84 : if (signe(x) <= 0) err_psi(x);
2326 77 : if (lgefint(x) > 3 || (n = itou(x)) > psi_n(prec2nbits(prec))) break;
2327 7 : av = avma;
2328 7 : y = der ? szeta(der + 1, prec): mpeuler(prec);
2329 7 : if (n > 1)
2330 : {
2331 0 : y = gsub(y, harmonic0(n - 1, stoi(der + 1)));
2332 0 : if (!odd(der)) y = gneg(y);
2333 0 : y = gmul(mpfact(der), y);
2334 0 : return gc_leaf(av, y);
2335 : }
2336 252 : case t_REAL: case t_COMPLEX: return cxpsi(x, der, prec);
2337 245 : case t_PADIC: return Qp_psi(x, der);
2338 231 : default:
2339 231 : av = avma; if (!(y = toser_i(x))) break;
2340 217 : if (!der) y = serpsi(y,prec);
2341 : else
2342 : {
2343 147 : y = zetahurwitz(stoi(der + 1), x, 0, prec2nbits(prec));
2344 147 : if(!odd(der)) y = gneg(y);
2345 147 : y = gmul(mpfact(der), y);
2346 : }
2347 210 : return gc_upto(av, y);
2348 : }
2349 84 : return trans_evalgen("psi",(void*)der,_gpsi_der,x,prec);
2350 : }
|
