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 &lt;stdio.h&gt;
#include &lt;stdlib.h&gt;
#include "mtwist.h"
#include &lt;math.h&gt;


int main(void) {
   int i;
   volatile double x;
   mt_seed();
   double u1;
   double u2;
   double w1;
   double w2;
   for(i = 0; i &lt; 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.
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)