Problem of Apollonius: Difference between revisions

From Rosetta Code
Content added Content deleted
Line 387: Line 387:


=={{header|Icon}} and {{header|Unicon}}==
=={{header|Icon}} and {{header|Unicon}}==
This is a translation of the Java version. [[Apollonius-unicon.png|thumb|Solution for Apollonius]]

<lang Icon>link graphics
<lang Icon>link graphics


Revision as of 05:27, 12 May 2011

Task
Problem of Apollonius
You are encouraged to solve this task according to the task description, using any language you may know.

Implement a solution to the Problem of Apollonius (description on wikipedia) which is the problem of finding the circle that is tangent to three specified circles. There is an algebraic solution which is pretty straightforward. The solutions to the example in the code are shown in the image (the red circle is "internally tangent" to all three black circles and the green circle is "externally tangent" to all three black circles).

Two solutions to the problem of apollonius

Ada

apollonius.ads: <lang Ada>package Apollonius is 

  type Point is record
     X, Y : Long_Float := 0.0;
  end record;

  type Circle is record
     Center : Point;
     Radius : Long_Float := 0.0;
  end record;

  type Tangentiality is (External, Internal);

  function Solve_CCC
    (Circle_1, Circle_2, Circle_3 : Circle;
     T1, T2, T3                   : Tangentiality := External)
     return                         Circle;

end Apollonius;</lang> apollonius.adb: <lang Ada>with Ada.Numerics.Generic_Elementary_Functions;

package body Apollonius is 

  package Math is new Ada.Numerics.Generic_Elementary_Functions
    (Long_Float);

  function Solve_CCC
    (Circle_1, Circle_2, Circle_3 : Circle;
     T1, T2, T3                   : Tangentiality := External)
     return                         Circle
  is
     S1 : Long_Float := 1.0;
     S2 : Long_Float := 1.0;
     S3 : Long_Float := 1.0;

     X1 : Long_Float renames Circle_1.Center.X;
     Y1 : Long_Float renames Circle_1.Center.Y;
     R1 : Long_Float renames Circle_1.Radius;

     X2 : Long_Float renames Circle_2.Center.X;
     Y2 : Long_Float renames Circle_2.Center.Y;
     R2 : Long_Float renames Circle_2.Radius;

     X3 : Long_Float renames Circle_3.Center.X;
     Y3 : Long_Float renames Circle_3.Center.Y;
     R3 : Long_Float renames Circle_3.Radius;
  begin
     if T1 = Internal then
        S1 := -S1;
     end if;
     if T2 = Internal then
        S2 := -S2;
     end if;
     if T3 = Internal then
        S3 := -S3;
     end if;

     declare
        V11 : constant Long_Float := 2.0 * X2 - 2.0 * X1;
        V12 : constant Long_Float := 2.0 * Y2 - 2.0 * Y1;
        V13 : constant Long_Float :=
          X1 * X1 - X2 * X2 + Y1 * Y1 - Y2 * Y2 - R1 * R1 + R2 * R2;
        V14 : constant Long_Float := 2.0 * S2 * R2 - 2.0 * S1 * R1;

        V21 : constant Long_Float := 2.0 * X3 - 2.0 * X2;
        V22 : constant Long_Float := 2.0 * Y3 - 2.0 * Y2;
        V23 : constant Long_Float :=
          X2 * X2 - X3 * X3 + Y2 * Y2 - Y3 * Y3 - R2 * R2 + R3 * R3;
        V24 : constant Long_Float := 2.0 * S3 * R3 - 2.0 * S2 * R2;

        W12 : constant Long_Float := V12 / V11;
        W13 : constant Long_Float := V13 / V11;
        W14 : constant Long_Float := V14 / V11;

        W22 : constant Long_Float := V22 / V21 - W12;
        W23 : constant Long_Float := V23 / V21 - W13;
        W24 : constant Long_Float := V24 / V21 - W14;

        P   : constant Long_Float := -W23 / W22;
        Q   : constant Long_Float := W24 / W22;
        M   : constant Long_Float := -W12 * P - W13;
        N   : constant Long_Float := W14 - W12 * Q;

        A   : constant Long_Float := N * N + Q * Q - 1.0;
        B   : constant Long_Float :=
          2.0 * M * N -
            2.0 * N * X1 +
              2.0 * P * Q -
                2.0 * Q * Y1 +
                  2.0 * S1 * R1;
        C   : constant Long_Float :=
          X1 * X1 +
            M * M -
              2.0 * M * X1 +
                P * P +
                  Y1 * Y1 -
                    2.0 * P * Y1 -
                      R1 * R1;

        D   : constant Long_Float := B * B - 4.0 * A * C;
        RS  : constant Long_Float := (-B - Math.Sqrt (D)) / (2.0 * A);
     begin
        return (Center => (X => M + N * RS, Y => P + Q * RS), Radius => RS);
     end;
  end Solve_CCC;

end Apollonius;</lang>

example test_apollonius.adb: <lang Ada>with Ada.Text_IO; with Apollonius;

procedure Test_Apollonius is 

  use Apollonius;
  package Long_Float_IO is new Ada.Text_IO.Float_IO (Long_Float);

  C1 : constant Circle := (Center => (X => 0.0, Y => 0.0), Radius => 1.0);
  C2 : constant Circle := (Center => (X => 4.0, Y => 0.0), Radius => 1.0);
  C3 : constant Circle := (Center => (X => 2.0, Y => 4.0), Radius => 2.0);

  R1 : Circle := Solve_CCC (C1, C2, C3, External, External, External);
  R2 : Circle := Solve_CCC (C1, C2, C3, Internal, Internal, Internal);

