マンデルブロ集合をAVXとOpenMPでできるだけ速く描く #C++

マンデルブロ集合

マンデルブロ集合とは、以下の漸化式で定義される複素数列 z_nにおいて、nを無限大にしたときに発散しない複素数c全体が作る集合です。

\begin{cases}
z_{n+1} = z_{n}^{2} + c \\\
z_0 = 0
\end{cases}

上の画像は、縦が実軸で下が正、横が虚軸で右が正とした複素数平面上にマンデルブロ集合をプロットしたものです。中央のカージオイドとそれに接する無数の円がマンデルブロ集合に属す複素数を表しています。

描画の方法

上の画像の各ピクセルの色は、上記の複素数列の絶対値が第何項でしきい値を超えたかに基づいて決定されています。これをエスケープタイムと呼び、実際にこの数列を計算してこれを求めるアルゴリズムをエスケープタイムアルゴリズムと呼びます。
ちなみに、複素数列のある項の絶対値が2.0以上になる場合は、必ず発散することが知られています。
また、計算機で無限大を扱うことはできないので、一定の回数（下記コードではiter_max回）計算したら打ち切って発散しないものとして扱います。

エスケープタイムアルゴリズムの実装と、それを使ってマンデルブロ集合のアスキーアートを得るプログラムを以下に示します。

#include <stdio.h>

// 複素数 {re, im} のエスケープタイムを求める
// 最大 iter_max 回反復する
int mandelbrot_escape_time(double re, double im, int iter_max) {
	double x = 0.0;
	double y = 0.0;
	for (int i = 0; i < iter_max; i++) {
		double tx = x * x - y * y;
		y = x * y * 2.0 + im;
		x = tx + re;
		if (x * x + y * y >= 4.0) {
			return i;
		}
	}
	return iter_max;
}

int main() {
	int height = 30;
	int width = 60;
	double re_min = -2.0;
	double re_max = 1.0;
	double re_range = re_max - re_min;
	double im_range = re_range * ((double)width / height) / 2;
	double im_min = -im_range / 2;
	double im_max = im_range / 2;
	int iter_max = '~'-' ';
	for (int r = 0; r < height; r++) {
		for (int c = 0; c < width; c++) {
			double re = re_min + re_range / height * r;
			double im = im_min + im_range / width * c;
			printf("%c", (' ' + mandelbrot_escape_time(re, im, iter_max)));
		}
		printf("\n");
	}
    return 0;
}

$ gcc mandelbrot.c -o mandelbrot_aa -O3
$ ./mandelbrot_aa
                  !!!"""""###$~$###"""""!!!                 
             !!!!"""""""#####&~&#####"""""""!!!!            
         !!!!!"""""""""#####%&~&%#####"""""""""!!!!!        
       !!!!!"""""""""""###$$%&~&%$$###"""""""""""!!!!!      
    !!!!!!!"""""""""""##$$$%&'~'&%$$$##"""""""""""!!!!!!!   
  !!!!!!!""""""""""""#$$$$%0+-~-+0%$$$$#""""""""""""!!!!!!! 
!!!!!!!!""""""""""""#$&&&&')-~~~-)'&&&&$#""""""""""""!!!!!!!
!!!!!!!""""""""""""#$%&(/.1~~~~~~~1./(&%$#""""""""""""!!!!!!
!!!!!!"""""""""""##$$$&'+7~~~~~~~~~7+'&$$$##"""""""""""!!!!!
!!!!!!""""""""""###$$$&)B~~~~~~~~~~~B)&$$$###""""""""""!!!!!
!!!!!"""""""""####$$$%&'*8~~~~~~~~~8*'&%$$$####"""""""""!!!!
!!!!!"""""""#####$$%%&&'(*.~~~~~~~.*('&&%%$$#####"""""""!!!!
!!!!"""""""#####$%%%'4*~7~~~~~~~~~~~7~*4'%%%$#####"""""""!!!
!!!!"""""######$&,+)+8~~~~~~~~~~~~~~~~~8+)+,&$######"""""!!!
!!!!""""#####$$%'+~~~~~~~~~~~~~~~~~~~~~~~~~+'%$$#####""""!!!
!!!!""""###$$$%&'i9~~~~~~~~~~~~~~~~~~~~~~~9i'&%$$$###""""!!!
!!!!"""##$$$%*'()~~~~~~~~~~~~~~~~~~~~~~~~~~~)('*%$$$##"""!!!
!!!!!""#2'&'(.~~~]~~~~~~~~~~~~~~~~~~~~~~~~~]~~~.('&'2#""!!!!
!!!!!""#$%(67~~~~j~~~~~~~~~~~~~~~~~~~~~~~~~j~~~~76(%$#""!!!!
!!!!!!"##$~&'/1A,6~~~~~~~~~~~~~~~~~~~~~~~~~6,A1/'&~$##"!!!!!
!!!!!!""###$$$%&(,~~~~~~~~~~~~~~~~~~~~~~~~~,(&%$$$###""!!!!!
!!!!!!!""####$$%%'+~~~~~~~~~~~~~~~~~~~~~~~+'%%$$####""!!!!!!
!!!!!!!!"""###$$%'.~~A~~~~~~~~+~~~~~~~~A~~.'%$$###"""!!!!!!!
!!!!!!!!!""""###$&.(&&(~.->;'&&&';>-.~(&&(.&$###""""!!!!!!!!
!!!!!!!!!!!"""""###$$$$$%$$$$$$$$$$$%$$$$$###"""""!!!!!!!!!!
!!!!!!!!!!!!"""""""""###################"""""""""!!!!!!!!!!!
!!!!!!!!!!!!!!"""""""""""""""""""""""""""""""""!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!"""""""""""""""""""""""""""!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!"""""""""""""""""""!!!!!!!!!!!!!!!!!!!!
$

