dynamic system, 1/x stationary distribution

#include <math.h>
#include <stdio.h>
#include <assert.h>

#define MAX_ITER 100000000
#define K 1024
#define B 7

typedef unsigned long long int uint64;
typedef unsigned long int uint32;
int main() {
  const int SHIFT = ((int)(log(B)/log(2.)));
  const int THRESH = (2 * K) << SHIFT;
  const double norm = 100. / MAX_ITER;

  uint64 histo[K] = { 0 };
  uint32 x =  2 * K - 1;

  int n;
  for (n = 0; n < MAX_ITER; ++n) {
    x = x * B + (rand() % B);
    const int shift = SHIFT + (x >= THRESH);
    x >>= shift;
    assert(x >= K && x < 2 * K);
    ++histo[x - K];
  }

  // print histo
  for (n = 0; n < K; ++n) {
    printf("%d %.5lf\n", n + K, norm * histo[n]);
  }
  return 0;
}