Sobol series works, need to change makefile to account for 2 functions that are LGPL

Makefile

@@ -1,7 +1,8 @@
 CC=gcc
-CFLAGS=-Wall -fPIC -O3 -ansi -pedantic-errors -pthread
+CFLAGS=-Wall -g -fPIC -O3 -pthread
-LDFLAGS=-lm
+CPPFLAGS=-pedantic -std=c++11
-PREFIX= /usr/local
+LDFLAGS=-lm -lstdc++
+PREFIX=/usr/local
 
 WINDOWS_CC=x86_64-w64-mingw32-gcc
 WINDOWS_CFLAGS= -Wall -O3 -ansi -pedantic-errors -pthread
@@ -9,7 +10,7 @@ WINDOWS_PLATFORM= windows
 
 .DEFAULT_GOAL := build/opt-pricer
 
-build/opt-pricer : src/opt-pricer.c build/depends/linux/gbm.o build/depends/linux/black_scholes.o build/depends/linux/utils.o
+build/opt-pricer : src/opt-pricer.c build/depends/linux/gbm.o build/depends/linux/black_scholes.o build/depends/linux/utils.o build/depends/linux/sobol.o build/depends/linux/asa241.o
 	@$(CC) $(CFLAGS) $^ $(LDFLAGS) -o $@
 
 build/depends/%/gbm.o : src/gbm_mc.c | folders
@@ -24,6 +25,12 @@ build/depends/%/utils.o : src/utils.c | folders
 build/depends/%/strptime.o : src/strptime.c | folders
 	@$(CC) $(CFLAGS) -c $^ $(LDFLAGS) -o $@
 
+build/depends/%/asa241.o : src/asa241.c | folders
+	@$(CC) $(CFLAGS) -c $^ $(LDFLAGS) -o $@
+
+build/depends/%/sobol.o : src/sobol.cpp | folders
+	@$(CC) $(CFLAGS) $(CPPFLAGS) -c $^ $(LDFLAGS) -o $@
+
 folders:
 	@mkdir -p build
 	@mkdir -p build/depends

src/asa241.c

