2025-03-26

偶奇で場合分けするやつ

二重階乗

$(2n)!! = 2^n n!$

$(2n-1)!! = \displaystyle \frac{(2n)!}{2^n n!}$

Wallis 積分

$\displaystyle \int_0^{\frac{\pi}{2}} \sin^{2n}x \, dx = \frac{(2n)!}{2^{2n+1} (n!)^2}\pi$

$\displaystyle \int_0^{\frac{\pi}{2}} \sin^{2n-1}x \, dx = \frac{2^{2n-2}((n-1)!)^2}{(2n-1)!}$

$\Gamma$ 関数

$\displaystyle \Gamma(n) = (n-1)!$

$\displaystyle \Gamma\left(n+\frac{1}{2}\right) = \frac{(2n)!}{4^nn!}\sqrt{\pi}$

$\displaystyle \Gamma\left(-n+\frac{1}{2}\right) = \frac{(-4)^nn!}{(2n)!}\sqrt{\pi}$

半径 $r$ の球の体積

$V_n(r)$ : 半径 $r$ の $n$ 次元球の体積

$V_{2n}(r) = \displaystyle\frac{\pi^n}{n!} r^{2n}$

$V_{2n-1}(r) = \displaystyle \frac{2^{2n-1}(n-1)!\pi^{n-1}}{(2n-1)!}r^{2n-1}$

Gauss 積分

$\displaystyle \int_0^\infty x^{2n}e^{-ax^2}\,dx = \frac{(2n)!}{2^{2n+1} n! a^{n+\frac{1}{2} }} \sqrt{\pi}$

$\displaystyle \int_0^\infty x^{2n-1}e^{-ax^2}\,dx = \frac{(n-1)!}{2a^n}$

2024-11-23

誘導部分グラフの数え上げ問題

注意

グラフは全て単純無向グラフとする。
グラフ $G$ の頂点数を $N$ 、辺数を $M$ と表す。
頂点 $i$ の次数を $d _ i$ と表す。
グラフの隣接行列の成分を $A _ {i ,j}$ と表す。

概要

以下の問題を考える。

$N$ 頂点 $M$ 辺のグラフ $G$ が与えられる。 $G$ に部分グラフとして含まれる $C _ 3$ の個数と $G$ の補グラフ $\overline{G}$ に部分グラフとして含まれる $C _ 3$ の個数の和を求めよ。

この問題は、以下のような考察により $O(N+M)$ 時間で解くことができることが知られている。

相異なる頂点の組 $(i,j,k)$ であって、以下の条件を満たすものを考える。

$G$ において、 $i, j$ 間に辺が存在し、 $i, k$ 間に辺が存在しない。

このような組の個数は、 $i$ を固定することで簡単に数えられる。具体的には、個数は $\sum _ {i=1} ^ n d _ i (n - 1 - d _ i)$ と表される。
相異なる $3$ 頂点の集合 $\{i, j, k\}$ について、 $\{i, j, k\}$ による $G$ の誘導部分グラフまたは $\overline{G}$ の誘導部分グラフが $C _ 3$ となる場合、個数に $0$ の寄与があり、そうでない場合、個数に $2$ の寄与がある。
よって、求める個数は
$\binom{N}{3} - \frac{1}{2} \sum _ {i=1}^N d _ i (N - 1 - d _ i)$ である。

$G$ に含まれる $C _ 3$ の個数を数える通常のアルゴリズムは $O(M \sqrt{M})$ 時間であり、 $O(N+M)$ 時間で簡単に解けるアルゴリズムは知られていない。にもかかわらず、補グラフとの個数の和は簡単に計算できるところが興味深い。本記事では、より一般に、グラフに含まれる特定の誘導部分グラフの個数を数え上げる問題の構造を分析し、このような現象が生じる理由を説明する。

誘導部分グラフ数え上げ問題

グラフ $G, H$ に対し、グラフ $G$ の相異なる $|V(H)|$ 個の頂点の組 $(v _ 1, v _ 2, \dots, v _ {|V(H)|})$ であって、それらの頂点による $G$ の誘導部分グラフが $H$ に同型であるものの個数を $count _ H(G)$ と書く。

