diff --git a/std/math.d b/std/math.d index 3b6d161acba..cf50a1d1ce4 100644 --- a/std/math.d +++ b/std/math.d @@ -6837,12 +6837,30 @@ Unqual!(Largest!(F, G)) pow(F, G)(F x, G y) @nogc @trusted pure nothrow { // Result is real only if y is an integer // Check for a non-zero fractional part - if (y > -1.0 / real.epsilon && y < 1.0 / real.epsilon) + enum maxOdd = pow(2.0L, real.mant_dig) - 1.0L; + static if (maxOdd > ulong.max) { - long w = cast(long)y; - if (w != y) + // Generic method, for any FP type + if (floor(y) != y) return sqrt(x); // Complex result -- create a NaN - if (w & 1) sign = -1.0; + + const hy = ldexp(y, -1); + 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; } @@ -6920,6 +6938,15 @@ Unqual!(Largest!(F, G)) pow(F, G)(F x, G y) @nogc @trusted pure nothrow assert(isIdentical(pow(0.0, 6.0), 0.0)); assert(isIdentical(pow(-0.0, 6.0), 0.0)); + // Issue #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(approxEqual(pow(two, three), 8.0)); assert(approxEqual(pow(two, -2.5), 0.1767767));