From 80f98830da730e3f99a8ef24b8e1e5e5abbb849d Mon Sep 17 00:00:00 2001 From: NunoSempere Date: Sat, 3 Jun 2023 03:52:04 -0600 Subject: [PATCH] start adding xorshift prng. --- C/out/samples | Bin 22456 -> 22496 bytes C/samples.c | 54 ++++++++++++++++++++++++++++++++------------------ 2 files changed, 35 insertions(+), 19 deletions(-) diff --git a/C/out/samples b/C/out/samples index 8c046ba8efbca3988f22ddf50caca31903b40c08..94914fa14d4fa901388ad12c0b85505dab22eb36 100755 GIT binary patch delta 1994 zcmZ`(e@v8R9DkpKlRE_OPCz__^A2utfW#dfN>kEPggZc`#1Q8c1DQ5gR5r=gl*H9d zX(!*>ve`OaTh!b`8;E#-0}Ya+32-!_CN~p%|-GyOlzDmN{>N_4;r6M`FE- z{#(VynB}DLJj;xiRQ5JJie}O+YZet-(o_FL6_bW_k=P;fSZu7G+AaBc=Wi!67@O^B zlLlj>U6;tVgN|6LYF~qelC4{@eUQyA2P{4f`w>+49o8@kPqWKeo>8qu6&91IVpj5w z#mJBheXI6YqgqWrHOwv5xnX_ni_QQ=EIGBqw-PA_u|snYQg+B*|CjP()mo(7U{q5s zdysNYBIQ|{Zk;#n?gZztE~>T`KJ9>`+Ve3Uj96l~=8gF0e zMbN7oi>C01b+PxukK2z-?|1g@691Q6tl;mQTNyBlQ2|51HgO)~Jen*?qn9|>;T-vc z^KRfMa6d2#3<8~h!44RVaqiSeQrkG^>w(s*oWB7K0#5;(u5tbgFaV5eq)sV>jg%xw zU2)Evfg#{-pbqU>oqX;4zy!;tWhfVOHEs6>tLB0hI7O0>BXdYesw9R{V$X+@Fmk+Sp zLR0QK{XB#ccQ8mRDhu^BIN34rcc`awnQj|)5!?;kT9OeC*atzmtG-dy z=_}Rl0(cXmI-LrEUS+vRW8Q3YT=fdhRoG!wJ0)WJX24j~6U3kk*HA>G4OcZqlh|G|H<#%L>`IIQdZPJtT-k-@ qGEdPVy=wYDrsWvi4^4x@EVmj3|1?yZ;r delta 1850 zcmZ`(YfO_@7(VBtbu9vIxmc=DTIfarnJtJjb4;oXT4Y)$3|O4vyu{1qFf-95%Y?`x z-6W~JBWki_XxtJLqqy#e8Z!$n8#BuWiKS4znMLCpqF#K33Io7B8FRET(C;%|`?c_$KW?Mm3wSXj4wHt`6$!r`822Vx3-f z~}?laTFkum1~qPF2exKQyW#7d{C&Z!+Wr71?G^$$&ZIxaH6(G)Y+8WZ;;bojyhLpPsMFql!N2uD;x)>QbnSz%hkR^ z{w8K%_WUkHCI?1@r~vi?8|C$g5Z`JT8@wU}!#ezj5G}w#U^_4Z+y`|1CB!hWCn|(f zVyyY95OqM?H6dDoJ;2jI|8*e-fkEJ?#8$9g827Q!aszcIATV${(1^M|1y%r$0LOvD z@*4da2QsugOHK^v1x9d7TY*7*XZL~Z4ho?S?tDH#d)#K;N{3-L(vaJ1@uS#rNG1oK zAlr@UC0>%2QHfHMRq7?+*O(PlP%@wUX+y~p-a;o!iX_b#ohr%I&jtAz>fT4D(s>yd z$Ambny2^*%0JEJ|m8NAJLSHgf#t%$K$-NXPUB#oF3y4*_IU$#)Y9M*}@NLlqW6qYi{S>DJyLT%+M zHN1mHJk@-Cs9?c8&d-y*@&da2pz@3*m@VhRSUhRq0NRksj4Vt;uT)!FtZ|aX6Nb4+ zW+fBr2)V0{ac-oM>TM>!O%Wx4lt(wUdCiss$WuPZhJQ_YJgc%_QeU;%5?0x;8sHm6 z#D1a?uh~7MvMx2bc(x>@=Thr)Nz`*pkcDDnDPruA-RZxUNVo_fKhaXl?)n4@>qDC|*cYj{C_ Nd9kvZPW-i;{{Rj&hkpP7 diff --git a/C/samples.c b/C/samples.c index a38f3532..ba77f3da 100644 --- a/C/samples.c +++ b/C/samples.c @@ -1,8 +1,8 @@ -#include #include +#include #include +#include #include -#include const float PI = 3.14159265358979323846; @@ -77,16 +77,32 @@ float split_array_sum(float** meta_array, int length, int divided_into) } +// Pseudo Random number generator + +uint32_t xorshift32(uint32_t* state) +{ + /* Algorithm "xor" from p. 4 of Marsaglia, "Xorshift RNGs" */ + uint32_t x = *state; + x ^= x << 13; + x ^= x >> 17; + x ^= x << 5; + return *state = x; +} + +inline float rand_xorshift32(uint32_t* state){ + return (float) xorshift32(state) / UINT32_MAX; +} + // Distribution & sampling functions -float rand_float(float to, unsigned int* seed) +float rand_float(float to, uint32_t* seed) { return ((float)rand_r(seed) / (float)RAND_MAX) * to; // See: for why to use rand_r: // rand() is not thread-safe, as it relies on (shared) hidden state. } -float ur_normal(unsigned int* seed) +float ur_normal(uint32_t* seed) { float u1 = rand_float(1.0, seed); float u2 = rand_float(1.0, seed); @@ -94,22 +110,22 @@ float ur_normal(unsigned int* seed) return z; } -float random_uniform(float from, float to, unsigned int* seed) +float random_uniform(float from, float to, uint32_t* seed) { return ((float)rand_r(seed) / (float)RAND_MAX) * (to - from) + from; } -float random_normal(float mean, float sigma, unsigned int* seed) +float random_normal(float mean, float sigma, uint32_t* seed) { return (mean + sigma * ur_normal(seed)); } -float random_lognormal(float logmean, float logsigma, unsigned int* seed) +float random_lognormal(float logmean, float logsigma, uint32_t* seed) { return expf(random_normal(logmean, logsigma, seed)); } -float random_to(float low, float high, unsigned int* seed) +float random_to(float low, float high, uint32_t* seed) { const float NORMAL95CONFIDENCE = 1.6448536269514722; float loglow = logf(low); @@ -120,7 +136,7 @@ float random_to(float low, float high, unsigned int* seed) } // Mixture function -void mixture(float (*samplers[])(unsigned int*), float* weights, int n_dists, float** results, int n_threads) +void mixture(float (*samplers[])(uint32_t*), float* weights, int n_dists, float** results, int n_threads) { // You can see a simpler version of this function in the git history // or in C-02-better-algorithm-one-thread/ @@ -135,10 +151,10 @@ void mixture(float (*samplers[])(unsigned int*), float* weights, int n_dists, fl float p1; int sample_index, i, split_array_length; - // unsigned int* seeds[n_threads]; - unsigned int** seeds = malloc(n_threads * sizeof(unsigned int*)); - for (unsigned int i = 0; i < n_threads; i++) { - seeds[i] = malloc(sizeof(unsigned int)); + // uint32_t* seeds[n_threads]; + uint32_t** seeds = malloc(n_threads * sizeof(uint32_t*)); + for (uint32_t i = 0; i < n_threads; i++) { + seeds[i] = malloc(sizeof(uint32_t)); *seeds[i] = i; } @@ -161,7 +177,7 @@ void mixture(float (*samplers[])(unsigned int*), float* weights, int n_dists, fl // free(normalized_weights); // free(cummulative_weights); free(cumsummed_normalized_weights); - for (unsigned int i = 0; i < n_threads; i++) { + for (uint32_t i = 0; i < n_threads; i++) { free(seeds[i]); } free(seeds); @@ -170,22 +186,22 @@ void mixture(float (*samplers[])(unsigned int*), float* weights, int n_dists, fl // Functions used for the BOTEC. // Their type has to be the same, as we will be passing them around. -float sample_0(unsigned int* seed) +float sample_0(uint32_t* seed) { return 0; } -float sample_1(unsigned int* seed) +float sample_1(uint32_t* seed) { return 1; } -float sample_few(unsigned int* seed) +float sample_few(uint32_t* seed) { return random_to(1, 3, seed); } -float sample_many(unsigned int* seed) +float sample_many(uint32_t* seed) { return random_to(2, 10, seed); } @@ -210,7 +226,7 @@ int main() // Generate mixture int n_dists = 4; float weights[] = { 1 - p_c, p_c / 2, p_c / 4, p_c / 4 }; - float (*samplers[])(unsigned int*) = { sample_0, sample_1, sample_few, sample_many }; + float (*samplers[])(uint32_t*) = { sample_0, sample_1, sample_few, sample_many }; mixture(samplers, weights, n_dists, dist_mixture, n_threads); printf("Sum(dist_mixture, N)/N = %f\n", split_array_sum(dist_mixture, N, n_threads) / N);