yukicoder No. 2530: Yellow Cards

問題解説 yukicoder

yukicoder.me

N \not\equiv 0 \pmod {998244353} のもとで、 O(\log N + \log K) 時間で解きます。

解法

[1] \Theta(N+K^2) 解法

dp[i][j] を、 i 回のイベントの直後に退場となった選手の人数がちょうど j 人であるような確率とします。

i 回のイベントの直後に退場となった選手の人数がちょうど j 人であるとき、今までの i 回のイエローカードの提示のうち 2j 回は選手を退場させるために使われているので、プレイエリアにいる選手のうちイエローカードを 1 回提示された選手の人数は i-2j 人となります。よって、以下のような漸化式が得られます。

  • dp[i+1][j] += \displaystyle\frac{N-(i-2j)}{N} dp[i][j]
  • dp[i+1][j+1] += \displaystyle\frac{i-2j}{N} dp[i][j]

また、初期値は dp[0][0] = 1、答えは \displaystyle\sum_{i=0}^{\infty} (N+i) \cdot dp[K][i] = N + \displaystyle\sum_{i=0}^\infty i \cdot dp[K][i] となります。(無限和の形で書いていますが、実際には加える要素数は O(K) 個となります。)

[2] \Theta(N+K) 解法

\displaystyle\sum_{j=0}^\infty dp[i][j] = 1 であることに注意します。 \displaystyle\sum_{j=0}^\infty j\cdot dp[i][j] の値に関心があるので、この値を dp2[i] とおきます。

dp の漸化式から、
\begin{align*} dp2[i+1] &= \sum_{j=0}^\infty \left(j\displaystyle\frac{N-(i-2j)}{N} + (j+1)\displaystyle\frac{i-2j}{N}\right)dp[i][j] \\ &= \sum_{j=0}^\infty \left(j+\displaystyle\frac{i-2j}{N}\right)dp[i][j] \\ &= \displaystyle\frac{N-2}{N}\sum_{j=0}^\infty j\cdot dp[i][j] + \frac{i}{N}\sum_{j=0}^\infty dp[i][j] \\ &= \displaystyle\frac{N-2}{N}dp2[i]+\frac{i}{N} \end{align*}

この漸化式に従い、 dp2[i+1] を dp2[i] から O(1) 時間で計算できます。

[3] \Theta(\log N + \log K) 解法

dp2 の漸化式を解きます。特殊解として dp2[i] = ai+b の形を仮定すると、

\begin{align*} a(i+1)+b &= \displaystyle\frac{N-2}{N}(ai+b) + \displaystyle\frac{i}{N} \\ \displaystyle\frac{2}{N}(ai+b)+a &= \displaystyle\frac{i}{N} \\ \displaystyle\frac{2a-1}{N}i + \left(a+\displaystyle\frac{2}{N}b\right) &= 0 \end{align*}

両辺を係数比較することで a = \displaystyle\frac{1}{2}, b = -\displaystyle\frac{N}{4} を得ます。よって、 dp2[i+1]-\left(\displaystyle\frac{1}{2}(i+1)-\displaystyle\frac{N}{4}\right) = \displaystyle\frac{N-2}{N}\left(dp2[i]-\left(\displaystyle\frac{1}{2}i-\displaystyle\frac{N}{4}\right)\right)

したがって、

dp2[K] = \displaystyle\frac{K}{2} - \displaystyle\frac{N}{4} + \displaystyle\frac{N}{4} \left(\displaystyle\frac{N-2}{N}\right)^K

なので、答えは

N + dp2[K] = \displaystyle\frac{K}{2} + \displaystyle\frac{3N}{4} + \displaystyle\frac{N}{4} \left(\displaystyle\frac{N-2}{N}\right)^K

となります。

これは \Theta(\log N + \log K) 時間で計算できます。

MMA Contest 016 I: regisys?

解法

[1] グリッドの問題への言い換え

各 i について、商品 i を購入できる一般人の人数を X _ i、商品 i を購入できる MMA 部員の人数を Y _ i とします。また、一般人の人数に 1 を加えた値を H、MMA 部員の人数に 1 を加えた値を W とおくと、以下のような問題に帰着されます。

H \times W のグリッドがある。 N 個の点があり、点 i はマス (X _ i, Y _ i) (0-indexed) にある。
i = 1, 2, \dots, H について以下の操作をそれぞれ行う。

  • x 座標が i 以上の点を 1 つ選んで取り除く。(以下操作 x _ i と呼ぶ)

i = 1, 2, \dots, W について以下の操作をそれぞれ行う。

  • y 座標が i 以上の点を 1 つ選んで取り除く。

残る点の最小個数を求めよ。(以下操作 y _ i と呼ぶ)

[2] Hall の結婚定理の利用

残る点の最小個数を求めるかわりに、最初にいくつかの点を消し、残った点の集合を U としたとき、操作により U に属する点をすべて取り除くことが可能であるために消す点の個数の最小値を求めることにします。

Hall の結婚定理により、 U に含まれる点をすべて取り除けることは、以下の条件と同値になります。

