ApollonianGasket.java

/**
 * Apollonian gasket.
 * Uses signed curvatures (outer negative) and chooses the correct ± sqrt branch by testing which candidate satisfies tangency best.
 */
public class ApollonianGasket {
    private static final int CANVAS = 900;
    private static final double MIN_RADIUS = 0.0005; // Lower = more recursion levels

    private static void drawCircle(double x, double y, double r) {
        StdDraw.circle(x, y, r);
    }

    // complex multiply, add, scale
    private static double[] cmul(double ar, double ai, double br, double bi) {
        return new double[]{ ar*br - ai*bi, ar*bi + ai*br };
    }
    private static double[] cadd(double ar, double ai, double br, double bi) {
        return new double[]{ ar+br, ai+bi };
    }
    private static double[] cscale(double ar, double ai, double s) {
        return new double[]{ ar*s, ai*s };
    }

    // principal complex sqrt
    private static double[] csqrt(double x, double y) {
        double mag = Math.hypot(x, y);
        double re = Math.sqrt((mag + x) / 2.0);
        double im = Math.signum(y) * Math.sqrt(Math.max(0.0, (mag - x) / 2.0));
        if (Double.isNaN(im)) im = 0.0;
        return new double[]{re, im};
    }

    // pick best of two z4 candidates by tangency error to parents
    private static double pickBestCenterAndDraw(
            double x1, double y1, double r1, double k1,
            double x2, double y2, double r2, double k2,
            double x3, double y3, double r3, double k3,
            double k4, double[] numerPlus, double[] numerMinus) {

        double r4 = 1.0 / Math.abs(k4);

        // candidate centers
        double[] z4p = { numerPlus[0] / k4, numerPlus[1] / k4 };
        double[] z4m = { numerMinus[0] / k4, numerMinus[1] / k4 };

        // error function: sum_i | |z4 - zi| - (ri + r4) |
        double errP = 0.0;
        double errM = 0.0;

        double[][] parents = new double[][] {
            {x1, y1, r1}, {x2, y2, r2}, {x3, y3, r3}
        };

        for (int i = 0; i < 3; i++) {
            double px = parents[i][0], py = parents[i][1], pr = parents[i][2];

            double dp = Math.hypot(z4p[0] - px, z4p[1] - py);
            errP += Math.abs(dp - (pr + r4));

            double dm = Math.hypot(z4m[0] - px, z4m[1] - py);
            errM += Math.abs(dm - (pr + r4));
        }

        double x4, y4;
        if (errP <= errM) {
            x4 = z4p[0]; y4 = z4p[1];
        } else {
            x4 = z4m[0]; y4 = z4m[1];
        }

        drawCircle(x4, y4, r4);
        return r4;
    }

