RandLib 2.0.0
Probability distributions, fitting, and heavy-tail tools for C++20
Loading...
Searching...
No Matches
math/gamma.hpp
Go to the documentation of this file.
1#pragma once
2
3#include <concepts>
4#include <type_traits>
5
6#include <randlib/math/constants.hpp>
7
8namespace randlib
9{
10
12
13namespace math
14{
15
22template <std::floating_point Real>
23Real factorial(Real n);
24
25template <std::integral Integer>
26long double factorial(Integer n);
27
33template <std::floating_point Real = long double>
34Real doubleFactorial(int n);
35
42template <std::floating_point Real = long double>
43Real lfact(size_t n);
44
50template <std::floating_point Real = long double>
51Real ldfact(size_t n);
52
59template <std::floating_point Real = long double>
60Real binom(size_t n, size_t k);
61
66template <std::floating_point Real = double>
67Real stirlingError(size_t n);
68
74template <std::floating_point Real = double>
75Real logFactorialRatio(size_t x, size_t y);
76
84template <std::floating_point Real>
85Real bd0(std::type_identity_t<Real> x, Real m);
86
91template <std::floating_point Real>
92Real digamma(Real x);
93
98template <std::floating_point Real>
99Real digammamLog(Real x);
100
105template <std::floating_point Real>
106Real trigamma(Real x);
107
117template <std::floating_point Real>
118Real lpgamma(std::type_identity_t<Real> a, Real x, std::type_identity_t<Real> logA, std::type_identity_t<Real> lgammaA);
119
128template <std::floating_point Real>
129Real lpgamma(std::type_identity_t<Real> a, Real x, std::type_identity_t<Real> logX);
130
137template <std::floating_point Real>
138Real lpgamma(std::type_identity_t<Real> a, Real x);
139
149template <std::floating_point Real>
150Real pgamma(std::type_identity_t<Real> a, Real x, std::type_identity_t<Real> logA, std::type_identity_t<Real> lgammaA);
151
160template <std::floating_point Real>
161Real pgamma(std::type_identity_t<Real> a, Real x, std::type_identity_t<Real> logX);
162
169template <std::floating_point Real>
170Real pgamma(std::type_identity_t<Real> a, Real x);
171
181template <std::floating_point Real>
182Real lqgamma(std::type_identity_t<Real> a, Real x, std::type_identity_t<Real> logA, std::type_identity_t<Real> lgammaA);
183
192template <std::floating_point Real>
193Real lqgamma(std::type_identity_t<Real> a, Real x, std::type_identity_t<Real> logX);
194
201template <std::floating_point Real>
202Real lqgamma(std::type_identity_t<Real> a, Real x);
203
213template <std::floating_point Real>
214Real qgamma(std::type_identity_t<Real> a, Real x, std::type_identity_t<Real> logA, std::type_identity_t<Real> lgammaA);
215
224template <std::floating_point Real>
225Real qgamma(std::type_identity_t<Real> a, Real x, std::type_identity_t<Real> logX);
226
233template <std::floating_point Real>
234Real qgamma(std::type_identity_t<Real> a, Real x);
235} // namespace math
236
237} // namespace randlib
238
239#include <randlib/math/functions.hpp>
240#include <randlib/math/beta.hpp>
241#include <randlib/math/quadrature.hpp>
242#include <randlib/math/roots.hpp>
243#include <randlib/math/optimize.hpp>
244#include <randlib/math/special_functions.hpp>
245
246namespace randlib
247{
248
249namespace math
250{
251
254constexpr int FACTORIAL_TABLESIZE = 175;
255
258constexpr long double FACTORIAL_TABLE[] = {
259 1.l,
260 3628800.l,
261 2432902008176640000.l,
262 265252859812191058636308480000000.l,
263 815915283247897734345611269596115894272000000000.l,
264 30414093201713378043612608166064768844377641568960512000000000000.l,
265 8320987112741390144276341183223364380754172606361245952449277696409600000000000000.l,
266 11978571669969891796072783721689098736458938142546425857555362864628009582789845319680000000000000000.l,
267 71569457046263802294811533723186532165584657342365752577109445058227039255480148842668944867280814080000000000000000000.l,
268 1485715964481761497309522733620825737885569961284688766942216863704985393094065876545992131370884059645617234469978112000000000000000000000.l,
269 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000.l,
270 15882455415227429404253703127090772871724410234473563207581748318444567162948183030959960131517678520479243672638179990208521148623422266876757623911219200000000000000000000000000.l,
271 6689502913449127057588118054090372586752746333138029810295671352301633557244962989366874165271984981308157637893214090552534408589408121859898481114389650005964960521256960000000000000000000000000000.l,
272 6466855489220473672507304395536485253155359447828049608975952322944781961185526165512707047229268452925683969240398027149120740074042105844737747799459310029635780991774612983803150965145600000000000000000000000000000000.l,
273 13462012475717524605876073858941615558355851148193967190051391468057460367090535696797920946629681836680869097041958983702264048370902871114013579941370766400374327741701139895604871545254810788060989321379840000000000000000000000000000000000.l,
274 57133839564458545904789328652610540031895535786011264182548375833179829124845398393126574488675311145377107878746854204162666250198684504466355949195922066574942592095735778929325357290444962472405416790722118445437122269675520000000000000000000000000000000000000.l,
275 471472363599206132240694321176194377951192623045460204976904578317542573467421580346978030238114995699562728104819596262106947389303901748942909887857509625114880781313585012959529941660203611234871833992565791817698209861793313332044734813700096000000000000000000000000000000000000000.l,
276 7257415615307998967396728211129263114716991681296451376543577798900561843401706157852350749242617459511490991237838520776666022565442753025328900773207510902400430280058295603966612599658257104398558294257568966313439612262571094946806711205568880457193340212661452800000000000000000000000000000000000000000.l,
277};
278
284template <std::floating_point Real>
285inline Real factorialForSmallValue(int n)
286{
287 constexpr int step = 10;
288 int residue = n % step;
289 if (residue <= step / 2) {
291 int nPrev = n - residue;
292 Real fact = static_cast<Real>(FACTORIAL_TABLE[nPrev / step]);
293 for (int i = 1; i <= residue; ++i)
294 fact *= nPrev + i;
295 return fact;
296 }
297
299 int nNext = n - residue + step;
300 Real denominator = 1;
301 for (int i = 0; i < step - residue; ++i)
302 denominator *= nNext - i;
303 return static_cast<Real>(FACTORIAL_TABLE[nNext / step]) / denominator;
304}
305
306template <std::floating_point Real>
307inline Real factorial(Real n)
308{
309 if (n < 0)
310 return 0.0;
311 return (n > FACTORIAL_TABLESIZE) ? std::tgamma(n + 1) : factorialForSmallValue<Real>(static_cast<int>(n));
312}
313
314template <std::integral Integer>
315inline long double factorial(Integer n)
316{
317 return factorial(static_cast<long double>(n));
318}
319
320template <std::floating_point Real>
321inline Real doubleFactorial(int n)
322{
323 if (n < 0)
324 return 0.0;
325 Real result = 1.0;
326 for (int factor = n; factor > 1; factor -= 2)
327 result *= factor;
328 return result;
329}
330
331constexpr long double LOGFACTORIAL_TABLE[] = {0.l,
332 0.l,
333 0.693147180559945286227l,
334 1.79175946922805495731l,
335 3.17805383034794575181l,
336 4.78749174278204581157l,
337 6.57925121201010121297l,
338 8.52516136106541466688l,
339 10.6046029027452508586l,
340 12.8018274800814690906l,
341 15.1044125730755158799l,
342 17.5023078458738865493l,
343 19.9872144956618846834l,
344 22.5521638531234209779l,
345 25.1912211827386798291l,
346 27.8992713838408903371l,
347 30.6718601060806754788l,
348 33.5050734501368836504l,
349 36.3954452080330526087l,
350 39.3398841871994946473l,
351 42.3356164607534850575l,
352 45.3801388984769147328l,
353 48.4711813518352201413l,
354 51.6066755677643769218l,
355 54.7847293981123115714l,
356 58.0036052229805250136l,
357 61.2617017610020013763l,
358 64.5575386270063233951l,
359 67.889743137181540078l,
360 71.2570389671680146648l,
361 74.6582363488301581356l,
362 78.0922235533153070719l,
363 81.5579594561150429399l,
364 85.0544670175815298307l,
365 88.5808275421976674124l,
366 92.1361756036870929165l,
367 95.719694542143201943l,
368 99.3306124547874276232l,
369 102.96819861451380973l,
370 106.63176026064346047l,
371 110.320639714757405159l,
372 114.034211781461706892l,
373 117.77188139974506953l,
374 121.533081515438624365l,
375 125.317271149356884052l,
376 129.123933639127216111l,
377 132.952575035616291643l,
378 136.802722637326382937l,
379 140.67392364823425055l,
380 144.565743946344895221l,
381 148.47776695177302031l,
382 152.409592584497374901l,
383 156.360836303078798437l,
384 160.331128216630901306l,
385 164.320112263195198921l,
386 168.32744544842768164l,
387 172.352797139162817075l,
388 176.395848406997345137l,
389 180.456291417543781108l,
390 184.533828861449478609l,
391 188.628173423671597675l,
392 192.739047287844897483l,
393 196.866181672889979382l,
394 201.009316399281516397l,
395 205.168199482641171016l,
396 209.342586752536817585l,
397 213.532241494563237438l,
398 217.73693411395424846l,
399 221.956441819130361637l,
400 226.190548323727625757l,
401 230.439043565776955802l,
402 234.701723442818263266l,
403 238.978389561834319466l,
404 243.268849002982733509l,
405 247.57291409618684952l,
406 251.890402209723191618l,
407 256.221135550009535109l,
408 260.564940971863222785l,
409 264.921649798552778066l,
410 269.29109765101981111l,
411 273.673124285693745605l,
412 278.067573440366174964l,
413 282.474292687630452292l,
414 286.893133295426991936l,
415 291.32395009427034438l,
416 295.766601350760595324l,
417 300.220948647014154176l,
418 304.686856765668721891l,
419 309.164193580146900331l,
420 313.652829949879048854l,
421 318.152639620209299665l,
422 322.663499126726208033l,
423 327.18528770377525916l,
424 331.717887196928472804l,
425 336.261181979198454428l,
426 340.815058870799020951l,
427 345.379407062266864159l,
428 349.954118040770254083l,
429 354.539085519440789085l,
430 359.13420536957545437l,
431 363.739375555563526632l,
432 368.354496072404685947l,
433 372.979468885689016133l,
434 377.614197873918612913l,
435 382.258588773060012045l,
436 386.912549123217502256l,
437 391.575988217329552299l,
438 396.248817051791547783l,
439 400.930948278915707306l,
440 405.622296161144902271l,
441 410.322776526937275321l,
442 415.032306728249579919l,
443 419.750805599544776214l,
444 424.478193418257092162l,
445 429.214391866651567398l,
446 433.959323995014813136l,
447 438.712914186121224702l,
448 443.475088120918996992l,
449 448.245772745384613245l,
450 453.024896238496125989l,
451 457.81238798127822065l,
452 462.608178526874951331l,
453 467.412199571608141468l,
454 472.224383926980635806l,
455 477.0446654925856933l,
456 481.87297922988790333l,
457 486.70926113683941594l,
458 491.553448223297948516l,
459 496.405478487217578731l,
460 501.265290891579240906l,
461 506.132825342034834648l,
462 511.008022665236012472l,
463 515.890824587822407921l,
464 520.781173716044122557l,
465 525.67901351599505233l,
466 530.584288294433576993l,
467 535.496943180169523657l,
468 540.41692410599773666l,
469 545.344177791154834267l,
470 550.278651724285509772l,
471 555.220294146894843834l,
472 560.169054037273099311l,
473 565.124881094874240262l,
474 570.087725725134191634l,
475 575.057539024710195008l,
476 580.034272767130801185l,
477 585.017879388839105559l,
478 590.008311975617743883l,
479 595.005524249382006019l,
480 600.009470555327425245l,
481 605.020105849423657673l,
482 610.037385686238621929l,
483 615.061266207084827329l,
484 620.091704128477431368l,
485 625.128656730890952531l,
486 630.172081847810204636l,
487 635.221937855059650246l,
488 640.278183660408103606l,
489 645.34077869343502698l,
490 650.409682895655237189l,
491 655.484856710889062015l,
492 660.566261075873512709l,
493 665.653857411105946085l,
494 670.747607611912712855l,
495 675.847474039736766827l,
496 680.953419513637527416l,
497 686.065407301994014233l,
498 691.183401114410685295l,
499 696.307365093814041757l,
500 701.437263808737156978l,
501 706.573062245787468783l,
502 711.714725802289990497l,
503 716.862220279103553366l,
504 722.015511873601212756l,
505 727.174567172815841332l,
506 732.339353146739199474l,
507 737.509837141777325087l,
508 742.685986874351215192l,
509 747.867770424643367733l,
510 753.055156230484158186l,
511 758.248113081374413014l,
512 763.446610112640087209l,
513 768.650616799716999594l,
514 773.860102952558349898l,
515 779.075038710167291356l,
516 784.295394535245691259l,
517 789.521141208958852076l,
518 794.752249825813464668l,
519 799.988691788643450309l,
520 805.230438803703009398l,
521 810.477462875863579939l,
522 815.729736303910158313l,
523 820.987231675938005537l,
524 826.249921864842804098l,
525 831.517780023906198039l,
526 836.790779582469781417l,
527 842.068894241700377279l,
528 847.352097970438421726l,
529 852.640365001132863654l,
530 857.933669825857350588l,
531 863.231987192405426867l,
532 868.535292100464630494l,
533 873.843559797865736982l,
534 879.156765776907491272l,
535 884.474885770751825476l,
536 889.797895749890244588l,
537 895.125771918679674855l,
538 900.458490711945046314l,
539 905.796028791646449463l,
540 911.138363043611207104l,
541 916.485470574328701332l,
542 921.837328707804772421l,
543 927.193914982476826481l,
544 932.555207148186127597l,
545 937.921183163207956568l,
546 943.291821191335770891l,
547 948.667099599019934431l,
548 954.046996952560334648l,
549 959.431492015349476787l,
550 964.820563745166055014l,
551 970.214191291518204707l,
552 975.612353993035981148l,
553 981.01503137490840345l,
554 986.422203146368360649l,
555 991.833849198223447274l,
556 997.249949600427953555l,
557 1002.67048459970021668l,
558 1008.0954346171816951l,
559 1013.52478024613617436l,
560 1018.9585022496903548l,
561 1024.39658155861343403l,
562 1029.83899926913522904l,
563 1035.28573664080158778l,
564 1040.73677509436743094l,
565 1046.19209620972492303l,
566 1051.65168172386916012l,
567 1057.11551352889478039l,
568 1062.58357367002986393l,
569 1068.05584434370143754l,
570 1073.53230789563281178l,
571 1079.01294681897479677l,
572 1084.49774375246556701l,
573 1089.98668147862213118l,
574 1095.47974292196272472l,
575 1100.97691114725603256l,
576 1106.47816935780065251l,
577 1111.9835008937332077l,
578 1117.49288923036101551l,
579 1123.00631797652590649l,
580 1128.52377087299055347l,
581 1134.04523179085276752l,
582 1139.57068472998480502l,
583 1145.10011381749609427l,
584 1150.63350330622370166l,
585 1156.17083757324212456l,
586 1161.71210111840059653l};
587
588constexpr size_t LOGFACTORIAL_TABLESIZE = 255;
589
590template <std::floating_point Real>
591inline Real lfact(size_t n)
592{
593 return (n > LOGFACTORIAL_TABLESIZE) ? std::lgamma(static_cast<Real>(n + 1))
594 : static_cast<Real>(LOGFACTORIAL_TABLE[n]);
595}
596
597template <std::floating_point Real>
598inline Real stirlingError(size_t n)
599{
600 if (n >= 10)
601 return lgammaCorrection(static_cast<Real>(n));
602 long double x = n;
603 return lfact(n) - (x + 0.5l) * std::log(x) + x - 0.5l * std::log(2 * randlib::consts::pi);
604}
605
606template <std::floating_point Real>
607inline Real logFactorialRatio(size_t x, size_t y)
608{
609 if (x <= LOGFACTORIAL_TABLESIZE && y <= LOGFACTORIAL_TABLESIZE)
610 return lfact<Real>(x) - lfact<Real>(y);
612 const Real difference = static_cast<Real>(x) - static_cast<Real>(y), yp1 = static_cast<Real>(y) + 1;
613 const Real logQuotient =
614 (2 * std::fabs(difference) < yp1) ? std::log1p(difference / yp1) : std::log((static_cast<Real>(x) + 1) / yp1);
615 return (yp1 - Real{0.5}) * logQuotient + difference * (std::log(static_cast<Real>(x) + 1) - 1) +
616 stirlingError<Real>(x + 1) - stirlingError<Real>(y + 1);
617}
618
619template <std::floating_point Real>
620inline Real bd0(std::type_identity_t<Real> x, Real m)
621{
622 if (std::fabs(x - m) < 0.1 * (x + m)) {
623 Real v = (x - m) / (x + m);
624 Real sum = (x - m) * v;
625 Real term = 2 * x * v;
626 Real vSq = v * v;
627 for (int j = 1; j < 1000; ++j) {
628 term *= vSq;
629 Real next = sum + term / (2 * j + 1);
630 if (next == sum)
631 return sum;
632 sum = next;
633 }
634 }
635 return x * std::log(x / m) + m - x;
636}
637
638template <std::floating_point Real>
639inline Real ldfact(size_t n)
640{
641 if (n & 1)
642 return lfact<Real>(n) - ldfact<Real>(n - 1);
643 size_t k = n >> 1;
644 return static_cast<Real>(k) * static_cast<Real>(randlib::consts::ln2) + lfact<Real>(k);
645}
646
647template <std::floating_point Real>
648inline Real binom(size_t n, size_t k)
649{
651 if (k == 0 || k == n)
652 return 1.0;
653 if (k == 1 || k == n - 1)
654 return n;
655 if (k == 2 || k == n - 2)
656 return 0.5 * n * (n - 1);
658 Real lfactn = lfact<Real>(n);
659 Real lfactk = lfact<Real>(k);
660 Real lfactnmk = lfact<Real>(n - k);
661 return std::exp(lfactn - lfactk - lfactnmk);
662}
663
668template <std::floating_point Real>
669inline Real digammamLogForLargeX(Real x)
670{
672 static constexpr long double taylorCoef[] = {0.00833333333333333333l, -0.00396825396825396825l,
673 0.00416666666666666666l, -0.00757575757575757576l,
674 0.02115090296908478727l, -0.08333333333333333333l};
675 static constexpr int bounds[] = {6000, 320, 75, 30, 20, 10, 7};
676
678 int degree = 0;
679 while (x < bounds[degree])
680 ++degree;
681
683 long double firstTerm = taylorCoef[5] / x - 0.5;
684 long double y = firstTerm / x;
685 for (int i = 0; i < degree; ++i)
686 y += taylorCoef[i] / std::pow(x, 2 * i + 4);
687 return y;
688}
689
690template <std::floating_point Real>
691inline Real digamma(Real x)
692{
694 if (x < 0.0) {
695 Real y = 1.0 - x;
696 Real z = randlib::consts::pi / std::tan(randlib::consts::pi * y);
697 return digamma(y) + z;
698 }
701 long double sum = 0;
702 while (x < 7.0) {
703 sum -= 1.0 / x;
704 ++x;
705 }
706 return sum + digammamLogForLargeX(x) + std::log(x);
707}
708
709template <std::floating_point Real>
710inline Real digammamLog(Real x)
711{
712 if (x < 0.0)
713 throw std::invalid_argument("Argument x of (digamma - log) function should be non-negative");
714 if (x == 0.0)
715 return -INFINITY;
716 return (x < 7.0) ? digamma(x) - std::log(x) : digammamLogForLargeX(x);
717}
718
719template <std::floating_point Real>
720inline Real trigamma(Real x)
721{
723 if (x < 0.0) {
724 Real z = randlib::consts::pi / std::sin(randlib::consts::pi * x);
725 return z * z - trigamma(1.0 - x);
726 }
728 static constexpr long double taylorCoef[] = {0.16666666666666666666l, -0.03333333333333333333l,
729 0.02380952380952380952l, -0.03333333333333333333l,
730 0.07575757575757575757l, -0.25311355311355311355l};
731 static constexpr int bounds[] = {1000, 140, 47, 25, 15, 10};
732
735 long double y = 0;
736 while (x < 10.0) {
737 y += 1.0 / (x * x);
738 ++x;
739 }
740
742 int degree = 1;
743 while (x < bounds[degree - 1])
744 ++degree;
745
747 long double firstTerm = 1.0 / x;
748 firstTerm += 0.5 / (x * x);
749 firstTerm += taylorCoef[0] / std::pow(x, 3);
750 y += firstTerm;
751 for (int i = 1; i < degree; ++i)
752 y += taylorCoef[i] / std::pow(x, 2 * i + 3);
753 return y;
754}
755
760template <std::floating_point Real>
761inline Real log1pmx(Real t)
762{
763 if (std::fabs(t) > 0.5)
764 return std::log1p(t) - t;
765 Real power = t, sum = 0.0;
766 for (int k = 2; k < 1000; ++k) {
767 power *= -t;
768 Real term = power / k;
769 sum += term;
770 if (std::fabs(term) <= 1e-17 * std::fabs(sum))
771 break;
772 }
773 return sum;
774}
775
781template <std::floating_point Real>
782inline Real logPrefactor(Real a, Real x, Real logX)
783{
784 if (a < 10.0)
785 return a * logX - x - std::lgamma(a);
786 Real t = (x - a) / a;
787 Real core = (std::fabs(t) <= 0.5) ? a * log1pmx(t) : a * (logX - std::log(a)) - (x - a);
788 return core + 0.5 * std::log(a / (2 * randlib::consts::pi)) - lgammaCorrection(a);
789}
790
794template <std::floating_point Real>
795inline Real lowerSeries(Real a, Real x)
796{
797 Real term = 1.0, sum = 1.0;
798 for (int n = 1; n < 1000000; ++n) {
799 term *= x / (a + n);
800 sum += term;
801 if (term <= 1e-17 * sum)
802 break;
803 }
804 return sum;
805}
806
811template <std::floating_point Real>
812inline Real upperContinuedFraction(Real a, Real x)
813{
814 constexpr Real tiny = std::same_as<Real, double>
815 ? Real{1e-300}
816 : std::numeric_limits<Real>::min() / std::numeric_limits<Real>::epsilon();
817 Real b = x + 1.0 - a, c = 1.0 / tiny, d = 1.0 / b, h = d;
818 for (int i = 1; i < 1000000; ++i) {
819 Real an = -i * (i - a);
820 b += 2.0;
821 d = an * d + b;
822 if (std::fabs(d) < tiny)
823 d = tiny;
824 c = b + an / c;
825 if (std::fabs(c) < tiny)
826 c = tiny;
827 d = 1.0 / d;
828 Real del = d * c;
829 h *= del;
830 if (std::fabs(del - 1.0) <= 1e-16)
831 break;
832 }
833 return h;
834}
835
836template <std::floating_point Real>
837struct IncompleteGamma
838{
839 Real p, q, logP, logQ;
840};
841
847template <std::floating_point Real>
848inline IncompleteGamma<Real> incompleteGamma(Real a, Real x, Real logX)
849{
850 if (x == 0.0)
851 return {0.0, 1.0, -INFINITY, 0.0};
852 if (a == 0.0 || std::isinf(x))
853 return {1.0, 0.0, 0.0, -INFINITY};
854 Real logPre = logPrefactor(a, x, logX);
855 if (x < a + 1.0) {
856 Real logP = logPre - std::log(a) + std::log(lowerSeries(a, x));
857 return {std::exp(logP), -std::expm1(logP), logP, math::log1mexp(logP)};
858 }
859 Real logQ = logPre + std::log(upperContinuedFraction(a, x));
860 return {-std::expm1(logQ), std::exp(logQ), math::log1mexp(logQ), logQ};
861}
862
863template <std::floating_point Real>
864inline void parametersVerification(Real a, Real x)
865{
866 if (!(x >= 0.0))
867 throw std::invalid_argument("Argument x of gamma function should be non-negative, but it's equal to " +
868 std::to_string(x));
869 if (!(a >= 0.0))
870 throw std::invalid_argument("Parameter a of gamma function should be non-negative, but it's equal to " +
871 std::to_string(a));
872}
873
874template <std::floating_point Real>
875inline Real lpgamma(std::type_identity_t<Real> a, Real x, std::type_identity_t<Real>, std::type_identity_t<Real>)
876{
877 parametersVerification(a, x);
878 return incompleteGamma(a, x, std::log(x)).logP;
879}
880
881template <std::floating_point Real>
882inline Real lpgamma(std::type_identity_t<Real> a, Real x, std::type_identity_t<Real> logX)
883{
884 parametersVerification(a, x);
885 return incompleteGamma(a, x, logX).logP;
886}
887
888template <std::floating_point Real>
889inline Real lpgamma(std::type_identity_t<Real> a, Real x)
890{
891 parametersVerification(a, x);
892 return incompleteGamma(a, x, std::log(x)).logP;
893}
894
895template <std::floating_point Real>
896inline Real pgamma(std::type_identity_t<Real> a, Real x, std::type_identity_t<Real>, std::type_identity_t<Real>)
897{
898 parametersVerification(a, x);
899 return incompleteGamma(a, x, std::log(x)).p;
900}
901
902template <std::floating_point Real>
903inline Real pgamma(std::type_identity_t<Real> a, Real x, std::type_identity_t<Real> logX)
904{
905 parametersVerification(a, x);
906 return incompleteGamma(a, x, logX).p;
907}
908
909template <std::floating_point Real>
910inline Real pgamma(std::type_identity_t<Real> a, Real x)
911{
912 parametersVerification(a, x);
913 return incompleteGamma(a, x, std::log(x)).p;
914}
915
916template <std::floating_point Real>
917inline Real lqgamma(std::type_identity_t<Real> a, Real x, std::type_identity_t<Real>, std::type_identity_t<Real>)
918{
919 parametersVerification(a, x);
920 return incompleteGamma(a, x, std::log(x)).logQ;
921}
922
923template <std::floating_point Real>
924inline Real lqgamma(std::type_identity_t<Real> a, Real x, std::type_identity_t<Real> logX)
925{
926 parametersVerification(a, x);
927 return incompleteGamma(a, x, logX).logQ;
928}
929
930template <std::floating_point Real>
931inline Real lqgamma(std::type_identity_t<Real> a, Real x)
932{
933 parametersVerification(a, x);
934 return incompleteGamma(a, x, std::log(x)).logQ;
935}
936
937template <std::floating_point Real>
938inline Real qgamma(std::type_identity_t<Real> a, Real x, std::type_identity_t<Real>, std::type_identity_t<Real>)
939{
940 parametersVerification(a, x);
941 return incompleteGamma(a, x, std::log(x)).q;
942}
943
944template <std::floating_point Real>
945inline Real qgamma(std::type_identity_t<Real> a, Real x, std::type_identity_t<Real> logX)
946{
947 parametersVerification(a, x);
948 return incompleteGamma(a, x, logX).q;
949}
950
951template <std::floating_point Real>
952inline Real qgamma(std::type_identity_t<Real> a, Real x)
953{
954 parametersVerification(a, x);
955 return incompleteGamma(a, x, std::log(x)).q;
956}
957
958} // namespace math
959
960} // namespace randlib
randlib::consts::ln2
constexpr long double ln2
Definition constants.hpp:18
randlib::consts::pi
constexpr long double pi
Definition constants.hpp:22
randlib::math::parametersVerification
void parametersVerification(Real a, Real x)
Definition math/gamma.hpp:864
randlib::math::qgamma
Real qgamma(std::type_identity_t< Real > a, Real x, std::type_identity_t< Real > logA, std::type_identity_t< Real > lgammaA)
Calculate upper regularized incomplete gamma function, accelerated by using precalculated values of l...
Definition math/gamma.hpp:938
randlib::math::factorial
Real factorial(Real n)
Calculate n! using table values for small n <= 255 and built-in gamma function for large n > 255.
Definition math/gamma.hpp:307
randlib::math::logPrefactor
Real logPrefactor(Real a, Real x, Real logX)
Definition math/gamma.hpp:782
randlib::math::digamma
Real digamma(Real x)
Definition math/gamma.hpp:691
randlib::math::upperContinuedFraction
Real upperContinuedFraction(Real a, Real x)
Modified Lentz evaluation of the continued fraction for Q(a, x), x ≥ a + 1.
Definition math/gamma.hpp:812
randlib::math::lfact
Real lfact(size_t n)
Calculate log(n!) using table values for small n <= 255 and built-in lgamma function for large n > 25...
Definition math/gamma.hpp:591
randlib::math::LOGFACTORIAL_TABLE
constexpr long double LOGFACTORIAL_TABLE[]
Definition math/gamma.hpp:331
randlib::math::ldfact
Real ldfact(size_t n)
ldoublefact
Definition math/gamma.hpp:639
randlib::math::log1mexp
Real log1mexp(Real x)
Definition functions.hpp:260
randlib::math::lpgamma
Real lpgamma(std::type_identity_t< Real > a, Real x, std::type_identity_t< Real > logA, std::type_identity_t< Real > lgammaA)
Calculate logarithm of lower incomplete gamma function, accelerated by using precalculated value of l...
Definition math/gamma.hpp:875
randlib::math::binom
Real binom(size_t n, size_t k)
Calculate binomial coefficient.
Definition math/gamma.hpp:648
randlib::math::bd0
Real bd0(std::type_identity_t< Real > x, Real m)
Deviance term of the saddle-point expansions of Poisson and binomial probabilities (C.
Definition math/gamma.hpp:620
randlib::math::lgammaCorrection
Real lgammaCorrection(Real x)
Definition special_functions.hpp:350
randlib::math::factorialForSmallValue
Real factorialForSmallValue(int n)
Get n! using product method with table.
Definition math/gamma.hpp:285
randlib::math::stirlingError
Real stirlingError(size_t n)
Definition math/gamma.hpp:598
randlib::math::doubleFactorial
Real doubleFactorial(int n)
Calculate n!!
Definition math/gamma.hpp:321
randlib::math::incompleteGamma
IncompleteGamma< Real > incompleteGamma(Real a, Real x, Real logX)
Regularized incomplete gamma ratios P(a, x) and Q(a, x) = 1 - P(a, x) and their logarithms.
Definition math/gamma.hpp:848
randlib::math::log1pmx
Real log1pmx(Real t)
Definition math/gamma.hpp:761
randlib::math::trigamma
Real trigamma(Real x)
Definition math/gamma.hpp:720
randlib::math::lowerSeries
Real lowerSeries(Real a, Real x)
Definition math/gamma.hpp:795
randlib::math::LOGFACTORIAL_TABLESIZE
constexpr size_t LOGFACTORIAL_TABLESIZE
Definition math/gamma.hpp:588
randlib::math::lqgamma
Real lqgamma(std::type_identity_t< Real > a, Real x, std::type_identity_t< Real > logA, std::type_identity_t< Real > lgammaA)
Calculate logarithm of upper incomplete gamma function, accelerated by using precalculated values of ...
Definition math/gamma.hpp:917
randlib::math::FACTORIAL_TABLESIZE
constexpr int FACTORIAL_TABLESIZE
Definition math/gamma.hpp:254
randlib::math::FACTORIAL_TABLE
constexpr long double FACTORIAL_TABLE[]
Definition math/gamma.hpp:258
randlib::math::pgamma
Real pgamma(std::type_identity_t< Real > a, Real x, std::type_identity_t< Real > logA, std::type_identity_t< Real > lgammaA)
Calculate lower regularized incomplete gamma function, accelerated by using precalculated value of lo...
Definition math/gamma.hpp:896
randlib::math::digammamLogForLargeX
Real digammamLogForLargeX(Real x)
Definition math/gamma.hpp:669
randlib::math::logFactorialRatio
Real logFactorialRatio(size_t x, size_t y)
Definition math/gamma.hpp:607
randlib::math::digammamLog
Real digammamLog(Real x)
Definition math/gamma.hpp:710
randlib
Definition bivariate_normal.hpp:10
special_functions.hpp
randlib::math::IncompleteGamma
Definition math/gamma.hpp:838
randlib::math::IncompleteGamma::logQ
Real logQ
Definition math/gamma.hpp:839
randlib::math::IncompleteGamma::logP
Real logP
Definition math/gamma.hpp:839
randlib::math::IncompleteGamma::q
Real q
Definition math/gamma.hpp:839
randlib::math::IncompleteGamma::p
Real p
Definition math/gamma.hpp:839