1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
|
// boost sinhc.hpp header file
// (C) Copyright Hubert Holin 2001.
// Distributed under the Boost Software License, Version 1.0. (See
// accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
// See http://www.boost.org for updates, documentation, and revision history.
#ifndef BOOST_SINHC_HPP
#define BOOST_SINHC_HPP
#ifdef _MSC_VER
#pragma once
#endif
#include <boost/math/tools/precision.hpp>
#include <boost/math/policies/error_handling.hpp>
#include <boost/math/special_functions/math_fwd.hpp>
#include <boost/math/special_functions/fpclassify.hpp>
#include <limits>
#include <string>
#include <stdexcept>
#include <cmath>
// These are the the "Hyperbolic Sinus Cardinal" functions.
namespace boost
{
namespace math
{
namespace detail
{
// This is the "Hyperbolic Sinus Cardinal" of index Pi.
template<typename T, typename Policy>
inline T sinhc_pi_imp(const T x, const Policy&)
{
using ::std::abs;
using ::std::sinh;
using ::std::sqrt;
static T const taylor_0_bound = tools::epsilon<T>();
static T const taylor_2_bound = sqrt(taylor_0_bound);
static T const taylor_n_bound = sqrt(taylor_2_bound);
if((boost::math::isinf)(x))
{
return policies::raise_overflow_error<T>("sinhc(%1%)", nullptr, Policy());
}
if (abs(x) >= taylor_n_bound)
{
return(sinh(x)/x);
}
else
{
// approximation by taylor series in x at 0 up to order 0
T result = static_cast<T>(1);
if (abs(x) >= taylor_0_bound)
{
T x2 = x*x;
// approximation by taylor series in x at 0 up to order 2
result += x2/static_cast<T>(6);
if (abs(x) >= taylor_2_bound)
{
// approximation by taylor series in x at 0 up to order 4
result += (x2*x2)/static_cast<T>(120);
}
}
return(result);
}
}
} // namespace detail
template <class T, class Policy>
inline typename tools::promote_args<T>::type sinhc_pi(T x, const Policy& pol)
{
typedef typename tools::promote_args<T>::type result_type;
return policies::checked_narrowing_cast<T, Policy>(detail::sinhc_pi_imp(static_cast<result_type>(x), pol), "sinhc(%1%)");
}
template <class T>
inline typename tools::promote_args<T>::type sinhc_pi(T x)
{
typedef typename tools::promote_args<T>::type result_type;
return sinhc_pi(static_cast<result_type>(x), policies::policy<>());
}
template<typename T, template<typename> class U>
inline U<T> sinhc_pi(const U<T> x)
{
using std::abs;
using std::sinh;
using std::sqrt;
using ::std::numeric_limits;
static T const taylor_0_bound = tools::epsilon<T>();
static T const taylor_2_bound = sqrt(taylor_0_bound);
static T const taylor_n_bound = sqrt(taylor_2_bound);
if (abs(x) >= taylor_n_bound)
{
return(sinh(x)/x);
}
else
{
// approximation by taylor series in x at 0 up to order 0
#ifdef __MWERKS__
U<T> result = static_cast<U<T> >(1);
#else
U<T> result = U<T>(1);
#endif
if (abs(x) >= taylor_0_bound)
{
U<T> x2 = x*x;
// approximation by taylor series in x at 0 up to order 2
result += x2/static_cast<T>(6);
if (abs(x) >= taylor_2_bound)
{
// approximation by taylor series in x at 0 up to order 4
result += (x2*x2)/static_cast<T>(120);
}
}
return(result);
}
}
}
}
#endif /* BOOST_SINHC_HPP */
|
