http://dl.acm.org/citation.cfm?id=1963487 の実装です

approximate_k_nn_graph.cc

#include <bits/stdc++.h>
using namespace std;
#define rep(i,n) for(int i=0;i<(n);++i)

using Pt = pair<float, float>;

map<pair<int, int>, float> memo;

float distance(vector<Pt>&ps, int i, int j) {
  if (i == j) return 0.0;
  if (i > j) return distance(ps, j, i);
  pair<int, int> key = { i, j };
  if (memo.count(key) > 0) return memo[key];
  auto&p = ps[i];
  auto&q = ps[j];
  float dx = p.first - q.first,
        dy = p.second - q.second;
  return memo[key] = dx * dx + dy * dy; 
}

using vi = vector<int>;
using vvi = vector<vi>;

vvi naiive(vector<Pt>&ps, int k) {
  const int n = ps.size();
  vvi ret(n);
  rep (i, n) {
    vector<pair<float, int>> ln;
    rep (j, n) {
      if (i == j) continue;
      ln.push_back({ distance(ps, i, j), j });
    }
    sort(begin(ln), end(ln));
    rep (j, k) {
      ret[i].push_back(ln[j].second);
    }
  }
  return ret;
}

vvi apprx(vector<Pt>&ps, int k) {
  int count_distance = 0;
  const int n = ps.size();
  using Item = pair<float, int>;
  vector<set<Item>> b(n);
  rep (i, n) {
    rep (j, k) {
      int l = i;
      while (i == l) l = rand() % n;
      float d = distance(ps, i, l); ++count_distance;
      b[i].insert({ d, l });
    }
  }
  rep (_, 13) { // (infinite-)loop
    vvi u(n);
    rep (i, n) {
      auto it = begin(b[i]);
      rep (j, k) {
        auto&p = *it; ++it;
        int l = p.second;
        u[i].push_back(l);
        u[l].push_back(i);
      }
    }
    int cx = 0;
    const int rho = 2;
    rep (i, n) {
      if (u[i].size() > k/rho) {
        random_shuffle(begin(u[i]), end(u[i]));
        u[i].resize(k/rho);
      }
      for (int j: u[i]) {
        for (int l: u[i]) {
          if (j == l) continue;
          float d = distance(ps, j, l); ++count_distance;
          Item itm1 = {d, l},
               itm2 = {d, j};
          if (b[j].count(itm1) == 0) b[j].insert(itm1); else ++cx;
          if (b[l].count(itm2) == 0) b[k].insert(itm2); else ++cx;
        }
      }
    }
    if (cx == 0) {
      cerr << "loop finished after " << _ << " loops" << endl;
      break;
    } else {
      cerr << "cx = " << cx << endl;
    }
  }
  vvi ret(n);
  rep (i, n) {
    auto it = begin(b[i]);
    rep (j, k) {
      auto&p = *it; ++it;
      ret[i].push_back(p.second);
    }
  }
  cerr << "#calc of distance = " << count_distance << endl;
  cerr << "scan rate = " << (float(count_distance) / float(n * (n - 1) / 2)) << endl;
  return ret;
}

int main() {

  vector<Pt> ps;
  {
    float x, y;
    while (cin >> x >> y) ps.push_back({ x, y });
  }
  const int n = ps.size();
  {
    const int k = 20;
    auto nng = naiive(ps, k);
    /*
    rep (i, n/10) {
      cout << "#" << i << ": ";
      rep (j, k) {
        cout << nng[i][j] << ' ';
      } cout << endl;
    }
    */

    auto anng = apprx(ps, k);

    /*
    rep (i, n/10) {
      cout << "#" << i << ": ";
      rep (j, k) {
        cout << anng[i][j] << ' ';
      } cout << endl;
    }
    */

    {
      float recall = 0.0;
      rep (i, n) {
        set<int> a;
        int c = 0;
        rep (j, k) a.insert(nng[i][j]);
        rep (j, k) if (a.count(anng[i][j])) ++c;
        recall += float(c) / float(k);
      }
      recall /= float(n);
      cerr << "recall = " << recall << endl;
    }
  }


  return 0;
}