9x9 数独の答えになる組み合わせは何通りか

# はじめに
この記事では、9x9 数独の答えになる組み合わせが何通りあるかを調べます。
D言語で並列処理を実装することで、処理時間を大幅に改善します。

# 数独の特徴
詳細な説明は、[前回の記事](https://qiita.com/devmynote/items/7e678379a014646267e8)を参照してください。

* 法則１：数字は交換が可能
* 法則２：同じブロック内の行は交換が可能
* 法則３：同じブロック内の列は交換が可能
* 法則４：ブロック行の交換が可能
* 法則５：ブロック列の交換が可能

# その１：シンプル版
まずは、処理時間は少しかかるけれど、わかりやすさに重点を置いたソースコードを実装しました。

## 処理時間短縮の工夫
以下の代表的な３パターンの組み合わせ総数を調べれば、答えにたどり着くことがわかります。
残りは、法則１～法則５を使うことで計算不要になります。
![sudoku99a.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/675901/aa60124d-1f3c-3483-fc16-a73cc9fd9f09.png)

青枠に囲まれた部分でパターン分けします。
左から、パターン１（456789型：黄色と緑が１列に並ぶ）、パターン２（457689型：緑の下に緑、黄色の下に黄色、緑の下に黄色がそれぞれ１つづつ来る）、パターン３（457986型：緑の下に黄色、黄色の下に緑が来る）とすると、以下の通りになります。
$[すべての組合せ]=(a+b\times6+c\times3)\times2_{※1}\times(3!)^2_{※2}\times9!_{※3}$
$a = [パターン１の組合せ総数]$
$b = [パターン２の組合せ総数]$
$c = [パターン３の組合せ総数]$
$※1：青枠の行について、上下並べ替えが可能$
$※2：青枠の緑と黄色のマーカーそれぞれについて、並べ替えが可能$
$※3：赤枠の数字について、交換が可能$
となります。

また、左下側のピンク色の枠について、$6! = 720$通り調べるところを代表的な$10$パターンについて調べ、あとは法則２、法則３を使って、$2\times3!\times3! = 72$倍することで、処理速度が短縮できます。
10パターンは、以下の通りに並べ替えられるものに絞れます。

* $A < D$
* $A < B < C$
* $D < E < F$

## ソースコード
さらに処理速度を短縮するために、D言語で並列処理を実装しました。
並列処理したい`foreach`ループで`parallel`を使うだけなので、驚くほど簡単です。
10個の並列処理を実行することで、処理時間をおよそ$\frac{1}{8}$に短縮できました。

``` d:sudoku99.d
import std.algorithm : remove;
import std.algorithm.iteration : filter, permutations;
import std.algorithm.sorting: nextPermutation;
import std.array : array;
import std.datetime;
import std.parallelism: parallel;
import std.stdio;

const NUM = 9;
const ROW = NUM;
const COL = NUM;
const BLOCK_ROW = 3;
const BLOCK_COL = NUM / BLOCK_ROW;
ubyte[COL][ROW] rect;
bool[NUM][ROW] yoko;
bool[NUM][COL] tate;
bool[NUM][BLOCK_COL][BLOCK_ROW] block;
long count;

bool push(int x, int y, int n)
{
	if ( yoko[y][n] == false
			&& tate[x][n] == false
			&& block[y / BLOCK_COL][x / BLOCK_ROW][n] == false ){
		rect[y][x] = cast(ubyte)(n + 1);
		yoko[y][n] = true;
		tate[x][n] = true;
		block[y / BLOCK_COL][x / BLOCK_ROW][n] = true;
		return ( true );
	}
	return ( false );
}

void pop(int x, int y, int n)
{
	rect[y][x] = 0;
	yoko[y][n] = false;
	tate[x][n] = false;
	block[y / BLOCK_COL][x / BLOCK_ROW][n] = false;
}

void next(int x, int y)
{
	if ( x == COL - 1 ){
		if ( y == ROW - 2 ){
			count++;
			return;
		} else {
			x = 0;
			y++;
		}
	} else {
		x++;
	}
	solve(x, y);
}

void solve(int x, int y)
{
	if ( rect[y][x] > 0 ){
		next(x, y);
		return;
	}
	for ( int n = 0; n < NUM; n++ ){
		if ( push(x, y, n) ){
			next(x, y);
			pop(x, y, n);
		}
	}
}

long countRect(ref ubyte[COL][ROW] r1)
{
	int[] left = [1, 2, 3, 4, 5, 6, 7, 8, 9];
	for ( int y = 0; y < 3; y++ ){
		left = remove!(a=> a == r1[y][0])(left);
	}
	auto lefts = left.permutations.filter!(
		a => a[0] < a[3] && a[0] < a[1] && a[1] < a[2] && a[3] < a[4] && a[4] < a[5]);
	ubyte[COL][ROW][] rectarr;
	foreach ( arr; lefts ){
		foreach ( y, n; arr.array ){
			r1[y + 3][0] = cast(ubyte)n;
		}
		rectarr ~= r1;
	}
	__gshared ulong total;
	total = 0;
	foreach ( r0; parallel(rectarr) ){
		count = 0;
		yoko = new bool[NUM][ROW];
		tate = new bool[NUM][COL];
		block = new bool[NUM][BLOCK_COL][BLOCK_ROW];
		rect = r0;
		for ( int y = 0; y < ROW; y++ ){
			for ( int x = 0; x < COL; x++ ){
				if ( rect[y][x] > 0 ){
					if ( push(x, y, rect[y][x] - 1) == false ){
						writeln("Illegal initial value");
					}
				}
			}
		}
		solve(0, 0);
		writefln("%(%(%d %)\n%)\n---%,3d", rect, count);
		synchronized {
			total += count;
		}
	}
	return ( total );
}

long solvePattern(ubyte val1, ubyte val2, ubyte val3, ubyte[] arr)
{
	long starttime = Clock.currStdTime();
	ulong subTotal;
	do {
		rect = [[1, 2, 3, 4, 5, 6, 7, 8, 9],
				[4, 5, 6, 0, 0, 0, 0, 0, 0],
				[7, 8, 9, 0, 0, 0, 0, 0, 0],
				[0, 0, 0, 0, 0, 0, 0, 0, 0],
				[0, 0, 0, 0, 0, 0, 0, 0, 0],
				[0, 0, 0, 0, 0, 0, 0, 0, 0],
				[0, 0, 0, 0, 0, 0, 0, 0, 0],
				[0, 0, 0, 0, 0, 0, 0, 0, 0],
				[0, 0, 0, 0, 0, 0, 0, 0, 0]];
		rect[1][2] = val1;
		rect[2][0] = val2;
		rect[2][2] = val3;
		rect[1][3..6] = arr[];
		subTotal += countRect(rect);
	} while( arr.nextPermutation );
	long time = (Clock.currStdTime() - starttime) / 10_000_000;
	writefln("sub total = %,3d", subTotal);
	writefln("time      = %,3d s", time);
	return ( subTotal );
}

void main()
{
	long total  = solvePattern(6, 7, 9, [7, 8, 9]);
		 total += solvePattern(7, 6, 9, [0, 8, 9]) * 6;
		 total += solvePattern(7, 9, 6, [0, 8, 9]) * 3;
	writefln("total     = %,3d * 2 * 36 * 72 * 362,880", total);
}
```

## コンパイル、実行結果（答え）
元のソースコード（4x4の答え算出版）を9x9用に定数を変更するだけだと、処理時間が1200年ほどかかる見込みでした。
ただ、これまで見てきた通りの見直しを行った結果、処理時間が33.3時間にまで短縮できました。

``` text:コマンドプロンプト
D:\Dev> ldc2 -m64 -release -O sudoku99.d

D:\Dev> sudoku99.exe

・・・省略・・・

1 2 3 4 5 6 7 8 9
4 5 6 9 8 7 0 0 0
7 8 9 0 0 0 0 0 0
2 0 0 0 0 0 0 0 0
3 0 0 0 0 0 0 0 0
6 0 0 0 0 0 0 0 0
5 0 0 0 0 0 0 0 0
8 0 0 0 0 0 0 0 0
9 0 0 0 0 0 0 0 0
---2,281,587,664
sub total = 127,491,348,608
time      = 17,159 s

・・・省略・・・

1 2 3 4 5 6 7 8 9
4 5 7 9 8 0 0 0 0
6 8 9 0 0 0 0 0 0
2 0 0 0 0 0 0 0 0
3 0 0 0 0 0 0 0 0
7 0 0 0 0 0 0 0 0
5 0 0 0 0 0 0 0 0
8 0 0 0 0 0 0 0 0
9 0 0 0 0 0 0 0 0
---6,825,771,168
sub total = 379,586,179,808
time      = 51,519 s

・・・省略・・・

1 2 3 4 5 6 7 8 9
4 5 7 9 8 0 0 0 0
9 8 6 0 0 0 0 0 0
2 0 0 0 0 0 0 0 0
5 0 0 0 0 0 0 0 0
6 0 0 0 0 0 0 0 0
3 0 0 0 0 0 0 0 0
7 0 0 0 0 0 0 0 0
8 0 0 0 0 0 0 0 0
---6,736,710,608
sub total = 380,379,290,944
time      = 51,149 s
total     = 3,546,146,300,288 * 2 * 36 * 72 * 362,880
```

答えの数値がlong型(64bit)の上限を超えるため、やや手抜きの出力結果になっていますが、正しい答えになっています！！
D言語で大きい数値を扱う場合、`BigInt`を使えばよいので、この記事の最後に補記しました。

# その２：速度重視版
シンプル版と比べてソースコードが長くなりますが、さらなる処理時間短縮に重点を置いたソースコードになります。

## 処理時間短縮の工夫
前回記事のコメント欄から[情報提供いただいた論文](http://kymst.net/sudoku_test/doc/SUDOKU_Index.pdf)を参考にすると、さらに処理時間が短縮できます。
シンプル版で調べた３パターンについて、上から3行目までを解析して、代表的なパターンを絞るという考え方です。
シンプル版だと、上の3行部分について$9072 = (1 + 3 + 3) * (3!)^4 $通り調べていることになりますが、速度重視版では$71$通りにまで絞れます。
これにより、処理時間をおよそ127倍（$=\frac{9072}{71}$）短縮することができました。

## ソースコード
`sudoku99x3.d`は上の3行を代表的な71パターンに絞り込む処理、`sudoku99x6.d`は下の6行の組合せを数える処理を実装しています。

``` d:sudoku99x3.d
import std.algorithm : sort;
import std.algorithm.mutation : swap;
import std.datetime;
import std.stdio;
import sudoku99x6;

ubyte[COL][ROW] rect;
bool[NUM][ROW] yoko;
bool[NUM][COL] tate;
bool[NUM][BLOCK_COL][BLOCK_ROW] block;
long count;
long[ubyte[COL][ROW]] cache1;
ubyte[COL][ROW][][long] cache2;

bool push(int x, int y, int n)
{
	if ( yoko[y][n] == false
			&& tate[x][n] == false
			&& block[y / BLOCK_COL][x / BLOCK_ROW][n] == false ){
		rect[y][x] = cast(ubyte)(n + 1);
		yoko[y][n] = true;
		tate[x][n] = true;
		block[y / BLOCK_COL][x / BLOCK_ROW][n] = true;
		return ( true );
	}
	return ( false );
}

void pop(int x, int y, int n)
{
	rect[y][x] = 0;
	yoko[y][n] = false;
	tate[x][n] = false;
	block[y / BLOCK_COL][x / BLOCK_ROW][n] = false;
}

void next3row(int x, int y)
{
	if ( x == COL - 1 ){
		if ( y == 2 ){
			cache1[rect] = count;
			cache2[count] ~= rect;
			count++;
			return;
		} else {
			x = 0;
			y++;
		}
	} else {
		x++;
	}
	solve3row(x, y);
}

void solve3row(int x, int y)
{
	if ( rect[y][x] > 0 ){
		next3row(x, y);
		return;
	}
	for ( int n = 0; n < NUM; n++ ){
		if ( push(x, y, n) ){
			next3row(x, y);
			pop(x, y, n);
		}
	}
}

void swapCol(ref ubyte[COL][ROW] r1, int x1, int x2)
{
	swap(r1[0][x1], r1[0][x2]);
	swap(r1[1][x1], r1[1][x2]);
	swap(r1[2][x1], r1[2][x2]);
}

void standardization(ref ubyte[COL][ROW] r1)
{
	int n = 0;
	for ( int x = 6; x < 9; x++ ){
		if ( r1[0][x] == r1[1][0] || r1[0][x] == r1[1][1] || r1[0][x] == r1[1][2] ){
			n++;
		}
	}
	if ( n >= 2 ){
		for ( int x = 3; x < 6; x++ ){
			swapCol(r1, x, x + 3);
		}
		n = 3 - n;
	}
	for ( int x = 0; x < 3; x++ ){
		ubyte v1 = r1[1][x];
		if ( r1[0][6] == v1 || r1[0][7] == v1 || r1[0][8] == v1 ){
			swapCol(r1, x, 2);
			for ( int x1 = 7; x1 < 9; x1++ ){
				if ( r1[0][x1] == v1 ){
					swapCol(r1, x1, 6);
				}
			}
		}
	}
	for ( int x = 3; x < 5; x++ ){
		ubyte v1 = r1[0][x];
		if ( r1[1][0] != v1 && r1[1][1] != v1 && r1[1][2] != v1 ){
			swapCol(r1, x, 5);
			break;
		}
	}
	if ( r1[2][1] == r1[0][5] ){
		swapCol(r1, 0, 1);
	}
	if ( r1[1][0] == r1[0][4] ){
		swapCol(r1, 3, 4);
	}
	if ( n == 0 && r1[2][0] == r1[0][8] ){
		swapCol(r1, 6, 8);
	}
	if ( r1[2][1] == r1[0][8] ){
		swapCol(r1, 7, 8);
	}
	ubyte[COL] v;
	for ( int x = 0; x < COL; x++ ){
		v[r1[0][x]] = cast(ubyte)(x + 1);
	}
	for ( int y = 0; y < 3; y++ ){
		for ( int x = 0; x < COL; x++ ){
			r1[y][x] = v[r1[y][x]];
		}
	}
}

void updateCache(ref ubyte[COL][ROW] r0, ref ubyte[COL][ROW] r1)
{
	void updateCacheInner(long before, long after){
		auto rx1 = cache2[before];
		foreach ( r0; rx1 ){
			cache1[r0] = after;
		}
		auto rx2 = cache2[after];
		rx2 ~= rx1;
		cache2[before] = rx2;
		cache2[after]  = rx2;
	}
	long* p0 = (r0 in cache1);
	long* p1 = (r1 in cache1);
	if ( p1 !is null ){
		if ( *p0 > *p1 ){
			updateCacheInner(*p0, *p1);
		} else if ( *p0 < *p1 ) {
			updateCacheInner(*p1, *p0);
		}
	} else {
//		writefln("---out of range\n%(%(%d %)\n%)", r1);
	}
}

void trans2x2()
{
	void trans2x2inner(ref ubyte[COL][ROW] r0, int x1, int x2, int y1, int y2){
		if ( r0[y1][x1] == r0[y2][x2] && r0[y1][x2] == r0[y2][x1] ){
			ubyte[COL][ROW] r1 = r0;
			swap(r1[y1][x1], r1[y2][x1]);
			swap(r1[y1][x2], r1[y2][x2]);
			standardization(r1);
			updateCache(r0, r1);
		}
	}
	auto keys = cache1.keys;
	foreach ( r0; keys.sort!("a < b") ){
		for ( int x1 = 0; x1 < 6; x1++ ){
			for ( int x2 = (x1 / 3 + 1) * 3; x2 < 9; x2++ ){
				trans2x2inner(r0, x1, x2, 0, 1);
				trans2x2inner(r0, x1, x2, 0, 2);
				trans2x2inner(r0, x1, x2, 1, 2);
			}
		}
	}
}

void trans3x2()
{
	void trans3x2inner(ref ubyte[COL][ROW] r0, int x1, int x2, int x3, int y1, int y2){
		if (  (r0[y1][x1] == r0[y2][x2] && r0[y1][x2] == r0[y2][x3] && r0[y1][x3] == r0[y2][x1])
			||(r0[y1][x1] == r0[y2][x3] && r0[y1][x2] == r0[y2][x1] && r0[y1][x3] == r0[y2][x2])){
			ubyte[COL][ROW] r1 = r0;
			swap(r1[y1][x1], r1[y2][x1]);
			swap(r1[y1][x2], r1[y2][x2]);
			swap(r1[y1][x3], r1[y2][x3]);
			standardization(r1);
			updateCache(r0, r1);
		}
	}
	auto keys = cache1.keys;
	foreach ( r0; keys.sort!("a < b") ){
		for ( int x1 = 0; x1 < 3; x1++ ){
			for ( int x2 = 3; x2 < 6; x2++ ){
				for ( int x3 = 6; x3 < 9; x3++ ){
					trans3x2inner(r0, x1, x2, x3, 0, 1);
					trans3x2inner(r0, x1, x2, x3, 0, 2);
					trans3x2inner(r0, x1, x2, x3, 1, 2);
				}
			}
		}
	}
}

void trans4x2()
{
	void trans4x2inner(ref ubyte[COL][ROW] r0, int x1, int x2, int x3, int x4, int y1, int y2){
		if (  (r0[y1][x1] == r0[y2][x4] && r0[y1][x2] == r0[y2][x3] && r0[y1][x3] == r0[y2][x1] && r0[y1][x4] == r0[y2][x2])
			||(r0[y1][x1] == r0[y2][x3] && r0[y1][x2] == r0[y2][x4] && r0[y1][x3] == r0[y2][x2] && r0[y1][x4] == r0[y2][x1])){
			ubyte[COL][ROW] r1 = r0;
			swap(r1[y1][x1], r1[y2][x1]);
			swap(r1[y1][x2], r1[y2][x2]);
			swap(r1[y1][x3], r1[y2][x3]);
			swap(r1[y1][x4], r1[y2][x4]);
			standardization(r1);
			updateCache(r0, r1);
		}
	}
	void trans4x2sub1(ref ubyte[COL][ROW] r0, int xx1, int xx2, int xx3){
		for ( int x1 = xx1; x1 < xx1 + 2; x1++ ){
			for ( int x2 = x1 + 1; x2 < xx1 + 3; x2++ ){
				for ( int x3 = xx2; x3 < xx2 + 3; x3++ ){
					for ( int x4 = xx3; x4 < xx3 + 3; x4++ ){
						trans4x2inner(r0, x1, x2, x3, x4, 0, 1);
						trans4x2inner(r0, x1, x2, x3, x4, 0, 2);
						trans4x2inner(r0, x1, x2, x3, x4, 1, 2);
					}
				}
			}
		}
	}
	void trans4x2sub2(ref ubyte[COL][ROW] r0, int xx1, int xx2){
		for ( int x1 = xx1; x1 < xx1 + 2; x1++ ){
			for ( int x2 = x1 + 1; x2 < xx1 + 3; x2++ ){
				for ( int x3 = xx2; x3 < xx2 + 2; x3++ ){
					for ( int x4 = x3 + 1; x4 < xx2 + 3; x4++ ){
						trans4x2inner(r0, x1, x2, x3, x4, 0, 1);
						trans4x2inner(r0, x1, x2, x3, x4, 0, 2);
						trans4x2inner(r0, x1, x2, x3, x4, 1, 2);
					}
				}
			}
		}
	}
	auto keys = cache1.keys;
	foreach ( r0; keys.sort!("a < b") ){
		trans4x2sub1(r0, 0, 3, 6);
		trans4x2sub1(r0, 3, 0, 6);
		trans4x2sub1(r0, 6, 0, 3);
		trans4x2sub2(r0, 0, 3);
		trans4x2sub2(r0, 0, 6);
		trans4x2sub2(r0, 0, 9);
	}
}

void trans2x3()
{
	void trans2x3inner(ref ubyte[COL][ROW] r0, int x1, int x2){
		if (  (r0[0][x1] == r0[1][x2] && r0[1][x1] == r0[2][x2] && r0[2][x1] == r0[0][x2])
			||(r0[0][x1] == r0[2][x2] && r0[1][x1] == r0[0][x2] && r0[2][x1] == r0[1][x2])){
			ubyte[COL][ROW] r1 = r0;
			swapCol(r1, x1, x2);
			standardization(r1);
			updateCache(r0, r1);
		}
	}
	auto keys = cache1.keys;
	foreach ( r0; keys.sort!("a < b") ){
		for ( int x1 = 0; x1 < 6; x1++ ){
			for ( int x2 = (x1 / 3 + 1) * 3; x2 < 9; x2++ ){
				trans2x3inner(r0, x1, x2);
			}
		}
	}
}

void transBlockColSwap()
{
	void transBlockColSwapInner(ref ubyte[COL][ROW] r0, int xx1, int xx2){
		ubyte[COL][ROW] r1 = r0;
		for ( int x1 = 0; x1 < 3; x1++ ){
			swapCol(r1, x1 + xx1, x1 + xx2);
		}
		for ( int x1 = 0; x1 < COL; x1++ ){
			swap(r1[1][x1], r1[2][x1]);
		}
		standardization(r1);
		updateCache(r0, r1);
	}
	auto keys = cache1.keys;
	foreach ( r0; keys.sort!("a < b") ){
		transBlockColSwapInner(r0, 0, 3);
		transBlockColSwapInner(r0, 0, 6);
		transBlockColSwapInner(r0, 3, 6);
	}
}

void transRowSwap()
{
	auto keys = cache1.keys;
	foreach ( r0; keys.sort!("a < b") ){
		ubyte[COL][ROW] r1 = r0;
		for ( int y1 = 0; y1 < 2; y1++ ){
			for ( int y2 = y1 + 1; y2 < 3; y2++ ){
				for ( int x1 = 0; x1 < 9; x1++ ){
					swap(r1[y1][x1], r1[y2][x1]);
				}
				standardization(r1);
				updateCache(r0, r1);
			}
		}
	}
}

void solvePattern(ubyte val1, ubyte val2, ubyte val3)
{
	rect = [[1, 2, 3, 4, 5, 6, 7, 8, 9],
			[4, 5, 6, 0, 0, 0, 0, 0, 0],
			[7, 8, 9, 0, 0, 0, 0, 0, 0],
			[0, 0, 0, 0, 0, 0, 0, 0, 0],
			[0, 0, 0, 0, 0, 0, 0, 0, 0],
			[0, 0, 0, 0, 0, 0, 0, 0, 0],
			[0, 0, 0, 0, 0, 0, 0, 0, 0],
			[0, 0, 0, 0, 0, 0, 0, 0, 0],
			[0, 0, 0, 0, 0, 0, 0, 0, 0]];
	rect[1][2] = val1;
	rect[2][0] = val2;
	rect[2][2] = val3;
	yoko = new bool[NUM][ROW];
	tate = new bool[NUM][COL];
	block = new bool[NUM][BLOCK_COL][BLOCK_ROW];
	for ( int y = 0; y < ROW; y++ ){
		for ( int x = 0; x < COL; x++ ){
			if ( rect[y][x] > 0 ){
				if ( push(x, y, rect[y][x] - 1) == false ){
					writeln("Illegal initial value");
					return;
				}
			}
		}
	}
	solve3row(0, 0);
}

void main()
{
	long starttime = Clock.currStdTime();
	solvePattern(6, 7, 9);
	solvePattern(7, 6, 9);
	solvePattern(7, 9, 6);
	trans2x2();
	trans3x2();
	trans4x2();
	trans2x3();
	transBlockColSwap();
	transRowSwap();
	long[] arr = new long[cache1.length];
	auto keys = cache1.keys;
	keys.sort!("a < b");
	foreach ( r0; keys ){
		long* p0 = (r0 in cache1);
		if ( r0[2][0] == 7 ){
			arr[*p0] += 1;
		} else if ( r0[2][0] == 6 ){
			arr[*p0] += 6;
		} else { // r0[2][0] == 9
			arr[*p0] += 3;
		}
	}
	int ind;
	long total;
	foreach ( id, num; arr ){
		if ( num > 0 ){
			ind++;
			writefln("%02s(%04s) x %s", ind, id, num);
			auto rx1 = cache2[id];
			foreach ( r1; rx1 ){
				long subTotal = countRect(r1);
				writefln("sub total = %,3d", subTotal);
				total += num * subTotal;
				long time = (Clock.currStdTime() - starttime) / 10_000_000;
				writefln("time   = %,3d s", time);
				break;
			}
		}
	}
	writefln("total = %,3d * 2 * 36 * 72 * 362,880", total);
}
```
``` d:sudoku99x6.d
module sudoku99x6;

import std.algorithm : remove;
import std.algorithm.iteration : filter, permutations;
import std.array : array;
import std.parallelism: parallel;
import std.stdio;

const NUM = 9;
const ROW = NUM;
const COL = NUM;
const BLOCK_ROW = 3;
const BLOCK_COL = NUM / BLOCK_ROW;
ubyte[COL][ROW] rect;
bool[NUM][ROW] yoko;
bool[NUM][COL] tate;
bool[NUM][BLOCK_COL][BLOCK_ROW] block;
long count;

bool push(int x, int y, int n)
{
	if ( yoko[y][n] == false
			&& tate[x][n] == false
			&& block[y / BLOCK_COL][x / BLOCK_ROW][n] == false ){
		rect[y][x] = cast(ubyte)(n + 1);
		yoko[y][n] = true;
		tate[x][n] = true;
		block[y / BLOCK_COL][x / BLOCK_ROW][n] = true;
		return ( true );
	}
	return ( false );
}

void pop(int x, int y, int n)
{
	rect[y][x] = 0;
	yoko[y][n] = false;
	tate[x][n] = false;
	block[y / BLOCK_COL][x / BLOCK_ROW][n] = false;
}

void next(int x, int y)
{
	if ( x == COL - 1 ){
		if ( y == ROW - 2 ){
			count++;
			return;
		} else {
			x = 0;
			y++;
		}
	} else {
		x++;
	}
	solve(x, y);
}

void solve(int x, int y)
{
	if ( rect[y][x] > 0 ){
		next(x, y);
		return;
	}
	for ( int n = 0; n < NUM; n++ ){
		if ( push(x, y, n) ){
			next(x, y);
			pop(x, y, n);
		}
	}
}

long countRect(ref ubyte[COL][ROW] r1)
{
	int[] left = [1, 2, 3, 4, 5, 6, 7, 8, 9];
	for ( int y = 0; y < 3; y++ ){
		left = remove!(a=> a == r1[y][0])(left);
	}
	auto lefts = left.permutations.filter!(
		a => a[0] < a[3] && a[0] < a[1] && a[1] < a[2] && a[3] < a[4] && a[4] < a[5]);
	ubyte[COL][ROW][] rectarr;
	foreach ( arr; lefts ){
		foreach ( y, n; arr.array ){
			r1[y + 3][0] = cast(ubyte)n;
		}
		rectarr ~= r1;
	}
	__gshared ulong total;
	total = 0;
	foreach( r0; parallel(rectarr) ){
		count = 0;
		yoko = new bool[NUM][ROW];
		tate = new bool[NUM][COL];
		block = new bool[NUM][BLOCK_COL][BLOCK_ROW];
		rect = r0;
		for ( int y = 0; y < ROW; y++ ){
			for ( int x = 0; x < COL; x++ ){
				if ( rect[y][x] > 0 ){
					if ( push(x, y, rect[y][x] - 1) == false ){
						writeln("Illegal initial value");
					}
				}
			}
		}
		solve(0, 0);
		writefln("%(%(%d %)\n%)\n---%,3d", rect, count);
		synchronized {
			total += count;
		}
	}
	return ( total );
}
```

## コンパイル、実行結果（答え）
処理時間がおよそ16分（958秒）にまで短縮できました。
参考にした[論文](http://kymst.net/sudoku_test/doc/SUDOKU_Index.pdf)では、2007年で処理時間が6時間との記述がありました。
当時と比べて、CPU、メモリ等のハードウェアやD言語を含むコンパイラ技術が進化したということでしょうね。

``` text:コマンドプロンプト
D:\Dev> ldc2 -m64 -release -O sudoku99x3.d sudoku99x6.d

D:\Dev> sudoku99x3.exe
01(0000) x 4
1 2 3 4 5 6 7 8 9
4 5 6 7 8 9 1 2 3
7 8 9 1 2 3 4 5 6
2 0 0 0 0 0 0 0 0
5 0 0 0 0 0 0 0 0
8 0 0 0 0 0 0 0 0
3 0 0 0 0 0 0 0 0
6 0 0 0 0 0 0 0 0
9 0 0 0 0 0 0 0 0
---6,696,000

・・・省略・・・

71(3793) x 54

・・・省略・・・

1 2 3 4 5 6 7 8 9
4 5 7 8 9 3 2 6 1
6 8 9 2 1 7 3 5 4
2 0 0 0 0 0 0 0 0
3 0 0 0 0 0 0 0 0
7 0 0 0 0 0 0 0 0
5 0 0 0 0 0 0 0 0
8 0 0 0 0 0 0 0 0
9 0 0 0 0 0 0 0 0
---11,186,176
sub total = 98,153,104
time   = 958 s
total = 3,546,146,300,288 * 2 * 36 * 72 * 362,880
```

# 補足：64ビットの上限を超える数値の出力
9x9 数独の答えになる組み合わせは、6,670,903,752,021,072,936,960通り（`= 3,546,146,300,288 * 2 * 36 * 72 * 362,880`）であることがわかりました。
この計算結果を「そのまま出力」と「カンマ区切り」の2パターンで実装してみました。

``` d:bigint1.d
import std.bigint;
import std.format;
import std.range;
import std.stdio;

void main(){
	long n = 3_546_146_300_288;
	BigInt bi = BigInt(n) * 2 * 36 * 72 * 362_880;
	writeln(bi);
	
	string s = bi.toDecimalString;
	writeln(format("%(%s,%)", bi.toDecimalString.retro.chunks(3)).retro);
}
```
``` text:コマンドプロンプト
D:\Dev> dmd -m64 bigint1.d

D:\Dev> bigint1
6670903752021072936960
6,670,903,752,021,072,936,960
```

# 関連記事

https://qiita.com/devmynote/items/017ca6c32dd1e45be4b9

https://qiita.com/devmynote/items/7e678379a014646267e8