161 lines
3.5 KiB
C
161 lines
3.5 KiB
C
/* dawsn.c
|
||
*
|
||
* Dawson's Integral
|
||
*
|
||
*
|
||
*
|
||
* SYNOPSIS:
|
||
*
|
||
* double x, y, dawsn();
|
||
*
|
||
* y = dawsn( x );
|
||
*
|
||
*
|
||
*
|
||
* DESCRIPTION:
|
||
*
|
||
* Approximates the integral
|
||
*
|
||
* x
|
||
* -
|
||
* 2 | | 2
|
||
* dawsn(x) = exp( -x ) | exp( t ) dt
|
||
* | |
|
||
* -
|
||
* 0
|
||
*
|
||
* Three different rational approximations are employed, for
|
||
* the intervals 0 to 3.25; 3.25 to 6.25; and 6.25 up.
|
||
*
|
||
*
|
||
* ACCURACY:
|
||
*
|
||
* Relative error:
|
||
* arithmetic domain # trials peak rms
|
||
* IEEE 0,10 10000 6.9e-16 1.0e-16
|
||
*
|
||
*
|
||
*/
|
||
|
||
/* dawsn.c */
|
||
|
||
|
||
/*
|
||
* Cephes Math Library Release 2.1: January, 1989
|
||
* Copyright 1984, 1987, 1989 by Stephen L. Moshier
|
||
* Direct inquiries to 30 Frost Street, Cambridge, MA 02140
|
||
*/
|
||
|
||
#include "mconf.h"
|
||
/* Dawson's integral, interval 0 to 3.25 */
|
||
static double AN[10] = {
|
||
1.13681498971755972054E-11,
|
||
8.49262267667473811108E-10,
|
||
1.94434204175553054283E-8,
|
||
9.53151741254484363489E-7,
|
||
3.07828309874913200438E-6,
|
||
3.52513368520288738649E-4,
|
||
-8.50149846724410912031E-4,
|
||
4.22618223005546594270E-2,
|
||
-9.17480371773452345351E-2,
|
||
9.99999999999999994612E-1,
|
||
};
|
||
|
||
static double AD[11] = {
|
||
2.40372073066762605484E-11,
|
||
1.48864681368493396752E-9,
|
||
5.21265281010541664570E-8,
|
||
1.27258478273186970203E-6,
|
||
2.32490249820789513991E-5,
|
||
3.25524741826057911661E-4,
|
||
3.48805814657162590916E-3,
|
||
2.79448531198828973716E-2,
|
||
1.58874241960120565368E-1,
|
||
5.74918629489320327824E-1,
|
||
1.00000000000000000539E0,
|
||
};
|
||
|
||
/* interval 3.25 to 6.25 */
|
||
static double BN[11] = {
|
||
5.08955156417900903354E-1,
|
||
-2.44754418142697847934E-1,
|
||
9.41512335303534411857E-2,
|
||
-2.18711255142039025206E-2,
|
||
3.66207612329569181322E-3,
|
||
-4.23209114460388756528E-4,
|
||
3.59641304793896631888E-5,
|
||
-2.14640351719968974225E-6,
|
||
9.10010780076391431042E-8,
|
||
-2.40274520828250956942E-9,
|
||
3.59233385440928410398E-11,
|
||
};
|
||
|
||
static double BD[10] = {
|
||
/* 1.00000000000000000000E0, */
|
||
-6.31839869873368190192E-1,
|
||
2.36706788228248691528E-1,
|
||
-5.31806367003223277662E-2,
|
||
8.48041718586295374409E-3,
|
||
-9.47996768486665330168E-4,
|
||
7.81025592944552338085E-5,
|
||
-4.55875153252442634831E-6,
|
||
1.89100358111421846170E-7,
|
||
-4.91324691331920606875E-9,
|
||
7.18466403235734541950E-11,
|
||
};
|
||
|
||
/* 6.25 to infinity */
|
||
static double CN[5] = {
|
||
-5.90592860534773254987E-1,
|
||
6.29235242724368800674E-1,
|
||
-1.72858975380388136411E-1,
|
||
1.64837047825189632310E-2,
|
||
-4.86827613020462700845E-4,
|
||
};
|
||
|
||
static double CD[5] = {
|
||
/* 1.00000000000000000000E0, */
|
||
-2.69820057197544900361E0,
|
||
1.73270799045947845857E0,
|
||
-3.93708582281939493482E-1,
|
||
3.44278924041233391079E-2,
|
||
-9.73655226040941223894E-4,
|
||
};
|
||
|
||
extern double MACHEP;
|
||
|
||
double dawsn(double xx)
|
||
{
|
||
double x, y;
|
||
int sign;
|
||
|
||
|
||
sign = 1;
|
||
if (xx < 0.0) {
|
||
sign = -1;
|
||
xx = -xx;
|
||
}
|
||
|
||
if (xx < 3.25) {
|
||
x = xx * xx;
|
||
y = xx * polevl(x, AN, 9) / polevl(x, AD, 10);
|
||
return (sign * y);
|
||
}
|
||
|
||
|
||
x = 1.0 / (xx * xx);
|
||
|
||
if (xx < 6.25) {
|
||
y = 1.0 / xx + x * polevl(x, BN, 10) / (p1evl(x, BD, 10) * xx);
|
||
return (sign * 0.5 * y);
|
||
}
|
||
|
||
|
||
if (xx > 1.0e9)
|
||
return ((sign * 0.5) / xx);
|
||
|
||
/* 6.25 to infinity */
|
||
y = 1.0 / xx + x * polevl(x, CN, 4) / (p1evl(x, CD, 5) * xx);
|
||
return (sign * 0.5 * y);
|
||
}
|