1 /* specfunc/debye.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 /* augmented to n=5 and 6 2005-11-08 by R. J. Mathar, http://www.strw.leidenuniv.nl/~mathar */
22
23 #include <config.h>
24 #include <gsl/gsl_math.h>
25 #include <gsl/gsl_errno.h>
26 #include <gsl/gsl_sf_debye.h>
27
28 #include "error.h"
29 #include "check.h"
30
31 #include "chebyshev.h"
32 #include "cheb_eval.c"
33
34 static double adeb1_data[17] = {
35 2.4006597190381410194,
36 0.1937213042189360089,
37 -0.62329124554895770e-02,
38 0.3511174770206480e-03,
39 -0.228222466701231e-04,
40 0.15805467875030e-05,
41 -0.1135378197072e-06,
42 0.83583361188e-08,
43 -0.6264424787e-09,
44 0.476033489e-10,
45 -0.36574154e-11,
46 0.2835431e-12,
47 -0.221473e-13,
48 0.17409e-14,
49 -0.1376e-15,
50 0.109e-16,
51 -0.9e-18
52 };
53 static cheb_series adeb1_cs = {
54 adeb1_data,
55 16,
56 -1.0, 1.0,
57 9
58 };
59
60 static double adeb2_data[18] = {
61 2.5943810232570770282,
62 0.2863357204530719834,
63 -0.102062656158046713e-01,
64 0.6049109775346844e-03,
65 -0.405257658950210e-04,
66 0.28633826328811e-05,
67 -0.2086394303065e-06,
68 0.155237875826e-07,
69 -0.11731280087e-08,
70 0.897358589e-10,
71 -0.69317614e-11,
72 0.5398057e-12,
73 -0.423241e-13,
74 0.33378e-14,
75 -0.2645e-15,
76 0.211e-16,
77 -0.17e-17,
78 0.1e-18
79 };
80 static cheb_series adeb2_cs = {
81 adeb2_data,
82 17,
83 -1.0, 1.0,
84 10
85 };
86
87 static double adeb3_data[17] = {
88 2.707737068327440945,
89 0.340068135211091751,
90 -0.12945150184440869e-01,
91 0.7963755380173816e-03,
92 -0.546360009590824e-04,
93 0.39243019598805e-05,
94 -0.2894032823539e-06,
95 0.217317613962e-07,
96 -0.16542099950e-08,
97 0.1272796189e-09,
98 -0.987963460e-11,
99 0.7725074e-12,
100 -0.607797e-13,
101 0.48076e-14,
102 -0.3820e-15,
103 0.305e-16,
104 -0.24e-17
105 };
106 static cheb_series adeb3_cs = {
107 adeb3_data,
108 16,
109 -1.0, 1.0,
110 10
111 };
112
113 static double adeb4_data[17] = {
114 2.781869415020523460,
115 0.374976783526892863,
116 -0.14940907399031583e-01,
117 0.945679811437042e-03,
118 -0.66132916138933e-04,
119 0.4815632982144e-05,
120 -0.3588083958759e-06,
121 0.271601187416e-07,
122 -0.20807099122e-08,
123 0.1609383869e-09,
124 -0.125470979e-10,
125 0.9847265e-12,
126 -0.777237e-13,
127 0.61648e-14,
128 -0.4911e-15,
129 0.393e-16,
130 -0.32e-17
131 };
132 static cheb_series adeb4_cs = {
133 adeb4_data,
134 16,
135 -1.0, 1.0,
136 10
137 };
138
139 static double adeb5_data[17] = {
140 2.8340269546834530149,
141 0.3994098857106266445,
142 -0.164566764773099646e-1,
143 0.10652138340664541e-2,
144 -0.756730374875418e-4,
145 0.55745985240273e-5,
146 -0.4190692330918e-6,
147 0.319456143678e-7,
148 -0.24613318171e-8,
149 0.1912801633e-9,
150 -0.149720049e-10,
151 0.11790312e-11,
152 -0.933329e-13,
153 0.74218e-14,
154 -0.5925e-15,
155 0.475e-16,
156 -0.39e-17
157 };
158 static cheb_series adeb5_cs = {
159 adeb5_data,
160 16,
161 -1.0, 1.0,
162 10
163 };
164
165 static double adeb6_data[17] = {
166 2.8726727134130122113,
167 0.4174375352339027746,
168 -0.176453849354067873e-1,
169 0.11629852733494556e-2,
170 -0.837118027357117e-4,
171 0.62283611596189e-5,
172 -0.4718644465636e-6,
173 0.361950397806e-7,
174 -0.28030368010e-8,
175 0.2187681983e-9,
176 -0.171857387e-10,
177 0.13575809e-11,
178 -0.1077580e-12,
179 0.85893e-14,
180 -0.6872e-15,
181 0.552e-16,
182 -0.44e-17
183 };
184 static cheb_series adeb6_cs = {
185 adeb6_data,
186 16,
187 -1.0, 1.0,
188 10
189 };
190
191
192 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
193
gsl_sf_debye_1_e(const double x,gsl_sf_result * result)194 int gsl_sf_debye_1_e(const double x, gsl_sf_result * result)
195 {
196 const double val_infinity = 1.64493406684822644;
197 const double xcut = -GSL_LOG_DBL_MIN;
198
199 /* CHECK_POINTER(result) */
200
201 if(x < 0.0) {
202 DOMAIN_ERROR(result);
203 }
204 else if(x < 2.0*GSL_SQRT_DBL_EPSILON) {
205 result->val = 1.0 - 0.25*x + x*x/36.0;
206 result->err = GSL_DBL_EPSILON * fabs(result->val);
207 return GSL_SUCCESS;
208 }
209 else if(x <= 4.0) {
210 const double t = x*x/8.0 - 1.0;
211 gsl_sf_result c;
212 cheb_eval_e(&adeb1_cs, t, &c);
213 result->val = c.val - 0.25 * x;
214 result->err = c.err + 0.25 * x * GSL_DBL_EPSILON;
215 return GSL_SUCCESS;
216 }
217 else if(x < -(M_LN2 + GSL_LOG_DBL_EPSILON)) {
218 const int nexp = floor(xcut/x);
219 const double ex = exp(-x);
220 double sum = 0.0;
221 double xk = nexp * x;
222 double rk = nexp;
223 int i;
224 for(i=nexp; i>=1; i--) {
225 sum *= ex;
226 sum += (1.0 + 1.0/xk)/rk;
227 rk -= 1.0;
228 xk -= x;
229 }
230 result->val = val_infinity/x - sum*ex;
231 result->err = GSL_DBL_EPSILON * fabs(result->val);
232 return GSL_SUCCESS;
233 }
234 else if(x < xcut) {
235 result->val = (val_infinity - exp(-x)*(x+1.0)) / x;
236 result->err = GSL_DBL_EPSILON * fabs(result->val);
237 return GSL_SUCCESS;
238 }
239 else {
240 result->val = val_infinity/x;
241 result->err = GSL_DBL_EPSILON * fabs(result->val);
242 return GSL_SUCCESS;
243 }
244 }
245
246
gsl_sf_debye_2_e(const double x,gsl_sf_result * result)247 int gsl_sf_debye_2_e(const double x, gsl_sf_result * result)
248 {
249 const double val_infinity = 4.80822761263837714;
250 const double xcut = -GSL_LOG_DBL_MIN;
251
252 /* CHECK_POINTER(result) */
253
254 if(x < 0.0) {
255 DOMAIN_ERROR(result);
256 }
257 else if(x < 2.0*M_SQRT2*GSL_SQRT_DBL_EPSILON) {
258 result->val = 1.0 - x/3.0 + x*x/24.0;
259 result->err = GSL_DBL_EPSILON * result->val;
260 return GSL_SUCCESS;
261 }
262 else if(x <= 4.0) {
263 const double t = x*x/8.0 - 1.0;
264 gsl_sf_result c;
265 cheb_eval_e(&adeb2_cs, t, &c);
266 result->val = c.val - x/3.0;
267 result->err = c.err + GSL_DBL_EPSILON * x/3.0;
268 return GSL_SUCCESS;
269 }
270 else if(x < -(M_LN2 + GSL_LOG_DBL_EPSILON)) {
271 const int nexp = floor(xcut/x);
272 const double ex = exp(-x);
273 double xk = nexp * x;
274 double rk = nexp;
275 double sum = 0.0;
276 int i;
277 for(i=nexp; i>=1; i--) {
278 sum *= ex;
279 sum += (1.0 + 2.0/xk + 2.0/(xk*xk)) / rk;
280 rk -= 1.0;
281 xk -= x;
282 }
283 result->val = val_infinity/(x*x) - 2.0 * sum * ex;
284 result->err = GSL_DBL_EPSILON * fabs(result->val);
285 return GSL_SUCCESS;
286 }
287 else if(x < xcut) {
288 const double x2 = x*x;
289 const double sum = 2.0 + 2.0*x + x2;
290 result->val = (val_infinity - 2.0 * sum * exp(-x)) / x2;
291 result->err = GSL_DBL_EPSILON * fabs(result->val);
292 return GSL_SUCCESS;
293 }
294 else {
295 result->val = (val_infinity/x)/x;
296 result->err = GSL_DBL_EPSILON * result->val;
297 CHECK_UNDERFLOW(result);
298 return GSL_SUCCESS;
299 }
300 }
301
302
gsl_sf_debye_3_e(const double x,gsl_sf_result * result)303 int gsl_sf_debye_3_e(const double x, gsl_sf_result * result)
304 {
305 const double val_infinity = 19.4818182068004875;
306 const double xcut = -GSL_LOG_DBL_MIN;
307
308 /* CHECK_POINTER(result) */
309
310 if(x < 0.0) {
311 DOMAIN_ERROR(result);
312 }
313 else if(x < 2.0*M_SQRT2*GSL_SQRT_DBL_EPSILON) {
314 result->val = 1.0 - 3.0*x/8.0 + x*x/20.0;
315 result->err = GSL_DBL_EPSILON * result->val;
316 return GSL_SUCCESS;
317 }
318 else if(x <= 4.0) {
319 const double t = x*x/8.0 - 1.0;
320 gsl_sf_result c;
321 cheb_eval_e(&adeb3_cs, t, &c);
322 result->val = c.val - 0.375*x;
323 result->err = c.err + GSL_DBL_EPSILON * 0.375*x;
324 return GSL_SUCCESS;
325 }
326 else if(x < -(M_LN2 + GSL_LOG_DBL_EPSILON)) {
327 const int nexp = floor(xcut/x);
328 const double ex = exp(-x);
329 double xk = nexp * x;
330 double rk = nexp;
331 double sum = 0.0;
332 int i;
333 for(i=nexp; i>=1; i--) {
334 double xk_inv = 1.0/xk;
335 sum *= ex;
336 sum += (((6.0*xk_inv + 6.0)*xk_inv + 3.0)*xk_inv + 1.0) / rk;
337 rk -= 1.0;
338 xk -= x;
339 }
340 result->val = val_infinity/(x*x*x) - 3.0 * sum * ex;
341 result->err = GSL_DBL_EPSILON * result->val;
342 return GSL_SUCCESS;
343 }
344 else if(x < xcut) {
345 const double x3 = x*x*x;
346 const double sum = 6.0 + 6.0*x + 3.0*x*x + x3;
347 result->val = (val_infinity - 3.0 * sum * exp(-x)) / x3;
348 result->err = GSL_DBL_EPSILON * result->val;
349 return GSL_SUCCESS;
350 }
351 else {
352 result->val = ((val_infinity/x)/x)/x;
353 result->err = GSL_DBL_EPSILON * result->val;
354 CHECK_UNDERFLOW(result);
355 return GSL_SUCCESS;
356 }
357 }
358
359
gsl_sf_debye_4_e(const double x,gsl_sf_result * result)360 int gsl_sf_debye_4_e(const double x, gsl_sf_result * result)
361 {
362 const double val_infinity = 99.5450644937635129;
363 const double xcut = -GSL_LOG_DBL_MIN;
364
365 /* CHECK_POINTER(result) */
366
367 if(x < 0.0) {
368 DOMAIN_ERROR(result);
369 }
370 else if(x < 2.0*M_SQRT2*GSL_SQRT_DBL_EPSILON) {
371 result->val = 1.0 - 2.0*x/5.0 + x*x/18.0;
372 result->err = GSL_DBL_EPSILON * result->val;
373 return GSL_SUCCESS;
374 }
375 else if(x <= 4.0) {
376 const double t = x*x/8.0 - 1.0;
377 gsl_sf_result c;
378 cheb_eval_e(&adeb4_cs, t, &c);
379 result->val = c.val - 2.0*x/5.0;
380 result->err = c.err + GSL_DBL_EPSILON * 2.0*x/5.0;
381 return GSL_SUCCESS;
382 }
383 else if(x < -(M_LN2 + GSL_LOG_DBL_EPSILON)) {
384 const int nexp = floor(xcut/x);
385 const double ex = exp(-x);
386 double xk = nexp * x;
387 double rk = nexp;
388 double sum = 0.0;
389 int i;
390 for(i=nexp; i>=1; i--) {
391 double xk_inv = 1.0/xk;
392 sum *= ex;
393 sum += ((((24.0*xk_inv + 24.0)*xk_inv + 12.0)*xk_inv + 4.0)*xk_inv + 1.0) / rk;
394 rk -= 1.0;
395 xk -= x;
396 }
397 result->val = val_infinity/(x*x*x*x) - 4.0 * sum * ex;
398 result->err = GSL_DBL_EPSILON * result->val;
399 return GSL_SUCCESS;
400 }
401 else if(x < xcut) {
402 const double x2 = x*x;
403 const double x4 = x2*x2;
404 const double sum = 24.0 + 24.0*x + 12.0*x2 + 4.0*x2*x + x4;
405 result->val = (val_infinity - 4.0 * sum * exp(-x)) / x4;
406 result->err = GSL_DBL_EPSILON * result->val;
407 return GSL_SUCCESS;
408 }
409 else {
410 result->val = (((val_infinity/x)/x)/x)/x;
411 result->err = GSL_DBL_EPSILON * result->val;
412 CHECK_UNDERFLOW(result);
413 return GSL_SUCCESS;
414 }
415 }
416
gsl_sf_debye_5_e(const double x,gsl_sf_result * result)417 int gsl_sf_debye_5_e(const double x, gsl_sf_result * result)
418 {
419 const double val_infinity = 610.405837190669483828710757875 ;
420 const double xcut = -GSL_LOG_DBL_MIN;
421
422 /* CHECK_POINTER(result) */
423
424 if(x < 0.0) {
425 DOMAIN_ERROR(result);
426 }
427 else if(x < 2.0*M_SQRT2*GSL_SQRT_DBL_EPSILON) {
428 result->val = 1.0 - 5.0*x/12.0 + 5.0*x*x/84.0;
429 result->err = GSL_DBL_EPSILON * result->val;
430 return GSL_SUCCESS;
431 }
432 else if(x <= 4.0) {
433 const double t = x*x/8.0 - 1.0;
434 gsl_sf_result c;
435 cheb_eval_e(&adeb5_cs, t, &c);
436 result->val = c.val - 5.0*x/12.0;
437 result->err = c.err + GSL_DBL_EPSILON * 5.0*x/12.0;
438 return GSL_SUCCESS;
439 }
440 else if(x < -(M_LN2 + GSL_LOG_DBL_EPSILON)) {
441 const int nexp = floor(xcut/x);
442 const double ex = exp(-x);
443 double xk = nexp * x;
444 double rk = nexp;
445 double sum = 0.0;
446 int i;
447 for(i=nexp; i>=1; i--) {
448 double xk_inv = 1.0/xk;
449 sum *= ex;
450 sum += (((((120.0*xk_inv + 120.0)*xk_inv + 60.0)*xk_inv + 20.0)*xk_inv + 5.0)*xk_inv+ 1.0) / rk;
451 rk -= 1.0;
452 xk -= x;
453 }
454 result->val = val_infinity/(x*x*x*x*x) - 5.0 * sum * ex;
455 result->err = GSL_DBL_EPSILON * result->val;
456 return GSL_SUCCESS;
457 }
458 else if(x < xcut) {
459 const double x2 = x*x;
460 const double x4 = x2*x2;
461 const double x5 = x4*x;
462 const double sum = 120.0 + 120.0*x + 60.0*x2 + 20.0*x2*x + 5.0*x4 + x5;
463 result->val = (val_infinity - 5.0 * sum * exp(-x)) / x5;
464 result->err = GSL_DBL_EPSILON * result->val;
465 return GSL_SUCCESS;
466 }
467 else {
468 result->val = ((((val_infinity/x)/x)/x)/x)/x;
469 result->err = GSL_DBL_EPSILON * result->val;
470 CHECK_UNDERFLOW(result);
471 return GSL_SUCCESS;
472 }
473 }
474
gsl_sf_debye_6_e(const double x,gsl_sf_result * result)475 int gsl_sf_debye_6_e(const double x, gsl_sf_result * result)
476 {
477 const double val_infinity = 4356.06887828990661194792541535 ;
478 const double xcut = -GSL_LOG_DBL_MIN;
479
480 /* CHECK_POINTER(result) */
481
482 if(x < 0.0) {
483 DOMAIN_ERROR(result);
484 }
485 else if(x < 2.0*M_SQRT2*GSL_SQRT_DBL_EPSILON) {
486 result->val = 1.0 - 3.0*x/7.0 + x*x/16.0;
487 result->err = GSL_DBL_EPSILON * result->val;
488 return GSL_SUCCESS;
489 }
490 else if(x <= 4.0) {
491 const double t = x*x/8.0 - 1.0;
492 gsl_sf_result c;
493 cheb_eval_e(&adeb6_cs, t, &c);
494 result->val = c.val - 3.0*x/7.0;
495 result->err = c.err + GSL_DBL_EPSILON * 3.0*x/7.0;
496 return GSL_SUCCESS;
497 }
498 else if(x < -(M_LN2 + GSL_LOG_DBL_EPSILON)) {
499 const int nexp = floor(xcut/x);
500 const double ex = exp(-x);
501 double xk = nexp * x;
502 double rk = nexp;
503 double sum = 0.0;
504 int i;
505 for(i=nexp; i>=1; i--) {
506 double xk_inv = 1.0/xk;
507 sum *= ex;
508 sum += ((((((720.0*xk_inv + 720.0)*xk_inv + 360.0)*xk_inv + 120.0)*xk_inv + 30.0)*xk_inv+ 6.0)*xk_inv+ 1.0) / rk;
509 rk -= 1.0;
510 xk -= x;
511 }
512 result->val = val_infinity/(x*x*x*x*x*x) - 6.0 * sum * ex;
513 result->err = GSL_DBL_EPSILON * result->val;
514 return GSL_SUCCESS;
515 }
516 else if(x < xcut) {
517 const double x2 = x*x;
518 const double x4 = x2*x2;
519 const double x6 = x4*x2;
520 const double sum = 720.0 + 720.0*x + 360.0*x2 + 120.0*x2*x + 30.0*x4 + 6.0*x4*x +x6 ;
521 result->val = (val_infinity - 6.0 * sum * exp(-x)) / x6;
522 result->err = GSL_DBL_EPSILON * result->val;
523 return GSL_SUCCESS;
524 }
525 else {
526 result->val = (((((val_infinity/x)/x)/x)/x)/x)/x ;
527 result->err = GSL_DBL_EPSILON * result->val;
528 CHECK_UNDERFLOW(result);
529 return GSL_SUCCESS;
530 }
531 }
532
533
534 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
535
536 #include "eval.h"
537
gsl_sf_debye_1(const double x)538 double gsl_sf_debye_1(const double x)
539 {
540 EVAL_RESULT(gsl_sf_debye_1_e(x, &result));
541 }
542
gsl_sf_debye_2(const double x)543 double gsl_sf_debye_2(const double x)
544 {
545 EVAL_RESULT(gsl_sf_debye_2_e(x, &result));
546 }
547
gsl_sf_debye_3(const double x)548 double gsl_sf_debye_3(const double x)
549 {
550 EVAL_RESULT(gsl_sf_debye_3_e(x, &result));
551 }
552
gsl_sf_debye_4(const double x)553 double gsl_sf_debye_4(const double x)
554 {
555 EVAL_RESULT(gsl_sf_debye_4_e(x, &result));
556 }
557
gsl_sf_debye_5(const double x)558 double gsl_sf_debye_5(const double x)
559 {
560 EVAL_RESULT(gsl_sf_debye_5_e(x, &result));
561 }
562
gsl_sf_debye_6(const double x)563 double gsl_sf_debye_6(const double x)
564 {
565 EVAL_RESULT(gsl_sf_debye_6_e(x, &result));
566 }
567