Add offcenter ellipsoid extent oracle.

* src/oracles.sc: Include gsl/gsl_poly.h.
(ellipsoid-extent-oracle-with-center): New function.
* include/oracles.h (ellipsoid_params): Add center member.
(ellipsoid_extent_oracle_with_center): Declare it.

1 files changed, 24 insertions, 0 deletions

diff --git a/src/oracles.sc b/src/oracles.sc
index 6f0ecbf..689ac9d 100644
--- a/src/oracles.sc
+++ b/src/oracles.sc
@@ -21,6 +21,7 @@
 
 (pre-include "math.h")
 (pre-include "gsl/gsl_blas.h")
+(pre-include "gsl/gsl_poly.h")
 (pre-include "gsl/gsl_randist.h")
 (pre-include "oracles.h")
 (pre-include "utils.h")
@@ -101,6 +102,29 @@
                             (gsl-vector-get (: params axes) i)))))
     (return (/ (sqrt k)))))
 
+(define (ellipsoid-extent-oracle-with-center r x -params) (double (const gsl-rng*) (const gsl-vector*) void*)
+  (let* ((params (const ellipsoid-params*) (convert-type -params ellipsoid-params*))
+         ;; a, b and c are the coefficients of a quadratic equation
+         ;; ax^2 + bx + c = 0.
+         (a double 0)
+         (b double 0)
+         (c double 0))
+    (for-i i (: (: params axes) size)
+      (set+ a (gsl-pow-2 (/ (gsl-vector-get x i)
+                            (gsl-vector-get (: params axes) i))))
+      (set+ b (/ (* (gsl-vector-get x i)
+                    (gsl-vector-get (: params center) i))
+                 (gsl-pow-2 (gsl-vector-get (: params axes) i))))
+      (set+ c (gsl-pow-2 (/ (gsl-vector-get (: params center) i)
+                            (gsl-vector-get (: params axes) i)))))
+    (set* b 2)
+    (set- c 1)
+    (declare k1 double
+             k2 double)
+    (gsl-poly-solve-quadratic a b c (address-of k1) (address-of k2))
+    ;; k2 > k1, k1 < 0, k2 > 0. k2 is our desired root.
+    (return k2)))
+
 (define (ellipsoid-true-volume dimension -params) (double (unsigned int) void*)
   (let* ((params (const ellipsoid-params*) (convert-type -params ellipsoid-params*))
          (vol double (volume-of-ball (: (: params axes) size))))