5 #include <gsl/gsl_sf_bessel.h>
6 #include <gsl/gsl_rng.h>
7 #include <gsl/gsl_randist.h>
9 static gsl_rng *rng; /* global generator */
11 void var_seed(uint64_t seed) {
12 gsl_rng_set(rng, seed);
15 void var_init(uint64_t seed) {
16 const gsl_rng_type * T;
20 rng = gsl_rng_alloc (T);
23 seed = (unsigned long)getpid();
29 unsigned int var_rhyper(unsigned int t, unsigned int n1, unsigned int n2) {
30 return gsl_ran_hypergeometric(rng, n1, n2, t);
33 /* must be two ints in each vector */
34 unsigned int var_rhyperv(unsigned int t, unsigned int *s, unsigned int *d) {
36 x = var_rhyper(t, s[0], s[1]);
37 d[0] = x; d[1] = t - x;
41 void var_rbinom(unsigned int t, unsigned int *n, unsigned int *d) {
42 d[0] = gsl_ran_binomial(rng, (double)n[0]/((double)n[0]+(double)n[1]), t);
46 /* TODO hardcoded to size of arrays == 2 */
47 void var_rmultinom(unsigned int t, unsigned int *n, unsigned int *d) {
51 p[0] = n[0]; p[1] = n[1];
52 gsl_ran_multinomial(rng, k, t, p, d);
55 /* pick one given n weighted probabilities */
56 int var_pick(unsigned int n, unsigned int *p) {
65 r = gsl_rng_uniform_int(rng, sum);
70 if (sum >= r) return i;