Mercurial > hg > octave-lyh
changeset 14655:43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
* oct-rand.cc (octave_rand:do_matrix): Remove method.
(octave_rand::do_float_scalar, octave_rand::do_float_nd_array,
octave_rand::do_float_vector): New methods.
* oct-rand.h (octave_rand:do_matrix, octave_rand::matrix): Remove methods.
(octave_rand::do_float_scalar, octave_rand::do_float_nd_array,
octave_rand::do_float_vector, octave_rand::float_scalar,
octave_rand::do_nd_array, octave_rand::float_vector): Declare new methods.
* randgamma.c (oct_fill_float_randg, oct_float_randg): New functions.
* randgamma.h (oct_fill_float_randg, oct_float_randg): Declare new functions.
* randpoisson.c (oct_fill_float_randp, oct_float_randp): New functions.
(poisson_cdf_lookup_float, poisson_rejection_float): New static functions.
* randpoisson.h (oct_fill_float_randp, oct_float_randp): Declare new functions.
* randmtzig.c (randu64) : Remove static function.
(ALLBITS): Remove compile flag logic.
(randu32): Change return type to float.
(oct_float_randu, oct_float_randn, oct_float_rande, oct_fill_float_randu,
oct_fill_float_randn, oct_fill_float_rande): New functions.
(create_ziggurat_float_tables): New static function.
* randmtzig.h (oct_float_randu, oct_float_randn, oct_float_rande,
oct_fill_float_randu, oct_fill_float_randn, oct_fill_float_rande):
Declare new functions.
* rand.cc (do_rand): Add logic to parse "double" or "single" as final argument.
author | David Bateman <dbateman@free.fr> |
---|---|
date | Sat, 19 May 2012 00:14:59 +0200 |
parents | 389f49a88656 |
children | 8d2ce821e38a |
files | liboctave/oct-rand.cc liboctave/oct-rand.h liboctave/randgamma.c liboctave/randgamma.h liboctave/randmtzig.c liboctave/randmtzig.h liboctave/randpoisson.c liboctave/randpoisson.h src/DLD-FUNCTIONS/rand.cc |
diffstat | 9 files changed, 747 insertions(+), 104 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/oct-rand.cc +++ b/liboctave/oct-rand.cc @@ -419,19 +419,120 @@ return retval; } -Matrix -octave_rand::do_matrix (octave_idx_type n, octave_idx_type m, double a) +float +octave_rand::do_float_scalar (float a) { - Matrix retval; + float retval = 0.0; + + if (use_old_generators) + { + double da = a; + double dretval; + switch (current_distribution) + { + case uniform_dist: + F77_FUNC (dgenunf, DGENUNF) (0.0, 1.0, dretval); + break; + + case normal_dist: + F77_FUNC (dgennor, DGENNOR) (0.0, 1.0, dretval); + break; + + case expon_dist: + F77_FUNC (dgenexp, DGENEXP) (1.0, dretval); + break; - if (n >= 0 && m >= 0) - { - retval.clear (n, m); + case poisson_dist: + if (da < 0.0 || xisnan(da) || xisinf(da)) + dretval = octave_NaN; + else + { + // workaround bug in ignpoi, by calling with different Mu + F77_FUNC (dignpoi, DIGNPOI) (da + 1, dretval); + F77_FUNC (dignpoi, DIGNPOI) (da, dretval); + } + break; - if (n > 0 && m > 0) - fill (retval.capacity(), retval.fortran_vec(), a); + case gamma_dist: + if (da <= 0.0 || xisnan(da) || xisinf(da)) + retval = octave_NaN; + else + F77_FUNC (dgengam, DGENGAM) (1.0, da, dretval); + break; + + default: + (*current_liboctave_error_handler) + ("rand: invalid distribution ID = %d", current_distribution); + break; + } + retval = dretval; } else + { + switch (current_distribution) + { + case uniform_dist: + retval = oct_float_randu (); + break; + + case normal_dist: + retval = oct_float_randn (); + break; + + case expon_dist: + retval = oct_float_rande (); + break; + + case poisson_dist: + // Keep poisson distribution in double precision for accuracy + retval = oct_randp (a); + break; + + case gamma_dist: + retval = oct_float_randg (a); + break; + + default: + (*current_liboctave_error_handler) + ("rand: invalid distribution ID = %d", current_distribution); + break; + } + + save_state (); + } + + return retval; +} + +Array<double> +octave_rand::do_vector (octave_idx_type n, double a) +{ + Array<double> retval; + + if (n > 0) + { + retval.clear (n, 1); + + fill (retval.capacity(), retval.fortran_vec(), a); + } + else if (n < 0) + (*current_liboctave_error_handler) ("rand: invalid negative argument"); + + return retval; +} + +Array<float> +octave_rand::do_float_vector (octave_idx_type n, float a) +{ + Array<float> retval; + + if (n > 0) + { + retval.clear (n, 1); + + fill (retval.capacity(), retval.fortran_vec(), a); + } + else if (n < 0) (*current_liboctave_error_handler) ("rand: invalid negative argument"); return retval; @@ -452,19 +553,17 @@ return retval; } -Array<double> -octave_rand::do_vector (octave_idx_type n, double a) +FloatNDArray +octave_rand::do_float_nd_array (const dim_vector& dims, float a) { - Array<double> retval; + FloatNDArray retval; - if (n > 0) + if (! dims.all_zero ()) { - retval.clear (n, 1); + retval.clear (dims); - fill (retval.capacity (), retval.fortran_vec (), a); + fill (retval.capacity(), retval.fortran_vec(), a); } - else if (n < 0) - (*current_liboctave_error_handler) ("rand: invalid negative argument"); return retval; } @@ -693,3 +792,94 @@ return; } + +void +octave_rand::fill (octave_idx_type len, float *v, float a) +{ + if (len < 1) + return; + + switch (current_distribution) + { + case uniform_dist: + if (use_old_generators) + { +#define RAND_FUNC(x) F77_FUNC (dgenunf, DGENUNF) (0.0, 1.0, x) + MAKE_RAND (len); +#undef RAND_FUNC + } + else + oct_fill_float_randu (len, v); + break; + + case normal_dist: + if (use_old_generators) + { +#define RAND_FUNC(x) F77_FUNC (dgennor, DGENNOR) (0.0, 1.0, x) + MAKE_RAND (len); +#undef RAND_FUNC + } + else + oct_fill_float_randn (len, v); + break; + + case expon_dist: + if (use_old_generators) + { +#define RAND_FUNC(x) F77_FUNC (dgenexp, DGENEXP) (1.0, x) + MAKE_RAND (len); +#undef RAND_FUNC + } + else + oct_fill_float_rande (len, v); + break; + + case poisson_dist: + if (use_old_generators) + { + double da = a; + if (da < 0.0 || xisnan(da) || xisinf(da)) +#define RAND_FUNC(x) x = octave_NaN; + MAKE_RAND (len); +#undef RAND_FUNC + else + { + // workaround bug in ignpoi, by calling with different Mu + double tmp; + F77_FUNC (dignpoi, DIGNPOI) (da + 1, tmp); +#define RAND_FUNC(x) F77_FUNC (dignpoi, DIGNPOI) (da, x) + MAKE_RAND (len); +#undef RAND_FUNC + } + } + else + oct_fill_float_randp (a, len, v); + break; + + case gamma_dist: + if (use_old_generators) + { + double da = a; + if (da <= 0.0 || xisnan(da) || xisinf(da)) +#define RAND_FUNC(x) x = octave_NaN; + MAKE_RAND (len); +#undef RAND_FUNC + else +#define RAND_FUNC(x) F77_FUNC (dgengam, DGENGAM) (1.0, da, x) + MAKE_RAND (len); +#undef RAND_FUNC + } + else + oct_fill_float_randg (a, len, v); + break; + + default: + (*current_liboctave_error_handler) + ("rand: invalid distribution ID = %d", current_distribution); + break; + } + + save_state (); + + return; +}
--- a/liboctave/oct-rand.h +++ b/liboctave/oct-rand.h @@ -27,8 +27,8 @@ #include <string> #include "dColVector.h" -#include "dMatrix.h" #include "dNDArray.h" +#include "fNDArray.h" #include "lo-ieee.h" class @@ -136,13 +136,24 @@ return instance_ok () ? instance->do_scalar (a) : octave_NaN; } - // Return a matrix of numbers from the sequence, filled in column - // major order. - static Matrix matrix (octave_idx_type r, octave_idx_type c, double a = 1.0) + // Return the next number from the sequence. + static float float_scalar (float a = 1.0) { - return instance_ok () ? instance->do_matrix (r, c, a) : Matrix (); + return instance_ok () ? instance->do_float_scalar (a) : octave_Float_NaN; } + // Return an array of numbers from the sequence. + static Array<double> vector (octave_idx_type n, double a = 1.0) + { + return instance_ok () ? instance->do_vector (n, a) : Array<double> (); + } + + // Return an array of numbers from the sequence. + static Array<float> float_vector (octave_idx_type n, float a = 1.0) + { + return instance_ok () ? instance->do_float_vector (n, a) : Array<float> (); + } + // Return an N-dimensional array of numbers from the sequence, // filled in column major order. static NDArray nd_array (const dim_vector& dims, double a = 1.0) @@ -150,10 +161,12 @@ return instance_ok () ? instance->do_nd_array (dims, a) : NDArray (); } - // Return an array of numbers from the sequence. - static Array<double> vector (octave_idx_type n, double a = 1.0) + + // Return an N-dimensional array of numbers from the sequence, + // filled in column major order. + static FloatNDArray float_nd_array (const dim_vector& dims, float a = 1.0) { - return instance_ok () ? instance->do_vector (n, a) : Array<double> (); + return instance_ok () ? instance->do_float_nd_array (dims, a) : FloatNDArray (); } private: @@ -220,16 +233,22 @@ // Return the next number from the sequence. double do_scalar (double a = 1.); - // Return a matrix of numbers from the sequence, filled in column - // major order. - Matrix do_matrix (octave_idx_type r, octave_idx_type c, double a = 1.); + // Return the next number from the sequence. + float do_float_scalar (float a = 1.); + + // Return an array of numbers from the sequence. + Array<double> do_vector (octave_idx_type n, double a = 1.); + + // Return an array of numbers from the sequence. + Array<float> do_float_vector (octave_idx_type n, float a = 1.); // Return an N-dimensional array of numbers from the sequence, // filled in column major order. NDArray do_nd_array (const dim_vector& dims, double a = 1.); - // Return an array of numbers from the sequence. - Array<double> do_vector (octave_idx_type n, double a = 1.); + // Return an N-dimensional array of numbers from the sequence, + // filled in column major order. + FloatNDArray do_float_nd_array (const dim_vector& dims, float a = 1.); // Some helper functions. @@ -248,6 +267,8 @@ void switch_to_generator (int dist); void fill (octave_idx_type len, double *v, double a); + + void fill (octave_idx_type len, float *v, float a); }; #endif
--- a/liboctave/randgamma.c +++ b/liboctave/randgamma.c @@ -141,3 +141,59 @@ oct_fill_randg(a,1,&ret); return ret; } + +#undef NAN +#undef RUNI +#undef RNOR +#undef REXP +#define NAN octave_Float_NaN +#define RUNI oct_float_randu() +#define RNOR oct_float_randn() +#define REXP oct_float_rande() + +void +oct_fill_float_randg (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./sqrt(9.*d); + + /* Handle invalid cases */ + if (a <= 0 || INFINITE(a)) + { + for (i=0; i < n; i++) + r[i] = 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 && log(u) >= 0.5*xsq + d*(1-v+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); + } +} + +float +oct_float_randg (float a) +{ + float ret; + oct_fill_float_randg(a,1,&ret); + return ret; +}
--- a/liboctave/randgamma.h +++ b/liboctave/randgamma.h @@ -32,6 +32,9 @@ extern OCTAVE_API double oct_randg (double a); extern OCTAVE_API void oct_fill_randg (double a, octave_idx_type n, double *p); +extern OCTAVE_API float oct_float_randg (float a); +extern OCTAVE_API void oct_fill_float_randg (float a, octave_idx_type n, float *p); + #ifdef __cplusplus } #endif
--- a/liboctave/randmtzig.c +++ b/liboctave/randmtzig.c @@ -64,10 +64,6 @@ http://www.math.keio.ac.jp/matumoto/emt.html email: matumoto@math.keio.ac.jp - * 2006-04-01 David Bateman - * * convert for use in octave, declaring static functions only used - * here and adding oct_ to functions visible externally - * * inverse sense of ALLBITS * 2004-01-19 Paul Kienzle * * comment out main * add init_by_entropy, get_state, set_state @@ -87,6 +83,15 @@ * * 2005-02-23 Paul Kienzle * * fix -DHAVE_X86_32 flag and add -DUSE_X86_32=0|1 for explicit control + * + * 2006-04-01 David Bateman + * * convert for use in octave, declaring static functions only used + * here and adding oct_ to functions visible externally + * * inverse sense of ALLBITS + * + * 2012-05-18 David Bateman + * * Remove randu64 and ALLBIT option + * * Add the single precision generators */ /* @@ -96,9 +101,6 @@ available. This is not necessary if your architecture has /dev/urandom defined. - Compile with -DALLBITS to disable 53-bit random numbers. This is about - 50% slower than using 32-bit random numbers. - Uses implicit -Di386 or explicit -DHAVE_X86_32 to determine if CPU=x86. You can force X86 behaviour with -DUSE_X86_32=1, or suppress it with -DUSE_X86_32=0. You should also consider -march=i686 or similar for @@ -126,21 +128,28 @@ static uint32_t randi32(void) returns 32-bit unsigned int static uint64_t randi53(void) returns 53-bit unsigned int static uint64_t randi54(void) returns 54-bit unsigned int - static uint64_t randi64(void) returns 64-bit unsigned int - static double randu32(void) returns 32-bit uniform in (0,1) + static float randu32(void) returns 32-bit uniform in (0,1) static double randu53(void) returns 53-bit uniform in (0,1) double oct_randu(void) returns M-bit uniform in (0,1) double oct_randn(void) returns M-bit standard normal double oct_rande(void) returns N-bit standard exponential + float oct_float_randu(void) returns M-bit uniform in (0,1) + float oct_float_randn(void) returns M-bit standard normal + float oct_float_rande(void) returns N-bit standard exponential + === Array generators === void oct_fill_randi32(octave_idx_type, uint32_t []) void oct_fill_randi64(octave_idx_type, uint64_t []) + void oct_fill_randu(octave_idx_type, double []) void oct_fill_randn(octave_idx_type, double []) void oct_fill_rande(octave_idx_type, double []) + void oct_fill_float_randu(octave_idx_type, float []) + void oct_fill_float_randn(octave_idx_type, float []) + void oct_fill_float_rande(octave_idx_type, float []) */ #if defined (HAVE_CONFIG_H) @@ -180,6 +189,7 @@ static int left = 1; static int initf = 0; static int initt = 1; +static int inittf = 1; /* initializes state[MT_N] with a seed */ void @@ -378,34 +388,14 @@ #endif } -#if 0 -// FIXME -- this doesn't seem to be used anywhere; should it be removed? -static uint64_t -randi64 (void) -{ - const uint32_t lo = randi32(); - const uint32_t hi = randi32(); -#if HAVE_X86_32 - uint64_t u; - uint32_t *p = (uint32_t *)&u; - p[0] = lo; - p[1] = hi; - return u; -#else - return (((uint64_t)hi<<32)|lo); -#endif -} -#endif - -#ifdef ALLBITS /* generates a random number on (0,1)-real-interval */ -static double +static float randu32 (void) { - return ((double)randi32() + 0.5) * (1.0/4294967296.0); + return ((float)randi32() + 0.5) * (1.0/4294967296.0); /* divided by 2^32 */ } -#else + /* generates a random number on (0,1) with 53-bit resolution */ static double randu53 (void) @@ -414,35 +404,22 @@ const uint32_t b=randi32()>>6; return (a*67108864.0+b+0.4) * (1.0/9007199254740992.0); } -#endif /* Determine mantissa for uniform doubles */ double oct_randu (void) { -#ifdef ALLBITS + return randu53 (); +} + +/* Determine mantissa for uniform floats */ +float +oct_float_randu (void) +{ return randu32 (); -#else - return randu53 (); -#endif } /* ===== Ziggurat normal and exponential generators ===== */ -#ifdef ALLBITS -# define ZIGINT uint32_t -# define EMANTISSA 4294967296.0 /* 32 bit mantissa */ -# define ERANDI randi32() /* 32 bits for mantissa */ -# define NMANTISSA 2147483648.0 /* 31 bit mantissa */ -# define NRANDI randi32() /* 31 bits for mantissa + 1 bit sign */ -# define RANDU randu32() -#else -# define ZIGINT uint64_t -# define EMANTISSA 9007199254740992.0 /* 53 bit mantissa */ -# define ERANDI randi53() /* 53 bits for mantissa */ -# define NMANTISSA EMANTISSA -# define NRANDI randi54() /* 53 bits for mantissa + 1 bit sign */ -# define RANDU randu53() -#endif #define ZIGGURAT_TABLE_SIZE 256 @@ -454,6 +431,13 @@ #define ZIGGURAT_EXP_INV_R 0.129918765548341586 #define EXP_SECTION_AREA 0.0039496598225815571993 +#define ZIGINT uint64_t +#define EMANTISSA 9007199254740992.0 /* 53 bit mantissa */ +#define ERANDI randi53() /* 53 bits for mantissa */ +#define NMANTISSA EMANTISSA +#define NRANDI randi54() /* 53 bits for mantissa + 1 bit sign */ +#define RANDU randu53() + static ZIGINT ki[ZIGGURAT_TABLE_SIZE]; static double wi[ZIGGURAT_TABLE_SIZE], fi[ZIGGURAT_TABLE_SIZE]; static ZIGINT ke[ZIGGURAT_TABLE_SIZE]; @@ -592,7 +576,6 @@ * Of course, different compilers and operating systems may * have something to do with this. */ -#if !defined(ALLBITS) # if HAVE_X86_32 /* 53-bit mantissa, 1-bit sign, x86 32-bit architecture */ double x; @@ -615,14 +598,6 @@ const double x = ( r&1 ? -rabs : rabs) * wi[idx]; # endif /* !HAVE_X86_32 */ if (rabs < (int64_t)ki[idx]) -#else /* ALLBITS */ - /* 32-bit mantissa */ - const uint32_t r = randi32(); - const uint32_t rabs = r&LMASK; - const int idx = (int)(r&0xFF); - const double x = ((int32_t)r) * wi[idx]; - if (rabs < ki[idx]) -#endif /* ALLBITS */ return x; /* 99.3% of the time we return here 1st try */ else if (idx == 0) { @@ -677,6 +652,173 @@ } } +#undef ZIGINT +#undef EMANTISSA +#undef ERANDI +#undef NMANTISSA +#undef NRANDI +#undef RANDU + +#define ZIGINT uint32_t +#define EMANTISSA 4294967296.0 /* 32 bit mantissa */ +#define ERANDI randi32() /* 32 bits for mantissa */ +#define NMANTISSA 2147483648.0 /* 31 bit mantissa */ +#define NRANDI randi32() /* 31 bits for mantissa + 1 bit sign */ +#define RANDU randu32() + +static ZIGINT fki[ZIGGURAT_TABLE_SIZE]; +static float fwi[ZIGGURAT_TABLE_SIZE], ffi[ZIGGURAT_TABLE_SIZE]; +static ZIGINT fke[ZIGGURAT_TABLE_SIZE]; +static float fwe[ZIGGURAT_TABLE_SIZE], ffe[ZIGGURAT_TABLE_SIZE]; + + +static void +create_ziggurat_float_tables (void) +{ + int i; + float x, x1; + + /* Ziggurat tables for the normal distribution */ + x1 = ZIGGURAT_NOR_R; + fwi[255] = x1 / NMANTISSA; + ffi[255] = exp (-0.5 * x1 * x1); + + /* Index zero is special for tail strip, where Marsaglia and Tsang + * defines this as + * k_0 = 2^31 * r * f(r) / v, w_0 = 0.5^31 * v / f(r), f_0 = 1, + * where v is the area of each strip of the ziggurat. + */ + fki[0] = (ZIGINT) (x1 * ffi[255] / NOR_SECTION_AREA * NMANTISSA); + fwi[0] = NOR_SECTION_AREA / ffi[255] / NMANTISSA; + ffi[0] = 1.; + + for (i = 254; i > 0; i--) + { + /* New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus + * need inverse operator of y = exp(-0.5*x*x) -> x = sqrt(-2*ln(y)) + */ + x = sqrt(-2. * log(NOR_SECTION_AREA / x1 + ffi[i+1])); + fki[i+1] = (ZIGINT)(x / x1 * NMANTISSA); + fwi[i] = x / NMANTISSA; + ffi[i] = exp (-0.5 * x * x); + x1 = x; + } + + fki[1] = 0; + + /* Zigurrat tables for the exponential distribution */ + x1 = ZIGGURAT_EXP_R; + fwe[255] = x1 / EMANTISSA; + ffe[255] = exp (-x1); + + /* Index zero is special for tail strip, where Marsaglia and Tsang + * defines this as + * k_0 = 2^32 * r * f(r) / v, w_0 = 0.5^32 * v / f(r), f_0 = 1, + * where v is the area of each strip of the ziggurat. + */ + fke[0] = (ZIGINT) (x1 * ffe[255] / EXP_SECTION_AREA * EMANTISSA); + fwe[0] = EXP_SECTION_AREA / ffe[255] / EMANTISSA; + ffe[0] = 1.; + + for (i = 254; i > 0; i--) + { + /* New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus + * need inverse operator of y = exp(-x) -> x = -ln(y) + */ + x = - log(EXP_SECTION_AREA / x1 + ffe[i+1]); + fke[i+1] = (ZIGINT)(x / x1 * EMANTISSA); + fwe[i] = x / EMANTISSA; + ffe[i] = exp (-x); + x1 = x; + } + fke[1] = 0; + + inittf = 0; +} + +/* + * Here is the guts of the algorithm. As Marsaglia and Tsang state the + * algorithm in their paper + * + * 1) Calculate a random signed integer j and let i be the index + * provided by the rightmost 8-bits of j + * 2) Set x = j * w_i. If j < k_i return x + * 3) If i = 0, then return x from the tail + * 4) If [f(x_{i-1}) - f(x_i)] * U < f(x) - f(x_i), return x + * 5) goto step 1 + * + * Where f is the functional form of the distribution, which for a normal + * distribution is exp(-0.5*x*x) + */ + +float +oct_float_randn (void) +{ + if (inittf) + create_ziggurat_float_tables(); + + while (1) + { + /* 32-bit mantissa */ + const uint32_t r = randi32(); + const uint32_t rabs = r&LMASK; + const int idx = (int)(r&0xFF); + const float x = ((int32_t)r) * fwi[idx]; + if (rabs < fki[idx]) + return x; /* 99.3% of the time we return here 1st try */ + else if (idx == 0) + { + /* As stated in Marsaglia and Tsang + * + * For the normal tail, the method of Marsaglia[5] provides: + * generate x = -ln(U_1)/r, y = -ln(U_2), until y+y > x*x, + * then return r+x. Except that r+x is always in the positive + * tail!!!! Any thing random might be used to determine the + * sign, but as we already have r we might as well use it + * + * [PAK] but not the bottom 8 bits, since they are all 0 here! + */ + float xx, yy; + do + { + xx = - ZIGGURAT_NOR_INV_R * log (RANDU); + yy = - log (RANDU); + } + while ( yy+yy <= xx*xx); + return (rabs&0x100 ? -ZIGGURAT_NOR_R-xx : ZIGGURAT_NOR_R+xx); + } + else if ((ffi[idx-1] - ffi[idx]) * RANDU + ffi[idx] < exp(-0.5*x*x)) + return x; + } +} + +float +oct_float_rande (void) +{ + if (inittf) + create_ziggurat_float_tables(); + + while (1) + { + ZIGINT ri = ERANDI; + const int idx = (int)(ri & 0xFF); + const float x = ri * fwe[idx]; + if (ri < fke[idx]) + return x; // 98.9% of the time we return here 1st try + else if (idx == 0) + { + /* As stated in Marsaglia and Tsang + * + * For the exponential tail, the method of Marsaglia[5] provides: + * x = r - ln(U); + */ + return ZIGGURAT_EXP_R - log(RANDU); + } + else if ((ffe[idx-1] - ffe[idx]) * RANDU + ffe[idx] < exp(-x)) + return x; + } +} + /* Array generators */ void oct_fill_randu (octave_idx_type n, double *p) @@ -701,3 +843,27 @@ for (i = 0; i < n; i++) p[i] = oct_rande(); } + +void +oct_fill_float_randu (octave_idx_type n, float *p) +{ + octave_idx_type i; + for (i = 0; i < n; i++) + p[i] = oct_float_randu(); +} + +void +oct_fill_float_randn (octave_idx_type n, float *p) +{ + octave_idx_type i; + for (i = 0; i < n; i++) + p[i] = oct_float_randn(); +} + +void +oct_fill_float_rande (octave_idx_type n, float *p) +{ + octave_idx_type i; + for (i = 0; i < n; i++) + p[i] = oct_float_rande(); +}
--- a/liboctave/randmtzig.h +++ b/liboctave/randmtzig.h @@ -82,11 +82,19 @@ extern OCTAVE_API double oct_randn (void); extern OCTAVE_API double oct_rande (void); +extern OCTAVE_API float oct_float_randu (void); +extern OCTAVE_API float oct_float_randn (void); +extern OCTAVE_API float oct_float_rande (void); + /* === Array generators === */ extern OCTAVE_API void oct_fill_randu (octave_idx_type n, double *p); extern OCTAVE_API void oct_fill_randn (octave_idx_type n, double *p); extern OCTAVE_API void oct_fill_rande (octave_idx_type n, double *p); +extern OCTAVE_API void oct_fill_float_randu (octave_idx_type n, float *p); +extern OCTAVE_API void oct_fill_float_randn (octave_idx_type n, float *p); +extern OCTAVE_API void oct_fill_float_rande (octave_idx_type n, float *p); + #ifdef __cplusplus } #endif
--- a/liboctave/randpoisson.c +++ b/liboctave/randpoisson.c @@ -368,6 +368,46 @@ } } +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 = (int)floor(lambda); + double P; + int tableidx; + size_t i = n; + + t[0] = P = exp(-lambda); + for (tableidx = 1; tableidx <= intlambda; tableidx++) { + P = P*lambda/(double)tableidx; + t[tableidx] = t[tableidx-1] + P; + } + + while (i-- > 0) { + double u = RUNI; + int k = (u > 0.458 ? intlambda : 0); + nextk: + if ( u <= t[k] ) { + p[i] = (float) k; + continue; + } + if (++k < tableidx) goto nextk; + + while (tableidx < TABLESIZE) { + P = P*lambda/(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] = (float)(tableidx-1); + } +} + /* From Press, et al., Numerical Recipes */ static void poisson_rejection (double lambda, double *p, size_t n) @@ -392,6 +432,30 @@ } } +/* From Press, et al., Numerical Recipes */ +static void +poisson_rejection_float (double lambda, float *p, size_t n) +{ + double sq = sqrt(2.0*lambda); + double alxm = log(lambda); + double g = lambda*alxm - LGAMMA(lambda+1.0); + size_t i; + + for (i = 0; i < n; i++) + { + double y, em, t; + do { + do { + y = tan(M_PI*RUNI); + em = sq * y + lambda; + } while (em < 0.0); + em = floor(em); + t = 0.9*(1.0+y*y)*exp(em*alxm-flogfak(em)-g); + } while (RUNI > t); + p[i] = em; + } +} + /* The cutoff of L <= 1e8 in the following two functions before using * the normal approximation is based on: * > L=1e8; x=floor(linspace(0,2*L,1000)); @@ -463,3 +527,68 @@ } return ret; } + +/* Generate a set of poisson numbers with the same distribution */ +void +oct_fill_float_randp (float FL, octave_idx_type n, float *p) +{ + double L = FL; + octave_idx_type i; + if (L < 0.0 || INFINITE(L)) + { + for (i=0; i<n; i++) + p[i] = 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 = sqrt(L); + for (i = 0; i < n; i++) + { + p[i] = floor(RNOR*sqrtL + L + 0.5); + if (p[i] < 0.0) + p[i] = 0.0; /* will probably never happen */ + } + } +} + +/* Generate one poisson variate */ +float +oct_float_randp (float FL) +{ + double L = FL; + float ret; + if (L < 0.0) ret = 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 *= RUNI; + } while (t > g); + ret = em; + } else if (L <= 1e8) { + /* numerical recipes */ + poisson_rejection_float(L, &ret, 1); + } else if (INFINITE(L)) { + /* FIXME R uses NaN, but the normal approx. suggests that as + * limit should be inf. Which is correct? */ + ret = NAN; + } else { + /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */ + ret = floor(RNOR*sqrt(L) + L + 0.5); + if (ret < 0.0) ret = 0.0; /* will probably never happen */ + } + return ret; +}
--- a/liboctave/randpoisson.h +++ b/liboctave/randpoisson.h @@ -32,6 +32,9 @@ extern OCTAVE_API double oct_randp (double L); extern OCTAVE_API void oct_fill_randp (double L, octave_idx_type n, double *p); +extern OCTAVE_API float oct_float_randp (float L); +extern OCTAVE_API void oct_fill_float_randp (float L, octave_idx_type n, float *p); + #ifdef __cplusplus } #endif
--- a/src/DLD-FUNCTIONS/rand.cc +++ b/src/DLD-FUNCTIONS/rand.cc @@ -1,3 +1,4 @@ + /* Copyright (C) 1996-2012 John W. Eaton @@ -60,6 +61,7 @@ NDArray a; int idx = 0; dim_vector dims; + bool is_single = false; unwind_protect frame; // Restore current distribution on any exit. @@ -68,6 +70,19 @@ octave_rand::distribution (distribution); + if (nargin > 0 && args(nargin-1).is_string()) + { + std::string s_arg = args(nargin-1).string_value (); + + if (s_arg == "single") + { + is_single = true; + nargin--; + } + else if (s_arg == "double") + nargin--; + } + if (additional_arg) { if (nargin == 0) @@ -298,27 +313,54 @@ dims.chop_trailing_singletons (); - if (additional_arg) + if (is_single) { - if (a.length() == 1) - return octave_rand::nd_array (dims, a(0)); - else + if (additional_arg) { - if (a.dims() != dims) + if (a.length() == 1) + return octave_rand::float_nd_array (dims, a(0)); + else { - error ("%s: mismatch in argument size", fcn); - return retval; + if (a.dims() != dims) + { + error ("%s: mismatch in argument size", fcn); + return retval; + } + octave_idx_type len = a.length (); + FloatNDArray m (dims); + float *v = m.fortran_vec (); + for (octave_idx_type i = 0; i < len; i++) + v[i] = octave_rand::float_scalar (a(i)); + return m; } - octave_idx_type len = a.length (); - NDArray m (dims); - double *v = m.fortran_vec (); - for (octave_idx_type i = 0; i < len; i++) - v[i] = octave_rand::scalar (a(i)); - return m; } + else + return octave_rand::float_nd_array (dims); } else - return octave_rand::nd_array (dims); + { + if (additional_arg) + { + if (a.length() == 1) + return octave_rand::nd_array (dims, a(0)); + else + { + if (a.dims() != dims) + { + error ("%s: mismatch in argument size", fcn); + return retval; + } + octave_idx_type len = a.length (); + NDArray m (dims); + double *v = m.fortran_vec (); + for (octave_idx_type i = 0; i < len; i++) + v[i] = octave_rand::scalar (a(i)); + return m; + } + } + else + return octave_rand::nd_array (dims); + } } DEFUN_DLD (rand, args, , @@ -332,6 +374,8 @@ @deftypefnx {Loadable Function} {@var{v} =} rand (\"seed\")\n\ @deftypefnx {Loadable Function} {} rand (\"seed\", @var{v})\n\ @deftypefnx {Loadable Function} {} rand (\"seed\", \"reset\")\n\ +@deftypefnx {Loadable Function} {} rand (@dots{}, \"single\")\n\ +@deftypefnx {Loadable Function} {} rand (@dots{}, \"double\")\n\ Return a matrix with random elements uniformly distributed on the\n\ interval (0, 1). The arguments are handled the same as the arguments\n\ for @code{eye}.\n\ @@ -402,6 +446,9 @@ \n\ The state or seed of the generator can be reset to a new random value\n\ using the \"reset\" keyword.\n\ +\n\ +The class of the value returned can be controlled by a trailing \"double\"\n\ +or \"single\" argument. These are the only valid classes.\n\ @seealso{randn, rande, randg, randp}\n\ @end deftypefn") { @@ -499,6 +546,8 @@ @deftypefnx {Loadable Function} {@var{v} =} randn (\"seed\")\n\ @deftypefnx {Loadable Function} {} randn (\"seed\", @var{v})\n\ @deftypefnx {Loadable Function} {} randn (\"seed\", \"reset\")\n\ +@deftypefnx {Loadable Function} {} randn (@dots{}, \"single\")\n\ +@deftypefnx {Loadable Function} {} randn (@dots{}, \"double\")\n\ Return a matrix with normally distributed random\n\ elements having zero mean and variance one. The arguments are\n\ handled the same as the arguments for @code{rand}.\n\ @@ -506,6 +555,9 @@ By default, @code{randn} uses the Marsaglia and Tsang ``Ziggurat technique''\n\ to transform from a uniform to a normal distribution.\n\ \n\ +The class of the value returned can be controlled by a trailing \"double\"\n\ +or \"single\" argument. These are the only valid classes.\n\ +\n\ Reference: G. Marsaglia and W.W. Tsang,\n\ @cite{Ziggurat Method for Generating Random Variables},\n\ J. Statistical Software, vol 5, 2000,\n\ @@ -565,12 +617,17 @@ @deftypefnx {Loadable Function} {@var{v} =} rande (\"seed\")\n\ @deftypefnx {Loadable Function} {} rande (\"seed\", @var{v})\n\ @deftypefnx {Loadable Function} {} rande (\"seed\", \"reset\")\n\ +@deftypefnx {Loadable Function} {} rande (@dots{}, \"single\")\n\ +@deftypefnx {Loadable Function} {} rande (@dots{}, \"double\")\n\ Return a matrix with exponentially distributed random elements. The\n\ arguments are handled the same as the arguments for @code{rand}.\n\ \n\ By default, @code{randn} uses the Marsaglia and Tsang ``Ziggurat technique''\n\ to transform from a uniform to an exponential distribution.\n\ \n\ +The class of the value returned can be controlled by a trailing \"double\"\n\ +or \"single\" argument. These are the only valid classes.\n\ +\n\ Reference: G. Marsaglia and W.W. Tsang,\n\ @cite{Ziggurat Method for Generating Random Variables},\n\ J. Statistical Software, vol 5, 2000,\n\ @@ -632,6 +689,8 @@ @deftypefnx {Loadable Function} {@var{v} =} randg (\"seed\")\n\ @deftypefnx {Loadable Function} {} randg (\"seed\", @var{v})\n\ @deftypefnx {Loadable Function} {} randg (\"seed\", \"reset\")\n\ +@deftypefnx {Loadable Function} {} randg (@dots{}, \"single\")\n\ +@deftypefnx {Loadable Function} {} randg (@dots{}, \"double\")\n\ Return a matrix with @code{gamma(@var{a},1)} distributed random elements.\n\ The arguments are handled the same as the arguments for @code{rand},\n\ except for the argument @var{a}.\n\ @@ -711,6 +770,9 @@ @end example\n\ \n\ @end table\n\ +\n\ +The class of the value returned can be controlled by a trailing \"double\"\n\ +or \"single\" argument. These are the only valid classes.\n\ @seealso{rand, randn, rande, randp}\n\ @end deftypefn") { @@ -898,6 +960,8 @@ @deftypefnx {Loadable Function} {@var{v} =} randp (\"seed\")\n\ @deftypefnx {Loadable Function} {} randp (\"seed\", @var{v})\n\ @deftypefnx {Loadable Function} {} randp (\"seed\", \"reset\")\n\ +@deftypefnx {Loadable Function} {} randp (@dots{}, \"single\")\n\ +@deftypefnx {Loadable Function} {} randp (@dots{}, \"double\")\n\ Return a matrix with Poisson distributed random elements with mean value\n\ parameter given by the first argument, @var{l}. The arguments\n\ are handled the same as the arguments for @code{rand}, except for the\n\ @@ -928,6 +992,9 @@ L. Montanet, et al., @cite{Review of Particle Properties}, Physical Review\n\ D 50 p1284, 1994.\n\ @end table\n\ +\n\ +The class of the value returned can be controlled by a trailing \"double\"\n\ +or \"single\" argument. These are the only valid classes.\n\ @seealso{rand, randn, rande, randg}\n\ @end deftypefn") {