Computer Physics Communications 184 (2013) 1889–1897

Contents lists available at SciVerse ScienceDirect

Computer Physics Communications

journal homepage: www.elsevier.com/locate/cpc

Finding the best portable congruential random number generators

Fatin Sezgin

a,∗

, Tevfik Metin Sezgin

b,1

a

Bilkent University, Ankara, Turkey

b

Koc University, Istanbul, Turkey

a r t i c l e i n f o

Article history:

Received 14 December 2012

Received in revised form

3 March 2013

Accepted 12 March 2013

Available online 22 March 2013

Keywords:

Lattice structure

Linear congruential generator

Portability

Random number

Spectral test

a b s t r a c t

Linear congruential random number generators must have large moduli to attain maximum periods, but

this creates integer overflow during calculations. Several methods have been suggested to remedy this

problem while obtaining portability. Approximate factoring is the most common method in portable

implementations, but there is no systematic technique for finding appropriate multipliers and an

exhaustive search is prohibitively expensive. We offer a very efficient method for finding all portable

multipliers of any given modulus value. Letting M = AB+C, the multiplier A gives a portable result if B−C

is positive. If it is negative, the portable multiplier can be defined as A =

⌊

M/B

⌋

. We also suggest a method

for discovering the most fertile search region for spectral top-quality multipliers in a two-dimensional

space. The method is extremely promising for best generator searches in very large moduli: 64-bit sizes

and above. As an application to an important and challenging problem, we examined the prime modulus

2

63

−25, suitable for 64-bit register size, and determined 12 high quality portable generators successfully

passing stringent spectral and empirical tests.

© 2013 Elsevier B.V. All rights reserved.

1. Introduction

Random numbers are essential tools in many applications, such

as simulation, education, cryptography, the arts, numerical analy-

sis, computer programming, VLSI testing, recreation and sampling.

Besides some physical and tabular sources, there are several deter-

ministic computational techniques to produce random sequences

of data, including congruential, shift register, lagged Fibonacci, in-

verse and cellular automata generators. Linear congruential gen-

erators have attracted the attention of many researchers because

they are efficient, easy to implement and their theory is well doc-

umented.

Mixed linear congruential generators produce a sequence of

integers {X

i

} defined by the recursion

X

n+1

= AX

n

+ c(mod M) (1)

with appropriately defined integer constants A, c, M and an initial

value X

0

. The special case of c = 0 has particular significance and

is called a multiplicative congruential generator. The maximum

possible period of these generators, M − 1, is limited by the

maximum integer size available in the computer and compiler. For

example, 32-bit computers can process integers up to 2

31

−1. This

∗

Corresponding author. Tel.: +90 312 290 5010; fax: +90 312 266 4607.

E-mail addresses: fatin@bilkent.edu.tr (F. Sezgin), mtsezgin@ku.edu.tr

(T.M. Sezgin).

1

Tel.: +90 212 338 1540; fax: +90 212 338 1548.

Mersenne prime is a popular modulus for multiplicative generators

for 32-bit arithmetic. Since the product AX

n

can produce very large

values, several methods have been suggested in the literature to

prevent integer overflow.

This study considers the portable implementation, relying on

the so-called approximate factoring method by using the constants

A, B and C, satisfying the expression M = AB + C such that B > C.

In the article, the term ‘‘portable’’ is used specifically to refer to

approximately factorable generators.

In Section 2 we offer a systematic search method for efficiently

determining all multipliers suitable for approximate factoring. Sec-

tion 3 gives a fruitful search interval for the best-quality con-

gruential generators with respect to a spectral test criterion. In

Section 4, we present application examples. In Section 5, we ap-

ply the method to the modulus M = 2

63

− 25, suitable for 64-bit

computers, and report 12 very good multipliers. In the last section,

we summarize the results and suggest some future applications.

2. Implementations for portable generators

Consider the multiplicative congruential generator X

n+1

=

48271X

n

(mod (2

31

− 1)). The calculations required by this gener-

