Extend all analytic functions in base for quaternions

[sethaxen]

All complex analytic functions f in Base can be extended to the quaternions with the following function:

using Quaternions

function quat_analytic(f, q::Quaternion)
    a = Quaternions.abs_imag(q)
    z = complex(q.s, a)
    fz = f(z)
    wr, wi = reim(fz)
    scale = wi / a
    if a > 0
        return Quaternion(wr, scale * q.v1, scale * q.v2, scale * q.v3)
    else  # quaternion may be real or complex
        return Quaternion(wr, oftype(scale, wi), zero(scale), zero(scale))
    end
end

The function is adapted from Theorem 5 of doi.org/10.1017/S0305004100055638. It can also be derived a completely different way by mapping a quaternion to a 2x2 skew-Hermitian matrix whose product preserves the Hamilton product, extending the analytic function to a matrix function, and then mapping back to the equivalent quaternions.

This function gives equivalent results to the current functions implementations here, but it's faster than all but log:

julia> using BenchmarkTools

julia> q = quatrand()
Quaternion{Float64}(-1.76527498999994, -0.9985320755133996, -0.606774830599748, 1.0149711723420227, false)

julia> @btime exp($q)
  31.384 ns (0 allocations: 0 bytes)
Quaternion{Float64}(0.003950512836446494, -0.11038429890473422, -0.0670768780605819, 0.11220158472113788, false)

julia> @btime $quat_analytic(exp, $q)
  28.991 ns (0 allocations: 0 bytes)
Quaternion{Float64}(0.003950512836446494, -0.11038429890473421, -0.06707687806058188, 0.11220158472113786, false)

julia> @btime log($q)
  74.247 ns (0 allocations: 0 bytes)
Quaternion{Float64}(0.8534278686043146, -1.5624466956559513, -0.9494470456448084, 1.5881696675557215, false)

julia> @btime $quat_analytic(log, $q)
  81.745 ns (0 allocations: 0 bytes)
Quaternion{Float64}(0.8534278686043146, -1.5624466956559513, -0.9494470456448083, 1.5881696675557215, false)

julia> @btime sqrt($q)
  131.847 ns (0 allocations: 0 bytes)
Quaternion{Float64}(0.539632108658829, -0.9251970550780371, -0.5622115704973448, 0.9404288181299799, false)

julia> @btime $quat_analytic(sqrt, $q)
  30.934 ns (0 allocations: 0 bytes)
Quaternion{Float64}(0.5396321086588289, -0.9251970550780371, -0.5622115704973447, 0.9404288181299798, false)

julia> @btime sin($q)
  158.434 ns (0 allocations: 0 bytes)
Quaternion{Float64}(-2.410402900676647, 0.2797836189575866, 0.1700152275131357, -0.28438977044324726, false)

julia> @btime $quat_analytic(sin, $q)
  48.684 ns (0 allocations: 0 bytes)
Quaternion{Float64}(-2.410402900676648, 0.2797836189575863, 0.17001522751313533, -0.2843897704432474, false)

julia> quat_analytic(tan, q) ≈ sin(q) / cos(q) ≈ cos(q) \ sin(q)
true

julia> quat_analytic(sin, quat_analytic(asin, q)) ≈ q
true

julia> quat_analytic(cosh, quat_analytic(acosh, q)) ≈ q
true

julia> quat_analytic(exp10, q) ≈ 10^q
true

julia> quat_analytic(inv, q) ≈ inv(q)
true