原始ピタゴラス数生成いろいろ

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

package {
    import flash.display.Sprite;
    import flash.text.TextField;
    import flash.utils.getTimer;
    public class Test extends Sprite {
        private var _tf : TextField;
  
        public function Test() {
            _tf = new TextField();
            _tf.width = 465;
            _tf.height = 465;
            addChild(_tf);
            
            var s : int;
            
            tr("斜辺が3000000までの原始ピタゴラス数の組の個数を数えてみたよ！");
            tr("ちなみに、これ以上に増やすとStern-Brocotの再帰版は死ぬよ！");
            tr();
            s = getTimer(); tr(testNaive(3000000)); tr((getTimer() - s) + "ms");tr();
            s = getTimer(); tr(testSternBrocot(3000000)); tr((getTimer() - s) + "ms");tr();
            s = getTimer(); tr(testSternBrocot2(3000000)); tr((getTimer() - s) + "ms");tr();
            s = getTimer(); tr(testMatrix(3000000)); tr((getTimer() - s) + "ms");tr();
        }
  
  		public function testNaive(N : int) : Number
  		{
  			tr("naive");
  			tr("互いに素なm>nでm,nの偶奇が一致しないときに、");
  			tr("(m*m-n*n, 2*m*n, m*m+n*n)は原始ピタゴラス数の組になる。");
  			var sqn : uint = Math.sqrt(N);
  			var ct : uint = 0;
  			for(var m : uint = 1;m <= sqn;m++){
  				for(var n : uint = 1 + (m & 1);n < m && m * m + n * n <= N;n += 2){
  					if(GCD(m, n) == 1){
//  						tr(m * m - n * n, 2 * m * n, m * m + n * n);
  						ct++;
  					}
  				}
  			}
  			return ct;
  		} 
  
		public function testSternBrocot(N : int) : Number
  		{
  			tr("Stern-Brocot木");
  			tr("有理数を列挙する木。これを使って、互いに素の判定を省略できる。");
  			return countWithFraction(0, 1, 1, 1, N);
  		}
  		
		public function testSternBrocot2(N : int) : Number
  		{
  			tr("Stern-Brocot木2");
  			tr("上の再帰をなくし、固定長Vectorにスタックしたバージョン。");
  			return countPythagoras(N);
        }
        
        private function testMatrix(N : uint) : Number
        {
        		tr("Matrix");
        		tr("ピタゴラス数を生み出す行列で生成。一番無駄がないはず。");
        		return countWithMatrix(3, 4, 5, N);
        }
        
        private static function GCD(a : uint, b : uint) : uint
        {
            while(b > 0){
                var c : uint = a;
                a = b;
                b = c % b;
            }
            return a;
        }
        
        private function countPythagoras(N : int) : Number
        {
			var sn : uint = 0;
			var sd : uint = 1;
			var en : uint = 1;
			var ed : uint = 1;
			var seg : Vector.<uint> = new Vector.<uint>(int(Math.sqrt(N) * 2));
        		var ret : Number = 0;
        		var p : uint = 0;
        		while(true){
        			var nn : uint = sn + en;
        			var nd : uint = sd + ed;
            		if(nn * nn + nd * nd < N){
	            		ret += (nn ^ nd) & 1;
					seg[p] = en;
					seg[p + 1] = ed;
					p += 2;
					en = nn; ed = nd;
        			}else{
        				if(p == 0)break;
        				sn = en;
        				sd = ed;
        				en = seg[p - 2];
        				ed = seg[p - 1];
        				p -= 2;
        			} 
        		}
        		return ret;
        }

        private function countWithFraction(sn : int, sd : int, en : int, ed : int, N : int) : Number
        {
        		var nn : int = sn + en;
        		var nd : int = sd + ed;
        		if(nn * nn + nd * nd > N)return 0;
        		return ((nn ^ nd) & 1) + countWithFraction(sn, sd, nn, nd, N) + countWithFraction(nn, nd, en, ed, N);
        }
        
        private function countWithMatrix(a : uint, b : uint, c : uint, N : uint) : Number
        {
        		if(c > N)return 0;
        		return 1 + countWithMatrix(a-2*b+2*c, 2*a-b+2*c, 2*a-2*b+3*c, N) +
        			countWithMatrix(-a+2*b+2*c, -2*a+b+2*c, -2*a+2*b+3*c, N) +
        			countWithMatrix(a+2*b+2*c, 2*a+b+2*c, 2*a+2*b+3*c, N);
        }
        
        private function tr(...o : Array) : void
        {
            _tf.appendText(o + "\n");
        }
    }
}