bessel_functions.cpp (16409B)
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> 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 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 366 EIGEN_DECLARE_TEST(bessel_functions) 367 { 368 CALL_SUBTEST_1(array_bessel_functions<ArrayXf>()); 369 CALL_SUBTEST_2(array_bessel_functions<ArrayXd>()); 370 }