1	// random number generation (out of line) -*- C++ -*-
2
3	// Copyright (C) 2009-2021 Free Software Foundation, Inc.
4	//
5	// This file is part of the GNU ISO C++ Library. This library is free
6	// software; you can redistribute it and/or modify it under the
7	// terms of the GNU General Public License as published by the
8	// Free Software Foundation; either version 3, or (at your option)
9	// any later version.
10
11	// This library is distributed in the hope that it will be useful,
12	// but WITHOUT ANY WARRANTY; without even the implied warranty of
13	// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14	// GNU General Public License for more details.
15
16	// Under Section 7 of GPL version 3, you are granted additional
17	// permissions described in the GCC Runtime Library Exception, version
18	// 3.1, as published by the Free Software Foundation.
19
20	// You should have received a copy of the GNU General Public License and
21	// a copy of the GCC Runtime Library Exception along with this program;
22	// see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
23	// <http://www.gnu.org/licenses/>.
24
25	/** @file bits/random.tcc
26	* This is an internal header file, included by other library headers.
27	* Do not attempt to use it directly. @headername{random}
28	*/
29
30	#ifndef _RANDOM_TCC
31	#define _RANDOM_TCC 1
32
33	#include <numeric> // std::accumulate and std::partial_sum
34
35	namespace std _GLIBCXX_VISIBILITY(default)
36	{
37	_GLIBCXX_BEGIN_NAMESPACE_VERSION
38
39	/// @cond undocumented
40	// (Further) implementation-space details.
41	namespace __detail
42	{
43	// General case for x = (ax + c) mod m -- use Schrage's algorithm
44	// to avoid integer overflow.
45	//
46	// Preconditions: a > 0, m > 0.
47	//
48	// Note: only works correctly for __m % __a < __m / __a.
49	template<typename _Tp, _Tp __m, _Tp __a, _Tp __c>
50	_Tp
51	_Mod<_Tp, __m, __a, __c, false, true>::
52	__calc(_Tp __x)
53	{
54	if (__a == 1)
55	__x %= __m;
56	else
57	{
58	static const _Tp __q = __m / __a;
59	static const _Tp __r = __m % __a;
60
61	_Tp __t1 = __a * (__x % __q);
62	_Tp __t2 = __r * (__x / __q);
63	if (__t1 >= __t2)
64	__x = __t1 - __t2;
65	else
66	__x = __m - __t2 + __t1;
67	}
68
69	if (__c != 0)
70	{
71	const _Tp __d = __m - __x;
72	if (__d > __c)
73	__x += __c;
74	else
75	__x = __c - __d;
76	}
77	return __x;
78	}
79
80	template<typename _InputIterator, typename _OutputIterator,
81	typename _Tp>
82	_OutputIterator
83	__normalize(_InputIterator __first, _InputIterator __last,
84	_OutputIterator __result, const _Tp& __factor)
85	{
86	for (; __first != __last; ++__first, ++__result)
87	*__result = *__first / __factor;
88	return __result;
89	}
90
91	} // namespace __detail
92	/// @endcond
93
94	template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
95	constexpr _UIntType
96	linear_congruential_engine<_UIntType, __a, __c, __m>::multiplier;
97
98	template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
99	constexpr _UIntType
100	linear_congruential_engine<_UIntType, __a, __c, __m>::increment;
101
102	template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
103	constexpr _UIntType
104	linear_congruential_engine<_UIntType, __a, __c, __m>::modulus;
105
106	template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
107	constexpr _UIntType
108	linear_congruential_engine<_UIntType, __a, __c, __m>::default_seed;
109
110	/**
111	* Seeds the LCR with integral value @p __s, adjusted so that the
112	* ring identity is never a member of the convergence set.
113	*/
114	template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
115	void
116	linear_congruential_engine<_UIntType, __a, __c, __m>::
117	seed(result_type __s)
118	{
119	if ((__detail::__mod<_UIntType, __m>(__c) == 0)
120	&& (__detail::__mod<_UIntType, __m>(__s) == 0))
121	_M_x = 1;
122	else
123	_M_x = __detail::__mod<_UIntType, __m>(__s);
124	}
125
126	/**
127	* Seeds the LCR engine with a value generated by @p __q.
128	*/
129	template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
130	template<typename _Sseq>
131	auto
132	linear_congruential_engine<_UIntType, __a, __c, __m>::
133	seed(_Sseq& __q)
134	-> _If_seed_seq<_Sseq>
135	{
136	const _UIntType __k0 = __m == 0 ? std::numeric_limits<_UIntType>::digits
137	: std::__lg(__m);
138	const _UIntType __k = (__k0 + 31) / 32;
139	uint_least32_t __arr[__k + 3];
140	__q.generate(__arr + 0, __arr + __k + 3);
141	_UIntType __factor = 1u;
142	_UIntType __sum = 0u;
143	for (size_t __j = 0; __j < __k; ++__j)
144	{
145	__sum += __arr[__j + 3] * __factor;
146	__factor *= __detail::_Shift<_UIntType, 32>::__value;
147	}
148	seed(__sum);
149	}
150
151	template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
152	typename _CharT, typename _Traits>
153	std::basic_ostream<_CharT, _Traits>&
154	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
155	const linear_congruential_engine<_UIntType,
156	__a, __c, __m>& __lcr)
157	{
158	using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
159
160	const typename __ios_base::fmtflags __flags = __os.flags();
161	const _CharT __fill = __os.fill();
162	__os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
163	__os.fill(__os.widen(' '));
164
165	__os << __lcr._M_x;
166
167	__os.flags(__flags);
168	__os.fill(__fill);
169	return __os;
170	}
171
172	template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
173	typename _CharT, typename _Traits>
174	std::basic_istream<_CharT, _Traits>&
175	operator>>(std::basic_istream<_CharT, _Traits>& __is,
176	linear_congruential_engine<_UIntType, __a, __c, __m>& __lcr)
177	{
178	using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
179
180	const typename __ios_base::fmtflags __flags = __is.flags();
181	__is.flags(__ios_base::dec);
182
183	__is >> __lcr._M_x;
184
185	__is.flags(__flags);
186	return __is;
187	}
188
189
190	template<typename _UIntType,
191	size_t __w, size_t __n, size_t __m, size_t __r,
192	_UIntType __a, size_t __u, _UIntType __d, size_t __s,
193	_UIntType __b, size_t __t, _UIntType __c, size_t __l,
194	_UIntType __f>
195	constexpr size_t
196	mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
197	__s, __b, __t, __c, __l, __f>::word_size;
198
199	template<typename _UIntType,
200	size_t __w, size_t __n, size_t __m, size_t __r,
201	_UIntType __a, size_t __u, _UIntType __d, size_t __s,
202	_UIntType __b, size_t __t, _UIntType __c, size_t __l,
203	_UIntType __f>
204	constexpr size_t
205	mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
206	__s, __b, __t, __c, __l, __f>::state_size;
207
208	template<typename _UIntType,
209	size_t __w, size_t __n, size_t __m, size_t __r,
210	_UIntType __a, size_t __u, _UIntType __d, size_t __s,
211	_UIntType __b, size_t __t, _UIntType __c, size_t __l,
212	_UIntType __f>
213	constexpr size_t
214	mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
215	__s, __b, __t, __c, __l, __f>::shift_size;
216
217	template<typename _UIntType,
218	size_t __w, size_t __n, size_t __m, size_t __r,
219	_UIntType __a, size_t __u, _UIntType __d, size_t __s,
220	_UIntType __b, size_t __t, _UIntType __c, size_t __l,
221	_UIntType __f>
222	constexpr size_t
223	mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
224	__s, __b, __t, __c, __l, __f>::mask_bits;
225
226	template<typename _UIntType,
227	size_t __w, size_t __n, size_t __m, size_t __r,
228	_UIntType __a, size_t __u, _UIntType __d, size_t __s,
229	_UIntType __b, size_t __t, _UIntType __c, size_t __l,
230	_UIntType __f>
231	constexpr _UIntType
232	mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
233	__s, __b, __t, __c, __l, __f>::xor_mask;
234
235	template<typename _UIntType,
236	size_t __w, size_t __n, size_t __m, size_t __r,
237	_UIntType __a, size_t __u, _UIntType __d, size_t __s,
238	_UIntType __b, size_t __t, _UIntType __c, size_t __l,
239	_UIntType __f>
240	constexpr size_t
241	mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
242	__s, __b, __t, __c, __l, __f>::tempering_u;
243
244	template<typename _UIntType,
245	size_t __w, size_t __n, size_t __m, size_t __r,
246	_UIntType __a, size_t __u, _UIntType __d, size_t __s,
247	_UIntType __b, size_t __t, _UIntType __c, size_t __l,
248	_UIntType __f>
249	constexpr _UIntType
250	mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
251	__s, __b, __t, __c, __l, __f>::tempering_d;
252
253	template<typename _UIntType,
254	size_t __w, size_t __n, size_t __m, size_t __r,
255	_UIntType __a, size_t __u, _UIntType __d, size_t __s,
256	_UIntType __b, size_t __t, _UIntType __c, size_t __l,
257	_UIntType __f>
258	constexpr size_t
259	mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
260	__s, __b, __t, __c, __l, __f>::tempering_s;
261
262	template<typename _UIntType,
263	size_t __w, size_t __n, size_t __m, size_t __r,
264	_UIntType __a, size_t __u, _UIntType __d, size_t __s,
265	_UIntType __b, size_t __t, _UIntType __c, size_t __l,
266	_UIntType __f>
267	constexpr _UIntType
268	mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
269	__s, __b, __t, __c, __l, __f>::tempering_b;
270
271	template<typename _UIntType,
272	size_t __w, size_t __n, size_t __m, size_t __r,
273	_UIntType __a, size_t __u, _UIntType __d, size_t __s,
274	_UIntType __b, size_t __t, _UIntType __c, size_t __l,
275	_UIntType __f>
276	constexpr size_t
277	mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
278	__s, __b, __t, __c, __l, __f>::tempering_t;
279
280	template<typename _UIntType,
281	size_t __w, size_t __n, size_t __m, size_t __r,
282	_UIntType __a, size_t __u, _UIntType __d, size_t __s,
283	_UIntType __b, size_t __t, _UIntType __c, size_t __l,
284	_UIntType __f>
285	constexpr _UIntType
286	mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
287	__s, __b, __t, __c, __l, __f>::tempering_c;
288
289	template<typename _UIntType,
290	size_t __w, size_t __n, size_t __m, size_t __r,
291	_UIntType __a, size_t __u, _UIntType __d, size_t __s,
292	_UIntType __b, size_t __t, _UIntType __c, size_t __l,
293	_UIntType __f>
294	constexpr size_t
295	mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
296	__s, __b, __t, __c, __l, __f>::tempering_l;
297
298	template<typename _UIntType,
299	size_t __w, size_t __n, size_t __m, size_t __r,
300	_UIntType __a, size_t __u, _UIntType __d, size_t __s,
301	_UIntType __b, size_t __t, _UIntType __c, size_t __l,
302	_UIntType __f>
303	constexpr _UIntType
304	mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
305	__s, __b, __t, __c, __l, __f>::
306	initialization_multiplier;
307
308	template<typename _UIntType,
309	size_t __w, size_t __n, size_t __m, size_t __r,
310	_UIntType __a, size_t __u, _UIntType __d, size_t __s,
311	_UIntType __b, size_t __t, _UIntType __c, size_t __l,
312	_UIntType __f>
313	constexpr _UIntType
314	mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
315	__s, __b, __t, __c, __l, __f>::default_seed;
316
317	template<typename _UIntType,
318	size_t __w, size_t __n, size_t __m, size_t __r,
319	_UIntType __a, size_t __u, _UIntType __d, size_t __s,
320	_UIntType __b, size_t __t, _UIntType __c, size_t __l,
321	_UIntType __f>
322	void
323	mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
324	__s, __b, __t, __c, __l, __f>::
325	seed(result_type __sd)
326	{
327	_M_x[0] = __detail::__mod<_UIntType,
328	__detail::_Shift<_UIntType, __w>::__value>(__sd);
329
330	for (size_t __i = 1; __i < state_size; ++__i)
331	{
332	_UIntType __x = _M_x[__i - 1];
333	__x ^= __x >> (__w - 2);
334	__x *= __f;
335	__x += __detail::__mod<_UIntType, __n>(__i);
336	_M_x[__i] = __detail::__mod<_UIntType,
337	__detail::_Shift<_UIntType, __w>::__value>(__x);
338	}
339	_M_p = state_size;
340	}
341
342	template<typename _UIntType,
343	size_t __w, size_t __n, size_t __m, size_t __r,
344	_UIntType __a, size_t __u, _UIntType __d, size_t __s,
345	_UIntType __b, size_t __t, _UIntType __c, size_t __l,
346	_UIntType __f>
347	template<typename _Sseq>
348	auto
349	mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
350	__s, __b, __t, __c, __l, __f>::
351	seed(_Sseq& __q)
352	-> _If_seed_seq<_Sseq>
353	{
354	const _UIntType __upper_mask = (~_UIntType()) << __r;
355	const size_t __k = (__w + 31) / 32;
356	uint_least32_t __arr[__n * __k];
357	__q.generate(__arr + 0, __arr + __n * __k);
358
359	bool __zero = true;
360	for (size_t __i = 0; __i < state_size; ++__i)
361	{
362	_UIntType __factor = 1u;
363	_UIntType __sum = 0u;
364	for (size_t __j = 0; __j < __k; ++__j)
365	{
366	__sum += __arr[__k * __i + __j] * __factor;
367	__factor *= __detail::_Shift<_UIntType, 32>::__value;
368	}
369	_M_x[__i] = __detail::__mod<_UIntType,
370	__detail::_Shift<_UIntType, __w>::__value>(__sum);
371
372	if (__zero)
373	{
374	if (__i == 0)
375	{
376	if ((_M_x[0] & __upper_mask) != 0u)
377	__zero = false;
378	}
379	else if (_M_x[__i] != 0u)
380	__zero = false;
381	}
382	}
383	if (__zero)
384	_M_x[0] = __detail::_Shift<_UIntType, __w - 1>::__value;
385	_M_p = state_size;
386	}
387
388	template<typename _UIntType, size_t __w,
389	size_t __n, size_t __m, size_t __r,
390	_UIntType __a, size_t __u, _UIntType __d, size_t __s,
391	_UIntType __b, size_t __t, _UIntType __c, size_t __l,
392	_UIntType __f>
393	void
394	mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
395	__s, __b, __t, __c, __l, __f>::
396	_M_gen_rand(void)
397	{
398	const _UIntType __upper_mask = (~_UIntType()) << __r;
399	const _UIntType __lower_mask = ~__upper_mask;
400
401	for (size_t __k = 0; __k < (__n - __m); ++__k)
402	{
403	_UIntType __y = ((_M_x[__k] & __upper_mask)
404	| (_M_x[__k + 1] & __lower_mask));
405	_M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1)
406	^ ((__y & 0x01) ? __a : 0));
407	}
408
409	for (size_t __k = (__n - __m); __k < (__n - 1); ++__k)
410	{
411	_UIntType __y = ((_M_x[__k] & __upper_mask)
412	| (_M_x[__k + 1] & __lower_mask));
413	_M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1)
414	^ ((__y & 0x01) ? __a : 0));
415	}
416
417	_UIntType __y = ((_M_x[__n - 1] & __upper_mask)
418	| (_M_x[0] & __lower_mask));
419	_M_x[__n - 1] = (_M_x[__m - 1] ^ (__y >> 1)
420	^ ((__y & 0x01) ? __a : 0));
421	_M_p = 0;
422	}
423
424	template<typename _UIntType, size_t __w,
425	size_t __n, size_t __m, size_t __r,
426	_UIntType __a, size_t __u, _UIntType __d, size_t __s,
427	_UIntType __b, size_t __t, _UIntType __c, size_t __l,
428	_UIntType __f>
429	void
430	mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
431	__s, __b, __t, __c, __l, __f>::
432	discard(unsigned long long __z)
433	{
434	while (__z > state_size - _M_p)
435	{
436	__z -= state_size - _M_p;
437	_M_gen_rand();
438	}
439	_M_p += __z;
440	}
441
442	template<typename _UIntType, size_t __w,
443	size_t __n, size_t __m, size_t __r,
444	_UIntType __a, size_t __u, _UIntType __d, size_t __s,
445	_UIntType __b, size_t __t, _UIntType __c, size_t __l,
446	_UIntType __f>
447	typename
448	mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
449	__s, __b, __t, __c, __l, __f>::result_type
450	mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
451	__s, __b, __t, __c, __l, __f>::
452	operator()()
453	{
454	// Reload the vector - cost is O(n) amortized over n calls.
455	if (_M_p >= state_size)
456	_M_gen_rand();
457
458	// Calculate o(x(i)).
459	result_type __z = _M_x[_M_p++];
460	__z ^= (__z >> __u) & __d;
461	__z ^= (__z << __s) & __b;
462	__z ^= (__z << __t) & __c;
463	__z ^= (__z >> __l);
464
465	return __z;
466	}
467
468	template<typename _UIntType, size_t __w,
469	size_t __n, size_t __m, size_t __r,
470	_UIntType __a, size_t __u, _UIntType __d, size_t __s,
471	_UIntType __b, size_t __t, _UIntType __c, size_t __l,
472	_UIntType __f, typename _CharT, typename _Traits>
473	std::basic_ostream<_CharT, _Traits>&
474	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
475	const mersenne_twister_engine<_UIntType, __w, __n, __m,
476	__r, __a, __u, __d, __s, __b, __t, __c, __l, __f>& __x)
477	{
478	using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
479
480	const typename __ios_base::fmtflags __flags = __os.flags();
481	const _CharT __fill = __os.fill();
482	const _CharT __space = __os.widen(' ');
483	__os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
484	__os.fill(__space);
485
486	for (size_t __i = 0; __i < __n; ++__i)
487	__os << __x._M_x[__i] << __space;
488	__os << __x._M_p;
489
490	__os.flags(__flags);
491	__os.fill(__fill);
492	return __os;
493	}
494
495	template<typename _UIntType, size_t __w,
496	size_t __n, size_t __m, size_t __r,
497	_UIntType __a, size_t __u, _UIntType __d, size_t __s,
498	_UIntType __b, size_t __t, _UIntType __c, size_t __l,
499	_UIntType __f, typename _CharT, typename _Traits>
500	std::basic_istream<_CharT, _Traits>&
501	operator>>(std::basic_istream<_CharT, _Traits>& __is,
502	mersenne_twister_engine<_UIntType, __w, __n, __m,
503	__r, __a, __u, __d, __s, __b, __t, __c, __l, __f>& __x)
504	{
505	using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
506
507	const typename __ios_base::fmtflags __flags = __is.flags();
508	__is.flags(__ios_base::dec | __ios_base::skipws);
509
510	for (size_t __i = 0; __i < __n; ++__i)
511	__is >> __x._M_x[__i];
512	__is >> __x._M_p;
513
514	__is.flags(__flags);
515	return __is;
516	}
517
518
519	template<typename _UIntType, size_t __w, size_t __s, size_t __r>
520	constexpr size_t
521	subtract_with_carry_engine<_UIntType, __w, __s, __r>::word_size;
522
523	template<typename _UIntType, size_t __w, size_t __s, size_t __r>
524	constexpr size_t
525	subtract_with_carry_engine<_UIntType, __w, __s, __r>::short_lag;
526
527	template<typename _UIntType, size_t __w, size_t __s, size_t __r>
528	constexpr size_t
529	subtract_with_carry_engine<_UIntType, __w, __s, __r>::long_lag;
530
531	template<typename _UIntType, size_t __w, size_t __s, size_t __r>
532	constexpr _UIntType
533	subtract_with_carry_engine<_UIntType, __w, __s, __r>::default_seed;
534
535	template<typename _UIntType, size_t __w, size_t __s, size_t __r>
536	void
537	subtract_with_carry_engine<_UIntType, __w, __s, __r>::
538	seed(result_type __value)
539	{
540	std::linear_congruential_engine<result_type, 40014u, 0u, 2147483563u>
541	__lcg(__value == 0u ? default_seed : __value);
542
543	const size_t __n = (__w + 31) / 32;
544
545	for (size_t __i = 0; __i < long_lag; ++__i)
546	{
547	_UIntType __sum = 0u;
548	_UIntType __factor = 1u;
549	for (size_t __j = 0; __j < __n; ++__j)
550	{
551	__sum += __detail::__mod<uint_least32_t,
552	__detail::_Shift<uint_least32_t, 32>::__value>
553	(__lcg()) * __factor;
554	__factor *= __detail::_Shift<_UIntType, 32>::__value;
555	}
556	_M_x[__i] = __detail::__mod<_UIntType,
557	__detail::_Shift<_UIntType, __w>::__value>(__sum);
558	}
559	_M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
560	_M_p = 0;
561	}
562
563	template<typename _UIntType, size_t __w, size_t __s, size_t __r>
564	template<typename _Sseq>
565	auto
566	subtract_with_carry_engine<_UIntType, __w, __s, __r>::
567	seed(_Sseq& __q)
568	-> _If_seed_seq<_Sseq>
569	{
570	const size_t __k = (__w + 31) / 32;
571	uint_least32_t __arr[__r * __k];
572	__q.generate(__arr + 0, __arr + __r * __k);
573
574	for (size_t __i = 0; __i < long_lag; ++__i)
575	{
576	_UIntType __sum = 0u;
577	_UIntType __factor = 1u;
578	for (size_t __j = 0; __j < __k; ++__j)
579	{
580	__sum += __arr[__k * __i + __j] * __factor;
581	__factor *= __detail::_Shift<_UIntType, 32>::__value;
582	}
583	_M_x[__i] = __detail::__mod<_UIntType,
584	__detail::_Shift<_UIntType, __w>::__value>(__sum);
585	}
586	_M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
587	_M_p = 0;
588	}
589
590	template<typename _UIntType, size_t __w, size_t __s, size_t __r>
591	typename subtract_with_carry_engine<_UIntType, __w, __s, __r>::
592	result_type
593	subtract_with_carry_engine<_UIntType, __w, __s, __r>::
594	operator()()
595	{
596	// Derive short lag index from current index.
597	long __ps = _M_p - short_lag;
598	if (__ps < 0)
599	__ps += long_lag;
600
601	// Calculate new x(i) without overflow or division.
602	// NB: Thanks to the requirements for _UIntType, _M_x[_M_p] + _M_carry
603	// cannot overflow.
604	_UIntType __xi;
605	if (_M_x[__ps] >= _M_x[_M_p] + _M_carry)
606	{
607	__xi = _M_x[__ps] - _M_x[_M_p] - _M_carry;
608	_M_carry = 0;
609	}
610	else
611	{
612	__xi = (__detail::_Shift<_UIntType, __w>::__value
613	- _M_x[_M_p] - _M_carry + _M_x[__ps]);
614	_M_carry = 1;
615	}
616	_M_x[_M_p] = __xi;
617
618	// Adjust current index to loop around in ring buffer.
619	if (++_M_p >= long_lag)
620	_M_p = 0;
621
622	return __xi;
623	}
624
625	template<typename _UIntType, size_t __w, size_t __s, size_t __r,
626	typename _CharT, typename _Traits>
627	std::basic_ostream<_CharT, _Traits>&
628	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
629	const subtract_with_carry_engine<_UIntType,
630	__w, __s, __r>& __x)
631	{
632	using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
633
634	const typename __ios_base::fmtflags __flags = __os.flags();
635	const _CharT __fill = __os.fill();
636	const _CharT __space = __os.widen(' ');
637	__os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
638	__os.fill(__space);
639
640	for (size_t __i = 0; __i < __r; ++__i)
641	__os << __x._M_x[__i] << __space;
642	__os << __x._M_carry << __space << __x._M_p;
643
644	__os.flags(__flags);
645	__os.fill(__fill);
646	return __os;
647	}
648
649	template<typename _UIntType, size_t __w, size_t __s, size_t __r,
650	typename _CharT, typename _Traits>
651	std::basic_istream<_CharT, _Traits>&
652	operator>>(std::basic_istream<_CharT, _Traits>& __is,
653	subtract_with_carry_engine<_UIntType, __w, __s, __r>& __x)
654	{
655	using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
656
657	const typename __ios_base::fmtflags __flags = __is.flags();
658	__is.flags(__ios_base::dec | __ios_base::skipws);
659
660	for (size_t __i = 0; __i < __r; ++__i)
661	__is >> __x._M_x[__i];
662	__is >> __x._M_carry;
663	__is >> __x._M_p;
664
665	__is.flags(__flags);
666	return __is;
667	}
668
669
670	template<typename _RandomNumberEngine, size_t __p, size_t __r>
671	constexpr size_t
672	discard_block_engine<_RandomNumberEngine, __p, __r>::block_size;
673
674	template<typename _RandomNumberEngine, size_t __p, size_t __r>
675	constexpr size_t
676	discard_block_engine<_RandomNumberEngine, __p, __r>::used_block;
677
678	template<typename _RandomNumberEngine, size_t __p, size_t __r>
679	typename discard_block_engine<_RandomNumberEngine,
680	__p, __r>::result_type
681	discard_block_engine<_RandomNumberEngine, __p, __r>::
682	operator()()
683	{
684	if (_M_n >= used_block)
685	{
686	_M_b.discard(block_size - _M_n);
687	_M_n = 0;
688	}
689	++_M_n;
690	return _M_b();
691	}
692
693	template<typename _RandomNumberEngine, size_t __p, size_t __r,
694	typename _CharT, typename _Traits>
695	std::basic_ostream<_CharT, _Traits>&
696	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
697	const discard_block_engine<_RandomNumberEngine,
698	__p, __r>& __x)
699	{
700	using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
701
702	const typename __ios_base::fmtflags __flags = __os.flags();
703	const _CharT __fill = __os.fill();
704	const _CharT __space = __os.widen(' ');
705	__os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
706	__os.fill(__space);
707
708	__os << __x.base() << __space << __x._M_n;
709
710	__os.flags(__flags);
711	__os.fill(__fill);
712	return __os;
713	}
714
715	template<typename _RandomNumberEngine, size_t __p, size_t __r,
716	typename _CharT, typename _Traits>
717	std::basic_istream<_CharT, _Traits>&
718	operator>>(std::basic_istream<_CharT, _Traits>& __is,
719	discard_block_engine<_RandomNumberEngine, __p, __r>& __x)
720	{
721	using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
722
723	const typename __ios_base::fmtflags __flags = __is.flags();
724	__is.flags(__ios_base::dec | __ios_base::skipws);
725
726	__is >> __x._M_b >> __x._M_n;
727
728	__is.flags(__flags);
729	return __is;
730	}
731
732
733	template<typename _RandomNumberEngine, size_t __w, typename _UIntType>
734	typename independent_bits_engine<_RandomNumberEngine, __w, _UIntType>::
735	result_type
736	independent_bits_engine<_RandomNumberEngine, __w, _UIntType>::
737	operator()()
738	{
739	typedef typename _RandomNumberEngine::result_type _Eresult_type;
740	const _Eresult_type __r
741	= (_M_b.max() - _M_b.min() < std::numeric_limits<_Eresult_type>::max()
742	? _M_b.max() - _M_b.min() + 1 : 0);
743	const unsigned __edig = std::numeric_limits<_Eresult_type>::digits;
744	const unsigned __m = __r ? std::__lg(__r) : __edig;
745
746	typedef typename std::common_type<_Eresult_type, result_type>::type
747	__ctype;
748	const unsigned __cdig = std::numeric_limits<__ctype>::digits;
749
750	unsigned __n, __n0;
751	__ctype __s0, __s1, __y0, __y1;
752
753	for (size_t __i = 0; __i < 2; ++__i)
754	{
755	__n = (__w + __m - 1) / __m + __i;
756	__n0 = __n - __w % __n;
757	const unsigned __w0 = __w / __n; // __w0 <= __m
758
759	__s0 = 0;
760	__s1 = 0;
761	if (__w0 < __cdig)
762	{
763	__s0 = __ctype(1) << __w0;
764	__s1 = __s0 << 1;
765	}
766
767	__y0 = 0;
768	__y1 = 0;
769	if (__r)
770	{
771	__y0 = __s0 * (__r / __s0);
772	if (__s1)
773	__y1 = __s1 * (__r / __s1);
774
775	if (__r - __y0 <= __y0 / __n)
776	break;
777	}
778	else
779	break;
780	}
781
782	result_type __sum = 0;
783	for (size_t __k = 0; __k < __n0; ++__k)
784	{
785	__ctype __u;
786	do
787	__u = _M_b() - _M_b.min();
788	while (__y0 && __u >= __y0);
789	__sum = __s0 * __sum + (__s0 ? __u % __s0 : __u);
790	}
791	for (size_t __k = __n0; __k < __n; ++__k)
792	{
793	__ctype __u;
794	do
795	__u = _M_b() - _M_b.min();
796	while (__y1 && __u >= __y1);
797	__sum = __s1 * __sum + (__s1 ? __u % __s1 : __u);
798	}
799	return __sum;
800	}
801
802
803	template<typename _RandomNumberEngine, size_t __k>
804	constexpr size_t
805	shuffle_order_engine<_RandomNumberEngine, __k>::table_size;
806
807	namespace __detail
808	{
809	// Determine whether an integer is representable as double.
810	template<typename _Tp>
811	constexpr bool
812	__representable_as_double(_Tp __x) noexcept
813	{
814	static_assert(numeric_limits<_Tp>::is_integer, "");
815	static_assert(!numeric_limits<_Tp>::is_signed, "");
816	// All integers <= 2^53 are representable.
817	return (__x <= (1ull << __DBL_MANT_DIG__))
818	// Between 2^53 and 2^54 only even numbers are representable.
819	|| (!(__x & 1) && __detail::__representable_as_double(__x >> 1));
820	}
821
822	// Determine whether x+1 is representable as double.
823	template<typename _Tp>
824	constexpr bool
825	__p1_representable_as_double(_Tp __x) noexcept
826	{
827	static_assert(numeric_limits<_Tp>::is_integer, "");
828	static_assert(!numeric_limits<_Tp>::is_signed, "");
829	return numeric_limits<_Tp>::digits < __DBL_MANT_DIG__
830	|| (bool(__x + 1u) // return false if x+1 wraps around to zero
831	&& __detail::__representable_as_double(__x + 1u));
832	}
833	}
834
835	template<typename _RandomNumberEngine, size_t __k>
836	typename shuffle_order_engine<_RandomNumberEngine, __k>::result_type
837	shuffle_order_engine<_RandomNumberEngine, __k>::
838	operator()()
839	{
840	constexpr result_type __range = max() - min();
841	size_t __j = __k;
842	const result_type __y = _M_y - min();
843	// Avoid using slower long double arithmetic if possible.
844	if _GLIBCXX17_CONSTEXPR (__detail::__p1_representable_as_double(__range))
845	__j *= __y / (__range + 1.0);
846	else
847	__j *= __y / (__range + 1.0L);
848	_M_y = _M_v[__j];
849	_M_v[__j] = _M_b();
850
851	return _M_y;
852	}
853
854	template<typename _RandomNumberEngine, size_t __k,
855	typename _CharT, typename _Traits>
856	std::basic_ostream<_CharT, _Traits>&
857	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
858	const shuffle_order_engine<_RandomNumberEngine, __k>& __x)
859	{
860	using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
861
862	const typename __ios_base::fmtflags __flags = __os.flags();
863	const _CharT __fill = __os.fill();
864	const _CharT __space = __os.widen(' ');
865	__os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
866	__os.fill(__space);
867
868	__os << __x.base();
869	for (size_t __i = 0; __i < __k; ++__i)
870	__os << __space << __x._M_v[__i];
871	__os << __space << __x._M_y;
872
873	__os.flags(__flags);
874	__os.fill(__fill);
875	return __os;
876	}
877
878	template<typename _RandomNumberEngine, size_t __k,
879	typename _CharT, typename _Traits>
880	std::basic_istream<_CharT, _Traits>&
881	operator>>(std::basic_istream<_CharT, _Traits>& __is,
882	shuffle_order_engine<_RandomNumberEngine, __k>& __x)
883	{
884	using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
885
886	const typename __ios_base::fmtflags __flags = __is.flags();
887	__is.flags(__ios_base::dec | __ios_base::skipws);
888
889	__is >> __x._M_b;
890	for (size_t __i = 0; __i < __k; ++__i)
891	__is >> __x._M_v[__i];
892	__is >> __x._M_y;
893
894	__is.flags(__flags);
895	return __is;
896	}
897
898
899	template<typename _IntType, typename _CharT, typename _Traits>
900	std::basic_ostream<_CharT, _Traits>&
901	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
902	const uniform_int_distribution<_IntType>& __x)
903	{
904	using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
905
906	const typename __ios_base::fmtflags __flags = __os.flags();
907	const _CharT __fill = __os.fill();
908	const _CharT __space = __os.widen(' ');
909	__os.flags(__ios_base::scientific | __ios_base::left);
910	__os.fill(__space);
911
912	__os << __x.a() << __space << __x.b();
913
914	__os.flags(__flags);
915	__os.fill(__fill);
916	return __os;
917	}
918
919	template<typename _IntType, typename _CharT, typename _Traits>
920	std::basic_istream<_CharT, _Traits>&
921	operator>>(std::basic_istream<_CharT, _Traits>& __is,
922	uniform_int_distribution<_IntType>& __x)
923	{
924	using param_type
925	= typename uniform_int_distribution<_IntType>::param_type;
926	using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
927
928	const typename __ios_base::fmtflags __flags = __is.flags();
929	__is.flags(__ios_base::dec | __ios_base::skipws);
930
931	_IntType __a, __b;
932	if (__is >> __a >> __b)
933	__x.param(param_type(__a, __b));
934
935	__is.flags(__flags);
936	return __is;
937	}
938
939
940	template<typename _RealType>
941	template<typename _ForwardIterator,
942	typename _UniformRandomNumberGenerator>
943	void
944	uniform_real_distribution<_RealType>::
945	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
946	_UniformRandomNumberGenerator& __urng,
947	const param_type& __p)
948	{
949	__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
950	__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
951	__aurng(__urng);
952	auto __range = __p.b() - __p.a();
953	while (__f != __t)
954	*__f++ = __aurng() * __range + __p.a();
955	}
956
957	template<typename _RealType, typename _CharT, typename _Traits>
958	std::basic_ostream<_CharT, _Traits>&
959	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
960	const uniform_real_distribution<_RealType>& __x)
961	{
962	using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
963
964	const typename __ios_base::fmtflags __flags = __os.flags();
965	const _CharT __fill = __os.fill();
966	const std::streamsize __precision = __os.precision();
967	const _CharT __space = __os.widen(' ');
968	__os.flags(__ios_base::scientific | __ios_base::left);
969	__os.fill(__space);
970	__os.precision(std::numeric_limits<_RealType>::max_digits10);
971
972	__os << __x.a() << __space << __x.b();
973
974	__os.flags(__flags);
975	__os.fill(__fill);
976	__os.precision(__precision);
977	return __os;
978	}
979
980	template<typename _RealType, typename _CharT, typename _Traits>
981	std::basic_istream<_CharT, _Traits>&
982	operator>>(std::basic_istream<_CharT, _Traits>& __is,
983	uniform_real_distribution<_RealType>& __x)
984	{
985	using param_type
986	= typename uniform_real_distribution<_RealType>::param_type;
987	using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
988
989	const typename __ios_base::fmtflags __flags = __is.flags();
990	__is.flags(__ios_base::skipws);
991
992	_RealType __a, __b;
993	if (__is >> __a >> __b)
994	__x.param(param_type(__a, __b));
995
996	__is.flags(__flags);
997	return __is;
998	}
999
1000
1001	template<typename _ForwardIterator,
1002	typename _UniformRandomNumberGenerator>
1003	void
1004	std::bernoulli_distribution::
1005	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1006	_UniformRandomNumberGenerator& __urng,
1007	const param_type& __p)
1008	{
1009	__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1010	__detail::_Adaptor<_UniformRandomNumberGenerator, double>
1011	__aurng(__urng);
1012	auto __limit = __p.p() * (__aurng.max() - __aurng.min());
1013
1014	while (__f != __t)
1015	*__f++ = (__aurng() - __aurng.min()) < __limit;
1016	}
1017
1018	template<typename _CharT, typename _Traits>
1019	std::basic_ostream<_CharT, _Traits>&
1020	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1021	const bernoulli_distribution& __x)
1022	{
1023	using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1024
1025	const typename __ios_base::fmtflags __flags = __os.flags();
1026	const _CharT __fill = __os.fill();
1027	const std::streamsize __precision = __os.precision();
1028	__os.flags(__ios_base::scientific | __ios_base::left);
1029	__os.fill(__os.widen(' '));
1030	__os.precision(std::numeric_limits<double>::max_digits10);
1031
1032	__os << __x.p();
1033
1034	__os.flags(__flags);
1035	__os.fill(__fill);
1036	__os.precision(__precision);
1037	return __os;
1038	}
1039
1040
1041	template<typename _IntType>
1042	template<typename _UniformRandomNumberGenerator>
1043	typename geometric_distribution<_IntType>::result_type
1044	geometric_distribution<_IntType>::
1045	operator()(_UniformRandomNumberGenerator& __urng,
1046	const param_type& __param)
1047	{
1048	// About the epsilon thing see this thread:
1049	// http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html
1050	const double __naf =
1051	(1 - std::numeric_limits<double>::epsilon()) / 2;
1052	// The largest _RealType convertible to _IntType.
1053	const double __thr =
1054	std::numeric_limits<_IntType>::max() + __naf;
1055	__detail::_Adaptor<_UniformRandomNumberGenerator, double>
1056	__aurng(__urng);
1057
1058	double __cand;
1059	do
1060	__cand = std::floor(std::log(1.0 - __aurng()) / __param._M_log_1_p);
1061	while (__cand >= __thr);
1062
1063	return result_type(__cand + __naf);
1064	}
1065
1066	template<typename _IntType>
1067	template<typename _ForwardIterator,
1068	typename _UniformRandomNumberGenerator>
1069	void
1070	geometric_distribution<_IntType>::
1071	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1072	_UniformRandomNumberGenerator& __urng,
1073	const param_type& __param)
1074	{
1075	__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1076	// About the epsilon thing see this thread:
1077	// http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html
1078	const double __naf =
1079	(1 - std::numeric_limits<double>::epsilon()) / 2;
1080	// The largest _RealType convertible to _IntType.
1081	const double __thr =
1082	std::numeric_limits<_IntType>::max() + __naf;
1083	__detail::_Adaptor<_UniformRandomNumberGenerator, double>
1084	__aurng(__urng);
1085
1086	while (__f != __t)
1087	{
1088	double __cand;
1089	do
1090	__cand = std::floor(std::log(1.0 - __aurng())
1091	/ __param._M_log_1_p);
1092	while (__cand >= __thr);
1093
1094	*__f++ = __cand + __naf;
1095	}
1096	}
1097
1098	template<typename _IntType,
1099	typename _CharT, typename _Traits>
1100	std::basic_ostream<_CharT, _Traits>&
1101	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1102	const geometric_distribution<_IntType>& __x)
1103	{
1104	using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1105
1106	const typename __ios_base::fmtflags __flags = __os.flags();
1107	const _CharT __fill = __os.fill();
1108	const std::streamsize __precision = __os.precision();
1109	__os.flags(__ios_base::scientific | __ios_base::left);
1110	__os.fill(__os.widen(' '));
1111	__os.precision(std::numeric_limits<double>::max_digits10);
1112
1113	__os << __x.p();
1114
1115	__os.flags(__flags);
1116	__os.fill(__fill);
1117	__os.precision(__precision);
1118	return __os;
1119	}
1120
1121	template<typename _IntType,
1122	typename _CharT, typename _Traits>
1123	std::basic_istream<_CharT, _Traits>&
1124	operator>>(std::basic_istream<_CharT, _Traits>& __is,
1125	geometric_distribution<_IntType>& __x)
1126	{
1127	using param_type = typename geometric_distribution<_IntType>::param_type;
1128	using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
1129
1130	const typename __ios_base::fmtflags __flags = __is.flags();
1131	__is.flags(__ios_base::skipws);
1132
1133	double __p;
1134	if (__is >> __p)
1135	__x.param(param_type(__p));
1136
1137	__is.flags(__flags);
1138	return __is;
1139	}
1140
1141	// This is Leger's algorithm, also in Devroye, Ch. X, Example 1.5.
1142	template<typename _IntType>
1143	template<typename _UniformRandomNumberGenerator>
1144	typename negative_binomial_distribution<_IntType>::result_type
1145	negative_binomial_distribution<_IntType>::
1146	operator()(_UniformRandomNumberGenerator& __urng)
1147	{
1148	const double __y = _M_gd(__urng);
1149
1150	// XXX Is the constructor too slow?
1151	std::poisson_distribution<result_type> __poisson(__y);
1152	return __poisson(__urng);
1153	}
1154
1155	template<typename _IntType>
1156	template<typename _UniformRandomNumberGenerator>
1157	typename negative_binomial_distribution<_IntType>::result_type
1158	negative_binomial_distribution<_IntType>::
1159	operator()(_UniformRandomNumberGenerator& __urng,
1160	const param_type& __p)
1161	{
1162	typedef typename std::gamma_distribution<double>::param_type
1163	param_type;
1164
1165	const double __y =
1166	_M_gd(__urng, param_type(__p.k(), (1.0 - __p.p()) / __p.p()));
1167
1168	std::poisson_distribution<result_type> __poisson(__y);
1169	return __poisson(__urng);
1170	}
1171
1172	template<typename _IntType>
1173	template<typename _ForwardIterator,
1174	typename _UniformRandomNumberGenerator>
1175	void
1176	negative_binomial_distribution<_IntType>::
1177	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1178	_UniformRandomNumberGenerator& __urng)
1179	{
1180	__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1181	while (__f != __t)
1182	{
1183	const double __y = _M_gd(__urng);
1184
1185	// XXX Is the constructor too slow?
1186	std::poisson_distribution<result_type> __poisson(__y);
1187	*__f++ = __poisson(__urng);
1188	}
1189	}
1190
1191	template<typename _IntType>
1192	template<typename _ForwardIterator,
1193	typename _UniformRandomNumberGenerator>
1194	void
1195	negative_binomial_distribution<_IntType>::
1196	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1197	_UniformRandomNumberGenerator& __urng,
1198	const param_type& __p)
1199	{
1200	__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1201	typename std::gamma_distribution<result_type>::param_type
1202	__p2(__p.k(), (1.0 - __p.p()) / __p.p());
1203
1204	while (__f != __t)
1205	{
1206	const double __y = _M_gd(__urng, __p2);
1207
1208	std::poisson_distribution<result_type> __poisson(__y);
1209	*__f++ = __poisson(__urng);
1210	}
1211	}
1212
1213	template<typename _IntType, typename _CharT, typename _Traits>
1214	std::basic_ostream<_CharT, _Traits>&
1215	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1216	const negative_binomial_distribution<_IntType>& __x)
1217	{
1218	using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1219
1220	const typename __ios_base::fmtflags __flags = __os.flags();
1221	const _CharT __fill = __os.fill();
1222	const std::streamsize __precision = __os.precision();
1223	const _CharT __space = __os.widen(' ');
1224	__os.flags(__ios_base::scientific | __ios_base::left);
1225	__os.fill(__os.widen(' '));
1226	__os.precision(std::numeric_limits<double>::max_digits10);
1227
1228	__os << __x.k() << __space << __x.p()
1229	<< __space << __x._M_gd;
1230
1231	__os.flags(__flags);
1232	__os.fill(__fill);
1233	__os.precision(__precision);
1234	return __os;
1235	}
1236
1237	template<typename _IntType, typename _CharT, typename _Traits>
1238	std::basic_istream<_CharT, _Traits>&
1239	operator>>(std::basic_istream<_CharT, _Traits>& __is,
1240	negative_binomial_distribution<_IntType>& __x)
1241	{
1242	using param_type
1243	= typename negative_binomial_distribution<_IntType>::param_type;
1244	using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
1245
1246	const typename __ios_base::fmtflags __flags = __is.flags();
1247	__is.flags(__ios_base::skipws);
1248
1249	_IntType __k;
1250	double __p;
1251	if (__is >> __k >> __p >> __x._M_gd)
1252	__x.param(param_type(__k, __p));
1253
1254	__is.flags(__flags);
1255	return __is;
1256	}
1257
1258
1259	template<typename _IntType>
1260	void
1261	poisson_distribution<_IntType>::param_type::
1262	_M_initialize()
1263	{
1264	#if _GLIBCXX_USE_C99_MATH_TR1
1265	if (_M_mean >= 12)
1266	{
1267	const double __m = std::floor(x: _M_mean);
1268	_M_lm_thr = std::log(x: _M_mean);
1269	_M_lfm = std::lgamma(__m + 1);
1270	_M_sm = std::sqrt(x: __m);
1271
1272	const double __pi_4 = 0.7853981633974483096156608458198757L;
1273	const double __dx = std::sqrt(x: 2 * __m * std::log(x: 32 * __m
1274	/ __pi_4));
1275	_M_d = std::round(x: std::max<double>(a: 6.0, b: std::min(a: __m, b: __dx)));
1276	const double __cx = 2 * __m + _M_d;
1277	_M_scx = std::sqrt(x: __cx / 2);
1278	_M_1cx = 1 / __cx;
1279
1280	_M_c2b = std::sqrt(x: __pi_4 * __cx) * std::exp(x: _M_1cx);
1281	_M_cb = 2 * __cx * std::exp(x: -_M_d * _M_1cx * (1 + _M_d / 2))
1282	/ _M_d;
1283	}
1284	else
1285	#endif
1286	_M_lm_thr = std::exp(x: -_M_mean);
1287	}
1288
1289	/**
1290	* A rejection algorithm when mean >= 12 and a simple method based
1291	* upon the multiplication of uniform random variates otherwise.
1292	* NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
1293	* is defined.
1294	*
1295	* Reference:
1296	* Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
1297	* New York, 1986, Ch. X, Sects. 3.3 & 3.4 (+ Errata!).
1298	*/
1299	template<typename _IntType>
1300	template<typename _UniformRandomNumberGenerator>
1301	typename poisson_distribution<_IntType>::result_type
1302	poisson_distribution<_IntType>::
1303	operator()(_UniformRandomNumberGenerator& __urng,
1304	const param_type& __param)
1305	{
1306	__detail::_Adaptor<_UniformRandomNumberGenerator, double>
1307	__aurng(__urng);
1308	#if _GLIBCXX_USE_C99_MATH_TR1
1309	if (__param.mean() >= 12)
1310	{
1311	double __x;
1312
1313	// See comments above...
1314	const double __naf =
1315	(1 - std::numeric_limits<double>::epsilon()) / 2;
1316	const double __thr =
1317	std::numeric_limits<_IntType>::max() + __naf;
1318
1319	const double __m = std::floor(__param.mean());
1320	// sqrt(pi / 2)
1321	const double __spi_2 = 1.2533141373155002512078826424055226L;
1322	const double __c1 = __param._M_sm * __spi_2;
1323	const double __c2 = __param._M_c2b + __c1;
1324	const double __c3 = __c2 + 1;
1325	const double __c4 = __c3 + 1;
1326	// 1 / 78
1327	const double __178 = 0.0128205128205128205128205128205128L;
1328	// e^(1 / 78)
1329	const double __e178 = 1.0129030479320018583185514777512983L;
1330	const double __c5 = __c4 + __e178;
1331	const double __c = __param._M_cb + __c5;
1332	const double __2cx = 2 * (2 * __m + __param._M_d);
1333
1334	bool __reject = true;
1335	do
1336	{
1337	const double __u = __c * __aurng();
1338	const double __e = -std::log(1.0 - __aurng());
1339
1340	double __w = 0.0;
1341
1342	if (__u <= __c1)
1343	{
1344	const double __n = _M_nd(__urng);
1345	const double __y = -std::abs(x: __n) * __param._M_sm - 1;
1346	__x = std::floor(x: __y);
1347	__w = -__n * __n / 2;
1348	if (__x < -__m)
1349	continue;
1350	}
1351	else if (__u <= __c2)
1352	{
1353	const double __n = _M_nd(__urng);
1354	const double __y = 1 + std::abs(x: __n) * __param._M_scx;
1355	__x = std::ceil(x: __y);
1356	__w = __y * (2 - __y) * __param._M_1cx;
1357	if (__x > __param._M_d)
1358	continue;
1359	}
1360	else if (__u <= __c3)
1361	// NB: This case not in the book, nor in the Errata,
1362	// but should be ok...
1363	__x = -1;
1364	else if (__u <= __c4)
1365	__x = 0;
1366	else if (__u <= __c5)
1367	{
1368	__x = 1;
1369	// Only in the Errata, see libstdc++/83237.
1370	__w = __178;
1371	}
1372	else
1373	{
1374	const double __v = -std::log(1.0 - __aurng());
1375	const double __y = __param._M_d
1376	+ __v * __2cx / __param._M_d;
1377	__x = std::ceil(x: __y);
1378	__w = -__param._M_d * __param._M_1cx * (1 + __y / 2);
1379	}
1380
1381	__reject = (__w - __e - __x * __param._M_lm_thr
1382	> __param._M_lfm - std::lgamma(__x + __m + 1));
1383
1384	__reject |= __x + __m >= __thr;
1385
1386	} while (__reject);
1387
1388	return result_type(__x + __m + __naf);
1389	}
1390	else
1391	#endif
1392	{
1393	_IntType __x = 0;
1394	double __prod = 1.0;
1395
1396	do
1397	{
1398	__prod *= __aurng();
1399	__x += 1;
1400	}
1401	while (__prod > __param._M_lm_thr);
1402
1403	return __x - 1;
1404	}
1405	}
1406
1407	template<typename _IntType>
1408	template<typename _ForwardIterator,
1409	typename _UniformRandomNumberGenerator>
1410	void
1411	poisson_distribution<_IntType>::
1412	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1413	_UniformRandomNumberGenerator& __urng,
1414	const param_type& __param)
1415	{
1416	__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1417	// We could duplicate everything from operator()...
1418	while (__f != __t)
1419	*__f++ = this->operator()(__urng, __param);
1420	}
1421
1422	template<typename _IntType,
1423	typename _CharT, typename _Traits>
1424	std::basic_ostream<_CharT, _Traits>&
1425	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1426	const poisson_distribution<_IntType>& __x)
1427	{
1428	using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1429
1430	const typename __ios_base::fmtflags __flags = __os.flags();
1431	const _CharT __fill = __os.fill();
1432	const std::streamsize __precision = __os.precision();
1433	const _CharT __space = __os.widen(' ');
1434	__os.flags(__ios_base::scientific | __ios_base::left);
1435	__os.fill(__space);
1436	__os.precision(std::numeric_limits<double>::max_digits10);
1437
1438	__os << __x.mean() << __space << __x._M_nd;
1439
1440	__os.flags(__flags);
1441	__os.fill(__fill);
1442	__os.precision(__precision);
1443	return __os;
1444	}
1445
1446	template<typename _IntType,
1447	typename _CharT, typename _Traits>
1448	std::basic_istream<_CharT, _Traits>&
1449	operator>>(std::basic_istream<_CharT, _Traits>& __is,
1450	poisson_distribution<_IntType>& __x)
1451	{
1452	using param_type = typename poisson_distribution<_IntType>::param_type;
1453	using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
1454
1455	const typename __ios_base::fmtflags __flags = __is.flags();
1456	__is.flags(__ios_base::skipws);
1457
1458	double __mean;
1459	if (__is >> __mean >> __x._M_nd)
1460	__x.param(param_type(__mean));
1461
1462	__is.flags(__flags);
1463	return __is;
1464	}
1465
1466
1467	template<typename _IntType>
1468	void
1469	binomial_distribution<_IntType>::param_type::
1470	_M_initialize()
1471	{
1472	const double __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p;
1473
1474	_M_easy = true;
1475
1476	#if _GLIBCXX_USE_C99_MATH_TR1
1477	if (_M_t * __p12 >= 8)
1478	{
1479	_M_easy = false;
1480	const double __np = std::floor(_M_t * __p12);
1481	const double __pa = __np / _M_t;
1482	const double __1p = 1 - __pa;
1483
1484	const double __pi_4 = 0.7853981633974483096156608458198757L;
1485	const double __d1x =
1486	std::sqrt(x: __np * __1p * std::log(x: 32 * __np
1487	/ (81 * __pi_4 * __1p)));
1488	_M_d1 = std::round(x: std::max<double>(a: 1.0, b: __d1x));
1489	const double __d2x =
1490	std::sqrt(__np * __1p * std::log(32 * _M_t * __1p
1491	/ (__pi_4 * __pa)));
1492	_M_d2 = std::round(x: std::max<double>(a: 1.0, b: __d2x));
1493
1494	// sqrt(pi / 2)
1495	const double __spi_2 = 1.2533141373155002512078826424055226L;
1496	_M_s1 = std::sqrt(x: __np * __1p) * (1 + _M_d1 / (4 * __np));
1497	_M_s2 = std::sqrt(x: __np * __1p) * (1 + _M_d2 / (4 * _M_t * __1p));
1498	_M_c = 2 * _M_d1 / __np;
1499	_M_a1 = std::exp(x: _M_c) * _M_s1 * __spi_2;
1500	const double __a12 = _M_a1 + _M_s2 * __spi_2;
1501	const double __s1s = _M_s1 * _M_s1;
1502	_M_a123 = __a12 + (std::exp(_M_d1 / (_M_t * __1p))
1503	* 2 * __s1s / _M_d1
1504	* std::exp(x: -_M_d1 * _M_d1 / (2 * __s1s)));
1505	const double __s2s = _M_s2 * _M_s2;
1506	_M_s = (_M_a123 + 2 * __s2s / _M_d2
1507	* std::exp(x: -_M_d2 * _M_d2 / (2 * __s2s)));
1508	_M_lf = (std::lgamma(__np + 1)
1509	+ std::lgamma(_M_t - __np + 1));
1510	_M_lp1p = std::log(x: __pa / __1p);
1511
1512	_M_q = -std::log(x: 1 - (__p12 - __pa) / __1p);
1513	}
1514	else
1515	#endif
1516	_M_q = -std::log(x: 1 - __p12);
1517	}
1518
1519	template<typename _IntType>
1520	template<typename _UniformRandomNumberGenerator>
1521	typename binomial_distribution<_IntType>::result_type
1522	binomial_distribution<_IntType>::
1523	_M_waiting(_UniformRandomNumberGenerator& __urng,
1524	_IntType __t, double __q)
1525	{
1526	_IntType __x = 0;
1527	double __sum = 0.0;
1528	__detail::_Adaptor<_UniformRandomNumberGenerator, double>
1529	__aurng(__urng);
1530
1531	do
1532	{
1533	if (__t == __x)
1534	return __x;
1535	const double __e = -std::log(1.0 - __aurng());
1536	__sum += __e / (__t - __x);
1537	__x += 1;
1538	}
1539	while (__sum <= __q);
1540
1541	return __x - 1;
1542	}
1543
1544	/**
1545	* A rejection algorithm when t * p >= 8 and a simple waiting time
1546	* method - the second in the referenced book - otherwise.
1547	* NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
1548	* is defined.
1549	*
1550	* Reference:
1551	* Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
1552	* New York, 1986, Ch. X, Sect. 4 (+ Errata!).
1553	*/
1554	template<typename _IntType>
1555	template<typename _UniformRandomNumberGenerator>
1556	typename binomial_distribution<_IntType>::result_type
1557	binomial_distribution<_IntType>::
1558	operator()(_UniformRandomNumberGenerator& __urng,
1559	const param_type& __param)
1560	{
1561	result_type __ret;
1562	const _IntType __t = __param.t();
1563	const double __p = __param.p();
1564	const double __p12 = __p <= 0.5 ? __p : 1.0 - __p;
1565	__detail::_Adaptor<_UniformRandomNumberGenerator, double>
1566	__aurng(__urng);
1567
1568	#if _GLIBCXX_USE_C99_MATH_TR1
1569	if (!__param._M_easy)
1570	{
1571	double __x;
1572
1573	// See comments above...
1574	const double __naf =
1575	(1 - std::numeric_limits<double>::epsilon()) / 2;
1576	const double __thr =
1577	std::numeric_limits<_IntType>::max() + __naf;
1578
1579	const double __np = std::floor(__t * __p12);
1580
1581	// sqrt(pi / 2)
1582	const double __spi_2 = 1.2533141373155002512078826424055226L;
1583	const double __a1 = __param._M_a1;
1584	const double __a12 = __a1 + __param._M_s2 * __spi_2;
1585	const double __a123 = __param._M_a123;
1586	const double __s1s = __param._M_s1 * __param._M_s1;
1587	const double __s2s = __param._M_s2 * __param._M_s2;
1588
1589	bool __reject;
1590	do
1591	{
1592	const double __u = __param._M_s * __aurng();
1593
1594	double __v;
1595
1596	if (__u <= __a1)
1597	{
1598	const double __n = _M_nd(__urng);
1599	const double __y = __param._M_s1 * std::abs(x: __n);
1600	__reject = __y >= __param._M_d1;
1601	if (!__reject)
1602	{
1603	const double __e = -std::log(1.0 - __aurng());
1604	__x = std::floor(x: __y);
1605	__v = -__e - __n * __n / 2 + __param._M_c;
1606	}
1607	}
1608	else if (__u <= __a12)
1609	{
1610	const double __n = _M_nd(__urng);
1611	const double __y = __param._M_s2 * std::abs(x: __n);
1612	__reject = __y >= __param._M_d2;
1613	if (!__reject)
1614	{
1615	const double __e = -std::log(1.0 - __aurng());
1616	__x = std::floor(x: -__y);
1617	__v = -__e - __n * __n / 2;
1618	}
1619	}
1620	else if (__u <= __a123)
1621	{
1622	const double __e1 = -std::log(1.0 - __aurng());
1623	const double __e2 = -std::log(1.0 - __aurng());
1624
1625	const double __y = __param._M_d1
1626	+ 2 * __s1s * __e1 / __param._M_d1;
1627	__x = std::floor(x: __y);
1628	__v = (-__e2 + __param._M_d1 * (1 / (__t - __np)
1629	-__y / (2 * __s1s)));
1630	__reject = false;
1631	}
1632	else
1633	{
1634	const double __e1 = -std::log(1.0 - __aurng());
1635	const double __e2 = -std::log(1.0 - __aurng());
1636
1637	const double __y = __param._M_d2
1638	+ 2 * __s2s * __e1 / __param._M_d2;
1639	__x = std::floor(x: -__y);
1640	__v = -__e2 - __param._M_d2 * __y / (2 * __s2s);
1641	__reject = false;
1642	}
1643
1644	__reject = __reject || __x < -__np || __x > __t - __np;
1645	if (!__reject)
1646	{
1647	const double __lfx =
1648	std::lgamma(__np + __x + 1)
1649	+ std::lgamma(__t - (__np + __x) + 1);
1650	__reject = __v > __param._M_lf - __lfx
1651	+ __x * __param._M_lp1p;
1652	}
1653
1654	__reject |= __x + __np >= __thr;
1655	}
1656	while (__reject);
1657
1658	__x += __np + __naf;
1659
1660	const _IntType __z = _M_waiting(__urng, __t - _IntType(__x),
1661	__param._M_q);
1662	__ret = _IntType(__x) + __z;
1663	}
1664	else
1665	#endif
1666	__ret = _M_waiting(__urng, __t, __param._M_q);
1667
1668	if (__p12 != __p)
1669	__ret = __t - __ret;
1670	return __ret;
1671	}
1672
1673	template<typename _IntType>
1674	template<typename _ForwardIterator,
1675	typename _UniformRandomNumberGenerator>
1676	void
1677	binomial_distribution<_IntType>::
1678	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1679	_UniformRandomNumberGenerator& __urng,
1680	const param_type& __param)
1681	{
1682	__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1683	// We could duplicate everything from operator()...
1684	while (__f != __t)
1685	*__f++ = this->operator()(__urng, __param);
1686	}
1687
1688	template<typename _IntType,
1689	typename _CharT, typename _Traits>
1690	std::basic_ostream<_CharT, _Traits>&
1691	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1692	const binomial_distribution<_IntType>& __x)
1693	{
1694	using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1695
1696	const typename __ios_base::fmtflags __flags = __os.flags();
1697	const _CharT __fill = __os.fill();
1698	const std::streamsize __precision = __os.precision();
1699	const _CharT __space = __os.widen(' ');
1700	__os.flags(__ios_base::scientific | __ios_base::left);
1701	__os.fill(__space);
1702	__os.precision(std::numeric_limits<double>::max_digits10);
1703
1704	__os << __x.t() << __space << __x.p()
1705	<< __space << __x._M_nd;
1706
1707	__os.flags(__flags);
1708	__os.fill(__fill);
1709	__os.precision(__precision);
1710	return __os;
1711	}
1712
1713	template<typename _IntType,
1714	typename _CharT, typename _Traits>
1715	std::basic_istream<_CharT, _Traits>&
1716	operator>>(std::basic_istream<_CharT, _Traits>& __is,
1717	binomial_distribution<_IntType>& __x)
1718	{
1719	using param_type = typename binomial_distribution<_IntType>::param_type;
1720	using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
1721
1722	const typename __ios_base::fmtflags __flags = __is.flags();
1723	__is.flags(__ios_base::dec | __ios_base::skipws);
1724
1725	_IntType __t;
1726	double __p;
1727	if (__is >> __t >> __p >> __x._M_nd)
1728	__x.param(param_type(__t, __p));
1729
1730	__is.flags(__flags);
1731	return __is;
1732	}
1733
1734
1735	template<typename _RealType>
1736	template<typename _ForwardIterator,
1737	typename _UniformRandomNumberGenerator>
1738	void
1739	std::exponential_distribution<_RealType>::
1740	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1741	_UniformRandomNumberGenerator& __urng,
1742	const param_type& __p)
1743	{
1744	__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1745	__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1746	__aurng(__urng);
1747	while (__f != __t)
1748	*__f++ = -std::log(result_type(1) - __aurng()) / __p.lambda();
1749	}
1750
1751	template<typename _RealType, typename _CharT, typename _Traits>
1752	std::basic_ostream<_CharT, _Traits>&
1753	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1754	const exponential_distribution<_RealType>& __x)
1755	{
1756	using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1757
1758	const typename __ios_base::fmtflags __flags = __os.flags();
1759	const _CharT __fill = __os.fill();
1760	const std::streamsize __precision = __os.precision();
1761	__os.flags(__ios_base::scientific | __ios_base::left);
1762	__os.fill(__os.widen(' '));
1763	__os.precision(std::numeric_limits<_RealType>::max_digits10);
1764
1765	__os << __x.lambda();
1766
1767	__os.flags(__flags);
1768	__os.fill(__fill);
1769	__os.precision(__precision);
1770	return __os;
1771	}
1772
1773	template<typename _RealType, typename _CharT, typename _Traits>
1774	std::basic_istream<_CharT, _Traits>&
1775	operator>>(std::basic_istream<_CharT, _Traits>& __is,
1776	exponential_distribution<_RealType>& __x)
1777	{
1778	using param_type
1779	= typename exponential_distribution<_RealType>::param_type;
1780	using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
1781
1782	const typename __ios_base::fmtflags __flags = __is.flags();
1783	__is.flags(__ios_base::dec | __ios_base::skipws);
1784
1785	_RealType __lambda;
1786	if (__is >> __lambda)
1787	__x.param(param_type(__lambda));
1788
1789	__is.flags(__flags);
1790	return __is;
1791	}
1792
1793
1794	/**
1795	* Polar method due to Marsaglia.
1796	*
1797	* Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
1798	* New York, 1986, Ch. V, Sect. 4.4.
1799	*/
1800	template<typename _RealType>
1801	template<typename _UniformRandomNumberGenerator>
1802	typename normal_distribution<_RealType>::result_type
1803	normal_distribution<_RealType>::
1804	operator()(_UniformRandomNumberGenerator& __urng,
1805	const param_type& __param)
1806	{
1807	result_type __ret;
1808	__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1809	__aurng(__urng);
1810
1811	if (_M_saved_available)
1812	{
1813	_M_saved_available = false;
1814	__ret = _M_saved;
1815	}
1816	else
1817	{
1818	result_type __x, __y, __r2;
1819	do
1820	{
1821	__x = result_type(2.0) * __aurng() - 1.0;
1822	__y = result_type(2.0) * __aurng() - 1.0;
1823	__r2 = __x * __x + __y * __y;
1824	}
1825	while (__r2 > 1.0 || __r2 == 0.0);
1826
1827	const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
1828	_M_saved = __x * __mult;
1829	_M_saved_available = true;
1830	__ret = __y * __mult;
1831	}
1832
1833	__ret = __ret * __param.stddev() + __param.mean();
1834	return __ret;
1835	}
1836
1837	template<typename _RealType>
1838	template<typename _ForwardIterator,
1839	typename _UniformRandomNumberGenerator>
1840	void
1841	normal_distribution<_RealType>::
1842	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1843	_UniformRandomNumberGenerator& __urng,
1844	const param_type& __param)
1845	{
1846	__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1847
1848	if (__f == __t)
1849	return;
1850
1851	if (_M_saved_available)
1852	{
1853	_M_saved_available = false;
1854	*__f++ = _M_saved * __param.stddev() + __param.mean();
1855
1856	if (__f == __t)
1857	return;
1858	}
1859
1860	__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1861	__aurng(__urng);
1862
1863	while (__f + 1 < __t)
1864	{
1865	result_type __x, __y, __r2;
1866	do
1867	{
1868	__x = result_type(2.0) * __aurng() - 1.0;
1869	__y = result_type(2.0) * __aurng() - 1.0;
1870	__r2 = __x * __x + __y * __y;
1871	}
1872	while (__r2 > 1.0 || __r2 == 0.0);
1873
1874	const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
1875	*__f++ = __y * __mult * __param.stddev() + __param.mean();
1876	*__f++ = __x * __mult * __param.stddev() + __param.mean();
1877	}
1878
1879	if (__f != __t)
1880	{
1881	result_type __x, __y, __r2;
1882	do
1883	{
1884	__x = result_type(2.0) * __aurng() - 1.0;
1885	__y = result_type(2.0) * __aurng() - 1.0;
1886	__r2 = __x * __x + __y * __y;
1887	}
1888	while (__r2 > 1.0 || __r2 == 0.0);
1889
1890	const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
1891	_M_saved = __x * __mult;
1892	_M_saved_available = true;
1893	*__f = __y * __mult * __param.stddev() + __param.mean();
1894	}
1895	}
1896
1897	template<typename _RealType>
1898	bool
1899	operator==(const std::normal_distribution<_RealType>& __d1,
1900	const std::normal_distribution<_RealType>& __d2)
1901	{
1902	if (__d1._M_param == __d2._M_param
1903	&& __d1._M_saved_available == __d2._M_saved_available)
1904	{
1905	if (__d1._M_saved_available
1906	&& __d1._M_saved == __d2._M_saved)
1907	return true;
1908	else if(!__d1._M_saved_available)
1909	return true;
1910	else
1911	return false;
1912	}
1913	else
1914	return false;
1915	}
1916
1917	template<typename _RealType, typename _CharT, typename _Traits>
1918	std::basic_ostream<_CharT, _Traits>&
1919	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1920	const normal_distribution<_RealType>& __x)
1921	{
1922	using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1923
1924	const typename __ios_base::fmtflags __flags = __os.flags();
1925	const _CharT __fill = __os.fill();
1926	const std::streamsize __precision = __os.precision();
1927	const _CharT __space = __os.widen(' ');
1928	__os.flags(__ios_base::scientific | __ios_base::left);
1929	__os.fill(__space);
1930	__os.precision(std::numeric_limits<_RealType>::max_digits10);
1931
1932	__os << __x.mean() << __space << __x.stddev()
1933	<< __space << __x._M_saved_available;
1934	if (__x._M_saved_available)
1935	__os << __space << __x._M_saved;
1936
1937	__os.flags(__flags);
1938	__os.fill(__fill);
1939	__os.precision(__precision);
1940	return __os;
1941	}
1942
1943	template<typename _RealType, typename _CharT, typename _Traits>
1944	std::basic_istream<_CharT, _Traits>&
1945	operator>>(std::basic_istream<_CharT, _Traits>& __is,
1946	normal_distribution<_RealType>& __x)
1947	{
1948	using param_type = typename normal_distribution<_RealType>::param_type;
1949	using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
1950
1951	const typename __ios_base::fmtflags __flags = __is.flags();
1952	__is.flags(__ios_base::dec | __ios_base::skipws);
1953
1954	double __mean, __stddev;
1955	bool __saved_avail;
1956	if (__is >> __mean >> __stddev >> __saved_avail)
1957	{
1958	if (!__saved_avail || (__is >> __x._M_saved))
1959	{
1960	__x._M_saved_available = __saved_avail;
1961	__x.param(param_type(__mean, __stddev));
1962	}
1963	}
1964
1965	__is.flags(__flags);
1966	return __is;
1967	}
1968
1969
1970	template<typename _RealType>
1971	template<typename _ForwardIterator,
1972	typename _UniformRandomNumberGenerator>
1973	void
1974	lognormal_distribution<_RealType>::
1975	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1976	_UniformRandomNumberGenerator& __urng,
1977	const param_type& __p)
1978	{
1979	__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
1980	while (__f != __t)
1981	*__f++ = std::exp(__p.s() * _M_nd(__urng) + __p.m());
1982	}
1983
1984	template<typename _RealType, typename _CharT, typename _Traits>
1985	std::basic_ostream<_CharT, _Traits>&
1986	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1987	const lognormal_distribution<_RealType>& __x)
1988	{
1989	using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
1990
1991	const typename __ios_base::fmtflags __flags = __os.flags();
1992	const _CharT __fill = __os.fill();
1993	const std::streamsize __precision = __os.precision();
1994	const _CharT __space = __os.widen(' ');
1995	__os.flags(__ios_base::scientific | __ios_base::left);
1996	__os.fill(__space);
1997	__os.precision(std::numeric_limits<_RealType>::max_digits10);
1998
1999	__os << __x.m() << __space << __x.s()
2000	<< __space << __x._M_nd;
2001
2002	__os.flags(__flags);
2003	__os.fill(__fill);
2004	__os.precision(__precision);
2005	return __os;
2006	}
2007
2008	template<typename _RealType, typename _CharT, typename _Traits>
2009	std::basic_istream<_CharT, _Traits>&
2010	operator>>(std::basic_istream<_CharT, _Traits>& __is,
2011	lognormal_distribution<_RealType>& __x)
2012	{
2013	using param_type
2014	= typename lognormal_distribution<_RealType>::param_type;
2015	using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2016
2017	const typename __ios_base::fmtflags __flags = __is.flags();
2018	__is.flags(__ios_base::dec | __ios_base::skipws);
2019
2020	_RealType __m, __s;
2021	if (__is >> __m >> __s >> __x._M_nd)
2022	__x.param(param_type(__m, __s));
2023
2024	__is.flags(__flags);
2025	return __is;
2026	}
2027
2028	template<typename _RealType>
2029	template<typename _ForwardIterator,
2030	typename _UniformRandomNumberGenerator>
2031	void
2032	std::chi_squared_distribution<_RealType>::
2033	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2034	_UniformRandomNumberGenerator& __urng)
2035	{
2036	__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2037	while (__f != __t)
2038	*__f++ = 2 * _M_gd(__urng);
2039	}
2040
2041	template<typename _RealType>
2042	template<typename _ForwardIterator,
2043	typename _UniformRandomNumberGenerator>
2044	void
2045	std::chi_squared_distribution<_RealType>::
2046	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2047	_UniformRandomNumberGenerator& __urng,
2048	const typename
2049	std::gamma_distribution<result_type>::param_type& __p)
2050	{
2051	__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2052	while (__f != __t)
2053	*__f++ = 2 * _M_gd(__urng, __p);
2054	}
2055
2056	template<typename _RealType, typename _CharT, typename _Traits>
2057	std::basic_ostream<_CharT, _Traits>&
2058	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2059	const chi_squared_distribution<_RealType>& __x)
2060	{
2061	using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2062
2063	const typename __ios_base::fmtflags __flags = __os.flags();
2064	const _CharT __fill = __os.fill();
2065	const std::streamsize __precision = __os.precision();
2066	const _CharT __space = __os.widen(' ');
2067	__os.flags(__ios_base::scientific | __ios_base::left);
2068	__os.fill(__space);
2069	__os.precision(std::numeric_limits<_RealType>::max_digits10);
2070
2071	__os << __x.n() << __space << __x._M_gd;
2072
2073	__os.flags(__flags);
2074	__os.fill(__fill);
2075	__os.precision(__precision);
2076	return __os;
2077	}
2078
2079	template<typename _RealType, typename _CharT, typename _Traits>
2080	std::basic_istream<_CharT, _Traits>&
2081	operator>>(std::basic_istream<_CharT, _Traits>& __is,
2082	chi_squared_distribution<_RealType>& __x)
2083	{
2084	using param_type
2085	= typename chi_squared_distribution<_RealType>::param_type;
2086	using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2087
2088	const typename __ios_base::fmtflags __flags = __is.flags();
2089	__is.flags(__ios_base::dec | __ios_base::skipws);
2090
2091	_RealType __n;
2092	if (__is >> __n >> __x._M_gd)
2093	__x.param(param_type(__n));
2094
2095	__is.flags(__flags);
2096	return __is;
2097	}
2098
2099
2100	template<typename _RealType>
2101	template<typename _UniformRandomNumberGenerator>
2102	typename cauchy_distribution<_RealType>::result_type
2103	cauchy_distribution<_RealType>::
2104	operator()(_UniformRandomNumberGenerator& __urng,
2105	const param_type& __p)
2106	{
2107	__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2108	__aurng(__urng);
2109	_RealType __u;
2110	do
2111	__u = __aurng();
2112	while (__u == 0.5);
2113
2114	const _RealType __pi = 3.1415926535897932384626433832795029L;
2115	return __p.a() + __p.b() * std::tan(__pi * __u);
2116	}
2117
2118	template<typename _RealType>
2119	template<typename _ForwardIterator,
2120	typename _UniformRandomNumberGenerator>
2121	void
2122	cauchy_distribution<_RealType>::
2123	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2124	_UniformRandomNumberGenerator& __urng,
2125	const param_type& __p)
2126	{
2127	__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2128	const _RealType __pi = 3.1415926535897932384626433832795029L;
2129	__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2130	__aurng(__urng);
2131	while (__f != __t)
2132	{
2133	_RealType __u;
2134	do
2135	__u = __aurng();
2136	while (__u == 0.5);
2137
2138	*__f++ = __p.a() + __p.b() * std::tan(__pi * __u);
2139	}
2140	}
2141
2142	template<typename _RealType, typename _CharT, typename _Traits>
2143	std::basic_ostream<_CharT, _Traits>&
2144	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2145	const cauchy_distribution<_RealType>& __x)
2146	{
2147	using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2148
2149	const typename __ios_base::fmtflags __flags = __os.flags();
2150	const _CharT __fill = __os.fill();
2151	const std::streamsize __precision = __os.precision();
2152	const _CharT __space = __os.widen(' ');
2153	__os.flags(__ios_base::scientific | __ios_base::left);
2154	__os.fill(__space);
2155	__os.precision(std::numeric_limits<_RealType>::max_digits10);
2156
2157	__os << __x.a() << __space << __x.b();
2158
2159	__os.flags(__flags);
2160	__os.fill(__fill);
2161	__os.precision(__precision);
2162	return __os;
2163	}
2164
2165	template<typename _RealType, typename _CharT, typename _Traits>
2166	std::basic_istream<_CharT, _Traits>&
2167	operator>>(std::basic_istream<_CharT, _Traits>& __is,
2168	cauchy_distribution<_RealType>& __x)
2169	{
2170	using param_type = typename cauchy_distribution<_RealType>::param_type;
2171	using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2172
2173	const typename __ios_base::fmtflags __flags = __is.flags();
2174	__is.flags(__ios_base::dec | __ios_base::skipws);
2175
2176	_RealType __a, __b;
2177	if (__is >> __a >> __b)
2178	__x.param(param_type(__a, __b));
2179
2180	__is.flags(__flags);
2181	return __is;
2182	}
2183
2184
2185	template<typename _RealType>
2186	template<typename _ForwardIterator,
2187	typename _UniformRandomNumberGenerator>
2188	void
2189	std::fisher_f_distribution<_RealType>::
2190	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2191	_UniformRandomNumberGenerator& __urng)
2192	{
2193	__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2194	while (__f != __t)
2195	*__f++ = ((_M_gd_x(__urng) * n()) / (_M_gd_y(__urng) * m()));
2196	}
2197
2198	template<typename _RealType>
2199	template<typename _ForwardIterator,
2200	typename _UniformRandomNumberGenerator>
2201	void
2202	std::fisher_f_distribution<_RealType>::
2203	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2204	_UniformRandomNumberGenerator& __urng,
2205	const param_type& __p)
2206	{
2207	__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2208	typedef typename std::gamma_distribution<result_type>::param_type
2209	param_type;
2210	param_type __p1(__p.m() / 2);
2211	param_type __p2(__p.n() / 2);
2212	while (__f != __t)
2213	*__f++ = ((_M_gd_x(__urng, __p1) * n())
2214	/ (_M_gd_y(__urng, __p2) * m()));
2215	}
2216
2217	template<typename _RealType, typename _CharT, typename _Traits>
2218	std::basic_ostream<_CharT, _Traits>&
2219	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2220	const fisher_f_distribution<_RealType>& __x)
2221	{
2222	using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2223
2224	const typename __ios_base::fmtflags __flags = __os.flags();
2225	const _CharT __fill = __os.fill();
2226	const std::streamsize __precision = __os.precision();
2227	const _CharT __space = __os.widen(' ');
2228	__os.flags(__ios_base::scientific | __ios_base::left);
2229	__os.fill(__space);
2230	__os.precision(std::numeric_limits<_RealType>::max_digits10);
2231
2232	__os << __x.m() << __space << __x.n()
2233	<< __space << __x._M_gd_x << __space << __x._M_gd_y;
2234
2235	__os.flags(__flags);
2236	__os.fill(__fill);
2237	__os.precision(__precision);
2238	return __os;
2239	}
2240
2241	template<typename _RealType, typename _CharT, typename _Traits>
2242	std::basic_istream<_CharT, _Traits>&
2243	operator>>(std::basic_istream<_CharT, _Traits>& __is,
2244	fisher_f_distribution<_RealType>& __x)
2245	{
2246	using param_type
2247	= typename fisher_f_distribution<_RealType>::param_type;
2248	using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2249
2250	const typename __ios_base::fmtflags __flags = __is.flags();
2251	__is.flags(__ios_base::dec | __ios_base::skipws);
2252
2253	_RealType __m, __n;
2254	if (__is >> __m >> __n >> __x._M_gd_x >> __x._M_gd_y)
2255	__x.param(param_type(__m, __n));
2256
2257	__is.flags(__flags);
2258	return __is;
2259	}
2260
2261
2262	template<typename _RealType>
2263	template<typename _ForwardIterator,
2264	typename _UniformRandomNumberGenerator>
2265	void
2266	std::student_t_distribution<_RealType>::
2267	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2268	_UniformRandomNumberGenerator& __urng)
2269	{
2270	__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2271	while (__f != __t)
2272	*__f++ = _M_nd(__urng) * std::sqrt(n() / _M_gd(__urng));
2273	}
2274
2275	template<typename _RealType>
2276	template<typename _ForwardIterator,
2277	typename _UniformRandomNumberGenerator>
2278	void
2279	std::student_t_distribution<_RealType>::
2280	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2281	_UniformRandomNumberGenerator& __urng,
2282	const param_type& __p)
2283	{
2284	__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2285	typename std::gamma_distribution<result_type>::param_type
2286	__p2(__p.n() / 2, 2);
2287	while (__f != __t)
2288	*__f++ = _M_nd(__urng) * std::sqrt(__p.n() / _M_gd(__urng, __p2));
2289	}
2290
2291	template<typename _RealType, typename _CharT, typename _Traits>
2292	std::basic_ostream<_CharT, _Traits>&
2293	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2294	const student_t_distribution<_RealType>& __x)
2295	{
2296	using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2297
2298	const typename __ios_base::fmtflags __flags = __os.flags();
2299	const _CharT __fill = __os.fill();
2300	const std::streamsize __precision = __os.precision();
2301	const _CharT __space = __os.widen(' ');
2302	__os.flags(__ios_base::scientific | __ios_base::left);
2303	__os.fill(__space);
2304	__os.precision(std::numeric_limits<_RealType>::max_digits10);
2305
2306	__os << __x.n() << __space << __x._M_nd << __space << __x._M_gd;
2307
2308	__os.flags(__flags);
2309	__os.fill(__fill);
2310	__os.precision(__precision);
2311	return __os;
2312	}
2313
2314	template<typename _RealType, typename _CharT, typename _Traits>
2315	std::basic_istream<_CharT, _Traits>&
2316	operator>>(std::basic_istream<_CharT, _Traits>& __is,
2317	student_t_distribution<_RealType>& __x)
2318	{
2319	using param_type
2320	= typename student_t_distribution<_RealType>::param_type;
2321	using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2322
2323	const typename __ios_base::fmtflags __flags = __is.flags();
2324	__is.flags(__ios_base::dec | __ios_base::skipws);
2325
2326	_RealType __n;
2327	if (__is >> __n >> __x._M_nd >> __x._M_gd)
2328	__x.param(param_type(__n));
2329
2330	__is.flags(__flags);
2331	return __is;
2332	}
2333
2334
2335	template<typename _RealType>
2336	void
2337	gamma_distribution<_RealType>::param_type::
2338	_M_initialize()
2339	{
2340	_M_malpha = _M_alpha < 1.0 ? _M_alpha + _RealType(1.0) : _M_alpha;
2341
2342	const _RealType __a1 = _M_malpha - _RealType(1.0) / _RealType(3.0);
2343	_M_a2 = _RealType(1.0) / std::sqrt(_RealType(9.0) * __a1);
2344	}
2345
2346	/**
2347	* Marsaglia, G. and Tsang, W. W.
2348	* "A Simple Method for Generating Gamma Variables"
2349	* ACM Transactions on Mathematical Software, 26, 3, 363-372, 2000.
2350	*/
2351	template<typename _RealType>
2352	template<typename _UniformRandomNumberGenerator>
2353	typename gamma_distribution<_RealType>::result_type
2354	gamma_distribution<_RealType>::
2355	operator()(_UniformRandomNumberGenerator& __urng,
2356	const param_type& __param)
2357	{
2358	__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2359	__aurng(__urng);
2360
2361	result_type __u, __v, __n;
2362	const result_type __a1 = (__param._M_malpha
2363	- _RealType(1.0) / _RealType(3.0));
2364
2365	do
2366	{
2367	do
2368	{
2369	__n = _M_nd(__urng);
2370	__v = result_type(1.0) + __param._M_a2 * __n;
2371	}
2372	while (__v <= 0.0);
2373
2374	__v = __v * __v * __v;
2375	__u = __aurng();
2376	}
2377	while (__u > result_type(1.0) - 0.0331 * __n * __n * __n * __n
2378	&& (std::log(__u) > (0.5 * __n * __n + __a1
2379	* (1.0 - __v + std::log(__v)))));
2380
2381	if (__param.alpha() == __param._M_malpha)
2382	return __a1 * __v * __param.beta();
2383	else
2384	{
2385	do
2386	__u = __aurng();
2387	while (__u == 0.0);
2388
2389	return (std::pow(__u, result_type(1.0) / __param.alpha())
2390	* __a1 * __v * __param.beta());
2391	}
2392	}
2393
2394	template<typename _RealType>
2395	template<typename _ForwardIterator,
2396	typename _UniformRandomNumberGenerator>
2397	void
2398	gamma_distribution<_RealType>::
2399	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2400	_UniformRandomNumberGenerator& __urng,
2401	const param_type& __param)
2402	{
2403	__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2404	__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2405	__aurng(__urng);
2406
2407	result_type __u, __v, __n;
2408	const result_type __a1 = (__param._M_malpha
2409	- _RealType(1.0) / _RealType(3.0));
2410
2411	if (__param.alpha() == __param._M_malpha)
2412	while (__f != __t)
2413	{
2414	do
2415	{
2416	do
2417	{
2418	__n = _M_nd(__urng);
2419	__v = result_type(1.0) + __param._M_a2 * __n;
2420	}
2421	while (__v <= 0.0);
2422
2423	__v = __v * __v * __v;
2424	__u = __aurng();
2425	}
2426	while (__u > result_type(1.0) - 0.0331 * __n * __n * __n * __n
2427	&& (std::log(__u) > (0.5 * __n * __n + __a1
2428	* (1.0 - __v + std::log(__v)))));
2429
2430	*__f++ = __a1 * __v * __param.beta();
2431	}
2432	else
2433	while (__f != __t)
2434	{
2435	do
2436	{
2437	do
2438	{
2439	__n = _M_nd(__urng);
2440	__v = result_type(1.0) + __param._M_a2 * __n;
2441	}
2442	while (__v <= 0.0);
2443
2444	__v = __v * __v * __v;
2445	__u = __aurng();
2446	}
2447	while (__u > result_type(1.0) - 0.0331 * __n * __n * __n * __n
2448	&& (std::log(__u) > (0.5 * __n * __n + __a1
2449	* (1.0 - __v + std::log(__v)))));
2450
2451	do
2452	__u = __aurng();
2453	while (__u == 0.0);
2454
2455	*__f++ = (std::pow(__u, result_type(1.0) / __param.alpha())
2456	* __a1 * __v * __param.beta());
2457	}
2458	}
2459
2460	template<typename _RealType, typename _CharT, typename _Traits>
2461	std::basic_ostream<_CharT, _Traits>&
2462	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2463	const gamma_distribution<_RealType>& __x)
2464	{
2465	using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2466
2467	const typename __ios_base::fmtflags __flags = __os.flags();
2468	const _CharT __fill = __os.fill();
2469	const std::streamsize __precision = __os.precision();
2470	const _CharT __space = __os.widen(' ');
2471	__os.flags(__ios_base::scientific | __ios_base::left);
2472	__os.fill(__space);
2473	__os.precision(std::numeric_limits<_RealType>::max_digits10);
2474
2475	__os << __x.alpha() << __space << __x.beta()
2476	<< __space << __x._M_nd;
2477
2478	__os.flags(__flags);
2479	__os.fill(__fill);
2480	__os.precision(__precision);
2481	return __os;
2482	}
2483
2484	template<typename _RealType, typename _CharT, typename _Traits>
2485	std::basic_istream<_CharT, _Traits>&
2486	operator>>(std::basic_istream<_CharT, _Traits>& __is,
2487	gamma_distribution<_RealType>& __x)
2488	{
2489	using param_type = typename gamma_distribution<_RealType>::param_type;
2490	using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2491
2492	const typename __ios_base::fmtflags __flags = __is.flags();
2493	__is.flags(__ios_base::dec | __ios_base::skipws);
2494
2495	_RealType __alpha_val, __beta_val;
2496	if (__is >> __alpha_val >> __beta_val >> __x._M_nd)
2497	__x.param(param_type(__alpha_val, __beta_val));
2498
2499	__is.flags(__flags);
2500	return __is;
2501	}
2502
2503
2504	template<typename _RealType>
2505	template<typename _UniformRandomNumberGenerator>
2506	typename weibull_distribution<_RealType>::result_type
2507	weibull_distribution<_RealType>::
2508	operator()(_UniformRandomNumberGenerator& __urng,
2509	const param_type& __p)
2510	{
2511	__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2512	__aurng(__urng);
2513	return __p.b() * std::pow(-std::log(result_type(1) - __aurng()),
2514	result_type(1) / __p.a());
2515	}
2516
2517	template<typename _RealType>
2518	template<typename _ForwardIterator,
2519	typename _UniformRandomNumberGenerator>
2520	void
2521	weibull_distribution<_RealType>::
2522	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2523	_UniformRandomNumberGenerator& __urng,
2524	const param_type& __p)
2525	{
2526	__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2527	__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2528	__aurng(__urng);
2529	auto __inv_a = result_type(1) / __p.a();
2530
2531	while (__f != __t)
2532	*__f++ = __p.b() * std::pow(-std::log(result_type(1) - __aurng()),
2533	__inv_a);
2534	}
2535
2536	template<typename _RealType, typename _CharT, typename _Traits>
2537	std::basic_ostream<_CharT, _Traits>&
2538	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2539	const weibull_distribution<_RealType>& __x)
2540	{
2541	using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2542
2543	const typename __ios_base::fmtflags __flags = __os.flags();
2544	const _CharT __fill = __os.fill();
2545	const std::streamsize __precision = __os.precision();
2546	const _CharT __space = __os.widen(' ');
2547	__os.flags(__ios_base::scientific | __ios_base::left);
2548	__os.fill(__space);
2549	__os.precision(std::numeric_limits<_RealType>::max_digits10);
2550
2551	__os << __x.a() << __space << __x.b();
2552
2553	__os.flags(__flags);
2554	__os.fill(__fill);
2555	__os.precision(__precision);
2556	return __os;
2557	}
2558
2559	template<typename _RealType, typename _CharT, typename _Traits>
2560	std::basic_istream<_CharT, _Traits>&
2561	operator>>(std::basic_istream<_CharT, _Traits>& __is,
2562	weibull_distribution<_RealType>& __x)
2563	{
2564	using param_type = typename weibull_distribution<_RealType>::param_type;
2565	using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2566
2567	const typename __ios_base::fmtflags __flags = __is.flags();
2568	__is.flags(__ios_base::dec | __ios_base::skipws);
2569
2570	_RealType __a, __b;
2571	if (__is >> __a >> __b)
2572	__x.param(param_type(__a, __b));
2573
2574	__is.flags(__flags);
2575	return __is;
2576	}
2577
2578
2579	template<typename _RealType>
2580	template<typename _UniformRandomNumberGenerator>
2581	typename extreme_value_distribution<_RealType>::result_type
2582	extreme_value_distribution<_RealType>::
2583	operator()(_UniformRandomNumberGenerator& __urng,
2584	const param_type& __p)
2585	{
2586	__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2587	__aurng(__urng);
2588	return __p.a() - __p.b() * std::log(-std::log(result_type(1)
2589	- __aurng()));
2590	}
2591
2592	template<typename _RealType>
2593	template<typename _ForwardIterator,
2594	typename _UniformRandomNumberGenerator>
2595	void
2596	extreme_value_distribution<_RealType>::
2597	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2598	_UniformRandomNumberGenerator& __urng,
2599	const param_type& __p)
2600	{
2601	__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2602	__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2603	__aurng(__urng);
2604
2605	while (__f != __t)
2606	*__f++ = __p.a() - __p.b() * std::log(-std::log(result_type(1)
2607	- __aurng()));
2608	}
2609
2610	template<typename _RealType, typename _CharT, typename _Traits>
2611	std::basic_ostream<_CharT, _Traits>&
2612	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2613	const extreme_value_distribution<_RealType>& __x)
2614	{
2615	using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2616
2617	const typename __ios_base::fmtflags __flags = __os.flags();
2618	const _CharT __fill = __os.fill();
2619	const std::streamsize __precision = __os.precision();
2620	const _CharT __space = __os.widen(' ');
2621	__os.flags(__ios_base::scientific | __ios_base::left);
2622	__os.fill(__space);
2623	__os.precision(std::numeric_limits<_RealType>::max_digits10);
2624
2625	__os << __x.a() << __space << __x.b();
2626
2627	__os.flags(__flags);
2628	__os.fill(__fill);
2629	__os.precision(__precision);
2630	return __os;
2631	}
2632
2633	template<typename _RealType, typename _CharT, typename _Traits>
2634	std::basic_istream<_CharT, _Traits>&
2635	operator>>(std::basic_istream<_CharT, _Traits>& __is,
2636	extreme_value_distribution<_RealType>& __x)
2637	{
2638	using param_type
2639	= typename extreme_value_distribution<_RealType>::param_type;
2640	using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2641
2642	const typename __ios_base::fmtflags __flags = __is.flags();
2643	__is.flags(__ios_base::dec | __ios_base::skipws);
2644
2645	_RealType __a, __b;
2646	if (__is >> __a >> __b)
2647	__x.param(param_type(__a, __b));
2648
2649	__is.flags(__flags);
2650	return __is;
2651	}
2652
2653
2654	template<typename _IntType>
2655	void
2656	discrete_distribution<_IntType>::param_type::
2657	_M_initialize()
2658	{
2659	if (_M_prob.size() < 2)
2660	{
2661	_M_prob.clear();
2662	return;
2663	}
2664
2665	const double __sum = std::accumulate(first: _M_prob.begin(),
2666	last: _M_prob.end(), init: 0.0);
2667	__glibcxx_assert(__sum > 0);
2668	// Now normalize the probabilites.
2669	__detail::__normalize(first: _M_prob.begin(), last: _M_prob.end(), result: _M_prob.begin(),
2670	factor: __sum);
2671	// Accumulate partial sums.
2672	_M_cp.reserve(n: _M_prob.size());
2673	std::partial_sum(first: _M_prob.begin(), last: _M_prob.end(),
2674	result: std::back_inserter(x&: _M_cp));
2675	// Make sure the last cumulative probability is one.
2676	_M_cp[_M_cp.size() - 1] = 1.0;
2677	}
2678
2679	template<typename _IntType>
2680	template<typename _Func>
2681	discrete_distribution<_IntType>::param_type::
2682	param_type(size_t __nw, double __xmin, double __xmax, _Func __fw)
2683	: _M_prob(), _M_cp()
2684	{
2685	const size_t __n = __nw == 0 ? 1 : __nw;
2686	const double __delta = (__xmax - __xmin) / __n;
2687
2688	_M_prob.reserve(__n);
2689	for (size_t __k = 0; __k < __nw; ++__k)
2690	_M_prob.push_back(__fw(__xmin + __k * __delta + 0.5 * __delta));
2691
2692	_M_initialize();
2693	}
2694
2695	template<typename _IntType>
2696	template<typename _UniformRandomNumberGenerator>
2697	typename discrete_distribution<_IntType>::result_type
2698	discrete_distribution<_IntType>::
2699	operator()(_UniformRandomNumberGenerator& __urng,
2700	const param_type& __param)
2701	{
2702	if (__param._M_cp.empty())
2703	return result_type(0);
2704
2705	__detail::_Adaptor<_UniformRandomNumberGenerator, double>
2706	__aurng(__urng);
2707
2708	const double __p = __aurng();
2709	auto __pos = std::lower_bound(__param._M_cp.begin(),
2710	__param._M_cp.end(), __p);
2711
2712	return __pos - __param._M_cp.begin();
2713	}
2714
2715	template<typename _IntType>
2716	template<typename _ForwardIterator,
2717	typename _UniformRandomNumberGenerator>
2718	void
2719	discrete_distribution<_IntType>::
2720	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2721	_UniformRandomNumberGenerator& __urng,
2722	const param_type& __param)
2723	{
2724	__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2725
2726	if (__param._M_cp.empty())
2727	{
2728	while (__f != __t)
2729	*__f++ = result_type(0);
2730	return;
2731	}
2732
2733	__detail::_Adaptor<_UniformRandomNumberGenerator, double>
2734	__aurng(__urng);
2735
2736	while (__f != __t)
2737	{
2738	const double __p = __aurng();
2739	auto __pos = std::lower_bound(__param._M_cp.begin(),
2740	__param._M_cp.end(), __p);
2741
2742	*__f++ = __pos - __param._M_cp.begin();
2743	}
2744	}
2745
2746	template<typename _IntType, typename _CharT, typename _Traits>
2747	std::basic_ostream<_CharT, _Traits>&
2748	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2749	const discrete_distribution<_IntType>& __x)
2750	{
2751	using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2752
2753	const typename __ios_base::fmtflags __flags = __os.flags();
2754	const _CharT __fill = __os.fill();
2755	const std::streamsize __precision = __os.precision();
2756	const _CharT __space = __os.widen(' ');
2757	__os.flags(__ios_base::scientific | __ios_base::left);
2758	__os.fill(__space);
2759	__os.precision(std::numeric_limits<double>::max_digits10);
2760
2761	std::vector<double> __prob = __x.probabilities();
2762	__os << __prob.size();
2763	for (auto __dit = __prob.begin(); __dit != __prob.end(); ++__dit)
2764	__os << __space << *__dit;
2765
2766	__os.flags(__flags);
2767	__os.fill(__fill);
2768	__os.precision(__precision);
2769	return __os;
2770	}
2771
2772	namespace __detail
2773	{
2774	template<typename _ValT, typename _CharT, typename _Traits>
2775	basic_istream<_CharT, _Traits>&
2776	__extract_params(basic_istream<_CharT, _Traits>& __is,
2777	vector<_ValT>& __vals, size_t __n)
2778	{
2779	__vals.reserve(__n);
2780	while (__n--)
2781	{
2782	_ValT __val;
2783	if (__is >> __val)
2784	__vals.push_back(__val);
2785	else
2786	break;
2787	}
2788	return __is;
2789	}
2790	} // namespace __detail
2791
2792	template<typename _IntType, typename _CharT, typename _Traits>
2793	std::basic_istream<_CharT, _Traits>&
2794	operator>>(std::basic_istream<_CharT, _Traits>& __is,
2795	discrete_distribution<_IntType>& __x)
2796	{
2797	using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
2798
2799	const typename __ios_base::fmtflags __flags = __is.flags();
2800	__is.flags(__ios_base::dec | __ios_base::skipws);
2801
2802	size_t __n;
2803	if (__is >> __n)
2804	{
2805	std::vector<double> __prob_vec;
2806	if (__detail::__extract_params(__is, __prob_vec, __n))
2807	__x.param({__prob_vec.begin(), __prob_vec.end()});
2808	}
2809
2810	__is.flags(__flags);
2811	return __is;
2812	}
2813
2814
2815	template<typename _RealType>
2816	void
2817	piecewise_constant_distribution<_RealType>::param_type::
2818	_M_initialize()
2819	{
2820	if (_M_int.size() < 2
2821	|| (_M_int.size() == 2
2822	&& _M_int[0] == _RealType(0)
2823	&& _M_int[1] == _RealType(1)))
2824	{
2825	_M_int.clear();
2826	_M_den.clear();
2827	return;
2828	}
2829
2830	const double __sum = std::accumulate(first: _M_den.begin(),
2831	last: _M_den.end(), init: 0.0);
2832	__glibcxx_assert(__sum > 0);
2833
2834	__detail::__normalize(first: _M_den.begin(), last: _M_den.end(), result: _M_den.begin(),
2835	factor: __sum);
2836
2837	_M_cp.reserve(n: _M_den.size());
2838	std::partial_sum(first: _M_den.begin(), last: _M_den.end(),
2839	result: std::back_inserter(x&: _M_cp));
2840
2841	// Make sure the last cumulative probability is one.
2842	_M_cp[_M_cp.size() - 1] = 1.0;
2843
2844	for (size_t __k = 0; __k < _M_den.size(); ++__k)
2845	_M_den[__k] /= _M_int[__k + 1] - _M_int[__k];
2846	}
2847
2848	template<typename _RealType>
2849	template<typename _InputIteratorB, typename _InputIteratorW>
2850	piecewise_constant_distribution<_RealType>::param_type::
2851	param_type(_InputIteratorB __bbegin,
2852	_InputIteratorB __bend,
2853	_InputIteratorW __wbegin)
2854	: _M_int(), _M_den(), _M_cp()
2855	{
2856	if (__bbegin != __bend)
2857	{
2858	for (;;)
2859	{
2860	_M_int.push_back(*__bbegin);
2861	++__bbegin;
2862	if (__bbegin == __bend)
2863	break;
2864
2865	_M_den.push_back(*__wbegin);
2866	++__wbegin;
2867	}
2868	}
2869
2870	_M_initialize();
2871	}
2872
2873	template<typename _RealType>
2874	template<typename _Func>
2875	piecewise_constant_distribution<_RealType>::param_type::
2876	param_type(initializer_list<_RealType> __bl, _Func __fw)
2877	: _M_int(), _M_den(), _M_cp()
2878	{
2879	_M_int.reserve(__bl.size());
2880	for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter)
2881	_M_int.push_back(*__biter);
2882
2883	_M_den.reserve(n: _M_int.size() - 1);
2884	for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
2885	_M_den.push_back(__fw(0.5 * (_M_int[__k + 1] + _M_int[__k])));
2886
2887	_M_initialize();
2888	}
2889
2890	template<typename _RealType>
2891	template<typename _Func>
2892	piecewise_constant_distribution<_RealType>::param_type::
2893	param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
2894	: _M_int(), _M_den(), _M_cp()
2895	{
2896	const size_t __n = __nw == 0 ? 1 : __nw;
2897	const _RealType __delta = (__xmax - __xmin) / __n;
2898
2899	_M_int.reserve(__n + 1);
2900	for (size_t __k = 0; __k <= __nw; ++__k)
2901	_M_int.push_back(__xmin + __k * __delta);
2902
2903	_M_den.reserve(__n);
2904	for (size_t __k = 0; __k < __nw; ++__k)
2905	_M_den.push_back(__fw(_M_int[__k] + 0.5 * __delta));
2906
2907	_M_initialize();
2908	}
2909
2910	template<typename _RealType>
2911	template<typename _UniformRandomNumberGenerator>
2912	typename piecewise_constant_distribution<_RealType>::result_type
2913	piecewise_constant_distribution<_RealType>::
2914	operator()(_UniformRandomNumberGenerator& __urng,
2915	const param_type& __param)
2916	{
2917	__detail::_Adaptor<_UniformRandomNumberGenerator, double>
2918	__aurng(__urng);
2919
2920	const double __p = __aurng();
2921	if (__param._M_cp.empty())
2922	return __p;
2923
2924	auto __pos = std::lower_bound(__param._M_cp.begin(),
2925	__param._M_cp.end(), __p);
2926	const size_t __i = __pos - __param._M_cp.begin();
2927
2928	const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
2929
2930	return __param._M_int[__i] + (__p - __pref) / __param._M_den[__i];
2931	}
2932
2933	template<typename _RealType>
2934	template<typename _ForwardIterator,
2935	typename _UniformRandomNumberGenerator>
2936	void
2937	piecewise_constant_distribution<_RealType>::
2938	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2939	_UniformRandomNumberGenerator& __urng,
2940	const param_type& __param)
2941	{
2942	__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
2943	__detail::_Adaptor<_UniformRandomNumberGenerator, double>
2944	__aurng(__urng);
2945
2946	if (__param._M_cp.empty())
2947	{
2948	while (__f != __t)
2949	*__f++ = __aurng();
2950	return;
2951	}
2952
2953	while (__f != __t)
2954	{
2955	const double __p = __aurng();
2956
2957	auto __pos = std::lower_bound(__param._M_cp.begin(),
2958	__param._M_cp.end(), __p);
2959	const size_t __i = __pos - __param._M_cp.begin();
2960
2961	const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
2962
2963	*__f++ = (__param._M_int[__i]
2964	+ (__p - __pref) / __param._M_den[__i]);
2965	}
2966	}
2967
2968	template<typename _RealType, typename _CharT, typename _Traits>
2969	std::basic_ostream<_CharT, _Traits>&
2970	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2971	const piecewise_constant_distribution<_RealType>& __x)
2972	{
2973	using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
2974
2975	const typename __ios_base::fmtflags __flags = __os.flags();
2976	const _CharT __fill = __os.fill();
2977	const std::streamsize __precision = __os.precision();
2978	const _CharT __space = __os.widen(' ');
2979	__os.flags(__ios_base::scientific | __ios_base::left);
2980	__os.fill(__space);
2981	__os.precision(std::numeric_limits<_RealType>::max_digits10);
2982
2983	std::vector<_RealType> __int = __x.intervals();
2984	__os << __int.size() - 1;
2985
2986	for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit)
2987	__os << __space << *__xit;
2988
2989	std::vector<double> __den = __x.densities();
2990	for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit)
2991	__os << __space << *__dit;
2992
2993	__os.flags(__flags);
2994	__os.fill(__fill);
2995	__os.precision(__precision);
2996	return __os;
2997	}
2998
2999	template<typename _RealType, typename _CharT, typename _Traits>
3000	std::basic_istream<_CharT, _Traits>&
3001	operator>>(std::basic_istream<_CharT, _Traits>& __is,
3002	piecewise_constant_distribution<_RealType>& __x)
3003	{
3004	using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
3005
3006	const typename __ios_base::fmtflags __flags = __is.flags();
3007	__is.flags(__ios_base::dec | __ios_base::skipws);
3008
3009	size_t __n;
3010	if (__is >> __n)
3011	{
3012	std::vector<_RealType> __int_vec;
3013	if (__detail::__extract_params(__is, __int_vec, __n + 1))
3014	{
3015	std::vector<double> __den_vec;
3016	if (__detail::__extract_params(__is, __den_vec, __n))
3017	{
3018	__x.param({ __int_vec.begin(), __int_vec.end(),
3019	__den_vec.begin() });
3020	}
3021	}
3022	}
3023
3024	__is.flags(__flags);
3025	return __is;
3026	}
3027
3028
3029	template<typename _RealType>
3030	void
3031	piecewise_linear_distribution<_RealType>::param_type::
3032	_M_initialize()
3033	{
3034	if (_M_int.size() < 2
3035	|| (_M_int.size() == 2
3036	&& _M_int[0] == _RealType(0)
3037	&& _M_int[1] == _RealType(1)
3038	&& _M_den[0] == _M_den[1]))
3039	{
3040	_M_int.clear();
3041	_M_den.clear();
3042	return;
3043	}
3044
3045	double __sum = 0.0;
3046	_M_cp.reserve(n: _M_int.size() - 1);
3047	_M_m.reserve(n: _M_int.size() - 1);
3048	for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
3049	{
3050	const _RealType __delta = _M_int[__k + 1] - _M_int[__k];
3051	__sum += 0.5 * (_M_den[__k + 1] + _M_den[__k]) * __delta;
3052	_M_cp.push_back(x: __sum);
3053	_M_m.push_back((_M_den[__k + 1] - _M_den[__k]) / __delta);
3054	}
3055	__glibcxx_assert(__sum > 0);
3056
3057	// Now normalize the densities...
3058	__detail::__normalize(first: _M_den.begin(), last: _M_den.end(), result: _M_den.begin(),
3059	factor: __sum);
3060	// ... and partial sums...
3061	__detail::__normalize(first: _M_cp.begin(), last: _M_cp.end(), result: _M_cp.begin(), factor: __sum);
3062	// ... and slopes.
3063	__detail::__normalize(first: _M_m.begin(), last: _M_m.end(), result: _M_m.begin(), factor: __sum);
3064
3065	// Make sure the last cumulative probablility is one.
3066	_M_cp[_M_cp.size() - 1] = 1.0;
3067	}
3068
3069	template<typename _RealType>
3070	template<typename _InputIteratorB, typename _InputIteratorW>
3071	piecewise_linear_distribution<_RealType>::param_type::
3072	param_type(_InputIteratorB __bbegin,
3073	_InputIteratorB __bend,
3074	_InputIteratorW __wbegin)
3075	: _M_int(), _M_den(), _M_cp(), _M_m()
3076	{
3077	for (; __bbegin != __bend; ++__bbegin, ++__wbegin)
3078	{
3079	_M_int.push_back(*__bbegin);
3080	_M_den.push_back(*__wbegin);
3081	}
3082
3083	_M_initialize();
3084	}
3085
3086	template<typename _RealType>
3087	template<typename _Func>
3088	piecewise_linear_distribution<_RealType>::param_type::
3089	param_type(initializer_list<_RealType> __bl, _Func __fw)
3090	: _M_int(), _M_den(), _M_cp(), _M_m()
3091	{
3092	_M_int.reserve(__bl.size());
3093	_M_den.reserve(n: __bl.size());
3094	for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter)
3095	{
3096	_M_int.push_back(*__biter);
3097	_M_den.push_back(__fw(*__biter));
3098	}
3099
3100	_M_initialize();
3101	}
3102
3103	template<typename _RealType>
3104	template<typename _Func>
3105	piecewise_linear_distribution<_RealType>::param_type::
3106	param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
3107	: _M_int(), _M_den(), _M_cp(), _M_m()
3108	{
3109	const size_t __n = __nw == 0 ? 1 : __nw;
3110	const _RealType __delta = (__xmax - __xmin) / __n;
3111
3112	_M_int.reserve(__n + 1);
3113	_M_den.reserve(n: __n + 1);
3114	for (size_t __k = 0; __k <= __nw; ++__k)
3115	{
3116	_M_int.push_back(__xmin + __k * __delta);
3117	_M_den.push_back(__fw(_M_int[__k] + __delta));
3118	}
3119
3120	_M_initialize();
3121	}
3122
3123	template<typename _RealType>
3124	template<typename _UniformRandomNumberGenerator>
3125	typename piecewise_linear_distribution<_RealType>::result_type
3126	piecewise_linear_distribution<_RealType>::
3127	operator()(_UniformRandomNumberGenerator& __urng,
3128	const param_type& __param)
3129	{
3130	__detail::_Adaptor<_UniformRandomNumberGenerator, double>
3131	__aurng(__urng);
3132
3133	const double __p = __aurng();
3134	if (__param._M_cp.empty())
3135	return __p;
3136
3137	auto __pos = std::lower_bound(__param._M_cp.begin(),
3138	__param._M_cp.end(), __p);
3139	const size_t __i = __pos - __param._M_cp.begin();
3140
3141	const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
3142
3143	const double __a = 0.5 * __param._M_m[__i];
3144	const double __b = __param._M_den[__i];
3145	const double __cm = __p - __pref;
3146
3147	_RealType __x = __param._M_int[__i];
3148	if (__a == 0)
3149	__x += __cm / __b;
3150	else
3151	{
3152	const double __d = __b * __b + 4.0 * __a * __cm;
3153	__x += 0.5 * (std::sqrt(x: __d) - __b) / __a;
3154	}
3155
3156	return __x;
3157	}
3158
3159	template<typename _RealType>
3160	template<typename _ForwardIterator,
3161	typename _UniformRandomNumberGenerator>
3162	void
3163	piecewise_linear_distribution<_RealType>::
3164	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
3165	_UniformRandomNumberGenerator& __urng,
3166	const param_type& __param)
3167	{
3168	__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
3169	// We could duplicate everything from operator()...
3170	while (__f != __t)
3171	*__f++ = this->operator()(__urng, __param);
3172	}
3173
3174	template<typename _RealType, typename _CharT, typename _Traits>
3175	std::basic_ostream<_CharT, _Traits>&
3176	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
3177	const piecewise_linear_distribution<_RealType>& __x)
3178	{
3179	using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
3180
3181	const typename __ios_base::fmtflags __flags = __os.flags();
3182	const _CharT __fill = __os.fill();
3183	const std::streamsize __precision = __os.precision();
3184	const _CharT __space = __os.widen(' ');
3185	__os.flags(__ios_base::scientific | __ios_base::left);
3186	__os.fill(__space);
3187	__os.precision(std::numeric_limits<_RealType>::max_digits10);
3188
3189	std::vector<_RealType> __int = __x.intervals();
3190	__os << __int.size() - 1;
3191
3192	for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit)
3193	__os << __space << *__xit;
3194
3195	std::vector<double> __den = __x.densities();
3196	for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit)
3197	__os << __space << *__dit;
3198
3199	__os.flags(__flags);
3200	__os.fill(__fill);
3201	__os.precision(__precision);
3202	return __os;
3203	}
3204
3205	template<typename _RealType, typename _CharT, typename _Traits>
3206	std::basic_istream<_CharT, _Traits>&
3207	operator>>(std::basic_istream<_CharT, _Traits>& __is,
3208	piecewise_linear_distribution<_RealType>& __x)
3209	{
3210	using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
3211
3212	const typename __ios_base::fmtflags __flags = __is.flags();
3213	__is.flags(__ios_base::dec | __ios_base::skipws);
3214
3215	size_t __n;
3216	if (__is >> __n)
3217	{
3218	vector<_RealType> __int_vec;
3219	if (__detail::__extract_params(__is, __int_vec, __n + 1))
3220	{
3221	vector<double> __den_vec;
3222	if (__detail::__extract_params(__is, __den_vec, __n + 1))
3223	{
3224	__x.param({ __int_vec.begin(), __int_vec.end(),
3225	__den_vec.begin() });
3226	}
3227	}
3228	}
3229	__is.flags(__flags);
3230	return __is;
3231	}
3232
3233
3234	template<typename _IntType, typename>
3235	seed_seq::seed_seq(std::initializer_list<_IntType> __il)
3236	{
3237	_M_v.reserve(n: __il.size());
3238	for (auto __iter = __il.begin(); __iter != __il.end(); ++__iter)
3239	_M_v.push_back(__detail::__mod<result_type,
3240	__detail::_Shift<result_type, 32>::__value>(*__iter));
3241	}
3242
3243	template<typename _InputIterator>
3244	seed_seq::seed_seq(_InputIterator __begin, _InputIterator __end)
3245	{
3246	if _GLIBCXX17_CONSTEXPR (__is_random_access_iter<_InputIterator>::value)
3247	_M_v.reserve(n: std::distance(__begin, __end));
3248
3249	for (_InputIterator __iter = __begin; __iter != __end; ++__iter)
3250	_M_v.push_back(__detail::__mod<result_type,
3251	__detail::_Shift<result_type, 32>::__value>(*__iter));
3252	}
3253
3254	template<typename _RandomAccessIterator>
3255	void
3256	seed_seq::generate(_RandomAccessIterator __begin,
3257	_RandomAccessIterator __end)
3258	{
3259	typedef typename iterator_traits<_RandomAccessIterator>::value_type
3260	_Type;
3261
3262	if (__begin == __end)
3263	return;
3264
3265	std::fill(__begin, __end, _Type(0x8b8b8b8bu));
3266
3267	const size_t __n = __end - __begin;
3268	const size_t __s = _M_v.size();
3269	const size_t __t = (__n >= 623) ? 11
3270	: (__n >= 68) ? 7
3271	: (__n >= 39) ? 5
3272	: (__n >= 7) ? 3
3273	: (__n - 1) / 2;
3274	const size_t __p = (__n - __t) / 2;
3275	const size_t __q = __p + __t;
3276	const size_t __m = std::max(a: size_t(__s + 1), b: __n);
3277
3278	#ifndef __UINT32_TYPE__
3279	struct _Up
3280	{
3281	_Up(uint_least32_t v) : _M_v(v & 0xffffffffu) { }
3282
3283	operator uint_least32_t() const { return _M_v; }
3284
3285	uint_least32_t _M_v;
3286	};
3287	using uint32_t = _Up;
3288	#endif
3289
3290	// k == 0, every element in [begin,end) equals 0x8b8b8b8bu
3291	{
3292	uint32_t __r1 = 1371501266u;
3293	uint32_t __r2 = __r1 + __s;
3294	__begin[__p] += __r1;
3295	__begin[__q] = (uint32_t)__begin[__q] + __r2;
3296	__begin[0] = __r2;
3297	}
3298
3299	for (size_t __k = 1; __k <= __s; ++__k)
3300	{
3301	const size_t __kn = __k % __n;
3302	const size_t __kpn = (__k + __p) % __n;
3303	const size_t __kqn = (__k + __q) % __n;
3304	uint32_t __arg = (__begin[__kn]
3305	^ __begin[__kpn]
3306	^ __begin[(__k - 1) % __n]);
3307	uint32_t __r1 = 1664525u * (__arg ^ (__arg >> 27));
3308	uint32_t __r2 = __r1 + (uint32_t)__kn + _M_v[__k - 1];
3309	__begin[__kpn] = (uint32_t)__begin[__kpn] + __r1;
3310	__begin[__kqn] = (uint32_t)__begin[__kqn] + __r2;
3311	__begin[__kn] = __r2;
3312	}
3313
3314	for (size_t __k = __s + 1; __k < __m; ++__k)
3315	{
3316	const size_t __kn = __k % __n;
3317	const size_t __kpn = (__k + __p) % __n;
3318	const size_t __kqn = (__k + __q) % __n;
3319	uint32_t __arg = (__begin[__kn]
3320	^ __begin[__kpn]
3321	^ __begin[(__k - 1) % __n]);
3322	uint32_t __r1 = 1664525u * (__arg ^ (__arg >> 27));
3323	uint32_t __r2 = __r1 + (uint32_t)__kn;
3324	__begin[__kpn] = (uint32_t)__begin[__kpn] + __r1;
3325	__begin[__kqn] = (uint32_t)__begin[__kqn] + __r2;
3326	__begin[__kn] = __r2;
3327	}
3328
3329	for (size_t __k = __m; __k < __m + __n; ++__k)
3330	{
3331	const size_t __kn = __k % __n;
3332	const size_t __kpn = (__k + __p) % __n;
3333	const size_t __kqn = (__k + __q) % __n;
3334	uint32_t __arg = (__begin[__kn]
3335	+ __begin[__kpn]
3336	+ __begin[(__k - 1) % __n]);
3337	uint32_t __r3 = 1566083941u * (__arg ^ (__arg >> 27));
3338	uint32_t __r4 = __r3 - __kn;
3339	__begin[__kpn] ^= __r3;
3340	__begin[__kqn] ^= __r4;
3341	__begin[__kn] = __r4;
3342	}
3343	}
3344
3345	template<typename _RealType, size_t __bits,
3346	typename _UniformRandomNumberGenerator>
3347	_RealType
3348	generate_canonical(_UniformRandomNumberGenerator& __urng)
3349	{
3350	static_assert(std::is_floating_point<_RealType>::value,
3351	"template argument must be a floating point type");
3352
3353	const size_t __b
3354	= std::min(a: static_cast<size_t>(std::numeric_limits<_RealType>::digits),
3355	b: __bits);
3356	const long double __r = static_cast<long double>(__urng.max())
3357	- static_cast<long double>(__urng.min()) + 1.0L;
3358	const size_t __log2r = std::log(x: __r) / std::log(x: 2.0L);
3359	const size_t __m = std::max<size_t>(a: 1UL,
3360	b: (__b + __log2r - 1UL) / __log2r);
3361	_RealType __ret;
3362	_RealType __sum = _RealType(0);
3363	_RealType __tmp = _RealType(1);
3364	for (size_t __k = __m; __k != 0; --__k)
3365	{
3366	__sum += _RealType(__urng() - __urng.min()) * __tmp;
3367	__tmp *= __r;
3368	}
3369	__ret = __sum / __tmp;
3370	if (__builtin_expect(__ret >= _RealType(1), 0))
3371	{
3372	#if _GLIBCXX_USE_C99_MATH_TR1
3373	__ret = std::nextafter(_RealType(1), _RealType(0));
3374	#else
3375	__ret = _RealType(1)
3376	- std::numeric_limits<_RealType>::epsilon() / _RealType(2);
3377	#endif
3378	}
3379	return __ret;
3380	}
3381
3382	_GLIBCXX_END_NAMESPACE_VERSION
3383	} // namespace
3384
3385	#endif
3386