ator may produce integers up to 1.473 × 2

46

and cause an integer

overflow in 32-bit arithmetic. There are several studies proposing

and reviewing techniques for the correct implementation of ran-

dom number generators, such as [1–8].

0010-4655/$ – see front matter © 2013 Elsevier B.V. All rights reserved.

http://dx.doi.org/10.1016/j.cpc.2013.03.013

1890 F. Sezgin, T.M. Sezgin / Computer Physics Communications 184 (2013) 1889–1897

2.1. Portable implementation by approximate factoring

Among many implementation methods, approximate factoring

has become highly popular for its effectiveness and simplicity. This

method relies on the decomposition of the modulus M as

M = AB + C (2)

where

B =

M

A

(3)

and C = M − AB, subject to the condition B > C, and with ⌊·⌋

denoting the floor function. It was first used by Wichmann and

Hill [9] in an informal remark. Sezgin [7] and Park and Miller [10]

independently proved the validity of the method. According to this

implementation, the code X

n+1

= AX

n

(mod M) can be written as:

K = int(X

n

/B)

X

n+1

= A ∗ (X

n

(mod B)) − C ∗ K

if (X

n+1

.lt.0)X

n+1

= X

n+1

+ M.

For example, since 2147483647 = 48271 ∗ 44488 + 3399, the

generator

X

n+1

= 48271X

n

(mod (2

31

− 1))

can be coded as:

k = ix/44488

ix = 48271 ∗ mod(ix, 44488) − 3399 ∗ k

if(ix.lt.0)ix = ix + 2147483647.

We state the following theorem for the condition of portability:

Theorem 1. For a modulus M, the multiplier A is portable if and only

if

M

A

>

M

A+1

.

Proof. Assume that A is portable. Then in (2) we have B > C . This

implies that B > M − AB. Inserting the value of B from (3) we get

M

A

> M − A

M

A

or

M

A

(A + 1) > M.

Dividing both sides by A +1 gives the asserted result. The opposite

is also true. Assume that

M

A

>

M

A+1

. Then the operations in

the reverse direction give

M

A

(A + 1) > M and this implies

BA + B > M = AB + C , giving B > C .

Corollary 1. All multipliers A ≤

√

M are portable.

Proof. It is easy to see that A =

√

M is portable. Because in this

case the above theorem is satisfied:

M

√

M

>

M

√

M + 1

.

On the other hand, when A <

√

M we will have B ≥ A. But since

A > C, we get B > C.

Corollary 2. There is no portable multiplier larger than

M

2

.

Table 1

Finding multipliers suitable for approximate factoring in modulus M = 103.

A B C B − C BX > X is Portable A =

13 7 12 −5 5 1 14 = 13 + 1

15 6 13 −7 7 2 17 = 15 + 2

18 5 13 −8 8 2 20 = 18 + 2

21 4 19 −15 15 4 25 = 21 + 4

26 3 25 −22 22 8 34 = 26 + 8

35 2 33 −31 31 16 51 = 35 + 16

Proof.

M

2

is portable, because in this case we will get the B value

B =

M

A

=

M

M

2

= 2.

Therefore the C value will be C = M − AB = M − 2

M

2

= 1

or 0 depending on the M being odd or even. We also see that

A =

M

2

+ k is not portable for positive integers k. Because if

k ≥ 1, we will get B = 1, and B > C implies C = 0. But this is not

the case since we get C = M − AB = M − A ≥ B.

We note that we always have C = M(mod A) < A, but most of

the time C > B. The portability condition is met by few multipliers.

For example, there are only 92679 multipliers satisfying this

condition for M = 2

31

−1 and this number constitutes 0.0043% of

all multiplier candidates. Hence, an exhaustive search for portable

generators is fruitless; some general rules must be determined

for an efficient search algorithm. In this respect, we present the

following remarks:

1. All multipliers ≤

√

M can be used for approximate factoring.

2. The last portable multiplier by approximate factoring is

