# Formatting an equation: X^Y to pow(X,Y)

I am trying to solve for the intersection of two circles with arbitrary centers and arbitrary radii.

(x-X1)^2+(y-Y1)^2==R1^2
(x-X2)^2+(y-Y2)^2==R2^2
Solve[% && %%,{x,y}]


It takes a while on my raspberry pi, so don't do it if you dont have to.

The output is around 47000 characters. I am going to format the answer and place it into a computer program written in OpenSCAD (think C code). OpenSCAD does not recognize the carret "^" operator, so I need every instance of X^Y in the output converted to pow(X,Y).

How do I do that?

I am working on mathematica via the command line, from a Windows 7 computer over my LAN. Assume I dont have access to the Mathematica IDE. Once I have the output I will then export the output to a text file, which I will then open in Windows, do a bit of editing, and then copy the formatted equations into OpenSCAD for processing.

Thanks.

• Do you mean for the second circle to be centered at {X1, Y2}, or is that a typo for {X2, Y2}? – 2012rcampion Apr 23 '15 at 15:27
• I guess you should be interested in this question Controlling measure zero sets of solutions with Manipulate. A case study – Artes Apr 23 '15 at 15:36
• If you're not restricted in introducing a few extra convenient substitutions, such as d=Sqrt[(x2-x1)^2+(y2-y1)^2] to somewhat shorten the calculations, the answer becomes easily expressible in just a few lines. This page can give some insights. Failing that, use ReplaceRepeated with the rule Power[a_,b_]->pow[a,b]. Mathematica will complain if you use the other (non-square) brackets, so some conversion to strings may be necessary. – LLlAMnYP Apr 23 '15 at 15:50
• may by you can use string:StringReplace[ToString[FullForm[x^y]], {"[" -> "(", "]" -> ")"}]? – Basheer Algohi Apr 23 '15 at 17:01
• There is always good ol' CForm. – wxffles Apr 23 '15 at 19:56

## 1 Answer

First I'll store the equations to a variable:

eqns = {(x - x1)^2 + (y - y1)^2 == r1^2, (x - x2)^2 + (y - y2)^2 == r2^2}


Note that I've switched to lowercase variable names; this is generally good practice to avoid conflicts with Mathematica's builtin identifiers, which are always uppercase.

Without loss of generality we can assume that the first circle has unit radius and is centered on the origin, and that the second circle's center lies on the x-axis. If this is not the case, we can always transform into a coordinate system where it is so.

This corresponds to the following substitutions:

soln = Solve[eqns /. {x1 | y1 | y2 -> 0, r1 -> 1, x2 -> d, r2 -> Sqrt[1 + d rr]}, {x, y}];
soln = FullSimplify @ soln
StringLength[ToString[soln, InputForm]]
(* 127 *)


The two solutions have the same x, and y which differ by only a sign. Therefore we only need to calculate one of the solutions. Let's take a closer look:

pt1 = {x, y} /. First[soln]

(* {(d - rr)/2, -(1/2) Sqrt[-(-2 + d - rr) (2 + d - rr)]} *)


It seems that y is actually a function of x. We can rewrite the expression like this:

x0 = (d - rr)/2;
pt2 = {x0, -Sqrt[(1 - x0) (1 + x0)]};
Reduce[pt2 == pt1] (* === True *)


Now all we need to do is calculate x0, and then translate back to the original coordinate system. For the transformation, we'll first rotate and scale to the correct size, then translate to the correct location. We can construct the rotation matrix from {d, 0} to offset coordinates {dx, dy} like so:

Simplify[ComplexExpand[RotationMatrix[{{d, 0}, {dx, dy}}]] /. {dx^2 + dy^2 -> d^2}, d > 0]
(* {{dx/d, -(dy/d)}, {dy/d, dx/d}} *)


We can then construct the full solution like this:

dx = (x2 - x1)/r1;
dy = (y2 - y1)/r1;
d = Sqrt[dx^2 + dy^2];
rr = ((r2/r1)^2 - 1) / d;

x0 = (d - rr)/2;
y0 = Sqrt[(1 - x0) (1 + x0)];

scale = r1 / d;

x1soln = (dx*x0 + dy*y0) * scale + x1;
y1soln = (dy*x0 - dx*y0) * scale + y1;

x2soln = (dx*x0 - dy*y0) * scale + x1;
y2soln = (dy*x0 + dx*y0) * scale + y1;


We can test the validity of the solutions by plugging them back into the original equations:

Simplify[eqns /. {{x -> x1soln, y -> y1soln}, {x -> x2soln, y -> y2soln}}]
(* {{True, True}, {True, True}} *)


