From 3a50c776c164f751bf2462ca99cac9aa45b2dfa5 Mon Sep 17 00:00:00 2001
From: kimikage <kimikage.ceo@gmail.com>
Date: Sun, 9 Aug 2020 17:35:17 +0900
Subject: [PATCH] Improve `isapprox`

This changes the association between absolute tolerance and relative tolerance
from additive (`+`) to `max` to match the behavior in `Base` of Julia v1.
This also fixes an overflow problem and a rounding problem.
This speeds up the operation with default tolerance settings.
---
 src/FixedPointNumbers.jl | 25 ++++++++++++++++++++-----
 test/fixed.jl            |  9 +++++++++
 test/normed.jl           |  8 ++++++++
 3 files changed, 37 insertions(+), 5 deletions(-)

diff --git a/src/FixedPointNumbers.jl b/src/FixedPointNumbers.jl
index 86c17e4f..3755280a 100644
--- a/src/FixedPointNumbers.jl
+++ b/src/FixedPointNumbers.jl
@@ -93,16 +93,31 @@ end
 """
     isapprox(x::FixedPoint, y::FixedPoint; rtol=0, atol=max(eps(x), eps(y)))
 
-For FixedPoint numbers, the default criterion is that `x` and `y` differ by no more than `eps`, the separation between adjacent fixed-point numbers.
+For FixedPoint numbers, the default criterion is that `x` and `y` differ by no
+more than `eps`, the separation between adjacent fixed-point numbers.
 """
-function isapprox(x::T, y::T; rtol=0, atol=max(eps(x), eps(y))) where {T <: FixedPoint}
-    maxdiff = T(atol+rtol*max(abs(x), abs(y)))
-    rx, ry, rd = reinterpret(x), reinterpret(y), reinterpret(maxdiff)
-    abs(signed(widen1(rx))-signed(widen1(ry))) <= rd
+@inline function isapprox(x::X, y::X; rtol=0, atol=eps(X)) where {X <: FixedPoint}
+    n, m = minmax(x, y) # m >= n
+    if rtol == zero(rtol)
+        _isapprox_atol(m, n, atol)
+    elseif atol == zero(atol)
+        _isapprox_rtol(m, n, rtol)
+    else
+        _isapprox_atol(m, n, atol) | _isapprox_rtol(m, n, rtol)
+    end
 end
 function isapprox(x::FixedPoint, y::FixedPoint; rtol=0, atol=max(eps(x), eps(y)))
     isapprox(promote(x, y)...; rtol=rtol, atol=atol)
 end
+function _isapprox_atol(m::X, n::X, atol::X) where {X <: FixedPoint}
+    unsigned(m.i - n.i) <= unsigned(max(atol, zero(X)).i)
+end
+function _isapprox_atol(m::X, n::X, atol) where {X <: FixedPoint}
+    unsigned(m.i - n.i) <= div(atol, eps(X))
+end
+function _isapprox_rtol(m::X, n::X, rtol) where {X <: FixedPoint}
+    unsigned(m.i - n.i) <= rtol * max(unsigned(abs(m.i)), unsigned(abs(n.i)))
+end
 
 # predicates
 isinteger(x::FixedPoint) = x == trunc(x) # TODO: use floor(x) when dropping support for Fixed{Int8,8}
diff --git a/test/fixed.jl b/test/fixed.jl
index 8f8ad8a6..9de30222 100644
--- a/test/fixed.jl
+++ b/test/fixed.jl
@@ -406,6 +406,15 @@ end
         @test all(x -> x + eps(F) ≈ x, xs)
         @test !any(x -> x - eps(F) ≈ x + eps(F), xs)
     end
+
+    @test isapprox(-0.5Q0f7, -1Q0f7, rtol=0.5, atol=0) # issue 209
+    @test isapprox(typemin(Q0f7), typemax(Q0f7), rtol=2.0)
+    @test !isapprox(zero(Q0f7), typemax(Q0f7), rtol=0.9)
+    @test isapprox(zero(Q0f7), eps(Q0f7), rtol=1e-6) # atol = eps(Q0f7)
+    @test !isapprox(eps(Q0f7), zero(Q0f7), rtol=1e-6, atol=1e-6)
+    @test !isapprox(1.0Q6f1, 1.5Q6f1, rtol=0.3, atol=0) # 1.5 * 0.3 < eps(Q6f1)
+
+    @test isapprox(eps(Q8f7), eps(Q0f7), rtol=1e-6)
 end
 
 @testset "clamp" begin
diff --git a/test/normed.jl b/test/normed.jl
index f3a1ff20..d88a7a12 100644
--- a/test/normed.jl
+++ b/test/normed.jl
@@ -431,6 +431,14 @@ end
         @test all(x -> x + eps(N) ≈ x, xs)
         @test !any(x -> x - eps(N) ≈ x + eps(N), xs)
     end
+
+    @test isapprox(typemin(N0f8), typemax(N0f8), rtol=1.0)
+    @test !isapprox(zero(N0f8), typemax(N0f8), rtol=0.9)
+    @test isapprox(zero(N0f8), eps(N0f8), rtol=1e-6) # atol = eps(N0f8)
+    @test !isapprox(eps(N0f8), zero(N0f8), rtol=1e-6, atol=1e-6)
+    @test !isapprox(0.66N6f2, 1.0N6f2, rtol=0.3, atol=0) # 1.0 * 0.3 < eps(N6f2)
+
+    @test isapprox(eps(N8f8), eps(N0f8), rtol=1e-6)
 end
 
 @testset "comparison" begin