begin 

  Ada.Text_IO.Put_Line ("R1:");
  Long_Float_IO.Put (R1.Center.X, Aft => 3, Exp => 0);
  Long_Float_IO.Put (R1.Center.Y, Aft => 3, Exp => 0);
  Long_Float_IO.Put (R1.Radius, Aft => 3, Exp => 0);
  Ada.Text_IO.New_Line;
  Ada.Text_IO.Put_Line ("R2:");
  Long_Float_IO.Put (R2.Center.X, Aft => 3, Exp => 0);
  Long_Float_IO.Put (R2.Center.Y, Aft => 3, Exp => 0);
  Long_Float_IO.Put (R2.Radius, Aft => 3, Exp => 0);
  Ada.Text_IO.New_Line;

end Test_Apollonius;</lang>

output: 

R1:
 2.000 2.100 3.900
R2:
 2.000 0.833 1.167

C

Translation of: Java

<lang c>#include <stdio.h>

  1. include <stdlib.h>
  2. include <math.h>
  1. define EXTERNALLY 1
  2. define INTERNALLY -1

typedef struct { 

   double x, y, r;

} Circle;


void circle_print(Circle c) { 

   printf("Circle[x=%.4lf, y=%.4lf, r=%.4lf]\n", c.x, c.y, c.r);

}

Circle solveApollonius(Circle c1, Circle c2, Circle c3, 

                       int s1, int s2, int s3) {
   const double v11 = 2 * c2.x - 2 * c1.x;
   const double v12 = 2 * c2.y - 2 * c1.y;
   const double v13 = c1.x*c1.x - c2.x*c2.x + c1.y*c1.y -
                      c2.y*c2.y - c1.r*c1.r + c2.r*c2.r;
   const double v14 = 2 * s2 * c2.r - 2 * s1 * c1.r;

   const double v21 = 2*c3.x - 2*c2.x;
   const double v22 = 2*c3.y - 2*c2.y;
   const double v23 = c2.x*c2.x - c3.x*c3.x + c2.y*c2.y -
                      c3.y*c3.y - c2.r*c2.r + c3.r*c3.r;
   const double v24 = 2 * s3 * c3.r - 2 * s2 * c2.r;

   const double w12 = v12 / v11;
   const double w13 = v13 / v11;
   const double w14 = v14 / v11;

   const double w22 = v22 / v21 - w12;
   const double w23 = v23 / v21 - w13;
   const double w24 = v24 / v21 - w14;

   const double P = -w23 / w22;
   const double Q = w24 / w22;
   const double M = -w12*P - w13;
   const double N = w14 - w12 * Q;

   const double a = N * N + Q * Q - 1;
   const double b = 2*M*N - 2*N*c1.x + 2*P*Q - 2*Q*c1.y + 2*s1*c1.r;
   const double c = c1.x*c1.x + M*M - 2*M*c1.x + P*P +
                    c1.y*c1.y - 2*P*c1.y - c1.r*c1.r;

   const double D = b*b -4*a*c;
   const double rs = (-b - sqrt(D)) / (2 * a);

   Circle result = { M + N * rs, P + Q * rs, rs };
   return result;

}

int main() { 

   Circle c1 = { 0.0, 0.0, 1.0 };
   Circle c2 = { 4.0, 0.0, 1.0 };
   Circle c3 = { 2.0, 4.0, 2.0 };

   Circle r1 = solveApollonius(c1, c2, c3,
                               EXTERNALLY, EXTERNALLY, EXTERNALLY);
   circle_print(r1);
   Circle r2 = solveApollonius(c1, c2, c3,
                               INTERNALLY, INTERNALLY, INTERNALLY);
   circle_print(r2);

   return 0;

}</lang> Output: 

Circle[x=2.0000, y=2.1000, r=3.9000]
Circle[x=2.0000, y=0.8333, r=1.1667]

D

Translated and improved from the Java version. <lang d>import std.stdio: writeln; import std.typecons: Tuple; import std.math: sqrt;


alias Tuple!(double, "x", double, "y", double, "r") Circle; enum Tangent { externally, internally }

