算法学习笔记:生成函数浅讲

生成函数浅讲

感觉这是一个非常牛逼的东西,写了点自己的感悟,可能讲得不是很清楚。

生成函数的定义就比较牛,将数列 \(\{a_i\}\) 写成一个函数 \(A(x)=\sum{a_ix^i}\) 的形式叫做普通生成函数。

此处的 \(x^i\) 没有实际意义,只是一个占位符。对于生成函数来说,绝大数多项式的运算法则都是可以用的。

比如:

\(A(x)+B(x) = \sum\limits{(a_i+b_i)x^i}\)

\(A(x)B(x)=\sum\limits_i\sum\limits_j(a_ib_j)x^{i+j}\)

一些常见的数列都可以写作生成函数的形式,而且往往都有较为简单的封闭形式:

\(\{1,0,0\dots\}\rightarrow 1\)

\(\{0,1,0\dots\}\rightarrow x\)

\(\{1,1,1\dots\}\rightarrow 1+x+x^2+\dots=\frac{1}{1-x}\)

\(\{1,0,1,0,1\dots\}\rightarrow 1+x^2+x^4\dots=\frac{1}{1-x^2}\)

我们用 \([x^n]\) 表示 \(x^n\) 的系数,即数列的第 \(n\) 位。

拓展一下它们的乘积有什么含义

\(x^k\) 相当于将数列向后平移 \(k\) 位。

\(\frac{1}{1-x}\) 相当于对数列做一次前缀和,而 \(\frac1{(1-x)^k}\) 就是求 \(k\) 次前缀和,把这个过程写出来就会发现很像一个旋转了一下的杨辉三角,第 \(i\) 位系数就是 \({i+k-1}\choose{k-1}\)

\(\frac1{1-x^p}\) 后第 \(i\) 位就是原数列第 \(i\) 位往前第 \(p\) 的整数倍的位的前缀和(听着很绕其实很简单,笔者试图用简洁的语言的描述但语文水平不允许),类似的 \(\frac1{(1-x^p)^k}\) 的数列只在 \(p\) 的整数倍的位上有值并且第 \(ip\) 位为 \({i+k-1}\choose{k-1}\)

然后我们就可以用生成函数做一些很牛的事情。

一道基础题:

P10780 食物

题意:在许多不同种类的食物中选出 \(n\) 个,每种食物的限制如下:

  • 承德汉堡:偶数个

  • 可乐:\(0\) 个或 \(1\)

  • 鸡腿:\(0\) 个,\(1\) 个或 \(2\)

  • 蜜桃多:奇数个

  • 鸡块:\(4\) 的倍数个

  • 包子:\(0\) 个,\(1\) 个,\(2\) 个或 \(3\)

  • 土豆片炒肉:不超过一个。

  • 面包:\(3\) 的倍数个

