libc random number generator flawed?

Question

Consider an algorithm to test the probability that a certain number is picked from a set of N unique numbers after a specific number of tries (for example, with N=2, what's the probability in Roulette (without 0) that it takes X tries for Black to win?).

The correct distribution for this is pow(1-1/N,X-1)*(1/N).

However, when I test this using the following code, there is always a deep ditch at X=31, independently from N, and independently from the seed.

Is this an intrinsic flaw that cannot be prevented due to the implementation specifics of the PRNG in use, is this a real bug, or am I overlooking something obvious?

// C

#include <sys/times.h>
#include <math.h>
#include <stdio.h>

int array[101];
void main(){

    int nsamples=10000000;
    double breakVal,diffVal;
    int i,cnt;

    // seed, but doesn't change anything
    struct tms time;
    srandom(times(&time));

    // sample
    for(i=0;i<nsamples;i++){
        cnt=1;
        do{
            if((random()%36)==0) // break if 0 is chosen
                break;
            cnt++;
        }while(cnt<100);
        array[cnt]++;
    }

    // show distribution
    for(i=1;i<100;i++){
        breakVal=array[i]/(double)nsamples; // normalize
        diffVal=breakVal-pow(1-1/36.,i-1)*1/36.; // difference to expected value
        printf("%d %.12g %.12g\n",i,breakVal,diffVal);
    }
}

Tested on an up-to-date Xubuntu 12.10 with libc6 package 2.15-0ubuntu20 and Intel Core i5-2500 SandyBridge, but I discovered this already a few years ago on an older Ubuntu machine.

I also tested this on Windows 7 using Unity3D/Mono (not sure which Mono version, though), and here the ditch happens at X=55 when using System.Random, while Unity's builtin Unity.Random has no visible ditch (at least not for X<100).

The distribution:

The differences:

Answer 1

This is due to glibc's random() function not being random enough. According to this page, for the random numbers returned by random(), we have:

oi = (oi-3 + oi-31) % 2^31

or:

oi = (oi-3 + oi-31 + 1) % 2^31.

Now take xi = oi % 36, and suppose the first equation above is the one used (this happens with a 50% chance for each number). Now if xi-31=0 and xi-3!=0, then the chance that xi=0 is less than 1/36. This is because 50% of the time oi-31 + oi-3 will be less than 2^31, and when that happens,

xi = oi % 36 = (oi-3 + oi-31) % 36 = oi-3 % 36 = xi-3,

which is nonzero. This causes the ditch you see 31 samples after a 0 sample.