Add 13 elementary functions

src/GaussianRandomVariables.jl

@@ -6,6 +6,7 @@ using ThickNumbers
 import Base: +, -, *, /, //, ^, inv
 import Base: abs, abs2, max, min, sqrt
 import Base: log, log2, log10, exp, exp2, exp10, sin, cos, sincos
+import Base: atan, tan, sinh, cosh, tanh, asinh, acosh, atanh, asin, acos, cbrt, expm1, log1p
 
 export GVar, ±
 export skewness, moment_error, distrust
@@ -464,4 +465,53 @@ end
 
 sincos(a::GVar) = (sin(a), cos(a))
 
+atan(a::GVar{<:AbstractFloat}) =
+    isempty(a) ? a : gmap(atan, x -> 1/(1 + x^2), x -> -2x/(1 + x^2)^2,
+                          x -> (6x^2 - 2)/(1 + x^2)^3, x -> 24x*(1 - x^2)/(1 + x^2)^4, a)
+
+tan(a::GVar{<:AbstractFloat}) =
+    isempty(a) ? a : gmap(tan,
+                          x -> 1 + tan(x)^2,
+                          x -> (t = tan(x); 2t*(1 + t^2)),
+                          x -> (t = tan(x); (1 + t^2)*(6t^2 + 2)),
+                          x -> (t = tan(x); t*(1 + t^2)*(24t^2 + 16)), a)
+
+sinh(a::GVar{<:AbstractFloat}) = isempty(a) ? a : gmap(sinh, cosh, sinh, cosh, sinh, a)
+cosh(a::GVar{<:AbstractFloat}) = isempty(a) ? a : gmap(cosh, sinh, cosh, sinh, cosh, a)
+
+tanh(a::GVar{<:AbstractFloat}) =
+    isempty(a) ? a : gmap(tanh,
+                          x -> 1 - tanh(x)^2,
+                          x -> (t = tanh(x); -2t*(1 - t^2)),
+                          x -> (t = tanh(x); (1 - t^2)*(6t^2 - 2)),
+                          x -> (t = tanh(x); t*(1 - t^2)*(16 - 24t^2)), a)
+
+asinh(a::GVar{<:AbstractFloat}) =
+    isempty(a) ? a : gmap(asinh, x -> 1/sqrt(1 + x^2), x -> -x/sqrt((1 + x^2)^3),
+                          x -> (2x^2 - 1)/sqrt((1 + x^2)^5), x -> (9x - 6x^3)/sqrt((1 + x^2)^7), a)
+
+acosh(a::GVar{<:AbstractFloat}) =
+    isempty(a) ? a : gmap(acosh, x -> 1/sqrt(x^2 - 1), x -> -x/sqrt((x^2 - 1)^3),
+                          x -> (2x^2 + 1)/sqrt((x^2 - 1)^5), x -> -(9x + 6x^3)/sqrt((x^2 - 1)^7), a)
+
+atanh(a::GVar{<:AbstractFloat}) =
+    isempty(a) ? a : gmap(atanh, x -> 1/(1 - x^2), x -> 2x/(1 - x^2)^2,
+                          x -> (6x^2 + 2)/(1 - x^2)^3, x -> 24x*(1 + x^2)/(1 - x^2)^4, a)
+
+asin(a::GVar{<:AbstractFloat}) =
+    isempty(a) ? a : gmap(asin, x -> 1/sqrt(1 - x^2), x -> x/sqrt((1 - x^2)^3),
+                          x -> (1 + 2x^2)/sqrt((1 - x^2)^5), x -> (9x + 6x^3)/sqrt((1 - x^2)^7), a)
+
+acos(a::GVar{<:AbstractFloat}) =
+    isempty(a) ? a : gmap(acos, x -> -1/sqrt(1 - x^2), x -> -x/sqrt((1 - x^2)^3),
+                          x -> -(1 + 2x^2)/sqrt((1 - x^2)^5), x -> -(9x + 6x^3)/sqrt((1 - x^2)^7), a)
+
+cbrt(a::GVar{<:AbstractFloat}) =
+    isempty(a) ? a : gmap(cbrt, x -> 1/(3*cbrt(x)^2), x -> -2/(9*cbrt(x)^5),
+                          x -> 10/(27*cbrt(x)^8), x -> -80/(81*cbrt(x)^11), a)
+
+# Shifts by an exact constant, so these inherit `exp`'s and `log`'s treatment.
+expm1(a::GVar{<:AbstractFloat}) = exp(a) - 1
+log1p(a::GVar{<:AbstractFloat}) = log(a + 1)
+
 end # module

test/runtests.jl

@@ -118,6 +118,43 @@ ispositive(x) = x > 0
         @test testscalar(log2, 3.0, 0.2; filter=ispositive)
     end
 
+    @testset "elementary functions" begin
+        # For a κ3 = 0, err = 0 input, `gmap`'s outputs are exact functions of the
+        # first four derivatives, so ForwardDiff validates the hand-written
+        # derivative closures directly.
+        D(f) = x -> ForwardDiff.derivative(f, x)
+        for (f, c) in ((atan, 0.7), (atan, -1.2), (tan, 0.7), (sinh, 0.8), (cosh, 0.8),
+                       (tanh, 0.6), (tanh, -0.5), (asinh, 0.9), (acosh, 1.7),
+                       (atanh, 0.4), (asin, 0.3), (acos, 0.3), (cbrt, 1.3), (cbrt, -2.0))
+            σ = 0.1
+            g = f(GVar(c, σ))
+            f1, f2, f3, f4 = D(f)(c), D(D(f))(c), D(D(D(f)))(c), D(D(D(D(f))))(c)
+            σ2 = σ^2
+            @test g.center ≈ f(c) + f2*σ2/2  rtol=1e-8
+            @test g.σ ≈ sqrt(f1^2*σ2 + f2^2*σ2^2/2 + f1*f3*σ2^2)  rtol=1e-8
+            @test g.κ3 ≈ 3*f1^2*f2*σ2^2 + f2^3*σ2^3  rtol=1e-8
+            @test g.err ≈ abs(f4)*σ2^2/8  rtol=1e-8
+        end
+
+        # Sampled moments confirm the propagation end to end.
+        @test testscalar(atan, 0.7, 0.1)
+        @test testscalar(tanh, 0.6, 0.1)
+        @test testscalar(sinh, 0.8, 0.1)
+        @test testscalar(cbrt, 1.3, 0.1)
+        @test testscalar(asin, 0.3, 0.05)
+
+        # Constant shifts are exact, so these coincide with `exp` and `log`.
+        a = 0.4 ± 0.1
+        @test expm1(a) ⩪ exp(a) - 1
+        @test rad(expm1(a)) == rad(exp(a))
+        @test log1p(a) ⩪ log(a + 1)
+        @test rad(log1p(a)) == rad(log(a + 1))
+
+        # Type is preserved for narrower floats.
+        @test atan(GVar(1.0f0, 0.5f0)) isa GVar{Float32}
+        @test cbrt(GVar(2.0f0, 0.2f0)) isa GVar{Float32}
+    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