Function random_wald

numpy/random/src/distributions/distributions.c:839–853

}

double random_wald(bitgen_t *bitgen_state, double mean, double scale) {
  double U, X, Y;
  double mu_2l;

  mu_2l = mean / (2 * scale);
  Y = random_standard_normal(bitgen_state);
  Y = mean * Y * Y;
  X = mean + mu_2l * (Y - sqrt(4 * scale * Y + Y * Y));
  U = next_double(bitgen_state);
  if (U <= mean / (mean + X)) {
    return X;
  } else {
    return mean * mean / X;
  }
}

double random_vonmises(bitgen_t *bitgen_state, double mu, double kappa) {
  double s;