/********************************/
/*     χ2分布の計算            */
/*          coded by Y.Suganuma */
/********************************/
import java.io.*;

public class Chi {

	static double p;   // α％値を計算するとき時α/100を設定
	static int dof;    // 自由度

	public static void main(String args[]) throws IOException
	{
		BufferedReader in = new BufferedReader(new InputStreamReader(System.in));
		double h, x, x1, x2, y;
		double z[] = new double [1];
		int sw, sw1[] = new int [1];

		System.out.print("自由度は？ ");
		dof = Integer.parseInt(in.readLine());
		System.out.println("目的とする結果は？ ");
		System.out.print("     =0 : 確率の計算（ P(X = x) 及び P(X < x) の値）\n");
		System.out.print("     =1 : ｐ％値（ P(X > u) = 0.01p となるuの値） ");
		sw = Integer.parseInt(in.readLine());

		if (sw == 0) {

			System.out.print("グラフ出力？（=1: yes,  =0: no） ");
			sw = Integer.parseInt(in.readLine());
					// 密度関数と分布関数の値
			if (sw == 0) {
				System.out.print("   データは？ ");
				x = Double.parseDouble(in.readLine());
				y = App.chi(x, z, dof);
				System.out.println("P(X = " + x + ") = " + z[0] + ",  P( X < " + x + ") = " + y +
                                   " (自由度 = " + dof + ")");
			}
					// グラフ出力
			else {
				String file1, file2;
				System.out.print("   密度関数のファイル名は？ ");
				file1 = in.readLine();
				System.out.print("   分布関数のファイル名は？ ");
				file2 = in.readLine();
				PrintStream out1 = new PrintStream(new FileOutputStream(file1));
				PrintStream out2 = new PrintStream(new FileOutputStream(file2));
				System.out.print("   データの下限は？ ");
				x1 = Double.parseDouble(in.readLine());
				System.out.print("   データの上限は？ ");
				x2 = Double.parseDouble(in.readLine());
				System.out.print("   刻み幅は？ ");
				h = Double.parseDouble(in.readLine());
				for (x = x1; x < x2+0.5*h; x += h) {
					y = App.chi(x, z, dof);
					out1.println(x + " " + z[0]);
					out2.println(x + " " + y);
				}
			}
		}
					// ％値
		else {
			System.out.print("％の値は？ ");
			x = Double.parseDouble(in.readLine());
			p = 0.01 * x;
			y = App.p_chi(sw1);
			System.out.println(x + "％値 = " + y + "  sw " + sw1[0] + " (自由度 = " + dof + ")");
		}
	}
}

-------------------------------------------------

/********************/
/*     関数値の計算 */
/********************/
class Kansu {
	private int sw;
					// コンストラクタ
	Kansu (int s) {sw = s;}
					// double型関数
	double snx(double x)
	{

		double y = 0.0, w[] = new double [1];

		switch (sw) {
						// 関数値（f(x)）の計算（正規分布）
			case 0:
				y = 1.0 - Chi.p - App.normal(x, w);
				break;
						// 関数値（f(x)）の計算（χ2分布）
			case 1:
				y = App.chi(x, w, Chi.dof) - 1.0 + Chi.p;
				break;
						// 関数の微分の計算（χ2分布）
			case 2:
				y = App.chi(x, w, Chi.dof);
				y = w[0];
				break;
		}

		return y;
	}
}

--------------------------------------------------

/****************************/
/*     科学技術系算用の手法 */
/****************************/
class App {

	/************************************************/
	/*     χ２乗分布の計算（P(X = ch), P(X < ch)） */
	/*          dd : P(X = ch)                      */
	/*          df : 自由度                         */
	/*          return : P(X < ch)                  */
	/************************************************/
	static double chi(double ch, double dd[], int df)
	{
		double pi = Math.PI;
		double chs, pis, pp, x, u, y[] = new double [1];
		int ia, i1;

		if (ch < 1.0e-10)
			ch = 1.0e-10;

		pis = Math.sqrt(pi);
		chs = Math.sqrt(ch);
		x   = 0.5 * ch;

		if(df%2 != 0) {
			u  = Math.sqrt(x) * Math.exp(-x) / pis;
			pp = 2.0 * (normal(chs, y) - 0.5);
			ia = 1;
		}

		else {
			u  = x * Math.exp(-x);
			pp = 1.0 - Math.exp(-x);
			ia = 2;
		}

		if (ia != df) {
			for (i1 = ia; i1 <= df-2; i1 += 2) {
				pp -= 2.0 * u / i1;
				u  *= (ch / i1);
			}
		}

		dd[0] = u / ch;

		return pp;
	}

	/**********************************************/
	/*     χ２乗分布のｐ％値（P(X > u) = 0.01p） */
	/*          ind : >=0 : normal（収束回数）    */
	/*                =-1 : 収束しなかった        */
	/**********************************************/
	static double p_chi(int ind[])
	{
		double po, x, xx = 0.0, x0, w;
					// 自由度＝１（正規分布を使用）
		if (Chi.dof == 1) {
			po     = Chi.p;
			Chi.p *= 0.5;
			x      = p_normal(ind);
			xx     = x * x;
			Chi.p  = po;
		}

		else {
					// 自由度＝２
			if (Chi.dof == 2)
				xx = -2.0 * Math.log(Chi.p);
					// 自由度＞２
			else {

				x  = p_normal(ind);   // 初期値計算のため

				if (ind[0] >= 0) {
						// ニュートン法
					w  = 2.0 / (9.0 * Chi.dof);
					x  = 1.0 - w + x * Math.sqrt(w);
					x0 = Math.pow(x, 3.0) * Chi.dof;

					Kansu kn1 = new Kansu(1);
					Kansu kn2 = new Kansu(2);

					xx = newton(x0, 1.0e-6, 1.0e-10, 100, ind, kn1, kn2);
				}
			}
		}

		return xx;
	}