/** Solves the Problem of Apollonius (finding a circle tangent to three other circles in the plane).

Params: 

 c1 = First circle of the problem.
 c2 = Second circle of the problem.
 c3 = Third circle of the problem.
 t1 = If the solution should be externally or internally tangent to c1.
 t2 = If the solution should be externally or internally tangent to c2.
 t3 = If the solution should be externally or internally tangent to c3.

Returns: The Circle that is tangent to c1, c2 and c3.

  • /

Circle solveApollonius(in Circle c1, in Circle c2, in Circle c3, 

                      in Tangent t1, in Tangent t2, in Tangent t3) {
   const double s1 = (t1 == Tangent.externally) ? 1.0 : -1.0;
   const double s2 = (t2 == Tangent.externally) ? 1.0 : -1.0;
   const double s3 = (t3 == Tangent.externally) ? 1.0 : -1.0;

   const double v11 = 2 * c2.x - 2 * c1.x;
   const double v12 = 2 * c2.y - 2 * c1.y;
   const double v13 = c1.x ^^ 2 - c2.x ^^ 2 +
                      c1.y ^^ 2 - c2.y ^^ 2 -
                      c1.r ^^ 2 + c2.r ^^ 2;
   const double v14 = 2 * s2 * c2.r - 2 * s1 * c1.r;

   const double v21 = 2 * c3.x - 2 * c2.x;
   const double v22 = 2 * c3.y - 2 * c2.y;
   const double v23 = c2.x ^^ 2 - c3.x ^^ 2 +
                      c2.y ^^ 2 - c3.y ^^ 2 -
                      c2.r ^^ 2 + c3.r ^^ 2;
   const double v24 = 2 * s3 * c3.r - 2 * s2 * c2.r;

   const double w12 = v12 / v11;
   const double w13 = v13 / v11;
   const double w14 = v14 / v11;

   const double w22 = v22 / v21 - w12;
   const double w23 = v23 / v21 - w13;
   const double w24 = v24 / v21 - w14;

   const double P = -w23 / w22;
   const double Q =  w24 / w22;
   const double M = -w12 * P - w13;
   const double N =  w14 - w12 * Q;

   const double a = N * N + Q ^^ 2 - 1;
   const double b = 2 * M * N - 2 * N * c1.x +
                    2 * P * Q - 2 * Q * c1.y +
                    2 * s1 * c1.r;
   const double c = c1.x ^^ 2 + M ^^ 2 - 2 * M * c1.x +
                    P ^^ 2 + c1.y ^^ 2 - 2 * P * c1.y - c1.r ^^ 2;

   // find a root of a quadratic equation.
   // This requires the circle centers not to be e.g. colinear
   const double D = b*b - 4*a*c;
   const double rs = (-b - sqrt(D)) / (2 * a);

   return Circle(M + N * rs, P + Q * rs, rs);

}


void main() { 

   auto c1 = Circle(0.0, 0.0, 1.0);
   auto c2 = Circle(4.0, 0.0, 1.0);
   auto c3 = Circle(2.0, 4.0, 2.0);

   alias Tangent.externally te;
   writeln(solveApollonius(c1, c2, c3, te, te, te));

   alias Tangent.internally ti;
   writeln(solveApollonius(c1, c2, c3, ti, ti, ti));

}</lang> Output: 

Tuple!(double,"x",double,"y",double,"r")(2, 2.1, 3.9)
Tuple!(double,"x",double,"y",double,"r")(2, 0.833333, 1.16667)

Fortran

Works with: Fortran version 90 and later

<lang fortran>program Apollonius 

 implicit none

 integer, parameter :: dp = selected_real_kind(15)

 type circle
   real(dp) :: x
   real(dp) :: y
   real(dp) :: radius
 end type
      
 type(circle) :: c1 , c2, c3, r

 c1 = circle(0.0, 0.0, 1.0)
 c2 = circle(4.0, 0.0, 1.0)
 c3 = circle(2.0, 4.0, 2.0)

 write(*, "(a,3f12.8))") "External tangent:", SolveApollonius(c1, c2, c3, 1, 1, 1)
 write(*, "(a,3f12.8))") "Internal tangent:", SolveApollonius(c1, c2, c3, -1, -1, -1)

contains

function SolveApollonius(c1, c2, c3, s1, s2, s3) result(res) 

 type(circle) :: res
 type(circle), intent(in) :: c1, c2, c3
 integer, intent(in) :: s1, s2, s3

 real(dp) :: x1, x2, x3, y1, y2, y3, r1, r2, r3
 real(dp) :: v11, v12, v13, v14
 real(dp) :: v21, v22, v23, v24
 real(dp) :: w12, w13, w14
 real(dp) :: w22, w23, w24
 real(dp) :: p, q, m, n, a, b, c, det
 
 x1 = c1%x; x2 = c2%x; x3 = c3%x
 y1 = c1%y; y2 = c2%y; y3 = c3%y
 r1 = c1%radius; r2 = c2%radius; r3 = c3%radius

 v11 = 2*x2 - 2*x1
 v12 = 2*y2 - 2*y1
 v13 = x1*x1 - x2*x2 + y1*y1 - y2*y2 - r1*r1 + r2*r2
 v14 = 2*s2*r2 - 2*s1*r1

 v21 = 2*x3 - 2*x2
 v22 = 2*y3 - 2*y2
 v23 = x2*x2 - x3*x3 + y2*y2 - y3*y3 - r2*r2 + r3*r3
 v24 = 2*s3*r3 - 2*s2*r2

 w12 = v12/v11
 w13 = v13/v11
 w14 = v14/v11

 w22 = v22/v21-w12
 w23 = v23/v21-w13
 w24 = v24/v21-w14

 p = -w23/w22
 q = w24/w22
 m = -w12*P - w13
 n = w14 - w12*q

 a = n*n + q*q - 1
 b = 2*m*n - 2*n*x1 + 2*p*q - 2*q*y1 + 2*s1*r1
 c = x1*x1 + m*m - 2*m*x1 + p*p + y1*y1 - 2*p*y1 - r1*r1

 det = b*b - 4*a*c
 res%radius = (-b-sqrt(det)) / (2*a)
 res%x = m + n*res%radius
 res%y = p + q*res%radius

end function end program</lang> Output 

External tangent:  2.00000000  2.10000000  3.90000000
Internal tangent:  2.00000000  0.83333333  1.16666667

Icon and Unicon

This is a translation of the Java version. thumb|Solution for Apollonius

<lang Icon>link graphics

record circle(x,y,r) global scale,xoffset,yoffset,yadjust

procedure main()

WOpen("size=400,400") | stop("Unable to open Window") scale := 28 xoffset := WAttrib("width") / 2 yoffset := ( yadjust := WAttrib("height")) / 2


WC(c1 := circle(0,0,1),"black") WC(c2 := circle(4,0,1),"black") WC(c3 := circle(2,4,2),"black") WC(c4 := Apollonius(c1,c2,c3,1,1,1),"green") #/ Expects "Circle[x=2.00,y=2.10,r=3.90]" (green circle in image) WC(c5 := Apollonius(c1,c2,c3,-1,-1,-1),"red") #/ Expects "Circle[x=2.00,y=0.83,r=1.17]" (red circle in image)


WAttrib("fg=blue") DrawLine( 0*scale+xoffset, yadjust-(-1*scale+yoffset), 0*scale+xoffset, yadjust-(4*scale+yoffset) ) DrawLine( -1*scale+xoffset, yadjust-(0*scale+yoffset), 4*scale+xoffset, yadjust-(0*scale+yoffset) ) WDone() end

procedure WC(c,fg) # write and plot circle WAttrib("fg="||fg) DrawCircle(c.x*scale+xoffset, yadjust-(c.y*scale+yoffset), c.r*scale) return write("Circle(x,y,r) := (",c.x,", ",c.y,", ",c.r,")") end

procedure Apollonius(c1,c2,c3,s1,s2,s3) # solve Apollonius 

 v11 := 2.*(c2.x - c1.x)
 v12 := 2.*(c2.y - c1.y)
 v13 := c1.x^2 - c2.x^2 + c1.y^2 - c2.y^2 - c1.r^2 + c2.r^2
 v14 := 2.*(s2*c2.r - s1*c1.r)

 v21 := 2.*(c3.x - c2.x)
 v22 := 2.*(c3.y - c2.y)
 v23 := c2.x^2 - c3.x^2 + c2.y^2 - c3.y^2 - c2.r^2 + c3.r^2
 v24 := 2.*(s3*c3.r - s2*c2.r)

 w12 := v12/v11
 w13 := v13/v11
 w14 := v14/v11

 w22 := v22/v21-w12
 w23 := v23/v21-w13
 w24 := v24/v21-w14

 P := -w23/w22
 Q := w24/w22
 M := -w12*P-w13
 N := w14 - w12*Q

 a := N*N + Q*Q - 1
 b := 2*M*N - 2*N*c1.x + 2*P*Q - 2*Q*c1.y + 2*s1*c1.r
 c := c1.x*c1.x + M*M - 2*M*c1.x + P*P + c1.y*c1.y - 2*P*c1.y - c1.r*c1.r

 #// Find a root of a quadratic equation. This requires the circle centers not to be e.g. colinear
 D := b*b-4*a*c
 rs := (-b-sqrt(D))/(2*a)
 xs := M + N * rs
 ys := P + Q * rs
 return circle(xs,ys,rs)

end</lang>

Output:

Circle(x,y,r) := (0, 0, 1)
Circle(x,y,r) := (4, 0, 1)
Circle(x,y,r) := (2, 4, 2)
Circle(x,y,r) := (2.0, 2.1, 3.9)
Circle(x,y,r) := (2.0, 0.8333333333333333, 1.166666666666667)

J

Solution <lang j>require 'math/misc/amoeba'

NB.*apollonius v solves Apollonius problems NB. y is Cx0 Cy0 R0 Cx1 Cy1 R1 Cx2 Cy2 R2 NB. x are radius scale factors to control which circles are included NB. in the common tangent circle. 1 to surround, _1 to exclude. NB. returns Cxs Cys Rs apollonius =: verb define"1 _ 

1 apollonius y

centers=. 2{."1 y
radii=. x * {:"1 y
goal=. 1e_20                               NB. goal simplex volume
dist=. radii + [: +/"1&.:*: centers -"1 ]  NB. distances to tangents
'soln err'=. ([: +/@:*:@, -/~@dist) f. amoeba goal centers
if. err > 10 * goal do.  return. end.    NB. no solution found
avg=. +/ % #
(, avg@dist) soln

)</lang>

Usage <lang j> ]rctst=: 0 0 1,4 0 1,:2 4 2 NB. Task circles 0 0 1 4 0 1 2 4 2 

  (_1 _1 _1 ,: 1 1 1) apollonius rctst  NB. internally & externally tangent solutions

2 0.83333333 1.1666667 2 2.1 3.9 </lang>

Java

Some lines in this example are too long (more than 80 characters). Please fix the code if it's possible and remove this message.

<lang Java>public class Circle { 

public double[] center;
public double radius;
public Circle(double[] center, double radius)
{
 this.center = center;
 this.radius = radius;
}
public String toString()
{
 return String.format("Circle[x=%.2f,y=%.2f,r=%.2f]",center[0],center[1],radius);
}

}

public class ApolloniusSolver { /** Solves the Problem of Apollonius (finding a circle tangent to three other circles in the plane). 

 * The method uses approximately 68 heavy operations (multiplication, division, square-roots). 
 * @param c1 One of the circles in the problem 
 * @param c2 One of the circles in the problem
 * @param c3 One of the circles in the problem
 * @param s1 An indication if the solution should be externally or internally tangent (+1/-1) to c1
 * @param s2 An indication if the solution should be externally or internally tangent (+1/-1) to c2
 * @param s3 An indication if the solution should be externally or internally tangent (+1/-1) to c3
 * @return The circle that is tangent to c1, c2 and c3. 
 */
