Beta-related functions. More...

Namespaces
namespace	detail

namespace	portable

Classes
struct	IncompleteGamma

struct	MinimizeResult
	The MinimizeResult struct Outcome of a multivariate minimization. More...

Functions
template<std::floating_point Real>
Real	logBeta (Real a, std::type_identity_t< Real > b)
	Calculate logarithm of beta function.

template<typename Left , typename Right > requires (std::is_arithmetic_v<Left> && std::is_arithmetic_v<Right> && (std::integral<Left> \|\| std::integral<Right> \|\| !std::same_as<Left, Right>))
auto	logBeta (Left a, Right b) -> std::common_type_t< double, Left, Right >

template<std::floating_point Real>
Real	ibeta (Real x, std::type_identity_t< Real > a, std::type_identity_t< Real > b, std::type_identity_t< Real > logBetaFun, std::type_identity_t< Real > logX, std::type_identity_t< Real > log1mX)
	Fast calculation of regularized beta function, using precalculated values.

template<std::floating_point Real>
Real	ibeta (Real x, std::type_identity_t< Real > a, std::type_identity_t< Real > b)
	Calculate regularized beta function.

template<std::floating_point Real>
Real	ibetaPowerSeries1 (Real x, Real a, Real b, Real logBetaFun, Real logX, Real log1mX)

template<std::floating_point Real>
Real	ibetacontinuedFraction1 (Real x, Real a, Real b, int number, Real logBetaFun, Real logX, Real log1mX)

template<std::floating_point Real>
Real	ibetaPowerSeries2 (Real x, Real a, Real b, Real logBetaFun, Real logX, Real log1mX)

template<std::floating_point Real>
Real	ibetaContinuedFraction2 (Real x, Real a, Real b, int number, Real logBetaFun, Real logX, Real log1mX)

template<std::floating_point Real>
Real	ibetaContinuedFractionLentz (Real x, Real a, Real b, Real logBetaFun, Real logX, Real log1mX)
	Continued fraction of I(x, a, b) (DiDonato and Morris, also Numerical Recipes betacf) by the modified Lentz method; it converges for x below (a + 1) / (a + b + 2) in far fewer terms than b for large shapes.

template<std::floating_point Real>
Real	ibetaLowerTail (Real x, Real a, Real b, Real logBetaFun, Real logX, Real log1mX)
	I(x, a, b) for x not above the mean a / (a + b) and shapes up to IBETA_LARGEST_FINITE_SHAPE.

template<int TABLE_SIZE>
constexpr std::array< LongDoublePair, TABLE_SIZE >	createZiggurat (const std::function< LongDoublePair(const LongDoublePair &)> &funPtr, LongDoublePair first, LongDoublePair second)

template<std::floating_point Real>
int	sign (Real x)

template<std::floating_point Real>
Real	atan (Real x)

template<std::floating_point Real>
Real	softplus (Real x)

template<std::floating_point Real>
Real	log1mexp (Real x)

template<std::floating_point Real>
Real	logexpm1l (Real x)

template<std::floating_point Real>
Real	log2mexp (Real x)

template<std::floating_point Real>
Real	erfinv (Real p)

template<std::floating_point Real>
Real	erfcinv (Real p)

template<std::floating_point Real>
Real	erfcx (Real x)

template<std::floating_point Real>
Real	xexpxsqerfc (Real x)
	xexpxsqerfc

template<std::floating_point Real>
Real	harmonicNumber (Real exponent, int number)

template<std::floating_point Real>
Real	logBesselI (std::type_identity_t< Real > nu, Real x)
	Calculates logarithm of modified Bessel function of the 1st kind.

template<std::floating_point Real>
Real	logBesselK (std::type_identity_t< Real > nu, Real x)
	Calculates logarithm of modified Bessel function of the 2nd kind.

template<std::floating_point Real>
Real	W0Lambert (Real x, std::type_identity_t< Real > epsilon=Real{1e-11})

template<std::floating_point Real>
std::complex< Real >	W0Lambert (std::complex< Real > z)

template<std::floating_point Real>
Real	Wm1Lambert (Real x, std::type_identity_t< Real > epsilon=Real{1e-11})

template<std::floating_point Real>
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)

template<std::floating_point Real>
Real	MarcumP (std::type_identity_t< Real > mu, Real x, std::type_identity_t< Real > y)

template<std::floating_point Real>
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)

template<std::floating_point Real>
Real	MarcumQ (std::type_identity_t< Real > mu, Real x, std::type_identity_t< Real > y)

template<std::floating_point Real>
Real	factorial (Real n)
	Calculate n! using table values for small n <= 255 and built-in gamma function for large n > 255.

template<std::integral Integer>
long double	factorial (Integer n)

template<std::floating_point Real = long double>
Real	doubleFactorial (int n)
	Calculate n!!

template<std::floating_point Real = long double>
Real	lfact (size_t n)
	Calculate log(n!) using table values for small n <= 255 and built-in lgamma function for large n > 255.

template<std::floating_point Real = long double>
Real	ldfact (size_t n)
	ldoublefact

template<std::floating_point Real = long double>
Real	binom (size_t n, size_t k)
	Calculate binomial coefficient.

template<std::floating_point Real = double>
Real	stirlingError (size_t n)

template<std::floating_point Real = double>
Real	logFactorialRatio (size_t x, size_t y)

