Add missing tests

Project.toml

@@ -1,6 +1,6 @@
 name = "Quaternions"
 uuid = "94ee1d12-ae83-5a48-8b1c-48b8ff168ae0"
-version = "0.5.1"
+version = "0.5.2"
 
 [deps]
 DualNumbers = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74"
@@ -12,7 +12,8 @@ DualNumbers = "0.5, 0.6"
 julia = "1"
 
 [extras]
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test"]
+test = ["ForwardDiff", "Test"]

src/DualQuaternion.jl

@@ -6,43 +6,43 @@ struct DualQuaternion{T<:Real} <: Number
   norm::Bool
 end
 
-DualQuaternion(q0::Quaternion, qe::Quaternion, n::Bool = false) =
-  DualQuaternion(promote(q0, qe)..., n)
+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}
+  return DualQuaternion{T}(
+    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,
+  )
+end
+function DualQuaternion{T}(q0::Quaternion) where {T<:Real}
+  return DualQuaternion{T}(convert(Quaternion{T}, q0), zero(Quaternion{T}), q0.norm)
+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)))
+end
+function DualQuaternion{T}(x::Real) where {T<:Real}
+  return DualQuaternion(convert(Quaternion{T}, x), zero(Quaternion{T}), abs(x) == one(x))
+end
 
-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(d1::Dual, d2::Dual, d3::Dual, d4::Dual, n::Bool) =
+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(a::Vector) = DualQuaternion(zero(Quaternion{typeof(a[1])}), Quaternion(a))
 
 const DualQuaternionF16 = DualQuaternion{Float16}
 const DualQuaternionF32 = DualQuaternion{Float32}
 const DualQuaternionF64 = DualQuaternion{Float64}
 
-convert(::Type{DualQuaternion{T}}, x::Real) where {T} = DualQuaternion(convert(T, x))
-
-convert(::Type{DualQuaternion{T}}, d::Dual) where {T} = DualQuaternion(convert(Dual{T}, d))
-
-convert(::Type{DualQuaternion{T}}, q::Quaternion) where {T} = DualQuaternion(convert(Quaternion{T}, q))
-
-convert(::Type{DualQuaternion{T}}, q::DualQuaternion{T}) where {T <: Real} = q
-
-convert(::Type{DualQuaternion{T}}, dq::DualQuaternion) where {T} =
-    DualQuaternion(convert(Quaternion{T}, dq.q0), convert(Quaternion{T}, dq.qe), dq.norm)
-
-promote_rule(::Type{DualQuaternion{T}}, ::Type{T}) where {T <: Real} = DualQuaternion{T}
-promote_rule(::Type{DualQuaternion}, ::Type{T}) where {T <: Real} = DualQuaternion
 promote_rule(::Type{DualQuaternion{T}}, ::Type{S}) where {T <: Real, S <: Real} = DualQuaternion{promote_type(T, S)}
-promote_rule(::Type{Quaternion{T}}, ::Type{DualQuaternion{S}}) where {T <: Real, S <: Real} = DualQuaternion{promote_type(T, S)}
+promote_rule(::Type{DualQuaternion{T}}, ::Type{Dual{S}}) where {T <: Real, S <: Real} = DualQuaternion{promote_type(T, S)}
+promote_rule(::Type{DualQuaternion{T}}, ::Type{Quaternion{S}}) where {T <: Real, S <: Real} = DualQuaternion{promote_type(T, S)}
 promote_rule(::Type{DualQuaternion{T}}, ::Type{DualQuaternion{S}}) where {T <: Real, S <: Real} = DualQuaternion{promote_type(T, S)}
 
 dualquat(q1, q2) = DualQuaternion(q1, q2)
@@ -97,16 +97,16 @@ function normalize(dq::DualQuaternion)
   end
 end
 
-function normalizea(dq::DualQuaternion)
+function normalizea(dq::DualQuaternion{T}) where {T}
   if (dq.norm)
