Improvements for drawing random numbers

Project.toml

@@ -1,10 +1,11 @@
 name = "Quaternions"
 uuid = "94ee1d12-ae83-5a48-8b1c-48b8ff168ae0"
-version = "0.4.2"
+version = "0.4.3"
 
 [deps]
 DualNumbers = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [compat]
 DualNumbers = "0.5, 0.6"

README.md

@@ -28,6 +28,8 @@ Implemented functions are:
     cos
     sqrt
     linpol  (interpolate between 2 normalized quaternions)
+    rand
+    randn
 
 [Dual quaternions](http://en.wikipedia.org/wiki/Dual_quaternion) are an extension, combining quaternions with
 [dual numbers](https://github.com/scidom/DualNumbers.jl).
@@ -53,6 +55,7 @@ further implemented here:
     exp
     log
     sqrt
+    rand
 
 [Octonions](http://en.wikipedia.org/wiki/Octonion) form the logical next step on the Complex-Quaternion path.
 They play a role, for instance, in the mathematical foundation of String theory.
@@ -70,3 +73,5 @@ They play a role, for instance, in the mathematical foundation of String theory.
     exp
     log
     sqrt
+    rand
+    randn

src/DualQuaternion.jl

@@ -169,3 +169,7 @@ end
 
 dualquatrand() = dualquat(quatrand(), quatrand())
 ndualquatrand() = normalize(dualquatrand())
+
+function rand(rng::AbstractRNG, ::Random.SamplerType{DualQuaternion{T}}) where {T<:Real}
+    return DualQuaternion{T}(rand(rng, Quaternion{T}), rand(rng, Quaternion{T}), false)
+end

src/Octonion.jl

@@ -159,3 +159,22 @@ function sqrt(o::Octonion)
 end
 
 octorand() = octo(randn(), randn(), randn(), randn(), randn(), randn(), randn(), randn())
+
+function rand(rng::AbstractRNG, ::Random.SamplerType{Octonion{T}}) where {T<:Real}
+  Octonion{T}(rand(rng, T), rand(rng, T), rand(rng, T), rand(rng, T),
+              rand(rng, T), rand(rng, T), rand(rng, T), rand(rng, T), false)
+end
+
+function randn(rng::AbstractRNG, ::Type{Octonion{T}}) where {T<:AbstractFloat}
+  Octonion{T}(
+      randn(rng, T) * INV_SQRT_EIGHT,
+      randn(rng, T) * INV_SQRT_EIGHT,
+      randn(rng, T) * INV_SQRT_EIGHT,
+      randn(rng, T) * INV_SQRT_EIGHT,
+      randn(rng, T) * INV_SQRT_EIGHT,
+      randn(rng, T) * INV_SQRT_EIGHT,
+      randn(rng, T) * INV_SQRT_EIGHT,
+      randn(rng, T) * INV_SQRT_EIGHT,
+      false,
+  )
+end

src/Quaternion.jl

@@ -184,8 +184,22 @@ function linpol(p::Quaternion, q::Quaternion, t::Real)
     end
 end
 
-quatrand()  = quat(randn(), randn(), randn(), randn())
-nquatrand() = normalize(quatrand())
+quatrand(rng = Random.GLOBAL_RNG)  = quat(randn(rng), randn(rng), randn(rng), randn(rng))
+nquatrand(rng = Random.GLOBAL_RNG) = normalize(quatrand(rng))
+
+function rand(rng::AbstractRNG, ::Random.SamplerType{Quaternion{T}}) where {T<:Real}
+    Quaternion{T}(rand(rng, T), rand(rng, T), rand(rng, T), rand(rng, T), false)
+end
+
+function randn(rng::AbstractRNG, ::Type{Quaternion{T}}) where {T<:AbstractFloat}
+    Quaternion{T}(
+        randn(rng, T) * 1//2,
+        randn(rng, T) * 1//2,
+        randn(rng, T) * 1//2,
+        randn(rng, T) * 1//2,
+        false,
+    )
+end
 
 ## Rotations
 

src/Quaternions.jl

@@ -5,7 +5,11 @@ module Quaternions
   import Base: +, -, *, /, ^
   import Base: abs, abs2, angle, conj, cos, exp, inv, isfinite, log, real, sin, sqrt
   import Base: convert, promote_rule, float
+  import Base: rand, randn
   import LinearAlgebra: norm, normalize
+  using Random
+
+  Base.@irrational INV_SQRT_EIGHT 0.3535533905932737622004 sqrt(big(0.125))
 
   include("Quaternion.jl")
   include("Octonion.jl")

test/test_Quaternion.jl

@@ -1,6 +1,7 @@
 
 using Quaternions: argq
 using LinearAlgebra
+using Random
 
 # creating random examples
 sample(QT::Type{Quaternion{T}}) where {T <: Integer} = QT(rand(-100:100, 4)..., false)
@@ -139,3 +140,83 @@ for _ in 1:100
         @test q ⊗ linpol(q1, q2, t) ≈ linpol(q ⊗ q1, q ⊗ q2, t)
     end
 end
+
+@testset "random quaternions" begin
+    @testset "quatrand" begin
+        rng = Random.MersenneTwister(42)
+        q1 = quatrand(rng)
+        @test q1 isa Quaternion
+        @test !q1.norm
+
+        q2 = quatrand()
+        @test q2 isa Quaternion
+        @test !q2.norm
+    end
+
+    @testset "nquatrand" begin
+        rng = Random.MersenneTwister(42)
+        q1 = nquatrand(rng)
+        @test q1 isa Quaternion
+        @test q1.norm
+
+        q2 = nquatrand()
+        @test q2 isa Quaternion
+        @test q2.norm
+    end
+
+    @testset "rand($H)" for H in (Quaternion, DualQuaternion, Octonion)
+        rng = Random.MersenneTwister(42)
+        q1 = rand(rng, H{Float64})
+        @test q1 isa H{Float64}
+        @test !q1.norm
+
+        q2 = rand(rng, H{Float32})
+        @test q2 isa H{Float32}
+        @test !q2.norm
+
+        qs = rand(rng, H{Float64}, 1000)
+        @test eltype(qs) === H{Float64}
+        @test length(qs) == 1000
+        xs = map(qs) do q
+            if q isa DualQuaternion
+                return [real(q.q0); Quaternions.imag(q.q0); real(q.qe); Quaternions.imag(q.qe)]
+            else
+                return [real(q); Quaternions.imag(q)]
+            end
+        end
+        xs_mean = sum(xs) / length(xs)
+        xs_var = sum(x -> abs2.(x .- xs_mean), xs) / (length(xs) - 1)
+        @test all(isapprox.(xs_mean, 0.5; atol=0.1))
+        @test all(isapprox.(xs_var, 1/12; atol=0.01))
+    end
+
+    @testset "randn($H)" for H in (Quaternion, Octonion)
+        rng = Random.MersenneTwister(42)
+        q1 = randn(rng, H{Float64})
+        @test q1 isa H{Float64}
+        @test !q1.norm
+
+        q2 = randn(rng, H{Float32})
+        @test q2 isa H{Float32}
+        @test !q2.norm
+
+        qs = randn(rng, H{Float64}, 10000)
+        @test eltype(qs) === H{Float64}
+        @test length(qs) == 10000
+        xs = map(qs) do q
+            if q isa DualQuaternion
+                return [real(q.q0); Quaternions.imag(q.q0); real(q.qe); Quaternions.imag(q.qe)]
+            else
+                return [real(q); Quaternions.imag(q)]
+            end
+        end
+        xs_mean = sum(xs) / length(xs)
+        xs_var = sum(x -> abs2.(x .- xs_mean), xs) / (length(xs) - 1)
+        @test all(isapprox.(xs_mean, 0; atol=0.1))
+        if H === Quaternion
+            @test all(isapprox.(xs_var, 1/4; atol=0.1))
+        else
+            @test all(isapprox.(xs_var, 1/8; atol=0.1))
+        end
+    end
+end