Close ThickNumbers interface gaps

Project.toml

@@ -8,6 +8,8 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 ThickNumbers = "b57aa878-5b76-4266-befc-f8e007760995"
 
 [compat]
+Aqua = "0.8"
+ForwardDiff = "0.10, 1"
 HypothesisTests = "0.11"
 SpecialFunctions = "2"
 StableRNGs = "1"
@@ -17,10 +19,12 @@ ThickNumbers = "1.1.2"
 julia = "1.10"
 
 [extras]
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 HypothesisTests = "09f84164-cd44-5f33-b23f-e6b0d136a0d5"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["HypothesisTests", "StableRNGs", "Statistics", "Test"]
+test = ["Aqua", "ForwardDiff", "HypothesisTests", "StableRNGs", "Statistics", "Test"]

src/GaussianRandomVariables.jl

@@ -5,7 +5,7 @@ using ThickNumbers
 
 import Base: +, -, *, /, //, ^, inv
 import Base: abs, abs2, max, min, sqrt
-import Base: log, exp, sin, cos, sincos
+import Base: log, log2, log10, exp, exp2, exp10, sin, cos, sincos
 
 export GVar, ±
 export skewness, moment_error, distrust
@@ -132,6 +132,13 @@ end
 
 ThickNumbers.loval(a::GVar) = a.center - a.σ
 ThickNumbers.hival(a::GVar) = a.center + a.σ
+# `center` and `σ` are the stored primitives, so report the midpoint, width and
+# radius from them directly. Deriving them from `loval`/`hival` round-trips
+# through `center ± σ` and loses an ulp, which can make `rad` fall below the `σ`
+# handed to `midrad` and so violate the `midrad` containment contract.
+ThickNumbers.mid(a::GVar) = a.center
+ThickNumbers.wid(a::GVar) = 2 * a.σ
+ThickNumbers.rad(a::GVar) = a.σ
 function ThickNumbers.lohi(::Type{G}, lo, hi) where G<:GVar
     center = (lo + hi)/2
     # Measure the radius from the *rounded* center, so a narrow span far from zero
@@ -329,6 +336,9 @@ function ^(a::GVar{T}, p::Integer) where T<:AbstractFloat
 end
 ^(a::GVar{<:Real}, p::Integer) = float(a)^p
 
+# Resolves the ambiguity between `^(::GVar, ::Real)` and Base's `^(::Number, ::Rational)`.
+^(a::GVar, p::Rational) = a^float(p)
+
 function ^(a::GVar, p::Real)
     isinteger(p) && return a^Int(p)
     isempty(a) && return a
@@ -406,6 +416,14 @@ end
 log(a::GVar{<:AbstractFloat}) =
     isempty(a) ? a : gmap(log, x -> 1/x, x -> -1/x^2, x -> 2/x^3, x -> -6/x^4, a)
 
+# Reduce the other bases to the natural one by rescaling, which is exact: the
+# argument scaling of `exp2`/`exp10` and the result scaling of `log2`/`log10`
+# both go through the constant-factor rules without adding truncation error.
+exp2(a::GVar{T}) where T<:AbstractFloat = exp(a * log(T(2)))
+exp10(a::GVar{T}) where T<:AbstractFloat = exp(a * log(T(10)))
+log2(a::GVar{T}) where T<:AbstractFloat = log(a) / log(T(2))
+log10(a::GVar{T}) where T<:AbstractFloat = log(a) / log(T(10))
+
 sqrt(a::GVar{<:AbstractFloat}) =
     isempty(a) ? a : gmap(sqrt, x -> 1/(2*sqrt(x)), x -> -1/(4*sqrt(x^3)),
                           x -> 3/(8*sqrt(x^5)), x -> -15/(16*sqrt(x^7)), a)

test/runtests.jl

@@ -3,7 +3,11 @@ using ThickNumbers
 using Statistics
 using HypothesisTests
 using StableRNGs
+using Aqua
+using ForwardDiff
 using Test
+push!(LOAD_PATH, pkgdir(ThickNumbers, "ThickNumbersInterfaceTests"))
+using ThickNumbersInterfaceTests
 
 # These comparisons are Monte Carlo, so an unseeded generator turns the suite into a
 # coin flip. Julia's own streams are not reproducible across releases, so seed a
@@ -37,6 +41,19 @@ end
 ispositive(x) = x > 0
 
 @testset "GaussianRandomVariables.jl" begin
