JJSierraM/squiggle.c
1
0
Fork 0
forked from personal/squiggle.c
Code Issues Pull Requests Packages Projects Releases Wiki Activity

rework error as a macro

Browse Source
This commit is contained in:
NunoSempere 2023-07-16 16:56:11 +02:00
parent c487bdfaf5
commit 84f94e2bea
2 changed files with 101 additions and 45 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
scratchpad/scratchpad
View File

Binary file not shown.

146
scratchpad/scratchpad.c
View File

@ -1,14 +1,26 @@
#include <float.h> // FLT_MAX, FLT_MIN #include <float.h> // FLT_MAX, FLT_MIN
#include <limits.h> // INT_MAX #include <limits.h> // INT_MAX
#include <math.h> // erf, sqrt #include <math.h> // erf, sqrt
#include <stdint.h> #include <stdint.h>
#include <stdio.h> #include <stdio.h>
#include <stdlib.h> #include <stdlib.h>
#include <time.h> #include <time.h>
#define EXIT_ON_ERROR 0 #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 { struct box {
int empty; int empty;
float content; float content;
@ -76,16 +88,7 @@ struct box inverse_cdf(float cdf(float), float p)
} }
if (!interval_found) { if (!interval_found) {
if (EXIT_ON_ERROR) { PROCESS_ERROR("Interval containing the target value not found, in function inverse_cdf");
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;
}
} else { } else {
int convergence_condition = 0; int convergence_condition = 0;
@ -114,22 +117,93 @@ struct box inverse_cdf(float cdf(float), float p)
result.content = low; result.content = low;
result.empty = 0; result.empty = 0;
} else { } else {
if (EXIT_ON_ERROR) { PROCESS_ERROR("Search process did not converge, in function inverse_cdf");
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;
}
} }
return result; 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 // Get random number between 0 and 1
uint32_t xorshift32(uint32_t* seed) uint32_t xorshift32(uint32_t* seed)
{ {
@ -211,16 +285,7 @@ struct box incbeta(float a, float b, float x)
struct box result; struct box result;
if (x < 0.0 || x > 1.0) { if (x < 0.0 || x > 1.0) {
if (EXIT_ON_ERROR) { PROCESS_ERROR("x out of bounds [0, 1], in function incbeta");
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;
}
} }
/*The continued fraction converges nicely for x < (a+1)/(a+b+2)*/ /*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) { PROCESS_ERROR("More loops needed, did not converge, in function incbeta");
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;
}
} }
struct box cdf_beta(float x) struct box cdf_beta(float x)
@ -412,6 +468,6 @@ int main()
} }
clock_t end_3 = clock(); clock_t end_3 = clock();
float time_spent_3 = (float)(end_3 - begin_3) / CLOCKS_PER_SEC; 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; return 0;
} }