⌊

M/2

⌋

.

3. Above

√

M, in some cases we may have B −C < 0. If that is the

case, we need an X value such that XB + (B − C) > 0. In order

to get a portable multiplier, this X value must be added to A.

Example. Let us take M = 103, a small modulus, for demonstra-

tive purposes. For values up to A = 10, all multipliers can be used in

approximate factoring. We observe that 11 and 12 are also satisfac-

tory. But A = 13 gives B = 7 and C = 12, causing a negative B −C

value: B − C = −5. We must have 7X > 5, getting X = 1. When

this occurs, 13 + 1 = 14 can be implemented by approximate fac-

toring. Following the same method, we get five more approximate

factoring generators, as summarized in Table 1.

The following theorem explains the method for finding the next

portable multiplier when B < C:

Theorem 2. When a multiplier

√

M < A < M/2 gives B < C, the

next portable generator can be obtained by using the multiplier

M

B

.

Proof. From (2) we get B−C = B−M +AB. If this value is negative,

we need to add an X value to A to get a positive B −C. It means that

we will have BX > −(B −C) = M −B −AB or BX > M −(A +1)B.

The X value will be

X >

M − (A + 1)B

B

.

We can take the smallest possible integer X as

X =

M − (A + 1)B

B

+ 1.

This can be simplified as:

X =

M

B

− (A + 1) + 1

=

M

B

− A.

F. Sezgin, T.M. Sezgin / Computer Physics Communications 184 (2013) 1889–1897 1891

Fig. 1. The distribution of B − C values for the modulus M = 19997 until A = 1000.

Therefore, when B − C < 0, the next portable multiplier must be

taken simply as A + X =

M

B

.

For example: With M = 103, A = 35;X =

103−36×2

2

+ 1 =

31

2

+ 1 = 16; therefore 35 + 16 = 51 is portable. Since B =

M

A

=

103

35

= 2, the same multiplier can be obtained simply as

M

B

=

103

2

= 51.

A FORTRAN program for finding portable multipliers is pre-

sented in Listing 1. It receives modulus values from the keyboard

and stops when 0 is given.

integer a, b

1 print*, "M=? Give 0 to stop"

read*, M

if (m.eq.0) stop

maxa=m/2

a=sqrt(float(m))

print*, "All multipliers less than or equal to ",a," are

portable"

print*, "Other portable multipliers are:"

2 a=a+1

if (a.gt.maxa) goto 1

b=m/a

if(b.gt.m-a*b) then

print*, a

else

a=m/b

print*, a

end if

goto 2

end

Listing 1. The FORTRAN program finding portable multipliers.

2.2. The distribution of B − C

B − C can be expressed as a function of A

B − C =

M

A

− M + A

M

A

.

It has a positive value until A ≤

√

M. After this limit negative

values start to be seen. For large values of A, B−C values are located

on straight lines with a slope

M

A

and the proportion of negative

points increases. The distribution of B − C is demonstrated for

modulus M = 19997 in Fig. 1. The B − C values are positive for

all multipliers A ≤

√

M = 141. The first negative value is seen

at A = 146, and as the multiplier gets larger, the proportion of

negative values increases, and the values points cluster on straight

lines. For large A, at each line there is only one positive value;

this situation explains the inefficiency of the exhaustive search for

portable multipliers. The figure shows multipliers until A = 1000.

The remaining multipliers have the same pattern of lines, only with

one positive point. Until

√

M, all multipliers are portable for all M.

After the square root, the number of portable multipliers declines.

The cumulative distribution of portable multipliers is depicted in

Fig. 2 for M = 2

31

− 1. Each point represents 10 000 multiplier

candidates. An inspection of the figure shows that 50% of the

portables are below A = 50000. Only 10% of the portables are

above A = 230000. The percentage of portable multipliers above

A = 500000 is lower than 5%. In the light of the discussion above,

an exhaustive search for large portable multipliers is impractical

for two reasons:

1. These multipliers are very scarce and searching for them is a

