• 伯努利数学习笔记&&Luogu P3711 仓鼠的数学题


    新科技

    Luogu P3711


    题意

    设$ S_{k,n}$表示$ displaystylesum_{i=0}^n i^k$

    求多项式$displaystylesum_{k=0}^n S_{k,x}a_k$的各项系数

    数组$ a$给定,$ n leq 100000$


    伯努利数 

    伯努利数$B$是一个数列,满足

    $$sum_{i=0}^n B_iinom{n+1}{i}=0$$

    可以用它来求自然数幂和

    $$ S_{k,n-1}=sum_{i=0}^{n-1}i^k=frac{1}{k+1}sum_{i=0}^kinom{k+1}{i}B_in^{k+1-i}$$

    如果已经得到了数列$ B$,求自然数幂和$S_{k,n}$是$ O(k)$的

    直接根据定义可以$ O(n^2)$递推伯努利数,考虑更快速的推法

    $$
    egin{aligned}
    sum_{i=0}^n B_iinom{n+1}{i}&=0\
    sum_{i=0}^{n-1} B_iinom{n}{i}&=0 (n>1)\
    B_n+sum_{i=0}^{n-1} B_iinom{n}{i}&=B_n (n>1)\
    B_n&=sum_{i=0}^nB_iinom{n}{i} (n>1)\
    frac{B_n}{n!}&=sum_{i=0}^nfrac{B_i}{i!(n-i)!}\
    end{aligned}
    $$

    设伯努利数的指数型生成函数为$ B$,伯努利数的第一项$ B_1=-frac{1}{2}$

    则有$B*e^x=B+x$

    整理得$B=frac{x}{e^x-1}=(frac{e^x-1}{x})^{-1}$

    直接多项式求逆即可

    时间复杂度$ O(n log n)$


    回到原题

    用伯努利数展开得

    $$
    egin{aligned}
    ans&=sum_{k=0}^na_k S_{k,x}\
    &=sum_{k=0}^na_k(x^k+frac{1}{k+1}sum_{i=0}^kinom{k+1}{i}B_ix^{k+1-i})\
    &=(sum_{k=0}^na_kx^k)+(sum_{k=0}^nk!sum_{i=0}^kfrac{B_i}{i!(k+1-i)!}x^{k+1-i})\
    ans[x^d]&=a_d+sum_{i=0}^{n+1}frac{B_i}{d!i!}(d+i-1)!\
    frac{ans[x^d]}{d!}&=a_d+sum_{i=0}^{n+1}frac{B_i}{i!}(d+i-1)!
    end{aligned}
    $$

    发现这是一个差卷积的形式

    按套路反转之后$ NTT$即可

    总复杂度仍是$ O(n log n)$


    代码 

    #include<ctime>
    #include<cmath>
    #include<cstdio>
    #include<cstring>
    #include<iostream>
    #include<algorithm>
    #include<queue>
    #include<vector>
    #define p 998244353
    #define rt register int
    #define ll long long
    #define ull unsigned long long
    using namespace std;
    inline ll read(){
        ll x=0;char zf=1;char ch=getchar();
        while(ch!='-'&&!isdigit(ch))ch=getchar();
        if(ch=='-')zf=-1,ch=getchar();
        while(isdigit(ch))x=x*10+ch-'0',ch=getchar();return x*zf;
    }
    void write(ll y){if(y<0)putchar('-'),y=-y;if(y>9)write(y/10);putchar(y%10+48);}
    void writeln(const ll y){write(y);putchar('
    ');}
    int k,m,n,x,y,z,cnt,ans;
    
    namespace Poly{
        #define poly vector<int>
        #define MAXN 524288
        int ksm(int x,int y=p-2){
            int ans=1;
            for(;y;y>>=1,x=1ll*x*x%p)if(y&1)ans=1ll*ans*x%p;
            return ans;
        }    
        void NTT(int n,poly &A,int fla){
            static ull F[MAXN],W[MAXN];A.resize(n);
            for(rt i=0,j=0;i<n;i++){
                F[i]=A[j];
                for(rt k=n>>1;(j^=k)<k;k>>=1);
            }
            for(rt i=1;i<n;i<<=1){
                const int w=W[1]=ksm(3,(p-1)/2/i);W[0]=1;
                for(rt k=2;k<i;k++)W[k]=1ll*W[k-1]*w%p;
                for(rt j=0;j<n;j+=i<<1){
                    for(rt k=0;k<i;k++){
                        const ull x=F[j+k],y=F[i+j+k]*W[k]%p;
                        F[j+k]=x+y,F[i+j+k]=x+p-y;
                    }
                }
            }
            for(rt i=0;i<n;i++)A[i]=F[i]%p;
            if(fla==-1){
                const int invn=ksm(n);
                reverse(A.begin()+1,A.end());
                for(rt i=0;i<n;i++)A[i]=1ll*A[i]*invn%p;
            }
        }
        poly Mul(poly x,poly y){
            int sz=x.size()+y.size()-1,lim=1;
            while(lim<=sz)lim<<=1;
            NTT(lim,x,1);NTT(lim,y,1);
            for(rt i=0;i<lim;i++)x[i]=1ll*x[i]*y[i]%p;
            NTT(lim,x,-1);x.resize(sz);return x;
        }
        poly Inv(poly a,int n=-1){
            if(n==-1)n=a.size();
            if(n==1)return {ksm(a[0])};
            poly c=Inv(a,n+1>>1),d(&a[0],&a[n]);
            int lim=1;while(lim<=n*2)lim<<=1;
            NTT(lim,c,1);NTT(lim,d,1);
            for(rt i=0;i<lim;i++)c[i]=1ll*c[i]*(2ll+p-1ll*d[i]*c[i]%p)%p;
            NTT(lim,c,-1);c.resize(n);return c;
        }
    }
    using namespace Poly;
    int inv[250010],jc[250010],njc[250010],a[250010];
    poly B;
    void init(int k){
        for(rt i=0;i<=1;i++)inv[i]=jc[i]=njc[i]=1;
        for(rt i=2;i<=k+2;i++){
            inv[i]=1ll*inv[p%i]*(p-p/i)%p;
            jc[i]=1ll*jc[i-1]*i%p;
            njc[i]=1ll*njc[i-1]*inv[i]%p;
        }
        B.resize(k+1);
        for(rt i=0;i<=k;i++)B[i]=njc[i+1];
        B=Inv(B);
        for(rt i=0;i<=k;i++)B[i]=1ll*B[i]*jc[i]%p;
    }
    int main(){
        n=read();init(n+2);
        for(rt i=0;i<=n;i++)a[i]=read();
        poly ans(n+2),C(n+1);
        for(rt i=0;i<=n;i++)B[i]=1ll*B[i]*njc[i]%p;
        for(rt i=0;i<=n;i++)C[i]=1ll*jc[i]*a[i]%p;
        reverse(&B[0],&B[n+1]);B.resize(n+1);C.resize(n+1);
        ans=Mul(B,C);
        for(rt i=0;i<=n+1;i++)ans[n+i-1]=1ll*ans[n+i-1]*njc[i]%p;
        for(rt i=0;i<=n;i++)(ans[n+i-1]+=a[i])%=p;write(a[0]),putchar(' ');
        for(rt i=1;i<=n+1;i++)write(ans[n+i-1]),putchar(' ');
        return 0;
    }
  • 相关阅读:
    Oracle数据库-建库、建表空间,建用户
    创建oracle数据库的表空间、用户、目录、导入导出文件等信息
    使用 UML 进行业务建模:理解业务用例与系统用例的相似和不同之处
    DataSet用法详细
    DataSet和DataTable详解
    经典SQL语句集锦
    JS脚本验证大全
    用Unity写一个12306验证器的恶搞图生成软件
    SQL Server中 sysobjects、syscolumns、systypes
    如何在oracle中导入导出(备份&恢复)dmp数据库文件
  • 原文地址:https://www.cnblogs.com/DreamlessDreams/p/10525098.html
Copyright © 2020-2023  润新知