FZU 2314

FZU 第十六届程序设计竞赛_重现赛 & FOJ Problem 2314 宝宝会求导

已知 (displaystyle f(x)={1over e^{-x}+1})

求 (x=0) 处泰勒展开第 (k) 项系数 (mod (10^9+7))

---

法一

已知该函数在 (x=0) 处泰勒展开式第 (k) 项系数为 (displaystyle {f^{(k)}(0)over k!})

由于 (({1over k!})mod (10^9+7)) 可以 (O(n)) 求出，故直接考虑求解 (f^{(k)}(0)mod (10^9+7))

不妨设 (displaystyle y=f(x)={1over e^{-x}+1}Rightarrow { ext dyover ext dx}=-{{ ext dover ext dx}(e^{-x}+1)over (e^{-x}+1)^2}={(e^{-x}+1)-1over (e^{-x}+1)^2}=y-y^2)

再次求导可得 (displaystyle { ext d^2 yover ext dx^2}={ ext dover ext dx}({ ext dyover ext dx})={ ext dyover ext dx}cdot { ext dover ext dy}({ ext dyover ext dx})={ ext dyover ext dx}cdot { ext dover ext dy}(y-y^2)=(1-2y)cdot { ext dyover ext dx}=(1-2y)(y-y^2)=y-3y^2+2y^3)

进而可观察出规律，若设 (displaystyle { ext d^n yover ext dx^n}=sum_{i=1}^{n+1}a_iy^i)

则 (displaystyle { ext d^{n+1} yover ext dx^{n+1}}={ ext dover ext dx}({ ext d^n yover ext dx^n})={ ext dyover ext dx}cdot ({ ext dover ext dy}sum_{i=1}^{n+1}a_iy^i))

由于后面这个显然是个多项式函数，前面的我们之前求导得到为 ((y-y^2))

( herefore displaystyle { ext d^{n+1} yover ext dx^{n+1}}=(y-y^2)(sum_{i=1}^{n+1}ia_iy^{i-1})=(1-y)(sum_{i=1}^{n+1}ia_iy^i))

也就是说，我们可以根据之 (n) 阶导的 (y) 多项式系数，可以 (O(n)) 地推出 ((n+1)) 阶导的 (y) 多项式系数

又 (ecause displaystyle y|_{x=0}={1over 2})

故我们可以 (O(n^2)) 地求解出问题结果

---

用 Ans[MAXN] 储存 (k) 阶导的答案， INV2[MAXN] 储存 (({1over 2})^k mod (10^9+7)) , Frac[MAXN] 储存 (({1over k!})mod (10^9+7)) ， F[MAXN] 储存当前阶导数的 (y) 多项式系数

可以先 (O(n)) 预处理出 (({1over 2})^k mod (10^9+7)) 与 (k!mod (10^9+7))

使用欧拉定理与快速幂 (({1over 1000!})equiv (1000!)^{(10^9+7)-1}(mod 10^9+7)) 求解出 (({1over 1000!})mod (10^9+7)) 。再反向扫回，求解出 (({1over k!})mod (10^9+7))

接下来通过递推算出 F[MAXN] 从而求解 Ans[MAXN]

由于从递推关系可得，式子从 (displaystyle sum_{i=0}^{n+1}a_iy^i) 变为 (displaystyle (1-y)(sum_{i=0}^{n+1}ia_iy^i))

故写出代码：

    for(int i=1;i<=1000;i++){
        for(int j=1;j<=i;j++) F[j]=j*F[j]%MOD;
        for(int j=i+1;j>=1;j--) F[j]=(F[j]-F[j-1]+MOD)%MOD;
    }

再加入刚刚求解的 INV2[MAXN],Frac[MAXN] 即可算出答案

    for(int i=1;i<=1000;i++){
        for(int j=1;j<=i;j++) F[j]=j*F[j]%MOD;
        for(int j=i+1;j>=1;j--) F[j]=(F[j]-F[j-1]+MOD)%MOD;
        for(int j=1;j<=i+1;j++) Ans[i]=(Ans[i]+F[j]*INV2[j]%MOD)%MOD;
        Ans[i]=Ans[i]*Frac[i]%MOD;
    }

最后主函数里面 (O(1)) 查询即可，总复杂度 (O(n^2+T),T) 为询问次数

---

【代码】

#include<iostream>
using namespace std;
typedef long long ll;
const ll MOD=1e9+7,inv2=MOD+1>>1,MAXN=1024;//inv2 为 2 的逆元
ll Ans[MAXN];
inline ll fpow(ll a,ll x){//快速幂
    ll ans=1;
    for(;x;x>>=1,a=a*a%MOD) if(x&1) ans=ans*a%MOD;
    return ans;
}
inline void pre(){
    ll F[MAXN]={0},INV2[MAXN],Frac[MAXN];
    INV2[0]=Frac[0]=1;//初始化
    for(int i=1;i<=1010;i++){//多算一点，防炸
        INV2[i]=INV2[i-1]*inv2%MOD;
        Frac[i]=i*Frac[i-1]%MOD;
    }
    Frac[1000]=fpow(Frac[1000],MOD-2);
    for(int i=999;i>=0;i--)
        Frac[i]=Frac[i+1]*(i+1)%MOD;

    Ans[0]=inv2;//初始化
    F[1]=1;//0 阶导时，y 的多项式为 y，即 1*y^1
    for(int i=1;i<=1000;i++){
        for(int j=1;j<=i;j++) F[j]=j*F[j]%MOD;
        for(int j=i+1;j>=1;j--) F[j]=(F[j]-F[j-1]+MOD)%MOD;
        for(int j=1;j<=i+1;j++) Ans[i]=(Ans[i]+F[j]*INV2[j]%MOD)%MOD;
        Ans[i]=Ans[i]*Frac[i]%MOD;
    }
}
int main(){
    pre();
    int N;
    while(cin>>N)
        cout<<Ans[N]<<endl;
    return 0;
}

