乱数調整で遊ぼう 別解

ろいしんさんの乱数調整で遊ぼうを中身を見ながら解析します。

1. 文字列を調べる

飛ばしてもよいです。バイナリエディタでは2バイトアラインされていないところの2バイト文字が見れなかったりするので念のため。
以下のRubyスクリプトで文字列を出力します

fname = ARGV[0]
File.binread(fname).scan(/(?:[\x81-\x9f\xe0-\xfc][\x40-\x7e\x80-\xfc]|[\x20-\x7e]){4,}/n).each do |matched|
  puts matched.encode("utf-8", "cp932", undef: :replace)
end

結果

おいしいみず
サイコソーダ
ミックスオレ
何をしますか?
おかね:%d円
しらべる -> 0
ポケモン -> 1
かいもの -> 2
________
| || ⊂⊃⊂⊃ ||
| || [][][][] ||
| ||
| || [][][][] ||
| || poke
 Θ ||
| || 口口口口 ||
|| ==== ||
    ̄ ̄ ̄ ̄ ̄ ̄
自動販売機が ある!
ほしい 飲み物は....
1. おいしいみず 200円
2. サイコソーダ 350円
3. ミックスオレ 400円
0. やめる
 おかね:%d円
...やっぱり やめた
YOU LOSE...
ペラップを しばいた!
もう一度遊びますか?
はい -> 1
いいえ -> 0
おかねが 足りない....
ガコン!
%sが でてきた!
当たりだ! もう一本
%sが でてきた!
YOU WIN!!
ベトベター
オニドリル
マグマッグ
ベトベトン
がんばりや
さみしがり
ゆうかん 
いじっぱり
やんちゃ 
ずぶとい 
すなお  
のんき  
わんぱく 
のうてんき
おくびょう
せっかち 
まじめ  
ようき  
むじゃき 
ひかえめ 
おっとり 
れいせい 
てれや  
うっかりや
おだやか 
おとなしい
なまいき 
しんちょう
きまぐれ 
するどいくちばし
きんのたま
モンスターボールを 買いました
おかね:%d円
野生の%sを つかまえた!
%s Lv.%d
性格:%s
個体値:%d-%d-%d-%d-%d-%d
なんと!
つかまえた%sが %sをもっていた!
%sを 売って
%d円 手に入れた
おかね:%d円
YOU LOSE...

2. とりあえず普通に遊ぶ

遊びます

3. IDAで逆アセンブルする

IDAで逆アセンブルします。

右下に文字列一覧があります。
試しに「当たりだ! もう一本」をクリックしましょう

「DATA XREF DATA XREF: sub_4011D0+9F↑o」というのがありますね。
これを使うと、この文字列を参照しているコードに移動できます。大変べんり。

無事、乱数を使ってあたりが出るかどうか判定しているところが見れました。
疑似乱数はx_{n+1} = (x_n * 41C64E6Dh + 6073h) % 2^32のようで、seedはdword_404374に格納されているようです。

初期seed決定も調べておきましょう。
dword_404374を選んで「Jump to xref to operand」を押すと、seedを読み書きしている部分が表示されます。

見つかりました。
これを見るとtime64()の戻り値 (1970年1月1日 00:00:00からの経過秒数)の下32bitを使っているようです。

ゲームの12/8からのバージョンではペラップが実際に音を鳴らすようになりました。
これを調べてみます。
Beepという文字列があるのでBeep() APIを使って音を鳴らしているのだと推定されます。呼んでいるところを見てみましょう。

(((seed >> 16) * 8192) >> 16) * 110 / 8192 + 440

によって周波数を定めているようです。

4. seed書き込みをトレースしながら実行

seed書き込みがあったらそのことを出力しながら実行してみます。

メニュー > Debugger > Breakpoints > Breakpoint listを開きます。

このように設定します。

そしてメニュー > Debugger > Start Processを押してみましょう。プログラムが実行されます。

このとき遊んでみると、seedが更新されるたびにTrace Windowにそのことが出力されます。

初期seedは0x5a2b5249であることがわかります。Rubyで確認してみると

irb(main):026:0> (Time.utc(1970,1,1) + 0x5a2b5249).getlocal
=> 2017-12-09 12:02:33 +0900

となって確かに起動時刻が使われています。

ペラップを鳴かせると1消費したり野生のポケモンを捕まえるといくつか消費したり、といったこともトレースされます。

長くなってきたのでこの辺りで。
これ以外の情報が気になる方はぜひ自分で解析してみてください。

ペラップの鳴き声からseedを特定するツールを作った

Pokémon RNG Advent Calendar 2017の8日目の記事です。

作りました。oupoとmizdraさんの二人の共同作です。

DSのポケモンにおいてペラップの鳴き声からseedを特定するためのツールです。要マイク。

使い方

1. ペラップの「おしゃべり」機能で880hzのサイン波を録音しておく (youtubeで880hzと検索すれば出ます)
2. このツールを開いた状態で、ペラップの鳴き声を鳴らす。するとツールの入力ボックスに自動的に入力される
3.検索したいモードなどを設定して検索を実行する
4. seed または消費数が見つかる

注意

マイクが遠くてペラップの鳴き声が拾われない、あるいは逆にペラップ以外のBGMの音も拾ってしまうという場合。このときはページ下部のmaxDecibelsを調整してください。-99から-1までの値を設定できます
鳴き声を取りこぼしてしまった場合は入力に「?」の行を入れておいてください。その行は任意の周波数にマッチするものとして処理されます。

開発

githubリポジトリを置いています。

プルリクエスト、イシュー投稿歓迎です。また、MITライセンスで公開しているのでご自由に。

今後の課題

  • Reactを全然使っていないのをリファクタリング
  • ペラップの音を拾う部分をもっと精度よく
  • 第5世代のパラメータを複数保存できるようにする
  • 第5世代のseed検索でキー入力に対応する

DSのROMのオーバーレイたちを展開するスクリプト

背景

DSはメインメモリが4MBしかなく、ゲームの全てのプログラムをメモリ上に載せることができません。そこで必要なときに必要なモジュールをメモリ上にロードしています(たとえば戦闘に入ったら戦闘用のプログラムモジュールをロードするなど)。
この各モジュールのことをオーバーレイと呼びます。

ROMファイルには各オーバーレイに対する
(ファイルID, 展開先のメモリアドレス)
の情報が入っています。

また、ファイルIDに対応する各ファイルのROM内のオフセットとサイズはROMのFAT領域に格納されています。

そこで各オーバーレイをファイルとして取り出し、圧縮を解凍し、逆アセンブルするスクリプトを書きました。

利点

