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