/**
* Copyright uwi ( http://wonderfl.net/user/uwi )
* MIT License ( http://www.opensource.org/licenses/mit-license.php )
* Downloaded from: http://wonderfl.net/c/eb8R
*/
package {
import flash.display.Sprite;
import flash.text.TextField;
import flash.utils.getTimer;
// @see http://projecteuler.net/index.php?section=problems&id=152
public class Euler152 extends Sprite {
private var _tf : TextField;
private var T2 : Array;
// 1/2を1/2^2,..,1/80^2から選んだ和であらわす組み合わせは何通りあるか。
public function Euler152() {
_tf = new TextField();
_tf.width = 465;
_tf.height = 465;
addChild(_tf);
var s : int = getTimer();
tr(solve(45));
// tr(solve(80));
var g : int = getTimer();
tr((g - s) + " ms");
}
// 分母が互いに素である分数同士を足すと、
// 和の分数には必ず互いの素因子が残ってしまうことを利用する。
//
// 上記の性質により、2~80のうち素因子が1個しか入っていない
// 40以上の素数の項は除外。
// 残りのうち、素因数の大きい方から見ていって、
// 分母からその素因数を消せる和の組み合わせをみつけていく。
// たとえば、1/1^2+1/2^2+1/3^3=49/36だが、
// これは素因数7を消せる和
// 1/7^2+1/14^2+1/21^2=1/36となる。
//
// このようにして2以外の素因数を分母から除いていく。
// 各項の最大素因数を見ているので、重複して除かれることはない。
// 最後に分母が2の累乗だけで構成された和を計算して1/2と
// 等しいかどうかチェックする。
private function solve(N : int) : Number
{
var primes : Array = sieveEratosthenes(N);
var i : uint, j : uint, k : uint, l : uint;
// 1/2^nの項
T2 = [];
for(i = 2;i < N;i *= 2){
T2.push([1, i*i]);
}
// 素因数を消せる和を構成できる分母(?)を並べる
var validss : Array = [[]];
for(i = 1;i < primes.length;i++){
var p : uint = primes[i];
if(p > N / 2)break;
var M : uint = N / p;
var validList : Array = new Array(M);
for(j = 1;j <= M;j++)validList[j] = 1;
for(k = i+1;k < primes.length && primes[k] <= M;k++){
for(l = primes[k];l <= M;l+=primes[k]){
validList[l] = 0;
}
}
var valids : Array = [];
for(j = 1;j <= M;j++){
if(validList[j] == 1){
valids.push(j);
}
}
tr(p, valids);
validss.push(valids);
}
// 素因数の大きい方から消せる和を構成していく。
return rec(validss, primes, validss.length - 1, [0, 1]);
}
private function rec(validss : Array, primes : Array, pos : int, sum : Array) : Number
{
var ret : Number = 0;
var ve : Array;
if(pos == 0){
// 素因数2のとき。和が1/2になるかチェック
for(var c : uint = 0;c < 1 << T2.length;c++){
var base : Array = sum;
for(var i : uint = 0;i < T2.length;i++){
if(((c >> i) & 1) != 0){
base = addFraction(base, T2[i]);
}
}
if(base[1] / base[0] == 2)ret++;
// if(base[0] == 1 && base[1] == 2)ret++;
}
return ret;
}
// 消せる和を列挙
var p : uint = primes[pos];
var ves : Array = enumValidExpressions(validss[pos], p, 0, [sum[0] * p * p, sum[1]]);
// tr(p, sum, ves, pos);
for each(ve in ves){
ret += rec(validss, primes, pos - 1, ve);
}
return ret;
}
// sumに、validsにある分母を用いて和をつくり、分母からpを消せる和を列挙
private function enumValidExpressions(valids : Array, p : uint, pos : uint, f : Array) : Array
{
if(f[0] / f[1] > p * p / 2)return []; // 和が1/2をこえたら意味がない
if(pos == valids.length){
// tr(f, ":", con);
return (f[0] % (p * p) == 0 && f[0] / f[1] <= p * p / 2) ? [[f[0] / p / p, f[1]]] : [];
}
var ve : Array = enumValidExpressions(valids, p, pos + 1, f);
return ve.concat(enumValidExpressions(valids, p, pos + 1, addFraction(f, [1, valids[pos] * valids[pos]]) ));
}
/*
private function enumValidExpressionsO(valids : Array, p : uint, pos : uint, f : Array, con : Array) : Array
{
if(pos == valids.length){
// tr(f, ":", con);
return (f[0] % (p * p) == 0 && f[0] / f[1] <= p * p / 2) ? [con.concat([f[0], f[1]])] : [];
}
var ve : Array = enumValidExpressionsO(valids, p, pos + 1, f, con);
return ve.concat(enumValidExpressionsO(valids, p, pos + 1, addFraction(f, [1, valids[pos] * valids[pos]]), con.concat(valids[pos])) );
}
*/
// 分数の足し算。 0:分子, 1:分母
private function addFraction(a : Array, b : Array) : Array
{
var d : int = a[1] * b[1];
var n : int = a[0] * b[1] + b[0] * a[1];
if(n == 0)return [0, 1];
var gcd : int = GCD(n, d);
return [n / gcd, d / gcd];
}
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 sieveEratosthenes(n : int) : Array
{
var nn : uint = ((n / 2 - 1) >>> 5) + 1;
var ar : Vector.<uint> = new Vector.<uint>(nn);
var i : uint, j : uint;
for(i = 0;i < nn - 1;i++)ar[i] = 0xffffffff;
ar[nn - 1] = (1 << ((n / 2 - 1) & 31)) - 1;
var sq : uint = (Math.sqrt(n) - 3) >>> 1;
for(var p : uint = 0;p <= sq;p++){
if(ar[p >>> 5] & (1 << (p & 31))){
var m : uint = (p << 1) + 3;
var m2 : uint = m << 1;
for(var mm : uint = m * m;mm <= n;mm += m2){
var ind : uint = (mm - 3) >>> 1;
ar[ind >>> 5] &= ~(1 << (ind & 31));
}
}
}
var ret : Array = [2];
for(i = 0;i < nn;i++){
for(j = 0;j <= 31;j++){
if(ar[i] & (1 << j))ret.push((((i << 5) | j) << 1) + 3);
}
}
return ret;
}
private function tr(...o : Array) : void
{
_tf.appendText(o + "\n");
_tf.scrollV = _tf.maxScrollV;
}
}
}
