diff --git a/src/DualQuaternion.jl b/src/DualQuaternion.jl index e39ba8a7..334b7064 100644 --- a/src/DualQuaternion.jl +++ b/src/DualQuaternion.jl @@ -1,39 +1,53 @@ -using DualNumbers - struct DualQuaternion{T<:Real} <: Number q0::Quaternion{T} qe::Quaternion{T} - norm::Bool end -DualQuaternion{T}(dq::DualQuaternion) where {T<:Real} = DualQuaternion{T}(dq.q0, dq.qe, dq.norm) -function DualQuaternion{T}(d1::Dual, d2::Dual, d3::Dual, d4::Dual, n::Bool=false) where {T<:Real} +Base.@deprecate( + DualQuaternion{T}(q0::Quaternion{T}, qe::Quaternion{T}, norm::Bool) where {T<:Real}, + DualQuaternion{T}(q0, qe) +) +Base.@deprecate( + DualQuaternion{T}(d1::Dual, d2::Dual, d3::Dual, d4::Dual, norm::Bool) where {T<:Real}, + DualQuaternion{T}(d1, d2, d3, d4) +) +Base.@deprecate( + DualQuaternion(q0::Quaternion, qe::Quaternion, norm::Bool), + DualQuaternion(q0, qe) +) +Base.@deprecate( + DualQuaternion(d1::Dual, d2::Dual, d3::Dual, d4::Dual, norm::Bool), + DualQuaternion(d1, d2, d3, d4) +) +Base.@deprecate dualquat(q0, qe, n) dualquat(q0, qe) +Base.@deprecate dualquat(d1, d2, d3, d4, n) dualquat(d1, d2, d3, d4) + +DualQuaternion{T}(dq::DualQuaternion) where {T<:Real} = DualQuaternion{T}(dq.q0, dq.qe) +function DualQuaternion{T}(d1::Dual, d2::Dual, d3::Dual, d4::Dual) where {T<:Real} return DualQuaternion{T}( - Quaternion(DualNumbers.value(d1), DualNumbers.value(d2), DualNumbers.value(d3), DualNumbers.value(d4), n), + Quaternion(DualNumbers.value(d1), DualNumbers.value(d2), DualNumbers.value(d3), DualNumbers.value(d4)), Quaternion(DualNumbers.epsilon(d1), DualNumbers.epsilon(d2), DualNumbers.epsilon(d3), DualNumbers.epsilon(d4)), - n, ) end function DualQuaternion{T}(q0::Quaternion) where {T<:Real} - return DualQuaternion{T}(convert(Quaternion{T}, q0), zero(Quaternion{T}), q0.norm) + return DualQuaternion{T}(convert(Quaternion{T}, q0), zero(Quaternion{T})) end function DualQuaternion{T}(d::Dual) where {T<:Real} return DualQuaternion( Quaternion{T}(DualNumbers.value(d)), - Quaternion{T}(DualNumbers.epsilon(d)), - (DualNumbers.value(d)==one(DualNumbers.value(d))) & iszero(DualNumbers.epsilon(d))) + Quaternion{T}(DualNumbers.epsilon(d))) end function DualQuaternion{T}(x::Real) where {T<:Real} - return DualQuaternion(convert(Quaternion{T}, x), zero(Quaternion{T}), abs(x) == one(x)) + return DualQuaternion(convert(Quaternion{T}, x), zero(Quaternion{T})) end -DualQuaternion(q0::Quaternion, qe::Quaternion) = DualQuaternion(promote(q0, qe)..., false) -DualQuaternion(d1::Dual, d2::Dual, d3::Dual, d4::Dual, n::Bool=false) = - DualQuaternion(Quaternion(DualNumbers.value(d1), DualNumbers.value(d2), DualNumbers.value(d3), DualNumbers.value(d4), n), - Quaternion(DualNumbers.epsilon(d1), DualNumbers.epsilon(d2), DualNumbers.epsilon(d3), DualNumbers.epsilon(d4)), n) -DualQuaternion(x::Real) = DualQuaternion(Quaternion(x), Quaternion(zero(x)), abs(x) == one(x)) -DualQuaternion(d::Dual) = DualQuaternion(Quaternion(DualNumbers.value(d)), Quaternion(DualNumbers.epsilon(d)), (DualNumbers.value(d)==one(DualNumbers.value(d))) & iszero(DualNumbers.epsilon(d))) -DualQuaternion(q::Quaternion) = DualQuaternion(q, zero(q), q.norm) +DualQuaternion(q0::Quaternion, qe::Quaternion) = DualQuaternion(promote(q0, qe)...) +DualQuaternion(d1::Dual, d2::Dual, d3::Dual, d4::Dual) = + DualQuaternion(Quaternion(DualNumbers.value(d1), DualNumbers.value(d2), DualNumbers.value(d3), DualNumbers.value(d4)), + Quaternion(DualNumbers.epsilon(d1), DualNumbers.epsilon(d2), DualNumbers.epsilon(d3), DualNumbers.epsilon(d4))) +DualQuaternion(x::Real) = DualQuaternion(Quaternion(x), Quaternion(zero(x))) +DualQuaternion(d::Dual) = DualQuaternion(Quaternion(DualNumbers.value(d)), Quaternion(DualNumbers.epsilon(d))) +DualQuaternion(q::Quaternion) = DualQuaternion(q, zero(q)) DualQuaternion(a::Vector) = DualQuaternion(zero(Quaternion{typeof(a[1])}), Quaternion(a)) const DualQuaternionF16 = DualQuaternion{Float16} @@ -46,11 +60,17 @@ promote_rule(::Type{DualQuaternion{T}}, ::Type{Quaternion{S}}) where {T <: Real, promote_rule(::Type{DualQuaternion{T}}, ::Type{DualQuaternion{S}}) where {T <: Real, S <: Real} = DualQuaternion{promote_type(T, S)} dualquat(q1, q2) = DualQuaternion(q1, q2) -dualquat(q1, q2, n) = DualQuaternion(q1, q2, n) dualquat(d1, d2, d3, d4) = DualQuaternion(d1, d2, d3, d4) -dualquat(d1, d2, d3, d4, n) = DualQuaternion(d1, d2, d3, d4, n) dualquat(x) = DualQuaternion(x) +function Base.getproperty(dq::DualQuaternion, k::Symbol) + if k === :norm + Base.depwarn("`dq.norm` is deprecated. Please use `isunit(dq)` instead.", Symbol("dq.norm")) + return isunit(dq) + end + return getfield(dq, k) +end + function show(io::IO, dq::DualQuaternion) show(io, dq.q0) print(io, " + ( ") @@ -69,58 +89,57 @@ Qe(dq::DualQuaternion) = dq.qe dual(dq.q0.v2, dq.qe.v2) / d, dual(dq.q0.v3, dq.qe.v3) / d) -abs2(dq::DualQuaternion) = dq.norm ? dual(one(dq.q0.s)) : +abs2(dq::DualQuaternion) = isunit(dq) ? dual(one(dq.q0.s)) : dual(abs2(dq.q0), 2.0 * (dq.q0.s * dq.qe.s + dq.q0.v1 * dq.qe.v1 + dq.q0.v2 * dq.qe.v2 + dq.q0.v3 * dq.qe.v3)) -abs(dq::DualQuaternion) = dq.norm ? dual(one(dq.q0.s)) : sqrt(abs2(dq)) +abs(dq::DualQuaternion) = isunit(dq) ? dual(one(dq.q0.s)) : sqrt(abs2(dq)) float(dq::DualQuaternion{T}) where T = convert(DualQuaternion{float(T)}, dq) -conj(dq::DualQuaternion) = DualQuaternion(conj(dq.q0), conj(dq.qe), dq.norm) -dconj(dq::DualQuaternion) = DualQuaternion(dq.q0, -dq.qe, dq.norm) +conj(dq::DualQuaternion) = DualQuaternion(conj(dq.q0), conj(dq.qe)) +dconj(dq::DualQuaternion) = DualQuaternion(dq.q0, -dq.qe) -inv(dq::DualQuaternion) = dq.norm ? conj(dq) : conj(dq) / abs2(dq) +inv(dq::DualQuaternion) = isunit(dq) ? conj(dq) : conj(dq) / abs2(dq) function normalize(dq::DualQuaternion) - if (dq.norm) + if (isunit(dq)) return dq end a = abs(dq) if abs(a) > 0 qa = dq / a - dualquat(qa.q0, qa.qe, true) + dualquat(qa.q0, qa.qe) else dq end end function normalizea(dq::DualQuaternion{T}) where {T} - if (dq.norm) + if (isunit(dq)) return (dq, one(Dual{T})) end a = abs(dq) if abs(a) > 0 qa = dq / a - dualquat(qa.q0, qa.qe, true), a + dualquat(qa.q0, qa.qe), a else dq, zero(Dual{T}) end end -(-)(dq::DualQuaternion) = DualQuaternion(-dq.q0, -dq.qe, dq.norm) +(-)(dq::DualQuaternion) = DualQuaternion(-dq.q0, -dq.qe) (+)(dq::DualQuaternion, dw::DualQuaternion) = DualQuaternion(dq.q0 + dw.q0, dq.qe + dw.qe) (-)(dq::DualQuaternion, dw::DualQuaternion) = DualQuaternion(dq.q0 - dw.q0, dq.qe - dw.qe) (*)(dq::DualQuaternion, dw::DualQuaternion) = DualQuaternion(dq.q0 * dw.q0, - dq.q0 * dw.qe + dq.qe * dw.q0, - dq.norm && dw.norm) + dq.q0 * dw.qe + dq.qe * dw.q0) (*)(dq::DualQuaternion, d::Dual) = (*)(Base.promote(dq, d)...) (*)(d::Dual, dq::DualQuaternion) = (*)(Base.promote(d, dq)...) (/)(dq::DualQuaternion, dw::DualQuaternion) = dq * inv(dw) -(==)(q::DualQuaternion, w::DualQuaternion) = (q.q0 == w.q0) & (q.qe == w.qe) # ignore .norm field +(==)(q::DualQuaternion, w::DualQuaternion) = (q.q0 == w.q0) & (q.qe == w.qe) function angleaxis(dq::DualQuaternion) tq = dq.qe * conj(dq.q0) @@ -156,7 +175,7 @@ function exp(dq::DualQuaternion) se = exp(se) dq = dualquat(quat(0.0, imag_part(dq.q0)...), quat(0.0, imag_part(dq.qe)...)) dq, th = normalizea(dq) - if dq.norm + if isunit(dq) dualquat(se) * (dualquat(cos(th)) + dq * dualquat(sin(th))) else dualquat(se) @@ -180,5 +199,5 @@ dualquatrand() = dualquat(quatrand(), quatrand()) ndualquatrand() = normalize(dualquatrand()) function rand(rng::AbstractRNG, ::Random.SamplerType{DualQuaternion{T}}) where {T<:Real} - return DualQuaternion{T}(rand(rng, Quaternion{T}), rand(rng, Quaternion{T}), false) + return DualQuaternion{T}(rand(rng, Quaternion{T}), rand(rng, Quaternion{T})) end diff --git a/src/Octonion.jl b/src/Octonion.jl index 52536e55..4a0c3941 100644 --- a/src/Octonion.jl +++ b/src/Octonion.jl @@ -7,20 +7,19 @@ struct Octonion{T<:Real} <: Number v5::T v6::T v7::T - norm::Bool end Octonion{T}(x::Real) where {T<:Real} = Octonion(convert(T, x)) Octonion{T}(x::Complex) where {T<:Real} = Octonion(convert(Complex{T}, x)) Octonion{T}(q::Quaternion) where {T<:Real} = Octonion(convert(Quaternion{T}, q)) Octonion{T}(o::Octonion) where {T<:Real} = - Octonion{T}(o.s, o.v1, o.v2, o.v3, o.v4, o.v5, o.v6, o.v7, o.norm) + Octonion{T}(o.s, o.v1, o.v2, o.v3, o.v4, o.v5, o.v6, o.v7) -Octonion(s::Real, v1::Real, v2::Real, v3::Real, v4::Real, v5::Real, v6::Real, v7::Real, n::Bool = false) = - Octonion(promote(s, v1, v2, v3, v4, v5, v6, v7)..., n) -Octonion(x::Real) = Octonion(x, zero(x), zero(x), zero(x), zero(x), zero(x), zero(x), zero(x), abs(x) == one(x)) -Octonion(z::Complex) = Octonion(z.re, z.im, zero(z.re), zero(z.re), zero(z.re), zero(z.re), zero(z.re), zero(z.re), abs(z) == one(z.re)) -Octonion(q::Quaternion) = Octonion(q.s, q.v1, q.v2, q.v3, zero(q.s), zero(q.s), zero(q.s), zero(q.s), q.norm) +Octonion(s::Real, v1::Real, v2::Real, v3::Real, v4::Real, v5::Real, v6::Real, v7::Real) = + Octonion(promote(s, v1, v2, v3, v4, v5, v6, v7)...) +Octonion(x::Real) = Octonion(x, zero(x), zero(x), zero(x), zero(x), zero(x), zero(x), zero(x)) +Octonion(z::Complex) = Octonion(z.re, z.im, zero(z.re), zero(z.re), zero(z.re), zero(z.re), zero(z.re), zero(z.re)) +Octonion(q::Quaternion) = Octonion(q.s, q.v1, q.v2, q.v3, zero(q.s), zero(q.s), zero(q.s), zero(q.s)) Octonion(s::Real, a::Vector) = Octonion(s, a[1], a[2], a[3], a[4], a[5], a[6], a[7]) Octonion(a::Vector) = Octonion(0, a[1], a[2], a[3], a[4], a[5], a[6], a[7]) @@ -34,10 +33,27 @@ promote_rule(::Type{Octonion{T}}, ::Type{Quaternion{S}}) where {T <: Real, S <: promote_rule(::Type{Octonion{T}}, ::Type{Octonion{S}}) where {T <: Real, S <: Real} = Octonion{promote_type(T, S)} octo(p, v1, v2, v3, v4, v5, v6, v7) = Octonion(p, v1, v2, v3, v4, v5, v6, v7) -octo(p, v1, v2, v3, v4, v5, v6, v7, n) = Octonion(p, v1, v2, v3, v4, v5, v6, v7, n) octo(x) = Octonion(x) octo(s, a) = Octonion(s, a) +Base.@deprecate( + Octonion{T}(s::Real, v1::Real, v2::Real, v3::Real, v4::Real, v5::Real, v6::Real, v7::Real, norm::Bool) where {T<:Real}, + Octonion{T}(s, v1, v2, v3, v4, v5, v6, v7) +) +Base.@deprecate( + Octonion(s::Real, v1::Real, v2::Real, v3::Real, v4::Real, v5::Real, v6::Real, v7::Real, norm::Bool), + Octonion(s, v1, v2, v3, v4, v5, v6, v7) +) +Base.@deprecate octo(p, v1, v2, v3, v4, v5, v6, v7, n) octo(p, v1, v2, v3, v4, v5, v6, v7) + +function Base.getproperty(o::Octonion, k::Symbol) + if k === :norm + Base.depwarn("`o.norm` is deprecated. Please use `isunit(o)` instead.", Symbol("o.norm")) + return isunit(o) + end + return getfield(o, k) +end + function show(io::IO, o::Octonion) pm(x) = x < 0 ? " - $(-x)" : " + $x" print(io, o.s, pm(o.v1), "im", pm(o.v2), "jm", pm(o.v3), "km", pm(o.v4), "ilm", pm(o.v5), "jlm", pm(o.v6), "klm", pm(o.v7), "lm") @@ -49,37 +65,37 @@ imag_part(o::Octonion) = (o.v1, o.v2, o.v3, o.v4, o.v5, o.v6, o.v7) (/)(o::Octonion, x::Real) = Octonion(o.s / x, o.v1 / x, o.v2 / x, o.v3 / x, o.v4 / x, o.v5 / x, o.v6 / x, o.v7 / x) -conj(o::Octonion) = Octonion(o.s, -o.v1, -o.v2, -o.v3, -o.v4, -o.v5, -o.v6, -o.v7, o.norm) +conj(o::Octonion) = Octonion(o.s, -o.v1, -o.v2, -o.v3, -o.v4, -o.v5, -o.v6, -o.v7) abs(o::Octonion) = sqrt(o.s * o.s + o.v1 * o.v1 + o.v2 * o.v2 + o.v3 * o.v3 + o.v4 * o.v4 + o.v5 * o.v5 + o.v6 * o.v6 + o.v7 * o.v7) float(q::Octonion{T}) where T = convert(Octonion{float(T)}, q) abs_imag(o::Octonion) = sqrt(o.v1 * o.v1 + o.v2 * o.v2 + o.v3 * o.v3 + o.v4 * o.v4 + o.v5 * o.v5 + o.v6 * o.v6 + o.v7 * o.v7) abs2(o::Octonion) = o.s * o.s + o.v1 * o.v1 + o.v2 * o.v2 + o.v3 * o.v3 + o.v4 * o.v4 + o.v5 * o.v5 + o.v6 * o.v6 + o.v7 * o.v7 -inv(o::Octonion) = o.norm ? conj(o) : conj(o) / abs2(o) +inv(o::Octonion) = isunit(o) ? conj(o) : conj(o) / abs2(o) isreal(o::Octonion) = iszero(o.v1) & iszero(o.v2) & iszero(o.v3) & iszero(o.v4) & iszero(o.v5) & iszero(o.v6) & iszero(o.v7) -isfinite(o::Octonion) = o.norm | (isfinite(real(o)) & isfinite(o.v1) & isfinite(o.v2) & isfinite(o.v3) & isfinite(o.v4) & isfinite(o.v5) & isfinite(o.v6) & isfinite(o.v7)) -iszero(o::Octonion) = ~o.norm & iszero(real(o)) & iszero(o.v1) & iszero(o.v2) & iszero(o.v3) & iszero(o.v4) & iszero(o.v5) & iszero(o.v6) & iszero(o.v7) +isfinite(o::Octonion) = isfinite(real(o)) & isfinite(o.v1) & isfinite(o.v2) & isfinite(o.v3) & isfinite(o.v4) & isfinite(o.v5) & isfinite(o.v6) & isfinite(o.v7) +iszero(o::Octonion) = iszero(real(o)) & iszero(o.v1) & iszero(o.v2) & iszero(o.v3) & iszero(o.v4) & iszero(o.v5) & iszero(o.v6) & iszero(o.v7) isnan(o::Octonion) = isnan(real(o)) | isnan(o.v1) | isnan(o.v2) | isnan(o.v3) | isnan(o.v4) | isnan(o.v5) | isnan(o.v6) | isnan(o.v7) -isinf(o::Octonion) = ~o.norm & (isinf(real(o)) | isinf(o.v1) | isinf(o.v2) | isinf(o.v3) | isinf(o.v4) | isinf(o.v5) | isinf(o.v6) | isinf(o.v7)) +isinf(o::Octonion) = isinf(real(o)) | isinf(o.v1) | isinf(o.v2) | isinf(o.v3) | isinf(o.v4) | isinf(o.v5) | isinf(o.v6) | isinf(o.v7) function normalize(o::Octonion) - if (o.norm) + if (isunit(o)) return o end o = o / abs(o) - Octonion(o.s, o.v1, o.v2, o.v3, o.v4, o.v5, o.v6, o.v7, true) + Octonion(o.s, o.v1, o.v2, o.v3, o.v4, o.v5, o.v6, o.v7) end function normalizea(o::Octonion) - if (o.norm) + if (isunit(o)) return (o, one(o.s)) end a = abs(o) o = o / a - (Octonion(o.s, o.v1, o.v2, o.v3, o.v4, o.v5, o.v6, o.v7, true), a) + (Octonion(o.s, o.v1, o.v2, o.v3, o.v4, o.v5, o.v6, o.v7), a) end -(-)(o::Octonion) = Octonion(-o.s, -o.v1, -o.v2, -o.v3, -o.v4, -o.v5, -o.v6, -o.v7, o.norm) +(-)(o::Octonion) = Octonion(-o.s, -o.v1, -o.v2, -o.v3, -o.v4, -o.v5, -o.v6, -o.v7) (+)(o::Octonion, w::Octonion) = Octonion(o.s + w.s, o.v1 + w.v1, @@ -121,7 +137,7 @@ end (/)(o::Octonion, w::Octonion) = o * inv(w) (==)(q::Octonion, w::Octonion) = (q.s == w.s) & (q.v1 == w.v1) & (q.v2 == w.v2) & (q.v3 == w.v3) & - (q.v4 == w.v4) & (q.v5 == w.v5) & (q.v6 == w.v6) & (q.v7 == w.v7) # ignore .norm field + (q.v4 == w.v4) & (q.v5 == w.v5) & (q.v6 == w.v6) & (q.v7 == w.v7) function exp(o::Octonion) s = o.s @@ -138,8 +154,7 @@ function exp(o::Octonion) scale * o.v4, scale * o.v5, scale * o.v6, - scale * o.v7, - abs(s) == 0) + scale * o.v7) end function log(o::Octonion) @@ -172,7 +187,7 @@ octorand() = octo(randn(), randn(), randn(), randn(), randn(), randn(), randn(), function rand(rng::AbstractRNG, ::Random.SamplerType{Octonion{T}}) where {T<:Real} Octonion{T}(rand(rng, T), rand(rng, T), rand(rng, T), rand(rng, T), - rand(rng, T), rand(rng, T), rand(rng, T), rand(rng, T), false) + rand(rng, T), rand(rng, T), rand(rng, T), rand(rng, T)) end function randn(rng::AbstractRNG, ::Type{Octonion{T}}) where {T<:AbstractFloat} @@ -185,6 +200,5 @@ function randn(rng::AbstractRNG, ::Type{Octonion{T}}) where {T<:AbstractFloat} randn(rng, T) * INV_SQRT_EIGHT, randn(rng, T) * INV_SQRT_EIGHT, randn(rng, T) * INV_SQRT_EIGHT, - false, ) end diff --git a/src/Quaternion.jl b/src/Quaternion.jl index 990bc3a0..6b171328 100644 --- a/src/Quaternion.jl +++ b/src/Quaternion.jl @@ -3,7 +3,6 @@ struct Quaternion{T<:Real} <: Number v1::T v2::T v3::T - norm::Bool end const QuaternionF16 = Quaternion{Float16} @@ -12,11 +11,11 @@ const QuaternionF64 = Quaternion{Float64} Quaternion{T}(x::Real) where {T<:Real} = Quaternion(convert(T, x)) Quaternion{T}(x::Complex) where {T<:Real} = Quaternion(convert(Complex{T}, x)) -Quaternion{T}(q::Quaternion) where {T<:Real} = Quaternion{T}(q.s, q.v1, q.v2, q.v3, q.norm) -Quaternion(s::Real, v1::Real, v2::Real, v3::Real, n::Bool = false) = - Quaternion(promote(s, v1, v2, v3)..., n) -Quaternion(x::Real) = Quaternion(x, zero(x), zero(x), zero(x), abs(x) == one(x)) -Quaternion(z::Complex) = Quaternion(z.re, z.im, zero(z.re), zero(z.re), abs(z) == one(z.re)) +Quaternion{T}(q::Quaternion) where {T<:Real} = Quaternion{T}(q.s, q.v1, q.v2, q.v3) +Quaternion(s::Real, v1::Real, v2::Real, v3::Real) = + Quaternion(promote(s, v1, v2, v3)...) +Quaternion(x::Real) = Quaternion(x, zero(x), zero(x), zero(x)) +Quaternion(z::Complex) = Quaternion(z.re, z.im, zero(z.re), zero(z.re)) Quaternion(s::Real, a::AbstractVector) = Quaternion(s, a[1], a[2], a[3]) Quaternion(a::AbstractVector) = Quaternion(0, a[1], a[2], a[3]) @@ -25,10 +24,27 @@ promote_rule(::Type{Quaternion{T}}, ::Type{Complex{S}}) where {T <: Real, S <: R promote_rule(::Type{Quaternion{T}}, ::Type{Quaternion{S}}) where {T <: Real, S <: Real} = Quaternion{promote_type(T, S)} quat(p, v1, v2, v3) = Quaternion(p, v1, v2, v3) -quat(p, v1, v2, v3, n) = Quaternion(p, v1, v2, v3, n) quat(x) = Quaternion(x) quat(s, a) = Quaternion(s, a) +Base.@deprecate( + Quaternion{T}(s::Real, v1::Real, v2::Real, v3::Real, norm::Bool) where {T<:Real}, + Quaternion{T}(s, v1, v2, v3) +) +Base.@deprecate( + Quaternion(s::Real, v1::Real, v2::Real, v3::Real, norm::Bool), + Quaternion(s, v1, v2, v3) +) +Base.@deprecate quat(p, v1, v2, v3, n) quat(p, v1, v2, v3) + +function Base.getproperty(q::Quaternion, k::Symbol) + if k === :norm + Base.depwarn("`q.norm` is deprecated. Please use `isunit(q)` instead.", Symbol("q.norm")) + return isunit(q) + end + return getfield(q, k) +end + function show(io::IO, q::Quaternion) pm(x) = x < 0 ? " - $(-x)" : " + $x" print(io, q.s, pm(q.v1), "im", pm(q.v2), "jm", pm(q.v3), "km") @@ -41,47 +57,47 @@ imag_part(q::Quaternion) = (q.v1, q.v2, q.v3) (/)(q::Quaternion, x::Real) = Quaternion(q.s / x, q.v1 / x, q.v2 / x, q.v3 / x) -conj(q::Quaternion) = Quaternion(q.s, -q.v1, -q.v2, -q.v3, q.norm) +conj(q::Quaternion) = Quaternion(q.s, -q.v1, -q.v2, -q.v3) abs(q::Quaternion) = sqrt(q.s * q.s + q.v1 * q.v1 + q.v2 * q.v2 + q.v3 * q.v3) float(q::Quaternion{T}) where T = convert(Quaternion{float(T)}, q) abs_imag(q::Quaternion) = sqrt(q.v1 * q.v1 + q.v2 * q.v2 + q.v3 * q.v3) abs2(q::Quaternion) = q.s * q.s + q.v1 * q.v1 + q.v2 * q.v2 + q.v3 * q.v3 -inv(q::Quaternion) = q.norm ? conj(q) : conj(q) / abs2(q) +inv(q::Quaternion) = isunit(q) ? conj(q) : conj(q) / abs2(q) isreal(q::Quaternion) = iszero(q.v1) & iszero(q.v2) & iszero(q.v3) -isfinite(q::Quaternion) = q.norm | (isfinite(q.s) & isfinite(q.v1) & isfinite(q.v2) & isfinite(q.v3)) -iszero(q::Quaternion) = ~q.norm & iszero(real(q)) & iszero(q.v1) & iszero(q.v2) & iszero(q.v3) +isfinite(q::Quaternion) = isfinite(q.s) & isfinite(q.v1) & isfinite(q.v2) & isfinite(q.v3) +iszero(q::Quaternion) = iszero(real(q)) & iszero(q.v1) & iszero(q.v2) & iszero(q.v3) isnan(q::Quaternion) = isnan(real(q)) | isnan(q.v1) | isnan(q.v2) | isnan(q.v3) -isinf(q::Quaternion) = ~q.norm & (isinf(q.s) | isinf(q.v1) | isinf(q.v2) | isinf(q.v3)) +isinf(q::Quaternion) = isinf(q.s) | isinf(q.v1) | isinf(q.v2) | isinf(q.v3) function normalize(q::Quaternion) - if (q.norm) + if (isunit(q)) return q end q = q / abs(q) - Quaternion(q.s, q.v1, q.v2, q.v3, true) + Quaternion(q.s, q.v1, q.v2, q.v3) end function normalizea(q::Quaternion) - if (q.norm) + if (isunit(q)) return (q, one(q.s)) end a = abs(q) q = q / a - (Quaternion(q.s, q.v1, q.v2, q.v3, true), a) + (Quaternion(q.s, q.v1, q.v2, q.v3), a) end function normalizeq(q::Quaternion) a = abs(q) if a > 0 q = q / a - Quaternion(q.s, q.v1, q.v2, q.v3, true) + Quaternion(q.s, q.v1, q.v2, q.v3) else - Quaternion(0.0, 1.0, 0.0, 0.0, true) + Quaternion(0.0, 1.0, 0.0, 0.0) end end -(-)(q::Quaternion) = Quaternion(-q.s, -q.v1, -q.v2, -q.v3, q.norm) +(-)(q::Quaternion) = Quaternion(-q.s, -q.v1, -q.v2, -q.v3) (+)(q::Quaternion, w::Quaternion) = Quaternion(q.s + w.s, q.v1 + w.v1, q.v2 + w.v2, q.v3 + w.v3) @@ -92,11 +108,10 @@ end (*)(q::Quaternion, w::Quaternion) = Quaternion(q.s * w.s - q.v1 * w.v1 - q.v2 * w.v2 - q.v3 * w.v3, q.s * w.v1 + q.v1 * w.s + q.v2 * w.v3 - q.v3 * w.v2, q.s * w.v2 - q.v1 * w.v3 + q.v2 * w.s + q.v3 * w.v1, - q.s * w.v3 + q.v1 * w.v2 - q.v2 * w.v1 + q.v3 * w.s, - q.norm && w.norm) + q.s * w.v3 + q.v1 * w.v2 - q.v2 * w.v1 + q.v3 * w.s) (/)(q::Quaternion, w::Quaternion) = q * inv(w) -(==)(q::Quaternion, w::Quaternion) = (q.s == w.s) & (q.v1 == w.v1) & (q.v2 == w.v2) & (q.v3 == w.v3) # ignore .norm field +(==)(q::Quaternion, w::Quaternion) = (q.s == w.s) & (q.v1 == w.v1) & (q.v2 == w.v2) & (q.v3 == w.v3) angleaxis(q::Quaternion) = angle(q), axis(q) @@ -140,20 +155,16 @@ function extend_analytic(f, q::Quaternion) w = f(z) wr, wi = reim(w) scale = wi / a - norm = _isexpfun(f) && iszero(s) if a > 0 - return Quaternion(wr, scale * q.v1, scale * q.v2, scale * q.v3, norm) + return Quaternion(wr, scale * q.v1, scale * q.v2, scale * q.v3) else # q == real(q), so f(real(q)) may be real or complex, i.e. wi may be nonzero. # we choose to embed complex numbers in the quaternions by identifying the first # imaginary quaternion basis with the complex imaginary basis. - return Quaternion(wr, oftype(scale, wi), zero(scale), zero(scale), norm) + return Quaternion(wr, oftype(scale, wi), zero(scale), zero(scale)) end end -_isexpfun(::Union{typeof(exp),typeof(exp2),typeof(exp10)}) = true -_isexpfun(::Any) = false - """ cis(q::Quaternion) @@ -245,14 +256,14 @@ function linpol(p::Quaternion, q::Quaternion, t::Real) Quaternion(sp * p.s + sq * q.s, sp * p.v1 + sq * q.v1, sp * p.v2 + sq * q.v2, - sp * p.v3 + sq * q.v3, true) + sp * p.v3 + sq * q.v3) end quatrand(rng = Random.GLOBAL_RNG) = quat(randn(rng), randn(rng), randn(rng), randn(rng)) nquatrand(rng = Random.GLOBAL_RNG) = normalize(quatrand(rng)) function rand(rng::AbstractRNG, ::Random.SamplerType{Quaternion{T}}) where {T<:Real} - Quaternion{T}(rand(rng, T), rand(rng, T), rand(rng, T), rand(rng, T), false) + Quaternion{T}(rand(rng, T), rand(rng, T), rand(rng, T), rand(rng, T)) end function randn(rng::AbstractRNG, ::Type{Quaternion{T}}) where {T<:AbstractFloat} @@ -261,7 +272,6 @@ function randn(rng::AbstractRNG, ::Type{Quaternion{T}}) where {T<:AbstractFloat} randn(rng, T) * 1//2, randn(rng, T) * 1//2, randn(rng, T) * 1//2, - false, ) end @@ -278,7 +288,7 @@ function qrotation(axis::AbstractVector{T}, theta) where {T <: Real} end s,c = sincos(theta / 2) scaleby = s / normaxis - Quaternion(c, scaleby * axis[1], scaleby * axis[2], scaleby * axis[3], true) + Quaternion(c, scaleby * axis[1], scaleby * axis[2], scaleby * axis[3]) end # Variant of the above where norm(rotvec) encodes theta @@ -289,7 +299,7 @@ function qrotation(rotvec::AbstractVector{T}) where {T <: Real} theta = norm(rotvec) s,c = sincos(theta / 2) scaleby = s / (iszero(theta) ? one(theta) : theta) - Quaternion(c, scaleby * rotvec[1], scaleby * rotvec[2], scaleby * rotvec[3], true) + Quaternion(c, scaleby * rotvec[1], scaleby * rotvec[2], scaleby * rotvec[3]) end function qrotation(dcm::AbstractMatrix{T}) where {T<:Real} diff --git a/src/Quaternions.jl b/src/Quaternions.jl index 6eeed774..e211fc0b 100644 --- a/src/Quaternions.jl +++ b/src/Quaternions.jl @@ -9,9 +9,22 @@ module Quaternions import LinearAlgebra: lyap, norm, normalize, sylvester using LinearAlgebra: cross, dot using Random + using DualNumbers Base.@irrational INV_SQRT_EIGHT 0.3535533905932737622004 sqrt(big(0.125)) + """ + isunit(x) + + Return `true` if for the hypercomplex number `x`, `abs(x) == 1`. + """ + function isunit end + isunit(x::Real) = isone(abs(x)) + isunit(x::Complex) = isone(abs(x)) + isunit(x::DualNumbers.Dual) = isunit(DualNumbers.value(x)) & iszero(DualNumbers.epsilon(x)) + # TODO: replace this with something more efficient + isunit(x::Number) = isone(abs(x)) + include("Quaternion.jl") include("Octonion.jl") include("DualQuaternion.jl") @@ -32,6 +45,7 @@ module Quaternions export octo export imag_part export dualquat + export isunit export angleaxis export angle export axis