数学板子

杨辉三角

long long C[1005][1005]{};
void pre()
{
    C[0][0]=1;
    for(int i=1;i<=n;i++)
    {
        C[i][0]=1;
        for(int j=1;j<=i;j++)
        {
		C[i][j]=c[i-1][j-1]+C[i-1][j];
         C[i][j]%=mod;
        }
	}
}

线性筛素数/积性函数

#define N 5005
int n,pri[N],phi[N],mob[N],d[N],num[N],s[N],sp[N],cnt;
bool np[N];//not prime
//pri质数 phi欧拉函数 mob莫比乌斯函数 约数个数d 约束和s 辅助函数num(最小质因子幂次) sp最小质因子等比数列和
void init()
{
    np[0]=np[1]=true;
	mob[1]=1;phi[1]=1;d[1]=1;
    for(int i=2;i<=n;i++)
    {
        if(!np[i])
        {
			pri[++cnt]=i;
            phi[i]=i-1;
            mob[i]=-1;
            d[i]=2;
            num[i]=1;
            s[i]=i+1;
            sp[i]=i+1;
        }
        for(int j=1;j<=cnt&&i*pri[j]<=n;j++)
        {
            np[i*pri[j]]==true;
            if(i%pri[j]==0)
            {
                phi[i*pri[j]]=phi[i]*pri[j];
                mob[i*pri[j]]=0;
                d[i*pri[j]]=d[i]/(num[i]+1)*(num[i]+2);
                num[i*pri[j]] = num[i]+1;
                s[i*pri[j]]=s[i]/sp[i]*(sp[i]*pri[j]+1);
                sp[i*pri[j]] = sp[i]*pri[j]+1;
                break;
		   }
           	    phi[i*pri[j]]=phi[i]*(pri[j]-1);
                mob[i*pri[j]]=-mob[i];
                d[i*pri[j]]=d[i]*2;
                num[i*pri[j]] = 1;
                s[i*pri[j]]=s[i]*(pri[j]+1);
                sp[i*pri[j]] = pri[j]+1;
		}
	}
}

杜教筛---phi和mob的前缀和

const long long N = 1e7
bool np[N+5];
vector<int> pri;
long long mob[N+5],smob[N+5],phi[N+5],sphi[N+5];
unordered_map<long long,long long> usmob,usphi;

void init(){
	mob[1]=1;phi[1]=1;
	for(int i=2;i<=N;++i){
		if(!np[i]){
			pri.push_back(i);
			mob[i]=-1;phi[i]=i-1;
		}
	for(int p:pri){
		if(p*i>N)break;
		np[p*i]=1;
		if(i%p==0){
		mob[p*i]=0;
		phi[p*i]=phi[i]*p;
		break;
		}
		mob[p*i]=mob[p]*mob[i];
		phi[p*i]=phi[p]*phi[i];
	}
	}
	for(int i=1;i<=N;++i){
	smob[i]=smob[i-1]+mob[i];
	sphi[i]=sphi[i-1]+phi[i];
	}
}
long long sum_mob(long long x){
	if(x<=N)return smob[x];
	if(usmob.count(x))return usmob[x];
	long long ans=1;
	for(long long l=2,r;l<=x;l=r+1){
	r=x/(x/l);
	ans-=(r-l+1)*sum_mob(x/l);
	}
	usmob[x]=ans;
	return ans;
}
long long sum_phi(long long x){
	if(x<=N)return sphi[x];
	if(usphi.count(x))return usphi[x];
	long long ans=x*(x+1)/2;
	for(long long l=2,r;l<=x;l=r+1){
		r=x/(x/l);
		ans-=(r-l+1)*sum_phi(x/l);
	}
	usphi[x]=ans;
	return ans;
}

Miller-Rabin素性测试和pollard-rho质因数分解

