sici.c
Go to the documentation of this file.
1 /* sici.c
2  *
3  * Sine and cosine integrals
4  *
5  *
6  *
7  * SYNOPSIS:
8  *
9  * double x, Ci, Si, sici();
10  *
11  * sici( x, &Si, &Ci );
12  *
13  *
14  * DESCRIPTION:
15  *
16  * Evaluates the integrals
17  *
18  * x
19  * -
20  * | cos t - 1
21  * Ci(x) = eul + ln x + | --------- dt,
22  * | t
23  * -
24  * 0
25  * x
26  * -
27  * | sin t
28  * Si(x) = | ----- dt
29  * | t
30  * -
31  * 0
32  *
33  * where eul = 0.57721566490153286061 is Euler's constant.
34  * The integrals are approximated by rational functions.
35  * For x > 8 auxiliary functions f(x) and g(x) are employed
36  * such that
37  *
38  * Ci(x) = f(x) sin(x) - g(x) cos(x)
39  * Si(x) = pi/2 - f(x) cos(x) - g(x) sin(x)
40  *
41  *
42  * ACCURACY:
43  * Test interval = [0,50].
44  * Absolute error, except relative when > 1:
45  * arithmetic function # trials peak rms
46  * IEEE Si 30000 4.4e-16 7.3e-17
47  * IEEE Ci 30000 6.9e-16 5.1e-17
48  */
49 
50 /*
51  * Cephes Math Library Release 2.1: January, 1989
52  * Copyright 1984, 1987, 1989 by Stephen L. Moshier
53  * Direct inquiries to 30 Frost Street, Cambridge, MA 02140
54  */
55 
56 #include "mconf.h"
57 
58 static double SN[] = {
59  -8.39167827910303881427E-11,
60  4.62591714427012837309E-8,
61  -9.75759303843632795789E-6,
62  9.76945438170435310816E-4,
63  -4.13470316229406538752E-2,
64  1.00000000000000000302E0,
65 };
66 
67 static double SD[] = {
68  2.03269266195951942049E-12,
69  1.27997891179943299903E-9,
70  4.41827842801218905784E-7,
71  9.96412122043875552487E-5,
72  1.42085239326149893930E-2,
73  9.99999999999999996984E-1,
74 };
75 
76 static double CN[] = {
77  2.02524002389102268789E-11,
78  -1.35249504915790756375E-8,
79  3.59325051419993077021E-6,
80  -4.74007206873407909465E-4,
81  2.89159652607555242092E-2,
82  -1.00000000000000000080E0,
83 };
84 
85 static double CD[] = {
86  4.07746040061880559506E-12,
87  3.06780997581887812692E-9,
88  1.23210355685883423679E-6,
89  3.17442024775032769882E-4,
90  5.10028056236446052392E-2,
91  4.00000000000000000080E0,
92 };
93 
94 static double FN4[] = {
95  4.23612862892216586994E0,
96  5.45937717161812843388E0,
97  1.62083287701538329132E0,
98  1.67006611831323023771E-1,
99  6.81020132472518137426E-3,
100  1.08936580650328664411E-4,
101  5.48900223421373614008E-7,
102 };
103 
104 static double FD4[] = {
105  /* 1.00000000000000000000E0, */
106  8.16496634205391016773E0,
107  7.30828822505564552187E0,
108  1.86792257950184183883E0,
109  1.78792052963149907262E-1,
110  7.01710668322789753610E-3,
111  1.10034357153915731354E-4,
112  5.48900252756255700982E-7,
113 };
114 
115 static double FN8[] = {
116  4.55880873470465315206E-1,
117  7.13715274100146711374E-1,
118  1.60300158222319456320E-1,
119  1.16064229408124407915E-2,
120  3.49556442447859055605E-4,
121  4.86215430826454749482E-6,
122  3.20092790091004902806E-8,
123  9.41779576128512936592E-11,
124  9.70507110881952024631E-14,
125 };
126 
127 static double FD8[] = {
128  /* 1.00000000000000000000E0, */
129  9.17463611873684053703E-1,
130  1.78685545332074536321E-1,
131  1.22253594771971293032E-2,
132  3.58696481881851580297E-4,
133  4.92435064317881464393E-6,
134  3.21956939101046018377E-8,
135  9.43720590350276732376E-11,
136  9.70507110881952025725E-14,
137 };
138 
139 static double GN4[] = {
140  8.71001698973114191777E-2,
141  6.11379109952219284151E-1,
142  3.97180296392337498885E-1,
143  7.48527737628469092119E-2,
144  5.38868681462177273157E-3,
145  1.61999794598934024525E-4,
146  1.97963874140963632189E-6,
147  7.82579040744090311069E-9,
148 };
149 
150 static double GD4[] = {
151  /* 1.00000000000000000000E0, */
152  1.64402202413355338886E0,
153  6.66296701268987968381E-1,
154  9.88771761277688796203E-2,
155  6.22396345441768420760E-3,
156  1.73221081474177119497E-4,
157  2.02659182086343991969E-6,
158  7.82579218933534490868E-9,
159 };
160 
161 static double GN8[] = {
162  6.97359953443276214934E-1,
163  3.30410979305632063225E-1,
164  3.84878767649974295920E-2,
165  1.71718239052347903558E-3,
166  3.48941165502279436777E-5,
167  3.47131167084116673800E-7,
168  1.70404452782044526189E-9,
169  3.85945925430276600453E-12,
170  3.14040098946363334640E-15,
171 };
172 
173 static double GD8[] = {
174  /* 1.00000000000000000000E0, */
175  1.68548898811011640017E0,
176  4.87852258695304967486E-1,
177  4.67913194259625806320E-2,
178  1.90284426674399523638E-3,
179  3.68475504442561108162E-5,
180  3.57043223443740838771E-7,
181  1.72693748966316146736E-9,
182  3.87830166023954706752E-12,
183  3.14040098946363335242E-15,
184 };
185 
186 extern double MACHEP;
187 
188 
189 int sici(double x, double *si, double *ci)
190 {
191  double z, c, s, f, g;
192  short sign;
193 
194  if (x < 0.0) {
195  sign = -1;
196  x = -x;
197  }
198  else
199  sign = 0;
200 
201 
202  if (x == 0.0) {
203  *si = 0.0;
204  *ci = -INFINITY;
205  return (0);
206  }
207 
208 
209  if (x > 1.0e9) {
210  if (cephes_isinf(x)) {
211  if (sign == -1) {
212  *si = -M_PI_2;
213  *ci = NAN;
214  }
215  else {
216  *si = M_PI_2;
217  *ci = 0;
218  }
219  return 0;
220  }
221  *si = M_PI_2 - cos(x) / x;
222  *ci = sin(x) / x;
223  }
224 
225 
226 
227  if (x > 4.0)
228  goto asympt;
229 
230  z = x * x;
231  s = x * polevl(z, SN, 5) / polevl(z, SD, 5);
232  c = z * polevl(z, CN, 5) / polevl(z, CD, 5);
233 
234  if (sign)
235  s = -s;
236  *si = s;
237  *ci = SCIPY_EULER + log(x) + c; /* real part if x < 0 */
238  return (0);
239 
240 
241 
242  /* The auxiliary functions are:
243  *
244  *
245  * *si = *si - M_PI_2;
246  * c = cos(x);
247  * s = sin(x);
248  *
249  * t = *ci * s - *si * c;
250  * a = *ci * c + *si * s;
251  *
252  * *si = t;
253  * *ci = -a;
254  */
255 
256 
257  asympt:
258 
259  s = sin(x);
260  c = cos(x);
261  z = 1.0 / (x * x);
262  if (x < 8.0) {
263  f = polevl(z, FN4, 6) / (x * p1evl(z, FD4, 7));
264  g = z * polevl(z, GN4, 7) / p1evl(z, GD4, 7);
265  }
266  else {
267  f = polevl(z, FN8, 8) / (x * p1evl(z, FD8, 8));
268  g = z * polevl(z, GN8, 8) / p1evl(z, GD8, 9);
269  }
270  *si = M_PI_2 - f * c - g * s;
271  if (sign)
272  *si = -(*si);
273  *ci = f * s - g * c;
274 
275  return (0);
276 }
s
RealScalar s
Definition: level1_cplx_impl.h:126
GD8
static double GD8[]
Definition: sici.c:173
ceres::sin
Jet< T, N > sin(const Jet< T, N > &f)
Definition: jet.h:439
c
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
CD
static double CD[]
Definition: sici.c:85
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
log
const EIGEN_DEVICE_FUNC LogReturnType log() const
Definition: ArrayCwiseUnaryOps.h:128
cephes_isinf
#define cephes_isinf(x)
Definition: mconf.h:100
FD8
static double FD8[]
Definition: sici.c:127
GN8
static double GN8[]
Definition: sici.c:161
CN
static double CN[]
Definition: sici.c:76
SCIPY_EULER
#define SCIPY_EULER
Definition: mconf.h:128
ceres::cos
Jet< T, N > cos(const Jet< T, N > &f)
Definition: jet.h:426
GN4
static double GN4[]
Definition: sici.c:139
p1evl
static double p1evl(double x, const double coef[], int N)
Definition: polevl.h:97
FN4
static double FN4[]
Definition: sici.c:94
polevl
static double polevl(double x, const double coef[], int N)
Definition: polevl.h:65
SD
static double SD[]
Definition: sici.c:67
pybind_wrapper_test_script.z
z
Definition: pybind_wrapper_test_script.py:61
GD4
static double GD4[]
Definition: sici.c:150
g
void g(const string &key, int i)
Definition: testBTree.cpp:41
M_PI_2
#define M_PI_2
Definition: mconf.h:118
tree::f
Point2(* f)(const Point3 &, OptionalJacobian< 2, 3 >)
Definition: testExpression.cpp:218
mconf.h
sici
int sici(double x, double *si, double *ci)
Definition: sici.c:189
FD4
static double FD4[]
Definition: sici.c:104
MACHEP
double MACHEP
Definition: const.c:54
SN
static double SN[]
Definition: sici.c:58
FN8
static double FN8[]
Definition: sici.c:115


gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:35:16