SRM 712 Div1 600 AverageVarianceSubtree

問題
ノードに重みがつけられた木が与えられる。
木のすべての部分木について重みの分散を計算し、この分散の平均値を求めよ。

解法
DFSしてノード上でDPする頻出テクニックを使う問題だが、分散を計算するために一工夫必要。すべての部分木の分散を効率よく計算するためには（ルート, 部分木のサイズ）ごとに"重みの二乗和"と"重みの和の二乗"を保持すればよい。

例として、ノード(a, b)からなる木とノード(c, d)からなる木をマージする作業を考えてみる。
木(a, b)の二乗和と木(c, d)の二乗和がわかっていたとすると、マージした木の二乗和はそのままこの二つを足せばいいだけなので簡単。

問題は、和の二乗の方で、
(a + b + c + d)^2 = a^2 + b^2 + c^2 + d^2 + 2ab + 2ac + 2ad + 2bc + 2bd + 2cd
(a + b)^2 = a^2 + 2ab + b^2
(c + d)^2 = c^2 + 2cd + d^2
となるので、単純に足すわけにはいかない。

和の二乗 = 二乗の和 + 2 * 異なる項の積和

となっているので、異なる項の積和をうまく部分問題から計算できればよいことが分かる。

木(a, b)の異なる項の積和 = ab
木(c, d)の異なる項の積和 = cd

木(a, b, c, d)の異なる項の積和 = ab + ac + ad + bc + bd + cd = ab + cd + (a + b)(c + d)

となるので、部分問題の異なる項の積和と部分問題の和を使えば計算できることが分かる。

以上をふまえて、

• sum[v][i] = ノードvを根とするサイズiの木の重み和
• sum_sq[v][i] = ノードvを根とするサイズiの木の重み二乗和
• sum_prd[v][i] = ノードvを根とするサイズiの木の重みの異なる項の積和

を葉から根方向にDPさせればよい。木をマージするときは、マージする相手方のサイズの分だけ自身が繰り返し現れるので、その分を掛けるのを忘れないようにする。

あとこの問題は数値計算の精度が厳しいらしく、計算誤差を小さくするため以下の工夫が必要。

• あらかじめ重みを中心化しておく。
• long doubleを使う。

ソースコード

using namespace std;

#define ALL(x) (x).begin(), (x).end()
#define EACH(itr,c) for(__typeof((c).begin()) itr=(c).begin(); itr!=(c).end(); itr++)  
#define FOR(i,b,e) for (int i=(int)(b); i<(int)(e); i++)
#define MP(x,y) make_pair(x,y)
#define REP(i,n) for(int i=0; i<(int)(n); i++)

int n;
long long sz[55][55];
double long sum[55][55];
double long sum_sq[55][55];
double long sum_prd[55][55];
double long w[55];
vector<int> child[55];

void dfs(int v) {
  REP (i, n+1) {
    sz[v][i] = 0;
    sum[v][i] = sum_sq[v][i] = sum_prd[v][i] = 0.0;
  }

  sz[v][1] = 1;
  sum[v][1] = w[v];
  sum_sq[v][1] = w[v] * w[v];

  for (auto &c : child[v]) {
    dfs(c);
    for (int i = n; i > 1; i--) {
      for (int j = 1; j < i; j++) {
        int k = i - j;
        sz[v][i] += sz[v][j] * sz[c][k];
        sum[v][i] += sum[v][j] * sz[c][k] + sum[c][k] * sz[v][j];
        sum_sq[v][i] += sum_sq[v][j] * sz[c][k] + sum_sq[c][k] * sz[v][j];
        sum_prd[v][i] += sum_prd[v][j] * sz[c][k] + sum_prd[c][k] * sz[v][j] + sum[v][j] * sum[c][k];
      }
    }
  }
}

class AverageVarianceSubtree {
  public:
  double average(vector<int> p, vector<int> weight) {
    n = weight.size();
    double long avg = 0.0;
    REP (i, n) avg += weight[i];
    avg /= n;
    REP (i, n) w[i] = weight[i] - avg;
    REP (i, n) child[i].clear();
    REP (i, p.size()) child[p[i]].push_back(i+1);

    dfs(0);

    long long tot = 0;
    long double ret = 0.0;

    REP (v, n) {
      for (int i = 1; i <= n; i++) {
        tot += sz[v][i];
        ret += sum_sq[v][i] / i - (sum_sq[v][i] + 2 * sum_prd[v][i]) / i / i;
      }
    }
    
    return ret / tot;
  }
};

