最新研究 Raspberry Pi 4のGPUプログラミング

$ git clone https://github.com/Idein/py-videocore6.git

$ sudo apt update
$ sudo apt install python3-pip
$ pip3 install --user --upgrade pip
$ pip3 install --user .

nop(sig=ldunifrf(r0))
L.loop_i
sub(r0, r0, 1, cond=’pushz’)
b(R.loop_i, cond=’na0’)
nop() # 遅延スロット
nop() # 遅延スロット
nop() # 遅延スロット

nop(sig=ldunifrf(r0))
with loop as loop_i:
    sub(r0, r0, 1, cond=’pushz’)
    loop_i.b(cond=’na0’)
    nop() # 遅延スロット
    nop() # 遅延スロット
    nop() # 遅延スロット

void *address = (void*)A;
uint32_t data = 0;
do {
    *(uint32_t)address = data;
    address = address + 4;
    data = data + 1;
    length = length - 1;
} while (length != 0);

length = length / 16;
void *address = (void*)A + 4 * element_index;
uint32_t data = element_index;
do {
    *(uint32_t*)address = data;
    address = (void*)address + 4 * 16;
    data = data + 16;
    length = length - 1;
} while (length != 0);

# np.array_equal (配列の要素の比較)のため
import numpy as np

# py-videocore6 ライブラリ
from videocore6.assembler import qpu
from videocore6.driver import Driver

# QPU で動作するメソッド

@qpu
def qpu_write_seq(asm):

    # 使用するレジスタに名前を付ける
    # 今回はアキュムレータを一時的な領域として，レジスタを値の格納に用いる
    reg_length  = rf0
    reg_address = rf1
    reg_stride  = rf2
    reg_data    = rf3

    # reg_length = length / 16
    nop(sig=ldunifrf(reg_length))
    shr(reg_length, reg_length, 4)

    # reg_address = address + 4 * element_index
    nop(sig=ldunifrf(reg_address))
    eidx(r0)
    shl(r0, r0, 2)
    add(reg_address, reg_address, r0)

    # それぞれのスレッドが 4 バイトを書き込み，1 QPU に 16 スレッドあるので，
    # ループするごとにアドレスに 64 を加える
    # reg_stride = 4 * 16 = 64
    shl(reg_stride, 4, 4)

    # reg_data = element_index
    eidx(reg_data)

    with loop as l:

        # 現在のアドレスにデータを書き込み，完了を待つ

        mov(tmud, reg_data)
        mov(tmua, reg_address)
        tmuwt()

        # アドレスに 64 を加える
        add(reg_address, reg_address, reg_stride)

        # データに 16 を加える
        # 即値は -16～15 しか使えないため，d + 16 = d - (-16)の関係を用いる
        sub(reg_data, reg_data, -16)

        # length から 1 を引き，0 になるか調べる
        sub(reg_length, reg_length, 1, cond='pushz')

        # 0 でない場合，ループの先頭にジャンプする
        l.b(cond='na0')

        # 遅延スロット 3 つ
        nop()
        nop()
        nop()

    # QPU プログラム実行終了のおまじない
    nop(sig=thrsw)
    nop(sig=thrsw)
    nop()
    nop()
    nop(sig=thrsw)
    nop()
    nop()
    nop()


# CPU で動作するメソッド

def write_seq():

    # 配列の長さ
    length = 256

    # py-videocore6 ドライバを初期化する
    with Driver()as drv:

        # QPU のプログラムをアセンブルする
        code = drv.program(qpu_write_seq)

        # 32 ビット整数の長さ length の配列を確保し，中身を 0 で埋める
        A = drv.alloc(length, dtype='uint32')
        A.fill(0)

        # uniform に length と A の先頭アドレスをセットする
        unif = drv.alloc(1024, dtype='uint32')
        unif[0] = length
        unif[1] = A.addresses()[0]

        # QPU プログラム実行前の配列 A の中身を表示する
        # すべて 0 のはず

        print('Before execution:')
        print(A)

        # QPU プログラムを実行する
        drv.execute(code, unif.addresses()[0])

        # QPU プログラム実行後の配列 A の中身を表示する

        # 0?length-1 の値が格納されているはず
        print('After execution:')
        print(A)

        # 本当に 0?length-1 の値が格納されているか確認する
        assert np.array_equal(A, np.arange(length))
        print('Test passed!')

