Add support for RealDot (#119)

* add support for RealDot

* add tests for realdot

* add compat table for RealDot

* Add `@inferred` in test for `realdot`

Co-authored-by: Seth Axen <seth@sethaxen.com>

* use realdot in slerp

* replace ≈ with == in realdot test

* remove muladd in realdot definition

* Add suggested comment

Co-authored-by: Seth Axen <seth@sethaxen.com>

* Use explicit RealDot.realdot

Co-authored-by: Seth Axen <seth@sethaxen.com>

* use RealDot.realdot in abs2

* add commutative law test for realdot

* add commutative law test for real∘dot

* bump version to v0.7.3

Co-authored-by: Seth Axen <seth@sethaxen.com>

Project.toml

@@ -1,13 +1,15 @@
 name = "Quaternions"
 uuid = "94ee1d12-ae83-5a48-8b1c-48b8ff168ae0"
-version = "0.7.2"
+version = "0.7.3"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+RealDot = "c1ae055f-0cd5-4b69-90a6-9a35b1a98df9"
 
 [compat]
 Aqua = "0.5"
+RealDot = "0.1"
 julia = "1"
 
 [extras]

src/Quaternion.jl

@@ -139,7 +139,7 @@ Base.conj(q::Quaternion) = Quaternion(q.s, -q.v1, -q.v2, -q.v3)
 Base.abs(q::Quaternion) = sqrt(abs2(q))
 Base.float(q::Quaternion{T}) where T = convert(Quaternion{float(T)}, q)
 abs_imag(q::Quaternion) = sqrt(q.v2 * q.v2 + (q.v1 * q.v1 + q.v3 * q.v3)) # ordered to match abs2
-Base.abs2(q::Quaternion) = (q.s * q.s + q.v2 * q.v2) + (q.v1 * q.v1 + q.v3 * q.v3)
+Base.abs2(q::Quaternion) = RealDot.realdot(q,q)
 Base.inv(q::Quaternion) = conj(q) / abs2(q)
 
 Base.isreal(q::Quaternion) = iszero(q.v1) & iszero(q.v2) & iszero(q.v3)
@@ -320,7 +320,7 @@ true
     iszero(qb0) && throw(DomainError(qb0, "The input quaternion must be non-zero."))
     qa = qa0 / abs(qa0)
     qb = qb0 / abs(qb0)
-    coshalftheta = qa.s * qb.s + qa.v1 * qb.v1 + qa.v2 * qb.v2 + qa.v3 * qb.v3
+    coshalftheta = RealDot.realdot(qa, qb)
 
     if coshalftheta < 0
         qb = -qb
@@ -386,6 +386,9 @@ LinearAlgebra.lyap(a::Quaternion, c::Quaternion) = lyap(promote(a, c)...)
 LinearAlgebra.lyap(a::Real, c::Quaternion) = c / -2a
 LinearAlgebra.lyap(a::Quaternion, c::Real) = c / -2real(a)
 
+## RealDot
+# ordering chosen so that real(p'q) == real(q'p) == realdot(p, q) == realdot(q, p), i.e. exact equality
+@inline RealDot.realdot(p::Quaternion, q::Quaternion) = (p.s * q.s + p.v2 * q.v2) + (p.v1 * q.v1 + p.v3 * q.v3)
 Base.widen(::Type{Quaternion{T}}) where {T} = Quaternion{widen(T)}
 
 Base.flipsign(x::Quaternion, y::Real) = ifelse(signbit(y), -x, x)

src/Quaternions.jl

@@ -2,6 +2,7 @@ module Quaternions
 
 using Random
 using LinearAlgebra
+using RealDot
 
 include("Quaternion.jl")
 

test/Quaternion.jl

@@ -580,6 +580,21 @@ end
         end
     end
 
+    @testset "RealDot with $T" for T in (Float32, Float64)
+        for _ in 1:10
+            q1 = randn(Quaternion{T})
+            q2 = randn(Quaternion{T})
+            # Check real∘dot is equal to realdot.
+            @test real(dot(q1,q2)) == @inferred(realdot(q1,q2))
+            # Check realdot is commutative.
+            @test realdot(q1,q2) == realdot(q2,q1)
+            # Check real∘dot is also commutative just in case.
+            @test real(dot(q1,q2)) == real(dot(q2,q1))
+            # Check the return type of realdot is correct.
+            @test realdot(q1,q2) isa T
+        end
+    end
+
     @testset "widen" begin
         @test widen(Quaternion{Int}) === Quaternion{Int128}
         @test widen(QuaternionF32) === QuaternionF64

test/runtests.jl

@@ -1,6 +1,7 @@
 using Test
 using Quaternions
 using Aqua
+using RealDot
 
 Aqua.test_all(Quaternions)