naoya_t@hatenablog

いわゆるチラシノウラであります

2013 TCO Marathon Round 3 "CirclesSeparation"

TCO13MMR3 (6/5〜19) に参戦しました。MM参加は通算6回目。

[6/26更新] 最終順位36位。\TシャツGET☆-(ノ゚Д゚)八(゚Д゚ )ノイエーイ/
http://community.topcoder.com/longcontest/stats/?module=ViewOverview&rd=15683
f:id:n4_t:20130626174825p:plain

回想

目標は

  • "LIMITED EDITION 2013 TOPCODER OPEN T-SHIRT"獲得*1
  • Ratingを黄色にする。(現在青1238)




ここに挙げられてる中では最弱ですが頑張りました


結局参加したのは最初の4日と最後2日だけだったけれど面白かった。
上位との差は焼きなまし技術・剛体シミュ技術等の差かな…

前半4日:他の円と重ならない場所をsegment treeを使って調べる関数を書けて満足していた

最初は全円の重心を中心に重くて大きいのから順にぐるぐる並べてみてた。同じ関数でスタート地点を個々の円の中心にしたらVisualizer上では疎らな感じになったけど70万弱→90万超まで上がった。

最盛期のスクショ
f:id:n4_t:20130619182553p:plain

10件のExampleだけ見ててもoverfittingするだけなので、パラメータ探索のために今回さくらVPSでseed1〜1000のスコアを取るやつを書いてみた。ローカル環境より微妙に速かった。

ちなみにTCOのサーバはローカルの6割程度の速度。

パラメータいじり倒してスコアはすこし伸びたけれど、伸び悩んで一旦諦める。

後半2日:他の円を蹴飛ばすことを覚えた

放置していたら60位ぐらいまで落ちていた。
なんとかTシャツ圏内に返り咲けないものか。

他の円を押し出したり押し返されたりするなんちゃって剛体シミュを書いてみた。
80万台前半まで落ちていたのを80万台後半まで取り戻した。

ビジュアライザ改造はoverlapがあっても表示するようにした+個々の円のamount of workを表示してみた程度。
f:id:n4_t:20130620041033p:plain

Masa-Yさんyowaさんの辺りまで追いついたけれど追い越せなかった。




最終的に出来たのは、重さとか面積とかを加味した順番で貪欲に、空いている場所で最良っぽい場所に置くのと他のやつを物理的に跳ね飛ばして置く(なんちゃって剛体シミュ)ケースと比較して良い方に置いていくだけのコード。

終了時点での暫定41位。自分より上は殆ど赤か黄色。
f:id:n4_t:20130620041047p:plain

最後の6時間ぐらいで30位あたりから41位まで降下した。赤の人たちはわざと点数を落とすコードを書いて潜伏していると聞いていたのでこれは想定内。
気になるのはTシャツボーダーライン。R1〜R3のランキングを重複排除してソートした感じだとTシャツボーダーは44位タイあたり。(微妙)

・・・

長い長いSystem Testの後、最終スコアが確定。8000点ほど上がって 891,999.81 点で、41位から36位へ上昇。Tシャツは確定。
(TLE等のために)0点になっているテストケースが3つあるらしい。(恐らく暫定で4つ→最終で3つ?)
ローカルでの結果からの推測では900000点台に乗るつもりだったので、zerosのケースを潰したら恐らく30位あたりまで上がれたのだろう。(←時間計測にExampleを活用すると良さそう)
http://community.topcoder.com/longcontest/stats/?module=ViewOverview&rd=15683
f:id:n4_t:20130626174825p:plain

\TシャツGET☆-(ノ゚Д゚)八(゚Д゚ )ノイエーイ/


レーティングも青1238→青1424 (+186)。
f:id:n4_t:20130626174808p:plain

その他


面白い。あとでじっくり読む

次回へのメモ

提出コード

最終コード貼っときます。コメントはあえてそのままで。

#include <iostream>
#include <sstream>
#include <cstdio>
#include <cmath>
#include <cctype>
#include <algorithm>
#include <string>
#include <vector>
#include <deque>
#include <stack>
#include <queue>
#include <list>
#include <map>
#include <set>

// #include "cout.h"
//#define KICK 1

