1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2015 Eugene Brevdo <[email protected]> 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 #ifndef EIGEN_BESSEL_FUNCTIONS_H 11 #define EIGEN_BESSEL_FUNCTIONS_H 12 13 namespace Eigen { 14 namespace internal { 15 16 // Parts of this code are based on the Cephes Math Library. 17 // 18 // Cephes Math Library Release 2.8: June, 2000 19 // Copyright 1984, 1987, 1992, 2000 by Stephen L. Moshier 20 // 21 // Permission has been kindly provided by the original author 22 // to incorporate the Cephes software into the Eigen codebase: 23 // 24 // From: Stephen Moshier 25 // To: Eugene Brevdo 26 // Subject: Re: Permission to wrap several cephes functions in Eigen 27 // 28 // Hello Eugene, 29 // 30 // Thank you for writing. 31 // 32 // If your licensing is similar to BSD, the formal way that has been 33 // handled is simply to add a statement to the effect that you are incorporating 34 // the Cephes software by permission of the author. 35 // 36 // Good luck with your project, 37 // Steve 38 39 40 /**************************************************************************** 41 * Implementation of Bessel function, based on Cephes * 42 ****************************************************************************/ 43 44 template <typename Scalar> 45 struct bessel_i0e_retval { 46 typedef Scalar type; 47 }; 48 49 template <typename T, typename ScalarType = typename unpacket_traits<T>::type> 50 struct generic_i0e { 51 EIGEN_DEVICE_FUNC rungeneric_i0e52 static EIGEN_STRONG_INLINE T run(const T&) { 53 EIGEN_STATIC_ASSERT((internal::is_same<T, T>::value == false), 54 THIS_TYPE_IS_NOT_SUPPORTED); 55 return ScalarType(0); 56 } 57 }; 58 59 template <typename T> 60 struct generic_i0e<T, float> { 61 EIGEN_DEVICE_FUNC 62 static EIGEN_STRONG_INLINE T run(const T& x) { 63 /* i0ef.c 64 * 65 * Modified Bessel function of order zero, 66 * exponentially scaled 67 * 68 * 69 * 70 * SYNOPSIS: 71 * 72 * float x, y, i0ef(); 73 * 74 * y = i0ef( x ); 75 * 76 * 77 * 78 * DESCRIPTION: 79 * 80 * Returns exponentially scaled modified Bessel function 81 * of order zero of the argument. 82 * 83 * The function is defined as i0e(x) = exp(-|x|) j0( ix ). 84 * 85 * 86 * 87 * ACCURACY: 88 * 89 * Relative error: 90 * arithmetic domain # trials peak rms 91 * IEEE 0,30 100000 3.7e-7 7.0e-8 92 * See i0f(). 93 * 94 */ 95 96 const float A[] = {-1.30002500998624804212E-8f, 6.04699502254191894932E-8f, 97 -2.67079385394061173391E-7f, 1.11738753912010371815E-6f, 98 -4.41673835845875056359E-6f, 1.64484480707288970893E-5f, 99 -5.75419501008210370398E-5f, 1.88502885095841655729E-4f, 100 -5.76375574538582365885E-4f, 1.63947561694133579842E-3f, 101 -4.32430999505057594430E-3f, 1.05464603945949983183E-2f, 102 -2.37374148058994688156E-2f, 4.93052842396707084878E-2f, 103 -9.49010970480476444210E-2f, 1.71620901522208775349E-1f, 104 -3.04682672343198398683E-1f, 6.76795274409476084995E-1f}; 105 106 const float B[] = {3.39623202570838634515E-9f, 2.26666899049817806459E-8f, 107 2.04891858946906374183E-7f, 2.89137052083475648297E-6f, 108 6.88975834691682398426E-5f, 3.36911647825569408990E-3f, 109 8.04490411014108831608E-1f}; 110 T y = pabs(x); 111 T y_le_eight = internal::pchebevl<T, 18>::run( 112 pmadd(pset1<T>(0.5f), y, pset1<T>(-2.0f)), A); 113 T y_gt_eight = pmul( 114 internal::pchebevl<T, 7>::run( 115 psub(pdiv(pset1<T>(32.0f), y), pset1<T>(2.0f)), B), 116 prsqrt(y)); 117 // TODO: Perhaps instead check whether all packet elements are in 118 // [-8, 8] and evaluate a branch based off of that. It's possible 119 // in practice most elements are in this region. 120 return pselect(pcmp_le(y, pset1<T>(8.0f)), y_le_eight, y_gt_eight); 121 } 122 }; 123 124 template <typename T> 125 struct generic_i0e<T, double> { 126 EIGEN_DEVICE_FUNC 127 static EIGEN_STRONG_INLINE T run(const T& x) { 128 /* i0e.c 129 * 130 * Modified Bessel function of order zero, 131 * exponentially scaled 132 * 133 * 134 * 135 * SYNOPSIS: 136 * 137 * double x, y, i0e(); 138 * 139 * y = i0e( x ); 140 * 141 * 142 * 143 * DESCRIPTION: 144 * 145 * Returns exponentially scaled modified Bessel function 146 * of order zero of the argument. 147 * 148 * The function is defined as i0e(x) = exp(-|x|) j0( ix ). 149 * 150 * 151 * 152 * ACCURACY: 153 * 154 * Relative error: 155 * arithmetic domain # trials peak rms 156 * IEEE 0,30 30000 5.4e-16 1.2e-16 157 * See i0(). 158 * 159 */ 160 161 const double A[] = {-4.41534164647933937950E-18, 3.33079451882223809783E-17, 162 -2.43127984654795469359E-16, 1.71539128555513303061E-15, 163 -1.16853328779934516808E-14, 7.67618549860493561688E-14, 164 -4.85644678311192946090E-13, 2.95505266312963983461E-12, 165 -1.72682629144155570723E-11, 9.67580903537323691224E-11, 166 -5.18979560163526290666E-10, 2.65982372468238665035E-9, 167 -1.30002500998624804212E-8, 6.04699502254191894932E-8, 168 -2.67079385394061173391E-7, 1.11738753912010371815E-6, 169 -4.41673835845875056359E-6, 1.64484480707288970893E-5, 170 -5.75419501008210370398E-5, 1.88502885095841655729E-4, 171 -5.76375574538582365885E-4, 1.63947561694133579842E-3, 172 -4.32430999505057594430E-3, 1.05464603945949983183E-2, 173 -2.37374148058994688156E-2, 4.93052842396707084878E-2, 174 -9.49010970480476444210E-2, 1.71620901522208775349E-1, 175 -3.04682672343198398683E-1, 6.76795274409476084995E-1}; 176 const double B[] = { 177 -7.23318048787475395456E-18, -4.83050448594418207126E-18, 178 4.46562142029675999901E-17, 3.46122286769746109310E-17, 179 -2.82762398051658348494E-16, -3.42548561967721913462E-16, 180 1.77256013305652638360E-15, 3.81168066935262242075E-15, 181 -9.55484669882830764870E-15, -4.15056934728722208663E-14, 182 1.54008621752140982691E-14, 3.85277838274214270114E-13, 183 7.18012445138366623367E-13, -1.79417853150680611778E-12, 184 -1.32158118404477131188E-11, -3.14991652796324136454E-11, 185 1.18891471078464383424E-11, 4.94060238822496958910E-10, 186 3.39623202570838634515E-9, 2.26666899049817806459E-8, 187 2.04891858946906374183E-7, 2.89137052083475648297E-6, 188 6.88975834691682398426E-5, 3.36911647825569408990E-3, 189 8.04490411014108831608E-1}; 190 T y = pabs(x); 191 T y_le_eight = internal::pchebevl<T, 30>::run( 192 pmadd(pset1<T>(0.5), y, pset1<T>(-2.0)), A); 193 T y_gt_eight = pmul( 194 internal::pchebevl<T, 25>::run( 195 psub(pdiv(pset1<T>(32.0), y), pset1<T>(2.0)), B), 196 prsqrt(y)); 197 // TODO: Perhaps instead check whether all packet elements are in 198 // [-8, 8] and evaluate a branch based off of that. It's possible 199 // in practice most elements are in this region. 200 return pselect(pcmp_le(y, pset1<T>(8.0)), y_le_eight, y_gt_eight); 201 } 202 }; 203 204 template <typename T> 205 struct bessel_i0e_impl { 206 EIGEN_DEVICE_FUNC 207 static EIGEN_STRONG_INLINE T run(const T x) { 208 return generic_i0e<T>::run(x); 209 } 210 }; 211 212 template <typename Scalar> 213 struct bessel_i0_retval { 214 typedef Scalar type; 215 }; 216 217 template <typename T, typename ScalarType = typename unpacket_traits<T>::type> 218 struct generic_i0 { 219 EIGEN_DEVICE_FUNC 220 static EIGEN_STRONG_INLINE T run(const T& x) { 221 return pmul( 222 pexp(pabs(x)), 223 generic_i0e<T, ScalarType>::run(x)); 224 } 225 }; 226 227 template <typename T> 228 struct bessel_i0_impl { 229 EIGEN_DEVICE_FUNC 230 static EIGEN_STRONG_INLINE T run(const T x) { 231 return generic_i0<T>::run(x); 232 } 233 }; 234 235 template <typename Scalar> 236 struct bessel_i1e_retval { 237 typedef Scalar type; 238 }; 239 240 template <typename T, typename ScalarType = typename unpacket_traits<T>::type > 241 struct generic_i1e { 242 EIGEN_DEVICE_FUNC 243 static EIGEN_STRONG_INLINE T run(const T&) { 244 EIGEN_STATIC_ASSERT((internal::is_same<T, T>::value == false), 245 THIS_TYPE_IS_NOT_SUPPORTED); 246 return ScalarType(0); 247 } 248 }; 249 250 template <typename T> 251 struct generic_i1e<T, float> { 252 EIGEN_DEVICE_FUNC 253 static EIGEN_STRONG_INLINE T run(const T& x) { 254 /* i1ef.c 255 * 256 * Modified Bessel function of order one, 257 * exponentially scaled 258 * 259 * 260 * 261 * SYNOPSIS: 262 * 263 * float x, y, i1ef(); 264 * 265 * y = i1ef( x ); 266 * 267 * 268 * 269 * DESCRIPTION: 270 * 271 * Returns exponentially scaled modified Bessel function 272 * of order one of the argument. 273 * 274 * The function is defined as i1(x) = -i exp(-|x|) j1( ix ). 275 * 276 * 277 * 278 * ACCURACY: 279 * 280 * Relative error: 281 * arithmetic domain # trials peak rms 282 * IEEE 0, 30 30000 1.5e-6 1.5e-7 283 * See i1(). 284 * 285 */ 286 const float A[] = {9.38153738649577178388E-9f, -4.44505912879632808065E-8f, 287 2.00329475355213526229E-7f, -8.56872026469545474066E-7f, 288 3.47025130813767847674E-6f, -1.32731636560394358279E-5f, 289 4.78156510755005422638E-5f, -1.61760815825896745588E-4f, 290 5.12285956168575772895E-4f, -1.51357245063125314899E-3f, 291 4.15642294431288815669E-3f, -1.05640848946261981558E-2f, 292 2.47264490306265168283E-2f, -5.29459812080949914269E-2f, 293 1.02643658689847095384E-1f, -1.76416518357834055153E-1f, 294 2.52587186443633654823E-1f}; 295 296 const float B[] = {-3.83538038596423702205E-9f, -2.63146884688951950684E-8f, 297 -2.51223623787020892529E-7f, -3.88256480887769039346E-6f, 298 -1.10588938762623716291E-4f, -9.76109749136146840777E-3f, 299 7.78576235018280120474E-1f}; 300 301 302 T y = pabs(x); 303 T y_le_eight = pmul(y, internal::pchebevl<T, 17>::run( 304 pmadd(pset1<T>(0.5f), y, pset1<T>(-2.0f)), A)); 305 T y_gt_eight = pmul( 306 internal::pchebevl<T, 7>::run( 307 psub(pdiv(pset1<T>(32.0f), y), 308 pset1<T>(2.0f)), B), 309 prsqrt(y)); 310 // TODO: Perhaps instead check whether all packet elements are in 311 // [-8, 8] and evaluate a branch based off of that. It's possible 312 // in practice most elements are in this region. 313 y = pselect(pcmp_le(y, pset1<T>(8.0f)), y_le_eight, y_gt_eight); 314 return pselect(pcmp_lt(x, pset1<T>(0.0f)), pnegate(y), y); 315 } 316 }; 317 318 template <typename T> 319 struct generic_i1e<T, double> { 320 EIGEN_DEVICE_FUNC 321 static EIGEN_STRONG_INLINE T run(const T& x) { 322 /* i1e.c 323 * 324 * Modified Bessel function of order one, 325 * exponentially scaled 326 * 327 * 328 * 329 * SYNOPSIS: 330 * 331 * double x, y, i1e(); 332 * 333 * y = i1e( x ); 334 * 335 * 336 * 337 * DESCRIPTION: 338 * 339 * Returns exponentially scaled modified Bessel function 340 * of order one of the argument. 341 * 342 * The function is defined as i1(x) = -i exp(-|x|) j1( ix ). 343 * 344 * 345 * 346 * ACCURACY: 347 * 348 * Relative error: 349 * arithmetic domain # trials peak rms 350 * IEEE 0, 30 30000 2.0e-15 2.0e-16 351 * See i1(). 352 * 353 */ 354 const double A[] = {2.77791411276104639959E-18, -2.11142121435816608115E-17, 355 1.55363195773620046921E-16, -1.10559694773538630805E-15, 356 7.60068429473540693410E-15, -5.04218550472791168711E-14, 357 3.22379336594557470981E-13, -1.98397439776494371520E-12, 358 1.17361862988909016308E-11, -6.66348972350202774223E-11, 359 3.62559028155211703701E-10, -1.88724975172282928790E-9, 360 9.38153738649577178388E-9, -4.44505912879632808065E-8, 361 2.00329475355213526229E-7, -8.56872026469545474066E-7, 362 3.47025130813767847674E-6, -1.32731636560394358279E-5, 363 4.78156510755005422638E-5, -1.61760815825896745588E-4, 364 5.12285956168575772895E-4, -1.51357245063125314899E-3, 365 4.15642294431288815669E-3, -1.05640848946261981558E-2, 366 2.47264490306265168283E-2, -5.29459812080949914269E-2, 367 1.02643658689847095384E-1, -1.76416518357834055153E-1, 368 2.52587186443633654823E-1}; 369 const double B[] = { 370 7.51729631084210481353E-18, 4.41434832307170791151E-18, 371 -4.65030536848935832153E-17, -3.20952592199342395980E-17, 372 2.96262899764595013876E-16, 3.30820231092092828324E-16, 373 -1.88035477551078244854E-15, -3.81440307243700780478E-15, 374 1.04202769841288027642E-14, 4.27244001671195135429E-14, 375 -2.10154184277266431302E-14, -4.08355111109219731823E-13, 376 -7.19855177624590851209E-13, 2.03562854414708950722E-12, 377 1.41258074366137813316E-11, 3.25260358301548823856E-11, 378 -1.89749581235054123450E-11, -5.58974346219658380687E-10, 379 -3.83538038596423702205E-9, -2.63146884688951950684E-8, 380 -2.51223623787020892529E-7, -3.88256480887769039346E-6, 381 -1.10588938762623716291E-4, -9.76109749136146840777E-3, 382 7.78576235018280120474E-1}; 383 T y = pabs(x); 384 T y_le_eight = pmul(y, internal::pchebevl<T, 29>::run( 385 pmadd(pset1<T>(0.5), y, pset1<T>(-2.0)), A)); 386 T y_gt_eight = pmul( 387 internal::pchebevl<T, 25>::run( 388 psub(pdiv(pset1<T>(32.0), y), 389 pset1<T>(2.0)), B), 390 prsqrt(y)); 391 // TODO: Perhaps instead check whether all packet elements are in 392 // [-8, 8] and evaluate a branch based off of that. It's possible 393 // in practice most elements are in this region. 394 y = pselect(pcmp_le(y, pset1<T>(8.0)), y_le_eight, y_gt_eight); 395 return pselect(pcmp_lt(x, pset1<T>(0.0)), pnegate(y), y); 396 } 397 }; 398 399 template <typename T> 400 struct bessel_i1e_impl { 401 EIGEN_DEVICE_FUNC 402 static EIGEN_STRONG_INLINE T run(const T x) { 403 return generic_i1e<T>::run(x); 404 } 405 }; 406 407 template <typename T> 408 struct bessel_i1_retval { 409 typedef T type; 410 }; 411 412 template <typename T, typename ScalarType = typename unpacket_traits<T>::type> 413 struct generic_i1 { 414 EIGEN_DEVICE_FUNC 415 static EIGEN_STRONG_INLINE T run(const T& x) { 416 return pmul( 417 pexp(pabs(x)), 418 generic_i1e<T, ScalarType>::run(x)); 419 } 420 }; 421 422 template <typename T> 423 struct bessel_i1_impl { 424 EIGEN_DEVICE_FUNC 425 static EIGEN_STRONG_INLINE T run(const T x) { 426 return generic_i1<T>::run(x); 427 } 428 }; 429 430 template <typename T> 431 struct bessel_k0e_retval { 432 typedef T type; 433 }; 434 435 template <typename T, typename ScalarType = typename unpacket_traits<T>::type> 436 struct generic_k0e { 437 EIGEN_DEVICE_FUNC 438 static EIGEN_STRONG_INLINE T run(const T&) { 439 EIGEN_STATIC_ASSERT((internal::is_same<T, T>::value == false), 440 THIS_TYPE_IS_NOT_SUPPORTED); 441 return ScalarType(0); 442 } 443 }; 444 445 template <typename T> 446 struct generic_k0e<T, float> { 447 EIGEN_DEVICE_FUNC 448 static EIGEN_STRONG_INLINE T run(const T& x) { 449 /* k0ef.c 450 * Modified Bessel function, third kind, order zero, 451 * exponentially scaled 452 * 453 * 454 * 455 * SYNOPSIS: 456 * 457 * float x, y, k0ef(); 458 * 459 * y = k0ef( x ); 460 * 461 * 462 * 463 * DESCRIPTION: 464 * 465 * Returns exponentially scaled modified Bessel function 466 * of the third kind of order zero of the argument. 467 * 468 * 469 * 470 * ACCURACY: 471 * 472 * Relative error: 473 * arithmetic domain # trials peak rms 474 * IEEE 0, 30 30000 8.1e-7 7.8e-8 475 * See k0(). 476 * 477 */ 478 479 const float A[] = {1.90451637722020886025E-9f, 2.53479107902614945675E-7f, 480 2.28621210311945178607E-5f, 1.26461541144692592338E-3f, 481 3.59799365153615016266E-2f, 3.44289899924628486886E-1f, 482 -5.35327393233902768720E-1f}; 483 484 const float B[] = {-1.69753450938905987466E-9f, 8.57403401741422608519E-9f, 485 -4.66048989768794782956E-8f, 2.76681363944501510342E-7f, 486 -1.83175552271911948767E-6f, 1.39498137188764993662E-5f, 487 -1.28495495816278026384E-4f, 1.56988388573005337491E-3f, 488 -3.14481013119645005427E-2f, 2.44030308206595545468E0f}; 489 const T MAXNUM = pset1<T>(NumTraits<float>::infinity()); 490 const T two = pset1<T>(2.0); 491 T x_le_two = internal::pchebevl<T, 7>::run( 492 pmadd(x, x, pset1<T>(-2.0)), A); 493 x_le_two = pmadd( 494 generic_i0<T, float>::run(x), pnegate( 495 plog(pmul(pset1<T>(0.5), x))), x_le_two); 496 x_le_two = pmul(pexp(x), x_le_two); 497 T x_gt_two = pmul( 498 internal::pchebevl<T, 10>::run( 499 psub(pdiv(pset1<T>(8.0), x), two), B), 500 prsqrt(x)); 501 return pselect( 502 pcmp_le(x, pset1<T>(0.0)), 503 MAXNUM, 504 pselect(pcmp_le(x, two), x_le_two, x_gt_two)); 505 } 506 }; 507 508 template <typename T> 509 struct generic_k0e<T, double> { 510 EIGEN_DEVICE_FUNC 511 static EIGEN_STRONG_INLINE T run(const T& x) { 512 /* k0e.c 513 * Modified Bessel function, third kind, order zero, 514 * exponentially scaled 515 * 516 * 517 * 518 * SYNOPSIS: 519 * 520 * double x, y, k0e(); 521 * 522 * y = k0e( x ); 523 * 524 * 525 * 526 * DESCRIPTION: 527 * 528 * Returns exponentially scaled modified Bessel function 529 * of the third kind of order zero of the argument. 530 * 531 * 532 * 533 * ACCURACY: 534 * 535 * Relative error: 536 * arithmetic domain # trials peak rms 537 * IEEE 0, 30 30000 1.4e-15 1.4e-16 538 * See k0(). 539 * 540 */ 541 542 const double A[] = { 543 1.37446543561352307156E-16, 544 4.25981614279661018399E-14, 545 1.03496952576338420167E-11, 546 1.90451637722020886025E-9, 547 2.53479107902614945675E-7, 548 2.28621210311945178607E-5, 549 1.26461541144692592338E-3, 550 3.59799365153615016266E-2, 551 3.44289899924628486886E-1, 552 -5.35327393233902768720E-1}; 553 const double B[] = { 554 5.30043377268626276149E-18, -1.64758043015242134646E-17, 555 5.21039150503902756861E-17, -1.67823109680541210385E-16, 556 5.51205597852431940784E-16, -1.84859337734377901440E-15, 557 6.34007647740507060557E-15, -2.22751332699166985548E-14, 558 8.03289077536357521100E-14, -2.98009692317273043925E-13, 559 1.14034058820847496303E-12, -4.51459788337394416547E-12, 560 1.85594911495471785253E-11, -7.95748924447710747776E-11, 561 3.57739728140030116597E-10, -1.69753450938905987466E-9, 562 8.57403401741422608519E-9, -4.66048989768794782956E-8, 563 2.76681363944501510342E-7, -1.83175552271911948767E-6, 564 1.39498137188764993662E-5, -1.28495495816278026384E-4, 565 1.56988388573005337491E-3, -3.14481013119645005427E-2, 566 2.44030308206595545468E0 567 }; 568 const T MAXNUM = pset1<T>(NumTraits<double>::infinity()); 569 const T two = pset1<T>(2.0); 570 T x_le_two = internal::pchebevl<T, 10>::run( 571 pmadd(x, x, pset1<T>(-2.0)), A); 572 x_le_two = pmadd( 573 generic_i0<T, double>::run(x), pmul( 574 pset1<T>(-1.0), plog(pmul(pset1<T>(0.5), x))), x_le_two); 575 x_le_two = pmul(pexp(x), x_le_two); 576 x_le_two = pselect(pcmp_le(x, pset1<T>(0.0)), MAXNUM, x_le_two); 577 T x_gt_two = pmul( 578 internal::pchebevl<T, 25>::run( 579 psub(pdiv(pset1<T>(8.0), x), two), B), 580 prsqrt(x)); 581 return pselect(pcmp_le(x, two), x_le_two, x_gt_two); 582 } 583 }; 584 585 template <typename T> 586 struct bessel_k0e_impl { 587 EIGEN_DEVICE_FUNC 588 static EIGEN_STRONG_INLINE T run(const T x) { 589 return generic_k0e<T>::run(x); 590 } 591 }; 592 593 template <typename T> 594 struct bessel_k0_retval { 595 typedef T type; 596 }; 597 598 template <typename T, typename ScalarType = typename unpacket_traits<T>::type> 599 struct generic_k0 { 600 EIGEN_DEVICE_FUNC 601 static EIGEN_STRONG_INLINE T run(const T&) { 602 EIGEN_STATIC_ASSERT((internal::is_same<T, T>::value == false), 603 THIS_TYPE_IS_NOT_SUPPORTED); 604 return ScalarType(0); 605 } 606 }; 607 608 template <typename T> 609 struct generic_k0<T, float> { 610 EIGEN_DEVICE_FUNC 611 static EIGEN_STRONG_INLINE T run(const T& x) { 612 /* k0f.c 613 * Modified Bessel function, third kind, order zero 614 * 615 * 616 * 617 * SYNOPSIS: 618 * 619 * float x, y, k0f(); 620 * 621 * y = k0f( x ); 622 * 623 * 624 * 625 * DESCRIPTION: 626 * 627 * Returns modified Bessel function of the third kind 628 * of order zero of the argument. 629 * 630 * The range is partitioned into the two intervals [0,8] and 631 * (8, infinity). Chebyshev polynomial expansions are employed 632 * in each interval. 633 * 634 * 635 * 636 * ACCURACY: 637 * 638 * Tested at 2000 random points between 0 and 8. Peak absolute 639 * error (relative when K0 > 1) was 1.46e-14; rms, 4.26e-15. 640 * Relative error: 641 * arithmetic domain # trials peak rms 642 * IEEE 0, 30 30000 7.8e-7 8.5e-8 643 * 644 * ERROR MESSAGES: 645 * 646 * message condition value returned 647 * K0 domain x <= 0 MAXNUM 648 * 649 */ 650 651 const float A[] = {1.90451637722020886025E-9f, 2.53479107902614945675E-7f, 652 2.28621210311945178607E-5f, 1.26461541144692592338E-3f, 653 3.59799365153615016266E-2f, 3.44289899924628486886E-1f, 654 -5.35327393233902768720E-1f}; 655 656 const float B[] = {-1.69753450938905987466E-9f, 8.57403401741422608519E-9f, 657 -4.66048989768794782956E-8f, 2.76681363944501510342E-7f, 658 -1.83175552271911948767E-6f, 1.39498137188764993662E-5f, 659 -1.28495495816278026384E-4f, 1.56988388573005337491E-3f, 660 -3.14481013119645005427E-2f, 2.44030308206595545468E0f}; 661 const T MAXNUM = pset1<T>(NumTraits<float>::infinity()); 662 const T two = pset1<T>(2.0); 663 T x_le_two = internal::pchebevl<T, 7>::run( 664 pmadd(x, x, pset1<T>(-2.0)), A); 665 x_le_two = pmadd( 666 generic_i0<T, float>::run(x), pnegate( 667 plog(pmul(pset1<T>(0.5), x))), x_le_two); 668 x_le_two = pselect(pcmp_le(x, pset1<T>(0.0)), MAXNUM, x_le_two); 669 T x_gt_two = pmul( 670 pmul( 671 pexp(pnegate(x)), 672 internal::pchebevl<T, 10>::run( 673 psub(pdiv(pset1<T>(8.0), x), two), B)), 674 prsqrt(x)); 675 return pselect(pcmp_le(x, two), x_le_two, x_gt_two); 676 } 677 }; 678 679 template <typename T> 680 struct generic_k0<T, double> { 681 EIGEN_DEVICE_FUNC 682 static EIGEN_STRONG_INLINE T run(const T& x) { 683 /* 684 * 685 * Modified Bessel function, third kind, order zero, 686 * exponentially scaled 687 * 688 * 689 * 690 * SYNOPSIS: 691 * 692 * double x, y, k0(); 693 * 694 * y = k0( x ); 695 * 696 * 697 * 698 * DESCRIPTION: 699 * 700 * Returns exponentially scaled modified Bessel function 701 * of the third kind of order zero of the argument. 702 * 703 * 704 * 705 * ACCURACY: 706 * 707 * Relative error: 708 * arithmetic domain # trials peak rms 709 * IEEE 0, 30 30000 1.4e-15 1.4e-16 710 * See k0(). 711 * 712 */ 713 const double A[] = { 714 1.37446543561352307156E-16, 715 4.25981614279661018399E-14, 716 1.03496952576338420167E-11, 717 1.90451637722020886025E-9, 718 2.53479107902614945675E-7, 719 2.28621210311945178607E-5, 720 1.26461541144692592338E-3, 721 3.59799365153615016266E-2, 722 3.44289899924628486886E-1, 723 -5.35327393233902768720E-1}; 724 const double B[] = { 725 5.30043377268626276149E-18, -1.64758043015242134646E-17, 726 5.21039150503902756861E-17, -1.67823109680541210385E-16, 727 5.51205597852431940784E-16, -1.84859337734377901440E-15, 728 6.34007647740507060557E-15, -2.22751332699166985548E-14, 729 8.03289077536357521100E-14, -2.98009692317273043925E-13, 730 1.14034058820847496303E-12, -4.51459788337394416547E-12, 731 1.85594911495471785253E-11, -7.95748924447710747776E-11, 732 3.57739728140030116597E-10, -1.69753450938905987466E-9, 733 8.57403401741422608519E-9, -4.66048989768794782956E-8, 734 2.76681363944501510342E-7, -1.83175552271911948767E-6, 735 1.39498137188764993662E-5, -1.28495495816278026384E-4, 736 1.56988388573005337491E-3, -3.14481013119645005427E-2, 737 2.44030308206595545468E0 738 }; 739 const T MAXNUM = pset1<T>(NumTraits<double>::infinity()); 740 const T two = pset1<T>(2.0); 741 T x_le_two = internal::pchebevl<T, 10>::run( 742 pmadd(x, x, pset1<T>(-2.0)), A); 743 x_le_two = pmadd( 744 generic_i0<T, double>::run(x), pnegate( 745 plog(pmul(pset1<T>(0.5), x))), x_le_two); 746 x_le_two = pselect(pcmp_le(x, pset1<T>(0.0)), MAXNUM, x_le_two); 747 T x_gt_two = pmul( 748 pmul( 749 pexp(-x), 750 internal::pchebevl<T, 25>::run( 751 psub(pdiv(pset1<T>(8.0), x), two), B)), 752 prsqrt(x)); 753 return pselect(pcmp_le(x, two), x_le_two, x_gt_two); 754 } 755 }; 756 757 template <typename T> 758 struct bessel_k0_impl { 759 EIGEN_DEVICE_FUNC 760 static EIGEN_STRONG_INLINE T run(const T x) { 761 return generic_k0<T>::run(x); 762 } 763 }; 764 765 template <typename T> 766 struct bessel_k1e_retval { 767 typedef T type; 768 }; 769 770 template <typename T, typename ScalarType = typename unpacket_traits<T>::type> 771 struct generic_k1e { 772 EIGEN_DEVICE_FUNC 773 static EIGEN_STRONG_INLINE T run(const T&) { 774 EIGEN_STATIC_ASSERT((internal::is_same<T, T>::value == false), 775 THIS_TYPE_IS_NOT_SUPPORTED); 776 return ScalarType(0); 777 } 778 }; 779 780 template <typename T> 781 struct generic_k1e<T, float> { 782 EIGEN_DEVICE_FUNC 783 static EIGEN_STRONG_INLINE T run(const T& x) { 784 /* k1ef.c 785 * 786 * Modified Bessel function, third kind, order one, 787 * exponentially scaled 788 * 789 * 790 * 791 * SYNOPSIS: 792 * 793 * float x, y, k1ef(); 794 * 795 * y = k1ef( x ); 796 * 797 * 798 * 799 * DESCRIPTION: 800 * 801 * Returns exponentially scaled modified Bessel function 802 * of the third kind of order one of the argument: 803 * 804 * k1e(x) = exp(x) * k1(x). 805 * 806 * 807 * 808 * ACCURACY: 809 * 810 * Relative error: 811 * arithmetic domain # trials peak rms 812 * IEEE 0, 30 30000 4.9e-7 6.7e-8 813 * See k1(). 814 * 815 */ 816 817 const float A[] = {-2.21338763073472585583E-8f, -2.43340614156596823496E-6f, 818 -1.73028895751305206302E-4f, -6.97572385963986435018E-3f, 819 -1.22611180822657148235E-1f, -3.53155960776544875667E-1f, 820 1.52530022733894777053E0f}; 821 const float B[] = {2.01504975519703286596E-9f, -1.03457624656780970260E-8f, 822 5.74108412545004946722E-8f, -3.50196060308781257119E-7f, 823 2.40648494783721712015E-6f, -1.93619797416608296024E-5f, 824 1.95215518471351631108E-4f, -2.85781685962277938680E-3f, 825 1.03923736576817238437E-1f, 2.72062619048444266945E0f}; 826 const T MAXNUM = pset1<T>(NumTraits<float>::infinity()); 827 const T two = pset1<T>(2.0); 828 T x_le_two = pdiv(internal::pchebevl<T, 7>::run( 829 pmadd(x, x, pset1<T>(-2.0)), A), x); 830 x_le_two = pmadd( 831 generic_i1<T, float>::run(x), plog(pmul(pset1<T>(0.5), x)), x_le_two); 832 x_le_two = pmul(x_le_two, pexp(x)); 833 x_le_two = pselect(pcmp_le(x, pset1<T>(0.0)), MAXNUM, x_le_two); 834 T x_gt_two = pmul( 835 internal::pchebevl<T, 10>::run( 836 psub(pdiv(pset1<T>(8.0), x), two), B), 837 prsqrt(x)); 838 return pselect(pcmp_le(x, two), x_le_two, x_gt_two); 839 } 840 }; 841 842 template <typename T> 843 struct generic_k1e<T, double> { 844 EIGEN_DEVICE_FUNC 845 static EIGEN_STRONG_INLINE T run(const T& x) { 846 /* k1e.c 847 * 848 * Modified Bessel function, third kind, order one, 849 * exponentially scaled 850 * 851 * 852 * 853 * SYNOPSIS: 854 * 855 * double x, y, k1e(); 856 * 857 * y = k1e( x ); 858 * 859 * 860 * 861 * DESCRIPTION: 862 * 863 * Returns exponentially scaled modified Bessel function 864 * of the third kind of order one of the argument: 865 * 866 * k1e(x) = exp(x) * k1(x). 867 * 868 * 869 * 870 * ACCURACY: 871 * 872 * Relative error: 873 * arithmetic domain # trials peak rms 874 * IEEE 0, 30 30000 7.8e-16 1.2e-16 875 * See k1(). 876 * 877 */ 878 const double A[] = {-7.02386347938628759343E-18, -2.42744985051936593393E-15, 879 -6.66690169419932900609E-13, -1.41148839263352776110E-10, 880 -2.21338763073472585583E-8, -2.43340614156596823496E-6, 881 -1.73028895751305206302E-4, -6.97572385963986435018E-3, 882 -1.22611180822657148235E-1, -3.53155960776544875667E-1, 883 1.52530022733894777053E0}; 884 const double B[] = {-5.75674448366501715755E-18, 1.79405087314755922667E-17, 885 -5.68946255844285935196E-17, 1.83809354436663880070E-16, 886 -6.05704724837331885336E-16, 2.03870316562433424052E-15, 887 -7.01983709041831346144E-15, 2.47715442448130437068E-14, 888 -8.97670518232499435011E-14, 3.34841966607842919884E-13, 889 -1.28917396095102890680E-12, 5.13963967348173025100E-12, 890 -2.12996783842756842877E-11, 9.21831518760500529508E-11, 891 -4.19035475934189648750E-10, 2.01504975519703286596E-9, 892 -1.03457624656780970260E-8, 5.74108412545004946722E-8, 893 -3.50196060308781257119E-7, 2.40648494783721712015E-6, 894 -1.93619797416608296024E-5, 1.95215518471351631108E-4, 895 -2.85781685962277938680E-3, 1.03923736576817238437E-1, 896 2.72062619048444266945E0}; 897 const T MAXNUM = pset1<T>(NumTraits<double>::infinity()); 898 const T two = pset1<T>(2.0); 899 T x_le_two = pdiv(internal::pchebevl<T, 11>::run( 900 pmadd(x, x, pset1<T>(-2.0)), A), x); 901 x_le_two = pmadd( 902 generic_i1<T, double>::run(x), plog(pmul(pset1<T>(0.5), x)), x_le_two); 903 x_le_two = pmul(x_le_two, pexp(x)); 904 x_le_two = pselect(pcmp_le(x, pset1<T>(0.0)), MAXNUM, x_le_two); 905 T x_gt_two = pmul( 906 internal::pchebevl<T, 25>::run( 907 psub(pdiv(pset1<T>(8.0), x), two), B), 908 prsqrt(x)); 909 return pselect(pcmp_le(x, two), x_le_two, x_gt_two); 910 } 911 }; 912 913 template <typename T> 914 struct bessel_k1e_impl { 915 EIGEN_DEVICE_FUNC 916 static EIGEN_STRONG_INLINE T run(const T x) { 917 return generic_k1e<T>::run(x); 918 } 919 }; 920 921 template <typename T> 922 struct bessel_k1_retval { 923 typedef T type; 924 }; 925 926 template <typename T, typename ScalarType = typename unpacket_traits<T>::type> 927 struct generic_k1 { 928 EIGEN_DEVICE_FUNC 929 static EIGEN_STRONG_INLINE T run(const T&) { 930 EIGEN_STATIC_ASSERT((internal::is_same<T, T>::value == false), 931 THIS_TYPE_IS_NOT_SUPPORTED); 932 return ScalarType(0); 933 } 934 }; 935 936 template <typename T> 937 struct generic_k1<T, float> { 938 EIGEN_DEVICE_FUNC 939 static EIGEN_STRONG_INLINE T run(const T& x) { 940 /* k1f.c 941 * Modified Bessel function, third kind, order one 942 * 943 * 944 * 945 * SYNOPSIS: 946 * 947 * float x, y, k1f(); 948 * 949 * y = k1f( x ); 950 * 951 * 952 * 953 * DESCRIPTION: 954 * 955 * Computes the modified Bessel function of the third kind 956 * of order one of the argument. 957 * 958 * The range is partitioned into the two intervals [0,2] and 959 * (2, infinity). Chebyshev polynomial expansions are employed 960 * in each interval. 961 * 962 * 963 * 964 * ACCURACY: 965 * 966 * Relative error: 967 * arithmetic domain # trials peak rms 968 * IEEE 0, 30 30000 4.6e-7 7.6e-8 969 * 970 * ERROR MESSAGES: 971 * 972 * message condition value returned 973 * k1 domain x <= 0 MAXNUM 974 * 975 */ 976 977 const float A[] = {-2.21338763073472585583E-8f, -2.43340614156596823496E-6f, 978 -1.73028895751305206302E-4f, -6.97572385963986435018E-3f, 979 -1.22611180822657148235E-1f, -3.53155960776544875667E-1f, 980 1.52530022733894777053E0f}; 981 const float B[] = {2.01504975519703286596E-9f, -1.03457624656780970260E-8f, 982 5.74108412545004946722E-8f, -3.50196060308781257119E-7f, 983 2.40648494783721712015E-6f, -1.93619797416608296024E-5f, 984 1.95215518471351631108E-4f, -2.85781685962277938680E-3f, 985 1.03923736576817238437E-1f, 2.72062619048444266945E0f}; 986 const T MAXNUM = pset1<T>(NumTraits<float>::infinity()); 987 const T two = pset1<T>(2.0); 988 T x_le_two = pdiv(internal::pchebevl<T, 7>::run( 989 pmadd(x, x, pset1<T>(-2.0)), A), x); 990 x_le_two = pmadd( 991 generic_i1<T, float>::run(x), plog(pmul(pset1<T>(0.5), x)), x_le_two); 992 x_le_two = pselect(pcmp_le(x, pset1<T>(0.0)), MAXNUM, x_le_two); 993 T x_gt_two = pmul( 994 pexp(pnegate(x)), 995 pmul( 996 internal::pchebevl<T, 10>::run( 997 psub(pdiv(pset1<T>(8.0), x), two), B), 998 prsqrt(x))); 999 return pselect(pcmp_le(x, two), x_le_two, x_gt_two); 1000 } 1001 }; 1002 1003 template <typename T> 1004 struct generic_k1<T, double> { 1005 EIGEN_DEVICE_FUNC 1006 static EIGEN_STRONG_INLINE T run(const T& x) { 1007 /* k1.c 1008 * Modified Bessel function, third kind, order one 1009 * 1010 * 1011 * 1012 * SYNOPSIS: 1013 * 1014 * float x, y, k1f(); 1015 * 1016 * y = k1f( x ); 1017 * 1018 * 1019 * 1020 * DESCRIPTION: 1021 * 1022 * Computes the modified Bessel function of the third kind 1023 * of order one of the argument. 1024 * 1025 * The range is partitioned into the two intervals [0,2] and 1026 * (2, infinity). Chebyshev polynomial expansions are employed 1027 * in each interval. 1028 * 1029 * 1030 * 1031 * ACCURACY: 1032 * 1033 * Relative error: 1034 * arithmetic domain # trials peak rms 1035 * IEEE 0, 30 30000 4.6e-7 7.6e-8 1036 * 1037 * ERROR MESSAGES: 1038 * 1039 * message condition value returned 1040 * k1 domain x <= 0 MAXNUM 1041 * 1042 */ 1043 const double A[] = {-7.02386347938628759343E-18, -2.42744985051936593393E-15, 1044 -6.66690169419932900609E-13, -1.41148839263352776110E-10, 1045 -2.21338763073472585583E-8, -2.43340614156596823496E-6, 1046 -1.73028895751305206302E-4, -6.97572385963986435018E-3, 1047 -1.22611180822657148235E-1, -3.53155960776544875667E-1, 1048 1.52530022733894777053E0}; 1049 const double B[] = {-5.75674448366501715755E-18, 1.79405087314755922667E-17, 1050 -5.68946255844285935196E-17, 1.83809354436663880070E-16, 1051 -6.05704724837331885336E-16, 2.03870316562433424052E-15, 1052 -7.01983709041831346144E-15, 2.47715442448130437068E-14, 1053 -8.97670518232499435011E-14, 3.34841966607842919884E-13, 1054 -1.28917396095102890680E-12, 5.13963967348173025100E-12, 1055 -2.12996783842756842877E-11, 9.21831518760500529508E-11, 1056 -4.19035475934189648750E-10, 2.01504975519703286596E-9, 1057 -1.03457624656780970260E-8, 5.74108412545004946722E-8, 1058 -3.50196060308781257119E-7, 2.40648494783721712015E-6, 1059 -1.93619797416608296024E-5, 1.95215518471351631108E-4, 1060 -2.85781685962277938680E-3, 1.03923736576817238437E-1, 1061 2.72062619048444266945E0}; 1062 const T MAXNUM = pset1<T>(NumTraits<double>::infinity()); 1063 const T two = pset1<T>(2.0); 1064 T x_le_two = pdiv(internal::pchebevl<T, 11>::run( 1065 pmadd(x, x, pset1<T>(-2.0)), A), x); 1066 x_le_two = pmadd( 1067 generic_i1<T, double>::run(x), plog(pmul(pset1<T>(0.5), x)), x_le_two); 1068 x_le_two = pselect(pcmp_le(x, pset1<T>(0.0)), MAXNUM, x_le_two); 1069 T x_gt_two = pmul( 1070 pexp(-x), 1071 pmul( 1072 internal::pchebevl<T, 25>::run( 1073 psub(pdiv(pset1<T>(8.0), x), two), B), 1074 prsqrt(x))); 1075 return pselect(pcmp_le(x, two), x_le_two, x_gt_two); 1076 } 1077 }; 1078 1079 template <typename T> 1080 struct bessel_k1_impl { 1081 EIGEN_DEVICE_FUNC 1082 static EIGEN_STRONG_INLINE T run(const T x) { 1083 return generic_k1<T>::run(x); 1084 } 1085 }; 1086 1087 template <typename T> 1088 struct bessel_j0_retval { 1089 typedef T type; 1090 }; 1091 1092 template <typename T, typename ScalarType = typename unpacket_traits<T>::type> 1093 struct generic_j0 { 1094 EIGEN_DEVICE_FUNC 1095 static EIGEN_STRONG_INLINE T run(const T&) { 1096 EIGEN_STATIC_ASSERT((internal::is_same<T, T>::value == false), 1097 THIS_TYPE_IS_NOT_SUPPORTED); 1098 return ScalarType(0); 1099 } 1100 }; 1101 1102 template <typename T> 1103 struct generic_j0<T, float> { 1104 EIGEN_DEVICE_FUNC 1105 static EIGEN_STRONG_INLINE T run(const T& x) { 1106 /* j0f.c 1107 * Bessel function of order zero 1108 * 1109 * 1110 * 1111 * SYNOPSIS: 1112 * 1113 * float x, y, j0f(); 1114 * 1115 * y = j0f( x ); 1116 * 1117 * 1118 * 1119 * DESCRIPTION: 1120 * 1121 * Returns Bessel function of order zero of the argument. 1122 * 1123 * The domain is divided into the intervals [0, 2] and 1124 * (2, infinity). In the first interval the following polynomial 1125 * approximation is used: 1126 * 1127 * 1128 * 2 2 2 1129 * (w - r ) (w - r ) (w - r ) P(w) 1130 * 1 2 3 1131 * 1132 * 2 1133 * where w = x and the three r's are zeros of the function. 1134 * 1135 * In the second interval, the modulus and phase are approximated 1136 * by polynomials of the form Modulus(x) = sqrt(1/x) Q(1/x) 1137 * and Phase(x) = x + 1/x R(1/x^2) - pi/4. The function is 1138 * 1139 * j0(x) = Modulus(x) cos( Phase(x) ). 1140 * 1141 * 1142 * 1143 * ACCURACY: 1144 * 1145 * Absolute error: 1146 * arithmetic domain # trials peak rms 1147 * IEEE 0, 2 100000 1.3e-7 3.6e-8 1148 * IEEE 2, 32 100000 1.9e-7 5.4e-8 1149 * 1150 */ 1151 1152 const float JP[] = {-6.068350350393235E-008f, 6.388945720783375E-006f, 1153 -3.969646342510940E-004f, 1.332913422519003E-002f, 1154 -1.729150680240724E-001f}; 1155 const float MO[] = {-6.838999669318810E-002f, 1.864949361379502E-001f, 1156 -2.145007480346739E-001f, 1.197549369473540E-001f, 1157 -3.560281861530129E-003f, -4.969382655296620E-002f, 1158 -3.355424622293709E-006f, 7.978845717621440E-001f}; 1159 const float PH[] = {3.242077816988247E+001f, -3.630592630518434E+001f, 1160 1.756221482109099E+001f, -4.974978466280903E+000f, 1161 1.001973420681837E+000f, -1.939906941791308E-001f, 1162 6.490598792654666E-002f, -1.249992184872738E-001f}; 1163 const T DR1 = pset1<T>(5.78318596294678452118f); 1164 const T NEG_PIO4F = pset1<T>(-0.7853981633974483096f); /* -pi / 4 */ 1165 T y = pabs(x); 1166 T z = pmul(y, y); 1167 T y_le_two = pselect( 1168 pcmp_lt(y, pset1<T>(1.0e-3f)), 1169 pmadd(z, pset1<T>(-0.25f), pset1<T>(1.0f)), 1170 pmul(psub(z, DR1), internal::ppolevl<T, 4>::run(z, JP))); 1171 T q = pdiv(pset1<T>(1.0f), y); 1172 T w = prsqrt(y); 1173 T p = pmul(w, internal::ppolevl<T, 7>::run(q, MO)); 1174 w = pmul(q, q); 1175 T yn = pmadd(q, internal::ppolevl<T, 7>::run(w, PH), NEG_PIO4F); 1176 T y_gt_two = pmul(p, pcos(padd(yn, y))); 1177 return pselect(pcmp_le(y, pset1<T>(2.0)), y_le_two, y_gt_two); 1178 } 1179 }; 1180 1181 template <typename T> 1182 struct generic_j0<T, double> { 1183 EIGEN_DEVICE_FUNC 1184 static EIGEN_STRONG_INLINE T run(const T& x) { 1185 /* j0.c 1186 * Bessel function of order zero 1187 * 1188 * 1189 * 1190 * SYNOPSIS: 1191 * 1192 * double x, y, j0(); 1193 * 1194 * y = j0( x ); 1195 * 1196 * 1197 * 1198 * DESCRIPTION: 1199 * 1200 * Returns Bessel function of order zero of the argument. 1201 * 1202 * The domain is divided into the intervals [0, 5] and 1203 * (5, infinity). In the first interval the following rational 1204 * approximation is used: 1205 * 1206 * 1207 * 2 2 1208 * (w - r ) (w - r ) P (w) / Q (w) 1209 * 1 2 3 8 1210 * 1211 * 2 1212 * where w = x and the two r's are zeros of the function. 1213 * 1214 * In the second interval, the Hankel asymptotic expansion 1215 * is employed with two rational functions of degree 6/6 1216 * and 7/7. 1217 * 1218 * 1219 * 1220 * ACCURACY: 1221 * 1222 * Absolute error: 1223 * arithmetic domain # trials peak rms 1224 * DEC 0, 30 10000 4.4e-17 6.3e-18 1225 * IEEE 0, 30 60000 4.2e-16 1.1e-16 1226 * 1227 */ 1228 const double PP[] = {7.96936729297347051624E-4, 8.28352392107440799803E-2, 1229 1.23953371646414299388E0, 5.44725003058768775090E0, 1230 8.74716500199817011941E0, 5.30324038235394892183E0, 1231 9.99999999999999997821E-1}; 1232 const double PQ[] = {9.24408810558863637013E-4, 8.56288474354474431428E-2, 1233 1.25352743901058953537E0, 5.47097740330417105182E0, 1234 8.76190883237069594232E0, 5.30605288235394617618E0, 1235 1.00000000000000000218E0}; 1236 const double QP[] = {-1.13663838898469149931E-2, -1.28252718670509318512E0, 1237 -1.95539544257735972385E1, -9.32060152123768231369E1, 1238 -1.77681167980488050595E2, -1.47077505154951170175E2, 1239 -5.14105326766599330220E1, -6.05014350600728481186E0}; 1240 const double QQ[] = {1.00000000000000000000E0, 6.43178256118178023184E1, 1241 8.56430025976980587198E2, 3.88240183605401609683E3, 1242 7.24046774195652478189E3, 5.93072701187316984827E3, 1243 2.06209331660327847417E3, 2.42005740240291393179E2}; 1244 const double RP[] = {-4.79443220978201773821E9, 1.95617491946556577543E12, 1245 -2.49248344360967716204E14, 9.70862251047306323952E15}; 1246 const double RQ[] = {1.00000000000000000000E0, 4.99563147152651017219E2, 1247 1.73785401676374683123E5, 4.84409658339962045305E7, 1248 1.11855537045356834862E10, 2.11277520115489217587E12, 1249 3.10518229857422583814E14, 3.18121955943204943306E16, 1250 1.71086294081043136091E18}; 1251 const T DR1 = pset1<T>(5.78318596294678452118E0); 1252 const T DR2 = pset1<T>(3.04712623436620863991E1); 1253 const T SQ2OPI = pset1<T>(7.9788456080286535587989E-1); /* sqrt(2 / pi) */ 1254 const T NEG_PIO4 = pset1<T>(-0.7853981633974483096); /* pi / 4 */ 1255 1256 T y = pabs(x); 1257 T z = pmul(y, y); 1258 T y_le_five = pselect( 1259 pcmp_lt(y, pset1<T>(1.0e-5)), 1260 pmadd(z, pset1<T>(-0.25), pset1<T>(1.0)), 1261 pmul(pmul(psub(z, DR1), psub(z, DR2)), 1262 pdiv(internal::ppolevl<T, 3>::run(z, RP), 1263 internal::ppolevl<T, 8>::run(z, RQ)))); 1264 T s = pdiv(pset1<T>(25.0), z); 1265 T p = pdiv( 1266 internal::ppolevl<T, 6>::run(s, PP), 1267 internal::ppolevl<T, 6>::run(s, PQ)); 1268 T q = pdiv( 1269 internal::ppolevl<T, 7>::run(s, QP), 1270 internal::ppolevl<T, 7>::run(s, QQ)); 1271 T yn = padd(y, NEG_PIO4); 1272 T w = pdiv(pset1<T>(-5.0), y); 1273 p = pmadd(p, pcos(yn), pmul(w, pmul(q, psin(yn)))); 1274 T y_gt_five = pmul(p, pmul(SQ2OPI, prsqrt(y))); 1275 return pselect(pcmp_le(y, pset1<T>(5.0)), y_le_five, y_gt_five); 1276 } 1277 }; 1278 1279 template <typename T> 1280 struct bessel_j0_impl { 1281 EIGEN_DEVICE_FUNC 1282 static EIGEN_STRONG_INLINE T run(const T x) { 1283 return generic_j0<T>::run(x); 1284 } 1285 }; 1286 1287 template <typename T> 1288 struct bessel_y0_retval { 1289 typedef T type; 1290 }; 1291 1292 template <typename T, typename ScalarType = typename unpacket_traits<T>::type> 1293 struct generic_y0 { 1294 EIGEN_DEVICE_FUNC 1295 static EIGEN_STRONG_INLINE T run(const T&) { 1296 EIGEN_STATIC_ASSERT((internal::is_same<T, T>::value == false), 1297 THIS_TYPE_IS_NOT_SUPPORTED); 1298 return ScalarType(0); 1299 } 1300 }; 1301 1302 template <typename T> 1303 struct generic_y0<T, float> { 1304 EIGEN_DEVICE_FUNC 1305 static EIGEN_STRONG_INLINE T run(const T& x) { 1306 /* j0f.c 1307 * Bessel function of the second kind, order zero 1308 * 1309 * 1310 * 1311 * SYNOPSIS: 1312 * 1313 * float x, y, y0f(); 1314 * 1315 * y = y0f( x ); 1316 * 1317 * 1318 * 1319 * DESCRIPTION: 1320 * 1321 * Returns Bessel function of the second kind, of order 1322 * zero, of the argument. 1323 * 1324 * The domain is divided into the intervals [0, 2] and 1325 * (2, infinity). In the first interval a rational approximation 1326 * R(x) is employed to compute 1327 * 1328 * 2 2 2 1329 * y0(x) = (w - r ) (w - r ) (w - r ) R(x) + 2/pi ln(x) j0(x). 1330 * 1 2 3 1331 * 1332 * Thus a call to j0() is required. The three zeros are removed 1333 * from R(x) to improve its numerical stability. 1334 * 1335 * In the second interval, the modulus and phase are approximated 1336 * by polynomials of the form Modulus(x) = sqrt(1/x) Q(1/x) 1337 * and Phase(x) = x + 1/x S(1/x^2) - pi/4. Then the function is 1338 * 1339 * y0(x) = Modulus(x) sin( Phase(x) ). 1340 * 1341 * 1342 * 1343 * 1344 * ACCURACY: 1345 * 1346 * Absolute error, when y0(x) < 1; else relative error: 1347 * 1348 * arithmetic domain # trials peak rms 1349 * IEEE 0, 2 100000 2.4e-7 3.4e-8 1350 * IEEE 2, 32 100000 1.8e-7 5.3e-8 1351 * 1352 */ 1353 1354 const float YP[] = {9.454583683980369E-008f, -9.413212653797057E-006f, 1355 5.344486707214273E-004f, -1.584289289821316E-002f, 1356 1.707584643733568E-001f}; 1357 const float MO[] = {-6.838999669318810E-002f, 1.864949361379502E-001f, 1358 -2.145007480346739E-001f, 1.197549369473540E-001f, 1359 -3.560281861530129E-003f, -4.969382655296620E-002f, 1360 -3.355424622293709E-006f, 7.978845717621440E-001f}; 1361 const float PH[] = {3.242077816988247E+001f, -3.630592630518434E+001f, 1362 1.756221482109099E+001f, -4.974978466280903E+000f, 1363 1.001973420681837E+000f, -1.939906941791308E-001f, 1364 6.490598792654666E-002f, -1.249992184872738E-001f}; 1365 const T YZ1 = pset1<T>(0.43221455686510834878f); 1366 const T TWOOPI = pset1<T>(0.636619772367581343075535f); /* 2 / pi */ 1367 const T NEG_PIO4F = pset1<T>(-0.7853981633974483096f); /* -pi / 4 */ 1368 const T NEG_MAXNUM = pset1<T>(-NumTraits<float>::infinity()); 1369 T z = pmul(x, x); 1370 T x_le_two = pmul(TWOOPI, pmul(plog(x), generic_j0<T, float>::run(x))); 1371 x_le_two = pmadd( 1372 psub(z, YZ1), internal::ppolevl<T, 4>::run(z, YP), x_le_two); 1373 x_le_two = pselect(pcmp_le(x, pset1<T>(0.0)), NEG_MAXNUM, x_le_two); 1374 T q = pdiv(pset1<T>(1.0), x); 1375 T w = prsqrt(x); 1376 T p = pmul(w, internal::ppolevl<T, 7>::run(q, MO)); 1377 T u = pmul(q, q); 1378 T xn = pmadd(q, internal::ppolevl<T, 7>::run(u, PH), NEG_PIO4F); 1379 T x_gt_two = pmul(p, psin(padd(xn, x))); 1380 return pselect(pcmp_le(x, pset1<T>(2.0)), x_le_two, x_gt_two); 1381 } 1382 }; 1383 1384 template <typename T> 1385 struct generic_y0<T, double> { 1386 EIGEN_DEVICE_FUNC 1387 static EIGEN_STRONG_INLINE T run(const T& x) { 1388 /* j0.c 1389 * Bessel function of the second kind, order zero 1390 * 1391 * 1392 * 1393 * SYNOPSIS: 1394 * 1395 * double x, y, y0(); 1396 * 1397 * y = y0( x ); 1398 * 1399 * 1400 * 1401 * DESCRIPTION: 1402 * 1403 * Returns Bessel function of the second kind, of order 1404 * zero, of the argument. 1405 * 1406 * The domain is divided into the intervals [0, 5] and 1407 * (5, infinity). In the first interval a rational approximation 1408 * R(x) is employed to compute 1409 * y0(x) = R(x) + 2 * log(x) * j0(x) / PI. 1410 * Thus a call to j0() is required. 1411 * 1412 * In the second interval, the Hankel asymptotic expansion 1413 * is employed with two rational functions of degree 6/6 1414 * and 7/7. 1415 * 1416 * 1417 * 1418 * ACCURACY: 1419 * 1420 * Absolute error, when y0(x) < 1; else relative error: 1421 * 1422 * arithmetic domain # trials peak rms 1423 * DEC 0, 30 9400 7.0e-17 7.9e-18 1424 * IEEE 0, 30 30000 1.3e-15 1.6e-16 1425 * 1426 */ 1427 const double PP[] = {7.96936729297347051624E-4, 8.28352392107440799803E-2, 1428 1.23953371646414299388E0, 5.44725003058768775090E0, 1429 8.74716500199817011941E0, 5.30324038235394892183E0, 1430 9.99999999999999997821E-1}; 1431 const double PQ[] = {9.24408810558863637013E-4, 8.56288474354474431428E-2, 1432 1.25352743901058953537E0, 5.47097740330417105182E0, 1433 8.76190883237069594232E0, 5.30605288235394617618E0, 1434 1.00000000000000000218E0}; 1435 const double QP[] = {-1.13663838898469149931E-2, -1.28252718670509318512E0, 1436 -1.95539544257735972385E1, -9.32060152123768231369E1, 1437 -1.77681167980488050595E2, -1.47077505154951170175E2, 1438 -5.14105326766599330220E1, -6.05014350600728481186E0}; 1439 const double QQ[] = {1.00000000000000000000E0, 6.43178256118178023184E1, 1440 8.56430025976980587198E2, 3.88240183605401609683E3, 1441 7.24046774195652478189E3, 5.93072701187316984827E3, 1442 2.06209331660327847417E3, 2.42005740240291393179E2}; 1443 const double YP[] = {1.55924367855235737965E4, -1.46639295903971606143E7, 1444 5.43526477051876500413E9, -9.82136065717911466409E11, 1445 8.75906394395366999549E13, -3.46628303384729719441E15, 1446 4.42733268572569800351E16, -1.84950800436986690637E16}; 1447 const double YQ[] = {1.00000000000000000000E0, 1.04128353664259848412E3, 1448 6.26107330137134956842E5, 2.68919633393814121987E8, 1449 8.64002487103935000337E10, 2.02979612750105546709E13, 1450 3.17157752842975028269E15, 2.50596256172653059228E17}; 1451 const T SQ2OPI = pset1<T>(7.9788456080286535587989E-1); /* sqrt(2 / pi) */ 1452 const T TWOOPI = pset1<T>(0.636619772367581343075535); /* 2 / pi */ 1453 const T NEG_PIO4 = pset1<T>(-0.7853981633974483096); /* -pi / 4 */ 1454 const T NEG_MAXNUM = pset1<T>(-NumTraits<double>::infinity()); 1455 1456 T z = pmul(x, x); 1457 T x_le_five = pdiv(internal::ppolevl<T, 7>::run(z, YP), 1458 internal::ppolevl<T, 7>::run(z, YQ)); 1459 x_le_five = pmadd( 1460 pmul(TWOOPI, plog(x)), generic_j0<T, double>::run(x), x_le_five); 1461 x_le_five = pselect(pcmp_le(x, pset1<T>(0.0)), NEG_MAXNUM, x_le_five); 1462 T s = pdiv(pset1<T>(25.0), z); 1463 T p = pdiv( 1464 internal::ppolevl<T, 6>::run(s, PP), 1465 internal::ppolevl<T, 6>::run(s, PQ)); 1466 T q = pdiv( 1467 internal::ppolevl<T, 7>::run(s, QP), 1468 internal::ppolevl<T, 7>::run(s, QQ)); 1469 T xn = padd(x, NEG_PIO4); 1470 T w = pdiv(pset1<T>(5.0), x); 1471 p = pmadd(p, psin(xn), pmul(w, pmul(q, pcos(xn)))); 1472 T x_gt_five = pmul(p, pmul(SQ2OPI, prsqrt(x))); 1473 return pselect(pcmp_le(x, pset1<T>(5.0)), x_le_five, x_gt_five); 1474 } 1475 }; 1476 1477 template <typename T> 1478 struct bessel_y0_impl { 1479 EIGEN_DEVICE_FUNC 1480 static EIGEN_STRONG_INLINE T run(const T x) { 1481 return generic_y0<T>::run(x); 1482 } 1483 }; 1484 1485 template <typename T> 1486 struct bessel_j1_retval { 1487 typedef T type; 1488 }; 1489 1490 template <typename T, typename ScalarType = typename unpacket_traits<T>::type> 1491 struct generic_j1 { 1492 EIGEN_DEVICE_FUNC 1493 static EIGEN_STRONG_INLINE T run(const T&) { 1494 EIGEN_STATIC_ASSERT((internal::is_same<T, T>::value == false), 1495 THIS_TYPE_IS_NOT_SUPPORTED); 1496 return ScalarType(0); 1497 } 1498 }; 1499 1500 template <typename T> 1501 struct generic_j1<T, float> { 1502 EIGEN_DEVICE_FUNC 1503 static EIGEN_STRONG_INLINE T run(const T& x) { 1504 /* j1f.c 1505 * Bessel function of order one 1506 * 1507 * 1508 * 1509 * SYNOPSIS: 1510 * 1511 * float x, y, j1f(); 1512 * 1513 * y = j1f( x ); 1514 * 1515 * 1516 * 1517 * DESCRIPTION: 1518 * 1519 * Returns Bessel function of order one of the argument. 1520 * 1521 * The domain is divided into the intervals [0, 2] and 1522 * (2, infinity). In the first interval a polynomial approximation 1523 * 2 1524 * (w - r ) x P(w) 1525 * 1 1526 * 2 1527 * is used, where w = x and r is the first zero of the function. 1528 * 1529 * In the second interval, the modulus and phase are approximated 1530 * by polynomials of the form Modulus(x) = sqrt(1/x) Q(1/x) 1531 * and Phase(x) = x + 1/x R(1/x^2) - 3pi/4. The function is 1532 * 1533 * j0(x) = Modulus(x) cos( Phase(x) ). 1534 * 1535 * 1536 * 1537 * ACCURACY: 1538 * 1539 * Absolute error: 1540 * arithmetic domain # trials peak rms 1541 * IEEE 0, 2 100000 1.2e-7 2.5e-8 1542 * IEEE 2, 32 100000 2.0e-7 5.3e-8 1543 * 1544 * 1545 */ 1546 1547 const float JP[] = {-4.878788132172128E-009f, 6.009061827883699E-007f, 1548 -4.541343896997497E-005f, 1.937383947804541E-003f, 1549 -3.405537384615824E-002f}; 1550 const float MO1[] = {6.913942741265801E-002f, -2.284801500053359E-001f, 1551 3.138238455499697E-001f, -2.102302420403875E-001f, 1552 5.435364690523026E-003f, 1.493389585089498E-001f, 1553 4.976029650847191E-006f, 7.978845453073848E-001f}; 1554 const float PH1[] = {-4.497014141919556E+001f, 5.073465654089319E+001f, 1555 -2.485774108720340E+001f, 7.222973196770240E+000f, 1556 -1.544842782180211E+000f, 3.503787691653334E-001f, 1557 -1.637986776941202E-001f, 3.749989509080821E-001f}; 1558 const T Z1 = pset1<T>(1.46819706421238932572E1f); 1559 const T NEG_THPIO4F = pset1<T>(-2.35619449019234492885f); /* -3*pi/4 */ 1560 1561 T y = pabs(x); 1562 T z = pmul(y, y); 1563 T y_le_two = pmul( 1564 psub(z, Z1), 1565 pmul(x, internal::ppolevl<T, 4>::run(z, JP))); 1566 T q = pdiv(pset1<T>(1.0f), y); 1567 T w = prsqrt(y); 1568 T p = pmul(w, internal::ppolevl<T, 7>::run(q, MO1)); 1569 w = pmul(q, q); 1570 T yn = pmadd(q, internal::ppolevl<T, 7>::run(w, PH1), NEG_THPIO4F); 1571 T y_gt_two = pmul(p, pcos(padd(yn, y))); 1572 // j1 is an odd function. This implementation differs from cephes to 1573 // take this fact in to account. Cephes returns -j1(x) for y > 2 range. 1574 y_gt_two = pselect( 1575 pcmp_lt(x, pset1<T>(0.0f)), pnegate(y_gt_two), y_gt_two); 1576 return pselect(pcmp_le(y, pset1<T>(2.0f)), y_le_two, y_gt_two); 1577 } 1578 }; 1579 1580 template <typename T> 1581 struct generic_j1<T, double> { 1582 EIGEN_DEVICE_FUNC 1583 static EIGEN_STRONG_INLINE T run(const T& x) { 1584 /* j1.c 1585 * Bessel function of order one 1586 * 1587 * 1588 * 1589 * SYNOPSIS: 1590 * 1591 * double x, y, j1(); 1592 * 1593 * y = j1( x ); 1594 * 1595 * 1596 * 1597 * DESCRIPTION: 1598 * 1599 * Returns Bessel function of order one of the argument. 1600 * 1601 * The domain is divided into the intervals [0, 8] and 1602 * (8, infinity). In the first interval a 24 term Chebyshev 1603 * expansion is used. In the second, the asymptotic 1604 * trigonometric representation is employed using two 1605 * rational functions of degree 5/5. 1606 * 1607 * 1608 * 1609 * ACCURACY: 1610 * 1611 * Absolute error: 1612 * arithmetic domain # trials peak rms 1613 * DEC 0, 30 10000 4.0e-17 1.1e-17 1614 * IEEE 0, 30 30000 2.6e-16 1.1e-16 1615 * 1616 */ 1617 const double PP[] = {7.62125616208173112003E-4, 7.31397056940917570436E-2, 1618 1.12719608129684925192E0, 5.11207951146807644818E0, 1619 8.42404590141772420927E0, 5.21451598682361504063E0, 1620 1.00000000000000000254E0}; 1621 const double PQ[] = {5.71323128072548699714E-4, 6.88455908754495404082E-2, 1622 1.10514232634061696926E0, 5.07386386128601488557E0, 1623 8.39985554327604159757E0, 5.20982848682361821619E0, 1624 9.99999999999999997461E-1}; 1625 const double QP[] = {5.10862594750176621635E-2, 4.98213872951233449420E0, 1626 7.58238284132545283818E1, 3.66779609360150777800E2, 1627 7.10856304998926107277E2, 5.97489612400613639965E2, 1628 2.11688757100572135698E2, 2.52070205858023719784E1}; 1629 const double QQ[] = {1.00000000000000000000E0, 7.42373277035675149943E1, 1630 1.05644886038262816351E3, 4.98641058337653607651E3, 1631 9.56231892404756170795E3, 7.99704160447350683650E3, 1632 2.82619278517639096600E3, 3.36093607810698293419E2}; 1633 const double RP[] = {-8.99971225705559398224E8, 4.52228297998194034323E11, 1634 -7.27494245221818276015E13, 3.68295732863852883286E15}; 1635 const double RQ[] = {1.00000000000000000000E0, 6.20836478118054335476E2, 1636 2.56987256757748830383E5, 8.35146791431949253037E7, 1637 2.21511595479792499675E10, 4.74914122079991414898E12, 1638 7.84369607876235854894E14, 8.95222336184627338078E16, 1639 5.32278620332680085395E18}; 1640 const T Z1 = pset1<T>(1.46819706421238932572E1); 1641 const T Z2 = pset1<T>(4.92184563216946036703E1); 1642 const T NEG_THPIO4 = pset1<T>(-2.35619449019234492885); /* -3*pi/4 */ 1643 const T SQ2OPI = pset1<T>(7.9788456080286535587989E-1); /* sqrt(2 / pi) */ 1644 T y = pabs(x); 1645 T z = pmul(y, y); 1646 T y_le_five = pdiv(internal::ppolevl<T, 3>::run(z, RP), 1647 internal::ppolevl<T, 8>::run(z, RQ)); 1648 y_le_five = pmul(pmul(pmul(y_le_five, x), psub(z, Z1)), psub(z, Z2)); 1649 T s = pdiv(pset1<T>(25.0), z); 1650 T p = pdiv( 1651 internal::ppolevl<T, 6>::run(s, PP), 1652 internal::ppolevl<T, 6>::run(s, PQ)); 1653 T q = pdiv( 1654 internal::ppolevl<T, 7>::run(s, QP), 1655 internal::ppolevl<T, 7>::run(s, QQ)); 1656 T yn = padd(y, NEG_THPIO4); 1657 T w = pdiv(pset1<T>(-5.0), y); 1658 p = pmadd(p, pcos(yn), pmul(w, pmul(q, psin(yn)))); 1659 T y_gt_five = pmul(p, pmul(SQ2OPI, prsqrt(y))); 1660 // j1 is an odd function. This implementation differs from cephes to 1661 // take this fact in to account. Cephes returns -j1(x) for y > 5 range. 1662 y_gt_five = pselect( 1663 pcmp_lt(x, pset1<T>(0.0)), pnegate(y_gt_five), y_gt_five); 1664 return pselect(pcmp_le(y, pset1<T>(5.0)), y_le_five, y_gt_five); 1665 } 1666 }; 1667 1668 template <typename T> 1669 struct bessel_j1_impl { 1670 EIGEN_DEVICE_FUNC 1671 static EIGEN_STRONG_INLINE T run(const T x) { 1672 return generic_j1<T>::run(x); 1673 } 1674 }; 1675 1676 template <typename T> 1677 struct bessel_y1_retval { 1678 typedef T type; 1679 }; 1680 1681 template <typename T, typename ScalarType = typename unpacket_traits<T>::type> 1682 struct generic_y1 { 1683 EIGEN_DEVICE_FUNC 1684 static EIGEN_STRONG_INLINE T run(const T&) { 1685 EIGEN_STATIC_ASSERT((internal::is_same<T, T>::value == false), 1686 THIS_TYPE_IS_NOT_SUPPORTED); 1687 return ScalarType(0); 1688 } 1689 }; 1690 1691 template <typename T> 1692 struct generic_y1<T, float> { 1693 EIGEN_DEVICE_FUNC 1694 static EIGEN_STRONG_INLINE T run(const T& x) { 1695 /* j1f.c 1696 * Bessel function of second kind of order one 1697 * 1698 * 1699 * 1700 * SYNOPSIS: 1701 * 1702 * double x, y, y1(); 1703 * 1704 * y = y1( x ); 1705 * 1706 * 1707 * 1708 * DESCRIPTION: 1709 * 1710 * Returns Bessel function of the second kind of order one 1711 * of the argument. 1712 * 1713 * The domain is divided into the intervals [0, 2] and 1714 * (2, infinity). In the first interval a rational approximation 1715 * R(x) is employed to compute 1716 * 1717 * 2 1718 * y0(x) = (w - r ) x R(x^2) + 2/pi (ln(x) j1(x) - 1/x) . 1719 * 1 1720 * 1721 * Thus a call to j1() is required. 1722 * 1723 * In the second interval, the modulus and phase are approximated 1724 * by polynomials of the form Modulus(x) = sqrt(1/x) Q(1/x) 1725 * and Phase(x) = x + 1/x S(1/x^2) - 3pi/4. Then the function is 1726 * 1727 * y0(x) = Modulus(x) sin( Phase(x) ). 1728 * 1729 * 1730 * 1731 * 1732 * ACCURACY: 1733 * 1734 * Absolute error: 1735 * arithmetic domain # trials peak rms 1736 * IEEE 0, 2 100000 2.2e-7 4.6e-8 1737 * IEEE 2, 32 100000 1.9e-7 5.3e-8 1738 * 1739 * (error criterion relative when |y1| > 1). 1740 * 1741 */ 1742 1743 const float YP[] = {8.061978323326852E-009f, -9.496460629917016E-007f, 1744 6.719543806674249E-005f, -2.641785726447862E-003f, 1745 4.202369946500099E-002f}; 1746 const float MO1[] = {6.913942741265801E-002f, -2.284801500053359E-001f, 1747 3.138238455499697E-001f, -2.102302420403875E-001f, 1748 5.435364690523026E-003f, 1.493389585089498E-001f, 1749 4.976029650847191E-006f, 7.978845453073848E-001f}; 1750 const float PH1[] = {-4.497014141919556E+001f, 5.073465654089319E+001f, 1751 -2.485774108720340E+001f, 7.222973196770240E+000f, 1752 -1.544842782180211E+000f, 3.503787691653334E-001f, 1753 -1.637986776941202E-001f, 3.749989509080821E-001f}; 1754 const T YO1 = pset1<T>(4.66539330185668857532f); 1755 const T NEG_THPIO4F = pset1<T>(-2.35619449019234492885f); /* -3*pi/4 */ 1756 const T TWOOPI = pset1<T>(0.636619772367581343075535f); /* 2/pi */ 1757 const T NEG_MAXNUM = pset1<T>(-NumTraits<float>::infinity()); 1758 1759 T z = pmul(x, x); 1760 T x_le_two = pmul(psub(z, YO1), internal::ppolevl<T, 4>::run(z, YP)); 1761 x_le_two = pmadd( 1762 x_le_two, x, 1763 pmul(TWOOPI, pmadd( 1764 generic_j1<T, float>::run(x), plog(x), 1765 pdiv(pset1<T>(-1.0f), x)))); 1766 x_le_two = pselect(pcmp_lt(x, pset1<T>(0.0f)), NEG_MAXNUM, x_le_two); 1767 1768 T q = pdiv(pset1<T>(1.0), x); 1769 T w = prsqrt(x); 1770 T p = pmul(w, internal::ppolevl<T, 7>::run(q, MO1)); 1771 w = pmul(q, q); 1772 T xn = pmadd(q, internal::ppolevl<T, 7>::run(w, PH1), NEG_THPIO4F); 1773 T x_gt_two = pmul(p, psin(padd(xn, x))); 1774 return pselect(pcmp_le(x, pset1<T>(2.0)), x_le_two, x_gt_two); 1775 } 1776 }; 1777 1778 template <typename T> 1779 struct generic_y1<T, double> { 1780 EIGEN_DEVICE_FUNC 1781 static EIGEN_STRONG_INLINE T run(const T& x) { 1782 /* j1.c 1783 * Bessel function of second kind of order one 1784 * 1785 * 1786 * 1787 * SYNOPSIS: 1788 * 1789 * double x, y, y1(); 1790 * 1791 * y = y1( x ); 1792 * 1793 * 1794 * 1795 * DESCRIPTION: 1796 * 1797 * Returns Bessel function of the second kind of order one 1798 * of the argument. 1799 * 1800 * The domain is divided into the intervals [0, 8] and 1801 * (8, infinity). In the first interval a 25 term Chebyshev 1802 * expansion is used, and a call to j1() is required. 1803 * In the second, the asymptotic trigonometric representation 1804 * is employed using two rational functions of degree 5/5. 1805 * 1806 * 1807 * 1808 * ACCURACY: 1809 * 1810 * Absolute error: 1811 * arithmetic domain # trials peak rms 1812 * DEC 0, 30 10000 8.6e-17 1.3e-17 1813 * IEEE 0, 30 30000 1.0e-15 1.3e-16 1814 * 1815 * (error criterion relative when |y1| > 1). 1816 * 1817 */ 1818 const double PP[] = {7.62125616208173112003E-4, 7.31397056940917570436E-2, 1819 1.12719608129684925192E0, 5.11207951146807644818E0, 1820 8.42404590141772420927E0, 5.21451598682361504063E0, 1821 1.00000000000000000254E0}; 1822 const double PQ[] = {5.71323128072548699714E-4, 6.88455908754495404082E-2, 1823 1.10514232634061696926E0, 5.07386386128601488557E0, 1824 8.39985554327604159757E0, 5.20982848682361821619E0, 1825 9.99999999999999997461E-1}; 1826 const double QP[] = {5.10862594750176621635E-2, 4.98213872951233449420E0, 1827 7.58238284132545283818E1, 3.66779609360150777800E2, 1828 7.10856304998926107277E2, 5.97489612400613639965E2, 1829 2.11688757100572135698E2, 2.52070205858023719784E1}; 1830 const double QQ[] = {1.00000000000000000000E0, 7.42373277035675149943E1, 1831 1.05644886038262816351E3, 4.98641058337653607651E3, 1832 9.56231892404756170795E3, 7.99704160447350683650E3, 1833 2.82619278517639096600E3, 3.36093607810698293419E2}; 1834 const double YP[] = {1.26320474790178026440E9, -6.47355876379160291031E11, 1835 1.14509511541823727583E14, -8.12770255501325109621E15, 1836 2.02439475713594898196E17, -7.78877196265950026825E17}; 1837 const double YQ[] = {1.00000000000000000000E0, 5.94301592346128195359E2, 1838 2.35564092943068577943E5, 7.34811944459721705660E7, 1839 1.87601316108706159478E10, 3.88231277496238566008E12, 1840 6.20557727146953693363E14, 6.87141087355300489866E16, 1841 3.97270608116560655612E18}; 1842 const T SQ2OPI = pset1<T>(.79788456080286535588); 1843 const T NEG_THPIO4 = pset1<T>(-2.35619449019234492885); /* -3*pi/4 */ 1844 const T TWOOPI = pset1<T>(0.636619772367581343075535); /* 2/pi */ 1845 const T NEG_MAXNUM = pset1<T>(-NumTraits<double>::infinity()); 1846 1847 T z = pmul(x, x); 1848 T x_le_five = pdiv(internal::ppolevl<T, 5>::run(z, YP), 1849 internal::ppolevl<T, 8>::run(z, YQ)); 1850 x_le_five = pmadd( 1851 x_le_five, x, pmul( 1852 TWOOPI, pmadd(generic_j1<T, double>::run(x), plog(x), 1853 pdiv(pset1<T>(-1.0), x)))); 1854 1855 x_le_five = pselect(pcmp_le(x, pset1<T>(0.0)), NEG_MAXNUM, x_le_five); 1856 T s = pdiv(pset1<T>(25.0), z); 1857 T p = pdiv( 1858 internal::ppolevl<T, 6>::run(s, PP), 1859 internal::ppolevl<T, 6>::run(s, PQ)); 1860 T q = pdiv( 1861 internal::ppolevl<T, 7>::run(s, QP), 1862 internal::ppolevl<T, 7>::run(s, QQ)); 1863 T xn = padd(x, NEG_THPIO4); 1864 T w = pdiv(pset1<T>(5.0), x); 1865 p = pmadd(p, psin(xn), pmul(w, pmul(q, pcos(xn)))); 1866 T x_gt_five = pmul(p, pmul(SQ2OPI, prsqrt(x))); 1867 return pselect(pcmp_le(x, pset1<T>(5.0)), x_le_five, x_gt_five); 1868 } 1869 }; 1870 1871 template <typename T> 1872 struct bessel_y1_impl { 1873 EIGEN_DEVICE_FUNC 1874 static EIGEN_STRONG_INLINE T run(const T x) { 1875 return generic_y1<T>::run(x); 1876 } 1877 }; 1878 1879 } // end namespace internal 1880 1881 namespace numext { 1882 1883 template <typename Scalar> 1884 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(bessel_i0, Scalar) 1885 bessel_i0(const Scalar& x) { 1886 return EIGEN_MATHFUNC_IMPL(bessel_i0, Scalar)::run(x); 1887 } 1888 1889 template <typename Scalar> 1890 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(bessel_i0e, Scalar) 1891 bessel_i0e(const Scalar& x) { 1892 return EIGEN_MATHFUNC_IMPL(bessel_i0e, Scalar)::run(x); 1893 } 1894 1895 template <typename Scalar> 1896 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(bessel_i1, Scalar) 1897 bessel_i1(const Scalar& x) { 1898 return EIGEN_MATHFUNC_IMPL(bessel_i1, Scalar)::run(x); 1899 } 1900 1901 template <typename Scalar> 1902 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(bessel_i1e, Scalar) 1903 bessel_i1e(const Scalar& x) { 1904 return EIGEN_MATHFUNC_IMPL(bessel_i1e, Scalar)::run(x); 1905 } 1906 1907 template <typename Scalar> 1908 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(bessel_k0, Scalar) 1909 bessel_k0(const Scalar& x) { 1910 return EIGEN_MATHFUNC_IMPL(bessel_k0, Scalar)::run(x); 1911 } 1912 1913 template <typename Scalar> 1914 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(bessel_k0e, Scalar) 1915 bessel_k0e(const Scalar& x) { 1916 return EIGEN_MATHFUNC_IMPL(bessel_k0e, Scalar)::run(x); 1917 } 1918 1919 template <typename Scalar> 1920 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(bessel_k1, Scalar) 1921 bessel_k1(const Scalar& x) { 1922 return EIGEN_MATHFUNC_IMPL(bessel_k1, Scalar)::run(x); 1923 } 1924 1925 template <typename Scalar> 1926 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(bessel_k1e, Scalar) 1927 bessel_k1e(const Scalar& x) { 1928 return EIGEN_MATHFUNC_IMPL(bessel_k1e, Scalar)::run(x); 1929 } 1930 1931 template <typename Scalar> 1932 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(bessel_j0, Scalar) 1933 bessel_j0(const Scalar& x) { 1934 return EIGEN_MATHFUNC_IMPL(bessel_j0, Scalar)::run(x); 1935 } 1936 1937 template <typename Scalar> 1938 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(bessel_y0, Scalar) 1939 bessel_y0(const Scalar& x) { 1940 return EIGEN_MATHFUNC_IMPL(bessel_y0, Scalar)::run(x); 1941 } 1942 1943 template <typename Scalar> 1944 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(bessel_j1, Scalar) 1945 bessel_j1(const Scalar& x) { 1946 return EIGEN_MATHFUNC_IMPL(bessel_j1, Scalar)::run(x); 1947 } 1948 1949 template <typename Scalar> 1950 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(bessel_y1, Scalar) 1951 bessel_y1(const Scalar& x) { 1952 return EIGEN_MATHFUNC_IMPL(bessel_y1, Scalar)::run(x); 1953 } 1954 1955 } // end namespace numext 1956 1957 } // end namespace Eigen 1958 1959 #endif // EIGEN_BESSEL_FUNCTIONS_H 1960