# Super fast sampling from a tricky random distribution

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.

• If you want to avoid compilers optimizing away the entire computation, the easiest way is to declare x as volatile double x. – Reimer Behrends May 22 '14 at 14:15
• @ReimerBehrends I am happy with that solution! Thank you. Now I just need some answers to the question :) – user9206 May 22 '14 at 15:46
• Incidentally, is the line w2 = -log(u2) correct? After all, u2 can be zero. – Reimer Behrends May 23 '14 at 9:41
• Looking at it more closely, I suspect that both u1 and u2 should be drawn from (0, 1) rather than [0,1). Also, 1/2 should probably be 1/2.0 in C (since / is integer division if both operands are ints). – Reimer Behrends May 23 '14 at 10:34
• @ReimerBehrends I will check the 1/2 possible bug. Thank you. For the other part the formula is taken from Table 1 of arxiv.org/pdf/0908.3961v2.pdf so if my implementation doesn't follow that then that is my mistake. I suppose mathematically the probability that u1 or u2 are exactly 0 is exactly 0. It's only with limited precision that you get a non-zero probability. – user9206 May 23 '14 at 13:20

# 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))

• I would like to remove the printing. The problem is that the whole code might be optimized out then! – user9206 May 17 '14 at 21:47
• @Lembik Write it to file? – qwr May 17 '14 at 22:47
• I added an if statement which should reduce the output size by a factor of about 10000. – user9206 May 18 '14 at 7:15
• @Lembik: if statements are slow. – Kyle Kanos May 18 '14 at 19:05
• @Lembik: The code needs to test the inequality each iteration, which is what makes it slow. A faster way to ensure it isn't optimized out is to just sum x and print it out at the end (or something similar). – Kyle Kanos May 18 '14 at 19:16

# 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

• Thank you! Do you know any faster RNGs which pass 15 diehard tests? – user9206 May 23 '14 at 6:50
• Certain variants of Xorshift should fit the bill, though I can't say that I'm much of an expert on PRNGs. – Reimer Behrends May 23 '14 at 9:36
• What is your current score using xorshift? Also, what do you think the prospects are for further speed ups? It seems sse might help...maybe. – user9206 May 25 '14 at 19:16
• The current score (on my machine, this can vary greatly by hardware) is about 2 (compared to compiling with -d:original). That's assuming that Xorshift is an acceptable PRNG. The limiting factor seems to be the need for at least 32 bits of precision. This makes it difficult to use cheaper approximations or to create a polynomial approximation of the formula as a whole. – Reimer Behrends May 26 '14 at 10:00
• The problem is that the distribution is dramatically skewed. Try plotting it and you will see. I am not sure how to capture this long tail without using a lot of bits. – user9206 May 26 '14 at 16:48

# 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.

• Thanks for this but your code takes 23 seconds compared to 20 seconds for mine on my machine. What timings are you getting and are you compiling mine with -03? – user9206 Jun 1 '14 at 6:56
• @Lembik I am running on i7 4770. Your code (including -O3): 21.8 sec, my code: 11.8 sec. I find it unlikely that gcc could beat numpy, because numpy operates efficiently on arrays and uses BLAS. Maybe try it again? I'm not sure what could have caused this difference. – qwr Jun 1 '14 at 21:46