1
#include "include/legacy-distributions.h"
2

3

4
static NPY_INLINE double legacy_double(aug_bitgen_t *aug_state) {
5 1
  return aug_state->bit_generator->next_double(aug_state->bit_generator->state);
6
}
7

8 1
double legacy_gauss(aug_bitgen_t *aug_state) {
9 1
  if (aug_state->has_gauss) {
10 1
    const double temp = aug_state->gauss;
11 1
    aug_state->has_gauss = false;
12 1
    aug_state->gauss = 0.0;
13 1
    return temp;
14
  } else {
15
    double f, x1, x2, r2;
16

17
    do {
18 1
      x1 = 2.0 * legacy_double(aug_state) - 1.0;
19 1
      x2 = 2.0 * legacy_double(aug_state) - 1.0;
20 1
      r2 = x1 * x1 + x2 * x2;
21 1
    } while (r2 >= 1.0 || r2 == 0.0);
22

23
    /* Polar method, a more efficient version of the Box-Muller approach. */
24 1
    f = sqrt(-2.0 * log(r2) / r2);
25
    /* Keep for next call */
26 1
    aug_state->gauss = f * x1;
27 1
    aug_state->has_gauss = true;
28 1
    return f * x2;
29
  }
30
}
31

32 1
double legacy_standard_exponential(aug_bitgen_t *aug_state) {
33
  /* We use -log(1-U) since U is [0, 1) */
34 1
  return -log(1.0 - legacy_double(aug_state));
35
}
36

37 1
double legacy_standard_gamma(aug_bitgen_t *aug_state, double shape) {
38
  double b, c;
39
  double U, V, X, Y;
40

41 1
  if (shape == 1.0) {
42 1
    return legacy_standard_exponential(aug_state);
43
  }
44 1
  else if (shape == 0.0) {
45
    return 0.0;
46 1
  } else if (shape < 1.0) {
47
    for (;;) {
48 1
      U = legacy_double(aug_state);
49 1
      V = legacy_standard_exponential(aug_state);
50 1
      if (U <= 1.0 - shape) {
51 1
        X = pow(U, 1. / shape);
52 1
        if (X <= V) {
53
          return X;
54
        }
55
      } else {
56 1
        Y = -log((1 - U) / shape);
57 1
        X = pow(1.0 - shape + shape * Y, 1. / shape);
58 1
        if (X <= (V + Y)) {
59
          return X;
60
        }
61
      }
62
    }
63
  } else {
64 1
    b = shape - 1. / 3.;
65 1
    c = 1. / sqrt(9 * b);
66
    for (;;) {
67
      do {
68 1
        X = legacy_gauss(aug_state);
69 1
        V = 1.0 + c * X;
70 1
      } while (V <= 0.0);
71

72 1
      V = V * V * V;
73 1
      U = legacy_double(aug_state);
74 1
      if (U < 1.0 - 0.0331 * (X * X) * (X * X))
75 1
        return (b * V);
76 1
      if (log(U) < 0.5 * X * X + b * (1. - V + log(V)))
77 1
        return (b * V);
78
    }
79
  }
80
}
81

82 1
double legacy_gamma(aug_bitgen_t *aug_state, double shape, double scale) {
83 1
  return scale * legacy_standard_gamma(aug_state, shape);
84
}
85

86 1
double legacy_pareto(aug_bitgen_t *aug_state, double a) {
87 1
  return exp(legacy_standard_exponential(aug_state) / a) - 1;
88
}
89

90 1
double legacy_weibull(aug_bitgen_t *aug_state, double a) {
91 1
  if (a == 0.0) {
92
    return 0.0;
93
  }
94 1
  return pow(legacy_standard_exponential(aug_state), 1. / a);
95
}
96

97 1
double legacy_power(aug_bitgen_t *aug_state, double a) {
98 1
  return pow(1 - exp(-legacy_standard_exponential(aug_state)), 1. / a);
99
}
100

101 1
double legacy_chisquare(aug_bitgen_t *aug_state, double df) {
102 1
  return 2.0 * legacy_standard_gamma(aug_state, df / 2.0);
103
}
104

