Marsaglia and Tsang

introduced the Ziggurat method. It provides a fast algorithm for generating both normally and exponentially distributed random numbers. The original paper also contains a corresponding implementation in the C language.

The listing in Figure shows this initial implementation. We have removed the code for generating exponentially distributed random numbers as well a comment header from the listing to keep the display more compact. The full listing is also included in the package for reference.

As can be seen from Figure , the Ziggurat algorithm is implemented in bare-bones C code using a number of macros, and two helper functions. The helper functions allow setting a seed, and deal with parameter updates needed in about 2.5% of cases. The actual core component—the function to draw a random number distributed according to thstandard normal distribution—is provided by the macro . Needless to say, using C macros is no longer cosidered . Possible side effects include inadvertent changes in globally visible variables, as well as possible bugs from the macro evaluation.

A more important concern is that this implementation uses types, and explicit bit mapping operations. The code was originally developed for 32-bit operating systems where and are typically four bytes (or 32 bits) wide. Hence the code does not produce correct results on a 64-bit operating system as (signed and unsigned) types are typically eight bytes (or 64 bits) wide (whereas is still 32 bits).

Our modified version introduced below overcomes both issues.

Leong, Zhang et al

show in a comment that the Ziggurat method as introduced by suffers from another weakness due to the generator (by Marsaglia). The authors show via a χ²-test that the generator has a short period of about 2³² − 1, or the four byte limit. Replacing it with the generator (also by Marsaglia) improves the performance beyond this limit.

#define MWC  ((znew<<16)+wnew )
#define SHR3 (jz=jsr, jsr^=(jsr<<13), \
        jsr^=(jsr>>17), jsr^=(jsr<<5),jz+jsr)
#define CONG (jcong=69069*jcong+1234567)
#define KISS ((MWC^CONG)+SHR3)

#define RNOR (hz=KISS, iz=hz&127, \
        (fabs(hz)<kn[iz]) ? hz*wn[iz] : nfix())

Following , Ziggurat code should use an improved uniform generator. Choices are either KISS as suggested initially, or another trusted (and fast) uniform generator such as the Mersenne Twister . Other Open Source implementations (such as the ones discussed below) frequently use the Mersenne Twister as the souce of uniformly distributed random numbers.

Voss’s implementation in GNU GSL

provided another Ziggurat implementation for use in the GNU Scientific Library or GSL . It uses the Mersenne Twister generator by which also avoids the issue identified by in which the originally-used uniform generator had too short a cycle.

notes two differences between his implementation and the original work by . First, he uses only 128 instead of 256 steps which reduces the memory requirements and computational cost at a possible (though presumably minor) loss of precision. Second, he uses an exponential distribution with tail rejection for the base strip, which is motivated by a simpler implementations. Both of these aspects could have implications for the statistical properties of the generator. Voss also appears to be unaware of the work by , yet sidesteps the issue raised by these author by relying on the Mersenne Twister generator.

Here, the implementation from the current GSL sources and file is used, and adapted to the class strucuture detailed in section .

Gretl

The Gretl econometrics program contains another Open Source implementation of the Ziggurat algorithm. The code credits the implementation by described above. review the Gretl implementation and performance of Ziggurat and find it to be satisfactory.

We use the implementation from the file from the current Gretl sources and adapt to the class strucuture detailed in section .

QuantLib

The QuantLib library for quantitative finance contains another open source implementation of Ziggurat. It is provided in the files and . As part of the experimental section, it is made available for further study and use, but not yet part of the default build. As before, we integrate this source into the class structure used here.

R Generators

Before comparing the speed of the different Ziggurat implementations, it is also illustrative to compare the different R generators. Figure provides a comparison.

We see that the Box-Muller generator is the slowest by some margin. However, both Kinderman-Ramage and Ahrens-Dieter are faster than the Inversion method chosen as the default in R. So even before considering Ziggurat generators, R users could reap a speed benefit simply by calling or .

Ziggurat Generators

All Ziggurat generators are significantly faster the the default generator in R which uses a inversion method. Among the Ziggurat generators, we notice that approaches using an external uniform number generator (GNU GSL, GNU Gretl, QuantLib) are all slower than our compact and self-contained implementation which is seen as the fastest method.

Standard Test for Uniform RNG draws

