1 /*
2 Copyright (C) 2015 Fredrik Johansson
3
4 This file is part of Arb.
5
6 Arb is free software: you can redistribute it and/or modify it under
7 the terms of the GNU Lesser General Public License (LGPL) as published
8 by the Free Software Foundation; either version 2.1 of the License, or
9 (at your option) any later version. See <http://www.gnu.org/licenses/>.
10 */
11
12 #include "acb_hypgeom.h"
13
14 void
acb_hypgeom_ei_asymp(acb_t res,const acb_t z,slong prec)15 acb_hypgeom_ei_asymp(acb_t res, const acb_t z, slong prec)
16 {
17 acb_t t, u;
18
19 acb_init(t);
20 acb_init(u);
21
22 acb_one(t);
23 acb_neg(u, z);
24
25 acb_hypgeom_u_asymp(u, t, t, u, -1, prec);
26 acb_div(u, u, z, prec);
27 acb_exp(t, z, prec);
28 acb_mul(u, u, t, prec);
29
30 if (arb_is_zero(acb_imagref(z)))
31 {
32 arb_zero(acb_imagref(u));
33 }
34 else if (arb_is_positive(acb_imagref(z)))
35 {
36 acb_const_pi(t, prec);
37 arb_add(acb_imagref(u), acb_imagref(u), acb_realref(t), prec);
38 }
39 else if (arb_is_negative(acb_imagref(z)))
40 {
41 acb_const_pi(t, prec);
42 arb_sub(acb_imagref(u), acb_imagref(u), acb_realref(t), prec);
43 }
44 else
45 {
46 /* add [-pi,pi] i */
47 acb_const_pi(t, prec);
48 arb_add_error(acb_imagref(u), acb_realref(t));
49 }
50
51 acb_swap(res, u);
52
53 acb_clear(t);
54 acb_clear(u);
55 }
56
57 /*
58 Ei(z) = z 2F2(1,1,2,2,z) + 0.5[log(z)-log(1/z)] + gamma
59 */
60 void
acb_hypgeom_ei_2f2(acb_t res,const acb_t z,slong prec)61 acb_hypgeom_ei_2f2(acb_t res, const acb_t z, slong prec)
62 {
63 acb_t a, t;
64 acb_struct b[2];
65
66 acb_init(a);
67 acb_init(b);
68 acb_init(b + 1);
69 acb_init(t);
70
71 acb_one(a);
72 acb_set_ui(b, 2);
73 acb_set_ui(b + 1, 2);
74
75 acb_hypgeom_pfq_direct(t, a, 1, b, 2, z, -1, prec);
76 acb_mul(t, t, z, prec);
77
78 arb_const_euler(acb_realref(a), prec);
79 arb_add(acb_realref(t), acb_realref(t), acb_realref(a), prec);
80
81 if (arb_is_zero(acb_imagref(z)))
82 {
83 if (arb_is_positive(acb_realref(z)))
84 {
85 acb_log(a, z, prec);
86 }
87 else
88 {
89 /* ok if overlapping zero -- will be indeterminate either way */
90 acb_neg(a, z);
91 acb_log(a, a, prec);
92 arb_zero(acb_imagref(a));
93 }
94
95 acb_add(t, t, a, prec);
96 }
97 else if (arb_is_nonzero(acb_imagref(z))
98 || arb_is_nonnegative(acb_realref(z))) /* not overlapping (-inf,0] */
99 {
100 acb_log(a, z, prec);
101 acb_add(t, t, a, prec);
102 }
103 else
104 {
105 acb_log(a, z, prec);
106 arb_zero(acb_imagref(a));
107
108 acb_const_pi(b, prec);
109 arb_add_error(acb_imagref(a), acb_realref(b));
110
111 acb_add(t, t, a, prec);
112 }
113
114 acb_swap(res, t);
115
116 acb_clear(a);
117 acb_clear(b);
118 acb_clear(b + 1);
119 acb_clear(t);
120 }
121
122 void
acb_hypgeom_ei(acb_t res,const acb_t z,slong prec)123 acb_hypgeom_ei(acb_t res, const acb_t z, slong prec)
124 {
125 if (acb_hypgeom_u_use_asymp(z, prec))
126 acb_hypgeom_ei_asymp(res, z, prec);
127 else
128 acb_hypgeom_ei_2f2(res, z, prec);
129 }
130
131