From 3910eacfd7d3ef876dfb63b21c874dbc4935e165 Mon Sep 17 00:00:00 2001
From: Tim Holy <tim.holy@gmail.com>
Date: Fri, 18 Apr 2025 11:17:51 -0500
Subject: [PATCH] exp: return the exact answer

Also makes the nonlinear tests more demanding.
---
 Project.toml                   |  4 +++-
 src/GaussianRandomVariables.jl |  3 +--
 test/runtests.jl               | 27 +++++++++++++++++----------
 3 files changed, 21 insertions(+), 13 deletions(-)

diff --git a/Project.toml b/Project.toml
index 6d9b061..5453a29 100644
--- a/Project.toml
+++ b/Project.toml
@@ -9,12 +9,14 @@ ThickNumbers = "b57aa878-5b76-4266-befc-f8e007760995"
 
 [compat]
 SpecialFunctions = "2"
+Statistics = "1"
 ThickNumbers = "1"
 julia = "1"
 
 [extras]
 HypothesisTests = "09f84164-cd44-5f33-b23f-e6b0d136a0d5"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["HypothesisTests", "Test"]
+test = ["HypothesisTests", "Statistics", "Test"]
diff --git a/src/GaussianRandomVariables.jl b/src/GaussianRandomVariables.jl
index ae24ef7..e075238 100644
--- a/src/GaussianRandomVariables.jl
+++ b/src/GaussianRandomVariables.jl
@@ -213,9 +213,8 @@ max(a::GVar, b::Real) = max(b, a)
 
 function Base.exp(a::GVar{<:AbstractFloat})
     isempty(a) && return a
-    return GVar(meanvar(exp, exp, exp, a.center, a.σ)...)
+    return GVar(exp(a.center + a.σ^2/2), exp(a.center + a.σ^2)*sqrt(1 - exp(-a.σ^2)))
 end
-Base.exp(a::GVar) = exp(float(a))
 
 function Base.log(a::GVar{<:AbstractFloat})
     isempty(a) && return a
diff --git a/test/runtests.jl b/test/runtests.jl
index 089f354..cac8b93 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,17 +1,24 @@
 using GaussianRandomVariables
 using ThickNumbers
+using Statistics
 using HypothesisTests
 using Test
 
-function testscalar(f, μ, σ; n=1000, pval = 0.001)
-    x = μ .+ σ .* randn(n)
+function testscalar(f, μ, σ; n=1000, pval = 0.001, filter=Returns(true), rtol=nothing)
+    x = Base.filter(filter, μ .+ σ .* randn(n))
+    @assert length(x) >= 0.9 * n
     y = f.(x)
     g = GVar(μ, σ)
     fg = f(g)
-    z = mid(fg) .+ rad(fg) .* randn(n)
-    return pvalue(EqualVarianceTTest(y, z)) > pval && pvalue(LeveneTest(y, z)) > pval
+    if rtol === nothing
+        z = mid(fg) .+ rad(fg) .* randn(n)
+        return pvalue(EqualVarianceTTest(y, z)) > pval && pvalue(LeveneTest(y, z)) > pval
+    end
+    return isapprox(mean(y), mid(fg); rtol) && isapprox(std(y), rad(fg); rtol)
 end
 
+ispositive(x) = x > 0
+
 @testset "GaussianRandomVariables.jl" begin
     x = 3 ± 1
     e = GVar(0, -1)          # empty
@@ -27,14 +34,14 @@ end
     @test rad(x - x) ≈ sqrt(2)
     @test mid(1/x) ≈ 1/mid(x) + rad(x)^2/mid(x)^3
     @test rad(1/x)^2 ≈ rad(x)^2 / mid(x)^4 + 2 * rad(x)^4 / mid(x)^6
-    @test sqrt(x) ⩪ exp(0.5 * log(x)) rtol=1e-3
+    @test sqrt(x) ⩪ exp(0.5 * log(x)) rtol=1e-2
     @test sqrt(x) ⩪ x^0.5
 
-    for μ in (2, 3, 10), σ in (0.25,)
+    for μ in (2, 3, 10), σ in (0.1, 0.25, 1)
         @test testscalar(x -> x^2, μ, σ)
-        @test testscalar(x -> x^1.8, μ, σ)
-        @test testscalar(exp, μ, σ)
-        @test testscalar(log, μ, σ)
-        @test testscalar(sqrt, μ, σ)
+        @test testscalar(x -> x^1.8, μ, σ; filter=ispositive)
+        @test testscalar(exp, μ, σ; rtol=0.03, n=10^6)
+        μ - 3*σ >= 0 &&  @test testscalar(log, μ, σ; filter=ispositive)
+        @test testscalar(sqrt, μ, σ; filter=ispositive)
     end
 end