Line data Source code
1 : #include "fd_rng.h"
2 : #include <math.h> /* FIXME: ELIMINATE DEPENDENCE AS THIS IS NOT GUARANTEED BIT LEVEL IDENTICAL BETWEEN PLATFORMS */
3 :
4 : float
5 3000 : fd_rng_float_robust( fd_rng_t * rng ) {
6 :
7 : /* Generate a float as though a continuum rand was generated in [1,2]
8 : and then round to nearest even rounded to the closest exactly
9 : representable float point value in [1,2]. (This can generate all
10 : such representable values in that range.) */
11 :
12 3000 : ulong u = (ulong)fd_rng_uint( rng ); /* 32-bit uniform rand */
13 3000 : ulong m = (u >> 9) + (u & 1UL); /* [0,2^23] trapezoidal rand */
14 3000 : float f = 1.f + ((float)m)*FLT_EPSILON; /* Exact */
15 :
16 : /* Now apply an exponent to yield the correct final result */
17 :
18 3000 : return ldexpf( f, -(int)fd_rng_coin_tosses( rng ) );
19 3000 : }
20 :
21 : float
22 58838621 : fd_rng_float_exp( fd_rng_t * rng ) {
23 58838621 : ulong u = 1UL + (fd_rng_ulong( rng ) >> 1); /* In (0,2^63] uniform */
24 58838621 : float f = ((float)u) * (1.f/(float)(1UL<<63)); /* In uniform (0,1] with reasonably low quant, exact */
25 58838621 : return -logf( f ); /* Convert to exp */
26 58838621 : }
27 :
28 : #if 0
29 : /* USE THIS WITH LOGF AT END REPLACED WITH FAST APPROX REJECTION IF NO
30 : LIBM. USE TRANSFORMATION ON CONTINUUM RAND IF WANTING NO QUANT
31 : ARTIFACTS. NEARLY EQUIVALENT CAN USE:
32 :
33 : F(x) = 1 - exp(-x) = 1 - exp2(-x*log2(e))
34 : -> y = x*log2(e), x = y*(1/log2(e))
35 : F(y) = 1 - exp2(-y)
36 : -> Y = -log2( U )
37 :
38 : where U is a uniform rand in [0,1]. But U can be generated by:
39 :
40 : U = (2^-E)( 1 + M )
41 :
42 : where E is IID integer rand with probability:
43 :
44 : P(E) = 2^-E for E>0
45 : = 0 otherwise
46 :
47 : and M is an IID uniform rand in [0,1]. Thus:
48 :
49 : -> Y = E - log2( 1 + M )
50 :
51 : yielding a variant with minimal tail artifacts. In single precision:
52 :
53 : e = fd_rng_coin_tosses( rng );
54 : m = fd_rng_uint( rng );
55 : y = e - lg1p( m / 2^32 );
56 : return y*(1/log2(e));
57 :
58 : and for double precision:
59 :
60 : e = fd_rng_coin_tosses( rng )
61 : m = fd_rng_ulong( rng );
62 : y = e - lg1p( m / 2^64 );
63 : return y*(1/log2(e));
64 :
65 : Eliminating floating point rep quant artifacts near the origin is
66 : much trickier due to the cancellations implied above (but also
67 : practically much less important ... they are no worse than the
68 : quantization artifacts of the typical generator). */
69 :
70 : float
71 : fd_rng_float_exp( fd_rng_t * rng ) {
72 :
73 : for(;;) {
74 :
75 : /* Generate the 7-bit ziggurat level and a 24-bit trapezoidal rand
76 : (e.g. "uniform" rand on [0,2^24] with the endpoints 0 and 2^24
77 : half as likely as the other points). */
78 :
79 : ulong u = (ulong)fd_rng_uint( rng );
80 :
81 : ulong l = (u >> 1) & (zig_level_cnt-1UL); /* uniform level (level_cnt must be a power-of-2 <= 128) */
82 : ulong m = (u >> 8) + (u & 1UL); /* ~2^24 bit trapezoidal rand */
83 :
84 : /* Compute a uniform rand as wide as the current ziggurat level
85 : (will be closed at both endpoints due to trapezoid rounding
86 : above). If willing to have an additional table, the
87 : multiplication by 2^-24 can be eliminated (technically though this
88 : is exact and can be done just be playing with with the exponent). */
89 :
90 : float x = zig_x[l+1UL]*((1.f/(float)(1<<24))*(float)m);
91 :
92 : /* If this is smaller than the level above this, we can immediately
93 : accept this rand. This occurs the vast majority of the time. */
94 :
95 : if( FD_LIKELY( x<=zig_x[l] ) ) break; /* Quick accept */
96 :
97 : /* We can't do a quick acceptance */
98 :
99 : if( FD_LIKELY( l<(zig_level_cnt-1UL) ) ) {
100 :
101 : /* We are picking a point uniformly in the l-th ziggurat strip. */
102 :
103 : float t = zig_y[l];
104 : float y = t + (zig_y[l+1UL]-t)*fd_rng_float_c0( rng );
105 : if( y < expf(-x) ) break; /* Fundamentally ~50/50 unpredictable */
106 :
107 : } else {
108 :
109 : /* We are picking a point uniformly in the ziggurat exponential
110 : tail. Note that x - R is an independent uniform rand between 0
111 : and 1. In finite precision arithmetic, (x-R)==0 will never
112 : occur due to the above acceptance check. Thus, it is safe to
113 : compute ln (x-R) to generate an exponentially distributed rand
114 : for the tail here. */
115 :
116 : x = zig_r - logf( x - zig_r );
117 : break;
118 :
119 : }
120 :
121 : /* Rejected */
122 : }
123 :
124 : return x;
125 : }
126 : #endif
127 :
128 : float
129 202551 : fd_rng_float_norm( fd_rng_t * rng ) {
130 :
131 : /* BEGIN AUTOGENERATED CODE *****************************************/
132 :
133 202551 : static float const zig_x[65] = {
134 202551 : 0.000000000000000000000000000000e+00f /* l= 0 */, 3.455520224239806172095417630130e-01f /* l= 1 */,
135 202551 : 4.621556810233836948395482607799e-01f /* l= 2 */, 5.450460069613645023634852793126e-01f /* l= 3 */,
136 202551 : 6.119536252363746848246563170282e-01f /* l= 4 */, 6.693178808991654190761917686547e-01f /* l= 5 */,
137 202551 : 7.202743599269989995133832427765e-01f /* l= 6 */, 7.666087419009196196129432565591e-01f /* l= 7 */,
138 202551 : 8.094445383821111665470331153482e-01f /* l= 8 */, 8.495396579691822718460002261676e-01f /* l= 9 */,
139 202551 : 8.874325556154150376969927394022e-01f /* l=10 */, 9.235214816642028361095943800319e-01f /* l=11 */,
140 202551 : 9.581106750231991023205452284728e-01f /* l=12 */, 9.914388611066755936521502357017e-01f /* l=13 */,
141 202551 : 1.023697658831656335945772817730e+00f /* l=14 */, 1.055043930218579828852511204307e+00f /* l=15 */,
142 202551 : 1.085608336108546385250818444579e+00f /* l=16 */, 1.115501429269408685757483667977e+00f /* l=17 */,
143 202551 : 1.144818099662757484029536325654e+00f /* l=18 */, 1.173640887904643574190903521082e+00f /* l=19 */,
144 202551 : 1.202042503654549259546786832420e+00f /* l=20 */, 1.230087774542608308986184340039e+00f /* l=21 */,
145 202551 : 1.257835180410570578152240628356e+00f /* l=22 */, 1.285338081362219465062293743962e+00f /* l=23 */,
146 202551 : 1.312645717221077709701594626868e+00f /* l=24 */, 1.339804034974004910067517382100e+00f /* l=25 */,
147 202551 : 1.366856386251642895349189821275e+00f /* l=26 */, 1.393844126728532207214715510357e+00f /* l=27 */,
148 202551 : 1.420807142148611732640484106582e+00f /* l=28 */, 1.447784320603175013198252174540e+00f /* l=29 */,
149 202551 : 1.474813987119431635871430463780e+00f /* l=30 */, 1.501934314168690078094073325765e+00f /* l=31 */,
150 202551 : 1.529183720117960190352200677832e+00f /* l=32 */, 1.556601266765526946100374472426e+00f /* l=33 */,
151 202551 : 1.584227066827416976904988055175e+00f /* l=34 */, 1.612102712541159403396816285348e+00f /* l=35 */,
152 202551 : 1.640271737439234119525742483514e+00f /* l=36 */, 1.668780124881012518812430089898e+00f /* l=37 */,
153 202551 : 1.697676879240338253859027295434e+00f /* l=38 */, 1.727014678920001637717554499041e+00f /* l=39 */,
154 202551 : 1.756850634895277701885175913876e+00f /* l=40 */, 1.787247184704316305928380181900e+00f /* l=41 */,
155 202551 : 1.818273160330322024778176848159e+00f /* l=42 */, 1.850005080182403822074664601072e+00f /* l=43 */,
156 202551 : 1.882528731753153053690014173682e+00f /* l=44 */, 1.915941134586538274016051519588e+00f /* l=45 */,
157 202551 : 1.950353006115250048047894682046e+00f /* l=46 */, 1.985891900706324235165341207665e+00f /* l=47 */,
158 202551 : 2.022706262852767138214066244828e+00f /* l=48 */, 2.060970741899825969985216023161e+00f /* l=49 */,
159 202551 : 2.100893279890740158706921580922e+00f /* l=50 */, 2.142724743933560004960706124599e+00f /* l=51 */,
160 202551 : 2.186772297656567175727290730514e+00f /* l=52 */, 2.233418418575520170057252533624e+00f /* l=53 */,
161 202551 : 2.283148713210882645120378131587e+00f /* l=54 */, 2.336593955757456149034331782666e+00f /* l=55 */,
162 202551 : 2.394596149760118465339014948157e+00f /* l=56 */, 2.458317361057721973189790776182e+00f /* l=57 */,
163 202551 : 2.529429815978553923233249078883e+00f /* l=58 */, 2.610473649022010540245164467166e+00f /* l=59 */,
164 202551 : 2.705599986069530120576243081842e+00f /* l=60 */, 2.822342497323976174845167053107e+00f /* l=61 */,
165 202551 : 2.976823798790561124028714035106e+00f /* l=62 */, 3.215929245508522913137711141118e+00f /* l=63 */,
166 202551 : 3.526881360327788386254538322007e+00f /* l=64 */
167 202551 : };
168 :
169 202551 : static float const zig_y[65] = {
170 202551 : 1.000000000000000000000000000000e+00f /* l= 0 */, 9.420441848916480126344477619149e-01f /* l= 1 */,
171 202551 : 8.987108451876446194939510037081e-01f /* l= 2 */, 8.619676182560859538047716432718e-01f /* l= 3 */,
172 202551 : 8.292416921465940136417270556191e-01f /* l= 4 */, 7.993205594629470858670725053052e-01f /* l= 5 */,
173 202551 : 7.715162251201959925242870874662e-01f /* l= 6 */, 7.453924046791949711084253327176e-01f /* l= 7 */,
174 202551 : 7.206510565419278137484770940802e-01f /* l= 8 */, 6.970774082324479442845238663651e-01f /* l= 9 */,
175 202551 : 6.745103421549214010562514620695e-01f /* l=10 */, 6.528251409771071645448889397834e-01f /* l=11 */,
176 202551 : 6.319228072029464656281239065549e-01f /* l=12 */, 6.117231258029589507463179287594e-01f /* l=13 */,
177 202551 : 5.921599774953054971587604327077e-01f /* l=14 */, 5.731780673128809060960613119828e-01f /* l=15 */,
178 202551 : 5.547305771308242055075976573164e-01f /* l=16 */, 5.367774409030761823035174384877e-01f /* l=17 */,
179 202551 : 5.192840512302357241190831071975e-01f /* l=18 */, 5.022202719018961842071223367068e-01f /* l=19 */,
180 202551 : 4.855596720803138878681283474581e-01f /* l=20 */, 4.692789240423391189252835808965e-01f /* l=21 */,
181 202551 : 4.533573236341009112014902027177e-01f /* l=22 */, 4.377764041761660086522593704483e-01f /* l=23 */,
182 202551 : 4.225196225029520580892660602812e-01f /* l=24 */, 4.075721013736326630086392874830e-01f /* l=25 */,
183 202551 : 3.929204164392402009802710699526e-01f /* l=26 */, 3.785524188002748156032031823237e-01f /* l=27 */,
184 202551 : 3.644570862756677649099863736115e-01f /* l=28 */, 3.506243980520670366787892163751e-01f /* l=29 */,
185 202551 : 3.370452285453843354473546511940e-01f /* l=30 */, 3.237112571905811243399165438861e-01f /* l=31 */,
186 202551 : 3.106148915554719103920833234156e-01f /* l=32 */, 2.977492017031598617858058342112e-01f /* l=33 */,
187 202551 : 2.851078641441273595351531267017e-01f /* l=34 */, 2.726851140512668792348013879767e-01f /* l=35 */,
188 202551 : 2.604757046803618125588395543213e-01f /* l=36 */, 2.484748731607814430951243836465e-01f /* l=37 */,
189 202551 : 2.366783120090016554051349020882e-01f /* l=38 */, 2.250821458811925040717320800621e-01f /* l=39 */,
190 202551 : 2.136829132292288532124100927656e-01f /* l=40 */, 2.024775526650522966412917846846e-01f /* l=41 */,
191 202551 : 1.914633939793006458976280803608e-01f /* l=42 */, 1.806381539102589138648434843870e-01f /* l=43 */,
192 202551 : 1.699999369289590334732324358735e-01f /* l=44 */, 1.595472415092518358135931233477e-01f /* l=45 */,
193 202551 : 1.492789726065822738132905095343e-01f /* l=46 */, 1.391944614029269361199806984142e-01f /* l=47 */,
194 202551 : 1.292934938280917421830808894390e-01f /* l=48 */, 1.195763500012655360156695570628e-01f /* l=49 */,
195 202551 : 1.100438576497440733127446653439e-01f /* l=50 */, 1.006974639150281902163835620612e-01f /* l=51 */,
196 202551 : 9.153933202201729600920110732631e-02f /* l=52 */, 8.257247254089309598817534793791e-02f /* l=53 */,
197 202551 : 7.380092428122808796957937671479e-02f /* l=54 */, 6.523000887995579990570474762657e-02f /* l=55 */,
198 202551 : 5.686669921543156531844320777935e-02f /* l=56 */, 4.872017206675432561731060171484e-02f /* l=57 */,
199 202551 : 4.080267659191996105732375653419e-02f /* l=58 */, 3.313098485527892895195507111383e-02f /* l=59 */,
200 202551 : 2.572902254561226234667566450942e-02f /* l=60 */, 1.863323273948142502092796546354e-02f /* l=61 */,
201 202551 : 1.190567663419300047354134047123e-02f /* l=62 */, 5.678316641776698130187294149412e-03f /* l=63 */,
202 202551 : 0.000000000000000000000000000000e+00f /* l=64 */
203 202551 : };
204 :
205 202551 : static ulong const zig_level_cnt = 64;
206 :
207 202551 : static float const zig_r = 3.215929245508522913137711141118e+00f;
208 202551 : static float const zig_rcp_r = 3.109521148192654732252473981369e-01f;
209 202551 : static float const zig_half_r = 1.607964622754261456568855570559e+00f;
210 202551 : static float const zig_tail = 2.852700996027780249580940719056e+00f;
211 :
212 : /* END AUTOGENERATED CODE *******************************************/
213 :
214 202551 : ulong s;
215 202551 : float x;
216 :
217 207225 : for(;;) {
218 :
219 : /* Generate the sign bit, up to 6-bit ziggurat level and a 24-bit
220 : trapezoidal rand (e.g. "uniform" rand on [0,2^24] with the
221 : endpoints 0 and 2^24 half as likely as the other points). */
222 :
223 207225 : ulong u = (ulong)fd_rng_uint( rng ); /* 32-bit rand */
224 :
225 207225 : /**/ s = (u >> 1) & 1UL; /* random sign */
226 207225 : ulong l = (u >> 2) & (zig_level_cnt-1UL); /* uniform level (zig_level_cnt must be power-of-2 <= 64) */
227 207225 : ulong m = (u >> 8) + (u & 1UL); /* 24-bit trapezoidal rand */
228 :
229 : /* Compute a uniform rand as wide as the current ziggurat level */
230 :
231 207225 : x = zig_x[l+1UL]*((1.f/16777216.f)*(float)m); /* Guaranteed in [0,x[l+1]] */
232 :
233 : /* If this is <= level above this, we can immediately accept this
234 : rand. This occurs the vast majority of the time. */
235 :
236 207225 : if( FD_LIKELY( x<=zig_x[l] ) ) break; /* Quick accept */
237 :
238 : /* We can't do a quick acceptance ... the y rand might matter */
239 :
240 10386 : float y = fd_rng_float_c( rng ); /* Guaranteed in [0,1] */
241 :
242 10386 : if( FD_LIKELY( l<(zig_level_cnt-1UL) ) ) {
243 :
244 : /* We are picking a point uniformly in the l-th ziggurat strip. */
245 :
246 10092 : y = zig_y[l]*(1.f-y) + zig_y[l+1UL]*y; /* Guaranteed in [y0,y1] */
247 :
248 10092 : } else {
249 :
250 : /* We are picking a point uniformly distributed in the ziggurat
251 : exponential tail.
252 :
253 : The simplest way to do this is to transform a new uniform
254 : rand into an exponential distribution that joins smoothly to
255 : the base of the Ziggurat.
256 :
257 : We theoretically can do slightly better by noting that x-r
258 : is an independent uniform rand in ( r==x[l], x[l+1] ]. Thus,
259 : instead of computing a new uniform rand for the tail x and
260 : transforming it to the desired exponential rand, we could get
261 : an exponential rand via:
262 :
263 : x' = r - (1/r) ln ( (x-r) / (x[l+1]-r) )
264 :
265 : which simplifies to:
266 :
267 : x = (r + (1/r) ln( x[l+1]-r) ) - (1/r) ln ( x-r )
268 :
269 : The first term is a compile time constant. The second term
270 : doesn't require a new rand but is otherwise the same
271 : transformation cost. x-r is guaranteed greater than 0 given
272 : the quick accept test above so the logf will not blow up. */
273 :
274 294 : x = zig_tail - zig_rcp_r*logf( x - zig_r );
275 294 : y *= expf( zig_r*(zig_half_r - x) );
276 :
277 294 : }
278 :
279 : /* FIXME: COULD GET RID OF 0.5 MULTIPLICATION BY USING PDF
280 : expf(-x^2), ADJUST TAIL COMP AND SCALING TO UNIT NORM AS PART OF
281 : THE SIGN INSERTION. */
282 :
283 10386 : if( y < expf((-0.5f)*(x*x)) ) break; /* Slow accept, approx 50/50 unpredictable */
284 :
285 : /* Rejected, try again */
286 10386 : }
287 :
288 202551 : return fd_float_if( (int)s, -x, x ); /* FIXME: LOAD AND MUL BY +/-1? LOAD AND COPYSIGNF BY +/-1? */
289 202551 : }
290 :
291 : #if FD_HAS_DOUBLE /* See float variants for details how these work */
292 :
293 : double
294 3000 : fd_rng_double_robust( fd_rng_t * rng ) {
295 3000 : ulong u = fd_rng_ulong( rng );
296 3000 : ulong m = (u >> 12) + (u & 1UL);
297 3000 : double d = 1. + ((double)m)*DBL_EPSILON;
298 3000 : return ldexp( d, -(int)fd_rng_coin_tosses( rng ) );
299 3000 : }
300 :
301 : double
302 219021 : fd_rng_double_exp( fd_rng_t * rng ) {
303 219021 : ulong u = 1UL + (fd_rng_ulong( rng ) >> 1);
304 219021 : double d = ((double)u) * (1./(double)(1UL<<63));
305 219021 : return -log( d );
306 219021 : }
307 :
308 : double
309 1112136 : fd_rng_double_norm( fd_rng_t * rng ) {
310 :
311 : /* BEGIN AUTOGENERATED CODE *****************************************/
312 :
313 1112136 : static double const zig_x[65] = {
314 1112136 : 0.000000000000000000000000000000e+00, 3.455520224239806172095417630130e-01,
315 1112136 : 4.621556810233836948395482607799e-01, 5.450460069613645023634852793126e-01,
316 1112136 : 6.119536252363746848246563170282e-01, 6.693178808991654190761917686547e-01,
317 1112136 : 7.202743599269989995133832427765e-01, 7.666087419009196196129432565591e-01,
318 1112136 : 8.094445383821111665470331153482e-01, 8.495396579691822718460002261676e-01,
319 1112136 : 8.874325556154150376969927394022e-01, 9.235214816642028361095943800319e-01,
320 1112136 : 9.581106750231991023205452284728e-01, 9.914388611066755936521502357017e-01,
321 1112136 : 1.023697658831656335945772817730e+00, 1.055043930218579828852511204307e+00,
322 1112136 : 1.085608336108546385250818444579e+00, 1.115501429269408685757483667977e+00,
323 1112136 : 1.144818099662757484029536325654e+00, 1.173640887904643574190903521082e+00,
324 1112136 : 1.202042503654549259546786832420e+00, 1.230087774542608308986184340039e+00,
325 1112136 : 1.257835180410570578152240628356e+00, 1.285338081362219465062293743962e+00,
326 1112136 : 1.312645717221077709701594626868e+00, 1.339804034974004910067517382100e+00,
327 1112136 : 1.366856386251642895349189821275e+00, 1.393844126728532207214715510357e+00,
328 1112136 : 1.420807142148611732640484106582e+00, 1.447784320603175013198252174540e+00,
329 1112136 : 1.474813987119431635871430463780e+00, 1.501934314168690078094073325765e+00,
330 1112136 : 1.529183720117960190352200677832e+00, 1.556601266765526946100374472426e+00,
331 1112136 : 1.584227066827416976904988055175e+00, 1.612102712541159403396816285348e+00,
332 1112136 : 1.640271737439234119525742483514e+00, 1.668780124881012518812430089898e+00,
333 1112136 : 1.697676879240338253859027295434e+00, 1.727014678920001637717554499041e+00,
334 1112136 : 1.756850634895277701885175913876e+00, 1.787247184704316305928380181900e+00,
335 1112136 : 1.818273160330322024778176848159e+00, 1.850005080182403822074664601072e+00,
336 1112136 : 1.882528731753153053690014173682e+00, 1.915941134586538274016051519588e+00,
337 1112136 : 1.950353006115250048047894682046e+00, 1.985891900706324235165341207665e+00,
338 1112136 : 2.022706262852767138214066244828e+00, 2.060970741899825969985216023161e+00,
339 1112136 : 2.100893279890740158706921580922e+00, 2.142724743933560004960706124599e+00,
340 1112136 : 2.186772297656567175727290730514e+00, 2.233418418575520170057252533624e+00,
341 1112136 : 2.283148713210882645120378131587e+00, 2.336593955757456149034331782666e+00,
342 1112136 : 2.394596149760118465339014948157e+00, 2.458317361057721973189790776182e+00,
343 1112136 : 2.529429815978553923233249078883e+00, 2.610473649022010540245164467166e+00,
344 1112136 : 2.705599986069530120576243081842e+00, 2.822342497323976174845167053107e+00,
345 1112136 : 2.976823798790561124028714035106e+00, 3.215929245508522913137711141118e+00,
346 1112136 : 3.526881360327788386254538322007e+00
347 1112136 : };
348 :
349 1112136 : static double const zig_y[65] = {
350 1112136 : 1.000000000000000000000000000000e+00, 9.420441848916480126344477619149e-01,
351 1112136 : 8.987108451876446194939510037081e-01, 8.619676182560859538047716432718e-01,
352 1112136 : 8.292416921465940136417270556191e-01, 7.993205594629470858670725053052e-01,
353 1112136 : 7.715162251201959925242870874662e-01, 7.453924046791949711084253327176e-01,
354 1112136 : 7.206510565419278137484770940802e-01, 6.970774082324479442845238663651e-01,
355 1112136 : 6.745103421549214010562514620695e-01, 6.528251409771071645448889397834e-01,
356 1112136 : 6.319228072029464656281239065549e-01, 6.117231258029589507463179287594e-01,
357 1112136 : 5.921599774953054971587604327077e-01, 5.731780673128809060960613119828e-01,
358 1112136 : 5.547305771308242055075976573164e-01, 5.367774409030761823035174384877e-01,
359 1112136 : 5.192840512302357241190831071975e-01, 5.022202719018961842071223367068e-01,
360 1112136 : 4.855596720803138878681283474581e-01, 4.692789240423391189252835808965e-01,
361 1112136 : 4.533573236341009112014902027177e-01, 4.377764041761660086522593704483e-01,
362 1112136 : 4.225196225029520580892660602812e-01, 4.075721013736326630086392874830e-01,
363 1112136 : 3.929204164392402009802710699526e-01, 3.785524188002748156032031823237e-01,
364 1112136 : 3.644570862756677649099863736115e-01, 3.506243980520670366787892163751e-01,
365 1112136 : 3.370452285453843354473546511940e-01, 3.237112571905811243399165438861e-01,
366 1112136 : 3.106148915554719103920833234156e-01, 2.977492017031598617858058342112e-01,
367 1112136 : 2.851078641441273595351531267017e-01, 2.726851140512668792348013879767e-01,
368 1112136 : 2.604757046803618125588395543213e-01, 2.484748731607814430951243836465e-01,
369 1112136 : 2.366783120090016554051349020882e-01, 2.250821458811925040717320800621e-01,
370 1112136 : 2.136829132292288532124100927656e-01, 2.024775526650522966412917846846e-01,
371 1112136 : 1.914633939793006458976280803608e-01, 1.806381539102589138648434843870e-01,
372 1112136 : 1.699999369289590334732324358735e-01, 1.595472415092518358135931233477e-01,
373 1112136 : 1.492789726065822738132905095343e-01, 1.391944614029269361199806984142e-01,
374 1112136 : 1.292934938280917421830808894390e-01, 1.195763500012655360156695570628e-01,
375 1112136 : 1.100438576497440733127446653439e-01, 1.006974639150281902163835620612e-01,
376 1112136 : 9.153933202201729600920110732631e-02, 8.257247254089309598817534793791e-02,
377 1112136 : 7.380092428122808796957937671479e-02, 6.523000887995579990570474762657e-02,
378 1112136 : 5.686669921543156531844320777935e-02, 4.872017206675432561731060171484e-02,
379 1112136 : 4.080267659191996105732375653419e-02, 3.313098485527892895195507111383e-02,
380 1112136 : 2.572902254561226234667566450942e-02, 1.863323273948142502092796546354e-02,
381 1112136 : 1.190567663419300047354134047123e-02, 5.678316641776698130187294149412e-03,
382 1112136 : 0.000000000000000000000000000000e+00
383 1112136 : };
384 :
385 1112136 : static ulong const zig_level_cnt = 64;
386 :
387 1112136 : static double const zig_r = 3.215929245508522913137711141118e+00;
388 1112136 : static double const zig_rcp_r = 3.109521148192654732252473981369e-01;
389 1112136 : static double const zig_half_r = 1.607964622754261456568855570559e+00;
390 1112136 : static double const zig_tail = 2.852700996027780249580940719056e+00;
391 :
392 : /* END AUTOGENERATED CODE *******************************************/
393 :
394 1112136 : ulong s;
395 1112136 : double x;
396 :
397 1137048 : for(;;) {
398 :
399 1137048 : ulong u = fd_rng_ulong( rng ); /* 64-bit rand */
400 :
401 1137048 : /**/ s = (u >> 1) & 1UL;
402 1137048 : ulong l = (u >> 2) & (zig_level_cnt-1UL); /* uniform level (zig_level_cnt must be power-of-2 <= 512) */
403 1137048 : ulong m = (u >> 11) + (u & 1UL); /* 53-bit trapezoidal rand */
404 :
405 1137048 : x = zig_x[l+1UL]*((1./9007199254740992.)*(double)m);
406 :
407 1137048 : if( FD_LIKELY( x<=zig_x[l] ) ) break;
408 :
409 57147 : double y = fd_rng_double_c( rng );
410 :
411 57147 : if( FD_LIKELY( l<(zig_level_cnt-1UL) ) ) {
412 55632 : y = zig_y[l]*(1.-y) + zig_y[l+1UL]*y;
413 55632 : } else {
414 1515 : x = zig_tail - zig_rcp_r*log( x - zig_r );
415 1515 : y *= exp( zig_r*(zig_half_r - x) );
416 1515 : }
417 :
418 57147 : if( y < exp((-0.5)*(x*x)) ) break;
419 57147 : }
420 :
421 1112136 : return fd_double_if( (int)s, -x, x );
422 1112136 : }
423 :
424 : #endif
425 :
|