template<std::floating_point Real>
Real	bd0 (std::type_identity_t< Real > x, Real m)
	Deviance term of the saddle-point expansions of Poisson and binomial probabilities (C.

template<std::floating_point Real>
Real	digamma (Real x)

template<std::floating_point Real>
Real	digammamLog (Real x)

template<std::floating_point Real>
Real	trigamma (Real x)

template<std::floating_point Real>
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 log(a) and log(Γ(a))

template<std::floating_point Real>
Real	lpgamma (std::type_identity_t< Real > a, Real x, std::type_identity_t< Real > logX)
	Calculate logarithm of lower incomplete gamma function, accelerated by using precalculated value of log(x)

template<std::floating_point Real>
Real	lpgamma (std::type_identity_t< Real > a, Real x)
	Calculate logarithm of lower incomplete gamma function.

template<std::floating_point Real>
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 log(a) and log(Γ(a))

template<std::floating_point Real>
Real	pgamma (std::type_identity_t< Real > a, Real x, std::type_identity_t< Real > logX)
	Calculate lower regularized incomplete gamma function, accelerated by using precalculated value of log(x)

template<std::floating_point Real>
Real	pgamma (std::type_identity_t< Real > a, Real x)
	Calculate lower regularized incomplete gamma function.

template<std::floating_point Real>
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 log(a) and log(Γ(a))

template<std::floating_point Real>
Real	lqgamma (std::type_identity_t< Real > a, Real x, std::type_identity_t< Real > logX)
	Calculate logarithm of upper incomplete gamma function, accelerated by using precalculated value of log(x)

template<std::floating_point Real>
Real	lqgamma (std::type_identity_t< Real > a, Real x)
	Calculate logarithm of upper incomplete gamma function.

template<std::floating_point Real>
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 log(a) and log(Γ(a))

template<std::floating_point Real>
Real	qgamma (std::type_identity_t< Real > a, Real x, std::type_identity_t< Real > logX)
	Calculate upper regularized incomplete gamma function, accelerated by using precalculated value of log(x)

template<std::floating_point Real>
Real	qgamma (std::type_identity_t< Real > a, Real x)
	Calculate upper regularized incomplete gamma function.

template<std::floating_point Real>
Real	factorialForSmallValue (int n)
	Get n! using product method with table.

template<std::floating_point Real>
Real	digammamLogForLargeX (Real x)

template<std::floating_point Real>
Real	log1pmx (Real t)

template<std::floating_point Real>
Real	logPrefactor (Real a, Real x, Real logX)

template<std::floating_point Real>
Real	lowerSeries (Real a, Real x)

template<std::floating_point Real>
Real	upperContinuedFraction (Real a, Real x)
	Modified Lentz evaluation of the continued fraction for Q(a, x), x ≥ a + 1.

template<std::floating_point Real>
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.

template<std::floating_point Real>
void	parametersVerification (Real a, Real x)

template<typename RealType >
double	parabolicMinimum (RealType a, RealType b, RealType c, double fa, double fb, double fc)

template<typename RealType >
bool	findMin (const std::function< double(RealType)> &funPtr, const Triplet< RealType > &abc, double &fx, RealType &root, double epsilon)
	Combined Brent's method.

template<typename RealType >
bool	findMin (const std::function< double(RealType)> &funPtr, RealType closePoint, RealType &root, long double epsilon=1e-8)
	Combined Brent's method.

template<std::floating_point Real = double, class F >
MinimizeResult< Real >	nelderMead (const F &f, std::vector< Real > start, Real tolerance=1e-10, int maxIterations=20000)
	Downhill simplex method of Nelder and Mead with the dimension-dependent coefficients of Gao and Han (2012).

template<std::floating_point Real = double, class F >
std::vector< Real >	numericalGradient (const F &f, const std::vector< Real > &x)

template<std::floating_point Real = double, class F >
std::vector< std::vector< Real > >	numericalHessian (const F &f, const std::vector< Real > &x)

template<std::floating_point Real = double, class F >
MinimizeResult< Real >	bfgs (const F &f, std::vector< Real > start, Real tolerance=1e-8, int maxIterations=1000)
	Quasi-Newton minimization of Broyden, Fletcher, Goldfarb and Shanno with central-difference gradients and a backtracking line search under the Armijo condition.

template<std::floating_point Real = double, class F >
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).

template<std::floating_point Real = double, class F >
Real	integrateTanhSinh (const F &f, std::type_identity_t< Real > a, std::type_identity_t< Real > b, std::type_identity_t< Real > tolerance=1e-11, int maxRefinements=12)
	Double-exponential tanh-sinh quadrature on a finite interval.

template<std::floating_point Real = double, class F >
std::complex< Real >	integrateOscillatory (const F &amplitude, std::type_identity_t< Real > frequency, std::type_identity_t< Real > lower=0, std::type_identity_t< Real > tolerance=1e-10, int maxCycles=512)
	QAWF-style Fourier integral on a semi-infinite interval.

template<std::floating_point Real = double, class F , class TailMass > requires std::invocable<const TailMass &, Real>
std::complex< Real >	integrateOscillatory (const F &amplitude, std::type_identity_t< Real > frequency, std::type_identity_t< Real > lower, const TailMass &tailMass, std::type_identity_t< Real > tolerance=1e-10, int maxCycles=512)
	QAWF-style Fourier integral of an amplitude with a known tail bound: the half-period walk stops once the remaining mass is below the tolerance, and is accelerated only beyond half of the mass.

template<std::floating_point Real, class F >
std::pair< Real, Real >	toms748 (const F &f, Real a, Real b, Real fa, Real fb, int maxIterations=200)
	Bracketing root finder of Alefeld, Potra and Shi (ACM TOMS algorithm 748): inverse cubic and Newton-quadratic interpolation steps, a double-length secant step and bisection whenever the bracket does not shrink by half; it converges like Brent's method at worst and faster on smooth functions.

template<std::floating_point Real = double, class F >
Real	findRoot (const F &f, std::type_identity_t< Real > a, std::type_identity_t< Real > b, int maxIterations=200)

template<typename RealType >
bool	findRootNewtonFirstOrder (const std::function< DoublePair(RealType)> &funPtr, RealType &root, long double funTol=1e-10, long double stepTol=1e-6)
	Newton's root-finding procedure, using first derivative.

template<typename RealType >
bool	findRootNewtonSecondOrder (const std::function< DoubleTriplet(RealType)> &funPtr, RealType &root, long double funTol=1e-10, long double stepTol=1e-6)
	Newton's root-finding procedure, using first and second derivatives.

bool	findRootNewtonFirstOrder2d (const std::function< DoublePair(DoublePair)> &funPtr, const std::function< std::tuple< DoublePair, DoublePair >(DoublePair)> &gradPtr, DoublePair &root, long double funTol=1e-10, long double stepTol=1e-6)
	Newton's root-finding procedure for 2 functions of 2 parameters, using Jacobian matrix.

template<std::floating_point Real>
Real	besselI (std::type_identity_t< Real > nu, Real x)

template<std::floating_point Real>
Real	besselK (std::type_identity_t< Real > nu, Real x)

template<std::floating_point Real>
Real	riemannZeta (Real s)

template<std::floating_point Real>
Real	lgammaCorrection (Real x)

Variables
constexpr int	FACTORIAL_TABLESIZE = 175

constexpr long double	FACTORIAL_TABLE []

constexpr long double	LOGFACTORIAL_TABLE []

constexpr size_t	LOGFACTORIAL_TABLESIZE = 255

Detailed Description

Beta-related functions.

Special functions that the C++17 standard library provides only optionally.

Gamma-related functions.

