184 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_Y1_HPP7#define BOOST_MATH_BESSEL_Y1_HPP8 9#ifdef _MSC_VER10#pragma once11#pragma warning(push)12#pragma warning(disable:4702) // Unreachable code (release mode only warning)13#endif14 15#include <boost/math/tools/config.hpp>16#include <boost/math/special_functions/detail/bessel_j1.hpp>17#include <boost/math/constants/constants.hpp>18#include <boost/math/tools/rational.hpp>19#include <boost/math/tools/big_constant.hpp>20#include <boost/math/policies/error_handling.hpp>21#include <boost/math/tools/assert.hpp>22 23#if defined(__GNUC__) && defined(BOOST_MATH_USE_FLOAT128)24//25// This is the only way we can avoid26// warning: non-standard suffix on floating constant [-Wpedantic]27// when building with -Wall -pedantic. Neither __extension__28// nor #pragma diagnostic ignored work :(29//30#pragma GCC system_header31#endif32 33// Bessel function of the second kind of order one34// x <= 8, minimax rational approximations on root-bracketing intervals35// x > 8, Hankel asymptotic expansion in Hart, Computer Approximations, 196836 37namespace boost { namespace math { namespace detail{38 39template <typename T, typename Policy>40BOOST_MATH_GPU_ENABLED T bessel_y1(T x, const Policy&);41 42template <typename T, typename Policy>43BOOST_MATH_GPU_ENABLED T bessel_y1(T x, const Policy&)44{45 BOOST_MATH_STATIC const T P1[] = {46 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.0535726612579544093e+13)),47 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.4708611716525426053e+12)),48 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.7595974497819597599e+11)),49 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 7.2144548214502560419e+09)),50 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -5.9157479997408395984e+07)),51 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.2157953222280260820e+05)),52 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.1714424660046133456e+02)),53 };54 BOOST_MATH_STATIC const T Q1[] = {55 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.0737873921079286084e+14)),56 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.1272286200406461981e+12)),57 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.7800352738690585613e+10)),58 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.2250435122182963220e+08)),59 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.8136470753052572164e+05)),60 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.2079908168393867438e+02)),61 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)),62 };63 BOOST_MATH_STATIC const T P2[] = {64 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.1514276357909013326e+19)),65 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -5.6808094574724204577e+18)),66 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.3638408497043134724e+16)),67 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.0686275289804744814e+15)),68 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -5.9530713129741981618e+13)),69 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.7453673962438488783e+11)),70 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.1957961912070617006e+09)),71 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.9153806858264202986e+06)),72 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.2337180442012953128e+03)),73 };74 BOOST_MATH_STATIC const T Q2[] = {75 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.3321844313316185697e+20)),76 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.6968198822857178911e+18)),77 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.0837179548112881950e+16)),78 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.1187010065856971027e+14)),79 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.0221766852960403645e+11)),80 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 6.3550318087088919566e+08)),81 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0453748201934079734e+06)),82 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.2855164849321609336e+03)),83 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)),84 };85 BOOST_MATH_STATIC const T PC[] = {86 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -4.4357578167941278571e+06)),87 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -9.9422465050776411957e+06)),88 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -6.6033732483649391093e+06)),89 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.5235293511811373833e+06)),90 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.0982405543459346727e+05)),91 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.6116166443246101165e+03)),92 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.0)),93 };94 BOOST_MATH_STATIC const T QC[] = {95 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -4.4357578167941278568e+06)),96 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -9.9341243899345856590e+06)),97 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -6.5853394797230870728e+06)),98 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.5118095066341608816e+06)),99 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.0726385991103820119e+05)),100 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.4550094401904961825e+03)),101 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)),102 };103 BOOST_MATH_STATIC const T PS[] = {104 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.3220913409857223519e+04)),105 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.5145160675335701966e+04)),106 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 6.6178836581270835179e+04)),107 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.8494262873223866797e+04)),108 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.7063754290207680021e+03)),109 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.5265133846636032186e+01)),110 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.0)),111 };112 BOOST_MATH_STATIC const T QS[] = {113 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 7.0871281941028743574e+05)),114 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.8194580422439972989e+06)),115 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.4194606696037208929e+06)),116 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.0029443582266975117e+05)),117 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.7890229745772202641e+04)),118 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.6383677696049909675e+02)),119 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)),120 };121 BOOST_MATH_STATIC const T x1 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.1971413260310170351e+00)),122 x2 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.4296810407941351328e+00)),123 x11 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.620e+02)),124 x12 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.8288260310170351490e-03)),125 x21 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.3900e+03)),126 x22 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -6.4592058648672279948e-06))127 ;128 T value, factor, r, rc, rs;129 130 BOOST_MATH_STD_USING131 using namespace boost::math::tools;132 using namespace boost::math::constants;133 134 BOOST_MATH_ASSERT(x > 0);135 136 if (x <= 4) // x in (0, 4]137 {138 T y = x * x;139 T z = 2 * log(x/x1) * bessel_j1(x) / pi<T>();140 r = evaluate_rational(P1, Q1, y);141 factor = (x + x1) * ((x - x11/256) - x12) / x;142 value = z + factor * r;143 }144 else if (x <= 8) // x in (4, 8]145 {146 T y = x * x;147 T z = 2 * log(x/x2) * bessel_j1(x) / pi<T>();148 r = evaluate_rational(P2, Q2, y);149 factor = (x + x2) * ((x - x21/256) - x22) / x;150 value = z + factor * r;151 }152 else // x in (8, \infty)153 {154 T y = 8 / x;155 T y2 = y * y;156 rc = evaluate_rational(PC, QC, y2);157 rs = evaluate_rational(PS, QS, y2);158 factor = 1 / (sqrt(x) * root_pi<T>());159 //160 // This code is really just:161 //162 // T z = x - 0.75f * pi<T>();163 // value = factor * (rc * sin(z) + y * rs * cos(z));164 //165 // But using the sin/cos addition rules, plus constants for sin/cos of 3PI/4166 // which then cancel out with corresponding terms in "factor".167 //168 T sx = sin(x);169 T cx = cos(x);170 value = factor * (y * rs * (sx - cx) - rc * (sx + cx));171 }172 173 return value;174}175 176}}} // namespaces177 178#ifdef _MSC_VER179#pragma warning(pop)180#endif181 182#endif // BOOST_MATH_BESSEL_Y1_HPP183 184