More accurate distance for projective spaces (#584)

* Improve accuracy for distance on projective spaces

* Add tests for small distances

* Increment patch number

* Don't test BigFloat

* Update test/manifolds/projective_space.jl

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

---------

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

Project.toml

@@ -1,7 +1,7 @@
 name = "Manifolds"
 uuid = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
 authors = ["Seth Axen <[email protected]>", "Mateusz Baran <[email protected]>", "Ronny Bergmann <[email protected]>", "Antoine Levitt <[email protected]>"]
-version = "0.8.53"
+version = "0.8.54"
 
 [deps]
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"

src/manifolds/ProjectiveSpace.jl

@@ -147,7 +147,16 @@ Note that this definition is similar to that of the [`AbstractSphere`](@ref).
 However, the absolute value ensures that all equivalent `p` and `q` have the same pairwise
 distance.
 """
-distance(::AbstractProjectiveSpace, p, q) = acos(min(abs(dot(p, q)), 1))
+function distance(::AbstractProjectiveSpace, p, q)
+    z = dot(p, q)
+    cosθ = abs(z)
+    T = float(real(Base.promote_eltype(p, q)))
+    # abs and relative error of acos is less than sqrt(eps(T))
+    cosθ < 1 - sqrt(eps(T)) / 8 && return acos(cosθ)
+    # improved accuracy for q close to p or -p
+    λ = sign(z)
+    return 2 * abs(atan(norm(p .* λ .- q), norm(p .* λ .+ q)))
+end
 
 function exp!(M::AbstractProjectiveSpace, q, p, X)
     θ = norm(M, p, X)

test/manifolds/projective_space.jl

@@ -305,4 +305,23 @@ include("../utils.jl")
         @test typeof(get_embedding(M)) === Euclidean{Tuple{2,2},ℂ}
         @test representation_size(M) == (2, 2)
     end
+
+    @testset "small distance tests" begin
+        @testset for fT in (Float32, Float64), T in (fT, Complex{fT}, Quaternion{fT})
+            𝔽 = T isa Complex ? ℂ : (T isa Quaternion ? ℍ : ℝ)
+            M = ProjectiveSpace(2, 𝔽)
+            rT = real(T)
+            atol = rtol = sqrt(eps(rT))
+            @testset for t in (zero(rT), eps(rT)^(1 // 4) / 8, eps(rT)^(1 // 4)),
+                λ in (one(T), (T <: Real ? -one(T) : sign(randn(T))))
+
+                p = project(M, randn(T, representation_size(M)))
+                X = project(M, p, randn(T, representation_size(M)))
+                X ./= norm(M, p, X)
+                project!(M, X, p, X)
+                @test distance(M, p, exp(M, p, t * X) * λ) ≈ t atol = atol rtol = rtol
+                @test distance(M, p, exp(M, p, (π - t) * X) * λ) ≈ t atol = atol rtol = rtol
+            end
+        end
+    end
 end