pfenerty
01-05-2006, 05:43 PM
Hello,
For the purpose of Monte Carlo probabilistic dart throwing (2 coordinates per dart throw, 1 rand per coordinate), I have a (combined, multiplicative) LCG that is parallelified, for embarrasingly parallel, distributed dart throwing, according to an algorithm that can be summarized as:
n dartServers each generate the same pseudo-random sequence, from the same seed, and then leapfrog one another in selecting pairs of rands from that sequence in order to generate, jointly, the complete sequence of pseudo-random dart throws. The servers accomplish this leapfrogging by replacing the sequential multiplier ('a') in the LCG relation:
X_sub_n+1 = ( a * X_sub_n ) mod m
with a leapfrog multiplier ('A') according to:
X_sub_n+k = ( A * X_sub_n ) mod m
where:
A = a^k
... with 'k' given by (2 * numberOfServers), since rands are needed in pairs.
So each dart throw requires, first, the next leapfrogged rand:
// rand for x coordinate
thisRand = (leapfrogMultiplier * previousRand) mod m
followed by its next sequential rand:
// rand for y coordinate
thisRand = (sequentialMultiplier * previousRand) mod m
This was seen to work correctly (i.e., the same sequence of rands produced by a single server is also produced jointly, by multiple servers). Some time later, I learned of E. L. Schrage's marvelous method for avoiding intermediate multiplication overflows, and sought to adapt my original algorithm accordingly.
But the typical implementation for Schrage's method cannot be adapted to this leapfrogging technique, it seems to me, because the "compensation" for 'k', that transforms 'a' into 'A', appears multiple times, in multiple terms, in Schrage's version (i.e., below, in 'A', 'Q', and 'R'):
// Schrage's next leapfrogged rand
x = A * (x - ( x / Q) * Q ) - R * (x / Q)
where:
Q = m / A
R = m mod A
This approach was implemented anyway, and confirmed NOT to work correctly, except when comparing the single server case, using Schrage's method, to any of the cases without Schrage (i.e., for the Schrage-version degenerate case A == a), which eliminates simple coding errors.
Is there a solution for this ... a leapfroggable implementation of Schrage's method? I'm afraid I don't have the background to derive it myself, if it exists.
Thanks so much!
Paul Fenerty
For the purpose of Monte Carlo probabilistic dart throwing (2 coordinates per dart throw, 1 rand per coordinate), I have a (combined, multiplicative) LCG that is parallelified, for embarrasingly parallel, distributed dart throwing, according to an algorithm that can be summarized as:
n dartServers each generate the same pseudo-random sequence, from the same seed, and then leapfrog one another in selecting pairs of rands from that sequence in order to generate, jointly, the complete sequence of pseudo-random dart throws. The servers accomplish this leapfrogging by replacing the sequential multiplier ('a') in the LCG relation:
X_sub_n+1 = ( a * X_sub_n ) mod m
with a leapfrog multiplier ('A') according to:
X_sub_n+k = ( A * X_sub_n ) mod m
where:
A = a^k
... with 'k' given by (2 * numberOfServers), since rands are needed in pairs.
So each dart throw requires, first, the next leapfrogged rand:
// rand for x coordinate
thisRand = (leapfrogMultiplier * previousRand) mod m
followed by its next sequential rand:
// rand for y coordinate
thisRand = (sequentialMultiplier * previousRand) mod m
This was seen to work correctly (i.e., the same sequence of rands produced by a single server is also produced jointly, by multiple servers). Some time later, I learned of E. L. Schrage's marvelous method for avoiding intermediate multiplication overflows, and sought to adapt my original algorithm accordingly.
But the typical implementation for Schrage's method cannot be adapted to this leapfrogging technique, it seems to me, because the "compensation" for 'k', that transforms 'a' into 'A', appears multiple times, in multiple terms, in Schrage's version (i.e., below, in 'A', 'Q', and 'R'):
// Schrage's next leapfrogged rand
x = A * (x - ( x / Q) * Q ) - R * (x / Q)
where:
Q = m / A
R = m mod A
This approach was implemented anyway, and confirmed NOT to work correctly, except when comparing the single server case, using Schrage's method, to any of the cases without Schrage (i.e., for the Schrage-version degenerate case A == a), which eliminates simple coding errors.
Is there a solution for this ... a leapfroggable implementation of Schrage's method? I'm afraid I don't have the background to derive it myself, if it exists.
Thanks so much!
Paul Fenerty