本日のC.A.D.

前回は S、T が ￢(0∈S) ∨ ￢(0∈T) を満たす場合でしたが，一般には

　　　∃s∃t(z=s*t ∧ s∈S ∧ t∈T)
　⇔　∃s∃t(z=s*t ∧ s∈S ∧ t∈T ∧ (s≠0 ∨ s=0))
　⇔　∃s∃t(z/s=t ∧ s∈S ∧ t∈T) ∨ ∃s∃t(z=0 ∧ s=0 ∧ s∈S ∧ t∈T)
　⇔　∃s(s∈S ∧ z/s∈T) ∨ (z=0 ∧ 0∈S ∧ T≠{})

のように選言となり，とくに，S={w;|w-c1|=r1}，T={w;|w-c2|=r2} の場合，この選言は次のように消去できます．

　　　∃s∃t(z=s*t ∧ |s-c1|=r1 ∧ |t-c2|=r2)　………　（３）
　⇔　∃s(|s-c1|=r1 ∧ |z/s-c2|=r2) ∨ (z=0 ∧ |c1|=r1 ∧ r2≧0)
　⇔　∃s(|s-c1|=r1 ∧ |z-c2*s|=r2*|s| ∧ s≠0) ∨ (z=0 ∧ |c1|=r1 ∧ r2≧0)
　⇔　∃s(|s-c1|=r1 ∧ |z-c2*s|=r2*|s| ∧ s≠0 ∧ r2≧0) ∨ ∃u(|u-c1|=r1 ∧ |z-c2*u|=r2*|u| ∧ u=0 ∧ r2≧0)
　⇔　∃s(|s-c1|=r1 ∧ |z-c2*s|=r2*|s| ∧ (s≠0 ∨ s=0) ∧ r2≧0)
　⇔　∃s(|s-c1|=r1 ∧ |z-c2*s|=r2*|s| ∧ r2≧0)　………　（４）

なお，（３）⇒（４）は z=s*t の代入で，（４）⇒（３）は s≠0，s=0 の場合にそれぞれ t=z/s，t=c2+r2 として得ることもできます．

さて，R1=r1^2，R2=r2^2 などとおき，（３），（４）を実変数化すると

f3 = Exists[{x1, y1, x2, y2}, y1*x2 + x1*y2 == y && x1*x2 - y1*y2 == x && (-a + x1)^2 + (-b + y1)^2 == R1 && (-c + x2)^2 + (-d + y2)^2 == R2];
f4 = Exists[{x1, y1}, (x1 - a)^2 + (y1 - b)^2 == R1 && (-(c*x1) + d*y1 + x)^2 + (-(c*y1) - d*x1 + y)^2 == R2*(x1^2 + y1^2) && 0 <= R2];

となり，f4 の 2 つの等式から x1^2 + y1^2 を消去した（根軸を表す）等式

(-(c*x1) + d*y1 + x)^2 + (-(c*y1) - d*x1 + y)^2 - R2*(x1^2 + y1^2) - (c^2 + d^2 - R2)*((x1 - a)^2 + (y1 - b)^2 - R1) == 0

の左辺を P+x1*Q+R*y1，つまり

PQR := {P->(-a^2)*c^2 - b^2*c^2 - a^2*d^2 - b^2*d^2 + c^2*R1 + d^2*R1 + a^2*R2 + b^2*R2 - R1*R2 + x^2 + y^2,
 Q->2*a*c^2 + 2*a*d^2 - 2*a*R2 - 2*c*x - 2*d*y,
 R->2*b*c^2 + 2*b*d^2 - 2*b*R2 + 2*d*x - 2*c*y};

とおくと，QEPCADB が

[]
(A,B,C,D,X,Y)
4
(EX)(EY)[X^2 + Y^2 = D /\ A X + B Y + C = 0].
finish
An equivalent quantifier-free formula:
D >= 0 /\ B^2 D + A^2 D - C^2 >= 0

と答えることにより，f4 は次のように QE できます．

Exists[{x1, y1}, (x1 - a)^2 + (y1 - b)^2 == R1 && P + Q*a + R*b + (x1 - a)*Q + R*(y1 - b) == 0 && 0 <= R2]
⇔
(P + Q*a + R*b)^2 - R1*(Q^2 + R^2) <=0  && 0 <= R1 && 0 <= R2
⇔ (using PQR)
((-a^2)*c^2 - b^2*c^2 - a^2*d^2 - b^2*d^2 + c^2*R1 + d^2*R1 + a^2*R2 + b^2*R2 - R1*R2 + x^2 + y^2 + b*(2*b*c^2 + 2*b*d^2 - 2*b*R2 + 2*d*x - 2*c*y) + a*(2*a*c^2 + 2*a*d^2 - 2*a*R2 - 2*c*x - 2*d*y))^2 - R1*((2*b*c^2 + 2*b*d^2 - 2*b*R2 + 2*d*x - 2*c*y)^2 + (2*a*c^2 + 2*a*d^2 - 2*a*R2 - 2*c*x - 2*d*y)^2) <= 0 && 0 <= R1 && 0 <= R2

こうした議論や置換を経ず，同じ結果を直接得ることは Mathematica でも難しく，stochastic な検算に留まるのが現状です．

F3[{x_, y_, a_, b_, c_, d_, R1_, R2_}] = Exists[{x1, y1, x2, y2}, y1*x2 + x1*y2 == y && x1*x2 - y1*y2 == x && (-a + x1)^2 + (-b + y1)^2 == R1 && (-c + x2)^2 + (-d + y2)^2 == R2];  
F5[{x_, y_, a_, b_, c_, d_, R1_, R2_}] = ((-a^2)*c^2 - b^2*c^2 - a^2*d^2 - b^2*d^2 + c^2*R1 + d^2*R1 + a^2*R2 + b^2*R2 - R1*R2 + x^2 + y^2 + b*(2*b*c^2 + 2*b*d^2 - 2*b*R2 + 2*d*x - 2*c*y) + a*(2*a*c^2 + 2*a*d^2 - 2*a*R2 - 2*c*x - 2*d*y))^2 - R1*((2*b*c^2 + 2*b*d^2 - 2*b*R2 + 2*d*x - 2*c*y)^2 + (2*a*c^2 + 2*a*d^2 - 2*a*R2 - 2*c*x - 2*d*y)^2) <= 0 && 0 <= R1 && 0 <= R2;  
F35[L_] := CylindricalDecomposition[Equivalent[F3[L], F5[L]]];

And@@Table[F35@(RandomSample[-10;;10^5,8]/10^5),10^4]//Timing
{174.248398, True}

https://github.com/chriswestbrown/tarski/raw/master/interpreter/ex/GenCircleCircle.txt