数学记录

数学记录(无码版)

多项式全家桶

多项式

定义

形如 \(\sum a_nx^n\) 的求和式若为有限项相加，那么称作多项式 记作\(f(x)=\sum_{n=0}^m a_nx^n\)

对于如 \(\sum_{n=0}^{\infty}a_nx^n\) 有一个非负整数次幂乘以一个常数，称作幂级数

运算

加减运算

两个多项式 \(F(x)=\sum_{n\geq0}a_nx^n\) 与 \(G(x)=\sum_{n\geq0}b_nx^n\) 有

\[F(x)\pm G(x)=\sum_{n\geq0}(a_n\pm b_n)x^n \]

卷积（乘法）

两个多项式 \(F(x)=\sum_{n\geq0}a_nx^n\) 与 \(G(x)=\sum_{n\geq0}b_nx^n\) 有

\[F(x)\times G(x)=\sum_{i\geq0}a_ix^i \sum_{j\geq0}b_jx^j\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ =\sum_{n\geq0}x^n\sum_{i\geq0}a_ib_{n-i} \]

FFT

单位根

对于 \(n\) 次单位根为 \(\omega_n\) 表示在复平面上以原点为圆心的半径为1的园，并从实部正半轴开始逆时针将园 \(n\) 等分后的原点与第一个 \(n\) 等分点构成的复数。

在弧度制下，该角弧度为 \(\frac{2\pi}{n}\) ，那么

\[\omega_n^k=\cos (k\times \frac{2\pi}{n}) +i\sin (k\times \frac{2\pi}{n}) \]

那么可以推导出：

\[(\omega_n^n)=(\omega_n^0)=1\\ \omega_{2n}^{2k}=\cos(2k\times\frac{2\pi}{2n})+i\sin(2k\times\frac{2\pi}{2n})=\omega_n^k\\ \omega_n^{\frac{n}{2}}=\cos(\frac{n}{2}\times\frac{2\pi}{n})+i\sin(\frac{n}{2}\times\frac{2\pi}{n})=\cos(\pi)+i\sin(\pi)=-1\\ \omega_n^{k+\frac n2}可以理解为在选择\frac{2k\pi}n 的基础上继续旋转\frac{2\pi \frac{n}2}{n}=\pi,即多旋转180^\circ,\\从向量的角度看，可以感性理解为其相乘，则\omega_n^{k+\frac n2}=-\omega_n^{k} \]

\[(\omega_n^k)^2=(\cos (k\times \frac{2\pi}{n}) +i\sin (k\times \frac{2\pi}{n}))\times(\cos (k\times \frac{2\pi}{n}) +i\sin (k\times \frac{2\pi}{n}))\\ =\cos^2(k\times \frac{2\pi}{n})-\sin^2(k\times \frac{2\pi}{n}) +i2\cos(k\times \frac{2\pi}{n})\sin(k\times \frac{2\pi}{n})\\ \\\because \begin{bmatrix} \cos \alpha&-\sin\alpha\\ \sin \alpha&\cos\alpha\\ \end{bmatrix} \begin{bmatrix} \cos \alpha\\ \sin \alpha \end{bmatrix}= \begin{bmatrix} \cos^2\alpha-\sin^2\alpha\\ 2\cos\alpha\sin\alpha\\ \end{bmatrix}\\ \therefore \cos 2\alpha=\cos^2\alpha-\sin^2\alpha\\ \sin 2\alpha=2\cos\alpha\sin\alpha\\ \therefore(\omega_n^k)^2=\cos(2k\times\frac{2\pi}{n})+i\sin(2k\times\frac{2\pi}{n})=\omega_n^{2k} \]

\[可以感性理解一下(\omega_n^i)(\omega_n^j)=(\omega_n^{i+j})，因为其可以看作向量旋转 \]

DFT

考虑将这 \(n\) 个 \(n\) 次单位根带入一个 \(n\) 次多项式，获得其点值表示法，复杂度肯定是 \(O(n^2)\) ，但是如果对下标进行奇偶性分离：

\[A(x)=\sum_{i=0}^{n-1} a_ix^i=a_0+a_1x^1+a_2x^2+\dots +a_{n-2}x^{n-2}+a_{n-1}x^{n-1}\\ =a_0+a_2x^2+\dots+a_{n-2}x^{n-2} \\+a_1x^1+a_3x^3+\dots+a_{n-1}x^{n-1}\\ 令A_1(x)=a_0+a_2x^1+\dots+a_{n-2}x^{\frac n2-1}\\A_2(x)=a_1+a_3x^1+\dots+a_{n-1}x^{\frac n2-1}\\ 那么A(x)=A_1(x^2)+xA_2(x^2)\\ 将其中一个单位根 \omega_n^k \ (k<\frac n2) 带入：\\ A(\omega_n^k)=A_1((\omega_n^k)^2)+\omega_n^kA_2((\omega_n^k)^2)\\ =A_1(\omega_n^2k)+\omega_n^kA_2(\omega_n^2k)\\ 若将\omega_n^{k+\frac n2}带入：\\ A`(-\omega_n^k)=A_1(\omega_n^2k)-\omega_n^kA_2(\omega_n^2k)\\ \]

我们可以发现，\(A\) 与 \(A`\) 只有一个常数项不同，所以在求 \(A\) 时可以 \(O(1)\) 算出第二个，所以 \(k\) 的取值范围时 \([0,\frac n 2-1]\) ，然后发现 \(A_1,A_2,A\) 的性质相同，所以可以递归转化，复杂度为 \(O(n\log_2 n)\)。

IDFT

现在已知 \((y_0,y_1,y_2,\dots,y_{n-1})\) 是 \((a_0,a_1,a_2,\dots,a_{n-1})\) 的点值表示法。

假设有个向量 \((c_0,c_1,c_2,\dots,c_{n-1})\) ，即多项式 \(B(x)=\sum\limits_{i=0}^{n-1}y_ix^i\) 在 \((\omega_n^0,\omega_n^{-1},\omega_n^{-2},\dots,\omega_n^{-(n-1)})\) 上的值。

那么：

\[\begin{aligned} c_k&=\sum_{i=0}^{n-1}y_i(\omega_n^{-k})^i\\ &=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_n^i)^j)(\omega_n^{-k})^i\\ &=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_n^i)^j(\omega_n^{-k})^i)\\ &=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j(\omega_n^i)^j(\omega_n^{-k})^i\\ &=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j(\omega_n^j)^i(\omega_n^{-k})^i\\ &=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j(\omega_n^{j-k})^i\\ &=\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}(\omega_n^{j-k})^i\\ \end{aligned} \]

\[设S(x,n)=\sum_{i=0}^{n-1}x^i\\ 当x=\omega_n^k时：\\ S(\omega_n^k,n)=1+(\omega_n^k)+(\omega_n^k)^2+\dots+(\omega_n^k)^{n-1}\\ (\omega_n^k)S(\omega_n^k,n)=(\omega_n^k)+(\omega_n^k)^2+(\omega_n^k)^3+\dots+(\omega_n^k)^n\\ 两式相减：\\ (\omega_n^k-1)S(\omega_n^k,n)=(\omega_n^k)^n-1\\ S(\omega_n^k,n)=\frac{(\omega_n^k)^n-1}{\omega_n^k-1}\\ \therefore S(\omega_n^k,n)=\frac{1-1}{\omega_n^k-1}\\ 但是在k=0时S(\omega_n^k,n)=n\\ 所以在原式c_k=\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}(\omega_n^{j-k})^i\\ 当i\neq k时，后一个\sum中值为0\\ 当i= k时，后一个\sum中值为n\\ 故c_k=na_k\\ a_k=\frac{c_k}{n} \]

总结

那么两个多项式相乘可以先将其分别转化为点值表示法，然后对于位置相乘得到结果的点值表示法，然后在跑一遍那个将多项式转化为点值的函数，将点值转化为多项式，最后依据 \(IDFT\) 除以即可。

#include <bits/stdc++.h>
#define pii pair<int,int>
#define ll long long
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
bool Mst;
const int Max=3e6+10;
const int mod=998244353;
const int inf=1e9+10;
const double PI=acos(-1.0);


inline int read(){
	int res=0,v=1;
	char c=getchar();
	while(c<'0'||c>'9'){v=(c=='-'?-1:1);c=getchar();}
	while(c>='0'&&c<='9'){res=(res<<3)+(res<<1)+(c^48);c=getchar();}
	return res*v;
}

struct comple{
	double a,b;
	comple(double a=0,double b=0):
		a(a),b(b){;}
	comple operator +(const comple &x)const{return comple(a+x.a,b+x.b);}
	comple operator -(const comple &x)const{return comple(a-x.a,b-x.b);}
	comple operator *(const comple &x)const{return comple(a*x.a-b*x.b,a*x.b+b*x.a);}
}; 

void Fast(int now,comple *a,int tag){
	if(now==1)return;
	comple a1[(now>>1)+10],a2[(now>>1)+10];
	for(int i=0;i<(now>>1);++i){
		a1[i]=a[i<<1];
		a2[i]=a[i<<1|1];
	}
	Fast(now>>1,a1,tag);
	Fast(now>>1,a2,tag);
	comple W=comple(cos(2.0*PI/now),tag*(sin(2.0*PI/now))),w=comple(1,0);
	for(int i=0;i<(now>>1);++i,w=w*W){
		comple z=w*a2[i];
		a[i]=a1[i]+z;
		a[i+(now>>1)]=a1[i]-z;
	}
}

comple a[Max],b[Max];

bool Med;
signed main(){
	int n=read();
	int m=read();
	for(int i=0;i<=n;++i)a[i].a=raed();
	for(int i=0;i<=m;++i)b[i].a=read();
	int pos=1;
	while(pos<=n+m)pos<<=1;
	Fast(pos,a,1);
	Fast(pos,b,1);
	for(int i=0;i<=pos;++i)a[i]=a[i]*b[i];
	Fast(pos,a,-1);
	for(int i=0;i<=n+m;++i)cout << (int)(a[i].a/pos+0.5)<< ' ';
	

	cerr<< "Time: "<<clock()/1000.0 << "s\n";
	cerr<< "Memory: " << (&Mst-&Med)/1000000.0 << "MB\n";
	return 0;
}

实际上对于上述算法还可以进行优化，就是不递归的写法。

因为我们发现一个数的位置是其元下标二进制翻转后的位置，所有可以预处理然后枚举长度合并。

#include <bits/stdc++.h>
#define pii pair<int,int>
#define ll long long
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
bool Mst;
const int Max=3e6+10;
const int mod=998244353;
const int inf=1e9+10;
const double PI=acos(-1.0);

inline int read(){
	int res=0,v=1;
	char c=getchar();
	while(c<'0'||c>'9'){v=(c=='-'?-1:1);c=getchar();}
	while(c>='0'&&c<='9'){res=(res<<3)+(res<<1)+(c^48);c=getchar();}
	return res*v;
}

struct comple {
	
	double a,b;
	comple (double a=0,double b=0):
		a(a),b(b){;}
	comple operator +(const comple &x)const{return comple(a+x.a,b+x.b);}
	comple operator -(const comple &x)const{return comple(a-x.a,b-x.b);}
	comple operator *(const comple &x)const{return comple(a*x.a-b*x.b,a*x.b+b*x.a);}
};

int rev[Max];

int limit;

void Fast(comple *a,int tag){
	for(int i=0;i<limit;++i)if(i<rev[i])swap(a[i],a[rev[i]]);
	for(int mid=1;mid<limit;mid<<=1){
		comple W=comple(cos(PI/mid),tag*sin(PI/mid));
		for(int R=mid<<1,j=0;j<limit;j+=R){
			comple w=comple(1,0);
			for(int l=0;l<mid;++l,w=w*W){
				comple x=a[l+j],y=a[l+j+mid]*w;
				a[l+j]=x+y;
				a[l+j+mid]=x-y;
			}
		}
	}
}

comple a[Max],b[Max];

bool Med;
signed main(){
	int n=read();
	int m=read();
	for(int i=0;i<=n;++i)a[i].a=read();
	for(int i=0;i<=m;++i)b[i].a=read();
	limit=1;
	int len=0;
	while(limit<=n+m)limit<<=1,++len;
	for(int i=0;i<=limit;++i)rev[i]=((rev[i>>1]>>1)|((i&1)<<(len-1)));
	Fast(a,1);
	Fast(b,1);
	for(int i=0;i<=limit;++i)a[i]=a[i]*b[i];
	Fast(a,-1);
	for(int i=0;i<=n+m;++i)cout << (int)(a[i].a/limit+0.5) << ' ';
	cout << "\n";

	cerr<< "Time: "<<clock()/1000.0 << "s\n";
	cerr<< "Memory: " << (&Mst-&Med)/1000000.0 << "MB\n";
	return 0;
}

NTT

原根

将最小的满足 \(a^l\equiv1\pmod p\) 的 \(l\) 称作 \(a\) 在模 \(p\) 意义下的阶。

如果 \(g\) 在模 \(p\) 意义下的阶是 \(\varphi(p)\)，那么称 \(g\) 是 \(p\) 的原根。

求原根可以这样做：如果 \(\gcd(g,p)=1\)，设 \(p_1,p_2 \dots p_k\) 为 \(p\) 的所有质因子，如果 \(g\) 是 \(p\) 的原根就满足对于任意 \(i\in[1,k],g^{\frac{\varphi(p)}{p_i}}\not\equiv1\pmod p\)

可以发现原根有一个非常优秀的性质就是对于在 \(1\) 到 \(\varphi(p)\) 中的任意一个幂次都是不一样的，而且是以 \(\varphi(p)\) 为长度的一个循环节，这与 \(\omega\) 的性质很像。而且当 \(p\) 为质数时，如果 \(\varphi(p)\) 有很多 \(2\) 这个因子就十分优秀，像 \(998244353=2^{23}\times 7\times 17+1\) 它的原根是 \(3\) 。

那么现在我们可以用 \(g^{\frac{\varphi(p)}{n}}\) 代替 \(n\) 次单位根就行。

#include <bits/stdc++.h>
#define pii pair<int,int>
#define ll long long
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
bool Mst;
const int Max=3e6+10;
const int mod=998244353;
const int inf=1e9+10;
const int g=3;
const int invg=332748118;

template <int mod>
struct modint{

	int val;

	static int norm(const int &x){return x<0?x+mod:x;}
	
	modint inv()const{
		int a=val,b=mod,u=1,v=0,t;
		while(b>0)t=a/b,swap(a-=t*b,b),swap(u-=t*v,v);
		return modint(u);
	}
	modint():val(0){}
	modint(const int &m):val(norm(m)){}
	modint(const ll &m):val(norm(m%mod)){}
	modint operator -()const{return modint(norm(-val));}
	bool operator ==(const modint &x){return val==x.val;}
	bool operator !=(const modint &x){return val!=x.val;}
	bool operator <=(const modint &x){return val<=x.val;}
	bool operator >=(const modint &x){return val>=x.val;}
	bool operator >(const modint &x){return val>x.val;}
	bool operator <(const modint &x){return val<x.val;}
	modint& operator *=(const modint &x){return val=static_cast<int>(1ll*val*x.val%mod),*this;}
	modint& operator +=(const modint &x){return val=(1ll*val+x.val)%mod,*this;}
	modint& operator -=(const modint &x){return val=norm(1ll*val-x.val),*this;}
	modint& operator >>=(const modint &x){return val>>=x.val,*this;}
	modint& operator <<=(const modint &x){return val<<=x.val,*this;}
	modint& operator ^=(const modint &x){return val^=x.val,*this;}
	modint operator >>(const modint &x){return modint(*this)>>=x;}
	modint operator <<(const modint &x){return modint(*this)<<=x;}
	modint& operator /=(const modint &x){return *this*=x.inv();}
	modint operator +(const modint &x){return modint(*this)+=x;}
	modint operator -(const modint &x){return modint(*this)-=x;}
	modint operator *(const modint &x){return modint(*this)*=x;}
	modint operator /(const modint &x){return modint(*this)/=x;}
	modint operator ^(const modint &x){return modint(*this)^=x;}
	friend std::ostream& operator<<(std::ostream& os,const modint &a){return os<<a.val;}
	friend std::istream& operator>>(std::istream& is,modint &a){return is>>a.val;}

};

typedef modint<998244353>m98;

inline int read(){
	int res=0,v=1;
	char c=getchar();
	while(c<'0'||c>'9'){v=(c=='-'?-1:1);c=getchar();}
	while(c>='0'&&c<='9'){res=(res<<3)+(res<<1)+(c^48);c=getchar();}
	return res*v;
}

m98 ksm(m98 a,int b){
	m98 ans=1;
	for(;b;b>>=1){
		if(b&1)ans=ans*a;
		a=a*a;
	}
	return ans;
}

m98 a[Max],b[Max];

int rev[Max],limit;


void Fast(m98 *a,int tag){
	for(int i=0;i<limit;++i)if(i<rev[i])swap(a[i],a[rev[i]]);
	for(int mid=1;mid<limit;mid<<=1){
		m98 W=ksm(tag==1?g:invg,(mod-1)/(mid<<1));
		for(int R=mid<<1,j=0;j<limit;j+=R){
			m98 w=1;
			for(int p=0;p<mid;++p,w=w*W){
				m98 x=a[p+j],y=w*a[p+j+mid];
				a[p+j]=x+y;
				a[p+j+mid]=x-y;
			}
		}
	}
	if(tag==-1){
		m98 inv=m98(limit).inv();
		for(int i=0;i<limit;++i)a[i]=a[i]*inv;
	}
}


bool Med;
signed main(){
	int n=read();
	int m=raed();
	for(int i=0;i<=n;++i)a[i].val=read();
	for(int i=0;i<=m;++i)b[i].val=raed();
	limit=1;
	int len=0;
	while(limit<=n+m)limit<<=1,++len;
	for(int i=0;i<=limit;++i)rev[i]=((rev[i>>1]>>1)|((i&1)<<(len-1)));
	Fast(a,1);
	Fast(b,1);
	for(int i=0;i<=limit;++i)a[i]=a[i]*b[i];
	Fast(a,-1);
	for(int i=0;i<=n+m;++i)cout << a[i] << ' ';
	cout << "\n";
	

	cerr<< "Time: "<<clock()/1000.0 << "s\n";
	cerr<< "Memory: " << (&Mst-&Med)/1000000.0 << "MB\n";
	return 0;
}
/*

*/

多项式牛顿迭代

给定多项式 \(g(x)\) 求在模 \(x^n\) 下的 \(f(x)\) 使其满足：

\[g(f(x))\equiv0\pmod {x^n} \]

通过倍增和泰勒展开(不会)得到：

