Skip to content

Generalized golden sequences, a form of low discrepancy sequence or quasi random numbers

License

Notifications You must be signed in to change notification settings

mschauer/GoldenSequences.jl

Repository files navigation

GoldenSequences.jl

Generalized golden sequences, a form of low discrepancy sequence or quasi random numbers See Martin Roberts: The Unreasonable Effectiveness of Quasirandom Sequences for background.

The d-dimensional sequence follows

x[i] = (x[i-1] .+ z) .% true, x[0] = x0

where

z = [ϕ[k]^(-i) for i in 1:d]

and ϕ[k] solves ϕ[k]^(d+1) = ϕ[k] + 1 (with ϕ[1] the golden mean.)

Golden sequence

julia> GoldenSequence(0.0)[1]
0.6180339887498949

Shifted golden sequence starting in 0.5

julia>  GoldenSequence(0.5)[0]
0.5

julia>  GoldenSequence(0.5)[1]
0.1180339887498949

GoldenSequence returns an infinite iterator:

julia> collect(take(GoldenSequence(0.0), 10))
10-element Array{Float64,1}:
 0.0                
 0.6180339887498949
 0.2360679774997898
 0.8541019662496847
 0.4721359549995796
 0.09016994374947451
 0.7082039324993694
 0.3262379212492643
 0.9442719099991592
 0.5623058987490541

Random colors: Low discrepancy series are good choice for (quasi-) random colors

using Colors
n = 20
c = map(x->RGB(x...), (take(GoldenSequence(3), n))) # perfect for random colors

Colors

2D golden sequence

julia>  GoldenSequence(2)[1]
(0.7548776662466927, 0.5698402909980532)

As low discrepancy series these number are well distributed (left), better than random numbers (right):

using Makie
n = 155
x = collect(Iterators.take(GoldenSequence(2), n))
p1 = scatter(x, markersize=0.02)
y = [(rand(),rand()) for i in 1:n]
p2 = scatter(y, markersize=0.02, color=:red)
vbox(p1, p2)

Quasi-random vs. random

Cartesian Golden Sequence

With a bit of effort, one can use Golden Sequences to generate spacefilling quasirandom sequences of cartesian indices. For example GoldenCartesianSequence((m[1], m[2])) will create a 2D Cartesian sequence corresponding to (approximate) samples of the Golden sequence in 1:m[1] x 1:m[2].

For that, GoldenCartesianSequence((m[1], ..., m[d])) will create full period linear congruential generators (LCG) x[i+1] = (x[i] + c[k]) % m[k] approximating phi[d]^[-k] by c[k]/m[k] such that c[k] and m[k] are coprime.

If m[1], ..., m[d] are coprime themselves, these LCG will have together period prod(m) and the sequence will be space filling, that is sort(collect(take(GoldenCartesianSequence(m), prod(m)))) == CartesianIndices(m)[:].

This means that if m[k] is the denominator of a good rational approximation c[k]//m[k] ≈ ϕ[d]^(-k), then the indices will be well distributed even for large i.

For example if m = (2819, 3508):

Quasi-random cartesian indices

The image shows the fraction of the first 0.00005 (red) and the first 0.002 indices (black) in the GoldenCartesianSequence((2819, 3508)).

In short, good m are coprime denominators of fractions given by the function rationalize

d = 2
m = @. denominator(rationalize(GoldenSequences.phis[d]^(-(1:d)), tol=0.0000001))
julia> m
2-element Array{Int64,1}:
 2819
 3508

julia> gcd(m[1], m[2])
1

There are some connections here to Knuth's multiplicative hashing method

hash(i) = mod(i*2654435769, 2^32)

where 2654435769 is approximately 2^32*ϕ.

Interface

GoldenSequence(n::Int) # Float64 n-dimensional golden sequence
GoldenSequence(x0::Number) # 1-d golden sequence shifted by `x0`
GoldenSequence(x0) # length(x)-d golden sequence shifted/starting in 'x0'

A flower

Flower petals grow in spots not covering older petals, the new spot is at an angle given by the golden sequence.

using Colors
using Makie
n = 20
c = map(x->RGB(x...), (take(GoldenSequence(3), n))) # perfect for random colors
x = collect(take(GoldenSequence(0.0), n))
petals = [(i*cos(2pi*x), i*sin(2pi*x)) for (i,x) in  enumerate(x)]
scatter(reverse(petals), color=c, markersize=10*(n:-1:1))

Flower petals

About

Generalized golden sequences, a form of low discrepancy sequence or quasi random numbers

Resources

License

Stars

Watchers

Forks

Packages

No packages published

Contributors 4

  •  
  •  
  •  
  •  

Languages