using System;

namespace alglib
{
    public class normaldistr
    {
        /*************************************************************************
        Error function

        The integral is

                                  x
                                   -
                        2         | |          2
          erf(x)  =  --------     |    exp( - t  ) dt.
                     sqrt(pi)   | |
                                 -
                                  0

        For 0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2); otherwise
        erf(x) = 1 - erfc(x).


        ACCURACY:

                             Relative error:
        arithmetic   domain     # trials      peak         rms
           IEEE      0,1         30000       3.7e-16     1.0e-16

        Cephes Math Library Release 2.8:  June, 2000
        Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
        *************************************************************************/
        public static double erf(double x)
        {
            double result = 0;
            double xsq = 0;
            double s = 0;
            double p = 0;
            double q = 0;

            s = Math.Sign(x);
            x = Math.Abs(x);
            if( (double)(x)<(double)(0.5) )
            {
                xsq = x*x;
                p = 0.007547728033418631287834;
                p = 0.288805137207594084924010+xsq*p;
                p = 14.3383842191748205576712+xsq*p;
                p = 38.0140318123903008244444+xsq*p;
                p = 3017.82788536507577809226+xsq*p;
                p = 7404.07142710151470082064+xsq*p;
                p = 80437.3630960840172832162+xsq*p;
                q = 0.0;
                q = 1.00000000000000000000000+xsq*q;
                q = 38.0190713951939403753468+xsq*q;
                q = 658.070155459240506326937+xsq*q;
                q = 6379.60017324428279487120+xsq*q;
                q = 34216.5257924628539769006+xsq*q;
                q = 80437.3630960840172826266+xsq*q;
                result = s*1.1283791670955125738961589031*x*p/q;
                return result;
            }
            if( (double)(x)>=(double)(10) )
            {
                result = s;
                return result;
            }
            result = s*(1-erfc(x));
            return result;
        }


        /*************************************************************************
        Complementary error function

         1 - erf(x) =

                                  inf.
                                    -
                         2         | |          2
          erfc(x)  =  --------     |    exp( - t  ) dt
                      sqrt(pi)   | |
                                  -
                                   x


        For small x, erfc(x) = 1 - erf(x); otherwise rational
        approximations are computed.


        ACCURACY:

                             Relative error:
        arithmetic   domain     # trials      peak         rms
           IEEE      0,26.6417   30000       5.7e-14     1.5e-14

        Cephes Math Library Release 2.8:  June, 2000
        Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
        *************************************************************************/
        public static double erfc(double x)
        {
            double result = 0;
            double p = 0;
            double q = 0;

            if( (double)(x)<(double)(0) )
            {
                result = 2-erfc(-x);
                return result;
            }
            if( (double)(x)<(double)(0.5) )
            {
                result = 1.0-erf(x);
                return result;
            }
            if( (double)(x)>=(double)(10) )
            {
                result = 0;
                return result;
            }
            p = 0.0;
            p = 0.5641877825507397413087057563+x*p;
            p = 9.675807882987265400604202961+x*p;
            p = 77.08161730368428609781633646+x*p;
            p = 368.5196154710010637133875746+x*p;
            p = 1143.262070703886173606073338+x*p;
            p = 2320.439590251635247384768711+x*p;
            p = 2898.0293292167655611275846+x*p;
            p = 1826.3348842295112592168999+x*p;
            q = 1.0;
            q = 17.14980943627607849376131193+x*q;
            q = 137.1255960500622202878443578+x*q;
            q = 661.7361207107653469211984771+x*q;
            q = 2094.384367789539593790281779+x*q;
            q = 4429.612803883682726711528526+x*q;
            q = 6089.5424232724435504633068+x*q;
            q = 4958.82756472114071495438422+x*q;
            q = 1826.3348842295112595576438+x*q;
            result = Math.Exp(-AP.Math.Sqr(x))*p/q;
            return result;
        }