\[f(x)\equiv f_0(x)-\frac{g(f_0(x))}{g'(f_0(x))}\pmod{x^n} \]

多项式求逆

设给定的多项式是 \(h(x)\)，那么有：

\[g(f(x))=\frac{1}{f(x)}-h(x)\equiv 0\pmod{x^n} \]

于是有：

\[f(x)\equiv f_0(x)-\frac{\frac{1}{f_0(x)}-h(x)}{-\frac{1}{f_0^2(x)}}\pmod{x^n}\\ \equiv 2f_0(x)-f_0^2(x)h(x)\pmod{x^n}\\ \]

而当 \(n=1\) 时，\([x^0]f_0(x)=[x^0]h(x)^{-1}\)

多项式开根

设给定的多项式是 \(h(x)\)，那么有：

\[g(f(x))\equiv f^2(x)-h(x)\pmod{x^n}\\ \]

于是有：

\[f(x)\equiv f_0(x)-\frac{f^2_0(x)-h(x)}{2f_0(x)}\pmod{x^n}\\ \equiv\frac{f_0^2(x)+h(x)}{2f_0(x)} \]

多项式对数函数

假设我们要求 \(G(x)\equiv\ln(F(x))\pmod{x^n}\) ，设函数 \(f(x)=\ln(x)\)，那么我们要求的是：

\[G(x)\equiv f(F(x))\pmod{x^n}\\ G'(x)=f(F(x))'\\ G'(x)=f'(F(x))F'(x)\\ \]

注意到 \(\ln(x)\) 的导数是 \(\frac{1}{x}\)，所以有：

\[G'\equiv \frac{F'(x)}{F(x)} \pmod{x^n} \]

因为我们只关心余数，所以可以求逆就行。

​

多项式指数函数

设给定的多项式是 \(h(x)\)，那么有：

\[g(f(x))\equiv\ln(f(x))-h(x)\pmod{x^n} \]

于是有：

\[f(x)\equiv f_0(x)-\frac{\ln(f_0(x))-h(x)}{\frac{1}{f_0(x)}}\pmod{x^n}\\ \equiv(f_0(x))(1-\ln(f_0(x))+h(x))\pmod{x^n} \]

多项式快速幂

设给定的多项式是 \(h(x)\)，那么求：

\[G(x)\equiv h(x)^k \pmod{x^n}\\ \ln(G(x))\equiv k\ln(h(x))\pmod{x^n}\\ G(x)\equiv e^{k\ln(h(x))}\pmod{x^n} \]

注意到要求 \(\ln\) 所以要保证第一位位 \(1\) ，所以对于部分不合法多项式应该左移 \(z\) 位到第一位不为 \(0\)，如果是别的数要先除以 \(a_0\)，但是最后要每一项乘以 \(a_0^{k\mod \varphi(p)}\) ，然后还要乘以一个 \(x^{zk\mod p}\)

Code

#include <bits/stdc++.h>
#define pii pair<int,int>
#define ll long long
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
bool Mst;
const int Max=6e5+10;
const int mod=998244353;
const int inf=1e9+10;

template <int mod>
struct modint{

	int val;

	static int norm(const int &x){return x<0?x+mod:x;}

	modint inv()const{
		int a=val,b=mod,u=1,v=0,t;
		while(b>0)t=a/b,swap(a-=t*b,b),swap(u-=t*v,v);
		return modint(u);
	}
	
	
	modint():val(0){}
	modint(const int &m):val(norm(m)){}
	modint(const ll &m):val(norm(m%mod)){}
	modint operator -()const{return modint(norm(-val));}
	bool operator ==(const modint &x){return val==x.val;}
	bool operator !=(const modint &x){return val!=x.val;}
	bool operator <=(const modint &x){return val<=x.val;}
	bool operator >=(const modint &x){return val>=x.val;}
	bool operator >(const modint &x){return val>x.val;}
	bool operator <(const modint &x){return val<x.val;}
	modint& operator *=(const modint &x){return val=static_cast<int>(1ll*val*x.val%mod),*this;}
	modint& operator <<=(const modint &x){return val=(1ll*val<<x.val)%mod,*this;}
	modint& operator +=(const modint &x){return val=(1ll*val+x.val)%mod,*this;}
	modint& operator -=(const modint &x){return val=norm(1ll*val-x.val),*this;}
	modint& operator >>=(const modint &x){return val>>=x.val,*this;}
	modint& operator ^=(const modint &x){return val^=x.val,*this;}
	modint operator >>(const modint &x){return modint(*this)>>=x;}
	modint operator <<(const modint &x){return modint(*this)<<=x;}
	modint& operator /=(const modint &x){return *this*=x.inv();}
	modint operator +(const modint &x){return modint(*this)+=x;}
	modint operator -(const modint &x){return modint(*this)-=x;}
	modint operator *(const modint &x){return modint(*this)*=x;}
	modint operator /(const modint &x){return modint(*this)/=x;}
	modint operator ^(const modint &x){return modint(*this)^=x;}
	friend std::ostream& operator<<(std::ostream& os,const modint &a){return os<<a.val;}
	friend std::istream& operator>>(std::istream& is,modint &a){return is>>a.val;}
};

typedef modint<998244353>m98;

inline int read(){
	int res=0,v=1;
	char c=getchar();
	while(c<'0'||c>'9'){v=(c=='-'?-1:1);c=getchar();}
	while(c>='0'&&c<='9'){res=(res<<3)+(res<<1)+(c^48);c=getchar();}
	return res*v;
}

m98 ksm(m98 a,int b){
	m98 ans=1;
	for(;b;b>>=1){
		if(b&1)ans=ans*a;
		a=a*a;
	}
	return ans;
}

m98 inv[Max],frac[Max];

void init(int len) {
	frac[0]=inv[0]=1; 
	for(int i=1;i<=len;++i)frac[i]=frac[i-1]*i;
	inv[len]=frac[len].inv();for(int i=len-1;i;--i)inv[i]=inv[i+1]*(i+1);
	for(int i=1;i<=len;++i)inv[i]*=frac[i-1];
}

int rev[Max];

void init(int lim,int len){for(int i=0;i<lim;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));}

struct Poly{  //导和积都到n，否则到n-1; 
	const int g=3;
	const int invg=332748118;
	
	m98 a[Max];
	void out(int n,int x=0){for(int i=0;i<n;++i)cout << a[i].val << " ";if(!x)cout << "\n";}
	void inte(int n){for(int i=n;i>=0;--i)a[i]=a[i-1]*inv[i];a[0]=0;}
	void ReadIn(int n){for(int i=0;i<n;++i)a[i].val=read();}
	void der(int n){for(int i=0;i<=n;++i)a[i]=a[i+1]*(i+1);}
	m98& operator [](const int x){return a[x];}
	
	void NTT(int tag,int lim){
		for(int i=0;i<lim;++i)if(i<rev[i])swap(a[i],a[rev[i]]);
		for(int mid=1;mid<lim;mid<<=1){
			m98 W=ksm(tag==1?g:invg,(mod-1)/(mid<<1));
			for(int R=mid<<1,j=0;j<lim;j+=R){
				m98 w=1;
				for(int p=0;p<mid;++p,w=w*W){
					m98 x=a[j+p],y=a[j+p+mid]*w;
					a[j+p]=x+y;
					a[j+p+mid]=x-y;
				}
			}
		}
		if(tag==-1){
			m98 inv=m98(lim).inv();
			for(int i=0;i<lim;++i)a[i]*=inv;
		}
	}
	
}a;

void GetInv(Poly &y,Poly &x,int n){    //第二个为参数，第一个为返回值
	static Poly b[2],c;
	int lim=2,len=1,bas=1;
	int pos=0;
	b[pos][0]=x[0].inv();
	while(bas<(n<<1)){
		init(lim,len);
		pos^=1;
		for(int i=0;i<bas;++i)b[pos][i]=b[pos^1][i]<<1;
		for(int i=0;i<=lim;++i)c[i]=x[i];
		for(int i=lim>>1;i<=lim;++i)c[i]=0;
		b[pos^1].NTT(1,lim);c.NTT(1,lim);
		for(int i=0;i<=lim;++i)b[pos^1][i]=b[pos^1][i]*b[pos^1][i]*c[i];
		b[pos^1].NTT(-1,lim);
		for(int i=0;i<bas;++i)b[pos][i]-=b[pos^1][i],b[pos^1][i]=0;
		bas<<=1,lim<<=1,++len;
	}
	for(int i=0;i<n;++i)y[i]=b[pos][i];
	for(int i=0;i<bas;++i)b[0][i]=b[1][i]=c[i]=0;
}

void GetSqrt(Poly &y,Poly &x,int n){
	
	static Poly t1,t2,s;
	s[0]=1;
	int lim=2,len=1,bas=1;
	while(bas<(n<<1)){
		for(int i=0;i<bas;++i)t1[i]=x[i];
		for(int i=0;i<bas;++i)t2[i]=s[i]<<1;
		GetInv(t2,t2,bas);init(lim,len);
		t2.NTT(1,lim);t1.NTT(1,lim);
		for(int i=0;i<=lim;++i)t1[i]=t1[i]*t2[i];
		t1.NTT(-1,lim);
		for(int i=0;i<bas;++i)s[i]=s[i]/2+t1[i];
		bas<<=1,lim<<=1,++len;
	}
	for(int i=0;i<n;++i)y[i]=s[i];
	for(int i=0;i<bas;++i)t1[i]=t2[i]=s[i]=0;
	
}

void GetLn(Poly &y,Poly&x,int n){
	static Poly t1,t2;
	for(int i=0;i<n;++i)t1[i]=t2[i]=x[i];
	t1.der(n);GetInv(t2,t2,n);
	int lim=1,len=0;
	while(lim<n<<1)lim<<=1,++len;
	init(lim,len);
	t1.NTT(1,lim),t2.NTT(1,lim) ;
	for(int i=0;i<lim;++i)t1[i]*=t2[i];
	t1.NTT(-1,lim);t1.inte(n) ;
	for(int i=0;i<n;++i)y[i]=t1[i];
	for(int i=0;i<lim;++i)t1[i]=t2[i]=0;
}


void GetExp(Poly &y,Poly &x,int n){
	static Poly s,tmp;s[0]=1;
	int bas=1,lim=2,len=1;
	while(bas<n<<1){
		for(int i=0;i<bas;++i)tmp[i]=s[i];
		GetLn(tmp,tmp,bas);
		for(int i=0;i<bas;++i)tmp[i]=x[i]-tmp[i];
		tmp[0]+=1;init(lim,len);tmp.NTT(1,lim),s.NTT(1,lim) ;
		for(int i=0;i<lim;++i)s[i]=s[i]*tmp[i];
		s.NTT(-1,lim) ;
		for(int i=bas;i<lim;++i)s[i]=0;
		lim<<=1,bas<<=1,++len;
	}
	for(int i=0;i<n;++i)y[i]=s[i];
	for(itn i=0;i<bas;++i)s[i]=tmp[i]=0;
}


int Si=0;
void expow(Poly &y,Poly &x,int n,m98 k,int k1){
	static Poly s;
	int z=n+1;
	for(int i=0;i<n;++i)if(x[i]!=0){z=i;break;}
	if(z==n+1||(z!=0&&Si)){
		for(int i=0;i<n;++i)y[i]=0;
		return;
	}
	for(itn i=z;i<n;++i)s[i-z]=x[i]/x[z];GetLn(s,s,n-z);
	for(int i=0;i<n-z;++i)s[i]*=k;GetExp(s,s,n-z);
	for(int i=0;i<n-z;++i)s[i]=s[i]*ksm(x[z],k1);
	for(int i=0;i<n;++i)y[i]=0;
	for(int i=0;1ll*z*k.val+i<n;++i)y[i+z*k.val]=s[i];
	for(int i=0;i<n;++i)s[i]=0; 
}


Poly b;

m98 k;
ll k1;

void In(){
	k=k1=0;
	char c=getchar();
	while(c<'0'||c>'9')c=getchar();
	while(c>='0'&&c<='9'){Si=Si|(k*10+(c^48)<k),k=(k<<3)+(k<<1)+(c^48),k1=(k1*10%(mod-1)+(c^48))%(mod-1),c=getchar();}
}

bool Med;
signed main(){
	init(300000);

	cerr<< "Time: "<<clock()/1000.0 << "s\n";
	cerr<< "Memory: " << (&Mst-&Med)/1000000.0 << "MB\n";
	return 0;
}

第二类斯特林数

定义

第二类斯特林数 \(\begin{Bmatrix}n\\m\end{Bmatrix}\) 也记作 \(S(n,m)\)，表示将 \(n\) 个区分的小球放入 \(m\) 个相同的盒子且不允许有空盒子的方案数。

递推公式

\[\begin{Bmatrix}n\\m\end{Bmatrix}=\begin{Bmatrix}n-1\\m-1\end{Bmatrix}+m\begin{Bmatrix}n-1\\m\end{Bmatrix} \]

通过组合意义来解释：

当加入一个球时可以将其单独放入一个空盒子中，方案数为 \(\begin{Bmatrix}n-1\\m-1\end{Bmatrix}\) ，

也可以将其放入一个非空的盒子中，方案数为 \(m\begin{Bmatrix}n-1\\m\end{Bmatrix}\)。

通过加法原理得到递推公式。

通向公式

\[\begin{Bmatrix}n\\m\end{Bmatrix}=\sum_{i=0}^m\frac{(-1)^{m-i}i^n}{i!(m-i)!} \]

可以使用容斥证明该式子：

1. 将 \(n\) 个区分的球放入 \(m\) 个区分的盒子中且允许有空盒子的方案为 \(A_m\)，
2. 将 \(n\) 个区分的球放入 \(m\) 个区分的盒子但不能有空盒子的方案是 \(B_m\)。

那么有 :

\[A_m=m^n=\sum_{i=0}^m\binom{m}{i}B_i \]

通过二项式反演可得：

\[\begin{aligned} B_m&=\sum_{i=0}^m(-1)^{m-i}\binom{m}{i}A_i\\&=\sum_{i=0}^m (-1)^{m-i}\frac{m!i^n}{i!(m-i)!} \end{aligned} \]

因为 \(B_m\) 正好 \(\begin{Bmatrix}n\\m\end{Bmatrix}\) 的 \(i!\) 倍，所以得到通向公式：

\[\begin{Bmatrix}n\\m\end{Bmatrix}=\sum_{i=0}^m\frac{(-1)^{m-i}i^n}{i!(m-i)!} \]

第二类斯特林数·行

发现其通向公式

\[\begin{Bmatrix}n\\i\end{Bmatrix}=\sum\limits_{j=0}^i\frac{(-1)^{i-j}j^n}{j!(i-j)!}\\ =\sum_{j=0}^i\frac{j^n}{j!}.\frac{(-1)^{i-j}}{(i-j)!} \]

这是一个卷积的形式，所以可以 NTT 求。

#include <bits/stdc++.h>
#define pii pair<int,int>
#define ll long long
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
bool Mst;
const int Max=1e6+10;
const int mod=167772161;
const int inf=1e9+10;

template <int mod>
struct modint{

	int val;

	static int norm(const int &x){return x<0?x+mod:x;}

	modint inv()const{
		int a=val,b=mod,u=1,v=0,t;
		while(b>0)t=a/b,swap(a-=t*b,b),swap(u-=t*v,v);
		return modint(u);
	}


	modint():val(0){}
	modint(const int &m):val(norm(m)){}
	modint(const ll &m):val(norm(m%mod)){}
	modint operator -()const{return modint(norm(-val));}
	bool operator ==(const modint &x){return val==x.val;}
	bool operator !=(const modint &x){return val!=x.val;}
	bool operator <=(const modint &x){return val<=x.val;}
	bool operator >=(const modint &x){return val>=x.val;}
	bool operator >(const modint &x){return val>x.val;}
	bool operator <(const modint &x){return val<x.val;}
	modint& operator *=(const modint &x){return val=static_cast<int>(1ll*val*x.val%mod),*this;}
	modint& operator <<=(const modint &x){return val=(1ll*val<<x.val)%mod,*this;}
	modint& operator +=(const modint &x){return val=(1ll*val+x.val)%mod,*this;}
	modint& operator -=(const modint &x){return val=norm(1ll*val-x.val),*this;}
	modint& operator >>=(const modint &x){return val>>=x.val,*this;}
	modint& operator ^=(const modint &x){return val^=x.val,*this;}
	modint operator >>(const modint &x){return modint(*this)>>=x;}
	modint operator <<(const modint &x){return modint(*this)<<=x;}
	modint& operator /=(const modint &x){return *this*=x.inv();}
	modint operator +(const modint &x){return modint(*this)+=x;}
	modint operator -(const modint &x){return modint(*this)-=x;}
	modint operator *(const modint &x){return modint(*this)*=x;}
	modint operator /(const modint &x){return modint(*this)/=x;}
	modint operator ^(const modint &x){return modint(*this)^=x;}
	friend std::ostream& operator<<(std::ostream& os,const modint &a){return os<<a.val;}
	friend std::istream& operator>>(std::istream& is,modint &a){return is>>a.val;}
};

typedef modint<167772161>m16;

inline int read(){
	int res=0,v=1;
	char c=getchar();
	while(c<'0'||c>'9'){v=(c=='-'?-1:1);c=getchar();}
	while(c>='0'&&c<='9'){res=(res<<3)+(res<<1)+(c^48);c=getchar();}
	return res*v;
}

int rev[Max];

void init(int lim,int len){
	for(int i=0;i<lim;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));
}

m16 ksm(m16 a,int b){
	m16 ans=1;
	for(;b;b>>=1){
		if(b&1)ans=ans*a;
		a=a*a;
	}
	return ans;
}

struct Poly{
	
	const m16 g=3;
	const m16 invg=(mod+1)/3;
	m16 a[Max];
	
	m16& operator [](const int x){return a[x];}
	
	void NTT(int tag,int lim){
		for(int i=0;i<lim;++i)if(i<rev[i])swap(a[i],a[rev[i]]);
		for(int mid=1;mid<lim;mid<<=1){
			m16 W=ksm(tag==1?g:invg,(mod-1)/(mid<<1));
			for(int R=mid<<1,j=0;j<lim;j+=R){
				m16 w=1;
				for(int p=0;p<mid;++p,w=w*W){
					m16 x=a[p+j],y=a[p+j+mid]*w;
					a[p+j]=x+y;
					a[p+j+mid]=x-y;
				}
			}
		}
		if(tag==-1){
			m16 inv=m16(lim).inv();
			for(int i=0;i<lim;++i)a[i]*=inv;
		}
	}
	
};


Poly a,b;

m16 frac[Max],inv[Max];

bool Med;
signed main(){
	int n=read();
	inv[0]=frac[0]=1;
	for(int i=1;i<=n;++i)frac[i]=frac[i-1]*i;
	inv[n]=frac[n].inv();
	for(int i=n-1;i>=1;--i)inv[i]=inv[i+1]*(i+1);
	for(int i=0;i<=n;++i){
		a[i]=ksm(i,n)*inv[i];
		b[i]=inv[i]*(i%2==0?1:-1);
	}
	int lim=1,len=0;
	while(lim<=n+n)lim<<=1,++len;
	init(lim,len);
	a.NTT(1,lim),b.NTT(1,lim) ;
	for(int i=0;i<lim;++i)a[i]*=b[i];
	a.NTT(-1,lim) ;
	for(int i=0;i<=n;++i)cotu << a[i] << ' ';
	cotu << "\n";


	cerr<< "Time: "<<clock()/1000.0 << "s\n";
	cerr<< "Memory: " << (&Mst-&Med)/1000000.0 << "MB\n";
	return 0;
}

第二类斯特林数·列