-    return (dq, one(dual))
+    return (dq, one(Dual{T}))
   end
   a = abs(dq)
   if abs(a) > 0
     qa = dq / a
     dualquat(qa.q0, qa.qe, true), a
   else
-    dq, zero(dual)
+    dq, zero(Dual{T})
   end
 end
 
@@ -117,6 +117,8 @@ end
 (*)(dq::DualQuaternion, dw::DualQuaternion) = DualQuaternion(dq.q0 * dw.q0,
                                                               dq.q0 * dw.qe + dq.qe * dw.q0,
                                                               dq.norm && dw.norm)
+(*)(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
 

src/Octonion.jl

@@ -10,6 +10,12 @@ struct Octonion{T<:Real} <: Number
   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(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))
@@ -22,18 +28,9 @@ const OctonionF16 = Octonion{Float16}
 const OctonionF32 = Octonion{Float32}
 const OctonionF64 = Octonion{Float64}
 
-convert(::Type{Octonion{T}}, x::Real) where {T} = Octonion(convert(T, x))
-convert(::Type{Octonion{T}}, z::Complex) where {T} = Octonion(convert(Complex{T}, z))
-convert(::Type{Octonion{T}}, q::Quaternion) where {T} = Octonion(convert(Quaternion{T}, q))
-convert(::Type{Octonion{T}}, o::Octonion{T}) where {T <: Real} = o
-convert(::Type{Octonion{T}}, o::Octonion) where {T} =
-  Octonion(convert(T, o.s), convert(T, o.v1), convert(T, o.v2), convert(T, o.v3), convert(T, o.v4), convert(T, o.v5), convert(T, o.v6), convert(T, o.v7), o.norm)
-
-promote_rule(::Type{Octonion{T}}, ::Type{T}) where {T <: Real} = Octonion{T}
-promote_rule(::Type{Octonion}, ::Type{T}) where {T <: Real} = Octonion
 promote_rule(::Type{Octonion{T}}, ::Type{S}) where {T <: Real, S <: Real} = Octonion{promote_type(T, S)}
-promote_rule(::Type{Complex{T}}, ::Type{Octonion{S}}) where {T <: Real, S <: Real} = Octonion{promote_type(T, S)}
-promote_rule(::Type{Quaternion{T}}, ::Type{Octonion{S}}) where {T <: Real, S <: Real} = Octonion{promote_type(T, S)}
+promote_rule(::Type{Octonion{T}}, ::Type{Complex{S}}) where {T <: Real, S <: Real} = Octonion{promote_type(T, S)}
+promote_rule(::Type{Octonion{T}}, ::Type{Quaternion{S}}) where {T <: Real, S <: Real} = Octonion{promote_type(T, 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)
@@ -54,9 +51,16 @@ imag(o::Octonion) = [o.v1, o.v2, o.v3, o.v4, o.v5, o.v6, o.v7]
 conj(o::Octonion) = Octonion(o.s, -o.v1, -o.v2, -o.v3, -o.v4, -o.v5, -o.v6, -o.v7, o.norm)
 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)
 
+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)
+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))
+
 function normalize(o::Octonion)
   if (o.norm)
     return o
@@ -122,7 +126,7 @@ function exp(o::Octonion)
   s = o.s
   se = exp(s)
   scale = se
-  th = abs(Octonion(imag(o)))
+  th = abs_imag(o)
   if th > 0
     scale *= sin(th) / th
   end
@@ -140,7 +144,7 @@ end
 function log(o::Octonion)
   o, a = normalizea(o)
   s = o.s
-  M = abs(Octonion(imag(o)))
+  M = abs_imag(o)
   th = atan(M, s)
   if M > 0
     M = th / M
@@ -153,7 +157,7 @@ function log(o::Octonion)
                      o.v6 * M,
                      o.v7 * M)
   else
