TrueType命令で擬似乱数: 線形合同法

TrueType命令で遊ぶシリーズ。

1. 初めてのTrueType命令: Windowsでは見えないフォントをつくる - にせねこメモ
2. フォントサイズに合わせて回転するフォントを作る(1) - にせねこメモ
3. PPEM・ポイントサイズを表示するフォント - にせねこメモ
4. TrueType命令で三角関数(sin, cos)を計算する - にせねこメモ
5. フォントサイズに合わせて回転するフォントを作る(2) - にせねこメモ

擬似乱数の生成

TrueType命令で擬似乱数が計算できれば夢が広がるのではないかなーと思ってやってみた。

擬似乱数は線形合同法を使うのが楽っぽい。実際に全然乱数ではないのだが、シミュレーションとかに使うのではないので、そこそこ乱数っぽく見えれば問題ないのではと思う。

線形合同法

線形合同法 - Wikipedia
Wikipediaに書いてあった、Park & Millerによる「最低基準」とされる次の式を計算する。

http://c-faq.com/lib/rand.htmlにこれを実装したCコードがあり、次のようになっている。(16807→48271に変更してある)

#define a 48271
#define m 2147483647
#define q (m / a)
#define r (m % a)

static long int seed = 1;

long int PMrand()
{
	long int hi = seed / q;
	long int lo = seed % q;
	long int test = a * lo - r * hi;
	if(test > 0)
		seed = test;
	else	seed = test + m;
	return seed;
}

ここで、

なので、

hi = seed / 44488;
lo = seed % 44488;
test = 48271 * lo - 3399 * hi;
if(test>0) seed = test;
else       seed = test+2147483647;

みたいな計算を実際にすればいいことがわかる。

実装

seedは正数とする。
seedをスタックトップにおいて呼び出すとseedから計算される次の(擬似)乱数が返ってくる関数0は次のようになった:

# func0: calculate random number (linear congruential generator)
# initial stack: ...|seed(>0)| 
# final stack:   ...|pseudo_random_number|

PUSHB_1 # func0
 0
FDEF     #seed|
DUP      #seed|seed|
PUSHW_2
 30000
 14488
ADD      #seed|seed|44488|
DUP      #seed|seed|44488|44488|
ROLL     #seed|44488|44488|seed|
SWAP     #seed|44488|seed|44488|
PUSHW_1  #seed|44488|seed|44488|64.0|
 4096
MUL      #seed|44488|seed|44488*64|
DIV      #seed|44488|seed/44488|
ROLL     #44488|seed/44488|seed|
ROLL     #seed/44488|seed|44488|
PUSHB_1
 3
CINDEX   #seed/44488|seed|44488|seed/44488|
PUSHW_1  #seed/44488|seed|44488|seed/44488|64.0|
 4096
MUL      #seed/44488|seed|44488|(seed/44488)*64|
MUL      #seed/44488|seed|44488*(seed/44488)|
SUB      #seed/44488|seed%44488| #seed%44488 = seed - 44488*(seed/44488)
PUSHW_3  #seed/44488|seed%44488|64.0|30000|18271|
 4096
 30000
 18271
ADD      #seed/44488|seed%44488|64.0|48271|
MUL      #seed/44488|seed%44488|48271*64|
MUL      #seed/44488|(seed%44488)*48271|
SWAP     #(seed%44488)*48271|seed/44488|
PUSHW_2  #(seed%44488)*48271|seed/44488|3399|64.0|
 3399
 4096
MUL      #(seed%44488)*48271|seed/44488|3399*64|
MUL      #(seed%44488)*48271|(seed/44488)*3399|
SUB      #test| # test = (seed%44488)*48271-(seed/44488)*3399
DUP      #test|test|
PUSHB_1
 0       #test|test|0|
GT       #test|test>0|
IF       #test|  #if test>0
  # do nothing here
ELSE
 PUSHW_7 #test|0xFF|256.0|0xFF|256.0|0xFF|256.0|0x7F|
  255   # 0xFF
  16384 # 256.0
  255   # 0xFF
  16384 # 256.0
  255   # 0xFF
  16384 # 256.0
  127   # 0x7F
 MUL     #test|0xFF|256.0|0xFF|256.0|0xFF|0x7F00|
 ADD     #test|0xFF|256.0|0xFF|256.0|0x7FFF|
 MUL     #test|0xFF|256.0|0xFF|0x7FFF00|
 ADD     #test|0xFF|256.0|0x7FFFFF|
 MUL     #test|0xFF|0x7FFFFF00|
 ADD     #test|0x7FFFFFFF| #2147483647 = 0x7FFFFFFF
 ADD     #test+2147483647|
EIF
ENDF

これを、例えば

PUSHB_1
 1
0 CALL

などとしてseed=1として呼び出すとスタックトップに48271が積まれる。

改良

