Computer Assisted Mathematics: Tools and Tactics for Solving Hard Problems

Daniel Lichtblau
danl@wolfram.com

Wolfram Research, Inc.
October, 2006

Abstract

In this talk I will present several problems that have caught my attention over the past few years. We will go over Mathematica formulations and solutions. Along the way we will meet with a branch-and-bound loop in its natural habitat, some rampaging Gröbner bases, a couple of tamed logic puzzles, and at least a dozen wild beasts (well...would you believe... a Mathematica developer wearing Halloween fangs?). You'll laugh. You'll cry. You'll write your holiday cards (be sure to spell my name correctly).

As the purpose is to illustrate a few of the many ways in which Mathematica can be used to advantage in tackling difficult problems, we will go into a bit of detail in selected examples. Do not let this deter you; there will be no exam, and it is the methods, not the problems, that are of importance. The examples are culled from problems I have seen on Usenet groups (primarily MathGroup), in articles, or have been asked in person.

() A simple math problem

At a workshop in 2004 I was posed the following: Find a positive integer such that, when multiplied by 9, you get the same result as by moving the least significant digit to the most significant. The person asking is editor-in-chief of the Journal of Symbolic Computation. He told me the (smallest) result had 44 digits, and it had taken him substantial computational effort to find it.

We code this as a simple integer programming problem, directly invoking the Mathematica function Reduce in order to solve it.

timesBy9OrRotateDigits[j_]:=Module[{n,m,k,res},res={ToRules[Reduce[{n==10*m+k,9*n==m+10^j*k,Element[{n,m,k},Integers],m≥0,1≤k≤9,n≥10^j},{n,m,k}]]};
If[res==={},{},{j+1,First[n/.res]}]]

We confirm below that the smallest solution has 44 digits. It can be shown to be the repeating block of the fraction 9/89.

Timing[Table[timesBy9OrRotateDigits[j],{j,45}]/.{}->Sequence[]]

His comment: "That's correct! But you cheated, using Mathematica:-)"

He, alas, had used another program...

() Algebraic replacement of expressions

This is a topic of some apparent interest insofar as it shows up with regularity on the Mathematica user's Usenet group comp.softsys.math.mathematica (aka MathGroup). One has a perhaps complicated expression, and wishes to replace subexpressions using perhaps a new "variable" to encapsulate them. The issue is that they might not be amenable to the familiar syntactic replacement rules of Mathematica.

For example, one might wish to replace (a^2) b with a new variable, and likewise have (a^4) (b^2) replaced by the square of that variable, etc. The particular example below comes from a recent MathGroup question [algebraic replacement].

expr=-((1+2*n)*((a^4*k^2+a^2*(-1+k^2*(q-z)^2)+2*(q-z)^2)*Cos[k*Sqrt[a^2+(q-z)^2]]-k*(a^2-2*(q-z)^2)*Sqrt[a^2+(q-z)^2]*Sin[k*Sqrt[a^2+(q-z)^2]])*Sin[((1+2*n)*Pi*z)/L])/(8*Pi*w*(a^2+(q-z)^2)^(5/2));

The desire was to replace (Sqrt[(a^2)+(((q-z))^2)]) with a new variable R.

Algebraic replacement...

I have written code to address this sort of question on a few occasions, and it seems to get a bit more complicated each time as the mission of the replacement gets more elaborate. Below is the current form I use. The key is to use the function PolynomialReduce to do a sort of "generalized division" in order to get an algebraic replacement.

