From 3ab59a0aa48f62cc95c4bb8bd308f59280bcf80f Mon Sep 17 00:00:00 2001
From: Seth Axen <seth@sethaxen.com>
Date: Wed, 7 Dec 2022 13:47:22 +0100
Subject: [PATCH] Fix log for negative quaternions (#115)

* Add log tests

* Add more general failing test

* Fix implementation of log for negative quat
---
 src/Quaternion.jl  | 17 ++++++++++++-----
 test/Quaternion.jl |  5 +++++
 2 files changed, 17 insertions(+), 5 deletions(-)

diff --git a/src/Quaternion.jl b/src/Quaternion.jl
index 9ec132e8..490a3b31 100644
--- a/src/Quaternion.jl
+++ b/src/Quaternion.jl
@@ -259,12 +259,19 @@ for f in (@static(VERSION ≥ v"1.6" ? (:sincos, :sincospi) : (:sincos,)))
     end
 end
 
+# this implementation is roughly 2x as fast as extend_analytic(log, q)
 function Base.log(q::Quaternion)
-    a = abs(q)
-    M = abs_imag(q)
-    theta = atan(M, q.s)
-    scale = theta / ifelse(iszero(M), oneunit(M), M)
-    return Quaternion(log(a), q.v1 * scale, q.v2 * scale, q.v3 * scale)
+    a = abs_imag(q)
+    theta = atan(a, q.s)
+    scale = theta / a
+    if a > 0
+        return Quaternion(log(abs(q)), scale * q.v1, scale * q.v2, scale * q.v3)
+    else
+        # q == real(q), so f(real(q)) may be real or complex.
+        # we choose to embed complex numbers in the quaternions by identifying the first
+        # imaginary quaternion basis with the complex imaginary basis.
+        return Quaternion(log(abs(q.s)), oftype(scale, theta), zero(scale), zero(scale))
+    end
 end
 
 Base.:^(q::Quaternion, w::Quaternion) = exp(w * log(q))
diff --git a/test/Quaternion.jl b/test/Quaternion.jl
index 51134e21..5895a433 100644
--- a/test/Quaternion.jl
+++ b/test/Quaternion.jl
@@ -385,6 +385,11 @@ end
         @testset "additional properties" begin
             @testset "log" begin
                 @test log(zero(QuaternionF64)) === Quaternion(-Inf, 0, 0, 0)
+                @test log(one(QuaternionF64)) === Quaternion(0.0, 0, 0, 0)
+                @test log(-one(QuaternionF64)) == Quaternion(0.0, π, 0, 0)
+                x = rand()
+                @test log(quat(x)) == quat(log(x))
+                @test log(quat(-x)) == Quaternion(reim(log(complex(-x)))..., 0, 0)
             end
 
             @testset "exp" begin