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