How many Sudoku puzzles exist?

Question

This is not a Sudoku solver, nor a Sudoku checker.

Your challenge is to write a function or script that, given as input the "block" size of a 2D Sudoku puzzle (which is 3 for the classic 9x9 board, 4 for a 16x16 board, etc.) will compute an approximation of the number of distinct puzzles (solutions) that exist for that size.

For example, given input 3, your program should print an approximation, to the desired precision, of the number 6,670,903,752,021,072,936,960 which is the known number of distinct 9x9 Sudoku puzzles, or 5,472,730,538 when taking into account the various symmetries. Your solution should state whether symmetries are counted or ignored.

"Desired precision" is left undefined: your program might run for a given time and then output its result, or compute it up to a given number of significant digits, or even run forever, printing better and better approximations. The point is that it should be possible to make it compute the result to any required precision, in a finite time. (So "42" is not an acceptable answer.) Restricting the precision of your result to the available machine floats is acceptable.

No accessing online resources, no storing the source code in the file name, etc.

---

PS: I know this is a Hard Problem (NP-complete if I'm not mistaken.) But this question is only asking for an approximate, statistical solution. For example you can try random configurations that satisfy one (or better two) constraints, compute how many of those exist and then check how often you get a puzzle that satisfies all three constraints. This will work in a decent time for small sizes (certainly for size=3 and possibly 4) but the algorithm should be generic enough to work for any size.

The best algorithm wins.

---

PS2: I changed from code-golf to code-challenge to better reflect the difficulty of the problem and encourage smarter solutions, over dumb but well-golfed ones. But since apparently "best algorithm" is unclear, let me try to define it properly.

Given enough time and disregarding constant factors (including CPU and intepreter speed), or equivalently, considering their asymptotical behaviour, which solution would converge to the exact result the fastest?

Answer 1

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