diff options
Diffstat (limited to 'libphobos/src/std/math/exponential.d')
| -rw-r--r-- | libphobos/src/std/math/exponential.d | 3439 |
1 files changed, 3439 insertions, 0 deletions
diff --git a/libphobos/src/std/math/exponential.d b/libphobos/src/std/math/exponential.d new file mode 100644 index 00000000000..a9ec930d90c --- /dev/null +++ b/libphobos/src/std/math/exponential.d @@ -0,0 +1,3439 @@ +// Written in the D programming language. + +/** +This is a submodule of $(MREF std, math). + +It contains several exponential and logarithm functions. + +Copyright: Copyright The D Language Foundation 2000 - 2011. + D implementations of exp, expm1, exp2, log, log10, log1p, and log2 + functions are based on the CEPHES math library, which is Copyright + (C) 2001 Stephen L. Moshier $(LT)steve@moshier.net$(GT) and are + incorporated herein by permission of the author. The author reserves + the right to distribute this material elsewhere under different + copying permissions. These modifications are distributed here under + the following terms: +License: $(HTTP www.boost.org/LICENSE_1_0.txt, Boost License 1.0). +Authors: $(HTTP digitalmars.com, Walter Bright), Don Clugston, + Conversion of CEPHES math library to D by Iain Buclaw and David Nadlinger +Source: $(PHOBOSSRC std/math/exponential.d) + +Macros: + TABLE_SV = <table border="1" cellpadding="4" cellspacing="0"> + <caption>Special Values</caption> + $0</table> + NAN = $(RED NAN) + PLUSMN = ± + INFIN = ∞ + PLUSMNINF = ±∞ + LT = < + GT = > + */ + +module std.math.exponential; + +import std.traits : isFloatingPoint, isIntegral, isSigned, isUnsigned, Largest, Unqual; + +static import core.math; +static import core.stdc.math; + +version (DigitalMars) +{ + version = INLINE_YL2X; // x87 has opcodes for these +} + +version (D_InlineAsm_X86) version = InlineAsm_X86_Any; +version (D_InlineAsm_X86_64) version = InlineAsm_X86_Any; + +version (InlineAsm_X86_Any) version = InlineAsm_X87; +version (InlineAsm_X87) +{ + static assert(real.mant_dig == 64); + version (CRuntime_Microsoft) version = InlineAsm_X87_MSVC; +} + +version (D_HardFloat) +{ + // FloatingPointControl.clearExceptions() depends on version IeeeFlagsSupport + version (IeeeFlagsSupport) version = FloatingPointControlSupport; +} + +/** + * Compute the value of x $(SUPERSCRIPT n), where n is an integer + */ +Unqual!F pow(F, G)(F x, G n) @nogc @trusted pure nothrow +if (isFloatingPoint!(F) && isIntegral!(G)) +{ + import std.traits : Unsigned; + + real p = 1.0, v = void; + Unsigned!(Unqual!G) m = n; + + if (n < 0) + { + if (n == -1) return 1 / x; + + m = cast(typeof(m))(0 - n); + v = p / x; + } + else + { + switch (n) + { + case 0: + return 1.0; + case 1: + return x; + case 2: + return x * x; + default: + } + + v = x; + } + + while (1) + { + if (m & 1) + p *= v; + m >>= 1; + if (!m) + break; + v *= v; + } + return p; +} + +/// +@safe pure nothrow @nogc unittest +{ + import std.math.operations : feqrel; + + assert(pow(2.0, 5) == 32.0); + assert(pow(1.5, 9).feqrel(38.4433) > 16); + assert(pow(real.nan, 2) is real.nan); + assert(pow(real.infinity, 2) == real.infinity); +} + +@safe pure nothrow @nogc unittest +{ + import std.math.operations : isClose, feqrel; + + // Make sure it instantiates and works properly on immutable values and + // with various integer and float types. + immutable real x = 46; + immutable float xf = x; + immutable double xd = x; + immutable uint one = 1; + immutable ushort two = 2; + immutable ubyte three = 3; + immutable ulong eight = 8; + + immutable int neg1 = -1; + immutable short neg2 = -2; + immutable byte neg3 = -3; + immutable long neg8 = -8; + + + assert(pow(x,0) == 1.0); + assert(pow(xd,one) == x); + assert(pow(xf,two) == x * x); + assert(pow(x,three) == x * x * x); + assert(pow(x,eight) == (x * x) * (x * x) * (x * x) * (x * x)); + + assert(pow(x, neg1) == 1 / x); + + assert(isClose(pow(xd, neg2), cast(double) (1 / (x * x)), 1e-15)); + assert(isClose(pow(xf, neg8), cast(float) (1 / ((x * x) * (x * x) * (x * x) * (x * x))), 1e-15)); + + assert(feqrel(pow(x, neg3), 1 / (x * x * x)) >= real.mant_dig - 1); +} + +@safe @nogc nothrow unittest +{ + import std.math.operations : isClose; + + assert(isClose(pow(2.0L, 10L), 1024, 1e-18)); +} + +// https://issues.dlang.org/show_bug.cgi?id=21601 +@safe @nogc nothrow pure unittest +{ + // When reals are large enough the results of pow(b, e) can be + // calculated correctly, if b is of type float or double and e is + // not too large. + static if (real.mant_dig >= 64) + { + // expected result: 3.790e-42 + assert(pow(-513645318757045764096.0f, -2) > 0.0); + + // expected result: 3.763915357831797e-309 + assert(pow(-1.6299717435255677e+154, -2) > 0.0); + } +} + +@safe @nogc nothrow unittest +{ + import std.math.operations : isClose; + import std.math.traits : isInfinity; + + static float f1 = 19100.0f; + static float f2 = 0.000012f; + + assert(isClose(pow(f1,9), 3.3829868e+38f)); + assert(isInfinity(pow(f1,10))); + assert(pow(f2,9) > 0.0f); + assert(isClose(pow(f2,10), 0.0f, 0.0, float.min_normal)); + + static double d1 = 21800.0; + static double d2 = 0.000012; + + assert(isClose(pow(d1,71), 1.0725339442974e+308)); + assert(isInfinity(pow(d1,72))); + assert(pow(d2,65) > 0.0f); + assert(isClose(pow(d2,66), 0.0, 0.0, double.min_normal)); + + static if (real.mant_dig == 64) // x87 + { + static real r1 = 21950.0L; + static real r2 = 0.000011L; + + assert(isClose(pow(r1,1136), 7.4066175654969242752260330529e+4931L)); + assert(isInfinity(pow(r1,1137))); + assert(pow(r2,998) > 0.0L); + assert(isClose(pow(r2,999), 0.0L, 0.0, real.min_normal)); + } +} + +@safe @nogc nothrow pure unittest +{ + import std.math.operations : isClose; + + enum f1 = 19100.0f; + enum f2 = 0.000012f; + + static assert(isClose(pow(f1,9), 3.3829868e+38f)); + static assert(pow(f1,10) > float.max); + static assert(pow(f2,9) > 0.0f); + static assert(isClose(pow(f2,10), 0.0f, 0.0, float.min_normal)); + + enum d1 = 21800.0; + enum d2 = 0.000012; + + static assert(isClose(pow(d1,71), 1.0725339442974e+308)); + static assert(pow(d1,72) > double.max); + static assert(pow(d2,65) > 0.0f); + static assert(isClose(pow(d2,66), 0.0, 0.0, double.min_normal)); + + static if (real.mant_dig == 64) // x87 + { + enum r1 = 21950.0L; + enum r2 = 0.000011L; + + static assert(isClose(pow(r1,1136), 7.4066175654969242752260330529e+4931L)); + static assert(pow(r1,1137) > real.max); + static assert(pow(r2,998) > 0.0L); + static assert(isClose(pow(r2,999), 0.0L, 0.0, real.min_normal)); + } +} + +/** + * Compute the power of two integral numbers. + * + * Params: + * x = base + * n = exponent + * + * Returns: + * x raised to the power of n. If n is negative the result is 1 / pow(x, -n), + * which is calculated as integer division with remainder. This may result in + * a division by zero error. + * + * If both x and n are 0, the result is 1. + * + * Throws: + * If x is 0 and n is negative, the result is the same as the result of a + * division by zero. + */ +typeof(Unqual!(F).init * Unqual!(G).init) pow(F, G)(F x, G n) @nogc @trusted pure nothrow +if (isIntegral!(F) && isIntegral!(G)) +{ + import std.traits : isSigned; + + typeof(return) p, v = void; + Unqual!G m = n; + + static if (isSigned!(F)) + { + if (x == -1) return cast(typeof(return)) (m & 1 ? -1 : 1); + } + static if (isSigned!(G)) + { + if (x == 0 && m <= -1) return x / 0; + } + if (x == 1) return 1; + static if (isSigned!(G)) + { + if (m < 0) return 0; + } + + switch (m) + { + case 0: + p = 1; + break; + + case 1: + p = x; + break; + + case 2: + p = x * x; + break; + + default: + v = x; + p = 1; + while (1) + { + if (m & 1) + p *= v; + m >>= 1; + if (!m) + break; + v *= v; + } + break; + } + return p; +} + +/// +@safe pure nothrow @nogc unittest +{ + assert(pow(2, 3) == 8); + assert(pow(3, 2) == 9); + + assert(pow(2, 10) == 1_024); + assert(pow(2, 20) == 1_048_576); + assert(pow(2, 30) == 1_073_741_824); + + assert(pow(0, 0) == 1); + + assert(pow(1, -5) == 1); + assert(pow(1, -6) == 1); + assert(pow(-1, -5) == -1); + assert(pow(-1, -6) == 1); + + assert(pow(-2, 5) == -32); + assert(pow(-2, -5) == 0); + assert(pow(cast(double) -2, -5) == -0.03125); +} + +@safe pure nothrow @nogc unittest +{ + immutable int one = 1; + immutable byte two = 2; + immutable ubyte three = 3; + immutable short four = 4; + immutable long ten = 10; + + assert(pow(two, three) == 8); + assert(pow(two, ten) == 1024); + assert(pow(one, ten) == 1); + assert(pow(ten, four) == 10_000); + assert(pow(four, 10) == 1_048_576); + assert(pow(three, four) == 81); +} + +// https://issues.dlang.org/show_bug.cgi?id=7006 +@safe pure nothrow @nogc unittest +{ + assert(pow(5, -1) == 0); + assert(pow(-5, -1) == 0); + assert(pow(5, -2) == 0); + assert(pow(-5, -2) == 0); + assert(pow(-1, int.min) == 1); + assert(pow(-2, int.min) == 0); + + assert(pow(4294967290UL,2) == 18446744022169944100UL); + assert(pow(0,uint.max) == 0); +} + +/**Computes integer to floating point powers.*/ +real pow(I, F)(I x, F y) @nogc @trusted pure nothrow +if (isIntegral!I && isFloatingPoint!F) +{ + return pow(cast(real) x, cast(Unqual!F) y); +} + +/// +@safe pure nothrow @nogc unittest +{ + assert(pow(2, 5.0) == 32.0); + assert(pow(7, 3.0) == 343.0); + assert(pow(2, real.nan) is real.nan); + assert(pow(2, real.infinity) == real.infinity); +} + +/** + * Calculates x$(SUPERSCRIPT y). + * + * $(TABLE_SV + * $(TR $(TH x) $(TH y) $(TH pow(x, y)) + * $(TH div 0) $(TH invalid?)) + * $(TR $(TD anything) $(TD $(PLUSMN)0.0) $(TD 1.0) + * $(TD no) $(TD no) ) + * $(TR $(TD |x| $(GT) 1) $(TD +$(INFIN)) $(TD +$(INFIN)) + * $(TD no) $(TD no) ) + * $(TR $(TD |x| $(LT) 1) $(TD +$(INFIN)) $(TD +0.0) + * $(TD no) $(TD no) ) + * $(TR $(TD |x| $(GT) 1) $(TD -$(INFIN)) $(TD +0.0) + * $(TD no) $(TD no) ) + * $(TR $(TD |x| $(LT) 1) $(TD -$(INFIN)) $(TD +$(INFIN)) + * $(TD no) $(TD no) ) + * $(TR $(TD +$(INFIN)) $(TD $(GT) 0.0) $(TD +$(INFIN)) + * $(TD no) $(TD no) ) + * $(TR $(TD +$(INFIN)) $(TD $(LT) 0.0) $(TD +0.0) + * $(TD no) $(TD no) ) + * $(TR $(TD -$(INFIN)) $(TD odd integer $(GT) 0.0) $(TD -$(INFIN)) + * $(TD no) $(TD no) ) + * $(TR $(TD -$(INFIN)) $(TD $(GT) 0.0, not odd integer) $(TD +$(INFIN)) + * $(TD no) $(TD no)) + * $(TR $(TD -$(INFIN)) $(TD odd integer $(LT) 0.0) $(TD -0.0) + * $(TD no) $(TD no) ) + * $(TR $(TD -$(INFIN)) $(TD $(LT) 0.0, not odd integer) $(TD +0.0) + * $(TD no) $(TD no) ) + * $(TR $(TD $(PLUSMN)1.0) $(TD $(PLUSMN)$(INFIN)) $(TD -$(NAN)) + * $(TD no) $(TD yes) ) + * $(TR $(TD $(LT) 0.0) $(TD finite, nonintegral) $(TD $(NAN)) + * $(TD no) $(TD yes)) + * $(TR $(TD $(PLUSMN)0.0) $(TD odd integer $(LT) 0.0) $(TD $(PLUSMNINF)) + * $(TD yes) $(TD no) ) + * $(TR $(TD $(PLUSMN)0.0) $(TD $(LT) 0.0, not odd integer) $(TD +$(INFIN)) + * $(TD yes) $(TD no)) + * $(TR $(TD $(PLUSMN)0.0) $(TD odd integer $(GT) 0.0) $(TD $(PLUSMN)0.0) + * $(TD no) $(TD no) ) + * $(TR $(TD $(PLUSMN)0.0) $(TD $(GT) 0.0, not odd integer) $(TD +0.0) + * $(TD no) $(TD no) ) + * ) + */ +Unqual!(Largest!(F, G)) pow(F, G)(F x, G y) @nogc @trusted pure nothrow +if (isFloatingPoint!(F) && isFloatingPoint!(G)) +{ + import core.math : fabs, sqrt; + import std.math.traits : isInfinity, isNaN, signbit; + + alias Float = typeof(return); + + static real impl(real x, real y) @nogc pure nothrow + { + // Special cases. + if (isNaN(y)) + return y; + if (isNaN(x) && y != 0.0) + return x; + + // Even if x is NaN. + if (y == 0.0) + return 1.0; + if (y == 1.0) + return x; + + if (isInfinity(y)) + { + if (isInfinity(x)) + { + if (!signbit(y) && !signbit(x)) + return F.infinity; + else + return F.nan; + } + else if (fabs(x) > 1) + { + if (signbit(y)) + return +0.0; + else + return F.infinity; + } + else if (fabs(x) == 1) + { + return F.nan; + } + else // < 1 + { + if (signbit(y)) + return F.infinity; + else + return +0.0; + } + } + if (isInfinity(x)) + { + if (signbit(x)) + { + long i = cast(long) y; + if (y > 0.0) + { + if (i == y && i & 1) + return -F.infinity; + else if (i == y) + return F.infinity; + else + return -F.nan; + } + else if (y < 0.0) + { + if (i == y && i & 1) + return -0.0; + else if (i == y) + return +0.0; + else + return F.nan; + } + } + else + { + if (y > 0.0) + return F.infinity; + else if (y < 0.0) + return +0.0; + } + } + + if (x == 0.0) + { + if (signbit(x)) + { + long i = cast(long) y; + if (y > 0.0) + { + if (i == y && i & 1) + return -0.0; + else + return +0.0; + } + else if (y < 0.0) + { + if (i == y && i & 1) + return -F.infinity; + else + return F.infinity; + } + } + else + { + if (y > 0.0) + return +0.0; + else if (y < 0.0) + return F.infinity; + } + } + if (x == 1.0) + return 1.0; + + if (y >= F.max) + { + if ((x > 0.0 && x < 1.0) || (x > -1.0 && x < 0.0)) + return 0.0; + if (x > 1.0 || x < -1.0) + return F.infinity; + } + if (y <= -F.max) + { + if ((x > 0.0 && x < 1.0) || (x > -1.0 && x < 0.0)) + return F.infinity; + if (x > 1.0 || x < -1.0) + return 0.0; + } + + if (x >= F.max) + { + if (y > 0.0) + return F.infinity; + else + return 0.0; + } + if (x <= -F.max) + { + long i = cast(long) y; + if (y > 0.0) + { + if (i == y && i & 1) + return -F.infinity; + else + return F.infinity; + } + else if (y < 0.0) + { + if (i == y && i & 1) + return -0.0; + else + return +0.0; + } + } + + // Integer power of x. + long iy = cast(long) y; + if (iy == y && fabs(y) < 32_768.0) + return pow(x, iy); + + real sign = 1.0; + if (x < 0) + { + // Result is real only if y is an integer + // Check for a non-zero fractional part + enum maxOdd = pow(2.0L, real.mant_dig) - 1.0L; + static if (maxOdd > ulong.max) + { + // Generic method, for any FP type + import std.math.rounding : floor; + if (floor(y) != y) + return sqrt(x); // Complex result -- create a NaN + + const hy = 0.5 * y; + if (floor(hy) != hy) + sign = -1.0; + } + else + { + // Much faster, if ulong has enough precision + const absY = fabs(y); + if (absY <= maxOdd) + { + const uy = cast(ulong) absY; + if (uy != absY) + return sqrt(x); // Complex result -- create a NaN + + if (uy & 1) + sign = -1.0; + } + } + x = -x; + } + version (INLINE_YL2X) + { + // If x > 0, x ^^ y == 2 ^^ ( y * log2(x) ) + // TODO: This is not accurate in practice. A fast and accurate + // (though complicated) method is described in: + // "An efficient rounding boundary test for pow(x, y) + // in double precision", C.Q. Lauter and V. Lefèvre, INRIA (2007). + return sign * exp2( core.math.yl2x(x, y) ); + } + else + { + // If x > 0, x ^^ y == 2 ^^ ( y * log2(x) ) + // TODO: This is not accurate in practice. A fast and accurate + // (though complicated) method is described in: + // "An efficient rounding boundary test for pow(x, y) + // in double precision", C.Q. Lauter and V. Lefèvre, INRIA (2007). + Float w = exp2(y * log2(x)); + return sign * w; + } + } + return impl(x, y); +} + +/// +@safe pure nothrow @nogc unittest +{ + import std.math.operations : isClose; + + assert(isClose(pow(2.0, 3.0), 8.0)); + assert(isClose(pow(1.5, 10.0), 57.6650390625)); + + // square root of 9 + assert(isClose(pow(9.0, 0.5), 3.0)); + // 10th root of 1024 + assert(isClose(pow(1024.0, 0.1), 2.0)); + + assert(isClose(pow(-4.0, 3.0), -64.0)); + + // reciprocal of 4 ^^ 2 + assert(isClose(pow(4.0, -2.0), 0.0625)); + // reciprocal of (-2) ^^ 3 + assert(isClose(pow(-2.0, -3.0), -0.125)); + + assert(isClose(pow(-2.5, 3.0), -15.625)); + // reciprocal of 2.5 ^^ 3 + assert(isClose(pow(2.5, -3.0), 0.064)); + // reciprocal of (-2.5) ^^ 3 + assert(isClose(pow(-2.5, -3.0), -0.064)); + + // reciprocal of square root of 4 + assert(isClose(pow(4.0, -0.5), 0.5)); + + // per definition + assert(isClose(pow(0.0, 0.0), 1.0)); +} + +/// +@safe pure nothrow @nogc unittest +{ + import std.math.operations : isClose; + + // the result is a complex number + // which cannot be represented as floating point number + import std.math.traits : isNaN; + assert(isNaN(pow(-2.5, -1.5))); + + // use the ^^-operator of std.complex instead + import std.complex : complex; + auto c1 = complex(-2.5, 0.0); + auto c2 = complex(-1.5, 0.0); + auto result = c1 ^^ c2; + // exact result apparently depends on `real` precision => increased tolerance + assert(isClose(result.re, -4.64705438e-17, 2e-4)); + assert(isClose(result.im, 2.52982e-1, 2e-4)); +} + +@safe pure nothrow @nogc unittest +{ + import std.math.traits : isNaN; + + assert(pow(1.5, real.infinity) == real.infinity); + assert(pow(0.5, real.infinity) == 0.0); + assert(pow(1.5, -real.infinity) == 0.0); + assert(pow(0.5, -real.infinity) == real.infinity); + assert(pow(real.infinity, 1.0) == real.infinity); + assert(pow(real.infinity, -1.0) == 0.0); + assert(pow(real.infinity, real.infinity) == real.infinity); + assert(pow(-real.infinity, 1.0) == -real.infinity); + assert(pow(-real.infinity, 2.0) == real.infinity); + assert(pow(-real.infinity, -1.0) == -0.0); + assert(pow(-real.infinity, -2.0) == 0.0); + assert(isNaN(pow(1.0, real.infinity))); + assert(pow(0.0, -1.0) == real.infinity); + assert(pow(real.nan, 0.0) == 1.0); + assert(isNaN(pow(real.nan, 3.0))); + assert(isNaN(pow(3.0, real.nan))); +} + +@safe @nogc nothrow unittest +{ + import std.math.operations : isClose; + + assert(isClose(pow(2.0L, 10.0L), 1024, 1e-18)); +} + +@safe pure nothrow @nogc unittest +{ + import std.math.operations : isClose; + import std.math.traits : isIdentical, isNaN; + import std.math.constants : PI; + + // Test all the special values. These unittests can be run on Windows + // by temporarily changing the version (linux) to version (all). + immutable float zero = 0; + immutable real one = 1; + immutable double two = 2; + immutable float three = 3; + immutable float fnan = float.nan; + immutable double dnan = double.nan; + immutable real rnan = real.nan; + immutable dinf = double.infinity; + immutable rninf = -real.infinity; + + assert(pow(fnan, zero) == 1); + assert(pow(dnan, zero) == 1); + assert(pow(rnan, zero) == 1); + + assert(pow(two, dinf) == double.infinity); + assert(isIdentical(pow(0.2f, dinf), +0.0)); + assert(pow(0.99999999L, rninf) == real.infinity); + assert(isIdentical(pow(1.000000001, rninf), +0.0)); + assert(pow(dinf, 0.001) == dinf); + assert(isIdentical(pow(dinf, -0.001), +0.0)); + assert(pow(rninf, 3.0L) == rninf); + assert(pow(rninf, 2.0L) == real.infinity); + assert(isIdentical(pow(rninf, -3.0), -0.0)); + assert(isIdentical(pow(rninf, -2.0), +0.0)); + + // @@@BUG@@@ somewhere + version (OSX) {} else assert(isNaN(pow(one, dinf))); + version (OSX) {} else assert(isNaN(pow(-one, dinf))); + assert(isNaN(pow(-0.2, PI))); + // boundary cases. Note that epsilon == 2^^-n for some n, + // so 1/epsilon == 2^^n is always even. + assert(pow(-1.0L, 1/real.epsilon - 1.0L) == -1.0L); + assert(pow(-1.0L, 1/real.epsilon) == 1.0L); + assert(isNaN(pow(-1.0L, 1/real.epsilon-0.5L))); + assert(isNaN(pow(-1.0L, -1/real.epsilon+0.5L))); + + assert(pow(0.0, -3.0) == double.infinity); + assert(pow(-0.0, -3.0) == -double.infinity); + assert(pow(0.0, -PI) == double.infinity); + assert(pow(-0.0, -PI) == double.infinity); + assert(isIdentical(pow(0.0, 5.0), 0.0)); + assert(isIdentical(pow(-0.0, 5.0), -0.0)); + assert(isIdentical(pow(0.0, 6.0), 0.0)); + assert(isIdentical(pow(-0.0, 6.0), 0.0)); + + // https://issues.dlang.org/show_bug.cgi?id=14786 fixed + immutable real maxOdd = pow(2.0L, real.mant_dig) - 1.0L; + assert(pow(-1.0L, maxOdd) == -1.0L); + assert(pow(-1.0L, -maxOdd) == -1.0L); + assert(pow(-1.0L, maxOdd + 1.0L) == 1.0L); + assert(pow(-1.0L, -maxOdd + 1.0L) == 1.0L); + assert(pow(-1.0L, maxOdd - 1.0L) == 1.0L); + assert(pow(-1.0L, -maxOdd - 1.0L) == 1.0L); + + // Now, actual numbers. + assert(isClose(pow(two, three), 8.0)); + assert(isClose(pow(two, -2.5), 0.1767766953)); + + // Test integer to float power. + immutable uint twoI = 2; + assert(isClose(pow(twoI, three), 8.0)); +} + +// https://issues.dlang.org/show_bug.cgi?id=20508 +@safe pure nothrow @nogc unittest +{ + import std.math.traits : isNaN; + + assert(isNaN(pow(-double.infinity, 0.5))); + + assert(isNaN(pow(-real.infinity, real.infinity))); + assert(isNaN(pow(-real.infinity, -real.infinity))); + assert(isNaN(pow(-real.infinity, 1.234))); + assert(isNaN(pow(-real.infinity, -0.751))); + assert(pow(-real.infinity, 0.0) == 1.0); +} + +/** Computes the value of a positive integer `x`, raised to the power `n`, modulo `m`. + * + * Params: + * x = base + * n = exponent + * m = modulus + * + * Returns: + * `x` to the power `n`, modulo `m`. + * The return type is the largest of `x`'s and `m`'s type. + * + * The function requires that all values have unsigned types. + */ +Unqual!(Largest!(F, H)) powmod(F, G, H)(F x, G n, H m) +if (isUnsigned!F && isUnsigned!G && isUnsigned!H) +{ + import std.meta : AliasSeq; + + alias T = Unqual!(Largest!(F, H)); + static if (T.sizeof <= 4) + { + alias DoubleT = AliasSeq!(void, ushort, uint, void, ulong)[T.sizeof]; + } + + static T mulmod(T a, T b, T c) + { + static if (T.sizeof == 8) + { + static T addmod(T a, T b, T c) + { + b = c - b; + if (a >= b) + return a - b; + else + return c - b + a; + } + + T result = 0, tmp; + + b %= c; + while (a > 0) + { + if (a & 1) + result = addmod(result, b, c); + + a >>= 1; + b = addmod(b, b, c); + } + + return result; + } + else + { + DoubleT result = cast(DoubleT) (cast(DoubleT) a * cast(DoubleT) b); + return result % c; + } + } + + T base = x, result = 1, modulus = m; + Unqual!G exponent = n; + + while (exponent > 0) + { + if (exponent & 1) + result = mulmod(result, base, modulus); + + base = mulmod(base, base, modulus); + exponent >>= 1; + } + + return result; +} + +/// +@safe pure nothrow @nogc unittest +{ + assert(powmod(1U, 10U, 3U) == 1); + assert(powmod(3U, 2U, 6U) == 3); + assert(powmod(5U, 5U, 15U) == 5); + assert(powmod(2U, 3U, 5U) == 3); + assert(powmod(2U, 4U, 5U) == 1); + assert(powmod(2U, 5U, 5U) == 2); +} + +@safe pure nothrow @nogc unittest +{ + ulong a = 18446744073709551615u, b = 20u, c = 18446744073709551610u; + assert(powmod(a, b, c) == 95367431640625u); + a = 100; b = 7919; c = 18446744073709551557u; + assert(powmod(a, b, c) == 18223853583554725198u); + a = 117; b = 7919; c = 18446744073709551557u; + assert(powmod(a, b, c) == 11493139548346411394u); + a = 134; b = 7919; c = 18446744073709551557u; + assert(powmod(a, b, c) == 10979163786734356774u); + a = 151; b = 7919; c = 18446744073709551557u; + assert(powmod(a, b, c) == 7023018419737782840u); + a = 168; b = 7919; c = 18446744073709551557u; + assert(powmod(a, b, c) == 58082701842386811u); + a = 185; b = 7919; c = 18446744073709551557u; + assert(powmod(a, b, c) == 17423478386299876798u); + a = 202; b = 7919; c = 18446744073709551557u; + assert(powmod(a, b, c) == 5522733478579799075u); + a = 219; b = 7919; c = 18446744073709551557u; + assert(powmod(a, b, c) == 15230218982491623487u); + a = 236; b = 7919; c = 18446744073709551557u; + assert(powmod(a, b, c) == 5198328724976436000u); + + a = 0; b = 7919; c = 18446744073709551557u; + assert(powmod(a, b, c) == 0); + a = 123; b = 0; c = 18446744073709551557u; + assert(powmod(a, b, c) == 1); + + immutable ulong a1 = 253, b1 = 7919, c1 = 18446744073709551557u; + assert(powmod(a1, b1, c1) == 3883707345459248860u); + + uint x = 100 ,y = 7919, z = 1844674407u; + assert(powmod(x, y, z) == 1613100340u); + x = 134; y = 7919; z = 1844674407u; + assert(powmod(x, y, z) == 734956622u); + x = 151; y = 7919; z = 1844674407u; + assert(powmod(x, y, z) == 1738696945u); + x = 168; y = 7919; z = 1844674407u; + assert(powmod(x, y, z) == 1247580927u); + x = 185; y = 7919; z = 1844674407u; + assert(powmod(x, y, z) == 1293855176u); + x = 202; y = 7919; z = 1844674407u; + assert(powmod(x, y, z) == 1566963682u); + x = 219; y = 7919; z = 1844674407u; + assert(powmod(x, y, z) == 181227807u); + x = 236; y = 7919; z = 1844674407u; + assert(powmod(x, y, z) == 217988321u); + x = 253; y = 7919; z = 1844674407u; + assert(powmod(x, y, z) == 1588843243u); + + x = 0; y = 7919; z = 184467u; + assert(powmod(x, y, z) == 0); + x = 123; y = 0; z = 1844674u; + assert(powmod(x, y, z) == 1); + + immutable ubyte x1 = 117; + immutable uint y1 = 7919; + immutable uint z1 = 1844674407u; + auto res = powmod(x1, y1, z1); + assert(is(typeof(res) == uint)); + assert(res == 9479781u); + + immutable ushort x2 = 123; + immutable uint y2 = 203; + immutable ubyte z2 = 113; + auto res2 = powmod(x2, y2, z2); + assert(is(typeof(res2) == ushort)); + assert(res2 == 42u); +} + +/** + * Calculates e$(SUPERSCRIPT x). + * + * $(TABLE_SV + * $(TR $(TH x) $(TH e$(SUPERSCRIPT x)) ) + * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) ) + * $(TR $(TD -$(INFIN)) $(TD +0.0) ) + * $(TR $(TD $(NAN)) $(TD $(NAN)) ) + * ) + */ +pragma(inline, true) +real exp(real x) @trusted pure nothrow @nogc // TODO: @safe +{ + import std.math.constants : LOG2E; + + version (InlineAsm_X87) + { + // e^^x = 2^^(LOG2E*x) + // (This is valid because the overflow & underflow limits for exp + // and exp2 are so similar). + if (!__ctfe) + return exp2Asm(LOG2E*x); + } + return expImpl(x); +} + +/// ditto +pragma(inline, true) +double exp(double x) @safe pure nothrow @nogc { return __ctfe ? cast(double) exp(cast(real) x) : expImpl(x); } + +/// ditto +pragma(inline, true) +float exp(float x) @safe pure nothrow @nogc { return __ctfe ? cast(float) exp(cast(real) x) : expImpl(x); } + +/// +@safe unittest +{ + import std.math.operations : feqrel; + import std.math.constants : E; + + assert(exp(0.0) == 1.0); + assert(exp(3.0).feqrel(E * E * E) > 16); +} + +private T expImpl(T)(T x) @safe pure nothrow @nogc +{ + import std.math : floatTraits, RealFormat; + import std.math.traits : isNaN; + import std.math.rounding : floor; + import std.math.algebraic : poly; + import std.math.constants : LOG2E; + + alias F = floatTraits!T; + static if (F.realFormat == RealFormat.ieeeSingle) + { + static immutable T[6] P = [ + 5.0000001201E-1, + 1.6666665459E-1, + 4.1665795894E-2, + 8.3334519073E-3, + 1.3981999507E-3, + 1.9875691500E-4, + ]; + + enum T C1 = 0.693359375; + enum T C2 = -2.12194440e-4; + + // Overflow and Underflow limits. + enum T OF = 88.72283905206835; + enum T UF = -103.278929903431851103; // ln(2^-149) + } + else static if (F.realFormat == RealFormat.ieeeDouble) + { + // Coefficients for exp(x) + static immutable T[3] P = [ + 9.99999999999999999910E-1L, + 3.02994407707441961300E-2L, + 1.26177193074810590878E-4L, + ]; + static immutable T[4] Q = [ + 2.00000000000000000009E0L, + 2.27265548208155028766E-1L, + 2.52448340349684104192E-3L, + 3.00198505138664455042E-6L, + ]; + + // C1 + C2 = LN2. + enum T C1 = 6.93145751953125E-1; + enum T C2 = 1.42860682030941723212E-6; + + // Overflow and Underflow limits. + enum T OF = 7.09782712893383996732E2; // ln((1-2^-53) * 2^1024) + enum T UF = -7.451332191019412076235E2; // ln(2^-1075) + } + else static if (F.realFormat == RealFormat.ieeeExtended || + F.realFormat == RealFormat.ieeeExtended53) + { + // Coefficients for exp(x) + static immutable T[3] P = [ + 9.9999999999999999991025E-1L, + 3.0299440770744196129956E-2L, + 1.2617719307481059087798E-4L, + ]; + static immutable T[4] Q = [ + 2.0000000000000000000897E0L, + 2.2726554820815502876593E-1L, + 2.5244834034968410419224E-3L, + 3.0019850513866445504159E-6L, + ]; + + // C1 + C2 = LN2. + enum T C1 = 6.9314575195312500000000E-1L; + enum T C2 = 1.4286068203094172321215E-6L; + + // Overflow and Underflow limits. + enum T OF = 1.1356523406294143949492E4L; // ln((1-2^-64) * 2^16384) + enum T UF = -1.13994985314888605586758E4L; // ln(2^-16446) + } + else static if (F.realFormat == RealFormat.ieeeQuadruple) + { + // Coefficients for exp(x) - 1 + static immutable T[5] P = [ + 9.999999999999999999999999999999999998502E-1L, + 3.508710990737834361215404761139478627390E-2L, + 2.708775201978218837374512615596512792224E-4L, + 6.141506007208645008909088812338454698548E-7L, + 3.279723985560247033712687707263393506266E-10L + ]; + static immutable T[6] Q = [ + 2.000000000000000000000000000000000000150E0, + 2.368408864814233538909747618894558968880E-1L, + 3.611828913847589925056132680618007270344E-3L, + 1.504792651814944826817779302637284053660E-5L, + 1.771372078166251484503904874657985291164E-8L, + 2.980756652081995192255342779918052538681E-12L + ]; + + // C1 + C2 = LN2. + enum T C1 = 6.93145751953125E-1L; + enum T C2 = 1.428606820309417232121458176568075500134E-6L; + + // Overflow and Underflow limits. + enum T OF = 1.135583025911358400418251384584930671458833e4L; + enum T UF = -1.143276959615573793352782661133116431383730e4L; + } + else + static assert(0, "Not implemented for this architecture"); + + // Special cases. + if (isNaN(x)) + return x; + if (x > OF) + return real.infinity; + if (x < UF) + return 0.0; + + // Express: e^^x = e^^g * 2^^n + // = e^^g * e^^(n * LOG2E) + // = e^^(g + n * LOG2E) + T xx = floor((cast(T) LOG2E) * x + cast(T) 0.5); + const int n = cast(int) xx; + x -= xx * C1; + x -= xx * C2; + + static if (F.realFormat == RealFormat.ieeeSingle) + { + xx = x * x; + x = poly(x, P) * xx + x + 1.0f; + } + else + { + // Rational approximation for exponential of the fractional part: + // e^^x = 1 + 2x P(x^^2) / (Q(x^^2) - P(x^^2)) + xx = x * x; + const T px = x * poly(xx, P); + x = px / (poly(xx, Q) - px); + x = (cast(T) 1.0) + (cast(T) 2.0) * x; + } + + // Scale by power of 2. + x = core.math.ldexp(x, n); + + return x; +} + +@safe @nogc nothrow unittest +{ + import std.math : floatTraits, RealFormat; + import std.math.operations : NaN, feqrel, isClose; + import std.math.constants : E; + import std.math.traits : isIdentical; + import std.math.algebraic : abs; + + version (IeeeFlagsSupport) import std.math.hardware : IeeeFlags, resetIeeeFlags, ieeeFlags; + version (FloatingPointControlSupport) + { + import std.math.hardware : FloatingPointControl; + + FloatingPointControl ctrl; + if (FloatingPointControl.hasExceptionTraps) + ctrl.disableExceptions(FloatingPointControl.allExceptions); + ctrl.rounding = FloatingPointControl.roundToNearest; + } + + static void testExp(T)() + { + enum realFormat = floatTraits!T.realFormat; + static if (realFormat == RealFormat.ieeeQuadruple) + { + static immutable T[2][] exptestpoints = + [ // x exp(x) + [ 1.0L, E ], + [ 0.5L, 0x1.a61298e1e069bc972dfefab6df34p+0L ], + [ 3.0L, E*E*E ], + [ 0x1.6p+13L, 0x1.6e509d45728655cdb4840542acb5p+16250L ], // near overflow + [ 0x1.7p+13L, T.infinity ], // close overflow + [ 0x1p+80L, T.infinity ], // far overflow + [ T.infinity, T.infinity ], + [-0x1.18p+13L, 0x1.5e4bf54b4807034ea97fef0059a6p-12927L ], // near underflow + [-0x1.625p+13L, 0x1.a6bd68a39d11fec3a250cd97f524p-16358L ], // ditto + [-0x1.62dafp+13L, 0x0.cb629e9813b80ed4d639e875be6cp-16382L ], // near underflow - subnormal + [-0x1.6549p+13L, 0x0.0000000000000000000000000001p-16382L ], // ditto + [-0x1.655p+13L, 0 ], // close underflow + [-0x1p+30L, 0 ], // far underflow + ]; + } + else static if (realFormat == RealFormat.ieeeExtended || + realFormat == RealFormat.ieeeExtended53) + { + static immutable T[2][] exptestpoints = + [ // x exp(x) + [ 1.0L, E ], + [ 0.5L, 0x1.a61298e1e069bc97p+0L ], + [ 3.0L, E*E*E ], + [ 0x1.1p+13L, 0x1.29aeffefc8ec645p+12557L ], // near overflow + [ 0x1.7p+13L, T.infinity ], // close overflow + [ 0x1p+80L, T.infinity ], // far overflow + [ T.infinity, T.infinity ], + [-0x1.18p+13L, 0x1.5e4bf54b4806db9p-12927L ], // near underflow + [-0x1.625p+13L, 0x1.a6bd68a39d11f35cp-16358L ], // ditto + [-0x1.62dafp+13L, 0x1.96c53d30277021dp-16383L ], // near underflow - subnormal + [-0x1.643p+13L, 0x1p-16444L ], // ditto + [-0x1.645p+13L, 0 ], // close underflow + [-0x1p+30L, 0 ], // far underflow + ]; + } + else static if (realFormat == RealFormat.ieeeDouble) + { + static immutable T[2][] exptestpoints = + [ // x, exp(x) + [ 1.0L, E ], + [ 0.5L, 0x1.a61298e1e069cp+0L ], + [ 3.0L, E*E*E ], + [ 0x1.6p+9L, 0x1.93bf4ec282efbp+1015L ], // near overflow + [ 0x1.7p+9L, T.infinity ], // close overflow + [ 0x1p+80L, T.infinity ], // far overflow + [ T.infinity, T.infinity ], + [-0x1.6p+9L, 0x1.44a3824e5285fp-1016L ], // near underflow + [-0x1.64p+9L, 0x0.06f84920bb2d4p-1022L ], // near underflow - subnormal + [-0x1.743p+9L, 0x0.0000000000001p-1022L ], // ditto + [-0x1.8p+9L, 0 ], // close underflow + [-0x1p+30L, 0 ], // far underflow + ]; + } + else static if (realFormat == RealFormat.ieeeSingle) + { + static immutable T[2][] exptestpoints = + [ // x, exp(x) + [ 1.0L, E ], + [ 0.5L, 0x1.a61299p+0L ], + [ 3.0L, E*E*E ], + [ 0x1.62p+6L, 0x1.99b988p+127L ], // near overflow + [ 0x1.7p+6L, T.infinity ], // close overflow + [ 0x1p+80L, T.infinity ], // far overflow + [ T.infinity, T.infinity ], + [-0x1.5cp+6L, 0x1.666d0ep-126L ], // near underflow + [-0x1.7p+6L, 0x0.026a42p-126L ], // near underflow - subnormal + [-0x1.9cp+6L, 0x0.000002p-126L ], // ditto + [-0x1.ap+6L, 0 ], // close underflow + [-0x1p+30L, 0 ], // far underflow + ]; + } + else + static assert(0, "No exp() tests for real type!"); + + const minEqualMantissaBits = T.mant_dig - 13; + T x; + version (IeeeFlagsSupport) IeeeFlags f; + foreach (ref pair; exptestpoints) + { + version (IeeeFlagsSupport) resetIeeeFlags(); + x = exp(pair[0]); + //printf("exp(%La) = %La, should be %La\n", cast(real) pair[0], cast(real) x, cast(real) pair[1]); + assert(feqrel(x, pair[1]) >= minEqualMantissaBits); + } + + // Ideally, exp(0) would not set the inexact flag. + // Unfortunately, fldl2e sets it! + // So it's not realistic to avoid setting it. + assert(exp(cast(T) 0.0) == 1.0); + + // NaN propagation. Doesn't set flags, bcos was already NaN. + version (IeeeFlagsSupport) + { + resetIeeeFlags(); + x = exp(T.nan); + f = ieeeFlags; + assert(isIdentical(abs(x), T.nan)); + assert(f.flags == 0); + + resetIeeeFlags(); + x = exp(-T.nan); + f = ieeeFlags; + assert(isIdentical(abs(x), T.nan)); + assert(f.flags == 0); + } + else + { + x = exp(T.nan); + assert(isIdentical(abs(x), T.nan)); + + x = exp(-T.nan); + assert(isIdentical(abs(x), T.nan)); + } + + x = exp(NaN(0x123)); + assert(isIdentical(x, NaN(0x123))); + } + + import std.meta : AliasSeq; + foreach (T; AliasSeq!(real, double, float)) + testExp!T(); + + // High resolution test (verified against GNU MPFR/Mathematica). + assert(exp(0.5L) == 0x1.A612_98E1_E069_BC97_2DFE_FAB6_DF34p+0L); + + assert(isClose(exp(3.0L), E * E * E, real.sizeof > double.sizeof ? 1e-15 : 1e-14)); +} + +/** + * Calculates the value of the natural logarithm base (e) + * raised to the power of x, minus 1. + * + * For very small x, expm1(x) is more accurate + * than exp(x)-1. + * + * $(TABLE_SV + * $(TR $(TH x) $(TH e$(SUPERSCRIPT x)-1) ) + * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) ) + * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) ) + * $(TR $(TD -$(INFIN)) $(TD -1.0) ) + * $(TR $(TD $(NAN)) $(TD $(NAN)) ) + * ) + */ +pragma(inline, true) +real expm1(real x) @trusted pure nothrow @nogc // TODO: @safe +{ + version (InlineAsm_X87) + { + if (!__ctfe) + return expm1Asm(x); + } + return expm1Impl(x); +} + +/// ditto +pragma(inline, true) +double expm1(double x) @safe pure nothrow @nogc +{ + return __ctfe ? cast(double) expm1(cast(real) x) : expm1Impl(x); +} + +/// ditto +pragma(inline, true) +float expm1(float x) @safe pure nothrow @nogc +{ + // no single-precision version in Cephes => use double precision + return __ctfe ? cast(float) expm1(cast(real) x) : cast(float) expm1Impl(cast(double) x); +} + +/// +@safe unittest +{ + import std.math.traits : isIdentical; + import std.math.operations : feqrel; + + assert(isIdentical(expm1(0.0), 0.0)); + assert(expm1(1.0).feqrel(1.71828) > 16); + assert(expm1(2.0).feqrel(6.3890) > 16); +} + +version (InlineAsm_X87) +private real expm1Asm(real x) @trusted pure nothrow @nogc +{ + version (X86) + { + enum PARAMSIZE = (real.sizeof+3)&(0xFFFF_FFFC); // always a multiple of 4 + asm pure nothrow @nogc + { + /* expm1() for x87 80-bit reals, IEEE754-2008 conformant. + * Author: Don Clugston. + * + * expm1(x) = 2^^(rndint(y))* 2^^(y-rndint(y)) - 1 where y = LN2*x. + * = 2rndy * 2ym1 + 2rndy - 1, where 2rndy = 2^^(rndint(y)) + * and 2ym1 = (2^^(y-rndint(y))-1). + * If 2rndy < 0.5*real.epsilon, result is -1. + * Implementation is otherwise the same as for exp2() + */ + naked; + fld real ptr [ESP+4] ; // x + mov AX, [ESP+4+8]; // AX = exponent and sign + sub ESP, 12+8; // Create scratch space on the stack + // [ESP,ESP+2] = scratchint + // [ESP+4..+6, +8..+10, +10] = scratchreal + // set scratchreal mantissa = 1.0 + mov dword ptr [ESP+8], 0; + mov dword ptr [ESP+8+4], 0x80000000; + and AX, 0x7FFF; // drop sign bit + cmp AX, 0x401D; // avoid InvalidException in fist + jae L_extreme; + fldl2e; + fmulp ST(1), ST; // y = x*log2(e) + fist dword ptr [ESP]; // scratchint = rndint(y) + fisub dword ptr [ESP]; // y - rndint(y) + // and now set scratchreal exponent + mov EAX, [ESP]; + add EAX, 0x3fff; + jle short L_largenegative; + cmp EAX,0x8000; + jge short L_largepositive; + mov [ESP+8+8],AX; + f2xm1; // 2ym1 = 2^^(y-rndint(y)) -1 + fld real ptr [ESP+8] ; // 2rndy = 2^^rndint(y) + fmul ST(1), ST; // ST=2rndy, ST(1)=2rndy*2ym1 + fld1; + fsubp ST(1), ST; // ST = 2rndy-1, ST(1) = 2rndy * 2ym1 - 1 + faddp ST(1), ST; // ST = 2rndy * 2ym1 + 2rndy - 1 + add ESP,12+8; + ret PARAMSIZE; + +L_extreme: // Extreme exponent. X is very large positive, very + // large negative, infinity, or NaN. + fxam; + fstsw AX; + test AX, 0x0400; // NaN_or_zero, but we already know x != 0 + jz L_was_nan; // if x is NaN, returns x + test AX, 0x0200; + jnz L_largenegative; +L_largepositive: + // Set scratchreal = real.max. + // squaring it will create infinity, and set overflow flag. + mov word ptr [ESP+8+8], 0x7FFE; + fstp ST(0); + fld real ptr [ESP+8]; // load scratchreal + fmul ST(0), ST; // square it, to create havoc! +L_was_nan: + add ESP,12+8; + ret PARAMSIZE; +L_largenegative: + fstp ST(0); + fld1; + fchs; // return -1. Underflow flag is not set. + add ESP,12+8; + ret PARAMSIZE; + } + } + else version (X86_64) + { + asm pure nothrow @nogc + { + naked; + } + version (Win64) + { + asm pure nothrow @nogc + { + fld real ptr [RCX]; // x + mov AX,[RCX+8]; // AX = exponent and sign + } + } + else + { + asm pure nothrow @nogc + { + fld real ptr [RSP+8]; // x + mov AX,[RSP+8+8]; // AX = exponent and sign + } + } + asm pure nothrow @nogc + { + /* expm1() for x87 80-bit reals, IEEE754-2008 conformant. + * Author: Don Clugston. + * + * expm1(x) = 2^(rndint(y))* 2^(y-rndint(y)) - 1 where y = LN2*x. + * = 2rndy * 2ym1 + 2rndy - 1, where 2rndy = 2^(rndint(y)) + * and 2ym1 = (2^(y-rndint(y))-1). + * If 2rndy < 0.5*real.epsilon, result is -1. + * Implementation is otherwise the same as for exp2() + */ + sub RSP, 24; // Create scratch space on the stack + // [RSP,RSP+2] = scratchint + // [RSP+4..+6, +8..+10, +10] = scratchreal + // set scratchreal mantissa = 1.0 + mov dword ptr [RSP+8], 0; + mov dword ptr [RSP+8+4], 0x80000000; + and AX, 0x7FFF; // drop sign bit + cmp AX, 0x401D; // avoid InvalidException in fist + jae L_extreme; + fldl2e; + fmul ; // y = x*log2(e) + fist dword ptr [RSP]; // scratchint = rndint(y) + fisub dword ptr [RSP]; // y - rndint(y) + // and now set scratchreal exponent + mov EAX, [RSP]; + add EAX, 0x3fff; + jle short L_largenegative; + cmp EAX,0x8000; + jge short L_largepositive; + mov [RSP+8+8],AX; + f2xm1; // 2^(y-rndint(y)) -1 + fld real ptr [RSP+8] ; // 2^rndint(y) + fmul ST(1), ST; + fld1; + fsubp ST(1), ST; + fadd; + add RSP,24; + ret; + +L_extreme: // Extreme exponent. X is very large positive, very + // large negative, infinity, or NaN. + fxam; + fstsw AX; + test AX, 0x0400; // NaN_or_zero, but we already know x != 0 + jz L_was_nan; // if x is NaN, returns x + test AX, 0x0200; + jnz L_largenegative; +L_largepositive: + // Set scratchreal = real.max. + // squaring it will create infinity, and set overflow flag. + mov word ptr [RSP+8+8], 0x7FFE; + fstp ST(0); + fld real ptr [RSP+8]; // load scratchreal + fmul ST(0), ST; // square it, to create havoc! +L_was_nan: + add RSP,24; + ret; + +L_largenegative: + fstp ST(0); + fld1; + fchs; // return -1. Underflow flag is not set. + add RSP,24; + ret; + } + } + else + static assert(0); +} + +private T expm1Impl(T)(T x) @safe pure nothrow @nogc +{ + import std.math : floatTraits, RealFormat; + import std.math.rounding : floor; + import std.math.algebraic : poly; + import std.math.constants : LN2; + + // Coefficients for exp(x) - 1 and overflow/underflow limits. + enum realFormat = floatTraits!T.realFormat; + static if (realFormat == RealFormat.ieeeQuadruple) + { + static immutable T[8] P = [ + 2.943520915569954073888921213330863757240E8L, + -5.722847283900608941516165725053359168840E7L, + 8.944630806357575461578107295909719817253E6L, + -7.212432713558031519943281748462837065308E5L, + 4.578962475841642634225390068461943438441E4L, + -1.716772506388927649032068540558788106762E3L, + 4.401308817383362136048032038528753151144E1L, + -4.888737542888633647784737721812546636240E-1L + ]; + + static immutable T[9] Q = [ + 1.766112549341972444333352727998584753865E9L, + -7.848989743695296475743081255027098295771E8L, + 1.615869009634292424463780387327037251069E8L, + -2.019684072836541751428967854947019415698E7L, + 1.682912729190313538934190635536631941751E6L, + -9.615511549171441430850103489315371768998E4L, + 3.697714952261803935521187272204485251835E3L, + -8.802340681794263968892934703309274564037E1L, + 1.0 + ]; + + enum T OF = 1.1356523406294143949491931077970764891253E4L; + enum T UF = -1.143276959615573793352782661133116431383730e4L; + } + else static if (realFormat == RealFormat.ieeeExtended) + { + static immutable T[5] P = [ + -1.586135578666346600772998894928250240826E4L, + 2.642771505685952966904660652518429479531E3L, + -3.423199068835684263987132888286791620673E2L, + 1.800826371455042224581246202420972737840E1L, + -5.238523121205561042771939008061958820811E-1L, + ]; + static immutable T[6] Q = [ + -9.516813471998079611319047060563358064497E4L, + 3.964866271411091674556850458227710004570E4L, + -7.207678383830091850230366618190187434796E3L, + 7.206038318724600171970199625081491823079E2L, + -4.002027679107076077238836622982900945173E1L, + 1.0 + ]; + + enum T OF = 1.1356523406294143949492E4L; + enum T UF = -4.5054566736396445112120088E1L; + } + else static if (realFormat == RealFormat.ieeeDouble) + { + static immutable T[3] P = [ + 9.9999999999999999991025E-1, + 3.0299440770744196129956E-2, + 1.2617719307481059087798E-4, + ]; + static immutable T[4] Q = [ + 2.0000000000000000000897E0, + 2.2726554820815502876593E-1, + 2.5244834034968410419224E-3, + 3.0019850513866445504159E-6, + ]; + } + else + static assert(0, "no coefficients for expm1()"); + + static if (realFormat == RealFormat.ieeeDouble) // special case for double precision + { + if (x < -0.5 || x > 0.5) + return exp(x) - 1.0; + if (x == 0.0) + return x; + + const T xx = x * x; + x = x * poly(xx, P); + x = x / (poly(xx, Q) - x); + return x + x; + } + else + { + // C1 + C2 = LN2. + enum T C1 = 6.9314575195312500000000E-1L; + enum T C2 = 1.428606820309417232121458176568075500134E-6L; + + // Special cases. + if (x > OF) + return real.infinity; + if (x == cast(T) 0.0) + return x; + if (x < UF) + return -1.0; + + // Express x = LN2 (n + remainder), remainder not exceeding 1/2. + int n = cast(int) floor((cast(T) 0.5) + x / cast(T) LN2); + x -= n * C1; + x -= n * C2; + + // Rational approximation: + // exp(x) - 1 = x + 0.5 x^^2 + x^^3 P(x) / Q(x) + T px = x * poly(x, P); + T qx = poly(x, Q); + const T xx = x * x; + qx = x + ((cast(T) 0.5) * xx + xx * px / qx); + + // We have qx = exp(remainder LN2) - 1, so: + // exp(x) - 1 = 2^^n (qx + 1) - 1 = 2^^n qx + 2^^n - 1. + px = core.math.ldexp(cast(T) 1.0, n); + x = px * qx + (px - cast(T) 1.0); + + return x; + } +} + +@safe @nogc nothrow unittest +{ + import std.math.traits : isNaN; + import std.math.operations : isClose, CommonDefaultFor; + + static void testExpm1(T)() + { + // NaN + assert(isNaN(expm1(cast(T) T.nan))); + + static immutable T[] xs = [ -2, -0.75, -0.3, 0.0, 0.1, 0.2, 0.5, 1.0 ]; + foreach (x; xs) + { + const T e = expm1(x); + const T r = exp(x) - 1; + + //printf("expm1(%Lg) = %Lg, should approximately be %Lg\n", cast(real) x, cast(real) e, cast(real) r); + assert(isClose(r, e, CommonDefaultFor!(T,T), CommonDefaultFor!(T,T))); + } + } + + import std.meta : AliasSeq; + foreach (T; AliasSeq!(real, double)) + testExpm1!T(); +} + +/** + * Calculates 2$(SUPERSCRIPT x). + * + * $(TABLE_SV + * $(TR $(TH x) $(TH exp2(x)) ) + * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) ) + * $(TR $(TD -$(INFIN)) $(TD +0.0) ) + * $(TR $(TD $(NAN)) $(TD $(NAN)) ) + * ) + */ +pragma(inline, true) +real exp2(real x) @nogc @trusted pure nothrow // TODO: @safe +{ + version (InlineAsm_X87) + { + if (!__ctfe) + return exp2Asm(x); + } + return exp2Impl(x); +} + +/// ditto +pragma(inline, true) +double exp2(double x) @nogc @safe pure nothrow { return __ctfe ? cast(double) exp2(cast(real) x) : exp2Impl(x); } + +/// ditto +pragma(inline, true) +float exp2(float x) @nogc @safe pure nothrow { return __ctfe ? cast(float) exp2(cast(real) x) : exp2Impl(x); } + +/// +@safe unittest +{ + import std.math.traits : isIdentical; + import std.math.operations : feqrel; + + assert(isIdentical(exp2(0.0), 1.0)); + assert(exp2(2.0).feqrel(4.0) > 16); + assert(exp2(8.0).feqrel(256.0) > 16); +} + +@safe unittest +{ + version (CRuntime_Microsoft) {} else // aexp2/exp2f/exp2l not implemented + { + assert( core.stdc.math.exp2f(0.0f) == 1 ); + assert( core.stdc.math.exp2 (0.0) == 1 ); + assert( core.stdc.math.exp2l(0.0L) == 1 ); + } +} + +version (InlineAsm_X87) +private real exp2Asm(real x) @nogc @trusted pure nothrow +{ + version (X86) + { + enum PARAMSIZE = (real.sizeof+3)&(0xFFFF_FFFC); // always a multiple of 4 + + asm pure nothrow @nogc + { + /* exp2() for x87 80-bit reals, IEEE754-2008 conformant. + * Author: Don Clugston. + * + * exp2(x) = 2^^(rndint(x))* 2^^(y-rndint(x)) + * The trick for high performance is to avoid the fscale(28cycles on core2), + * frndint(19 cycles), leaving f2xm1(19 cycles) as the only slow instruction. + * + * We can do frndint by using fist. BUT we can't use it for huge numbers, + * because it will set the Invalid Operation flag if overflow or NaN occurs. + * Fortunately, whenever this happens the result would be zero or infinity. + * + * We can perform fscale by directly poking into the exponent. BUT this doesn't + * work for the (very rare) cases where the result is subnormal. So we fall back + * to the slow method in that case. + */ + naked; + fld real ptr [ESP+4] ; // x + mov AX, [ESP+4+8]; // AX = exponent and sign + sub ESP, 12+8; // Create scratch space on the stack + // [ESP,ESP+2] = scratchint + // [ESP+4..+6, +8..+10, +10] = scratchreal + // set scratchreal mantissa = 1.0 + mov dword ptr [ESP+8], 0; + mov dword ptr [ESP+8+4], 0x80000000; + and AX, 0x7FFF; // drop sign bit + cmp AX, 0x401D; // avoid InvalidException in fist + jae L_extreme; + fist dword ptr [ESP]; // scratchint = rndint(x) + fisub dword ptr [ESP]; // x - rndint(x) + // and now set scratchreal exponent + mov EAX, [ESP]; + add EAX, 0x3fff; + jle short L_subnormal; + cmp EAX,0x8000; + jge short L_overflow; + mov [ESP+8+8],AX; +L_normal: + f2xm1; + fld1; + faddp ST(1), ST; // 2^^(x-rndint(x)) + fld real ptr [ESP+8] ; // 2^^rndint(x) + add ESP,12+8; + fmulp ST(1), ST; + ret PARAMSIZE; + +L_subnormal: + // Result will be subnormal. + // In this rare case, the simple poking method doesn't work. + // The speed doesn't matter, so use the slow fscale method. + fild dword ptr [ESP]; // scratchint + fld1; + fscale; + fstp real ptr [ESP+8]; // scratchreal = 2^^scratchint + fstp ST(0); // drop scratchint + jmp L_normal; + +L_extreme: // Extreme exponent. X is very large positive, very + // large negative, infinity, or NaN. + fxam; + fstsw AX; + test AX, 0x0400; // NaN_or_zero, but we already know x != 0 + jz L_was_nan; // if x is NaN, returns x + // set scratchreal = real.min_normal + // squaring it will return 0, setting underflow flag + mov word ptr [ESP+8+8], 1; + test AX, 0x0200; + jnz L_waslargenegative; +L_overflow: + // Set scratchreal = real.max. + // squaring it will create infinity, and set overflow flag. + mov word ptr [ESP+8+8], 0x7FFE; +L_waslargenegative: + fstp ST(0); + fld real ptr [ESP+8]; // load scratchreal + fmul ST(0), ST; // square it, to create havoc! +L_was_nan: + add ESP,12+8; + ret PARAMSIZE; + } + } + else version (X86_64) + { + asm pure nothrow @nogc + { + naked; + } + version (Win64) + { + asm pure nothrow @nogc + { + fld real ptr [RCX]; // x + mov AX,[RCX+8]; // AX = exponent and sign + } + } + else + { + asm pure nothrow @nogc + { + fld real ptr [RSP+8]; // x + mov AX,[RSP+8+8]; // AX = exponent and sign + } + } + asm pure nothrow @nogc + { + /* exp2() for x87 80-bit reals, IEEE754-2008 conformant. + * Author: Don Clugston. + * + * exp2(x) = 2^(rndint(x))* 2^(y-rndint(x)) + * The trick for high performance is to avoid the fscale(28cycles on core2), + * frndint(19 cycles), leaving f2xm1(19 cycles) as the only slow instruction. + * + * We can do frndint by using fist. BUT we can't use it for huge numbers, + * because it will set the Invalid Operation flag is overflow or NaN occurs. + * Fortunately, whenever this happens the result would be zero or infinity. + * + * We can perform fscale by directly poking into the exponent. BUT this doesn't + * work for the (very rare) cases where the result is subnormal. So we fall back + * to the slow method in that case. + */ + sub RSP, 24; // Create scratch space on the stack + // [RSP,RSP+2] = scratchint + // [RSP+4..+6, +8..+10, +10] = scratchreal + // set scratchreal mantissa = 1.0 + mov dword ptr [RSP+8], 0; + mov dword ptr [RSP+8+4], 0x80000000; + and AX, 0x7FFF; // drop sign bit + cmp AX, 0x401D; // avoid InvalidException in fist + jae L_extreme; + fist dword ptr [RSP]; // scratchint = rndint(x) + fisub dword ptr [RSP]; // x - rndint(x) + // and now set scratchreal exponent + mov EAX, [RSP]; + add EAX, 0x3fff; + jle short L_subnormal; + cmp EAX,0x8000; + jge short L_overflow; + mov [RSP+8+8],AX; +L_normal: + f2xm1; + fld1; + fadd; // 2^(x-rndint(x)) + fld real ptr [RSP+8] ; // 2^rndint(x) + add RSP,24; + fmulp ST(1), ST; + ret; + +L_subnormal: + // Result will be subnormal. + // In this rare case, the simple poking method doesn't work. + // The speed doesn't matter, so use the slow fscale method. + fild dword ptr [RSP]; // scratchint + fld1; + fscale; + fstp real ptr [RSP+8]; // scratchreal = 2^scratchint + fstp ST(0); // drop scratchint + jmp L_normal; + +L_extreme: // Extreme exponent. X is very large positive, very + // large negative, infinity, or NaN. + fxam; + fstsw AX; + test AX, 0x0400; // NaN_or_zero, but we already know x != 0 + jz L_was_nan; // if x is NaN, returns x + // set scratchreal = real.min + // squaring it will return 0, setting underflow flag + mov word ptr [RSP+8+8], 1; + test AX, 0x0200; + jnz L_waslargenegative; +L_overflow: + // Set scratchreal = real.max. + // squaring it will create infinity, and set overflow flag. + mov word ptr [RSP+8+8], 0x7FFE; +L_waslargenegative: + fstp ST(0); + fld real ptr [RSP+8]; // load scratchreal + fmul ST(0), ST; // square it, to create havoc! +L_was_nan: + add RSP,24; + ret; + } + } + else + static assert(0); +} + +private T exp2Impl(T)(T x) @nogc @safe pure nothrow +{ + import std.math : floatTraits, RealFormat; + import std.math.traits : isNaN; + import std.math.rounding : floor; + import std.math.algebraic : poly; + + // Coefficients for exp2(x) + enum realFormat = floatTraits!T.realFormat; + static if (realFormat == RealFormat.ieeeQuadruple) + { + static immutable T[5] P = [ + 9.079594442980146270952372234833529694788E12L, + 1.530625323728429161131811299626419117557E11L, + 5.677513871931844661829755443994214173883E8L, + 6.185032670011643762127954396427045467506E5L, + 1.587171580015525194694938306936721666031E2L + ]; + + static immutable T[6] Q = [ + 2.619817175234089411411070339065679229869E13L, + 1.490560994263653042761789432690793026977E12L, + 1.092141473886177435056423606755843616331E10L, + 2.186249607051644894762167991800811827835E7L, + 1.236602014442099053716561665053645270207E4L, + 1.0 + ]; + } + else static if (realFormat == RealFormat.ieeeExtended) + { + static immutable T[3] P = [ + 2.0803843631901852422887E6L, + 3.0286971917562792508623E4L, + 6.0614853552242266094567E1L, + ]; + static immutable T[4] Q = [ + 6.0027204078348487957118E6L, + 3.2772515434906797273099E5L, + 1.7492876999891839021063E3L, + 1.0000000000000000000000E0L, + ]; + } + else static if (realFormat == RealFormat.ieeeDouble) + { + static immutable T[3] P = [ + 1.51390680115615096133E3L, + 2.02020656693165307700E1L, + 2.30933477057345225087E-2L, + ]; + static immutable T[3] Q = [ + 4.36821166879210612817E3L, + 2.33184211722314911771E2L, + 1.00000000000000000000E0L, + ]; + } + else static if (realFormat == RealFormat.ieeeSingle) + { + static immutable T[6] P = [ + 6.931472028550421E-001L, + 2.402264791363012E-001L, + 5.550332471162809E-002L, + 9.618437357674640E-003L, + 1.339887440266574E-003L, + 1.535336188319500E-004L, + ]; + } + else + static assert(0, "no coefficients for exp2()"); + + // Overflow and Underflow limits. + enum T OF = T.max_exp; + enum T UF = T.min_exp - 1; + + // Special cases. + if (isNaN(x)) + return x; + if (x > OF) + return real.infinity; + if (x < UF) + return 0.0; + + static if (realFormat == RealFormat.ieeeSingle) // special case for single precision + { + // The following is necessary because range reduction blows up. + if (x == 0.0f) + return 1.0f; + + // Separate into integer and fractional parts. + const T i = floor(x); + int n = cast(int) i; + x -= i; + if (x > 0.5f) + { + n += 1; + x -= 1.0f; + } + + // Rational approximation: + // exp2(x) = 1.0 + x P(x) + x = 1.0f + x * poly(x, P); + } + else + { + // Separate into integer and fractional parts. + const T i = floor(x + cast(T) 0.5); + int n = cast(int) i; + x -= i; + + // Rational approximation: + // exp2(x) = 1.0 + 2x P(x^^2) / (Q(x^^2) - P(x^^2)) + const T xx = x * x; + const T px = x * poly(xx, P); + x = px / (poly(xx, Q) - px); + x = (cast(T) 1.0) + (cast(T) 2.0) * x; + } + + // Scale by power of 2. + x = core.math.ldexp(x, n); + + return x; +} + +@safe @nogc nothrow unittest +{ + import std.math.operations : feqrel, NaN, isClose; + import std.math.traits : isIdentical; + import std.math.constants : SQRT2; + + assert(feqrel(exp2(0.5L), SQRT2) >= real.mant_dig -1); + assert(exp2(8.0L) == 256.0); + assert(exp2(-9.0L)== 1.0L/512.0); + + static void testExp2(T)() + { + // NaN + const T specialNaN = NaN(0x0123L); + assert(isIdentical(exp2(specialNaN), specialNaN)); + + // over-/underflow + enum T OF = T.max_exp; + enum T UF = T.min_exp - T.mant_dig; + assert(isIdentical(exp2(OF + 1), cast(T) T.infinity)); + assert(isIdentical(exp2(UF - 1), cast(T) 0.0)); + + static immutable T[2][] vals = + [ + // x, exp2(x) + [ 0.0, 1.0 ], + [ -0.0, 1.0 ], + [ 0.5, SQRT2 ], + [ 8.0, 256.0 ], + [ -9.0, 1.0 / 512 ], + ]; + + foreach (ref val; vals) + { + const T x = val[0]; + const T r = val[1]; + const T e = exp2(x); + + //printf("exp2(%Lg) = %Lg, should be %Lg\n", cast(real) x, cast(real) e, cast(real) r); + assert(isClose(r, e)); + } + } + + import std.meta : AliasSeq; + foreach (T; AliasSeq!(real, double, float)) + testExp2!T(); +} + +/********************************************************************* + * Separate floating point value into significand and exponent. + * + * Returns: + * Calculate and return $(I x) and $(I exp) such that + * value =$(I x)*2$(SUPERSCRIPT exp) and + * .5 $(LT)= |$(I x)| $(LT) 1.0 + * + * $(I x) has same sign as value. + * + * $(TABLE_SV + * $(TR $(TH value) $(TH returns) $(TH exp)) + * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD 0)) + * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD int.max)) + * $(TR $(TD -$(INFIN)) $(TD -$(INFIN)) $(TD int.min)) + * $(TR $(TD $(PLUSMN)$(NAN)) $(TD $(PLUSMN)$(NAN)) $(TD int.min)) + * ) + */ +T frexp(T)(const T value, out int exp) @trusted pure nothrow @nogc +if (isFloatingPoint!T) +{ + import std.math : floatTraits, RealFormat, MANTISSA_MSB, MANTISSA_LSB; + import std.math.traits : isSubnormal; + + if (__ctfe) + { + // Handle special cases. + if (value == 0) { exp = 0; return value; } + if (value == T.infinity) { exp = int.max; return value; } + if (value == -T.infinity || value != value) { exp = int.min; return value; } + // Handle ordinary cases. + // In CTFE there is no performance advantage for having separate + // paths for different floating point types. + T absValue = value < 0 ? -value : value; + int expCount; + static if (T.mant_dig > double.mant_dig) + { + for (; absValue >= 0x1.0p+1024L; absValue *= 0x1.0p-1024L) + expCount += 1024; + for (; absValue < 0x1.0p-1021L; absValue *= 0x1.0p+1021L) + expCount -= 1021; + } + const double dval = cast(double) absValue; + int dexp = cast(int) (((*cast(const long*) &dval) >>> 52) & 0x7FF) + double.min_exp - 2; + dexp++; + expCount += dexp; + absValue *= 2.0 ^^ -dexp; + // If the original value was subnormal or if it was a real + // then absValue can still be outside the [0.5, 1.0) range. + if (absValue < 0.5) + { + assert(T.mant_dig > double.mant_dig || isSubnormal(value)); + do + { + absValue += absValue; + expCount--; + } while (absValue < 0.5); + } + else + { + assert(absValue < 1 || T.mant_dig > double.mant_dig); + for (; absValue >= 1; absValue *= T(0.5)) + expCount++; + } + exp = expCount; + return value < 0 ? -absValue : absValue; + } + + Unqual!T vf = value; + ushort* vu = cast(ushort*)&vf; + static if (is(immutable T == immutable float)) + int* vi = cast(int*)&vf; + else + long* vl = cast(long*)&vf; + int ex; + alias F = floatTraits!T; + + ex = vu[F.EXPPOS_SHORT] & F.EXPMASK; + static if (F.realFormat == RealFormat.ieeeExtended || + F.realFormat == RealFormat.ieeeExtended53) + { + if (ex) + { // If exponent is non-zero + if (ex == F.EXPMASK) // infinity or NaN + { + if (*vl & 0x7FFF_FFFF_FFFF_FFFF) // NaN + { + *vl |= 0xC000_0000_0000_0000; // convert NaNS to NaNQ + exp = int.min; + } + else if (vu[F.EXPPOS_SHORT] & 0x8000) // negative infinity + exp = int.min; + else // positive infinity + exp = int.max; + + } + else + { + exp = ex - F.EXPBIAS; + vu[F.EXPPOS_SHORT] = (0x8000 & vu[F.EXPPOS_SHORT]) | 0x3FFE; + } + } + else if (!*vl) + { + // vf is +-0.0 + exp = 0; + } + else + { + // subnormal + + vf *= F.RECIP_EPSILON; + ex = vu[F.EXPPOS_SHORT] & F.EXPMASK; + exp = ex - F.EXPBIAS - T.mant_dig + 1; + vu[F.EXPPOS_SHORT] = ((-1 - F.EXPMASK) & vu[F.EXPPOS_SHORT]) | 0x3FFE; + } + return vf; + } + else static if (F.realFormat == RealFormat.ieeeQuadruple) + { + if (ex) // If exponent is non-zero + { + if (ex == F.EXPMASK) + { + // infinity or NaN + if (vl[MANTISSA_LSB] | + (vl[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) // NaN + { + // convert NaNS to NaNQ + vl[MANTISSA_MSB] |= 0x0000_8000_0000_0000; + exp = int.min; + } + else if (vu[F.EXPPOS_SHORT] & 0x8000) // negative infinity + exp = int.min; + else // positive infinity + exp = int.max; + } + else + { + exp = ex - F.EXPBIAS; + vu[F.EXPPOS_SHORT] = F.EXPBIAS | (0x8000 & vu[F.EXPPOS_SHORT]); + } + } + else if ((vl[MANTISSA_LSB] | + (vl[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) == 0) + { + // vf is +-0.0 + exp = 0; + } + else + { + // subnormal + vf *= F.RECIP_EPSILON; + ex = vu[F.EXPPOS_SHORT] & F.EXPMASK; + exp = ex - F.EXPBIAS - T.mant_dig + 1; + vu[F.EXPPOS_SHORT] = F.EXPBIAS | (0x8000 & vu[F.EXPPOS_SHORT]); + } + return vf; + } + else static if (F.realFormat == RealFormat.ieeeDouble) + { + if (ex) // If exponent is non-zero + { + if (ex == F.EXPMASK) // infinity or NaN + { + if (*vl == 0x7FF0_0000_0000_0000) // positive infinity + { + exp = int.max; + } + else if (*vl == 0xFFF0_0000_0000_0000) // negative infinity + exp = int.min; + else + { // NaN + *vl |= 0x0008_0000_0000_0000; // convert NaNS to NaNQ + exp = int.min; + } + } + else + { + exp = (ex - F.EXPBIAS) >> 4; + vu[F.EXPPOS_SHORT] = cast(ushort)((0x800F & vu[F.EXPPOS_SHORT]) | 0x3FE0); + } + } + else if (!(*vl & 0x7FFF_FFFF_FFFF_FFFF)) + { + // vf is +-0.0 + exp = 0; + } + else + { + // subnormal + vf *= F.RECIP_EPSILON; + ex = vu[F.EXPPOS_SHORT] & F.EXPMASK; + exp = ((ex - F.EXPBIAS) >> 4) - T.mant_dig + 1; + vu[F.EXPPOS_SHORT] = + cast(ushort)(((-1 - F.EXPMASK) & vu[F.EXPPOS_SHORT]) | 0x3FE0); + } + return vf; + } + else static if (F.realFormat == RealFormat.ieeeSingle) + { + if (ex) // If exponent is non-zero + { + if (ex == F.EXPMASK) // infinity or NaN + { + if (*vi == 0x7F80_0000) // positive infinity + { + exp = int.max; + } + else if (*vi == 0xFF80_0000) // negative infinity + exp = int.min; + else + { // NaN + *vi |= 0x0040_0000; // convert NaNS to NaNQ + exp = int.min; + } + } + else + { + exp = (ex - F.EXPBIAS) >> 7; + vu[F.EXPPOS_SHORT] = cast(ushort)((0x807F & vu[F.EXPPOS_SHORT]) | 0x3F00); + } + } + else if (!(*vi & 0x7FFF_FFFF)) + { + // vf is +-0.0 + exp = 0; + } + else + { + // subnormal + vf *= F.RECIP_EPSILON; + ex = vu[F.EXPPOS_SHORT] & F.EXPMASK; + exp = ((ex - F.EXPBIAS) >> 7) - T.mant_dig + 1; + vu[F.EXPPOS_SHORT] = + cast(ushort)(((-1 - F.EXPMASK) & vu[F.EXPPOS_SHORT]) | 0x3F00); + } + return vf; + } + else // static if (F.realFormat == RealFormat.ibmExtended) + { + assert(0, "frexp not implemented"); + } +} + +/// +@safe unittest +{ + import std.math.operations : isClose; + + int exp; + real mantissa = frexp(123.456L, exp); + + assert(isClose(mantissa * pow(2.0L, cast(real) exp), 123.456L)); + + assert(frexp(-real.nan, exp) && exp == int.min); + assert(frexp(real.nan, exp) && exp == int.min); + assert(frexp(-real.infinity, exp) == -real.infinity && exp == int.min); + assert(frexp(real.infinity, exp) == real.infinity && exp == int.max); + assert(frexp(-0.0, exp) == -0.0 && exp == 0); + assert(frexp(0.0, exp) == 0.0 && exp == 0); +} + +@safe @nogc nothrow unittest +{ + import std.math.operations : isClose; + + int exp; + real mantissa = frexp(123.456L, exp); + + // check if values are equal to 19 decimal digits of precision + assert(isClose(mantissa * pow(2.0L, cast(real) exp), 123.456L, 1e-18)); +} + +@safe unittest +{ + import std.math : floatTraits, RealFormat; + import std.math.traits : isIdentical; + import std.meta : AliasSeq; + import std.typecons : tuple, Tuple; + + static foreach (T; AliasSeq!(real, double, float)) + {{ + Tuple!(T, T, int)[] vals = [ // x,frexp,exp + tuple(T(0.0), T( 0.0 ), 0), + tuple(T(-0.0), T( -0.0), 0), + tuple(T(1.0), T( .5 ), 1), + tuple(T(-1.0), T( -.5 ), 1), + tuple(T(2.0), T( .5 ), 2), + tuple(T(float.min_normal/2.0f), T(.5), -126), + tuple(T.infinity, T.infinity, int.max), + tuple(-T.infinity, -T.infinity, int.min), + tuple(T.nan, T.nan, int.min), + tuple(-T.nan, -T.nan, int.min), + + // https://issues.dlang.org/show_bug.cgi?id=16026: + tuple(3 * (T.min_normal * T.epsilon), T( .75), (T.min_exp - T.mant_dig) + 2) + ]; + + foreach (elem; vals) + { + T x = elem[0]; + T e = elem[1]; + int exp = elem[2]; + int eptr; + T v = frexp(x, eptr); + assert(isIdentical(e, v)); + assert(exp == eptr); + } + + static if (floatTraits!(T).realFormat == RealFormat.ieeeExtended) + { + static T[3][] extendedvals = [ // x,frexp,exp + [0x1.a5f1c2eb3fe4efp+73L, 0x1.A5F1C2EB3FE4EFp-1L, 74], // normal + [0x1.fa01712e8f0471ap-1064L, 0x1.fa01712e8f0471ap-1L, -1063], + [T.min_normal, .5, -16381], + [T.min_normal/2.0L, .5, -16382] // subnormal + ]; + foreach (elem; extendedvals) + { + T x = elem[0]; + T e = elem[1]; + int exp = cast(int) elem[2]; + int eptr; + T v = frexp(x, eptr); + assert(isIdentical(e, v)); + assert(exp == eptr); + } + } + }} + + // CTFE + alias CtfeFrexpResult= Tuple!(real, int); + static CtfeFrexpResult ctfeFrexp(T)(const T value) + { + int exp; + auto significand = frexp(value, exp); + return CtfeFrexpResult(significand, exp); + } + static foreach (T; AliasSeq!(real, double, float)) + {{ + enum Tuple!(T, T, int)[] vals = [ // x,frexp,exp + tuple(T(0.0), T( 0.0 ), 0), + tuple(T(-0.0), T( -0.0), 0), + tuple(T(1.0), T( .5 ), 1), + tuple(T(-1.0), T( -.5 ), 1), + tuple(T(2.0), T( .5 ), 2), + tuple(T(float.min_normal/2.0f), T(.5), -126), + tuple(T.infinity, T.infinity, int.max), + tuple(-T.infinity, -T.infinity, int.min), + tuple(T.nan, T.nan, int.min), + tuple(-T.nan, -T.nan, int.min), + + // https://issues.dlang.org/show_bug.cgi?id=16026: + tuple(3 * (T.min_normal * T.epsilon), T( .75), (T.min_exp - T.mant_dig) + 2) + ]; + + static foreach (elem; vals) + { + static assert(ctfeFrexp(elem[0]) is CtfeFrexpResult(elem[1], elem[2])); + } + + static if (floatTraits!(T).realFormat == RealFormat.ieeeExtended) + { + enum T[3][] extendedvals = [ // x,frexp,exp + [0x1.a5f1c2eb3fe4efp+73L, 0x1.A5F1C2EB3FE4EFp-1L, 74], // normal + [0x1.fa01712e8f0471ap-1064L, 0x1.fa01712e8f0471ap-1L, -1063], + [T.min_normal, .5, -16381], + [T.min_normal/2.0L, .5, -16382] // subnormal + ]; + static foreach (elem; extendedvals) + { + static assert(ctfeFrexp(elem[0]) is CtfeFrexpResult(elem[1], cast(int) elem[2])); + } + } + }} +} + +@safe unittest +{ + import std.meta : AliasSeq; + void foo() { + static foreach (T; AliasSeq!(real, double, float)) + {{ + int exp; + const T a = 1; + immutable T b = 2; + auto c = frexp(a, exp); + auto d = frexp(b, exp); + }} + } +} + +/****************************************** + * Extracts the exponent of x as a signed integral value. + * + * If x is not a special value, the result is the same as + * $(D cast(int) logb(x)). + * + * $(TABLE_SV + * $(TR $(TH x) $(TH ilogb(x)) $(TH Range error?)) + * $(TR $(TD 0) $(TD FP_ILOGB0) $(TD yes)) + * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD int.max) $(TD no)) + * $(TR $(TD $(NAN)) $(TD FP_ILOGBNAN) $(TD no)) + * ) + */ +int ilogb(T)(const T x) @trusted pure nothrow @nogc +if (isFloatingPoint!T) +{ + import std.math : floatTraits, RealFormat, MANTISSA_MSB, MANTISSA_LSB; + + import core.bitop : bsr; + alias F = floatTraits!T; + + union floatBits + { + T rv; + ushort[T.sizeof/2] vu; + uint[T.sizeof/4] vui; + static if (T.sizeof >= 8) + ulong[T.sizeof/8] vul; + } + floatBits y = void; + y.rv = x; + + int ex = y.vu[F.EXPPOS_SHORT] & F.EXPMASK; + static if (F.realFormat == RealFormat.ieeeExtended || + F.realFormat == RealFormat.ieeeExtended53) + { + if (ex) + { + // If exponent is non-zero + if (ex == F.EXPMASK) // infinity or NaN + { + if (y.vul[0] & 0x7FFF_FFFF_FFFF_FFFF) // NaN + return FP_ILOGBNAN; + else // +-infinity + return int.max; + } + else + { + return ex - F.EXPBIAS - 1; + } + } + else if (!y.vul[0]) + { + // vf is +-0.0 + return FP_ILOGB0; + } + else + { + // subnormal + return ex - F.EXPBIAS - T.mant_dig + 1 + bsr(y.vul[0]); + } + } + else static if (F.realFormat == RealFormat.ieeeQuadruple) + { + if (ex) // If exponent is non-zero + { + if (ex == F.EXPMASK) + { + // infinity or NaN + if (y.vul[MANTISSA_LSB] | ( y.vul[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) // NaN + return FP_ILOGBNAN; + else // +- infinity + return int.max; + } + else + { + return ex - F.EXPBIAS - 1; + } + } + else if ((y.vul[MANTISSA_LSB] | (y.vul[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) == 0) + { + // vf is +-0.0 + return FP_ILOGB0; + } + else + { + // subnormal + const ulong msb = y.vul[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF; + const ulong lsb = y.vul[MANTISSA_LSB]; + if (msb) + return ex - F.EXPBIAS - T.mant_dig + 1 + bsr(msb) + 64; + else + return ex - F.EXPBIAS - T.mant_dig + 1 + bsr(lsb); + } + } + else static if (F.realFormat == RealFormat.ieeeDouble) + { + if (ex) // If exponent is non-zero + { + if (ex == F.EXPMASK) // infinity or NaN + { + if ((y.vul[0] & 0x7FFF_FFFF_FFFF_FFFF) == 0x7FF0_0000_0000_0000) // +- infinity + return int.max; + else // NaN + return FP_ILOGBNAN; + } + else + { + return ((ex - F.EXPBIAS) >> 4) - 1; + } + } + else if (!(y.vul[0] & 0x7FFF_FFFF_FFFF_FFFF)) + { + // vf is +-0.0 + return FP_ILOGB0; + } + else + { + // subnormal + enum MANTISSAMASK_64 = ((cast(ulong) F.MANTISSAMASK_INT) << 32) | 0xFFFF_FFFF; + return ((ex - F.EXPBIAS) >> 4) - T.mant_dig + 1 + bsr(y.vul[0] & MANTISSAMASK_64); + } + } + else static if (F.realFormat == RealFormat.ieeeSingle) + { + if (ex) // If exponent is non-zero + { + if (ex == F.EXPMASK) // infinity or NaN + { + if ((y.vui[0] & 0x7FFF_FFFF) == 0x7F80_0000) // +- infinity + return int.max; + else // NaN + return FP_ILOGBNAN; + } + else + { + return ((ex - F.EXPBIAS) >> 7) - 1; + } + } + else if (!(y.vui[0] & 0x7FFF_FFFF)) + { + // vf is +-0.0 + return FP_ILOGB0; + } + else + { + // subnormal + const uint mantissa = y.vui[0] & F.MANTISSAMASK_INT; + return ((ex - F.EXPBIAS) >> 7) - T.mant_dig + 1 + bsr(mantissa); + } + } + else // static if (F.realFormat == RealFormat.ibmExtended) + { + assert(0, "ilogb not implemented"); + } +} +/// ditto +int ilogb(T)(const T x) @safe pure nothrow @nogc +if (isIntegral!T && isUnsigned!T) +{ + import core.bitop : bsr; + if (x == 0) + return FP_ILOGB0; + else + { + static assert(T.sizeof <= ulong.sizeof, "integer size too large for the current ilogb implementation"); + return bsr(x); + } +} +/// ditto +int ilogb(T)(const T x) @safe pure nothrow @nogc +if (isIntegral!T && isSigned!T) +{ + import std.traits : Unsigned; + // Note: abs(x) can not be used because the return type is not Unsigned and + // the return value would be wrong for x == int.min + Unsigned!T absx = x >= 0 ? x : -x; + return ilogb(absx); +} + +/// +@safe pure unittest +{ + assert(ilogb(1) == 0); + assert(ilogb(3) == 1); + assert(ilogb(3.0) == 1); + assert(ilogb(100_000_000) == 26); + + assert(ilogb(0) == FP_ILOGB0); + assert(ilogb(0.0) == FP_ILOGB0); + assert(ilogb(double.nan) == FP_ILOGBNAN); + assert(ilogb(double.infinity) == int.max); +} + +/** +Special return values of $(LREF ilogb). + */ +alias FP_ILOGB0 = core.stdc.math.FP_ILOGB0; +/// ditto +alias FP_ILOGBNAN = core.stdc.math.FP_ILOGBNAN; + +/// +@safe pure unittest +{ + assert(ilogb(0) == FP_ILOGB0); + assert(ilogb(0.0) == FP_ILOGB0); + assert(ilogb(double.nan) == FP_ILOGBNAN); +} + +@safe nothrow @nogc unittest +{ + import std.math : floatTraits, RealFormat; + import std.math.operations : nextUp; + import std.meta : AliasSeq; + import std.typecons : Tuple; + static foreach (F; AliasSeq!(float, double, real)) + {{ + alias T = Tuple!(F, int); + T[13] vals = // x, ilogb(x) + [ + T( F.nan , FP_ILOGBNAN ), + T( -F.nan , FP_ILOGBNAN ), + T( F.infinity, int.max ), + T( -F.infinity, int.max ), + T( 0.0 , FP_ILOGB0 ), + T( -0.0 , FP_ILOGB0 ), + T( 2.0 , 1 ), + T( 2.0001 , 1 ), + T( 1.9999 , 0 ), + T( 0.5 , -1 ), + T( 123.123 , 6 ), + T( -123.123 , 6 ), + T( 0.123 , -4 ), + ]; + + foreach (elem; vals) + { + assert(ilogb(elem[0]) == elem[1]); + } + }} + + // min_normal and subnormals + assert(ilogb(-float.min_normal) == -126); + assert(ilogb(nextUp(-float.min_normal)) == -127); + assert(ilogb(nextUp(-float(0.0))) == -149); + assert(ilogb(-double.min_normal) == -1022); + assert(ilogb(nextUp(-double.min_normal)) == -1023); + assert(ilogb(nextUp(-double(0.0))) == -1074); + static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended) + { + assert(ilogb(-real.min_normal) == -16382); + assert(ilogb(nextUp(-real.min_normal)) == -16383); + assert(ilogb(nextUp(-real(0.0))) == -16445); + } + else static if (floatTraits!(real).realFormat == RealFormat.ieeeDouble) + { + assert(ilogb(-real.min_normal) == -1022); + assert(ilogb(nextUp(-real.min_normal)) == -1023); + assert(ilogb(nextUp(-real(0.0))) == -1074); + } + + // test integer types + assert(ilogb(0) == FP_ILOGB0); + assert(ilogb(int.max) == 30); + assert(ilogb(int.min) == 31); + assert(ilogb(uint.max) == 31); + assert(ilogb(long.max) == 62); + assert(ilogb(long.min) == 63); + assert(ilogb(ulong.max) == 63); +} + +/******************************************* + * Compute n * 2$(SUPERSCRIPT exp) + * References: frexp + */ + +pragma(inline, true) +real ldexp(real n, int exp) @safe pure nothrow @nogc { return core.math.ldexp(n, exp); } +///ditto +pragma(inline, true) +double ldexp(double n, int exp) @safe pure nothrow @nogc { return core.math.ldexp(n, exp); } +///ditto +pragma(inline, true) +float ldexp(float n, int exp) @safe pure nothrow @nogc { return core.math.ldexp(n, exp); } + +/// +@nogc @safe pure nothrow unittest +{ + import std.meta : AliasSeq; + static foreach (T; AliasSeq!(float, double, real)) + {{ + T r; + + r = ldexp(3.0L, 3); + assert(r == 24); + + r = ldexp(cast(T) 3.0, cast(int) 3); + assert(r == 24); + + T n = 3.0; + int exp = 3; + r = ldexp(n, exp); + assert(r == 24); + }} +} + +@safe pure nothrow @nogc unittest +{ + import std.math : floatTraits, RealFormat; + + static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended || + floatTraits!(real).realFormat == RealFormat.ieeeExtended53 || + floatTraits!(real).realFormat == RealFormat.ieeeQuadruple) + { + assert(ldexp(1.0L, -16384) == 0x1p-16384L); + assert(ldexp(1.0L, -16382) == 0x1p-16382L); + int x; + real n = frexp(0x1p-16384L, x); + assert(n == 0.5L); + assert(x==-16383); + assert(ldexp(n, x)==0x1p-16384L); + } + else static if (floatTraits!(real).realFormat == RealFormat.ieeeDouble) + { + assert(ldexp(1.0L, -1024) == 0x1p-1024L); + assert(ldexp(1.0L, -1022) == 0x1p-1022L); + int x; + real n = frexp(0x1p-1024L, x); + assert(n == 0.5L); + assert(x==-1023); + assert(ldexp(n, x)==0x1p-1024L); + } + else static assert(false, "Floating point type real not supported"); +} + +/* workaround https://issues.dlang.org/show_bug.cgi?id=14718 + float parsing depends on platform strtold +@safe pure nothrow @nogc unittest +{ + assert(ldexp(1.0, -1024) == 0x1p-1024); + assert(ldexp(1.0, -1022) == 0x1p-1022); + int x; + double n = frexp(0x1p-1024, x); + assert(n == 0.5); + assert(x==-1023); + assert(ldexp(n, x)==0x1p-1024); +} + +@safe pure nothrow @nogc unittest +{ + assert(ldexp(1.0f, -128) == 0x1p-128f); + assert(ldexp(1.0f, -126) == 0x1p-126f); + int x; + float n = frexp(0x1p-128f, x); + assert(n == 0.5f); + assert(x==-127); + assert(ldexp(n, x)==0x1p-128f); +} +*/ + +@safe @nogc nothrow unittest +{ + import std.math.operations : isClose; + + static real[3][] vals = // value,exp,ldexp + [ + [ 0, 0, 0], + [ 1, 0, 1], + [ -1, 0, -1], + [ 1, 1, 2], + [ 123, 10, 125952], + [ real.max, int.max, real.infinity], + [ real.max, -int.max, 0], + [ real.min_normal, -int.max, 0], + ]; + int i; + + for (i = 0; i < vals.length; i++) + { + real x = vals[i][0]; + int exp = cast(int) vals[i][1]; + real z = vals[i][2]; + real l = ldexp(x, exp); + + assert(isClose(z, l, 1e-6)); + } + + real function(real, int) pldexp = &ldexp; + assert(pldexp != null); +} + +private +{ + import std.math : floatTraits, RealFormat; + + version (INLINE_YL2X) {} else + { + static if (floatTraits!real.realFormat == RealFormat.ieeeQuadruple) + { + // Coefficients for log(1 + x) = x - x**2/2 + x**3 P(x)/Q(x) + static immutable real[13] logCoeffsP = [ + 1.313572404063446165910279910527789794488E4L, + 7.771154681358524243729929227226708890930E4L, + 2.014652742082537582487669938141683759923E5L, + 3.007007295140399532324943111654767187848E5L, + 2.854829159639697837788887080758954924001E5L, + 1.797628303815655343403735250238293741397E5L, + 7.594356839258970405033155585486712125861E4L, + 2.128857716871515081352991964243375186031E4L, + 3.824952356185897735160588078446136783779E3L, + 4.114517881637811823002128927449878962058E2L, + 2.321125933898420063925789532045674660756E1L, + 4.998469661968096229986658302195402690910E-1L, + 1.538612243596254322971797716843006400388E-6L + ]; + static immutable real[13] logCoeffsQ = [ + 3.940717212190338497730839731583397586124E4L, + 2.626900195321832660448791748036714883242E5L, + 7.777690340007566932935753241556479363645E5L, + 1.347518538384329112529391120390701166528E6L, + 1.514882452993549494932585972882995548426E6L, + 1.158019977462989115839826904108208787040E6L, + 6.132189329546557743179177159925690841200E5L, + 2.248234257620569139969141618556349415120E5L, + 5.605842085972455027590989944010492125825E4L, + 9.147150349299596453976674231612674085381E3L, + 9.104928120962988414618126155557301584078E2L, + 4.839208193348159620282142911143429644326E1L, + 1.0 + ]; + + // Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2) + // where z = 2(x-1)/(x+1) + static immutable real[6] logCoeffsR = [ + 1.418134209872192732479751274970992665513E5L, + -8.977257995689735303686582344659576526998E4L, + 2.048819892795278657810231591630928516206E4L, + -2.024301798136027039250415126250455056397E3L, + 8.057002716646055371965756206836056074715E1L, + -8.828896441624934385266096344596648080902E-1L + ]; + static immutable real[7] logCoeffsS = [ + 1.701761051846631278975701529965589676574E6L, + -1.332535117259762928288745111081235577029E6L, + 4.001557694070773974936904547424676279307E5L, + -5.748542087379434595104154610899551484314E4L, + 3.998526750980007367835804959888064681098E3L, + -1.186359407982897997337150403816839480438E2L, + 1.0 + ]; + } + else + { + // Coefficients for log(1 + x) = x - x**2/2 + x**3 P(x)/Q(x) + static immutable real[7] logCoeffsP = [ + 2.0039553499201281259648E1L, + 5.7112963590585538103336E1L, + 6.0949667980987787057556E1L, + 2.9911919328553073277375E1L, + 6.5787325942061044846969E0L, + 4.9854102823193375972212E-1L, + 4.5270000862445199635215E-5L, + ]; + static immutable real[7] logCoeffsQ = [ + 6.0118660497603843919306E1L, + 2.1642788614495947685003E2L, + 3.0909872225312059774938E2L, + 2.2176239823732856465394E2L, + 8.3047565967967209469434E1L, + 1.5062909083469192043167E1L, + 1.0000000000000000000000E0L, + ]; + + // Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2) + // where z = 2(x-1)/(x+1) + static immutable real[4] logCoeffsR = [ + -3.5717684488096787370998E1L, + 1.0777257190312272158094E1L, + -7.1990767473014147232598E-1L, + 1.9757429581415468984296E-3L, + ]; + static immutable real[4] logCoeffsS = [ + -4.2861221385716144629696E2L, + 1.9361891836232102174846E2L, + -2.6201045551331104417768E1L, + 1.0000000000000000000000E0L, + ]; + } + } +} + +/************************************** + * Calculate the natural logarithm of x. + * + * $(TABLE_SV + * $(TR $(TH x) $(TH log(x)) $(TH divide by 0?) $(TH invalid?)) + * $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes) $(TD no)) + * $(TR $(TD $(LT)0.0) $(TD $(NAN)) $(TD no) $(TD yes)) + * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no) $(TD no)) + * ) + */ +real log(real x) @safe pure nothrow @nogc +{ + import std.math.constants : LN2, LOG2, SQRT1_2; + import std.math.traits : isInfinity, isNaN, signbit; + import std.math.algebraic : poly; + + version (INLINE_YL2X) + return core.math.yl2x(x, LN2); + else + { + // C1 + C2 = LN2. + enum real C1 = 6.93145751953125E-1L; + enum real C2 = 1.428606820309417232121458176568075500134E-6L; + + // Special cases. + if (isNaN(x)) + return x; + if (isInfinity(x) && !signbit(x)) + return x; + if (x == 0.0) + return -real.infinity; + if (x < 0.0) + return real.nan; + + // Separate mantissa from exponent. + // Note, frexp is used so that denormal numbers will be handled properly. + real y, z; + int exp; + + x = frexp(x, exp); + + // Logarithm using log(x) = z + z^^3 R(z) / S(z), + // where z = 2(x - 1)/(x + 1) + if ((exp > 2) || (exp < -2)) + { + if (x < SQRT1_2) + { // 2(2x - 1)/(2x + 1) + exp -= 1; + z = x - 0.5; + y = 0.5 * z + 0.5; + } + else + { // 2(x - 1)/(x + 1) + z = x - 0.5; + z -= 0.5; + y = 0.5 * x + 0.5; + } + x = z / y; + z = x * x; + z = x * (z * poly(z, logCoeffsR) / poly(z, logCoeffsS)); + z += exp * C2; + z += x; + z += exp * C1; + + return z; + } + + // Logarithm using log(1 + x) = x - .5x^^2 + x^^3 P(x) / Q(x) + if (x < SQRT1_2) + { + exp -= 1; + x = 2.0 * x - 1.0; + } + else + { + x = x - 1.0; + } + z = x * x; + y = x * (z * poly(x, logCoeffsP) / poly(x, logCoeffsQ)); + y += exp * C2; + z = y - 0.5 * z; + + // Note, the sum of above terms does not exceed x/4, + // so it contributes at most about 1/4 lsb to the error. + z += x; + z += exp * C1; + + return z; + } +} + +/// +@safe pure nothrow @nogc unittest +{ + import std.math.operations : feqrel; + import std.math.constants : E; + + assert(feqrel(log(E), 1) >= real.mant_dig - 1); +} + +/************************************** + * Calculate the base-10 logarithm of x. + * + * $(TABLE_SV + * $(TR $(TH x) $(TH log10(x)) $(TH divide by 0?) $(TH invalid?)) + * $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes) $(TD no)) + * $(TR $(TD $(LT)0.0) $(TD $(NAN)) $(TD no) $(TD yes)) + * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no) $(TD no)) + * ) + */ +real log10(real x) @safe pure nothrow @nogc +{ + import std.math.constants : LOG2, LN2, SQRT1_2; + import std.math.algebraic : poly; + import std.math.traits : isNaN, isInfinity, signbit; + + version (INLINE_YL2X) + return core.math.yl2x(x, LOG2); + else + { + // log10(2) split into two parts. + enum real L102A = 0.3125L; + enum real L102B = -1.14700043360188047862611052755069732318101185E-2L; + + // log10(e) split into two parts. + enum real L10EA = 0.5L; + enum real L10EB = -6.570551809674817234887108108339491770560299E-2L; + + // Special cases are the same as for log. + if (isNaN(x)) + return x; + if (isInfinity(x) && !signbit(x)) + return x; + if (x == 0.0) + return -real.infinity; + if (x < 0.0) + return real.nan; + + // Separate mantissa from exponent. + // Note, frexp is used so that denormal numbers will be handled properly. + real y, z; + int exp; + + x = frexp(x, exp); + + // Logarithm using log(x) = z + z^^3 R(z) / S(z), + // where z = 2(x - 1)/(x + 1) + if ((exp > 2) || (exp < -2)) + { + if (x < SQRT1_2) + { // 2(2x - 1)/(2x + 1) + exp -= 1; + z = x - 0.5; + y = 0.5 * z + 0.5; + } + else + { // 2(x - 1)/(x + 1) + z = x - 0.5; + z -= 0.5; + y = 0.5 * x + 0.5; + } + x = z / y; + z = x * x; + y = x * (z * poly(z, logCoeffsR) / poly(z, logCoeffsS)); + goto Ldone; + } + + // Logarithm using log(1 + x) = x - .5x^^2 + x^^3 P(x) / Q(x) + if (x < SQRT1_2) + { + exp -= 1; + x = 2.0 * x - 1.0; + } + else + x = x - 1.0; + + z = x * x; + y = x * (z * poly(x, logCoeffsP) / poly(x, logCoeffsQ)); + y = y - 0.5 * z; + + // Multiply log of fraction by log10(e) and base 2 exponent by log10(2). + // This sequence of operations is critical and it may be horribly + // defeated by some compiler optimizers. + Ldone: + z = y * L10EB; + z += x * L10EB; + z += exp * L102B; + z += y * L10EA; + z += x * L10EA; + z += exp * L102A; + + return z; + } +} + +/// +@safe pure nothrow @nogc unittest +{ + import std.math.algebraic : fabs; + + assert(fabs(log10(1000) - 3) < .000001); +} + +/** + * Calculates the natural logarithm of 1 + x. + * + * For very small x, log1p(x) will be more accurate than + * log(1 + x). + * + * $(TABLE_SV + * $(TR $(TH x) $(TH log1p(x)) $(TH divide by 0?) $(TH invalid?)) + * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no) $(TD no)) + * $(TR $(TD -1.0) $(TD -$(INFIN)) $(TD yes) $(TD no)) + * $(TR $(TD $(LT)-1.0) $(TD -$(NAN)) $(TD no) $(TD yes)) + * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no) $(TD no)) + * ) + */ +real log1p(real x) @safe pure nothrow @nogc +{ + import std.math.traits : isNaN, isInfinity, signbit; + import std.math.constants : LN2; + + version (INLINE_YL2X) + { + // On x87, yl2xp1 is valid if and only if -0.5 <= lg(x) <= 0.5, + // ie if -0.29 <= x <= 0.414 + return (core.math.fabs(x) <= 0.25) ? core.math.yl2xp1(x, LN2) : core.math.yl2x(x+1, LN2); + } + else + { + // Special cases. + if (isNaN(x) || x == 0.0) + return x; + if (isInfinity(x) && !signbit(x)) + return x; + if (x == -1.0) + return -real.infinity; + if (x < -1.0) + return real.nan; + + return log(x + 1.0); + } +} + +/// +@safe pure unittest +{ + import std.math.traits : isIdentical, isNaN; + import std.math.operations : feqrel; + + assert(isIdentical(log1p(0.0), 0.0)); + assert(log1p(1.0).feqrel(0.69314) > 16); + + assert(log1p(-1.0) == -real.infinity); + assert(isNaN(log1p(-2.0))); + assert(log1p(real.nan) is real.nan); + assert(log1p(-real.nan) is -real.nan); + assert(log1p(real.infinity) == real.infinity); +} + +/*************************************** + * Calculates the base-2 logarithm of x: + * $(SUB log, 2)x + * + * $(TABLE_SV + * $(TR $(TH x) $(TH log2(x)) $(TH divide by 0?) $(TH invalid?)) + * $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes) $(TD no) ) + * $(TR $(TD $(LT)0.0) $(TD $(NAN)) $(TD no) $(TD yes) ) + * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no) $(TD no) ) + * ) + */ +real log2(real x) @safe pure nothrow @nogc +{ + import std.math.traits : isNaN, isInfinity, signbit; + import std.math.constants : SQRT1_2, LOG2E; + import std.math.algebraic : poly; + + version (INLINE_YL2X) + return core.math.yl2x(x, 1.0L); + else + { + // Special cases are the same as for log. + if (isNaN(x)) + return x; + if (isInfinity(x) && !signbit(x)) + return x; + if (x == 0.0) + return -real.infinity; + if (x < 0.0) + return real.nan; + + // Separate mantissa from exponent. + // Note, frexp is used so that denormal numbers will be handled properly. + real y, z; + int exp; + + x = frexp(x, exp); + + // Logarithm using log(x) = z + z^^3 R(z) / S(z), + // where z = 2(x - 1)/(x + 1) + if ((exp > 2) || (exp < -2)) + { + if (x < SQRT1_2) + { // 2(2x - 1)/(2x + 1) + exp -= 1; + z = x - 0.5; + y = 0.5 * z + 0.5; + } + else + { // 2(x - 1)/(x + 1) + z = x - 0.5; + z -= 0.5; + y = 0.5 * x + 0.5; + } + x = z / y; + z = x * x; + y = x * (z * poly(z, logCoeffsR) / poly(z, logCoeffsS)); + goto Ldone; + } + + // Logarithm using log(1 + x) = x - .5x^^2 + x^^3 P(x) / Q(x) + if (x < SQRT1_2) + { + exp -= 1; + x = 2.0 * x - 1.0; + } + else + x = x - 1.0; + + z = x * x; + y = x * (z * poly(x, logCoeffsP) / poly(x, logCoeffsQ)); + y = y - 0.5 * z; + + // Multiply log of fraction by log10(e) and base 2 exponent by log10(2). + // This sequence of operations is critical and it may be horribly + // defeated by some compiler optimizers. + Ldone: + z = y * (LOG2E - 1.0); + z += x * (LOG2E - 1.0); + z += y; + z += x; + z += exp; + + return z; + } +} + +/// +@safe unittest +{ + import std.math.operations : isClose; + + assert(isClose(log2(1024.0L), 10)); +} + +@safe @nogc nothrow unittest +{ + import std.math.operations : isClose; + + // check if values are equal to 19 decimal digits of precision + assert(isClose(log2(1024.0L), 10, 1e-18)); +} + +/***************************************** + * Extracts the exponent of x as a signed integral value. + * + * If x is subnormal, it is treated as if it were normalized. + * For a positive, finite x: + * + * 1 $(LT)= $(I x) * FLT_RADIX$(SUPERSCRIPT -logb(x)) $(LT) FLT_RADIX + * + * $(TABLE_SV + * $(TR $(TH x) $(TH logb(x)) $(TH divide by 0?) ) + * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD +$(INFIN)) $(TD no)) + * $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes) ) + * ) + */ +real logb(real x) @trusted nothrow @nogc +{ + version (InlineAsm_X87_MSVC) + { + version (X86_64) + { + asm pure nothrow @nogc + { + naked ; + fld real ptr [RCX] ; + fxtract ; + fstp ST(0) ; + ret ; + } + } + else + { + asm pure nothrow @nogc + { + fld x ; + fxtract ; + fstp ST(0) ; + } + } + } + else + return core.stdc.math.logbl(x); +} + +/// +@safe @nogc nothrow unittest +{ + assert(logb(1.0) == 0); + assert(logb(100.0) == 6); + + assert(logb(0.0) == -real.infinity); + assert(logb(real.infinity) == real.infinity); + assert(logb(-real.infinity) == real.infinity); +} + +/************************************* + * Efficiently calculates x * 2$(SUPERSCRIPT n). + * + * scalbn handles underflow and overflow in + * the same fashion as the basic arithmetic operators. + * + * $(TABLE_SV + * $(TR $(TH x) $(TH scalb(x))) + * $(TR $(TD $(PLUSMNINF)) $(TD $(PLUSMNINF)) ) + * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) ) + * ) + */ +pragma(inline, true) +real scalbn(real x, int n) @safe pure nothrow @nogc { return _scalbn(x,n); } + +/// ditto +pragma(inline, true) +double scalbn(double x, int n) @safe pure nothrow @nogc { return _scalbn(x,n); } + +/// ditto +pragma(inline, true) +float scalbn(float x, int n) @safe pure nothrow @nogc { return _scalbn(x,n); } + +/// +@safe pure nothrow @nogc unittest +{ + assert(scalbn(0x1.2345678abcdefp0L, 999) == 0x1.2345678abcdefp999L); + assert(scalbn(-real.infinity, 5) == -real.infinity); + assert(scalbn(2.0,10) == 2048.0); + assert(scalbn(2048.0f,-10) == 2.0f); +} + +pragma(inline, true) +private F _scalbn(F)(F x, int n) +{ + import std.math.traits : isInfinity; + + if (__ctfe) + { + // Handle special cases. + if (x == F(0.0) || isInfinity(x)) + return x; + } + return core.math.ldexp(x, n); +} + +@safe pure nothrow @nogc unittest +{ + // CTFE-able test + static assert(scalbn(0x1.2345678abcdefp0L, 999) == 0x1.2345678abcdefp999L); + static assert(scalbn(-real.infinity, 5) == -real.infinity); + // Test with large exponent delta n where the result is in bounds but 2.0L ^^ n is not. + enum initialExponent = real.min_exp + 2, resultExponent = real.max_exp - 2; + enum int n = resultExponent - initialExponent; + enum real x = 0x1.2345678abcdefp0L * (2.0L ^^ initialExponent); + enum staticResult = scalbn(x, n); + static assert(staticResult == 0x1.2345678abcdefp0L * (2.0L ^^ resultExponent)); + assert(scalbn(x, n) == staticResult); +} + |