7K4B blog

猫でも分かる何か

形式的冪級数 LIB::FormalPowerSeries

使い方

 f=9+8z+7z^2 g=9+8z^2+7z^3 に対して
 fg=81+72z+135z^2+127z^3+112z^4+49z^5 を計算する

#include <bits/stdc++.h>
#include <atcoder/all>
using namespace std;
using namespace atcoder;
const long long mint1=1000000007;
const long long mint9=998244353; 

using mint=static_modint<mint9>;
using fps=FormalPowerSeries<mint>;
using sfps=vector<pair<int,mint>>;

int main(void){
	fps f={9,8,7},g1={9,0,8,7};
	f.resize(6);
	fps ans1=f*g1; // 畳み込みを使った演算
	cout<<ans1[5].val()<<endl;
	sfps g2={{0,9},{1,0},{2,8},{3,7}};
	fps ans2=f*g2; // 疎な多項式の演算
	cout<<ans2[5].val()<<endl;
	return 0;
}
49
49

3行の using を書いてからメイン処理を書く。fは3次式で積が5次式なので、fを5次にリサイズする必要がある(リサイズしなかった場合、3次以上の情報は捨てられる)。畳み込みを使う場合はfps型、疎な(項の少ない)多項式にはsfps型を使うと高速に計算できる。なお atcoder ライブラリを利用するので<atcoder/all>が必要。

逆元を使った高速化テク

 \displaystyle (1-x) \sum_{k=0}^{\infty} x^k = 1 より  \displaystyle \sum_{k=0}^{\infty} x^k=\frac{1}{1-x}
