diff --git a/Project.toml b/Project.toml index 48ef108..2dd21cd 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "IsApprox" uuid = "28f27b66-4bd8-47e7-9110-e2746eb8bed7" authors = ["John Lapeyre "] -version = "0.1.2" +version = "0.1.1" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/src/IsApprox.jl b/src/IsApprox.jl index d64ec23..772e19c 100644 --- a/src/IsApprox.jl +++ b/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") diff --git a/src/core.jl b/src/core.jl index f492b30..db8a179 100644 --- a/src/core.jl +++ b/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 diff --git a/src/other_applications.jl b/src/other_applications.jl index e371133..82f947d 100644 --- a/src/other_applications.jl +++ b/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) """