Skip to content

Commit bd3a3bf

Browse files
LRFLEWLRFLEW
committed
linear_congruential_engine: add using more precision to prevent overflow
1 parent 8c5c4d9 commit bd3a3bf

File tree

6 files changed

+166
-47
lines changed

6 files changed

+166
-47
lines changed

‎libcxx/include/__random/linear_congruential_engine.h

Lines changed: 68 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -30,28 +30,45 @@ template <unsigned long long __a,
3030
unsigned long long __c,
3131
unsigned long long __m,
3232
unsigned long long _Mp,
33-
bool _MightOverflow = (__a != 0 && __m != 0 && __m - 1 > (_Mp - __c) / __a),
34-
bool _OverflowOK = ((__m & (__m - 1)) == 0ull), // m = 2^n
35-
bool _SchrageOK = (__a != 0 && __m != 0 && __m % __a <= __m / __a)> // r <= q
33+
bool _HasOverflow = (__a != 0ull && (__m & (__m - 1ull)) != 0ull), // a != 0, m != 0, m != 2^n
34+
bool _Full = (!_HasOverflow || __m - 1ull <= (_Mp - __c) / __a), // (a * x + c) % m works
35+
bool _Part = (!_HasOverflow || __m - 1ull <= _Mp / __a), // (a * x) % m works
36+
bool _Schrage = (_HasOverflow && __m % __a <= __m / __a)> // r <= q
3637
struct __lce_alg_picker {
37-
static_assert(!_MightOverflow || _OverflowOK || _SchrageOK,
38-
"The current values of a, c, and m cannot generate a number "
39-
"within bounds of linear_congruential_engine.");
38+
static _LIBCPP_CONSTEXPR const int __mode = _Full ? 0 : _Part ? 1 : _Schrage ? 2 : 3;
4039

41-
static _LIBCPP_CONSTEXPR const bool __use_schrage = _MightOverflow && !_OverflowOK && _SchrageOK;
40+
#ifdef _LIBCPP_HAS_NO_INT128
41+
static_assert(_Mp != (unsigned long long)(~0) || _Full || _Part || _Schrage,
42+
"The current values for a, c, and m are not currently supported on platforms without __int128");
43+
#endif
4244
};
4345

4446
template <unsigned long long __a,
4547
unsigned long long __c,
4648
unsigned long long __m,
4749
unsigned long long _Mp,
48-
bool _UseSchrage = __lce_alg_picker<__a, __c, __m, _Mp>::__use_schrage>
50+
int _Mode = __lce_alg_picker<__a, __c, __m, _Mp>::__mode>
4951
struct __lce_ta;
5052

5153
// 64
5254

