Project Euler 276

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

package {
    import flash.display.Sprite;
    import flash.text.TextField;
    import flash.utils.getTimer;
    // @see http://projecteuler.net/index.php?section=problems&id=276
    public class Euler276 extends Sprite {
        private var _tf : TextField;
  
        public function Euler276() {
            _tf = new TextField();
            _tf.width = 465;
            _tf.height = 465;
            addChild(_tf);
            
            var s : int = getTimer();
            tr(solve(1000000));
//            tr(solve(10000000));
//            test2();
            var g : int = getTimer();
            tr((g - s) + " ms");
        }
        
        // 解説は長いので下に
        private function solve(N : uint) : Array
        {
        		var V : Object = {};
        		V[0] = [0, 0];
        		V[1] = [0, 0];
        		V[2] = [0, 0];
        		
        		var p : uint = N / 3;
        		var sum : Array = [0, 0];
        		for(var i : uint = 3;i <= N;i++){
        			sum = N2.add(sum, [all(i), 0]);
        			if(uint(N / p) == i){
	        			V[i] = sum;
	        			while(uint(N / p) == i && p >= 1)p--;
        			}
        		}
        		
        		p = 0;
        		for(i = N / 3;i >= 1;i--){
        			var k : uint = N / i;
        			if(uint(N / i) > p){
	        			for(var j : uint = k / 3;j >= 2;j--){
	        				V[k] = N2.sub(V[k], V[uint(k / j)]);
//	        				V[k] -= V[uint(k / j)];
	        			}
	        			p = uint(N / i);
        			}
        		}
        		return V[N];
        }
        
        private function all(n : uint) : Number
        {
        		return S(n-1)-S(n-uint(n/3)-1)-uint(n/2)*uint(n/4)+uint(n/4)*uint(n/4)-uint(n/3)*(uint(n/3)-1)/2;
        }
        
        private function S(n : uint) : Number
        {
        		return n&1 ? (n*n-1)/4 : n*n/4;
        }
        
        // a<=b<=cのもとで三角不等式はa+b>cだけとなる。
        // P(N):a<=b<=cかつa+b>cかつa+b+c<N, 
        // Q(N):P(N)かつgcd(a,b,c)=1,
        // S(N)={(a,b,c)|Q(N)}, f(N)=#S(N)とおいておく。
        // 最終的にはf(10^7)を求めれば良い。
        //
        // さて、集合A(N)={(a,b,c)|P(N)}の要素で、
        // gcd(a,b,c)=kとなるものの個数はf([N/k])となる。
        // k=1,...,[N/3]で重複なくA(N)が覆えるため、
        // #A(N)=Σ_[k=1~[N/3]] f([N/k])
        // f(N)=#A(N) - Σ_[k=2~[N/3]] f([N/k])
        // となるので、Nについて漸化式を小さい方から求めていけば良い。
        // 
        // #A(N)の値を求めるために、
        // 集合t(n)={(a,b,c)|a<=b<=cかつa+b>cかつa+b+c=n}の要素を数えることを考える。
        // これは、a-b空間上で(0,n/2),(n/3,n/3),(n/4,n/4)を結ぶ三角形の周と内部の格子点の数に等しい。
        // ただし、a+b=N/2上の格子点はカウントしない。
        // s(n)=Σ_[k=1~n] [k/2]とおくと、
        // t(n) = s(n-1)-s(n-[n/3]-1)-[n/2]*[n/4]+[n/4]^2-[n/3]([n/3]-1)/2
        // と表される。
        // (ちなみにs(n)=(n^2-1)/4 (n:odd), n^2/4 (n:even))
        // これより、#A(N)=Σ_{k=1~N-1} t(k) として求められる。
        //
        // M以下のすべてのNについてf(N)を求めて行くと、O(M√M)かかる。
        // f(M)を求めるだけなら、{[M/k]|kは自然数}についてのみ求めていけば良い。
        // こうするとO(M)となる。
        
        private function test2() : void
        {
        		var a : uint, b : uint, c : uint;
        		var ct : uint = 0;
        		for(a = 1;a <= 34;a++){
        			for(b = a;b <= 68;b++){
        				var g : uint = GCD(a, b);
        				for(c = b;a + b + c <= 97;c++){
        					if(a + b > c && GCD(g, c) == 1){
//        						tr(a,b,c);
        						ct++;
        					}
        				}
        			}
        		}
        		tr("test2 : " + ct);
        }
        
        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 test() : void
        	{
        		var a : uint, b : uint, c : uint;
        		var ct : Array = new Array(11);
        		for(var i : uint = 0;i < 11;i++)ct[i] = 0;
        		
        		for(a = 1;a <= 10;a++){
        			for(b = a;b <= 10;b++){
        				for(c = b;c <= 10;c++){
        					if(a + b > c){
        						if(a + b + c <= 10){
        							for(i = a+b+c;i<=10;i++){
        								ct[i]++;
        							}
//        							ct[a+b+c]++;
        						}
        					}
        				}
        			}
        		}
        		tr(ct);
        	}
        
        private function tr(...o : Array) : void
        {
            _tf.appendText(o + "\n");
            _tf.scrollV = _tf.maxScrollV;
        }
        
    }
}

class N2
{
	public static const N : Number = 1E15;
	
	public static function add(a : Array, b : Array) : Array
	{
		var r0 : Number = a[0] + b[0];
		var rp : int = Math.floor(r0 / N);
		r0 -= rp * N;
		var r1 : Number = a[1] + b[1] + rp;
		return [r0, r1];
	}
	
	public static function sub(a : Array, b : Array) : Array
	{
		var r0 : Number = a[0] - b[0];
		var rp : int = Math.floor(r0 / N);
		r0 -= rp * N;
		var r1 : Number = a[1] - b[1] + rp;
		return [r0, r1];
	}
}