Classes
struct	Segment

Functions
template<typename RealType >
bool	areClose (RealType a, RealType b, RealType eps=1e-6)

template<std::floating_point Real>
Real	erfinvChebyshevSeries (Real x, long double t, const long double *array, int size)

template<std::floating_point Real>
Real	erfinvAux1 (Real beta)

template<std::floating_point Real>
Real	erfinvAux2 (Real beta)

template<std::floating_point Real>
Real	erfinvAux3 (Real beta)

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

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

template<std::floating_point Real>
Real	WLambert (Real x, Real w0, Real epsilon)

template<std::floating_point Real>
Real	MarcumPSeries (Real mu, Real x, Real y, Real logX, Real logY)

template<std::floating_point Real>
Real	MarcumQSeries (Real mu, Real x, Real y, Real logX, Real logY)

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

template<std::floating_point Real>
Real	MarcumPAsymptoticForLargeXY (Real mu, Real x, Real y, Real sqrtX, Real sqrtY)

template<std::floating_point Real>
Real	MarcumPForMuSmallerThanOne (Real mu, Real x, Real y, Real logX, Real logY)

template<std::floating_point Real>
Real	MarcumQIntegrand (Real theta, Real xi, Real sqrt1pXiSq, Real mu, Real y)

template<std::floating_point Real>
Real	MarcumQIntergralRepresentation (Real mu, Real x, Real y, Real sqrtX, Real sqrtY)

template<std::floating_point Real>
std::pair< std::complex< Real >, Real >	accelerateOscillatorySums (const std::vector< std::complex< Real > > &partialSums)

template<std::floating_point Real, class F >
Segment< Real >	gaussKronrod21 (const F &f, Real a, Real b, Real &absEstimate)
	21-point Gauss-Kronrod rule with its embedded 10-point Gauss rule (QUADPACK qk21)

template<std::floating_point Real, class F >
Real	integrateFinite (const F &f, Real a, Real b, Real tolerance, int maxSubdivisions)
	Global adaptive scheme (QUADPACK QAG) on a finite interval.

template<std::floating_point Real, class F >
Real	evaluateTanhSinh (const F &f, Real x, Real endpointDistance)

template<std::floating_point Real, class F >
std::pair< Real, Real >	tanhSinhLevel (const F &f, Real a, Real b, Real step)

template<std::floating_point Real, class F , class TailMass >
std::complex< Real >	integrateHalfPeriods (const F &amplitude, Real frequency, Real lower, Real tolerance, int maxCycles, const TailMass &tailMass)

template<std::floating_point Real>
bool	bracketResolved (Real a, Real b)
	\|b - a\| no longer resolvable at the relative precision of Real

template<std::floating_point Real>
Real	safeDivide (Real numerator, Real denominator, Real fallback)

template<std::floating_point Real>
Real	secantStep (Real a, Real b, Real fa, Real fb)

template<std::floating_point Real>
Real	newtonQuadraticStep (Real a, Real b, Real d, Real fa, Real fb, Real fd, int steps)

template<std::floating_point Real>
Real	inverseCubicStep (Real a, Real b, Real d, Real e, Real fa, Real fb, Real fd, Real fe)

template<std::floating_point Real>
void	reciprocalGammaTerms (Real mu, Real &gam1, Real &gam2, Real &gampl, Real &gammi)
	Values needed by Temme's series for \|μ\| ≤ 1/2: gam1 = (1/Γ(1-μ) - 1/Γ(1+μ)) / (2μ), gam2 = (1/Γ(1-μ) + 1/Γ(1+μ)) / 2, gampl = 1/Γ(1+μ), gammi = 1/Γ(1-μ).

template<std::floating_point Real>
void	besselIK (Real nu, Real x, Real &ri, Real &rk)
	Temme's method with Steed's continued fraction (Numerical Recipes, "bessik")

template<std::floating_point Real>
void	checkBesselArguments (Real nu, Real x)

