using System;
namespace Volpe.Cafe.MonteCarlo.Model
{
    [Serializable]
    public class Random
    {
        #region 
        public Random(int seed)
        {
            this._seed = seed;
            this._ran2iff = 0;
            this._ir = new int[97];
            this._betaCon = new double[3];
            this._instance = 255;
        }
        #endregion
        #region 
        public double RanUnif()
        {
            int seed = this._seed;    
            int m = 714025;
            int ia = 1366;
            int ic = 150899;
            int iy = 0;
            double rm = 0.0000014005112D;
            double result;
            do
            {
                if (seed < 0 || this._ran2iff == 0)
                {
                    this._ran2iff = 1;
                    seed = Math.Abs((ic - seed) % m);
                    for (int i = 0; i < 97; i++)
                    {
                        seed = (ia * seed + ic) % m;
                        this._ir[i] = seed;
                    }
                    seed = (ia * seed + ic) % m;
                    iy = seed;
                }
                int j = 1 + (int)(96 * iy / m);
                if (j > 96 || j < 1) { throw new IndexOutOfRangeException("j=" + j + ", seed=" + this._seed); }
                iy = this._ir[j];
                result = iy * rm;
                seed = (ia * seed + ic) % m;
                this._ir[j] = seed;
            } while (Math.Abs(result - 0) < double.Epsilon);
            this._seed = seed;    
            return result;
        }
        public double RanBeta(double aShape, double bShape)
        {
            if (this._betaCon[0] <= 0)
            {
                this._betaCon[0] = (aShape < bShape) ? aShape : bShape;
                this._betaCon[0] = (this._betaCon[0] <= 1) ? 1 / this._betaCon[0] :
                    this._betaCon[0] = Math.Sqrt((aShape + bShape - 2) / (2 * aShape * bShape - aShape - bShape));
                this._betaCon[1] = aShape + bShape;
                this._betaCon[2] = aShape + 1 / this._betaCon[0];
            }
            double u1, u2, v, w, t;
            do
            {
                u1 = this.RanUnif();
                u2 = this.RanUnif();
                v  = this._betaCon[0] * Math.Log(u1 / (1 - u1));
                w  = aShape * Math.Exp(v);
                t  = Math.Log(this._betaCon[1] / (bShape + w)) + this._betaCon[2];
                v  = this._betaCon[0] * t;
            } while (this._betaCon[1] * t * v - 1.13862944 < Math.Log(u1 * u1 * u2));
            return w / (bShape + w);
        }
        public double RanBeta(double aShape, double bShape, double low, double high)
        {
            double value = this.RanBeta(aShape, bShape);
            return low + (high - low) * value;
        }
        public double RanNormal(double mean, double stdDev)
        {
            double rCode;
            if (this._instance < 254)
            {
                rCode = this._instance;
                this._instance = 255;
            }
            else
            {
                this._instance = Math.Sqrt(-2 * Math.Log(this.RanUnif()));
                double u2 = 6.2831852 * this.RanUnif();
                rCode = this._instance * Math.Cos(u2);      
                this._instance *= Math.Sin(u2);             
            }
            return mean + stdDev * rCode;
        }
        public double RanNormalConstrained(double mean, double stdDev, int retryIfOutOfBounds, double minBound, double maxBound,
            bool inclusiveBounds)
        {
            double rnd = this.RanNormal(mean, stdDev);
            if (retryIfOutOfBounds > 0)
            {
                if ((!double.IsNaN(minBound) && ((!inclusiveBounds && rnd <= minBound) || rnd < minBound)) ||
                    (!double.IsNaN(maxBound) && ((!inclusiveBounds && rnd >= maxBound) || rnd > maxBound)))
                {
                    return this.RanNormalConstrained(mean, stdDev, retryIfOutOfBounds - 1, minBound, maxBound, inclusiveBounds);
                }
            }
            return rnd;
        }
        #endregion
        #region 
        private int _seed;
        private int _ran2iff;        
        private int[] _ir;            
        private double[] _betaCon;    
        private double _instance;    
        #endregion
    }
}

