約数集合でのゼータ変換・メビウス変換的なやつと畳み込み

昨日の CodeChef の July Lunchtime 2018 に参加したときに思いついたのでメモ.

※何も見ずに書いてるからもしかしたら間違いだらけかもしれないので注意です.

競プロで使うゼータ変換(?)は, 各集合について, その上位または下位集合の和を求めるやつです.
式で書くと

g(S) = Σ f(T) (S ⊆ T)　または　g(S) = Σ f(T) (T ⊆ S)

みたいな感じ.
で, これは集合をビットで管理してるときにも使えるけど, 約数集合とかでも使えるというのが今回の話.

上位集合をまとめる方のゼータ変換を使って以下を行います.

g(x) = Σ f(y) (y は x の倍数)

これで各値の個数 f から, 各値の倍数の個数 g を求めることができます.
実装は単純で, 前処理で各値の約数を列挙しておいて, 小さい値から約数側に足していきます.
ちなみに大きい値からやると重複が起こるのでダメ.

for (int x = 2; x <= MAX; x++) {
    for (int y : divs[x]) { // divs は最初に列挙しておく
        if (y != x) {
            f[y] += f[x];
        }
    }
}

こっちはゼータ変換の逆変換なので, 全操作を逆順にやっていくだけです.
大きい値からやりましょう.

for (int x = MAX; x >= 2; x--) {
    for (int y : divs[x]) {
        if (y != x) {
            f[y] -= f[x];
        }
    }
}

普通のやつだと, 二つの列をゼータ変換してからかけ合わせてメビウス変換すると, 上位集合なら積集合の個数, 下位集合なら和集合の個数みたいなのが求まります.
くわしくはこの記事を参照.
色々な畳み込み - kazuma8128’s blog

それを今回のに適用すると, 積集合に相当するものは GCD になります. (下位集合でやれば LCM になりそう)
式で書くとこう.

h(z) = Σ f(x) * g(y) (gcd(x, y) = z)

これをすると, 二つの列の任意のペアの gcd を取るときの, 各 GCD の値に対して個数が求まります.
(ところでこういうのも畳み込みっていえるのかな？よく知らない)

計算量はゼータ変換, メビウス変換の部分がボトルネックで O(X^(3/2)) (Xは最大値) になりますが, 一回前処理しとけば X までの値の約数の個数の和程度しかかからないので結構早くなります.

これを利用すれば以下の問題が解けます.
ただしこの畳み込みをさらにまとめて計算しないといけないので少し難しいですが.

ビットの方の高速ゼータ変換みたいな感じでこれも高速化できたりするのかな.
まあこれでも前処理すれば十分早いけど.

(追記 9/27)
よく考えたら前処理も変換も計算量は O(NlogN) でした. 調和級数になってるので.