Super fast sampling from a tricky random distribution

Question

This is a micro-optimization challenge. The task is to sample from the maximally skewed stable distribution F(x;1,-1,Pi/2,0). See Table 1 of http://arxiv.org/pdf/0908.3961v2.pdf . It is easiest to describe how to do this using the following C code which uses the Mersenne Twister RNG C code from here .

#include <stdio.h>
#include <stdlib.h>
#include "mtwist.h"
#include <math.h>


int main(void) {
   int i;
   volatile double x;
   mt_seed();
   double u1;
   double u2;
   double w1;
   double w2;
   for(i = 0; i < 100000000; ++i) {
     u1 = mt_drand();
     u2 = mt_drand();
     w1 = M_PI*(u1-1/2.0);
     w2 = -log(u2);
     x = tan(w1)*(M_PI_2-w1)+log(w2*cos(w1)/(M_PI_2-w1));
   }
   return EXIT_SUCCESS;
}

I compile with

gcc -Wall -O3 random.c mtwist-1.5/mtwist.c -o random -lm

The running time is 20 seconds taking 5,000,000 iterations per second.

There are two obvious ways of making this fast. The first is to choose a very fast uniform random number generator. If you choose to use anything but the Mersenne twister then you must run it through the Diehard tests and check the P-values provided. Use the code here http://www.stat.fsu.edu/pub/diehard/ to do this. Your uniform random number generator must pass at least 15 of these tests.

The second way is to find fast micro-optimizations for the transcendental functions, perhaps using some CPU specific features. For example in my case these features of the AMD FX-8350 CPU.

Rules The code must maintain at least 32 bits of accuracy throughout. The code should be compilable on ubuntu using easy to install free software only. Please provide full compile and run instructions. You should also provide a diehard log for any non Mersenne twister RNGs as was done for Build a random number generator that passes the Diehard tests . Note the very simple RNG that passes 15 tests in that link.

Scores You can report your score as the time taken by the code above divided by the time taken by your code on your own machine. For running your code you can use any compiler and compilation options you like. To test the code above you should use gcc 4.8 or later with -O3. If any submission take less than one second I will increase N by factors of 10 until the fastest submission takes longer than 1 second.

Test machine My test machine was a standard 8GB RAM ubuntu install on an AMD FX-8350 Processor. Your machine can be anything that supports gcc 4.8 or later.

The Intel optimization manual says

If there is no critical need to evaluate the transcendental functions using the extended precision of 80 bits, applications should consider an alternate, software-based approach, such as a look-up-table-based algorithm using interpolation techniques. It is possible to improve transcendental performance with these techniques by choosing the desired numeric precision and the size of the look-up table, and by taking advantage of the parallelism of the SSE and the SSE2 instructions.

According to this very helpful table,fcos has latency 154 and fptan has latency 166-231 on my AMD FX 8350.

Answer 1

GMPY2

Abuse of GMP. @OP the slowest part of the code is printing the numbers, if leaving them in memory is fine the code would be much much faster.

import math
import random
import gmpy2
from gmpy2 import mpz, mpfr

gmpy2.get_context().precision = 11 # 32 bits of accuracy

N = 100000000
random_state = gmpy2.random_state()
pi = math.pi
half_pi = pi / 2
size = 0
for i in range(N):
    u1 = gmpy2.mpfr_random(random_state)
    u2 = gmpy2.mpfr_random(random_state) 
    w1 = gmpy2.mul(pi, gmpy2.sub(u1,1/2))
    w2 = -gmpy2.log(u2)
    a = gmpy2.mul(gmpy2.tan(w1), gmpy2.sub(half_pi, w1))
    b = gmpy2.div(gmpy2.log(gmpy2.mul(w2, gmpy2.cos(w1))), gmpy2.sub(half_pi, w1))
    x = gmpy2.add(a, b)

Answer 2

Nimrod (Score = 1.71, 4.8 vs. 8.2 seconds)

Observe that x->[[cos x, -sin x], [sin x, cos x]] is a group isomorphism between [0, 2*pi) (with respect to addition modulo 2*pi) and the rotation matrices (as a subgroup of SL(2, R)).