waste of time.

2. As will be seen in Section 3, the spectral quality of large portable

multipliers is extremely poor.

3. Selecting good multipliers

After finding multipliers suitable for implementation by ap-

proximate factoring, the next step is to determine the full-period

condition and the quality of these multipliers. Random number

generators must be subjected to several theoretical and empirical

tests before they can be used in serious applications. In this respect,

the spectral test is a very reliable theoretical tool to distinguish be-

tween bad and good congruential generators. This test is explained

in detail by Knuth [11]. Letting 0 < A < M and A be relatively

prime to M, the test determines the values of

ν

t

= min

t

i=1

s

i

2

(4)

for 2 ≤ t ≤ T , given that

t

i=1

s

i

A

i−1

≡ 0 (Mod M) (5)

where s

i

are integers 0 ≤ s

i

< M and (s

1

, . . . , s

t

) = (0, . . . , 0).

1892 F. Sezgin, T.M. Sezgin / Computer Physics Communications 184 (2013) 1889–1897

Fig. 2. Cumulative distribution of portable multiplier proportions for modulus M = 2

31

− 1.

Table 2

The upper and lower limits of the best multipliers for the modulus M = 2

31

− 1.

Percentile Lower limit of A Upper limit of A Lower limit of S

2

Best 20% 42374 = 0.914

√

M 71370 = 1.540

√

M 0.851

Best 10% 44963 = 0.970

√

M 67157 = 1.449

√

M 0.903

Best 5% 46207 = 0.997

√

M 65197 = 1.407

√

M 0.928

Best 1% 48008 = 1.036

√

M 53790 = 1.161

√

M 0.964

Letting d

t

= 1/ν

t

represent the maximum distance between

adjacent hyper-planes determined by the points of the lattice in

a t-dimensional space, and d

∗

t

represent the lower bound of this

distance, a figure of merit called the normalized spectral test value

is defined as

S

t

= d

∗

t

/d

t

. (6)

For our portable multipliers, if A is small, the identity (5) can be

defined as A − A = 0 (Mod M), implying s

1

= 1 and s

2

= −A. By

(4) these values will produce the following minimum ν

2

value:

ν

2

= min

2

i=1

s

i

2

=

1 + A

2

.

Since d

∗

2

is the constant 1/M

1/2

(4/3)

1/4

, from (6), the normal-

ized spectral value for t = 2 will be:

S

2

=

√

1 + A

2

√

M(4/3)

1/4

.

This is approximately 0.93A/

√

M and explains the steady in-

crease in the quality of small multipliers when they approach the

square root of the modulus. S

2

will reach its maximum value S

2

= 1

when A =

√

M/0.93 = 1.075

√

M. When the multiplier exceeds

√

M by a large margin, this situation produces very small ν

2

values.

Because a very large A creates a smaller B value, and since C is even

smaller, taking the identity (5) as AB + C = M (Mod M) = 0, the

generator produces a very small ν

2

2

= B

2

+ C

2

value. The search

for the best two-dimensional lattice structure can be restricted

to a narrow interval. Investigating the distribution of S

2

values in

Fig. 5 for portable multipliers 30000 < A < 68000 in a generator

with M = 2

31

− 1, we can notice that the peak S

2

= 0.9985 is

reached at A = 49883 = 1.076

√

M. Sezgin [12] presented per-

centiles of normalized spectral test for dimensions t = 2, . . . , 8.

These percentiles remain unchanged across modulus values. Ta-

ble 2 presents the intervals for the best multiplier percentile limits

for t = 2 in M = 2

31

− 1.

4. Some numerical examples

In Section 2, the distribution of B−C was described, and a graph

was presented in Fig. 1 for a small modulus. Here, in Fig. 3, we

demonstrate the pattern of large multipliers for a larger modulus,

M = 2

31

− 1, from A = 300000 to A = 310000. There are only

231 non-negative values (portable multipliers) in this range. The

