special_functions/detail/bessel_j1.hpp

0001 //  Copyright (c) 2006 Xiaogang Zhang
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_MATH_BESSEL_J1_HPP
0007 #define BOOST_MATH_BESSEL_J1_HPP
0008
0009 #ifdef _MSC_VER
0010 #pragma once
0011 #endif
0012
0013 #include <boost/math/tools/config.hpp>
0014 #include <boost/math/constants/constants.hpp>
0015 #include <boost/math/tools/rational.hpp>
0016 #include <boost/math/tools/big_constant.hpp>
0017 #include <boost/math/tools/assert.hpp>
0018
0019 #if defined(__GNUC__) && defined(BOOST_MATH_USE_FLOAT128)
0020 //
0021 // This is the only way we can avoid
0022 // warning: non-standard suffix on floating constant [-Wpedantic]
0023 // when building with -Wall -pedantic.  Neither __extension__
0024 // nor #pragma diagnostic ignored work :(
0025 //
0026 #pragma GCC system_header
0027 #endif
0028
0029 // Bessel function of the first kind of order one
0030 // x <= 8, minimax rational approximations on root-bracketing intervals
0031 // x > 8, Hankel asymptotic expansion in Hart, Computer Approximations, 1968
0032
0033 namespace boost { namespace math{  namespace detail{
0034
0035 template <typename T>
0036 BOOST_MATH_GPU_ENABLED T bessel_j1(T x);
0037
0038 template <class T>
0039 struct bessel_j1_initializer
0040 {
0041    struct init
0042    {
0043       BOOST_MATH_GPU_ENABLED init()
0044       {
0045          do_init();
0046       }
0047       BOOST_MATH_GPU_ENABLED static void do_init()
0048       {
0049          bessel_j1(T(1));
0050       }
0051       BOOST_MATH_GPU_ENABLED void force_instantiate()const{}
0052    };
0053    BOOST_MATH_STATIC const init initializer;
0054    BOOST_MATH_GPU_ENABLED static void force_instantiate()
0055    {
0056       #ifndef BOOST_MATH_HAS_GPU_SUPPORT
0057       initializer.force_instantiate();
0058       #endif
0059    }
0060 };
0061
0062 template <class T>
0063 const typename bessel_j1_initializer<T>::init bessel_j1_initializer<T>::initializer;
0064
0065 template <typename T>
0066 BOOST_MATH_GPU_ENABLED T bessel_j1(T x)
0067 {
0068     bessel_j1_initializer<T>::force_instantiate();
0069
0070     BOOST_MATH_STATIC const T P1[] = {
0071          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.4258509801366645672e+11)),
0072          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 6.6781041261492395835e+09)),
0073          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.1548696764841276794e+08)),
0074          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 9.8062904098958257677e+05)),
0075          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -4.4615792982775076130e+03)),
0076          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0650724020080236441e+01)),
0077          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.0767857011487300348e-02))
0078     };
0079     BOOST_MATH_STATIC const T Q1[] = {
0080          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.1868604460820175290e+12)),
0081          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.2091902282580133541e+10)),
0082          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.0228375140097033958e+08)),
0083          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.9117614494174794095e+05)),
0084          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0742272239517380498e+03)),
0085          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)),
0086          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.0))
0087     };
0088     BOOST_MATH_STATIC const T P2[] = {
0089          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.7527881995806511112e+16)),
0090          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.6608531731299018674e+15)),
0091          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.6658018905416665164e+13)),
0092          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.5580665670910619166e+11)),
0093          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.8113931269860667829e+09)),
0094          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.0793266148011179143e+06)),
0095          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -7.5023342220781607561e+03)),
0096          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.6179191852758252278e+00))
0097     };
0098     BOOST_MATH_STATIC const T Q2[] = {
0099          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.7253905888447681194e+18)),
0100          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.7128800897135812012e+16)),
0101          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.4899346165481429307e+13)),
0102          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.7622777286244082666e+11)),
0103          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 6.4872502899596389593e+08)),
0104          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.1267125065029138050e+06)),
0105          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.3886978985861357615e+03)),
0106          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0))
0107     };
0108     BOOST_MATH_STATIC const T PC[] = {
0109         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -4.4357578167941278571e+06)),
0110         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -9.9422465050776411957e+06)),
0111         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -6.6033732483649391093e+06)),
0112         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.5235293511811373833e+06)),
0113         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.0982405543459346727e+05)),
0114         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.6116166443246101165e+03)),
0115         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.0))
0116     };
0117     BOOST_MATH_STATIC const T QC[] = {
0118         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -4.4357578167941278568e+06)),
0119         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -9.9341243899345856590e+06)),
0120         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -6.5853394797230870728e+06)),
0121         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.5118095066341608816e+06)),
0122         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.0726385991103820119e+05)),
0123         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.4550094401904961825e+03)),
0124         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0))
0125     };
0126     BOOST_MATH_STATIC const T PS[] = {
0127          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.3220913409857223519e+04)),
0128          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.5145160675335701966e+04)),
0129          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 6.6178836581270835179e+04)),
0130          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.8494262873223866797e+04)),
0131          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.7063754290207680021e+03)),
0132          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.5265133846636032186e+01)),
0133          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.0))
0134     };
0135     BOOST_MATH_STATIC const T QS[] = {
0136          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 7.0871281941028743574e+05)),
0137          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.8194580422439972989e+06)),
0138          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.4194606696037208929e+06)),
0139          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.0029443582266975117e+05)),
0140          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.7890229745772202641e+04)),
0141          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.6383677696049909675e+02)),
0142          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0))
0143     };
0144
0145     BOOST_MATH_STATIC const T x1  =  static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.8317059702075123156e+00));
0146     BOOST_MATH_STATIC const T x2  =  static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 7.0155866698156187535e+00));
0147     BOOST_MATH_STATIC const T x11 =  static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 9.810e+02));
0148     BOOST_MATH_STATIC const T x12 =  static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.2527979248768438556e-04));
0149     BOOST_MATH_STATIC const T x21 =  static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.7960e+03));
0150     BOOST_MATH_STATIC const T x22 =  static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.8330184381246462950e-05));
0151
0152     T value, factor, r, rc, rs, w;
0153
0154     BOOST_MATH_STD_USING
0155     using namespace boost::math::tools;
0156     using namespace boost::math::constants;
0157
0158     w = abs(x);
0159     if (x == 0)
0160     {
0161         return static_cast<T>(0);
0162     }
0163     if (w <= 4)                       // w in (0, 4]
0164     {
0165         T y = x * x;
0166         BOOST_MATH_ASSERT(sizeof(P1) == sizeof(Q1));
0167         r = evaluate_rational(P1, Q1, y);
0168         factor = w * (w + x1) * ((w - x11/256) - x12);
0169         value = factor * r;
0170     }
0171     else if (w <= 8)                  // w in (4, 8]
0172     {
0173         T y = x * x;
0174         BOOST_MATH_ASSERT(sizeof(P2) == sizeof(Q2));
0175         r = evaluate_rational(P2, Q2, y);
0176         factor = w * (w + x2) * ((w - x21/256) - x22);
0177         value = factor * r;
0178     }
0179     else                                // w in (8, \infty)
0180     {
0181         T y = 8 / w;
0182         T y2 = y * y;
0183         BOOST_MATH_ASSERT(sizeof(PC) == sizeof(QC));
0184         BOOST_MATH_ASSERT(sizeof(PS) == sizeof(QS));
0185         rc = evaluate_rational(PC, QC, y2);
0186         rs = evaluate_rational(PS, QS, y2);
0187         factor = 1 / (sqrt(w) * constants::root_pi<T>());
0188         //
0189         // What follows is really just:
0190         //
0191         // T z = w - 0.75f * pi<T>();
0192         // value = factor * (rc * cos(z) - y * rs * sin(z));
0193         //
0194         // but using the sin/cos addition rules plus constants
0195         // for the values of sin/cos of 3PI/4 which then cancel
0196         // out with corresponding terms in "factor".
0197         //
0198         T sx = sin(x);
0199         T cx = cos(x);
0200         value = factor * (rc * (sx - cx) + y * rs * (sx + cx));
0201     }
0202
0203     BOOST_MATH_ASSERT(x >= 0);  // Negative values handled by the caller.
0204
0205     return value;
0206 }
0207
0208 }}} // namespaces
0209
0210 #endif // BOOST_MATH_BESSEL_J1_HPP
0211