Finally had a chance to look at this and yes, there is a bug, but the reason is interesting. The random number generator creates 32-bit numbers internally. The std::uniform_real_distribution creates 64-bit doubles. So for each draw of a double, it need to ask the generator for 2 32-but numbers. So to discard `ndiscard` doubles, the randomDouble function should have discarded `2*ndiscard` random numbers.
Forum