先假设每个盒子也要区分，考虑对于一个盒子计算，其生成函数是 \(G(x)=\sum\limits_{i\geq 1}\frac{x^i}{i!}=(e^x-1)\)，其第 \(i\) 项表示在这个盒子放入了 \(i\) 个球。 这样对于 \(m\) 个盒子，\(F(x)=G(x)^m\)，那么如果有 \(n\) 个球，答案就是 \([x^n]G(x)^m\) 由于目前的盒子是区分的，所以要乘以一个 \(m!\)，故有 \(n\) 个球的答案是 \(m![x^n]G(x)^m\)。这样可以多项式快速幂求，都是发现多项式 \(G(x)\) 的首项是 \(0\) 所以要多项式平移。这里使用的是 EGF 其卷积形式如下 \(\sum\limits_{i\geq 0}x^i\sum\limits_{j=0}^i\frac{1}{i!(i-j)!}a_ib_j=\sum\limits_{i\geq 0} \frac{x^i}{i!}\sum\limits_{j=0}^i\binom{i}{j}a_ib_j\) 就是考虑了 \(a\),\(b\) 合并的顺序，而最后要乘一个 \(i!\) 去消掉 EGF 的 \(i!\) 这个系数。

#include <bits/stdc++.h>
#define pii pair<int,int>
#define ll long long
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
bool Mst;
const int Max=2e6+10;
const int mod=167772161;
const int inf=1e9+10;

template <int mod>
struct modint{

	int val;

	static int norm(const int &x){return x<0?x+mod:x;}

	modint inv()const{
		int a=val,b=mod,u=1,v=0,t;
		while(b>0)t=a/b,swap(a-=t*b,b),swap(u-=t*v,v);
		return modint(u);
	}
	
	
	modint():val(0){}
	modint(const int &m):val(norm(m)){}
	modint(const ll &m):val(norm(m%mod)){}
	modint operator -()const{return modint(norm(-val));}
	bool operator ==(const modint &x){return val==x.val;}
	bool operator !=(const modint &x){return val!=x.val;}
	bool operator <=(const modint &x){return val<=x.val;}
	bool operator >=(const modint &x){return val>=x.val;}
	bool operator >(const modint &x){return val>x.val;}
	bool operator <(const modint &x){return val<x.val;}
	modint& operator *=(const modint &x){return val=static_cast<int>(1ll*val*x.val%mod),*this;}
	modint& operator <<=(const modint &x){return val=(1ll*val<<x.val)%mod,*this;}
	modint& operator +=(const modint &x){return val=(1ll*val+x.val)%mod,*this;}
	modint& operator -=(const modint &x){return val=norm(1ll*val-x.val),*this;}
	modint& operator >>=(const modint &x){return val>>=x.val,*this;}
	modint& operator ^=(const modint &x){return val^=x.val,*this;}
	modint operator >>(const modint &x){return modint(*this)>>=x;}
	modint operator <<(const modint &x){return modint(*this)<<=x;}
	modint& operator /=(const modint &x){return *this*=x.inv();}
	modint operator +(const modint &x){return modint(*this)+=x;}
	modint operator -(const modint &x){return modint(*this)-=x;}
	modint operator *(const modint &x){return modint(*this)*=x;}
	modint operator /(const modint &x){return modint(*this)/=x;}
	modint operator ^(const modint &x){return modint(*this)^=x;}
	friend std::ostream& operator<<(std::ostream& os,const modint &a){return os<<a.val;}
	friend std::istream& operator>>(std::istream& is,modint &a){return is>>a.val;}
};

typedef modint<167772161>m16;

inline int read(){
	int res=0,v=1;
	char c=getchar();
	while(c<'0'||c>'9'){v=(c=='-'?-1:1);c=getchar();}
	while(c>='0'&&c<='9'){res=(res<<3)+(res<<1)+(c^48);c=getchar();}
	return res*v;
}

m16 ksm(m16 a,int b){
	m16 ans=1;
	for(;b;b>>=1){
		if(b&1)ans=ans*a;
		a=a*a;
	}
	return ans;
}

m16 inv[Max],frac[Max],Inv[Max];

void init(int len) {
	frac[0]=inv[0]=1; 
	for(int i=1;i<=len;++i)frac[i]=frac[i-1]*i;
	Inv[len]=frac[len].inv();for(int i=len-1;i;--i)Inv[i]=Inv[i+1]*(i+1);
	for(int i=1;i<=len;++i)inv[i]=Inv[i]*frac[i-1];
}

int rev[Max];

void init(int lim,int len){for(int i=0;i<lim;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));}

struct Poly{  //导和积都到n，否则到n-1; 
	const m16 g=3;
	const m16 invg=(mod+1)/3;
	
	m16 a[Max];
	void out(int n,int x=0){for(int i=0;i<n;++i)cout << a[i].val << " ";if(!x)cout << "\n";}
	void inte(int n){for(int i=n;i>=0;--i)a[i]=a[i-1]*inv[i];a[0]=0;}
	void ReadIn(int n){for(int i=0;i<n;++i)a[i].val=read();}
	void der(int n){for(int i=0;i<=n;++i)a[i]=a[i+1]*(i+1);}
	m16& operator [](const int x){return a[x];}
	
	void NTT(int tag,int lim){
		for(int i=0;i<lim;++i)if(i<rev[i])swap(a[i],a[rev[i]]);
		for(int mid=1;mid<lim;mid<<=1){
			m16 W=ksm(tag==1?g:invg,(mod-1)/(mid<<1));
			for(int R=mid<<1,j=0;j<lim;j+=R){
				m16 w=1;
				for(int p=0;p<mid;++p,w=w*W){
					m16 x=a[j+p],y=a[j+p+mid]*w;
					a[j+p]=x+y;
					a[j+p+mid]=x-y;
				}
			}
		}
		if(tag==-1){
			m16 inv=m16(lim).inv();
			for(int i=0;i<lim;++i)a[i]*=inv;
		}
	}
	
}a;

void GetInv(Poly &y,Poly &x,int n){    //第二个为参数，第一个为返回值
	static Poly b[2],c;
	int lim=2,len=1,bas=1;
	int pos=0;
	b[pos][0]=x[0].inv();
	while(bas<=(n<<1)){
		init(lim,len);
		pos^=1;
		for(int i=0;i<bas;++i)b[pos][i]=b[pos^1][i]<<1;
		for(int i=0;i<=lim;++i)c[i]=x[i];
		for(int i=lim>>1;i<=lim;++i)c[i]=0;
		b[pos^1].NTT(1,lim);c.NTT(1,lim);
		for(int i=0;i<=lim;++i)b[pos^1][i]=b[pos^1][i]*b[pos^1][i]*c[i];
		b[pos^1].NTT(-1,lim);
		for(int i=0;i<bas;++i)b[pos][i]-=b[pos^1][i],b[pos^1][i]=0;
		bas<<=1,lim<<=1,++len;
	}
	for(int i=0;i<n;++i)y[i]=b[pos][i];
	for(int i=0;i<bas;++i)b[0][i]=b[1][i]=c[i]=0;
}

void GetLn(Poly &y,Poly&x,int n){
	static Poly t1,t2;
	for(int i=0;i<n;++i)t1[i]=t2[i]=x[i];
	t1.der(n);GetInv(t2,t2,n);
	int lim=1,len=0;
	while(lim<=n<<1)lim<<=1,++len;
	init(lim,len);
	t1.NTT(1,lim),t2.NTT(1,lim) ;
	for(int i=0;i<lim;++i)t1[i]*=t2[i];
	t1.NTT(-1,lim);t1.inte(n) ;
	for(int i=0;i<n;++i)y[i]=t1[i];
	for(int i=0;i<lim;++i)t1[i]=t2[i]=0;
}


void GetExp(Poly &y,Poly &x,int n){
	static Poly s,tmp;s[0]=1;
	int bas=1,lim=2,len=1;
	while(bas<=n<<1){
		for(int i=0;i<bas;++i)tmp[i]=s[i];
		GetLn(tmp,tmp,bas);
		for(int i=0;i<bas;++i)tmp[i]=x[i]-tmp[i];
		tmp[0]+=1;init(lim,len);tmp.NTT(1,lim),s.NTT(1,lim) ;
		for(int i=0;i<lim;++i)s[i]=s[i]*tmp[i];
		s.NTT(-1,lim) ;
		for(int i=bas;i<lim;++i)s[i]=0;
		lim<<=1,bas<<=1,++len;
	}
	for(int i=0;i<n;++i)y[i]=s[i];
	for(itn i=0;i<bas;++i)s[i]=tmp[i]=0;
}

void expow(Poly &y,Poly &x,int n,int k){
	static Poly s;
	for(itn i=0;i<=n;++i)s[i]=x[i];GetLn(s,s,n);
	for(int i=0;i<=n;++i)s[i]*=k;GetExp(s,s,n);
	for(int i=0;i<=n;++i)y[i]=s[i];
	for(int i=0;i<=n;++i)s[i]=0;
}


Poly b;

bool Med;
signed main(){
	int n=read()+1;
	int k=read();
	init(n<<3);
	for(int i=1;i<n;++i)b[i-1]=Inv[i];
	expow(b,b,n,k);
	for(int i=0;i<n;++i){
		if(i<k)cout << 0 << " ";
		else cout << (b[i-k]*Inv[k]*frac[i])<< ' ';
	}
	

	cerr<< "Time: "<<clock()/1000.0 << "s\n";
	cerr<< "Memory: " << (&Mst-&Med)/1000000.0 << "MB\n";
	return 0;
}
/*

第一类斯特林数

定义

第一类斯特林数 \(\begin{bmatrix}n\\k\end{bmatrix}\) 表示将 \(n\) 个区分的元素划分为 \(k\) 个互不区分的非空轮换中的方案数。

递推公式

\[\begin{bmatrix}n\\m\end{bmatrix}=\begin{bmatrix}n-1\\m-1\end{bmatrix}+(n-1)\begin{bmatrix}n-1\\m\end{bmatrix} \]

因为加入一个元素时有以下两种情况：

1. 加在之前一个轮换中，此时它可以加在 \(n-1\) 个元素的右边，故有 \((n-1)\begin{bmatrix}n-1\\m\end{bmatrix}\) 种方案
2. 新开一个轮换，故有 \(\begin{bmatrix}n-1\\m-1\end{bmatrix}\) 种方案

第一类斯特林数·行

考虑像第二类斯特林数

第一类斯特林数·列

可考虑对于每个轮换生成函数为 \(G(x)=\sum\limits_{i\geq 0}\frac{x^i}i\) 所有可以直接多项式快速幂求 \(F(x)=G(x)^m\)。

#include <bits/stdc++.h>
#define pii pair<int,int>
#define ll long long
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
bool Mst;
const int Max=2e6+10;
const int mod=167772161;
const int inf=1e9+10;

template <int mod>
struct modint{

	int val;

	static int norm(const int &x){return x<0?x+mod:x;}

	modint inv()const{
		int a=val,b=mod,u=1,v=0,t;
		while(b>0)t=a/b,swap(a-=t*b,b),swap(u-=t*v,v);
		return modint(u);
	}
	
	
	modint():val(0){}
	modint(const int &m):val(norm(m)){}
	modint(const ll &m):val(norm(m%mod)){}
	modint operator -()const{return modint(norm(-val));}
	bool operator ==(const modint &x){return val==x.val;}
	bool operator !=(const modint &x){return val!=x.val;}
	bool operator <=(const modint &x){return val<=x.val;}
	bool operator >=(const modint &x){return val>=x.val;}
	bool operator >(const modint &x){return val>x.val;}
	bool operator <(const modint &x){return val<x.val;}
	modint& operator *=(const modint &x){return val=static_cast<int>(1ll*val*x.val%mod),*this;}
	modint& operator <<=(const modint &x){return val=(1ll*val<<x.val)%mod,*this;}
	modint& operator +=(const modint &x){return val=(1ll*val+x.val)%mod,*this;}
	modint& operator -=(const modint &x){return val=norm(1ll*val-x.val),*this;}
	modint& operator >>=(const modint &x){return val>>=x.val,*this;}
	modint& operator ^=(const modint &x){return val^=x.val,*this;}
	modint operator >>(const modint &x){return modint(*this)>>=x;}
	modint operator <<(const modint &x){return modint(*this)<<=x;}
	modint& operator /=(const modint &x){return *this*=x.inv();}
	modint operator +(const modint &x){return modint(*this)+=x;}
	modint operator -(const modint &x){return modint(*this)-=x;}
	modint operator *(const modint &x){return modint(*this)*=x;}
	modint operator /(const modint &x){return modint(*this)/=x;}
	modint operator ^(const modint &x){return modint(*this)^=x;}
	friend std::ostream& operator<<(std::ostream& os,const modint &a){return os<<a.val;}
	friend std::istream& operator>>(std::istream& is,modint &a){return is>>a.val;}
};

typedef modint<167772161>m16;

inline int read(){
	int res=0,v=1;
	char c=getchar();
	while(c<'0'||c>'9'){v=(c=='-'?-1:1);c=getchar();}
	while(c>='0'&&c<='9'){res=(res<<3)+(res<<1)+(c^48);c=getchar();}
	return res*v;
}

m16 ksm(m16 a,int b){
	m16 ans=1;
	for(;b;b>>=1){
		if(b&1)ans=ans*a;
		a=a*a;
	}
	return ans;
}

m16 inv[Max],frac[Max],Inv[Max];

void init(int len) {
	frac[0]=inv[0]=1; 
	for(int i=1;i<=len;++i)frac[i]=frac[i-1]*i;
	Inv[len]=frac[len].inv();for(int i=len-1;i;--i)Inv[i]=Inv[i+1]*(i+1);
	for(int i=1;i<=len;++i)inv[i]=Inv[i]*frac[i-1];
}

int rev[Max];

void init(int lim,int len){for(int i=0;i<lim;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));}

struct Poly{  //导和积都到n，否则到n-1; 
	const m16 g=3;
	const m16 invg=(mod+1)/3;
	
	m16 a[Max];
	void out(int n,int x=0){for(int i=0;i<n;++i)cout << a[i].val << " ";if(!x)cout << "\n";}
	void inte(int n){for(int i=n;i>=0;--i)a[i]=a[i-1]*inv[i];a[0]=0;}
	void ReadIn(int n){for(int i=0;i<n;++i)a[i].val=read();}
	void der(int n){for(int i=0;i<=n;++i)a[i]=a[i+1]*(i+1);a[n]=0;}
	m16& operator [](const int x){return a[x];}
	
	void NTT(int tag,int lim){
		for(int i=0;i<lim;++i)if(i<rev[i])swap(a[i],a[rev[i]]);
		for(int mid=1;mid<lim;mid<<=1){
			m16 W=ksm(tag==1?g:invg,(mod-1)/(mid<<1));
			for(int R=mid<<1,j=0;j<lim;j+=R){
				m16 w=1;
				for(int p=0;p<mid;++p,w=w*W){
					m16 x=a[j+p],y=a[j+p+mid]*w;
					a[j+p]=x+y;
					a[j+p+mid]=x-y;
				}
			}
		}
		if(tag==-1){
			m16 inv=m16(lim).inv();
			for(int i=0;i<lim;++i)a[i]*=inv;
		}
	}
	
}a;

void GetInv(Poly &y,Poly &x,int n){    //第二个为参数，第一个为返回值
	static Poly b[2],c;
	int lim=2,len=1,bas=1;
	int pos=0;
	b[pos][0]=x[0].inv();
	while(bas<(n<<1)){
		init(lim,len);
		pos^=1;
		for(int i=0;i<bas;++i)b[pos][i]=b[pos^1][i]<<1;
		for(int i=0;i<=lim;++i)c[i]=x[i];
		for(int i=lim>>1;i<=lim;++i)c[i]=0;
		b[pos^1].NTT(1,lim);c.NTT(1,lim);
		for(int i=0;i<=lim;++i)b[pos^1][i]=b[pos^1][i]*b[pos^1][i]*c[i];
		b[pos^1].NTT(-1,lim);
		for(int i=0;i<bas;++i)b[pos][i]-=b[pos^1][i],b[pos^1][i]=0;
		bas<<=1,lim<<=1,++len;
	}
	for(int i=0;i<n;++i)y[i]=b[pos][i];
	for(int i=0;i<bas;++i)b[0][i]=b[1][i]=c[i]=0;
}

void GetLn(Poly &y,Poly&x,int n){
	static Poly t1,t2;
	for(int i=0;i<n;++i)t1[i]=t2[i]=x[i];
	t1.der(n);GetInv(t2,t2,n);
	int lim=1,len=0;
	while(lim<n<<1)lim<<=1,++len;
	init(lim,len);
	t1.NTT(1,lim),t2.NTT(1,lim) ;
	for(int i=0;i<lim;++i)t1[i]*=t2[i];
	t1.NTT(-1,lim);t1.inte(n) ;
	for(int i=0;i<n;++i)y[i]=t1[i];
	for(int i=0;i<lim;++i)t1[i]=t2[i]=0;
}


void GetExp(Poly &y,Poly &x,int n){
	static Poly s,tmp;s[0]=1;
	int bas=1,lim=2,len=1;
	while(bas<n<<1){
		for(int i=0;i<bas;++i)tmp[i]=s[i];
		GetLn(tmp,tmp,bas);
		for(int i=0;i<bas;++i)tmp[i]=x[i]-tmp[i];
		tmp[0]+=1;init(lim,len);tmp.NTT(1,lim),s.NTT(1,lim) ;
		for(int i=0;i<lim;++i)s[i]=s[i]*tmp[i];
		s.NTT(-1,lim) ;
		for(int i=bas;i<lim;++i)s[i]=0;
		lim<<=1,bas<<=1,++len;
	}
	for(int i=0;i<n;++i)y[i]=s[i];
	for(itn i=0;i<bas;++i)s[i]=tmp[i]=0;
}

void expow(Poly &y,Poly &x,int n,int k){
	static Poly s;
	for(itn i=0;i<n;++i)s[i]=x[i];GetLn(s,s,n);
	for(int i=0;i<n;++i)s[i]*=k;GetExp(s,s,n);
	for(int i=0;i<n;++i)y[i]=s[i];
	for(int i=0;i<n;++i)s[i]=0;
}


Poly b;

bool Med;
signed main(){
	int n=read()+1;
	int k=read();
	init(n);
	for(int i=0;i<n-1;++i)b[i]=inv[i+1];
	expow(b,b,n,k);
	m16 z=frac[k].inv();
	for(int i=0;i<n;++i){
		if(i<k)cout << 0 << " ";
		else cout << (b[i-k]*z*frac[i]).val<< ' ';
	}
	

	cerr<< "Time: "<<clock()/1000.0 << "s\n";
	cerr<< "Memory: " << (&Mst-&Med)/1000000.0 << "MB\n";
	return 0;
}
/*

*/

P5401 [CTS2019] 珍珠

设 \(cnt_i\) 表示颜色为 \(i\) 的珍珠个数。那么对于一个合法序列满足：

\[\begin{aligned} \sum_{i=1}^{D} \lfloor\frac{cnt_i}{2}\rfloor&\geq m\\ \sum_{i=1}^D\frac{cnt_i-(cnt_i\mod 2)}{2}&\geq m\\ \sum_{i=1}^Dcnt_i-(cnt_i\mod2)&\geq2m\\ n-\sum_{i=1}^Dcnt_i\mod2&\geq 2m\\ \sum_{i=1}^D cnt_i \mod2&\leq n-2m \end{aligned} \]

然后可以特判两种情况，分别是 \(n<2m\) 和 \(2m\leq n-D\) ，答案分别为 0 和 \(D^n\) 。

对于其他情况，我们发现答案只与每种珍珠出现的奇偶性有关，那么我们设 \(odd_i\) 表示恰好有 \(i\) 个奇数的方案数，那么 \(ans=\sum\limits_{i=0}^{n-2m} odd_i\) 。

