Merge pull request #83 from stephentyrone/pown

Initial pass at pow(Self, Int) for Float and Double

Author: stephentyrone

Sources/Real/ScalarConformances.swift

@@ -50,11 +50,21 @@ extension Float: Real {
   }
 
   @_transparent public static func pow(_ x: Float, _ n: Int) -> Float {
-    // TODO: this implementation is not quite correct, because n may be
-    // rounded in conversion to Float. This only effects very extreme cases,
-    // so we'll leave it alone for now; however, it gets the sign wrong if
-    // it rounds an odd number to an even number, so we should fix it soon.
-    return libm_powf(x, Float(n))
+    // If n is exactly representable as Float, we can just call powf:
+    if let y = Float(exactly: n) {
+      return libm_powf(x, y)
+    }
+    // Otherwise, n is too large to losslessly represent as Float.
+    // The range of "interesting" n is -1488522191 ... 1744361944; outside
+    // of this range, all x != 1 overflow or underflow, so only the parity
+    // of x matters. We don't really care about the specific range at all,
+    // only that the bounds fit exactly into two Floats. Mask the low 24
+    // bits of n, get pow with that exponent (this contains the parity),
+    // then get pow with the rest (this may round, but if it does we've
+    // saturated anyway, so it doesn't matter).
+    let low = n & 0xffffff
+    let high = n - low
+    return libm_powf(x, Float(low)) * libm_powf(x, Float(high))
   }
 
   @_transparent public static func root(_ x: Float, _ n: Int) -> Float {
@@ -129,11 +139,22 @@ extension Double: Real {
   }
 
   @_transparent public static func pow(_ x: Double, _ n: Int) -> Double {
-    // TODO: this implementation is not quite correct, because n may be
-    // rounded in conversion to Double. This only effects very extreme cases,
-    // so we'll leave it alone for now; however, it gets the sign wrong if
-    // it rounds an odd number to an even number, so we should fix it soon.
-    return libm_pow(x, Double(n))
+    // If n is exactly representable as Double, we can just call pow:
+    // Note that all calls on a 32b platform go down this path.
+    if let y = Double(exactly: n) {
+      return libm_pow(x, y)
+    }
+    // Otherwise, n is too large to losslessly represent as Double, so we
+    // just split it into two parts, high and low. This is always exact,
+    // so the only source of error is pow itself and the multiplication.
+    //
+    // mask constant is spelled in this funny way because if we just anded
+    // with the hex value, we'd get a compile error on 32b platforms, even
+    // though this whole branch is dead code on 32b.
+    let mask = Int(truncatingIfNeeded: 0x1fffffffffffff as UInt64)
+    let low = n & mask
+    let high = n - low
+    return libm_pow(x, Double(low)) * libm_pow(x, Double(high))
   }
 
   @_transparent public static func root(_ x: Double, _ n: Int) -> Double {
@@ -191,10 +212,9 @@ extension Float80: Real {
   }
 
   @_transparent public static func pow(_ x: Float80, _ n: Int) -> Float80 {
-    // TODO: this implementation is not quite correct, because n may be
-    // rounded in conversion to Float80. This only effects very extreme cases,
-    // so we'll leave it alone for now; however, it gets the sign wrong if
-    // it rounds an odd number to an even number, so we should fix it soon.
+    // Every Int value is exactly representable as Float80, so we don't need
+    // to do anything fancy--unlike Float and Double, we can just call the
+    // libm pow function.
     return libm_powl(x, Float80(n))
   }
 

Tests/RealTests/RealTests.swift

@@ -1,4 +1,4 @@
-//===--- ElementaryFunctionTests.swift ------------------------*- swift -*-===//
+//===--- ElementaryFunctionChecks.swift ------------------------*- swift -*-===//
 //
 // This source file is part of the Swift Numerics open source project
 //
@@ -22,28 +22,42 @@ func sanityCheck<T>(_ expected: TestLiteralType, _ actual: T,
                     ulps allowed: T = 16,
                     file: StaticString = #file, line: UInt = #line)
   where T: BinaryFloatingPoint {
+  // Shortcut relative-error check if we got the sign wrong; it's OK to
+  // underflow to zero a little bit early, but we don't want to allow going
+  // right through zero to the other side.
+  XCTAssert(actual.sign == expected.sign, "\(actual) != \(expected) as \(T.self)", file: file, line: line)
   // Default tolerance is 16 ulps; It's OK to relax this as needed for new
-  // platforms, as these tests are *not* intended to validate the math
+  // platforms, as these checks are *not* intended to validate the math
   // library--they are only intended to check that the Swift bindings are
   // calling the right functions in the math library. It's important, however
-  // not to relax the tolerance beyond a few hundred ulps, because these tests
+  // not to relax the tolerance beyond a few hundred ulps, because these checks
   // need to detect errors where the *wrong function* is being called; e.g.
   // we need to flag an implentation that inadvertently called the C hypotf
   // function instead of hypot. This is especially important because the C
   // shims that we're calling through will allow silent type conversions.
   if actual == T(expected) || actual.isNaN && expected.isNaN {
     return
   }
-  //  Compute error in ulp, compare to tolerance.
-  let absoluteError = T(abs(TestLiteralType(actual) - expected))
-  let ulpError = absoluteError / T(expected).ulp
-  XCTAssert(ulpError <= allowed, "\(actual) != \(expected) as \(T.self)" +
-            "\n  \(ulpError)-ulp error exceeds \(allowed)-ulp tolerance.",
-            file: file, line: line)
+  // Special-case where expected or observed is infinity.
+  // Artificially knock everything down a binade, treat actual infinity as
+  // the base of the next binade up.
+  if actual.isInfinite || T(expected).isInfinite {
+    let scaledExpected = TestLiteralType(signOf: expected,
+      magnitudeOf: expected.isInfinite ? TestLiteralType.greatestFiniteMagnitude.binade : 0.5 * expected
+    )
+    let scaledActual = T(signOf: actual,
+      magnitudeOf: actual.isInfinite ? T.greatestFiniteMagnitude.binade : 0.5 * actual
+    )
+    return sanityCheck(scaledExpected, scaledActual, ulps: allowed, file: file, line: line)
+  }
+  // Compute error in ulp, compare to tolerance.
+  let absoluteError = T(abs(TestLiteralType(actual) - expected)).magnitude
+  let ulpError = absoluteError / max(T(expected).magnitude, T.leastNormalMagnitude).ulp
+  XCTAssert(ulpError <= allowed, "\(actual) != \(expected) as \(T.self)\n\(ulpError)-ulp error exceeds \(allowed)-ulp tolerance.", file: file, line: line)
 }
 
 internal extension ElementaryFunctions where Self: BinaryFloatingPoint {
-  static func elementaryFunctionTests() {
+  static func elementaryFunctionChecks() {
     sanityCheck(1.1863995522992575361931268186727044683, Self.acos(0.375))
     sanityCheck(0.3843967744956390830381948729670469737, Self.asin(0.375))
     sanityCheck(0.3587706702705722203959200639264604997, Self.atan(0.375))
@@ -69,7 +83,7 @@ internal extension ElementaryFunctions where Self: BinaryFloatingPoint {
 }
 
 internal extension Real where Self: BinaryFloatingPoint {
-  static func realFunctionTests() {
+  static func realFunctionChecks() {
     sanityCheck(1.2968395546510096659337541177924511598, Self.exp2(0.375))
     sanityCheck(2.3713737056616552616517527574788898386, Self.exp10(0.375))
     sanityCheck(-1.415037499278843818546261056052183491, Self.log2(0.375))
@@ -85,25 +99,121 @@ internal extension Real where Self: BinaryFloatingPoint {
     XCTAssertEqual(.minus, Self.signGamma(-2.375))
     #endif
   }
+  
+  static func testPownCommon() {
+    // If x is -1, then the result is ±1 with sign chosen by parity of n.
+    // Simply converting n to Real will flip parity when n is large, so
+    // first check that we get those cases right.
+    XCTAssertEqual(Self.pow(-1,  0),  1)
+    XCTAssertEqual(Self.pow(-1,  1), -1)
+    XCTAssertEqual(Self.pow(-1, -1), -1)
+    XCTAssertEqual(Self.pow(-1,  2),  1)
+    XCTAssertEqual(Self.pow(-1, -2),  1)
+    XCTAssertEqual(Self.pow(-1,  Int.max - 1), 1)
+    XCTAssertEqual(Self.pow(-1, -Int.max + 1), 1)
+    XCTAssertEqual(Self.pow(-1,  Int.max), -1)
+    XCTAssertEqual(Self.pow(-1, -Int.max), -1)
+    XCTAssertEqual(Self.pow(-1,  Int.min),  1)
+  }
+}
+
+extension Float {
+  static func testPown() {
+    testPownCommon()
+    let u = Float(1).nextUp
+    let d = Float(1).nextDown
+    // Smallest exponents not exactly representable as Float.
+    sanityCheck(-7.3890560989306677280287919329569359, Float.pow(-u, 0x1000001))
+    sanityCheck(-0.3678794082804575860056608283059288, Float.pow(-d, 0x1000001))
+    // Exponents close to overflow boundary.
+    sanityCheck(-3.4028231352500001570898203463449749e38, Float.pow(-u, 744261161))
+    sanityCheck( 3.4028235408981285772043562848249166e38, Float.pow(-u, 744261162))
+    sanityCheck(-3.4028239465463053543440887892352174e38, Float.pow(-u, 744261163))
+    sanityCheck( 3.4028233551634475284795244782720072e38, Float.pow(-d, -1488522190))
+    sanityCheck(-3.4028235579875369356575053576685267e38, Float.pow(-d, -1488522191))
+    sanityCheck( 3.4028237608116384320940078199368685e38, Float.pow(-d, -1488522192))
+    // Exponents close to underflow boundary.
+    sanityCheck( 7.0064936491761438872280296737844625e-46, Float.pow(-u, -872181048))
+    sanityCheck(-7.0064928139371132951305928725186420e-46, Float.pow(-u, -872181049))
+    sanityCheck( 7.0064919786981822712727285793333389e-46, Float.pow(-u, -872181050))
+    sanityCheck(-7.0064924138100205091278464932003585e-46, Float.pow(-d, 1744361943))
+    sanityCheck( 7.0064919961905290625123586120258840e-46, Float.pow(-d, 1744361944))
+    sanityCheck(-7.0064915785710625079583096856510544e-46, Float.pow(-d, 1744361945))
+    // Just hammer max/min exponents, these always saturate, but this will reveal
+    // errors in some implementations that one could try.
+    sanityCheck( .infinity, Self.pow(-u,  Int.max - 1))
+    sanityCheck( 0.0,       Self.pow(-d,  Int.max - 1))
+    sanityCheck( 0.0,       Self.pow(-u, -Int.max + 1))
+    sanityCheck( .infinity, Self.pow(-d, -Int.max + 1))
+    sanityCheck(-.infinity, Self.pow(-u,  Int.max))
+    sanityCheck(-0.0,       Self.pow(-d,  Int.max))
+    sanityCheck(-0.0,       Self.pow(-u, -Int.max))
+    sanityCheck(-.infinity, Self.pow(-d, -Int.max))
+    sanityCheck( 0.0,       Self.pow(-u,  Int.min))
+    sanityCheck( .infinity, Self.pow(-d,  Int.min))
+  }
 }
 
-final class ElementaryFunctionTests: XCTestCase {
+extension Double {
+  static func testPown() {
+    testPownCommon()
+    // Following tests only make sense (and are only necessary) on 64b platforms.
+#if arch(arm64) || arch(x86_64)
+    let u: Double = 1.nextUp
+    let d: Double = 1.nextDown
+    // Smallest exponent not exactly representable as Double.
+    sanityCheck(-7.3890560989306502272304274605750685, Double.pow(-u, 0x20000000000001))
+    sanityCheck(-0.1353352832366126918939994949724833, Double.pow(-u, -0x20000000000001))
+    sanityCheck(-0.3678794411714422603312898889458068, Double.pow(-d, 0x20000000000001))
+    sanityCheck(-2.7182818284590456880451484776630468, Double.pow(-d, -0x20000000000001))
+    // Exponents close to overflow boundary.
+    sanityCheck( 1.7976931348623151738531864721534215e308, Double.pow(-u, 3196577161300664268))
+    sanityCheck(-1.7976931348623155730212483790972209e308, Double.pow(-u, 3196577161300664269))
+    sanityCheck( 1.7976931348623159721893102860411089e308, Double.pow(-u, 3196577161300664270))
+    sanityCheck( 1.7976931348623157075547244136070910e308, Double.pow(-d, -6393154322601327474))
+    sanityCheck(-1.7976931348623159071387553670790721e308, Double.pow(-d, -6393154322601327475))
+    sanityCheck( 1.7976931348623161067227863205510754e308, Double.pow(-d, -6393154322601327476))
+    // Exponents close to underflow boundary.
+    sanityCheck( 2.4703282292062334560337346683707907e-324, Double.pow(-u, -3355781687888880946))
+    sanityCheck(-2.4703282292062329075106789791206172e-324, Double.pow(-u, -3355781687888880947))
+    sanityCheck( 2.4703282292062323589876232898705654e-324, Double.pow(-u, -3355781687888880948))
+    sanityCheck(-2.4703282292062332640976590913373022e-324, Double.pow(-d, 6711563375777760775))
+    sanityCheck( 2.4703282292062329898361312467121758e-324, Double.pow(-d, 6711563375777760776))
+    sanityCheck(-2.4703282292062327155746034020870799e-324, Double.pow(-d, 6711563375777760777))
+    // Just hammer max/min exponents, these always saturate, but this will reveal
+    // errors in some implementations that one could try.
+    sanityCheck( .infinity, Self.pow(-u,  Int.max - 1))
+    sanityCheck( 0.0,       Self.pow(-d,  Int.max - 1))
+    sanityCheck( 0.0,       Self.pow(-u, -Int.max + 1))
+    sanityCheck( .infinity, Self.pow(-d, -Int.max + 1))
+    sanityCheck(-.infinity, Self.pow(-u,  Int.max))
+    sanityCheck(-0.0,       Self.pow(-d,  Int.max))
+    sanityCheck(-0.0,       Self.pow(-u, -Int.max))
+    sanityCheck(-.infinity, Self.pow(-d, -Int.max))
+    sanityCheck( 0.0,       Self.pow(-u,  Int.min))
+    sanityCheck( .infinity, Self.pow(-d,  Int.min))
+#endif
+  }
+}
+
+final class ElementaryFunctionChecks: XCTestCase {
 
   func testFloat() {
-    Float.elementaryFunctionTests()
-    Float.realFunctionTests()
+    Float.elementaryFunctionChecks()
+    Float.realFunctionChecks()
+    Float.testPown()
   }
 
   func testDouble() {
-    Double.elementaryFunctionTests()
-    Double.realFunctionTests()
+    Double.elementaryFunctionChecks()
+    Double.realFunctionChecks()
+    Double.testPown()
   }
 
-  
   #if (arch(i386) || arch(x86_64)) && !os(Windows) && !os(Android)
   func testFloat80() {
-    Float80.elementaryFunctionTests()
-    Float80.realFunctionTests()
+    Float80.elementaryFunctionChecks()
+    Float80.realFunctionChecks()
   }
   #endif
 }