A revision of the subtract-with-borrow random number generators

TL;DR

This paper reimplements the popular subtract-with-borrow random number generator RANLUX as a large-integer linear congruential generator, achieving faster performance and improved seeding capabilities on modern computers.

Contribution

It introduces a novel, efficient modular multiplication method for RANLUX, enabling faster random number generation and unique sequence seeding.

Findings

Significant speed improvement over traditional RANLUX

Comparable performance to other high-quality generators

Fast state skipping for unique sequence seeding

Abstract

The most popular and widely used subtract-with-borrow generator, also known as RANLUX, is reimplemented as a linear congruential generator using large integer arithmetic with the modulus size of 576 bits. Modern computers, as well as the specific structure of the modulus inferred from RANLUX, allow for the development of a fast modular multiplication -- the core of the procedure. This was previously believed to be slow and have too high cost in terms of computing resources. Our tests show a significant gain in generation speed which is comparable with other fast, high quality random number generators. An additional feature is the fast skipping of generator states leading to a seeding scheme which guarantees the uniqueness of random number sequences.

Tables

Table 1. Table 1: The multiplier A ≡ a p mod m 𝐴 modulo superscript 𝑎 𝑝 𝑚 A\equiv a^{p}\bmod m which corresponds to the RANLUX strategy to waste p − 24 𝑝 24 p-24 random numbers before next 24 numbers will be delivered to a user as well as the associated luxury levels.

luxury level	$p$	$a^{p}modm$
0	24	\seqsplitfffffffffffffffffffffffefffffffffffffffffffffffffffffffffffffffffffffffffffffffffffe000000000000000000000001000000000000000000000000000000000000
1	48	\seqsplit000000000000000000000002ffffffffffffffffffffffff000000000000000000000000000000000001fffffffffffffffffffffffc000000000000000000000001000000000001
2	97	\seqsplitffffff000000000008000000000009fffffffffffefffffffffff1000000000000000000000006ffffff000004fffffffffff6ffffffffffec000000000001000000000015000001
3	223	\seqsplit00028b000000000bba00000000026cfffffffff8e4fffffffff96000000000027b0000000007d0fffffffffe25ffffffffeef0fffffffffa0a000000000942000000000ba6000000
4	389	\seqsplit0df0600000002ee0020000000b9242ffffffdf6604ffffffe4ab160000000d92ab0000001e93f2fffffff593cfffffffb9c8a6ffffffe525740000002c38960000002ecac9000000
$-$	1024	\seqsplite1754cefa19deea6f58651c8ac11b437ba841c49eca3003ff0ef508f058cfdab6105ca16980e6a3ab12a823219e1cd0007281433953609f1cc9c5ca19cf7f0c6d3899b14b7c5ee90
$-$	2048	\seqsplitb48c187cf5b22097492edfcc0cc8e753ff74e54107684ed2256c3d3c662ea36c20b2ca60cb78c5096d8a15a13bee7cb0e64dcb31c48228ec4cec2c78af55c101ed7faa90747aaad9

Table 2. Table 2: CPU clocks needed to multiply two 576 bit numbers for generic functions from GMP (mul_basecase_*) and the author’s 9 × 9 9 9 9\times 9 limbs multiplications for various CPU families.

	Core2	Haswell	Skylake
mul_basecase_core2	370.9	218.5	197.1
mul_basecase_coreihwl	n.a.	205.6	203.7
mul_basecase_coreibwl	n.a.	n.a.	163.0
mul9x9	534.7	198.8	192.0
mul9x9mulx	n.a.	162.1	154.5
mul9x9mulxadox	n.a.	n.a.	119.4

Table 3. Table 3: CPU clocks needed for multiplication of two 576 bit numbers modulo 2 576 − 2 240 + 1 superscript 2 576 superscript 2 240 1 2^{576}-2^{240}+1 .

	Haswell	Skylake
mulmod9x9mulx	201.2	191.6
mulmod9x9mulxadox	n.a.	155.4

Table 4. Table 4: CPU clocks needed to sum 10 9 superscript 10 9 10^{9} random numbers at Haswell CPU normalized to a single number. All tests were compiled by gcc v6.2.1 with -O3 optimization flag.

	double	float
