約数ゼータ・メビウス変換のメモ

Last updated Apr 18, 2022

• Algorithm

以前僕が書いた記事を復習するときにとても読みにくかったので、もう一度書き直します。

倍数ゼータ・メビウス変換についても書いたので、合わせて読むといいかも。同じ問題を別視点で考えてます。

関数$f(n)$に対する約数ゼータ変換は以下の定義です。$d|n$は「$d$は$n$の約数」という意味です(dはdivisorのd)。約数の足し合わせですね。

$$ F(n) = \sum_{ d | n } f(d) $$

約数メビウス変換は、約数ゼータ変換の逆向きの操作です。つまり、$F(n)$から$f(n)$を求める操作です。

ちなみに、メビウス変換はメビウスの反転公式から メビウス関数を$\mu(n)$とすれば、

$$ f(n) = \sum_{d | n} F(d) \mu(\frac{n}{d}) $$

です。メビウス関数については蟻本にも書いてあります。

高速ゼータ変換の約数版 O(N log(log(N))) - noshi91のメモがはやいです！

ある$n, m$について、$F(\gcd(n, m))$を考えてみます。$f$は$F$をメビウス変換した関数です。

$$ F(\gcd(n, m)) = \sum_{d | gcd(n, m)} f(d) = \sum_{d | n, d | m} f(d) $$

最後の$\sum$の条件が2つあるので、これをどうにか$d | n$だけにしてみます。

$$ c(d) = \begin{cases} 1 & d | m \\\ 0 & \text{otherwise} \end{cases} $$

とすれば以下のようになります。

$$ F(\gcd(n, m)) = \sum_{d | n} f(d) c(d) $$

こうすると例えば、$\sum_i F(\gcd(n, m_i))$の計算は

$$ \sum_i F(\gcd(n, m_i)) = \sum_{d | n} f(d) (\sum_i c_i(d)) $$

とまとめて計算できます。

AGC038 C - LCMs

$$ \begin{align} & \sum_{i = 0}^{N - 1} \sum_{j = 0}^{i - 1} \operatorname{lcm}(A_i, A_j) \\\ = & \sum_{i = 0}^{N - 1} \sum_{j = 0}^{i - 1} A_i \frac{A_j}{\gcd(A_i, A_j)} \\\ \end{align} $$

$F(n) = \frac{1}{n}$として、その約数メビウス変換を$f(n)$としましょう。すると、

$$ \begin{align} = & \sum_{i = 0}^{N - 1} A_i \sum_{j = 0}^{i - 1} F(\gcd(A_i, A_j)) A_j \\\ = & \sum_{i = 0}^{N - 1} A_i \sum_{j = 0}^{i - 1} (\sum_{d | \gcd(A_i, A_j)} f(d) A_j) \\\ = & \sum_{i = 0}^{N - 1} A_i \sum_{j = 0}^{i - 1} (\sum_{d | A_i, d | A_j} f(d) A_j) \\\ \end{align} $$

ここで、上の$c(d)$のような係数テクを導入します。

$$ s_j(d) = \begin{cases} A_j & d | A_j \\\ 0 & \text{otherwise} \end{cases} $$

とすれば、$d | A_j$の条件を取り除けて

$$ \begin{align} = & \sum_{i = 0}^{N - 1} A_i \sum_{j = 0}^{i - 1} (\sum_{d | A_i} f(d) s_j(d)) \\\ = & \sum_{i = 0}^{N - 1} A_i \sum_{d | A_i} f(d) (\sum_{j = 0}^{i - 1} s_j(d)) \\\ \end{align} $$

$\sum_{j = 0}^{i - 1} s_j(d)$は$i$を進めながら更新できます。計算量は、$f$を求める約数メビウス変換が$O(A \log\log A)$、約数の数は$O(\sqrt A)$で抑えられるので、全体で$O(A \log\log A + N \sqrt A)$です。

提出コード

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

const i64 MAX = 1e6 + 1;
 
