From c31a57815be4d3ef5f6fd28cc9a34b3be55bfe56 Mon Sep 17 00:00:00 2001 From: Arun Isaac Date: Fri, 5 Feb 2021 15:00:27 +0530 Subject: 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. --- src/utils.c | 160 ------------------------------------------------------------ 1 file changed, 160 deletions(-) delete mode 100644 src/utils.c (limited to 'src/utils.c') 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 -#include -#include -#include -#include -#include "utils.h" - -// This function exists so that guile code can access the value of -// M_PI as provided by . -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 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 -- cgit v1.2.3