ルジャンドル記号の計算

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

package {
    import flash.display.Sprite;
    import flash.text.TextField;
    import flash.utils.getTimer;
    
    // @see http://ja.wikipedia.org/wiki/%E5%B9%B3%E6%96%B9%E5%89%B0%E4%BD%99%E3%81%AE%E7%9B%B8%E4%BA%92%E6%B3%95%E5%89%87
    // @see http://en.wikipedia.org/wiki/Quadratic_reciprocity
    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 = getTimer();
            
            for(var i : uint = 1;i < 23;i++){
                tr("(" + i + "|23)=" + calcLegendre(i, 23));
            }
            
            var g : int = getTimer();
            tr((g - s) + " ms");
        }
        
        private function calcLegendre(a : uint, b : uint) : int
        {
            if(GCD(a, b) != 1)return 0;
            
            // 分母を分解
            var ret : uint = 1;
            for(var i : uint = 2;i * i <= b;i++){
                for(var ct : uint = 0;b % i == 0;ct++){
                    b /= i;
                }
                if(ct % 2 == 1){
                    ret *= calcLegendreCore(a, i);
                }
            }
            if(b != 1)ret *= calcLegendreCore(a, b);
            return ret;
        }
        
        // (a,b)=1, b:odd prime
        private function calcLegendreCore(a : uint, b : uint) : int
        {
//            tr(a, b);
            if(a == 1){
                return 1;
            }
            
            // Jacobi symbol
            if(a == 2){
                return (b % 8 == 1 || b % 8 == 7) ? 1 : -1;
            }
            if(a == b - 1){
                return b % 4 == 1 ? 1 : -1;
            }
            
            var ret : int = 1;
            if(a < b){
                // Quadratic Reciprocity
                if(a % 2 == 1 && b % 2 == 1 && isPrime(a) && isPrime(b)){
                    var c : uint = a; a = b; b = c;
                    if(a % 4 == 3 && b % 4 == 3){
                        ret = -ret;
                    }
                }
            }
            a %= b;
            
            for(var i : uint = 2;i * i <= a;i++){
                for(var ct : uint = 0;a % i == 0;ct++){
                    a /= i;
                }
                if(ct % 2 == 1){
                    ret *= calcLegendreCore(i, b);
                }
            }
            if(a != 1)ret *= calcLegendreCore(a, b);
            return ret;
        }
        
        private var _cache : Object = {};
        
        private function isPrime(n : Number) : Boolean
        {
            if(n <= 1)return false;
            if(n % 2 == 0)return n == 2;
            if(_cache[n])return _cache[n];
            
            var sq : int = Math.sqrt(n);
            for(var i : int = 3;i <= sq;i+=2){
                if(n % i == 0){
                    _cache[n] = false;
                    return false;
                }
            }
            _cache[n] = true;
            return true;
        }
        
	private static function GCD(a : uint, b : uint) : uint
	{
	    return b == 0 ? a : GCD(b, a % b);
	}

        private function tr(...o : Array) : void
        {
            _tf.appendText(o + "\n");
        }
    }
}