]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // (C) Copyright John Maddock 2006. |
2 | // Use, modification and distribution are subject to the | |
3 | // Boost Software License, Version 1.0. (See accompanying file | |
4 | // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) | |
5 | ||
6 | #ifndef BOOST_MATH_SF_TRIGAMMA_HPP | |
7 | #define BOOST_MATH_SF_TRIGAMMA_HPP | |
8 | ||
9 | #ifdef _MSC_VER | |
10 | #pragma once | |
11 | #endif | |
12 | ||
13 | #include <boost/math/special_functions/math_fwd.hpp> | |
14 | #include <boost/math/tools/rational.hpp> | |
15 | #include <boost/math/tools/series.hpp> | |
16 | #include <boost/math/tools/promotion.hpp> | |
17 | #include <boost/math/policies/error_handling.hpp> | |
18 | #include <boost/math/constants/constants.hpp> | |
19 | #include <boost/mpl/comparison.hpp> | |
20 | #include <boost/math/tools/big_constant.hpp> | |
21 | #include <boost/math/special_functions/polygamma.hpp> | |
22 | ||
23 | namespace boost{ | |
24 | namespace math{ | |
25 | namespace detail{ | |
26 | ||
27 | template<class T, class Policy> | |
28 | T polygamma_imp(const int n, T x, const Policy &pol); | |
29 | ||
30 | template <class T, class Policy> | |
31 | T trigamma_prec(T x, const mpl::int_<53>*, const Policy&) | |
32 | { | |
33 | // Max error in interpolated form: 3.736e-017 | |
34 | static const T offset = BOOST_MATH_BIG_CONSTANT(T, 53, 2.1093254089355469); | |
35 | static const T P_1_2[] = { | |
36 | BOOST_MATH_BIG_CONSTANT(T, 53, -1.1093280605946045), | |
37 | BOOST_MATH_BIG_CONSTANT(T, 53, -3.8310674472619321), | |
38 | BOOST_MATH_BIG_CONSTANT(T, 53, -3.3703848401898283), | |
39 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.28080574467981213), | |
40 | BOOST_MATH_BIG_CONSTANT(T, 53, 1.6638069578676164), | |
41 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.64468386819102836), | |
42 | }; | |
43 | static const T Q_1_2[] = { | |
44 | BOOST_MATH_BIG_CONSTANT(T, 53, 1.0), | |
45 | BOOST_MATH_BIG_CONSTANT(T, 53, 3.4535389668541151), | |
46 | BOOST_MATH_BIG_CONSTANT(T, 53, 4.5208926987851437), | |
47 | BOOST_MATH_BIG_CONSTANT(T, 53, 2.7012734178351534), | |
48 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.64468798399785611), | |
49 | BOOST_MATH_BIG_CONSTANT(T, 53, -0.20314516859987728e-6), | |
50 | }; | |
51 | // Max error in interpolated form: 1.159e-017 | |
52 | static const T P_2_4[] = { | |
53 | BOOST_MATH_BIG_CONSTANT(T, 53, -0.13803835004508849e-7), | |
54 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.50000049158540261), | |
55 | BOOST_MATH_BIG_CONSTANT(T, 53, 1.6077979838469348), | |
56 | BOOST_MATH_BIG_CONSTANT(T, 53, 2.5645435828098254), | |
57 | BOOST_MATH_BIG_CONSTANT(T, 53, 2.0534873203680393), | |
58 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.74566981111565923), | |
59 | }; | |
60 | static const T Q_2_4[] = { | |
61 | BOOST_MATH_BIG_CONSTANT(T, 53, 1.0), | |
62 | BOOST_MATH_BIG_CONSTANT(T, 53, 2.8822787662376169), | |
63 | BOOST_MATH_BIG_CONSTANT(T, 53, 4.1681660554090917), | |
64 | BOOST_MATH_BIG_CONSTANT(T, 53, 2.7853527819234466), | |
65 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.74967671848044792), | |
66 | BOOST_MATH_BIG_CONSTANT(T, 53, -0.00057069112416246805), | |
67 | }; | |
68 | // Maximum Deviation Found: 6.896e-018 | |
69 | // Expected Error Term : -6.895e-018 | |
70 | // Maximum Relative Change in Control Points : 8.497e-004 | |
71 | static const T P_4_inf[] = { | |
72 | static_cast<T>(0.68947581948701249e-17L), | |
73 | static_cast<T>(0.49999999999998975L), | |
74 | static_cast<T>(1.0177274392923795L), | |
75 | static_cast<T>(2.498208511343429L), | |
76 | static_cast<T>(2.1921221359427595L), | |
77 | static_cast<T>(1.5897035272532764L), | |
78 | static_cast<T>(0.40154388356961734L), | |
79 | }; | |
80 | static const T Q_4_inf[] = { | |
81 | static_cast<T>(1.0L), | |
82 | static_cast<T>(1.7021215452463932L), | |
83 | static_cast<T>(4.4290431747556469L), | |
84 | static_cast<T>(2.9745631894384922L), | |
85 | static_cast<T>(2.3013614809773616L), | |
86 | static_cast<T>(0.28360399799075752L), | |
87 | static_cast<T>(0.022892987908906897L), | |
88 | }; | |
89 | ||
90 | if(x <= 2) | |
91 | { | |
92 | return (offset + boost::math::tools::evaluate_polynomial(P_1_2, x) / tools::evaluate_polynomial(Q_1_2, x)) / (x * x); | |
93 | } | |
94 | else if(x <= 4) | |
95 | { | |
96 | T y = 1 / x; | |
97 | return (1 + tools::evaluate_polynomial(P_2_4, y) / tools::evaluate_polynomial(Q_2_4, y)) / x; | |
98 | } | |
99 | T y = 1 / x; | |
100 | return (1 + tools::evaluate_polynomial(P_4_inf, y) / tools::evaluate_polynomial(Q_4_inf, y)) / x; | |
101 | } | |
102 | ||
103 | template <class T, class Policy> | |
104 | T trigamma_prec(T x, const mpl::int_<64>*, const Policy&) | |
105 | { | |
106 | // Max error in interpolated form: 1.178e-020 | |
107 | static const T offset_1_2 = BOOST_MATH_BIG_CONSTANT(T, 64, 2.109325408935546875); | |
108 | static const T P_1_2[] = { | |
109 | BOOST_MATH_BIG_CONSTANT(T, 64, -1.10932535608960258341), | |
110 | BOOST_MATH_BIG_CONSTANT(T, 64, -4.18793841543017129052), | |
111 | BOOST_MATH_BIG_CONSTANT(T, 64, -4.63865531898487734531), | |
112 | BOOST_MATH_BIG_CONSTANT(T, 64, -0.919832884430500908047), | |
113 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.68074038333180423012), | |
114 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.21172611429185622377), | |
115 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.259635673503366427284), | |
116 | }; | |
117 | static const T Q_1_2[] = { | |
118 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.0), | |
119 | BOOST_MATH_BIG_CONSTANT(T, 64, 3.77521119359546982995), | |
120 | BOOST_MATH_BIG_CONSTANT(T, 64, 5.664338024578956321), | |
121 | BOOST_MATH_BIG_CONSTANT(T, 64, 4.25995134879278028361), | |
122 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.62956638448940402182), | |
123 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.259635512844691089868), | |
124 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.629642219810618032207e-8), | |
125 | }; | |
126 | // Max error in interpolated form: 3.912e-020 | |
127 | static const T P_2_8[] = { | |
128 | BOOST_MATH_BIG_CONSTANT(T, 64, -0.387540035162952880976e-11), | |
129 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.500000000276430504), | |
130 | BOOST_MATH_BIG_CONSTANT(T, 64, 3.21926880986360957306), | |
131 | BOOST_MATH_BIG_CONSTANT(T, 64, 10.2550347708483445775), | |
132 | BOOST_MATH_BIG_CONSTANT(T, 64, 18.9002075150709144043), | |
133 | BOOST_MATH_BIG_CONSTANT(T, 64, 21.0357215832399705625), | |
134 | BOOST_MATH_BIG_CONSTANT(T, 64, 13.4346512182925923978), | |
135 | BOOST_MATH_BIG_CONSTANT(T, 64, 3.98656291026448279118), | |
136 | }; | |
137 | static const T Q_2_8[] = { | |
138 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.0), | |
139 | BOOST_MATH_BIG_CONSTANT(T, 64, 6.10520430478613667724), | |
140 | BOOST_MATH_BIG_CONSTANT(T, 64, 18.475001060603645512), | |
141 | BOOST_MATH_BIG_CONSTANT(T, 64, 31.7087534567758405638), | |
142 | BOOST_MATH_BIG_CONSTANT(T, 64, 31.908814523890465398), | |
143 | BOOST_MATH_BIG_CONSTANT(T, 64, 17.4175479039227084798), | |
144 | BOOST_MATH_BIG_CONSTANT(T, 64, 3.98749106958394941276), | |
145 | BOOST_MATH_BIG_CONSTANT(T, 64, -0.000115917322224411128566), | |
146 | }; | |
147 | // Maximum Deviation Found: 2.635e-020 | |
148 | // Expected Error Term : 2.635e-020 | |
149 | // Maximum Relative Change in Control Points : 1.791e-003 | |
150 | static const T P_8_inf[] = { | |
151 | BOOST_MATH_BIG_CONSTANT(T, 64, -0.263527875092466899848e-19), | |
152 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.500000000000000058145), | |
153 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.0730121433777364138677), | |
154 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.94505878379957149534), | |
155 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.0517092358874932620529), | |
156 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.07995383547483921121), | |
157 | }; | |
158 | static const T Q_8_inf[] = { | |
159 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.0), | |
160 | BOOST_MATH_BIG_CONSTANT(T, 64, -0.187309046577818095504), | |
161 | BOOST_MATH_BIG_CONSTANT(T, 64, 3.95255391645238842975), | |
162 | BOOST_MATH_BIG_CONSTANT(T, 64, -1.14743283327078949087), | |
163 | BOOST_MATH_BIG_CONSTANT(T, 64, 2.52989799376344914499), | |
164 | BOOST_MATH_BIG_CONSTANT(T, 64, -0.627414303172402506396), | |
165 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.141554248216425512536), | |
166 | }; | |
167 | ||
168 | if(x <= 2) | |
169 | { | |
170 | return (offset_1_2 + boost::math::tools::evaluate_polynomial(P_1_2, x) / tools::evaluate_polynomial(Q_1_2, x)) / (x * x); | |
171 | } | |
172 | else if(x <= 8) | |
173 | { | |
174 | T y = 1 / x; | |
175 | return (1 + tools::evaluate_polynomial(P_2_8, y) / tools::evaluate_polynomial(Q_2_8, y)) / x; | |
176 | } | |
177 | T y = 1 / x; | |
178 | return (1 + tools::evaluate_polynomial(P_8_inf, y) / tools::evaluate_polynomial(Q_8_inf, y)) / x; | |
179 | } | |
180 | ||
181 | template <class T, class Policy> | |
182 | T trigamma_prec(T x, const mpl::int_<113>*, const Policy&) | |
183 | { | |
184 | // Max error in interpolated form: 1.916e-035 | |
185 | ||
186 | static const T P_1_2[] = { | |
187 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.999999999999999082554457936871832533), | |
188 | BOOST_MATH_BIG_CONSTANT(T, 113, -4.71237311120865266379041700054847734), | |
189 | BOOST_MATH_BIG_CONSTANT(T, 113, -7.94125711970499027763789342500817316), | |
190 | BOOST_MATH_BIG_CONSTANT(T, 113, -5.74657746697664735258222071695644535), | |
191 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.404213349456398905981223965160595687), | |
192 | BOOST_MATH_BIG_CONSTANT(T, 113, 2.47877781178642876561595890095758896), | |
193 | BOOST_MATH_BIG_CONSTANT(T, 113, 2.07714151702455125992166949812126433), | |
194 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.858877899162360138844032265418028567), | |
195 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.20499222604410032375789018837922397), | |
196 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0272103140348194747360175268778415049), | |
197 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0015764849020876949848954081173520686), | |
198 | }; | |
199 | static const T Q_1_2[] = { | |
200 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.0), | |
201 | BOOST_MATH_BIG_CONSTANT(T, 113, 4.71237311120863419878375031457715223), | |
202 | BOOST_MATH_BIG_CONSTANT(T, 113, 9.58619118655339853449127952145877467), | |
203 | BOOST_MATH_BIG_CONSTANT(T, 113, 11.0940067269829372437561421279054968), | |
204 | BOOST_MATH_BIG_CONSTANT(T, 113, 8.09075424749327792073276309969037885), | |
205 | BOOST_MATH_BIG_CONSTANT(T, 113, 3.87705890159891405185343806884451286), | |
206 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.22758678701914477836330837816976782), | |
207 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.249092040606385004109672077814668716), | |
208 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0295750413900655597027079600025569048), | |
209 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.00157648490200498142247694709728858139), | |
210 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.161264050344059471721062360645432809e-14), | |
211 | }; | |
212 | ||
213 | // Max error in interpolated form: 8.958e-035 | |
214 | static const T P_2_4[] = { | |
215 | BOOST_MATH_BIG_CONSTANT(T, 113, -2.55843734739907925764326773972215085), | |
216 | BOOST_MATH_BIG_CONSTANT(T, 113, -12.2830208240542011967952466273455887), | |
217 | BOOST_MATH_BIG_CONSTANT(T, 113, -23.9195022162767993526575786066414403), | |
218 | BOOST_MATH_BIG_CONSTANT(T, 113, -24.9256431504823483094158828285470862), | |
219 | BOOST_MATH_BIG_CONSTANT(T, 113, -14.7979122765478779075108064826412285), | |
220 | BOOST_MATH_BIG_CONSTANT(T, 113, -4.46654453928610666393276765059122272), | |
221 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.0191439033405649675717082465687845002), | |
222 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.515412052554351265708917209749037352), | |
223 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.195378348786064304378247325360320038), | |
224 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0334761282624174313035014426794245393), | |
225 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.002373665205942206348500250056602687), | |
226 | }; | |
227 | static const T Q_2_4[] = { | |
228 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.0), | |
229 | BOOST_MATH_BIG_CONSTANT(T, 113, 4.80098558454419907830670928248659245), | |
230 | BOOST_MATH_BIG_CONSTANT(T, 113, 9.99220727843170133895059300223445265), | |
231 | BOOST_MATH_BIG_CONSTANT(T, 113, 11.8896146167631330735386697123464976), | |
232 | BOOST_MATH_BIG_CONSTANT(T, 113, 8.96613256683809091593793565879092581), | |
233 | BOOST_MATH_BIG_CONSTANT(T, 113, 4.47254136149624110878909334574485751), | |
234 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.48600982028196527372434773913633152), | |
235 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.319570735766764237068541501137990078), | |
236 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0407358345787680953107374215319322066), | |
237 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.00237366520593271641375755486420859837), | |
238 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.239554887903526152679337256236302116e-15), | |
239 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.294749244740618656265237072002026314e-17), | |
240 | }; | |
241 | ||
242 | static const T y_offset_2_4 = BOOST_MATH_BIG_CONSTANT(T, 113, 3.558437347412109375); | |
243 | ||
244 | // Max error in interpolated form: 4.319e-035 | |
245 | static const T P_4_8[] = { | |
246 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.166626112697021464248967707021688845e-16), | |
247 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.499999999999997739552090249208808197), | |
248 | BOOST_MATH_BIG_CONSTANT(T, 113, 6.40270945019053817915772473771553187), | |
249 | BOOST_MATH_BIG_CONSTANT(T, 113, 41.3833374155000608013677627389343329), | |
250 | BOOST_MATH_BIG_CONSTANT(T, 113, 166.803341854562809335667241074035245), | |
251 | BOOST_MATH_BIG_CONSTANT(T, 113, 453.39964786925369319960722793414521), | |
252 | BOOST_MATH_BIG_CONSTANT(T, 113, 851.153712317697055375935433362983944), | |
253 | BOOST_MATH_BIG_CONSTANT(T, 113, 1097.70657567285059133109286478004458), | |
254 | BOOST_MATH_BIG_CONSTANT(T, 113, 938.431232478455316020076349367632922), | |
255 | BOOST_MATH_BIG_CONSTANT(T, 113, 487.268001604651932322080970189930074), | |
256 | BOOST_MATH_BIG_CONSTANT(T, 113, 119.953445242335730062471193124820659), | |
257 | }; | |
258 | static const T Q_4_8[] = { | |
259 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.0), | |
260 | BOOST_MATH_BIG_CONSTANT(T, 113, 12.4720855670474488978638945855932398), | |
261 | BOOST_MATH_BIG_CONSTANT(T, 113, 78.6093129753298570701376952709727391), | |
262 | BOOST_MATH_BIG_CONSTANT(T, 113, 307.470246050318322489781182863190127), | |
263 | BOOST_MATH_BIG_CONSTANT(T, 113, 805.140686101151538537565264188630079), | |
264 | BOOST_MATH_BIG_CONSTANT(T, 113, 1439.12019760292146454787601409644413), | |
265 | BOOST_MATH_BIG_CONSTANT(T, 113, 1735.6105285756048831268586001383127), | |
266 | BOOST_MATH_BIG_CONSTANT(T, 113, 1348.32500712856328019355198611280536), | |
267 | BOOST_MATH_BIG_CONSTANT(T, 113, 607.225985860570846699704222144650563), | |
268 | BOOST_MATH_BIG_CONSTANT(T, 113, 119.952317857277045332558673164517227), | |
269 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.000140165918355036060868680809129436084), | |
270 | }; | |
271 | ||
272 | // Maximum Deviation Found: 2.867e-035 | |
273 | // Expected Error Term : 2.866e-035 | |
274 | // Maximum Relative Change in Control Points : 2.662e-004 | |
275 | static const T P_8_16[] = { | |
276 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.184828315274146610610872315609837439e-19), | |
277 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.500000000000000004122475157735807738), | |
278 | BOOST_MATH_BIG_CONSTANT(T, 113, 3.02533865247313349284875558880415875), | |
279 | BOOST_MATH_BIG_CONSTANT(T, 113, 13.5995927517457371243039532492642734), | |
280 | BOOST_MATH_BIG_CONSTANT(T, 113, 35.3132224283087906757037999452941588), | |
281 | BOOST_MATH_BIG_CONSTANT(T, 113, 67.1639424550714159157603179911505619), | |
282 | BOOST_MATH_BIG_CONSTANT(T, 113, 83.5767733658513967581959839367419891), | |
283 | BOOST_MATH_BIG_CONSTANT(T, 113, 71.073491212235705900866411319363501), | |
284 | BOOST_MATH_BIG_CONSTANT(T, 113, 35.8621515614725564575893663483998663), | |
285 | BOOST_MATH_BIG_CONSTANT(T, 113, 8.72152231639983491987779743154333318), | |
286 | }; | |
287 | static const T Q_8_16[] = { | |
288 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.0), | |
289 | BOOST_MATH_BIG_CONSTANT(T, 113, 5.71734397161293452310624822415866372), | |
290 | BOOST_MATH_BIG_CONSTANT(T, 113, 25.293404179620438179337103263274815), | |
291 | BOOST_MATH_BIG_CONSTANT(T, 113, 62.2619767967468199111077640625328469), | |
292 | BOOST_MATH_BIG_CONSTANT(T, 113, 113.955048909238993473389714972250235), | |
293 | BOOST_MATH_BIG_CONSTANT(T, 113, 130.807138328938966981862203944329408), | |
294 | BOOST_MATH_BIG_CONSTANT(T, 113, 102.423146902337654110717764213057753), | |
295 | BOOST_MATH_BIG_CONSTANT(T, 113, 44.0424772805245202514468199602123565), | |
296 | BOOST_MATH_BIG_CONSTANT(T, 113, 8.89898032477904072082994913461386099), | |
297 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.0296627336872039988632793863671456398), | |
298 | }; | |
299 | // Maximum Deviation Found: 1.079e-035 | |
300 | // Expected Error Term : -1.079e-035 | |
301 | // Maximum Relative Change in Control Points : 7.884e-003 | |
302 | static const T P_16_inf[] = { | |
303 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0), | |
304 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.500000000000000000000000000000087317), | |
305 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.345625669885456215194494735902663968), | |
306 | BOOST_MATH_BIG_CONSTANT(T, 113, 9.62895499360842232127552650044647769), | |
307 | BOOST_MATH_BIG_CONSTANT(T, 113, 3.5936085382439026269301003761320812), | |
308 | BOOST_MATH_BIG_CONSTANT(T, 113, 49.459599118438883265036646019410669), | |
309 | BOOST_MATH_BIG_CONSTANT(T, 113, 7.77519237321893917784735690560496607), | |
310 | BOOST_MATH_BIG_CONSTANT(T, 113, 74.4536074488178075948642351179304121), | |
311 | BOOST_MATH_BIG_CONSTANT(T, 113, 2.75209340397069050436806159297952699), | |
312 | BOOST_MATH_BIG_CONSTANT(T, 113, 23.9292359711471667884504840186561598), | |
313 | }; | |
314 | static const T Q_16_inf[] = { | |
315 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.0), | |
316 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.357918006437579097055656138920742037), | |
317 | BOOST_MATH_BIG_CONSTANT(T, 113, 19.1386039850709849435325005484512944), | |
318 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.874349081464143606016221431763364517), | |
319 | BOOST_MATH_BIG_CONSTANT(T, 113, 98.6516097434855572678195488061432509), | |
320 | BOOST_MATH_BIG_CONSTANT(T, 113, -16.1051972833382893468655223662534306), | |
321 | BOOST_MATH_BIG_CONSTANT(T, 113, 154.316860216253720989145047141653727), | |
322 | BOOST_MATH_BIG_CONSTANT(T, 113, -40.2026880424378986053105969312264534), | |
323 | BOOST_MATH_BIG_CONSTANT(T, 113, 60.1679136674264778074736441126810223), | |
324 | BOOST_MATH_BIG_CONSTANT(T, 113, -13.3414844622256422644504472438320114), | |
325 | BOOST_MATH_BIG_CONSTANT(T, 113, 2.53795636200649908779512969030363442), | |
326 | }; | |
327 | ||
328 | if(x <= 2) | |
329 | { | |
330 | return (2 + boost::math::tools::evaluate_polynomial(P_1_2, x) / tools::evaluate_polynomial(Q_1_2, x)) / (x * x); | |
331 | } | |
332 | else if(x <= 4) | |
333 | { | |
334 | return (y_offset_2_4 + boost::math::tools::evaluate_polynomial(P_2_4, x) / tools::evaluate_polynomial(Q_2_4, x)) / (x * x); | |
335 | } | |
336 | else if(x <= 8) | |
337 | { | |
338 | T y = 1 / x; | |
339 | return (1 + tools::evaluate_polynomial(P_4_8, y) / tools::evaluate_polynomial(Q_4_8, y)) / x; | |
340 | } | |
341 | else if(x <= 16) | |
342 | { | |
343 | T y = 1 / x; | |
344 | return (1 + tools::evaluate_polynomial(P_8_16, y) / tools::evaluate_polynomial(Q_8_16, y)) / x; | |
345 | } | |
346 | T y = 1 / x; | |
347 | return (1 + tools::evaluate_polynomial(P_16_inf, y) / tools::evaluate_polynomial(Q_16_inf, y)) / x; | |
348 | } | |
349 | ||
350 | template <class T, class Tag, class Policy> | |
351 | T trigamma_imp(T x, const Tag* t, const Policy& pol) | |
352 | { | |
353 | // | |
354 | // This handles reflection of negative arguments, and all our | |
355 | // error handling, then forwards to the T-specific approximation. | |
356 | // | |
357 | BOOST_MATH_STD_USING // ADL of std functions. | |
358 | ||
359 | T result = 0; | |
360 | // | |
361 | // Check for negative arguments and use reflection: | |
362 | // | |
363 | if(x <= 0) | |
364 | { | |
365 | // Reflect: | |
366 | T z = 1 - x; | |
367 | // Argument reduction for tan: | |
368 | if(floor(x) == x) | |
369 | { | |
370 | return policies::raise_pole_error<T>("boost::math::trigamma<%1%>(%1%)", 0, (1-x), pol); | |
371 | } | |
372 | T s = fabs(x) < fabs(z) ? boost::math::sin_pi(x, pol) : boost::math::sin_pi(z, pol); | |
373 | return -trigamma_imp(z, t, pol) + boost::math::pow<2>(constants::pi<T>()) / (s * s); | |
374 | } | |
375 | if(x < 1) | |
376 | { | |
377 | result = 1 / (x * x); | |
378 | x += 1; | |
379 | } | |
380 | return result + trigamma_prec(x, t, pol); | |
381 | } | |
382 | ||
383 | template <class T, class Policy> | |
384 | T trigamma_imp(T x, const mpl::int_<0>*, const Policy& pol) | |
385 | { | |
386 | return polygamma_imp(1, x, pol); | |
387 | } | |
388 | // | |
389 | // Initializer: ensure all our constants are initialized prior to the first call of main: | |
390 | // | |
391 | template <class T, class Policy> | |
392 | struct trigamma_initializer | |
393 | { | |
394 | struct init | |
395 | { | |
396 | init() | |
397 | { | |
398 | typedef typename policies::precision<T, Policy>::type precision_type; | |
399 | do_init(mpl::bool_<precision_type::value && (precision_type::value <= 113)>()); | |
400 | } | |
401 | void do_init(const mpl::true_&) | |
402 | { | |
403 | boost::math::trigamma(T(2.5), Policy()); | |
404 | } | |
405 | void do_init(const mpl::false_&){} | |
406 | void force_instantiate()const{} | |
407 | }; | |
408 | static const init initializer; | |
409 | static void force_instantiate() | |
410 | { | |
411 | initializer.force_instantiate(); | |
412 | } | |
413 | }; | |
414 | ||
415 | template <class T, class Policy> | |
416 | const typename trigamma_initializer<T, Policy>::init trigamma_initializer<T, Policy>::initializer; | |
417 | ||
418 | } // namespace detail | |
419 | ||
420 | template <class T, class Policy> | |
421 | inline typename tools::promote_args<T>::type | |
422 | trigamma(T x, const Policy&) | |
423 | { | |
424 | typedef typename tools::promote_args<T>::type result_type; | |
425 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
426 | typedef typename policies::precision<T, Policy>::type precision_type; | |
427 | typedef typename mpl::if_< | |
428 | mpl::or_< | |
429 | mpl::less_equal<precision_type, mpl::int_<0> >, | |
430 | mpl::greater<precision_type, mpl::int_<114> > | |
431 | >, | |
432 | mpl::int_<0>, | |
433 | typename mpl::if_< | |
434 | mpl::less<precision_type, mpl::int_<54> >, | |
435 | mpl::int_<53>, | |
436 | typename mpl::if_< | |
437 | mpl::less<precision_type, mpl::int_<65> >, | |
438 | mpl::int_<64>, | |
439 | mpl::int_<113> | |
440 | >::type | |
441 | >::type | |
442 | >::type tag_type; | |
443 | ||
444 | typedef typename policies::normalise< | |
445 | Policy, | |
446 | policies::promote_float<false>, | |
447 | policies::promote_double<false>, | |
448 | policies::discrete_quantile<>, | |
449 | policies::assert_undefined<> >::type forwarding_policy; | |
450 | ||
451 | // Force initialization of constants: | |
452 | detail::trigamma_initializer<value_type, forwarding_policy>::force_instantiate(); | |
453 | ||
454 | return policies::checked_narrowing_cast<result_type, Policy>(detail::trigamma_imp( | |
455 | static_cast<value_type>(x), | |
456 | static_cast<const tag_type*>(0), forwarding_policy()), "boost::math::trigamma<%1%>(%1%)"); | |
457 | } | |
458 | ||
459 | template <class T> | |
460 | inline typename tools::promote_args<T>::type | |
461 | trigamma(T x) | |
462 | { | |
463 | return trigamma(x, policies::policy<>()); | |
464 | } | |
465 | ||
466 | } // namespace math | |
467 | } // namespace boost | |
468 | #endif | |
469 |