每种食物都是以「个」为单位,只要总数加起来是 \(n\) 就算一种方案。对于给出的 \(n\) ,你需要计算出方案数,对 \(10007\) 取模。(\(n\le 10^{500}\)

\(sol\):首先对于每种食物构造一个数列,数列的第 \(i\) 位表示总数为 \(i\) 时这种食物有多少选法,显然只有 \(0/1\) ,如下:

\(\{1,0,1,0,1\dots\}\rightarrow 1+x^2+x^4\dots=\frac1{1-x^2}\)

\(\{1,1,0,0,0\dots\}\rightarrow 1+x\)

\(\{1,1,1,0,0\dots\}\rightarrow 1+x+x^2\)

\(\{0,1,0,1,0\dots\}\rightarrow x+x^3+x^5\dots=\frac{x}{1-x^2}\)

\(\{1,0,0,0,1,0\dots\}\rightarrow 1+x^4+x^8\dots=\frac1{1-x^4}\)

\(\{1,1,1,1,0\dots\}\rightarrow 1+x+x^2+x^3\)

\(\{1,1,0,0,0\dots\}\rightarrow 1+x\)

\(\{1,0,0,1,0\dots\}\rightarrow 1+x^3+x^6\dots=\frac1{1-x^3}\)

然后把它们乘起来就把数量上的加法转化为了系数上的相加,所以总数为 \(n\) 的答案就变成了 \(x^n\) 的系数。结果化简一下就是 \(\frac x{(1-x)^4}\) ,然后结合上面总结的 \(\frac x{(1-x)^k}\) 的规律,可以得到答案就是 \({n+2\choose3}\),边读入边取模即可。

进阶应用

这道题中的生成函数是比较显然的,但是很多题目需要我们通过生成函数转化一些其他做法,并且最终生成函数的通项的求法也会有很多种,下面看一些不同解法的题。

[国家集训队] 整数的lqp拆分 (特征方程)

定义 \(F_0=0,F_1=1,F_n=F_{n-1}+F_{n-2} (n>1)\) (其实就是斐波那契数列)

\(\sum\prod\limits_{i=1}^m F_{a_i}\)

\(m>0\)

\(a_1,a_2...a_m>0\)
\(a_1+a_2+...+a_m=n\)
由于答案可能非常大,所以要对 \(10^9 + 7\) 取模。(\(1\le n \le 10^{10000}\)

\(sol\) :首先考虑DP,设 \(f_{i,j}\) 表示固定 \(m\)\(i\)\(n\)\(j\) 时的答案。然后就会发现一个非常显然的暴力转移:

\[f_{i,j}=\sum\limits_{k=0}^jf_{i-1,k}\times F_{j-k} \]

按理来说这里的 \(k\) 只能从 \(1\) 枚举到 \(j-1\) 但是 \(F_0=0,f_{i,0}=0\) 所以对答案是没有影响的。最终答案就应该是 \(\sum\limits_{i=1}^nf_{i,n}\)

我们把 \(f_{i-1}\)\(f_i\) 的递推看成一次卷积操作,就会发现 \(f_i=F^i\) ,那么答案就是:

\[[x^n]\sum\limits_{i=1}^nF^i \]

看着似乎不太可做,但是我们发现当 \(i>n\) 时,\(F^i\) 的第 \(n\) 项显然为 \(0\) 所以可以改写成:

\[\begin{aligned} [x^n]\sum\limits_{i=1}^{\infty} F^i\\ =[x^n]\frac F{1-F} \end{aligned} \]

把斐波那契数列的生成函数 \(F=\frac{x}{1-x-x^2}\) 代进去就是 \(\frac x{1-2x-x^2}\) ,这个生成函数可以用特征方程来解出来,这里详细讲一下。

分子就是平移一位是好处理的,而分母转换一下就是:

\[f_0=1,f_1=2,f_{n+1}=2f_n+f_{n-1} \]

我们考虑把后面的式子写成一种等比数列的形式:

\[f_{n+1}-pf_n=q(f_n-pf_{n-1}) \]

然后我们可以求出这个等比数列的通项:

\[\begin{aligned} f_{n+1}-pf_n=q^n(f_1-pf_0)\\ f_{n+1}-pf_n=q^n(2-p) \end{aligned} \]

\(p\)\(q\) 应该满足条件:

\(pq=-1,p+q=2\)

所以 \(p\)\(q\) 是方程 \(x^2-2x-1=0\) 的两个根 \(1+\sqrt2\)\(1-\sqrt2\)

同时这里也可以看出 \(p\)\(q\) 是等价的(可以交换),再结合上面的式子可以得到一个方程组:

\[\begin{aligned} f_{n+1}-qf_n=p^n(2-q)\\ f_{n+1}-pf_n=q^n(2-p) \end{aligned} \]

然后我们把 \(f_n\) 解出来,再把 \(p\)\(q\) 代进去:

\[f_n=\frac{(1+\sqrt2)^{n+1}-(1-\sqrt2)^{n+1}}{2\sqrt2} \]

因为这是平移过的,所以最后答案是

\[f_{n-1}=\frac{(1+\sqrt2)^n-(1-\sqrt2)^n}{2\sqrt2} \]

Two Snuke (拉格朗日插值)

题意:对于一个十元组 \((a_1,a_2,b_1,b_2,c_1,c_2,d_1,d_2,e_1,e_2)\),定义它合法当且仅当满足下列条件:

  • \(0\le a_1<a_2\)

  • \(0\le b_1<b_2\)

  • \(0\le c_1<c_2\)

  • \(0\le d_1<d_2\)

  • \(0\le e_1<e_2\)

  • \(a_1+a_2+b_1+b_2+c_1+c_2+d_1,d_2+e_1+e_2\le n\)

    现在有 \(T\) 组数据,每次对于给定的 \(n\) 求所有合法十元组 \((a_2-a_1)(b_2-b_1)(c_2-c_1)(d_2-d_1)(e_2-e_1)\) 之和对 \(10^9+7\) 取模的结果。

    \(1\le T\le 100,1\le n\le 10^9\)

这道题似乎有比较简单的组合意义解法,但是我觉得用生成函数做这道题更帅一些,其实是不会组合意义做法。之前写过题解,这里就不讲了。

code

代码
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod=1e9+7;
ll fac[100]={1},sum[100],n=16;
inline ll rd()
{
	char c;ll f=1;
	while(!isdigit(c=getchar()))if(c=='-')f=-1;
	ll x=c-'0';
	while(isdigit(c=getchar()))x=x*10+(c^48);
	return x*f;
}
inline ll qp(ll x,ll y)
{
	ll res=1;
	while(y)
	{
		if(y&1) (res*=x)%=mod;
		(x*=x)%=mod,y>>=1;
	}
	return res;
}
ll C(ll n,ll m)
{
	ll s=1;
	for(int i=1;i<=m;i++) (s*=(n-i+1))%=mod;
	return s*qp(fac[m],mod-2)%mod;
}
ll ans(ll k)
{
	if(k<5) return 0;k-=5;
	n=min(16ll,k/2);
	for(int j=0;j<=n;j++)
		sum[j]=C(j+4,4)*C(k-2*j+10,10)%mod;
	k/=2;ll s=0;
	for(int j=1;j<=n;j++) (sum[j]+=sum[j-1])%=mod;
	if(k<=n) return (sum[k]+mod)%mod;
	for(int i=0;i<=n;i++)
	{
		ll s1=sum[i],s2=1;
		for(int j=0;j<=n;j++) if(j!=i)
			s1=s1*(k-j)%mod,s2=s2*(i-j)%mod;
		(s+=s1*qp(s2,mod-2)%mod)%=mod;
	}
	return (s+mod)%mod;
}
int main()
{
	for(int i=1;i<=n;i++) fac[i]=fac[i-1]*i%mod;
	for(ll t=rd();t--;) printf("%lld\n",ans(rd()));
	return 0;
}

Card Deck Score (待定系数)

题意:有一些卡牌,每张卡牌上有一个数字,具体的,有 \(b_i\) 张卡牌上的数字为 \(a_i\)

求出拿走其中 \(m\) 张卡牌的贡献之和。贡献为这些卡牌的乘积。对于本质相同的卡牌组合,只算一次。

  • \(n\leq 16,m\leq 10^{18},b_i\leq 10^{17},1\leq a_i<mod\)

\(sol\)

首先对于每种卡牌构造生成函数:

\[\sum\limits_{j=0}^{b_i} a_i^jx^j=\frac{1-a_i^{b_i+1}x^{b_i+1}}{1-a_ix} \]

最终答案的就是:

\[\begin{aligned} [x^m]\prod\limits_{i=1}^n \frac{1-a_i^{b_i+1}x^{b_i+1}}{1-a_ix}\\ =[x^m]\frac{\prod\limits_{i=1}^n1-a_i^{b_i+1}x^{b_i+1}}{\prod\limits_{i=1}^n1-a_ix} \end{aligned} \]

这里因为 \(n\) 很小,只有 \(16\) 所以可以\(2^n\)暴力把分子部分展开。假设展开有 \(w\) 项,第 \(i\) 项为 \(c_ix^{d_i}\) ,则累加每一项的贡献得到的答案就是:

\[\begin{aligned} \sum\limits_{i=1}^wc_i\times[x^{n-d_i}]\frac{1}{\prod\limits_{i=1}^n1-a_ix} \end{aligned} \]

现在问题变成了 \(2^n\) 次询问分母某一位的值。处理分母可以考虑待定系数,如果可以写成求和的形式就好做多了,于是我们令:

\[\begin{aligned} \prod\limits_{i=1}^n\frac1{1-a_ix}=\sum\limits_{i=1}^n \frac{p_i}{1-a_ix}\\ \sum\limits_{i=1}^n \frac{p_i}{1-a_ix}\prod\limits_{j=1}^n 1-a_jx=1\\ \sum\limits_{i=1}^n p_i\prod\limits_{j=1,j\ne i}^n 1-a_jx=1\\ \end{aligned} \]

一般来说生成函数中的 \(x\) 是没有实际意义的,但是这里我们可以考虑把它当作一个普通的函数,而普通函数有一个优势,可以代值求待定系数。所以我们代入 \(x=\frac1{a_k}\) 那么当 \(j=k\) 时,式子中 \(1-a_jx=0\) ,也就是说只有 \(i=k\) 这一项被保留了下来:

\[\begin{aligned} p_i\prod\limits_{j=1,j\ne i}^n 1-\frac{a_j}{a_i}=1\\ p_i=\frac1{\prod\limits_{j=1,j\ne i}^n 1-\frac{a_j}{a_i}} \end{aligned} \]

于是我们可以求出每一项的待定系数,现在再看一下分母的形式就非常好求了:

\[\begin{aligned} [x^k]\sum\limits_{i=1}^n \frac{p_i}{1-a_ix}\\ \sum\limits_{i=1}^n[x^k]\frac{p_i}{1-a_ix}\\ =\sum\limits_{i=1}^np_i\times a_i^k \end{aligned} \]

至此,我们得到了一个单次 \(O(n)\) 查询某一位的做法,最终复杂度:\(O(n2^n)\)​。

code

代码
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod=998244353,N=1e6+5;
inline ll rd()
{
	char c;ll f=1;
	while(!isdigit(c=getchar()))if(c=='-')f=-1;
	ll x=c-'0';
	while(isdigit(c=getchar()))x=x*10+(c^48);
	return x*f;
}
ll qp(ll x,ll y)
{
    ll res=1;
    while(y)
    {
        if(y&1)res=res*x%mod;
        x=x*x%mod;
        y>>=1;
    }
    return res;
}
ll n,m,a[20],in,b[20],x[20],p[20],nx[N],P[N],w=1;
int main()
{
	n=rd(),m=rd();
    for(int i=1;i<=n;i++)
		a[i]=rd(),b[i]=rd(),x[i]=-qp(a[i],b[i]+1);
    nx[w]=1,P[w]=0;
    for(int i=1;i<=n;i++)
    {
    	int t=w;
        for(int j=1;j<=t;j++)
        {
            nx[++w]=nx[j]*x[i]%mod;
            P[w]=P[j]+b[i]+1;
        }
	}
    for(int i=1;i<=n;i++)
    {
        p[i]=1,in=qp(a[i],mod-2);
        for(int j=1;j<=n;j++) if(i!=j)
            p[i]=p[i]*qp(1-in*a[j]%mod+mod,mod-2)%mod;
    }
    ll ans=0;
    for(int i=1;i<=w;i++)
    {
        if(P[i]>m) continue;
        ll now=m-P[i];
        ll res=0;
        for(int j=1;j<=n;j++)
            (res+=(p[j]*qp(a[j],now)%mod))%=mod;
        ans=(ans+nx[i]*res%mod+mod)%mod;
    }
    cout<<ans;
    return 0;
}

一种暴力做法

上面的题用了一些比较巧妙的方法将生成函数求了出来,但总有一些毒瘤题,很难找到相适应的方法,这时候通常可以根据题意列出方程,但是最终得到的很有可能是一个极其丑陋的二次方程,解出的形式看上去就很难做,我们就有一种方法将它转化为更简洁的形式。

首先,我们有方程 \(G^2+F_1G+F_2=0\) ,其中 \(G,F_1,F_2\) 均为连续可导函数 让我们来证明这个式子: \((4F_2+F_1^2)G'+(F_1F_1'-2F_2')G+2F_1'F_2-F_1F_2'=0\)​。

先对原式求导得:

\[\begin{aligned} 2GG'+F_1G'+G'F_1+F_2'=0\\ G'=-\frac{F_1G'+F_2'}{2G+F_1}\\ G'=-(\frac{F_1'}2+\frac{F_2'-\frac{F_1F_2'}2}{2G+F_1}) \end{aligned} \]

接下来,我们尝试把分母上的 \(2G+F_1\) 变到分子上(相当于消掉二次项),于是将原式变形:

\[\begin{aligned} \frac{G^2+F_1G+F_2}{2G+F_1}=0\\ \frac G2+\frac{F_1}4+\frac{F_2'-\frac{F_1F_1'}2}{2G+F_1}=0\\ \frac 1{2G+F_1}=-\frac{\frac G2+\frac{F_1}4}{F_2'-\frac{F_1F_1'}2}\\ =-\frac{2G+F_1}{4F_2'-2F_1F_1} \end{aligned} \]

再将 \(\frac 1{2G+F_1}\) 代回上面的式子得:

\[\begin{aligned} G'=-\frac{F_1'}2+\frac{(F_2'-\frac{F_1F_2'}2 )(2G+F_1)}{4F_2'-2F_1F_1}\\ (4F_2'-2F_1F_1)G'=-2F_2F_1'+\frac{F_1F_2'}2-2GF_2'-GF_1F_1'+F_1F_2'-\frac{F_1F_2'}2\\ (4F_2+F_1^2)G'+(F_1F_1'-2F_2')G+2F_1'F_2-F_1F_2'=0 \end{aligned} \]

现在,让我们推广一下,当 \(F_1G^2+F_2G+F_3=0\) 时,其中 \(G,F_1,F_2,F_3\) 均为连续可导函数,是否会有类似的关系。

我们将两边同时除以 \(F_1\) 得到

\[\begin{aligned} G^2+\frac{F_2}{F_1}G+\frac{F_3}{F_1}=0 \end{aligned} \]

由刚证明的定理得:

\[(\frac{4F_3}{F_1}+\frac{F_2^2}{F_1^2})G'+(\frac{F_2}{F_1}\times\frac{F_1F_2'-F_1'F_2}{F_1^2}-2\times\frac{F_1F_3'-F_1'F_3}{F_1^2})G+2\frac{F_1F_2'-F_1'F_2}{F_1^2}\times \frac{F_3}{F_1}-\frac{F_2}{F_1}\frac{F_1F_3'-F_1'F_3}{F_1^2}=0 \]

化简一下可以得到

\[(4F_3F_1^2+F_2^2F_1)G'+(F_2F_1F_2'-F_2^2F_1-2F_1^2F_3'+2F_1F_1'F_3)G+2F_1F_2'F_3-2F_1'F_2F_3-F_2F_1F_3+F_1'F_2F_3=0 \]

最终我们得到了这个优美的式子。

下面用一道例题讲讲这个看似抽象的定理具体的应用:

直接自然溢出啥事没有 加强版

题面:给定一个正整数 \(n\),问有多少个长度为 \(n\) 的字符串,满足这个字符串是一个「程序片段」。

具体定义如下:

单个分号 ; 是一个「语句」。

空串 是一个「程序片段」。

如果字符串 A 是「程序片段」,字符串 B 是「语句」,则 AB 是「程序片段」。

如果字符串 A 是「程序片段」,则 {A} 是「语句块」。

如果字符串 A 是「语句块」,则 A 是「语句」,[]A[]()A 都是「函数」。

如果字符串 A 是「函数」,则 (A) 是「函数」,AA() 都是「值」。

如果字符串 A 是「值」,则 (A) 是「值」,A; 是「语句」。

注意:AB 并不代表 BA

\(n\le 10^7\)

sol:

这个题各种字符串的关系极其复杂,于是我们用生成函数来表示,\(A,B,C,D,E\) 分别表示语句,程序片段,语句块,函数,值。

\[\begin{aligned} \begin{equation} \begin{cases} B=AB+1\\ C=x^2B\\ D=x^2C+x^4C+x^2D\\ E=x^2E+D\\ A=x+C+xE \end{cases} \end{equation} \end{aligned} \]

其中 \(B\) 就是我们需要的生成函数,我们试着解一下:

将2式代入3式得:

\[\begin{aligned} D=x^4B+x^6B+x^2D\\ D=\frac{x^4+x^6}{1-x^2}B \end{aligned} \]

再代入4式得:

\[\begin{aligned} E=x^2E+\frac{x^4+x^6}{1-x^2}B\\ E=\frac{x^4+x^6}{1-2x^2+x^4}B \end{aligned} \]

然后代入5式得:

\[\begin{aligned} A=x+x^2B+\frac{x^5+x^7}{1-2x^2+x^4}B \end{aligned} \]

最后代入1式得:

\[\begin{aligned} B=xB+x^2B^2+\frac{x^5+x^7}{1-2x^2+x^4}B^2+1\\ (x-1)B+\frac{x^7+x^6+x^5-2x^4+x^2}{1-2x^2+x^4}B^2+1=0\\ (x^7+x^6+x^5-2x^4+x^2)B^2+(x^5-x^4-2x^3+x-1)B+x^4-2x^2+1=0 \end{aligned} \]

然后我们发现这就是上面讲的那个定理的形式,所以我们接下来就可以大力推式子了!!!

这里我们令:

\[\begin{aligned} F_1=x^7+x^6+x^5-2x^4+x^2\\ F_2=x^5-x^4-2x^3+x-1\\ F_3=x^4-2x^2+1 \end{aligned} \]

代入进去经过化简得到

\[\begin{aligned} (4x^{18}+7x^{17}+5x^{16}-20x^{15}-33x^{14}+5x^{13}+55x^{12}+42x^{11}-67x^{10}-63x^9+59x^8+40x^7-31x^6-13x^5+9x^4+2x^3-x^2)B'+(6x^{17}+8x^{16}-4x^{15}-30x^{14}-63x^{13}+29x^{12}+111x^{11}+27x^{10}-94x^9-78x^8+70x^7+64x^6-39x^5-23x^4+15x^3+3x^2-2x)B+(-x^{15}+3x^{14}+11x^{13}-15x^{12}-30x^{11}+22x^{10}+40x^9-6x^8-37x^7-9x^6+27x^5+5x^4-12x^3+2x)=0\\ \end{aligned} \]

约掉一个 \(x\) 之后可以得到:

\[(4x^{17}+7x^{16}+5x^{15}-20x^{14}-33x^{13}+5x^{12}+55x^{11}+42x^{10}-67x^9-63x^8+59x^7+40x^6-31x^5-13x^4+9x^3+2x^2-x)B'+(6x^{16}+8x^{15}-4x^{14}-30x^{13}-63x^{13}+29x^{11}+111x^{10}+27x^9-94x^8-78x^7+70x^6+64x^5-39x^4-23x^3+15x^2+3x-2x)B+(-x^{14}+3x^{13}+11x^{12}-15x^{11}-30x^{10}+22x^9+40x^8-6x^7-37x^6-9x^5+27x^4+5x^3-12x^2+2x)=0 \]

代入第 \(n\) 项,,令答案为 \(F_n\) 我们就得到了 \(O(n)\) 的递推了。

\[\begin{aligned} (n+2)F_n+2(n-1)F_{n-1}+9(n-2)F_{n-2}-13(n-3)F_{n-3}-31(n-4)F_{n-4}+40(n-5)F_{n-5}+59(n-6)F_{n-6}-63(n-7)F_{n-}-67(n-8)F_{n-8}+42(n-9)F_{n-9}+55(n-10)F_{n-10}+5(n-11)F_{n-11}-33(n-12)F_{n-12}-20(n-13)F_{n-13}+5(n-14)F_{n-14}+7(n-15)F_{n-15}+4(n-16)F_{n-16}+3F_{n-1}+15F_{n-2}-23F_{n-3}-39F_{n-4}+64F_{n-5}+70F_{n-6}-78F_{n-7}-94F_{n-8}+27F_{n-9}+111F_{n-10}+29F_{n-11}-63F_{n-12}-30F_{n-13}-4F_{n-14}+8F_{n-15}+6F_{n-16}-[n=14]+3[n=13]+11[n=12]-15[n=11]-30[n=10]+22[n=9]+40[n=8]-6[n=7]-37[n=6]-9[n=5]+27[n=4]+5[n=3]-12[n=2]=0 \end{aligned} \]

code

代码
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mod=998244353,N=1e7+5;
ll f[N],ans,sum,s1[N],s2,n,h[15];
ll qp(ll x,ll y)
{
	ll res=1;
	while(y)
	{
		if(y&1) (res*=x)%=mod;
		y>>=1,(x*=x)%=mod;
	}
	return res;
}
ll F(int i){return(i>=0)?f[i]:0;}
ll H(int i){return(i<=14)?h[i]:0;}
int main()
{
	cin>>n; sum=s1[0]=s2=f[0]=1;
	for(int i=1;i<=n+2;i++) sum=sum*i%mod,s1[i]=s1[i-1]*i%mod; 
	sum=qp(sum,mod-2);
	for(int i=n+2;i>=1;i--) s1[i]=sum*s1[i-1]%mod*s2%mod,s2=s2*i%mod;
	h[14]=-1,h[13]=3,h[12]=11,h[11]=-15,h[10]=-30,h[9]=22,h[8]=40,h[7]=-6,h[6]=-37,h[5]=-9,h[4]=27,h[3]=5,h[2]=-12;
	for(int i=1;i<=n;i++) f[i]=((2*(i-1)*F(i-1)+9*(i-2)*F(i-2)-13*(i-3)*F(i-3)-31*(i-4)*F(i-4)+40*(i-5)*F(i-5)+59*(i-6)*F(i-6)-63*(i-7)*F(i-7)-67*(i-8)*F(i-8)+42*(i-9)*F(i-9)+55*(i-10)*F(i-10)+5*(i-11)*F(i-11)-33*(i-12)*F(i-12)-20*(i-13)*F(i-13)+5*(i-14)*F(i-14)+7*(i-15)*F(i-15)+4*(i-16)*F(i-16)+3*F(i-1)+15*F(i-2)-23*F(i-3)-39*F(i-4)+64*F(i-5)+70*F(i-6)-78*F(i-7)-94*F(i-8)+27*F(i-9)+111*F(i-10)+29*F(i-11)-63*F(i-12)-30*F(i-13)-4*F(i-14)+8*F(i-15)+6*F(i-16)+H(i))%mod+mod)*s1[i+2]%mod;
	cout<<f[n];
	return 0;
}
posted @ 2025-01-16 20:33  Re_Star  阅读(144)  评论(4)    收藏  举报