グラフを引数とする実数値関数 $f$ が以下の条件を満たすとき、 $f$ は 誘導部分グラフ数え上げ関数 であるという。

あるグラフ $H _ 1, H _ 2, \dots, H _ k$ と実数 $c _ 1, c _ 2, \dots, c _ k$ が存在し、任意のグラフ $G$ について $f(G) = \sum _ {i=1} ^ k c _ i \cdot count _ {H _ i} (G)$ と表される。

さらに、 $H _ i$ がすべて $n$ 頂点のグラフであるとき、 $f$ は $n$ 頂点誘導部分グラフ数え上げ関数 であるという。

すべての誘導部分グラフ数え上げ関数の集合を $V$ とおき、その部分集合である、 $n$ 頂点誘導部分グラフ数え上げ関数の集合を $V _ n$ とおく。 $V$ および $V _ n$ は、実数値関数に関する自然な加法とスカラー倍について $\mathbb{R}$ 上のベクトル空間をなす。 $n$ 頂点のグラフが同型を除いて $H _ 1, \dots, H _ {a _ n}$ の $a _ n$ 個であるとき、 $count _ {H _ i} (H _ j) = \delta _ {i,j}$ であるから、 $V _ n$ は $(count _ {H _ i}) _ i$ を基底とする $a _ n$ 次元ベクトル空間となる。

さらに、任意の関数 $F(N,M)$ について、 $O(F(N,M))$ 時間で計算できる誘導部分グラフ数え上げ関数の集合は $V$ の部分空間をなす。なぜならば、いくつかの関数 $f _ 1, \dots, f _ k$ がそれぞれ $O(F(N,M))$ 時間で計算できるならば、それらの線形結合も $O(F(N,M))$ 時間で計算できるからである。

$3$ 頂点の場合

$O(N+M)$ 時間で計算できる誘導部分グラフ数え上げ関数のなす $V _ 3$ の部分空間を $W$ とおく。 $W$ を特定することが目標である。

$3$ 頂点のグラフは同型を除いて $4$ つである。辺が $0, 1, 2, 3$ 本であるものをそれぞれ $H _ 1, H _ 2, H _ 3, H _ 4$ とおく。 $V _ 3$ の元を、基底 $(count _ {H _ i}) _ i$ による成分表示で $\mathbb{R} ^ 4$ の元と同一視する。

以下、 $\sum _ {i,j,k}$ と書いたら、暗黙に $i, j, k$ は相異なることを課すとする。

$(c_1, c_2, c_3, c_4) \in \mathbb{R}^4$ とし、 $f(n)$ を $f(n) = c _ {n+1} \, (n = 0,1,2,3)$ なる関数とすると、 $\sum _ {i, j, k} f(A _ {i,j} + A _ {j,k} + A _ {k,i})$ は $V _ 3$ に属し、その成分表示は $(c _ 1, c _ 2, c _ 3, c _ 4)$ である。

$f(n) = 1$ とおくと、 $\sum _ {i, j, k} f(A _ {i,j} + A _ {j,k} + A _ {k,i}) = N(N-1)(N-2)$ である。よって、 $N(N-1)(N-2) \in W$ で、その成分表示は $(1, 1, 1, 1)$ である。

$f(n) = n$ とおくと、 $i, j, k$ の対称性に注意し、 $\sum _ {i, j, k} f(A _ {i,j} + A _ {j,k} + A _ {k,i}) = 3 \sum _ {i, j, k} A _ {i,j} = 6M(N-2)$ である。よって、 $6M(N-2) \in W$ で、その成分表示は $(0, 1, 2, 3)$ である。

