Project Euler 160

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

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

        // N!が素因数2を十分な数含んでいるとする。
        
        // 求める値は、N!の0を除いた下位5桁なので、
        // N!から10をできるだけ抜く、つまり、
        // N!のなかの素因数5の数をMとしてN!/10^Mを100000でわった余りを求めればよい。
        // まずP=N!/5^Mについて求める。Pは5と互いに素なので、
        // (A)100000の中の素因数5だけで構成された3125(=5^5)でPをわった余りを求めて、
        // (B)それを残りの32(=2^5)で割れば(32の、3125についての剰余群の逆元をかければ)よい。
        // N!/10^Mは、上の値からさらに2^Mで割った値なので、(B)の^5を^(M+5)にすればよい。
        
        // (A)は、N!の5n番目の項抜きのを3125でわった余りを求めて、
        // それに [N/5]!の5n番目の項抜きの・・と繰り返して求めてすべてかければよい。
        // (B)は、オイラーの定理より、3125と互いに素なaについて、
        // a^φ(3125)=1 (mod 3125)なので、それを使えばよい。
        
        // 最後に2^5をかけて答えを得る。
        private function solve(N : Number) : int
        {
            var u : Number;
            
            // 素因数5の個数を求める
            var n5 : Number = 0;
            for(u = N;u > 0;u = Math.floor(u / 5)){
                n5 += Math.floor(u / 5);
            }
            tr(n5);
            
            // N!の5n番目の項抜きの3125を法とした余りのテーブルをつくる
            var modf : Array = new Array(3125);
            var mul : int = 1;
            for(var i : int = 0;i < 3125;i++){
                if(i % 5 != 0) mul = (mul * i) % 3125;
                modf[i] = mul;
            }
            
            var ret : int = 1;
            // (A)
            for(u = N;u > 0;u = Math.floor(u / 5)){
                var q : int = u / 3125;
                var r : int = u - q * 3125;
                // overflowするかもしれない
                ret = (ret * exmod(modf[3124], q, 3125) * modf[r]) % 3125;
            }
            
            // (B)
            var t3125 : int = 3125 * 4 / 5;
            var n5mod : int = (n5 + 5) - Math.floor((n5 + 5) / t3125) * t3125;
            ret = (ret * exmod(2, t3125 - n5mod, 3125)) % 3125;
            
            return ret * 32;
        }
        
        // a^n(mod p) ただし a<2^16
        private function exmod(a : int, n : int, p : int) : int
        {
            var mul : int = a;
            var ret : int = 1;
            for(var i : int = n;i > 0;i >>= 1){
                if((i & 1) == 1){
                    ret *= mul;
                    ret %= p;
                }
                mul *= mul;
                mul %= p;
            }
            return ret;
        }

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