1 /* specfunc/sinint.c
2 *
3 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 3 of the License, or (at
8 * your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 */
19
20 /* Author: G. Jungman */
21
22 #include <config.h>
23 #include <gsl/gsl_math.h>
24 #include <gsl/gsl_errno.h>
25 #include <gsl/gsl_sf_trig.h>
26 #include <gsl/gsl_sf_expint.h>
27
28 #include "error.h"
29
30 #include "chebyshev.h"
31 #include "cheb_eval.c"
32
33 /*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/
34
35 /* based on SLATEC r9sifg.f, W. Fullerton */
36
37 /*
38 series for f1 on the interval 2.00000e-02 to 6.25000e-02
39 with weighted error 2.82e-17
40 log weighted error 16.55
41 significant figures required 15.36
42 decimal places required 17.20
43 */
44 static double f1_data[20] = {
45 -0.1191081969051363610,
46 -0.0247823144996236248,
47 0.0011910281453357821,
48 -0.0000927027714388562,
49 0.0000093373141568271,
50 -0.0000011058287820557,
51 0.0000001464772071460,
52 -0.0000000210694496288,
53 0.0000000032293492367,
54 -0.0000000005206529618,
55 0.0000000000874878885,
56 -0.0000000000152176187,
57 0.0000000000027257192,
58 -0.0000000000005007053,
59 0.0000000000000940241,
60 -0.0000000000000180014,
61 0.0000000000000035063,
62 -0.0000000000000006935,
63 0.0000000000000001391,
64 -0.0000000000000000282
65 };
66 static cheb_series f1_cs = {
67 f1_data,
68 19,
69 -1, 1,
70 10
71 };
72
73 /*
74
75 series for f2 on the interval 0.00000e+00 to 2.00000e-02
76 with weighted error 4.32e-17
77 log weighted error 16.36
78 significant figures required 14.75
79 decimal places required 17.10
80 */
81 static double f2_data[29] = {
82 -0.0348409253897013234,
83 -0.0166842205677959686,
84 0.0006752901241237738,
85 -0.0000535066622544701,
86 0.0000062693421779007,
87 -0.0000009526638801991,
88 0.0000001745629224251,
89 -0.0000000368795403065,
90 0.0000000087202677705,
91 -0.0000000022601970392,
92 0.0000000006324624977,
93 -0.0000000001888911889,
94 0.0000000000596774674,
95 -0.0000000000198044313,
96 0.0000000000068641396,
97 -0.0000000000024731020,
98 0.0000000000009226360,
99 -0.0000000000003552364,
100 0.0000000000001407606,
101 -0.0000000000000572623,
102 0.0000000000000238654,
103 -0.0000000000000101714,
104 0.0000000000000044259,
105 -0.0000000000000019634,
106 0.0000000000000008868,
107 -0.0000000000000004074,
108 0.0000000000000001901,
109 -0.0000000000000000900,
110 0.0000000000000000432
111 };
112 static cheb_series f2_cs = {
113 f2_data,
114 28,
115 -1, 1,
116 14
117 };
118
119 /*
120
121 series for g1 on the interval 2.00000e-02 to 6.25000e-02
122 with weighted error 5.48e-17
123 log weighted error 16.26
124 significant figures required 15.47
125 decimal places required 16.92
126 */
127 static double g1_data[21] = {
128 -0.3040578798253495954,
129 -0.0566890984597120588,
130 0.0039046158173275644,
131 -0.0003746075959202261,
132 0.0000435431556559844,
133 -0.0000057417294453025,
134 0.0000008282552104503,
135 -0.0000001278245892595,
136 0.0000000207978352949,
137 -0.0000000035313205922,
138 0.0000000006210824236,
139 -0.0000000001125215474,
140 0.0000000000209088918,
141 -0.0000000000039715832,
142 0.0000000000007690431,
143 -0.0000000000001514697,
144 0.0000000000000302892,
145 -0.0000000000000061400,
146 0.0000000000000012601,
147 -0.0000000000000002615,
148 0.0000000000000000548
149 };
150 static cheb_series g1_cs = {
151 g1_data,
152 20,
153 -1, 1,
154 13
155 };
156
157 /*
158
159 series for g2 on the interval 0.00000e+00 to 2.00000e-02
160 with weighted error 5.01e-17
161 log weighted error 16.30
162 significant figures required 15.12
163 decimal places required 17.07
164 */
165 static double g2_data[34] = {
166 -0.0967329367532432218,
167 -0.0452077907957459871,
168 0.0028190005352706523,
169 -0.0002899167740759160,
170 0.0000407444664601121,
171 -0.0000071056382192354,
172 0.0000014534723163019,
173 -0.0000003364116512503,
174 0.0000000859774367886,
175 -0.0000000238437656302,
176 0.0000000070831906340,
177 -0.0000000022318068154,
178 0.0000000007401087359,
179 -0.0000000002567171162,
180 0.0000000000926707021,
181 -0.0000000000346693311,
182 0.0000000000133950573,
183 -0.0000000000053290754,
184 0.0000000000021775312,
185 -0.0000000000009118621,
186 0.0000000000003905864,
187 -0.0000000000001708459,
188 0.0000000000000762015,
189 -0.0000000000000346151,
190 0.0000000000000159996,
191 -0.0000000000000075213,
192 0.0000000000000035970,
193 -0.0000000000000017530,
194 0.0000000000000008738,
195 -0.0000000000000004487,
196 0.0000000000000002397,
197 -0.0000000000000001347,
198 0.0000000000000000801,
199 -0.0000000000000000501
200 };
201 static cheb_series g2_cs = {
202 g2_data,
203 33,
204 -1, 1,
205 20
206 };
207
208
209 /* x >= 4.0 */
fg_asymp(const double x,gsl_sf_result * f,gsl_sf_result * g)210 static void fg_asymp(const double x, gsl_sf_result * f, gsl_sf_result * g)
211 {
212 /*
213 xbig = sqrt (1.0/r1mach(3))
214 xmaxf = exp (amin1(-alog(r1mach(1)), alog(r1mach(2))) - 0.01)
215 xmaxg = 1.0/sqrt(r1mach(1))
216 xbnd = sqrt(50.0)
217 */
218 const double xbig = 1.0/GSL_SQRT_DBL_EPSILON;
219 const double xmaxf = 1.0/GSL_DBL_MIN;
220 const double xmaxg = 1.0/GSL_SQRT_DBL_MIN;
221 const double xbnd = 7.07106781187;
222
223 const double x2 = x*x;
224
225 if(x <= xbnd) {
226 gsl_sf_result result_c1;
227 gsl_sf_result result_c2;
228 cheb_eval_e(&f1_cs, (1.0/x2-0.04125)/0.02125, &result_c1);
229 cheb_eval_e(&g1_cs, (1.0/x2-0.04125)/0.02125, &result_c2);
230 f->val = (1.0 + result_c1.val)/x;
231 g->val = (1.0 + result_c2.val)/x2;
232 f->err = result_c1.err/x + 2.0 * GSL_DBL_EPSILON * fabs(f->val);
233 g->err = result_c2.err/x2 + 2.0 * GSL_DBL_EPSILON * fabs(g->val);
234 }
235 else if(x <= xbig) {
236 gsl_sf_result result_c1;
237 gsl_sf_result result_c2;
238 cheb_eval_e(&f2_cs, 100.0/x2-1.0, &result_c1);
239 cheb_eval_e(&g2_cs, 100.0/x2-1.0, &result_c2);
240 f->val = (1.0 + result_c1.val)/x;
241 g->val = (1.0 + result_c2.val)/x2;
242 f->err = result_c1.err/x + 2.0 * GSL_DBL_EPSILON * fabs(f->val);
243 g->err = result_c2.err/x2 + 2.0 * GSL_DBL_EPSILON * fabs(g->val);
244 }
245 else {
246 f->val = (x < xmaxf ? 1.0/x : 0.0);
247 g->val = (x < xmaxg ? 1.0/x2 : 0.0);
248 f->err = 2.0 * GSL_DBL_EPSILON * fabs(f->val);
249 g->err = 2.0 * GSL_DBL_EPSILON * fabs(g->val);
250 }
251
252 return;
253 }
254
255
256 /* based on SLATEC si.f, W. Fullerton
257
258 series for si on the interval 0.00000e+00 to 1.60000e+01
259 with weighted error 1.22e-17
260 log weighted error 16.91
261 significant figures required 16.37
262 decimal places required 17.45
263 */
264
265 static double si_data[12] = {
266 -0.1315646598184841929,
267 -0.2776578526973601892,
268 0.0354414054866659180,
269 -0.0025631631447933978,
270 0.0001162365390497009,
271 -0.0000035904327241606,
272 0.0000000802342123706,
273 -0.0000000013562997693,
274 0.0000000000179440722,
275 -0.0000000000001908387,
276 0.0000000000000016670,
277 -0.0000000000000000122
278 };
279
280 static cheb_series si_cs = {
281 si_data,
282 11,
283 -1, 1,
284 9
285 };
286
287 /*
288 series for ci on the interval 0.00000e+00 to 1.60000e+01
289 with weighted error 1.94e-18
290 log weighted error 17.71
291 significant figures required 17.74
292 decimal places required 18.27
293 */
294 static double ci_data[13] = {
295 -0.34004281856055363156,
296 -1.03302166401177456807,
297 0.19388222659917082877,
298 -0.01918260436019865894,
299 0.00110789252584784967,
300 -0.00004157234558247209,
301 0.00000109278524300229,
302 -0.00000002123285954183,
303 0.00000000031733482164,
304 -0.00000000000376141548,
305 0.00000000000003622653,
306 -0.00000000000000028912,
307 0.00000000000000000194
308 };
309 static cheb_series ci_cs = {
310 ci_data,
311 12,
312 -1, 1,
313 9
314 };
315
316
317 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
318
gsl_sf_Si_e(const double x,gsl_sf_result * result)319 int gsl_sf_Si_e(const double x, gsl_sf_result * result)
320 {
321 double ax = fabs(x);
322
323 /* CHECK_POINTER(result) */
324
325 if(ax < GSL_SQRT_DBL_EPSILON) {
326 result->val = x;
327 result->err = 0.0;
328 return GSL_SUCCESS;
329 }
330 else if(ax <= 4.0) {
331 gsl_sf_result result_c;
332 cheb_eval_e(&si_cs, (x*x-8.0)*0.125, &result_c);
333 result->val = x * (0.75 + result_c.val);
334 result->err = ax * result_c.err;
335 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
336 return GSL_SUCCESS;
337 }
338 else {
339 /* Note there is no loss of precision
340 * here bcause of the leading constant.
341 */
342 gsl_sf_result f;
343 gsl_sf_result g;
344 fg_asymp(ax, &f, &g);
345 result->val = 0.5 * M_PI - f.val*cos(ax) - g.val*sin(ax);
346 result->err = f.err + g.err;
347 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
348 if(x < 0.0) result->val = -result->val;
349 return GSL_SUCCESS;
350 }
351 }
352
353
gsl_sf_Ci_e(const double x,gsl_sf_result * result)354 int gsl_sf_Ci_e(const double x, gsl_sf_result * result)
355 {
356 /* CHECK_POINTER(result) */
357
358 if(x <= 0.0) {
359 DOMAIN_ERROR(result);
360 }
361 else if(x <= 4.0) {
362 const double lx = log(x);
363 const double y = (x*x-8.0)*0.125;
364 gsl_sf_result result_c;
365 cheb_eval_e(&ci_cs, y, &result_c);
366 result->val = lx - 0.5 + result_c.val;
367 result->err = 2.0 * GSL_DBL_EPSILON * (fabs(lx) + 0.5) + result_c.err;
368 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
369 return GSL_SUCCESS;
370 }
371 else {
372 gsl_sf_result sin_result;
373 gsl_sf_result cos_result;
374 int stat_sin = gsl_sf_sin_e(x, &sin_result);
375 int stat_cos = gsl_sf_cos_e(x, &cos_result);
376 gsl_sf_result f;
377 gsl_sf_result g;
378 fg_asymp(x, &f, &g);
379 result->val = f.val*sin_result.val - g.val*cos_result.val;
380 result->err = fabs(f.err*sin_result.val);
381 result->err += fabs(g.err*cos_result.val);
382 result->err += fabs(f.val*sin_result.err);
383 result->err += fabs(g.val*cos_result.err);
384 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
385 return GSL_ERROR_SELECT_2(stat_sin, stat_cos);
386 }
387 }
388
389
390 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
391
392 #include "eval.h"
393
gsl_sf_Si(const double x)394 double gsl_sf_Si(const double x)
395 {
396 EVAL_RESULT(gsl_sf_Si_e(x, &result));
397 }
398
gsl_sf_Ci(const double x)399 double gsl_sf_Ci(const double x)
400 {
401 EVAL_RESULT(gsl_sf_Ci_e(x, &result));
402 }
403