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

<?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;
/*
          １次の場合
*/
	if ($n == 1) {
		if (abs($a[0]) < $eps)
			$ind = 1;
		else {
			$rl[$k] = -$a[1] / $a[0];
			$im[$k] = 0.0;
		}
	}
/*
          ２次の場合
*/
	else if ($n == 2) {
					// １次式
		if (abs($a[0]) < $eps) {
			if (abs($a[1]) < $eps)
				$ind = 1;
			else {
				$rl[$k] = -$a[2] / $a[1];
				$im[$k] = 0.0;
			}
		}
					// ２次式
		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;
			}
		}
	}
					// ３次以上の場合
	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) {
						// ２次方程式を解く
			$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;
}

?>