        /*************************************************************************
        Normal distribution function

        Returns the area under the Gaussian probability density
        function, integrated from minus infinity to x:

                                   x
                                    -
                          1        | |          2
           ndtr(x)  = ---------    |    exp( - t /2 ) dt
                      sqrt(2pi)  | |
                                  -
                                 -inf.

                    =  ( 1 + erf(z) ) / 2
                    =  erfc(z) / 2

        where z = x/sqrt(2). Computation is via the functions
        erf and erfc.


        ACCURACY:

                             Relative error:
        arithmetic   domain     # trials      peak         rms
           IEEE     -13,0        30000       3.4e-14     6.7e-15

        Cephes Math Library Release 2.8:  June, 2000
        Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
        *************************************************************************/
        public static double normaldistribution(double x)
        {
            double result = 0;

            result = 0.5*(erf(x/1.41421356237309504880)+1);
            return result;
        }


        /*************************************************************************
        Inverse of the error function

        Cephes Math Library Release 2.8:  June, 2000
        Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
        *************************************************************************/
        public static double inverf(double e)
        {
            double result = 0;

            result = invnormaldistribution(0.5*(e+1))/Math.Sqrt(2);
            return result;
        }


        /*************************************************************************
        Inverse of Normal distribution function

        Returns the argument, x, for which the area under the
        Gaussian probability density function (integrated from
        minus infinity to x) is equal to y.


        For small arguments 0 < y < exp(-2), the program computes
        z = sqrt( -2.0 * log(y) );  then the approximation is
        x = z - log(z)/z  - (1/z) P(1/z) / Q(1/z).
        There are two rational functions P/Q, one for 0 < y < exp(-32)
        and the other for y up to exp(-2).  For larger arguments,
        w = y - 0.5, and  x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2)).

        ACCURACY:

                             Relative error:
        arithmetic   domain        # trials      peak         rms
           IEEE     0.125, 1        20000       7.2e-16     1.3e-16
           IEEE     3e-308, 0.135   50000       4.6e-16     9.8e-17

        Cephes Math Library Release 2.8:  June, 2000
        Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
        *************************************************************************/
        public static double invnormaldistribution(double y0)
        {
            double result = 0;
            double expm2 = 0;
            double s2pi = 0;
            double x = 0;
            double y = 0;
            double z = 0;
            double y2 = 0;
            double x0 = 0;
            double x1 = 0;
            int code = 0;
            double p0 = 0;
            double q0 = 0;
            double p1 = 0;
            double q1 = 0;
            double p2 = 0;
            double q2 = 0;

            expm2 = 0.13533528323661269189;
            s2pi = 2.50662827463100050242;
            if( (double)(y0)<=(double)(0) )
            {
                result = -AP.Math.MaxRealNumber;
                return result;
            }
            if( (double)(y0)>=(double)(1) )
            {
                result = AP.Math.MaxRealNumber;
                return result;
            }
            code = 1;
            y = y0;
            if( (double)(y)>(double)(1.0-expm2) )
            {
                y = 1.0-y;
                code = 0;
            }
            if( (double)(y)>(double)(expm2) )
            {
                y = y-0.5;
                y2 = y*y;
                p0 = -59.9633501014107895267;
                p0 = 98.0010754185999661536+y2*p0;
                p0 = -56.6762857469070293439+y2*p0;
                p0 = 13.9312609387279679503+y2*p0;
                p0 = -1.23916583867381258016+y2*p0;
                q0 = 1;
                q0 = 1.95448858338141759834+y2*q0;
                q0 = 4.67627912898881538453+y2*q0;
                q0 = 86.3602421390890590575+y2*q0;
                q0 = -225.462687854119370527+y2*q0;
                q0 = 200.260212380060660359+y2*q0;
                q0 = -82.0372256168333339912+y2*q0;
                q0 = 15.9056225126211695515+y2*q0;
                q0 = -1.18331621121330003142+y2*q0;
                x = y+y*y2*p0/q0;
                x = x*s2pi;
                result = x;
                return result;
            }
            x = Math.Sqrt(-(2.0*Math.Log(y)));
            x0 = x-Math.Log(x)/x;
            z = 1.0/x;
            if( (double)(x)<(double)(8.0) )
            {
                p1 = 4.05544892305962419923;
                p1 = 31.5251094599893866154+z*p1;
                p1 = 57.1628192246421288162+z*p1;
                p1 = 44.0805073893200834700+z*p1;
                p1 = 14.6849561928858024014+z*p1;
                p1 = 2.18663306850790267539+z*p1;
                p1 = -(1.40256079171354495875*0.1)+z*p1;
                p1 = -(3.50424626827848203418*0.01)+z*p1;
                p1 = -(8.57456785154685413611*0.0001)+z*p1;
                q1 = 1;
                q1 = 15.7799883256466749731+z*q1;
                q1 = 45.3907635128879210584+z*q1;
                q1 = 41.3172038254672030440+z*q1;
                q1 = 15.0425385692907503408+z*q1;
                q1 = 2.50464946208309415979+z*q1;
                q1 = -(1.42182922854787788574*0.1)+z*q1;
                q1 = -(3.80806407691578277194*0.01)+z*q1;
                q1 = -(9.33259480895457427372*0.0001)+z*q1;
                x1 = z*p1/q1;
            }
            else
            {
                p2 = 3.23774891776946035970;
                p2 = 6.91522889068984211695+z*p2;
                p2 = 3.93881025292474443415+z*p2;
                p2 = 1.33303460815807542389+z*p2;
                p2 = 2.01485389549179081538*0.1+z*p2;
                p2 = 1.23716634817820021358*0.01+z*p2;
                p2 = 3.01581553508235416007*0.0001+z*p2;
                p2 = 2.65806974686737550832*0.000001+z*p2;
                p2 = 6.23974539184983293730*0.000000001+z*p2;
                q2 = 1;
                q2 = 6.02427039364742014255+z*q2;
                q2 = 3.67983563856160859403+z*q2;
                q2 = 1.37702099489081330271+z*q2;
                q2 = 2.16236993594496635890*0.1+z*q2;
                q2 = 1.34204006088543189037*0.01+z*q2;
                q2 = 3.28014464682127739104*0.0001+z*q2;
                q2 = 2.89247864745380683936*0.000001+z*q2;
                q2 = 6.79019408009981274425*0.000000001+z*q2;
                x1 = z*p2/q2;
            }
            x = x0-x1;
            if( code!=0 )
            {
                x = -x;
            }
            result = x;
            return result;
        }
    }
}

