Defer to libm cbrt for root(x,3).

This lets us get exact cases right only platforms where cbrt(x) has sub-ulp accuracy. It doesn't get all the other non-power-of-two roots quite right, and is still dependent on the host system, but 2 and 3 are the most common roots people take by far, so it's a good first step.

Sources/NumericsShims/include/NumericsShims.h

@@ -81,6 +81,10 @@ HEADER_SHIM float libm_powf(float x, float y) {
   return __builtin_powf(x, y);
 }
 
+HEADER_SHIM float libm_cbrtf(float x) {
+  return __builtin_cbrtf(x);
+}
+
 HEADER_SHIM float libm_atan2f(float y, float x) {
   return __builtin_atan2f(y, x);
 }
@@ -204,6 +208,10 @@ HEADER_SHIM double libm_pow(double x, double y) {
   return __builtin_pow(x, y);
 }
 
+HEADER_SHIM double libm_cbrt(double x) {
+  return __builtin_cbrt(x);
+}
+
 HEADER_SHIM double libm_atan2(double y, double x) {
   return __builtin_atan2(y, x);
 }
@@ -322,6 +330,10 @@ HEADER_SHIM long double libm_powl(long double x, long double y) {
   return __builtin_powl(x, y);
 }
 
+HEADER_SHIM long double libm_cbrtl(long double x) {
+  return __builtin_cbrtl(x);
+}
+
 HEADER_SHIM long double libm_atan2l(long double y, long double x) {
   return __builtin_atan2l(y, x);
 }

Sources/Real/Real.swift

@@ -39,15 +39,6 @@ extension Real {
     return pow(10, x)
   }
 
-  @_transparent
-  public static func root(_ x: Self, _ n: Int) -> Self {
-    guard x >= 0 || n % 2 != 0 else { return .nan }
-    // TODO: this implementation is not quite correct, because n may be
-    // rounded in conversion to Self. This only affects very extreme cases,
-    // so we'll leave it alone for now.
-    return Self(signOf: x, magnitudeOf: pow(x.magnitude, 1/Self(n)))
-  }
-
   #if !os(Windows)
   public static func signGamma(_ x: Self) -> FloatingPointSign {
     // Gamma is strictly positive for x >= 0.

Sources/Real/ScalarConformances.swift

@@ -57,6 +57,16 @@ extension Float: Real {
     return libm_powf(x, Float(n))
   }
 
+  @_transparent public static func root(_ x: Float, _ n: Int) -> Float {
+    guard x >= 0 || n % 2 != 0 else { return .nan }
+    // Workaround the issue mentioned below for the specific case of n = 3
+    // where we can fallback on cbrt.
+    if n == 3 { return libm_cbrtf(x) }
+    // TODO: this implementation is not quite correct, because either n or
+    // 1/n may be not be representable as Double.
+    return Float(signOf: x, magnitudeOf: libm_powf(x.magnitude, 1/Float(n)))
+  }
+
   @_transparent public static func atan2(y: Float, x: Float) -> Float {
     return libm_atan2f(y, x)
   }
@@ -126,6 +136,16 @@ extension Double: Real {
     return libm_pow(x, Double(n))
   }
 
+  @_transparent public static func root(_ x: Double, _ n: Int) -> Double {
+    guard x >= 0 || n % 2 != 0 else { return .nan }
+    // Workaround the issue mentioned below for the specific case of n = 3
+    // where we can fallback on cbrt.
+    if n == 3 { return libm_cbrt(x) }
+    // TODO: this implementation is not quite correct, because either n or
+    // 1/n may be not be representable as Double.
+    return Double(signOf: x, magnitudeOf: libm_pow(x.magnitude, 1/Double(n)))
+  }
+
   @_transparent public static func atan2(y: Double, x: Double) -> Double {
     return libm_atan2(y, x)
   }
@@ -178,6 +198,16 @@ extension Float80: Real {
     return libm_powl(x, Float80(n))
   }
 
+  @_transparent public static func root(_ x: Float80, _ n: Int) -> Float80 {
+    guard x >= 0 || n % 2 != 0 else { return .nan }
+    // Workaround the issue mentioned below for the specific case of n = 3
+    // where we can fallback on cbrt.
+    if n == 3 { return libm_cbrtl(x) }
+    // TODO: this implementation is not quite correct, because either n or
+    // 1/n may be not be representable as Double.
+    return Float80(signOf: x, magnitudeOf: libm_powl(x.magnitude, 1/Float80(n)))
+  }
+
   @_transparent public static func atan2(y: Float80, x: Float80) -> Float80 {
     return libm_atan2l(y, x)
   }

Tests/RealTests/RealTests.swift

@@ -61,6 +61,7 @@ internal extension ElementaryFunctions where Self: BinaryFloatingPoint {
     sanityCheck(-0.980829253011726236856451127452003999, Self.log(0.375))
     sanityCheck(0.3184537311185346158102472135905995955, Self.log(onePlus: 0.375))
     sanityCheck(-0.7211247851537041911608191553900547941, Self.root(-0.375, 3))
+    XCTAssertEqual(-10, Self.root(-1000, 3))
     sanityCheck(0.6123724356957945245493210186764728479, Self.sqrt(0.375))
     sanityCheck(0.54171335479545025876069682133938570, Self.pow(0.375, 0.625))
     sanityCheck(-0.052734375, Self.pow(-0.375, 3))