def main():
    write_seq()
if __name__ == '__main__':
    main()

$ python3 write_seq.py

with Driver()as drv:

code = drv.program(qpu_write_seq)

A = drv.alloc(length, dtype='uint32')

A.fill(0)

print(A)

unif = drv.alloc(1024, dtype='uint32')
unif[0] = length
unif[1] = A.addresses()[0]

drv.execute(code, unif.addresses()[0])

@qpu

red_1cpu_9}]
reg_length  = rf0
reg_address = rf1
reg_stride  = rf2
reg_data    = rf3

# reg_length = length / 16
nop(sig=ldunifrf(reg_length))
shr(reg_length, reg_length, 4)

# reg_address = address + 4 * element_index
nop(sig=ldunifrf(reg_address))
eidx(r0)
shl(r0, r0, 2)
add(reg_address, reg_address, r0)

# それぞれのスレッドが 4 バイトを書き込み，1 QPU に 16 スレッドあるので，
# ループするごとにアドレスに 64 を加える
# reg_stride = 4 * 16 = 64
shl(reg_stride, 4, 4)

# reg_data = element_index
eidx(reg_data)

with loop as l:

mov(tmud, reg_data)
mov(tmua, reg_address)
tmuwt()

add(reg_address, reg_address, reg_stride)

sub(reg_data, reg_data, -16)

sub(reg_length, reg_length, 1, cond='pushz')

l.b(cond='na0')

nop()
nop()
nop()

nop(sig=thrsw)
nop(sig=thrsw)
nop()
nop()
nop(sig=thrsw)
nop()
nop()
nop()

# np.array_equal (配列の要素の比較)のため
import numpy as np

# py-videocore6 ライブラリ
from videocore6.assembler import qpu
from videocore6.driver import Driver

# QPU で動作するメソッド
@qpu
def qpu_write_multi(asm):

    reg_qpu_num = rf0

    # QPU 番号 (qpu_num)は (thread_index >> 2)& 0b1111 となる
    tidx(r0)
    shr(r0, r0, 2)
    band(reg_qpu_num, r0, 0b1111)

    # element_index + 16 * reg_qpu_num を書き込む
    eidx(r0)
    shl(r1, reg_qpu_num, 4)
    add(tmud, r0, r1)

    # 書き込み先は A + 4 * element_index + 64 * reg_qpu_num
    nop(sig=ldunifrf(r0))
    eidx(r1)
    shl(r1, r1, 2)
    add(r0, r0, r1)
    shl(r1, reg_qpu_num, 6)
    add(tmua, r0, r1)

    # 書き込みの完了を待つ
    tmuwt()

    # QPU プログラム実行終了のおまじない
    nop(sig=thrsw)
    nop(sig=thrsw)
    nop()
    nop()
    nop(sig=thrsw)
    nop()
    nop()
    nop()

# CPU で動作するメソッド
def write_multi():

    # 8 つの QPU を使用する
    nthreads = 8

    # py-videocore6 ドライバを初期化する
    with Driver()as drv:

        # QPU のプログラムをアセンブルする
        code = drv.program(qpu_write_multi)

        # 32 ビット整数の長さ 16 * nthreads の配列を確保し，
        # 中身を 0 で埋める
        A = drv.alloc((nthreads, 16), dtype='uint32')
        A.fill(0)

        # uniform に A の先頭アドレスをセットする
        unif = drv.alloc(1024, dtype='uint32')
        unif[0] = A.addresses()[0, 0]

        # QPU プログラム実行前の配列 A の中身を表示する

        # すべて 0 のはず
        print('Before execution:')
        print(A)

        # QPU プログラムを実行する
        drv.execute(code, unif.addresses()[0], thread=nthreads)

        # QPU プログラム実行後の配列 A の中身を表示する
        print('After execution:')
        print(A)

        # A に 0?16*nthreads-1 の値が格納されているか確認する
        assert np.array_equal(A.reshape(-1), range(16 * nthreads))
        print('Test passed!')