#include<bits/stdc++.h>
using namespace std;
typedef longlong ll;
ll T,maxm;
ll qpow(ll x,ll exp,ll mod){
	ll re=1;
	while(exp){
	if(exp&1)re=(__int128)re*x%mod;
	x=(__int128)x*x%mod;
	exp>>=1;
	}
 	return re;
}
bool is_prime(llp){
	if(p<3)return p==2;
	if(!(p&1))return false;
	const static ll bse[]={2,325,9375,28178,450775,9780504,1795265022};
	ll d=p-1,r=0;
	while(!(d&1)){d>>=1;r++;}
	for(ll a:bse){
	ll v=qpow(a,d,p);
	if(v<=1||v==p-1)continue;
	for(inti=1;i<=r;++i){
	v=(__int128)v*v%p;
	if(v==p-1&&i!=r){v=1;break;}
	if(v<=1)return false;
	}
	if(v!=1)return false;
	}
	return true;
}
ll pollard_rho(ll p){
	if(p==4) return 2;
	if(is_prime(p))return p;
	while(1){
		ll c=1ll*rand()%(p-1)+1;
		auto f=[=](ll x){
		return((__int128)x*x+c)%p;
		};
		ll t=0,r=0,mul=1,bd=1,tmp;
		while(1){
			for(ll stp=1;stp<=bd;++stp){
				t=f(t);
				mul=(__int128)mul*abs(t-r)%p;
        		if(!(stp%127)){
				tmp=__gcd(mul,p);
					if(tmp>1)return tmp;
				}
			}
			tmp=__gcd(mul,p);
			if(tmp>1)return tmp;
			bd<<=1;
			r=t;mul=1;
		}
	}
}
void findfac(ll n){
	if(n==1||n<=maxm)return;
	if(is_prime(n)){maxm=max(n,maxm);return;}
	ll v=pollard_rho(n);
	while(n%v==0)n/=v;
	findfac(v);
	findfac(n);
}
        

exgcd/扩展欧几里得

long long exgcd(long long a, long long b, long long& x, long long& y)
{
	if (!b) return a;
	long long d = exgcd(b, a % b, y, x);
	y -= a / b * x;
	return d;
}

扩展中国剩余定理

#include<bits/stdc++.h>
using  namespace std;
typedef longlong ll;
typedef __int128 lll;
const int N=1e5+5;
int n;
ll a[N],m[N];
ll exgcd(ll a,ll b,ll& x,ll& y){
	if(b==0){
		x=1;y=0;
	return a;
	}
	ll re=exgcd(b,a%b,x,y),tmp;
	tmp=x;x=y;y=tmp-(a/b)*y;
	return re;
}
ll exCRT(int n,ll* a,ll* m){
	ll ans=a[1],mod=m[1];
	for(inti=2;i<=n;++i){//两两合并
		ll x,y,d=exgcd(mod,m[i],x,y),tmp;
		if((a[i]-ans)%d)return-1;
		lll px=(lll)(a[i]-ans)/d*x;
		tmp=m[i]/d;
		x=(px%tmp+tmp)%tmp;
		ans=((lll)x*mod+ans)%(mod/d*m[i]);
		mod=mod/d*m[i];
     }
	return ans;
}

快速幂

long long qpow(long long a,long long b)
{
    long long ans=1;
    while(b!=0)
    {
       if(b&1) ans*=a,ans%=mod;
        a*=a;
        a%=mod;
        b>>=1;
	}
    return ans;
}

逆元

费马小定理（结合快速幂）

long long get_inv1(long long x)
{
    return qpow(x,mod-2)%mod;
}

扩展欧几里得

long long get_inv2(long long a,long long b)
{
    //求a的逆元
    long long d,x,y;
    d=exgcd(a,b,x,y);
    return d==1?(x+b)%b:-1;
}

线性求逆元

long long inv[100005]{};
void pre()
{
inv[1] = 1;
for (int i = 2; i <= n; ++i) {
  inv[i] = (long long)(mod - mod / i) * inv[mod % i] % mod;
}
}

线性求任意n数逆元

