180 lines · plain
1// Copyright (c) 2006 Xiaogang Zhang2// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0. (See accompanying file4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)5 6#ifndef BOOST_MATH_BESSEL_J1_HPP7#define BOOST_MATH_BESSEL_J1_HPP8 9#ifdef _MSC_VER10#pragma once11#endif12 13#include <boost/math/tools/config.hpp>14#include <boost/math/constants/constants.hpp>15#include <boost/math/tools/rational.hpp>16#include <boost/math/tools/big_constant.hpp>17#include <boost/math/tools/assert.hpp>18 19#if defined(__GNUC__) && defined(BOOST_MATH_USE_FLOAT128)20//21// This is the only way we can avoid22// warning: non-standard suffix on floating constant [-Wpedantic]23// when building with -Wall -pedantic. Neither __extension__24// nor #pragma diagnostic ignored work :(25//26#pragma GCC system_header27#endif28 29// Bessel function of the first kind of order one30// x <= 8, minimax rational approximations on root-bracketing intervals31// x > 8, Hankel asymptotic expansion in Hart, Computer Approximations, 196832 33namespace boost { namespace math{ namespace detail{34 35template <typename T>36BOOST_MATH_GPU_ENABLED T bessel_j1(T x)37{38 BOOST_MATH_STATIC const T P1[] = {39 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.4258509801366645672e+11)),40 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 6.6781041261492395835e+09)),41 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.1548696764841276794e+08)),42 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 9.8062904098958257677e+05)),43 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -4.4615792982775076130e+03)),44 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0650724020080236441e+01)),45 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.0767857011487300348e-02))46 };47 BOOST_MATH_STATIC const T Q1[] = {48 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.1868604460820175290e+12)),49 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.2091902282580133541e+10)),50 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.0228375140097033958e+08)),51 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.9117614494174794095e+05)),52 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0742272239517380498e+03)),53 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)),54 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.0))55 };56 BOOST_MATH_STATIC const T P2[] = {57 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.7527881995806511112e+16)),58 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.6608531731299018674e+15)),59 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.6658018905416665164e+13)),60 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.5580665670910619166e+11)),61 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.8113931269860667829e+09)),62 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.0793266148011179143e+06)),63 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -7.5023342220781607561e+03)),64 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.6179191852758252278e+00))65 };66 BOOST_MATH_STATIC const T Q2[] = {67 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.7253905888447681194e+18)),68 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.7128800897135812012e+16)),69 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.4899346165481429307e+13)),70 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.7622777286244082666e+11)),71 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 6.4872502899596389593e+08)),72 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.1267125065029138050e+06)),73 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.3886978985861357615e+03)),74 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0))75 };76 BOOST_MATH_STATIC const T PC[] = {77 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -4.4357578167941278571e+06)),78 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -9.9422465050776411957e+06)),79 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -6.6033732483649391093e+06)),80 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.5235293511811373833e+06)),81 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.0982405543459346727e+05)),82 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.6116166443246101165e+03)),83 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.0))84 };85 BOOST_MATH_STATIC const T QC[] = {86 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -4.4357578167941278568e+06)),87 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -9.9341243899345856590e+06)),88 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -6.5853394797230870728e+06)),89 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.5118095066341608816e+06)),90 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.0726385991103820119e+05)),91 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.4550094401904961825e+03)),92 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0))93 };94 BOOST_MATH_STATIC const T PS[] = {95 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.3220913409857223519e+04)),96 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.5145160675335701966e+04)),97 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 6.6178836581270835179e+04)),98 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.8494262873223866797e+04)),99 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.7063754290207680021e+03)),100 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.5265133846636032186e+01)),101 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.0))102 };103 BOOST_MATH_STATIC const T QS[] = {104 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 7.0871281941028743574e+05)),105 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.8194580422439972989e+06)),106 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.4194606696037208929e+06)),107 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.0029443582266975117e+05)),108 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.7890229745772202641e+04)),109 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.6383677696049909675e+02)),110 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0))111 };112 113 BOOST_MATH_STATIC const T x1 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.8317059702075123156e+00));114 BOOST_MATH_STATIC const T x2 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 7.0155866698156187535e+00));115 BOOST_MATH_STATIC const T x11 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 9.810e+02));116 BOOST_MATH_STATIC const T x12 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.2527979248768438556e-04));117 BOOST_MATH_STATIC const T x21 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.7960e+03));118 BOOST_MATH_STATIC const T x22 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.8330184381246462950e-05));119 120 T value, factor, r, rc, rs, w;121 122 BOOST_MATH_STD_USING123 using namespace boost::math::tools;124 using namespace boost::math::constants;125 126 w = abs(x);127 if (x == 0)128 {129 return static_cast<T>(0);130 }131 if (w <= 4) // w in (0, 4]132 {133 T y = x * x;134 BOOST_MATH_ASSERT(sizeof(P1) == sizeof(Q1));135 r = evaluate_rational(P1, Q1, y);136 factor = w * (w + x1) * ((w - x11/256) - x12);137 value = factor * r;138 }139 else if (w <= 8) // w in (4, 8]140 {141 T y = x * x;142 BOOST_MATH_ASSERT(sizeof(P2) == sizeof(Q2));143 r = evaluate_rational(P2, Q2, y);144 factor = w * (w + x2) * ((w - x21/256) - x22);145 value = factor * r;146 }147 else // w in (8, \infty)148 {149 T y = 8 / w;150 T y2 = y * y;151 BOOST_MATH_ASSERT(sizeof(PC) == sizeof(QC));152 BOOST_MATH_ASSERT(sizeof(PS) == sizeof(QS));153 rc = evaluate_rational(PC, QC, y2);154 rs = evaluate_rational(PS, QS, y2);155 factor = 1 / (sqrt(w) * constants::root_pi<T>());156 //157 // What follows is really just:158 //159 // T z = w - 0.75f * pi<T>();160 // value = factor * (rc * cos(z) - y * rs * sin(z));161 //162 // but using the sin/cos addition rules plus constants163 // for the values of sin/cos of 3PI/4 which then cancel164 // out with corresponding terms in "factor".165 //166 T sx = sin(x);167 T cx = cos(x);168 value = factor * (rc * (sx - cx) + y * rs * (sx + cx));169 }170 171 BOOST_MATH_ASSERT(x >= 0); // Negative values handled by the caller.172 173 return value;174}175 176}}} // namespaces177 178#endif // BOOST_MATH_BESSEL_J1_HPP179 180