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         189 :             {
      14         189 :                 if (arg < static_cast<T>(3.75)) {
      15         169 :                     const T ratio = arg / static_cast<T>(3.75);
      16         169 :                     const T y = ratio * ratio;
      17             : 
      18         169 :                     const T p0 = static_cast<T>(0.45813e-2);
      19         169 :                     const T p1 = static_cast<T>(0.360768e-1) + y * p0;
      20         169 :                     const T p2 = static_cast<T>(0.2659732) + y * p1;
      21         169 :                     const T p3 = static_cast<T>(1.2067492) + y * p2;
      22         169 :                     const T p4 = static_cast<T>(3.0899424) + y * p3;
      23         169 :                     const T p5 = static_cast<T>(3.5156229) + y * p4;
      24         169 :                     const T p6 = static_cast<T>(1.0) + y * p5;
      25             : 
      26         169 :                     return std::log(p6);
      27         169 :                 } 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          20 :                     const T y = static_cast<T>(3.75) / arg;
      32             : 
      33          20 :                     const T p0 = static_cast<T>(0.392377e-2);
      34          20 :                     const T p1 = static_cast<T>(-0.1647633e-1) + y * p0;
      35          20 :                     const T p2 = static_cast<T>(0.2635537e-1) + y * p1;
      36          20 :                     const T p3 = static_cast<T>(-0.2057706e-1) + y * p2;
      37          20 :                     const T p4 = static_cast<T>(0.916281e-2) + y * p3;
      38          20 :                     const T p5 = static_cast<T>(-0.157565e-2) + y * p4;
      39          20 :                     const T p6 = static_cast<T>(0.225319e-2) + y * p5;
      40          20 :                     const T p7 = static_cast<T>(0.1328592e-1) + y * p6;
      41          20 :                     const T p8 = static_cast<T>(0.39894228) + y * p7;
      42             : 
      43          20 :                     return arg - std::sqrt(arg) + std::log(p8);
      44          20 :                 }
      45         189 :             }
      46             :         };
      47             : 
      48             :         struct BesselFn_1_0 {
      49             :             template <typename T>
      50             :             __host__ __device__ constexpr T operator()(const T& arg) const noexcept
      51         346 :             {
      52         346 :                 if (arg < static_cast<T>(3.75)) {
      53         288 :                     T bessel_0_p6 = 0;
      54         288 :                     const T ratio = arg / static_cast<T>(3.75);
      55         288 :                     const T y = ratio * ratio;
      56             : 
      57         288 :                     {
      58         288 :                         const T p0 = static_cast<T>(0.45813e-2);
      59         288 :                         const T p1 = static_cast<T>(0.360768e-1) + y * p0;
      60         288 :                         const T p2 = static_cast<T>(0.2659732) + y * p1;
      61         288 :                         const T p3 = static_cast<T>(1.2067492) + y * p2;
      62         288 :                         const T p4 = static_cast<T>(3.0899424) + y * p3;
      63         288 :                         const T p5 = static_cast<T>(3.5156229) + y * p4;
      64         288 :                         bessel_0_p6 = static_cast<T>(1.0) + y * p5;
      65         288 :                     }
      66             : 
      67         288 :                     T bessel_1_p6 = 0;
      68         288 :                     {
      69         288 :                         const T p0 = static_cast<T>(0.32411e-3);
      70         288 :                         const T p1 = static_cast<T>(0.301532e-2) + y * p0;
      71         288 :                         const T p2 = static_cast<T>(0.2658733e-1) + y * p1;
      72         288 :                         const T p3 = static_cast<T>(0.15084934) + y * p2;
      73         288 :                         const T p4 = static_cast<T>(0.51498869) + y * p3;
      74         288 :                         const T p5 = static_cast<T>(0.87890594) + y * p4;
      75         288 :                         bessel_1_p6 = arg * (static_cast<T>(0.5) + y * p5);
      76         288 :                     }
      77             : 
      78         288 :                     return bessel_1_p6 / bessel_0_p6;
      79         288 :                 } 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          58 :                     const T y = static_cast<T>(3.75) / arg;
      84             : 
      85          58 :                     T bessel_0_p8;
      86          58 :                     {
      87          58 :                         const T p0 = static_cast<T>(0.392377e-2);
      88          58 :                         const T p1 = static_cast<T>(-0.1647633e-1) + y * p0;
      89          58 :                         const T p2 = static_cast<T>(0.2635537e-1) + y * p1;
      90          58 :                         const T p3 = static_cast<T>(-0.2057706e-1) + y * p2;
      91          58 :                         const T p4 = static_cast<T>(0.916281e-2) + y * p3;
      92          58 :                         const T p5 = static_cast<T>(-0.157565e-2) + y * p4;
      93          58 :                         const T p6 = static_cast<T>(0.225319e-2) + y * p5;
      94          58 :                         const T p7 = static_cast<T>(0.1328592e-1) + y * p6;
      95          58 :                         bessel_0_p8 = static_cast<T>(0.39894228) + y * p7;
      96          58 :                     }
      97             : 
      98          58 :                     T bessel_1_p8;
      99          58 :                     {
     100          58 :                         const T p0 = static_cast<T>(0.420059e-2);
     101          58 :                         const T p1 = static_cast<T>(0.1787654e-1) - y * p0;
     102          58 :                         const T p2 = static_cast<T>(-0.2895312e-1) + y * p1;
     103          58 :                         const T p3 = static_cast<T>(0.2282967e-1) + y * p2;
     104          58 :                         const T p4 = static_cast<T>(-0.1031555e-1) + y * p3;
     105          58 :                         const T p5 = static_cast<T>(0.163801e-2) + y * p4;
     106          58 :                         const T p6 = static_cast<T>(-0.362018e-2) + y * p5;
     107          58 :                         const T p7 = static_cast<T>(-0.3988024e-1) + y * p6;
     108          58 :                         bessel_1_p8 = static_cast<T>(0.39894228) + y * p7;
     109          58 :                     }
     110             : 
     111          58 :                     return bessel_1_p8 / bessel_0_p8;
     112          58 :                 }
     113         346 :             }
     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