Language: C++

#include <algorithm>
#include <cassert>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <vector>


#define EPSILON 1e-9


class point_t;


class vector_t {
public:
  double x;
  double y;
  double z;
  
public:
  vector_t(double x = 0, double y = 0, double z = 0) : x(x), y(y), z(z) {};
  
  vector_t(const point_t& p);

public:
  double operator[](int i) const {
    return (i == 0) ? x : (i == 1) ? y : z;
  };
  
  double& operator[](int i) {
    return (i == 0) ? x : (i == 1) ? y : z;
  };

public:
  vector_t operator-() const {
    return vector_t(-x, -y, -z);
  };
  
  vector_t operator+(const vector_t& v) const {
    return vector_t(x + v.x, y + v.y, z + v.z);
  };
  
  double operator*(const vector_t& v) const {
    return x * v.x + y * v.y + z * v.z;
  };
  
  vector_t operator*(double a) const {
    return vector_t(x * a, y * a, z * a);
  };
  
  vector_t operator/(double a) const {
    return vector_t(x / a, y / a, z / a);
  };
  
  vector_t operator^(const vector_t& v) const {
    return vector_t(y * v.z - z * v.y,
                    z * v.x - x * v.z,
                    x * v.y - y * v.x);
  };
    
  vector_t operator-=(const vector_t& v) {
    x -= v.x;
    y -= v.y;
    z -= v.z;
    
    return *this;
  };

public:
  double length() const {
    return sqrt(*this * *this);
  };

  vector_t normal() const {
    return *this / length();
  };
};


class point_t {
public:
  double x;
  double y;
  double z;
  double w;

public:
  point_t(double x = 0, double y = 0, double z = 0, double w = 1) : x(x), y(y), z(z), w(w) {};

  point_t(const vector_t& v) : x(v.x), y(v.y), z(v.z), w(1) {};

public:
  double operator[](int i) const {
    return (i == 0) ? x : (i == 1) ? y : (i == 2) ? z : w;
  };
  
  double& operator[](int i) {
    return (i == 0) ? x : (i == 1) ? y : (i == 2) ? z : w;
  };
  
public:
  point_t operator+(const vector_t& v) const {
    return point_t(x + v.x, y + v.y, z + v.z);
  };

  vector_t operator-(const point_t& p) const {
    return vector_t(x - p.x, y - p.y, z - p.z);
  };
};


vector_t::vector_t(const point_t& p) : x(p.x / p.w), y(p.y / p.w), z(p.z / p.w)
{
}


class matrix_t
{
public:
  matrix_t() {
    identitize();
  };
  
public:
  const double* operator[](int i) const {
    return matrix_[i];
  };

  double* operator[](int i) {
    return matrix_[i];
  };

public:
  void identitize() {
    memset(matrix_, 0, sizeof(matrix_));

    matrix_[0][0] =
      matrix_[1][1] =
      matrix_[2][2] =
      matrix_[3][3] = 1;
  };

public:
  matrix_t operator*(const matrix_t& m) const {
    matrix_t r;

    for (int i = 0; i < 4; i ++)
      for (int j = 0; j < 4; j ++)
        r[i][j] =
          matrix_[i][0] * m[0][j] +
          matrix_[i][1] * m[1][j] +
          matrix_[i][2] * m[2][j] +
          matrix_[i][3] * m[3][j];

    return r;
  };

  vector_t operator*(const vector_t& v) const {
    vector_t r;

    for (int i = 0; i < 3; i ++)
      r[i] =
        matrix_[i][0] * v.x +
        matrix_[i][1] * v.y +
        matrix_[i][2] * v.z +
        matrix_[i][3];
    
    return r;
  };

  point_t operator*(const point_t& p) const {
    point_t r;

    for (int i = 0; i < 4; i ++)
      r[i] =
        matrix_[i][0] * p.x +
        matrix_[i][1] * p.y +
        matrix_[i][2] * p.z +
        matrix_[i][3] * p.w;
    
    return r;
  };

public:
  static matrix_t translation(double x = 0, double y = 0, double z = 0) {
    matrix_t m;

    m[0][3] = x;
    m[1][3] = y;
    m[2][3] = z;

    return m;
  };

  static matrix_t rotation(double theta = 0) {
    matrix_t m;
    
    m[0][0] =  cos(theta);
    m[0][1] = -sin(theta);
    m[1][0] =  sin(theta);
    m[1][1] =  cos(theta);
    
    return m;
  };

  static matrix_t scale(double s = 1) {
    matrix_t m;

    m[0][0] = m[1][1] = m[2][2] = s;

    return m;
  };

private:
  double matrix_[4][4];
};


