完全なTRNGではなく、個人が趣味でTRNGを実装評価しているだけです。
C言語で疑似乱数を使うとき、昔からよく見るコードがある。

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

int main(void) {
    srand((unsigned)time(NULL));

    printf("%d\n", rand());
    return 0;
}

C言語の入門本や授業などで最初に習うと思う。

time(NULL)は秒単位の時刻であり、探索空間が狭い。
プログラムが実行された時刻がだいたい分かれば、seedを総当たりできる。
しかもrand()自体も暗号用途を想定した乱数生成器ではない。

そこで今回は、FPGA上にリングオシレータ TRNGを作って、物理現象由来のエントロピーを取り出してみる。

真正乱数生成

True Random Number Generatorといい、頭文字を取ってTRNGである。

TRNGは物理現象をエントロピー源にする。
イメージしやすいものはサイコロである。
サイコロは投げるという物理的な動作によって1から6までの数値を得ることができる。

コンピュータでトランジスタの熱雑音やクロック差などを用いることもあるようだ。

FPGAで比較的作りやすいものとしては、リングオシレータがある。

リングオシレータの原理

詳しいことは ACRiルームのブログ記事 が非常にわかりやすかった。

簡単に要約すると、NOT回路を奇数個繋いでループにすることで、回路を一周すると信号が反転して、発振するようになる。

このとき、配線遅延などが発生するため、きれいな発振とはならず、ジッタが生じる。
この出力をサンプリングすることで乱数的な出力を得られるという仕組みのようだ。

実装

PYNQ-Z2ボード いつものFPGA検証ボード PYNQ-Z2 のPL側にリングオシレータを合成し、AXI-StreamでAXI DMAへ流してPS側で色々試してみる。

なお、いつもは高位合成(HLS)を用いて実装を行ってきたが、リングオシレータは同期回路ではないため、同期回路を生成するHLSとは非常に相性が悪い。
仕組みとしてもNOTを奇数個接続するものと、簡易的であるため、HDLで記述した。
(あと、少し前まではAIのHDLの精度はpythonやjsと比較すると残念なモノであったが、最近はそれなりに動くものが出てくる。)

リングオシレータ

`timescale 1ns / 1ps
`default_nettype none

(* DONT_TOUCH = "true", KEEP_HIERARCHY = "yes" *)
module ro_cell #(
    parameter integer STAGES = 5
)(
    input  wire enable,
    output wire ro_out
);

    (* KEEP = "true", DONT_TOUCH = "true" *)
    wire [STAGES-1:0] ring;

    (* KEEP = "true", DONT_TOUCH = "true" *)
    assign ring[0] = enable ? ~ring[STAGES-1] : 1'b0;

    genvar i;
    generate
        for (i = 1; i < STAGES; i = i + 1) begin : gen_inv
            (* KEEP = "true", DONT_TOUCH = "true" *)
            assign ring[i] = ~ring[i-1];
        end
    endgenerate

    (* KEEP = "true", DONT_TOUCH = "true" *)
    assign ro_out = ring[STAGES-1];

endmodule

