Project Euler - Problem 247
「続きを読む」の使い方を覚えたのでトップページがシャキっとした．
本来こんなことしてる暇はないが追い詰められるほど感情が無になりやるべきことをやりたくなくなる病気にかかっているため Problem 247 をやります．
Squares under a hyperbola
これはインデックス (3, 3) までの図です
rectangles = {};
genRect[left_, below_, x_, y_] :=
Module[
{edge = N[(-(x + y) + Sqrt[(x - y)^2 + 4])/2]},
If[Or[left > 3, below > 3], Return[]];
AppendTo[
rectangles,
{EdgeForm[Thin], Transparent, Rectangle[{x, y}, {x + edge, y + edge}]}
];
genRect[left + 1, below, x + edge, y];
genRect[left, below + 1, x, y + edge];
];
genRect[0, 0, 1, 0];
Show[
Plot[
1/x, {x, 1, 3}, PlotRange -> {0, 1},
PlotStyle -> {Thin, Black},
AspectRatio -> 0.5, ImageSize -> 700
],
Graphics[rectangles]
]
これは画像の生成呪文です
はてなのシンタックスハイライトめちゃくちゃ多言語対応してるっぽいけど Mathematica 先生はなかった(たぶん)
フラクタル構造のようなものがみえているので再帰で解けそう
正方形ですが，
- 左下の座標が (x, y)
- 右上が双曲線に接している
- 一辺 a
とすると
$\displaystyle a = \frac{ -(x + y) + \sqrt{ (x - y)^{2} + 4 } }{ 2 }$
なのでこれを使う．
インデックス (l, b)，左下の座標 (x, y) の正方形を (l, b, x, y) で表すと，
- 上に (l, b + 1, x, y + a)
- 右に (l + 1, b, x + a, y)
の正方形が生えます．
これで再帰を書ける．方針としては
- インデックス (3, 3) の正方形の中で最も番号の大きい( = 辺の長さが短い)正方形の辺の長さ d を求める
- 辺の長さ d 以下の正方形を全てカウントする
以上 2 ステップに分けてやった．
d を求める
double d = 1; void rec(int l, int b, double x, double y) { if (l > 3 || b > 3) return; double a = (-(x + y) + sqrt((x - y) * (x - y) + 4)) / 2; if (a < d) d = a; rec(l + 1, b, x + a, y); rec(l, b + 1, x, y + a); } int main() { rec(0, 0, 1, 0); return 0; }
これはそのままですね
辺の長さ d 以下の正方形を全てカウントする
double d; int counter = 0; void rec2(double x, double y) { double a = (-(x + y) + sqrt((x - y) * (x - y) + 4)) / 2; if (a < d) return; counter++; rec2(x, y + a); rec2(x + a, y); } int main() { rec2(1, 0); cout << counter << endl; return 0; }
これもそのままですね，と言いたいが再帰回数がべらぼうに多くスタックオーバーフローするので再帰を使わずに書かなければいけない．
typedef pair<double, double> Pdd; double d; int counter = 0; void rec2_stack() { double x, y, a; stack<Pdd> st; st.push(Pdd(1, 0)); while (!st.empty()) { x = st.top().first; y = st.top().second; a = (-(x + y) + sqrt((x - y) * (x - y) + 4)) / 2; st.pop(); if (a < d) continue; counter++; st.push(Pdd(x, y + a)); st.push(Pdd(x + a, y)); } } int main() { rec2_stack(); cout << counter << endl; return 0; }
スタックでなんとかしたやつ．これで 10 秒くらい．
ただよく観察すると b = 0 の正方形，つまり x 軸に接している正方形の数がバカみたいに多く，b > 1 の正方形は l < 300 ぐらいの領域にしか存在しないことがわかるので，b = 0 とそれ以外で分けてカウントするとだいぶ速くなる．スタックも使わなくていい．
double d; int counter = 0; // b = 0 void below0() { double a = 1; double x = 0; while (1) { x += a; a = (-x + sqrt(x * x + 4)) / 2; if (a < d) break; counter++; } } // それ以外 void rec3(int l, int b, double x, double y) { double a = (-(x + y) + sqrt((x - y) * (x - y) + 4)) / 2; if (a < d || l > 300) return; if (b != 0) counter++; rec3(l + 1, b, x + a, y); rec3(l, b + 1, x, y + a); } int main() { below0(); rec3(0, 0, 1, 0); cout << counter << endl; return 0; }
これで 0.02 秒くらいになった．めでたしめでたし
コード全体
#include <iostream> #include <time.h> #include <stack> using namespace std; typedef long long int lint; typedef pair<double, double> Pdd; double d = 1; int counter = 0; // d を求める void rec(int l, int b, double x, double y) { if (l > 3 || b > 3) return; double a = (-(x + y) + sqrt((x - y) * (x - y) + 4)) / 2; if (a < d) d = a; rec(l + 1, b, x + a, y); rec(l, b + 1, x, y + a); } // stack overflow void rec2(double x, double y) { double a = (-(x + y) + sqrt((x - y) * (x - y) + 4)) / 2; if (a < d) return; counter++; rec2(x, y + a); rec2(x + a, y); } // stack version void rec2_stack() { double x, y, a; stack<Pdd> st; st.push(Pdd(1, 0)); while (!st.empty()) { x = st.top().first; y = st.top().second; a = (-(x + y) + sqrt((x - y) * (x - y) + 4)) / 2; st.pop(); if (a < d) continue; counter++; st.push(Pdd(x, y + a)); st.push(Pdd(x + a, y)); } } // speedy void below0() { double a = 1; double x = 0; while (1) { x += a; a = (-x + sqrt(x * x + 4)) / 2; if (a < d) break; counter++; } } void rec3(int l, int b, double x, double y) { double a = (-(x + y) + sqrt((x - y) * (x - y) + 4)) / 2; if (a < d || l > 300) return; if (b != 0) counter++; rec3(l + 1, b, x + a, y); rec3(l, b + 1, x, y + a); } void PE0247() { rec(0, 0, 1, 0); //rec2_stack(); rec3(0, 0, 1, 0); below0(); cout << counter << endl; } int main() { clock_t start, end; start = clock(); PE0247(); end = clock(); printf("%d msec.\n", end - start); return 0; }