Test for random number generators are often focussed on the case of uniform generators which are the most common type of generators. As detailed for example in , a test proceeds as follow:

Here, we can apply this test for first converting the N(0, 1) distributed values produced by the given Ziggurat implementation to U(0, 1) distributed values by using the inverse of the normal distribution. We are then ready to simulate and test. Figure below displays Q-Q plots for the empirical distribution against the uniform, and displays the p-values of a Kolmogorov-Smirnov as well as a Wilcoxon test.

We see that five of the six generators pass the test. In the case of the original generator, we can see the departure from the expected diagonal clearly once we draw more than 4.2 × 10⁹ numbers (which is the limit of the representation of an unsigned four-byte nunber). However, only the Kolmogorov-Smirnov test can formally reject; the Wilcoxon test appear to lack sufficient power in this setting. The QuantLib generator is seen as suspicious which p-value just below a conventional rejection level.

Normal Test for Normal RNG draws

We can propose a simpler variant of the test outlined in the previous section. As the random numbers we are drawing are following a N(0, 1) distribution, the sum of their values follows a N(0, n) distribution. This allows us to skip one inversion step:

Results, shown in Figure , are qualitatively similar to the result discussed above. Kolmogorov-Smirnov rejects for the generator. However, we note that the Wilcoxon test now has a lower p-value—we would now reject at conventional test levels. The QuantLib implementation is now rejected by the Wilcoxon test but not the Kolmogorov-Smirnov.

χ² test

As can be seen in Figure , the original proposal by fails as was shown by . All other tests pass again.

C++ Implementation

Preceding work by and also contained implementations in the C language. These versions were implemented in just a few lines, and used idioms common to C programmers such as macros and global variables.

C++ programming style permits encapsulation in order to avoid possible collisions and side-effects. Moreover, by using a modest amount of object-oriented programming we can use a class structure with a common base class to express commonalities between the implementations. The following code segment shows the virtual base class used here.

#ifndef RcppZiggurat__Zigg_h
#define RcppZiggurat__Zigg_h

#include <cmath>
#include <stdint.h>         // or cstdint (C++11)

namespace Ziggurat {

   class Zigg {
   public:
      virtual ~Zigg() {};
      virtual void setSeed(const uint32_t s) = 0;
      // no getSeed() as GSL has none
      virtual double norm() = 0;
   };
}

#endif

As shown in the preceding code segment, we provide two user-accessible functions to obtain a normal random deviate, and to set the seed, respectively. The actual implementation uses portable types such as , an unsigned 32-bit integer provided by the C header file , which provides correct and identical results on both 32-bit and 64-bit operating systems.

Each actual implementation can then encapsulate its state variable as a private variable inaccessible to other functions. Such a small core to each class also makes it feasible to provide a Ziggurat generator in each thread in a parallel execution framework.

Our Ziggurat implementation has no external dependencies and can therefore be used in other projects. The testing framework used for this note has a single dependency on the GNU GSL as the generator by is used via its GSL implementations. The generator and testing framework in the corresponding R package have a build-dependency on R, and are of course accessed by R. But the generator discussed here could equally well be used in standalone programs or with other scripting languages.

In the package, the Rcpp Attributes feature of the C++/R integration package are used to access instances of the corresponding generator class.

#include <Rcpp.h>

#include <Ziggurat.h>

static Ziggurat::Ziggurat::Ziggurat zigg;

// [[Rcpp::export]]
Rcpp::NumericVector zrnorm(int n) {
    Rcpp::NumericVector x(n);
    for (int i=0; i<n; i++) {
        x[i] = zigg.norm();
    }
    return x;
}

// [[Rcpp::export]]
void zsetseed(unsigned long int s) {
    zigg.setSeed(s);
    return;
}

In this particular reference implementation, we have chosen a namespace for the entire project. Within this namespace, we opted to provided an extra namespace layer for each generator as some of these generators still use global variables—which are therefore shielded in their own namespace. For example, for the generator, we use . Next is the name of the actual class—which in the case of the reference implementation shown above is once again leading to the triple use of the term. Actual deployment of the Ziggurat generator (without comparison to other implementations and concerns about interaction between variables, particularly for the older implementations having global variables) can of course be used with a single namespace.

The remainder shows how two functions and are provided via the attribute as described in the Rcpp Attributes vignette of the package .