replacementFunction[expr_,rep_,vars_]:=With[{num=Numerator[expr],den=Denominator[expr],hed=Head[expr]},If[PolynomialQ[num,vars]&&PolynomialQ[den,vars],PolynomialReduce[num,rep,vars][[2]]/PolynomialReduce[den,rep,vars][[2]],If[Head[hed]===Symbol&&MemberQ[Attributes[hed],NumericFunction],Map[replacementFunction[#,rep,vars]&,expr],expr]]
]

There is a price to pay for getting the power of algebraic reduction. We really need to use polynomials in our replacements, and that means we cannot work directly with the desired radical. This is unfortunate in that it will leave us with radicals of powers, but we can fix that using e.g. PowerExpand or (as a post-processing step) a syntactic replacement rule.

replacementFunction[expr,a^2+(q-z)^2-R^2,{a,q,z}]/.(R^aa_)^(bb_):>R^(aa*bb)

Had we multiple replacements, perhaps with interdependencies, then we might want to form a Gröbner basis prior to use of PolynomialReduce. That's another topic.

() Gröbner bases

I do not propose to explain in any detail what these are. For our purposes they are to be regarded as an algebraic form of "rewriting rules". That is to say, we have polynomials and we want to rewrite according to some (polynomial) equivalence relations, perhaps eliminating some variables along the way.

Various tactics are used to extend this idea to handling rational functions and radicals. We will see examples of this in later sections.

() Unraveling a parametrized nested radical

This next problem was posed last year on MathGroup [parametrized root].

expr=(-1+(Sqrt[2])-(Sqrt[(1+2 (Sqrt[2])-s) (-1+s)])+s) (-1-(Sqrt[2])+(Sqrt[2]) s-(Sqrt[-1-2 (Sqrt[2])+2 (2+(Sqrt[2])-s) s]))-(-1-(Sqrt[2])+(Sqrt[2]) s+(Sqrt[-1-2 (Sqrt[2])+2 (2+(Sqrt[2])-s) s])) (-2-2 (Sqrt[2])+2 s+2 √(1-(1/4) (((-1-(Sqrt[2])+(Sqrt[(1+2 (Sqrt[2])-s) (-1+s)])+s))^2)));

This is an algebraic function and we wish to find the exact solution where it vanishes around s≈2.707. A brief look shows that, at the root, it "goes complex" as the parameter s increases.

Plot[{Re[expr],Im[expr]}, {s,2,3}]

Nested radical...

One way is to create a Gröbner basis (in effect, describing the vanishing set expr==0) and let it do all the hard work of sorting out the algebraic dependencies.

gb=GroebnerBasis[expr,s];

In general we might hope there is a "reasonable" element in the basis, in particular one that is explicitly polynomial in s. And indeed it turns out the first element satisfies that description.

poly=First[gb]

Nested radical...

We now want to figure out a minimal polynomial component that vanishes near 2.707. First we find the vanishing point more precisely.

acc=100;
sroot=FindRoot[expr==0,{s,2.7},AccuracyGoal->acc,PrecisionGoal->acc,WorkingPrecision->(1.2*acc)];

Now we factor our polynomial and find the factor that vanishes at sroot.

fax=Drop[Map[First,FactorList[poly,Extension->Automatic]],1];
minpoly=First[Select[fax,Chop[#/.sroot,10^(-acc)]===0&]]

We now use Solve to find our value as a Root object.

First[Select[s/.Solve[minpoly==0,s],Abs[#-(s/.sroot)]<10^(-acc/5)&]]

Nested radical...

We now use some number theory functionality, buried under the hood and based on lattice reduction, to deduce the same result from the high rpecision approximating value. We assume the result will have a minimal polynomial of degree no larger than 20. This formerly required an add-on package.

RootApproximant[s/.sroot,20]

It turns out that for this particular example, if we let RootApproximant guess the degree it will come up with the right result.

RootApproximant[s/.sroot]

() Implicitizing a rational surface

We present an example of a parametrized rational surface, wherein we use a Gröbner basis to find the implicit form. The example comes from fairly recent literature [parametrized surface]. Why do we want to implicitize it? Well, there might be many reasons, but one is that it actually is a bit easier to visualize because the plot routines seem to handle that form better. I think this is because the curve where the parametric form denominator vanishes is problematic, and moreover there are several base points (where numerators ALSO vanish), near which we are likely to confuse the plotting routines.

We regard the parametric form as a system of polynomials x,y,z each parametrized by s,t but we also have a common denominator w likewise parametrized.

w=(t^3)+(t^2)-t+(s^2) t-1-s+(s^2)+(s^3);
polys={2 (t^3)+4 (t^2)+2 t+4 s t+(s^2) t+2+3 s+(s^2),-2 s (t^2)-2 t-s t+2+s-2 (s^2)-(s^3),2 (t^2)-3 s (t^2)-2 t-3 s t-2 (s^2) t-2 s-3 (s^2)-(s^3)};

Implicitizing...

So let's have a look at this thing. I confess it is a bit tricky to get the parametric plot to do anything nice with it, and I am using a variant with settings rigged by Michael Trott MMAL ("Mathematica Master At Large").

L=2.2;
ParametricPlot3D[Evaluate[polys/w],{s,-L,L},{t,-L,L},Mesh->False,PlotPoints->40,
PlotRange->5,Exclusions->{w==0, Offset -> 0.01, MaxRecursion -> 3},
ColorFunction -> (Hue[#5]&)]

Implicitizing...

Pretty, but it's a bit difficult to figure out. We will implicitize. How? What one might try is to form a Gröbner basis of polynomials {w x-x(s,t),w y-y(s,t),w z-z(s,t)} where w is an explicit polynomial in the parameters, (x,y) are distinct variables, and (x(s,t),y(s,t)) are again explicit polynomials. The goal would be to eliminate the parameters (s,t). This alone does not work due to presence of base points. One repairs the situation by adding a polynomial w r-1 in a new "reciprocal" variable r. This relation forces the denominator to be nonzero. We add this new variable to the list of those we seek to eliminate.

vars={x,y,z};
elims={s,t,r};
Timing[implicit=First[GroebnerBasis[Append[w*vars-polys,w*r-1],vars,elims,MonomialOrder->EliminationOrder]]]

Implicitizing...

Nothing to it. Now we do a contour plot.

l=20;
ContourPlot3D[Evaluate[imp==0],{x, -l, l},{y,-l,l},{z,-l,l}]

Much nicer, I think.

() Puzzle: Johnny's ideal woman

This puzzle appeared in a university publication in 1998. It was asked in MathGroup in 2003 [Johnny's ideal woman].

Johnny's ideal woman is brunette, blue eyed, slender, and tall. He knows four women: Adele, Betty, Carol, and Doris. Only one of the four women has all four characteristics that Johnny requires.
1. Only three of the women are both brunette and slender.
2. Only two of the women are both brunette and tall.
3. Only two of the women are both slender and tall.
4. Only one of the women is both brunette and blue eyed.
5. Adele and Betty have the same color eyes.
6. Betty and Carol have the same color hair.
7. Carol and Doris have different builds.
8. Doris and Adele are the same height.

At IMS 2004 I showed a method that had been presented in comp.soft-sys.math.mathematica by Lars Rassmusson. It involves setting up and solving a system of equations (which, under the hood, uses a Gröbner basis).

Here I will show a logic programming approach. It is simple, yet effective. Also of interest is the contrast to a method we will cover later. The method below is similar to that of a couple of responses that appeared in the forum. I think it is generally more effective, however, in that it makes use of intermediate pruning steps so that set sizes do not get very large. Another place this approach was used to advantage was by Michael Trott (MMAL) in solving the notorious "Who owns the zebra" puzzle [Trott TMJ article].

Johnny's ideal woman...

The idea is to enlarge, via enumeration of possibilities, and then prune the search space iteratively as new variables are taken into account. The enumeration steps consist of adding all possible values of a new variable to each partial solution. These are then pruned using available constraints in an effort to keep the size of the partial solution set reasonable. We will illustrate the method on this example. We begin by showing the power set of eye combinations that can go with the four women (we will not show this for the remaining features).

gals={Adele,Betty,Carol,Doris};
features={blueeye,slender,brunette,tall};

eyes=Flatten[Outer[List,Sequence@@Table[{blueeye,Not[blueeye]},{4}]],3]

Johnny's ideal woman...

We form the lists of possibilities combining women with eye color.

galsAndEyes=Flatten[Outer[Thread[And[##]]&,{gals},eyes,1],1]

Now we use the relevant constraint to prune this list.

galsAndEyes=Cases[galsAndEyes,{AorB_&&a_,___,BorA_&&a_,___}/;MatchQ[AorB,Adele|Betty]&&MatchQ[BorA,Adele|Betty]]

Johnny's ideal woman...

We add an attribute for girth.

widths=Flatten[Outer[List,Sequence@@Table[{slender,Not[slender]},{4}]],3];
galsAndEyesAndWidth=Flatten[Outer[Thread[And[##]]&,galsAndEyes,widths,1],1];

This time we can prune with more than one constraint.

galsAndEyesAndWidth=Cases[galsAndEyesAndWidth,{___,And[CorD_,_,a_],___,And[DorC_,_,b_],___}/;MatchQ[CorD,Carol|Doris]&&MatchQ[DorC,Carol|Doris]&&b==Not[a]];
galsAndEyesAndWidth=Select[galsAndEyesAndWidth,Length[Select[#,Function[x,MatchQ[x,_&&blueeye&&slender]]]]==3&];

Now add hair coloring.

hairs=Flatten[Outer[List,Sequence@@Table[{brunette,Not[brunette]},{4}]],3];
galsAndEyesAndWidthAndHair=Flatten[Outer[Thread[And[##]]&,galsAndEyesAndWidth,hairs,1],1];

Again we prune.

galsAndEyesAndWidthAndHair=Cases[galsAndEyesAndWidthAndHair,{___,And[BorC_,_,_,a_],___,And[CorB_,_,_,a_],___}/;MatchQ[BorC,Betty|Carol]&&MatchQ[CorB,Betty|Carol]];
galsAndEyesAndWidthAndHair=Select[galsAndEyesAndWidthAndHair,Length[Select[#,Function[x,MatchQ[x,_&&blueeye&&_&&brunette]]]]==1&];

Johnny's ideal woman...

Finally add height.

heights=Flatten[Outer[List,Sequence@@Table[{tall,Not[tall]},{4}]],3];
galsAndEyesAndWidthAndHairAndHeight=Flatten[Outer[Thread[And[##]]&,galsAndEyesAndWidthAndHair,heights,1],1];
galsAndEyesAndWidthAndHairAndHeight=Cases[galsAndEyesAndWidthAndHairAndHeight,{___,And[AorD_,_,_,a_],___,And[DorA_,_,_,a_],___}/;MatchQ[AorD,Adele|Doris]&&MatchQ[DorA,Adele|Doris]];

This time we prune with three constraints. When finished, we will have only one ideal woman (who, as it happens, repeats herself. Possibly Johnny was not aware of this trait.)

galsAndEyesAndWidthAndHairAndHeight=Flatten[Outer[Thread[And[##]]&,galsAndEyesAndWidthAndHair,heights,1],1];
galsAndEyesAndWidthAndHairAndHeight=Cases[galsAndEyesAndWidthAndHairAndHeight,{___,And[AorD_,_,_,a_],___,And[DorA_,_,_,a_],___}/;MatchQ[AorD,Adele|Doris]&&MatchQ[DorA,Adele|Doris]];
galsAndEyesAndWidthAndHairAndHeight=Select[galsAndEyesAndWidthAndHairAndHeight,Length[Select[#,Function[x,MatchQ[x,_&&blueeye&&_&&_&&tall]]]]==2&];
galsAndEyesAndWidthAndHairAndHeight=Select[galsAndEyesAndWidthAndHairAndHeight,Length[Select[#,Function[x,MatchQ[x,__&&brunette&&tall]]]]==2&]

We see that only Adele has the requisite features.

() The envelope, please (a curvy boundary)

Or: How Gröbner bases saved me a lot of work, after I already did it

This is a problem I first encountered around eight years ago, and have spoken about on a few occasions [boundary curve]. It originated in a doctoral dissertation in number theory. The investigation focused on a trigonometric map between planar regions, and one question concerned finding the envelope curve(s) bounding the range.

trigpoly[n_,a_,b_]=Cos[n a]+Cos[n b]+Cos[n (a-b)];

I will work with a variant that is computationally more challenging than the original one.

x[a_,b_]=trigpoly[1,a,b];
y[a_,b_]=trigpoly[5,a,b];

Curvy boundary...

We plot the images of around 16000 random points from the square [0,2π]×[0,2π].

pointlist[m_]:=Table[2*Pi*{Random[],Random[]},{2^m}]
points2k=pointlist[14];
data[{point__}]:={x[point],y[point]};
datalist2k=Map[data,points2k];
ListPlot[datalist2k,PlotRange->All]

Curvy boundary...

There are several ways to go about finding the boundary curve. The one I used was to cast as an optimization problem using Lagrange multipliers (one might alternatively look for a vanishing Jacobian). The idea is for fixed x coordinate, find corresponding y values that are minimal/maximal. As we work with trig functions we will make them polynomial variables and augment with the usual identity (sin^2)+(cos^2)-1. There are alternatives but this is what I found to work well.

parameters={a,b};
mainvars={x,y};
xpoly=x-x[a,b];
ypoly=y-y[a,b];
xypolys={xpoly,ypoly};
grad[expr_,a_,b_]:={(∂_a)expr,(∂_b)expr}
trigsubs={Cos[a]->(Subscript[c, a]),Sin[a]->(Subscript[s, a]),Cos[b]->(Subscript[c, b]),Sin[b]->(Subscript[s, b])};
gradientpolys=grad[xpoly,a,b]-λ grad[ypoly,a,b];
trigidentities={(Subscript[c, a]^2)+(Subscript[s, a]^2)-1,(Subscript[c, b]^2)+(Subscript[s, b]^2)-1};
elimvars={λ,(Subscript[c, a]),(Subscript[c, b]),(Subscript[s, a]),(Subscript[s, b])};
polys=TrigExpand[Join[trigidentities,xypolys,gradientpolys]]/.trigsubs

Curvy boundary...

We have our system of polynomials. We now compute an elimination ideal, getting rid of the trigonometric variables and the Lagrange multiplier.

Originally I could not use purely symbolic GroebnerBasis to find the envelope curves (it hung). The computation involved the Mathematica functions NMinimize, FindRoot, NDSolve, NullSpace, RowReduce, and Rationalize. With a dash of LatticeReduce for extra measure. Altogether a fascinating foray into symbolic-numeric computation.

But it turns out Mathematica's GroebnerBasis can (now) handle this, if asked just right. This involves using the Gröbner walk method. It also involved repairing some trouble in that code, which is why it was not handling the problem three years ago when I first tried it.

Timing[implicitpoly=Factor[First[GroebnerBasis[polys,mainvars,elimvars,Sort->True,MonomialOrder->EliminationOrder,Method->"GroebnerWalk"]]]]

Curvy boundary...

Let's have a look at it.

Needs["Graphics`"]

ImplicitPlot[Evaluate[implicitpoly==0],{x,-1.5,3.2},{y,-2,3.2},PlotPoints->200,AspectRatio->1,PlotStyle->{{Thickness[.01],GrayLevel[.5],Dashing[{.15,.04}]},{Automatic}}]

() Independent vertices in a graph

This was discussed not long ago on the Usenet news group sci.math.symbolic [independent vertices]. We wish to find a maximally independent set of vertices in a graph. This is an NP-complete problem from combinatorial optimization. Undaunted, we proceed as below. First we need code to create a graph with a set of independent vertices of some given size (with low probability, there could be a larger set). The code is adapted from some that Maxim Rytin showed in the news group thread.

shuffle[n_Integer,m_Integer]:=Module[{res=Range[n]},Do[rand=Random[Integer,{j,n}];
res[[{j,rand}]]=res[[{rand,j}]];,{j,m}];
Take[res,m]]

Now we create a graph with 500 vertices and an independent set of 100 vertices. We show the set as it is what we hopre to reproduce.

{n,m}={500,100};
density=1/4;
SeedRandom[1111];
indepset=Sort[shuffle[n,m]]
mat=Array[Boole[#1<#2&&Random[]<density&&Complement[{##},indepset]=!={}]&,{n,n}];

Independent vertices...

The variable mat is the familiar matrix of edges, where each row and column represents a vertex.

We will recover this independent set using quadratic programming and an appropriate objective function. I show how to do this with FindMinimum. One uses the (new) QuadraticProgramming method as it seems to do well with a nonconvex problem having linear constraints. Some of the ideas behind this approach, using quadratic programming, are discussed in a technical report [Karmarker Resende Ramakrishnan report].

The key is to come up with an appropriate objective function. Well, before that, we need an appropriate formulation of the problem. We will use a variable for each vertex and constrain it to take on values between zero and one. A value of zero means the variable is not in the independent set, while one means it is, so we will further encourage them to take on one or the other of the two boundary values by adept use of the objective function.

This objective will have three components. The first is the product of matrix values above the diagonal times the corresponding vertex variables. A typical monomial looks like mat[[i,j]]*x[i]*x[j] and a value of one means we regard the two vertices as being in the independent set, but really they cannot be because there is a vertex between them. We could, I suppose, have quadratic constraints that force all these to be zero or at least "small", but I think it is better to keep the constraints simple and put our effort into the objective function.

Independent vertices...

The second component to this function will be the negated total of the variables. Without a component such as this, they'd all be zero, which would give us no vertices in our independent set (bad). So we want the objective function to try to push some variables to be large rather than small. We will weight this term by the reciprocal of the expected size of the set (this is cheating, because we know the size, but in practice one could take reasonable guesses as to what might be good weights). We use this smallish weighting because otherwise the objective function might be encouraged to give a result with too many elements in the "independent" set, in particular containing some that gave contributions near one from the monomials in the first component discussed above.

The third component to this function is similar to the second one. It is simply the sum of squares of all variables, similarly weighted.

If we have values of one for all vertex variables that are in the independent set, and values of zero for the rest, then we get values of 0, -1, and -1 respectively for the three objective function components, for a total objective function value of -2 (this, as noted above, assumes we got the weight "correct"). If other vertex variables take on positive values then the first component increases by more than the others decrease. Offhand I do not know if it is possible to decrease some of the one values, and increase some of the zero values, in such a way that the total objective function decreases. Maybe. But this is only a heuristic method, albeit one that seems to work well.

Independent vertices...

independentVerticesQP[mat_,n_,m_]:=Module[
{x,vars,poly1,poly2,poly3,poly,sol},
vars=Array[x,n];
poly1=Sum[mat[[i,j]]*x[i]*x[j],{i,n-1},{j,i+1,n}];
poly2=Total[vars];
poly3=vars.vars;
poly=poly1-poly2/m-poly3/m;
sol=FindMinimum[{poly,Flatten[{(0≤#≤1&)/@vars}]},vars,Method->QuadraticProgramming];
{sol[[1]],Flatten[Position[Round[Chop[vars/.sol[[2]]]],_?(#=!=0&),Heads->False]]}]

Timing[iset=independentVerticesQP[mat,n,m]]

Lo and behold, we have recovered our independent set.

Independent vertices...

We will illustrate the independent set on a smaller example (so the picture is not overcrowded with vertices and edges).

{n,m}={100,20};
density=1/8;
SeedRandom[1111];
indepset=Sort[shuffle[n,m]]
mat=Array[Boole[#1<#2&&Random[]<density&&Complement[{##},indepset]=!={}]&,{n,n}];

Timing[indies=independentVerticesQP[mat,n,m]]

Well, we found a set of independent vertices a bit larger than the one we started with. This is no doubt due to the low density of the graph. We will show this with blue vertices in the set and red ones outside.We'll have  light green edges connecting vertices not in the set, red edges connecting vertices in the set to vertices outside, and blue edges connecting pairs in the set (the salient feature being we do not want to see any blue edges).

iset=indies[[2]];
GraphPlot[mat,Method->"CircularEmbedding",VertexRenderingFunction->Function[{c,v},{If[MemberQ[iset,v],Blue,Red],PointSize[Medium],Point[c]}],VertexLabeling->Automatic,EdgeRenderingFunction->Function[{r,v,l},ii=Intersection[iset,v];{If[ii=={},Green,If[Length[ii]==2,Blue,Red]],Line[r]}]]

() Depicting cylinders containing five given points

Over the past several years I have worked on various aspects of the following problem: given five points in space, find all right circular cylinders containing them. A bit of thought indicates that, generically, one needs four values to define the axial line, and a fifth necessary value is the radius, so having an equation to be satisfied by each point means five is the number we expect will give us finitely (but not zero) actual cylinders. Well, we might need to go to complex space, but... Anyway, I do not propose to go into a lengthy discussion on this topic (I have done so in past, e.g.[cylinders through five points]).

One configuration that is of interest is that where the five points are the vertices of two regular tetrahedra glued together along a common face. We orient so that this face lies in the horizontal plane, and thus the last two points are aligned vertically In this case there happen to be six real cylinders containing all five points, and it can be shown that this is the largest (other than the nongeneric cases where there are infinitely many, e.g. the points are all collinear).

Some general theory indicates that if we raise the top point vertically, we continue to have six cylinders until we hit a particular value at which...poof...they all disappear (which is to say, solutions to the radial and axial parameters become complex valued). I always found this interesting and wanted to visualize it. With some (okay, considerable) help from Jeff Bryant, we now have a Manipulate example that will show this.

Depicting cylinders...

First we require a small amount of code to solve for the cylinder parameters. While I do not propose to explain this in any detail, the idea is to represent the axial line as offset+t vec where vec={a,c,1} is a (generic) direction vector for the cylinder axis, and offset={b,d,0} is a point on the axis. To obtain an equation we project our point onto the axis, find the distance between the point and this projected value, and insist that equal the (unknown) radius r. Again, five points yield five equations relating above cylinder parameters (a,b,c,d,r). We use NSolve once again as our computational engine.

perp[vec1_,vec_,offset_]:=vec1-offset-Projection[vec1-offset,vec,Dot]

solveCylinders[pts_List,vec_,offset_,prec_:Automatic]:=Module[{exprs,k,perps},perps=Table[perp[pts[[k]],vec,offset],{k,5}];
exprs=(Numerator[Together[#1.#1-rsqr]]&)/@perps;
NSolve[exprs,{a,b,c,d,rsqr},WorkingPrecision->prec]]

Next we put this into a manipulation.

Depicting cylinders...

Manipulate[dpoints={{1,0,0},{-1/2,Sqrt[3]/2,0},{-1/2,-Sqrt[3]/2,0},{0,0,k*Sqrt[2]},{0,0,-Sqrt[2]}};
vec={a,c,1};
offset={b,d,0};
solns=solveCylinders[dpoints,vec,offset];
top[j_]:=(2*vec+offset)/.solns[[j]];
bot[j_]:=(-2*vec+offset)/.solns[[j]];
{p1,p2,p3,p4,p5}=dpoints;
Graphics3D[{ControlActive[Opacity[1],Opacity[.6]],Red,Sphere[#,0.04]&/@{p1,p2,p3},Green,Sphere[p4,0.04],Blue,Sphere[p5,0.04],White,Polygon[{{p1,p2,p4},{p2,p3,p4},{p3,p1,p4}}],Polygon[{{p1,p2,p5},{p2,p3,p5},{p3,p1,p5}}],Yellow,Cylinder[{top[1],bot[1]},0.70],Blue,Cylinder[{top[2],bot[2]},0.70],
Orange,Cylinder[{top[3],bot[3]},0.70],
Green,Cylinder[{top[4],bot[4]},0.70],
Red,Cylinder[{top[5],bot[5]},0.70],Purple,Cylinder[{top[6],bot[6]},0.70]
},Axes->True,AxesLabel->{"x","y",None},PlotRange->{{-3,3},{-3,3},{-3,3}}],{k,.26,5/4.+.1},TrackedSymbols:>{k}]

() A logic puzzle with 0-1 inequality constraints

This next example illustrates a type of popular logic puzzle often found in supermarket and airport publications. This one is called "Dinner on the run" [Jean Hannagan puzzle]; I cribbed it from my son's puzzle magazine. A description of the problem is as follows. Five men order five sandwiches from Max's Deli, such that each has a distinct fillings, toppings, spreads, and breads. Further clues are as below, and we are to deduce the components of each person's sandwich. We follow customary rules of interpretation of such puzzles, e.g. "Neither the sandwich with cheese nor the one with mayo was on white bread" means, among other things, that the cheese and mayo sandwiches are themselves distinct.

(1) One filling is salami. One topping is tomatoes. One spread is butter. One bread is pumpernickel. Frank placed the order.
(2) Regardless of order, Max refuses to put mustard or ketchup on tuna salad, or onions on turkey.
(3) Tim does not like onions.
(4) The roast beef was not on rye, and did not contain pickles.
(5) The tuna salad was not on seven-grain bread.
(6) Neither tuna salad nor sandwhich with lettuce was spread with mayo.
(7) The turkey, which was on whole wheat, did not have mustard.
(8) The sandwich with lettuce did not have mustard.
(9) Neither Nick nor the person requesting ketchup and onions used white bread.
(10) The ham sandwich had neither pickes nor onions.
(11) The sandwich with pickles did not have mayo.
(12) Jim did not order relish.
(13) The sandwich with lettuce did not have relish.
(14) The sandwich with cheese did not have mustard. Moreover it did not go to Tim.
(15) Neither the sandwich with cheese nor the one with mayo was on white bread.
(16) Elmer does not eat rye bread.
(17) Tim did not order his with whole wheat bread.
(18) The sandwich with whole wheat did not have mayo.
(19) Nick did not order tuna salad.
(20) Jim did not order seven-grain bread.
(21) Neither Jim nor Nick used lettuce.

Logic puzzle...

We first use the information just to figure out the five items in each of five categories. The men are Frank, Jim, Nick, Tim, and Elmer. The fillings are salami, tuna salad, turkey, ham, and roast beef. The toppings are tomatoes, lettuce, onions, cheese, and pickles. The spreads are butter, ketchup, mustard, mayo, and relish. The breads are pumpernickel, seven-grain, white, rye, and wheat. Note that this preprocessing step is in itself a bit of a challenge. (Is relish a topping or spread? Is cheese a filling or topping?)

As there are a scant (((5!))^4) possibilities (around 200 million, that is) we could, with some work, do a brute force search. We will use a constraint satisfaction approach via linear programming as we believe it is instructive to demonstrate handling of inequations in this setting.

Conceptually, we arbitrarily number the men Frank=1,...,Elmer=5. We create variables for each of the categories: filling, top, spread, and bread. For each such variable we create subvariables, one per sandwich number. These must be zero or one, and they sum to one (that is, exactly one has value one). For example, we will have the equation

(∑_(j=1)^5)mayo[j]=1

Moreover we have transposes of the above type of equation. That is, we know the sum of all subvariables for a particular sandwich and category must be also be one. For example, we have

mayo[1]+ketchup[1]+mustard[1]+relish[1]+butter[1]=1

We furthermore have equations relating the basic variables to their corresponding subvariables, such as

mayo=mayo[1]+2 mayo[2]+3 mayo[3]+4 mayo[4]+5 mayo[5]

Logic puzzle...

For efficiency we only work with the knapsack variables, using tactics equivalent to these above equations at the end to put the solution into a reasonable form.

Finally we have all the inequations and the handful of equations implied by the 21 rules above. These we handle as follows. From rule 7, say, we know that turkey=wheat. This means we equate all the corresponding subvariables. From rule 11 we know that pickles≠mayo. Some reflection shows that this equivalent to the set of subvariable inequalities mayo[j]+pickles[j]≤1 for 1≤j≤5. When we equate a pair of items, e.g. turkey and wheat, we realize that by equating each pair of their five subvariables we can remove one. Similarly observe that, for example, Jim did not take lettuce, so we have an equation for the variable representing lettuce used by Jim (it is zero). Again we can remove it from further consideration after we use it to simplify whatever constraints contained that variable. Use of such equations as means to remove variables and constraints tends to make the computation much faster. In this case it gave a speed improvement by a factor of around 10.

As all constraints are linear, we will solve this by treating it as a linear programming problem. Curiously we are not actually optimizing anything; it is an example of a "constraint satisfaction" problem. This gives us some freedom to make up objective functions, and in general that can be used to improve efficiency, but we will forego that in the interest of simplicity.

Logic puzzle...

The idea of the branching method is to solve what are termed relaxations of the problem wherein integrality is not enforced but inequality constraints are in use. For any solution one looks for noninteger parts. If all are integer the solution is fine. Otherwise one takes, say, the first coordinate that is not an integer and spawns a pair of new problems, one constraining the corresponding variable to be less or equal to the floor of the value in the solution, and the other constraining it to be greater or equal to the ceiling of that value. This is the "branching". We do not utilize the "bounding" part in this example, but in cases where there is an actual objective function of interest, one uses it to remove subproblems for which we know there is no hope of improving on the current "best" solution (once a viable solution has been found). Note that this approach is also useful for what are referred to as mixed problems, wherein some but not all variables may be constrained to take on integer values. One simply branches only on those so constrained.

The following code will quit after finding one valid solution. We are trusting that the author used adequate conditions to ensure that a solution exists and is unique (which is in fact the case). It is not hard to modify the code to allow for other possibilities. Note that in contrast to the method I presented yesterday [ILP talk notebook] for handling Frobenius instance problems, we now have many integer variables but constrained to take on values in a small set. In this case it is neither necessary nor even desirable to preprocess with lattice reduction, or to base branching decisions thereon.

Logic puzzle...

findSandwiches[]:=Module[{consumers,fills,tops,spreads,breads,frank,jim,nick,tim,elmer,salami,tuna,turkey,ham,roastbeef,tomatoes,lettuce,onions,cheese,pickles,butter,ketchup,mustard,mayo,relish,pumpernickel,sevengrain,white,rye,wheat,fillvars,topvars,spreadvars,breadvars,varlists,all01vars,pv,prodvars,constraints1,constraints2,constraints3,constraints4,constraints,vars,program,stack,soln,soln01,badpos,counter=0,eps=10^(-6),names,suffix,partsoln,psvars},consumers={frank,jim,nick,tim,elmer};
fills={salami,tuna,turkey,ham,roastbeef};
tops={tomatoes,lettuce,onions,cheese,pickles};
spreads={butter,ketchup,mustard,mayo,relish};
breads={pumpernickel,sevengrain,white,rye,wheat};
fillvars=Outer[#1[#2]&,fills,consumers];
topvars=Outer[#1[#2]&,tops,consumers];
spreadvars=Outer[#1[#2]&,spreads,consumers];
breadvars=Outer[#1[#2]&,breads,consumers];
varlists={fillvars,topvars,spreadvars,breadvars};
all01vars=Flatten[varlists];
constraints1=Map[0≤#≤1&,all01vars];
constraints2=Map[Apply[Plus,#]==1&,{fillvars,topvars,spreadvars,breadvars},{2}];
constraints3=Map[Apply[Plus,#]==1&,Map[Transpose,{fillvars,topvars,spreadvars,breadvars}],{2}];
ineq[v1_,v2_,v3__]:={ineq[v1,v2],ineq[v1,v3],ineq[v2,v3]};
ineq[v1_,v2_]:=Map[v1[#]+v2[#]≤1&,consumers];
eq[v1_,v2_]:=Map[v1[#]==v2[#]&,consumers];
constraints4={ineq[mustard,tuna],ineq[ketchup,tuna],ineq[onions,turkey],ineq[roastbeef,rye],ineq[roastbeef,pickles],ineq[tuna,sevengrain],ineq[tuna,lettuce,mayo],ineq[turkey,mustard],ineq[mustard,lettuce],ineq[ketchup,white],ineq[ham,pickles],ineq[ham,onions],ineq[pickles,mayo],ineq[lettuce,relish],ineq[cheese,mustard],ineq[cheese,mayo,white],ineq[mayo,wheat]};
constraints=Flatten[Join[constraints1,constraints2,constraints3,constraints4]];
partsoln=Flatten[{eq[turkey,wheat],eq[ketchup,onions],onions[tim]==0,rye[elmer]==0,wheat[tim]==0,white[nick]==0,onions[nick]==0,relish[jim]==0,cheese[tim]==0,tuna[nick]==0,sevengrain[jim]==0,lettuce[jim]==0,lettuce[nick]==0}];
psvars=Map[First,partsoln];
partsoln=partsoln/.Equal->Rule;
constraints=constraints//.partsoln/.{True->Sequence[],((aa_/;Head[aa]=!=Plus)≤1)->Sequence[]};
all01vars=Complement[all01vars,psvars];
program=constraints;
stack={program,{}};
While[stack=!={},counter++;program=stack[[1]];stack=stack[[2]];
Internal`DeactivateMessages[soln=NMinimize[{1,program},all01vars]];
If[!FreeQ[soln,Indeterminate],Continue[]];
soln=Chop[soln[[2]],eps];
soln01=all01vars/.soln;
badpos=Position[soln01,(aa_/;Chop[aa-Round[aa],eps]=!=0),{1},1,Heads->False];
If[badpos=={},partsoln=partsoln/.soln;
soln=Join[soln,partsoln];
soln=Split[Sort[(soln/.aa_?NumberQ:>Round[aa])/.{HoldPattern[_->0]:>Sequence[],HoldPattern[hh_[bb_]->1]:>bb==hh,HoldPattern[Rule[_,bb_/;Not[NumberQ[bb]]]]:>Sequence[]}],#1[[1]]===#2[[1]]&];
names=Map[#[[1,1]]&,soln];
soln=MapThread[Prepend,{soln/.aa_==bb_:>bb,names}];
soln=Map[ToString,soln,{2}];
suffix=soln[[1,1]];
suffix=StringDrop[suffix,StringPosition[suffix,"$"][[1,1]]-1];
soln=Map[StringReplace[#,suffix->""]&,soln,{2}];
Break[],badpos=badpos[[1,1]];
stack={Append[program,all01vars[[badpos]]==0],stack};
stack={Append[program,all01vars[[badpos]]==1],stack};];];
{counter,soln}]

Logic puzzle...

After all that work, solving is a breeze. Who'da thunk it?

SetOptions[LinearProgramming,Method->InteriorPoint];

Timing[orders=findSandwiches[]]

Similar methods can be used for solving, say, sudoku puzzles.

() Summary

Well, here we are. As promised, we touched upon Gröbner bases, logic puzzles, branch-and-bound, and more. The methods presented apply to a vast array of problems in mathematical optimization and elsewhere.

Should I survive the expected tomato attack, I hope to return to problems from MathGroup and elsewhere in a future talk, but perhaps with more of a number theory and calculus focus.

() Acknowledgements

I thank Jeff Bryant, Brett Champion, Yifan Hu, Michael Trott, and Robby Villegas for assisting in various parts of the preparation of these slides. Mostly small matters, such as gathering the shards of my fractured sanity.

() Selected references

[Algebraic replacement]
http://forums.wolfram.com/mathgroup/archive/2006/Aug/msg00283.html

[Parametrized nested radical]
http://forums.wolfram.com/mathgroup/archive/2005/May/msg00502.html

[Parametrized rational surface]
A. Chtcherba and D. Kapur. Conditions for determinantal formula for resultant of a polynomial system. Proceedings of the 2006 International Symposium on Symbolic and Algebraic Computation (ISSAC 2006), 55-62. ACM Press 2006.
See also: J. Zheng, T. Sederberg, E-W Chionh, and D Cox. Implicitizing rational surfaces with base points using the method of moving surfaces. In: Topics in Algebraic Geometry and Geometric Modeling. Contemporary Mathematics volume 334, 151-168. AMS 2003.

[Johnny's ideal woman]
http://forums.wolfram.com/mathgroup/archive/2003/Sep/msg00219.html
http://forums.wolfram.com/mathgroup/archive/2003/Sep/msg00258.html
The original statement of the problem is from Georgia College & State University BITS & BYTES 3:3, February 1998.
http://www.gcsu.edu/acad_affairs/coll_artsci/mathcomp_sci/bits/V3N3Feb98.html
This work has also appeared in The Mathematica Journal 9(2):289-291, "Tricks of the Trade" column, edited by P. Abbott.

[Who own's the zebra?]
M. Trott. Solving puzzles with Mathematica. In "Trott's Corner" column, The Mathematica Journal 7(3):291-307 (puzzle 3). 1999.

[Trigonometric boundary curve]
D. Lichtblau. Computing curves bounding trigonometric planar maps: symbolic and hybrid methods. Proceedings of the 2004 Workshop on Automated Deduction in Geometry (ADG 2004). Lecture Notes in Computer Science 3763, 70-91. 2006.
http://www.springerlink.com/content/2267x4v48wpg1680/?p=0e3d41cd8764421c9518b74e42de8892&pi=12

[Independent vertices]
http://groups.google.com/group/sci.math.symbolic/msg/c06cbe2ba0a3ace7?dmode=source

[Karmarker, Resende, Ramakrishnan report]
An interior point approach to the maximum independent set problem in dense random graphs. Mathematical Sciences Research Center, AT&T Bell Laboratories, Murray Hill, NJ. 1989.

[Cylinders through five points]
http://webs.uvigo.es/adg2006/talks.html
Daniel Lichtblau. Cylinders Through Five Points: Complex and Real Enumerative Geometry.

[Hannagan's "Dinner on the Run"]
J. Hannagan. Dinner on the run. Dell Variety Puzzles and Word Games 135:62. May 2004.

[Integer linear programming intro]
D. Lichtblau. An Introduction to Integer Linear Programming in Mathematica. Notebook for presentation at Wolfram Technology Conference 2006.
http://library.wolfram.com/infocenter/Conferences/6449