From 5b28f01d77aaf02dc3717c986b8b90f8e99b7f5a Mon Sep 17 00:00:00 2001 From: Nobuyoshi Nakada Date: Mon, 4 May 2020 00:00:27 +0900 Subject: Make int-pair-to-real conversion more portable And utilize more bits even if DBL_MANT_DIG > 53. --- random.c | 21 ++++++++++++++++----- 1 file changed, 16 insertions(+), 5 deletions(-) diff --git a/random.c b/random.c index e5676656cc..de2767a93d 100644 --- a/random.c +++ b/random.c @@ -14,6 +14,7 @@ #include #include #include +#include #include #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; -- cgit v1.2.3