Bairstow's method

1変数高次方程式の解を求める
精度はfindRootsの第2引数で調整
@see http://www.sist.ac.jp/~suganuma/kougi/other_lecture/SE/num/num.htm#3.1

by uwi 09 May 2010

Embed

/**
 * Copyright uwi ( http://wonderfl.net/user/uwi )
 * MIT License ( http://www.opensource.org/licenses/mit-license.php )
 * Downloaded from: http://wonderfl.net/c/4xGn
 */

package {
    import flash.display.Sprite;
    import flash.events.*;
    import com.bit101.components.*;
    
    // 1変数高次方程式の解を求める
    // 精度はfindRootsの第2引数で調整
    // @see http://www.sist.ac.jp/~suganuma/kougi/other_lecture/SE/num/num.htm#3.1
    public class RootFinder extends Sprite {
    		private var _text : Text;
		private var _label : Text;
		private var _rootsRe : Text;
		private var _rootsAll : Text;
		private var _solve : PushButton;
    		    	
        public function RootFinder() {
            _text = new Text(this, 30, 60, "1 2 3 4");
            _text.width = 400;
            _text.height = 20;
            _text.textField.restrict = "-0-9 +\.";
            _text.textField.addEventListener(Event.CHANGE, onTextChange);
            
            _label = new Text(this, 30, 100, "");
            _label.editable = false;
            _label.width = 400;
            _label.height = 40;
            
            var lr : Label = new Label(this, 30, 160, "real roots");
            _rootsRe = new Text(this, 90, 160, "");
            _rootsRe.editable = false;
            _rootsRe.width = 340;
            _rootsRe.height = 70;
            
            var la : Label = new Label(this, 30, 230, "all roots");
            _rootsAll = new Text(this, 90, 230, "");
            _rootsAll.editable = false;
            _rootsAll.width = 340;
            _rootsAll.height = 70;
            
            onTextChange(null);
            
            _solve = new PushButton(this, 130, 300, "Solve it!", onSolve);
            _solve.width = 200;
            _solve.height = 100;
        }
        
        private function onTextChange(e : Event = null) : void
        {
        		var a : Array = _text.text.split(" ");
        		var text : String = ""; 
        		var first : Boolean = true;
        		for(var i : uint = 0;i < a.length;i++){
        			var deg : int = a.length - 1 - i;
        			a[i] = parseFloat(a[i]);
        			if(isNaN(a[i]) || a[i] == 0)continue;
        			if(first){
        				first = false;
        			}else{
        				if(a[i] >= 0)text += "+";
        			}
        			if(deg == 0){
        				text += a[i];
        			}else if(deg == 1){
        				text += (a[i] == 1 ? "" : (a[i] == -1 ? "-" : a[i])) + "x";
        			}else{
        				text += (a[i] == 1 ? "" : (a[i] == -1 ? "-" : a[i])) + "x^" + deg;
        			}
        		}
        		_label.text = text + " = 0";
        }
        
        private function onSolve(e : MouseEvent) : void
        {
        		var c : Array = []
        		for each(var el : String in _text.text.split(" ")){
        			var num : Number = parseFloat(el);
        			c.push(isNaN(num) ? 0 : num);
        		}
        		
        		// real roots
        		var retre : Array = findRoots(c, 0.0001, true);
        		retre.sort(Array.NUMERIC);
        		_rootsRe.text = "x = " + (retre.length ? retre.join(",") : "None!");
        		
        		// all roots
        		var retall : Array = findRoots(c, 0.0001, false);
        		retall.sortOn(["re", "im"], Array.NUMERIC);
        		var strs : Array = [];
        		for each(var rt : Object in retall){
        			strs.push(
        				(rt.re||"") + (rt.im ? ((rt.im>=0 && rt.re!=0?"+":"") + rt.im + "i") : "")
        				);
        		}
        		_rootsAll.text = "x = " + (retall.length ? strs.join(",") : "None!");
        }
        
        private static function findRoots(a : Array, eps : Number, realOnly : Boolean) : Array
        {
        		cutHead(a);
        		var p : Number = a.length >= 1 ? a[1]/a[0] : 0;
        		var q : Number = a.length >= 2 ? a[2]/a[0] : 0;
        		var ret : Array = [];
        		while(true){
        			var n : uint = a.length;
        			if(n == 3){
        				var sqa : Number = a[1] * a[1] - 4 * a[0] * a[2];
        				if(realOnly){
	        				if(sqa >= 0){
	        					ret.push((-a[1]+Math.sqrt(sqa))/(2*a[0]), (-a[1]-Math.sqrt(sqa))/(2*a[0]));
	        				}
        				}else{
        					if(sqa >= 0){
	        					ret.push({re:(-a[1]+Math.sqrt(sqa))/(2*a[0])}, {re:(-a[1]-Math.sqrt(sqa))/(2*a[0])});
	        				}else{
	        					ret.push(
	        					{re:-a[1]/(2*a[0]), im:Math.sqrt(-sqa)/(2*a[0])},
	        					{re:-a[1]/(2*a[0]), im:-Math.sqrt(-sqa)/(2*a[0])}
        						);
        					}
        				}
        				break;
    				}else if(n == 2){
    					if(realOnly){
	        				ret.push(-a[1]/a[0]);
    					}else{
    						ret.push({re:-a[1]/a[0]});
    					}
        				break;
        			}else if(n <= 1){
        				break;
        			}
        				
        			var b : Array = [];
        			var c : Array = [];
        			for(var i : uint = 0;i < n;i++){
        				b.push(a[i] - p*(b[i-1]||0) - q*(b[i-2]||0));
        				c.push(b[i] - p*(c[i-1]||0) - q*(c[i-2]||0));
        			}
    				var D : Number = c[n-3]*c[n-3]-c[n-4]*(c[n-2]-b[n-2]);
    				var dp : Number = (b[n-2]*c[n-3]-b[n-1]*c[n-4])/D;
    				var dq : Number = (b[n-1]*c[n-3]-b[n-2]*(c[n-2]-b[n-2]))/D;
    				
        			if(Math.abs(dp) < eps && Math.abs(dq) < eps){
        				// x^2+px+q=0
        				var sq : Number = p * p - 4 * q;
        				if(realOnly){
	        				if(sq >= 0){
	        					ret.push((-p+Math.sqrt(sq))/2, (-p-Math.sqrt(sq))/2);
	        				}
        				}else{
        					if(sq >= 0){
	        					ret.push({re:(-p+Math.sqrt(sq))/2}, {re:(-p-Math.sqrt(sq))/2});
	        				}else{
	        					ret.push(
	        					{re:-p/2, im:Math.sqrt(-sq)/2},
	        					{re:-p/2, im:-Math.sqrt(-sq)/2}
        						);
        					}
        				}
        				a = b.slice(0, n - 2);
        			}else{
        				p += dp;
        				q += dq;
        			}
        		} 
    			return ret;
        }
        
        private static function cutHead(a : Array) : Array
        {
        		while(a.length > 0 && a[0] == 0)a.shift();
        		return a;
        }
    }
}