RandLib 2.0.0
Probability distributions, fitting, and heavy-tail tools for C++20
Loading...
Searching...
No Matches
functions.hpp
Go to the documentation of this file.
1#pragma once
2
3#include <array>
4#include <cmath>
5#include <concepts>
6#include <functional>
7#include <complex>
8#include <limits>
9#include <vector>
10#include <utility>
11#include <type_traits>
12#include <algorithm>
13#include <numeric>
14
15#include <randlib/math/constants.hpp>
16
17namespace randlib
18{
19
20namespace math
21{
22
23template <int TABLE_SIZE>
24constexpr std::array<LongDoublePair, TABLE_SIZE>
25createZiggurat(const std::function<LongDoublePair(const LongDoublePair &)> &funPtr, LongDoublePair first,
26 LongDoublePair second);
27
28namespace detail
29{
30
37template <typename RealType>
38bool areClose(RealType a, RealType b, RealType eps = 1e-6)
39{
40 if (a == b)
41 return true;
42 RealType fa = std::fabs(a);
43 RealType fb = std::fabs(b);
44 return std::fabs(b - a) < eps * std::max(fa, fb);
45}
46
47} // namespace detail
48
53template <std::floating_point Real>
54int sign(Real x);
55
60template <std::floating_point Real>
61Real atan(Real x);
62
67template <std::floating_point Real>
68Real softplus(Real x);
69
74template <std::floating_point Real>
75Real log1mexp(Real x);
76
81template <std::floating_point Real>
82Real logexpm1l(Real x);
83
88template <std::floating_point Real>
89Real log2mexp(Real x);
90
95template <std::floating_point Real>
96Real erfinv(Real p);
97
102template <std::floating_point Real>
103Real erfcinv(Real p);
104
109template <std::floating_point Real>
110Real erfcx(Real x);
111
117template <std::floating_point Real>
118Real xexpxsqerfc(Real x);
119
125template <std::floating_point Real>
126Real harmonicNumber(Real exponent, int number);
127
134template <std::floating_point Real>
135Real logBesselI(std::type_identity_t<Real> nu, Real x);
136
143template <std::floating_point Real>
144Real logBesselK(std::type_identity_t<Real> nu, Real x);
145
151template <std::floating_point Real>
152Real W0Lambert(Real x, std::type_identity_t<Real> epsilon = Real{1e-11});
153
158template <std::floating_point Real>
159std::complex<Real> W0Lambert(std::complex<Real> z);
160
166template <std::floating_point Real>
167Real Wm1Lambert(Real x, std::type_identity_t<Real> epsilon = Real{1e-11});
168
179template <std::floating_point Real>
180Real MarcumP(std::type_identity_t<Real> mu, Real x, std::type_identity_t<Real> y, std::type_identity_t<Real> sqrtX,
181 std::type_identity_t<Real> sqrtY, std::type_identity_t<Real> logX, std::type_identity_t<Real> logY);
182
189template <std::floating_point Real>
190Real MarcumP(std::type_identity_t<Real> mu, Real x, std::type_identity_t<Real> y);
191
202template <std::floating_point Real>
203Real MarcumQ(std::type_identity_t<Real> mu, Real x, std::type_identity_t<Real> y, std::type_identity_t<Real> sqrtX,
204 std::type_identity_t<Real> sqrtY, std::type_identity_t<Real> logX, std::type_identity_t<Real> logY);
205
212template <std::floating_point Real>
213Real MarcumQ(std::type_identity_t<Real> mu, Real x, std::type_identity_t<Real> y);
214} // namespace math
215
216} // namespace randlib
217
218#include <randlib/math/gamma.hpp>
219#include <randlib/math/special_functions.hpp>
220#include <randlib/math/beta.hpp>
221#include <randlib/math/quadrature.hpp>
222#include <randlib/math/roots.hpp>
223#include <randlib/math/optimize.hpp>
224
225namespace randlib
226{
227
228namespace math
229{
230
231template <std::floating_point Real>
232inline int sign(Real x)
233{
234 return (x > 0) ? 1 : ((x < 0) ? -1 : 0);
235}
236
237template <std::floating_point Real>
238inline Real atan(Real x)
239{
244 if (x == 0.0)
245 return 0.0;
246 if (x > 1.0)
247 return randlib::consts::half_pi - std::atan(1.0 / x);
248 return (x < -1.0) ? -randlib::consts::half_pi - std::atan(1.0 / x) : std::atan(x);
249}
250
251template <std::floating_point Real>
252inline Real softplus(Real x)
253{
254 if (x < 20.0)
255 return std::log1p(std::exp(x));
256 return (x < 35.0) ? x + std::exp(-x) : x;
257}
258
259template <std::floating_point Real>
260inline Real log1mexp(Real x)
261{
262 return (x < -randlib::consts::ln2) ? std::log1p(-std::exp(x)) : std::log(-std::expm1(x));
263}
264
265template <std::floating_point Real>
266inline Real logexpm1l(Real x)
267{
268 if (x < 20.0)
269 return std::log(std::expm1(x));
270 return (x < 35.0) ? x - std::exp(-x) : x;
271}
272
273template <std::floating_point Real>
274inline Real log2mexp(Real x)
275{
276 return std::log1p(-std::expm1(x));
277}
278
279namespace detail
280{
281
282template <std::floating_point Real>
283inline Real erfinvChebyshevSeries(Real x, long double t, const long double *array, int size)
284{
286 long double Tn = t, Tnm1 = 1, sum = 0.0;
287 for (int i = 1; i != size; ++i) {
288 sum += array[i] * Tn;
290 long double temp = Tn;
291 Tn *= 2 * t;
292 Tn -= Tnm1;
293 Tnm1 = temp;
294 }
295 return static_cast<Real>(static_cast<long double>(x) * (array[0] + sum));
296}
297
298template <std::floating_point Real>
299inline Real erfinvAux1(Real beta)
300{
302 static constexpr long double D5 = -9.199992358830151031278420l, D6 = 2.794990820124599493768426l;
303 static constexpr int muSize = 26;
304 static constexpr long double mu[muSize] = {.9885750640661893136460358l,
305 .0108577051845994776160281l,
306 -.0017511651027627952594825l,
307 .0000211969932065633437984l,
308 .0000156648714042435087911l,
309 -.05190416869103124261e-5l,
310 -.00371357897426717780e-5l,
311 .00012174308662357429e-5l,
312 -.00001768115526613442e-5l,
313 -.119372182556161e-10l,
314 .003802505358299e-10l,
315 -.000660188322362e-10l,
316 -.000087917055170e-10l,
317 -.3506869329e-15l,
318 -.0697221497e-15l,
319 -.0109567941e-15l,
320 -.0011536390e-15l,
321 -.0000263938e-15l,
322 .05341e-20l,
323 -.22610e-20l,
324 .09552e-20l,
325 -.05250e-20l,
326 .02487e-20l,
327 -.01134e-20l,
328 .00420e-20l};
329 return erfinvChebyshevSeries(beta, D5 / std::sqrt(beta) + D6, mu, muSize);
330}
331
332template <std::floating_point Real>
333inline Real erfinvAux2(Real beta)
334{
336 static constexpr long double D3 = -0.5594576313298323225436913l, D4 = 2.287915716263357638965891l;
337 static constexpr int deltaSize = 38;
338 static constexpr long double delta[deltaSize] = {.9566797090204925274526373l,
339 -.0231070043090649036999908l,
340 -.0043742360975084077333218l,
341 -.0005765034226511854809364l,
342 -.0000109610223070923931242l,
343 .0000251085470246442787982l,
344 .0000105623360679477511955l,
345 .27544123300306391503e-5l,
346 .04324844983283380689e-5l,
347 -.00205303366552086916e-5l,
348 -.00438915366654316784e-5l,
349 -.00176840095080881795e-5l,
350 -.00039912890280463420e-5l,
351 -.00001869324124559212e-5l,
352 .00002729227396746077e-5l,
353 .00001328172131565497e-5l,
354 .318342484482286e-10l,
355 .016700607751926e-10l,
356 -.020364649611537e-10l,
357 -.009648468127965e-10l,
358 -.002195672778128e-10l,
359 -.000095689813014e-10l,
360 .000137032572230e-10l,
361 .000062538505417e-10l,
362 .000014584615266e-10l,
363 .1078123993e-15l,
364 -.0709229988e-15l,
365 -.0391411775e-15l,
366 -.0111659209e-15l,
367 -.0015770366e-15l,
368 .0002853149e-15l,
369 .0002716662e-15l,
370 .0000176835e-15l,
371 .09828e-20l,
372 .20464e-20l,
373 .08020e-20l,
374 .01650e-20l};
375 return erfinvChebyshevSeries(beta, D3 * beta + D4, delta, deltaSize);
376}
377
378template <std::floating_point Real>
379inline Real erfinvAux3(Real beta)
380{
382 static constexpr long double D1 = -1.548813042373261659512742l, D2 = 2.565490123147816151928163l;
383 static constexpr int lambdaSize = 27;
384 static constexpr long double lambda[lambdaSize] = {.9121588034175537733059200l,
385 -.0162662818676636958546661l,
386 .0004335564729494453650589l,
387 .0002144385700744592065205l,
388 .26257510757648130176e-5l,
389 -.30210910501037969912e-5l,
390 -.00124060618367572157e-5l,
391 .00624066092999917380e-5l,
392 -.00005401247900957858e-5l,
393 -.00014232078975315910e-5l,
394 .343840281955305e-10l,
395 .335848703900138e-10l,
396 -.014584288516512e-10l,
397 -.008102174258833e-10l,
398 .000525324085874e-10l,
399 .000197115408612e-10l,
400 -.000017494333828e-10l,
401 -.4800596619e-15l,
402 .0557302987e-15l,
403 .0116326054e-15l,
404 -.0017262489e-15l,
405 -.0002784973e-15l,
406 .0000524481e-15l,
407 .65270e-20l,
408 -.15707e-20l,
409 -.01475e-20l,
410 .00450e-20l};
411 return erfinvChebyshevSeries(beta, D1 * beta + D2, lambda, lambdaSize);
412}
413
414template <std::floating_point Real>
415inline Real erfinvAux4(Real p)
416{
418 static constexpr int xiSize = 39;
419 static constexpr long double xi[xiSize] = {.9928853766189408231495800l,
420 .1204675161431044864647846l,
421 .0160781993420999447267039l,
422 .0026867044371623158279591l,
423 .0004996347302357262947170l,
424 .0000988982185991204409911l,
425 .0000203918127639944337340l,
426 .43272716177354218758e-5l,
427 .09380814128593406758e-5l,
428 .02067347208683427411e-5l,
429 .00461596991054300078e-5l,
430 .00104166797027146217e-5l,
431 .00023715009995921222e-5l,
432 .00005439284068471390e-5l,
433 .00001255489864097987e-5l,
434 .291381803663201e-10l,
435 .067949421808797e-10l,
436 .015912343331569e-10l,
437 .003740250585245e-10l,
438 .000882087762421e-10l,
439 .000208650897725e-10l,
440 .000049488041039e-10l,
441 .000011766394740e-10l,
442 .2803855725e-15l,
443 .0669506638e-15l,
444 .0160165495e-15l,
445 .0038382583e-15l,
446 .0009212851e-15l,
447 .0002214615e-15l,
448 .0000533091e-15l,
449 .0000128488e-15l,
450 .31006e-20l,
451 .07491e-20l,
452 .01812e-20l,
453 .00439e-20l,
454 .00106e-20l,
455 .00026e-20l,
456 .00006e-20l,
457 .00002e-20l};
458 return erfinvChebyshevSeries(p, p * p / 0.32 - 1.0, xi, xiSize);
459}
460
461} // namespace detail
462
463template <std::floating_point Real>
464inline Real erfinv(Real p)
465{
467 if (p < 0.0)
468 return -erfinv(-p);
469 if (p > 1.0)
470 throw std::invalid_argument("Argument p should be in interval [-1, 1]");
471 if (p == 1.0)
472 return INFINITY;
473 if (p == 0.0)
474 return 0.0;
475 if (p < 0.8)
476 return detail::erfinvAux4(p);
478 Real beta = std::sqrt(-std::log1p(-p * p));
479 if (p < 0.9975)
480 return detail::erfinvAux3(beta);
481 return (1.0 - p < 5e-16) ? detail::erfinvAux1(beta) : detail::erfinvAux2(beta);
482}
483
484template <std::floating_point Real>
485inline Real erfcinv(Real p)
486{
488 if (p > 1.0)
489 return -erfcinv(2.0 - p);
490 if (p == 0.0)
491 return INFINITY;
492 if (p == 1.0)
493 return 0.0;
494 if (p > 0.2)
495 return detail::erfinvAux4(1.0 - p);
496 Real pSq = p * p, p2 = 2 * p;
497 Real beta = std::sqrt(-std::log(p2 - pSq));
498 if (p > 0.0025)
499 return detail::erfinvAux3(beta);
500 return (p > 5e-16) ? detail::erfinvAux2(beta) : detail::erfinvAux1(beta);
501}
502
503template <std::floating_point Real>
504inline Real xexpxsqerfc(Real x)
505{
506 static constexpr int MAX_X = 10;
507 static constexpr int N = 10;
508 if (x < MAX_X) {
509 Real y = x * x;
510 y += std::log(std::erfc(x));
511 return x * std::exp(y);
512 }
513 Real log2xSq = randlib::consts::ln2 + 2 * std::log(x);
514 Real sum = 0.0;
515 for (int n = 1; n != N; ++n) {
516 Real add = math::ldfact(2 * n - 1);
517 add -= n * log2xSq;
518 add = std::exp(add);
519 sum += (n & 1) ? -add : add;
520 }
521 return (1.0 + sum) / randlib::consts::sqrtpi;
522}
523
524template <std::floating_point Real>
525inline Real erfcx(Real x)
526{
527 if (x < 1.0)
528 return std::exp(x * x) * std::erfc(x);
529 return xexpxsqerfc(x) / x;
530}
531
532template <std::floating_point Real>
533inline Real harmonicNumber(Real exponent, int number)
534{
535 if (number < 1)
536 return 0;
537 if (exponent == 1)
538 return static_cast<Real>(randlib::consts::egamma) + digamma(static_cast<Real>(number + 1));
539 if (exponent == 2)
540 return static_cast<Real>(randlib::consts::pi_sq) / 6 - trigamma(static_cast<Real>(number + 1));
541 Real res = 1.0;
542 for (int i = 2; i <= number; ++i)
543 res += std::pow(i, -exponent);
544 return res;
545}
546
547namespace detail
548{
549
556template <std::floating_point Real>
557inline Real logHankelSeries(Real nu, Real x, int sign)
558{
559 const Real mu = 4 * nu * nu;
560 Real term = 1, sum = 1;
561 for (int k = 1; k < 64; ++k) {
562 const Real odd = 2 * k - 1;
563 const Real next = term * sign * (mu - odd * odd) / (8 * k * x);
564 if (!(std::fabs(next) < std::fabs(term)))
565 break;
566 term = next;
567 sum += term;
568 if (std::fabs(term) <= std::numeric_limits<Real>::epsilon() * std::fabs(sum))
569 break;
570 }
571 return std::log(sum);
572}
573
574} // namespace detail
575
576template <std::floating_point Real>
577inline Real logBesselI(std::type_identity_t<Real> nu, Real x)
578{
579 if (x < 0) {
580 Real roundNu = std::round(nu);
581 bool nuIsInt = detail::areClose(nu, roundNu);
582 if (nuIsInt) {
583 int nuInt = roundNu;
584 return (nuInt % 2) ? NAN : logBesselI(nu, -x);
585 }
586 return -INFINITY;
587 }
588
589 if (x == 0) {
590 if (nu == 0)
591 return 0.0;
592 Real roundNu = std::round(nu);
593 bool nuIsInt = detail::areClose(nu, roundNu);
594 return (nu > 0 || nuIsInt) ? -INFINITY : INFINITY;
595 }
596
597 if (std::fabs(nu) == 0.5) {
599 Real y = x - 0.5 * (randlib::consts::ln2 + randlib::consts::lnpi + std::log(x));
600 y += (nu > 0) ? math::log1mexp(-2 * x) : math::softplus(-2 * x);
601 return y;
602 }
603
604 if (nu < 0) {
606 Real besseli = math::besselI(-nu, x);
607 Real sinPiAbsNu = std::sin(-randlib::consts::pi * nu);
608 Real y = 0;
609 if (sinPiAbsNu == 0 || detail::areClose(nu, std::round(nu)))
610 y = besseli;
611 else {
612 Real besselk = math::besselK(-nu, x);
613 y = besseli + randlib::consts::two_inv_pi * sinPiAbsNu * besselk;
614 }
615 return (y <= 0) ? -INFINITY : std::log(y);
616 }
617
618 Real besseli = math::besselI(nu, x);
619 if (std::isfinite(besseli) && !(besseli > Real{1e-250})) {
621 const Real quarterSquare = Real{0.25} * x * x;
622 Real term = 1, sum = 1;
623 for (int j = 1; j < 100000; ++j) {
624 term *= quarterSquare / (j * (nu + j));
625 sum += term;
626 if (term <= std::numeric_limits<Real>::epsilon() * sum)
627 break;
628 }
629 return nu * std::log(Real{0.5} * x) - std::lgamma(nu + 1) + std::log(sum);
630 }
631 return std::isfinite(besseli) ? std::log(besseli)
632 : x - 0.5 * (randlib::consts::ln2 + randlib::consts::lnpi + std::log(x)) +
633 detail::logHankelSeries<Real>(nu, x, -1);
634}
635
636template <std::floating_point Real>
637inline Real logBesselK(std::type_identity_t<Real> nu, Real x)
638{
640 nu = std::fabs(nu);
641
642 if (x == 0.0)
643 return INFINITY;
644
645 Real besselk = 0;
646 if (nu == 0.5 || (besselk = math::besselK(nu, x)) == 0)
647 return 0.5 * (randlib::consts::lnpi - randlib::consts::ln2 - std::log(x)) - x +
648 detail::logHankelSeries<Real>(nu, x, 1);
649
650 if (!std::isfinite(besselk))
651 return (nu == 0) ? std::log(-std::log(x))
652 : std::lgamma(nu) - static_cast<Real>(randlib::consts::ln2) - nu * std::log(Real{0.5} * x);
653
654 return std::log(besselk);
655}
656
663namespace detail
664{
665
666template <std::floating_point Real>
667inline Real WLambert(Real x, Real w0, Real epsilon)
668{
669 Real w = w0;
670 Real step = 0;
671 do {
672 Real ew = std::exp(w);
673 Real wew = w * ew;
674 Real numerator1 = wew - x;
675 Real wp1 = w + 1;
676 Real denominator1 = ew * wp1;
677 Real numerator2 = (w + 2) * numerator1;
678 Real denominator2 = 2 * wp1;
679 step = numerator2 / denominator2;
680 step = numerator1 / (denominator1 - step);
681 w -= step;
682 } while (std::fabs(step) > epsilon);
683 return w;
684}
685
686} // namespace detail
687
688template <std::floating_point Real>
689inline Real W0Lambert(Real x, std::type_identity_t<Real> epsilon)
690{
691 Real w = 0;
692 if (x < -randlib::consts::inv_e)
693 throw std::invalid_argument("Argument x should be greater than -1/e, but it's equal to " + std::to_string(x));
694 if (x > 10) {
695 Real logX = std::log(x);
696 Real loglogX = std::log(logX);
697 w = logX - loglogX;
698 }
699 return detail::WLambert(x, w, epsilon);
700}
701
702template <std::floating_point Real>
703inline std::complex<Real> W0Lambert(std::complex<Real> z)
704{
705 using Complex = std::complex<Real>;
706 if (z == Complex(0))
707 return 0;
709 const Complex logShifted = std::log(Real(1) + z);
710 Complex w = logShifted * (Real(1) - std::log(Real(1) + logShifted) / (Real(2) + logShifted));
711 for (int iteration = 0; iteration < 64; ++iteration) {
712 const Complex exponential = std::exp(w), residual = w * exponential - z, wp1 = w + Real(1);
713 const Complex step = residual / (exponential * wp1 - (w + Real(2)) * residual / (Real(2) * wp1));
714 w -= step;
715 if (std::abs(step) <= 4 * std::numeric_limits<Real>::epsilon() * std::abs(w))
716 break;
717 }
718 return w;
719}
720
721template <std::floating_point Real>
722inline Real Wm1Lambert(Real x, std::type_identity_t<Real> epsilon)
723{
724 Real w = -2;
725 if (x < -randlib::consts::inv_e || x > 0)
726 throw std::invalid_argument(
727 "Argument x should be greater than -1/e and smaller or equal to 0, but it's equal to " + std::to_string(x));
728 if (x > -0.1) {
729 Real logmX = std::log(-x);
730 Real logmlogmX = std::log(-logmX);
731 w = logmX - logmlogmX;
732 }
733 return detail::WLambert(x, w, epsilon);
734}
735
736namespace detail
737{
738
747template <std::floating_point Real>
748inline Real MarcumPSeries(Real mu, Real x, Real y, Real logX, Real logY)
749{
751 static constexpr Real ln2piEps = -35.0;
752 Real lgammamu = std::lgamma(mu);
753 Real C = lgammamu - ln2piEps + mu;
754
757 Real root = std::max(Real{0.5} * (mu * mu + 4 * x * y - mu), Real{1});
758 Real logXY = logX + logY;
759 if (!math::findRootNewtonSecondOrder<Real>(
760 [C, mu, logXY](Real n) {
761 Real npmu = n + mu;
762 Real logn = std::log(n), lognpmu = std::log(npmu);
763 Real first = logn - 2 - logXY;
764 first *= n;
765 first += npmu * lognpmu;
766 first -= C;
767 Real second = logn + lognpmu - logXY;
768 Real third = 1.0 / n + 1.0 / npmu;
769 return Triplet<Real>(first, second, third);
770 },
771 root))
773 throw std::runtime_error("Marcum P function: failure in numerical procedure");
774
776 Real sum = 0.0;
778 int n0 = static_cast<int>(std::max(std::ceil(root), Real{5}));
779 Real mpn0 = mu + n0;
780 Real P = pgamma(mpn0, y, logY);
781 Real diffP = (mpn0 - 1) * logY - y - std::lgamma(mpn0);
782 diffP = std::exp(diffP);
783 for (int n = n0; n > 0; --n) {
784 Real term = n * logX - x - lfact(n);
785 Real mupnm1 = mu + n - 1;
786 term = std::exp(term) * P;
787 sum += term;
788 if (n % 5 == 0) {
791 P = pgamma(mupnm1, y, logY);
792 diffP = (mupnm1 - 1) * logY - y - std::lgamma(mupnm1);
793 diffP = std::exp(diffP);
794 } else {
796 P += diffP;
797 diffP *= mupnm1 / y;
798 }
799 }
801 Real lastTerm = lpgamma(mu, y, logY) - x;
802 lastTerm = std::exp(lastTerm);
803 sum += lastTerm;
804 return sum;
805}
806
807template <std::floating_point Real>
808inline Real MarcumQSeries(Real mu, Real x, Real y, Real logX, Real logY)
809{
810 if (x <= 0)
811 return qgamma(mu, y, logY);
812 const Real saddle = Real{0.5} * (std::sqrt(mu * mu + 4 * x * y) - mu);
813 const int center = static_cast<int>(std::floor(std::max(x, saddle)));
814 Real shape = mu + center;
815 Real weight = std::exp(-x + center * logX - std::lgamma(center + 1.0));
816 Real q = qgamma(shape, y, logY);
817 Real delta = std::exp(shape * logY - y - std::lgamma(shape + 1));
818 Real sum = weight * q;
819
820 Real lowerWeight = weight;
821 Real previousTerm = sum;
822 int decreasing = 0;
823 for (int n = center; n > 0; --n) {
824 lowerWeight *= n / x;
825 const Real lowerQ = qgamma(mu + n - 1, y, logY);
826 const Real term = lowerWeight * lowerQ;
827 sum += term;
828 decreasing = (term <= previousTerm) ? decreasing + 1 : 0;
829 if (decreasing >= 8 && term <= 4 * std::numeric_limits<Real>::epsilon() * sum)
830 break;
831 previousTerm = term;
832 }
833
834 decreasing = 0;
835 previousTerm = weight * q;
836 for (int n = center + 1; n < 1000000; ++n) {
837 weight *= x / n;
838 q += delta;
839 shape += 1;
840 delta *= y / shape;
841 const Real term = weight * std::min(Real{1}, q);
842 sum += term;
843 decreasing = (term <= previousTerm) ? decreasing + 1 : 0;
844 if (decreasing >= 8 && term <= 4 * std::numeric_limits<Real>::epsilon() * sum)
845 break;
846 previousTerm = term;
847 }
848 return std::clamp(sum, Real{0}, Real{1});
849}
850
853template <std::floating_point Real>
854inline Real MarcumPSaddleSeries(Real mu, Real x, Real y, Real logX, Real logY)
855{
856 if (x <= 0)
857 return pgamma(mu, y, logY);
858 const Real saddle = Real{0.5} * (std::sqrt(mu * mu + 4 * x * y) - mu);
859 const int center = static_cast<int>(std::floor(std::max(x, saddle)));
860 const Real shape = mu + center;
861 const Real weight = std::exp(-x + center * logX - std::lgamma(center + 1.0));
862 const Real p = pgamma(shape, y, logY);
863 Real sum = weight * p;
864
865 Real lowerWeight = weight, lowerP = p;
866 Real delta = std::exp((shape - 1) * logY - y - std::lgamma(shape));
867 Real previousTerm = sum;
868 int decreasing = 0;
869 for (int n = center; n > 0; --n) {
870 lowerWeight *= n / x;
871 lowerP += delta;
872 delta *= (mu + n - 1) / y;
873 const Real term = lowerWeight * std::min(Real{1}, lowerP);
874 sum += term;
875 decreasing = (term <= previousTerm) ? decreasing + 1 : 0;
876 if (decreasing >= 8 && term <= 4 * std::numeric_limits<Real>::epsilon() * sum)
877 break;
878 previousTerm = term;
879 }
880
881 Real upperWeight = weight;
882 previousTerm = weight * p;
883 decreasing = 0;
884 for (int n = center + 1; n < 1000000; ++n) {
885 upperWeight *= x / n;
886 const Real term = upperWeight * pgamma(mu + n, y, logY);
887 sum += term;
888 decreasing = (term <= previousTerm) ? decreasing + 1 : 0;
889 if (decreasing >= 8 && term <= 4 * std::numeric_limits<Real>::epsilon() * sum)
890 break;
891 previousTerm = term;
892 }
893 return std::clamp(sum, Real{0}, Real{1});
894}
895
904template <std::floating_point Real>
905inline Real MarcumPAsymptoticForLargeXY(Real mu, Real x, Real y, Real sqrtX, Real sqrtY)
906{
907 Real xi = 2 * sqrtX * sqrtY;
908 Real sigma = (y + x) / xi - 1;
909 Real sigmaXi = sigma * xi;
910 Real rho = sqrtY / sqrtX;
911 Real aux = std::erfc(sqrtX - sqrtY);
912 Real Phi = (sigma == 0) ? Real{0} : sign(x - y) * std::sqrt(static_cast<Real>(randlib::consts::pi) / sigma) * aux;
913 Real Psi0 = aux / std::sqrt(rho);
914 Real logXi = randlib::consts::ln2 + 0.5 * std::log(x * y);
916 int n0 = static_cast<int>(std::max(std::ceil(sigmaXi), Real{7}));
917 Real sum = 0.0;
918 Real A1 = 1.0, A2 = 1.0;
919 for (int n = 1; n <= n0; ++n) {
921 Real nmHalf = n - 0.5;
922 Phi *= -sigma;
923 Real temp = -sigmaXi - nmHalf * logXi;
924 Phi += std::exp(temp);
925 Phi /= nmHalf;
927 Real coef1 = (nmHalf - mu) * (nmHalf + mu);
928 Real coef2 = (nmHalf - mu + 1) * (nmHalf + mu - 1);
929 Real denom = -2 * n;
930 A1 *= coef1 / denom;
931 A2 *= coef2 / denom;
933 Real Psi = Phi * (A2 - A1 / rho);
934 sum += (n & 1) ? Psi : -Psi;
935 }
936 sum /= randlib::consts::sqrt2pi;
937 sum += Psi0;
938 sum *= 0.5 * std::pow(rho, mu);
939 return sum;
940}
941
950template <std::floating_point Real>
951inline Real MarcumPForMuSmallerThanOne(Real mu, Real x, Real y, Real logX, Real logY)
952{
953 // TODO: check Krishnamoorthy paper for alternative representation
954
959 Real aux = x + mu * static_cast<Real>(randlib::consts::ln2) + std::lgamma(mu);
960 Real log2x = randlib::consts::ln2 + logX;
961 Real I = (randlib::consts::ln2 + logY) * mu - aux;
962 I = std::exp(I) / mu;
963
964 Real mum1 = mu - 1.0;
965 I += math::integrate(
966 [x, log2x, mum1, aux](Real t) {
967 if (t <= 0)
968 return Real{0};
970 Real log2T = std::log(2 * t);
971 Real exponent = mum1 * log2T;
972 exponent -= aux;
973
975 Real logBessel = math::logBesselI(mum1, 2 * std::sqrt(x * t));
976 Real z = mum1 * (log2T - log2x);
977 Real log2F = 0.5 * z - t - x + logBessel;
978
980 return std::exp(log2F) - 2 * std::exp(exponent);
981 },
982 0, y);
983 return I;
984}
985
986template <std::floating_point Real>
987inline Real MarcumQIntegrand(Real theta, Real xi, Real sqrt1pXiSq, Real mu, Real y)
988{
989 Real sinTheta = std::sin(theta);
990 Real theta_sinTheta = theta / sinTheta;
991 Real rho = std::hypot(theta_sinTheta, xi);
992 Real theta_sinThetapRho = theta_sinTheta + rho;
993 Real r = 0.5 * theta_sinThetapRho / y;
994 Real cosTheta = std::cos(theta);
995 Real psi = cosTheta * rho - sqrt1pXiSq;
996 Real frac = theta_sinThetapRho / (1.0 + sqrt1pXiSq);
997 psi -= std::log(frac);
998 Real numerator = (sinTheta - theta * cosTheta) / (sinTheta * rho);
999 numerator += cosTheta - r;
1000 numerator *= r;
1001 Real denominator = r * r - 2.0 * r * cosTheta + 1.0;
1002 Real f = numerator / denominator;
1003 return std::exp(mu * psi) * f;
1004}
1005
1006template <std::floating_point Real>
1007inline Real MarcumQIntergralRepresentation(Real mu, Real x, Real y, Real sqrtX, Real sqrtY)
1008{
1009 Real xi = 2.0 * sqrtX * sqrtY / mu;
1010 Real sqrt1pXiSq = std::sqrt(1.0 + xi * xi);
1011 Real yPrime = y / mu;
1012 Real s0 = 0.5 * (1.0 + sqrt1pXiSq) / yPrime;
1013 Real phi = x / s0 + y * s0 - std::log(s0) * mu;
1014 auto integrand = [=](Real theta) { return MarcumQIntegrand(theta, xi, sqrt1pXiSq, mu, yPrime); };
1015 Real integral =
1016 math::integrate(integrand, -static_cast<Real>(randlib::consts::pi), static_cast<Real>(randlib::consts::pi));
1017 return 0.5 * std::exp(-x - y + phi) / randlib::consts::pi * integral;
1018}
1019
1020} // namespace detail
1021
1022template <std::floating_point Real>
1023inline Real MarcumP(std::type_identity_t<Real> mu, Real x, std::type_identity_t<Real> y,
1024 std::type_identity_t<Real> sqrtX, std::type_identity_t<Real> sqrtY, std::type_identity_t<Real> logX,
1025 std::type_identity_t<Real> logY)
1026{
1027 if (x < 0.0 || y <= 0.0)
1028 return 0.0;
1029 if (std::isinf(y))
1030 return 1.0;
1031 if (std::isinf(x))
1032 return 0.0;
1033
1034 if (x < 30)
1035 return detail::MarcumPSeries(mu, x, y, logX, logY);
1036
1038 const Real saddle = Real{0.5} * (std::hypot(mu, 2 * sqrtX * sqrtY) - mu);
1039 if (std::isfinite(saddle) && std::max(x, saddle) < std::numeric_limits<int>::max() - 1)
1040 return (y < mu + x) ? detail::MarcumPSaddleSeries(mu, x, y, logX, logY)
1041 : 1 - detail::MarcumQSeries(mu, x, y, logX, logY);
1042
1043 Real xi = 2 * sqrtX * sqrtY;
1044 if (xi > 30 && mu * mu < 2 * xi)
1045 return detail::MarcumPAsymptoticForLargeXY(mu, x, y, sqrtX, sqrtY);
1046
1047 Real temp = std::sqrt(4 * x + 2 * mu);
1048 Real f1 = x + mu - temp, f2 = x + mu + temp;
1049 if (y > f1 && y < f2)
1050 // IF mu > 135
1051 return 1.0 - detail::MarcumQIntergralRepresentation(mu, x, y, sqrtX, sqrtY);
1052
1053 // TODO: implement the rest techniques
1054
1055 Real mum1 = mu - 1;
1056 return math::integrate(
1057 [mum1, logX, x](Real t) {
1058 if (t < 0.0)
1059 return Real{0};
1060 if (t == 0.0)
1061 return (mum1 == 0) ? Real{0.5} * std::exp(-x) : Real{0};
1062 Real logBesseli = math::logBesselI(mum1, 2 * std::sqrt(x * t));
1063 Real z = 0.5 * mum1 * (std::log(t) - logX);
1064 Real h = t + x;
1065 return std::exp(logBesseli + z - h);
1066 },
1067 0, y);
1068}
1069
1070template <std::floating_point Real>
1071inline Real MarcumP(std::type_identity_t<Real> mu, Real x, std::type_identity_t<Real> y)
1072{
1073 Real sqrtX = std::sqrt(x), sqrtY = std::sqrt(y);
1074 Real logX = std::log(x), logY = std::log(y);
1075 return MarcumP(mu, x, y, sqrtX, sqrtY, logX, logY);
1076}
1077
1078template <std::floating_point Real>
1079inline Real MarcumQ(std::type_identity_t<Real> mu, Real x, std::type_identity_t<Real> y,
1080 std::type_identity_t<Real> sqrtX, std::type_identity_t<Real> sqrtY, std::type_identity_t<Real> logX,
1081 std::type_identity_t<Real> logY)
1082{
1083 if (x < 0 || y <= 0)
1084 return 1;
1085 if (std::isinf(y))
1086 return 0;
1087 if (std::isinf(x))
1088 return 1;
1089 const Real saddle = Real{0.5} * (std::hypot(mu, 2 * sqrtX * sqrtY) - mu);
1090 if (!std::isfinite(saddle) || std::max(x, saddle) >= std::numeric_limits<int>::max() - 1)
1091 return 1 - MarcumP(mu, x, y, sqrtX, sqrtY, logX, logY);
1094 if (x >= 30 && y < mu + x)
1095 return 1 - detail::MarcumPSaddleSeries(mu, x, y, logX, logY);
1096 return detail::MarcumQSeries(mu, x, y, logX, logY);
1097}
1098
1099template <std::floating_point Real>
1100inline Real MarcumQ(std::type_identity_t<Real> mu, Real x, std::type_identity_t<Real> y)
1101{
1102 Real sqrtX = std::sqrt(x), sqrtY = std::sqrt(y);
1103 Real logX = std::log(x), logY = std::log(y);
1104 return MarcumQ(mu, x, y, sqrtX, sqrtY, logX, logY);
1105}
1106
1107} // namespace math
1108
1109} // namespace randlib
randlib::consts::ln2
constexpr long double ln2
Definition constants.hpp:18
randlib::consts::egamma
constexpr long double egamma
Definition constants.hpp:36
randlib::consts::half_pi
constexpr long double half_pi
Definition constants.hpp:23
randlib::consts::sqrtpi
constexpr long double sqrtpi
Definition constants.hpp:31
randlib::consts::pi
constexpr long double pi
Definition constants.hpp:22
randlib::consts::two_inv_pi
constexpr long double two_inv_pi
Definition constants.hpp:25
randlib::consts::pi_sq
constexpr long double pi_sq
Definition constants.hpp:39
randlib::consts::sqrt2pi
constexpr long double sqrt2pi
Definition constants.hpp:32
randlib::consts::inv_e
constexpr long double inv_e
Definition constants.hpp:35
randlib::consts::lnpi
constexpr long double lnpi
Definition constants.hpp:21
randlib::detail::standardNormal
Real standardNormal(G &g)
Ziggurat algorithm of Marsaglia and Tsang with 256 layers and an explicit sign bit.
Definition variates.hpp:45
randlib::math::detail::MarcumQIntergralRepresentation
Real MarcumQIntergralRepresentation(Real mu, Real x, Real y, Real sqrtX, Real sqrtY)
Definition functions.hpp:1007
randlib::math::detail::erfinvChebyshevSeries
Real erfinvChebyshevSeries(Real x, long double t, const long double *array, int size)
Definition functions.hpp:283
randlib::math::detail::WLambert
Real WLambert(Real x, Real w0, Real epsilon)
Definition functions.hpp:667
randlib::math::detail::areClose
bool areClose(RealType a, RealType b, RealType eps=1e-6)
Definition functions.hpp:38
randlib::math::detail::erfinvAux1
Real erfinvAux1(Real beta)
Definition functions.hpp:299
randlib::math::detail::MarcumQIntegrand
Real MarcumQIntegrand(Real theta, Real xi, Real sqrt1pXiSq, Real mu, Real y)
Definition functions.hpp:987
randlib::math::detail::erfinvAux4
Real erfinvAux4(Real p)
Definition functions.hpp:415
randlib::math::detail::MarcumPSeries
Real MarcumPSeries(Real mu, Real x, Real y, Real logX, Real logY)
Definition functions.hpp:748
randlib::math::detail::erfinvAux2
Real erfinvAux2(Real beta)
Definition functions.hpp:333
randlib::math::detail::logHankelSeries
Real logHankelSeries(Real nu, Real x, int sign)
Definition functions.hpp:557
randlib::math::detail::MarcumQSeries
Real MarcumQSeries(Real mu, Real x, Real y, Real logX, Real logY)
Definition functions.hpp:808
randlib::math::detail::erfinvAux3
Real erfinvAux3(Real beta)
Definition functions.hpp:379
randlib::math::detail::MarcumPSaddleSeries
Real MarcumPSaddleSeries(Real mu, Real x, Real y, Real logX, Real logY)
Σ e^{-x} x^n / n! P(μ + n, y), summed outwards from the dominant Poisson term.
Definition functions.hpp:854
randlib::math::detail::MarcumPAsymptoticForLargeXY
Real MarcumPAsymptoticForLargeXY(Real mu, Real x, Real y, Real sqrtX, Real sqrtY)
Definition functions.hpp:905
randlib::math::detail::MarcumPForMuSmallerThanOne
Real MarcumPForMuSmallerThanOne(Real mu, Real x, Real y, Real logX, Real logY)
Definition functions.hpp:951
randlib::math::besselI
Real besselI(std::type_identity_t< Real > nu, Real x)
Definition special_functions.hpp:320
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::erfinv
Real erfinv(Real p)
Definition functions.hpp:464
randlib::math::log2mexp
Real log2mexp(Real x)
Definition functions.hpp:274
randlib::math::digamma
Real digamma(Real x)
Definition math/gamma.hpp:691
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::createZiggurat
constexpr std::array< LongDoublePair, TABLE_SIZE > createZiggurat(const std::function< LongDoublePair(const LongDoublePair &)> &funPtr, LongDoublePair first, LongDoublePair second)
randlib::math::ldfact
Real ldfact(size_t n)
ldoublefact
Definition math/gamma.hpp:639
randlib::math::sign
int sign(Real x)
Definition functions.hpp:232
randlib::math::log1mexp
Real log1mexp(Real x)
Definition functions.hpp:260
randlib::math::erfcx
Real erfcx(Real x)
Definition functions.hpp:525
randlib::math::harmonicNumber
Real harmonicNumber(Real exponent, int number)
Definition functions.hpp:533
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::W0Lambert
Real W0Lambert(Real x, std::type_identity_t< Real > epsilon=Real{1e-11})
Definition functions.hpp:689
randlib::math::atan
Real atan(Real x)
Definition functions.hpp:238
randlib::math::Wm1Lambert
Real Wm1Lambert(Real x, std::type_identity_t< Real > epsilon=Real{1e-11})
Definition functions.hpp:722
randlib::math::softplus
Real softplus(Real x)
Definition functions.hpp:252
randlib::math::logBesselI
Real logBesselI(std::type_identity_t< Real > nu, Real x)
Calculates logarithm of modified Bessel function of the 1st kind.
Definition functions.hpp:577
randlib::math::trigamma
Real trigamma(Real x)
Definition math/gamma.hpp:720
randlib::math::MarcumP
Real MarcumP(std::type_identity_t< Real > mu, Real x, std::type_identity_t< Real > y, std::type_identity_t< Real > sqrtX, std::type_identity_t< Real > sqrtY, std::type_identity_t< Real > logX, std::type_identity_t< Real > logY)
Definition functions.hpp:1023
randlib::math::MarcumQ
Real MarcumQ(std::type_identity_t< Real > mu, Real x, std::type_identity_t< Real > y, std::type_identity_t< Real > sqrtX, std::type_identity_t< Real > sqrtY, std::type_identity_t< Real > logX, std::type_identity_t< Real > logY)
Definition functions.hpp:1079
randlib::math::xexpxsqerfc
Real xexpxsqerfc(Real x)
xexpxsqerfc
Definition functions.hpp:504
randlib::math::logBesselK
Real logBesselK(std::type_identity_t< Real > nu, Real x)
Calculates logarithm of modified Bessel function of the 2nd kind.
Definition functions.hpp:637
randlib::math::integrate
Real integrate(const F &f, std::type_identity_t< Real > a, std::type_identity_t< Real > b, std::type_identity_t< Real > tolerance=1e-11, int maxSubdivisions=512)
Global adaptive Gauss-Kronrod quadrature (QUADPACK QAG).
Definition quadrature.hpp:226
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::logexpm1l
Real logexpm1l(Real x)
Definition functions.hpp:266
randlib::math::erfcinv
Real erfcinv(Real p)
Definition functions.hpp:485
randlib::math::besselK
Real besselK(std::type_identity_t< Real > nu, Real x)
Definition special_functions.hpp:330
randlib
Definition bivariate_normal.hpp:10
randlib::Triplet
std::tuple< T, T, T > Triplet
Definition fwd.hpp:36
randlib::LongDoublePair
Pair< long double > LongDoublePair
Definition fwd.hpp:40
special_functions.hpp