Geometric() could be 5x faster.

[LilithHafner]

julia> Geometric()
Geometric{Float64}(p=0.5)

julia> @b Geometric() rand
10.111 ns

julia> @b Geometric() rand(_, 1000)
8.625 μs (3 allocs: 7.875 KiB)

julia> struct OneHalf <: Real end

julia> Base.rand(rng::Random.AbstractRNG, ::Geometric{OneHalf}) = leading_zeros(rand(rng, UInt64))

julia> Geometric2() = Geometric{OneHalf}(OneHalf())
Geometric2 (generic function with 1 method)

julia> Geometric2()
Geometric{OneHalf}(p=OneHalf())

julia> @b Geometric2() rand
2.283 ns

julia> @b Geometric2() rand(_, 1000)
1.117 μs (3 allocs: 7.875 KiB)

julia> using StatsBase # Correctness sanity check

julia> countmap(rand(Geometric(), 1000))
Dict{Int64, Int64} with 9 entries:
  0 => 506
  4 => 30
  5 => 9
  6 => 8
  2 => 137
  7 => 3
  8 => 1
  3 => 56
  1 => 250

julia> countmap(rand(Geometric2(), 1000))
Dict{Int64, Int64} with 11 entries:
  5  => 13
  8  => 2
  1  => 263
  0  => 479
  6  => 11
  9  => 2
  3  => 67
  7  => 5
  4  => 28
  15 => 1
  2  => 129

This is accurate to within 2^-64 (i.e. one would need to take approximately 2^64 samples to demonstrate that this implementation is not faithful to the distribution it claims to sample from; 55261506563602523656 to demonstrate at p=0.05).

It would come at the downside of making Geometric() a different type than Geometric(0.5).

[devmotion]

AFAICT this is the same idea as outlined in #1883, isn't it?

[LilithHafner]

Yes. Same idea. One thing that's new here is the idea to use leading_zeros.