已知\(a_{i} = \sum_{j=1}^{i} \{^{i} _{j} \}b_{j}\), 給出\(a_{1} 到 a_{n}\) ;html
求\(b_{l} 到 b_{r}\)在\(1e9+7\)的意義下取模的值;c++
\(1 \le l \le r \le n \le 10^5\)優化
\(r-l \le 100\)spa
\(0 \le a_{i} \lt 10^9 + 7\)code
斯特林反演(http://www.javashuo.com/article/p-qionilex-kr.html):htm
因爲\(r-l\)比較小,能夠直接暴力求,只須要快速求出\([^n_1] - [^n_n]\)blog
\([_m^n] = [x^m] x^{\overline n}\) , 其中\(x^{\overline n} = \Pi_{i=1}^n(x+i-1)\)get
直接作是\(n log ^2n\)的,假設在計算\(F_n(x) = x^{\overline n}\)時已經計算了\(F_\frac{n}{2}(x)\),記\(m = \frac{n}{2} , F_{m} = \sum_{i=0}^{m}a_{i} x^{i}\)it
只須要計算\(F'_{m} (x) = F_{m}(x+m)\)和\(F_{m}\)相乘;class
\[ F'_m(x) = \sum_{i=0}^{m} a_{i} (x+m)^i \\ = \sum_{i=0}^{m} a _{i}\sum_{j=0}^{i}(^i_j)x^j m^{i-j} \\ = \sum_{i=0}^{m} \sum_{j=0}^{i} a_{i} \frac{i!}{j!(i-j)!}x^{j}m^{i-j} \\ = \sum_{i=0}^{m} \frac{x^i}{i!} \sum_{j=0}^{m-i} a_{i+j}(i+j)! \ \frac{m^{j}}{j!} \]
記$A(x)=\sum_{i=0}^{m} a_{i} i! x^{i}, B(x) = \sum_{i=0}^{m}\frac{m^{i}}{i!} x^{m-i} $
相乘以後取第\(m+i\)項能夠算後面的東西
但是不是\(ntt\)模數,直接取模會爆\(long \ long\)須要用到拆係數\(ntt\);
將\(A(x)*B(x)\)分解成\((k * A1(x) + A2(x)) * (k * B1(x) + B2(x)) , k 通常取\sqrt{n} (2^{15} ) 左右\)
暴力是7次\(fft\),能夠優化到須要4次\(fft\)
兩次\(fft\)合併成一次的作法,須要計算\(dft(A(x))和dft(B(x))\)
記: \[\begin{align} P(x) = A(x) + iB(x)\\ Q(x) = A(x) -iB(x) \end{align}\]
考慮直接作\(P\)爲\(P’\),考慮\(Q \ dft\)以後的結果\(Q'\)
\[ \begin{align} Q'[k] &= A(w_{n}^{k}) - iB(w_{n}^{k}) \\ &= \sum_{j=0}^{n-1} (a_{j}-ib_{j})w^{jk}_{n} \\ &= \sum_{j=0}^{n-1}(a_{j}-ib_{j})(cos(\frac{2\pi jk}{n} ) + isin(\frac{2\pi jk}{n}) ) \\ &= \sum_{j=0}^{n-1}(a_{j}cos(\frac{2\pi jk}{n}) + b_{j}sin(\frac{2\pi jk}{n}))- i(b_{j}cos(\frac{2\pi jk}{n}) - a_{j}sin(\frac{2\pi jk}{n})) \\ &= conj \sum_{j=0}^{n-1} (a_{j}cos(\frac{-2\pi jk}{n}) - b_{j}sin(\frac{-2\pi jk}{n} ) ) +i( b_{j}cos(\frac{-2\pi jk}{n}) + a_{j}sin(\frac{-2\pi jk}{n})) \\ &= conj \sum_{j=0}^{n-1} (a_{j}+ib_{j})(cos(\frac{-2\pi jk}{n})+isin(\frac{-2\pi jk}{n})) \\ &= conj \ \sum_{j=0}^{n-1} (a_{j}+ib_{j})w^{-jk}_{n} \\ &= conj \ P'[n-k] \end{align} \]
因此由\(P'\)直接獲得\(Q'\):
即:\[ \begin{align} A'(x) = \frac{P'(x) + Q'(x)}{2} \\ B'(x) = \frac{P'(x) - Q'(x)}{2i}\end{align} g \]
結果也存在實部和虛部裏直接\(idft\)回去就行了;
注意精度;
#include<bits/stdc++.h> #define ld long double #define ll long long using namespace std; const int N=400010,mod=1e9+7; const ld pi=acos(-1); int n,l,r,rev[N],fac[N],inv[N],p[N],a[N],b[N]; char gc(){ static char*p1,*p2,s[1000000]; if(p1==p2)p2=(p1=s)+fread(s,1,1000000,stdin); return(p1==p2)?EOF:*p1++; } int rd(){ int x=0;char c=gc(); while(c<'0'||c>'9')c=gc(); while(c>='0'&&c<='9')x=(x<<1)+(x<<3)+c-'0',c=gc(); return x; } int pw(int x,int y){ int re=1; while(y){ if(y&1)re=(ll)re*x%mod; y>>=1,x=(ll)x*x%mod; } return re; } struct C{ ld x,y; C(ld _x=0,ld _y=0):x(_x),y(_y){}; C operator +(const C&A)const{return C(x+A.x,y+A.y);} C operator -(const C&A)const{return C(x-A.x,y-A.y);} C operator *(const C&A)const{return C(x*A.x-y*A.y,x*A.y+y*A.x);} C operator /(const ld&A)const{return C(x/A,y/A);} C operator !()const{return C(x,-y);} }; void fft(C*A,int len,int f){ for(int i=0;i<len;++i)if(i<rev[i])swap(A[i],A[rev[i]]); for(int i=1;i<len;i<<=1){ C wn=C(cos(pi/i),f*sin(pi/i)); for(int j=0;j<len;j+=i<<1){ C w=C(1,0); for(int k=0;k<i;++k,w=w*wn){ C x=A[j+k],y=w*A[j+k+i]; A[j+k]=x+y,A[j+k+i]=x-y; } } } if(!~f){for(int i=0;i<len;++i){A[i]=A[i]/len;}} } void mul(int*A,int*B,int n){ static C t1[N],t2[N],t3[N],t4[N]; int L,len; for(L=0,len=1;len<(n+1)<<1;++L,len<<=1); for(int i=0;i<len;++i){rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));} for(int i=0;i<=n;++i)t1[i]=C(A[i]>>15,A[i]&0x7fff),t2[i]=C(B[i]>>15,B[i]&0x7fff); for(int i=n+1;i<len;++i)t1[i]=t2[i]=C(0,0); fft(t1,len,1),fft(t2,len,1); for(int i=0;i<len;++i){ C x1=t1[i],y1=!t1[(len-i)&(len-1)]; C x2=t2[i],y2=!t2[(len-i)&(len-1)]; t3[i]=(x1+y1)*x2*C(0.5,0); t4[i]=(x1-y1)*x2*C(0,-0.5); // t3[i]=((x1+y1)*(x2+y2)+(x1+y1)*(x2-y2))*C(0.25,0); // t4[i]=((x1-y1)*(x2+y2)+(x1-y1)*(x2-y2))*C(0,-0.25); } fft(t3,len,-1),fft(t4,len,-1); for(int i=0;i<=n<<1;++i){ A[i]=((((ll)(t3[i].x+0.1)%mod)<<30)+(((ll)(t3[i].y+t4[i].x+0.1)%mod)<<15)+(ll)(t4[i].y+0.1))%mod; } } void solve(int*A,int n){ if(n==0){A[0]=1;return;} if(n==1){A[1]=1;return;} if(n&1){ solve(A,n-1); for(int i=n;i;--i)A[i]=(A[i-1]+(ll)A[i]*(n-1))%mod; A[0]=0; return; } static int t1[N],t2[N]; int m=n>>1; solve(A,m); for(int i=0;i<=m;++i)t1[i]=(ll)A[i]*fac[i]%mod; for(int i=0,t=1;i<=m;++i,t=(ll)t*m%mod)t2[m-i]=(ll)t*inv[i]%mod; mul(t1,t2,m); for(int i=0;i<=m;++i){t1[i]=(ll)inv[i]*t1[m+i]%mod;} mul(A,t1,m); } int main(){ freopen("gift.in","r",stdin); freopen("gift.out","w",stdout); n=rd();l=rd();r=rd(); for(int i=1;i<=n;++i)p[i]=rd(); for(int i=fac[0]=inv[0]=1;i<=n;++i){ fac[i]=(ll)fac[i-1]*i%mod; inv[i]=pw(fac[i],mod-2); } solve(a,l-1); for(int i=l-1;i<=r;++i){ for(int j=1;j<=i;++j){ int t=(ll)a[j]*p[j]%mod; if((i-j)&1)b[i]=(b[i]-t+mod)%mod; else b[i]=(b[i]+t)%mod; } for(int j=i+1;j;--j){a[j]=(a[j-1]+(ll)a[j]*i)%mod;}a[0]=0; // b[i]=(b[i]-b[i-1]+mod)%mod; } for(int i=l;i<=r;++i)printf("%d ",(b[i]-b[i-1]+mod)%mod); return 0; }