看到恰好可以考虑容斥，定义 \(f_i\) 表示钦定 \(i\) 个是奇数，其余的不进行限制的方案数，那么 \(f_i=\sum\limits_{j=i}^{D} \binom{j}{i} odd_j\)

由二项式反演可以得到：

\[\begin{aligned} odd_i&=\sum_{j=i}^D (-1)^{j-i}\binom{j}{i}f(j)\\ &=\sum_{j=i}^D(-1)^{j-i}\frac{j!}{i!(j-i)!} f(j)\\ &=\frac{1}{i!} \sum_{j=i}^D \frac{(-1)^{j-i}}{(j-i)!}j!f(j)\\ \end{aligned} \]

通过EGF的卷积定义：

\[\begin{aligned} (F*G)[k]&=\sum_{i+j=k} \binom{k}{i}F[i]G[j]\\ \frac{(F*G)[k]}{k!}&=\sum_{i+j=k}\frac{F[i]}{i!}\frac{G[j]}{j!} \end{aligned} \]

我们可以发现上面式子中：

\[\sum_{j=i}^D \frac{(-1)^{j-i}}{(j-i)!}j!f(j) \]

是个卷积形式。

考虑先求 \(f(i)\) ，因为奇数位是1的EGF是：

\[\begin{aligned} \sum\limits_{i=0}\frac{[i\%2==1]}{i!}x^i&=\frac{e^x-e^{-x}}{2}\\ \therefore f_k&=\binom{D}{k}n^k(e^x)^{D-k}\\ &=\binom{D}{k}n![x^n]\frac{1}{2^k}(e^x-e^{-x})^k(e^x)^{D-k}\\ &=\binom{D}{k}\frac{n!}{2^k}[x^n](e^x-e^{-x})^k(e^x)^{D-k}\\ \end{aligned} \]

由二项式定理 \((x+y)^k=\sum\limits_{i=0}^k\binom{k}{i}x^iy^{k-i}\) 可知:

\[\begin{aligned} f_k&=\binom Dk \frac{n!}{2^k}[x^n]\sum_{i=0}^k\binom kie^{ix}(-e^x)^{k-i}(e^x)^{D-k}\\ &=\binom Dk \frac{n!}{2^k}\sum_{i=0}^k\binom ki(-1)^{k-i}e^{ix-(k-i)x+(D-k)x}[x^n]\\ &=\binom Dk \frac{n!}{2^k}\sum_{i=0}^k\binom ki(-1)^{k-i}[x^n]e^{(D-2(k-i))x}\\ \because [x^n]e^{ax}&=\frac{a^n}{n!}\\ \therefore f_k&=\binom Dk\frac{n!}{2^k}\sum_{i=0}^k\binom ki(-1)^{k-i}\frac{(D-2(k-j))^n}{n!}\\ &=\frac{D!}{k!(D-k)!}\frac{n!}{2^k}\sum_{i=0}^k\frac{k!}{i!(i-k)!}(-1)^{k-i}\frac{(D-2(k-j))^n}{n!}\\ 令 j&=k-i\\ \therefore f_k&=\frac{D!}{k!(D-k)!}\frac{n!}{2^k}\sum_{j=0}^k\frac{k!}{j!(k-j)!}(-1)^j\frac{(D-2j)^n}{n!}\\ &=\frac{D!}{(D-k)!2^k}\sum_{i=0}^k\frac{(-1)^i(D-2i)^n}{i!}\frac{1}{(k-i)!} \end{aligned} \]

显然后面求和号中是一个EGF的卷积形式，于是我们可以求出 \(f_k\) 。

#include <bits/stdc++.h>
#define int long long
#define reaD read
#define raed read
#define haed head
#define cotu cout
using namespace std;
const int Max=3e6+10;
const int mod=998244353;
const int g=3;
const int gi=332748118;



inline int read(){
	int res=0,v=1;
	char c=getchar();
	while(c<'0'||c>'9'){v=(c=='-'?-1:1);c=getchar();}
	while(c>='0'&&c<='9'){res=(res<<3)+(res<<1)+(c^48);c=getchar();}
	return res*v;
}

int a[Max],b[Max];


namespace Ntt{
	
	int ksm(int a,int b){
		int ans=1;
		for(;b;b>>=1){
			if(b&1)ans=ans*a%mod;
			a=a*a%mod;
		}
		return ans;
	}

	int p[Max],lim,pos;
	
	void rever(int lim,int pos){
		for(int i=1;i<lim;++i){
			p[i]=(p[i>>1]>>1)|((i&1)<<(pos-1));
		}
	}

	void NTT(int *a,int tag){
		rever(lim,pos);
		for(int i=0;i<lim;++i){
			if(i<p[i])swap(a[i],a[p[i]]);
		}
		for(int mid=1;mid<lim;mid<<=1){
			int W=ksm(tag==1?g:gi,(mod-1)/(mid<<1))%mod;
			for(int r=mid<<1,j=0;j<lim;j+=r){
				int w=1;
				for(int k=0;k<mid;++k,w=w*W%mod){
					int q=w*a[j+k+mid]%mod;
					int l=a[j+k];
					a[j+k]=(l+q)%mod;
					a[j+k+mid]=(l-q+mod)%mod;
				}
			}
		}
		if(tag==-1){
			int inv=ksm(lim,mod-2) ;
			for(int i=0;i<=lim;++i)a[i]=a[i]*inv%mod;
		}
	}
	
	void Solve(int *a,int *b){
		NTT(a,1);
		NTT(b,1);
		for(int i=0;i<=lim;++i)a[i]=a[i]*b[i]%mod;
		NTT(a,-1);
	}
}

using namespace Ntt;

int inv[Max],frac[Max];

void init(int n){
	frac[0]=1;
	for(int i=1;i<=n;++i){
		frac[i]=frac[i-1]*i%mod;
		inv[i]=ksm(frac[i],mod-2);
	}
	inv[0]=1;
}

signed main(){
	int d,n,m;
	d=read(),n=read(),m=read();
	init(200000);
	if(n<2*m){
		cout << 0 << "\n";
		return 0;
	}
	if(2*m<=n-d){
		cout << ksm(d,n) %mod << "\n";
		return 0;
	}
	
	//处理f 
	for(int i=0;i<=d;++i){
		a[i]=((i%2==0?1:-1)*ksm(d-2*i,n)%mod*inv[i]%mod+mod)%mod;;
	}
	for(int i=0;i<=d;++i){
		b[i]=inv[i];
	}
	lim=1;
	pos=0;
	while(lim<=2*d)lim<<=1,++pos;
	rever(lim,pos);
	Solve(a,b);
	for(int i=0;i<=d;++i)a[i]=a[i]*frac[d]%mod*inv[d-i]%mod*ksm(ksm(2,i)%mod,mod-2)%mod*frac[i]%mod;
	for(int i=d+1;i<=lim;++i)a[i]=b[i]=0;
	
	//处理odd 
	for(int i=0;i<=d;++i){
		b[d-i]=(((i%2==0)?inv[i]:-inv[i])%mod+mod)%mod;
	}
	Solve(a,b);
	int ans=0;
	for(int i=0;i<=n-2*m;++i){
		ans=(ans+a[i+d]*inv[i]%mod)%mod;
	}
	printf("%lld\n",ans);
	return 0;
}

P4980 【模板】Polya 定理

\[\begin{aligned} ans&=\frac1n\sum_{i=1}^nn^{\gcd(n,i)}\\ &=\frac1n\sum_{d|n}n^d\sum_{i=1}^{n}[\gcd(n,i)==d]\\ &=\frac1n\sum_{d|n}n^d\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}[\gcd(\frac{n}{d},i)==1]\\ &=\frac1n\sum_{d|n}n^d\varphi(\frac{n}{d}) \end{aligned} \]

#include <bits/stdc++.h>
#define int long long
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
const int Max=1e5+10;
const int mod=1e9+7;
const int inf=1e9+10;



inline int read(){
	int res=0,v=1;
	char c=getchar();
	while(c<'0'||c>'9'){v=(c=='-'?-1:1);c=getchar();}
	while(c>='0'&&c<='9'){res=(res<<3)+(res<<1)+(c^48);c=getchar();}
	return res*v;
}

int ksm(int a,int b){
	int ans=1;
	for(;b;b>>=1){
		if(b&1)ans=ans*a%mod;
		a=a*a%mod;
	}
	return ans;
}

int phi(int x){
	int ans=x;
	for(int i=2;i*i<=x;++i){
		if(x%i==0){
			ans-=ans/i;
			while(x%i==0)x/=i;
		}
	}
	if(x>1)ans-=ans/x;
	return ans;
}

int GetAns(int n,int d){
	return ksm(n,d)*phi(n/d)%mod;
}

int Get(){
	int n=read();
	int ans=0;
	for(int i=1;i*i<=n;++i){
		if(n%i==0){
			ans=(ans+GetAns(n,i))%mod;
			if(i*i!=n)ans=(ans+GetAns(n,n/i))%mod;
		}
	}
	return ans*ksm(n,mod-2)%mod;
}

signed main(){
	int T=read();
	while(T--)cout << Get() << "\n";

	return 0;
}

P10322 高洁（Purity）

\[\sum_{i=1}^ni^{v(i)+1}\\ d=\prod_{i=1}p_i^{t_i}\\ 设t_{max}为\max(t_1,t_2\dots) \\ 令：f(x)=\sum_{i=1}^ni^{x+1}[v(i)==x]\\ \therefore ans=\sum_{i=1}^{t_{max}}f(i)\\ 设d_t=\prod_{i=1}p_{i}^{\lceil\frac{t_i}{k}\rceil}\\ \therefore d|i^k等价于d_k|i\\ 并且[v(i)==k] 等价于[d_{k-1}\nmid i]\&\&[d_k|i]\\ \therefore f(k)= \begin{cases} \sum\limits_{i=1}^ni[d_{t_{max}}\nmid i]&&&k=0\\ \sum\limits_{i=1}^ni^2[d_1|i]&&&k=1\\ \sum\limits_{i=1}^ni^{k+1}[d_{k-1}\nmid i][d_k|i]&&&k>1 \end{cases}\\ \]

当 \(k=0\) 时：

\[\begin{aligned} f(k)&=\sum_{i=1}^ni[1-d_{t_{max}}|i]\\ &=\sum_{i=1}^ni-\sum_{i=1}^ni[d_{t_{max}}|i]\\ &=\sum_{i=1}^ni-d_{t_{max}}\sum_{i=1}^{\lfloor\frac{n}{d_{t_{max}}}\rfloor}i \end{aligned} \]

当 \(k=1\) 时：
转枚举 \(i\) 变为枚举 \(i\) 的因子

\[\begin{aligned} f(k)&=\sum_{i=1}^{\lfloor\frac{n}{d_1}\rfloor}(id_1)^2\\ &=d_1^2\sum_{i=1}^{\lfloor\frac{n}{d_1}\rfloor}(i)^2\\ \end{aligned} \]

当 \(k>1\) 时：

\[\begin{aligned} f(k)&=\sum_{i=1}^{\lfloor\frac{n}{d_k}\rfloor}(id_k)^{k+1}[\frac{d_{k-1}}{d_k}\nmid i]\\ &=d_k^{k+1}\sum_{i=1}^{\lfloor\frac{n}{d_k}\rfloor}i^{k+1}[1-d_{k-1}|i]\\ &=d_{k}^{k+1}(\sum_{i=1}^{\lfloor\frac{n}{d_k}\rfloor}i^{k+1}-\sum_{i=1}^{\lfloor\frac{n}{d_k}\rfloor}i^{k+1}[\frac{d_{k-1}}{d_k}|i])\\ &=d_{k}^{k+1}(\sum_{i=1}^{\lfloor\frac{n}{d_k}\rfloor}i^{k+1}-\sum_{i=1}^{\lfloor\frac{n}{d_{k-1}}\rfloor}(i\frac{d_{k-1}}{d_k})^{k+1})\\ &=d_k^{k+1}\sum_{i=1}^{\lfloor\frac{n}{d_k}\rfloor}i^{k+1}-d_{k-1}^{k+1}\sum_{i=1}^{\lfloor\frac{n}{d_{k-1}}\rfloor}i^{k+1} \end{aligned} \]

设 \(sum(n,k)=\sum\limits_{i=1}^ni^k\)

\[f(k)= \begin{cases} sum(n,1)-d_{t_{max}}sum(\lfloor\frac{n}{d_{t_{max}}}\rfloor,1)&&k=0\\ d_1^2sum(\lfloor\frac{n}{d_1}\rfloor,2)&&k=1\\ d_k^{k+1}sum(\lfloor\frac{n}{d_k}\rfloor,k+1)-d_{k-1}^{k+1}sum(\lfloor\frac{n}{d_{k-1}}\rfloor,k+1)&&k>1 \end{cases} \]

#include <bits/stdc++.h>
#define int __int128
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
const int Max=1e5+10;
const int mod=998244353;
const int inf=1e9+10;
namespace Fread {
	const int SIZE=1<<21;
	char buf[SIZE],*S,*T;
	inline char getchar() {
	    if(S==T){
	        T=(S=buf)+fread(buf,1,SIZE,stdin);
	        if(S==T)return '\n';
	    }
	    return *S++;
	}
}
namespace Fwrite {
	const int SIZE=1<<21;
	char buf[SIZE],*S=buf,*T=buf+SIZE;
	inline void flush(){
	    fwrite(buf,1,S-buf,stdout);
	    S=buf;
	}
	inline void putchar(char c){
	    *S++=c;
	    if(S==T)flush();
	}
	struct POPOSSIBLE{
	    ~POPOSSIBLE(){flush();}
	}ztr;
}
#define getchar Fread :: getchar
#define putchar Fwrite :: putchar
namespace Fastio{
	struct Reader{
	    template<typename T>
    	Reader& operator >> (T& x) {
        	char c=getchar();
        	T f=1;
        	while(c<'0'||c>'9'){
            	if(c=='-')f=-1;
            	c=getchar();
        	}
        	x=0;
        	while(c>='0'&&c<='9'){
	            x=x*10+(c-'0');
	            c=getchar();
	        }
	        x*=f;
	        return *this;
    	}
	    Reader& operator >> (char& c) {
	        c=getchar();
	        while(c==' '||c=='\n')c=getchar();
	        return *this;
	    }
	    Reader& operator >> (char* str) {
	        int len=0;
	        char c=getchar();
	        while(c==' '||c=='\n')c=getchar();
	        while(c!=' '&&c!='\n'&&c!='\r'){
	            str[len++]=c;
            	c=getchar();
        	}
			str[len]='\0';
			return *this;
	    }
	    Reader(){}
	}cin;
	struct Writer{
	    template<typename T>
	    Writer& operator << (T x) {
	        if(x==0){putchar('0');return *this;}
	        if(x<0){putchar('-');x=-x;}
	        static int sta[45];
	        int top=0;
	        while(x){sta[++top]=x%10;x/=10;}
	        while(top){putchar(sta[top]+'0');--top;}
	        return *this;
    	}
    	Writer& operator << (char c) {
        	putchar(c);
        	return *this;
	    }
    	Writer& operator << (char* str) {
        	int cur=0;
        	while(str[cur])putchar(str[cur++]);
        	return *this;
    	}
    	Writer& operator << (const char* str) {
    	    int cur=0;
    	    while(str[cur])putchar(str[cur++]);
			return *this;
    	}
    	Writer(){}
	}cout;
}
#define cin Fastio :: cin
#define cout Fastio :: cout
#define endl Fastio :: endl


inline int read(){
	int x;
	cin>> x;
	return x;
}
int ksm(int a,int b){
	int ans=1;
	for(;b;b>>=1){
		if(b&1)ans=ans*a%mod;
		a=a*a%mod;
	}
	return ans%mod;
}

int P[Max],q[Max],s[Max];
int n,d,t,a[Max];



int sum(int n,int k){
	P[0]=1;
	q[0]=n%mod;
	for(int i=1;i<=k+1;++i){
		P[i]=P[i-1]*i%mod;
		q[i]=q[i-1]*(n-i)%mod;
	}
	s[k+2]=1;
	for(int i=k+1;i;--i){
		s[i]=s[i+1]*(n-i)%mod;
	}
	int ans=0;
	int res=0;
	for(int i=1;i<=k+1;++i){
		res=(res+ksm(i,k)%mod)%mod;
		ans=(ans+res*q[i-1]%mod*s[i+1]%mod*ksm(P[i]*((k-i+1)%2==0?1:-1)*P[k+1-i]%mod,mod-2)%mod)%mod;
	}
	return (ans%mod+mod)%mod;
}

int f(int k){
	if(k==0)return ((sum(n,1)-a[t]*sum(n/a[t],1)%mod)%mod+mod)%mod;
	if(k==1)return a[1]*a[1]%mod*sum(n/a[1],2)%mod;
	return ((ksm(a[k],k+1)*sum(n/a[k],k+1)%mod-ksm(a[k-1],k+1)*sum(n/a[k-1],k+1)%mod)%mod+mod)%mod;
}

vector<int>p;
vector<int>b;
void work(){
	n=read();
	d=read();
	if(d==1){
		cout << sum(n,2) << "\n";
		return;
	}
	t=0;
	p.clear();
	b.clear();
	for(int i=2;i*i<=d;++i){
		if(d%i==0){
			p.push_back(i);     //质因数分解 
			b.push_back(1);
			int z=b.size()-1;
			d/=i;
			while(d%i==0){
				b[z]++;
				d/=i;
			}
			t=max(t,b[z]);
		}
	}
	if(d!=1){
		p.push_back(d);
		b.push_back(1);
		int z=b.size()-1;
		t=max(t,b[z]);
	}
	int z=p.size();
	for(int i=1;i<=t;++i){
		a[i]=1;
		for(int j=0;j<z;++j){
			a[i]=a[i]*(ksm(p[j],ceil(1.0*b[j]/i)))%mod;  //计算d_k 
		}
	} 
	int ans=0;
	for(int i=0;i<=t;++i){
		ans=(ans+f(i)%mod)%mod;
	}
	cout << ans << "\n";
}

signed main(){
	int t=read();
	while(t--)work();
	return 0;
}

不知道为什么要开 __int128 才过，不知道哪里爆 long long 了。

P6271 [湖北省队互测2014] 一个人的数论

设 \(F(x)=\sum\limits_{i=1}^ni^k\)，\(G(x)=\sum\limits_{i=1}^ni^k[\gcd(i,n)==1]\)

\[\begin{aligned} \because F(x)&=\sum_{d|x}d^kG(\frac{x}{d})\\ \end{aligned} \]

由莫反可知：

\[\begin{aligned} G(x)&=\sum_{d|x}d^k\mu(d)F(\frac{x}{d})\\ \end{aligned} \]

因为F(x)是应该x+1次多项式，即:

\[\begin{aligned} \\F(x)&=\sum_{i=0}^xf_ix^i\\ \therefore G(x)&=\sum_{d|x}d^k\mu(d)\sum_{i=0}^{k+1}f_i(\frac{x}{d})^i\\ &=\sum_{i=0}^{k+1}f_ix^i\sum_{d|x}d^{k-i}\mu(d)\\ \end{aligned} \]

由于 \(\sum_{d|x}d^{k-i}\mu(d)\) 为积性函数

