Skip to content

<random>: Replace inaccurate _XLgamma() with UCRT lgamma() #5891

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
Stephan T. Lavavej (StephanTLavavej) merged 1 commit into from
Nov 27, 2025
Merged

<random>: Replace inaccurate _XLgamma() with UCRT lgamma() #5891

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 4 additions & 8 deletions stl/inc/random
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,10 +65,6 @@ _INLINE_VAR constexpr long double _Exp1 = 2.7182818284590452353602874713526625
_INLINE_VAR constexpr long double _Two32 = 4294967296.0L;
_INLINE_VAR constexpr long double _Two31 = 2147483648.0L;

extern "C++" _CRTIMP2_PURE float __CLRCALL_PURE_OR_CDECL _XLgamma(float) noexcept;
extern "C++" _CRTIMP2_PURE double __CLRCALL_PURE_OR_CDECL _XLgamma(double) noexcept;
extern "C++" _CRTIMP2_PURE long double __CLRCALL_PURE_OR_CDECL _XLgamma(long double) noexcept;

_INLINE_VAR constexpr int _Nwords = 4;

template <class _Elem, class _Traits>
Expand Down Expand Up @@ -1926,7 +1922,7 @@ public:
}
_Sqrt = _CSTD sqrt(2.0 * _Mean0);
_Logm = _CSTD log(_Mean0);
_Gx1 = _Mean0 * _Logm - _XLgamma(_Mean0 + 1.0);
_Gx1 = _Mean0 * _Logm - _STD lgamma(_Mean0 + 1.0);
}

_Ty1 _Mean = 0.0;
Expand Down Expand Up @@ -2030,7 +2026,7 @@ private:
}

if (_Nrand_impl<_Ty1>(_Eng)
<= 0.9 * (1.0 + _Yx * _Yx) * _CSTD exp(_Res * _Par0._Logm - _XLgamma(_Res + 1.0) - _Par0._Gx1)) {
<= 0.9 * (1.0 + _Yx * _Yx) * _CSTD exp(_Res * _Par0._Logm - _STD lgamma(_Res + 1.0) - _Par0._Gx1)) {
return _Res;
}
}
Expand Down Expand Up @@ -2086,7 +2082,7 @@ public:
_Px = _Px0;
_Pp = _Px < 0.5 ? _Px : (1.0 - _Px);
_Mean = _Tx * _Pp;
_Gx1 = _XLgamma(_Tx + 1.0);
_Gx1 = _STD lgamma(_Tx + 1.0);
_Sqrt = _CSTD sqrt(2 * _Mean * (1 - _Pp));
_Logp = _CSTD log(_Pp);
_Logp1 = _CSTD log(1.0 - _Pp);
Expand Down Expand Up @@ -2232,7 +2228,7 @@ private:
}
if (_Nrand_impl<_Ty1>(_Eng)
<= 1.2 * _Par0._Sqrt * (1.0 + _Yx * _Yx)
* _CSTD exp(_Par0._Gx1 - _XLgamma(_Res + 1.0) - _XLgamma(_Par0._Tx - _Res + 1.0)
* _CSTD exp(_Par0._Gx1 - _STD lgamma(_Res + 1.0) - _STD lgamma(_Par0._Tx - _Res + 1.0)
+ _Res * _Par0._Logp + (_Par0._Tx - _Res) * _Par0._Logp1)) {
break;
}
Expand Down
60 changes: 9 additions & 51 deletions stl/src/xlgamma.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,63 +1,21 @@
// Copyright (c) Microsoft Corporation.
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception

// moderately accurate lgamma function for random

// TRANSITION, ABI: This should be superseded by the CRT's lgamma().

#include <cmath>

// #include <random>
_STD_BEGIN
_CRTIMP2_PURE float __CLRCALL_PURE_OR_CDECL _XLgamma(float) noexcept;
_CRTIMP2_PURE double __CLRCALL_PURE_OR_CDECL _XLgamma(double) noexcept;
_CRTIMP2_PURE long double __CLRCALL_PURE_OR_CDECL _XLgamma(long double) noexcept;

float __CLRCALL_PURE_OR_CDECL _XLgamma(float x) noexcept { // moderately accurate log gamma
static const float coeff[6] = {76.18009172947146F, -86.50532032941677F, 24.01409824083091F, -1.23173972450155F,
0.1208650973866179E-2F, -0.5395239384953E-5F};

float val0 = x + 5.5F;
val0 -= (x + 0.5F) * _STD log(val0);
float val1 = 1.000000000190015F;
float y = x + 1.0F;

for (int i = 0; i < 6; ++i, y += 1.0) {
val1 += coeff[i] / y;
}

return -val0 + _STD log(2.5066282746310005F * val1 / x);
// TRANSITION, ABI: preserved for binary compatibility
_CRTIMP2_PURE float __CLRCALL_PURE_OR_CDECL _XLgamma(float x) noexcept {
return _STD lgamma(x);
}

double __CLRCALL_PURE_OR_CDECL _XLgamma(double x) noexcept { // moderately accurate log gamma
static const double coeff[6] = {76.18009172947146, -86.50532032941677, 24.01409824083091, -1.23173972450155,
0.1208650973866179E-2, -0.5395239384953E-5};

double val0 = x + 5.5;
val0 -= (x + 0.5) * _STD log(val0);
double val1 = 1.000000000190015;
double y = x + 1;

for (int i = 0; i < 6; ++i, y += 1.0) {
val1 += coeff[i] / y;
}

return -val0 + _STD log(2.5066282746310005 * val1 / x);
// TRANSITION, ABI: preserved for binary compatibility
_CRTIMP2_PURE double __CLRCALL_PURE_OR_CDECL _XLgamma(double x) noexcept {
return _STD lgamma(x);
}

long double __CLRCALL_PURE_OR_CDECL _XLgamma(long double x) noexcept { // moderately accurate log gamma
static const long double coeff[6] = {76.18009172947146, -86.50532032941677, 24.01409824083091, -1.23173972450155,
0.1208650973866179E-2, -0.5395239384953E-5};

long double val0 = x + 5.5;
val0 -= (x + 0.5) * _STD log(val0);
long double val1 = 1.000000000190015;
long double y = x + 1;

for (int i = 0; i < 6; ++i, y += 1.0) {
val1 += coeff[i] / y;
}

return -val0 + _STD log(2.5066282746310005 * val1 / x);
// TRANSITION, ABI: preserved for binary compatibility
_CRTIMP2_PURE long double __CLRCALL_PURE_OR_CDECL _XLgamma(long double x) noexcept {
return _STD lgamma(x);
}
_STD_END