The 14 March is
Pi Day (π Day) and is a celebration of Pi (3.14 in the
mm.dd date notation). There are a lot of
activities for this day. One of these activities is the Pi Day Sudoku.
Almost exactly one year ago, I blogged about the
Pi Day Sudoku 2009 problem in these two posts:
The problem is an extended version of a Sudoku problem:
Rules: Fill in the grid so that each row, column, and jigsaw region contains 1-9 exactly once and π [Pi] three times.
(Click to enlarge the picture)
Since it was a competition (closed June 1, 2009), I didn't published any models of this problem when blogging about the problem. Now, a year after, it seems to be a good time to publish them. I then implemented two version, one in MiniZinc, and one in Comet::
Both models use the same approach (details differs however). It is the same as for the plain Sudoku, with two exceptions:
- The common approach in plain Sudoku is to use the global constraint
alldifferent for stating that the values in the rows, columns, and regions should be different. Since there should be 3 occurrences of π in each row, column and region, this approach don't work. As mentioned in Solving Pi Day Sudoku 2009 with the global cardinality constraint, my first approach was a home made version of the constraint alldifferent except 0 and Pi but it was too slow. Instead (and via a suggestion of Mikael Zayenz Lagerkvist) I changed to the global constraint global_cardinality, which was much faster.
- The other difference to the standard approach is that the regions are not 3x3 (or MxM), but "jigsaw regions". It was not especially hard to make this representation (though boring to type them in). The constraint for checking the regions are (in MiniZinc):
/\ % the regions
forall(i in 1..num_regions) (
check([x[regions[j,1],regions[j,2]] | j in 1+(i-1)*n..1+((i-1)*n)+11 ])
)
Here I will not publish the solution to the puzzle since I gather that there are a lot of people out there that will solve it by there own. And there is at least one valid solution out there.
Summary
This was a fun problem, especially since I learned some new things by implement the models. As a constraint programming challenge is was quite harder than this year puzzle:
Pi Day 2010:
Rules: Fill in the grid so that each row, column, and block contains 1-9 exactly once. This puzzle only has 18 clues! That is conjectured to be the least number of clues that a unique-solution rotationally symmetric puzzle can have. To celebrate Pi Day, the given clues are the first 18 digits of π = 3.14159265358979323...
[Yes, I have implemented a MiniZinc model for this problem as well; it is a standard Sudoku problem. No, I will not publish the model or a solution until the deadline, June 1, 2010.]
For more about Pi Day Sudoku 2010, see the blog 360:
Pi Day Sudoku is back.
Also, see the following models that implements Killer Sudoku, which use the same approach as Pi Day Sudoku 2009:
The MiniZinc code
Here is the MiniZinc model (
sudoku_pi.mzn), slightly edited.
include "globals.mzn";
int: n = 12;
int: X = 0; % the unknown
int: P = -1; % π (Pi)
predicate check(array[int] of var int: x) =
global_cardinality(x, occ) % :: domain
;
array[1..n, 1..n] of var -1..9: x :: is_output;
array[-1..9] of 0..3: occ = array1d(-1..9, [3, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1]);
array[1..11] of 0..3: occ2 = [3, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1];
% solve satisfy;
solve :: int_search([x[i,j] | i,j in 1..n], first_fail, indomain_min, complete) satisfy;
constraint
% copy the hints
forall(i in 1..n, j in 1..n) (
x[i,j] != 0
/\
if dat[i,j] != X then
x[i,j] = dat[i,j]
else
true
endif
)
/\ % rows
forall(i in 1..n) (
check([x[i,j] | j in 1..n])
)
/\ % columns
forall(j in 1..n) (
check([x[i,j] | i in 1..n])
)
/\ % the regions
forall(i in 1..num_regions) (
check([x[regions[j,1],regions[j,2]] | j in 1+(i-1)*n..1+((i-1)*n)+11 ])
)
;
output
[ show(occ) ++ "\n"] ++
[
if j = 1 then "\n" else " " endif ++
show(x[i,j])
| i in 1..n, j in 1..n
] ++ ["\n"];
% data
array[1..n,1..n] of int: dat = array2d(1..n, 1..n,
[
4,9,7, P,5,X,X,X,X, X,X,X,
X,P,X, 8,X,X,9,6,1, 5,2,X,
X,8,X, 1,X,X,X,P,X, 7,X,X,
X,X,X, X,X,X,X,P,X, 4,X,X,
5,3,9, 6,X,X,X,X,X, X,X,X,
9,4,X, P,P,P,7,X,X, X,X,X,
X,X,X, X,X,6,2,5,P, X,7,4,
X,X,X, X,X,X,X,X,P, P,3,8,
X,7,8, 4,6,9,X,X,X, X,X,X,
X,X,3, X,P,X,X,4,7, 1,6,9,
X,X,4, X,1,X,X,X,6, X,P,X,
X,X,X, X,X,X,X,X,4, X,5,X
]);
% The regions
int: num_regions = 12;
array[1..num_regions*12, 1..2] of int: regions = array2d(1..num_regions*12, 1..2,
[
% Upper left dark green
1,1 , 1,2 , 1,3 ,
2,1 , 2,2 , 2,3 ,
3,1 , 3,2 ,
4,1 , 4,2 ,
5,1 , 5,2 ,
% Upper mid light dark green
1,4 , 1,5 , 1,6 , 1,7 , 1,8 , 1,9 ,
2,4 , 2,5 , 2,6 , 2,7 , 2,8 , 2,9 ,
% Upper right green
1,10 , 1,11 , 1,12 ,
2,10 , 2,11 , 2,12 ,
3,11 , 3,12 ,
4,11 , 4,12 ,
5,11 , 5,12 ,
% Mid upper left "blue"
3,3 , 3,4 , 3,5 , 3,6 ,
4,3 , 4,4 , 4,5 , 4,6 ,
5,3 , 5,4 , 5,5 , 5,6 ,
% Mid Upper right blue
3,7 , 3,8 , 3,9 , 3,10 ,
4,7 , 4,8 , 4,9 , 4,10 ,
5,7 , 5,8 , 5,9 , 5,10 ,
% Mid left green
6,1 , 6,2 , 6,3 ,
7,1 , 7,2 , 7,3 ,
8,1 , 8,2 , 8,3 ,
9,1 , 9,2 , 9,3 ,
% Mid left blue
6,4 , 6,5 ,
7,4 , 7,5 ,
8,4 , 8,5 ,
9,4 , 9,5 ,
10,4 , 10,5 ,
11,4 , 11,5 ,
% Mid mid green
6,6 , 6,7 ,
7,6 , 7,7 ,
8,6 , 8,7 ,
9,6 , 9,7 ,
10,6 , 10,7 ,
11,6 , 11,7 ,
% Mid right blue
6,8 , 6,9 ,
7,8 , 7,9 ,
8,8 , 8,9 ,
9,8 , 9,9 ,
10,8 , 10,9 ,
11,8 , 11,9 ,
% Mid right green
6,10 , 6,11 , 6,12 ,
7,10 , 7,11 , 7,12 ,
8,10 , 8,11 , 8,12 ,
9,10 , 9,11 , 9,12 ,
% Lower left dark green
10,1 , 10,2 , 10,3 ,
11,1 , 11,2 , 11,3 ,
12,1 , 12,2 , 12,3 , 12,4 , 12,5 , 12,6 ,
% Lower right dark green
10,10 , 10,11 , 10,12 ,
11,10 , 11,11 , 11,12 ,
12,7 , 12,8 , 12,9 , 12,10 , 12,11 , 12,12
]);
The Comet code
The Comet model (
sudoku_pi.co) use the same principle as the MiniZinc model. However, the representation of the regions are different where I instead of a matrix use a more object oriented approach with two
tuple for the structures. For some reasons that I have forgot now, I didn't create a function
check in this Comet model, instead stated the
cardinality constraints directly.
import cotfd;
int t0 = System.getCPUTime();
int n = 12;
int P = -1; // Pi
int X = 0; // unknown
range R = -1..9;
set{int} V = {-1,1,2,3,4,5,6,7,8,9};
// regions where 1..9 is alldiff + 3 Pi
tuple Pos {
int row;
int col;
}
tuple Region {
set{Pos} p;
}
int num_regions = 12;
Region regions[1..num_regions] =
[
// Upper left dark green
Region({Pos(1,1),Pos(1,2),Pos(1,3),
Pos(2,1),Pos(2,2),Pos(2,3),
Pos(3,1),Pos(3,2),
Pos(4,1),Pos(4,2),
Pos(5,1), Pos(5,2)}),
// Upper mid light dark green
Region({Pos(1,4), Pos(1,5), Pos(1,6), Pos(1,7), Pos(1,8), Pos(1,9),
Pos(2,4), Pos(2,5), Pos(2,6), Pos(2,7), Pos(2,8), Pos(2,9)}),
// Upper right green
Region({Pos(1,10), Pos(1,11), Pos(1,12),
Pos(2,10), Pos(2,11), Pos(2,12),
Pos(3,11), Pos(3,12),
Pos(4,11), Pos(4,12),
Pos(5,11), Pos(5,12) }),
// Mid upper left "blue"
Region({Pos(3,3), Pos(3,4),Pos(3,5), Pos(3,6),
Pos(4,3), Pos(4,4),Pos(4,5), Pos(4,6),
Pos(5,3), Pos(5,4),Pos(5,5), Pos(5,6)}),
// Mid Upper right blue
Region({Pos(3,7), Pos(3,8), Pos(3,9), Pos(3,10),
Pos(4,7), Pos(4,8), Pos(4,9), Pos(4,10),
Pos(5,7), Pos(5,8), Pos(5,9), Pos(5,10)}),
// Mid left green
Region({Pos(6,1), Pos(6,2),Pos(6,3),
Pos(7,1), Pos(7,2),Pos(7,3),
Pos(8,1), Pos(8,2),Pos(8,3),
Pos(9,1), Pos(9,2),Pos(9,3)}),
// Mid left blue
Region({Pos(6,4),Pos(6,5),
Pos(7,4),Pos(7,5),
Pos(8,4),Pos(8,5),
Pos(9,4),Pos(9,5),
Pos(10,4),Pos(10,5),
Pos(11,4),Pos(11,5)}),
// Mid mid green
Region({Pos(6,6),Pos(6,7),
Pos(7,6),Pos(7,7),
Pos(8,6),Pos(8,7),
Pos(9,6),Pos(9,7),
Pos(10,6),Pos(10,7),
Pos(11,6),Pos(11,7)}),
// Mid right blue
Region({Pos(6,8), Pos(6,9),
Pos(7,8), Pos(7,9),
Pos(8,8), Pos(8,9),
Pos(9,8), Pos(9,9),
Pos(10,8), Pos(10,9),
Pos(11,8), Pos(11,9)}),
// Mid right green
Region({Pos(6,10), Pos(6,11), Pos(6,12),
Pos(7,10), Pos(7,11), Pos(7,12),
Pos(8,10), Pos(8,11), Pos(8,12),
Pos(9,10), Pos(9,11), Pos(9,12)}),
// Lower left dark green
Region({Pos(10,1),Pos(10,2), Pos(10,3),
Pos(11,1),Pos(11,2), Pos(11,3),
Pos(12,1),Pos(12,2), Pos(12,3),Pos(12,4),Pos(12,5), Pos(12,6)}),
// Lower right dark green
Region({Pos(10,10), Pos(10,11),Pos(10,12),
Pos(11,10), Pos(11,11),Pos(11,12),
Pos(12,7),Pos(12,8), Pos(12,9),Pos(12,10),Pos(12,11), Pos(12,12)})
];
// the hints
int data[1..n,1..n] =
[
[4,9,7, P,5,X,X,X,X, X,X,X],
[X,P,X, 8,X,X,9,6,1, 5,2,X],
[X,8,X, 1,X,X,X,P,X, 7,X,X],
[X,X,X, X,X,X,X,P,X, 4,X,X],
[5,3,9, 6,X,X,X,X,X, X,X,X],
[9,4,X, P,P,P,7,X,X, X,X,X],
[X,X,X, X,X,6,2,5,P, X,7,4],
[X,X,X, X,X,X,X,X,P, P,3,8],
[X,7,8, 4,6,9,X,X,X, X,X,X],
[X,X,3, X,P,X,X,4,7, 1,6,9],
[X,X,4, X,1,X,X,X,6, X,P,X],
[X,X,X, X,X,X,X,X,4, X,5,X]
];
Integer num_solutions(0);
Solver m();
var{int} x[1..n,1..n](m, R);
// int occ[-1..9] = [3,0,1,1,1,1,1,1,1,1,1];
var{int} occ[-1..9](m,0..3);
int occ_count[-1..9] = [3,0,1,1,1,1,1,1,1,1,1];
exploreall {
// get the hints
forall(i in 1..n) {
forall(j in 1..n) {
int c = data[i,j];
if (c == P) {
cout << "P";
} else {
cout << data[i,j];
}
if (c != 0) {
m.post(x[i,j] == data[i,j]);
}
cout << " ";
}
cout << endl;
}
forall(i in 1..n, j in 1..n) {
m.post(x[i,j] != 0);
}
forall(i in -1..9) {
m.post(occ[i] == occ_count[i]);
}
// rows
forall(i in 1..n) {
m.post(cardinality(occ, all(j in 1..n) x[i,j]));
}
// columns
forall(j in 1..n) {
m.post(cardinality(occ, all(i in 1..n) x[i,j]));
}
// regions
forall(r in 1..num_regions) {
m.post(cardinality(occ, all(i in regions[r].p) x[i.row,i.col]));
}
} using {
// reversing i and j gives faster solution
forall(i in 1..n, j in 1..n: !x[i,j].bound()) {
tryall(v in V : x[i,j].memberOf(v)) by(v)
m.label(x[i,j], v);
onFailure
m.diff(x[i,j], v);
}
int t1 = System.getCPUTime();
cout << "time: " << (t1-t0) << endl;
cout << "#choices = " << m.getNChoice() << endl;
cout << "#fail = " << m.getNFail() << endl;
cout << "#propag = " << m.getNPropag() << endl;
num_solutions++;
forall(i in 1..n) {
forall(j in 1..n) {
int c = x[i,j];
if (c == P) {
cout << "P";
} else {
cout << x[i,j];
}
cout << " ";
}
cout << endl;
}
}
cout << "\nnum_solutions: " << num_solutions << endl;
cout << endl << endl;
int t1 = System.getCPUTime();
cout << "time: " << (t1-t0) << endl;
cout << "#choices = " << m.getNChoice() << endl;
cout << "#fail = " << m.getNFail() << endl;
cout << "#propag = " << m.getNPropag() << endl;
