update

Author: KumaCS

.verify-helper/docs/index.md

@@ -1,2 +0,0 @@
-[![Actions Status](https://github.com/KumaCS/library-cpp/workflows/verify/badge.svg)](https://github.com/KumaCS/library-cpp/actions)
-[![GitHub Pages](https://img.shields.io/static/v1?label=GitHub+Pages&message=+&color=brightgreen&logo=github)](https://KumaCS.github.io/library-cpp/)

data-structure/binary-indexed-tree.hpp

@@ -5,7 +5,7 @@ struct BinaryIndexedTree {
   int size;
   vector<T> data;
   BinaryIndexedTree(int n) : size(n), data(n) {}
-  T add(int p, T x) {
+  void add(int p, T x) {
     for (p++; p <= size; p += p & -p) data[p - 1] += x;
   }
   T sum(int p) {
@@ -15,9 +15,9 @@ struct BinaryIndexedTree {
   }
   T sum(int l, int r) { return sum(r) - sum(l); }
 
-  template <bool (*f)(S)>
+  template <bool (*f)(T)>
   int lower_bound() const {
-    return lower_bound([](S x) { return f(x); });
+    return lower_bound([](T x) { return f(x); });
   }
   template <class F>
   int lower_bound(F f) const {
@@ -26,10 +26,8 @@ struct BinaryIndexedTree {
     while (r < size) r = r << 1;
     T s = 0;
     for (int len = r; len > 0; len >>= 1) {
-      if (x + len < size && !f(s + data[x + len])) {
-        s += data[x + len];
-        x += len;
-      }
+      if (x + len < size && !f(s + data[x + len]))
+        s += data[x += len];
     }
     return x + 1;
   }

data-structure/integer-set.hpp

@@ -1,6 +1,6 @@
 #pragma once
 
-template <class T = unsigned int, class U = unsigned long long>
+template <class T = int, class U = unsigned long long>
 class IntegerSet {
  private:
   const static T B = 6, W = 64, MASK = W - 1;
@@ -24,11 +24,12 @@ class IntegerSet {
     while (n > 0);
     start = vector<T>(1, 0);
     start.reserve(len.size() + 1);
-    for (const auto v : len) start.push_back(start.back() + v);
+    for (auto v : len) start.push_back(start.back() + v);
     data = vector<U>(start.back());
   }
   IntegerSet(string s = "") {
-    T n = size = s.size();
+    size = s.size();
+    T n = size;
     vector<T> len;
     do len.push_back((n >>= B) + 1);
     while (n > 0);
@@ -58,11 +59,11 @@ class IntegerSet {
     }
   }
   bool contains(T x) const { return (data[x >> B] >> (x & MASK)) & 1; }
-  optional<T> min(T x) const {
+  T min(T x) const {
     T d = 0, i = x;
     while (true) {
-      if (d + 1 >= start.size()) return nullopt;
-      if ((i >> B) >= start[d + 1] - start[d]) return nullopt;
+      if (d + 1 >= start.size()) return -1;
+      if ((i >> B) >= start[d + 1] - start[d]) return -1;
       U m = data[start[d] + (i >> B)] & ((~U(0)) << (i & MASK));
       if (m == 0) {
         d++;
@@ -77,11 +78,11 @@ class IntegerSet {
     }
     return i;
   }
-  optional<T> max(T x) const {
+  T max(T x) const {
     T d = 0, i = x;
     while (true) {
-      if (i >= size) return nullopt;
-      if (d >= data.size()) return nullopt;
+      if (i < 0) return -1;
+      if (d >= data.size()) return -1;
       U m = data[start[d] + (i >> B)] & ~((~U(1)) << (i & MASK));
       if (m == 0) {
         d++;

docs/flow/max-flow.md

@@ -8,7 +8,7 @@ push-relabel algorithm により $n$ 頂点 $m$ 辺のとき $O(n^2m)$ 時間．
   - $\sum_{v\in V}f(v,u)-\sum_{v\in V}f(u,v)\geq 0$
 - 頂点 $u$ の excess flow とは $e(u)=\sum_{v\in V}f(v,u)-\sum_{v\in V}f(u,v)$ のこと．
 - 頂点 $u$ が overflow しているとは $e(u)\gt 0$ であること．
-- height function $h:V\to\mathbb{N}$ は $h(s)=|V|,h(t)=0$ およびすべての残余辺 $(u,v)\in E_f$ に対し $h(u)\leq h(v)+1$ を満たすもの．
+- height function $h:V\to\mathbb{N}$ は $h(s)=\lvert V\rvert,h(t)=0$ およびすべての残余辺 $(u,v)\in E_f$ に対し $h(u)\leq h(v)+1$ を満たすもの．
 
 preflow を更新する push，height function を更新する relabel を組み合わせたアルゴリズム．
 

docs/fps/polynomial-interpolation.md

@@ -0,0 +1,15 @@
+$n$ 次多項式 $f$ で $0\leq i\leq n$ に対し $f(x_i)=y_i$ を満たすものを $O(n(\log n)^2)$ 時間で求める．
+
+## アルゴリズム
+
+$g(x)=\prod_{i=0}^{n}(x-x_i)$ とすれば Lagrange 補間より以下が成り立つ．
+
+$$f(x)=\sum_{i=0}^{n}\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)$ 時間で計算できる．
+
+$g(x)$ の計算過程はそれ以降にも流用できる．

docs/fps/prefix-sum-of-polynomial.md

@@ -1 +1,30 @@
-多項式 $f$ について $g(n)=\sum_{i=0}^{n-1}f(i)$ を満たす多項式 $g$ を求める．
+多項式 $f$ について $g(n)=\sum_{i=0}^{n}f(i)$ を満たす多項式 $g$ を求める．
+
+$d=\deg f$ として $O(d\log d)$ 時間．
+
+## アルゴリズム
+
+Bernoulli 数 $B_i=\left[\frac{x^i}{i!}\right]\dfrac{xe^x}{e^x-1}$ を用いると以下が成り立つ．
+
+$$\begin{align*}
+\sum_{i=1}^{n}i^k
+&=k![x^k]\sum_{i=1}^{n}e^{ix}\\
+&=k![x^k]\frac{(e^{nx}-1)e^x}{e^x-1}\\
+&=k![x^k]\frac{e^{nx}-1}{x}\cdot\frac{xe^x}{e^x-1}\\
+&=k!\sum_{i=0}^{k}\frac{n^{i+1}}{(i+1)!}\cdot\frac{B_{k-i}}{(k-i)!}
+\end{align*}$$
+
+これは Faulhaber の公式として知られている．
+
+$f(x)=\sum_{i=0}^{d-1}f_ix^i$ とすると，上の結果から $1\leq l\leq d$ に対して
+
+$$\begin{align*}
+[x^l]g(x)
+&=\frac{1}{l!}\sum_{k=l-1}^{d-1}f_kk!\cdot\frac{B_{k-(l-1)}}{(k-(l-1))!}\\
+&=\frac{1}{l!}[x^{d-l}]\left(\sum_{i=0}^{d-1}f_ii!\cdot x^{d-1-i}\right)\frac{xe^x}{e^x-1}
+\end{align*}$$
+
+が成り立つ．また $[x^0]g(x)=f_0$ であるから，結局 $g$ は $O(d\log d)$ 時間で計算できる．
+
+- https://codeforces.com/blog/entry/98563
+- https://www.codechef.com/problems/SERSUM

docs/fps/product-of-polynomials.md

@@ -0,0 +1,7 @@
+多項式 $f_1,f_2,\dots,f_n$ の総積を，次数の総和を $n$ として $O(n(\log n)^2)$ 時間で計算する．
+
+分割統治すればよい．
+
+## FFT 回数削減
+
+それぞれを最低限の長さで FFT し，ダブリングをしていけばよい．

fps/fps-rational.hpp

@@ -0,0 +1,37 @@
+#pragma once
+#include "fps/formal-power-series.hpp"
+
+template <class mint>
+struct FPSRational {
+  using F = FormalPowerSeries<mint>;
+  using R = FPSRational;
+  F num, den;
+  R& operator+=(const R& r) {
+    num *= r.den;
+    num += den * r.num;
+    den *= r.den;
+    return *this;
+  }
+  R& operator-=(const R& r) {
+    num *= r.den;
+    num -= den * r.num;
+    den *= r.den;
+    return *this;
+  }
+  R& operator*=(const R& r) {
+    num *= r.num;
+    den *= r.den;
+    return *this;
+  }
+  R& operator/=(const R& r) {
+    num *= r.den;
+    den *= r.num;
+    return *this;
+  }
+  R operator+(const R& r) const { return R(*this) += r; }
+  R operator-(const R& r) const { return R(*this) -= r; }
+  R operator*(const R& r) const { return R(*this) *= r; }
+  R operator/(const R& r) const { return R(*this) /= r; }
+  R inv() const { return {den, num}; }
+  F approx(int deg) const { return (den * num.inv(deg)).pre(deg); }
+};

fps/polynomial-interpolation.hpp

@@ -0,0 +1,24 @@
+#pragma once
+#include "fps/formal-power-series.hpp"
+
+template <class mint>
+FormalPowerSeries<mint> PolynomialInterpolation(const vector<mint>& x, const vector<mint>& y) {
+  using fps = FormalPowerSeries<mint>;
+  assert(x.size() == y.size());
+  int n = x.size();
+  if (n == 0) return {y[0]};
+  vector<fps> prod(2 * n);
+  for (int i = 0; i < n; i++) prod[i + n] = {-x[i], 1};
+  for (int i = n - 1; i > 0; i--) prod[i] = prod[i * 2] * prod[i * 2 + 1];
+  vector<fps> fs(2 * n);
+  fs[1] = prod[1].diff();
+  for (int i = 2; i < 2 * n; i++) fs[i] = fs[i / 2] % prod[i];
+  for (int i = n; i < n * 2; i++) fs[i] = {y[i - n] / fs[i][0]};
+  for (int i = n - 1; i > 0; i--) fs[i] = fs[(i << 1) | 0] * prod[(i << 1) | 1] + fs[(i << 1) | 1] * prod[(i << 1) | 0];
+  return fs[1];
+}
+
+/**
+ * @brief 多項式補間
+ * @docs docs/fps/polynomial-interpolation.md
+ */

fps/prefix-sum-of-polynomial.hpp

@@ -1,21 +1,27 @@
 #pragma once
+#include "modint/factorial.hpp"
 #include "fps/formal-power-series.hpp"
 
-// sum_{i=0}^{n-1}r^i*poly(i)
-// f[i]=poly(i)
+// g(n)=sum_{i=0}^{n}f(i)
 template <class mint>
-mint SumOfExpPoly(long long n, mint r, vector<mint>& f) {
+FormalPowerSeries<mint> PrefixSumOfPolynomial(FormalPowerSeries<mint> f) {
+  if (f.empty()) return {};
+  using fact = Factorial<mint>;
+  mint c = f[0];
   int d = f.size();
-  mint p = 1;
-  vector<mint> g(d + 1, 0);
-  for (int i = 0; i < d; i++) {
-    g[i + 1] = g[i] + f[i] * p;
-    p *= r;
-  }
-
-  return s;
+  fact::reserve(d);
+  for (int i = 0; i < d; i++) f[i] *= fact::fact(i);
+  reverse(f.begin(), f.end());
+  FormalPowerSeries<mint> g(d);
+  for (int i = 0; i < d; i++) g[i] = fact::fact_inv(i + 1) * (i & 1 ? -1 : 1);
+  f *= g.inv();
+  f.resize(d);
+  f.push_back(c);
+  reverse(f.begin(), f.end());
+  for (int i = 1; i <= d; i++) f[i] *= fact::fact_inv(i);
+  return f;
 }
 /**
- * @brief Prefix Sum of Polynomial
+ * @brief 多項式の Prefix Sum
  * @docs docs/fps/prefix-sum-of-polynomial.md
  */