aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorArun Isaac2021-06-30 14:49:03 +0530
committerArun Isaac2021-06-30 14:49:03 +0530
commit0c7e6308c128a0fab0e484d974ca661b26cb459a (patch)
treef4bfff9334ed9826512f115bf86eabc5a0e61cfd /src
parentec527081eaaec124c070288befc0815b139ee232 (diff)
downloadnsmc-0c7e6308c128a0fab0e484d974ca661b26cb459a.tar.gz
nsmc-0c7e6308c128a0fab0e484d974ca661b26cb459a.tar.lz
nsmc-0c7e6308c128a0fab0e484d974ca661b26cb459a.zip
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.
Diffstat (limited to 'src')
-rw-r--r--src/oracles.sc24
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))))