<random>: binomial_distribution is broken

[statementreply]

Describe the bug

binomial_distribution doesn't produce the correct random number distribution.

Command-line test case

G:\Temp>type devcom_606865.cpp
#include <cstdio>
#include <random>

using namespace std;

int main() {
    constexpr int num = 1000;
    constexpr double prob = 0.001238;
    constexpr int iters = 100'000;
    constexpr int seed = 12345;

    mt19937 mt_rand(seed);
    binomial_distribution<int> distribution(num, prob);

    double sum = 0;
    for (int i = 0; i < iters; i++) {
        sum += distribution(mt_rand);
    }

    printf("Actual average            = %f\n", sum / iters);
    printf("Mathematical expectiation = %f\n", prob * num);
}

G:\Temp>cl /EHsc /W4 /WX devcom_606865.cpp
用于 x64 的 Microsoft (R) C/C++ 优化编译器 19.27.29009.1 版
版权所有(C) Microsoft Corporation。保留所有权利。

devcom_606865.cpp
Microsoft (R) Incremental Linker Version 14.27.29009.1
Copyright (C) Microsoft Corporation.  All rights reserved.

/out:devcom_606865.exe
devcom_606865.obj

G:\Temp>.\devcom_606865.exe
Actual average            = 0.996960
Mathematical expectiation = 1.238000

Expected behavior

The average value should be close to the mathematical expectation.

STL version

https://github.com/microsoft/STL/commit/5be7d49c243979231dff35919d3a3d813d75d714

Additional context

It appears that the code in <random> does not exactly follow the referenced Numerical Recipes function.

The <random> header code casts the result of the (_Par0._Sqrt * _Yx + _Par0._Mean) calculation to an int, before performing the comparison:

STL/stl/inc/random

Lines 2497 to 2502 in 5be7d49

	for (;;) { // generate a tentative value
	_Yx = static_cast<_Ty1>(_CSTD tan(_Pi * _NRAND(_Eng, _Ty1)));
	_Res = static_cast<_Ty>(_Par0._Sqrt * _Yx + _Par0._Mean);
	if (_Ty{0} <= _Res && _Res <= _Par0._Tx) {
	break;
	}

The Numerical Recipes bnldev function calculates the floor after performing the comparison (the variable em is declared as a double):

do {
	angle=PI*ran1(idum);
	y=tan(angle);
	em=sq*y+am;
} while (em < 0.0 || em >= (en+1.0));
em=floor(em);

When the original issue reporter modified the <random> header as follows:

for (;;) { // generate a tentative value
    _Yx = (_Ty1)_CSTD tan(_Pi * _NRAND(_Eng, _Ty1));
    double val = (_Par0._Sqrt * _Yx + _Par0._Mean);
    _Res = (_Ty)(val);
    if (static_cast<_Ty>(0) <= val && val <= _Par0._Tx)
        break;
}

(sic, unfortunately still incorrect) and recompile, she/he got a more reasonable result.

Also tracked by DevCom-606865 and Microsoft-internal VSO-923008 / AB#923008.

[statementreply]

It appears that there's a similar bug in poisson_distribution.