public static Circle solveApollonius(Circle c1, Circle c2, Circle c3, int s1, int s2, int s3)
{
 float x1 = c1.center[0];
 float y1 = c1.center[1];
 float r1 = c1.radius;
 float x2 = c2.center[0];
 float y2 = c2.center[1];
 float r2 = c2.radius;
 float x3 = c3.center[0];
 float y3 = c3.center[1];
 float r3 = c3.radius;

 //Currently optimized for fewest multiplications. Should be optimized for readability
 float v11 = 2*x2 - 2*x1;
 float v12 = 2*y2 - 2*y1;
 float v13 = x1*x1 - x2*x2 + y1*y1 - y2*y2 - r1*r1 + r2*r2;
 float v14 = 2*s2*r2 - 2*s1*r1;

 float v21 = 2*x3 - 2*x2;
 float v22 = 2*y3 - 2*y2;
 float v23 = x2*x2 - x3*x3 + y2*y2 - y3*y3 - r2*r2 + r3*r3;
 float v24 = 2*s3*r3 - 2*s2*r2;

 float w12 = v12/v11;
 float w13 = v13/v11;
 float w14 = v14/v11;

 float w22 = v22/v21-w12;
 float w23 = v23/v21-w13;
 float w24 = v24/v21-w14;

 float P = -w23/w22;
 float Q = w24/w22;
 float M = -w12*P-w13;
 float N = w14 - w12*Q;

 float a = N*N + Q*Q - 1;
 float b = 2*M*N - 2*N*x1 + 2*P*Q - 2*Q*y1 + 2*s1*r1;
 float c = x1*x1 + M*M - 2*M*x1 + P*P + y1*y1 - 2*P*y1 - r1*r1;

 // Find a root of a quadratic equation. This requires the circle centers not to be e.g. colinear
 float D = b*b-4*a*c;
 float rs = (-b-Math.sqrt(D))/(2*a);
 float xs = M + N * rs;
 float ys = P + Q * rs;
 return new Circle(new double[]{xs,ys}, rs);
}
public static void main(final String[] args)
{
 Circle c1 = new Circle(new double[]{0,0}, 1);
 Circle c2 = new Circle(new double[]{4,0}, 1);
 Circle c3 = new Circle(new double[]{2,4}, 2);
 System.out.println(solveApollonius(c1,c2,c3,1,1,1));    // Expects "Circle[x=2.00,y=2.10,r=3.90]" (green circle in image)
 System.out.println(solveApollonius(c1,c2,c3,-1,-1,-1)); // Expects "Circle[x=2.00,y=0.83,r=1.17]" (red circle in image)
}

}</lang>