今回は、描画領域内のすべてのピクセルに対するエスケープタイムを求めるプログラムを作成します。

実験環境

生成範囲等 (冒頭の画像と同じ)
- width: 1920
- height: 1080
- real: -1.5 ~ 0.5
- imag: -1.7777... ~ 1.7777...
- 最大反復回数: 100000
ソフトウェア / ハードウェア
- 言語: C++20
- コンパイラ: gcc version 12.2.0 (Debian 12.2.0-14)
- 最適化オプション: -O3 -mtune=native -march=native -DNDEBUG -fopenmp
- 計算機: ThinkPad X13 Gen3
  - AMD Ryzen7 PRO 6850U (8 cores, 16 threads)
    - 命令セット: AES, AMD-V, AVX, AVX2, FMA3, MMX(+), SHA, SSE, SSE2, SSE3, SSE4.1, SSE4.2, SSE4A, SSSE3, x86-64
  - 32GB LPDDR5

美麗なマンデルブロ集合の画像を得るには、エスケープタイムを色に変換する必要がありますが、その変換は非本質的なのでそれにかかる時間は除外します。また、メモリの確保にかかる時間も除外します。

ナイーブなエスケープタイムアルゴリズム

高速化前

上のCをC++に直し、上記の設定で実行すると、計算時間は95秒でした。

SIMD

複数の命令に対して同じ命令を実行するアーキテクチャをSIMDと呼びます。x86_64のAVX以降のSIMD拡張命令セットは256bitのレジスタを持ちます。これを用いてdouble型を4並列にして計算します。
これをシングルコアで実行したところ、24.3秒(3.9倍速) になりました。4並列になっているので、その効果は十分出ていると言えるでしょう。

OpenMP

このアルゴリズムはピクセルごとが論理的に独立しており、メモリアクセスの局所性も高く、いわゆるEmbarrassingly parallel なアルゴリズムです。
そこでOpenMPを使用してピクセルごとに並列化します。

これだけで、計算時間は14秒(6.7倍速) になりました。
8コアのマシンで実行しているので妥当な範囲と言えると思います。

SIMD + OpenMP

もちろん、SIMD化したものをマルチスレッド実行することもできます。同様にピクセルごとの並列化を行ったところ3.9秒かかりました。

OpenMPだけのときと比べて3.59倍、SIMDだけのときと比べて6.2倍、全体で24.36倍高速です。

境界追跡アルゴリズム

単純なエスケープタイムアルゴリズムで冒頭の画像を生成したとき、発散しないピクセルで実行した反復回数（上記コードのmandelbrot_escape_timeのfor文の実行回数）は439億1430万回、発散するピクセルでは1408万9203回となりました。つまり、全体の99.96%の反復はマンデルブロ集合に属すピクセルの計算で行われています。
そこで、発散するピクセルの計算を可能な限り省略すべく、境界追跡アルゴリズムを導入します。

アルゴリズム

このアルゴリズムは、エスケープタイムが変化する境界部分のみを計算します。以下のようなアルゴリズムです。

描画領域の外周のエスケープタイムを計算する
外周において、隣り合うピクセル同士のエスケープタイムが異なる境界を探し、境界に接しているピクセルをQUEUEにpushする。
QUEUEが空でないとき、以下を繰り返す。
1. QUEUEからピクセルをpopする。これをpと呼ぶ。
2. pのエスケープタイムを計算する。
3. pの8近傍のうち、エスケープタイムが計算されているものと境界が存在するかをチェックする。
4. pの8近傍のうち、QUEUEに追加されておらず、エスケープタイムが計算されていないもののうち、境界に接しているものをQUEUEにpushする。

このアルゴリズムを実行すると、以下の画像で黒い部分が計算されます。

発散する部分の反復回数の合計は10億3420万回となり、439億1430万回のうち97.64%を削っています。

Single Thread & Scaler

シングルスレッドかつスカラーで上記のアルゴリズムを実行すると、3.0秒かかりました。
ナイーブなアルゴリズムと比べると31.7倍、SIMDとOpenMPで高速化したものと比べても1.3倍高速になっています。

SIMD

