initial attempt on bc

buggy because wrong base for log, but it's a start

bc/estimate.bc

@@ -0,0 +1,34 @@
+p_a = 0.8
+p_b = 0.5
+p_c = p_a * p_b
+
+weights[0] = 1 - p_c
+weights[1] = p_c / 2
+weights[2] = p_c / 4
+weights[3] = p_c / 4
+
+/* We'll have to define the mixture manually */
+define mixture(){
+  p = sample_unit_uniform()
+  if(p <= weights[0]){
+    return 0
+  } 
+  if(p <= (weights[0] + weights[1])){
+    return 1
+  }
+  if(p<= (weights[0] + weights[1] + weights[2])){
+    return sample_to(1, 3)
+  }
+  return sample_to(2, 10)
+}
+
+/* n_samples = 1000000 */
+n_samples = 10000
+sum=0
+for(i=0; i < n_samples; i++){
+  /* samples[i] = mixture() */
+  sum += mixture()
+}
+sum/n_samples
+
+halt

bc/extra/beta.bc

@@ -0,0 +1,51 @@
+define sample_gamma(alpha){
+  /*
+  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
+  Note that the Wikipedia page for the gamma distribution includes a scaling parameter
+  k or beta
+  https://en.wikipedia.org/wiki/Gamma_distribution
+  such that gamma_k(alpha, k) = k * gamma(alpha)
+  or gamma_beta(alpha, beta) = gamma(alpha) / beta
+  So far I have not needed to use this, and thus the second parameter is by default 1.
+  */
+
+  if (alpha >= 1) {
+    d = alpha - (1/3);
+    c = 1.0 / sqrt(9.0 * d);
+    while (1) {
+      v=-1
+      while(v<=0) {
+        x = sample_unit_normal();
+        v = 1 + c * x;
+      }
+      v = v * v * v;
+      u = sample_unit_uniform();
+      if (u < (1 - (0.0331 * (x * x * x * x)))) { /* Condition 1 */
+        /* 
+        the 0.0331 doesn't inspire much confidence
+        however, this isn't the whole story
+        by knowing that Condition 1 implies condition 2
+        we realize that this is just a way of making the algorithm faster
+        i.e., of not using the logarithms
+        */
+        return d * v;
+      }
+      if (log(u, 2) < ((0.5 * (x * x)) + (d * (1 - v + log(v, 2))))) { /* Condition 2 */
+        return d * v;
+      }
+    }
+  } else {
+    return sample_gamma(1 + alpha) * p(sample_unit_uniform(), 1 / alpha);
+    /* see note in p. 371 of https://dl.acm.org/doi/pdf/10.1145/358407.358414 */
+  }
+}
+
+define sample_beta(a, b)
+{
+    /* See: https://en.wikipedia.org/wiki/Gamma_distribution#Related_distributions */
+    gamma_a = sample_gamma(a);
+    gamma_b = sample_gamma(b);
+    return gamma_a / (gamma_a + gamma_b);
+}

bc/makefile

@@ -0,0 +1,5 @@
+compute:
+	ghbc -l squiggle.bc estimate.bc
+
+time:
+	time ghbc -l squiggle.bc estimate.bc

bc/notes.md

@@ -6,7 +6,7 @@ https://pubs.opengroup.org/onlinepubs/9699919799.2018edition/utilities/bc.html
 ## gh-bc
 
 To build
-./configure
+./configure.sh -O3
 make
 sudo cp bin/bc /usr/bin/ghbc
 

bc/squiggle.bc

@@ -1,4 +1,4 @@
-scale = 32
+scale = 8
 /* seed = 12345678910 */
 pi = 4 * atan(1)
 normal90confidence=1.6448536269514727148638489079916
@@ -70,55 +70,3 @@ define sample_to(low, high){
   loghigh = log(high, 2)
   return e(sample_normal_from_90_confidence_interval(loglow, loghigh))
 }
-
-define sample_gamma(alpha){
-  /*
-  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
-  Note that the Wikipedia page for the gamma distribution includes a scaling parameter
-  k or beta
-  https://en.wikipedia.org/wiki/Gamma_distribution
-  such that gamma_k(alpha, k) = k * gamma(alpha)
-  or gamma_beta(alpha, beta) = gamma(alpha) / beta
-  So far I have not needed to use this, and thus the second parameter is by default 1.
-  */
-
-  if (alpha >= 1) {
-    d = alpha - (1/3);
-    c = 1.0 / sqrt(9.0 * d);
-    while (1) {
-      v=-1
-      while(v<=0) {
-        x = sample_unit_normal();
-        v = 1 + c * x;
-      }
-      v = v * v * v;
-      u = sample_unit_uniform();
-      if (u < (1 - (0.0331 * (x * x * x * x)))) { /* Condition 1 */
-        /* 
-        the 0.0331 doesn't inspire much confidence
-        however, this isn't the whole story
-        by knowing that Condition 1 implies condition 2
-        we realize that this is just a way of making the algorithm faster
-        i.e., of not using the logarithms
-        */
-        return d * v;
-      }
-      if (log(u, 2) < ((0.5 * (x * x)) + (d * (1 - v + log(v, 2))))) { /* Condition 2 */
-        return d * v;
-      }
-    }
-  } else {
-    return sample_gamma(1 + alpha) * p(sample_unit_uniform(), 1 / alpha);
-    /* see note in p. 371 of https://dl.acm.org/doi/pdf/10.1145/358407.358414 */
-  }
-}
-
-define sample_beta(a, b)
-{
-    /* See: https://en.wikipedia.org/wiki/Gamma_distribution#Related_distributions */
-    gamma_a = sample_gamma(a);
-    gamma_b = sample_gamma(b);
-    return gamma_a / (gamma_a + gamma_b);
-}