C++
What i will present here is an algorithm, illustrated with an example for a 3x3 case. It could theoretically be extended to the NxN case, but that would need a much more powerful computer and/or some ingenious tweaks. I will mention some improvements as I go through.
Before going further, let's note the symmetries of the Sudoku grid, i.e. the transformations which lead to another grid in a trivial way. For block size 3, the symmetries are as follows:
Horizontal symmetry
**The N=3 sudoku is said to consist of 3 "bands" of 3 "rows" each**
permute the three bands: 3! permutations = 6
permute the rows in each band: 3 bands, 3! permutations each =(3!)^3=216
Vertical symmetry
**The N=3 sudoku is said to consist of 3 "stacks" of 3 "columns" each.**
the count is the same as for horizontal.
Note that horizontal and vertical reflections of the grid can be achieved by a combination of these, so they do not need to be counted.
There is one more spatial symmetry to be considered, which is transposing, which is a factor of 2
. This gives the total spatial symmetry of
2*(N!*(N!)^N)^2 = 2*(6*216)^2=3359232 spatial symmetries for the case N=3.
Then there is another, very important symmetry, called relabelling.
Relabelling gives a further (N^2)!=9!=362880 symmetries for the case N=3. So the total
number of symmetries is 362880*3359232=1218998108160.
The total number of solutions cannot be found simply by multiplying the number of symmetry-unique solutions by this number, because there are a number (less than 1%) of automorphic solutions. That means that for these special solutions there is a symmetry operation that maps them to themselves, or multiple symmetry operations that map them to the same other solution.
To estimate the number of solutions, I approach the problem in 4 steps:
1.Fill an array r[362880][12]
with all possible permutations of the numbers 0 to 8. (this is programming, and it is in C, so we are not going to use 1 to 9.) If you are astute you will notice that the second subscript is 12 not 9. This is because, while doing this, bearing in mind that we are going to consider this to be a "row" we also calculate three more integers r[9,10,11] == 1<<a | 1<<b | 1<<c
where 9,10,11 refer to the first, second and third stack and a,b,c are the three numbers present in each stack for that row.
2.Fill an array b
with all possible solutions of a band of 3 rows. To keep this reasonably small, only include those solutions where the top row is 012,345,678. I do this by brute force, by generating all possible middle rows and ANDing r[0][10,11,12]
with r[i][10,11,12]
. Any positive value means there are two identical numbers in the same square and the band is invalid. When there is a valid combination for the first two rows, I search the 3rd (bottom) row with the same technique.
I dimensioned the array as b[2000000][9] but the program only finds 1306368 solutions. I did not know how many there were, so I left the array dimension like that. This is actually only half the possible solutions for a single band (verified on wikipedia), because I only scan the 3rd row from the current value for i
upwards. The remaining half of the solutions can be found trivially by exchanging the 2nd and 3rd rows.
The way the information is stored in array b
is a little confusing at first. instead of using each integer to store the numbers 0..8
found in a given position, here each integer considers one of the numbers 0..8
and indicates in which columns it can be found. thus b[x][7]==100100001
would indicate that for solution x the number 7 is found in columns 0,5 and 8 (from right to left.) The reason for this representation is that we need to generate the rest of the possibilities for the band by relabelling, and this representation makes it convenient to do this.
The two steps above comprise the setup and take about a minute (possibly less if I removed the unnecessary data output. The two steps below are the actual search.)
3 Search randomly for solutions for the first two bands that do not clash (i.e. do not have the same number twice in a given column. We pick a random solution for band 1, assuming always permutation 0, and a random solution for band 2 with a random permutation. A result is normally found in less than 9999 tries (first stage hit rate in the thousands range) and takes a fraction of a second. By permutation, I mean that for the second band we take a solution from b[][] where the first row is always 012,345,678 and relabel it so that any possible sequence of numbers on the first row is possible.
4 When a hit is found in step 3, search for a solution for the third band which does not clash with the other two. We do not want to make just one try, otherwise the processing time for step 3 would be wasted. On the other hand we do not want to put an inordinate amount of effort into this.
Just for fun, last night I did it the dumbest way possible, but it was still interesting (because it nothing for ages, then found large numbers of solutions in bursts.) It took all night to get one datapoint, even with the little hack (!z)
I did to abort the last k
loop as soon as we know this is not a valid solution (which makes it run nearly 9 times faster.) It found 1186585 solutions for the complete grid
after searching all 362880 relabellings of all 1306368 canonical solutions for the last block, a total of 474054819840 possibilities. That's a hit rate of 1 in 400000 for the second stage. I will try again soon with a random search rather than a scan. It should give a reasonable answer in just a few million tries, which should take only a few seconds.
The overall answer should be (362880*(1306368*2))^3*hit rate=8.5E35*hit rate. By back calculating from the number in the question, I expect a hit rate of 1/1.2E14. What I've got so far with my single datapoint is 1/(400000*1000) which is out by a factor of about a million. This could be an anomaly of chance, an error in my program, or an error in my math. I won't know which it is until I run a few more tests.
I'll leave this here for tonight. The text is a bit scrappy, I will tidy it up soon and hopefully add some more results, and maybe a few words on how to make it faster and how to extend the concept to N=4. I don't think I'll be making too many more changes to my program, though :-)
Ah.. the program:
#include "stdafx.h"
#define _CRT_RAND_S
#include <algorithm>
#include <time.h>
unsigned int n[] = { 0,1,2,3,4,5,6,7,8 }, r[362880][12], b[2000000][9],i,j,k,l,u,v,w,x,y,z;
int main () {
//Run through all possible permutations of n[] and load them into r[][]
i=0;
do {
r[i][9] = r[i][10] = r[i][11]=0;
for (l = 0; l < 9; l++){
r[i][l] = n[l];
r[i][9 + l / 3] |= 1 << n[l];
}
if((i+1)%5040==0) printf("%d%d%d %d%d%d %d%d%d %o %o %o %o \n"
,r[i][0],r[i][1],r[i][2],r[i][3],r[i][4],r[i][5],r[i][6],r[i][7],r[i][8],r[i][9],r[i][10],r[i][11],r[i][9]+r[i][10]+r[i][11]);
i++;
} while ( std::next_permutation(n,n+9) );
//Initialise b[][]
for (l = 0; l<2000000; l++) for (k = 0; k<9; k++) b[l][k]=0;
//fill b[][] with all solutions of the first band, where row0 ={0,1,2,3,4,5,6,7,8} and row1<row2
l=0;
for (i = 0; i<362880; i++)
if (!(r[0][9] & r[i][9] | r[0][10] & r[i][10] | r[0][11] & r[i][11])){printf("%d %d \n",i,l);
for (j=i; j<362880;j++)
if(!(r[0][9]&r[j][9] | r[0][10]&r[j][10] | r[0][11]&r[j][11] | r[j][9]&r[i][9] | r[j][10]&r[i][10] | r[j][11]&r[i][11] )){
for (k = 0; k < 9; k++){
b[l][r[0][k]]|=1<<k;
b[l][r[i][k]]|=1<<k;
b[l][r[j][k]]|=1<<k;
}
l++;
}
// printf("%d%d%d %d%d%d %d%d%d %o %o %o %o \n"
// ,r[i][0],r[i][1],r[i][2],r[i][3],r[i][4],r[i][5],r[i][6],r[i][7],r[i][8],r[i][9],r[i][10],r[i][11],r[i][9]+r[i][10]+r[i][11]);
// printf("%d%d%d %d%d%d %d%d%d %o %o %o %o \n"
// ,r[j][0],r[j][1],r[j][2],r[j][3],r[j][4],r[j][5],r[j][6],r[j][7],r[j][8],r[j][9],r[j][10],r[j][11],r[j][9]+r[j][10]+r[j][11]);
// printf("%d %d %o %o %o %o %o %o %o %o %o \n",i,l,b[l][0],b[l][1],b[l][2],b[l][3],b[l][4],b[l][5],b[l][6],b[l][7],b[l][8]);
}
// find a random solution for the first 2 bands
l=0;
do{
rand_s(&u); u /= INT_MIN / -653184; //1st band selection
rand_s(&v); v /= INT_MIN / -181440; //2nd band permutation
rand_s(&w); w /= INT_MIN / -653184; //2nd band selection
z = 0;
for (k = 0; k < 9; k++) z |= b[u][k] & b[w][r[v][k]];
l++;
} while (z);
printf("finished random after %d tries \n",l);
printf("found solution with top band %d permutation 0, and middle band %d permutation %d \n",u,w,v);
getchar();
// scan all possibilities for the last band
l=0;
for (i = 0; i < 362880; i++) for (j = 0; j < 1306368; j++){
z=0;
for(k=0;(k<9)&&(!z);k++) z|= b[u][k] & b[j][r[i][k]] | b[j][r[i][k]] & b[w][r[v][k]];
if (!z){ l++; printf("solution %d : i= %d j=%d",l,i,j); }
}
printf("finished bottom band scan at %d millisec \n", clock()); getchar();
}