![プログラミング]()
![Polygon]()
プログラミング
4200000000Hzでゴリ押す！！！(2)
今回はアルゴリズムとCPUのパワーでパズルを全力で解きます。
Exponential Idle（Android版 / iOS版）というクリッカー系放置ゲーム中のミニゲームである、矢印パズルを攻略します。
こちらのnoteの記事にめっちゃ影響を受けて書くことにしました。
と言っても、初級と中級は私でも頑張れば解けたのと、上級も先行研究があるので、エキスパートに絞って話を進めていきます（まあ上級も同じプログラムをちょちょっといじれば簡単に解けますが）。訂正とお詫び前回のサムネイルの一部に「♪ラデツキー行進曲（1.5倍速）」とありましたが、正しくは1.543倍速です。謹んでお詫び申し上げます。ルール説明（エキスパート）図1：矢印パズル（エキスパート）のルール六角格子状に並べられたタイルをタップすると、そのタイルを含む周囲（最大）7個のタイルがそれぞれ時計回りに60°回転します。
それを繰り返してすべてのタイルの向きを上向きに揃える、というルールです。目標盤面を入力として受け取り、各タイルのタップ回数を出力する。基本の考え方1回タップすることに60°回転するので、6回同じ場所をタップすれば盤面の状態は変わらず、加法巡回群をなすような6元集合、$\mathbb{Z}/6\mathbb{Z}$ を持ってくる。解法案1参考文献[3]と同様に、$\mathbb{Z}/6\mathbb{Z}$ 上の37次の線形方程式を解くことに帰着させる方法です。が。
聡明な読者の方ならお気づきでしょう。$\mathbb{Z}/6\mathbb{Z}$ は体でしょうか。2,3,4は乗法逆元を持たないので体ではないですね。
すると線形方程式を解くときに除法が使えないということになり、掃き出し法は使えないことになります（ピボットを1にできないだけで掃き出せないわけではないだろうが、一般的なアルゴリズムをそのまま持ってくるのは厳しそう）。2掃き出しで解けない（難しい）と言っても、元は有限個しかないので $6^{37}=6.2 \times 10^{28}$ 通りすべてを試してしまえという考え方。猿。
いや無理やけど。もし1パターンを1nsで検証できたとしても6.2e19 s掛かります。実際は1パターンの検証で約7x37回の足し算をすることになるのでもっとかかります。
imos法で最適化すればどうにかなるようなレベルの話ではないです。やるだけ無駄。3大本命。DFSです。DFS自体は木の全探索なのですが、探索の必要がない部分を枝刈りすればだいぶ計算量が減ります。どれぐらい絞れるかは後述。
今回はこれを採用。プログラムこれが解を探索するプログラムです。
入力の数字は、ディスプレイモードを数字にして表示された数字ということになっています（内部的には-1して持っているが）。#include <iostream>
#include <vector>
#include <tuple>
#include <algorithm>
using index_t = std::tuple<int, int>;
using board_t = std::vector<std::vector<int>>;
const int mod = 6;
const size_t board_row = 7, board_col = 7;
void dfs(board_t &board, board_t &hands, index_t index);
void flip_around(board_t &board, index_t index, int amount);
index_t next(index_t index);
bool is_in_board(index_t index);
void print_board(const board_t board);
const index_t end_index = next({ 6, 6 });
int main() {
board_t board(board_row, std::vector<int>(board_col));
board_t hands(board_row, std::vector<int>(board_col));
for (index_t i = { 0, 0 }; i != end_index; i = next(i)) {
const auto [r, c] = i;
std::cout << r << " " << c << ": " << std::flush;
std::cin >> board[r][c];
board[r][c]--;
}
dfs(board, hands, { 0, 0 });
return 0;
}
void dfs(board_t &board, board_t &hands, index_t index) {
if (index == end_index) {
print_board(hands);
exit(0);
return;
}
auto [r, c] = index;
// 右端辺での枝刈り
if (c == 6 && r > 3 && board[r - 1][c - 1] != board[r - 1][c]) {
return;
}
// 下端辺での枝刈り
if (r == 6 && c > 3 && board[r - 1][c - 1] != board[r][c - 1]) {
return;
}
// 上端、左端辺以外での枝刈り
if (r != 0 && c != 0) {
// 左上マスを0にするように動かす
const int hand = (mod - board[r - 1][c - 1]) % mod;
hands[r][c] = hand;
flip_around(board, index, hand);
dfs(board, hands, next(index));
flip_around(board, index, mod - hand);
} else {
for (int hand = 0; hand < mod; hand++) {
hands[r][c] = hand;
flip_around(board, index, hand);
dfs(board, hands, next(index));
flip_around(board, index, mod - hand);
}
}
}
void flip_around(board_t &board, index_t index, int amount) {
// 正体不明だがちゃんと最初だけ初期化処理が走る
static std::array<index_t, 7u> dpos{
dpos[0] = { -1, -1 },
dpos[1] = { -1, 0 },
dpos[2] = { 0, -1 },
dpos[3] = { 0, 0 },
dpos[4] = { 0, 1 },
dpos[5] = { 1, 0 },
dpos[6] = { 1, 1 }
};
const auto [r, c] = index;
for (const auto [dr, dc] : dpos) {
if (!is_in_board({ r + dr, c + dc })) continue;
board[r + dr][c + dc] += amount;
board[r + dr][c + dc] %= mod;
}
}
index_t next(index_t index) {
auto [r, c] = index;
if (c == 6 || c - r >= 3) {
r++;
c = std::max(r - 3, 0);
} else {
c++;
}
return { r, c };
}
bool is_in_board(index_t index) {
auto [r, c] = index;
return std::abs(r - 3) <= 3 && std::abs(c - 3) <= 3 && std::abs(c - r) <= 3;
}
void print_board(const board_t board) {
for (size_t i = 0; i < board_row; i++) {
for (size_t j = 0; j < board_col; j++) {
if (is_in_board({ i, j })) {
std::cout << board[i][j] << " ";
} else {
std::cout << " ";
}
}
std::cout << "\n";
}
std::cout << std::endl;
}
こんなに const とか & とか書くんだったら最初からRustでやればよかったなと若干の後悔。
ぼちぼち解説していきます。
まずここ。using index_t = std::tuple<int, int>;
using board_t = std::vector<std::vector<int>>;index_t というエイリアスに std::tuple<int, int> を当てています。
index_t index = {row, column}; となっているならば、index は、row 行 column 列のパネルの位置を表します。
board_t には std::vector<std::vector<int>> を当てていますが、
board_t board(R, std::vector<int>(C)); としたならば、board は、R 行 C 列のパネルの二次元配列を表します。（入力と）パネルの現状態と出力する手数に使います。
図2：indexの指し方
次はここ。index_t next(index_t index) {
auto [r, c] = index;
if (c == 6 || c - r >= 3) {
r++;
c = std::max(r - 3, 0);
} else {
c++;
}
return { r, c };
}next関数はその名の通り、次に探索する（タップ回数を決定する）点を決定する関数で、DFSでの探索順に直接関わってきます。
定義はまあ見りゃわかりますね（ifの条件に c - r >= 3 とか書いてあるが、これは条件をいっぱい書きたくなかったからうまいこと纏めてあるだけ）。
2行目の auto [r, c] = index; ってのは構造化束縛というC++17以降で追加された機能です。分割代入だと思って大体問題ありません。
一応探索順を図に示しておくとこうです。図3：DFSの探索順この順序にしているのは、単純に分かりやすいという以外にもう一点理由があって、枝刈りがしやすいという点。
枝刈りの方法はdfs関数を見ればわかります。
void dfs(board_t &board, board_t &hands, index_t index) {
if (index == end_index) {
print_board(hands);
exit(0);
return;
}
auto [r, c] = index;
// 右端辺での枝刈り
if (c == 6 && r > 3 && board[r - 1][c - 1] != board[r - 1][c]) {
return;
}
// 下端辺での枝刈り
if (r == 6 && c > 3 && board[r - 1][c - 1] != board[r][c - 1]) {
return;
}
// 上端、左端辺以外での枝刈り
if (r != 0 && c != 0) {
// 左上マスを0にするように動かす
const int hand = (mod - board[r - 1][c - 1]) % mod;
hands[r][c] = hand;
flip_around(board, index, hand);
dfs(board, hands, next(index));
flip_around(board, index, mod - hand);
} else {
for (int hand = 0; hand < mod; hand++) {
hands[r][c] = hand;
flip_around(board, index, hand);
dfs(board, hands, next(index));
flip_around(board, index, mod - hand);
}
}
}私もそこまでアホではないし、この関数は深い再帰を伴うので引数の大部分は参照で受けるようにしました。
部のC#erの皆さんが卒倒するかもしれませんが、vector<T>は値型です。というか、C/C++に値型・参照型という分け方はありません。
明示しない限り常にすべて値渡しになるので、ちゃんと参照渡しになるよう指定してやります。
さておき、枝刈りの話。これも図に示したほうが早いでしょう。図4：dfs関数の枝刈り図3に示した探索順から分かるように、27,32,36のノードでは、そのノードを操作した後は左上と真上の変化がないので、その2箇所の状態は必ず同じでなければならず、そうでなければ後退して探索し直すことになります。34,35,36でも左上と左下で同じことが起こります。
また図3から分かるように、ある時点で探索しているパネルの左上にパネルがあったとき、探索しているパネルを動かした後に再び左上のパネルが変わることはないので、左上は0にならなくてはなりません。よって左上のパネルを元にただ一通りに決まります。逆に言えば、0～5の中から自由に選べるのは左上にパネルが存在しない0,1,2,3,4,9,15のたった7個のパネルのタップ回数のみです。したがって、検証の必要なパターン数は $6^7 \approx 10^{5.45}$ 通りだけになります（ただし、関数呼び出し自体が重いことと、これらの点の不正性が分かるまでに再帰の深いところに潜ることがあるのでやはりある程度計算に時間はかかる）。
こんな感じで殆どの点で枝刈りが行えることがわかりますね。
if (index == end_index) {
print_board(hands);
exit(0);
return;
}で、最後36番が終わればここの分岐にたどり着いて終わりです（なんでこの人exitとreturn両方書いてるんだろう）。
でも図4を見てもらえばわかりますが、36番ノードだけチェックが入ってないんですよね。
これでいいのかと言われれば良くて、解が存在するならば、36番ノードだけが非0であることはありません。
流石に解が存在しない問題は生成されないだろうということでここは一つ。
flip_around(board, index, hand);
dfs(board, hands, next(index));
flip_around(board, index, mod - hand);ところでこの部分ですが、再帰DFSにおける常套手段です。
ある変換（ここではある1パネルをn回タップすること）が可逆であるならば、dfs呼び出しから帰ってきた後に逆変換（ここでは同じパネルを6-n回タップすること）を行えばうまいことDFSになります。是非覚えておいてください。
以上。閑話void flip_around(board_t &board, index_t index, int amount) {
// 正体不明だがちゃんと最初だけ初期化処理が走る
static std::array<index_t, 7u> dpos{
dpos[0] = { -1, -1 },
dpos[1] = { -1, 0 },
dpos[2] = { 0, -1 },
dpos[3] = { 0, 0 },
dpos[4] = { 0, 1 },
dpos[5] = { 1, 0 },
dpos[6] = { 1, 1 }
};
...
}この部分、何ですか。自分で書いたけどどういう理屈で初期化されているのかわかりません。C++プロの方、教えて下さい。
aggregateならメンバ名を指定して初期化する方法があるけど、その延長と考えていいんですかね。参考文献[1]Conic Games「Exponential Idle - Google Play のアプリ」<https://play.google.com/store/apps/details?id=com.conicgames.exponentialidle>
[2]Gilles-Philippe Paille「「Exponential Idle」をApp Storeで」<https://apps.apple.com/jp/app/exponential-idle/id1538487382>
[3]ちゃそ「Exponential Idle #2 矢印パズル攻略法（と二元体GF(2)上の線形方程式について）」<https://note.com/so_ra_64/n/n9c2eb6a5ef6f>
[4]ますたー。/繰り上げP「焼き甜花ちゃんシリーズの作り方＆素材」<https://www.nicovideo.jp/watch/sm37861810>
[5]いもす研「いもす法」<https://imoz.jp/algorithms/imos_method.html>
![>]()
内容を見る
![プログラミング]()
![Polygon]()
プログラミング
4200000000Hzでゴリ押す！！！(1)
何でもかんでもアルゴリズムとマシンパワー（主にCPU）でゴリ押します。そういう企画です。
「AtCoder水下位競プロerである私ならこういう思考をするぞ」っていうのを書きます（個人の意見ですが）。
2021年 東京大学(理系) 数学 問4(4)初回は東大の入試問題をやります。この問題をパワーでゴリ押すのは \(O(1)\) 番煎じになりますが、気にしないことにします。
実はこれ春休みの帰省中に数人で集まって解いた問題の一つで、僕だけパワーで解決しようとしてました。
問題\({}_{2021}C_{37}\bmod 4\) を求めよ。
本来はいろいろと前置きがあって、それを利用して引数を落として計算するんですが、そんなの私には分かりません。
さて、この問題をどうやって解いてやりましょうか。
解法1\({}_{2021}C_{37}\) を求める
で、問題は \({}_{2021}C_{37}\) がでかすぎることです。どの程度か簡単に見積もってみましょう。
$$
\begin{aligned}
{}_{2021}C_{37}
&=\frac{2021 \times \dots \times 1985}{37 \times \dots \times 1} \\
&>\left(\frac{1985}{37}\right)^{37} \\
&>2^{212.58}
\end{aligned}
$$
困りました。これでは128ビット符号付き整数(→GCCの__int128、C++ Boostのint128_t、Rustのi128など)でも収まらず、直接計算するのは厳しいです。
十分な長さの固定長整数型や多倍長整数を利用すれば出来ないこともないので一応やってみましょう。
Pythonのint型は多倍長なので何食わぬ顔で書けば実装できます。
n = 2021
r = 37
MOD = 4
ans = 1
for i in range(min(r, n - r)):
ans = ans * (n - i) // (1 + i)
print(f'{n}C{r} = {ans}')
print(f'{n}C{r} % {MOD} = {ans % MOD}')
実行結果は次のようになりました。
2021C37 = 10549453167137272405699426118445044729600037976073052501105621045645733312800275
2021C37 % 4 = 3
な、成程……！
$${}_{2021}C_{37}=10549453167137272405699426118445044729600037976073052501105621045645733312800275$$
なんだ……！そして \(10549453167137272405699426118445044729600037976073052501105621045645733312800275\) を \(4\) で割ったあまりは \(3\) なんだ……！
多倍長整数での \(N \times M\) および \(N / M\) \((N\gg M)\) が \(O(\log N \log\log N)\) で計算できる（例えばFFTとNewton-Raphson法で実装）という仮定のもとで、
この方法による \({}_{n}C_{r} (r \leq n/2)\) の計算量は大きめに見積もると、
$$
\begin{aligned}&
O(\log {}_{n}C_{1} \log\log {}_{n}C_{1})
+ \dots + O(\log {}_{n}C_{r} \log\log {}_{n}C_{r}) \\
=& O(\log n \log\log n) + \dots + O(\log n^r \log\log n^r) \\
=& O(r^2 \log n \log\log n)
\end{aligned}
$$
となります（この評価は割とガバガバですが）。
解法2\({}_{n}C_{r}=\dfrac{{}_{n}P_{r}}{r!}\) とオイラーの定理を利用する
\({}_{2021}C_{37} \bmod 4\) を計算しようと思ったらこの方法が一番に思いつきました。
オイラーの定理は、競技プログラミング界隈では常識であるところのフェルマーの小定理の一般化です。
この定理は次のようなものです。
\(m\) 以下の正の整数のうち \(m\) に互いに素であるものの個数を \(\varphi(m)\) とする（言い換えれば、\(\varphi(m)=|\{n ; \gcd(m,n)=1, 1 \leq n \leq m\}|\)）。
\(a\) と \(m\) が互いに素であるとき、
$$a^{\varphi(m)} \equiv 1 \pmod{m}$$
これの何が嬉しいって、分子 \(a\) と分母 \(b\) を \(\bmod m\) で別々に計算できるので、オーバーフローさせずに固定長整数型で計算できます。
そして分子に \(\bmod m\) における分母の逆元 \(b^{-1}\) を掛ける事によって求めることが出来ます。できるはずでした。
しかしよーく考えてみましょう。\({}_{2021}P_{37},37! \pmod{4}\) っていくつになるでしょうか？
……そうです。\(0\) になります。\(0\) に何掛けたって \(4\) で割ったあまりが \(1\) になることはありません。
すなわち、\(b \equiv 0\) の逆元 \(b^{-1}\) なんてものは存在しないので、この方法は使えません。
\({}_{2021}C_{37} \bmod (10^9+7)\) であれば求められるんですけどね……
解法2'ただこれを回避する方法があります。
$${}_{n}C_{r}=\frac{a}{b},a={}_{n}P_{r},b=r!$$
とします。ここで、
$$
\begin{aligned}
a &=2^{\alpha_0} \times 3^{\alpha_1}
\times 5^{\alpha_2} \times 7^{\alpha_3} \times \dots \\
b &=2^{\beta_0} \times 3^{\beta_1}
\times 5^{\beta_2} \times 7^{\beta_3} \times \dots
\end{aligned}
$$
のように素因数分解し、約分することで
$$
{}_{n}C_{r}=\frac{a}{b}
=2^{\alpha_0-\beta_0} \times 3^{\alpha_1-\beta_1}
\times 5^{\alpha_2-\beta_2} \times 7^{\alpha_3-\beta_3} \times \dots
$$
として計算することができます（\({}_{n}C_{r}\) は整数であるから、任意の \(i\) について \(\alpha_i \geq \beta_i\) が成り立つ）。
実際にこれを実装してみました。
from collections import defaultdict
n = 2021
r = 37
MOD = 4
spf = list(range(n + 1))
for i in range(2, n + 1):
if spf[i] != i:
continue
for j in range(i * 2, n + 1, i):
if spf[j] == j:
spf[j] = i
pf = defaultdict(lambda:0)
for i in range(min(r, n - r)):
m = n - i
while m > 1:
pf[spf[m]] += 1
m //= spf[m]
m = 1 + i
while m > 1:
pf[spf[m]] -= 1
m //= spf[m]
ans = 1
for p, v in pf.items():
ans *= pow(p, v, MOD)
ans %= MOD
print(f'{n}C{r} % {MOD} = {ans}')
結果は当然解法1と同じく3です。
各自然数を \(O(\sqrt{n})\) で素因数分解して全体 \(O(r \sqrt{n})\) で求める方法もありますが、\(r \fallingdotseq n/2\) の場合も想定して、エラトステネスの篩と同様の前計算 \(O(n \log\log n)\) でSPF（Smallest Prime Factor）配列を求め、一つの自然数の素因数分解は \(O(\log n)\) で、全体 \(O(n \log\log n + r \log n)\) でできるようにしました。
全体の計算量は、\(O(n \log\log n + r \log n)\) です。
東大の問題の場合 \(r\) が小さいので微妙ですが、\(r \fallingdotseq n/2\) の場合はこちらのほうが早そうです。
解法3パスカルの三角形を使う
今回の大本命です。
解法2では計算途中で \(4 \equiv 0\) で掛けたり割ったりするから計算できなかったのが悪いので、いっそのこと足し算だけで求めてしまおうという考え。
$$
{}_{n}C_{r}=\begin{cases}
0 &(n<r) \\
1 &(r=0) \\
{}_{n-1}C_{r} + {}_{n-1}C_{r-1} &(\text{otherwise})
\end{cases}
$$
この漸化式を使えば解けます。軽く証明してみましょう。
$$
\begin{aligned}
{}_{n-1}C_{r} + {}_{n-1}C_{r-1}
&=\frac{(n-1)!}{(n-r-1)!r!}+\frac{(n-1)!}{(n-r)!(r-1)!} \\
&=\frac{n-r}{n-r}\frac{(n-1)!}{(n-r-1)!r!}
+\frac{r}{r}\frac{(n-1)!}{(n-r)!(r-1)!} \\
&=(n-r)\frac{(n-1)!}{(n-r)!r!}+r\frac{(n-1)!}{(n-r)!r!} \\
&=\frac{n!}{(n-r)!r!}={}_{n}C_{r} \\
\end{aligned}
$$
ではこれを実装します。
実装1: 再帰関数
漸化式と言えば再帰ですよね。見たまんま実装できるので良いと思います。
n = 2021
r = 37
MOD = 4
def comb(n, r, mod):
if n < r:
return 0
elif r == 0:
return 1
else:
return (comb(n - 1, r, mod) + comb(n - 1, r - 1, mod)) % mod
print(f'{n}C{r} % {MOD} = {comb(n, r, MOD)}')
でもちょっと待ってください。これ、comb関数の呼び出し回数どうなると思います？
$$
\operatorname{call}(n,r)=\begin{cases}
1 &(n<r) \\
1 &(r=0) \\
\operatorname{call}(n-1,r) + \operatorname{call}(n-1,r-1) + 1 &(\text{otherwise})
\end{cases}
$$
として定義した \(\operatorname{call}(n, r)\) 回になるんですよ。つまり少なくとも \({}_{n}C_{r}\) より大きい回数呼び出されます。
大きく見積もれば \(O(2^n)\)。\(n=2021\) で、実行していい計算量ではありません。
実装2: メモ化再帰
上の実装の問題点は、\({}_{n}C_{r}\) を計算するために、\({}_{n-x}C_{r-y}\) の計算を \({}_{x}C_{y}\) 回もしてしまっているということです。
ここで、メモ化という手法を取ることによって、各 \((x,y)\) に対して \({}_{n-x}C_{r-y}\) を \(1\) 回のみ計算すれば良いようにすることができます。
from sys import setrecursionlimit
setrecursionlimit(3000)
n = 2021
r = 37
MOD = 4
memo = [[None] * (r + 1) for _ in range(n + 1)]
def comb(n, r, mod):
if n < r:
return 0
elif r == 0:
return 1
elif memo[n][r] != None:
return memo[n][r]
else:
memo[n][r] = (comb(n - 1, r, mod) + comb(n - 1, r - 1, mod)) % mod
return memo[n][r]
print(f'{n}C{r} % {MOD} = {comb(n, r, MOD)}')
コードを見ればそう難しい話ではないですね。単に二次元配列に計算済みの値を保存して、計算済みならば保存した値を返すようにしているだけです。
計算量は \(O(nr)\) ですが、空間計算量も \(O(nr)\) なのが痛いところです。
実装3: 動的計画法
動的計画法とは、問題を部分問題に分割し、自明な部分問題から遷移していって元の問題を解く手法です。
メモ化再帰もこれに含まれますが、単純再帰は含まれません（動的計画法の要件に計算結果の記録が含まれるため）。
ここでは、ボトムアップに（すなわち、自明な部分問題から順に）求める方法を指します。
パスカルの三角形の \(0\) 行目から \(n\) 行目まで順番に埋めていくイメージですね。
n = 2021
r = 37
MOD = 4
dp = [[0] * (r + 1) for _ in range(n + 1)]
dp[0][0] = 1
for i in range(1, n + 1):
dp[i][0] = 1
for j in range(1, r + 1):
dp[i][j] = (dp[i - 1][j] + dp[i - 1][j - 1]) % MOD
print(f'{n}C{r} % {MOD} = {dp[n][r]}')
これもやはり \(O(nr)\) です。
実装3': SIMD
普通に順次計算すると1クロックで1命令で、計算が終わってまた1命令出すんですが、その1命令に複数のオペランドに対する同一の演算を詰め込むことができます。これがSIMDです。
ベクトルの足し算や内積などが一番よく利用されている例だと思います。
さて、PythonでもSIMDをすることができます。
Numpyというライブラリを利用することで、配列同士の和を取ったりでき、自動でSIMDを使うようになります。内部的に PADDB 命令とかを出してるんでしょうね。
import numpy as np
n = 2021
r = 37
MOD = 4
dp = np.zeros((n + 1, r + 1), np.int8)
dp[0, 0] = 1
for i in range(1, n + 1):
dp[i, 0] = 1
dp[i, 1:] = (dp[i - 1, 1:r + 1] + dp[i - 1, 0:r]) % MOD
print(f'{n}C{r} % {MOD} = {dp[n, r]}')
定数倍の高速化なので、計算量は変わらず \(O(nr)\) です。なお、この条件 \((n,r)=(2021,37)\) だとこっちのほうが遅くなりましたが、\((n,r)=(20000,1000)\) とかにするとめちゃめちゃ高速化できます。
まとめ答えは \(3\) です。
\(n,r\) がともに大きくなった場合は解法2'が最も高速で、\(O(n \log\log n + r \log n)\) で解けます。
高速素因数分解は覚えとけ。
![>]()
内容を見る