Variables
constexpr double	IBETA_LARGEST_FINITE_SHAPE = 4096
	Beyond this size of a shape the finite continued fractions of ibetaLowerTail accumulate rounding errors, and its power series converge at the rate x close to 1.

constexpr long double	minimumTolerance = 1e-21L

template<std::floating_point Real>
constexpr Real	machineEpsilon = std::numeric_limits<Real>::epsilon()

constexpr int	maxSeriesIterations = 100000

Detailed Description

Parameters

x
w0
epsilon

Returns

Function Documentation

◆ accelerateOscillatorySums()

template<std::floating_point Real>

std::pair< std::complex< Real >, Real > randlib::math::detail::accelerateOscillatorySums ( const std::vector< std::complex< Real > > & partialSums )

◆ areClose()

template<typename RealType >

bool randlib::math::detail::areClose	(	RealType	a,
		RealType	b,
		RealType	eps = `1e-6`
	)

Parameters

a
b
eps

Returns: |a - b| < eps * max(a, b)

◆ besselIK()

template<std::floating_point Real>

void randlib::math::detail::besselIK	(	Real	nu,
		Real	x,
		Real &	ri,
		Real &	rk
	)

inline

Temme's method with Steed's continued fraction (Numerical Recipes, "bessik")

Parameters

nu	order ν ≥ 0
x	argument x > 0
ri	I_ν(x)
rk	K_ν(x)

CF1 for I'_ν / I_ν by the modified Lentz method

Downward recurrence for I to order μ

Temme's series for K_μ and K_{μ+1}

Steed's CF2 for K_μ and K_{μ+1}

Wronskian gives I_μ, then scale the recurrence and recur K upwards

◆ bracketResolved()

template<std::floating_point Real>

bool randlib::math::detail::bracketResolved	(	Real	a,
		Real	b
	)

|b - a| no longer resolvable at the relative precision of Real

◆ checkBesselArguments()

template<std::floating_point Real>

void randlib::math::detail::checkBesselArguments	(	Real	nu,
		Real	x
	)

inline

◆ erfinvAux1()

template<std::floating_point Real>

Real randlib::math::detail::erfinvAux1 ( Real beta )

inline

|1 - p| < 5e-16

◆ erfinvAux2()

template<std::floating_point Real>

Real randlib::math::detail::erfinvAux2 ( Real beta )

inline

5e-16 < |1 - p| < 0.025

◆ erfinvAux3()

template<std::floating_point Real>

Real randlib::math::detail::erfinvAux3 ( Real beta )

inline

0.8 < p < 0.975

◆ erfinvAux4()

template<std::floating_point Real>

Real randlib::math::detail::erfinvAux4 ( Real p )

inline

0 < p < 0.8

◆ erfinvChebyshevSeries()

template<std::floating_point Real>

Real randlib::math::detail::erfinvChebyshevSeries	(	Real	x,
		long double	t,
		const long double *	array,
		int	size
	)

inline

We approximate inverse erf via Chebyshev polynomials

calculate next Chebyshev polynomial

◆ evaluateTanhSinh()

template<std::floating_point Real, class F >

Real randlib::math::detail::evaluateTanhSinh	(	const F &	f,
		Real	x,
		Real	endpointDistance
	)

◆ gaussKronrod21()

template<std::floating_point Real, class F >

Segment< Real > randlib::math::detail::gaussKronrod21	(	const F &	f,
		Real	a,
		Real	b,
		Real &	absEstimate
	)

21-point Gauss-Kronrod rule with its embedded 10-point Gauss rule (QUADPACK qk21)

Parameters

f	integrand
a	lower boundary
b	upper boundary
absEstimate	accumulates the integral of \|f\| over [a, b]

Returns: segment with the Kronrod estimate and |Kronrod - Gauss| as its error

◆ integrateFinite()

template<std::floating_point Real, class F >

Real randlib::math::detail::integrateFinite	(	const F &	f,
		Real	a,
		Real	b,
		Real	tolerance,
		int	maxSubdivisions
	)

Global adaptive scheme (QUADPACK QAG) on a finite interval.