using namespace std;
#define all(c)  (c).begin(),(c).end()
#define tr(i,c)  for(typeof((c).begin()) i=(c).begin(); i!=(c).end(); i++)
#define rep(var,n)  for(int var=0;var<(n);var++)
#define repN(var,n)  for(int var=1;var<=(n);var++)
#define found(s,e)  ((s).find(e)!=(s).end())
#define cons(a,b) make_pair(a,b)
typedef pair<int,int> ii;

extern int param;

#define PI 3.14159265358979
#define eps 1.0e-9

#ifdef USE_LONG_DOUBLE
typedef long double Double;
# define Hypot hypotl
# define Sin sinl
# define Cos cosl
# define Atan2 atan2l
# define Abs fabsl
# define Sqrt sqrtl
# define Pow powl
# define Log2 log2l
# define Log10 log10l
# define Log logl
#else
typedef double Double;
# define Hypot hypot
# define Sin sin
# define Cos cos
# define Atan2 atan2
# define Abs fabs
# define Sqrt sqrt
# define Pow pow
# define Log2 log2
# define Log10 log10
# define Log log
#endif

clock_t T0;
#define TIME_LIMIT (5.555 * 0.65)
inline double T() { return (double)(clock() - T0)/CLOCKS_PER_SEC; }
inline double rT() { return 1.0 - T()/TIME_LIMIT; }

Double rndf() { return (Double)rand() / RAND_MAX; }

typedef Double position;
typedef int contents;
struct interval {
  position low, high;
  interval(position low, position high) :
      low(low), high(high) { }
};

struct segment_tree {
  vector<position> pos;
  struct node {
    vector<contents> values;
    position B, E, M;
    node *left, *right;
  } *root;
  template <class IN>
  segment_tree(IN begin, IN end) : pos(begin, end) {
    root = build_tree(0, pos.size()-1);
  }
  node *build_tree(int i, int j) {
    int m = (i+j)/2;
    node *t = new node;
    t->B = pos[i]; t->E = pos[j]; t->M = pos[m];
    t->left  = (i+1 < j ? build_tree(i, m) : NULL);
    t->right = (i+1 < j ? build_tree(m, j) : NULL);
    return t;
  }
  void insert(const interval& I, contents c) { insert(root, I, c); }
  void insert(node *v, const interval& I, contents c) {
    if (I.low <= v->B && v->E <= I.high) {
      v->values.push_back( c );
    } else {
      if (I.low  < v->M) insert(v->left , I, c);
      if (I.high > v->M) insert(v->right, I, c);
    }
  }
  template <class OUT>
  void query(position p, OUT out) { query(root, p, out); }
  template <class OUT>
  void query(node *t, position p, OUT out) {
    if (!t || p < t->B || p >= t->E) return;
    copy(t->values.begin(), t->values.end(), out);
    if (p < t->M) query(t->left, p, out);   // 半開区間の探索になる
    else          query(t->right, p, out);
  }
};



class CirclesSeparation {
 public:
  vector<double> minimumWork(vector<double> x, vector<double> y, vector<double> r, vector<double> m);
};

class Circle {
 public:
  int id;
  Double x, y, r, m;
  Double px, py;

 public:
  Circle(int _id, double _x, double _y, double _r, double _m) {
    id = _id;
    x = (Double)_x; y = (Double)_y; r = (Double)_r; m = (Double)_m;
    px = x; py = y;
  }
};

int N;
Double gx, gy;
Double totalS;

vector<Circle*> cs;
vector<Double> s, px, py;

vector<Double> ex, ey, er, em;
Double escore;
vector<int> ei;
vector<Double> kx, ky;

inline Double rad2deg(Double rad) {
  Double deg = rad / PI * 180.0;
  return (deg < 0) ? deg + 360 : deg;
}

