Mercurial > octave
changeset 25434:859ad1f0b85e
use templates to eliminate duplicate code in random number generators
* randgamma.h, randgamma.cc: Use templates to eliminate nearly duplicate
specializations.
* randpoisson.h, randpoisson.cc: Likewise.
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Sat, 17 Mar 2018 11:59:05 -0500 |
| parents | 49e0447413ad |
| children | a52e6fb674b1 |
| files | liboctave/numeric/randgamma.cc liboctave/numeric/randgamma.h liboctave/numeric/randpoisson.cc liboctave/numeric/randpoisson.h |
| diffstat | 4 files changed, 51 insertions(+), 244 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/numeric/randgamma.cc Sat Mar 17 07:02:30 2018 -0500 +++ b/liboctave/numeric/randgamma.cc Sat Mar 17 11:59:05 2018 -0500 @@ -89,37 +89,31 @@ namespace octave { - -#define INFINITE lo_ieee_isinf -#define RUNI rand_uniform<double> () -#define RNOR rand_normal<double> () -#define REXP rand_exponential<double> () - - template <> void rand_gamma<double> (double a, octave_idx_type n, double *r) + template <typename T> void rand_gamma (T a, octave_idx_type n, T *r) { octave_idx_type i; /* If a < 1, start by generating gamma (1+a) */ - const double d = (a < 1. ? 1.+a : a) - 1./3.; - const double c = 1./std::sqrt (9.*d); + const T d = (a < 1. ? 1.+a : a) - 1./3.; + const T c = 1./std::sqrt (9.*d); /* Handle invalid cases */ - if (a <= 0 || INFINITE(a)) + if (a <= 0 || lo_ieee_isinf (a)) { for (i=0; i < n; i++) - r[i] = numeric_limits<double>::NaN (); + r[i] = numeric_limits<T>::NaN (); return; } for (i=0; i < n; i++) { - double x, xsq, v, u; + T x, xsq, v, u; restart: - x = RNOR; + x = rand_normal<T> (); v = (1+c*x); v *= v*v; if (v <= 0) goto restart; /* rare, so don't bother moving up */ - u = RUNI; + u = rand_uniform<T> (); xsq = x*x; if (u >= 1.-0.0331*xsq*xsq && std::log (u) >= 0.5*xsq + d*(1-v+std::log (v))) goto restart; @@ -130,55 +124,10 @@ /* Use gamma(a) = gamma(1+a)*U^(1/a) */ /* Given REXP = -log(U) then U^(1/a) = exp(-REXP/a) */ for (i = 0; i < n; i++) - r[i] *= exp (-REXP/a); + r[i] *= exp (-rand_exponential<T> () / a); } } -#undef RUNI -#undef RNOR -#undef REXP - -#define RUNI rand_uniform<float> () -#define RNOR rand_normal<float> () -#define REXP rand_exponential<float> () - - template <> void rand_gamma<float> (float a, octave_idx_type n, float *r) - { - octave_idx_type i; - /* If a < 1, start by generating gamma(1+a) */ - const float d = (a < 1. ? 1.+a : a) - 1./3.; - const float c = 1./std::sqrt (9.*d); - - /* Handle invalid cases */ - if (a <= 0 || INFINITE(a)) - { - for (i=0; i < n; i++) - r[i] = numeric_limits<float>::NaN (); - return; - } - - for (i=0; i < n; i++) - { - float x, xsq, v, u; - frestart: - x = RNOR; - v = (1+c*x); - v *= v*v; - if (v <= 0) - goto frestart; /* rare, so don't bother moving up */ - u = RUNI; - xsq = x*x; - if (u >= 1.-0.0331*xsq*xsq && std::log (u) >= 0.5*xsq + d*(1-v+std::log (v))) - goto frestart; - r[i] = d*v; - } - if (a < 1) - { - /* Use gamma(a) = gamma(1+a)*U^(1/a) */ - /* Given REXP = -log(U) then U^(1/a) = exp(-REXP/a) */ - for (i = 0; i < n; i++) - r[i] *= exp (-REXP/a); - } - } - + template void rand_gamma (double, octave_idx_type, double *); + template void rand_gamma (float, octave_idx_type, float *); }
--- a/liboctave/numeric/randgamma.h Sat Mar 17 07:02:30 2018 -0500 +++ b/liboctave/numeric/randgamma.h Sat Mar 17 11:59:05 2018 -0500 @@ -34,12 +34,6 @@ void rand_gamma (T a, octave_idx_type n, T *p); - template <> void - rand_gamma<double> (double a, octave_idx_type n, double *p); - - template <> void - rand_gamma<float> (float a, octave_idx_type n, float *p); - template <typename T> T rand_gamma (T a)
--- a/liboctave/numeric/randpoisson.cc Sat Mar 17 07:02:30 2018 -0500 +++ b/liboctave/numeric/randpoisson.cc Sat Mar 17 11:59:05 2018 -0500 @@ -264,19 +264,23 @@ /* The remainder of the file is by Paul Kienzle */ + /* Table size is predicated on the maximum value of lambda + * we want to store in the table, and the maximum value of + * returned by the uniform random number generator on [0,1). + * With lambda==10 and u_max = 1 - 1/(2^32+1), we + * have poisson_pdf(lambda,36) < 1-u_max. If instead our + * generator uses more bits of mantissa or returns a value + * in the range [0,1], then for lambda==10 we need a table + * size of 46 instead. For long doubles, the table size + * will need to be longer still. */ +#define TABLESIZE 46 + /* Given uniform u, find x such that CDF(L,x)==u. Return x. */ - static void poisson_cdf_lookup (double lambda, double *p, size_t n) + + template <typename T> + static void + poisson_cdf_lookup (double lambda, T *p, size_t n) { - /* Table size is predicated on the maximum value of lambda - * we want to store in the table, and the maximum value of - * returned by the uniform random number generator on [0,1). - * With lambda==10 and u_max = 1 - 1/(2^32+1), we - * have poisson_pdf(lambda,36) < 1-u_max. If instead our - * generator uses more bits of mantissa or returns a value - * in the range [0,1], then for lambda==10 we need a table - * size of 46 instead. For long doubles, the table size - * will need to be longer still. */ -#define TABLESIZE 46 double t[TABLESIZE]; /* Precompute the table for the u up to and including 0.458. @@ -311,7 +315,7 @@ nextk: if (u <= t[k]) { - p[i] = static_cast<double> (k); + p[i] = static_cast<T> (k); continue; } if (++k < tableidx) @@ -333,81 +337,14 @@ /* We are assuming that the table size is big enough here. * This should be true even if rand_uniform is returning values in * the range [0,1] rather than [0,1). */ - p[i] = static_cast<double> (tableidx-1); - } - } - - static void poisson_cdf_lookup_float (double lambda, float *p, size_t n) - { - double t[TABLESIZE]; - - /* Precompute the table for the u up to and including 0.458. - * We will almost certainly need it. */ - int intlambda = static_cast<int> (std::floor (lambda)); - double P; - int tableidx; - size_t i = n; - - t[0] = P = exp (-lambda); - for (tableidx = 1; tableidx <= intlambda; tableidx++) - { - P = P*lambda/static_cast<double> (tableidx); - t[tableidx] = t[tableidx-1] + P; - } - - while (i-- > 0) - { - double u = octave::rand_uniform<double> (); - int k = (u > 0.458 ? intlambda : 0); - nextk: - if (u <= t[k]) - { - p[i] = static_cast<float> (k); - continue; - } - if (++k < tableidx) - goto nextk; - - while (tableidx < TABLESIZE) - { - P = P*lambda/static_cast<double> (tableidx); - t[tableidx] = t[tableidx-1] + P; - if (t[tableidx] == t[tableidx-1]) t[tableidx] = 1.0; - tableidx++; - if (u <= t[tableidx-1]) break; - } - - p[i] = static_cast<float> (tableidx-1); + p[i] = static_cast<T> (tableidx-1); } } /* From Press, et al., Numerical Recipes */ - static void poisson_rejection (double lambda, double *p, size_t n) - { - double sq = std::sqrt (2.0*lambda); - double alxm = std::log (lambda); - double g = lambda*alxm - xlgamma (lambda+1.0); - size_t i; - - for (i = 0; i < n; i++) - { - double y, em, t; - do - { - do - { - y = tan (M_PI*octave::rand_uniform<double> ()); - em = sq * y + lambda; - } while (em < 0.0); - em = std::floor (em); - t = 0.9*(1.0+y*y)*exp (em*alxm-flogfak (em)-g); - } while (octave::rand_uniform<double> () > t); - p[i] = em; - } - } - - /* From Press, et al., Numerical Recipes */ - static void poisson_rejection_float (double lambda, float *p, size_t n) + template <typename T> + static void + poisson_rejection (double lambda, T *p, size_t n) { double sq = std::sqrt (2.0*lambda); double alxm = std::log (lambda); @@ -442,17 +379,18 @@ * and large enough. */ /* Generate a set of poisson numbers with the same distribution */ - template <> void rand_poisson<double> (double L, octave_idx_type n, double *p) + template <typename T> void rand_poisson (T L_arg, octave_idx_type n, T *p) { + double L = L_arg; octave_idx_type i; if (L < 0.0 || lo_ieee_isinf (L)) { for (i=0; i<n; i++) - p[i] = numeric_limits<double>::NaN (); + p[i] = numeric_limits<T>::NaN (); } else if (L <= 10.0) { - poisson_cdf_lookup (L, p, n); + poisson_cdf_lookup<T> (L, p, n); } else if (L <= 1e8) { @@ -465,18 +403,22 @@ const double sqrtL = std::sqrt (L); for (i = 0; i < n; i++) { - p[i] = std::floor (rand_normal<double> () * sqrtL + L + 0.5); + p[i] = std::floor (rand_normal<T> () * sqrtL + L + 0.5); if (p[i] < 0.0) p[i] = 0.0; /* will probably never happen */ } } } + template void rand_poisson<double> (double, octave_idx_type, double *); + template void rand_poisson<float> (float, octave_idx_type, float *); + /* Generate one poisson variate */ - template <> double rand_poisson<double> (double L) + template <typename T> T rand_poisson (T L_arg) { - double ret; - if (L < 0.0) ret = numeric_limits<double>::NaN (); + double L = L_arg; + T ret; + if (L < 0.0) ret = numeric_limits<T>::NaN (); else if (L <= 12.0) { /* From Press, et al. Numerical recipes */ @@ -486,98 +428,30 @@ do { ++em; - t *= rand_uniform<double> (); + t *= rand_uniform<T> (); } while (t > g); ret = em; } else if (L <= 1e8) { /* numerical recipes */ - poisson_rejection (L, &ret, 1); + poisson_rejection<T> (L, &ret, 1); } else if (lo_ieee_isinf (L)) { /* FIXME: R uses NaN, but the normal approximation suggests that * limit should be Inf. Which is correct? */ - ret = numeric_limits<double>::NaN (); + ret = numeric_limits<T>::NaN (); } else { /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */ - ret = std::floor (rand_normal<double> () * std::sqrt (L) + L + 0.5); + ret = std::floor (rand_normal<T> () * std::sqrt (L) + L + 0.5); if (ret < 0.0) ret = 0.0; /* will probably never happen */ } return ret; } - /* Generate a set of poisson numbers with the same distribution */ - template <> void rand_poisson<float> (float FL, octave_idx_type n, float *p) - { - double L = FL; - octave_idx_type i; - if (L < 0.0 || lo_ieee_isinf (L)) - { - for (i=0; i<n; i++) - p[i] = numeric_limits<double>::NaN (); - } - else if (L <= 10.0) - { - poisson_cdf_lookup_float (L, p, n); - } - else if (L <= 1e8) - { - for (i=0; i<n; i++) - p[i] = pprsc (L); - } - else - { - /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */ - const double sqrtL = std::sqrt (L); - for (i = 0; i < n; i++) - { - p[i] = std::floor (rand_normal<float> () * sqrtL + L + 0.5); - if (p[i] < 0.0) - p[i] = 0.0; /* will probably never happen */ - } - } - } - - /* Generate one poisson variate */ - template <> float rand_poisson<float> (float FL) - { - double L = FL; - float ret; - if (L < 0.0) ret = numeric_limits<float>::NaN (); - else if (L <= 12.0) - { - /* From Press, et al. Numerical recipes */ - double g = exp (-L); - int em = -1; - double t = 1.0; - do - { - ++em; - t *= rand_uniform<float> (); - } while (t > g); - ret = em; - } - else if (L <= 1e8) - { - /* numerical recipes */ - poisson_rejection_float (L, &ret, 1); - } - else if (lo_ieee_isinf (L)) - { - /* FIXME: R uses NaN, but the normal approximation suggests that - * limit should be Inf. Which is correct? */ - ret = numeric_limits<float>::NaN (); - } - else - { - /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */ - ret = std::floor (rand_normal<float> () * std::sqrt (L) + L + 0.5); - if (ret < 0.0) ret = 0.0; /* will probably never happen */ - } - return ret; - } + template double rand_poisson<double> (double); + template float rand_poisson<float> (float); }
--- a/liboctave/numeric/randpoisson.h Sat Mar 17 07:02:30 2018 -0500 +++ b/liboctave/numeric/randpoisson.h Sat Mar 17 11:59:05 2018 -0500 @@ -30,19 +30,9 @@ namespace octave { - template <typename T> T rand_poisson (T L); - - template <> double rand_poisson<double> (double L); - template <> float rand_poisson<float> (float L); + template <typename T> void rand_poisson (T L, octave_idx_type n, T *p); - template <typename T> void - rand_poisson (T L, octave_idx_type n, T *p); - - template <> void - rand_poisson<double> (double L, octave_idx_type n, double *p); - - template <> void - rand_poisson<float> (float L, octave_idx_type n, float *p); + template <typename T> T rand_poisson (T L); } OCTAVE_DEPRECATED (4.4, "use 'octave::rand_poisson<double>' instead")