FZU 2314

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

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

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

---

法一

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

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

不妨设 \(\displaystyle y=f(x)={1\over e^{-x}+1}\Rightarrow {\text dy\over \text dx}=-{{\text d\over \text dx}(e^{-x}+1)\over (e^{-x}+1)^2}={(e^{-x}+1)-1\over (e^{-x}+1)^2}=y-y^2\)

再次求导可得 \(\displaystyle {\text d^2 y\over \text dx^2}={\text d\over \text dx}({\text dy\over \text dx})={\text dy\over \text dx}\cdot {\text d\over \text dy}({\text dy\over \text dx})={\text dy\over \text dx}\cdot {\text d\over \text dy}(y-y^2)=(1-2y)\cdot {\text dy\over \text dx}=(1-2y)(y-y^2)=y-3y^2+2y^3\)

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

则 \(\displaystyle {\text d^{n+1} y\over \text dx^{n+1}}={\text d\over \text dx}({\text d^n y\over \text dx^n})={\text dy\over \text dx}\cdot ({\text d\over \text dy}\sum_{i=1}^{n+1}a_iy^i)\)

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

\(\therefore \displaystyle {\text d^{n+1} y\over \text 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\) 多项式系数

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

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

---

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

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

使用欧拉定理与快速幂 \(({1\over 1000!})\equiv (1000!)^{(10^9+7)-1}(\mod 10^9+7)\) 求解出 \(({1\over 1000!})\mod (10^9+7)\) 。再反向扫回，求解出 \(({1\over 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 {1\over n!}x^n\)

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

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

则多项式系数 \(F[n]\equiv ({(-1)^n\over 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 0\pmod {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}+1\pmod{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]<<"\n";

    cout.flush();
    return 0;
}

三模 NTT 版