libstdc++ and MSVC implement the C++17 mathematical special functions, libc++ does not. These wrappers forward to std:: when it is available and fall back to RandLib's own implementations otherwise (or always, with RANDLIB_PORTABLE_SPECIAL_FUNCTIONS defined).

Function Documentation

◆ atan()

template<std::floating_point Real>

Real randlib::math::atan ( Real x )

inline

Parameters

x

Returns: safe atan(x)

For small absolute values we use standard technique Otherwise we use relation atan(x) = +/-π/2 - atan(1/x) to avoid numeric problems

◆ bd0()

template<std::floating_point Real>

Real randlib::math::bd0	(	std::type_identity_t< Real >	x,
		Real	m
	)

inline

Deviance term of the saddle-point expansions of Poisson and binomial probabilities (C.

Loader, "Fast and accurate computation of binomial probabilities", 2000)

Parameters

x	non-negative number
m	positive number

Returns: x log(x / m) + m - x, without cancellation for x close to m

◆ besselI()

template<std::floating_point Real>

Real randlib::math::besselI	(	std::type_identity_t< Real >	nu,
		Real	x
	)

inline

Parameters

nu	order, ν ≥ 0
x	argument, x ≥ 0

Returns: modified Bessel function of the first kind I_ν(x)

◆ besselK()

template<std::floating_point Real>

Real randlib::math::besselK	(	std::type_identity_t< Real >	nu,
		Real	x
	)

inline

Parameters

nu	order, ν ≥ 0
x	argument, x > 0

Returns: modified Bessel function of the second kind K_ν(x)

◆ bfgs()

template<std::floating_point Real = double, class F >

MinimizeResult< Real > randlib::math::bfgs	(	const F &	f,
		std::vector< Real >	start,
		Real	tolerance = `1e-8`,
		int	maxIterations = `1000`
	)

Quasi-Newton minimization of Broyden, Fletcher, Goldfarb and Shanno with central-difference gradients and a backtracking line search under the Armijo condition.

Parameters

f	smooth function of a std::span<const Real>
start	starting point
tolerance	the iteration stops when the largest gradient component is below this tolerance times max(1, \|f\|), or when a step no longer changes f
maxIterations	iteration limit

H <- (I - ρ s yᵀ) H (I - ρ y sᵀ) + ρ s sᵀ

◆ binom()

template<std::floating_point Real = long double>

Real randlib::math::binom	(	size_t	n,
		size_t	k
	)

inline

Calculate binomial coefficient.

Parameters

n	non-negative integer
k	non-negative integer, k < n

Returns: C(n, k)

first check trivial cases

next run general procedure

first check trivial cases

next run general procedure

◆ createZiggurat()

template<int TABLE_SIZE>

constexpr std::array< LongDoublePair, TABLE_SIZE > randlib::math::createZiggurat	(	const std::function< LongDoublePair(const LongDoublePair &)> &	funPtr,
		LongDoublePair	first,
		LongDoublePair	second
	)

constexpr

◆ digamma()

template<std::floating_point Real>

Real randlib::math::digamma ( Real x )

inline

Parameters

x	input parameter

Returns: digamma(x) = ln(Γ(x))' = d(Γ(x))/Γ(x)

Negative argument

shift to minimum value, for which series expansion is applicable

Negative argument

shift to minimum value, for which series expansion is applicable

◆ digammamLog()

template<std::floating_point Real>

Real randlib::math::digammamLog ( Real x )

inline

Parameters

x	input parameter

Returns: digamma(x) - log(x)

◆ digammamLogForLargeX()

template<std::floating_point Real>

Real randlib::math::digammamLogForLargeX ( Real x )

inline

Parameters

x

Returns: digamma(x) - log(x) for x > 7

set up tables

choose the amount of terms in series

apply

◆ doubleFactorial()

template<std::floating_point Real = long double>

Real randlib::math::doubleFactorial ( int n )

inline

Calculate n!!

Parameters

n	non-negative integer number

Returns: n!!

◆ erfcinv()

template<std::floating_point Real>

Real randlib::math::erfcinv ( Real p )

inline

Parameters

p

Returns: inverse complementary error function: such x that erfc(x) = p

Consider special cases

◆ erfcx()

template<std::floating_point Real>

Real randlib::math::erfcx ( Real x )

inline

Parameters

x

Returns: scaled complementary error function exp(x^2) erfc(x)

◆ erfinv()

template<std::floating_point Real>

Real randlib::math::erfinv ( Real p )

inline

Parameters

p

Returns: inverse error function: such x that erf(x) = p

Consider special cases

Handle tails

Consider special cases

Handle tails

◆ factorial() [1/2]

template<std::integral Integer>

long double randlib::math::factorial ( Integer n )

inline

◆ factorial() [2/2]

template<std::floating_point Real>

Real randlib::math::factorial ( Real n )

inline

Calculate n! using table values for small n <= 255 and built-in gamma function for large n > 255.

Parameters

n	non-negative integer number

Returns: n!

◆ factorialForSmallValue()

template<std::floating_point Real>

Real randlib::math::factorialForSmallValue ( int n )

inline

Get n! using product method with table.

Parameters

n	non-negative integer number

Returns: n!

go up

go down

◆ findMin() [1/2]

template<typename RealType >

bool randlib::math::findMin	(	const std::function< double(RealType)> &	funPtr,
		const Triplet< RealType > &	abc,
		double &	fx,
		RealType &	root,
		double	epsilon
	)

Combined Brent's method.

Parameters

funPtr
abc	lower boundary / middle / upper boundary
fx	funPtr(b)
root	such x that funPtr(x) is min
epsilon	tolerance

Returns: true if success

The tolerance is relative beyond |x| = 1: an absolute one cannot be met once it is below the spacing of doubles near the minimum, and the iteration would never end

accept u

use golden ratio instead of parabolic approximation

golden ratio [x, c]

golden ratio [a, x]

setting the closest distance between u and x

◆ findMin() [2/2]

template<typename RealType >

bool randlib::math::findMin	(	const std::function< double(RealType)> &	funPtr,
		RealType	closePoint,
		RealType &	root,
		long double	epsilon = `1e-8`
	)

Combined Brent's method.

Parameters

funPtr
closePoint	point that is nearby minimum
root	such x that funPtr(x) is min
epsilon	tolerance

Returns: true if success

we go to the right

parabolic interpolation

go to the left

parabolic interpolation

◆ findRoot()

template<std::floating_point Real = double, class F >

Real randlib::math::findRoot	(	const F &	f,
		std::type_identity_t< Real >	a,
		std::type_identity_t< Real >	b,
		int	maxIterations = `200`
	)