long long a[100005],s[100005],sv[100005]{},inv[100005];
void pre()
{
s[0] = 1;
for (int i = 1; i <= n; ++i) s[i] = s[i - 1] * a[i] % mod;
sv[n] = qpow(s[n], mod - 2);
// 当然这里也可能用 exgcd 来求逆元.
for (int i = n; i >= 1; --i) sv[i - 1] = sv[i] * a[i] % mod;
for (int i = 1; i <= n; ++i) inv[i] = sv[i] * s[i - 1] % mod;
}

欧拉函数(单个数的)

#include <cmath>

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

欧拉函数（多个数）

constexpr int N = 1E7;
constexpr int P = 1000003;

bool isprime[N + 1];
int phi[N + 1];
std::vector<int> primes;

std::fill(isprime + 2, isprime + N + 1, true);
phi[1] = 1;
for (int i = 2; i <= N; i++) {
    if (isprime[i]) {
        primes.push_back(i);
        phi[i] = i - 1;
    }
    for (auto p : primes) {
        if (i * p > N) {
            break;
        }
        isprime[i * p] = false;
        if (i % p == 0) {
            phi[i * p] = phi[i] * p;
            break;
        }
        phi[i * p] = phi[i] * (p - 1);
    }
}

矩阵快速幂

struct mat {
  LL a[sz][sz];

  mat() { memset(a, 0, sizeof a); }

  mat operator-(const mat& T) const {
    mat res;
    for (int i = 0; i < sz; ++i)
      for (int j = 0; j < sz; ++j) {
        res.a[i][j] = (a[i][j] - T.a[i][j]) % MOD;
      }
    return res;
  }

  mat operator+(const mat& T) const {
    mat res;
    for (int i = 0; i < sz; ++i)
      for (int j = 0; j < sz; ++j) {
        res.a[i][j] = (a[i][j] + T.a[i][j]) % MOD;
      }
    return res;
  }

  mat operator*(const mat& T) const {
    mat res;
    int r;
    for (int i = 0; i < sz; ++i)
      for (int k = 0; k < sz; ++k) {
        r = a[i][k];
        for (int j = 0; j < sz; ++j)
          res.a[i][j] += T.a[k][j] * r, res.a[i][j] %= MOD;
      }
    return res;
  }

  mat operator^(LL x) const {
    mat res, bas;
    for (int i = 0; i < sz; ++i) res.a[i][i] = 1;
    for (int i = 0; i < sz; ++i)
      for (int j = 0; j < sz; ++j) bas.a[i][j] = a[i][j] % MOD;
    while (x) {
      if (x & 1) res = res * bas;
      bas = bas * bas;
      x >>= 1;
    }
    return res;
  }
};

高斯消元

double ans[N],a[N][N];
int n;
int gauss(){
	int dim=0;
	for(int i=1;i<=n;i++){
		int r=dim+1;
		int t=r;
		for(int j=t+1;j<=n;j++)
			if(fabs(a[t][i])<fabs(a[j][i]))t=j;
			if(fabs(a[t][i])<eps)continue;//把非0元素所在行交换到当前行
			if(r!=t)swap(a[r],a[t]);//第r行第一项变成1
			double tmp=a[r][i];
			for(int j=i;j<=n+1;j++) a[r][j]/=tmp;//变成上三角用第i行去消掉其他所有行的第c列
		for(int j=r+1;j<=n;j++){
		tmp=a[j][i];
		if(fabs(tmp)<eps)continue;
		for(int k=i;k<=n+1;k++) a[j][k]-=a[r][k]*tmp;
		}
		dim++;
	}
	if(dim<n){
    	for(inti=dim+1;i<=n;i++) if(fabs(a[i][n+1])>eps) return -1;//无解
		return dim;//无穷多解
	}//唯一解
	ans[n]=a[n][n+1];
	for(inti=n-1;i>=1;i--){
	ans[i]=a[i][n+1];
	for(intj=i+1;j<=n;j++)
	ans[i]-=a[i][j]*ans[j];
	}
return dim;
}
    

