Theziggurat algorithm is analgorithm forpseudo-random number sampling. Belonging to the class ofrejection sampling algorithms, it relies on an underlying source of uniformly-distributed random numbers, typically from apseudo-random number generator, as well as precomputed tables. The algorithm is used to generate values from amonotonically decreasingprobability distribution. It can also be applied tosymmetricunimodal distributions, such as thenormal distribution, by choosing a value from one half of the distribution and then randomly choosing which half the value is considered to have been drawn from. It was developed byGeorge Marsaglia and others in the 1960s.

A typical value produced by the algorithm only requires the generation of one random floating-point value and one random table index, followed by one table lookup, one multiply operation and one comparison. Sometimes (2.5% of the time, in the case of a normal orexponential distribution when using typical table sizes)^{[citation needed]} more computations are required. Nevertheless, the algorithm is computationally much faster^{[citation needed]} than the two most commonly used methods of generating normally distributed random numbers, theMarsaglia polar method and theBox–Muller transform, which require at least one logarithm and onesquare root calculation for each pair of generated values. However, since the ziggurat algorithm is more complex to implement it is best used when large quantities of random numbers are required.

The termziggurat algorithm dates from Marsaglia's paper with Wai Wan Tsang in 2000; it is so named because it is conceptually based on covering the probability distribution with rectangular segments stacked in decreasing order of size, resulting in a figure that resembles aziggurat.

The ziggurat algorithm is a rejection sampling algorithm; it randomly generates a point in a distribution slightly larger than the desired distribution, then tests whether the generated point is inside the desired distribution. If not, it tries again. Given a random point underneath a probability density curve, itsx coordinate is a random number with the desired distribution.

The distribution the ziggurat algorithm chooses from is made up ofn equal-area regions;${\displaystyle n-1}$ rectangles that cover the bulk of the desired distribution, on top of a non-rectangular base that includes the tail of the distribution.

Given a monotone decreasing probability density function${\displaystyle f(x)}$, defined for all${\displaystyle x\ge 0}$ , the base of the ziggurat is defined as all points inside the distribution and below${\displaystyle y_1=f(x_1)}$. This consists of a rectangular region from${\displaystyle (0,0)}$ to${\displaystyle (x_1,y_1)}$, and the (typically infinite) tail of the distribution, where${\displaystyle x\ge x_1(\text{and }y\le y_1)}$.

This layer (call it layer 0) has areaA. On top of this, add a rectangular layer of width${\displaystyle x_1}$ and height${\displaystyle A/x_1}$, so it also has area${\displaystyle A}$. The top of this layer is at height${\displaystyle y_2=y_1+A/x_1}$ and intersects the density function at a point${\displaystyle (x_2,y_2)}$, where${\displaystyle y_2=f(x_2)}$. This layer includes every point in the density function between${\displaystyle y_1}$ and${\displaystyle y_2}$, but (unlike the base layer) also includes points such as${\displaystyle (x_1,y_2)}$ which are not in the desired distribution.

Further layers are then stacked on top. To use a precomputed table of size${\displaystyle n}$ (n = 256 is typical), one chooses${\displaystyle x_1}$such that${\displaystyle x_n=0}$, meaning that the top box, layer${\displaystyle n-1}$, reaches the distribution's peak at${\displaystyle (0,f(0))}$ exactly.

Layer${\displaystyle i}$ extends vertically in range${\displaystyle [y_i,y_{i+1}]}$, and can be divided into two regions horizontally: the (generally larger) portion in range${\displaystyle [0,x_{i+1}]}$ which is entirely contained within the desired distribution, and the (small) portion in range${\displaystyle [x_{i+1},x]}$, which is only partially contained.

Ignoring for a moment the problem of layer 0, and given uniform random variables${\displaystyle U_0}$ and${\displaystyle U_1\in [0,1)}$, the ziggurat algorithm can be described as:

1. Choose a random layer.
2. Let${\displaystyle x=U_0{x}_i}$.
3. If, return${\displaystyle x}$.
4. Let${\displaystyle y=y_i+U_1(y_{i+1}-y_i)}$.
5. Compute${\displaystyle f(x)}$. If, return${\displaystyle x}$.
6. Otherwise, choose new random numbers and go back to step 1.