Double try_(int L, Double th, Circle& ci) { // [0..L] とかぶらない場所を探す
  vector<pair<Double,Double> > segs;

  Double cx = ci.px, cy = ci.py;
  Double c = Cos(th), s = Sin(th);

  for (int i=0; i<=L; ++i) {
    // (gx, gy) + r(cos th, sin th)
    // with ex[i], ey[i], er[i]
    // (ex_i, ey_i) + r(cos th+90, sin th+90) = (gx, gy) + R(cos th, sin th)
    // (ex_i, ey_i) + r(-sin th, cos th) = (gx, gy) + R(cos th, sin th)
    // ex_i - r sin th = gx + R cos th
    // ey_i + r cos th = gy + R sin th
    // R cos th = -((ey_i - gy) - r sin th)
    // R sin th = ((ex_i - gx) - r cos th)

    Double dx = ex[i] - cx, dy = ey[i] - cy;
    Double R = dx*c + dy*s, r = dx*s - dy*c;

    Double r_ = Abs(r);
    Double rsum = er[i] + ci.r + eps;
    if (rsum > r_) {
      Double delta = Sqrt(rsum*rsum - r_*r_) + eps;
      Double lo = R - delta, hi = R + delta;
      segs.push_back(cons(lo, hi));
    }
  }

  int W = segs.size();
  if (W == 0) return 0;

  // 端点を取る (-> sort + uniq)
  set<position> spos;
  tr(it, segs) {
    spos.insert(it->first);
    spos.insert(it->second);
  }

  segment_tree st(all(spos)); // O(n log n)
  rep(i, W) {
    st.insert(interval(segs[i].first, segs[i].second), i);
  }

  Double t = 0; // eps;
  while (true) { // 平均何周回ってるんだろう
    vector<contents> out(L+1, -1);
    st.query((position)t, out.begin());
    Double kmax = t;

    if (out.size() > 0) {
      tr(it, out) {
        if (*it == -1) break;
        kmax = max(kmax, segs[*it].second);
      }
    }
    if (kmax == t) break; // 見つからなかった → t=0に置いていい
    t = kmax;
  }
  return t;
}

int TLIMIT, T360, TREM, Tpick_max, Tdiv;

int queueing_cnt;

vector<int> fall(int L, Double div) {
  priority_queue<pair<Double, int> > pq;

  rep(i, L) {
    int id = ei[i];
    Double dx = kx[i] - cs[id]->x,
           dy = ky[i] - cs[id]->y;
    Double d = Hypot(dx, dy);
    if (d == 0) continue;

    Double cost = cs[id]->m * d;
    Double e_x = dx / d, e_y = dy / d;

    Double g = cost / div;

    // 133.979
    // g /= totalS; // 133.678 772.41
    g /= Sqrt(totalS); // 133.876 771.274
    /// g *= totalS; // 134.036

    kx[i] -= g * e_x;
    ky[i] -= g * e_y;

    pq.push(cons(cs[id]->m, i));
  }

  vector<int> res;
  while (!pq.empty()) { res.push_back(pq.top().second); pq.pop(); }
  return res;
}

