Migrate C source to SC.

sph-sc is a scheme-like S-expression syntax for C. It elements much of
the pain and repetition involved in writing C syntax.

* src/extent-sampling.c, src/gaussian-nd-random.c, src/nd-random.c,
src/oracles.c, src/utils.c: Delete files.
* src/extent-sampling.sc, src/gaussian-nd-random.sc,
src/macros/macros.sc, src/nd-random.sc, src/oracles.sc, src/utils.sc:
New files.
* CMakeLists.txt: Generate C source files from SC source files.

12 files changed, 568 insertions, 583 deletions

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 94264c7..fc3a296 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -2,8 +2,24 @@ cmake_minimum_required(VERSION 3.10)
 project(extent-sampling VERSION 0.1)
 
 find_package(GSL REQUIRED)
+# TODO: Make indent optional
+find_program(INDENT NAMES indent REQUIRED)
+find_program(SC NAMES sc REQUIRED)
+
+# Generate C source files from SC source files.
+file(GLOB SC_SOURCES "src/*.sc")
+foreach(sc_source ${SC_SOURCES})
+  get_filename_component(source_basename ${sc_source} NAME_WE)
+  set(c_file ${source_basename}.c)
+  add_custom_command(
+    OUTPUT ${c_file}
+    COMMAND ${SC} ${sc_source} ${c_file}
+    COMMAND ${INDENT} ${c_file}
+    DEPENDS ${sc_source}
+    VERBATIM)
+  list(APPEND C_SOURCES ${c_file})
+endforeach()
 
 include_directories("include")
-file(GLOB SOURCES "src/*.c")
 
