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

代数方程式（ベアストウ法）

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

　　プログラムは，実係数代数方程式 (x + 1)(x - 2)(x - 3)(x2 + x + 1) = 0 の解を，ベアストウ法で解いた例です．

  1. C++ 

    /************************************/
/* 代数方程式の解（ベアストウ法）   */
/*      例：(x+1)(x-2)(x-3)(x2+x+1) */
/*           =x5-3x4-2x3+3x2+7x+6=0 */
/*      coded by Y.Suganuma         */
/************************************/
#include <stdio.h>

int Bairstow(int, int, double, double, double,
             double *, double *, double *, double *, double *);

int main()
{
	double *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    = 5;
	a    = new double [n+1];
	b    = new double [n+1];
	c    = new double [n+1];
	rl   = new double [n];
	im   = new double [n];

	a[0] = 1.0;
	a[1] = -3.0;
	a[2] = -2.0;
	a[3] = 3.0;
	a[4] = 7.0;
	a[5] = 6.0;
					// 計算
	ind = Bairstow(n, ct, eps, p0, q0, a, b, c, rl, im);
					// 出力
	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;

	return 0;
}

/*************************************************/
/* 実係数代数方程式の解（ベアストウ法）          */
/*      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;
/*
          １次の場合
*/
	if (n == 1) {
		if (fabs(a[0]) < eps)
			ind = 1;
		else {
			rl[0] = -a[1] / a[0];
			im[0] = 0.0;
		}
	}
