diff options
| author | Arun Isaac | 2021-06-30 14:49:03 +0530 |
|---|---|---|
| committer | Arun Isaac | 2021-06-30 14:49:03 +0530 |
| commit | 0c7e6308c128a0fab0e484d974ca661b26cb459a (patch) | |
| tree | f4bfff9334ed9826512f115bf86eabc5a0e61cfd /src | |
| parent | ec527081eaaec124c070288befc0815b139ee232 (diff) | |
| download | nsmc-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.sc | 24 |
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)))) |