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