	/*****************************************************/
	/*     標準正規分布N(0,1)の計算（P(X = x), P(X < x)）*/
	/*          w : P(X = x)                             */
	/*          return : P(X < x)                        */
	/*****************************************************/
	static double normal(double x, double w[])
	{
		double y, z, P;
	/*
	     確率密度関数（定義式）
	*/
		w[0] = Math.exp(-0.5 * x * x) / Math.sqrt(2.0*Math.PI);
	/*
	     確率分布関数（近似式を使用）
	*/
		y = 0.70710678118654 * Math.abs(x);
		z = 1.0 + y * (0.0705230784 + y * (0.0422820123 +
            y * (0.0092705272 + y * (0.0001520143 + y * (0.0002765672 +
            y * 0.0000430638)))));
		P = 1.0 - Math.pow(z, -16.0);

		if (x < 0.0)
			P = 0.5 - 0.5 * P;
		else
			P = 0.5 + 0.5 * P;

		return P;
	}

	/**********************************************************************/
	/*     標準正規分布N(0,1)のｐ％値（P(X > u) = 0.01p）（二分法を使用） */
	/*          ind : >= 0 : normal（収束回数）                           */
	/*                = -1 : 収束しなかった                               */
	/**********************************************************************/
	static double p_normal(int ind[])
	{
		double u;
		int sw[] = new int [1];

		Kansu kn = new Kansu(0);

		u        = bisection(-7.0, 7.0, 1.0e-6, 1.0e-10, 100, sw, kn);
		ind[0]   = sw[0];

		return u;
	}

	/*********************************************************/
	/*     Newton法による非線形方程式(f(x)=0)の解            */
	/*          x1 : 初期値                                  */
	/*          eps1 : 終了条件１（｜x(k+1)-x(k)｜＜eps1）   */
	/*          eps2 : 終了条件２（｜f(x(k))｜＜eps2）       */
	/*          max : 最大試行回数                           */
	/*          ind : 実際の試行回数                         */
	/*                （負の時は解を得ることができなかった） */
	/*          kn1 : 関数を計算するクラスオブジェクト       */
	/*          kn2 : 関数の微分を計算するクラスオブジェクト */
	/*          return : 解                                  */
	/*********************************************************/
	static double newton(double x1, double eps1, double eps2, int max,
                         int ind[], Kansu kn1, Kansu kn2)
	{
		double g, dg, x;
		int sw;

		x      = x1;
		ind[0] = 0;
		sw     = 0;

		while (sw == 0 && ind[0] >= 0) {

			ind[0]++;
			sw = 1;
			g  = kn1.snx(x1);

			if (Math.abs(g) > eps2) {
				if (ind[0] <= max) {
					dg = kn2.snx(x1);
					if (Math.abs(dg) > eps2) {
						x = x1 - g / dg;
						if (Math.abs(x-x1) > eps1 && Math.abs(x-x1) > eps1*Math.abs(x)) {
							x1 = x;
							sw = 0;
						}
					}
					else
						ind[0] = -1;
				}
				else
					ind[0] = -1;
			}
		}

		return x;
	}

	/*************************************************************/
	/*     二分法による非線形方程式(f(x)=0)の解                  */
	/*          x1,x2 : 初期値                                   */
	/*          eps1 : 終了条件１（｜x(k+1)-x(k)｜＜eps1）       */
	/*          eps2 : 終了条件２（｜f(x(k))｜＜eps2）           */
	/*          max : 最大試行回数                               */
	/*          ind : 実際の試行回数                             */
	/*                （負の時は解を得ることができなかった）     */
	/*          kn : 関数値を計算するクラスオブジェクト          */
	/*          return : 解                                      */
	/*************************************************************/
	static double bisection(double x1, double x2, double eps1, double eps2, int max, int ind[], Kansu kn)
	{
		double f0, f1, f2, x0 = 0.0;
		int sw;

		f1 = kn.snx(x1);
		f2 = kn.snx(x2);

		if (f1*f2 > 0.0)
			ind[0] = -1;

		else {
			ind[0] = 0;
			if (f1*f2 == 0.0)
				x0 = (f1 == 0.0) ? x1 : x2;
			else {
				sw = 0;
				while (sw == 0 && ind[0] >= 0) {
					sw      = 1;
					ind[0] += 1;
					x0      = 0.5 * (x1 + x2);
					f0      = kn.snx(x0);
					if (Math.abs(f0) > eps2) {
						if (ind[0] <= max) {
							if (Math.abs(x1-x2) > eps1 && Math.abs(x1-x2) > eps1*Math.abs(x2)) {
								sw = 0;
								if (f0*f1 < 0.0) {
									x2 = x0;
									f2 = f0;
								}
								else {
									x1 = x0;
									f1 = f0;
								}
							}
						}
						else
							ind[0] = -1;
					}
				}
			}
		}

		return x0;
	}
}