diff options
Diffstat (limited to 'libm/double/dawsn.c')
| -rw-r--r-- | libm/double/dawsn.c | 392 |
1 files changed, 392 insertions, 0 deletions
diff --git a/libm/double/dawsn.c b/libm/double/dawsn.c new file mode 100644 index 000000000..4f8d27a0c --- /dev/null +++ b/libm/double/dawsn.c @@ -0,0 +1,392 @@ +/* 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 + * DEC 0,10 6000 7.4e-17 1.4e-17 + * + * + */ + +/* dawsn.c */ + + +/* +Cephes Math Library Release 2.8: June, 2000 +Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier +*/ + +#include <math.h> +/* Dawson's integral, interval 0 to 3.25 */ +#ifdef UNK +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, +}; +#endif +#ifdef DEC +static unsigned short AN[40] = { +0027107,0176630,0075752,0107612, +0030551,0070604,0166707,0127727, +0031647,0002210,0117120,0056376, +0033177,0156026,0141275,0140627, +0033516,0112200,0037035,0165515, +0035270,0150613,0016423,0105634, +0135536,0156227,0023515,0044413, +0037055,0015273,0105147,0064025, +0137273,0163145,0014460,0166465, +0040200,0000000,0000000,0000000, +}; +static unsigned short AD[44] = { +0027323,0067372,0115566,0131320, +0030714,0114432,0074206,0006637, +0032137,0160671,0044203,0026344, +0033252,0146656,0020247,0100231, +0034303,0003346,0123260,0022433, +0035252,0125460,0173041,0155415, +0036144,0113747,0125203,0124617, +0036744,0166232,0143671,0133670, +0037442,0127755,0162625,0000100, +0040023,0026736,0003604,0106265, +0040200,0000000,0000000,0000000, +}; +#endif +#ifdef IBMPC +static unsigned short AN[40] = { +0x51f1,0x0f7d,0xffb3,0x3da8, +0xf5fb,0x9db8,0x2e30,0x3e0d, +0x0ba0,0x13ca,0xe091,0x3e54, +0xb833,0xd857,0xfb82,0x3eaf, +0xbd6a,0x07c3,0xd290,0x3ec9, +0x7174,0x63a2,0x1a31,0x3f37, +0xa921,0xe4e9,0xdb92,0xbf4b, +0xed03,0x714c,0xa357,0x3fa5, +0x1da7,0xa326,0x7ccc,0xbfb7, +0x0000,0x0000,0x0000,0x3ff0, +}; +static unsigned short AD[44] = { +0xd65a,0x536e,0x6ddf,0x3dba, +0xc1b4,0x4f10,0x9323,0x3e19, +0x659c,0x2910,0xfc37,0x3e6b, +0xf013,0xc414,0x59b5,0x3eb5, +0x04a3,0xd4d6,0x60dc,0x3ef8, +0x3b62,0x1ec4,0x5566,0x3f35, +0x7532,0xf550,0x92fc,0x3f6c, +0x36f7,0x58f7,0x9d93,0x3f9c, +0xa008,0xbcb2,0x55fd,0x3fc4, +0x9197,0xc0f0,0x65bb,0x3fe2, +0x0000,0x0000,0x0000,0x3ff0, +}; +#endif +#ifdef MIEEE +static unsigned short AN[40] = { +0x3da8,0xffb3,0x0f7d,0x51f1, +0x3e0d,0x2e30,0x9db8,0xf5fb, +0x3e54,0xe091,0x13ca,0x0ba0, +0x3eaf,0xfb82,0xd857,0xb833, +0x3ec9,0xd290,0x07c3,0xbd6a, +0x3f37,0x1a31,0x63a2,0x7174, +0xbf4b,0xdb92,0xe4e9,0xa921, +0x3fa5,0xa357,0x714c,0xed03, +0xbfb7,0x7ccc,0xa326,0x1da7, +0x3ff0,0x0000,0x0000,0x0000, +}; +static unsigned short AD[44] = { +0x3dba,0x6ddf,0x536e,0xd65a, +0x3e19,0x9323,0x4f10,0xc1b4, +0x3e6b,0xfc37,0x2910,0x659c, +0x3eb5,0x59b5,0xc414,0xf013, +0x3ef8,0x60dc,0xd4d6,0x04a3, +0x3f35,0x5566,0x1ec4,0x3b62, +0x3f6c,0x92fc,0xf550,0x7532, +0x3f9c,0x9d93,0x58f7,0x36f7, +0x3fc4,0x55fd,0xbcb2,0xa008, +0x3fe2,0x65bb,0xc0f0,0x9197, +0x3ff0,0x0000,0x0000,0x0000, +}; +#endif + +/* interval 3.25 to 6.25 */ +#ifdef UNK +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, +}; +#endif +#ifdef DEC +static unsigned short BN[44] = { +0040002,0045342,0113762,0004360, +0137572,0120346,0172745,0144046, +0037300,0151134,0123440,0117047, +0136663,0025423,0014755,0046026, +0036157,0177561,0027535,0046744, +0135335,0161052,0071243,0146535, +0034426,0154060,0164506,0135625, +0133420,0005356,0100017,0151334, +0032303,0066137,0024013,0046212, +0131045,0016612,0066270,0047574, +0027435,0177025,0060625,0116363, +}; +static unsigned short BD[40] = { +/*0040200,0000000,0000000,0000000,*/ +0140041,0140101,0174552,0037073, +0037562,0061503,0124271,0160756, +0137131,0151760,0073210,0110534, +0036412,0170562,0117017,0155377, +0135570,0101374,0074056,0037276, +0034643,0145376,0001516,0060636, +0133630,0173540,0121344,0155231, +0032513,0005602,0134516,0007144, +0131250,0150540,0075747,0105341, +0027635,0177020,0012465,0125402, +}; +#endif +#ifdef IBMPC +static unsigned short BN[44] = { +0x411e,0x52fe,0x495c,0x3fe0, +0xb905,0xdebc,0x541c,0xbfcf, +0x13c5,0x94e4,0x1a4b,0x3fb8, +0xa983,0x633d,0x6562,0xbf96, +0xa9bd,0x25eb,0xffee,0x3f6d, +0x79ac,0x4e54,0xbc45,0xbf3b, +0xd773,0x1d28,0xdb06,0x3f02, +0xfa5b,0xd001,0x015d,0xbec2, +0x6991,0xe501,0x6d8b,0x3e78, +0x09f0,0x4d97,0xa3b1,0xbe24, +0xb39e,0xac32,0xbfc2,0x3dc3, +}; +static unsigned short BD[40] = { +/*0x0000,0x0000,0x0000,0x3ff0,*/ +0x47c7,0x3f2d,0x3808,0xbfe4, +0x3c3e,0x7517,0x4c68,0x3fce, +0x122b,0x0ed1,0x3a7e,0xbfab, +0xfb60,0x53c1,0x5e2e,0x3f81, +0xc7d8,0x8f05,0x105f,0xbf4f, +0xcc34,0xc069,0x795f,0x3f14, +0x9b53,0x145c,0x1eec,0xbed3, +0xc1cd,0x5729,0x6170,0x3e89, +0xf15c,0x0f7c,0x1a2c,0xbe35, +0xb560,0x02a6,0xbfc2,0x3dd3, +}; +#endif +#ifdef MIEEE +static unsigned short BN[44] = { +0x3fe0,0x495c,0x52fe,0x411e, +0xbfcf,0x541c,0xdebc,0xb905, +0x3fb8,0x1a4b,0x94e4,0x13c5, +0xbf96,0x6562,0x633d,0xa983, +0x3f6d,0xffee,0x25eb,0xa9bd, +0xbf3b,0xbc45,0x4e54,0x79ac, +0x3f02,0xdb06,0x1d28,0xd773, +0xbec2,0x015d,0xd001,0xfa5b, +0x3e78,0x6d8b,0xe501,0x6991, +0xbe24,0xa3b1,0x4d97,0x09f0, +0x3dc3,0xbfc2,0xac32,0xb39e, +}; +static unsigned short BD[40] = { +/*0x3ff0,0x0000,0x0000,0x0000,*/ +0xbfe4,0x3808,0x3f2d,0x47c7, +0x3fce,0x4c68,0x7517,0x3c3e, +0xbfab,0x3a7e,0x0ed1,0x122b, +0x3f81,0x5e2e,0x53c1,0xfb60, +0xbf4f,0x105f,0x8f05,0xc7d8, +0x3f14,0x795f,0xc069,0xcc34, +0xbed3,0x1eec,0x145c,0x9b53, +0x3e89,0x6170,0x5729,0xc1cd, +0xbe35,0x1a2c,0x0f7c,0xf15c, +0x3dd3,0xbfc2,0x02a6,0xb560, +}; +#endif + +/* 6.25 to infinity */ +#ifdef UNK +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, +}; +#endif +#ifdef DEC +static unsigned short CN[20] = { +0140027,0030427,0176477,0074402, +0040041,0012617,0112375,0162657, +0137461,0000761,0074120,0135160, +0036607,0004325,0117246,0115525, +0135377,0036345,0064750,0047732, +}; +static unsigned short CD[20] = { +/*0040200,0000000,0000000,0000000,*/ +0140454,0127521,0071653,0133415, +0040335,0144540,0016105,0045241, +0137711,0112053,0155034,0062237, +0037015,0002102,0177442,0074546, +0135577,0036345,0064750,0052152, +}; +#endif +#ifdef IBMPC +static unsigned short CN[20] = { +0xef20,0xffa7,0xe622,0xbfe2, +0xbcb6,0xf29f,0x22b1,0x3fe4, +0x174e,0x2f0a,0x203e,0xbfc6, +0xd36b,0xb3d4,0xe11a,0x3f90, +0x09fb,0xad3d,0xe79c,0xbf3f, +}; +static unsigned short CD[20] = { +/*0x0000,0x0000,0x0000,0x3ff0,*/ +0x76e2,0x2e75,0x95ea,0xc005, +0xa954,0x0388,0xb92c,0x3ffb, +0x8c94,0x7b43,0x3285,0xbfd9, +0x4f2d,0x5fe4,0xa088,0x3fa1, +0x0a8d,0xad3d,0xe79c,0xbf4f, +}; +#endif +#ifdef MIEEE +static unsigned short CN[20] = { +0xbfe2,0xe622,0xffa7,0xef20, +0x3fe4,0x22b1,0xf29f,0xbcb6, +0xbfc6,0x203e,0x2f0a,0x174e, +0x3f90,0xe11a,0xb3d4,0xd36b, +0xbf3f,0xe79c,0xad3d,0x09fb, +}; +static unsigned short CD[20] = { +/*0x3ff0,0x0000,0x0000,0x0000,*/ +0xc005,0x95ea,0x2e75,0x76e2, +0x3ffb,0xb92c,0x0388,0xa954, +0xbfd9,0x3285,0x7b43,0x8c94, +0x3fa1,0xa088,0x5fe4,0x4f2d, +0xbf4f,0xe79c,0xad3d,0x0a8d, +}; +#endif + +#ifdef ANSIPROT +extern double chbevl ( double, void *, int ); +extern double sqrt ( double ); +extern double fabs ( double ); +extern double polevl ( double, void *, int ); +extern double p1evl ( double, void *, int ); +#else +double chbevl(), sqrt(), fabs(), polevl(), p1evl(); +#endif +extern double PI, MACHEP; + +double dawsn( xx ) +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 ); +} |