class point_circle_intersect : public std::unary_function<point_t, bool>
{
public:
  point_circle_intersect(const point_t& o, double r) : o_(o), rr_(r * r + EPSILON) {};

public:
  bool operator()(const point_t& p) const {
    vector_t v(p - o_);

    return v * v < rr_;
  };

private:
  point_t o_;
  
  double rr_;
};


class circle_circle_intersect
{
public:
  circle_circle_intersect(const point_t& p, double r) : p_(p), r_(r) {};

public:
  bool operator()(const point_t& p, double r) const {
    vector_t v(p - p_);

    r += r_;

    return v * v - r * r < EPSILON;
  };

private:
  point_t p_;

  double r_;
};


int n;

std::vector<point_t> points;
std::vector<double>  rs;


double normalize(double a)
{
  while (a < 0)
    a += 2 * M_PI;

  while (a > 2 * M_PI)
    a -= 2 * M_PI;

  return a;
}

double dfs(int i, double s, double e, int j)
{
  if (j == n)
    return rs[i] * (e - s);

  if (i == j) {
    return dfs(i, s, e, j + 1);
  }
  else if (circle_circle_intersect(points[i], rs[i])(points[j], rs[j])) {
    vector_t v1(points[i]), v2(points[j]);
    double   r1(rs    [i]), r2(rs    [j]);

    double a = -(v1.x - v2.x) / (v1.y - v2.y);
    double b = (v2 * v2 - v1 * v1 + r1 * r1 - r2 * r2) / 2 / (- v1.y + v2.y);

    double A = a * a + 1;
    double B = a * (b - v1.y) - v1.x;
    double C = b * (b - 2 * v1.y) + v1 * v1 - r1 * r1;

    double D = sqrt(B * B - A * C);
 
    point_t p1((- B + D) / A, a * (- B + D) / A + b);
    point_t p2((- B - D) / A, a * (- B - D) / A + b);

    double ss = normalize(atan2(p1.y - v1.y, p1.x - v1.x));
    double ee = normalize(atan2(p2.y - v1.y, p2.x - v1.x));

    if (ss > ee)
      std::swap(ss, ee);

    double m = (ss + ee) / 2;

    if (point_circle_intersect(v2, r2)(matrix_t::translation(v1.x, v1.y) *
                                       matrix_t::scale(r1)               *
                                       matrix_t::rotation(m)             *
                                       vector_t(1))) {
      ;
    }
    else {
      std::swap(ss, ee);

      ss = normalize(ss);
      ee = normalize(ee);
    }

    if (ss < ee) {
      if (ss <= s && e <= ee) {
        return 0;
      }
      else if (ee <= s || e <= ss) {
        return dfs(i, s, e, j + 1);
      }
      else if (s <= ss && ee <= e) {
        return dfs(i, s, ss, j + 1) + dfs(i, ee, e, j + 1);
      }
      else if (ss <= e) {
        return dfs(i, ss, e, j + 1);
      }
      else if (ee <= e) {
        return dfs(i, s, ee, j + 1);
      }
      assert(false);
    }
    else {
      return dfs(i, std::max(s, ee), std::min(e, ss), j + 1);
    }
  }
  else {
    return dfs(i, s, e, j + 1);
  }
}


int main(int argc, char** argv)
{
  while (1) {
    std::cin >> n;

    if (n == 0)
      break;
    
    points.resize(n);
    rs    .resize(n);
    
    for (int i = 0; i < n; i ++)
      std::cin >> points[i].x >> points[i].y >> rs[i];

    double a = 0;

    for (int i = 0; i < n; i ++)
      a += dfs(i, 0, 2 * M_PI, 0);

    std::cout << std::setprecision(14) << a << std::endl;
  }
  
  return 0;
}

Compiler output:
a.cc: In member function 'void matrix_t::identitize()':
a.cc:136:5: error: 'memset' was not declared in this scope
  136 |     memset(matrix_, 0, sizeof(matrix_));
      |     ^~~~~~
a.cc:7:1: note: 'memset' is defined in header '<cstring>'; this is probably fixable by adding '#include <cstring>'
    6 | #include <vector>
  +++ |+#include <cstring>
    7 | 
a.cc: At global scope:
a.cc:220:44: warning: 'template<class _Arg, class _Result> struct std::unary_function' is deprecated [-Wdeprecated-declarations]
  220 | class point_circle_intersect : public std::unary_function<point_t, bool>
      |                                            ^~~~~~~~~~~~~~
In file included from /usr/include/c++/14/string:49,
                 from /usr/include/c++/14/bits/locale_classes.h:40,
                 from /usr/include/c++/14/bits/ios_base.h:41,
                 from /usr/include/c++/14/iomanip:42,
                 from a.cc:4:
/usr/include/c++/14/bits/stl_function.h:117:12: note: declared here
  117 |     struct unary_function
      |            ^~~~~~~~~~~~~~