任意模数行列式

//P7112【模板】行列式求值https://www.luogu.com.cn/problem/P7112
ll n,p[N][N],mod;
ll det(){
	assert(mod!=0);
	ll ans=1;
	for(int i=1;i<=n;++i){
		for(int j=i+1;j<=n;++j) 
            while(p[j][i]!=0){//gcdstep辗转相减
	ll t=p[i][i]/p[j][i];
	if(t) for(int k=i;k<=n;++k) p[i][k]=(p[i][k]-p[j][k]*t)%mod;
	swap(p[i],p[j]);
	ans*=-1;
	}
	ans=ans*p[i][i]%mod;
	if(!ans)return 0;
}
return(ans+mod)%mod;
}

抑或方程组

bitset<1010> p[2010];
//p[1~n]：增广矩阵，0位置为常数
//n为未知数个数，m为方程个数，返回方程组的解（多解/无解返回一个空的vector）
vector<bool>GaussElimination(intn,intm){
	//循环消去第i个元
	for(inti=1;i<=n;i++){
	int cur=i;
	while(cur<=m&&!p[cur].test(i))cur++;
	//第i个元的所有系数均为0，有多解
	if(cur>m)returnstd::vector<bool>(0);
	if(cur!=i)swap(p[cur],p[i]);
	for(int j=1;j<=m;j++)
	if(i!=j&&p[j].test(i))p[j]^=p[i];
    }
    vector<bool> ans(n+1,0);
    for(int i=1;i<=n;i++) ans[i]=p[i].test(0);
    return ans;
}

线性基

#include<bits/stdc++.h>
#define ll long long
using namespace std;
structL_B{
	ll d[61],p[61];
	ll cnt;
	L_B(){
		memset(d,0,sizeof(d));
		memset(p,0,sizeof(p));
		cnt=0;
	}
	bool insert(llval){//插入是否插入成功
	for(ll i=60;i>=0;i--)
		if(val&(1LL<<i)){
			if(!d[i]){
			d[i]=val;
			break;
			}
		val^=d[i];
		}
  	return val>0;
	}
	ll query_max(){//取若干个数求异或最大值
		ll ret=0;
		for(ll i=60;i>=0;i--)
			if((ret^d[i])>ret)
			ret^=d[i];
		return ret;
	}
	ll query_min(){//取若干个数求异或最小值
		for(ll i=0;i<=60;i++)
		if(d[i])
		return d[i];
    	return 0;
	}
	void rebuild(){//重构线性基
		for(ll i=60;i>=0;i--)
			for(ll j=i-1;j>=0;j--)
				if(d[i]&(1LL<<j)) d[i]^=d[j];
				for(ll i=0;i<=60;i++)
					if(d[i]) p[cnt++]=d[i];
	}
	ll kthquery(ll k){//查询第k大，之前需要rebuild
		ll ret=0;
		if(k>=(1LL<<cnt))return -1;
		for(ll i=60;i>=0;i--)
		if(k&(1LL<<i))
		ret^=p[i];
	return ret;
	}
}lb;
L_B merge(const L_B& n1,const L_B& n2){//将线性基n2插入线形基n1
L_B ret=n1;
for(ll i=60;i>=0;i--) if(n2.d[i])
ret.insert(n1.d[i]);
return ret;
}

第一类斯特林数（求n 选 i）

