题解 贝尔数

传送门

因为crt板子忘了只有暴力分……

发现模数很奇怪，分解后发现等于 \(31*35*41*43*47\)
又发现下面给了个对任意质数成立的递推式
当 \(p=2\) 时发现就是个简单的矩阵加速了
于是想到对上面5个模数都跑一遍，然后crt合并

• 当模数很奇怪时，考虑质因数分解一下
• 当一个式子对任意质数成立时，考虑选几个质数求一下然后crt合并

crt板子：

ll crt(int n, ll* r, ll* m) {
	ll M=1, ans=0;
	for (int i=1; i<=n; ++i) M*=m[i];
	for (int i=1; i<=n; ++i) {
		ll w=M/m[i], x, y;
		exgcd(w, m[i], x, y);
		ans=(ans+r[i]*w*x)%M;
	}
	return (ans+M)%M;
}

Code:

#include <bits/stdc++.h>
using namespace std;
#define INF 0x3f3f3f3f
#define N 100010
#define ll long long
#define reg register int
//#define int long long

char buf[1<<21], *p1=buf, *p2=buf;
#define getchar() (p1==p2&&(p2=(p1=buf)+fread(buf, 1, 1<<21, stdin)), p1==p2?EOF:*p1++)
inline int read() {
	int ans=0, f=1; char c=getchar();
	while (!isdigit(c)) {if (c=='-') f=-f; c=getchar();}
	while (isdigit(c)) {ans=(ans<<3)+(ans<<1)+(c^48); c=getchar();}
	return ans*f;
}

int n;
const ll mod=95041567;
inline void md(ll& a, ll b) {a+=b; a=a>=mod?a-mod:a;}
inline ll md(ll a) {return a>=mod?a-mod:a;}

namespace force{
	ll C[1010][1010], f[1010];
	void init() {
		C[1][0]=C[1][1]=1;
		for (int i=2; i<=1000; ++i) {
			C[i][0]=1;
			for (int j=1; j<=i; ++j) 
				C[i][j]=md(C[i-1][j]+C[i-1][j-1]);
		}
		f[0]=f[1]=1;
		for (int i=2; i<=1000; ++i)
			for (int j=0; j<i; ++j)
				md(f[i], C[i-1][j]*f[j]%mod);
	}
	void solve() {printf("%lld\n", f[n]);}
}

namespace task1{
	const ll div[]={31, 37, 41, 43, 47};
	ll r[10], mod2;
	inline void md2(ll& a, ll b) {a+=b; a=a>=mod2?a-mod2:a;}
	struct matrix{
		int n, m;
		ll a[100][100];
		matrix() {}
		matrix(int x, int y) {n=x; m=y; memset(a, 0, sizeof(a));}
		inline void resize(int x, int y) {n=x; m=y; memset(a, 0, sizeof(a));}
		inline void put() {for (int i=1; i<=n; ++i) {for (int j=1; j<=m; ++j) cout<<a[i][j]<<' '; cout<<endl;}}
		inline ll* operator [] (int t) {return a[t];}
		inline matrix operator * (matrix b) {
			matrix c(n, b.m);
			for (reg i=1; i<=n; ++i)
				for (reg k=1; k<=m; ++k) 
					for (reg j=1; j<=b.m; ++j)
						md2(c[i][j], a[i][k]*b[k][j]%mod2);
			return c;
		}
	}mat, t;
	matrix qpow(matrix a, ll b) {
		matrix ans=a; --b;
		while (b) {
			if (b&1) ans=ans*a;
			a=a*a; b>>=1;
		}
		return ans;
	}
	void exgcd(ll a, ll b, ll& x, ll& y) {
		if (!b) {x=1; y=0; return ;}
		exgcd(b, a%b, y, x);
		y-=a/b*x;
	}
	ll crt() {
		ll M=mod, s=0;
		for (int i=0; i<5; ++i) {
			ll w=M/div[i], x, y;
			exgcd(w, div[i], x, y);
			s=(s+r[i]*w*x)%mod;
		}
		return (s+mod)%mod;
	}
	void solve() {
		for (int k=0; k<5; ++k) {
			mod2=div[k];
			mat.resize(1, 50);
			for (int i=1; i<=50; ++i) mat[1][i]=force::f[i]%mod2;
			// cout<<"mat: "; mat.put();
			t.resize(50, 50);
			for (int i=1; i<50; ++i) t[i+1][i]=1;
			t[50-div[k]+1][50]=t[50-div[k]+2][50]=1;
			// cout<<"---t---"<<endl; t.put();
			// cout<<"div: "<<div[k]<<endl;
			t=qpow(t, n-1);
			mat=mat*t;
			// cout<<"mat2: "; mat.put();
			r[k]=mat[1][1];
		}
		// cout<<"r: "; for (int i=0; i<5; ++i) cout<<r[i]<<' '; cout<<endl;
		printf("%lld\n", crt());
	}
}

signed main()
{
	freopen("bell.in", "r", stdin);
	freopen("bell.out", "w", stdout);

	int T=read();
	force::init();
	while (T--) {
		n=read();
		// force::solve();
		if (n>1) task1::solve();
		else force::solve();
		// force::solve();
	}

	return 0;
}