From 0f5f68f21f6a38697bd4fa863e5684b74cd5f80d Mon Sep 17 00:00:00 2001
From: NunoSempere <nuno.sempere@protonmail.com>
Date: Wed, 7 Feb 2024 19:44:02 +0100
Subject: [PATCH] add possible gotcha to folk wisdom file

---
 FOLK_WISDOM.md          |  59 ++++++++++++++
 scratchpad/scratchpad   | Bin 26992 -> 26832 bytes
 scratchpad/scratchpad.c | 166 +++++++---------------------------------
 3 files changed, 86 insertions(+), 139 deletions(-)

diff --git a/FOLK_WISDOM.md b/FOLK_WISDOM.md
index 6a29f9f..2dc9537 100644
--- a/FOLK_WISDOM.md
+++ b/FOLK_WISDOM.md
@@ -269,6 +269,65 @@ The last commit that has code to sample from arbitrary cdfs is `8f6919fa2...`. Y
 
 - Even though the C standard is ambiguous about this, this code assumes that doubles are 64 bit precision (otherwise the xorshift code should be different).
 
+### Consider this code
+
+Consider sampling from the unit uniform in this manner:
+
+```
+#include "../squiggle.h"
+#include "../squiggle_more.h"
+#include <math.h>
+#include <stdint.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+int main()
+{
+    // Replicate <https://arxiv.org/pdf/1806.02404.pdf>, and in particular the red line in page 11.
+    // Could also be interesting to just produce and save many samples.
+
+    // set randomness seed
+    uint64_t* seed = malloc(sizeof(uint64_t));
+    *seed = UINT64_MAX/64; // xorshift can't start with a seed of 0
+    
+    int n_samples = 100*MILLION;
+    int p_sixteenth = 0;
+    int p_eighth = 0;
+    int p_quarter = 0;
+    int p_half = 0;
+    double sample;
+    for(int i=0; i<n_samples; i++){
+        sample = sample_unit_uniform(seed);
+        // printf("%lf\n", sample);
+        if (sample < 1.0/16.0){
+            p_sixteenth++;
+            p_eighth++;
+            p_quarter++;
+            p_half++;
+        } else if(sample < 0.125){
+            p_eighth++;
+            p_quarter++;
+            p_half++;
+        } else if(sample < 0.25){
+            p_quarter++;
+            p_half++;
+        } else if(sample < 0.5){
+            p_half++;
+        }else{
+            // printf("Sample > 0.5\n");
+        }
+    }
+    printf("p_16th: %lf; p_eighth; %lf; p_quarter: %lf; p_half: %lf", ((double)p_sixteenth)/n_samples, (double)p_eighth/n_samples, (double)p_quarter/n_samples, (double)p_half/n_samples);
+}
+```
+
+What will be printed out? In particular, consider that not all floating points have the same density. 
+
+<details style="border-style: solid; border-width: 2px; padding-top: 20px; padding-left: 20px; padding-right: 20px; margin-bottom: 20px;">
+<summary>Click on the arrow to see the answer</summary>
+The p_eighth will be ~0.125, p_quarter will be ~0.25, p_half will be ~0.5. This is because these random numbers are produced by generating ints and then dividing by the maximum int. There may be additional gotchas here, but this at least ensures that intervals of the same length in [0,1] will have the same number of samples. It's just that those that are smaller will be represented with more precision.
+</details>
+
 ## Related projects
 
 - [Squiggle](https://www.squiggle-language.com/)
diff --git a/scratchpad/scratchpad b/scratchpad/scratchpad
index f3010c5a74d179998e74509811bd7296b6a10383..964534600bfc838b207052aa971df797c335c3f2 100755
GIT binary patch
delta 9724
zcma)C4_H*kwZHem;)=3&QIJ0hEO@~M1QsJv5{;}bx+oe^@GtU5V1)z(VKK2auI(l^
z`w<i~ZDLZJPy1fmB(F`A*w-4{CX4@L)PT_hqDF!-xrhp4RDve^elzzjqA%}#8NPeZ
z%$YN1=A1e6o0+@xp5j|t`KHjorqLWr84BU_+vEQ<+~lv6OxxR2{4XUag-Z!!6u(x`
z`jB7px1w8SD%=OD#4CJ>ydr!YB(&`GdFyvM8{gRZZtLHEpC0e$S(Q&V`lN;Mzo25m
z7`#I8w3eN$O7CkJN8S-a`7g+qJ|U#WH<TECO?(M0@hjN-IN9x+!oNz^`KFK--%x%x
z<G#-1O`1@C4a*fF_jgjPN#P5~cFi2kztKSqG>;DaD{_HcIA<tA21!f0GzFtmO<ZO5
z8m?|#t(~i>eu}HBSW8y;853&i?Iq>KwOmPYWo2~<S5;P3QnP`pQc899(z4oG&h}Jo
zg}sdY!*90Zq2$Tjob0TbvkE6Ar%tBo1^m&s5B_Lg;CxZmfchb))cOwA(H0!XgHK(k
z$78}%A6_w>gY-p}=i2^~ABQ$QDAU9v|0Jm1!HOJb52FqIxB<$<^J~^X!}&ogEh`EI
zrjRx<uA*=FLLppDkdiSyz&ELUZIM#%mG_p3UOSTaAraaTei%v9E)MWISp)eH@{D#4
z?@JEhttA1&L-^r)Vuyd^_pZKFaF~~VK$V%7WTxIS=Vhi|HRoogewCMcb@4yBFK=qB
z&MR6u*W>YIxsF)(g3{!Z<r>yK#@1wZ9_G!?LwFAIAkJ>yJ~`Xvkq^S!K2e|aY!b)W
z#1#3dPf<$Sq?CHsKRu+41C8pB+r*x;cuwK@5Tdb_Oti<{CnlV?Hl}?)k>fs~^U`M)
z6Hf1?^GI_i=<9HwSa1%{K|H7N9GCy`gU4fa9@M6~SFSSWn^&0&%!LPRVw8mqYC5M3
zqksJ<I_zYOyDJkm4GILv=q=TaaUT_5cJqh2CZJY~dKR5>95IX;vm|Z;$5|T#{u)Qg
z^w-gyCpadX<GyV|F$)=Mqh=;0>HYR#YlG&8dB{KKwh#Aoh*6uVngc3@d>E5Y8<kC7
z7&BHlGmpG6CL-Y~T5$D=c~<B2cWKt@x?~?|ZJ7Rd6y&bE9*<m#P7iAi%_F^IBvKmk
zknrD$<dqOf*p^5>2zf#npGXozqx@qa+9xj82#F*&bo?+1CUVhTvMDr2n0}YEhRzjI
z6G%kZPJzUee}={QZyQkLtK-QxVdIA_iiiA(0TLKqApCBCtP6iyC^0g@!tWHp8Ac`;
zV}#)Ez9Z(bVSf8jWosg)@EJnYcVze21pjsL7iK14`vh`k?8?M(thoGGpBR`e>EW2C
zpTWGl`>bLDCVcsRd2GL;{7xKMIZhgO9)6aW^^+~*mJ3P!<jZj}V=8V_34OMtg%V%M
zCK{XNyxSx!B04yM&Pr+v_!e5bJm7jPrfpqF5GjdB7v7E~uSJaa|F2t04-K*8qX;Rc
z3L>(s@GMd&TCujvd?mz@VP#jqZ4wx{K={ioQXDxxRF*hS{N-k#OP`n@(1{X;qUn^#
z&d507O^F<gj26~N<eSJyaSkfw)i;&SW=bR|DkeUTHK0Anzg221lAak8N|yjSKiVfQ
z3lJqz7<DZsYCKiHyp39NBZe*~<0r_&zUBMzXhO!5YD37_uW!&kU^3N2Ri9xDeZ+|v
zvfB_TJQG8j4LMOOph0mvx@eZ?eXqKawu+b@LsFw>4Qq?Wm@oE`pGPN*S_3Uz&oU{@
zl~!JQgZv|UL&P+QK+Q-=o`e9ix_Z3c)@!1p_soxJ5QNFS<a6mig~xix=2(|7D}s!O
zTOfRLm8^_gDfnFjb2_;hXVp8ddOX>#&TQ8;tLunde3i_f@LBXPXv8-L7muau5X;|<
zj7+Y0JPZTmH?ELOW2~^hn^YQ2LQ^++&A3uL(5+Sjr{wkBL>nJ3-03C{#a|T0TqdCj
z@%}5J1!kt*2q*Ir;)FY)q%OfKyc0??l9R~ggpyIW!#M7SxG8u}I1@~|5Jn0UqlGIW
z<d=yPh1bH!+bI1u1oGkJYT{A{3tCYhyb$fOT>Y*-pcT+)x#j{kR_f5tLfL#BwTZdM
zY+`cz`iBtX6M-xW$jK11iFqwF@F8Pt?-Td8`606(4FGBRAKAnu&3)p;{%qIbt|zcG
zZOL97x*<jqG9`FS={Fva+4()cDb3uN`9rqrhP5$i&lvLl6oc^oH>6|ASjPtRBJYIQ
zN`Or@3MwXcX~7$!Zor&N2bSf!n+H2ix^6B8^k={f&1X0xP`^xoIn=BVP==Y!7-<Ja
z+6mLnt)JpG?JI~(KOW31zibm7c`~&Rg2r!QUoW-q5L(@aeP6>qz_vbi6EuPk2QL2-
z^a|?HH2M7cG^+I=#&s3!yL}dnrk1JB`xcC^sb(2kzxT75#oT_xbPrTsqAqh?&USr`
zuyiOGQ7Y->E(R-lWWALZ$pG2A^2bUIM5&!-IvT6E=>RqU3Ooh~Kf&h8c$wL471y?*
z`+LR@&MPy%6?--2Rk0gU^8XnAvZD4U8UAxc?N2pawc8<u%6$r2mP(=b49ALx;cbIB
zN)4Zj(Xf>a1k`Iz<7^7c&!|)RQ!H;slR?WB!yW&}l>bRl`;#npDr$ea<wjUhDw*Vg
zTOJ49Cy>8qIW`4Y-oelwI|yW#2lS0nLVPYp1)DE2R~vf~1#G4f>bY2y$!*FCMv)#(
z^s(J&276m+9J*mIP<nqmjc_b?wt^4X9O9Cr=<qvl|8{ArT!y|~*RoyTFlU?Ovc92Z
zqkINsdhTLy2mH%~{fp7xM#%`>NB7e68eIz(X)3sO$%=QSlz}&<81H-<*7{t@gjL}%
zy5&D!WCx;!^f+|Db`L~sO@NUB=KwEmZSlboaH)BTInSJHc72;i7r)i{18N#jV|Cqm
z-<N`N$$;yO90GM&u0z<k#b-Z8E&B4YHqO(U8!ZwxgA3M1Uy2Pd8~c$%+dWo{qB6YI
z*{lV>w@+MUXB#(r7cyCEdUO3*2!RDZ_x{YkUQ0K@#z5*|>g;^Pb7LkZGvScSA_d_%
zLG@=#*fZ6G-7pBTZM%SN;S=5IheV6TAq82TN3_-*W~b#jj<=?NSAW{MQ3?{oEqAbk
z(`k{>!>Q)Wg0;ERH=76WP=z<dYR>9(>uqB7QJZMPgr3DfFMw9J+QjY_YkEsPj;z(k
zTxd|&u61_aMF}D{1j<=7xd7UI6!TBz>pu5IEjTtvoNYlZXgA86+!?nH`PkE}CCAJS
z*^ONKQ8ByOxt~wJ{7m$Z@B*wQr_2qxzv0qb#BA*U`m`$CD&PpaWN`LR61O$ul-vOS
z;0SMppNuv~6SWSKmDE)|9(IDp*%m-)Jb}<Z12^Z=*;saL{V1w10t!)oQ~cEd%=kWf
ze|Qf{IkOAHHtsjOKJ<pA36gu^L+sT_(~y^Ee(CX8(o3bFy*M{=`{)%$bY!+5f`d?<
z>ds76ckpls&S|j!L77|S5RU_q0Sl$=xwD%hl_{{guCcCC-Dc+zUcS^pwcn|iEv`PO
zJS?hrf>gJB3{_c(K0X`08|`beTnD@Slry~pbBS>uV!?vt5{P9XKM8r*8Zisu1sZz?
z!vjSX-r3SrtX=?o^2Mhf#fXlY-=i1u+wIQp{LSLyA2#Orc3sDr|K2*nA<>%c>bZB@
z=z8@6I$(Rf12-EeP5V&Cx|rUIt^gVR6f)3<@gg9jgF?m;bczit{~UI(1KrOnQLhi!
z&W1)0#%jPDx(xJw(-1H$89t_texFVcHutyS5jvD1<G6Ui!DP?j&a)54?m^j=knK8*
zQJNbvxqVc@bMC!(b7nCslk3`oBOki>18p{bNVAP-*1c1~4f9fr@k8?xdRp+phU?5?
zU<Jhv*sP%gav&4&12UjVFEPqa4hW{)PvGOr_YOrgN|XL0fcdA5`iBNwr`JEY!R2CS
z&p!?S(35Y#az$yxZiP-y2So#{Hl}z$QS87%qxVYVmsI;~MSDBSAF&>wy%XN2nwL?{
zxRgxnx41g3t}j@`XypUF>KR!($o7>eqRVfv8ep6e2QYE#Y9OT9dDmaunuMIS<a-90
z4*AA;_!4VT%K)<nbr{Jsr~q0U;eoq9)7R_Jee9bfX%#juaM=kFngOS&xDn&3E9T=i
z(g1v9FiY2WW9|V8iBN4obl{+fq3FAV)T<@E+yeO=T5I}%eFAy|KB8MJ?f$2}!W+uP
zEp9Ad3XjfaAmc!NJCM;JxUd4O4ZvG*%Wc5B>$24aJnM0*5NYL)c@a{tA|z1J29dMy
z{RL18gWdZTgarN_y$>QGQw{&yFu7+xhdgLC+IQoy`%{oIo%sOJ*f4VIAXuK+6!YWk
z5mraZsfM7Jyl^=^YXtN0Mk&<UZ`>44&4OD3<{~cOW#qlv3b>V$9)3!SL>+_1h>>Uv
z&Ux`yzlPtASfpvg=!*G(Vx*a^LbURg8giv+nXc=Zu6}$ZaHOIU;3@@#p&7og>{S+Q
z(kSrc&~x<L1BL+@#9<gPh%KDKU|7>o1U$scSI<%e?4}6V36&=SODaIhr%(mlHy|!7
zF1sWFlK+NktIJ`*F<|i8pxSGf^cbDRHNZN>R7+k=FGnMW>A-mm2G<!}g7NCbia9A)
zu>q7yh9LmY@O=iFQGCC}^vIczvtl?}od0{=vl}Mvq8$N)hgcWEf@r9Phj}@Gb#Q}P
zco`O6!?!fAg?*n7!FZGXVT$p|3dRHC+48{bu?w)e^-Bip4yVI+3*U7Iw^VgmSg?fX
z+C%um{?xM;409UmgUw2IW5CPk)Z)61?q5VMyD@exh-^pqiuV%Ou}if?e(vg~a4{-_
z!)-$ysq+9Y&xUmR`TENZxejS6<3n8x^kB0}M4h!^Bk<sUtgdTkm^%%hQuZYb!Q4mb
zJl-!F%2WgdQqiHP-h>{yfvUzs@^7FGW(GSN(HxniSWwhpd9urB*HhDQgq+yB7i`Xa
z`wO%42Y-wBc(1u}wqIuBCSRb+w62G+4K(mWw+Atc@7h5>Dh;{FCJkU2U{Cwei|mBl
z^OzfUlkMJyEu|3Sy>Cl-7-Q~=^pH!*bA^C*lA2N={H2X-Pnk5O{WE3)JDn&WdG^bP
z6y8YSn6VX#aR@93cpcwl5u+`LuVA|$8AyqFco)vd3=pxH@+mapOOrgdO@(q?ez4i2
z3W&*$KG4=I-~Eg%GL27In%y}4pD}RCk}TJ!=DaM|O?o80n};R!XRsG%$A3=`=W_L#
zJ&#Xar+G0$TV7L``mnwH5kqXHZK9#3u&iQDxxM^RZ*g6Hajm_qcCfm<xYDK;;QM>&
z?HNkmt<n~ib};(vXD^#cM@Fdr-%qR=z9`2(LlN$sLA2(Hi_k9@YvWdND`OLq6HsD=
zi!FJ?z!_po{khJUp68M}pMQ-r#wOIM)$yoIRC#KnI?7{-EHFRfFt1<3M(V%VUaDja
zX|b_NF4g4CO&a8-4f3W8^3n%+>4Q9_3yxb*U$wfd)=+ITtgcvFXQ-;Lvl~{I8P?R6
z6~ki&hrPUbtsynpe~4gRaaB!a**)2EI2)pal8Wm2WmVO+8w~8%f?|6`^;(0ix>o7L
zR$Nh8RtlXo-5bD<=?xrIY*OXP$L6q%S6`v&y%{Ht({x3hCRwFf8B8-?Kqs0UO4*X3
z%-jHGZ00xz{i%i?trni}c$OhwbJFA447wR~2WZpB9?#pL?iP<H1N7}yj|aa_aE4PJ
zPY|f-6OYFPx(n0-+H%_CDFWp_^>{XdW?-f40Cgki{{?EmK63@M26Utk=ipk9nTY*8
z1Lu<r(9NK^pxhT8PdVjZ@OYjDt-)cPUJv!S_Km~orRFQtgKoa$@ic)peeLnkPx5+v
zPdW~|3sj3Ou?b(~GC<uu9*-MTk6W)1E42nR!_R;bBeM()6Fx%JfHr|X56a<dOp70R
zVpd%E({MrxIGkk|&cPff;hBu5_apM=tT=~!+~Z+<T8D&dJ7>q!LHY2RH{laI_X3_a
z_}WPMvvvAy!W>=bR-a6rq0u)}XM90p(V4dR&C*SE4$IPI2!Hq2O`WAPnRUjQI>Ss|
z=uDk{rcOHtJ2n#S^Bm5$-0K+MYR2L?d>!L@@%UoGG_(_|&LF(w&l=GNBf;xO-IK`E
zMCA?O<zixE6?v1;s>sv2Y2cM(vVW|O2F%7!$Vp4Mu<JCrY>5%>o*}x-WMN`Ec`P$h
z$ZsbFnURm}#3oQi`$D2a7>D`q2amdtCd#8RBw;_<z&OgVXgoE2PWEMn3x7LHKFypf
z>^M(iX3rI_d_mTNa2;gd?6}8dT%(53k;9!Y2v%L_7N2=KlhCOP&1Bx9>VJ*MX+Z?7
zSJh7ieI2}Mog^-6l7o)qD-@eLJ)Vo8LsZSt8JxZrU6IdvA6@x;U5$W=FzZ6iO!Y>H
z(qObA_>w`n$LWk>p0@fxtkHMo;EbB-j9PS?d7n5R%&Fx7CiO68^*E;WL^US%BouzE
z(`Jqkp4Lq@-&cLFb~O|@hi&nDLDT5F)n}W)m7>p`=(xWNS6LdgETm@B0nEOC04${D
z&;g9o5Po*BPv(exN*$Rrm}toQVs&(2l`KV`tJWDxbXvwA2fh)Pj0gFX!LI>7kIK_1
zv+A@|Pd4~_z4BJf#B80`0)8p@+;ue+*f_H=PEwmSN&9Eo1**vPg`@YF{40C{2Y-R2
z8j<+zyGN7k`I6A~2&tW45Ur)gBB_KsN@op8!BhOgORmo^jNXC$mnNJHSaPcwmCd!0
zbqi+j3&}eRqWJJVtqY=rFt1{@1tUKAVKkgX=LGu<=N5v7@NQNb!R;V9Ios%%I&+fw
z781B9gzqQui?;CwvU^c5Ik{+VmK&eVh9-Rf=LE!j?*n|?qZ!}(8>_bsbsdl2tBs_1
zaf(nAL|$B6AY@3SZ}A_|<14un_)5~88zr=05V>(-i?9+#L46bL2Z>=RXd$_t8=|qJ
z)<Ay-4|@xX%fR6`gCZ(*FEvrQu-`*Yk6L}vy&9T(gVs9eml)bd3k3+IG?Zz2gz|+z
zuI&+}?0uh7<@NYUgQj{_{+!D1h*kJ2RQ^?!&#8J?<UwR4#*BUwxLxT=0~UJ7CMp$N
z4L%*y#J~zPoun=e3CmJzd(~KCu!!1xTFd{Oy|>bL>hD#_76Yo09#seQXYi>z9(2}w
zwAG7CwZRLGN&|*N$eu$3#n0Wy1qxxDwn&lizWZ6XBgploA&zhaFZHM1pwv!Qg(i6U
zF$!ZMbkX#%%Fj^U;Pp?A$}fsn<h*FN+RIn{S*BETwT#b^BTEAvTU41{Axfp!g)gi8
zrb!B)p$xKL<Iy85O+^ad`+gT=WbvnlH`_N+xkuHLfs;2)yy{xF!e7C4P)tVh8gp80
zu+YHR+y;ek7r>_-KNy<bs(cZ?hSS7$SkPY5yfh@jA0SDGRg=IfIfnAcQuH207epHU
zTtX!}ROdaf2y*EvzlbhSB!=^#k9%8cioH3$p71L|!hWH);M8dJZnaG$bw!BgO|Nv?
zzf`$hQA{!CotDEYzc*6hXDWm{1wQ*9mWXSrM26ZjN0s=V@~JiaMjBo*0h$N?*F}|!
z2A|rxIcU$Ql?!;XF+YmDl|PI4t}>B@tHzRNR%MfOtMc%FH+<n}l3EZl3jZ9FYRbs)
z{%^~q<TNs;AcMSGkeFCpxu$G&ZE;~uZFOmViM<e_YpZMVFEVvpacyn!hQc*v_QI;N
z(u(4>q^IEFu?n-MwqmWlFx^yGQc+l5QD?7SQ(Ig`EQQMght}8Gi|uvf&BDm!q1-`v
zuI%Xw`+cMG{xi<;!Y9d`qG=lTA39`ne()Yw(HeeFPjNKQ8%bcv3_gS8l%xcD;oeIF
z$#!tV)qTQCO{5jvg(}xeYe-;eN}zX7@zOSuQ<^f|yYuYW<0_rQ3-v~Fsw~;x3yjz>
z{B5%V-#NCF!2k8%Q)Ig>WmNwxWv;x`t=V(Rww8}D6e)cDLX|2+d6QP((WIulh~G?F
Z%cl%?uTa=&B(P%2@SBB7DI6w7{SUge)xZD%

delta 9173
zcmb7K3s_XwwLWuT<e|(NkT(Jc7+@eoz-YuIUS|;Gpy=lV4Jrr-gBTGI2V+v((IDzX
zX)rrCHBE1>=FyrqU!$?@HGO695pB{YsI_Xe^_4USpNR1hV{-1l_Bo@*=KF5%<~w`s
zz4qE`t-aRTYwvTI-UFuogQo7-a5-Aw#i0XX+s3S)g~*|LPu_SeH}nm=-g(2$B23L<
zzsR03T~Cn<4U+p{nMuc&v6sY;BY*Vag<E;=o*G)0eQU&zU(WkS*{E&tY(>DRp~5Oo
zW(XG!xIS9Uwu{3|@3PYYF{VCVzeY8)2IiTn(8^^y17l4`*pk3p){E9JIE{<f?=X`&
zmvzq%WFyTvthY9htuRNLD%eJIiFrQ+rGj@P3mG(PsIYWNn_FmYc~V%quvu^|Z(W?v
z+6LK2RaY-{*VZknu4`OW-B7!*NoZ|%*EQC*35#8e>sp@@7S}d4wbTh9Yi?~@*z9f~
zO15XvgbDd0g_5$;!tvE3hTm;Hgtj1vLIC~*Jb{84kKP+31feDJT0qV5Y*J8yCrG#i
zYAYV^#uo<(7`cTu6%*Qv0BMZM3w(T#V1~5`7rAg@m;f18qu(CCFofj#`MKk$FADL7
zckET1jS#E@L5M*k3QyjB1HNc>E+{)L&`$?=s#@<rK_P@a7(9b*3VzHKzzzqOn1Wbh
zNDLbmGGvD%#5aifMC&`&deJlUJrXqIcaB0yP&zIs>XVwLHAt6_E8e}zPPC~6DBfbB
zp<7vR?t4%W6z^pa_y(l#N?SpLA~o+<M?j+51xk>7bcf_mAStrhiNvQ$MDm1A@{2np
z9-T(kX)c&SQ=-#!>olB2ESCb_9YE@&wLxe}kDK2FM*Y<{T=L>OY;~@AI_ZnV^@>h1
z{XP<pRi{}>G#}|SqjefCIoL0GAjdCRA=1}CIv`o=ZSsILyB=$l?EBF#vc;k<x8*^#
zC<%PfgONb92hT1%MK&w1gHiW1NS<8tc5s%~oQ3-bq>9AX(CavFk~*J<g!d>Wc2dRt
zGa$|5<4V`43lGrp?Yq_EUt{U}iflYwxQj}`!Eg~*UfqNcf$hfJ!RyxHEBE-+r}g24
z`-ZvAeQKpX8g<{O;jsVEM>S}pi~N?$Wqq_yALaJsF4+%lUZ({pk_VA-3mXx3Y9h;2
zG}&S~C&QdFNS_5MM)u-CH1RqQb&uaQZG0S&X81B5V!?0hGbBRFy(z8P3JdaiRs`-v
z;Px)707YWK10Ft5X@ffu&B5IAzw=ujc78kCZ~(|*WO%@A16y^$ZJ6DAS$*#-z>Ymi
zzHe^T?glBjh>I>K%8J9_i5BO%06fqCoR32VgiggtjvxOv81=<sv|+fEyc&9f`6S>d
zwhuwL`VUB;&imBaf==H^Y>>Ou4c)AJX!NXEsc=D8*vS+EV`ITVi+a9P5WeU}_gyrU
zF7rDC#@uvA!KlsD?)14sv{O>@N)o%<fS~ts`cco6vfD!w#NtvG7nzt{45ZNhz)YoM
zY$`Ps?`Q5I%KEVx=%}-AX_{IOq3@d$OIb~%O_`~9zeYlVPvMzh#A=V)2_o2$TPHlw
z2HY-niH;j$BWu2ZP=__GRek4x*4M~ZL_R7OXR&u9ZDMj3I~n<exWmpSMJ0#M|B~3v
zHa~A??x?h&c02Uw{*t{CHBGGgl3j_KDE4HsvC%&f(=yrd=+w}UZ|L+|1`CQw3p$el
z`mb-W{Fr&-A8xRXG0VlP>74Mtuj_;#r*p!0(n0u_>#R05ZqOMF!QxY6#))0mSx;<s
zXck5h7GYWG>{{&H2i~$#clHNZcA3oz3orjZPnL?5?E`+7I`^8+e7%i*FU}UU81(A8
zYwV|Sv&G5Rm?=Ir`pqwh!dk|IlAEbV{m(DhnD~^aU1=~-ey?s*unT(=2(~dYl6Az7
z5glplkMU`t=dbD#64KcDcw6e%hzvEB^F;PRenL&*8EAv4<+*(ie!=n+rig#J%2p+$
z#aclvJ+o4zRaAzxp$ApWZK>?{2^r$S6!vLCig+W11r1A*jv&_5R|j-q7#cY&HS0IP
z=@Qfr^%0)%u$Wkw;UI?gLv+~a6y_OrDYZoA)JJJ09cUvW4xx?p@my*CnarL`j7iA2
zLNaKB>$^~m;s|TvwPe<lm?Vyp+3Cb-$<HJcH!?d!^r(}N8Mzs0R-!AJO-mXd)Cd}N
zXg_;7DLZ0qKQ_~9PK8Oq>XZHKNYYb@^+}w~5Suze-{gL4Tg|0OI~FIe7sc|6>{iO_
zV)F&IIn^s3if0d`O%Vgmv+t$N6(2qi<X9GEQ><^F)3h@0nKJJs#k)^kb&f5z9Z5mP
z!Cs8I5-SM#K4h5KIB2xH&q(&_u|DQXPZtMNwjn)F^r`HR>2sy4sxe#&KEF|!JtIp@
z=wmZ7&WaO0X9byAp&#OwfMt&Ah+&H|GsLuL#xfQ0SQPtvW_QF|+~uxHKZyDtKoC9H
z63smJ6frW2y=u=9hexws^ma#rHlBrLJ??Qn?wsYE=`5c)0~@$P@vg81D|O$LdT%Q|
zp9gSj$lT6$XIG_7QlI)v(~3JkeGeBQ!6}uVz|BE*0mo0TQtw`;_hiS}el-;+XO9**
z*?XzXd#<#rljs$1pp9om#anK(sy_sc;#rTTQrE9|`<0Ge4yET*vSR875EXBsPkA)b
z)_&00^`ur>w@>k&Q@W;@For8inJrK8B1rOmQnCXl9{S#WMW&mI_iswqkPnC;cFex-
zKdPr`5~!4UZ(v4|vL0JP@w#mVTCcL+0|*4Lpuj>0DIM<zE04NuvF#r!0nn|q<9$cj
zqa$tYKXHyZ@q+?qSCLIq>JDWcE$dow7|tycO-fhgUW7B92b_7C_i|~+zF@6)z$t~7
z+$<aORr}|_sc)a9?0y-?CXN@GoHRM|VR*?YW&fF6B}DRFaRxgAw*~vr2C4j$25I>D
zWe*J~QjQN(kArONmI$PQt@d<V8NZ4g`!*S7UOtZ5r-|dN2>AIy55oJcn|y-|u@!mG
zDBer{Fwgl?eI6so$G`;I)V*go8PD2Q<mrw*6?3a7cRMScRnB?NYX8A;RZ6<0?ELXd
zRNfp%r?XGl*}%&s-re6f8@mlwSuL;=1>?o9&#;s+^Thu;!&Z$MF}m;cT{-^K;joSy
zLTAMRBqSt{d12OY3gZgO5>dCU?A(~thxc<7G6d3FWe*Qdt2$U$oaX}^IvApp2s1wb
zw;uJvX;%11ny0+f`==SostdtLdG*Trz;qlV2+P&Bn>dW}eVqr}U5fWCHbO8qz;A&Q
zoKy3CWnJb@-0t7~0*F6SrFc)O8%`m@_XXOBsMLQaEA>{|tQB*V^*YCT<sn@}dycYd
zPq2t|Q06@~*IBt+@ec4EgPf2F*6L+&;|SFZE2t~cB?~*D@NMi&?^WpJfJYPONqk#J
z{SveEHQuw1B3l9Nht*JtM-uL!xRlWTg<9I70S{fgU=&<aktfN5xE`r6&~AtHqrKW8
z^tjRVamwSM1p{dfz0f#TUzAn{p@9jAQES*k(B-!z3T8eF2mKkr`l0l)Pb%KE+9!4X
zm@6^-fb{Y%sfc4d>wwC?$n}9}yWu?aajQlRQP5JVgc6A2y9vFLD|%B}`w3^6yU(xI
zIo`LQ0k30Ff?_`G6$>fEm4Rp|vv;+xJp!s7lRSP+RzHo~3waOox%cShjW#0P;b;#R
z+;N@MxgA5;2o7>x6?}L1ItFoWQZ@I=vPy`@w#$u@VU%SN7#Wp>7BE+(wQgMboIA&w
z;D21(GNEAp!7}gXaDE5cY6CXP`DZDQVI8OOAo?rL_T@uI<C=#4!Sk;+bYm5?c34Y#
zHiCvbhoClQArbJIr~7U|p6d9A5kpAsMihCjBKDW_<K7w-7e;onhArgwqF^ZNYC={;
zHT9mDT5W5@gnVQl{k{te;bLv?K|+JaaFl9?)~~STffK*OAY81S&}%4{LgzDxcj@`1
z+{3BewyB`jQw+vCuS=cZAgX|?v{j(STLEhL0LJ^fv}PB~kCrN}DbQmNUap|{Gw&R7
zC;n#dfob;NjlC#NI=uI8>;(Yh(IZXQri^b4s0$<l4$COxis?Im<FO>UlhQRAytC~d
zB`bMIUWVy&zmnFfx{4>IHP?B|-b*eaHT!0AC*owFC3u!SO`b$E>>KE=hskZUfEke1
zz6!VDEU~px0&4U^9K{4rLp=8&?NBlbJ2x14q|1B@fTtrh<b<YG%;CxBUsBTf%)2q<
z`P2E;4c*@v_i!F%E%Nm}q5lGf{;E9`GDh?}>CMM}3J=mV*7qcIHUz)}C#7c_2+>)C
z+Mpv_5bdeHNghs&J>XwsJJ#5j=r_?_{L_Rt|0{kI;V*)V+IXDg>Y46Nrc}H~@FsE%
z=Cj`u!HDd9pLzo+)uUk3c`vCK=&pxdj@WD&Vk@P41hzhoJK$pPEWdp)cD)_zrcJNj
z{N(>Ba4iHH@*cwr4)MkZh<{(;y88t3cP#G>EKdHH#SObFD;#TEA870JwGaFDauCMx
zi0yVOy$Mw9)=!()Aw*3^#hms?_F2}@=sm}ft2(|m{a}={tN5EToIG6*c7`xZwk#%n
z!t$~c$|r>&XWa(bT)aIx(Sd#Fa?sZZ+_!x4cv;ujuaIM%vr4^xbIvUFUUQZ^aiD?a
z1)yh?Wzkz>u=@6q9nU?o#Js~1>JXoI1be2`?tXijOHNOh-Ho+w`AL^t+vbv+TikMO
zQ%ke!5jnl7ArzrG6CV#3$ctMrSX)=u-d5{&$xST_mbTW`x#YT*#rSww>t5K>yfo+k
z<^^^_^vqxXZu_E$*^VgQ5PgKi0pZ|4=pVOPd+w!Yw<P%V@d)4CpTC;@%AnDFeC6ul
zPMyx-)8U<ZrN6~m9kKkV!RNbie|Ta+$n<HU4z|gW6Au*~Ir740Ien>YPp{7s(i?JQ
zR6-Stc0N*4+U$0<Evs$9s+!%k3!7=(?namFu5DZ3a?2QL*A3Lr(%xLp88x&w*U{p^
zbeXGdsjIrKzJXO2#>u+eTjO-4`rpGgL-|+-+Q%NYvv8P)z8}-&Ncc$8R-!Fqd@Eoc
zpbyZ|t7&v;68@@bjez}tM+tve)9gqW+W;Q}^!*K~2{7-7rZoZ99M!aSfZd3d*8zJG
zdp^LaNMWA>Zab-IkpTiuN`&-Cd@QN~EC=ibYy`C8ezXFx8?c-3xcI#VDBx}4CgA3O
zYFbes`go1%2DIXBYbT(L_p?)gn=fcu4!-L30y+TmE<!=TZopQ+e!z8R1sYw_w9P;W
zmo@DKpaburR(#i#@piouup96#KmqUgM*&*_Z&M#%4GQqVvK4S;5E~eu5#N1K({eE)
z6a{hb+t8#B)64(NMigauTK@vhG;9>S1_*9Ejf5u~8)X>a20Ye7nl^!O#TM%tvBV;~
zEmo%`SjYkDEj$kR>{nzGK0d(`+YwM|sR_sr2(`qH=cFY#I`l8rC+M3qF`Q2r6uhbm
z{dk%Vu}ekq(KYbvBEqB4kPLl(|5uh&9541B28d=;ivyEBg_lQ>rivv`{4CTGJbths
zqw5AP_$X^DP8I8qu@{Sni!YpD?-wVDdrq)_7AMUK{JW;PKu3ivFbok;0~{HL3gPmA
z+W_2C9HVPd1l&#F9wnIpi6)CY8|#O7JT({_c9PATFi{jvv!4T$pJCSl8qc!Ak_>U`
zr|gN6BoC3k03#J3Zl1yvVS>f_oLFdy{b4|%MeYhLwAeS8otC^cgPfLvj-Vn-w<(}P
zv=jj4wAd$DWRMnGtc4_()I0@VO`mI8CfZcUiC+M>4!Fk&(PEM7Kuq`{fe7&4?-;Nl
zKw}ZOGiicikwwK4OcNdhz6N;Ch-^H;5?la$Iq+M7-$wLUzy`61yj^1195~OiHE^P3
zTVNZY6a6`Gei7H#IatV}4Mw|T7<PxlSUsP-CctTtI|54(FDp!;78(4UeAcAFO^~q%
zbL{64@Wz}jve<@zB8$8RQtTb(iI&agfXyb$R&%jso4E-P!p&k;sKpK<r{88y&tS1M
z)S#L;_`90LgY6JF5<71bwq)a-6_!|H*I<PXQMq?TL|X&e6aNOS4BG4sAeOPkYl$)6
zC=+2<XDHtvFcFz|8)%v?v7eTV2rWXF{$t_0#O9Y4hK$5US25eXBR_OSfR+AjgW|(C
z7b>QWTZXc-$u@CQ0c)E)FQp$7sZe&HvY$353Y}q8GT4>L)hX)`NmP(Kc-e0>47A{5
zOQwu7DQx?cWK;Bx15=X4I6q@!1_uH7Z&EZ%nHCigA}D|{rv1D(SXjrVO?!@XIx}s!
zNoL{GV@$0qYx;90E88+XihVSFV(C<41%8GJXo!gKmG|L?8+iYU{q$relMpXhS?!Em
zQHWvB&zL95nQUOjPa*N8@=VhwtfxF#?1mEM8Kz2hr939T8M`<FLR#@LjtZrFD&I5k
z@@T!!U4`ZkX_omljLP_j;?TH6(8~fot<~ss{`x(li#uR&+4iuG=aGcwr)<^C7}I&K
z!B>77S}zLEWq^+mUJMsDQ8e6BBH5N%F&@fHq<23hW)!;3QhCI{`zR<;<{0=%emq$Z
z#h)sw47@c@$Mao|<_ZJ9jh`L3?&!n9%Q`+=n2I%1NjC(%X-tswbv)nMX#Sp=9#0PM
zHKLUMSxMz2EyDB-)sw?vLR!kUG<}RGRx}fUw|<9d?FN0`13Eo@7NAl9Ja<ZB30mP|
zoDf{Y2kunl+z4`;d4a#>d75Y;UvTJTNHV;5iGmEAs`vfJ!i#>qO~>b8-BjK%@Q(Xj
z{c8h1)i9*rz`Z&?T96I@x^zA~nR6vQl)_1cmGSzJ-_&O4NO=F&=jSm<{)xwm4G>$S
z)BEcfgTB@9G0)f-OfmTT|Nn1jD5UH39zn*xf~c4b3y&}u^&3ItU&{mopO>c7`;W_M
zz>`AO8XfPis|@-Yn4QY6pa_-u9Iv<Wow0e&;8=swLH7$3ep&*qXG><s#Qnq=?H<hs
zg>^a*w$iBSw`>b&-ZyAY#p;9p3&mk}05q40a<U0Ik{jCawcizk4Sa8kKDXH5A7kKa
z9?}u~L_#wYc>0~e(8Tcug=}z~X3U&Uc%m?A6#Jp`{hVTUyfbHtDdF2cE^s*fKP$kG
z379rlj#yf|xV6dEM)&FJ`r76Nt~QobS;T%>J!pu5Y-yk0<YKjz`NZp6zh}71Z6w=R
zIV{}I3zzunF?rQ>3)x$hlOyiAu~*mY%x+fBV)LtF*oLZc5&Y*8cR7t1K8k%%m6Pq~
zMS9I|SuWHzEpW|mtF3NrYpHLqb64ZdsJW#Lzo1ykQszCJtdkSAy?J2+1_Th+xt1=i
zcD1#&w6V5%6&62BQV+TvpO+NqTE5WDGOHI!_bLf-Y-9CU7FsoU$HnUTrX5pjWs}Lq
zmdqb#va{{;bHn{7l^^2B12)<?;{0$p%dg7~R}8Elu4GGq^&f|R=wsXKazp%w=DHny
zbtNWoR0bR28XoG;*Es8*a81C;v)z>&e*aHp%GecGZba*NJuLiC8MI?WL$fK|zhA6V
l_3>b~q&AuLE~qh`V)>1uL*%(Sei>VW-itMQFCNQE{$HXb(~kfE

diff --git a/scratchpad/scratchpad.c b/scratchpad/scratchpad.c
index da52936..b173ad3 100644
--- a/scratchpad/scratchpad.c
+++ b/scratchpad/scratchpad.c
@@ -1,14 +1,10 @@
 #include "../squiggle.h"
-// #include "../squiggle_more.h"
+#include "../squiggle_more.h"
 #include <math.h>
 #include <stdint.h>
 #include <stdio.h>
 #include <stdlib.h>
 
-double sample_loguniform(double a, double b, uint64_t* seed){
-    return exp(sample_uniform(log(a), log(b), seed));
-}
-
 int main()
 {
     // Replicate <https://arxiv.org/pdf/1806.02404.pdf>, and in particular the red line in page 11.
@@ -18,140 +14,32 @@ int main()
     uint64_t* seed = malloc(sizeof(uint64_t));
     *seed = UINT64_MAX/64; // xorshift can't start with a seed of 0
     
-    double sample_fermi_naive(uint64_t* seed){
-        double rate_of_star_formation = sample_loguniform(1,100, seed);
-        double fraction_of_stars_with_planets = sample_loguniform(0.1, 1, seed);
-        double number_of_habitable_planets_per_star_system = sample_loguniform(0.1, 1, seed);
-        double rate_of_life_formation_in_habitable_planets = sample_lognormal(1, 50, seed);
-        double fraction_of_habitable_planets_in_which_any_life_appears = -expm1(-rate_of_life_formation_in_habitable_planets);
-        // double fraction_of_habitable_planets_in_which_any_life_appears = 1-exp(-rate_of_life_formation_in_habitable_planets);
-        // but with more precision
-        double fraction_of_planets_with_life_in_which_intelligent_life_appears = sample_loguniform(0.001, 1, seed);
-        double fraction_of_intelligent_planets_which_are_detectable_as_such = sample_loguniform(0.01, 1, seed);
-        double longevity_of_detectable_civilizations = sample_loguniform(100, 10000000000, seed);
-
-        // printf(" rate_of_star_formation = %lf\n", rate_of_star_formation);
-        // printf(" fraction_of_stars_with_planets = %lf\n", fraction_of_stars_with_planets);
-        // printf(" number_of_habitable_planets_per_star_system = %lf\n", number_of_habitable_planets_per_star_system);
-        // printf(" rate_of_life_formation_in_habitable_planets = %.16lf\n", rate_of_life_formation_in_habitable_planets);
-        // printf(" fraction_of_habitable_planets_in_which_any_life_appears = %lf\n", fraction_of_habitable_planets_in_which_any_life_appears);
-        // printf(" fraction_of_planets_with_life_in_which_intelligent_life_appears = %lf\n", fraction_of_planets_with_life_in_which_intelligent_life_appears);
-        // printf(" fraction_of_intelligent_planets_which_are_detectable_as_such = %lf\n", fraction_of_intelligent_planets_which_are_detectable_as_such);
-        // printf(" longevity_of_detectable_civilizations = %lf\n", longevity_of_detectable_civilizations);
-
-        // Expected number of civilizations in the Milky way;
-        // see footnote 3 (p. 5)
-        double n = rate_of_star_formation *
-          fraction_of_stars_with_planets *
-          number_of_habitable_planets_per_star_system * 
-          fraction_of_habitable_planets_in_which_any_life_appears * 
-          fraction_of_planets_with_life_in_which_intelligent_life_appears * 
-          fraction_of_intelligent_planets_which_are_detectable_as_such *
-          longevity_of_detectable_civilizations;
-
-        return n;
-    }
-
-    double sample_fermi_paradox_naive(uint64_t* seed){
-        double n = sample_fermi_naive(seed);
-        return ((n > 1) ? 1 : 0);
-    }
-
-    double n = 1000000;
-    double naive_fermi_proportion = 0;
-    for(int i=0; i<n; i++){
-        double result = sample_fermi_paradox_naive(seed);
-        // printf("result: %lf\n", result);
-        naive_fermi_proportion+=result;
-    }
-    printf("Naïve %% that we are not alone: %lf\n", naive_fermi_proportion/n);
-
-
-    // Thinking in log space
-    double sample_fermi_logspace(uint64_t* seed){
-        double log_rate_of_star_formation = sample_uniform(log(1), log(100), seed);
-        double log_fraction_of_stars_with_planets = sample_uniform(log(0.1), log(1), seed);
-        double log_number_of_habitable_planets_per_star_system = sample_uniform(log(0.1), log(1), seed);
-        double log_fraction_of_planets_with_life_in_which_intelligent_life_appears = sample_uniform(log(0.001), log(1), seed);
-        double log_fraction_of_intelligent_planets_which_are_detectable_as_such = sample_uniform(log(0.01), log(1), seed);
-        double log_longevity_of_detectable_civilizations = sample_uniform(log(100), log(10000000000), seed);
-
-        // printf(" log_rate_of_star_formation = %lf\n", log_rate_of_star_formation);
-        // printf(" log_fraction_of_stars_with_planets = %lf\n", log_fraction_of_stars_with_planets);
-        // printf(" log_number_of_habitable_planets_per_star_system = %lf\n", log_number_of_habitable_planets_per_star_system);
-        // printf(" log_fraction_of_planets_with_life_in_which_intelligent_life_appears = %lf\n", log_fraction_of_planets_with_life_in_which_intelligent_life_appears);
-        // printf(" log_fraction_of_intelligent_planets_which_are_detectable_as_such = %lf\n", log_fraction_of_intelligent_planets_which_are_detectable_as_such);
-        // printf(" log_longevity_of_detectable_civilizations = %lf\n", log_longevity_of_detectable_civilizations);
-
-        double log_n1 = 
-            log_rate_of_star_formation +
-            log_fraction_of_stars_with_planets +
-            log_number_of_habitable_planets_per_star_system + 
-            log_fraction_of_planets_with_life_in_which_intelligent_life_appears + 
-            log_fraction_of_intelligent_planets_which_are_detectable_as_such +
-            log_longevity_of_detectable_civilizations;
-        // printf("first part of calculation: %lf\n", log_n1);
-
-        /* Consider fraction_of_habitable_planets_in_which_any_life_appears separately.
-        Imprecisely, we could do:
-
-        double rate_of_life_formation_in_habitable_planets = sample_lognormal(1, 50, seed);
-        double fraction_of_habitable_planets_in_which_any_life_appears = 1- exp(-rate_of_life_formation_in_habitable_planets);
-        double log_fraction_of_habitable_planets_in_which_any_life_appears = log(1-fraction_of_habitable_planets_in_which_any_life_appears);
-        double n = exp(log_n1) * fraction_of_habitable_planets_in_which_any_life_appears;
-        // or: 
-        double n2 = exp(log_n1 + log(fraction_of_habitable_planets_in_which_any_life_appears))
-
-        However, we lose all precision here.
-
-        Now, say
-        a = underlying normal 
-        b = rate_of_life_formation_in_habitable_planets = exp(underlying normal)
-        c = 1 - exp(-b) = fraction_of_habitable_planets_in_which_any_life_appears
-        d = log(c)
-
-        Now, is there some way we can d more efficiently/precisely?
-        Turns out there is!
-
-        Looking at the Taylor expansion for c = 1 - exp(-b), it's b - b^2/2 + b^3/6 - x^b/24, etc.
-        // https://www.wolframalpha.com/input?i=1-exp%28-x%29
-        When b ~ 0 (as is often the case), this is close to b.
-
-        But now, if b ~ 0
-        c ~ b
-        and d = log(c) ~ log(b) = log(exp(a)) = a
-        */
-        double log_rate_of_life_formation_in_habitable_planets = sample_normal(1, 50, seed);
-        // printf("log_rate_of_life_formation_in_habitable_planets: %lf\n", log_rate_of_life_formation_in_habitable_planets);
-
-        double log_fraction_of_habitable_planets_in_which_any_life_appears;
-        if(log_rate_of_life_formation_in_habitable_planets < -32){
-            log_fraction_of_habitable_planets_in_which_any_life_appears = log_rate_of_life_formation_in_habitable_planets;
-        } else{
-            double rate_of_life_formation_in_habitable_planets = exp(log_rate_of_life_formation_in_habitable_planets);
-            double fraction_of_habitable_planets_in_which_any_life_appears = -expm1(-rate_of_life_formation_in_habitable_planets);
-            log_fraction_of_habitable_planets_in_which_any_life_appears = log(fraction_of_habitable_planets_in_which_any_life_appears);
+    int n_samples = 100*MILLION;
+    int p_sixteenth = 0;
+    int p_eighth = 0;
+    int p_quarter = 0;
+    int p_half = 0;
+    double sample;
+    for(int i=0; i<n_samples; i++){
+        sample = sample_unit_uniform(seed);
+        // printf("%lf\n", sample);
+        if (sample < 1.0/16.0){
+            p_sixteenth++;
+            p_eighth++;
+            p_quarter++;
+            p_half++;
+        } else if(sample < 0.125){
+            p_eighth++;
+            p_quarter++;
+            p_half++;
+        } else if(sample < 0.25){
+            p_quarter++;
+            p_half++;
+        } else if(sample < 0.5){
+            p_half++;
+        }else{
+            // printf("Sample > 0.5\n");
         }
-        // printf(" log_fraction_of_habitable_planets_in_which_any_life_appears: %lf\n", log_fraction_of_habitable_planets_in_which_any_life_appears);
-
-        double log_n = log_n1 + log_fraction_of_habitable_planets_in_which_any_life_appears;
-
-        return log_n;
     }
-
-    double sample_fermi_paradox_logspace(uint64_t* seed){
-        double n = sample_fermi_logspace(seed);
-        return ((n > 0) ? 1 : 0);
-    }
-
-    double logspace_fermi_proportion = 0;
-    for(int i=0; i<n; i++){
-        double result = sample_fermi_paradox_logspace(seed);
-        // printf("result: %lf\n", result);
-        logspace_fermi_proportion+=result;
-    }
-    printf("Using more accurate logspace computations, %% that we are not alone: %lf\n", logspace_fermi_proportion/n);
-    double result2;
-
-    free(seed);
+    printf("p_16th: %lf; p_eighth; %lf; p_quarter: %lf; p_half: %lf", ((double)p_sixteenth)/n_samples, (double)p_eighth/n_samples, (double)p_quarter/n_samples, (double)p_half/n_samples);
 }