def main():
    write_multi()
if __name__ == '__main__':
    main()

$ python3 write_multi.py

A = drv.alloc((nthreads, 16), dtype='uint32')

tidx(r0)
shr(r0, r0, 2)
band(reg_qpu_num, r0, 0b1111)

eidx(r0)
shl(r1, reg_qpu_num, 4)
add(tmud, r0, r1)

nop(sig=ldunifrf(r0))
eidx(r1)
shl(r1, r1, 2)
add(r0, r0, r1)
shl(r1, reg_qpu_num, 6)
add(tmua, r0, r1)

$ sudo apt update
$ sudo apt install libstb-dev

void mandelbrot_naive(const int height, const int width, const int maxiter,
        uint32_t *img, const float ymin, const float ymax, const float xmin,
        const float xmax)
{
    const float dy = (ymax - ymin)/ height;
    const float dx = (xmax - xmin)/ width;

    float cy = ymin;
    for (int i = 0; i < height; ++i, cy += dy){

        float cx = xmin;
        for (int j = 0; j < width; ++j, cx += dx){
            float y = 0.f, x = 0.f;
            int k;
            for (k = maxiter - 1; k >= 0; --k){
                const float y2 = 2.f * x * y + cy;
                const float x2 = x * x - y * y + cx;
                if (y2 * y2 + x2 * x2 > 4.f)
                    break;
                y = y2;
                x = x2;
            }
            /*
             * k < 0 つまり maxiter 回の反復を行なっても発散しなかった場合，

             * n → ∞ で発散しない，つまりマンデルブロ集合とみなす．
             * マンデルブロ集合を白 (ALPHA = B = G = R = 255)，

             * そうでないピクセルを黒 (ALPHA = 255, B = G = R = 0)にする．
             */
            *img++ = (k < 0)? 0xffffffff : 0xff000000;
        }
    }
}

#include 
#include 
#include 
#include 
#include 

#define STB_IMAGE_WRITE_IMPLEMENTATION
#include 

static
double getsec(void)
{
    struct timespec t;
    clock_gettime(CLOCK_MONOTONIC, &t);
    return t.tv_sec + t.tv_nsec * 1e-9;
}