/*************************************************************************
AP library
Copyright (c) 2003-2009 Sergey Bochkanov (ALGLIB project).

>>> LICENSE >>>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation (www.fsf.org); either version 2 of the
License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

A copy of the GNU General Public License is available at
http://www.fsf.org/licensing/licenses

>>> END OF LICENSE >>>
*************************************************************************/
namespace AP
{
    /********************************************************************
    Class defining a complex number with double precision.
    ********************************************************************/
    public struct Complex
    {
        public double x;
        public double y;

        public Complex(double _x)
        {
            x = _x;
            y = 0;
        }
        public Complex(double _x, double _y)
        {
            x = _x;
            y = _y;
        }
        public static implicit operator Complex(double _x)
        {
            return new Complex(_x);
        }
        public static bool operator==(Complex lhs, Complex rhs)
        {
            return ((double)lhs.x==(double)rhs.x) & ((double)lhs.y==(double)rhs.y);
        }
        public static bool operator!=(Complex lhs, Complex rhs)
        {
            return ((double)lhs.x!=(double)rhs.x) | ((double)lhs.y!=(double)rhs.y);
        }
        public static Complex operator+(Complex lhs)
        {
            return lhs;
        }
        public static Complex operator-(Complex lhs)
        {
            return new Complex(-lhs.x,-lhs.y);
        }
        public static Complex operator+(Complex lhs, Complex rhs)
        {
            return new Complex(lhs.x+rhs.x,lhs.y+rhs.y);
        }
        public static Complex operator-(Complex lhs, Complex rhs)
        {
            return new Complex(lhs.x-rhs.x,lhs.y-rhs.y);
        }
        public static Complex operator*(Complex lhs, Complex rhs)
        { 
            return new Complex(lhs.x*rhs.x-lhs.y*rhs.y, lhs.x*rhs.y+lhs.y*rhs.x);
        }
        public static Complex operator/(Complex lhs, Complex rhs)
        {
            Complex result;
            double e;
            double f;
            if( System.Math.Abs(rhs.y)<System.Math.Abs(rhs.x) )
            {
                e = rhs.y/rhs.x;
                f = rhs.x+rhs.y*e;
                result.x = (lhs.x+lhs.y*e)/f;
                result.y = (lhs.y-lhs.x*e)/f;
            }
            else
            {
                e = rhs.x/rhs.y;
                f = rhs.y+rhs.x*e;
                result.x = (lhs.y+lhs.x*e)/f;
                result.y = (-lhs.x+lhs.y*e)/f;
            }
            return result;
        }
		public override int GetHashCode() 
		{ 
			return x.GetHashCode() ^ y.GetHashCode(); 
		}
		public override bool Equals(object obj) 
		{ 
			if( obj is byte)
				return Equals(new Complex((byte)obj));
			if( obj is sbyte)
				return Equals(new Complex((sbyte)obj));
			if( obj is short)
				return Equals(new Complex((short)obj));
			if( obj is ushort)
				return Equals(new Complex((ushort)obj));
			if( obj is int)
				return Equals(new Complex((int)obj));
			if( obj is uint)
				return Equals(new Complex((uint)obj));
			if( obj is long)
				return Equals(new Complex((long)obj));
			if( obj is ulong)
				return Equals(new Complex((ulong)obj));
			if( obj is float)
				return Equals(new Complex((float)obj));
			if( obj is double)
				return Equals(new Complex((double)obj));
			if( obj is decimal)
				return Equals(new Complex((double)(decimal)obj));
			return base.Equals(obj); 
		}	
	}    
    
