zipf.cpp

#include <cassert>
#include <cstdio>
#include <cmath>
#include <random>

int zipf(double alpha, int n, std::mt19937_64 &rng)
{
  static bool first = true;      // Static first time flag
  static double c = 0;          // Normalization constant
  static double *sum_probs;     // Pre-calculated sum of probabilities
  double z;                     // Uniform random number (0 < z < 1)
  int zipf_value;               // Computed exponential value to be returned
  int    i;                     // Loop counter
  int low, high, mid;           // Binary-search bounds

  // Compute normalization constant on first call only
  if (first == true)
  {
    for (i=1; i<=n; i++)
      c = c + (1.0 / pow((double) i, alpha));
    c = 1.0 / c;

    sum_probs = new double[(n+1)*sizeof(*sum_probs)];
    sum_probs[0] = 0;
    for (i=1; i<=n; i++) {
      sum_probs[i] = sum_probs[i-1] + c / pow((double) i, alpha);
    }
    first = false;
  }

  // Pull a uniform random number (0 < z < 1)
  do
  {
    std::uniform_real_distribution<double> unif(0, 1);
    z = unif(rng);
  }
  while ((z == 0) || (z == 1));

  // Map z to the value
  low = 1, high = n;
  do {
    mid = floor((low+high)/2);
    if (sum_probs[mid] >= z && sum_probs[mid-1] < z) {
      zipf_value = mid;
      break;
    } else if (sum_probs[mid] >= z) {
      high = mid-1;
    } else {
      low = mid+1;
    }
  } while (low <= high);

  // Assert that zipf_value is between 1 and N
  assert((zipf_value >=1) && (zipf_value <= n));

  return(zipf_value);
}

int main() {
    std::mt19937_64 rng;
    std::seed_seq seq{1,2,3,4,5};
    rng.seed(seq);

    int N = 10;
    double alpha = 0.75;
    for (int i = 0; i < 40; i++) {
        printf("%d ", zipf(alpha, N, rng));
    }

    return 0;
}