From 4deb2044d675927730e97093001022166f8b554e Mon Sep 17 00:00:00 2001 From: NunoSempere Date: Sun, 16 Jul 2023 12:09:41 +0200 Subject: [PATCH] add more examples, tweak code a bit --- scratchpad/scratchpad | Bin 16912 -> 17144 bytes scratchpad/scratchpad.c | 52 ++++++++++++++++++++++++++++------------ 2 files changed, 37 insertions(+), 15 deletions(-) diff --git a/scratchpad/scratchpad b/scratchpad/scratchpad index e5dacaac26b8e9e1f7d06204af388e18fdc34aff..79ca1c2495f43680c09afbcd5b86d52c251a03dd 100755 GIT binary patch literal 17144 zcmeHOeQ;FQb-&UgB!jeX6~S>kq9Y3`DkvpcRr>b!0wTBdLZDCt#GB@eycgauRfkSOWRRwl_I z6sE#u@}yTnRmI!s#%oMDrre(G&uz8~8{Mw_XKK|WSB7d$y`j~)biazr?<>lB#FE~y z(i>KKOa;p;%JHO_(63SXW2bJkQE|)vLfG{-DLp%FQSwZw4=VSc*F*o4{x&GRF}E}$ z?DT*M3#Q!P??O)vt&;d`!oz%@>TkAj=vQ%MO5;tXt~at{)rz{_P;GA{n%q;nr*T#7 zsukWu%)4B+lV8*a<vp)?i|Eq4uRy8{~~a= zDU$QQA6!29ZL{DH&w{T4?#9nvb_2M{4ib;?>7(ITB?9f8_XI-Wcz9@V^y8~T44~y`gh%Wkqy}hw6 z*z4L840K1L!QRNDWS(@C`MzKzD%N$hH?IvW_uk=Mk*zQHu7txhk~F#wP>)d!h9n7@mDtSr6P3hCPRXmv14i3S? zDa}>&3PArz#@(&`>L>oeDhJqZEjsQWTv>+daXOA;30^o^V^@im-8YfGw6 zO#8K4)6?36CBH|;FOsXN5&y`9`(Qfx9USqSWTGSWZ)PAVns940(`a^*JVI&JkbY>Z zZ;Nm1DEuF&P=3ItcTsOhXUgdiNHcRS`94^61o(xwQ29DlPE)=$yx*QWz0IFG=O6gx zz*WVyLR=aGdkZT5!IGPZ+<95ALaP7E?WiNwuQA58 zrTfvF7*!PVA5r7r$|)kZT=AzS{UaZ&@sC_~_={fgpT43mfPo(>17*fILW?|Ree!Sr zYHAkA+cx?KR)2%4{?vKB+<$oWAA!ky^oe0)LU8ewl0ooATcMkyKd~E)${(_)c}`Ja zwYrhr<&=zP%PFv1M@Wy_De1JjQ@bnBSn(1XA-QX=<{2679)d0g>=ZH5Pm(uJkP91PLFnXdpUW~ADv%roJtjkA%p&wte3R1rH0PB zWvitHS4kLD_D5#@VBq6FK%lnPeAhnou<{q+wYp2Y`m@Ja zjbnV7j_y@7Wxn~boWCGHrD2R*rr73D68JL`_^L$!2a5=~LN(PBtvLDQ_+bEc4WQ$G z+3}`d%1D*(9w3L0(RoT9BmXro#Q~wzo-C&9J@^S(po4hFJdV>7?*Wn53ry0f-yn}X z&{#5G8g^BbN(vv#O?w`5pl~=C8ka$%-&;fJ4-QyM%>u11qj^S~v}aBlXS6yfawFVW zfO-cF%Ik2k>A?Hi0lJ*Yp>b6W1tAwfuR`VBi2<2&>Usd=Cl_9{oZ@mEad_I356N-dNdHOsEmQy8#M)f_8$Lc< ziVMcMiiQ)Zv5`!rB;zFwW2sWiKQ!MQE~T_HWyaLCg*eE_D@v*qF(|nX-qZ0xe(q<> zlXM0Tby9OD{sL3I#d&i+^nImcaP?o+;6nL8d+O)DjqRxqeVxA4C%8?heID@r9k+Z+ zb702!HR3e#sY74d@V1(7J5o~}sSjFG|L!vu{HuRpw8($wFOu&em}TZ=%eT!Juuh5O zFQpdGt0YN$p%fxeh(I9%g$NWP@c%gi^!}+U)E!7hBi*rhU!XqFAljq5!tq4d9qV>4 zOSpr&ds(+TlDG?Go|H{Ik_^VfAt_6`Q78s`CEFU0$Kvit)SU?S_4kJ3?(Sr?OOM2& z90RtE&d7K)rn|d>y%IYK{tR7e`px{LHqv`_CSX~Uk6Qtj)RKJhEb1Yd_U+OdVGLF z?|!KiJ$j!g+T$u(TwYp6Paz1WZyFlWcO40|l+vp|f!h+X0>7aT&?f;-SDUl?ceMGt z%l3;ki|)Ghj+&bYCjBk=9fO=yBoKvg9lxW%n?Sjc%v1OsLw`#YNZ2rbF9W+xVWj`p z_)$!j6Xtii!p<^rJF0XKU541Kp=0woT?54}&U(iKMf04lR%e-ie)03DwK>b$u-sVz zkx%fezGN6brr5MPU5AQ&&gucjI%kvPi%_a=ce;GevgY}HsJF>_%Y5=_1>}z+);FL| zg)l#`zXK*?rLZl)&H}rGgisps>j!q}eZ!ap6)-wj-0rM?%Heam56$s8Jp(1=@^k3j zv)1W`JIzj4v$M=MzXjrd4xLl57e<|Q|A2I#cCaiA2(5Op^C zob|BT>~uG?A?YjS)FkZSDMXT==ysNh zopzeAC|6Y8fTE6G2>&( zG4?yHS+;jOQ2NJ_5ZvE}vN`W4=W zTTJ}Fx9vdIR_6(f1J4}z`*f=n2flOV@6U7bIpV8%?UZEmGFN`C;C#%*O9kg`E?y>h zJah4D1m{66K2LDI<>K?Rd6|pjEX^}64xFv-ytq^J=fyR#KQHbQ>AZNQ;CUd|AAGUx zwjyTTAwHM=9^cyPz|)0IH(C*vugcCV)>enOPVjt^i+^6^%V&qE&Yu6d@|eAvtcY26 zh(*Gc7w6e+Mu+CjIpU%*gP1jZz@ zcwdtuRacGUKLalkSC8X)CC}Fx?tz}Si9C5V9|w-v;i_@J0=$^6YPNj75%{$Lo74qO zO~0Z+;@sN|Z3gb1VY9XGOZ@7*+6LSu7G9NCVc^tn|8=Gb`AiGd+Xyed`uuuC$!}K| zoI2I>K8a`D$X1?H@@bVnx2YC-CPn_(^Z!NQ?yUAq74SfZdb} zN$z{VDQ=T0&dXUQ(G?HsT|NE5khe=j%qQysJT>=r#i9wlKOXD%){CyVp3sx{(xywy zY-9=OeSt38%96l7lTa+Mvp2Q_kJ&?dES?AilY2y0jGnRMxw*HoepN$$A=;}F2?XQu z;NCzus>k<=?s%{-90(=*`u0M_!UZ6y=ZFRZZTI=^X$`dA+d`XXas@&n(DLBDzI)o& z=C;VaHh{3BDbVUy7W^&uiNLy!_07JH!1}hf`&&B$oxbLdRm%#Gc7Vo0qlF?fuapjd3QPGZePl^V0)7zyotSidTfIRwOO-^pOL+D09&(LU$xL}2rkbTJ*NM_)oYp?9 z&-c?z-DtBGUhgnQ`|Ygyo7DX^Q+gj|Eqr}otQ{4rKHr}+b+e$NcKxpa$JaM8_Mg|6 zOu2ux%**aqWwSw4Xl=>*e810hQVEie6lXquby*M*V06f=&+7uFN6}_2WW!n>0Yhs; z*5~yDQ(m{S|18J!G3e7;lX+feFx{X!;~2C5?86``6jRpc^$63n66F3_pO60;rQfXN zcsARI4Q$BuW)oxz=sZF2POHA3F+4PUv z^!FN9O+cte(A2Kat!^(iYe_Ail8ca<)*>xoE)3(bq&-5a6HH$v4 zH=nHqh{AFvH^VdjDHeSdYj51z;4vtiF1FHLPxQx1-HkSfA-! zC|mVKlaj7CDfh!RSdaN~C=D%S0FKek({VVZo6YKLjgd14jN^L_dMVpjC ww+bAsHz*|R56?gJT#kL0it69?FPW7WZN|A@maRl3pZ>zTOo>G{1sf~=3-x3AR{#J2 delta 3620 zcmZ`+eQZPN)_uUtj zRJ+!@=bm$Z=iGD8{q*O3($E2E$d}XTBD_R(Le{T3vg8-kR<0bG3}v0L>e3Tdnj@*P z^(LJ)C{c7?0?EKPQmFGIk*8_9u&YpLDbFn3q!!0EobjYF?_Jg*nLSrxZ0w#}n^%ps2r&Wz{@+0R|u z==9ED33+Hg5j|&w;SU}&k%%`%|Z|Uwy?xd=975&-D zTsmn9rI=?EeuS8LEy6?RCLVyuBuV-^>iLxq%#vph(FHXp`yzmkDwX^uhEWYHce4_WZZCX~Wmt;Z`gOTHy zhxAQUeG_tDFG84-70-K2cQze3!zQUszWZUDyzN?Iu%tyMDKgH2Z<9mkrscXnrVnYF z#@y}I-F=_AdyBj4DPt~}wQ#d-GI7{i$8ja|L12HWEg^54;2zUacq*Dm^nPGFkGMW@ z!ci_n+&H})EjbS+^iyP#ggh96<*R5U2CX|7Df>L>L8mHf z_e~`2+N7Mf58;#Zf&yR7@g_mUXf>Uvf*q!3%7hy*r?_3#ogIRP4+*s=~Z`rYSV_qz++^8 zC1jD`WHSTBGsdYa9>kou_5>{=ES_E0#-O5I1uDt4@t{Kv-vb*DTB7pX7qqK#nA<5e zK?lS@%Ll)W)JG58kq=a&0bxtfqR-#m)W4DMsVUi&1TJ)&B_B*p&^ImD{u7;d1|3Mo zxPzIY=V03Zg{EOKaDVtvw^b@`w&_ur)EDpee zp|EgJ!+G6fES^80?>)g#`oavHk05v&5F)Er<##6JwdW3CO`hjQY#_4wsr=44InKqu zWVq0Iy;nVAmaBSX1RdXjH-sRE)AAcf88Un!rFS479OIUs@}Y&ey#*<_VI;%V5_H^p zfEgB30kb_pn;ib{t{ya#2qW|*@c8(Ag+2e7L&z>SbcOfPlEa;VyipDZSjsE;m|xm@ zzHKPQ*fx3Z<;2;^03<`O02E;)gCRMD*lR!2X9u1!8Mn6c7^28n?Vmn7);1)o7vw|r zOxU-X)uG3*Z3&7Y=pG=DiEXYP_;YMi^}xf} z#@N7>V#i|xzR3z{zIipfhqCY;yQ^tKpw~dhLDQg@L4Eg- zcr3JQpg#bA@jmt{(}O-?`YBCQpg;CN)7n984>fI6qA%uuBexbNVL(Vcliv$LG*>Vl zuA}rv`HNCMw77~%QRd=H;!_}B#ss1;%$xXF5V>EG4WoS%vZX@CyuZiV2yPEMR@XX} zgQl3%*K4kFt~cM3vYfsboJ!1TTQet~eKUCM{*NTS?j?OM7 zHF`u!I3J~D;g_Xx`ffNpb}j6a(8c{_UMY&WJ!U(p1@%h(oSQ=u^t(k>QX}0DolLsv zN$7?dXK|5Z9(tol1=LmSmC~p^W)B&mF>s@(J!UT%r`y4e(1GHe^k(tgYBo+2F8?9% zGGp@h4Wak_oEw|WCw5_yIg4+O44^%_2AHIsXxI0z~c#3&vV zLQI_Wda>7kf@|lJngT;EDAT6cfkt3EyRQ!*l?X;U%jZV|KfA9xG zJHo@TctY|00F|4`O>9b^#@8`#MXw60ZeeBYtO3yri9T)h7K^9p1hg)bpNMDL4CIvX zwqAHM_U<``$M{Gjj7ON$sp|zzc*+e?@qq+gdF#YSwe%f2Ml>WMG>3t>ouTwfv70^@S>TF_Nj7K>-4Y4s8hhNJ#<|E&0*Ri8EVB0tyMB5lQaW=| z7$-)8`pZ+)R}rA&<%+#gW_wV(gDPVV$gPp) #define VERBOSE 1 -// to do: reuse more informative printing from build-your-own-lisp - // Errors // https://mccue.dev/pages/7-27-22-c-errors // Huh, first time I've really needed a struct @@ -17,6 +15,7 @@ // options: // - exit // - pass structs +// to do: reuse more informative printing from build-your-own-lisp? struct box { int empty; float content; @@ -93,12 +92,10 @@ struct box inverse_cdf(float cdf(float), float p){ int convergence_condition = 0; int count = 0; while(!convergence_condition && (count < (INT_MAX/2) )){ - if(VERBOSE){ - printf("while loop\n"); - } float mid = (high + low)/2; int mid_not_new = (mid == low) || (mid == high); - if(VERBOSE){ + if(0){ + printf("while loop\n"); printf("low: %f, high: %f\n", low, high); printf("mid: %f\n", mid); } @@ -161,32 +158,57 @@ struct box sampler(float cdf(float), uint32_t* seed){ return result; } -// to-do: integrals => beta distribution +// ~~[-] to-do: integrals => beta distribution~~ +// => instead, use this incomplete beta function implementation, based on continuous fractions: +// +// // main with an example int main(){ // Uniform: - struct box result = inverse_cdf(cdf_uniform_0_1, 0.5); - if(result.empty){ - printf("Inverse not calculated\n"); + struct box result_1 = inverse_cdf(cdf_uniform_0_1, 0.5); + char* name_1 = "cdf_uniform_0_1"; + if(result_1.empty){ + printf("Inverse for %s not calculated\n", name_1); exit(1); }else{ - printf("Inverse of the cdf at %f is: %f\n", 0.5, result.content); + printf("Inverse of %s at %f is: %f\n", name_1, 0.5, result_1.content); } // Squared cdf - struct box result2 = inverse_cdf(cdf_squared_0_1, 0.5); - if(result2.empty){ - printf("Inverse not calculated\n"); + struct box result_2 = inverse_cdf(cdf_squared_0_1, 0.5); + char* name_2 = "cdf_squared_0_1"; + if(result_2.empty){ + printf("Inverse for %s not calculated\n", name_2); exit(1); }else{ - printf("Inverse of the cdf at %f is: %f\n", 0.5, result2.content); + printf("Inverse of %s at %f is: %f\n", name_2, 0.5, result_2.content); + } + + // Normal cdf + struct box result_3 = inverse_cdf(cdf_normal_0_1, 0.5); + char* name_3 = "cdf_normal_0_1"; + if(result_3.empty){ + printf("Inverse for %s not calculated\n", name_3); + exit(1); + }else{ + printf("Inverse of %s at %f is: %f\n", name_3, 0.5, result_3.content); } // set randomness seed uint32_t* seed = malloc(sizeof(uint32_t)); *seed = 1000; // xorshift can't start with 0 + printf("\n\nGetting some samples from %s:\n", name_3); + int n=100; + for(int i=0;i