MUMPS

Translation of: Java

<lang MUMPS>APOLLONIUS(CIR1,CIR2,CIR3,S1,S2,S3) 

;Circles are passed in as strings with three parts with a "^" separator in the order x^y^r
;The three circles are CIR1, CIR2, and CIR3
;The S1, S2, and S3 parameters determine if the solution will be internally or externally
;tangent to the circle. (+1 external, -1 internal)
;CIRR is the circle returned in the same format as the input circles
;
;Xn, Yn, and Rn are the values for a circle n - following the precedents from the
;other examples because doing $Pieces would make this confusing to read
NEW X1,X2,X3,Y1,Y2,Y3,R1,R2,R3,RS,V11,V12,V13,V14,V21,V22,V23,V24,W12,W13,W14,W22,W23,W24,P,M,N,Q,A,B,C,D
NEW CIRR
SET X1=$PIECE(CIR1,"^",1),X2=$PIECE(CIR2,"^",1),X3=$PIECE(CIR3,"^",1)
SET Y1=$PIECE(CIR1,"^",2),Y2=$PIECE(CIR2,"^",2),Y3=$PIECE(CIR3,"^",2)
SET R1=$PIECE(CIR1,"^",3),R2=$PIECE(CIR2,"^",3),R3=$PIECE(CIR3,"^",3)
SET V11=(2*X2)-(2*X1)
SET V12=(2*Y2)-(2*Y1)
SET V13=(X1*X1)-(X2*X2)+(Y1*Y1)-(Y2*Y2)-(R1*R1)+(R2*R2)
SET V14=(2*S2*R2)-(2*S1*R1)
SET V21=(2*X3)-(2*X2)
SET V22=(2*Y3)-(2*Y2)
SET V23=(X2*X2)-(X3*X3)+(Y2*Y2)-(Y3*Y3)-(R2*R2)+(R3*R3)
SET V24=(2*S3*R3)-(2*S2*R2)
SET W12=V12/V11
SET W13=V13/V11
SET W14=V14/V11
SET W22=(V22/V21)-W12 ;Parentheses for insurance - MUMPS evaluates left to right
SET W23=(V23/V21)-W13
SET W24=(V24/V21)-W14
SET P=-W23/W22
SET Q=W24/W22
SET M=-(W12*P)-W13
SET N=W14-(W12*Q)
SET A=(N*N)+(Q*Q)-1
SET B=(2*M*N)-(2*N*X1)+(2*P*Q)-(2*Q*Y1)+(2*S1*R1)
SET C=(X1*X1)+(M*M)+(2*M*X1)+(P*P)+(Y1*Y1)-(2*P*Y1)-(R1*R1)
SET D=(B*B)-(4*A*C)
SET RS=(-B-(D**.5))/(2*A)
SET $PIECE(CIRR,"^",1)=M+(N*RS)
SET $PIECE(CIRR,"^",2)=P+(Q*RS)
SET $PIECE(CIRR,"^",3)=RS
KILL X1,X2,X3,Y1,Y2,Y3,R1,R2,R3,RS,V11,V12,V13,V14,V21,V22,V23,V24,W12,W13,W14,W22,W23,W24,P,M,N,Q,A,B,C,D
QUIT CIRR</lang>