/*
          ２次の場合
*/
	else if (n == 2) {
					// １次式
		if (fabs(a[0]) < eps) {
			if (fabs(a[1]) < eps)
				ind = 1;
			else {
				rl[0] = -a[2] / a[1];
				im[0] = 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[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;
			}
		}
	}
					// ３次以上の場合
	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) {
						// ２次方程式を解く
			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 

    /************************************/
/* 代数方程式の解（ベアストウ法）   */
/*      例：(x+1)(x-2)(x-3)(x2+x+1) */
/*           =x5-3x4-2x3+3x2+7x+6=0 */
/*      coded by Y.Suganuma         */
/************************************/
import java.io.*;

public class Test {
	public static void main(String args[]) throws IOException
	{
		double 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    = 5;
		a    = new double [n+1];
		b    = new double [n+1];
		c    = new double [n+1];
		rl   = new double [n];
		im   = new double [n];

		a[0] = 1.0;
		a[1] = -3.0;
		a[2] = -2.0;
		a[3] = 3.0;
		a[4] = 7.0;
		a[5] = 6.0;
					// 計算
		ind = Bairstow(n, ct, eps, p0, q0, a, b, c, rl, im, 0);
					// 出力
		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 : 結果の実部と虚部                */
	/*      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;
	}
}

  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 aa  = (document.getElementById("coe").value).split(/ {1,}/);
			for (let i1= 0; i1 < aa.length; i1++)
				a[i1]  = parseFloat(aa[i1]);
			let b   = new Array();
			let c   = new Array();
			let rl  = new Array();
			let im  = Array();

			let ind = bairstow(n, ct, eps, p0, q0, a, b, c, rl, im, 0);
					// 出力
			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 : 結果の実部と虚部                */
		/*      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;
		/*
		          １次の場合
		*/
			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;
		}
	</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="bairstow.gif"></P>
	</DL>

	<DIV STYLE="text-align:center">
		次数：<INPUT ID="order" STYLE="font-size: 100%" TYPE="text" SIZE="2" VALUE="5"> 
		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>
		係数（次数が高い順）：<INPUT ID="coe" STYLE="font-size: 100%;" TYPE="text" SIZE="50" VALUE="1 -3 -2 3 7 6"><BR><BR>
		<TEXTAREA ID="ans" COLS="50" ROWS="15" STYLE="font-size: 100%"></TEXTAREA>
	</DIV>

</BODY>

</HTML>

  4. PHP 

    <?php

/************************************/
/* 代数方程式の解（ベアストウ法）   */
/*      例：(x+1)(x-2)(x-3)(x2+x+1) */
/*           =x5-3x4-2x3+3x2+7x+6=0 */
/*      coded by Y.Suganuma         */
/************************************/
					// データの設定
	$ct  = 1000;
	$eps = 1.0e-10;
	$p0  = 0.0;
	$q0  = 0.0;
	$n   = 5;
	$a   = array($n+1);
	$b   = array($n+1);
	$c   = array($n+1);
	$rl  = array($n);
	$im  = array($n);

	$a[0] = 1.0;
	$a[1] = -3.0;
	$a[2] = -2.0;
	$a[3] = 3.0;
	$a[4] = 7.0;
	$a[5] = 6.0;
					// 計算
	$ind = Bairstow($n, $ct, $eps, $p0, $q0, $a, $b, $c, $rl, $im, 0);
					// 出力
	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 : 結果の実部と虚部                */
/*      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;
}

?>

  5. Ruby 

    #***********************************/
# 代数方程式の解（ベアストウ法）   */
#      例：(x+1)(x-2)(x-3)(x2+x+1) */
#           =x5-3x4-2x3+3x2+7x+6=0 */
#      coded by Y.Suganuma         */
#***********************************/

#************************************************/
# 実係数代数方程式の解（ベアストウ法）          */
#      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
#
#          １次の場合
#
	if n == 1
		if a[0].abs() < eps
			ind = 1
		else
			rl[k] = -a[1] / a[0]
			im[k] = 0.0
		end
#
#          ２次の場合
#
	elsif n == 2
					# １次式
		if a[0].abs() < eps
			if a[1].abs() < eps
				ind = 1
			else
				rl[k] = -a[2] / a[1]
				im[k] = 0.0
			end
					# ２次式
		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
					# ３次以上の場合
	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
						# ２次方程式を解く
			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    = 5
a    = [1.0, -3.0, -2.0, 3.0, 7.0, 6.0]
b    = Array.new(n+1)
c    = Array.new(n+1)
rl   = Array.new(n)
im   = Array.new(n)
				# 計算
ind = Bairstow(n, ct, eps, p0, q0, a, b, c, rl, im, 0)
				# 出力
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 : 結果
#      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
			# １次の場合
	if n == 1 :
		if abs(a[0]) < eps :
			ind = 1
		else :
			r[k] = complex(-a[1] / a[0], 0)
			# ２次の場合
	elif n == 2 :
				# １次式
		if abs(a[0]) < eps :
			if abs(a[1]) < eps :
				ind = 1
			else :
				r[k] = complex(-a[2] / a[1], 0)
				# ２次式
		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)
				# ３次以上の場合
	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 :
					# ２次方程式を解く
			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

############################################
# 代数方程式の解（ベアストウ法）
#      例：(x+1)(x-2)(x-3)(x2+x+1)
#           =x5-3x4-2x3+3x2+7x+6=0
#      coded by Y.Suganuma
############################################
			# データの設定
ct  = 1000
eps = 1.0e-10
p0  = 0.0
q0  = 0.0
n   = 5
a   = np.array([1.0, -3.0, -2.0, 3.0, 7.0, 6.0])
b   = np.empty(n+1, np.float)
c   = np.empty(n+1, np.float)
r   = np.empty(n, np.complex)
			# 計算
ind = Bairstow(n, ct, eps, p0, q0, a, b, c, r, 0)
			# 出力
if ind > 0 :
	print("収束しませんでした！")
else :
	for i1 in range(0, n) :
		print("   " + str(r[i1]))

  7. C# 

    /************************************/
/* 代数方程式の解（ベアストウ法）   */
/*      例：(x+1)(x-2)(x-3)(x2+x+1) */
/*           =x5-3x4-2x3+3x2+7x+6=0 */
/*      coded by Y.Suganuma         */
/************************************/
using System;

class Program
{
	static void Main()
	{
		Test1 ts = new Test1();
	}
}

class Test1
{
	public Test1()
	{
					// データの設定
		int ct      = 1000;
		int n       = 5;
		double eps  = 1.0e-10;
		double p0   = 0.0;
		double q0   = 0.0;
		double[] a  = {1.0, -3.0, -2.0, 3.0, 7.0, 6.0};
		double[] b  = new double [n+1];
		double[] c  = new double [n+1];
		double[] rl = new double [n];
		double[] im = new double [n];

					// 計算
		int ind = Bairstow(n, ct, eps, p0, q0, a, b, c, rl, im, 0);
					// 出力
		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 : 結果の実部と虚部                */
	/*      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;
	/*
	          １次の場合
	*/
		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) {
			double D;
						// １次式
			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;
			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) {
							// ２次方程式を解く
				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 

    ''''''''''''''''''''''''''''''''''''
' 代数方程式の解（ベアストウ法）   '
'      例：(x+1)(x-2)(x-3)(x2+x+1) '
'           =x5-3x4-2x3+3x2+7x+6=0 '
'      coded by Y.Suganuma         '
''''''''''''''''''''''''''''''''''''
Module Test
	Sub Main()
					' データの設定
		Dim ct As Integer = 1000
		Dim n As Integer  = 5
		Dim eps As Double = 1.0e-10
		Dim p0 As Double  = 0.0
		Dim q0 As Double  = 0.0
		Dim a() As Double = {1.0, -3.0, -2.0, 3.0, 7.0, 6.0}
		Dim b(n+1) As Double
		Dim c(n+1) As Double
		Dim rl(n) As Double
		Dim im(n) As Double

					' 計算
		Dim ind As Integer = Bairstow(n, ct, eps, p0, q0, a, b, c, rl, im, 0)
					' 出力
		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 : 結果の実部と虚部                '
	'      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
				'
				' １次の場合
				'
		If n = 1
			If Math.Abs(a(0)) < eps
				ind = 1
			Else
				rl(p) = -a(1) / a(0)
				im(p) = 0.0
			End If
				'
				' ２次の場合
				'
		ElseIf n = 2
			Dim D As Double
						' １次式
			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
						' ２次式
			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
						' ３次以上の場合
		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
							' ２次方程式を解く
				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

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