はじめに

マンデルブロ集合の描画を普通に実装していると、浮動小数点表現の限界により意外と簡単に拡大率の上限に達する。例えばFloat64でも倍率は1E14あたりが限界になるが、複雑で面白い画像はもっと深い場所にあったりする。
ここで、摂動論(Perturbation Theory)の考え方を導入すると、この限界を大幅に伸ばすことができる。本記事ではその理論や実装方法、さらに高速化の手法として2023年時点で主流であるBLAも紹介する。

本記事は以下の記事の続編。

実験用コード

摂動論

摂動論についてWikipediaの説明から引用する。

考えている問題Aを、厳密に解ける問題Bに小さな変更（摂動）が加えられた問題であるとみなす。
問題Aの近似解は、問題Bの厳密解に、摂動が加わったことによって生じる小さな補正（摂動項）を加えたものであると考える。
ここで求めるべき摂動項は、問題Bの厳密解の組み合わせ、すなわち一次結合の形で表現出来ると考え、その係数を与えられた条件から順次求める。

本記事で対象にするマンデルブロ集合描画の計算は、複素平面上の注目する座標においてその周辺の区間含めて発散するかどうかを調べることであり、拡大率が大きくなることは周辺区間がどんどん微小になっていくことになる。というわけで摂動論がうまく適用できれば嬉しいのだが、実際に摂動項だけ独立で計算できることが知られている。
以下、数式解説のスクリーンショット。¹

実装

実装する場合は、以下のような手順になる。

参照点に関してのみ、事前にその軌跡を高精度で計算して全て記録
実際の描画に必要な摂動項の部分は、事前計算した軌跡を参照しながら通常の精度(Float64またはFloat32)で計算
参照点の値と摂動項を足し合わせて軌跡を算出し、その発散を判定

参照点の計算は高精度で行う必要があるが、軌跡の記録自体は通常の精度で行なっても問題ない(むしろそうしないと速度的なメリットがほぼなくなってしまう)。
ここで、参照点が発散してしまうと摂動項部の計算が無駄になるので、発散しない点を選ぶ必要がある。ということだったのだが、割と最近(2021年あたり？)Rebasingという手法が開発され、厳密に選ばなくても計算が破綻しないように実装できるようになった。

基本のみ

以下、最小限に実装したつもりのコード。

import numpy as np
from numba import njit, prange
import math