前述したとおり、AVX2では4並列で計算できるのでQUEUEから4つ一気にpopします。また、同じピクセルをpushすることがないよう、QUEUEに追加したピクセルにはフラグを立てておくことにします。

SIMD化の結果、1.45秒になりました。SIMD化前と比べて2.06倍速になっています。

理論値の4倍に大きく届かない理由は、ナイーブなアルゴリズムの場合は発散しない4ピクセルをパックするケースが比較的多かったのに対して、このアルゴリズムではそのようなケースが減ったことにより、早期に発散するピクセルの無駄な計算が発生することによるものだと思われます。

OpenMP

ナイーブなアルゴリズムではすべてのピクセルが独立に処理できましたが、このアルゴリズムの場合は既存の境界の周辺以外は計算できないため、独立ではありません。また、SIMDの場合は、4つ一気にpopするだけでしたが、この方法はマルチスレッドではスケールしません。なぜなら、queue(std::deque)の操作はatomicに行われる必要があり、その部分で並列度を下げざるを得ないからです。

そこで、描画領域を縦横に分割し、それの一つ一つを各実行スレッドに割り当て、同様のアルゴリズムを実行します。このアルゴリズムではある場所のエスケープタイムを計算したあと、その8近傍しか参照しないため、外周を先に計算することでそれが番兵の役割を果たします。

以下の画像は、描画領域を縦に16分割して計算したものです。16個の短冊状の枠ごとに、外周のエスケープタイムを計算したことがわかると思います。

この方法で並列化を行ったところ、0.65秒となり、ついに1秒を切りました。
並列化前と比べると、4.16倍の高速化となりました。8コアのマシンにしては遅いですが、すべてのコアに平等にタスクが分配されていないことが支配的要因だと思われます。

SIMD + OpenMP

すべての高速化を施すと、0.31秒になりました。
シングルコアかつスカラーによる実装と比べると9.7倍速、OpenMPのみと比べると2.1倍速、SIMDのみと比べると4.7倍速になっています。

他のベンチマーク

下のように適当にズームするとき、各フレームの生成にかかる時間を計測しました。

結果は下のグラフのようになりました。青がシングルかつスカラー、赤がAVX、黄がAVXとOpenMP、緑がAVXのみの効果、水がAVXとOpenMPの効果、橙がAVXと比べたAVXとOpenMPの効果です。
AVXとOpenMPにより、最大で11.4倍程度、平均8.4倍程度の高速化を実現しました。

全体の傾向として、ズームするほどAVXの効果が向上しています。マンデルブロ集合の縁での変化が緩やかになったことで、前述した無駄な計算が緩和されたことによると考えられます。

まとめ

マンデルブロ集合の描画のナイーブなアルゴリズムに対してSIMD化およびマルチスレッド化を行い、24倍程度高速化しました。また、境界探索アルゴリズムを導入し、ナイーブなアルゴリズムに比べて31.7倍程度高速化しました。さらに、境界探索アルゴリズムをSIMD化およびマルチスレッド化し、合計で約300倍程度の高速化を実現しました。

私の計算機で実現できる、一枚絵の計算速度はこれで限界に近いように思います。ここから早くするには以下の方法が考えられます。

インラインアセンブリの利用
- エスケープタイムを求める関数でインラインアセンブリを用いると多少速くなるかもしれません。
GPUの利用
- ただ、ナイーブなアルゴリズムをRTX3060に移植しただけでは勝てないことがわかっています。
AVX512の利用
- 8並列になるだけではなく、様々な機能追加によって多少速くなるかもしれません。
ベクトル内の偏りの調整
- ベクトル内の反復回数の偏りによってSIMDの利用効率が下がることが分かっているので、これをどうにかすると多少速くなる可能性があります。
- AVX512を使う場合にはこれが重要になりそうです。
並列度を上げやすいアルゴリズム
- 今回挙げた境界追跡アルゴリズムの並列化方法は、並列度を上げようとすると外周の計算が増えるため、並列度と計算量がトレードオフの関係にあります。
- このため、GPUのようなメニーコア環境とはあまり相性が良いとは言えません。
- より並列度を上げやすい方法があれば、より大規模な計算機も活用できるでしょう。
FPGAの利用
- 単に速いエスケープタイムルーチンが組めるだけでなく、64ビットを上回る精度で計算できます。
- 乗算器の設計が専門分野の一部なのでいつかやってみたいと思っています。

おまけ: Buddhabrot

ブッダブロはマンデルブロ集合に属さない複素数が作る複素数列z_nの軌跡の確率分布です。上の画像の明るい部分は複素数列がよく通り、逆に暗いところは通らないことを表します。

ブッダブロの描画も、マンデルブロ集合に属す複素数の計算をいかに省けるかが勝負になります。
上の画像は、境界追跡アルゴリズムでマンデルブロ集合を計算しておき、ランダムに生成した複素数が属すピクセルの四隅がマンデルブロ集合に属すときはその複素数もマンデルブロ集合に属すものとして扱っています。

マンデルブロ集合をAVXとOpenMPでできるだけ速く描く