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 | |
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__ |
20 | // nor #pragma diagnostic ignored work :( |
21 | // |
22 | #pragma GCC system_header |
23 | #endif |
24 | |
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> |
39 | T lgamma_small_imp(T z, T zm1, T zm2, const boost::integral_constant<int, 64>&, const Policy& /* l */, const Lanczos&) |
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 | // |
121 | // If z is less than 1 use recurrence to shift to |
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> |
228 | T lgamma_small_imp(T z, T zm1, T zm2, const boost::integral_constant<int, 113>&, const Policy& /* l */, const Lanczos&) |
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 | // |
319 | // If z is less than 1 use recurrence to shift to |
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> |
485 | T lgamma_small_imp(T z, T zm1, T zm2, const boost::integral_constant<int, 0>&, const Policy& pol, const Lanczos&) |
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 | |
533 | |