SRM 531 Div2 600 NoRepeatPlaylist

問題概要
スマホにN曲の歌が入っている。
これらの曲を組み合わせてP曲のプレイリスト を作りたい。ただしプレイリストは以下の条件を満たす必要がある。
1) すべての曲は最低でも1回はプレイされなければならない
2) 同じ曲をプレイする場合は、最低でも間に別の曲をM曲プレイしなければならない
プレイリストの構成方法のパターン数を求めよ。

解法
よく分からななかったのでとりあえずDPしてみた。

const long long MOD = 1e9 + 7;
int N, M, P;
long long dp[124][124][124];

long long solve(int pos, int cnt, int ng) {
  if (pos == P)
    return cnt == N;
  if (cnt > N)
    return 0;
  if (ng >= N)
    return 0;
  
  long long &ret = dp[pos][cnt][ng];
  if (ret == -1) {
    ret = 0;
    ret += (cnt-ng) * solve(pos+1, cnt, min(ng+1, M));
    ret += (N-cnt) * solve(pos+1, cnt+1, min(ng+1, M));
    ret %= MOD;
  }
  return ret;
}

class NoRepeatPlaylist {
 public:
  int numPlaylists(int N, int M, int P) {
    ::N = N;
    ::M = M;
    ::P = P;
    memset(dp, -1, sizeof(dp));
    return solve(0, 0, 0);
  }
};

より高速な解法
(1)の条件がなければ簡単に計算できることに注目して包除原理する。
N曲以下を使って(2)を満たすパターン数を数える。
これには、N-1曲しか使ってないパターンがあるのでそれを引く....
という感じに足して引いてを繰り返す。

const long long MOD = 1e9 + 7;
long long comb[128][128];

class NoRepeatPlaylist {
 public:
  int numPlaylists(int N, int M, int P) {
    comb[0][0] = 1;
    for (int i = 1; i < 128; i++) {
      comb[i][0] = comb[i][i] = 1;
      for (int j = 1; j < i; j++)
        comb[i][j] = (comb[i-1][j-1] + comb[i-1][j]) % MOD;
    }

    int sign = 1;
    long long ret = 0;
    for (int i = N; i > 0; i--) {
      long long pat = 1;
      for (int j = 0; j < P; j++)
        pat = pat * max(i-j, i-M) % MOD;
      ret += sign * pat * comb[N][i];
      ret = (ret % MOD + MOD) % MOD;
      sign = -sign;
    }
    return ret;
  }
};

SRM 507 Div1 500 CubePacking

問題
辺の長さが1の立方体がNs個、辺の長さがLの立方体がNb個ある。
これらの立方体を出来るだけ体積が小さい直方体に格納したい。
直方体の体積の最小値を求めよ。

解法
解の下限はNb*L*L*L+Nsである。
解の上限はNb*L*L*L+Ns+L*L-1である。（L*Lの底面を持つ直方体に立方体を詰め込む）
ということで、 解の候補をすべて列挙しても高々L*L程度のループをまわせばよい。

あとは、それぞれの解の候補に対して、その体積になるような直方体を全列挙すればよい。これは時間がかかりそうだが約数の個数は思ったほど多くないため十分に間に合う。

以下に整数の約数の個数がどの程度かを示す。
この結果はN周辺の数をサンプリングし、それらの約数の個数の平均値を取ったものである。

N	約数の個数
106	14.9
107	17.3
108	19.6
109	21.9

ざっくりと高々log(N)程度と覚えておくと良さそう（あくまでも平均値であることに注意）。

ソースコード

using namespace std;

#define ALL(x) (x).begin(), (x).end()
#define EACH(itr,c) for(__typeof((c).begin()) itr=(c).begin(); itr!=(c).end(); itr++)  
#define FOR(i,b,e) for (int i=(int)(b); i<(int)(e); i++)
#define MP(x,y) make_pair(x,y)
#define REP(i,n) for(int i=0; i<(int)(n); i++)

vector<long long> divisor(long long n) {
  vector<long long> ret;
  for (long long i = 1; i * i <= n; i++) {
    if (n % i == 0) {
      ret.push_back(i);
      if (i != n / i)
        ret.push_back(n / i);
    }
  }
  return ret;
}

class CubePacking {
  public:
  int getMinimumVolume(int Ns, int Nb, int L) {
    int s = L*L*L*Nb+Ns;
    for (int x = s;; x++) {
      auto ds = divisor(x);
      REP (i, ds.size()) REP (j, i+1) {
        int a = ds[i];
        int b = ds[j];
        int c = x / a / b;
        if (a * b * c != x)
          continue;
        if ((a/L) * (b/L) * (c/L) >= Nb)
          return x;
      }
    }
    return -1;
  }
};