@njit('void(uint8[:,:,:], float64[:,:],  float64, int64, float64)', parallel=True)
def draw_image(image,  ref,  pixel_size,  max_iters, bailout):
    LOG2 = math.log(2)
    height, width = image.shape[0], image.shape[1]
    max_ref_iteration = ref.shape[0]
    white = np.array([255,255,255], np.float32)
    red = np.array([255,0,0], np.float32)
    for y in prange(height):
        dc_y = pixel_size * (y - height//2)
        for x in range(width):
            dc_x = pixel_size * (x - width//2)

            dz_x, dz_y = 0,0
            iteration = 0
            ref_iteration =0
            while iteration < max_iters:
                ref_x, ref_y = ref[ref_iteration][0], ref[ref_iteration][1]
                zx, zy = ref_x+dz_x, ref_y+dz_y
                radius2 = zx*zx + zy*zy
                if radius2 > bailout: break

                # Rebasing
                dradius2 = dz_x*dz_x + dz_y*dz_y
                if radius2<dradius2 or ref_iteration==max_ref_iteration-1:
                    dz_x,dz_y = zx, zy
                    ref_iteration = 0
                    ref_x, ref_y = ref[0][0], ref[0][1]

                dz_x_1 = 2*(ref_x*dz_x-ref_y*dz_y) + (dz_x*dz_x - dz_y*dz_y)  + dc_x
                dz_y_1 = 2*(ref_x*dz_y+ref_y*dz_x) + 2*dz_x*dz_y + dc_y
                dz_x, dz_y = dz_x_1, dz_y_1

                iteration += 1
                ref_iteration+=1

            if iteration<max_iters:
                color = white if iteration<max_ref_iteration else red
                image[height-1-y, x] = color*((iteration/1000)%1)
    return

from decimal import Decimal, getcontext
def calc_references( center_real, center_imag, magnification, max_iters, bailout, dtype):
    getcontext().prec = int(Decimal(magnification).log10())+5
    Z=[]
    radius2 = 0.0
    iteration = 0
    c_x, c_y = Decimal(center_real), Decimal(center_imag)
    z_x, z_y = 0, 0
    while iteration < max_iters and radius2 <= bailout:
        Z.append((z_x,z_y))
        z_x, z_y = (z_x*z_x-z_y*z_y) + c_x, (2*z_x*z_y) + c_y
        radius2 = (z_x*z_x) + (z_y*z_y)
        iteration += 1
    return np.array(Z,dtype)

def calc_pixel_size(width, magnification,dtype ):
    return dtype(Decimal(4.0/width) / Decimal(magnification))


def draw(width, height, center_real, center_imag, magnification, max_iters, bailout):
    ref = calc_references( center_real, center_imag, magnification, max_iters, bailout, np.float64)
    pixel_size = calc_pixel_size(width, magnification, np.float64)

    image = np.zeros((height, width,3), dtype = np.uint8)
    draw_image(image, ref, pixel_size, max_iters, bailout)
    return image

width, height = 512, 512
center_real = "-1.768610493014677074503175653270226520239677907588665494812359766720721863405719772"
center_imag = "0.001266613503868717702066411192242601576193940560471409839484404683951701639387836214"
magnification = '1e72'

max_iters = 10000
bailout=4.0

#draw
image = draw(width, height, center_real,center_imag, magnification, max_iters, bailout)

座標や拡大率に関しては、文字列で渡さないと丸められてしまうので注意。
参照点は画像の中央として、事前にDecimalを使って高精度で計算している。
以下、これを使って描画した画像。

画像中で赤い部分はRebasingが効いている箇所。これは画面中央の参照点より発散が遅い場所で、実際に対応する参照点は計算していないことを示す。
ちなみに普通に着色すると以下のような画像になる。

倍率を1E72としており、前述の1E14という制限を大きく超えている。この実装での制限は1E300程度になる。1E14というのは仮数部のbit数による制限で、一方1E300というのは指数部の制限による。
これ以上を実現するためには、指数部だけ別に保持しながら演算を行う手法がよく使われるようだ。

マンデルブロ集合の描画で"深い場所"というのは発散までの回数が大きい(10万をこえるのがざら)ということなので、CPUのみだと計算時間的に厳しくなってくる。
これについてはCUDA対応で速くすることが考えられるが、実際に対応した場合との描画速度比較は次節でまとめて行う。

BLA

摂動論では摂動部分を冪級数展開して高次の成分を無視するという近似計算で使うことが多いようだが、マンデルブロ集合の描画でも高次成分を無視することで計算を一気にスキップさせるというテクニックが用いられる。最新の描画アプリだと、CPUのみでも場合によっては上記のCUDA対応の実装よりも短時間で描画が終わったりする。

かつては摂動論導入と同時に提案されたSeries Approximationが主流だったが、2023年現在はBivariate Linear Approximation(BLA)が取って代わっているようだ。(これも先ほどのRebasingの提唱者と同じ人が開発した模様)

以下、BLAの数式解説のスクリーンショット。²

上の数式で言うと、mが現在のステップ数として、条件が合えばそこからlステップ分スキップしてnに行けることを示している。

実装面ではこのサイトの説明が比較的わかりやすい。

まず事前に係数と閾値のテーブルを以下のように作成する。

参照点での軌跡を高精度で計算
参照点の各ステップ(iteration)で1ステップ先に行ける係数と、適用条件の閾値を求める
上で計算した係数を２回連続で適用させて、2つ先に行ける係数と、適用条件の閾値を計算
これを繰り返して、4、8、16、・・・と先に行ける係数と条件を限界まで計算

ここでの係数計算は通常の精度で可能で、かつ一ステップ目の計算以降は過去の計算結果を利用できるので、テーブル作成は高速に実行できる。

発散判定の実装は、ここで作成したテーブルを参照しながら適用条件を満たす場合に随時スキップするようにする。

実際のコーディングに関しては、こちらのコードが大変参考になった。
以下、BLAの計算と描画のコード。

from typing import NamedTuple
class BLA(NamedTuple):
    A: float
    B: float
    r: float

@njit('float64[:,:,:](float64,  float64[:,:] )')
def calc_bla_coef(pixel_size, ref):
    def cabs(x): 
        return math.sqrt(x[0]*x[0] + x[1]*x[1])
    def cmul(x,y,coef=1.0): 
        return ( coef*(x[0]*y[0]-x[1]*y[1]), coef*(x[0]*y[1]+x[1]*y[0]) )
    def cadd(x,y): 
        return (x[0]+y[0], x[1]+y[1])
    blas = []
    bla = []
    A0, B0 = (1.0,0.0), (0.0,0.0)
    epsilon = 2**-24 # 1e-6 # 2**-53 # 
    for m in range(1, len(ref)):
        zm = (ref[m,0],ref[m,1])
        A = cmul(zm,A0, 2.0)
        B = cadd( cmul(zm,B0,2.0) , (1.0,0.0))
        xA = cabs(A)
        r = max(0,epsilon*xA - pixel_size/xA)
        b1 = BLA(A, B, r)
        bla.append(b1)
    blas.append(bla)

    src = 0
    msrc = len(ref)-1
    while(msrc>1):
        mdst = (msrc+1)//2
        bla = []
        blas.append(bla)
        for m in range(mdst):
            m2   = m*2
            m2_1 = m2 + 1
            if m2_1<msrc:
                x2   = blas[src][m2]
                x2_1 = blas[src][m2_1]
                Al2 = cmul( x2_1.A, x2.A )
                Bl2 = cadd( cmul( x2_1.A, x2.B ) , x2_1.B)
                xA = cabs(x2.A)
                xB = cabs(x2.B)
                r = min(x2.r, max(0, (x2_1.r - xB*pixel_size)/xA ))
                bla.append( BLA(Al2, Bl2, r) ) 
            else:
                bla.append( blas[src][m2] ) 

        msrc = (msrc+1)//2
        src+=1
    buf = np.zeros((len(blas), len(blas[0]), 5),  np.float64)
    for i,bla in enumerate(blas):
        # print(i,bla[0])
        for j in range(len(bla)):
            buf[i,j,0] = bla[j].A[0]
            buf[i,j,1] = bla[j].A[1]
            buf[i,j,2] = bla[j].B[0]
            buf[i,j,3] = bla[j].B[1]
            buf[i,j,4] = bla[j].r
    return buf

@njit('void(uint8[:,:,:], float64[:,:],  float64[:,:,:], float64, int64, float64, float64)', parallel=True)
def bla_draw_image(image,  ref, blas,  pixel_size,  max_iters, bailout, density ):
    LOG2 = math.log(2)
    height, width = image.shape[0], image.shape[1]
    max_ref_iteration = ref.shape[0]

    for y in prange(height):
        dc_y = pixel_size * (y - height//2)
        for x in range(width):
            dc_x = pixel_size * (x - width//2)
            dz_x, dz_y = 0.0,0.0
            tmp_color = np.zeros((3,), np.float32)

            iteration = 0
            ref_iteration =0
            while iteration < max_iters:
                ref_x, ref_y = ref[ref_iteration][0], ref[ref_iteration][1]
                zx, zy = ref_x+dz_x, ref_y+dz_y
                radius2 = zx*zx + zy*zy
                if radius2 > bailout: break

                # Rebasing
                dradius2 = dz_x*dz_x + dz_y*dz_y
                if radius2<dradius2 or ref_iteration==max_ref_iteration-1:
                    dz_x,dz_y = zx, zy
                    ref_iteration = 0
                    ref_x, ref_y = ref[0][0], ref[0][1]

                # BLA
                m = ref_iteration
                l = 0
                if m!=0:
                    r = math.sqrt(dradius2)
                    ix = m-1
                    for level in range(blas.shape[0]):
                        ixm = ix * (2**level) + 1 
                        if m==ixm and r<blas[level,ix,4]:
                            b = blas[level,ix]
                            l = 2**level
                        else:
                            break
                        ix = ix//2
                n = m+l
                if l!=0 and n < max_ref_iteration:
                    A_x, A_y, B_x, B_y = b[0],b[1],b[2],b[3]
                    dz_x_1 = (dz_x*A_x - dz_y*A_y) + (dc_x*B_x - dc_y*B_y)
                    dz_y_1 = (dz_x*A_y + dz_y*A_x) + (dc_x*B_y + dc_y*B_x)
                    ref_iteration += l
                    iteration += l
                else:
                    dz_x_1 = 2*(ref_x*dz_x-ref_y*dz_y) + (dz_x*dz_x - dz_y*dz_y)  + dc_x
                    dz_y_1 = 2*(ref_x*dz_y+ref_y*dz_x) + 2*dz_x*dz_y + dc_y
                    ref_iteration+=1
                    iteration += 1
                dz_x, dz_y = dz_x_1, dz_y_1


            if iteration<max_iters:
                # smoothing
                log_zn = math.log(radius2) / 2
                nu = math.log(log_zn / LOG2) / LOG2
                alpha = iteration + 1 - nu
                # alpha = iteration

                alpha *= 0.05*density
                tmp_color = get_color(alpha)
            image[height-1-y, x] = tmp_color
    return

描画速度比較

計測は上記のカラー画像を生成するのに要した時間で行う。座標等の設定は以下の通り。

width, height = 512, 512
center_real = "-1.768610493014677074503175653270226520239677907588665494812359766720721863405719772"
center_imag = "0.001266613503868717702066411192242601576193940560471409839484404683951701639387836214"
magnification = '1e72'
sampling=2

事前に行う参照点の計算時間は共通かつ無視できる程度なので、以下の比較には含めない。

PTはPerturbation Theoryの略。CoalbのGPUはT4使用。
ColabのCPUは現在の一般的なPC用のCPUに比べて低性能なので、筆者所有のM1 Mac Miniでの結果も記載する。

環境	手法	時間(秒)
Colab(CPU)	PT	49.916
Colab(CPU)	PT＋BLA	8.055
Colab(GPU)	PT	1.998
Colab(GPU)	PT+BLA	0.303
M1 Mac(CPU)	PT	9.227
M1 Mac(CPU)	PT+BLA	0.776

GPUを使ったPT+BLAが最も速いが、M1でのPT+BLAと大差ない。GPUとしてはT4よりP100の方がFloat64性能が大幅に高いので、そちらと比較するともう少し差が出る。ただしM1より高性能なCPUを使えば、当然もっと速くなるので差が縮まる。

BLAで高速化できる程度は座標等によって異なり、場合によってはもっと速くなることもある。

事前に行う参照点の計算時間は今回の比較で省略したが、高精度で行うこともあり数秒から数分かかる場合もある。ここが一番時間がかかる処理になりがちなので、多くの描画アプリでは結果をキャッシュして再利用できるようにしているようだ。一点にズームしていくような動画を作成する場合は、最初だけ計算すればあとはそれをずっと使い続ければ良い。

まとめ

マンデルブロ集合描画の限界を伸ばし描画速度を上げる、摂動論及びBLAを実装して効果を確認した。
アルゴリズムだけでここまで性能向上するのには驚かされた。

最後に、4x4スーパーサンプリングで生成した画像をサンプルとしていくつか掲載しておく。