実際に乱数として使うのであれば、seedはスタックに積まずにstorage areaに格納することにするのが良いのではないかと思う。
prepで初期seedをstorage areaに書きこみ、乱数関数ではstorage areaからseedを読み込み、次の乱数を計算して値を返す(この乱数は次のseedとなるためstorage areaに書きこまれる)。
ここでは、storage areaの0番をseedの記録に使用することにする。
このような関数1を先ほどの関数を弄って作成した。

# func0: calculate random number (linear congruential generator)
# initial stack: ...| 
# final stack:   ...|pseudo_random_number|

PUSHB_1 # func1
 1
FDEF     #|
PUSHB_1
 0   # storage area number
RS       #seed|  # seed = StorageArea[0]
DUP      #seed|seed|
PUSHW_2
 30000
 14488
ADD      #seed|seed|44488|
DUP      #seed|seed|44488|44488|
ROLL     #seed|44488|44488|seed|
SWAP     #seed|44488|seed|44488|
PUSHW_1  #seed|44488|seed|44488|64.0|
 4096
MUL      #seed|44488|seed|44488*64|
DIV      #seed|44488|seed/44488|
ROLL     #44488|seed/44488|seed|
ROLL     #seed/44488|seed|44488|
PUSHB_1
 3
CINDEX   #seed/44488|seed|44488|seed/44488|
PUSHW_1  #seed/44488|seed|44488|seed/44488|64.0|
 4096
MUL      #seed/44488|seed|44488|(seed/44488)*64|
MUL      #seed/44488|seed|44488*(seed/44488)|
SUB      #seed/44488|seed%44488| #seed%44488 = seed - 44488*(seed/44488)
PUSHW_3  #seed/44488|seed%44488|64.0|30000|18271|
 4096
 30000
 18271
ADD      #seed/44488|seed%44488|64.0|48271|
MUL      #seed/44488|seed%44488|48271*64|
MUL      #seed/44488|(seed%44488)*48271|
SWAP     #(seed%44488)*48271|seed/44488|
PUSHW_2  #(seed%44488)*48271|seed/44488|3399|64.0|
 3399
 4096
MUL      #(seed%44488)*48271|seed/44488|3399*64|
MUL      #(seed%44488)*48271|(seed/44488)*3399|
SUB      #test| # test = (seed%44488)*48271-(seed/44488)*3399
DUP      #test|test|
PUSHB_1
 0       #test|test|0|
GT       #test|test>0|
IF       #test|  #if test>0
  #do nothing here
ELSE
 PUSHW_7 #test|0xFF|256.0|0xFF|256.0|0xFF|256.0|0x7F|
  255   # 0xFF
  16384 # 256.0
  255   # 0xFF
  16384 # 256.0
  255   # 0xFF
  16384 # 256.0
  127   # 0x7F
 MUL     #test|0xFF|256.0|0xFF|256.0|0xFF|0x7F00|
 ADD     #test|0xFF|256.0|0xFF|256.0|0x7FFF|
 MUL     #test|0xFF|256.0|0xFF|0x7FFF00|
 ADD     #test|0xFF|256.0|0x7FFFFF|
 MUL     #test|0xFF|0x7FFFFF00|
 ADD     #test|0x7FFFFFFF| #2147483647 = 0x7FFFFFFF
 ADD     #test+2147483647|
EIF
DUP      #seed|seed|
PUSHB_1  #seed|seed|0|
 0
SWAP     #seed|0|seed|
WS
ENDF

これを'prep'で、適当な値で初期化する。

PUSHB_1
 0
PUSHB_1
 1  #初期seed (>0)
WS      # StorageArea[0] = 1

すると、関数1を呼び出すと、スタックに乱数が積まれていくようになるはずである。
試しに、グリフから

10 1 LOOPCALL

と呼び出してみると、スタックに乱数が積まれていく。

右側中段のスタック表示ウィンドウを見ると、上で記載したC言語のコードの出力と同じになっていることが確かめられる。

実際に使う段には、最大値が大きいので、適当な数で割ってスケールして使う必要があると思う。

コード(コピペ用)

'fpgm'

PUSHB_1
 1
FDEF
PUSHB_1
 0
RS
DUP
PUSHW_2
 30000
 14488
ADD
DUP
ROLL
SWAP
PUSHW_1
 4096
MUL
DIV
ROLL
ROLL
PUSHB_1
 3
CINDEX
PUSHW_1
 4096
MUL
MUL
SUB
PUSHW_3
 4096
 30000
 18271
ADD
MUL
MUL
SWAP
PUSHW_2
 3399
 4096
MUL
MUL
SUB
DUP
PUSHB_1
 0
GT
IF
ELSE
PUSHW_7
 255
 16384
 255
 16384
 255
 16384
 127
MUL
ADD
MUL
ADD
MUL
ADD
ADD
EIF
DUP
PUSHB_1
 0
SWAP
WS
ENDF

'prep'

PUSHB_1
 0
PUSHB_1
 1
WS

参考

• 線形合同法 - Wikipedia
• http://c-faq.com/lib/rand.html