64 bit Modular Exponentiation in Fortran

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.

Using PRNGs in Parallel

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 O(n) complexity. An LCG PRNG can be “leapfrogged” efficiently using a recurrence relation [1]. The computationally expensive component of this calculation is the calculation of the modular exponent of the LCG parameter a.

Fast Calculation Of Integer Powers

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 O(log(n)). 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.

Fortran Native Implementation of Modular Exponentiation

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

[1] http://nayuki.eigenstate.org/page/fast-skipping-in-a-linear-congruential-generator

Share This

Tweet this!   Share on LinkedIn   Share on Facebook

Leave a Reply

Your email address will not be published. Required fields are marked *