Double kick(int M, Double x, Double y) {
  // M個の先行者たち + M+1番目
  rep(i, M) {
    kx[i] = ex[i]; ky[i] = ey[i];
  }

  // fprintf(stderr, "kicking. ei[M] = %d, ci.id = %d¥n", ei[M], ci.id);
  kx[M] = x; ky[M] = y;

  queue<int> q;
  set<int> qs;
  q.push(M); qs.insert(M);

  //vector<int> ids = fall(M+1, 50);
  vector<int> ids = fall(M+1, 36);
  tr(it, ids) {
    int i = *it;
    if (!found(qs,i)) { q.push(i); qs.insert(i); }
  }

  while (!q.empty()) {
    int t = q.front(); q.pop();
    qs.erase(t);
    ++queueing_cnt;
#ifdef KICK
    fprintf(stderr, "(M=%d) Q%d: %d¥n", M, queueing_cnt, t);
#endif
    rep(i, M+1) {
      if (i == t) continue;
      
      Double dx = kx[i] - kx[t], dy = ky[i] - ky[t];
      Double d = Hypot(dx, dy);
#ifdef KICK
      if (d == 0) {
        fprintf(stderr, "assert fails: d = 0¥n");
      }
#endif
      Double rr = er[i] + er[t];
      if (d >= rr + eps) continue;

      // Double ux = dx / d, uy = dy / d; // u: unit vector
      // Double rest = rr+eps*5000 - d;
      Double rest = rr+eps*10000 - d;
      // Double ux = dx/d, uy = dy/d;
      Double vx = dx * (rest / d), vy = dy * (rest / d);

      //Double ipi = 0, ipt = 0;
      // Double ri = em[t]*(1.0 + di*1e-3), rt = em[i]*(1.0 + dt*1e-3);

      Double ri = em[t], rt = em[i];
      //ri = em[t] * Hypot(kx[t]+vx/5 - cs[t]->x, ky[t]+vy/5 - cs[t]->y);
      //rt = em[i] * Hypot(kx[i]-vx/5 - cs[i]->x, ky[i]-vy/5 - cs[i]->y);

      //fprintf(stderr, "ri: %g + %g, rt: %g + %g¥n", em[t], ipi, em[i], ipt);
      // ri = Pow(ri, 1.5); rt = Pow(rt, 1.5);
      
      Double li = rt / (ri+rt), lt = ri / (ri+rt);
#ifdef KICK
      fprintf(stderr, "  %d hits %d; rest = %g, rate = (%g, %g)¥n", t, i, rest, li, lt);
#endif
      // i側をriで、m側をrmで
      kx[i] += vx * li;
      ky[i] += vy * li;

      kx[t] -= vx * lt;
      ky[t] -= vy * lt;

      if (!found(qs,i)) { q.push(i); qs.insert(i); }
      if (!found(qs,t)) { q.push(t); qs.insert(t); }
    }
  }

  // kscore
  Double kscore = 0;
  rep(i, M+1) {
    int id = ei[i];
    Circle* ci = cs[id];
    kscore += ci->m * Hypot(ci->x - kx[i], ci->y - ky[i]);
  }
  return kscore;
}

Double put(int M, Circle& ci) {
  // M+1番目のcircleを置くよ
  if (M == 0) {
    ex[M] = ci.px;
    ey[M] = ci.py;
    er[M] = ci.r;
    em[M] = ci.m;
    ei[M] = ci.id;
    return 0;
  }
  // TLIMIT >= 164
  // maxで164
  Double dd = PI*2 / T360;

  bool has_zero = false;

  Double GD_RATIO = 0; //0.001;

  // 1段 (*T)
  priority_queue<pair<Double,pair<Double,Double> > > pq; // (cost, (u, th))
  for (Double dh=-PI; dh<PI; dh+=dd) {
    Double u = try_(M-1, dh, ci); // M * O(?)

    Double _x = ci.px + u*Cos(dh), _y = ci.py + u*Sin(dh);
    Double gd = Hypot(_x - gx, _y - gy);

    Double cost = ci.m * u + gd*GD_RATIO;
    pq.push(cons(-cost, cons(u, dh)));
    if (u == 0) { has_zero = true; break; }
  }

  Double min_u, min_th;

  if (has_zero) {
    min_u = pq.top().second.first;
    min_th = pq.top().second.second;
  } else {
    // dhの左右に、隣りのdhとの中点までをdivで均等割する
    Double d_ = dd/(Tdiv*2+1);

    // 2段 (pick_max * div * 2)
    priority_queue<pair<Double,pair<Double,Double> > > pq2; // (cost, (u, th))
    int c = 0;

    while (!pq.empty()) {
      pq2.push(pq.top());

      Double u = pq.top().second.first,
            dh = pq.top().second.second;
      pq.pop();

      if (u == 0) break;

      for (int k=-Tdiv; k<=Tdiv; ++k) {
        if (k == 0) continue;
        Double dh2 = dh + d_*k;
        Double u = try_(M-1, dh2, ci);

        Double _x = ci.px + u*Cos(dh2), _y = ci.py + u*Sin(dh2);
        Double gd = Hypot(_x - gx, _y - gy);

        Double cost = ci.m * u + gd*GD_RATIO;
        pq2.push(cons(-cost, cons(u, dh2)));
      }

      ++c;
      if (c == Tpick_max) break;
    }

    min_u = pq2.top().second.first;
    min_th = pq2.top().second.second;
  }

  ex[M] = ci.px + min_u*Cos(min_th);
  ey[M] = ci.py + min_u*Sin(min_th);
  er[M] = ci.r;
  em[M] = ci.m;
  ei[M] = ci.id;
  // fprintf(stderr, "putting. [%d] = %d (%g %g | %g)¥n", M, ci.id, (double)ex[M], (double)ey[M], (double)er[M]);

  return ci.m * min_u;
}