There are some optimizations you can do to make the computation faster. I won't go over this in too much detail, since it's all done by hand, but here's my best result:

dx = (x2 - x1);
dy = (y2 - y1);
d2 = dx^2 + dy^2;
rr = (r2 - r1) * (r2 + r1) / d2;

x0 = (1 - rr)/2;
y0 = Sqrt[r1*r1/d2 - x0*x0];

x1soln = x1 + dx*x0 + dy*y0;
x2soln = x1 + dx*x0 - dy*y0;

y1soln = y1 + dy*x0 - dx*y0;
y2soln = y1 + dy*x0 + dx*y0;

• I saved some expensive divide operations by propagating the /r1 through the whole expression until it cancelled out with a *r1.
• I saved a multiplication by computing rr with (r2-r1)*(r2+r1) instead of r2*r2 - r1*r1.
• I traded the Sqrt in d for an extra multiplication in y0.
• Some more cancellations occurred and scale turned into 1, saving four multiplications.
• By rearranging the final sums, hopefully your compiler will save the products from one expression to the next, saving four multiplications.

Overall we save:

• One square root (expensive)
• Three divisions
• Two multiplications (up to six with compiler optimization)
• (up to two additions with compiler optimization)

## Update

You can try to get Mathematica to handle this automatically using an experimental/undocumented function. We can start exactly as you did:

soln = {x, y} /. Solve[{(x - x1)^2 + (y - y1)^2 == r1^2,
(x - x2)^2 + (y - y2)^2 == r2^2}, {x, y}]

soln = FullSimplify[soln];


Note that the second step takes around 350x longer (12.5 s vs 30 ms) than the first on my laptop (i7, 1.9 GHz), so it may take a long time on the RPi (ARM11, 0.7 GHz).

We can then use ExperimentalOptimizeExpression to try and apply 'compiler optimization':

compiled = ExperimentalOptimizeExpression[soln];
expression = Map[Hold, compiled, {2}][[1, 2]];
symbols = Union@Cases[expression, s_Symbol /; (Context[s] == "Compile"), Infinity];
exprString =
ToString[expression /.
Thread[symbols -> (Symbol["temp" <> IntegerString[#]] &) /@
Range[Length[symbols]]], CForm]
StringReplace[exprString, {RegularExpression["^Hold$\$$|$\$$$"] -> "", RegularExpression[",(?![^(]\\))"] -> ";\n", RegularExpression["(\\d)\\.(\\D)"] -> "$1\$2"}]


This outputs some C-like code:

temp1 = -x2;
temp2 = x1 + temp1;
temp3 = Power(temp2,2);
temp4 = -y2;
temp5 = y1 + temp4;
temp6 = Power(temp5,2);
temp7 = temp3 + temp6;
temp20 = -temp3;
temp21 = -temp6;
temp8 = 1/temp7;
temp11 = Power(r1,2);
temp9 = Power(r2,2);
temp30 = Power(r1,4);
temp31 = -temp9;
temp32 = temp31 + temp3 + temp6;
temp33 = Power(temp32,2);
temp34 = temp9 + temp3 + temp6;
temp35 = -2*temp11*temp34;
temp36 = temp30 + temp33 + temp35;
temp37 = -(temp36*temp6);
temp38 = Sqrt(temp37);
// ...


It goes on for about 45 lines. You can add whatever extra string replacement rules you need to get the correct formatting (e.g. "Power" -> "Math.pow").

Note that there are a lot of silly things going on here, like temp31 = -temp9. Why not just make temp9 = -Power(r2,-2), or do temp32 = (-temp32) + ...? I think that OptimizeExpression optimizes for speed above all, so it's not guaranteed to be compact, which is maybe not what you want.

• heres a case where I suspect google-ing 'circle intersection algorithm' would have saved a lot of time.. – george2079 Apr 24 '15 at 15:17
• @george2079 Probably, but I like to pontificate =) – 2012rcampion Apr 24 '15 at 15:23
• Thank you @2012rcampion. I'm trying to implement the solution now. You are right, that sqrt substitution I think was the key. The Solve[] took 0.18 seconds, even on my Pi. Thanks for teaching me about optimization (in general), I think I'm going to need to understand a bunch of this if I'm going to continue to work on a Pi. – DrXenocide Apr 25 '15 at 9:30
• @2012rcampion: in equation blocks 6 and 7 is dy=x2-x1? or dy=y2-y1? – DrXenocide Apr 25 '15 at 9:33
• @DrX Yes, my bad. Let me double check the rest of my code real quick... Ok, should be fixed now. – 2012rcampion Apr 25 '15 at 14:07