From 6e228dcc6bfdab03cf0259ac1f7b45431960ccaf Mon Sep 17 00:00:00 2001
From: NunoSempere <nuno.sempere@protonmail.com>
Date: Sun, 23 Jul 2023 13:02:56 +0200
Subject: [PATCH] replace all floats (32 bits) with doubles (64 bits)

to fix bug after switching xorshift32 => xorshift64
---
 README.md                                    |  16 +--
 examples/01_one_sample/example               | Bin 22288 -> 22290 bytes
 examples/01_one_sample/example.c             |  20 ++--
 examples/02_many_samples/example             | Bin 22328 -> 22330 bytes
 examples/02_many_samples/example.c           |  20 ++--
 examples/03_gcc_nested_function/example      | Bin 22344 -> 22346 bytes
 examples/03_gcc_nested_function/example.c    |  20 ++--
 examples/04_sample_from_cdf_simple/example   | Bin 22432 -> 22435 bytes
 examples/04_sample_from_cdf_simple/example.c |  40 +++----
 examples/05_sample_from_cdf_beta/example     | Bin 26552 -> 26554 bytes
 examples/05_sample_from_cdf_beta/example.c   |  24 ++--
 examples/06_gamma_beta/example               | Bin 22160 -> 22162 bytes
 examples/06_gamma_beta/example.c             |  10 +-
 squiggle.c                                   | 116 +++++++++----------
 squiggle.h                                   |  36 +++---
 test/test                                    | Bin 22312 -> 22314 bytes
 test/test.c                                  |  45 +++----
 17 files changed, 175 insertions(+), 172 deletions(-)

diff --git a/README.md b/README.md
index 4185179..45a3e34 100644
--- a/README.md
+++ b/README.md
@@ -68,7 +68,7 @@ This library provides two approaches:
 ```C
 struct box {
     int empty;
-    float content;
+    double content;
     char* error_msg;
 };
 ```
@@ -131,9 +131,9 @@ int main(){
     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);
-    float b = 2 * a;
-    float c = b / a;
+    double a = sample_to(1, 10, seed);
+    double b = 2 * a;
+    double c = b / a;
 
     printf("a: %f, b: %f, c: %f\n", a, b, c);
     // a: 0.607162, b: 1.214325, c: 0.500000
@@ -153,7 +153,7 @@ vs
 #include <stdlib.h>
 #include <stdio.h>
 