int main() {
  i64 N;
  cin >> N;
  vector<i64> A(N);
  rep(i,0,N) cin >> A[i];
 
  vector<fp> f(MAX);
  rep(i,1,MAX) {
    f[i] = fp(i).inv();
  }
  inverse_divisor_transform(f);
  
  vector<fp> sum(MAX);
  fp ans;
  for(auto a: A) {
    fp res;
    for(i64 d = 1; d * d <= a; d++) {
      if(a % d == 0) {
        res += f[d] * sum[d];
        sum[d] += fp(a);
        i64 dd = a / d;
        if(dd != d) {
          res += f[dd] * sum[dd];
          sum[dd] += fp(a);
        }
      }
    }
    ans += res * fp(a);
  }
  cout << ans << endl;
}

約数列挙にosa_k法を使うとめっちゃはやくなります。

提出コード

ABC248 G - GCD cost on the treeのパスグラフ版から考えて、木に応用しようとしたけど失敗しちゃったのでパスグラフの解法だけでもメモっておきます…

$N$要素からなる数列$A_0, \cdots, A_{N - 1}$が与えられて以下を求める問題です。

$$ \sum_{i = 0}^{N - 1} \sum_{j = 0}^{i - 1} (i - j + 1) \gcd(A_j, \cdots, A_i) $$

式変形していきます。

$$ \begin{align} & = \sum_{i = 0}^{N - 1} \sum_{j = 0}^{i - 1} (i + 1) \gcd(A_j, \cdots, A_i) \\\ & - \sum_{i = 0}^{N - 1} \sum_{j = 0}^{i - 1} j \gcd(A_j, \cdots, A_i) \end{align} $$

$F(n) = n$として、その約数メビウス変換を$f(n)$とします。まずは前半部分を計算しましょう。 $$ \begin{align} & \sum_{i = 0}^{N - 1} \sum_{j = 0}^{i - 1} (i + 1) \gcd(A_j, \cdots, A_i) \\\ = & \sum_{i = 0}^{N - 1} (i + 1) \sum_{j = 0}^{i - 1} F(\gcd(A_j, \cdots, A_i)) \\\ = & \sum_{i = 0}^{N - 1} (i + 1) \sum_{j = 0}^{i - 1} (\sum_{d | A_j, \cdots d | A_i} f(d)) \\\ \end{align} $$

先にあったような係数テクをします。

$$ s_j(d) = \begin{cases} 1 & d | A_j, \cdots d | A_{i - 1} \\\ 0 & \text{otherwise} \end{cases} $$

とすると、

$$ \begin{align} = & \sum_{i = 0}^{N - 1} (i + 1) \sum_{j = 0}^{i - 1} (\sum_{d | A_i} f(d) s_j(d)) \\\ = & \sum_{i = 0}^{N - 1} (i + 1) \sum_{d | A_i} f(d) (\sum_{j = 0}^{i - 1} s_j(d)) \\\ \end{align} $$

$\sum_{j = 0}^{i - 1} s_j(d)$は、$i$を進めながら更新できます。$S_{i}(d) = \sum_{j = 0}^{i - 1} s_j(d)$としたとき、

$$ S_{i + 1}(d) = \begin{cases} 1 + S_i(d) & d | A_i \\\ 0 & \text{otherwise} \end{cases} $$

とすればよいです。実装方法としては、$S_{i + 1}(d) \neq 0$となる$d$の数が$O(\sqrt A)$で抑えられることを活かす方法と、最後に$S(d)$を更新したときの$i$を記録しておく方法がありそうです。

後半部分も計算します。

$$ \begin{align} & \sum_{i = 0}^{N - 1} \sum_{j = 0}^{i - 1} j \gcd(A_j, \cdots, A_i) \\\ = & \sum_{i = 0}^{N - 1} \sum_{j = 0}^{i - 1} j F(\gcd(A_j, \cdots, A_i)) \\\ = & \sum_{i = 0}^{N - 1} \sum_{j = 0}^{i - 1} (\sum_{d | A_j, \cdots d | A_i} f(d) j) \\\ \end{align} $$