55+
#ifndef _LIBCPP_HAS_NO_INT128
56+
template <unsigned long long _Ap, unsigned long long _Cp, unsigned long long _Mp>
57+
struct __lce_ta<_Ap, _Cp, _Mp, (unsigned long long)(~0), 3> {
58+
typedef unsigned long long result_type;
59+
_LIBCPP_HIDE_FROM_ABI static result_type next(result_type __xp) {
60+
__extension__ typedef unsigned __int128 calc_type;
61+
const calc_type __a = static_cast<calc_type>(_Ap);
62+
const calc_type __c = static_cast<calc_type>(_Cp);
63+
const calc_type __m = static_cast<calc_type>(_Mp);
64+
const calc_type __x = static_cast<calc_type>(__xp);
65+
return static_cast<result_type>((__a * __x + __c) % __m);
66+
}
67+
};
68+
#endif
69+
5370
template <unsigned long long __a, unsigned long long __c, unsigned long long __m>
54-
struct __lce_ta<__a, __c, __m, (unsigned long long)(~0), true> {
71+
struct __lce_ta<__a, __c, __m, (unsigned long long)(~0), 2> {
5572
typedef unsigned long long result_type;
5673
_LIBCPP_HIDE_FROM_ABI static result_type next(result_type __x) {
5774
// Schrage's algorithm
@@ -66,7 +83,7 @@ struct __lce_ta<__a, __c, __m, (unsigned long long)(~0), true> {
6683
};
6784

6885
template <unsigned long long __a, unsigned long long __m>
69-
struct __lce_ta<__a, 0, __m, (unsigned long long)(~0), true> {
86+
struct __lce_ta<__a, 0ull, __m, (unsigned long long)(~0), 2> {
7087
typedef unsigned long long result_type;
7188
_LIBCPP_HIDE_FROM_ABI static result_type next(result_type __x) {
7289
// Schrage's algorithm
@@ -80,21 +97,40 @@ struct __lce_ta<__a, 0, __m, (unsigned long long)(~0), true> {
8097
};
8198

8299
template <unsigned long long __a, unsigned long long __c, unsigned long long __m>
83-
struct __lce_ta<__a, __c, __m, (unsigned long long)(~0), false> {
100+
struct __lce_ta<__a, __c, __m, (unsigned long long)(~0), 1> {
101+
typedef unsigned long long result_type;
102+
_LIBCPP_HIDE_FROM_ABI static result_type next(result_type __x) {
103+
// Use (((a*x) % m) + c) % m
104+
__x = (__a * __x) % __m;
105+
__x += __c - (__x >= __m - __c) * __m;
106+
return __x;
107+
}
108+
};
109+
110+
template <unsigned long long __a, unsigned long long __c, unsigned long long __m>
111+
struct __lce_ta<__a, __c, __m, (unsigned long long)(~0), 0> {
84112
typedef unsigned long long result_type;
85113
_LIBCPP_HIDE_FROM_ABI static result_type next(result_type __x) { return (__a * __x + __c) % __m; }
86114
};
87115

88116
template <unsigned long long __a, unsigned long long __c>
89-
struct __lce_ta<__a, __c, 0, (unsigned long long)(~0), false> {
117+
struct __lce_ta<__a, __c, 0ull, (unsigned long long)(~0), 0> {
90118
typedef unsigned long long result_type;
91119
_LIBCPP_HIDE_FROM_ABI static result_type next(result_type __x) { return __a * __x + __c; }
92120
};
93121

94122
// 32
95123

124+
template <unsigned long long __a, unsigned long long __c, unsigned long long __m>
125+
struct __lce_ta<__a, __c, __m, unsigned(~0), 3> {
126+
typedef unsigned result_type;
127+
_LIBCPP_HIDE_FROM_ABI static result_type next(result_type __x) {
128+
return static_cast<result_type>(__lce_ta<__a, __c, __m, (unsigned long long)(~0)>::next(__x));
129+
}
130+
};
131+
96132
template <unsigned long long _Ap, unsigned long long _Cp, unsigned long long _Mp>
97-
struct __lce_ta<_Ap, _Cp, _Mp, unsigned(~0), true> {
133+
struct __lce_ta<_Ap, _Cp, _Mp, unsigned(~0), 2> {
98134
typedef unsigned result_type;
99135
_LIBCPP_HIDE_FROM_ABI static result_type next(result_type __x) {
100136
const result_type __a = static_cast<result_type>(_Ap);
@@ -112,7 +148,7 @@ struct __lce_ta<_Ap, _Cp, _Mp, unsigned(~0), true> {
112148
};
113149

114150
template <unsigned long long _Ap, unsigned long long _Mp>
115-
struct __lce_ta<_Ap, 0, _Mp, unsigned(~0), true> {
151+
struct __lce_ta<_Ap, 0ull, _Mp, unsigned(~0), 2> {
116152
typedef unsigned result_type;
117153
_LIBCPP_HIDE_FROM_ABI static result_type next(result_type __x) {
118154
const result_type __a = static_cast<result_type>(_Ap);
@@ -128,7 +164,21 @@ struct __lce_ta<_Ap, 0, _Mp, unsigned(~0), true> {
128164
};
129165

130166
template <unsigned long long _Ap, unsigned long long _Cp, unsigned long long _Mp>
131-
struct __lce_ta<_Ap, _Cp, _Mp, unsigned(~0), false> {
167+
struct __lce_ta<_Ap, _Cp, _Mp, unsigned(~0), 1> {
168+
typedef unsigned result_type;
169+
_LIBCPP_HIDE_FROM_ABI static result_type next(result_type __x) {
170+
const result_type __a = static_cast<result_type>(_Ap);
171+
const result_type __c = static_cast<result_type>(_Cp);
172+
const result_type __m = static_cast<result_type>(_Mp);
173+
// Use (((a*x) % m) + c) % m
174+
__x = (__a * __x) % __m;
175+
__x += __c - (__x >= __m - __c) * __m;
176+
return __x;
177+
}
178+
};
179+
180+
template <unsigned long long _Ap, unsigned long long _Cp, unsigned long long _Mp>
181+
struct __lce_ta<_Ap, _Cp, _Mp, unsigned(~0), 0> {
132182
typedef unsigned result_type;
133183
_LIBCPP_HIDE_FROM_ABI static result_type next(result_type __x) {
134184
const result_type __a = static_cast<result_type>(_Ap);
@@ -139,7 +189,7 @@ struct __lce_ta<_Ap, _Cp, _Mp, unsigned(~0), false> {
139189
};
140190

141191
template <unsigned long long _Ap, unsigned long long _Cp>
142-
struct __lce_ta<_Ap, _Cp, 0, unsigned(~0), false> {
192+
struct __lce_ta<_Ap, _Cp, 0ull, unsigned(~0), 0> {
143193
typedef unsigned result_type;
144194
_LIBCPP_HIDE_FROM_ABI static result_type next(result_type __x) {
145195
const result_type __a = static_cast<result_type>(_Ap);
@@ -150,8 +200,8 @@ struct __lce_ta<_Ap, _Cp, 0, unsigned(~0), false> {
150200

151201
// 16
152202

153-
template <unsigned long long __a, unsigned long long __c, unsigned long long __m, bool __b>
154-
struct __lce_ta<__a, __c, __m, (unsigned short)(~0), __b> {
203+
template <unsigned long long __a, unsigned long long __c, unsigned long long __m, int __mode>
204+
struct __lce_ta<__a, __c, __m, (unsigned short)(~0), __mode> {
155205
typedef unsigned short result_type;
156206
_LIBCPP_HIDE_FROM_ABI static result_type next(result_type __x) {
157207
return static_cast<result_type>(__lce_ta<__a, __c, __m, unsigned(~0)>::next(__x));

‎libcxx/test/std/numerics/rand/rand.eng/rand.eng.lcong/alg.pass.cpp

Lines changed: 50 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -38,12 +38,12 @@ int main(int, char**)
3838

3939
// m might overflow. The overflow is not OK and result will be in bounds
4040
// so we should use Schrage's algorithm
41-
typedef std::linear_congruential_engine<T, (1ull << 32), 0, (1ull << 63) + 1> E2;
41+
typedef std::linear_congruential_engine<T, 0x100000000ull, 0, (1ull << 63) + 1ull> E2;
4242
E2 e2;
4343
// make sure Schrage's algorithm is used (it would be 0s after the first otherwise)
44-
assert(e2() == (1ull << 32));
44+
assert(e2() == 0x100000000ull);
4545
assert(e2() == (1ull << 63) - 1ull);
46-
assert(e2() == (1ull << 63) - (1ull << 33) + 1ull);
46+
assert(e2() == (1ull << 63) - 0x1ffffffffull);
4747
// make sure result is in bounds
4848
assert(e2() < (1ull << 63) + 1);
4949
assert(e2() < (1ull << 63) + 1);
@@ -56,29 +56,62 @@ int main(int, char**)
5656
typedef std::linear_congruential_engine<T, 0x18000001ull, 0x12347ull, (3ull << 56)> E3;
5757
E3 e3;
5858
// make sure Schrage's algorithm is used
59-
assert(e3() == 402727752ull);
60-
assert(e3() == 162159612030764687ull);
61-
assert(e3() == 108176466184989142ull);
59+
assert(e3() == 0x18012348ull);
60+
assert(e3() == 0x2401b4ed802468full);
61+
assert(e3() == 0x18051ec400369d6ull);
6262
// make sure result is in bounds
6363
assert(e3() < (3ull << 56));
6464
assert(e3() < (3ull << 56));
6565
assert(e3() < (3ull << 56));
6666
assert(e3() < (3ull << 56));
6767
assert(e3() < (3ull << 56));
6868

69-
// m will not overflow so we should not use Schrage's algorithm
70-
typedef std::linear_congruential_engine<T, 1ull, 1, (1ull << 48)> E4;
69+
// 32-bit case:
70+
// m might overflow. The overflow is not OK, result will be in bounds,
71+
// and Schrage's algorithm is incompatible here. Need to use 64 bit arithmetic.
72+
typedef std::linear_congruential_engine<unsigned, 0x10009u, 0u, 0x7fffffffu> E4;
7173
E4 e4;
74+
// make sure enough precision is used
75+
assert(e4() == 0x10009u);
76+
assert(e4() == 0x120053u);
77+
assert(e4() == 0xf5030fu);
78+
// make sure result is in bounds
79+
assert(e4() < 0x7fffffffu);
80+
assert(e4() < 0x7fffffffu);
81+
assert(e4() < 0x7fffffffu);
82+
assert(e4() < 0x7fffffffu);
83+
assert(e4() < 0x7fffffffu);
84+
85+
#ifndef _LIBCPP_HAS_NO_INT128
86+
// m might overflow. The overflow is not OK, result will be in bounds,
87+
// and Schrage's algorithm is incompatible here. Need to use 128 bit arithmetic.
88+
typedef std::linear_congruential_engine<T, 0x100000001ull, 0ull, (1ull << 61) - 1ull> E5;
89+
E5 e5;
90+
// make sure enough precision is used
91+
assert(e5() == 0x100000001ull);
92+
assert(e5() == 0x200000009ull);
93+
assert(e5() == 0xb00000019ull);
94+
// make sure result is in bounds
95+
assert(e5() < (1ull << 61) - 1ull);
96+
assert(e5() < (1ull << 61) - 1ull);
97+
assert(e5() < (1ull << 61) - 1ull);
98+
assert(e5() < (1ull << 61) - 1ull);
99+
assert(e5() < (1ull << 61) - 1ull);
100+
#endif
101+
102+
// m will not overflow so we should not use Schrage's algorithm
103+
typedef std::linear_congruential_engine<T, 1ull, 1, (1ull << 48)> E6;
104+
E6 e6;
72105
// make sure the correct algorithm was used
73-
assert(e4() == 2ull);
74-
assert(e4() == 3ull);
75-
assert(e4() == 4ull);
106+
assert(e6() == 2ull);
107+
assert(e6() == 3ull);
108+
assert(e6() == 4ull);
76109
// make sure result is in bounds
77-
assert(e4() < (1ull << 48));
78-
assert(e4() < (1ull << 48));
79-
assert(e4() < (1ull << 48));
80-
assert(e4() < (1ull << 48));
81-
assert(e4() < (1ull << 48));
110+
assert(e6() < (1ull << 48));
111+
assert(e6() < (1ull << 48));
112+
assert(e6() < (1ull << 48));
113+
assert(e6() < (1ull << 48));
114+
assert(e6() < (1ull << 48));
82115

83116
return 0;
84-
}
117+
}

‎libcxx/test/std/numerics/rand/rand.eng/rand.eng.lcong/assign.pass.cpp

Lines changed: 12 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -61,24 +61,33 @@ test()
6161
test1<T, A, 0, M>();
6262
test1<T, A, M - 2, M>();
6363
test1<T, A, M - 1, M>();
64+
}
65+
66+
template <class T>
67+
void test_ext() {
68+
const T M(static_cast<T>(-1));
6469

65-
/*
6670
// Cases where m is odd and m % a > m / a (not implemented)
6771
test1<T, M - 2, 0, M>();
6872
test1<T, M - 2, M - 2, M>();
6973
test1<T, M - 2, M - 1, M>();
7074
test1<T, M - 1, 0, M>();
7175
test1<T, M - 1, M - 2, M>();
7276
test1<T, M - 1, M - 1, M>();
73-
*/
7477
}
7578

7679
int main(int, char**)
7780
{
7881
test<unsigned short>();
82+
test_ext<unsigned short>();
7983
test<unsigned int>();
84+
test_ext<unsigned int>();
8085
test<unsigned long>();
86+
test_ext<unsigned long>();
8187
test<unsigned long long>();
88+
#ifndef _LIBCPP_HAS_NO_INT128
89+
test_ext<unsigned long long>();
90+
#endif
8291

83-
return 0;
92+
return 0;
8493
}

‎libcxx/test/std/numerics/rand/rand.eng/rand.eng.lcong/copy.pass.cpp

Lines changed: 12 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -60,24 +60,33 @@ test()
6060
test1<T, A, 0, M>();
6161
test1<T, A, M - 2, M>();
6262
test1<T, A, M - 1, M>();
63+
}
64+
65+
template <class T>
66+
void test_ext() {
67+
const T M(static_cast<T>(-1));
6368

64-
/*
6569
// Cases where m is odd and m % a > m / a (not implemented)
6670
test1<T, M - 2, 0, M>();
6771
test1<T, M - 2, M - 2, M>();
6872
test1<T, M - 2, M - 1, M>();
6973
test1<T, M - 1, 0, M>();
7074
test1<T, M - 1, M - 2, M>();
7175
test1<T, M - 1, M - 1, M>();
72-
*/
7376
}
7477

7578
int main(int, char**)
7679
{
7780
test<unsigned short>();
81+
test_ext<unsigned short>();
7882
test<unsigned int>();
83+
test_ext<unsigned int>();
7984
test<unsigned long>();
85+
test_ext<unsigned long>();
8086
test<unsigned long long>();
87+
#ifndef _LIBCPP_HAS_NO_INT128
88+
test_ext<unsigned long long>();
89+
#endif
8190

82-
return 0;
91+
return 0;
8392
}

‎libcxx/test/std/numerics/rand/rand.eng/rand.eng.lcong/default.pass.cpp

Lines changed: 12 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -58,24 +58,33 @@ test()
5858
test1<T, A, 0, M>();
5959
test1<T, A, M - 2, M>();
6060
test1<T, A, M - 1, M>();
61+
}
62+
63+
template <class T>
64+
void test_ext() {
65+
const T M(static_cast<T>(-1));
6166

62-
/*
6367
// Cases where m is odd and m % a > m / a (not implemented)
6468
test1<T, M - 2, 0, M>();
6569
test1<T, M - 2, M - 2, M>();
6670
test1<T, M - 2, M - 1, M>();
6771
test1<T, M - 1, 0, M>();
6872
test1<T, M - 1, M - 2, M>();
6973
test1<T, M - 1, M - 1, M>();
70-
*/
7174
}
7275

7376
int main(int, char**)
7477
{
7578
test<unsigned short>();
79+
test_ext<unsigned short>();
7680
test<unsigned int>();
81+
test_ext<unsigned int>();
7782
test<unsigned long>();
83+
test_ext<unsigned long>();
7884
test<unsigned long long>();
85+
#ifndef _LIBCPP_HAS_NO_INT128
86+
test_ext<unsigned long long>();
87+
#endif
7988

80-
return 0;
89+
return 0;
8190
}

0 commit comments

Comments
 (0)