Implement sin and cos

src/GaussianRandomVariables.jl

@@ -5,7 +5,7 @@ using ThickNumbers
 
 import Base: +, -, *, /, //, ^, inv
 import Base: abs, abs2, max, min, sqrt
-import Base: log, exp
+import Base: log, exp, sin, cos, sincos
 
 export GVar, ±
 export skewness, moment_error, distrust
@@ -410,4 +410,40 @@ 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)
 
+# sin and cos of a Gaussian are not Gaussian, but their first three moments are
+# exact in closed form from the characteristic function E[exp(itx)] = exp(itc - t²σ²/2).
+# With s = sin c, k = cos c and the Gaussian damping factors d1 = exp(-σ²/2),
+# d2 = exp(-2σ²), d9 = exp(-9σ²/2):
+#
+#   E[sin x] = s·d1                       E[cos x] = k·d1
+#   E[sin²x] = (1 - cos2c·d2)/2           E[cos²x] = (1 + cos2c·d2)/2
+#   E[sin³x] = (3s·d1 - sin3c·d9)/4       E[cos³x] = (3k·d1 + cos3c·d9)/4
+#
+# The mean is exact in σ, so nothing is charged to `err`; as for `exp`, only the
+# input's own skew and center error propagate, to leading order through the local
+# derivatives (f' = k, f'' = -s, f''' = -k for sin; negated/rotated for cos).
+function sin(a::GVar{<:AbstractFloat})
+    isempty(a) && return a
+    c, σ, κ3, err = a.center, a.σ, a.κ3, a.err
+    d1, d2, d9 = exp(-σ^2/2), exp(-2σ^2), exp(-9σ^2/2)
+    s, k = sin(c), cos(c)
+    μ = s*d1
+    v = (1 - cos(2c)*d2)/2 - μ^2
+    e3 = (3s*d1 - sin(3c)*d9)/4
+    return assemble(μ - k*κ3/6, v - s*k*κ3, (e3 - 3μ*(v + μ^2) + 2μ^3) + k^3*κ3, abs(k*d1)*err)
+end
+
+function cos(a::GVar{<:AbstractFloat})
+    isempty(a) && return a
+    c, σ, κ3, err = a.center, a.σ, a.κ3, a.err
+    d1, d2, d9 = exp(-σ^2/2), exp(-2σ^2), exp(-9σ^2/2)
+    s, k = sin(c), cos(c)
+    μ = k*d1
+    v = (1 + cos(2c)*d2)/2 - μ^2
+    e3 = (3k*d1 + cos(3c)*d9)/4
+    return assemble(μ + s*κ3/6, v + s*k*κ3, (e3 - 3μ*(v + μ^2) + 2μ^3) - s^3*κ3, abs(s*d1)*err)
+end
+
+sincos(a::GVar) = (sin(a), cos(a))
+
 end # module

test/runtests.jl

@@ -193,6 +193,52 @@ ispositive(x) = x > 0
         @test abs2(3 ± 0.5) ⩪ (3 ± 0.5)^2
     end
 
+    # sin and cos of a Gaussian are not Gaussian, but their first three moments are
+    # exact, so mid, rad and κ3 match closed forms and nothing is charged to `err`.
+    @testset "sin and cos" begin
+        c, σ = 0.7, 0.5
+        d1, d2 = exp(-σ^2/2), exp(-2σ^2)
+        s = sin(GVar(c, σ))
+        @test mid(s) ≈ sin(c)*d1
+        @test rad(s)^2 ≈ (1 - cos(2c)*d2)/2 - (sin(c)*d1)^2
+        @test moment_error(s) == 0
+        k = cos(GVar(c, σ))
+        @test mid(k) ≈ cos(c)*d1
+        @test rad(k)^2 ≈ (1 + cos(2c)*d2)/2 - (cos(c)*d1)^2
+        @test moment_error(k) == 0
+
+        # A zero-mean input is symmetric: sin has zero mean and no skew, while cos
+        # peaks at 1 and folds its mass downward, so it is skewed toward low values.
+        @test mid(sin(0 ± σ)) == 0
+        @test sin(0 ± σ).κ3 == 0
+        @test mid(cos(0 ± σ)) ≈ exp(-σ^2/2)
+        @test skewness(cos(0 ± σ)) < 0
+
+        # A vanishing spread is an ordinary number.
+        @test sin(GVar(c, 0.0)) ⩪ GVar(sin(c), 0.0)
+        @test cos(GVar(c, 0.0)) ⩪ GVar(cos(c), 0.0)
+
+        # `sincos` returns the pair; the empty set propagates through all three.
+        sc = sincos(GVar(c, σ))
+        @test sc[1] ⩪ sin(GVar(c, σ)) && sc[2] ⩪ cos(GVar(c, σ))
+        e = GVar(0.0, -1.0)
+        @test isempty(sin(e)) && isempty(cos(e))
+        @test all(isempty, sincos(e))
+
+        @test @inferred(sin(GVar(1.0, 0.5))) isa GVar{Float64}
+        @test @inferred(sincos(GVar(1.0, 0.5))) isa Tuple{GVar{Float64},GVar{Float64}}
+
+        # Sampled moments confirm the closed forms across a range of spreads.
+        for (μ, σ) in ((0.6, 0.2), (1.0, 0.5), (2.0, 0.8), (0.9, 1.0))
+            @test testscalar(sin, μ, σ; rtol=0.02, n=10^6)
+            @test testscalar(cos, μ, σ; rtol=0.02, n=10^6)
+        end
+        for (μ, σ) in ((1.0, 0.5), (2.0, 0.8), (0.9, 1.0))
+            @test testskew(sin, μ, σ)
+            @test testskew(cos, μ, σ)
+        end
+    end
+
     @testset "min and max" begin
         a, b = GVar(3.0, 1.0), GVar(3.5, 1.0)    # spans [2,4] and [2.5,4.5]
         @test loval(min(a, b)) ≈ 2.0 && hival(min(a, b)) ≈ 4.0