Parameters

f	continuous function with f(a) and f(b) of opposite signs
a	lower end of the bracket
b	upper end of the bracket
maxIterations	maximum number of refinement steps

Returns: root of f in [a, b], to full precision (TOMS 748)

◆ findRootNewtonFirstOrder()

template<typename RealType >

bool randlib::math::findRootNewtonFirstOrder	(	const std::function< DoublePair(RealType)> &	funPtr,
		RealType &	root,
		long double	funTol = `1e-10`,
		long double	stepTol = `1e-6`
	)

Newton's root-finding procedure, using first derivative.

Parameters

funPtr	mapping x \|-> (f(x), f'(x))
root	starting point in input and such x that f(x) = 0 in output
funTol	function tolerance
stepTol	step tolerance

Returns: true if success, false otherwise

Sanity check

Check convergence criteria

◆ findRootNewtonFirstOrder2d()

bool randlib::math::findRootNewtonFirstOrder2d	(	const std::function< DoublePair(DoublePair)> &	funPtr,
		const std::function< std::tuple< DoublePair, DoublePair >(DoublePair)> &	gradPtr,
		DoublePair &	root,
		long double	funTol = `1e-10`,
		long double	stepTol = `1e-6`
	)

inline

Newton's root-finding procedure for 2 functions of 2 parameters, using Jacobian matrix.

Parameters

funPtr	mapping x, y \|-> (f(x, y), g(x, y))
gradPtr	x, y \|-> (f_x(x, y), f_y(x, y), g_x(x, y)), g_y(x, y))
root	starting point in input and such x that f(x) = 0 in output
funTol	function tolerance
stepTol	step tolerance

Returns: true if success, false otherwise

Sanity check

Check convergence criteria

Sanity check

Check convergence criteria

◆ findRootNewtonSecondOrder()

template<typename RealType >

bool randlib::math::findRootNewtonSecondOrder	(	const std::function< DoubleTriplet(RealType)> &	funPtr,
		RealType &	root,
		long double	funTol = `1e-10`,
		long double	stepTol = `1e-6`
	)

Newton's root-finding procedure, using first and second derivatives.

Parameters

funPtr	mapping x \|-> (f(x), f'(x), f''(x))
root	starting point in input and such x that f(x) = 0 in output
funTol	function tolerance
stepTol	step tolerance

Returns: true if success, false otherwise

Sanity check

Check convergence criteria

◆ harmonicNumber()

template<std::floating_point Real>

Real randlib::math::harmonicNumber	(	Real	exponent,
		int	number
	)

inline

Parameters

exponent
number

Returns: sum_{i=1}^{number} i^{-exponent}

◆ ibeta() [1/2]

template<std::floating_point Real>

Real randlib::math::ibeta	(	Real	x,
		std::type_identity_t< Real >	a,
		std::type_identity_t< Real >	b
	)

inline

Calculate regularized beta function.

Parameters

x
a
b

Returns: I(x, a, b) = B(x, a, b) / B(a, b)

Check special and incorrect values of parameters

◆ ibeta() [2/2]

template<std::floating_point Real>

Real randlib::math::ibeta	(	Real	x,
		std::type_identity_t< Real >	a,
		std::type_identity_t< Real >	b,
		std::type_identity_t< Real >	logBetaFun,
		std::type_identity_t< Real >	logX,
		std::type_identity_t< Real >	log1mX
	)

inline

Fast calculation of regularized beta function, using precalculated values.

Parameters

x
a
b
logBetaFun	log(B(a, b))
logX	log(x)
log1mX	log(1-x)

Returns: I(x, a, b) = B(x, a, b) / B(a, b)

Check special and incorrect values of parameters

x may be a rounded complement such as 1 - t for tiny t, so 1 - x is recovered from log(1 - x)

For large shapes the continued fraction converges below c = (a + 1) / (a + b + 2). For a tiny x the complementary fraction at 1 - x would lose its digits, and the direct one still converges up to about four standard deviations above c, so the side whose argument is tiny is preferred there

If x is greater than mean of Beta distribution, convergence of complementary distribution function is faster. Close to 1 the series below converge only at the rate x, so the argument is never moved there: a tiny argument beyond a tinier mean is evaluated directly, and one within this distance of 1 under the mean (for a tiny second shape) through the complement, losing at most ε / I(x, a, b)

Check special and incorrect values of parameters

x may be a rounded complement such as 1 - t for tiny t, so 1 - x is recovered from log(1 - x)

For large shapes the continued fraction converges below c = (a + 1) / (a + b + 2). For a tiny x the complementary fraction at 1 - x would lose its digits, and the direct one still converges up to about four standard deviations above c, so the side whose argument is tiny is preferred there

If x is greater than mean of Beta distribution, convergence of complementary distribution function is faster. Close to 1 the series below converge only at the rate x, so the argument is never moved there: a tiny argument beyond a tinier mean is evaluated directly, and one within this distance of 1 under the mean (for a tiny second shape) through the complement, losing at most ε / I(x, a, b)

◆ ibetacontinuedFraction1()

template<std::floating_point Real>

Real randlib::math::ibetacontinuedFraction1	(	Real	x,
		Real	a,
		Real	b,
		int	number,
		Real	logBetaFun,
		Real	logX,
		Real	log1mX
	)

inline

avoid underflow

◆ ibetaContinuedFraction2()

template<std::floating_point Real>

Real randlib::math::ibetaContinuedFraction2	(	Real	x,
		Real	a,
		Real	b,
		int	number,
		Real	logBetaFun,
		Real	logX,
		Real	log1mX
	)

inline

◆ ibetaContinuedFractionLentz()

template<std::floating_point Real>

Real randlib::math::ibetaContinuedFractionLentz	(	Real	x,
		Real	a,
		Real	b,
		Real	logBetaFun,
		Real	logX,
		Real	log1mX
	)

inline

Continued fraction of I(x, a, b) (DiDonato and Morris, also Numerical Recipes betacf) by the modified Lentz method; it converges for x below (a + 1) / (a + b + 2) in far fewer terms than b for large shapes.

◆ ibetaLowerTail()

template<std::floating_point Real>

Real randlib::math::ibetaLowerTail	(	Real	x,
		Real	a,
		Real	b,
		Real	logBetaFun,
		Real	logX,
		Real	log1mX
	)

inline

I(x, a, b) for x not above the mean a / (a + b) and shapes up to IBETA_LARGEST_FINITE_SHAPE.

We truncate after 0.45 instead of 0.5 here in order to avoid x being too close to 0.5

equivalent continued fraction #1

equivalent continued fraction #2

◆ ibetaPowerSeries1()

template<std::floating_point Real>

Real randlib::math::ibetaPowerSeries1	(	Real	x,
		Real	a,
		Real	b,
		Real	logBetaFun,
		Real	logX,
		Real	log1mX
	)

inline

avoid underflow

◆ ibetaPowerSeries2()

template<std::floating_point Real>

Real randlib::math::ibetaPowerSeries2	(	Real	x,
		Real	a,
		Real	b,
		Real	logBetaFun,
		Real	logX,
		Real	log1mX
	)

inline

◆ incompleteGamma()

template<std::floating_point Real>

IncompleteGamma< Real > randlib::math::incompleteGamma	(	Real	a,
		Real	x,
		Real	logX
	)

inline

Regularized incomplete gamma ratios P(a, x) and Q(a, x) = 1 - P(a, x) and their logarithms.

The smaller ratio is always computed directly (series for x < a + 1, continued fraction otherwise), never by subtraction, so both tails keep their relative accuracy.

◆ integrate()

template<std::floating_point Real = double, class F >

Real randlib::math::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).