dummy	9.1	9.1
std::minstd_rand	35.1	20.2
std::mt19937_64	36.0	37.0
std::ranlux24_base	47.2	26.0
std::ranlux48_base	24.4	26.1
std::ranlux24	387.0	197.7
std::ranlux48	640.5	638.4
gsl_ranlxs0		84.2
gsl_ranlxs1		125.9
gsl_ranlxs2		215.5
gsl_ranlxd1	213.8
gsl_ranlxd2 ($p=397$)	394.0
gsl_ranlux ($p=223$)		185.0
gsl_ranlux389 ($p=389$)		315.4
TRandom1 ($p=389$)		362.7
RANLUX ($p=389$)		378.3
ranlxs (array, SSE, $p=397$)		50.4
ranlxd (array, SSE, $p=397$)	95.0
RANLUX++	29.2	20.0
RANLUX++ (array)	24.7	15.7

Table 5. Table 5: CPU clocks needed to sum 10 9 superscript 10 9 10^{9} single precision random numbers of the subtract-with-borrow algorithm with the skipping at Haswell CPU normalized to a single number. The skipping mode corresponds to the highest luxury level 4 of RANLUX. All tests were compiled by gcc v6.2.1 with -O3 optimization flag.

scalar	scalar(asm)	SSE2	AVX2
40.3	28.5	12.0	7.95

Equations

$$x_{i+1}=(a\cdot x_i+c)\mathrm{mod}m,$$