\[\begin{aligned} \because n&=\prod_{p_i|n}p_i^{a_i} \\ \therefore G(x)&=\sum_{i=0}^{k+1}f_ix^i\prod_{p_j|x}\sum_{z=0}^{a_i}\mu(p_j^z)p_j^{z(k-i)}\\ \therefore \mu(p^k)&= \begin{cases} 1&&&k=0\\ -1&&&k=1\\ 0&&&k>1 \end{cases}\\ \therefore G(x)&=\sum_{i=0}^{k+1}f_ix^i\prod_{p_j|x}(1-p_j^{k-i})\\ \therefore G(n)&=\sum_{i=0}^{k+1}f_in^i\prod_{p_j|n}(1-p_j^{k-i})\\ \end{aligned} \]

那么只需要插出 \(F(x)\) 的多项式就行。

#include <bits/stdc++.h>
#define int long long
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
const int Max=1e5+10;
const int mod=1e9+7;
const int inf=1e9+10;



inline int read(){
	int res=0,v=1;
	char c=getchar();
	while(c<'0'||c>'9'){v=(c=='-'?-1:1);c=getchar();}
	while(c>='0'&&c<='9'){res=(res<<3)+(res<<1)+(c^48);c=getchar();}
	return res*v;
}

int ksm(int a,int b){
	int ans=1;
	for(;b;b>>=1){
		if(b&1)ans=ans*a%mod;
		a=a*a%mod;
	}
	return ans;
}

int a[Max],f[Max];
int len;
void Mul(int c,int b){
	++len;
	for(int i=len;i>=1;--i){
		a[i]=(a[i]*b+a[i-1]*c)%mod;
	}
	a[0]=a[0]*b%mod;
}

void mul(int v){
	for(int i=0;i<=len;++i){
		a[i]=a[i]*v%mod;
	}
}

void add(){
	for(int i=0;i<=len;++i){
		f[i]=(f[i]+a[i])%mod;
		a[i]=0;
	}
	len=0;
	a[0]=1;
}

void init(int k){
	a[0]=1;
	len=0;
	int val=0;
	for(int i=1;i<=k+2;++i){
		val=(val+ksm(i,k))%mod;
		int m=ksm(val,mod-2);
		for(int j=1;j<=k+2;++j){
			if(i==j)continue;
			Mul(1,mod-j);
			m=m*(i+mod-j)%mod;
		}
		mul(ksm(m,mod-2));
		add();
	}
}

int p[Max];

signed main(){
	int k=read(),w=read();
	int n=1;
	for(int i=1;i<=w;++i){
		p[i]=read();
		n=n*ksm(p[i],read())%mod;
	}
	init(k);
//	for(int i=1;i<=k+2;++i)cout<< f[i] << ' ';
//	cout << "\n";
	int ans=0;
	int in=n;
	for(int i=1;i<=k+2;++i){
		int pos=in*f[i]%mod;
		in=in*n%mod;
		for(int j=1;j<=w;++j){
			pos=pos*(1+mod-ksm(p[j],k-i-1+mod)%mod)%mod;
		}
		ans=((ans+pos)%mod+mod)%mod;
	}
	cout << ans << "\n";
	return 0;
}

P10269 实力派

第一个问题：

\[\begin{aligned} Ans_1&=\sum_{T\subseteq S\&\&|T|=k}[\gcd_{x\in T}x=1]\\ &=\sum_{T\subseteq S\&\&|T|=k}\sum_{i|\gcd_{x\in T}x}\mu(i)\\ &=\sum_{i=1}^n\mu(i)C_{cnt_i}^k\\ \end{aligned} \]

\(cnt_i\) 表示 \(i\) 的倍数的出现次数，\(n\) 为值域。

第二个问题：

\[考虑容斥，如果一个数作为\gcd那么其他选出来的数必定是其倍数，\\并且不是其倍数的倍数\\ 例如：当i作为\gcd时，ans_i=i(C_{cnt_i}^k-C_{cnt_{2i}}^k-C_{cnt_{3i}}^k-C_{cnt_{5i}}^k+C_{cnt_{6i}}^k\dots)\\ 可以发现ans_i=i\sum_{j=1}\mu(j)C_{cnt_j}^k\\ \begin{aligned} \therefore Ans_2&=\sum_{i=1}^nans_i=\sum_{i=1}^ni\sum_{j=1}^{\lfloor\frac{n}{i}\rfloor}\mu(j)C_{cnt_{ij}}^k\\ &=\sum_{i=1}^nC_{cnt_i}^k\sum_{j|i}\mu(j)\frac{i}{j}\\ &=\sum_{i=1}^nC_{cnt_i}^k\sum_{j|i}\mu(j)id(\frac ij)\\ \because \mu\ * id &=\varphi\\ \therefore Ans_2&=\sum_{i=1}^nC_{cnt_i}^k\varphi(i)\\\ \end{aligned} \]

现在这样是 \(O(nm)\) 的暴力，肯定过不去，考虑优化。可以发现这两个式子只与 \(cnt_i\) 的值有关，那么我们可以将式子改写为：

\[\begin{aligned} Ans_1&=\sum_{i=1}^vC_{i}^k\sum_{j=1}^n\mu(j)[cnt_j==i]\\ Ans_2&=\sum_{i=1}^vC_{i}^k\sum_{j=1}^n\varphi(j)[cnt_j==i]\\ v&=\max_{i=1}^n\ cnt_i \end{aligned} \]

可以发现两个式子后面一个求和号中的值是不随 \(k\) 变化而变化的，那么我们可以花费 \(O(n)\) 的代价算出来其和为 \(S1_i\),\(S2_i\)，故：

\[Ans_1=\sum_{i=1}^vC_i^kS1_i\\ Ans_2=\sum_{i=1}^vC_i^kS2_i\\ \]

由于 \(cnt_i\) 的种类不会过多，大约是 \(O(\sqrt n)\) 种，可以暴力了。

#include <bits/stdc++.h>
#define ll long long
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
const int Max=1e6+10;
int mod=998244353;
const int inf=1e9+10;

namespace fast_IO {
#define IOSIZE 1000000
	char ibuf[IOSIZE], obuf[IOSIZE], *p1 = ibuf, *p2 = ibuf, *p3 = obuf;
#define getchar() ((p1==p2)and(p2=(p1=ibuf)+fread(ibuf,1,IOSIZE,stdin),p1==p2)?(EOF):(*p1++))
#define putchar(x) ((p3==obuf+IOSIZE)&&(fwrite(obuf,p3-obuf,1,stdout),p3=obuf),*p3++=x)
#define isdigit(ch) (ch>47&&ch<58)
#define isspace(ch) (ch<33)
	template<typename T> inline T read() { T s = 0; int w = 1; char ch; while (ch = getchar(), !isdigit(ch) and (ch != EOF)) if (ch == '-') w = -1; if (ch == EOF) return false; while (isdigit(ch)) s = s * 10 + ch - 48, ch = getchar(); return s * w; }
	template<typename T> inline bool read(T &s) { s = 0; int w = 1; char ch; while (ch = getchar(), !isdigit(ch) and (ch != EOF)) if (ch == '-') w = -1; if (ch == EOF) return false; while (isdigit(ch)) s = s * 10 + ch - 48, ch = getchar(); return s *= w, true; }
	template<typename T> inline void print(T x) { if (x < 0) putchar('-'), x = -x; if (x > 9) print(x / 10); putchar(x % 10 + 48); }
	inline bool read(char &s) { while (s = getchar(), isspace(s)); return true; }
	inline bool read(char *s) { char ch; while (ch = getchar(), isspace(ch)); if (ch == EOF) return false; while (!isspace(ch)) *s++ = ch, ch = getchar(); *s = '\000'; return true; }
	inline void print(char x) { putchar(x); }
	inline void print(char *x) { while (*x) putchar(*x++); }
	inline void print(const char *x) { for (int i = 0; x[i]; i++) putchar(x[i]); }
	inline bool read(std::string& s) { s = ""; char ch; while (ch = getchar(), isspace(ch)); if (ch == EOF) return false; while (!isspace(ch)) s += ch, ch = getchar(); return true; }
	inline void print(std::string x) { for (int i = 0, n = x.size(); i < n; i++) putchar(x[i]); }
	inline bool read(bool &b) { char ch; while(ch=getchar(), isspace(ch)); b=ch^48; return true; }
	inline void print(bool b) { putchar(b+48); }
	template<typename T, typename... T1> inline int read(T& a, T1&... other) { return read(a) + read(other...); }
	template<typename T, typename... T1> inline void print(T a, T1... other) { print(a), print(other...); }
	struct Fast_IO { ~Fast_IO() { fwrite(obuf, p3 - obuf, 1, stdout); } } io;
	template<typename T> Fast_IO& operator >> (Fast_IO &io, T &b) { return read(b), io; }
	template<typename T> Fast_IO& operator << (Fast_IO &io, T b) { return print(b), io; }
#define cout io
#define cin io
#define endl '\n'
} using namespace fast_IO;


inline int read(){
	int x;
	cin>> x;
	return x;
}
int n,m;

ll ksm(ll a,int b){
	ll ans=1;
	for(;b;b>>=1){
		if(b&1)ans=ans*a%mod;
		a=a*a%mod;
	}
	return ans;
}


int pri[Max],phi[Max],flag[Max],pos,mu[Max];

void pre(int n){
	mu[1]=phi[1]=1;
	for(int i=2;i<=n;++i){
		if(!flag[i]){
			pri[++pos]=i;
			phi[i]=i-1;
			mu[i]=-1;
		}
		for(int j=1;j<=pos&&pri[j]*i<=n;++j){
			flag[i*pri[j]]=1;
			if(i%pri[j]==0){
				phi[i*pri[j]]=phi[i]*pri[j];
				break;
			}
			mu[i*pri[j]]-=mu[i];
			phi[i*pri[j]]=phi[i]*(pri[j]-1);
		}
	}
}

int a[Max];

ll p[Max];
ll inv[Max];
ll C(int n,int m){
	return p[n]*inv[m]%mod*inv[n-m]%mod;//ksm(p[m],mod-2)%mod*ksm(p[n-m],mod-2)%mod;
}

ll val[Max];
ll ans1[Max],ans2[Max];

void add(ll &a,ll b){
	a=a+b;
	a+=(a<0?mod:0);
	a-=((a>=mod)?mod:0);
}

signed main(){
	n=read();
	m=read();
	mod=read();
	int v=0;
	p[0]=1;
	for(int i=1;i<=n;++i){
		a[i]=read();
		v=max(v,a[i]);
		p[i]=p[i-1]*i%mod;
		val[a[i]]++;
	}
	pre(max(n,v));
	for(int i=1;i<=v;++i){
		for(int j=2*i;j<=v;j+=i) {
			val[i]+=val[j];
		}
	}
	inv[n]=ksm(p[n],mod-2);
	for(int i=n-1;i>=0;--i){
		inv[i]=inv[i+1]*(i+1)%mod;
	}
	for(int i=1;i<=max(n,v);++i)mu[i]=mu[i]+(mu[i]<0?mod:0);
	for(int i=1;i<=v;++i){
		for(int k=1;k<=val[i];++k){
			add(ans2[k],C(val[i],k)*phi[i]%mod);
			add(ans1[k],C(val[i],k)*mu[i]%mod);
		}
	}
	for(int i=1;i<=m;++i){
		int x=read();
		cotu << ans1[x] << ' ' << ans2[x] << "\n";
	}

	return 0;
}

P7360 「JZOI-1」红包

\[\prod_{i_1=1}^n\prod_{i_2=1}^n\dots\prod_{i_k=i}^nlcm(i_1,i_2\dotsi_k)\mod 998244353\\ 考虑对每个质数计算其贡献：\\ a_{p,j}表示在j质因数分解以后的质因子p的个数\\（后文可能会省略p，若省略则默认为p_i） \\ ans=\prod_{p_i|n}p_i^{\sum\limits_{j_1=1}^n\sum\limits_{j_2=1}^n\dots\sum\limits_{j_k=1}^n\max(a_{p_i,j_1},a_{p_i,j_2},\dots a_{p_i,j_k})}\\ 尝试对p_i计算\sum\limits_{j_1=1}^n\sum\limits_{j_2=1}^n\dots\sum\limits_{j_k=1}^n\max(a_{p_i,j_1},a_{p_i,j_2},\dots a_{p_i,j_k})\\ 假设:t=\max(a_{j_1},a_{j_2},\dots a_{j_k})\\ 那么满足\max(a_{j_1},a_{j_2},\dots a_{j_k})<t时的j的可供取值的个数是： n-b_t\\等等，好像有锅\\ b_{p，j}表示在[1,n]中拥有只是j个质因子p的数的个数（省略同a\\ 那么存在k^{(n-b_t)}种取值使得\max(a_{j_1},a_{j_2}\dots a_{j_k})<t\\ 那么是不是可以对于每个p暴力枚举t，做差得到\max(a_{j_1},a_{j_2}\dots a_{j_k})=t的个数。\\ ans=\prod_{p_i|n}p_i^{\sum\limits_{t=1}^{p_i^t\leq n}t\times(k^{n-b_{t+1}}-k^{n-b_i})} \\ 问题关键在于求b_{p,j},这个我们可以花费 O(n\log n)的时间对每个数质因数分解，\\然后求前缀和就行\\ \]

现在问题在于暴力（或前缀和）的复杂度是多少（二者复杂度相同），可以列出计算式：

\[\sum_{p\in prime\&\&p<n}\log_pn \]

通过打表发现在 \(1e6\) 时这个值大约是 \(75000\)

但是好像复杂度有问题（我也不会算），主要是还要加快速幂，先搁一下

P4128 [SHOI2006] 有色图

因为有 \(Burnside\) 引理：

\[|X/G|=\frac{1}{|G|}\sum_{g\in G}|X^g|\\ \]

\(|X/G|\) 是群 \(G\) 作用在集合 \(X\) 的轨道。
\(X^g\) 是 \(g\) 作用在 \(X\) 上的不变的集合。

那么我们要求的就是 \(|X/G|\)。

由于对顶点重排是一个置换群，每一种置换是有若干轮换构成的，我们可以发现有若干边是等价的，即他们的颜色必须相同，而对于每个等价类都有 \(m\) 种选择，那么如果有 \(k\) 个等价类，那么就会存在 \(m^k\) 种染色方式，使其在经过若干次该操作后会恢复成原状。（注意：一个等价类并不等同于一个轮换，因为是经过若干操作后变化原状，并不是一次操作。）

现在问题的关键在于求解 \(k\)，如果我们将一个置换拆解为 \(len\) 个轮换，第 \(i\) 个轮换长度是 \(b_i\) 。

先考虑在同一个轮换中的贡献：我们可以发现，在同一个轮换的同一个等价类中的边长度是一样的，因为这样才能在若干次操作后恢复原因。由于在一个轮换中我们选择的是点，那么他们之间两两连边也应该取出，那么就构成了一个完全图。而在一个有 \(a\) 个点的完全图中长度不同的边有 \(\lfloor\frac{a}{2}\rfloor\) 中，于是这 \(len\) 个轮换的贡献是 \(\sum\limits_{i=1}^{len}\lfloor\frac{b_i}{2}\rfloor\) 。

这是在同一个轮换中，但是在一个轮换中我们只取出其内部连边，然而两两轮换中的连边还没有分类。那么如何统计两两轮换间的贡献呢？

我们可以发现在只两个轮换 \(i\) 和 \(j\) 时 ，只需要经过 \(\operatorname{lcm}(b_i,b_j)\) 次操作就会恢复原样，那么其等价类长度为 \(\operatorname{lcm}(b_i,b_j)\)，由于有 \(b_ib_j\) 条边，那么会产生 \(\frac{b_ib_j}{\operatorname{lcm}(b_i,b_j)}=\gcd(b_i,b_j)\) 个等价类，那么两两轮换中的贡献是：\(\sum\limits_{i=1}^{len}\sum\limits_{j=1}^{i-1}\gcd(b_i,b_j)\)。

那么可以得到 \(k=\sum\limits_{i=1}^{len}\lfloor\frac{b_i}{2}\rfloor+\sum\limits_{i=1}^{len}\sum\limits_{j=1}^{i-1}\gcd(b_i,b_j)\)，但是新的问题又出现了，就是该群的置换个数即 \(|G|\) 是 \(n!\) 的，我们没法去直接枚举所有置换。但是我们发现答案好像只与 \(b\) 数组有关系，那么我们是不是只用去找的 \(b\) 数组的所有情况然后去重就行了呢？

说做就做：将这 \(n\) 个数随便放置的方案是 \(n!\) 种，但是对于一个长度为 \(b_i\) 的轮换会重复计算 \(b_i!\) 种，那么总共有 \(\frac{n!}{\prod\limits_ib_i!}\) ，这是考虑不同轮换之间的方案，而对于一个轮换内部，其实是一个园排列，那么对于一个长度是 \(b_i\) 的轮换会存在 \((b_i-1)!\) 种不同方案，那么总方案就是 \(\prod\limits_i(b_i-1)!\) 种。有一开始的式子我们发现答案与轮换的顺序没有关系，那么我们还要去除因为顺序所带来的重复，假设 \(b_i\) 出现了 \(cnt_i\) 次，那么重复次数为 \(cnt_i!\)，所以全部重复次数为 \(\prod\limits_icnt_i!\)。

那么我们将最后式子写出来：

\[\begin{aligned} k&=\sum\limits_{i=1}^{len}\lfloor\frac{b_i}{2}\rfloor+\sum\limits_{i=1}^{len}\sum\limits_{j=1}^{i-1}\gcd(b_i,b_j)\\ |X/G|&=\frac{1}{|G|}\sum_{g\in G}\frac{m^{k}\times n!\times \prod(b_i-1)!}{\prod (b_i!cnt_i!)}\\ \end{aligned} \]

由于 \(|G|\) 是 \(n!\)，那么整理这个式子得到：

\[\begin{aligned} k&=\sum\limits_{i=1}^{len}\lfloor\frac{b_i}{2}\rfloor+\sum\limits_{i=1}^{len}\sum\limits_{j=1}^{i-1}\gcd(b_i,b_j)\\ |X/G|&=\sum_{g\in G}\frac{m^{k}}{\prod (b_icnt_i!)}\\ \end{aligned} \]

我们发现 \(b\) 是 \(n\) 的一种分拆形式，并且 \(n\) 本身不大，那么我们可以暴力枚举，计算贡献。

#include <bits/stdc++.h>
#define int long long
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
const int Max=1e2+10;
int mod=998244353;
const int inf=1e9+10;



inline int read(){
	int res=0,v=1;
	char c=getchar();
	while(c<'0'||c>'9'){v=(c=='-'?-1:1);c=getchar();}
	while(c>='0'&&c<='9'){res=(res<<3)+(res<<1)+(c^48);c=getchar();}
	return res*v;
}

int frac[Max],b[Max];
int n,ans=0,m;

void pre(){
	frac[0]=1;
	for(int i=1;i<=n;++i){
		frac[i]=frac[i-1]*i%mod;
	}
}

int ksm(int a,int b){
	int ans=1;
	for(;b;b>>=1){
		if(b&1)ans=ans*a%mod;
		a=a*a%mod;
	}
	return ans;
}


