squiggle.c

Self-contained Monte Carlo estimation in C99
Log | Files | Refs | README

commit 68e7730f24c33268a151dd197534f8c28ed376e0
parent 8f69dd1e588db16d8f59a523dbe2ebc448683663
Author: NunoSempere <nuno.sempere@protonmail.com>
Date:   Sun, 16 Jul 2023 21:08:05 +0200

reformat squiggle.c, remake examples.

Diffstat:
Mexamples/01_one_sample/example | 0
Mexamples/02_many_samples/example | 0
Mexamples/03_gcc_nested_function/example | 0
Mscratchpad/scratchpad | 0
Mscratchpad/scratchpad.c | 2+-
Msquiggle.c | 69+++++++++++++++++++++++++++++++++++----------------------------------
6 files changed, 36 insertions(+), 35 deletions(-)

diff --git a/examples/01_one_sample/example b/examples/01_one_sample/example Binary files differ. diff --git a/examples/02_many_samples/example b/examples/02_many_samples/example Binary files differ. diff --git a/examples/03_gcc_nested_function/example b/examples/03_gcc_nested_function/example Binary files differ. diff --git a/scratchpad/scratchpad b/scratchpad/scratchpad Binary files differ. diff --git a/scratchpad/scratchpad.c b/scratchpad/scratchpad.c @@ -4,7 +4,7 @@ #include <stdint.h> #include <stdio.h> #include <stdlib.h> -#include <sys/types.h> +// #include <sys/types.h> #include <time.h> #define EXIT_ON_ERROR 0 diff --git a/squiggle.c b/squiggle.c @@ -3,29 +3,29 @@ #include <stdlib.h> // PI constant -const float PI = M_PI;// 3.14159265358979323846; +const float PI = M_PI; // 3.14159265358979323846; // Pseudo Random number generator -uint32_t xorshift32 -(uint32_t* seed) +uint32_t xorshift32(uint32_t* seed) { - // Algorithm "xor" from p. 4 of Marsaglia, "Xorshift RNGs" - // See <https://stackoverflow.com/questions/53886131/how-does-xorshift32-works> - // https://en.wikipedia.org/wiki/Xorshift - // Also some drama: <https://www.pcg-random.org/posts/on-vignas-pcg-critique.html>, <https://prng.di.unimi.it/> - - uint32_t x = *seed; - x ^= x << 13; - x ^= x >> 17; - x ^= x << 5; - return *seed = x; + // Algorithm "xor" from p. 4 of Marsaglia, "Xorshift RNGs" + // See <https://stackoverflow.com/questions/53886131/how-does-xorshift32-works> + // https://en.wikipedia.org/wiki/Xorshift + // Also some drama: <https://www.pcg-random.org/posts/on-vignas-pcg-critique.html>, <https://prng.di.unimi.it/> + + uint32_t x = *seed; + x ^= x << 13; + x ^= x >> 17; + x ^= x << 5; + return *seed = x; } // Distribution & sampling functions -float rand_0_to_1(uint32_t* seed){ - return ((float) xorshift32(seed)) / ((float) UINT32_MAX); +float rand_0_to_1(uint32_t* seed) +{ + return ((float)xorshift32(seed)) / ((float)UINT32_MAX); } float rand_float(float max, uint32_t* seed) @@ -33,7 +33,7 @@ float rand_float(float max, uint32_t* seed) return rand_0_to_1(seed) * max; } -float ur_normal(uint32_t* seed) +float unit_normal(uint32_t* seed) { float u1 = rand_0_to_1(seed); float u2 = rand_0_to_1(seed); @@ -48,7 +48,7 @@ float random_uniform(float from, float to, uint32_t* seed) float random_normal(float mean, float sigma, uint32_t* seed) { - return (mean + sigma * ur_normal(seed)); + return (mean + sigma * unit_normal(seed)); } float random_lognormal(float logmean, float logsigma, uint32_t* seed) @@ -90,24 +90,25 @@ float mixture(float (*samplers[])(uint32_t*), float* weights, int n_dists, uint3 // You can see a simpler version of this function in the git history // or in C-02-better-algorithm-one-thread/ float sum_weights = array_sum(weights, n_dists); - float* cumsummed_normalized_weights = (float*) malloc(n_dists * sizeof(float)); - cumsummed_normalized_weights[0] = weights[0]/sum_weights; + float* cumsummed_normalized_weights = (float*)malloc(n_dists * sizeof(float)); + cumsummed_normalized_weights[0] = weights[0] / sum_weights; for (int i = 1; i < n_dists; i++) { - cumsummed_normalized_weights[i] = cumsummed_normalized_weights[i - 1] + weights[i]/sum_weights; + cumsummed_normalized_weights[i] = cumsummed_normalized_weights[i - 1] + weights[i] / sum_weights; + } + + float result; + int result_set_flag = 0; + float p = random_uniform(0, 1, seed); + for (int k = 0; k < n_dists; k++) { + if (p < cumsummed_normalized_weights[k]) { + result = samplers[k](seed); + result_set_flag = 1; + break; + } } + if (result_set_flag == 0) + result = samplers[n_dists - 1](seed); - float result; - int result_set_flag = 0; - float p = random_uniform(0, 1, seed); - for (int k = 0; k < n_dists; k++) { - if (p < cumsummed_normalized_weights[k]) { - result = samplers[k](seed); - result_set_flag = 1; - break; - } - } - if(result_set_flag == 0) result = samplers[n_dists-1](seed); - - free(cumsummed_normalized_weights); - return result; + free(cumsummed_normalized_weights); + return result; }