关于数独最终合法结果的计数

The number of valid Sudoku solution grids for the standard 9×9 grid was calculated by Bertram Felgenhauer in 2005 to be 6,670,903,752,021,072,936,960 ,This number is equal to 9! × 72^2 × 2^7 × 27,704,267,971, the last factor of which is prime.
上面问题是说合法的数独最终结果总共有6,670,903,752,021,072,936,960 个。

有什么好的算法可以得到6,670,903,752,021,072,936,960 ？或者27,704,267,971，722 各自怎么来 

其实通过google搜索27,704,267,971就可以搜索到一些相关的文章了。

首先，其中系数9！×72^2是比较容易理解的：

9!就是说明总是将9个数字进行置换，将左上角的3*3个方格中的数字变成正好顺序排列了。
而72^2是说分别对于行和列做相同的标准化操作。比如对于列
中间的3列可以任意列交换使得它们第一行中的数据单调增（有3!=6个置换）
右边的3列可以任意列交换使得它们第一行中的数据单调增（有3!=6个置换）
然后中间3列和右边3列也可以交换，使得第4列中第一个数小于第7列中第一个数。（2中置换）
总共有6*6*2=72中置换。

经过上面的标准化后，每个标准化后的局面可以代表9!×72^2个局面了。

然后我们的程序需要搜索出2^7 × 27,704,267,971个标准局面就可以了。

我的程序，分成两部分，第一部分，对所有的局面根据上面3行的内容进行分类（也就是不管余下6行的内容），可以分出101个类。

/*Classification of valid sudoku result*/
#include <stdio.h>
#include <stdlib.h>

#ifdef WIN32
typedef  __int64 longlong;
#define LLD "%I64d"
#else
typedef long long longlong;
#define LLD "%lld"
#endif

#define L0 ((longlong)0)
#define L1 ((longlong)1)
#define MAX_COUNT 36288
#define MAX_CLASS 416
short stack[MAX_COUNT][14];
int stack_count;
int board[9][9];
longlong count;
int class_count;
int class_number;
typedef int template;
#define TEMPLATE_COUNT 5184
#define MASK1 ((1<<12)-1)

template temps[100][9][TEMPLATE_COUNT];
int temps_count[100][9];
int class_id[MAX_COUNT];
int class_elem[MAX_CLASS];
int class_weight[MAX_CLASS];

signed char is_index[9][9];
signed char is_value[9][9];

int index_value[900];
int index_pos[900];
int index_count;
int sort[10][6]={
{0,1,2,3,4,5},
{0,1,3,2,4,5},
{0,1,4,2,3,5},
{0,1,5,2,3,4},
{0,2,3,1,4,5},
{0,2,4,1,3,5},
{0,2,5,1,3,4},
{0,3,4,1,2,5},
{0,3,5,1,2,4},
{0,4,5,1,2,3}
};
int row[6]={4,5,6,7,8,9};
int column[6]={2,3,5,6,8,9};
int sorder[9]={8,3,1,4,2,7,5,6,0};
int class[100][9];

//Decode the compressed representation of the first 3 columns in statck[stack_pointer]
void set_board(int stack_pointer){
    int i,j;
    for(i=0;i<3;i++)for(j=0;j<3;j++){
board[i][j]=i*3+j+1;
    }
    for(i=3;i<9;i++){
board[0][i]=row[sort[stack[stack_pointer][0]][i-3]];
    }
    for(i=1;i<10;i++){
int v=stack[stack_pointer][i];
for(j=0;j<12;j++){
    if(v&(1<<j)){
int u=j/6+1;
int v=j%6+3;
board[u][v]=i;
    }
}
    }
}

int cmp_stack(const void *p,const void *q){
    const short *sp;
    const short *sq;
    int i;
    sp=(const short *)p;
    sq=(const short *)q;
    for(i=10;i<14;i++){
if(sp[i]!=sq[i])
    return sp[i]-sq[i];
    }
    for(i=0;i<10;i++){
if(sp[i]!=sq[i])
    return sp[i]-sq[i];
    }
    return 0;
}

//Classify all 9*9 elements into 2 class.
//The left top 3*3 elements and all elements in first row and first column (total 21 elements)
//  are classified as index elements and their value will be assigned in advance
//1 2 3 * * * * * *
//4 5 6
//7 8 9
//A
//B
//C
//D
//E
//F
//The array is_index gives the index of each index elements
//The array is_value gives the index of each value elements (Non-index elemetns)
//The number of first 9 index elements and the last 6 index elements are filled in advance
//For the other 6 index elements in first column, it is required that for the 2 triple
//  index elements vector, all numbers should be in increasing order and the first element
//  of the first triple is smaller than the first element in the secodn triple
//  That's A<B<C, D<E<F and A<D. So there're total 10 possible filling methods
//  Array row[sort[][]] gives the 10 possible fillings.
//For each number from 1 to 9, it could fill 1 or 2 elements of the first 15 index elements.
//  The index_value[] array record all possible patterns
//  The array class[][] record the pattern id of each number in each possible filling method
//For each class, the later searching algorithm will search for the number in value elements
 
void generate_all_index(){
   int i,j;
   int model[21];
   for(i=0;i<3;i++)for(j=0;j<3;j++){
      is_index[i][j]=i*3+j;
      is_value[i][j]=-1;
   }
   for(i=3;i<9;i++){
      is_index[i][0]=12+i;
      is_value[i][0]=-1;
   }
   for(j=3;j<9;j++){
      is_index[0][j]=6+j;
      is_value[0][j]=-1;
   }
   for(i=1;i<3;i++)for(j=3;j<9;j++){
      is_index[i][j]=-1;
      is_value[i][j]=(i-1)*6+(j-3);
   }
   for(i=3;i<9;i++)for(j=1;j<9;j++){
      is_index[i][j]=-1;
      is_value[i][j]=12+(i-3)*8+j-1;
   }
   for(i=0;i<9;i++)model[i]=i+1;