Measuring the Philippine Coastline.

// Traverse grid cells intersected by segment a -> b (Amanatides-Woo DDA)
static void traverse(hset *h, pt a, pt b, double s) {
    long ix = (long)floor(a.x / s), iy = (long)floor(a.y / s);
    long ix1 = (long)floor(b.x / s), iy1 = (long)floor(b.y / s);
    hs_insert(h, cell_key(ix, iy));
    
    double dx = b.x - a.x, dy = b.y - a.y;
    int stepx = dx > 0 ? 1 : -1, stepy = dy > 0 ? 1 : -1;
    double tmx = dx != 0 ? ((stepx > 0 ? (ix + 1) * s - a.x : ix * s - a.x) / dx) : 1e300;
    double tmy = dy != 0 ? ((stepy > 0 ? (iy + 1) * s - a.y : iy * s - a.y) / dy) : 1e300;
    double tdx = dx != 0 ? s / fabs(dx) : 1e300;
    double tdy = dy != 0 ? s / fabs(dy) : 1e300;
    
    // March across cell boundaries until reaching endpoint
    while (ix != ix1 || iy != iy1) {
        if (tmx < tmy) { if (tmx >= 1.0) break; ix += stepx; tmx += tdx; }
        else           { if (tmy >= 1.0) break; iy += stepy; tmy += tdy; }
        hs_insert(h, cell_key(ix, iy));
    }
}

// Compass divider walk: step along shoreline using chord radius r
static double yardstick(double r) {
    pt a = P[0];
    int i = 0;
    double steps = 0;
    while (i < n) {
        // Find first vertex outside current radius
        int j = i + 1;
        while (j <= n && dist(P[j % n], a) < r) j++;
        if (j > n) { steps += dist(a, P[0]) / r; break; }
        
        // Solve exact circle-segment intersection
        pt p = (j - 1 == i) ? a : P[j - 1];
        pt q = P[j % n];
        double Ax = p.x - a.x, Ay = p.y - a.y;
        double Bx = q.x - p.x, By = q.y - p.y;
        double B2 = Bx*Bx + By*By, AB = Ax*Bx + Ay*By, C = Ax*Ax + Ay*Ay - r*r;
        double disc = fmax(AB*AB - B2*C, 0.0);
        double t = B2 > 0 ? (-AB + sqrt(disc)) / B2 : 0;
        
        a.x = p.x + fmin(fmax(t, 0.0), 1.0) * Bx;
        a.y = p.y + fmin(fmax(t, 0.0), 1.0) * By;
        steps += 1.0;
        i = j - 1;
    }
    return steps;
}

// Bending yardstick: swing radius r; pivot if chord cuts through land
static double bend_walk(double r) {
    pt a = P[0];
    int i = 0; double steps = 0;
    while (i < n) {
        pt hit; int j; double rem;
        if (!arc_hit(a, r, a, i, &hit, &j, &rem)) { steps += rem / r; break; }
        
        // Check if intermediate points lie on both sides of the chord
        double cx = hit.x - a.x, cy = hit.y - a.y, clen = hypot(cx, cy);
        int pos = 0, neg = 0, hinge = -1; double maxdev = 0;
        for (int v = i + 1; v <= j; v++) {
            pt w = P[v % n];
            double s = (cx * (w.y - a.y) - cy * (w.x - a.x)) / (clen > 0 ? clen : 1);
            if (s > 1e-12 * r) pos = 1;
            if (s < -1e-12 * r) neg = 1;
            if (fabs(s) > maxdev) { maxdev = fabs(s); hinge = v; }
        }
        // If blocked: bend at the furthest vertex and spend remaining length
        if (pos && neg && hinge > i && hinge <= j) {
            double used = dist(a, P[hinge % n]);
            double left = r - used;
            if (left > 1e-12 * r) {
                pt hit2; int j2;
                if (arc_hit(P[hinge % n], left, P[hinge % n], hinge, &hit2, &j2, &rem)) {
                    a = hit2; i = j2; steps += 1.0; continue;
                }
            }
        }
        a = hit; i = j; steps += 1.0;
    }
    return steps;
}