---

法二

感谢神仙学弟 Stump 提出的方法

考虑 (displaystyle e^x=sum_{i=0}^infty {1over n!}x^n)

故 (displaystyle e^{-x}=sum_{i=0}^infty {(-1)^nover n!}x^n)

由于询问的最高次数 (kleq 1000) 故我们直接取 (F(x)equiv (e^{-x}+1)pmod {x^{1024}})

则多项式系数 (F[n]equiv ({(-1)^nover n!}+[n==0])pmod {10^9+7})

最后只要令多项式 (G(x)equiv F^{-1}(x)pmod {x^{1024}})

暴力求解出多项式 (G(x)) 即可得出答案：对于每次的询问 (k) ，输出 (G[k]) 。

而求逆可以考虑 (O(n^2)) 暴力方法：

(G(x)equiv F^{-1}(x)pmod x) 时，(G[0]equiv (F[0])^{-1}pmod x)

否则求 (G(x)equiv F^{-1}(x)pmod {x^n}) 时，先递归下去，求解出 (G(x)equiv F^{-1}(x)pmod {x^{n-1}})

而后根据式子 (displaystyle (G*F)[n-1]=sum_{i=0}^{n-1}F[i]G[n-1-i]=F[0]G[n-1]+sum_{i=1}^{n-1}F[i]G[n-1-i]equiv 0pmod {10^9+7})

故 (displaystyle G[n-1]equiv -(F[0])^{-1}cdot sum_{i=1}^{n-1}F[i]G[n-1-i]equiv -G[0]cdot sum_{i=1}^{n-1}F[i]G[n-1-i]pmod {10^9+7})

前半段求 (F(x)) 用预处理阶乘的方法，时间复杂度 (T(n)=O(n)+O(log n)+O(n)=O(n))

后半段求 (G(x)) ，时间复杂度为 (displaystyle T(n)=log n+sum_{i=2}^n sum_{j=1}^i 1=O(n^2))

故总复杂度为 (O(n^2+T))

#include<iostream>
using namespace std;
typedef long long ll;
const ll MOD=1e9+7,MAXN=1024;
inline ll fpow(ll a,ll x) { ll ans=1; for(;x;x>>=1,a=a*a%MOD) if(x&1) ans=ans*a%MOD; return ans; }
inline ll inv(ll a) { return fpow(a,MOD-2); }
inline void inv(ll f[MAXN],ll g[MAXN],ll mod){
    if(mod==1){
        g[0]=inv(f[0]);
        return ;
    }
    inv(f,g,mod-1);
    for(int i=1;i<mod;i++) g[mod-1]=(g[mod-1]+f[i]*g[mod-1-i]%MOD)%MOD;
    g[mod-1]=(MOD-g[mod-1]*g[0]%MOD)%MOD;
}
ll Frac[MAXN],F[MAXN],G[MAXN];
inline void pre(){
    ios::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);

    Frac[0]=1;
    for(int i=1;i<MAXN;i++) Frac[i]=Frac[i-1]*i%MOD;
    F[MAXN-1]=inv(Frac[MAXN-1]);
    for(int i=MAXN-1;i>0;i--) F[i-1]=F[i]*i%MOD;
    for(int i=1;i<MAXN;i+=2) F[i]=MOD-F[i];
    F[0]++;
    inv(F,G,1024);
}
int main(){
    pre();
    int K;
    while(cin>>K)
        cout<<G[K]<<endl;
    return 0;
}

---

法三

由于已知 (e^{-x}+1pmod{x^{1024}}) ，对其进行任意模数的多项式求逆即可

MTT 版

#include<iostream>
#include<cmath>
using namespace std;
typedef long long ll;
typedef long double db;
typedef pair<db,db> pdd;
#define fi first
#define se second
inline pdd operator + (const pdd &a,const pdd &b) { return pdd(a.fi+b.fi,a.se+b.se); }
inline pdd operator - (const pdd &a,const pdd &b) { return pdd(a.fi-b.fi,a.se-b.se); }
inline pdd operator * (const pdd &a,const pdd &b) { return pdd(a.fi*b.fi-a.se*b.se, a.fi*b.se+a.se*b.fi); }
inline pdd operator / (const pdd &a,db b) { return pdd(a.fi/b,a.se/b); }
inline pdd operator ~ (const pdd &p) { return pdd(p.fi,-p.se); }
inline pdd operator / (const pdd &a,const pdd &b) { return a*(~b)/(b*(~b)).fi; }