int work(int len){
	int pos=0;
	int u=1;
	for(int i=1;i<=len;++i){
		pos=pos+(b[i]/2);        //一个环 
		u=u*b[i]%mod;
		for(int j=1;j<i;++j){
			pos=pos+__gcd(b[i],b[j]);//环交贡献 
		}
	}
	for(int i=1;i<=len;){
		int j=i+1;
		while(b[j]==b[i]&&j<=len)++j;    //去重 
		u=u*frac[j-i]%mod;
		i=j;
	}
	return ksm(u%mod,mod-2)*ksm(m,pos)%mod;
}
//注意pos的值不能取模 


void dfs(int now,int last,int sum){  //分拆 
	if(!sum){
		ans=(ans+work(now-1))%mod;
		return;
	}
	for(int i=last;i<=sum;++i){
		b[now]=i;
		dfs(now+1,i,sum-i);
	}
}

signed main(){
	n=read();
	m=read();
	mod=read();
	pre();
	dfs(1,1,n);
	cout << ans << "\n";
	return 0;
}

然后对于 P4727 [HNOI2009] 图的同构计数 可以将问题转换为每条边染不染色，即 \(m=2\) 的情况求解。

P4619 [SDOI2018] 旧试题

令 \((i,j)\) 表示 \([\gcd(i,j)==1]\)

\[\sum_{i=1}^A\sum_{j=1}^B\sum_{k=1}^Cd(ijk)\\ \]

先考虑 \(C=1\) 时：

有结论 \(d(ij)=\sum\limits_{x|i}\sum\limits_{y|i}(x,y)\) 。因为我们可以对于每个质因子单独考虑，对于质因子 \(P\) 它在 \(i\) 中出现了 \(p_i\) 次，在 \(j\) 中出现 \(p_j\) 次，那么它会在 \(ij\) 中出现 \(p_i+p_j\) 次。我们知道 \(d(i)=\prod\limits_{P\in pri\&\&P|i} (p_i+1)\) 因为对于一个 \(P\) 你可以选择的次数是在 \([0,p_i]\) 中，根据乘法原理可知。那么质因子 \(P\) 在 \(d(i,j)\) 中的贡献也是 \(p_i+p_j+1\)。而当 \([\gcd(x,y)==1]\) 时，对于质因子 \(P\) 它要么不出现，要么只在 \(x\) 或 \(y\) 中出现（废话）。这样它的贡献是 \(p_i+p_j+1\) 现在考虑是不是乘法原理，因为对于每个 \(P\) 是相互独立的，假设在 \(x\) 中出现了 \(a_x\) 次，那么是不会影响另外质因子 \(P^{'}\) 。所以可以证明。

那么将这个结论推广到 \(C\neq 1\) 时，那么有 \(d(ijk)=\sum\limits_{x|i}\sum\limits_{y|i}\sum\limits_{z|k}(x,y)(y,z)(x,k)\) 。

那么：

\[\begin{aligned} ans&=\sum_{i=1}^A\sum_{j=1}^B\sum_{k=1}^C\sum_{x|i}\sum_{y|i}\sum_{z|k}(x,y)(y,z)(x,k)\\ &=\sum_{i=1}^A\sum_{j=1}^B\sum_{k=1}^C(x,y)(y,x)(x,k)\frac{A}{x} \end{aligned} \]

List Generation

考虑答案为从点 \((0,0)\) 到点 \((n,m)\) 的只向上或右走的路径长度之和，考虑枚举转折点，为了不重，将转折点分为两类：

1. 从左边来然后向上走
2. 从下边来然后向右走

于是考虑选 \(i\) 转折点 1 的方案有 \(\binom ni\binom mi\) ，然后路径上一共有 \(n+m-1\) 个点，由于已经钦定了 \(i\) 个点了，那么剩下 \(s=n+m-i-1\) 个点可选可不选，对于这 \(s\) 个点的贡献是：

\[\begin{aligned} \sum_{j=0}^s\binom sj(i+j+2)&=2^s(i+2)+\sum_{j=0}^s\binom sj j\\ &=2^s(i+2)+\frac{\sum_{j=0}^s\binom sj s}{2}\\ &=2^s(i+2)+s2^{s-1}\\ \therefore ans&=\sum_{i=0}^{\min(n,m)}\binom ni \binom mi(2^s(i+2)+s2^{s-1}) \end{aligned} \]

#include <bits/stdc++.h>
#define pii pair<int,int>
#define int long long
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
bool Mst;
const int Max=1e7+10;
const int mod=1e9+7;
const int inf=1e9+10;


inline int read(){
	int res=0,v=1;
	char c=getchar();
	while(c<'0'||c>'9'){v=(c=='-'?-1:1);c=getchar();}
	while(c>='0'&&c<='9'){res=(res<<3)+(res<<1)+(c^48);c=getchar();}
	return res*v;
}

int frc[Max],pr[Max],inv[Max];

int ksm(int a,int b){
	int ans=1;
	for(;b;b>>=1){
		if(b&1)ans=ans*a%mod;
		a=a*a%mod;
	}
	return ans;
}

void pre(){
	pr[0]=1;
	for(int i=1;i<=1e7;++i){
		pr[i]=pr[i-1]*2%mod;
	}
	inv[0]=frc[0]=1;
	for(int i=1;i<=5e6+2;++i){
		frc[i]=frc[i-1]*i%mod;
	}
	int z=5e6+2;
	inv[z]=ksm(frc[z],mod-2);
	for(int i=z-1;i>=1;--i){
		inv[i]=inv[i+1]*(i+1)%mod;
	}
}

int C(int n,int m){  //n选m 
	if(n<0||m<0||n<m)return 0;
	return frc[n]*inv[m]%mod*inv[n-m]%mod;
}

void work(){
	int n=read();
	int m=read();
	int ans=0;
	for(int i=0;i<=min(n,m);++i){
		int s=n+m-i-1;
		ans=(ans+(C(n,i)*C(m,i)%mod*(pr[s]*(i+2)%mod+(s?s*pr[s-1]%mod:0))%mod))%mod;
	}
	cout << ans << "\n";
}

bool Med;
signed main(){
	pre();
	int T=read();
	while(T--)work();


	cerr<< "Time: "<<clock()/1000.0 << "s\n";
	cerr<< "Memory: " << (&Mst-&Med)/1000000.0 << "MB\n";
	return 0;
}

P6031 CF1278F Cards 加强版

因为每次洗牌是独立的，所以一次抽中王牌的概率 \(p=\frac{1}{m}\)，所以有：

\[ans=\sum_{i=0}^n\binom nip^i(1-p)^{n-i}i^k \]

用第二类斯特林数拆幂次：

因为 \(m^n\) 等于将 \(n\) 个区分的小球放入 \(m\) 个区分的盒子中，允许有空的方案数，用第二类斯特林数表示就是：

\[m^n=\sum_{i=0}^{\min(n,m)} \binom{m}{i}i!\begin{Bmatrix}n\\i\end{Bmatrix} \]

代入原式可得：

\[\begin{aligned} ans&=\sum_{i=0}^n\binom ni p^i(1-p)^{n-i}\sum_{j=0}^{k} \binom ij j!\begin{Bmatrix}k\\j\end{Bmatrix}\\ &=\sum_{j=0}^{k}j!\begin{Bmatrix}k\\j\end{Bmatrix}\sum_{i=0}^n\binom ni \binom ij p^i(1-p)^{n-i}\\ \end{aligned} \]

考虑第二个求和号中的内容：

\[\begin{aligned} &=\sum_{i=0}^n\binom ni\binom ijp^i(1-p)^{n-i}\\ &=\sum_{i=0}^n \frac{n!i!}{i!(n-i)!j!(i-j)!}p^i(1-p)^{n-i}\\ &=\frac{1}{j!}\sum_{i=0}^n\frac{n!}{(n-i)!(i-j)!}p^i(1-p)^{n-i}\\ &=\frac{1}{j!}\sum_{i=0}^n\frac{(n-j)!}{(n-i)!(i-j)!}\frac{n!}{(n-j)!}p^i(1-p)^{n-i}\\ &=\frac{n!}{(n-j)!j!}\sum_{i=0}^n\binom{n-j}{n-i}p^i(1-p)^{n-i}\\ &=\binom nj\sum_{i=0}^n\binom{n-j}{i-j}p^i(1-p)^{n-i}\\ &=p^j\binom nj \sum_{i=0}^{n-j}\binom{n-j}{i}p^i(1-p)^{n-i-j}\\ &=p^j\binom{n}{j}\\ \therefore ans&=\sum_{j=0}^kj!\begin{Bmatrix}k\\j\end{Bmatrix}p^j\binom nj \end{aligned} \]

\[\because m^n=\sum_{i=0}^m\binom mi i!\begin{Bmatrix}n\\i\end{Bmatrix} \]

通过二项式反演得到：

\[\begin{Bmatrix}n\\m\end{Bmatrix}=\frac{1}{m!}\sum_{i=0}^m(-1)^{m-i}\binom mi i^n \]

代回原式得到：

\[\begin{aligned} ans&=\sum_{j=0}^k\frac{1}{j!}\sum_{i=0}^j(-1)^{i-j}\binom ji i^kj!\binom njp^j\\ &=\sum_{i=0}^k(-1)^ii^k\sum_{j=i}^k\binom nj\binom ji(-p)^j\\ &=\sum_{i=0}^k(-1)^ii^k\sum_{j=i}^k\frac{n!j!}{j!i!(n-j)!(j-i)!}(-p)^j\\ &=\sum_{i=0}^k(-1)^ii^k\binom ni\sum_{j=i}^k \binom {n-i}{j-i}(-p)^j\\ &=\sum_{i=0}^ki^kp^i\binom ni\sum_{j=0}^{k-i}\binom {n-i}{j}(-p)^j\\ \end{aligned} \]

设 \(S(i)=\sum\limits_{j=0}^{k-i}\binom {n-i}j(-p)^j\)，裂开组合数：

\[\begin{aligned} S(i)&=\sum_{j=0}^{k-i}\left( \binom{n-i-1}{j}+\binom{n-i-1}{j-1}\right)(-p)^j\\ &=S(i+1)+\binom{n-i-1}{k-i}(-p)^{k-i}+\sum_{j=0}^{k-i}\binom{n-i-1}{j-1}(-p)^j\\ \end{aligned} \]

令 \(s=S(i+1)+\binom{n-i-1}{k-i}(-p)^{k-i}\),

\[\begin{aligned} S(i)&=s+(-p)\sum_{j=0}^{k-i}\binom{n-i-1}{j-1}(-p)^{j-1}\\ &=s+(-p)\sum_{j=0}^{k-i-1}\binom{n-i-1}j(-p)^j\\ \therefore S(i)&=c+(-p)S(i+1)=(1-p)S(i+1)+\binom{n-i-1}{k-i}(-1)^{k-i}\\(S(k)&=1)\\ \therefore ans&=\sum_{i=0}^ki^k\binom ni p^i S(i) \end{aligned} \]

#include <bits/stdc++.h>
#define pii pair<int,int>
#define int long long
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
bool Mst;
const int Max=1e7+10;
const int mod=998244353;
const int inf=1e9+10;


inline int read(){
	int res=0,v=1;
	char c=getchar();
	while(c<'0'||c>'9'){v=(c=='-'?-1:1);c=getchar();}
	while(c>='0'&&c<='9'){res=(res<<3)+(res<<1)+(c^48);c=getchar();}
	return res*v;
}

int ksm(int a,int b) {
	int ans=1;
	for(;b;b>>=1){
		if(b&1)ans=ans*a%mod;
		a=a*a%mod;
	}
	return ans;
}

int pri[Max],idk[Max],cnt;
int inv[Max];

int k;

void pre(int n){
	inv[1]=idk[1]=1;
	for(int i=2;i<=n;++i){
		if(!idk[i]){
			pri[++cnt]=i;
			idk[i]=ksm(i,k);
			inv[i]=ksm(i,mod-2);
		}
		for(int j=2;j<=cnt&&i*pri[j]<=n;++j){
			idk[i*pri[j]]=idk[i]*idk[pri[j]]%mod;
			inv[i*pri[j]]=inv[i]*inv[pri[j]]%mod;
			if(i%pri[j]==0)break;
		}
	}
}

int n,m;

int frac[Max],Inv[Max];

void init(){
	frac[0]=Inv[0]=1;
	for(int i=1;i<=n;++i)frac[i]=frac[i-1]*i%mod;
	Inv[n]=ksm(frac[n],mod-2);
	for(int i=n-1;i>=1;--i)Inv[i]=Inv[i+1]*(i+1)%mod;
}

int C(int n,int m){
	return frac[n]*Inv[m]%mod*Inv[n-m]%mod;
}

void Work1(){  //n<=k
	int ans=0;
	init();
	pre(n);
	int p=ksm(m,mod-2);
	int Q=((1-p)%mod+mod)%mod;
	int Z=ksm(Q,n);
	Q=ksm(Q,mod-2);
	int z=1;
	for(int i=0;i<=n;++i){
		ans=(ans+C(n,i)*z%mod*Z%mod*idk[i]%mod)%mod;
		z=z*p%mod;
		Z=Z*Q%mod;
	}
	cout << ans << "\n";
}

int S[Max];

void Work2(){
	pre(k);
	int p=ksm(m,mod-2);
	int c=1,x=1;
	S[k]=1;
	for(int i=k-1;i;--i){
		c=c*(n-i-1)%mod*inv[k-i]%mod;
		x=((-1ll*x*p)%mod+mod)%mod;
		S[i]=(S[i+1]*((1-p)+mod)%mod+c*x)%mod;
	}
	c=x=1;
	int ans=0;
	for(int i=1;i<=k;++i){
		x=x*p%mod;
		c=c*(n-i+1)%mod*inv[i]%mod;
		ans=(ans+S[i]*c%mod*idk[i]%mod*x%mod)%mod;
	}	
	cout << ans << "\n";
}

bool Med;
signed main(){
	n=read();
	m=read();
	k=read();
	if(n<=k)Work1();
	else Work2();


	cerr<< "Time: "<<clock()/1000.0 << "s\n";
	cerr<< "Memory: " << (&Mst-&Med)/1000000.0 << "MB\n";
	return 0;
}

P4389 付公主的背包

首先对于一个体积为 \(V\) 物品可以很轻易的求出其生成函数：

\[F_V(x)=\sum_{i\geq 0}[i\%V=0]x^i=\frac{1}{1-x^V} \]

所以最后的答案是 ：

\[[x^m]\prod_{i=1}^nF_{a_i}(x) \]

但是现在是 \(O(n^2\log n)\) 的，但是因为有如下操作：

\[G(x)*F(x)=e^{\ln(G(x))+\ln(F(x))} \]

所以考虑对所有生成函数求 \(\ln\)，设：

\[\begin{aligned} G(x)&=\ln(\frac{1}{1-x^V})\\ &=\ln1-\ln(1-x^V)\\ &=-\ln(1-x^V)\\ &=-\ln(F(x))\\ G'(x)&=\frac{F'(x)}{F(x)}\\ &=\frac{-Vx^{V-1}}{1-x^V}\\ &=-\sum_{i\geq0}Vx^{iV+V-1}\\ G(x)&=-\sum_{i\geq 0}\frac{x^{iV}}{i} \end{aligned} \]

最后求 \(\exp\) 就行。

#include <bits/stdc++.h>
#define pii pair<int,int>
#define ll long long
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
bool Mst;
const int Max=5e5+10;
const int mod=998244353;
const int inf=1e9+10;

template <int mod>
struct modint{

	int val;

	static int norm(const int &x){return x<0?x+mod:x;}

	modint inv()const{
		int a=val,b=mod,u=1,v=0,t;
		while(b>0)t=a/b,swap(a-=t*b,b),swap(u-=t*v,v);
		return modint(u);
	}


	modint():val(0){}
	modint(const int &m):val(norm(m)){}
	modint(const ll &m):val(norm(m%mod)){}
	modint operator -()const{return modint(norm(-val));}
	bool operator ==(const modint &x){return val==x.val;}
	bool operator !=(const modint &x){return val!=x.val;}
	bool operator <=(const modint &x){return val<=x.val;}
	bool operator >=(const modint &x){return val>=x.val;}
	bool operator >(const modint &x){return val>x.val;}
	bool operator <(const modint &x){return val<x.val;}
	modint& operator *=(const modint &x){return val=static_cast<int>(1ll*val*x.val%mod),*this;}
	modint& operator <<=(const modint &x){return val=(1ll*val<<x.val)%mod,*this;}
	modint& operator +=(const modint &x){return val=(1ll*val+x.val)%mod,*this;}
	modint& operator -=(const modint &x){return val=norm(1ll*val-x.val),*this;}
	modint& operator >>=(const modint &x){return val>>=x.val,*this;}
	modint& operator ^=(const modint &x){return val^=x.val,*this;}
	modint operator >>(const modint &x){return modint(*this)>>=x;}
	modint operator <<(const modint &x){return modint(*this)<<=x;}
	modint& operator /=(const modint &x){return *this*=x.inv();}
	modint operator +(const modint &x){return modint(*this)+=x;}
	modint operator -(const modint &x){return modint(*this)-=x;}
	modint operator *(const modint &x){return modint(*this)*=x;}
	modint operator /(const modint &x){return modint(*this)/=x;}
	modint operator ^(const modint &x){return modint(*this)^=x;}
	friend std::ostream& operator<<(std::ostream& os,const modint &a){return os<<a.val;}
	friend std::istream& operator>>(std::istream& is,modint &a){return is>>a.val;}
};

typedef modint<998244353>m98;

inline int read(){
	int res=0,v=1;
	char c=getchar();
	while(c<'0'||c>'9'){v=(c=='-'?-1:1);c=getchar();}
	while(c>='0'&&c<='9'){res=(res<<3)+(res<<1)+(c^48);c=getchar();}
	return res*v;
}

m98 ksm(m98 a,int b){
	m98 ans=1;
	for(;b;b>>=1){
		if(b&1)ans=ans*a;
		a=a*a;
	}
	return ans;
}

m98 inv[Max],frac[Max];

void init(int len){
	inv[0]=frac[0]=1;
	for(int i=1;i<=len;++i)frac[i]=frac[i-1]*i;
	inv[len]=frac[len].inv();for(int i=len-1;i>=1;--i)inv[i]=inv[i+1]*(i+1);
	for(int i=1;i<=len;++i)inv[i]*=frac[i-1];
}

int rev[Max];

void init(int lim,int len){for(int i=0;i<lim;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));}

struct Poly{
	const int g=3;
	const int invg=(mod+1)/3;
	m98 a[Max];
	
	m98& operator [](const int x){return a[x];}
	void out(int n,int x=0){for(int i=0;i<n;++i)cout << a[i].val << ' ';if(x==0)cout << "\n";}
	void inte(int n){for(int i=n-1;i>0;--i)a[i]=a[i-1]*inv[i];a[0]=0;}
	void der(int n){for(int i=0;i<n;++i)a[i]=a[i+1]*(i+1);a[n]=0;}
	void ReadIn(int n){for(int i=0;i<n;++i)a[i].val=read();}
	
	void NTT(int tag,int lim){
		for(int i=0;i<lim;++i)if(i<rev[i])swap(a[i],a[rev[i]]);
		for(int mid=1;mid<lim;mid<<=1){
			m98 W=ksm(tag==1?g:invg,(mod-1)/(mid<<1));
			for(int R=mid<<1,j=0;j<lim;j+=R){
				m98 w=1;
				for(int p=0;p<mid;++p,w=w*W){
					m98 x=a[p+j],y=a[p+j+mid]*w;
					a[p+j]=x+y;
					a[p+j+mid]=x-y;
				}
			}
		}
		if(tag==-1){
			m98 inv=m98(lim).inv();
			for(int i=0;i<lim;++i)a[i]*=inv;
		}
	}
};

