move the math functions to a separate file; remove dependency on GSL due to license

src/black_scholes.c

@@ -1,15 +1,11 @@
 #define _XOPEN_SOURCE
 
 #include "black_scholes.h"
+#include "utils.h"
 #include <cblas.h>
 #include <math.h>
 #include <time.h>
 
-double normalcdf(double value)
-{
-	return 0.5 * erfc(-value * M_SQRT1_2);
-}
-
 double bsm(double spot, double rfr, double vol, double strike,
 		struct tm expiry, struct tm value, int type)
 {

src/gbm_mc.c

@@ -1,12 +1,10 @@
 #define _XOPEN_SOURCE
 #define max(X,Y) (((X) > (Y)) ? (X) : (Y))
 
-#include <gsl/gsl_statistics_double.h>
-#include <gsl/gsl_rng.h>
-#include <gsl/gsl_randist.h>
+#include "gbm_mc.h"
+#include "utils.h"
 #include <math.h>
 #include <time.h>
-#include "gbm_mc.h"
 
 double gbm_simulation(double spot, double rfr, double vol, double tte, double rand)
 {
@@ -21,32 +19,23 @@ double gbm_simulation(double spot, double rfr, double vol, double tte, double ra
 double gbm(double spot, double rfr, double vol, double strike, struct tm expiry,
 		struct tm value, int type, int sims)
 {
-	const gsl_rng_type *T;
 	double tte, expiry_date, value_date, level, price, rand;
-	double *results;
-	gsl_rng *r;
+	double results = 0;
 	int i;
 
-	results = (double *)malloc(sizeof(double) * sims);
+	if (sims < 1) sims = 1;
 
 	expiry_date = mktime(&expiry);
 	value_date = mktime(&value);
 	tte = difftime(expiry_date, value_date) / (60 * 60 * 24 * 365);
 
-	/* GSL RNG setup */
-	gsl_rng_env_setup();
-	T = gsl_rng_default;
-	r = gsl_rng_alloc(T);
-
 	for (i=0; i<sims; i++) {
-		rand = gsl_ran_gaussian(r, 1);
+		rand = gaussrand();
 		level = gbm_simulation(spot, rfr, vol, tte, rand);
-		results[i] = max((level - strike) * type, 0);
+		results += max((level - strike) * type, 0);
 	}
 
-	price = gsl_stats_mean(results, 1, sims) * 1 / pow((1 + rfr), tte);
+	price = results / sims * 1 / pow((1 + rfr), tte);
 
-	gsl_rng_free(r);
-	free(results);
 	return price;
 }

src/gbm_mc.h

@@ -1,6 +1,8 @@
 #ifndef GBM_MC
 #define GBM_MC
 
+#include <time.h>
+
 double gbm_simulation(double spot, double rfr, double vol, double tte, double rand);
 
 double gbm(double spot, double rfr, double vol, double strike, struct tm expiry, struct tm value, int type, int sims);

src/utils.c

@@ -0,0 +1,52 @@
+#define M_SQRT1_2  0.707106781186547524401
+
+#include <math.h>
+#include <stdlib.h>
+
+double normalcdf(double z)
+{
+	/* https://www.johndcook.com/blog/cpp_phi/ */
+	double a1, a2, a3, a4, a5, p, t, y;
+	int sign;
+
+	a1 =  0.254829592;
+	a2 = -0.284496736;
+	a3 =  1.421413741;
+	a4 = -1.453152027;
+	a5 =  1.061405429;
+	p  =  0.3275911;
+
+	sign = 1;
+	if (z < 0)
+		sign = -1;
+	z = fabs(z) / pow(2, 0.5);
+	t = 1 / (1 + p * z);
+	y = 1 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * exp(-z * z);
+	return 0.5 * (1 + sign * y);
+}
+
+double gaussrand()
+{
+	static double V1, V2, S;
+	static int phase = 0;
+	double X;
+
+	if (phase == 0) {
+		do {
+			double U1 = (double)rand() / RAND_MAX;
+			double U2 = (double)rand() / RAND_MAX;
+
+			V1 = 2 * U1 - 1;
+			V2 = 2 * U2 - 1;
+			S = V1 * V1 + V2 * V2;
+		} while (S >= 1 || S == 0);
+
+		X = V1 * sqrt(-2 * log(S) / S);
+	} else {
+		X = V2 * sqrt(-2 * log(S) / S);
+	}
+
+	phase = 1 - phase;
+
+	return X;
+}

src/utils.h

@@ -0,0 +1,8 @@
+#ifndef UTILS_H
+#define UTILS_H
+
+double normalcdf(double value);
+
+double gaussrand();
+
+#endif