Hatena::ブログ(Diary)

oupoの日記

2017-06-16

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

21:48

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

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

f:id:oupo:20170616213707p:image

f:id:oupo:20170616213708p:image

f:id:oupo:20170616213706p:image

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

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()

2017-01-08

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

14:28

背景

乱数調整のツールを作っていると次のようなコードを書きたくなります。

しかし、このコードは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のデストラクタが再帰的に呼ばれ、再帰の回数が大きくなりすぎてスタックオーバーフローするためです。

2016-12-13

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

00:01

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は次の連立方程式を解けば得られます。

f:id:oupo:20161212215141p:image

あとは、この方程式グレブナー基底を求めるなり、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))

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

2016-12-10

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

18:12

疑似乱数TinyMTの逆算について考えます。Page not found · GitHub Pagesの続きともいえるような内容です。

上の記事は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を算出する行列を計算しました。

さきさんのツールと結果が一致します。

2016-12-07

32bit LCGのseedを手計算で求めよう (半分ネタ)

11:54

 X_{n+1} = (a X_n + b) ¥bmod 2^{32}

     a = ¥mathrm{0x41c64e6d},¥;b = ¥mathrm{0x6073}

で与えられる線形合同法疑似乱数において、その上位16bit

Y_n = ¥lfloor X_n / 2^{16} ¥rfloor

が観測できるという状況を考えます。

そのとき $X_0$ の値を求めるにはどうすればよいでしょうか。

$Y_0$ を $\alpha$ とおくとき $X_0$ は

 X_0 = 2^{16} ¥alpha + x,¥;0 ¥le x < 2^{16}

とかけます。よって $2^{16}$ 通りの $x$ すべてをしらみつぶして、 $Y_1$ の値が一致するものを調べればよさそうです。$Y_1$ の値が一致するものは複数個ある可能性がある場合があるのでその場合は $Y_2$ の値を調べます。これで疑似乱数seedを特定できるでしょう。

しかし、$2^{16}$ 通り調べるのは手計算ではできないですね。

そこで手計算できる方法を考えます。

まず、

 X_{n+2^{16}} = (¥mathrm{0xdfa40001} X_n + ¥mathrm{0x47310000}) ¥bmod 2^{32}

であることに注意します。

特に係数がそれぞれ $2^{16}$ の倍数に1足した値と$2^{16}$ の倍数になっているので、それらを$2^{16} A + 1$, $2^{16} B$とおいてみます。

すると上の $X_0$ の候補の値について

 ¥begin{align*} X_{2^{16}} &= (2^{16} A + 1) (2^{16} ¥alpha + x) + 2^{16} B ¥¥ &= 2^{16}(Ax+¥alpha+B) + x¥end{align}

が成り立ちます。ここで $Y_{2^{16}}$ を $\beta$ とおけば

 Ax+¥alpha+B ¥equiv ¥beta ¥pmod{2^{16}}

なのでこれを解くことにより $x$ を得ます!

ただし $A$ は4の倍数になっているので、得られるのは $x \bmod 2^{14}$ です。したがってこの方法で得られるseedは4通りあるので、その後 $Y_1$ などを使って一通りに絞らないといけません。

上の方法をまとめると

1) 乱数を1つ生成し、その値(16bit)を得る ($\alpha$ とする)

2) 乱数を2^{16}-1 = 65535個進める

3) 乱数を1つ生成し、その値(16bit)を得る ($\beta$ とする)

4) $\alpha$ と $\beta$ からseedが4通りに絞れる

となります。

さて、よく見ると 0xdfa40001 は 2^{16} の倍数 + 1 であるだけでなく、 2^{18} の倍数 + 1 です。ここに気付くと上の方法は改良できます。

2^{16}個先を考えるのではなく2^{14}個先を考えるのです。

 X_{n+2^{14}} = (¥mathrm{0xf7e90001} X_n + ¥mathrm{0xf1cc4000}) ¥bmod 2^{32}

