]>
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_SPECIAL_FUNCTIONS_DETAIL_LGAMMA_SMALL | |
7 | #define BOOST_MATH_SPECIAL_FUNCTIONS_DETAIL_LGAMMA_SMALL | |
8 | ||
9 | #ifdef _MSC_VER | |
10 | #pragma once | |
11 | #endif | |
12 | ||
13 | #include <boost/math/tools/big_constant.hpp> | |
14 | ||
92f5a8d4 TL |
15 | #if defined(__GNUC__) && defined(BOOST_MATH_USE_FLOAT128) |
16 | // | |
17 | // This is the only way we can avoid | |
18 | // warning: non-standard suffix on floating constant [-Wpedantic] | |
19 | // when building with -Wall -pedantic. Neither __extension__ | |
f67539c2 | 20 | // nor #pragma diagnostic ignored work :( |
92f5a8d4 TL |
21 | // |
22 | #pragma GCC system_header | |
23 | #endif | |
24 | ||
7c673cae FG |
25 | namespace boost{ namespace math{ namespace detail{ |
26 | ||
27 | // | |
28 | // These need forward declaring to keep GCC happy: | |
29 | // | |
30 | template <class T, class Policy, class Lanczos> | |
31 | T gamma_imp(T z, const Policy& pol, const Lanczos& l); | |
32 | template <class T, class Policy> | |
33 | T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos& l); | |
34 | ||
35 | // | |
36 | // lgamma for small arguments: | |
37 | // | |
38 | template <class T, class Policy, class Lanczos> | |
f67539c2 | 39 | T lgamma_small_imp(T z, T zm1, T zm2, const boost::integral_constant<int, 64>&, const Policy& /* l */, const Lanczos&) |
7c673cae FG |
40 | { |
41 | // This version uses rational approximations for small | |
42 | // values of z accurate enough for 64-bit mantissas | |
43 | // (80-bit long doubles), works well for 53-bit doubles as well. | |
44 | // Lanczos is only used to select the Lanczos function. | |
45 | ||
46 | BOOST_MATH_STD_USING // for ADL of std names | |
47 | T result = 0; | |
48 | if(z < tools::epsilon<T>()) | |
49 | { | |
50 | result = -log(z); | |
51 | } | |
52 | else if((zm1 == 0) || (zm2 == 0)) | |
53 | { | |
54 | // nothing to do, result is zero.... | |
55 | } | |
56 | else if(z > 2) | |
57 | { | |
58 | // | |
59 | // Begin by performing argument reduction until | |
60 | // z is in [2,3): | |
61 | // | |
62 | if(z >= 3) | |
63 | { | |
64 | do | |
65 | { | |
66 | z -= 1; | |
67 | zm2 -= 1; | |
68 | result += log(z); | |
69 | }while(z >= 3); | |
70 | // Update zm2, we need it below: | |
71 | zm2 = z - 2; | |
72 | } | |
73 | ||
74 | // | |
75 | // Use the following form: | |
76 | // | |
77 | // lgamma(z) = (z-2)(z+1)(Y + R(z-2)) | |
78 | // | |
79 | // where R(z-2) is a rational approximation optimised for | |
80 | // low absolute error - as long as it's absolute error | |
81 | // is small compared to the constant Y - then any rounding | |
82 | // error in it's computation will get wiped out. | |
83 | // | |
84 | // R(z-2) has the following properties: | |
85 | // | |
86 | // At double: Max error found: 4.231e-18 | |
87 | // At long double: Max error found: 1.987e-21 | |
88 | // Maximum Deviation Found (approximation error): 5.900e-24 | |
89 | // | |
90 | static const T P[] = { | |
91 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.180355685678449379109e-1)), | |
92 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.25126649619989678683e-1)), | |
93 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.494103151567532234274e-1)), | |
94 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.172491608709613993966e-1)), | |
95 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.259453563205438108893e-3)), | |
96 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.541009869215204396339e-3)), | |
97 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.324588649825948492091e-4)) | |
98 | }; | |
99 | static const T Q[] = { | |
100 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.1e1)), | |
101 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.196202987197795200688e1)), | |
102 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.148019669424231326694e1)), | |
103 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.541391432071720958364e0)), | |
104 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.988504251128010129477e-1)), | |
105 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.82130967464889339326e-2)), | |
106 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.224936291922115757597e-3)), | |
107 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.223352763208617092964e-6)) | |
108 | }; | |
109 | ||
110 | static const float Y = 0.158963680267333984375e0f; | |
111 | ||
112 | T r = zm2 * (z + 1); | |
113 | T R = tools::evaluate_polynomial(P, zm2); | |
114 | R /= tools::evaluate_polynomial(Q, zm2); | |
115 | ||
116 | result += r * Y + r * R; | |
117 | } | |
118 | else | |
119 | { | |
120 | // | |
f67539c2 | 121 | // If z is less than 1 use recurrence to shift to |
7c673cae FG |
122 | // z in the interval [1,2]: |
123 | // | |
124 | if(z < 1) | |
125 | { | |
126 | result += -log(z); | |
127 | zm2 = zm1; | |
128 | zm1 = z; | |
129 | z += 1; | |
130 | } | |
131 | // | |
132 | // Two approximations, on for z in [1,1.5] and | |
133 | // one for z in [1.5,2]: | |
134 | // | |
135 | if(z <= 1.5) | |
136 | { | |
137 | // | |
138 | // Use the following form: | |
139 | // | |
140 | // lgamma(z) = (z-1)(z-2)(Y + R(z-1)) | |
141 | // | |
142 | // where R(z-1) is a rational approximation optimised for | |
143 | // low absolute error - as long as it's absolute error | |
144 | // is small compared to the constant Y - then any rounding | |
145 | // error in it's computation will get wiped out. | |
146 | // | |
147 | // R(z-1) has the following properties: | |
148 | // | |
149 | // At double precision: Max error found: 1.230011e-17 | |
150 | // At 80-bit long double precision: Max error found: 5.631355e-21 | |
151 | // Maximum Deviation Found: 3.139e-021 | |
152 | // Expected Error Term: 3.139e-021 | |
153 | ||
154 | // | |
155 | static const float Y = 0.52815341949462890625f; | |
156 | ||
157 | static const T P[] = { | |
158 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.490622454069039543534e-1)), | |
159 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.969117530159521214579e-1)), | |
160 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.414983358359495381969e0)), | |
161 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.406567124211938417342e0)), | |
162 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.158413586390692192217e0)), | |
163 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.240149820648571559892e-1)), | |
164 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.100346687696279557415e-2)) | |
165 | }; | |
166 | static const T Q[] = { | |
167 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.1e1)), | |
168 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.302349829846463038743e1)), | |
169 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.348739585360723852576e1)), | |
170 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.191415588274426679201e1)), | |
171 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.507137738614363510846e0)), | |
172 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.577039722690451849648e-1)), | |
173 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.195768102601107189171e-2)) | |
174 | }; | |
175 | ||
176 | T r = tools::evaluate_polynomial(P, zm1) / tools::evaluate_polynomial(Q, zm1); | |
177 | T prefix = zm1 * zm2; | |
178 | ||
179 | result += prefix * Y + prefix * r; | |
180 | } | |
181 | else | |
182 | { | |
183 | // | |
184 | // Use the following form: | |
185 | // | |
186 | // lgamma(z) = (2-z)(1-z)(Y + R(2-z)) | |
187 | // | |
188 | // where R(2-z) is a rational approximation optimised for | |
189 | // low absolute error - as long as it's absolute error | |
190 | // is small compared to the constant Y - then any rounding | |
191 | // error in it's computation will get wiped out. | |
192 | // | |
193 | // R(2-z) has the following properties: | |
194 | // | |
195 | // At double precision, max error found: 1.797565e-17 | |
196 | // At 80-bit long double precision, max error found: 9.306419e-21 | |
197 | // Maximum Deviation Found: 2.151e-021 | |
198 | // Expected Error Term: 2.150e-021 | |
199 | // | |
200 | static const float Y = 0.452017307281494140625f; | |
201 | ||
202 | static const T P[] = { | |
203 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.292329721830270012337e-1)), | |
204 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.144216267757192309184e0)), | |
205 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.142440390738631274135e0)), | |
206 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.542809694055053558157e-1)), | |
207 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.850535976868336437746e-2)), | |
208 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.431171342679297331241e-3)) | |
209 | }; | |
210 | static const T Q[] = { | |
211 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.1e1)), | |
212 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.150169356054485044494e1)), | |
213 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.846973248876495016101e0)), | |
214 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.220095151814995745555e0)), | |
215 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.25582797155975869989e-1)), | |
216 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.100666795539143372762e-2)), | |
217 | static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.827193521891290553639e-6)) | |
218 | }; | |
219 | T r = zm2 * zm1; | |
220 | T R = tools::evaluate_polynomial(P, T(-zm2)) / tools::evaluate_polynomial(Q, T(-zm2)); | |
221 | ||
222 | result += r * Y + r * R; | |
223 | } | |
224 | } | |
225 | return result; | |
226 | } | |
227 | template <class T, class Policy, class Lanczos> | |
f67539c2 | 228 | T lgamma_small_imp(T z, T zm1, T zm2, const boost::integral_constant<int, 113>&, const Policy& /* l */, const Lanczos&) |
7c673cae FG |
229 | { |
230 | // | |
231 | // This version uses rational approximations for small | |
232 | // values of z accurate enough for 113-bit mantissas | |
233 | // (128-bit long doubles). | |
234 | // | |
235 | BOOST_MATH_STD_USING // for ADL of std names | |
236 | T result = 0; | |
237 | if(z < tools::epsilon<T>()) | |
238 | { | |
239 | result = -log(z); | |
240 | BOOST_MATH_INSTRUMENT_CODE(result); | |
241 | } | |
242 | else if((zm1 == 0) || (zm2 == 0)) | |
243 | { | |
244 | // nothing to do, result is zero.... | |
245 | } | |
246 | else if(z > 2) | |
247 | { | |
248 | // | |
249 | // Begin by performing argument reduction until | |
250 | // z is in [2,3): | |
251 | // | |
252 | if(z >= 3) | |
253 | { | |
254 | do | |
255 | { | |
256 | z -= 1; | |
257 | result += log(z); | |
258 | }while(z >= 3); | |
259 | zm2 = z - 2; | |
260 | } | |
261 | BOOST_MATH_INSTRUMENT_CODE(zm2); | |
262 | BOOST_MATH_INSTRUMENT_CODE(z); | |
263 | BOOST_MATH_INSTRUMENT_CODE(result); | |
264 | ||
265 | // | |
266 | // Use the following form: | |
267 | // | |
268 | // lgamma(z) = (z-2)(z+1)(Y + R(z-2)) | |
269 | // | |
270 | // where R(z-2) is a rational approximation optimised for | |
271 | // low absolute error - as long as it's absolute error | |
272 | // is small compared to the constant Y - then any rounding | |
273 | // error in it's computation will get wiped out. | |
274 | // | |
275 | // Maximum Deviation Found (approximation error) 3.73e-37 | |
276 | ||
277 | static const T P[] = { | |
278 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.018035568567844937910504030027467476655), | |
279 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.013841458273109517271750705401202404195), | |
280 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.062031842739486600078866923383017722399), | |
281 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.052518418329052161202007865149435256093), | |
282 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.01881718142472784129191838493267755758), | |
283 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0025104830367021839316463675028524702846), | |
284 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.00021043176101831873281848891452678568311), | |
285 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.00010249622350908722793327719494037981166), | |
286 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.11381479670982006841716879074288176994e-4), | |
287 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.49999811718089980992888533630523892389e-6), | |
288 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.70529798686542184668416911331718963364e-8) | |
289 | }; | |
290 | static const T Q[] = { | |
291 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.0), | |
292 | BOOST_MATH_BIG_CONSTANT(T, 113, 2.5877485070422317542808137697939233685), | |
293 | BOOST_MATH_BIG_CONSTANT(T, 113, 2.8797959228352591788629602533153837126), | |
294 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.8030885955284082026405495275461180977), | |
295 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.69774331297747390169238306148355428436), | |
296 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.17261566063277623942044077039756583802), | |
297 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.02729301254544230229429621192443000121), | |
298 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0026776425891195270663133581960016620433), | |
299 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.00015244249160486584591370355730402168106), | |
300 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.43997034032479866020546814475414346627e-5), | |
301 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.46295080708455613044541885534408170934e-7), | |
302 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.93326638207459533682980757982834180952e-11), | |
303 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.42316456553164995177177407325292867513e-13) | |
304 | }; | |
305 | ||
306 | T R = tools::evaluate_polynomial(P, zm2); | |
307 | R /= tools::evaluate_polynomial(Q, zm2); | |
308 | ||
309 | static const float Y = 0.158963680267333984375F; | |
310 | ||
311 | T r = zm2 * (z + 1); | |
312 | ||
313 | result += r * Y + r * R; | |
314 | BOOST_MATH_INSTRUMENT_CODE(result); | |
315 | } | |
316 | else | |
317 | { | |
318 | // | |
f67539c2 | 319 | // If z is less than 1 use recurrence to shift to |
7c673cae FG |
320 | // z in the interval [1,2]: |
321 | // | |
322 | if(z < 1) | |
323 | { | |
324 | result += -log(z); | |
325 | zm2 = zm1; | |
326 | zm1 = z; | |
327 | z += 1; | |
328 | } | |
329 | BOOST_MATH_INSTRUMENT_CODE(result); | |
330 | BOOST_MATH_INSTRUMENT_CODE(z); | |
331 | BOOST_MATH_INSTRUMENT_CODE(zm2); | |
332 | // | |
333 | // Three approximations, on for z in [1,1.35], [1.35,1.625] and [1.625,1] | |
334 | // | |
335 | if(z <= 1.35) | |
336 | { | |
337 | // | |
338 | // Use the following form: | |
339 | // | |
340 | // lgamma(z) = (z-1)(z-2)(Y + R(z-1)) | |
341 | // | |
342 | // where R(z-1) is a rational approximation optimised for | |
343 | // low absolute error - as long as it's absolute error | |
344 | // is small compared to the constant Y - then any rounding | |
345 | // error in it's computation will get wiped out. | |
346 | // | |
347 | // R(z-1) has the following properties: | |
348 | // | |
349 | // Maximum Deviation Found (approximation error) 1.659e-36 | |
350 | // Expected Error Term (theoretical error) 1.343e-36 | |
351 | // Max error found at 128-bit long double precision 1.007e-35 | |
352 | // | |
353 | static const float Y = 0.54076099395751953125f; | |
354 | ||
355 | static const T P[] = { | |
356 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.036454670944013329356512090082402429697), | |
357 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.066235835556476033710068679907798799959), | |
358 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.67492399795577182387312206593595565371), | |
359 | BOOST_MATH_BIG_CONSTANT(T, 113, -1.4345555263962411429855341651960000166), | |
360 | BOOST_MATH_BIG_CONSTANT(T, 113, -1.4894319559821365820516771951249649563), | |
361 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.87210277668067964629483299712322411566), | |
362 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.29602090537771744401524080430529369136), | |
363 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.0561832587517836908929331992218879676), | |
364 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.0053236785487328044334381502530383140443), | |
365 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.00018629360291358130461736386077971890789), | |
366 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.10164985672213178500790406939467614498e-6), | |
367 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.13680157145361387405588201461036338274e-8) | |
368 | }; | |
369 | static const T Q[] = { | |
370 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.0), | |
371 | BOOST_MATH_BIG_CONSTANT(T, 113, 4.9106336261005990534095838574132225599), | |
372 | BOOST_MATH_BIG_CONSTANT(T, 113, 10.258804800866438510889341082793078432), | |
373 | BOOST_MATH_BIG_CONSTANT(T, 113, 11.88588976846826108836629960537466889), | |
374 | BOOST_MATH_BIG_CONSTANT(T, 113, 8.3455000546999704314454891036700998428), | |
375 | BOOST_MATH_BIG_CONSTANT(T, 113, 3.6428823682421746343233362007194282703), | |
376 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.97465989807254572142266753052776132252), | |
377 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.15121052897097822172763084966793352524), | |
378 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.012017363555383555123769849654484594893), | |
379 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0003583032812720649835431669893011257277) | |
380 | }; | |
381 | ||
382 | T r = tools::evaluate_polynomial(P, zm1) / tools::evaluate_polynomial(Q, zm1); | |
383 | T prefix = zm1 * zm2; | |
384 | ||
385 | result += prefix * Y + prefix * r; | |
386 | BOOST_MATH_INSTRUMENT_CODE(result); | |
387 | } | |
388 | else if(z <= 1.625) | |
389 | { | |
390 | // | |
391 | // Use the following form: | |
392 | // | |
393 | // lgamma(z) = (2-z)(1-z)(Y + R(2-z)) | |
394 | // | |
395 | // where R(2-z) is a rational approximation optimised for | |
396 | // low absolute error - as long as it's absolute error | |
397 | // is small compared to the constant Y - then any rounding | |
398 | // error in it's computation will get wiped out. | |
399 | // | |
400 | // R(2-z) has the following properties: | |
401 | // | |
402 | // Max error found at 128-bit long double precision 9.634e-36 | |
403 | // Maximum Deviation Found (approximation error) 1.538e-37 | |
404 | // Expected Error Term (theoretical error) 2.350e-38 | |
405 | // | |
406 | static const float Y = 0.483787059783935546875f; | |
407 | ||
408 | static const T P[] = { | |
409 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.017977422421608624353488126610933005432), | |
410 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.18484528905298309555089509029244135703), | |
411 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.40401251514859546989565001431430884082), | |
412 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.40277179799147356461954182877921388182), | |
413 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.21993421441282936476709677700477598816), | |
414 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.069595742223850248095697771331107571011), | |
415 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.012681481427699686635516772923547347328), | |
416 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0012489322866834830413292771335113136034), | |
417 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.57058739515423112045108068834668269608e-4), | |
418 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.8207548771933585614380644961342925976e-6) | |
419 | }; | |
420 | static const T Q[] = { | |
421 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.0), | |
422 | BOOST_MATH_BIG_CONSTANT(T, 113, -2.9629552288944259229543137757200262073), | |
423 | BOOST_MATH_BIG_CONSTANT(T, 113, 3.7118380799042118987185957298964772755), | |
424 | BOOST_MATH_BIG_CONSTANT(T, 113, -2.5569815272165399297600586376727357187), | |
425 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.0546764918220835097855665680632153367), | |
426 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.26574021300894401276478730940980810831), | |
427 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.03996289731752081380552901986471233462), | |
428 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.0033398680924544836817826046380586480873), | |
429 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.00013288854760548251757651556792598235735), | |
430 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.17194794958274081373243161848194745111e-5) | |
431 | }; | |
432 | T r = zm2 * zm1; | |
433 | T R = tools::evaluate_polynomial(P, T(0.625 - zm1)) / tools::evaluate_polynomial(Q, T(0.625 - zm1)); | |
434 | ||
435 | result += r * Y + r * R; | |
436 | BOOST_MATH_INSTRUMENT_CODE(result); | |
437 | } | |
438 | else | |
439 | { | |
440 | // | |
441 | // Same form as above. | |
442 | // | |
443 | // Max error found (at 128-bit long double precision) 1.831e-35 | |
444 | // Maximum Deviation Found (approximation error) 8.588e-36 | |
445 | // Expected Error Term (theoretical error) 1.458e-36 | |
446 | // | |
447 | static const float Y = 0.443811893463134765625f; | |
448 | ||
449 | static const T P[] = { | |
450 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.021027558364667626231512090082402429494), | |
451 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.15128811104498736604523586803722368377), | |
452 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.26249631480066246699388544451126410278), | |
453 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.21148748610533489823742352180628489742), | |
454 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.093964130697489071999873506148104370633), | |
455 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.024292059227009051652542804957550866827), | |
456 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.0036284453226534839926304745756906117066), | |
457 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0002939230129315195346843036254392485984), | |
458 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.11088589183158123733132268042570710338e-4), | |
459 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.13240510580220763969511741896361984162e-6) | |
460 | }; | |
461 | static const T Q[] = { | |
462 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.0), | |
463 | BOOST_MATH_BIG_CONSTANT(T, 113, -2.4240003754444040525462170802796471996), | |
464 | BOOST_MATH_BIG_CONSTANT(T, 113, 2.4868383476933178722203278602342786002), | |
465 | BOOST_MATH_BIG_CONSTANT(T, 113, -1.4047068395206343375520721509193698547), | |
466 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.47583809087867443858344765659065773369), | |
467 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.09865724264554556400463655444270700132), | |
468 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.012238223514176587501074150988445109735), | |
469 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.00084625068418239194670614419707491797097), | |
470 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.2796574430456237061420839429225710602e-4), | |
471 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.30202973883316730694433702165188835331e-6) | |
472 | }; | |
473 | // (2 - x) * (1 - x) * (c + R(2 - x)) | |
474 | T r = zm2 * zm1; | |
475 | T R = tools::evaluate_polynomial(P, T(-zm2)) / tools::evaluate_polynomial(Q, T(-zm2)); | |
476 | ||
477 | result += r * Y + r * R; | |
478 | BOOST_MATH_INSTRUMENT_CODE(result); | |
479 | } | |
480 | } | |
481 | BOOST_MATH_INSTRUMENT_CODE(result); | |
482 | return result; | |
483 | } | |
484 | template <class T, class Policy, class Lanczos> | |
f67539c2 | 485 | T lgamma_small_imp(T z, T zm1, T zm2, const boost::integral_constant<int, 0>&, const Policy& pol, const Lanczos&) |
7c673cae FG |
486 | { |
487 | // | |
488 | // No rational approximations are available because either | |
489 | // T has no numeric_limits support (so we can't tell how | |
490 | // many digits it has), or T has more digits than we know | |
491 | // what to do with.... we do have a Lanczos approximation | |
492 | // though, and that can be used to keep errors under control. | |
493 | // | |
494 | BOOST_MATH_STD_USING // for ADL of std names | |
495 | T result = 0; | |
496 | if(z < tools::epsilon<T>()) | |
497 | { | |
498 | result = -log(z); | |
499 | } | |
500 | else if(z < 0.5) | |
501 | { | |
502 | // taking the log of tgamma reduces the error, no danger of overflow here: | |
503 | result = log(gamma_imp(z, pol, Lanczos())); | |
504 | } | |
505 | else if(z >= 3) | |
506 | { | |
507 | // taking the log of tgamma reduces the error, no danger of overflow here: | |
508 | result = log(gamma_imp(z, pol, Lanczos())); | |
509 | } | |
510 | else if(z >= 1.5) | |
511 | { | |
512 | // special case near 2: | |
513 | T dz = zm2; | |
514 | result = dz * log((z + Lanczos::g() - T(0.5)) / boost::math::constants::e<T>()); | |
515 | result += boost::math::log1p(dz / (Lanczos::g() + T(1.5)), pol) * T(1.5); | |
516 | result += boost::math::log1p(Lanczos::lanczos_sum_near_2(dz), pol); | |
517 | } | |
518 | else | |
519 | { | |
520 | // special case near 1: | |
521 | T dz = zm1; | |
522 | result = dz * log((z + Lanczos::g() - T(0.5)) / boost::math::constants::e<T>()); | |
523 | result += boost::math::log1p(dz / (Lanczos::g() + T(0.5)), pol) / 2; | |
524 | result += boost::math::log1p(Lanczos::lanczos_sum_near_1(dz), pol); | |
525 | } | |
526 | return result; | |
527 | } | |
528 | ||
529 | }}} // namespaces | |
530 | ||
531 | #endif // BOOST_MATH_SPECIAL_FUNCTIONS_DETAIL_LGAMMA_SMALL | |
532 |