FPGA上で動くBrainfuckCPUを作った

これはBrainf*ck Advent Calendar 2019 - Adventar25日目の記事です。

Brainfuck命令セットを実行するCPUをVerilog HDLで記述し、FPGAで動かしてみた記録になります。

まず、ハードウェアを意識しつつ命令セットシミュレータ(というかインタプリタ)を作成しました。 ハードウェアを意識しつつ、というのは、インタプリタのメインループの中には一サイクルで完結しそうなものしか書かない、くらいの意味です。 Brainfuckの場合、「対応する[に飛ぶ」のような命令がありますが、これは一サイクルでは終わらないでしょうから、そういったものをどう処理するかを考えながら書くということになります。

以下にC++で書かれた、ハードウェアを意識したBrainfuck命令セットシミュレータを載せます。 実装の都合上、Brainfuck,命令はサポートしていません。

#include <iostream>
#include <array>

struct State {
    size_t ip;
    size_t dp;
    int depth;
};
std::array<unsigned char, 30000> ram;
const char* rom =
#include "mandelbrot.bf"
;

State update( const State& state ) {
    State new_state;
    if( state.depth == 0 ) {
        if( rom[state.ip] == 0x5B && ram[state.dp] == 0 ) { new_state.depth = 1; }
        else if( rom[state.ip] == 0x5D && ram[state.dp] != 0 ) { new_state.depth = -1; }
        else { new_state.depth = 0; }
    } else {
        new_state.depth = state.depth + (rom[state.ip] == 0x5B) - (rom[state.ip] == 0x5D);
    }

    new_state.ip = new_state.depth < 0 ? state.ip - 1 : state.ip + 1;

    new_state.dp = rom[state.ip] == 0x3E && state.depth == 0 ? state.dp + 1 :
                   rom[state.ip] == 0x3C && state.depth == 0 ? state.dp - 1
                                                             : state.dp;

    const bool we = (rom[state.ip] == 0x2B || rom[state.ip] == 0x2D) && state.depth == 0;
    if (we) ram[state.dp] = rom[state.ip] == 0x2D ? ram[state.dp] - 1 : ram[state.dp] + 1;

    if( rom[state.ip] == 0x2E && state.depth == 0 ) { putchar(ram[state.dp]); }

    return new_state;
}

int main() {
    for( State state { 0, 0, 0 }; rom[state.ip]; ) {
            state = update( state );
    }
}

これをほとんどそのままVerilog HDLの記述にすると、Brainfuck CPUの完成です。

module cpu(uart_wr, uart_data, clk, rstn);
    output wire uart_wr;
    output wire[7:0] uart_data;
    input  wire clk;
    input  wire rstn;

// state
    reg[14:0] ip;
    reg[14:0] dp;
    reg[14:0] depth;
    reg sync_rstn;

    wire[14:0] next_ip;
    wire[14:0] next_dp;
    wire[14:0] next_depth;

// temporary
    wire executing;

// module
    wire[7:0] insn;
    wire[7:0] rdata;
    wire we;
    wire[7:0] wdata;
    insnmem insnmem0(next_ip, insn, clk);
    datamem datamem0(next_dp, rdata, we, wdata, clk);

// logic
    assign executing = depth == 0;

    assign next_depth = !sync_rstn ? 0 :
                        executing ? insn == 8'h5B && rdata == 0 ?  1 :
                                    insn == 8'h5D && rdata != 0 ? -1
                                                                :  0
                                  : depth + (insn == 8'h5B) - (insn == 8'h5D);

    assign next_ip = !sync_rstn ? 0 :
                     next_depth[14] ? ip - 1 : ip + 1;

    assign next_dp = !sync_rstn ? 0 :
                     insn == 8'h3E && executing ? dp + 1 :
                     insn == 8'h3C && executing ? dp - 1
                                                : dp;

    assign we = (insn == 8'h2B || insn == 8'h2D) && executing;
    assign wdata = insn == 8'h2D ? rdata - 1 : rdata + 1;

    assign uart_wr = insn == 8'h2E && executing;
    assign uart_data = rdata;

// flip-flop
    always @(posedge clk) begin
        ip <= next_ip;
        dp <= next_dp;
        depth <= next_depth;
        sync_rstn <= rstn;
    end
endmodule

module insnmem(next_ip, insn, clk);
    input  wire[14:0] next_ip;
    output reg  [7:0] insn;
    input  wire       clk;

    reg[7:0] mem[0:32767];
    initial $readmemb("D:/vivado/brainfuck/project_1/mandelbrot.bin", mem);

    always @(posedge clk) begin
       insn <= mem[next_ip];
    end
endmodule

module datamem(addr, rdata, we, wdata, clk);
    input  wire[14:0] addr;
    output reg  [7:0] rdata;
    input  wire       we;
    input  wire [7:0] wdata;
    input  wire       clk;

    reg[7:0] mem[0:32767];
    reg[31:0] i; initial for(i=0; i<32768;i=i+1) mem[i] = 0;

    always @(posedge clk) begin
        if (we) begin
            rdata <= wdata;
            mem[addr] <= wdata;
        end else begin
            rdata <= mem[addr];
        end
    end
endmodule

const教みたいな記述方法になっているので普通のVerilog HDLコードを見慣れている方には奇異なコードに映るかもしれません。

以下、このVerilog HDLコードを解説していきます。

まず、このCPUから出ていく信号線には、uart_wruart_dataがあります。これは、外にUART通信を行うモジュールを作る予定で、そこに送るデータになっています。 このCPUに入ってくる信号線はclkrstnです。それぞれ、クロック信号、負論理リセット信号(1になっている状態がリセット、0になっている状態が通常状態という意味で、負論理を意味するnが末尾についています)になっています。

reg変数はこのCPUのステートです。ipdpdepthは、もとのC++ソースコードにあったものと同じ役割です。 その他、同期化された負論理リセット信号であるsync_rstn用のフリップフロップを加えました。

cpuモジュールの残りの部分はC++コードとほぼ同一なので説明は不要でしょう。 リセットがかかっている間はipdpdepth0にすること、クロックに同期してCPUのステートを(一度に)更新すること、などがHDL特有の点でしょうか。

insnmemモジュールは、命令メモリモジュールです。 クロック同期でデータが読み出せるごく普通のROMになっています。 実行すべきプログラムはCPUを論理合成するときに読み込みます。 一つのプログラムしか動かせないCPUが出来上がるということです。

datamemモジュールは、データメモリモジュールです。 データの書き込みも読み出しもクロック同期のRAMになっています。 rdata <= wdataのような一見変に見える記述は、書き込んだデータを即読み出せるようにするためです。

特筆すべきは、アドレス線が一本しかないことです。 通常、読み書きできるRAMには、読み出し用のアドレス線と書き込み用のアドレス線の二つがあるはずです。 Brainfuckの命令セットには「読んで更新して書き込む (read-modify-write)」命令が存在し、複雑に見えますが、アドレス線自体は一本で十分なのです。 また、命令を読み出さなくてもどの番地を読み出せばよいかが決まるため、何も工夫しなくてもシングルサイクルプロセッサ(一クロックで一命令を処理するプロセッサ)が作れます。

次に、後回しにしていたUART通信モジュールを作りました。

module top(sysclk, cpu_resetn, uart_rx_out);
    input  wire sysclk;
    input  wire cpu_resetn;
    output wire uart_rx_out;

    wire uart_wr;
    wire[7:0] uart_data;

    cpu cpu0(uart_wr, uart_data, sysclk, cpu_resetn);
    uart uart0(uart_rx_out, uart_wr, uart_data, sysclk, cpu_resetn);
endmodule

module uart(uart_tx, uart_wr, uart_data, clk, rstn);
    output reg  uart_tx;
    input  wire uart_wr;
    input  wire[7:0] uart_data;
    input  wire clk;
    input  wire rstn;

    reg[28:0] div_clk;
    reg[10:0] shifter;
    always @(posedge clk or negedge rstn) begin
        if (!rstn) begin
            uart_tx <= 1;
            shifter <= 0;
            div_clk <= 0;
        end else begin
            if (shifter != 0 && !div_clk[28]) begin
                uart_tx <= shifter[0];
                shifter <= shifter >> 1;
            end else if (shifter == 0 && uart_wr) begin
                shifter <= { 2'b11, uart_data[7:0], 1'b0 };
            end
            div_clk <= div_clk[28] ? div_clk + 115200 : div_clk + 115200 - 100000000;
        end
    end
endmodule

UART通信のボーレートは115200Hzを使用することにします。 このCPUは後述するように100MHzで動くので、そのクロックを分周することで115200Hzのクロックを作り出します。 UARTのプロトコルに従うと、送るべきデータは、順に0、8bitデータの下位から上位に向かって1bitずつ、1、1、の11bitです。


さて、このVerilog HDLコードを論理合成し、FPGAに書き込みます。 論理合成~配置配線はVivado 2019.2で行いました。 FPGAは、Xilinx社Artix-7シリーズのNexys Videoというものを使いました。 Vivadoのタイミング解析によれば、合成結果はちょうど100MHzで動作するらしいです。

FPGAmandelbrot.bf*1が内蔵されたCPUを書き込み、実際に動かし、UARTで送られてくる信号をTeraTerm*2で見てみると、確かに動作していることがわかります。


FPGAでCPUを動かせばさぞ速いだろうと思っていたのですが、さすがにBrainfuckくらい簡単な命令セットだと、数GHzで動くCPU上でシミュレーションしたほうが高速です。 なんとかもっと速いBrainfuckCPUを作れないでしょうか?

ここで普通に考えるのは、パイプライン化されたCPUを作るというものです。 上記HDLコードのクリティカルパス(動作周波数を決定する、最も遅延の長いパス)は、命令メモリ読み出し→デコード→next_depth計算(≒分岐判定)→次サイクルIP決定→IP配送、になっています。

しかし、Brainfuckコードには大量の分岐命令が存在するため、安易にパイプライン化をするとストール(パイプラインに仕事を供給できないサイクル)が頻発し、性能が上がらないという事態になります。 この問題に対処するために分岐予測を取り入れた場合、クリティカルパスは、分岐予測器メモリ読み出し→何らかの比較→予測された次サイクルIP決定→IP配送、となりますが、このパスも結構時間がかかり、それほどクロック周波数は上がらないのです。 よって、パイプライン化による性能向上は限られています。

それよりも、BrainfuckCPUの実行効率を落としている、「分岐の飛び先が命令の情報だけでわからない」という部分を解決するハードウェアを追加したほうが、高速化に貢献します。 この目的には、通常のCPUにも搭載されている、Branch Target Buffer (BTB)とまったく同じものが使えます。 なお、分岐方向は予測によるのではなく、実際の実行によって得られる方向を採用します。

BTBを実装したBrainfuckCPUのコードが以下になります。

module top(sysclk, cpu_resetn, uart_rx_out);
    input  wire sysclk;
    input  wire cpu_resetn;
    output wire uart_rx_out;

    wire uart_busy;
    wire uart_wr;
    wire[7:0] uart_data;

    cpu cpu0(uart_busy, uart_wr, uart_data, sysclk, cpu_resetn);
    uart uart0(uart_rx_out, uart_busy, uart_wr, uart_data, sysclk, cpu_resetn);
endmodule

module uart(uart_tx, busy, uart_wr, uart_data, clk, rstn);
    output reg  uart_tx;
    output wire busy;
    input  wire uart_wr;
    input  wire[7:0] uart_data;
    input  wire clk;
    input  wire rstn;

    reg[28:0] div_clk;
    reg[10:0] shifter;

    assign busy = shifter != 0;

    always @(posedge clk or negedge rstn) begin
        if (!rstn) begin
            uart_tx <= 1;
            shifter <= 0;
            div_clk <= 0;
        end else begin
            if (busy && !div_clk[28]) begin
                uart_tx <= shifter[0];
                shifter <= shifter >> 1;
            end else if (!busy && uart_wr) begin
                shifter <= { 2'b11, uart_data[7:0], 1'b0 };
            end
            div_clk <= div_clk[28] ? div_clk + 115200 : div_clk + 115200 - 100000000;
        end
    end
endmodule

module cpu(uart_busy, uart_wr, uart_data, clk, rstn);
    input  wire uart_busy;
    output wire uart_wr;
    output wire[7:0] uart_data;
    input  wire clk;
    input  wire rstn;

// state
    reg[14:0] ip;
    reg[14:0] dp;
    reg[14:0] depth;
    reg[14:0] jump_ip;
    reg sync_rstn;

    wire[14:0] next_ip;
    wire[14:0] next_dp;
    wire[14:0] next_depth;
    wire[14:0] next_jump_ip;

// temporary
    wire executing;
    wire taken;
    wire[14:0] next_depth_temp;

// module
    wire[7:0] insn;
    wire[7:0] rdata;
    wire we;
    wire[7:0] wdata;
    wire btb_hit;
    wire[14:0] btb_target;
    wire btb_we;
    wire[14:0] btb_write_ip;
    wire[14:0] btb_write_target;
    insnmem insnmem0(next_ip, insn, clk);
    datamem datamem0(next_dp, rdata, we, wdata, clk);
    btb btb0(next_ip, btb_hit, btb_target, btb_we, btb_write_ip, btb_write_target, clk);

// logic
    assign executing = depth == 0;

    assign next_depth_temp = executing ? insn == 8'h5B & rdata == 0 ?  1 :
                                         insn == 8'h5D & rdata != 0 ? -1
                                                                    :  0
                                       : depth + (insn == 8'h5B) - (insn == 8'h5D);

    assign taken = executing && next_depth_temp != 0;

    assign next_depth = !sync_rstn ? 0 :
                        taken && btb_hit ? 0 : next_depth_temp;

    assign next_ip = !sync_rstn ? 0 :
                     taken && btb_hit ? btb_target :
                     uart_busy && uart_wr ? ip :
                     next_depth_temp[14] ? ip - 1 : ip + 1;

    assign next_dp = !sync_rstn ? 0 :
                     insn == 8'h3E && executing ? dp + 1 :
                     insn == 8'h3C && executing ? dp - 1
                                                : dp;

    assign next_jump_ip = !sync_rstn ? 0 :
                          taken ? ip : jump_ip;

    assign we = (insn == 8'h2B || insn == 8'h2D) && executing;
    assign wdata = insn == 8'h2D ? rdata - 1 : rdata + 1;

    assign btb_we = next_depth_temp == 0 && !executing;
    assign btb_write_ip = jump_ip;
    assign btb_write_target = next_ip;

    assign uart_wr = insn == 8'h2E && executing;
    assign uart_data = rdata;

// flip-flop
    always @(posedge clk) begin
        ip <= next_ip;
        dp <= next_dp;
        depth <= next_depth;
        jump_ip <= next_jump_ip;
        sync_rstn <= rstn;
    end
endmodule

module btb(next_ip, hit, target, we, write_ip, write_target, clk);
    input  wire[14:0] next_ip;
    output wire       hit;
    output wire[14:0] target;
    input  wire       we;
    input  wire[14:0] write_ip;
    input  wire[14:0] write_target;
    input  wire       clk;

    reg before_we;
    reg[14:0] ip;

    wire[9:0] partial_addr;
    wire[20:0] rdata;
    wire[20:0] wdata;
    btbmem btbmem0(partial_addr, rdata, we, wdata, clk);

    assign hit    = !before_we && rdata[20] && rdata[19:15] == ip[14:10];
    assign target = rdata[14:0];

    assign partial_addr = we ? write_ip[9:0] : next_ip[9:0];
    assign wdata  = { 1'b1, write_ip[14:10], write_target };

    always @(posedge clk) begin
        before_we <= we;
        ip <= next_ip;
    end
endmodule

module btbmem(addr, rdata, we, wdata, clk);
    input  wire [9:0] addr;
    output reg [20:0] rdata;
    input  wire       we;
    input  wire[20:0] wdata;
    input  wire       clk;

    reg[20:0] mem[0:1023];
    integer i; initial for(i=0; i<1024;i=i+1) mem[i] = 0;

    always @(posedge clk) begin
        if (we) begin
            rdata <= wdata;
            mem[addr] <= wdata;
        end else begin
            rdata <= mem[addr];
        end
    end
endmodule

module insnmem(next_ip, insn, clk);
    input  wire[14:0] next_ip;
    output reg  [7:0] insn;
    input  wire       clk;

    reg[7:0] mem[0:32767];
    initial $readmemb("D:/vivado/brainfuck/project_1/mandelbrot.bin", mem);

    always @(posedge clk) begin
       insn <= mem[next_ip];
    end
endmodule

module datamem(addr, rdata, we, wdata, clk);
    input  wire[14:0] addr;
    output reg  [7:0] rdata;
    input  wire       we;
    input  wire [7:0] wdata;
    input  wire       clk;

    reg[7:0] mem[0:32767];
    integer i; initial for(i=0; i<32768;i=i+1) mem[i] = 0;

    always @(posedge clk) begin
        if (we) begin
            rdata <= wdata;
            mem[addr] <= wdata;
        end else begin
            rdata <= mem[addr];
        end
    end
endmodule

BTBは1024エントリ、ダイレクトマップのものを作りました。エントリ数は、件のFPGAに乗っているBrockRAM一つ(36Kbit)に収まるように決定しました。 ポート数をケチるため*3、書き込みポートと読み出しポートを共有化し、BTBに書き込むサイクルの次サイクルではBTBからの読み出しができない(常にBTBミスとなる)ものとしました。BTBへの書き込みは頻繁には起こらない(BTBミスしたときのみBTBに書き込む)ため、これはそれほど問題になりません。 BTBのエントリは、上位ビットから、validビット(1bit)、タグ(5bit)、分岐先(15bit)、の計21bitです。

validビットが0であるか、タグが一致しなかった場合は、BTBミスとなり、元のCPUとまったく同じ動作になります。 BTBヒットであっても、直ちにその値が次のサイクルのIPになるわけではありません。データメモリを読み出した値によっては分岐不成立となり、IP+1が次のサイクルのIPになります。 BTBヒットかつ分岐成立の場合のみ、記録されている分岐先が次のサイクルのIPとなります(taken && btb_hit ? btb_target :)。

BTBを実装しても、動作周波数は変わらず、100MHzになりました。 なお、Brainfuck命令の実行が速くなりすぎ、UART送信にかかる時間よりも短い時間間隔で二つの.命令が実行されてしまうことがあり、送信データに欠けが発生しました。 そのため、UART送信中に.命令が実行されそうになった場合はCPU全体を止めるというuart_busy && uart_wr ? ip :を追加しました。


実際にFPGAに書き込んで動作させてみると、結構早くなった印象を受けました。

C++で書かれたシミュレータで確認したところ、最初のBrainfuckCPUは10521107970命令を31892362997サイクルで実行していました。 BTBを追加したBrainfuckCPUでは、(UART送信に伴うストールを除けば)10521107970命令を10824019551サイクルで実行できるようになりました。 約三倍の高速化です!

まとめ

FPGA上で動く、BTBつきのシングルサイクルBrainfuckCPU(100MHz動作)を作った話でした。

今後の課題

UART受信回路を載せて、プログラムを後から書き込めたり、,命令を実行できたりするCPUを作ったりしていきたいです。

*1:リセットがかかってもデータメモリが0クリアされない問題があるため、元のプログラムの先頭にデータメモリの0番地~307番地(mandelbrot.bfで使用する範囲)を0クリアするコードを追加したものになっています。

*2:端末の改行コード設定をAUTOにしました

*3:クリティカルパスは、命令メモリ読み出し→デコード→next_depth_temp計算→次サイクルIP決定→BTBへの書き込み、となっているため、このポート数削減による微妙な遅延短縮が動作周波数に直結します

VivadoのBlockRAM推論がおかしい

東大のプロセッサ実験のwiki *1を見ると、読み出しアドレスをクロックに同期してregに書き込めば、(write-firstの)BlockRAMに推論されると書かれています。

そこで、1read/1write・読み書き両方ともクロック同期のRAMを、以下のように書いてみます。

module ram(raddr, rdata, we, waddr, wdata, clk);
    input  wire[9:0] raddr;
    output wire[7:0] rdata;
    input  wire      we;
    input  wire[9:0] waddr;
    input  wire[7:0] wdata;
    input  wire      clk;

    reg[7:0] mem[0:1023];
    reg[31:0] i; initial for(i=0; i<1024;i=i+1) mem[i] = 0;

    reg[9:0] raddr_reg;
    always @(posedge clk) begin
        if (we) mem[waddr] <= wdata;
        raddr_reg <= raddr;
    end
    assign rdata = mem[raddr_reg];
endmodule

これをVivadoで合成してみると、確かにBlockRAMに推論されます。

しかし、verilogシミュレーションとタイミングシミュレーションの結果(≒FPGAに焼いた時の結果)がまるで違います!いったいどういうことなのでしょうか。

観察

まず、このRAMを使うモジュールを定義します。書き込むアドレスは0番地だけですが、愚直にそう書いてしまうとうまくいきません。 最適化で他の番地が消えてしまうからです。よってあらゆる番地にアクセスできそうな雰囲気を出しておく必要があります(ptr+1の部分)。

module top(sysclk, cpu_resetn, sw, led);
    input  wire sysclk;
    input  wire cpu_resetn;
    input  wire[2:0] sw;
    output wire[7:0] led;

    reg[9:0] ptr;
    wire[9:0] next_ptr;
    assign next_ptr = sw[2] ? ptr + 1 : 0;

    ram ram0(
        .raddr(next_ptr),
        .rdata(led),
        .we(sw[0]),
        .waddr(ptr),
        .wdata(led + sw[1]),
        .clk(sysclk)
    );

    always @(posedge sysclk or negedge cpu_resetn) begin
        if (!cpu_resetn) 
            ptr <= 0;
        else
            ptr <= next_ptr;
    end    
endmodule

次に、テストベンチを書きます。20サイクルの間、RAMの0番地から読みだした値+1を書き込むということを繰り返します。

`timescale 1ns / 1ps

module test(led);
    output [7:0] led;
    
    reg      sysclk;
    reg      cpu_resetn;
    reg[2:0] sw;
    
    top top0(sysclk, cpu_resetn, sw, led);

    initial begin
        sysclk <= 0;
        cpu_resetn <= 1;
        sw <= 0;
        #100 cpu_resetn <= 0;
        # 99 cpu_resetn <= 1;
        #100 sw <= 3;
        #200 $finish;
    end

    always #5 sysclk <= ~sysclk;
endmodule

これを実行してみて、RAMの0番地に何が書かれているかを確認してみます。すると、verilogシミュレーションでは20、タイミングシミュレーションでは10になっています。

verilogシミュレーションの結果は、verilog HDLの記法上もそうであるように、ちゃんとwrite-firstになっていることを示唆しています。つまり、書き込んだ値が即読み出せるということです。

一方、タイミングシミュレーションでは、どうもread-firstになっているようです。書き込んだ値は、直後のサイクルではまだ見えないということになります。

解決法

Xilinx社のマニュアルでは、write-firstなBlockRAMに推論してほしい時の書き方として、以下の書き方が推奨されていました。

always @(posedge clk) begin
    if(we) begin
        do <= data;
        mem[address] <= data;
    end else
        do <= mem[address];
end

実際、この書き方をした場合は、おかしな問題は発生しません。

結論

Xilinx社のマニュアルに記載してある通りに書けば、おかしなことは発生しません。

しかし、そうでない書き方をした時に不完全な推論が行われ、verilog HDLの記述どおりに合成されないのは、不親切とかいう次元を超えて明らかにVivadoのバグだと思います。

とはいっても、RAMに書いた直後に読みだすというパターンは普通は存在しないため、問題につながることは少ないようです。 キャッシュメモリをBlockRAMで書いたりすると、書いた*2直後に読み出すということは十分あり得そうです(12/24 指摘を受け追記)。 そもそも、BlockRAMを書くつもりでなくても、クロック同期なメモリっぽいものを書いてしまうと勝手にBlockRAM化されてわけのわからないバグにつながるということも考えられます。

1read1writeのRAMの正しい作り方(12/24 追記)

「解決法」のところで書いたものはreadポートとwriteポートが共有されているものになっています。1read1writeのRAMを作るにはどうしたらいいか?という疑問があったので作ってみました。

まず、以下のBlockRAMのマニュアルの20ページ目を見てみると、たとえwrite-firstを指定したとしても、readポートで指定した番地とwriteポートで指定した番地が同じ場合、readポート側に何が読みだされるかは保証されないと記載されています。 https://www.xilinx.com/support/documentation/user_guides/ug473_7Series_Memory_Resources.pdf

よって、このBlockRAMプリミティブを使うだけでは絶対に1read1writeのRAMを作ることはできません。

そこで、readアドレスとwriteアドレスが一致している場合のフォワーディング回路を自前で用意する必要がありそうです(自動で推論してほしいですね)。

以下のように書くと、ちゃんとBlockRAMに推論され、最初の記述と論理的に同じ動作になり、合成してもこの問題は発生しなくなります。

module ram(raddr, rdata, we, waddr, wdata, clk);
    input  wire[9:0] raddr;
    output wire[7:0] rdata;
    input  wire      we;
    input  wire[9:0] waddr;
    input  wire[7:0] wdata;
    input  wire      clk;

    reg[7:0] mem[0:1023];
    reg[31:0] i; initial for(i=0; i<1024;i=i+1) mem[i] = 0;

    reg[9:0] last_raddr, last_waddr;
    reg[7:0] rdata1, rdata2;
    always @(posedge clk) begin
        if (we) begin
            mem[waddr] <= wdata;
            rdata1 <= wdata;
        end else begin
            rdata1 <= mem[waddr];
        end
        rdata2 <= mem[raddr];
        last_raddr <= raddr;
        last_waddr <= waddr;
    end
    assign rdata = last_raddr == last_waddr ? rdata1 : rdata2;
endmodule

writeポートの方は、「解決法」のところで書いたような、Xilinx社マニュアル推奨のwrite-firstな読み出しになる書き方をしています。 よってそこから読み出した結果rdata1は直前の書き込み値last_wdataの代わりとして使えます。

これによる無駄なフリップフロップ削減効果に加え、そもそもwrite-firstにしたほうが高速と書かれているので、たぶんこれが一番いい書き方なのだと思います。

*1:12/26追記:見れなくなっていました。アーカイブRAMの書き方 - マイクロプロセッサの設計と実装

*2:キャッシュメモリへのwriteアクセス(直後に読み出すことは少なさそう)ではなく、キャッシュラインの置換時のRAM書き込み

分岐命令の意味論一覧

機械語には複数種類の分岐命令があり、しかも一つの分岐命令でも異なる意味論で使ったりします。 どういう意味論を持っているか、まとめてみました。

以下では、無条件分岐命令か、あるいは条件分岐命令で条件が成り立った場合の意味論についてまとめます。 条件分岐命令で条件が不成立だった場合(単に次の命令に進むだけでいつもとおなじ)は除きます。

一つの関数内で完結するもの

無条件分岐

特定の番地にプログラムの制御を移す、もっとも単純な分岐命令です。 一つの関数内で完結することから、分岐元と分岐先のアドレスの差は小さいことが多く、相対アドレス指定が行われることが多いです。

ジャンプテーブルを使った多方向分岐

C言語switch文や、関数型言語のパターンマッチ等を効率的に実装する場合などに出現します。 汎用レジスタに書かれている値への分岐命令として実現されます。 原理的には相対アドレス指定でも問題ないのですが、そうなっている命令セットはあまり見かけません。

関数間にまたがるもの

関数呼び出し

特定の番地にプログラムの制御を移しますが、呼び出し先ルーチンが終了した場合にこの命令の次の命令から再開する、という機能も兼ね備えた分岐命令です。 一命令に複数の機能がありますが、非常に頻出であるため、RISCであっても採用される命令です。 関数呼び出し命令をどのように実装するかは、命令セットによっていろいろな方法が模索されましたが、どの方法も欠点があります。 RISC-Vでは、呼び出し先関数のアドレス指定に、相対アドレス指定を採用しています。これだと遠いアドレスにある関数が呼び出せないように思えますが、相対アドレスの絶対値が大きい場合には、二命令かけて相対アドレスを指定する、という方法をとっています。

間接関数呼び出し

関数呼び出しとほぼ同じですが、呼び出す関数が動的に決まるというものです。例えば、関数ポインタを使った関数呼び出しや、C++の仮想関数呼び出しなどです。

関数からの復帰

関数を終了し、関数呼び出し命令の次の命令に制御を移す命令です。 「関数呼び出し命令の次の命令」のアドレスは、絶対アドレスで保存されているはずです。

末尾呼び出し

特定の番地にプログラムの制御を移しますが、呼び出し先ルーチンが終了した場合にこの関数に戻るのではなく、この関数を呼び出した命令の次の命令から再開するような分岐命令です。 やたらに複雑な命令のように見えますが、(呼び出し規約を無視すれば)無条件分岐命令とまったく同じ動作になっています。 呼び出し先関数のアドレスの指定方法は、関数呼び出しと同じになります。

コルーチン呼び出し

コルーチンは、関数の一般化です。 コルーチンと関数とで異なるのは、以下の二点です。 関数の場合、常に最後まで実行されたのち、呼び出し元に帰ってきます。一方、コルーチンの場合、その途中の任意の場所で実行を中断し、呼び出し元に帰ってくる可能性があります。 関数の場合、次回呼び出し時もまた先頭から実行します。一方、コルーチンの場合、次回呼び出し時は、その中断したところから再開します。

こんな複雑なことを実現する分岐命令も、RISC-Vには存在します。中断したアドレスは、必然的に絶対アドレスが使われます。

なお、コルーチンからの復帰というのは、別のコルーチンの呼び出しに相当するため、ここに含まれます。

機械語を意識して表にまとめてみる

表にまとめてみると、以下のようになります。参考のため、RISC-Vの分岐命令(jjaljrjalrtailcallret)を併記しました。

ただし、以下でRASというのは、return address stackの略で、関数呼び出しに関連する命令の分岐予測を高精度に行うハードウェアのことです。 関数呼び出しに関連するアドレスだけが入れられる、汎用でないレジスタだと思ってもおおよそ間違っていません。

一つの関数内で完結する分岐 関数をまたぐ分岐で、RASに書き込まない 関数をまたぐ分岐で、RASに書き込む
無条件分岐(j 末尾呼び出し(j 関数呼び出し(jal 分岐先が近く、固定
N/A 末尾呼び出し(tail 関数呼び出し(call 分岐先が遠く、固定
ジャンプテーブル分岐(jr 間接末尾呼び出し(jr 間接関数呼び出し(jalr 分岐先が不定
N/A 関数からの復帰(ret コルーチン呼び出し/復帰(jalr RASを読み出し、そこに分岐

継続をあらわに扱った場合

先の表は、機械語・ハードウェアを強く意識したものになっていました。継続を使うと、抽象的なままで考えることができます。

関数からの復帰は、継続の末尾呼び出しなので、間接末尾呼び出しと同じはずです(実際、retjrの特殊な場合です)。 また、ジャンプテーブルも、複数の継続の中から一つの継続を選んで呼び出す、と解釈できるため、やはり間接末尾呼び出しと同じになっています。

継続渡しスタイルに変換してしまえば、無条件分岐と末尾呼び出しの区別はなくなります。また、関数はコルーチンの退化したバージョンに過ぎません。

これに加え、分岐先が遠い・近いというのは完全に機械語の制約であり、本質的なものではないことに注目すれば、分岐命令は本質的に以下のように分類されるということがわかります。

  • 静的に決まっている継続の末尾呼び出し(j
  • 動的に決まる継続の末尾呼び出し(jr
  • 現在の継続を取り出しつつ、静的に決まっている継続を呼び出し(jal
  • 現在の継続を取り出しつつ、動的に決まる継続を呼び出し(jalr

Brainfuck コードを最適化コンパイルする~C++メタプログラミングで~

この記事は、Brainf*ck Advent Calendar 2019 - Adventarの八日目の記事です。

Brainfuckインタプリタで実行するとかなり遅いです。

高速に実行する一つの簡単な方法に、「C言語に変換する」という方法があります。 Wikipediaにも書いてありますが、以下の簡単な文字列置換によって、C言語に変換できます。

  • >++ptr;
  • <--ptr;
  • +++*ptr;
  • ---*ptr;
  • ,*ptr=getchar();
  • .putchar(*ptr);
  • [while(*ptr){
  • ]}

あとは、適当な“おまじない”をくっつければコンパイルできるC言語ソースコードになります。

#include <stdio.h>

int main() {
    char data[30000] = {};
    char* ptr = data;
    【ここに変換したコード本体】
}

さて、この方法では変換という手順が入るのが気になります。何とか変換せずに、“おまじない”をくっつけるだけでコンパイルできるようにできないでしょうか?

それ自体は簡単で、単にインタプリタソースコードとくっつければ、コンパイルできるC言語ソースコードになります(第1二村射影)。 しかし、これは結局インタプリタと同程度の速度でしか動作せず、高速化したいという当初の目的を達成できていません。

「変換せず、“おまじない”をくっつけるだけでコンパイル可能」と「コンパイルして高速なバイナリを得たい」の両立はできないのでしょうか?

C++メタプログラミングを駆使して、これを実現します。

前提

使ったメタプログラミングテクニック

コンパイル時に値を計算するのではなく、プログラムコードを計算すればいいだけなので、困難な点は少ないです。

完成品

#include <cstdio>
#include <string>
#include <utility>

using Byte = int;

constexpr std::size_t find_open( const char* bf_code, std::size_t begin, std::size_t end ) {
        for( std::size_t i = begin; i != end; ++i ) {
                if( bf_code[i] == '[' ) { return i; }
        }
        return end;
}

constexpr std::size_t find_close( const char* bf_code, std::size_t begin, std::size_t end ) {
        std::size_t depth = 0;
        for( std::size_t i = begin; i != end; ++i ) {
                if( bf_code[i] == '[' ) { ++depth; }
                if( bf_code[i] == ']' ) { --depth; }
                if( depth == 0 ) { return i; }
        }
        return end;
}

template<char C>
Byte* action( Byte* ptr ) {
        static_assert( C != ']', "Unmatched ']'" );
        return ptr;
}

template<> Byte* action<'+'>( Byte* ptr ) { ++*ptr; return ptr; }
template<> Byte* action<'-'>( Byte* ptr ) { --*ptr; return ptr; }
template<> Byte* action<'>'>( Byte* ptr ) { return ptr+1; }
template<> Byte* action<'<'>( Byte* ptr ) { return ptr-1; }
template<> Byte* action<','>( Byte* ptr ) { *ptr = std::getchar(); return ptr; }
template<> Byte* action<'.'>( Byte* ptr ) { std::putchar(*ptr); return ptr; }

template<std::size_t Begin, class Code, std::size_t... Indeces>
Byte* ExecuteSimple( Byte* ptr, Code bf_code, std::index_sequence<Indeces...> ) {
        ((ptr = action<bf_code()[Begin+Indeces]>( ptr )), ...);
        return ptr;
}

template<std::size_t Begin, std::size_t End, class Code>
Byte* ExecuteRange( Byte* ptr, Code bf_code ) {
        constexpr std::size_t Open = find_open( bf_code(), Begin, End );
        ptr = ExecuteSimple<Begin>( ptr, bf_code, std::make_index_sequence<Open-Begin>() );
        if constexpr( Open != End ) {
                constexpr std::size_t Close = find_close( bf_code(), Open, End );
                static_assert( Close != End, "Unmatched '['" );

                while( *ptr ) {
                        ptr = ExecuteRange<Open+1, Close>( ptr, bf_code );
                }
                ptr = ExecuteRange<Close+1, End>( ptr, bf_code );
        }
        return ptr;
}

int main() {
        constexpr auto bf_code = []{return
#include "mandelbrot.bf"
        ;};
        constexpr std::size_t CodeSize = std::char_traits<char>::length( bf_code() );
        Byte data[30000] = {};
        Byte* ptr = data;
        ExecuteRange<0, CodeSize>( ptr, bf_code );
}

解説

using Byte = int;

Brainfuckデータメモリ一要素を表す型の定義です。後述しますが、charとするよりもintとしたほうが高速です。ただし、ラップアラウンドを前提としているBrainfuckプログラムは正しく動作しなくなります。

find_openfind_close

find_open(bf_code, begin, end)は、Brainfuck命令列bf_codeを受け取り、そのbegin文字目*1以上end文字目未満の間で最も近い位置にある[を探し、その位置を返します。見つからなかった場合は、endを返します。

find_close(bf_code, begin, end)は、Brainfuck命令列bf_codeを受け取り、そのbegin文字目以上end文字目未満の間で最も近い位置にある[と対応の取れた]の位置を返します。該当しない場合、endを返します。

以上の二つの関数は、自明に副作用がない純粋な関数であり、また単純型のみで型付けできるため、constexpr関数として実装しています。

action

Brainfuck命令をテンプレート引数として受け取り、その動作をする関数に実体化されるテンプレートです。 全ての動作は、現在のデータメモリを指し示すポインタを受け取り、目的の動作を行い、更新されたデータメモリを指し示すポインタを返す関数として実装します。

6つの命令+-><,.について特殊化することでディスパッチを実現しています。 []の命令動作は、ローカルには完結しないため、ここでは取り扱いません(]が来た場合、対応が取れていない]であるため、コンパイルエラーにするstatic_assertを入れています)。

ExecuteSimple

ExecuteSimpleはテンプレートで、ExecuteSimple<Begin, Code, Indeces...>は、Brainfuck命令列を返す関数の型Code、その中で[]を含まない非負整数区間の始点を指定するBegin区間と同じ要素数の初項0公差1の数列を指定するIndecesを受け取り、その一連の動作をする関数に実体化されます。 テンプレート引数を推論する都合上、受け取る順番が変則的です。

Indecesは、0, 1, 2, 3, 4, 5のような非負整数列になっており、((【Indecesを使った式】), ...);の形で使うことで、【Indecesが0の時の式】, 【Indecesが1の時の式】, ……, 【Indecesが5の時の式】;のように展開されます。 このように展開されることは、ループアンローリングを強制していることに相当します。このループアンローリングは、コンパイラの最適化を助けます。

ExecuteRange

ExecuteRangeはテンプレートで、ExecuteRange<Begin, End, Code>は、Brainfuck命令列を返す関数の型Code、その中で実行すべき非負整数区間*2を指定するBeginEndを受け取り、その動作をする関数に実体化されるテンプレートです。 このテンプレートの実体化の際、再帰的にこのテンプレートを、異なるテンプレート引数で実体化することがあります。

まず、find_openで実行すべき区間の中に[があるかを確認します。もし[がなければ、ExecuteSimpleを呼び出すだけで仕事は完了です。

もし、[が含まれていれば、対応する]find_closeで探します。これは正しいBrainfuckプログラムであればあるはずです(もしなければ、対応が取れていない[であり、コンパイルエラーにします)。

[]で囲まれている区間は、そのBrainfuckプログラムを、while文の中で再帰的にExecuteRangeを呼び出すことで変換します。 また、残りのBrainfuckプログラムも、再帰的にExecuteRangeに渡し、変換します。

if constexprを使うことにより、テンプレートの再帰的展開が無限ループに陥らないようにすることができます(使わない場合、実行されないif文の中身も展開されてしまって無限ループになります)。

main

ゼロ初期化されたデータメモリの確保、Brainfuckデータメモリを指し示すポインタの確保、変換されたプログラムの呼び出しを行うだけです。

bf_codeがへんてこな書き方になっているのは以下の理由です。 actionExecuteSimpleExecuteRangeはいずれもテンプレート引数により動作を変化させる必要があります。 そのため、異なるBrainfuck命令列には、異なる型がつかないとうまくいきません。 C++ラムダ式は、すべて異なる型がつくことになっているので、定数を返すラムダ式を書くことで、これを実現できます。 ExecuteSimpleExecuteRangeが、Brainfuck命令列ではなく、それを返す関数を受け取るような関数になっているのもこれが原因です。

速度測定

  • ベンチマークとしてErik Bosman氏作mandelbrot.bを使用
  • 環境:Windows10 (1903), WSL (Microsoft 4.4.0-18362.476-Microsoft 4.4.35), Intel(R) Core(TM) i7-6700HQ CPU @ 2.60GHz(実際には、Intel(R) Turbo Boost Technology 2.0により3.50GHzで稼働しているはずです)
  • C++コンパイラ:g++-9 (Ubuntu 9.2.1-17ubuntu1~16.04) 9.2.1 20191102
  • C++コンパイル方法:g++-9 -std=c++17 -O3 -march=native
  • 実行速度の計測:time ./a.out >/dev/nullを20回実行し、最も速かったものを採用

何も工夫しないインタプリタの場合

Brainfuckのインタプリタを書いた - よーるで、ナイーブなインタプリタと、g++-9以降向けに最適化した*3ナイーブなインタプリタを実装しました。 これらを使った場合、前者では約40秒、後者では約27秒かかりました。

C言語に変換した場合

最初に述べた方法でC言語に変換したソースコードコンパイルすると、その実行時間は0.953秒でした。

データメモリ用の配列をcharではなくintで宣言した場合、実行時間はわずかに縮まり、0.906秒ほどになりました。

今回作った、C++メタプログラミングによる方法

実行速度は0.844秒まで高速化されました。

Byteとしてcharを使った場合、0.875秒とほんのわずかですが劣ります。これは、コンパイラcharを取り扱うのが苦手であることに起因します。 x86は1Byteレジスタを持つため本来は4Byteレジスタとの間で転送する必要はないのですが、コンパイラがそれを読み切れず、不要な変換命令を挟んでしまうことが原因のようです。

ちなみに、clang++-9を用いた場合は、0.953秒でした。

ところで、言うまでもないことですが、32倍もの高速化が行われる理由は、そのほとんどがC++コンパイラの最適化能力にあります。 今回書いたコードは、Brainfuckのコードをちょっとパースして、C++コンパイラにわかりやすい形に変換したにすぎず、実際に最適化を行っているわけではありません。

なぜC言語に変換するよりC++メタプログラミングを使ったほうが高速なのか

謎です。

C++メタプログラミング版では、ループ内部が関数化されている影響で生存解析等の最適化が行いやすくなった結果のようにも思われますが、確定的なことはよくわかりません。

もっと高速化することは可能か

Brainfuckプログラムに頻出のパターンをとらえ、それを高速化するようなコンパイラ(あるいはメタプログラミングを駆使した“おまじない”)を書くことは可能です。

例えば、Brainfuck のプログラムを C に変換してコンパイルしてみる実験 - Qiitaで行われているようなC言語変換をしてからコンパイルしたバイナリを実行した場合、その実行時間は0.609秒と、今回作った方法を上回ります。

Brainfuckプログラムに頻出のイディオムを、今回作った方法で最適化できているか確認してみます。 結論から言うと、簡単なテストプログラムの場合、全部最適化できてしまいました。 上記記事のC言語変換ではできているが、コンパイラの汎用的な最適化では見抜けないというパターンは見つけられませんでした。

[-]

ゼロ代入イディオムです。今回作った方法では、ちゃんと最適化されます。

入力Brainfuckプログラム,[-].とした場合、コンパイル結果は以下のようになりました。

        mov     rdi, QWORD PTR stdin[rip] # 第一引数=stdin
        call    _IO_getc
        mov     rsi, QWORD PTR stdout[rip] # 第二引数=stdout
        xor     edi, edi # 第一引数
        call    _IO_putc

xor edi, edix86のゼロ代入イディオムで、確かに最適化されていることがわかります。

[->+<]

破壊的加算イディオムです。-+の位置が違ったり、<>の個数が異なったり入れ替わっていたりするような場合もありますが、すべて同じです。,>,[-<+>]<.のようなごく単純な例では、多少見抜いて最適化されるようです。

        mov     rdi, QWORD PTR stdin[rip]
        call    _IO_getc
        mov     rdi, QWORD PTR stdin[rip]
        mov     ebp, eax # 一回目のgetcharで返ってきた値をebpに保存
        call    _IO_getc
        test    eax, eax # 二回目のgetcharで返ってきた値が0であるかチェック
        mov     edx, ebp # 一回目のgetcharで返ってきた値をedxにコピー
        lea     ebp, [rbp+0+rax] # 二回目のgetcharで返ってきた値とebpに保存した値を足し算
        cmove   ebp, edx # 二回目のgetcharで返ってきた値が0の時はedx、そうでないときは加算結果
        mov     rsi, QWORD PTR stdout[rip]
        mov     edi, ebp # ……をputcharの引数に渡す
        call    _IO_putc

getcで返ってきた値が0の場合でも足し算した値をそのまま使えばいいのに、そこは見抜けなかったようです。

ちなみに、clangでコンパイルした場合、この部分を以下のように完全に最適化してくれます。

        mov     rdi, qword ptr [rip + stdin]
        call    _IO_getc
        mov     ebx, eax # 一回目のgetcharで返ってきた値をebxに保存
        mov     rdi, qword ptr [rip + stdin]
        call    _IO_getc
        lea     edi, [rax + rbx] # 二回目のgetcharで返ってきた値とebxに保存した値を足し算し、
                                 # putcharの引数に渡す
        mov     rsi, qword ptr [rip + stdout]
        call    _IO_putc

gccは、この手の「結局同じになるのだから分岐しなくていい」というのを見抜けませんが、clangは見抜いてくれます。 総合的な最適化能力はgccのほうが上のようですが……。

まとめ

C++メタプログラミングを使い、Brainfuckソースコードにちょっと“おまじない”をくっつけたソースコードコンパイルするだけで高速なバイナリを得る方法を示しました。

それによって得られるバイナリは、従来の「C言語に変換する」手法を超える最適化がかかったものになっていました。

また、この手の別言語コンパイラを流用した最適化を行う場合、Brainfuckイディオムを特別扱いせずとも単純な例であればコンパイラ最適化が働くこと、 データメモリ用の配列を4Byteなどコンパイラが得意としている語長にすると更なる高速化が行われるということがわかりました。


明日は、roodni さんが1993年のbrainfuckについて書かれるそうです。1993年ってちょっと特別ですね。

*1:0-origin

*2:左閉右開半区間

*3:分岐予測器の予測的中率が高くなるようなコードが出力されるように書いたコードになっています。

容量上限付き可変長配列の空間計算量

boost::container::static_vectorのような動作をする、容量の上限Nコンパイル時に決められた値であり、実行時に配列の要素数sizeを(上限までの範囲で)変更できるようなデータ構造を考えます。 これを実現するには最低何bitの空間が必要でしょうか?

要素一つ当たり必要なビット数を、Kbitだとします。

 N=0 の場合

これは、自明に0bitで実現可能です。

 K=0 の場合

これは、自明に \lfloor \log_2(N+1) \rfloorbitです。配列の要素数size0である場合、1である場合、……Nである場合、の N+1通りを区別する必要があるためです。

それ以外の場合で、 2^K \ge Nの場合

この場合、 NK + 1bitで実現可能です。

まず、この値が下限値であることを示します。 要素一つ当たりに必要なビット数がKbitで、最大N要素入りうるので NKbitは確実に必要です。 N要素入っている状態は \left(2^K\right)^N = 2^{NK}状態あり、それ以外に0要素入っている状態が少なくともありますから、最低でも 2^{NK}+1通りを区別する必要があります。  NKbitでは 2^{NK}通りの状態しか区別できませんので、最低でも NK+1bit必要ということになります。

次に、実現する方法を示します。 まず、一エントリがKbitであり、N要素ある配列を普通に用意します。これには NKbitが費やされます。残りの1bitは、「sizeNと等しい」ことを表すフラグにします。 容量をケチるため、フラグの状態がtrueであるかfalseであるかによって、データ構造を変化させます。具体的には、以下のようにします。

  • フラグがtrueである、つまりsizeNと等しい場合、用意した配列にsize個のKbitデータを詰め込む
  • フラグがfalseである、つまりsizeNより小さい場合、用意した配列にsize個のKbitデータを詰め込む。配列の末尾にはsizeを書き込む

フラグがtrueである場合、問題が発生しないことはすぐにわかります。フラグがfalseである場合について少し補足します。

  • 「配列の末尾には」となっていますが、sizeNより小さいので、配列の末尾一エントリは使われていません。よって、必ず使えます。
  • sizeを書き込む」となっていますが、sizeとしてあり得るのは01、……N-1のN通りなので、(二進数として書き込むことで)Kbitで足ります。

それ以外の場合

もし、要素数取得関数size()の時間計算量が定数である必要がある*1なら、以下のケチりテクニックは使えません。 なくてもよいなら、上記と同じ NK + 1bitで実現する方法があります。ただし、size()の時間計算量が定数でなくなるということからもわかるように、かなり非効率な実装となります。

まず、先と同じ議論により、 NK+1bitが下限値であることはすぐにわかります。

実現する方法も、先の方法とおおよそ同じです。フラグがfalseの場合の議論が成り立っていないので、以下のように変更します。

  • フラグがfalseである、つまりsizeNより小さい場合、用意した配列にsize個のKbitデータを詰め込む。使っていないエントリの先頭にtrueを、他のエントリにはfalseを書き込む

少し補足します。

  • 「使っていないエントリの先頭に」となっていますが、sizeNより小さいので、使っていないエントリは常に存在します。
  • trueを」「falseを書き込む」となっていますが、 K>0であるので、そのエントリはtruefalseを区別できるビット幅を持っていることが保証されます。ただし、タグなしunionという、かなり危険なコードを書く必要があります。
  • 「ほかのエントリには」となっていますが、そのようなエントリがなければ何もする必要はありません。

このように書き込んだデータ構造を解釈するには、後ろから解釈する必要があります。前から解釈することは不可能です。

後ろから解釈する、というのはつまり、フラグ→配列の末尾→配列の末尾のその前→……という順番で確認していくということです。 このように確認を進めていき、最初にtrueが書き込まれていた部分までに並んでいたfalseの数が、N-sizeに対応しています。 つまり、sizeの情報が一進数で書き込まれているということになります。

別の解釈もできます。フラグを「trueなら配列にはN個の要素が書き込まれている、falseなら配列には『容量の上限がN-1の"この仕組み"のデータ構造』が書き込まれている」という切り替えフラグだと解釈すると、再帰的なデータ構造になっていることが確かめられます。基底ケース(N=0)ではfalseがあり得ないので常にtrueとなっていることとも整合します。

参考:要素としてあり得る状態が 2^K通り未満である場合

この場合、きっかり NKbitで実現可能です。その実例は、ヌル文字(NUL)終端文字列であり、文字列以外にも一般化可能です。

このように書き込んだデータ構造は、前から解釈することになります(一種の接頭符号になっています)。

size()の時間計算量が定数でないという特徴は似ていますが、前から解釈できる特徴により、コンパイル時に容量上限を決めなくてよいという利点があります。

*1:C++STLコンテナの要件。size関数を持たないコンテナも存在します(例えばstd::forward_list)

Microsoft(R) Excel(R) が有効数字を減らしてくる現象について

まずA1セルに1を入れ、A2セルには=A1*2と記入します。次に、A2セルのフィルハンドル(セルの右下にカーソルを合わせたときに出現する十字型のカーソル)を下にドラッグし、100セルほどコピーします。

このようにすると、2Nという数列を計算することができます。

しかし、A列の書式設定を「数値」にしてみると、途中から一の位が0になっていることがわかります。 2Nという数列のいかなる項も、一の位が0になることはありませんから、これはおかしいです。

内部で浮動小数点数演算をしているため、誤差が生じているのだ、といった考え方もあると思います。 しかし、近年の計算機が使う浮動小数点数の基数はほぼまちがいなく2であり、2Nはオーバーフローしない限りは正確に表せるはずです。 実際、後で示すように、Excelが使用している浮動小数点数は、ほぼ確実にIEEE 754 binary64(二進倍精度浮動小数点数)です。

実験

どこから正しく表示できないのか

まず、この現象が発生し始めるのはどこかを調べます。 すると、=10^15+1という数式を入力したセルは、その表示が1000000000000000となり、この現象が発生していることがわかります。 一方、=10^15-1という数式を入力したセルは、その表示が999999999999999となるため、この現象は発生していないことがわかります。

つまり、「十進表記に変換する際、その数が正確に表せるかとは無関係に、上から15桁は十進表記し、残りの桁は0で埋めて桁合わせする」というアルゴリズムになっているようです。

倍精度浮動小数点数と仮定して、どこから誤差が生じるのか

ところで、次の浮動小数点数との間隔が1を超えるようになるような最小の正の浮動小数点数は、253です。 253+1は、浮動小数点数では正確に表せません。 逆に言うと、絶対値が253未満の正の浮動小数点数について、その数より1大きい数は、浮動小数点数で表せます。

253は十進法で表すと9007199254740992ですが、Excelで計算させると最後の2が0に変わってしまいます。

内部では、実は正確に保持している

1000000000000000を入力したセルに1を足し続けるといずれ1000000000000010になります。 つまり、見た目では変な表記になっていますが、内部的には正確な値を保持しているようです。

実は、グラフを用いると、これを”可視化”することができます。

f:id:lpha_z:20191117234956p:plain

横軸をX、縦軸をYとしたとき、Y=253+Xであるような点をいくつか打ちました。

Xが負の領域では、きれいに一直線に並んでおり、正確な値が計算されていることがわかります。ただし、縦軸の目盛数値は、最終桁が0に変わってしまっています。

Xが正の領域では、がたがたとして変な感じになっています。これは、以下の二点によるものと考えるとすべての説明がつきます。

  • 253より大きな、浮動小数点数で表せる数は、253+2, 253+4, 253+6, ……のように続くので、正確に表せないときがある
  • 演算結果浮動小数点数で正確に表せない場合、最も近い浮動小数点数に丸める。ただし、そういった浮動小数点数の候補が二個ある場合、仮数部の末尾が0であるようなものを選ぶ(いわゆる偶数丸め)

やはり、ExcelIEEE 754 binary64を内部的に使っていると考えるのが妥当なようです。

強制的に15桁しか表示しないのは、十進浮動小数点数で表した時の仮数部が1である場合に精度がほぼ1桁分減ってしまうのでやめてほしいのですが……。