Infinite limits are mapped onto (0, 1] by x = a + (1 - t) / t, as in QUADPACK QAGI; the rule never evaluates the endpoints.

Parameters

f	integrand
a	lower limit, possibly -∞
b	upper limit, possibly +∞
tolerance	relative error tolerance
maxSubdivisions	maximal number of bisections

Returns: ∫_a^b f(x) dx

◆ integrateOscillatory() [1/2]

template<std::floating_point Real = double, class F , class TailMass >
requires std::invocable<const TailMass &, Real>

std::complex< Real > randlib::math::integrateOscillatory	(	const F &	amplitude,
		std::type_identity_t< Real >	frequency,
		std::type_identity_t< Real >	lower,
		const TailMass &	tailMass,
		std::type_identity_t< Real >	tolerance = `1e-10`,
		int	maxCycles = `512`
	)

QAWF-style Fourier integral of an amplitude with a known tail bound: the half-period walk stops once the remaining mass is below the tolerance, and is accelerated only beyond half of the mass.

Parameters

amplitude	non-oscillatory real amplitude
frequency	positive angular frequency
lower	finite lower limit
tailMass	bound on the integral of \|amplitude\| from its argument to infinity
tolerance	absolute convergence tolerance
maxCycles	maximum number of half-periods

Returns: integral of amplitude(x) exp(i frequency x) from lower to infinity

◆ integrateOscillatory() [2/2]

template<std::floating_point Real = double, class F >

std::complex< Real > randlib::math::integrateOscillatory	(	const F &	amplitude,
		std::type_identity_t< Real >	frequency,
		std::type_identity_t< Real >	lower = `0`,
		std::type_identity_t< Real >	tolerance = `1e-10`,
		int	maxCycles = `512`
	)

QAWF-style Fourier integral on a semi-infinite interval.

The interval is divided into half-periods and the resulting alternating sequence is accelerated with Wynn's epsilon algorithm.

Parameters

amplitude	non-oscillatory real amplitude
frequency	positive angular frequency
lower	finite lower limit
tolerance	absolute convergence tolerance
maxCycles	maximum number of half-periods

Returns: integral of amplitude(x) exp(i frequency x) from lower to infinity

◆ integrateTanhSinh()

template<std::floating_point Real = double, class F >

Real randlib::math::integrateTanhSinh	(	const F &	f,
		std::type_identity_t< Real >	a,
		std::type_identity_t< Real >	b,
		std::type_identity_t< Real >	tolerance = `1e-11`,
		int	maxRefinements = `12`
	)

Double-exponential tanh-sinh quadrature on a finite interval.

This rule is intended for integrable endpoint singularities. A callable that accepts (x, endpointDistance) receives a negative distance a - x near the left endpoint and a positive distance b - x near the right endpoint, allowing cancellation-free evaluation when x rounds to an endpoint.

Parameters

f	integrand accepting either `f(x)` or `f(x, endpointDistance)`
a	finite lower limit
b	finite upper limit
tolerance	error tolerance relative to the estimated integral of `\|f\|`
maxRefinements	maximum number of step halvings

Returns: ∫_a^b f(x) dx

◆ ldfact()

template<std::floating_point Real = long double>

Real randlib::math::ldfact ( size_t n )

inline

ldoublefact

Parameters

n

Returns: log(n!!)

◆ lfact()

template<std::floating_point Real = long double>

Real randlib::math::lfact ( size_t n )

inline

Calculate log(n!) using table values for small n <= 255 and built-in lgamma function for large n > 255.

Parameters

n	non-negative integer number

Returns: log(n!)

◆ lgammaCorrection()

template<std::floating_point Real>

Real randlib::math::lgammaCorrection ( Real x )

inline

Parameters

x	argument, x ≥ 10

Returns: log Γ(x) - ((x - 1/2) log x - x + log √(2π)), the remainder of Stirling's formula

◆ log1mexp()

template<std::floating_point Real>

Real randlib::math::log1mexp ( Real x )

inline

Parameters

x

Returns: log(1 - exp(x))

◆ log1pmx()

template<std::floating_point Real>

Real randlib::math::log1pmx ( Real t )

inline

Parameters

t

Returns: log(1 + t) - t, accurate also for small |t|

◆ log2mexp()

template<std::floating_point Real>

Real randlib::math::log2mexp ( Real x )

inline

Parameters

x

Returns: log(2 - exp(x))

◆ logBesselI()

template<std::floating_point Real>

Real randlib::math::logBesselI	(	std::type_identity_t< Real >	nu,
		Real	x
	)

inline

Calculates logarithm of modified Bessel function of the 1st kind.

Parameters

nu
x

Returns: log(I_ν(x))

log(sinh(x)) or log(cosh(x))

I(ν, x) = I(-ν, x) - 2 / π sin(πν) K(-ν, x), and sin(πν) = -sin(π|ν|) for ν < 0

Large order: I_ν(x) = (x/2)^ν / Γ(ν + 1) Σ_j (x²/4)^j / (j! (ν + 1)_j) underflows, its logarithm does not

log(sinh(x)) or log(cosh(x))

I(ν, x) = I(-ν, x) - 2 / π sin(πν) K(-ν, x), and sin(πν) = -sin(π|ν|) for ν < 0

Large order: I_ν(x) = (x/2)^ν / Γ(ν + 1) Σ_j (x²/4)^j / (j! (ν + 1)_j) underflows, its logarithm does not