$$Y_{i+1}=\{\begin{array}{ll}(y_2,\dots ,y_r,\Delta ,0),& \text{if}\Delta \ge 0\\ (y_2,\dots ,y_r,\Delta +b,1),& \text{otherwise}\end{array}$$

$$x_i=x(Y_i)\equiv \sum _{j=1}^r{y}_j{b}^{j-1}-\sum _{j=1}^s{y}_{r-s+j}{b}^{j-1}+k$$

$$\underset{p\text{times}}{\underbrace{a\cdot (a\cdot (\dots )\mathrm{mod}m)\mathrm{mod}m}}=(a^p\mathrm{mod}m)\cdot x\mathrm{mod}m=A\cdot x\mathrm{mod}m,$$

$$x_s=(a^n{)}^s\cdot x_0\mathrm{mod}m,$$

Full text

††journal: Computer Physics Communications

PROGRAM SUMMARY/NEW VERSION PROGRAM SUMMARY

Program Title: RANLUX++

Licensing provisions: GPLv3

Programming language: C++, C, Assembler

1 Introduction

The well known Linear Congruential Generator (LCG) is a recurrent sequence of numbers calculated as follows:

[TABLE]

where $x_{0}$ is the initial state or seed, $a$ – the multiplier, $c$ – the increment and $m$ – the modulus. The particular choice of the parameters $a$, $c$ and $m$ with period $q$ – the minimal number when $x_{q}=x_{0}$, can be found in the literature [1]. Commonly used LCGs are limited to $m\leq 2^{64}$, and have poor statistical properties. Thus they are not used for Monte-Carlo physical simulations.

This situation can be mitigated when $m$ reaches several hundreds or even thousand bits. The cost of the increased range of $m$ is to deal with arbitrary precision integer arithmetic which was believed to be prohibitively expensive for practical purposes. In the last two decades there has been tremendous progress in modern central processor units (CPU) especially for personal computers (PC) which can be employed for long arithmetic.

We have explored the possibility to use the long arithmetic in LCG to improve the quality of generated random numbers and found that, despite a substantial increase in calculations, the time to generate a single random number is not proportionally risen. In fact for some parameters, the computational time decreased compared to ordinary LCGs with machine word modulus size.

2 Subtract-with-borrow generator

At this point no specific constraints on $a$, $c$ and $m$ parameters of LCG have been applied. As a good starting point we choose the subtract-with-borrow generator first introduced in [2] and the intimate connection with LCG has been shown as a part of the period calculation. The algorithm has been extensively studied in [3] to improve statistical quality of generated numbers. Based on this study the generator RANLUX [4] was developed and now it is widely used in physics simulations as well as in other fields where random numbers with high statistical quality are required. However the current method employed by RANLUX to achieve the high quality makes it one of the slowest generators on the market.

The definition of the subtract-with-borrow generator is the following: let $b$ some integer greater than 1 also called the base and vector $Y=(y_{1},\ldots,y_{r},k)$ with the length $r+1$, where $0\leq y_{i}<b$ and $k$ or the carry equals 0 or 1. Then define a recursive transformation of the vector $Y_{i}$ with the rule:

[TABLE]

where $\Delta=y_{r-s+1}-y_{1}-k$ and $r$ and $s$ also called the lags. As shown in the work [5], this recursion is equivalent to LCG with the modulus $m=b^{r}-b^{s}+1$, the multiplier $a=m-(m-1)/b$ and $c=0$ with the relation:

[TABLE]

In the RANLUX generator the lags $r=24$ and $s=10$ with the base $b=2^{24}$ are chosen among other suggested parameters in [2], and thus the modulus $m=b^{24}-b^{10}+1$ is a prime number and the multiplier $a=m-(m-1)/b=b^{24}-b^{23}-b^{10}+b^{9}+1$. With those parameters the period $q$ is equal to $(m-1)/48$.

Due to the selected base $b$ the natural choice to keep the generator state is a vector of length 24 composed of 24-bit numbers. This implementation uses the properties of the modulus $m$ to avoid long arithmetic calculations, and a single step equivalent to one modular multiplication $(a\cdot x\bmod m)$ that requires only subtraction of two 24-bit numbers and carry propagation. In the original FORTRAN implementation, 24-bit numbers were stored as floats to avoid at that time, a high cost integer-to-float conversion.

2.1 Remainder

The simple structure of the modulus $m$ allows us to calculate the remainder using only additions, subtractions and bit shifts. The modulus $m$ and thus the generator state $x$ have size of $24\cdot 24=576$ bits and fits into 9 64-bit machine words. The result of the product $z=a\cdot x$ fits into 18 64-bit machine words which can be represented as a 48 element array of 24-bit numbers: $z=[z_{0},z_{1},\ldots,z_{46},z_{47}]$. The number $r$ obtained by the procedure shown in Algorithm 1 is congruent to $z\bmod m$ and $r<b^{24}$. Note the product $c\cdot m$ is also only bit shifting due to the simple structure of $m$. The calculation of $c$ is a sum of carry bits of each arithmetic operation.

2.2 Skipping

Examining the result of a single step of Eq. 1 one can note that the main part of the number $x_{i}$ is preserved in its successor $x_{i+1}$ which is just rotated by 24 bits. This strong correlation is the reason of the poor statistical quality of the original subtract-with-borrow generator [2]. The bright idea developed in [3] is to apply the transformation (2) many times to break correlations between nearby states before using the state for actual physical simulation. The drawback of this method is obvious – all intermediate states have to be explicitly calculated even if they are not needed. Despite the single step being simple with small resource consumption, good statistical quality requires several hundred steps thus in total, the skipping requires a lot of time. This is a luxury to spend resources and not use the results. Thus so-called luxury levels were introduced as aliases for how many generatated numbers have to be wasted.

Using Eq. 1 we can efficienty skip numbers since all $p$ recurrent steps collapse to a single multiplication:

[TABLE]

where the factor $A\equiv(a^{p}\bmod m)$ is precomputed and thus the cost to calculate the next state with or without skipping is the same. Any state in the entire period $q=(m-1)/48\approx 10^{171}$ can be calculated in no more than $2\times\log_{2}(q)\approx 1140$ long multiplications using fast exponentiation by squaring which takes order of tens of $\mu$sec on modern CPUs.

In the Table 2.2 the precomputed values of $A\equiv(a^{p}\bmod m)$ where the values of $p$ is taken from [4] are shown for illustrative purposes. In the initial rows, long chains of 0 or 1 in binary representation are clearly visible and this can be interepreted such that for each bit of the state $x_{i+p}$ only a few bits of the state $x_{i}$ contributes. Even at the highest luxury level 4 there are still some patterns observable and a demading user maybe not be completely satisfied. For such user the two last rows would be more attractive especially since it is for free! Such chaotic multipiers mean that if any single bit of the state $x_{i}$ is changed in the next step the altered state will be absolutely different from the unaltered one.

With explicit long multiplication, there is no need to keep the multiplier $A$ as a power of $a$, it can be adjusted to get the full period, $m-1$. As an example the number $(a^{2048}+13\bmod m)$ is a primitive root modulo $m$ and with this multiplier all numbers in the range $1\ldots m-1$ will appear in the sequence only once with any initial $x_{0}$ from the same range.