#include<bits/stdc++.h>
typedef long long LL;
const int N=550050;
const int mod=167772161;
LL pow_mod(LL a,LL b){
	LL ans=1;
	for(;b;b>>=1,a=a*a%mod)if(b&1)ans=ans*a%mod;
return ans;
}
int L,rev[N];
LL w[N],inv[N],fac[N],ifac[N];
void Init(int n){
	L=1;
	while(L<=n)L<<=1;
	for(int i=1;i<L;++i)
	rev[i]=(rev[i>>1]>>1)|((i&1)*L/2);	
	LL wn=pow_mod(3,(mod-1)/L);
	w[L>>1]=1;
	for(int i=L>>1;i<L;++i)w[i+1]=w[i]*wn%mod;
	for(int i=(L>>1)-1;i;--i)w[i]=w[i<<1];
}
void DFT(LL* A,int len){
	int k=__builtin_ctz(L)-__builtin_ctz(len);
	for(int i=1;i<len;++i){
	int j=rev[i]>>k;
	if(j>i)std::swap(A[i],A[j]);
	}
	for(int h=1;h<len;h<<=1)
		for(int i=0;i<len;i+=(h<<1))
			for(int j=0;j<h;++j){
				LL t=A[i+j+h]*w[j+h]%mod;
				A[i+j+h]=A[i+j]-t;
                A[i+j]+=t;
			}
		for(int i=0;i<len;++i)A[i]%=mod;
}
void IDFT(LL* A,int len){
	std::reverse(A+1,A+len);
	DFT(A,len);
	int v=mod-(mod-1)/len;
	for(int i=0;i<len;++i)A[i]=A[i]*v%mod;
}
void offset(const LL* f,int n,LL c,LL* g){
//g(x)=f(x+c)
//g[i]=1/i!sum_{j=i}^nj!f[j]c^(j-i)/(j-i)!
	static LL tA[N],tB[N];
	int l=1;while(l<=n+n)l<<=1;
	for(int i=0;i<n;++i)tA[n-i-1]=f[i]*fac[i]%mod;
	LL pc=1;
	for(int i=0;i<n;++i,pc=pc*c%mod)tB[i]=pc*ifac[i]%mod;
	for(int i=n;i<l;++i)tA[i]=tB[i]=0;
	DFT(tA,l);DFT(tB,l);
	for(int i=0;i<l;++i)tA[i]=tA[i]*tB[i]%mod;
	IDFT(tA,l);
	for(int i=0;i<n;++i)
	g[i]=tA[n-i-1]*ifac[i]%mod;
}
void Solve(int n,LL* f){
	if(n==0)return void(f[0]=1);
	static LL tA[N],tB[N];
	int m=n/2;
	Solve(m,f);
	int l=1; while(l<=n)l<<=1;
	offset(f,m+1,m,tA);
	for(int i=0;i<=m;++i)tB[i]=f[i];
	for(int i=m+1;i<l;++i)tA[i]=tB[i]=0;
	DFT(tA,l);DFT(tB,l);
	for(int i=0;i<l;++i)tA[i]=tA[i]*tB[i]%mod;
	IDFT(tA,l);
	if(n&1)
	for(int i=0;i<=n;++i)
	f[i]=((i?tA[i-1]:0)+(n-1)*tA[i])%mod;
	else
	for(int i=0;i<=n;++i) f[i]=tA[i];
}
LL f[N];
int main(){
	int n;
	scanf("%d",&n);
	Init(n*2);
    inv[1]=1;
	for(int i=2;i<=n;++i)inv[i]=-(mod/i)*inv[mod%i]%mod;
	fac[0]=ifac[0]=1;
	for(int i=1;i<=n;++i){
		fac[i]=fac[i-1]*i%mod;
		ifac[i]=ifac[i-1]*inv[i]%mod;
	}
	Solve(n,f);
	for(int i=0;i<=n;++i)
		printf("%lld",(f[i]+mod)%mod);
	return 0;
}

第一类斯特林数（i选k）

