ziggurat

README

ziggurat

The Ziggurat Algorithm is an extremely fast algorithm for the random sampling from arbitrary probability distributions. Just pass your distribution as input, and get a random number generator as output.

What

The goal of random number generation is to sample from the area under the probability distribution curve. The Ziggurat algorithm does this by covering the probability distribution up into equal areas with a large number of rectangles, which just barely exceed the area of the probability distribution, then samples from a random point in a randomly selected rectangle. There is a very small chance that the sampled point lies outside the probability distribution, in which case we simply retry (with the same rectangle). With enough rectangles, you can construct them such that the probability of landing outside the probability distribution is extremely small, and in the vast majority of cases, sampling a random point from the distribution is equivalent to picking a random point in rectangle.

This code enables the automated construction of these extremely fast random number generators for unimodal, univariate probability distributions. This code is closely integrated with gonum, so you can supply most gonum distributions as input and get your fast(er) random number generator as output. You can also write your own distributions as long as they fulfill the ziggurat.Distribution interface.

For fastest results, I recommend using xorshift as your random source.

Why

Fast random number generation is very important for large-scale Monte Carlo simulations, commonly used in scientific computing and statistics. Random number generation is often the bottleneck in these simulations, and I needed this code to run my own simulations at as large a scale as possible.

For the normal distribution and exponentiatial distribution, the Go standard library already uses pre-built hand-optimized Ziggurat algorithms to generate fast random numbers, so this library will not outperform them; however this library enables the creation of similarly-performing algorithms for arbitrary probability distributions (e.g. Triangle, Gamma, and Students' T distributions), and is much faster than existing libraries in those cases.

How

package main

import (
	"github.com/argusdusty/ziggurat"
	"github.com/vpxyz/xorshift/xoroshiro128plus"
	"gonum.org/v1/gonum/stat/distuv"
)

func main() {
	src := xoroshiro128plus.NewSource(1)
	distribution := distuv.UnitNormal // Swap this for most gonum univariate distributions
	rng := ziggurat.ToZiggurat (Symmetric)(distribution, src) // The normal distribution is symmetric, so we can use the more efficient symmetric ziggurat
	randomNormalValue := rng.Rand()
}

Note that gonum 1.16.0 is required due to the use of math/rand/v2.

Benchmarks

ziggurat>go test -run="^$" -bench=.
goos: windows
goarch: amd64
pkg: github.com/argusdusty/ziggurat
cpu: AMD Ryzen 9 5900X 12-Core Processor

Distribution	Algorithm	RNG	Time
Beta (alpha=4, beta=4)	Ziggurat	Default	21.84ns/op
Beta (alpha=4, beta=4)	Ziggurat	xoroshiro128+	11.45ns/op
Beta (alpha=4, beta=4)	Ziggurat (Symmetric)	Default	9.347ns/op
Beta (alpha=4, beta=4)	Ziggurat (Symmetric)	xoroshiro128+	4.557ns/op
Beta (alpha=4, beta=4)	Gonum	Default	48.69ns/op
Beta (alpha=4, beta=4)	Gonum	xoroshiro128+	38.19ns/op
Gamma (alpha=1)	Ziggurat	Default	8.775ns/op
Gamma (alpha=1)	Ziggurat	xoroshiro128+	3.927ns/op
Gamma (alpha=1)	Gonum	Default	11.64ns/op
Gamma (alpha=1)	Gonum	xoroshiro128+	8.195ns/op
Half-Normal	Ziggurat	Default	8.431ns/op
Half-Normal	Ziggurat	xoroshiro128+	3.560ns/op
Normal	Ziggurat	Default	21.22ns/op
Normal	Ziggurat	xoroshiro128+	10.58ns/op
Normal	Ziggurat (Symmetric)	Default	8.611ns/op
Normal	Ziggurat (Symmetric)	xoroshiro128+	3.868ns/op
Normal	Gonum	Default	11.24ns/op
Normal	Gonum	xoroshiro128+	6.978ns/op
Normal	Stdlib	Default	8.816ns/op
Normal	Stdlib	xoroshiro128+	3.975ns/op
Student's t (dof=5)	Ziggurat	Default	22.55ns/op
Student's t (dof=5)	Ziggurat	xoroshiro128+	11.80ns/op
Student's t (dof=5)	Ziggurat (Symmetric)	Default	9.703ns/op
Student's t (dof=5)	Ziggurat (Symmetric)	xoroshiro128+	5.016ns/op
Student's t (dof=5)	Gonum	Default	42.80ns/op
Student's t (dof=5)	Gonum	xoroshiro128+	35.68ns/op
Triangle (a=0, b=1, c=0)	Ziggurat	Default	8.416ns/op
Triangle (a=0, b=1, c=0)	Ziggurat	xoroshiro128+	3.572ns/op
Triangle (a=0, b=1, c=0)	Gonum	Default	17.85ns/op
Triangle (a=0, b=1, c=0)	Gonum	xoroshiro128+	15.66ns/op

Documentation

Index

• Constants
• func ToSymmetricZiggurat(distribution Distribution, src rand.Source) distuv.Rander
• func ToZiggurat(distribution Distribution, src rand.Source) distuv.Rander
• type Distribution

Constants

const (
	ZIGGURAT_BIT_LENGTH = 10 // 10 is the largest feasible number here, otherwise we start to eat into floating point accuracy
	ZIGGURAT_N          = 1 << ZIGGURAT_BIT_LENGTH
)

Functions

func ToSymmetricZiggurat

func ToSymmetricZiggurat(distribution Distribution, src rand.Source) distuv.Rander

func ToZiggurat

func ToZiggurat(distribution Distribution, src rand.Source) distuv.Rander

Types

type Distribution

type Distribution interface {
	Mode() float64              // The x value for which the PDF is maximized.
	Prob(x float64) float64     // The PDF function for this distribution.
	Survival(x float64) float64 // The survival function (1-CDF) for this distribution.
	Quantile(p float64) float64 // The quantile function (integral(CDF)) for this distribution.
}