From 1e2283de433f6a3dff85b8d16ad741f37651eb81 Mon Sep 17 00:00:00 2001 From: Huon Wilson Date: Thu, 31 Oct 2013 23:27:26 +1100 Subject: [PATCH] std::rand: correct an off-by-one in the Ziggurat code. The code was using (in the notation of Doornik 2005) `f(x_{i+1}) - f(x_{i+2})` rather than `f(x_i) - f(x_{i+1})`. This corrects that, and removes the F_DIFF tables which caused this problem in the first place. They `F_DIFF` tables are a micro-optimisation (in theory, they could easily be a micro-pessimisation): that `if` gets hit about 1% of the time for Exp/Normal, and the rest of the condition involves RNG calls and a floating point `exp`, so it is unlikely that saving a single FP subtraction will be very useful (especially as more tables means more memory reads and higher cache pressure, as well as taking up space in the binary (although only ~2k in this case)). Closes #10084. Notably, unlike that issue suggests, this wasn't a problem with the Exp tables. It affected Normal too, but since it is symmetric, there was no bias in the mean (as the bias was equal on the positive and negative sides and so cancelled out) but it was visible as a variance slightly lower than it should be. --- src/etc/ziggurat_tables.py | 7 +- src/libstd/rand/distributions.rs | 7 +- src/libstd/rand/ziggurat_tables.rs | 132 ----------------------------- 3 files changed, 5 insertions(+), 141 deletions(-) diff --git a/src/etc/ziggurat_tables.py b/src/etc/ziggurat_tables.py index 51c0da39bd5..263c76a88dd 100755 --- a/src/etc/ziggurat_tables.py +++ b/src/etc/ziggurat_tables.py @@ -19,7 +19,7 @@ from math import exp, sqrt, log, floor import random # The order should match the return value of `tables` -TABLE_NAMES = ['X', 'F', 'F_DIFF'] +TABLE_NAMES = ['X', 'F'] # The actual length of the table is 1 more, to stop # index-out-of-bounds errors. This should match the bitwise operation @@ -43,13 +43,10 @@ def tables(r, v, f, f_inv): # cache the f's fvec = [0]*(TABLE_LEN+1) - fdiff = [0]*(TABLE_LEN+1) for i in range(TABLE_LEN+1): fvec[i] = f(xvec[i]) - if i > 0: - fdiff[i] = fvec[i] - fvec[i-1] - return xvec, fvec, fdiff + return xvec, fvec # Distributions # N(0, 1) diff --git a/src/libstd/rand/distributions.rs b/src/libstd/rand/distributions.rs index e7bcf8ce5d3..a6fa1a2bd87 100644 --- a/src/libstd/rand/distributions.rs +++ b/src/libstd/rand/distributions.rs @@ -207,7 +207,6 @@ fn ziggurat(rng: &mut R, symmetric: bool, X: ziggurat_tables::ZigTable, F: ziggurat_tables::ZigTable, - F_DIFF: ziggurat_tables::ZigTable, pdf: &'static fn(f64) -> f64, zero_case: &'static fn(&mut R, f64) -> f64) -> f64 { static SCALE: f64 = (1u64 << 53) as f64; @@ -237,7 +236,7 @@ fn ziggurat(rng: &mut R, return zero_case(rng, u); } // algebraically equivalent to f1 + DRanU()*(f0 - f1) < 1 - if F[i + 1] + F_DIFF[i + 1] * rng.gen() < pdf(x) { + if F[i + 1] + (F[i] - F[i + 1]) * rng.gen() < pdf(x) { return x; } } @@ -288,7 +287,7 @@ impl Rand for StandardNormal { rng, true, // this is symmetric &ziggurat_tables::ZIG_NORM_X, - &ziggurat_tables::ZIG_NORM_F, &ziggurat_tables::ZIG_NORM_F_DIFF, + &ziggurat_tables::ZIG_NORM_F, pdf, zero_case)) } } @@ -366,7 +365,7 @@ impl Rand for Exp1 { Exp1(ziggurat(rng, false, &ziggurat_tables::ZIG_EXP_X, - &ziggurat_tables::ZIG_EXP_F, &ziggurat_tables::ZIG_EXP_F_DIFF, + &ziggurat_tables::ZIG_EXP_F, pdf, zero_case)) } } diff --git a/src/libstd/rand/ziggurat_tables.rs b/src/libstd/rand/ziggurat_tables.rs index aca2457cac4..049ef3dbb59 100644 --- a/src/libstd/rand/ziggurat_tables.rs +++ b/src/libstd/rand/ziggurat_tables.rs @@ -145,72 +145,6 @@ pub static ZIG_NORM_F: [f64, .. 257] = 0.887984660763399880, 0.898095921906304051, 0.908726440060562912, 0.919991505048360247, 0.932060075968990209, 0.945198953453078028, 0.959879091812415930, 0.977101701282731328, 1.000000000000000000]; -pub static ZIG_NORM_F_DIFF: [f64, .. 257] = - [0.000000000000000000, 0.000782818165911943, 0.001348786815607765, 0.001428899847265509, - 0.001484430705892882, 0.001528472172127356, 0.001565707298030807, 0.001598388670308183, - 0.001627786418212004, 0.001654692743837703, 0.001679637706201265, 0.001702994844613767, - 0.001725038123187510, 0.001745974954326004, 0.001765966477270568, 0.001785140598493315, - 0.001803600702759419, 0.001821431661060659, 0.001838704088536796, 0.001855477433793579, - 0.001871802266665008, 0.001887722003144375, 0.001903274226858077, 0.001918491715965767, - 0.001933403251421835, 0.001948034260540625, 0.001962407334827158, 0.001976542650643127, - 0.001990458313945481, 0.002004170645086643, 0.002017694415851860, 0.002031043048104267, - 0.002044228781321551, 0.002057262814738517, 0.002070155428613822, 0.002082916088226049, - 0.002095553533492583, 0.002108075856553551, 0.002120490569226280, 0.002132804661891696, - 0.002145024655099026, 0.002157156644953973, 0.002169206343177243, 0.002181179112575302, - 0.002193079998548175, 0.002204913757158977, 0.002216684880213121, 0.002228397617726446, - 0.002240055998106505, 0.002251663846325885, 0.002263224800326716, 0.002274742325862292, - 0.002286219729956393, 0.002297660173134250, 0.002309066680560787, 0.002320442152205823, - 0.002331789372137141, 0.002343111017035562, 0.002354409664009627, 0.002365687797781804, - 0.002376947817308683, 0.002388192041889739, 0.002399422716815966, 0.002410642018598946, - 0.002421852059823287, 0.002433054893654529, 0.002444252518034679, 0.002455446879594508, - 0.002466639877306970, 0.002477833365903986, 0.002489029159078809, 0.002500229032490808, - 0.002511434726590794, 0.002522647949281448, 0.002533870378427505, 0.002545103664226889, - 0.002556349431455662, 0.002567609281597438, 0.002578884794865288, 0.002590177532127119, - 0.002601489036740262, 0.002612820836305291, 0.002624174444343735, 0.002635551361907296, - 0.002646953079123743, 0.002658381076686089, 0.002669836827288052, 0.002681321797012387, - 0.002692837446676144, 0.002704385233135737, 0.002715966610556786, 0.002727583031652520, - 0.002739235948893221, 0.002750926815690169, 0.002762657087557796, 0.002774428223256353, - 0.002786241685917290, 0.002798098944155558, 0.002810001473169871, 0.002821950755833219, - 0.002833948283778004, 0.002845995558475284, 0.002858094092312607, 0.002870245409671041, - 0.002882451048004164, 0.002894712558920987, 0.002907031509275432, 0.002919409482262880, - 0.002931848078526783, 0.002944348917277934, 0.002956913637427061, 0.002969543898733384, - 0.002982241382970874, 0.002995007795115689, 0.003007844864553855, 0.003020754346314269, - 0.003033738022328147, 0.003046797702715820, 0.003059935227105459, 0.003073152465984053, - 0.003086451322084072, 0.003099833731808721, 0.003113301666695822, 0.003126857134927052, - 0.003140502182881588, 0.003154238896738770, 0.003168069404132778, 0.003181995875862154, - 0.003196020527657495, 0.003210145622009941, 0.003224373470066433, 0.003238706433592253, - 0.003253146927007733, 0.003267697419501892, 0.003282360437226572, 0.003297138565578506, - 0.003312034451571411, 0.003327050806304299, 0.003342190407532641, 0.003357456102345890, - 0.003372850809960137, 0.003388377524629727, 0.003404039318688046, 0.003419839345721265, - 0.003435780843885239, 0.003451867139373843, 0.003468101650046629, 0.003484487889225119, - 0.003501029469670069, 0.003517730107746697, 0.003534593627793237, 0.003551623966702611, - 0.003568825178730639, 0.003586201440546166, 0.003603757056536316, 0.003621496464384588, - 0.003639424240937217, 0.003657545108379068, 0.003675863940735269, 0.003694385770723563, - 0.003713115796977806, 0.003732059391668707, 0.003751222108547281, 0.003770609691440940, - 0.003790228083232539, 0.003810083435355216, 0.003830182117840641, 0.003850530729957835, - 0.003871136111486317, 0.003892005354668437, 0.003913145816891062, 0.003934565134149914, - 0.003956271235355358, 0.003978272357543333, 0.004000577062061084, 0.004023194251800533, - 0.004046133189565926, 0.004069403517661885, 0.004093015278800460, 0.004116978938436600, - 0.004141305408647655, 0.004166006073685835, 0.004191092817346642, 0.004216578052307351, - 0.004242474751606884, 0.004268796482457593, 0.004295557442594244, 0.004322772499391836, - 0.004350457232007221, 0.004378627976825644, 0.004407301876525049, 0.004436496933105327, - 0.004466232065271192, 0.004496527170598785, 0.004527403192966406, 0.004558882195791591, - 0.004590987441673855, 0.004623743479123199, 0.004657176237135574, 0.004691313128472929, - 0.004726183162616859, 0.004761817069491636, 0.004798247435199299, 0.004835508851176451, - 0.004873638078381815, 0.004912674228345848, 0.004952658963181422, 0.004993636716962402, - 0.005035654941235035, 0.005078764377854039, 0.005123019362831771, 0.005168478165478940, - 0.005215203367812893, 0.005263262290042703, 0.005312727468930079, 0.005363677197016692, - 0.005416196132139284, 0.005470375988385734, 0.005526316321746716, 0.005584125426278286, - 0.005643921359735682, 0.005705833121505521, 0.005770002010457520, 0.005836583196307310, - 0.005905747545561058, 0.005977683752542928, 0.006052600837980204, 0.006130731092920838, - 0.006212333565464245, 0.006297698213369562, 0.006387150879090475, 0.006481059288027780, - 0.006579840329791975, 0.006683968961788356, 0.006793989182803495, 0.006910527673723577, - 0.007034310911336661, 0.007166186857056056, 0.007307152748134871, 0.007458391141830445, - 0.007621317291194862, 0.007797642342679434, 0.007989459040836144, 0.008199360125510702, - 0.008430605346682607, 0.008687362737884952, 0.008975066840784529, 0.009300967772353674, - 0.009675004947253041, 0.010111261142904171, 0.010630518154258861, 0.011265064987797335, - 0.012068570920629962, 0.013138877484087819, 0.014680138359337902, 0.017222609470315398, - 0.022898298717268672]; pub static ZIG_EXP_R: f64 = 7.697117470131050077; pub static ZIG_EXP_X: [f64, .. 257] = [8.697117470131052741, 7.697117470131050077, 6.941033629377212577, 6.478378493832569696, @@ -344,69 +278,3 @@ pub static ZIG_EXP_F: [f64, .. 257] = 0.775956852040116218, 0.791527636972496285, 0.808421651523009044, 0.826993296643051101, 0.847785500623990496, 0.871704332381204705, 0.900469929925747703, 0.938143680862176477, 1.000000000000000000]; -pub static ZIG_EXP_F_DIFF: [f64, .. 257] = - [0.000000000000000000, 0.000287067661533533, 0.000513134928485678, 0.000569030497974398, - 0.000609667963417335, 0.000642831049855169, 0.000671465984262828, 0.000697030342996893, - 0.000720360862708599, 0.000741986223663093, 0.000762263730113694, 0.000781447246315807, - 0.000799724254382053, 0.000817237547791934, 0.000834098656693235, 0.000850396538527769, - 0.000866203416804620, 0.000881578828420777, 0.000896572504999613, 0.000911226471926952, - 0.000925576608509206, 0.000939653828282008, 0.000953484986066785, 0.000967093584871414, - 0.000980500333784669, 0.000993723593313716, 0.001006779734568374, 0.001019683431705467, - 0.001032447902101660, 0.001045085105172934, 0.001057605908173612, 0.001070020225402434, - 0.001082337135821582, 0.001094564983022843, 0.001106711460658764, 0.001118783685829211, - 0.001130788262427001, 0.001142731336065933, 0.001154618641914802, 0.001166455546523074, - 0.001178247084534012, 0.001189997991027938, 0.001201712730115490, 0.001213395520299268, - 0.001225050357040701, 0.001236681032901414, 0.001248291155571943, 0.001259884164055092, - 0.001271463343231895, 0.001283031837006378, 0.001294592660197942, 0.001306148709326875, - 0.001317702772419903, 0.001329257537945404, 0.001340815602974395, 0.001352379480650950, - 0.001363951607045839, 0.001375534347457789, 0.001387130002219621, 0.001398740812059381, - 0.001410368963061376, 0.001422016591266340, 0.001433685786946429, 0.001445378598586011, - 0.001457097036596827, 0.001468843076792140, 0.001480618663643060, 0.001492425713336909, - 0.001504266116655995, 0.001516141741693663, 0.001528054436422108, 0.001540006031125918, - 0.001551998340713470, 0.001564033166917514, 0.001576112300394977, 0.001588237522735750, - 0.001600410608388780, 0.001612633326513305, 0.001624907442762655, 0.001637234721007311, - 0.001649616925003372, 0.001662055820012304, 0.001674553174376953, 0.001687110761059388, - 0.001699730359144919, 0.001712413755316500, 0.001725162745304071, 0.001737979135312442, - 0.001750864743431488, 0.001763821401032123, 0.001776850954151601, 0.001789955264870927, - 0.001803136212688003, 0.001816395695889220, 0.001829735632922019, 0.001843157963772116, - 0.001856664651347151, 0.001870257682870316, 0.001883939071285826, 0.001897710856679738, - 0.001911575107717528, 0.001925533923102574, 0.001939589433056721, 0.001953743800826108, - 0.001967999224215228, 0.001982357937151347, 0.001996822211282223, 0.002011394357609747, - 0.002026076728162574, 0.002040871717710169, 0.002055781765521847, 0.002070809357173103, - 0.002085957026402963, 0.002101227357025226, 0.002116622984897121, 0.002132146599948981, - 0.002147800948277823, 0.002163588834309782, 0.002179513123034188, 0.002195576742314159, - 0.002211782685277469, 0.002228134012792427, 0.002244633856033434, 0.002261285419141418, - 0.002278091981983449, 0.002295056903017983, 0.002312183622271174, 0.002329475664429648, - 0.002346936642057179, 0.002364570258941101, 0.002382380313575932, 0.002400370702791893, - 0.002418545425535629, 0.002436908586812392, 0.002455464401797752, 0.002474217200128692, - 0.002493171430384328, 0.002512331664766249, 0.002531702603989994, 0.002551289082400404, - 0.002571096073321844, 0.002591128694658967, 0.002611392214760672, 0.002631892058563845, - 0.002652633814032662, 0.002673623238910738, 0.002694866267805934, 0.002716369019626269, - 0.002738137805389534, 0.002760179136428037, 0.002782499733014893, 0.002805106533435520, - 0.002828006703534697, 0.002851207646767162, 0.002874717014785921, 0.002898542718600849, - 0.002922692940346749, 0.002947176145699226, 0.002972001096982591, 0.002997176867015228, - 0.003022712853742948, 0.003048618795714386, 0.003074904788455568, 0.003101581301807876, - 0.003128659198296080, 0.003156149752600867, 0.003184064672214937, 0.003212416119368622, - 0.003241216734320596, 0.003270479660111680, 0.003300218568896729, 0.003330447689969929, - 0.003361181839619420, 0.003392436452949343, 0.003424227617828290, 0.003456572111131984, - 0.003489487437467131, 0.003522991870580083, 0.003557104497672658, 0.003591845266868621, - 0.003627235038102472, 0.003663295637722386, 0.003700049917134574, 0.003737521815846301, - 0.003775736429304177, 0.003814720081962375, 0.003854500406067995, 0.003895106426696382, - 0.003936568653631844, 0.003978919180756157, 0.004022191793678687, 0.004066422086428989, - 0.004111647588127876, 0.004157907900659452, 0.004205244848493050, 0.004253702641940915, - 0.004303328055299205, 0.004354170621502118, 0.004406282845128784, 0.004459720435841752, - 0.004514542564613699, 0.004570812145417769, 0.004628596145424491, 0.004687965927177740, - 0.004748997626717266, 0.004811772572194672, 0.004876377748206484, 0.004942906311860507, - 0.005011458167522187, 0.005082140608288488, 0.005155069033533799, 0.005230367753417398, - 0.005308170893076836, 0.005388623411430704, 0.005471882252147620, 0.005558117647517014, - 0.005647514599798176, 0.005740274569295156, 0.005836617404105682, 0.005936783553485037, - 0.006041036615386131, 0.006149666279423593, 0.006262991739818591, 0.006381365669577810, - 0.006505178868201678, 0.006634865721946159, 0.006770910649812723, 0.006913855752425535, - 0.007064309938019209, 0.007222959874423007, 0.007390583214465396, 0.007568064673498798, - 0.007756415714389786, 0.007956798835585532, 0.008170557788458321, 0.008399255510700199, - 0.008644722212900025, 0.008909116987305010, 0.009195007664428712, 0.009505475652925033, - 0.009844255532840629, 0.010215923852312625, 0.010626158965710175, 0.011082105722287849, - 0.011592898788496009, 0.012170432837851575, 0.012830529553771619, 0.013594766864701180, - 0.014493463190219380, 0.015570784932380066, 0.016894014550512759, 0.018571645120042057, - 0.020792203980939394, 0.023918831757214210, 0.028765597544542998, 0.037673750936428774, - 0.061856319137823523];