形式的 Dirichlet 級数
(number-theory/formal-dirichlet-series.hpp)

• View this file on GitHub
• Last update: 2025-11-20 21:02:27+09:00
• Include: #include "number-theory/formal-dirichlet-series.hpp"

形式的 Dirichlet 級数ライブラリ．

prefix の具体的な値を管理する．

形式的 Dirichlet 級数

$R$ を環とする．

$R$ の数列 $A,B:\mathbb{N}\to R$ の Dirichlet 積 $A*B:\mathbb{N}\to R$ を以下で定める． \((A*B)(k)=\sum_{ij=k}A(i)B(j)\)

$R^{\mathbb{N}}$ に対して加算を成分毎の加算，乗算を Dirichlet 積とすると環になる．これを形式的 Dirichlet 級数環と呼ぶことにする． 以下特に断りのない場合積の記号 $*$ は省略する．

Dirichlet 級数の記法を流用し $A\in R^{\mathbb{N}}$ を以下のようにも表す： \(\sum_{n=1}^{\infty}\frac{A(n)}{n^{-s}}.\) また $[n^{-s}]A=A(n)$ とする．

この記法を用いれば加法の単位元は $0$，乗法の単位元は $1$ である． また Dirichlet 級数 $A$ の逆元が存在する必要十分条件は $[1^{-s}]A$ が逆元を持つことである．

乗算/除算

$A,B$ の前 $N$ 項がわかっているとき，$AB$ の前 $N$ 項が $O(N\log N)$ 時間で計算できる．

exp/log

