using System;

//Качественный генератор нормальных распределений
namespace Entropy.Maths.RandomNumbers
{
    public class Ziggurat : INormalRandom
    {
        private const double LAST_POINT = 3.6541528853610088;
        private const double AREA = 0.00492867323399;
        private double[] _x = new double[257];
        private double[] _r = new double[256];
        private IUniformRandom _random;

        public Ziggurat(IUniformRandom random)
        {
            _random = random;
            this.Init();
        }

        public double Next
        {
            get
            {
                while (true)
                {
                    double u = 2.0 * _random.Next - 1.0;
                    int i = (int)(_random.Next * 256);
                    if (Math.Abs(u) < _r)
                        return u * _x;
                    if (i == 0)
                        return this.GetFromTail(u < 0);
                    double x = u * _x;
                    double f0 = Math.Exp(-0.5 * (_x * _x - x * x));
                    double f1 = Math.Exp(-0.5 * (_x[i + 1] * _x[i + 1] - x * x));
                    if (f1 + _random.Next * (f0 - f1) < 1.0)
                        return x;
                }
            }
        }

        private void Init()
        {
            double r = LAST_POINT;
            double f = Math.Exp(-0.5 * r * r);

            _x[0] = AREA / f;
            _x[1] = r;
            _x[256] = 0;

            for (int i = 2; i < 256; ++i)
            {
                _x = Math.Sqrt(-2.0 * Math.Log(AREA / _x[i - 1] + f));
                f = Math.Exp(-0.5 * _x * _x);
            }

            for (int i = 0; i < 256; ++i)
                _r = _x[i + 1] / _x;
        }

        private double GetFromTail(bool isNegative)
        {
            double x, y;
            do
            {
                x = Math.Log(_random.Next) / LAST_POINT;
                y = Math.Log(_random.Next);
            }
            while (-2.0 * y < x * x);
            return isNegative ? x - LAST_POINT : LAST_POINT - x;
        }
    }
}