add

Author: KumaCS

docs/fps/count-subset-sum.md

@@ -0,0 +1,11 @@
+整数列 $A=(A_0,A_1,\dots,A_{M-1})$ が与えられる．
+$k=0,1,\dots,N$ に対して $\sum_{i\in I}A_i=k$ を満たす $I\subseteq\{0,1,\dots,M-1\}$ の個数を $O(N\log N)$ 時間で列挙する．
+
+## アルゴリズム
+
+$\prod_{i}(1+x^{A_i})$ を $N$ 次まで計算できればよい．
+以下のように変形する．
+
+$$\prod_{i}(1+x^{A_i})=\exp\left(\log\prod_{i}(1+x^{A_i})\right)=\exp\left(\sum_{i}\log(1+x^{A_i})\right)$$
+
+$\log(1+x^{A_i})$ は $A_i$ の倍数次の項だけをもつので，$\exp$ の中身は調和級数の議論で $O(N\log N)$ 時間で計算できる．

docs/fps/polynomial-gcd.md

@@ -0,0 +1,99 @@
+多項式 $f(x),g(x)$ に対する $\gcd(f(x),g(x))$ を $n=\max(\deg f,\deg g)$ として $O(n(\log n)^2)$ 時間で計算する．
+
+また $f(x)h(x)\equiv 1\pmod{g(x)}$ となる多項式 $h(x)$ も同じ計算量で求められる．
+
+## Half-GCD 法
+
+Euclid の互除法をそのまま適用すると，一回の反復で次数が $1$ しか減少しない場合に $\Omega(n^2\log n)$ 時間かかる．
+解決するのが Half-GCD 法．
+
+- [https://codeforces.com/blog/entry/101850]
+- [https://dl.acm.org/doi/10.1145/800125.804045]
+
+$f,g$ の GCD は互除法で求められる．
+
+- $f_0=f,f_1=g$ とする（$\deg f\geq\deg g$）．
+- $i=1,2,\dots$ について順に $f_{i-1}$ を $f_i$ で割った商を $q_i$ とし，$f_{i+1}=f_i-f_{i+1}q_i$ と定めていく．
+- $f_{i+1}=0$ となったとき $\operatorname{GCD}(f,g)=f_i$ である．
+
+過程は行列で表示できる：$M_i=\begin{pmatrix}0&1\\1&-q_i\end{pmatrix}$ とすると $\begin{pmatrix}f_i\\f_{i+1}\end{pmatrix}=M_i\cdots M_1\begin{pmatrix}f_0\\f_1\end{pmatrix}$ である．
+
+$\deg f_t \gt \frac{1}{2}\deg f\geq \deg f_{t+1}$ となる $t$ について $\operatorname{HGCD}(f,g)=M_t\cdots M_1$ とおく．
+
+HGCD は再帰的に計算できる（証明は後述）．ただし $\deg f\geq\deg g$ とし，`f >> k` は $f\operatorname{div}x^k$ を意味する．
+
+```
+def HGCD(f, g):
+  if deg(g) <= deg(f) / 2: return I
+  k = floor(deg(f) / 2) + 1
+  M = HGCD(f >> k, g >> k)
+  (f, g)^T = M * (f, g)^T
+  if deg(g) <= deg(f) / 2: return M
+  (通常の互除法の 1 step を行う)
+  if deg(g) <= deg(f) / 2: return M
+  l = 2 * k - deg(f)
+  return HGCD(f >> l, g >> l) * M
+```
+
+2 回の HGCD ではどちらも次数が半分以下になっている．
+また多項式除算を 1 回行っている．
+従って $T(n)=2T(n/2)+n\log n$ を解いて $T(n)=n(\log n)^2$ であるから計算量は $O(n(\log n)^2)$ 時間である．
+
+
+HGCD を用いれば以下のアルゴリズムで GCD が求められる．
+ループごとに $f$ の次数が半減するので $T(n)+T(n/2)+T(n/4)+\cdots=O(n(\log n)^2)$ 時間．
+
+```
+def GCD(f, g):
+  while g != 0:
+    (f, g) = HGCD(f, g) * (f, g)
+    if g == 0: break
+    (f, g) = (g, f % g)
+  return f
+```
+
+
+### 正当性
+
+形式的 Laurent 級数（有限個の負冪の項を許す）として考えると，$f$ を $g$ で割った商 $q$ は
+
+$$q(x^{-1})=\frac{f(x^{-1})}{g(x^{-1})}\bmod x$$
+
+と表せる．
+
+
+
+> **補題.** 多項式 $f,g\ (\deg f-\deg g=k)$ に対し $f$ を $g$ で割った商 $q$ は $f,g$ のそれぞれ上から $k+1$ 次の係数にのみ依存する（特に $\deg f$ の値によらない）．
+
+<details>
+<summary>証明</summary>
+
+両辺に $x^k=x^{\deg f-\deg g}$ をかけると
+$$\begin{align*}
+x^kq(x^{-1})
+&=\frac{x^{\deg f}f(x^{-1})}{x^{\deg g}g(x^{-1})}\bmod x^{k+1}\\
+&=\frac{x^{\deg f}f(x^{-1})\bmod x^{k+1}}{x^{\deg g}g(x^{-1})\bmod x^{k+1}}\bmod x^{k+1}
+\end{align*}$$
+であるのでよい．
+(証明終)
+
+</details>
+
+---
+
+> **定理.** $(f,g)$ に対し互除法を適用して得られる行列を $M_1,M_2,\dots$ とする．非負整数 $k$ に対し，ある $t$ が存在して $\operatorname{HGCD}(f\operatorname{div}x^k,g\operatorname{div}x^k)=M_t\cdots M_1$ となる．
+
+<details>
+<summary>証明</summary>
+
+$n=\deg f$ とする．
+
+HGCD 側は $\deg f'_i\gt\frac{1}{2}(n-k)$ の間操作を行う．
+これをみたす $i$ に対して以下を示すことができればよい．
+- $\deg f_i=\deg f'_i+k$
+- $f_i$ と $f'_i$ のそれぞれ上 $\deg f_{i-1}-\lfloor\frac{n+k}{2}\rfloor+1=\deg f'_{i-1}-\lfloor\frac{n-k}{2}\rfloor+1$ 次が一致する
+
+(TODO: 証明する，実験では成り立ちそう & これより次数評価は緩められなさそう)
+(証明終)
+
+</details>

docs/fps/polynomial-interpolation.md

@@ -1,15 +1,49 @@
-$n$ 次多項式 $f$ で $0\leq i\leq n$ に対し $f(x_i)=y_i$ を満たすものを $O(n(\log n)^2)$ 時間で求める．
+高々 $n-1$ 次の多項式 $f$ で $0\leq i\lt n$ に対し $f(x_i)=y_i$ を満たすものを $O(n(\log n)^2)$ 時間で求める．
+
+また評価点が等比級数ならば $O(n\log n)$ 時間で求められる．
 
 ## アルゴリズム
 
-$g(x)=\prod_{i=0}^{n}(x-x_i)$ とすれば Lagrange 補間より以下が成り立つ．
+$g(x)=\prod_{i=0}^{n-1}(x-x_i)$ とすれば Lagrange 補間より以下が成り立つ．
 
-$$f(x)=\sum_{i=0}^{n}\frac{y_ig(x)}{g'(x_i)(x-x_i)}$$
+$$f(x)=\sum_{i=0}^{n-1}\frac{y_ig(x)}{g'(x_i)(x-x_i)}$$
 
 以下のようにして $O(n(\log n)^2)$ 時間で計算できる．
 
 - $g(x)$ は分割統治により $O(n(\log n)^2)$ 時間で計算できる．
 - $g'(x_i)$ は多点評価により $O(n(\log n)^2)$ 時間で列挙できる．
-- $\sum_{i=0}^{n}\frac{y_i}{g'(x_i)(x-x_i)}$ は有理数の和として分割統治により $O(n(\log n)^2)$ 時間で計算できる．
+- $\sum_{i=0}^{n-1}\frac{y_i}{g'(x_i)(x-x_i)}$ は有理数の和として分割統治により $O(n(\log n)^2)$ 時間で計算できる．
+
+$g(x)$ の計算過程はそれ以降にも流用できる．
+
+## 等比級数の場合
+
+- [https://judge.yosupo.jp/problem/polynomial_interpolation_on_geometric_sequence]
+- [https://noshi91.github.io/algorithm-encyclopedia/polynomial-interpolation-geometric#noredirect]
+
+$i=0,\dots,n-1$ について $f(ar^i)$ がわかっているとする．ただし $0\leq i\lt j\lt n$ ならば $ar^i\neq ar^j$．
+
+$f$ は以下のように表示できる．
+$$f(x)=\sum_{i=0}^{n-1}y_i\prod_{j\neq i}\frac{x-ar^j}{a(r^i-r^j)}$$
+
+$g(x)=x^{n-1}f(x^{-1})$ とすると $g$ は以下のように表示できる．
+$$g(x)=\sum_{i=0}^{n-1}y_i\prod_{j\neq i}\frac{1-ar^jx}{a(r^i-r^j)}$$
+
+分母は以下のように変形できる．
+$$\begin{align*}
+\prod_{j\neq i}a(r^i-r^j)
+&=\prod_{j=0}^{i-1}a(r^i-r^j)\prod_{j=i+1}^{n-1}a(r^i-r^j)\\
+&=\prod_{j=1}^{i}r^{i-j}a(r^j-1)\prod_{j=1}^{n-1-i}ar^i(1-r^j)\\
+&=(-1)^{i}r^{i(i-1)/2+i(n-1-i)}\prod_{j=1}^{i}a(1-r^j)\prod_{j=1}^{n-1-i}a(1-r^j)
+\end{align*}$$
+
+よって $w_i=y_i\prod_{j\neq i}\frac{1}{a(r^i-r^j)}$ は全ての $i$ に対し $O(n)$ 時間で求められるので以下を計算できればよい．
+$$g(x)=\sum_{i=0}^{n-1}w_i\prod_{j\neq i}(1-ar^jx)=\left(\prod_{i=0}^{n-1}(1-ar^ix)\right)\left(\sum_{i=0}^{n-1}\frac{w_i}{1-ar^ix}\right)$$
+
+$F(x)=\prod_{i=0}^{n-1}(1-ar^ix)$ とおくと $(1-ax)F(rx)=F(x)(1-ar^nx)$ より $F_i=[x^i]F(x)$ について
+$$F_0=1,\quad (1-r^i)F_i=a(r^n-r^{i-1})F_{i-1}$$
+が成り立ち，$O(n)$ 時間で $F$ は求められる．
 
-$g(x)$ の計算過程はそれ以降にも流用できる．
+$G(x)=\sum_{i=0}^{n-1}\frac{w_i}{1-ar^ix}$ とおくと
+$$G(x)=\sum_{i=0}^{n-1}\sum_{j=0}^{\infty}w_i(ar^ix)^j=\sum_{j=0}^{\infty}a^j\left(\sum_{i=0}^{n-1}w_ir^{ij}\right)x^j$$
+である．$G$ は $n-1$ 次まで求められればよいので，多項式 $\sum_{i=0}^{n-1}w_ix^i$ を $x=r^0,r^1,\dots,r^{n-1}$ で多点評価すればよい．

docs/set/power-projection-of-set-power-series.md

@@ -0,0 +1,26 @@
+$a$ を集合冪級数とする．
+$\sum_{s}w_s[x^s]a^i$ を $i=0,\dots,M-1$ について列挙する．
+
+$O(N^22^N+NM)$ 時間．
+
+- [https://maspypy.com/%e9%9b%86%e5%90%88%e3%81%b9%e3%81%8d%e7%b4%9a%e6%95%b0%e9%96%a2%e9%80%a3-3-%e5%a4%9a%e9%a0%85%e5%bc%8f%e3%81%a8%e3%81%ae%e5%90%88%e6%88%90%e3%81%ae%e8%bb%a2%e7%bd%ae]
+
+## 導入
+
+$a$ を fix したとき，合成 $f(a)$ の計算は $\sum_{i}f_i[x^s]a^i$ の $s$ についての列挙．
+
+この転置問題は $\sum_{s}w_s[x^s]a^i$ の $i$ についての列挙．
+
+多項式 $f$ との合成は以下の要領で計算していた．
+
+- $g(x)=f(x+c)$ とする
+- $g$ を EGF にする
+- $g$ と $a-c$ を合成する
+
+$g(x)=f(x+c)$ のとき $g_i=\sum_{j\geq i}\binom{j}{i}c^{j-i}f_j$ である．この転置は難しくない．
+
+また $g$ を EGF にするには $x^k$ の係数を $k!$ 倍すればよいので，特に転置は同じものである．
+
+$[x^0]a=0$ のときの EGF との合成の転置を考える．
+subset convolution $a\star b$ は $b$ を固定すれば $a$ について線形写像．
+その転置アルゴリズムは $a\star b=\operatorname{rev}(\operatorname{rev}(a)\star b)$ である．

fps/count-subset-sum.hpp

@@ -0,0 +1,27 @@
+#pragma once
+#include "modint/factorial.hpp"
+#include "fps/formal-power-series.hpp"
+
+// #{ I | sum[i in I]a[i] = k } for k=0,1,...,n
+template <class mint>
+FormalPowerSeries<mint> CountSubsetSum(vector<int> a, int n) {
+  using fact = Factorial<mint>;
+  using fps = FormalPowerSeries<mint>;
+  vector<int> c(n + 1, 0);
+  for (auto v : a)
+    if (0 <= v && v <= n) c[v]++;
+  FormalPowerSeries<mint> f(n + 1, 0);
+  for (int i = 1; i <= n; i++) {
+    if (c[i] == 0) continue;
+    mint s = c[i];
+    for (int j = 1; i * j <= n; j++) {
+      f[i * j] += s * fact::inv(j);
+      s = -s;
+    }
+  }
+  return f.exp();
+}
+/**
+ * @brief Count Subset Sum
+ * @docs docs/fps/count-subset-sum.md
+ */

fps/formal-power-series.hpp

@@ -4,31 +4,37 @@ template <class mint>
 struct FormalPowerSeries : vector<mint> {
   using vector<mint>::vector;
   using FPS = FormalPowerSeries;
-  FPS &operator+=(const FPS &r) {
+  FormalPowerSeries(const vector<mint>& r) : vector<mint>(r) {}
+  FormalPowerSeries(vector<mint>&& r) : vector<mint>(std::move(r)) {}
+  FPS& operator=(const vector<mint>& r) {
+    vector<mint>::operator=(r);
+    return *this;
+  }
+  FPS& operator+=(const FPS& r) {
     if (r.size() > this->size()) this->resize(r.size());
     for (int i = 0; i < (int)r.size(); i++) (*this)[i] += r[i];
     return *this;
   }
-  FPS &operator+=(const mint &r) {
+  FPS& operator+=(const mint& r) {
     if (this->empty()) this->resize(1);
     (*this)[0] += r;
     return *this;
   }
-  FPS &operator-=(const FPS &r) {
+  FPS& operator-=(const FPS& r) {
     if (r.size() > this->size()) this->resize(r.size());
     for (int i = 0; i < (int)r.size(); i++) (*this)[i] -= r[i];
     return *this;
   }
-  FPS &operator-=(const mint &r) {
+  FPS& operator-=(const mint& r) {
     if (this->empty()) this->resize(1);
     (*this)[0] -= r;
     return *this;
   }
-  FPS &operator*=(const mint &v) {
+  FPS& operator*=(const mint& v) {
     for (int k = 0; k < (int)this->size(); k++) (*this)[k] *= v;
     return *this;
   }
-  FPS &operator/=(const FPS &r) {
+  FPS& operator/=(const FPS& r) {
     if (this->size() < r.size()) {
       this->clear();
       return *this;
@@ -38,7 +44,7 @@ struct FormalPowerSeries : vector<mint> {
       FPS f(*this), g(r);
       g.shrink();
       mint coeff = g.at(g.size() - 1).inv();
-      for (auto &x : g) x *= coeff;
+      for (auto& x : g) x *= coeff;
       int deg = (int)f.size() - (int)g.size() + 1;
       int gs = g.size();
       FPS quo(deg);
@@ -52,19 +58,19 @@ struct FormalPowerSeries : vector<mint> {
     }
     return *this = ((*this).rev().pre(n) * r.rev().inv(n)).pre(n).rev();
   }
-  FPS &operator%=(const FPS &r) {
+  FPS& operator%=(const FPS& r) {
     *this -= *this / r * r;
     shrink();
     return *this;
   }
-  FPS operator+(const FPS &r) const { return FPS(*this) += r; }
-  FPS operator+(const mint &v) const { return FPS(*this) += v; }
-  FPS operator-(const FPS &r) const { return FPS(*this) -= r; }
-  FPS operator-(const mint &v) const { return FPS(*this) -= v; }
-  FPS operator*(const FPS &r) const { return FPS(*this) *= r; }
-  FPS operator*(const mint &v) const { return FPS(*this) *= v; }
-  FPS operator/(const FPS &r) const { return FPS(*this) /= r; }
-  FPS operator%(const FPS &r) const { return FPS(*this) %= r; }
+  FPS operator+(const FPS& r) const { return FPS(*this) += r; }
+  FPS operator+(const mint& v) const { return FPS(*this) += v; }
+  FPS operator-(const FPS& r) const { return FPS(*this) -= r; }
+  FPS operator-(const mint& v) const { return FPS(*this) -= v; }
+  FPS operator*(const FPS& r) const { return FPS(*this) *= r; }
+  FPS operator*(const mint& v) const { return FPS(*this) *= v; }
+  FPS operator/(const FPS& r) const { return FPS(*this) /= r; }
+  FPS operator%(const FPS& r) const { return FPS(*this) %= r; }
   FPS operator-() const {
     FPS ret(this->size());
     for (int i = 0; i < (int)this->size(); i++) ret[i] = -(*this)[i];
@@ -130,7 +136,7 @@ struct FormalPowerSeries : vector<mint> {
   }
   mint eval(mint x) const {
     mint r = 0, w = 1;
-    for (auto &v : *this) r += w * v, w *= x;
+    for (auto& v : *this) r += w * v, w *= x;
     return r;
   }
   FPS log(int deg = -1) const {
@@ -160,10 +166,10 @@ struct FormalPowerSeries : vector<mint> {
     return FPS(deg, mint(0));
   }
 
-  static void *ntt_ptr;
+  static void* ntt_ptr;
   static void set_ntt();
-  FPS &operator*=(const FPS &r);
-  FPS middle_product(const FPS &r) const;
+  FPS& operator*=(const FPS& r);
+  FPS middle_product(const FPS& r) const;
   void ntt();
   void intt();
   void ntt_doubling();
@@ -172,4 +178,4 @@ struct FormalPowerSeries : vector<mint> {
   FPS exp(int deg = -1) const;
 };
 template <typename mint>
-void *FormalPowerSeries<mint>::ntt_ptr = nullptr;
+void* FormalPowerSeries<mint>::ntt_ptr = nullptr;

fps/multipoint-evaluation.hpp

@@ -47,6 +47,6 @@ vector<mint> MultipointEvaluationGeometric(FormalPowerSeries<mint> f, mint a, mi
 }
 
 /**
- * @brief Multipoint Evaluation
+ * @brief 多点評価
  * @docs docs/fps/multipoint-evaluation.md
  */