move to xorshift64. Better precision.

README.md

@@ -128,7 +128,7 @@ In squiggle.c, this ambiguity doesn't exist, at the cost of much greater overhea
 
 int main(){
     // set randomness seed
-    uint32_t* seed = malloc(sizeof(uint32_t));
+    uint64_t* seed = malloc(sizeof(uint64_t));
     *seed = 1000; // xorshift can't start with a seed of 0
 
     float a = sample_to(1, 10, seed);
@@ -153,7 +153,7 @@ vs
 #include <stdlib.h>
 #include <stdio.h>
 
-float draw_xyz(uint32_t* seed){
+float draw_xyz(uint64_t* seed){
     // function could also be placed inside main with gcc nested functions extension.
     return sample_to(1, 20, seed);
 }
@@ -161,7 +161,7 @@ float draw_xyz(uint32_t* seed){
 
 int main(){
     // set randomness seed
-    uint32_t* seed = malloc(sizeof(uint32_t));
+    uint64_t* seed = malloc(sizeof(uint64_t));
     *seed = 1000; // xorshift can't start with a seed of 0
 
     float a = draw_xyz(seed);

examples/01_one_sample/example.c

@@ -4,29 +4,29 @@
 #include <stdio.h>
 
 // Estimate functions
-float sample_0(uint32_t* seed)
+float sample_0(uint64_t* seed)
 {
     return 0;
 }
 
-float sample_1(uint32_t* seed)
+float sample_1(uint64_t* seed)
 {
     return 1;
 }
 
-float sample_few(uint32_t* seed)
+float sample_few(uint64_t* seed)
 {
     return sample_to(1, 3, seed);
 }
 
-float sample_many(uint32_t* seed)
+float sample_many(uint64_t* seed)
 {
     return sample_to(2, 10, seed);
 }
 
 int main(){
     // set randomness seed
-		uint32_t* seed = malloc(sizeof(uint32_t));
+		uint64_t* seed = malloc(sizeof(uint64_t));
 		*seed = 1000; // xorshift can't start with 0
 
     float p_a = 0.8;
@@ -35,7 +35,7 @@ int main(){
 
     int n_dists = 4;
     float weights[] = { 1 - p_c, p_c / 2, p_c / 4, p_c / 4 };
-    float (*samplers[])(uint32_t*) = { sample_0, sample_1, sample_few, sample_many };
+    float (*samplers[])(uint64_t*) = { sample_0, sample_1, sample_few, sample_many };
 
     float result_one = sample_mixture(samplers, weights, n_dists, seed);
 		printf("result_one: %f\n", result_one);

examples/02_many_samples/example.c

@@ -4,29 +4,29 @@
 #include "../../squiggle.h"
 
 // Estimate functions
-float sample_0(uint32_t* seed)
+float sample_0(uint64_t* seed)
 {
     return 0;
 }
 
-float sample_1(uint32_t* seed)
+float sample_1(uint64_t* seed)
 {
     return 1;
 }
 
-float sample_few(uint32_t* seed)
+float sample_few(uint64_t* seed)
 {
     return sample_to(1, 3, seed);
 }
 
-float sample_many(uint32_t* seed)
+float sample_many(uint64_t* seed)
 {
     return sample_to(2, 10, seed);
 }
 
 int main(){
     // set randomness seed
-		uint32_t* seed = malloc(sizeof(uint32_t));
+		uint64_t* seed = malloc(sizeof(uint64_t));
 		*seed = 1000; // xorshift can't start with 0
 
     float p_a = 0.8;
@@ -35,7 +35,7 @@ int main(){
 
     int n_dists = 4;
     float weights[] = { 1 - p_c, p_c / 2, p_c / 4, p_c / 4 };
-    float (*samplers[])(uint32_t*) = { sample_0, sample_1, sample_few, sample_many };
+    float (*samplers[])(uint64_t*) = { sample_0, sample_1, sample_few, sample_many };
 
 		int n_samples = 1000000;
 		float* result_many = (float *) malloc(n_samples * sizeof(float));

examples/03_gcc_nested_function/example.c

@@ -5,7 +5,7 @@
 
 int main(){
     // set randomness seed
-		uint32_t* seed = malloc(sizeof(uint32_t));
+		uint64_t* seed = malloc(sizeof(uint64_t));
 		*seed = 1000; // xorshift can't start with 0
 
     float p_a = 0.8;
@@ -14,12 +14,12 @@ int main(){
 
     int n_dists = 4;
 		
-		float sample_0(uint32_t* seed){ return 0; }
-		float sample_1(uint32_t* seed) { return 1; }
-		float sample_few(uint32_t* seed){ return sample_to(1, 3, seed); }
-		float sample_many(uint32_t* seed){ return sample_to(2, 10, seed); } 
+		float sample_0(uint64_t* seed){ return 0; }
+		float sample_1(uint64_t* seed) { return 1; }
+		float sample_few(uint64_t* seed){ return sample_to(1, 3, seed); }
+		float sample_many(uint64_t* seed){ return sample_to(2, 10, seed); } 
 		
-    float (*samplers[])(uint32_t*) = { sample_0, sample_1, sample_few, sample_many };
+    float (*samplers[])(uint64_t*) = { sample_0, sample_1, sample_few, sample_many };
 		float weights[] = { 1 - p_c, p_c / 2, p_c / 4, p_c / 4 };
 
 		int n_samples = 1000000;

examples/04_sample_from_cdf_simple/example.c

@@ -49,7 +49,7 @@ void test_inverse_cdf_float(char* cdf_name, float cdf_float(float))
     }
 }
 
-void test_and_time_sampler_float(char* cdf_name, float cdf_float(float), uint32_t* seed)
+void test_and_time_sampler_float(char* cdf_name, float cdf_float(float), uint64_t* seed)
 {
     printf("\nGetting some samples from %s:\n", cdf_name);
     clock_t begin = clock();
@@ -75,7 +75,7 @@ int main()
 
     // Testing samplers
     // set randomness seed
-    uint32_t* seed = malloc(sizeof(uint32_t));
+    uint64_t* seed = malloc(sizeof(uint64_t));
     *seed = 1000; // xorshift can't start with 0
 
     // Test float sampler

examples/05_sample_from_cdf_beta/example.c

@@ -131,7 +131,7 @@ void test_inverse_cdf_box(char* cdf_name, struct box cdf_box(float))
     }
 }
 
-void test_and_time_sampler_box(char* cdf_name, struct box cdf_box(float), uint32_t* seed)
+void test_and_time_sampler_box(char* cdf_name, struct box cdf_box(float), uint64_t* seed)
 {
     printf("\nGetting some samples from %s:\n", cdf_name);
     clock_t begin = clock();
@@ -154,7 +154,7 @@ int main()
     test_inverse_cdf_box("cdf_beta", cdf_beta);
 
     // Test box sampler
-    uint32_t* seed = malloc(sizeof(uint32_t));
+    uint64_t* seed = malloc(sizeof(uint64_t));
     *seed = 1000; // xorshift can't start with 0
     test_and_time_sampler_box("cdf_beta", cdf_beta, seed);
     // Ok, this is slower than python!!

examples/06_gamma_beta/example.c

@@ -8,7 +8,7 @@
 int main()
 {
     // set randomness seed
-    uint32_t* seed = malloc(sizeof(uint32_t));
+    uint64_t* seed = malloc(sizeof(uint64_t));
     *seed = 1000; // xorshift can't start with 0
 
     int n =  1000 * 1000;

squiggle.c

@@ -14,14 +14,14 @@
 const float PI = 3.14159265358979323846; // M_PI in gcc gnu99
 
 // Pseudo Random number generator
-uint32_t xorshift32(uint32_t* seed)
+uint64_t xorshift32(uint32_t* seed)
 {
     // Algorithm "xor" from p. 4 of Marsaglia, "Xorshift RNGs"
-    // See <https://stackoverflow.com/questions/53886131/how-does-xorshift32-works>
+    // See <https://stackoverflow.com/questions/53886131/how-does-xorshift64-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;
+    uint64_t x = *seed;
     x ^= x << 13;
     x ^= x >> 17;
     x ^= x << 5;
@@ -31,7 +31,7 @@ uint32_t xorshift32(uint32_t* seed)
 uint64_t xorshift64(uint64_t* seed)
 {
     // Algorithm "xor" from p. 4 of Marsaglia, "Xorshift RNGs"
-    // See <https://stackoverflow.com/questions/53886131/how-does-xorshift32-works>
+    // See <https://stackoverflow.com/questions/53886131/how-does-xorshift64-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/>
 
@@ -44,13 +44,13 @@ uint64_t xorshift64(uint64_t* seed)
 
 // Distribution & sampling functions
 // Unit distributions
-float sample_unit_uniform(uint32_t* seed)
+float sample_unit_uniform(uint64_t* seed)
 {
     // samples uniform from [0,1] interval.
-    return ((float)xorshift32(seed)) / ((float)UINT32_MAX);
+    return ((float)xorshift64(seed)) / ((float)UINT64_MAX);
 }
 
-float sample_unit_normal(uint32_t* seed)
+float sample_unit_normal(uint64_t* seed)
 {
     // See: <https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform>
     float u1 = sample_unit_uniform(seed);
@@ -60,22 +60,22 @@ float sample_unit_normal(uint32_t* seed)
 }
 
 // Composite distributions
-float sample_uniform(float start, float end, uint32_t* seed)
+float sample_uniform(float start, float end, uint64_t* seed)
 {
     return sample_unit_uniform(seed) * (end - start) + start;
 }
 
-float sample_normal(float mean, float sigma, uint32_t* seed)
+float sample_normal(float mean, float sigma, uint64_t* seed)
 {
     return (mean + sigma * sample_unit_normal(seed));
 }
 
-float sample_lognormal(float logmean, float logsigma, uint32_t* seed)
+float sample_lognormal(float logmean, float logsigma, uint64_t* seed)
 {
     return expf(sample_normal(logmean, logsigma, seed));
 }
 
-float sample_to(float low, float high, uint32_t* seed)
+float sample_to(float low, float high, uint64_t* seed)
 {
     // Given a (positive) 90% confidence interval,
     // returns a sample from a lognormal
@@ -88,7 +88,7 @@ float sample_to(float low, float high, uint32_t* seed)
     return sample_lognormal(logmean, logsigma, seed);
 }
 
-float sample_gamma(float alpha, uint32_t* seed)
+float sample_gamma(float alpha, uint64_t* seed)
 {
 
     // A Simple Method for Generating Gamma Variables, Marsaglia and Wan Tsang, 2001
@@ -125,7 +125,7 @@ float sample_gamma(float alpha, uint32_t* seed)
     }
 }
 
-float sample_beta(float a, float b, uint32_t* seed)
+float sample_beta(float a, float b, uint64_t* seed)
 {
     float gamma_a = sample_gamma(a, seed);
     float gamma_b = sample_gamma(b, seed);
@@ -168,7 +168,7 @@ float array_std(float* array, int length)
 }
 
 // Mixture function
-float sample_mixture(float (*samplers[])(uint32_t*), float* weights, int n_dists, uint32_t* seed)
+float sample_mixture(float (*samplers[])(uint64_t*), float* weights, int n_dists, uint64_t* seed)
 {
     // You can see a simpler version of this function in the git history
     // or in C-02-better-algorithm-one-thread/
@@ -364,13 +364,13 @@ struct box inverse_cdf_box(struct box cdf_box(float), float p)
 }
 
 // Sampler based on inverse cdf and randomness function
-struct box sampler_cdf_box(struct box cdf(float), uint32_t* seed)
+struct box sampler_cdf_box(struct box cdf(float), uint64_t* seed)
 {
     float p = sample_unit_uniform(seed);
     struct box result = inverse_cdf_box(cdf, p);
     return result;
 }
-struct box sampler_cdf_float(float cdf(float), uint32_t* seed)
+struct box sampler_cdf_float(float cdf(float), uint64_t* seed)
 {
     float p = sample_unit_uniform(seed);
     struct box result = inverse_cdf_float(cdf, p);
@@ -378,7 +378,7 @@ struct box sampler_cdf_float(float cdf(float), uint32_t* seed)
 }
 
 /* Could also define other variations, e.g.,
-float sampler_danger(struct box cdf(float), uint32_t* seed)
+float sampler_danger(struct box cdf(float), uint64_t* seed)
 {
     float p = sample_unit_uniform(seed);
     struct box result = inverse_cdf_box(cdf, p);

squiggle.h

@@ -1,24 +1,24 @@
 #ifndef SQUIGGLEC
 #define SQUIGGLEC
 
-// uint32_t header
+// uint64_t header
 #include <stdint.h>
 
 // Pseudo Random number generator
-uint32_t xorshift32(uint32_t* seed);
+uint64_t xorshift64(uint64_t* seed);
 
 // Basic distribution sampling functions
-float sample_unit_uniform(uint32_t* seed);
-float sample_unit_normal(uint32_t* seed);
+float sample_unit_uniform(uint64_t* seed);
+float sample_unit_normal(uint64_t* seed);
 
 // Composite distribution sampling functions
-float sample_uniform(float start, float end, uint32_t* seed);
-float sample_normal(float mean, float sigma, uint32_t* seed);
-float sample_lognormal(float logmean, float logsigma, uint32_t* seed);
-float sample_to(float low, float high, uint32_t* seed);
+float sample_uniform(float start, float end, uint64_t* seed);
+float sample_normal(float mean, float sigma, uint64_t* seed);
+float sample_lognormal(float logmean, float logsigma, uint64_t* seed);
+float sample_to(float low, float high, uint64_t* seed);
 
-float sample_gamma(float alpha, uint32_t* seed);
-float sample_beta(float a, float b, uint32_t* seed);
+float sample_gamma(float alpha, uint64_t* seed);
+float sample_beta(float a, float b, uint64_t* seed);
 
 // Array helpers
 float array_sum(float* array, int length);
@@ -27,7 +27,7 @@ float array_mean(float* array, int length);
 float array_std(float* array, int length);
 
 // Mixture function
-float sample_mixture(float (*samplers[])(uint32_t*), float* weights, int n_dists, uint32_t* seed);
+float sample_mixture(float (*samplers[])(uint64_t*), float* weights, int n_dists, uint64_t* seed);
 
 // Box
 struct box {
@@ -47,7 +47,7 @@ struct box inverse_cdf_float(float cdf(float), float p);
 struct box inverse_cdf_box(struct box cdf_box(float), float p);
 
 // Samplers from cdf
-struct box sampler_cdf_float(float cdf(float), uint32_t* seed);
-struct box sampler_cdf_box(struct box cdf(float), uint32_t* seed);
+struct box sampler_cdf_float(float cdf(float), uint64_t* seed);
+struct box sampler_cdf_box(struct box cdf(float), uint64_t* seed);
 
 #endif

test/test.c

@@ -4,9 +4,9 @@
 #include <stdlib.h>
 #include <stdio.h>
 
-#define N 10 * 1000 * 1000
+#define N  1000 * 1000
 
-void test_unit_uniform(uint32_t* seed){
+void test_unit_uniform(uint64_t* seed){
 	float* unit_uniform_array = malloc(sizeof(float) * N);
 	
 	for(int i=0; i<N; i++){
@@ -40,7 +40,7 @@ void test_unit_uniform(uint32_t* seed){
 
 }
 
-void test_uniform(float start, float end, uint32_t* seed){
+void test_uniform(float start, float end, uint64_t* seed){
 	float* uniform_array = malloc(sizeof(float) * N);
 	
 	for(int i=0; i<N; i++){
@@ -76,7 +76,7 @@ void test_uniform(float start, float end, uint32_t* seed){
 
 int main(){
     // set randomness seed
-    uint32_t* seed = malloc(sizeof(uint32_t));
+    uint64_t* seed = malloc(sizeof(uint64_t));
     *seed = 1000; // xorshift can't start with a seed of 0
    
 		test_unit_uniform(seed);