#ifndef CHISQRDISTRI_CK
#define CHISQRDISTRI_OK

#include "adapIntgrtn.h"

using namespace adapInt;

namespace ChiSqrD {

	int df;

	double Chi_kernel( double x ) {
		double v = Math::Pow( x, df/2.0 - 1.0 ) * Math::Exp(-x/2.0);
		return v;
	}

	double gamma( double r ) {
		if ( r > 1.25 ) {
			return (r - 1.0) * gamma(r - 1.0);
		} else {
			if (r > 0.75) {
				return 1.0;
			} else {
				return Math::Sqrt(Math::PI);
			};
		};
	}

	double Normal_Dist( double x ) {
		return Math::Exp( -x * x / 2.0 );
	}

	double Cum_Chi_distribution( double chi ) {
		if ( chi <= 0.0 ) {
			return 0.0;
		} else {
			if ( df <= 1 ) {
				AdaptiveGL gl;
				return gl.integral( -Math::Sqrt(chi), Math::Sqrt(chi), Normal_Dist ) / Math::Sqrt(2.0 * Math::PI);
			} else {
				AdaptiveGL gl;
				double v = gl.integral( 0.0, chi, Chi_kernel );
				return v / (gamma(df/2.0) * Math::Pow(2.0, df/2.0));
			};
		};
	}


	//     Bisection法による根の計算    

	void Bisection( double (*f)(double),
					double s,					 //   fが増加関数のとき s > 0.0
					double c,					 //   根に対する関数値
					double L_b, double U_b,		 //   根の存在区間の下限と上限
					double & Root,				 //   根  f(Root) = c
					double acc ){
		if ( (f(L_b) - c) * (f(U_b) - c) > 0.0 ){
			throw  gcnew Exception("上限と下限の組み合わせが不適切です。");
		}
		do {
			double m = 0.5 * (L_b + U_b);
			double v = f(m);
			if ( s * (v - c) > 0.0 ) {
				U_b = m;
			} else {
				L_b = m;
			}
		} while ( (Math::Abs(U_b - L_b) >= acc * Math::Abs(U_b))
					&&
				  ( (Math::Abs(U_b) + Math::Abs(L_b)) >= 1.0e-14)
				);
		Root = 0.5 * (L_b + U_b);
	}




}

#endif
