はじめに

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

数独の特徴

詳細な説明は、前回の記事を参照してください。

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

その１：シンプル版

まずは、処理時間は少しかかるけれど、わかりやすさに重点を置いたソースコードを実装しました。

処理時間短縮の工夫

以下の代表的な３パターンの組み合わせ総数を調べれば、答えにたどり着くことがわかります。
残りは、法則１～法則５を使うことで計算不要になります。

青枠に囲まれた部分でパターン分けします。
左から、パターン１（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}$に短縮できました。

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時間にまで短縮できました。

コマンドプロンプト

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を使えばよいので、この記事の最後に補記しました。

その２：速度重視版

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

処理時間短縮の工夫

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

ソースコード

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

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);
}

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秒）にまで短縮できました。
参考にした論文では、2007年で処理時間が6時間との記述がありました。
当時と比べて、CPU、メモリ等のハードウェアやD言語を含むコンパイラ技術が進化したということでしょうね。

コマンドプロンプト

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パターンで実装してみました。

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);
}

コマンドプロンプト

D:\Dev> dmd -m64 bigint1.d

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

関連記事