problem_of_apollonius.sf

#!/usr/bin/ruby

class Circle(x,y,r) {
   method to_s { "Circle(#{x}, #{y}, #{r})" };
}

func solve_apollonius(Array c, Array s) -> Circle {

    var 𝑣11 = (2*c[1].x - 2*c[0].x);
    var 𝑣12 = (2*c[1].y - 2*c[0].y);
    var 𝑣13 = (c[0].x**2 - c[1].x**2 + c[0].y**2 - c[1].y**2 - c[0].r**2 + c[1].r**2);
    var 𝑣14 = (2*s[1]*c[1].r - 2*s[0]*c[0].r);

    var 𝑣21 = (2*c[2].x - 2*c[1].x);
    var 𝑣22 = (2*c[2].y - 2*c[1].y);
    var 𝑣23 = (c[1].x**2 - c[2].x**2 + c[1].y**2 - c[2].y**2 - c[1].r**2 + c[2].r**2);
    var 𝑣24 = (2*s[2]*c[2].r - 2*s[1]*c[1].r);

    var 𝑤12 = (𝑣12 / 𝑣11);
    var 𝑤13 = (𝑣13 / 𝑣11);
    var 𝑤14 = (𝑣14 / 𝑣11);

    var 𝑤22 = (𝑣22/𝑣21 - 𝑤12);
    var 𝑤23 = (𝑣23/𝑣21 - 𝑤13);
    var 𝑤24 = (𝑣24/𝑣21 - 𝑤14);

    var 𝑃 = (-𝑤23 / 𝑤22);
    var 𝑄 = (𝑤24 / 𝑤22);
    var 𝑀 = ((-𝑤12)*𝑃 - 𝑤13);
    var 𝑁 = (𝑤14 - 𝑤12*𝑄);

    var 𝑎 = (𝑁**2 + 𝑄**2 - 1);
    var 𝑏 = (2*𝑀*𝑁 - 2*𝑁*c[0].x + 2*𝑃*𝑄 - 2*𝑄*c[0].y + 2*s[0]*c[0].r);
    var 𝑐 = (c[0].x**2 + 𝑀**2 - 2*𝑀*c[0].x + 𝑃**2 + c[0].y**2 - 2*𝑃*c[0].y - c[0].r**2);

    var 𝐷 = (𝑏**2 - 4*𝑎*𝑐);
    var rs = ((-𝑏 - 𝐷.sqrt.rat_approx) / 2*𝑎);

    var xs = (𝑀 + 𝑁*rs);
    var ys = (𝑃 + 𝑄*rs);

    Circle(xs, ys, rs);
}

var c = [Circle(0, 0, 1), Circle(4, 0, 1), Circle(2, 4, 2)];

var a = solve_apollonius(c, %n<1 1 1>);
var b = solve_apollonius(c, %n<-1 -1 -1>);

say a;
say b;

a.x == 2     || die "error";
a.y == 21/10 || die "error";
a.r == 39/10 || die "error";
b.x == 2     || die "error";
b.r == 7/6   || die "error";