In use:

USER>WRITE C1
0^0^1
USER>WRITE C2
4^0^1
USER>WRITE C3
2^4^2
USER>WRITE $$APOLLONIUS^ROSETTA(C1,C2,C3,1,1,1)
2^2.1^3.9
USER>WRITE $$APOLLONIUS^ROSETTA(C1,C2,C3,-1,-1,-1)
2^.833333333333333333^1.166666666666666667

OCaml

Translation of: C

<lang ocaml>type point = { x:float; y:float } type circle = { 

 center: point;
 radius: float;

}

let new_circle ~x ~y ~r = 

 { center = { x=x; y=y };
   radius = r }

let print_circle ~c = 

 Printf.printf "Circle(x=%.2f, y=%.2f, r=%.2f)\n"
   c.center.x c.center.y c.radius

let defxyr c = 

 (c.center.x,
  c.center.y,
  c.radius)

let solve_apollonius ~c1 ~c2 ~c3 

                    ~s1 ~s2 ~s3 =
 let ( * ) = ( *. ) in
 let ( / ) = ( /. ) in
 let ( + ) = ( +. ) in
 let ( - ) = ( -. ) in

 let x1, y1, r1 = defxyr c1
 and x2, y2, r2 = defxyr c2
 and x3, y3, r3 = defxyr c3 in

 let v11 = 2.0 * x2 - 2.0 * x1
 and v12 = 2.0 * y2 - 2.0 * y1
 and v13 = x1*x1 - x2*x2 + y1*y1 - y2*y2 - r1*r1 + r2*r2
 and v14 = (2.0 * s2 * r2) - (2.0 * s1 * r1)

 and v21 = 2.0 * x3 - 2.0 * x2
 and v22 = 2.0 * y3 - 2.0 * y2
 and v23 = x2*x2 - x3*x3 + y2*y2 - y3*y3 - r2*r2 + r3*r3
 and v24 = (2.0 * s3 * r3) - (2.0 * s2 * r2) in

 let w12 = v12 / v11
 and w13 = v13 / v11
 and w14 = v14 / v11 in

 let w22 = v22 / v21 - w12
 and w23 = v23 / v21 - w13
 and w24 = v24 / v21 - w14 in

 let p = -. w23 / w22
 and q = w24 / w22 in
 let m = -. w12 * p - w13
 and n = w14 - w12 * q in

 let a = n*n + q*q - 1.0
 and b = 2.0*m*n - 2.0*n*x1 + 2.0*p*q - 2.0*q*y1 + 2.0*s1*r1
 and c = x1*x1 + m*m - 2.0*m*x1 + p*p + y1*y1 - 2.0*p*y1 - r1*r1 in

 let d = b * b - 4.0 * a * c in
 let rs = (-. b - (sqrt d)) / (2.0 * a) in

 let xs = m + n * rs
 and ys = p + q * rs in

 (new_circle xs ys rs)

let () = 

 let c1 = new_circle 0.0 0.0 1.0
 and c2 = new_circle 4.0 0.0 1.0
 and c3 = new_circle 2.0 4.0 2.0 in

 let r1 = solve_apollonius c1 c2 c3 1.0 1.0 1.0 in
 print_circle r1;

 let r2 = solve_apollonius c1 c2 c3 (-1.) (-1.) (-1.) in
 print_circle r2;
</lang>


PureBasic

Translation of: Java

<lang PureBasic>Structure Circle 

 XPos.f
 YPos.f
 Radius.f

EndStructure

Procedure ApolloniusSolver(*c1.Circle,*c2.Circle,*c3.Circle, s1, s2, s3) 

 Define.f  ; This tells the compiler that all non-specified new variables
           ; should be of float type (.f).
 x1=*c1\XPos:  y1=*c1\YPos:  r1=*c1\Radius
 x2=*c2\XPos:  y2=*c2\YPos:  r2=*c2\Radius  
 x3=*c3\XPos:  y3=*c3\YPos:  r3=*c3\Radius
 
 v11 = 2*x2 - 2*x1
 v12 = 2*y2 - 2*y1
 v13 = x1*x1 - x2*x2 + y1*y1 - y2*y2 - r1*r1 + r2*r2
 v14 = 2*s2*r2 - 2*s1*r1
 
 v21 = 2*x3 - 2*x2
 v22 = 2*y3 - 2*y2
 v23 = x2*x2 - x3*x3 + y2*y2 - y3*y3 - r2*r2 + r3*r3
 v24 = 2*s3*r3 - 2*s2*r2
 
 w12 = v12/v11
 w13 = v13/v11
 w14 = v14/v11
 
 w22 = v22/v21-w12
 w23 = v23/v21-w13
 w24 = v24/v21-w14
 
 P = -w23/w22
 Q =  w24/w22
 M = -w12*P-w13
 N =  w14-w12*Q
 
 a = N*N + Q*Q - 1
 b = 2*M*N - 2*N*x1 + 2*P*Q - 2*Q*y1 + 2*s1*r1
 c = x1*x1 + M*M - 2*M*x1 + P*P + y1*y1 - 2*P*y1 - r1*r1
 
 D= b*b - 4*a*c
 
 Define *result.Circle=AllocateMemory(SizeOf(Circle))
 ; Allocate memory for a returned Structure of type Circle.
 ; This memory should be freed later but if not, PureBasic’s
 ; internal framework will do so when the program shuts down.
 If *result
   *result\Radius=(-b-Sqr(D))/(2*a)
   *result\XPos  =M+N * *result\Radius
   *result\YPos  =P+Q * *result\Radius
 EndIf
 ProcedureReturn *result ; Sending back a pointer

