下面构造生成函数\(F(x)=\sum_{i=0}^{\infty} dp(i)x^i\)
根据递推式,可以列出\(F(x)=AxF(x)+BF^2(x)+kx\),然后移项化简后可得\(F(x)=\frac{1-Ax\pm\sqrt{(1-Ax)^2-4kBx}}{2B}\).由于没有\(dp(0)\),所以\(F(x)\)没有常数项,那么有用的根只有\(F(x)=\frac{1-Ax-\sqrt{(1-Ax)^2-4kBx}}{2B}\)
发现生成函数里只有那个根号是要考虑的,所以考虑求出\(\sqrt{(1-Ax)^2-4kBx}=\sqrt{A^2x^2-(2A+4kB)x+1}\).先令\(a=A^2,b=-(2A+4kB),c=1\),写成\(G(x)=\sqrt{ax^2+bx+c}\),然后设\(H(x)=ax^2+bx+c\),则有\(G(x)=H^{\frac{1}{2}}(x)\),然后对两边求导,得
\(G'(x)=\frac{1}{2}H^{-\frac{1}{2}}(x)H'(x)\)
\(G'(x)H(x)=\frac{1}{2}H^{\frac{1}{2}}(x)H'(x)\)
\(G'(x)H(x)=\frac{1}{2}G(x)H'(x)\)
又因为有\(G(x)=\sum_{i=0}^{\infty}g_ix^i,G'(x)=\sum_{i=0}^{\infty}g_{i+1}(i+1)x^i\),对比\(G'(x)H(x)=\frac{1}{2}G(x)H'(x)\)两边系数后可得\(g_{i+1}=\frac{b(\frac{1}{2}-i)g_i+a(2-i)g_{i-1}}{c(i+1)}\),带入\(a,b,c\)后得\(g_i=\frac{-(2A+4kB)(\frac{3}{2}-i)g_{i-1}+A^2(3-i)g_{i-2}}{i}\),然后可以得到\(G(x)\),进而得到\(F(x)\),最后预处理\(\sum_{i=1}^{n} dp(i)^2\)就可以单次\(O(1)\)求出答案
\(\tiny\text{这方法去年就有洛谷日报了,一年了还不会/kk}\)
#include<bits/stdc++.h>
#define LL long long
using namespace std;
const int N=1e7+10,mod=1e9+7;
int rd()
{
int x=0,w=1;char ch=0;
while(ch<'0'||ch>'9'){if(ch=='-') w=-1;ch=getchar();}
while(ch>='0'&&ch<='9'){x=x*10+(ch^48);ch=getchar();}
return x*w;
}
int fpow(int a,int b){int an=1;while(b){if(b&1) an=1ll*an*a%mod;a=1ll*a*a%mod,b>>=1;}return an;}
int ginv(int a){return fpow(a,mod-2);}
int n,k,a,b,inv[N],f[N];
int main()
{
//////////////
inv[0]=inv[1]=1;
for(int i=2;i<=N-5;++i) inv[i]=(mod-1ll*mod/i*inv[mod%i]%mod)%mod;
n=rd(),k=rd(),a=rd(),b=rd();
f[0]=1,f[1]=(1ll*mod+mod-a-2ll*k*b%mod)%mod;
int yyb=2ll*f[1]%mod,inv2=ginv(2),i2t3=3ll*inv2%mod;
for(int i=2;i<=n;++i) f[i]=1ll*(1ll*yyb*(i2t3-i+mod)%mod*f[i-1]+1ll*a*a%mod*(3-i+mod)%mod*f[i-2])%mod*inv[i]%mod;
f[0]-=1,f[1]=(f[1]+a)%mod;
int invb=ginv((b+b)%mod);
for(int i=0;i<=n;++i) f[i]=1ll*(mod-f[i])*invb%mod;
for(int i=1;i<=n;++i) f[i]=(1ll*f[i]*f[i]+f[i-1])%mod;
int q=rd();
while(q--)
{
int l=rd(),r=rd();
printf("%d\n",(f[r]-f[l-1]+mod)%mod);
}
return 0;
}原文地址：https://www.cnblogs.com/smyjr/p/12264611.html