+    @testset "Aqua" begin
+        Aqua.test_all(GaussianRandomVariables)
+    end
+
+    @testset "ThickNumbers interface" begin
+        ThickNumbersInterfaceTests.test_reserved()
+        ThickNumbersInterfaceTests.test_required(GVar{Float64})
+        ThickNumbersInterfaceTests.test_required(GVar, [Float32, Float64])
+        ThickNumbersInterfaceTests.test_optional(GVar{Float64})
+        ThickNumbersInterfaceTests.test_optional(GVar, [Float32, Float64])
+        ThickNumbersInterfaceTests.test_FPTNviolations(GVar(1.0, 2.0))
+    end
+
     @testset "arithmetic" begin
         x = 3 ± 1
         e = GVar(0, -1)          # empty
@@ -59,6 +76,10 @@ ispositive(x) = x > 0
         @test x^0 ⩪ oneunit(x)
         @test x^1 ⩪ x
         @test 1/(1/x) ⩪ x rtol=0.1
+        @test x^(1//2) ⩪ x^0.5      # a Rational exponent reduces to the real one
+        @test x^(4//2) ⩪ x^2        # an integer-valued Rational takes the exact path
+        # `mid`, `wid` and `rad` report the stored primitives exactly.
+        @test mid(x) == 3.0 && wid(x) == 2.0 && rad(x) == 1.0
     end
 
     @testset "distributions" begin
@@ -76,6 +97,27 @@ ispositive(x) = x > 0
         end
     end
 
+    @testset "log and exp bases" begin
+        a = 3.0 ± 0.2
+        # `exp2`/`exp10`/`log2`/`log10` are rescalings of the natural-base
+        # functions, and the rescaling is exact, so the spans coincide.
+        @test exp2(a) ⩪ exp(a * log(2.0))
+        @test exp10(a) ⩪ exp(a * log(10.0))
+        @test log2(a) ⩪ log(a) / log(2.0)
+        @test log10(a) ⩪ log(a) / log(10.0)
+        # Zero-spread inputs reduce to the ordinary functions.
+        @test exp2(GVar(3.0, 0.0)) ⩪ GVar(8.0, 0.0)
+        @test exp10(GVar(2.0, 0.0)) ⩪ GVar(100.0, 0.0)
+        @test log2(GVar(8.0, 0.0)) ⩪ GVar(3.0, 0.0)
+        @test log10(GVar(100.0, 0.0)) ⩪ GVar(2.0, 0.0)
+        # Type is preserved for narrower floats.
+        @test exp2(GVar(1.0f0, 0.5f0)) isa GVar{Float32}
+        @test log2(GVar(3.0f0, 0.2f0)) isa GVar{Float32}
+        # Sampled moments confirm the propagation.
+        @test testscalar(exp2, 1.0, 0.2; rtol=0.03, n=10^6)
+        @test testscalar(log2, 3.0, 0.2; filter=ispositive)
+    end
+
     # x^p is a polynomial, so the Gaussian moments terminate: mean, variance and
     # third cumulant are all exact.
     @testset "exact moments of integer powers" begin
@@ -344,4 +386,24 @@ ispositive(x) = x > 0
         @test valuetype(GVar(1.0, 2.0)) === Float64
         @test convert(GVar{Float64}, 3) ⩪ GVar(3.0, 0.0)
     end
+
+    # Restricting the scalar-times-`GVar` methods to `BaseReals` keeps a
+    # `ForwardDiff.Dual` from matching them, so differentiation carries `GVar`
+    # values through in the `Dual` slots. Each analytic derivative is written as
+    # the same sequence of `GVar` operations ForwardDiff performs, so the spans
+    # match exactly (note `w * w`, not `w^2`: ForwardDiff differentiates one order
+    # at a time, and `GVar` treats repeated factors as independent).
+    @testset "ForwardDiff" begin
+        x, w = 2.0, 0.0 ± 0.5
+        ϕ(t) = sin(x + (1 + t) * w)
+        ϕ′(t) = cos(x + (1 + t) * w) * w
+        ϕ′′(t) = -sin(x + (1 + t) * w) * w * w
+        ϕ′′′(t) = -cos(x + (1 + t) * w) * w * w * w
+        dϕ(t) = ForwardDiff.derivative(ϕ, t)
+        ddϕ(t) = ForwardDiff.derivative(dϕ, t)
+        dddϕ(t) = ForwardDiff.derivative(ddϕ, t)
+        @test ϕ′(0) ≐ dϕ(0)
+        @test ϕ′′(0) ≐ ddϕ(0)
+        @test ϕ′′′(0) ≐ dddϕ(0)
+    end
 end