1 /** 2 * Implementation of the xoroshiro128+ uniform random number generator. 3 * 4 * Authors: 5 * $(LINK2 http://braingam.es/, Joseph Rushton Wakeling) 6 * 7 * Copyright: 8 * Written in 2016 by Joseph Rushton Wakeling. 9 * 10 * License: 11 * $(LINK2 https://creativecommons.org/publicdomain/zero/1.0/legalcode, Creative Commons CC0) 12 * (public domain) 13 */ 14 module dxorshift.xoroshiro128plus; 15 16 /** 17 * The xoroshiro128+ uniform random number generator is the 18 * fastest full-period generator passing the BigCrush test 19 * suite without systematic failures. 20 * 21 * Due to the relatively short period (2 ^^ 128 - 1) it is 22 * acceptable only for applications with a mild amount of 23 * parallelism; for applications requiring many parallel 24 * random sequences, `Xorshift1024star` is recommended 25 * instead. 26 * 27 * Besides passing BigCrush, this generator passes the 28 * PractRand test suite up to (and including) 16TB, with 29 * the exception of binary rank tests, which fail due to 30 * the lowest bit being a linear feedback shift register 31 * (LFSR). All other bits pass all tests. 32 * 33 * Credits: This code is ported from the public-domain 34 * reference implementation by David Blackman 35 * and Sebastiano Vigna, available online at 36 * $(LINK http://xorshift.di.unimi.it/xoroshiro128plus.c) 37 */ 38 struct Xoroshiro128plus 39 { 40 public: 41 /// Marks this range as a uniform random number generator 42 enum bool isUniformRandom = true; 43 44 /// Smallest generated value (0) 45 enum ulong min = ulong.min; 46 47 /// Largest generated value 48 enum ulong max = ulong.max; 49 50 /* Copy-by-value is disabled to avoid unintended 51 * duplication of random sequences; use the `dup` 52 * property if you really wish to copy the state 53 * of the RNG. 54 */ 55 @disable this(this); 56 57 // RNG can only be initialized with a seed 58 @disable this(); 59 60 /** 61 * Constructor (RNG instances can only be initialized 62 * with a specified seed). 63 */ 64 this(ulong s) @nogc @safe nothrow pure 65 { 66 this.seed(s); 67 } 68 69 /// ditto 70 this(ulong s0, ulong s1) @nogc @safe nothrow pure 71 { 72 this.seed(s0, s1); 73 } 74 75 /// ditto 76 this(ulong[2] s) @nogc @safe nothrow pure 77 { 78 this.seed(s); 79 } 80 81 /// ditto 82 this(R)(R range) 83 if (isInputRange!R && is(Unqual!(ElementType!R) == ulong)) 84 { 85 this.seed(range); 86 } 87 88 /// Range primitives 89 enum bool empty = false; 90 91 /// ditto 92 ulong front() @nogc @property @safe const nothrow pure 93 { 94 return this.state[0] + this.state[1]; 95 } 96 97 /// ditto 98 void popFront() @nogc @safe nothrow pure 99 in 100 { 101 assert(this.state != [0uL, 0uL]); 102 } 103 body 104 { 105 immutable ulong s1 = this.state[1] ^ this.state[0]; 106 this.state[0] = rotateLeft(this.state[0], 55) ^ s1 ^ (s1 << 14); 107 this.state[1] = rotateLeft(s1, 36); 108 } 109 110 /** 111 * Jump function, equivalent to 2 ^^ 64 calls to 112 * `popFront()`; can be used to generate 2 ^^ 64 113 * non-overlapping subsequences for 114 * parallel computation. 115 */ 116 void jump() @nogc @safe nothrow pure 117 { 118 enum ulong[2] jump_ = [0xbeac0467eba5facb, 0xd86b048b86aa9922]; 119 enum jumpSize = jump_.sizeof / (*(jump_.ptr)).sizeof; 120 ulong s0 = 0; 121 ulong s1 = 0; 122 123 for (size_t i = 0; i < jumpSize; ++i) 124 { 125 for (int b = 0; b < 64; ++b) 126 { 127 if (jump_[i] & 1uL << b) 128 { 129 s0 ^= this.state[0]; 130 s1 ^= this.state[1]; 131 } 132 popFront(); 133 } 134 } 135 136 this.state[0] = s0; 137 this.state[1] = s1; 138 } 139 140 /** 141 * Provides a copy of this RNG instance with identical 142 * internal state. 143 * 144 * This property is provided in preference to `save` 145 * so as to allow the user to duplicate RNG state 146 * when explicitly desired, but without risking 147 * unintended copies by functions that save forward 148 * ranges provided as input. 149 */ 150 typeof(this) dup() @nogc @property @safe nothrow pure 151 { 152 return typeof(this)(this); 153 } 154 155 /// (Re)seeds the generator. 156 void seed(ulong s) @nogc @safe nothrow pure 157 { 158 import dxorshift.splitmix64 : SplitMix64; 159 this.seed(SplitMix64(s)); 160 } 161 162 /// ditto 163 void seed(ulong s0, ulong s1) @nogc @safe nothrow pure 164 in 165 { 166 // seeds are not both 0 167 assert(!(!s0 && !s1)); 168 } 169 body 170 { 171 this.state[0] = s0; 172 this.state[1] = s1; 173 popFront(); 174 } 175 176 /// ditto 177 void seed(ulong[2] s) @nogc @safe nothrow pure 178 in 179 { 180 assert(!(!s[0] && !s[1])); 181 } 182 body 183 { 184 seed(s[0], s[1]); 185 } 186 187 /// ditto 188 void seed(R)(R range) 189 if (isInputRange!R && is(Unqual!(ElementType!R) == ulong)) 190 { 191 foreach (ref s; this.state) 192 { 193 import std.range.primitives : empty, front, popFront; 194 assert(!range.empty, "Insufficient elements to populate RNG state!"); 195 s = range.front; 196 range.popFront(); 197 } 198 assert(this.state != typeof(this.state).init, "Seed elements cannot all be zero!"); 199 this.popFront(); 200 } 201 202 private: 203 // 128 bits of state 204 ulong[2] state; 205 206 // Helper constructor used to implement `dup` 207 this(ref typeof(this) that) @nogc @safe nothrow pure 208 { 209 this.state[] = that.state[]; 210 } 211 212 // Simulated rotate operation used in `popFront()` 213 static ulong rotateLeft(ulong x, int k) @nogc @safe nothrow pure 214 in 215 { 216 assert(0 <= k); 217 assert(k <= 64); 218 } 219 body 220 { 221 return (x << k) | (x >> (64 - k)); 222 } 223 224 import std.range.primitives : ElementType, isInputRange; 225 import std.traits : Unqual; 226 } 227 228 /// 229 unittest 230 { 231 import std.array : array; 232 import std.random : isUniformRNG, randomSample, uniform; 233 import std.range : iota, take; 234 import dxorshift.xoroshiro128plus; 235 236 // xoroshiro128+ generators must be initialized 237 // with a specified seed 238 auto gen = Xoroshiro128plus(123456); 239 240 // verify it is indeed a uniform RNG as defined 241 // in the standard library, whether accessed 242 // directly or via a pointer 243 static assert(isUniformRNG!(typeof(gen))); 244 static assert(isUniformRNG!(typeof(&gen))); 245 246 // since the postblit is disabled, we must 247 // pass a pointer to any functionality that 248 // would otherwise copy the RNG by value 249 assert((&gen).take(2).array == [14854895758870614632uL, 250 2102156639392820999uL]); 251 252 // this means, of course, that we must guarantee 253 // the lifetime of the pointer is valid for the 254 // lifetime of any functionality that uses it 255 auto sample = iota(100).randomSample(10, &gen).array; 256 257 // however, we can pass the RNG as-is to any 258 // functionality that takes it by ref and does 259 // not try to copy it by value 260 auto val = uniform!"(]"(3.5, 4.0, gen); 261 262 // in circumstances where we really want to 263 // copy the RNG state, we can use `dup` 264 auto gen2 = gen.dup; 265 assert((&gen).take(9).array == (&gen2).take(9).array); 266 } 267 268 unittest 269 { 270 import std.array : array; 271 import std.random : isUniformRNG, isSeedable; 272 import std.range : take; 273 274 import dxorshift.splitmix64 : SplitMix64; 275 276 static assert(isUniformRNG!Xoroshiro128plus); 277 278 static assert(isSeedable!Xoroshiro128plus); 279 static assert(isSeedable!(Xoroshiro128plus, ulong[2])); 280 static assert(isSeedable!(Xoroshiro128plus, ulong[])); 281 static assert(isSeedable!(Xoroshiro128plus, ulong)); 282 static assert(isSeedable!(Xoroshiro128plus, SplitMix64)); 283 284 // output comparisons to reference implementation, 285 // using constructor, seeding, and duplication 286 { 287 auto gen = Xoroshiro128plus(123, 456); 288 assert((&gen).take(10).array == [4431571926312075699uL, 16834163345174162131uL, 289 4468099366319113814uL, 167286530559998105uL, 290 18053350147704165166uL, 9927833068668020801uL, 291 14561249726909464726uL, 316732314664215799uL, 292 8800682043873892537uL, 8955312909536390945uL]); 293 294 gen.seed(123, 456); 295 auto gen2 = gen.dup; 296 assert((&gen).take(10).array == [4431571926312075699uL, 16834163345174162131uL, 297 4468099366319113814uL, 167286530559998105uL, 298 18053350147704165166uL, 9927833068668020801uL, 299 14561249726909464726uL, 316732314664215799uL, 300 8800682043873892537uL, 8955312909536390945uL]); 301 302 assert((&gen2).take(10).array == [4431571926312075699uL, 16834163345174162131uL, 303 4468099366319113814uL, 167286530559998105uL, 304 18053350147704165166uL, 9927833068668020801uL, 305 14561249726909464726uL, 316732314664215799uL, 306 8800682043873892537uL, 8955312909536390945uL]); 307 } 308 309 { 310 ulong[2] s = [12345uL, 67890uL]; 311 312 auto gen = Xoroshiro128plus(s); 313 assert((&gen).take(10).array == [2059148541540170003uL, 6156794878792115187uL, 314 559523256690861310uL, 15314907387785043984uL, 315 4915457426679953335uL, 5462571845969584332uL, 316 9658602537074702831uL, 17979359875003347608uL, 317 5174518315773110499uL, 2532010971873184518uL]); 318 gen.seed(s); 319 auto gen2 = gen.dup; 320 assert((&gen).take(10).array == [2059148541540170003uL, 6156794878792115187uL, 321 559523256690861310uL, 15314907387785043984uL, 322 4915457426679953335uL, 5462571845969584332uL, 323 9658602537074702831uL, 17979359875003347608uL, 324 5174518315773110499uL, 2532010971873184518uL]); 325 326 assert((&gen2).take(10).array == [2059148541540170003uL, 6156794878792115187uL, 327 559523256690861310uL, 15314907387785043984uL, 328 4915457426679953335uL, 5462571845969584332uL, 329 9658602537074702831uL, 17979359875003347608uL, 330 5174518315773110499uL, 2532010971873184518uL]); 331 } 332 333 { 334 auto s = 123456; 335 336 auto gen = Xoroshiro128plus(s); 337 assert((&gen).take(10).array == [14854895758870614632uL, 2102156639392820999uL, 338 13092495043793465900uL, 8397221095866455920uL, 339 6262852887298792196uL, 16202237309782713452uL, 340 14544835201639844962uL, 7120381903468495472uL, 341 4724551740662753335uL, 3230748688607015409uL]); 342 343 gen.seed(s); 344 auto gen2 = gen.dup; 345 assert((&gen).take(10).array == [14854895758870614632uL, 2102156639392820999uL, 346 13092495043793465900uL, 8397221095866455920uL, 347 6262852887298792196uL, 16202237309782713452uL, 348 14544835201639844962uL, 7120381903468495472uL, 349 4724551740662753335uL, 3230748688607015409uL]); 350 351 assert((&gen2).take(10).array == [14854895758870614632uL, 2102156639392820999uL, 352 13092495043793465900uL, 8397221095866455920uL, 353 6262852887298792196uL, 16202237309782713452uL, 354 14544835201639844962uL, 7120381903468495472uL, 355 4724551740662753335uL, 3230748688607015409uL]); 356 } 357 358 { 359 auto s = SplitMix64(123456); 360 361 auto gen = Xoroshiro128plus(&s); 362 assert((&gen).take(10).array == [14854895758870614632uL, 2102156639392820999uL, 363 13092495043793465900uL, 8397221095866455920uL, 364 6262852887298792196uL, 16202237309782713452uL, 365 14544835201639844962uL, 7120381903468495472uL, 366 4724551740662753335uL, 3230748688607015409uL]); 367 368 s.seed(123456); 369 gen.seed(&s); 370 auto gen2 = gen.dup; 371 assert((&gen).take(10).array == [14854895758870614632uL, 2102156639392820999uL, 372 13092495043793465900uL, 8397221095866455920uL, 373 6262852887298792196uL, 16202237309782713452uL, 374 14544835201639844962uL, 7120381903468495472uL, 375 4724551740662753335uL, 3230748688607015409uL]); 376 377 assert((&gen2).take(10).array == [14854895758870614632uL, 2102156639392820999uL, 378 13092495043793465900uL, 8397221095866455920uL, 379 6262852887298792196uL, 16202237309782713452uL, 380 14544835201639844962uL, 7120381903468495472uL, 381 4724551740662753335uL, 3230748688607015409uL]); 382 } 383 384 // compare jump to reference implementation, 385 // using constructor, seeding and duplication 386 { 387 auto gen = Xoroshiro128plus(123, 456); 388 gen.jump(); 389 assert((&gen).take(10).array == [3109680772672824672uL, 11190329315615627403uL, 390 2415690012231644097uL, 2347094600162878539uL, 391 18099586205610688946uL, 7375268959557732117uL, 392 12413671816612458655uL, 14565394836119542025uL, 393 5936088160154203578uL, 12124177863926731024uL]); 394 395 gen.seed(123, 456); 396 auto gen2 = gen.dup; 397 gen.jump(); 398 assert((&gen).take(10).array == [3109680772672824672uL, 11190329315615627403uL, 399 2415690012231644097uL, 2347094600162878539uL, 400 18099586205610688946uL, 7375268959557732117uL, 401 12413671816612458655uL, 14565394836119542025uL, 402 5936088160154203578uL, 12124177863926731024uL]); 403 404 gen2.jump(); 405 assert((&gen2).take(10).array == [3109680772672824672uL, 11190329315615627403uL, 406 2415690012231644097uL, 2347094600162878539uL, 407 18099586205610688946uL, 7375268959557732117uL, 408 12413671816612458655uL, 14565394836119542025uL, 409 5936088160154203578uL, 12124177863926731024uL]); 410 } 411 }