Merge pull request #235 from stephentyrone/complex-pow

Fix edge-cases for Complex.pow(.zero, <Int>), and Real.pow(zero, zero)

Author: stephentyrone

Sources/ComplexModule/Complex+ElementaryFunctions.swift

@@ -420,7 +420,9 @@ extension Complex: ElementaryFunctions {
 
   @inlinable
   public static func pow(_ z: Complex, _ n: Int) -> Complex {
-    if z.isZero { return .zero }
+    if z.isZero {
+      return n < 0 ? .infinity : n == 0 ? .one : .zero
+    }
     // TODO: this implementation is not quite correct, because n may be
     // rounded in conversion to RealType. This only effects very extreme
     // cases, so we'll leave it alone for now.

Sources/RealModule/Double+Real.swift

@@ -107,14 +107,14 @@ extension Double: Real {
     libm_exp2(x)
   }
 
-  #if os(macOS) || os(iOS) || os(tvOS) || os(watchOS)
+#if os(macOS) || os(iOS) || os(tvOS) || os(watchOS)
   @_transparent
   public static func exp10(_ x: Double) -> Double {
     libm_exp10(x)
   }
-  #endif
+#endif
 
-  #if os(macOS) && arch(x86_64)
+#if os(macOS) && arch(x86_64)
   // Workaround for macOS bug (<rdar://problem/56844150>) where hypot can
   // overflow for values very close to the overflow boundary of the naive
   // algorithm. Since this is only for macOS, we can just unconditionally
@@ -125,12 +125,12 @@ extension Double: Real {
     let y80 = Float80(y)
     return Double(Float80.sqrt(x80*x80 + y80*y80))
   }
-  #else
+#else
   @_transparent
   public static func hypot(_ x: Double, _ y: Double) -> Double {
     libm_hypot(x, y)
   }
-  #endif
+#endif
 
   @_transparent
   public static func gamma(_ x: Double) -> Double {
@@ -150,6 +150,7 @@ extension Double: Real {
   @_transparent
   public static func pow(_ x: Double, _ y: Double) -> Double {
     guard x >= 0 else { return .nan }
+    if x == 0 && y == 0 { return .nan }
     return libm_pow(x, y)
   }
 
@@ -211,13 +212,13 @@ extension Double: Real {
     libm_atan2(y, x)
   }
 
-  #if !os(Windows)
+#if !os(Windows)
   @_transparent
   public static func logGamma(_ x: Double) -> Double {
     var dontCare: Int32 = 0
     return libm_lgamma(x, &dontCare)
   }
-  #endif
+#endif
 
   @_transparent
   public static func _mulAdd(_ a: Double, _ b: Double, _ c: Double) -> Double {

Sources/RealModule/ElementaryFunctions.swift

@@ -228,14 +228,24 @@ public protocol ElementaryFunctions: AdditiveArithmetic {
 
   /// exp(y * log(x)) computed with additional internal precision.
   ///
-  /// See also `sqrt()` and `root()`.
+  /// The edge-cases of this function are defined based on the behavior of the
+  /// expression `exp(y log x)`, matching IEEE 754's `powr` operation.
+  /// In particular, this means that if `x` and `y` are both zero, `pow(x,y)`
+  /// is `nan` for real types and `infinity` for complex types, rather than 1.
+  ///
+  /// There is also a `pow(_:Self,_:Int)` overload, whose behavior is defined
+  /// in terms of repeated multiplication, and hence returns 1 for this case.
   ///
+  /// See also `sqrt()` and `root()`.
   static func pow(_ x: Self, _ y: Self) -> Self
 
   /// `x` raised to the nth power.
   ///
-  /// See also `sqrt()` and `root()`.
+  /// The edge-cases of this function are defined in terms of repeated
+  /// multiplication or division, rather than exp(n log x). In particular,
+  /// `Float.pow(0, 0)` is 1.
   ///
+  /// See also `sqrt()` and `root()`.
   static func pow(_ x: Self, _ n: Int) -> Self
 
   /// The [square root][wiki] of `x`.
@@ -248,6 +258,5 @@ public protocol ElementaryFunctions: AdditiveArithmetic {
   /// The nth root of `x`.
   ///
   /// See also `pow()` and `sqrt()`.
-  ///
   static func root(_ x: Self, _ n: Int) -> Self
 }

Sources/RealModule/Float+Real.swift

@@ -137,6 +137,7 @@ extension Float: Real {
   @_transparent
   public static func pow(_ x: Float, _ y: Float) -> Float {
     guard x >= 0 else { return .nan }
+    if x == 0 && y == 0 { return .nan }
     return libm_powf(x, y)
   }
 

Sources/RealModule/Float80+Real.swift

@@ -132,6 +132,7 @@ extension Float80: Real {
   @_transparent
   public static func pow(_ x: Float80, _ y: Float80) -> Float80 {
     guard x >= 0 else { return .nan }
+    if x == 0 && y == 0 { return .nan }
     return libm_powl(x, y)
   }
 

Tests/ComplexTests/ElementaryFunctionTests.swift

@@ -400,6 +400,18 @@ final class ElementaryFunctionTests: XCTestCase {
     }
   }
 
+  func testPowR<T: Real & FixedWidthFloatingPoint>(_ type: T.Type) {
+    XCTAssertEqual(Complex<T>.pow(.zero, -.one),  .infinity)
+    XCTAssertEqual(Complex<T>.pow(.zero,  .zero), .infinity)
+    XCTAssertEqual(Complex<T>.pow(.zero, +.one),  .infinity)
+  }
+  
+  func testPowN<T: Real & FixedWidthFloatingPoint>(_ type: T.Type) {
+    XCTAssertEqual(Complex<T>.pow(.zero, -1), .infinity)
+    XCTAssertEqual(Complex<T>.pow(.zero,  0), .one)
+    XCTAssertEqual(Complex<T>.pow(.zero, +1), .zero)
+  }
+  
   func testFloat() {
     testExp(Float.self)
     testExpMinusOne(Float.self)
@@ -411,6 +423,8 @@ final class ElementaryFunctionTests: XCTestCase {
     testAcosh(Float.self)
     testAsinh(Float.self)
     testAtanh(Float.self)
+    testPowR(Float.self)
+    testPowN(Float.self)
   }
 
   func testDouble() {
@@ -424,9 +438,11 @@ final class ElementaryFunctionTests: XCTestCase {
     testAcosh(Double.self)
     testAsinh(Double.self)
     testAtanh(Double.self)
+    testPowR(Double.self)
+    testPowN(Double.self)
   }
 
-  #if (arch(i386) || arch(x86_64)) && !os(Windows) && !os(Android)
+#if (arch(i386) || arch(x86_64)) && !os(Windows) && !os(Android)
   func testFloat80() {
     testExp(Float80.self)
     testExpMinusOne(Float80.self)
@@ -438,6 +454,8 @@ final class ElementaryFunctionTests: XCTestCase {
     testAcosh(Float80.self)
     testAsinh(Float80.self)
     testAtanh(Float80.self)
+    testPowR(Float80.self)
+    testPowN(Float80.self)
   }
-  #endif
+#endif
 }

Tests/RealTests/ElementaryFunctionChecks.swift

@@ -120,39 +120,59 @@ internal extension Real where Self: BinaryFloatingPoint {
     assertClose(0.4041169094348222983238250859191217675, Self.erf(0.375))
     assertClose(0.5958830905651777016761749140808782324, Self.erfc(0.375))
     assertClose(2.3704361844166009086464735041766525098, Self.gamma(0.375))
-    #if !os(Windows)
+#if !os(Windows)
     assertClose( -0.11775527074107877445136203331798850, Self.logGamma(1.375))
     XCTAssertEqual(.plus,  Self.signGamma(1.375))
     XCTAssertEqual(.minus, Self.signGamma(-2.375))
-    #endif
+#endif
+  }
+}
+
+extension Real {
+  static func powZeroChecks() {
+    // pow(_:Self,_:Self) is defined by exp(y log(x)) and has edge-cases to
+    // match. In particular, if x is zero, log(x) is -infinity, so pow(0,0)
+    // is exp(0 * -infinity) = exp(nan) = nan.
+    XCTAssertEqual(pow(0, -1 as Self), infinity)
+    XCTAssert(pow(0, 0 as Self).isNaN)
+    XCTAssertEqual(pow(0,  1 as Self), zero)
+    // pow(_:Self,_:Int) is defined by repeated multiplication or division,
+    // and hence pow(0, 0) is 1.
+    XCTAssertEqual(pow(0, -1), infinity)
+    XCTAssertEqual(pow(0,  0), 1)
+    XCTAssertEqual(pow(0,  1), zero)
   }
 }
 
 final class ElementaryFunctionChecks: XCTestCase {
 
-  #if !((os(macOS) || targetEnvironment(macCatalyst)) && arch(x86_64))
+#if !((os(macOS) || targetEnvironment(macCatalyst)) && arch(x86_64))
   func testFloat16() {
     if #available(macOS 11.0, iOS 14.0, watchOS 14.0, tvOS 7.0, *) {
       Float16.elementaryFunctionChecks()
       Float16.realFunctionChecks()
+      Float16.powZeroChecks()
     }
   }
-  #endif
+#endif
 
   func testFloat() {
     Float.elementaryFunctionChecks()
     Float.realFunctionChecks()
+    Float.powZeroChecks()
   }
 
   func testDouble() {
     Double.elementaryFunctionChecks()
     Double.realFunctionChecks()
+    Double.powZeroChecks()
   }
 
-  #if (arch(i386) || arch(x86_64)) && !os(Windows) && !os(Android)
+#if (arch(i386) || arch(x86_64)) && !os(Windows) && !os(Android)
   func testFloat80() {
     Float80.elementaryFunctionChecks()
     Float80.realFunctionChecks()
+    Float80.powZeroChecks()
   }
-  #endif
+#endif
 }