boost/math/special_functions/detail/bessel_j1.hpp

   1 //  Copyright (c) 2006 Xiaogang Zhang
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5
   6 #ifndef BOOST_MATH_BESSEL_J1_HPP
   7 #define BOOST_MATH_BESSEL_J1_HPP
   8
   9 #ifdef _MSC_VER
  10 #pragma once
  11 #endif
  12
  13 #include <boost/math/constants/constants.hpp>
  14 #include <boost/math/tools/rational.hpp>
  15 #include <boost/math/tools/big_constant.hpp>
  16 #include <boost/assert.hpp>
  17
  18 // Bessel function of the first kind of order one
  19 // x <= 8, minimax rational approximations on root-bracketing intervals
  20 // x > 8, Hankel asymptotic expansion in Hart, Computer Approximations, 1968
  21
  22 namespace boost { namespace math{  namespace detail{
  23
  24 template <typename T>
  25 T bessel_j1(T x);
  26
  27 template <class T>
  28 struct bessel_j1_initializer
  29 {
  30    struct init
  31    {
  32       init()
  33       {
  34          do_init();
  35       }
  36       static void do_init()
  37       {
  38          bessel_j1(T(1));
  39       }
  40       void force_instantiate()const{}
  41    };
  42    static const init initializer;
  43    static void force_instantiate()
  44    {
  45       initializer.force_instantiate();
  46    }
  47 };
  48
  49 template <class T>
  50 const typename bessel_j1_initializer<T>::init bessel_j1_initializer<T>::initializer;
  51
  52 template <typename T>
  53 T bessel_j1(T x)
  54 {
  55     bessel_j1_initializer<T>::force_instantiate();
  56
  57     static const T P1[] = {
  58          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.4258509801366645672e+11)),
  59          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 6.6781041261492395835e+09)),
  60          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.1548696764841276794e+08)),
  61          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 9.8062904098958257677e+05)),
  62          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -4.4615792982775076130e+03)),
  63          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0650724020080236441e+01)),
  64          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.0767857011487300348e-02))
  65     };
  66     static const T Q1[] = {
  67          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.1868604460820175290e+12)),
  68          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.2091902282580133541e+10)),
  69          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.0228375140097033958e+08)),
  70          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.9117614494174794095e+05)),
  71          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0742272239517380498e+03)),
  72          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)),
  73          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.0))
  74     };
  75     static const T P2[] = {
  76          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.7527881995806511112e+16)),
  77          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.6608531731299018674e+15)),
  78          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.6658018905416665164e+13)),
  79          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.5580665670910619166e+11)),
  80          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.8113931269860667829e+09)),
  81          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.0793266148011179143e+06)),
  82          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -7.5023342220781607561e+03)),
  83          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.6179191852758252278e+00))
  84     };
  85     static const T Q2[] = {
  86          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.7253905888447681194e+18)),
  87          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.7128800897135812012e+16)),
  88          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.4899346165481429307e+13)),
  89          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.7622777286244082666e+11)),
  90          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 6.4872502899596389593e+08)),
  91          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.1267125065029138050e+06)),
  92          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.3886978985861357615e+03)),
  93          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0))
  94     };
  95     static const T PC[] = {
  96         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -4.4357578167941278571e+06)),
  97         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -9.9422465050776411957e+06)),
  98         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -6.6033732483649391093e+06)),
  99         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.5235293511811373833e+06)),
 100         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.0982405543459346727e+05)),
 101         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.6116166443246101165e+03)),
 102         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.0))
 103     };
 104     static const T QC[] = {
 105         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -4.4357578167941278568e+06)),
 106         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -9.9341243899345856590e+06)),
 107         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -6.5853394797230870728e+06)),
 108         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.5118095066341608816e+06)),
 109         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.0726385991103820119e+05)),
 110         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.4550094401904961825e+03)),
 111         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0))
 112     };
 113     static const T PS[] = {
 114          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.3220913409857223519e+04)),
 115          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.5145160675335701966e+04)),
 116          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 6.6178836581270835179e+04)),
 117          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.8494262873223866797e+04)),
 118          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.7063754290207680021e+03)),
 119          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.5265133846636032186e+01)),
 120          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.0))
 121     };
 122     static const T QS[] = {
 123          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 7.0871281941028743574e+05)),
 124          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.8194580422439972989e+06)),
 125          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.4194606696037208929e+06)),
 126          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.0029443582266975117e+05)),
 127          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.7890229745772202641e+04)),
 128          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.6383677696049909675e+02)),
 129          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0))
 130     };
 131     static const T x1  =  static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.8317059702075123156e+00)),
 132                    x2  =  static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 7.0155866698156187535e+00)),
 133                    x11 =  static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 9.810e+02)),
 134                    x12 =  static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.2527979248768438556e-04)),
 135                    x21 =  static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.7960e+03)),
 136                    x22 =  static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.8330184381246462950e-05));
 137
 138     T value, factor, r, rc, rs, w;
 139
 140     BOOST_MATH_STD_USING
 141     using namespace boost::math::tools;
 142     using namespace boost::math::constants;
 143
 144     w = abs(x);
 145     if (x == 0)
 146     {
 147         return static_cast<T>(0);
 148     }
 149     if (w <= 4)                       // w in (0, 4]
 150     {
 151         T y = x * x;
 152         BOOST_ASSERT(sizeof(P1) == sizeof(Q1));
 153         r = evaluate_rational(P1, Q1, y);
 154         factor = w * (w + x1) * ((w - x11/256) - x12);
 155         value = factor * r;
 156     }
 157     else if (w <= 8)                  // w in (4, 8]
 158     {
 159         T y = x * x;
 160         BOOST_ASSERT(sizeof(P2) == sizeof(Q2));
 161         r = evaluate_rational(P2, Q2, y);
 162         factor = w * (w + x2) * ((w - x21/256) - x22);
 163         value = factor * r;
 164     }
 165     else                                // w in (8, \infty)
 166     {
 167         T y = 8 / w;
 168         T y2 = y * y;
 169         BOOST_ASSERT(sizeof(PC) == sizeof(QC));
 170         BOOST_ASSERT(sizeof(PS) == sizeof(QS));
 171         rc = evaluate_rational(PC, QC, y2);
 172         rs = evaluate_rational(PS, QS, y2);
 173         factor = 1 / (sqrt(w) * constants::root_pi<T>());
 174         //
 175         // What follows is really just:
 176         //
 177         // T z = w - 0.75f * pi<T>();
 178         // value = factor * (rc * cos(z) - y * rs * sin(z));
 179         //
 180         // but using the sin/cos addition rules plus constants
 181         // for the values of sin/cos of 3PI/4 which then cancel
 182         // out with corresponding terms in "factor".
 183         //
 184         T sx = sin(x);
 185         T cx = cos(x);
 186         value = factor * (rc * (sx - cx) + y * rs * (sx + cx));
 187     }
 188
 189     if (x < 0)
 190     {
 191         value *= -1;                 // odd function
 192     }
 193     return value;
 194 }
 195
 196 }}} // namespaces
 197
 198 #endif // BOOST_MATH_BESSEL_J1_HPP
 199