diff --git a/scratchpad/makefile b/scratchpad/makefile index cc2fe53..cf6d6b5 100644 --- a/scratchpad/makefile +++ b/scratchpad/makefile @@ -10,7 +10,7 @@ CC=gcc # required for nested functions # Main file SRC=scratchpad.c -OUTPUT=scratchpad +OUTPUT=./scratchpad ## Dependencies MATH=-lm diff --git a/scratchpad/scratchpad b/scratchpad/scratchpad index 79ca1c2..1516d8f 100755 Binary files a/scratchpad/scratchpad and b/scratchpad/scratchpad differ diff --git a/scratchpad/scratchpad.c b/scratchpad/scratchpad.c index c6bf9e2..58e5d41 100644 --- a/scratchpad/scratchpad.c +++ b/scratchpad/scratchpad.c @@ -4,6 +4,7 @@ #include // FLT_MAX, FLT_MIN #include #include // erf, sqrt +#include #define EXIT_ON_ERROR 0 @@ -41,9 +42,9 @@ float cdf_squared_0_1(float x) float cdf_normal_0_1(float x) { - float mean = 0; - float std = 1; - return 0.5 * (1 + erf((x - mean) / (std * sqrt(2)))); // erf from math.h + // float mean = 0; + // float std = 1; + return 0.5 * (1 + erf((x - 0) / (1 * sqrt(2)))); // erf from math.h } // Inverse cdf @@ -94,13 +95,9 @@ struct box inverse_cdf(float cdf(float), float p) while (!convergence_condition && (count < (INT_MAX / 2))) { float mid = (high + low) / 2; int mid_not_new = (mid == low) || (mid == high); - if (0) { - printf("while loop\n"); - printf("low: %f, high: %f\n", low, high); - printf("mid: %f\n", mid); - } - + // float width = high - low; if (mid_not_new) { + // if ((width < 1e-8) || mid_not_new){ convergence_condition = 1; } else { float mid_sign = cdf(mid) - p; @@ -166,6 +163,16 @@ struct box sampler(float cdf(float), uint32_t* seed) return result; } +// For comparison, raw sampler +const float PI = 3.14159265358979323846; +float sampler_normal_0_1(uint32_t* seed) +{ + float u1 = rand_0_to_1(seed); + float u2 = rand_0_to_1(seed); + float z = sqrtf(-2.0 * log(u1)) * sin(2 * PI * u2); + return z; +} + // to do: add beta. // for the cdf, use this incomplete beta function implementation, based on continuous fractions: // @@ -208,16 +215,32 @@ int main() // set randomness seed uint32_t* seed = malloc(sizeof(uint32_t)); *seed = 1000; // xorshift can't start with 0 + int n = 1000000; printf("\n\nGetting some samples from %s:\n", name_3); - int n = 100; - for (int i = 0; i < n; i++) { + clock_t begin = clock(); + for (int i = 0; i < n; i++) { struct box sample = sampler(cdf_normal_0_1, seed); if (sample.empty) { printf("Error in sampler function"); } else { - printf("%f\n", sample.content); + // printf("%f\n", sample.content); } } + clock_t end = clock(); + double time_spent = (double)(end - begin) / CLOCKS_PER_SEC; + printf("Time spent: %f", time_spent); + + // Get some normal samples using the previous method. + clock_t begin_2 = clock(); + printf("\n\nGetting some samples from sampler_normal_0_1\n"); + for (int i = 0; i < n; i++) { + float normal_sample = sampler_normal_0_1(seed); + // printf("%f\n", normal_sample); + } + clock_t end_2 = clock(); + double time_spent_2 = (double)(end_2 - begin_2) / CLOCKS_PER_SEC; + printf("Time spent: %f", time_spent_2); + return 0; }