This is a MELG-64 library for Fortran-2008.
This library is based on S. Harase and T. Kimoto's implementation.

The following contents were originally present in the original README.

The 64-bit Maximally Equidistributed F₂-Linear Generators with Mersenne Prime Period (MELG-64) are 64-bit Mersenne-Twister-type pseudorandom number generators developed between 2014 and 2017, and the corresponding paper was published on ACM TOMS in 2018.

S. Harase and T. Kimoto, "Implementing 64-bit maximally equidistributed F₂-linear generators with Mersenne prime period", ACM Transactions on Mathematical Software, Volume 44, Issue 3, April 2018, Article No. 30, 11 pp. Artcle

CPUs and operating systems are moving from 32 to 64 bits, and hence it is important to have good 64-bit pseudorandom number generators (PRNGs) designed to fully exploit these word lengths.

The 32-bit Mersenne Twister (MT) MT19937 (Matsumoto and Nishimura, 1998) is one of the most widely used PRNGs, but it is not completely optimized in terms of high-dimensional uniformity, which is a theoretical criterion of PRNGs. The 32-bit WELL generators (Panneton et. al., 2006) was developed in order to overcome this weakness.

However, for 64-bit PRNGs, MT19937-64 (Nishimura, 2000) and SFMT19937 using SIMD (Saito and Matsumoto, 2008), etc., have been proposed, but there exists no 64-bit MT-type long-period linear PRNG completely optimized for high-dimensional uniformity, such as a variant of WELL generators.

In this page, we introduce 64-bit maximally equidistributed F₂-linear PRNGs (MELG-64) that are optimal in this respect and have speeds equivalent to 64-bit Mersenne Twisters.

MELG19937-64 has the following properties:

• Very long period 2¹⁹⁹³⁷-1 ≈ 10⁶⁰⁰⁰;
• High-dimensional uniformity completely optimized;
• Fast generation competitive with MT19937-64;
• Memory size requiring only 312 words (similarly to MT19937-64).

We provide the codes for MELG-64 with various period lengths from 2⁶⁰⁷-1 to 2⁴⁴⁴⁹⁷-1. The jump-ahead algorithm is also implemented in order to obtain disjoint streams in parallel computing. (The default skip size is 2²⁵⁶.)

Please click the "Code" button in the upper right of the content pane, and clone this page:

git clone https://github.com/sharase/melg-64.git

or download the zip file. If you want to use MELG19937-64, for example,

cd melg-64
cd melg19937-64
gcc melg19937-64.c -o melg19937-64 -O3 -Wall
./melg19937-64

Before using, please initialize the state by the function init_by_array64(init, length) in a similar way as Mersenne Twisters.

The high-dimensional uniformity is a theoretical criterion for PRNGs, which is assessed via the dimension of equidistribution with v-bit accuracy as follows:

Let

x₀, x₁, ..., x_P-1, x_P = x₀, ...

be a unsigned w-bit binary integer sequence with period P, where w is the word size of the intended machine. Let trunc_v(x_i)
denote the number formed by the v most significant bits of x_i. Consider the kv-bit vectors for the entire period:

(trunc_v(x_i), trunc_v(x_i+1), ..., trunc_v(x_i+k-1)),   i = 0, ..., P-1 .

A pseudorandom sequence x_i of w-bit integers of period P is said to be k-dimensionally equidistributed with v-bit accuracy if each of the 2^kv possible combinations of bits occurs the same number of times over the whole period P, except for the all-zero combination that occurs once less often.

The largest value of k with this property is called the dimension of equidistribution with v-bit accuracy, denoted by k(v).

This definition is based on the assumption that the higher digits are large numbers. In particular, the dimension of equidistribution ensures that the output values with the v most significant bits are uniformly distributed up to dimension k(v). Thus, as a criterion of uniformity, larger values of k(v) for each 1 ≤ v ≤ w is desirable.

Now we have a trivial upper bound

k(v) ≤ ⌊ log₂(P+1) / v ⌋

for each v = 1, 2, ...,w. Define the sum of the gaps

Δ := ∑ ⌊ log₂(P+1) / v ⌋ - k(v)),

where the sum is over all 1 ≤ v ≤ w.

If Δ = 0, the generator is said to be maximally equidistributed (ME).

The aim of our study is to design maximally equidistributed F₂-linear PRNGs with similar speed as 64-bit Mersenne Twisters.

We compare the following MT-type PRNGs corresponding to 64-bit output sequences:

• MELG19937-64: the 64-bit integer output of our proposed generator;
• MT19937-64: the 64-bit integer output of the 64-bit Mersenne Twister (downloaded from here);
• MT19937-64 (ID3): the 64-bit integer output of a 64-bit Mersenne Twister based on a five-term recursion (ID3) (Nishimura 2000);
• SFMT19937-64 (without SIMD): the 64-bit integer output of the SIMD-oriented Fast Mersenne Twister SFMT19937 without SIMD (Saito and Matsumoto 2008);
• SFMT19937-64 (with SIMD): the 64-bit integer output of the foregoing with SIMD (Saito and Matsumoto 2008).

We measure the CPU time (in seconds) taken to generate 10⁹ 64-bit unsigned integers. The following table summarizes the timing and the figures of merit Δ and N₁. See the remark below for the definition of N₁.

• CPU time (Intel): Intel Core i7-3770 (3.40GHz) Linux gcc compiler with -O3
• CPU time (AMD): AMD Phenom II X6 1045T (2.70 GHz) Linux gcc compiler with -O3

SFMT19937 is very fast but Δ for SFMT19937 is large. (In fact, the SFMT generators are optimized under the assumption that one will mainly be using 32-bit output sequences. For double-precision floating-point numbers, dSFMT is faster than SFMT and is also improved from the viewpoint of the dimensions of equidistribution with v-bit accuracy.)

Remark:

For MT-type PRNGs, each bit of the sequence obeys a linear feedback shift register generator with the characteristic polynomial P(x) over the two-element field F₂. Let N₁ be the number of nonzero coefficients of P(x). In addition to the excellent equidistribution, as a secondary criterion, N₁ should be large enough. This criterion implies that the PRNG avoids a long-lasting impact for poor initialization, such as 0-excess states (Panneton et al, 2006).

For more information, please see the presentation slides here.

• S. Harase and T. Kimoto, "Implementing 64-bit maximally equidistributed F₂-linear generators with Mersenne prime period", ACM Transactions on Mathematical Software, Volume 44, Issue 3, April 2018, Article No. 30, 11 Pages. Artcle
• S. Harase, "Conversion of Mersenne Twister to double-precision floating-point numbers", Mathematics and Computers in Simulation, Volume 161, July 2019, Pages 76-83. Article

日本語の発表スライドもこちらに公開しました。Youtubeの解説動画も作成しました。

This work was partially supported by JSPS KAKENHI Grant Numbers JP18K18016, JP26730015, JP26310211, JP15K13460, JP12J07985.