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
25namespace boost{ namespace math{ namespace detail{
26
27//
28// These need forward declaring to keep GCC happy:
29//
30template <class T, class Policy, class Lanczos>
31T gamma_imp(T z, const Policy& pol, const Lanczos& l);
32template <class T, class Policy>
33T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos& l);
34
35//
36// lgamma for small arguments:
37//
38template <class T, class Policy, class Lanczos>
39T 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}
227template <class T, class Policy, class Lanczos>
228T 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}
484template <class T, class Policy, class Lanczos>
485T 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

source code of include/boost/math/special_functions/detail/lgamma_small.hpp