前回記事の続き。
【C#】正規乱数 (1):Box-Mullar 法, 中心極限定理を利用した方法 - てくメモ
逆関数法
近似によって求まった、正規分布の累積分布関数の逆関数を用いる手法がある。
そのひとつが以下のリンク。
An algorithm for computing the inverse normal cumulative distribution function(※ リンク切れ)
└ Wayback Machine アーカイブ
C# に書き下したものは以下。
大量のマジックナンバーがいかにも逆という空気を感じさせる。
using static System.Math; // ※ コード中のコメントは参考元のもの public static double InvNorm(double p) { const double a1 = -3.969683028665376e+01; const double a2 = 2.209460984245205e+02; const double a3 = -2.759285104469687e+02; const double a4 = 1.383577518672690e+02; const double a5 = -3.066479806614716e+01; const double a6 = 2.506628277459239e+00; const double b1 = -5.447609879822406e+01; const double b2 = 1.615858368580409e+02; const double b3 = -1.556989798598866e+02; const double b4 = 6.680131188771972e+01; const double b5 = -1.328068155288572e+01; const double c1 = -7.784894002430293e-03; const double c2 = -3.223964580411365e-01; const double c3 = -2.400758277161838e+00; const double c4 = -2.549732539343734e+00; const double c5 = 4.374664141464968e+00; const double c6 = 2.938163982698783e+00; const double d1 = 7.784695709041462e-03; const double d2 = 3.224671290700398e-01; const double d3 = 2.445134137142996e+00; const double d4 = 3.754408661907416e+00; const double pLow = 0.02425; const double pHigh = 1 - pLow; if (p < pLow) { // Rational approximation for lower region. var q = Sqrt(-2 * Log(p)); return (((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6) / ((((d1 * q + d2) * q + d3) * q + d4) * q + 1); } else if (pHigh < p) { // Rational approximation for upper region: var q = Sqrt(-2 * Log(1 - p)); return -(((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6) / ((((d1 * q + d2) * q + d3) * q + d4) * q + 1); } else { // Rational approximation for central region: var q = p - 0.5; var r = q * q; return (((((a1 * r + a2) * r + a3) * r + a4) * r + a5) * r + a6) * q / (((((b1 * r + b2) * r + b3) * r + b4) * r + b5) * r + 1); } }
累積分布関数の逆関数なので、確率変数 p を与えることでパーセンタイルを得ることができる。
このことについて、例として定番ないわゆる受験偏差値を考える。
試しに偏差値75とWeb検索したところ上位0.6%と出てきたので、p として0.994を与えて、境界値が75になっているか確認してみる。偏差値は平均50・標準偏差10になるよう調整されているので、それを考慮する。
Console.WriteLine(RandomUtil.InvNorm(0.994) * 10 + 50); // 75.12144330576888
さて、例では p を指定して与えたけれど、これが [0, 1) の一様乱数ならどうなるか。
正規乱数となる。
Ziggurat 法
事前計算で用意したルックアップテーブルを用いる非常に高速な手法。
前回冒頭でも提示した次の記事において、アルゴリズムの内容・実装が詳述されている。
詳説 Ziggurat 法 ~ 正規分布・指数分布乱数の高速生成 - 屋根裏工房改
アルゴリズムの内容そのものはリンク先参照として、高速たるゆえんが分かる部分を引用。
(※ 原コードではメソッド中にローカル関数としてFallback
の実装があるが省略)
// ※ コード中のコメントは参考元のもの [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]; } // 以降の処理は相対的に重いため別関数に分け、上の高速パスのインライン化を試みる return Fallback(rng, u1); }
リンク先の解説によれば、約 98.5% はFallback
前に値を返すということ。
その 98.5% では、ルックアップテーブルから値を引いて一様乱数と掛けるだけ。
いかにも速い……
比較
N = 65536 でループ。
- Box-Mullar 法については一度に生成する両方を使った
- Ziggurat 法は引用リンク先の改良法を利用したが、用いる一様乱数については他との並びで置き換えた
- Ziggurat 法で事前に格納する値については、引用リンク先のものをそのまま使用した
比較用コード折りたたみ(注:非常に長い)
[Params(1 << 16)] public int N = 1 << 16; private NormProvider norm = new(); private sealed class NormProvider { private readonly Random rand = new(); private readonly Koku koku = new(5); private sealed class Koku { private readonly double a, sd; private readonly int ni; public double N { get => ni + a; init { if (value < 1d) throw new ArgumentOutOfRangeException("N < 1"); double n = value; ni = (int)n; a = n - ni; var sigma = (ni / 3d) + (a * a / 3d); sd = Sqrt(sigma); } } public Koku(int n) => N = n; public double Next(Random rand) { double s = 0d; for (int i = 0; i < ni; i++) { s += rand.NextStandardUniform(); } s += rand.NextStandardUniform() * a; return s / sd; } } public NormProvider() { } public (double Z1, double Z2) NextBoxMullar() { double x = rand.NextDouble(); double y = rand.NextDouble(); double r = Sqrt(-2 * Log(1 - x)); (double sin, double cos) = SinCos(y * (PI * 2)); return (r * sin, r * cos); } public double NextRand12() { return rand.NextDouble() + rand.NextDouble() + rand.NextDouble() + rand.NextDouble() + rand.NextDouble() + rand.NextDouble() + rand.NextDouble() + rand.NextDouble() + rand.NextDouble() + rand.NextDouble() + rand.NextDouble() + rand.NextDouble() - 6d; } public double NextKoku() => koku.Next(rand); public double NextInvNorm() => InvNorm(rand.NextDouble()); private static double InvNorm(double p) { const double a1 = -3.969683028665376e+01; const double a2 = 2.209460984245205e+02; const double a3 = -2.759285104469687e+02; const double a4 = 1.383577518672690e+02; const double a5 = -3.066479806614716e+01; const double a6 = 2.506628277459239e+00; const double b1 = -5.447609879822406e+01; const double b2 = 1.615858368580409e+02; const double b3 = -1.556989798598866e+02; const double b4 = 6.680131188771972e+01; const double b5 = -1.328068155288572e+01; const double c1 = -7.784894002430293e-03; const double c2 = -3.223964580411365e-01; const double c3 = -2.400758277161838e+00; const double c4 = -2.549732539343734e+00; const double c5 = 4.374664141464968e+00; const double c6 = 2.938163982698783e+00; const double d1 = 7.784695709041462e-03; const double d2 = 3.224671290700398e-01; const double d3 = 2.445134137142996e+00; const double d4 = 3.754408661907416e+00; const double pLow = 0.02425; const double pHigh = 1 - pLow; if (p < pLow) { // Rational approximation for lower region. var q = Sqrt(-2 * Log(p)); return (((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6) / ((((d1 * q + d2) * q + d3) * q + d4) * q + 1); } else if (pHigh < p) { // Rational approximation for upper region: var q = Sqrt(-2 * Log(1 - p)); return -(((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6) / ((((d1 * q + d2) * q + d3) * q + d4) * q + 1); } else { // Rational approximation for central region: var q = p - 0.5; var r = q * q; return (((((a1 * r + a2) * r + a3) * r + a4) * r + a5) * r + a6) * q / (((((b1 * r + b2) * r + b3) * r + b4) * r + b5) * r + 1); } } public double NextZiggurat() => ModifiedZigguratNormal(rand); 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(Random rng) { const int TableSize = 256; // 四角形のレイヤー数 const int IMax = 253; long u1 = (long)rng.NextUInt64(); int tableX = (int)(u1 & (TableSize - 1)); // 四角形のレイヤーなら x 座標を返す if (tableX < IMax) { return u1 * ModifiedZigguratNormalX[tableX]; } // 以降の処理は相対的に重いため別関数に分け、上の高速パスのインライン化を試みる [MethodImpl(MethodImplOptions.NoInlining)] static double Fallback(Random 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.NextUInt64(); 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.NextUInt64() >> 1; double tx; for (; ; ux = rng.NextUInt64() >> 1, uy = rng.NextUInt64() >> 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.NextUInt64() >> 1; double tx; for (; ; ux = rng.NextUInt64() >> 1, uy = rng.NextUInt64() >> 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.NextUInt64() >> 1; double tx; for (; ; ux = rng.NextUInt64() >> 1, uy = rng.NextUInt64() >> 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 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(Random rng) { const int TableSize = 256; const int IMax = 252; ulong u1 = rng.NextUInt64(); int tableX = (int)(u1 & (TableSize - 1)); if (tableX < IMax) { return u1 * ModifiedZigguratExponentialX[tableX]; } [MethodImpl(MethodImplOptions.NoInlining)] static double Fallback(Random 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.NextUInt64(); 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.NextUInt64(); int tableX = (int)(u1 & (TableSize - 1)); if (tableX < IMax) { return tail + u1 * ModifiedZigguratExponentialX[tableX]; } } else { // へこみレイヤー ulong ux = u1; ulong uy = rng.NextUInt64(); double tx; for (; ; ux = rng.NextUInt64(), uy = rng.NextUInt64()) { // 三角形の右上; 上下左右反転させる 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); } } [Benchmark] public double BoxMullar() { double result = 0; for (int i = 0; i < N; i += 2) { var (r1, r2) = norm.NextBoxMullar(); result += r1; result += r2; } return result; } [Benchmark] public double Rand12() { double result = 0; for (int i = 0; i < N; i++) { result += norm.NextRand12(); } return result; } [Benchmark] public double Koku() { double result = 0; for (int i = 0; i < N; i++) { result += norm.NextKoku(); } return result; } [Benchmark] public double InvNorm() { double result = 0; for (int i = 0; i < N; i++) { result += norm.NextInvNorm(); } return result; } [Benchmark] public double Ziggurat() { double result = 0; for (int i = 0; i < N; i++) { result += norm.NextZiggurat(); } return result; }
- 注
NextUInt64
はulong
の一様乱数を生成する別途の拡張メソッドNextStandardUniform
は[-1.0 - 1.0) の一様乱数を生成する別途の拡張メソッド
Method | N | Mean | Error | StdDev | ---------- |------ |-----------:|-----------:|---------:| BoxMullar | 65536 | 1,713.4 μs | 340.3 μs | 18.65 μs | Rand12 | 65536 | 4,265.1 μs | 758.7 μs | 41.58 μs | Koku | 65536 | 1,589.1 μs | 1,045.9 μs | 57.33 μs | InvNorm | 65536 | 993.3 μs | 376.2 μs | 20.62 μs | Ziggurat | 65536 | 390.8 μs | 211.7 μs | 11.60 μs |
一様乱数12個の手法はそれなりに他から後れ、もうひとつの簡便な手法であるコクのある乱数も Box-Mullar 法に対し微差に速い程度で、一様乱数を何度も使うのはそれなりに負担となるのが感じられる。
そして、Ziggurat 法はやはり速い。
数万回で数百マイクロ秒からミリ秒の単位なので常に速い方法が要求されるわけではないようにみえるけれど、ここに速さを求める場合には際立って速い Ziggurat 法は有力だと感じた。