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