lib/libm/dawsn.c

   1 /*                                                      dawsn.c
   2  *
   3  *      Dawson's Integral
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * double x, y, dawsn();
  10  *
  11  * y = dawsn( x );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  * Approximates the integral
  18  *
  19  *                             x
  20  *                             -
  21  *                      2     | |        2
  22  *  dawsn(x)  =  exp( -x  )   |    exp( t  ) dt
  23  *                          | |
  24  *                           -
  25  *                           0
  26  *
  27  * Three different rational approximations are employed, for
  28  * the intervals 0 to 3.25; 3.25 to 6.25; and 6.25 up.
  29  *
  30  *
  31  * ACCURACY:
  32  *
  33  *                      Relative error:
  34  * arithmetic   domain     # trials      peak         rms
  35  *    IEEE      0,10        10000       6.9e-16     1.0e-16
  36  *    DEC       0,10         6000       7.4e-17     1.4e-17
  37  *
  38  *
  39  */
  40 \f
  41 /*                                                      dawsn.c */
  42
  43
  44 /*
  45 Cephes Math Library Release 2.8:  June, 2000
  46 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
  47 */
  48
  49 #include "mconf.h"
  50 /* Dawson's integral, interval 0 to 3.25 */
  51 #ifdef UNK
  52 static double AN[10] = {
  53  1.13681498971755972054E-11,
  54  8.49262267667473811108E-10,
  55  1.94434204175553054283E-8,
  56  9.53151741254484363489E-7,
  57  3.07828309874913200438E-6,
  58  3.52513368520288738649E-4,
  59 -8.50149846724410912031E-4,
  60  4.22618223005546594270E-2,
  61 -9.17480371773452345351E-2,
  62  9.99999999999999994612E-1,
  63 };
  64 static double AD[11] = {
  65  2.40372073066762605484E-11,
  66  1.48864681368493396752E-9,
  67  5.21265281010541664570E-8,
  68  1.27258478273186970203E-6,
  69  2.32490249820789513991E-5,
  70  3.25524741826057911661E-4,
  71  3.48805814657162590916E-3,
  72  2.79448531198828973716E-2,
  73  1.58874241960120565368E-1,
  74  5.74918629489320327824E-1,
  75  1.00000000000000000539E0,
  76 };
  77 #endif
  78 #ifdef DEC
  79 static unsigned short AN[40] = {
  80 0027107,0176630,0075752,0107612,
  81 0030551,0070604,0166707,0127727,
  82 0031647,0002210,0117120,0056376,
  83 0033177,0156026,0141275,0140627,
  84 0033516,0112200,0037035,0165515,
  85 0035270,0150613,0016423,0105634,
  86 0135536,0156227,0023515,0044413,
  87 0037055,0015273,0105147,0064025,
  88 0137273,0163145,0014460,0166465,
  89 0040200,0000000,0000000,0000000,
  90 };
  91 static unsigned short AD[44] = {
  92 0027323,0067372,0115566,0131320,
  93 0030714,0114432,0074206,0006637,
  94 0032137,0160671,0044203,0026344,
  95 0033252,0146656,0020247,0100231,
  96 0034303,0003346,0123260,0022433,
  97 0035252,0125460,0173041,0155415,
  98 0036144,0113747,0125203,0124617,
  99 0036744,0166232,0143671,0133670,
 100 0037442,0127755,0162625,0000100,
 101 0040023,0026736,0003604,0106265,
 102 0040200,0000000,0000000,0000000,
 103 };
 104 #endif
 105 #ifdef IBMPC
 106 static unsigned short AN[40] = {
 107 0x51f1,0x0f7d,0xffb3,0x3da8,
 108 0xf5fb,0x9db8,0x2e30,0x3e0d,
 109 0x0ba0,0x13ca,0xe091,0x3e54,
 110 0xb833,0xd857,0xfb82,0x3eaf,
 111 0xbd6a,0x07c3,0xd290,0x3ec9,
 112 0x7174,0x63a2,0x1a31,0x3f37,
 113 0xa921,0xe4e9,0xdb92,0xbf4b,
 114 0xed03,0x714c,0xa357,0x3fa5,
 115 0x1da7,0xa326,0x7ccc,0xbfb7,
 116 0x0000,0x0000,0x0000,0x3ff0,
 117 };
 118 static unsigned short AD[44] = {
 119 0xd65a,0x536e,0x6ddf,0x3dba,
 120 0xc1b4,0x4f10,0x9323,0x3e19,
 121 0x659c,0x2910,0xfc37,0x3e6b,
 122 0xf013,0xc414,0x59b5,0x3eb5,
 123 0x04a3,0xd4d6,0x60dc,0x3ef8,
 124 0x3b62,0x1ec4,0x5566,0x3f35,
 125 0x7532,0xf550,0x92fc,0x3f6c,
 126 0x36f7,0x58f7,0x9d93,0x3f9c,
 127 0xa008,0xbcb2,0x55fd,0x3fc4,
 128 0x9197,0xc0f0,0x65bb,0x3fe2,
 129 0x0000,0x0000,0x0000,0x3ff0,
 130 };
 131 #endif
 132 #ifdef MIEEE
 133 static unsigned short AN[40] = {
 134 0x3da8,0xffb3,0x0f7d,0x51f1,
 135 0x3e0d,0x2e30,0x9db8,0xf5fb,
 136 0x3e54,0xe091,0x13ca,0x0ba0,
 137 0x3eaf,0xfb82,0xd857,0xb833,
 138 0x3ec9,0xd290,0x07c3,0xbd6a,
 139 0x3f37,0x1a31,0x63a2,0x7174,
 140 0xbf4b,0xdb92,0xe4e9,0xa921,
 141 0x3fa5,0xa357,0x714c,0xed03,
 142 0xbfb7,0x7ccc,0xa326,0x1da7,
 143 0x3ff0,0x0000,0x0000,0x0000,
 144 };
 145 static unsigned short AD[44] = {
 146 0x3dba,0x6ddf,0x536e,0xd65a,
 147 0x3e19,0x9323,0x4f10,0xc1b4,
 148 0x3e6b,0xfc37,0x2910,0x659c,
 149 0x3eb5,0x59b5,0xc414,0xf013,
 150 0x3ef8,0x60dc,0xd4d6,0x04a3,
 151 0x3f35,0x5566,0x1ec4,0x3b62,
 152 0x3f6c,0x92fc,0xf550,0x7532,
 153 0x3f9c,0x9d93,0x58f7,0x36f7,
 154 0x3fc4,0x55fd,0xbcb2,0xa008,
 155 0x3fe2,0x65bb,0xc0f0,0x9197,
 156 0x3ff0,0x0000,0x0000,0x0000,
 157 };
 158 #endif
 159
 160 /* interval 3.25 to 6.25 */
 161 #ifdef UNK
 162 static double BN[11] = {
 163  5.08955156417900903354E-1,
 164 -2.44754418142697847934E-1,
 165  9.41512335303534411857E-2,
 166 -2.18711255142039025206E-2,
 167  3.66207612329569181322E-3,
 168 -4.23209114460388756528E-4,
 169  3.59641304793896631888E-5,
 170 -2.14640351719968974225E-6,
 171  9.10010780076391431042E-8,
 172 -2.40274520828250956942E-9,
 173  3.59233385440928410398E-11,
 174 };
 175 static double BD[10] = {
 176 /*  1.00000000000000000000E0,*/
 177 -6.31839869873368190192E-1,
 178  2.36706788228248691528E-1,
 179 -5.31806367003223277662E-2,
 180  8.48041718586295374409E-3,
 181 -9.47996768486665330168E-4,
 182  7.81025592944552338085E-5,
 183 -4.55875153252442634831E-6,
 184  1.89100358111421846170E-7,
 185 -4.91324691331920606875E-9,
 186  7.18466403235734541950E-11,
 187 };
 188 #endif
 189 #ifdef DEC
 190 static unsigned short BN[44] = {
 191 0040002,0045342,0113762,0004360,
 192 0137572,0120346,0172745,0144046,
 193 0037300,0151134,0123440,0117047,
 194 0136663,0025423,0014755,0046026,
 195 0036157,0177561,0027535,0046744,
 196 0135335,0161052,0071243,0146535,
 197 0034426,0154060,0164506,0135625,
 198 0133420,0005356,0100017,0151334,
 199 0032303,0066137,0024013,0046212,
 200 0131045,0016612,0066270,0047574,
 201 0027435,0177025,0060625,0116363,
 202 };
 203 static unsigned short BD[40] = {
 204 /*0040200,0000000,0000000,0000000,*/
 205 0140041,0140101,0174552,0037073,
 206 0037562,0061503,0124271,0160756,
 207 0137131,0151760,0073210,0110534,
 208 0036412,0170562,0117017,0155377,
 209 0135570,0101374,0074056,0037276,
 210 0034643,0145376,0001516,0060636,
 211 0133630,0173540,0121344,0155231,
 212 0032513,0005602,0134516,0007144,
 213 0131250,0150540,0075747,0105341,
 214 0027635,0177020,0012465,0125402,
 215 };
 216 #endif
 217 #ifdef IBMPC
 218 static unsigned short BN[44] = {
 219 0x411e,0x52fe,0x495c,0x3fe0,
 220 0xb905,0xdebc,0x541c,0xbfcf,
 221 0x13c5,0x94e4,0x1a4b,0x3fb8,
 222 0xa983,0x633d,0x6562,0xbf96,
 223 0xa9bd,0x25eb,0xffee,0x3f6d,
 224 0x79ac,0x4e54,0xbc45,0xbf3b,
 225 0xd773,0x1d28,0xdb06,0x3f02,
 226 0xfa5b,0xd001,0x015d,0xbec2,
 227 0x6991,0xe501,0x6d8b,0x3e78,
 228 0x09f0,0x4d97,0xa3b1,0xbe24,
 229 0xb39e,0xac32,0xbfc2,0x3dc3,
 230 };
 231 static unsigned short BD[40] = {
 232 /*0x0000,0x0000,0x0000,0x3ff0,*/
 233 0x47c7,0x3f2d,0x3808,0xbfe4,
 234 0x3c3e,0x7517,0x4c68,0x3fce,
 235 0x122b,0x0ed1,0x3a7e,0xbfab,
 236 0xfb60,0x53c1,0x5e2e,0x3f81,
 237 0xc7d8,0x8f05,0x105f,0xbf4f,
 238 0xcc34,0xc069,0x795f,0x3f14,
 239 0x9b53,0x145c,0x1eec,0xbed3,
 240 0xc1cd,0x5729,0x6170,0x3e89,
 241 0xf15c,0x0f7c,0x1a2c,0xbe35,
 242 0xb560,0x02a6,0xbfc2,0x3dd3,
 243 };
 244 #endif
 245 #ifdef MIEEE
 246 static unsigned short BN[44] = {
 247 0x3fe0,0x495c,0x52fe,0x411e,
 248 0xbfcf,0x541c,0xdebc,0xb905,
 249 0x3fb8,0x1a4b,0x94e4,0x13c5,
 250 0xbf96,0x6562,0x633d,0xa983,
 251 0x3f6d,0xffee,0x25eb,0xa9bd,
 252 0xbf3b,0xbc45,0x4e54,0x79ac,
 253 0x3f02,0xdb06,0x1d28,0xd773,
 254 0xbec2,0x015d,0xd001,0xfa5b,
 255 0x3e78,0x6d8b,0xe501,0x6991,
 256 0xbe24,0xa3b1,0x4d97,0x09f0,
 257 0x3dc3,0xbfc2,0xac32,0xb39e,
 258 };
 259 static unsigned short BD[40] = {
 260 /*0x3ff0,0x0000,0x0000,0x0000,*/
 261 0xbfe4,0x3808,0x3f2d,0x47c7,
 262 0x3fce,0x4c68,0x7517,0x3c3e,
 263 0xbfab,0x3a7e,0x0ed1,0x122b,
 264 0x3f81,0x5e2e,0x53c1,0xfb60,
 265 0xbf4f,0x105f,0x8f05,0xc7d8,
 266 0x3f14,0x795f,0xc069,0xcc34,
 267 0xbed3,0x1eec,0x145c,0x9b53,
 268 0x3e89,0x6170,0x5729,0xc1cd,
 269 0xbe35,0x1a2c,0x0f7c,0xf15c,
 270 0x3dd3,0xbfc2,0x02a6,0xb560,
 271 };
 272 #endif
 273
 274 /* 6.25 to infinity */
 275 #ifdef UNK
 276 static double CN[5] = {
 277 -5.90592860534773254987E-1,
 278  6.29235242724368800674E-1,
 279 -1.72858975380388136411E-1,
 280  1.64837047825189632310E-2,
 281 -4.86827613020462700845E-4,
 282 };
 283 static double CD[5] = {
 284 /* 1.00000000000000000000E0,*/
 285 -2.69820057197544900361E0,
 286  1.73270799045947845857E0,
 287 -3.93708582281939493482E-1,
 288  3.44278924041233391079E-2,
 289 -9.73655226040941223894E-4,
 290 };
 291 #endif
 292 #ifdef DEC
 293 static unsigned short CN[20] = {
 294 0140027,0030427,0176477,0074402,
 295 0040041,0012617,0112375,0162657,
 296 0137461,0000761,0074120,0135160,
 297 0036607,0004325,0117246,0115525,
 298 0135377,0036345,0064750,0047732,
 299 };
 300 static unsigned short CD[20] = {
 301 /*0040200,0000000,0000000,0000000,*/
 302 0140454,0127521,0071653,0133415,
 303 0040335,0144540,0016105,0045241,
 304 0137711,0112053,0155034,0062237,
 305 0037015,0002102,0177442,0074546,
 306 0135577,0036345,0064750,0052152,
 307 };
 308 #endif
 309 #ifdef IBMPC
 310 static unsigned short CN[20] = {
 311 0xef20,0xffa7,0xe622,0xbfe2,
 312 0xbcb6,0xf29f,0x22b1,0x3fe4,
 313 0x174e,0x2f0a,0x203e,0xbfc6,
 314 0xd36b,0xb3d4,0xe11a,0x3f90,
 315 0x09fb,0xad3d,0xe79c,0xbf3f,
 316 };
 317 static unsigned short CD[20] = {
 318 /*0x0000,0x0000,0x0000,0x3ff0,*/
 319 0x76e2,0x2e75,0x95ea,0xc005,
 320 0xa954,0x0388,0xb92c,0x3ffb,
 321 0x8c94,0x7b43,0x3285,0xbfd9,
 322 0x4f2d,0x5fe4,0xa088,0x3fa1,
 323 0x0a8d,0xad3d,0xe79c,0xbf4f,
 324 };
 325 #endif
 326 #ifdef MIEEE
 327 static unsigned short CN[20] = {
 328 0xbfe2,0xe622,0xffa7,0xef20,
 329 0x3fe4,0x22b1,0xf29f,0xbcb6,
 330 0xbfc6,0x203e,0x2f0a,0x174e,
 331 0x3f90,0xe11a,0xb3d4,0xd36b,
 332 0xbf3f,0xe79c,0xad3d,0x09fb,
 333 };
 334 static unsigned short CD[20] = {
 335 /*0x3ff0,0x0000,0x0000,0x0000,*/
 336 0xc005,0x95ea,0x2e75,0x76e2,
 337 0x3ffb,0xb92c,0x0388,0xa954,
 338 0xbfd9,0x3285,0x7b43,0x8c94,
 339 0x3fa1,0xa088,0x5fe4,0x4f2d,
 340 0xbf4f,0xe79c,0xad3d,0x0a8d,
 341 };
 342 #endif
 343
 344 #ifdef ANSIPROT
 345 extern double chbevl ( double, void *, int );
 346 extern double sqrt ( double );
 347 extern double fabs ( double );
 348 extern double polevl ( double, void *, int );
 349 extern double p1evl ( double, void *, int );
 350 #else
 351 double chbevl(), sqrt(), fabs(), polevl(), p1evl();
 352 #endif
 353 extern double PI, MACHEP;
 354
 355 double dawsn( xx )
 356 double xx;
 357 {
 358 double x, y;
 359 int sign;
 360
 361
 362 sign = 1;
 363 if( xx < 0.0 )
 364         {
 365         sign = -1;
 366         xx = -xx;
 367         }
 368
 369 if( xx < 3.25 )
 370 {
 371 x = xx*xx;
 372 y = xx * polevl( x, AN, 9 )/polevl( x, AD, 10 );
 373 return( sign * y );
 374 }
 375
 376
 377 x = 1.0/(xx*xx);
 378
 379 if( xx < 6.25 )
 380         {
 381         y = 1.0/xx + x * polevl( x, BN, 10) / (p1evl( x, BD, 10) * xx);
 382         return( sign * 0.5 * y );
 383         }
 384
 385
 386 if( xx > 1.0e9 )
 387         return( (sign * 0.5)/xx );
 388
 389 /* 6.25 to infinity */
 390 y = 1.0/xx + x * polevl( x, CN, 4) / (p1evl( x, CD, 5) * xx);
 391 return( sign * 0.5 * y );
 392 }