	/********************************************************************
	AP math namespace
	********************************************************************/
	public struct rcommstate
	{
		public int stage;
		public int[] ia;
		public bool[] ba;
		public double[] ra;
		public AP.Complex[] ca;
	};

	/********************************************************************
    AP math namespace
    ********************************************************************/
    public class Math
    {
        //public static System.Random RndObject = new System.Random(System.DateTime.Now.Millisecond);
        public static System.Random RndObject = new System.Random(System.DateTime.Now.Millisecond + 1000*System.DateTime.Now.Second + 60*1000*System.DateTime.Now.Minute);

        public const double MachineEpsilon = 5E-16;
        public const double MaxRealNumber = 1E300;
        public const double MinRealNumber = 1E-300;
        
        public static double RandomReal()
        {
            double r = 0;
            lock(RndObject){ r = RndObject.NextDouble(); }
            return r;
        }
        public static int RandomInteger(int N)
        {
            int r = 0;
            lock(RndObject){ r = RndObject.Next(N); }
            return r;
        }
        public static double Sqr(double X)
        {
            return X*X;
        }        
        public static double AbsComplex(Complex z)
        {
            double w;
            double xabs;
            double yabs;
            double v;
    
            xabs = System.Math.Abs(z.x);
            yabs = System.Math.Abs(z.y);
            w = xabs>yabs ? xabs : yabs;
            v = xabs<yabs ? xabs : yabs; 
            if( v==0 )
                return w;
            else
            {
                double t = v/w;
                return w*System.Math.Sqrt(1+t*t);
            }
        }
        public static Complex Conj(Complex z)
        {
            return new Complex(z.x, -z.y); 
        }    
        public static Complex CSqr(Complex z)
        {
            return new Complex(z.x*z.x-z.y*z.y, 2*z.x*z.y); 
        }

    }
}
 

 
using System;