Always bisect the segment with the largest error estimate until the total error is below the tolerance relative to the integral, or to a few ulps of the integral of |f| when the integrand cancels

QUADPACK's roundoff detection: bisections that stop reducing the error mean that rounding errors of the integrand dominate, and further subdivision cannot reach the tolerance

Re-sum to avoid drift from the running updates

◆ integrateHalfPeriods()

template<std::floating_point Real, class F , class TailMass >

std::complex< Real > randlib::math::detail::integrateHalfPeriods	(	const F &	amplitude,
		Real	frequency,
		Real	lower,
		Real	tolerance,
		int	maxCycles,
		const TailMass &	tailMass
	)

A known bound on the remaining ∫|amplitude| replaces the negligible-cycle heuristic, which mistakes empty half-periods below the bulk of the amplitude for a converged tail

◆ inverseCubicStep()

template<std::floating_point Real>

Real randlib::math::detail::inverseCubicStep	(	Real	a,
		Real	b,
		Real	d,
		Real	e,
		Real	fa,
		Real	fb,
		Real	fd,
		Real	fe
	)

◆ logHankelSeries()

template<std::floating_point Real>

Real randlib::math::detail::logHankelSeries	(	Real	nu,
		Real	x,
		int	sign
	)

inline

Parameters

nu	order
x	large argument
sign	-1 for I and +1 for K

Returns: log Σ_k sign^k a_k(ν) / x^k of the Hankel expansions of I and K, truncated at the smallest term

◆ MarcumPAsymptoticForLargeXY()

template<std::floating_point Real>

Real randlib::math::detail::MarcumPAsymptoticForLargeXY	(	Real	mu,
		Real	x,
		Real	y,
		Real	sqrtX,
		Real	sqrtY
	)

inline

Parameters

mu
x
y
sqrtX	√x
sqrtY	√y

Returns: asymptotic expansion for Marcum-P function for large x*y

sanity check

change φ

calculate A(μ-1) and A(μ)

compute term ψ and add it to the sum

◆ MarcumPForMuSmallerThanOne()

template<std::floating_point Real>

Real randlib::math::detail::MarcumPForMuSmallerThanOne	(	Real	mu,
		Real	x,
		Real	y,
		Real	logX,
		Real	logY
	)

inline

Parameters

mu
x
y
logX	precomputed log(x)
logY	precomputed log(y)

Returns

in this case we use numerical integration, however we have singularity point at 0, so we get rid of it by subtracting the function which has the same behaviour at this point

Calculate log of leveling factor

Calculate log(2f(t))

Return difference f(t) - factor

◆ MarcumPSaddleSeries()

template<std::floating_point Real>

Real randlib::math::detail::MarcumPSaddleSeries	(	Real	mu,
		Real	x,
		Real	y,
		Real	logX,
		Real	logY
	)

inline

Σ e^{-x} x^n / n! P(μ + n, y), summed outwards from the dominant Poisson term.

Below that term the lower incomplete gamma function grows by the stable recurrence P(a - 1, y) = P(a, y) + y^{a-1} e^{-y} / Γ(a)

◆ MarcumPSeries()

template<std::floating_point Real>

Real randlib::math::detail::MarcumPSeries	(	Real	mu,
		Real	x,
		Real	y,
		Real	logX,
		Real	logY
	)

inline

Parameters

mu
x
y
logX	log(x)
logY	log(y)

Returns: series expansion for Marcum-P function

~log(2πε) for ε = 1e-16

solving equation f(n) = 0 to find first negleted term

unexpected return

series expansion

sanity check

every 5 iterations we recalculate P and diffP in order to achieve enough accuracy

otherwise we use recurrent relations

add the last 0-term

◆ MarcumQIntegrand()

template<std::floating_point Real>

Real randlib::math::detail::MarcumQIntegrand	(	Real	theta,
		Real	xi,
		Real	sqrt1pXiSq,
		Real	mu,
		Real	y
	)

inline

◆ MarcumQIntergralRepresentation()

template<std::floating_point Real>

