#include <math.h>
#include <assert.h>

#ifdef UNIT_TEST_NORMAL_DISTRIBUTION
#include <iostream>
#include <limits>
#include <vector>
#endif

#include "NormalDistribution.h"

NormalDistribution::NormalDistribution(int numberOfPoints)
{
  assert(numberOfPoints >= 7);
  _max = numberOfPoints - 1;
  _mean = _max / 2;
  _stdDev = sqrt(numberOfPoints)/2;
  _points = new double[numberOfPoints];
  // The following statement worked fine where N was 1000.  But you couldn't
  // go much higher before you exceeded the bounds of a double precision
  // number.  If N was 10000 then you'd have all 0's.
  //_points[0] = pow(0.5, _max);
  // Instead we start with 1 and divide by 2 each step along the way.
  _points[0] = 1.0;
  for (int i = 1; i < numberOfPoints; i++)
    {
      _points[i] = 0.0;
    }
  for (int row = 1; row < numberOfPoints; row++)
    {
      for (int i = row; i; i--)
	{
	  _points[i] = (_points[i] + _points[i - 1]) / 2.0;
	}
      _points[0] /= 2.0;
    }
  for (int i = 1; i < numberOfPoints; i++)
  {
    _points[i] += _points[i - 1];
  }
  #ifdef UNIT_TEST_NORMAL_DISTRIBUTION
  for (int i = 0; i < numberOfPoints; i++)
    {
      //std::cout<<i<<":  "<<_points[i]<<"\n";
    }
  #endif
}

NormalDistribution::~NormalDistribution()
{
  delete [] _points;
}

double NormalDistribution::probability(double from, double to) const
{
  return probability(to) - probability(from);
}

double NormalDistribution::probability(double to) const
{
  if (isnan(to))
    // We get this occassionally.  I'm not sure why.  Without this test we
    // would dump core.  Now if we get NaN as an input, we return NaN as an
    // output.
    return to;
  double _location = to * _stdDev + _mean;
  if (_location <= 0)
    {
      _location = 0;
    }
  else if (_location > _max)
    {
      _location = _max;
    }
  int lower = (int)floor(_location);
  double higherPart = _location - lower;
  if ((lower == _max) || (higherPart == 0))
    {
      return _points[lower];
    }
  return _points[lower] * (1 - higherPart) + _points[lower + 1] * higherPart;
}

#ifdef UNIT_TEST_NORMAL_DISTRIBUTION

NormalDistribution d(10000);

void showProbability(double from, double to)
{
  std::cout<<"probability("<<from<<","<<to<<")=";
  std::flush(std::cout);
  // The following line has been known to dump core.
  std::cout<<d.probability(from, to)<<std::endl;
};

typedef std::numeric_limits<double> limits;
typedef std::numeric_limits<long double> llimits;

