15th October 2014
Many applications depend on streams of pseudo random numbers from a pseudo random number generator (PRNG). For example, Monte Carlo based simulations and machine learning algorithms which include some random component require the output of a PRNG. Often these PRNGs are seeded explicitly so that we can reproduce identical output for a particular set of parameters or input data set.
If we want to accelerate the processing of the these algorithms through parallelisation and we also wish to preserve the repeatability of the output, we require a mechanism for reproducing the same stream of random numbers that would have been encountered by the serial version. This is often achieved through a technique known as PRNG “leapfrog” or “skipahead”. Simply running the generator through the required number of random elements is too costly, with approx complexity. An LCG PRNG can be “leapfrogged” efficiently using a recurrence relation . The computationally expensive component of this calculation is the calculation of the modular exponent of the LCG parameter .
However, integer powers can be calculated very efficiently using binary exponentiation, which is the same technique used in elliptic-curve cryptography. This reduces the complexity to very approximately . As the skipahead relation uses a modular exponent, we can use modular exponentiation which prevents the intermediate results growing very large and would require arbitrary precision integers.
In fact, the modular exponent is required frequently enough that the Python pow() function takes an optional third parameter to specify the modulus of the output.
It was recently necessary for me to use this technique in a Fortran application, without big integer support.
To solve this I came up with the implementation below. The modipow function calculates the integer exponent of the parameter base, raised to the power of ex, modulus m. It works for the full range of signed 64 bit integers.
function modipow(base,ex,m) integer(8) :: modipow integer(8),intent(in) :: base,ex,m integer(8) :: res,base_t,ex_t res = 1 base_t = mod(base,m) ex_t = ex do while (ex_t.gt.0) if (AND(ex_t,1)) then res = res * base_t res = mod(res,m) endif ex_t = ishft(ex_t,-1) base_t = base_t * base_t base_t = mod(base_t,m) enddo modipow = mod(res,m) end function modipow