6670903752021072936960
据说在2003年这个数字被计算出来了。
我们也可以计算一下呦！
方法呢？当然是记忆化搜索。
为了简化状态，我们每个3*3一起填。
记录每个列、每个行当前剩下的数字的集合。开一个数组来存储搜索到当前状态的方法数。
当然，状态太多了，数组是开不下的。于是，我们要充分利用对称性来减少状态数量。
第一部分的对称性来自数字的对称。1~9是可以相互替换的。这样，我们可以把状态的个数减少很多。
第二部分的对称来自行和列的对称。注意，我们把第1、2、3行相互交换是无所谓的，第1、2、3列之间也是可以整列交换的。于是，状态的个数还可以减少6的4次方。
最后的对称来自整3行的对称。在填完第3~6行第1~3列之后，我们可以把4~6列和7~9列整体交换。这样，状态个数再次除以2。
状态总数缩小到了几百万。可以存下了。
写代码出来是一个不容易的工作。
经过很努力的优化，我的程序可以在4分钟内出结果了。
#include <cstdio> #include <algorithm> #include <string.h> #include <ctime> using namespace std; const int Max3=280; const int NMax=Max3*Max3*Max3; const long long Mod=0; long long dp1[NMax],dp2[NMax]; int repr1[NMax][4],repr2[NMax][4]; int A[12]; long long cur; int bitcount[512]; int comb[280]; const int AMax=50; int Ans[50]; int idcnt=0; int getid(const int *a,int y){ idcnt++; int perm[9],n=0; for (int i=0;i<9;i++)perm[i]=-1; int b[12]; for (int i=0;i<12;i++){ if (bitcount[a[i]]<=3)b[i]=a[i]; else b[i]=a[i]^511; } int val[9]; for (int i=0;i<9;i++)val[i]=i; for (int i=0;i<12;i++){ int s=b[i]; if (s){ int j=__builtin_ctz(s); val[j]|=((1<<20)>>i); s&=s-1; j=__builtin_ctz(s); val[j]|=((1<<20)>>i); s&=s-1; j=__builtin_ctz(s); val[j]|=((1<<20)>>i); //for (int j=0;j<9;j++)if (s&(1<<j))val[j]|=(1<<(20-i)); } } if (b[0]){ for (int i=0;i<3;i++){ int s=b[i]; int ps[3]; ps[0]=__builtin_ctz(s); s&=s-1; ps[1]=__builtin_ctz(s); s&=s-1; ps[2]=__builtin_ctz(s); //for (int j=0;j<9;j++) //if (s&(1<<j))ps[pc++]=j; if (val[ps[0]]<val[ps[1]])swap(ps[0],ps[1]); if (val[ps[0]]<val[ps[2]])swap(ps[0],ps[2]); if (val[ps[1]]<val[ps[2]])swap(ps[1],ps[2]); perm[ps[0]]=n; perm[ps[1]]=n+1; perm[ps[2]]=n+2; n+=3; } for (int i=3;i<12;i++){ int t=0,s=b[i]; if (s){ int j=__builtin_ctz(s); t|=(1<<perm[j]); s&=s-1; j=__builtin_ctz(s); t|=(1<<perm[j]); s&=s-1; j=__builtin_ctz(s); t|=(1<<perm[j]); } //for (int j=0;j<9;j++)if (s&(1<<j)){ //if (perm[j]==-1){ //perm[j]=n++; //} //t|=(1<<perm[j]); //} b[i]=t; } }else{ for (int i=0;i<12;i++){ int t=0,s=b[i]; if (s){ int j=__builtin_ctz(s); if (perm[j]==-1)perm[j]=n++; t|=(1<<perm[j]); s&=s-1; j=__builtin_ctz(s); if (perm[j]==-1)perm[j]=n++; t|=(1<<perm[j]); s&=s-1; j=__builtin_ctz(s); if (perm[j]==-1)perm[j]=n++; t|=(1<<perm[j]); } //for (int j=0;j<9;j++)if (s&(1<<j)){ //if (perm[j]==-1){ //perm[j]=n++; //} //t|=(1<<perm[j]); //} b[i]=t; } } int r=0; for (int i=3;i<12;i+=3){ if (b[i]>b[i+1])swap(b[i],b[i+1]); if (b[i]>b[i+2])swap(b[i],b[i+2]); if (b[i+1]>b[i+2])swap(b[i+1],b[i+2]); } if (y==0){ if (b[6]>b[9]){ swap(b[6],b[9]); swap(b[7],b[10]); swap(b[8],b[11]); } } for (int i=3;i<12;i+=3){ int u=b[i],v=b[i+1],w=b[i+2]; if (u==0 && v==0 && w==0)r=r*Max3; else{ int t=(u<<18)|(v<<9)|w; int p=0,q=Max3; while (q-p>1){ int m=(q+p)>>1; if (comb[m]<=t)p=m; else q=m; } r=r*Max3+p; } } for (int i=0;i<4;i++){ int t=0; for (int j=0;j<3;j++) t|=a[i*3+j]<<(j*9); repr1[r][i]=t; } return r; } void dfs(int j,int x,int y,int u){ if (y==3)x++,y=0; if (x==3){ int B[12]; for (int i=0;i<j*3;i++)B[i]=A[i]; for (int i=j*3;i<j*3+3;i++)B[i]=A[i+3]; for (int i=j*3+3;i<j*3+6;i++)B[i]=A[i-3]; for (int i=j*3+6;i<12;i++)B[i]=A[i]; int k=getid(B,j); dp1[k]+=cur; if (Mod && dp1[k]>=Mod) dp1[k]-=Mod; }else{ for (int i=0;i<9;i++)if (!((1<<i)&u) && !((1<<i)&(A[j*3+x]|A[j*3+3+y]))){ A[j*3+x]|=(1<<i); A[j*3+3+y]|=(1<<i); dfs(j,x,y+1,u|(1<<i)); A[j*3+x]^=(1<<i); A[j*3+3+y]^=(1<<i); } } } int ways[3],n1; void dfs1(int a,int u){ if (a==3){ if (ways[0]<ways[1] && ways[1]<ways[2]){ comb[n1]=(ways[0]<<18)|(ways[1]<<9)|ways[2]; n1++; } return; } for (int i=0;i<9;i++)if (!((1<<i)&u)){ for (int j=0;j<i;j++)if (!((1<<j)&u)){ for (int k=0;k<j;k++)if (!((1<<k)&u)){ ways[a]=(1<<i)|(1<<j)|(1<<k); dfs1(a+1,u|ways[a]); } } } } void showret(){ int l=AMax-1; while (l && Ans[l]==0)l--; for (int i=l;i>=0;i--)printf("%d",Ans[i]); puts(""); } int main() { bitcount[0]=0; for (int i=1;i<512;i++)bitcount[i]=bitcount[i>>1]+(i&1); n1=0; dfs1(0,0); sort(comb,comb+n1); memset(dp1,0,sizeof(dp1)); for (int i=0;i<12;i++)A[i]=0; dp1[getid(A,-1)]=1; memset(Ans,0,sizeof(Ans)); for (int i=0;i<3;i++){ for (int j=0;j<3;j++){ printf("i=%d j=%d\n",i,j); memcpy(dp2,dp1,sizeof(dp2)); memcpy(repr2,repr1,sizeof(repr2)); memset(dp1,0,sizeof(dp1)); int cnt=0; for (int k=0;k<NMax;k++)if (dp2[k]){ cur=dp2[k]; if (i==0 && j==1)cur/=362880; for (int l=0;l<12;l++) A[l]=(repr2[k][l/3]>>(9*(l%3)))&511; if (j==0){ for (int l=11;l>=3;l--)A[l]=A[l-3]; A[0]=A[1]=A[2]=0; } if (k%1000==0){ printf("on %d got ",k); for (int i=0;i<12;i++){ printf("("); for (int j=0;j<9;j++) if (A[i]&(1<<j))printf("%d",j); printf(")"); } printf(" with %lld\n",cur); } dfs(j,0,0,0); cnt++; } printf("cnt=%d\n",cnt); } } for (int i=0;i<362880;i++){ long long t=dp1[0]; for (int j=0;j<AMax;j++){ t+=Ans[j]; Ans[j]=t%10; t/=10; } } showret(); return 0; }