先にあったような係数テクをします。

$$ t_j(d) = \begin{cases} j & d | A_j, \cdots d | A_{i - 1} \\\ 0 & \text{otherwise} \end{cases} $$

とすると、

$$ \begin{align} = & \sum_{i = 0}^{N - 1} (i + 1) \sum_{j = 0}^{i - 1} (\sum_{d | A_i} f(d) t_j(d)) \\\ = & \sum_{i = 0}^{N - 1} (i + 1) \sum_{d | A_i} f(d) (\sum_{j = 0}^{i - 1} t_j(d)) \\\ \end{align} $$

$T_{i} = \sum_{j = 0}^{i - 1} t_j(d)$も同様に、$i$を進めながら更新できます。

$$ T_{i + 1}(d) = \begin{cases} i + T_i(d) & d | A_i \\\ 0 & \text{otherwise} \end{cases} $$

これでパスグラフに関しては解けているはず…？下に検証コードを載せておきます。

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100

const i64 MAX = 1e5 + 1;

fp solve(i64 N, vector<i64> A) {
  vector<fp> f(MAX);
  rep(i,0,MAX) f[i] = fp(i);
  inverse_divisor_transform(f);

  fp left;
  {
  vector<pair<int, fp>> sum(MAX, { -1, fp(0) });
  rep(i,0,N) {
    fp res;
    for(i64 d = 1; d * d <= A[i]; d++) {
      if(A[i] % d == 0) {
        auto& [b, s] = sum[d];
        res += (b + 1 == i ? s : fp(0)) * f[d];
        if(b + 1 != i) {
          s = 0;
        }
        b = i;
        s += fp(1);
        if(A[i] / d != d) {
          i64 dd = A[i] / d;
          auto& [b, s] = sum[dd];
          res += (b + 1 == i ? s : fp(0)) * f[dd];
          if(b + 1 != i) {
            s = 0;
          }
          b = i;
          s += fp(1);
        }
      }
    }
    left += res * fp(i + 1);
  }
  }
  {
  vector<pair<int, fp>> sum(MAX, { -1, fp(0) });
  rep(i,0,N) {
    fp res;
    for(i64 d = 1; d * d <= A[i]; d++) {
      if(A[i] % d == 0) {
        auto& [b, s] = sum[d];
        res += (b + 1 == i ? s : fp(0)) * f[d];
        if(b + 1 != i) {
          s = 0;
        }
        b = i;
        s += fp(i);
        if(A[i] / d != d) {
          i64 dd = A[i] / d;
          auto& [b, s] = sum[dd];
          res += (b + 1 == i ? s : fp(0)) * f[dd];
          if(b + 1 != i) {
            s = 0;
          }
          b = i;
          s += fp(i);
        }
      }
    }
    left -= res;
  }
  }
  return left;
}

i64 gcd(i64 a, i64 b) {
  if(b == 0) return a;
  return gcd(b, a % b);
}

fp naive(i64 N, std::vector<i64> A) {
  fp ans;
  rep(i,0,N) {
    i64 G = A[i];
    rep(j,i + 1,N) {
      G = gcd(G, A[j]);
      ans += fp(G) * fp(j - i + 1);
    }
  }
  return ans;
}

int main() {
  rep(s,0,1000) {
    i64 N = 1000;
    mt19937 mt(s);
    uniform_int_distribution<> dist(1, (int)1e5);
    vector<i64> A(N);
    rep(i,0,N) A[i] = dist(mt);
    cout << "seed: " << s << endl;
    //cout << N << endl;
    //cout << A << endl;
    fp so, na;
    cout << "solve: " << (so = solve(N, A)) << endl;
    cout << "naive: " << (na = naive(N, A)) << endl;
    assert(so == na);
  }
}

木バージョンが解けないので、がんばる