• [https://atcoder.jp/contests/abc428/editorial/14241]

形式的 Dirichlet 級数 $A$ について $[1^{-s}]A$ が適当な条件を満たすとき形式的冪級数 $\exp x,\log x$ との合成によって $\exp A,\log A$ を考えることができる．

$A$ から $\exp A$ や $\log A$ を高速に計算できないだろうか．

ここで $n$ の重複を込めた素因数の個数を $\Omega(n)$ として，演算子 $\mathfrak{D}$ を \(\mathfrak{D}\left(\sum_{n\geq 1}\frac{a_n}{n^s}\right)=\sum_{n\geq 1}\frac{a_n\Omega(n)}{n^s}\)

によって定める．これは微分に対応する演算になっており，以下の性質を満たす．

• $\mathfrak{D}(A+B)=(\mathfrak{D} A)+(\mathfrak{D} B)$．
• $\mathfrak{D}(AB)=(\mathfrak{D} A)B+A(\mathfrak{D} B)$．
• $\mathfrak{D}(A^k)=k(\mathfrak{D} A)A^{k-1}$．
• 形式的冪級数 $f$ に対し，合成 $f(A)$ が定義されるならば $\mathfrak{D}(f(A))=f’(A)(\mathfrak{D} A)$．

逆に上二つの性質を満たす演算子は $\mathfrak{D}(p^{-s})$ を定めれば決まる． $\mathfrak{D}(p^{-s})=p^{-s}$ とすれば上記の $\mathfrak{D}$ が得られる．

従って以下が成り立つ． \(\mathfrak{D}\exp A=(\mathfrak{D} A)\exp A,\quad \mathfrak{D}\log A=\frac{\mathfrak{D} A}{A}\)

これを用いれば $A$ の前 $N$ 項がわかっているとき $\exp,\log$ の前 $N$ 項は $O(N\log N)$ 時間で計算できる．

Code

#pragma once

template <class mint>
struct FormalDirichletSeries : vector<mint> {
  using vector<mint>::vector;
  using FDS = FormalDirichletSeries;
  FDS& operator+=(const mint& v) {
    if (1 < this->size()) (*this)[1] += v;
    return *this;
  }
  FDS& operator+=(const FDS& r) {
    if (r.size() > this->size()) this->resize(r.size());
    for (int i = 1; i < (int)r.size(); i++) (*this)[i] += r[i];
    return *this;
  }
  FDS& operator-=(const mint& v) {
    if (1 < this->size()) (*this)[1] -= v;
    return *this;
  }
  FDS& operator-=(const FDS& r) {
    if (r.size() > this->size()) this->resize(r.size());
    for (int i = 1; i < (int)r.size(); i++) (*this)[i] -= r[i];
    return *this;
  }
  FDS& operator*=(const mint& v) {
    for (int i = 1; i < (int)this->size(); i++) (*this)[i] *= v;
    return *this;
  }
  FDS& operator*=(const FDS& r) {
    FDS p(this->size(), 0);
    for (int i = 1; i < (int)r.size(); i++)
      for (int j = 1; i * j < (int)p.size(); j++)
        p[i * j] += r[i] * (*this)[j];
    return *this = p;
  }
  FDS& operator/=(const FDS& r) {
    mint c = r[1].inv();
    for (int i = 1; i < (int)this->size(); i++) {
      (*this)[i] *= c;
      for (int j = 2; j < (int)r.size() && i * j < (int)this->size(); j++)
        (*this)[i * j] -= (*this)[i] * r[j];
    }
    return *this;
  }
  FDS operator+(const mint& v) const { return FDS(*this) += v; }
  FDS operator+(const FDS& r) const { return FDS(*this) += r; }
  FDS operator-(const mint& v) const { return FDS(*this) -= v; }
  FDS operator-(const FDS& r) const { return FDS(*this) -= r; }
  FDS operator*(const mint& v) const { return FDS(*this) *= v; }
  FDS operator*(const FDS& r) const { return FDS(*this) *= r; }
  FDS operator/(const FDS& r) const { return FDS(*this) /= r; }
  FDS operator-() const {
    FDS ret(this->size());
    for (int i = 1; i < (int)this->size(); i++) ret[i] = -(*this)[i];
    return ret;
  }

  static vector<mint> inv;
  static void init_inv() {
    if (inv.empty()) {
      inv.assign(33, 0);
      inv[1] = mint(1);
      auto mod = mint::get_mod();
      for (int i = 2; i < (int)inv.size(); i++) inv[i] = (-inv[mod % i]) * (mod / i);
    }
  }
  FDS integral() const {
    int n = this->size();
    if (n <= 1) return FDS(n);
    calc_lpf(n);
    FDS ret(n);
    init_inv();
    for (int i = 2; i < ret.size(); i++) ret[i] = (*this)[i] * inv[npf[i]];
    return ret;
  }
  FDS diff() const {
    calc_lpf(this->size());
    FDS ret(this->size());
    for (int i = 1; i < (int)this->size(); i++) ret[i] = (*this)[i] * npf[i];
    return ret;
  }
  FDS log(int n = -1) const {
    if (n == -1) n = (int)this->size();
    if (n <= 1) return FDS(n);
    assert((*this)[1] == mint(1));
    FDS f = *this;
    f.resize(n);
    return (f.diff() / f).integral();
  }
  FDS exp(int n = -1) const {
    if (n == -1) n = (int)this->size();
    if (n <= 1) return FDS(n);
    assert((*this)[1] == mint(0));
    auto df = this->diff();
    FDS ret(n);
    if (1 < n) ret[1] = 1;
    init_inv();
    for (int i = 1; i < n; i++) {
      if (i > 1) ret[i] *= inv[npf[i]];
      for (int j = 2; j < df.size() && i * j < n; j++) ret[i * j] += df[j] * ret[i];
    }
    return ret;
  }

  static vector<int> npf;
  static vector<int> lpf;
  static void calc_lpf(int n) {
    if (lpf.empty()) {
      lpf = {0, 1, 2, 3, 2, 5, 2, 7};
      npf = {0, 0, 1, 1, 2, 1, 2, 1};
    }
    int l = lpf.size(), r = l;
    while (r <= n) r <<= 1;
    if (l == r) return;
    lpf.resize(r);
    npf.resize(r);
    for (int i = 2; i < l; i++) {
      if (lpf[i] != i) continue;
      for (int j = i * 2; j < r; j += i)
        if (lpf[j] == 0) lpf[j] = i;
    }
    for (int i = l; i < r; i++) {
      if (lpf[i] != 0) continue;
      lpf[i] = i;
      for (int j = i * 2; j < r; j += i)
        if (lpf[j] == 0) lpf[j] = i;
    }
    for (int i = l; i < r; i++)
      npf[i] = npf[i / lpf[i]] + 1;
  }
};
template <typename mint>
vector<mint> FormalDirichletSeries<mint>::inv = {};
template <typename mint>
vector<int> FormalDirichletSeries<mint>::npf = {};
template <typename mint>
vector<int> FormalDirichletSeries<mint>::lpf = {};
/**
 * @brief 形式的 Dirichlet 級数
 * @docs docs/number-theory/formal-dirichlet-series.md
 */

#line 2 "number-theory/formal-dirichlet-series.hpp"

template <class mint>
struct FormalDirichletSeries : vector<mint> {
  using vector<mint>::vector;
  using FDS = FormalDirichletSeries;
  FDS& operator+=(const mint& v) {
    if (1 < this->size()) (*this)[1] += v;
    return *this;
  }
  FDS& operator+=(const FDS& r) {
    if (r.size() > this->size()) this->resize(r.size());
    for (int i = 1; i < (int)r.size(); i++) (*this)[i] += r[i];
    return *this;
  }
  FDS& operator-=(const mint& v) {
    if (1 < this->size()) (*this)[1] -= v;
    return *this;
  }
  FDS& operator-=(const FDS& r) {
    if (r.size() > this->size()) this->resize(r.size());
    for (int i = 1; i < (int)r.size(); i++) (*this)[i] -= r[i];
    return *this;
  }
  FDS& operator*=(const mint& v) {
    for (int i = 1; i < (int)this->size(); i++) (*this)[i] *= v;
    return *this;
  }
  FDS& operator*=(const FDS& r) {
    FDS p(this->size(), 0);
    for (int i = 1; i < (int)r.size(); i++)
      for (int j = 1; i * j < (int)p.size(); j++)
        p[i * j] += r[i] * (*this)[j];
    return *this = p;
  }
  FDS& operator/=(const FDS& r) {
    mint c = r[1].inv();
    for (int i = 1; i < (int)this->size(); i++) {
      (*this)[i] *= c;
      for (int j = 2; j < (int)r.size() && i * j < (int)this->size(); j++)
        (*this)[i * j] -= (*this)[i] * r[j];
    }
    return *this;
  }
  FDS operator+(const mint& v) const { return FDS(*this) += v; }
  FDS operator+(const FDS& r) const { return FDS(*this) += r; }
  FDS operator-(const mint& v) const { return FDS(*this) -= v; }
  FDS operator-(const FDS& r) const { return FDS(*this) -= r; }
  FDS operator*(const mint& v) const { return FDS(*this) *= v; }
  FDS operator*(const FDS& r) const { return FDS(*this) *= r; }
  FDS operator/(const FDS& r) const { return FDS(*this) /= r; }
  FDS operator-() const {
    FDS ret(this->size());
    for (int i = 1; i < (int)this->size(); i++) ret[i] = -(*this)[i];
    return ret;
  }

  static vector<mint> inv;
  static void init_inv() {
    if (inv.empty()) {
      inv.assign(33, 0);
      inv[1] = mint(1);
      auto mod = mint::get_mod();
      for (int i = 2; i < (int)inv.size(); i++) inv[i] = (-inv[mod % i]) * (mod / i);
    }
  }
  FDS integral() const {
    int n = this->size();
    if (n <= 1) return FDS(n);
    calc_lpf(n);
    FDS ret(n);
    init_inv();
    for (int i = 2; i < ret.size(); i++) ret[i] = (*this)[i] * inv[npf[i]];
    return ret;
  }
  FDS diff() const {
    calc_lpf(this->size());
    FDS ret(this->size());
    for (int i = 1; i < (int)this->size(); i++) ret[i] = (*this)[i] * npf[i];
    return ret;
  }
  FDS log(int n = -1) const {
    if (n == -1) n = (int)this->size();
    if (n <= 1) return FDS(n);
    assert((*this)[1] == mint(1));
    FDS f = *this;
    f.resize(n);
    return (f.diff() / f).integral();
  }
  FDS exp(int n = -1) const {
    if (n == -1) n = (int)this->size();
    if (n <= 1) return FDS(n);
    assert((*this)[1] == mint(0));
    auto df = this->diff();
    FDS ret(n);
    if (1 < n) ret[1] = 1;
    init_inv();
    for (int i = 1; i < n; i++) {
      if (i > 1) ret[i] *= inv[npf[i]];
      for (int j = 2; j < df.size() && i * j < n; j++) ret[i * j] += df[j] * ret[i];
    }
    return ret;
  }

  static vector<int> npf;
  static vector<int> lpf;
  static void calc_lpf(int n) {
    if (lpf.empty()) {
      lpf = {0, 1, 2, 3, 2, 5, 2, 7};
      npf = {0, 0, 1, 1, 2, 1, 2, 1};
    }
    int l = lpf.size(), r = l;
    while (r <= n) r <<= 1;
    if (l == r) return;
    lpf.resize(r);
    npf.resize(r);
    for (int i = 2; i < l; i++) {
      if (lpf[i] != i) continue;
      for (int j = i * 2; j < r; j += i)
        if (lpf[j] == 0) lpf[j] = i;
    }
    for (int i = l; i < r; i++) {
      if (lpf[i] != 0) continue;
      lpf[i] = i;
      for (int j = i * 2; j < r; j += i)
        if (lpf[j] == 0) lpf[j] = i;
    }
    for (int i = l; i < r; i++)
      npf[i] = npf[i / lpf[i]] + 1;
  }
};
template <typename mint>
vector<mint> FormalDirichletSeries<mint>::inv = {};
template <typename mint>
vector<int> FormalDirichletSeries<mint>::npf = {};
template <typename mint>
vector<int> FormalDirichletSeries<mint>::lpf = {};
/**
 * @brief 形式的 Dirichlet 級数
 * @docs docs/number-theory/formal-dirichlet-series.md
 */

Back to top page