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