情報学部 菅沼ホーム 目次 索引

行列の固有値(フレーム法+ベアストウ法)

    1. A. C++
    2. B. Java
    3. C. JavaScript
    4. D. PHP
    5. E. Ruby
    6. F. Python
    7. G. C#
    8. H. VB

  プログラムは,行列の固有値をフレーム法とベアストウ法によって求めるためのものです.

  1. C++

    /********************************************/
    /* 行列の固有値(フレーム法+ベアストウ法) */
    /*      coded by Y.Suganuma                 */
    /********************************************/
    #include <stdio.h>
    
    int Bairstow(int, int, double, double, double,
                 double *, double *, double *, double *, double *);
    int Frame(int, int, double, double, double, double *, double *, double *,
              double *, double *, double **, double **, double **);
    
    int main()
    {
    	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];
    	H1   = new double * [n];
    	H2   = new double * [n];
    	for (i1 = 0; i1 < n; i1++) {
    		A[i1]  = new double [n];
    		H1[i1] = new double [n];
    		H2[i1] = new double [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)
    		printf("収束しませんでした!\n");
    	else {
    		for (i1 = 0; i1 < n; i1++)
    			printf("   %f  i %f\n", rl[i1], im[i1]);
    	}
    
    	delete [] a;
    	delete [] b;
    	delete [] c;
    	delete [] rl;
    	delete [] im;
    	for (i1 = 0; i1 < n; i1++) {
    		delete [] A[i1];
    		delete [] H1[i1];
    		delete [] H2[i1];
    	}
    	delete [] A;
    	delete [] H1;
    	delete [] H2;
    
    	return 0;
    }
    
    /*************************************************/
    /* 行列の固有値(フレーム法+ベアストウ法)      */
    /*      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                      */
    /*************************************************/
    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);
    
    	return ind;
    }
    
    /*************************************************/
    /* 実係数代数方程式の解(ベアストウ法)          */
    /*      n : 次数                                 */
    /*      ct : 最大繰り返し回数                    */
    /*      eps : 収束判定条件                       */
    /*      p0, q0 : x2+px+qにおけるp,qの初期値      */
    /*      a : 係数(最高次から与え,値は変化する) */
    /*      b, c : 作業域((n+1)次の配列)           */
    /*      rl, im : 結果の実部と虚部                */
    /*      return : =0 : 正常                       */
    /*               =1 : 収束せず                   */
    /*      coded by Y.Suganuma                      */
    /*************************************************/
    #include <math.h>
    
    int Bairstow(int n, int ct, double eps, double p0, double q0, 
                 double *a, double *b, double *c, double *rl, double *im)
    {
    	double D, dp, dq, p1 = p0, p2 = 0.0, q1 = q0, q2 = 0.0;
    	int i1, ind = 0, count = 0;
    /*
              1次の場合
    */
    	if (n == 1) {
    		if (fabs(a[0]) < eps)
    			ind = 1;
    		else {
    			rl[0] = -a[1] / a[0];
    			im[0] = 0.0;
    		}
    	}
    /*
              2次の場合
    */
    	else if (n == 2) {
    					// 1次式
    		if (fabs(a[0]) < eps) {
    			if (fabs(a[1]) < eps)
    				ind = 1;
    			else {
    				rl[0] = -a[2] / a[1];
    				im[0] = 0.0;
    			}
    		}
    					// 2次式
    		else {
    			D = a[1] * a[1] - 4.0 * a[0] * a[2];
    			if (D < 0.0) {   // 虚数
    				D      = sqrt(-D);
    				a[0]  *= 2.0;
    				rl[0]  = -a[1] / a[0];
    				rl[1]  = -a[1] / a[0];
    				im[0]  = D / a[0];
    				im[1]  = -im[0];
    			}
    			else {           // 実数
    				D      = sqrt(D);
    				a[0]   = 1.0 / (2.0 * a[0]);
    				rl[0]  = a[0] * (-a[1] + D);
    				rl[1]  = a[0] * (-a[1] - D);
    				im[0]  = 0.0;
    				im[1]  = 0.0;
    			}
    		}
    	}
    					// 3次以上の場合
    	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 (fabs(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 (fabs(dp) < eps && fabs(dq) < eps)
    					ind = 0;
    				else {
    					count++;
    					p1 = p2;
    					q1 = q2;
    				}
    			}
    		}
    
    		if (ind == 0) {
    						// 2次方程式を解く
    			D = p2 * p2 - 4.0 * q2;
    			if (D < 0.0) {   // 虚数
    				D      = sqrt(-D);
    				rl[0]  = -0.5 * p2;
    				rl[1]  = -0.5 * p2;
    				im[0]  = 0.5 * D;
    				im[1]  = -im[0];
    			}
    			else {           // 実数
    				D      = sqrt(D);
    				rl[0]  = 0.5 * (-p2 + D);
    				rl[1]  = 0.5 * (-p2 - D);
    				im[0]  = 0.0;
    				im[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[2], &im[2]);
    		}
    	}
    
    	return ind;
    }
    			

  2. Java

    /********************************************/
    /* 行列の固有値(フレーム法+ベアストウ法) */
    /*      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;
    	/*
    	          1次の場合
    	*/
    		if (n == 1) {
    			if (Math.abs(a[0]) < eps)
    				ind = 1;
    			else {
    				rl[p] = -a[1] / a[0];
    				im[p] = 0.0;
    			}
    		}
    	/*
    	          2次の場合
    	*/
    		else if (n == 2) {
    						// 1次式
    			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;
    				}
    			}
    						// 2次式
    			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;
    				}
    			}
    		}
    						// 3次以上の場合
    		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) {
    							// 2次方程式を解く
    				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;
    	}
    }
    			

  3. JavaScript

      ここをクリックすると,任意のデータに対して,画面上で解を得ることができます.
    <!DOCTYPE HTML>
    
    <HTML>
    
    <HEAD>
    
    	<TITLE>行列の固有値(フレーム法+ベアストウ法)</TITLE>
    	<META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=utf-8">
    	<SCRIPT TYPE="text/javascript">
    		function main()
    		{
    					// データの設定
    			let ct  = parseInt(document.getElementById("trial").value);
    			let eps = 1.0e-10;
    			let p0  = parseFloat(document.getElementById("p0").value);
    			let q0  = parseFloat(document.getElementById("q0").value);
    			let n   = parseInt(document.getElementById("order").value);
    			let a   = new Array();
    			let b   = new Array();
    			let c   = new Array();
    			let rl  = new Array();
    			let im  = Array();
    			let A   = new Array();
    			for (let i1 = 0; i1 < n; i1++)
    				A[i1] = new Array();
    			let s = (document.getElementById("ar").value).split(/ {1,}|\n{1,}/);
    			let k = 0;
    			for (let i1 = 0; i1 < n; i1++) {
    				for (let i2 = 0; i2 < n; i2++) {
    					A[i1][i2] = parseFloat(s[k]);
    					k++;
    				}
    			}
    			let H1 = new Array();
    			for (let i1 = 0; i1 < n; i1++)
    				H1[i1] = new Array();
    			let H2 = new Array();
    			for (let i1 = 0; i1 < n; i1++)
    				H2[i1] = new Array();
    
    			ind = frame(n, ct, eps, p0, q0, a, b, c, rl, im, A, H1, H2);
    					// 出力
    			if (ind > 0)
    				document.getElementById("ans").value = "解を求めることができません!";
    			else {
    				let str = "";
    				for (let i1 = 0; i1 < n; i1++)
    					str = str + rl[i1] + "  i " + im[i1] + "\n";
    				document.getElementById("ans").value = str;
    			}
    		}
    
    		/*************************************************/
    		/* 行列の固有値(フレーム法+ベアストウ法)      */
    		/*      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                      */
    		/*************************************************/
    		function frame(n, ct, eps, p0, q0, a, b, c, rl, im, A, H1, H2)
    		{
    			let i1;
    			let i2;
    			let i3;
    			let i4;
    			let ind;
    
    			a[0] = 1.0;
    						// b1の計算
    			a[1] = 0.0;
    			for (i1 = 0; i1 < n; i1++)
    				a[1] -= A[i1][i1];
    						// b2の計算
    			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;
    						// b3からbnの計算
    			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                      */
    		/*************************************************/
    		function bairstow(n, ct, eps, p0, q0, a, b, c, rl, im, p)
    		{
    			let D;
    			let dp;
    			let dq;
    			let p1 = p0;
    			let p2 = 0.0;
    			let q1 = q0;
    			let q2 = 0.0;
    			let i1;
    			let ind = 0;
    			let count = 0;
    		/*
    		          1次の場合
    		*/
    			if (n == 1) {
    				if (Math.abs(a[0]) < eps)
    					ind = 1;
    				else {
    					rl[p] = -a[1] / a[0];
    					im[p] = 0.0;
    				}
    			}
    		/*
    		          2次の場合
    		*/
    			else if (n == 2) {
    							// 1次式
    				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;
    					}
    				}
    							// 2次式
    				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;
    					}
    				}
    			}
    						// 3次以上の場合
    			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) {
    								// 2次方程式を解く
    					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;
    		}
    	</SCRIPT>
    
    </HEAD>
    
    <BODY STYLE="font-size: 130%; background-color: #eeffee;">
    
    	<H2 STYLE="text-align:center"><B>行列の固有値(フレーム法+ベアストウ法)</B></H2>
    
    	<DL>
    		<DT>  テキストフィールドおよびテキストエリアには,例として,以下に示す行列の固有値を求める場合に対する値が設定されています.他の問題を実行する場合は,それらを適切に修正してください.
    		<P STYLE="text-align:center"><IMG SRC="eigen.gif"></P>
    	</DL>
    
    	<DIV STYLE="text-align:center">
    		次数:<INPUT ID="order" STYLE="font-size: 100%" TYPE="text" SIZE="2" VALUE="3"> 
    		p0:<INPUT ID="p0" STYLE="font-size: 100%;" TYPE="text" SIZE="2" VALUE="0"> 
    		q0:<INPUT ID="q0" STYLE="font-size: 100%;" TYPE="text" SIZE="2" VALUE="0"> 
    		最大繰り返し回数:<INPUT ID="trial" STYLE="font-size: 100%;" TYPE="text" SIZE="4" VALUE="1000"> 
    		<BUTTON STYLE="font-size: 100%; background-color: pink" onClick="main()">OK</BUTTON><BR><BR>
    		<TEXTAREA ID="ar" COLS="50" ROWS="15" STYLE="font-size: 100%">7 2 1
    2 1 -4
    1 -4 2</TEXTAREA><BR><BR>
    		<TEXTAREA ID="ans" COLS="50" ROWS="15" STYLE="font-size: 100%"></TEXTAREA>
    	</DIV>
    
    </BODY>
    
    </HTML>
    			

  4. PHP

    <?php
    
    /********************************************/
    /* 行列の固有値(フレーム法+ベアストウ法) */
    /*      coded by Y.Suganuma                 */
    /********************************************/
    
    					// データの設定
    	$ct  = 1000;
    	$eps = 1.0e-10;
    	$p0  = 0.0;
    	$q0  = 0.0;
    	$n   = 3;
    	$a   = array($n+1);
    	$b   = array($n+1);
    	$c   = array($n+1);
    	$rl  = array($n);
    	$im  = array($n);
    	$A   = array($n);
    	$H1  = array($n);
    	$H2  = array($n);
    	for ($i1 = 0; $i1 < $n; $i1++) {
    		$A[$i1]  = array($n);
    		$H1[$i1] = array($n);
    		$H2[$i1] = array($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)
    		printf("収束しませんでした!\n");
    	else {
    		for ($i1 = 0; $i1 < $n; $i1++)
    			printf("   %f  i %f\n", $rl[$i1], $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                      */
    /*************************************************/
    function Frame($n, $ct, $eps, $p0, $q0, $a, $b, $c, &$rl, &$im, $A, $H1, $H2)
    {
    	$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 : 結果の実部と虚部                */
    /*      k : 結果を設定する配列の位置             */
    /*      return : =0 : 正常                       */
    /*               =1 : 収束せず                   */
    /*      coded by Y.Suganuma                      */
    /*************************************************/
    function Bairstow($n, $ct, $eps, $p0, $q0, $a, $b, $c, &$rl, &$im, $k)
    {
    	
    	$p1    = $p0;
    	$p2    = 0.0;
    	$q1    = $q0;
    	$q2    = 0.0;
    	$ind   = 0;
    	$count = 0;
    /*
              1次の場合
    */
    	if ($n == 1) {
    		if (abs($a[0]) < $eps)
    			$ind = 1;
    		else {
    			$rl[$k] = -$a[1] / $a[0];
    			$im[$k] = 0.0;
    		}
    	}
    /*
              2次の場合
    */
    	else if ($n == 2) {
    					// 1次式
    		if (abs($a[0]) < $eps) {
    			if (abs($a[1]) < $eps)
    				$ind = 1;
    			else {
    				$rl[$k] = -$a[2] / $a[1];
    				$im[$k] = 0.0;
    			}
    		}
    					// 2次式
    		else {
    			$D = $a[1] * $a[1] - 4.0 * $a[0] * $a[2];
    			if ($D < 0.0) {   // 虚数
    				$D         = sqrt(-$D);
    				$a[0]     *= 2.0;
    				$rl[$k]    = -$a[1] / $a[0];
    				$rl[$k+1]  = -$a[1] / $a[0];
    				$im[$k]    = $D / $a[0];
    				$im[$k+1]  = -$im[$k];
    			}
    			else {           // 実数
    				$D        = sqrt($D);
    				$a[0]     = 1.0 / (2.0 * $a[0]);
    				$rl[$k]   = $a[0] * (-$a[1] + $D);
    				$rl[$k+1] = $a[0] * (-$a[1] - $D);
    				$im[$k]   = 0.0;
    				$im[$k+1] = 0.0;
    			}
    		}
    	}
    					// 3次以上の場合
    	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 (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 (abs($dp) < $eps && abs($dq) < $eps)
    					$ind = 0;
    				else {
    					$count++;
    					$p1 = $p2;
    					$q1 = $q2;
    				}
    			}
    		}
    
    		if ($ind == 0) {
    						// 2次方程式を解く
    			$D = $p2 * $p2 - 4.0 * $q2;
    			if ($D < 0.0) {   // 虚数
    				$D        = sqrt(-$D);
    				$rl[$k]   = -0.5 * $p2;
    				$rl[$k+1] = -0.5 * $p2;
    				$im[$k]   = 0.5 * $D;
    				$im[$k+1] = -$im[$k];
    			}
    			else {           // 実数
    				$D        = sqrt($D);
    				$rl[$k]   = 0.5 * (-$p2 + $D);
    				$rl[$k+1] = 0.5 * (-$p2 - $D);
    				$im[$k]   = 0.0;
    				$im[$k+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, $k+2);
    		}
    	}
    
    	return $ind;
    }
    
    ?>
    			

  5. Ruby

    #*******************************************/
    # 行列の固有値(フレーム法+ベアストウ法) */
    #      coded by Y.Suganuma                 */
    #*******************************************/
    
    #************************************************/
    # 行列の固有値(フレーム法+ベアストウ法)      */
    #      n : 次数                                 */
    #      ct : 最大繰り返し回数                    */
    #      eps : 収束判定条件                       */
    #      p0, q0 : x2+px+qにおけるp,qの初期値      */
    #      a : 係数(最高次から与え,値は変化する) */
    #      b, c : 作業域((n+1)次の配列)           */
    #      rl, im : 結果の実部と虚部                */
    #      aa : 行列                                */
    #      h1, h2 : 作業域(nxnの行列)             */
    #      return : =0 : 正常                       */
    #               =1 : 収束せず                   */
    #      coded by Y.Suganuma                      */
    #************************************************/
    def Frame(n, ct, eps, p0, q0, a, b, c, rl, im, aa, h1, h2)
    	a[0] = 1.0
    					# a1の計算
    	a[1] = 0.0
    	for i1 in 0 ... n
    		a[1] -= aa[i1][i1]
    	end
    					# a2の計算
    	for i1 in 0 ... n
    		for i2 in 0 ... n
    			h1[i1][i2] = aa[i1][i2]
    		end
    		h1[i1][i1] += a[1]
    	end
    
    	a[2] = 0.0
    	for i1 in 0 ... n
    		for i2 in 0 ... n
    			a[2] -= aa[i1][i2] * h1[i2][i1]
    		end
    	end
    	a[2] *= 0.5
    					# a3からanの計算
    	for i1 in 3 ... n+1
    
    		for i2 in 0 ... n
    			for i3 in 0 ... n
    				h2[i2][i3] = 0.0
    				for i4 in 0 ... n
    					h2[i2][i3] += aa[i2][i4] * h1[i4][i3]
    				end
    			end
    			h2[i2][i2] += a[i1-1]
    		end
    
    		a[i1] = 0.0
    		for i2 in 0 ... n
    			for i3 in 0 ... n
    				a[i1] -= aa[i2][i3] * h2[i3][i2]
    			end
    		end
    		a[i1] /= i1
    
    		for i2 in 0 ... n
    			for i3 in 0 ... n
    				h1[i2][i3] = h2[i2][i3]
    			end
    		end
    	end
    					# ベアストウ法
    	ind = Bairstow(n, ct, eps, p0, q0, a, b, c, rl, im, 0)
    
    	return ind
    end
    
    #************************************************/
    # 実係数代数方程式の解(ベアストウ法)          */
    #      n : 次数                                 */
    #      ct : 最大繰り返し回数                    */
    #      eps : 収束判定条件                       */
    #      p0, q0 : x2+px+qにおけるp,qの初期値      */
    #      a : 係数(最高次から与え,値は変化する) */
    #      b,c : 作業域((n+1)次の配列)            */
    #      rl, im : 結果の実部と虚部                */
    #      k : 結果の位置                           */
    #      return : =0 : 正常                       */
    #               =1 : 収束せず                   */
    #      coded by Y.Suganuma                      */
    #************************************************/
    def Bairstow(n, ct, eps, p0, q0, a, b, c, rl, im, k)
    			# 初期設定
    	p1    = p0
    	p2    = 0.0
    	q1    = q0
    	q2    = 0.0
    	ind   = 0
    	count = 0
    #
    #          1次の場合
    #
    	if n == 1
    		if a[0].abs() < eps
    			ind = 1
    		else
    			rl[k] = -a[1] / a[0]
    			im[k] = 0.0
    		end
    #
    #          2次の場合
    #
    	elsif n == 2
    					# 1次式
    		if a[0].abs() < eps
    			if a[1].abs() < eps
    				ind = 1
    			else
    				rl[k] = -a[2] / a[1]
    				im[k] = 0.0
    			end
    					# 2次式
    		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[k]    = -a[1] / a[0]
    				rl[k+1]  = -a[1] / a[0]
    				im[k]    = d / a[0]
    				im[k+1]  = -im[0]
    			else           # 実数
    				d       = Math.sqrt(d)
    				a[0]    = 1.0 / (2.0 * a[0])
    				rl[k]   = a[0] * (-a[1] + d)
    				rl[k+1] = a[0] * (-a[1] - d)
    				im[k]   = 0.0
    				im[k+1] = 0.0
    			end
    		end
    					# 3次以上の場合
    	else
    						# 因数分解
    		ind = 1
    		while ind > 0 && count <= ct
    			for i1 in 0 ... n+1
    				if i1 == 0
    					b[i1] = a[i1]
    				elsif i1 == 1
    					b[i1] = a[i1] - p1 * b[i1-1]
    				else
    					b[i1] = a[i1] - p1 * b[i1-1] - q1 * b[i1-2]
    				end
    			end
    			for i1 in 0 ... n+1
    				if i1 == 0
    					c[i1] = b[i1]
    				elsif i1 == 1
    					c[i1] = b[i1] - p1 * c[i1-1]
    				else
    					c[i1] = b[i1] - p1 * c[i1-1] - q1 * c[i1-2]
    				end
    			end
    			d = c[n-2] * c[n-2] - c[n-3] * (c[n-1] - b[n-1])
    			if d.abs() < 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 dp.abs() < eps && dq.abs() < eps
    					ind = 0
    				else
    					count += 1
    					p1 = p2
    					q1 = q2
    				end
    			end
    		end
    
    		if ind == 0
    						# 2次方程式を解く
    			d = p2 * p2 - 4.0 * q2
    			if d < 0.0   # 虚数
    				d       = Math.sqrt(-d)
    				rl[k]   = -0.5 * p2
    				rl[k+1] = -0.5 * p2
    				im[k]   = 0.5 * d
    				im[k+1] = -im[k]
    			else           # 実数
    				d       = Math.sqrt(d)
    				rl[k]   = 0.5 * (-p2 + d)
    				rl[k+1] = 0.5 * (-p2 - d)
    				im[k]   = 0.0
    				im[k+1] = 0.0
    			end
    						# 残りの方程式を解く
    			n -= 2
    			for i1 in 0 ... n+1
    				a[i1] = b[i1]
    			end
    			ind = Bairstow(n, ct, eps, p0, q0, a, b, c, rl, im, k+2)
    		end
    	end
    
    	return ind
    end
    
    					# データの設定
    ct   = 1000
    eps  = 1.0e-10
    p0   = 0.0
    q0   = 0.0
    n    = 3
    a    = Array.new(n+1)
    b    = Array.new(n+1)
    c    = Array.new(n+1)
    rl   = Array.new(n)
    im   = Array.new(n)
    aa   = Array.new(n)
    h1   = Array.new(n)
    h2   = Array.new(n)
    for i1 in 0 ... n
    	aa[i1] = Array.new(n)
    	h1[i1] = Array.new(n)
    	h2[i1] = Array.new(n)
    end
    
    aa[0][0] = 7.0
    aa[0][1] = 2.0
    aa[0][2] = 1.0
    
    aa[1][0] = 2.0
    aa[1][1] = 1.0
    aa[1][2] = -4.0
    
    aa[2][0] = 1.0
    aa[2][1] = -4.0
    aa[2][2] = 2.0
    				# 計算
    ind = Frame(n, ct, eps, p0, q0, a, b, c, rl, im, aa, h1, h2)
    				# 出力
    if ind > 0
    	printf("収束しませんでした!\n")
    else
    	for i1 in 0 ... n
    		printf("   %f  i %f\n", rl[i1], im[i1])
    	end
    end
    			

  6. Python

    # -*- coding: UTF-8 -*-
    import numpy as np
    from math import *
    
    ############################################
    # 行列の固有値(フレーム法+ベアストウ法)
    #      n : 次数
    #      ct : 最大繰り返し回数
    #      eps : 収束判定条件
    #      p0, q0 : x2+px+qにおけるp,qの初期値
    #      a : 係数(最高次から与え,値は変化する)
    #      b, c : 作業域((n+1)次の配列)
    #      r : 結果
    #      A : 行列
    #      H1, H2 : 作業域(nxnの行列)
    #      return : =0 : 正常
    #               =1 : 収束せず 
    #      coded by Y.Suganuma
    ############################################
    
    def Frame(n, ct, eps, p0, q0, a, b, c, r, A, H1, H2) :
    
    	a[0] = 1.0
    			# a1の計算
    	a[1] = 0.0
    	for i1 in range(0, n) :
    		a[1] -= A[i1][i1]
    			# a2の計算
    	for i1 in range(0, n) :
    		for i2 in range(0, n) :
    			H1[i1][i2] = A[i1][i2]
    		H1[i1][i1] += a[1]
    
    	a[2] = 0.0
    	for i1 in range(0, n) :
    		for i2 in range(0, n) :
    			a[2] -= A[i1][i2] * H1[i2][i1]
    	a[2] *= 0.5
    			# a3からanの計算
    	for i1 in range(3, n+1) :
    
    		for i2 in range(0, n) :
    			for i3 in range(0, n) :
    				H2[i2][i3] = 0.0
    				for i4 in range(0, n) :
    					H2[i2][i3] += A[i2][i4] * H1[i4][i3]
    			H2[i2][i2] += a[i1-1]
    
    		a[i1] = 0.0
    		for i2 in range(0, n) :
    			for i3 in range(0, n) :
    				a[i1] -= A[i2][i3] * H2[i3][i2]
    		a[i1] /= i1
    
    		for i2 in range(0, n) :
    			for i3 in range(0, n) :
    				H1[i2][i3] = H2[i2][i3]
    			# ベアストウ法
    	ind = Bairstow(n, ct, eps, p0, q0, a, b, c, r, 0)
    
    	return ind
    
    ############################################
    # 実係数代数方程式の解(ベアストウ法)
    #      n : 次数
    #      ct : 最大繰り返し回数
    #      eps : 収束判定条件
    #      p0, q0 : x2+px+qにおけるp,qの初期値
    #      a : 係数(最高次から与え,値は変化する)
    #      b,c : 作業域((n+1)次の配列)
    #      r : 結果
    #      k : 結果の位置
    #      return : =0 : 正常
    #               =1 : 収束せず
    #      coded by Y.Suganuma
    ############################################
    
    def Bairstow(n, ct, eps, p0, q0, a, b, c, r, k) :
    
    	p1    = p0
    	p2    = 0.0
    	q1    = q0
    	q2    = 0.0
    	ind   = 0
    	count = 0
    			# 1次の場合
    	if n == 1 :
    		if abs(a[0]) < eps :
    			ind = 1
    		else :
    			r[k] = complex(-a[1] / a[0], 0)
    			# 2次の場合
    	elif n == 2 :
    				# 1次式
    		if abs(a[0]) < eps :
    			if abs(a[1]) < eps :
    				ind = 1
    			else :
    				r[k] = complex(-a[2] / a[1], 0)
    				# 2次式
    		else :
    			D = a[1] * a[1] - 4.0 * a[0] * a[2]
    			if D < 0.0 :   # 虚数
    				D      = sqrt(-D)
    				a[0]  *= 2.0
    				r[k]   = complex(-a[1] / a[0], D / a[0])
    				r[k+1] = complex(-a[1] / a[0], -D / a[0])
    			else :   # 実数
    				D      = sqrt(D)
    				a[0]   = 1.0 / (2.0 * a[0])
    				r[k]   = complex(a[0] * (-a[1] + D), 0)
    				r[k+1] = complex(a[0] * (-a[1] - D), 0)
    				# 3次以上の場合
    	else :
    					# 因数分解
    		ind = 1
    		while ind > 0 and count <= ct :
    			for i1 in range(0, n+1) :
    				if i1 == 0 :
    					b[i1] = a[i1]
    				elif 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 in range(0, n+1) :
    				if i1 == 0 :
    					c[i1] = b[i1]
    				elif 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 fabs(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 abs(dp) < eps and fabs(dq) < eps :
    					ind = 0
    				else :
    					count += 1
    					p1 = p2
    					q1 = q2
    
    		if ind == 0 :
    					# 2次方程式を解く
    			D = p2 * p2 - 4.0 * q2
    			if D < 0.0 :   # 虚数
    				D      = sqrt(-D)
    				r[k]   = complex(-0.5 * p2, 0.5 * D)
    				r[k+1] = complex(-0.5 * p2, -0.5 * D)
    			else :   # 実数
    				D      = sqrt(D)
    				r[k]   = complex(0.5 * (-p2 + D), 0)
    				r[k+1] = complex(0.5 * (-p2 - D), 0)
    					# 残りの方程式を解く
    			n -= 2
    			for i1 in range(0, n+1) :
    				a[i1] = b[i1]
    			ind = Bairstow(n, ct, eps, p0, q0, a, b, c, r, k+2)
    
    	return ind
    
    ############################################
    # 行列の固有値(フレーム法+ベアストウ法) */
    #      coded by Y.Suganuma                 */
    ############################################
    			# データの設定
    ct  = 1000
    eps = 1.0e-10
    p0  = 0.0
    q0  = 0.0
    n   = 3
    a   = np.empty(n+1, np.float)
    b   = np.empty(n+1, np.float)
    c   = np.empty(n+1, np.float)
    r   = np.empty(n, np.complex)
    A   = np.array([[7, 2, 1],[2, 1, -4],[1, -4, 2]], np.float)
    H1  = np.empty((n, n), np.float)
    H2  = np.empty((n, n), np.float)
    			# 計算
    ind = Frame(n, ct, eps, p0, q0, a, b, c, r, A, H1, H2)
    			# 出力
    if ind > 0 :
    	print("収束しませんでした!")
    else :
    	for i1 in range(0, n) :
    		print("   " + str(r[i1]))
    			

  7. C#

    /********************************************/
    /* 行列の固有値(フレーム法+ベアストウ法) */
    /*      coded by Y.Suganuma                 */
    /********************************************/
    using System;
    
    class Program
    {
    	static void Main()
    	{
    		Test1 ts = new Test1();
    	}
    }
    
    class Test1
    {
    	public Test1()
    	{
    					// データの設定
    		int ct       = 1000;
    		int n        = 3;
    		double eps   = 1.0e-10;
    		double p0    = 0.0;
    		double q0    = 0.0;
    		double[] a   = new double [n+1];
    		double[] b   = new double [n+1];
    		double[] c   = new double [n+1];
    		double[] rl  = new double [n];
    		double[] im  = new double [n];
    		double[][] A = new double [][] {
    			new double[] {7.0, 2.0, 1.0},
    			new double[] {2.0, 1.0, -4.0},
    			new double[] {1.0, -4.0, 2.0}
    		};
    		double[][] H1 = new double [n][];
    		double[][] H2 = new double [n][];
    		for (int i1 = 0; i1 < n; i1++) {
    			H1[i1] = new double [n];
    			H2[i1] = new double [n];
    		}
    					// 計算
    		int ind = Frame(n, ct, eps, p0, q0, a, b, c, rl, im, A, H1, H2);
    					// 出力
    		if (ind > 0)
    			Console.WriteLine("収束しませんでした!");
    		else {
    			for (int i1 = 0; i1 < n; i1++)
    				Console.WriteLine("   " + 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                      */
    	/*************************************************/
    	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)
    	{
    		a[0] = 1.0;
    					// a1の計算
    		a[1] = 0.0;
    		for (int i1 = 0; i1 < n; i1++)
    			a[1] -= A[i1][i1];
    					// a2の計算
    		for (int i1 = 0; i1 < n; i1++) {
    			for (int i2 = 0; i2 < n; i2++)
    				H1[i1][i2] = A[i1][i2];
    			H1[i1][i1] += a[1];
    		}
    
    		a[2] = 0.0;
    		for (int i1 = 0; i1 < n; i1++) {
    			for (int i2 = 0; i2 < n; i2++)
    				a[2] -= A[i1][i2] * H1[i2][i1];
    		}
    		a[2] *= 0.5;
    					// a3からanの計算
    		for (int i1 = 3; i1 <= n; i1++) {
    
    			for (int i2 = 0; i2 < n; i2++) {
    				for (int i3 = 0; i3 < n; i3++) {
    					H2[i2][i3] = 0.0;
    					for (int 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 (int i2 = 0; i2 < n; i2++) {
    				for (int i3 = 0; i3 < n; i3++)
    					a[i1] -= A[i2][i3] * H2[i3][i2];
    			}
    			a[i1] /= i1;
    
    			for (int i2 = 0; i2 < n; i2++) {
    				for (int i3 = 0; i3 < n; i3++)
    					H1[i2][i3] = H2[i2][i3];
    			}
    		}
    					// ベアストウ法
    		int 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)
    	{
    		int ind = 0;
    	/*
    	          1次の場合
    	*/
    		if (n == 1) {
    			if (Math.Abs(a[0]) < eps)
    				ind = 1;
    			else {
    				rl[p] = -a[1] / a[0];
    				im[p] = 0.0;
    			}
    		}
    	/*
    	          2次の場合
    	*/
    		else if (n == 2) {
    			double D;
    						// 1次式
    			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;
    				}
    			}
    						// 2次式
    			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;
    				}
    			}
    		}
    						// 3次以上の場合
    		else {
    							// 因数分解
    			ind = 1;
    			int count = 0;
    			double D, dp, dq, p1 = p0, p2 = 0.0, q1 = q0, q2 = 0.0;
    			while (ind > 0 && count <= ct) {
    				for (int 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 (int 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) {
    							// 2次方程式を解く
    				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 (int 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;
    	}
    }
    			

  8. VB

    ''''''''''''''''''''''''''''''''''''''''''''
    ' 行列の固有値(フレーム法+ベアストウ法) '
    '      coded by Y.Suganuma                 '
    ''''''''''''''''''''''''''''''''''''''''''''
    Module Test
    	Sub Main()
    					' データの設定
    		Dim ct As Integer = 1000
    		Dim n As Integer  = 3
    		Dim eps As Double = 1.0e-10
    		Dim p0 As Double  = 0.0
    		Dim q0 As Double  = 0.0
    		Dim a(n+1) As Double
    		Dim b(n+1) As Double
    		Dim c(n+1) As Double
    		Dim rl(n) As Double
    		Dim im(n) As Double
    		Dim AA(,) As Double = {{7.0, 2.0, 1.0}, {2.0, 1.0, -4.0}, {1.0, -4.0, 2.0}}
    		Dim H1(n,n) As Double
    		Dim H2(n,n) As Double
    					' 計算
    		Dim ind As Integer = Frame(n, ct, eps, p0, q0, a, b, c, rl, im, AA, H1, H2)
    					' 出力
    		If ind > 0
    			Console.WriteLine("収束しませんでした!")
    		Else
    			For i1 As Integer = 0 To n-1
    				Console.WriteLine("   " & rl(i1) & "  i " & im(i1))
    			Next
    		End If
    
    	End Sub
    
    	'''''''''''''''''''''''''''''''''''''''''''''''''
    	' 行列の固有値(フレーム法+ベアストウ法)      '
    	'      n : 次数                                 '
    	'      ct : 最大繰り返し回数                    '
    	'      eps : 収束判定条件                       '
    	'      p0, q0 : x2+px+qにおけるp,qの初期値      '
    	'      a : 係数(最高次から与え,値は変化する) '
    	'      b, c : 作業域((n+1)次の配列)           '
    	'      rl, im : 結果の実部と虚部                '
    	'      AA : 行列                                '
    	'      H1, H2 : 作業域(nxnの行列)             '
    	'      return : =0 : 正常                       '
    	'               =1 : 収束せず                   '
    	'      coded by Y.Suganuma                      '
    	'''''''''''''''''''''''''''''''''''''''''''''''''
    	Function Frame(n As Integer, ct As Integer, eps As Double, p0 As Double,
    	               q0 As Double, a() As Double, b() As Double, c() As Double,
    	               rl() As Double, im() As Double, AA(,) As Double, H1(,) As Double,
    	               H2(,) As Double)
    
    		a(0) = 1.0
    					' a1の計算
    		a(1) = 0.0
    		For i1 As Integer = 0 To n-1
    			a(1) -= AA(i1,i1)
    		Next
    					' a2の計算
    		For i1 As Integer = 0 To n-1
    			For i2 As Integer = 0 To n-1
    				H1(i1,i2) = AA(i1,i2)
    			Next
    			H1(i1,i1) += a(1)
    		Next
    
    		a(2) = 0.0
    		For i1 As Integer = 0 To n-1
    			For i2 As Integer = 0 To n-1
    				a(2) -= AA(i1,i2) * H1(i2,i1)
    			Next
    		Next
    		a(2) *= 0.5
    					' a3からanの計算
    		For i1 As Integer = 3 To n
    
    			For i2 As Integer = 0 To n-1
    				For i3 As Integer = 0 To n-1
    					H2(i2,i3) = 0.0
    					For i4 As Integer = 0 To n-1
    						H2(i2,i3) += AA(i2,i4) * H1(i4,i3)
    					Next
    				Next
    				H2(i2,i2) += a(i1-1)
    			Next
    
    			a(i1) = 0.0
    			For i2 As Integer = 0 To n-1
    				For i3 As Integer = 0 To n-1
    					a(i1) -= AA(i2,i3) * H2(i3,i2)
    				Next
    			Next
    			a(i1) /= i1
    
    			For i2 As Integer = 0 To n-1
    				For i3 As Integer = 0 To n-1
    					H1(i2,i3) = H2(i2,i3)
    				Next
    			Next
    		Next
    					' ベアストウ法
    		Dim ind As Integer = Bairstow(n, ct, eps, p0, q0, a, b, c, rl, im, 0)
    
    		Return ind
    
    	End Function
    
    	'''''''''''''''''''''''''''''''''''''''''''''''''
    	' 実係数代数方程式の解(ベアストウ法)          '
    	'      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                      '
    	'''''''''''''''''''''''''''''''''''''''''''''''''
    	Function Bairstow(n As Integer, ct As Integer, eps As Double, p0 As Double,
    	                  q0 As Double, a() As Double, b() As Double, c() As Double,
    	                  rl() As Double, im() As Double, p As Integer)
    
    		Dim ind As Integer = 0
    				'
    				' 1次の場合
    				'
    		If n = 1
    			If Math.Abs(a(0)) < eps
    				ind = 1
    			Else
    				rl(p) = -a(1) / a(0)
    				im(p) = 0.0
    			End If
    				'
    				' 2次の場合
    				'
    		ElseIf n = 2
    			Dim D As Double
    						' 1次式
    			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
    				End If
    						' 2次式
    			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
    				End If
    			End If
    						' 3次以上の場合
    		Else
    							' 因数分解
    			ind = 1
    			Dim count As Integer = 0
    			Dim D As Double
    			Dim dp As Double
    			Dim dq As Double
    			Dim p1 As Double = p0
    			Dim p2 As Double = 0.0
    			Dim q1 As Double = q0
    			Dim q2 As Double = 0.0
    
    			Do While ind > 0 and count <= ct
    				For i1 As Integer = 0 To n
    					If i1 = 0
    						b(i1) = a(i1)
    					ElseIf i1 = 1
    						b(i1) = a(i1) - p1 * b(i1-1)
    					Else
    						b(i1) = a(i1) - p1 * b(i1-1) - q1 * b(i1-2)
    					End If
    				Next
    				For i1 As Integer = 0 To n
    					If i1 = 0
    						c(i1) = b(i1)
    					ElseIf i1 = 1
    						c(i1) = b(i1) - p1 * c(i1-1)
    					Else
    						c(i1) = b(i1) - p1 * c(i1-1) - q1 * c(i1-2)
    					End If
    				Next
    				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 and Math.Abs(dq) < eps
    						ind = 0
    					Else
    						count += 1
    						p1 = p2
    						q1 = q2
    					End If
    				End If
    			Loop
    	
    			If ind = 0
    							' 2次方程式を解く
    				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
    				End If
    							' 残りの方程式を解く
    				n -= 2
    				For i1 As Integer = 0 To n
    					a(i1) = b(i1)
    				Next
    				ind = Bairstow(n, ct, eps, p0, q0, a, b, c, rl, im, p+2)
    			End If
    		End If
    	
    		Return ind
    
    	End Function
    
    End Module
    			

情報学部 菅沼ホーム 目次 索引