-float draw_xyz(uint64_t* seed){
+double draw_xyz(uint64_t* seed){
     // function could also be placed inside main with gcc nested functions extension.
     return sample_to(1, 20, seed);
 }
@@ -164,9 +164,9 @@ int main(){
     uint64_t* seed = malloc(sizeof(uint64_t));
     *seed = 1000; // xorshift can't start with a seed of 0
 
-    float a = draw_xyz(seed);
-    float b = 2 * draw_xyz(seed);
-    float c = b / a;
+    double a = draw_xyz(seed);
+    double b = 2 * draw_xyz(seed);
+    double c = b / a;
 
     printf("a: %f, b: %f, c: %f\n", a, b, c);
     // a: 0.522484, b: 10.283501, c: 19.681936
diff --git a/examples/01_one_sample/example b/examples/01_one_sample/example
index 5038416938686b8b98b2981d837dfe1f044c2ceb..0395c5c0b5d39ae55b9019189cd89cf8a91aedf4 100755
GIT binary patch
delta 32
jcmbQRj&agD#tq><Y$^GrNja&LV|=bdnVT(qZ*T$t*MAKG

delta 30
icmbQVj&Z^|#tq><tZ6y<i6xU`eXc`To2`6rZ~_3nJ`6Yj

diff --git a/examples/01_one_sample/example.c b/examples/01_one_sample/example.c
index f1687ae..c0668a8 100644
--- a/examples/01_one_sample/example.c
+++ b/examples/01_one_sample/example.c
@@ -4,22 +4,22 @@
 #include <stdio.h>
 
 // Estimate functions
-float sample_0(uint64_t* seed)
+double sample_0(uint64_t* seed)
 {
     return 0;
 }
 
-float sample_1(uint64_t* seed)
+double sample_1(uint64_t* seed)
 {
     return 1;
 }
 
-float sample_few(uint64_t* seed)
+double sample_few(uint64_t* seed)
 {
     return sample_to(1, 3, seed);
 }
 
-float sample_many(uint64_t* seed)
+double sample_many(uint64_t* seed)
 {
     return sample_to(2, 10, seed);
 }
@@ -29,15 +29,15 @@ int main(){
 		uint64_t* seed = malloc(sizeof(uint64_t));
 		*seed = 1000; // xorshift can't start with 0
 
-    float p_a = 0.8;
-    float p_b = 0.5;
-    float p_c = p_a * p_b;
+    double p_a = 0.8;
+    double p_b = 0.5;
+    double p_c = p_a * p_b;
 
     int n_dists = 4;
-    float weights[] = { 1 - p_c, p_c / 2, p_c / 4, p_c / 4 };
-    float (*samplers[])(uint64_t*) = { sample_0, sample_1, sample_few, sample_many };
+    double weights[] = { 1 - p_c, p_c / 2, p_c / 4, p_c / 4 };
+    double (*samplers[])(uint64_t*) = { sample_0, sample_1, sample_few, sample_many };
 
-    float result_one = sample_mixture(samplers, weights, n_dists, seed);
+    double result_one = sample_mixture(samplers, weights, n_dists, seed);
 		printf("result_one: %f\n", result_one);
 		free(seed);
 }
diff --git a/examples/02_many_samples/example b/examples/02_many_samples/example
index 4b781cf33257928a6f489624af77bf0b03772ee7..67494bdb3a3696d3bc52171ef3d0d5e3dc93fe85 100755
GIT binary patch
delta 32
jcmdn7j&avI#tn@=Y$^GrNja&L+kCD=nVXY*uW<qZ=sOL~

delta 30
icmdnBj&a92#tn@=tZ6y<i6xWUeXc`Tn^SzRaRLC%+zkW(

diff --git a/examples/02_many_samples/example.c b/examples/02_many_samples/example.c
index 7f3faa0..18ec1eb 100644
--- a/examples/02_many_samples/example.c
+++ b/examples/02_many_samples/example.c
@@ -4,22 +4,22 @@
 #include "../../squiggle.h"
 
 // Estimate functions
-float sample_0(uint64_t* seed)
+double sample_0(uint64_t* seed)
 {
     return 0;
 }
 
-float sample_1(uint64_t* seed)
+double sample_1(uint64_t* seed)
 {
     return 1;
 }
 
-float sample_few(uint64_t* seed)
+double sample_few(uint64_t* seed)
 {
     return sample_to(1, 3, seed);
 }
 
-float sample_many(uint64_t* seed)
+double sample_many(uint64_t* seed)
 {
     return sample_to(2, 10, seed);
 }
@@ -29,16 +29,16 @@ int main(){
 		uint64_t* seed = malloc(sizeof(uint64_t));
 		*seed = 1000; // xorshift can't start with 0
 
-    float p_a = 0.8;
-    float p_b = 0.5;
-    float p_c = p_a * p_b;
+    double p_a = 0.8;
+    double p_b = 0.5;
+    double p_c = p_a * p_b;
 
     int n_dists = 4;
-    float weights[] = { 1 - p_c, p_c / 2, p_c / 4, p_c / 4 };
-    float (*samplers[])(uint64_t*) = { sample_0, sample_1, sample_few, sample_many };
+    double weights[] = { 1 - p_c, p_c / 2, p_c / 4, p_c / 4 };
+    double (*samplers[])(uint64_t*) = { sample_0, sample_1, sample_few, sample_many };
 
 		int n_samples = 1000000;
-		float* result_many = (float *) malloc(n_samples * sizeof(float));
+		double* result_many = (double *) malloc(n_samples * sizeof(double));
 		for(int i=0; i<n_samples; i++){
       result_many[i] = sample_mixture(samplers, weights, n_dists, seed);
 		}
diff --git a/examples/03_gcc_nested_function/example b/examples/03_gcc_nested_function/example
index 7c4bfc34ac18e1eaf0c433e810b0edfd293ab236..f669b125a2ba8a56dc561dbe0ff5ce84c488f779 100755
GIT binary patch
delta 32
jcmX@Hj`7qw#tmD2*i!OKlX6lg@A5efWp3{B-Ny+4^ZyQ6

delta 30
icmX@Lj`74g#tmD2SkrRy6H6xV_Bjn<ZSL{i#|Z$~&<&6P

diff --git a/examples/03_gcc_nested_function/example.c b/examples/03_gcc_nested_function/example.c
index 62c8e53..72eb925 100644
--- a/examples/03_gcc_nested_function/example.c
+++ b/examples/03_gcc_nested_function/example.c
@@ -8,22 +8,22 @@ int main(){
 		uint64_t* seed = malloc(sizeof(uint64_t));
 		*seed = 1000; // xorshift can't start with 0
 
-    float p_a = 0.8;
-    float p_b = 0.5;
-    float p_c = p_a * p_b;
+    double p_a = 0.8;
+    double p_b = 0.5;
+    double p_c = p_a * p_b;
 
     int n_dists = 4;
 		
-		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); } 
+		double sample_0(uint64_t* seed){ return 0; }
+		double sample_1(uint64_t* seed) { return 1; }
+		double sample_few(uint64_t* seed){ return sample_to(1, 3, seed); }
+		double sample_many(uint64_t* seed){ return sample_to(2, 10, seed); } 
 		
-    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 };
+    double (*samplers[])(uint64_t*) = { sample_0, sample_1, sample_few, sample_many };
+		double weights[] = { 1 - p_c, p_c / 2, p_c / 4, p_c / 4 };
 
 		int n_samples = 1000000;
-		float* result_many = (float *) malloc(n_samples * sizeof(float));
+		double* result_many = (double *) malloc(n_samples * sizeof(double));
 		for(int i=0; i<n_samples; i++){
       result_many[i] = sample_mixture(samplers, weights, n_dists, seed);
 		}
diff --git a/examples/04_sample_from_cdf_simple/example b/examples/04_sample_from_cdf_simple/example
index 84ed2c2917fd2899bf36904916402e6ee4bfdadf..e272595a5ff9678fcb1ab8c34d49782feacc82f9 100755
GIT binary patch
delta 43
ocmZ3mo^kPd#tpZ8*i!OKlX6lg-}hMxV|MyFL7AH`_+H=y0HC=Nv;Y7A

delta 40
ncmZ3yo^io?#tpZ8SkrRy6H6vP@L39Fb^AI&Seq~VUf=`(Yx5C^

diff --git a/examples/04_sample_from_cdf_simple/example.c b/examples/04_sample_from_cdf_simple/example.c
index 0d21ec5..7734adb 100644
--- a/examples/04_sample_from_cdf_simple/example.c
+++ b/examples/04_sample_from_cdf_simple/example.c
@@ -8,7 +8,7 @@
 #define NUM_SAMPLES 1000000
 
 // Example cdf
-float cdf_uniform_0_1(float x)
+double cdf_uniform_0_1(double x)
 {
     if (x < 0) {
         return 0;
@@ -19,7 +19,7 @@ float cdf_uniform_0_1(float x)
     }
 }
 
-float cdf_squared_0_1(float x)
+double cdf_squared_0_1(double x)
 {
     if (x < 0) {
         return 0;
@@ -30,17 +30,17 @@ float cdf_squared_0_1(float x)
     }
 }
 
-float cdf_normal_0_1(float x)
+double cdf_normal_0_1(double x)
 {
-    float mean = 0;
-    float std = 1;
+    double mean = 0;
+    double std = 1;
     return 0.5 * (1 + erf((x - mean) / (std * sqrt(2)))); // erf from math.h
 }
 
 // Some testers
-void test_inverse_cdf_float(char* cdf_name, float cdf_float(float))
+void test_inverse_cdf_double(char* cdf_name, double cdf_double(double))
 {
-    struct box result = inverse_cdf_float(cdf_float, 0.5);
+    struct box result = inverse_cdf_double(cdf_double, 0.5);
     if (result.empty) {
         printf("Inverse for %s not calculated\n", cdf_name);
         exit(1);
@@ -49,12 +49,12 @@ 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), uint64_t* seed)
+void test_and_time_sampler_double(char* cdf_name, double cdf_double(double), uint64_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_cdf_float(cdf_float, seed);
+        struct box sample = sampler_cdf_double(cdf_double, seed);
         if (sample.empty) {
             printf("Error in sampler function for %s", cdf_name);
         } else {
@@ -62,26 +62,26 @@ void test_and_time_sampler_float(char* cdf_name, float cdf_float(float), uint64_
         }
     }
     clock_t end = clock();
-    float time_spent = (float)(end - begin) / CLOCKS_PER_SEC;
+    double time_spent = (double)(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 double
+    test_inverse_cdf_double("cdf_uniform_0_1", cdf_uniform_0_1);
+    test_inverse_cdf_double("cdf_squared_0_1", cdf_squared_0_1);
+    test_inverse_cdf_double("cdf_normal_0_1", cdf_normal_0_1);
 
     // Testing samplers
     // set randomness seed
     uint64_t* seed = malloc(sizeof(uint64_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);
+    // Test double sampler
+    test_and_time_sampler_double("cdf_uniform_0_1", cdf_uniform_0_1, seed);
+    test_and_time_sampler_double("cdf_squared_0_1", cdf_squared_0_1, seed);
+    test_and_time_sampler_double("cdf_normal_0_1", cdf_normal_0_1, seed);
 
     // Get some normal samples using a previous approach
     printf("\nGetting some samples from sample_unit_normal\n");
@@ -89,12 +89,12 @@ int main()
     clock_t begin_2 = clock();
 
     for (int i = 0; i < NUM_SAMPLES; i++) {
-        float normal_sample = sample_unit_normal(seed);
+        double normal_sample = sample_unit_normal(seed);
         // printf("%f\n", normal_sample);
     }
 
     clock_t end_2 = clock();
-    float time_spent_2 = (float)(end_2 - begin_2) / CLOCKS_PER_SEC;
+    double time_spent_2 = (double)(end_2 - begin_2) / CLOCKS_PER_SEC;
     printf("Time spent: %f\n", time_spent_2);
 
     free(seed);
diff --git a/examples/05_sample_from_cdf_beta/example b/examples/05_sample_from_cdf_beta/example
index 853558f01f31952818063528b379fc2e387ec7a4..b1afd6f63c6f9d2951ab11adf7b77c25dea1ce6b 100755
GIT binary patch
delta 32
jcmdmSo^jWC#tm$-Y$^GrNja&Lxno~InVWya9^eE3;$02g

delta 30
icmdmWo^i){#tm$-tZ6y<i6xVHVqZX5n}5X~-~<53rwtPT

diff --git a/examples/05_sample_from_cdf_beta/example.c b/examples/05_sample_from_cdf_beta/example.c
index 1328ba2..bd7d28a 100644
--- a/examples/05_sample_from_cdf_beta/example.c
+++ b/examples/05_sample_from_cdf_beta/example.c
@@ -10,11 +10,11 @@
 #define TINY_BETA 1.0e-30
 
 // Incomplete beta function
-struct box incbeta(float a, float b, float x)
+struct box incbeta(double a, double b, double x)
 {
     // Descended from <https://github.com/codeplea/incbeta/blob/master/incbeta.c>,
     // <https://codeplea.com/incomplete-beta-function-c>
-    // but modified to return a box struct and floats instead of doubles.
+    // but modified to return a box struct and doubles instead of doubles.
     // [ ] to do: add attribution in README
     // Original code under this license:
     /*
@@ -60,17 +60,17 @@ struct box incbeta(float a, float b, float x)
     }
 
     /*Find the first part before the continued fraction.*/
-    const float lbeta_ab = lgamma(a) + lgamma(b) - lgamma(a + b);
-    const float front = exp(log(x) * a + log(1.0 - x) * b - lbeta_ab) / a;
+    const double lbeta_ab = lgamma(a) + lgamma(b) - lgamma(a + b);
+    const double front = exp(log(x) * a + log(1.0 - x) * b - lbeta_ab) / a;
 
     /*Use Lentz's algorithm to evaluate the continued fraction.*/
-    float f = 1.0, c = 1.0, d = 0.0;
+    double f = 1.0, c = 1.0, d = 0.0;
 
     int i, m;
     for (i = 0; i <= 200; ++i) {
         m = i / 2;
 
-        float numerator;
+        double numerator;
         if (i == 0) {
             numerator = 1.0; /*First numerator is 1.0.*/
         } else if (i % 2 == 0) {
@@ -89,7 +89,7 @@ struct box incbeta(float a, float b, float x)
         if (fabs(c) < TINY_BETA)
             c = TINY_BETA;
 
-        const float cd = c * d;
+        const double cd = c * d;
         f *= cd;
 
         /*Check for stop.*/
@@ -105,7 +105,7 @@ struct box incbeta(float a, float b, float x)
     return PROCESS_ERROR("More loops needed, did not converge, in function incbeta");
 }
 
-struct box cdf_beta(float x)
+struct box cdf_beta(double x)
 {
     if (x < 0) {
         struct box result = { .empty = 0, .content = 0 };
@@ -114,13 +114,13 @@ struct box cdf_beta(float x)
         struct box result = { .empty = 0, .content = 1 };
         return result;
     } else {
-        float successes = 1, failures = (2023 - 1945);
+        double successes = 1, failures = (2023 - 1945);
         return incbeta(successes, failures, x);
     }
 }
 
 // Some testers
-void test_inverse_cdf_box(char* cdf_name, struct box cdf_box(float))
+void test_inverse_cdf_box(char* cdf_name, struct box cdf_box(double))
 {
     struct box result = inverse_cdf_box(cdf_box, 0.5);
     if (result.empty) {
@@ -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), uint64_t* seed)
+void test_and_time_sampler_box(char* cdf_name, struct box cdf_box(double), uint64_t* seed)
 {
     printf("\nGetting some samples from %s:\n", cdf_name);
     clock_t begin = clock();
@@ -144,7 +144,7 @@ void test_and_time_sampler_box(char* cdf_name, struct box cdf_box(float), uint64
         }
     }
     clock_t end = clock();
-    float time_spent = (float)(end - begin) / CLOCKS_PER_SEC;
+    double time_spent = (double)(end - begin) / CLOCKS_PER_SEC;
     printf("Time spent: %f\n", time_spent);
 }
 
diff --git a/examples/06_gamma_beta/example b/examples/06_gamma_beta/example
index 20513bc6d08ac3b043606cf8d7d4ed0f2bed6c34..11e0a92f61ef2f8a590d5830b1bb4f0b90484951 100755
GIT binary patch
delta 32
jcmbQRmT}Tr#touAY$^GrNja&LrF>39nVVnvEaL<K&bbWk

delta 30
icmbQVmT|&b#touAtZ6y<i6xVzeNIDIo8R~>;{*V<P7Eag

diff --git a/examples/06_gamma_beta/example.c b/examples/06_gamma_beta/example.c
index 8c11342..3a7d422 100644
--- a/examples/06_gamma_beta/example.c
+++ b/examples/06_gamma_beta/example.c
@@ -14,15 +14,15 @@ int main()
     int n =  1000 * 1000;
 		/*
 		for (int i = 0; i < n; i++) {
-        float gamma_0 = sample_gamma(0.0, seed);
+        double gamma_0 = sample_gamma(0.0, seed);
         // printf("sample_gamma(0.0): %f\n", gamma_0);
     }
 		printf("\n");
 		*/
 
-		float* gamma_1_array = malloc(sizeof(float) * n);
+		double* gamma_1_array = malloc(sizeof(double) * n);
     for (int i = 0; i < n; i++) {
-        float gamma_1 = sample_gamma(1.0, seed);
+        double gamma_1 = sample_gamma(1.0, seed);
         // printf("sample_gamma(1.0): %f\n", gamma_1);
 				gamma_1_array[i] = gamma_1;
     }
@@ -30,9 +30,9 @@ int main()
 		free(gamma_1_array);
 		printf("\n");
 		
-		float* beta_1_2_array = malloc(sizeof(float) * n);
+		double* beta_1_2_array = malloc(sizeof(double) * n);
 		for (int i = 0; i < n; i++) {
-        float beta_1_2 = sample_beta(1, 2.0, seed);
+        double beta_1_2 = sample_beta(1, 2.0, seed);
         // printf("sample_beta(1.0, 2.0): %f\n", beta_1_2);
 				beta_1_2_array[i] = beta_1_2;
     }
diff --git a/squiggle.c b/squiggle.c
index 3d3f4e6..8832c36 100644
--- a/squiggle.c
+++ b/squiggle.c
@@ -1,4 +1,4 @@
-#include <float.h>
+#include <double.h>
 #include <limits.h>
 #include <math.h>
 #include <stdint.h>
@@ -11,7 +11,7 @@
 #define EXIT_ON_ERROR 0
 #define PROCESS_ERROR(error_msg) process_error(error_msg, EXIT_ON_ERROR, __FILE__, __LINE__)
 
-const float PI = 3.14159265358979323846; // M_PI in gcc gnu99
+const double PI = 3.14159265358979323846; // M_PI in gcc gnu99
 
 // Pseudo Random number generator
 uint64_t xorshift32(uint32_t* seed)
@@ -44,58 +44,58 @@ uint64_t xorshift64(uint64_t* seed)
 
 // Distribution & sampling functions
 // Unit distributions
-float sample_unit_uniform(uint64_t* seed)
+double sample_unit_uniform(uint64_t* seed)
 {
     // samples uniform from [0,1] interval.
-    return ((float)xorshift64(seed)) / ((float)UINT64_MAX);
+    return ((double)xorshift64(seed)) / ((double)UINT64_MAX);
 }
 
-float sample_unit_normal(uint64_t* seed)
+double sample_unit_normal(uint64_t* seed)
 {
     // See: <https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform>
-    float u1 = sample_unit_uniform(seed);
-    float u2 = sample_unit_uniform(seed);
-    float z = sqrtf(-2.0 * log(u1)) * sin(2 * PI * u2);
+    double u1 = sample_unit_uniform(seed);
+    double u2 = sample_unit_uniform(seed);
+    double z = sqrtf(-2.0 * log(u1)) * sin(2 * PI * u2);
     return z;
 }
 
 // Composite distributions
-float sample_uniform(float start, float end, uint64_t* seed)
+double sample_uniform(double start, double end, uint64_t* seed)
 {
     return sample_unit_uniform(seed) * (end - start) + start;
 }
 
-float sample_normal(float mean, float sigma, uint64_t* seed)
+double sample_normal(double mean, double sigma, uint64_t* seed)
 {
     return (mean + sigma * sample_unit_normal(seed));
 }
 
-float sample_lognormal(float logmean, float logsigma, uint64_t* seed)
+double sample_lognormal(double logmean, double logsigma, uint64_t* seed)
 {
     return expf(sample_normal(logmean, logsigma, seed));
 }
 
-float sample_to(float low, float high, uint64_t* seed)
+double sample_to(double low, double high, uint64_t* seed)
 {
     // Given a (positive) 90% confidence interval,
     // returns a sample from a lognormal
     // with a matching 90% c.i.
-    const float NORMAL95CONFIDENCE = 1.6448536269514722;
-    float loglow = logf(low);
-    float loghigh = logf(high);
-    float logmean = (loglow + loghigh) / 2;
-    float logsigma = (loghigh - loglow) / (2.0 * NORMAL95CONFIDENCE);
+    const double NORMAL95CONFIDENCE = 1.6448536269514722;
+    double loglow = logf(low);
+    double loghigh = logf(high);
+    double logmean = (loglow + loghigh) / 2;
+    double logsigma = (loghigh - loglow) / (2.0 * NORMAL95CONFIDENCE);
     return sample_lognormal(logmean, logsigma, seed);
 }
 
-float sample_gamma(float alpha, uint64_t* seed)
+double sample_gamma(double alpha, uint64_t* seed)
 {
 
     // A Simple Method for Generating Gamma Variables, Marsaglia and Wan Tsang, 2001
     // https://dl.acm.org/doi/pdf/10.1145/358407.358414
     // see also the references/ folder
     if (alpha >= 1) {
-        float d, c, x, v, u;
+        double d, c, x, v, u;
         d = alpha - 1.0 / 3.0;
         c = 1.0 / sqrt(9.0 * d);
         while (1) {
@@ -125,24 +125,24 @@ float sample_gamma(float alpha, uint64_t* seed)
     }
 }
 
-float sample_beta(float a, float b, uint64_t* seed)
+double sample_beta(double a, double b, uint64_t* seed)
 {
-    float gamma_a = sample_gamma(a, seed);
-    float gamma_b = sample_gamma(b, seed);
+    double gamma_a = sample_gamma(a, seed);
+    double gamma_b = sample_gamma(b, seed);
     return gamma_a / (gamma_a + gamma_b);
 }
 
 // Array helpers
-float array_sum(float* array, int length)
+double array_sum(double* array, int length)
 {
-    float sum = 0.0;
+    double sum = 0.0;
     for (int i = 0; i < length; i++) {
         sum += array[i];
     }
     return sum;
 }
 
-void array_cumsum(float* array_to_sum, float* array_cumsummed, int length)
+void array_cumsum(double* array_to_sum, double* array_cumsummed, int length)
 {
     array_cumsummed[0] = array_to_sum[0];
     for (int i = 1; i < length; i++) {
@@ -150,16 +150,16 @@ void array_cumsum(float* array_to_sum, float* array_cumsummed, int length)
     }
 }
 
-float array_mean(float* array, int length)
+double array_mean(double* array, int length)
 {
-    float sum = array_sum(array, length);
+    double sum = array_sum(array, length);
     return sum / length;
 }
 
-float array_std(float* array, int length)
+double array_std(double* array, int length)
 {
-    float mean = array_mean(array, length);
-    float std = 0.0;
+    double mean = array_mean(array, length);
+    double std = 0.0;
     for (int i = 0; i < length; i++) {
         std += (array[i] - mean);
 				std *= std;
@@ -169,20 +169,20 @@ float array_std(float* array, int length)
 }
 
 // Mixture function
-float sample_mixture(float (*samplers[])(uint64_t*), float* weights, int n_dists, uint64_t* seed)
+double sample_mixture(double (*samplers[])(uint64_t*), double* 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/
-    float sum_weights = array_sum(weights, n_dists);
-    float* cumsummed_normalized_weights = (float*)malloc(n_dists * sizeof(float));
+    double sum_weights = array_sum(weights, n_dists);
+    double* cumsummed_normalized_weights = (double*)malloc(n_dists * sizeof(double));
     cumsummed_normalized_weights[0] = weights[0] / sum_weights;
     for (int i = 1; i < n_dists; i++) {
         cumsummed_normalized_weights[i] = cumsummed_normalized_weights[i - 1] + weights[i] / sum_weights;
     }
 
-    float result;
+    double result;
     int result_set_flag = 0;
-    float p = sample_uniform(0, 1, seed);
+    double p = sample_uniform(0, 1, seed);
     for (int k = 0; k < n_dists; k++) {
         if (p < cumsummed_normalized_weights[k]) {
             result = samplers[k](seed);
@@ -200,7 +200,7 @@ float sample_mixture(float (*samplers[])(uint64_t*), float* weights, int n_dists
 // Sample from an arbitrary cdf
 struct box {
     int empty;
-    float content;
+    double content;
     char* error_msg;
 };
 
@@ -219,13 +219,13 @@ struct box process_error(const char* error_msg, int should_exit, char* file, int
 
 // Inverse cdf at point
 // Two versions of this function:
-//   - raw, dealing with cdfs that return floats
-//     - input: cdf: float => float, p
+//   - raw, dealing with cdfs that return doubles
+//     - input: cdf: double => double, p
 //     - output: Box(number|error)
 //   - box, dealing with cdfs that return a box.
-//     - input: cdf: float => Box(number|error), p
+//     - input: cdf: double => Box(number|error), p
 //     - output: Box(number|error)
-struct box inverse_cdf_float(float cdf(float), float p)
+struct box inverse_cdf_double(double cdf(double), double p)
 {
     // given a cdf: [-Inf, Inf] => [0,1]
     // returns a box with either
@@ -233,8 +233,8 @@ struct box inverse_cdf_float(float cdf(float), float 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;
+    double low = -1.0;
+    double high = 1.0;
 
     // 1. Make sure that cdf(low) < p < cdf(high)
     int interval_found = 0;
@@ -260,14 +260,14 @@ struct box inverse_cdf_float(float cdf(float), float p)
         int convergence_condition = 0;
         int count = 0;
         while (!convergence_condition && (count < (INT_MAX / 2))) {
-            float mid = (high + low) / 2;
+            double mid = (high + low) / 2;
             int mid_not_new = (mid == low) || (mid == high);
-            // float width = high - low;
+            // double width = high - low;
             // if ((width < 1e-8) || mid_not_new){
             if (mid_not_new) {
                 convergence_condition = 1;
             } else {
-                float mid_sign = cdf(mid) - p;
+                double mid_sign = cdf(mid) - p;
                 if (mid_sign < 0) {
                     low = mid;
                 } else if (mid_sign > 0) {
@@ -288,7 +288,7 @@ struct box inverse_cdf_float(float cdf(float), float p)
     }
 }
 
-struct box inverse_cdf_box(struct box cdf_box(float), float p)
+struct box inverse_cdf_box(struct box cdf_box(double), double p)
 {
     // given a cdf: [-Inf, Inf] => Box([0,1])
     // returns a box with either
@@ -296,8 +296,8 @@ struct box inverse_cdf_box(struct box cdf_box(float), float 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;
+    double low = -1.0;
+    double high = 1.0;
 
     // 1. Make sure that cdf(low) < p < cdf(high)
     int interval_found = 0;
@@ -332,9 +332,9 @@ struct box inverse_cdf_box(struct box cdf_box(float), float p)
         int convergence_condition = 0;
         int count = 0;
         while (!convergence_condition && (count < (INT_MAX / 2))) {
-            float mid = (high + low) / 2;
+            double mid = (high + low) / 2;
             int mid_not_new = (mid == low) || (mid == high);
-            // float width = high - low;
+            // double width = high - low;
             if (mid_not_new) {
                 // if ((width < 1e-8) || mid_not_new){
                 convergence_condition = 1;
@@ -343,7 +343,7 @@ struct box inverse_cdf_box(struct box cdf_box(float), float p)
                 if (cdf_mid.empty) {
                     return PROCESS_ERROR(cdf_mid.error_msg);
                 }
-                float mid_sign = cdf_mid.content - p;
+                double mid_sign = cdf_mid.content - p;
                 if (mid_sign < 0) {
                     low = mid;
                 } else if (mid_sign > 0) {
@@ -365,23 +365,23 @@ 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), uint64_t* seed)
+struct box sampler_cdf_box(struct box cdf(double), uint64_t* seed)
 {
-    float p = sample_unit_uniform(seed);
+    double p = sample_unit_uniform(seed);
     struct box result = inverse_cdf_box(cdf, p);
     return result;
 }
-struct box sampler_cdf_float(float cdf(float), uint64_t* seed)
+struct box sampler_cdf_double(double cdf(double), uint64_t* seed)
 {
-    float p = sample_unit_uniform(seed);
-    struct box result = inverse_cdf_float(cdf, p);
+    double p = sample_unit_uniform(seed);
+    struct box result = inverse_cdf_double(cdf, p);
     return result;
 }
 
 /* Could also define other variations, e.g.,
-float sampler_danger(struct box cdf(float), uint64_t* seed)
+double sampler_danger(struct box cdf(double), uint64_t* seed)
 {
-    float p = sample_unit_uniform(seed);
+    double p = sample_unit_uniform(seed);
     struct box result = inverse_cdf_box(cdf, p);
 		if(result.empty){
 			exit(1);
diff --git a/squiggle.h b/squiggle.h
index aa3399b..6e2106b 100644
--- a/squiggle.h
+++ b/squiggle.h
@@ -8,31 +8,31 @@
 uint64_t xorshift64(uint64_t* seed);
 
 // Basic distribution sampling functions
-float sample_unit_uniform(uint64_t* seed);
-float sample_unit_normal(uint64_t* seed);
+double sample_unit_uniform(uint64_t* seed);
+double sample_unit_normal(uint64_t* seed);
 
 // Composite distribution sampling functions
-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);
+double sample_uniform(double start, double end, uint64_t* seed);
+double sample_normal(double mean, double sigma, uint64_t* seed);
+double sample_lognormal(double logmean, double logsigma, uint64_t* seed);
+double sample_to(double low, double high, uint64_t* seed);
 
-float sample_gamma(float alpha, uint64_t* seed);
-float sample_beta(float a, float b, uint64_t* seed);
+double sample_gamma(double alpha, uint64_t* seed);
+double sample_beta(double a, double b, uint64_t* seed);
 
 // Array helpers
-float array_sum(float* array, int length);
-void array_cumsum(float* array_to_sum, float* array_cumsummed, int length);
-float array_mean(float* array, int length);
-float array_std(float* array, int length);
+double array_sum(double* array, int length);
+void array_cumsum(double* array_to_sum, double* array_cumsummed, int length);
+double array_mean(double* array, int length);
+double array_std(double* array, int length);
 
 // Mixture function
-float sample_mixture(float (*samplers[])(uint64_t*), float* weights, int n_dists, uint64_t* seed);
+double sample_mixture(double (*samplers[])(uint64_t*), double* weights, int n_dists, uint64_t* seed);
 
 // Box
 struct box {
     int empty;
-    float content;
+    double content;
     char* error_msg;
 };
 
@@ -43,11 +43,11 @@ struct box {
 struct box process_error(const char* error_msg, int should_exit, char* file, int line);
 
 // Inverse cdf
-struct box inverse_cdf_float(float cdf(float), float p);
-struct box inverse_cdf_box(struct box cdf_box(float), float p);
+struct box inverse_cdf_double(double cdf(double), double p);
+struct box inverse_cdf_box(struct box cdf_box(double), double p);
 
 // Samplers from cdf
-struct box sampler_cdf_float(float cdf(float), uint64_t* seed);
-struct box sampler_cdf_box(struct box cdf(float), uint64_t* seed);
+struct box sampler_cdf_double(double cdf(double), uint64_t* seed);
+struct box sampler_cdf_box(struct box cdf(double), uint64_t* seed);
 
 #endif
diff --git a/test/test b/test/test
index 708e1fffaaafad773cdd48fedd1ba858f1b6889f..f6a849cc61c3ce0206fae4dda622ad76896d46e7 100755
GIT binary patch
delta 3950
zcmYjU3s6+o8NPSfRaV)(OJre}=ZcqjL0}gn8d6bKK<)~rID%P}sL3FchmP7nOhyt=
zEV{Ks%-=~(OfBs=X-%Vpp_#&@i3`MclA_II;`qV@>qRsm7$4Di`<;8&$DKL*pYy-Y
zf4={m|D65%KK{dfe2ar`+wL4PuSn#yZv{p2Zu4v5`slTGAMfTPJx}Pv2b2Pa^HDQn
z7Csdhx2EXp&QF#sxK!BK_R++h?R8<Ik#C{DjE&(-+dIdu9ak?mUa={EQdLzx%qd|U
zCkGD8frH8tp5sOZn`JafqoyHHIWVL&j@9|6_3@|>Ss+G0(eyPK$_iatCPoPE?pLjX
zbObhn)UCX4(9L(>!$gf=*s_Ap_EcvINt<ARueybTk42X^#X+mtB4{L6Mo_BVkjrt`
zTjb{4RCMHKQ!*ntw`%P2R@1qe=>IlYV^;fwqz2}wCq$utn+cn<M`+4~c>WWTCnTnR
z0v1<0rCe^zeVcXVz{Q#fd28-F7$|#wQB|c1Qywxk)3ynj^(}y4A>eT+c0}z99Fzk+
zCC#mjjlkJ#3=z4p$EZ{p87h(acS~}9sJ;Ao6+|Tf>kWb9?hgd1Z&avEnTkZ-2Q5+*
zl~|jGXNeerSTADX2jHF)BxRIlMp>Rz6PS2$*9cO<`qjeF$qr<RF{4B{S}7WpHxW5f
zP$`B(9Av6zXhT#~{Q<axUJqIs)MQ=F2s4}IIGhP(T@VGU1-EH8hOZU_wd1f+cHPw+
zT0YQbw-t!SM^s;Ds?(u*i>bZ>RV37LU>t&4@~pSg`|QB48Lv4vzKtm?4Eb6Dv%z+-
zM;JGkJa&8iB5%1b@XG;)sJE<~VO^?W6OP{!(Wb0}tuJs0`!R$&G7w5I)&LFs7Kq`1
zTFZj}einQ=XKK%dfWN7nuxJoxZ;yYFjhljWYgC|r+odU`R!v%^IRk`nnOz4>>|HGu
zPj#UVs%a4k(N4aAQlm4%s}V91cAl0-d-z4PEm|u0ELKCiN^}6!N-<qST~ea&s47yj
zQk0qJ9ay{;v?!Dv3W0x6X!;c3!Cg1pITEq15pDf)b3(MHzZi=I3%|-gDJ@2--+6y7
z<PvjD)8_IhH9xDWkd()ENl*@))3%b_c+l9v_7rpq?<xzRQj&AA_6(9%`bg5eO9wu<
zt)*?56Z!cN{lPQH)@=vJVd|BlLwN^0Z6ywdd-e{*FmMMaLrCH%81qMdE%_ZOd@;T0
zSk0fNA03jl`j)C9mZKv=g{j*~`xESx<Ft%x!VGNf${SSbTxPg>llD7({MU5Hxr=wx
zd$E(PpWb+A7~{|k`(th6zKaPN{)g_xuH`dmP2B5-^U<_3_HpaWKRq-#9R0v#OLWMj
zo}P|R8~e?54Iy&J<?&}J7(Xd$4>(v(iD*Q+jiMlRDeutD_=)xoJ93iMt8D_vDlCIh
z97kD>DDo!c@K4kF1j%x0^x>r0)Rr)L{7+EQ)^=r-zE4Qyy%dp{&(EjQ#N_(tArm%s
zqZBew&~Bhzh4zsVq0&+bDd06R3jWCvVevJ{egt?25k!rcuI0g1EoAmqV|R`S3yun#
z-6#gkDpR*Yf|jCKEXMARtP`UK_9F426H{%$t{9btv~DcMbewJ^Cf0v6Ar$R1Bicf)
zvP`vu2^T?BEFn`F6&CcvS(O+rY;Hy9z{>QxW>&!}39|UIZp#(~Bo~*ky;_08gQFP*
zKSwbg>Wtr13@gEXOW623w%fOc$pNCIP4zA43KDRx=#>NAa@Y43ppKr#=(;>sAAujK
z^}-NLlqhUD1SgnLTTzg$=()$(p*@8z8jJqF=$=4327K)6AO)LkMawJb1Mw3me7)=H
zSHoeYC@J3#>$!whgs7i3zzi@^iX2J@r0|pA(QL6@91!-x%PCM)SX5tWMDTaRe}U+N
zkJIQSq(!os{5_ZnhgWFYg5#nAgX#eIa0WewB|Wsw87HRT=3&^r4HHBZ0c`zUN`zjs
z*~pf#j&@Bf6jF2rR{a6~KWJp)tdft0RCQnxbP;WYs9ePm!A6J!pP_={e2xl<KbOma
zU)W_2=eTW@9K9~HizhGy9*YlqKlyHZU2;Ef`3CSdKdm<2HV928XqA|jE=ZRJ=@_C?
zghR*5p^K4VKPDVHiMajPqJEa;ZaJs7_A45m<f$JOmVOa9DoDrW=7Md;hjlwfZai$n
zU5J$X;jiMToO8JL612Y*GEWH-Qnujax;VTVQ0rNWPs7DvAM0cxt387o*`18}5i}ke
z%El`iaOkpoNH@#FifLoYw}#z=v|+NBc0Mli+0_5|EPftMnq0}hK@F3u`Kfei^2Q1I
zxVbv7qe7MhE`83Za)IiksWb1H(XALzRjrIFtI<<_R<K9dvjym|Gy&Cxa^MaEl30S6
zD^Sa9x-7k9u$gI5>K49<&Zib<8sP2KNC>m<;_)W9M}^J*(%(W-C$C*5_EivM-B^O6
zTTA}5_2vZ_!@xeD&ZN1bJle>AT4<ULb>=8JoNQCx;t$cEr?gJIItE@hcj+!&oj%>=
z3UGCot7v)pn(^6SYh!yKolZ|kdFsKq!^Cm8y$6zUT~)ULNmr*|_pyr-GUl20Lq1T7
z9Y)I?3H(-iE5m2@GTtyQ6vi{enMlouO)17Z(7VVN_?F%6KY7c%fxG2$bH+@t7nY>5
z9T=_*wa?1z<|%0QeC)|G#H3_qRj0ccqO+g+ie_<~%Zt`U^3=Yms}A1#c>PH0Y(6gR
zSJ$C|R|GeFMpY#QH+)uAOF_AFs`>&OgBlI=(i5@uZT(O|wdw<%32M6tRnP{|D$ut7
zscHkL>yoPe0dx%#+Xfm$mcIlo9aPm}P!1=_g%Vo|n#bN<C<A4nd7v+_@ikSg19hPW
z?f~t@ZPFjcrEzO;9XaqGS$-1}gZAR`41(qXu5UowKow9IWlYQGgS2+qWakd73(vs^
zA4)n`gVsX(rxn*1pmzjqW}m9=WueWr+TJqEu{xT==2>08H+rp7lgVrKG>$2-Ht=Dk
zhH+NUT&v`@x}Y=1YMWyXFB}gu7n%)d^<%b8=yQcuTa&?Obs2tRwdGr}%^*jc2AziQ
zRW+3{aPzItHmAM+aIEQ6LZ<_~=XD(pj|FZGcu_ylzU)}L1oU5iMAQn=COVrP$9K?e
z{9YmZ^hsHMV854Pi9i|+^Q?Jc`w+X$YYm?>-mnWjpEZ0wBFCeQOUJ^$qZg+qrH13!
z?+N*4;Vy?{H+a(?;njfGOZ%r!NT|E2s_~ee3*&gJ%df3c#&9<@8tB6GVnhEhWzLw&
zm(p`HX7J(k?u^x`;p|EAh`T-nS@j5yW#l?|LCSiv$f!U<*YYN`2hE*Ea|<pz_LY)1
zcQ)^D-;le^$ouJn*Uft<a!w2{QR<u_XCJPUM?btdI0K-yQC2?wqLDAZ2!GeM=P%=H
z=vw|d*uGd$1pfYlB6CCNKjB{D3pblR%oN{V+FpnW=|rKMA0~b-^WvQA=G(|Ox5(@f
zb**@MckVO%8oCCFhb9)qI7?G>$zOp8`;LG<#*$16ied~l2fb93Y3R$K_lh!|t$3BP
zPiBbXn*tD+ZYy0Va`PNT$}!Hi(8}Dp3RfJ`=_9WkgPAL2H{VJDIf>syd*sc0A7%N9
zoDIpk`u)t>jk;$KZTH0hi9J3ypVxlM_km$d)#}>kSFWH<C1*q3_Uie`V^HP;3*Baa
irtZg2T?>n_LUWnh{P6Fww`_)43T+greM8yb4F3n+M_v5@

delta 3955
zcmZ8k4^&iT7JqM;i5X?)4Z<)pATZcd1_A~O6kEv|0eK`E#87}_TE@{bwH0eSwSiK=
z9HO9G+h)5i?%CYj?G#13E!yS?YhtUZY`a<XXENK~hz0@`ruFvsy#eIxJLk=P-~IR9
z`@8qu@6DfH<zIf4Z@2OO4*Rg-xdcwVo})<KWoXfdRxRDdyo--?-w;|(Q}Ss(kKU2I
zv21n3%in$3z4ydj$u+@_j=P@gel$$f@}0CN+{&-*?g@WjR(<(_^{|#2&)ekZ6-AMb
zaB>*dz9W+FuzU~CaYK3`k!u{{>%3!v+4y-P$ARhwK_%a?Tpq5jjYmT3F(I-*j9?1;
z!H{!8Qm+^xxOyMlFF5;QBRG5I-)K1IbNee+=H%(0uyNx*Y+3%*cl6Z?i3J=-U2|e8
z_Qqk8Hls@=dASHrv5YR^xbb$WEzh3Egt|G5ES{*dNUgf}A(Ni)*BO+6kl4UL>6@g`
z4>Dm}&KZh}n#;dLxlsuzFM-8vzpg@R%-hKxCEr(d5mIYj6BhCdKPrm69$UVv`z6&!
zWz=(qJg!^?N{4uheTOAqUrAdhV+)Cov)vdXQe&T1&ek$iBJ;15EdFx)xp&;l1g->m
zwh0_JwGf<vF`;ID0t!0?TBIn-sb&?=5-|d?{u_xOL<u;-DUVTFv@v%mo{6gujiMC8
zc&-(O2mXe7F=n&~N1H^gycv<B1T~@_Vn0(oMh`?s*S`jL(CfqPg<5V%O$#%#@c^6&
z-a`=iYXz6?Su7u10n`q_Mt<sw>d^Qrb$4@tsGXs@k*O|(s?1dHhAJ!N!C44u^&OrX
z&mF;~jMtVox05O4Px@K`v!QOV&oFK}IV_g?3Xiu!^8I*l3h7D}89jbWM6<jRmj8#k
z8~K<a!*7A8pF%%#Ci<UE0WJrPEPZp%6LV))VIpKZCR7&C+ZKmqgGtp|uQ~;U%%YAi
zS}bv^U?(`wqiREf_ZuK}VN|F(JSwb7UqVSS>H1_;3Pt;nZi~s}S5b3}lm82K#l)nW
zE-MPEut}7dvv%aG?kNi1cPH^aEHu5w+NUjNKaIsCnZQOqqwHnU%iWY>b#C4j%gn0R
zizruhjhKyJub3%s`$5IKawgs<E~$8D+Bps`Hi<U*I5=tqHWtANxC-9LUiP-UyP+h%
z3h_`+ZN^j^*S1r@`aIu3zp%~c@1xzeE&OFN*`21Q2}MCq9Yx~0Gj=)^Z=ou?an@cq
zK|Fr>x3txMr^Yxzr|iZ2Hxw7Uo4<wj$0nNwzq@L<$fg=<<IJ-rtdoYrlo0m-zmnSH
z9@Z!^bRc$tsrB2dCSS))o4gz|Y4SVTJU2Cbd|Xk4L^s_M=jMN;z}#!2kAj0NOGGXD
zMk@-=L-JqA8lPl2Z2`hrC)J*UtjZV;29?>OX<htcelxuo?=)T-ySixwos3VOZGn=S
zTf{hN5>ohel$Mas-%nc;=GE_vMn~+_;QeEu;2y+10PX%!p~m<qqySwPqSYT^;~Fsy
z(m}x2hsxB5j#?q(WG%+xs8DuPXx;`-%;Hv9GH>WDyy?}aRZC}ZADaY;OyVDtgw$Rg
zlF*G&LSjPw2T_wz){m+Qu-4QKf(dg!Mvagmj|pXG;b^_67n)^^Utpx$qncH*5rFDZ
z$aVl5070}DI<N|O;oWk=4-%{D#e1UL1}nieA?$pH^=xaH<fD$H1RAs@(D#Ybjni#}
ziOj1TM8Xx@p;i($_sfMB(Rr<5*ySjEtGiV9BeDpQo=!@t7g}VvW$Cjdk$YAA-E-JJ
zx=s+dL43|$zz0j!Ay_~PXqx5sM_8X%=^le7fHbZNQHZcocdF^fc>Cu=+u(zhz>Alo
z2S%Y_sw2eFdiAA)7|5$*=S9`0IXx=m3)4pRWrNCTlsSm~25?*AB3aZ<41`K}t||~6
zctUYq*jH%E$C05Ti3R}P7?GK)aX7hgWb`%$c$Tho>YQZg&0tx7444OlVR%JZm3|J-
zDn8&Il8=Q3Xas6zj=!BAODYs}Aq92-pp30FkaSbYlf#M<G-1RejR;XbjwRxb5Q7iH
z9$hphB>cdJEe<q{Ph=V)U&5S%Q+D?FhQY(|VQ<fw%iT5e&ho~yC`;QarSXzRX!?Ry
zim497(<?a73C{O~ch6M_7LWmWsFYy_a+r=>ySN_R?K>(s|1~D8IVQCgJguEPlCe`b
ztESX=M2nfK9zDR#57c@A2QaNZDJ?#-{X-ajAY^=uUC_cuA4-68Ag6W++ztg8c3Eih
zzThS{)9@vO{+LYZmh+L*IpI`drtMTZ|FouVh{}>Zbl(Dr*U_5`ZsKEUY(WhsP-${4
zuc5b+cSc3w2<sWgP%rU){SKq@>r~>re#!d=Jf;K{MLnA3rFfFxlhH%djjv)`F8MAa
z7$++ea|LP?XutDb%_n+VkkZ26O0T7?$QZ!4+=G!2<{ZP5;2IN}-wJ&WIeR$uM6u*S
zkoDbaOzT2&raodw!WssaczPw(5zVQq%0Z!NDb$&x(Qpc+KE*#nTdwO&`ujET+I%SF
z($OP@T;_lqa+yb&X<KJ!gRQPDUG(R)`1y;c*DX4Z!!aJ5hvuklN8dVnZV&m`Oc&B_
zHtc{rSc*QQ23tJ8n%1Nj8zLER1Sb>YX@<0PB0Y9~96ks=w-x*T<Ebe2ee3ahd{-)@
zw)6<F%S%#NAL`|U-KLCQ{*BBm4Gso&ExqwL-NbPY5AH*ZLDKbsW$c0D^lp`KX7~Mv
ze|qDQ(axoOT-eXX!NjKp7x+X`oCq}VsiKsEj(~1p>(3PB3QyHH#MU>Qg$gEN7w8hu
z0B9*_=~s%f9<<@ViqZf&0{Sed8%6E{ZAZnAfjWj1WdziWcgum%y%y9BYR1Iyf|i1A
zU~3G#M?u|~ufGNz!3lLXj7#O(am3m1vC{B8HU=HRIqU~r3+O%oJqap<x+y+8pYNh|
z*~#`6<b`j<=?g~%SBJZS_GGUxv9H*E+)DzAQm3JJvhAxH_^>$5ER%D&$>A~CmYK}U
zO!~sv&~f6{1CdYQaTlJsLX){kv)p8B3@bKeYIc}x`6jc+g!V!TxLJ>s(z=|e)J<m;
zC55qZF-Z+1E#TdyayZTgx)r>$XJ~g$tYzX$MY(T=s1qVD9nXp5TPTR%w@J6~+AJr~
zy#RDr7>$~nP42KEgv|mdoUPf7r(%<SCE~_c7?%drWvW`3m|{DxDC}Hh7Yk)AByWP3
zt#U%VI`9ItXJJ(QlNS`lh26QZDBk3#RCANC`~x&P>D0m%8tp|&SadzVmC6<^;+^DM
zv?ayKK161?3n0kHXYdLj8wW2yNx4N@J-&44qC2YFZ`h|bbmIJD*IM%AE#+IfcjVov
zMb=+<T)dkims$B-N?A5+AHmr&^TIa=r$Igr%F4$-H7d?8!oRkT{5$!_>0<svu-#Zt
z#5d6E1x1FnlYb2_5MS7Aa5Ga}b<|Oa4e3;&i$6*Hax0%n_T?_Vor;$i8A`>FmWy6k
zemgv0grt;`imdjA`60>gfe5?$LFXWnOe!m~YDTPdZ&8NEoI$&aGVBfbl4Td0DeSrq
zATZr7`l85%{URl+-9MR`E2P4$nAGu;N3tTCDv8DClM?w|v`=d0{ghQ)Wbd38QlF}<
z1up^fTuL3qR_N?2cJYqxkBawdbhVqdtgE9ZN<Ixeb#GWX?;8FdYAwg#<cyG)7V0Z6
b0%o(lF2mLT$obwyhTO?E$?b0R{!#NkJn&y_

diff --git a/test/test.c b/test/test.c
index 7e472f6..2c8f909 100644
--- a/test/test.c
+++ b/test/test.c
@@ -4,22 +4,22 @@
 #include <stdlib.h>
 #include <stdio.h>
 
-#define N  1000 * 1000
+#define N  100
 
 void test_unit_uniform(uint64_t* seed){
-	float* unit_uniform_array = malloc(sizeof(float) * N);
+	double* unit_uniform_array = malloc(sizeof(double) * N);
 	
 	for(int i=0; i<N; i++){
 		unit_uniform_array[i] = sample_unit_uniform(seed);
 	}
 	
-	float mean = array_mean(unit_uniform_array, N);
-	float expected_mean = 0.5;
-	float delta_mean = mean - expected_mean;
+	double mean = array_mean(unit_uniform_array, N);
+	double expected_mean = 0.5;
+	double delta_mean = mean - expected_mean;
 
-	float std = array_std(unit_uniform_array, N);
-	float expected_std = sqrt(1.0/12.0);
-	float delta_std = std - expected_std;
+	double std = array_std(unit_uniform_array, N);
+	double expected_std = sqrt(1.0/12.0);
+	double delta_std = std - expected_std;
 	
 	printf("Mean of unit uniform: %f, vs expected mean: %f, delta: %f\n", mean, expected_mean, delta_mean);
 	printf("Std of unit uniform: %f, vs expected std: %f, delta: %f\n", std, expected_std, delta_std);
@@ -40,23 +40,22 @@ void test_unit_uniform(uint64_t* seed){
 
 }
 
-void test_uniform(float start, float end, uint64_t* seed){
-	float* uniform_array = malloc(sizeof(float) * N);
+void test_uniform(double start, double end, uint64_t* seed){
+	double* uniform_array = malloc(sizeof(double) * N);
 	
 	for(int i=0; i<N; i++){
 		uniform_array[i] = sample_uniform(start, end, seed);
 	}
 	
-	float mean = array_mean(uniform_array, N);
-	float expected_mean = (start + end) / 2; 
-	float delta_mean = mean - expected_mean;
+	double mean = array_mean(uniform_array, N);
+	double expected_mean = (start + end) / 2; 
+	double delta_mean = mean - expected_mean;
 	
-	float std = array_std(uniform_array, N);
-	float expected_std = sqrt(1.0/12.0) * fabs(end-start);
-	float delta_std = std - expected_std;
+	double std = array_std(uniform_array, N);
+	double expected_std = sqrt(1.0/12.0) * fabs(end-start);
+	double delta_std = std - expected_std;
 	
-
-	float width = fabs(end - start);
+	double width = fabs(end - start);
 	if(fabs(delta_mean) > width * 1.0/1000.0){
 		printf("[-] Mean test for [%.1f, %.1f] uniform NOT passed.\n", start, end);
 		printf("Mean of [%.1f, %.1f] uniform: %f, vs expected mean: %f, delta: %f\n", start, end, mean, expected_mean, mean - expected_mean);
@@ -67,6 +66,10 @@ void test_uniform(float start, float end, uint64_t* seed){
 	if(fabs(delta_std) > width * 1.0/1000.0){
 		printf("[-] Std test for [%.1f, %.1f] uniform NOT passed.\n", start, end);
 		printf("Std of [%.1f, %.1f] uniform: %f, vs expected std: %f, delta: %f\n", start, end, std, expected_std, std - expected_std);
+		for(int i=0; i<N; i++){
+
+			printf("%.1f, ", uniform_array[i]);
+		}
 	}else {
 		printf("[x] Std test for unit uniform PASSED\n");
 	}
@@ -81,9 +84,9 @@ int main(){
    
 		test_unit_uniform(seed);
 
-		for(int i=0; i<100; i++){
-			float start = sample_uniform(-10, 10, seed);
-			float end = sample_uniform(-10, 10, seed);
+		for(int i=0; i<1; i++){
+			double start = sample_uniform(-10, 10, seed);
+			double end = sample_uniform(-10, 10, seed);
 			if ( end > start){
 				test_uniform(start, end, seed);
 			}