diff --git a/base/librandom.jl b/base/librandom.jl index 4e169026d448b..d7d07488693e8 100644 --- a/base/librandom.jl +++ b/base/librandom.jl @@ -61,7 +61,7 @@ function dsfmt_init_by_array(s::DSFMT_state, seed::Vector{Uint32}) end function dsfmt_gv_init_by_array(seed::Vector{Uint32}) - ccall((:dsfmt_gv_init_by_array,:librandom), + ccall((:dsfmt_gv_init_by_array,:librandom), Void, (Ptr{Uint32}, Int32), seed, length(seed)) @@ -136,7 +136,6 @@ for (genrand_fill, gv_genrand_fill, genrand_fill_name, gv_genrand_fill_name) in end return A end - end end @@ -153,51 +152,541 @@ function dsfmt_gv_genrand_uint32() ()) end -## randmtzig - -function randmtzig_randn(s::DSFMT_state) - ccall((:randmtzig_randn,:librandom), - Float64, - (Ptr{Void},), - s.val) -end +# randmtzig +## Tables for normal variates +const ki = + Uint64[0x0007799ec012f7b3, 0,0x0006045f4c7de363,0x0006d1aa7d5ec0a6, + 0x000728fb3f60f778,0x0007592af4e9fbc0,0x000777a5c0bf655d,0x00078ca3857d2256, + 0x00079bf6b0ffe58c,0x0007a7a34ab092ae,0x0007b0d2f20dd1cb,0x0007b83d3aa9cb52, + 0x0007be597614224e,0x0007c3788631abea,0x0007c7d32bc192ef,0x0007cb9263a6e86d, + 0x0007ced483edfa84,0x0007d1b07ac0fd39,0x0007d437ef2da5fc,0x0007d678b069aa6e, + 0x0007d87db38c5c87,0x0007da4fc6a9ba63,0x0007dbf611b37f3c,0x0007dd7674d0f286, + 0x0007ded5ce8205f7,0x0007e018307fb62c,0x0007e141081bd124,0x0007e2533d712de9, + 0x0007e3514bbd7718,0x0007e43d54944b52,0x0007e5192f25ef43,0x0007e5e67481118d, + 0x0007e6a6897c1ce2,0x0007e75aa6c7f64d,0x0007e803df8ee498,0x0007e8a326eb6273, + 0x0007e93954717a28,0x0007e9c727f8648f,0x0007ea4d4cc85a3d,0x0007eacc5c4907a9, + 0x0007eb44e0474cf6,0x0007ebb754e47419,0x0007ec242a3d8474,0x0007ec8bc5d69646, + 0x0007ecee83d3d6e9,0x0007ed4cb8082f46,0x0007eda6aee01710,0x0007edfcae2dfe68, + 0x0007ee4ef5dccd3f,0x0007ee9dc08c394f,0x0007eee9441a17c8,0x0007ef31b21b4fb2, + 0x0007ef773846a8a7,0x0007efba00d35a17,0x0007effa32ccf69f,0x0007f037f25e1279, + 0x0007f0736112d12c,0x0007f0ac9e145c26,0x0007f0e3c65e1fcc,0x0007f118f4ed8e55, + 0x0007f14c42ed0dc9,0x0007f17dc7daa0c3,0x0007f1ad99aac6a5,0x0007f1dbcce80015, + 0x0007f20874cf56bf,0x0007f233a36a3b9a,0x0007f25d69a604ae,0x0007f285d7694a92, + 0x0007f2acfba75e3c,0x0007f2d2e4720909,0x0007f2f79f09c345,0x0007f31b37ec883c, + 0x0007f33dbae36abc,0x0007f35f330f08d5,0x0007f37faaf2fa79,0x0007f39f2c805381, + 0x0007f3bdc11f4f1d,0x0007f3db71b83851,0x0007f3f846bba121,0x0007f4144829f846, + 0x0007f42f7d9a8b9e,0x0007f449ee420432,0x0007f463a0f8675e,0x0007f47c9c3ea77b, + 0x0007f494e643cd8e,0x0007f4ac84e9c475,0x0007f4c37dc9cd51,0x0007f4d9d638a433, + 0x0007f4ef934a5b6b,0x0007f504b9d5f33e,0x0007f5194e78b352,0x0007f52d55994a96, + 0x0007f540d36aba0d,0x0007f553cbef0e78,0x0007f56642f9ec90,0x0007f5783c32f31e, + 0x0007f589bb17f609,0x0007f59ac2ff1525,0x0007f5ab5718b15a,0x0007f5bb7a71427d, + 0x0007f5cb2ff3100a,0x0007f5da7a67cebe,0x0007f5e95c7a24e8,0x0007f5f7d8b7171f, + 0x0007f605f18f5ef4,0x0007f613a958ad0a,0x0007f621024ed7ea,0x0007f62dfe94f8cb, + 0x0007f63aa036777a,0x0007f646e928065a,0x0007f652db488f88,0x0007f65e786213ff, + 0x0007f669c22a7d8b,0x0007f674ba44645a,0x0007f67f623fc8dc,0x0007f689bb9ac294, + 0x0007f693c7c22482,0x0007f69d881217a6,0x0007f6a6fdd6ac37,0x0007f6b02a4c61ee, + 0x0007f6b90ea0a7f4,0x0007f6c1abf254c1,0x0007f6ca03521664,0x0007f6d215c2db82, + 0x0007f6d9e43a355a,0x0007f6e16fa0b329,0x0007f6e8b8d23729,0x0007f6efc09e4569, + 0x0007f6f687c84cc0,0x0007f6fd0f07ea0a,0x0007f703570925e2,0x0007f709606cad03, + 0x0007f70f2bc80370,0x0007f714b9a5b292,0x0007f71a0a85725e,0x0007f71f1edc4d9e, + 0x0007f723f714c17a,0x0007f728938ed843,0x0007f72cf4a03fa1,0x0007f7311a945a17, + 0x0007f73505ac4bf8,0x0007f738b61f03bd,0x0007f73c2c193dc1,0x0007f73f67bd835d, + 0x0007f74269242559,0x0007f745305b31a1,0x0007f747bd666429,0x0007f74a103f12ed, + 0x0007f74c28d414f6,0x0007f74e0709a42e,0x0007f74faab939fa,0x0007f75113b16657, + 0x0007f75241b5a156,0x0007f753347e16b9,0x0007f753ebb76b7d,0x0007f75467027d05, + 0x0007f754a5f4199e,0x0007f754a814b208,0x0007f7546ce003af,0x0007f753f3c4bb2a, + 0x0007f7533c240e92,0x0007f75245514f41,0x0007f7510e91726d,0x0007f74f971a9012, + 0x0007f74dde135797,0x0007f74be2927972,0x0007f749a39e051d,0x0007f747202aba8b, + 0x0007f744571b4e3d,0x0007f741473f9efe,0x0007f73def53dc44,0x0007f73a4dff9c00, + 0x0007f73661d4deaf,0x0007f732294f0040,0x0007f72da2d19444,0x0007f728cca72bdb, + 0x0007f723a5000367,0x0007f71e29f09628,0x0007f7185970156c,0x0007f7123156c102, + 0x0007f70baf5c1e2d,0x0007f704d1150a24,0x0007f6fd93f1a4e6,0x0007f6f5f53b10b6, + 0x0007f6edf211023f,0x0007f6e587671cea,0x0007f6dcb202167a,0x0007f6d36e749c65, + 0x0007f6c9b91bf4c6,0x0007f6bf8e1c541b,0x0007f6b4e95ce015,0x0007f6a9c6835700, + 0x0007f69e20ef5211,0x0007f691f3b517ec,0x0007f6853997f322,0x0007f677ed03ff1a, + 0x0007f66a08075bdc,0x0007f65b844ab75a,0x0007f64c5b091860,0x0007f63c8506d4bc, + 0x0007f62bfa8798ff,0x0007f61ab34364b1,0x0007f608a65a599a,0x0007f5f5ca4737e8, + 0x0007f5e214d05b49,0x0007f5cd7af7066e,0x0007f5b7f0e4c2a1,0x0007f5a169d68fcf, + 0x0007f589d80596a6,0x0007f5712c8d0174,0x0007f557574c912b,0x0007f53c46c77193, + 0x0007f51fe7feb9f3,0x0007f5022646ecfb,0x0007f4e2eb17ab1d,0x0007f4c21dd4a3d2, + 0x0007f49fa38ea394,0x0007f47b5ebb62eb,0x0007f4552ee27473,0x0007f42cf03d58f5, + 0x0007f4027b4854a0,0x0007f3d5a44119e0,0x0007f3a63a8fb553,0x0007f37408155101, + 0x0007f33ed05b55ec,0x0007f3064f9c183f,0x0007f2ca399c7ba1,0x0007f28a384bb940, + 0x0007f245ea1b7a2c,0x0007f1fcdffe8f1c,0x0007f1ae9af758ce,0x0007f15a8917f27f, + 0x0007f10001ccaaac,0x0007f09e413c418a,0x0007f034627733d8,0x0007efc15815b8d5, + 0x0007ef43e2bf7f55,0x0007eeba84e31dff,0x0007ee237294df89,0x0007ed7c7c170141, + 0x0007ecc2f0d95d3b,0x0007ebf377a46782,0x0007eb09d6deb286,0x0007ea00a4f17809, + 0x0007e8d0d3da63d6,0x0007e771023b0fcf,0x0007e5d46c2f08d9,0x0007e3e937669691, + 0x0007e195978f1176,0x0007deb2c0e05c1d,0x0007db0362002a1a,0x0007d6202c15143a, + 0x0007cf4b8f00a2cc,0x0007c4fd24520efe,0x0007b362fbf81816,0x00078d2d25998e25] +const wi = + [1.7367254121602630e-15,9.5586603514556339e-17,1.2708704834810623e-16, + 1.4909740962495474e-16,1.6658733631586268e-16,1.8136120810119029e-16, + 1.9429720153135588e-16,2.0589500628482093e-16,2.1646860576895422e-16, + 2.2622940392218116e-16,2.3532718914045892e-16,2.4387234557428771e-16, + 2.5194879829274225e-16,2.5962199772528103e-16,2.6694407473648285e-16, + 2.7395729685142446e-16,2.8069646002484804e-16,2.8719058904113930e-16, + 2.9346417484728883e-16,2.9953809336782113e-16,3.0543030007192440e-16, + 3.1115636338921572e-16,3.1672988018581815e-16,3.2216280350549905e-16, + 3.2746570407939751e-16,3.3264798116841710e-16,3.3771803417353232e-16, + 3.4268340353119356e-16,3.4755088731729758e-16,3.5232663846002031e-16, + 3.5701624633953494e-16,3.6162480571598339e-16,3.6615697529653540e-16, + 3.7061702777236077e-16,3.7500889278747798e-16,3.7933619401549554e-16, + 3.8360228129677279e-16,3.8781025861250247e-16,3.9196300853257678e-16, + 3.9606321366256378e-16,4.0011337552546690e-16,4.0411583124143332e-16, + 4.0807276830960448e-16,4.1198623774807442e-16,4.1585816580828064e-16, + 4.1969036444740733e-16,4.2348454071520708e-16,4.2724230518899761e-16, + 4.3096517957162941e-16,4.3465460355128760e-16,4.3831194100854571e-16, + 4.4193848564470665e-16,4.4553546609579137e-16,4.4910405058828750e-16, + 4.5264535118571397e-16,4.5616042766900381e-16,4.5965029108849407e-16, + 4.6311590702081647e-16,4.6655819856008752e-16,4.6997804906941950e-16, + 4.7337630471583237e-16,4.7675377680908526e-16,4.8011124396270155e-16, + 4.8344945409350080e-16,4.8676912627422087e-16,4.9007095245229938e-16, + 4.9335559904654139e-16,4.9662370843221783e-16,4.9987590032409088e-16, + 5.0311277306593187e-16,5.0633490483427195e-16,5.0954285476338923e-16, + 5.1273716399787966e-16,5.1591835667857364e-16,5.1908694086703434e-16, + 5.2224340941340417e-16,5.2538824077194543e-16,5.2852189976823820e-16, + 5.3164483832166176e-16,5.3475749612647295e-16,5.3786030129452348e-16, + 5.4095367096239933e-16,5.4403801186554671e-16,5.4711372088173611e-16, + 5.5018118554603362e-16,5.5324078453927836e-16,5.5629288815190902e-16, + 5.5933785872484621e-16,5.6237605106900435e-16,5.6540781286489604e-16, + 5.6843348504368141e-16,5.7145340215092040e-16,5.7446789269419609e-16, + 5.7747727947569648e-16,5.8048187991076857e-16,5.8348200633338921e-16, + 5.8647796628943653e-16,5.8947006281858718e-16,5.9245859472561339e-16, + 5.9544385684180598e-16,5.9842614027720281e-16,6.0140573266426640e-16, + 6.0438291839361250e-16,6.0735797884236057e-16,6.1033119259564394e-16, + 6.1330283566179110e-16,6.1627318168165963e-16,6.1924250213258470e-16, + 6.2221106652737879e-16,6.2517914260879998e-16,6.2814699653988953e-16, + 6.3111489309056042e-16,6.3408309582080600e-16,6.3705186726088149e-16, + 6.4002146908880247e-16,6.4299216230548961e-16,6.4596420740788321e-16, + 6.4893786456033965e-16,6.5191339376461587e-16,6.5489105502874154e-16, + 6.5787110853507413e-16,6.6085381480782587e-16,6.6383943488035057e-16, + 6.6682823046247459e-16,6.6982046410815579e-16,6.7281639938375311e-16, + 6.7581630103719006e-16,6.7882043516829803e-16,6.8182906940062540e-16, + 6.8484247305500383e-16,6.8786091732516637e-16,6.9088467545571690e-16, + 6.9391402292275690e-16,6.9694923761748294e-16,6.9999060003307640e-16, + 7.0303839345521508e-16,7.0609290415654822e-16,7.0915442159548734e-16, + 7.1222323861967788e-16,7.1529965167453030e-16,7.1838396101720629e-16, + 7.2147647093647067e-16,7.2457748997883870e-16,7.2768733118146927e-16, + 7.3080631231227429e-16,7.3393475611774048e-16,7.3707299057898310e-16, + 7.4022134917657997e-16,7.4338017116476479e-16,7.4654980185558890e-16, + 7.4973059291369793e-16,7.5292290266240584e-16,7.5612709640179217e-16, + 7.5934354673958895e-16,7.6257263393567558e-16,7.6581474626104873e-16, + 7.6907028037219191e-16,7.7233964170182985e-16,7.7562324486711744e-16, + 7.7892151409638524e-16,7.8223488367564108e-16,7.8556379841610841e-16, + 7.8890871414417552e-16,7.9227009821522709e-16,7.9564843005293662e-16, + 7.9904420171571300e-16,8.0245791849212591e-16,8.0589009952726568e-16, + 8.0934127848215009e-16,8.1281200422845008e-16,8.1630284158098775e-16, + 8.1981437207065329e-16,8.2334719476060504e-16,8.2690192710884700e-16, + 8.3047920588053737e-16,8.3407968811366288e-16,8.3770405214202216e-16, + 8.4135299867980282e-16,8.4502725197240968e-16,8.4872756101861549e-16, + 8.5245470086955962e-16,8.5620947401062333e-16,8.5999271183276646e-16, + 8.6380527620052589e-16,8.6764806112455816e-16,8.7152199454736980e-16, + 8.7542804025171749e-16,8.7936719990210427e-16,8.8334051523084080e-16, + 8.8734907038131345e-16,8.9139399442240861e-16,8.9547646404950677e-16, + 8.9959770648910994e-16,9.0375900262601175e-16,9.0796169037400680e-16, + 9.1220716831348461e-16,9.1649689962191353e-16,9.2083241632623076e-16, + 9.2521532390956933e-16,9.2964730630864167e-16,9.3413013134252651e-16, + 9.3866565661866598e-16,9.4325583596767065e-16,9.4790272646517382e-16, + 9.5260849610662787e-16,9.5737543220974496e-16,9.6220595062948384e-16, + 9.6710260588230542e-16,9.7206810229016259e-16,9.7710530627072088e-16, + 9.8221725991905411e-16,9.8740719604806711e-16,9.9267855488079765e-16, + 9.9803500261836449e-16,1.0034804521436181e-15,1.0090190861637457e-15, + 1.0146553831467086e-15,1.0203941464683124e-15,1.0262405372613567e-15, + 1.0322001115486456e-15,1.0382788623515399e-15,1.0444832676000471e-15, + 1.0508203448355195e-15,1.0572977139009890e-15,1.0639236690676801e-15, + 1.0707072623632994e-15,1.0776584002668106e-15,1.0847879564403425e-15, + 1.0921079038149563e-15,1.0996314701785628e-15,1.1073733224935752e-15, + 1.1153497865853155e-15,1.1235791107110833e-15,1.1320817840164846e-15, + 1.1408809242582780e-15,1.1500027537839792e-15,1.1594771891449189e-15, + 1.1693385786910960e-15,1.1796266352955801e-15,1.1903876299282890e-15, + 1.2016759392543819e-15,1.2135560818666897e-15,1.2261054417450561e-15, + 1.2394179789163251e-15,1.2536093926602567e-15,1.2688244814255010e-15, + 1.2852479319096109e-15,1.3031206634689985e-15,1.3227655770195326e-15, + 1.3446300925011171e-15,1.3693606835128518e-15,1.3979436672775240e-15, + 1.4319989869661328e-15,1.4744848603597596e-15,1.5317872741611144e-15, + 1.6227698675312968e-15] +const fi = + [1.0000000000000000e+00,9.7710170126767082e-01,9.5987909180010600e-01, + 9.4519895344229909e-01,9.3206007595922991e-01,9.1999150503934646e-01, + 9.0872644005213032e-01,8.9809592189834297e-01,8.8798466075583282e-01, + 8.7830965580891684e-01,8.6900868803685649e-01,8.6003362119633109e-01, + 8.5134625845867751e-01,8.4291565311220373e-01,8.3471629298688299e-01, + 8.2672683394622093e-01,8.1892919160370192e-01,8.1130787431265572e-01, + 8.0384948317096383e-01,7.9654233042295841e-01,7.8937614356602404e-01, + 7.8234183265480195e-01,7.7543130498118662e-01,7.6863731579848571e-01, + 7.6195334683679483e-01,7.5537350650709567e-01,7.4889244721915638e-01, + 7.4250529634015061e-01,7.3620759812686210e-01,7.2999526456147568e-01, + 7.2386453346862967e-01,7.1781193263072152e-01,7.1183424887824798e-01, + 7.0592850133275376e-01,7.0009191813651117e-01,6.9432191612611627e-01, + 6.8861608300467136e-01,6.8297216164499430e-01,6.7738803621877308e-01, + 6.7186171989708166e-01,6.6639134390874977e-01,6.6097514777666277e-01, + 6.5561147057969693e-01,6.5029874311081637e-01,6.4503548082082196e-01, + 6.3982027745305614e-01,6.3465179928762327e-01,6.2952877992483625e-01, + 6.2445001554702606e-01,6.1941436060583399e-01,6.1442072388891344e-01, + 6.0946806492577310e-01,6.0455539069746733e-01,5.9968175261912482e-01, + 5.9484624376798689e-01,5.9004799633282545e-01,5.8528617926337090e-01, + 5.8055999610079034e-01,5.7586868297235316e-01,5.7121150673525267e-01, + 5.6658776325616389e-01,5.6199677581452390e-01,5.5743789361876550e-01, + 5.5291049042583185e-01,5.4841396325526537e-01,5.4394773119002582e-01, + 5.3951123425695158e-01,5.3510393238045717e-01,5.3072530440366150e-01, + 5.2637484717168403e-01,5.2205207467232140e-01,5.1775651722975591e-01, + 5.1348772074732651e-01,5.0924524599574761e-01,5.0502866794346790e-01, + 5.0083757512614835e-01,4.9667156905248933e-01,4.9253026364386815e-01, + 4.8841328470545758e-01,4.8432026942668288e-01,4.8025086590904642e-01, + 4.7620473271950547e-01,4.7218153846772976e-01,4.6818096140569321e-01, + 4.6420268904817391e-01,4.6024641781284248e-01,4.5631185267871610e-01, + 4.5239870686184824e-01,4.4850670150720273e-01,4.4463556539573912e-01, + 4.4078503466580377e-01,4.3695485254798533e-01,4.3314476911265209e-01, + 4.2935454102944126e-01,4.2558393133802180e-01,4.2183270922949573e-01, + 4.1810064983784795e-01,4.1438753404089090e-01,4.1069314827018799e-01, + 4.0701728432947315e-01,4.0335973922111429e-01,3.9972031498019700e-01, + 3.9609881851583223e-01,3.9249506145931540e-01,3.8890886001878855e-01, + 3.8534003484007706e-01,3.8178841087339344e-01,3.7825381724561896e-01, + 3.7473608713789086e-01,3.7123505766823922e-01,3.6775056977903225e-01, + 3.6428246812900372e-01,3.6083060098964775e-01,3.5739482014578022e-01, + 3.5397498080007656e-01,3.5057094148140588e-01,3.4718256395679348e-01, + 3.4380971314685055e-01,3.4045225704452164e-01,3.3711006663700588e-01, + 3.3378301583071823e-01,3.3047098137916342e-01,3.2717384281360129e-01, + 3.2389148237639104e-01,3.2062378495690530e-01,3.1737063802991350e-01, + 3.1413193159633707e-01,3.1090755812628634e-01,3.0769741250429189e-01, + 3.0450139197664983e-01,3.0131939610080288e-01,2.9815132669668531e-01, + 2.9499708779996164e-01,2.9185658561709499e-01,2.8872972848218270e-01, + 2.8561642681550159e-01,2.8251659308370741e-01,2.7943014176163772e-01, + 2.7635698929566810e-01,2.7329705406857691e-01,2.7025025636587519e-01, + 2.6721651834356114e-01,2.6419576399726080e-01,2.6118791913272082e-01, + 2.5819291133761890e-01,2.5521066995466168e-01,2.5224112605594190e-01, + 2.4928421241852824e-01,2.4633986350126363e-01,2.4340801542275012e-01, + 2.4048860594050039e-01,2.3758157443123795e-01,2.3468686187232990e-01, + 2.3180441082433859e-01,2.2893416541468023e-01,2.2607607132238020e-01, + 2.2323007576391746e-01,2.2039612748015194e-01,2.1757417672433113e-01, + 2.1476417525117358e-01,2.1196607630703015e-01,2.0917983462112499e-01, + 2.0640540639788071e-01,2.0364274931033485e-01,2.0089182249465656e-01, + 1.9815258654577511e-01,1.9542500351413428e-01,1.9270903690358912e-01, + 1.9000465167046496e-01,1.8731181422380025e-01,1.8463049242679927e-01, + 1.8196065559952254e-01,1.7930227452284767e-01,1.7665532144373500e-01, + 1.7401977008183875e-01,1.7139559563750595e-01,1.6878277480121151e-01, + 1.6618128576448205e-01,1.6359110823236570e-01,1.6101222343751107e-01, + 1.5844461415592431e-01,1.5588826472447920e-01,1.5334316106026283e-01, + 1.5080929068184568e-01,1.4828664273257453e-01,1.4577520800599403e-01, + 1.4327497897351341e-01,1.4078594981444470e-01,1.3830811644855071e-01, + 1.3584147657125373e-01,1.3338602969166913e-01,1.3094177717364430e-01, + 1.2850872227999952e-01,1.2608687022018586e-01,1.2367622820159654e-01, + 1.2127680548479021e-01,1.1888861344290998e-01,1.1651166562561080e-01, + 1.1414597782783835e-01,1.1179156816383801e-01,1.0944845714681163e-01, + 1.0711666777468364e-01,1.0479622562248690e-01,1.0248715894193508e-01, + 1.0018949876880981e-01,9.7903279038862284e-02,9.5628536713008819e-02, + 9.3365311912690860e-02,9.1113648066373634e-02,8.8873592068275789e-02, + 8.6645194450557961e-02,8.4428509570353374e-02,8.2223595813202863e-02, + 8.0030515814663056e-02,7.7849336702096039e-02,7.5680130358927067e-02, + 7.3522973713981268e-02,7.1377949058890375e-02,6.9245144397006769e-02, + 6.7124653827788497e-02,6.5016577971242842e-02,6.2921024437758113e-02, + 6.0838108349539864e-02,5.8767952920933758e-02,5.6710690106202902e-02, + 5.4666461324888914e-02,5.2635418276792176e-02,5.0617723860947761e-02, + 4.8613553215868521e-02,4.6623094901930368e-02,4.4646552251294443e-02, + 4.2684144916474431e-02,4.0736110655940933e-02,3.8802707404526113e-02, + 3.6884215688567284e-02,3.4980941461716084e-02,3.3093219458578522e-02, + 3.1221417191920245e-02,2.9365939758133314e-02,2.7527235669603082e-02, + 2.5705804008548896e-02,2.3902203305795882e-02,2.2117062707308864e-02, + 2.0351096230044517e-02,1.8605121275724643e-02,1.6880083152543166e-02, + 1.5177088307935325e-02,1.3497450601739880e-02,1.1842757857907888e-02, + 1.0214971439701471e-02,8.6165827693987316e-03,7.0508754713732268e-03, + 5.5224032992509968e-03,4.0379725933630305e-03,2.6090727461021627e-03, + 1.2602859304985975e-03] -function randmtzig_fill_randn!(s::DSFMT_state, A) - ccall((:randmtzig_fill_randn,:librandom), - Void, - (Ptr{Void}, Ptr{Float64}, Int), - s.val, A, length(A)) - return A -end +## Tables for exponential variates +const ke = + [0x000e290a13924be4,0 ,0x0009beadebce18c0,0x000c377ac71f9e08, + 0x000d4ddb99075857,0x000de893fb8ca23e,0x000e4a8e87c4328e,0x000e8dff16ae1cba, + 0x000ebf2deab58c5a,0x000ee49a6e8b9639,0x000f0204efd64ee5,0x000f19bdb8ea3c1c, + 0x000f2d458bbe5bd2,0x000f3da104b78236,0x000f4b86d784571f,0x000f577ad8a7784f, + 0x000f61de83da32ac,0x000f6afb7843cce7,0x000f730a57372b44,0x000f7a37651b0e68, + 0x000f80a5bb6eea52,0x000f867189d3cb5c,0x000f8bb1b4f8fbbd,0x000f9079062292b9, + 0x000f94d70ca8d43a,0x000f98d8c7dcaa99,0x000f9c8928abe083,0x000f9ff175b734a6, + 0x000fa319996bc47e,0x000fa6085f8e9d08,0x000fa8c3a62e1991,0x000fab5084e1f660, + 0x000fadb36c84cccb,0x000faff041086846,0x000fb20a6ea22bb9,0x000fb404fb42cb3d, + 0x000fb5e295158174,0x000fb7a59e99727a,0x000fb95038c8789d,0x000fbae44ba684ec, + 0x000fbc638d822e60,0x000fbdcf89209ffb,0x000fbf29a303cfc5,0x000fc0731df1089d, + 0x000fc1ad1ed6c8b1,0x000fc2d8b02b5c8a,0x000fc3f6c4d92131,0x000fc5083ac9ba7d, + 0x000fc60ddd1e9cd7,0x000fc7086622e825,0x000fc7f881009f0c,0x000fc8decb41ac71, + 0x000fc9bbd623d7ec,0x000fca9027c5b26e,0x000fcb5c3c319c49,0x000fcc20864b4449, + 0x000fccdd70a35d41,0x000fcd935e34bf80,0x000fce42ab0db8bd,0x000fceebace7ec02, + 0x000fcf8eb3b0d0e7,0x000fd02c0a049b60,0x000fd0c3f59d199d,0x000fd156b7b5e27e, + 0x000fd1e48d670342,0x000fd26daff73552,0x000fd2f2552684bf,0x000fd372af7233c2, + 0x000fd3eeee528f62,0x000fd4673e73543b,0x000fd4dbc9e72ff8,0x000fd54cb856dc2c, + 0x000fd5ba2f2c4119,0x000fd62451ba02c3,0x000fd68b415fcff5,0x000fd6ef1dabc161, + 0x000fd75004790eb7,0x000fd7ae120c583f,0x000fd809612dbd0a,0x000fd8620b40effa, + 0x000fd8b8285b78fe,0x000fd90bcf594b1d,0x000fd95d15efd426,0x000fd9ac10bfa70c, + 0x000fd9f8d364df06,0x000fda437086566c,0x000fda8bf9e3c9ff,0x000fdad28062fed5, + 0x000fdb17141bff2d,0x000fdb59c4648085,0x000fdb9a9fda83cd,0x000fdbd9b46e3ed4, + 0x000fdc170f6b5d05,0x000fdc52bd81a3fb,0x000fdc8ccacd07ba,0x000fdcc542dd3902, + 0x000fdcfc30bcb794,0x000fdd319ef77143,0x000fdd6597a0f60c,0x000fdd98245a48a3, + 0x000fddc94e575271,0x000fddf91e64014f,0x000fde279ce914cb,0x000fde54d1f0a06b, + 0x000fde80c52a47d0,0x000fdeab7def394e,0x000fded50345eb36,0x000fdefd5be59fa1, + 0x000fdf248e39b26f,0x000fdf4aa064b4b0,0x000fdf6f98435894,0x000fdf937b6f30bb, + 0x000fdfb64f414572,0x000fdfd818d48262,0x000fdff8dd07fed8,0x000fe018a08122c5, + 0x000fe03767adaa5a,0x000fe05536c58a14,0x000fe07211ccb4c5,0x000fe08dfc94c532, + 0x000fe0a8fabe8ca2,0x000fe0c30fbb87a6,0x000fe0dc3ecf3a5a,0x000fe0f48b107522, + 0x000fe10bf76a82ef,0x000fe122869e4200,0x000fe1383b4327e1,0x000fe14d17c83188, + 0x000fe1611e74c023,0x000fe1745169635a,0x000fe186b2a09177,0x000fe19843ef4e08, + 0x000fe1a90705bf64,0x000fe1b8fd6fb37c,0x000fe1c828951444,0x000fe1d689ba4bfd, + 0x000fe1e4220099a5,0x000fe1f0f26655a0,0x000fe1fcfbc726d4,0x000fe2083edc2831, + 0x000fe212bc3bfeb4,0x000fe21c745adfe4,0x000fe225678a8895,0x000fe22d95fa23f4, + 0x000fe234ffb62282,0x000fe23ba4a800d9,0x000fe2418495fddd,0x000fe2469f22bffb, + 0x000fe24af3cce90e,0x000fe24e81ee9859,0x000fe25148bcda1a,0x000fe253474703fe, + 0x000fe2547c75fdc6,0x000fe254e70b7550,0x000fe25485a0fd1a,0x000fe25356a71450, + 0x000fe2515864173b,0x000fe24e88f316f2,0x000fe24ae64296fb,0x000fe2466e132f61, + 0x000fe2411df611bd,0x000fe23af34b6f73,0x000fe233eb40bf42,0x000fe22c02cee01c, + 0x000fe22336b81711,0x000fe2198385e5cd,0x000fe20ee586b707,0x000fe20358cb5dfc, + 0x000fe1f6d92465b1,0x000fe1e9621f2c9e,0x000fe1daef02c8da,0x000fe1cb7accb0a6, + 0x000fe1bb002d22ca,0x000fe1a9798349b9,0x000fe196e0d9140d,0x000fe1832fdebc44, + 0x000fe16e5fe5f932,0x000fe15869dccfcf,0x000fe1414647fe78,0x000fe128ed3cf8b2, + 0x000fe10f565b69cf,0x000fe0f478c633ab,0x000fe0d84b1bdd9e,0x000fe0bac36e6688, + 0x000fe09bd73a6b5c,0x000fe07b7b5d920b,0x000fe059a40c26d2,0x000fe03644c5d7f9, + 0x000fe011504979b3,0x000fdfeab887b95d,0x000fdfc26e94a448,0x000fdf986297e306, + 0x000fdf6c83bb8663,0x000fdf3ec0193eee,0x000fdf0f04a5d30a,0x000fdedd3d1aa204, + 0x000fdea953dcfc13,0x000fde7331e3100e,0x000fde3abe9626f3,0x000fddffdfb1dbd5, + 0x000fddc2791ff351,0x000fdd826cd068c7,0x000fdd3f9a8d3857,0x000fdcf9dfc95b0d, + 0x000fdcb1176a55fe,0x000fdc65198ba50c,0x000fdc15bb3b2daa,0x000fdbc2ce2dc4ae, + 0x000fdb6c206aaaca,0x000fdb117becb4a2,0x000fdab2a6379bf1,0x000fda4f5fdfb4e9, + 0x000fd9e76401f3a4,0x000fd97a67a9ce20,0x000fd9081922142a,0x000fd8901f2d4b02, + 0x000fd812182170e1,0x000fd78d98e23cd4,0x000fd7022bb3f083,0x000fd66f4edf96ba, + 0x000fd5d473200306,0x000fd530f9ccff94,0x000fd48432b7b351,0x000fd3cd59a8469f, + 0x000fd30b9368f90a,0x000fd23dea45f500,0x000fd16349e2e04b,0x000fd07a7a3ef98b, + 0x000fcf8219b5df06,0x000fce7895bcfcdf,0x000fcd5c220ad5e3,0x000fcc2aadbc17dd, + 0x000fcae1d5e81fbd,0x000fc97ed4e778fa,0x000fc7fe6d4d720f,0x000fc65ccf39c2fc, + 0x000fc4957623cb04,0x000fc2a2fc826dc8,0x000fc07ee19b01ce,0x000fbe213c1cf493, + 0x000fbb8051ac1566,0x000fb890078d120e,0x000fb5411a5b9a96,0x000fb18000547134, + 0x000fad334827f1e2,0x000fa839276708b9,0x000fa263b32e37ee,0x000f9b72d1c52cd1, + 0x000f930a1a281a05,0x000f889f023d820a,0x000f7b577d2be5f4,0x000f69c650c40a8f, + 0x000f51530f0916d9,0x000f2cb0e3c5933e,0x000eeefb15d605d9,0x000e6da6ecf27460] -function randmtzig_gv_randn() - ccall((:randmtzig_gv_randn,:librandom), - Float64, - ()) -end +const we = + [1.9311480126418366e-15,1.4178028487910829e-17,2.3278824993382448e-17, + 3.0487830247064320e-17,3.6665697714474878e-17,4.2179302189289733e-17, + 4.7222561556862764e-17,5.1911915446217879e-17,5.6323471083955047e-17, + 6.0510082606427647e-17,6.4510165096727506e-17,6.8352646803700541e-17, + 7.2059939574689050e-17,7.5649815537392981e-17,7.9136643961951065e-17, + 8.2532235563518929e-17,8.5846436168850513e-17,8.9087554865647428e-17, + 9.2262679629663719e-17,9.5377914505292719e-17,9.8438560874559257e-17, + 1.0144925809006294e-16,1.0441409405585343e-16,1.0733669323436384e-16, + 1.1022028745670189e-16,1.1306777346479334e-16,1.1588176009705533e-16, + 1.1866460730417886e-16,1.2141845865694359e-16,1.2414526862326387e-16, + 1.2684682560606153e-16,1.2952477151912284e-16,1.3218061851538810e-16, + 1.3481576335745444e-16,1.3743149982367625e-16,1.4002902946807859e-16, + 1.4260947099321287e-16,1.4517386844829297e-16,1.4772319842763584e-16, + 1.5025837641447456e-16,1.5278026239101652e-16,1.5528966581595696e-16, + 1.5778735005459581e-16,1.6027403633350909e-16,1.6275040728083524e-16, + 1.6521711010420076e-16,1.6767475945078279e-16,1.7012393998770646e-16, + 1.7256520873568226e-16,1.7499909718432365e-16,1.7742611321380505e-16, + 1.7984674284430714e-16,1.8226145183195818e-16,1.8467068712763576e-16, + 1.8707487821298258e-16,1.8947443832625899e-16,1.9186976558915995e-16, + 1.9426124404443042e-16,1.9664924461299023e-16,1.9903412597830144e-16, + 2.0141623540485899e-16,2.0379590949693882e-16,2.0617347490308439e-16, + 2.0854924897123771e-16,2.1092354035891528e-16,2.1329664960238294e-16, + 2.1566886964838970e-16,2.1804048635167009e-16,2.2041177894111562e-16, + 2.2278302045723950e-16,2.2515447816331350e-16,2.2752641393233694e-16, + 2.2989908461180186e-16,2.3227274236804366e-16,2.3464763501180916e-16, + 2.3702400630653389e-16,2.3940209626069303e-16,2.4178214140547710e-16, + 2.4416437505894123e-16,2.4654902757768304e-16,2.4893632659702250e-16, + 2.5132649726057970e-16,2.5371976244007951e-16,2.5611634294614988e-16, + 2.5851645773082391e-16,2.6092032408240577e-16,2.6332815781331452e-16, + 2.6574017344147618e-16,2.6815658436579989e-16,2.7057760303623509e-16, + 2.7300344111887955e-16,2.7543430965657619e-16,2.7787041922541278e-16, + 2.8031198008751431e-16,2.8275920234049704e-16,2.8521229606393309e-16, + 2.8767147146315804e-16,2.9013693901073754e-16,2.9260890958589514e-16, + 2.9508759461219033e-16,2.9757320619372521e-16,3.0006595725014739e-16, + 3.0256606165070789e-16,3.0507373434762511e-16,3.0758919150899939e-16, + 3.1011265065151543e-16,3.1264433077316750e-16,3.1518445248623523e-16, + 3.1773323815073683e-16,3.2029091200858335e-16,3.2285770031865573e-16, + 3.2543383149302610e-16,3.2801953623454359e-16,3.3061504767600738e-16, + 3.3322060152114841e-16,3.3583643618764577e-16,3.3846279295240445e-16, + 3.4109991609932597e-16,3.4374805306980633e-16,3.4640745461620167e-16, + 3.4907837495850680e-16,3.5176107194449828e-16,3.5445580721360130e-16, + 3.5716284636474652e-16,3.5988245912849274e-16,3.6261491954370031e-16, + 3.6536050613905045e-16,3.6811950211971757e-16,3.7089219555951389e-16, + 3.7367887959883854e-16,3.7647985264877841e-16,3.7929541860172334e-16, + 3.8212588704887531e-16,3.8497157350504876e-16,3.8783279964117988e-16, + 3.9070989352498183e-16,3.9360318987020748e-16,3.9651303029500381e-16, + 3.9943976358986842e-16,4.0238374599574693e-16,4.0534534149283966e-16, + 4.0832492210071775e-16,4.1132286819038357e-16,4.1433956880894741e-16, + 4.1737542201763194e-16,4.2043083524385856e-16,4.2350622564821518e-16, + 4.2660202050715582e-16,4.2971865761233266e-16,4.3285658568752094e-16, + 4.3601626482415681e-16,4.3919816693657415e-16,4.4240277623809919e-16, + 4.4563058973923611e-16,4.4888211776926172e-16,4.5215788452263475e-16, + 4.5545842863172421e-16,4.5878430376746227e-16,4.6213607926964266e-16, + 4.6551434080870692e-16,4.6891969108099157e-16,4.7235275053955480e-16, + 4.7581415816285534e-16,4.7930457226372470e-16,4.8282467134125866e-16, + 4.8637515497845119e-16,4.8995674478861404e-16,4.9357018541385775e-16, + 4.9721624557917034e-16,5.0089571920591141e-16,5.0460942658884340e-16, + 5.0835821564116245e-16,5.1214296321235415e-16,5.1596457648410618e-16, + 5.1982399444994938e-16,5.2372218948478484e-16,5.2766016901098856e-16, + 5.3163897726836902e-16,5.3565969719590503e-16,5.3972345243389779e-16, + 5.4383140945596370e-16,5.4798477984116296e-16,5.5218482269752343e-16, + 5.5643284724928722e-16,5.6073021560139669e-16,5.6507834569605064e-16, + 5.6947871447763482e-16,5.7393286128396354e-16,5.7844239148359912e-16, + 5.8300898038105864e-16,5.8763437741400573e-16,5.9232041066909314e-16, + 5.9706899174600906e-16,6.0188212100252363e-16,6.0676189321700068e-16, + 6.1171050370897217e-16,6.1673025496306200e-16,6.2182356380685327e-16, + 6.2699296919933262e-16,6.3224114069342115e-16,6.3757088764394262e-16, + 6.4298516924135947e-16,6.4848710546189033e-16,6.5407998903644809e-16, + 6.5976729855445663e-16,6.6555271283433428e-16,6.7144012671064882e-16, + 6.7743366840910103e-16,6.8353771870512740e-16,6.8975693209068478e-16, + 6.9609626020748846e-16,7.0256097784459588e-16,7.0915671184495837e-16, + 7.1588947332085531e-16,7.2276569364381212e-16,7.2979226475290851e-16, + 7.3697658441912426e-16,7.4432660721604146e-16,7.5185090208325131e-16, + 7.5955871753377488e-16,7.6746005575784274e-16,7.7556575712157906e-16, + 7.8388759686228577e-16,7.9243839615735500e-16,8.0123215021130834e-16, + 8.1028417659131464e-16,8.1961128778061250e-16,8.2923199285818092e-16, + 8.3916673441467979e-16,8.4943816836487701e-16,8.6007149633349414e-16, + 8.7109486293879040e-16,8.8253983380721398e-16,8.9444197485198646e-16, + 9.0684155971316690e-16,9.1978444098118649e-16,9.3332313294229516e-16, + 9.4751817065249841e-16,9.6243983456584759e-16,9.7817036547844198e-16, + 9.9480684723838795e-16,1.0124650144288319e-15,1.0312843657756166e-15, + 1.0514351604044550e-15,1.0731281954224043e-15,1.0966288068517408e-15, + 1.1222774909350319e-15,1.1505212963006663e-15,1.1819635283304206e-15, + 1.2174462832361815e-15,1.2581958069755114e-15,1.3060984107128082e-15, + 1.3642786158057857e-15,1.4384889932178723e-15,1.5412190700064194e-15, + 1.7091034077168055e-15] +const fe = + [1.0000000000000000e+00,9.3814368086217470e-01,9.0046992992574648e-01, + 8.7170433238120359e-01,8.4778550062398961e-01,8.2699329664305032e-01, + 8.0842165152300838e-01,7.9152763697249562e-01,7.7595685204011555e-01, + 7.6146338884989628e-01,7.4786862198519510e-01,7.3503809243142348e-01, + 7.2286765959357202e-01,7.1127476080507601e-01,7.0019265508278816e-01, + 6.8956649611707799e-01,6.7935057226476536e-01,6.6950631673192473e-01, + 6.6000084107899970e-01,6.5080583341457110e-01,6.4189671642726609e-01, + 6.3325199421436607e-01,6.2485273870366598e-01,6.1668218091520766e-01, + 6.0872538207962201e-01,6.0096896636523223e-01,5.9340090169173343e-01, + 5.8601031847726803e-01,5.7878735860284503e-01,5.7172304866482582e-01, + 5.6480919291240017e-01,5.5803828226258745e-01,5.5140341654064129e-01, + 5.4489823767243961e-01,5.3851687200286191e-01,5.3225388026304332e-01, + 5.2610421398361973e-01,5.2006317736823360e-01,5.1412639381474856e-01, + 5.0828977641064288e-01,5.0254950184134772e-01,4.9690198724154955e-01, + 4.9134386959403253e-01,4.8587198734188491e-01,4.8048336393045421e-01, + 4.7517519303737737e-01,4.6994482528395998e-01,4.6478975625042618e-01, + 4.5970761564213769e-01,4.5469615747461550e-01,4.4975325116275500e-01, + 4.4487687341454851e-01,4.4006510084235390e-01,4.3531610321563657e-01, + 4.3062813728845883e-01,4.2599954114303434e-01,4.2142872899761658e-01, + 4.1691418643300288e-01,4.1245446599716118e-01,4.0804818315203240e-01, + 4.0369401253053028e-01,3.9939068447523107e-01,3.9513698183329016e-01, + 3.9093173698479711e-01,3.8677382908413765e-01,3.8266218149600983e-01, + 3.7859575940958079e-01,3.7457356761590216e-01,3.7059464843514600e-01, + 3.6665807978151416e-01,3.6276297335481777e-01,3.5890847294874978e-01, + 3.5509375286678746e-01,3.5131801643748334e-01,3.4758049462163698e-01, + 3.4388044470450241e-01,3.4021714906678002e-01,3.3658991402867761e-01, + 3.3299806876180899e-01,3.2944096426413633e-01,3.2591797239355619e-01, + 3.2242848495608917e-01,3.1897191284495724e-01,3.1554768522712895e-01, + 3.1215524877417955e-01,3.0879406693456019e-01,3.0546361924459026e-01, + 3.0216340067569353e-01,2.9889292101558179e-01,2.9565170428126120e-01, + 2.9243928816189257e-01,2.8925522348967775e-01,2.8609907373707683e-01, + 2.8297041453878075e-01,2.7986883323697292e-01,2.7679392844851736e-01, + 2.7374530965280297e-01,2.7072259679906002e-01,2.6772541993204479e-01, + 2.6475341883506220e-01,2.6180624268936298e-01,2.5888354974901623e-01, + 2.5598500703041538e-01,2.5311029001562946e-01,2.5025908236886230e-01, + 2.4743107566532763e-01,2.4462596913189211e-01,2.4184346939887721e-01, + 2.3908329026244918e-01,2.3634515245705964e-01,2.3362878343743335e-01, + 2.3093391716962741e-01,2.2826029393071670e-01,2.2560766011668407e-01, + 2.2297576805812019e-01,2.2036437584335949e-01,2.1777324714870053e-01, + 2.1520215107537868e-01,2.1265086199297828e-01,2.1011915938898826e-01, + 2.0760682772422204e-01,2.0511365629383771e-01,2.0263943909370902e-01, + 2.0018397469191127e-01,1.9774706610509887e-01,1.9532852067956322e-01, + 1.9292814997677135e-01,1.9054576966319539e-01,1.8818119940425432e-01, + 1.8583426276219711e-01,1.8350478709776746e-01,1.8119260347549629e-01, + 1.7889754657247831e-01,1.7661945459049488e-01,1.7435816917135349e-01, + 1.7211353531532006e-01,1.6988540130252766e-01,1.6767361861725019e-01, + 1.6547804187493600e-01,1.6329852875190182e-01,1.6113493991759203e-01, + 1.5898713896931421e-01,1.5685499236936523e-01,1.5473836938446808e-01, + 1.5263714202744286e-01,1.5055118500103989e-01,1.4848037564386679e-01, + 1.4642459387834494e-01,1.4438372216063478e-01,1.4235764543247220e-01, + 1.4034625107486245e-01,1.3834942886358020e-01,1.3636707092642886e-01, + 1.3439907170221363e-01,1.3244532790138752e-01,1.3050573846833077e-01, + 1.2858020454522817e-01,1.2666862943751067e-01,1.2477091858083096e-01, + 1.2288697950954514e-01,1.2101672182667483e-01,1.1916005717532768e-01, + 1.1731689921155557e-01,1.1548716357863353e-01,1.1367076788274431e-01, + 1.1186763167005630e-01,1.1007767640518538e-01,1.0830082545103380e-01, + 1.0653700405000166e-01,1.0478613930657017e-01,1.0304816017125772e-01, + 1.0132299742595363e-01,9.9610583670637132e-02,9.7910853311492199e-02, + 9.6223742550432798e-02,9.4549189376055859e-02,9.2887133556043541e-02, + 9.1237516631040155e-02,8.9600281910032858e-02,8.7975374467270218e-02, + 8.6362741140756913e-02,8.4762330532368119e-02,8.3174093009632383e-02, + 8.1597980709237419e-02,8.0033947542319905e-02,7.8481949201606421e-02, + 7.6941943170480503e-02,7.5413888734058410e-02,7.3897746992364746e-02, + 7.2393480875708738e-02,7.0901055162371829e-02,6.9420436498728755e-02, + 6.7951593421936601e-02,6.6494496385339774e-02,6.5049117786753749e-02, + 6.3615431999807334e-02,6.2193415408540995e-02,6.0783046445479633e-02, + 5.9384305633420266e-02,5.7997175631200659e-02,5.6621641283742877e-02, + 5.5257689676697037e-02,5.3905310196046087e-02,5.2564494593071692e-02, + 5.1235237055126281e-02,4.9917534282706372e-02,4.8611385573379497e-02, + 4.7316792913181548e-02,4.6033761076175170e-02,4.4762297732943282e-02, + 4.3502413568888183e-02,4.2254122413316234e-02,4.1017441380414819e-02, + 3.9792391023374125e-02,3.8578995503074857e-02,3.7377282772959361e-02, + 3.6187284781931423e-02,3.5009037697397410e-02,3.3842582150874330e-02, + 3.2687963508959535e-02,3.1545232172893609e-02,3.0414443910466604e-02, + 2.9295660224637393e-02,2.8188948763978636e-02,2.7094383780955800e-02, + 2.6012046645134217e-02,2.4942026419731783e-02,2.3884420511558171e-02, + 2.2839335406385240e-02,2.1806887504283581e-02,2.0787204072578117e-02, + 1.9780424338009743e-02,1.8786700744696030e-02,1.7806200410911362e-02, + 1.6839106826039948e-02,1.5885621839973163e-02,1.4945968011691148e-02, + 1.4020391403181938e-02,1.3109164931254991e-02,1.2212592426255381e-02, + 1.1331013597834597e-02,1.0464810181029979e-02,9.6144136425022099e-03, + 8.7803149858089753e-03,7.9630774380170400e-03,7.1633531836349839e-03, + 6.3819059373191791e-03,5.6196422072054830e-03,4.8776559835423923e-03, + 4.1572951208337953e-03,3.4602647778369040e-03,2.7887987935740761e-03, + 2.1459677437189063e-03,1.5362997803015724e-03,9.6726928232717454e-04, + 4.5413435384149677e-04] -function randmtzig_gv_fill_randn!(A) - ccall((:randmtzig_gv_fill_randn,:librandom), - Void, - (Ptr{Float64}, Int), - A, length(A)) - return A -end +ziggurat_nor_r = 3.6541528853610087963519472518 +ziggurat_nor_inv_r = inv(ziggurat_nor_r) +ziggurat_exp_r = 7.6971174701310497140446280481 -function randmtzig_exprnd() - ccall((:randmtzig_exprnd,:librandom), - Float64, - ()) -end +randi() = reinterpret(Uint64,dsfmt_gv_genrand_close1_open2()) & 0x000fffffffffffff +randi(state::DSFMT_state) = reinterpret(Uint64,dsfmt_genrand_close1_open2(state)) & 0x000fffffffffffff +for (lhs, rhs) in (([], []), + ([:(state::DSFMT_state)], [:state])) + @eval begin + function randmtzig_randn($(lhs...)) + @inbounds begin + while true + r = randi($(rhs...)) + rabs = int64(r>>1) # One bit for the sign + idx = rabs & 0xFF + x = (r&1 != 0x000000000 ? -rabs : rabs)*wi[idx+1] + if rabs < ki[idx+1] + return x # 99.3% of the time we return here 1st try + elseif idx == 0 + while true + xx = -$(ziggurat_nor_inv_r)*log(rand($(rhs...))) + yy = -log(rand($(rhs...))) + if yy+yy > xx*xx + return (rabs & 0x100) != 0x000000000 ? -$(ziggurat_nor_r)-xx : $(ziggurat_nor_r)+xx + end + end + elseif (fi[idx] - fi[idx+1])*rand($(rhs...)) + fi[idx+1] < exp(-0.5*x*x) + return x # return from the triangular area + end + end + end + end -function randmtzig_fill_exprnd!(A) - ccall((:randmtzig_fill_exprnd,:librandom), - Void, - (Ptr{Float64}, Int), - A, length(A)) - return A + function randmtzig_exprnd($(lhs...)) + @inbounds begin + while true + ri = randi($(rhs...)) + idx = ri & 0xFF + x = ri*we[idx+1] + if ri < ke[idx+1] + return x # 98.9% of the time we return here 1st try + elseif idx == 0 + return $(ziggurat_exp_r) - log(rand($(rhs...))) + elseif (fe[idx] - fe[idx+1])*rand($(rhs...)) + fe[idx+1] < exp(-x) + return x # return from the triangular area + end + end + end + end + end end - ## Windows entropy @windows_only begin diff --git a/base/random.jl b/base/random.jl index edada50dfc763..fbdd647b40ad8 100644 --- a/base/random.jl +++ b/base/random.jl @@ -213,10 +213,10 @@ randbool!(B::BitArray) = rand!(B) # The Ziggurat Method for generating random variables - Marsaglia and Tsang # Paper and reference code: http://www.jstatsoft.org/v05/i08/ -randn() = randmtzig_gv_randn() +randn() = randmtzig_randn() randn(rng::MersenneTwister) = randmtzig_randn(rng.state) -randn!(A::Array{Float64}) = randmtzig_gv_fill_randn!(A) -randn!(rng::MersenneTwister, A::Array{Float64}) = randmtzig_fill_randn!(rng.state, A) +randn!(A::Array{Float64}) = (for i = 1:length(A);A[i] = randmtzig_randn();end;A) +randn!(rng::MersenneTwister, A::Array{Float64}) = (for i = 1:length(A);A[i] = randmtzig_randn(rng.state);end;A) randn(dims::Dims) = randn!(Array(Float64, dims)) randn(dims::Int...) = randn!(Array(Float64, dims...)) diff --git a/base/sysimg.jl b/base/sysimg.jl index 98a4570a704c4..0cc475085094e 100644 --- a/base/sysimg.jl +++ b/base/sysimg.jl @@ -149,12 +149,6 @@ end # reduction along dims include("reducedim.jl") # macros in this file relies on string.jl -# random number generation and statistics -include("statistics.jl") -include("librandom.jl") -include("random.jl") -importall .Random - # basic data structures include("ordering.jl") importall .Order @@ -165,6 +159,27 @@ include("sort.jl") importall .Sort include("combinatorics.jl") +# rounding utilities +include("rounding.jl") +importall .Rounding + +# BigInts and BigFloats +include("gmp.jl") +importall .GMP +include("mpfr.jl") +importall .MPFR +big(n::Integer) = convert(BigInt,n) +big(x::FloatingPoint) = convert(BigFloat,x) +big(q::Rational) = big(num(q))//big(den(q)) +big(z::Complex) = complex(big(real(z)),big(imag(z))) +@vectorize_1arg Number big + +# random number generation and statistics +include("statistics.jl") +include("librandom.jl") +include("random.jl") +importall .Random + # distributed arrays and memory-mapped arrays include("darray.jl") include("mmap.jl") @@ -209,21 +224,6 @@ include("fftw.jl") include("dsp.jl") importall .DSP -# rounding utilities -include("rounding.jl") -importall .Rounding - -# BigInts and BigFloats -include("gmp.jl") -importall .GMP -include("mpfr.jl") -importall .MPFR -big(n::Integer) = convert(BigInt,n) -big(x::FloatingPoint) = convert(BigFloat,x) -big(q::Rational) = big(num(q))//big(den(q)) -big(z::Complex) = complex(big(real(z)),big(imag(z))) -@vectorize_1arg Number big - # (s)printf macros include("printf.jl") importall .Printf diff --git a/deps/Makefile b/deps/Makefile index 487c160ebec13..7518a41990653 100644 --- a/deps/Makefile +++ b/deps/Makefile @@ -614,7 +614,7 @@ install-openspecfun: $(OPENSPECFUN_OBJ_TARGET) ## LIBRANDOM ## LIBRANDOM_OBJ_TARGET = $(build_shlibdir)/librandom.$(SHLIB_EXT) -LIBRANDOM_OBJ_SOURCE = random/librandom.$(SHLIB_EXT) +LIBRANDOM_OBJ_SOURCE = librandom.$(SHLIB_EXT) LIBRANDOM_CFLAGS = $(CFLAGS) -DNDEBUG -DDSFMT_MEXP=19937 $(fPIC) -DDSFMT_DO_NOT_USE_OLD_NAMES ifneq ($(USEMSVC), 1) @@ -627,29 +627,27 @@ ifeq ($(ARCH), x86_64) LIBRANDOM_CFLAGS += -msse2 -DHAVE_SSE2 endif -random/dsfmt-$(DSFMT_VER).tar.gz: +dsfmt-$(DSFMT_VER).tar.gz: $(JLDOWNLOAD) $@ http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/dSFMT-src-$(DSFMT_VER).tar.gz touch -c $@ -random/dsfmt-$(DSFMT_VER)/config.status: random/dsfmt-$(DSFMT_VER).tar.gz - cd random && \ +dsfmt-$(DSFMT_VER)/config.status: dsfmt-$(DSFMT_VER).tar.gz mkdir -p dsfmt-$(DSFMT_VER) && \ $(TAR) -C dsfmt-$(DSFMT_VER) --strip-components 1 -xf dsfmt-$(DSFMT_VER).tar.gz && \ cd dsfmt-$(DSFMT_VER) && patch < ../dSFMT.h.patch && patch < ../dSFMT.c.patch echo 1 > $@ -$(LIBRANDOM_OBJ_SOURCE): random/jl_random.c random/randmtzig.c random/dsfmt-$(DSFMT_VER)/config.status - cd random && \ - $(CC) $(CPPFLAGS) $(LIBRANDOM_CFLAGS) $(LDFLAGS) jl_random.c -o librandom.$(SHLIB_EXT) && \ +$(LIBRANDOM_OBJ_SOURCE): dsfmt-$(DSFMT_VER)/config.status dsfmt-$(DSFMT_VER)/dSFMT.c + $(CC) $(CPPFLAGS) $(LIBRANDOM_CFLAGS) $(LDFLAGS) dsfmt-$(DSFMT_VER)/dSFMT.c -o librandom.$(SHLIB_EXT) && \ $(INSTALL_NAME_CMD)librandom.$(SHLIB_EXT) librandom.$(SHLIB_EXT) $(LIBRANDOM_OBJ_TARGET): $(LIBRANDOM_OBJ_SOURCE) cp $< $@ clean-random: - -rm -f random/librandom.$(SHLIB_EXT) + -rm -f librandom.$(SHLIB_EXT) distclean-random: clean-random - -cd random && rm -rf *.tar.gz dsfmt-$(DSFMT_VER) + -rf *.tar.gz dsfmt-$(DSFMT_VER) -get-random: random/dsfmt-$(DSFMT_VER).tar.gz -configure-random: random/dsfmt-$(DSFMT_VER)/config.status +get-random: dsfmt-$(DSFMT_VER).tar.gz +configure-random: dsfmt-$(DSFMT_VER)/config.status compile-random: $(LIBRANDOM_OBJ_SOURCE) check-random: compile-random install-random: $(LIBRANDOM_OBJ_TARGET) diff --git a/deps/random/dSFMT.c.patch b/deps/dSFMT.c.patch similarity index 100% rename from deps/random/dSFMT.c.patch rename to deps/dSFMT.c.patch diff --git a/deps/random/dSFMT.h.patch b/deps/dSFMT.h.patch similarity index 100% rename from deps/random/dSFMT.h.patch rename to deps/dSFMT.h.patch diff --git a/deps/random/.gitignore b/deps/random/.gitignore deleted file mode 100644 index 9baf7f91a9395..0000000000000 --- a/deps/random/.gitignore +++ /dev/null @@ -1,5 +0,0 @@ -/dsfmt-* -/mt19937-64* -/libMT* -/librandom.so -/librandom.dylib diff --git a/deps/random/jl_random.c b/deps/random/jl_random.c deleted file mode 100644 index 3159d8d669267..0000000000000 --- a/deps/random/jl_random.c +++ /dev/null @@ -1,6 +0,0 @@ -/* - random numbers -*/ - -#include "dsfmt-2.2/dSFMT.c" -#include "randmtzig.c" diff --git a/deps/random/randmtzig.c b/deps/random/randmtzig.c deleted file mode 100644 index 218d534352c83..0000000000000 --- a/deps/random/randmtzig.c +++ /dev/null @@ -1,853 +0,0 @@ -/* - Parts Copyright (C) 2012 Viral B. Shah - All rights reserved. - - Modifications made for julia to support dsfmt and only __LP64__ systems. - 52-bits of randomness are used from the mantissa of random double precision - numbers generated by dsfmt. - */ - -/* - A C-program for MT19937, with initialization improved 2002/2/10. - Coded by Takuji Nishimura and Makoto Matsumoto. - This is a faster version by taking Shawn Cokus's optimization, - Matthe Bellew's simplification, Isaku Wada's real version. - David Bateman added normal and exponential distributions following - Marsaglia and Tang's Ziggurat algorithm. - - Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, - Copyright (C) 2004, David Bateman - All rights reserved. - - Redistribution and use in source and binary forms, with or without - modification, are permitted provided that the following conditions - are met: - - 1. Redistributions of source code must retain the above copyright - notice, this list of conditions and the following disclaimer. - - 2. Redistributions in binary form must reproduce the above copyright - notice, this list of conditions and the following disclaimer in the - documentation and/or other materials provided with the distribution. - - 3. The names of its contributors may not be used to endorse or promote - products derived from this software without specific prior written - permission. - - THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS - "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT - LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR - A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER - OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, - EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, - PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR - PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF - LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING - NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS - SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - - - Any feedback is very welcome. - http://www.math.keio.ac.jp/matumoto/emt.html - email: matumoto@math.keio.ac.jp - -*/ - -#include -#include -#include -#include -#ifndef _MSC_VER -#include -#endif - -#ifdef STANDALONE -#include -#include - -#define DSFMT_DO_NOT_USE_OLD_NAMES -#include "dsfmt-2.2/dSFMT.c" -#endif - -typedef ptrdiff_t randmtzig_idx_type; -typedef signed char randmtzig_int8_t; -typedef unsigned char randmtzig_uint8_t; -typedef short randmtzig_int16_t; -typedef unsigned short randmtzig_uint16_t; -typedef int randmtzig_int32_t; -typedef unsigned int randmtzig_uint32_t; -typedef long long randmtzig_int64_t; -typedef unsigned long long randmtzig_uint64_t; - -/* Declarations */ - -extern double randmtzig_randn (dsfmt_t *dsfmt); -extern void randmtzig_fill_randn (dsfmt_t *dsfmt, double *p, randmtzig_idx_type n); -extern double randmtzig_gv_randn (void); -extern void randmtzig_gv_fill_randn (double *p, randmtzig_idx_type n); -extern double randmtzig_exprnd (void); -extern void randmtzig_fill_exprnd (double *p, randmtzig_idx_type n); - -/* ===== Uniform generators ===== */ - -inline static randmtzig_uint64_t gv_randi (void) -{ - double r = dsfmt_gv_genrand_close1_open2(); - return *((uint64_t *) &r) & 0x000fffffffffffff; -} - -/* generates a random number on (0,1) with 53-bit resolution */ -inline static double gv_randu (void) -{ - return dsfmt_gv_genrand_open_open(); -} - -inline static randmtzig_uint64_t randi (dsfmt_t *dsfmt) -{ - double r = dsfmt_genrand_close1_open2(dsfmt); - return *((uint64_t *) &r) & 0x000fffffffffffff; -} - -/* generates a random number on (0,1) with 53-bit resolution */ -inline static double randu (dsfmt_t *dsfmt) -{ - return dsfmt_genrand_open_open(dsfmt); -} - -/* ===== Ziggurat normal and exponential generators ===== */ -# define ZIGINT randmtzig_uint64_t -# define EMANTISSA 4503599627370496 /* 52 bit mantissa */ -# define ERANDI gv_randi() /* 52 bits for mantissa */ -# define NMANTISSA 2251799813685248 -# define NRANDI gv_randi() /* 51 bits for mantissa + 1 bit sign */ -# define RANDU gv_randu() - -#define ZIGGURAT_TABLE_SIZE 256 - -#define ZIGGURAT_NOR_R 3.6541528853610088 -#define ZIGGURAT_NOR_INV_R 0.27366123732975828 -#define NOR_SECTION_AREA 0.00492867323399 - -#define ZIGGURAT_EXP_R 7.69711747013104972 -#define ZIGGURAT_EXP_INV_R 0.129918765548341586 -#define EXP_SECTION_AREA 0.0039496598225815571993 - - -/* -This code is based on the paper Marsaglia and Tsang, "The ziggurat method -for generating random variables", Journ. Statistical Software. Code was -presented in this paper for a Ziggurat of 127 levels and using a 32 bit -integer random number generator. This version of the code, uses the -Mersenne Twister as the integer generator and uses 256 levels in the -Ziggurat. This has several advantages. - - 1) As Marsaglia and Tsang themselves states, the more levels the few - times the expensive tail algorithm must be called - 2) The cycle time of the generator is determined by the integer - generator, thus the use of a Mersenne Twister for the core random - generator makes this cycle extremely long. - 3) The license on the original code was unclear, thus rewriting the code - from the article means we are free of copyright issues. - 4) Compile flag for full 53-bit random mantissa. - -It should be stated that the authors made my life easier, by the fact that -the algorithm developed in the text of the article is for a 256 level -ziggurat, even if the code itself isn't... - -One modification to the algorithm developed in the article, is that it is -assumed that 0 <= x < Inf, and "unsigned long"s are used, thus resulting in -terms like 2^32 in the code. As the normal distribution is defined between --Inf < x < Inf, we effectively only have 31 bit integers plus a sign. Thus -in Marsaglia and Tsang, terms like 2^32 become 2^31. We use NMANTISSA for -this term. The exponential distribution is one sided so we use the -full 32 bits. We use EMANTISSA for this term. - -It appears that I'm slightly slower than the code in the article, this -is partially due to a better generator of random integers than they -use. But might also be that the case of rapid return was optimized by -inlining the relevant code with a #define. As the basic Mersenne -Twister is only 25% faster than this code I suspect that the main -reason is just the use of the Mersenne Twister and not the inlining, -so I'm not going to try and optimize further. -*/ - - -// void randmtzig_create_ziggurat_tables (void) -// { -// int i; -// double x, x1; - -// /* Ziggurat tables for the normal distribution */ -// x1 = ZIGGURAT_NOR_R; -// wi[255] = x1 / NMANTISSA; -// fi[255] = exp (-0.5 * x1 * x1); - -// /* Index zero is special for tail strip, where Marsaglia and Tsang -// * defines this as -// * k_0 = 2^31 * r * f(r) / v, w_0 = 0.5^31 * v / f(r), f_0 = 1, -// * where v is the area of each strip of the ziggurat. -// */ -// ki[0] = (ZIGINT) (x1 * fi[255] / NOR_SECTION_AREA * NMANTISSA); -// wi[0] = NOR_SECTION_AREA / fi[255] / NMANTISSA; -// fi[0] = 1.; - -// for (i = 254; i > 0; i--) -// { -// /* New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus -// * need inverse operator of y = exp(-0.5*x*x) -> x = sqrt(-2*ln(y)) -// */ -// x = sqrt(-2. * log(NOR_SECTION_AREA / x1 + fi[i+1])); -// ki[i+1] = (ZIGINT)(x / x1 * NMANTISSA); -// wi[i] = x / NMANTISSA; -// fi[i] = exp (-0.5 * x * x); -// x1 = x; -// } - -// ki[1] = 0; - -// /* Zigurrat tables for the exponential distribution */ -// x1 = ZIGGURAT_EXP_R; -// we[255] = x1 / EMANTISSA; -// fe[255] = exp (-x1); - -// /* Index zero is special for tail strip, where Marsaglia and Tsang -// * defines this as -// * k_0 = 2^32 * r * f(r) / v, w_0 = 0.5^32 * v / f(r), f_0 = 1, -// * where v is the area of each strip of the ziggurat. -// */ -// ke[0] = (ZIGINT) (x1 * fe[255] / EXP_SECTION_AREA * EMANTISSA); -// we[0] = EXP_SECTION_AREA / fe[255] / EMANTISSA; -// fe[0] = 1.; - -// for (i = 254; i > 0; i--) -// { -// /* New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus -// * need inverse operator of y = exp(-x) -> x = -ln(y) -// */ -// x = - log(EXP_SECTION_AREA / x1 + fe[i+1]); -// ke[i+1] = (ZIGINT)(x / x1 * EMANTISSA); -// we[i] = x / EMANTISSA; -// fe[i] = exp (-x); -// x1 = x; -// } -// ke[1] = 0; -// } - -// Tables for randn -static ZIGINT ki[ZIGGURAT_TABLE_SIZE] = - {2104047571230236, 0,1693657211688499,1919380038164751, - 2015384402142420,2068365869415708,2101878624030987,2124958784087614, - 2141808670783638,2154644611559370,2164744887580145,2172897953690771, - 2179616279367521,2185247251864556,2190034623104318,2194154434518163, - 2197736978772008,2200880740889623,2203661538008543,2206138681107245, - 2208359231804928,2210361007256700,2212174742387166,2213825672703393, - 2215334711001466,2216719334486539,2217994262138197,2219171977964129, - 2220263139537873,2221276900116549,2222221164932202,2223102796828387, - 2223927782546019,2224701368169460,2225428170203747,2226112267247709, - 2226757276104752,2227366415327922,2227942558554233,2228488279492093, - 2229005890046815,2229497472774805,2229964908626691,2230409900758245, - 2230833995044249,2231238597815812,2231624991249884,2231994346765634, - 2232347736722468,2232686144665663,2233010474325699,2233321557544631, - 2233620161275830,2233906993781039,2234182710130112,2234447917093281, - 2234703177502812,2234949014149981,2235185913274123,2235414327692697, - 2235634679614740,2235847363174420,2236052746716668,2236251174862705, - 2236442970379808,2236628435876608,2236807855342616,2236981495548416, - 2237149607321006,2237312426707072,2237470176035519,2237623064889274, - 2237771290995262,2237915041040474,2238054491421185,2238189808931596, - 2238321151397547,2238448668260322,2238572501115061,2238692784207837, - 2238809644895031,2238923204068302,2239033576548092,2239140871448347, - 2239245192514865,2239346638439450,2239445303151863,2239541276091355, - 2239634642459413,2239725483455210,2239813876495104,2239899895417414, - 2239983610673598,2240065089506859,2240144396119109,2240221591827156, - 2240296735208897,2240369882240222,2240441086423317,2240510398906937, - 2240577868599239,2240643542273660,2240707464668327,2240769678579424, - 2240830224948918,2240889142947021,2240946470049710,2241002242111632, - 2241056493434688,2241109256832545,2241160563691345,2241210444026824, - 2241258926538069,2241306038658085,2241351806601384,2241396255408737, - 2241439408989263,2241481290159988,2241521920683014,2241561321300414, - 2241599511766981,2241636510880914,2241672336512567,2241707005631317, - 2241740534330669,2241772937851645,2241804230604542,2241834426189118, - 2241863537413270,2241891576310240,2241918554154426,2241944481475803, - 2241969368073032,2241993223025259,2242016054702647,2242037870775672, - 2242058678223188,2242078483339294,2242097291739004,2242115108362739, - 2242131937479636,2242147782689690,2242162646924702,2242176532448058, - 2242189440853303,2242201373061504,2242212329317384,2242222309184204, - 2242231311537365,2242239334556685,2242246375717338,2242252431779384, - 2242257498775863,2242261571999386,2242264645987166,2242266714504423, - 2242267770526080,2242267806216682,2242266812908434,2242264781077261, - 2242261700316790,2242257559310117,2242252345799249,2242246046552055, - 2242238647326588,2242230132832599,2242220486690050,2242209691384432, - 2242197728218658,2242184577261284,2242170217290794,2242154625735654, - 2242137778609814,2242119650443302,2242100214207531,2242079441234882, - 2242057301132111,2242033761687055,2242008788768083,2241982346215658, - 2241954395725333,2241924896721420,2241893806220494,2241861078683807, - 2241826665857576,2241790516600019,2241752576693859,2241712788642894, - 2241671091451056,2241627420382213,2241581706698751,2241533877376746, - 2241483854795259,2241431556397014,2241376894317324,2241319774977796, - 2241260098640839,2241197758920517,2241132642244683,2241064627262631, - 2240993584191722,2240919374095516,2240841848084869,2240760846432212, - 2240676197587764,2240587717084761,2240495206318733,2240398451183547, - 2240297220544145,2240191264522592,2240080312570135,2239964071293311, - 2239842221996510,2239714417896679,2239580280957705,2239439398282173, - 2239291317986176,2239135544468183,2238971532964959,2238798683265249, - 2238616332424332,2238423746288075,2238220109591870,2238004514345197, - 2237775946143192,2237533267957802,2237275200846732,2237000300869931, - 2236706931309079,2236393229029127,2236057063479481,2235695986373225, - 2235307169458838,2234887326941556,2234432617919425,2233938522519742, - 2233399683022654,2232809697779175,2232160850599794,2231443750584617, - 2230646845562145,2229755753817960,2228752329126507,2227613325162477, - 2226308442121145,2224797391720369,2223025347823800,2220915633329775, - 2218357446086993,2215184158448627,2211132412537323,2205758503851011, - 2198248265654920,2186916352102052,2167562552481677,2125549880839429}; - -static double wi[ZIGGURAT_TABLE_SIZE] = - {17367254121656703e-31,9558660348275145e-32,12708704832820278e-32, - 14909740960986864e-32,16658733630346416e-32,18136120809053487e-32, - 1942972015219358e-31,20589500627632916e-32,21646860576118966e-32, - 2262294039150043e-31,23532718913376864e-32,24387234556800803e-32, - 25194879828681465e-32,2596219977196592e-31,26694407473112964e-32, - 2739572968463095e-31,280696460019946e-30,28719058903642897e-32, - 29346417484275224e-32,29953809336344285e-32,30543030006769113e-32, - 3111563633851158e-31,3167298801818414e-31,3221628035016365e-31, - 32746570407564125e-32,33264798116476e-29,337718034169968e-30, - 34268340352771636e-32,34755088731390227e-32,3523266384567022e-31, - 3570162463362898e-31,3616248057128073e-31,36615697529342477e-32, - 3706170277693123e-31,37500889278448874e-32,3793361940125627e-31, - 38360228129389374e-32,3878102586096749e-31,3919630085297984e-31, - 39606321365983254e-32,40011337552278087e-32,4041158312387907e-31, - 4080727683070036e-31,4119862377455137e-31,41585816580575855e-32, - 41969036444492247e-32,4234845407127582e-31,42724230518658345e-32, - 43096517956924877e-32,4346546035489394e-31,4383119410062289e-31, - 4419384856424202e-31,4455354660935343e-31,4491040505860591e-31, - 4526453511835132e-31,45616042766683e-29,4596502910863464e-31, - 4631159070186941e-31,4665581985579899e-31,469978049067346e-30, - 4733763047137822e-31,4767537768070579e-31,4801112439606964e-31, - 4834494540915173e-31,4867691262722585e-31,4900709524503576e-31, - 4933555990446197e-31,4966237084303158e-31,499875900322208e-30, - 5031127730640677e-31,5063349048324261e-31,5095428547615612e-31, - 5127371639960692e-31,5159183566767805e-31,5190869408652579e-31, - 5222434094116442e-31,52538824077020155e-32,5285218997665102e-31, - 5316448383199491e-31,5347574961247755e-31,5378603012928409e-31, - 5409536709607314e-31,5440380118638932e-31,5471137208800966e-31, - 550181185544408e-30,5532407845376661e-31,5562928881503102e-31, - 5593378587232605e-31,5623760510674315e-31,5654078128633358e-31, - 5684334850421336e-31,5714534021493849e-31,5744678926926726e-31, - 5774772794741848e-31,5804818799092685e-31,5834820063319006e-31, - 5864779662879593e-31,589470062817121e-30,5924585947241581e-31, - 5954438568403615e-31,598426140275769e-30,601405732662843e-30, - 6043829183921996e-31,6073579788409578e-31,6103311925942512e-31, - 6133028356604082e-31,6162731816802865e-31,6192425021312213e-31, - 6222110665260248e-31,6251791426074554e-31,6281469965385542e-31, - 6311148930892342e-31,6340830958194888e-31,6370518672595733e-31, - 640021469087503e-30,6429921623041988e-31,645964207406601e-30, - 648937864559066e-30,6519133937633505e-31,6548910550274845e-31, - 6578711085338253e-31,6608538148065851e-31,6638394348791179e-31, - 6668282304612498e-31,6698204641069389e-31,6728163993825439e-31, - 6758163010359885e-31,6788204351671041e-31,681829069399439e-30, - 6848424730538249e-31,6878609173239948e-31,6908846754545526e-31, - 6939140229215998e-31,696949237616333e-30,6999906000319335e-31, - 7030383934540792e-31,7060929041554193e-31,7091544215943653e-31, - 7122232386185626e-31,7152996516734219e-31,7183839610161045e-31, - 7214764709353755e-31,7245774899777502e-31,72768733118038725e-32, - 7308063123111988e-31,7339347561166714e-31,7370729905779203e-31, - 7402213491755235e-31,7433801711637146e-31,7465498018545449e-31, - 7497305929126601e-31,7529229026613742e-31,7561270964007667e-31, - 7593435467385694e-31,7625726339346621e-31,7658147462600412e-31, - 7690702803711903e-31,7723396417008341e-31,7756232448661274e-31, - 778921514095401e-30,7822348836746627e-31,7855637984151357e-31, - 7889087141432085e-31,7922700982142658e-31,7956484300519808e-31, - 7990442017147628e-31,8024579184911813e-31,8058900995263265e-31, - 8093412784812165e-31,812812004227522e-30,8163028415800651e-31, - 8198143720697359e-31,8233471947596931e-31,8269019271079405e-31, - 8304792058796362e-31,834079688112767e-30,8377040521411316e-31, - 8413529986789175e-31,8450272519715296e-31,8487275610177406e-31, - 85245470086869e-29,8562094740097588e-31,8599927118319072e-31, - 86380527619967175e-32,8676480611237092e-31,8715219945465259e-31, - 8754280402508787e-31,8793671999012706e-31,8833405152300122e-31, - 88734907038049e-29,8913939944215902e-31,8954764640486935e-31, - 8995977064883017e-31,9037590026252085e-31,9079616903732087e-31, - 9122071683126914e-31,9164968996211253e-31,9208324163254476e-31, - 9252153239087913e-31,9296473063078686e-31,9341301313417584e-31, - 938665656617903e-30,9432558359669126e-31,9479027264644209e-31, - 95260849610588e-29,957375432209002e-30,962205950628746e-30, - 9671026058815726e-31,972068102289435e-30,9771053062699983e-31, - 9822172599183368e-31,9874071960473548e-31,9926785548800904e-31, - 9980350026176626e-31,10034804521429213e-31,10090190861630543e-31, - 10146553831460223e-31,10203941464676316e-31,1026240537260681e-30, - 10322001115479755e-31,10382788623508751e-31,10444832675993878e-31, - 10508203448348659e-31,1057297713900341e-30,10639236690670377e-31, - 10707072623626628e-31,107765840026618e-29,10847879564397177e-31, - 10921079038143372e-31,109963147017795e-29,11073733224929686e-31, - 11153497865847152e-31,11235791107104895e-31,11320817840158973e-31, - 11408809242576976e-31,1150002753783406e-30,11594771891443527e-31, - 11693385786905373e-31,1179626635295029e-30,11903876299277459e-31, - 1201675939253847e-30,12135560818661637e-31,12261054417445396e-31, - 12394179789158183e-31,12536093926597603e-31,1268824481425016e-30, - 12852479319091384e-31,13031206634685398e-31,13227655770190893e-31, - 13446300925006917e-31,13693606835124475e-31,13979436672771461e-31, - 14319989869657897e-31,14744848603594667e-31,1531787274160907e-30, - 16227698675312968e-31}; - -static double fi[ZIGGURAT_TABLE_SIZE] = - {1.,.9771017012827331,.9598790918124174,.9451989534530794, - .9320600759689914,.9199915050483614,.9087264400605639,.898095921906305, - .8879846607634008,.8783096558161477,.869008688043794,.8600336212030095, - .8513462584651245,.842915653118442,.8347162929929313,.8267268339520951, - .8189291916094156,.8113078743182208,.8038494831763903,.7965423304282554, - .7893761435711993,.7823418326598627,.775431304986139,.7686373158033355, - .7619533468415471,.7553735065117552,.7488924472237273,.7425052963446368, - .7362075981312672,.729995264565803,.7238645334728822,.717811932634902, - .711834248882359,.7059285013367979,.7000919181404905,.694321916130033, - .6886160830085275,.6829721616487918,.6773880362225135,.6718617199007669, - .6663913439123812,.6609751477802419,.6556114705832252,.650298743114295, - .6450354808242524,.6398202774564395,.6346517992909606,.6295287799281287, - .6244500155502747,.6194143606090396,.6144207238920772,.6094680649288958, - .6045553907005499,.599681752622168,.5948462437709915,.5900479963357922, - .5852861792663006,.5805599961036837,.5758686829752109,.5712115067380753, - .5665877632589521,.5619967758172782,.5574378936214867,.5529104904285204, - .5484139632579217,.5439477311926505,.5395112342595453,.5351039323830201, - .5307253044061945,.5263748471741873,.5220520746747954,.5177565172322012, - .5134877207497434,.5092452459981365,.5050286679458292,.5008375751284826, - .49667156905479676,.4925302636461491,.4884132847077125,.48432026942891204, - .48025086591125016,.4762047327216842,.4721815384698837,.46818096140782267, - .4642026890502793,.460246417814924,.45631185268077407,.4523987068638829, - .4485067015092144,.4446355653977281,.4407850346677702,.43695485254992955, - .43314476911457434,.42935454103134185,.42558393133990086,.4218327092313535, - .41810064983968476,.4143875340427069,.41069314827198344,.40701728433124823, - .4033597392228692,.399720314981932,.3960988185175474,.39249506146101104, - .3889088600204649,.38534003484173424,.3817884108750316,.37825381724723833, - .37473608713949164,.37123505766982157,.3677505697805964,.36428246813054976, - .36083060099117586,.3573948201472906,.35397498080156936,.35057094148288126, - .34718256395825153,.3438097131482915,.3404522570459456,.33711006663841303, - .33378301583210873,.3304709813805373,.32717384281495887,.32389148237773235, - .3206237849582305,.3173706380312227,.3141319315976305,.310907558127564, - .307697412505554,.30450139197789644,.30131939610203423,.29815132669790145, - .2949970878011627,.2918565856182811,.28872972848335393,.28561642681665805, - .28251659308484933,.27943014176276515,.2763569892967811,.27329705406967564, - .2702502563669598,.26721651834463167,.2641957639983174,.2611879191337636, - .25819291133864797,.2552106699556771,.25224112605694377,.2492842124195167, - .24633986350223877,.24340801542371202,.24048860594144916,.23758157443217368, - .2346868618732527,.23180441082524855,.22893416541557743,.22607607132326474, - .22323007576478943,.22039612748101145,.21757417672517823,.2147641752520084, - .21196607630785277,.20917983462193548,.20640540639867916,.20364274931112133, - .20089182249543117,.19815258654653795,.1954250035148854,.19270903690432864, - .19000465167119293,.18731181422451676,.18463049242750437,.18196065560021638, - .17930227452353026,.17665532144440646,.17401977008249914,.17139559563815535, - .16878277480185,.1661812857651097,.1635911082329826,.16101222343811727, - .1584446141565199,.15588826472506426,.15334316106083742,.15080929068240986, - .1482866427331284,.14577520800653765,.14327497897404687,.14078594981496803, - .138308116449064,.13584147657175705,.13338602969216254,.13094177717412792, - .12850872228047336,.12608687022065,.12367622820205106,.12127680548523516, - .11888861344334545,.11651166562603685,.11414597782825504,.1117915681642454, - .10944845714720981,.10711666777507266,.10479622562286683,.10248715894230599, - .10018949876917177,.09790327903921535,.09562853671335306,.09336531191302634, - .09111364806670041,.08887359206859394,.08664519445086755,.08442850957065445, - .0822235958134955,.08003051581494733,.07784933670237201,.07568013035919481, - .07352297371424082,.07137794905914183,.06924514439725017,.06712465382802392, - .06501657797147035,.06292102443797778,.06083810834975175,.05876795292113793, - .056710690106399425,.05466646132507786,.05263541827697361,.05061772386112175, - .04861355321603513,.04662309490208967,.044646552251446515,.042684144916619336, - .04073611065607874,.0388027074046569,.03688421568869112,.034980941461833046, - .033093219458688684,.031221417192023686,.02936593975823011,.027527235669693315, - .02570580400863265,.023902203305873237,.022117062707379908,.020351096230109344, - .018605121275783343,.016880083152595836,.015177088307982065,.013497450601780796, - .0118427578579431,.0102149714397311,.008616582769422912,.007050875471392109, - .005522403299264755,.0040379725933718715,.002609072746106362,.0012602859304985975}; - -// Tables for exprnd -static ZIGINT ke[ZIGGURAT_TABLE_SIZE] = - {3985772928715748, 0,2742928985168065,3438700186803721, - 3744780257810519,3914896975372863,4022625697542798,4096776410635450, - 4150853606149210,4192001604687417,4224344877584101,4250427292531740, - 4271901371161554,4289886428824118,4305167164135199,4318309783140431, - 4329732973408940,4339752937704679,4348612900760388,4356502988721768, - 4363573953227346,4369946852445020,4375720012348349,4380974119031481, - 4385776001930298,4390181484145305,4394237557465219,4397984061535398, - 4401454994146430,4404679543790856,4407682910787985,4410486965794400, - 4413110782053579,4415571068741702,4417882526198713,4420058138987325, - 4422109419110772,4424046609003130,4425878851844253,4427614335173868, - 4429260412563040,4430823707156475,4432310200160197,4433725306767517, - 4435073941555377,4436360575016074,4437589282595121,4438763787369085, - 4439887497305303,4440963537889317,4441994780778252,4442983869033585, - 4443933239400428,4444845142028910,4445721657973833,4446564714759241, - 4447376100252993,4448157475061632,4448910383626429,4449636264176642, - 4450336457674983,4451012215872352,4451664708573597,4452295030203006, - 4452904205747010,4453493196141906,4454062903166143,4454614173889474, - 4455147804725090,4455664545125435,4456165100957688,4456650137590828, - 4457120282722585,4457576128971459,4458018236256245,4458447133983073, - 4458863323057847,4459267277740095,4459659447352586,4460040257859578, - 4460410113325310,4460769397263133,4461118473884710,4461457689257740, - 4461787372379910,4462107836175980,4462419378424319,4462722282618581, - 4463016818769709,4463303244152965,4463581804004301,4463852732169940, - 4464116251712773,4464372575478779,4464621906626490,4464864439122178, - 4465100358203284,4465329840812355,4465553056003596,4465770165323939, - 4465981323170417,4466186677125455,4466386368271563,4466580531486827, - 4466769295722448,4466952784263502,4467131114974006,4467304400527265, - 4467472748622447,4467636262188208,4467795039574164,4467949174730939, - 4468098757379442,4468243873170018,4468384603832024,4468521027314373, - 4468653217917530,4468781246417428,4468905180181701,4469025083278642, - 4469141016579234,4469253037852582,4469361201855066,4469465560413474, - 4469566162502383,4469663054316032,4469756279334881,4469845878387080, - 4469931889704995,4470014348976986,4470093289394551,4470168741694984, - 4470240734199652,4470309292847996,4470374441227332,4470436200598525, - 4470494589917605,4470549625853344,4470601322800852,4470649692891185, - 4470694745996980,4470736489734116,4470774929459349,4470810068263924, - 4470841906963074,4470870444081369,4470895675833821,4470917596102651, - 4470936196409614,4470951465883737,4470963391224346,4470971956659198, - 4470977143897542,4470978932077904,4470977297710362,4470972214613072, - 4470963653842747,4470951583618802,4470935969240827,4470916772999009, - 4470893954077117,4470867468447603,4470837268758338,4470803304210460, - 4470765520426769,4470723859310029,4470678258890503,4470628653161980, - 4470574971905457,4470517140499614,4470455079717082,4470388705505446, - 4470317928751818,4470242655029689,4470162784326669,4470078210751556, - 4469988822219058,4469894500110287,4469795118907000,4469690545797298, - 4469580640250319,4469465253557163,4469344228335006,4469217397991048, - 4469084586142556,4468945605988875,4468800259630802,4468648337332217, - 4468489616718259,4468323861903709,4468150822544456,4467970232804102, - 4467781810226787,4467585254506222,4467380246139658,4467166444954116, - 4466943488490515,4466710990229518,4466468537640691,4466215690034133, - 4465951976190801,4465676891744455,4465389896284247,4465090410142477, - 4464777810826750,4464451429049612,4464110544301482,4463754379904174, - 4463382097472202,4462992790697122,4462585478355953,4462159096427753, - 4461712489182116,4461244399078944,4460753455289386,4460238160612098, - 4459696876515553,4459127805983956,4458528973779075,4457898203649722, - 4457233091920646,4456530976767892,4455788902331217,4455003576616607, - 4454171321891082,4453288015951104,4452349022232651,4451349106194827, - 4450282334707462,4449141954247903,4447920242480611,4446608326137821, - 4445195955871677,4443671225661690,4442020220072463,4440226566619900, - 4438270861888260,4436129927556552,4433775834104270,4431174602388627, - 4428284451100006,4425053392146958,4421415870372502,4417287970124084, - 4412560416174562,4407088078325945,4400673742272494,4393042098597073, - 4383796248451589,4372341169422858,4357740343059956,4338425130125967, - 4311541827049177,4271262897902398,4203411844498905,4061213381260384}; - -static double we[ZIGGURAT_TABLE_SIZE] = - {19311480126418366e-31,1417802848791084e-32,23278824993382457e-33, - 30487830247064326e-33,3666569771447489e-32,4217930218928974e-32, - 4722256155686277e-32,51911915446217885e-33,5632347108395505e-32, - 6051008260642765e-32,645101650967275e-31,6835264680370054e-32, - 7205993957468906e-32,7564981553739299e-32,7913664396195108e-32, - 8253223556351894e-32,8584643616885051e-32,8908755486564743e-32, - 9226267962966373e-32,9537791450529272e-32,9843856087455926e-32, - 10144925809006294e-32,10441409405585343e-32,10733669323436384e-32, - 1102202874567019e-31,11306777346479334e-32,11588176009705533e-32, - 11866460730417886e-32,1214184586569436e-31,12414526862326387e-32, - 12684682560606153e-32,12952477151912284e-32,1321806185153881e-31, - 13481576335745447e-32,13743149982367625e-32,14002902946807862e-32, - 14260947099321287e-32,14517386844829297e-32,14772319842763584e-32, - 15025837641447456e-32,15278026239101652e-32,15528966581595696e-32, - 1577873500545958e-31,1602740363335091e-31,16275040728083524e-32, - 16521711010420076e-32,16767475945078279e-32,17012393998770646e-32, - 17256520873568226e-32,17499909718432365e-32,17742611321380505e-32, - 17984674284430714e-32,18226145183195818e-32,18467068712763576e-32, - 18707487821298258e-32,18947443832625902e-32,19186976558915997e-32, - 19426124404443042e-32,19664924461299023e-32,19903412597830144e-32, - 20141623540485899e-32,20379590949693882e-32,2061734749030844e-31, - 2085492489712377e-31,21092354035891528e-32,21329664960238294e-32, - 21566886964838972e-32,2180404863516701e-31,22041177894111562e-32, - 2227830204572395e-31,2251544781633135e-31,22752641393233694e-32, - 22989908461180186e-32,23227274236804366e-32,23464763501180916e-32, - 2370240063065339e-31,23940209626069303e-32,2417821414054771e-31, - 24416437505894123e-32,24654902757768304e-32,2489363265970225e-31, - 2513264972605797e-31,2537197624400795e-31,2561163429461499e-31, - 2585164577308239e-31,26092032408240577e-32,2633281578133145e-31, - 2657401734414762e-31,2681565843657999e-31,2705776030362351e-31, - 27300344111887955e-32,27543430965657624e-32,2778704192254128e-31, - 2803119800875143e-31,28275920234049704e-32,2852122960639331e-31, - 28767147146315804e-32,29013693901073754e-32,29260890958589514e-32, - 29508759461219033e-32,2975732061937252e-31,3000659572501474e-31, - 3025660616507079e-31,3050737343476251e-31,3075891915089994e-31, - 31011265065151543e-32,3126443307731675e-31,31518445248623523e-32, - 31773323815073683e-32,32029091200858335e-32,32285770031865573e-32, - 3254338314930261e-31,3280195362345436e-31,3306150476760074e-31, - 3332206015211484e-31,33583643618764577e-32,33846279295240445e-32, - 34109991609932597e-32,34374805306980633e-32,34640745461620167e-32, - 3490783749585068e-31,3517610719444983e-31,3544558072136013e-31, - 3571628463647465e-31,35988245912849274e-32,3626149195437003e-31, - 36536050613905045e-32,36811950211971757e-32,3708921955595139e-31, - 37367887959883854e-32,3764798526487784e-31,37929541860172334e-32, - 3821258870488753e-31,38497157350504876e-32,3878327996411799e-31, - 39070989352498183e-32,3936031898702075e-31,3965130302950038e-31, - 3994397635898684e-31,40238374599574693e-32,40534534149283966e-32, - 4083249221007178e-31,41132286819038357e-32,4143395688089474e-31, - 417375422017632e-30,42043083524385856e-32,4235062256482152e-31, - 4266020205071558e-31,42971865761233266e-32,43285658568752094e-32, - 4360162648241568e-31,43919816693657415e-32,4424027762380992e-31, - 4456305897392361e-31,4488821177692617e-31,4521578845226347e-31, - 4554584286317242e-31,4587843037674623e-31,4621360792696427e-31, - 4655143408087069e-31,4689196910809916e-31,4723527505395548e-31, - 4758141581628553e-31,4793045722637247e-31,4828246713412587e-31, - 4863751549784512e-31,489956744788614e-30,4935701854138577e-31, - 4972162455791703e-31,5008957192059114e-31,5046094265888434e-31, - 5083582156411624e-31,5121429632123542e-31,5159645764841062e-31, - 5198239944499494e-31,5237221894847848e-31,5276601690109886e-31, - 531638977268369e-30,535659697195905e-30,5397234524338979e-31, - 5438314094559637e-31,547984779841163e-30,5521848226975234e-31, - 5564328472492872e-31,5607302156013967e-31,5650783456960506e-31, - 5694787144776348e-31,5739328612839635e-31,5784423914835991e-31, - 5830089803810586e-31,5876343774140057e-31,5923204106690931e-31, - 5970689917460091e-31,6018821210025236e-31,6067618932170007e-31, - 6117105037089722e-31,616730254963062e-30,6218235638068533e-31, - 6269929691993326e-31,6322411406934211e-31,6375708876439426e-31, - 6429851692413595e-31,6484871054618903e-31,6540799890364481e-31, - 6597672985544566e-31,6655527128343343e-31,6714401267106488e-31, - 677433668409101e-30,6835377187051274e-31,6897569320906848e-31, - 6960962602074885e-31,7025609778445959e-31,7091567118449584e-31, - 7158894733208553e-31,7227656936438121e-31,7297922647529085e-31, - 7369765844191243e-31,7443266072160415e-31,7518509020832513e-31, - 7595587175337749e-31,7674600557578427e-31,7755657571215791e-31, - 7838875968622858e-31,792438396157355e-30,8012321502113083e-31, - 8102841765913146e-31,8196112877806125e-31,8292319928581809e-31, - 8391667344146798e-31,849438168364877e-30,8600714963334941e-31, - 8710948629387904e-31,882539833807214e-30,8944419748519865e-31, - 9068415597131669e-31,9197844409811865e-31,9333231329422952e-31, - 9475181706524984e-31,9624398345658476e-31,978170365478442e-30, - 994806847238388e-30,1012465014428832e-30,10312843657756166e-31, - 1051435160404455e-30,10731281954224043e-31,10966288068517408e-31, - 1122277490935032e-30,11505212963006663e-31,11819635283304206e-31, - 12174462832361815e-31,12581958069755114e-31,13060984107128082e-31, - 13642786158057857e-31,14384889932178723e-31,15412190700064194e-31, - 17091034077168055e-31}; - -static double fe[ZIGGURAT_TABLE_SIZE] = - { 1.0,.9381436808621746,.9004699299257464,.8717043323812036, - .8477855006239896,.8269932966430503,.8084216515230084,.7915276369724956, - .7759568520401156,.7614633888498963,.7478686219851951,.7350380924314235, - .722867659593572,.711274760805076,.7001926550827882,.689566496117078, - .6793505722647654,.6695063167319247,.6600008410789997,.650805833414571, - .6418967164272661,.6332519942143661,.6248527387036659,.6166821809152077, - .608725382079622,.6009689663652322,.5934009016917334,.586010318477268, - .578787358602845,.5717230486648258,.5648091929124002,.5580382822625874, - .5514034165406413,.5448982376724396,.5385168720028619,.5322538802630432, - .5261042139836197,.5200631773682336,.5141263938147486,.5082897764106429, - .5025495018413477,.49690198724154955,.49134386959403253,.4858719873418849, - .4804833639304542,.4751751930373774,.46994482528396,.4647897562504262, - .4597076156421377,.45469615747461545,.449753251162755,.4448768734145485, - .4400651008423539,.4353161032156366,.43062813728845883,.42599954114303434, - .4214287289976166,.4169141864330029,.4124544659971612,.4080481831520324, - .4036940125305303,.3993906844752311,.39513698183329016,.3909317369847971, - .38677382908413765,.38266218149600983,.3785957594095808,.37457356761590216, - .370594648435146,.36665807978151416,.3627629733548178,.3589084729487498, - .35509375286678746,.35131801643748334,.347580494621637,.3438804447045024, - .34021714906678,.33658991402867755,.332998068761809,.3294409642641363, - .3259179723935562,.3224284849560891,.31897191284495724,.31554768522712895, - .31215524877417955,.3087940669345602,.30546361924459026,.3021634006756935, - .2988929210155818,.2956517042812612,.2924392881618926,.28925522348967775, - .2860990737370768,.28297041453878075,.27986883323697287,.27679392844851736, - .27374530965280297,.27072259679906,.2677254199320448,.2647534188350622, - .2618062426893629,.25888354974901623,.2559850070304154,.25311029001562946, - .2502590823688623,.24743107566532763,.2446259691318921,.24184346939887721, - .23908329026244918,.23634515245705964,.23362878343743335,.2309339171696274, - .2282602939307167,.22560766011668407,.22297576805812017,.2203643758433595, - .21777324714870053,.21520215107537868,.21265086199297828,.21011915938898826, - .20760682772422204,.2051136562938377,.20263943909370902,.20018397469191127, - .19774706610509887,.19532852067956322,.19292814997677132,.1905457696631954, - .1881811994042543,.1858342627621971,.18350478709776746,.1811926034754963, - .1788975465724783,.17661945459049488,.1743581691713535,.17211353531532006, - .16988540130252766,.1676736186172502,.165478041874936,.16329852875190182, - .16113493991759203,.1589871389693142,.15685499236936523,.15473836938446808, - .15263714202744286,.1505511850010399,.1484803756438668,.14642459387834494, - .14438372216063478,.1423576454324722,.14034625107486245,.1383494288635802, - .13636707092642886,.13439907170221363,.13244532790138752,.13050573846833077, - .12858020454522817,.12666862943751067,.12477091858083096,.12288697950954514, - .12101672182667483,.11916005717532768,.11731689921155557,.11548716357863353, - .11367076788274431,.1118676316700563,.11007767640518538,.1083008254510338, - .10653700405000166,.10478613930657017,.10304816017125772,.10132299742595363, - .09961058367063713,.0979108533114922,.0962237425504328,.09454918937605586, - .09288713355604354,.09123751663104016,.08960028191003286,.08797537446727022, - .08636274114075691,.08476233053236812,.08317409300963238,.08159798070923742, - .0800339475423199,.07848194920160642,.0769419431704805,.07541388873405841, - .07389774699236475,.07239348087570874,.07090105516237183,.06942043649872875, - .0679515934219366,.06649449638533977,.06504911778675375,.06361543199980733, - .062193415408540995,.06078304644547963,.059384305633420266,.05799717563120066, - .05662164128374288,.05525768967669704,.05390531019604609,.05256449459307169, - .05123523705512628,.04991753428270637,.0486113855733795,.04731679291318155, - .04603376107617517,.04476229773294328,.04350241356888818,.042254122413316234, - .04101744138041482,.039792391023374125,.03857899550307486,.03737728277295936, - .03618728478193142,.03500903769739741,.03384258215087433,.032687963508959535, - .03154523217289361,.030414443910466604,.029295660224637393,.028188948763978636, - .0270943837809558,.026012046645134217,.024942026419731783,.02388442051155817, - .02283933540638524,.02180688750428358,.020787204072578117,.019780424338009743, - .01878670074469603,.01780620041091136,.016839106826039948,.015885621839973163, - .014945968011691148,.014020391403181938,.013109164931254991,.012212592426255381, - .011331013597834597,.010464810181029979,.00961441364250221,.008780314985808975, - .00796307743801704,.007163353183634984,.006381905937319179,.005619642207205483, - .004877655983542392,.004157295120833795,.003460264777836904,.002788798793574076, - .0021459677437189063,.0015362997803015724,.0009672692823271745,.00045413435384149677}; - - -/* - * Here is the guts of the algorithm. As Marsaglia and Tsang state the - * algorithm in their paper - * - * 1) Calculate a random signed integer j and let i be the index - * provided by the rightmost 8-bits of j - * 2) Set x = j * w_i. If j < k_i return x - * 3) If i = 0, then return x from the tail - * 4) If [f(x_{i-1}) - f(x_i)] * U < f(x) - f(x_i), return x - * 5) goto step 1 - * - * Where f is the functional form of the distribution, which for a normal - * distribution is exp(-0.5*x*x) - */ - -/* NOTE: This is identical to randmtzig_gv_randn() below except for the random number generation */ -double randmtzig_randn (dsfmt_t *dsfmt) -{ - while (1) - { - /* arbitrary mantissa (selected by randi, with 1 bit for sign) */ - const randmtzig_uint64_t r = randi(dsfmt); - const randmtzig_int64_t rabs=r>>1; - const int idx = (int)(rabs&0xFF); - const double x = ( r&1 ? -rabs : rabs) * wi[idx]; - - if (rabs < (randmtzig_int64_t)ki[idx]) { - return x; /* 99.3% of the time we return here 1st try */ - } else if (idx == 0) { - /* As stated in Marsaglia and Tsang - * - * For the normal tail, the method of Marsaglia[5] provides: - * generate x = -ln(U_1)/r, y = -ln(U_2), until y+y > x*x, - * then return r+x. Except that r+x is always in the positive - * tail!!!! Any thing random might be used to determine the - * sign, but as we already have r we might as well use it - * - * [PAK] but not the bottom 8 bits, since they are all 0 here! - */ - double xx, yy; - do { - xx = - ZIGGURAT_NOR_INV_R * log (randu(dsfmt)); - yy = - log (randu(dsfmt)); - } - while ( yy+yy <= xx*xx); - return (rabs&0x100 ? -ZIGGURAT_NOR_R-xx : ZIGGURAT_NOR_R+xx); - } else if ((fi[idx-1] - fi[idx]) * randu(dsfmt) + fi[idx] < exp(-0.5*x*x)) { - return x; - } - - } -} - -/* NOTE: This is identical to randmtzig_randn() above except for the random number generation */ -double randmtzig_gv_randn (void) -{ - while (1) - { - /* arbitrary mantissa (selected by NRANDI, with 1 bit for sign) */ - const randmtzig_uint64_t r = NRANDI; - const randmtzig_int64_t rabs=r>>1; - const int idx = (int)(rabs&0xFF); - const double x = ( r&1 ? -rabs : rabs) * wi[idx]; - - if (rabs < (randmtzig_int64_t)ki[idx]) { - return x; /* 99.3% of the time we return here 1st try */ - } else if (idx == 0) { - /* As stated in Marsaglia and Tsang - * - * For the normal tail, the method of Marsaglia[5] provides: - * generate x = -ln(U_1)/r, y = -ln(U_2), until y+y > x*x, - * then return r+x. Except that r+x is always in the positive - * tail!!!! Any thing random might be used to determine the - * sign, but as we already have r we might as well use it - * - * [PAK] but not the bottom 8 bits, since they are all 0 here! - */ - double xx, yy; - do { - xx = - ZIGGURAT_NOR_INV_R * log (RANDU); - yy = - log (RANDU); - } - while ( yy+yy <= xx*xx); - return (rabs&0x100 ? -ZIGGURAT_NOR_R-xx : ZIGGURAT_NOR_R+xx); - } else if ((fi[idx-1] - fi[idx]) * RANDU + fi[idx] < exp(-0.5*x*x)) { - return x; - } - - } -} - -double randmtzig_exprnd (void) -{ - while (1) - { - ZIGINT ri = ERANDI; - const int idx = (int)(ri & 0xFF); - const double x = ri * we[idx]; - if (ri < ke[idx]) - return x; // 98.9% of the time we return here 1st try - else if (idx == 0) - { - /* As stated in Marsaglia and Tsang - * - * For the exponential tail, the method of Marsaglia[5] provides: - * x = r - ln(U); - */ - return ZIGGURAT_EXP_R - log(RANDU); - } - else if ((fe[idx-1] - fe[idx]) * RANDU + fe[idx] < exp(-x)) - return x; - } -} - -void randmtzig_gv_fill_randn (double *p, randmtzig_idx_type n) -{ - randmtzig_idx_type i; - for (i = 0; i < n; i++) - p[i] = randmtzig_gv_randn(); -} - -void randmtzig_fill_randn (dsfmt_t *dsfmt, double *p, randmtzig_idx_type n) -{ - randmtzig_idx_type i; - for (i = 0; i < n; i++) - p[i] = randmtzig_randn(dsfmt); -} - - -void randmtzig_fill_exprnd (double *p, randmtzig_idx_type n) -{ - randmtzig_idx_type i; - for (i = 0; i < n; i++) - p[i] = randmtzig_exprnd(); -} - - -#ifdef STANDALONE - -int main(int ac, char *av[]) { - if (ac == 1) { - printf("Usage: randmtzig \n"); - return (-1); - } - - int n = atoi(av[1]); - time_t t1; - - dsfmt_gv_init_gen_rand(0); - - double *p; posix_memalign((void **)&p, 16, n*sizeof(double)); - uint32_t *u; posix_memalign((void **)&u, 16, 2*n*sizeof(uint32_t)); - - t1 = clock(); - dsfmt_gv_fill_array_close_open(p, n); - printf("Uniform fill (n): %f\n", (clock() - t1) / (double) CLOCKS_PER_SEC); - - t1 = clock(); - for (int i = 0; i < n; i++) p[i] = dsfmt_gv_genrand_close_open(); - printf("Uniform (n): %f\n", (clock() - t1) / (double) CLOCKS_PER_SEC); - - t1 = clock(); - for (int i = 0; i < 2*n; i++) u[i] = dsfmt_gv_genrand_uint32(); - printf("Uniform 32-bit ints (2*n): %f\n", (clock() - t1) / (double) CLOCKS_PER_SEC); - - memset((void *)p, 0, n*sizeof(double)); - t1 = clock(); - for (int i = 0; i < n; i++) p[i] = randmtzig_gv_randn(); - printf("Normal (n): %f\n", (clock() - t1) / (double) CLOCKS_PER_SEC); - for (int i = 0; i < 10; i++) printf("%lf\n", p[i]); - - return 0; -} - -#endif diff --git a/test/random.jl b/test/random.jl index dbc3d4bd5ba3e..4327211128f05 100644 --- a/test/random.jl +++ b/test/random.jl @@ -16,11 +16,11 @@ @test rand(MersenneTwister(5294967296)) == 0.3498809918210497 # randn -@test randn(MersenneTwister(42)) == -0.5560268761438381 +@test randn(MersenneTwister(42)) == -0.5560268761463861 A = zeros(2, 2) randn!(MersenneTwister(42), A) -@test A == [-0.5560268761438381 0.02715533800914659; - -0.4443833571072952 -0.29948409035852047] +@test A == [-0.5560268761463861 0.027155338009193845; + -0.444383357109696 -0.29948409035891055] for T in (Int8, Uint8, Int16, Uint16, Int32, Uint32, Int64, Uint64, Int128, Uint128, Char, Float16, Float32, Float64, Rational{Int}) @@ -38,3 +38,86 @@ if sizeof(Int32) < sizeof(Int) @test typeof(r) == Int32 @test -1 <= r <= typemax(Int32) end + +# Test ziggurat tables +ziggurat_table_size = 256 +nmantissa = 2^51 # one bit for the sign +ziggurat_nor_r = BigFloat("3.65415288536100879635194725185604664812733315920964488827246397029393565706474") +nor_section_area = ziggurat_nor_r*exp(-ziggurat_nor_r^2/2) + erfc(ziggurat_nor_r/sqrt(BigFloat(2)))*sqrt(big(π)/2) +emantissa = 2^52 +ziggurat_exp_r = BigFloat("7.69711747013104971404462804811408952334296818528283253278834867283241051210533") +exp_section_area = (ziggurat_exp_r + 1)*exp(-ziggurat_exp_r) + +const ki = Array(Uint64, ziggurat_table_size) +const wi = Array(Float64, ziggurat_table_size) +const fi = Array(Float64, ziggurat_table_size) +# Tables for exponential variates +const ke = Array(Uint64, ziggurat_table_size) +const we = Array(Float64, ziggurat_table_size) +const fe = Array(Float64, ziggurat_table_size) +function randmtzig_fill_ziggurat_tables() # Operates on the global arrays + wib = big(wi) + fib = big(fi) + web = big(we) + feb = big(fe) + # Ziggurat tables for the normal distribution + x1 = ziggurat_nor_r + wib[256] = x1/nmantissa + fib[256] = exp(-0.5*x1*x1) + # Index zero is special for tail strip, where Marsaglia and Tsang + # defines this as + # k_0 = 2^31 * r * f(r) / v, w_0 = 0.5^31 * v / f(r), f_0 = 1, + # where v is the area of each strip of the ziggurat. + ki[1] = uint(itrunc(x1*fib[256]/nor_section_area*nmantissa)) + wib[1] = nor_section_area/fib[256]/nmantissa + fib[1] = one(BigFloat) + + for i = 255:-1:2 + # New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus + # need inverse operator of y = exp(-0.5*x*x) -> x = sqrt(-2*ln(y)) + x = sqrt(-2.0*log(nor_section_area/x1 + fib[i+1])) + ki[i+1] = uint64(itrunc(x/x1*nmantissa)) + wib[i] = x/nmantissa + fib[i] = exp(-0.5*x*x) + x1 = x + end + + ki[2] = uint64(0) + + # Zigurrat tables for the exponential distribution + x1 = ziggurat_exp_r + web[256] = x1/emantissa + feb[256] = exp(-x1) + + # Index zero is special for tail strip, where Marsaglia and Tsang + # defines this as + # k_0 = 2^32 * r * f(r) / v, w_0 = 0.5^32 * v / f(r), f_0 = 1, + # where v is the area of each strip of the ziggurat. + ke[1] = uint64(itrunc(x1*feb[256]/exp_section_area*emantissa)) + web[1] = exp_section_area/feb[256]/emantissa + feb[1] = one(BigFloat) + + for i = 255:-1:2 + # New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus + # need inverse operator of y = exp(-x) -> x = -ln(y) + x = -log(exp_section_area/x1 + feb[i+1]) + ke[i+1] = uint64(itrunc(x/x1*emantissa)) + web[i] = x/emantissa + feb[i] = exp(-x) + x1 = x + end + ke[2] = zero(Uint64) + + wi[:] = float64(wib) + fi[:] = float64(fib) + we[:] = float64(web) + fe[:] = float64(feb) + return nothing +end +randmtzig_fill_ziggurat_tables() +@test all(ki == Base.LibRandom.ki) +@test all(wi == Base.LibRandom.wi) +@test all(fi == Base.LibRandom.fi) +@test all(ke == Base.LibRandom.ke) +@test all(we == Base.LibRandom.we) +@test all(fe == Base.LibRandom.fe)