From d5fcbe966bc501db8bf6a3809633f0b82e6ae547 Mon Sep 17 00:00:00 2001 From: Pieter Wuille Date: Tue, 18 Jun 2024 12:39:56 -0400 Subject: [PATCH] random: improve precision of MakeExponentiallyDistributed --- src/random.cpp | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/src/random.cpp b/src/random.cpp index 530dfff1ed7..cb5c127e0da 100644 --- a/src/random.cpp +++ b/src/random.cpp @@ -775,5 +775,14 @@ void RandomInit() double MakeExponentiallyDistributed(uint64_t uniform) noexcept { - return -std::log1p((uniform >> 16) * -0.0000000000000035527136788 /* -1/2^48 */); + // To convert uniform into an exponentially-distributed double, we use two steps: + // - Convert uniform into a uniformly-distributed double in range [0, 1), use the expression + // ((uniform >> 11) * 0x1.0p-53), as described in https://prng.di.unimi.it/ under + // "Generating uniform doubles in the unit interval". Call this value x. + // - Given an x in uniformly distributed in [0, 1), we find an exponentially distributed value + // by applying the quantile function to it. For the exponential distribution with mean 1 this + // is F(x) = -log(1 - x). + // + // Combining the two, and using log1p(x) = log(1 + x), we obtain the following: + return -std::log1p((uniform >> 11) * -0x1.0p-53); }