行列の固有値（フレーム法＋ベアストウ法）

/********************************************/
/* 行列の固有値（フレーム法＋ベアストウ法） */
/*      coded by Y.Suganuma                 */
/********************************************/
import java.io.*;

public class Test {
	public static void main(String args[]) throws IOException
	{
		double A[][], H1[][], H2[][], a[], b[], c[], rl[], im[], p0, q0, eps;
		int i1, ind, ct, n;
					// データの設定
		ct   = 1000;
		eps  = 1.0e-10;
		p0   = 0.0;
		q0   = 0.0;
		n    = 3;
		a    = new double [n+1];
		b    = new double [n+1];
		c    = new double [n+1];
		rl   = new double [n];
		im   = new double [n];
		A    = new double [n][n];
		H1   = new double [n][n];
		H2   = new double [n][n];

		A[0][0] = 7.0;
		A[0][1] = 2.0;
		A[0][2] = 1.0;

		A[1][0] = 2.0;
		A[1][1] = 1.0;
		A[1][2] = -4.0;

		A[2][0] = 1.0;
		A[2][1] = -4.0;
		A[2][2] = 2.0;
					// 計算
		ind = Frame(n, ct, eps, p0, q0, a, b, c, rl, im, A, H1, H2);
					// 出力
		if (ind > 0)
			System.out.println("収束しませんでした！");
		else {
			for (i1 = 0; i1 < n; i1++)
				System.out.println("   " + rl[i1] + "  i " + im[i1]);
		}
	}

	/*************************************************/
	/* 行列の固有値（フレーム法＋ベアストウ法）      */
	/*      n : 次数                                 */
	/*      ct : 最大繰り返し回数                    */
	/*      eps : 収束判定条件                       */
	/*      p0, q0 : x2+px+qにおけるp,qの初期値      */
	/*      a : 係数（最高次から与え，値は変化する） */
	/*      b, c : 作業域（(n+1)次の配列）           */
	/*      rl, im : 結果の実部と虚部                */
	/*      A : 行列                                 */
	/*      H1, H2 : 作業域（nxnの行列）             */
	/*      return : =0 : 正常                       */
	/*               =1 : 収束せず                   */
	/*      coded by Y.Suganuma                      */
	/*************************************************/
	static int Frame(int n, int ct, double eps, double p0, double q0, double a[], double b[],
           double c[], double rl[], double im[], double A[][], double H1[][], double H2[][])
	{
		int i1, i2, i3, i4, ind;

		a[0] = 1.0;
					// a1の計算
		a[1] = 0.0;
		for (i1 = 0; i1 < n; i1++)
			a[1] -= A[i1][i1];
					// a2の計算
		for (i1 = 0; i1 < n; i1++) {
			for (i2 = 0; i2 < n; i2++)
				H1[i1][i2] = A[i1][i2];
			H1[i1][i1] += a[1];
		}

		a[2] = 0.0;
		for (i1 = 0; i1 < n; i1++) {
			for (i2 = 0; i2 < n; i2++)
				a[2] -= A[i1][i2] * H1[i2][i1];
		}
		a[2] *= 0.5;
					// a3からanの計算
		for (i1 = 3; i1 <= n; i1++) {

			for (i2 = 0; i2 < n; i2++) {
				for (i3 = 0; i3 < n; i3++) {
					H2[i2][i3] = 0.0;
					for (i4 = 0; i4 < n; i4++)
						H2[i2][i3] += A[i2][i4] * H1[i4][i3];
				}
				H2[i2][i2] += a[i1-1];
			}

			a[i1] = 0.0;
			for (i2 = 0; i2 < n; i2++) {
				for (i3 = 0; i3 < n; i3++)
					a[i1] -= A[i2][i3] * H2[i3][i2];
			}
			a[i1] /= i1;

			for (i2 = 0; i2 < n; i2++) {
				for (i3 = 0; i3 < n; i3++)
					H1[i2][i3] = H2[i2][i3];
			}
		}
					// ベアストウ法
		ind = Bairstow(n, ct, eps, p0, q0, a, b, c, rl, im, 0);

		return ind;
	}