void GetInv(Poly &y,Poly &x,int n){
	Poly b[2],c;
	b[0][0]=x[0].inv();
	int pos=0;
	int lim=2,bas=1,len=1;
	while(bas<(n<<1)){
		init(lim,len);
		pos^=1;
		for(int i=0;i<bas;++i)b[pos][i]=b[pos^1][i]<<1;
		for(int i=0;i<lim;++i)c[i]=x[i];
		for(int i=lim>>1;i<=lim;++i)c[i]=0;
		b[pos^1].NTT(1,lim),c.NTT(1,lim) ;
		for(int i=0;i<lim;++i)b[pos^1][i]*=b[pos^1][i]*c[i];
		b[pos^1].NTT(-1,lim) ;
		for(int i=0;i<bas;++i)b[pos][i]-=b[pos^1][i];
		bas<<=1,lim<<=1,++len;
	}
	for(int i=0;i<n;++i)y[i]=b[pos][i];
	for(itn i=0;i<bas;++i)b[0][i]=b[1][i]=c[i]=0;
}

void GetLn(Poly &y,Poly &x,int n){
	Poly t1,t2;
	for(int i=0;i<n;++i)t1[i]=t2[i]=x[i];
	t1.der(n);GetInv(t2,t2,n);
	int lim=1,len=0;
	while(lim<n<<1)lim<<=1,++len;
	init(lim,len);
	t1.NTT(1,lim) ; t2.NTT(1,lim) ;
	for(int i=0;i<lim;++i)t1[i]*=t2[i];
	t1.NTT(-1,lim),t1.inte(n) ;
	for(int i=0;i<n;++i)y[i]=t1[i];
	for(int i=0;i<n;++i)t1[i]=t2[i]=0;
}

void GetExp(Poly &x,Poly &y,int n){
	Poly s,tmp;
	s[0]=1;
	int lim=2,bas=1,len=1;
	while(bas<n<<1){
		for(int i=0;i<bas;++i)tmp[i]=s[i];
		GetLn(tmp,tmp,bas);
		for(int i=0;i<bas;++i)tmp[i]=x[i]-tmp[i];
		tmp[0]+=1;init(lim,len);tmp.NTT(1,lim),s.NTT(1,lim) ;
		for(int i=0;i<lim;++i)s[i]*=tmp[i];
		s.NTT(-1,lim);
		for(int i=bas;i<lim;++i)s[i]=0;
		lim<<=1,bas<<=1,++len;
	}
	for(int i=0;i<n;++i)y[i]=s[i];
	for(int i=0;i<bas;++i)tmp[i]=s[i]=0;
}

int cnt[Max];

Poly x;

bool Med;
signed main(){
	int n=read();
	int m=read();
	init(m+10);
	for(int i=1;i<=n;++i)cnt[read()]++;
	for(int i=1;i<=m;++i){
		for(int j=1;j*i<=m;++j){
			x[i*j]+=inv[j]*cnt[i];
		}
	}
	GetExp(x,x,m+1);
	for(int i=1;i<=m;++i)cout << x[i].val << "\n";

	cerr<< "Time: "<<clock()/1000.0 << "s\n";
	cerr<< "Memory: " << (&Mst-&Med)/1000000.0 << "MB\n";
	return 0;
}
/*

*/

P4841 [集训队作业2013] 城市规划

发现联通这个限制比较强，考虑转换，设 \(g_i\) 表示 \(i\) 个点的简单无向图个数，显然有 \(g_i=2^{\binom i2}\)，设 \(f_i\) 表示 \(i\) 个点的无向连通图个数。这样，我们考虑第一个点所在连通块大小，有如下转移方程：\(g_n=\sum\limits_{i=1}^n\binom{n-1}{i-1}f_{i}g_{n-i}\)，推一下式子得到：

\[g_n=\sum\limits_{i=1}^n\binom{n-1}{i-1}f_{i}g_{n-i} =\sum_{i=1}^n\frac{(n-1)!}{(i-1)!(n-i)!}f_{i}g_{n-i}\\ \frac{g_n}{(n-1)!}=\sum_{i=1}^n\frac{f_{i}}{(i-1)!}.\frac{g_{n-i}}{(n-i)!} \]

设：

\[F=\sum_{i\geq0}\frac{f_i}{i!}x^i\\ H=\sum_{i\geq1}\frac{2^{\binom i2}}{(i-1)!}x^i\\ G=\sum_{i\geq0}\frac{2^{\binom i2}}{i!}x^i \\ \because H=F*G\\ \therefore F=H*G^{-1} \]

#include <bits/stdc++.h>
#define pii pair<int,int>
#define ll long long
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
bool Mst;
const int Max=6e5+10;
const int mod=1004535809;
const int inf=1e9+10;

template <int mod>
struct modint{

	int val;

	static int norm(const int &x){return x<0?x+mod:x;}

	modint inv()const{
		int a=val,b=mod,u=1,v=0,t;
		while(b>0)t=a/b,swap(a-=t*b,b),swap(u-=t*v,v);
		return modint(u);
	}
	
	
	modint():val(0){}
	modint(const int &m):val(norm(m)){}
	modint(const ll &m):val(norm(m%mod)){}
	modint operator -()const{return modint(norm(-val));}
	bool operator ==(const modint &x){return val==x.val;}
	bool operator !=(const modint &x){return val!=x.val;}
	bool operator <=(const modint &x){return val<=x.val;}
	bool operator >=(const modint &x){return val>=x.val;}
	bool operator >(const modint &x){return val>x.val;}
	bool operator <(const modint &x){return val<x.val;}
	modint& operator *=(const modint &x){return val=static_cast<int>(1ll*val*x.val%mod),*this;}
	modint& operator <<=(const modint &x){return val=(1ll*val<<x.val)%mod,*this;}
	modint& operator +=(const modint &x){return val=(1ll*val+x.val)%mod,*this;}
	modint& operator -=(const modint &x){return val=norm(1ll*val-x.val),*this;}
	modint& operator >>=(const modint &x){return val>>=x.val,*this;}
	modint& operator ^=(const modint &x){return val^=x.val,*this;}
	modint operator >>(const modint &x){return modint(*this)>>=x;}
	modint operator <<(const modint &x){return modint(*this)<<=x;}
	modint& operator /=(const modint &x){return *this*=x.inv();}
	modint operator +(const modint &x){return modint(*this)+=x;}
	modint operator -(const modint &x){return modint(*this)-=x;}
	modint operator *(const modint &x){return modint(*this)*=x;}
	modint operator /(const modint &x){return modint(*this)/=x;}
	modint operator ^(const modint &x){return modint(*this)^=x;}
	friend std::ostream& operator<<(std::ostream& os,const modint &a){return os<<a.val;}
	friend std::istream& operator>>(std::istream& is,modint &a){return is>>a.val;}
};

typedef modint<1004535809>m19;

inline int read(){
	int res=0,v=1;
	char c=getchar();
	while(c<'0'||c>'9'){v=(c=='-'?-1:1);c=getchar();}
	while(c>='0'&&c<='9'){res=(res<<3)+(res<<1)+(c^48);c=getchar();}
	return res*v;
}

m19 ksm(m19 a,int b){
	m19 ans=1;
	for(;b;b>>=1){
		if(b&1)ans=ans*a;
		a=a*a;
	}
	return ans;
}

m19 inv[Max],frac[Max],Inv[Max];

void init(int len) {
	Inv[0]=frac[0]=inv[0]=1; 
	for(int i=1;i<=len;++i)frac[i]=frac[i-1]*i;
	Inv[len]=frac[len].inv();for(int i=len-1;i;--i)Inv[i]=Inv[i+1]*(i+1);
	for(int i=1;i<=len;++i)inv[i]=Inv[i]*frac[i-1];
}

int rev[Max];

void init(int lim,int len){for(int i=0;i<lim;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));}

struct Poly{  //导和积都到n，否则到n-1; 
	const int g=3;
	const int invg=(mod+1)/3;
	
	m19 a[Max];
	void out(int n,int x=0){for(int i=0;i<n;++i)cout << a[i].val << " ";if(!x)cout << "\n";}
	void inte(int n){for(int i=n;i>0;--i)a[i]=a[i-1]*inv[i];a[0]=0;}
	void ReadIn(int n){for(int i=0;i<n;++i)a[i].val=read();}
	void der(int n){for(int i=0;i<=n;++i)a[i]=a[i+1]*(i+1);}
	m19& operator [](const int x){return a[x];}
	
	void NTT(int tag,int lim){
		for(int i=0;i<lim;++i)if(i<rev[i])swap(a[i],a[rev[i]]);
		for(int mid=1;mid<lim;mid<<=1){
			m19 W=ksm(tag==1?g:invg,(mod-1)/(mid<<1));
			for(int R=mid<<1,j=0;j<lim;j+=R){
				m19 w=1;
				for(int p=0;p<mid;++p,w=w*W){
					m19 x=a[j+p],y=a[j+p+mid]*w;
					a[j+p]=x+y;
					a[j+p+mid]=x-y;
				}
			}
		}
		if(tag==-1){
			m19 inv=m19(lim).inv();
			for(int i=0;i<lim;++i)a[i]*=inv;
		}
	}
	
}a;

void GetInv(Poly &y,Poly &x,int n){    //第二个为参数，第一个为返回值
	static Poly b[2],c;
	int lim=2,len=1,bas=1;
	int pos=0;
	b[pos][0]=x[0].inv();
	while(bas<(n<<1)){
		init(lim,len);
		pos^=1;
		for(int i=0;i<bas;++i)b[pos][i]=b[pos^1][i]<<1;
		for(int i=0;i<=lim;++i)c[i]=x[i];
		for(int i=lim>>1;i<=lim;++i)c[i]=0;
		b[pos^1].NTT(1,lim);c.NTT(1,lim);
		for(int i=0;i<=lim;++i)b[pos^1][i]=b[pos^1][i]*b[pos^1][i]*c[i];
		b[pos^1].NTT(-1,lim);
		for(int i=0;i<bas;++i)b[pos][i]-=b[pos^1][i],b[pos^1][i]=0;
		bas<<=1,lim<<=1,++len;
	}
	for(int i=0;i<n;++i)y[i]=b[pos][i];
	for(int i=0;i<bas;++i)b[0][i]=b[1][i]=c[i]=0;
}

void GetSqrt(Poly &y,Poly &x,int n){
	
	static Poly t1,t2,s;
	s[0]=1;
	int lim=2,len=1,bas=1;
	while(bas<(n<<1)){
		for(int i=0;i<bas;++i)t1[i]=x[i];
		for(int i=0;i<bas;++i)t2[i]=s[i]<<1;
		GetInv(t2,t2,bas);init(lim,len);
		t2.NTT(1,lim);t1.NTT(1,lim);
		for(int i=0;i<=lim;++i)t1[i]=t1[i]*t2[i];
		t1.NTT(-1,lim);
		for(int i=0;i<bas;++i)s[i]=s[i]/2+t1[i];
		bas<<=1,lim<<=1,++len;
	}
	for(int i=0;i<n;++i)y[i]=s[i];
	for(int i=0;i<bas;++i)t1[i]=t2[i]=s[i]=0;
	
}

void GetLn(Poly &y,Poly&x,int n){
	static Poly t1,t2;
	for(int i=0;i<n;++i)t1[i]=t2[i]=x[i];
	t1.der(n);GetInv(t2,t2,n);
	int lim=1,len=0;
	while(lim<n<<1)lim<<=1,++len;
	init(lim,len);
	t1.NTT(1,lim),t2.NTT(1,lim) ;
	for(int i=0;i<lim;++i)t1[i]*=t2[i];
	t1.NTT(-1,lim);t1.inte(n) ;
	for(int i=0;i<n;++i)y[i]=t1[i];
	for(int i=0;i<lim;++i)t1[i]=t2[i]=0;
}


void GetExp(Poly &y,Poly &x,int n){
	static Poly s,tmp;s[0]=1;
	int bas=1,lim=2,len=1;
	while(bas<n<<1){
		for(int i=0;i<bas;++i)tmp[i]=s[i];
		GetLn(tmp,tmp,bas);
		for(int i=0;i<bas;++i)tmp[i]=x[i]-tmp[i];
		tmp[0]+=1;init(lim,len);tmp.NTT(1,lim),s.NTT(1,lim) ;
		for(int i=0;i<lim;++i)s[i]=s[i]*tmp[i];
		s.NTT(-1,lim) ;
		for(int i=bas;i<lim;++i)s[i]=0;
		lim<<=1,bas<<=1,++len;
	}
	for(int i=0;i<n;++i)y[i]=s[i];
	for(itn i=0;i<bas;++i)s[i]=tmp[i]=0;
}


int Si=0;
void expow(Poly &y,Poly &x,int n,m19 k,int k1){
	static Poly s;
	int z=n+1;
	for(int i=0;i<n;++i)if(x[i]!=0){z=i;break;}
	if(z==n+1||(z!=0&&Si)){
		for(int i=0;i<n;++i)y[i]=0;
		return;
	}
	for(itn i=z;i<n;++i)s[i-z]=x[i]/x[z];GetLn(s,s,n-z);
	for(int i=0;i<n-z;++i)s[i]*=k;GetExp(s,s,n-z);
	for(int i=0;i<n-z;++i)s[i]=s[i]*ksm(x[z],k1);
	for(int i=0;i<n;++i)y[i]=0;
	for(int i=0;1ll*z*k.val+i<n;++i)y[i+z*k.val]=s[i];
	for(int i=0;i<n;++i)s[i]=0; 
}


Poly b;

m19 k;
ll k1;

void In(){
	k=k1=0;
	char c=getchar();
	while(c<'0'||c>'9')c=getchar();
	while(c>='0'&&c<='9'){Si=Si|(k*10+(c^48)<k),k=(k<<3)+(k<<1)+(c^48),k1=(k1*10%(mod-1)+(c^48))%(mod-1),c=getchar();}
}

Poly H,G,F;

bool Med;
signed main(){
	int n=read();
	init(n<<2);
	for(int i=1;i<=n;++i)H[i]=ksm(2,(1ll*i*(i-1)/2)%(mod-1))*Inv[i-1];
	for(int i=0;i<=n;++i)G[i]=ksm(2,(1ll*i*(i-1)/2)%(mod-1))*Inv[i];
	GetInv(G,G,n+1);
	int lim=1,len=0;
	while(lim<=n<<1)lim<<=1,++len;
	init(lim,len);
	H.NTT(1,lim),G.NTT(1,lim);
	for(int i=0;i<lim;++i)F[i]=G[i]*H[i];
	F.NTT(-1,lim);
	cout << F[n]*frac[n-1] << "\n";

	cerr<< "Time: "<<clock()/1000.0 << "s\n";
	cerr<< "Memory: " << (&Mst-&Med)/1000000.0 << "MB\n";
	return 0;
}

P5748 集合划分计数

考虑先对单个集合生成函数 \(F(x)=e^x-1=\sum\limits_{i\geq1}\frac{x^i}{i!}\) 然后现在令答案的生成函数是 \(G(x)\)，显然有：

\[G(x)=\sum_{i=0}^n\frac{F^i(x)}{i!}=e^{F(x)} \]

#include <bits/stdc++.h>
#define pii pair<int,int>
#define ll long long
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
bool Mst;
const int Max=6e5+10;
const int mod=998244353;
const int inf=1e9+10;

template <int mod>
struct modint{

	int val;

	static int norm(const int &x){return x<0?x+mod:x;}

	modint inv()const{
		int a=val,b=mod,u=1,v=0,t;
		while(b>0)t=a/b,swap(a-=t*b,b),swap(u-=t*v,v);
		return modint(u);
	}
	
	
	modint():val(0){}
	modint(const int &m):val(norm(m)){}
	modint(const ll &m):val(norm(m%mod)){}
	modint operator -()const{return modint(norm(-val));}
	bool operator ==(const modint &x){return val==x.val;}
	bool operator !=(const modint &x){return val!=x.val;}
	bool operator <=(const modint &x){return val<=x.val;}
	bool operator >=(const modint &x){return val>=x.val;}
	bool operator >(const modint &x){return val>x.val;}
	bool operator <(const modint &x){return val<x.val;}
	modint& operator *=(const modint &x){return val=static_cast<int>(1ll*val*x.val%mod),*this;}
	modint& operator <<=(const modint &x){return val=(1ll*val<<x.val)%mod,*this;}
	modint& operator +=(const modint &x){return val=(1ll*val+x.val)%mod,*this;}
	modint& operator -=(const modint &x){return val=norm(1ll*val-x.val),*this;}
	modint& operator >>=(const modint &x){return val>>=x.val,*this;}
	modint& operator ^=(const modint &x){return val^=x.val,*this;}
	modint operator >>(const modint &x){return modint(*this)>>=x;}
	modint operator <<(const modint &x){return modint(*this)<<=x;}
	modint& operator /=(const modint &x){return *this*=x.inv();}
	modint operator +(const modint &x){return modint(*this)+=x;}
	modint operator -(const modint &x){return modint(*this)-=x;}
	modint operator *(const modint &x){return modint(*this)*=x;}
	modint operator /(const modint &x){return modint(*this)/=x;}
	modint operator ^(const modint &x){return modint(*this)^=x;}
	friend std::ostream& operator<<(std::ostream& os,const modint &a){return os<<a.val;}
	friend std::istream& operator>>(std::istream& is,modint &a){return is>>a.val;}
};

typedef modint<998244353>m98;

inline int read(){
	int res=0,v=1;
	char c=getchar();
	while(c<'0'||c>'9'){v=(c=='-'?-1:1);c=getchar();}
	while(c>='0'&&c<='9'){res=(res<<3)+(res<<1)+(c^48);c=getchar();}
	return res*v;
}

m98 ksm(m98 a,int b){
	m98 ans=1;
	for(;b;b>>=1){
		if(b&1)ans=ans*a;
		a=a*a;
	}
	return ans;
}

m98 inv[Max],frac[Max],Inv[Max];

void init(int len) {
	Inv[0]=frac[0]=inv[0]=1; 
	for(int i=1;i<=len;++i)frac[i]=frac[i-1]*i;
	Inv[len]=frac[len].inv();for(int i=len-1;i;--i)Inv[i]=Inv[i+1]*(i+1);
	for(int i=1;i<=len;++i)inv[i]=Inv[i]*frac[i-1];
}

int rev[Max];

void init(int lim,int len){for(int i=0;i<lim;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));}

struct Poly{  //导和积都到n，否则到n-1; 
	const int g=3;
	const int invg=(mod+1)/3;
	
	m98 a[Max];
	void out(int n,int x=0){for(int i=0;i<n;++i)cout << a[i].val << " ";if(!x)cout << "\n";}
	void inte(int n){for(int i=n;i>0;--i)a[i]=a[i-1]*inv[i];a[0]=0;}
	void ReadIn(int n){for(int i=0;i<n;++i)a[i].val=read();}
	void der(int n){for(int i=0;i<=n;++i)a[i]=a[i+1]*(i+1);}
	m98& operator [](const int x){return a[x];}
	