const int LimBit=12, M=1<<LimBit<<1;
const db pi=acos(-1),eps=1e-9;
int a[M], b[M];
struct FFT{
    int N,na,nb,Rev[M+M],*rev, p, m, InV[M];
    pdd W[2][M+M],*w[2];
    inline ll fpow(ll a,ll x) const { ll ans=1; for(;x;x>>=1,a=a*a%p) if(x&1) ans=ans*a%p; return ans; }
    inline ll inv(ll a) const { return fpow(a,p-2); }
    inline int add(int a,int b) const { return (a+b>=p)?(a+b-p):(a+b); }
    inline int dis(int a,int b) const { return (a-b<0)?(a-b+p):(a-b); }

    FFT(int p_):p(p_){
        InV[1]=1;
        for(int i=2;i<p&&i<M;++i)
            InV[i]=dis(p,(ll)p/i*InV[p%i]%p);

        m=sqrt(p)+1;
        w[0]=W[0]; w[1]=W[1]; rev=Rev;
        for(int d=LimBit;d>=0;--d){
            int N=1<<d;
            for(int i=0;!(i>>d);++i){
                rev[i]=(rev[i>>1]>>1)|((i&1)<<d-1);
                w[0][i]=w[1][i]=pdd( cos(2*pi*i/N) , sin(2*pi*i/N) );
                w[1][i].se=-w[1][i].se;
            }
            w[0]+=N; w[1]+=N; rev+=N;
        }
    }
    inline void work(){
        w[0]=W[0]; w[1]=W[1]; rev=Rev;
        for(int d=LimBit;1<<d>N;--d)
            w[0]+=1<<d, w[1]+=1<<d, rev+=1<<d;
    }
    inline void fft(pdd *a,int f){
        static pdd x;
        for(int i=0;i<N;++i) if(i<rev[i]) swap(a[i],a[rev[i]]);
        for(int i=1;i<N;i<<=1)
            for(int j=0,t=N/(i<<1);j<N;j+=i<<1)
                for(int k=0,l=0;k<i;++k,l+=t)
                    x=w[f][l]*a[j+k+i], a[j+k+i]=a[j+k]-x, a[j+k]=a[j+k]+x;
        if(f) for(int i=0;i<N;++i) a[i]=a[i]/N;
    }
    inline void doit(int *a,int *b,int na,int nb){
        static pdd x, y, p0[M], p1[M];
        for(N=1;N<na+nb-1;N<<=1);
        for(int i=na;i<N;++i) a[i]=0;
        for(int i=nb;i<N;++i) b[i]=0;
        for(int i=0;i<N;++i)
            p1[i]=pdd(a[i]/m,a[i]%m),p0[i]=pdd(b[i]/m,b[i]%m);
        work(); fft(p1,0); fft(p0,0);
        for(int i=0,j;i+i<=N;++i){
            j=(N-1)&(N-i);
            x=p0[i], y=p0[j];
            p0[i]=(x-(~y))*p1[i]/pdd(0,2);
            p1[i]=(x+(~y))*p1[i]/pdd(2,0);
            if(i==j) continue;
            p0[j]=(y-(~x))*p1[j]/pdd(0,2);
            p1[j]=(y+(~x))*p1[j]/pdd(2,0);
        }
        fft(p1,1); fft(p0,1);
        for(int i=0;i<N;++i)
            a[i]=((((ll)(p1[i].fi+0.5+eps)*m%p+(ll)(p1[i].se+0.5+eps)+(ll)(p0[i].fi+0.5+eps))%p*m%p+(ll)(p0[i].se+0.5+eps))%p+p)%p;
    }
    inline void get_inv(int *f,int *g,int n){
        for(int i=1;i<n;++i) g[i]=0; g[0]=inv(f[0]);
        for(int l=1;l<n;l<<=1){
            for(int i=l;i<l<<1;++i) g[i]=0;
            for(int i=0;i<l<<1;++i) a[i]=f[i];
            doit(a,g,l<<1,l<<1);
            for(int i=0;i<l<<1;++i) a[i]=dis(0,a[i]); a[0]=add(a[0],2);
            doit(g,a,l<<1,l<<1);
        }
        for(int i=n;i<N;++i) g[i]=0;
    }
}*fft;

int n,p=1e9+7,f[M],g[M];
inline void pre(){
    f[0]=1;
    for(int i=1;i<1024;++i)
        f[i]=(ll)f[i-1]*i%p;
    f[1023]=fft->inv(f[1023]);
    for(int i=1023;i;--i)
        f[i-1]=(ll)f[i]*i%p;
    for(int i=1;i<1024;i+=2)
        f[i]=fft->dis(0,f[i]);
    ++f[0];
}
int main(){
    ios::sync_with_stdio(0);
    cin.tie(0); cout.tie(0);
    fft=new FFT(p);
    pre();

    fft->get_inv(f,g,1024);
    while(cin>>n) cout<<g[n]<<"
";

    cout.flush();
    return 0;
}

三模 NTT 版