ORANGE pico の乱数のひみつ

[以前の記事](https://qiita.com/mikecat_mixc/items/ac1905163f9f04c7e29a)で、[ORANGE pico](http://www.picosoft.co.jp/orange/) の RAM の先頭に乱数に関係すると予想されるデータが格納されていることを紹介した。
今回は、このデータと乱数の関係を観察し、ORANGE pico の乱数の性質を調べる。

# RAM の先頭のデータと乱数

以下のプログラムにより、`rnd` 関数により乱数の値を得ながら、RAMの先頭16バイトのデータを出力する。
RAMのデータは、乱数の値の取得前に1回、そして取得するたびに1回ずつ出力する。
乱数に影響を与えないよう、乱数ではなく固定の数値列を使用してメモリー配列の位置を取得する。

```
10 memarrayaddr=&H03C8:dim array(3)=[&HA436,&H63DB,&HF632,&H4AE6]
20 for i=0 to 3
30 mpoke i*2+0,array(i)
40 mpoke i*2+1,array(i)>>8
50 next
60 offset=0
70 offset=offset-1
80 if array(offset)<>array(0) then goto 70
90 if array(offset+1)<>array(1) then goto 70
100 if array(offset+2)<>array(2) || array(offset+3)<>array(3) then goto 70
110 offset=offset-memarrayaddr/2
120 uartput 1,"---------- "
130 for i=0 to 7
140 uartput 1,format$(" %02X %02X",array(offset+i) & &HFF, (array(offset+i)>>8) & &HFF)
150 next
160 uartput 1,chr$(&H0A)
170 for c=2 to 20
180 uartput 1,format$("&H%08X ",rnd(&H10000000))
190 for i=0 to 7
200 uartput 1,format$(" %02X %02X",array(offset+i) & &HFF, (array(offset+i)>>8) & &HFF)
210 next
220 uartput 1,chr$(&H0A)
230 next
```

実行すると、UARTに以下の出力が得られた。

```
----------  33 13 49 05 B5 3B 12 1F E5 55 9A 15 15 CD 5B 07
&H0CA345EA  EA 45 A3 DC 33 13 49 05 B5 3B 12 1F E5 55 9A 15
&H0B5116E6  E6 16 51 1B EA 45 A3 DC 33 13 49 05 B5 3B 12 1F
&H051049AA  AA 49 10 95 E6 16 51 1B EA 45 A3 DC 33 13 49 05
&H088D00B0  B0 00 8D D8 AA 49 10 95 E6 16 51 1B EA 45 A3 DC
&H0EC7825E  5E 82 C7 1E B0 00 8D D8 AA 49 10 95 E6 16 51 1B
&H0DB24146  46 41 B2 8D 5E 82 C7 1E B0 00 8D D8 AA 49 10 95
&H0AF81443  43 14 F8 9A 46 41 B2 8D 5E 82 C7 1E B0 00 8D D8
&H0AC00F2C  2C 0F C0 2A 43 14 F8 9A 46 41 B2 8D 5E 82 C7 1E
&H0837AD58  58 AD 37 08 2C 0F C0 2A 43 14 F8 9A 46 41 B2 8D
&H07906569  69 65 90 17 58 AD 37 08 2C 0F C0 2A 43 14 F8 9A
&H0D9031D4  D4 31 90 4D 69 65 90 17 58 AD 37 08 2C 0F C0 2A
&H0703EE25  25 EE 03 67 D4 31 90 4D 69 65 90 17 58 AD 37 08
&H02EBD2F0  F0 D2 EB D2 25 EE 03 67 D4 31 90 4D 69 65 90 17
&H06C45EE9  E9 5E C4 46 F0 D2 EB D2 25 EE 03 67 D4 31 90 4D
&H0A16D974  74 D9 16 8A E9 5E C4 46 F0 D2 EB D2 25 EE 03 67
&H021C7CD5  D5 7C 1C F2 74 D9 16 8A E9 5E C4 46 F0 D2 EB D2
&H0EEC4C34  34 4C EC 7E D5 7C 1C F2 74 D9 16 8A E9 5E C4 46
&H0ABB6616  16 66 BB 1A 34 4C EC 7E D5 7C 1C F2 74 D9 16 8A
&H065AC14C  4C C1 5A 26 16 66 BB 1A 34 4C EC 7E D5 7C 1C F2
```

続けてもう一度実行すると、以下の出力が得られた。

```
----------  4C C1 5A 26 16 66 BB 1A 34 4C EC 7E D5 7C 1C F2
&H07B1EB86  86 EB B1 37 4C C1 5A 26 16 66 BB 1A 34 4C EC 7E
&H0B208CA8  A8 8C 20 2B 86 EB B1 37 4C C1 5A 26 16 66 BB 1A
&H0A6AD40C  0C D4 6A EA A8 8C 20 2B 86 EB B1 37 4C C1 5A 26
&H0ACA38AC  AC 38 CA 1A 0C D4 6A EA A8 8C 20 2B 86 EB B1 37
&H029F0DA8  A8 0D 9F A2 AC 38 CA 1A 0C D4 6A EA A8 8C 20 2B
&H0DF5909F  9F 90 F5 8D A8 0D 9F A2 AC 38 CA 1A 0C D4 6A EA
&H0183FF99  99 FF 83 31 9F 90 F5 8D A8 0D 9F A2 AC 38 CA 1A
&H0AC7AE5D  5D AE C7 7A 99 FF 83 31 9F 90 F5 8D A8 0D 9F A2
&H006F1EE0  E0 1E 6F 20 5D AE C7 7A 99 FF 83 31 9F 90 F5 8D
&H013F031A  1A 03 3F 01 E0 1E 6F 20 5D AE C7 7A 99 FF 83 31
&H0F6E4B93  93 4B 6E 2F 1A 03 3F 01 E0 1E 6F 20 5D AE C7 7A
&H089CBD65  65 BD 9C 68 93 4B 6E 2F 1A 03 3F 01 E0 1E 6F 20
&H005C3688  88 36 5C 30 65 BD 9C 68 93 4B 6E 2F 1A 03 3F 01
&H0982C44A  4A C4 82 C9 88 36 5C 30 65 BD 9C 68 93 4B 6E 2F
&H04ED3C3A  3A 3C ED 94 4A C4 82 C9 88 36 5C 30 65 BD 9C 68
&H0917CC57  57 CC 17 19 3A 3C ED 94 4A C4 82 C9 88 36 5C 30
&H082E518B  8B 51 2E C8 57 CC 17 19 3A 3C ED 94 4A C4 82 C9
&H07517C50  50 7C 51 17 8B 51 2E C8 57 CC 17 19 3A 3C ED 94
&H0AA09E6C  6C 9E A0 EA 50 7C 51 17 8B 51 2E C8 57 CC 17 19
```

ここから、以下のことがわかる。

* RAM先頭のデータの最初の4バイトが、ほぼそのまま乱数の値として出てきている。
  (最大値として `&H10000000` を指定しているので十六進数の最上位桁は `0` になっているが、それ以外は一致している)
* 乱数を生成するごとに、RAM先頭のデータの最初の12バイトが、4バイトずらしてコピーされている。

この「状態の一部をそのまま乱数として出力する」「乱数生成で状態の大部分をコピーする」というのは、[Xorshift](https://ja.wikipedia.org/wiki/Xorshift) にみられる特徴である。
さらに、RAM先頭のデータの初期値を4バイトずつリトルエンディアンで解釈すると、以下の値になっている。

```
0x05491333  0x1F123BB5  0x159A55E5  0x075BCD15
```

さらにこれらの値を十進数に変換すると、以下の値になっている。

```
88675123  521288629  362436069  123456789
```

これらの値は、論文 [Xorshift RNGs](https://www.jstatsoft.org/v08/i14/paper) に掲載されている `xor128` の状態の初期値と一致し、左からそれぞれ `w`、`z`、`y`、`x` に対応している。

そこで、乱数の生成に Xorshift が使用されていると仮定し、以下のC言語のプログラムで乱数生成を試した。

```c
#include <stdio.h>

static unsigned int x = 123456789;
static unsigned int y = 362436069;
static unsigned int z = 521288629;
static unsigned int w = 88675123;

unsigned int randint(void) {
	unsigned int t;
	t = (x ^ (x << 11));
	x = y; y = z; z = w;
	w = (w ^ (w >> 19)) ^ (t ^ (t >> 8));
	return w;
}

int main(void) {
	int i;
	for (i = 1; i < 20; i++) {
		printf("%08X\n", randint() % 0x10000000);
	}
	return 0;
}
```

このプログラムを[実行する](https://ideone.com/gbigtw)と、以下の出力が得られた。

```
0CA345EA
0B5116E6
051049AA
088D00B0
0EC7825E
0DB24146
0AF81443
0AC00F2C
0837AD58
07906569
0D9031D4
0703EE25
02EBD2F0
06C45EE9
0A16D974
021C7CD5
0EEC4C34
0ABB6616
065AC14C
```

ORANGE pico の出力と、この範囲ではピッタリ一致している。
したがって、ORANGE pico の乱数生成には Xorshift が使用されていると推測できる。

また、ORANGE pico をリセットして先ほどのプログラムを再び実行すると、再び同じ結果が得られた。
もう一度リセットして実行すると、また同じ結果が得られた。
よって、ORANGE pico で生成される乱数は、RAMに格納されている状態を直接書き換えなければリセットしてから `rnd()` を呼び出した回数 (と `rnd()` に渡す最大値) によって決まり、同じ回数・同じ最大値であれば常に同じ値になると推測できる。
(乱数を初期化 (系列を設定) するコマンドや関数は、今のところ見つかっていない)

# 中途半端な最大値と乱数

最大値 `&H10000000` を指定した場合は、ORANGE pico の乱数と、C言語の Xorshift で `randint` の戻り値を単純に最大値で割った剰余で生成した乱数が一致するようだった。
この最大値は、`randint` が返す値として考えられる値の種類数 `0x100000000` を割り切る。
しかし、種類数を割り切れない最大値を指定するとどうなるだろうか？
この場合、単純に剰余をとってしまうと、同じ乱数の値に変換される生の `randint` の戻り地の数が変わってしまい、乱数が一様分布でなくなってしまう。

`0x100000000` (`4294967296`) を割った余りがなるべく大きくなる値を探すと、`0x56000000` が見つかった。(余りは種類数の約1/3の `1409286144`)
そこで、以下のプログラムで、これを最大値として乱数を生成してみた。
最大値を変更しただけでなく、乱数を十六進数で出力する意味が薄くなるので十進数での出力に変更した。

```
10 memarrayaddr=&H03C8:dim array(3)=[&HA436,&H63DB,&HF632,&H4AE6]
20 for i=0 to 3
30 mpoke i*2+0,array(i)
40 mpoke i*2+1,array(i)>>8
50 next
60 offset=0
70 offset=offset-1
80 if array(offset)<>array(0) then goto 70
90 if array(offset+1)<>array(1) then goto 70
100 if array(offset+2)<>array(2) || array(offset+3)<>array(3) then goto 70
110 offset=offset-memarrayaddr/2
120 uartput 1,"---------- "
130 for i=0 to 7
140 uartput 1,format$(" %02X %02X",array(offset+i) & &HFF, (array(offset+i)>>8) & &HFF)
150 next
160 uartput 1,chr$(&H0A)
170 for c=2 to 20
180 uartput 1,format$("%10d ",rnd(&H56000000))
190 for i=0 to 7
200 uartput 1,format$(" %02X %02X",array(offset+i) & &HFF, (array(offset+i)>>8) & &HFF)
210 next
220 uartput 1,chr$(&H0A)
230 next
```

このプログラムを続けて2回実行すると、以下の結果が得られた。

```
----------  33 13 49 05 B5 3B 12 1F E5 55 9A 15 15 CD 5B 07
 816006634  EA 45 A3 DC 33 13 49 05 B5 3B 12 1F E5 55 9A 15
 458299110  E6 16 51 1B EA 45 A3 DC 33 13 49 05 B5 3B 12 1F
1058032042  AA 49 10 95 E6 16 51 1B EA 45 A3 DC 33 13 49 05
 747438256  B0 00 8D D8 AA 49 10 95 E6 16 51 1B EA 45 A3 DC
 516391518  5E 82 C7 1E B0 00 8D D8 AA 49 10 95 E6 16 51 1B
 934428998  46 41 B2 8D 5E 82 C7 1E B0 00 8D D8 AA 49 10 95
1157108803  43 14 F8 9A 46 41 B2 8D 5E 82 C7 1E B0 00 8D D8
 717229868  2C 0F C0 2A 43 14 F8 9A 46 41 B2 8D 5E 82 C7 1E
 137866584  58 AD 37 08 2C 0F C0 2A 43 14 F8 9A 46 41 B2 8D
 395339113  69 65 90 17 58 AD 37 08 2C 0F C0 2A 43 14 F8 9A
1301295572  D4 31 90 4D 69 65 90 17 58 AD 37 08 2C 0F C0 2A
 285470245  25 EE 03 67 D4 31 90 4D 69 65 90 17 58 AD 37 08
 652989168  F0 D2 EB D2 25 EE 03 67 D4 31 90 4D 69 65 90 17
1187274473  E9 5E C4 46 F0 D2 EB D2 25 EE 03 67 D4 31 90 4D
 873912692  74 D9 16 8A E9 5E C4 46 F0 D2 EB D2 25 EE 03 67
1176272085  D5 7C 1C F2 74 D9 16 8A E9 5E C4 46 F0 D2 EB D2
 686574644  34 4C EC 7E D5 7C 1C F2 74 D9 16 8A E9 5E C4 46
 448488982  16 66 BB 1A 34 4C EC 7E D5 7C 1C F2 74 D9 16 8A
 643481932  4C C1 5A 26 16 66 BB 1A 34 4C EC 7E D5 7C 1C F2
```

```
----------  4C C1 5A 26 16 66 BB 1A 34 4C EC 7E D5 7C 1C F2
 934407046  86 EB B1 37 4C C1 5A 26 16 66 BB 1A 34 4C EC 7E
 723553448  A8 8C 20 2B 86 EB B1 37 4C C1 5A 26 16 66 BB 1A
1047188492  0C D4 6A EA A8 8C 20 2B 86 EB B1 37 4C C1 5A 26
 449460396  AC 38 CA 1A 0C D4 6A EA A8 8C 20 2B 86 EB B1 37
1285492136  A8 0D 9F A2 AC 38 CA 1A 0C D4 6A EA A8 8C 20 2B
 938840223  9F 90 F5 8D A8 0D 9F A2 AC 38 CA 1A 0C D4 6A EA
 830734233  99 FF 83 31 9F 90 F5 8D A8 0D 9F A2 AC 38 CA 1A
 617066077  5D AE C7 7A 99 FF 83 31 9F 90 F5 8D A8 0D 9F A2
 544153312  E0 1E 6F 20 5D AE C7 7A 99 FF 83 31 9F 90 F5 8D
  20906778  1A 03 3F 01 E0 1E 6F 20 5D AE C7 7A 99 FF 83 31
 795757459  93 4B 6E 2F 1A 03 3F 01 E0 1E 6F 20 5D AE C7 7A
 312261989  65 BD 9C 68 93 4B 6E 2F 1A 03 3F 01 E0 1E 6F 20
 811349640  88 36 5C 30 65 BD 9C 68 93 4B 6E 2F 1A 03 3F 01
 495109194  4A C4 82 C9 88 36 5C 30 65 BD 9C 68 93 4B 6E 2F
1055734842  3A 3C ED 94 4A C4 82 C9 88 36 5C 30 65 BD 9C 68
 420990039  57 CC 17 19 3A 3C ED 94 4A C4 82 C9 88 36 5C 30
 472797579  8B 51 2E C8 57 CC 17 19 3A 3C ED 94 4A C4 82 C9
 391216208  50 7C 51 17 8B 51 2E C8 57 CC 17 19 3A 3C ED 94
1050713708  6C 9E A0 EA 50 7C 51 17 8B 51 2E C8 57 CC 17 19
```

乱数の最大値として `&H10000000` を用いたプログラムの出力と、この範囲の乱数の状態は一致しているようである。
乱数の偏りをなくすために悪い値が出たら再生成をする処理をすると、`rnd` 関数の呼び出し1回で行う乱数生成の回数が変わり、乱数の状態がズレることが予想される。
しかし、乱数の状態のズレはこの範囲では観測されなかった。

さらに、以下のC言語のプログラムで、単純な剰余で乱数を生成してみる。

```c
#include <stdio.h>

static unsigned int x = 123456789;
static unsigned int y = 362436069;
static unsigned int z = 521288629;
static unsigned int w = 88675123;

unsigned int randint(void) {
	unsigned int t;
	t = (x ^ (x << 11));
	x = y; y = z; z = w;
	w = (w ^ (w >> 19)) ^ (t ^ (t >> 8));
	return w;
}

int main(void) {
	int i;
	for (i = 0; i < 38; i++) {
		printf("%10u\n", randint() % 0x56000000);
	}
	return 0;
}
```

以下の[実行結果](https://ideone.com/HjTJAn)が得られた。

```
 816006634
 458299110
1058032042
 747438256
 516391518
 934428998
1157108803
 717229868
 137866584
 395339113
1301295572
 285470245
 652989168
1187274473
 873912692
1176272085
 686574644
 448488982
 643481932
 934407046
 723553448
1047188492
 449460396
1285492136
 938840223
 830734233
 617066077
 544153312
  20906778
 795757459
 312261989
 811349640
 495109194
1055734842
 420990039
 472797579
 391216208
1050713708
```

ORANGE pico により生成した乱数と、この範囲ではピッタリ一致している。
よって、ORANGE pico は、乱数の偏りをなくす工夫はしない単純な剰余で乱数を生成していると推測できる。

# さらに多くの乱数を生成し、分布を確認する

再び最大値 `&H56000000` を用い、乱数を生成する。
さらに、生成した乱数を16個のブロックに分け、それぞれのブロックに入った乱数を数える。
今回は `160000` 個の乱数を生成するので、等確率であれば各ブロックに `10000` 個程度の乱数が入るはずである。
なお、今回は出力が少ないため、`uartput` ではなく `print` を使用している。

```
10 max=&H56000000
20 blocknum=16
30 trynum=160000
40 dim count(blocknum-1,1)
50 for i=0 to blocknum-1
60 count(i,0)=0:count(i,1)=0
70 next
80 blocksize=max/blocknum
90 for i=1 to trynum
100 block=rnd(max) / blocksize
110 count(block,0)=count(block,0)+1
120 if count(block,0)>9999 then count(block,0)=0:count(block,1)=count(block,1)+1
120 next
130 for i=0 to blocknum-1
140 print format$("%2d ",i);
150 if count(i,1)>0 then print format$("%4d%04d",count(i,1),count(i,0)) else print format$("%8d",count(i,0))
160 next
```

リセットしてから実行すると、以下の出力が得られた。

```
 0    10054
 1    10116
 2    10013
 3     9974
 4    10115
 5    10004
 6    10222
 7    10005
 8     9923
 9    10168
10    10011
11    10119
12    10086
13    10021
14    10217
15     8952
```

最後のブロックに入った乱数が目立って少なく、偏りの存在を示唆している。

同様の処理をC言語でも行ってみた。

```c
#include <stdio.h>

static unsigned int x = 123456789;
static unsigned int y = 362436069;
static unsigned int z = 521288629;
static unsigned int w = 88675123;

unsigned int randint(void) {
	unsigned int t;
	t = (x ^ (x << 11));
	x = y; y = z; z = w;
	w = (w ^ (w >> 19)) ^ (t ^ (t >> 8));
	return w;
}

int main(void) {
	unsigned int max = 0x56000000u;
	unsigned int blocksize = max / 16;
	int count[16] = {0};
	int i;
	for (i = 0; i < 160000; i++) {
#if 1
		unsigned int value = randint() % max;
#else
		unsigned int value;
		int ok = 0;
		unsigned int r = (0xffffffffu % max + 1) % max;
		do {
			unsigned int x = randint();
			if (r == 0 || x < 0xffffffff - (r - 1)) { ok = 1; value = x % max; }
		} while (!ok);
#endif
		count[value / blocksize]++;
	}
	for (i = 0; i < 16; i++) {
		printf("%2d %8d\n", i, count[i]);
	}
	return 0;
}
```

[実行する](https://ideone.com/sZDTZk)と、以下の出力が得られた。

```
 0    10054
 1    10116
 2    10013
 3     9974
 4    10115
 5    10004
 6    10222
 7    10005
 8     9923
 9    10168
10    10011
11    10119
12    10086
13    10021
14    10217
15     8952
```

今回も ORANGE pico と出力がピッタリ一致した。
よって、やはり ORANGE pico でも単純な剰余によって乱数を生成していると推測できる。

`#if 1` を `#if 0` に変え、`randint` の返り値が中途半端な部分だった場合は生成をやり直す処理を入れて[実行する](https://ideone.com/QnwktD)と、以下の出力が得られた。

```
 0     9842
 1    10070
 2    10015
 3     9920
 4    10081
 5     9930
 6    10118
 7    10081
 8    10010
 9    10027
10     9874
11    10075
12     9987
13     9944
14     9969
15    10057
```
すべての値が 10000±200 の範囲に収まっており、偏りが減っている。
よって、偏りを減らす計算方法はあるにもかかわらず、使用していない、といえる。

# 結論

ORANGE pico の `rnd` 関数は、Xorshift (`xor128`) を用いて32ビットの乱数を1個生成し、その値を `rnd` 関数に渡された引数で割った余りを返していると推測できる。
Xorshift の状態は ORANGE pico のリセット時に論文に載っている値に初期化され、(RAM上の状態を直接書き換えなければ) リセットから同じ回数、同じ最大値で `rnd` 関数を呼び出すと常に同じ値を返すと推測できる。