#include<bits/stdc++.h>
using namespace std;
#define Int register int
#define mod 167772161
#define MAXN 531072
#define Gi 3
long long quick_pow(long long a,long long b ,long long c){
	int res=1;
	while(b){
	if(b&1)res=1ll*res*a%c;
	a=1ll*a*a%c;
	b>>=1;
	}
	return res;
}
int limit=1,l,r[MAXN];
void NTT(int*a,int type){
	for(Int i=0;i<limit;++i)if(i<r[i]) swap(a[i],a[r[i]]);
	for(Int mid=1;mid<limit;mid<<=1){
		int Wn=quick_pow(Gi,(mod-1)/(mid<<1),mod);
		if(type==-1)Wn=quick_pow(Wn,mod-2,mod);
		for(Int R=mid<<1,j=0;j<limit;j+=R){
			for(Int k=0,w=1;k<mid;++k,w=1ll*w*Wn%mod){
				int x=a[j+k],y=1ll*w*a[j+k+mid]%mod;
				a[j+k]=(x+y)%mod,a[j+k+mid]=(x+mod-y)%mod;
			}
        }
	}
	if(type==1)return;
	int Inv=quick_pow(limit,mod-2,mod);
	for(Int i=0;i<limit;++i)a[i]=1ll*a[i]*Inv%mod;
}
int c[MAXN];
void Solve(int len,int* a,int* b){
	if(len==1) return b[0]=quick_pow(a[0],mod-2,mod),void();
	Solve((len+1)>>1,a,b);
	limit=1,l=0;
	while(limit<(len<<1))limit<<=1,l++;
	for(Int i=0;i<limit;++i)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
	for(Int i=0;i<len;++i)c[i]=a[i];
	for(Int i=len;i<limit;++i)c[i]=0;
	NTT(c,1);NTT(b,1);
	for(Int i=0;i<limit;++i)b[i]=1ll*b[i]*(2+mod-1ll*c[i]*b[i]%mod)%mod;
	NTT(b,-1);
	for(Int i=len;i<limit;++i)b[i]=0;
}
void deravitive(int* a,int n){
	for(Int i=1;i<=n;++i)a[i-1]=1ll*a[i]*i%mod;
	a[n]=0;
}
void inter(int* a,int n){
	for(Int i=n;i>=1;--i)a[i]=1ll*a[i-1]*quick_pow(i,mod-2,mod)%mod;
	a[0]=0;
}
int b[MAXN];
void Ln(int* a,int n){
	memset(b,0,sizeof(b));
	Solve(n,a,b);deravitive(a,n);
	while(limit<=n)limit<<=1,l++;
	for(Int i=0;i<limit;++i)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
	NTT(a,1),NTT(b,1);
	for(Int i=0;i<limit;++i)a[i]=1ll*a[i]*b[i]%mod;
	NTT(a,-1);
	inter(a,n);
	for(Int i=n+1;i<limit;++i)a[i]=0;
}
int F0[MAXN];
void Exp(int* a,int* B,int n){
	if(n==1)return B[0]=1,void();
	Exp(a,B,(n+1)>>1);
	for(Int i=0;i<limit;++i)F0[i]=B[i];
	Ln(F0,n);
	F0[0]=(a[0]+1+mod-F0[0])%mod;
    for(Int i=1;i<n;++i)F0[i]=(a[i]+mod-F0[i])%mod;
	NTT(F0,1);NTT(B,1);
	for(Int i=0;i<limit;++i)B[i]=1ll*F0[i]*B[i]%mod;
	NTT(B,-1);
	for(Int i=n;i<limit;++i)B[i]=0;
}
int read(){
	int x=0;
    char c=getchar();
    int f=1;
	while(c<'0'||c>'9'){if(c=='-')f=-f;c=getchar();}
	while(c>='0'&&c<='9'){x=(int)((int)(x<<3)%mod+(int)(x<<1)%mod+c-'0')%mod;c=getchar();}
	return x*f;
}
void write(int x){
	if(x<0){x=-x;putchar('-');}
	if(x>9)write(x/10);
	putchar(x%10+'0');
}
int n,k;
int fac[MAXN],A[MAXN],B[MAXN];
signed main(){
	n=read(),k=read();
	for(Int i=0;i<n;++i)A[i]=quick_pow(i+1,mod-2,mod);
	Ln(A,n);
	for(Int i=0;i<n;++i)A[i]=1ll*A[i]*k%mod;
	Exp(A,B,n);fac[0]=1;
	for(Int i=1;i<=max(n,k);++i)fac[i]=1ll*fac[i-1]*i%mod;
	for(Int i=n;i>=k;--i)B[i]=B[i-k];
	for(Int i=0;i<k;++i)B[i]=0;int Inv=quick_pow(fac[k],mod-2,mod);
	for(Int i=0;i<=n;++i)write(1ll*B[i]*fac[i]%mod*Inv%mod),
	putchar(' ');
	putchar('\n');
return 0;
}
            

