Mercurial > octave
changeset 27855:9405e2be91d0
Restrict rand ("single") to range (0, 1) (bug #41742).
* NEWS: Announce change in rand function.
* randmtzig.cc (randu24): Rename function from randu32.
New algorithm uses masking to reduce randi32() to 24-bits.
Test for value of 0 to exclude left endpoint of range.
Divide result by 2^24 which guarantees right endpoint of 1 is excluded.
* randmtzig.cc: Rename all instances of randu32 with randu24.
| author | Rik <rik@octave.org> |
|---|---|
| date | Thu, 19 Dec 2019 11:57:54 -0800 |
| parents | a57e88065f89 |
| children | b3a03dbde858 |
| files | NEWS liboctave/numeric/randmtzig.cc |
| diffstat | 2 files changed, 22 insertions(+), 6 deletions(-) [+] |
line wrap: on
line diff
--- a/NEWS Thu Dec 19 11:29:22 2019 -0800 +++ b/NEWS Thu Dec 19 11:57:54 2019 -0800 @@ -18,6 +18,15 @@ endpoints are symmetric. This is more intuitive and also compatible with recent changes made in Matlab R2019b. +- The underlying algorithm of the `rand` function has been changed. + For single precision outputs, the algorithm has been fixed so that it + produces values strictly in the range (0, 1). Previously, it could + occasionally generate the right endpoint value of 1 (See bug #41742). + In addition, the new implementation uses a uniform interval between + floating point values in the range (0, 1) rather than targeting a + uniform density (# of random integers / length along real number + line). + - The `edit` function option `"editinplace"` now defaults to `true` and the option `"home"` now defaults to the empty matrix `[]`. Files will no longer be copied to the user's HOME directory for editing. The old
--- a/liboctave/numeric/randmtzig.cc Thu Dec 19 11:29:22 2019 -0800 +++ b/liboctave/numeric/randmtzig.cc Thu Dec 19 11:57:54 2019 -0800 @@ -138,7 +138,7 @@ 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 float randu32 (void) returns 32-bit uniform in (0,1) + static float randu24 (void) returns 24-bit uniform in (0,1) static double randu53 (void) returns 53-bit uniform in (0,1) double rand_uniform (void) returns M-bit uniform in (0,1) @@ -377,10 +377,17 @@ } /* generates a random number on (0,1)-real-interval */ - static float randu32 (void) + static float randu24 (void) { - return (static_cast<float> (randi32 ()) + 0.5) * (1.0/4294967296.0); - /* divided by 2^32 */ + uint32_t i; + + do + { + i = randi32 () & static_cast<uint32_t> (0xFFFFFF); + } + while (i == 0); + + return i * (1.0f / 16777216.0f); } /* generates a random number on (0,1) with 53-bit resolution */ @@ -404,7 +411,7 @@ float rand_uniform<float> (void) { - return randu32 (); + return randu24 (); } /* ===== Ziggurat normal and exponential generators ===== */ @@ -665,7 +672,7 @@ #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() +#define RANDU randu24() static ZIGINT fki[ZIGGURAT_TABLE_SIZE]; static float fwi[ZIGGURAT_TABLE_SIZE], ffi[ZIGGURAT_TABLE_SIZE];