@@ -0,0 +1,579 @@
+# include <stdlib.h>
+# include <stdio.h>
+# include <math.h>
+# include <time.h>
+
+# include "asa241.h"
+
+/******************************************************************************/
+
+void normal_01_cdf_values ( int *n_data, double *x, double *fx )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+    NORMAL_01_CDF_VALUES returns some values of the Normal 01 CDF.
+
+  Discussion:
+
+    In Mathematica, the function can be evaluated by:
+
+      Needs["Statistics`ContinuousDistributions`"]
+      dist = NormalDistribution [ 0, 1 ]
+      CDF [ dist, x ]
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    28 August 2004
+
+  Author:
+
+    John Burkardt
+
+  Reference:
+
+    Milton Abramowitz, Irene Stegun,
+    Handbook of Mathematical Functions,
+    National Bureau of Standards, 1964,
+    ISBN: 0-486-61272-4,
+    LC: QA47.A34.
+
+    Stephen Wolfram,
+    The Mathematica Book,
+    Fourth Edition,
+    Cambridge University Press, 1999,
+    ISBN: 0-521-64314-7,
+    LC: QA76.95.W65.
+
+  Parameters:
+
+    Input/output, int *N_DATA.  The user sets N_DATA to 0 before the
+    first call.  On each call, the routine increments N_DATA by 1, and
+    returns the corresponding data; when there is no more data, the
+    output value of N_DATA will be 0 again.
+
+    Output, double *X, the argument of the function.
+
+    Output, double *FX, the value of the function.
+*/
+{
+# define N_MAX 17
+
+  double fx_vec[N_MAX] = { 
+     0.5000000000000000E+00,  
+     0.5398278372770290E+00,  
+     0.5792597094391030E+00,  
+     0.6179114221889526E+00,  
+     0.6554217416103242E+00,  
+     0.6914624612740131E+00,  
+     0.7257468822499270E+00,  
+     0.7580363477769270E+00,  
+     0.7881446014166033E+00,  
+     0.8159398746532405E+00,  
+     0.8413447460685429E+00,  
+     0.9331927987311419E+00,  
+     0.9772498680518208E+00,  
+     0.9937903346742239E+00,  
+     0.9986501019683699E+00,  
+     0.9997673709209645E+00,  
+     0.9999683287581669E+00 };
+
+  double x_vec[N_MAX] = { 
+     0.0000000000000000E+00,    
+     0.1000000000000000E+00,  
+     0.2000000000000000E+00,  
+     0.3000000000000000E+00,  
+     0.4000000000000000E+00,  
+     0.5000000000000000E+00,  
+     0.6000000000000000E+00,  
+     0.7000000000000000E+00,  
+     0.8000000000000000E+00,  
+     0.9000000000000000E+00,  
+     0.1000000000000000E+01,  
+     0.1500000000000000E+01,  
+     0.2000000000000000E+01,  
+     0.2500000000000000E+01,  
+     0.3000000000000000E+01,  
+     0.3500000000000000E+01,  
+     0.4000000000000000E+01 };
+
+  if ( *n_data < 0 )
+  {
+    *n_data = 0;
+  }
+
+  *n_data = *n_data + 1;
+
+  if ( N_MAX < *n_data )
+  {
+    *n_data = 0;
+    *x = 0.0;
+    *fx = 0.0;
+  }
+  else
+  {
+    *x = x_vec[*n_data-1];
+    *fx = fx_vec[*n_data-1];
+  }
+
+  return;
+# undef N_MAX
+}
+/******************************************************************************/
+
+float r4_huge ( )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+    R4_HUGE returns a "huge" R4.
+
+  Discussion:
+
+    The value returned by this function is NOT required to be the
+    maximum representable R4.  This value varies from machine to machine,
+    from compiler to compiler, and may cause problems when being printed.
+    We simply want a "very large" but non-infinite number.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    14 February 2009
+
+  Author:
+
+    John Burkardt
+
+  Parameters:
+
+    Output, float R4_HUGE, a "huge" R4 value.
+*/
+{
+  float value;
+
+  value = 1.0E+30;
+
+  return value;
+}
+/******************************************************************************/
+
+float r4_normal_01_cdf_inverse ( float p )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+    R4_NORMAL_01_CDF_INVERSE inverts the standard normal CDF.
+
+  Discussion:
+
+    The result is accurate to about 1 part in 10**7.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    27 December 2004
+
+  Author:
+
+    Original FORTRAN77 version by Michael Wichura.
+    C version by John Burkardt.
+
+  Reference:
+
+    Michael Wichura,
+    The Percentage Points of the Normal Distribution,
+    Algorithm AS 241,
+    Applied Statistics,
+    Volume 37, Number 3, pages 477-484, 1988.
+
+  Parameters:
+
+    Input, float P, the value of the cumulative probability densitity function.
+    0 < P < 1.  If P is outside this range, an "infinite" result is returned.
+
+    Output, float R4_NORMAL_01_CDF_INVERSE, the normal deviate value with the 
+    property that the probability of a standard normal deviate being less than or
+    equal to this value is P.
+*/
+{
+  static float a[4] = { 3.3871327179, 50.434271938, 159.29113202, 59.109374720 };
+  static float b[4] = { 1.0, 17.895169469, 78.757757664, 67.187563600 };
+  static float c[4] = { 1.4234372777, 2.7568153900, 1.3067284816, 0.17023821103 };
+  static float const1 = 0.180625;
+  static float const2 = 1.6;
+  static float d[3] = { 1.0, 0.73700164250, 0.12021132975 };
+  static float e[4] = { 6.6579051150, 3.0812263860, 0.42868294337, 0.017337203997 };
+  static float f[3] = { 1.0, 0.24197894225, 0.012258202635 };
+  float q;
+  float r;
+  static float split1 = 0.425;
+  static float split2 = 5.0;
+  float value;
+
+  if ( p <= 0.0 )
+  {
+    value = - r4_huge ( );
+    return value;
+  }
+
+  if ( 1.0 <= p )
+  {
+    value = r4_huge ( );
+    return value;
+  }
+
+  q = p - 0.5;
+
+  if ( fabs ( q ) <= split1 )
+  {
+    r = const1 - q * q;
+    value = q * r4poly_value ( 4, a, r ) / r4poly_value ( 4, b, r );
+  }
+  else
+  {
+    if ( q < 0.0 )
+    {
+      r = p;
+    }
+    else
+    {
+      r = 1.0 - p;
+    }
+
+    if ( r <= 0.0 )
+    {
+      value = - 1.0;
+      exit ( 1 );
+    }
+
+    r = sqrt ( -log ( r ) );
+
+    if ( r <= split2 )
+    {
+      r = r - const2;
+      value = r4poly_value ( 4, c, r ) / r4poly_value ( 3, d, r );
+    }
+    else
+    {
+      r = r - split2;
+      value = r4poly_value ( 4, e, r ) / r4poly_value ( 3, f, r );
+    }
+
+    if ( q < 0.0 )
+    {
+      value = - value;
+    }
+
+  }
+
+  return value;
+}
+/******************************************************************************/
+
+float r4poly_value ( int n, float a[], float x )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+    R4POLY_VALUE evaluates a real polynomial.
+
+  Discussion:
+
+    For sanity's sake, the value of N indicates the NUMBER of 
+    coefficients, or more precisely, the ORDER of the polynomial,
+    rather than the DEGREE of the polynomial.  The two quantities
+    differ by 1, but cause a great deal of confusion.
+
+    Given N and A, the form of the polynomial is:
+
+      p(x) = a[0] + a[1] * x + ... + a[n-2] * x^(n-2) + a[n-1] * x^(n-1)
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    13 August 2004
+
+  Author:
+
+    John Burkardt
+
+  Parameters:
+
+    Input, int N, the order of the polynomial.
+
+    Input, float A[N], the coefficients of the polynomial.
+    A[0] is the constant term.
+
+    Input, float X, the point at which the polynomial is to be evaluated.
+
+    Output, float R4POLY_VALUE, the value of the polynomial at X.
+*/
+{
+  int i;
+  float value;
+
+  value = 0.0;
+
+  for ( i = n - 1; 0 <= i; i-- )
+  {
+    value = value * x + a[i];
+  }
+
+  return value;
+}
+/******************************************************************************/
+
+double r8_huge ( )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+    R8_HUGE returns a "huge" R8.
+
+  Discussion:
+
+    The value returned by this function is NOT required to be the
+    maximum representable R8.  This value varies from machine to machine,
+    from compiler to compiler, and may cause problems when being printed.
+    We simply want a "very large" but non-infinite number.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    06 October 2007
+
+  Author:
+
+    John Burkardt
+
+  Parameters:
+
+    Output, double R8_HUGE, a "huge" R8 value.
+*/
+{
+  double value;
+
+  value = 1.0E+30;
+
+  return value;
+}
+/******************************************************************************/
+
+double r8_normal_01_cdf_inverse ( double p )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+    R8_NORMAL_01_CDF_INVERSE inverts the standard normal CDF.
+
+  Discussion:
+
+    The result is accurate to about 1 part in 10^16.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    19 March 2010
+
+  Author:
+
+    Original FORTRAN77 version by Michael Wichura.
+    C version by John Burkardt.
+
+  Reference:
+
+    Michael Wichura,
+    The Percentage Points of the Normal Distribution,
+    Algorithm AS 241,
+    Applied Statistics,
+    Volume 37, Number 3, pages 477-484, 1988.
+
+  Parameters:
+
+    Input, double P, the value of the cumulative probability 
+    densitity function.  0 < P < 1.  If P is outside this range, an "infinite"
+    value is returned.
+
+    Output, double R8_NORMAL_01_CDF_INVERSE, the normal deviate value 
+    with the property that the probability of a standard normal deviate being 
+    less than or equal to this value is P.
+*/
+{
+  static double a[8] = { 
+    3.3871328727963666080,     1.3314166789178437745e+2,
+    1.9715909503065514427e+3,  1.3731693765509461125e+4,
+    4.5921953931549871457e+4,  6.7265770927008700853e+4,
+    3.3430575583588128105e+4,  2.5090809287301226727e+3 };
+  static double b[8] = {
+    1.0,                       4.2313330701600911252e+1,
+    6.8718700749205790830e+2,  5.3941960214247511077e+3,
+    2.1213794301586595867e+4,  3.9307895800092710610e+4,
+    2.8729085735721942674e+4,  5.2264952788528545610e+3 };
+  static double c[8] = {
+    1.42343711074968357734,     4.63033784615654529590,
+    5.76949722146069140550,     3.64784832476320460504,
+    1.27045825245236838258,     2.41780725177450611770e-1,
+    2.27238449892691845833e-2,  7.74545014278341407640e-4 };
+  static double const1 = 0.180625;
+  static double const2 = 1.6;
+  static double d[8] = {
+    1.0,                        2.05319162663775882187,
+    1.67638483018380384940,     6.89767334985100004550e-1,
+    1.48103976427480074590e-1,  1.51986665636164571966e-2,
+    5.47593808499534494600e-4,  1.05075007164441684324e-9 };
+  static double e[8] = {
+    6.65790464350110377720,     5.46378491116411436990,
+    1.78482653991729133580,     2.96560571828504891230e-1,
+    2.65321895265761230930e-2,  1.24266094738807843860e-3,
+    2.71155556874348757815e-5,  2.01033439929228813265e-7 };
+  static double f[8] = {
+    1.0,                        5.99832206555887937690e-1,
+    1.36929880922735805310e-1,  1.48753612908506148525e-2,
+    7.86869131145613259100e-4,  1.84631831751005468180e-5,
+    1.42151175831644588870e-7,  2.04426310338993978564e-15 };
+  double q;
+  double r;
+  static double split1 = 0.425;
+  static double split2 = 5.0;
+  double value;
+
+  if ( p <= 0.0 )
+  {
+    value = - r8_huge ( );
+    return value;
+  }
+
+  if ( 1.0 <= p )
+  {
+    value = r8_huge ( );
+    return value;
+  }
+
+  q = p - 0.5;
+
+  if ( fabs ( q ) <= split1 )
+  {
+    r = const1 - q * q;
+    value = q * r8poly_value ( 8, a, r ) / r8poly_value ( 8, b, r );
+  }
+  else
+  {
+    if ( q < 0.0 )
+    {
+      r = p;
+    }
+    else
+    {
+      r = 1.0 - p;
+    }
+
+    if ( r <= 0.0 )
+    {
+      value = - 1.0;
+      exit ( 1 );
+    }
+
+    r = sqrt ( -log ( r ) );
+
+    if ( r <= split2 )
+    {
+      r = r - const2;
+      value = r8poly_value ( 8, c, r ) / r8poly_value ( 8, d, r ); 
+     }
+     else
+     {
+       r = r - split2;
+       value = r8poly_value ( 8, e, r ) / r8poly_value ( 8, f, r );
+    }
+
+    if ( q < 0.0 )
+    {
+      value = - value;
+    }
+
+  }
+
+  return value;
+}
+/******************************************************************************/
+
+double r8poly_value ( int n, double a[], double x )
+
+/******************************************************************************/
+/*
+  Purpose:
+
+    R8POLY_VALUE evaluates a double precision polynomial.
+
+  Discussion:
+
+    For sanity's sake, the value of N indicates the NUMBER of 
+    coefficients, or more precisely, the ORDER of the polynomial,
+    rather than the DEGREE of the polynomial.  The two quantities
+    differ by 1, but cause a great deal of confusion.
+
+    Given N and A, the form of the polynomial is:
+
+      p(x) = a[0] + a[1] * x + ... + a[n-2] * x^(n-2) + a[n-1] * x^(n-1)
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license. 
+
+  Modified:
+
+    19 March 2010
+
+  Author:
+
+    John Burkardt
+
+  Parameters:
+
+    Input, int N, the order of the polynomial.
+
+    Input, double A[N], the coefficients of the polynomial.
+    A[0] is the constant term.
+
+    Input, double X, the point at which the polynomial is to be evaluated.
+
+    Output, double R8POLY_VALUE, the value of the polynomial at X.
+*/
+{
+  int i;
+  double value;
+
+  value = 0.0;
+
+  for ( i = n - 1; 0 <= i; i-- )
+  {
+    value = value * x + a[i];
+  }
+
+  return value;
+}

