LiveArchive 7227 Equilibrium State

問題
https://icpcarchive.ecs.baylor.edu/external/72/7230.pdf

XY平面に、M個のばねとN個の物体が釣り合って、K個の定点で固定されている。物体の重さはなく、バネのたるみもない。
ばね定数と、定点の座標が与えられるので、物体の座標を順に求めよ。
バネの両端点は、{定点, 物体}, あるいは {物体, 物体} となっている。
N(物体の数) <= 1000

解法
本番では通らない解法だと思うが、テストケースが弱いのかLiveArchiveでは通ってしまった。


各物体ごとに力の釣り合い式を立てる。未知数の個数=式の数を確認する。
その後はガウスジョルダンでACしてしまった。O(N^3)でテストケース有りなので普通に間に合うはずがないのでは?

力の釣り合いをX, Y成分ごとに1次元ベクトルで考えて方程式を立式すれば良い。
各々の物体について、繋がれたばねから掛かる力をフックの定理を適用して立式する。物体の座標をxとして他の物体や定点の座標を引くと、有向の伸びの長さが求まるので、w_1 * (x - x_1) + w_2 * (x - x_2) + ... = 0 という力の釣り合い式が求まる。
式を展開して、変数ごとにまとめる。定点と繋がれたばねから掛かる力は定数となるので右辺に持っていく。

N個の物体それぞれに立式し、未知数の数もN個であることを確認して、左辺を表すの行列に変数の係数を、右辺を表すベクトルに定数を埋める。
あとはX, Yそれぞれの成分について連立方程式を解く。

typedef double result_type;

double const EPS = 1e-8;  // for floating point number

typedef vector<result_type> vec;
typedef vector<vec> mat;

vec gauss_jordan(mat const& A, vec const& b) {
  int n = A.size();
  mat B(n, vec(n+1));
  rep(i, n) rep(j, n) B[i][j] = A[i][j];
  rep(i, n) B[i][n] = b[i];

  rep(i, n) {
    int pivot = i;
    REP(j, i, n) {
      if(abs(B[j][i]) > abs(B[pivot][i])) { pivot = j; }
    }
    swap(B[i], B[pivot]);

    if(abs(B[i][i]) < EPS) { return vec(); }

    REP(j, i+1, n+1) B[i][j] /= B[i][i];
    rep(j, n) {
      if(i != j) {
        REP(k, i+1, n+1) B[j][k] -= B[j][i] * B[i][k];
      }
    }
  }

  vec x(n);
  rep(i, n) x[i] = B[i][n];
  return x;
}

int K, M, N;
struct point { double x, y; };
vector<point> fpts;
vector<vector<pair<int, int> > > G;

int main() {

  int T; cin >> T;
  rep(_, T) {
    cout << "Test case number : " << _ + 1 << endl;
    cin >> K >> M >> N;
    fpts.clear(), fpts.resize(K);
    rep(i, K) {
      cin >> fpts[i].x >> fpts[i].y;
    }

    G.clear();
    G.resize(N);

    rep(i, M) {
      double w; cin >> w;
      int a, b; cin >> a >> b;
      if(a > 0) a--;
      if(b > 0) b--;
      if(a >= 0) {
        G[a].push_back(make_pair(b, w));
      }
      if(b >= 0) {
        G[b].push_back(make_pair(a, w));
      }
    }

    // i行目: 物体iについての各成分ごとの立式
    //  j列目: 物体jから受ける力
    //  B[i]: fixed pointsから受ける力
    mat Ax(N), Ay(N); rep(i, N) Ax[i].resize(N), Ay[i].resize(N);
    vec Bx(N), By(N);

    rep(target, N) {
      rep(i, G[target].size()) {
        int other = G[target][i].first;
        double w = G[target][i].second;
        if(other >= 0) {
          Ax[target][target] += w;
          Ax[target][other] -= w;
          Ay[target][target] += w;
          Ay[target][other] -= w;
        }
        else {
          Ax[target][target] += w;
          Bx[target] += w * fpts[-other-1].x;
          Ay[target][target] += w;
          By[target] += w * fpts[-other-1].y;
        }
      }
    }

    vec rx = gauss_jordan(Ax, Bx);
    vec ry = gauss_jordan(Ay, By);

    rep(i, N) {
      printf("%d %.5f %.5f\n", i + 1, rx[i], ry[i]);
    }

  }
  
  return 0;
}