From: lewis@bway.net (Robert H. Lewis) Subject: Apollonius Sphere and Circle Again Date: Mon, 26 Mar 2001 12:18:35 -0500 Newsgroups: sci.math.symbolic Summary: Appolonius (kissing-sphere) problem as a CAS/Grobner basis challenge Greetings, A few weeks ago I posted here a query about the "Appolonius Sphere problem" (also spelled "Apollonius"). Briefly, the Appolonius Circle Problem is a well known problem of geometry. The problem is to find the circle(s) tangent to three given circles in the plane. I could not find any solution on the web to the obvious generalization: find a sphere tangent to four given spheres. The problem was brought to my attention by Stephen Bridgett of the U.K. who says it has potential medical applications. He wrote to me: > I'm trying to find the exact equation for the > coordinates of the centre (sx,sy,sz) and radius (sr) > of the sphere (in some cases there will be two such > spheres), which tangentially touches four other > spheres with centres at (ax,ay,az), (bx,by,bz), > (cx,cy,cz), (dx,dy,dz) and radii of ar,br,cr,dr > respectively. This means solving the following > equations to find sx,sy,sz and sr: > > (sx-ax)^2 + (sy-ay)^2 + (sz-az)^2 = (sr+ar)^2 > (sx-bx)^2 + (sy-by)^2 + (sz-bz)^2 = (sr+br)^2 > (sx-cx)^2 + (sy-cy)^2 + (sz-cz)^2 = (sr+cr)^2 > (sx-dx)^2 + (sy-dy)^2 + (sz-dz)^2 = (sr+dr)^2 > > I've used Mathematica to solve the 2D case of > a circle touching three other circles, but so > far Mathematica has failed with the above 3D > case. > I presented this as a challenge to computer algebra systems -- and their users! Can you solve the above system of equations for sx, sy, sz, and sr? My system Fermat (web page below) solved it rather easily. It found a quadratic equation for sx in terms of the given parameters ax, ay, az, bx, ... , dr. Similarly for sy, sz, sr. Secondly and easier, can your system solve the equations after substituting in ax=ay=az=0? This is just a translation so is certainly valid mathematically. One could also plug in az=1 if necessary, but let's not, just for fun! ____________________________________ So, the "full problem" solves the equations as given above, and the "reduced" or "simplified" problem first plugs in ax=ay=az=0. I have received a number of interesting replies: (1) Allan Wilks writes that the problem is completely solved theoretically in a recent preprint by Mallows, Lagaris, and him. Since he didn't post it here I won't give the url. However, the problem solved in that paper is much more restricted -- they assume that ALL FOUR circles (five spheres, ...) are tangent to each other, and call it "the Descartes circle theorem." Here we do not require that the original circles or spheres be tangent to each other. (2) Mark D. (Danny) Rintoul has a clever idea to simplify the equations: "First, subtract the 4th equation from the first 3. This will give you a set of 3 equations that are linear in the 4 unknowns (sx,sy,sz,and sr). Now, solve for sx, sy, and sz in terms of the knowns *and* sr (assuming it is known). Then using this equation for the vector S in terms of s_r, calculate what s_r needs to be to satisfy the fourth equation." (3) Henry Cejtin writes: "A fine test problem. I have a simple Groebner basis library, written in Scheme. Nothing fancy. It does various monomial orderings. I used an elimination-ordering and to get the equation for sr (in terms of the a's, b's, c's and d's) took just over 1 hour, and it was using about 50 meg of RAM. Doing it for the special case where ax, ay and az are all = 0 took 46 seconds and used a bit under 8 meg of RAM. Solving for sx took 1.75 hours CPU time and used 60 meg RAM. In the ax = ay = az = 0 special case it took 300 seconds CPU time and used under 13 meg of RAM. All of these times are on an 400 MHz Pentium II machine." These are by far the best times I have seen for Groebner basis methods. (4) Stephen Bridgett: "I ran Mathematica with the simplified version of the equations (ie. ax=ay=az=0) with the Solve command. I left it running all night, and this morning it had a solution for 'sx'. The Timing[] function in Mathematica for computing sy returned a time of 4.6 hours for the simplified case of ax=ay=az=0. I also found a function called MaxMemoryUsed[], which returned 41,191,720 bytes. [It could not do the full problem.] (5) Robert Israel writes: In Maple 6 (Sun SPARC Solaris) The full problem, with(Groebner): eqs:= {(sx-ax)^2 + (sy-ay)^2 + (sz-az)^2 - (sr+ar)^2,    (sx-bx)^2 + (sy-By)^2 + (sz-bz)^2 - (sr+br)^2,    (sx-cx)^2 + (sy-cy)^2 + (sz-cz)^2 - (sr+cr)^2,    (sx-dx)^2 + (sy-dy)^2 + (sz-dz)^2 - (sr+dr)^2}: (I couldn't use "by" as a variable because this is a reserved word in Maple) G:= gbasis(eqs, plex(sx,sy,sz,sr)): I stopped it without reaching a solution after bytes used=647,253,504, alloc=103265824, time=2526.36 [secs] > Simplified problem ax=ay=az=0: > eqs:= subs(ax=0,ay=0,az=0,eqs):   G:= gbasis(eqs, plex(sx,sy,sz,sr)): bytes used=574,205,204, alloc=63034088, time=839.40 seconds. The Groebner basis G consists of 4 polynomials: the first is a quadratic in sr with 2457 terms, the others each have 48 terms and are of the form a sr + b sz + c, a sr + b sy + c and a sr + b sx + c respectively.  Of course, given G we could get a Groebner basis for the original problem by substituting sx - ax for sx, bx - ax for bx, etc.  After expanding, the quadratic in sr then has 18366 terms, and the other three polynomials have 144 each. The problem could be made much more tractable: in addition to [ax,ay,az]=0 you could take by=bz=0, cz=0, and ar = 1 (i.e. choose the origin as the centre of one sphere, the unit of length as the radius of this sphere, the x axis as the direction to the second sphere, and the xy plane as the plane containing the three centres).  Then Maple takes less than one second of CPU time to come up with the Groebner basis. Conclusion: always use symmetry as much as possible to eliminate unnecessary variables. (6) My own solution (Robert H. Lewis): I used my system Fermat (see URL below) with the Dixon-KSY resultant method [Kapur-Saxena-Yang paper in the 1996 ISSAC conference.] The reduced problem took 285 seconds, 45 meg of RAM to get sx, working over Z, the integers as ground ring. Working mod a prime (44449) reduces it to 136 seconds and 37 meg of RAM. The full problem for sx takes 2200 seconds (almost all of which is the computation of the Content of the original 1.3 million term equation) and about 390 meg of RAM. Working mod p=44449 reduces this to 205 seconds. This is all on a Macintosh blue and white G3, 400 mhz. Finding sr is much easier, over Z, for the reduced problem: it takes 19 seconds and about 2 meg of RAM. But the full problem takes 425 seconds over Z/44449, and I ran out of memory (approx. 400 meg) trying it over Z. Thanks to everyone who contributed! Robert H. Lewis Mathematics Department Fordham University Bronx NY 10458 www.bway.net/~lewis/ -----= Posted via Newsfeeds.Com, Uncensored Usenet News =----- http://www.newsfeeds.com - The #1 Newsgroup Service in the World! -----== Over 80,000 Newsgroups - 16 Different Servers! =-----