noshi91のメモ

データ構造のある風景

添え字 gcd での畳み込みで AGC038-C を解く

概要

数列  a, b に対し、その gcd 畳み込み*1  c を以下のように定めます

 \displaystyle c_{k} := \sum_{\gcd(i, j) = k} a_{i} * b_{j}

この畳み込みの計算方法と、応用として AGC038-C を解く手順を示します。



数列の変換で要素毎の乗算に変換

数列から数列への関数  f を導入します。  f は倍数の和を取る変換で、以下のように定義します。

 \displaystyle f(x)_{k} := \sum_{k \mid i} x_{i}

高速ゼータ変換を知っている人なら、上位集合に対する高速ゼータ変換の類似物と言えば分かりやすいかもしれません。 実は、 a b の gcd 畳み込み  c に対して以下の式が成立します。

 f(c)_{i} = f(a)_{i} * f(b)_{i}

よって  f の逆変換  f^{-1} が存在すれば

 c = f^{-1}(f(a) * f(b)) ( * は要素毎の乗算)

によって  c を計算することが出来ます。 このような変換を考えることは、 \gcd に限らず  | + による他種の畳み込みにも共通しています。

変換の概説

 \gcd(i, j) = k が成立しているとき、 k \mid i かつ  k \mid j です。 条件を緩めて、 k \mid i かつ  k \mid j なる  i,j について  \sum a_{i} * b_{j} ならば容易に計算することが出来ます。 これは  i,j が独立なので  f(a)_{k} * f(b)_{k} となるためです。
このような  i,j について  \gcd(i,j) = k は必ずしも成立しませんが、 k \mid \gcd(i,j) が成立します。 これをよく見ると、 f(c)_{k} に等しいことが分かります。



 f,f^{-1} \Theta(N \log N) で計算する


 f

 k について  k \mid i なる  i を全て足しても、調和級数により計算量は  \Theta(N \log N) となります。

void f(std::vector<int> &a) {
  int n = a.size();
  for (int k = 1; k < n; ++k) {
    for (int i = k * 2; i < n; i += k) {
      a[k] += a[i];
    }
  }
}

 f^{-1} を考える都合により、in-place な実装例を示しました。
値を重複して足し合わせないよう、 k のループは昇順に回す必要があります。

 f^{-1}

 fアルゴリズムをそのまま逆順になぞればよいです。
同様に、 k の昇降に気を付けてください。

void inv_f(std::vector<int> &a) {
  int n = a.size();
  for (int k = n - 1; k >= 1; --k) {
    for (int i = k * 2; i < n; i += k) {
      a[k] -= a[i];
    }
  }
}

 i のループは並列なのでどちらに回しても大丈夫です。



AGC038-C を解く

数列  x x_{i} := i * (i が A に含まれる個数) として、 x x を gcd 畳み込みしたものを  y とします。

 \displaystyle y_{k} = \sum_{\gcd(A_{i},A_{j}) = k} A_{i} * A_{j}

であるため、基本的には

 \displaystyle \sum \frac{y_{k}}{k}

を計算すればよいです。
ただし  i \lt j に反する  i=j,i \gt j が足されていることに注意して、補正します。

以上を実装した提出がこちらです。

Submission #7659496 - AtCoder Grand Contest 038

先頭に貼ってある modint 構造体については以下の記事を参照してください。

modint 構造体を使ってみませんか? (C++) - noshi91のメモ



関連する文献

ゼータ変換・メビウス変換を理解する - Qiita

gcd に限らず、ゼータ変換とメビウス変換全般を説明した丁寧で分かりやすい記事です。



更なる高速化  \Theta(N (\log(\log(N)))

高速ゼータ変換の約数版 - noshi91のメモ

こちらのコードなどを参照してください。

*1:正式名称が分かっていません、ご存知の方がいらっしゃいましたら教えてください