Merge pull request #145 from stephentyrone/mulAdd-and-cosm1

Added mulAdd and cosMinusOne operations for Real.

Author: stephentyrone

Sources/RealModule/Double+Real.swift

@@ -218,4 +218,9 @@ extension Double: Real {
     return libm_lgamma(x, &dontCare)
   }
   #endif
+  
+  @_transparent
+  public static func _mulAdd(_ a: Double, _ b: Double, _ c: Double) -> Double {
+    _numerics_muladd(a, b, c)
+  }
 }

Sources/RealModule/Float+Real.swift

@@ -191,4 +191,9 @@ extension Float: Real {
     return libm_lgammaf(x, &dontCare)
   }
   #endif
+  
+  @_transparent
+  public static func _mulAdd(_ a: Float, _ b: Float, _ c: Float) -> Float {
+    _numerics_muladdf(a, b, c)
+  }
 }

Sources/RealModule/Float16+Real.swift

@@ -172,5 +172,10 @@ extension Float16: Real {
     Float16(.logGamma(Float(x)))
   }
   #endif
+  
+  @_transparent
+  public static func _mulAdd(_ a: Float16, _ b: Float16, _ c: Float16) -> Float16 {
+    _numerics_muladdf16(a, b, c)
+  }
 }
 #endif

Sources/RealModule/Real.swift

@@ -40,6 +40,17 @@ extension Real {
     return pow(10, x)
   }
 
+  /// cos(x) - 1, computed in such a way as to maintain accuracy for small x.
+  ///
+  /// See also:
+  /// -
+  /// - `ElementaryFunctions.expMinusOne()`
+  @_transparent
+  public static func cosMinusOne(_ x: Self) -> Self {
+    let sinxOver2 = sin(x/2)
+    return -2*sinxOver2*sinxOver2
+  }
+  
   #if !os(Windows)
   public static func signGamma(_ x: Self) -> FloatingPointSign {
     // Gamma is strictly positive for x >= 0.
@@ -76,6 +87,11 @@ extension Real {
   }
   #endif
 
+  @_transparent
+  public static func _mulAdd(_ a: Self, _ b: Self, _ c: Self) -> Self {
+    a*b + c
+  }
+  
   @_transparent
   public static func sqrt(_ x: Self) -> Self {
     return x.squareRoot()

Sources/RealModule/RealFunctions.swift

@@ -112,4 +112,9 @@ public protocol RealFunctions: ElementaryFunctions {
   /// - `logGamma()`
   static func signGamma(_ x: Self) -> FloatingPointSign
 #endif
+  
+  /// a*b + c, computed _either_ with an FMA or with separate multiply and add.
+  ///
+  /// Whichever is faster should be chosen by the compiler statically.
+  static func _mulAdd(_ a: Self, _ b: Self, _ c: Self) -> Self
 }

Sources/_NumericsShims/include/_NumericsShims.h

@@ -382,6 +382,28 @@ HEADER_SHIM long double libm_lgammal(long double x, int *signp) {
 }
 #endif
 
+// MARK: - fast mul-add inlines
+/// a*b + c evaluated _either_ as two operations or fma, whichever is faster.
+HEADER_SHIM _Float16 _numerics_muladdf16(_Float16 a, _Float16 b, _Float16 c) {
+#pragma STDC FP_CONTRACT ON
+  return a*b + c;
+}
+
+/// a*b + c evaluated _either_ as two operations or fma, whichever is faster.
+HEADER_SHIM float _numerics_muladdf(float a, float b, float c) {
+#pragma STDC FP_CONTRACT ON
+  return a*b + c;
+}
+
+/// a*b + c evaluated _either_ as two operations or fma, whichever is faster.
+HEADER_SHIM double _numerics_muladd(double a, double b, double c) {
+#pragma STDC FP_CONTRACT ON
+  return a*b + c;
+}
+
+// No long-double muladd operation, because no one has built an FMA for it
+// (except for Itanium, which Swift doesn't support).
+
 // MARK: - shims to import C complex operations for timing purposes
 // Clang doesn't provide complex arithmetic on Windows (because MSVC
 // doesn't), so we can't define these there, or we'll get link errors.