ElementaryFunctions compliance for Complex

Hi everybody,

I needed exponential and trigonometric functions for complex numbers, so I wrote code to add ElementaryFunctions compliance to Complex for Double, Float, and Float80. The code compiles and runs (on MacOS).

I've played around with it a bit, and I can get results, although the results for Float80 are wrong and I'm not sure how to fix that.

For example, I get these values for `cosh(Complex(8,0)

Float: (1490.4791, 0.0)
Double: (1490.479161252178, 0.0)
Float80: (2.7686593833135738328, 0.0)

I don't have have real implementations of expMinusOne or log(onePlus: ), as they aren't included in the CLang builtins. For these functions, I put in dummy code until I can implement a working solution for. For example:

  public static func expMinusOne(_ x: Complex<Double>) -> Complex<Double> {
    exp(x) - 1
  }

The code is available at GitHub - HaydenMcCabe/swift-numerics at ComplexElementaryFunctions

Here are some notes on what I've done so far:

_NumericsShims.h

I added additional inline functions to _NumericsShims.h in the style of the other functions, although because they need to use structs to pass function arguments, there is additional code in each function to cast types. Like the other functions, they only utilize functions that are included in the Clang builtins library. There are several functions for each type (real, imag, proj, conj, arg, abs), that are bridged to Swift, but not utilized.

The name of the struct CComplex was changed to CComplexDouble for consistency with the newly added structs CComplexFloat and CComplexFloat80. CComplex was not used outside of this file.

There are two existing functions, libm_cdiv and libm_cmul, which are not utilized, and can likely be removed, as ComplexModule implements this directly on the Complex type.

Complex+ElementaryFunctions.swift

This file defines the conditional conformance of Complex when its associated RealType conforms to the new protocol ComplexFunctionReal. The format is as follows:

extension Complex: ElementaryFunctions where RealType: ComplexFunctionReal {
    public static func exp(_ x: Complex<RealType>) -> Complex<RealType> {
        RealType.exp(x)
    }

    /// Additional function definitions
}

ComplexFunctionReal.swift

I don't care for the name, but I couldn't think of anything better to name it for now.

This protocol defines the functions of ElementaryFunctions with complex arguments and return types that conforming types must provide. It is defined as follows:

public protocol ComplexFunctionReal: Real {
    static func exp(_ x: Complex<Self>) -> Complex<Self>
    static func expMinusOne(_ x: Complex<Self>) -> Complex<Self>
    static func cosh(_ x: Complex<Self>) -> Complex<Self>
    static func sinh(_ x: Complex<Self>) -> Complex<Self>
    static func tanh(_ x: Complex<Self>) -> Complex<Self>
    static func cos(_ x: Complex<Self>) -> Complex<Self>
    static func sin(_ x: Complex<Self>) -> Complex<Self>
    static func tan(_ x: Complex<Self>) -> Complex<Self>
    static func log(_ x: Complex<Self>) -> Complex<Self>
    static func log(onePlus x: Complex<Self>) -> Complex<Self>
    static func acosh(_ x: Complex<Self>) -> Complex<Self>
    static func asinh(_ x: Complex<Self>) -> Complex<Self>
    static func atanh(_ x: Complex<Self>) -> Complex<Self>
    static func acos(_ x: Complex<Self>) -> Complex<Self>
    static func asin(_ x: Complex<Self>) -> Complex<Self>
    static func atan(_ x: Complex<Self>) -> Complex<Self>
    static func pow(_ x: Complex<Self>, _ y: Complex<Self>) -> Complex<Self>
    static func pow(_ x: Complex<Self>, _ n: Int) -> Complex<Self>
    static func sqrt(_ x: Complex<Self>) -> Complex<Self>
    static func root(_ x: Complex<Self>, _ n: Int) -> Complex<Self>   
}

Double+ComplexFunctionReal.swift
Float+ComplexFunctionReal.swift
Float80+ComplexFunctionReal.swift

These functions bridge the functions required by ComplexFunctionsReal, in the same style as the corresponding files in RealModule that provide compliance to Real. They simply cast between the Complex Swift struct to the corresponding C struct and call the functions defined in the _NumericsShims module as in the following code:

extension Double: ComplexFunctionReal {
  @_transparent
  public static func exp(_ x: Complex<Double>) -> Complex<Double> {
    let a = CComplexDouble(real: x.real, imag: x.imaginary)
    let b = libm_cexp(a)
    return Complex<Double>(b.real, b.imag)
  }
 
  //  Additional functions
}

Package.swift

The _NumericsShims target was added as a dependency of ComplexModule

DifferentiableTests.swift

I added a file to the ComplexTests test suite, although I haven't implemented the tests. There is one test defined, but it's only there to verify that tests could execute and the exp function can be called on Complex types. I'm not sure what would be required of this test suite, as I'm not versed in the issues that arise from complex floating-point types.

Please let me know if you have suggestions regarding the issues I mentioned. I'd like to see this functionality in the Numerics package if I can get these worked out.

Hi @HaydenMcCabe! I've implemented these functions in native Swift (and it can be done generically) here. These match the C/C++ choice of branch cuts too. They'd be trivial to adapt for Swift Numerics (just be aware of some of the differences in design choices with regard to magnitude, etc.) and don't require a shim.

2 Likes

I might’ve missed it, but I don’t see equivalents of expm1 and log1p for Complex in your Numeric Annex.

I’m not looking at a numerical recipes reference, but just with pen and paper I’ve derived the following possible implementations:

extension Complex {
  static func expMinusOne(_ z: Self) -> Self {
    let (x, y) = (z.x, z.y)
    let a = RealType.expMinusOne(x)
    let c = RealType.cos(y)
    let s = RealType.sin(y)
    let u = a*c - s*s / (1 + c)
    let v = a*s + s
    return Self(u, v)
  }
  
  static func logOnePlus(_ z: Self) -> Self {
    let (x, y) = (z.x, z.y)
    let a = x*x + y*y + 2*x
    let u = RealType.logOnePlus(a) / 2
    let v = RealType.atan2(y: y, x: x+1)
    return Self(u, v)
  }
}

They’ll need auditing to properly handle non-finite values, as well as potential overflow, but the core calculations should be correct and maintain proper precision.

1 Like

Wow. That’s awesome. Thanks, I’ll take a look at that.

Thanks. I’ll try looking at how this works.