Implement sylvester and lyap for quaternions (#65)

* Implement iszero for Quaternion

* Implement sylvester and lyap

* Simplify implementation

* Hit other branch

* Test for more random numbers

* Define isnan

* Match real/complex behavior for returning NaNs

* Apply suggestions from code review

Co-authored-by: Yuto Horikawa <hyrodium@gmail.com>

* Add promotion methods

* Add methods for combinations of real arguments

* Implement isreal

* Simplify when one is real

* Use bitwise operations

* Add isinf

* Test isfinite/isinf

* Test sylvester/lyap more thoroughly

* Avoid huge fractions

* Skip test when lyap did not promote

* Call sylvester

* Move sylvester/lyap tests to new test file

* Increment version number

Co-authored-by: Yuto Horikawa <hyrodium@gmail.com>

Author: sethaxen

Project.toml

@@ -1,6 +1,6 @@
 name = "Quaternions"
 uuid = "94ee1d12-ae83-5a48-8b1c-48b8ff168ae0"
-version = "0.5.0"
+version = "0.5.1"
 
 [deps]
 DualNumbers = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74"

src/Quaternion.jl

@@ -55,7 +55,11 @@ 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)
 
-isfinite(q::Quaternion) = q.norm ? true : (isfinite(q.s) && isfinite(q.v1) && isfinite(q.v2) && isfinite(q.v3))
+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)
+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))
 
 function normalize(q::Quaternion)
     if (q.norm)
@@ -381,3 +385,38 @@ function slerp(qa::Quaternion{T}, qb::Quaternion{T}, t::T) where {T}
         qa.v3 * ratio_a + qm.v3 * ratio_b,
     )
 end
+
+function sylvester(a::Quaternion{T}, b::Quaternion{T}, c::Quaternion{T}) where {T<:Real}
+    isreal(a) && return sylvester(real(a), b, c)
+    isreal(b) && return sylvester(a, real(b), c)
+    abs2a = abs2(a)
+    abs2b = abs2(b)
+    if abs2a > abs2b
+        inva = conj(a) / abs2a
+        d1 = -2real(b) - a - inva * abs2b
+        x = d1 \ (c + inva * c * conj(b))
+    else
+        invb = conj(b) / abs2b
+        d2 = -2real(a) - b - invb * abs2a
+        x = (c + conj(a) * c * invb) / d2
+    end
+    return x
+end
+sylvester(a::Quaternion, b::Quaternion, c::Quaternion) = sylvester(promote(a, b, c)...)
+sylvester(a::Quaternion, b::Quaternion, c::Real) = sylvester(promote(a, b, c)...)
+# if either a or b commute with x, use a simpler expression
+sylvester(a::Real, b::Real, c::Quaternion) = c / -(a + b)
+sylvester(a::Real, b::Quaternion, c::Quaternion) = c / -(a + b)
+sylvester(a::Quaternion, b::Real, c::Quaternion) = -(a + b) \ c
+sylvester(a::Real, b::Quaternion, c::Real) = -c / (a + b)
+sylvester(a::Quaternion, b::Real, c::Real) = (a + b) \ -c
+
+function lyap(a::Quaternion{T}, c::Quaternion{T}) where {T<:Real}
+    # if a commutes with c, use a simpler expression
+    (isreal(a) || isreal(c)) && return c / -2real(a)
+    return (c + a \ c * a) / -4real(a)
+end
+lyap(a::Quaternion, c::Quaternion) = lyap(promote(a, c)...)
+# if a commutes with c, use a simpler expression
+lyap(a::Real, c::Quaternion) = c / -2a
+lyap(a::Quaternion, c::Real) = c / -2real(a)

src/Quaternions.jl

@@ -3,10 +3,10 @@ __precompile__()
 module Quaternions
 
   import Base: +, -, *, /, ^, ==
-  import Base: abs, abs2, angle, conj, cos, exp, inv, isfinite, log, real, sin, sqrt
+  import Base: abs, abs2, angle, conj, cos, exp, inv, isreal, isfinite, isinf, iszero, isnan, log, real, sin, sqrt
   import Base: convert, promote_rule, float
   import Base: rand, randn
-  import LinearAlgebra: norm, normalize
+  import LinearAlgebra: lyap, norm, normalize, sylvester
   using Random
 
   Base.@irrational INV_SQRT_EIGHT 0.3535533905932737622004 sqrt(big(0.125))

test/Quaternion.jl

@@ -439,4 +439,68 @@ Base.:(/)(a::MyReal, b::Real) = a.val / b
             end
         end
     end
+
+    @testset "sylvester/lyap" begin
+        Ts = (Float64, QuaternionF64)
+        Ttrips = [(Ta, Tb, Tc) for Ta in Ts for Tb in Ts for Tc in Ts]
+        Ttrips = filter(x -> any(y -> y <: Quaternion, x), Ttrips)
+        @testset "($Ta, $Tb, $Tc)" for (Ta, Tb, Tc) in Ttrips 
+            for _ in 1:100
+                a = randn(Ta)
+                b = randn(Tb)
+                c = randn(Tc)
+                x = @inferred sylvester(a, b, c)
+                @test a * x + x * b ≈ -c
+                x = @inferred sylvester(b, a, c)
+                @test b * x + x * a ≈ -c
+                @test iszero(sylvester(a, b, zero(c)))
+                @test sylvester(a, zero(b), c) ≈ a \ -c
+                @test sylvester(zero(a), b, c) ≈ -c / b
+                @test iszero(sylvester(zero(a), b, zero(c)))
+                @test iszero(sylvester(a, zero(b), zero(c)))
+                @test iszero(sylvester(a, b, zero(c)))
+                # @test isnan(sylvester(zero(a), zero(b), c))
+    
+                @test @inferred(lyap(a, c)) ≈ sylvester(a, a', c)
+                @test @inferred(lyap(b, c)) ≈ sylvester(b, b', c)
+                @test iszero(lyap(a, zero(c)))
+            end
+            @testset "nan/inf return same as for complex" begin
+                Tza, Tzb, Tzc = map(T -> T <: Quaternion ? complex(real(T)) : T, (Ta, Tb, Tc))
+                a, b = zero(Ta), zero(Tb)
+                za, zb = zero(Tza), zero(Tzb)
+                @testset for f in (one, zero, randn)
+                    x = sylvester(a, b, f(Tc))
+                    zx = sylvester(za, zb, f(Tzc))
+                    if isinf(zx)
+                        @test isinf(x)
+                    elseif isnan(zx)
+                        @test isnan(x)
+                    end
+                    if VERSION ≥ v"1.7"
+                        x = lyap(a, f(Tc))
+                        zx = lyap(za, f(Tzc))
+                        if isinf(zx)
+                            @test isinf(x)
+                        elseif isnan(zx)
+                            @test isnan(x)
+                        end
+                    end
+                end
+            end    
+        end
+        @testset "rational" begin
+            a = Quaternion(1, 2, 3, 4)
+            b = Quaternion(1//2, 2//2, 3//2, 4//2)
+            c = Quaternion(-1//2, 2//2, -4//2, -3//2)
+            @test @inferred(sylvester(a, b, c)) isa Quaternion{Rational{Int}}
+            @test @inferred(lyap(a, c)) isa Quaternion{Rational{Int}}
+            a = randn(QuaternionF32)
+            @test @inferred(sylvester(a, b, c)) isa QuaternionF32
+            @test @inferred(lyap(a, c)) isa QuaternionF32
+            null = zero(Quaternion{Rational{Int}})
+            @test_throws DivideError sylvester(null, null, null)
+            @test_throws DivideError lyap(null, null)
+        end
+    end    
 end