diff --git a/scratchpad/scratchpad b/scratchpad/scratchpad index 462f9cd..5bb58d1 100755 Binary files a/scratchpad/scratchpad and b/scratchpad/scratchpad differ diff --git a/scratchpad/scratchpad.c b/scratchpad/scratchpad.c index 20e5330..5c9372b 100644 --- a/scratchpad/scratchpad.c +++ b/scratchpad/scratchpad.c @@ -1,14 +1,26 @@ -#include // FLT_MAX, FLT_MIN -#include // INT_MAX -#include // erf, sqrt +#include // FLT_MAX, FLT_MIN +#include // INT_MAX +#include // erf, sqrt #include #include #include #include #define EXIT_ON_ERROR 0 +#define MAX_ERROR_LENGTH 500 +#define PROCESS_ERROR(...) \ + do { \ + if (EXIT_ON_ERROR) { \ + printf("@, in %s (%d)", __FILE__, __LINE__); \ + exit(1); \ + } else { \ + char error_msg[MAX_ERROR_LENGTH]; \ + snprintf(error_msg, MAX_ERROR_LENGTH, "@, in %s (%d)", __FILE__, __LINE__); \ + struct box error = { .empty = 1, .error_msg = error_msg }; \ + return error; \ + } \ + } while (0) -// Errors struct box { int empty; float content; @@ -76,16 +88,7 @@ struct box inverse_cdf(float cdf(float), float p) } if (!interval_found) { - if (EXIT_ON_ERROR) { - printf("Interval containing the target value not found, in function inverse_cdf, in %s (%d)", __FILE__, __LINE__); - exit(1); - } else { - char error_msg[200]; - snprintf(error_msg, 200, "Interval containing the target value not found in function inverse_cdf, in %s (%d)", __FILE__, __LINE__); - result.empty = 1; - result.error_msg = error_msg; - return result; - } + PROCESS_ERROR("Interval containing the target value not found, in function inverse_cdf"); } else { int convergence_condition = 0; @@ -114,22 +117,93 @@ struct box inverse_cdf(float cdf(float), float p) result.content = low; result.empty = 0; } else { - if (EXIT_ON_ERROR) { - printf("Search process did not converge, in function inverse_cdf, in %s (%d)", __FILE__, __LINE__); - exit(1); - } else { - char error_msg[200]; - snprintf(error_msg, 200, "Search process did not converge, in function inverse_cdf, in %s (%d)", __FILE__, __LINE__); - result.empty = 1; - result.error_msg = error_msg; - return result; - } + PROCESS_ERROR("Search process did not converge, in function inverse_cdf"); } return result; } } +// Inverse cdf at point, but this time taking a struct box. +struct box inverse_cdf_box(struct box cdf_box(float), float p) +{ + // given a cdf: [-Inf, Inf] => Box([0,1]) + // returns a box with either + // x such that cdf(x) = p + // or an error + // if EXIT_ON_ERROR is set to 1, it exits instead of providing an error + + struct box result; + float low = -1.0; + float high = 1.0; + + // 1. Make sure that cdf(low) < p < cdf(high) + int interval_found = 0; + while ((!interval_found) && (low > -FLT_MAX / 4) && (high < FLT_MAX / 4)) { + // ^ Using FLT_MIN and FLT_MAX is overkill + // but it's also the *correct* thing to do. + struct box cdf_low = cdf_box(low); + if(cdf_low.empty){ + PROCESS_ERROR(cdf_low.error_msg); + } + + struct box cdf_high=cdf_box(high); + if(cdf_high.empty){ + PROCESS_ERROR(cdf_low.error_msg); + } + + int low_condition = (cdf_low.content < p); + int high_condition = (p < cdf_high.content); + if (low_condition && high_condition) { + interval_found = 1; + } else if (!low_condition) { + low = low * 2; + } else if (!high_condition) { + high = high * 2; + } + } + + if (!interval_found) { + PROCESS_ERROR("Interval containing the target value not found, in function inverse_cdf"); + } else { + + int convergence_condition = 0; + int count = 0; + while (!convergence_condition && (count < (INT_MAX / 2))) { + float mid = (high + low) / 2; + int mid_not_new = (mid == low) || (mid == high); + // float width = high - low; + if (mid_not_new) { + // if ((width < 1e-8) || mid_not_new){ + convergence_condition = 1; + } else { + struct box cdf_mid = cdf_box(mid); + if(cdf_mid.empty){ + PROCESS_ERROR(cdf_mid.error_msg); + } + float mid_sign = cdf_mid.content - p; + if (mid_sign < 0) { + low = mid; + } else if (mid_sign > 0) { + high = mid; + } else if (mid_sign == 0) { + low = mid; + high = mid; + } + } + } + + if (convergence_condition) { + result.content = low; + result.empty = 0; + return result; + } else { + PROCESS_ERROR("Search process did not converge, in function inverse_cdf"); + } + + } +} + // Get random number between 0 and 1 uint32_t xorshift32(uint32_t* seed) { @@ -211,16 +285,7 @@ struct box incbeta(float a, float b, float x) struct box result; if (x < 0.0 || x > 1.0) { - if (EXIT_ON_ERROR) { - printf("x = %f, x out of bounds [0, 1], in function incbeta, in %s (%d)", __FILE__, __LINE__); - exit(1); - } else { - char error_msg[200]; - snprintf(error_msg, 200, "x = %f, x out of bounds [0, 1], in function incbeta, in %s (%d)", x, __FILE__, __LINE__); - result.empty = 1; - result.error_msg = error_msg; - return result; - } + PROCESS_ERROR("x out of bounds [0, 1], in function incbeta"); } /*The continued fraction converges nicely for x < (a+1)/(a+b+2)*/ @@ -276,16 +341,7 @@ struct box incbeta(float a, float b, float x) } } - if (EXIT_ON_ERROR) { - printf("More loops needed, did not converge, in function incbeta, in %s (%d)", __FILE__, __LINE__); - exit(1); - } else { - char error_msg[200]; - snprintf(error_msg, 200, "More loops needed, did not converge, in function incbeta, in %s (%d)", __FILE__, __LINE__); - result.empty = 1; - result.error_msg = error_msg; - return result; - } + PROCESS_ERROR("More loops needed, did not converge, in function incbeta"); } struct box cdf_beta(float x) @@ -412,6 +468,6 @@ int main() } clock_t end_3 = clock(); float time_spent_3 = (float)(end_3 - begin_3) / CLOCKS_PER_SEC; - printf("Time spent: %f", time_spent_3); + printf("Time spent: %f\n", time_spent_3); return 0; }