Mercurial > hg > vamp-build-and-test

annotate DEPENDENCIES/generic/include/boost/math/special_functions/detail/bessel_y1.hpp @ 133:4acb5d8d80b6 tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Don't fail environmental check if README.md exists (but .txt and no-suffix don't)
author Chris Cannam
date Tue, 30 Jul 2019 12:25:44 +0100
parents 2665513ce2d3
children
rev   line source
Chris@16 1 // Copyright (c) 2006 Xiaogang Zhang
Chris@16 2 // Use, modification and distribution are subject to the
Chris@16 3 // Boost Software License, Version 1.0. (See accompanying file
Chris@16 4 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@16 5
Chris@16 6 #ifndef BOOST_MATH_BESSEL_Y1_HPP
Chris@16 7 #define BOOST_MATH_BESSEL_Y1_HPP
Chris@16 8
Chris@16 9 #ifdef _MSC_VER
Chris@16 10 #pragma once
Chris@16 11 #endif
Chris@16 12
Chris@16 13 #include <boost/math/special_functions/detail/bessel_j1.hpp>
Chris@16 14 #include <boost/math/constants/constants.hpp>
Chris@16 15 #include <boost/math/tools/rational.hpp>
Chris@16 16 #include <boost/math/tools/big_constant.hpp>
Chris@16 17 #include <boost/math/policies/error_handling.hpp>
Chris@16 18 #include <boost/assert.hpp>
Chris@16 19
Chris@16 20 // Bessel function of the second kind of order one
Chris@16 21 // x <= 8, minimax rational approximations on root-bracketing intervals
Chris@16 22 // x > 8, Hankel asymptotic expansion in Hart, Computer Approximations, 1968
Chris@16 23
Chris@16 24 namespace boost { namespace math { namespace detail{
Chris@16 25
Chris@16 26 template <typename T, typename Policy>
Chris@16 27 T bessel_y1(T x, const Policy&);
Chris@16 28
Chris@16 29 template <class T, class Policy>
Chris@16 30 struct bessel_y1_initializer
Chris@16 31 {
Chris@16 32 struct init
Chris@16 33 {
Chris@16 34 init()
Chris@16 35 {
Chris@16 36 do_init();
Chris@16 37 }
Chris@16 38 static void do_init()
Chris@16 39 {
Chris@16 40 bessel_y1(T(1), Policy());
Chris@16 41 }
Chris@16 42 void force_instantiate()const{}
Chris@16 43 };
Chris@16 44 static const init initializer;
Chris@16 45 static void force_instantiate()
Chris@16 46 {
Chris@16 47 initializer.force_instantiate();
Chris@16 48 }
Chris@16 49 };
Chris@16 50
Chris@16 51 template <class T, class Policy>
Chris@16 52 const typename bessel_y1_initializer<T, Policy>::init bessel_y1_initializer<T, Policy>::initializer;
Chris@16 53
Chris@16 54 template <typename T, typename Policy>
Chris@16 55 T bessel_y1(T x, const Policy& pol)
Chris@16 56 {
Chris@16 57 bessel_y1_initializer<T, Policy>::force_instantiate();
Chris@16 58
Chris@16 59 static const T P1[] = {
Chris@16 60 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.0535726612579544093e+13)),
Chris@16 61 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.4708611716525426053e+12)),
Chris@16 62 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.7595974497819597599e+11)),
Chris@16 63 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 7.2144548214502560419e+09)),
Chris@16 64 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -5.9157479997408395984e+07)),
Chris@16 65 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.2157953222280260820e+05)),
Chris@16 66 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.1714424660046133456e+02)),
Chris@16 67 };
Chris@16 68 static const T Q1[] = {
Chris@16 69 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.0737873921079286084e+14)),
Chris@16 70 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.1272286200406461981e+12)),
Chris@16 71 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.7800352738690585613e+10)),
Chris@16 72 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.2250435122182963220e+08)),
Chris@16 73 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.8136470753052572164e+05)),
Chris@16 74 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.2079908168393867438e+02)),
Chris@16 75 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)),
Chris@16 76 };
Chris@16 77 static const T P2[] = {
Chris@16 78 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.1514276357909013326e+19)),
Chris@16 79 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -5.6808094574724204577e+18)),
Chris@16 80 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.3638408497043134724e+16)),
Chris@16 81 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.0686275289804744814e+15)),
Chris@16 82 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -5.9530713129741981618e+13)),
Chris@16 83 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.7453673962438488783e+11)),
Chris@16 84 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.1957961912070617006e+09)),
Chris@16 85 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.9153806858264202986e+06)),
Chris@16 86 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.2337180442012953128e+03)),
Chris@16 87 };
Chris@16 88 static const T Q2[] = {
Chris@16 89 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.3321844313316185697e+20)),
Chris@16 90 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.6968198822857178911e+18)),
Chris@16 91 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.0837179548112881950e+16)),
Chris@16 92 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.1187010065856971027e+14)),
Chris@16 93 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.0221766852960403645e+11)),
Chris@16 94 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 6.3550318087088919566e+08)),
Chris@16 95 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0453748201934079734e+06)),
Chris@16 96 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.2855164849321609336e+03)),
Chris@16 97 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)),
Chris@16 98 };
Chris@16 99 static const T PC[] = {
Chris@16 100 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -4.4357578167941278571e+06)),
Chris@16 101 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -9.9422465050776411957e+06)),
Chris@16 102 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -6.6033732483649391093e+06)),
Chris@16 103 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.5235293511811373833e+06)),
Chris@16 104 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.0982405543459346727e+05)),
Chris@16 105 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.6116166443246101165e+03)),
Chris@16 106 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.0)),
Chris@16 107 };
Chris@16 108 static const T QC[] = {
Chris@16 109 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -4.4357578167941278568e+06)),
Chris@16 110 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -9.9341243899345856590e+06)),
Chris@16 111 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -6.5853394797230870728e+06)),
Chris@16 112 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.5118095066341608816e+06)),
Chris@16 113 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.0726385991103820119e+05)),
Chris@16 114 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.4550094401904961825e+03)),
Chris@16 115 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)),
Chris@16 116 };
Chris@16 117 static const T PS[] = {
Chris@16 118 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.3220913409857223519e+04)),
Chris@16 119 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.5145160675335701966e+04)),
Chris@16 120 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 6.6178836581270835179e+04)),
Chris@16 121 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.8494262873223866797e+04)),
Chris@16 122 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.7063754290207680021e+03)),
Chris@16 123 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.5265133846636032186e+01)),
Chris@16 124 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.0)),
Chris@16 125 };
Chris@16 126 static const T QS[] = {
Chris@16 127 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 7.0871281941028743574e+05)),
Chris@16 128 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.8194580422439972989e+06)),
Chris@16 129 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.4194606696037208929e+06)),
Chris@16 130 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.0029443582266975117e+05)),
Chris@16 131 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.7890229745772202641e+04)),
Chris@16 132 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.6383677696049909675e+02)),
Chris@16 133 static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)),
Chris@16 134 };
Chris@16 135 static const T x1 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.1971413260310170351e+00)),
Chris@16 136 x2 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.4296810407941351328e+00)),
Chris@16 137 x11 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.620e+02)),
Chris@16 138 x12 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.8288260310170351490e-03)),
Chris@16 139 x21 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.3900e+03)),
Chris@16 140 x22 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -6.4592058648672279948e-06))
Chris@16 141 ;
Chris@16 142 T value, factor, r, rc, rs;
Chris@16 143
Chris@16 144 BOOST_MATH_STD_USING
Chris@16 145 using namespace boost::math::tools;
Chris@16 146 using namespace boost::math::constants;
Chris@16 147
Chris@16 148 if (x <= 0)
Chris@16 149 {
Chris@16 150 return policies::raise_domain_error<T>("bost::math::bessel_y1<%1%>(%1%,%1%)",
Chris@16 151 "Got x == %1%, but x must be > 0, complex result not supported.", x, pol);
Chris@16 152 }
Chris@16 153 if (x <= 4) // x in (0, 4]
Chris@16 154 {
Chris@16 155 T y = x * x;
Chris@16 156 T z = 2 * log(x/x1) * bessel_j1(x) / pi<T>();
Chris@16 157 r = evaluate_rational(P1, Q1, y);
Chris@16 158 factor = (x + x1) * ((x - x11/256) - x12) / x;
Chris@16 159 value = z + factor * r;
Chris@16 160 }
Chris@16 161 else if (x <= 8) // x in (4, 8]
Chris@16 162 {
Chris@16 163 T y = x * x;
Chris@16 164 T z = 2 * log(x/x2) * bessel_j1(x) / pi<T>();
Chris@16 165 r = evaluate_rational(P2, Q2, y);
Chris@16 166 factor = (x + x2) * ((x - x21/256) - x22) / x;
Chris@16 167 value = z + factor * r;
Chris@16 168 }
Chris@16 169 else // x in (8, \infty)
Chris@16 170 {
Chris@16 171 T y = 8 / x;
Chris@16 172 T y2 = y * y;
Chris@16 173 rc = evaluate_rational(PC, QC, y2);
Chris@16 174 rs = evaluate_rational(PS, QS, y2);
Chris@16 175 factor = 1 / (sqrt(x) * root_pi<T>());
Chris@16 176 //
Chris@16 177 // This code is really just:
Chris@16 178 //
Chris@16 179 // T z = x - 0.75f * pi<T>();
Chris@16 180 // value = factor * (rc * sin(z) + y * rs * cos(z));
Chris@16 181 //
Chris@16 182 // But using the sin/cos addition rules, plus constants for sin/cos of 3PI/4
Chris@16 183 // which then cancel out with corresponding terms in "factor".
Chris@16 184 //
Chris@16 185 T sx = sin(x);
Chris@16 186 T cx = cos(x);
Chris@16 187 value = factor * (y * rs * (sx - cx) - rc * (sx + cx));
Chris@16 188 }
Chris@16 189
Chris@16 190 return value;
Chris@16 191 }
Chris@16 192
Chris@16 193 }}} // namespaces
Chris@16 194
Chris@16 195 #endif // BOOST_MATH_BESSEL_Y1_HPP
Chris@16 196