lcm 畳み込み
(convolution/lcm.hpp)
- View this file on GitHub
- Last update: 2025-10-10 17:35:46+09:00
- Include:
#include "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)$ を返す.
- $c_0$ は未定義.
- $1\leq k\leq N$ に対し $c_k=\sum_{\gcd(i,j)=k}$.
$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
*/