◆ logBesselK()

template<std::floating_point Real>

Real randlib::math::logBesselK	(	std::type_identity_t< Real >	nu,
		Real	x
	)

inline

Calculates logarithm of modified Bessel function of the 2nd kind.

Parameters

nu
x

Returns: log(K_ν(x))

K(-ν, x) = K(ν, x)

◆ logBeta() [1/2]

template<typename Left , typename Right >
requires (std::is_arithmetic_v<Left> && std::is_arithmetic_v<Right> && (std::integral<Left> || std::integral<Right> || !std::same_as<Left, Right>))

auto randlib::math::logBeta	(	Left	a,
		Right	b
	)		-> std::common_type_t<double, Left, Right>

inline

◆ logBeta() [2/2]

template<std::floating_point Real>

Real randlib::math::logBeta	(	Real	a,
		std::type_identity_t< Real >	b
	)

inline

Calculate logarithm of beta function.

Parameters

a	positive parameter
b	positive parameter

Returns: log(B(a, b)) = log(Γ(a)) + log(Γ(b)) - log(Γ(a + b))

Stirling corrections avoid both overflow of Γ and cancellation between large log Γ values

◆ logexpm1l()

template<std::floating_point Real>

Real randlib::math::logexpm1l ( Real x )

inline

Parameters

x

Returns: log(exp(x) - 1)

◆ logFactorialRatio()

template<std::floating_point Real = double>

Real randlib::math::logFactorialRatio	(	size_t	x,
		size_t	y
	)

inline

Parameters

x	non-negative integer number
y	non-negative integer number

Returns: log(x! / y!), through Stirling's series instead of a difference of large log-factorials

log(x!) = (x + 1/2) log(x + 1) - (x + 1) + log(2π) / 2 + stirlingError(x + 1)

◆ logPrefactor()

template<std::floating_point Real>

Real randlib::math::logPrefactor	(	Real	a,
		Real	x,
		Real	logX
	)

inline

Returns: log(x^a e^{-x} / Γ(a)); for a ≥ 10 it is formed as a·log1pmx(x/a - 1) + log(a / 2π) / 2 - (Stirling remainder of log Γ(a)), which keeps full relative accuracy when x is close to a large a

◆ lowerSeries()

template<std::floating_point Real>

Real randlib::math::lowerSeries	(	Real	a,
		Real	x
	)

inline

Returns: Σ_{n≥0} x^n / ((a+1)(a+2)...(a+n)), so that P(a, x) = x^a e^{-x} / Γ(a+1) times it

◆ lpgamma() [1/3]

template<std::floating_point Real>

Real randlib::math::lpgamma	(	std::type_identity_t< Real >	a,
		Real	x
	)

inline

Calculate logarithm of lower incomplete gamma function.

Parameters

a	non-negative parameter
x	non-negative parameter

Returns: log(P(a, x))

◆ lpgamma() [2/3]

template<std::floating_point Real>

Real randlib::math::lpgamma	(	std::type_identity_t< Real >	a,
		Real	x,
		std::type_identity_t< Real >	logA,
		std::type_identity_t< Real >	lgammaA
	)

inline

Calculate logarithm of lower incomplete gamma function, accelerated by using precalculated value of log(a) and log(Γ(a))

Parameters

a	non-negative parameter
x	non-negative parameter
logA	log(a)
lgammaA	log(Γ(a))

Returns: log(P(a, x))

◆ lpgamma() [3/3]

template<std::floating_point Real>

Real randlib::math::lpgamma	(	std::type_identity_t< Real >	a,
		Real	x,
		std::type_identity_t< Real >	logX
	)

inline

Calculate logarithm of lower incomplete gamma function, accelerated by using precalculated value of log(x)

Parameters

a	non-negative parameter
x	non-negative parameter
logX	log(x)

Returns: log(P(a, x))

◆ lqgamma() [1/3]

template<std::floating_point Real>

Real randlib::math::lqgamma	(	std::type_identity_t< Real >	a,
		Real	x
	)

inline

Calculate logarithm of upper incomplete gamma function.

Parameters

a	non-negative parameter
x	non-negative parameter

Returns: log(Q(a, x))

◆ lqgamma() [2/3]

template<std::floating_point Real>

Real randlib::math::lqgamma	(	std::type_identity_t< Real >	a,
		Real	x,
		std::type_identity_t< Real >	logA,
		std::type_identity_t< Real >	lgammaA
	)

inline

Calculate logarithm of upper incomplete gamma function, accelerated by using precalculated values of log(a) and log(Γ(a))

Parameters

a	non-negative parameter
x	non-negative parameter
logA	log(a)
lgammaA	log(Γ(a))

Returns: log(Q(a, x))

◆ lqgamma() [3/3]

template<std::floating_point Real>

Real randlib::math::lqgamma	(	std::type_identity_t< Real >	a,
		Real	x,
		std::type_identity_t< Real >	logX
	)

inline

Calculate logarithm of upper incomplete gamma function, accelerated by using precalculated value of log(x)

Parameters

a	non-negative parameter
x	non-negative parameter
logX	log(x)

Returns: log(Q(a, x))

◆ MarcumP() [1/2]

template<std::floating_point Real>

Real randlib::math::MarcumP	(	std::type_identity_t< Real >	mu,
		Real	x,
		std::type_identity_t< Real >	y
	)

inline

Parameters

mu
x
y

Returns: 1 - Marcum Q-function

◆ MarcumP() [2/2]

template<std::floating_point Real>

Real randlib::math::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
	)

inline

Parameters

mu
x
y
sqrtX	√x
sqrtY	√y
logX	log(x)
logY	log(y)

Returns: 1 - Marcum Q-function

Saddle-centred Poisson mixtures: the lower tail directly below the mean μ + x, the complement above it

◆ MarcumQ() [1/2]

template<std::floating_point Real>

Real randlib::math::MarcumQ	(	std::type_identity_t< Real >	mu,
		Real	x,
		std::type_identity_t< Real >	y
	)

inline

Parameters

mu
x
y

Returns: Marcum Q-function

◆ MarcumQ() [2/2]

template<std::floating_point Real>

Real randlib::math::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
	)

inline

Parameters

mu
x
y
sqrtX	√x
sqrtY	√y
logX	log(x)
logY	log(y)

Returns: Marcum Q-function

Below the mean μ + x the complement of the lower-tail mixture is accurate where Q is close to 1, and it agrees with MarcumP there

◆ nelderMead()

template<std::floating_point Real = double, class F >

