From c7d9f87ad7dada50a314ee6aa159935dded998b7 Mon Sep 17 00:00:00 2001
From: NunoSempere <nuno.sempere@protonmail.com>
Date: Sun, 16 Jul 2023 11:08:59 +0200
Subject: [PATCH] add sampler, add normal cdf

---
 scratchpad/scratchpad.c | 42 ++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 41 insertions(+), 1 deletion(-)

diff --git a/scratchpad/scratchpad.c b/scratchpad/scratchpad.c
index 9ed24e9..120814c 100644
--- a/scratchpad/scratchpad.c
+++ b/scratchpad/scratchpad.c
@@ -3,6 +3,7 @@
 #include <stdio.h>
 #include <float.h> // FLT_MAX, FLT_MIN
 #include <limits.h> // INT_MAX
+#include <math.h>
 #define VERBOSE 1
 
 // to do: reuse more informative printing from build-your-own-lisp
@@ -42,6 +43,12 @@ float cdf_squared_0_1(float x){
 	}
 }
 
+float cdf_normal_0_1(float x){
+	float mean = 0;
+	float std = 1;
+	return 0.5 * ( 1 + erf((x-mean)/(std * sqrt(2)) ));
+}
+
 // Inverse cdf
 struct box inverse_cdf(float cdf(float), float p){
 	// given a cdf: [-Inf, Inf] => [0,1]
@@ -124,8 +131,37 @@ struct box inverse_cdf(float cdf(float), float p){
 
 }
 
+// Get random number between 0 and 1
+uint32_t xorshift32
+(uint32_t* seed)
+{
+	// Algorithm "xor" from p. 4 of Marsaglia, "Xorshift RNGs"
+  // See <https://stackoverflow.com/questions/53886131/how-does-xorshift32-works>
+  // https://en.wikipedia.org/wiki/Xorshift
+	// Also some drama: <https://www.pcg-random.org/posts/on-vignas-pcg-critique.html>, <https://prng.di.unimi.it/>
+
+	uint32_t x = *seed;
+	x ^= x << 13;
+	x ^= x >> 17;
+	x ^= x << 5;
+	return *seed = x;
+}
+
+// Distribution & sampling functions
+
+float rand_0_to_1(uint32_t* seed){
+	return ((float) xorshift32(seed)) / ((float) UINT32_MAX);
+}
+
 // sampler based on inverse cdf
-// to-do: integrals
+struct box sampler(float cdf(float), uint32_t* seed){
+	struct box result;
+	float p = rand_0_to_1(seed);
+	result = inverse_cdf(cdf, p);
+	return result;
+}
+
+// to-do: integrals => beta distribution
 
 // main with an example
 int main(){
@@ -147,6 +183,10 @@ int main(){
 	}else{
 		printf("Inverse of the cdf at %f is: %f\n", 0.5, result2.content);
 	}
+
+	// set randomness seed
+	uint32_t* seed = malloc(sizeof(uint32_t));
+	*seed = 1000; // xorshift can't start with 0
 	
 	return 0;
 }