SRM 502 Div1 500 TheProgrammingContestDivOne

問題概要
プロコンでは問題がn問出題される。
それぞれの問題に対して、

• point[i] 問題iを解いたときにもらえる得点
• minus[i] 問題iの単位時間あたりの得点の減衰
• time[i] 問題iを解くのに必要な時間

が与えられる。問題を解くのに使える時間はTである。
最適な戦略を用いたときに得られる得点の最大値を求めよ。

解法
この問題の面白いところは、

• 問題を解く順番の最適化
• 解くべき問題集合の最適化

という2つの最適化を考えなければいけないところ。

もし解く順序が分かっていれば0-1ナップサック問題に帰着できるので、順序を考える。

i=0,1,2,..,n-1の順で解くことを考えてみる。

このとき、

minus[k+1] * time[k] < minus[k] * time[k+1]

を満たさなければ、kとk+1の順序は逆にしたほうがいいことがわかる。

つまり問題を解く最適な順序は、

time[i] / minus[i]

が小さい順ということが分かる。

よって上の順でソートして、0-1ナップサック問題を解けばよい。

ソースコード

#define ALL(x) (x).begin(), (x).end()
#define EACH(itr,c) for(__typeof((c).begin()) itr=(c).begin(); itr!=(c).end(); itr++)  
#define FOR(i,b,e) for (int i=(int)(b); i<(int)(e); i++)
#define MP(x,y) make_pair(x,y)
#define REP(i,n) for(int i=0; i<(int)(n); i++)

long long dp[100000+1];
class TheProgrammingContestDivOne {
public:
  int find(int T, vector<int> point, vector<int> minus, vector<int> time) {
    int n = point.size();

    // optimize order
    vector<pair<double, int> > vs;
    REP (i, n) vs.push_back(make_pair((double)time[i]/minus[i], i));
    sort(vs.begin(), vs.end());
    
    // optimize set
    memset(dp, 0, sizeof(dp));
    for (auto &pr: vs) {
      int i = pr.second;
      for (long long j = T; j >= time[i]; j--) {
        dp[j] = max(dp[j], dp[j-time[i]] + point[i] - minus[i] * j);
      }
    }
    return *max_element(dp, dp+T+1);
  }
};

SRM 657 Div2 1000 PolynomialRemainder

問題

a, b, cが与えられる。

a x^2 + b x + c = 0 (mod 1,000,000,000)

となるようなxを求めよ。

解法

10^9で割って0になるので、

2^9で割って0、5^9で割って0となるようなxを求めて、中国剰余定理をすればよい。

実装

using namespace std;

#define ALL(x) (x).begin(), (x).end()
#define EACH(itr,c) for(__typeof((c).begin()) itr=(c).begin(); itr!=(c).end(); itr++)  
#define FOR(i,b,e) for (int i=(int)(b); i<(int)(e); i++)
#define MP(x,y) make_pair(x,y)
#define REP(i,n) for(int i=0; i<(int)(n); i++)

long long extgcd(long long a, long long b, long long &x, long long &y) {
  long long d = a;
  if (b != 0) {
    d = extgcd(b, a%b, y, x);
    y -= a/b * x;
  } else {
    x = 1;
    y = 0;
  }
  return d;
}

pair<long long, long long> crm(long long a, long long n, long long b, long long m) {
  a %= n;
  b %= m;
  long long p, q;
  extgcd(n, m, p, q);
  long long mod = n * m;
  long long x = (a * q  * m) % mod + (b * p * n) % mod;
  x = (x % mod + mod) % mod;
  return make_pair(x, mod);
}

long long calc(long long a, long long b, long long c, long long x, long long mod) {
  long long ret = ((a % mod * x) % mod) * x % mod;
  ret = (ret + b * x) % mod;
  ret = (ret + c) % mod;
  return ret;
}

class PolynomialRemainder {
public:
  int findRoot(int a, int b, int c) {
    long long s, t;

    int n = 1;
    REP (i, 9) n *= 2;
    for (s = 0; s < n; s++) {
      if (calc(a,b,c,s,n) == 0)
        break;
    }
    if (s == n)
      return -1;

    int m = 1;
    REP (i, 9) m *= 5;
    for (t = 0; t < m; t++) {
      if (calc(a,b,c,t,m) == 0)
        break;
    }
    if (t == m)
      return -1;

    cout << s << " " << t << endl;
    auto ret = crm(s, n, t, m);
    return ret.first;
  }
};

