## Using Maple 15 to Solve a Peg Board Puzzle Game

Maple 15 In high school I was briefly fascinated by a triangular "jump all but one" game, commonly found at Cracker Barrel restaurants.  The basic premise is that any peg can "jump" over an adjacent peg to occupy the empty hole next to the jumped peg.  The jumped peg is then removed.  The goal is to continue jumping pegs until there is only one left.

The instructions on the face of the Cracker Barrel version of this game say, "LEAVE ONLY ONE -- YOU'RE A GENIUS".  Wanting to claim the right to call myself a genius, unlike ordinary kids, who might just play the game a few times, I sat down on my Turbo-XT and started writing BASIC code.  The algorithm I came up with ran a bit slow, so I directed output to my printer and let it run over night.  In the morning the program was still chugging along.  I advanced the paper feed on the dot-matrix lineprinter -- the kind that used continuous feed paper with perforated edges and holes on each side.   Into view came 3 solutions represented by a string of numbers.   A quick check verified that I was now a genius.

With the addition of Maple 15's new parallel multi-process features, and beefed up Grid Computing package, I was recently thinking about big examples that could show off how fast Maple can run on modern multi-core computers.  My mind turned back to this toy game, and I wondered just how long it would take to find a solution in Maple?

It turns out that a brute force algorithm for finding solutions is fairly easy to write.  Here, I generalized it to a triangle with N-pegs on each side.

Peg := proc( N )
local peg, i, j, r, A;

peg := proc( A, N, path, sol )
local i, j, p;

p := [op(path),A];
for i from 1 to N do
for j from 1 to N-i+1 do
if A[i,j] then
if i+2 <= N-j+1 and A[i+1,j] and not A[i+2,j] then
peg( subsop([i,j]=false,[i+1,j]=false,[i+2,j]=true,A), N, p, sol );
end if;
if i-2 > 0 and A[i-1,j] and not A[i-2,j] then
peg( subsop([i,j]=false,[i-1,j]=false,[i-2,j]=true,A), N, p, sol );
end if;

if i-2 > 0 and j+2 <= N-i+3 and A[i-1,j+1] and not A[i-2,j+2] then
peg( subsop([i,j]=false,[i-1,j+1]=false,[i-2,j+2]=true,A), N, p, sol );
end if;
if j-2 > 0 and i+2 <= N-j+3 and A[i+1,j-1] and not A[i+2,j-2] then
peg( subsop([i,j]=false,[i+1,j-1]=false,[i+2,j-2]=true,A), N, p, sol );
end if;

if j+2 <= N-i+1 and A[i,j+1] and not A[i,j+2] then
peg( subsop([i,j]=false,[i,j+1]=false,[i,j+2]=true,A), N, p, sol );
end if;
if j-2 > 0 and A[i,j-1] and not A[i,j-2] then
peg( subsop([i,j]=false,[i,j-1]=false,[i,j-2]=true,A), N, p, sol );
end if;

end if;
end do;
end do;
if numelems(p) = N*(N+1)/2-1 then
sol(numelems(sol)+1) := p;
end if;
end proc:

A := [seq([seq(true,j=1..N-i+1)],i=1..N)]; #fill all holes
A[1,1] := false; #pick first spot as empty
r := Array(1..0);
peg(A,N,[],r);
r;
end proc: The Cracker Barrel version had 5 pegs on each edge.   Running this in serial, Maple finds 29760 solutions in just under a minute.  That's the number of solutions when starting with a corner peg empty. > t := time[real]():
 > R := Peg(5);
 > time[real]() - t;  (1) Each entry in the array,  R, stores an encoding of the solution.  So you don't have to decode it, below is a simple procedure to animate any given solution.  Click the plot generated and press the play button to see how to proceed through a given solution. PegPlot := proc( moves, N )
local pegframe, i;

pegframe := proc( A, N )
local i, j, p1, p2;

#plot board
p1 := plots:-pointplot( [
seq(seq([(i-N)/2+j,N-i],j=1..N-i+1),i=1..N)
],symbol=circle,symbolsize=31,axes=none);

#plot pegs
p2 := plots:-pointplot( [
seq(seq(`if`(A[i,j],[(i-N)/2+j,N-i],NULL),j=1..N-i+1),i=1..N)
],symbol=solidcircle,symbolsize=30,color=blue,axes=none);

