Skip to content

Instantly share code, notes, and snippets.

@andanteyk
Created April 30, 2023 07:14
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save andanteyk/d496afcc71427243e552aadf40680890 to your computer and use it in GitHub Desktop.
Save andanteyk/d496afcc71427243e552aadf40680890 to your computer and use it in GitHub Desktop.
Ziggurat / Modified Ziggurat method
// This code is licensed under the terms of the MIT license.
// For details, see https://gist.github.com/andanteyk/d08ab296665b3fc68df58beff3ea39cb .
// How to build:
// > dotnet new console
// > dotnet add package BenchmarkDotNet
// > dotnet run -c Release
using BenchmarkDotNet.Attributes;
using BenchmarkDotNet.Configs;
using RngLab.Rng.Generators;
using System.Runtime.CompilerServices;
using System.Numerics;
using System.Runtime.InteropServices;
using System.Security.Cryptography;
using BenchmarkDotNet.Running;
using RngLab.Benchmark.Distributions;
namespace RngLab.Rng.Generators
{
/// <summary>
/// 擬似乱数生成器の基底クラス
/// </summary>
public abstract class RandomBase
{
public abstract ulong Next();
}
/// <summary>
/// Seiran128 擬似乱数生成器
/// </summary>
public sealed class Seiran : RandomBase
{
private ulong _state0, _state1;
public Seiran()
{
Span<ulong> buffer = stackalloc ulong[2];
do
{
RandomNumberGenerator.Fill(MemoryMarshal.AsBytes(buffer));
} while (buffer[0] == 0 && buffer[1] == 0);
_state0 = buffer[0];
_state1 = buffer[1];
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public override ulong Next()
{
(ulong s0, ulong s1) = (_state0, _state1);
ulong result = BitOperations.RotateLeft((s0 + s1) * 9, 29) + s0;
_state0 = s0 ^ BitOperations.RotateLeft(s1, 29);
_state1 = s0 ^ s1 << 9;
return result;
}
}
}
namespace RngLab.Benchmark.Distributions
{
[MemoryDiagnoser]
[DisassemblyDiagnoser(printSource: true)]
[GroupBenchmarksBy(BenchmarkLogicalGroupRule.ByCategory)]
[CategoriesColumn]
public class Ziggurat
{
public Seiran Rng = new();
[Benchmark(Baseline = true), BenchmarkCategory("Normal")] public (double, double) BoxMuller() => BoxMuller(Rng);
[Benchmark, BenchmarkCategory("Normal")] public (double, double) PolarMethod() => PolarMethod(Rng);
[Benchmark, BenchmarkCategory("Normal")] public double ZigguratNormal() => ZigguratNormal(Rng);
[Benchmark, BenchmarkCategory("Normal")] public double ModifiedZigguratNormal() => ModifiedZigguratNormal(Rng);
[Benchmark, BenchmarkCategory("Normal")] public double ModifiedZigguratNormalSeparate() => ModifiedZigguratNormalSeparate(Rng);
[Benchmark, BenchmarkCategory("Normal")] public double ModifiedZigguratNormalSimple() => ModifiedZigguratNormalSimple(Rng);
[Benchmark, BenchmarkCategory("Normal")] public double ModifiedZigguratNormal128() => ModifiedZigguratNormal128(Rng);
[Benchmark(Baseline = true), BenchmarkCategory("Exponential")] public double Inverse() => Inverse(Rng);
[Benchmark, BenchmarkCategory("Exponential")] public double ZigguratExponential() => ZigguratExponential(Rng);
[Benchmark, BenchmarkCategory("Exponential")] public double ModifiedZigguratExponential() => ModifiedZigguratExponential(Rng);
/// <summary>
/// すべてのテーブルを生成する
/// 適当にブレークポイントを掛けてデータを集めてください
/// </summary>
public static void CreateTableAll()
{
var nor128 = CreateZigguratTable(128,
x => Math.Exp(-x * x / 2),
x => Math.Sqrt(2 * Math.PI) * 0.5 * (1 + Erf(x / Math.Sqrt(2))),
y => Math.Sqrt(-2 * Math.Log(y))
);
var nor256 = CreateZigguratTable(256,
x => Math.Exp(-x * x / 2),
x => Math.Sqrt(2 * Math.PI) * 0.5 * (1 + Erf(x / Math.Sqrt(2))),
y => Math.Sqrt(-2 * Math.Log(y))
);
var exp128 = CreateZigguratTable(128,
x => Math.Exp(-x),
x => 1 - Math.Exp(-x),
y => -Math.Log(y)
);
var exp256 = CreateZigguratTable(256,
x => Math.Exp(-x),
x => 1 - Math.Exp(-x),
y => -Math.Log(y)
);
var modnor128 = CreateModifiedZigguratTable(128, 128 - 3, true,
x => Math.Exp(-x * x / 2),
x => Math.Sqrt(2 * Math.PI) * 0.5 * (1 + Erf(x / Math.Sqrt(2))),
y => Math.Sqrt(-2 * Math.Log(y)),
x => -x * Math.Exp(-x * x / 2)
);
var modnor256 = CreateModifiedZigguratTable(256, 256 - 3, true,
x => Math.Exp(-x * x / 2),
x => Math.Sqrt(2 * Math.PI) * 0.5 * (1 + Erf(x / Math.Sqrt(2))),
y => Math.Sqrt(-2 * Math.Log(y)),
x => -x * Math.Exp(-x * x / 2)
);
var modexp128 = CreateModifiedZigguratTable(128, 128 - 4, false,
x => Math.Exp(-x),
x => 1 - Math.Exp(-x),
y => -Math.Log(y),
x => -Math.Exp(-x)
);
var modexp256 = CreateModifiedZigguratTable(256, 256 - 4, false,
x => Math.Exp(-x),
x => 1 - Math.Exp(-x),
y => -Math.Log(y),
x => -Math.Exp(-x)
);
}
/// <summary>
/// 誤差関数
/// </summary>
public static double Erf(double real)
{
if (double.IsInfinity(real))
return Math.Sign(real); // ±∞ => ±1
double erf = 2 / Math.Sqrt(Math.PI);
double sum = 0;
double factorial = 1;
for (int i = 0; true; i++)
{
double element = (((i & 1) != 0 ? -1 : 1) * Math.Pow(real, 2 * i + 1)) /
(factorial * (2 * i + 1));
if (sum + element == sum)
break;
sum += element;
factorial *= i + 1;
}
return sum * erf;
}
/// <summary>
/// Ziggurat 法のテーブル生成
/// </summary>
/// <param name="tableSize">128 or 256</param>
/// <param name="probabilityDensityFunction">確率密度関数</param>
/// <param name="cumulativeDistributionFunction">累積分布関数</param>
/// <param name="inverseProbabilityDensityFunction">確率密度関数の逆関数</param>
/// <returns>
/// * 各レイヤー右上角の X 座標
/// * 各レイヤー右上角の Y 座標
/// * 各レイヤーの (確率密度関数と重ならない四角形の面積 / 全体の面積)
/// </returns>
public static (double[] rectangleX, double[] rectangleY, double[] insideRectangleRatio)
CreateZigguratTable(int tableSize,
Func<double, double> probabilityDensityFunction,
Func<double, double> cumulativeDistributionFunction,
Func<double, double> inverseProbabilityDensityFunction)
{
double CalculateArea(double x0)
{
return x0 * probabilityDensityFunction(x0) + cumulativeDistributionFunction(double.PositiveInfinity) - cumulativeDistributionFunction(x0);
}
// 一番下のレイヤーの右上角の x 座標 ref. dn
double x0;
// 1 つのレイヤーの面積 ref. vn
double area;
{
// 誤差関数
double CalculateError(double x0)
{
double area = CalculateArea(x0);
double xi = x0;
for (int i = 1; i < tableSize; i++)
{
xi = inverseProbabilityDensityFunction(area / xi + probabilityDensityFunction(xi));
// 計算失敗時に負になる場合がある
if (xi < 0)
return double.NegativeInfinity;
}
// 面積が足りない場合などに NaN になる場合がある
return double.IsNaN(xi) ? double.NegativeInfinity : xi;
}
// 探索を開始する最小値/最大値 (ヒューリスティック; 計算に失敗した場合は適宜変更する)
double left = 2;
double right = 10;
// 二分法で x0 を求める
while (true)
{
double mid = left * 0.5 + right * 0.5;
if (left == mid && mid < right)
mid = Math.BitIncrement(left);
else if (left < mid && mid == right)
mid = Math.BitDecrement(right);
if (left == mid || right == mid)
break;
double leftError = CalculateError(left);
double midError = CalculateError(mid);
double rightError = CalculateError(right);
if (Math.Sign(leftError) != Math.Sign(midError))
right = mid;
else if (Math.Sign(rightError) != Math.Sign(midError))
left = mid;
else
throw new InvalidOperationException("something wrong. should change left/right");
}
x0 = right;
area = CalculateArea(x0);
}
// 各レイヤー右上角の X 座標 ref. wn
var rectangleX = new double[tableSize];
// 各レイヤー右上角の Y 座標 ref. fn
var rectangleY = new double[tableSize];
// 各レイヤーの (確率密度関数と重ならない四角形の面積 / 全体の面積) ref. kn
var insideRectangleRatio = new double[tableSize];
// rectangleX の設定
// [0]: 一番下のレイヤーを長方形と見立てたときの横幅 (insideRectangleRatio[0] の設定のため)
// [1]: 一番下のレイヤー内の長方形の右端
rectangleX[0] = area / probabilityDensityFunction(x0);
rectangleX[1] = x0;
for (int i = 2; i < rectangleX.Length; i++)
{
double xi = inverseProbabilityDensityFunction(area / rectangleX[i - 1] + probabilityDensityFunction(rectangleX[i - 1]));
rectangleX[i] = xi;
if (xi < 0)
throw new InvalidOperationException("negative x; parameter may be wrong");
}
// rectangleY の設定
// [^1]: 頂点の y 座標
for (int i = 0; i < rectangleY.Length - 1; i++)
{
rectangleY[i] = probabilityDensityFunction(rectangleX[i + 1]);
}
rectangleY[^1] = probabilityDensityFunction(0); // == 1;
// insideRectangleRatio の設定
// [^1]: 一番上のレイヤーはすべて範囲外なので 0
for (int i = 0; i < insideRectangleRatio.Length - 1; i++)
{
insideRectangleRatio[i] = rectangleX[i + 1] / rectangleX[i];
}
insideRectangleRatio[^1] = 0;
return (rectangleX, rectangleY, insideRectangleRatio);
}
/// <summary>
/// 改良 Ziggurat 法のテーブル生成
/// </summary>
/// <param name="tableSize">128 or 256</param>
/// <param name="iMax">正規分布は <paramref name="tableSize"/>-3, 指数分布は <paramref name="tableSize"/>-4 </param>
/// <param name="isSymmetric">正規分布は true, 指数分布は false</param>
/// <param name="probabilityDensityFunction">確率密度関数</param>
/// <param name="cumulativeDistributionFunction">累積分布関数</param>
/// <param name="inversePdf">確率密度関数の逆関数</param>
/// <param name="derivativePdf">確率密度関数の傾き (微分)</param>
/// <returns>
/// * X に係数 2^-n を掛けたもの; n=63(正規分布) or 64(指数分布)
/// * Y に係数 2^-n を掛けたもの
/// * 上位 56 bit が 1 つめの重み, 下位 8 bit が 2 つめのインデックスである Walker's Alias method のテーブル
/// * 変曲点 (確率密度関数のへこみとふくらみの境界) のレイヤーのインデックス
/// * 最大のふくらみの比
/// * 最大のへこみの比
/// </returns>
public static (double[] x, double[] y, ulong[] aliasBox, int inflectionIndex, ulong maxConvexRate, ulong maxConcaveRate)
CreateModifiedZigguratTable(int tableSize, int iMax, bool isSymmetric,
Func<double, double> probabilityDensityFunction,
Func<double, double> cumulativeDistributionFunction,
Func<double, double> inversePdf,
Func<double, double> derivativePdf)
{
// 1 レイヤーの面積
var area = (cumulativeDistributionFunction(double.PositiveInfinity) - cumulativeDistributionFunction(0)) / tableSize;
// 二分法
static double Bisection(Func<double, double> errorFunction, double initialLeft, double initialRight)
{
double left = initialLeft;
double right = initialRight;
while (true)
{
double mid = left * 0.5 + right * 0.5;
if (left == mid && mid < right)
mid = Math.BitIncrement(left);
else if (left < mid && mid == right)
mid = Math.BitDecrement(right);
if (left == mid || right == mid)
break;
double leftError = errorFunction(left);
double midError = errorFunction(mid);
double rightError = errorFunction(right);
if (Math.Sign(leftError) != Math.Sign(midError))
right = mid;
else if (Math.Sign(rightError) != Math.Sign(midError))
left = mid;
else
throw new InvalidOperationException("something wrong. should change left/right");
}
return left;
}
// xtop (x[iMax]) を二分法で求める
double xtop;
{
// 探索の最小/最大値。(ヒューリスティック; 動かないときは適宜変更する)
double xtopLeft = 0.01;
double xtopRight = 0.5 * (isSymmetric ? 1.0 : 0.4);
// 探索の誤差関数。最後のレイヤーの x 座標まで計算可能なら負の、不可能なら正の値を返す
double CalculateError(double xtop)
{
for (int i = iMax - 1; i >= 0; i--)
{
xtop = inversePdf(probabilityDensityFunction(xtop) - area / xtop);
if (double.IsNaN(xtop))
return i + 1;
}
return -xtop;
}
xtop = Bisection(CalculateError, xtopLeft, xtopRight);
}
// x, y の設定
var x = new double[iMax + 1];
var y = new double[iMax + 1];
x[^1] = 0;
y[^1] = probabilityDensityFunction(0);
x[^2] = xtop;
y[^2] = probabilityDensityFunction(xtop);
for (int i = x.Length - 3; i >= 0; i--)
{
x[i] = inversePdf(probabilityDensityFunction(x[i + 1]) - area / x[i + 1]);
y[i] = probabilityDensityFunction(x[i]);
// verification
double actual = x[i + 1] * (y[i + 1] - y[i]);
if (Math.Abs(actual - area) > 1.0 / (1ul << 52) || double.IsNaN(actual))
throw new InvalidOperationException("area mismatch");
}
// 各レイヤー端の面積の計算
var edgeArea = new double[iMax + 1];
edgeArea[0] = cumulativeDistributionFunction(double.PositiveInfinity) - cumulativeDistributionFunction(x[0]);
for (int i = 1; i < edgeArea.Length; i++)
{
edgeArea[i] = (cumulativeDistributionFunction(x[i - 1]) - cumulativeDistributionFunction(x[i]))
- y[i - 1] * (x[i - 1] - x[i]);
}
// verification
{
double actual = edgeArea.Sum();
if (Math.Abs(actual - area * (tableSize - iMax)) > 1.0 / (1ul << 50))
throw new InvalidOperationException("edgeArea mismatch");
}
// Walker's Alias method によるテーブル作成
var firstWeights = new double[tableSize];
var secondIndexes = Enumerable.Range(0, tableSize).ToArray();
double edgeAverage = (tableSize - iMax) * area / tableSize;
{
for (int i = 0; i < edgeArea.Length; i++)
firstWeights[i] = edgeArea[i];
// 小さい順に並べ替える
var sorted = Enumerable.Range(0, tableSize).ToArray();
Array.Sort(sorted, Comparer<int>.Create((a, b) => firstWeights[a].CompareTo(firstWeights[b])));
// 最大の要素から最小の要素に分配していく
int small = 0;
int large = firstWeights.Length - 1;
while (small < large)
{
secondIndexes[sorted[small]] = sorted[large];
firstWeights[sorted[large]] -= edgeAverage - firstWeights[sorted[small]];
switch (firstWeights[sorted[large]].CompareTo(edgeAverage))
{
case > 0:
small++;
break;
case 0:
small++;
firstWeights[sorted[large]] = edgeAverage;
secondIndexes[sorted[large]] = sorted[large];
large--;
break;
case < 0:
(sorted[small], sorted[large]) = (sorted[large], sorted[small]);
large--;
break;
}
}
}
// firstWeights 正規化と検証
for (int i = 0; i < firstWeights.Length; i++)
{
if (firstWeights[i] - edgeAverage > 1.0 / (1ul << 52))
throw new InvalidOperationException("weight mismatch");
firstWeights[i] /= edgeAverage;
}
// 変曲点・ふくらみ/へこみの最大値の探索
int inflectionIndex = tableSize;
double maxConvexRatio = 0;
double maxConcaveRatio = 0;
for (int i = 1; i < x.Length; i++)
{
// > 0 ならふくらみ(convex), < 0 ならへこみ(concave)
double GradientDifference(double t)
{
double pdfGradient = derivativePdf(t);
double lineGradient = (y[i] - y[i - 1]) / (x[i] - x[i - 1]);
return pdfGradient - lineGradient;
}
double maxGradientPoint = Bisection(GradientDifference, x[i], x[i - 1]);
if (inflectionIndex == tableSize && GradientDifference(maxGradientPoint) > 0)
{
inflectionIndex = i;
}
{
double lineGradient = (y[i] - y[i - 1]) / (x[i] - x[i - 1]);
double diff = probabilityDensityFunction(maxGradientPoint) -
(lineGradient * maxGradientPoint + y[i] - lineGradient * x[i]);
double ratio = diff / (y[i] - y[i - 1]);
if (i < inflectionIndex)
maxConcaveRatio = Math.Max(-ratio, maxConcaveRatio);
else
maxConvexRatio = Math.Max(ratio, maxConvexRatio);
}
}
// x, y のスケール
for (int i = 0; i < x.Length; i++)
{
x[i] *= (isSymmetric ? 1.0 : 0.5) / (1ul << 63);
y[i] *= (isSymmetric ? 1.0 : 0.5) / (1ul << 63);
}
// Walker's Alias method のテーブルを一つにまとめる
var aliasBox = new ulong[tableSize];
for (int i = 0; i < aliasBox.Length; i++)
{
ulong ulongWeight = (ulong)(firstWeights[i] * (2.0 * (1ul << 63)));
// verification
if ((ulongWeight & ((ulong)tableSize - 1)) != 0)
throw new InvalidOperationException("weight may lose precision");
aliasBox[i] = ulongWeight ^ (ulong)secondIndexes[i];
}
return (x, y, aliasBox, inflectionIndex, (ulong)(maxConvexRatio * (2.0 * (1ul << 63))), (ulong)(maxConcaveRatio * (2.0 * (1ul << 63))));
}
/// <summary>
/// https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform
/// </summary>
public static (double x, double y) BoxMuller<TRng>(TRng rng)
where TRng : notnull, RandomBase
{
// [0, 1)
double u0 = (rng.Next() >> 11) * (1.0 / (1ul << 53));
double u1 = (rng.Next() >> 11) * (1.0 / (1ul << 53));
double r = Math.Sqrt(-2 * Math.Log(1 - u0));
(double sin, double cos) = Math.SinCos(u1 * (Math.PI * 2));
return (r * cos, r * sin);
}
/// <summary>
/// https://en.wikipedia.org/wiki/Marsaglia_polar_method
/// </summary>
public static (double x, double y) PolarMethod<TRng>(TRng rng)
where TRng : notnull, RandomBase
{
double x, y, s;
do
{
// [-1, 1)
x = ((long)rng.Next() >> 10) * (1.0 / (1ul << 53));
y = ((long)rng.Next() >> 10) * (1.0 / (1ul << 53));
s = x * x + y * y;
} while (s >= 1 || s == 0);
double r = Math.Sqrt(-2 * Math.Log(s) / s);
return (x * r, y * r);
}
private static readonly double[] ZigguratNormalX = new double[] { 3.9107579595237043, 3.6541528853612593, 3.4492782985617532, 3.3202447338401821, 3.2245750520481811, 3.1478892895183974, 3.0835261320025542, 3.027837791770017, 2.9786032798822775, 2.9343668672093322, 2.894121053613866, 2.8571387308736864, 2.8228773968269127, 2.7909211740024045, 2.7609440052804706, 2.7326853590445026, 2.70593365612356, 2.6805146432862492, 2.6562830375772535, 2.6331163936320987, 2.6109105184893453, 2.5895759867088137, 2.5690354526823764, 2.5492215503253211, 2.5300752321603972, 2.5115444416272426, 2.4935830412716, 2.4761499396710813, 2.4592083743352684, 2.4427253182009325, 2.42667098493772, 2.4110184139016977, 2.3957431197825105, 2.3808227951726733, 2.3662370567178832, 2.3519672273797418, 2.33799614879713, 2.3243080188717387, 2.3108882506019826, 2.297723348903479, 2.2848008027251119, 2.2721089902290061, 2.2596370951744165, 2.2473750329480229, 2.2353133849305591, 2.2234433400931528, 2.2117566428848074, 2.2002455466119271, 2.1889027716270157, 2.1777214677409522, 2.1666951803549721, 2.1558178198774054, 2.1450836340485613, 2.1344871828466938, 2.1240233156902049, 2.1136871506873387, 2.1034740557155671, 2.0933796311394861, 2.0833996939990032, 2.0735302635194461, 2.0637675478124398, 2.0541079316513642, 2.0445479652182481, 2.03508435373034, 2.02571394786458, 2.0164337349069341, 2.0072408305612632, 1.9981324713591586, 1.9891060076181815, 1.9801588969012243, 1.9712886979344113, 1.9624930649451195, 1.9537697423854079, 1.945116560009444, 1.936531428276465, 1.9280123340534405, 1.9195573365939675, 1.9111645637720374, 1.9028322085512177, 1.8945585256714979, 1.8863418285375806, 1.8781804862937981, 1.8700729210720737, 1.8620176054004858, 1.8540130597610183, 1.8460578502860066, 1.8381505865836323, 1.8302899196835873, 1.8224745400947211, 1.8147031759671228, 1.8069745913516659, 1.79928758455057, 1.7916409865530174, 1.7840336595503012, 1.7764644955253877, 1.7689324149121384, 1.7614363653197851, 1.7539753203185515, 1.7465482782826074, 1.7391542612868018, 1.7317923140538587, 1.7244615029489456, 1.717160915018729, 1.7098896570722131, 1.7026468548008398, 1.6954316519354835, 1.6882432094381228, 1.6810807047261067, 1.6739433309270635, 1.6668302961626096, 1.6597408228591322, 1.6526741470840112, 1.6456295179057432, 1.6386061967765142, 1.6316034569358457, 1.6246205828340128, 1.6176568695739995, 1.6107116223708198, 1.60378415602709, 1.5968737944237896, 1.5899798700251984, 1.583101723397043, 1.5762387027369265, 1.56939016341615, 1.5625554675320774, 1.5557339834702153, 1.5489250854752188, 1.5421281532300539, 1.5353425714425726, 1.5285677294387776, 1.5218030207620699, 1.5150478427777931, 1.5083015962823969, 1.5015636851165561, 1.4948335157815931, 1.4881104970585544, 1.4813940396293015, 1.474683555698977, 1.4679784586192084, 1.4612781625114117, 1.4545820818895538, 1.4478896312817273, 1.441200224849883, 1.4345132760070591, 1.427828197031431, 1.4211443986764922, 1.4144612897766626, 1.4077782768475986, 1.4010947636804589, 1.3944101509293574, 1.3877238356912009, 1.3810352110770892, 1.3743436657744088, 1.3676485835987275, 1.3609493430345434, 1.3542453167639046, 1.3475358711818661, 1.3408203658976923, 1.3340981532206579, 1.3273685776292334, 1.3206309752223737, 1.3138846731515479, 1.3071289890320688, 1.300363230332185, 1.293586693738306, 1.2867986644946126, 1.2799984157151976, 1.2731852076667469, 1.2663582870196313, 1.2595168860651274, 1.2526602218963219, 1.2457874955500639, 1.238897891107136, 1.231990574747597, 1.2250646937580041, 1.2181193754869677, 1.2111537262451981, 1.2041668301458937, 1.1971577478809672, 1.1901255154282315, 1.1830691426842401, 1.1759876120170198, 1.1688798767324153, 1.1617448594472084, 1.15458145036154, 1.1473885054224768, 1.1401648443697949, 1.1329092486541936, 1.1256204592172099, 1.1182971741210386, 1.1109380460152867, 1.1035416794263686, 1.0961066278537686, 1.0886313906557457, 1.081114409705189, 1.0735540657942411, 1.0659486747639475, 1.058296483332521, 1.050595664592797, 1.0428443131460381, 1.0350404398353523, 1.0271819660375807, 1.0192667174674432, 1.0112924174419795, 1.0032566795466822, 0.99515699963712634, 0.98699074710112478, 0.978755155296315, 0.97044731106634363, 0.96206414322518963, 0.95360240988326583, 0.94505868447037722, 0.93642934028882008, 0.92771053340427945, 0.91889818365190556, 0.90998795349907058, 0.90097522446361256, 0.89185507073537251, 0.88262222958763836, 0.87327106809137733, 0.863795545555871, 0.85418917101077385, 0.84444495491181415, 0.83455535408909487, 0.82451220875506026, 0.8143066701380417, 0.803929116992859, 0.79336905884357589, 0.78261502331025445, 0.77165442422766217, 0.76047340643327954, 0.74905666202106924, 0.73738721143763775, 0.72544614091343629, 0.71321228519451418, 0.70066184111046281, 0.68776789279955486, 0.67449982284118892, 0.66082257424845559, 0.6466957148991842, 0.63207223639042209, 0.61689699001230169, 0.60110461776075474, 0.58461676611138036, 0.56733825705909191, 0.54915170233275123, 0.52990972066750952, 0.509423329608476, 0.487443966146143, 0.463634336798436, 0.43751840221625116, 0.40838913462146925, 0.37512133288940913, 0.33573751922784406, 0.28617459180977611, 0.2152418960131707 };
private static readonly double[] ZigguratNormalY = new double[] { 0.0012602859304974439, 0.002609072746099267, 0.00403797259335825, 0.0055224032992442374, 0.0070508754713644239, 0.0086165827693878132, 0.010214971439688376, 0.011842757857892563, 0.013497450601722269, 0.0151770883079154, 0.016880083152520892, 0.018605121275699969, 0.020351096230017414, 0.022117062707279287, 0.0239022033057638, 0.025705804008514285, 0.027527235669565889, 0.029365939758093519, 0.031221417191877831, 0.033093219458533447, 0.034980941461668344, 0.036884215688516825, 0.038802707404472905, 0.040736110655884963, 0.042684144916415659, 0.04464655225123286, 0.046623094901865926, 0.048613553215801165, 0.050617723860877491, 0.052635418276718977, 0.054666461324812718, 0.056710690106123673, 0.058767952920851532, 0.060838108349454613, 0.062921024437669754, 0.065016577971151415, 0.06712465382769392, 0.069245144396909028, 0.071377949058789442, 0.073522973713877157, 0.075680130358819722, 0.077849336701985419, 0.080030515814549091, 0.082223595813085623, 0.084428509570232818, 0.086645194450434088, 0.088873592068148557, 0.091113648066243016, 0.093365311912556842, 0.095628536712871318, 0.097903279038721328, 0.10018949876866533, 0.10248715894178702, 0.10479622562233526, 0.10711666777452844, 0.10944845714665281, 0.11179156816367555, 0.11414597782767222, 0.11651166562544095, 0.11888861344273637, 0.12127680548461285, 0.12367622820141533, 0.12608687022000081, 0.12850872227981067, 0.1309417771734516, 0.13338602969147251, 0.13584147657105317, 0.13830811644834623, 0.14078594981423628, 0.14327497897330105, 0.14577520800577762, 0.148286642732354, 0.15080929068162108, 0.15334316106003412, 0.15588826472424633, 0.15844461415568725, 0.16101222343726984, 0.16359110823212028, 0.16618128576423236, 0.16878277480095755, 0.17139559563724771, 0.17401977008157618, 0.17665532144346804, 0.1793022745225763, 0.18196065559924679, 0.18463049242651905, 0.18731181422351556, 0.19000465167017574, 0.19270903690329533, 0.1954250035138359, 0.19815258654547216, 0.20089182249434889, 0.20364274931002255, 0.20640540639756363, 0.20917983462080317, 0.21196607630670355, 0.2147641752508421, 0.21757417672399487, 0.22039612747981074, 0.22323007576357126, 0.22607607132202903, 0.2289341654143241, 0.2318044108239774, 0.23468686187196361, 0.23758157443086661, 0.24048860594012389, 0.24340801542236842, 0.24633986350087667, 0.249284212418136, 0.25224112605554438, 0.25521066995425878, 0.25819291133721067, 0.26118791913230727, 0.26419576399684175, 0.26721651834313659, 0.27025025636544514, 0.27329705406814125, 0.2763569892952269, 0.27943014176119091, 0.28251659308325483, 0.28561642681504318, 0.28872972848171846, 0.29185658561662492, 0.29499708779948569, 0.29815132669620342, 0.301319396100315, 0.30450139197615589, 0.3076974125037919, 0.31090755812578014, 0.3141319315958247, 0.31737063802939486, 0.32062378495638039, 0.3238914823758598, 0.32717384281306355, 0.33047098137861913, 0.33378301583016734, 0.33711006663644844, 0.3404522570439576, 0.34380971314627984, 0.347182563956216, 0.35057094148082157, 0.35397498079948536, 0.357394820145182, 0.3608306009890424, 0.36428246812839121, 0.36775056977841253, 0.37123505766761228, 0.3747360871372567, 0.37825381724497742, 0.38178841087274445, 0.38534003483942053, 0.38890886001812447, 0.39249506145864371, 0.39609881851515277, 0.39972031497950977, 0.40335973922041923, 0.40701728432877016, 0.41069314826947695, 0.41438753404017176, 0.41810064983712059, 0.42183270922876004, 0.42558393133727773, 0.42935454102868881, 0.43314476911189109, 0.43695485254721572, 0.44078503466502544, 0.444635565394952, 0.44850670150640665, 0.45239870686104311, 0.45631185267790186, 0.46024641781201886, 0.46420268904734086, 0.4681809614048506, 0.4721815384668776, 0.47620473271864355, 0.48025086590817456, 0.48432026942580109, 0.48841328470456563, 0.49253026364296582, 0.49667156905157667, 0.50083757512522509, 0.505028667942534, 0.50924524599480292, 0.51348772074637083, 0.5177565172287889, 0.52205207467134307, 0.52637484717069427, 0.53072530440266008, 0.53510393237944376, 0.53951123425592618, 0.54394773118898809, 0.54841396325421521, 0.552910490424769, 0.55743789361768981, 0.5619967758134351, 0.56658776325506177, 0.5712115067341369, 0.57586868297122362, 0.58055999609964681, 0.58528617926221294, 0.59004799633165272, 0.59484624376679918, 0.59968175261792189, 0.60455539069624886, 0.60946806492453864, 0.61442072388766289, 0.61941436060456678, 0.624450015545742, 0.62952877992353484, 0.63465179928630411, 0.63982027745171888, 0.645035480819466, 0.65029874310944136, 0.65561147057830238, 0.66097514777524813, 0.66639134390731447, 0.67186171989562515, 0.67738803621729471, 0.6829721616434935, 0.68861608300314736, 0.69432191612456839, 0.70009191813493854, 0.7059285013311557, 0.71183424887662339, 0.7178119326290695, 0.72386453346694934, 0.72999526455976582, 0.73620759812512149, 0.74250529633837792, 0.74889244721735027, 0.75537350650525459, 0.76195334683491711, 0.76863731579656958, 0.77543130497923, 0.78234183265280255, 0.78937614356397945, 0.79654233042086586, 0.8038494831688201, 0.81130787431045737, 0.81892919160144473, 0.82672683394390045, 0.83471629298449379, 0.84291565310973948, 0.85134625845613077, 0.860033621193693, 0.86900868803411613, 0.87830965580606057, 0.88798466075284344, 0.89809592189519782, 0.908726440048799, 0.91999150503578564, 0.93206007595537421, 0.94519895343804106, 0.95987909179524333, 0.97710170126172169, 1 };
private static readonly double[] ZigguratNormalInsideRectangleRatio = new double[] { 0.93438482339784146, 0.9439337670790281, 0.96259114123224732, 0.97118595481323156, 0.97621833534900215, 0.97955355109528319, 0.98194004595758977, 0.9837393825978511, 0.98514860539779781, 0.9862846687490916, 0.98722157019175127, 0.98800851576559712, 0.98867955694412057, 0.98925904142289189, 0.98976486079329329, 0.99021047087166425, 0.99060619510024306, 0.99096009202199753, 0.9912785484011204, 0.99156669443232515, 0.99182870051292471, 0.99206799331942253, 0.99228741575505808, 0.99248934712540759, 0.9926757946571686, 0.992848464053455, 0.99300881450026679, 0.99315810199358789, 0.99329741379121994, 0.99342769604769288, 0.99354977616117823, 0.9936643809806216, 0.99377215174421885, 0.993873656416444, 0.99396939993918676, 0.9940598327987008, 0.994145358223779, 0.99422633826463747, 0.9943030989512045, 0.99437593469007746, 0.99444511202859875, 0.99451087289024254, 0.994573437366301, 0.99463300613354155, 0.99468976255525132, 0.99474387451318824, 0.994795496009965, 0.99484476857486315, 0.99489182250075348, 0.9949367779354179, 0.9949797458469517, 0.99502082887994, 0.99506012211660577, 0.99509771375505118, 0.99513368571497651, 0.99516811417978746, 0.99520107008277459, 0.995232619544, 0.99526282426363755, 0.99529174187675662, 0.99531942627389669, 0.99534592789122311, 0.99537129397358115, 0.99539556881335944, 0.99541879396771271, 0.99544100845639982, 0.995462248942217, 0.99548254989578477, 0.99550194374624068, 0.99552046101921721, 0.99553813046333184, 0.99555497916628022, 0.99557103266150548, 0.99558631502631523, 0.99560084897222323, 0.9956146559282133, 0.9956277561175525, 0.995640168628716, 0.99565191148092913, 0.995663001684783, 0.99567345529833806, 0.99568328747908408, 0.99569251253209423, 0.99570114395468035, 0.99570919447782258, 0.99571667610462522, 0.9957236001460279, 0.99572997725397661, 0.99573581745224582, 0.9957411301650817, 0.99574592424382358, 0.99575020799164649, 0.99575398918655456, 0.99575727510274581, 0.99576007253045506, 0.99576238779437731, 0.99576422677075871, 0.99576559490324223, 0.99576649721754262, 0.99576693833501806, 0.99576692248520682, 0.9957664535173838, 0.99576553491118991, 0.9957641697863856, 0.99576236091176984, 0.99576011071330739, 0.99575742128149891, 0.99575429437802865, 0.99575073144171822, 0.99574673359381649, 0.99574230164264588, 0.99573743608762955, 0.99573213712271824, 0.99572640463923279, 0.99572023822813738, 0.99571363718175621, 0.99570660049494386, 0.99569912686571838, 0.99569121469536226, 0.99568286208799683, 0.99567406684963478, 0.99566482648670873, 0.99565513820407725, 0.99564499890250668, 0.99563440517562163, 0.99562335330632257, 0.99561183926265906, 0.9955998586931516, 0.99558740692155145, 0.99557447894102047, 0.99556106940772537, 0.99554717263381787, 0.99553278257979161, 0.99551789284618941, 0.995502496664638, 0.99548658688818548, 0.99547015598091393, 0.99545319600679316, 0.99543569861774472, 0.99541765504087887, 0.99539905606486445, 0.9953798920253869, 0.99536015278965118, 0.99533982773987573, 0.995318905755724, 0.99529737519561445, 0.99527522387684486, 0.99525243905446115, 0.9952290073987986, 0.99520491497161, 0.99518014720069981, 0.99515468885296454, 0.99512852400574148, 0.99510163601635426, 0.99507400748973462, 0.9950456202439959, 0.99501645527381477, 0.99498649271147244, 0.99495571178538955, 0.9949240907759821, 0.99489160696863876, 0.994858236603617, 0.99482395482262709, 0.99478873561185832, 0.99475255174117716, 0.9947153746992069, 0.99467717462396721, 0.99463792022872566, 0.99459757872268117, 0.99455611572606506, 0.9945134951792004, 0.99446967924502583, 0.99442462820453559, 0.99437830034453645, 0.9943306518370626, 0.99428163660972479, 0.9942312062061972, 0.99417930963596024, 0.994125893212333, 0.99407090037772, 0.99401427151488952, 0.99395594374296725, 0.993895850696692, 0.99383392228731093, 0.9937700844433186, 0.99370425882903357, 0.99363636253877685, 0.99356630776415478, 0.99349400143165123, 0.99341934480739791, 0.99334223306560854, 0.99326255481671988, 0.99318019159079463, 0.99309501727115479, 0.9930068974725722, 0.99291568885758263, 0.99282123838362, 0.99272338247267478, 0.99262194609401588, 0.9925167417491656, 0.99240756834677712, 0.99229420995320716, 0.99217643440249148, 0.99205399174689413, 0.99192661252630221, 0.99179400583126154, 0.99165585713037485, 0.99151182582792508, 0.99136154251183151, 0.99120460584513936, 0.99104057904600462, 0.990868985891185, 0.99068930616606743, 0.990500970469704, 0.990303354265623, 0.990095771047515, 0.98987746446227787, 0.98964759920003453, 0.98940525041995608, 0.9891493914298386, 0.98887887927354934, 0.98859243779989714, 0.98828863768420283, 0.987965872743101, 0.98762233171485958, 0.987255964459405, 0.98686444124727868, 0.98644510343143876, 0.985994903297113, 0.98551033021605039, 0.98498731932552863, 0.98442113771213846, 0.98380624136276962, 0.98313609373740962, 0.98240293339783336, 0.98159747319695245, 0.9807085063183768, 0.97972238371372611, 0.97862231119247456, 0.97738738919734169, 0.97599127836276789, 0.97440030911606157, 0.97257074532084742, 0.97044472541008764, 0.96794407128365678, 0.96496053534296011, 0.96133984665684658, 0.95685442306062907, 0.95115412026544921, 0.94367126739894869, 0.93342161736002993, 0.91853896464974361, 0.89501046672497309, 0.85237596461647558, 0.75213489308039172, 0 };
public static double ZigguratNormal<TRng>(TRng rng)
where TRng : notnull, RandomBase
{
const int TableSize = 256;
while (true)
{
long u1 = (long)rng.Next();
// [-1, 1)
double u1Double = (u1 >> 10) * (1.0 / (1ul << 53));
int rectangleIndex = (int)(u1 & (TableSize - 1));
// 確率密度関数の範囲内にある四角形なら、そのまま返す
if (Math.Abs(u1Double) < ZigguratNormalInsideRectangleRatio[rectangleIndex])
return u1Double * ZigguratNormalX[rectangleIndex];
if (rectangleIndex == 0)
{
// 一番下のレイヤーの尾の部分
double s, t;
double x0 = ZigguratNormalX[1];
do
{
s = -Math.Log(1 - (rng.Next() >> 11) * (1.0 / (1ul << 53))) / x0;
t = -Math.Log(1 - (rng.Next() >> 11) * (1.0 / (1ul << 53)));
} while (s * s > t + t);
return (x0 + s) * (u1 < 0 ? -1 : 1);
}
else
{
// それ以外の端の部分
double x = u1Double * ZigguratNormalX[rectangleIndex];
double y = (rng.Next() >> 11) * (1.0 / (1ul << 53));
if (ZigguratNormalY[rectangleIndex - 1] + y * (ZigguratNormalY[rectangleIndex] - ZigguratNormalY[rectangleIndex - 1]) <= Math.Exp(-0.5 * x * x))
return x;
}
}
}
private static readonly double[] ModifiedZigguratNormalX = { 3.942166282539769E-19, 3.7204945004118776E-19, 3.5827024480628514E-19, 3.4807476236540124E-19, 3.3990177171882035E-19, 3.3303778360340052E-19, 3.2709438817617477E-19, 3.2183577132495033E-19, 3.1710758541840374E-19, 3.1280307407034012E-19, 3.088452065580397E-19, 3.0517650624107304E-19, 3.0175290292584557E-19, 2.9853983440705282E-19, 2.9550967462801758E-19, 2.9263997988491624E-19, 2.8991225869977438E-19, 2.8731108780226257E-19, 2.84823463271013E-19, 2.8243831535194356E-19, 2.8014613964727E-19, 2.7793871261807768E-19, 2.7580886921411178E-19, 2.7375032698308729E-19, 2.7175754543391022E-19, 2.698256124753846E-19, 2.6795015188771481E-19, 2.6612724730440009E-19, 2.6435337927976614E-19, 2.6262537282028419E-19, 2.6094035335224123E-19, 2.5929570954330983E-19, 2.5768906173214707E-19, 2.5611823497719588E-19, 2.5458123593393342E-19, 2.5307623292372444E-19, 2.5160153867798386E-19, 2.5015559533646177E-19, 2.4873696135403144E-19, 2.4734430003079192E-19, 2.4597636942892717E-19, 2.446320134791244E-19, 2.4331015411139196E-19, 2.4200978427132945E-19, 2.4072996170445874E-19, 2.3946980340903342E-19, 2.3822848067252665E-19, 2.3700521461931791E-19, 2.3579927220741325E-19, 2.3460996262069967E-19, 2.334366340105445E-19, 2.3227867054673835E-19, 2.311354897430376E-19, 2.3000654002704233E-19, 2.28891298527976E-19, 2.2778926905921892E-19, 2.2669998027527317E-19, 2.2562298398527411E-19, 2.2455785360727256E-19, 2.2350418274933907E-19, 2.2246158390513289E-19, 2.2142968725296249E-19, 2.2040813954857555E-19, 2.1939660310297606E-19, 2.1839475483749618E-19, 2.1740228540916853E-19, 2.1641889840016519E-19, 2.1544430956570613E-19, 2.1447824613540345E-19, 2.1352044616350571E-19, 2.1257065792395107E-19, 2.1162863934653125E-19, 2.1069415749082023E-19, 2.0976698805483467E-19, 2.0884691491567363E-19, 2.0793372969963634E-19, 2.0702723137954107E-19, 2.0612722589717129E-19, 2.0523352580895635E-19, 2.0434594995315797E-19, 2.0346432313698151E-19, 2.0258847584216418E-19, 2.0171824394771313E-19, 2.0085346846857531E-19, 1.9999399530912018E-19, 1.9913967503040587E-19, 1.9829036263028147E-19, 1.9744591733545175E-19, 1.9660620240469855E-19, 1.9577108494251483E-19, 1.9494043572246305E-19, 1.9411412901962158E-19, 1.9329204245152932E-19, 1.9247405682708166E-19, 1.9166005600287069E-19, 1.9084992674649821E-19, 1.9004355860642335E-19, 1.892408437879372E-19, 1.8844167703488431E-19, 1.8764595551677744E-19, 1.8685357872097447E-19, 1.8606444834960931E-19, 1.852784682209879E-19, 1.8449554417517925E-19, 1.8371558398354866E-19, 1.8293849726199562E-19, 1.8216419538767388E-19, 1.8139259141898443E-19, 1.8062360001864449E-19, 1.7985713737964738E-19, 1.790931211539384E-19, 1.7833147038364195E-19, 1.7757210543468423E-19, 1.768149479326639E-19, 1.7605992070083135E-19, 1.7530694770004407E-19, 1.7455595397057215E-19, 1.7380686557563472E-19, 1.7305960954655261E-19, 1.7231411382940904E-19, 1.7157030723311375E-19, 1.7082811937877138E-19, 1.7008748065025788E-19, 1.6934832214591352E-19, 1.6861057563126349E-19, 1.6787417349268046E-19, 1.6713904869190636E-19, 1.6640513472135291E-19, 1.656723655601024E-19, 1.6494067563053264E-19, 1.6420999975549112E-19, 1.6348027311594529E-19, 1.6275143120903658E-19, 1.6202340980646725E-19, 1.6129614491314931E-19, 1.6056957272604589E-19, 1.5984362959313479E-19, 1.5911825197242491E-19, 1.5839337639095552E-19, 1.57668939403708E-19, 1.5694487755235887E-19, 1.5622112732380259E-19, 1.5549762510837067E-19, 1.5477430715767269E-19, 1.5405110954198328E-19, 1.5332796810709686E-19, 1.5260481843056971E-19, 1.5188159577726681E-19, 1.5115823505412761E-19, 1.5043467076406196E-19, 1.4971083695888395E-19, 1.4898666719118714E-19, 1.4826209446506113E-19, 1.4753705118554368E-19, 1.468114691066983E-19, 1.4608527927820112E-19, 1.4535841199031451E-19, 1.4463079671711862E-19, 1.4390236205786415E-19, 1.4317303567630175E-19, 1.4244274423783481E-19, 1.4171141334433217E-19, 1.4097896746642792E-19, 1.4024532987312287E-19, 1.3951042255849034E-19, 1.3877416616527576E-19, 1.3803647990516385E-19, 1.3729728147547174E-19, 1.3655648697200824E-19, 1.3581401079782068E-19, 1.35069765567529E-19, 1.3432366200692418E-19, 1.3357560884748263E-19, 1.3282551271542044E-19, 1.3207327801488085E-19, 1.3131880680481522E-19, 1.3056199866908071E-19, 1.2980275057923783E-19, 1.2904095674948603E-19, 1.2827650848312722E-19, 1.2750929400989208E-19, 1.2673919831340477E-19, 1.2596610294799505E-19, 1.2518988584399369E-19, 1.2441042110056516E-19, 1.2362757876504158E-19, 1.2284122459762067E-19, 1.2205121982017847E-19, 1.2125742084782237E-19, 1.2045967900166969E-19, 1.1965784020118015E-19, 1.188517446341955E-19, 1.1804122640264086E-19, 1.1722611314162059E-19, 1.1640622560939106E-19, 1.1558137724540872E-19, 1.1475137369333182E-19, 1.1391601228549045E-19, 1.1307508148492589E-19, 1.1222836028063025E-19, 1.1137561753107903E-19, 1.1051661125053528E-19, 1.0965108783189755E-19, 1.0877878119905372E-19, 1.0789941188076654E-19, 1.0701268599703639E-19, 1.0611829414763283E-19, 1.0521591019102927E-19, 1.043051899002755E-19, 1.033857694803547E-19, 1.0245726392923698E-19, 1.0151926522209309E-19, 1.0057134029488234E-19, 9.96130287996728E-20, 9.864384059945989E-20, 9.76632529647558E-20, 9.6670707427623454E-20, 9.5665606240866658E-20, 9.46473083804332E-20, 9.3615125017323484E-20, 9.256831437088727E-20, 9.1506075837638762E-20, 9.04275432677257E-20, 8.9331777233763668E-20, 8.8217756102327859E-20, 8.7084365674892307E-20, 8.5930387109612138E-20, 8.4754482764244337E-20, 8.3555179508462331E-20, 8.2330848933585352E-20, 8.1079683729129829E-20, 7.9799669284133852E-20, 7.8488549286072733E-20, 7.714378370093468E-20, 7.5762496979467554E-20, 7.4341413578485329E-20, 7.2876776807378419E-20, 7.1364245443525362E-20, 6.9798760240761054E-20, 6.8174368944799042E-20, 6.6483992986198527E-20, 6.4719110345162767E-20, 6.2869314813103686E-20, 6.0921687548281251E-20, 5.8859873575576818E-20, 5.6662675116090981E-20, 5.4301813630894559E-20, 5.17381717444942E-20, 4.8915031722398539E-20, 4.5744741890755295E-20, 4.2078802568583392E-20, 3.7625986722404749E-20, 3.16285898058819E-20, 0 };
private static readonly double[] ModifiedZigguratNormalY = { 1.4598410796621224E-22, 3.0066613427945036E-22, 4.6129728815105761E-22, 6.2663350049236694E-22, 7.9594524761883951E-22, 9.6874655021707428E-22, 1.1446877002379678E-21, 1.3235036304379412E-21, 1.5049857692053373E-21, 1.688965300071954E-21, 1.8753025382711875E-21, 2.0638798423695436E-21, 2.2545966913644952E-21, 2.4473661518802054E-21, 2.6421122727763804E-21, 2.8387681187880171E-21, 3.0372742567457558E-21, 3.2375775699986863E-21, 3.4396303157949074E-21, 3.6433893657998091E-21, 3.8488155868912591E-21, 4.0558733309493069E-21, 4.2645300104283891E-21, 4.4747557422305367E-21, 4.686523046535586E-21, 4.8998065902775521E-21, 5.1145829672105745E-21, 5.3308305082046459E-21, 5.5485291167032029E-21, 5.767660125269071E-21, 5.9882061699178717E-21, 6.2101510795442477E-21, 6.4334797782257449E-21, 6.6581781985714183E-21, 6.8842332045893437E-21, 7.1116325227957336E-21, 7.3403646804903347E-21, 7.5704189502886659E-21, 7.8017853001379984E-21, 8.0344543481570242E-21, 8.2684173217333329E-21, 8.5036660203915217E-21, 8.7401927820109687E-21, 8.9779904520282081E-21, 9.217052355306156E-21, 9.4573722703928955E-21, 9.69894440592696E-21, 9.9417633789758589E-21, 1.0185824195119829E-20, 1.0431122230114781E-20, 1.0677653212987408E-20, 1.0925413210432012E-20, 1.1174398612392903E-20, 1.1424606118728722E-20, 1.1676032726866311E-20, 1.1928675720361041E-20, 1.2182532658289385E-20, 1.2437601365406799E-20, 1.2693879923010687E-20, 1.2951366660454158E-20, 1.3210060147261467E-20, 1.3469959185800735E-20, 1.3731062804473641E-20, 1.3993370251385593E-20, 1.4256880988463133E-20, 1.4521594685988372E-20, 1.4787511217522905E-20, 1.505463065519617E-20, 1.5322953265335221E-20, 1.5592479504415048E-20, 1.5863210015310325E-20, 1.6135145623830982E-20, 1.6408287335525598E-20, 1.6682636332737935E-20, 1.6958193971903124E-20, 1.7234961781071113E-20, 1.7512941457646084E-20, 1.779213486633149E-20, 1.807254403727107E-20, 1.8354171164377274E-20, 1.8637018603838942E-20, 1.8921088872801E-20, 1.9206384648209465E-20, 1.9492908765815639E-20, 1.9780664219333854E-20, 2.006965415974783E-20, 2.0359881894760853E-20, 2.0651350888385693E-20, 2.0944064760670542E-20, 2.1238027287557472E-20, 2.1533242400870493E-20, 2.1829714188430486E-20, 2.2127446894294606E-20, 2.2426444919118282E-20, 2.2726712820637813E-20, 2.3028255314272291E-20, 2.3331077273843579E-20, 2.3635183732413304E-20, 2.394057988323637E-20, 2.4247271080830292E-20, 2.4555262842160345E-20, 2.486456084794038E-20, 2.5175170944049631E-20, 2.5487099143065938E-20, 2.5800351625916006E-20, 2.61149347436437E-20, 2.6430855019297335E-20, 2.6748119149937423E-20, 2.7066734008766265E-20, 2.7386706647381211E-20, 2.7708044298153576E-20, 2.8030754376735287E-20, 2.8354844484695771E-20, 2.8680322412291655E-20, 2.9007196141372144E-20, 2.9335473848423237E-20, 2.9665163907754E-20, 2.999627489482863E-20, 3.0328815589748068E-20, 3.0662794980885293E-20, 3.0998222268678766E-20, 3.1335106869588609E-20, 3.1673458420220558E-20, 3.2013286781622988E-20, 3.2354602043762618E-20, 3.2697414530184806E-20, 3.304173480286495E-20, 3.3387573667257349E-20, 3.3734942177548944E-20, 3.4083851642125214E-20, 3.4434313629256261E-20, 3.4786339973011394E-20, 3.5139942779411176E-20, 3.5495134432826177E-20, 3.585192760263246E-20, 3.6210335250134172E-20, 3.657037063576439E-20, 3.6932047326575888E-20, 3.7295379204034264E-20, 3.7660380472126407E-20, 3.802706566579829E-20, 3.8395449659736661E-20, 3.8765547677510173E-20, 3.9137375301086418E-20, 3.9510948480742178E-20, 3.9886283545385442E-20, 4.0263397213308578E-20, 4.0642306603393553E-20, 4.1023029246790973E-20, 4.140558309909645E-20, 4.1789986553048823E-20, 4.2176258451776819E-20, 4.2564418102621753E-20, 4.2954485291566191E-20, 4.3346480298300112E-20, 4.3740423911958146E-20, 4.4136337447563716E-20, 4.4534242763218286E-20, 4.493416227807625E-20, 4.5336118991149031E-20, 4.5740136500984466E-20, 4.6146239026271279E-20, 4.6554451427421133E-20, 4.6964799229185082E-20, 4.7377308644364926E-20, 4.7792006598684163E-20, 4.8208920756888113E-20, 4.8628079550147814E-20, 4.9049512204847647E-20, 4.9473248772842596E-20, 4.9899320163277668E-20, 5.0327758176068971E-20, 5.0758595537153414E-20, 5.11918659356227E-20, 5.1627604062866077E-20, 5.2065845653856434E-20, 5.2506627530725224E-20, 5.2949987648783478E-20, 5.339596514515945E-20, 5.3844600390237606E-20, 5.42959350420994E-20, 5.47500121041839E-20, 5.5206875986405109E-20, 5.5666572569983857E-20, 5.6129149276275828E-20, 5.6594655139902512E-20, 5.70631408865206E-20, 5.7534659015596954E-20, 5.8009263888591254E-20, 5.8487011822987619E-20, 5.8967961192659828E-20, 5.9452172535103507E-20, 5.9939708666122642E-20, 6.0430634802618953E-20, 6.0925018694200555E-20, 6.1422930764402872E-20, 6.1924444262401555E-20, 6.2429635426193963E-20, 6.2938583658336226E-20, 6.3451371715447575E-20, 6.3968085912834963E-20, 6.4488816345752724E-20, 6.5013657128995346E-20, 6.5542706656731714E-20, 6.6076067884730729E-20, 6.6613848637404208E-20, 6.7156161942412992E-20, 6.770312639595058E-20, 6.825486656224642E-20, 6.8811513411327837E-20, 6.9373204799659693E-20, 6.9940085998959121E-20, 7.0512310279279515E-20, 7.1090039553397179E-20, 7.1673445090644808E-20, 7.22627083096558E-20, 7.2858021661057338E-20, 7.3459589613035812E-20, 7.4067629754967565E-20, 7.4682374037052829E-20, 7.5304070167226678E-20, 7.5932983190698559E-20, 7.6569397282483766E-20, 7.721361778948769E-20, 7.7865973566417028E-20, 7.8526819659456767E-20, 7.9196540403850572E-20, 7.987555301703798E-20, 8.0564311788901642E-20, 8.1263312996426188E-20, 8.1973100703706316E-20, 8.2694273652634046E-20, 8.3427493508836792E-20, 8.4173494807453416E-20, 8.4933097052832066E-20, 8.57072195782309E-20, 8.64968999859307E-20, 8.7303317295655339E-20, 8.8127821378859516E-20, 8.8971970928196666E-20, 8.9837583239314076E-20, 9.0726800697869543E-20, 9.1642181484063544E-20, 9.2586826406702765E-20, 9.3564561480278864E-20, 9.4580210012636175E-20, 9.564001555085037E-20, 9.6752334770503142E-20, 9.7928851697808831E-20, 9.9186905857531331E-20, 1.0055456271343398E-19, 1.0208407377305566E-19, 1.0390360993240711E-19, 1.0842021724855044E-19 };
private static readonly ulong[] ModifiedZigguratNormalAliasBox = { 0xc6d1bcb1d7c5c8fb, 0xc58e65fed90c18f7, 0xbdfe727f7b5100f4, 0xdf559b5c3fd230f1, 0xfd37a81820b170ef, 0xda845732a9b27003, 0xc18bb1c5c69448f3, 0xaea2e0fe0c50f0f5, 0x9fc6849420a6c8f6, 0x93c2f42cbec160f7, 0x89d49e0c64c588f8, 0x8178f3a67d05e0f9, 0x7a554ee0b04510f9, 0x7428aefe15e7f4fa, 0x6ec330fc6a9140fb, 0x6a00bca9d20520fb, 0x65c592cf1a7894fb, 0x61fc00beb1018800, 0x5e92cd0eec7fc800, 0x5b7c1cb2569b60fc, 0x58aca85e4418b0fc, 0x561b28c74e477cfc, 0x53bfe910a63250fc, 0x5194745ca09590fd, 0x4f9356ee453d70fd, 0x4db7eda0d49b90fd, 0x4bfe400809666cfd, 0x4a62e28fed3b58fd, 0x48e2deb6dce1f0fd, 0x477b9ff6890330fd, 0x462ae44f1e19e8fd, 0x44eeafaec37ed8fd, 0x43c54198bbb850fd, 0x42ad0c96e746f4fd, 0x41a4af1f120398fd, 0x40aaed9f2142d0fd, 0x3fbead7e1c2650fd, 0x3edef0e17972b0fd, 0x3e0ad317806c02fd, 0x3d418587825728fd, 0x3c824d1301ff92fd, 0x3bcc7fd4ba1b34fd, 0x3b1f832e92fdb8fd, 0x3a7aca1973b674fd, 0x39ddd3aedfad24fd, 0x394829e0df9906fd, 0x38b9605bcc8b1afd, 0x3831138b23cd28fd, 0x37aee7bbd833a6fd, 0x37328858ffdef4fd, 0x36bba73f647326fd, 0x3649fc239c9b0cfd, 0x35dd440a3ca3e4fd, 0x357540cd98fcdcfd, 0x3511b8b1b6f68cfd, 0x34b27603192b44fd, 0x345746bed152d2fd, 0x33fffc451e8a3afd, 0x33ac6b11ce6374fd, 0x335c6a7d782972fd, 0x330fd4834a7bf0fd, 0x32c6858cd99e9efd, 0x32805c436b10a6fd, 0x323d39647de4dafd, 0x31fcff9b1cbc2efd, 0x31bf935c0afec2fd, 0x3184dac54d663cfd, 0x314cbd80aea4b6fd, 0x311724a8b10300fd, 0x30e3faaf5c30f4fd, 0x30b32b47f4eed4fd, 0x3084a3519d16c8fd, 0x305850c4757e66fd, 0x302e229fe376a2fd, 0x300608da4510a6fd, 0x2fdff451a1855afd, 0x2fbbd6be63f1eefd, 0x2f99a2a64bd204fd, 0x2f794b507d5784fd, 0x2f5ac4babd80a4fd, 0x2f3e038f49a680fd, 0x2f22fd1b9810c4fd, 0x2f09a74772a6d0fd, 0x2ef1f88d34f902fd, 0x2edbe7f1c515c6fd, 0x2ec76cfe738d88fd, 0x2eb47fba0c19dafd, 0x2ea318a28e748efd, 0x2e9330a894d4cafd, 0x2e84c1291fc1cafd, 0x2e77c3e97a4dcefd, 0x2e6c3312903184fd, 0x2e62092c7d15f2fd, 0x2e59411b71ded4fd, 0x2e51d61b49fe3efd, 0x2e4bc3bccd1166fd, 0x2e4705e25f1ef2fd, 0x2e4398bd39fd94fd, 0x2e4178cab1778afd, 0x2e40a2d1fb9f2cfd, 0x2e4113e198c4c8fd, 0x2e42c94de1bd5afd, 0x2e45c0ae0c8020fd, 0x2e49f7dbb2dc20fd, 0x2e4f6cf035308efd, 0x2e561e43fc4bbefd, 0x2e5e0a6cbd8db0fd, 0x2e67303c4fe530fd, 0x2e718ebfe48e86fd, 0x2e7d253f029ffafd, 0x2e89f33a94dd24fd, 0x2e97f86c0ae26afd, 0x2ea734c573754cfd, 0x2eb7a8703ea7dafd, 0x2ec953cd280342fd, 0x2edc377436e01cfd, 0x2ef0543416c832fd, 0x2f05ab123d28a2fd, 0x2f1c3d4b2da8d4fd, 0x2f340c5207034afd, 0x2f4d19d1325b6efd, 0x2f6767aadb615afd, 0x2f82f7f87920a0fd, 0x2f9fcd0c79737afd, 0x2fbde971df0bb2fd, 0x2fdd4fedcdd7bafd, 0x2ffe037f772576fd, 0x30200761cc76dafd, 0x30435f0c42b708fd, 0x30680e33921a30fd, 0x308e18cbac79aafd, 0x30b5830850abb8fd, 0x30de515f4a0566fd, 0x310888894f49a6fd, 0x31342d84491134fd, 0x31614594c1269afd, 0x318fd6486c0d4cfd, 0x31bfe57779f8b4fd, 0x31f17947e4e340fd, 0x3224982ee3906cfd, 0x325948f43c8910fd, 0x328f92b4cc68f2fd, 0x32c77ce56f2f8efd, 0x33010f562b3bd6fd, 0x333c5235b88b54fd, 0x33794e1501d94afd, 0x33b80beac83996fd, 0x33f89517e43358fd, 0x343af36b92a066fd, 0x347f3127c9566cfd, 0x34c55906684954fd, 0x350d763e286330fd, 0x3557948851df6efd, 0x35a3c02639f2e2fd, 0x35f205e7cfddc8fd, 0x3642733226acc4fd, 0x3695160663fbc8fd, 0x36e9fd09463002fd, 0x3741378b39d708fd, 0x379ad590b4f468fd, 0x37f6e7db2e059cfd, 0x38557ff35d637cfd, 0x38b6b03277c0eefd, 0x391a8bce450fdefd, 0x398126e4252d20fd, 0x39ea9685db8464fd, 0x3a56f0c70e33d8fd, 0x3ac64ccb521042fd, 0x3b38c2d5fabeb2fd, 0x3bae6c5a0bd7defd, 0x3c27640c403414fd, 0x3ca3c5f5b1b4f2fd, 0x3d23af886fd6bcfd, 0x3da73fb5118a5afd, 0x3e2e97023a62d6fd, 0x3eb9d7a61501d0fd, 0x3f4925a1214e68fd, 0x3fdca6dc15b5f4fd, 0x40748346fbe7a8fd, 0x4110e4fb856bb8fd, 0x41b1f8627a85d0fd, 0x4257ec5b004f08fd, 0x4302f26661d574fd, 0x43b33ed7116e40fd, 0x4469090393bba4fd, 0x45248b7deb8264fd, 0x45e6045012d9e0fd, 0x46adb53da6aed4fd, 0x477be40bd0fa20fd, 0x4850dacf591c1cfd, 0x492ce842a0fd70fd, 0x4a1060233adb10fd, 0x4afb9b98ffb8dcfd, 0x4beef9a70adef0fd, 0x4ceadfa81a5098fd, 0x4defb9d72c7494fd, 0x4efdfbe72fe4d0fd, 0x501621a9608c34fd, 0x5138afc6a7a9f0fd, 0x5266348c2347b8fc, 0x539f48cf4aad78fc, 0x54e490eb262980fc, 0x5636bddb69fe9cfc, 0x57968e77cf0548fc, 0x5904d0d6975054fc, 0x5a8263d9629960fc, 0x5c1038ebfe230400, 0x5daf55fbfd149000, 0x5f60d7b2ef00d400, 0x6125f3fa86025800, 0x62fffcdb2c4b1800, 0x64f063bf243f3800, 0x66f8bd2da9e8d8fb, 0x691ac5123b4ba8fb, 0x6b5863a76bef24fb, 0x6db3b324abb4d8fb, 0x702f0650b62f28fa, 0x72ccf025b46f0cfa, 0x75904cbbd8328cfa, 0x787c4bbbb09d30fa, 0x7b947ca99e2e4cf9, 0x7edcdd6e5388c4f9, 0x8259eb9b2bd398f8, 0x8610b907d7b2c0f8, 0x8a070493b6c2f001, 0x8e435809cf16f001, 0x92cd2c7476d78801, 0x97ad168ac8e140f7, 0x9cecfd6d8baee0f7, 0xa2985e97c21140f6, 0xa8bca2e6cb4158f5, 0xaf6989f6a2190002, 0xb6b1b2fd462b6002, 0xbeab4d13c9c930f4, 0xc770fcdd4cdad0f3, 0xd1230b71f7ef48f2, 0xdbe8fb696740b003, 0xe7f3aeac31d898f1, 0xf5805d6770e7f0f0, 0xfffffffffd71c0ef, 0xfdeb966ea72228ef, 0xf1fe2559fa5188f0, 0xf7bdcf2ab32e5003, 0xe4368349dbba40f2, 0xfcccff52b1bd70f3, 0xcf0a381011797802, 0xce3286a74622e0f5, 0xae75adb6e22950f6, 0xd722c4a089ebd801, 0xc9dc5f7b7d2d80f8, 0xd7c9bdcbf93010f9, 0xcfef53b6a4d3c8fa, 0xb786b6325d47c000, 0xb06b5777e76e98fc, 0x00000000000000fd, 0x00000000000000fd };
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static double ModifiedZigguratNormal<TRng>(TRng rng)
where TRng : notnull, RandomBase
{
const int TableSize = 256;
// 四角形のレイヤー数
const int IMax = 253;
long u1 = (long)rng.Next();
int tableX = (int)(u1 & (TableSize - 1));
// 四角形のレイヤーなら x 座標を返す
if (tableX < IMax)
{
return u1 * ModifiedZigguratNormalX[tableX];
}
// 以降の処理は相対的に重いため別関数に分け、上の高速パスのインライン化を試みる
[MethodImpl(MethodImplOptions.NoInlining)]
static double Fallback(TRng rng, long u1)
{
const int InflectionIndex = 204;
const ulong MaxConvexRate = 0x3efb83be6450cc00;
const ulong MaxConcaveRate = 0x151b6b6b7cd81f00;
[MethodImpl(MethodImplOptions.AggressiveInlining)]
static double SampleX(double x, int tableY) => Math.FusedMultiplyAdd(x, ModifiedZigguratNormalX[tableY - 1] - ModifiedZigguratNormalX[tableY], ModifiedZigguratNormalX[tableY] * (1ul << 63));
[MethodImpl(MethodImplOptions.AggressiveInlining)]
static double SampleY(double y, int tableY) => Math.FusedMultiplyAdd(y, ModifiedZigguratNormalY[tableY - 1] - ModifiedZigguratNormalY[tableY], ModifiedZigguratNormalY[tableY] * (1ul << 63));
// Walker's Alias method のテーブルを引く
ulong u2 = rng.Next();
int tableAlias = (int)(u2 & (TableSize - 1));
int tableY = (u2 >> 8) < (ModifiedZigguratNormalAliasBox[tableAlias] >> 8) ? tableAlias : (int)(ModifiedZigguratNormalAliasBox[tableAlias] & (TableSize - 1));
if (tableY == 0)
{
// 最下層レイヤー: 尾から生成する
double x0 = ModifiedZigguratNormalX[0] * (1ul << 63);
double s, t;
do
{
s = ModifiedZigguratExponential(rng) / x0;
t = ModifiedZigguratExponential(rng);
} while (s * s > t + t);
return (x0 + s) * (u1 < 0 ? -1 : 1);
}
else if (tableY < InflectionIndex)
{
// へこみレイヤー
ulong ux = (ulong)u1 << 1 >> 1;
ulong uy = rng.Next() >> 1;
double tx;
for (; ; ux = rng.Next() >> 1, uy = rng.Next() >> 1)
{
// 三角形の右上; 上下左右反転させる
if (uy < ux)
{
ux = (1ul << 63) - 1 - ux;
uy = (1ul << 63) - 1 - uy;
}
// 三角形の左下; 採択
if (uy >= ux + MaxConcaveRate)
{
tx = SampleX(ux, tableY);
break;
}
// 確率密度関数より下なら採択
tx = SampleX(ux, tableY);
double ty = SampleY(uy, tableY);
if (ty <= Math.Exp(-0.5 * tx * tx))
{
break;
}
}
return tx * (u1 < 0 ? -1 : 1);
}
else if (tableY == InflectionIndex)
{
// へこみ+ふくらみレイヤー
ulong ux = (ulong)u1 << 1 >> 1;
ulong uy = rng.Next() >> 1;
double tx;
for (; ; ux = rng.Next() >> 1, uy = rng.Next() >> 1)
{
// 三角形の左下; 採択
if (uy >= ux + MaxConcaveRate)
{
tx = SampleX(ux, tableY);
break;
}
// 三角形の右上; 棄却
if (uy + MaxConvexRate <= ux)
{
continue;
}
// 確率密度関数より下なら採択
tx = SampleX(ux, tableY);
double ty = SampleY(uy, tableY);
if (ty <= Math.Exp(-0.5 * tx * tx))
{
break;
}
}
return tx * (u1 < 0 ? -1 : 1);
}
else
{
// ふくらみレイヤー
ulong ux = (ulong)u1 << 1 >> 1;
ulong uy = rng.Next() >> 1;
double tx;
for (; ; ux = rng.Next() >> 1, uy = rng.Next() >> 1)
{
// 三角形の左下; 採択
if (uy >= ux)
{
tx = SampleX(ux, tableY);
break;
}
// 三角形の右上; 棄却
if (uy + MaxConvexRate <= ux)
{
continue;
}
// 確率密度関数より下なら採択
tx = SampleX(ux, tableY);
double ty = SampleY(uy, tableY);
if (ty <= Math.Exp(-tx * tx / 2))
{
break;
}
}
return tx * (u1 < 0 ? -1 : 1);
}
}
return Fallback(rng, u1);
}
private static readonly ulong[] ModifiedZigguratNormalAliasThreshold = { 0xc6d1bcb1d7c5c800, 0xc58e65fed90c1800, 0xbdfe727f7b510000, 0xdf559b5c3fd23000, 0xfd37a81820b17000, 0xda845732a9b27000, 0xc18bb1c5c6944800, 0xaea2e0fe0c50f000, 0x9fc6849420a6c800, 0x93c2f42cbec16000, 0x89d49e0c64c58800, 0x8178f3a67d05e000, 0x7a554ee0b0451000, 0x7428aefe15e7f400, 0x6ec330fc6a914000, 0x6a00bca9d2052000, 0x65c592cf1a789400, 0x61fc00beb1018800, 0x5e92cd0eec7fc800, 0x5b7c1cb2569b6000, 0x58aca85e4418b000, 0x561b28c74e477c00, 0x53bfe910a6325000, 0x5194745ca0959000, 0x4f9356ee453d7000, 0x4db7eda0d49b9000, 0x4bfe400809666c00, 0x4a62e28fed3b5800, 0x48e2deb6dce1f000, 0x477b9ff689033000, 0x462ae44f1e19e800, 0x44eeafaec37ed800, 0x43c54198bbb85000, 0x42ad0c96e746f400, 0x41a4af1f12039800, 0x40aaed9f2142d000, 0x3fbead7e1c265000, 0x3edef0e17972b000, 0x3e0ad317806c0200, 0x3d41858782572800, 0x3c824d1301ff9200, 0x3bcc7fd4ba1b3400, 0x3b1f832e92fdb800, 0x3a7aca1973b67400, 0x39ddd3aedfad2400, 0x394829e0df990600, 0x38b9605bcc8b1a00, 0x3831138b23cd2800, 0x37aee7bbd833a600, 0x37328858ffdef400, 0x36bba73f64732600, 0x3649fc239c9b0c00, 0x35dd440a3ca3e400, 0x357540cd98fcdc00, 0x3511b8b1b6f68c00, 0x34b27603192b4400, 0x345746bed152d200, 0x33fffc451e8a3a00, 0x33ac6b11ce637400, 0x335c6a7d78297200, 0x330fd4834a7bf000, 0x32c6858cd99e9e00, 0x32805c436b10a600, 0x323d39647de4da00, 0x31fcff9b1cbc2e00, 0x31bf935c0afec200, 0x3184dac54d663c00, 0x314cbd80aea4b600, 0x311724a8b1030000, 0x30e3faaf5c30f400, 0x30b32b47f4eed400, 0x3084a3519d16c800, 0x305850c4757e6600, 0x302e229fe376a200, 0x300608da4510a600, 0x2fdff451a1855a00, 0x2fbbd6be63f1ee00, 0x2f99a2a64bd20400, 0x2f794b507d578400, 0x2f5ac4babd80a400, 0x2f3e038f49a68000, 0x2f22fd1b9810c400, 0x2f09a74772a6d000, 0x2ef1f88d34f90200, 0x2edbe7f1c515c600, 0x2ec76cfe738d8800, 0x2eb47fba0c19da00, 0x2ea318a28e748e00, 0x2e9330a894d4ca00, 0x2e84c1291fc1ca00, 0x2e77c3e97a4dce00, 0x2e6c331290318400, 0x2e62092c7d15f200, 0x2e59411b71ded400, 0x2e51d61b49fe3e00, 0x2e4bc3bccd116600, 0x2e4705e25f1ef200, 0x2e4398bd39fd9400, 0x2e4178cab1778a00, 0x2e40a2d1fb9f2c00, 0x2e4113e198c4c800, 0x2e42c94de1bd5a00, 0x2e45c0ae0c802000, 0x2e49f7dbb2dc2000, 0x2e4f6cf035308e00, 0x2e561e43fc4bbe00, 0x2e5e0a6cbd8db000, 0x2e67303c4fe53000, 0x2e718ebfe48e8600, 0x2e7d253f029ffa00, 0x2e89f33a94dd2400, 0x2e97f86c0ae26a00, 0x2ea734c573754c00, 0x2eb7a8703ea7da00, 0x2ec953cd28034200, 0x2edc377436e01c00, 0x2ef0543416c83200, 0x2f05ab123d28a200, 0x2f1c3d4b2da8d400, 0x2f340c5207034a00, 0x2f4d19d1325b6e00, 0x2f6767aadb615a00, 0x2f82f7f87920a000, 0x2f9fcd0c79737a00, 0x2fbde971df0bb200, 0x2fdd4fedcdd7ba00, 0x2ffe037f77257600, 0x30200761cc76da00, 0x30435f0c42b70800, 0x30680e33921a3000, 0x308e18cbac79aa00, 0x30b5830850abb800, 0x30de515f4a056600, 0x310888894f49a600, 0x31342d8449113400, 0x31614594c1269a00, 0x318fd6486c0d4c00, 0x31bfe57779f8b400, 0x31f17947e4e34000, 0x3224982ee3906c00, 0x325948f43c891000, 0x328f92b4cc68f200, 0x32c77ce56f2f8e00, 0x33010f562b3bd600, 0x333c5235b88b5400, 0x33794e1501d94a00, 0x33b80beac8399600, 0x33f89517e4335800, 0x343af36b92a06600, 0x347f3127c9566c00, 0x34c5590668495400, 0x350d763e28633000, 0x3557948851df6e00, 0x35a3c02639f2e200, 0x35f205e7cfddc800, 0x3642733226acc400, 0x3695160663fbc800, 0x36e9fd0946300200, 0x3741378b39d70800, 0x379ad590b4f46800, 0x37f6e7db2e059c00, 0x38557ff35d637c00, 0x38b6b03277c0ee00, 0x391a8bce450fde00, 0x398126e4252d2000, 0x39ea9685db846400, 0x3a56f0c70e33d800, 0x3ac64ccb52104200, 0x3b38c2d5fabeb200, 0x3bae6c5a0bd7de00, 0x3c27640c40341400, 0x3ca3c5f5b1b4f200, 0x3d23af886fd6bc00, 0x3da73fb5118a5a00, 0x3e2e97023a62d600, 0x3eb9d7a61501d000, 0x3f4925a1214e6800, 0x3fdca6dc15b5f400, 0x40748346fbe7a800, 0x4110e4fb856bb800, 0x41b1f8627a85d000, 0x4257ec5b004f0800, 0x4302f26661d57400, 0x43b33ed7116e4000, 0x4469090393bba400, 0x45248b7deb826400, 0x45e6045012d9e000, 0x46adb53da6aed400, 0x477be40bd0fa2000, 0x4850dacf591c1c00, 0x492ce842a0fd7000, 0x4a1060233adb1000, 0x4afb9b98ffb8dc00, 0x4beef9a70adef000, 0x4ceadfa81a509800, 0x4defb9d72c749400, 0x4efdfbe72fe4d000, 0x501621a9608c3400, 0x5138afc6a7a9f000, 0x5266348c2347b800, 0x539f48cf4aad7800, 0x54e490eb26298000, 0x5636bddb69fe9c00, 0x57968e77cf054800, 0x5904d0d697505400, 0x5a8263d962996000, 0x5c1038ebfe230400, 0x5daf55fbfd149000, 0x5f60d7b2ef00d400, 0x6125f3fa86025800, 0x62fffcdb2c4b1800, 0x64f063bf243f3800, 0x66f8bd2da9e8d800, 0x691ac5123b4ba800, 0x6b5863a76bef2400, 0x6db3b324abb4d800, 0x702f0650b62f2800, 0x72ccf025b46f0c00, 0x75904cbbd8328c00, 0x787c4bbbb09d3000, 0x7b947ca99e2e4c00, 0x7edcdd6e5388c400, 0x8259eb9b2bd39800, 0x8610b907d7b2c000, 0x8a070493b6c2f000, 0x8e435809cf16f000, 0x92cd2c7476d78800, 0x97ad168ac8e14000, 0x9cecfd6d8baee000, 0xa2985e97c2114000, 0xa8bca2e6cb415800, 0xaf6989f6a2190000, 0xb6b1b2fd462b6000, 0xbeab4d13c9c93000, 0xc770fcdd4cdad000, 0xd1230b71f7ef4800, 0xdbe8fb696740b000, 0xe7f3aeac31d89800, 0xf5805d6770e7f000, 0xfffffffffd71c000, 0xfdeb966ea7222800, 0xf1fe2559fa518800, 0xf7bdcf2ab32e5000, 0xe4368349dbba4000, 0xfcccff52b1bd7000, 0xcf0a381011797800, 0xce3286a74622e000, 0xae75adb6e2295000, 0xd722c4a089ebd800, 0xc9dc5f7b7d2d8000, 0xd7c9bdcbf9301000, 0xcfef53b6a4d3c800, 0xb786b6325d47c000, 0xb06b5777e76e9800, 0x0000000000000000, 0x0000000000000000 };
private static readonly byte[] ModifiedZigguratNormalAliasIndex = { 0xfb, 0xf7, 0xf4, 0xf1, 0xef, 0x03, 0xf3, 0xf5, 0xf6, 0xf7, 0xf8, 0xf9, 0xf9, 0xfa, 0xfb, 0xfb, 0xfb, 0x00, 0x00, 0xfc, 0xfc, 0xfc, 0xfc, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfc, 0xfc, 0xfc, 0xfc, 0xfc, 0xfc, 0xfc, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0xfb, 0xfb, 0xfb, 0xfb, 0xfa, 0xfa, 0xfa, 0xfa, 0xf9, 0xf9, 0xf8, 0xf8, 0x01, 0x01, 0x01, 0xf7, 0xf7, 0xf6, 0xf5, 0x02, 0x02, 0xf4, 0xf3, 0xf2, 0x03, 0xf1, 0xf0, 0xef, 0xef, 0xf0, 0x03, 0xf2, 0xf3, 0x02, 0xf5, 0xf6, 0x01, 0xf8, 0xf9, 0xfa, 0x00, 0xfc, 0xfd, 0xfd };
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static double ModifiedZigguratNormalSeparate<TRng>(TRng rng)
where TRng : notnull, RandomBase
{
const int TableSize = 256;
const int IMax = 253;
long u1 = (long)rng.Next();
int tableX = (int)(u1 & (TableSize - 1));
if (tableX < IMax)
{
return u1 * ModifiedZigguratNormalX[tableX];
}
[MethodImpl(MethodImplOptions.NoInlining)]
static double Fallback(TRng rng, long u1)
{
const int InflectionIndex = 204;
const ulong MaxConvexRate = 0x3efb83be6450cc00;
const ulong MaxConcaveRate = 0x151b6b6b7cd81f00;
[MethodImpl(MethodImplOptions.AggressiveInlining)]
static double SampleX(double x, int tableY) => Math.FusedMultiplyAdd(x, ModifiedZigguratNormalX[tableY - 1] - ModifiedZigguratNormalX[tableY], ModifiedZigguratNormalX[tableY] * (1ul << 63));
[MethodImpl(MethodImplOptions.AggressiveInlining)]
static double SampleY(double y, int tableY) => Math.FusedMultiplyAdd(y, ModifiedZigguratNormalY[tableY - 1] - ModifiedZigguratNormalY[tableY], ModifiedZigguratNormalY[tableY] * (1ul << 63));
ulong u2 = rng.Next();
int tableAlias = (int)(u2 & (TableSize - 1));
int tableY = u2 < ModifiedZigguratNormalAliasThreshold[tableAlias] ? tableAlias : (ModifiedZigguratNormalAliasIndex[tableAlias] & (TableSize - 1));
if (tableY == 0)
{
// 最下層レイヤー: 尾から生成する
double x0 = ModifiedZigguratNormalX[0] * (1ul << 63);
double s, t;
do
{
s = ModifiedZigguratExponential(rng) / x0;
t = ModifiedZigguratExponential(rng);
} while (s * s > t + t);
return (x0 + s) * (u1 < 0 ? -1 : 1);
}
else if (tableY < InflectionIndex)
{
// へこみレイヤー
ulong ux = (ulong)u1 << 1 >> 1;
ulong uy = rng.Next() >> 1;
double tx;
for (; ; ux = rng.Next() >> 1, uy = rng.Next() >> 1)
{
// 三角形の右上; 上下左右反転させる
if (uy < ux)
{
ux = (1ul << 63) - 1 - ux;
uy = (1ul << 63) - 1 - uy;
}
// 三角形の左下; 採択
if (uy >= ux + MaxConcaveRate)
{
tx = SampleX(ux, tableY);
break;
}
// 確率密度関数より下なら採択
tx = SampleX(ux, tableY);
double ty = SampleY(uy, tableY);
if (ty <= Math.Exp(-0.5 * tx * tx))
{
break;
}
}
return tx * (u1 < 0 ? -1 : 1);
}
else if (tableY == InflectionIndex)
{
// へこみ+ふくらみレイヤー
ulong ux = (ulong)u1 << 1 >> 1;
ulong uy = rng.Next() >> 1;
double tx;
for (; ; ux = rng.Next() >> 1, uy = rng.Next() >> 1)
{
// 三角形の左下; 採択
if (uy >= ux + MaxConcaveRate)
{
tx = SampleX(ux, tableY);
break;
}
// 三角形の右上; 棄却
if (uy + MaxConvexRate <= ux)
{
continue;
}
// 確率密度関数より下なら採択
tx = SampleX(ux, tableY);
double ty = SampleY(uy, tableY);
if (ty <= Math.Exp(-0.5 * tx * tx))
{
break;
}
}
return tx * (u1 < 0 ? -1 : 1);
}
else
{
// ふくらみレイヤー
ulong ux = (ulong)u1 << 1 >> 1;
ulong uy = rng.Next() >> 1;
double tx;
for (; ; ux = rng.Next() >> 1, uy = rng.Next() >> 1)
{
// 三角形の左下; 採択
if (uy >= ux)
{
tx = SampleX(ux, tableY);
break;
}
// 三角形の右上; 棄却
if (uy + MaxConvexRate <= ux)
{
continue;
}
// 確率密度関数より下なら採択
tx = SampleX(ux, tableY);
double ty = SampleY(uy, tableY);
if (ty <= Math.Exp(-tx * tx / 2))
{
break;
}
}
return tx * (u1 < 0 ? -1 : 1);
}
}
return Fallback(rng, u1);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static double ModifiedZigguratNormalSimple<TRng>(TRng rng)
where TRng : notnull, RandomBase
{
const int TableSize = 256;
const int IMax = 253;
long u1 = (long)rng.Next();
int tableX = (int)(u1 & (TableSize - 1));
if (tableX < IMax)
{
return u1 * ModifiedZigguratNormalX[tableX];
}
[MethodImpl(MethodImplOptions.NoInlining)]
static double Fallback(TRng rng, long u1)
{
const ulong MaxConvexRate = 0x3efb83be6450cc00;
const ulong MaxConcaveRate = 0x151b6b6b7cd81f00;
[MethodImpl(MethodImplOptions.AggressiveInlining)]
static double SampleX(double x, int tableY) => Math.FusedMultiplyAdd(x, ModifiedZigguratNormalX[tableY - 1] - ModifiedZigguratNormalX[tableY], ModifiedZigguratNormalX[tableY] * (1ul << 63));
[MethodImpl(MethodImplOptions.AggressiveInlining)]
static double SampleY(double y, int tableY) => Math.FusedMultiplyAdd(y, ModifiedZigguratNormalY[tableY - 1] - ModifiedZigguratNormalY[tableY], ModifiedZigguratNormalY[tableY] * (1ul << 63));
ulong u2 = rng.Next();
int tableAlias = (int)(u2 & (TableSize - 1));
int tableY = (u2 >> 8) < (ModifiedZigguratNormalAliasBox[tableAlias] >> 8) ? tableAlias : (int)(ModifiedZigguratNormalAliasBox[tableAlias] & (TableSize - 1));
if (tableY == 0)
{
// 最下層レイヤー: 尾から生成する
double x0 = ModifiedZigguratNormalX[0] * (1ul << 63);
double s, t;
do
{
s = ModifiedZigguratExponential(rng) / x0;
t = ModifiedZigguratExponential(rng);
} while (s * s > t + t);
return (x0 + s) * (u1 < 0 ? -1 : 1);
}
else
{
// へこみ+ふくらみレイヤー
ulong ux = (ulong)u1 << 1 >> 1;
ulong uy = rng.Next() >> 1;
double tx;
for (; ; ux = rng.Next() >> 1, uy = rng.Next() >> 1)
{
// 三角形の左下; 採択
if (uy >= ux + MaxConcaveRate)
{
tx = SampleX(ux, tableY);
break;
}
// 三角形の右上; 棄却
if (uy + MaxConvexRate <= ux)
{
continue;
}
// 確率密度関数より下なら採択
tx = SampleX(ux, tableY);
double ty = SampleY(uy, tableY);
if (ty <= Math.Exp(-0.5 * tx * tx))
{
break;
}
}
return tx * (u1 < 0 ? -1 : 1);
}
}
return Fallback(rng, u1);
}
private static readonly double[] ModifiedZigguratNormal128X = { 3.7106324284260855E-19, 3.4727570413107269E-19, 3.323444774437097E-19, 3.2121224908146307E-19, 3.1222980429682549E-19, 3.0464149991809007E-19, 2.9803511076420591E-19, 2.921599105695414E-19, 2.8685152781782515E-19, 2.819961073506187E-19, 2.7751138553662935E-19, 2.7333590333967138E-19, 2.6942247998686363E-19, 2.6573407039455179E-19, 2.6224102984818063E-19, 2.5891924716988487E-19, 2.5574883427861844E-19, 2.5271318375598605E-19, 2.4979827662160572E-19, 2.46992164374329E-19, 2.4428457500930542E-19, 2.4166660891160468E-19, 2.3913050101329783E-19, 2.3666943255096012E-19, 2.3427738046333343E-19, 2.3194899571109774E-19, 2.2967950407422031E-19, 2.27464624601461E-19, 2.2530050205613187E-19, 2.2318365055805198E-19, 2.2111090625550372E-19, 2.1907938733570644E-19, 2.1708646004152319E-19, 2.1512970963653542E-19, 2.1320691547215583E-19, 2.1131602947487912E-19, 2.0945515750059417E-19, 2.0762254310455146E-19, 2.0581655335638063E-19, 2.0403566639418661E-19, 2.0227846046377408E-19, 2.0054360423116331E-19, 1.9882984819084119E-19, 1.9713601702024109E-19, 1.9546100275400929E-19, 1.938037586706696E-19, 1.9216329380010643E-19, 1.9053866797345862E-19, 1.8892898734803279E-19, 1.8733340034909345E-19, 1.8575109397817676E-19, 1.84181290444157E-19, 1.8262324407887132E-19, 1.8107623850384662E-19, 1.7953958401870402E-19, 1.7801261518525667E-19, 1.7649468858425159E-19, 1.7498518072421474E-19, 1.7348348608400093E-19, 1.7198901527247684E-19, 1.7050119329032078E-19, 1.6901945788023964E-19, 1.6754325795301234E-19, 1.6607205207769148E-19, 1.6460530702505223E-19, 1.6314249635398298E-19, 1.61683099030978E-19, 1.6022659807322798E-19, 1.587724792060147E-19, 1.573202295252039E-19, 1.5586933615559802E-19, 1.544192848957532E-19, 1.5296955883958015E-19, 1.5151963696462781E-19, 1.5006899267637876E-19, 1.48617092297157E-19, 1.4716339348733594E-19, 1.4570734358542019E-19, 1.4424837785222545E-19, 1.42785917602761E-19, 1.4131936820748498E-19, 1.3984811694229824E-19, 1.3837153066389848E-19, 1.3688895328385253E-19, 1.3539970301085422E-19, 1.3390306932599166E-19, 1.3239830965029308E-19, 1.3088464565715265E-19, 1.2936125917421152E-19, 1.278272876095657E-19, 1.2628181882539626E-19, 1.2472388536775736E-19, 1.2315245794366204E-19, 1.2156643801492993E-19, 1.1996464935141011E-19, 1.1834582835272523E-19, 1.1670861290570431E-19, 1.1505152949164969E-19, 1.1337297819012335E-19, 1.1167121513942944E-19, 1.099443319021019E-19, 1.0819023103773537E-19, 1.0640659699317903E-19, 1.0459086116407605E-19, 1.0274015963700992E-19, 1.008512816517032E-19, 9.8920606173870359E-20, 9.69440230601868E-20, 9.4916834002649629E-20, 9.28336265649053E-20, 9.068811185467488E-20, 8.84729121995528E-20, 8.6179278738656915E-20, 8.3796708594045952E-20, 8.1312414505672258E-20, 7.8710571352560578E-20, 7.5971213692188272E-20, 7.3068565203215115E-20, 6.9968397782173119E-20, 6.66236322285616E-20, 6.2966505034176623E-20, 5.8893338600748748E-20, 5.4231068392380193E-20, 4.8648631701179604E-20, 1.1815783248660808E-20, 0 };
private static readonly double[] ModifiedZigguratNormal128Y = { 3.1018627649684215E-22, 6.4161955637510225E-22, 9.8794312560512861E-22, 1.3462692043002265E-21, 1.7149038575531662E-21, 2.092720818384637E-21, 2.4789126511088548E-21, 2.8728706223370715E-21, 3.2741190539483217E-21, 3.6822761974052209E-21, 4.09702936300823E-21, 4.5181183024263428E-21, 4.9453236553271337E-21, 5.3784586517967231E-21, 5.8173629899501654E-21, 6.2618982143503573E-21, 6.7119441587115478E-21, 7.1673961613792455E-21, 7.6281628536032465E-21, 8.0941643801446191E-21, 8.5653309515002216E-21, 9.0416016541867231E-21, 9.522923464472197E-21, 1.000925042440781E-20, 1.0500542948741835E-20, 1.0996767238436776E-20, 1.1497894781820016E-20, 1.2003901928398776E-20, 1.2514769523418681E-20, 1.3030482593593847E-20, 1.3551030076263881E-20, 1.4076404585668654E-20, 1.4606602211168553E-20, 1.5141622343146861E-20, 1.5681467523061883E-20, 1.6226143314710132E-20, 1.6775658194246422E-20, 1.7330023456905948E-20, 1.7889253138703671E-20, 1.8453363951662365E-20, 1.9022375231352259E-20, 1.9596308895721431E-20, 2.0175189414363279E-20, 2.0759043787511009E-20, 2.1347901534173894E-20, 2.1941794688938658E-20, 2.2540757807055911E-20, 2.3144827977517197E-20, 2.3754044843906076E-20, 2.4368450632877551E-20, 2.4988090190185977E-20, 2.5613011024243715E-20, 2.6243263357251694E-20, 2.6878900184000486E-20, 2.7519977338496977E-20, 2.81665535686277E-20, 2.8818690619127371E-20, 2.9476453323179411E-20, 3.0139909703036017E-20, 3.0809131080109693E-20, 3.1484192195055722E-20, 3.2165171338437835E-20, 3.2852150492648107E-20, 3.3545215485837615E-20, 3.4244456158707755E-20, 3.4949966545115331E-20, 3.566184506755843E-20, 3.6380194748736775E-20, 3.7105123440521248E-20, 3.7836744071825718E-20, 3.8575174917051726E-20, 3.9320539886976922E-20, 4.007296884418466E-20, 4.083259794538864E-20, 4.159957001329878E-20, 4.2374034941007474E-20, 4.3156150132256428E-20, 4.3946080981381527E-20, 4.4744001397236093E-20, 4.5550094375973711E-20, 4.636455262824389E-20, 4.7187579267134457E-20, 4.8019388564104405E-20, 4.8860206781214487E-20, 4.97102730892105E-20, 5.0569840582484182E-20, 5.1439177403673269E-20, 5.2318567992724419E-20, 5.3208314477699412E-20, 5.4108738227546543E-20, 5.5020181590595773E-20, 5.5943009846810014E-20, 5.6877613407015433E-20, 5.7824410298671146E-20, 5.8783848985521609E-20, 5.9756411578088429E-20, 6.0742617503911476E-20, 6.1743027721409269E-20, 6.2758249580088645E-20, 6.37889424537937E-20, 6.4835824304371336E-20, 6.5899679372783664E-20, 6.6981367246417778E-20, 6.8081833619494419E-20, 6.920212315426687E-20, 7.0343394973099032E-20, 7.1506941478666208E-20, 7.2694211431046214E-20, 7.3906838536149409E-20, 7.5146677265941156E-20, 7.6415848310395542E-20, 7.7716797072911523E-20, 7.90523701634132E-20, 8.0425917258201741E-20, 8.18414295922275E-20, 8.3303732859302562E-20, 8.4818763619618413E-20, 8.6393978971419646E-20, 8.8038989137302267E-20, 8.9766585213831557E-20, 9.1594520964429774E-20, 9.3548879952675868E-20, 9.567125609340419E-20, 9.80371751669185E-20, 1.077782751194112E-19, 1.0842021724855044E-19 };
private static readonly ulong[] ModifiedZigguratNormal128AliasBox = { 0xd2b57fdbdb75d07b, 0xfb2212e9fcf5c878, 0xe58eccc587e88078, 0xb49989b3d043807a, 0x9786bb871cea787d, 0x84151812b4e9e87b, 0x7615dc89a6221c00, 0x6b7d90dc5ec4807c, 0x632c46853100787c, 0x5c75ce5e9235647c, 0x56ec890b8b05ac7c, 0x52472e51ffe1887c, 0x4e52cc819833cc7c, 0x4aead17032fb247c, 0x47f4459001fc087c, 0x455ad1136e0b1c7c, 0x430ecdd5b1fe387c, 0x4103fd57a79a347c, 0x3f30a6610665ce7c, 0x3d8cf610fbf38e7c, 0x3c128df58cb6d27c, 0x3abc30e71a3d107c, 0x39858558bc9e5c7c, 0x386ae6e5c8a6e07c, 0x376942e15f69ac7c, 0x367dfcf6debf0e7c, 0x35a6d9ca574d447c, 0x34e1ee1b848e887c, 0x342d9158fcc1487c, 0x338852d7e48ee67c, 0x32f0f1198714b27c, 0x326652ad7a76ae7c, 0x31e78059c6cb267c, 0x3173a04725eb2a7c, 0x3109f1fe24dc967c, 0x30a9cb0d21693e7c, 0x30529437d282bc7c, 0x3003c717a507ce7c, 0x2fbcec18ad96b07c, 0x2f7d98c41f02bc7c, 0x2f456e4a706a5c7c, 0x2f1418439a8a0a7c, 0x2ee94b9cc7b1367c, 0x2ec4c5ac10a7e47c, 0x2ea64b652b19c47c, 0x2e8da8a9fdac9e7c, 0x2e7aafb345d2c47c, 0x2e6d388df84e447c, 0x2e6520aace715c7c, 0x2e624a7d96db647c, 0x2e649d2a914f227c, 0x2e6c043ff9589e7c, 0x2e786f7ad17b827c, 0x2e89d295735e647c, 0x2ea0251f4be7ac7c, 0x2ebb625cceb1747c, 0x2edb892f3e32a87c, 0x2f009c037094ca7c, 0x2f2aa0c7d8b64e7c, 0x2f59a0e8b8dda47c, 0x2f8da9535a06a67c, 0x2fc6ca7f15d3427c, 0x3005187d2e41227c, 0x3048ab0f2ce7b67c, 0x30919dc3c3f4b67c, 0x30e0101ae6076a7c, 0x313425b12f36a47c, 0x318e067328116e7c, 0x31edded92a53f47c, 0x3253e02c4c43827c, 0x32c040d55d74227c, 0x33333cb72ee3b47c, 0x33ad1594dc735e7c, 0x342e1385c418287c, 0x34b685791bd2a07c, 0x3546c1ca4a43127c, 0x35df26e8f839be7c, 0x36801c17458e627c, 0x372a1240c3cdb67c, 0x37dd84ee96d7cc7c, 0x389afb5c464f5a7c, 0x396309b2b583147c, 0x3a36526ecdc29c7c, 0x3b1587fb608f007c, 0x3c016e86049ca47c, 0x3cfade198583ee7c, 0x3e02c5099a700a7c, 0x3f1a2abe4e15267c, 0x404232f104174c7c, 0x417c217018a3047c, 0x42c95e831616207c, 0x442b7c1053edd07c, 0x45a43baa51300c7c, 0x473595b72ad4547c, 0x48e1c1eece6f407c, 0x4aab417db9ddbc7c, 0x4c94eb2e0ef5d87c, 0x4ea1fa13b509787c, 0x50d61f5cdb5ff07c, 0x533598167fd0ec7c, 0x55c547f6274a307c, 0x588ada905afad47c, 0x5b8cecdc21778c7c, 0x5ed3418b3b39607c, 0x626703aa7fc2407c, 0x66531c50b262bc7c, 0x6aa4a20b9d5acc7c, 0x6f6b6b9126f49000, 0x74bad37afa7c9800, 0x7aaac15ad1935800, 0x815916b1957a007b, 0x88ebaec41b86307b, 0x91933b6d9ffa607d, 0x9b8f77433e6d407a, 0xa73576cddf06f87a, 0xb4f978040fde2079, 0xc57eae2e090f1801, 0xd9b1ce7564419078, 0xf2f80cd49743c877, 0xfffffffffc905877, 0xf9803610f3fef077, 0xe85e73ad7805c801, 0xb6ccadba631ed079, 0xdc889e670729e07d, 0x9a2a8a4b9eabc000, 0xe538d252165e007a, 0x000000000000007c, 0x000000000000007c };
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static double ModifiedZigguratNormal128<TRng>(TRng rng)
where TRng : notnull, RandomBase
{
const int TableSize = 128;
const int IMax = 125;
long u1 = (long)rng.Next();
int tableX = (int)(u1 & (TableSize - 1));
if (tableX < IMax)
{
return u1 * ModifiedZigguratNormal128X[tableX];
}
[MethodImpl(MethodImplOptions.NoInlining)]
static double Fallback(TRng rng, long u1)
{
const int InflectionIndex = 102;
const ulong MaxConvexRate = 0x3fdb883b8b97e800;
const ulong MaxConcaveRate = 0x14f98b67d5e73700;
[MethodImpl(MethodImplOptions.AggressiveInlining)]
static double SampleX(double x, int tableY) => Math.FusedMultiplyAdd(x, ModifiedZigguratNormal128X[tableY - 1] - ModifiedZigguratNormal128X[tableY], ModifiedZigguratNormal128X[tableY] * (1ul << 63));
[MethodImpl(MethodImplOptions.AggressiveInlining)]
static double SampleY(double y, int tableY) => Math.FusedMultiplyAdd(y, ModifiedZigguratNormal128Y[tableY - 1] - ModifiedZigguratNormal128Y[tableY], ModifiedZigguratNormal128Y[tableY] * (1ul << 63));
ulong u2 = rng.Next();
int tableAlias = (int)(u2 & (TableSize - 1));
int tableY = (u2 >> 8) < (ModifiedZigguratNormal128AliasBox[tableAlias] >> 8) ? tableAlias : (int)(ModifiedZigguratNormal128AliasBox[tableAlias] & (TableSize - 1));
if (tableY == 0)
{
// 最下層レイヤー: 尾から生成する
double x0 = ModifiedZigguratNormal128X[0] * (1ul << 63);
double s, t;
do
{
s = ModifiedZigguratExponential(rng) / x0;
t = ModifiedZigguratExponential(rng);
} while (s * s > t + t);
return (x0 + s) * (u1 < 0 ? -1 : 1);
}
else if (tableY < InflectionIndex)
{
// へこみレイヤー
ulong ux = (ulong)u1 << 1 >> 1;
ulong uy = rng.Next() >> 1;
double tx;
for (; ; ux = rng.Next() >> 1, uy = rng.Next() >> 1)
{
// 三角形の右上; 上下左右反転させる
if (uy < ux)
{
ux = (1ul << 63) - 1 - ux;
uy = (1ul << 63) - 1 - uy;
}
// 三角形の左下; 採択
if (uy >= ux + MaxConcaveRate)
{
tx = SampleX(ux, tableY);
break;
}
// 確率密度関数より下なら採択
tx = SampleX(ux, tableY);
double ty = SampleY(uy, tableY);
if (ty <= Math.Exp(-0.5 * tx * tx))
{
break;
}
}
return tx * (u1 < 0 ? -1 : 1);
}
else if (tableY == InflectionIndex)
{
// へこみ+ふくらみレイヤー
ulong ux = (ulong)u1 << 1 >> 1;
ulong uy = rng.Next() >> 1;
double tx;
for (; ; ux = rng.Next() >> 1, uy = rng.Next() >> 1)
{
// 三角形の左下; 採択
if (uy >= ux + MaxConcaveRate)
{
tx = SampleX(ux, tableY);
break;
}
// 三角形の右上; 棄却
if (uy + MaxConvexRate <= ux)
{
continue;
}
// 確率密度関数より下なら採択
tx = SampleX(ux, tableY);
double ty = SampleY(uy, tableY);
if (ty <= Math.Exp(-0.5 * tx * tx))
{
break;
}
}
return tx * (u1 < 0 ? -1 : 1);
}
else
{
// ふくらみレイヤー
ulong ux = (ulong)u1 << 1 >> 1;
ulong uy = rng.Next() >> 1;
double tx;
for (; ; ux = rng.Next() >> 1, uy = rng.Next() >> 1)
{
// 三角形の左下; 採択
if (uy >= ux)
{
tx = SampleX(ux, tableY);
break;
}
// 三角形の右上; 棄却
if (uy + MaxConvexRate <= ux)
{
continue;
}
// 確率密度関数より下なら採択
tx = SampleX(ux, tableY);
double ty = SampleY(uy, tableY);
if (ty <= Math.Exp(-tx * tx / 2))
{
break;
}
}
return tx * (u1 < 0 ? -1 : 1);
}
}
return Fallback(rng, u1);
}
public static double Inverse<T>(T rng)
where T : notnull, RandomBase
{
return -Math.Log(1 - ((rng.Next() >> 11) * (1.0 / (1ul << 53))));
}
private static readonly double[] ZigguratExponentialX = new double[] { 8.6971174701308573, 7.6971174701310288, 6.9410336293772019, 6.4783784938325635, 6.1441646657724682, 5.8821443157953963, 5.666410167454031, 5.4828906275260607, 5.3230905057543971, 5.1814872813015, 5.0542884899813041, 4.9387770859012514, 4.8329397410251129, 4.735242996601742, 4.6444918854200861, 4.5597370617073523, 4.4802117465284228, 4.4052876934735741, 4.3344436803172739, 4.2672424802773667, 4.2033137137351853, 4.1423408656640524, 4.0840513104082987, 4.0282085446479377, 3.9746060666737897, 3.9230625001354911, 3.8734176703995105, 3.8255294185223381, 3.7792709924116692, 3.7345288940397987, 3.69120109023742, 3.6491955157608551, 3.6084288131289108, 3.5688252656483388, 3.5303158891293451, 3.4928376547740614, 3.456332821132762, 3.4207483572511217, 3.3860354424603032, 3.3521490309001116, 3.3190474709707503, 3.2866921715990705, 3.2550473085704517, 3.2240795652862659, 3.1937579032122425, 3.1640533580259751, 3.1349388580844422, 3.1063890623398258, 3.078380215254092, 3.0508900166154569, 3.0238975044556784, 2.9973829495161324, 2.9713277599210914, 2.9457143948950475, 2.9205262865127426, 2.8957477686001436, 2.8713640120155381, 2.8473609656351906, 2.8237253024500375, 2.80044437025074, 2.7775061464397588, 2.7548991965623473, 2.7326126361947027, 2.7106360958679314, 2.6889596887418059, 2.6675739807732688, 2.6464699631518109, 2.6256390267977903, 2.6050729387408373, 2.5847638202141425, 2.5647041263169075, 2.5448866271118722, 2.5253043900378302, 2.5059507635285962, 2.4868193617402121, 2.4679040502973675, 2.4491989329782524, 2.4306983392644224, 2.4123968126888733, 2.3942890999214606, 2.3763701405361428, 2.3586350574093395, 2.341079147703037, 2.3236978743901986, 2.306486858283582, 2.2894418705322717, 2.272558825553157, 2.2558337743672214, 2.2392628983129113, 2.222842503111039, 2.2065690132576661, 2.1904389667232222, 2.1744490099377769, 2.1585958930438882, 2.1428764653998442, 2.1272876713173705, 2.1118265460190444, 2.0964902118017172, 2.0812758743932274, 2.0661808194905782, 2.0512024094685875, 2.0363380802487723, 2.0215853383189288, 2.0069417578945212, 1.9924049782135793, 1.9779727009573631, 1.963642687789551, 1.9494127580071874, 1.935280786297054, 1.9212447005915307, 1.9073024800183902, 1.8934521529393107, 1.8796917950722138, 1.8660195276928306, 1.852433515911178, 1.8389319670188822, 1.825513128903522, 1.8121752885263929, 1.7989167704602931, 1.785735935484128, 1.7726311792313076, 1.7596009308890768, 1.7466436519460764, 1.7337578349855738, 1.7209420025219375, 1.70819470587806, 1.6955145241015401, 1.6829000629175561, 1.6703499537164543, 1.657862852574175, 1.6454374393037259, 1.6330724165359938, 1.6207665088282603, 1.6085184617988608, 1.5963270412864858, 1.5841910325326913, 1.5721092393862321, 1.5600804835278905, 1.5481036037145159, 1.5361774550410345, 1.5243009082192287, 1.5124728488721195, 1.5006921768428192, 1.4889578055167483, 1.4772686611561361, 1.4656236822457476, 1.4540218188487959, 1.4424620319720149, 1.4309432929388821, 1.4194645827699857, 1.4080248915695381, 1.3966232179170446, 1.3852585682631247, 1.3739299563284932, 1.3626364025050894, 1.3513769332583379, 1.3401505805295075, 1.3289563811371192, 1.3177933761763276, 1.3066606104151768, 1.2955571316866037, 1.2844819902750155, 1.273434238296244, 1.2624129290696182, 1.2514171164808554, 1.2404458543344095, 1.2294981956938522, 1.2185731922087935, 1.2076698934267645, 1.1967873460884062, 1.1859245934042053, 1.1750806743109148, 1.164254622705682, 1.1534454666557779, 1.1426522275816759, 1.1318739194110818, 1.121109547701334, 1.1103581087274148, 1.0996185885326011, 1.0888899619385506, 1.0781711915113761, 1.0674612264799714, 1.0567590016025552, 1.046063435977048, 1.0353734317905323, 1.0246878730026212, 1.0140056239571005, 1.0033255279157007, 0.99264640550727978, 0.98196705308506638, 0.97128624098390726, 0.96060271166867039, 0.94991517776408, 0.93922231995526639, 0.9285227847472145, 0.91781518207004831, 0.90709808271569436, 0.896370015589894, 0.88562946476175564, 0.87487486629102917, 0.86410460481100859, 0.85331700984237746, 0.84251035181037282, 0.83168283773427742, 0.820832606554416, 0.80995772405742261, 0.79905617735549139, 0.78812586886949687, 0.777164609759134, 0.766170112735439, 0.75513998418198647, 0.74407171550051243, 0.73296267358436962, 0.72181009030876064, 0.71061105090965948, 0.69936248110323651, 0.68806113277375247, 0.67670356802952725, 0.66528614139268238, 0.65380497984766961, 0.642255960424541, 0.63063468493349506, 0.61893645139488085, 0.607156221620305, 0.59528858429150766, 0.58332771274877426, 0.57126731653259311, 0.55910058551154529, 0.54682012516331513, 0.53441788123717027, 0.52188505159213994, 0.50921198244365928, 0.49638804551867588, 0.48340149165346652, 0.47023927508217384, 0.45688684093142512, 0.44332786607355745, 0.42954394022541581, 0.41551416960036158, 0.40121467889628309, 0.386617977941125, 0.37169214532992267, 0.35639976025839937, 0.34069648106485489, 0.32452911701691545, 0.30783295467493849, 0.29052795549123694, 0.27251318547847153, 0.25365836338591918, 0.23379048305968231, 0.21267151063097472, 0.18995868962244072, 0.16512762256419711, 0.13730498094002377, 0.10483850756583221, 0.063852163815019666 };
private static readonly double[] ZigguratExponentialY = new double[] { 0.00045413435384150625, 0.00096726928232718462, 0.0015362997803015821, 0.0021459677437189167, 0.0027887987935740857, 0.0034602647778369132, 0.004157295120833804, 0.0048776559835424, 0.0056196422072054891, 0.0063819059373191834, 0.0071633531836349847, 0.0079630774380170365, 0.0087803149858089683, 0.0096144136425022047, 0.010464810181029972, 0.01133101359783459, 0.012212592426255367, 0.013109164931254979, 0.014020391403181929, 0.014945968011691135, 0.015885621839973142, 0.016839106826039927, 0.017806200410911341, 0.018786700744696006, 0.019780424338009712, 0.020787204072578086, 0.021806887504283553, 0.022839335406385209, 0.023884420511558143, 0.024942026419731752, 0.026012046645134186, 0.027094383780955765, 0.028188948763978594, 0.029295660224637358, 0.03041444391046657, 0.031545232172893567, 0.032687963508959493, 0.033842582150874281, 0.035009037697397355, 0.036187284781931367, 0.037377282772959312, 0.0385789955030748, 0.039792391023374063, 0.04101744138041475, 0.042254122413316164, 0.043502413568888121, 0.044762297732943226, 0.0460337610761751, 0.047316792913181478, 0.048611385573379413, 0.049917534282706288, 0.051235237055126191, 0.052564494593071595, 0.053905310196045983, 0.055257689676696933, 0.056621641283742766, 0.057997175631200555, 0.059384305633420148, 0.060783046445479529, 0.0621934154085409, 0.063615431999807209, 0.065049117786753624, 0.066494496385339649, 0.06795159342193649, 0.06942043649872863, 0.070901055162371718, 0.072393480875708627, 0.073897746992364607, 0.075413888734058271, 0.076941943170480351, 0.078481949201606255, 0.080033947542319725, 0.081597980709237239, 0.083174093009632175, 0.08476233053236791, 0.086362741140756691, 0.087975374467270009, 0.08960028191003265, 0.091237516631039961, 0.092887133556043361, 0.094549189376055651, 0.096223742550432576, 0.097910853311491991, 0.09961058367063691, 0.10132299742595341, 0.10304816017125748, 0.10478613930656992, 0.1065370040500014, 0.10830082545103352, 0.11007767640518512, 0.11186763167005603, 0.11367076788274404, 0.11548716357863324, 0.11731689921155526, 0.11916005717532738, 0.12101672182667453, 0.12288697950954483, 0.12477091858083066, 0.12666862943751034, 0.12858020454522784, 0.13050573846833044, 0.13244532790138716, 0.13439907170221324, 0.13636707092642847, 0.13834942886357982, 0.14034625107486204, 0.14235764543247181, 0.14438372216063433, 0.1464245938783445, 0.14848037564386635, 0.15055118500103945, 0.15263714202744241, 0.15473836938446761, 0.15685499236936479, 0.15898713896931377, 0.16113493991759159, 0.16329852875190137, 0.16547804187493556, 0.16767361861724975, 0.16988540130252722, 0.17211353531531962, 0.17435816917135305, 0.17661945459049444, 0.17889754657247786, 0.18119260347549584, 0.18350478709776702, 0.18583426276219667, 0.18818119940425385, 0.19054576966319492, 0.19292814997677085, 0.19532852067956272, 0.19774706610509835, 0.20018397469191074, 0.20263943909370846, 0.20511365629383715, 0.20760682772422148, 0.2101191593889877, 0.21265086199297772, 0.2152021510753781, 0.21777324714869994, 0.22036437584335888, 0.22297576805811958, 0.22560766011668346, 0.22826029393071612, 0.23093391716962683, 0.23362878343743271, 0.236345152457059, 0.23908329026244851, 0.24184346939887655, 0.24462596913189141, 0.24743107566532693, 0.25025908236886157, 0.25311029001562874, 0.2559850070304146, 0.25888354974901545, 0.2618062426893622, 0.26475341883506143, 0.267725419932044, 0.27072259679905925, 0.27374530965280219, 0.27679392844851652, 0.27986883323697209, 0.28297041453877997, 0.286099073737076, 0.28925522348967686, 0.29243928816189169, 0.29565170428126025, 0.29889292101558085, 0.30216340067569264, 0.30546361924458931, 0.30879406693455924, 0.31215524877417866, 0.31554768522712795, 0.31897191284495624, 0.32242848495608811, 0.32591797239355508, 0.32944096426413522, 0.33299806876180782, 0.33658991402867644, 0.34021714906677891, 0.3438804447045013, 0.34758049462163582, 0.35131801643748212, 0.35509375286678624, 0.35890847294874856, 0.36276297335481655, 0.36665807978151294, 0.37059464843514478, 0.37457356761590094, 0.37859575940957957, 0.38266218149600856, 0.38677382908413638, 0.39093173698479577, 0.39513698183328883, 0.39939068447522974, 0.40369401253052889, 0.408048183152031, 0.41245446599715979, 0.41691418643300149, 0.42142872899761519, 0.4259995411430329, 0.43062813728845734, 0.43531610321563508, 0.44006510084235234, 0.44487687341454696, 0.44975325116275339, 0.45469615747461384, 0.459707615642136, 0.46478975625042451, 0.46994482528395831, 0.47517519303737565, 0.48048336393045249, 0.48587198734188314, 0.4913438695940307, 0.49690198724154766, 0.50254950184134572, 0.50828977641064088, 0.51412639381474656, 0.52006317736823149, 0.52610421398361762, 0.5322538802630411, 0.53851687200285969, 0.54489823767243739, 0.55140341654063907, 0.55803828226258523, 0.564809192912398, 0.57172304866482349, 0.57878735860284269, 0.58601031847726559, 0.59340090169173088, 0.60096896636522967, 0.60872538207961946, 0.616682180915205, 0.6248527387036632, 0.63325199421436329, 0.6418967164272632, 0.6508058334145681, 0.6600008410789967, 0.66950631673192162, 0.67935057226476214, 0.68956649611707466, 0.70019265508278472, 0.71127476080507235, 0.72286765959356813, 0.73503809243141938, 0.74786862198519077, 0.76146338884989162, 0.77595685204011067, 0.79152763697249029, 0.80842165152300249, 0.82699329664304377, 0.84778550062398217, 0.87170433238119494, 0.9004699299257356, 0.93814368086215949, 1 };
private static readonly double[] ZigguratExponentialInsideRectangleRatio = new double[] { 0.88501937527759034, 0.90177052075821362, 0.93334492234895694, 0.94841088269568263, 0.95735460160488228, 0.96332389401564811, 0.96761273284061833, 0.970854767561948, 0.97339830605926758, 0.97545129720201773, 0.97714586250685509, 0.9785701312217, 0.97978523431731257, 0.98083496216629562, 0.98175154014612565, 0.982559232231441, 0.98327667144016018, 0.98391841394121538, 0.98449600340983368, 0.98501871716040235, 0.98549410007825688, 0.9859283537627439, 0.98632662483499056, 0.98669322171877949, 0.98703177983587354, 0.98734538903362712, 0.98763669297965662, 0.98790796748635723, 0.98816118281496423, 0.98839805366842159, 0.98862007962999832, 0.98882857811923941, 0.98902471143770987, 0.98920950910943584, 0.9893838864474741, 0.98954866007259057, 0.98970456095429549, 0.98985224542540939, 0.98999230452958, 0.99012527198992906, 0.9902516310312921, 0.99037182024466175, 0.9904862386477, 0.99059525006749294, 0.99069918694952086, 0.9907983536789351, 0.99089302948572933, 0.9909834709936054, 0.99106991446267267, 0.99115257776820054, 0.991231662151089, 0.99130735377031254, 0.99137982508307188, 0.99144923607463231, 0.99151573535666215, 0.99157946115023932, 0.99164054216744923, 0.99169909840360526, 0.99175524185050945, 0.99180907713980859, 0.99186070212431776, 0.99191020840419342, 0.99195768180396959, 0.99200320280572929, 0.99204684694304812, 0.99208868515978688, 0.99212878413733741, 0.99216720659349922, 0.99220401155581028, 0.99223925461183016, 0.99227298813859888, 0.99230526151325438, 0.9923361213065709, 0.99236561146099878, 0.99239377345461544, 0.99242064645225525, 0.99244626744495057, 0.99247067137870859, 0.99249389127353815, 0.992515958333559, 0.99253690204893674, 0.99255675029032031, 0.992575529396393, 0.99259326425509031, 0.99260997837898168, 0.992625693975279, 0.99264043201087893, 0.99265421227281769, 0.99266705342448014, 0.99267897305787656, 0.99268998774226858, 0.99270011306940664, 0.99270936369561391, 0.99271775338093615, 0.99272529502555107, 0.992732000703623, 0.99273788169476451, 0.99274294851326073, 0.99274721093519136, 0.99275067802357986, 0.99275335815168719, 0.992755259024552, 0.99275638769888164, 0.99275675060137647, 0.99275635354557468, 0.99275520174728593, 0.992753299838688, 0.992750651881144, 0.99274726137679492, 0.992743131278984, 0.99273826400155152, 0.992732661427047, 0.99272632491388946, 0.99271925530251526, 0.99271145292053586, 0.992702917586933, 0.99269364861531262, 0.99268364481623561, 0.99267290449863754, 0.9926614254703523, 0.99264920503774434, 0.99263624000445749, 0.99262252666928052, 0.99260806082312936, 0.99259283774514362, 0.99257685219788694, 0.99256009842164672, 0.992542570127816, 0.99252426049134423, 0.99250516214223872, 0.99248526715609164, 0.99246456704361208, 0.99244305273913, 0.99242071458804337, 0.99239754233317123, 0.99237352509997256, 0.9923486513805887, 0.99232290901666131, 0.99229628518087176, 0.99226876635714678, 0.99224033831946779, 0.99221098610921488, 0.99218069401097442, 0.99214944552672835, 0.9921172233483414, 0.992084009328251, 0.99204978444826, 0.99201452878632423, 0.99197822148121506, 0.99194084069493038, 0.99190236357271777, 0.991862766200558, 0.9918220235599513, 0.99178010947982875, 0.9917369965854046, 0.99169265624376, 0.99164705850594392, 0.9916001720453459, 0.99155196409208446, 0.99150240036313253, 0.991451444987864, 0.99139906042870551, 0.991345207396519, 0.99128984476033011, 0.9912329294509743, 0.991174416358196, 0.99111425822069354, 0.99105240550855578, 0.99098880629749, 0.99092340613417673, 0.99085614789203169, 0.99078697161658169, 0.99071581435959077, 0.99064261000097809, 0.99056728905748914, 0.99048977847696273, 0.99041000141693258, 0.99032787700616487, 0.99024332008758831, 0.99015624094091947, 0.99006654498309177, 0.98997413244440968, 0.98987889801810292, 0.989780730480718, 0.98967951228048034, 0.98957511909044193, 0.98946741932286209, 0.98935627360084655, 0.98924153418280314, 0.98912304433473308, 0.9890006376447642, 0.98887413727364459, 0.98874335513410516, 0.98860809099110691, 0.988468131473929, 0.98832324898986468, 0.988173200527905, 0.98801772633919471, 0.9878565484792, 0.9876893691943851, 0.98751586913370226, 0.98733570536230042, 0.987148509151455, 0.98695388351475388, 0.98675140045588272, 0.98654059788784931, 0.98632097617695813, 0.98609199425710281, 0.985853065250742, 0.98560355152190426, 0.98534275907338709, 0.98506993118442876, 0.98478424116596541, 0.98448478408730655, 0.9841705672997576, 0.98384049954802222, 0.98349337841764151, 0.98312787581408667, 0.98274252110380866, 0.98233568146602468, 0.98190553890171117, 0.98145006321710981, 0.98096698013499284, 0.98045373347714948, 0.979907440091483, 0.979324835846818, 0.97870221056073026, 0.97803532912240787, 0.97731933527055037, 0.97654863341021259, 0.97571674239409945, 0.97481611319623207, 0.9738378996382967, 0.97277166744713273, 0.97160502140444271, 0.97032312239453888, 0.96890805450552342, 0.96733798498544366, 0.96558603352123551, 0.96361872652488012, 0.96139384751146961, 0.95885738974131551, 0.95593914209661879, 0.95254613726150628, 0.94855265223826835, 0.94378444893278224, 0.93799298941024356, 0.93081134015791733, 0.9216746490790465, 0.9096670995657411, 0.89320233377217584, 0.86928175221887716, 0.83150825287660968, 0.76354482443449556, 0.60905258284914432, 0 };
public static double ZigguratExponential<TRng>(TRng rng)
where TRng : notnull, RandomBase
{
const int TableSize = 256;
double tail = 0;
while (true)
{
ulong u1 = rng.Next();
double u1Double = (u1 >> 11) * (1.0 / (1ul << 53));
int rectangleIndex = (int)(u1 & (TableSize - 1));
// 確率密度関数の範囲内にある四角形なら、そのまま返す
if (u1Double < ZigguratExponentialInsideRectangleRatio[rectangleIndex])
return tail + u1Double * ZigguratExponentialX[rectangleIndex];
if (rectangleIndex == 0)
{
// 一番下のレイヤー(尾の部分)
double x0 = ZigguratExponentialX[1];
tail += x0;
}
else
{
// それ以外
double x = u1Double * ZigguratExponentialX[rectangleIndex];
double y = (rng.Next() >> 11) * (1.0 / (1ul << 53));
if (ZigguratExponentialY[rectangleIndex - 1] + y * (ZigguratExponentialY[rectangleIndex] - ZigguratExponentialY[rectangleIndex - 1]) <= Math.Exp(-x))
return tail + x;
}
}
}
private static readonly double[] ModifiedZigguratExponentialX = { 4.10331203376756E-19, 3.6986866175804254E-19, 3.4566656688958066E-19, 3.2823679410482584E-19, 3.1456269979909514E-19, 3.03286120648027E-19, 2.9367638051868817E-19, 2.8529425264268634E-19, 2.7785472845811339E-19, 2.7116219451872107E-19, 2.6507648848254508E-19, 2.5949369628854142E-19, 2.5433461308999266E-19, 2.4953746469398331E-19, 2.4505312947224859E-19, 2.408418950532468E-19, 2.3687119326822434E-19, 2.3311397903598484E-19, 2.2954754508892096E-19, 2.2615263895329144E-19, 2.2291279408177043E-19, 2.1981381563184253E-19, 2.1684337983553293E-19, 2.1399071809234915E-19, 2.11246365135325E-19, 2.0860195626732108E-19, 2.0605006261232858E-19, 2.0358405612934665E-19, 2.0119799815503497E-19, 1.9888654671438722E-19, 1.9664487892667293E-19, 1.9446862564655205E-19, 1.9235381609360234E-19, 1.902968306909011E-19, 1.8829436069272401E-19, 1.8634337346015125E-19, 1.8444108246124117E-19, 1.8258492124400068E-19, 1.807725207664377E-19, 1.7900168957659048E-19, 1.7727039642266747E-19, 1.7557675494392153E-19, 1.739190101501551E-19, 1.7229552644453697E-19, 1.7070477698281687E-19, 1.691453341937061E-19, 1.6761586131144529E-19, 1.6611510479342966E-19, 1.6464188751402265E-19, 1.6319510264101051E-19, 1.6177370811405434E-19, 1.6037672165540421E-19, 1.5900321625239331E-19, 1.5765231605910447E-19, 1.5632319267132589E-19, 1.5501506173467129E-19, 1.5372717985068675E-19, 1.5245884175002801E-19, 1.5120937770547304E-19, 1.4997815116072294E-19, 1.4876455655371315E-19, 1.4756801731556631E-19, 1.4638798402842153E-19, 1.45223932727213E-19, 1.4407536333208373E-19, 1.4294179819953481E-19, 1.4182278078165824E-19, 1.4071787438389916E-19, 1.3962666101276579E-19, 1.3854874030576456E-19, 1.3748372853660133E-19, 1.3643125768936716E-19, 1.3539097459603044E-19, 1.3436254013209542E-19, 1.3334562846576738E-19, 1.3233992635639463E-19, 1.3134513249834234E-19, 1.3036095690679895E-19, 1.2938712034232588E-19, 1.28423353771241E-19, 1.2746939785917756E-19, 1.2652500249538757E-19, 1.2558992634556372E-19, 1.246639364311392E-19, 1.2374680773319347E-19, 1.2283832281924351E-19, 1.2193827149133938E-19, 1.210464504540078E-19, 1.2016266300070284E-19, 1.192867187175259E-19, 1.1841843320307338E-19, 1.1755762780335641E-19, 1.1670412936081657E-19, 1.1585776997653414E-19, 1.1501838678479183E-19, 1.1418582173921758E-19, 1.1335992140978604E-19, 1.1254053679000983E-19, 1.117275231136981E-19, 1.1092073968070402E-19, 1.1012004969112224E-19, 1.0932532008743433E-19, 1.085364214041337E-19, 1.077532276243935E-19, 1.06975616043369E-19, 1.0620346713775331E-19, 1.0543666444122949E-19, 1.0467509442548528E-19, 1.0391864638647763E-19, 1.0316721233565366E-19, 1.0242068689585317E-19, 1.0167896720163442E-19, 1.0094195280378056E-19, 1.0020954557775858E-19, 9.9481649635916368E-20, 9.8758171243215533E-20, 9.8039018736309816E-20, 9.7324102445789393E-20, 9.6613334621421653E-20, 9.5906629360228306E-20, 9.5203902537247469E-20, 9.4505071738837578E-20, 9.3810056198387454E-20, 9.31187767343039E-20, 9.2431155690154966E-20, 9.1747116876852892E-20, 9.1066585516766538E-20, 9.03894881896586E-20, 8.9715752780347442E-20, 8.9045308427998322E-20, 8.83780854769529E-20, 8.7714015429009779E-20, 8.7053030897072713E-20, 8.6395065560086262E-20, 8.5740054119181881E-20, 8.5087932254960365E-20, 8.44386365858392E-20, 8.3792104627395537E-20, 8.3148274752638179E-20, 8.2507086153143368E-20, 8.1868478800991455E-20, 8.1232393411442874E-20, 8.0598771406293152E-20, 7.9967554877848161E-20, 7.9338686553461617E-20, 7.8712109760577806E-20, 7.808776839222306E-20, 7.7465606872890166E-20, 7.6845570124760031E-20, 7.6227603534205167E-20, 7.5611652918519362E-20, 7.4997664492817877E-20, 7.438558483705183E-20, 7.3775360863079942E-20, 7.31669397817399E-20, 7.2560269069860586E-20, 7.1955296437155025E-20, 7.13519697929326E-20, 7.0750237212566977E-20, 7.015004690365451E-20, 6.95513471717952E-20, 6.8954086385926047E-20, 6.8358212943133333E-20, 6.7763675232867276E-20, 6.7170421600478691E-20, 6.6578400309993487E-20, 6.59875595060358E-20, 6.5397847174806131E-20, 6.4809211104014845E-20, 6.4221598841665558E-20, 6.3634957653576154E-20, 6.3049234479517678E-20, 6.2464375887843187E-20, 6.188032802846975E-20, 6.129703658406669E-20, 6.0714446719292248E-20, 6.0132503027908849E-20, 5.9551149477593744E-20, 5.8970329352247147E-20, 5.83899851915836E-20, 5.7810058727774453E-20, 5.723049081888938E-20, 5.6651221378862847E-20, 5.607218930368675E-20, 5.5493332393503676E-20, 5.4914587270244651E-20, 5.4335889290421791E-20, 5.3757172452648771E-20, 5.3178369299420046E-20, 5.2599410812633142E-20, 5.2020226302285734E-20, 5.1440743287720515E-20, 5.0860887370724853E-20, 5.0280582099717818E-20, 4.9699748824173369E-20, 4.9118306538333753E-20, 4.8536171713160077E-20, 4.7953258115345353E-20, 4.7369476612077128E-20, 4.6784734960079542E-20, 4.6198937577284764E-20, 4.561198529527826E-20, 4.5023775090426467E-20, 4.4434199791323838E-20, 4.3843147759883749E-20, 4.3250502543035532E-20, 4.2656142491570618E-20, 4.2059940342192631E-20, 4.1461762758256734E-20, 4.0861469824017271E-20, 4.0258914486419624E-20, 3.9653941937549631E-20, 3.9046388929762231E-20, 3.8436083014214539E-20, 3.7822841691982579E-20, 3.7206471465089576E-20, 3.658676677254669E-20, 3.5963508793815561E-20, 3.5336464098833411E-20, 3.47053831197502E-20, 3.4069998414628236E-20, 3.3430022687305141E-20, 3.2785146520105071E-20, 3.2135035766684288E-20, 3.14793285404618E-20, 3.0817631719071592E-20, 3.0149516866075868E-20, 2.9474515446425115E-20, 2.8792113179942989E-20, 2.8101743334798716E-20, 2.7402778706749675E-20, 2.6694521954501665E-20, 2.5976193858994976E-20, 2.5246918933169189E-20, 2.4505707611314757E-20, 2.3751433966683067E-20, 2.298280750063274E-20, 2.2198336948998792E-20, 2.1396283151074303E-20, 2.0574596636715019E-20, 1.9730833381303158E-20, 1.886203856570086E-20, 1.79645820431022E-20, 1.7033918345550304E-20, 1.6064223820782044E-20, 1.5047823458200018E-20, 1.3974234727799188E-20, 1.2828456524359345E-20, 1.1587604878401979E-20, 1.0213347614125671E-20, 8.63088516510677E-21, 1.8691277361760884E-21, 0 };
private static readonly double[] ModifiedZigguratExponentialY = { 2.7976027475557515E-23, 5.9012549913509879E-23, 9.2222116933672047E-23, 1.2719515233348412E-22, 1.6368847155753889E-22, 2.0153866066352558E-22, 2.4062739159746744E-22, 2.8086457448290812E-22, 3.2217910270220888E-22, 3.6451331171730957E-22, 4.0781944228160054E-22, 4.520572684174019E-22, 4.9719244243198682E-22, 5.4319530229844609E-22, 5.9003998877305356E-22, 6.3770377674155041E-22, 6.8616655881885455E-22, 7.3541043971875075E-22, 7.8541941287201191E-22, 8.3617909921871784E-22, 8.8767653375151537E-22, 9.39899989255219E-22, 9.928388293916143E-22, 1.0464833852026511E-21, 1.1008248504979007E-21, 1.1558551926152977E-21, 1.2115670758062617E-21, 1.2679537950710331E-21, 1.3250092187085155E-21, 1.3827277381830083E-21, 1.4411042241734185E-21, 1.5001339878773739E-21, 1.5598127468065068E-21, 1.6201365944400758E-21, 1.6811019732093431E-21, 1.74270565037044E-21, 1.8049446963929617E-21, 1.8678164655485767E-21, 1.931318578430991E-21, 1.9954489061776306E-21, 2.0602055561959357E-21, 2.125586859224434E-21, 2.1915913575816748E-21, 2.2582177944755211E-21, 2.3254651042617279E-21, 2.3933324035547874E-21, 2.4618189831059857E-21, 2.5309243003739369E-21, 2.6006479727217234E-21, 2.6709897711824342E-21, 2.7419496147415344E-21, 2.8135275650903038E-21, 2.8857238218095836E-21, 2.9585387179475211E-21, 3.0319727159588386E-21, 3.1060264039765707E-21, 3.1807004923902052E-21, 3.2559958107068082E-21, 3.331913304674072E-21, 3.4084540336463018E-21, 3.4856191681762061E-21, 3.5634099878170283E-21, 3.6418278791210032E-21, 3.7208743338214959E-21, 3.8005509471873057E-21, 3.8808594165387593E-21, 3.9618015399161166E-21, 4.0433792148917293E-21, 4.1255944375181539E-21, 4.2084493014051493E-21, 4.2919459969191414E-21, 4.37608681049931E-21, 4.4608741240850228E-21, 4.5463104146498117E-21, 4.6323982538375496E-21, 4.7191403076969E-21, 4.8065393365105021E-21, 4.894598194715693E-21, 4.983319830913927E-21, 5.0727072879663022E-21, 5.1627637031729616E-21, 5.25349230853432E-21, 5.344896431092389E-21, 5.4369794933506531E-21, 5.5297450137711825E-21, 5.623196607347895E-21, 5.7173379862550425E-21, 5.8121729605702182E-21, 5.9077054390713131E-21, 6.0039394301070831E-21, 6.1008790425410963E-21, 6.1985284867690043E-21, 6.2968920758092644E-21, 6.3959742264675574E-21, 6.4957794605752805E-21, 6.5963124063026928E-21, 6.6975777995473823E-21, 6.7995804853988689E-21, 6.9023254196803455E-21, 7.0058176705686239E-21, 7.1100624202935653E-21, 7.2150649669183374E-21, 7.3208307262020838E-21, 7.427365233546625E-21, 7.5346741460290269E-21, 7.64276324452201E-21, 7.7516384359042979E-21, 7.8613057553631861E-21, 7.9717713687917565E-21, 8.0830415752833376E-21, 8.1951228097259628E-21, 8.308021645499782E-21, 8.4217447972805243E-21, 8.5362991239523415E-21, 8.6516916316335227E-21, 8.7679294768187883E-21, 8.8850199696421055E-21, 9.00297057726413E-21, 9.1217889273886856E-21, 9.24148281191289E-21, 9.3620601907158E-21, 9.4835291955907153E-21, 9.6058981343265841E-21, 9.7291754949442317E-21, 9.8533699500934251E-21, 9.9784903616171675E-21, 1.010454578528994E-20, 1.0231545475736935E-20, 1.0359498891541785E-20, 1.0488415700550663E-20, 1.0618305785381053E-20, 1.0749179249143974E-20, 1.0881046421388921E-20, 1.1013917864281284E-20, 1.1147804379022596E-20, 1.1282717012524507E-20, 1.1418667064347987E-20, 1.155566609391999E-20, 1.1693725928040416E-20, 1.1832858668693039E-20, 1.197307670117479E-20, 1.2114392702558688E-20, 1.2256819650506589E-20, 1.2400370832448864E-20, 1.2545059855149203E-20, 1.2690900654673779E-20, 1.2837907506785232E-20, 1.2986095037783149E-20, 1.3135478235814109E-20, 1.3286072462675743E-20, 1.3437893466140902E-20, 1.3590957392829556E-20, 1.3745280801657969E-20, 1.3900880677896509E-20, 1.4057774447869568E-20, 1.4215979994333246E-20, 1.43755156725689E-20, 1.4536400327233133E-20, 1.4698653310007725E-20, 1.486229449809581E-20, 1.5027344313614035E-20, 1.5193823743933803E-20, 1.5361754363028524E-20, 1.5531158353887937E-20, 1.570205853206498E-20, 1.5874478370425466E-20, 1.604844202517616E-20, 1.6223974363252439E-20, 1.6401100991152985E-20, 1.6579848285315668E-20, 1.6760243424136097E-20, 1.6942314421738426E-20, 1.7126090163616655E-20, 1.7311600444274307E-20, 1.7498876007000826E-20, 1.7687948585934524E-20, 1.7878850950574509E-20, 1.8071616952917891E-20, 1.8266281577413691E-20, 1.8462880993941777E-20, 1.8661452614043478E-20, 1.8862035150651046E-20, 1.9064668681585509E-20, 1.9269394717117605E-20, 1.9476256271913919E-20, 1.9685297941721188E-20, 1.9896565985175711E-20, 2.011010841116287E-20, 2.0325975072194052E-20, 2.054421776431546E-20, 2.0764890334116344E-20, 2.0988048793463279E-20, 2.1213751442653714E-20, 2.144205900275679E-20, 2.1673034757993714E-20, 2.1906744709105117E-20, 2.2143257738760407E-20, 2.2382645790186164E-20, 2.2624984060329141E-20, 2.2870351209027197E-20, 2.3118829585841496E-20, 2.3370505476409172E-20, 2.3625469370411695E-20, 2.3883816253525597E-20, 2.4145645926034935E-20, 2.441106335114639E-20, 2.4680179036466914E-20, 2.4953109452591E-20, 2.522997749331276E-20, 2.5510912982642653E-20, 2.579605323458912E-20, 2.6085543672584608E-20, 2.6379538516522633E-20, 2.667820154666292E-20, 2.6981706955199747E-20, 2.7290240298129617E-20, 2.760399956226778E-20, 2.7923196364936903E-20, 2.8248057307096873E-20, 2.8578825504645344E-20, 2.8915762327478304E-20, 2.9259149381897149E-20, 2.9609290779395838E-20, 2.9966515744169338E-20, 3.0331181623398431E-20, 3.0703677379217488E-20, 3.108442766024987E-20, 3.1473897575051851E-20, 3.1872598321607185E-20, 3.228109386876898E-20, 3.2700008940944537E-20, 3.313003863165466E-20, 3.35719600725733E-20, 3.4026646723650843E-20, 3.4495086044066494E-20, 3.4978401579282243E-20, 3.5477880897439209E-20, 3.5995011394472534E-20, 3.6531526869552723E-20, 3.7089469133133434E-20, 3.767127106708655E-20, 3.8279871085571472E-20, 3.8918874931706407E-20, 3.9592791337014738E-20, 4.0307387768676638E-20, 4.1070251384909018E-20, 4.1891722989140241E-20, 4.2786564624839063E-20, 4.3777229834795032E-20, 4.490119402885342E-20, 4.6231235710575519E-20, 5.2372836842253038E-20, 5.4210108624275222E-20 };
private static readonly ulong[] ModifiedZigguratExponentialAliasBox = { 0xd6701309607c00fa, 0xb706635304f880f7, 0xe707a22b605bc0f4, 0xbb6690d043b840f2, 0xf7f0e20d4e6640ef, 0xf861ac912e35c0ed, 0xfc733e1db4d4c0eb, 0xea7c6f626bb600ed, 0xd6530fc94de280ef, 0xc606ebbf928680f0, 0xb88e65112bba80f2, 0xad3835f893548003, 0xa3894f25d92c00f4, 0x9b2975b2fba90002, 0x93d7ac60134580f5, 0x8d62f853f5de00f6, 0x87a5b38769c500f7, 0x82826dc545e100f7, 0x7de1c92824e08001, 0x79b0fa0dfb6380f8, 0x75e0b41799ab00f8, 0x72646198d58f80f9, 0x6f318ee74f9080f9, 0x6c3f7a79c27800fa, 0x6986bf9802fe00fa, 0x670114812b3700fa, 0x64a91705bf5e0000, 0x627a2400936a0000, 0x60703714cfdd0000, 0x5e87d0c784300000, 0x5cbde1897c370000, 0x5b0fb89f5ef300fc, 0x597af618aa1c00fc, 0x57fd7f361fa200fc, 0x569574c494a800fc, 0x55412b0c19c600fc, 0x53ff2307c87b00fb, 0x52ce04aae61d00fb, 0x51ac9a033e0500fb, 0x5099cb12df2700fb, 0x4f949a42401e00fb, 0x4e9c2151ca6c00fb, 0x4daf8eb65fa700fb, 0x4cce23501edc00fb, 0x4bf7306d94af00fb, 0x4b2a160fe10200fb, 0x4a6641665f8200fb, 0x49ab2b79c02900fb, 0x48f858000d6900fb, 0x484d5453dc2f00fb, 0x47a9b689eb5c00fb, 0x470d1ca1578900fb, 0x46772bcaa7a700fb, 0x45e78fc30cfc00fb, 0x455dfa4133f600fb, 0x44da22716a1d00fb, 0x445bc47f6c1d00fb, 0x43e2a12c2a3c00fb, 0x436e7d6e025600fb, 0x42ff221a76ac00fb, 0x42945b981cb100fb, 0x422df997f96c00fb, 0x41cbced566df00fb, 0x416db0dbf94f00fb, 0x411377d273ad00fb, 0x40bcfe4aa2cd00fb, 0x406a2115621200fb, 0x401abf1a58d000fb, 0x3fceb93355b300fb, 0x3f85f20a882200fb, 0x3f404dfbb77d00fb, 0x3efdb2f7ead100fb, 0x3ebe086b4eb700fb, 0x3e8137254cfd00fb, 0x3e472942828d00fb, 0x3e0fca184dd500fb, 0x3ddb062229e100fb, 0x3da8caf036b700fb, 0x3d790717509200fb, 0x3d4baa2226ba00fb, 0x3d20a4838a6e00fb, 0x3cf7e789b8f400fb, 0x3cd165528b6100fb, 0x3cad10c0885500fb, 0x3c8add70c62800fb, 0x3c6abfb16baf00fb, 0x3c4cac7900e800fb, 0x3c30995e3f4700fb, 0x3c167c90877b00fb, 0x3bfe4cd0d0fe00fb, 0x3be8016b2aa700fb, 0x3bd39230969400fb, 0x3bc0f7716a7100fb, 0x3bb029f812b300fb, 0x3ba123041a5f00fb, 0x3b93dc45b02800fb, 0x3b884fd96bd900fb, 0x3b7e78445a4900fb, 0x3b765070640100fb, 0x3b6fd3a8e46e00fb, 0x3b6afd979e3d00fb, 0x3b67ca41c8f100fb, 0x3b663605700e00fb, 0x3b663d970db500fb, 0x3b67ddff35dc00fb, 0x3b6b1498a5ce00fb, 0x3b6fdf0e54a800fb, 0x3b763b59c45000fb, 0x3b7e27c18cfc00fb, 0x3b87a2d7edea00fb, 0x3b92ab79a32a00fb, 0x3b9f40ccd1d000fb, 0x3bad6240241200fb, 0x3bbd0f8a075200fb, 0x3bce48a7fd3800fb, 0x3be10dde233400fb, 0x3bf55fb6da5c00fb, 0x3c0b3f02826000fb, 0x3c22acd7549400fb, 0x3c3baa9178a600fb, 0x3c5639d320a400fb, 0x3c725c84c0ea00fb, 0x3c9014d583c000fb, 0x3caf653bc01a00fb, 0x3cd05075a4ce00fb, 0x3cf2d989f35000fb, 0x3d1703c8e9ae00fb, 0x3d3cd2cd4b1000fb, 0x3d644a7d78bc00fb, 0x3d8d6f0cd39200fb, 0x3db844fd17c800fb, 0x3de4d11ff60a00fb, 0x3e131898d0e600fb, 0x3e4320de958c00fb, 0x3e74efbdccb200fb, 0x3ea88b5ac4ee00fb, 0x3eddfa33fd8600fb, 0x3f154324a35a00fb, 0x3f4e6d675ce600fb, 0x3f898099348600fb, 0x3fc684bcada400fb, 0x4005823d350200fb, 0x404681f2a8f400fb, 0x40898d25245600fb, 0x40cead91243000fb, 0x4115ed6bbef200fb, 0x415f576762da00fb, 0x41aaf6b8aa3000fb, 0x41f8d71b961c00fb, 0x424904d92de200fb, 0x429b8ccd588600fb, 0x42f07c6d351600fb, 0x4347e1cdc8fe00fb, 0x43a1cbab21f800fb, 0x43fe496feeae00fb, 0x445d6b3d957000fb, 0x44bf41f4c4e400fb, 0x4523df3ea8ca00fb, 0x458b5596ac7e00fb, 0x45f5b854e36400fb, 0x46631bb91c4e00fb, 0x46d394f6bcdc00fb, 0x47473a416e3c00fb, 0x47be22da94ec00fb, 0x4838671fba7000fb, 0x48b6209a002a00fb, 0x49376a0ea4be00fb, 0x49bc5f90997000fb, 0x4a451e938a8600fb, 0x4ad1c6000bb400fb, 0x4b627649699800fb, 0x4bf75184fef000fb, 0x4c907b834ce200fb, 0x4d2e19eaf70c00fb, 0x4dd05455bee000fb, 0x4e77546fd2dc00fb, 0x4f234619799800fb, 0x4fd4578b6d8a00fb, 0x508ab97e16ca00fb, 0x51469f540fb600fb, 0x52083f47eafe00fb, 0x52cfd29e0ba000fb, 0x539d95da7f5c00fb, 0x5471c8fb8e1a00fb, 0x554cafb93d9000fc, 0x562e91ca95aa00fc, 0x5717bb30e5b000fc, 0x58087c89f87800fc, 0x59012b69cb7000fc, 0x5a0222bca07000fc, 0x5b0bc3325da800fc, 0x5c1e73b4533000fc, 0x5d3aa1e671700000, 0x5e60c2b562e00000, 0x5f9152f31ac80000, 0x60ccd8035e740000, 0x6213e09a7b7c0000, 0x6367059046fc0000, 0x64c6eaca1ab40000, 0x6634403e83ec00fa, 0x67afc316678800fa, 0x693a3eef3f3c00fa, 0x6ad48f430f4800fa, 0x6c7fa0fb666400fa, 0x6e3c743590e800f9, 0x700c1e3edef400f9, 0x71efcbd1b44400f9, 0x73e8c39d0cac00f9, 0x75f86921c3c400f8, 0x78203ff3c96c00f8, 0x7a61ef6ee04000f8, 0x7cbf46f2a9740001, 0x7f3a42bc810c0001, 0x81d5117b70fc0001, 0x84921abea46c00f7, 0x87740667b2a400f7, 0x8a7dc550a9e000f6, 0x8db29b62b57c00f6, 0x91162b665fa400f5, 0x94ac84e89cb800f5, 0x987a34a5a9680002, 0x9c845807104800f4, 0xa0d0b466b73800f4, 0xa565d2f88b8400f3, 0xaa4b227a9e2800f3, 0xaf89201e3d480003, 0xb529898ac11400f2, 0xbb379a6f82dc00f1, 0xc1c058e4de2000f0, 0xc8d2f4f36d940004, 0xd081411ddbd800ef, 0xd8e04befb33800ee, 0xe20925a810280005, 0xec19e192689000ed, 0xf736e941b9c800ec, 0xfffffffffffcc0eb, 0xf84e32105543c006, 0xefc4d596acd3c0ec, 0xf04f089008e8c005, 0xe51b75a76848c0ee, 0xd9f6ddd4a2914004, 0xf4c648847f92c0f0, 0xc13cd60fe90ec0f1, 0xc2e92e87e11bc003, 0xb4ec402cc0ffc0f3, 0xa05bdee4bcda0002, 0xbf28b27170a080f5, 0xa64cb7a7255e80f6, 0xa99e6bb6d07f8001, 0xe7ae9dfdc03100f8, 0xf8596f91c66500f9, 0xd717854fe71f00fc, 0xd8d13e002ab80000, 0x00000000000000fb, 0x00000000000000fb, 0x00000000000000fb };
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static double ModifiedZigguratExponential<TRng>(TRng rng)
where TRng : notnull, RandomBase
{
const int TableSize = 256;
const int IMax = 252;
ulong u1 = rng.Next();
int tableX = (int)(u1 & (TableSize - 1));
if (tableX < IMax)
{
return u1 * ModifiedZigguratExponentialX[tableX];
}
[MethodImpl(MethodImplOptions.NoInlining)]
static double Fallback(TRng rng, ulong u1)
{
const ulong MaxConcaveRate = 0x17b3cab860ee7800;
[MethodImpl(MethodImplOptions.AggressiveInlining)]
static double SampleX(double x, int tableY) => Math.FusedMultiplyAdd(x, ModifiedZigguratExponentialX[tableY - 1] - ModifiedZigguratExponentialX[tableY], ModifiedZigguratExponentialX[tableY] * (2.0 * (1ul << 63)));
[MethodImpl(MethodImplOptions.AggressiveInlining)]
static double SampleY(double y, int tableY) => Math.FusedMultiplyAdd(y, ModifiedZigguratExponentialY[tableY - 1] - ModifiedZigguratExponentialY[tableY], ModifiedZigguratExponentialY[tableY] * (2.0 * (1ul << 63)));
double tail = 0;
while (true)
{
ulong u2 = rng.Next();
int tableAlias = (int)(u2 & (TableSize - 1));
int tableY = (u2 >> 8) < (ModifiedZigguratExponentialAliasBox[tableAlias] >> 8) ? tableAlias : (int)(ModifiedZigguratExponentialAliasBox[tableAlias] & (TableSize - 1));
if (tableY == 0)
{
// 最下層レイヤー: 尾から生成する
// 指数分布の尾は指数分布であり再帰的に適用可能なので、 tail を増やして再試行する
double x0 = ModifiedZigguratExponentialX[0] * (2.0 * (1ul << 63));
tail += x0;
u1 = rng.Next();
int tableX = (int)(u1 & (TableSize - 1));
if (tableX < IMax)
{
return tail + u1 * ModifiedZigguratExponentialX[tableX];
}
}
else
{
// へこみレイヤー
ulong ux = u1;
ulong uy = rng.Next();
double tx;
for (; ; ux = rng.Next(), uy = rng.Next())
{
// 三角形の右上; 上下左右反転させる
if (uy < ux)
{
ux = ~ux;
uy = ~uy;
}
// 三角形の左下; 採択
if (uy >= ux + MaxConcaveRate)
{
tx = SampleX(ux, tableY);
break;
}
// 確率密度関数より下なら採択
tx = SampleX(ux, tableY);
double ty = SampleY(uy, tableY);
if (ty <= Math.Exp(-tx))
{
break;
}
}
return tail + tx;
}
}
}
return Fallback(rng, u1);
}
}
}
namespace RngLab
{
public static class Program
{
public static void Main(string[] args)
{
var summary = BenchmarkRunner.Run<Ziggurat>();
}
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment