diff --git a/src/generator/congruential32.spec.ts b/src/generator/congruential32.spec.ts index 2344fbb..a6911ab 100644 --- a/src/generator/congruential32.spec.ts +++ b/src/generator/congruential32.spec.ts @@ -128,12 +128,43 @@ describe('congruential32', () => { ].map((v) => v | 0), ); }); + it('Should produce the right sequence after jump for seed=42', () => { + const g = congruential32(42); + g.jump(); + const data = []; + for (let idx = 0; idx !== 100; ++idx) { + const v = g.next(); + data.push(v); + } + // The following values correspond to the values that would have been extracted from `g` + // if instead of `g.jump()` we called `for (let i = 0; i !== 2**16; ++i) { g.next(); }`. + expect(data).toEqual( + [ + 2465055980, 3883212298, 4021781768, 738993281, 3215657551, 984983963, 2896152994, 1422049114, 1152760575, + 2899723742, 2082487871, 398492622, 76255798, 2215279246, 3044661771, 2239891048, 2928515811, 2990029361, + 2305593392, 1500586776, 3294838601, 2335870109, 833774352, 7718367, 285290200, 1060774088, 1896682433, + 2661497019, 2001778250, 3610235348, 1510287817, 570761941, 510977537, 3514690393, 842682139, 3499042025, + 2712543904, 481540341, 3898195024, 2670247668, 1239841767, 241610001, 850833705, 1071541028, 1993336276, + 3560123125, 3430856005, 4113308155, 3182374902, 4242096142, 663756062, 2557311941, 599938835, 1467766910, + 2297401161, 1634419517, 977439997, 2264808165, 635274734, 487214458, 3255669425, 2273302381, 3079637125, + 2702596895, 485245105, 3288217115, 461445586, 1798588098, 846485922, 974131284, 4134376994, 2028790940, + 244515625, 163930541, 1877278150, 1188919153, 1605435649, 1872886775, 3162597328, 2231559874, 444075616, + 1187015306, 1287526668, 774513201, 3912888716, 285946065, 2523757377, 276618804, 754296248, 2745085974, + 1322021033, 4113569277, 905852707, 4114890895, 3409788845, 894566056, 4198542707, 756173224, 513315557, + 699881382, + ].map((v) => v | 0), + ); + }); it('Should return the same sequence given same seeds', () => fc.assert(p.sameSeedSameSequences(congruential32))); it('Should return the same sequence when built from state', () => fc.assert(p.clonedFromStateSameSequences(congruential32, congruential32FromState))); it('Should return the same sequence if called twice', () => fc.assert(p.sameSequencesIfCallTwice(congruential32))); it('Should generate values between -2**31 and 2**31 -1', () => fc.assert(p.valuesInRange(congruential32))); + it('Should not depend on ordering between jump and next', () => fc.assert(p.noOrderNextJump(congruential32))); it('Should impact itself with next', () => fc.assert(p.changeSelfWithNext(congruential32))); + it('Should impact itself with jump', () => fc.assert(p.changeSelfWithJump(congruential32))); it('Should not impact clones when impacting itself on next', () => fc.assert(p.noChangeOnClonedWithNext(congruential32))); + it('Should not impact clones when impacting itself on jump', () => + fc.assert(p.noChangeOnClonedWithJump(congruential32))); }); diff --git a/src/generator/congruential32.ts b/src/generator/congruential32.ts index b4fdc3f..7a124d0 100644 --- a/src/generator/congruential32.ts +++ b/src/generator/congruential32.ts @@ -1,4 +1,4 @@ -import type { RandomGenerator } from '../types/RandomGenerator'; +import type { JumpableRandomGenerator } from '../types/JumpableRandomGenerator'; // Inspired from java.util.Random implementation // http://grepcode.com/file/repository.grepcode.com/java/root/jdk/openjdk/6-b14/java/util/Random.java#Random.next%28int%29 @@ -8,7 +8,19 @@ const INCREMENT: number = 0x00269ec3; const MASK: number = 0xffffffff; const MASK_2: number = (1 << 31) - 1; -class LinearCongruential32 implements RandomGenerator { +// Jump constants: equivalent to 2^16 calls to next() (= 3 * 2^16 LCG steps). +// Precomputed using the LCG jump-ahead algorithm from: +// F. Brown, "Random Number Generation with Arbitrary Strides", Trans. Am. Nucl. Soc. (Nov. 1994) +// Also described at: https://www.pcg-random.org/posts/bounded-rands.html +// +// For an LCG s_{n+1} = (a * s_n + c) mod m, jumping by k steps yields: +// s_{n+k} = (A * s_n + C) mod m +// where A = a^k mod m and C = c * (a^{k-1} + ... + 1) mod m. +// With a=0x000343fd, c=0x00269ec3, m=2^32, and k=3*2^16=196608: +const JUMP_MULTIPLIER: number = 0x76dc0001; +const JUMP_INCREMENT: number = 0x369b0000; + +class LinearCongruential32 implements JumpableRandomGenerator { constructor(private seed: number) {} clone(): LinearCongruential32 { return new LinearCongruential32(this.seed); @@ -27,6 +39,12 @@ class LinearCongruential32 implements RandomGenerator { const vnext = v3 + ((v2 + (v1 << 15)) << 15); return vnext | 0; } + jump(): void { + // equivalent to 2^16 calls to next() + // can be used to generate 2^16 non-overlapping subsequences for the full 2^32 period + // Math.imul is required because seed * JUMP_MULTIPLIER can exceed 2^53 + this.seed = (Math.imul(this.seed, JUMP_MULTIPLIER) + JUMP_INCREMENT) & MASK; + } getState(): readonly number[] { return [this.seed]; } @@ -39,7 +57,7 @@ function computeValueFromNextSeed(nextseed: number) { return (nextseed & MASK_2) >> 16; } -export function congruential32FromState(state: readonly number[]): RandomGenerator { +export function congruential32FromState(state: readonly number[]): JumpableRandomGenerator { const valid = state.length === 1; if (!valid) { throw new Error('The state must have been produced by a congruential32 RandomGenerator'); @@ -47,6 +65,6 @@ export function congruential32FromState(state: readonly number[]): RandomGenerat return new LinearCongruential32(state[0]); } -export function congruential32(seed: number): RandomGenerator { +export function congruential32(seed: number): JumpableRandomGenerator { return new LinearCongruential32(seed); }