Step 1 amounts to choosing a low-resolutiony coordinate. Step 3 tests if thex coordinate is clearly within the desired density function without knowing more about the y coordinate. If it is not, step 4 chooses a high-resolution y coordinate, and step 5 does the rejection test.

With closely spaced layers, the algorithm terminates at step 3 a very large fraction of the time. For the top layer${\displaystyle n-1}$, however, this test always fails, because${\displaystyle x_n=0}$.

Layer 0 can also be divided into a central region and an edge, but the edge is an infinite tail. To use the same algorithm to check if the point is in the central region, generate a fictitious${\displaystyle x_0=A/y_1}$. This will generate points with with the correct frequency, and in the rare case that layer 0 is selected and${\displaystyle x\ge x_1}$, use a special fallback algorithm to select a point at random from the tail. Because the fallback algorithm is used less than one time in a thousand, speed is not essential.

Thus, the full ziggurat algorithm for one-sided distributions is:

1. Choose a random layer.
2. Let${\displaystyle x=U_0{x}_i}$.
3. If, return${\displaystyle x}$.
4. If${\displaystyle i=0}$, generate a point from the tail using the fallback algorithm.
5. Let${\displaystyle y=y_i+U_1(y_{i+1}-y_i)}$.
6. Compute${\displaystyle f(x)}$. If, return${\displaystyle x}$.
7. Otherwise, choose new random numbers and go back to step 1.

For a two-sided distribution, the result must be negated 50% of the time. This can often be done conveniently by choosing${\displaystyle U_0\in [-1,1]}$ and, in step 3, testing if.

Because the ziggurat algorithm only generatesmost outputs very rapidly, and requires a fallback algorithm whenever${\displaystyle x>x_1}$ it is always more complex than a more direct implementation. The specific fallback algorithm depends on the distribution.

For an exponential distribution, the tail looks just like the body of the distribution. One way is to fall back to the most elementary algorithmE = −ln(U₁) and letx = x₁ − ln(U₁). Another is to call the ziggurat algorithmrecursively and addx₁ to the result.

For a normal distribution, Marsaglia suggests a compact algorithm:

1. Letx = −ln(U₁)/x₁.
2. Lety = −ln(U₂).
3. If 2y >x², returnx + x₁.
4. Otherwise, go back to step 1.

Sincex₁ ≈ 3.5 for typical table sizes, the test in step 3 is almost always successful. Since −ln(U₁) is an exponentially distributed variate, an implementation of the exponential distribution may be used.

The algorithm can be performed efficiently with precomputed tables ofx_i andy_i =f(x_i), but there are some modifications to make it even faster:

• Nothing in the ziggurat algorithm depends on the probability distribution function being normalized (integral under the curve equal to 1), removingnormalizing constants can speed up the computation off(x).
• Most uniform random number generators are based on integer random number generators which return an integer in the range [0,  2³² − 1]. A table of 2⁻³²x_i lets you use such numbers directly forU₀.
• When computing two-sided distributions using a two-sidedU₀ as described earlier, the random integer can be interpreted as a signed number in the range [−2³¹, 2³¹ − 1], and a scale factor of 2⁻³¹ can be used.
• Rather than comparingU₀x_i tox_i +1 in step 3, it is possible to precomputex_i +1/x_i and compareU₀ with that directly. IfU₀ is an integer random number generator, these limits may be premultiplied by 2³² (or 2³¹, as appropriate) so an integer comparison can be used.
• With the above two changes, the table of unmodifiedx_i values is no longer needed and may be deleted.
• When generatingIEEE 754 single-precisionfloating point values, which only have a 24-bit mantissa (including the implicit leading 1), the least-significant bits of a 32-bit integer random number are not used. These bits may be used to select the layer number. (See the references below for a detailed discussion of this.)
• The first three steps may be put into aninline function, which can call an out-of-line implementation of the less frequently needed steps.

It is possible to store the entire table precomputed, or just include the valuesn,y₁,A, and an implementation off⁻¹(y) in thesource code, and compute the remaining values when initializing the random number generator.

As previously described, you can findx_i =f⁻¹(y_i) andy_i +1 = y_i + A/x_i. Repeatn − 1 times for the layers of the ziggurat. At the end, you should havey_n = f(0). There will be someround-off error, but it is a usefulsanity test to see that it is acceptably small.

