#include <stdint.h>
#include <stdio.h>
#include <unistd.h>

#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

static gsl_rng *rng;  /* global generator */

void var_seed(uint64_t seed) {
	gsl_rng_set(rng, seed);
}

void var_init(uint64_t seed) {
	const gsl_rng_type * T;

	gsl_rng_env_setup();
	T = gsl_rng_default;
	rng = gsl_rng_alloc (T);

	if (seed == 0) {
		seed = (unsigned long)getpid();
	}

	var_seed(seed);
}

unsigned int var_rhyper(unsigned int t, unsigned int n1, unsigned int n2) {
	return gsl_ran_hypergeometric(rng, n1, n2, t);
}

/* must be two ints in each vector */
unsigned int var_rhyperv(unsigned int t, unsigned int *s, unsigned int *d) {
	int x;
	x = var_rhyper(t, s[0], s[1]);
	d[0] = x; d[1] = t - x;
	return x;
}

void var_rbinom(unsigned int t, unsigned int *n, unsigned int *d) {
	d[0] = gsl_ran_binomial(rng, (double)n[0]/((double)n[0]+(double)n[1]), t);
	d[1] = t - d[0];
}

/* TODO hardcoded to size of arrays == 2 */
void var_rmultinom(unsigned int t, unsigned int *n, unsigned int *d) {
	double p[2];
	size_t k = 2;

	p[0] = n[0]; p[1] = n[1];
	gsl_ran_multinomial(rng, k, t, p, d);
}

/* pick one given n weighted probabilities */
int var_pick(unsigned int n, unsigned int *p) {
	long sum;
	int i;
	unsigned long r;

	sum = 0;
	for (i=0;i<n;i++) {
		sum += p[i];
	}
	r = gsl_rng_uniform_int(rng, sum);

	sum = 0;
	for (i=0;i<n;i++) {
		sum += p[i];
		if (sum >= r) return i;
	}
	return i;
}