Double scoring(vector<Circle*>& cs, int L) {
  Double score = 0;
  rep(i, L) {
    //if (emk[i]) continue;
    int id = ei[i];
    Circle* ci = cs[id];
    score += ci->m * Hypot(ci->x - ex[i], ci->y - ey[i]);
  }
  return score;
}

vector<double> subWork(vector<double>& x, vector<double>& y, vector<double>& r, vector<double>& m, int t) {
  // fprintf(stderr, "PARAM = %d¥n", param);
  priority_queue<pair<Double, int> > pq;

  rep(i, N) {
    Double ox = px[i] - x[i], oy = py[i] - y[i];
    Double ofs = Hypot(ox,oy);
    if (N > 400) ofs = 0;
    
    Double dx = (Double)px[i] - gx, dy = (Double)py[i] - gy;
    Double dist = Hypot(dx,dy); // 重心からの距離

    Double mw = 0;
    rep(j, N) {
      if (j == i) continue;
      Double d = Hypot(x[j]-x[i], y[j]-y[i]);
      Double rr = r[j] + r[i];
      Double move = (rr - d)/rr;
      //if (move > 0) mw += move;
      if (move > 0) mw += 1;
    }

    Double score;

    // m小さいの, r大きいのは外に飛ばしたい
    // r小さいのを優先した方が満足数が増える?
    // 重心からの距離:どうだろう

    //score = -s[i]*80 + m[i]*2 - dist; // 145.01
    //score = -s[i]*120 + m[i]*2 - dist; //
    //score = -s[i]*200 + m[i]*2 - dist; //
    //score = -s[i]*80 + m[i]*1.5 - dist; // 146.182

    Double WT;
    if (N < 100) {
      WT = 50;
    }
    else if (N < 150) WT = 55;
    else if (N < 200) WT = 75;
    else if (N < 250) WT = 90;
    else if (N < 300) WT = 100;
    else if (N < 350) WT = 110;
    else if (N < 400) WT = 120;
    else if (N < 450) WT = 130;
    else if (N <= 500) WT = 140;
    rep(j, t) {
      WT *= 0.6;
    }

    Double MW_MAG = 0.0005;
    Double ns = totalS * N;
    if (ns < 100.0) MW_MAG = 0;
    else if (ns < 200.0) MW_MAG = 0.0005;
    else if (ns < 300.0) MW_MAG = 0.0033;
    else MW_MAG = 0.01;

    int M_MAG = 200;

    M_MAG = 3425 / sqrt(N);
    /*
    // score = -ofs*(1.0*3)/m[i] + -s[i]*WT + m[i]*(0.01*M_MAG) - dist - 0.0005*mw; // 143.714
    score = -ofs*(1.0*3)/m[i] + -s[i]*WT + m[i]*(0.01*M_MAG) - dist - 0.0005*mw; // 136.266
    score = -ofs*(1.0*2)/m[i] + -s[i]*WT + m[i]*(0.01*M_MAG) - dist - 0.0005*mw; // 136.166
    score = -ofs*(1.0*2)/m[i] + -s[i]*WT + m[i]*(0.01*M_MAG) - dist; // 136.186
    //score = -ofs*(1.0*2)/m[i] + -s[i]*WT*0 + m[i]*(0.01*M_MAG) - dist; // faible
    //score = -ofs*(1.0*2)/m[i] + -s[i]*WT*2 + m[i]*(0.01*M_MAG) - dist; // 137.49
    score = -ofs*(1.0*2)/m[i] + -s[i]*WT + m[i]*(0.015*M_MAG) - dist; // 138.564

    score = -ofs*1.5/m[i] + -s[i]*WT + m[i]*(0.01*M_MAG) - dist; // 134.465
    score = -ofs*2.0/m[i] + -s[i]*WT + m[i]*(0.01*M_MAG) - dist - 0.0005*mw; // 134.357
    score = -ofs*2.0/m[i] + -s[i]*WT + m[i]*(0.01*M_MAG) - dist; // 134.534
    */
    score = -ofs*2.0/m[i] + -s[i]*WT + m[i]*(0.01*M_MAG) - dist - 0.001*mw; // 133.55
    
    // score = m[i];

    pq.push(cons(score, i));
  }

  // cerr << "maxR = " << maxR << endl; 0.06〜0.14あたり

  ex.resize(N);
  ey.resize(N);
  er.resize(N);
  em.resize(N);
  ei.resize(N);
  escore = 0.0;

  kx.resize(N);
  ky.resize(N);

  Double lN = Log10(N);

  // TLIMIT = (int)(600000.0 / lN / lN / N);
  TLIMIT = (int)(400000.0 / lN / lN / N);
  if (TLIMIT < 96) TLIMIT = 96;

  //T360 = min(180, TLIMIT-56);
  Tpick_max = 4;
  Tdiv = 7;
  if (N > 250) Tdiv = 5;
  T360 = min(180, TLIMIT-2*Tdiv*Tpick_max);
  // fprintf(stderr, "TLIMIT=%d, T360=%d, div=%d, pick_max=%d¥n", TLIMIT, T360, Tdiv, Tpick_max);

  vector<int> cids; // 暫定スコアのいい順にidが入っている
  while (!pq.empty()) {
    cids.push_back(pq.top().second);
    pq.pop();
  }

  queueing_cnt = 0;

  rep(i, N) {
    int id = cids[i];
    escore += put(i, *cs[id]); // ex[i], ey[i] に位置

    Double cx = cs[id]->x, cy = cs[id]->y;
    // if (i >= 350) continue;

    if (rT() < 0.25) continue;

    priority_queue<pair<Double,pair<vector<Double>,vector<Double> > > > pq;

    Double kscore = kick(i, cx, cy);
    if (kscore < escore * 0.99999) {
      pq.push(cons(-kscore, cons(kx, ky)));
    }

    #define B_DIV 10

    // Double r_[] = { 0.0333, 0.0666, 0.1 };
    rep(b, B_DIV) { // 134.891 779.358
    // rep(b, 16) { / 134.715
    // rep(b, 8) { // 135.066
    // rep(b, 12) { // 134.337
      Double th = (1.0/B_DIV*b - 0.5)*PI*2;
      // Double th = (rndf() - 0.5)*PI*2;

      rep(c, 1) {
        // Double r = rndf(); // 135.68
        // Double r = 0.01 * rndf(); // 137.37
        // Double r = 0.05 * rndf(); // 134.695
        //Double r = 0.1 * rndf(); // 134.218
        // Double r = 0.01667 * (1+c); // r_[c]; // 0.0666;
        Double r = 0.0666;
        //Double r = 0.15 * rndf(); // 134.083
        // Double r = 0.2 * rndf(); // 134.305
        Double rx = r * Cos(th), ry = r * Sin(th);
        Double kscore = kick(i, cx+rx, cy+ry);
        if (kscore < escore * 0.99999) {
          pq.push(cons(-kscore, cons(kx, ky)));
        }
      }
    }

    while (!pq.empty()) {
      Double kscore = -pq.top().first;
      // if (kscore > escore * 0.99999) break;
      // if (kscore > escore * 0.99999) break;

      #ifdef KICK
      fprintf(stderr, "[i=%d] escore:%g kscore:%g¥n", i, escore, kscore);
      #endif
      rep(j, i+1) {
        ex[j] = pq.top().second.first[j];
        ey[j] = pq.top().second.second[j];
      }
      escore = kscore;
      pq.pop();
      break;
    }

    while (!pq.empty()) pq.pop();
  }


#ifdef KICK
  fprintf(stderr, "queueing count: %d / %d¥n", queueing_cnt, N);
#endif
    // ret
  vector<double> ans(N*2);

  {
    vector<Double> fx(N), fy(N);
    rep(i, N) {
      int id = ei[i];
      fx[id] = ex[i];
      fy[id] = ey[i];
    }

    rep(i, N) {
      ans[i*2] = (double)fx[i]; // - mdx;
      ans[i*2+1] = (double)fy[i]; // - mdy;
    }
  }

  return ans;
}



