1 //******************************************************************************
2 //
3 // File: Mcg1Random.java
4 // Package: edu.rit.util
5 // Unit: Class edu.rit.util.Mcg1Random
6 //
7 // This Java source file is copyright (C) 2008 by Alan Kaminsky. All rights
8 // reserved. For further information, contact the author, Alan Kaminsky, at
9 // ark@cs.rit.edu.
10 //
11 // This Java source file is part of the Parallel Java Library ("PJ"). PJ is free
12 // software; you can redistribute it and/or modify it under the terms of the GNU
13 // General Public License as published by the Free Software Foundation; either
14 // version 3 of the License, or (at your option) any later version.
15 //
16 // PJ is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
18 // A PARTICULAR PURPOSE. See the GNU General Public License for more details.
19 //
20 // Linking this library statically or dynamically with other modules is making a
21 // combined work based on this library. Thus, the terms and conditions of the GNU
22 // General Public License cover the whole combination.
23 //
24 // As a special exception, the copyright holders of this library give you
25 // permission to link this library with independent modules to produce an
26 // executable, regardless of the license terms of these independent modules, and
27 // to copy and distribute the resulting executable under terms of your choice,
28 // provided that you also meet, for each linked independent module, the terms
29 // and conditions of the license of that module. An independent module is a module
30 // which is not derived from or based on this library. If you modify this library,
31 // you may extend this exception to your version of the library, but you are not
32 // obligated to do so. If you do not wish to do so, delete this exception
33 // statement from your version.
34 //
35 // A copy of the GNU General Public License is provided in the file gpl.txt. You
36 // may also obtain a copy of the GNU General Public License on the World Wide
37 // Web at http://www.gnu.org/licenses/gpl.html.
38 //
39 //******************************************************************************
40 package edu.rit.util;
41
42 import java.io.Serial;
43
44 /**
45 * Class Mcg1Random provides a default pseudorandom number generator (PRNG)
46 * designed for use in parallel scientific programming. To create an instance of
47 * class Mcg1Random, either use the <code>Mcg1Random()</code> constructor, or use
48 * the static <code>getInstance(long,String)</code> method in class {@linkplain
49 * Random}.
50 * <P>
51 * Class Mcg1Random uses L'Ecuyer's 63-bit multiplicative congruential
52 * generator:
53 * <PRE>
54 * seed := seed * A (mod M);
55 * </PRE> with <I>A</I> = 2307085864 and <I>M</I> = 2<SUP>63</SUP>-25. For
56 * further information, see P. L'Ecuyer, F. Blouin, and R. Couture, A search for
57 * good multiple recursive random number generators, <I>ACM Transactions on
58 * Modeling and Computer Simulation,</I> 3(2):87-98, April 1993.
59 *
60 * @author Alan Kaminsky
61 * @version 01-Mar-2008
62 */
63 public class Mcg1Random
64 extends Random {
65
66 @Serial
67 private static final long serialVersionUID = 1L;
68
69 // Hidden data members.
70 // Multiplicative congruential generator parameters.
71 private static final long A = 2307085864L;
72 private static final long M = 9223372036854775783L;
73
74 // Table of powers of A (mod M).
75 // powtable[i] = A^(2^i) (mod M), i = 0, 1, 2, ..., 62.
76 private static final long[] powtable = new long[]{
77 /* 0*/2307085864L,
78 /* 1*/ 5322645183868626496L,
79 /* 2*/ 983401115462215297L,
80 /* 3*/ 3556108090190705823L,
81 /* 4*/ 7990665143195102590L,
82 /* 5*/ 2110036525984475599L,
83 /* 6*/ 7043012601020815633L,
84 /* 7*/ 8705155707092105232L,
85 /* 8*/ 3648485552813098205L,
86 /* 9*/ 3168429798853819517L,
87 /*10*/ 7370936612916750461L,
88 /*11*/ 7860663018156131952L,
89 /*12*/ 3001105880121306407L,
90 /*13*/ 2701734581708584636L,
91 /*14*/ 44173215984149523L,
92 /*15*/ 4386281867185367357L,
93 /*16*/ 6179163218358095360L,
94 /*17*/ 7483044026478026567L,
95 /*18*/ 3475714592143337300L,
96 /*19*/ 1764426730688581302L,
97 /*20*/ 3750657437672096664L,
98 /*21*/ 622726075290379426L,
99 /*22*/ 5708473958970181660L,
100 /*23*/ 4021546582722653103L,
101 /*24*/ 2336213934427760687L,
102 /*25*/ 1250271094601288883L,
103 /*26*/ 3574383011208782094L,
104 /*27*/ 8396902035548884488L,
105 /*28*/ 8461483610275050157L,
106 /*29*/ 4570169555765982077L,
107 /*30*/ 8905831846701231221L,
108 /*31*/ 8735916407118983196L,
109 /*32*/ 2440495732904503112L,
110 /*33*/ 1885457269016286005L,
111 /*34*/ 4972446378304258072L,
112 /*35*/ 5086882142287647560L,
113 /*36*/ 7606891628733932672L,
114 /*37*/ 1492990033908793408L,
115 /*38*/ 9099993837175275499L,
116 /*39*/ 164616137930049276L,
117 /*40*/ 5117944347055477320L,
118 /*41*/ 3732738446422589684L,
119 /*42*/ 577797231373159603L,
120 /*43*/ 2884327325873197522L,
121 /*44*/ 4833803989390835826L,
122 /*45*/ 7647846260763424785L,
123 /*46*/ 4871120313232679781L,
124 /*47*/ 2522743552130321382L,
125 /*48*/ 2285147082121189109L,
126 /*49*/ 3702619298913044713L,
127 /*50*/ 7517285182136659617L,
128 /*51*/ 1501022168611987834L,
129 /*52*/ 4083684657803873370L,
130 /*53*/ 1174110446001111617L,
131 /*54*/ 82581059520186299L,
132 /*55*/ 1334190853588951475L,
133 /*56*/ 3130709730706025384L,
134 /*57*/ 8886205968707213290L,
135 /*58*/ 993283284549990895L,
136 /*59*/ 3258516944203296282L,
137 /*60*/ 4273233140749644635L,
138 /*61*/ 7682756089153477585L,
139 /*62*/ 8243539608199123644L,};
140
141 // Seed for this PRNG.
142 private long seed;
143
144 // 128 bytes of extra padding to avert cache interference.
145 private transient long p0, p1, p2, p3, p4, p5, p6, p7;
146 private transient long p8, p9, pa, pb, pc, pd, pe, pf;
147
148 // Exported constructors.
149 /**
150 * Construct a new PRNG with the given seed. The seed must not be 0.
151 *
152 * @param seed Seed.
153 * @exception IllegalArgumentException (unchecked exception) Thrown if
154 * <code>seed</code> = 0.
155 */
156 public Mcg1Random(long seed) {
157 setSeed(seed);
158 }
159
160 // Exported operations.
161 /**
162 * {@inheritDoc}
163 *
164 * Set this PRNG's seed. The seed must not be 0.
165 * @exception IllegalArgumentException (unchecked exception) Thrown if
166 * <code>seed</code> = 0.
167 */
168 public void setSeed(long seed) {
169 if (seed == 0L) {
170 throw new IllegalArgumentException("Mcg1Random.setSeed(): seed = 0 illegal");
171 }
172
173 // Make sure seed is nonnegative.
174 this.seed = seed & 0x7FFFFFFFFFFFFFFFL;
175 }
176
177 // Hidden operations.
178 /**
179 * Return the next 64-bit pseudorandom value in this PRNG's sequence.
180 *
181 * @return Pseudorandom value.
182 */
183 protected long next() {
184 // Multiply seed (a 64-bit number) by A (a 32-bit number) yielding x (a
185 // 96-bit number). Bits 63-0 of x are stored in x_63_0. Bits 95-32 of x
186 // are stored in x_95_32. Note that these overlap.
187 long tmp_63_0 = (seed & 0x00000000FFFFFFFFL) * A;
188 long tmp_95_32 = (seed >>> 32) * A;
189 long x_63_0 = (tmp_95_32 << 32) + tmp_63_0;
190 long x_95_32 = tmp_95_32 + (tmp_63_0 >>> 32);
191
192 // Compute x mod M, where M = 2^63-25. For the algorithm, see the
193 // Handbook of Applied Cryptography, Section 14.3.4.
194 // q = int (x / 2^63)
195 long q = x_95_32 >>> 31;
196
197 // r = x mod 2^63
198 long r = x_63_0 & 0x7FFFFFFFFFFFFFFFL;
199
200 // r = r + (q * 25) mod 2^63
201 r += q * 25L;
202
203 // If there was a carry into the high-order bit of r, or if r >= M,
204 // subtract M.
205 if (r < 0L || r >= M) {
206 r -= M;
207 }
208
209 // r = x mod M becomes the new seed.
210 seed = r;
211
212 // Since seed is only in the range 0 .. 2^63 - 1, left-shift yielding a
213 // number in the range -2^63 .. 2^63 - 1.
214 return seed << 1;
215 }
216
217 /**
218 * {@inheritDoc}
219 *
220 * Return the 64-bit pseudorandom value the given number of positions ahead
221 * in this PRNG's sequence.
222 */
223 protected long next(long skip) {
224 // Compute seed * A^skip (mod M).
225 int i = 0;
226 while (skip != 0L) {
227 if ((skip & 1L) != 0L) {
228 seed = modMultiply(powtable[i], seed);
229 }
230 skip >>>= 1;
231 ++i;
232 }
233
234 // Since seed is only in the range 0 .. 2^63 - 1, left-shift yielding a
235 // number in the range -2^63 .. 2^63 - 1.
236 return seed << 1;
237 }
238
239 /**
240 * Returns a * b (mod M). a and b are assumed to be in the range 0 ..
241 * 2^63-1.
242 *
243 * @param a First number to multiply.
244 * @param b Second number to multiply.
245 *
246 * @return a * b (mod M).
247 */
248 private static long modMultiply(long a,
249 long b) {
250 // Let a = s*2^32 + t, b = u*2^32 + v, where s, t, u, and v are 32-bit
251 // numbers. Form the four 64-bit products tv, sv, tu, and su. Add these
252 // up in the appropriate combinations to get x = a * b, where x =
253 // x_127_64*2^64 + x_63_0:
254 // +--------+--------+
255 // | s | t | = a
256 // +--------+--------+
257 // +--------+--------+
258 // * | u | v | = b
259 // +--------+--------+
260 // ----------------------------------------
261 // +-----------------+
262 // | tv |
263 // +-----------------+
264 // +-----------------+
265 // | sv |
266 // +-----------------+
267 // +-----------------+
268 // | tu |
269 // +-----------------+
270 // +-----------------+
271 // + | su |
272 // +-----------------+
273 // ----------------------------------------
274 // +-----------------+-----------------+
275 // | x_127_64 | x_63_0 | = x = a * b
276 // +-----------------+-----------------+
277 long s = a >>> 32;
278 long t = a & 0xFFFFFFFFL;
279 long u = b >>> 32;
280 long v = b & 0xFFFFFFFFL;
281 long tv = t * v;
282 long sv = s * v;
283 long tu = t * u;
284 long su = s * u;
285 long tmp = (tv >>> 32) + (sv & 0xFFFFFFFFL) + (tu & 0xFFFFFFFFL);
286 long x_63_0 = (tv & 0xFFFFFFFFL) + (tmp << 32);
287 long x_127_64 = (tmp >>> 32) + (sv >>> 32) + (tu >>> 32) + su;
288
289 // Compute x mod M, where M = 2^63-25. For the algorithm, see the
290 // Handbook of Applied Cryptography, Section 14.3.4.
291 // q = int (x / 2^63)
292 long q = (x_127_64 << 1) | (x_63_0 >>> 63);
293
294 // r = x mod 2^63
295 long r = x_63_0 & 0x7FFFFFFFFFFFFFFFL;
296
297 while (q > 0L) {
298 // qc = q * 25
299 // Multiply q (a 64-bit number) by c (a 32-bit number) yielding qc
300 // (a 96-bit number). Bits 63-0 of qc are stored in qc_63_0. Bits
301 // 95-32 of qc are stored in qc_95_32. Note that these overlap.
302 long tmp_63_0 = (q & 0xFFFFFFFFL) * 25L;
303 long tmp_95_32 = (q >>> 32) * 25L;
304 long qc_63_0 = (tmp_95_32 << 32) + tmp_63_0;
305 long qc_95_32 = tmp_95_32 + (tmp_63_0 >>> 32);
306
307 // q = int (qc / 2^63)
308 q = qc_95_32 >>> 31;
309
310 // r = r + (qc mod 2^63)
311 r += qc_63_0 & 0x7FFFFFFFFFFFFFFFL;
312
313 // If there was a carry into the high-order bit of r, or if r >= M,
314 // subtract M.
315 if (r < 0L || r >= M) {
316 r -= M;
317 }
318 }
319
320 // Return r = x mod M.
321 return r;
322 }
323
324 }