ndsdis2 (http://i486.mods.jp/old/nds/ndshack.html)はメインメモリを全て逆アセンブルする仕組みでした。
これではプログラム以外が格納されているメモリも逆アセンブルするので無駄があります。
またメモリダンプしたそのときに載っているオーバーレイしか逆アセンブルできないという欠点があります。
今回のスクリプトはすべてのプログラムデータを無駄なく逆アセンブルできます。

必要なもの

スクリプト

nitrodis.rb

# NDSファイルの中のoverlayたちを展開して逆アセンブルする

# Usage: ruby nitrodis.rb <filename>.nds


require "fileutils"

OBJDUMP = 'C:/devkitPro/devkitARM/bin/arm-none-eabi-objdump.exe'
SIZEOF_OVR = 32

Fat = Struct.new(:start, :end)
Ovr = Struct.new(:id, :ramaddr, :ramsize, :bsssize, :start, :end, :file_id, :reserved)

def read_u32(f, ofs)
  f.seek ofs
  f.read(4).unpack("V")[0]
end

def read_u32_list(f, ofs, len)
  f.seek ofs
  f.read(4 * len).unpack("V*")
end

def read_cstr(f, ofs, len)
  f.seek ofs
  f.read(len).unpack("Z*")[0]
end

def read(f, ofs, bytes)
  f.seek ofs
  f.read(bytes)
end

def do_file(fname, bin, addr)
  File.binwrite(fname+".orig", bin)
  system "./DSDecmp", "-d", "-f", "lzovl", fname+".orig", fname
  destfname = File.exist?(fname) ? fname : fname + ".orig"
  system OBJDUMP, "-b", "binary", "-m", "arm",
         "-Mforce-thumb", "--adjust-vma=#{addr}",
         "-D", destfname, {1 => destfname + ".dis"}
end

fname = ARGV[0]
open(fname, "rb") do |f|
  romname = read_cstr(f, 0, 12).gsub(" ", "")
  arm9exe_start = read_u32(f, 0x20)
  arm9exe_size = read_u32(f, 0x2c)

  fat_off = read_u32(f, 0x48)
  fat_size = read_u32(f, 0x4c)
  num_files = fat_size / 8
  arm9_overlay_off = read_u32(f, 0x50)
  arm9_overlay_size = read_u32(f, 0x54)
  num_overlay_9 = arm9_overlay_size / SIZEOF_OVR
  fats = []
  ovrs = []
  num_files.times do |i|
    s, e = read_u32_list(f, fat_off + i * 8, 2)
    fats << Fat.new(s, e)
  end
  dirname = "overlay-#{romname}"
  #FileUtils.rm_r dirname
  FileUtils.mkdir_p dirname

  bin = read(f, arm9exe_start, arm9exe_size)
  do_file "#{dirname}/arm9exe", bin, 0x02000000

  num_overlay_9.times do |i|
    ovr = Ovr.new(*read_u32_list(f, arm9_overlay_off + i * SIZEOF_OVR, 8))
    ovrs << ovr
    fat = fats[ovr.file_id]
    #puts "%d %08x %d %d %08x %08x %d %d %08x" % [ovr.file_id, ovr.ramaddr, ovr.ramsize, ovr.bsssize, ovr.start, ovr.end, fat.start, fat.end - fat.start, ovr.reserved]
    #puts "%d %08x %d" % [ovr.file_id, ovr.ramaddr, ovr.ramsize]
    bin = read(f, fat.start, fat.end - fat.start)
    ovrfname = "#{dirname}/#{ovr.file_id}"
    do_file ovrfname, bin, ovr.ramaddr
  end
end

DeSmuMEの3Dレンダラの行列をLuaからいじれるようにするhack

(2015年にやったものを発掘して記事にしているのでもう細部を忘れています)

ゲーム内の処理とか一切触らずにこんな風に斜めから映したり、街の全体像を映したりできます。



(建物の一部が映ってなかったり、キャラクターの位置がずれているのはすでに修正済みのバグだった気がします)

EmuLuaにemu.extramatrixmul(), emu.extramatrixproj()という関数を用意しています。
それを使ってたとえばこんな風に書きます。

function translate(x, y, z)
	emu.extramatrixmul({1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, x, y, z, 1})
end

function rotateY(theta)
	local c = math.cos(theta)
	local s = math.sin(theta)
	emu.extramatrixmul({c, 0, s, 0, 0, 1, 0, 0,-s, 0, c, 0, 0, 0, 0, 1})
end

function rotateX(theta)
	local c = math.cos(theta)
	local s = math.sin(theta)
	emu.extramatrixmul({1, 0, 0, 0, 0, c,-s, 0, 0, s, c, 0, 0, 0, 0, 1})
end

function scale(k)
	emu.extramatrixmul({k, 0, 0, 0, 0, k, 0, 0, 0, 0, k, 0, 0, 0, 0, 1})
end

function proj(scalar, x, y, z)
	emu.extramatrixproj({scalar, 0, 0, 0, 0, scalar, 0, 0, 0, 0, 1, 0, x, y, z, 1})
end

emu.registerbefore(function ()
	emu.extramatrixinit()
	translate(-200, 0, 0)
	rotateY(0.5)
	scale(0.6)
	proj(1/2,0,0, 0)
end)

DeSmuMEのソースコードのdiffです。
元のなったDeSmuMEのリビジョンは忘れました(

diff -ru trunk-copy/desmume/src/gfx3d.cpp 3d-hack/desmume/src/gfx3d.cpp
--- trunk-copy/desmume/src/gfx3d.cpp	2014-09-02 11:36:48.539893300 +0900
+++ 3d-hack/desmume/src/gfx3d.cpp	2015-12-17 15:36:03.263874500 +0900
@@ -54,6 +54,7 @@
 #include "GPU_OSD.h"
 #endif
 
+//#define DEBUG_MATRIX
 
 /*
 thoughts on flush timing:
@@ -320,6 +321,10 @@
 static CACHE_ALIGN s32		mtxCurrent [4][16];
 static CACHE_ALIGN s32		mtxTemporal[16];
 static u32 mode = 0;
+s32 extraMatrix[16];
+s32 extraMatrixProj[16];
+
+static const char * const mode_name[] = {"projection", "coordinate", "directional", "texture"};
 
 // Indexes for matrix loading/multiplication
 static u8 ML4x4ind = 0;
@@ -534,6 +539,8 @@
 	MatrixInit (mtxCurrent[2]);
 	MatrixInit (mtxCurrent[3]);
 	MatrixInit (mtxTemporal);
+	MatrixInit (extraMatrix);
+	MatrixInit (extraMatrixProj);
 
 	MatrixStackInit(&mtxStack[0]);
 	MatrixStackInit(&mtxStack[1]);
@@ -560,7 +567,7 @@
 
 	memset(gfx3d_convertedScreen,0,sizeof(gfx3d_convertedScreen));
 
-	gfx3d.state.clearDepth = DS_DEPTH15TO24(0x7FFF);
+	gfx3d.state.clearDepth = 0xffffffffu;
 	
 	clInd2 = 0;
 	isSwapBuffers = FALSE;
@@ -587,6 +594,8 @@
 	return fx32_shiftdown(fx32_mul(a[0],b[0]) + fx32_mul(a[1],b[1]) + fx32_mul(a[2],b[2]));
 }
 
+static void printMatrix4x4(s32 *m);
+
 #define SUBMITVERTEX(ii, nn) polylist->list[polylist->count].vertIndexes[ii] = tempVertInfo.map[nn];
 //Submit a vertex to the GE
 static void SetVertex()
@@ -623,8 +632,11 @@
 	//so that we only have to multiply one matrix here
 	//(we could lazy cache the concatenated clip matrix and only generate it
 	//when we need to)
-	MatrixMultVec4x4_M2(mtxCurrent[0], coordTransformed);
-
+	MatrixMultVec4x4(mtxCurrent[1],coordTransformed);
+	MatrixMultVec4x4(extraMatrix,coordTransformed);
+	MatrixMultVec4x4(mtxCurrent[0],coordTransformed);
+	MatrixMultVec4x4(extraMatrixProj,coordTransformed);
+	
 	//printf("%f %f %f\n",s16coord[0]/4096.0f,s16coord[1]/4096.0f,s16coord[2]/4096.0f);
 	//printf("x %f %f %f %f\n",mtxCurrent[0][0]/4096.0f,mtxCurrent[0][1]/4096.0f,mtxCurrent[0][2]/4096.0f,mtxCurrent[0][3]/4096.0f);
 	//printf(" = %f %f %f %f\n",coordTransformed[0]/4096.0f,coordTransformed[1]/4096.0f,coordTransformed[2]/4096.0f,coordTransformed[3]/4096.0f);
@@ -645,6 +657,18 @@
 		printf("wtf\n");
 	}
 	VERT &vert = vertlist->list[vertIndex];
+	
+	//coordTransformed[0] *= 0.5;
+	//coordTransformed[1] *= 0.5;
+#ifdef DEBUG_MATRIX
+	printf("projection matrix = ");
+	printMatrix4x4(mtxCurrent[0]);
+	puts("");
+	printf("coordinate matrix = ");
+	printMatrix4x4(mtxCurrent[1]);
+	puts("");
+	printf("setVertex: %d (%.1f,%.1f,%.1f) -> (%.1f,%.1f,%.1f,%.1f)\n", vertIndex, coord[0] / 4096.0f, coord[1] / 4096.0f, coord[2] / 4096.0f, coordTransformed[0] / 4096.0f, coordTransformed[1] / 4096.0f, coordTransformed[2] / 4096.0f, coordTransformed[3] / 4096.0f);
+#endif
 
 	//printf("%f %f %f\n",coordTransformed[0],coordTransformed[1],coordTransformed[2]);
 	//if(coordTransformed[1] > 20) 
@@ -839,6 +863,17 @@
 
 
 //===============================================================================
+static void printMatrix4x4(s32 *m)
+{
+	printf("[[%.1f,%.1f,%.1f,%.1f],[%.1f,%.1f,%.1f,%.1f],[%.1f,%.1f,%.1f,%.1f],[%.1f,%.1f,%.1f,%.1f]]", m[0] / 4096.0, m[1] / 4096.0, m[2] / 4096.0, m[3] / 4096.0, m[4] / 4096.0, m[5] / 4096.0, m[6] / 4096.0, m[7] / 4096.0, m[8] / 4096.0, m[9] / 4096.0, m[10] / 4096.0, m[11] / 4096.0, m[12] / 4096.0, m[13] / 4096.0, m[14] / 4096.0, m[15] / 4096.0);
+}
+
+static void printMatrix4x3(s32 *m)
+{
+	printf("[[%.1f,%.1f,%.1f],[%.1f,%.1f,%.1f],[%.1f,%.1f,%.1f],[%.1f,%.1f,%.1f]]", m[0] / 4096.0, m[1] / 4096.0, m[2] / 4096.0, m[4] / 4096.0, m[5] / 4096.0, m[6] / 4096.0, m[8] / 4096.0, m[9] / 4096.0, m[10] / 4096.0, m[12] / 4096.0, m[13] / 4096.0, m[14] / 4096.0);
+}
+
+
 static void gfx3d_glMatrixMode(u32 v)
 {
 	mode = (v&3);
@@ -850,13 +885,22 @@
 {
 	//this command always works on both pos and vector when either pos or pos-vector are the current mtx mode
 	short mymode = (mode==1?2:mode);
-
+	
+#ifdef DEBUG_MATRIX
+	printf("pushMatrix: %s (cur = ", mode_name[mymode]);
+	printMatrix4x4(mtxCurrent[mymode]);
+	printf(")\n");
 	MatrixStackPushMatrix(&mtxStack[mymode], mtxCurrent[mymode]);
+#endif
 
 	GFX_DELAY(17);
 
-	if(mymode==2)
+	if(mymode==2) {
+#ifdef DEBUG_MATRIX
+		printf("pushMatrix: %s\n", mode_name[1]);
+#endif
 		MatrixStackPushMatrix (&mtxStack[1], mtxCurrent[1]);
+	}
 }
 
 static void gfx3d_glPopMatrix(s32 i)
@@ -877,11 +921,18 @@
 	//please note that our ability to skip treating this as signed is dependent on the modular addressing later. if that ever changes, we need to change this back.
 
 	MatrixStackPopMatrix(mtxCurrent[mymode], &mtxStack[mymode], i);
+#ifdef DEBUG_MATRIX
+	printf("popMatrix: %s %d\n", mode_name[mymode], i);
+#endif
 
 	GFX_DELAY(36);
 
-	if (mymode == 2)
+	if (mymode == 2) {
+#ifdef DEBUG_MATRIX
+		printf("popMatrix: %s %d\n", mode_name[1], i);
+#endif
 		MatrixStackPopMatrix(mtxCurrent[1], &mtxStack[1], i);
+	}
 }
 
 static void gfx3d_glStoreMatrix(u32 v)
@@ -903,12 +954,23 @@
 	if(v==31)
 		MMU_new.gxstat.se = 1;
 
+#ifdef DEBUG_MATRIX
+	printf("storeMatrix: %s %d (mtx = ", mode_name[mymode], v);
+	printMatrix4x4(mtxCurrent[mymode]);
+	printf(")\n");
+#endif
 	MatrixStackLoadMatrix (&mtxStack[mymode], v, mtxCurrent[mymode]);
 
 	GFX_DELAY(17);
 
-	if(mymode==2)
+	if(mymode==2) {
+#ifdef DEBUG_MATRIX
+		printf("storeMatrix: %s %d (mtx = ", mode_name[1], v);
+		printMatrix4x4(mtxCurrent[1]);
+		printf(")\n");
+#endif
 		MatrixStackLoadMatrix (&mtxStack[1], v, mtxCurrent[1]);
+	}
 }
 
 static void gfx3d_glRestoreMatrix(u32 v)
@@ -930,23 +992,36 @@
 	if(v==31)
 		MMU_new.gxstat.se = 1;
 
-
+#ifdef DEBUG_MATRIX
+	printf("restoreMatrix: %s %d\n", mode_name[mymode], v);
+#endif
 	MatrixCopy (mtxCurrent[mymode], MatrixStackGetPos(&mtxStack[mymode], v));
 
 	GFX_DELAY(36);
 
-	if (mymode == 2)
+	if (mymode == 2) {
+#ifdef DEBUG_MATRIX
+		printf("restoreMatrix: %s %d\n", mode_name[1], v);
+#endif
 		MatrixCopy (mtxCurrent[1], MatrixStackGetPos(&mtxStack[1], v));
+	}
 }
 
 static void gfx3d_glLoadIdentity()
 {
+#ifdef DEBUG_MATRIX
+	printf("loadIdentity: %s\n", mode_name[mode]);
+#endif
 	MatrixIdentity (mtxCurrent[mode]);
 
 	GFX_DELAY(19);
 
-	if (mode == 2)
+	if (mode == 2) {
+#ifdef DEBUG_MATRIX
+		printf("loadIdentity: %s\n", mode_name[1]);
+#endif
 		MatrixIdentity (mtxCurrent[1]);
+	}
 
 	//printf("identity: %d to: \n",mode); MatrixPrint(mtxCurrent[1]);
 }
@@ -959,13 +1034,22 @@
 	if(ML4x4ind<16) return FALSE;
 	ML4x4ind = 0;
 
+#ifdef DEBUG_MATRIX
+	printf("loadMatrix4x4: %s ", mode_name[mode]);
+	printMatrix4x4(mtxCurrent[mode]);
+	puts("");
+#endif
+
 	GFX_DELAY(19);
 
 	//vector_fix2float<4>(mtxCurrent[mode], 4096.f);
 
-	if (mode == 2)
+	if (mode == 2) {
+#ifdef DEBUG_MATRIX
+		printf("copy from %s to %s\n", mode_name[2], mode_name[1]);
+#endif
 		MatrixCopy (mtxCurrent[1], mtxCurrent[2]);
-
+	}
 	//printf("load4x4: matrix %d to: \n",mode); MatrixPrint(mtxCurrent[1]);
 	return TRUE;
 }
@@ -985,10 +1069,19 @@
 	mtxCurrent[mode][3] = mtxCurrent[mode][7] = mtxCurrent[mode][11] = 0;
 	mtxCurrent[mode][15] = (1<<12);
 
+#ifdef DEBUG_MATRIX
+	printf("loadMatrix4x3: %s ", mode_name[mode]);
+	printMatrix4x3(mtxCurrent[mode]);
+	puts("");
+#endif
 	GFX_DELAY(30);
 
-	if (mode == 2)
+	if (mode == 2) {
+#ifdef DEBUG_MATRIX
+		printf("copy from %s to %s\n", mode_name[2], mode_name[1]);
+#endif
 		MatrixCopy (mtxCurrent[1], mtxCurrent[2]);
+	}
 	//printf("load4x3: matrix %d to: \n",mode); MatrixPrint(mtxCurrent[1]);
 	return TRUE;
 }
@@ -1005,10 +1098,18 @@
 
 	//vector_fix2float<4>(mtxTemporal, 4096.f);
 
+#ifdef DEBUG_MATRIX
+	printf("multMatrix4x4: %s ", mode_name[mode]);
+	printMatrix4x4(mtxTemporal);
+	puts("");
+#endif
 	MatrixMultiply (mtxCurrent[mode], mtxTemporal);
 
 	if (mode == 2)
 	{
+#ifdef DEBUG_MATRIX
+		printf("copy from %s to %s\n", mode_name[2], mode_name[1]);
+#endif
 		MatrixMultiply (mtxCurrent[1], mtxTemporal);
 		GFX_DELAY_M2(30);
 	}
@@ -1036,10 +1137,18 @@
 	mtxTemporal[3] = mtxTemporal[7] = mtxTemporal[11] = 0;
 	mtxTemporal[15] = 1<<12;
 
+#ifdef DEBUG_MATRIX
+	printf("multMatrix4x3: %s ", mode_name[mode]);
+	printMatrix4x3(mtxTemporal);
+	puts("");
+#endif
 	MatrixMultiply (mtxCurrent[mode], mtxTemporal);
 
 	if (mode == 2)
 	{
+#ifdef DEBUG_MATRIX
+		printf("copy from %s to %s\n", mode_name[2], mode_name[1]);
+#endif
 		MatrixMultiply (mtxCurrent[1], mtxTemporal);
 		GFX_DELAY_M2(30);
 	}
@@ -1070,10 +1179,16 @@
 	mtxTemporal[15] = 1<<12;
 	mtxTemporal[12] = mtxTemporal[13] = mtxTemporal[14] = 0;
 
+#ifdef DEBUG_MATRIX
+	printf("multMatrix3x3: %s [[%.1f,%.1f,%.1f],[%.1f,%.1f,%.1f],[%.1f,%.1f,%.1f]]\n", mode_name[mode], mtxCurrent[mode][0] / 4096.0, mtxCurrent[mode][1] / 4096.0, mtxCurrent[mode][2] / 4096.0, mtxCurrent[mode][4] / 4096.0, mtxCurrent[mode][5] / 4096.0, mtxCurrent[mode][6] / 4096.0, mtxCurrent[mode][8] / 4096.0, mtxCurrent[mode][9] / 4096.0, mtxCurrent[mode][10] / 4096.0);
+#endif
 	MatrixMultiply (mtxCurrent[mode], mtxTemporal);
 
 	if (mode == 2)
 	{
+#ifdef DEBUG_MATRIX
+		printf("copy from %s to %s\n", mode_name[2], mode_name[1]);
+#endif
 		MatrixMultiply (mtxCurrent[1], mtxTemporal);
 		GFX_DELAY_M2(30);
 	}
@@ -1095,6 +1210,9 @@
 	if(scaleind<3) return FALSE;
 	scaleind = 0;
 
+#ifdef DEBUG_MATRIX
+	printf("scaleMatrix: %s [%.1f,%.1f,%.1f]\n", mode_name[(mode==2?1:mode)], scale[0] / 4096.0, scale[1] / 4096.0, scale[2] / 4096.0);
+#endif
 	MatrixScale (mtxCurrent[(mode==2?1:mode)], scale);
 	//printf("scale: matrix %d to: \n",mode); MatrixPrint(mtxCurrent[1]);
 
@@ -1117,12 +1235,18 @@
 	if(transind<3) return FALSE;
 	transind = 0;
 
+#ifdef DEBUG_MATRIX
+	printf("translateMatrix: %s [%.1f,%.1f,%.1f]\n", mode_name[mode], trans[0] / 4096.0, trans[1] / 4096.0, trans[2] / 4096.0);
+#endif
 	MatrixTranslate (mtxCurrent[mode], trans);
 
 	GFX_DELAY(22);
 
 	if (mode == 2)
 	{
+#ifdef DEBUG_MATRIX
+		printf("translateMatrix: %s [%.1f,%.1f,%.1f]\n", mode_name[1], trans[0] / 4096.0, trans[1] / 4096.0, trans[2] / 4096.0);
+#endif
 		MatrixTranslate (mtxCurrent[1], trans);
 		GFX_DELAY_M2(30);
 	}
@@ -1538,11 +1662,15 @@
 		//but change it all to floating point and do it that way instead
 		CACHE_ALIGN float temp1[16] = {mtxCurrent[1][0]/4096.0f,mtxCurrent[1][1]/4096.0f,mtxCurrent[1][2]/4096.0f,mtxCurrent[1][3]/4096.0f,mtxCurrent[1][4]/4096.0f,mtxCurrent[1][5]/4096.0f,mtxCurrent[1][6]/4096.0f,mtxCurrent[1][7]/4096.0f,mtxCurrent[1][8]/4096.0f,mtxCurrent[1][9]/4096.0f,mtxCurrent[1][10]/4096.0f,mtxCurrent[1][11]/4096.0f,mtxCurrent[1][12]/4096.0f,mtxCurrent[1][13]/4096.0f,mtxCurrent[1][14]/4096.0f,mtxCurrent[1][15]/4096.0f};
 		CACHE_ALIGN float temp0[16] = {mtxCurrent[0][0]/4096.0f,mtxCurrent[0][1]/4096.0f,mtxCurrent[0][2]/4096.0f,mtxCurrent[0][3]/4096.0f,mtxCurrent[0][4]/4096.0f,mtxCurrent[0][5]/4096.0f,mtxCurrent[0][6]/4096.0f,mtxCurrent[0][7]/4096.0f,mtxCurrent[0][8]/4096.0f,mtxCurrent[0][9]/4096.0f,mtxCurrent[0][10]/4096.0f,mtxCurrent[0][11]/4096.0f,mtxCurrent[0][12]/4096.0f,mtxCurrent[0][13]/4096.0f,mtxCurrent[0][14]/4096.0f,mtxCurrent[0][15]/4096.0f};
+		CACHE_ALIGN float extra[16] = {extraMatrix[0]/4096.0f,extraMatrix[1]/4096.0f,extraMatrix[2]/4096.0f,extraMatrix[3]/4096.0f,extraMatrix[4]/4096.0f,extraMatrix[5]/4096.0f,extraMatrix[6]/4096.0f,extraMatrix[7]/4096.0f,extraMatrix[8]/4096.0f,extraMatrix[9]/4096.0f,extraMatrix[10]/4096.0f,extraMatrix[11]/4096.0f,extraMatrix[12]/4096.0f,extraMatrix[13]/4096.0f,extraMatrix[14]/4096.0f,extraMatrix[15]/4096.0f};
+		CACHE_ALIGN float extra2[16] = {extraMatrixProj[0]/4096.0f,extraMatrixProj[1]/4096.0f,extraMatrixProj[2]/4096.0f,extraMatrixProj[3]/4096.0f,extraMatrixProj[4]/4096.0f,extraMatrixProj[5]/4096.0f,extraMatrixProj[6]/4096.0f,extraMatrixProj[7]/4096.0f,extraMatrixProj[8]/4096.0f,extraMatrixProj[9]/4096.0f,extraMatrixProj[10]/4096.0f,extraMatrixProj[11]/4096.0f,extraMatrixProj[12]/4096.0f,extraMatrixProj[13]/4096.0f,extraMatrixProj[14]/4096.0f,extraMatrixProj[15]/4096.0f};
 
 		DS_ALIGN(16) VERT_POS4f vert = { verts[i].x, verts[i].y, verts[i].z, verts[i].w };
-
+		
 		_NOSSE_MatrixMultVec4x4(temp1,verts[i].coord);
+		_NOSSE_MatrixMultVec4x4(extra,verts[i].coord);
 		_NOSSE_MatrixMultVec4x4(temp0,verts[i].coord);
+		_NOSSE_MatrixMultVec4x4(extra2,verts[i].coord);
 	}
 
 	//clip each poly
@@ -1654,7 +1782,11 @@
 
 void gfx3d_glClearDepth(u32 v)
 {
-	gfx3d.state.clearDepth = DS_DEPTH15TO24(v);
+	if (v == 0x7fff) {
+		gfx3d.state.clearDepth = 0xffffffffu;
+	} else {
+		gfx3d.state.clearDepth = DS_DEPTH15TO24(v);
+	}
 }
 
 // Ignored for now
diff -ru trunk-copy/desmume/src/lua-engine.cpp 3d-hack/desmume/src/lua-engine.cpp
--- trunk-copy/desmume/src/lua-engine.cpp	2014-09-02 11:36:48.390793300 +0900
+++ 3d-hack/desmume/src/lua-engine.cpp	2015-12-17 15:35:16.475717600 +0900
@@ -3788,6 +3788,42 @@
 	return 1;
 }
 
+extern s32 extraMatrix[16];
+extern s32 extraMatrixProj[16];
+
+DEFINE_LUA_FUNCTION(emu_extramatrixinit, "")
+{
+	MatrixIdentity(extraMatrix);
+	MatrixIdentity(extraMatrixProj);
+	return 0;
+}
+
+DEFINE_LUA_FUNCTION(emu_extramatrixmul, "matrix")
+{
+	s32 tmp[16];
+	luaL_checktype(L, 1, LUA_TTABLE);
+	for (int i = 0; i < 16; i++) {
+		lua_rawgeti(L, 1, i + 1);
+		tmp[i] = (s32)(lua_tonumber(L, -1) * 4096);
+		lua_pop(L, 1);
+	}
+	MatrixMultiply(extraMatrix, tmp);
+	return 0;
+}
+
+DEFINE_LUA_FUNCTION(emu_extramatrixproj, "matrix")
+{
+	s32 tmp[16];
+	luaL_checktype(L, 1, LUA_TTABLE);
+	for (int i = 0; i < 16; i++) {
+		lua_rawgeti(L, 1, i + 1);
+		tmp[i] = (s32)(lua_tonumber(L, -1) * 4096);
+		lua_pop(L, 1);
+	}
+	MatrixMultiply(extraMatrixProj, tmp);
+	return 0;
+}
+
 // TODO
 /*
 DEFINE_LUA_FUNCTION(emu_loadrom, "filename")
@@ -4541,6 +4577,9 @@
 	{"registermenustart", emu_registermenustart},
 	// alternative names
 //	{"openrom", emu_loadrom},
+	{"extramatrixinit", emu_extramatrixinit},
+	{"extramatrixmul", emu_extramatrixmul},
+	{"extramatrixproj", emu_extramatrixproj},
 	{NULL, NULL}
 };
 static const struct luaL_reg guilib [] =
diff -ru trunk-copy/desmume/src/rasterize.cpp 3d-hack/desmume/src/rasterize.cpp
--- trunk-copy/desmume/src/rasterize.cpp	2014-09-02 11:36:37.489425300 +0900
+++ 3d-hack/desmume/src/rasterize.cpp	2015-12-17 15:01:54.567041000 +0900
@@ -869,7 +869,8 @@
 			left->Step();
 			right->Step();
 
-			if(!RENDERER && _debug_thisPoly)
+			//if(!RENDERER && _debug_thisPoly)
+			if (0)
 			{
 				//debug drawing
 				bool top = (horizontal&&first);
@@ -1034,6 +1035,7 @@
 	template<bool SLI>
 	FORCEINLINE void mainLoop(SoftRasterizerEngine* const engine)
 	{
+		INFO("mainLoop() polyCount=%d\n", engine->clippedPolyCounter);
 		this->engine = engine;
 		lastTexKey = NULL;
 
@@ -1236,7 +1238,7 @@
 	//I am not sure whether it is right, though. previously this was cleared to 0, as a guess,
 	//but in spiderman2 some fires with polyid 0 try to render on top of the background
 	clearFragment.polyid.translucent = kUnsetTranslucentPolyID; 
-	clearFragment.depth = gfx3d.renderState.clearDepth;
+	clearFragment.depth = 0xffffffffu; //gfx3d.renderState.clearDepth;
 	clearFragment.stencil = 0;
 	clearFragment.isTranslucentPoly = 0;
 	clearFragment.fogged = BIT15(gfx3d.renderState.clearColor);
@@ -1686,17 +1688,18 @@
 
 	softRastHasNewData = true;
 	
-	if (rasterizerCores > 1)
+	/*if (rasterizerCores > 1)
 	{
 		for(unsigned int i = 0; i < rasterizerCores; i++)
 		{
 			rasterizerUnitTask[i].execute(&execRasterizerUnit, (void *)i);
 		}
 	}
-	else
+	else*/
 	{
 		rasterizerUnit[0].mainLoop<false>(&mainSoftRasterizer);
 	}
+	//driver->EMU_PauseEmulation(true);
 }
 
static void SoftRastRenderFinish()

計算結果をキャッシュする&コピーコストのかからないMT

背景

乱数調整のツールを作っていると次のようなコードを書きたくなります。
しかし、このコードはshow_ivs()を呼び出すたびにMTがコピーされるので遅いです。

struct MT {
	uint32_t table[624];
	int index;
};

void init_genrand(MT *mt, uint32_t seed);
uint32_t genrand_int32(MT *mt);

void show_ivs(MT mt) {
	for (int i = 0; i < 6; i ++) {
		printf("%02d ", (int)((uint64_t)(genrand_int32(&mt) * 32) >> 32));
	}
	printf("\n");
}

void list(uint32_t seed) {
	MT mt;
	init_genrand(&mt, seed);
	for (int i = 0; i < 100; i ++) {
		show_ivs(mt);
		genrand_int32(&mt);
	}
}

本題

そこでコピーのコストがかからず、また一度計算した値はキャッシュするMTを作ってみました。

しかし、連結リストを使っているので遅い気がします。
配列に乱数値を貯めておき、それを使うという方法の方がずっといいかもしれません。

#include <memory>
#include <cstdint>
#include <cstdio>

#define N 624
#define M 397

namespace MT {

// デバッグ用
int objects_count = 0;
int calc_count = 0;

const uint32_t MATRIX_A = 0x9908b0dfUL;
const uint32_t UPPER_MASK = 0x80000000UL;
const uint32_t LOWER_MASK = 0x7fffffffUL;

class State;

class Node {
	friend State;
private:
	uint32_t m_val;
	std::shared_ptr<Node> m_next;
	std::shared_ptr<Node> m_prev_624;
	std::shared_ptr<Node> m_prev_623;
	std::shared_ptr<Node> m_prev_227;
public:
	Node(uint32_t val) : m_val(val) {
		objects_count++;
	}
	Node(std::shared_ptr<Node> prev_624, std::shared_ptr<Node> prev_623, std::shared_ptr<Node> prev_227)
		: m_prev_624(prev_624), m_prev_623(prev_623), m_prev_227(prev_227) {
		calc_value();
		objects_count++;
	}

	Node(Node &other) : m_val(other.m_val), m_next(other.m_next),
		m_prev_624(other.m_prev_624), m_prev_623(other.m_prev_623), m_prev_227(other.m_prev_227) {
		
		objects_count++;
	}

	~Node() {
		objects_count--;
	}

	std::shared_ptr<Node> next() {
		if (m_next != nullptr) {
			return m_next;
		} else {
			m_next = std::make_shared<Node>(Node(m_prev_624->next(), m_prev_623->next(), m_prev_227->next()));
			m_prev_624 = m_prev_623 = m_prev_227 = nullptr;
			return m_next;
		}
	}

private:
	void calc_value()
	{
		uint32_t y = (m_prev_624->m_val & UPPER_MASK) | (m_prev_623->m_val & LOWER_MASK);
		m_val = m_prev_227->m_val ^ (y >> 1) ^ ((y & 1) != 0 ? MATRIX_A : 0);
		calc_count++;
	}
};

class State {
	std::shared_ptr<Node> m_node;
public:
	State(uint32_t seed) {
		uint32_t vals[N];
		std::shared_ptr<Node> nodes[N];
		vals[0] = seed;
		for (int i = 1; i < N; i++)
		{
			vals[i] = 1812433253 * (vals[i - 1] ^ (vals[i - 1] >> 30)) + i;
		}
		for (int i = 0; i < N; i++)
		{
			nodes[i] = std::make_shared<Node>(Node(vals[i]));
		}
		for (int i = 1; i < N; i++)
		{
			nodes[i - 1]->m_next = nodes[i];
		}
		m_node = std::make_shared<Node>(Node(nodes[0], nodes[1], nodes[M]));
		nodes[N-1]->m_next = m_node;
	}


	uint32_t rand()
	{
		uint32_t value = m_node->m_val;
		m_node = m_node->next();
		return temper(value);
	}
private:
	uint32_t temper(uint32_t y)
	{
		y ^= (y >> 11);
		y ^= (y << 7) & 0x9d2c5680;
		y ^= (y << 15) & 0xefc60000;
		y ^= (y >> 18);
		return y;
	}

};
}

void show_ivs(MT::State mt) {
	for (int i = 0; i < 6; i++) {
		printf("%02d ", (int)(((uint64_t)mt.rand() * 32) >> 32));
	}
	printf("\n");
}

int main() {
	MT::State mt(0);
	for (int i = 0; i < 10000; i++) {
		show_ivs(mt);
		mt.rand();
	}
	printf("objects_count=%d\n", MT::objects_count);
	printf("calc_count=%d\n", MT::calc_count);
	return 0;
}

calc_count(次の乱数値を計算する回数)は10006回となって最小限におさえられています。
またobjects_count(現在メモリに存在するノードの数)は625となって不要なノードは消されていることが分かります。

問題点

MT::State mt(0);

の行の次に

MT::State mt2 = mt;

を追加するとmain()関数終了時にSEGVしてしまいます。
これはNodeのデストラクタが再帰的に呼ばれ、再帰の回数が大きくなりすぎてスタックオーバーフローするためです。

63変数連立2次方程式が解ければTinyMTのstateは連続した64個の乱数値の下位2bitから得られる

Pokémon RNG Advent Calendar 2016 13日目の記事です。

疑似乱数TinyMTにおいて得られた乱数値から乱数生成器の状態を逆算する話です。

SM孵化乱数列を計算する : ただの雑記byさきによってTinyMTのstateは連続した127個の乱数値の下位1bitがあれば得られることが示されました。

これを連続した64個の乱数値の下位2bitから得られるともっと楽できそうです。

これは63変数63本のF_2係数連立2次方程式が解ければ可能です。
たとえば、連続する64個の下位2bitが01, 00, 00, …, 00であるstateは次の連立方程式を解けば得られます。

あとは、この方程式のグレブナー基底を求めるなり、SAT-based Constraint Solverにかけるなりして解けばいいですね!
どちらも試しましたが、変数の個数の指数のオーダーの時間がかかって63変数は無理そうでした!

ソースコードは以下の通り。SageMathCloudで動かすことができます。

mat1 = 0x8f7011ee
mat2 = 0xfc78ff1f
tmat = 0x3793fdff

F2 = GF(2)

# 整数を32次元F2ベクトルに変換する
def int_to_f2(x):
    return vector([F2((x >> i) & 1) for i in range(32)])

# 整数列をF2ベクトルに変換する
def ints_to_f2(xs):
    return vector([F2((xs[floor(i / 32)] >> (i % 32)) & 1) for i in range(32*len(xs))])

# 32次元F2ベクトルを整数に変換する
def f2_to_int(vec):
    x = 0
    for i in range(32):
        x |= Integer(vec[i]) << i
    return x

# F2ベクトルを整数列に変換する
def f2_to_ints(vec):
    return [f2_to_int(vec[i*32:i*32+32]) for i in range(floor(len(vec) / 32))]

# 行列を横につなげる
def join_matrix_horizontal(mats):
    m = mats[0].nrows()
    n = mats[0].ncols()
    return matrix([[mats[floor(j / n)][i, j % n] for j in range(n*len(mats))] for i in range(m)])

# 行列を縦につなげる
def join_matrix_vertical(mats):
    m = mats[0].nrows()
    n = mats[0].ncols()
    return matrix([[mats[floor(i / m)][i % m, j] for j in range(n)] for i in range(m*len(mats))])

# 行列からそのi番目の行だけを得る
def nth_row(mat, i):
    m = mat.nrows()
    n = mat.ncols()
    return matrix([[mat[i, j] for j in range(n)]])

# 128次元ベクトルを32次元ベクトルを4つつなげたものと思ってそのi番目をとってくるという変換の行列
def elem_matrix(i):
    mats = [Mat(GF(2), 32, 32).zero() for j in range(4)]
    mats[i] = Mat(GF(2), 32, 32).identity_matrix()
    return join_matrix_horizontal(mats)

# 真なら1、偽なら0を返す
def xi(x):
    if x:
        return 1
    else:
        return 0

# 32次元ベクトルに対して右シフトを表す行列
def rshift_matrix(k):
    return matrix([[F2(xi(j - i == k)) for j in range(32)] for i in range(32)])

# 32次元ベクトルに対して左シフトを表す行列
def lshift_matrix(k):
    return rshift_matrix(-k)

# 32次元ベクトルに対して第kビットを得るという変換の行列
def get_bit_matrix(k):
    return matrix([[F2(xi(j == k)) for j in range(32)]])

get_0bit = get_bit_matrix(0)

# n次元ベクトルをn×1行列に変換する
def vect_to_matrix(vec):
    return matrix([[vec[i]] for i in range(len(vec))])

# 32次元ベクトルに対してaとandをとるという変換を表す行列
def mask_matrix(a):
    return matrix([[F2(a[i] * xi(i == j)) for j in range(32)] for i in range(32)])

# 128次元ベクトルから調律した32次元ベクトルを返す変換を表す行列
def temper_matrix():
    t0 = elem_matrix(3)
    t1 = elem_matrix(0)
    t2 = rshift_matrix(8) * elem_matrix(2)
    return t0 + (t1 + t2) + vect_to_matrix(int_to_f2(tmat)) * get_0bit * (t1 + t2)

# 128次元ベクトルから調律に使われるt1の値を返す変換を表す行列
def temper_t1_matrix():
    return elem_matrix(0)

# 128次元ベクトルから調律に使われるt2の値を返す変換を表す行列
def temper_t2_matrix():
    return rshift_matrix(8) * elem_matrix(2)

temper = temper_matrix()
temper_t1 = temper_t1_matrix()
temper_t2 = temper_t2_matrix()

# 128次元のstateから次のstateを得る変換を表す行列
def next_state_matrix():
    y = elem_matrix(3)
    x = mask_matrix(int_to_f2(0x7fffffff)) * elem_matrix(0) + elem_matrix(1) + elem_matrix(2)
    x += lshift_matrix(1) * x
    y += rshift_matrix(1) * y + x
    st0 = elem_matrix(1)
    st1 = elem_matrix(2) + vect_to_matrix(int_to_f2(mat1)) * get_0bit * y
    st2 = x + lshift_matrix(10) * y + vect_to_matrix(int_to_f2(mat2)) * get_0bit * y
    st3 = y
    return join_matrix_vertical([st0, st1, st2, st3])

next_state = next_state_matrix()

# 128次元ベクトルを4つの32次元ベクトルをつなげたものとみたとき第0要素の最上位bitを無視して127bitを得るという変換の行列
def from_128bit_to_127bit_matrix():
    m = [[F2(0) for j in range(128)] for i in range(127)]
    for i in range(31):
        m[i][i] = F2(1)
    for i in range(31, 127):
        m[i][i+1] = F2(1)
    return matrix(m)

# 逆に127bitから1bit水増しして128次元ベクトルを返すという変換を表す行列
def from_127bit_to_128bit_matrix():
    m = [[F2(0) for j in range(127)] for i in range(128)]
    for i in range(31):
        m[i][i] = F2(1)
    for i in range(31, 127):
        m[i+1][i] = F2(1)
    return matrix(m)

from_128bit_to_127bit = from_128bit_to_127bit_matrix()
from_127bit_to_128bit = from_127bit_to_128bit_matrix()

# 127次元のstateから次のstateを得る変換を表す行列
next_state_127 = from_128bit_to_127bit * next_state * from_127bit_to_128bit

# 127次元のstateから連続する127個の乱数のそれぞれ最下位bitを得る変換の行列
def from_state_to_random_stream_matrix():
    mats = [None for i in range(127)]
    step = Mat(GF(2), 127, 127).identity_matrix()
    A = get_0bit * temper * from_127bit_to_128bit
    for i in range(127):
        mats[i] = A * step
        step *= next_state_127
    return join_matrix_vertical(mats)

# stateから調律後の第1bitの連続した127個分を得る
def from_state_to_tempered_bit1_stream_matrix():
    mats = [None for i in range(127)]
    step = Mat(GF(2), 127, 127).identity_matrix()
    A = get_bit_matrix(1) * temper * from_127bit_to_128bit
    for i in range(127):
        mats[i] = A * step
        step *= next_state_127
    return join_matrix_vertical(mats)

# stateから調律に使われるt1の第0bitの連続した127個分を得る
def from_state_to_t1_stream_matrix():
    mats = [None for i in range(127)]
    step = Mat(GF(2), 127, 127).identity_matrix()
    A = get_0bit * temper_t1 * from_127bit_to_128bit
    for i in range(127):
        mats[i] = A * step
        step *= next_state_127
    return join_matrix_vertical(mats)

# stateから調律に使われるt2の第0bitの連続した127個分を得る
def from_state_to_t2_stream_matrix():
    mats = [None for i in range(127)]
    step = Mat(GF(2), 127, 127).identity_matrix()
    A = get_0bit * temper_t2 * from_127bit_to_128bit
    for i in range(127):
        mats[i] = A * step
        step *= next_state_127
    return join_matrix_vertical(mats)

from_state_to_random_stream = from_state_to_random_stream_matrix()

# 127個の連続した第0bitからstateを得る行列 (逆行列を計算)
from_random_stream_to_state = from_state_to_random_stream.inverse()

from_state_to_tempered_bit1_stream = from_state_to_tempered_bit1_stream_matrix()
from_state_to_t1_stream = from_state_to_t1_stream_matrix()
from_state_to_t2_stream = from_state_to_t2_stream_matrix()

# 連立一次方程式の解のパラメータから解を得る変換の行列
def from_param_to_state_matrix(solution0, basis):
    m = [[0 for j in range(len(basis)+1)] for i in range(127)]
    for i in range(127):
        for j in range(len(basis)+1):
            if j == 0:
                m[i][j] = solution0[i]
            else:
                m[i][j] = basis[j-1][i]
    return matrix(m)

# 係数がベクトルaの成分の一次式を文字列で返す
def poly(a):
    poly = []
    for i in range(len(a)):
        if a[i] != 0:
            if i == 0:
                poly.append("1")
            else:
                poly.append("x"+str(i-1))
    if poly != []:
        return "+".join(poly)
    else:
        return "0"

# a, b, cから決まる方程式を文字列で返す
def equation(a, b, c):
    return poly(a) + "+(" + poly(b) + ")*(" + poly(c) + ")=0"

# 乱数列の最下位bitが1,0,…,0,…(最初に1、続いて0が63個)であるようなstateを求めるという一次方程式を解き、解空間(solution1, basis)を得る
I = matrix([[xi(i == j) for j in range(127)] for i in range(64)])
vec = vector([xi(i == 0) for i in range(64)])
solution0 = (I * from_random_stream_to_state).solve_right(vec)
basis = kernel(((I * from_random_stream_to_state).transpose())).basis()

from_param_to_state = from_param_to_state_matrix(solution1, basis)

for i in range(len(basis)):
    a = list(nth_row(from_state_to_tempered_bit1_stream, i) * from_param_to_state)[0]
    b = list(nth_row(from_state_to_t1_stream, i) * from_param_to_state)[0]
    c = list(nth_row(from_state_to_t2_stream, i) * from_param_to_state)[0]
    print(equation(a, b, c))

得られた連立方程式は以下にアップロードしてあります。

出力関数が線形な場合のTinyMTにおいて連続する4つの乱数値からstateを得る

疑似乱数TinyMTの逆算について考えます。https://odanpoyo.github.io/2016/12/10/Tiny-Mersenne-twister%E3%81%AE%E9%80%86%E7%AE%97/の続きともいえるような内容です。
上の記事はstateから一つ前のstateを計算するという話でしたが、この記事では連続する4つの乱数値からstateを得ることを考えます。

TinyMTは何も設定しなければ出力関数が非線形ですが、LINEARITY_CHECKをONにすると線形になります。その場合、連続する4つの乱数値からstateを得るのはただの(F_2係数の)連立一次方程式を解くだけなので簡単にできます。

SageMath Cloudで動かせるプログラムを作りました。
動かす場合、メニューのKernel→Change Kernel→SageMath 7.4と選ぶ必要があるかもしれません。

なお、TinyMTの実装はTinyMT/tinymt32.h at master · MersenneTwister-Lab/TinyMT · GitHubのtinymt32_next_stateとtinymt32_temperを参照。

LINEARITY_CHECKがOFFの場合、つまり出力関数が非線形な場合は乱数値からstateを得るのは現実的な実行時間じゃ無理だろうなと思っていました。しかしそれをさき (@water_blow)さんがくつがえしてくれました。

そうなんです。LINEARITY_CHECKがOFFのとき以下のように算術加算が使われていて非線形なので逆算は不可能に見えます。しかし下位1bitを見ればxorしていることと変わらないので上のツイートのように逆算が可能になります。

	inline static uint32_t tinymt32_temper(tinymt32_t * random) {
		uint32_t t0, t1;
		t0 = random->status[3];
#if defined(LINEARITY_CHECK)
		t1 = random->status[0]
			^ (random->status[2] >> TINYMT32_SH8);
#else
		t1 = random->status[0]
			+ (random->status[2] >> TINYMT32_SH8);
#endif
		t0 ^= t1;
		t0 ^= -((int32_t)(t1 & 1)) & random->tmat;
		return t0;
	}

全く気が付きませんでした。すごいですね。さきさんがこの件について記事を書いてくれることを期待してみます。

追記(2016/12/11)
さきさんが書いてくれました。

追記(2016/12/11)
こちらでも連続した127個の乱数の下位1bitからstateを算出する行列を計算しました。
さきさんのツールと結果が一致します。

筆者: oupo (連絡先: oupo.nejiki@gmail.com)