From 84f94e2beae20780f8e4c76edd40d3fea3df7a6f Mon Sep 17 00:00:00 2001 From: NunoSempere Date: Sun, 16 Jul 2023 16:56:11 +0200 Subject: [PATCH] rework error as a macro --- scratchpad/scratchpad | Bin 21816 -> 21904 bytes scratchpad/scratchpad.c | 146 +++++++++++++++++++++++++++------------- 2 files changed, 101 insertions(+), 45 deletions(-) diff --git a/scratchpad/scratchpad b/scratchpad/scratchpad index 462f9cdd009246e5d78607b0d738dd492b59d6b5..5bb58d1d9675814f65bbbc0eac64058c4b697f1c 100755 GIT binary patch delta 6559 zcmai24R}=5nZ7g0U=or!GmtPd`JKoRCb(oWLfX(8n+XHCk+}RNYXSjfP10a#j1sVQ zu>(*3956O|Y*+NzU25IMEoj4rD%Gl!1Ta`z0{9c5N~(f)3`$7&i58rFzjN;mXU7xcg~#yzcO^}F?86iTW3h(DO8En|5*0vI8(CPGbQtBhMZW-6t!#V zDlvErnU*UfNCQ0&@`^w5!%1!4iqs(He$&e&vFXUQ6tYuxaB#kgj(j#EX z-HNl8oj)SC=f$T>QZ1_=liQOnaNo4Vz)@v&BoZNY8%eeo+ycAG9k3$r$SVYkj9rb= zjL5A1FjP;ZNz&^E08&lm?M9y*958}6%2pj!rf5cc;7_T9Nmv^}o2+f0U@buAem2pZ z2DBW!>*fx89f`b-8A-n1RpYHZi9QsoAcRTv==yT-kT>*-0{coUBud-oW?@K4M!O_< zm}P8*t`r=Ai%|Dqg10v~z;B;!4g}XoqcWQl8fBHHBm%q(uH<01+#PlTaPxhbmNx`P zBdOZFm2{Z&HJvkTdQR!O7KxBWjGxoTOS6Ynq=4`gn8N%4ir6SyYdQi@h*C%yfU2@F z32y{G61$H1Sa~(3E!1C;|Eu+s#9exbFrw(!{}<$ z?M~6cF>$*0g8A3e{Z9%2P4{`IDl%S}Z6=sB+iYPjqnTDl;)T;w^f`xU&UuhhmRTmO zMq`awi5a?Ur3ePZO6*nM`-=)qq$ZXl8A5>e*J!>Xe1s;JfHWPnctfk#_^Nzsd~17a zb6AG4S>k`Prw5M17};el`Iz?5_F9(V7YQGx*_Q^t5MH!3U3Ex#W?bHlScd-%FQpzl z3ic}Db3E8>D*Jk{yJ@Y$E7hbl&WrHh*oZal;TT*5kkJwv*Rtd0xDL!CF$V-{frvokqShHGSo)DC)vq{N{&%J#QcS%rm;VfpQ81 zRW3Vnvv&kV;$;jei>`yI;n&kE>mlNkg z%x^Ekzc+Nx>X)hbyhc4ac)GMb;e{mZjds(15aor-ydkY{(ca6L;dmxV9kf(9U4nm+ zx1DoiRut?fa1hd}bqESb@smQ%FQt9om^g#4O@G+9J&7ynIg=yuM_Qh32i-=VlaX3H z2LLu(DHblY5DV7!Ut%BpqY1Pl%P|2CqO@P+k7ay6P<>HUP5=aOIe2a51mI<8*5^xIBP_EeI%oAp<-$SUFC&?}3D^9W2#b6N5WN-*)Yo1_W%`D7~3T+KwnIi8% z1Qc)ZTS9WMm$kI0f@LJwD^^*Rs}gcx)2a9zgy^7#<_1@y+&zRNyFP^9@0sTiYh8~G zC2&25@G=^D7(W+cT26~p zLuTXZK+7nAg9`WnBqo4lv|tCYmQ)QK@(;zozHgLOmAVp{=gS8+G0#!fx`kW?+yff0 zFl~PkO29mhD4szrVr|`6*zwAut#2Kn9NJvVe$wwfE_V+)8uc!OY| zT&%)5@9rM;v5IS`uyY$QID$?uMnAzQ(rndfl+Ga}iFu1FA9d$YE%OYr*4_MarZdpP z=t{R~9VM>!INb(;pz4iUKl@3D7IUarJC{VfmQh8jb0U|%A?_(=5%)|GCFT6YxKmcp zTqk5jLX?zGX+ber-P3IVyyRf=}C%ZWGTJBBy7EK25i{=toN0++s^ zm}8&X_2TO8RR&e2s^fRrVfLV~;Z}ZwAsm-yfevcbF1Ju~*5Xp!Kft{q3XlG zFb&Vx&be6H_aR+MNN{Z6ys0Q%;V>_|2?RtuG3XX+MplcvXDDT3;nV znA|pSnCC6lx>ra9hFDvQzt-EpNhw`b}9sJYhjv?-zylP6kWkxg6s?S{$kOYF$v0ltZoNTNMn8yht1_`N zh7ZELW{0^Ch3}?r&cZ`%JAHz%mQC~)si`l5kWd|C7lu3OJ1F>`5(Bg9H^!rADGr4% z%E~kvV|ZOw>UZM-Iea&6quLeTcSmV^!7DQn(|;)mzUwP5310G*`GQx=<@ULI!CqOK zeT}NUvh-}%>siaEJf3S!DEueq`g-TAjZXKh>TGVw-DE9W|HB&R#)g{u#_u?1)grgL zs=W<{0G(vG3`&*{fdh0H)PnVuPDdhc(0tq!9eApk<()pmm_b zpj$zAVLo(%4}KPjjDa?M9*I<7b98`i01bg|1(i^1pD;=}Qa>7PK`n5Y zivY;v7<3otAgIFcn>Rc8$&Vuu4=hOj8*qe5k_*52_&Gn||2HqoG4KHl5MO(6g_Iic z>*p8d*#dUV1E1sKA*o&Xxxl+ga3qPO_v2Rxp7=LOHd|`UF~Y%H;MalwEaCT=EltK^ zb7yo^lOfJ*@tI@%Q;q*?J`i1E4n?mMv@=?gZh>et#p7%wa*L62b&=W9YV?`yO;NW& zGU_g9+hGg_7Eit3+z0)B=o1Y=U#h~-PX-X!7IG}~{iI(^`a7ZD7rh8dMN{i0n1TzH zU8mVvL2?0Au;jXqL{XrUZow&nKV-y1JX*z6@I=-k@{LQ(gVDzK%?cn7N7sRtfU#DMVA0g5FO6~9 z7)+=R3NdRJUzayMFdv*M^y+Yf+l`!17hwf0f^&QkqW)}{yb1&z$tDdVw+hgvrx3sA zz}xVrNaW`PZ?W0(sIkOsZ;vW6J0FSmncb}>pES zv!%!!vt+6tdgoxIZq&au7WAz~hZHkuRzb~5 zgJwmqmODuvsH9vK4DA_uXcH0ETrD&2P}$Ne0@EM0(XATd8g2AL8a15f-eOx|r7ojo z>3A%3z7`|BI$ZO@o>G_5)xAA|IW?G8Fa@Fjkw*y_kdn`T$x z6(W*8yLtRVvCZMvEc9zg^!isSM%R)qA~`uf?%~U2TUxWGJB~T*)bgA- zJ-bjwjxS$~GU8a_vFh7Ir~7&33eU7Yah{=AX;hBtdTPbjhJ;t+dt3=O@Ot)YiHoNJDKc3TjHVP>CynbVeKbLlGGVy zopnM$Yf1zA`%$&WK9nU%4Xk-s?J0Z*!?v9lACx3*F|?|?E$4tEJY1iern1Q42q*5OZdtxeFN zEk!`k_NbxlTmMm;kYD~{0RH`2izBu6WxDoR07CI2Q2f#qF!+1&a)H#g00F){jD=Tl zr8ZlVDnf^mo%JAUX1ppwM;dy0!SsxvFO=yTRoT4oC|hnm3{LHbNSzw$SNkIhg59~= zOo@e>oTOE}Gh&)Sj4#Mv(v_zZxU^=MD>(;n4%Y3FQYV(W7jYc^fRi@vwFj zl3HNGPTm&%5#GWjtqe$(Q!4RO)4h|PVO^G+qaxsxQdcsWn{-$cXsBblBOi@`Hv=|C zV1u~_S^n@STYPksrTASo*8NjtiWc|hLTov-314PG7-+CydK4-KVTAwB3`g+n6%i@x z%tKjdTQMW^vBd7pK;DgLuMP7{F~5Ec5v>?6MEPfAqayUOHWSm{&=6_D?;vWp(%0xEuSZYkWalI1neegi9ATczta~qt znnZ|*c@R14`YaZuZ0AR8rxyd}J`X;#*DlWo zVDrz_u3ed@VL&YUIJni_dgBPF3-nIHIv<8hYR6&fGv4rumA=}bZ>4Y5KrftYMbbPU zvaSwFf+wp)XsUjXb?pKI^rD5pFwd3EU)QtoG3!2p%z~1aC@Al+?lL$Ege8j*?&Da* zfNKWcwJ1+_zcxEj@M8uhP_O}moxg<=y@xex5!$`9MWD9&FYw|;bWheDM0v8DcKCrN zvE8MpJ1)d~!!JNkT3_|6L4e=V%im*ZbP_P11<)G*Np)HH4923Q(IL(PRb3N#N?2)! zh!iAusYf&S2{_v1ate26z+K_JP1}?poDK-5YsqPdfRA;WuX-MiY!REgm`WuaSuZBS z0L%X!`_?ecJI~SEeGc`&+$X6dB101h&4^4;2z?=P2#mBO7AghiJ|Y~;6#f|QVO-_1 zX+!9>;46Gql#C51;T^mLHbAk!)6BNjMt;MQlB+Dr%p#NRlr1Vp@TL z+G8#uaO+sBUEOW|q!8oIaBD*Jw2|$B$$BbE6a0nK7V!(I8S)#2Ji(DQ?+r*uZ8t4b zcN9HGLp5}!HCf$V^b1V1r#?ay0T?}E{@+4g=~$Sv6aJ5~mHki)m_3f<9%Y`vE!6WY zYaRRvN1mezJlZV~9bwhxwU|&h2km`F);ijdIM6MXLXXWpmP-220^^||>nvuVn&*M~ zZ>|MKC{Sx5IjA*YwOJE-Pz3R0?A=t}nFPyHc|zJk`M8~bA?;?T24~SzlEv+Wb!N~upq|1*Y1vD^W*44Uxso?O0N)Gi zFlV08%^RrRV5pX=vSwfer{#K(w8G1xe!xFM@!E8@nnT1;K@2l4^YOIV#g}K0-s-;y zUz4}&Ph$8WM@JE1!{niYJ(fL= zi&yAl+`YyS6(ze=7_wQ41%yk<*QK1FNna};p3d*d*d>2i&zdj?cQo(c_9+P-> zX679y;m76_%EbwckPFuRLo}fOKY#~pB7kuK?SV_YBXjN4s!Oyzu?KSQGQG=R$;nH2 zK1Gs3H$;jixJk)1b@4fhKcR7$ifp<0e!g9qlVkyhI1}23`7=tkseu1c$qXXT$Ysq* zM6O@~uZhtKgf}b>y#17wQq*1$Uig^j80-EnP^HWrVeT<1^AS+9egNvIk5x~Aa2nY6 zpn4D`o3|eX+&=q9KpX(_eN^cMtf6(g>oZv5mm8R8h;{vvZ=R_HW$<=}+Zw4`d0)qc za}btm`x{ExH^OTFY32MEW48gyT_v3Y==Hr;oT+eFT6dDt9y$ec1EMHU1K&ZT5G5fF zSMA{sqR}tw=WOVtuUg&xpo?PFi&=#@>-;?3L^}>sPF4&w#0*jha0bVjDhkk#V)eU< z@kSNmO}!T6CiLVzAMm?p%?sAU6cPYJtghXMwMreS%appDp44Cq z)`T%5>|t$YOrL@^dtIz^4brLZF8nrtH%x15d7k@jbfzB%f~JMv!+ zZM-hmnLK(fc5vN)o)^xYDL^llLglOZ{^vsP?okMEd@3?&2HzDiem~hYC`7PC+X^r=*JJ*7}W`>zb4` zEzPZUjm?eCYn9e@^-62q#$FHLZ;;%@}X4-?*v% z%WDD+^NEnVNpa@}=0Yi+txA!S3;b54rL9$IX;8jG@|%GLjpHX03s{_OkFXaVrY080TS4finO2iyQSjHP)Ju;3zA z60j2R3Sc{+jr0Jz0DDod#j=zwDQM)ofM^HYM+|^A?A!vJCe^(E<~iwOXQELLEJ>wP z(N-kMh2KK_HoVQ{1+&xK@Lb@Dh6l}FYQ<0Cz6FlpR``u?9Ptp>WBBa^?k0iJxt}LI zaN_>~)oQD^#tQ>~0Dcg7a!c}AZ9C*UtnG5#4pXAl=Cj6^Cdp4(cgYpjUb#-dcJSQ> z)+va$I1-KCCQGQIGOMji_F0`f;_d`z+-gWWp$rPjlQvo(f&62TC)i7N$Vp0UF?IpI z55+?`O0pj~kqOKR*96=;;1&ojNy7bDY6tFN;0_T^4O$iXnM7-RS(1e9DLo3JdLbI^ zn1Tr$+aWKtJ}b+wSr5vKtzmf$;GICM(gVCRDejUClpmxzzzUFbBkw^^FOH^r$w7suq6I@43EJ8zf?TEPg!J|niBhGrU18pdQ5vB zpU#_#Typ%JzK4qLnF7e)_2rw6@stvWymbMeSMp|x4PQ3&h}8w0&Hq;7sPkeZBaSdG z!JmF8#e^Te4qtBIW2R#=wHTPBbjen!V~{ZMe3d_byCdx(j7eL}af2xB<*(iDFbPi* zUIfyFYQBvpl+H1|!o8)Aw6_J$AypfHLB-_G@H$ZG-b1@N_6Jvtn#A{%I=o?|O_p%a zrpo!$GZU^&Bjg+%iN8_Ar19j3eKZC(W)g&&Zon{MGV%N(X{g7zR^j4`LN_a}yy zu(8=t`{KEY@N$Y2Z)n%ClJp!&rM+rMDTOS;Ryp@i553B&w;LFIFAqCEyX&U~a-_mZ*zeL{EIgkUo~qUso6QWmlBP^|ku9o8??P zA6cBA7^?#8$9YT2u^;bS;z@~>YYgt=2bOrIZ8a>%;Go<$vg8w!WdlBx>7hmO@QJFM a5}u{&9-g0a@3LDG?)B@lbbjEAIsXFhLcs_C 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; }