1 /* shichi.c
2 *
3 * Hyperbolic sine and cosine integrals
4 *
5 *
6 *
7 * SYNOPSIS:
8 *
9 * double x, Chi, Shi, shichi();
10 *
11 * shichi( x, &Chi, &Shi );
12 *
13 *
14 * DESCRIPTION:
15 *
16 * Approximates the integrals
17 *
18 * x
19 * -
20 * | | cosh t - 1
21 * Chi(x) = eul + ln x + | ----------- dt,
22 * | | t
23 * -
24 * 0
25 *
26 * x
27 * -
28 * | | sinh t
29 * Shi(x) = | ------ dt
30 * | | t
31 * -
32 * 0
33 *
34 * where eul = 0.57721566490153286061 is Euler's constant.
35 * The integrals are evaluated by power series for x < 8
36 * and by Chebyshev expansions for x between 8 and 88.
37 * For large x, both functions approach exp(x)/2x.
38 * Arguments greater than 88 in magnitude return NPY_INFINITY.
39 *
40 *
41 * ACCURACY:
42 *
43 * Test interval 0 to 88.
44 * Relative error:
45 * arithmetic function # trials peak rms
46 * IEEE Shi 30000 6.9e-16 1.6e-16
47 * Absolute error, except relative when |Chi| > 1:
48 * IEEE Chi 30000 8.4e-16 1.4e-16
49 */
50
51 /*
52 * Cephes Math Library Release 2.0: April, 1987
53 * Copyright 1984, 1987 by Stephen L. Moshier
54 * Direct inquiries to 30 Frost Street, Cambridge, MA 02140
55 */
56
57
58 #include "mconf.h"
59
60 /* x exp(-x) shi(x), inverted interval 8 to 18 */
61 static double S1[] = {
62 1.83889230173399459482E-17,
63 -9.55485532279655569575E-17,
64 2.04326105980879882648E-16,
65 1.09896949074905343022E-15,
66 -1.31313534344092599234E-14,
67 5.93976226264314278932E-14,
68 -3.47197010497749154755E-14,
69 -1.40059764613117131000E-12,
70 9.49044626224223543299E-12,
71 -1.61596181145435454033E-11,
72 -1.77899784436430310321E-10,
73 1.35455469767246947469E-9,
74 -1.03257121792819495123E-9,
75 -3.56699611114982536845E-8,
76 1.44818877384267342057E-7,
77 7.82018215184051295296E-7,
78 -5.39919118403805073710E-6,
79 -3.12458202168959833422E-5,
80 8.90136741950727517826E-5,
81 2.02558474743846862168E-3,
82 2.96064440855633256972E-2,
83 1.11847751047257036625E0
84 };
85
86 /* x exp(-x) shi(x), inverted interval 18 to 88 */
87 static double S2[] = {
88 -1.05311574154850938805E-17,
89 2.62446095596355225821E-17,
90 8.82090135625368160657E-17,
91 -3.38459811878103047136E-16,
92 -8.30608026366935789136E-16,
93 3.93397875437050071776E-15,
94 1.01765565969729044505E-14,
95 -4.21128170307640802703E-14,
96 -1.60818204519802480035E-13,
97 3.34714954175994481761E-13,
98 2.72600352129153073807E-12,
99 1.66894954752839083608E-12,
100 -3.49278141024730899554E-11,
101 -1.58580661666482709598E-10,
102 -1.79289437183355633342E-10,
103 1.76281629144264523277E-9,
104 1.69050228879421288846E-8,
105 1.25391771228487041649E-7,
106 1.16229947068677338732E-6,
107 1.61038260117376323993E-5,
108 3.49810375601053973070E-4,
109 1.28478065259647610779E-2,
110 1.03665722588798326712E0
111 };
112
113 /* x exp(-x) chin(x), inverted interval 8 to 18 */
114 static double C1[] = {
115 -8.12435385225864036372E-18,
116 2.17586413290339214377E-17,
117 5.22624394924072204667E-17,
118 -9.48812110591690559363E-16,
119 5.35546311647465209166E-15,
120 -1.21009970113732918701E-14,
121 -6.00865178553447437951E-14,
122 7.16339649156028587775E-13,
123 -2.93496072607599856104E-12,
124 -1.40359438136491256904E-12,
125 8.76302288609054966081E-11,
126 -4.40092476213282340617E-10,
127 -1.87992075640569295479E-10,
128 1.31458150989474594064E-8,
129 -4.75513930924765465590E-8,
130 -2.21775018801848880741E-7,
131 1.94635531373272490962E-6,
132 4.33505889257316408893E-6,
133 -6.13387001076494349496E-5,
134 -3.13085477492997465138E-4,
135 4.97164789823116062801E-4,
136 2.64347496031374526641E-2,
137 1.11446150876699213025E0
138 };
139
140 /* x exp(-x) chin(x), inverted interval 18 to 88 */
141 static double C2[] = {
142 8.06913408255155572081E-18,
143 -2.08074168180148170312E-17,
144 -5.98111329658272336816E-17,
145 2.68533951085945765591E-16,
146 4.52313941698904694774E-16,
147 -3.10734917335299464535E-15,
148 -4.42823207332531972288E-15,
149 3.49639695410806959872E-14,
150 6.63406731718911586609E-14,
151 -3.71902448093119218395E-13,
152 -1.27135418132338309016E-12,
153 2.74851141935315395333E-12,
154 2.33781843985453438400E-11,
155 2.71436006377612442764E-11,
156 -2.56600180000355990529E-10,
157 -1.61021375163803438552E-9,
158 -4.72543064876271773512E-9,
159 -3.00095178028681682282E-9,
160 7.79387474390914922337E-8,
161 1.06942765566401507066E-6,
162 1.59503164802313196374E-5,
163 3.49592575153777996871E-4,
164 1.28475387530065247392E-2,
165 1.03665693917934275131E0
166 };
167
168 static double hyp3f0(double a1, double a2, double a3, double z);
169
170 /* Sine and cosine integrals */
171
172 extern double MACHEP;
173
shichi(x,si,ci)174 int shichi(x, si, ci)
175 double x;
176 double *si, *ci;
177 {
178 double k, z, c, s, a, b;
179 short sign;
180
181 if (x < 0.0) {
182 sign = -1;
183 x = -x;
184 }
185 else
186 sign = 0;
187
188
189 if (x == 0.0) {
190 *si = 0.0;
191 *ci = -NPY_INFINITY;
192 return (0);
193 }
194
195 if (x >= 8.0)
196 goto chb;
197
198 if (x >= 88.0)
199 goto asymp;
200
201 z = x * x;
202
203 /* Direct power series expansion */
204 a = 1.0;
205 s = 1.0;
206 c = 0.0;
207 k = 2.0;
208
209 do {
210 a *= z / k;
211 c += a / k;
212 k += 1.0;
213 a /= k;
214 s += a / k;
215 k += 1.0;
216 }
217 while (fabs(a / s) > MACHEP);
218
219 s *= x;
220 goto done;
221
222
223 chb:
224 /* Chebyshev series expansions */
225 if (x < 18.0) {
226 a = (576.0 / x - 52.0) / 10.0;
227 k = exp(x) / x;
228 s = k * chbevl(a, S1, 22);
229 c = k * chbevl(a, C1, 23);
230 goto done;
231 }
232
233 if (x <= 88.0) {
234 a = (6336.0 / x - 212.0) / 70.0;
235 k = exp(x) / x;
236 s = k * chbevl(a, S2, 23);
237 c = k * chbevl(a, C2, 24);
238 goto done;
239 }
240
241 asymp:
242 if (x > 1000) {
243 *si = NPY_INFINITY;
244 *ci = NPY_INFINITY;
245 }
246 else {
247 /* Asymptotic expansions
248 * http://functions.wolfram.com/GammaBetaErf/CoshIntegral/06/02/
249 * http://functions.wolfram.com/GammaBetaErf/SinhIntegral/06/02/0001/
250 */
251 a = hyp3f0(0.5, 1, 1, 4.0/(x*x));
252 b = hyp3f0(1, 1, 1.5, 4.0/(x*x));
253 *si = cosh(x)/x * a + sinh(x)/(x*x) * b;
254 *ci = sinh(x)/x * a + cosh(x)/(x*x) * b;
255 }
256 if (sign) {
257 *si = -*si;
258 }
259 return 0;
260
261 done:
262 if (sign)
263 s = -s;
264
265 *si = s;
266
267 *ci = NPY_EULER + log(x) + c;
268 return (0);
269 }
270
271
272 /*
273 * Evaluate 3F0(a1, a2, a3; z)
274 *
275 * The series is only asymptotic, so this requires z large enough.
276 */
hyp3f0(double a1,double a2,double a3,double z)277 static double hyp3f0(double a1, double a2, double a3, double z)
278 {
279 int n, maxiter;
280 double err, sum, term, m;
281
282 m = pow(z, -1.0/3);
283 if (m < 50) {
284 maxiter = m;
285 }
286 else {
287 maxiter = 50;
288 }
289
290 term = 1.0;
291 sum = term;
292 for (n = 0; n < maxiter; ++n) {
293 term *= (a1 + n) * (a2 + n) * (a3 + n) * z / (n + 1);
294 sum += term;
295 if (fabs(term) < 1e-13 * fabs(sum) || term == 0) {
296 break;
297 }
298 }
299
300 err = fabs(term);
301
302 if (err > 1e-13 * fabs(sum)) {
303 return NPY_NAN;
304 }
305
306 return sum;
307 }
308