MinimizeResult< Real > randlib::math::nelderMead	(	const F &	f,
		std::vector< Real >	start,
		Real	tolerance = `1e-10`,
		int	maxIterations = `20000`
	)

Downhill simplex method of Nelder and Mead with the dimension-dependent coefficients of Gao and Han (2012).

The initial simplex steps 5% along each coordinate (0.00025 for zero coordinates), as in MATLAB's fminsearch.

Parameters

f	function of a std::span<const Real>
start	starting point
tolerance	the iteration stops when the function values and the coordinates of the simplex agree to this relative tolerance
maxIterations	iteration limit

◆ numericalGradient()

template<std::floating_point Real = double, class F >

std::vector< Real > randlib::math::numericalGradient	(	const F &	f,
		const std::vector< Real > &	x
	)

Parameters

f	function of a std::span<const Real>
x	point

Returns: central-difference gradient with steps proportional to the cube root of the machine epsilon

◆ numericalHessian()

template<std::floating_point Real = double, class F >

std::vector< std::vector< Real > > randlib::math::numericalHessian	(	const F &	f,
		const std::vector< Real > &	x
	)

Parameters

f	function of a std::span<const Real>
x	point

Returns: symmetric central-difference Hessian with steps proportional to the fourth root of the machine epsilon

◆ parabolicMinimum()

template<typename RealType >

double randlib::math::parabolicMinimum	(	RealType	a,
		RealType	b,
		RealType	c,
		double	fa,
		double	fb,
		double	fc
	)

Parameters

a	point to the left of b
b	middle point
c	point to the right of b
fa	f(a)
fb	f(b)
fc	f(c)

Returns: minimum of interpolated parabola

◆ parametersVerification()

template<std::floating_point Real>

void randlib::math::parametersVerification	(	Real	a,
		Real	x
	)

inline

◆ pgamma() [1/3]

template<std::floating_point Real>

Real randlib::math::pgamma	(	std::type_identity_t< Real >	a,
		Real	x
	)

inline

Calculate lower regularized incomplete gamma function.

Parameters

a	non-negative parameter
x	non-negative parameter

Returns: P(a, x)

◆ pgamma() [2/3]

template<std::floating_point Real>

Real randlib::math::pgamma	(	std::type_identity_t< Real >	a,
		Real	x,
		std::type_identity_t< Real >	logA,
		std::type_identity_t< Real >	lgammaA
	)

inline

Calculate lower regularized incomplete gamma function, accelerated by using precalculated value of log(a) and log(Γ(a))

Parameters

a	non-negative parameter
x	non-negative parameter
logA	log(a)
lgammaA	log(Γ(a))

Returns: P(a, x)

◆ pgamma() [3/3]

template<std::floating_point Real>

Real randlib::math::pgamma	(	std::type_identity_t< Real >	a,
		Real	x,
		std::type_identity_t< Real >	logX
	)

inline

Calculate lower regularized incomplete gamma function, accelerated by using precalculated value of log(x)

Parameters

a	non-negative parameter
x	non-negative parameter
logX	log(x)

Returns: P(a, x)

◆ qgamma() [1/3]

template<std::floating_point Real>

Real randlib::math::qgamma	(	std::type_identity_t< Real >	a,
		Real	x
	)

inline

Calculate upper regularized incomplete gamma function.

Parameters

a	non-negative parameter
x	non-negative parameter

Returns: Q(a, x)

◆ qgamma() [2/3]

template<std::floating_point Real>

Real randlib::math::qgamma	(	std::type_identity_t< Real >	a,
		Real	x,
		std::type_identity_t< Real >	logA,
		std::type_identity_t< Real >	lgammaA
	)

inline

Calculate upper regularized incomplete gamma function, accelerated by using precalculated values of log(a) and log(Γ(a))

Parameters

a	non-negative parameter
x	non-negative parameter
logA	log(a)
lgammaA	log(Γ(a))

Returns: Q(a, x)

◆ qgamma() [3/3]

template<std::floating_point Real>

Real randlib::math::qgamma	(	std::type_identity_t< Real >	a,
		Real	x,
		std::type_identity_t< Real >	logX
	)

inline

Calculate upper regularized incomplete gamma function, accelerated by using precalculated value of log(x)

Parameters

a	non-negative parameter
x	non-negative parameter
logX	log(x)

Returns: Q(a, x)

◆ riemannZeta()

template<std::floating_point Real>

Real randlib::math::riemannZeta ( Real s )

inline

Parameters

s	argument, s > 1

Returns: Riemann zeta function ζ(s)

◆ sign()

template<std::floating_point Real>

int randlib::math::sign ( Real x )

inline

Parameters

x

Returns: sign of x

◆ softplus()

template<std::floating_point Real>

Real randlib::math::softplus ( Real x )

inline

Parameters

x

Returns: log(1 + exp(x))

◆ stirlingError()

template<std::floating_point Real = double>

Real randlib::math::stirlingError ( size_t n )

inline

Parameters

n	positive integer number

Returns: log(n!) - (n + 1/2) log(n) + n - log(2π)/2, the error of Stirling's formula

◆ toms748()

template<std::floating_point Real, class F >

std::pair< Real, Real > randlib::math::toms748	(	const F &	f,
		Real	a,
		Real	b,
		Real	fa,
		Real	fb,
		int	maxIterations = `200`
	)

Bracketing root finder of Alefeld, Potra and Shi (ACM TOMS algorithm 748): inverse cubic and Newton-quadratic interpolation steps, a double-length secant step and bisection whenever the bracket does not shrink by half; it converges like Brent's method at worst and faster on smooth functions.

Parameters

f	continuous function with f(a) and f(b) of opposite signs
a	lower end of the bracket
b	upper end of the bracket
fa	f(a)
fb	f(b)
maxIterations	maximal number of evaluations of f

Returns: the final bracket [a, b] around the root; a == b if f vanishes exactly

replaces [a, b] by the sub-bracket around c and keeps the dropped end in (d, fd)

double-length secant step from the end with the smaller residual

bisect when the bracket did not shrink by half

◆ trigamma()

template<std::floating_point Real>

Real randlib::math::trigamma ( Real x )

inline

Parameters

x	input parameter

Returns: trigamma(x) = digamma'(x)

Negative argument

set up tables

shift to minimum value, for which series expansion is applicable

choose the amount of terms in series

apply

Negative argument

set up tables

shift to minimum value, for which series expansion is applicable

choose the amount of terms in series

apply

◆ upperContinuedFraction()

template<std::floating_point Real>

Real randlib::math::upperContinuedFraction	(	Real	a,
		Real	x
	)