Real randlib::math::detail::MarcumQIntergralRepresentation	(	Real	mu,
		Real	x,
		Real	y,
		Real	sqrtX,
		Real	sqrtY
	)

inline

◆ MarcumQSeries()

template<std::floating_point Real>

Real randlib::math::detail::MarcumQSeries	(	Real	mu,
		Real	x,
		Real	y,
		Real	logX,
		Real	logY
	)

inline

◆ newtonQuadraticStep()

template<std::floating_point Real>

Real randlib::math::detail::newtonQuadraticStep	(	Real	a,
		Real	b,
		Real	d,
		Real	fa,
		Real	fb,
		Real	fd,
		int	steps
	)

◆ reciprocalGammaTerms()

template<std::floating_point Real>

void randlib::math::detail::reciprocalGammaTerms	(	Real	mu,
		Real &	gam1,
		Real &	gam2,
		Real &	gampl,
		Real &	gammi
	)

inline

Values needed by Temme's series for |μ| ≤ 1/2: gam1 = (1/Γ(1-μ) - 1/Γ(1+μ)) / (2μ), gam2 = (1/Γ(1-μ) + 1/Γ(1+μ)) / 2, gampl = 1/Γ(1+μ), gammi = 1/Γ(1-μ).

Odd part of the Taylor series 1/Γ(1+z) = Σ a_{k+1} z^k

◆ safeDivide()

template<std::floating_point Real>

Real randlib::math::detail::safeDivide	(	Real	numerator,
		Real	denominator,
		Real	fallback
	)

◆ secantStep()

template<std::floating_point Real>

Real randlib::math::detail::secantStep	(	Real	a,
		Real	b,
		Real	fa,
		Real	fb
	)

◆ tanhSinhLevel()

template<std::floating_point Real, class F >

std::pair< Real, Real > randlib::math::detail::tanhSinhLevel	(	const F &	f,
		Real	a,
		Real	b,
		Real	step
	)

◆ WLambert()

template<std::floating_point Real>

Real randlib::math::detail::WLambert	(	Real	x,
		Real	w0,
		Real	epsilon
	)

inline

Variable Documentation

◆ IBETA_LARGEST_FINITE_SHAPE

constexpr double randlib::math::detail::IBETA_LARGEST_FINITE_SHAPE = 4096

inlineconstexpr

Beyond this size of a shape the finite continued fractions of ibetaLowerTail accumulate rounding errors, and its power series converge at the rate x close to 1.

◆ machineEpsilon

template<std::floating_point Real>

constexpr Real randlib::math::detail::machineEpsilon = std::numeric_limits<Real>::epsilon()

inlineconstexpr

◆ maxSeriesIterations

constexpr int randlib::math::detail::maxSeriesIterations = 100000

constexpr

◆ minimumTolerance

constexpr long double randlib::math::detail::minimumTolerance = 1e-21L

inlineconstexpr

Classes

Functions

Variables

Detailed Description

Function Documentation

◆ accelerateOscillatorySums()

◆ areClose()

◆ besselIK()

◆ bracketResolved()

◆ checkBesselArguments()

◆ erfinvAux1()

◆ erfinvAux2()

◆ erfinvAux3()

◆ erfinvAux4()

◆ erfinvChebyshevSeries()

◆ evaluateTanhSinh()

◆ gaussKronrod21()

◆ integrateFinite()

◆ integrateHalfPeriods()

◆ inverseCubicStep()

◆ logHankelSeries()

◆ MarcumPAsymptoticForLargeXY()

◆ MarcumPForMuSmallerThanOne()

◆ MarcumPSaddleSeries()

◆ MarcumPSeries()

◆ MarcumQIntegrand()

◆ MarcumQIntergralRepresentation()

◆ MarcumQSeries()

◆ newtonQuadraticStep()

◆ reciprocalGammaTerms()

◆ safeDivide()

◆ secantStep()

◆ tanhSinhLevel()

◆ WLambert()

Variable Documentation

◆ IBETA_LARGEST_FINITE_SHAPE

◆ machineEpsilon

◆ maxSeriesIterations

◆ minimumTolerance