`default_nettype wire

Vivadoでは同期回路を作ることが想定されているためか、こんな異常な処理を行う処理はエラーが表示された。

DRC LUTLP-1: Combinatorial Loop Alert

これを回避するためには、制約ファイルで許可する必要がある。

set_property ALLOW_COMBINATORIAL_LOOPS true [get_nets -hierarchical -regexp {.*u_ro_cell/ring\[4\]}]

(色々記述したが、Vivadoによって制約ファイルが最適化された)
これを入れないと、配置配線までは進んでもbitstream生成直前で止まってしまった。

AXI DMAでPYNQへ流す

TRNGから出力されたビット列は32bit長ワードにしてAXI4-Streamで出力した。

m_axis_tdata[31:0] : 32bitのraw entropy word
m_axis_tvalid      : tdataが有効なとき1
m_axis_tready      : AXI DMA側が受け取り可能なとき1
m_axis_tlast       : 1フレームの終端

今回のRTLでは、1024ワードごとにTLASTを出すようにした。 つまり1回のDMA転送単位は4096Bである。

AXI-Stream出力

AXI-Stream側のトップは次のような構成にした。

module ro_trng_entropy_axis #(
    parameter integer NUM_RO          = 31,
    parameter integer STAGES          = 5,
    parameter integer STARTUP_DISCARD = 1024,
    parameter integer FRAME_WORDS     = 1024
)(
    input  wire        aclk,
    input  wire        aresetn,

    input  wire        enable,

    output reg  [31:0] m_axis_tdata,
    output reg         m_axis_tvalid,
    input  wire        m_axis_tready,
    output reg         m_axis_tlast,

    output reg  [31:0] raw_bit_count,
    output reg  [31:0] raw_word_count,
    output reg  [31:0] dropped_bit_count,
    output reg  [31:0] frame_count
);

NUM_ROはリングオシレータの個数、STAGESは各リングオシレータの反転段数である。
今回はNUM_RO=31, STAGES=5で試した。

STARTUP_DISCARD=1024として、enable直後の1024個は安定性を重視して捨てる。

また、FRAME_WORDS=1024として、1024ワードごとにm_axis_tlastを出すようにした。

ブロックダイヤグラム

Vivadoでブロックダイヤグラムを作成し、 作成したRTLを右クリックして、ブロックダイヤグラムに持っていくと配置できる。

ブロックダイヤグラム

これを合成・配置配線・ビットストリーム生成をする。

PYNQで実験

PYNQ側では、AXI DMAの受信チャンネルを使ってDDRへ書き込ませる。

from pynq import Overlay, allocate
import numpy as np

ol = Overlay("ro_trng.bit")
dma = ol.axi_dma_0

buf = allocate(shape=(1024,), dtype=np.uint32)

dma.recvchannel.transfer(buf)
dma.recvchannel.wait()

data = np.array(buf, dtype=np.uint32)

print(data[:16])
print("unique:", len(np.unique(data)))
print("mean:", data.mean())

実行すると、次のような値が取れた。

[3471455669  203510213   21924752 1044509769 1509687426 1237344356
 2261647426  801325290  220333620  503688473   67371577 3557916881
 3356116487  799413274  307927570 2074307515]
unique: 1024
mean: 1819063105.256836

1024ワードすべてが異なる値になっており、少なくとも同じ値を繰り返しているわけではなさそう。

FRAME_WORDSとDMA転送長の罠

最初、262,144ワードのバッファを用意して受信したところ、ほとんどが0になった。

一瞬TRNGが壊れているように見えるが、これはDMA転送長とTLASTの扱いを勘違いしていた。

今回のRTLではFRAME_WORDS=1024なので、1024ワードごとにTLASTが出る。 AXI DMAはTLASTを受け取ると、そのパケットの転送が完了したと判断する。

つまり、262,144ワードのバッファを指定しても、最初の1024ワードだけが書き込まれ、残りは未書き込みのままだったようだ。

大量に取得したい場合は、1024ワードずつ複数回DMAする。

FRAME_WORDS = 1024
N_FRAMES = 256

buf = allocate(shape=(FRAME_WORDS,), dtype=np.uint32)
chunks = []

for i in range(N_FRAMES):
    dma.recvchannel.transfer(buf)
    dma.recvchannel.wait()
    chunks.append(np.array(buf, dtype=np.uint32).copy())

data = np.concatenate(chunks)

print("frames:", N_FRAMES)
print("total words:", len(data))
print("unique words:", len(np.unique(data)))

エントロピービットの分布を見る

取得したuint32列をbit列に変換し、1の比率を調べる。

raw_bytes = data.view(np.uint8)
raw_bits = np.unpackbits(raw_bytes)

ones = raw_bits.sum()
total = raw_bits.size

print("bits:", total)
print("ones:", int(ones))
print("zeros:", int(total - ones))
print("one ratio:", ones / total)

256フレーム、合計262,144ワード(8,288,608bit)で評価した結果は次のようになった。

frames: 256
total words: 262144
frame one ratio mean: 0.4284254312515259
frame one ratio min: 0.420806884765625
frame one ratio max: 0.43475341796875
frame one ratio std: 0.002683141734286474
total one ratio: 0.4284254312515259
unique words: 262125

unique wordsは262144ワード中262125が異なる値だった。

一方で、bitの1比率は約0.428だった。 理想的な乱数なら0.5付近になるはずなので、エントロピーの生のビット列には明確な偏りがある。

詳細に考察するために、最初の16ワードを2進数で表示した。

for i, x in enumerate(data[:16]):
    print(f"{i:04d}: {int(x):08x}  {int(x):032b}")
0000: e7ab7228  11100111101010110111001000101000
0001: b83a60c9  10111000001110100110000011001001
0002: d6c0eb83  11010110110000001110101110000011
0003: 01244927  00000001001001000100100100100111
0004: fc91a000  11111100100100011010000000000000
0005: f8685e27  11111000011010000101111000100111
0006: 84b06724  10000100101100000110011100100100
0007: be7b8202  10111110011110111000001000000010
0008: 11440728  00010001010001000000011100101000
0009: 05500394  00000101010100000000001110010100
0010: 220312c9  00100010000000110001001011001001
0011: 08c20124  00001000110000100000000100100100
0012: 2678968a  00100110011110001001011010001010
0013: 10393676  00010000001110010011011001110110
0014: 022438d4  00000010001001000011100011010100
0015: ede67a04  11101101111001100111101000000100

ぱっと見、0が多い気がする。

ヒストグラムも出力してみた。 エントロピービット列のヒストグラム

この時点でTRNGによる乱数(っぽい値)の取得はできているので、終わりにしてもいいが、 乱数として使えるレベルにしてみたい。

2bitごとでの分布を見てみる

00,01,10,11の2ビットごとで観測してみたいと思う。

bits = np.unpackbits(data.view(np.uint8))
pairs = bits[:len(bits)//2*2].reshape(-1, 2)

pair00 = np.count_nonzero((pairs[:,0] == 0) & (pairs[:,1] == 0))
pair01 = np.count_nonzero((pairs[:,0] == 0) & (pairs[:,1] == 1))
pair10 = np.count_nonzero((pairs[:,0] == 1) & (pairs[:,1] == 0))
pair11 = np.count_nonzero((pairs[:,0] == 1) & (pairs[:,1] == 1))

print("00:", pair00)
print("01:", pair01)
print("10:", pair10)
print("11:", pair11)
00: 1382325
01: 1014563
10: 1015502
11: 781914

2ビットごとの頻度

これより、00が多く、11が少ないので、やはりエントロピービット列は0に偏っている。
しかし、0110はいい感じにちらばっているため、乱数として利用できそうである。

ノイマン抽出

ノイマン抽出 では、2bitずつ見て次のように処理する。

  • 00 破棄
  • 11 破棄
  • 01 => 0
  • 10 => 1
valid01 = (pairs[:, 0] == 0) & (pairs[:, 1] == 1)
valid10 = (pairs[:, 0] == 1) & (pairs[:, 1] == 0)

mask = valid01 | valid10

# 01 -> 0, 10 -> 1
vn_bits = pairs[mask, 0]

print("raw bits:", len(bits))
print("vn bits:", len(vn_bits))
print("vn output rate per raw bit:", len(vn_bits) / len(bits))
print("vn one ratio:", vn_bits.mean())

結果

raw bits: 8388608
vn bits: 2030065
vn output rate per raw bit: 0.24200260639190674
vn one ratio: 0.500231273382872

ノイマン抽出後のヒストグラム エントロピービット列では1の比率が約0.428だったが、抽出後は約0.50023になった。

出力効率は約24.2%である。 つまり、エントロピービット列の約4分の1が抽出後の乱数として残る。

まとめ

今回、PYNQ-Z2のPL上にリングオシレータTRNGを作り、AXI DMAでPS側へエントロピービット列を流してみた。

そのままではエントロピーは0に偏っていたが、0110の頻度はほぼ同じで、古典的なノイマン抽出を行い、効率は悪いが乱数っぽくなった。

今回作ったものは、分類としては簡易的なMURO型TRNGに近いと思う。
複数のリングオシレータをXORで混合し、その出力をZynqのFCLKでサンプリングしている。
ただし、ACRiの記事で紹介されているような厳密なMURO型とは異なり、RO由来の分周クロックなどは使っていない簡易的な実装である。

今後、NIST SP800-90やAIS31検定を試してみる。