Improve pow and powi correctness and compile-time support

dub.sdl

@@ -10,6 +10,11 @@ buildType "unittest" {
     versions "mir_core_test"
 }
 
+buildType "unittest-cov" {
+    buildOptions "unittests" "coverage" "debugMode" "debugInfo"
+    versions "mir_core_test"
+}
+
 buildType "unittest-verbose" {
     buildOptions "unittests" "debugMode" "debugInfo"
     versions "mir_core_test"

source/mir/math/common.d

@@ -3,12 +3,59 @@ Common floating point math functions.
 
 This module has generic LLVM-oriented API compatible with all D compilers.
 
+$(D pow) follows libm special-value conventions: a base of 1 yields 1 even
+for a NaN exponent, and either sign of 1 raised to an infinite exponent
+yields 1. For zero or infinite bases, a negative result requires a negative
+base and an odd integer exponent; the sign of zero is preserved accordingly.
+
 License: $(HTTP www.apache.org/licenses/LICENSE-2.0, Apache-2.0)
 Authors:   Ilia Ki, Phobos Team
 +/
 module mir.math.common;
 
 import mir.internal.utility: isComplex, isComplexOf, isFloatingPoint;
+private import std.traits: Unqual;
+
+// Normalize special values before using Phobos, whose runtime pow can differ
+// from LLVM/libm and the compiler's CTFE builtin for these inputs.
+private bool powSpecial(T)(in T x, in T power, out Unqual!T result)
+    @safe pure nothrow @nogc
+    if (isFloatingPoint!T)
+{
+    import std.math: signbit;
+
+    if (power == 0 || x == 1)
+    {
+        result = 1;
+        return true;
+    }
+    if (power == T.infinity || power == -T.infinity)
+    {
+        auto magnitude = x < 0 ? -x : x;
+        if (magnitude == 1)
+            result = 1;
+        else if (magnitude > 1)
+            result = power > 0 ? T.infinity : 0;
+        else if (magnitude < 1)
+            result = power > 0 ? 0 : T.infinity;
+        else
+            return false; // NaN base: let the backend propagate it.
+        return true;
+    }
+    if ((x == 0 || x == T.infinity || x == -T.infinity)
+        && (power > 0 || power < 0))
+    {
+        // Positive powers keep zero at zero and infinity at infinity;
+        // negative powers exchange them.
+        result = (x == 0) == (power > 0) ? 0 : T.infinity;
+        // Avoid integer casts: real can represent odd integers beyond long.max.
+        auto remainder = power % 2;
+        if (signbit(x) && (remainder == 1 || remainder == -1))
+            result = -result;
+        return true;
+    }
+    return false;
+}
 
 version(LDC)
 {
@@ -126,7 +173,11 @@ enum
 
 version(LDC)
 {
-    import ldc.intrinsics: LLVM_version;
+    static import ldc.intrinsics;
+    static if (__traits(hasMember, ldc.intrinsics, "LLVM_major"))
+        private enum llvmMajor = ldc.intrinsics.LLVM_major;
+    else
+        private enum llvmMajor = ldc.intrinsics.LLVM_version / 100;
     nothrow @nogc pure @safe:
 
     pragma(LDC_intrinsic, "llvm.sqrt.f#")
@@ -141,26 +192,58 @@ version(LDC)
     ///
     T cos(T)(in T val) if (isFloatingPoint!T);
 
-    static if (LLVM_version >= 1300)
+    static if (llvmMajor >= 13)
         pragma(LDC_intrinsic, "llvm.powi.f#.i32")
-        ///
-        T powi(T)(in T val, int power) if (isFloatingPoint!T);
+        private T llvmPowi(T)(in T val, int power) if (isFloatingPoint!T);
     else 
         pragma(LDC_intrinsic, "llvm.powi.f#")
-        ///
-        T powi(T)(in T val, int power) if (isFloatingPoint!T);
+        private T llvmPowi(T)(in T val, int power) if (isFloatingPoint!T);
 
-    version(mir_core_test)
-    unittest
+    version (CRuntime_Microsoft)
+    // Keep LLVM from speculating both pow intrinsics before selecting a result.
+    pragma(inline, false)
+    private T widePowi(T)(T val, int power) if (isFloatingPoint!T)
     {
-        assert(powi(3.0, int(2)) == 9);
-        float f = 3;
-        assert(powi(f, int(2)) == 9);
+        return cast(T) llvmPowi(cast(double) val, power);
+    }
+
+    ///
+    T powi(T)(in T val, int power) if (isFloatingPoint!T)
+    {
+        // LLVM power intrinsics cannot be evaluated by CTFE.
+        if (__ctfe)
+        {
+            import std.math: pow;
+            return cast(T) pow(val, power);
+        }
+        version (CRuntime_Microsoft)
+        static if (T.mant_dig == float.mant_dig)
+        {
+            // LLVM lowers float powi to powf on this target, converting the
+            // exponent to float. Beyond 2^^24 this can lose integer parity.
+            // Double represents every int exactly, including int.min.
+            if (power < -0x1000000 || power > 0x1000000)
+                return widePowi(val, power);
+        }
+        return llvmPowi(val, power);
     }
 
     pragma(LDC_intrinsic, "llvm.pow.f#")
+    private T llvmPow(T)(in T val, in T power) if (isFloatingPoint!T);
+
     ///
-    T pow(T)(in T val, in T power) if (isFloatingPoint!T);
+    T pow(T)(in T val, in T power) if (isFloatingPoint!T)
+    {
+        if (__ctfe)
+        {
+            import std.math: pow;
+            Unqual!T result;
+            if (powSpecial!T(val, power, result))
+                return result;
+            return cast(T) pow(val, power);
+        }
+        return llvmPow(val, power);
+    }
 
     pragma(LDC_intrinsic, "llvm.exp.f#")
     ///
@@ -359,7 +442,14 @@ else static if (__VERSION__ >= 2082) // DMD 2.082 onward.
     ///
     T cos(T)(in T x) if (isFloatingPoint!T) { return std.math.cos(x); }
     ///
-    T pow(T)(in T x, in T power) if (isFloatingPoint!T) { alias y = power; mixin(mixinCMath2!`pow`); }
+    T pow(T)(in T x, in T power) if (isFloatingPoint!T)
+    {
+        Unqual!T result;
+        if (powSpecial!T(x, power, result))
+            return result;
+        alias y = power;
+        mixin(mixinCMath2!`pow`);
+    }
     ///
     T powi(T)(in T x, int power) if (isFloatingPoint!T) { alias y = power; mixin(mixinCMath2!`pow`); }
     ///
@@ -488,7 +578,14 @@ else // DMD version prior to 2.082
     ///
     T cos(T)(in T x) if (isFloatingPoint!T) { mixin(mixinCMath!`cos`); }
     ///
-    T pow(T)(in T x, in T power) if (isFloatingPoint!T) { alias y = power; mixin(mixinCMath2!`pow`); }
+    T pow(T)(in T x, in T power) if (isFloatingPoint!T)
+    {
+        Unqual!T result;
+        if (powSpecial!T(x, power, result))
+            return result;
+        alias y = power;
+        mixin(mixinCMath2!`pow`);
+    }
     ///
     T powi(T)(in T x, int power) if (isFloatingPoint!T) { alias y = power; mixin(mixinCMath2!`pow`); }
     ///
@@ -573,6 +670,156 @@ version (mir_core_test)
     assert(feqrel(pow(2.0L, -0.5L), cos(PI / 4)) >= real.mant_dig - 1);
 }
 
+version (mir_core_test)
+@nogc nothrow pure @safe unittest
+{
+    import std.meta: AliasSeq;
+
+    static foreach (T; AliasSeq!(float, double, real))
+    {{
+        // Exercise the CTFE path separately from the runtime implementation.
+        enum root = pow(T(4), T(0.5));
+        enum reciprocal = powi(T(2), -3);
+        static assert(is(typeof(root) == T));
+        static assert(is(typeof(reciprocal) == T));
+        static assert(root == 2 && reciprocal == T(0.125));
+        static assert(powi(T(-1), int.min) == 1);
+        static assert(powi(T(-1), int.max) == -1);
+        static assert(pow(T(0), T(0)) == 1);
+        static assert(powi(T(0), 0) == 1);
+
+        T base = 2;
+        assert(pow(base, T(3)) == 8);
+        assert(powi(base, -3) == reciprocal);
+        assert(pow(T(4), T(0.5)) == root);
+        assert(powi(T(-1), int.min) == 1);
+        assert(powi(T(-1), 31) == -1);
+
+        // Preserve explicit instantiation and function-address use.
+        auto power = &pow!T;
+        auto integerPower = &powi!T;
+        assert(power(base, T(3)) == 8);
+        assert(integerPower(base, -3) == reciprocal);
+    }}
+}
+
+version (mir_core_test)
+@nogc nothrow pure @safe unittest
+{
+    import std.meta: AliasSeq;
+    import std.math: signbit, isNaN;
+
+    // Keep the exponent unknown to the optimizer so these tests exercise
+    // runtime intrinsic lowering rather than constant folding.
+    pragma(inline, false)
+    static T runtimePowi(T)(T base, int exponent)
+    {
+        return powi(base, exponent);
+    }
+
+    static foreach (T; AliasSeq!(float, double, real))
+    {{
+        static assert(powi(T.nan, 0) == 1);
+        static assert(signbit(powi(-T(0), 3)));
+        static assert(powi(-T(0), -3) == -T.infinity);
+        static assert(signbit(powi(-T.infinity, -3)));
+
+        foreach (exponent; [int.min + 1, -0x1000001, -0xffffff,
+                           0xffffff, 0x1000001, int.max])
+        {
+            assert(runtimePowi(T(-1), exponent) == -1);
+            auto zeroPower = runtimePowi(-T(0), exponent);
+            auto infinityPower = runtimePowi(-T.infinity, exponent);
+            assert(signbit(zeroPower));
+            assert(signbit(infinityPower));
+            assert(zeroPower == (exponent < 0 ? -T.infinity : -T(0)));
+            assert(infinityPower == (exponent < 0 ? -T(0) : -T.infinity));
+        }
+        foreach (exponent; [int.min, -0x1000000, 0, 0x1000000, int.max - 1])
+            assert(runtimePowi(T(-1), exponent) == 1);
+        assert(runtimePowi(T.nan, 0) == 1);
+        assert(isNaN(runtimePowi(T.nan, 3)));
+        assert(runtimePowi(T(2), -3) == T(0.125));
+    }}
+}
+
+version (mir_core_test)
+@nogc nothrow pure @safe unittest
+{
+    import std.meta: AliasSeq;
+    import std.math: signbit, isNaN;
+
+    pragma(inline, false)
+    static T runtimePow(T)(T base, T exponent)
+    {
+        return pow(base, exponent);
+    }
+
+    static foreach (T; AliasSeq!(float, double, real))
+    {{
+        // Explicit expected values: do not use CTFE as the runtime oracle.
+        enum T[3][] specialCases = [
+            [T(1), T.nan, T(1)],
+            [T(1), T.infinity, T(1)],
+            [T(1), -T.infinity, T(1)],
+            [T(-1), T.infinity, T(1)],
+            [T(-1), -T.infinity, T(1)],
+            [T.nan, T(0), T(1)],
+            [T.nan, T(2), T.nan],
+            [T.nan, T.infinity, T.nan],
+            [T(2), T.nan, T.nan],
+            [-T.infinity, T.nan, T.nan],
+            [-T.infinity, -T.infinity, T(0)],
+            [-T.infinity, T.infinity, T.infinity],
+            [-T.infinity, T(-0.5), T(0)],
+            [-T.infinity, T(0.5), T.infinity],
+            [-T.infinity, T(-3), -T(0)],
+            [-T.infinity, T(3), -T.infinity],
+            [-T.infinity, T(-2), T(0)],
+            [-T.infinity, T(2), T.infinity],
+            [T.infinity, -T.infinity, T(0)],
+            [T.infinity, T.infinity, T.infinity],
+            [-T(0), T(-3), -T.infinity],
+            [-T(0), T(3), -T(0)],
+            [-T(0), T(-0.5), T.infinity],
+            [-T(0), T(0.5), T(0)],
+            [-T(0), T.infinity, T(0)],
+            [-T(0), -T.infinity, T.infinity],
+            [T(0), T(-3), T.infinity],
+            [T(0), T(3), T(0)],
+            [T(0.5), T.infinity, T(0)],
+            [T(-0.5), -T.infinity, T.infinity],
+            [T(-2), T.infinity, T.infinity],
+            [T(2), -T.infinity, T(0)],
+            [T(-2), T(0.5), T.nan]];
+
+        static foreach (entry; specialCases)
+        {{
+            enum expected = entry[2];
+            enum compileTime = pow(entry[0], entry[1]);
+            static if (isNaN(expected))
+                static assert(isNaN(compileTime));
+            else
+                static assert(compileTime == expected
+                    && signbit(compileTime) == signbit(expected));
+
+            auto actual = runtimePow(entry[0], entry[1]);
+            if (isNaN(expected))
+                assert(isNaN(actual));
+            else
+                assert(actual == expected && signbit(actual) == signbit(expected));
+        }}
+
+        static if (T.mant_dig >= 64)
+        {
+            enum T largeOdd = 0x1p63L + 1;
+            static assert(signbit(pow(-T(0), largeOdd)));
+            assert(signbit(runtimePow(-T(0), largeOdd)));
+            assert(runtimePow(-T.infinity, largeOdd) == -T.infinity);
+        }
+    }}
+}
+
 /// Overload for cdouble, cfloat and creal
 @optmath auto fabs(T)(in T x)
     if (isComplex!T)