diff options
author | Nobuyoshi Nakada <nobu@ruby-lang.org> | 2020-05-04 00:00:27 +0900 |
---|---|---|
committer | Nobuyoshi Nakada <nobu@ruby-lang.org> | 2020-05-04 00:00:27 +0900 |
commit | 5b28f01d77aaf02dc3717c986b8b90f8e99b7f5a (patch) | |
tree | e08a2cbb37912e675c254c7caea656dfe68988db /random.c | |
parent | 5aaa75e7c1f4b7912c10ffdcb1cac581e20eda39 (diff) |
Make int-pair-to-real conversion more portable
And utilize more bits even if DBL_MANT_DIG > 53.
Diffstat (limited to 'random.c')
-rw-r--r-- | random.c | 21 |
1 files changed, 16 insertions, 5 deletions
@@ -14,6 +14,7 @@ #include <errno.h> #include <limits.h> #include <math.h> +#include <float.h> #include <time.h> #ifdef HAVE_UNISTD_H @@ -75,12 +76,22 @@ genrand_real(struct MT *mt) return int_pair_to_real_exclusive(a, b); } +static const double dbl_reduce_scale = /* 2**(-DBL_MANT_DIG) */ + (1.0 + / (double)(DBL_MANT_DIG > 2*31 ? (1ul<<31) : 1.0) + / (double)(DBL_MANT_DIG > 1*31 ? (1ul<<31) : 1.0) + / (double)(1ul<<(DBL_MANT_DIG%31))); + static double int_pair_to_real_exclusive(uint32_t a, uint32_t b) { - a >>= 5; - b >>= 6; - return(a*67108864.0+b)*(1.0/9007199254740992.0); + static const int a_shift = DBL_MANT_DIG < 64 ? + (64-DBL_MANT_DIG)/2 : 0; + static const int b_shift = DBL_MANT_DIG < 64 ? + (64-DBL_MANT_DIG)-a_shift : 0; + a >>= a_shift; + b >>= b_shift; + return (a*(double)(1ul<<(32-b_shift))+b)*dbl_reduce_scale; } /* generates a random number on [0,1] with 53-bit resolution*/ @@ -148,7 +159,7 @@ static double int_pair_to_real_inclusive(uint32_t a, uint32_t b) { double r; - enum {dig = 53}; + enum {dig = DBL_MANT_DIG}; enum {dig_u = dig-32, dig_r64 = 64-dig, bmask = ~(~0u<<(dig_r64))}; #if defined HAVE_UINT128_T const uint128_t m = ((uint128_t)1 << dig) | 1; @@ -163,7 +174,7 @@ int_pair_to_real_inclusive(uint32_t a, uint32_t b) b = (b >> dig_r64) + (((a >> dig_u) + (b & bmask)) >> dig_r64); r = (double)a * (1 << dig_u) + b; #endif - return ldexp(r, -dig); + return r * dbl_reduce_scale; } VALUE rb_cRandom; |