EndProcedure

If OpenConsole() 

 Define.Circle c1, c2, c3
 Define *c.Circle  ; '*c' is defined as a pointer to a circle-structure.
 c1\Radius=1
 c2\XPos=4:  c2\Radius=1
 c3\XPos=2:  c3\YPos=4:  c3\Radius=2
 
 *c=ApolloniusSolver(@c1, @c2, @c3, 1, 1, 1)
 If *c ; Verify that *c got allocated
   PrintN("Circle [x="+StrF(*c\XPos,2)+", y="+StrF(*c\YPos,2)+", r="+StrF(*c\Radius,2)+"]")
   FreeMemory(*c)  ; We are done with *c for the first calculation
 EndIf        
 
 *c=ApolloniusSolver(@c1, @c2, @c3,-1,-1,-1)
 If *c
   PrintN("Circle [x="+StrF(*c\XPos,2)+", y="+StrF(*c\YPos,2)+", r="+StrF(*c\Radius,2)+"]")
   FreeMemory(*c) 
 EndIf
 Print("Press ENTER to exit"): Input()

EndIf</lang> 

Circle [x=2.00, y=2.10, r=3.90]
Circle [x=2.00, y=0.83, r=1.17]
Press ENTER to exit

Python

Translation of: Java

. Although a Circle class is defined, the solveApollonius function is defined in such a way that any three valued tuple or list could be used instead of c1, c2, and c3. The function calls near the end use instances of the Circle class, whereas the docstring shows how the same can be achieved using simple tuples. (And also serves as a simple doctest)

<lang python> from collections import namedtuple import math

Circle = namedtuple('Circle', 'x, y, r')

def solveApollonius(c1, c2, c3, s1, s2, s3): 

   
   >>> solveApollonius((0, 0, 1), (4, 0, 1), (2, 4, 2), 1,1,1)
   Circle(x=2.0, y=2.1, r=3.9)
   >>> solveApollonius((0, 0, 1), (4, 0, 1), (2, 4, 2), -1,-1,-1)
   Circle(x=2.0, y=0.8333333333333333, r=1.1666666666666667) 
   
   x1, y1, r1 = c1
   x2, y2, r2 = c2
   x3, y3, r3 = c3

   v11 = 2*x2 - 2*x1
   v12 = 2*y2 - 2*y1
   v13 = x1*x1 - x2*x2 + y1*y1 - y2*y2 - r1*r1 + r2*r2
   v14 = 2*s2*r2 - 2*s1*r1

   v21 = 2*x3 - 2*x2
   v22 = 2*y3 - 2*y2
   v23 = x2*x2 - x3*x3 + y2*y2 - y3*y3 - r2*r2 + r3*r3
   v24 = 2*s3*r3 - 2*s2*r2

   w12 = v12/v11
   w13 = v13/v11
   w14 = v14/v11

   w22 = v22/v21-w12
   w23 = v23/v21-w13
   w24 = v24/v21-w14

   P = -w23/w22
   Q = w24/w22
   M = -w12*P-w13
   N = w14 - w12*Q

   a = N*N + Q*Q - 1
   b = 2*M*N - 2*N*x1 + 2*P*Q - 2*Q*y1 + 2*s1*r1
   c = x1*x1 + M*M - 2*M*x1 + P*P + y1*y1 - 2*P*y1 - r1*r1

   # Find a root of a quadratic equation. This requires the circle centers not to be e.g. colinear
   D = b*b-4*a*c
   rs = (-b-math.sqrt(D))/(2*a)

   xs = M+N*rs
   ys = P+Q*rs

   return Circle(xs, ys, rs)

if __name__ == '__main__': 

   c1, c2, c3 = Circle(0, 0, 1), Circle(4, 0, 1), Circle(2, 4, 2)
   print(solveApollonius(c1, c2, c3, 1, 1, 1))    #Expects "Circle[x=2.00,y=2.10,r=3.90]" (green circle in image)
   print(solveApollonius(c1, c2, c3, -1, -1, -1)) #Expects "Circle[x=2.00,y=0.83,r=1.17]" (red circle in image)</lang>

Sample Output 

Circle(x=2.0, y=2.1, r=3.9)
Circle(x=2.0, y=0.8333333333333333, r=1.1666666666666667)

=begin Problem of Apollonius http://rosettacode.org/wiki/Problem_of_Apollonius

Ruby

Translation of: Java