src/asa241.h

@@ -0,0 +1,7 @@
+void normal_01_cdf_values ( int *n_data, double *x, double *fx );
+float r4_huge ( );
+float r4_normal_01_cdf_inverse ( float p );
+float r4poly_value ( int n, float a[], float x );
+double r8_huge ( );
+double r8_normal_01_cdf_inverse ( double p );
+double r8poly_value ( int n, double a[], double x );

src/gbm_mc.c

@@ -2,7 +2,9 @@
 #define max(X, Y) (((X) > (Y)) ? (X) : (Y))
 
 #include "gbm_mc.h"
+#include "sobol.h"
 #include "utils.h"
+#include "asa241.h"
 
 #include <math.h>
 #include <pthread.h>
@@ -103,23 +105,30 @@ void *run_simulations(void *opt_ptr)
 void gbm(struct Option *opt)
 {
 	int i = 0, j = 0;
+	int seed;
+	int nrandoms;
 	pthread_t *threads;
 	struct Option *options;
 	double **randoms;
 
+	seed = (unsigned)time(NULL);
 	options = malloc(sizeof(struct Option) * NUM_THREADS);
+	nrandoms = opt->sims / NUM_THREADS;
 	randoms = malloc(sizeof(double *) * NUM_THREADS);
 
 	/* Create 2D array */
 	for(i=0; i<NUM_THREADS; i++) {
 		randoms[i] = malloc(sizeof(double) * opt->sims / NUM_THREADS);
+		randoms[i] = i8_sobol_generate(1, nrandoms, (seed + nrandoms * i));
+		for(j=0; j<nrandoms; j++) {
+			randoms[i][j] = r8_normal_01_cdf_inverse(randoms[i][j]);
+		}
 	}
 
 	/* Fill 2D array with normal randoms
 	 * The purpose is to give the Option structs the random
 	 * numbers they will need for simulations, as opposed to
 	 * having the individual threads do so (rand() is not thread-safe)
-	 */
 	j = 0;
 	i = 0;
 	for(i = 0; i < NUM_THREADS;){
@@ -130,7 +139,7 @@ void gbm(struct Option *opt)
 		}
 		j++;
 	}
