noshi91のメモ

データ構造のある風景

高速ゼータ変換の約数版 O(N log(log(N)))

添え字 gcd での畳み込みで AGC038-C を解く - noshi91のメモ
高速化前の実装と、実用例などが書かれています

概要

kazuma8128.hatenablog.com  高速ゼータ変換は添え字を集合とみなして各添え字について部分集合の和を計算します。これと同じことを約数でも出来るという話が kazuma8128 さんのブログで言及されていて、面白いと思ったので実装と計算量について書きます。
 高速ゼータ変換と高速メビウス変換とその上位集合版があるように、これも約数と倍数と (逆元があるとき) その逆変換があります。これらを組み合わせることで  \gcd  \rm lcm で畳み込みみたいなのが出来ます。詳しくはリンク先の記事を参照してください。

普通に列挙  O ( N \log ( N ) )

 各添え字について、その約数を普通に全部足します。この計算量は一見すると非自明ですが、 k を約数に持つ数は  N / k 個なので調和級数 O ( N \log ( N ) ) となります。
 実装はエラトステネスの篩の素数以外も全部見る版みたいにするとやりやすいです。in_place で実装する場合、足した値が重複しないようにループの昇降に気を付けましょう。以下に約数の実装を載せます、倍数の場合や逆変換も似た感じです。

template <class T> void divisor_transform(std::vector<T> &a) {
  int n = a.size();
  for (int i = n -1; i > 0; --i) {
    for (int j = i; (j += i) < n;) {
      a[j] += a[i];
    }
  }
}


高速ゼータ変換みたいにする  O ( N \log ( \log ( N ) ) )

 本題。先ほどの実装は普通に全列挙しているので、高速ゼータ変換で言えば  3 ^ n で全列挙する方法になります。これを各素因数について高速ゼータ変換みたいに独立に処理すると高速化できます。
 高速ゼータ変換では、各bitについて「そのbit以外が等しい添え字同士で演算する」を行うと計算が出来ている、という具合でした。同様に、各素数について「その素数以外の素因数が数も含めて等しい数同士で演算する」をします。ただし素因数を持つ数はbitと違い0/1の2値ではないので、少し計算を変える必要があります。素数  p について処理するときを考えてみましょう。 p 以外の素因数が等しいような添え字の集合は  \lbrace c , c p , c p ^ 2 , c p ^ 3 \ldots \rbrace となっています。この中だけで求める変換を行えばよいので、 p の指数が小さいほうから累積和を取って行けばよいことが分かります。
 結局これを整理すると  N 以下の素数  p について、 a _ k a _ {k p} に加算するという操作を  k の小さい方から行えばよいことが分かります。 k の範囲は  k p \lt N までであるため、これはまさにエラトステネスの篩そのものであることが分かります。よって計算量は  O ( N \log ( \log ( N ) ) ) となります。

template <class T> void divisor_transform(std::vector<T> &a) {
  int n = a.size();
  std::vector<bool> sieve(n, true);
  for (int p = 2; p < n; ++p) {
    if (sieve[p]) {
      for (int k = 1; k * p < n; ++k) {
        sieve[k * p] = false;
        a[k * p] += a[k];
      }
    }
  }
}


GCD Sum を解いてみる

 高速ゼータ変換同様、 \gcd \rm lcm での畳み込みみたいなのが出来ます。
Solution: 22095570 | CodeChef

諸注意

・上記実装では読みやすさを意識して vector<bool> や int 等を使っていますが、多分遅いです。
・a.size() == 0 等のコーナーケースを処理していない実装です
・a[0] については無視していますが、GCD Sum 等の問題ではちゃんと処理に入れたほうが都合がよいです。