j0.c
Go to the documentation of this file.
1 /* j0.c
2  *
3  * Bessel function of order zero
4  *
5  *
6  *
7  * SYNOPSIS:
8  *
9  * double x, y, j0();
10  *
11  * y = j0( x );
12  *
13  *
14  *
15  * DESCRIPTION:
16  *
17  * Returns Bessel function of order zero of the argument.
18  *
19  * The domain is divided into the intervals [0, 5] and
20  * (5, infinity). In the first interval the following rational
21  * approximation is used:
22  *
23  *
24  * 2 2
25  * (w - r ) (w - r ) P (w) / Q (w)
26  * 1 2 3 8
27  *
28  * 2
29  * where w = x and the two r's are zeros of the function.
30  *
31  * In the second interval, the Hankel asymptotic expansion
32  * is employed with two rational functions of degree 6/6
33  * and 7/7.
34  *
35  *
36  *
37  * ACCURACY:
38  *
39  * Absolute error:
40  * arithmetic domain # trials peak rms
41  * IEEE 0, 30 60000 4.2e-16 1.1e-16
42  *
43  */
44  /* y0.c
45  *
46  * Bessel function of the second kind, order zero
47  *
48  *
49  *
50  * SYNOPSIS:
51  *
52  * double x, y, y0();
53  *
54  * y = y0( x );
55  *
56  *
57  *
58  * DESCRIPTION:
59  *
60  * Returns Bessel function of the second kind, of order
61  * zero, of the argument.
62  *
63  * The domain is divided into the intervals [0, 5] and
64  * (5, infinity). In the first interval a rational approximation
65  * R(x) is employed to compute
66  * y0(x) = R(x) + 2 * log(x) * j0(x) / M_PI.
67  * Thus a call to j0() is required.
68  *
69  * In the second interval, the Hankel asymptotic expansion
70  * is employed with two rational functions of degree 6/6
71  * and 7/7.
72  *
73  *
74  *
75  * ACCURACY:
76  *
77  * Absolute error, when y0(x) < 1; else relative error:
78  *
79  * arithmetic domain # trials peak rms
80  * IEEE 0, 30 30000 1.3e-15 1.6e-16
81  *
82  */
83 
84 /*
85  * Cephes Math Library Release 2.8: June, 2000
86  * Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
87  */
88 
89 /* Note: all coefficients satisfy the relative error criterion
90  * except YP, YQ which are designed for absolute error. */
91 
92 #include "mconf.h"
93 
94 static double PP[7] = {
95  7.96936729297347051624E-4,
96  8.28352392107440799803E-2,
97  1.23953371646414299388E0,
98  5.44725003058768775090E0,
99  8.74716500199817011941E0,
100  5.30324038235394892183E0,
101  9.99999999999999997821E-1,
102 };
103 
104 static double PQ[7] = {
105  9.24408810558863637013E-4,
106  8.56288474354474431428E-2,
107  1.25352743901058953537E0,
108  5.47097740330417105182E0,
109  8.76190883237069594232E0,
110  5.30605288235394617618E0,
111  1.00000000000000000218E0,
112 };
113 
114 static double QP[8] = {
115  -1.13663838898469149931E-2,
116  -1.28252718670509318512E0,
117  -1.95539544257735972385E1,
118  -9.32060152123768231369E1,
119  -1.77681167980488050595E2,
120  -1.47077505154951170175E2,
121  -5.14105326766599330220E1,
122  -6.05014350600728481186E0,
123 };
124 
125 static double QQ[7] = {
126  /* 1.00000000000000000000E0, */
127  6.43178256118178023184E1,
128  8.56430025976980587198E2,
129  3.88240183605401609683E3,
130  7.24046774195652478189E3,
131  5.93072701187316984827E3,
132  2.06209331660327847417E3,
133  2.42005740240291393179E2,
134 };
135 
136 static double YP[8] = {
137  1.55924367855235737965E4,
138  -1.46639295903971606143E7,
139  5.43526477051876500413E9,
140  -9.82136065717911466409E11,
141  8.75906394395366999549E13,
142  -3.46628303384729719441E15,
143  4.42733268572569800351E16,
144  -1.84950800436986690637E16,
145 };
146 
147 static double YQ[7] = {
148  /* 1.00000000000000000000E0, */
149  1.04128353664259848412E3,
150  6.26107330137134956842E5,
151  2.68919633393814121987E8,
152  8.64002487103935000337E10,
153  2.02979612750105546709E13,
154  3.17157752842975028269E15,
155  2.50596256172653059228E17,
156 };
157 
158 /* 5.783185962946784521175995758455807035071 */
159 static double DR1 = 5.78318596294678452118E0;
160 
161 /* 30.47126234366208639907816317502275584842 */
162 static double DR2 = 3.04712623436620863991E1;
163 
164 static double RP[4] = {
165  -4.79443220978201773821E9,
166  1.95617491946556577543E12,
167  -2.49248344360967716204E14,
168  9.70862251047306323952E15,
169 };
170 
171 static double RQ[8] = {
172  /* 1.00000000000000000000E0, */
173  4.99563147152651017219E2,
174  1.73785401676374683123E5,
175  4.84409658339962045305E7,
176  1.11855537045356834862E10,
177  2.11277520115489217587E12,
178  3.10518229857422583814E14,
179  3.18121955943204943306E16,
180  1.71086294081043136091E18,
181 };
182 
183 extern double SQ2OPI;
184 
185 double j0(double x)
186 {
187  double w, z, p, q, xn;
188 
189  if (x < 0)
190  x = -x;
191 
192  if (x <= 5.0) {
193  z = x * x;
194  if (x < 1.0e-5)
195  return (1.0 - z / 4.0);
196 
197  p = (z - DR1) * (z - DR2);
198  p = p * polevl(z, RP, 3) / p1evl(z, RQ, 8);
199  return (p);
200  }
201 
202  w = 5.0 / x;
203  q = 25.0 / (x * x);
204  p = polevl(q, PP, 6) / polevl(q, PQ, 6);
205  q = polevl(q, QP, 7) / p1evl(q, QQ, 7);
206  xn = x - M_PI_4;
207  p = p * cos(xn) - w * q * sin(xn);
208  return (p * SQ2OPI / sqrt(x));
209 }
210 
211 /* y0() 2 */
212 /* Bessel function of second kind, order zero */
213 
214 /* Rational approximation coefficients YP[], YQ[] are used here.
215  * The function computed is y0(x) - 2 * log(x) * j0(x) / M_PI,
216  * whose value at x = 0 is 2 * ( log(0.5) + EUL ) / M_PI
217  * = 0.073804295108687225.
218  */
219 
220 double y0(double x)
221 {
222  double w, z, p, q, xn;
223 
224  if (x <= 5.0) {
225  if (x == 0.0) {
227  return -INFINITY;
228  }
229  else if (x < 0.0) {
230  sf_error("y0", SF_ERROR_DOMAIN, NULL);
231  return NAN;
232  }
233  z = x * x;
234  w = polevl(z, YP, 7) / p1evl(z, YQ, 7);
235  w += M_2_PI * log(x) * j0(x);
236  return (w);
237  }
238 
239  w = 5.0 / x;
240  z = 25.0 / (x * x);
241  p = polevl(z, PP, 6) / polevl(z, PQ, 6);
242  q = polevl(z, QP, 7) / p1evl(z, QQ, 7);
243  xn = x - M_PI_4;
244  p = p * sin(xn) + w * q * cos(xn);
245  return (p * SQ2OPI / sqrt(x));
246 }
w
RowVector3d w
Definition: Matrix_resize_int.cpp:3
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
M_2_PI
#define M_2_PI
Definition: mconf.h:121
ceres::sin
Jet< T, N > sin(const Jet< T, N > &f)
Definition: jet.h:439
DR1
static double DR1
Definition: j0.c:159
PP
static double PP[7]
Definition: j0.c:94
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
log
const EIGEN_DEVICE_FUNC LogReturnType log() const
Definition: ArrayCwiseUnaryOps.h:128
QQ
static double QQ[7]
Definition: j0.c:125
y0
double y0(double x)
Definition: j0.c:220
ceres::cos
Jet< T, N > cos(const Jet< T, N > &f)
Definition: jet.h:426
p1evl
static double p1evl(double x, const double coef[], int N)
Definition: polevl.h:97
DR2
static double DR2
Definition: j0.c:162
PQ
static double PQ[7]
Definition: j0.c:104
Eigen::numext::q
EIGEN_DEVICE_FUNC const Scalar & q
Definition: SpecialFunctionsImpl.h:1984
QP
static double QP[8]
Definition: j0.c:114
M_PI_4
#define M_PI_4
Definition: mconf.h:119
polevl
static double polevl(double x, const double coef[], int N)
Definition: polevl.h:65
YP
static double YP[8]
Definition: j0.c:136
pybind_wrapper_test_script.z
z
Definition: pybind_wrapper_test_script.py:61
RQ
static double RQ[8]
Definition: j0.c:171
j0
double j0(double x)
Definition: j0.c:185
YQ
static double YQ[7]
Definition: j0.c:147
mconf.h
p
float * p
Definition: Tutorial_Map_using.cpp:9
sf_error
void sf_error(const char *func_name, sf_error_t code, const char *fmt,...)
Definition: sf_error.c:41
RP
static double RP[4]
Definition: j0.c:164
SF_ERROR_SINGULAR
@ SF_ERROR_SINGULAR
Definition: sf_error.h:10
NULL
#define NULL
Definition: ccolamd.c:609
SF_ERROR_DOMAIN
@ SF_ERROR_DOMAIN
Definition: sf_error.h:16
ceres::sqrt
Jet< T, N > sqrt(const Jet< T, N > &f)
Definition: jet.h:418
SQ2OPI
double SQ2OPI
Definition: const.c:65


gtsam
Author(s):
autogenerated on Thu Jun 13 2024 03:02:47