Add UpToPhase

Project.toml

@@ -1,7 +1,7 @@
 name = "IsApprox"
 uuid = "28f27b66-4bd8-47e7-9110-e2746eb8bed7"
 authors = ["John Lapeyre <[email protected]>"]
-version = "0.1.2"
+version = "0.1.1"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/IsApprox.jl

@@ -1,6 +1,6 @@
 module IsApprox
 
-export AbstractApprox, Equal, EachApprox, Approx
+export AbstractApprox, Equal, EachApprox, Approx, UpToPhase
 export ispossemidef, isunitary, isinvolution, isidempotent, isnormal, commutes, anticommutes
 
 include("core.jl")

src/core.jl

@@ -71,3 +71,52 @@ function Base.isapprox(a::EachApprox, A::AbstractArray, B::AbstractArray)
     end
     return true
 end
+
+"""
+    UpToPhase <: AbstractApprox
+
+Demands that each pair of elements are approximately equal up to a phase,
+that is a number whose absolute value is one.
+`kw` are keyword pairs that are forwarded to `isapprox`.
+"""
+struct UpToPhase{T} <: AbstractApprox
+    kw::T
+end
+UpToPhase(; kws...) = UpToPhase(kws)
+
+function Base.isapprox(a::UpToPhase, x::Number, y::Number)
+    aa = Approx(;a.kw...)
+    if isapprox(aa, x, zero(x))
+        return isapprox(aa, y, zero(y))
+    elseif isapprox(aa, y, zero(y))
+        return isapprox(aa, x, zero(x))
+    end
+    return isunitary(x / y, aa)
+end
+
+function Base.isapprox(_app::UpToPhase,  A::AbstractArray, B::AbstractArray)
+    n1 = length(A)
+    n2 = length(B)
+    if n1 != n2
+        return false
+    end
+    app = Approx(;_app.kw...)
+    seen_non_zero_flag = false
+    z = zero(eltype(A)) # TODO use promotion
+    for (a, b) in zip(A, B)
+        if iszero(a, app)
+            !iszero(b, app) && return false
+        elseif iszero(b, app)
+            !iszero(a, app) && return false
+        else
+            if ! seen_non_zero_flag
+                z = a/b
+                isunitary(z, app) || return false
+                seen_non_zero_flag = true
+            else
+                isapprox(app, a/b, z) || return false
+            end
+        end
+    end
+    return true
+end

src/other_applications.jl

@@ -99,7 +99,9 @@ Return `true` if `m` is unitary. If `m` is real, this tests orthogonality.
 isunitary(m::AbstractMatrix, approx_test::AbstractApprox=Equal()) =
     _isunitary(m, approx_test, LinearAlgebra.dot, _identity)
 
-isunitary(x::Number, approx_test::AbstractApprox=Equal()) = isone(abs(x), approx_test)
+# Careful, because we use abs2, the requested precision will probably be wrong. But, abs2 is much faster.
+# Maybe there is way to pass this information.
+isunitary(x::Number, approx_test::AbstractApprox=Equal()) = isone(abs2(x), approx_test)
 isunitary(J::LinearAlgebra.UniformScaling, approx_test::AbstractApprox=Equal()) = isunitary(J.λ, approx_test)
 
 """