Friday, November 16, 2007

Enigma 1469

Here is how Mathematica can be used to solve New Scientist Enigma number 1469 in the 17 November 2007 issue.

WARNING: See the comments for details. There is an error in the solution given below that I need to fix. In working on a fix I discovered a bug in the HamiltonianCycle function in Mathematica's Combinatorica package (this bug has been acknowledged by the author of the Combinatorica package) which destroyed any chance of my fix working fully correctly, although the correct solution to Enigma 1469 could be seen lurking in amongst the erroneous clutter. I will abandon my solution to Enigma 1469 because it has absorbed too much of my time already, and the discovery of the HamiltonianCycle bug can stand as a testimonial to the cleverness of Enigma 1469. I suppose that I ought to leave the rest of this posting intact otherwise the comments would be orphaned, and also there are some useful techniques illustrated in my erroneous solution.

Load the Combinatorica package.

Derive the formula for triangular numbers.
1/2 n (1 + n)

Compute a list of all triangular numbers less than 1000.
triangularnumbers = Table[triangularformula, {n,44}]
{1, 3, 6, 10, 15, 21, 28, 36, 45, 55, 66, 78, 91, 105, 120, 136, 153, 171, 190, 210, 231, 253, 276, 300, 325, 351, 378, 406, 435, 465, 496, 528, 561, 595, 630, 666, 703, 741, 780, 820, 861, 903, 946, 990}

Convert this to the corresponding list of lists of digits.
triangulardigits = IntegerDigits[triangularnumbers]
{{1}, {3}, {6}, {1, 0}, {1, 5}, {2, 1}, {2, 8}, {3, 6}, {4, 5}, {5, 5}, {6, 6}, {7, 8}, {9, 1}, {1, 0, 5}, {1, 2, 0}, {1, 3, 6}, {1, 5, 3}, {1, 7, 1}, {1, 9, 0}, {2, 1, 0}, {2, 3, 1}, {2, 5, 3}, {2, 7, 6}, {3, 0, 0}, {3, 2, 5}, {3, 5, 1}, {3, 7, 8}, {4, 0, 6}, {4, 3, 5}, {4, 6, 5}, {4, 9, 6}, {5, 2, 8}, {5, 6, 1}, {5, 9, 5}, {6, 3, 0}, {6, 6, 6}, {7, 0, 3}, {7, 4, 1}, {7, 8, 0}, {8, 2, 0}, {8, 6, 1}, {9, 0, 3}, {9, 4, 6}, {9, 9, 0}}