inline

Modified Lentz evaluation of the continued fraction for Q(a, x), x ≥ a + 1.

Returns: h such that Q(a, x) = x^a e^{-x} / Γ(a) · h

◆ W0Lambert() [1/2]

template<std::floating_point Real>

Real randlib::math::W0Lambert	(	Real	x,
		std::type_identity_t< Real >	epsilon = `Real{1e-11}`
	)

inline

Parameters

x
epsilon

Returns: W0 branch of Lambert W function

◆ W0Lambert() [2/2]

template<std::floating_point Real>

std::complex< Real > randlib::math::W0Lambert ( std::complex< Real > z )

inline

Parameters

z	complex argument off the branch cut (-∞, -1/e]

Returns: principal branch W0 of the Lambert W function: such w that w e^w = z

Winitzki's approximation, refined by Halley's iteration

◆ Wm1Lambert()

template<std::floating_point Real>

Real randlib::math::Wm1Lambert	(	Real	x,
		std::type_identity_t< Real >	epsilon = `Real{1e-11}`
	)

inline

Parameters

x
epsilon

Returns: W-1 branch of Lambert W function

◆ xexpxsqerfc()

template<std::floating_point Real>

Real randlib::math::xexpxsqerfc ( Real x )

inline

xexpxsqerfc

Parameters

x

Returns: x * exp(x^2) * erfc(x)

Variable Documentation

◆ FACTORIAL_TABLE

constexpr long double randlib::math::FACTORIAL_TABLE[]

constexpr

Initial value:

= {
    1.l,
    3628800.l,
    2432902008176640000.l,
    265252859812191058636308480000000.l,
    815915283247897734345611269596115894272000000000.l,
    30414093201713378043612608166064768844377641568960512000000000000.l,
    8320987112741390144276341183223364380754172606361245952449277696409600000000000000.l,
    11978571669969891796072783721689098736458938142546425857555362864628009582789845319680000000000000000.l,
    71569457046263802294811533723186532165584657342365752577109445058227039255480148842668944867280814080000000000000000000.l,
    1485715964481761497309522733620825737885569961284688766942216863704985393094065876545992131370884059645617234469978112000000000000000000000.l,
    93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000.l,
    15882455415227429404253703127090772871724410234473563207581748318444567162948183030959960131517678520479243672638179990208521148623422266876757623911219200000000000000000000000000.l,
    6689502913449127057588118054090372586752746333138029810295671352301633557244962989366874165271984981308157637893214090552534408589408121859898481114389650005964960521256960000000000000000000000000000.l,
    6466855489220473672507304395536485253155359447828049608975952322944781961185526165512707047229268452925683969240398027149120740074042105844737747799459310029635780991774612983803150965145600000000000000000000000000000000.l,
    13462012475717524605876073858941615558355851148193967190051391468057460367090535696797920946629681836680869097041958983702264048370902871114013579941370766400374327741701139895604871545254810788060989321379840000000000000000000000000000000000.l,
    57133839564458545904789328652610540031895535786011264182548375833179829124845398393126574488675311145377107878746854204162666250198684504466355949195922066574942592095735778929325357290444962472405416790722118445437122269675520000000000000000000000000000000000000.l,
    471472363599206132240694321176194377951192623045460204976904578317542573467421580346978030238114995699562728104819596262106947389303901748942909887857509625114880781313585012959529941660203611234871833992565791817698209861793313332044734813700096000000000000000000000000000000000000000.l,
    7257415615307998967396728211129263114716991681296451376543577798900561843401706157852350749242617459511490991237838520776666022565442753025328900773207510902400430280058295603966612599658257104398558294257568966313439612262571094946806711205568880457193340212661452800000000000000000000000000000000000000000.l,
}

◆ FACTORIAL_TABLESIZE

constexpr int randlib::math::FACTORIAL_TABLESIZE = 175

constexpr

◆ LOGFACTORIAL_TABLE

constexpr long double randlib::math::LOGFACTORIAL_TABLE[]

constexpr

◆ LOGFACTORIAL_TABLESIZE

constexpr size_t randlib::math::LOGFACTORIAL_TABLESIZE = 255

constexpr

Namespaces

Classes

Functions

Variables

Detailed Description

Function Documentation

◆ atan()

◆ bd0()

◆ besselI()

◆ besselK()

◆ bfgs()

◆ binom()

◆ createZiggurat()

◆ digamma()

◆ digammamLog()

◆ digammamLogForLargeX()

◆ doubleFactorial()

◆ erfcinv()

◆ erfcx()

◆ erfinv()

◆ factorial() [1/2]

◆ factorial() [2/2]

◆ factorialForSmallValue()

◆ findMin() [1/2]

◆ findMin() [2/2]

◆ findRoot()

◆ findRootNewtonFirstOrder()

◆ findRootNewtonFirstOrder2d()

◆ findRootNewtonSecondOrder()

◆ harmonicNumber()

◆ ibeta() [1/2]

◆ ibeta() [2/2]

◆ ibetacontinuedFraction1()

◆ ibetaContinuedFraction2()

◆ ibetaContinuedFractionLentz()

◆ ibetaLowerTail()

◆ ibetaPowerSeries1()

◆ ibetaPowerSeries2()

◆ incompleteGamma()

◆ integrate()

◆ integrateOscillatory() [1/2]

◆ integrateOscillatory() [2/2]

◆ integrateTanhSinh()

◆ ldfact()

◆ lfact()

◆ lgammaCorrection()

◆ log1mexp()

◆ log1pmx()

◆ log2mexp()

◆ logBesselI()

◆ logBesselK()

◆ logBeta() [1/2]

◆ logBeta() [2/2]

◆ logexpm1l()

◆ logFactorialRatio()

◆ logPrefactor()

◆ lowerSeries()

◆ lpgamma() [1/3]

◆ lpgamma() [2/3]

◆ lpgamma() [3/3]

◆ lqgamma() [1/3]

◆ lqgamma() [2/3]

◆ lqgamma() [3/3]

◆ MarcumP() [1/2]

◆ MarcumP() [2/2]

◆ MarcumQ() [1/2]

◆ MarcumQ() [2/2]

◆ nelderMead()

◆ numericalGradient()

◆ numericalHessian()

◆ parabolicMinimum()

◆ parametersVerification()

◆ pgamma() [1/3]

◆ pgamma() [2/3]

◆ pgamma() [3/3]

◆ qgamma() [1/3]

◆ qgamma() [2/3]

◆ qgamma() [3/3]

◆ riemannZeta()

◆ sign()