-
+	*/
 	/* Set the number of simulations on a per-thread basis
 	 */
 	opt->sims = opt->sims / NUM_THREADS;

src/sobol.h

@@ -0,0 +1,39 @@
+#include <stdint.h>
+
+#ifdef __cplusplus
+extern "C"
+{
+#endif
+
+int i4_bit_hi1 ( int n );
+int i4_bit_lo0 ( int n );
+int i4_max ( int i1, int i2 );
+int i4_min ( int i1, int i2 );
+void i4_sobol ( int dim_num, int *seed, float quasi[ ] );
+float *i4_sobol_generate ( int m, int n, int skip );
+int i4_uniform ( int b, int c, int *seed );
+
+int i8_bit_hi1 ( long long int n );
+int i8_bit_lo0 ( long long int n );
+long long int i8_max ( long long int i1, long long int i2 );
+long long int i8_min ( long long int i1, long long int i2 );
+void i8_sobol ( int dim_num, long long int *seed, double quasi[ ] );
+double *i8_sobol_generate ( int m, int n, int skip );
+long long int i8_uniform ( long long int b, long long int c, int *seed );
+
+float r4_abs ( float x );
+int r4_nint ( float x );
+float r4_uniform_01 ( int *seed );
+
+double r8_abs ( double x );
+int r8_nint ( double x );
+double r8_uniform_01 ( int *seed );
+
+/* void r8mat_write ( string output_filename, int m, int n, double table[] ); */
+
+int tau_sobol ( int dim_num );
+void timestamp ( );
+
+#ifdef __cplusplus
+}
+#endif