	/*************************************************/
	/* 実係数代数方程式の解（ベアストウ法）          */
	/*      n : 次数                                 */
	/*      ct : 最大繰り返し回数                    */
	/*      eps : 収束判定条件                       */
	/*      p0, q0 : x2+px+qにおけるp,qの初期値      */
	/*      a : 係数（最高次から与え，値は変化する） */
	/*      b,c : 作業域（(n+1)次の配列）            */
	/*      rl, im : 結果の実部と虚部                */
	/*      p : 答えの位置                           */
	/*      return : =0 : 正常                       */
	/*               =1 : 収束せず                   */
	/*      coded by Y.Suganuma                      */
	/*************************************************/
	static int Bairstow(int n, int ct, double eps, double p0, double q0, 
	                    double a[], double b[], double c[], double rl[], double im[], int p)
	{
		double D, dp, dq, p1 = p0, p2 = 0.0, q1 = q0, q2 = 0.0;
		int i1, ind = 0, count = 0;
	/*
	          １次の場合
	*/
		if (n == 1) {
			if (Math.abs(a[0]) < eps)
				ind = 1;
			else {
				rl[p] = -a[1] / a[0];
				im[p] = 0.0;
			}
		}
	/*
	          ２次の場合
	*/
		else if (n == 2) {
						// １次式
			if (Math.abs(a[0]) < eps) {
				if (Math.abs(a[1]) < eps)
					ind = 1;
				else {
					rl[p] = -a[2] / a[1];
					im[p] = 0.0;
				}
			}
						// ２次式
			else {
				D = a[1] * a[1] - 4.0 * a[0] * a[2];
				if (D < 0.0) {   // 虚数
					D        = Math.sqrt(-D);
					a[0]    *= 2.0;
					rl[p]    = -a[1] / a[0];
					rl[p+1]  = -a[1] / a[0];
					im[p]    = D / a[0];
					im[p+1]  = -im[p];
				}
				else {           // 実数
					D        = Math.sqrt(D);
					a[0]     = 1.0 / (2.0 * a[0]);
					rl[p]    = a[0] * (-a[1] + D);
					rl[p+1]  = a[0] * (-a[1] - D);
					im[p]    = 0.0;
					im[p+1]  = 0.0;
				}
			}
		}
						// ３次以上の場合
		else {
							// 因数分解
			ind = 1;
			while (ind > 0 && count <= ct) {
				for (i1 = 0; i1 <= n; i1++) {
					if (i1 == 0)
						b[i1] = a[i1];
					else if (i1 == 1)
						b[i1] = a[i1] - p1 * b[i1-1];
					else
						b[i1] = a[i1] - p1 * b[i1-1] - q1 * b[i1-2];
				}
				for (i1 = 0; i1 <= n; i1++) {
					if (i1 == 0)
						c[i1] = b[i1];
					else if (i1 == 1)
						c[i1] = b[i1] - p1 * c[i1-1];
					else
						c[i1] = b[i1] - p1 * c[i1-1] - q1 * c[i1-2];
				}
				D = c[n-2] * c[n-2] - c[n-3] * (c[n-1] - b[n-1]);
				if (Math.abs(D) < eps)
					return ind;
				else {
					dp = (b[n-1] * c[n-2] - b[n] * c[n-3]) / D;
					dq = (b[n] * c[n-2] - b[n-1] * (c[n-1] - b[n-1])) / D;
					p2 = p1 + dp;
					q2 = q1 + dq;
					if (Math.abs(dp) < eps && Math.abs(dq) < eps)
						ind = 0;
					else {
						count++;
						p1 = p2;
						q1 = q2;
					}
				}
			}
	
			if (ind == 0) {
							// ２次方程式を解く
				D = p2 * p2 - 4.0 * q2;
				if (D < 0.0) {   // 虚数
					D        = Math.sqrt(-D);
					rl[p]    = -0.5 * p2;
					rl[p+1]  = -0.5 * p2;
					im[p]    = 0.5 * D;
					im[p+1]  = -im[p];
				}
				else {           // 実数
					D        = Math.sqrt(D);
					rl[p]    = 0.5 * (-p2 + D);
					rl[p+1]  = 0.5 * (-p2 - D);
					im[p]    = 0.0;
					im[p+1]  = 0.0;
				}
							// 残りの方程式を解く
				n -= 2;
				for (i1 = 0; i1 <= n; i1++)
					a[i1] = b[i1];
				ind = Bairstow(n, ct, eps, p0, q0, a, b, c, rl, im, p+2);
			}
		}
	
		return ind;
	}
}