/**
* Copyright uwi ( http://wonderfl.net/user/uwi )
* MIT License ( http://www.opensource.org/licenses/mit-license.php )
* Downloaded from: http://wonderfl.net/c/nbKK
*/
package {
import flash.display.Sprite;
import flash.text.TextField;
import flash.utils.getTimer;
// @see http://projecteuler.net/index.php?section=problems&id=
public class Euler extends Sprite {
private var _tf : TextField;
// Euler Problem 258の予定
public function Euler() {
_tf = new TextField();
_tf.width = 465;
_tf.height = 465;
addChild(_tf);
var s : int = getTimer();
tr(solve());
// tr(exGCD(2008,2002));
var g : int = getTimer();
tr((g - s) + " ms");
}
private function solve() : int
{
var N : Number = Math.pow(10, 18);
var M0 : Number = 500000000000000; // 係数パスカル三角形の左端の長さ
var M1 : Number = Math.floor(N / 1999); // 係数パスカル三角形の右端の長さ
var n2000 : Number = M1 / 2000; // 2000列の個数
var mod : int = (n2000 - Math.floor(n2000)) * 2000; // テーブルの周回数の余りのユニットの余りの長さ
n2000 = Math.floor(n2000);
tr(M1, n2000);
var ret : int = 0;
var mul : int = 1;
var p : int;
var sum : int;
var i : int, j : int;
for each(p in [2, 5]){
var modp : int = modd(n2000, p);
var t : Array = cmodtable(p);
sum = 0;
sum += t[(modp - 1) * p + 0];
for(j = 1;j < mod;j++){
sum += t[modp * p + (j % p)];
}
sum %= p;
tr(p, modp, sum);
ret = chinese(ret, sum, mul, p);
mul *= p;
}
tr(ret);
// 20092010 = 2 * 5 * 859 * 2339
for each(p in [859, 2339]){
var rep : Number = n2000 / p; // テーブルの周回数
var rep2 : int = modd(n2000, p);
var t : Array = cmodtable(p);
// ユニットをつくる
// 1周をつくる
var cir : int = 0;
var d2 : int = 2000 / p;
var m2 : int = 2000 % p;
for(i = 0;i < p;i++){
var mi : int = i * p;
var start : int = (i * 2000) % p;
cir += exmod(2, i, p) * d2;
for(j = 0;j < m2;j++){
cir += t[mi + ((start + j) % p)];
}
cir %= p;
}
var sum : int = (cir * modd(Math.floor(rep), p)) % p;
tr(p, Math.floor(rep), rep2, mod, cir, sum);
d2 = 2000 / p;
m2 = 2000 % p;
for(i = 0;i < rep2;i++){
mi = i * p;
start = (i * 2000) % p;
sum += exmod(2, i, p) * d2;
for(j = 0;j < m2;j++){
sum += t[mi + ((start + j) % p)];
}
sum %= p;
}
// sum = sum + t[0] - t[(i - 1) * p + (((i - 1) * 2000) % p)];
// ユニットの余りを作る
d2 = mod / p;
m2 = mod % p;
mi = (i % p) * p;
start = (i * 2000) % p;
sum += t[(mi - p) + start];
sum += exmod(2, i % p, p) * d2;
for(j = 1;j < m2;j++){
sum += t[mi + ((start + j) % p)];
}
sum %= p;
tr(p, sum);
ret = chinese(ret, sum, mul, p);
mul *= p;
}
return ret;
}
private function chinese(a1 : int, a2 : int, m1 : int, m2 : int) : Number
{
var e : Array = exGCD(m1, m2);
if(e[0] != 1)return 0;
for(var ret : Number = a1 + (a2 - a1) * e[1] * m1;ret < 0;ret += m1 * m2);
return ret;
}
// (a,b)=1のときの、ax+by=1の解
private function exGCD(a : int, b : int) : Array
{
return rec(a, b, 1, 0, 0, 1);
}
private function rec(a : int, b : int, p : int, q : int, r : int, s : int) : Array
{
if(b == 0)return [a, p, r];
var c : int = a / b;
return rec(b, a % b, q, p - c * q, s, r - c * s);
}
// a : prime
// C(n, k) % aのテーブル(n<kのときは0)
private function cmodtable(a : int) : Array
{
var invs : Array = new Array(a);
for(var i : int = 1;i < a;i++){
invs[i] = inv(i, a);
}
var ret : Array = new Array(a * a);
for(var r : int = 0;r < a;r++){
var v : int = 1;
for(var c : int = 0;c < a;c++){
ret[r * a + c] = c > r ? 0 : v;
v = (v * (r - c)) % a;
v = (v * invs[c + 1]) % a;
}
}
return ret;
}
// p : prime
private function inv(n : int, p : int) : int
{
// p-2乗する
return exmod(n, p - 2, p);
}
// 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 modd(a : Number, b : int) : int
{
var r : Number = a / b;
return int(Math.round((r - Math.floor(r)) * b));
}
private function tr(...o : Array) : void
{
_tf.appendText(o + "\n");
}
}
}
