この記事は, 2024 夏のブログリレー 1 日目の記事です.

はじめに

https://atcoder.jp/contests/arc136/tasks/arc136_d を解いたときに新しく勉強したことがあったので, 忘れないように書いておきます.
※本記事の後半にはこの問題の解法を述べる部分があります.

ゼータ・メビウス変換について

まず, ゼータ・メビウス変換の基本的な知識をまとめます.

ゼータ変換

数列の添字をまとめた添字集合 $J$ を考えます. このときに, この集合で半順序関係 $⪯$ が成立して, かつ関数 $f(x)$ が任意の $x\in J$ で定義されている時, $f$ をゼータ変換した関数 $F$ を次のように定義します.

$$F(n)=\sum _{x⪯n}f(x)$$

半順序関係

ある集合 $S$ について, 次の 3 つの性質を満たしているとき, 関係 $⪯$ は 半順序関係であるといいます.

1. 反射律 : 任意の $a\in S$ について, $a⪯a$
2. 反対称律 : 任意の $a,b\in S$ について, $a⪯b\wedge b⪯a\leftrightarrow a=b$
3. 推移律 : 任意の $a,b,c\in S$ について, $a⪯b\wedge b⪯c⟹a⪯c$

例 : $\mathbb{N}$ 上で, 関係 $\le$ は半順序関係である.

ゼータ変換は累積和を取る

上の定義を見ると, 競技プログラミング的には $F$ は累積和そのものであるとわかります. 実際, $J=\{x|1\le x\le n(n\in \mathbb{N})\}$ において, 関係 $\le$ を考えると, このとき数列 $\{F(1),F(2),\cdots ,F(n)\}$ は数列 $\{f(1),f(2),\cdots ,f(n)\}$ の累積和をとったものであることが自然にわかると思います.

メビウス変換

ゼータ変換と同様の条件下で, $F$ をメビウス変換した関数$f$を次のように定義します.

$$f(n)=\sum _{x⪯n}\mu (x,n)F(x)$$

ここで, $\mu (x,n)$ はメビウス関数であり, 次の性質を満たします.

$$\sum _{m\le x\le n}\mu (x,n)=\delta (m,n)$$

$\delta (m,n)$ はクロネッカーのデルタ.

メビウス変換は差分を取る

直感的には, ゼータ変換の逆のことを行っていると捉えることができます. すると先のゼータ変換では累積和を取ったのに対し, このメビウス変換では差分を取っていると考えられます. 実際, $\mu (x,y)$ を $x=y$ のとき $1$ , $x=y-1$ のとき $-1$ , そうでないとき $0$ をとる関数としてメビウス変換を行えば, 先の累積和から元の数列を復元することができます.

ゼータ変換の高速化

今回は, 次のようなゼータ変換を高速化することを考えます.

ある有限集合 $U=\{1,2,\cdots ,n\}$ に対し, 写像 $f:2^U\to \mathbb{R}$ が定められているとする. この時, $f$ を次のようにゼータ変換した関数 $F$ を求める.

$$\forall S\subseteq U,F(S)=\sum _{T\subseteq S}f(T)$$

$O(4^n)$ 解

for (int S = 0; S < (1 << n); S++) {
    for (int T = 0; T < (1 << n); T++) {
        if ((S & T) == T) F[S] += f[T];
    }
}

一番簡単ですぐに思いつくのは, 愚直解だと思います. これなら確実に求められますが, 計算量は $(2^n{)}^2=4^n$ で少し重めです.

$O(3^n)$ 解

これを高速化するために, まずは $T$ が見る集合の数を減らしたいという気持ちになります. 今のところ $T$ は $2^U$ のすべての要素を調べていますが, 例えば $S=\{1,3,4\}$ のときに $T=U$ のケースは明らかに見る必要がないです. このような無駄を削ることで, $T$ が $S$ の部分集合のみを調べるようにすると, 高速化が出来ます.

for (int S = 0; S < (1 << n); S++) {
    for (int T = S;; T = (T - 1) & S) {
        F[S] += f[T];
    }
}

サイズ $k$ の $U$ の部分集合の数は ${}_n{C}_k$ 個であり, それぞれについて部分集合は $2^k$ 個あるので, 計算量は ${\sum }_{k=0}^n{}_n{C}_k{2}^k=(1+2{)}^n=3^n$ です.

$O(n{2}^n)$ 解 (高速ゼータ変換)

これをさらに高速化するために, dp のようなことをします. これによって計算量 $O(n{2}^n)$ で行うゼータ変換を高速ゼータ変換といいます.

$$dp(S,i):=i以下の要素がSの部分集合で,他はSと一致している時のfの総和$$

とします. 例えば, $dp(10011,2)=f(10000)+f(10001)+f(10010)+f(10011)$ のようになります. つまり, 下 $k$ 桁以外は $S$ と一致して, 下 $k$ 桁については任意のケタ $i$ について $i=1⟹i\in S$ を満たすようなものすべてに対する $f$ の総和を求めるということになります.

次の 4 式を見比べてみましょう.

• $dp(10011,2)=f(10000)+f(10001)+f(10010)+f(10011)$
• $dp(10001,2)=f(10000)+f(10001)$
• $dp(10011,1)=f(10010)+f(10011)$
• $dp(10001,1)=f(10000)+f(10001)$

すると, 次のことが分かります.

• $dp(10001,2)=dp(10001,1)$
• $dp(10011,2)=dp(10001,1)+dp(10011,1)$

