BiComplex numbers in Julia

BiComplexes.jl

module BiComplexes

using AutoHashEquals

export BiComplex, derivative

@auto_hash_equals struct BiComplex{T} <: Number
  re::Complex{T}
  im::Complex{T}
end
const imim  = BiComplex(0 + 0im, 1 + 0im)

Base.isequal(a::BiComplex, b::Complex) = (a.re == b.re) && (a.im == b.im)
Base.real(b::BiComplex) = b.re
Base.imag(b::BiComplex) = b.im

Base.conj(a::BiComplex) = BiComplex(a.re, -a.im)
Base.abs2(a::BiComplex) = abs2(a.re) + abs2(a.im)
Base.abs(a::BiComplex) = sqrt(abs2(a.re) + abs2(a.im))
inv(a::BiComplex) = conj(a) / (a.re^2 + a.im^2)

Base.:+(a::BiComplex, b::BiComplex) = BiComplex(a.re + b.re, a.im + b.im)
Base.:-(a::BiComplex, b::BiComplex) = BiComplex(a.re - b.re, a.im - b.im)
Base.:+(a::BiComplex, b::Number) = BiComplex(a.re + b, a.im)
Base.:-(a::BiComplex, b::Number) = BiComplex(a.re - b, a.im)
Base.:-(a::Number, b::BiComplex) = BiComplex(a - b.re, -b.im)
Base.:-(a::BiComplex) = BiComplex(-a.re, -a.im)

Base.angle(a::BiComplex) = atan(a.im / a.re)

function Base.convert(::Type{BiComplex{T}}, x::Real) where T
  return BiComplex(T(x) + 0im, zero(T) + 0im)
end
function Base.convert(::Type{BiComplex{T}}, x::Complex) where T
  return BiComplex(T(x), zero(T))
end
function Base.promote_rule(::Type{BiComplex{T}}, ::Type{U}) where {T, U<:Real}
  return BiComplex{promote(T, U)}
end

function Base.:*(a::BiComplex, b::BiComplex)
  return BiComplex(real(a) * real(b) - imag(a) * imag(b),
                   real(a) * imag(b) + imag(a) * real(b))
end
Base.:*(a::BiComplex, b::Number) = BiComplex(real(a) * b, imag(a) * b)
Base.:*(a::Number, b::BiComplex) = b * a

Base.:/(a::BiComplex, b::BiComplex) = a * inv(b)
Base.:/(a::Number, b::BiComplex) = a * inv(b)
Base.:/(a::BiComplex, b::Number) = BiComplex(a.re / b, a.im / b)

Base.exp(a::BiComplex) = exp(a.re) * BiComplex(cos(a.im), sin(a.im))
Base.cos(a::BiComplex) = (exp(imim * a) + exp(-imim * a)) / 2
Base.sin(a::BiComplex) = (exp(imim * a) - exp(-imim * a)) / 2 / imim

#αβ(a::Complex, b::Complex) = [1 -im; 1 im] * [a; b]
αβ(a::Complex, b::Complex) = (a - im * b, a + im * b)
αβ(a::BiComplex) = αβ(a.re, a.im)

reim(α::Complex, β::Complex) = ((α + β) / 2, (im * α - im * β) / 2)

const e⃗2 = BiComplex(1 + 0im, 0 + im)
const e⃗ᵀ2 = BiComplex(1 + 0im, 0 - im)
const e⃗ = e⃗2 / 2
const e⃗ᵀ = e⃗ᵀ2 / 2

function _op(op::T, a::BiComplex) where {T}
  α, β = αβ(a) 
  return (op(α) * e⃗2 + op(β) * e⃗ᵀ2) / 2
end
Base.:^(a::BiComplex, x::Float64) = _op(z -> z^x, a)
Base.sqrt(a::BiComplex) = _op(sqrt, a)
Base.cbrt(a::BiComplex) = _op(cbrt, a)
Base.acos(a::BiComplex) = _op(acos, a)
Base.asin(a::BiComplex) = _op(asin, a)
Base.log(a::BiComplex) = _op(log, a)

#function complexmodulus(a::BiComplex)
#  return sqrt(a.re^2 + a.im^2)
#end

function derivative(f::T, x::Complex{U}, h=sqrt(eps(U))) where {T,U<:Real}
  return imag(f(BiComplex(x, Complex(h, 0)))) / h
end

end


using .BiComplexes
using Test

@testset "+ - * /" begin
@testset "+ - * /" begin
  a = BiComplex(randn(ComplexF64), randn(ComplexF64))
  @test a - a == BiComplex(0 + 0im, 0 + 0im)
  @test a + a == 2 * a
  @test isapprox(a / a, BiComplex(1.0 + 0im, 0.0 + 0im), rtol=4eps(), atol=0)
  @test isapprox((a * a) / a / a, BiComplex(1.0 + 0im, 0.0 + 0im), rtol=4eps(), atol=0)
end

@testset "exp sin cos log" begin
  a = BiComplex(randn(ComplexF64), randn(ComplexF64))
  unit = BiComplex(1+0im, 0+0im)
  @test isapprox(sin(a)^2 + cos(a)^2, unit, rtol=sqrt(eps()), atol=10eps())
  @test isapprox(asin(sin(a)), a, rtol=sqrt(eps()), atol=10eps())
  @test isapprox(acos(cos(a)), a, rtol=sqrt(eps()), atol=10eps())
  @test isapprox(sin(asin(a)), a, rtol=sqrt(eps()), atol=10eps())
  @test isapprox(cos(acos(a)), a, rtol=sqrt(eps()), atol=10eps())
  @test isapprox(log(exp(a)), a, rtol=sqrt(eps()), atol=10eps())
  @test isapprox(exp(log(a)), a, rtol=sqrt(eps()), atol=10eps())
end

@testset "derivative" begin
  f(x) = sin(x) * exp(-x^2) / x^3
  g(x) = (cos(x) / x^3 - 2 * x * sin(x) / x^3 - 3 * sin(x) / x^4 ) * exp(-x^2)
  
  z = 1.0 + im
  @test isapprox(BiComplexes.derivative(f, z), g(z), rtol=sqrt(eps()),
                 atol=10eps())
  z = randn(ComplexF64) 
  @test isapprox(BiComplexes.derivative(f, z), g(z), rtol=sqrt(eps()),
                 atol=10eps())
end
end