From 487de4a731f4af5726d61e6996fdaa12d8bc9e0b Mon Sep 17 00:00:00 2001 From: NunoSempere Date: Sun, 16 Jul 2023 17:33:37 +0200 Subject: [PATCH] readd inverse cdf for raw floats. rework examples. --- scratchpad/scratchpad | Bin 21904 -> 22056 bytes scratchpad/scratchpad.c | 208 ++++++++++++++++++++++++++++------------ 2 files changed, 149 insertions(+), 59 deletions(-) diff --git a/scratchpad/scratchpad b/scratchpad/scratchpad index 5bb58d1d9675814f65bbbc0eac64058c4b697f1c..b1224e6724193197c269ebb434be8725b123a481 100755 GIT binary patch literal 22056 zcmeHP3v^S*nI2mHe9C$oL~}@I0#`$SdA?^vdkk* z6XzI16(H+wUum1OrF%{jPPb>5-A%XjaT?hMVrY^Y@>n*I#*jy)5CRc-IBCe;@1Hv( zUtOf6?e3nlXOHI`>7V)L|Nr^tpMU1wnJeA^1t3(FG z!xpkqd{1XrvoRzzaKdrBYG_ObPLgZ>1U|;RF2A1dGKH!d;)MgKIRk$kW2rM z^575V!EeulAI*b*0k|C>bE*T7%N}1I{3m(vEqU;89(;Qqd;xGfKIT*lAeWuL27Zm1 znaYUA77D)w%vW8v(iaFugPX$bv0$`r<+A41NU+You{p?mzD+Hy5np@EAC38ZI%`hG z;MO+Qyvg6v;%CjRn^^l@(O8HDqaoHFj!;{TqAg!rG#rVAe2qCa>XTf|N!{^XgAQsfnH_aS5FAI<>077%u zgHP8tzl&)#mYo(@oR`P2H$=JPsxYo7MeIG4DUm!)IJv$|4u?x=6|g`lVLHEs^HVu? zCCL}EQfd@QUI*m*YZuoYRjl(89}?H!l3HC(;<9{VwT{!8Kz)~K$bh5aOiCK?@fnn{ z3kH0W0UtKtwEoD{B=E5SM3;C~n*lejUt0~h90L;QG~oQY%(XD~4Fhgm?|Te5wJp;_ z23-1x>h>CNTDxT0XTT>)5Olu*$AD(i0RwK$pjhk;xMIKu47iLDsvb1px#qn=#b|ca z1EU@o^}whHMm_Mq?tzQS)gP!`Z5TZBs3x5{t!Y{hFVl@cI&~~lrW=5Csy9=n3x7KGFPSpku%uHz&Xnl}Af5VVrc4+9 zbn5<0nJ)P0)Loe}UFg#(f2K?~I_cC+nKIqrq*LBZnQj2mDR-t!7yfkWbD1(-@YAVD zvW(uo_+{#CCFA`?d5Nh!-&CG$D$g*Lr`gyR6CJ@1&x=S<~SP30F& z<>!oL?^<=w{A)2t>YgQkb2Ve$x|%&Ru0m$*nqQME5g%pAF5KmSZK}>@M{TlAhZgD3 zG=lEXp+_)6P~Y8qotup|k`52^9KWC8A%)cEV*_=V2# zYj7mbl17dm)6t(3`W$a(iHjPDf23*2Z{w>veyC?x6$>fn6Y;*@Go|i9kk==1a;8`( zrS87?vEI}q_fe43#qOgQRb^J+MP*&_0_cjb#eh$$bd-<|s!34IRMcElZc8LrqY=~` zPXk;AFbz;~4k)g{i^@$$AX(#>sOkhMH z-oVk{0(E4Jp$*NpruT;16^D`5TGM{78*oPxQOXCll%lfqaa4kH^BBe1+aKzI zEMj7S3Yqy#>n{y^T3KKG8?YxI{gbApmHOgepg^%kl$K-~f&S2l`_bMS#W@_RQqI3p zsW_7tmHPKV-?~_-+y+0KKN9jP&V$IDRPTSI7SkuK*k(kb46VJ8u2gK(*6M}d2h?3f zGo0Kr+Z`%6-3f&AJ&uKGuTtL+TTa6`XevLU?1+O6J)8FGRR{p(yrwvRcTs6Ls+7O4 zIFBhiXJJMw<-b#QehwM5e;A*ql67z)3GIiF8t;!ER(sFdyAvTL(yw&WVI@?Ya32pT z_nuIDhwSRE-*7Vm>fGzOT}t;Ov{KxEkh*buKL%MuAoURaoz(15(W6NLhX;n&Icn0% zqy~o_-!+ahs*gEr$sc_HBhYQ+VfZVkq{-=3ZhI{xw8R|O0}VJTlI`GyyCBWkAmHpx zuGc|baiyaSAe@7SoXNSULLe5Qe-H!0#q1$fi43NdhG)Rqm$Z^R41kqfH|#^c0^fCx z61aNZlW=vN!>0EHwd)Q7K(smgpaTt31n8mCKxl9Bx{O9Uh778US6vA|?VA5J5Onf5 zHCMVrFPuXY>kqq0uO>Ztg=`x`SJ+qA>tT-=&OqJ zjm@NW{zGMm9Oq#Q+LxjC5I)1AvxsH{GL%1o7#vV`&cIUHl_*$p1hJ@WUx8A*kK&w$ z_1Fh#WfqK>bq2mWNdbInUqQz$(0Uj_b`V+-$j+Bh(eN@ndDJ9Y`DEY*#Xwf`l~nx$vb38-RHA0jPk>TQq{LTp>#h$jIN`~5Bro#03BaX)7RWFx zKoj98TE}@UOoVa6%FdTaRrwRrmSHM7Pb*GL!um60$RieF8^g_3d3y8L(a1@qyflI~Hcvup^P7-U z$}cKAe}q{^0c|&gP}8$eHWEb$p-VWALqJ<_%t-+~{_U)QZqdhCah_3}N6C??F9OM2 ze}Ei0d;iFMQ(WW*xJ`M$LXzE=Xr`0!3%dN7)L9_Ow=Yrvdf*krj~c(U#t=WdVF`Kd z+ngePCg+HsvlZMAN_k4zaRb5}6R4TPt-|U!77jj(2vA9Hjr<5O3c7&zQrYv7YZk4Hoxw_BeD+#aar9Ng*M zQ-y65r?3bYn1&jP()YJ?DCNIbcHF~lI)M@4V@bDtBA)7x!-X?pqq4J?ug;xd#f_>6 zpQdThe0DP}mwDs#pGhD)PH~~b^1|C1hSna%^?T)kH%WhYlD?@L))2=~>iJyljzhbA z`Qt&02z4jqQoIcJ*NEI|kmQ=Tvch&f@0iD6zm@F&_mTcwN99M7cS6o^OG)v#YplL~ z>WdHdCMRLTG}t>t+o!R(f7R=0xy-r6;(L<#Sj=_42k6)N$9z0LaUVb*;$u?h=pr;S zn?!qmhWu7tz6G~Km&s%Hm+cx08;+12h?glh5_j$`BJTlh9*f0?Wv2QHL0}OA&0j#3b{?Sz8yrywJ@S}Pgs@;j?PbHe9dKz!Nz40UMt-vKZ zUwacq+It-5D#rF|rDnaoae8};C7Psq-rlY5g!h)3J@day3f5M~FL~>#<8Q;ncL89v zd`^vjzzxrwPn0AtzKK=!L3~xmpMamlIY3SPk2%LVv$yxahRz;({1#MWOWyHFawc*VNV58E94FPx92P3q7rz5`+O6@MqTBDFp{AX- zcfHq_S_xUc-e5?0G}gr5<$5?XSbLMR^oEq~MVP{8znL}s9t<4DL;7bkTF&y%^8oVr z=b3ZxPXu3onSbo)4Ed)HA(qcSooF@L4%NDU#vy0)&$BpwGy3N_u;=s7Y{-6sf5z(# zW&D#lj1cQM`KM;jl1A)iy|+}ypC<>s^i_}YGz)dJA$P4Iy%*R$r79H zJ}K_UYU1gd_?wmSKYF#vFR5Mq1?r*~J6W4z)&Gd=)-7qjLLQ9sX!A zkSRw{>2J;yHwI&V-7B&oNK;FZ)|kD~-`v>I?2iQlcnCdyMKBf%M>g5pTd7rlOIvfW z-5!dzwt#){cvcmSLK2qiRZ)AWBhnZPw?^!`M(O*H(IX^%1{)bWsZ@LD40yxRe<(eg zuD*-0hnJ%KIP~41Y4_p4=DvGTExcg>Su+>)Z*@fu%)Dv-Eh_7S#urN5KQt7_}JkiE=V8};W2#5fDeJnL_E9kvB76}FwPTU zd-0hH?28=H>z{un~qWpAfYU zz=}br9)RiunNPlZnB)-)-5L2xl3zjc--rA_Q5A%iPG~ZyLK8aNBy`<8!dE?#>68b^%Q<$|Roj`FKgX{q82qLfd zshav!mFZI@^+}~Zt+09`)`ec)qosx(@p1VL&28v0{+^puj;2u$jCx?y1EU@o^}zo( z56J&-4bhDx65SS1GSe0vW)lVF7fj^wc!ei*hEDObE1u;4#PIt(o=OD8&&zm{{}-cs zVoI|9osYCuD(}U$n3B8}?-LpMpQ-^-PrvV_bfqvD&uVz0`yfjA2}xEGVyt$x4%6)h zrCK2`+<<@7=?TA=;;BLyDDgItk#;8ppZw37Y+v3BNM5xXZRUYv<@qb>{n?oT; zm5q;QGJdo12~5U$HjY^@$H!1%!R)u^#I4MhV;n6^VfLIjl;DRe=43%I$aTe7XCdYa zBckUm_@}x2KfbZnf?aYC-`pBe`mg$w%>P-9H5T?6CfBEI{0f$9eps+0lj}(~KUTw9 zBckUm>?&r=(->+sN_D{YF2wCOz^xgCbAg!FNFp|L3#+G0s1kH7wBO zjo-QY+s7HETsDr+a6Olg<0-)}uUk_8d7imEt_sm{tPYor`}0`13+XC}kT#|=@P7)R zM_ec6b&bxsa)dss%!A*O2j2wTZfsEGyEs3_S5(R~hYhgP&MuG7?{NOh6+J5S?7vp8 z<9mK&dx6`a#~jZCdGuUEc;V&q`eYvd_c(526~g{&*iE83;` zK%}Op6lz87Y?^o)mfp9Xel7CydTJbMe#pC&g#aH9Dgl-=|TuVZN7{)~kTv z_lR{vj>|+`=SYvuuIrKaKUWCcF7&tQ)$BUpq~9DrbbTc~8^jGGKRbksN8tNErz?