Because w1 is effectively of the form pi*k/2^32, where k is a signed 32-bit integer, we can split k in two 16-bit words, h and l and always represent w1 as the sum of pi*h/2^16+pi*l/2^32 and tabulate rotation matrices for both addends in 2*sizeof(double)*2^16 bytes each. We can then derive sin(w1) and cos(w1) from the product of the two rotation matrices per the isomorphism above (and tan(w1) as sin(w1)/cos(w1)).

The code below is in Nimrod, but should be directly portable to C, with the necessary renaming to integer and float types (float64 in Nimrod corresponds to double in C, and int corresponds to intptr_t, i.e. 64-bit on most modern architectures) and essentially no change in runtime. It currently does not print anything, but uses a volatile variable as described in my earlier comment to avoid the code being optimized away.

Further speedups may come from looking at how to make the computation of the natural logarithm cheaper (and, of course, replace the Mersenne Twister with a different RNG, but that's not very interesting).

Compile with nimrod cc -d:release random.nim.

Compile with -d:validate instead to check floating point errors with respect to the normal sin/cos calculations.

Update: Added an option to use an Xorshift RNG, enabled by compiling with -d:xorshift. This reduces runtime by approximately a quarter. Also fixed RNG generation to avoid u1 == 0 or u2 == 0.

import math, mersenne, unsigned

# It is sufficient to represent the rotation matrix
# [ [ cos x, -sin x], [sin x, cos x] ]]
# by just the values in the left column/bottom row.

type RotMat = tuple[sine, cosine: float64]

var tab1, tab2: array[65536, RotMat]

for i in 0..65535:
  tab1[i] = (sin(pi*float64(i-32768)/float64(65536)),
             cos(pi*float64(i-32768)/float64(65536)))
  tab2[i] = (sin(pi*float64(i)/float64(65536*65536)),
             cos(pi*float64(i)/float64(65536*65536)))

proc `*`(x, y: RotMat): RotMat {.inline.} =
  result.cosine = x.cosine*y.cosine-x.sine*y.sine
  result.sine = x.cosine*y.sine+x.sine*y.cosine

proc sincos(x: int): RotMat =
  let r1 = tab1[(x shr 16) and 65535]
  let r2 = tab2[x and 65535]
  return r1*r2

when defined(xorshift):
  type RNGState = object
    x, y, z, w: uint32

  proc getNum(s: var RNGState): int =
    let t = s.x xor (s.x shl 11)
    s.x = s.y; s.y = s.z; s.z = s.w
    s.w = s.w xor (s.w shr 19) xor t xor (t shr 8)
    result = int(s.w)

  var rng = RNGState(x: 1)
else:
  var rng = newMersenneTwister(1)

const maxRand = 1 shl 32

for i in 1..100_000_000:
  var t1, t2: int
  while true:
    t1 = rng.getNum; if t1 != 0: break
  let u1 = t1/maxRand
  while true:
    t2 = rng.getNum; if t2 != 0: break
  let u2 = t2/maxRand
  let w1 = pi*(u1-1/2)
  let w2 = -ln(u2)
  when defined(original):
    let (sw1, cw1) = (sin(w1), cos(w1))
  else:
    let (sw1, cw1) = sincos(t1)
  var x {.volatile.} = sw1/cw1*(pi/2-w1)+ln(w2*cw1/(pi/2-w1))
  when defined(validate):
    var x2 = sin(w1)/cos(w1)*(pi/2-w1)+ln(w2*cos(w1)/(pi/2-w1))
    echo sw1-sin(w1)
    echo cw1-cos(w1)
    echo x-x2

Answer 3

Numpy (score = 1.8)

In response to my previous code not being fast enough (I admit it was very slow), I used numpy which is extremely fast at working with large arrays of numbers and uses low level number libraries (BLAS). This is almost a direct translation of your code, but already is noticeably faster and has the bonus of being very readable.

import math
from numpy import *

pi = math.pi
half_pi = pi / 2

for i in range(10000):
    u1 = random.rand(10000)
    u2 = random.rand(10000)

    w1 = pi * (u1 - 1/2)
    w2 = -log(u2)

    a = half_pi - w1
    x = tan(w1)*a + log(w2*cos(w1)/a)

print("done")

I am running on i7 4770 with 8 gb RAM. Your code (including -O3 in GCC): 21.8 sec, my code: 11.8 sec.