When actually filling in the table values, just assume thatx_n = 0 andy_n = f(0), and accept the slight difference in layern − 1's area as rounding error.

Given an initial (guess at)x₁, you need a way to compute the areat of the tail for whichx >x₁. For the exponential distribution, this is juste^−x₁, while for the normal distribution, assuming you are using the unnormalizedf(x) =e^−x2/2, this is√π/2 erfc(x/√2). For more awkward distributions,numerical integration may be required.

With this in hand, fromx₁, you can findy₁ =f(x₁), the areat in the tail, and the area of the base layerA =x₁y₁ + t.

Then compute the seriesy_i andx_i as above. Ify_i >f(0) for anyi <n, then the initial estimatex₁ was too low, leading to too large an areaA. Ify_n <f(0), then the initial estimatex₁ was too high.

Given this, use aroot-finding algorithm (such as thebisection method) to find the valuex₁ which producesy_n−1 as close tof(0) as possible. Alternatively, look for the value which makes the area of the topmost layer,x_n−1(f(0) − y_n−1), as close to the desired valueA as possible. This saves one evaluation off⁻¹(x) and is actually the condition of greatest interest.

Christopher D. McFarland has proposed a further-optimized version.^[1] This applies three algorithmic changes, at the expense of slightly larger tables.

First, the common case considers only the rectangular portions, from (0,y_i −1) to (x_i, y_i) The odd-shaped regions to the right of these (mostly almost triangular, plus the tail) are handled separately. This simplifies and speeds up the algorithm'sfast path.

Second, the exact area of the odd-shaped regions is used; they arenot rounded up to include the entire rectangle to (x_i −1, y_i). This increases the probability that the fast path will be used.

One major consequence of this is that the number of layers is slightly less thann. Even though the area of the odd-shaped portions is taken exactly, the total adds up to more than one layer's worth. The area per layer is adjusted so that the number of rectangular layers is an integer. If the initial 0 ≤ i < n exceeds the number of rectangular layers, phase 2 proceeds.

If the value sought lies in any of the odd-shaped regions, thealias method is used to choose one, based on its true area. This is a small amount of additional work, and requires additional alias tables, but chooses one of the layers' right-hand sides.

The chosen odd-shaped region is rejection sampled, but if a sample is rejected, the algorithm doesnot return to the beginning. The true area of each odd-shaped region was used to choose a layer, so the rejection-sampling loop stays in that layer until a point is chosen.

Third, the almost-triangular shape of most odd-shaped portions is exploited, although this must be divided into three cases depending on thesecond derivative of the probability distribution function in the selected layer.

If the function is convex (as the exponential distribution is everywhere, and the normal distribution is for |x| > 1), then the function is strictly contained within the lower triangle. Two unit uniform deviatesU₁ andU₂ are chosen, and before they are scaled to the rectangle enclosing the odd-shaped region, their sum is tested. IfU₁ + U₂ > 1, the point is in the upper triangle and can be reflected to (1−U₁, 1−U₂). Then, ifU₁ + U₂ < 1−ε, for some suitable toleranceε, the point is definitely below the curve and can immediately be accepted. Only for points very close to the diagonal is it necessary to compute the distribution functionf(x) to perform an exact rejection test. (The toleranceε should in theory depend on the layer, but a single maximum value can be used on all layers with little loss.)

If the function is concave (as the normal distribution is for |x| < 1), it includes a small portion of the upper triangle so reflection is impossible, but points whose normalized coordinates satisfyU₁ +U₂ ≤ 1 can be immediately accepted, and points for whichU₁ +U₂ > 1+ε can be immediately rejected.

In the one layer which straddles |x| = 1, the normal distribution has aninflection point, and the exact rejection test must be applied if 1−ε <U₁ + U₂ < 1+ε.

The tail is handled as in the original Ziggurat algorithm, and can be thought of as a fourth case for the shape of the odd-shaped region to the right.

• Pseudorandom number generators
• Non-uniform random numbers
• Statistical algorithms

• Articles with short description
• Short description matches Wikidata
• All articles with unsourced statements
• Articles with unsourced statements from September 2011
• Articles with unsourced statements from December 2021