int main(void)
{
    const int height = 1024, width = 1024;
    const int maxiter = 16384;
    const float ymin = -1.5f, ymax = 1.5f;
    const float xmin = -2.f, xmax = 1.f;

    uint32_t *img;
    int ret = posix_memalign((void**)&img, 4096,
            height * width * sizeof(*img));
    assert(!ret);

    const double start = getsec();
    mandelbrot_naive(height, width, maxiter, img, ymin, ymax, xmin, xmax);
    const double end = getsec();
    printf("Time: 

    stbi_write_png("mandelbrot.png", width, height, sizeof(*img), img,
            sizeof(*img)* width);

    return 0;
}

const float dy = (ymax - ymin)/ height;
const float dx = (xmax - xmin)/ width;

float cy = ymin;
for (int i = 0; i < height; ++i, cy += dy){
    float cx = xmin;
    for (int j = 0; j < width; ++j, cx += dx){

import math
import time

from PIL import Image
import numpy as np

from videocore6 import pack_unpack
from videocore6.assembler import qpu
from videocore6.driver import Driver

@qpu
def qpu_mandelbrot(asm, *, nqpus, unroll_shift):

    params = [
            'height', 'width',
            'nloops_i', 'nloops_j', 'maxiter',
            'addr',
            'ymin', 'xmin',
            'dy', 'dx',
    ]

    values = [
            'i', 'j', 'k',
            'cx', 'cy',
            'x', 'y',
            'color_black',
            'stride',
    ]

    g = globals()
    for i, reg in enumerate(params + values):
        g[f'reg_{reg}'] = g[f'rf{i}']

    for param in params:
        nop(sig = ldunifrf(g[f'reg_{param}']))

    # reg_addr += 4 * element_index
    eidx(r0)
    shl(r0, r0, 2)
    add(reg_addr, reg_addr, r0)

    # reg_maxiter = reg_maxiter - 1
    shr(reg_maxiter, reg_maxiter, unroll_shift)
    sub(reg_maxiter, reg_maxiter, 1)

    if nqpus == 1:
        mov(reg_stride, 0)
    elif nqpus == 8:
        # stride = 4 * width * 7
        umul24(r0, reg_width, 7)
        shl(reg_stride, r0, 2)
        # r0 = qpu_num
        # reg_addr += 4 * width * qpu_num
        tidx(r0)
        shr(r0, r0, 2)
        band(r0, r0, 0b1111)
        umul24(r1, r0, reg_width)
        shl(r1, r1, 2)
        add(reg_addr, reg_addr, r1)
        # reg_ymin += reg_dy / 8 * qpu_num
        itof(r0, r0)
        fmul(r1, reg_dy, 1. / 8)
        fmul(r0, r1, r0)
        fadd(reg_ymin, reg_ymin, r0)
    else:
        raise Exception('nqpus must be 1 or 8')

    # reg_xmin += reg_dx / 16 * element_index
    eidx(r0)
    itof(r0, r0)
    fmul(r1, reg_dx, 1. / 16)
    fmul(r0, r0, r1)
    fadd(reg_xmin, reg_xmin, r0)

    # reg_color_black = 0xff000000
    mov(r0, 0xffffffff)
    shl(r0, r0, 12)
    shl(reg_color_black, r0, 12)

    # i = nloops_i, nloops_i - 1, ..., 1
    mov(reg_i, reg_nloops_i)

    # cy = ymin, ymin + dy, ..., ymin + dy * (height - 1)
    mov(reg_cy, reg_ymin)
    with loop as ly:

        # j = nloops_j, nloops_j - 1, ..., 1
        # cx = xmin, xmin + dx * 16, ..., xmin + dx * 16 * (width / 16 - 1)
        mov(reg_j, reg_nloops_j).mov(reg_cx, reg_xmin)
        # z = x + i y = 0
        mov(reg_x, 0).mov(reg_y, 0)
        mov(reg_k, reg_maxiter)
        with loop as lx:

            # k = maxiter - 1, maxiter - 2, ..., 0
            with loop as lit:

                for i in range(1 << unroll_shift):
                    # r0 = x * x + cx
                    # r1 = y * y
                    # r2 = 2. * x
                    fadd(r2, reg_x, reg_x).fmul(r0, reg_x, reg_x)
                    fadd(r0, r0, reg_cx).fmul(r1, reg_y, reg_y)
                    # x' = x * x - y * y + cx
                    # r0 = 2. * x * y
                    fsub(reg_x, r0, r1).fmul(r0, r2, reg_y)
                    # y' = 2 * x * y + cy
                    # r1 = x' * x'
                    fadd(reg_y, r0, reg_cy).fmul(r1, reg_x, reg_x)
                # r0 = 4 - x' * x' - y' * y'
                # k = k - 1
                fsub(r0, 4., r1).fmul(r1, reg_y, reg_y)
                fsub(r0, r0, r1).sub(reg_k, reg_k, 1)

                # If either of r0 and reg_k is negative, then exit the loop.
                bor(null, r0, reg_k, cond='pushn')
                lit.b(cond='anyna')
                mov(null, r0, cond='pushn')
                mov(r0, reg_color_black, cond='ifa').mov(r0, 0xffffffff, cond='ifna')
                mov(r1, reg_addr)

            sub(reg_j, reg_j, 1, cond='pushz').mov(tmud, r0)
            fadd(reg_cx, reg_cx, reg_dx).mov(tmua, r1)
            lx.b(cond='anyna')
            mov(reg_x, 0).mov(reg_y, 0)
            shl(r0, 4, 4).mov(reg_k, reg_maxiter)
            tmuwt().add(reg_addr, reg_addr, r0)

        sub(reg_i, reg_i, 1, cond = 'pushz')
        ly.b(cond = 'anyna')
        fadd(reg_cy, reg_cy, reg_dy)
        add(reg_addr, reg_addr, reg_stride)
        nop()

    nop(sig = thrsw)
    nop(sig = thrsw)
    nop()
    nop()
    nop(sig = thrsw)
    nop()
    nop()
    nop()

def mandelbrot(*, height, width, maxiter, ymin, ymax, xmin, xmax, ydiv, xdiv,
        unroll_shift, nqpus):

    assert height 
    assert width 
    nloops_i = height // ydiv
    nloops_j = width // xdiv
    dy = (ymax - ymin)/ height * ydiv
    dx = (xmax - xmin)/ width * xdiv

    assert maxiter 

    with Driver()as drv:

        code = drv.program(qpu_mandelbrot, nqpus=nqpus,
                unroll_shift=unroll_shift)
        #drv.dump_code(code); exit()

        A = drv.alloc((height, width), dtype = 'uint32')

        unif = drv.alloc(1024, dtype = 'uint32')

        for i, v in enumerate([height, width, nloops_i, nloops_j, maxiter,
                A.addresses()[0, 0],
                *pack_unpack('f', 'I', [ymin, xmin, dy, dx])]):
            unif[i] = v

        start = time.time()
        drv.execute(code, unif.addresses()[0], thread=nqpus)
        end = time.time()

        print(f'Time: {end - start} sec')

        Image.fromarray(A, mode = 'RGBA').save('mandel.png')


def main():

    nqpus = 8

    unroll_shift = 6

    maxiter = 16384
    height, width = 1024, 1024
    xmin, xmax = -2., 1.
    ymin, ymax = -1.5, 1.5

    print(f'Running with {nqpus} QPUs')
    print(f'unroll = {1 << unroll_shift}')
    print(f'maxiter = {maxiter}')
    print(f'height = {height}, width = {width}')
    print(f'xmin = {xmin}, xmax = {xmax}')
    print(f'ymin = {ymin}, ymax = {ymax}')

    mandelbrot(height=height, width=width, maxiter=maxiter, ymin=ymin,
            ymax=ymax, xmin=xmin, xmax=xmax, ydiv=nqpus, xdiv=16,
            unroll_shift=unroll_shift, nqpus=nqpus)

if __name__ == '__main__':
    main()

#include 
#include 
void mandelbrot_neon_naive(const int height, const int width, const int maxiter,
        uint32_t *img, const float ymin, const float ymax, const float xmin,
        const float xmax)
{

#pragma omp parallel default(none)\
        firstprivate(height, width, img, ymin, ymax, xmin, xmax)
    {

        const int nthreads = omp_get_num_threads();
        const int thread_num = omp_get_thread_num();
        assert(height 

        printf("

        const int nlanes = 4;
        assert(width 

        const float32x4_t dy = vdupq_n_f32((ymax - ymin)/ height * nthreads);
        const float32x4_t dx = vdupq_n_f32((xmax - xmin)/ width * nlanes);

        float32x4_t cx_init = vdupq_n_f32(0.f);
        for (int lane = 0; lane < nlanes; ++lane)

            cx_init = vsetq_lane_f32((xmax - xmin)/ width * lane, cx_init,
                    lane);

        cx_init = vaddq_f32(cx_init, vdupq_n_f32(xmin));

        img += width * thread_num;

        float32x4_t cy = vdupq_n_f32(
                ymin + (ymax - ymin)/ height * thread_num);
        for (int i = 0; i < height / nthreads; ++i){
            float32x4_t cx = cx_init;
            for (int j = 0; j < width / nlanes; ++j){
                uint32x4_t is_working;
                float32x4_t y = vdupq_n_f32(0.f), x = vdupq_n_f32(0.f);

                for (int k = maxiter - 1; k >= 0; --k){
                    float32x4_t s = vmlaq_f32(cx, x, x);
                    float32x4_t t = vaddq_f32(x, x);
                    x = vmlsq_f32(s, y, y);
                    y = vmlaq_f32(cy, t, y);

                    float32x4_t u = vmlaq_f32(vmulq_f32(y, y), x, x);
                    /* is_working = (x * x + y * y <= 4.f)? 0xffffffff : 0 */
                    is_working = vcleq_f32(u, vdupq_n_f32(4.f));

                    uint32x2_t v = vorr_u32(vget_low_u32(is_working),
                            vget_high_u32(is_working));
                    if (!vget_lane_u64(vreinterpret_u64_u32(v), 0))
                        break;
                }

                /*
                 * is_working は，各ピクセルにおいて発散した場合は 0，

                 * 発散しなかった場合は 0xffffffff になっている．
                 * 発散したピクセルは 0xff000000 (ALPHA = 255, B = G = R = 0)，
                 * 発散しなかったピクセルは 0xffffffff (ALPHA = B = G = R = 255)
                 * にしたいので，is_working と 0xff000000 のビットごとの論理和を

                 * 取ったものを採用すればよい．
                 */

                vst1q_u32(img, vorrq_u32(is_working, vdupq_n_u32(0xff000000)));
                cx = vaddq_f32(cx, dx);
                img += nlanes;
            }
            cy = vaddq_f32(cy, dy);
            img += width * (nthreads - 1);
        }
    }
}

void mandelbrot_neon_unroll4(const int height, const int width, const int maxiter,
        uint32_t *img, const float ymin, const float ymax, const float xmin,
        const float xmax)
{

#pragma omp parallel default(none)\
        firstprivate(height, width, img, ymin, ymax, xmin, xmax)
    {
        const int nthreads = omp_get_num_threads();
        const int thread_num = omp_get_thread_num();
        assert(height 

        printf("

        /*
         * ARM の NEON は 128 ビット幅のレジスタを扱えるため，
         * 32 ビットの整数および浮動小数点数を同時に 4 つ処理することができる
         */

        const int nlanes = 4;
        assert(width 

        const float32x4_t dy = vdupq_n_f32((ymax - ymin)/ height * nthreads);
        const float32x4_t dx = vdupq_n_f32((xmax - xmin)/ width * nlanes * 4);

        float32x4_t cx0_init = vdupq_n_f32(0.f);
        for (int lane = 0; lane < nlanes; ++lane)
            cx0_init = vsetq_lane_f32(
                    (xmax - xmin)/ width * (nlanes * 0 + lane)+ xmin,
                    cx0_init, lane);

        float32x4_t cx1_init = vdupq_n_f32(0.f);
        for (int lane = 0; lane < nlanes; ++lane)
            cx1_init = vsetq_lane_f32(
                    (xmax - xmin)/ width * (nlanes * 1 + lane)+ xmin,
                    cx1_init, lane);

        float32x4_t cx2_init = vdupq_n_f32(0.f);
        for (int lane = 0; lane < nlanes; ++lane)
            cx2_init = vsetq_lane_f32(
                    (xmax - xmin)/ width * (nlanes * 2 + lane)+ xmin,
                    cx2_init, lane);

        float32x4_t cx3_init = vdupq_n_f32(0.f);
        for (int lane = 0; lane < nlanes; ++lane)
            cx3_init = vsetq_lane_f32(
                    (xmax - xmin)/ width * (nlanes * 3 + lane)+ xmin,
                    cx3_init, lane);

        img += width * thread_num;

        float32x4_t cy = vdupq_n_f32(
                ymin + (ymax - ymin)/ height * thread_num);
        for (int i = 0; i < height / nthreads; ++i){
            float32x4_t cx0 = cx0_init;
            float32x4_t cx1 = cx1_init;
            float32x4_t cx2 = cx2_init;
            float32x4_t cx3 = cx3_init;
            for (int j = 0; j < width / nlanes / 4; ++j){
                uint32x4_t is_working0, is_working1, is_working2, is_working3;
                float32x4_t y0 = vdupq_n_f32(0.f), x0 = vdupq_n_f32(0.f);
                float32x4_t y1 = vdupq_n_f32(0.f), x1 = vdupq_n_f32(0.f);
                float32x4_t y2 = vdupq_n_f32(0.f), x2 = vdupq_n_f32(0.f);
                float32x4_t y3 = vdupq_n_f32(0.f), x3 = vdupq_n_f32(0.f);
                for (int k = (maxiter >> 4)- 1; k >= 0; --k){

#pragma GCC unroll 65534
                    for (int l = (1 << 4)- 1; l >= 0; --l){
                        float32x4_t s0 = vmlaq_f32(cx0, x0, x0);
                        float32x4_t s1 = vmlaq_f32(cx1, x1, x1);
                        float32x4_t s2 = vmlaq_f32(cx2, x2, x2);
                        float32x4_t s3 = vmlaq_f32(cx3, x3, x3);
                        float32x4_t t0 = vaddq_f32(x0, x0);
                        float32x4_t t1 = vaddq_f32(x1, x1);
                        float32x4_t t2 = vaddq_f32(x2, x2);
                        float32x4_t t3 = vaddq_f32(x3, x3);
                        x0 = vmlsq_f32(s0, y0, y0);
                        x1 = vmlsq_f32(s1, y1, y1);
                        x2 = vmlsq_f32(s2, y2, y2);
                        x3 = vmlsq_f32(s3, y3, y3);
                        y0 = vmlaq_f32(cy, t0, y0);
                        y1 = vmlaq_f32(cy, t1, y1);
                        y2 = vmlaq_f32(cy, t2, y2);
                        y3 = vmlaq_f32(cy, t3, y3);
                    }
                    float32x4_t s0 = vmlaq_f32(vmulq_f32(y0, y0), x0, x0);
                    float32x4_t s1 = vmlaq_f32(vmulq_f32(y1, y1), x1, x1);
                    float32x4_t s2 = vmlaq_f32(vmulq_f32(y2, y2), x2, x2);
                    float32x4_t s3 = vmlaq_f32(vmulq_f32(y3, y3), x3, x3);
                    /* is_working = (x * x + y * y <= 4.f)? 0xffffffff : 0 */
                    is_working0 = vcleq_f32(s0, vdupq_n_f32(4.f));
                    is_working1 = vcleq_f32(s1, vdupq_n_f32(4.f));
                    is_working2 = vcleq_f32(s2, vdupq_n_f32(4.f));
                    is_working3 = vcleq_f32(s3, vdupq_n_f32(4.f));
                    uint32x4_t is_working = vorrq_u32(
                            vorrq_u32(
                                vorrq_u32(is_working0, is_working1),
                                is_working2),
                            is_working3);
                    uint32x2_t t = vorr_u32(vget_low_u32(is_working),
                            vget_high_u32(is_working));
                    if (!vget_lane_u64(vreinterpret_u64_u32(t), 0))
                        break;
                }
                /*
                 * is_working は，各ピクセルにおいて発散した場合は 0，

                 * 発散しなかった場合は 0xffffffff になっている．
                 * 発散したピクセルは 0xff000000 (ALPHA = 255, B = G = R = 0)，
                 * 発散しなかったピクセルは 0xffffffff (ALPHA = B = G = R = 255)
                 * にしたいので，is_working と 0xff000000 のビットごとの論理和を

                 * 取ったものを採用すればよい．
                 */
                vst1q_u32(img + nlanes * 0,
                        vorrq_u32(is_working0, vdupq_n_u32(0xff000000)));
                vst1q_u32(img + nlanes * 1,
                        vorrq_u32(is_working1, vdupq_n_u32(0xff000000)));
                vst1q_u32(img + nlanes * 2,
                        vorrq_u32(is_working2, vdupq_n_u32(0xff000000)));
                vst1q_u32(img + nlanes * 3,
                        vorrq_u32(is_working3, vdupq_n_u32(0xff000000)));
                img += nlanes * 4;
                cx0 = vaddq_f32(cx0, dx);
                cx1 = vaddq_f32(cx1, dx);
                cx2 = vaddq_f32(cx2, dx);
                cx3 = vaddq_f32(cx3, dx);
            }
            img += width * (nthreads - 1);
            cy = vaddq_f32(cy, dy);
        }
    }
}