WvBW167U9&>pG0zn}L&`?QR{Fe%_Xc|A9RC!=wk}Ywp)R;C5q! zB0rVuv9V4u;QWjdb$!5v%d({YAow@fby&{pH#wfg$_pP!{UWaBNS3UJC*epVJ<@Qs z;q@D?3f35nwZ}T}3$8{s^1cpVti^}dWkiDX8jV1!Z&P#YMm!P;#9E{6K7Yqn*4Ro9 zH}KfQwXkA=J2w-(P9yB|N2C62K0JVlZo^A7{4GIWprfT_8$=8xA2?%KtUllJo4qTm zd{wI|=^Y{23<2h=T(`=*vU*u|4S%HwAiSu?SEULKYURz$x1whCQg4lK_44Iws_J}o z-la8Fq?Nzrqp`h%n~JjWogf~MeBB4V!h_!bF@ns&t3L1|l5F7P-gzRQ+k~21{V^Y& zI`MauWXs2b?eIY);EUm*l}}1+Y~9Ly!GIqRq4*0-vSrLKIT@*)$D*XG!WV1Bb1u9; zWu$7(A||yPo+HttF#fFTvX`G2B++m4BQpdDjMLwh0vBez+ohSmy~Q9u7aW21R$mir z!>eFOm+^6&uX;5U1;P`&+zh@5i2NSYOCozH6s|08N8h)8}a$@ZX6lq z^wON1Z^FrbJx#U-Ih3%T2NvF%V_;aseQX$J=4#*867z2ajYW0ZBy)I+PB7ZWTzIk> zbZv@sxZ0wvZNX@4n*rL`5ylIr!T|w!msZb?`8P4n-sEp@Vy?ip2%6AoEUMSg6Iwjz z%qsX$7Y#Q1NkHVEak`=E7Rmg7rx*mVZ5^%TfV8xNPzy_mz^$a}ZI@d1apHwfNAH zgOsQ1BPB`gyfpG9&l@B{x6MZRdLb{VxJ{7k*h6HeYLmR&&q`{~X1`4SR@CD^ORds= zdA=j5Y+s)5nC1fy(U z%FF(LPslG7d~%;JsXYH8d6F^9{}>qgMC!+_nVtj{w_O6vey0)u*!+Nc$!FxJh1~Ye*{3Q6;~5 z{v8tXa{M#r9~HVp!W6LP_MbAz%X1P*=~=%qncF{!N@_y3FV9mX9Te+{)FI=?~IN>Wbx z@gIRv|1%Ke0aKCv7axyByZ`_I literal 21904 zcmeHPeRNyJl^DVh6LgiCzuGul4DU@HnN-# zmIZShI~8D>Za;v22Fez;KpR@xrKK&horENX;v|#?(t-&DB%FXzfEr5B{_eaP`RT=` zdwTYtJ?l9~x-)m~y>sWzotZoHq<8Dm>Seh(IZP!JyOdF><4g@xE-MDEkQsn-R>}(T zJ)51&CZU`sa9Vx222iV|!@1coU*b)G#8*t2G2l5GnklF>BuISYCGQzpje^Xom5I+x zSy?th^9+qnK~R;PMmdJhSoA*_2Y=T% z_-n_(?;i*MPr%#oF{TfN z8GhGHZnlW|s%ozA1%i>_hEQ8H7^%79vZj`Bu*Sc>Dad@j4b3fKUt824iTZpR>6nDU z&8@6ygTJ}i&zf2`u(lf`(FPWbG_bZ%nA)nBZTVUwp>VXpSKs($7Tg?)vex!!8=652 z^({3a9iB!)rLXAlV&9b&+jV&TdVY@%Pi?ETLx-QDkbn>A z@U&*Bv{Q#KP)NWhba)JCCdGAlO9so>ULBt6@B=#hG#!2-J7Iwd3rtvG!U7W(_`hU< zllH*6mAo%grT`(8{92|=fk7&HTc%6_K`MD;rc8lP zD(TOZDd0&Zuga7uP)H>`nKA_ksiZqorobSTJWrL;yAzjE?<$$-AIhr@z{4Zno}zbWG4>RVv^Xj^a7W)s zzV^AqwXh;dJBt{5$}T}OiEAL#d&k0Zph)(|7=Izt%fovh>5!77NnZFM%G90^EQC&y zLi^BsVhuz>S_*lRcDq5^H}%q_+N*&ZmG;CZ>GV_NbDrM>@SV5^WvEdHB6P7&@q1&t z-SHN2n;#+UTv~b7KVhWd9!jg$|A=}e=Jg5Dk#LrKq(X_fYTQY_QA#j0JKrt*;#Xq$jcK^gAY%bkZ$ z@~ifpk)FN})sr#EP+fB=M)eX$p&OVSs0V2nuKG_z`JMI-Wu5b?&md97H439`7oJk%n$>5bFWxlVj2EPTZfNcVH-z|BA0f z7&scMwL+K7SZ&kLqhr-eL!Dfu&7 zM*Bnihx@G?Dr zauX`0J|quwgocE>dN;xOt>)*?rjfs{m-wh2mHW7B@1{aj7XBA1LAh}fcl7o*YzHr< z=3Xjf)_dB&$#YUz6n+Ep#Bykq;zRv<~>TY76~C%^G$#O?;Uv;I)nkX+Ka_{Ci{iv3W*QRNmy$x6>_Xib~L=kJWtfF zMbDlU3q|+6nBUx$;N8WPaU7-I9w=t)*N_iAYF|S~pW)8^{5x^Jv~Mex{-2?OJAco= z)5j|X;yxhUN6461ZR+gU*4FlF>kipsCxoFCu%CnP8hdfVH6rGru+4WwC9OZC+&8rn zM38PLc^x{>T_c-Xx$9}@1!mq|rK!msFG9IjN!tqq7byvypk)8nc>6S=u~cv?#V8H@ zmI3!cZ||VZd&ix^Y60*30%1npH4{R;x9y@KxR%DBJ6`3E1JwJW@6q~=+(GUt`Yo=D zL?;S1{u^3e_!!FEb#h~Crtih*dq4I1*QiP?{i~SBHPkx{q*WtZ+eZZV5`kSu03%R9 zP-ibrF2s};euaVA*olTKMZ>j6w8l4edVA{ZHgdRMy$8+R@f)uFH6Z&PmAtW6s(T85 zQiR2=$NV@z;wK-Z(`&?r<&C|xF8Q&L-T`_HBiR*0!*Bl?VY?ooZ|eMuTK`{7ejAv? zBaoo;SP^d`go*tPe(QkIAY|4z@T>~3aJzu$Zo&9|iiv}+RNG-La+vfXFgDR1gzm!# zV4hAs>F6)#ckfOg;x0&2m2iXU$b*@8A4ea`x1QizU!%@A$Q_ENqkX;eKF^(pMxFO| z(QA0K<{Rz*<~|{Gm@XnI_$M{5`6^l<#o&Ym7&uyyMJ-f>awB}(DuH5D2Fv`(8_{S*5UV6uY@AMgCWN zgu4naLK_>Wp|yo}Oti{&of3Wo=9e`7qsbxiqW~Jx=F_S3|XRkp+Jy_snIP+{}B|JkoX3=EAEZ; z^R6z94#NwV6YikV)d@PFDdK$yL=-=Rv&lZ)+dG8lJRHa8_uRFccZaCoE~xd8@B=si zEuB_#7Yr=jam_C792j-Q2Z_*mSa}2Q`W6rw@54HPyylPaK%uub3xqyGe^7|Y+jWlA zjrN$F2wT|0&yo^R`=F`L~1Qf%yUeQ6{m zhcM~Nfg(=BMU-`_Dg*4KNq z(A^gspc4fFoG1ojld<1x^|X9zGy?z#*a;wD9}ppAY%)db_v1RcU!$M&>GVV$1$Asv z^3jj7w9(qb@?Lz$WXRY@axl_G%ZOU&XlOsyCp6Hs$F{YmuQxf(y+3xK_f@ldKh^Dz zokg(;sxX_e@D?PJ!xWvYF)xNtYD7c>7!)U*=PKH0>g@TO=Al{PbZgAd!rUYx!N_zn zcme`jG{p2GuqhUSO=2em4?V3ju8gNvp_Me|`%0{5Su@)Uggee%AMtIEV{}edV%982 zYaZ?z;M?v4s**dB+%Y25MD45~YJo$L_zpB(i#QJZiHFyoLi{wh6_qIz`ceBjYB3Yn zzJkDL2z|QQUW&bUcLR3~@@{f%4?2R+ZcLc{*zu5kCb0reB$VjGBZsbdyPsyg`yl4k z2+(y;H&pUl;uJ^>Em!9O@~@S&A2-@C6UXt}nRv1Rzda8iG4mb0-%>=5Y zZ5J8q+WzgLgO#)~c@are+snn$P>p=RZo@Y%7N~BrTNCF&GHH)I3~)jO9{DiN+7l46 zTcf%J*d}{P;{CVD2e1%g=&;~|!O@=>)R3CsN_#PoyHFjlI}*P@6~gp4&;zgmBiL*W zxKS8L@w(?h+Lzd_@k2==8tOrp68OlApspr$v3J8~Y7y@MyBj&2<{dc9V5b#6gaM1A zMA-hsV@e#B+jkC9&&(I0eOyKO0Nn8uFoj!@)q19Bz3?k|@CDKav_P>nTyM>Y>@ghLu?rP7 zI7DpICs1+k7lq<*<-lE1f2CKVAc{Q4`Gfs z<^GA$+r#hf=dR;?+g<1ZcfH9SPxI~vkjJtX=iLV}mRpbTt@~+Y9q(xo0$j0hwP>UE zBzL?3uUd0}JCCZ6p}kEPJhnFoKdt~dcfEym5-SCF{Yf`?sVI z`gR*I-1Q9aene0V@oh7O)~|5KU-`DD;Wb3nLR5}J=pNQ2g!))|5SP(`0r6PkF`TMF zrG;jwp^@r}lPv_XwXA^HyY{068f}d9ux>QEBjlQn-N}uBqyB~w(Bk$~@)#sk3f+d^U6g>5#+g@HM&KG5K64~H6BBF(-MpPN;MHwGhZL0d}$k@%yw3ma^qwu=F# z2->zA+x?MXK+qCz7>xcVQCu2{v_wdDo4>iWDHyRew1?}Xp_VYK2{i|8ZLPs@lmr=y zh^Z*4hM@_y9Bzr)>itdi?M?n@FaVk&u9C0EYps z2qD=k>9hlI0I<%)SnKQQ^Z?M;A4{jpf$spU0o)1L2sjMb0k{JJ6FmeNd@G$E0c<^y zPS+y*-VWFV7zgYCWH@kt18@LvC!p+m@YI>AUtB*I7W833Mapc2vCiBB~*k8oa#HE78fIz5A8Bk&KP@62PdcIGa(?94ND=1j3zJ(m2+sk#4Yc_OdM63_Ds zxHFHjOF)`OHb#he3GpqpSi5sQmf}v+7eQ#c9?Zp%1^#7IBbL2j9{@YSIM_uWD~TUt z)s^9|B>v^Ze=qp=<}C$N#neWfC^Vs=8>OtPN9n3Yv*P#EbOyEJ*o)srh(vZFzc>dt zF*Ya@?it`pftxR=U`xP3-~zz)5zcF|<^E)fC4cEuR*vGws4K=H+DjclzdLg;w+!ax zzGz86;o-bSz$zftY38e#YWgsDimZbKO@2B0Y{&iKH7K)y#agF{7voWfx<=GJNWLb< z<7yfYTV^~=ujUvBg!)0gJP6+{Bbi9$_&fw$)1TAnZ;`IcEY`bnt1QJmrV5Mg&ODFB z(QWovN;)T1TFOlin5I}tF0(j17F&g-xWZzsu;gDpwGwQ{A!9G(29YNjUnLoLnU-2? zJ?OCGPIINDq3Rw288W2~3AanQL&BXB#w8q( za8SZw2{kGI^KbRP!{wJuwh%_bu~EoFxmK?l~I%@5`ot zViu2KU=z#D{7+pk#l%cZ{r@~0pT}+~-`Vt&nDXmv{A8y5KO3LVlwV}y&tS@* zvhh=x^4o0uRHpnp8;{vO#<-X;>uqD=@rW2djxZ(@tJ<-cD(La() z^cXE~VxMK||M>b^6HdE(#UWLXrT?WCW&ZD~uQ9PRnOc{!@n^BI=7kBXj9N#s>9Go~ z)MK=~iOpcvG4X1Z9Oa?qC=Wjz8O3NhSQuUiUk16TmK4}bpy8s=>s&Se>9L7SP%TAw1 zi$VVx*ne2=n`)nr3bbour~Bc*k0VFiae&_FpAU~i|2u() z|Hutp+}%Kr%Lt!)`uwG*%VYKH4Jl{GY^@V6*_YHf^ttF4+0AOx2>T>`k;Iovd#byO za|K=!=|GvJ-!3cCSr0ct+T#l!9!)swTf%M8)<{dM)6K@lm$3Rsv@P0>Uli7}(HC*}qRl?M z2O}J$_hSTFd>fit*5esRAledX^ZDC1v-%c#z=7uy&eD=a?y-sJJsBaNKN9h8@!`Qo zWD9GE_?v^iKznoZ7O?0_K2S!pNPWI#S9`8l>RWncCB6J3n;^h^m20l_Tv2sdc8z#f z2vB&_jBlw|O7K=*&3wzNS5$baeJhqNTeY;tSL3OuUP@BM>pSY(+J&kp>tFU!Uanr| zK`-dQ`$R?&^Y9*zPkVI^eFqqVEvn^=6h*S(B#(@?IgqBp~g;%GFb1;u+|vTJ-kyp=||A-%3<%y-b} zM$w07W>8E2qL^%ncwdbUu~^te_Yjde+qN`E{p$gv5e+t~9Nu3OjI=T*o(%?_8^Y~Q zn%u!ibc+tPzCDCjI)wrf<*BHe8})Bsg1ph+*2tWJEnzgF!DvLQp=X}-pi+>@f)90( zV3VI1WUjR-%A8_coX9#iw4g|twqQMTMuVGCq*)K@7UA&DV56MpjR8ngB~6Al;WZAG zM^pah5QJ&$7+i*BthpHr$muKMKZcEiI|@n}X_|K~@V|l_z)sK6Ip4{H0iz zDJf_ZrO_{So>v$OL-qW%l3zi)B9u_=UkXj}xSn5~8x^!ELJ5uh6c^wxM=eVK9kP7| zRr~7x#Hha?{B(Ds_|-Wy1If5Zq$6b&zq;oL0Y^3l;;>96P%yajc1m!it*{HtIGF(WA(EB!SrV?D@_O%=a7Co342gsOeTulj$t zO8HWx?duG;xY2y2OQZ%$w$bfB?*luMC15>*T8R+dd2Pc?$l@z&|M4Rj^1I+Q6^GK4sun_Xi59drC!bwBK&Yuf|`+ z;R;$a6`BO4x03g)fnVK6DEPP{l+f7z5d*)vw@}b4=d%){{703KAwzcqicFnb)%ydD z^s21jyWr(IesvGBQXZrf9Y{n{Wd-THOunnySNA3AeUNHm(fp&N_>>(V29Eqc!@q*4 z@IJ{U8hl6tnQz6f;3P2X`7`$~M+s1Ts$2kO!m0MvIlx-3Q5fl|Eq$6M+b_{0 zkyrfce5tAxzaAP~OpunBYaBbI0(4%aHk3YUTATr#(p}@wz6ZS;e9lm*+EuhP-3%n6 P>ZjEj#S8<3p~(IlzLmNH diff --git a/scratchpad/scratchpad.c b/scratchpad/scratchpad.c index 8e070c7..8280379 100644 --- a/scratchpad/scratchpad.c +++ b/scratchpad/scratchpad.c @@ -4,6 +4,7 @@ #include #include #include +#include #include #define EXIT_ON_ERROR 0 @@ -20,6 +21,7 @@ return error; \ } \ } while (0) +#define NUM_SAMPLES 10 struct box { int empty; @@ -41,7 +43,6 @@ float cdf_uniform_0_1(float x) float cdf_squared_0_1(float x) { - float result; if (x < 0) { return 0; } else if (x > 1) { @@ -174,7 +175,76 @@ struct box cdf_beta(float x) } // Inverse cdf at point -struct box inverse_cdf(struct box cdf_box(float), float p) +// Two versions of this function: +// - raw, dealing with cdfs that return floats +// - box, dealing with cdfs that return a box. + +// Inverse cdf +struct box inverse_cdf_float(float cdf(float), float p) +{ + // given a cdf: [-Inf, Inf] => [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 + + 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. + + int low_condition = (cdf(low) < p); + int high_condition = (p < cdf(high)); + 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 ((width < 1e-8) || mid_not_new){ + if (mid_not_new) { + convergence_condition = 1; + } else { + float mid_sign = cdf(mid) - 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) { + struct box result = {.empty = 0, .content = low}; + return result; + } else { + PROCESS_ERROR("Search process did not converge, in function inverse_cdf"); + } + + } +} + +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 @@ -277,10 +347,16 @@ float rand_0_to_1(uint32_t* seed) } // Sampler based on inverse cdf and randomness function -struct box sampler(struct box cdf(float), uint32_t* seed) +struct box sampler_box_cdf(struct box cdf(float), uint32_t* seed) { float p = rand_0_to_1(seed); - struct box result = inverse_cdf(cdf, p); + struct box result = inverse_cdf_box(cdf, p); + return result; +} +struct box sampler_float_cdf(float cdf(float), uint32_t* seed) +{ + float p = rand_0_to_1(seed); + struct box result = inverse_cdf_float(cdf, p); return result; } @@ -294,83 +370,97 @@ float sampler_normal_0_1(uint32_t* seed) return z; } -int main() -{ - - // Get the inverse cdf of a [0,1] uniform distribution at 0.5 - 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); +// Some testers +void test_inverse_cdf_float(char* cdf_name, float cdf_float(float)){ + struct box result = inverse_cdf_float(cdf_float, 0.5); + if (result.empty) { + printf("Inverse for %s not calculated\n", cdf_name); exit(1); } else { - printf("Inverse of %s at %f is: %f\n", name_1, 0.5, result_1.content); + printf("Inverse of %s at %f is: %f\n", cdf_name, 0.5, result.content); } - // Get the inverse cdf of a [0,1] squared distribution at 0.5 - 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); +} +void test_inverse_cdf_box(char* cdf_name, struct box cdf_box(float)){ + struct box result = inverse_cdf_box(cdf_box, 0.5); + if (result.empty) { + printf("Inverse for %s not calculated\n", cdf_name); exit(1); } else { - printf("Inverse of %s at %f is: %f\n", name_2, 0.5, result_2.content); + printf("Inverse of %s at %f is: %f\n", cdf_name, 0.5, result.content); } - // Get the inverse of a normal(0,1) cdf distribution - 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); - } +} - // Use the sampler on a normal(0,1) - // set randomness seed - uint32_t* seed = malloc(sizeof(uint32_t)); - *seed = 1000; // xorshift can't start with 0 - int n = 100; - - printf("\n\nGetting some samples from %s:\n", name_3); +void test_and_time_sampler_float(char* cdf_name, float cdf_float(float), uint32_t* seed){ + printf("\nGetting some samples from %s:\n", cdf_name); clock_t begin = clock(); - for (int i = 0; i < n; i++) { - struct box sample = sampler(cdf_normal_0_1, seed); + for (int i = 0; i < NUM_SAMPLES; i++) { + struct box sample = sampler_float_cdf(cdf_float, seed); if (sample.empty) { - printf("Error in sampler function"); + printf("Error in sampler function for %s", cdf_name); } else { printf("%f\n", sample.content); } } clock_t end = clock(); float time_spent = (float)(end - begin) / CLOCKS_PER_SEC; - printf("Time spent: %f", time_spent); + printf("Time spent: %f\n", 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(); - float time_spent_2 = (float)(end_2 - begin_2) / CLOCKS_PER_SEC; - printf("Time spent: %f", time_spent_2); - - // Get some beta samples - clock_t begin_3 = clock(); - printf("\n\nGetting some samples from box sampler_dangerous_beta\n"); - for (int i = 0; i < n; i++) { - struct box sample = sampler(cdf_dangerous_beta, seed); +void test_and_time_sampler_box(char* cdf_name, struct box cdf_box(float), uint32_t* seed){ + printf("\nGetting some samples from %s:\n", cdf_name); + clock_t begin = clock(); + for (int i = 0; i < NUM_SAMPLES; i++) { + struct box sample = sampler_box_cdf(cdf_box, seed); if (sample.empty) { - printf("Error in sampler function"); + printf("Error in sampler function for %s", cdf_name); } else { printf("%f\n", sample.content); } } - clock_t end_3 = clock(); - float time_spent_3 = (float)(end_3 - begin_3) / CLOCKS_PER_SEC; - printf("Time spent: %f\n", time_spent_3); + clock_t end = clock(); + float time_spent = (float)(end - begin) / CLOCKS_PER_SEC; + printf("Time spent: %f\n", time_spent); +} + +int main() +{ + // Test inverse cdf float + test_inverse_cdf_float("cdf_uniform_0_1", cdf_uniform_0_1); + test_inverse_cdf_float("cdf_squared_0_1", cdf_squared_0_1); + test_inverse_cdf_float("cdf_normal_0_1", cdf_normal_0_1); + + // Test inverse cdf box + test_inverse_cdf_box("cdf_beta", cdf_beta); + + // Testing samplers + // set randomness seed + uint32_t* seed = malloc(sizeof(uint32_t)); + *seed = 1000; // xorshift can't start with 0 + + // Test float sampler + test_and_time_sampler_float("cdf_uniform_0_1", cdf_uniform_0_1, seed); + test_and_time_sampler_float("cdf_squared_0_1", cdf_squared_0_1, seed); + test_and_time_sampler_float("cdf_normal_0_1", cdf_normal_0_1, seed); + + // Get some normal samples using a previous approach + printf("\nGetting some samples from sampler_normal_0_1\n"); + + clock_t begin_2 = clock(); + + for (int i = 0; i < NUM_SAMPLES; i++) { + float normal_sample = sampler_normal_0_1(seed); + printf("%f\n", normal_sample); + } + + clock_t end_2 = clock(); + float time_spent_2 = (float)(end_2 - begin_2) / CLOCKS_PER_SEC; + printf("Time spent: %f\n", time_spent_2); + + // Test box sampler + test_and_time_sampler_box("cdf_beta", cdf_beta, seed); + + free(seed); return 0; }