namespace alglib
{
    public class normaldistr
    {
        /*************************************************************************
        Error function

        The integral is

                                  x
                                   -
                        2         | |          2
          erf(x)  =  --------     |    exp( - t  ) dt.
                     sqrt(pi)   | |
                                 -
                                  0

        For 0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2); otherwise
        erf(x) = 1 - erfc(x).


        ACCURACY:

                             Relative error:
        arithmetic   domain     # trials      peak         rms
           IEEE      0,1         30000       3.7e-16     1.0e-16

        Cephes Math Library Release 2.8:  June, 2000
        Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
        *************************************************************************/
        public static double erf(double x)
        {
            double result = 0;
            double xsq = 0;
            double s = 0;
            double p = 0;
            double q = 0;

            s = Math.Sign(x);
            x = Math.Abs(x);
            if( (double)(x)<(double)(0.5) )
            {
                xsq = x*x;
                p = 0.007547728033418631287834;
                p = 0.288805137207594084924010+xsq*p;
                p = 14.3383842191748205576712+xsq*p;
                p = 38.0140318123903008244444+xsq*p;
                p = 3017.82788536507577809226+xsq*p;
                p = 7404.07142710151470082064+xsq*p;
                p = 80437.3630960840172832162+xsq*p;
                q = 0.0;
                q = 1.00000000000000000000000+xsq*q;
                q = 38.0190713951939403753468+xsq*q;
                q = 658.070155459240506326937+xsq*q;
                q = 6379.60017324428279487120+xsq*q;
                q = 34216.5257924628539769006+xsq*q;
                q = 80437.3630960840172826266+xsq*q;
                result = s*1.1283791670955125738961589031*x*p/q;
                return result;
            }
            if( (double)(x)>=(double)(10) )
            {
                result = s;
                return result;
            }
            result = s*(1-erfc(x));
            return result;
        }


        /*************************************************************************
        Complementary error function

         1 - erf(x) =

                                  inf.
                                    -
                         2         | |          2
          erfc(x)  =  --------     |    exp( - t  ) dt
                      sqrt(pi)   | |
                                  -
                                   x


        For small x, erfc(x) = 1 - erf(x); otherwise rational
        approximations are computed.


        ACCURACY:

                             Relative error:
        arithmetic   domain     # trials      peak         rms
           IEEE      0,26.6417   30000       5.7e-14     1.5e-14

        Cephes Math Library Release 2.8:  June, 2000
        Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
        *************************************************************************/
        public static double erfc(double x)
        {
            double result = 0;
            double p = 0;
            double q = 0;

            if( (double)(x)<(double)(0) )
            {
                result = 2-erfc(-x);
                return result;
            }
            if( (double)(x)<(double)(0.5) )
            {
                result = 1.0-erf(x);
                return result;
            }
            if( (double)(x)>=(double)(10) )
            {
                result = 0;
                return result;
            }
            p = 0.0;
            p = 0.5641877825507397413087057563+x*p;
            p = 9.675807882987265400604202961+x*p;
            p = 77.08161730368428609781633646+x*p;
            p = 368.5196154710010637133875746+x*p;
            p = 1143.262070703886173606073338+x*p;
            p = 2320.439590251635247384768711+x*p;
            p = 2898.0293292167655611275846+x*p;
            p = 1826.3348842295112592168999+x*p;
            q = 1.0;
            q = 17.14980943627607849376131193+x*q;
            q = 137.1255960500622202878443578+x*q;
            q = 661.7361207107653469211984771+x*q;
            q = 2094.384367789539593790281779+x*q;
            q = 4429.612803883682726711528526+x*q;
            q = 6089.5424232724435504633068+x*q;
            q = 4958.82756472114071495438422+x*q;
            q = 1826.3348842295112595576438+x*q;
            result = Math.Exp(-AP.Math.Sqr(x))*p/q;
            return result;
        }


