Line data    Source code

       1             : #pragma once
       2             : 
       3             : #include "CUDADefines.h"
       4             : #include <cmath>
       5             : 
       6             : namespace elsa::fn
       7             : {
       8             :     namespace detail
       9             :     {
      10             :         struct BesselFn_log_0 {
      11             :             template <typename T>
      12             :             __host__ __device__ constexpr T operator()(const T& arg) const noexcept
      13         215 :             {
      14         215 :                 if (arg < static_cast<T>(3.75)) {
      15         164 :                     const T ratio = arg / static_cast<T>(3.75);
      16         164 :                     const T y = ratio * ratio;
      17             : 
      18         164 :                     const T p0 = static_cast<T>(0.45813e-2);
      19         164 :                     const T p1 = static_cast<T>(0.360768e-1) + y * p0;
      20         164 :                     const T p2 = static_cast<T>(0.2659732) + y * p1;
      21         164 :                     const T p3 = static_cast<T>(1.2067492) + y * p2;
      22         164 :                     const T p4 = static_cast<T>(3.0899424) + y * p3;
      23         164 :                     const T p5 = static_cast<T>(3.5156229) + y * p4;
      24         164 :                     const T p6 = static_cast<T>(1.0) + y * p5;
      25             : 
      26         164 :                     return std::log(p6);
      27         164 :                 } else {
      28             :                     // see Numerical Recipes in C - 2nd Edition
      29             :                     // by W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery
      30             :                     // p.237
      31          51 :                     const T y = static_cast<T>(3.75) / arg;
      32             : 
      33          51 :                     const T p0 = static_cast<T>(0.392377e-2);
      34          51 :                     const T p1 = static_cast<T>(-0.1647633e-1) + y * p0;
      35          51 :                     const T p2 = static_cast<T>(0.2635537e-1) + y * p1;
      36          51 :                     const T p3 = static_cast<T>(-0.2057706e-1) + y * p2;
      37          51 :                     const T p4 = static_cast<T>(0.916281e-2) + y * p3;
      38          51 :                     const T p5 = static_cast<T>(-0.157565e-2) + y * p4;
      39          51 :                     const T p6 = static_cast<T>(0.225319e-2) + y * p5;
      40          51 :                     const T p7 = static_cast<T>(0.1328592e-1) + y * p6;
      41          51 :                     const T p8 = static_cast<T>(0.39894228) + y * p7;
      42             : 
      43          51 :                     return arg - std::sqrt(arg) + std::log(p8);
      44          51 :                 }
      45         215 :             }
      46             :         };
      47             : 
      48             :         struct BesselFn_1_0 {
      49             :             template <typename T>
      50             :             __host__ __device__ constexpr T operator()(const T& arg) const noexcept
      51         399 :             {
      52         399 :                 if (arg < static_cast<T>(3.75)) {
      53         323 :                     T bessel_0_p6 = 0;
      54         323 :                     const T ratio = arg / static_cast<T>(3.75);
      55         323 :                     const T y = ratio * ratio;
      56             : 
      57         323 :                     {
      58         323 :                         const T p0 = static_cast<T>(0.45813e-2);
      59         323 :                         const T p1 = static_cast<T>(0.360768e-1) + y * p0;
      60         323 :                         const T p2 = static_cast<T>(0.2659732) + y * p1;
      61         323 :                         const T p3 = static_cast<T>(1.2067492) + y * p2;
      62         323 :                         const T p4 = static_cast<T>(3.0899424) + y * p3;
      63         323 :                         const T p5 = static_cast<T>(3.5156229) + y * p4;
      64         323 :                         bessel_0_p6 = static_cast<T>(1.0) + y * p5;
      65         323 :                     }
      66             : 
      67         323 :                     T bessel_1_p6 = 0;
      68         323 :                     {
      69         323 :                         const T p0 = static_cast<T>(0.32411e-3);
      70         323 :                         const T p1 = static_cast<T>(0.301532e-2) + y * p0;
      71         323 :                         const T p2 = static_cast<T>(0.2658733e-1) + y * p1;
      72         323 :                         const T p3 = static_cast<T>(0.15084934) + y * p2;
      73         323 :                         const T p4 = static_cast<T>(0.51498869) + y * p3;
      74         323 :                         const T p5 = static_cast<T>(0.87890594) + y * p4;
      75         323 :                         bessel_1_p6 = arg * (static_cast<T>(0.5) + y * p5);
      76         323 :                     }
      77             : 
      78         323 :                     return bessel_1_p6 / bessel_0_p6;
      79         323 :                 } else {
      80             :                     // see Numerical Recipes in C - 2nd Edition
      81             :                     // by W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery
      82             :                     // p.238
      83          76 :                     const T y = static_cast<T>(3.75) / arg;
      84             : 
      85          76 :                     T bessel_0_p8;
      86          76 :                     {
      87          76 :                         const T p0 = static_cast<T>(0.392377e-2);
      88          76 :                         const T p1 = static_cast<T>(-0.1647633e-1) + y * p0;
      89          76 :                         const T p2 = static_cast<T>(0.2635537e-1) + y * p1;
      90          76 :                         const T p3 = static_cast<T>(-0.2057706e-1) + y * p2;
      91          76 :                         const T p4 = static_cast<T>(0.916281e-2) + y * p3;
      92          76 :                         const T p5 = static_cast<T>(-0.157565e-2) + y * p4;
      93          76 :                         const T p6 = static_cast<T>(0.225319e-2) + y * p5;
      94          76 :                         const T p7 = static_cast<T>(0.1328592e-1) + y * p6;
      95          76 :                         bessel_0_p8 = static_cast<T>(0.39894228) + y * p7;
      96          76 :                     }
      97             : 
      98          76 :                     T bessel_1_p8;
      99          76 :                     {
     100          76 :                         const T p0 = static_cast<T>(0.420059e-2);
     101          76 :                         const T p1 = static_cast<T>(0.1787654e-1) - y * p0;
     102          76 :                         const T p2 = static_cast<T>(-0.2895312e-1) + y * p1;
     103          76 :                         const T p3 = static_cast<T>(0.2282967e-1) + y * p2;
     104          76 :                         const T p4 = static_cast<T>(-0.1031555e-1) + y * p3;
     105          76 :                         const T p5 = static_cast<T>(0.163801e-2) + y * p4;
     106          76 :                         const T p6 = static_cast<T>(-0.362018e-2) + y * p5;
     107          76 :                         const T p7 = static_cast<T>(-0.3988024e-1) + y * p6;
     108          76 :                         bessel_1_p8 = static_cast<T>(0.39894228) + y * p7;
     109          76 :                     }
     110             : 
     111          76 :                     return bessel_1_p8 / bessel_0_p8;
     112          76 :                 }
     113         399 :             }
     114             :         };
     115             :     } // namespace detail
     116             : 
     117             :     static constexpr __device__ detail::BesselFn_log_0 bessel_log_0;
     118             :     static constexpr __device__ detail::BesselFn_1_0 bessel_1_0;
     119             : } // namespace elsa::fn