ここに出てくる係数をあらためて $2^{16} A + 1 = 0xf7e90001, 2^{14} B = 0xf1cc4000$ とおきます。

すると $X_0$ の候補の値について

 ¥begin{align*} X_{2^{14}} &= (2^{16} A + 1) (2^{16} ¥alpha + x) + 2^{14} B ¥¥ &= 2^{14}(4Ax+4¥alpha+B) + x¥end{align}

が成り立ちます。ここで $X_{2^{14}}$ を $2^{16} \beta + y$ とおけば

 ¥begin{align*} 2^{14}(4Ax+4¥alpha+B)+x = 2^{16} ¥beta + y.¥end{align}

$x = 2^{14} x' + x'', y = 2^{14} y' + y'', 0 \ge x', y' < 4, 0 \ge x'', y'' < 2^{14}$ とおけば

 ¥begin{align*} 2^{14}(4Ax+4¥alpha+B+x’)+x’’ = 2^{14} (4 ¥beta + y’) + y’’.¥end{align}

上位18bitをとることにより

 4Ax+4¥alpha+B + x’ ¥equiv 4 ¥beta + y’ ¥pmod{2^{18}}

$y' - x' = \delta$とけば

 4Ax+4¥alpha+B ¥equiv 4 ¥beta + ¥delta ¥pmod{2^{18}} (1)

     -4 < ¥delta < 4

この式に適する$\delta$を選び、$x$について解くと解の候補が得られます。$x$と$\delta$の整合性をチェックし、解になっているか確かめられます。

上の方法をまとめると

1) 乱数を1つ生成し、その値(16bit)を得る ($\alpha$ とする)

2) 乱数を2^{14}-1 = 16383個進める

3) 乱数を1つ生成し、その値(16bit)を得る ($\beta$ とする)

4) $\alpha$ と $\beta$ からseedが最大2通り(最小0通り)に決定される

となります。

試しにやってみましょう。$\alpha = 11111, \beta = 22222$としてみます。

式(1)は

 253860x+44444+247601 ¥equiv 88888 + ¥delta ¥pmod{2^{18}}

です。整理すると

 253860x ¥equiv -203157 + ¥delta ¥pmod{2^{18}}.

$-203157 + \delta$は4の倍数にならざるを得ないことから$\delta = 1, -3$です。

$\delta = 1$のとき:

 253860x ¥equiv -203156 ¥pmod{2^{18}}.

両辺を4で割り、

 63465x ¥equiv -50789 ¥pmod{2^{16}}.

63465の2^{16}を法とした逆元はユークリッドの互除法により20569と求まるので

 x ¥equiv 20569 ¥cdot (-50789) ¥pmod{2^{16}}.

計算すると

 x ¥equiv 30435 ¥pmod{2^{16}}

と求まります。

このとき$x' = 1$なので$\delta = 1$と整合します ($\delta = y' - x'$をみたす$0 \ge y' < 4$があります)。

$\delta = -3$のとき:

同様の計算で

 x ¥equiv 51004 ¥pmod{2^{16}}

と求まります。

このとき$x' = 3$なので$\delta = 1$と整合しません ($\delta = y' - x'$をみたす$0 \ge y' < 4$がありません)。

よって解は30435の一つです。

ゆえにseed 2^{16} ¥cdot 11111 + 30435 = 728200931 = ¥mathrm{0x2b6776e3}です。

2^{14}回乱数を進めるなんて2^{16}通り機械でしらみつぶしするよりずっと大変じゃないかということで、半分ネタ記事でした。

追記 (2016/12/8)

2^{14}の方で必ず解が一通りになるという勘違いをしていてその間違いを残したまま公開していました。現在は修正済みです。

外野席外野席 2016/12/07 20:42 >>2^{14}回乱数を進めるなんて2^{16}通り機械でしらみつぶしするよりずっと大変じゃないか

~14 にしても ~16にしても 通常の手計算で対処すべき範疇をかなり逸脱しているかと
あえて申し上げます
ろいしんさんのスタイルは私も好きですが、あなたがそれに倣うことはお勧めいたしがたい