105 1
double legacy_noncentral_chisquare(aug_bitgen_t *aug_state, double df,
106
                                   double nonc) {
107
  double out;
108 1
  if (nonc == 0) {
109 1
    return legacy_chisquare(aug_state, df);
110
  }
111 1
  if (1 < df) {
112 1
    const double Chi2 = legacy_chisquare(aug_state, df - 1);
113 1
    const double n = legacy_gauss(aug_state) + sqrt(nonc);
114 1
    return Chi2 + n * n;
115
  } else {
116 1
    const long i = random_poisson(aug_state->bit_generator, nonc / 2.0);
117 1
    out = legacy_chisquare(aug_state, df + 2 * i);
118
    /* Insert nan guard here to avoid changing the stream */
119 1
    if (npy_isnan(nonc)){
120
      return NPY_NAN;
121
    } else {
122 1
    return out;
123
    }
124
  }
125
}
126

127 1
double legacy_noncentral_f(aug_bitgen_t *aug_state, double dfnum, double dfden,
128
                           double nonc) {
129 1
  double t = legacy_noncentral_chisquare(aug_state, dfnum, nonc) * dfden;
130 1
  return t / (legacy_chisquare(aug_state, dfden) * dfnum);
131
}
132

133 1
double legacy_wald(aug_bitgen_t *aug_state, double mean, double scale) {
134
  double U, X, Y;
135
  double mu_2l;
136

137 1
  mu_2l = mean / (2 * scale);
138 1
  Y = legacy_gauss(aug_state);
139 1
  Y = mean * Y * Y;
140 1
  X = mean + mu_2l * (Y - sqrt(4 * scale * Y + Y * Y));
141 1
  U = legacy_double(aug_state);
142 1
  if (U <= mean / (mean + X)) {
143
    return X;
144
  } else {
145 1
    return mean * mean / X;
146
  }
147
}
148

149 1
double legacy_normal(aug_bitgen_t *aug_state, double loc, double scale) {
150 1
  return loc + scale * legacy_gauss(aug_state);
151
}
152

153 1
double legacy_lognormal(aug_bitgen_t *aug_state, double mean, double sigma) {
154 1
  return exp(legacy_normal(aug_state, mean, sigma));
155
}
156

157 1
double legacy_standard_t(aug_bitgen_t *aug_state, double df) {
158
  double num, denom;
159

160 1
  num = legacy_gauss(aug_state);
161 1
  denom = legacy_standard_gamma(aug_state, df / 2);
162 1
  return sqrt(df / 2) * num / sqrt(denom);
163
}
164

165 1
int64_t legacy_negative_binomial(aug_bitgen_t *aug_state, double n, double p) {
166 1
  double Y = legacy_gamma(aug_state, n, (1 - p) / p);
167 1
  return (int64_t)random_poisson(aug_state->bit_generator, Y);
168
}
169

170 1
double legacy_standard_cauchy(aug_bitgen_t *aug_state) {
171 1
  return legacy_gauss(aug_state) / legacy_gauss(aug_state);
172
}
173

174 1
double legacy_beta(aug_bitgen_t *aug_state, double a, double b) {
175
  double Ga, Gb;
176

177 1
  if ((a <= 1.0) && (b <= 1.0)) {
178
    double U, V, X, Y;
179
    /* Use Johnk's algorithm */
180

181
    while (1) {
182 1
      U = legacy_double(aug_state);
183 1
      V = legacy_double(aug_state);
184 1
      X = pow(U, 1.0 / a);
185 1
      Y = pow(V, 1.0 / b);
186

187 1
      if ((X + Y) <= 1.0) {
188 1
        if (X + Y > 0) {
189 1
          return X / (X + Y);
190
        } else {
191 1
          double logX = log(U) / a;
192 1
          double logY = log(V) / b;
193 1
          double logM = logX > logY ? logX : logY;
194 1
          logX -= logM;
195 1
          logY -= logM;
196

197 1
          return exp(logX - log(exp(logX) + exp(logY)));
198
        }
199
      }
200
    }
201
  } else {
202 1
    Ga = legacy_standard_gamma(aug_state, a);
203 1
    Gb = legacy_standard_gamma(aug_state, b);
204 1
    return Ga / (Ga + Gb);
205
  }
206
}
207

