はじめに
競技プログラミング、特にヒューリスティックなコンテストですと、そこまで大きくないグリッド上の最短距離を求めたくなる時が多いかと思います。
本記事では会社の競プロ部の活動の一環として、グリッド上の2点間の最短距離を求める手法についてまとめます。取り扱うのは幅優先探索(Breadth Fast Search)およびDijkstra法についてでして、A*法などは出てきませんのであしからず。
なお本記事ではC++を利用し、以下のusing宣言を仮定します。すみません手癖です。
template<typename T>
using vec = std::vector<T>;
template<typename T>
using vec2 = vec<vec<T>>;
using i2 = std::pair<int,int>;
using clk = std::chrono::steady_clock;
またコンパイルはgcc version 14.2.0 (Homebrew GCC 14.2.0_1) を用いて以下のように行います
g++ -O3 main.cpp -o ./out --std=gnu++20
動作環境はMacBookでチップがAir Apple M3、メモリが16GBです。
移動コストがすべて1の場合
サイズが$L$のグリッドが与えられたとします。グリッドは0と1で表され、1の部分には侵入できないとします。
# L = 16
# [0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0]
# [0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0]
# [0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0]
# [0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0]
# [0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0]
# [0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0]
# [0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0]
# [0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0]
# [1 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0]
# [0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0]
# [0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0]
# [0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0]
# [0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0]
# [0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0]
# [0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0]
# [0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0]
2次元でのBFS
グリッド上の移動コストがすべて1の場合、最短距離は幅優先探索(BFS)で求めることが出来ます。グリッドが2次元vectorで表されたとして、以下のように書けます(雑ですみません……)。
uint32_t bfs2d(const vec2<int>& field, const i2& s, const i2& t){
constexpr uint32_t inf = 1<<30;
static const vec<i2> dirs = {{1,0},{0,1},{-1,0},{0,-1}};
static std::deque<i2> dq;
dq.clear();
dq.emplace_back(s);
static vec2<uint32_t> dist(L, vec<uint32_t>(L, inf));
rep(i,L) std::ranges::fill(dist[i], inf);
const auto& [sx,sy] = s;
const auto& [tx,ty] = t;
dist[sx][sy] = 0;
while(!dq.empty()){
const auto [x,y] = dq.front();
dq.pop_front();
if(x == tx && y == ty) break;
for(const auto& [dx,dy] : dirs){
int nx = x + dx;
int ny = y + dy;
if(nx < 0 || L <= nx || ny < 0 || L <= ny) continue;
if(field[nx][ny]) continue;
if(dist[x][y] + 1 < dist[nx][ny]){
dist[nx][ny] = dist[x][y] + 1;
dq.emplace_back(nx,ny);
}
}
}
return dist[tx][ty];
}
1次元でのBFS
さて、一般に2次元配列をそのまま使うのはアクセス速度などの関係で遅いことが多く、高速化のために1次元化してから探索するテクニックがあります。例えば下記のようにグラフ構造で表して
vec2<int> make1dGraph(vec2<int>& field){
vec2<int> result(L * L);
static const vec<i2> dirs = {{1,0},{0,1},{-1,0},{0,-1}};
rep(x,size){
rep(y,size){
int v = L * x + y;
for(const auto& [dx,dy] : dirs){
int nx = x + dx;
int ny = y + dy;
if(nx < 0 || L <= nx || ny < 0 || L <= ny) continue;
if(field[nx][ny] == 1) continue;
int nv = L * nx + ny;
result[v].push_back(nv);
}
}
}
return result;
}
探索します。
uint32_t bfs1d(const vec2<int>& graph, const int s, const int t){
constexpr uint32_t inf = 1<<30;
static std::deque<int> dq;
dq.clear();
dq.push_back(s);
static vec<uint32_t> dist(L*L);
std::ranges::fill(dist, inf);
dist[s] = 0;
while(!dq.empty()){
const auto v = dq.front();
dq.pop_front();
if(v==t) break;
for(const auto nv : graph[v]){
if(dist[v] + 1 < dist[nv]){
dist[nv] = dist[v] + 1;
dq.push_back(nv);
}
}
}
return dist[t];
}
ビッドボードを使ったBFS
その他、ビットボードを使うテクニックがあります。AtCoder社様が主催するAtCoder Heuristic Constest 046に関するこちらの解説スライドを実装の参考にさせて頂きました。
まずグリッドをビットマップで表します。
constexpr size_t BIT_SIZE = L * L;
std::bitset<BIT_SIZE> makeBitBoard(const vec2<int>& field){
std::bitset<BIT_SIZE> result;
rep(x, L){
rep(y, L){
if(field[x][y] == 1){
result.set(x * L + y);
}
}
}
return result;
}
実際は通れるグリッドを知りたいのでフリップしたものを使います。
auto bit_board = makeBitBoard(field);
auto free_bit_board = bit_board;
free_bit_board.flip();
解説スライドにあります左右のガードも用意します。
auto right_guard = std::bitset<BIT_SIZE>();
auto left_guard = std::bitset<BIT_SIZE>();
rep(x,L){
rep(y,L){
int v = L*x + y;
if(y < L-1) right_guard.set(v);
if(0 < y) left_guard.set(v);
}
}
そしてこれらを使ってbfsします。
uint32_t bfsBit(const std::bitset<BIT_SIZE>& free_bit_board, int s, int t){
if(s==t) return 0;
static auto current = std::bitset<BIT_SIZE>();
current.reset();
current[s] = true;
rep(d,BIT_SIZE){
auto next = current;
next |= (current & right_guard)<<1; // right move
next |= (current & left_guard)>>1; // left move
next |= (current >> size); // upper move
next |= (current << size); // lower move
next &= free_bit_board;
if(current == next){
break;
}
if(next[t]) return d+1;
std::swap(current, next);
}
return 1<<30;
}
速度調査
さて、では速度調査をします。速度調査では$M=10000$として$M$個の始点と終点のペア$(start, end)$を用意し、それらすべての最短距離を下記のように計算します。
long long result = 0;
auto start_time = clk::now();
for(const auto& [start, end] : queries){
result += bfs2d(field, start, end);
}
auto end_time = clk::now();
結果は以下のようになりました(単位:秒)。グリッドのサイズが小さい($L\leq 32$)時のビットボードの速さが際立ちますね。また2D版と1D版の差はサイズが大きくなるほどに縮まるようでした。
| 関数\L | 16 | 32 | 64 | 128 |
|---|---|---|---|---|
| bfs2d | 0.011450 | 0.036547 | 0.140537 | 0.509499 |
| bfs1d | 0.006311 | 0.026202 | 0.116278 | 0.444023 |
| bfsBit | 0.001316 | 0.005523 | 0.108464 | 0.761697 |
移動コストが非負の場合
次に移動コストが非負の場合を考えます。グリッド間のパスにコストを置いても良いですが、今回は簡単のためにグリッドそれぞれにコスト$c_{ij}$を置いてそこに移動するとそのコストが掛かるとします。グリッド全体は以下のような感じで、壁は無しとします。
L=16, コスト最大: 10
# [7 2 9 9 8 2 2 1 9 4 2 8 1 7 3 6]
# [1 5 5 9 9 9 9 6 4 3 8 7 6 6 0 7]
# [9 4 2 0 2 8 7 3 4 8 4 3 3 5 0 7]
# [5 3 3 6 1 0 2 2 0 7 3 5 8 4 9 0]
# [1 8 2 2 1 9 1 4 0 7 9 0 8 2 7 0]
# [6 1 9 0 7 1 1 5 0 6 1 3 7 9 1 6]
# [5 8 9 1 4 8 6 8 8 5 6 3 5 1 0 3]
# [2 5 5 8 5 2 8 8 2 0 1 4 3 8 1 6]
# [7 1 1 6 8 6 5 1 4 5 5 3 2 4 7 4]
# [4 5 3 7 6 4 2 7 6 5 3 7 9 9 5 4]
# [2 9 4 9 4 0 2 5 9 8 7 9 3 5 7 4]
# [1 4 8 9 2 6 2 9 6 3 5 4 8 1 1 1]
# [2 3 0 1 5 8 5 0 5 8 3 2 2 3 0 8]
# [9 3 5 7 6 7 2 8 6 8 4 5 6 7 2 9]
# [6 5 7 8 7 6 3 6 2 9 7 7 8 1 2 5]
# [5 7 7 9 3 3 3 0 2 3 9 9 8 7 1 7]
Dijkstra法
移動コストが非負の場合は、Dijkstra法が使えます。std::priority_queueを用いた実装が以下になります。
uint32_t Dijkstra(const vec2<i2>& graph, const int s, const int t){
constexpr uint32_t inf = 1<<30;
std::priority_queue<i2, vec<i2>, std::greater<i2>> hq;
hq.emplace(0, s);
static vec<uint32_t> dist(L*L);
std::ranges::fill(dist, inf);
dist[s] = 0;
while(!hq.empty()){
const auto [c,v] = hq.top(); // ここを参照にすると不具合するので注意
hq.pop();
if(v==t) break;
if(dist[v] < c) continue;
for(const auto [nv,dc] : graph[v]){
if(dist[v] + dc < dist[nv]){
dist[nv] = dist[v] + dc;
hq.emplace(dist[nv], nv);
}
}
}
return dist[t];
}
Radix Heap使用版
Dijkstra法の速度は現時点での最短頂点をどう取り出すかが肝でして、std::priority_queueよりその操作が高速なデータ構造を使うことでより速いDijksta法が実現できます。今回はこちらのサイトで紹介されているRadix Heapの実装をほぼそのまま拝借させて頂きました![]()
// from https://ei1333.github.io/algorithm/radix-heap.html (Thanks!)
template< typename T >
struct RadixHeap
{
using uint = unsigned;
std::vector< std::pair< uint, T > > v[33];
uint size, last;
RadixHeap() : size(0), last(0) {}
bool empty() const { return size == 0; }
inline int getbit(int a)
{
return a ? 32 - __builtin_clz(a) : 0;
}
void push(uint key, const T &value)
{
++size;
v[getbit(key ^ last)].emplace_back(key, value);
}
std::pair< uint, T > pop()
{
if(v[0].empty()) {
int idx = 1;
while(v[idx].empty()) ++idx;
last = std::min_element(std::begin(v[idx]), std::end(v[idx]))->first;
for(auto &p : v[idx]) v[getbit(p.first ^ last)].emplace_back(p);
v[idx].clear();
}
--size;
auto ret = v[0].back();
v[0].pop_back();
return ret;
}
void clear(){
size = 0; last = 0;
rep(i,33) v[i].clear();
}
};
uint32_t DijkstraByRadixHeap(const vec2<i2>& graph, const int s, const int t){
constexpr uint32_t inf = 1<<30;
static RadixHeap<int> hq;
hq.clear();
hq.push(0, s);
static vec<uint32_t> dist(L*L);
std::ranges::fill(dist, inf);
dist[s] = 0;
while(!hq.empty()){
const auto [c,v] = hq.pop();
if(v==t) break;
if(dist[v] < c) continue;
for(const auto [nv,dc] : graph[v]){
if(dist[v] + dc < dist[nv]){
dist[nv] = dist[v] + dc;
hq.push(dist[nv], nv);
}
}
}
return dist[t];
}
0...kBFS利用版
一般に非負コストの最短距離ですとDijkstra法ですが、例えば移動コストが0 or 1の場合は01BFSと呼ばれる
手法でより高速に最短距離を求めることが出来ます。距離0と1に対応するキューを用意する感じです(実際はひとつのDequeで表現されることが多いです)。
同様の考えで、最短距離がさほど大きくならないことが保証出来るならば、それぞれの距離に応じたキューみたいなものを用意しておけば高速に探索できるのではという話です。正式名称が分からないので本記事ではとりあえず0..kBFSと書いておきます。
(どの距離がキューに溜まっているかを高速に知るためにまたRadixHeap使ってます)
// ref: https://yosupo.hatenablog.com/entry/2015/04/03/224649 (Thanks!)
struct RadixHeapInt {
using uint = unsigned int;
vec<uint> v[33];
uint last, sz;
RadixHeapInt() {
last = sz = 0;
}
int bsr(uint x){
if (x==0) return -1;
return 31 - __builtin_clz(x);
}
void push(uint x){
sz++;
v[bsr(x^last)+1].push_back(x);
}
uint top() { // 追加
if (!v[0].size()) {
int i = 1;
while (!v[i].size()) i++;
uint new_last = *std::min_element(v[i].begin(), v[i].end());
for(uint x : v[i]){
v[bsr(x^new_last)+1].push_back(x);
}
last = new_last;
v[i].clear();
}
return last;
}
uint pop() {
if (!v[0].size()) {
int i = 1;
while (!v[i].size()) i++;
uint new_last = *std::min_element(v[i].begin(), v[i].end());
for(uint x : v[i]){
v[bsr(x^new_last)+1].push_back(x);
}
last = new_last;
v[i].clear();
}
sz--;
v[0].pop_back();
return last;
}
void clear(){
for(auto& sub : v) sub.clear();
last = sz = 0;
}
};
struct KQueue{ // queue?
vec2<int> qs;
int max_size;
RadixHeapInt index_set;
KQueue(int max_size):max_size(max_size){
qs = vec2<int>(max_size);
}
void push(int c, int s){
if(qs[c].empty()){
index_set.push(c);
}
qs[c].push_back(s);
}
i2 pop(){
int c = index_set.top();
int res = qs[c].back();
qs[c].pop_back();
if(qs[c].empty()){
index_set.pop();
}
return {c,res};
}
bool empty(){
return index_set.sz == 0;
}
void clear(){
while(index_set.sz > 0){
int c = index_set.pop();
qs[c].clear();
}
index_set.clear();
}
};
uint32_t DijkstraByKBFS(const vec2<i2>& graph, const int s, const int t){
constexpr uint32_t inf = 1<<30;
static auto kq = KQueue(500);
kq.clear();
kq.push(0, s);
static vec<uint32_t> dist(L*L);
std::ranges::fill(dist, inf);
dist[s] = 0;
while(!kq.empty()){
const auto [c,v] = kq.pop();
if(v==t) break;
if(dist[v] < c) continue;
for(const auto [nv,dc] : graph[v]){
if(dist[v] + dc < dist[nv]){
dist[nv] = dist[v] + dc;
kq.push(dist[nv], nv);
}
}
}
return dist[t];
}
速度調査
では速度調査をします。速度調査の仕方は先程と同様です($M=10000$)。ただしコストの最大値$MC$に対して、キュー(?)はとりあえず$L*MC$本用意しています。
| 関数\(L,M) | (16,10) | (32,10) | (64,10) | (64,30) |
|---|---|---|---|---|
| priority_queue | 0.01892 | 0.215918 | 1.09289 | 1.049 |
| RadixHeap版 | 0.010245 | 0.077246 | 0.417854 | 0.47432 |
| 0..kBFS版 | 0.008369 | 0.054789 | 0.316994 | 0.365121 |
こう見るとなんか0..kBFSが速い感じがしますね……理屈は正直考えれてないです。実装起因も大きいかもですね。
まとめ
以上、グリッド上の最短距離の求め方の紹介と、それぞれの速度を調べてみました。ビットボードやRadix Heapなどは今まで試せてなかったので良い機会になりました。解説記事を書いてくださっていた方に感謝です。
今回は終点を明示的に指定しましたが"特定の条件を満たすマスを見つけるまで"みたいなのにも使えるかと思いますので、ぜひまた活用し高速化していきましょう。
ここまでお読み頂きありがとうございました!