SRM 653 Div2 250 CountryGroup

問題概要

椅子にn人の人が座っている。同じ国から来た人は必ず隣に座っているらしい。

左から順に、あなたの国からは何人の人がここにいますか？と聞いていく。

これに対する回答がn個与えられるので、矛盾しないかどうかと矛盾しない場合は何カ国の人がここにいるかを求めよ。

解法

stackを使うと綺麗に書ける。

using namespace std;

#define ALL(x) (x).begin(), (x).end()
#define EACH(itr,c) for(__typeof((c).begin()) itr=(c).begin(); itr!=(c).end(); itr++)  
#define FOR(i,b,e) for (int i=(int)(b); i<(int)(e); i++)
#define MP(x,y) make_pair(x,y)
#define REP(i,n) for(int i=0; i<(int)(n); i++)

class CountryGroup {
  public:
  int solve(vector<int> a) {
    int r = 0;
    while (a.size()) {
      int x = a.back();
      REP (i, x) {
        if (a.empty() || a.back() != x)
          return -1;
        a.pop_back();
      }
      ++r;
    }
    return r;
  }
};

SRM 571 Div2 1000 MagicMoleculeEasy

問題概要
グラフG(V, E)が与えられる。
各ノードには得点p[i]が割り振られている。
Vの中からk個のノードを選びたい。ただし、ノードv-w間にエッジが存在するとき、v, wのいずれかは必ず選ばなければならない。
選ばれたノードの得点の和の最大値を求めよ。

解法
あるエッジ(v, w)に注目すると、
vが選ばれなかった場合、wは必ず選ばなければならない。
wが選ばれなかった場合、vは必ず選ばなければならない。
という2パターンの分岐が発生する。どちらかを選ぶとkは1つ減る。

つまり全列挙を行っても高々2^k程度のパターンしかなく全列挙が可能。
ノードを選ぶ・選ばないではなく、エッジで結ばれるノードのどちらを選ぶかで状態を考えるという発想力が必要な問題。

今回の問題のように、入力のサイズ（=今回はグラフのサイズ）とは関係のないパラメータkに対してのみ指数時間かかるアルゴリズムをFPT（Fixed Parameter Tractable）アルゴリズムと呼ぶらしい。

実装

using namespace std;

#define ALL(x) (x).begin(), (x).end()
#define EACH(itr,c) for(__typeof((c).begin()) itr=(c).begin(); itr!=(c).end(); itr++)  
#define FOR(i,b,e) for (int i=(int)(b); i<(int)(e); i++)
#define MP(x,y) make_pair(x,y)
#define REP(i,n) for(int i=0; i<(int)(n); i++)

int n;
bool conn[50][50];
int point[50];
vector<pair<int, int> > es;

int solve(int k, set<int> &s) {
  int p;
  for (p = 0; p < es.size(); p++) {
    if (!s.count(es[p].first) && !s.count(es[p].second))
      break;
  }

  if (p == es.size()) {
    int ret = 0;
    vector<int> ps;
    REP (i, n) {
      if (s.count(i))
        ret += point[i];
      else
        ps.push_back(point[i]);
    }
    sort(ps.rbegin(), ps.rend());
    REP (i, k) ret += ps[i];
    return ret;
  }

  if (es.size() - p > k * n)
    return -1;

  int ret = -1;
  int x = es[p].first;
  int y = es[p].second;

  s.insert(x);
  ret = max(ret, solve(k-1, s));
  s.erase(x);

  s.insert(y);
  ret = max(ret, solve(k-1, s));
  s.erase(y);
  
  return ret;
}

class MagicMoleculeEasy {
public:
  int maxMagicPower(vector<int> magicPower, vector<string> magicBond, int k) {
    n = magicPower.size();
    memset(conn, 0, sizeof(conn));
    REP (i, n) REP (j, n) conn[i][j] = magicBond[i][j] == 'Y';
    REP (i, n) point[i] = magicPower[i];
    es.clear();
    REP (i, n) REP (j, i) if (magicBond[i][j] == 'Y') es.push_back(make_pair(i, j));
    set<int> s;
    return solve(k, s);
  }
};

SRM 652 Div2 1000 NoRightTurnDiv2

問題概要
二次元座標上にN個の点が与えられる。
このN個の点すべてをある順番に従って訪れたい。以下の条件を満たすように移動するとき、どのような順で点を訪れればよいか求めよ。

