diff options
Diffstat (limited to 'gsl-1.9/randist/gauss.c')
| -rw-r--r-- | gsl-1.9/randist/gauss.c | 143 |
1 files changed, 143 insertions, 0 deletions
diff --git a/gsl-1.9/randist/gauss.c b/gsl-1.9/randist/gauss.c new file mode 100644 index 0000000..2ceb804 --- /dev/null +++ b/gsl-1.9/randist/gauss.c @@ -0,0 +1,143 @@ +/* randist/gauss.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2006 James Theiler, Brian Gough + * Copyright (C) 2006 Charles Karney + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#include <config.h> +#include <math.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_rng.h> +#include <gsl/gsl_randist.h> + +/* Of the two methods provided below, I think the Polar method is more + * efficient, but only when you are actually producing two random + * deviates. We don't produce two, because then we'd have to save one + * in a static variable for the next call, and that would screws up + * re-entrant or threaded code, so we only produce one. This makes + * the Ratio method suddenly more appealing. + * + * [Added by Charles Karney] We use Leva's implementation of the Ratio + * method which avoids calling log() nearly all the time and makes the + * Ratio method faster than the Polar method (when it produces just one + * result per call). Timing per call (gcc -O2 on 866MHz Pentium, + * average over 10^8 calls) + * + * Polar method: 660 ns + * Ratio method: 368 ns + * + */ + +/* Polar (Box-Mueller) method; See Knuth v2, 3rd ed, p122 */ + +double +gsl_ran_gaussian (const gsl_rng * r, const double sigma) +{ + double x, y, r2; + + do + { + /* choose x,y in uniform square (-1,-1) to (+1,+1) */ + x = -1 + 2 * gsl_rng_uniform_pos (r); + y = -1 + 2 * gsl_rng_uniform_pos (r); + + /* see if it is in the unit circle */ + r2 = x * x + y * y; + } + while (r2 > 1.0 || r2 == 0); + + /* Box-Muller transform */ + return sigma * y * sqrt (-2.0 * log (r2) / r2); +} + +/* Ratio method (Kinderman-Monahan); see Knuth v2, 3rd ed, p130. + * K+M, ACM Trans Math Software 3 (1977) 257-260. + * + * [Added by Charles Karney] This is an implementation of Leva's + * modifications to the original K+M method; see: + * J. L. Leva, ACM Trans Math Software 18 (1992) 449-453 and 454-455. */ + +double +gsl_ran_gaussian_ratio_method (const gsl_rng * r, const double sigma) +{ + double u, v, x, y, Q; + const double s = 0.449871; /* Constants from Leva */ + const double t = -0.386595; + const double a = 0.19600; + const double b = 0.25472; + const double r1 = 0.27597; + const double r2 = 0.27846; + + do /* This loop is executed 1.369 times on average */ + { + /* Generate a point P = (u, v) uniform in a rectangle enclosing + the K+M region v^2 <= - 4 u^2 log(u). */ + + /* u in (0, 1] to avoid singularity at u = 0 */ + u = 1 - gsl_rng_uniform (r); + + /* v is in the asymmetric interval [-0.5, 0.5). However v = -0.5 + is rejected in the last part of the while clause. The + resulting normal deviate is strictly symmetric about 0 + (provided that v is symmetric once v = -0.5 is excluded). */ + v = gsl_rng_uniform (r) - 0.5; + + /* Constant 1.7156 > sqrt(8/e) (for accuracy); but not by too + much (for efficiency). */ + v *= 1.7156; + + /* Compute Leva's quadratic form Q */ + x = u - s; + y = fabs (v) - t; + Q = x * x + y * (a * y - b * x); + + /* Accept P if Q < r1 (Leva) */ + /* Reject P if Q > r2 (Leva) */ + /* Accept if v^2 <= -4 u^2 log(u) (K+M) */ + /* This final test is executed 0.012 times on average. */ + } + while (Q >= r1 && (Q > r2 || v * v > -4 * u * u * log (u))); + + return sigma * (v / u); /* Return slope */ +} + +double +gsl_ran_gaussian_pdf (const double x, const double sigma) +{ + double u = x / fabs (sigma); + double p = (1 / (sqrt (2 * M_PI) * fabs (sigma))) * exp (-u * u / 2); + return p; +} + +double +gsl_ran_ugaussian (const gsl_rng * r) +{ + return gsl_ran_gaussian (r, 1.0); +} + +double +gsl_ran_ugaussian_ratio_method (const gsl_rng * r) +{ + return gsl_ran_gaussian_ratio_method (r, 1.0); +} + +double +gsl_ran_ugaussian_pdf (const double x) +{ + return gsl_ran_gaussian_pdf (x, 1.0); +} + |