208 1
double legacy_f(aug_bitgen_t *aug_state, double dfnum, double dfden) {
209 1
  return ((legacy_chisquare(aug_state, dfnum) * dfden) /
210 1
          (legacy_chisquare(aug_state, dfden) * dfnum));
211
}
212

213 1
double legacy_exponential(aug_bitgen_t *aug_state, double scale) {
214 1
  return scale * legacy_standard_exponential(aug_state);
215
}
216

217

218 1
static RAND_INT_TYPE legacy_random_binomial_original(bitgen_t *bitgen_state,
219
                                                     double p,
220
                                                     RAND_INT_TYPE n,
221
                                                     binomial_t *binomial) {
222
  double q;
223

224 1
  if (p <= 0.5) {
225 1
    if (p * n <= 30.0) {
226 1
      return random_binomial_inversion(bitgen_state, n, p, binomial);
227
    } else {
228 1
      return random_binomial_btpe(bitgen_state, n, p, binomial);
229
    }
230
  } else {
231 1
    q = 1.0 - p;
232 1
    if (q * n <= 30.0) {
233 1
      return n - random_binomial_inversion(bitgen_state, n, q, binomial);
234
    } else {
235 1
      return n - random_binomial_btpe(bitgen_state, n, q, binomial);
236
    }
237
  }
238
}
239

240

241 1
int64_t legacy_random_binomial(bitgen_t *bitgen_state, double p,
242
                               int64_t n, binomial_t *binomial) {
243 1
  return (int64_t) legacy_random_binomial_original(bitgen_state, p,
244
                                                   (RAND_INT_TYPE) n,
245
                                                   binomial);
246
}
247

248

249 1
static RAND_INT_TYPE random_hypergeometric_hyp(bitgen_t *bitgen_state,
250
                                               RAND_INT_TYPE good,
251
                                               RAND_INT_TYPE bad,
252
                                               RAND_INT_TYPE sample) {
253
  RAND_INT_TYPE d1, k, z;
254
  double d2, u, y;
255

256 1
  d1 = bad + good - sample;
257 1
  d2 = (double)MIN(bad, good);
258

259 1
  y = d2;
260 1
  k = sample;
261 1
  while (y > 0.0) {
262 1
    u = next_double(bitgen_state);
263 1
    y -= (RAND_INT_TYPE)floor(u + y / (d1 + k));
264 1
    k--;
265 1
    if (k == 0)
266
      break;
267
  }
268 1
  z = (RAND_INT_TYPE)(d2 - y);
269 1
  if (good > bad)
270 1
    z = sample - z;
271 1
  return z;
272
}
273