また  \displaystyle \sum_{k=0}^{\infty} f^k = \frac{1}{1-f} (ただし  f は定数項を持たない多項式

 \displaystyle f= \sum_{k=1}^{K-1} x^k = \frac{x-x^K}{1-x}, g= \sum_{k=1}^{\infty} x^k = \frac{x}{1-x} とすると
 h=f+fgf+fgfgf+...=\frac{f}{1-gf}=\frac{x-x^2-x^K+x^{K+1}}{1+2x+x^{K+1}} となり  [x^N]h が答え
https://atcoder.jp/contests/tdpc/submissions/45972918

メンバ関数

  • F &operator*=(const F &g): 畳み込みで高速化したn次式とm次式の積  O(k \log k) (k=n+m)
  • F &operator*=(vector<pair<int,T>> g): 疎な多項式の積  O(nm)(n,mは項数)
  • F &operator/=(const F &g): 畳み込みで高速化したn次式とm次式の商  O(k \log k) (k=n+m)
  • F &operator/=(vector<pair<int,T>> g): 疎な多項式の商  O(nm)(n,mは項数)
#include <bits/stdc++.h>
#include <atcoder/all>
using namespace std;
using namespace atcoder;

namespace LIB{
#define rep2(i,m,n) for(int i=(m);i<(n);++i)
#define rep1(i,n) rep2(i,0,n)
#define drep2(i,m,n) for(int i=(m)-1;i>=(n);--i)
#define drep1(i,n) drep2(i,n,0)

	template<class T>struct FormalPowerSeries:vector<T>{
		using vector<T>::vector;
		using vector<T>::operator=;
		using F=FormalPowerSeries;

		F operator-() const {
			F res(*this);
			for (auto &e : res) e = -e;
			return res;
		}
		F &operator*=(const T &g) {
			for (auto &e : *this) e *= g;
			return *this;
		}
		F &operator/=(const T &g) {
			assert(g != T(0));
			*this *= g.inv();
			return *this;
		}
		F &operator+=(const F &g) {
			int n = (*this).size(), m = g.size();
			rep1(i, min(n, m)) (*this)[i] += g[i];
			return *this;
		}
		F &operator-=(const F &g) {
			int n = (*this).size(), m = g.size();
			rep1(i, min(n, m)) (*this)[i] -= g[i];
			return *this;
		}
		F &operator<<=(const int d) {
			int n = (*this).size();
			(*this).insert((*this).begin(), d, 0);
			(*this).resize(n);
			return *this;
		}
		F &operator>>=(const int d) {
			int n = (*this).size();
			(*this).erase((*this).begin(), (*this).begin() + min(n, d));
			(*this).resize(n);
			return *this;
		}
		F inv(int d = -1) const {
			int n = (*this).size();
			assert(n != 0 && (*this)[0] != 0);
			if (d == -1) d = n;
			assert(d > 0);
			F res{(*this)[0].inv()};
			while (res.size() < d) {
				int m = size(res);
				F f(begin(*this), begin(*this) + min(n, 2*m));
				F r(res);
				f.resize(2*m), internal::butterfly(f);
				r.resize(2*m), internal::butterfly(r);
				rep1(i, 2*m) f[i] *= r[i];
				internal::butterfly_inv(f);
				f.erase(f.begin(), f.begin() + m);
				f.resize(2*m), internal::butterfly(f);
				rep1(i, 2*m) f[i] *= r[i];
				internal::butterfly_inv(f);
				T iz = T(2*m).inv(); iz *= -iz;
				rep1(i, m) f[i] *= iz;
				res.insert(res.end(), f.begin(), f.begin() + m);
			}
			return {res.begin(), res.begin() + d};
		}
		F &operator*=(const F &g) {
			int n = (int)(*this).size();
			*this = convolution(*this, g);
			(*this).resize(n);
			return *this;
		}
		F &operator/=(const F &g) {
		   int n = (*this).size();
		   *this = convolution(*this, g.inv(n));
		   (*this).resize(n);
		   return *this;
		}
		F &operator*=(vector<pair<int, T>> g) {
			int n = (int)(*this).size();
			auto [d, c] = g.front();
			if (d == 0) g.erase(g.begin());
			else c = 0;
			drep1(i, n) {
				(*this)[i] *= c;
				for (auto &[j, b] : g) {
					if (j > i) break;
					(*this)[i] += (*this)[i-j] * b;
				}
			}
			return *this;
		}
		F &operator/=(vector<pair<int, T>> g) {
			int n = (*this).size();
			auto [d, c] = g.front();
			assert(d == 0 && c != T(0));
			T ic = c.inv();
			g.erase(g.begin());
			rep1(i, n) {
				for (auto &[j, b] : g) {
					if (j > i) break;
					(*this)[i] -= (*this)[i-j] * b;
				}
				(*this)[i] *= ic;
			}
			return *this;
		}
		void multiply(const int d, const T c) { 
			int n = (*this).size();
			if (c == T(1)) drep1(i, n-d) (*this)[i+d] += (*this)[i];
			else if (c == T(-1)) drep1(i, n-d) (*this)[i+d] -= (*this)[i];
			else drep1(i, n-d) (*this)[i+d] += (*this)[i] * c;
		}
		void divide(const int d, const T c) {
			int n = (*this).size();
			if (c == T(1)) rep1(i, n-d) (*this)[i+d] -= (*this)[i];
			else if (c == T(-1)) rep1(i, n-d) (*this)[i+d] += (*this)[i];
			else rep1(i, n-d) (*this)[i+d] -= (*this)[i] * c;
		}
		T eval(const T &a) const {
			T x(1), res(0);
			for (auto e : *this) res += e * x, x *= a;
			return res;
		}
		F operator*(const T &g) const { return F(*this) *= g; }
		F operator/(const T &g) const { return F(*this) /= g; }
		F operator+(const F &g) const { return F(*this) += g; }
		F operator-(const F &g) const { return F(*this) -= g; }
		F operator<<(const int d) const { return F(*this) <<= d; }
		F operator>>(const int d) const { return F(*this) >>= d; }
		F operator*(const F &g) const { return F(*this) *= g; }
		F operator/(const F &g) const { return F(*this) /= g; }
		F operator*(vector<pair<int, T>> g) const { return F(*this) *= g; }
		F operator/(vector<pair<int, T>> g) const { return F(*this) /= g; }
	};

#undef rep2
#undef rep1
#undef drep2
#undef drep1
}