vector<double>
CirclesSeparation::minimumWork(vector<double> x, vector<double> y, vector<double> r, vector<double> m) {
  T0 = clock();

  N = x.size();

  gx = gy = 0;
  Double mA = 0;
  Double maxR = 0;
  totalS = 0;

  s.resize(N);
  cs.resize(N);
  px.resize(N);
  py.resize(N);

  // priority_queue<pair<Double, Circle*> > pq;
  rep(i, N) {
    r[i] += 1e-5; //eps;
    cs[i] = new Circle(i, x[i], y[i], r[i], m[i]);
    maxR = max(maxR, (Double)r[i]);

    gx += x[i] * m[i];
    gy += y[i] * m[i];
    mA += m[i];

    s[i] = (Double)r[i] * r[i];
    totalS += s[i];

    px[i] = (Double)x[i] + (rndf() - 0.5)*0.0001;
    py[i] = (Double)y[i] + (rndf() - 0.5)*0.0001;
  }
  gx /= mA; gy /= mA; // 重心

  Double G = 1, H = 1, E = 0.000001;

  vector<Double> fx(N,0), fy(N,0);
  // while (rep) {
  int REPEAT = 1000;
  if (N > 400) REPEAT = 0;
  rep(z, REPEAT) {
  // rep(z, 1000) {
    //
    fx.assign(N, 0);
    fy.assign(N, 0);

    for (int i=0; i<N-1; ++i) {
      for (int j=i+1; j<N; ++j) {
        // #i の位置: px[i], py[i]
        // #i の重心: cs[i].x, cs[i].y (不変)
        Double dx = px[j] - px[i], dy = py[j] - py[i];
        Double vw = r[i] + r[j];
        Double ru = Hypot(dx, dy);
        if (ru > vw) continue;

        if (ru < 0.0001) {
          ru = vw * 0.1;
          dx = rndf(); dy = rndf();
          Double r_ = Hypot(dx, dy);
          dx = dx/r_*ru, dy = dy/r_*ru;
        }
        Double ex = dx/ru, ey = dy/ru;
        Double r2 = ru*ru;
        // iにかかる分
        Double f = G / r2;
        if (f > vw) f = vw;

        Double fi = f + m[j] / (m[i] + m[j]);
        fx[i] += -ex * fi;
        fy[i] += -ey * fi;
        // jにかかる分
        Double fj = f / m[i] / (m[i] + m[j]);
        fx[j] += ex * fj;
        fy[j] += ey * fj;
      }
    }

    rep(i, N) {
      Double dx = px[i] - cs[i]->x, dy = py[i] - cs[i]->y;
      Double ru = Hypot(dx, dy);
      if (ru == 0) continue;
      Double ex = dx / ru, ey = dy / ru;
      Double r2 = ru*ru;
      Double f = H * m[i] * m[i] / r2;
      if (f > ru) f = ru;
      fx[i] += -ex * f;
      fy[i] += -ey * f;
    }

    rep(i, N) {
      // fprintf(stderr, "%d + (%g, %g)¥n", i, fx[i], fy[i]);
      px[i] += E * fx[i];
      py[i] += E * fy[i];
    }
  }

#if 0
  rep(i, N) {
    cs[i]->p(px[i], py[i]);
  }
#endif

  priority_queue<pair<double,vector<double> > > ans_pq;

  double r1 = rT();
  int t = 0;
  do {
    // Double t_0 = T();
    vector<double> ret = subWork(x, y, r, m, t++);

    double cost = 0;
    rep(i, N) {
      double xi = ret[i*2], yi = ret[i*2+1];
      cost += m[i] * hypot(xi-x[i], yi-y[i]);
    }
    ans_pq.push(cons(-cost, ret));

    // fprintf(stderr, "[N=%d][%g %g] TURN %d: SCORED %g¥n", N, t_0, T(), t, cost);
  } while (rT() > r1/2);

  // ending
  rep(i, N) delete cs[i];

  // fprintf(stderr, "FINAL %g¥n", -ans_pq.top().first);

  return ans_pq.top().second;
}

*1:Round1〜3での最高順位の上位100位タイまで。http://community.topcoder.com/tco13/marathon/marathon-rules/ 参照