| 554 | } |
| 555 | |
| 556 | double rk_gauss(rk_state *state) { |
| 557 | if (state->has_gauss) { |
| 558 | const double tmp = state->gauss; |
| 559 | state->gauss = 0; |
| 560 | state->has_gauss = 0; |
| 561 | return tmp; |
| 562 | } else { |
| 563 | double f, x1, x2, r2; |
| 564 | |
| 565 | do { |
| 566 | x1 = 2.0 * rk_double(state) - 1.0; |
| 567 | x2 = 2.0 * rk_double(state) - 1.0; |
| 568 | r2 = x1 * x1 + x2 * x2; |
| 569 | } while (r2 >= 1.0 || r2 == 0.0); |
| 570 | |
| 571 | /* Polar method, a more efficient version of the Box-Muller approach. */ |
| 572 | f = sqrt(-2.0 * log(r2) / r2); |
| 573 | /* Keep for next call */ |
| 574 | state->gauss = f * x1; |
| 575 | state->has_gauss = 1; |
| 576 | return f * x2; |
| 577 | } |
| 578 | } |