diff --git a/std/internal/math/gammafunction.d b/std/internal/math/gammafunction.d index 468f7eff8f5..d4e51f5653a 100644 --- a/std/internal/math/gammafunction.d +++ b/std/internal/math/gammafunction.d @@ -22,6 +22,9 @@ module std.internal.math.gammafunction; import std.internal.math.errorfunction; import std.math; import core.math : fabs, sin, sqrt; +import std.algorithm.iteration : fold; +import std.algorithm.searching : any; +import std.range : only; pure: nothrow: @@ -77,6 +80,34 @@ immutable real[8] logGammaDenominator = [ -0x1.00f95ced9e5f54eep+9L, 1.0L ]; + +/* Given a set of real values where at least one of them is NaN, it returns the + * NaN with the largest payload. This implements the NaN handling policy + * followed by D operators that accept multiple floating point arguments. + * + * It preserves the sign of the NaN. Also, when multiple NaNs have the largest + * payload, it returns the leftmost one. This means that + * largestNaNPayload(-NaN(1), NaN(1)) returns -NaN(1), and + * largestNaNPayload(NaN(1), -NaN(1)) returns NaN(1). + * + * When none of the provided values are NaN, the result is undefined. + */ +real largestNaNPayload(real first, real[] rest ...) +{ + return fold!((a, b) => cmp(abs(a), abs(b)) >= 0 ? a : b)(rest, first); +} + +@safe unittest +{ + assert(largestNaNPayload(NaN(0)) is NaN(0)); + assert(largestNaNPayload(NaN(1), NaN(0)) is NaN(1)); + assert(largestNaNPayload(NaN(2), NaN(3), NaN(1)) is NaN(3)); + assert(largestNaNPayload(-10.0L, -real.nan, 1.0L) is -real.nan); + assert(largestNaNPayload(-NaN(1), NaN(1)) is -NaN(1)); + assert(largestNaNPayload(NaN(1), -NaN(1)) is NaN(1)); +} + + /* * Helper function: Gamma function computed by Stirling's formula. * @@ -1012,14 +1043,7 @@ enum real BETA_BIGINV = 1.084202172485504434007e-19L; */ real betaIncomplete(real aa, real bb, real xx ) { - // If any parameters are NaN, return the NaN with the largest payload. - if (isNaN(aa) || isNaN(bb) || isNaN(xx)) - { - // With cmp, - // -NaN(larger) < -NaN(smaller) < -inf < inf < NaN(smaller) < NaN(larger). - const largerParam = cmp(abs(aa), abs(bb)) >= 0 ? aa : bb; - return cmp(abs(xx), abs(largerParam)) >= 0 ? xx : largerParam; - } + if (only(aa, bb, xx).any!isNaN) return largestNaNPayload(xx, aa, bb); // domain errors if (signbit(aa) == 1 || signbit(bb) == 1) return real.nan; @@ -1762,14 +1786,6 @@ real betaDistPowerSeries(real a, real b, real x ) /*************************************** * Incomplete gamma integral and its complement * - * These functions are defined by - * - * gammaIncomplete = ( $(INTEGRATE 0, x) $(POWER e, -t) $(POWER t, a-1) dt )/ $(GAMMA)(a) - * - * gammaIncompleteCompl(a,x) = 1 - gammaIncomplete(a,x) - * = ($(INTEGRATE x, ∞) $(POWER e, -t) $(POWER t, a-1) dt )/ $(GAMMA)(a) - * - * In this implementation both arguments must be positive. * The integral is evaluated by either a power series or * continued fraction expansion, depending on the relative * values of a and x. @@ -1777,22 +1793,26 @@ real betaDistPowerSeries(real a, real b, real x ) real gammaIncomplete(real a, real x ) in { - assert(x >= 0); - assert(a > 0); + if (!any!isNaN(only(a, x))) + { + assert(x >= 0); + assert(signbit(a) == 0); + } } do { - /* left tail of incomplete gamma function: - * - * inf. k - * a -x - x - * x e > ---------- - * - - - * k=0 | (a+k+1) - * - */ - if (x == 0) - return 0.0L; + // pass x first, so that if x and a are NaNs with the same payload but with + // opposite signs, return x. + if (any!isNaN(only(a, x))) return largestNaNPayload(x, a); + + // domain violation + if (signbit(a) == 1 || x < 0.0L) return real.nan; + + // edge cases + if (x == 0.0L) return 0.0L; + if (x is real.infinity) return 1.0L; + if (a is +0.0L) return 1.0L; + if (a is real.infinity) return 0.0L; if ( (x > 1.0L) && (x > a ) ) return 1.0L - gammaIncompleteCompl(a,x); @@ -1823,13 +1843,27 @@ do real gammaIncompleteCompl(real a, real x ) in { - assert(x >= 0); - assert(a > 0); + if (!any!isNaN(only(a, x))) + { + assert(x >= 0); + assert(signbit(a) == 0); + } } do { - if (x == 0) - return 1.0L; + // pass x first, so that if x and a are NaNs with the same payload but with + // opposite signs, return x. + if (any!isNaN(only(a, x))) return largestNaNPayload(x, a); + + // domain violation + if (signbit(a) == 1 || x < 0.0L) return real.nan; + + // edge cases + if (x == 0.0L) return 1.0L; + if (x is real.infinity) return 0.0L; + if (a is +0.0L) return 0.0L; + if (a is real.infinity) return 1.0L; + if ( (x < 1.0L) || (x < a) ) return 1.0L - gammaIncomplete(a,x); @@ -2039,6 +2073,18 @@ ihalve: @safe unittest { +assert(gammaIncomplete(NaN(4), -NaN(4)) is -NaN(4)); +assert(!isNaN(gammaIncomplete(+0.0L, 1.0L))); +assert(!isNaN(gammaIncomplete(1.0L, -0.0L))); +assert(gammaIncomplete(+0.0L, 0.0L) == 0.0L); +assert(gammaIncomplete(real.infinity, real.infinity) == 1.0L); + +assert(gammaIncompleteCompl(NaN(4), -NaN(4)) is -NaN(4)); +assert(!isNaN(gammaIncompleteCompl(+0.0L, 1.0L))); +assert(!isNaN(gammaIncompleteCompl(1.0L, -0.0L))); +assert(gammaIncompleteCompl(+0.0L, 0.0L) == 1.0L); +assert(gammaIncompleteCompl(real.infinity, real.infinity) == 0.0L); + //Values from Excel's GammaInv(1-p, x, 1) assert(fabs(gammaIncompleteComplInv(1, 0.5L) - 0.693147188044814L) < 0.00000005L); assert(fabs(gammaIncompleteComplInv(12, 0.99L) - 5.42818075054289L) < 0.00000005L); diff --git a/std/mathspecial.d b/std/mathspecial.d index 8889cbf6c67..aa538c393f7 100644 --- a/std/mathspecial.d +++ b/std/mathspecial.d @@ -29,6 +29,7 @@ * NAN = $(RED NAN) * SUP = $0 * GAMMA = Γ + * LGAMMA =γ * PSI = Ψ * THETA = θ * INTEGRAL = ∫ @@ -47,9 +48,11 @@ * LT = < * LE = ≤ * GT = > + * GE = ≥ * SQRT = √ * HALF = ½ * COMPLEX = ℂ + * NOBR = $1 * * Copyright: Based on the CEPHES math library, which is * Copyright (C) 1994 Stephen L. Moshier (moshier@world.std.com). @@ -396,43 +399,112 @@ real betaIncompleteInverse(real a, real b, real y ) return std.internal.math.gammafunction.betaIncompleteInv(a, b, y); } -/** Incomplete gamma integral and its complement +/** Regularized lower incomplete gamma function P(a,x) * - * These functions are defined by + * Mathematically, P(a,x) = $(LGAMMA)(a,x)/$(GAMMA)(a), where $(LGAMMA)(a,x) is the lower incomplete + * gamma function, $(LGAMMA)(a,x) = $(INTEGRATE 0, x)$(POWER t, a-1)$(POWER e, -t)dt, a $(GT) 0, and + * x $(GE) 0. * - * gammaIncomplete = ( $(INTEGRATE 0, x) $(POWER e, -t) $(POWER t, a-1) dt )/ $(GAMMA)(a) + * Params: + * a = the shape parameter, must be positive + * x = the fraction of integration completion from below, must be non-negative + * + * Returns: + * It returns P(a,x), an element of [0,1]. * - * gammaIncompleteCompl(a,x) = 1 - gammaIncomplete(a,x) - * = ($(INTEGRATE x, $(INFIN)) $(POWER e, -t) $(POWER t, a-1) dt )/ $(GAMMA)(a) + * $(TABLE_SV + * $(TR $(TH a) $(TH x) $(TH gammaIncomplete(a, x)) ) + * $(TR $(TD negative) $(TD) $(TD $(NAN)) ) + * $(TR $(TD) $(TD $(LT) 0) $(TD $(NAN)) ) + * $(TR $(TD positive) $(TD 0) $(TD 0) ) + * $(TR $(TD positive) $(TD $(INFIN)) $(TD 1) ) + * $(TR $(TD +0) $(TD $(GT) 0) $(TD 1) ) + * $(TR $(TD $(INFIN)) $(TD (0, $(INFIN))) $(TD 0) ) + * ) * - * In this implementation both arguments must be positive. - * The integral is evaluated by either a power series or - * continued fraction expansion, depending on the relative - * values of a and x. + * See_Also: $(LREF gamma) and $(LREF gammaIncompleteCompl) */ real gammaIncomplete(real a, real x ) in { - assert(x >= 0); - assert(a > 0); + // allow NaN input to pass through so that it can be addressed by the + // internal NaN payload propagation logic + if (!isNaN(a) && !isNaN(x)) + { + assert(signbit(a) == 0, "a must be positive"); + assert(x >= 0, "x must be non-negative"); + } } +out(p; isNaN(p) || (p >= 0 && p <= 1)) do { return std.internal.math.gammafunction.gammaIncomplete(a, x); } -/** ditto */ -real gammaIncompleteCompl(real a, real x ) +/// +@safe unittest +{ + assert(gammaIncomplete(1, 1) == 1 - 1/E); + assert(gammaIncomplete(1, 0) == 0); + assert(gammaIncomplete(1, real.infinity) == 1); + assert(gammaIncomplete(+0., 1) == 1); + assert(gammaIncomplete(real.infinity, 1) == 0); +} + +/** Regularized upper incomplete gamma function Q(a,x) + * + * Mathematically, Q(a,x) = $(GAMMA)(a,x)/$(GAMMA)(a), where $(GAMMA)(a,x) is the upper incomplete + * gamma function, $(GAMMA)(a,x) = $(INTEGRATE x, $(INFIN))$(POWER t, a-1)$(POWER e, -t)dt, + * $(NOBR a $(GT) 0), and x $(GE) 0. Note that P(a,x) + Q(a,x) = 1 or Q(a,x) = 1 - P(a,x), so Q is + * the complement of P. + * + * Params: + * a = the shape parameter, must be positive + * x = the fraction of integration completion from above, must be non-negative + * + * Returns: + * It returns Q(a,x), an element of [0,1]. + * + * $(TABLE_SV + * $(TR $(TH a) $(TH x) $(TH gammaIncompleteCompl(a, x)) ) + * $(TR $(TD negative) $(TD) $(TD $(NAN)) ) + * $(TR $(TD) $(TD $(LT) 0) $(TD $(NAN)) ) + * $(TR $(TD positive) $(TD 0) $(TD 1) ) + * $(TR $(TD positive) $(TD $(INFIN)) $(TD 0) ) + * $(TR $(TD +0) $(TD $(GT) 0) $(TD 0) ) + * $(TR $(TD $(INFIN)) $(TD (0, $(INFIN))) $(TD 1) ) + * ) + * + * See_Also: $(LREF gamma) and $(LREF gammaIncomplete) + */ + real gammaIncompleteCompl(real a, real x ) in { - assert(x >= 0); - assert(a > 0); + // allow NaN input to pass through so that it can be addressed by the + // internal NaN payload propagation logic + if (!isNaN(a) && !isNaN(x)) + { + assert(signbit(a) == 0, "a must be positive"); + assert(x >= 0, "x must be non-negative"); + } } +out(q; isNaN(q) || (q >= 0 && q <= 1)) do { return std.internal.math.gammafunction.gammaIncompleteCompl(a, x); } +/// +@safe unittest +{ + assert(isClose(gammaIncompleteCompl(2, 1), 2/E)); + assert(gammaIncompleteCompl(1, 0) == 1); + assert(gammaIncompleteCompl(1, real.infinity) == 0); + assert(gammaIncompleteCompl(+0., 1) == 0); + assert(gammaIncompleteCompl(real.infinity, 1) == 1); + assert(isClose(gammaIncompleteCompl(1, 2), 1-gammaIncomplete(1, 2))); +} + /** Inverse of complemented incomplete gamma integral * * Given a and p, the function finds x such that