$f(n) = n^2$ とおくと、 $i, j, k$ の対称性に注意し、 $\sum _ {i, j, k} f(A _ {i,j} + A _ {j,k} + A _ {k,i}) = 3 \sum _ {i, j, k} A _ {i,j}^2 + 6 \sum _ {i, j, k} A _ {i,j} A _ {i,k} = 6M(N-2) + \sum _ {i=1} ^ N d _ i (d _ i - 1)$ である。よって、 $6M(N-2) + 6\sum _ {i=1} ^ N d _ i (d _ i - 1) \in W$ で、その成分表示は $(0, 1, 4, 9)$ である。 $6M(N-2)$ を引き $2$ で割ると、 $3\sum _ {i=1} ^ N d _ i (d _ i - 1) \in W$ で、その成分表示は $(0, 0, 1, 3)$ である。

$(1,1,1,1), (0,1,2,3), (0,0,1,3)$ は $1$ 次独立であるから、 $W$ は少なくとも $3$ 次元である。

グラフに含まれる $C _ 3$ を数え上げる問題、すなわち関数 $(1,0,0,0)$ を計算する問題が $O(N+M)$ 時間で解けないと仮定する。このとき、 $W \subsetneq V _ 3$ であるから、 $\dim W < \dim V _ 3 = 4$ である。よって、 $W$ は $(1,1,1,1), (0,1,2,3), (0,0,1,3)$ により張られる部分空間に一致する。これは、 $W = \{(c _ 1, c _ 2, c _ 3, c _ 4) \mid c _ 1 - 3 c _ 2 + 3 c _ 3 - c _ 4 = 0\}$ とも表される。

例 1

本記事冒頭の問題を解く方法を導出する。 $3$ 頂点に順序をつけて考え、最後に $6$ で割ればよい。 $(1, 0, 0, 1)$ は $1 - 3 \cdot 0 + 3 \cdot 0 - 1 = 0$ を満たすので、 $W$ の元である。よって、 $(1,0,0,1)$ は $O(N+M)$ で計算できる。

$(1,0,0,1)$ を $W$ の基底 $(1,1,1,1), (0,1,2,3), (0,0,1,3)$ の線形結合で表すと、 $(1,0,0,1) = (1,1,1,1) - (0,1,2,3) + (0,0,1,3)$ となる。よって、 $(1,0,0,1)$ を計算する計算式は $N(N-1)(N-2) - 6M(N-2) + 3 \sum _ {i=1} ^ N d _ i (d _ i - 1)$ である。これは確かに冒頭で求めた式の $6$ 倍と一致する。

例 2

以下の問題を考える。

相異なる頂点 $i, j, k$ であって、 $\{i, j, k\}$ により誘導される誘導部分グラフがちょうど $2$ 本の辺を持つものの個数を求めよ。

$(0, 0, 1, 0)$ を計算すればよい。 $1 \cdot 0 - 3 \cdot 0 + 3 \cdot 1 - 0 \cdot 0 = 3 \neq 0$ より、 $(0,0,1,0) \notin W$ である。よって、(上の仮定の下で) この問題は $O(N+M)$ 時間で解くことができず、 $C _ 3$ の数え上げアルゴリズムが必要である。

$(1,1,1,1), (0,1,2,3), (0,0,1,3), (1,0,0,0)$ は $V _ 3$ の基底をなす。 $(0,0,1,0)$ をこれらの線形結合で表すと、 $(0,0,1,0) = -3(1,1,1,1)+3(0,1,2,3)-2(0,0,1,3)+3(1,0,0,0)$ となる。よって、 $(0,0,1,0)$ を計算する計算式は $-3N(N-1)(N-2) + 3 \cdot 6M(N-2) - 2 \cdot 3\sum _ i d _ i (d _ i - 1) + 3count _ {C _ 3} G$ であり、 $O(M\sqrt{M})$ 時間で計算できる。

2023-11-04

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)$ 時間で計算できます。

2023-08-12

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;
}

2022-12-26

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;
}

2022-12-26

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;
  }
}

2022-12-08

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++;
  }
}

二重階乗

Wallis 積分

関数

半径 の球の体積

Gauss 積分

注意

概要

誘導部分グラフ数え上げ問題

頂点の場合

例 1

例 2

解法

[1] 解法

[2] 解法

[3] 解法

解法