dawsn.c
Go to the documentation of this file.
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  *
37  *
38  */
39 
40 /* dawsn.c */
41 
42 
43 /*
44  * Cephes Math Library Release 2.1: January, 1989
45  * Copyright 1984, 1987, 1989 by Stephen L. Moshier
46  * Direct inquiries to 30 Frost Street, Cambridge, MA 02140
47  */
48 
49 #include "mconf.h"
50 /* Dawson's integral, interval 0 to 3.25 */
51 static double AN[10] = {
52  1.13681498971755972054E-11,
53  8.49262267667473811108E-10,
54  1.94434204175553054283E-8,
55  9.53151741254484363489E-7,
56  3.07828309874913200438E-6,
57  3.52513368520288738649E-4,
58  -8.50149846724410912031E-4,
59  4.22618223005546594270E-2,
60  -9.17480371773452345351E-2,
61  9.99999999999999994612E-1,
62 };
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 
78 /* interval 3.25 to 6.25 */
79 static double BN[11] = {
80  5.08955156417900903354E-1,
81  -2.44754418142697847934E-1,
82  9.41512335303534411857E-2,
83  -2.18711255142039025206E-2,
84  3.66207612329569181322E-3,
85  -4.23209114460388756528E-4,
86  3.59641304793896631888E-5,
87  -2.14640351719968974225E-6,
88  9.10010780076391431042E-8,
89  -2.40274520828250956942E-9,
90  3.59233385440928410398E-11,
91 };
92 
93 static double BD[10] = {
94  /* 1.00000000000000000000E0, */
95  -6.31839869873368190192E-1,
96  2.36706788228248691528E-1,
97  -5.31806367003223277662E-2,
98  8.48041718586295374409E-3,
99  -9.47996768486665330168E-4,
100  7.81025592944552338085E-5,
101  -4.55875153252442634831E-6,
102  1.89100358111421846170E-7,
103  -4.91324691331920606875E-9,
104  7.18466403235734541950E-11,
105 };
106 
107 /* 6.25 to infinity */
108 static double CN[5] = {
109  -5.90592860534773254987E-1,
110  6.29235242724368800674E-1,
111  -1.72858975380388136411E-1,
112  1.64837047825189632310E-2,
113  -4.86827613020462700845E-4,
114 };
115 
116 static double CD[5] = {
117  /* 1.00000000000000000000E0, */
118  -2.69820057197544900361E0,
119  1.73270799045947845857E0,
120  -3.93708582281939493482E-1,
121  3.44278924041233391079E-2,
122  -9.73655226040941223894E-4,
123 };
124 
125 extern double MACHEP;
126 
127 double dawsn(double xx)
128 {
129  double x, y;
130  int sign;
131 
132 
133  sign = 1;
134  if (xx < 0.0) {
135  sign = -1;
136  xx = -xx;
137  }
138 
139  if (xx < 3.25) {
140  x = xx * xx;
141  y = xx * polevl(x, AN, 9) / polevl(x, AD, 10);
142  return (sign * y);
143  }
144 
145 
146  x = 1.0 / (xx * xx);
147 
148  if (xx < 6.25) {
149  y = 1.0 / xx + x * polevl(x, BN, 10) / (p1evl(x, BD, 10) * xx);
150  return (sign * 0.5 * y);
151  }
152 
153 
154  if (xx > 1.0e9)
155  return ((sign * 0.5) / xx);
156 
157  /* 6.25 to infinity */
158  y = 1.0 / xx + x * polevl(x, CN, 4) / (p1evl(x, CD, 5) * xx);
159  return (sign * 0.5 * y);
160 }
BN
static double BN[11]
Definition: dawsn.c:79
x
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
Definition: gnuplot_common_settings.hh:12
sign
const EIGEN_DEVICE_FUNC SignReturnType sign() const
Definition: ArrayCwiseUnaryOps.h:219
CN
static double CN[5]
Definition: dawsn.c:108
dawsn
double dawsn(double xx)
Definition: dawsn.c:127
p1evl
static double p1evl(double x, const double coef[], int N)
Definition: polevl.h:97
CD
static double CD[5]
Definition: dawsn.c:116
polevl
static double polevl(double x, const double coef[], int N)
Definition: polevl.h:65
MACHEP
double MACHEP
Definition: const.c:54
BD
static double BD[10]
Definition: dawsn.c:93
y
Scalar * y
Definition: level1_cplx_impl.h:124
mconf.h
AN
static double AN[10]
Definition: dawsn.c:51
AD
static double AD[11]
Definition: dawsn.c:64


gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:32:22