1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2016 Gael Guennebaud <gael.guennebaud@inria.fr>
5 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9
10 #include "main.h"
11 #include "../Eigen/SpecialFunctions"
12
13 template<typename X, typename Y>
verify_component_wise(const X & x,const Y & y)14 void verify_component_wise(const X& x, const Y& y)
15 {
16 for(Index i=0; i<x.size(); ++i)
17 {
18 if((numext::isfinite)(y(i))) {
19 VERIFY_IS_APPROX( x(i), y(i) );
20 }
21 else if((numext::isnan)(y(i)))
22 VERIFY((numext::isnan)(x(i)));
23 else
24 VERIFY_IS_EQUAL( x(i), y(i) );
25 }
26 }
27
array_bessel_functions()28 template<typename ArrayType> void array_bessel_functions()
29 {
30 // Test Bessel function i0. Reference results obtained with SciPy.
31 {
32 ArrayType x(21);
33 ArrayType expected(21);
34 ArrayType res(21);
35
36 x << -20.0, -18.0, -16.0, -14.0, -12.0, -10.0, -8.0, -6.0, -4.0, -2.0, 0.0,
37 2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0;
38
39 expected << 4.35582826e+07, 6.21841242e+06, 8.93446228e+05, 1.29418563e+05,
40 1.89489253e+04, 2.81571663e+03, 4.27564116e+02, 6.72344070e+01,
41 1.13019220e+01, 2.27958530e+00, 1.00000000e+00, 2.27958530e+00,
42 1.13019220e+01, 6.72344070e+01, 4.27564116e+02, 2.81571663e+03,
43 1.89489253e+04, 1.29418563e+05, 8.93446228e+05, 6.21841242e+06,
44 4.35582826e+07;
45
46 CALL_SUBTEST(res = bessel_i0(x);
47 verify_component_wise(res, expected););
48 }
49
50 // Test Bessel function i0e. Reference results obtained with SciPy.
51 {
52 ArrayType x(21);
53 ArrayType expected(21);
54 ArrayType res(21);
55
56 x << -20.0, -18.0, -16.0, -14.0, -12.0, -10.0, -8.0, -6.0, -4.0, -2.0, 0.0,
57 2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0;
58
59 expected << 0.0897803118848, 0.0947062952128, 0.100544127361,
60 0.107615251671, 0.116426221213, 0.127833337163, 0.143431781857,
61 0.16665743264, 0.207001921224, 0.308508322554, 1.0, 0.308508322554,
62 0.207001921224, 0.16665743264, 0.143431781857, 0.127833337163,
63 0.116426221213, 0.107615251671, 0.100544127361, 0.0947062952128,
64 0.0897803118848;
65
66 CALL_SUBTEST(res = bessel_i0e(x);
67 verify_component_wise(res, expected););
68 }
69
70 // Test Bessel function i1. Reference results obtained with SciPy.
71 {
72 ArrayType x(21);
73 ArrayType expected(21);
74 ArrayType res(21);
75
76 x << -20.0, -18.0, -16.0, -14.0, -12.0, -10.0, -8.0, -6.0, -4.0, -2.0, 0.0,
77 2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0;
78
79 expected << -4.24549734e+07, -6.04313324e+06, -8.65059436e+05, -1.24707259e+05,
80 -1.81413488e+04, -2.67098830e+03, -3.99873137e+02, -6.13419368e+01,
81 -9.75946515e+00, -1.59063685e+00, 0.00000000e+00, 1.59063685e+00,
82 9.75946515e+00, 6.13419368e+01, 3.99873137e+02, 2.67098830e+03,
83 1.81413488e+04, 1.24707259e+05, 8.65059436e+05, 6.04313324e+06,
84 4.24549734e+07;
85
86 CALL_SUBTEST(res = bessel_i1(x);
87 verify_component_wise(res, expected););
88 }
89
90 // Test Bessel function i1e. Reference results obtained with SciPy.
91 {
92 ArrayType x(21);
93 ArrayType expected(21);
94 ArrayType res(21);
95
96 x << -20.0, -18.0, -16.0, -14.0, -12.0, -10.0, -8.0, -6.0, -4.0, -2.0, 0.0,
97 2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0;
98
99 expected << -0.0875062221833, -0.092036796872, -0.0973496147565,
100 -0.103697667463, -0.11146429929, -0.121262681384, -0.134142493293,
101 -0.152051459309, -0.178750839502, -0.215269289249, 0.0, 0.215269289249,
102 0.178750839502, 0.152051459309, 0.134142493293, 0.121262681384,
103 0.11146429929, 0.103697667463, 0.0973496147565, 0.092036796872,
104 0.0875062221833;
105
106 CALL_SUBTEST(res = bessel_i1e(x);
107 verify_component_wise(res, expected););
108 }
109
110 // Test Bessel function j0. Reference results obtained with SciPy.
111 {
112 ArrayType x(77);
113 ArrayType expected(77);
114 ArrayType res(77);
115
116 x << -38., -37., -36., -35., -34., -33., -32., -31., -30.,
117 -29., -28., -27., -26., -25., -24., -23., -22., -21., -20., -19.,
118 -18., -17., -16., -15., -14., -13., -12., -11., -10., -9., -8.,
119 -7., -6., -5., -4., -3., -2., -1., 0., 1., 2., 3.,
120 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14.,
121 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25.,
122 26., 27., 28., 29., 30., 31., 32., 33., 34., 35., 36.,
123 37., 38.;
124
125 expected << 0.11433274, 0.01086237, -0.10556738,
126 -0.12684568, -0.03042119, 0.09727067, 0.13807901, 0.05120815,
127 -0.08636798, -0.14784876, -0.07315701, 0.07274192, 0.15599932,
128 0.09626678, -0.05623027, -0.16241278, -0.12065148, 0.03657907,
129 0.16702466, 0.14662944, -0.01335581, -0.16985425, -0.17489907,
130 -0.01422447, 0.17107348, 0.2069261 , 0.04768931, -0.1711903 ,
131 -0.24593576, -0.09033361, 0.17165081, 0.30007927, 0.15064526,
132 -0.17759677, -0.39714981, -0.26005195, 0.22389078, 0.76519769,
133 1. , 0.76519769, 0.22389078, -0.26005195, -0.39714981,
134 -0.17759677, 0.15064526, 0.30007927, 0.17165081, -0.09033361,
135 -0.24593576, -0.1711903 , 0.04768931, 0.2069261 , 0.17107348,
136 -0.01422447, -0.17489907, -0.16985425, -0.01335581, 0.14662944,
137 0.16702466, 0.03657907, -0.12065148, -0.16241278, -0.05623027,
138 0.09626678, 0.15599932, 0.07274192, -0.07315701, -0.14784876,
139 -0.08636798, 0.05120815, 0.13807901, 0.09727067, -0.03042119,
140 -0.12684568, -0.10556738, 0.01086237, 0.11433274;
141
142 CALL_SUBTEST(res = bessel_j0(x);
143 verify_component_wise(res, expected););
144 }
145
146 // Test Bessel function j1. Reference results obtained with SciPy.
147 {
148 ArrayType x(81);
149 ArrayType expected(81);
150 ArrayType res(81);
151
152 x << -40., -39., -38., -37., -36., -35., -34., -33., -32., -31., -30.,
153 -29., -28., -27., -26., -25., -24., -23., -22., -21., -20., -19.,
154 -18., -17., -16., -15., -14., -13., -12., -11., -10., -9., -8.,
155 -7., -6., -5., -4., -3., -2., -1., 0., 1., 2., 3.,
156 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14.,
157 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25.,
158 26., 27., 28., 29., 30., 31., 32., 33., 34., 35., 36.,
159 37., 38., 39., 40.;
160
161 expected << -0.12603832, -0.0640561 , 0.05916189, 0.13058004, 0.08232981,
162 -0.04399094, -0.13297118, -0.10061965, 0.02658903, 0.13302432,
163 0.11875106, -0.0069342 , -0.13055149, -0.13658472, -0.01504573,
164 0.12535025, 0.15403807, 0.03951932, -0.11717779, -0.17112027,
165 -0.06683312, 0.10570143, 0.18799489, 0.09766849, -0.09039718,
166 -0.20510404, -0.13337515, 0.07031805, 0.2234471 , 0.1767853 ,
167 -0.04347275, -0.24531179, -0.23463635, 0.00468282, 0.27668386,
168 0.32757914, 0.06604333, -0.33905896, -0.57672481, -0.44005059,
169 0. , 0.44005059, 0.57672481, 0.33905896, -0.06604333,
170 -0.32757914, -0.27668386, -0.00468282, 0.23463635, 0.24531179,
171 0.04347275, -0.1767853 , -0.2234471 , -0.07031805, 0.13337515,
172 0.20510404, 0.09039718, -0.09766849, -0.18799489, -0.10570143,
173 0.06683312, 0.17112027, 0.11717779, -0.03951932, -0.15403807,
174 -0.12535025, 0.01504573, 0.13658472, 0.13055149, 0.0069342 ,
175 -0.11875106, -0.13302432, -0.02658903, 0.10061965, 0.13297118,
176 0.04399094, -0.08232981, -0.13058004, -0.05916189, 0.0640561 ,
177 0.12603832;
178
179 CALL_SUBTEST(res = bessel_j1(x);
180 verify_component_wise(res, expected););
181 }
182 // Test Bessel function k0e. Reference results obtained with SciPy.
183 {
184 ArrayType x(42);
185 ArrayType expected(42);
186 ArrayType res(42);
187
188 x << 0.25, 0.5, 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12.,
189 13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25.,
190 26., 27., 28., 29., 30., 31., 32., 33., 34., 35., 36., 37., 38.,
191 39., 40.;
192
193 expected << 1.97933385, 1.52410939, 1.14446308, 0.84156822,
194 0.6977616 , 0.60929767, 0.54780756, 0.50186313, 0.4658451 ,
195 0.43662302, 0.41229555, 0.39163193, 0.3737955 , 0.35819488,
196 0.34439865, 0.33208364, 0.32100235, 0.31096159, 0.30180802,
197 0.29341821, 0.28569149, 0.27854488, 0.2719092 , 0.26572635,
198 0.25994703, 0.25452917, 0.2494366 , 0.24463801, 0.24010616,
199 0.23581722, 0.23175022, 0.22788667, 0.22421014, 0.22070602,
200 0.21736123, 0.21416406, 0.21110397, 0.20817141, 0.20535778,
201 0.20265524, 0.20005668, 0.19755558;
202
203 CALL_SUBTEST(res = bessel_k0e(x);
204 verify_component_wise(res, expected););
205 }
206
207 // Test Bessel function k0. Reference results obtained with SciPy.
208 {
209 ArrayType x(42);
210 ArrayType expected(42);
211 ArrayType res(42);
212
213 x << 0.25, 0.5, 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12.,
214 13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25.,
215 26., 27., 28., 29., 30., 31., 32., 33., 34., 35., 36., 37., 38.,
216 39., 40.;
217
218 expected << 1.54150675, 0.92441907, 4.21024438e-01, 1.13893873e-01,
219 3.47395044e-02, 1.11596761e-02, 3.69109833e-03, 1.24399433e-03,
220 4.24795742e-04, 1.46470705e-04, 5.08813130e-05, 1.77800623e-05,
221 6.24302055e-06, 2.20082540e-06, 7.78454386e-07, 2.76137082e-07,
222 9.81953648e-08, 3.49941166e-08, 1.24946640e-08, 4.46875334e-09,
223 1.60067129e-09, 5.74123782e-10, 2.06176797e-10, 7.41235161e-11,
224 2.66754511e-11, 9.60881878e-12, 3.46416156e-12, 1.24987740e-12,
225 4.51286453e-13, 1.63053459e-13, 5.89495073e-14, 2.13247750e-14,
226 7.71838266e-15, 2.79505752e-15, 1.01266123e-15, 3.67057597e-16,
227 1.33103515e-16, 4.82858338e-17, 1.75232770e-17, 6.36161716e-18,
228 2.31029936e-18, 8.39286110e-19;
229
230 CALL_SUBTEST(res = bessel_k0(x);
231 verify_component_wise(res, expected););
232 }
233
234 // Test Bessel function k0e. Reference results obtained with SciPy.
235 {
236 ArrayType x(42);
237 ArrayType expected(42);
238 ArrayType res(42);
239
240 x << 0.25, 0.5, 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12.,
241 13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25.,
242 26., 27., 28., 29., 30., 31., 32., 33., 34., 35., 36., 37., 38.,
243 39., 40.;
244
245 expected << 1.97933385, 1.52410939, 1.14446308, 0.84156822,
246 0.6977616 , 0.60929767, 0.54780756, 0.50186313,
247 0.4658451 , 0.43662302, 0.41229555, 0.39163193,
248 0.3737955 , 0.35819488, 0.34439865, 0.33208364,
249 0.32100235, 0.31096159, 0.30180802, 0.29341821,
250 0.28569149, 0.27854488, 0.2719092 , 0.26572635,
251 0.25994703, 0.25452917, 0.2494366 , 0.24463801,
252 0.24010616, 0.23581722, 0.23175022, 0.22788667,
253 0.22421014, 0.22070602, 0.21736123, 0.21416406,
254 0.21110397, 0.20817141, 0.20535778, 0.20265524,
255 0.20005668, 0.19755558;
256
257 CALL_SUBTEST(res = bessel_k0e(x);
258 verify_component_wise(res, expected););
259 }
260
261 // Test Bessel function k1. Reference results obtained with SciPy.
262 {
263 ArrayType x(42);
264 ArrayType expected(42);
265 ArrayType res(42);
266
267 x << 0.25, 0.5, 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12.,
268 13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25.,
269 26., 27., 28., 29., 30., 31., 32., 33., 34., 35., 36., 37., 38.,
270 39., 40.;
271
272 expected << 3.74702597, 1.65644112, 6.01907230e-01, 1.39865882e-01,
273 4.01564311e-02, 1.24834989e-02, 4.04461345e-03, 1.34391972e-03,
274 4.54182487e-04, 1.55369212e-04, 5.36370164e-05, 1.86487735e-05,
275 6.52086067e-06, 2.29075746e-06, 8.07858841e-07, 2.85834365e-07,
276 1.01417294e-07, 3.60715712e-08, 1.28570417e-08, 4.59124963e-09,
277 1.64226697e-09, 5.88305797e-10, 2.11029922e-10, 7.57898116e-11,
278 2.72493059e-11, 9.80699893e-12, 3.53277807e-12, 1.27369078e-12,
279 4.59568940e-13, 1.65940011e-13, 5.99574032e-14, 2.16773200e-14,
280 7.84189960e-15, 2.83839927e-15, 1.02789171e-15, 3.72416929e-16,
281 1.34991783e-16, 4.89519373e-17, 1.77585196e-17, 6.44478588e-18,
282 2.33973340e-18, 8.49713195e-19;
283
284 CALL_SUBTEST(res = bessel_k1(x);
285 verify_component_wise(res, expected););
286 }
287
288 // Test Bessel function k1e. Reference results obtained with SciPy.
289 {
290 ArrayType x(42);
291 ArrayType expected(42);
292 ArrayType res(42);
293
294 x << 0.25, 0.5, 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12.,
295 13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25.,
296 26., 27., 28., 29., 30., 31., 32., 33., 34., 35., 36., 37., 38.,
297 39., 40.;
298
299 expected << 4.81127659, 2.73100971, 1.63615349, 1.03347685,
300 0.80656348, 0.68157595, 0.60027386, 0.54217591,
301 0.49807158, 0.46314909, 0.43462525, 0.41076657,
302 0.39043094, 0.37283175, 0.35740757, 0.34374563,
303 0.33153489, 0.32053597, 0.31056123, 0.30146131,
304 0.29311559, 0.2854255 , 0.27830958, 0.27169987,
305 0.26553913, 0.25977879, 0.25437733, 0.249299 ,
306 0.24451285, 0.23999191, 0.2357126 , 0.23165413,
307 0.22779816, 0.22412841, 0.22063036, 0.21729103,
308 0.21409878, 0.21104314, 0.20811462, 0.20530466,
309 0.20260547, 0.20000997;
310
311 CALL_SUBTEST(res = bessel_k1e(x);
312 verify_component_wise(res, expected););
313 }
314
315 // Test Bessel function y0. Reference results obtained with SciPy.
316 {
317 ArrayType x(42);
318 ArrayType expected(42);
319 ArrayType res(42);
320
321 x << 0.25, 0.5, 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12.,
322 13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25.,
323 26., 27., 28., 29., 30., 31., 32., 33., 34., 35., 36., 37., 38.,
324 39., 40.;
325
326 expected << -0.93157302, -0.44451873, 0.08825696, 0.51037567, 0.37685001,
327 -0.01694074, -0.30851763, -0.28819468, -0.02594974, 0.22352149,
328 0.2499367 , 0.05567117, -0.16884732, -0.22523731, -0.07820786,
329 0.12719257, 0.2054643 , 0.095811 , -0.0926372 , -0.18755216,
330 -0.10951969, 0.0626406 , 0.17020176, 0.1198876 , -0.03598179,
331 -0.15283403, -0.12724943, 0.01204463, 0.13521498, 0.13183647,
332 0.00948116, -0.11729573, -0.13383266, -0.02874248, 0.09913483,
333 0.13340405, 0.04579799, -0.08085609, -0.13071488, -0.06066076,
334 0.06262353, 0.12593642;
335
336 CALL_SUBTEST(res = bessel_y0(x);
337 verify_component_wise(res, expected););
338 }
339
340 // Test Bessel function y1. Reference results obtained with SciPy.
341 {
342 ArrayType x(42);
343 ArrayType expected(42);
344 ArrayType res(42);
345
346 x << 0.25, 0.5, 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12.,
347 13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25.,
348 26., 27., 28., 29., 30., 31., 32., 33., 34., 35., 36., 37., 38.,
349 39., 40.;
350
351 expected << -2.70410523, -1.47147239, -0.78121282, -0.10703243,
352 0.32467442, 0.39792571, 0.14786314, -0.17501034, -0.30266724,
353 -0.15806046, 0.10431458, 0.24901542, 0.16370554, -0.05709922,
354 -0.21008141, -0.16664484, 0.02107363, 0.17797517, 0.16720504,
355 0.00815513, -0.14956011, -0.16551161, -0.03253926, 0.12340586,
356 0.1616692 , 0.05305978, -0.09882996, -0.15579655, -0.07025124,
357 0.07552213, 0.14803412, 0.08442557, -0.05337283, -0.13854483,
358 -0.09578012, 0.03238588, 0.12751273, 0.10445477, -0.01262946,
359 -0.11514066, -0.11056411, -0.00579351;
360
361 CALL_SUBTEST(res = bessel_y1(x);
362 verify_component_wise(res, expected););
363 }
364 }
365
EIGEN_DECLARE_TEST(bessel_functions)366 EIGEN_DECLARE_TEST(bessel_functions)
367 {
368 CALL_SUBTEST_1(array_bessel_functions<ArrayXf>());
369 CALL_SUBTEST_2(array_bessel_functions<ArrayXd>());
370 }
371