<lang ruby>class Circle 

 def initialize(x, y, r)
   @x, @y, @r = [x, y, r].map(&:to_f)
 end
 attr_reader :x, :y, :r

 def self.apollonius(c1, c2, c3, s1=1, s2=1, s3=1)
   x1, y1, r1 = [c1.x, c1.y, c1.r]
   x2, y2, r2 = [c2.x, c2.y, c2.r]
   x3, y3, r3 = [c3.x, c3.y, c3.r]

   v11 = 2*x2 - 2*x1
   v12 = 2*y2 - 2*y1
   v13 = x1**2 - x2**2 + y1**2 - y2**2 - r1**2 + r2**2
   v14 = 2*s2*r2 - 2*s1*r1

   v21 = 2*x3 - 2*x2
   v22 = 2*y3 - 2*y2
   v23 = x2**2 - x3**2 + y2**2 - y3**2 - r2**2 + r3**2
   v24 = 2*s3*r3 - 2*s2*r2

   w12 = v12/v11
   w13 = v13/v11
   w14 = v14/v11

   w22 = v22/v21 - w12
   w23 = v23/v21 - w13
   w24 = v24/v21 - w14

   p = -w23/w22
   q = w24/w22
   m = -w12*p - w13
   n = w14 - w12*q

   a = n**2 + q**2 - 1
   b = 2*m*n - 2*n*x1 + 2*p*q - 2*q*y1 + 2*s1*r1
   c = x1**2 + m**2 - 2*m*x1 + p**2 + y1**2 - 2*p*y1 - r1**2

   d = b**2 - 4*a*c
   rs = (-b - Math.sqrt(d)) / (2*a)
   xs = m + n*rs
   ys = p + q*rs

   return self.new(xs, ys, rs)
 end

end


p c1 = Circle.new(0, 0, 1) p c2 = Circle.new(2, 4, 2) p c3 = Circle.new(4, 0, 1)

p Circle.apollonius(c1, c2, c3) p Circle.apollonius(c1, c2, c3, -1, -1, -1)</lang>

outputs 

#<Circle:0x1012e058 @x=0.0, @y=0.0, @r=1.0>
#<Circle:0x1012deb4 @x=2.0, @y=4.0, @r=2.0>
#<Circle:0x1012dd48 @x=4.0, @y=0.0, @r=1.0>
#<Circle:0x1012d0b4 @x=2.0, @y=2.1, @r=3.9>
#<Circle:0x1012c474 @x=2.0, @y=0.833333333333333, @r=1.16666666666667>

Tcl

Translation of: Java
Works with: Tcl version 8.6

or

Library: TclOO

<lang tcl>package require TclOO; # Just so we can make a circle class

oo::class create circle { 

   variable X Y Radius
   constructor {x y radius} {

namespace import ::tcl::mathfunc::double set X [double $x]; set Y [double $y]; set Radius [double $radius] 

   }
   method values {} {list $X $Y $Radius}
   method format {} {

format "Circle\[o=(%.2f,%.2f),r=%.2f\]" $X $Y $Radius 

   }

}

proc solveApollonius {c1 c2 c3 {s1 1} {s2 1} {s3 1}} { 

   if {abs($s1)!=1||abs($s2)!=1||abs($s3)!=1} {

error "wrong sign; must be 1 or -1" 

   }

   lassign [$c1 values] x1 y1 r1
   lassign [$c2 values] x2 y2 r2
   lassign [$c3 values] x3 y3 r3

   set v11 [expr {2*($x2 - $x1)}]
   set v12 [expr {2*($y2 - $y1)}]
   set v13 [expr {$x1**2 - $x2**2 + $y1**2 - $y2**2 - $r1**2 + $r2**2}]
   set v14 [expr {2*($s2*$r2 - $s1*$r1)}]

   set v21 [expr {2*($x3 - $x2)}]
   set v22 [expr {2*($y3 - $y2)}]
   set v23 [expr {$x2**2 - $x3**2 + $y2**2 - $y3**2 - $r2**2 + $r3**2}]
   set v24 [expr {2*($s3*$r3 - $s2*$r2)}]

   set w12 [expr {$v12 / $v11}]
   set w13 [expr {$v13 / $v11}]
   set w14 [expr {$v14 / $v11}]

   set w22 [expr {$v22 / $v21 - $w12}]
   set w23 [expr {$v23 / $v21 - $w13}]
   set w24 [expr {$v24 / $v21 - $w14}]

   set P [expr {-$w23 / $w22}]
   set Q [expr {$w24 / $w22}]
   set M [expr {-$w12 * $P - $w13}]
   set N [expr {$w14 - $w12 * $Q}]

   set a [expr {$N**2 + $Q**2 - 1}]
   set b [expr {2*($M*$N - $N*$x1 + $P*$Q - $Q*$y1 + $s1*$r1)}]
   set c [expr {($x1-$M)**2 + ($y1-$P)**2 - $r1**2}]

   set rs [expr {(-$b - sqrt($b**2 - 4*$a*$c)) / (2*$a)}]
   set xs [expr {$M + $N*$rs}]
   set ys [expr {$P + $Q*$rs}]

   return [circle new $xs $ys $rs]

}</lang> Demonstration code: <lang tcl>set c1 [circle new 0 0 1] set c2 [circle new 4 0 1] set c3 [circle new 2 4 2] set sA [solveApollonius $c1 $c2 $c3] set sB [solveApollonius $c1 $c2 $c3 -1 -1 -1] puts [$sA format] puts [$sB format]</lang> Output: 

Circle[o=(2.00,2.10),r=3.90]
Circle[o=(2.00,0.83),r=1.17]

Note that the Tcl code uses the ** (exponentiation) operator to shorten and simplify some operations, and that the circle class is forcing the interpretation of every circle's coordinates as double-precision floating-point numbers.

Retrieved from "https://rosettacode.org/wiki/Problem_of_Apollonius?oldid=154754"