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


gtsam
Author(s):
autogenerated on Sat Nov 16 2024 04:01:55