U の任意の部分集合 U' について、 U' に属する点を取り除くことが可能な操作の集合を \Gamma(U') とおくと、 |U'| \leq |\Gamma(U')| である。

U' が空の場合は明らかに成立するので、以下 U' は空でないとします。 U' に属する点の x 座標の最大値を X _ {\text{max}}、 y 座標の最大値を Y _ {\text{max}} とおくと、 \Gamma(U') = \{ x _ 1, x _ 2, \dots, x _ {X _ {\text{max}}}, y _ 1, y _ 2, \dots, y _ {Y _ {\text{max}}} \} となるので、 |\Gamma(U')| = X _ {\text{max}} + Y _ {\text{max}} です。よって |\Gamma(U')| は X _ {\text{max}} と Y _ {\text{max}} にしかよらないので、ある X _ {\text{max}} と Y _ {\text{max}} について x \leq X _ {\text{max}}, y \leq Y _ {\text{max}} に含まれる点をすべて含むような部分集合 U' のみ考えればよいです。

したがって、 0 \leq i \leq H, 0 \leq j \leq W に対し、 x \leq i, y \leq j に含まれる点の個数を A _ {ij} とおくと、点をすべて消せることは A _ {ij} \leq i + j がすべての i, j について成り立つことと同値になります。(点が存在しない x または y 座標についても考えていますが、点が存在するところまで i, j を小さくした条件から従います。) また、 B _ {ij} = A _ {ij} - i - j とおきます。

[3] B _ {ij} の性質

B _ {ij} は anti-monge になる、つまり、 i _ 1 < i _ 2, j _ 1 < j _ 2 に対し、 B _ {i _ 1 j _ 1} + B _ {i _ 2 j _ 2} \geq B _ {i _ 1 j _ 2} + B _ {i _ 2 j _ 1} が成り立つことを示します。

定義から、 B _ {i _ 1 j _ 1} + B _ {i _ 2 j _ 2} - B _ {i _ 1 j _ 2} - B _ {i _ 2 j _ 1} = A _ {i _ 1 j _ 1} + A _ {i _ 2 j _ 2} - A _ {i _ 1 j _ 2} - A _ {i _ 2 j _ 1} です。これは領域 i _ 1 \leq x < i _ 2, j _ 1 \leq y < j _ 2 に含まれる点の個数となるので 0 以上です。

消す点の個数の最小値は、実は B _ {ij} の最大値に等しくなります。これを帰納法により示します。

  • B _ {ij} の最大値は 0 以上である。 B _ {ij} の最大値が 0 の場合、点を消さなくても操作により点をすべて取り除ける。
  • B _ {ij} の最大値が m であるとする。 B _ {ij} = m となるような (i, j) のうち、 i の最小値を i'、 j の最小値を j' とおくと、 B _ {i'j'} = m であることを示す。ある i' \leq i, j' \leq j であって B _ {ij'} = B _ {i'j} = m であるようなものが存在する。 i = i' または j = j' の場合は自明なのでそうでないとすると、B の anti-monge 性より B _ {ij} + B _ {i'j'} \geq 2m で、 B の最大値は m より B _ {ij} = B _ {i'j'} = m となる。 B _ {i'j'} = m より、 x \leq i', y \leq j' の範囲に i' + j' + m > 0 個の点が存在する。このうち 1 つを選んで消すと、 B _ {ij} の最大値が m - 1 の場合に帰着する。また、それ以外の点を消すと B _ {ij} の最大値は m のままであるから、消す点の個数の最小値は m 個である。

[4] データ構造による高速化

B _ {ij} の最大値を直接求めると \Omega(HW) 時間がかかるため、これをデータ構造を用いて高速化します。区間 add・区間 max の遅延セグメント木を用いて平面走査すると O ( ( H+N) \log W + W) 時間で求めることができます。時間計算量は全体で O ( ( H+N+W) ( \log H + \log W)) 時間となります。

実装例

#include <bits/stdc++.h>
using namespace std;
const int INF = 10000000;
template <typename T>
struct lazy_segment_tree{
  int N;
  vector<T> ST;
  vector<T> lazy;
  lazy_segment_tree(vector<int> &A){
    int n = A.size();
    N = 1;
    while (N < n){
      N *= 2;
    }
    ST = vector<T>(N * 2 - 1, -INF);
    lazy = vector<T>(N * 2 - 1, 0);
    for (int i = 0; i < n; i++){
      ST[N - 1 + i] = A[i];
    }
    for (int i = N - 2; i >= 0; i--){
      ST[i] = max(ST[i * 2 + 1], ST[i * 2 + 2]);
    }
  }
  void eval(int i){
    if (i < N - 1){
      lazy[i * 2 + 1] += lazy[i];
      lazy[i * 2 + 2] += lazy[i];
    }
    ST[i] += lazy[i];
    lazy[i] = 0;
  }
  void range_add(int L, int R, T x, int i, int l, int r){
    eval(i);
    if (R <= l || r <= L){
      return;
    } else if (L <= l && r <= R){
      lazy[i] += x;
      eval(i);
    } else {
      int m = (l + r) / 2;
      range_add(L, R, x, i * 2 + 1, l, m);
      range_add(L, R, x, i * 2 + 2, m, r);
      ST[i] = max(ST[i * 2 + 1], ST[i * 2 + 2]);
    }
  }
  void range_add(int L, int R, T x){
    range_add(L, R, x, 0, 0, N);
  }
  T all(){
    eval(0);
    return ST[0];
  }
};
int main(){
  int N, M;
  cin >> N >> M;
  vector<int> A(N);
  for (int i = 0; i < N; i++){
    cin >> A[i];
  }
  vector<int> B(N);
  for (int i = 0; i < N; i++){
    cin >> B[i];
  }
  vector<int> T(M), C(M);
  for (int i = 0; i < M; i++){
    cin >> T[i] >> C[i];
  }
  vector<vector<int>> C2(2);
  for (int i = 0; i < M; i++){
    C2[T[i]].push_back(C[i]);
  }
  for (int i = 0; i < 2; i++){
    sort(C2[i].begin(), C2[i].end());
  }
  int H = C2[0].size() + 1;
  int W = C2[1].size() + 1;
  vector<int> x(N), y(N);
  for (int i = 0; i < N; i++){
    x[i] = C2[0].end() - lower_bound(C2[0].begin(), C2[0].end(), A[i]);
    y[i] = C2[1].end() - lower_bound(C2[1].begin(), C2[1].end(), B[i]);
  }
  vector<vector<int>> add(H);
  for (int i = 0; i < N; i++){
    add[x[i]].push_back(y[i]);
  }
  vector<int> a(W);
  for (int i = 0; i < W; i++){
    a[i] = -i;
  }
  lazy_segment_tree<int> ST(a);
  int ans = 0;
  for (int i = 0; i < H; i++){
    for (int j : add[i]){
      ST.range_add(j, W, 1);
    }
    ans = max(ans, ST.all());
    ST.range_add(0, W, -1);
  }
  cout << ans << endl;
}

AGC060 C: Large Heap

問題解説 AtCoder

atcoder.jp

解法

[1] 条件の言い換え

対称性により、不等号の向きをすべて反転させても求める確率は変わりません。後の計算がやりやすくなるので、不等号の向きをすべて反転させて考えます。

長さ 2^N-1 の実数列 X = (x_1, x_2, \dots, x_{2^N-1}) が以下の条件を満たすとき、 X は [0,1] 上の列であると呼ぶことにします。

  • 0 \leq x_i \leq 1 \ (1 \leq i \leq 2^N-1)

[0,1] 上の列 X = (x_1, x_2, \dots, x_{2^N-1}) がさらに以下の条件をすべて満たすとき、 X は [0,1] 上のヒープ的な列であると呼ぶことにします。

  • x_i > x_{2i} \ (1 \leq i \leq 2^{N-1}-1)
  • x_i > x_{2i+1} \ (1 \leq i \leq 2^{N-1}-1)

このとき、求める確率は [0,1] 上のヒープ的な列 X を一様ランダムに 1 つ選んだときに x_U > x_V となる確率に一致します。これは以下のようにしてわかります。

  • [0,1] 上の列 (x_1, x_2, \dots, x_{2^N-1}) がヒープ的であるか、また x_U > x_V であるかは、 x_i の間の大小関係のみに依存する。
  • [0,1] 上の列を一様ランダムに 1 つ選んだとき、等しい要素が存在する確率は 0 であるから無視できる。また、各要素が相異なるとき、要素間の大小関係としてありうるものは (2^N-1)! 通りあるが、それらは等確率で現れる。

さらに、以下のような根付き木 T を考えます。

  • 2^N-1 個の頂点を持つ。
  • 根は頂点 1 である。
  • 1 \leq i \leq 2^{N-1}-1 について、頂点 i は頂点 2i, 2i+1 の 2 つの子を持つ。

このとき、 [0,1] 上の列 (x_1, x_2, \dots, x_{2^N-1}) がヒープ的であることは、木 T の各頂点 i に x_i を書いたとき、頂点 u が頂点 v の親であるような木 T の頂点の組 (u,v) すべてについて頂点 u に書かれた値が頂点 v に書かれた値より大きいことと同値になります。

以上により、順列に関する問題を、木の各頂点に区間 [0,1] に含まれる実数を書く問題に言い換えることができました。

[0,1] 上の列 X = (x_1, x_2, \dots, x_{2^N-1}) を一様ランダムに選んだとき、 X がヒープ的である確率を p_1(N), X がヒープ的でかつ x_U > x_V である確率を p_2(N, A, B) とおくと、答えは \displaystyle\frac{p_2(N, A, B)}{p_1(N)} となります。よって、 p_1(N), p_2(N, A, B) をそれぞれ求めればよいです。

[2] p_1(N) の計算

p_1(N) を求めるときは、実数で考えるのではなく、順列のままで考えると求めやすいです。

以下の事実が典型として知られています。

N 頂点の根付き木がある。長さ N の順列 P_1, P_2, \dots, P_N を一様ランダムに 1 つ選んだとき、根以外の全ての頂点 v に対し、その親を p とすると P_p > P_v となる確率は、頂点 i の部分木の頂点数を \text{size}_i とおくと \displaystyle \frac{1}{\prod_{i=1}^N \text{size}_i} である。

木 T は深さ N-1 の完全二分木です。 i = 0, 1, \dots, N-1 に対し、深さ i の頂点は 2^i 個存在し、それらの頂点の部分木の頂点数はすべて 2^{N-i}-1 です。したがって、 \displaystyle p_1(N) = \frac{1}{\prod_{i=0}^{N-1} \left(2^{N-i}-1\right)^{2^i}} となります。

もしくは、最終的に用いる値が \displaystyle\frac{1}{p_1(N)} であることから、漸化式 \displaystyle\frac{1}{p_1(1)}=1, \frac{1}{p_1(N)}=\left(\frac{1}{p_1(N-1)}\right)^2\left(2^N-1\right) を用いると簡単に計算することができます。

[3] p_2(N, A, B) の計算

頂点 U, V に書く値については特別な制約があるため、最初に固定しておきます。具体的には、頂点 U には必ず a を、頂点 V には必ず b を書くことにします。ただし、 0 \leq b \lt a \leq 1 とします。

以下のような木 DP を考えます。

  • \text{dp}[v](a,b,x): 木 T の頂点 v の部分木の各頂点に区間 [0,1] に含まれる実数を独立に一様ランダムに書くが、頂点 U に数を書く場合には必ず a を書き、頂点 V に数を書く場合には必ず b を書くとき、頂点 v に書いた数が x 以下であり、かつ頂点 v の部分木で大小関係の条件が満たされている確率

ただし、頂点 v の部分木で大小関係の条件が満たされているとは、頂点 v の部分木に含まれる各頂点 w に対し頂点 w に書いた数が頂点 w の子に書いた数より大きくなることと定義します。また、 x は 0 \leq x \leq 1 の範囲のみ考えるとします。

DP の遷移を考えます。

頂点 v が葉で、 U または V である場合

頂点 v が葉なので、大小関係に関する制約はありません。

v = U のとき、頂点 v に書かれる数は必ず a なので、頂点 v に書かれる数が x 以下になる確率は x \lt a のとき 0, x \geq a のとき 1 となります。よって、 \text{dp}[v](a,b,x) = \begin{cases} 0 & (x \lt a) \\ 1 & (x \geq a) \end{cases} となります。

同様に、 v = V のときは \text{dp}[v](a,b,x) = \begin{cases} 0 & (x \lt b) \\ 1 & (x \geq b) \end{cases} となります。

頂点 v が葉で、 U でも V でもない場合

上と同じく、大小関係に関する制約はないので、頂点 v に書かれる数は区間 [0,1] から一様ランダムに選ばれます。よって、頂点 v に書かれる数が x 以下になる確率は x となるので、 \text{dp}[v](a,b,x)=x となります。

頂点 v が葉でなく、 U または V である場合

v = U のとき、頂点 v に書かれる数は必ず a です。よって、 x \lt a のとき明らかに \text{dp}[v](a,b,x)=0 となります。

x \geq a とします。頂点 v の部分木で大小関係の条件が満たされることは、以下の条件がともに成り立つことと同値になります。

  • 頂点 2v に書かれた数が a 以下であり、かつ頂点 2v の部分木で大小関係の条件が満たされている。
  • 頂点 2v+1 に書かれた数が a 以下であり、かつ頂点 2v+1 の部分木で大小関係の条件が満たされている。

よって、頂点 v の部分木で大小関係の条件が満たされる確率は \text{dp}[2v](a,b,a) \cdot \text{dp}[2v+1](a,b,a) です。

以上より、遷移は \text{dp}[v](a,b,x) = \begin{cases} 0 & (x \lt a) \\ \text{dp}[2v](a,b,a) \cdot \text{dp}[2v+1](a,b,a) & (x \geq a) \end{cases} となります。

同様に、 v = V の場合の遷移は \text{dp}[v](a,b,x) = \begin{cases} 0 & (x \lt b) \\ \text{dp}[2v](a,b,b) \cdot \text{dp}[2v+1](a,b,b) & (x \geq b) \end{cases} となります。

頂点 v が葉でなく、 U でも V でもない場合

頂点 v に書く数を t とします。頂点 v の部分木で大小関係の条件が満たされることは、以下の条件がともに成り立つことと同値になります。

  • 頂点 2v に書かれた数が t 以下であり、かつ頂点 2v の部分木で大小関係の条件が満たされている。
  • 頂点 2v+1 に書かれた数が t 以下であり、かつ頂点 2v+1 の部分木で大小関係の条件が満たされている。

これが成り立つ確率は \text{dp}[2v](a,b,t) \cdot \text{dp}[2v+1](a,b,t) です。

t は区間 [0,1] から一様ランダムに選ばれますが、 t \leq x である必要があるので、遷移は \displaystyle\text{dp}[v](a,b,x) = \int_0^x \text{dp}[2v](a,b,t) \cdot \text{dp}[2v+1](a,b,t) \ dt となります。

\text{dp}[1](a,b,1) は、木 T の頂点 U に a を、頂点 V に頂点 b を、他の各頂点に区間 [0,1] に含まれる実数を独立に一様ランダムに選んで書くとき、各辺の大小関係の条件が満たされる確率となります。よって、 \displaystyle p_2(N, A, B) = \int_0^1 \int_0^a \text{dp}[1](a,b,1) \ db da となります。

[4] 計算例

例として、サンプル 3 の結果を計算する過程を示します。

  • \displaystyle\frac{1}{p_1(4)} = 3^4 \cdot 7^2 \cdot 15 = 59535 である。
  • 頂点 8 は葉で頂点 U であるので、 \text{dp}[8](a,b,x) = \begin{cases} 1 & (x \geq a) \\ 0 & (x \lt a) \end{cases} である。以後、このような場合に \text{dp}[8](a,b,x) = 1 \ (x \geq a) と省略する。
  • 頂点 9 は葉で頂点 U でも頂点 V でもないので、 \text{dp}[9](a,b,x)=x である。頂点 10, 11, 12, 13, 14, 15 も同様である。
  • \displaystyle\text{dp}[4](a,b,x) = \int_0^x \text{dp}[8](a,b,t) \cdot \text{dp}[9](a,b,t) \ dt = \int_a^x 1 \cdot t \ dt = \frac{1}{2}\left(x^2-a^2\right) \ (x \geq a) である。
  • \displaystyle\text{dp}[5](a,b,x) = \int_0^x \text{dp}[10](a,b,t) \cdot \text{dp}[11](a,b,t) \ dt = \int_0^x t \cdot t \ dt = \frac{1}{3}x^3 である。頂点 6 も同様である。
  • \text{dp}[7](a,b,x) = \text{dp}[14](a,b,b) \cdot \text{dp}[15](a,b,b) = b^2 \ (x \geq b) である。
  • \displaystyle\text{dp}[2](a,b,x) = \int_0^x \text{dp}[4](a,b,t) \cdot \text{dp}[5](a,b,t) \ dt = \int_a^x \frac{1}{2}\left(t^2-a^2\right) \cdot \frac{1}{3}t^3 \ dt = \frac{1}{72}\left(2x^6 - 3x^4a^2+ a^6\right) \ (x \geq a) である。
  • \displaystyle\text{dp}[3](a,b,x) = \int_0^x \text{dp}[6](a,b,t) \cdot \text{dp}[7](a,b,t) \ dt = \int_b^x \frac{1}{3}t^3 \cdot b^2 \ dt = \frac{1}{12}\left(b^2x^4 - b^6\right) \ (x \geq b) である。
  • \displaystyle\text{dp}[1](a,b,1) = \int_a^1 \frac{1}{72}\left(2t^6 - 3t^4a^2 + a^6\right) \cdot \frac{1}{12}\left(b^2t^4 - b^6\right) \ dt である。
  • \displaystyle p_2(N,A,B) = \int_0^1 \int_0^a \int_a^1 \frac{1}{72}\left(2t^6 - 3t^4a^2 + a^6\right) \cdot \frac{1}{12}\left(b^2t^4 - b^6\right) \ dt db da である。これを計算すると \displaystyle\frac{89}{38102400} となる。
  • 答えは \displaystyle\frac{89}{38102400} \cdot 59535 = \frac{89}{640} である。

[5] 計算量の削減

頂点は 2^N-1 個あるため、すべての頂点について \text{dp} を計算することはできません。しかし、木 T には構造が同じ部分が多数存在することを利用すると、 \text{dp} を計算する必要のある頂点を O(N) 個に減らすことができます。各頂点の \text{dp} は a, b, x の多項式となりますが、それらの特徴と計算手順は以下のようになります。ただし、 P = 998244353 とします。

  • v を U の先祖でも V の先祖でもない頂点とすると、 \text{dp}[v] は頂点 v の高さのみに依存する x のみの単項式である。 \text{dp}[v] は O(N \log P) 時間の前計算のもと、 O(1) 時間で得ることができる。
  • v を U の先祖で V の先祖でない頂点とする。また、頂点 v, U の距離を d とおくと、 \text{dp}[v] は a, x の d+1 項の多項式であり、次数はすべて等しい。 d>0 のとき、 v の子で部分木に U を含むほうを w とおくと、 \text{dp}[v] は \text{dp}[w] から O(d \log P) 時間で計算できる。
  • v を V の先祖で U の先祖でない頂点とする。また、頂点 v, V の距離を d とおくと、 \text{dp}[v] は b, x の d+1 項の多項式であり、次数はすべて等しい。 d>0 のとき、 v の子で部分木に V を含むほうを w とおくと、 \text{dp}[v] は \text{dp}[w] から O(d \log P) 時間で計算できる。

これらを用いると、 \text{dp}[2], \text{dp}[3] を O(N^2 \log P) 時間で計算することができます。また、 \text{dp}[2] は a, x の O(N) 項の多項式、 \text{dp}[3] は b, x の O(N) 項の多項式なので、 \text{dp}[2](a,b,x) \cdot \text{dp}[3](a,b,x) は O(N^2) 項の多項式となります。

\displaystyle p_2(N,A,B) = \int_0^1 \int_0^a \int_a^1 \text{dp}[2](a,b,x) \cdot \text{dp}[3](a,b,x) \ dx db da なので、積分を項ごとに分解すると、 \displaystyle \int_0^1 \int_0^a \int_a^1 x^{e_1}a^{e_2}b^{e_3} \ dx db da を計算することに帰着されます。

0 \leq a \leq 1, 0 \leq b \leq a, a \leq x \leq 1 \Leftrightarrow 0 \leq x \leq 1, 0 \leq a \leq x, 0 \leq b \leq a なので、積分の順序を交換して
\displaystyle \int_0^1 \int_0^a \int_a^1 x^{e_1}a^{e_2}b^{e_3} \ dx db da = \int_0^1 \int_0^x \int_0^a x^{e_1}a^{e_2}b^{e_3} \ db da dx = \frac{1}{(e_3+1)(e_2+e_3+1)(e_1+e_2+e_3+1)} となります。これで p_2(N,A,B) を O(N^2 \log P) 時間で求めることができたので、答えも O(N^2 \log P) 時間で求めることができます。

実装

atcoder.jp

#include <bits/stdc++.h>
using namespace std;
const long long MOD = 998244353;
long long modpow(long long a, long long b){
  long long ans = 1;
  while (b > 0){
    if (b % 2 == 1){
      ans *= a;
      ans %= MOD;
    }
    a *= a;
    a %= MOD;
    b /= 2;
  }
  return ans;
}
long long modinv(long long a){
  a %= MOD;
  return modpow(a, MOD - 2);
}
vector<long long> POW2;
struct poly{
  vector<tuple<long long, long long, long long>> P;
  poly(): P(){
  }
  int size(){
    return P.size();
  }
  void multiply_monomial(poly f){
    int N = P.size();
    for (int i = 0; i < N; i++){
      get<0>(P[i]) *= get<0>(f.P[0]);
      get<0>(P[i]) %= MOD;
      get<1>(P[i]) += get<1>(f.P[0]);
      get<1>(P[i]) %= MOD;
      get<2>(P[i]) += get<2>(f.P[0]);
      get<2>(P[i]) %= MOD;
    }
  }
  void integrate(){
    int N = P.size();
    for (int i = 0; i < N; i++){
      get<1>(P[i])++;
      get<1>(P[i]) %= MOD;
      get<0>(P[i]) *= modinv(get<1>(P[i]));
      get<0>(P[i]) %= MOD;
    }
  }
  void integrate_a(){
    int N = P.size();
    long long sum = 0;
    long long deg = (get<1>(P[0]) + get<2>(P[0]) + 1) % MOD;
    for (int i = 0; i < N; i++){
      get<1>(P[i])++;
      get<1>(P[i]) %= MOD;
      get<0>(P[i]) *= modinv(get<1>(P[i]));
      get<0>(P[i]) %= MOD;
      sum += MOD - get<0>(P[i]);
      sum %= MOD;
    }
    P.push_back(make_tuple(sum, 0, deg));
  }
};
int main(){
  int N, A, B;
  cin >> N >> A >> B;
  POW2.resize(N + 1);
  POW2[0] = 1;
  for (int i = 0; i < N; i++){
    POW2[i + 1] = POW2[i] * 2 % MOD;
  }
  vector<poly> f(N);
  f[N - 1].P.push_back(make_tuple(1, 1, 0));
  for (int i = N - 2; i >= 0; i--){
    f[i] = f[i + 1];
    f[i].multiply_monomial(f[i]);
    f[i].integrate();
  }
  poly ga;
  if (A == N - 1){
    ga.P.push_back(make_tuple(1, 0, 0));
  } else {
    ga = f[A + 1];
    ga.multiply_monomial(ga);
    swap(get<1>(ga.P[0]), get<2>(ga.P[0]));
  }
  for (int i = A; i >= 2; i--){
    ga.multiply_monomial(f[i]);
    ga.integrate_a();
  }
  poly gb;
  if (B == N - 1){
    gb.P.push_back(make_tuple(1, 0, 0));
  } else {
    gb = f[B + 1];
    gb.multiply_monomial(gb);
    swap(get<1>(gb.P[0]), get<2>(gb.P[0]));
  }
  for (int i = B; i >= 2; i--){
    gb.multiply_monomial(f[i]);
    gb.integrate_a();
  }
  int N1 = ga.size();
  int N2 = gb.size();
  long long p2 = 0;
  vector<long long> T(N2);
  for (int i = 0; i < N2; i++){
    T[i] = modinv(get<2>(gb.P[i]) + 1);
  }
  for (int i = 0; i < N1; i++){
    for (int j = 0; j < N2; j++){
      long long c = get<0>(ga.P[i]) * get<0>(gb.P[j]) % MOD;
      long long x = (get<1>(ga.P[i]) + get<1>(gb.P[j])) % MOD;
      long long a = get<2>(ga.P[i]);
      long long b = get<2>(gb.P[j]);
      p2 += T[j] * modinv(a + b + 2) % MOD * modinv(x + a + b + 3) % MOD * c % MOD;
      p2 %= MOD;
    }
  }
  long long p1inv = 1;
  for (int i = 2; i <= N; i++){
    p1inv *= p1inv;
    p1inv %= MOD;
    p1inv *= POW2[i] - 1;
    p1inv %= MOD;
  }
  cout << p2 * p1inv % MOD << endl;
}

AOJ 2154: Shore Erosion

この記事は帰ってきた AOJ-ICPC Advent Calendar 2022 22 日目の記事です。

https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=2154

解法

与えられた多角形の内部を D、頂点を反時計回りに P _ 0, P _ 1, \dots, P _ {N-1} とします。また、 P _ N = P _ 0 であるとします。

辺 P _ i P _ {i+1} からのマンハッタン距離が R 以下である点の集合を D _ i とすると、島のある点が浸食される条件は、その点がある i に対し D _ i に含まれることです。よって、浸食された後の島は D \setminus (D _ 0 \cup D _ 1 \cup D _ {n-1}) となります。

領域 D _ i は四角形もしくは六角形になります。具体的には、点 A ( A _ x , A _ y ) , B ( B _ x , B _ y ) があるとき、線分 AB からのマンハッタン距離が R 以下である領域は以下のように得られます。

  • A _ x + A _ y \leq B _ x + B _ y と仮定しても一般性を失わないので、そのように仮定する。
  • A _ x - A _ y > B _ x - B _ y の場合、 6 点 (A _ x - R, A _ y), (B _ x - R, B _ y), (B _ x, B _ y + R), (B _ x + R, B _ y), (A _ x + R, A _ y), (A _ x, A _ y - R) を頂点とする六角形である。
  • そうでない場合、 6 点 (A _ x - R, A _ y), (A _ x, A _ y + R), (B _ x, B _ y + R), (B _ x + R, B _ y), (B _ x, B _ y - R), (A _ x, A _ y - R) を頂点とする六角形である。

Slab Decomposition で平面を D, D _ 0, D _ 1, \dots, D _ {N-1} の辺で分割し、それぞれの領域が浸食された後の島であるかどうかを調べます。ビット演算を用いて、それぞれの領域がどれに含まれるかを計算すると得られます。その後、島である領域と島でない領域が隣接する部分それぞれについて、長さを足し合わせればよいです。

実装

judge.u-aizu.ac.jp

#include <bits/stdc++.h>
using namespace std;
const double eps = 0.000001;
int sign(double x){
  if (x > eps){
    return 1;
  } else if (x < -eps){
    return -1;
  } else {
    return 0;
  }
}
struct point{
  double x, y;
  point(){
  }
  point(double x, double y): x(x), y(y){
  }
  point operator +(point P){
    return point(x + P.x, y + P.y);
  }
  point operator -(point P){
    return point(x - P.x, y - P.y);
  }
  point operator *(double k){
    return point(x * k, y * k);
  }
  point operator /(double k){
    return point(x / k, y / k);
  }
};
double abs(point P){
  return sqrt(P.x * P.x + P.y * P.y);
}
double dist(point P, point Q){
  return abs(Q - P);
}
point rotate(point P, double t){
  return point(P.x * cos(t) - P.y * sin(t), P.x * sin(t) + P.y * cos(t));
}
double cross(point P, point Q){
  return P.x * Q.y - P.y * Q.x;
}
struct line{
  point A, B;
  line(){
  }
  line(point A, point B): A(A), B(B){
  }
};
point vec(line L){
  return L.B - L.A;
}
bool is_parallel(line L1, line L2){
  return abs(cross(vec(L1), vec(L2))) < eps;
}
bool segment_line_intersect(line L1, line L2){
  return sign(cross(L1.A - L2.A, vec(L2))) * sign(cross(L1.B - L2.A, vec(L2))) <= 0;
}
bool segment_intersect(line L1, line L2){
  return segment_line_intersect(L1, L2) && segment_line_intersect(L2, L1);
}
point line_intersection(line L1, line L2){
  return L1.A + vec(L1) * cross(L2.A - L1.A, vec(L2)) / cross(vec(L1), vec(L2));
}
vector<double> unique_eps(vector<double> A){
  int N = A.size();
  sort(A.begin(), A.end());
  vector<double> A2;
  A2.push_back(A[0]);
  for (int i = 1; i < N; i++){
    if (A[i] - A[i - 1] > eps){
      A2.push_back(A[i]);
    }
  }
  return A2;
}
template <typename T>
struct slab_decomposition{
  function<T(T, T)> f;
  T e;
  vector<pair<line, T>> S;
  vector<double> X;
  int N;
  vector<vector<double>> Y;
  vector<int> P;
  vector<vector<tuple<int, int, T>>> E;
  slab_decomposition(function<T(T, T)> f, T e): f(f), e(e){
  }
  void add_segment(line s, T x){
    S.push_back(make_pair(s, x));
  }
  vector<vector<pair<T, int>>> build(){
    int M = S.size();
    for (int i = 0; i < M; i++){
      if (S[i].first.A.x > S[i].first.B.x){
        swap(S[i].first.A, S[i].first.B);
      }
    }
    for (int i = 0; i < M; i++){
      X.push_back(S[i].first.A.x);
      X.push_back(S[i].first.B.x);
    }
    for (int i = 0; i < M; i++){
      for (int j = i + 1; j < M; j++){
        if (!is_parallel(S[i].first, S[j].first) && segment_intersect(S[i].first, S[j].first)){
          X.push_back(line_intersection(S[i].first, S[j].first).x);
        }
      }
    }
    X = unique_eps(X);
    N = X.size();
    vector<int> L(M), R(M);
    for (int i = 0; i < M; i++){
      L[i] = lower_bound(X.begin(), X.end(), S[i].first.A.x - eps) - X.begin();
      R[i] = lower_bound(X.begin(), X.end(), S[i].first.B.x - eps) - X.begin();
    }
    vector<vector<double>> A(M, vector<double>(N));
    Y.resize(N);
    for (int i = 0; i < M; i++){
      for (int j = L[i]; j <= R[i]; j++){
        A[i][j] = line_intersection(S[i].first, line(point(X[j], 0), point(X[j], 1))).y;
        Y[j].push_back(A[i][j]);
      }
    }
    for (int i = 0; i < N; i++){
      Y[i] = unique_eps(Y[i]);
    }
    vector<vector<int>> B(M, vector<int>(N));
    for (int i = 0; i < M; i++){
      for (int j = L[i]; j <= R[i]; j++){
        B[i][j] = lower_bound(Y[j].begin(), Y[j].end(), A[i][j] - eps) - Y[j].begin();
      }
    }
    E.resize(N - 1);
    for (int i = 0; i < N - 1; i++){
      E[i].push_back(make_tuple(-1, -1, e));
    }
    for (int i = 0; i < M; i++){
      for (int j = L[i]; j < R[i]; j++){
        E[j].push_back(make_tuple(B[i][j], B[i][j + 1], S[i].second));
      }
    }
    for (int i = 0; i < N - 1; i++){
      E[i].push_back(make_tuple(Y[i].size(), Y[i + 1].size(), e));
    }
    for (int i = 0; i < N - 1; i++){
      sort(E[i].begin(), E[i].end());
      if (!E[i].empty()){
        vector<tuple<int, int, T>> E2;
        int cnt = E[i].size();
        E2.push_back(E[i][0]);
        for (int j = 1; j < cnt; j++){
          if (get<0>(E[i][j - 1]) == get<0>(E[i][j]) && get<1>(E[i][j - 1]) == get<1>(E[i][j])){
            get<2>(E2.back()) = f(get<2>(E2.back()), get<2>(E[i][j]));
          } else {
            E2.push_back(E[i][j]);
          }
        }
        E[i] = E2;
      }
    }
    P.resize(N + 2);
    P[0] = 0;
    P[1] = 1;
    for (int i = 0; i < N - 1; i++){
      P[i + 2] = P[i + 1] + E[i].size() - 1;
    }
    P[N + 1] = P[N] + 1;
    vector<vector<int>> idL(N), idR(N);
    for (int i = 0; i < N; i++){
      idL[i] = vector<int>(Y[i].size() + 1);
      idR[i] = vector<int>(Y[i].size() + 1);
    }
    idL[0] = vector<int>(Y[0].size() + 1, 0);
    for (int i = 0; i < N - 1; i++){
      int cnt = E[i].size();
      for (int j = 0; j < cnt - 1; j++){
        for (int k = get<0>(E[i][j]) + 1; k <= get<0>(E[i][j + 1]); k++){
          idR[i][k] = P[i + 1] + j;
        }
        for (int k = get<1>(E[i][j]) + 1; k <= get<1>(E[i][j + 1]); k++){
          idL[i + 1][k] = P[i + 1] + j;
        }
      }
    }
    idR[N - 1] = vector<int>(Y[N - 1].size() + 1, P[N]);
    int V = P[N + 1];
    vector<vector<pair<T, int>>> G(V);
    for (int i = 0; i <= N; i++){
      for (int j = P[i]; j < P[i + 1] - 1; j++){
        G[j].push_back(make_pair(get<2>(E[i - 1][j - P[i] + 1]), j + 1));
        G[j + 1].push_back(make_pair(get<2>(E[i - 1][j - P[i] + 1]), j));
      }
    }
    for (int i = 0; i < N; i++){
      for (int j = 0; j <= Y[i].size(); j++){
        G[idL[i][j]].push_back(make_pair(e, idR[i][j]));
        G[idR[i][j]].push_back(make_pair(e, idL[i][j]));
      }
    }
    return G;
  }
  int find(point p){
    int x = lower_bound(X.begin(), X.end(), p.x - eps) - X.begin();
    if (x == 0){
      return 0;
    } else if (x == N){
      return P[N];
    } else {
      int ans = P[x];
      for (int i = 1; i < E[x - 1].size() - 1; i++){
        point A(X[x - 1], Y[x - 1][get<0>(E[x - 1][i])]);
        point B(X[x], Y[x][get<1>(E[x - 1][i])]);
        if (cross(B - A, p - A) > 0){
          ans++;
        }
      }
      return ans;
    }
  }
  vector<double> area(){
    int V = P[N + 1];
    vector<double> ans(V, -1);
    for (int i = 0; i < N - 1; i++){
      for(int j = 1; j < E[i].size() - 2; j++){
        double l = Y[i][get<0>(E[i][j + 1])] - Y[i][get<0>(E[i][j])];
        double r = Y[i + 1][get<1>(E[i][j + 1])] - Y[i + 1][get<1>(E[i][j])];
        ans[P[i + 1] + j] = (l + r) * (X[i + 1] - X[i]) / 2;
      }
    }
    return ans;
  }
  vector<double> length(){
    int V = P[N + 1];
    vector<double> ans(V, -1);
    for (int i = 0; i < N - 1; i++){
      for(int j = 1; j < E[i].size() - 1; j++){
        double l = Y[i][get<0>(E[i][j])];
        double r = Y[i + 1][get<1>(E[i][j])];
        ans[P[i + 1] + j - 1] = hypot(X[i + 1] - X[i], r - l);
      }
    }
    return ans;
  }
};
int main(){
  cout << fixed << setprecision(20);
  while (true){
    int N, R;
    cin >> N >> R;
    if (N == 0 && R == 0){
      break;
    }
    vector<point> P(N + 1);
    for (int i = 0; i < N; i++){
      cin >> P[i].x >> P[i].y;
    }
    P[N] = P[0];
    vector<vector<point>> P2(N, vector<point>(6));
    for (int i = 0; i < N; i++){
      point A = P[i], B = P[i + 1];
      if (A.x + A.y > B.x + B.y){
        swap(A, B);
      }
      if (A.x - A.y > B.x - B.y){
        P2[i][0] = point(A.x - R, A.y);
        P2[i][1] = point(B.x - R, B.y);
        P2[i][2] = point(B.x, B.y + R);
        P2[i][3] = point(B.x + R, B.y);
        P2[i][4] = point(A.x + R, A.y);
        P2[i][5] = point(A.x, A.y - R);
      } else {
        P2[i][0] = point(A.x - R, A.y);
        P2[i][1] = point(A.x, A.y + R);
        P2[i][2] = point(B.x, B.y + R);
        P2[i][3] = point(B.x + R, B.y);
        P2[i][4] = point(B.x, B.y - R);
        P2[i][5] = point(A.x, A.y - R);
      }
    }
    for (int i = 0; i <= N; i++){
      P[i] = rotate(P[i], 1);
    }
    for (int i = 0; i < N; i++){
      for (int j = 0; j < 6; j++){
        P2[i][j] = rotate(P2[i][j], 1);
      }
    }
    slab_decomposition<__int128_t> D(bit_xor<__int128_t>(), 0);
    for (int i = 0; i < N; i++){
      D.add_segment(line(P[i], P[i + 1]), 1);
    }
    for (int i = 0; i < N; i++){
      for (int j = 0; j < 6; j++){
        D.add_segment(line(P2[i][j], P2[i][(j + 1) % 6]), (__int128_t) 1 << (i + 1));
      }
    }
    vector<vector<pair<__int128_t, int>>> E = D.build();
    int V = E.size();
    vector<__int128_t> A(V, -1);
    A[0] = 0;
    queue<int> Q;
    Q.push(0);
    while (!Q.empty()){
      int v = Q.front();
      Q.pop();
      for (pair<__int128_t, int> e : E[v]){
        int w = e.second;
        if (A[w] == -1){
          A[w] = A[v] ^ e.first;
          Q.push(w);
        }
      }
    }
    vector<double> L = D.length();
    double ans = 0;
    for (int i = 0; i < V - 1; i++){
      if (A[i] != 1 && A[i + 1] == 1 || A[i] == 1 && A[i + 1] != 1){
        ans += L[i];
      }
    }
    cout << ans << endl;
  }
}

AOJ 2394: Longest Lane

この記事は帰ってきた AOJ-ICPC Advent Calendar 2022 8 日目の記事です。

https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=2394

問題概要

N 頂点の凸と限らない多角形が与えられる。この多角形の内部または周上に完全に含まれる線分の長さの最大値を求めよ。

解法

多角形の頂点を反時計回りに P _ 0, P _ 1, \dots, P _ {N-1} とします。

実は、考える線分の候補は多角形の 2 頂点を通る線分のみとしてよいです。その 2 頂点を P _ i, P _ j とすると、以下を調べる必要があります。

  • 線分 P _ i P _ j が多角形の内部に含まれるか
  • P _ i から、 P _ j と逆方向にどこまで伸ばせるか
  • P _ j から、 P _ i と逆方向にどこまで伸ばせるか

多角形のそれぞれの辺と頂点について、線分を伸ばすことをどのように妨げるか調べます。辺については、その辺が線分 P _ i P _ j と端点以外で交わるかどうか調べ、交わる場合は交点と P _ i, P _ j の位置関係を調べればよいです。

頂点については、前後の辺も含めて位置関係を調べなくてはならないため、複雑になります。結局、以下の判定に帰着します。

  • 4 点 A, B, C, P が与えられる。半直線 AB を点 A を中心として AC まで反時計回りに回転させたとき、半直線が通過した領域に点 P が与えられるか判定する。

これは B - A と C - A の外積の符号で場合分けすることで簡単な判定に帰着します。以下の実装では、関数 point_in_angle でこの判定を行っています。

実装

judge.u-aizu.ac.jp

#include <bits/stdc++.h>
using namespace std;
const double eps = 0.000001;
const double INF = 100000000;
int sign(double x){
  if (x > eps){
    return 1;
  } else if (x < -eps){
    return -1;
  } else {
    return 0;
  }
};
struct point{
  double x, y;
  point(){
  }
  point(double x, double y): x(x), y(y){
  }
  point operator +(point P){
    return point(x + P.x, y + P.y);
  }
  point operator -(point P){
    return point(x - P.x, y - P.y);
  }
};
double abs(point P){
  return sqrt(P.x * P.x + P.y * P.y);
}
double dist(point P, point Q){
  return abs(Q - P);
}
double dot(point P, point Q){
  return P.x * Q.x + P.y * Q.y;
}
double cross(point P, point Q){
  return P.x * Q.y - P.y * Q.x;
}
bool point_in_angle(point A, point B, point C, point P){
  B = B - A;
  C = C - A;
  P = P - A;
  if (abs(cross(B, P)) < eps && dot(B, P) > 0){
    return true;
  }
  if (abs(cross(C, P)) < eps && dot(C, P) > 0){
    return true;
  }
  if (abs(cross(B, C)) < eps){
    if (dot(B, C) > 0){
      return false;
    } else {
      return cross(B, P) > 0;
    }
  }
  if (cross(B, C) > 0){
    return cross(B, P) > 0 && cross(P, C) > 0;
  } else {
    return cross(B, P) > 0 || cross(P, C) > 0;
  }
}
struct line{
  point A, B;
  line(point A, point B): A(A), B(B){
  }
};
point vec(line L){
  return L.B - L.A;
}
bool segment_line_intersect(line L1, line L2){
  return sign(cross(L2.A - L1.A, vec(L1))) * sign(cross(L2.B - L1.A, vec(L1))) < 0;
}
bool segment_intersect(line L1, line L2){
  return segment_line_intersect(L1, L2) && segment_line_intersect(L2, L1);
}
int main(){
  cout << fixed << setprecision(10);
  int T = 0;
  while (true){
    int N;
    cin >> N;
    if (N == 0){
      break;
    }
    vector<point> P(N * 2);
    for (int i = 0; i < N; i++){
      cin >> P[i].x >> P[i].y;
    }
    for (int i = 0; i < N; i++){
      P[N + i] = P[i];
    }
    double ans = 0;
    for (int i = 1; i <= N; i++){
      for (int j = i + 1; j < i + N - 1; j++){
        double mn = -INF, mx = INF;
        bool ok = true;
        if (!point_in_angle(P[i], P[i + 1], P[i - 1], P[j])){
          ok = false;
        }
        if (!point_in_angle(P[i], P[i + 1], P[i - 1], P[i] + (P[i] - P[j]))){
          mn = 0;
        }
        if (!point_in_angle(P[j], P[j + 1], P[j - 1], P[i])){
          ok = false;
        }
        if (!point_in_angle(P[j], P[j + 1], P[j - 1], P[j] + (P[j] - P[i]))){
          mx = 1;
        }
        for (int k = 0; k < N; k++){
          if (k != i - 1 && k != i % N && k != (j - 1) % N && k != j % N){
            if (segment_line_intersect(line(P[i], P[j]), line(P[k], P[k + 1]))){
              double r = cross(P[k] - P[i], P[k + 1] - P[k]) / cross(P[j] - P[i], P[k + 1] - P[k]);
              if (r < eps){
                mn = max(mn, r);
              } else if (r > 1 - eps){
                mx = min(mx, r);
              } else {
                ok = false;
              }
            }
          }
        }
        for (int k = 0; k < N; k++){
          if (k != i % N && k != j % N){
            if (abs(cross(P[k] - P[i], P[j] - P[i])) < eps){
              if (!point_in_angle(P[k], P[k + 1], P[k + N - 1], P[k] + P[i] - P[j]) || !point_in_angle(P[k], P[k + 1], P[k + N - 1], P[k] + P[j] - P[i])){
                if (dot(P[k] - P[i], P[j] - P[i]) < 0){
                  mn = max(mn, -abs(P[k] - P[i]) / abs(P[j] - P[i]));
                } else if (dot(P[k] - P[j], P[j] - P[i]) > 0){
                  mx = min(mx, abs(P[k] - P[i]) / abs(P[j] - P[i]));
                } else {
                  ok = false;
                }
              }
            }
          }
        }
        if (ok){
          ans = max(ans, (mx - mn) * dist(P[i], P[j]));
        }
      }
    }
    cout << "Case " << T + 1 << ": " << ans << endl;
    T++;
  }
}

AOJ 1394: Fair Chocolate-Cutting

この記事は帰ってきた AOJ-ICPC Advent Calendar 2022 1 日目の記事です。

https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=1394

問題概要

n 頂点の凸多角形が与えられる。この多角形の辺上に端点があり、多角形の面積を 2 等分する線分の長さとしてありうる最大値と最小値を求めよ。

解法

多角形の頂点を反時計回りに P _ 0, P _ 1, \dots, P _ {n-1}、面積を S とします。また、 P _ N = P _ 0 であるとします。

線分の端点を A, B (ただし A は辺 P _ i P _ {i+1} 上、 B は辺 P _ j P _ {j+1} 上の点) とします。多角形 P _ {i+1} P _ {i+2} \dots P _ j の面積を S _ 1、多角形 P _ {j+1} P _ {j+2} \dots P _ i の面積を S _ 2 とおくと、明らかに S _ 1, S _ 2 \leq \frac{S}{2} が必要です。

辺 P _ i P _ {i+1} と辺 P _ j P _ {j+1} が平行かどうかで場合分けをします。

辺 P _ i P _ {i+1} と辺 P _ j P _ {j+1} が平行でない場合

直線 P _ i P _ {i+1} と直線 P _ j P _ {j+1} の交点を Q とします。このとき、 Q が半直線 P _ i P _ {i+1} 上にある場合のみ考えるとしてよいです。(そうでない場合、 i と j を逆にして考える)

QP _ {i+1} = a _ 1, QA = a, QP _ i = a _ 2, QP _ j = b _ 1, QB = b, QP _ {j+1} = b _ 2, \angle P _ {i+1} Q P _ j = \theta とおくと、 a, b の取りうる値の範囲は a _ 1 \leq a \leq a _ 2, b _ 1 \leq b \leq b _ 2 となります。また、多角形 AP _ {i+1} P _ {i+2} \dots P _ j B の面積は S _ 1 + \frac{1}{2} (ab - a _ 1 b _ 1) \sin \theta となります。これが \frac{1}{2}S と等しくなればよいので、式変形して

ab = a _ 1 b _ 1 + \frac{S - 2S _ 1}{\sin \theta}

となります。右辺を p とおくと、 a = \frac{p}{b}, b = \frac{p}{a} です。よって、 a _ \text{min} = \max(a _ 1, \frac{p}{b _ 2}), a _ \text{max} = \min(a _ 2, \frac{p}{b _ 1}) とおくと、 a の取りうる値の範囲は a _ \text{min} \leq a \leq a _ \text{max} となります。

線分 AB の長さは余弦定理より \sqrt{a ^ 2 + b ^ 2 - 2ab \cos \theta} です。よって a ^ 2 + b ^ 2 - 2ab \cos \theta の最大値と最小値を求めればよいですが、これは (a + b) ^ 2 - 2p (1 + \cos\theta) に等しいので、 a + b = a + \frac{p}{a} の最大値と最小値を考えればよいです。

a > 0 において、 a + \frac{p}{a} は下に凸な関数で a = \sqrt{p} で極小値を取ります。従って、最大値を与える a の候補は a _ \text{min}, a _ \text{max}、最小値を与える a の候補は a _ \text{min}, a _ \text{max}, \sqrt{p} となります。

辺 P _ i P _ {i+1} と辺 P _ j P _ {j+1} が平行な場合

P _ {i+1} A = a, P _ {i+1} P _ i = a _ 0, P _ j B = b, P _ j P _ {j+1} = b _ 0 とおくと、と、 a, b の取りうる値の範囲は 0 \leq a \leq a _ 0, 0 \leq b \leq b _ 0 となります。また、直線 P _ i P _ {i +1}, P _ j P _ {j+1} の距離を d とおくと、多角形 AP _ {i+1} P _ {i+2} \dots P _ j B の面積は S _ 1 + \frac{1}{2} (a+b)d となります。これが \frac{1}{2}S と等しくなればよいので、式変形して

a+b = \frac{S - 2S _ 1}{d}

となります。右辺を s とおくと、 a = s - b, b = s - a です。よって、 a _ \text{min} = \max(0, s - b _ 0), a _ \text{max} = \min(a _ 0, s) とおくと、 a の取りうる値の範囲は a _ \text{min} \leq a \leq a _ \text{max} となります。

幾何的な考察から、最大値を与える a の候補は a _ \text{min}, a _ \text{max}、最小値の候補は a = a _ \text{min}, a _ \text{max} と AB \perp P _ i P _ {i+1} となるような a です。 AB \perp P _ i P _ {i+1} となるような a が存在するかは、内積の計算で判定することができます。

実装

judge.u-aizu.ac.jp

#include <bits/stdc++.h>
using namespace std;
const long double eps = 0.000001;
const long double INF = 100000;
struct point{
  long double x, y;
  point(){
  }
  point(long double x, long double y): x(x), y(y){
  }
  point operator +(point P){
    return point(x + P.x, y + P.y);
  }
  point operator -(point P){
    return point(x - P.x, y - P.y);
  }
  point operator *(long double k){
    return point(x * k, y * k);
  }
  point operator /(long double k){
    return point(x / k, y / k);
  }
};
long double abs(point P){
  return sqrt(P.x * P.x + P.y * P.y);
}
long double dist(point P, point Q){
  return abs(Q - P);
}
point unit(point P){
  return P / abs(P);
}
long double dot(point P, point Q){
  return P.x * Q.x + P.y * Q.y;
}
long double cross(point P, point Q){
  return P.x * Q.y - P.y * Q.x;
}
long double angle_cosine(point P, point Q){
  return dot(P, Q) / (abs(P) * abs(Q));
}
long double angle_sine(point P, point Q){
  return abs(cross(P, Q)) / (abs(P) * abs(Q));
}
long double angle_one_minus_cosine(point P, point Q){
  return 1 - angle_cosine(P, Q);
}
long double area(point A, point B, point C){
  return abs(cross(B - A, C - A)) / 2;
}
struct line{
  point A, B;
  line(){
  }
  line(point A, point B): A(A), B(B){
  }
};
point vec(line L){
  return L.B - L.A;
}
double point_line_distance(point P, line L){
  return abs(cross(P - L.A, vec(L))) / abs(vec(L));
}
point line_intersection(line L1, line L2){
  return L1.A + vec(L1) * cross(L2.A - L1.A, vec(L2)) / cross(vec(L1), vec(L2));
}
long double get_dist(long double a, long double b, long double one_minus_cost){
  return sqrt(max(pow(a - b, 2) + 2 * a * b * one_minus_cost, (long double) 0));
}
int main(){
  cout << fixed << setprecision(20);
  int n;
  cin >> n;
  vector<point> P(n + 1);
  for (int i = 0; i < n; i++){
    cin >> P[i].x >> P[i].y;
  }
  P[n] = P[0];
  long double S = 0;
  for (int i = 0; i < n; i++){
    S += cross(P[i], P[i + 1]) / 2;
  }
  long double mn = INF, mx = 0;
  for (int i = 0; i < n; i++){
    long double S1 = 0, S2 = S;
    for (int p = 1; p < n; p++){
      int j = (i + p) % n;
      S2 -= area(P[i], P[j], P[j + 1]);
      if (p >= 3){
        int j2 = (j - 1 + n) % n;
        S1 += area(P[j], P[j2], P[i + 1]);
      }
      if (S1 <= S / 2 + eps && S2 <= S / 2 + eps){
        long double C = cross(P[j + 1] - P[j], P[i] - P[i + 1]);
        if (C > eps){
          long double sint = angle_sine(P[j + 1] - P[j], P[i] - P[i + 1]);
          long double one_minus_cost = angle_one_minus_cosine(P[j + 1] - P[j], P[i] - P[i + 1]);
          point A = line_intersection(line(P[i], P[i + 1]), line(P[j], P[j + 1]));
          long double a1 = dist(A, P[i + 1]);
          long double a2 = dist(A, P[i]);
          long double b1 = dist(A, P[j]);
          long double b2 = dist(A, P[j + 1]);
          long double pr = a1 * b1 + (S - S1 * 2) / sint;
          long double xmn = max(a1, pr / max(b2, eps));
          long double xmx = min(a2, pr / max(b1, eps));
          if (xmn <= xmx + eps){
            long double d1 = get_dist(xmn, pr / xmn, one_minus_cost);
            mn = min(mn, d1);
            mx = max(mx, d1);
            long double d2 = get_dist(xmx, pr / xmx, one_minus_cost);
            mn = min(mn, d2);
            mx = max(mx, d2);
            long double x = sqrt(pr);
            if (xmn - eps <= x && x <= xmx + eps){
              long double d3 = get_dist(x, x, one_minus_cost);
              mn = min(mn, d3);
            }
          }
        } else if (C > -eps){
          long double a = dist(P[i], P[i + 1]);
          long double b = dist(P[j], P[j + 1]);
          long double h = point_line_distance(P[i + 1], line(P[j], P[j + 1]));
          long double s = (S - S1 * 2) / h;
          if (-eps <= s && s <= a + b + eps){
            long double xmn = max((long double) 0, s - b);
            long double xmx = min(a, s);
            point P1 = P[i + 1] + unit(P[i] - P[i + 1]) * xmn;
            point Q1 = P[j] + unit(P[j + 1] - P[j]) * (s - xmn);
            long double d1 = dist(P1, Q1);
            mn = min(mn, d1);
            mx = max(mx, d1);
            point P2 = P[i + 1] + unit(P[i] - P[i + 1]) * xmx;
            point Q2 = P[j] + unit(P[j + 1] - P[j]) * (s - xmx);
            long double d2 = dist(P2, Q2);
            mn = min(mn, d2);
            mx = max(mx, d2);
            if (dot(P[i] - P1, Q1 - P1) > -eps && dot(P[i] - P2, Q2 - P2) < eps){
              mn = min(mn, h);
            }
          }
        }
      }
    }
  }
  cout << mn << endl;
  cout << mx << endl;
}

Pinely Round 1 (Div. 1 + Div. 2) D. Carry Bit 別解

問題解説 Codeforces

https://codeforces.com/contest/1761/problem/D

問題概要

非負整数 x, y に対し 2 進数で x+y を計算したときの繰り上がりの回数を f(x, y) とおく。 整数の組 a, b \in [ 0, 2 ^ n) であって、 f(a, b) = k であるようなものの個数を 10 ^ 9+7 で割った余りを求めよ。

解法

a _ {n,k}, b _ {n,k} を以下のように定義します。

  • 整数 a, b \in [ 0, 2 ^ n) であって、 f(a, b) = k で 2^{n-1} の位から 2 ^ n の位への繰り上がりがあるものの個数を a _ {n, k} とする。
  • 整数 a, b \in [ 0, 2 ^ n) であって、 f(a, b) = k で 2^{n-1} の位から 2 ^ n の位への繰り上がりがないものの個数を b _ {n, k} とする。

また、 F(x,y) = \sum _ {n, k} a _ {n, k} x ^ n y ^ k, G(x,y) = \sum _ {n, k} b _ {n, k} x ^ n y ^ k とします。

このとき、 F(x,y), G(x,y) について以下が成り立ちます。

  • F(x,y) = 3xF(x,y) + xG(x,y) + 1
  • G(x,y) = xyF(x,y) + 3xyG(x,y)

(1-3xy)G(x,y)=xyF(x,y) なので、上の式の両辺を (1-3xy) 倍しこれを代入すると

(1-3xy)F(x,y) = 3x(1-3xy)F(x,y) + x ^ 2yF(x,y) + (1-3xy)

F(x,y) = \frac{1-3xy}{1-3x-3xy+8x ^ 2y}

よって G(x,y) = \frac{xy}{1-3x-3xy+8x ^ 2y} なので、

F(x,y)+G(x,y) = \frac{1-2xy}{1-3x-3xy+8x ^ 2y} = (1-2xy) \sum _ {i=0} ^ {\infty} (3x+3xy-8x ^ 2y) ^ i となります。

求めたい答えは [ x ^ ny ^ k ] (F(x,y)+G(x,y)) = [ x ^ n y ^ k ] (1-2xy) \sum _ {i = 0} ^ {\infty} (3x+3xy-8x ^ 2y) ^ i なので、 n, k について [ x ^ n y ^ k ] \sum _ {i = 0} ^ {\infty} (3x+3xy-8x ^ 2y) ^ i を求めることを考えます。

多項定理を用いて展開して

[ x ^ n y ^ k ] \sum _ {i = 0} ^ {\infty} (3x+3xy-8x ^ 2y) ^ i = [ x ^ n y ^ k ] \sum _ {a,b,c} \frac{(a+b+c)!}{a!b!c!} (3x) ^ a (3xy) ^ b (-8x ^ 2y) ^ c となります。

次数を比較すると、 a+b+2c=n, b+c=k となるような a, b, c のみ考えればよいです。 c を固定すると a, b は一意に定まります。この a, b, c について \frac{(a+b+c)!}{a!b!c!} 3 ^ {a+b} (-8) ^ c を計算し、総和を求めればよいです。

実装

https://codeforces.com/contest/1761/submission/186767977