Project Euler 227

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

package {
    import flash.display.Sprite;
    import flash.text.TextField;
    import flash.utils.getTimer;
    // @see http://projecteuler.net/index.php?section=problems&id=227
    public class Euler227 extends Sprite {
        private var _tf : TextField;
  
        public function Euler227() {
            _tf = new TextField();
            _tf.width = 465;
            _tf.height = 465;
            addChild(_tf);
            
            var s : int = getTimer();
            tr(solve(50)); 
				
            var g : int = getTimer();
            tr((g - s) + " ms");
        }

        // 距離の離散マルコフ過程で、距離0の点からの遷移確率を0にしたものを考える。
        // 遷移行列は、距離が変わらない確率18/36, ±1だけ増える確率8/36, ±2だけ増える確率1/36
        // を基本とした51*51の五重対角行列Aになる。
        // 期待値は、初期状態をvとして、
        // Ex=0v+1Av+2A^2v+・・
        // と表される。これを計算すると、Ex=(E-A)^-2 Av
        // となるのでこれを計算すればよい。(E-A)^-2を左辺に移行すれば連立方程式になるので、
        // それを解くことにする。
        // この問題の場合のvは51次ベクトルで、最後の項だけ1で他は0.
        private function solve(N : int) : Number
        {
        		var M : int = N + 1;
        		var A : Array = new Array(M * M);
        		var i : uint;
        		for(i = 0;i < A.length;i++)A[i] = 0;
        		/*
        		for(i = 2;i <= N-2;i++){
        			A[(i-2)*M+i] = 1 / 36;
        			A[(i-1)*M+i] = 8 / 36;
        			A[(i  )*M+i] = 18 / 36;
        			A[(i+1)*M+i] = 8 / 36;
        			A[(i+2)*M+i] = 1 / 36;
        		}
    			A[(M-2-2)*M+M-2] = 1 / 36;
    			A[(M-2-1)*M+M-2] = 8 / 36;
    			A[(M-2  )*M+M-2] = 19 / 36;
    			A[(M-2+1)*M+M-2] = 8 / 36;
    			
    			A[(M-1-2)*M+M-1] = 2 / 36;
    			A[(M-1-1)*M+M-1] = 16 / 36;
    			A[(M-1  )*M+M-1] = 18 / 36;
    			
    			A[(1-1)*M+1] = 8 / 36;
    			A[(1  )*M+1] = 19 / 36;
    			A[(1+1)*M+1] = 8 / 36;
    			A[(1+2)*M+1] = 1 / 36;
    			*/
        		for(i = 2;i <= N-2;i++){
        			A[(i-2)*M+i] = 1;
        			A[(i-1)*M+i] = 8;
        			A[(i  )*M+i] = 18;
        			A[(i+1)*M+i] = 8;
        			A[(i+2)*M+i] = 1;
        		}
    			A[(M-2-2)*M+M-2] = 1;
    			A[(M-2-1)*M+M-2] = 8;
    			A[(M-2  )*M+M-2] = 19;
    			A[(M-2+1)*M+M-2] = 8;
    			
    			A[(M-1-2)*M+M-1] = 2;
    			A[(M-1-1)*M+M-1] = 16;
    			A[(M-1  )*M+M-1] = 18;
    			
    			A[(1-1)*M+1] = 8;
    			A[(1  )*M+1] = 19;
    			A[(1+1)*M+1] = 8;
    			A[(1+2)*M+1] = 1;
    			
    			A[(0  )*M+0] = 0;
    			
    			// A-E
    			var AE : Array = A.concat();
    			for(i = 0;i < M;i++){
    				AE[i*M+i] -= 36;
    			}
    			
    			// (A-E)^2
    			var AE2 : Array = mulMatMat(AE, AE, M, M, M);
    			/*
    			for(i = 0;i < M;i++){
    				tr(AE2.slice(i*M, (i+1)*M));
    			}
    			*/ 
    			
    			var p : Array = new Array(M);
    			for(i = 0;i < M;i++)p[i] = 0;
    			p[M-1] = 1;
    			p = mulMatVec(A, p, M, M);
    			
    			var ex : Array = gaussEliminate(AE2, p);
    			
    			// 有効数字10桁ってかいてあるのに11桁目が安定しないので正解だすのに苦労・・
    			return ex[0] * 36;
        }
        
        private static function mulMatVec(a : Array, b : Array, r : uint, c : uint) : Array
        {
        		var ret : Array = new Array(r);
        		for(var i : uint = 0;i < r;i++){
        			var sum : Number = 0;
        			for(var j : uint = 0;j < c;j++){
        				sum += a[i*c+j]*b[j];
        			}
        			ret[i] = sum;
        		}
        		return ret;
        }

        private static function mulMatMat(a : Array, b : Array, ar : uint, ac : uint, bc : uint) : Array
        {
        		var ret : Array = new Array(ar * bc);
        		for(var i : uint = 0;i < ar;i++){
        			for(var j : uint = 0;j < bc;j++){
	        			var sum : Number = 0;
	        			for(var k : uint = 0;k < ac;k++){
	        				sum += a[i*ac+k] * b[k*bc+j];
	        			}
	        			ret[i*bc+j] = sum;
        			}
        		}
       		return ret;
        }

        private function gaussEliminate(ab : Array, bb : Array) : Array
        {
        		var a : Array = ab.concat();
        		var b : Array = bb.concat();
            var n : int = b.length;
            if(n == 0 || a.length != n * n)return null;
            
            var i : int, j : int, k : int;
            for(i = 0;i < n;i++){
                var pmax : Number = -1;
                var p : int = -1;
                for(j = i;j < n;j++){
                    var aji : Number = a[j * n + i];
                    if(Math.abs(aji) > pmax){
                        pmax = Math.abs(aji);
                        p = j;
                    }
                }
                if(p == -1)return null;
                
                if(p != i){
                    var d : Number;
                    for(j = 0;j < n;j++){
                        d = a[i * n + j];
                        a[i * n + j] = a[p * n + j];
                        a[p * n + j] = d;
                    }
                    d = b[i];
                    b[i] = b[p];
                    b[p] = d;
                }
                
                for(j = i + 1;j < n;j++){
                    var m : Number = a[j * n + i] / a[i * n + i];
                    a[j * n + i] = 0;
                    for(k = i + 1;k < n;k++){
                        a[j * n + k] -= a[i * n + k] * m;
                    }
                    b[j] -= b[i] * m;
                }
            }
                
            if(a[n * n - 1] == 0.0)return null;
            var x : Array = new Array(n);
            x[n - 1] = b[n - 1] / a[n * n - 1];
            for(i = n - 2;i >= 0;i--){
                x[i] = b[i];
                for(j = i+1;j < n;j++){
                    x[i] -= a[i * n + j] * x[j];
                }
                x[i] /= a[i * n + i];
            }
            tr(x); 
            return x;
        }
        
        private function tr(...o : Array) : void
        {
            _tf.appendText(o + "\n");
            _tf.scrollV = _tf.maxScrollV;
        }
    }
}