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

public class F {

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

	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("自由度１は？ ");
		dof1 = Integer.parseInt(in.readLine());
		System.out.print("自由度２は？ ");
		dof2 = 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.f(x, z, dof1, dof2);
				System.out.println("P(X = " + x + ") = " + z[0] + ",  P( X < " + x + ") = " + y +
                                   " (自由度 = " + dof1 + " " + dof2 + ")");
			}
					// グラフ出力
			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.f(x, z, dof1, dof2);
					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_f(sw1);
			System.out.println(x + "％値 = " + y + "  sw " + sw1[0] +
                               " (自由度 = " + dof1 + " " + dof2 + ")");
		}
	}
}

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

/********************/
/*     関数値の計算 */
/********************/
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 - F.p - App.normal(x, w);
				break;
						// 関数値（f(x)）の計算（Ｆ分布）
			case 1:
				y = App.f(x, w, F.dof1, F.dof2) - 1.0 + F.p;
				break;
						// 関数の微分の計算（Ｆ分布）
			case 2:
				y = App.f(x, w, F.dof1, F.dof2);
				y = w[0];
				break;
		}

		return y;
	}
}

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

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

	/********************************************/
	/*     ｆ分布の計算（P(X = ff), P(X < ff)） */
	/*          dd : P(X = ff)                  */
	/*          df1,df2 : 自由度                */
	/*          return : P(X < ff)              */
	/********************************************/
	static double f(double ff, double dd[], int df1, int df2)
	{
		double pi = Math.PI;
		double pp, u, x;
		int ia, ib, i1;

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

		x  = ff * df1 / (ff * df1 + df2);

		if(df1%2 == 0) {
			if (df2%2 == 0) {
				u  = x * (1.0 - x);
				pp = x;
				ia = 2;
				ib = 2;
			}
			else {
				u  = 0.5 * x * Math.sqrt(1.0-x);
				pp = 1.0 - Math.sqrt(1.0-x);
				ia = 2;
				ib = 1;
			}
		}

		else {
			if (df2%2 == 0) {
				u  = 0.5 * Math.sqrt(x) * (1.0 - x);
				pp = Math.sqrt(x);
				ia = 1;
				ib = 2;
			}
			else {
				u  = Math.sqrt(x*(1.0-x)) / pi;
				pp = 1.0 - 2.0 * Math.atan2(Math.sqrt(1.0-x), Math.sqrt(x)) / pi;
				ia = 1;
				ib = 1;
			}
		}

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

		if (ib != df2) {
			for (i1 = ib; i1 <= df2-2; i1 += 2) {
				pp += 2.0 * u / i1;
				u  *= (1.0 - x) * (i1 + df1) / i1;
			}
		}

		dd[0] = u / ff;

		return pp;
	}

	/********************************************/
	/*     ｆ分布のｐ％値（P(X > u) = 0.01p）   */
	/*          ind : >= 0 : normal（収束回数） */
	/*                = -1 : 収束しなかった     */
	/********************************************/
	static double p_f(int ind[])
	{
		double a = 0.0, a1 = 0.0, b = 0.0, b1 = 0.0, df1 = 0.0, df2 = 0.0,
               e = 0.0, ff = 0.0, f0, x, y1, y2, yq = 0.0;
		int sw = 0, MAX = 340;
	/*
	     初期値計算の準備
	          while文は，大きな自由度によるガンマ関数の
	          オーバーフローを避けるため
	*/
		while (sw >= 0) {

			df1 = 0.5 * (F.dof1 - sw);
			df2 = 0.5 * F.dof2;
			a   = 2.0 / (9.0 * (F.dof1 - sw));
			a1  = 1.0 - a;
			b   = 2.0 / (9.0 * F.dof2);
			b1  = 1.0 - b;

			yq  = p_normal(ind);

			e   = b1 * b1 - b * yq * yq;

			if (e > 0.8 || (F.dof1+F.dof2-sw) <= MAX)
				sw = -1;
			else {
				sw += 1;
				if ((F.dof1-sw) == 0)
					sw = -2;
			}
		}

		if (sw == -2)
			ind[0] = -1;
	/*
	     f0 : 初期値
	*/
		else {
			if (e > 0.8) {
				x  = (a1 * b1 + yq * Math.sqrt(a1*a1*b+a*e)) / e;
				f0 = Math.pow(x, 3.0);
			}
			else {
				y1 = Math.pow((double)F.dof2, df2-1.0);
				y2 = Math.pow((double)F.dof1, df2);
				x  = gamma(df1+df2, ind) / gamma(df1, ind) / gamma(df2, ind) *
                     2.0 * y1 / y2 / F.p;
				f0 = Math.pow(x, 2.0/F.dof2);
			}
	/*
	     ニュートン法
	*/
			Kansu kn1 = new Kansu(1);
			Kansu kn2 = new Kansu(2);

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

		return ff;
	}

	/********************************************/
	/*     Γ（ｘ）の計算（ガンマ関数，近似式） */
	/*          ier : =0 : normal               */
	/*                =-1 : x=-n (n=0,1,2,･･･)  */
	/*          return : 結果                   */
	/********************************************/
	static double gamma(double x, int ier[])
	{
		double err, g, s, t, v, w, y;
		int k;

		ier[0] = 0;

		if (x > 5.0) {
			v = 1.0 / x;
			s = ((((((-0.000592166437354 * v + 0.0000697281375837) * v +
                0.00078403922172) * v - 0.000229472093621) * v -
                0.00268132716049) * v + 0.00347222222222) * v +
                0.0833333333333) * v + 1.0;
			g = 2.506628274631001 * Math.exp(-x) * Math.pow(x,x-0.5) * s;
		}

		else {

			err = 1.0e-20;
			w   = x;
			t   = 1.0;

			if (x < 1.5) {

				if (x < err) {
					k = (int)x;
					y = (double)k - x;
					if (Math.abs(y) < err || Math.abs(1.0-y) < err)
						ier[0] = -1;
				}

				if (ier[0] == 0) {
					while (w < 1.5) {
						t /= w;
						w += 1.0;
					}
				}
			}

			else {
				if (w > 2.5) {
					while (w > 2.5) {
						w -= 1.0;
						t *= w;
					}
				}
			}

			w -= 2.0;
			g  = (((((((0.0021385778 * w - 0.0034961289) * w + 
                 0.0122995771) * w - 0.00012513767) * w + 0.0740648982) * w +
                 0.0815652323) * w + 0.411849671) * w + 0.422784604) * w +
                 0.999999926;
			g *= t;
		}

		return g;
	}

	/*****************************************************/
	/*     標準正規分布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;
	}
}