-add_library(extentsampling SHARED ${SOURCES})
+add_library(extentsampling SHARED ${C_SOURCES})
diff --git a/src/extent-sampling.c b/src/extent-sampling.c
deleted file mode 100644
index eff3b3f..0000000
--- a/src/extent-sampling.c
+++ /dev/null
@@ -1,161 +0,0 @@
-/*
-  This file contains all the core extent sampling routines. From time
-  to time as research progresses, I keep refactoring this code to be
-  clearer and more efficient. This mostly involves throwing out old
-  code even when such old code may be required in the future. I
-  figured that keeping old code around is not worth the trouble
-  because it impedes clarity of thought.
-*/
-
-#include <math.h>
-#include <stddef.h>
-#include <gsl/gsl_blas.h>
-#include <gsl/gsl_integration.h>
-#include <gsl/gsl_rstat.h>
-#include <gsl/gsl_randist.h>
-#include "extent-sampling.h"
-#include "nd-random.h"
-#include "utils.h"
-
-#define VOLUME_MINIMUM_SAMPLES 100
-#define INTEGRAL_MINIMUM_SAMPLES 100
-#define CONFIDENCE_INTERVAL_FACTOR 1.96
-
-#define INTEGRATION_INTERVALS 1000
-
-#define SOLID_ANGLE_START 0.0
-#define SOLID_ANGLE_LAST 0.5
-#define SOLID_ANGLE_STEPS 100
-
-double volume (extent_oracle_t extent_oracle, double true_volume,
-               const gsl_rng* r, unsigned int dimension, double rtol,
-               gsl_rstat_workspace* stats)
-{
-  gsl_vector* x = gsl_vector_alloc(dimension);
-  double volume, error;
-  double vn = ln_volume_of_ball(dimension);
-  do {
-    random_direction_vector(r, x);
-    double extent = extent_oracle(x);
-    double volume_from_sample = exp(vn + dimension*log(extent));
-    gsl_rstat_add(volume_from_sample, stats);
-    volume = gsl_rstat_mean(stats);
-    error = CONFIDENCE_INTERVAL_FACTOR * gsl_rstat_sd_mean(stats) / volume;
-  } while ((error > rtol) || (rerror(volume, true_volume) > rtol) || (gsl_rstat_n(stats) < VOLUME_MINIMUM_SAMPLES));
-  gsl_vector_free(x);
-  return volume;
-}
-
-double volume_window (extent_oracle_t extent_oracle, double true_volume,
-                      const gsl_rng* r, unsigned int dimension, double rtol,
-                      unsigned int* number_of_samples)
-{
-  gsl_rstat_workspace* stats = gsl_rstat_alloc();
-  gsl_vector* x = gsl_vector_alloc(dimension);
-  double volume, error;
-  double vn = ln_volume_of_ball(dimension);
-  // This is the window length used in Volume.m of Lovasz-Vempala's
-  // code.
-  int window_length = 4*dimension*dimension + 500;
-  int accurate_estimates = 0;
-  do {
-    random_direction_vector(r, x);
-    double extent = extent_oracle(x);
-    double volume_from_sample = exp(vn + dimension*log(extent));
-    gsl_rstat_add(volume_from_sample, stats);
-    volume = gsl_rstat_mean(stats);
-    error = rerror(volume, true_volume);
-    if (error < rtol) accurate_estimates++;
-    else accurate_estimates = 0;
-  } while (accurate_estimates < window_length);
-  *number_of_samples = gsl_rstat_n(stats);
-  gsl_vector_free(x);
-  gsl_rstat_free(stats);
-  return volume;
-}
-
-double integral_per_direction
-(integrand_t integrand, const gsl_vector* direction, const gsl_rng* r, unsigned int n,
- double radius, double rtol, int* neval)
-{
-  double f (double r, void* params)
-  {
-    (*neval)++;
-    return gsl_pow_uint(r, n-1) * integrand(r, direction);
-  }
-
-  gsl_function gsl_f = {&f};
-  double result, error;
-  gsl_integration_workspace *w = gsl_integration_workspace_alloc(INTEGRATION_INTERVALS);
-  gsl_integration_qag(&gsl_f, 0, radius, 0, rtol, INTEGRATION_INTERVALS, GSL_INTEG_GAUSS15, w, &result, &error);
-  gsl_integration_workspace_free(w);
-  return surface_area_of_ball(n) * result;
-}
-
-double integral
-(integrand_t integrand, extent_oracle_t extent_oracle, double true_integral,
- const gsl_rng* r, unsigned int dimension, double rtol,
- gsl_rstat_workspace* stats)
-{
-  gsl_vector* x = gsl_vector_alloc(dimension);
-  double integral, error;
-  do {
-    int neval = 0;
-    random_direction_vector(r, x);
-    double extent = extent_oracle(x);
-    double integral_from_sample = integral_per_direction(integrand, x, r, dimension, extent, rtol, &neval);
-    gsl_rstat_add(integral_from_sample, stats);
-    integral = gsl_rstat_mean(stats);
-    error = CONFIDENCE_INTERVAL_FACTOR * gsl_rstat_sd_mean(stats) / integral;
-  } while ((error > rtol) || (rerror(integral, true_integral) > rtol) || (gsl_rstat_n(stats) < INTEGRAL_MINIMUM_SAMPLES));
-  if (error > rtol)
-    GSL_ERROR_VAL("integral failed to converge", GSL_ETOL, integral);
-  gsl_vector_free(x);
-  return integral;
-}
-
-double volume_cone (extent_oracle_t extent_oracle, const gsl_rng* r,
-                    const gsl_vector* mean, double omega_min, double omega_max,
-                    unsigned int number_of_samples, double* variance)
-{
-  unsigned int n = mean->size;
-  gsl_vector* x = gsl_vector_alloc(n);
-  double theta_min = solid_angle_to_planar_angle(omega_min * surface_area_of_ball(n), n);
-  double theta_max = solid_angle_to_planar_angle(omega_max * surface_area_of_ball(n), n);
-  double vn = ln_volume_of_ball(n);
-  gsl_rstat_workspace* stats = gsl_rstat_alloc();
-  for (int i=0; i<number_of_samples; i++) {
-    hollow_cone_random_vector(r, mean, theta_min, theta_max, x);
-    double extent = extent_oracle(x);
-    double volume_from_sample = exp(vn + n*log(extent));
-    gsl_rstat_add(volume_from_sample, stats);
-  }
-
-  if (variance) *variance = gsl_rstat_variance(stats);
-
-  gsl_rstat_free(stats);
-  gsl_vector_free(x);
-  return gsl_rstat_mean(stats) * (omega_max - omega_min);
-}
-
-double volume_experiment
-(extent_oracle_t extent_oracle, const gsl_rng* r,
- const gsl_vector* mean, unsigned int samples_per_cone,
- double solid_angle_factor, double solid_angle_threshold_exponent_factor,
- unsigned int* number_of_samples)
-{
-  int n = mean->size;
-  double volume = 0;
-  double omega_max = 0.5;
-  int N = 0;
-  do {
-    double omega_min = omega_max / solid_angle_factor;
-    volume += volume_cone(extent_oracle, r, mean, omega_min, omega_max, samples_per_cone, NULL);
-    N += samples_per_cone;
-    omega_max = omega_min;
-  } while (omega_max > pow(2, -solid_angle_threshold_exponent_factor*n));
-  volume += volume_cone(extent_oracle, r, mean, 0, omega_max, samples_per_cone, NULL);
-  N += samples_per_cone;
-  if (number_of_samples) *number_of_samples = N;
-  return volume;
-}
diff --git a/src/extent-sampling.sc b/src/extent-sampling.sc
new file mode 100644
index 0000000..c6220e6
--- /dev/null
+++ b/src/extent-sampling.sc
@@ -0,0 +1,155 @@
+(sc-include "macros/macros")
+
+(pre-include "stddef.h")
+(pre-include "gsl/gsl_integration.h")
+(pre-include "gsl/gsl_rstat.h")
+(pre-include "gsl/gsl_vector.h")
+(pre-include "nd-random.h")
+(pre-include "utils.h")
+(pre-include "extent-sampling.h")
+
+(sc-define-syntax (with-vector var size body ...)
+  (with-alloc var gsl-vector*
+              (gsl-vector-alloc size)
+              gsl-vector-free
+              body ...))
+
+(sc-define-syntax (with-rstats var body ...)
+  (with-alloc var gsl-rstat-workspace*
+              (gsl-rstat-alloc) gsl-rstat-free
+              body ...))
+
+(sc-define-syntax (with-integration-workspace var intervals body ...)
+  (with-alloc var gsl-integration-workspace*
+              (gsl-integration-workspace-alloc intervals)
+              gsl-integration-workspace-free
+              body ...))
+
+(pre-define CONFIDENCE-INTERVAL-FACTOR 1.96)
+
+(pre-let* (VOLUME-MINIMUM-SAMPLES 100)
+  (define (volume extent-oracle true-volume r dimension rtol stats)
+    (double extent-oracle-t double (const gsl-rng*) (unsigned int) double gsl-rstat-workspace*)
+    (declare volume double)
+    (let* ((vn double (ln-volume-of-ball dimension)))
+      (with-vector x dimension
+        (do-while (or (> (/ (* CONFIDENCE-INTERVAL-FACTOR
+                               (gsl-rstat-sd-mean stats))
+                            volume)
+                         rtol)
+                      (> (rerror volume true-volume) rtol)
+                      (< (gsl-rstat-n stats) VOLUME-MINIMUM-SAMPLES))
+          (random-direction-vector r x)
+          (gsl-rstat-add (exp (+ vn (* dimension (log (extent-oracle x)))))
+                         stats)
+          (set volume (gsl-rstat-mean stats)))))
+    (return volume)))
+
+(define (volume-window extent-oracle true-volume r dimension rtol number-of-samples)
+  (double extent-oracle-t double (const gsl-rng*) (unsigned int) double (unsigned int*))
+  (declare volume double)
+  (let* ((vn double (ln-volume-of-ball dimension))
+         ;; This is the window length used in Volume.m of
+         ;; Lovasz-Vempala's code.
+         (window-length int (+ (* 4 dimension dimension) 500))
+         (accurate-estimates int 0))
+    (with-rstats stats
+      (with-vector x dimension
+        (do-while (< accurate-estimates window-length)
+          (random-direction-vector r x)
+          (gsl-rstat-add (exp (+ vn (* dimension (log (extent-oracle x)))))
+                         stats)
+          (set volume (gsl-rstat-mean stats))
+          (cond
+           ((< (rerror volume true-volume) rtol)
+            (set+ accurate-estimates 1))
+           (else (set accurate-estimates 0)))))
+      (set (pointer-get number-of-samples)
+           (gsl-rstat-n stats))))
+  (return volume))
+
+(pre-let* (INTEGRATION-INTERVALS 1000)
+  (define (integral-per-direction integrand direction r n radius rtol neval)
+    (double integrand-t (const gsl-vector*) (const gsl-rng*) (unsigned int) double double int*)
+    (define (f r params) (double double void*)
+      (set+ (pointer-get neval) 1)
+      (return (* (gsl-pow-uint r (- n 1))
+                 (integrand r direction))))
+    (declare gsl-f (struct-variable gsl-function (function (address-of f)))
+             result double
+             error double)
+    (with-integration-workspace w INTEGRATION-INTERVALS
+      (gsl-integration-qag (address-of gsl-f)
+                           0
+                           radius
+                           0
+                           rtol
+                           INTEGRATION-INTERVALS
+                           GSL-INTEG-GAUSS15
+                           w
+                           (address-of result)
+                           (address-of error)))
+    (return (* (surface-area-of-ball n)
+               result))))
+
+(pre-let* (INTEGRAL-MINIMUM-SAMPLES 100)
+  (define (integral integrand extent-oracle true-integral r dimension rtol stats)
+    (double integrand-t extent-oracle-t double (const gsl-rng*) (unsigned int) double gsl-rstat-workspace*)
+    (declare integral double
+             error double)
+    (with-vector x dimension
+      (do-while (or (> error rtol)
+                    (> (rerror integral true-integral) rtol)
+                    (< (gsl-rstat-n stats) INTEGRAL-MINIMUM-SAMPLES))
+        ;; TODO: If this neval is not used anywhere, remove it.
+        (let* ((neval int 0))
+          (random-direction-vector r x)
+          (gsl-rstat-add (integral-per-direction integrand x r dimension
+                                                 (extent-oracle x) rtol (address-of neval))
+                         stats))
+        (set integral (gsl-rstat-mean stats))
+        (set error (/ (* CONFIDENCE-INTERVAL-FACTOR (gsl-rstat-sd-mean stats))
+                      integral))))
+    (cond
+     ((> error rtol)
+      (GSL-ERROR-VAL "Integral failed to converge" GSL-ETOL integral)))
+    (return integral)))
+
+(define (volume-cone extent-oracle r mean omega-min omega-max number-of-samples variance)
+  (double extent-oracle-t (const gsl-rng*) (const gsl-vector*) double double (unsigned int) double*)
+  (declare volume double)
+  (let* ((dimension (unsigned int) (: mean size))
+         (vn double (ln-volume-of-ball dimension))
+         (theta-min double (solid-angle->planar-angle (* omega-min (surface-area-of-ball dimension)) dimension))
+         (theta-max double (solid-angle->planar-angle (* omega-max (surface-area-of-ball dimension)) dimension)))
+    (with-rstats stats
+      (with-vector x dimension
+        (for-i i number-of-samples
+          (hollow-cone-random-vector r mean theta-min theta-max x)
+          (gsl-rstat-add (exp (+ vn (* dimension (log (extent-oracle x)))))
+                         stats)))
+      (cond
+       (variance (set (pointer-get variance) (gsl-rstat-variance stats))))
+      (set volume (* (gsl-rstat-mean stats)
+                     (- omega-max omega-min)))))
+  (return volume))
+
+(define (volume-experiment extent-oracle r mean samples-per-cone solid-angle-factor
+                           solid-angle-threshold-exponent-factor number-of-samples)
+  (double extent-oracle-t (const gsl-rng*) (const gsl-vector*)
+          (unsigned int) double double (unsigned int*))
+  (define
+    volume double 0
+    omega-max double 0.5
+    N int 0)
+  (let* ((dimension int (: mean size)))
+    (do-while (> omega-max (pow 2 (- (* dimension solid-angle-threshold-exponent-factor))))
+      (let* ((omega-min double (/ omega-max solid-angle-factor)))
+        (set+ volume (volume-cone extent-oracle r mean omega-min omega-max samples-per-cone NULL))
+        (set+ N samples-per-cone)
+        (set omega-max omega-min))))
+  (set+ volume (volume-cone extent-oracle r mean 0 omega-max samples-per-cone NULL))
+  (set+ N samples-per-cone)
+  (cond
+   (number-of-samples (set (pointer-get number-of-samples) N)))
+  (return volume))
diff --git a/src/gaussian-nd-random.c b/src/gaussian-nd-random.c
deleted file mode 100644
index 16abc94..0000000
--- a/src/gaussian-nd-random.c
+++ /dev/null
@@ -1,51 +0,0 @@
-#include <gsl/gsl_blas.h>
-#include <gsl/gsl_randist.h>
-#include <gsl/gsl_sf_gamma.h>
-#include "nd-random.h"
-#include "utils.h"
-
-#define INTEGRATION_EPSABS 1e-3
-#define INTEGRATION_EPSREL 1e-3
-
-unsigned int shifted_gaussian_random_vector
-(const gsl_rng* r, const gsl_vector* mean, double theta_max, double standard_deviation, gsl_vector* x)
-{
-  for (int i=0; i<x->size; i++)
-    gsl_vector_set(x, i, gsl_ran_gaussian(r, standard_deviation));
-  gsl_vector_add(x, mean);
-  gsl_blas_dscal(1/gsl_blas_dnrm2(x), x);
-  if (angle_between_vectors(x, mean) > theta_max)
-    return 1 + shifted_gaussian_random_vector(r, mean, theta_max, standard_deviation, x);
-  else return 1;
-}
-
-double shifted_gaussian_pdf (double theta, double mean, double theta_max,
-                             double standard_deviation, unsigned int dimension, gsl_integration_workspace* w)
-{
-  // NOTE: Due to the rejection sampling, the returned result is only
-  // proportional to the actual PDF.
-  if (theta > theta_max) return 0;
-  
-  double projection = mean*cos(theta);
-  double integrand (double r, void* params)
-  {
-    return gsl_pow_uint(r, dimension - 1)
-      * gaussian_pdf((r - projection) / standard_deviation);
-  }
-  
-  gsl_function gsl_integrand = {&integrand};
-  double result, error;
-  /* gsl_integration_qagiu(&gsl_integrand, 0, INTEGRATION_EPSABS, INTEGRATION_EPSREL, */
-  /*                       w->limit, w, &result, &error); */
-  double rmax = 0.5 * (projection + sqrt(gsl_pow_2(projection) + 4*(dimension-1)*gsl_pow_2(standard_deviation)));
-  double halfwidth = 3 * standard_deviation;
-  gsl_integration_qag(&gsl_integrand,
-                      rmax - halfwidth < 0 ? 0 : rmax - halfwidth,
-                      rmax + halfwidth,
-                      INTEGRATION_EPSABS, INTEGRATION_EPSREL,
-                      w->limit, GSL_INTEG_GAUSS41, w, &result, &error);
-  return result
-    * gaussian_pdf(sqrt(gsl_pow_2(mean) - gsl_pow_2(projection)) / standard_deviation)
-    / gsl_pow_uint(M_SQRT2*sqrt(M_PI), dimension - 2)
-    / gsl_pow_uint(standard_deviation, dimension);
-}
diff --git a/src/gaussian-nd-random.sc b/src/gaussian-nd-random.sc
new file mode 100644
index 0000000..bfde566
--- /dev/null
+++ b/src/gaussian-nd-random.sc
@@ -0,0 +1,72 @@
+(sc-include "macros/macros")
+
+(pre-include "gsl/gsl_blas.h")
+(pre-include "gsl/gsl_integration.h")
+(pre-include "gsl/gsl_randist.h")
+(pre-include "utils.h")
+
+(pre-define INTEGRATION-EPSABS 1e-3)
+(pre-define INTEGRATION-EPSREL 1e-3)
+
+(define (shifted-gaussian-random-vector r mean theta-max standard-deviation x)
+  ((unsigned int) (const gsl-rng*) (const gsl-vector*) double double gsl-vector*)
+  "Return a random vector generated by the shifted Gaussian method."
+  (for-i i (: x size)
+         (gsl-vector-set x i (gsl-ran-gaussian r standard-deviation)))
+  (gsl-vector-add x mean)
+  (gsl-blas-dscal (/ 1 (gsl-blas-dnrm2 x)) x)
+  (cond
+   ((> (angle-between-vectors x mean) theta-max)
+    (return (+ 1 (shifted-gaussian-random-vector r mean theta-max standard-deviation x))))
+   (else (return 1))))
+
+(define (shifted-gaussian-pdf theta mean theta-max standard-deviation dimension w)
+  (double double double double double (unsigned int) gsl-integration-workspace*)
+  "Return the unnormalized probability density function at THETA for a
+random vector generated using the shifted Gaussian method."
+  (cond
+   ((> theta theta-max) (return 0))
+   (else
+    (define projection double (* mean (cos theta)))
+    (define (integrand r params) (double double void*)
+      (return (* (gsl-pow-uint r (- dimension 1))
+                 (gaussian-pdf (/ (- r projection)
+                                  standard-deviation)))))
+
+    (declare gsl-integrand (struct-variable gsl-function (function (address-of integrand)))
+             result double
+             error double)
+    
+    ;; (gsl-integration-qagiu (address-of gsl-integrand)
+    ;;                        0
+    ;;                        INTEGRATION-EPSABS
+    ;;                        INTEGRATION-EPSREL
+    ;;                        (: w limit)
+    ;;                        w
+    ;;                        (address-of result)
+    ;;                        (address-of error))
+    (let* ((rmax double (/ (+ projection
+                              (sqrt (+ (gsl-pow-2 projection)
+                                       (* 4
+                                          (- dimension 1)
+                                          (gsl-pow-2 standard-deviation)))))
+                           2))
+           (half-width double (* 3 standard-deviation)))
+      (gsl-integration-qag (address-of gsl-integrand)
+                           (if* (< (- rmax half-width) 0)
+                                0
+                                (- rmax half-width))
+                           (+ rmax half-width)
+                           INTEGRATION-EPSABS
+                           INTEGRATION-EPSREL
+                           (: w limit)
+                           GSL-INTEG-GAUSS41
+                           w
+                           (address-of result)
+                           (address-of error)))
+    (return (/ (* result
+                  (gaussian-pdf (/ (sqrt (- (gsl-pow-2 mean)
+                                            (gsl-pow-2 projection)))
+                                   standard-deviation)))
+               (gsl-pow-uint (sqrt (* 2 M-PI)) (- dimension 2))
+               (gsl-pow-uint standard-deviation dimension))))))
diff --git a/src/macros/macros.sc b/src/macros/macros.sc
new file mode 100644
index 0000000..6ffc7ec
--- /dev/null
+++ b/src/macros/macros.sc
@@ -0,0 +1,12 @@
+(sc-define-syntax (for-i index limit body ...)
+  (for ((define index int 0)
+        (< index limit)
+        (set+ index 1))
+    body ...))
+
+(sc-define-syntax (with-alloc var type allocate free body ...)
+  (let* ()
+    (declare var type)
+    (set var allocate)
+    body ...
+    (free var)))
diff --git a/src/nd-random.c b/src/nd-random.c
deleted file mode 100644
index 25abc7a..0000000
--- a/src/nd-random.c
+++ /dev/null
@@ -1,101 +0,0 @@
-#include <float.h>
-#include <math.h>
-#include <gsl/gsl_blas.h>
-#include <gsl/gsl_randist.h>
-#include <gsl/gsl_sf_gamma.h>
-#include "nd-random.h"
-#include "utils.h"
-
-static double beta_inc_unnormalized (double a, double b, double x);
-static double incomplete_wallis_integral (double theta, unsigned int m);
-
-void random_direction_vector (const gsl_rng* r, gsl_vector* x)
-{
-  gsl_ran_dir_nd(r, x->size, x->data);
-}
-
-static void rotate_from_nth_canonical (gsl_vector* x, const gsl_vector* orient)
-{
-  const size_t n = x->size;
-  double xn = gsl_vector_get(x, n - 1);
-  double mun = gsl_vector_get(orient, n - 1);
-  gsl_vector_const_view orient_sub = gsl_vector_const_subvector(orient, 0, n - 1);
-  double b = gsl_blas_dnrm2(&orient_sub.vector);
-  double a = (dot_product(orient, x) - xn*mun) / b;
-  double c = mun, s = sqrt(1 - gsl_pow_2(c));
-  gsl_blas_daxpy((xn*s + a*(c - 1))/b, orient, x);
-  gsl_vector_set(x, n - 1, gsl_vector_get(x, n - 1)
-                 + xn*(c - 1) - a*s
-                 - mun*(xn*s + a*(c - 1))/b);
-}
-
-/* TODO: There is an edge case when mean is the (n-1)th canonical
-   basis vector. Fix it.
-*/
-void hollow_cone_random_vector (const gsl_rng* r, const gsl_vector* mean, double theta_min, double theta_max, gsl_vector* x)
-{
-  unsigned int n = x->size;
-  gsl_ran_dir_nd(r, n - 1, x->data);
-
-  // Generate random vector around the nth canonical basis vector
-  double omega_min = planar_angle_to_solid_angle(theta_min, n);
-  double omega_max = planar_angle_to_solid_angle(theta_max, n);
-  gsl_vector_scale(x, cos(theta_max) * tan(solid_angle_to_planar_angle(gsl_ran_flat(r, omega_min, omega_max), n)));
-  gsl_vector_set(x, n - 1, cos(theta_max));
-  gsl_vector_scale(x, 1.0/gsl_blas_dnrm2(x));
-
-  // Rotate to arbitrary basis
-  rotate_from_nth_canonical(x, mean);
-}
-
-void subsampling_random_vector (const gsl_rng* r, const gsl_vector* mean, double theta_max, gsl_vector* x)
-{
-  hollow_cone_random_vector(r, mean, 0, theta_max, x);
-}
-
-static double beta_inc_unnormalized (double a, double b, double x)
-{
-  return gsl_sf_beta_inc(a, b, x) * gsl_sf_beta(a, b);
-}
-
-static double incomplete_wallis_integral (double theta, unsigned int m)
-{
-  /**
-     @param theta 0 < theta < pi
-     @param m
-     @return \int_0^\theta \sin^m x dx
-   **/
-  if ((theta < 0) || (theta > M_PI))
-    GSL_ERROR("Incomplete Wallis integral only allows theta in [0,pi]", GSL_EDOM);
-  if (theta < M_PI_2)
-    return 0.5 * beta_inc_unnormalized((m+1)/2.0, 0.5, gsl_pow_2(sin(theta)));
-  else
-    return 0.5 * (gsl_sf_beta((m+1)/2.0, 0.5) + beta_inc_unnormalized(0.5, (m+1)/2.0, gsl_pow_2(cos(theta))));
-}
-
-double planar_angle_to_solid_angle (double planar_angle, unsigned int dimension)
-{
-  return 2 * pow(M_PI, (dimension - 1)/2.0)
-    * incomplete_wallis_integral(planar_angle, dimension - 2)
-    / gsl_sf_gamma((dimension - 1)/2.0);
-}
-
-double solid_angle_to_planar_angle (double solid_angle, unsigned int dimension)
-{
-  double f (double planar_angle, void* params)
-  {
-    return planar_angle_to_solid_angle(planar_angle, dimension) - solid_angle;
-  }
-
-  gsl_function gsl_f = {&f};
-  /* if (fabs(GSL_FN_EVAL(&gsl_f, 0))/surface_area_of_ball(dimension) < BISECTION_FUNCTION_EPSABS) */
-  /*   return 0; */
-  /* else if (fabs(GSL_FN_EVAL(&gsl_f, M_PI))/surface_area_of_ball(dimension) < BISECTION_FUNCTION_EPSABS) */
-  /*   return M_PI; */
-  /* else return bisection(&gsl_f, 0, M_PI); */
-  if (solid_angle == 0)
-    return 0;
-  else if (solid_angle == surface_area_of_ball(dimension))
-    return M_PI;
-  else return bisection(&gsl_f, 0, M_PI);
-}
diff --git a/src/nd-random.sc b/src/nd-random.sc
new file mode 100644
index 0000000..3d96113
--- /dev/null
+++ b/src/nd-random.sc
@@ -0,0 +1,98 @@
+(pre-include "math.h")
+(pre-include "gsl/gsl_blas.h")
+(pre-include "gsl/gsl_randist.h")
+(pre-include "gsl/gsl_sf_gamma.h")
+(pre-include "gsl/gsl_vector.h")
+(pre-include "utils.h")
+
+(define (random-direction-vector r x)
+  (void (const gsl-rng*) gsl-vector*)
+  "Generate a random vector distributed uniformly on the unit
+sphere. Write the result to X."
+  (gsl-ran-dir-nd r (: x size) (: x data)))
+
+(define (rotate-from-nth-canonical x orient)
+  ((static void) gsl-vector* (const gsl-vector*))
+  (let* ((n size-t (: x size))
+         (xn double (gsl-vector-get x (- n 1)))
+         (mun double (gsl-vector-get orient (- n 1)))
+         (orient-sub gsl-vector-const-view
+                     (gsl-vector-const-subvector orient 0 (- n 1)))
+         (b double (gsl-blas-dnrm2 (address-of (struct-get orient-sub vector))))
+         (a double (/ (- (dot-product orient x)
+                         (* xn mun))
+                      b))
+         (s double (sqrt (- 1 (gsl-pow-2 mun)))))
+    (gsl-blas-daxpy (/ (+ (* xn s)
+                          (* a (- mun 1)))
+                       b)
+                    orient
+                    x)
+    (gsl-vector-set x
+                    (- n 1)
+                    (+ (gsl-vector-get x (- n 1))
+                       (* xn (- mun 1))
+                       (- (* a s))
+                       (- (/ (* mun (+ (* xn s)
+                                       (* a (- mun 1))))
+                             b))))))
+
+(define (beta-inc-unnormalized a b x) ((static double) double double double)
+  (return (* (gsl-sf-beta-inc a b x)
+             (gsl-sf-beta a b))))
+
+(define (incomplete-wallis-integral theta m) ((static double) double (unsigned int))
+  "Return the incomplete Wallis integral \\int_0^\\theta \\sin^m x
+dx. THETA should be in [0,pi]."
+  (cond
+   ((or (< theta 0) (> theta M-PI))
+    (GSL-ERROR "Incomplete Wallis integral only allows theta in [0,pi]" GSL-EDOM))
+   ((< theta M-PI-2)
+    (return (/ (beta-inc-unnormalized (* 0.5 (+ m 1)) 0.5 (gsl-pow-2 (sin theta)))
+               2)))
+   (else
+    (return (/ (+ (gsl-sf-beta (* 0.5 (+ m 1)) 0.5)
+                  (beta-inc-unnormalized 0.5 (* 0.5 (+ m 1)) (gsl-pow-2 (cos theta))))
+               2)))))
+
+(define (planar-angle->solid-angle planar-angle dimension) ((static double) double (unsigned int))
+  (return (/ (* 2
+                (pow M-PI (* 0.5 (- dimension 1)))
+                (incomplete-wallis-integral planar-angle (- dimension 2)))
+             (gsl-sf-gamma (* 0.5 (- dimension 1))))))
+
+(define (solid-angle->planar-angle solid-angle dimension) ((static double) double (unsigned int))
+  (define (f planar-angle params) (double double void*)
+    (return (- (planar-angle->solid-angle planar-angle dimension)
+               solid-angle)))
+  
+  (declare gsl-f (struct-variable gsl-function (function (address-of f))))
+  ;; TODO: Equality comparisons to floating point values may be
+  ;; problematic. Fix it.
+  (cond
+   ((= solid-angle 0) (return 0))
+   ((= solid-angle (surface-area-of-ball dimension)) (return M-PI))
+   (else (return (bisection (address-of gsl-f) 0 M-PI)))))
+
+;; TODO: There is an edge case when mean is the (n-1)th canonical
+;; basis vector. Fix it.
+(define (hollow-cone-random-vector r mean theta-min theta-max x)
+  (void (const gsl-rng*) (const gsl-vector*) double double gsl-vector*)
+  ;; Generate random vector around the nth canonical basis vector.
+  (let* ((n size-t (: x size)))
+    (gsl-ran-dir-nd r (- n 1) (: x data))
+    (gsl-vector-scale x (* (cos theta-max)
+                           (tan (solid-angle->planar-angle
+                                 (gsl-ran-flat r
+                                               (planar-angle->solid-angle theta-min n)
+                                               (planar-angle->solid-angle theta-max n))
+                                 n))))
+    (gsl-vector-set x (- n 1) (cos theta-max)))
+  (gsl-vector-scale x (/ 1 (gsl-blas-dnrm2 x)))
+
+  ;; Rotate to arbitrary basis.
+  (rotate-from-nth-canonical x mean))
+
+(define (subsampling-random-vector r mean theta-max x)
+  (void (const gsl-rng*) (const gsl-vector*) double gsl-vector*)
+  (hollow-cone-random-vector r mean 0 theta-max x))
diff --git a/src/oracles.c b/src/oracles.c
deleted file mode 100644
index 2b440eb..0000000
--- a/src/oracles.c
+++ /dev/null
@@ -1,108 +0,0 @@
-#include <math.h>
-#include <gsl/gsl_math.h>
-#include <gsl/gsl_randist.h>
-#include "oracles.h"
-#include "utils.h"
-
-double bernoulli_extent_generator
-(const gsl_rng* r, double p, double r0, double r1)
-{
-  return gsl_ran_bernoulli(r, p) ? r1 : r0;
-}
-
-double bernoulli_true_volume (double p, double r0, double r1, unsigned int dimension)
-{
-  return volume_of_ball(dimension)
-    * (p * gsl_pow_uint(r1, dimension)
-       + (1 - p) * gsl_pow_uint(r0, dimension));
-}
-
-double uniform_extent_generator (const gsl_rng* r, double a, double b)
-{
-  return gsl_ran_flat(r, a, b);
-}
-
-// TODO: Verify the accuracy of this function for non-trivial a, b
-double uniform_true_volume (double a, double b, unsigned int dimension)
-{
-  return exp(ln_volume_of_ball(dimension) + dimension*log(b) - log(dimension + 1))
-    - exp(ln_volume_of_ball(dimension) + dimension*log(a) - log(dimension + 1));
-}
-
-double beta_extent_generator
-(const gsl_rng* r, double alpha, double beta)
-{
-  return gsl_ran_beta(r, alpha, beta);
-}
-
-double beta_true_volume (double alpha, double beta, unsigned int dimension)
-{
-  double vol = volume_of_ball(dimension);
-  for (unsigned int r=0; r<dimension; r++)
-    vol *= (alpha + r) / (alpha + beta + r);
-  return vol;
-}
-
-static double infinity_norm (const gsl_vector* x)
-{
-  double max = fabs(gsl_vector_get(x, 0));
-  for (int i=1; i<x->size; i++)
-    max = GSL_MAX(max, fabs(gsl_vector_get(x, i)));
-  return max;
-}
-
-double cube_extent_oracle (const gsl_vector* x, double edge)
-{
-  return edge / 2 / infinity_norm(x);
-}
-
-double cube_extent_oracle_with_center
-(const gsl_vector* x, const gsl_vector* center, double edge)
-{
-  double min = (0.5*edge - GSL_SIGN(gsl_vector_get(x, 0))*gsl_vector_get(center, 0))
-    / fabs(gsl_vector_get(x, 0));
-  for (int i=1; i<center->size; i++)
-    min = GSL_MIN(min, (0.5*edge - GSL_SIGN(gsl_vector_get(x, i))*gsl_vector_get(center, i))
-                  / fabs(gsl_vector_get(x, i)));
-  return min;
-}
-
-double cube_true_volume (double edge, unsigned int dimension)
-{
-  return gsl_pow_uint(edge, dimension);
-}
-
-double cube_maximum_extent (double edge, unsigned int dimension)
-{
-  return edge * sqrt(dimension) / 2;
-}
-
-double ellipsoid_extent_oracle (const gsl_vector* x, const gsl_vector* axes)
-{
-  double k = 0;
-  for (int i=0; i<axes->size; i++)
-    k += gsl_pow_2(gsl_vector_get(x, i) / gsl_vector_get(axes, i));
-  return 1/sqrt(k);
-}
-
-double ellipsoid_true_volume (const gsl_vector* axes)
-{
-  unsigned int dimension = axes->size;
-  double vol = volume_of_ball(dimension);
-  for (int i=0; i<dimension; i++)
-    vol *= gsl_vector_get(axes, i);
-  return vol;
-}
-
-double spheroid_extent_oracle (const gsl_vector* x, double eccentricity)
-{
-  unsigned int dimension = x->size;
-  gsl_vector_const_view xsub = gsl_vector_const_subvector(x, 1, dimension - 1);
-  return 1/sqrt(gsl_pow_2(gsl_blas_dnrm2(&xsub.vector))
-                + gsl_pow_2(gsl_vector_get(x, 0) / eccentricity));
-}
-
-double spheroid_true_volume (double eccentricity, unsigned int dimension)
-{
-  return volume_of_ball(dimension) * eccentricity;
-}
diff --git a/src/oracles.sc b/src/oracles.sc
new file mode 100644
index 0000000..9680479
--- /dev/null
+++ b/src/oracles.sc
@@ -0,0 +1,90 @@
+(sc-include "macros/macros")
+
+(pre-include "math.h")
+(pre-include "gsl/gsl_blas.h")
+(pre-include "gsl/gsl_randist.h")
+(pre-include "utils.h")
+
+(define (bernoulli-extent-generator r p r0 r1) (double (const gsl-rng*) double double double)
+  (return (if* (gsl-ran-bernoulli r p) r1 r0)))
+
+(define (bernoulli-true-volume p r0 r1 dimension) (double double double double (unsigned int))
+  (return (* (volume-of-ball dimension)
+             (+ (* p (gsl-pow-uint r1 dimension))
+                (* (- 1 p) (gsl-pow-uint r0 dimension))))))
+
+(define (uniform-extent-generator r a b) (double (const gsl-rng*) double double)
+  (return (gsl-ran-flat r a b)))
+
+;; TODO: Verify the accuracy of this function for non-trivial a, b.
+(define (uniform-true-volume a b dimension) (double double double (unsigned int))
+  (return (- (exp (+ (ln-volume-of-ball dimension)
+                     (* dimension (log b))
+                     (- (log (+ dimension 1)))))
+             (exp (+ (ln-volume-of-ball dimension)
+                     (* dimension (log a))
+                     (- (log (+ dimension 1))))))))
+
+(define (beta-extent-generator r alpha beta) (double (const gsl-rng*) double double)
+  (return (gsl-ran-beta r alpha beta)))
+
+(define (beta-true-volume alpha beta dimension) (double double double (unsigned int))
+  (let* ((vol double (volume-of-ball dimension)))
+    (for-i r dimension
+      (set* vol (/ (+ alpha r)
+                   (+ alpha beta r))))
+    (return vol)))
+
+(define (infinity-norm x) ((static double) (const gsl-vector*))
+  (let* ((max double (fabs (gsl-vector-get x 0))))
+    ;; TODO: Start this loop from i = 1, not i = 0. That would be
+    ;; slightly faster.
+    (for-i i (: x size)
+      (set max (GSL-MAX max (fabs (gsl-vector-get x i)))))
+    (return max)))
+
+(define (cube-extent-oracle x edge) (double (const gsl-vector*) double)
+  (return (/ edge 2 (infinity-norm x))))
+
+(sc-define-syntax (compute-cube-extent-oracle-minimizand i)
+  (/ (- (/ edge 2)
+        (* (GSL-SIGN (gsl-vector-get x i))
+           (gsl-vector-get center i)))
+     (fabs (gsl-vector-get x i))))
+
+(define (cube-extent-oracle-with-center x center edge) (double (const gsl-vector*) (const gsl-vector*) double)
+  (let* ((min double (compute-cube-extent-oracle-minimizand 0)))
+    ;; TODO: Start this loop from i = 1, not i = 0. That would be
+    ;; slightly faster.
+    (for-i i (: center size)
+      (set min (GSL-MIN min (compute-cube-extent-oracle-minimizand i))))
+    (return min)))
+
+(define (cube-true-volume edge dimension) (double double (unsigned int))
+  (return (gsl-pow-uint edge dimension)))
+
+(define (cube-maximum-extent edge dimension) (double double (unsigned int))
+  (return (/ (* edge (sqrt dimension))
+             2)))
+
+(define (ellipsoid-extent-oracle x axes) (double (const gsl-vector*) (const gsl-vector*))
+  (let* ((k double 0))
+    (for-i i (: axes size)
+      (set+ k (gsl-pow-2 (/ (gsl-vector-get x i)
+                            (gsl-vector-get axes i)))))
+    (return (/ (sqrt k)))))
+
+(define (ellipsoid-true-volume axes) (double (const gsl-vector*))
+  (let* ((vol double (volume-of-ball (: axes size))))
+    (for-i i (: axes size)
+      (set* vol (gsl-vector-get axes i)))
+    (return vol)))
+
+(define (spheroid-extent-oracle x eccentricity) (double (const gsl-vector*) double)
+  (let* ((xsub gsl-vector-const-view
+               (gsl-vector-const-subvector x 1 (- (: x size) 1))))
+    (return (/ (sqrt (+ (gsl-pow-2 (gsl-blas-dnrm2 (address-of (struct-get xsub vector))))
+                        (gsl-pow-2 (/ (gsl-vector-get x 0) eccentricity))))))))
+
+(define (spheroid-true-volume eccentricity dimension) (double double (unsigned int))
+  (return (* (volume-of-ball dimension) eccentricity)))
diff --git a/src/utils.c b/src/utils.c
deleted file mode 100644
index b22e07b..0000000
--- a/src/utils.c
+++ /dev/null
@@ -1,160 +0,0 @@
-#include <gsl/gsl_blas.h>
-#include <gsl/gsl_math.h>
-#include <gsl/gsl_randist.h>
-#include <gsl/gsl_sf.h>
-#include <gsl/gsl_sf_gamma.h>
-#include "utils.h"
-
-// This function exists so that guile code can access the value of
-// M_PI as provided by <gsl_math.h>.
-double pi ()
-{
-  /**
-   @return Value of pi
-  **/
-  return M_PI;
-}
-
-double volume_of_ball (unsigned int dimension)
-{
-  /**
-     @param dimension Dimension of the ball (2 for a circle, 3 for a sphere, etc.)
-     @return Volume of the unit ball
-   **/
-  return exp(ln_volume_of_ball(dimension));
-}
-
-double ln_volume_of_ball (unsigned int dimension)
-{
-  /**
-     @param dimension Dimension of the ball (2 for a circle, 3 for a sphere, etc.)
-     @return Natural logarithm of the volume of the unit ball
-   **/
-  return (dimension/2.0)*M_LNPI - gsl_sf_lngamma(1 + dimension/2.0);
-}
-
-double surface_area_of_ball (unsigned int dimension)
-{
-  /**
-     @param dimension Dimension of the ball (2 for a circle, 3 for a sphere, etc.)
-     @return Surface area of the unit ball
-   **/
-  return dimension * volume_of_ball(dimension);
-}
-
-double ln_surface_area_of_ball (unsigned int dimension)
-{
-  /**
-     @param dimension Dimension of the ball (2 for a circle, 3 for a sphere, etc.)
-     @return Natural logarithm of the surface area of the unit ball
-   **/
-  return log(dimension) + ln_volume_of_ball(dimension);
-}
-
-double lower_incomplete_gamma (double s, double x)
-{
-  return gsl_sf_gamma(s) * gsl_sf_gamma_inc_P(s, x);
-}
-
-double angle_between_vectors (const gsl_vector* x, const gsl_vector* y)
-{
-  /**
-     @param x Vector 1
-     @param y Vector 2
-     @return Angle between the two vectors in the range [0, pi]
-  **/
-  double dot_product;
-  gsl_blas_ddot(x, y, &dot_product);
-  // TODO: Is this a valid floating point comparison?
-  if (dot_product != 0)
-    return acos(dot_product / gsl_blas_dnrm2(x) / gsl_blas_dnrm2(y));
-  else return 0;
-}
-
-double dot_product (const gsl_vector* x, const gsl_vector* y)
-{
-  /**
-     @param x Vector 1
-     @param y Vector 2
-     @return Dot product between the two vectors
-  **/
-  double result;
-  gsl_blas_ddot(x, y, &result);
-  return result;
-}
-
-double gaussian_pdf (double x)
-{
-  /**
-     @param x
-     @return exp(-x^2/2) / sqrt(2*pi)
-   **/
-  return gsl_ran_gaussian_pdf(x, 1.0);
-}
-
-double gaussian_cdf (double x)
-{
-  /**
-     @param x
-     @return int_{-inf}^x gaussian_pdf(t) dt
-   **/
-  return 0.5*(1 + gsl_sf_erf(x/M_SQRT2));
-}
-
-double rerror (double approx, double exact)
-{
-  /**
-     @param approx Approximate value
-     @param exact Exact value
-     @return Relative error
-  **/
-  return fabs(1 - approx/exact);
-}
-
-#define BISECTION_EPSABS 0
-#define BISECTION_EPSREL 1e-6
-#define BISECTION_MAX_ITERATIONS 1000
-#define BISECTION_FUNCTION_EPSABS 1e-9
-double bisection (gsl_function* f, double a, double b)
-{
-  const gsl_root_fsolver_type* T = gsl_root_fsolver_bisection;
-  gsl_root_fsolver* s = gsl_root_fsolver_alloc(T);
-  double solution;
-  int status;
-
-  void error_handler (const char* reason, const char* file, int line, int gsl_errno)
-  {
-    fprintf(stderr, "Bisection error handler invoked.\n");
-    fprintf(stderr, "f(%g) = %g\n", a, GSL_FN_EVAL(f, a));
-    fprintf(stderr, "f(%g) = %g\n", b, GSL_FN_EVAL(f, b));
-    fprintf(stderr, "gsl: %s:%d: ERROR: %s\n", file, line, reason);
-    abort();
-  }
-
-  gsl_error_handler_t* old_handler = gsl_set_error_handler(error_handler);
-  gsl_root_fsolver_set(s, f, a, b);
-  do {
-    status = gsl_root_fsolver_iterate (s);
-    solution = gsl_root_fsolver_root (s);
-    a = gsl_root_fsolver_x_lower (s);
-    b = gsl_root_fsolver_x_upper (s);
-    status = gsl_root_test_interval (a, b, BISECTION_EPSABS, BISECTION_EPSREL);
-  } while (status == GSL_CONTINUE);
-  gsl_root_fsolver_free(s);
-  gsl_set_error_handler(old_handler);
-  return solution;
-}
-
-double bisection_rlimit (gsl_function* f, double a, double b)
-{
-  int sign = SIGNUM(GSL_FN_EVAL(f, a));
-  for (int i=0; i<BISECTION_MAX_ITERATIONS; i++)
-    if (sign*GSL_FN_EVAL(f, b) > 0) b *= 2;
-    else return b;
-  fprintf(stderr, "Bisection bracketing failed\n");
-  exit(1);
-}
-#undef BISECTION_EPSABS
-#undef BISECTION_EPSREL
-#undef BISECTION_MAX_ITERATIONS
-#undef BISECTION_FUNCTION_EPSABS
diff --git a/src/utils.sc b/src/utils.sc
new file mode 100644
index 0000000..751cd77
--- /dev/null
+++ b/src/utils.sc
@@ -0,0 +1,123 @@
+(sc-include "macros/macros")
+
+(pre-include "gsl/gsl_math.h")
+(pre-include "gsl/gsl_blas.h")
+(pre-include "gsl/gsl_randist.h")
+(pre-include "gsl/gsl_roots.h")
+(pre-include "gsl/gsl_sf.h")
+
+(pre-define (SIGNUM x)
+  (if* (< x 0) -1 1))
+
+;; This function exists so that guile code can access the value of
+;; M_PI as provided by <gsl_math.h>.
+(define (pi) (double)
+  "Return the value of pi."
+  (return M-PI))
+
+(define (ln-volume-of-ball dimension) (double (unsigned int))
+  "Return the natural logarithm of the volume of the unit
+ball. DIMENSION of 2 corresponds to a circle, 3 corresponds to a
+sphere, etc."
+  (return (- (* 0.5 dimension M-LNPI)
+             (gsl-sf-lngamma (+ 1 (* 0.5 dimension))))))
+
+(define (volume-of-ball dimension) (double (unsigned int))
+  "Return the volume of the unit ball of DIMENSION. DIMENSION of 2
+corresponds to a circle, 3 corresponds to a sphere, etc."
+  (return (exp (ln-volume-of-ball dimension))))
+
+(define (ln-surface-area-of-ball dimension) (double (unsigned int))
+  "Return the natural logarithm of the surface area of the unit
+ball. DIMENSION of 2 corresponds to a circle, 3 corresponds to a
+sphere, etc."
+  (return (+ (log dimension)
+             (ln-volume-of-ball dimension))))
+
+(define (surface-area-of-ball dimension) (double (unsigned int))
+  "Return the surface area of the unit ball of DIMENSION. DIMENSION of
+2 corresponds to a circle, 3 corresponds to a sphere, etc."
+  (return (* dimension (volume-of-ball dimension))))
+
+(define (lower-incomplete-gamma s x) (double double double)
+  (return (* (gsl-sf-gamma s)
+             (gsl-sf-gamma-inc-P s x))))
+
+(define (angle-between-vectors x y) (double (const gsl-vector*) (const gsl-vector*))
+  "Return the angle between vectors X and Y. The returned value is in
+the range [0,pi]."
+  (declare dot-product double)
+  (gsl-blas-ddot x y (address-of dot-product))
+  ;; TODO: Is this a valid floating point comparison?
+  (return (if* (= dot_product 0)
+               0
+               (acos (/ dot-product
+                        (gsl-blas-dnrm2 x)
+                        (gsl-blas-dnrm2 y))))))
+
+(define (dot-product x y) (double (const gsl-vector*) (const gsl-vector*))
+  "Return the dot product of vectors X and Y."
+  (declare result double)
+  (gsl-blas-ddot x y (address-of result))
+  (return result))
+
+(define (gaussian-pdf x) (double double)
+  "Return exp(-x^2/2) / sqrt(2*pi)"
+  (return (gsl-ran-gaussian-pdf x 1)))
+
+(define (gaussian-cdf x) (double double)
+  "Return \\int_{-\\inf}^x gaussian_pdf(t) dt."
+  (return (* 0.5 (+ 1 (gsl-sf-erf (/ x M-SQRT2))))))
+
+(define (rerror approx exact) (double double double)
+  "Return the relative error between approximate value APPROX and
+exact value EXACT."
+  (return (fabs (- 1 (/ approx exact)))))
+
+(sc-define-syntax* (with-root-fsolver solver solver-type function a b body ...)
+  (let ((solver-type (sc-gensym)))
+    `(begin
+       (let* ((,solver-type (const gsl-root-fsolver-type*) ,solver-type))
+         (with-alloc solver gsl-root-fsolver*
+                     (gsl-root-fsolver-alloc ,solver-type)
+                     gsl-root-fsolver-free
+                     (gsl-root-fsolver-set solver ,function ,a ,b)
+                     ,@body)))))
+
+(sc-define-syntax* (with-error-handler handler body ...)
+  (let ((old-handler (sc-gensym)))
+    `(begin
+       (let* ((,old-handler gsl-error-handler-t* (gsl-set-error-handler ,handler)))
+         ,@body
+         (gsl-set-error-handler ,old-handler)))))
+
+(pre-let* (BISECTION-EPSABS 0 BISECTION-EPSREL 1e-6 BISECTION-MAX-ITERATIONS 1000)
+  (define (bisection f a b) (double gsl-function* double double)
+    (declare solution double)
+    (define (error-handler reason file line gsl-errno) (void (const char*) (const char*) int int)
+      (fprintf stderr "Bisection error handler invoked.\n")
+      (fprintf stderr "f(%g) = %g\n" a (GSL_FN_EVAL f a))
+      (fprintf stderr "f(%g) = %g\n" b (GSL_FN_EVAL f b))
+      (fprintf stderr "gsl: %s:%d: ERROR: %s\n" file line reason)
+      (abort))
+
+    (with-root-fsolver solver gsl-root-fsolver-bisection f a b
+      (with-error-handler error-handler
+        (do-while (= (gsl-root-test-interval (gsl-root-fsolver-x-lower solver)
+                                             (gsl-root-fsolver-x-upper solver)
+                                             BISECTION-EPSABS
+                                             BISECTION-EPSREL)
+                     GSL-CONTINUE)
+          (gsl-root-fsolver-iterate solver)))
+      (set solution (gsl-root-fsolver-root solver)))
+    (return solution))
+
+  (define (bisection-rlimit f a b) (double gsl-function* double double)
+    (let* ((sign int (SIGNUM (GSL_FN_EVAL f a))))
+      (for-i i BISECTION-MAX-ITERATIONS
+        (cond
+         ((> (* sign (GSL_FN_EVAL f b)) 0)
+          (set* b 2))
+         (else (return b)))))
+    (fprintf stderr "Bisection bracketing failed\n")
+    (abort)))