    private static void recurse(
            double x1, double y1, double r1, double k1,
            double x2, double y2, double r2, double k2,
            double x3, double y3, double r3, double k3) {

        // curvature k4 using + branch for small inner circle
        double sumk = k1 + k2 + k3;
        double under = k1*k2 + k2*k3 + k3*k1;
        if (under < 0 && Math.abs(under) < 1e-16) under = 0;
        double k4 = sumk + 2.0 * Math.sqrt(Math.max(0.0, under));
        double r4 = 1.0 / Math.abs(k4);
        if (r4 < MIN_RADIUS) return;

        // complex zi
        double[] z1 = {x1, y1}, z2 = {x2, y2}, z3 = {x3, y3};

        double[] z1z2 = cmul(z1[0], z1[1], z2[0], z2[1]);
        double[] z2z3 = cmul(z2[0], z2[1], z3[0], z3[1]);
        double[] z3z1 = cmul(z3[0], z3[1], z1[0], z1[1]);

        double[] S = new double[] {
            k1*k2 * z1z2[0] + k2*k3 * z2z3[0] + k3*k1 * z3z1[0],
            k1*k2 * z1z2[1] + k2*k3 * z2z3[1] + k3*k1 * z3z1[1]
        };

        double[] sqrtS = csqrt(S[0], S[1]);

        // sum k_i z_i
        double[] sumkz = new double[] {
            k1*z1[0] + k2*z2[0] + k3*z3[0],
            k1*z1[1] + k2*z2[1] + k3*z3[1]
        };

        // two numerators: +2*sqrtS and -2*sqrtS
        double[] numerPlus  = cadd(sumkz[0], sumkz[1],  2.0*sqrtS[0],  2.0*sqrtS[1]);
        double[] numerMinus = cadd(sumkz[0], sumkz[1], -2.0*sqrtS[0], -2.0*sqrtS[1]);

        // pick the candidate that best satisfies tangency, draw it, get radius
        double rChosen = pickBestCenterAndDraw(
            x1,y1,r1,k1, x2,y2,r2,k2, x3,y3,r3,k3, k4, numerPlus, numerMinus);

        // compute k4 signed: follow sign of k4 computed earlier (outer negative handled by callers)
        double signedK4 = k4;

        // recurse on triples with the chosen child — we need the child's center to pass to further recursion,
        // recompute chosen center quickly here
        double[] z4p = { numerPlus[0] / k4, numerPlus[1] / k4 };
        double[] z4m = { numerMinus[0] / k4, numerMinus[1] / k4 };

        // choose same way as before (recompute errors)
        double errP = 0, errM = 0;
        double[][] parents = new double[][] {
            {x1,y1,r1}, {x2,y2,r2}, {x3,y3,r3}
        };
        for (int i=0;i<3;i++){
            double px=parents[i][0], py=parents[i][1], pr=parents[i][2];
            errP += Math.abs(Math.hypot(z4p[0]-px, z4p[1]-py) - (pr + rChosen));
            errM += Math.abs(Math.hypot(z4m[0]-px, z4m[1]-py) - (pr + rChosen));
        }
        double x4 = (errP <= errM) ? z4p[0] : z4m[0];
        double y4 = (errP <= errM) ? z4p[1] : z4m[1];

        // recurse
        recurse(x1,y1,r1,k1, x2,y2,r2,k2, x4,y4,rChosen,signedK4);
        recurse(x1,y1,r1,k1, x3,y3,r3,k3, x4,y4,rChosen,signedK4);
        recurse(x2,y2,r2,k2, x3,y3,r3,k3, x4,y4,rChosen,signedK4);
    }

    public static void main(String[] args) {
        StdDraw.setCanvasSize(CANVAS, CANVAS);
        StdDraw.setXscale(0,1);
        StdDraw.setYscale(0,1);

        double R_outer = 0.48;
        double ox = 0.5, oy = 0.5;
        drawCircle(ox, oy, R_outer);
        double k_outer = -1.0 / R_outer; // outer negative

        double r_inner = R_outer / (1.0 + 2.0 / Math.sqrt(3.0));
        double k_inner = 1.0 / r_inner;
        double d = R_outer - r_inner;

        double angleTop = Math.toRadians(90);
        double angleBL  = Math.toRadians(210);
        double angleBR  = Math.toRadians(330);

        double cx1 = ox + d * Math.cos(angleTop);
        double cy1 = oy + d * Math.sin(angleTop);

        double cx2 = ox + d * Math.cos(angleBL);
        double cy2 = oy + d * Math.sin(angleBL);

        double cx3 = ox + d * Math.cos(angleBR);
        double cy3 = oy + d * Math.sin(angleBR);

        drawCircle(cx1, cy1, r_inner);
        drawCircle(cx2, cy2, r_inner);
        drawCircle(cx3, cy3, r_inner);

        // triple of inner circles (all positive curvature)
        recurse(cx1, cy1, r_inner, k_inner, cx2, cy2, r_inner, k_inner, cx3, cy3, r_inner, k_inner);

        // triples that include outer circle (pass outer negative curvature)
        recurse(cx1, cy1, r_inner, k_inner, cx2, cy2, r_inner, k_inner, ox, oy, R_outer, k_outer);
        recurse(cx1, cy1, r_inner, k_inner, cx3, cy3, r_inner, k_inner, ox, oy, R_outer, k_outer);
        recurse(cx2, cy2, r_inner, k_inner, cx3, cy3, r_inner, k_inner, ox, oy, R_outer, k_outer);
    }
}