すなわち, この dp は, $i\notin S$ のとき $dp(S,i)=dp(S,i-1),$ $i\in S$ のとき $dp(S,i)=dp(S,i-1)+dp(S\setminus \{i\},i-1)$ と遷移すればよさそうです. こうすると, $0\le i\le n$ であることと, $n$ 個の要素からなる集合の部分集合は $2^n$ 個であることから, 計算量 $O(n{2}^n)$ でゼータ変換を行うことができます.

for (int i = 0; i < (1 << n); i++) {
    dp[i][0] = f(i);
}

for (int i = 1; i <= n; i++) {
    for (int j = 0; j < (1 << n); j++) {
        if (j & (1 << (i - 1))) {
            dp[j][i] = dp[j][i - 1] + dp[j - (1 << (i - 1))][i - 1];
        }
        else dp[j][i] = dp[j][i - 1];
    }
}

ところで, この遷移を観察すると, 任意の $S$ について $dp(S,i)$ は必ず $dp(\ast ,i-1)$ の形のところから遷移しているので, dp 配列を 1 次元にして使い回せるとメモリを節約できて嬉しいです.
添字のループの順番に気を付けながら実装すると, 次のように dp 配列の次元を落とすことができます. 1 次元にできたので, 名前は dp ではなく上の問題設定に従って F になりました.

for (int i = 0; i < (1 << n); i++) {
    F[i] = f(i);
}

for (int i = 0; i < n; i++) {
    for (int j = 0; j < (1 << n); j++) {
        if (j & (1 << i)) {
            F[j] += F[j - (1 << i)];
        }
    }
}

このようにして, $O(n{2}^n)$ でゼータ変換を行うことができました.

例

これを用いて冒頭に挙げた問題を解いていきます.

問題概要

長さ $n$ の整数列 $A=\{A_1,A_2,\cdots ,A_n\}$ がある. このうち, $1\le i<j\le n$ のペア $(i,j)$ であって, $A_i+A_j$ の計算をする際に繰り上がりが生じないものの個数を求めよ.

$2\le n\le 10^6,0\le A_i\le 10^6(1\le i\le n)$

考察

この問題は 10 進数ですが, まずは 2 進数で同様の問題が与えられたときにどう解くかを考えてみます.
2 進数の整数がそれぞれ $a=(a_1{a}_2\cdots a_m{)}_{(2)},b=(b_1{b}_2\cdots b_m{)}_{(2)}$ と表されるとします. この時, $a+b$ で繰り上げが発生しないためには, $a_i=b_i=1$ を満たす $i$ が存在しなければよいです. 即ち, $a_i\&b_i=0$ を満たせばよいです.

ここで, 次のような dp を考えてみましょう.

$$dp(i):=iより任意のbitの値が小さいものの数$$

この dp を考えると, ある整数 $a$ と足し合わせたときに繰り上がりが発生しないものの個数は, dp のサイズを $N=2^d$ として $dp(N-a-1)$ となります.
遷移ですが, ある bit が立っていなかったら, その bit が立っている数へ遷移するという形になります.

for (int i = 0; i < n; i++) {
    for (int j = 0; j < (1 << n); j++) {
        if (!(j & (1 << i))) {
            dp[j + (1 << i)] += dp[j];
        }
    }
}

このように, 高速ゼータ変換と全く同じような考え方で解くことができるとわかります.

さて, 元のように 10 進数に戻して考えてみます.
10 進数の整数がそれぞれ $a=(a_1{a}_2\cdots a_m{)}_{(10)},b=(b_1{b}_2\cdots b_m{)}_{(10)}$ と表されるとすると, 繰り上がりが発生しない条件は, $a_i+b_i\ge 10$ を満たす $i$ が存在しないことです. よって, 各桁について, $a_i<9$ であれば $a_i+1<10$ より $a\to a^{\prime }=(a_1{a}_2\cdots (a_i+1)\cdots a_m)$ の遷移が存在します.
この遷移をすべての $0\le a\le 10^n-1$ に対して行う, という操作をすべてのケタに対して行えばよいです.

dp 配列の初期化については, $A_i=A_i+0$ というペアだと考え, $\forall i\in \{1,2,\cdots ,n\},dp(A_i)=1$ と初期化しておきます (他はすべて $0$ ) . ただ, $A_i$ が重複することがあるので, その際はその数だけ足します. 問題の制約より, $A_i$ は高々 6 桁なので, dp 配列のサイズは $10^6$ でよいでしょう.

vector<int>dp(1000000, 0);
for (int i = 0; i < n; i++) {
    dp[A[i]]++;
}

int digit = 1;
for (int i = 0; i < 6; i++) {
    for (int j = 0; j < 1000000; j++) {
        if ((j / digit) % 10 < 9) {
            dp[j + digit] += dp[j];
        }
    }
    digit *= 10;
}

最終的に出力する答えは, $\frac{1}{2}{\sum }_{i=1}^{n}dp(999999-A[i])$ となります. ただし, $A[i]+A[i]$ をするときに繰り上がりが発生しない場合, その分を答えから引く必要があるので注意してください.

int ans = 0;
for (int i = 0; i < n; i++) {
    ans += dp[999999 - A[i]];
}

// a+b で繰り上がりが発生しなければ true, 発生すれば false を返す関数 check(a, b) を作っておくと便利

for (int i = 0; i < n; i++) {
    if (check(A[i], A[i])) ans--;
}
cout << ans / 2 << endl;

あとは頑張ってコードを書くと, AC を得られます.
https://atcoder.jp/contests/arc136/submissions/56672429

最後に

調べたらもっと深いところまで色々出てきたのでそのへんも勉強していきたいですね.
明日は kashiwade さん, tatyam さん, alt--er さんです.