274
/* D1 = 2*sqrt(2/e) */
275
/* D2 = 3 - 2*sqrt(3/e) */
276
#define D1 1.7155277699214135
277
#define D2 0.8989161620588988
278 1
static RAND_INT_TYPE random_hypergeometric_hrua(bitgen_t *bitgen_state,
279
                                                RAND_INT_TYPE good,
280
                                                RAND_INT_TYPE bad,
281
                                                RAND_INT_TYPE sample) {
282
  RAND_INT_TYPE mingoodbad, maxgoodbad, popsize, m, d9;
283
  double d4, d5, d6, d7, d8, d10, d11;
284
  RAND_INT_TYPE Z;
285
  double T, W, X, Y;
286

287 1
  mingoodbad = MIN(good, bad);
288 1
  popsize = good + bad;
289 1
  maxgoodbad = MAX(good, bad);
290 1
  m = MIN(sample, popsize - sample);
291 1
  d4 = ((double)mingoodbad) / popsize;
292 1
  d5 = 1.0 - d4;
293 1
  d6 = m * d4 + 0.5;
294 1
  d7 = sqrt((double)(popsize - m) * sample * d4 * d5 / (popsize - 1) + 0.5);
295 1
  d8 = D1 * d7 + D2;
296 1
  d9 = (RAND_INT_TYPE)floor((double)(m + 1) * (mingoodbad + 1) / (popsize + 2));
297 1
  d10 = (random_loggam(d9 + 1) + random_loggam(mingoodbad - d9 + 1) +
298 1
         random_loggam(m - d9 + 1) + random_loggam(maxgoodbad - m + d9 + 1));
299 1
  d11 = MIN(MIN(m, mingoodbad) + 1.0, floor(d6 + 16 * d7));
300
  /* 16 for 16-decimal-digit precision in D1 and D2 */
301

302
  while (1) {
303 1
    X = next_double(bitgen_state);
304 1
    Y = next_double(bitgen_state);
305 1
    W = d6 + d8 * (Y - 0.5) / X;
306

307
    /* fast rejection: */
308 1
    if ((W < 0.0) || (W >= d11))
309 1
      continue;
310

311 1
    Z = (RAND_INT_TYPE)floor(W);
312 1
    T = d10 - (random_loggam(Z + 1) + random_loggam(mingoodbad - Z + 1) +
313 1
               random_loggam(m - Z + 1) + random_loggam(maxgoodbad - m + Z + 1));
314

315
    /* fast acceptance: */
316 1
    if ((X * (4.0 - X) - 3.0) <= T)
317
      break;
318

319
    /* fast rejection: */
320 1
    if (X * (X - T) >= 1)
321 1
      continue;
322
    /* log(0.0) is ok here, since always accept */
323 1
    if (2.0 * log(X) <= T)
324
      break; /* acceptance */
325
  }
326

327
  /* this is a correction to HRUA* by Ivan Frohne in rv.py */
328 1
  if (good > bad)
329 1
    Z = m - Z;
330

331
  /* another fix from rv.py to allow sample to exceed popsize/2 */
332 1
  if (m < sample)
333 1
    Z = good - Z;
334

335 1
  return Z;
336
}
337
#undef D1
338
#undef D2
339

340 1
static RAND_INT_TYPE random_hypergeometric_original(bitgen_t *bitgen_state,
341
                                                    RAND_INT_TYPE good,
342
                                                    RAND_INT_TYPE bad,
343
                                                    RAND_INT_TYPE sample)
344
{
345 1
  if (sample > 10) {
346 1
    return random_hypergeometric_hrua(bitgen_state, good, bad, sample);
347 1
  } else if (sample > 0) {
348 1
    return random_hypergeometric_hyp(bitgen_state, good, bad, sample);
349
  } else {
350
    return 0;
351
  }
352
}
353

354

355
/*
356
 * This is a wrapper function that matches the expected template. In the legacy
357
 * generator, all int types are long, so this accepts int64 and then converts
358
 * them to longs. These values must be in bounds for long and this is checked
359
 * outside this function
360
 *
361
 * The remaining are included for the return type only
362
 */
363 1
int64_t legacy_random_hypergeometric(bitgen_t *bitgen_state, int64_t good,
364
                                     int64_t bad, int64_t sample) {
365 1
  return (int64_t)random_hypergeometric_original(bitgen_state,
366
                                                 (RAND_INT_TYPE)good,
367
                                                 (RAND_INT_TYPE)bad,
368
                                                 (RAND_INT_TYPE)sample);
369
}
370

371

372 1
int64_t legacy_random_logseries(bitgen_t *bitgen_state, double p) {
373 1
  return (int64_t)random_logseries(bitgen_state, p);
374
}
375

376 1
 int64_t legacy_random_poisson(bitgen_t *bitgen_state, double lam) {
377 1
  return (int64_t)random_poisson(bitgen_state, lam);
378
}
379

380 1
 int64_t legacy_random_zipf(bitgen_t *bitgen_state, double a) {
381 1
  return (int64_t)random_zipf(bitgen_state, a);
382
}
383

384 1
 int64_t legacy_random_geometric(bitgen_t *bitgen_state, double p) {
385 1
  return (int64_t)random_geometric(bitgen_state, p);
386
}
387

388 1
 void legacy_random_multinomial(bitgen_t *bitgen_state, RAND_INT_TYPE n,
389
                               RAND_INT_TYPE *mnix, double *pix, npy_intp d,
390
                               binomial_t *binomial) {
391 1
  return random_multinomial(bitgen_state, n, mnix, pix, d, binomial);
392
}

Read our documentation on viewing source code .

Loading