Skip to the content.

:heavy_check_mark: lcm 畳み込み
(convolution/lcm.hpp)

GCD 畳み込み

長さ $\red{N+1}$ の vector $(a_0,a_1,\dots,a_N),(b_0,b_1,\dots,b_N)$ を受け取り,以下を満たす長さ $N+1$ の vector $(c_0,c_1,\dots,c_N)$ を返す.

$a$ を倍数ゼータ変換した列を $A$ とする.つまり $A_i=\sum_{i\mid j}a_j$ とする.$B,C$ も同様とする.

このとき $k\mid\gcd(i,j)\iff k\mid i\land k\mid j$ を用いると

\[\begin{align*} C_k &=\sum_{k\mid l}c_l\\ &=\sum_{k\mid\gcd(i,j)}a_ib_j\\ &=\sum_{k\mid i}\sum_{k\mid j}a_ib_j\\ &=A_kB_k \end{align*}\]

となる.従って $c$ は $a,b$ を倍数ゼータ変換し,各点積をとって倍数メビウス変換すれば得られる.

LCM 畳み込み

倍数ゼータ変換の代わりに約数ゼータ変換を用いれば LCM について畳み込みできる.

lcm の場合,$\operatorname{lcm}(i,j)\leq N$ とは限らない.

Depends on

Verified with

Code

#pragma once

#include "../number-theory/divisor-multiple-transform.hpp"

template <class T>
vector<T> LcmConvolution(vector<T> a, vector<T> b) {
  assert(a.size() == b.size());
  if (a.empty()) return {};
  DivisorZetaTransform(a);
  DivisorZetaTransform(b);
  for (int i = 0; i < a.size(); i++) a[i] *= b[i];
  DivisorMobiusTransform(a);
  return a;
}

/**
 * @brief lcm 畳み込み
 * @docs docs/convolution/gcd-lcm.md
 */
#line 2 "convolution/lcm.hpp"

#line 2 "number-theory/divisor-multiple-transform.hpp"

template <class T>
void DivisorZetaTransform(vector<T>& a) {
  if (a.empty()) return;
  int n = a.size();
  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];
      }
    }
  }
}

template <class T>
void DivisorMobiusTransform(vector<T>& a) {
  if (a.empty()) return;
  int n = a.size();
  vector<bool> sieve(n, true);
  for (int p = 2; p < n; p++) {
    if (sieve[p]) {
      for (int k = (n - 1) / p; k > 0; k--) {
        sieve[k * p] = false;
        a[k * p] -= a[k];
      }
    }
  }
}

template <class T>
void MultipleZetaTransform(vector<T>& a) {
  if (a.empty()) return;
  int n = a.size();
  vector<bool> sieve(n, true);
  for (int p = 2; p < n; p++) {
    if (sieve[p]) {
      for (int k = (n - 1) / p; k > 0; k--) {
        sieve[k * p] = false;
        a[k] += a[k * p];
      }
    }
  }
}

template <class T>
void MultipleMobiusTransform(vector<T>& a) {
  if (a.empty()) return;
  int n = a.size();
  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] -= a[k * p];
      }
    }
  }
}

/**
 * @brief 約数・倍数変換
 * @docs docs/number-theory/divisor-multiple-transform.md
 */
#line 4 "convolution/lcm.hpp"

template <class T>
vector<T> LcmConvolution(vector<T> a, vector<T> b) {
  assert(a.size() == b.size());
  if (a.empty()) return {};
  DivisorZetaTransform(a);
  DivisorZetaTransform(b);
  for (int i = 0; i < a.size(); i++) a[i] *= b[i];
  DivisorMobiusTransform(a);
  return a;
}

/**
 * @brief lcm 畳み込み
 * @docs docs/convolution/gcd-lcm.md
 */
Back to top page