Function random_standard_normal_f

numpy/random/src/distributions/distributions.c:179–211

}

float random_standard_normal_f(bitgen_t *bitgen_state) {
  uint32_t r;
  int sign;
  uint32_t rabs;
  int idx;
  float x, xx, yy;
  for (;;) {
    /* r = n23sb8 */
    r = next_uint32(bitgen_state);
    idx = r & 0xff;
    sign = (r >> 8) & 0x1;
    rabs = (r >> 9) & 0x0007fffff;
    x = rabs * wi_float[idx];
    if (sign & 0x1)
      x = -x;
    if (rabs < ki_float[idx])
      return x; /* # 99.3% of the time return here */
    if (idx == 0) {
      for (;;) {
        /* Switch to 1.0 - U to avoid log(0.0), see GH 13361 */
        xx = -ziggurat_nor_inv_r_f * npy_log1pf(-next_float(bitgen_state));
        yy = -npy_log1pf(-next_float(bitgen_state));
        if (yy + yy > xx * xx)
          return ((rabs >> 8) & 0x1) ? -(ziggurat_nor_r_f + xx)
                                     : ziggurat_nor_r_f + xx;
      }
    } else {
      if (((fi_float[idx - 1] - fi_float[idx]) * next_float(bitgen_state) +
           fi_float[idx]) < exp(-0.5 * x * x))
        return x;
    }
  }
}

void random_standard_normal_fill_f(bitgen_t *bitgen_state, npy_intp cnt, float *out) {
  npy_intp i;