int main(int argc, char *argv[])
{
  std::cout<<"probability(0)="<<d.probability(0)<<"\n";
  const std::vector< double > smallish
  {
    -1e+99, limits::lowest(), -limits::infinity(),
      limits::signaling_NaN(), -limits::signaling_NaN()
  };
  for (double small : smallish)
    showProbability(small, 0);
  for (double small : smallish)
    showProbability(small, -4);
  std::cout<<"probability(-4,-3)="<<d.probability(-4, -3)<<"\n";
  std::cout<<"probability(-3,-2)="<<d.probability(-3, -2)<<"\n";
  std::cout<<"probability(-2,-1)="<<d.probability(-2, -1)<<"\n";
  std::cout<<"probability(-1,0)="<<d.probability(-1, 0)<<"\n";
  std::cout<<"probability(0,1)="<<d.probability(0, 1)<<"\n";
  std::cout<<"probability(1,2)="<<d.probability(1, 2)<<"\n";
  std::cout<<"probability(2,3)="<<d.probability(2, 3)<<"\n";
  std::cout<<"probability(3,4)="<<d.probability(3, 4)<<"\n";
  const std::vector< double > bigish
  {
    1e+99, limits::max(), limits::infinity(),
      limits::signaling_NaN(), -limits::signaling_NaN()
  };
  for (double big : bigish)
    showProbability(4, big);
  for (double big : bigish)
    showProbability(0, big);
  std::cout<<"probability(0,0)="<<d.probability(0, 0)<<"\n";
  std::cout<<"probability(-1,1)="<<d.probability(-1, 1)<<"\n";
  std::cout<<"probability(-2,2)="<<d.probability(-2, 2)<<"\n";
  std::cout<<"probability(-3,3)="<<d.probability(-3, 3)<<"\n";
  std::cout<<"probability(-4,4)="<<d.probability(-4, 4)<<"\n";
  for (double big : bigish)
    for (double small : smallish)
      showProbability(small, big);
  std::cout<<"std::numeric_limits<double>\n"
	   <<"limits::max()="<<limits::max()<<"\n"
	   <<"limits::min()="<<limits::min()<<"\n"
	   <<"limits::lowest()="<<limits::lowest()<<"\n"
	   <<"limits::epsilon()="<<limits::epsilon()<<"\n"
	   <<"limits::digits="<<limits::digits<<"\n"
	   <<"limits::digits10="<<limits::digits10<<"\n"
	   <<"limits::max_exponent="<<limits::max_exponent<<"\n"
	   <<"limits::max_exponent10="<<limits::max_exponent10<<"\n"
	   <<"limits::min_exponent="<<limits::min_exponent<<"\n"
	   <<"limits::min_exponent10="<<limits::min_exponent10<<"\n"
	   <<"limits::has_infinity="<<limits::has_infinity<<"\n"
	   <<"limits::has_quiet_NaN="<<limits::has_quiet_NaN<<"\n"
	   <<"limits::has_signaling_NaN="<<limits::has_signaling_NaN<<"\n"
	   <<"limits::infinity()="<<limits::infinity()<<"\n"
	   <<"-limits::infinity()="<<-limits::infinity()<<"\n"
	   <<"limits::quiet_NaN()="<<limits::quiet_NaN()<<"\n"
	   <<"-limits::quiet_NaN()="<<-limits::quiet_NaN()<<"\n"
	   <<"limits::signaling_NaN()="<<limits::signaling_NaN()<<"\n"
	   <<"-limits::signaling_NaN()="<<-limits::signaling_NaN()<<"\n";
  std::cout<<"std::numeric_limits<long double>\n"
	   <<"llimits::max()="<<llimits::max()<<"\n"
	   <<"llimits::min()="<<llimits::min()<<"\n"
	   <<"llimits::lowest()="<<llimits::lowest()<<"\n"
	   <<"llimits::epsilon()="<<llimits::epsilon()<<"\n"
	   <<"llimits::digits="<<llimits::digits<<"\n"
	   <<"llimits::digits10="<<llimits::digits10<<"\n"
	   <<"llimits::max_exponent="<<llimits::max_exponent<<"\n"
	   <<"llimits::max_exponent10="<<llimits::max_exponent10<<"\n"
	   <<"llimits::min_exponent="<<llimits::min_exponent<<"\n"
	   <<"llimits::min_exponent10="<<llimits::min_exponent10<<"\n"
	   <<"llimits::has_infinity="<<llimits::has_infinity<<"\n"
	   <<"llimits::has_quiet_NaN="<<llimits::has_quiet_NaN<<"\n"
	   <<"llimits::has_signaling_NaN="<<llimits::has_signaling_NaN<<"\n"
	   <<"llimits::infinity()="<<llimits::infinity()<<"\n"
	   <<"-llimits::infinity()="<<-llimits::infinity()<<"\n"
	   <<"llimits::quiet_NaN()="<<llimits::quiet_NaN()<<"\n"
	   <<"-llimits::quiet_NaN()="<<-llimits::quiet_NaN()<<"\n"
	   <<"llimits::signaling_NaN()="<<llimits::signaling_NaN()<<"\n"
	   <<"-llimits::signaling_NaN()="<<-llimits::signaling_NaN()<<"\n";
  return 0;
}

#endif