Fast normal random number generators on vector processors

We consider pseudo-random number generators suitable for vector processors. In particular, we describe vectorised implementations of the Box-Muller and Polar methods, and show that they give good performance on the Fujitsu VP2200. We also consider some other popular methods, e.g. the Ratio method of Kinderman and Monahan (1977) (as improved by Leva (1992)), and the method of Von Neumann and Forsythe, and show why they are unlikely to be competitive with the Polar method on vector processors.

💡 Research Summary

The paper addresses the problem of generating normally distributed pseudo‑random numbers efficiently on vector processors, a task that is central to Monte‑Carlo simulations, statistical modelling, and many scientific computing applications. The authors begin by reviewing the classic algorithms for normal variate generation: the Box‑Muller transform, the Polar (or Marsaglia‑Bray) method, the Ratio‑of‑Uniforms technique introduced by Kinderman and Monahan (1977) and later refined by Leva (1992), and the Von Neumann‑Forsythe rejection method. While these algorithms are well‑studied on scalar architectures, their performance characteristics change dramatically on SIMD‑oriented vector machines because of issues such as pipeline stalls, vector length preservation, and memory bandwidth utilization.

The core contribution of the paper is a detailed vectorised implementation of the Box‑Muller and Polar methods on the Fujitsu VP2200, a high‑performance vector supercomputer. For Box‑Muller, the authors note that each generated pair requires a logarithm and two trigonometric evaluations (sine and cosine). On the VP2200 these functions are not hardware‑accelerated, so the authors employ polynomial approximations (Taylor series) and pre‑computed lookup tables to keep the operations vectorisable. Nevertheless, the need for both log and trig evaluations introduces a substantial computational cost and creates divergent control flow when the vector registers contain values that fall outside the domain of the approximations.

The Polar method, by contrast, eliminates the explicit trigonometric calls. It generates two uniform variates u and v, computes s = u² + v², and accepts the pair only if s < 1. Accepted pairs are transformed using r = sqrt(-2 ln s / s) and the final normal variates are x = u r, y = v r. The authors show how to handle the rejection step efficiently in a vector context: after each pass, a mask identifies the elements that failed the s < 1 test; a “compaction” operation packs the surviving elements together, and new uniform variates are generated only for the rejected slots. This approach preserves the vector length for the majority of the computation and minimizes branch divergence, leading to high utilisation of the VP2200’s vector registers (often exceeding 90 % occupancy).

The paper then analyses the Ratio‑of‑Uniforms method and its Leva improvement. Both rely on evaluating a piecewise rational approximation of the inverse cumulative distribution function and involve multiple conditional branches to select the appropriate approximation region. When vectorised, each element of a vector may follow a different branch, causing the vector pipeline to stall and dramatically reducing throughput. Empirical measurements on the VP2200 confirm that the Ratio method achieves only about 45 % vector register utilisation and suffers from irregular memory access patterns due to the need to fetch different coefficients for each element.

Similarly, the Von Neumann‑Forsythe method, which repeatedly rescales a uniform variate until it falls within a target interval, incurs frequent conditional checks and loop iterations. The authors demonstrate that on a vector processor the method’s performance is limited by the same divergence problems, yielding a throughput comparable to the Ratio method but still far below that of the Polar approach.

Performance experiments are presented for each algorithm, using a workload of 10⁸ normal variates. The Polar method consistently outperforms the Box‑Muller implementation by roughly 30 % in wall‑clock time, primarily because it avoids the costly trigonometric evaluations and benefits from the efficient rejection‑masking scheme. Compared with the Ratio and Von Neumann‑Forsythe methods, the Polar algorithm is 2–3 times faster and uses significantly less memory bandwidth, making it the clear winner for vectorised normal variate generation on the VP2200.

Implementation details are discussed in depth. The logarithm required by the Polar method is approximated with a fifth‑order polynomial when hardware support is absent. The square‑root operation is performed using the processor’s native vector sqrt instruction. The authors also describe how they vectorise the underlying uniform generator (a linear congruential generator) to avoid pipeline stalls at the seed‑generation stage. The paper highlights the importance of “vector mask” and “compaction” primitives, which are specific to the VP2200 instruction set but have analogues on other SIMD architectures (e.g., AVX‑512 compress/store masks). By leveraging these primitives, the authors maintain a high degree of parallelism even in the presence of rejection sampling.

In the concluding section, the authors distil three design principles for high‑performance normal random number generation on vector processors: (1) minimise non‑linear function calls (log, trig) by using approximations or lookup tables; (2) handle rejection or conditional branches with vector masks and compaction to keep the vector length stable; and (3) organise memory accesses to be contiguous and regular, thereby reducing bandwidth pressure. The Polar method satisfies all three criteria and, as demonstrated on the Fujitsu VP2200, provides the best overall performance. The authors suggest that the same techniques can be transferred to modern SIMD‑rich platforms such as GPUs and many‑core CPUs, offering a practical roadmap for developers of scientific libraries that require fast, high‑quality normal random variates.