第二类斯特林数 n选i

#include<bits/stdc++.h>
using namespace std;
#define int long long
inline void read(int& x){
	char c=getchar();x=0;int f=1;
	while(c>'9'||c<'0'){if(c=='-')f=-1;c=getchar();}
	while(c<='9'&&c>='0')x=(x<<1)+(x<<3)+c-'0',c=getchar();
	x=x*f;
}
const int p=167772161ll,w=3,N=2e6+10;
inline int qpow(int a,int b){
	int k=1ll;
	while(b){
		if(b&1)k=k*a%p;
	a=a*a%p;
	b=b>>1;
	}
return k;
}
int inv[N],n,f[N],g[N],lim,len,rev[N];
inline int upmod(int x){
return (x%p+p)%p;
}
inline void ntt(int* a,int f){
	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){
		int wn=qpow(w,(((p-1)/(mid<<1)*f)+p-1));
		for(int j=0;j<lim;j+=(mid<<1)){
			int g=1;
			for(int k=0;k<mid;k++,g=g*wn%p){
			int x=a[k+j],y=g*a[k+j+mid]%p;
			a[k+j]=upmod(x+y);
			a[k+j+mid]=upmod(x-y+p);
			}
		}
	}
	if(f==-1){
	int Inv=qpow(lim,(p-2));
	for(int i=0;i<lim;i++)a[i]=a[i]*Inv%p;
	}
}
signed main(){
	read(n);n++;
	inv[0]=1;
	for(int i=1;i<n;i++)inv[i]=inv[i-1]*i%p;
	for(int i=1;i<n;i++)inv[i]=qpow(inv[i],p-2);
	for(int i=0;i<n;i++){
	f[i]=(i&1?(p-inv[i]):inv[i]);
	g[i]=qpow(i,n-1)*inv[i]%p;
     }
	lim=1,len=0;
	while(lim<=(n<<1))len++,lim<<=1;
	for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));
	ntt(f,1);ntt(g,1);
	for(int i=0;i<lim;i++)f[i]=f[i]*g[i]%p;
	ntt(f,-1);
	for(int i=0;i<n;i++)printf("%lld",f[i]);
}

第二类斯特林数（i选k）

