193 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_J0_HPP7#define BOOST_MATH_BESSEL_J0_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 zero30// 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_j0(T x);37 38template <typename T>39BOOST_MATH_GPU_ENABLED T bessel_j0(T x)40{41#ifdef BOOST_MATH_INSTRUMENT42 static bool b = false;43 if (!b)44 {45 std::cout << "bessel_j0 called with " << typeid(x).name() << std::endl;46 std::cout << "double = " << typeid(double).name() << std::endl;47 std::cout << "long double = " << typeid(long double).name() << std::endl;48 b = true;49 }50#endif51 52 BOOST_MATH_STATIC const T P1[] = {53 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -4.1298668500990866786e+11)),54 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.7282507878605942706e+10)),55 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -6.2140700423540120665e+08)),56 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 6.6302997904833794242e+06)),57 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.6629814655107086448e+04)),58 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0344222815443188943e+02)),59 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.2117036164593528341e-01))60 };61 BOOST_MATH_STATIC const T Q1[] = {62 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.3883787996332290397e+12)),63 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.6328198300859648632e+10)),64 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.3985097372263433271e+08)),65 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.5612696224219938200e+05)),66 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 9.3614022392337710626e+02)),67 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)),68 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.0))69 };70 BOOST_MATH_STATIC const T P2[] = {71 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.8319397969392084011e+03)),72 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.2254078161378989535e+04)),73 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -7.2879702464464618998e+03)),74 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0341910641583726701e+04)),75 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.1725046279757103576e+04)),76 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.4176707025325087628e+03)),77 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 7.4321196680624245801e+02)),78 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.8591703355916499363e+01))79 };80 BOOST_MATH_STATIC const T Q2[] = {81 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.5783478026152301072e+05)),82 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.4599102262586308984e+05)),83 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -8.4055062591169562211e+04)),84 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.8680990008359188352e+04)),85 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.9458766545509337327e+03)),86 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.3307310774649071172e+02)),87 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.5258076240801555057e+01)),88 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0))89 };90 BOOST_MATH_STATIC const T PC[] = {91 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.2779090197304684302e+04)),92 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.1345386639580765797e+04)),93 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.1170523380864944322e+04)),94 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.4806486443249270347e+03)),95 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.5376201909008354296e+02)),96 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.8961548424210455236e-01))97 };98 BOOST_MATH_STATIC const T QC[] = {99 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.2779090197304684318e+04)),100 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.1370412495510416640e+04)),101 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.1215350561880115730e+04)),102 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.5028735138235608207e+03)),103 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.5711159858080893649e+02)),104 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0))105 };106 BOOST_MATH_STATIC const T PS[] = {107 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -8.9226600200800094098e+01)),108 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.8591953644342993800e+02)),109 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.1183429920482737611e+02)),110 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.2300261666214198472e+01)),111 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.2441026745835638459e+00)),112 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -8.8033303048680751817e-03))113 };114 BOOST_MATH_STATIC const T QS[] = {115 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.7105024128512061905e+03)),116 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.1951131543434613647e+04)),117 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 7.2642780169211018836e+03)),118 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.4887231232283756582e+03)),119 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 9.0593769594993125859e+01)),120 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0))121 };122 123 BOOST_MATH_STATIC const T x1 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.4048255576957727686e+00));124 BOOST_MATH_STATIC const T x2 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.5200781102863106496e+00));125 BOOST_MATH_STATIC const T x11 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 6.160e+02));126 BOOST_MATH_STATIC const T x12 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.42444230422723137837e-03));127 BOOST_MATH_STATIC const T x21 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.4130e+03));128 BOOST_MATH_STATIC const T x22 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.46860286310649596604e-04));129 130 T value, factor, r, rc, rs;131 132 BOOST_MATH_STD_USING133 using namespace boost::math::tools;134 using namespace boost::math::constants;135 136 BOOST_MATH_ASSERT(x >= 0); // reflection handled elsewhere.137 138 if (x == 0)139 {140 return static_cast<T>(1);141 }142 if (x <= 4) // x in (0, 4]143 {144 T y = x * x;145 BOOST_MATH_ASSERT(sizeof(P1) == sizeof(Q1));146 r = evaluate_rational(P1, Q1, y);147 factor = (x + x1) * ((x - x11/256) - x12);148 value = factor * r;149 }150 else if (x <= 8.0) // x in (4, 8]151 {152 T y = 1 - (x * x)/64;153 BOOST_MATH_ASSERT(sizeof(P2) == sizeof(Q2));154 r = evaluate_rational(P2, Q2, y);155 factor = (x + x2) * ((x - x21/256) - x22);156 value = factor * r;157 }158 else // x in (8, \infty)159 {160 T y = 8 / x;161 T y2 = y * y;162 BOOST_MATH_ASSERT(sizeof(PC) == sizeof(QC));163 BOOST_MATH_ASSERT(sizeof(PS) == sizeof(QS));164 rc = evaluate_rational(PC, QC, y2);165 rs = evaluate_rational(PS, QS, y2);166 factor = constants::one_div_root_pi<T>() / sqrt(x);167 //168 // What follows is really just:169 //170 // T z = x - pi/4;171 // value = factor * (rc * cos(z) - y * rs * sin(z));172 //173 // But using the addition formulae for sin and cos, plus174 // the special values for sin/cos of pi/4.175 //176 T sx = sin(x);177 T cx = cos(x);178 BOOST_MATH_INSTRUMENT_VARIABLE(rc);179 BOOST_MATH_INSTRUMENT_VARIABLE(rs);180 BOOST_MATH_INSTRUMENT_VARIABLE(factor);181 BOOST_MATH_INSTRUMENT_VARIABLE(sx);182 BOOST_MATH_INSTRUMENT_VARIABLE(cx);183 value = factor * (rc * (cx + sx) - y * rs * (sx - cx));184 }185 186 return value;187}188 189}}} // namespaces190 191#endif // BOOST_MATH_BESSEL_J0_HPP192 193