	void NTT(int tag,int lim){
		for(int i=0;i<lim;++i)if(i<rev[i])swap(a[i],a[rev[i]]);
		for(int mid=1;mid<lim;mid<<=1){
			m98 W=ksm(tag==1?g:invg,(mod-1)/(mid<<1));
			for(int R=mid<<1,j=0;j<lim;j+=R){
				m98 w=1;
				for(int p=0;p<mid;++p,w=w*W){
					m98 x=a[j+p],y=a[j+p+mid]*w;
					a[j+p]=x+y;
					a[j+p+mid]=x-y;
				}
			}
		}
		if(tag==-1){
			m98 inv=m98(lim).inv();
			for(int i=0;i<lim;++i)a[i]*=inv;
		}
	}
	
}a;

void GetInv(Poly &y,Poly &x,int n){    //第二个为参数，第一个为返回值
	static Poly b[2],c;
	int lim=2,len=1,bas=1;
	int pos=0;
	b[pos][0]=x[0].inv();
	while(bas<(n<<1)){
		init(lim,len);
		pos^=1;
		for(int i=0;i<bas;++i)b[pos][i]=b[pos^1][i]<<1;
		for(int i=0;i<=lim;++i)c[i]=x[i];
		for(int i=lim>>1;i<=lim;++i)c[i]=0;
		b[pos^1].NTT(1,lim);c.NTT(1,lim);
		for(int i=0;i<=lim;++i)b[pos^1][i]=b[pos^1][i]*b[pos^1][i]*c[i];
		b[pos^1].NTT(-1,lim);
		for(int i=0;i<bas;++i)b[pos][i]-=b[pos^1][i],b[pos^1][i]=0;
		bas<<=1,lim<<=1,++len;
	}
	for(int i=0;i<n;++i)y[i]=b[pos][i];
	for(int i=0;i<bas;++i)b[0][i]=b[1][i]=c[i]=0;
}

void GetSqrt(Poly &y,Poly &x,int n){
	
	static Poly t1,t2,s;
	s[0]=1;
	int lim=2,len=1,bas=1;
	while(bas<(n<<1)){
		for(int i=0;i<bas;++i)t1[i]=x[i];
		for(int i=0;i<bas;++i)t2[i]=s[i]<<1;
		GetInv(t2,t2,bas);init(lim,len);
		t2.NTT(1,lim);t1.NTT(1,lim);
		for(int i=0;i<=lim;++i)t1[i]=t1[i]*t2[i];
		t1.NTT(-1,lim);
		for(int i=0;i<bas;++i)s[i]=s[i]/2+t1[i];
		bas<<=1,lim<<=1,++len;
	}
	for(int i=0;i<n;++i)y[i]=s[i];
	for(int i=0;i<bas;++i)t1[i]=t2[i]=s[i]=0;
	
}

void GetLn(Poly &y,Poly&x,int n){
	static Poly t1,t2;
	for(int i=0;i<n;++i)t1[i]=t2[i]=x[i];
	t1.der(n);GetInv(t2,t2,n);
	int lim=1,len=0;
	while(lim<n<<1)lim<<=1,++len;
	init(lim,len);
	t1.NTT(1,lim),t2.NTT(1,lim) ;
	for(int i=0;i<lim;++i)t1[i]*=t2[i];
	t1.NTT(-1,lim);t1.inte(n) ;
	for(int i=0;i<n;++i)y[i]=t1[i];
	for(int i=0;i<lim;++i)t1[i]=t2[i]=0;
}


void GetExp(Poly &y,Poly &x,int n){
	static Poly s,tmp;s[0]=1;
	int bas=1,lim=2,len=1;
	while(bas<n<<1){
		for(int i=0;i<bas;++i)tmp[i]=s[i];
		GetLn(tmp,tmp,bas);
		for(int i=0;i<bas;++i)tmp[i]=x[i]-tmp[i];
		tmp[0]+=1;init(lim,len);tmp.NTT(1,lim),s.NTT(1,lim) ;
		for(int i=0;i<lim;++i)s[i]=s[i]*tmp[i];
		s.NTT(-1,lim) ;
		for(int i=bas;i<lim;++i)s[i]=0;
		lim<<=1,bas<<=1,++len;
	}
	for(int i=0;i<n;++i)y[i]=s[i];
	for(itn i=0;i<bas;++i)s[i]=tmp[i]=0;
}


int Si=0;
void expow(Poly &y,Poly &x,int n,m98 k,int k1){
	static Poly s;
	int z=n+1;
	for(int i=0;i<n;++i)if(x[i]!=0){z=i;break;}
	if(z==n+1||(z!=0&&Si)){
		for(int i=0;i<n;++i)y[i]=0;
		return;
	}
	for(itn i=z;i<n;++i)s[i-z]=x[i]/x[z];GetLn(s,s,n-z);
	for(int i=0;i<n-z;++i)s[i]*=k;GetExp(s,s,n-z);
	for(int i=0;i<n-z;++i)s[i]=s[i]*ksm(x[z],k1);
	for(int i=0;i<n;++i)y[i]=0;
	for(int i=0;1ll*z*k.val+i<n;++i)y[i+z*k.val]=s[i];
	for(int i=0;i<n;++i)s[i]=0; 
}


Poly b;

m98 k;
ll k1;

void In(){
	k=k1=0;
	char c=getchar();
	while(c<'0'||c>'9')c=getchar();
	while(c>='0'&&c<='9'){Si=Si|(k*10+(c^48)<k),k=(k<<3)+(k<<1)+(c^48),k1=(k1*10%(mod-1)+(c^48))%(mod-1),c=getchar();}
}

Poly F;

bool Med;
signed main(){
	int n=1e5;
	init(200000);
	for(int i=1;i<=n;++i)F[i]=Inv[i];
	GetExp(F,F,n+1);
	m98 p=1;
	for(int i=1;i<=n;++i,p=p*i)F[i]*=p;
	int T=read();
	while(T--)cout << F[read()].val << "\n";
	
	cerr<< "Time: "<<clock()/1000.0 << "s\n";
	cerr<< "Memory: " << (&Mst-&Med)/1000000.0 << "MB\n";
	return 0;
}

P5824 十二重计数法

球与盒子互不相同

无限制

显然是 \(m^n\)。

每个盒子最多一个球

发现 \(n>m\) 时答案为 0，从组合上讲答案是 \(n!\binom mn\)。

每个盒子至少一个球

从 EGF 的角度出发，每个盒子的 EGF 是 \((e^x-1)\)，故答案是 \([x^n](e^x-1)^m\)，然后直接多项式快速幂？

球要区分，盒子不区分

无限制

是第二类斯特林数，答案为 \(\sum\limits_{i=1}^m\begin{Bmatrix}n\\i\end{Bmatrix}\)。

每个盒子最多一个球

答案为 \([n\leq m]\)

每个盒子至少一个球

第二类斯特林数 \(\begin{Bmatrix}n\\m\end{Bmatrix}\)。

球不区分，盒子要区分

无限制

考虑插板法 \(\binom{n+m-1}{m-1}\)。

每个盒子最多一个球

答案为选 \(n\) 个有球的盒子 \(\binom mn\)。

每个盒子至少一个球

先给每个盒子放一个，然后插板法 \(\binom{n-1}{m-1}\)。

球与盒子都不区分

无限制

发现这个就是P4389 付公主的背包，只不过重量只有 \(1\) 。

每个盒子最多一个球

答案是 \([n\leq m]\)。

每个盒子至少一个球

同样，先保住下限，然后跑无限制。

#include <bits/stdc++.h>
#define pii pair<int,int>
#define ll long long
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
bool Mst;
const int Max=1e6+10;
const int mod=998244353;
const int inf=1e9+10;

template <int mod>
struct modint{

	int val;

	static int norm(const int &x){return x<0?x+mod:x;}

	modint inv()const{
		int a=val,b=mod,u=1,v=0,t;
		while(b>0)t=a/b,swap(a-=t*b,b),swap(u-=t*v,v);
		return modint(u);
	}
	
	
	modint():val(0){}
	modint(const int &m):val(norm(m)){}
	modint(const ll &m):val(norm(m%mod)){}
	modint operator -()const{return modint(norm(-val));}
	bool operator ==(const modint &x){return val==x.val;}
	bool operator !=(const modint &x){return val!=x.val;}
	bool operator <=(const modint &x){return val<=x.val;}
	bool operator >=(const modint &x){return val>=x.val;}
	bool operator >(const modint &x){return val>x.val;}
	bool operator <(const modint &x){return val<x.val;}
	modint& operator *=(const modint &x){return val=static_cast<int>(1ll*val*x.val%mod),*this;}
	modint& operator <<=(const modint &x){return val=(1ll*val<<x.val)%mod,*this;}
	modint& operator +=(const modint &x){return val=(1ll*val+x.val)%mod,*this;}
	modint& operator -=(const modint &x){return val=norm(1ll*val-x.val),*this;}
	modint& operator >>=(const modint &x){return val>>=x.val,*this;}
	modint& operator ^=(const modint &x){return val^=x.val,*this;}
	modint operator >>(const modint &x){return modint(*this)>>=x;}
	modint operator <<(const modint &x){return modint(*this)<<=x;}
	modint& operator /=(const modint &x){return *this*=x.inv();}
	modint operator +(const modint &x){return modint(*this)+=x;}
	modint operator -(const modint &x){return modint(*this)-=x;}
	modint operator *(const modint &x){return modint(*this)*=x;}
	modint operator /(const modint &x){return modint(*this)/=x;}
	modint operator ^(const modint &x){return modint(*this)^=x;}
	friend std::ostream& operator<<(std::ostream& os,const modint &a){return os<<a.val;}
	friend std::istream& operator>>(std::istream& is,modint &a){return is>>a.val;}
};

typedef modint<998244353>m98;

inline int read(){
	int res=0,v=1;
	char c=getchar();
	while(c<'0'||c>'9'){v=(c=='-'?-1:1);c=getchar();}
	while(c>='0'&&c<='9'){res=(res<<3)+(res<<1)+(c^48);c=getchar();}
	return res*v;
}

m98 ksm(m98 a,int b){
	m98 ans=1;
	for(;b;b>>=1){
		if(b&1)ans=ans*a;
		a=a*a;
	}
	return ans;
}

m98 inv[Max],frac[Max],Inv[Max];

void init(int len) {
	Inv[0]=frac[0]=inv[0]=1; 
	for(int i=1;i<=len;++i)frac[i]=frac[i-1]*i;
	Inv[len]=frac[len].inv();for(int i=len-1;i;--i)Inv[i]=Inv[i+1]*(i+1);
	for(int i=1;i<=len;++i)inv[i]=Inv[i]*frac[i-1];
}

m98 C(int n,int m){
	return frac[n]*Inv[m]*Inv[n-m];
}

int rev[Max];

void init(int lim,int len){for(int i=0;i<lim;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));}

struct Poly{  //导和积都到n，否则到n-1; 
	const int g=3;
	const int invg=(mod+1)/3;
	
	m98 a[Max];
	void out(int n,int x=0){for(int i=0;i<n;++i)cout << a[i].val << " ";if(!x)cout << "\n";}
	void inte(int n){for(int i=n;i>0;--i)a[i]=a[i-1]*inv[i];a[0]=0;}
	void ReadIn(int n){for(int i=0;i<n;++i)a[i].val=read();}
	void der(int n){for(int i=0;i<=n;++i)a[i]=a[i+1]*(i+1);}
	m98& operator [](const int x){return a[x];}
	
	void NTT(int tag,int lim){
		for(int i=0;i<lim;++i)if(i<rev[i])swap(a[i],a[rev[i]]);
		for(int mid=1;mid<lim;mid<<=1){
			m98 W=ksm(tag==1?g:invg,(mod-1)/(mid<<1));
			for(int R=mid<<1,j=0;j<lim;j+=R){
				m98 w=1;
				for(int p=0;p<mid;++p,w=w*W){
					m98 x=a[j+p],y=a[j+p+mid]*w;
					a[j+p]=x+y;
					a[j+p+mid]=x-y;
				}
			}
		}
		if(tag==-1){
			m98 inv=m98(lim).inv();
			for(int i=0;i<lim;++i)a[i]*=inv;
		}
	}
	
};

void GetInv(Poly &y,Poly &x,int n){    //第二个为参数，第一个为返回值
	static Poly b[2],c;
	int lim=2,len=1,bas=1;
	int pos=0;
	b[pos][0]=x[0].inv();
	while(bas<(n<<1)){
		init(lim,len);
		pos^=1;
		for(int i=0;i<bas;++i)b[pos][i]=b[pos^1][i]<<1;
		for(int i=0;i<=lim;++i)c[i]=x[i];
		for(int i=lim>>1;i<=lim;++i)c[i]=0;
		b[pos^1].NTT(1,lim);c.NTT(1,lim);
		for(int i=0;i<=lim;++i)b[pos^1][i]=b[pos^1][i]*b[pos^1][i]*c[i];
		b[pos^1].NTT(-1,lim);
		for(int i=0;i<bas;++i)b[pos][i]-=b[pos^1][i],b[pos^1][i]=0;
		bas<<=1,lim<<=1,++len;
	}
	for(int i=0;i<n;++i)y[i]=b[pos][i];
	for(int i=0;i<bas;++i)b[0][i]=b[1][i]=c[i]=0;
}

void GetSqrt(Poly &y,Poly &x,int n){
	
	static Poly t1,t2,s;
	s[0]=1;
	int lim=2,len=1,bas=1;
	while(bas<(n<<1)){
		for(int i=0;i<bas;++i)t1[i]=x[i];
		for(int i=0;i<bas;++i)t2[i]=s[i]<<1;
		GetInv(t2,t2,bas);init(lim,len);
		t2.NTT(1,lim);t1.NTT(1,lim);
		for(int i=0;i<=lim;++i)t1[i]=t1[i]*t2[i];
		t1.NTT(-1,lim);
		for(int i=0;i<bas;++i)s[i]=s[i]/2+t1[i];
		bas<<=1,lim<<=1,++len;
	}
	for(int i=0;i<n;++i)y[i]=s[i];
	for(int i=0;i<bas;++i)t1[i]=t2[i]=s[i]=0;
	
}

void GetLn(Poly &y,Poly&x,int n){
	static Poly t1,t2;
	for(int i=0;i<n;++i)t1[i]=t2[i]=x[i];
	t1.der(n);GetInv(t2,t2,n);
	int lim=1,len=0;
	while(lim<n<<1)lim<<=1,++len;
	init(lim,len);
	t1.NTT(1,lim),t2.NTT(1,lim) ;
	for(int i=0;i<lim;++i)t1[i]*=t2[i];
	t1.NTT(-1,lim);t1.inte(n) ;
	for(int i=0;i<n;++i)y[i]=t1[i];
	for(int i=0;i<lim;++i)t1[i]=t2[i]=0;
}

void GetExp(Poly &y,Poly &x,int n){
	static Poly s,tmp;s[0]=1;
	int bas=1,lim=2,len=1;
	while(bas<n<<1){
		for(int i=0;i<bas;++i)tmp[i]=s[i];
		GetLn(tmp,tmp,bas);
		for(int i=0;i<bas;++i)tmp[i]=x[i]-tmp[i];
		tmp[0]+=1;init(lim,len);tmp.NTT(1,lim),s.NTT(1,lim) ;
		for(int i=0;i<lim;++i)s[i]=s[i]*tmp[i];
		s.NTT(-1,lim) ;
		for(int i=bas;i<lim;++i)s[i]=0;
		lim<<=1,bas<<=1,++len;
	}
	for(int i=0;i<n;++i)y[i]=s[i];
	for(itn i=0;i<bas;++i)s[i]=tmp[i]=0;
}

int Si=0;
void expow(Poly &y,Poly &x,int n,m98 k,int k1){
	static Poly s;
	int z=n+1;
	for(int i=0;i<n;++i)if(x[i]!=0){z=i;break;}
	if(z==n+1||(z!=0&&Si)){
		for(int i=0;i<n;++i)y[i]=0;
		return;
	}
	for(itn i=z;i<n;++i)s[i-z]=x[i]/x[z];GetLn(s,s,n-z);
	for(int i=0;i<n-z;++i)s[i]*=k;GetExp(s,s,n-z);
	for(int i=0;i<n-z;++i)s[i]=s[i]*ksm(x[z],k1);
	for(int i=0;i<n;++i)y[i]=0;
	for(int i=0;1ll*z*k.val+i<n;++i)y[i+z*k.val]=s[i];
	for(int i=0;i<n;++i)s[i]=0; 
}



m98 k;
ll k1;

void In(){
	k=k1=0;
	char c=getchar();
	while(c<'0'||c>'9')c=getchar();
	while(c>='0'&&c<='9'){Si=Si|(k*10+(c^48)<k),k=(k<<3)+(k<<1)+(c^48),k1=(k1*10%(mod-1)+(c^48))%(mod-1),c=getchar();}
}


void Solve1(int n,int m){cout << ksm(m,n)<< "\n";}
void Solve2(int n,int m){if(n>m)cout << "0\n";else cout << frac[n]*C(m,n) << "\n";}
void Solve3(int n,int m){
	Poly b;b[0]=0;
	for(int i=1;i<=n;++i)b[i]=Inv[i];
	expow(b,b,n+1,m,m);
	cout << b[n]*frac[n] << '\n';
}

Poly a,b;
void Solve4(int n,int m){
	for(int i=0;i<=m;++i){
		a[i]=ksm(i,n)*Inv[i];
		b[i]=Inv[i]*(i%2==0?1:-1);
	}
	int lim=1,len=0;
	while(lim<=m+m)lim<<=1,++len;
	init(lim,len);
	a.NTT(1,lim),b.NTT(1,lim) ;
	for(int i=0;i<lim;++i)a[i]*=b[i];
	a.NTT(-1,lim);
	m98 ans=0;
	for(int i=1;i<=m;++i)ans+=a[i];
	cout << ans << "\n";
}

void Solve5(int n,int m){cout <<(n<=m)<< "\n";}
void Solve6(int n,int m){cout << a[m]<< "\n";}
void Solve7(int n,int m){cout << C(n+m-1,m-1)<< "\n";}
void Solve8(int n,int m){cout << C(m,n)<< "\n";}
void Solve9(int n,int m){if(n<m)cout << 0<< "\n";else cout <<  C(n-1,m-1)<< "\n";}

Poly x;

void Solve10(int n,int m){
	for(int i=1;i<=m;++i){
		for(int j=1;j*i<=n;++j){
			x[i*j]-=inv[j];
		}
	}
	GetExp(x,x,n+1);
	GetInv(x,x,n+1);
	cout << x[n] << "\n";
}

void Solve11(int n,int m){cout << (n<=m) << "\n";}
void Solve12(int n,int m){
	if(n<m)cout << 0<< "\n";
	else cout << x[n-m]<< "\n";
}


bool Med;
signed main(){
	int n=read(),m=read();
	init(500000);
	Solve1(n,m);Solve2(n,m);Solve3(n,m);
	Solve4(n,m);Solve5(n,m);Solve6(n,m);
	Solve7(n,m);Solve8(n,m);Solve9(n,m);
	Solve10(n,m);Solve11(n,m);Solve12(n,m);
	cerr<< "Time: "<<clock()/1000.0 << "s\n";
	cerr<< "Memory: " << (&Mst-&Med)/1000000.0 << "MB\n";
	return 0;
}