• 点から点への移動するときは二点を結ぶ直線上を移動しなければならない
• 移動したパスが交差してはならない
• 右にターンすることはできない（パスは反時計周りにならないといけない）

解法
最も右側の点から始める。最も右側の点が複数個ある場合は、その中で最も下側にあるものを選ぶ。あとは反時計まわりになるように、外側の点からgreedyに訪問すればいい。反時計まわりの判定、最も外側の判定は外積を使って行う。

実装

using namespace std;

#define ALL(x) (x).begin(), (x).end()
#define EACH(itr,c) for(__typeof((c).begin()) itr=(c).begin(); itr!=(c).end(); itr++)  
#define FOR(i,b,e) for (int i=(int)(b); i<(int)(e); i++)
#define MP(x,y) make_pair(x,y)
#define REP(i,n) for(int i=0; i<(int)(n); i++)

int vis[50];

long long cross(long long x1, long long y1, long long x2, long long y2) {
  return x1 * y2 - x2 * y1;
}

class NoRightTurnDiv2 {
  public:
  vector<int> findPath(vector<int> x, vector<int> y) {
    memset(vis, 0, sizeof(vis));
    vector<int> ret;
    int n = x.size();
    int p2 = 0;
    FOR (i, 1, n) {
      if (make_pair(x[i], -y[i]) > make_pair(x[p2], -y[p2]))
        p2 = i;
    }
    x.push_back(x[p2]);
    y.push_back(y[p2]-1000);
    vis[p2] = true;
    ret.push_back(p2);
    int p1 = n;
    
    REP (i, n-1) {
      int p3 = -1;
      REP (j, n) {
        if (vis[j])
          continue;
        if (p3 == -1 ||
            (cross(x[p2]-x[p1], x[j]-x[p2], y[p2]-y[p1], y[j]-y[p2]) >= 0 &&
             cross(x[j]-x[p2], x[p3]-x[p2], y[j]-y[p2], y[p3]-y[p2]) >= 0)) {
          p3 = j;
        }
      }
      vis[p3] = true;
      ret.push_back(p3);
      p1 = p2;
      p2 = p3;
    }
    return ret;
  }
};

SRM 647 Div2 1000 BuildingTowers

問題概要
N個のビルがあり、ラベルが1からNまで振られている。
以下の条件を満たすとき、ビルNの高さの最大値を求めよ。

• ビル1の高さは0
• すべてのビルの高さは非負整数
• 隣り合うビルの高さの差はK以下
• ビルx[i]の高さはt[i]以下

解法
与えられたx[]のビルについて、最大の高さを求める。
このとき、前から見たときと、後ろから見たときに、3つ目の条件を満たすように最大の高さを計算しないといけない。
x[]のビルについて決まったら、x[i]とx[i+1]の間の区間に立てられるビルの最大の高さを求める。これをすべてのiについて行う。

後ろから見たときのチェックは見落としがちなので、注意が必要。

実装

using namespace std;

#define ALL(x) (x).begin(), (x).end()
#define EACH(itr,c) for(__typeof((c).begin()) itr=(c).begin(); itr!=(c).end(); itr++)  
#define FOR(i,b,e) for (int i=(int)(b); i<(int)(e); i++)
#define MP(x,y) make_pair(x,y)
#define REP(i,n) for(int i=0; i<(int)(n); i++)

const long long oo = 1LL<<60;

class BuildingTowers {
  public:
  long long maxHeight(int N, int K, vector<int> x_, vector<int> t_) {
    vector<long long> x(ALL(x_));
    vector<long long> t(ALL(t_));

    if (x.size() && x[0] == 1) t[0] = 0;
    else {
      x.insert(x.begin(), 1);
      t.insert(t.begin(), 0);
    }

    int m = x.size();
    for (int i = 1; i < m; i++)
      t[i] = min(t[i], t[i-1] + K * (x[i] - x[i-1]));
    for (int i = m-2; i >= 0; i--)
      t[i] = min(t[i], t[i+1] + K * (x[i+1] - x[i]));
    
    long long ret = 0;
    REP (i, m-1) {
      long long z = (t[i+1] - t[i] + K * (x[i] + x[i+1])) / (2 * K);
      if (x[i] <= z && z <= x[i+1])
        ret = max(ret, t[i] + K * (z - x[i]));
      ++z;
      if (x[i] <= z && z <= x[i+1])
        ret = max(ret, t[i+1] + K * (x[i+1] - z));
    }
    ret = max(ret, t[m-1] + K * (N - x[m-1]));
    return ret;
  }
};