#include<algorithm>
#include<cstdio>
#define mod 167772161
#define G 3
#define Maxn 270000
using namespace std;
int n,k,r[Maxn<<2];
long long invn,invG;
long long fac[Maxn],inv[Maxn];
long long powM(long long a,long long t=mod-2){
	long long ans=1,buf=a;
	while(t){	if(t&1)ans=(ans*buf)%mod;
	buf=(buf*buf)%mod;
	t>>=1;
	}
   return ans;
}
void NTT(long long* f,bool op,int n){
	for(int i=0;i<n;i++)
		if(r[i]<i)swap(f[r[i]],f[i]);
	for(int len=1;len<n;len<<=1){
		int w=powM(op==1?G:invG,(mod-1)/len/2);
		for(int p=0;p<n;p+=len+len){
			long long buf=1;
			for(int i=p;i<p+len;i++){
				int sav=f[i+len]*buf%mod;
                 f[i+len]=f[i]-sav;
				if(f[i+len]<0)f[i+len]+=mod;
				f[i]=f[i]+sav;
				if(f[i]>=mod)f[i]-=mod;
				buf=buf*w%mod;
			}//F(x)=FL(x^2)+x*FR(x^2)
			//F(W^k)=FL(w^k)+W^k*FR(w^k)
		//F(W^{k+n/2})=FL(w^k)-W^k*FR(w^k)
		}
	}
}
long long g[Maxn<<2];
void rev(long long* f,int len){
	for(int i=0;i<len;i++)g[i]=f[i];
	for(int i=0;i<len;i++)f[len-i-1]=g[i];
}
//令f=f*g(modx^lim)
void times(long long* f,long long* gg,int len,int lim){
	int m=len+len,n;
	for(int i=0;i<len;i++)g[i]=gg[i];
	for(n=1;n<m;n<<=1);invn=powM(n);
	for(int i=len;i<n;i++)g[i]=0;
	for(int i=0;i<n;i++)
	r[i]=(r[i>>1]>>1)|((i&1)?n>>1:0);
	NTT(f,1,n);NTT(g,1,n);
	for(int i=0;i<n;++i)f[i]=(f[i]*g[i])%mod;
	NTT(f,0,n);
	for(int i=0;i<lim;++i)f[i]=(f[i]*invn)%mod;
	for(int i=lim;i<n;++i)f[i]=0;
}
void Init(int lim){
	inv[1]=inv[0]=fac[0]=1;
	for(int i=1;i<=lim;i++)fac[i]=fac[i-1]*i%mod;
	for(int i=2;i<=lim;i++)
	inv[i]=inv[mod%i]*(mod-mod/i)%mod;
	for(int i=2;i<=lim;i++)inv[i]=inv[i-1]*inv[i]%mod;
	for(int i=1;i<=lim;i++)inv[i]=powM(fac[i]);
}
long long p[Maxn<<2];//求出F(x-c)
void fminus(long long* s,long long* f,int len,int c){
	c=mod-c;
	for(int i=0;i<len;i++)
	p[len-i-1]=f[i]*fac[i]%mod;
	long long buf=1;
	for(int i=0;i<len;i++,buf=buf*c%mod)
	s[i]=buf*inv[i]%mod;
	times(p,s,len,len);
	for(int i=0;i<len;i++)s[len-i-1]=p[i]*inv[len-i-1]%mod;
	for(int i=len;i<len+len;i++)s[i]=0;
}
long long f[Maxn<<2],s[Maxn<<2];
void solve(long long* f,int n){
	if(n==1){f[0]=0;f[1]=1;}
	else if(n&1){
		solve(f,n-1);f[n]=0;//再乘上(x-n+1)就好了
		for(int i=n;i>0;i--)
		f[i]=(f[i-1]+(mod-n+1)*f[i])%mod;
		f[0]=f[0]*(mod-n+1)%mod;
	}
	else{
		solve(f,n/2);//S(x)=F(x+n/2)
		fminus(s,f,n/2+1,n/2);
		times(f,s,n/2+1,n+1);
	}
}
void invp(long long* f,int len){
	for(int i=0;i<k+1;i++)s[i]=p[i]=0;//注意清空
	long long* r=s,*rr=p;
	int n=1;for(;n<len;n<<=1);
	rr[0]=powM(f[0]);
	for(int len=2;len<=n;len<<=1){
		for(int i=0;i<len;i++)
		r[i]=rr[i]*2%mod;
		times(rr,rr,len/2,len);
		times(rr,f,len,len);
		for(int i=0;i<len;i++)
		rr[i]=(r[i]-rr[i]+mod)%mod;
	}
    for(int i=0;i<len;i++) f[i]=rr[i];
}
int main(){
	scanf("%d%d",&n,&k);
	if(k>n){
		for(int i=0;i<=n;i++)printf("0");
	return 0;
	}
    invG=powM(G);
	Init(k);solve(f,k+1);
	for(int i=0;i<k+1;i++)f[i]=f[i+1];
	rev(f,k+1);
	for(int i=n-k+1;i<k+1;i++)f[i]=0;
	for(int i=k+1;i<n-k+1;i++)f[i]=0;
	invp(f,n-k+1);
	for(int i=0;i<k;i++)printf("0");
	for(int i=0;i<n-k+1;i++)printf("%lld",f[i]);
    return 0;
}