        /*************************************************************************
        Normal distribution function

        Returns the area under the Gaussian probability density
        function, integrated from minus infinity to x:

                                   x
                                    -
                          1        | |          2
           ndtr(x)  = ---------    |    exp( - t /2 ) dt
                      sqrt(2pi)  | |
                                  -
                                 -inf.

                    =  ( 1 + erf(z) ) / 2
                    =  erfc(z) / 2

        where z = x/sqrt(2). Computation is via the functions
        erf and erfc.


        ACCURACY:

                             Relative error:
        arithmetic   domain     # trials      peak         rms
           IEEE     -13,0        30000       3.4e-14     6.7e-15

        Cephes Math Library Release 2.8:  June, 2000
        Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
        *************************************************************************/
        public static double normaldistribution(double x)
        {
            double result = 0;

            result = 0.5*(erf(x/1.41421356237309504880)+1);
            return result;
        }


        /*************************************************************************
        Inverse of the error function

        Cephes Math Library Release 2.8:  June, 2000
        Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
        *************************************************************************/
        public static double inverf(double e)
        {
            double result = 0;

            result = invnormaldistribution(0.5*(e+1))/Math.Sqrt(2);
            return result;
        }


        /*************************************************************************
        Inverse of Normal distribution function

        Returns the argument, x, for which the area under the
        Gaussian probability density function (integrated from
        minus infinity to x) is equal to y.


        For small arguments 0 < y < exp(-2), the program computes
        z = sqrt( -2.0 * log(y) );  then the approximation is
        x = z - log(z)/z  - (1/z) P(1/z) / Q(1/z).
        There are two rational functions P/Q, one for 0 < y < exp(-32)
        and the other for y up to exp(-2).  For larger arguments,
        w = y - 0.5, and  x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2)).

        ACCURACY:

                             Relative error:
        arithmetic   domain        # trials      peak         rms
           IEEE     0.125, 1        20000       7.2e-16     1.3e-16
           IEEE     3e-308, 0.135   50000       4.6e-16     9.8e-17

        Cephes Math Library Release 2.8:  June, 2000
        Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
        *************************************************************************/
        public static double invnormaldistribution(double y0)
        {
            double result = 0;
            double expm2 = 0;
            double s2pi = 0;
            double x = 0;
            double y = 0;
            double z = 0;
            double y2 = 0;
            double x0 = 0;
            double x1 = 0;
            int code = 0;
            double p0 = 0;
            double q0 = 0;
            double p1 = 0;
            double q1 = 0;
            double p2 = 0;
            double q2 = 0;

            expm2 = 0.13533528323661269189;
            s2pi = 2.50662827463100050242;
            if( (double)(y0)<=(double)(0) )
            {
                result = -AP.Math.MaxRealNumber;
                return result;
            }
            if( (double)(y0)>=(double)(1) )
            {
                result = AP.Math.MaxRealNumber;
                return result;
            }
            code = 1;
            y = y0;
            if( (double)(y)>(double)(1.0-expm2) )
            {
                y = 1.0-y;
                code = 0;
            }
            if( (double)(y)>(double)(expm2) )
            {
                y = y-0.5;
                y2 = y*y;
                p0 = -59.9633501014107895267;
                p0 = 98.0010754185999661536+y2*p0;
                p0 = -56.6762857469070293439+y2*p0;
                p0 = 13.9312609387279679503+y2*p0;
                p0 = -1.23916583867381258016+y2*p0;
                q0 = 1;
                q0 = 1.95448858338141759834+y2*q0;
                q0 = 4.67627912898881538453+y2*q0;
                q0 = 86.3602421390890590575+y2*q0;
                q0 = -225.462687854119370527+y2*q0;
                q0 = 200.260212380060660359+y2*q0;
                q0 = -82.0372256168333339912+y2*q0;
                q0 = 15.9056225126211695515+y2*q0;
                q0 = -1.18331621121330003142+y2*q0;
                x = y+y*y2*p0/q0;
                x = x*s2pi;
                result = x;
                return result;
            }
            x = Math.Sqrt(-(2.0*Math.Log(y)));
            x0 = x-Math.Log(x)/x;
            z = 1.0/x;
            if( (double)(x)<(double)(8.0) )
            {
                p1 = 4.05544892305962419923;
                p1 = 31.5251094599893866154+z*p1;
                p1 = 57.1628192246421288162+z*p1;
                p1 = 44.0805073893200834700+z*p1;
                p1 = 14.6849561928858024014+z*p1;
                p1 = 2.18663306850790267539+z*p1;
                p1 = -(1.40256079171354495875*0.1)+z*p1;
                p1 = -(3.50424626827848203418*0.01)+z*p1;
                p1 = -(8.57456785154685413611*0.0001)+z*p1;
                q1 = 1;
                q1 = 15.7799883256466749731+z*q1;
                q1 = 45.3907635128879210584+z*q1;
                q1 = 41.3172038254672030440+z*q1;
                q1 = 15.0425385692907503408+z*q1;
                q1 = 2.50464946208309415979+z*q1;
                q1 = -(1.42182922854787788574*0.1)+z*q1;
                q1 = -(3.80806407691578277194*0.01)+z*q1;
                q1 = -(9.33259480895457427372*0.0001)+z*q1;
                x1 = z*p1/q1;
            }
            else
            {
                p2 = 3.23774891776946035970;
                p2 = 6.91522889068984211695+z*p2;
                p2 = 3.93881025292474443415+z*p2;
                p2 = 1.33303460815807542389+z*p2;
                p2 = 2.01485389549179081538*0.1+z*p2;
                p2 = 1.23716634817820021358*0.01+z*p2;
                p2 = 3.01581553508235416007*0.0001+z*p2;
                p2 = 2.65806974686737550832*0.000001+z*p2;
                p2 = 6.23974539184983293730*0.000000001+z*p2;
                q2 = 1;
                q2 = 6.02427039364742014255+z*q2;
                q2 = 3.67983563856160859403+z*q2;
                q2 = 1.37702099489081330271+z*q2;
                q2 = 2.16236993594496635890*0.1+z*q2;
                q2 = 1.34204006088543189037*0.01+z*q2;
                q2 = 3.28014464682127739104*0.0001+z*q2;
                q2 = 2.89247864745380683936*0.000001+z*q2;
                q2 = 6.79019408009981274425*0.000000001+z*q2;
                x1 = z*p2/q2;
            }
            x = x0-x1;
            if( code!=0 )
            {
                x = -x;
            }
            result = x;
            return result;
        }
    }
}