Construct the adjacency matrix of allowed adjacencies between triangular numbers.
adjacencies = Outer[(If[Last[#1]==First[#2], 1, 0])&, triangulardigits, triangulardigits, 1];

Convert the adjacency matrix to a set of rules that allow me to add a "tooltip" to each node of the corresponding the graph so that it can be explored by hovering the mouse over the graph: (1) extract the sparse array form of the adjacency matrix, (2) convert the sparse array into a list of rules for which graph nodes are linked by edges, (3) add a tooltip to each graph node.
adjacencies2 = SparseArray[adjacencies];
adjacencies3 = ((adjacencies2//ArrayRules//Most) /. HoldPattern[{i_,j_}->1] :> i->j);
adjacencies4 = adjacencies3 /. x_?IntegerQ :> Tooltip[x, triangularnumbers[[x]]];

Display the graph with tooltips.
GraphPlot[adjacencies4, DirectedEdges->True]
[Output omitted]

Convert the adjacency matrix to a graph. I need this representation in order to use graph manipulation algorithms.
g = FromAdjacencyMatrix[adjacencies, Type->Directed];

Extract the longest cycle in the graph, and extract the edges within this cycle.
longestcyclenodes = ExtractCycles[g]//Sort//Last//Most;
longestcycleedges = Partition[longestcyclenodes,2,1,{1,1}];

Display the graph with the longest cycle highlighted.
GraphPlot[adjacencies4, EdgeRenderingFunction -> (If[MemberQ[longestcycleedges,#2], {Red,Arrow[#1]}, {Opacity[0.3],LightGray,Arrow[#1]}]&)]

Convert the longest cycle to the corresponding a list of triangular numbers, and thence the sequence of digits encountered when going around the cycle.
longestcyclenumbers = triangularnumbers[[longestcyclenodes]];
longestcycledigits = longestcyclenumbers//IntegerDigits//Flatten

[Output omitted]

Compute the number of digits encountered when going around the longest cycle.
[Solution omitted]


Roman said...

good afternoon,

I solved the enigma in Matlab (don't have mathematica), but I'm not 100% sure if I didn't miss something along the way. The solution I got was 21 digits, and there are actually quite a few solutions with 21 digits. Do you reach the same solution?

Thanks very much!

Roman said...

I found a small oversight in the code when I had another look just now, I now come to a solution of 30 digits. Is this more like it? Thanks very much for your help.

Stephen Luttrell said...

There is only one longest cycle, and it is shorter than what you have found.

Roman said...

but, for instance, this works:
{1,15,55,595,528,861,153,3,351,105,561 ,171}
A sequence of 12 numbers (30 digits), and there are 287 more solutions with 30 digits like this one. I don't think that the numbers have to be monotonically ascending, do they?

Stephen Luttrell said...

Yes, you are right. Thanks for pointing that out.

I have now looked more closely at the documentation for ExtractCycles, and it says:

ExtractCycles[g] gives a maximal list of edge-disjoint cycles in graph g.

I had carelessly overlooked the "edge-disjoint" part of this definition, so the cycles that I extracted are not guaranteed to include the absolutely longest cycle, and this case is clearly one where ExtractCycles does not give the result I had hoped that it would.

In Mathematica there doesn't seem to be a single function that computes the longest cycle (the graph "circumference"), but there is one for computing the shortest cycle (the graph "girth"). I will therefore have to write a "patch" to fix the Mathematica solution as soon as I have time.

Again, thanks for pointing this out.

Stephen Luttrell said...

I have now had a chance to look at this again, and I have a solution based on an approach where I carefully remove nodes from the graph until the function HamiltonianCycle can find a complete cycle in the remaining graph.

However, I found what appears to be a bug in HamiltonianCycle, so I will hold off posting this solution until I can be sure that it will work OK.

Stephen Luttrell said...

I posted a comment entitled Possible bug in HamiltonianCycle to MathGroup, and got a response from Daniel Lichtblau of Wolfram Research who explained that the behaviour that I had observed was an "undocumented feature" of the implementation of HamiltonianCycle.

Merde, merde, and thrice merde! This is not what I wanted to hear.

I have now emailed the author of the Combinatorica package which contains the implementation of the HamiltonianCycle function, to ask for confirmation that the behaviour of HamiltonianCycle that I have observed is what he intended, i.e. there is a part of the algorithm that treats directed graphs as if they are undirected. I'm sure that the answer to this question should be "no", and I would like to see what is going to be done about it.

Stephen Luttrell said...

...and the response from the author of the Combinatorica package re the bug (or is it an "undocumented feature") in HamiltonianCycle is:

"It does seem a fix worth making."

At least the problem has now been reported and acknowledged, but I guess that its solution will not be immediately available.

But that leaves me with a problem, because my use of Mathematica to solve New Scientist Enigma problems is designed to show off how you can solve non-trivial problems using only a few lines of Mathematica code.

Maybe I could rework my solution to Enigma 1469 to avoid these problems, but I think that I will call a halt to this solution because it has served its (rather unexpected) purpose of highlighting a bug (OK, "undocumented feature") in HamiltonianCycle.

My thanks go to commenter "roman" for firing a shot over my bow which led to this discovery.

Roman said...

Glad to be of help. I've never worked with Mathematica, I just wrote a short impromptu program in Matlab which is probably not very robust or fast, and based on some assumptions. I can send you the code if you like... I'm curious to see what the new scientist deems the solution to be.

Stephen Luttrell said...

Using the buggy version of HamiltonianCycle I obtained a 30-digit solution like the one you found, and also a 32-digit solution. It was when I looked closely at the 32-digit solution that I realised that it broke the rules a little, and on further investigation I traced this to the inability of HamiltonianCycle to correctly compute cycles in directed graphs.

Some of my Matlab-using colleagues will find this all to be hugely amusing, because I am less than polite about Matlab's capabilities compared to those of Mathematica!

There is no need to send me your Matlab solution. My reasons for solving these New Scientist Enigmas are described here.