Created
March 2, 2021 20:03
-
-
Save jwscook/e096706547684ef091dd8c3b98b611cb to your computer and use it in GitHub Desktop.
BiComplex numbers in Julia
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment