CPUを用いたFizzBuzzの高速化について

2022年4月1日に、κeenさんがFizzBuzzの高速化について技術記事を投稿なさっていました。

https://keens.github.io/blog/2022/04/01/gpuwomochiitafizzbuzznokousokukanotorikumi/

この記事はあくまでジョーク記事という扱いですが、筆者もFizzBuzzの高速化に興味が湧いたので、実装してみることにしました。

## FizzBuzz問題とは

FizzBuzzはプログラミングのお題として採り上げられる問題の1つです。
https://ja.wikipedia.org/wiki/Fizz_Buzz

この記事では次のような問題だとして扱うことにします。

- 整数を表すNが引数として与えられる
- 1からNまで順番に「その数が3で割り切れるならfizzを」「その数が5で割り切れるならbuzzを」「その数が3でも5でも割り切れるならfizzbuzzを」「どれにも当てはまらないならその数自身を」末尾に改行文字（LF）を付与して出力する
- 出力先は標準出力とする

これをひとまず実装するだけならそれほど難しくありません。例えばPythonで書いてみます。

```python:fizzbuzz.py
import sys

N = int(sys.argv[1])
for i in range(1, N + 1):
    if i % 15 == 0:
        print("fizzbuzz")
    elif i % 3 == 0:
        print("fizz")
    elif i % 5 == 0:
        print("buzz")
    else:
        print(i)
```

実行例です。

```console
$ python fizzbuzz.py 20
1
2
fizz
4
buzz
fizz
7
8
fizz
buzz
11
fizz
13
14
fizzbuzz
16
17
fizz
19
buzz
```

今回はこのfizzbuzz問題の高速化を考えます。以降は特に、ベンチマークの条件として「1から1億までのfizzbuzzの出力」に絞って考えることにしました。1億を選んだ理由ですが、実行時間が数秒程度で長すぎず短すぎず、出力も現実的なサイズだったので、取り扱いやすかったからです。

## 評価環境
これからfizzbuzzのコードを書き、性能の測定・改善をしていくのですが、ベンチマーク等は以下の環境で行いました。

| 項目 | 内容 |
|-----|-----|
| ハードウェア | AWSのc6i.largeインスタンス |
| OS | Ubuntu 22.04 LTS |
| コンパイラ | GCC 11.3.0 |
| コンパイルオプション | `--std=c++17 -O2 -g3` |

## 実装

実装はC++で行いました。ただし、実装時間の関係上、色々な箇所で手抜きを行っています。あまりきれいなコードではないと思いますのでご注意ください。

- 全体的にC++17に準拠していますが、一部に環境依存のもの（例えばシステムコール`write`をそのまま呼んでいます）や、プラクティスとして妥当か微妙な部分（例えばC由来の関数を呼んでいる等）があるかもしれません
- エラーハンドリングは全体的に適当です
- 記事中のコードスニペットの中にはincludeを省略している箇所があります
- ベンチマークとして1億までのfizzbuzzを考えることにしました。この条件のもとで手抜きしている箇所があります
  - 例えば、十進数で9桁までの整数しか扱わないと実装上仮定している箇所があります

### 最初の実装

チューニングを頑張る前に、とりあえず普通にfizzbuzzを書いてみることにしました。

```cpp
int main(int argc, char *argv[]) {
  assert(argc == 2);
  const std::string arg_num{argv[1]};
  const auto N = std::stoll(arg_num);
  for (long long i = 1; i <= N; i++) {
    if (i % 15 == 0) {
      puts("fizzbuzz");
    } else if (i % 3 == 0) {
      puts("fizz");
    } else if (i % 5 == 0) {
      puts("buzz");
    } else {
      printf("%lld\n", i);
    }
  }
}
```

ざっとベンチマークを取ることで雰囲気を把握します。

```console
$ /usr/bin/time -p ./main_naive 100000000 >/dev/null
real 3.51
user 3.49
sys 0.01
```

1億までのfizzbuzzは約700MiB[^100m-fizzbuzz]になります。これを3.51sで処理しているので出力のスループットとしては199MiB/sほどです。あまり速くはないですね。sysがおよそ0.3sと比較的小さいので`write`（バッファから標準出力への書き出し）はあまり問題になっておらず、単純にプログラムのどこかが遅いようです。

[^100m-fizzbuzz]: 小数点第二位を四捨五入すると700.1MiBのようです

さらに、後のことを考えて「引数`n0`から`num`個のfizzbuzz列を生成し、`dst`が指す先に書く」バージョンの実装も準備しておきます。

```cpp
size_t write_fizzbuzz_block_naive(char *dst, int64_t n0, int64_t num) {
  constexpr std::string_view fizz{"fizz\n"};
  constexpr std::string_view buzz{"buzz\n"};
  constexpr std::string_view fizzbuzz{"fizzbuzz\n"};
  const char *dst0 = dst;
  for (int64_t i = n0; i < n0 + num; i++) {
    if (i % 15 == 0) {
      fizzbuzz.copy(dst, fizzbuzz.size());
      dst += fizzbuzz.size();
    } else if (i % 3 == 0) {
      fizz.copy(dst, fizz.size());
      dst += fizz.size();
    } else if (i % 5 == 0) {
      buzz.copy(dst, buzz.size());
      dst += buzz.size();
    } else {
      dst += sprintf(dst, "%PRId64\n", i);
    }
  }
  return dst - dst0;
}
```

この実装もそれほど速くないようです。

一般に、コードのどの部分が遅いか調査するには、プロファイラを当てることになると思います。しかし、今回の場合は可能性がありそうな場所は限られます。というか(s)printfぐらいしかないですね。たとえばsprintfをコードから消すことで、実行が大幅に高速になることは簡単に調べられます。

### 整数の文字列化

さて、sprintfがなぜ遅いかは調査が難しい（とても高機能な関数なので……）ので、今回の用途に特化した関数を自作することにします。

結局やりたいことは整数を文字列に変換することですが、もちろん、この手法はよく研究されています。Daniel Lemireの説明する素朴な技法は、要するに10進数でいう2桁ごとの値を計算し、事前定義したテーブルから引く方法です。

https://lemire.me/blog/2022/03/28/converting-integers-to-decimal-strings-faster-with-avx-512/

今回は、先の考察にあるように、出力させる桁数が予め確定させられる点も実装に活用します。具体的には、テンプレートの引数で桁数を渡してやり、計算する桁数を定数としてコンパイラの最適化を効かせます。

```cpp
template <int d> void to_string_tree_table(uint64_t x, char *out) {
  static const char table[200] = {
      0x30, 0x30, 0x30, 0x31, 0x30, 0x32, 0x30, 0x33, 0x30, 0x34, 0x30, 0x35,
      0x30, 0x36, 0x30, 0x37, 0x30, 0x38, 0x30, 0x39, 0x31, 0x30, 0x31, 0x31,
      0x31, 0x32, 0x31, 0x33, 0x31, 0x34, 0x31, 0x35, 0x31, 0x36, 0x31, 0x37,
      0x31, 0x38, 0x31, 0x39, 0x32, 0x30, 0x32, 0x31, 0x32, 0x32, 0x32, 0x33,
      0x32, 0x34, 0x32, 0x35, 0x32, 0x36, 0x32, 0x37, 0x32, 0x38, 0x32, 0x39,
      0x33, 0x30, 0x33, 0x31, 0x33, 0x32, 0x33, 0x33, 0x33, 0x34, 0x33, 0x35,
      0x33, 0x36, 0x33, 0x37, 0x33, 0x38, 0x33, 0x39, 0x34, 0x30, 0x34, 0x31,
      0x34, 0x32, 0x34, 0x33, 0x34, 0x34, 0x34, 0x35, 0x34, 0x36, 0x34, 0x37,
      0x34, 0x38, 0x34, 0x39, 0x35, 0x30, 0x35, 0x31, 0x35, 0x32, 0x35, 0x33,
      0x35, 0x34, 0x35, 0x35, 0x35, 0x36, 0x35, 0x37, 0x35, 0x38, 0x35, 0x39,
      0x36, 0x30, 0x36, 0x31, 0x36, 0x32, 0x36, 0x33, 0x36, 0x34, 0x36, 0x35,
      0x36, 0x36, 0x36, 0x37, 0x36, 0x38, 0x36, 0x39, 0x37, 0x30, 0x37, 0x31,
      0x37, 0x32, 0x37, 0x33, 0x37, 0x34, 0x37, 0x35, 0x37, 0x36, 0x37, 0x37,
      0x37, 0x38, 0x37, 0x39, 0x38, 0x30, 0x38, 0x31, 0x38, 0x32, 0x38, 0x33,
      0x38, 0x34, 0x38, 0x35, 0x38, 0x36, 0x38, 0x37, 0x38, 0x38, 0x38, 0x39,
      0x39, 0x30, 0x39, 0x31, 0x39, 0x32, 0x39, 0x33, 0x39, 0x34, 0x39, 0x35,
      0x39, 0x36, 0x39, 0x37, 0x39, 0x38, 0x39, 0x39,
  };
  uint64_t top = x / 100000000;
  uint64_t bottom = x % 100000000;
  uint64_t toptop = top / 10000;
  uint64_t topbottom = top % 10000;
  uint64_t bottomtop = bottom / 10000;
  uint64_t bottombottom = bottom % 10000;
  uint64_t toptoptop = toptop / 100;
  uint64_t toptopbottom = toptop % 100;
  uint64_t topbottomtop = topbottom / 100;
  uint64_t topbottombottom = topbottom % 100;
  uint64_t bottomtoptop = bottomtop / 100;
  uint64_t bottomtopbottom = bottomtop % 100;
  uint64_t bottombottomtop = bottombottom / 100;
  uint64_t bottombottombottom = bottombottom % 100;

  static_assert(d <= 10); // とりあえず今回必要な桁数に限定する
  switch (d) {
  case 1:
    memcpy(out, &table[2 * bottombottombottom] + 1, 1);
    break;
  case 2:
    memcpy(out, &table[2 * bottombottombottom], 2);
    break;
  case 3:
    memcpy(out, &table[2 * bottombottomtop] + 1, 1);
    memcpy(out + 1, &table[2 * bottombottombottom], 2);
    break;
  case 4:
    memcpy(out, &table[2 * bottombottomtop], 2);
    memcpy(out + 2, &table[2 * bottombottombottom], 2);
    break;
  case 5:
    memcpy(out, &table[2 * bottomtopbottom] + 1, 1);
    memcpy(out + 1, &table[2 * bottombottomtop], 2);
    memcpy(out + 3, &table[2 * bottombottombottom], 2);
    break;
  case 6:
    memcpy(out, &table[2 * bottomtopbottom], 2);
    memcpy(out + 2, &table[2 * bottombottomtop], 2);
    memcpy(out + 4, &table[2 * bottombottombottom], 2);
    break;
  case 7:
    memcpy(out, &table[2 * bottomtoptop] + 1, 1);
    memcpy(out + 1, &table[2 * bottomtopbottom], 2);
    memcpy(out + 3, &table[2 * bottombottomtop], 2);
    memcpy(out + 5, &table[2 * bottombottombottom], 2);
    break;
  case 8:
    memcpy(out, &table[2 * bottomtoptop], 2);
    memcpy(out + 2, &table[2 * bottomtopbottom], 2);
    memcpy(out + 4, &table[2 * bottombottomtop], 2);
    memcpy(out + 6, &table[2 * bottombottombottom], 2);
    break;
  case 9:
    memcpy(out, &table[2 * topbottombottom] + 1, 1);
    memcpy(out + 1, &table[2 * bottomtoptop], 2);
    memcpy(out + 3, &table[2 * bottomtopbottom], 2);
    memcpy(out + 5, &table[2 * bottombottomtop], 2);
    memcpy(out + 7, &table[2 * bottombottombottom], 2);
    break;
  case 10:
    memcpy(out, &table[2 * topbottombottom], 2);
    memcpy(out + 2, &table[2 * bottomtoptop], 2);
    memcpy(out + 4, &table[2 * bottomtopbottom], 2);
    memcpy(out + 6, &table[2 * bottombottomtop], 2);
    memcpy(out + 8, &table[2 * bottombottombottom], 2);
    break;
  }
}

template <int d> size_t write_number(uint64_t x, char *out) {
  to_string_tree_table<d>(x, out);
  out[d] = '\n';
  return d + 1;
}
```

今回の用途では、整数を文字列化して出力するだけでなく末尾に改行文字の付与が必要なので、それをする`write_number`も定義しています。

この実装とsprintfで性能比較のマイクロベンチマークを書いてみます。桁数を10に固定した上で、固定の値の文字列化をひたすら計算させます。

```cpp:bench_i2s.cpp
int main() {
  constexpr auto N = 20 * 1000 * 1000ll;
  const auto x = 1234567890ll;
  char out[16] = {'\0'};
  // sprintf
  time([&]() -> void {
    for (int i = 0; i < N; i++) {
      sprintf(out, "%lld\n", x);
    }
  });
  // 専用実装
  time([&]() -> void {
    for (int i = 0; i < N; i++) {
      write_number<10>(x, out);
    }
  });
}
```

ただし、コードの特定区間で実行時間を取るために、以下の`time`関数を用いました。

```cpp
#include <chrono>
#include <functional>
#include <stdio.h>

inline void time(std::function<void()> f) {
  using std::chrono::duration;
  using std::chrono::high_resolution_clock;
  using std::chrono::milliseconds;

  auto t1 = high_resolution_clock::now();
  f();
  auto t2 = high_resolution_clock::now();
  duration<double, std::milli> ms_double = t2 - t1;
  fprintf(stderr, "%f[ms]\n", ms_double.count());
}
```

実行例です。

```console
$ g++ --std=c++17 -O2 -o bench bench_i2s.cpp
$ ./bench 
1049.527316[ms]
17.328950[ms]
```

今回の環境では60倍ほど高速化できているようです。

### 元の関数を書き換える

整数を文字列に変換する実装を得たので、ようやく元の関数`write_fizzbuzz_block_naive`を書き換える準備が整いました。改良版である`write_fizzbuzz_block_opt`を以下に提示します。

```cpp
size_t write_fizzbuzz_block_opt(char *dst, const long long n0, const long long num) {
  const char *dst0 = dst;
  const auto n1 = n0 + num - 1;
  assert(n1 <= 999'999'999ll); // これ以上は未実装
  // 4桁目までは行数としても少ないので、基本の実装（write_fizzbuzz_block_naive）に流す
  if (n0 <= 9'999ll) {
    const auto x = std::max(0ll, n0);
    const auto y = std::min(9'999ll, n1);
    dst += write_fizzbuzz_block_naive(dst, x, y - x + 1);
  }
  // 5桁部分
  if (n0 <= 99'999ll && 10'000ll <= n1) {
    const auto x = std::max(10'000ll, n0);
    const auto y = std::min(99'999ll, n1);
    dst += write_fizzbuzz_block_d<5>(dst, x, y - x + 1);
  }
  // 6桁部分
  if (n0 <= 999'999ll && 100'000ll <= n1) {
    const auto x = std::max(100'000ll, n0);
    const auto y = std::min(999'999ll, n1);
    dst += write_fizzbuzz_block_d<6>(dst, x, y - x + 1);
  }
  // 7桁部分
  if (n0 <= 9'999'999ll && 1'000'000ll <= n1) {
    const auto x = std::max(1'000'000ll, n0);
    const auto y = std::min(9'999'999ll, n1);
    dst += write_fizzbuzz_block_d<7>(dst, x, y - x + 1);
  }
  // 8桁部分
  if (n0 <= 99'999'999ll && 10'000'000ll <= n1) {
    const auto x = std::max(10'000'000ll, n0);
    const auto y = std::min(99'999'999ll, n1);
    dst += write_fizzbuzz_block_d<8>(dst, x, y - x + 1);
  }
  // 9桁部分
  if (n0 <= 999'999'999ll && 100'000'000ll <= n1) {
    const auto x = std::max(100'000'000ll, n0);
    const auto y = std::min(999'999'999ll, n1);
    dst += write_fizzbuzz_block_d<9>(dst, x, y - x + 1);
  }
  return dst - dst0;
}
```

やっていることは単純なのですが見た目が複雑なので順番に考えていきます。

まず引数は「引数`n0`から`num`個のfizzbuzz列を生成し、`dst`が指す先に書く」という意味でした。出力するべき区間の末尾を`n1`と書くことにすると、`n1 = n0 + num -1`です。例えば1から200個のfizzbuzzを出力するとき、末尾は1+200-1 = 200だと考えると分かりやすいかもしれません[^fizzbuzz-200]。

[^fizzbuzz-200]: fizzbuzzのルールにより、実際の出力は200ではなく"buzz"になります。

つまり区間[n0, n1]に含まれる整数をfizzbuzzルールに従って出力することを考えればよいのですが、前述の通り出力桁数は予め固定にしなければなりません。今回の問題の範囲では、桁数は1から9までを考えれば十分なので、その範囲で各桁について場合分けを行っています。

区間[n0, n1]が桁数dを含んでいるかは、その桁数の区間（例えば3桁だと[100, 999]）と重なり判定をしていると見なすと分かりやすいです。

https://koseki.hatenablog.com/entry/20111021/range

各桁それぞれの場合において、その桁数の数字は区間中のどの一部分にあたるかを計算しているのが`x, y`の定義部分です。これは要するに桁数の区間と共通部分を取ればよいのでmin, max関数を使って計算しています。

### 改善その1: sprintfの差し替え

ここまでで改良した実装を使ってベンチマークを取ります。

```cpp:main_opt1.cpp
int main(int argc, char *argv[]) {
  assert(argc == 2);
  const std::string arg_num{argv[1]};
  const auto N = std::stoll(arg_num);
  auto buff = new std::array<char, 512 * 1000 * 1000>();
  auto s = write_fizzbuzz_block_opt(buff->data(), 1, N);
  write(1, buff->data(), s);
}
```

実行例です。

```console
$ /usr/bin/time -p ./main_opt1 100000000 >/dev/null 
real 0.62
user 0.36
sys 0.26
```

数倍程度の高速化が見られますが、userとsysを比較してみたとき相対的にsysの大きさが目立ちます。そこで、次節ではIO周りについて改善の余地がないか考えることにしました。


### 改善その2: 出力の分割
ここまでの実装では、1からNまでのfizzbuzz列を全てメインメモリに出力してからシステムコール`write`による書き出しを行っていました。

しかし、fizzbuzzの出力において全体をメモリ上に保持しておく必要性はないので、メモリの無駄使いです。無駄なだけならまだいいとしても、より細かな単位で出力したほうが参照の局所性も改善するはずで、性能向上が期待できるでしょう。

というわけで、1からNまでのN要素を、ほどほどの区間`BLOCK_SIZE`で分割した上で、個々の区間ごとに`write`していくことにします。

```cpp
int main(int argc, char *argv[]) {
  assert(argc == 2);
  const std::string arg_num{argv[1]};
  const auto N = std::stoll(arg_num);
  auto buff = new std::array<char, 32 * 1024 * 1024>();
  constexpr auto BLOCK_SIZE = 100 * 1000ll;
  int64_t i;
  // BLOCK_SIZEごとにfizzbuzz計算をして、結果をwriteに渡している
  for (i = 1; i + BLOCK_SIZE <= N + 1; i += BLOCK_SIZE) {
    auto s = write_fizzbuzz_block_opt(buff->data(), i, BLOCK_SIZE);
    write(1, buff->data(), s);
  }
  // loop remainder
  auto s = write_fizzbuzz_block_opt(buff->data(), i, N - i + 1);
  write(1, buff->data(), s);
}
```

実行例です。

```console
$ /usr/bin/time -p ./main_opt2 100000000 >/dev/null 
real 0.34
user 0.33
sys 0.01
```

期待通りsysが改善している様子があります。

## 結果のまとめ

ここまでの結果をまとめます。

| 実装    | 実行時間（秒） | スループット（MiB/s） |
|--------|--------------|------------|
| 最初の実装 | 3.51 | 199 |
| 改善その1 | 0.62 | 1130 |
| 改善その1 + 改善その2 | 0.34 | 2060 |

最初は素朴な実装から始まりましたが、2つの改良を入れることで性能は約10倍ほど改善し、スループットは2GiB/sほどにまでなりました。

## さらなる高速化は可能か？

いくつかの可能性は思いつくので、今回の試みがアプローチの全てではないように思われます。

- BLOCK_SIZEなどいくつかのパラメータをより適切な値にチューニングする
  - 今回は雑に値を決め打ちしたので最適値からは遠い可能性がある
- マルチスレッド化
- アルゴリズム的な工夫の余地もまだあるかもしれない
  - 整数の文字列変換アルゴリズムは、先の記事に記載があるようにAVX-512対応が可能らしい？
  - FizzBuzzは15の倍数で出力パターンがループする点を上手く活用できないだろうか？

## おわりに

FizzBuzzはプログラマの間でおなじみの有名問題ですが、これの高速化（？）に取り組んだ事例はκeenさんのものしか知らなかったので、今回試しに取り組んでみました。

執筆時間の関係で各項目について深堀りできていない点があります。例えば、sprintfをとりあえずで差し替えましたが、本来ちゃんとライブラリの中の実装まで追うことで挙動把握したいものです。また、元ネタの記事はGPUなことを踏まえると、CPU実装でもマルチスレッド化でさらなる高速化が見込めるのかは気になる点ですが検証できていません。

これら不完全な点はありますが、今回試行錯誤した範囲の工夫でもある程度の性能向上を実現できたので、結果にはそれなりに満足しています。