plots[display]([p2,p1]);
end proc;

plots:-display([seq(pegframe(moves[i],N),i=1..numelems(moves))],insequence=true);

end proc:

> PegPlot(R,5);  The algorithm for finding solutions is brute-force.  It searches through every possible combination of moves in a recursive way.  The outlook for parallelizing the code seemed pretty good.  An initial version which just split each of the N  starting points off onto another process immediately cut the speed in half.  This was a pretty good result for doing almost no work.  Looking at the balance of work being done, I saw that this split wasn't balanced.  One node was doing half the work, while the other three cores were long finished.  After some tweaking I came up with the following: ParallelPeg := proc( N )

local peg, i, j, r, holes, nholes, me, numnodes, A;

peg := proc( A, N, path, sol, lev )
local i, j, p, k1, k2, h, split;

p := [op(path),A];

if lev = 6 then
me := Grid:-MyNode();
numnodes := Grid:-NumNodes();
nholes := numelems(holes);
split := holes[(me)*floor(nholes/numnodes) + 1
.. `if`(me=numnodes-1,nholes,(me+1)*floor(nholes/numnodes))];
#printf("Node %d level %d split: %La\n",Grid:-MyNode(),lev,split); >
else
split := holes;
end if;

for h in split do
i,j := op(h);
if A[i,j] then
if i+2 <= N-j+1 and A[i+1,j] and not A[i+2,j] then
peg( subsop([i,j]=false,[i+1,j]=false,[i+2,j]=true,A), N, p, sol, lev+1 );
end if;
if i-2 > 0 and A[i-1,j] and not A[i-2,j] then
peg( subsop([i,j]=false,[i-1,j]=false,[i-2,j]=true,A), N, p, sol, lev+1 );
end if;

if i-2 > 0 and j+2 <= N-i+3 and A[i-1,j+1] and not A[i-2,j+2] then
peg( subsop([i,j]=false,[i-1,j+1]=false,[i-2,j+2]=true,A), N, p, sol, lev+1 );
end if;
if j-2 > 0 and i+2 <= N-j+3 and A[i+1,j-1] and not A[i+2,j-2] then
peg( subsop([i,j]=false,[i+1,j-1]=false,[i+2,j-2]=true,A), N, p, sol, lev+1 );
end if;

if j+2 <= N-i+1 and A[i,j+1] and not A[i,j+2] then
peg( subsop([i,j]=false,[i,j+1]=false,[i,j+2]=true,A), N, p, sol, lev+1 );
end if;
if j-2 > 0 and A[i,j-1] and not A[i,j-2] then
peg( subsop([i,j]=false,[i,j-1]=false,[i,j-2]=true,A), N, p, sol, lev+1 );
end if;

end if;
end do;
if numelems(p) = N*(N+1)/2-1 then
sol(numelems(sol)+1) := p;
end if;
NULL;
end proc:

printf("Node %d called with N = %a\n",Grid:-MyNode(),N);
holes := [seq(seq([i,j],j=1..N-i+1),i=1..N)];

A := [seq([seq(true,j=1..N-i+1)],i=1..N)]; #fill all holes
A[1,1] := false; #pick first spot as empty
r := Array(1..0);
peg(A,N,[],r,1);

printf("Node %d found %d solutions\n",Grid:-MyNode(),numelems(r));
if Grid:-MyNode() = 0 then
else
Grid:-Send(0,r);
end if;
r;
end proc:

 > t := time[real]():
 > R2 := Grid:-Launch(ParallelPeg,5,numnodes=5);
 > time[real]()-t;
 Node 3 called with N = 5 Node 2 called with N = 5 Node 1 called with N = 5 Node 4 called with N = 5 Node 0 called with N = 5 Node 3 found 3456 solutions Node 2 found 2925 solutions Node 1 found 7309 solutions Node 0 found 7969 solutions Node 4 found 8101 solutions  (2)

I ran this using 5 processes computing on a 4-core cpu.  The balance of work was better, with no node doing more than about 30% of the work for the 5-peg triangle.  The run time went from 47 seconds down to 19 seconds.  The important thing to note was that the code stayed mostly the same.

I haven't had the patience to run a 6-peg triangle using this code yet.  Maybe I should tie the output to my printer?  