<random>: Implement modified ziggurat algorithm for normal_distribution

[statementreply]

This PR implements modified ziggurat algorithm for normal_distribution, based on McFarland, C. D. (2015). A modified ziggurat algorithm for generating exponentially and normally distributed pseudorandom numbers.

The main differences between this PR and the paper:

• The layers are moved upwards to reduce the area of non-rectangular regions, creating a bottommost layer of which the size of the rectangular region is smaller than 1/N.
• Different criteria are used for the always sampled and always rejected regions of ziggurat overhangs.

Benchmark results for x64 on Intel Core Ultra 7 265KF:

Benchmark	Time	Speed Up
BM_Generator<std::mt19937_64>	2.40 ns
BM_Generator<b_r::mt19937_64>	1.60 ns
BM_Distribution<std::mt19937_64, std_old::normal_distribution<double>>	7.00 ns	1.00
BM_Distribution<std::mt19937_64, std_new::normal_distribution<double>>	3.15 ns	2.22
BM_Distribution<std::mt19937_64, boost_r::normal_distribution<double>>	3.94 ns	1.78
BM_Distribution<b_r::mt19937_64, boost_r::normal_distribution<double>>	3.57 ns	1.96
BM_Generator<std::mt19937>	2.43 ns
BM_Generator<b_r::mt19937>	1.40 ns
BM_Distribution<std::mt19937, std_old::normal_distribution<float>>	6.77 ns	1.00
BM_Distribution<std::mt19937, std_new::normal_distribution<float>>	3.14 ns	2.16
BM_Distribution<std::mt19937, boost_r::normal_distribution<float>>	3.98 ns	1.70
BM_Distribution<b_r::mt19937, boost_r::normal_distribution<float>>	3.14 ns	2.16

Removed tests\GH_004618_normal_distribution_avoids_resets which became moot.

Fixes #1003.

[frederick-vs-ja]

Can we remove _Signed is it's always true?

Suggested change

	template <class _Ty, class _Uty, bool _Signed, int _Lx>
	template <class _Ty, class _Uty, int _Lx>

[statementreply]

It's reserved for exponential_distribution, where _Signed would be false.

[frederick-vs-ja]

Suggested change

	using _Xtype = conditional_t<_Signed, make_signed_t<_Uty>, _Uty>;
	using _Xtype = make_signed_t<_Uty>;

Also, would it be better to rename it to _Sint_type?

[frederick-vs-ja]

Let's avoid numeric_limits.

Suggested change

	static constexpr _Ty _Width_scale = (_Signed ? _Ty{1} : _Ty{2}) * _Ty{(numeric_limits<_Uty>::max)() / 2u + 1u};
	static constexpr _Ty _Width_scale = _Ty{static_cast<_Uty>(-1) / 2u + 1u};

[MattStephanson]

In common cases (e.g., CVTSI2SD with ties-to-even) this conversion to floating point is slightly biased toward values where the significand ends in 0. It's also correlated with _Layer because a tie only happens when the lower eight or eleven bits of _Xbits are 1 followed by all 0s. Occurs below.

I'm unsure how much effect this has, but if it doesn't affect the performance much, maybe it would be cleaner to mask or shift out the lower bits as they're used.

[Edit: I shouldn't say that round-to-even bias happens only with those value of the lower bits. I think it's sufficient but not necessary.]

[MattStephanson]

The upper bits of _Xbits could be used on the first iteration. Also for the tail, below.

[MattStephanson]

Marsaglia's tail method might be faster if we don't have a fast exponential generator, as it only needs sqrt when returning, instead of an additional log or exp. It can be found in Luc Devroye's book, section IX.1.2, or this DTIC report.