proportion of portables is 2.31%. As the multipliers increase, this

proportion will fall sharply.

For large modulus values like M = 2

48

, 2

63

or 2

128

, only partial

searches are conducted in the random number literature, because,

as mentioned in Section 3, conducting an exhaustive search for

good portable generators is a tedious task. Even with restrictions

on the search, it requires a great amount of computer time. Begin-

ning in 2006 and ending in 2009, Tang and Chang [13] conducted

an exhaustive search for good 64-bit congruential random number

generators under some restrictions on an IBM computer running

Linux and compiled with gcc-02. After this effort, they found only

35 portable multipliers having a spectral performance value above

S

T

> 0.760. Fig. 4 depicts A/

√

M values for these best generators.

It is interesting to observe that these generators could have been

found in a very short time, because all the multipliers are extremely

close to

√

M.

We conducted a search for M = 2

31

−1. Fig. 5 shows the distri-

bution of S

2

values for portable multipliers 30000 < A < 68000.

Actually, the graph starts as a straight line from A = 2 and con-

tinues near to

√

M in this manner. In later parts, there is an inter-

esting fractal structure and a steady decrease of the spectral value

with increasing A. The remaining part is depicted in Fig. 6 until

A = 2000000 exhibiting the sharp decline in the spectral value.

An example for a larger modulus is depicted in Fig. 7: M =

4503599090499601 in a narrow interval near the square root for

S

2

> 0.91 values. Results agree with our remark in Section 3. The

search gave the spectral maximum A as 72179138, which is also

1.076

√

M.

As proven in [12,14], the squared figure of merit defined in

(4), ν

t

2

, changes as a polynomial of degree 2(t − 1) and this

causes very chaotic behavior in higher dimensions. Therefore, in

determining the best multipliers, one should conduct the search

in two-dimensional spectral values first. Fig. 8 shows S

2

, S

3

and

S

4

values for a 63-bit modulus: M = 9223372012704246017.

Although its square root, A = 3037000496, is very good for a

two-dimensional space and gives S

2

= 0.931 there, it causes a

very small spectral value in dimensions above two. This fact can

be observed from the graphs of S

3

and S

4

in Fig. 8. The situation is

similar for higher dimensions. But fortunately, the bad lattice range

is very short. In the worst case, within a range of length about 2900,

F. Sezgin, T.M. Sezgin / Computer Physics Communications 184 (2013) 1889–1897 1893

Fig. 3. The distribution of B − C for multipliers from A = 300000 to A = 310000 in modulus M = 2

31

− 1.

Fig. 4. Distribution of A/

√

M values for the 35 best multipliers found by Tang and Chang for 64-bit random number generators.

Fig. 5. Distribution of S

2

values for portable multipliers 30000 < A < 68000 in a generator with M = 2

31

− 1.

S

3

reaches values larger than 0.9. In higher dimensions, the range of

bad lattices is negligible. As a result, we suggest searching around

the square root of the modulus, where the spectral test value of the

two-dimensional space reaches its maximum value and provides

the most fertile search area.

5. An application to find 64-bit congruential generators

L’Ecuyer and Simard [15] point out that linear congruential

generators with moduli up to 2

61

fail several tests in the strin-

gent test package TestU01. It would be interesting to investigate

the quality of larger moduli by this package. Since 64-bit reg-

isters can allow the representation of integer arithmetic up to

2

63

− 1, we searched the best portable random number genera-

tors for the prime modulus 2

63

− 25. All multipliers determined

in this study produce full period. If the multiplier A is a primi-

tive root of the modulus M, the generator will have full period

when M is prime. This is a necessary and sufficient condition. In

our case the factorization of M − 1 is 2

63

− 26 = 2 × 3

4

× 17 ×

23 ×319279 ×456065899. Therefore, A is a full-period multiplier

in cases where, for all prime factors q, A

M−1/q

= 1 Mod(M). In

this search, we determined portable full-period multipliers in the

top 1% region of a two-dimensional space and investigated their