同じ rendering ループを 2 通りで実装した: メインスレッド上の vanilla JavaScript と、まったく同じループを書いた GLSL fragment shader。浅いズームでは GPU が 30–100 倍速い。だが あるスケールを超えると、GPU は CPU と明らかに異なる — そして間違った — 像を返し始める。本記事はその「あるスケール」がどこで、なぜ起きるかの話。
📦 GitHub: https://github.com/sen-ltd/webgl-mandelbrot
🧮 Demo: https://sen.ltd/portfolio/webgl-mandelbrot/
スクリーンショットはデモの float32 breaks プリセット: 同じビューポートを 512×512、500 iterations で 2 度 render している。左パネルが JavaScript CPU レンダラ、右が WebGL GPU レンダラ、完全に同一の escape-time カーネル を走らせている。違って見えるのは、片方が動いていて片方が壊れているから。
前提
Mandelbrot 集合は z ← z² + c を z = 0 から反復することで定義される。|z| が永遠に有界なら c は集合に属する。実装上は |z|² > 4 (証明済の脱出半径) になるか反復回数が上限に達するまで回す。脱出時の反復回数で色付け。
内側ループは tiny:
// JavaScript、1 ピクセル
export function escapeCount(cx, cy, maxIter) {
let x = 0, y = 0;
let x2 = 0, y2 = 0;
for (let i = 0; i < maxIter; i++) {
if (x2 + y2 > 4) return i;
y = 2 * x * y + cy;
x = x2 - y2 + cx;
x2 = x * x;
y2 = y * y;
}
return maxIter;
}
GLSL 版は同じループを 1 単語ずつ移植したもの:
// Fragment shader、1 ピクセル 1 invocation、並列
float x = 0.0, y = 0.0;
float x2 = 0.0, y2 = 0.0;
int iter = MAX_ITER;
for (int i = 0; i < MAX_ITER; i++) {
if (x2 + y2 > 4.0) { iter = i; break; }
y = 2.0 * x * y + c.y;
x = x2 - y2 + c.x;
x2 = x * x;
y2 = y * y;
}
1 つだけ罠: GLSL ES 1.00 は for ループの上限を コンパイル時定数 にしか許さない。なので MAX_ITER は shader への #define として焼き付け、UI で max-iter スライダを動かすと shader が再コンパイルされる。これが 2 実装の唯一の構造的差異。
速度差
私のラップトップで 512×512、500 iterations、seahorse valley にズームした条件:
| Render | ms | 相対 |
|---|---|---|
| CPU | ~230 | 1× |
| GPU | ~2.5 | ~90× |
驚くほどではない: CPU は 1 スレッドで 262,144 ピクセル × 最大 500 iterations を直列に。GPU は warp/wavefront 幅 (32 か 64) × ドライバが割り当てた compute unit 数だけ並列にこなす。
書いておくべき副題: WebGL で GPU 時間を honest に測るのは案外厄介。gl.drawArrays は非同期で、戻った瞬間には何も render されていない。gl.finish() はドライバ依存で信頼できない。動く のは draw 直後に gl.readPixels(0, 0, 1, 1, …) を呼ぶこと: 1 ピクセル読み戻すには pipeline が flush される必要があるので、これが事実上の barrier になる。デモは 10 回ループして平均を取っている。表示される値は「render 呼出 → ピクセル読み出し可能」の wall clock。
ここから面白くなる
Mandelbrot 集合は self-similar。永遠にズームでき、各スケールで新しい filament・spiral・ミニ Mandelbrot が現れる。ズームというのは要するにビューポートのスケールパラメータを縮めること、つまりピクセル座標を複素平面にマップする以下の式で:
vec2 c = uCenter + vec2(uv.x * uScale * aspect, uv.y * uScale);
uv.x は [-0.5, 0.5]、uScale はビューの複素平面上の高さ、aspect は非正方 canvas の補正。1000 倍ズームインなら uScale を 1/1000 にする。
各ピクセルの c は uCenter + 0.49 × uScale × aspect のような値。uScale が極端に小さくなると、これは 大きい数 + すごく小さい数 になる。そして c は float。GLSL ES 1.00 の highp float は IEEE 754 単精度: 符号 1 + 指数 8 + 仮数 23 ビット。相対精度は 2⁻²³ ≈ 1.19 × 10⁻⁷。
seahorse valley の中心は (-0.744, 0.132) 付近、なので uCenter + … の絶対精度は 0.75 × 1.19 × 10⁻⁷ ≈ 9 × 10⁻⁸。ピクセル間距離 uScale / canvasHeight がこれを下回った瞬間、隣り合うピクセルは同じ複素数に丸められる。複数 (時に 12 個前後) のご近所ピクセルが同一の反復回数を出し、画像はもうフラクタルではなく ブロックモザイク になる。
これが上のスクリーンショット。uScale = 5 × 10⁻⁶、512 ピクセル canvas で、ピクセル間隔は ~10⁻⁸ — float32 解像度より一桁細かい。CPU は同じ数学を double (仮数 52 bit、相対精度 ~2.2 × 10⁻¹⁶) で回しているので何も気付かない。GPU は JPEG quality 0 を通したような絵を返す。
「精度を上げればいいのでは」
簡単には行かない。GLSL ES 1.00 に double はない。WebGL 2 / GLSL ES 3.00 のコア仕様にもまだない。GL_ARB_gpu_shader_fp64 拡張を使う desktop OpenGL 4.0+ で初めて double が使える。ブラウザを target にするなら float で頑張るしかない。
deep-zoom Mandelbrot コミュニティが実際に使う回避策は 2 系統:
-
double-single / float-float 算術。各座標を 2 つの float の組
(hi, lo)で表現し、和が望む値になるようにする。加算と乗算を Knuth 流の compensated 演算列に展開する。ALU コスト ~2-4 倍、bandwidth ~2 倍を払って、~14 桁の精度を買う。shader はかなり複雑になり、-ffast-math系の最適化を許すと compensated 演算が壊される。 -
Perturbation method。リファレンスピクセル 1 点の orbit を CPU 側
doubleで 1 回計算しておく。他のピクセルは reference からの 小さい delta をfloatで計算する。delta は 0 付近なので、float の相対精度は座標 ~1.0 ではなく距離 ~10⁻¹⁴ に対して効く。glitch detection (reference orbit が近傍と発散する) の数値罠があり、対処 (Zhuoran のアルゴリズム、BLA) は論文 1 本分の労力。本格的な fractal ソフト (Kalles Fraktaler 等) はこれを採用している。
このデモはどちらもやらない。壁を見せる のがデモの目的であって、迂回することではない。
移植可能な学び
Mandelbrot は道具。本題は精度の話: shader 内で uCenter + uv * uScale を計算していて、uCenter と uScale * 何か小さい値 の両方が結果に効くなら、deep zoom でまったく同じバグが待っている。tile rendering、Google Maps 風 mercator 投影、原点から遠いカメラの deferred shading — 「2 つの大きい数の差の下位ビットに有用な情報がある」ような場面では、単精度が同じ形で破綻する。
修正の形もだいたい同じ: 大 + 小ではなく小 + 小 に演算を組み替える。CPU 側 (double) でカメラ/座標の origin を近傍に translate しておき、それ以降の演算をすべて origin 相対 (float) で GPU 側に投げる。これは perturbation Mandelbrot のコアアイデアそのものであり、CAD ソフトが world 座標を double で持って vertex shader には float delta だけ送る理由でもある。
デモが他にやっていること
- CPU/GPU パネルでパン/ズームがロックされているので、ズームしてどこで分岐するかが直接見える
- 6 プリセット が浅い → 深いズームを順に踏ませる。最後が GPU が壊れる地点
- Max-iter スライダ は shader を都度再コンパイル (定数上限ループの制約のため)
-
honest GPU タイマ は
readPixelsflush で 10 フレーム平均
JavaScript ~160 行 + GLSL ~45 行、bundler なし依存なし。
触る
プリセットを左から右にクリック。Spiral までは 2 パネル区別不能。Deep で右に微弱な縦縞が出る。float32 breaks で右はブロックモザイク。そこが仮数 24 bit を踏み越えた瞬間。
🛠 本記事は SEN 合同会社 が公開している小さな開発者ツール群の 1 つ。他は portfolio 一覧 から。