下面构造生成函数(F(x)=sum_{i=0}^{infty} dp(i)x^i)
根据递推式,可以列出(F(x)=AxF(x)+BF^2(x)+kx),然后移项化简后可得(F(x)=frac{1-Axpmsqrt{(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))求出答案
( iny ext{这方法去年就有洛谷日报了,一年了还不会/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
",(f[r]-f[l-1]+mod)%mod);
}
return 0;
}