-    return Octonion(log(a), th, 0.0, 0.0)
+    return Octonion(complex(log(a), ifelse(iszero(a), zero(th), th)))
   end
 end
 

src/Quaternion.jl

@@ -10,26 +10,18 @@ const QuaternionF16 = Quaternion{Float16}
 const QuaternionF32 = Quaternion{Float32}
 const QuaternionF64 = Quaternion{Float64}
 
-(::Type{Quaternion{T}})(x::Real) where {T<:Real} = Quaternion(convert(T, x))
-(::Type{Quaternion{T}})(q::Quaternion{T}) where {T<:Real} = q
-(::Type{Quaternion{T}})(q::Quaternion) where {T<:Real} = Quaternion{T}(q.s, q.v1, q.v2, q.v3, q.norm)
+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(s::Real, a::AbstractVector) = Quaternion(s, a[1], a[2], a[3])
 Quaternion(a::AbstractVector) = Quaternion(0, a[1], a[2], a[3])
 
-convert(::Type{Quaternion{T}}, x::Real) where {T} = Quaternion(convert(T, x))
-convert(::Type{Quaternion{T}}, z::Complex) where {T} = Quaternion(convert(Complex{T}, z))
-convert(::Type{Quaternion{T}}, q::Quaternion{T}) where {T <: Real} = q
-convert(::Type{Quaternion{T}}, q::Quaternion) where {T} =
-    Quaternion(convert(T, q.s), convert(T, q.v1), convert(T, q.v2), convert(T, q.v3), q.norm)
-
-promote_rule(::Type{Quaternion{T}}, ::Type{T}) where {T <: Real} = Quaternion{T}
-promote_rule(::Type{Quaternion}, ::Type{T}) where {T <: Real} = Quaternion
 promote_rule(::Type{Quaternion{T}}, ::Type{S}) where {T <: Real, S <: Real} = Quaternion{promote_type(T, S)}
-promote_rule(::Type{Complex{T}}, ::Type{Quaternion{S}}) where {T <: Real, S <: Real} = Quaternion{promote_type(T, S)}
+promote_rule(::Type{Quaternion{T}}, ::Type{Complex{S}}) where {T <: Real, S <: Real} = Quaternion{promote_type(T, S)}
 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)
@@ -240,32 +232,19 @@ function linpol(p::Quaternion, q::Quaternion, t::Real)
         q = qm
     end
     c = p.s * q.s + p.v1 * q.v1 + p.v2 * q.v2 + p.v3 * q.v3
-    if c > - 1.0
-        if c < 1.0
-            o = acos(c)
-            s = sin(o)
-            sp = sin((1 - t) * o) / s
-            sq = sin(t * o) / s
-        else
-            sp = 1 - t
-            sq = t
-        end
-        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)
+    if c < 1.0
+        o = acos(c)
+        s = sin(o)
+        sp = sin((1 - t) * o) / s
+        sq = sin(t * o) / s
     else
-        s  =  p.v3
-        v1 = -p.v2
-        v2 =  p.v1
-        v3 = -p.s
-        sp = sin((0.5 - t) * pi)
-        sq = sin(t * pi)
-        Quaternion(s,
-                   sp * p.v1 + sq * v1,
-                   sp * p.v2 + sq * v2,
-                   sp * p.v3 + sq * v3, true)
+        sp = 1 - t
+        sq = t
     end
+    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)
 end
 
 quatrand(rng = Random.GLOBAL_RNG)  = quat(randn(rng), randn(rng), randn(rng), randn(rng))

src/Quaternions.jl

@@ -7,6 +7,7 @@ module Quaternions
   import Base: convert, promote_rule, float
   import Base: rand, randn
   import LinearAlgebra: lyap, norm, normalize, sylvester
+  using LinearAlgebra: cross, dot
   using Random
 
   Base.@irrational INV_SQRT_EIGHT 0.3535533905932737622004 sqrt(big(0.125))