這題首先我們可以推出遞推式：

我們先把加數個數大於等於 \(2\) 的限制去掉，最後再減回去即可。

\[f_0=1 \]

\[f_n=\sum\limits_{i=1}^{n} j \cdot f_{i-j} \]

暴力程式碼：

#include<bits/stdc++.h>
using namespace std;
template<typename T>inline void read(T& FF){
	FF=0;T RR=1;char CH=getchar();
	for(;!isdigit(CH);CH=getchar())if(CH=='-')RR=-1;
	for(;isdigit(CH);CH=getchar())FF=(FF<<1)+(FF<<3)+(CH^48);
	FF*=RR;
}
int f[100010];
int main(){int T;read(T);while(T--){
	int n;read(n);
	f[0]=1;
	for(int i=1;i<=n;i++){
		f[i]=0;
		for(int j=1;j<=i;j++)
			f[i]+=j*f[i-j];
	}
	cout<<f[n]-n<<endl;}
	return 0;
}
 

然後我們看這個東西怎麼矩乘。

我們來做差。

\[\begin{aligned}f_n-f_{n-1}&=\sum\limits_{i=1}^{n} j \cdot f_{i-j}-\sum\limits_{i=1}^{n-1} j \cdot f_{i-j-1}\\&=f_{i-1}+2\cdot f_{i-2}+\cdots + i\cdot f_{0}-[f_{i-2}+2\cdot f_{i-3}+\cdots + (i-1)\cdot f_{0}]\\&=f_{i-1}+f_{i-2}+\cdots+f_{0}\end{aligned} \]

我們把這個東西叫做 \(g_{i-1}\)

那麼，現在矩陣乘法要做的事情就是利用 \(f_{i-1}\) 和 \(g_{i-1}\) 求出 \(f_i\)，我們發現 \(f_i=f_{i-1}+g_{i-1}\)，故可得矩陣：

\[\begin{bmatrix}1&1\\1&2\\\end{bmatrix}\begin{bmatrix}f_{i-1}\\g_{i-1}\end{bmatrix}=\begin{bmatrix}f_{i}\\g_{i}\end{bmatrix} \]

然後我們可以得到：

\[\begin{bmatrix}0\\1\end{bmatrix}\begin{bmatrix}1&1\\1&2\end{bmatrix}^n \]

用矩乘算一下即可。

\(code\)：

#include <bits/stdc++.h>
#define int long long
using namespace std;
template<typename T>inline void read(T &FF){
	T RR=1;FF=0;char CH=getchar();
	for(;!isdigit(CH);CH=getchar())if(CH=='-')RR=-1;
	for(;isdigit(CH);CH=getchar())FF=(FF<<1)+(FF<<3)+(CH^48);
	FF*=RR;
}
const int N=5,MOD=1000000007;
struct Matrix{
	int a[N][N],n,m;
	Matrix(){
		memset(a,0,sizeof(a));
	}
	void build(){
		for(int i=1;i<=2;i++)a[i][i]=1;
	}
}a,ans;
Matrix operator*(const Matrix x,const Matrix y){
	Matrix ans;
	ans.n=x.m;
	ans.m=y.n;
	for(int k=1;k<=x.m;k++)
		for(int i=1;i<=x.n;i++)
			for(int j=1;j<=y.m;j++)
				ans.a[i][j]=(ans.a[i][j]+1ll*x.a[i][k]*y.a[k][j]%MOD)%MOD;
	return ans;
}
Matrix pw(Matrix a,int k){
	Matrix ans;ans.n=ans.m=2;ans.build();
	for(;k;k>>=1,a=a*a)
		if(k&1)ans=ans*a;
	return ans;
}
signed main(){
	int T;read(T);
	while(T--){
		int n;read(n);
		Matrix a,b;
		a.n=2;a.m=2;
		a.a[1][1]=a.a[1][2]=a.a[2][1]=1;a.a[2][2]=2;
		b.n=2;b.m=1;
		b.a[1][1]=0;b.a[2][1]=1;
		a=pw(a,n);
		a=a*b;
		cout<<(a.a[1][1]-n+MOD)%MOD<<endl;
	}
	return 0;
}