Problem of Apollonius: Difference between revisions

From Rosetta Code
Content added Content deleted
(→‎{{header|Python}}: Add comment.)
m (→‎{{header|Java}}: lines too long)
Line 98: Line 98:


=={{header|Java}}==
=={{header|Java}}==
{{lines too long}}
<lang Java>public class Circle{
<lang Java>public class Circle{
public double[] center;
public double[] center;

Revision as of 21:00, 14 August 2010

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

C

Translation of: Java

<lang c>#include <stdio.h>

  1. include <stdlib.h>
  2. include <math.h>

typedef struct {

 struct { double x, y; } center;
 double radius;

} circle;


circle *circle_new(double x, double y, double rs) {

 circle *r = malloc(sizeof(circle));
 if ( r == NULL ) return NULL;
 r->center.x = x;
 r->center.y = y;
 r->radius = rs;
 return r;

}

void circle_print(circle *c) {

 printf("Circle[x=%.2lf,y=%.2lf,r=%.2lf]\n",
        c->center.x, c->center.y, c->radius);

}

  1. define DEFXYR(N) double x##N = c##N->center.x; \
                 double y##N = c##N->center.y; \

double r##N = c##N->radius; circle *solveApollonius(circle *c1, circle *c2, circle *c3,

                       int s1, int s2, int s3)

{

   DEFXYR(1); DEFXYR(2); DEFXYR(3);
   double v11 = 2*x2 - 2*x1;
   double v12 = 2*y2 - 2*y1;
   double v13 = x1*x1 - x2*x2 + y1*y1 - y2*y2 - r1*r1 + r2*r2;
   double v14 = 2*s2*r2 - 2*s1*r1;

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

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

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

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

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

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

   double xs = M+N*rs;
   double ys = P+Q*rs;
   circle *r = malloc(sizeof(circle));
   r->center.x = xs;
   r->center.y = ys;
   r->radius = rs;
   return r;

}

int main() {

 circle *c1 = circle_new(0.0,0.0,1.0);
 circle *c2 = circle_new(4.0,0.0,1.0);
 circle *c3 = circle_new(2.0,4.0,2.0);
 
 circle *r1 = solveApollonius(c1, c2, c3, 1, 1, 1);
 circle_print(r1);
 circle *r2 = solveApollonius(c1, c2, c3, -1, -1, -1);
 circle_print(r2);
 
 free(c1); free(c2); free(c3);
 free(r1); free(r2);

}</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(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>

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)

Tcl

Translation of: Java

<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.