てきとーな日記

2010-01-02

[]Σ電卓(2) 18:38

いろいろ改良したらかなりまともになったのでとりあえず完成

六角形の個数を求めるやつだと、

$x1$y1$x2$y2$x3$y3$x4$y4$x5$y5$x6$y6[0<=x1,y1,x2,y2,x3,y3,x4,y4,x5,y5,x6,y6<=N][x1+y1,x2+y2,x3+y3,x4+y4,x5+y5,x6+y6<=N][x1<x2][y1=y2][y2<y3][x2=x3][x3-y3>x4-y4][x3+y3=x4+y4][x4>x5][y4=y5][y5>y6][x5=x6][x6-y6<x1-y1][x6+y6=x1+y1]@N>=1

こんな入力を与えると、

long res = N * (12 + N * (4 + N * (-15 + N * (-5 + N * (3 + N))))) / 720;
if ((1 + N) % 2 == 0) res += (45 + N * (-48 + N * (-52 + N * (60 + N * (5 + N * (-12 + 2 * N)))))) / 2880;
if (N % 2 == 0) res += N * (-48 + N * (-52 + N * (60 + N * (5 + N * (-12 + 2 * N))))) / 2880;
return res;

こんな出力を吐いてくれるようになった。

オプションでmodとったり、BigIntegerにしたりも可能。

SRM454の方は、20行くらいあって長過ぎorz

もっと頑張れそうな気もするけどめんどいからいいやw

そのうちどっかで公開するかも。

ひっそりと公開しました http://ow.ly/SuLM

パスワードはてきと

2009-12-23

[]Σ電卓 00:33

目的

TopCoderとかでいちいちΣ計算するのがだるかったのでプログラムで自動化を試みた。

ただΣ計算するだけではあまり意味が無いので、条件式を自動で場合分けして計算するものを考えたところ、式は多項式に限定し、条件式は線形不等式と定数を法とする合同式に限定すれば、サイズ小さければ解けそう…ということになったので作成開始。

最近のでいうとSRM455の550やSRM454の1000などを解けるようなのを目指す。

方針

[条件式]で条件が真の時1、偽の時0を返す関数を表すことにし、

¥sum_{x=a}^b f(x)

$x[a<=x][x<=b]f(x)

と書く。($xはxが整数全体を動くことを表す)

条件式は線形不等式[A<=0]と、線形合同式[a|A](Aがaで割り切れるか)に限定(aが小さいことを前提としているので合同式は線形でなくても良い気もする)

あとは、てきとーに構文木を作った上で、再帰的に

modの除去

$x [c | ax + B] f(x)

= [(a, c) | B] $x [c / (a, c) | a / (a,c) x + B / (a, c)] f(x)

= [(a, c) | B] $x f( (cx - mB) / (a, c) ) #mはa/(a,c)のc/(a,c)での逆元

場合分け

[ax <= A][bx <= B]

= [bA <= aB][ax <= A] + [bA > aB][bx <= B]

Σ計算

$x[A <= ax][bx <= B]p(x)

= [bA <= aB]($j[0 <= j < b][b | B - j]P( (B - j) / b ) - $i[0 <= i < a][a | A + i]P( (A + i) / a - 1) ) #P(x)はΣ_{i=1}^x p(i)

を行っていけば最終的に定数ループだけになって解けるはず

ただし、条件式が多かったり係数が大きいと、場合分けと定数ループで項数が指数爆発して大変なことに…

実装

最初めんどくさいの構文解析くらいだろと思ってOCamlでやったが、途中で合同式が必須なことに気づき、大変なことに…

Javaにしようかとも思ったけど結局OCamlでそのまま書いた

それっぽいのが出来上がってSRM455の550をクリア

これが計算に使ったもとの式

$n[3<=n][n<=N](N-n+1)(N-n+2)$a$b$c[1<=a][1<=b][1<=c][a+b<n][b+c<n][c+a<n]

計算結果は出力手抜きしたままなので項の数が70近くあって大変なことに…

long f(long N) {
	return (3 * cond(-1 * N + 7 * 1 <= 0) * cond((1 * N + 1 * 1) % 2 == 0) * pow(N, 6)
			 + -9 * cond(-1 * N + 7 * 1 <= 0) * cond((1 * N + 1 * 1) % 2 == 0) * pow(N, 5)
			 + -15 * cond(-1 * N + 7 * 1 <= 0) * cond((1 * N + 1 * 1) % 2 == 0) * pow(N, 4)
			 + 90 * cond(-1 * N + 7 * 1 <= 0) * cond((1 * N + 1 * 1) % 2 == 0) * pow(N, 3)
			 + -123 * cond(-1 * N + 7 * 1 <= 0) * cond((1 * N + 1 * 1) % 2 == 0) * pow(N, 2)
			 ...
			 + 5760 * cond(-1 * N + 3 * 1 <= 0)) % M * pow(2880 * 2, M - 2);

てきとーに同じ条件分岐でくくったりすればもう少しはましになりそう

さらに、この式思いつかなくてもだれでも解けるように頂点6個の座標*2の変数12個にしてゴリゴリ六角形条件を書き連ねた場合

$x1$y1$x2$y2$x3$y3$x4$y4$x5$y5$x6$y6

[0<=x1][0<=y1][x1+y1<=n]

[0<=x2][0<=y2][x2+y2<=n]

[0<=x3][0<=y3][x3+y3<=n]

[0<=x4][0<=y4][x4+y4<=n]

[0<=x5][0<=y5][x5+y5<=n]

[0<=x6][0<=y6][x6+y6<=n]

[x1<x2][y1=y2]

[x2-y2=x3-y3][x2+y2<x3+y3]

[x3=x4][y3<y4]

[x4>x5][y4=y5]

[x5-y5=x6-y6][x5+y5>x6+y6]

[x6=x1][y6>y1]

こんな感じ?(未検証)

も試してみたところ項数が指数爆発して停止しないorz

手抜きして等号は不等式2つに置換してるので、とりあえず手動で等号除去して変数6つにしたら解けた!

ただし項数が300近くて、しかもオーバーフローしたのであってるのかは知らない

改良

とりあえず等号の処理を特別扱いして、多倍長演算にし、出力まともにすればコンテストで出る問題くらいは解けるようになるのかな…

というかこれくらいすでにどっかにあってもいい気がするのだが…(Mathematica先生は本当にやってくれないのだろうか…)