静岡理工科大学 菅沼ホーム 目次 索引

代数方程式(ベアストウ法)

    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;
    /*
              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

      アプレット版では,任意の代数方程式に対して,画面上で解を得ることができます.
    /************************************/
    /* 代数方程式の解(ベアストウ法)   */
    /*      例:(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;
    	/*
    	          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;
    	}
    }
    			
    アプレット版
    /****************************/
    /* 代数方程式(ベアストウ) */
    /*      coded by Y.Suganuma */
    /****************************/
    import java.awt.*;
    import java.awt.event.*;
    import java.applet.*;
    import java.util.StringTokenizer;
    
    public class Bairstow extends Applet implements ActionListener {
    	TextField order, coef, _p0, _q0, max;
    	TextArea ta;
    	Button bt;
    	Font f = new Font("TimesRoman", Font.BOLD, 20);
    
    	/************/
    	/* 初期設定 */
    	/************/
    	public void init()
    	{
    					// レイアウト,背景色,フォント
    		setLayout(new BorderLayout(5, 5));
    		setBackground(new Color(225, 255, 225));
    		setFont(f);
    					// 上のパネル
    		Panel pn1 = new Panel();
    		add(pn1, BorderLayout.NORTH);
    
    		pn1.add(new Label("次数:"));
    		order = new TextField("5", 3);
    		pn1.add(order);
    
    		pn1.add(new Label(" p0:"));
    		_p0 = new TextField("0", 3);
    		pn1.add(_p0);
    
    		pn1.add(new Label(" q0:"));
    		_q0 = new TextField("0", 3);
    		pn1.add(_q0);
    
    		pn1.add(new Label(" 最大繰り返し回数:"));
    		max = new TextField("1000", 5);
    		pn1.add(max);
    
    		pn1.add(new Label("  "));
    		bt = new Button("OK");
    		bt.setBackground(Color.pink);
    		bt.addActionListener(this);
    		pn1.add(bt);
    					// 中央のパネル(問題)
    		Panel pn2 = new Panel();
    		add(pn2, BorderLayout.CENTER);
    
    		pn2.add(new Label("係数(次数が高い順):"));
    		coef = new TextField("1 -3 -2 3 7 6", 50);
    		pn2.add(coef);
    					// 下のパネル(答え)
    		Panel pn3 = new Panel();
    		add(pn3, BorderLayout.SOUTH);
    
    		ta = new TextArea(10, 30);
    		pn3.add(ta);
    	}
    
    	/******************************/
    	/* ボタンが押されたときの処理 */
    	/******************************/
    	public void actionPerformed(ActionEvent e)
    	{
    		if (e.getSource() == bt) {
    			double a[], b[], c[], rl[], im[], p0, q0, eps;
    			int ind, ct, n;
    					// 問題の設定
    			eps  = 1.0e-10;
    			n    = Integer.parseInt(order.getText());
    			ct   = Integer.parseInt(max.getText());
    			p0   = Double.parseDouble(_p0.getText());
    			q0   = Double.parseDouble(_q0.getText());
    			a    = new double [n+1];
    			b    = new double [n+1];
    			c    = new double [n+1];
    			rl   = new double [n];
    			im   = new double [n];
    
    			String s = coef.getText();
    			StringTokenizer str = new StringTokenizer(s, " \n");
    			int k = 0;
    			while (str.hasMoreTokens()) {
    				a[k] = Double.parseDouble(str.nextToken());
    				k++;
    			}
    					// 計算
    			ind = bairstow(n, ct, eps, p0, q0, a, b, c, rl, im, 0);
    							// 結果
    			if (ind > 0) {
    				ta.setForeground(Color.red);
    				ta.setText(" 解を求めることができません");
    			}
    			else {
    				ta.setForeground(Color.black);
    				ta.setText("");
    				for (int i1 = 0; i1 < n; i1++)
    					ta.append(" " + String.format("%.5f",rl[i1]) + " i " + String.format("%.5f",im[i1]) + "\n");
    			}
    		}
    	}
    
    	/*************************************************/
    	/* 実係数代数方程式の解(ベアストウ法)          */
    	/*      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                      */
    	/*************************************************/
    	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 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;
    		/*
    		          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="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;
    /*
              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

    #***********************************/
    # 代数方程式の解(ベアストウ法)   */
    #      例:(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
    #
    #          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    = 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
    			# 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
    
    ############################################
    # 代数方程式の解(ベアストウ法)
    #      例:(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;
    	/*
    	          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

    ''''''''''''''''''''''''''''''''''''
    ' 代数方程式の解(ベアストウ法)   '
    '      例:(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
    				'
    				' 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
    			

静岡理工科大学 菅沼ホーム 目次 索引