/*************************************************************************
AP library
Copyright (c) 2003-2009 Sergey Bochkanov (ALGLIB project).

>>> LICENSE >>>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation (www.fsf.org); either version 2 of the
License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

A copy of the GNU General Public License is available at
http://www.fsf.org/licensing/licenses

>>> END OF LICENSE >>>
*************************************************************************/
namespace AP
{
    /********************************************************************
    Class defining a complex number with double precision.
    ********************************************************************/
    public struct Complex
    {
        public double x;
        public double y;

        public Complex(double _x)
        {
            x = _x;
            y = 0;
        }
        public Complex(double _x, double _y)
        {
            x = _x;
            y = _y;
        }
        public static implicit operator Complex(double _x)
        {
            return new Complex(_x);
        }
        public static bool operator==(Complex lhs, Complex rhs)
        {
            return ((double)lhs.x==(double)rhs.x) & ((double)lhs.y==(double)rhs.y);
        }
        public static bool operator!=(Complex lhs, Complex rhs)
        {
            return ((double)lhs.x!=(double)rhs.x) | ((double)lhs.y!=(double)rhs.y);
        }
        public static Complex operator+(Complex lhs)
        {
            return lhs;
        }
        public static Complex operator-(Complex lhs)
        {
            return new Complex(-lhs.x,-lhs.y);
        }
        public static Complex operator+(Complex lhs, Complex rhs)
        {
            return new Complex(lhs.x+rhs.x,lhs.y+rhs.y);
        }
        public static Complex operator-(Complex lhs, Complex rhs)
        {
            return new Complex(lhs.x-rhs.x,lhs.y-rhs.y);
        }
        public static Complex operator*(Complex lhs, Complex rhs)
        { 
            return new Complex(lhs.x*rhs.x-lhs.y*rhs.y, lhs.x*rhs.y+lhs.y*rhs.x);
        }
        public static Complex operator/(Complex lhs, Complex rhs)
        {
            Complex result;
            double e;
            double f;
            if( System.Math.Abs(rhs.y)<System.Math.Abs(rhs.x) )
            {
                e = rhs.y/rhs.x;
                f = rhs.x+rhs.y*e;
                result.x = (lhs.x+lhs.y*e)/f;
                result.y = (lhs.y-lhs.x*e)/f;
            }
            else
            {
                e = rhs.x/rhs.y;
                f = rhs.y+rhs.x*e;
                result.x = (lhs.y+lhs.x*e)/f;
                result.y = (-lhs.x+lhs.y*e)/f;
            }
            return result;
        }
		public override int GetHashCode() 
		{ 
			return x.GetHashCode() ^ y.GetHashCode(); 
		}
		public override bool Equals(object obj) 
		{ 
			if( obj is byte)
				return Equals(new Complex((byte)obj));
			if( obj is sbyte)
				return Equals(new Complex((sbyte)obj));
			if( obj is short)
				return Equals(new Complex((short)obj));
			if( obj is ushort)
				return Equals(new Complex((ushort)obj));
			if( obj is int)
				return Equals(new Complex((int)obj));
			if( obj is uint)
				return Equals(new Complex((uint)obj));
			if( obj is long)
				return Equals(new Complex((long)obj));
			if( obj is ulong)
				return Equals(new Complex((ulong)obj));
			if( obj is float)
				return Equals(new Complex((float)obj));
			if( obj is double)
				return Equals(new Complex((double)obj));
			if( obj is decimal)
				return Equals(new Complex((double)(decimal)obj));
			return base.Equals(obj); 
		}	
	}    
    
	/********************************************************************
	AP math namespace
	********************************************************************/
	public struct rcommstate
	{
		public int stage;
		public int[] ia;
		public bool[] ba;
		public double[] ra;
		public AP.Complex[] ca;
	};

	/********************************************************************
    AP math namespace
    ********************************************************************/
    public class Math
    {
        //public static System.Random RndObject = new System.Random(System.DateTime.Now.Millisecond);
        public static System.Random RndObject = new System.Random(System.DateTime.Now.Millisecond + 1000*System.DateTime.Now.Second + 60*1000*System.DateTime.Now.Minute);

        public const double MachineEpsilon = 5E-16;
        public const double MaxRealNumber = 1E300;
        public const double MinRealNumber = 1E-300;
        
        public static double RandomReal()
        {
            double r = 0;
            lock(RndObject){ r = RndObject.NextDouble(); }
            return r;
        }
        public static int RandomInteger(int N)
        {
            int r = 0;
            lock(RndObject){ r = RndObject.Next(N); }
            return r;
        }
        public static double Sqr(double X)
        {
            return X*X;
        }        
        public static double AbsComplex(Complex z)
        {
            double w;
            double xabs;
            double yabs;
            double v;
    
            xabs = System.Math.Abs(z.x);
            yabs = System.Math.Abs(z.y);
            w = xabs>yabs ? xabs : yabs;
            v = xabs<yabs ? xabs : yabs; 
            if( v==0 )
                return w;
            else
            {
                double t = v/w;
                return w*System.Math.Sqrt(1+t*t);
            }
        }
        public static Complex Conj(Complex z)
        {
            return new Complex(z.x, -z.y); 
        }    
        public static Complex CSqr(Complex z)
        {
            return new Complex(z.x*z.x-z.y*z.y, 2*z.x*z.y); 
        }

    }
}