• FFT,NTT 专题


    学习傅里叶的基本性质及其代码,可以参考大神理解

    还有 ACdream 的博客

    贴一下NTT的模板:

    using namespace std;
    typedef long long ll;
    
    int n;
    const ll MOD=998244353;
    const int N = 400005;
    const int g=3;
    int len;
    ll A[N];
    long long a[N],b[N],wn[30];
    long long q_pow(long long x,long long y,long long P)
    {
        long long ans=1;
        while(y>0)
        {
            if(y&1)ans=ans*x%P;
            x=x*x%P;
            y>>=1;
        }
        return ans;
    }
    void init()
    {
        for(int i=0;i<21;i++)
        {
            int t=1<<i;
            wn[i]=q_pow(g,(MOD-1)/t,MOD);
        }
    }
    ///雷德算法,2^M=len,将第i位的数与“i的二进制反转之后的位”的数交换
    void rader(long long F[],int len)
    {
        int j=len/2;///模拟二进制反转进位的的位置
        for(int i=1;i<len-1;i++)
        {
            if(i<j)swap(F[i],F[j]);///该出手时就出手
            int k=len/2;
            while(j>=k)
            {
                j-=k;
                k>>=1;
            }
            if(j<k)j+=k;
        }
    }
    void NTT(long long F[],int len,int t)
    {
        int id=0;
        rader(F,len);
        for(int h=2;h<=len;h<<=1)
        {
            id++;
            for(int j=0;j<len;j+=h)
            {
                long long E=1;
                for(int k=j;k<j+h/2;k++)
                {
                    long long u=F[k];
                    long long v=(E*F[k+h/2])%MOD;
                    F[k]=(u+v)%MOD;
                    F[k+h/2]=((u-v)%MOD+MOD)%MOD;
                    E=(E*wn[id])%MOD;
                }
            }
        }
        if(t==-1)
        {
            for(int i=1;i<len/2;i++)swap(F[i],F[len-i]);
            long long ni=q_pow(len,MOD-2,MOD);
            for(int i=0;i<len;i++)F[i]=(F[i]%MOD*ni)%MOD;
        }
    }
    
    void work()///卷积,点乘,插值
    {
        NTT(a,len,1);
        NTT(b,len,1);
        for(int i=0;i<len;i++)
            a[i]=(a[i]*b[i])%MOD;
        NTT(a,len,-1);
    }
    int main()
    {
        #ifndef ONLINE_JUDGE
        freopen("in.txt","r",stdin);
        #endif
        init();
        int T_T;
        scanf("%d",&T_T);
        for(int kase=1;kase<=T_T;kase++)
        {
            sc(n);
            len=1;
            while(len<=2*n)
                len<<=1;
        /**
           这部分就是你对公式的变形并把他转化为卷积形式,分别存在a[],b[]里
        */
        work();
      }
      return 0;
    }

    对于傅里叶运用的要点,要认真对待对下面卷积公式的理解

    训练:

    传送门:hdu 5829 Rikka with Subset

    题意:给你一个数组A[],然后计算F[k],F[k]指A[]所有子集中,先对每个子集前k大的数的和,然后再对把所以子集所求值求和

    思路:因为我也是初学,所以我找的一个大神的博客学习的http://blog.csdn.net/cqu_hyx/article/details/52194696

    对于其中的几点,我做一点补充:

    1. f[k]计算的是第k大的数对答案的贡献,而比他大的数的贡献没有计算;不过只需要累加分f[1]..f[k-1]就是前k大的数的和了
    2. 对比以上卷积公式,我们得到的是   a[i]=2n-i/i!   ,  b[i]=A[i]*(i-1)! ;如果我们b[0]=A[1]*0!,那么我们reverse一下,b[i]=b[n-i]了;在我们推导的f[k]公式里,f[k]是与a[i]*b[i+k]相关的,且i最大为(n-k),经过我们对b数组reverse,f[k]是与a[i]*b[(n-k)-i]相关的,这刚好和上面的卷积公式一致,所以我们可以放心ntt了。
    3. 由卷积公式可知,y[n-k]是a[i]*b[(n-k)-i]的卷积结果,而这个结果存在a[]里;并且f[k]也等于a[i]*b[(n-k)-i]的卷积,所以f[k]与y[n-k]即a[n-k]相关
    /**************************************************************
        Problem:hdu 5829 Rikka with Subset
        User: youmi
        Language: C++
        Result: Accepted
        Time:2605MS
        Memory:11804K
    ****************************************************************/
    //#pragma comment(linker, "/STACK:1024000000,1024000000")
    //#include<bits/stdc++.h>
    #include <iostream>
    #include <cstdio>
    #include <cstring>
    #include <algorithm>
    #include <map>
    #include <stack>
    #include <set>
    #include <sstream>
    #include <cmath>
    #include <queue>
    #include <deque>
    #include <string>
    #include <vector>
    #define zeros(a) memset(a,0,sizeof(a))
    #define ones(a) memset(a,-1,sizeof(a))
    #define sc(a) scanf("%d",&a)
    #define sc2(a,b) scanf("%d%d",&a,&b)
    #define sc3(a,b,c) scanf("%d%d%d",&a,&b,&c)
    #define scs(a) scanf("%s",a)
    #define sclld(a) scanf("%I64d",&a)
    #define pt(a) printf("%d
    ",a)
    #define ptlld(a) printf("%I64d
    ",a)
    #define rep(i,from,to) for(int i=from;i<=to;i++)
    #define irep(i,to,from) for(int i=to;i>=from;i--)
    #define Max(a,b) ((a)>(b)?(a):(b))
    #define Min(a,b) ((a)<(b)?(a):(b))
    #define lson (step<<1)
    #define rson (lson+1)
    #define eps 1e-6
    #define oo 0x3fffffff
    #define TEST cout<<"*************************"<<endl
    const double pi=4*atan(1.0);
    
    using namespace std;
    typedef long long ll;
    
    int n;
    const ll MOD=998244353;
    const int N = 400005;
    const int g=3;
    int len;
    ll inv2;
    ll A[N];
    ll inv[N],fac[N],f[N],bit[N];
    long long a[N],b[N],wn[30];
    long long q_pow(long long x,long long y,long long P)
    {
        long long ans=1;
        while(y>0)
        {
            if(y&1)ans=ans*x%P;
            x=x*x%P;
            y>>=1;
        }
        return ans;
    }
    void init()
    {
        for(int i=0;i<21;i++)
        {
            int t=1<<i;
            wn[i]=q_pow(g,(MOD-1)/t,MOD);
        }
        inv[0]=inv[1]=1;
        rep(i,2,1e5)
            inv[i]=(inv[i-1]*q_pow(i,MOD-2,MOD))%MOD;
        fac[0]=fac[1]=1;
        rep(i,2,1e5)
            fac[i]=(fac[i-1]*i)%MOD;
        bit[0]=1;
        rep(i,1,1e5)
            bit[i]=bit[i-1]*2%MOD;
        inv2=q_pow(2,MOD-2,MOD);
    }
    ///雷德算法,2^M=len,将第i位的数与“i的二进制反转之后的位”的数交换
    void rader(long long F[],int len)
    {
        int j=len/2;///模拟二进制反转进位的的位置
        for(int i=1;i<len-1;i++)
        {
            if(i<j)swap(F[i],F[j]);///该出手时就出手
            int k=len/2;
            while(j>=k)
            {
                j-=k;
                k>>=1;
            }
            if(j<k)j+=k;
        }
    }
    void NTT(long long F[],int len,int t)
    {
        int id=0;
        rader(F,len);
        for(int h=2;h<=len;h<<=1)
        {
            id++;
            for(int j=0;j<len;j+=h)
            {
                long long E=1;
                for(int k=j;k<j+h/2;k++)
                {
                    long long u=F[k];
                    long long v=(E*F[k+h/2])%MOD;
                    F[k]=(u+v)%MOD;
                    F[k+h/2]=((u-v)%MOD+MOD)%MOD;
                    E=(E*wn[id])%MOD;
                }
            }
        }
        if(t==-1)
        {
            for(int i=1;i<len/2;i++)swap(F[i],F[len-i]);
            long long ni=q_pow(len,MOD-2,MOD);
            for(int i=0;i<len;i++)F[i]=(F[i]%MOD*ni)%MOD;
        }
    }
    
    void work()///卷积,点乘,插值
    {
        NTT(a,len,1);
        NTT(b,len,1);
        for(int i=0;i<len;i++)
            a[i]=(a[i]*b[i])%MOD;
        NTT(a,len,-1);
    }
    int main()
    {
        #ifndef ONLINE_JUDGE
        freopen("in.txt","r",stdin);
        #endif
        init();
        int T_T;
        scanf("%d",&T_T);
        for(int kase=1;kase<=T_T;kase++)
        {
            sc(n);
            len=1;
            while(len<=2*n)
                len<<=1;
            rep(i,1,n)
                sclld(A[i]);
            sort(A+1,A+1+n,greater<ll>());
            zeros(a);zeros(b);
            rep(i,0,n-1)
                a[i]=(bit[n-i]*inv[i])%MOD;
            rep(i,0,n-1)
                b[i]=(A[i+1]*fac[i])%MOD;
            reverse(b,b+n);
            work();
            ll r=inv2;
            f[0]=0;
            rep(i,1,n)
            {
                f[i]=a[n-i]*inv[i-1]%MOD*r%MOD;
                r=r*inv2%MOD;
                f[i]=(f[i]+f[i-1])%MOD;
                printf("%I64d ",f[i]);
            }
            puts("");
        }
    }
    View Code
    不为失败找借口,只为成功找方法
  • 相关阅读:
    大型网站技术架构-阅读笔记1
    如何发挥一个字节的极限,存储大量内容
    利用easyui创建一个简单的登录页面
    linux tomcat 快捷操作
    linux 安装jdk
    Linux-查看服务器的信息
    HTTP协议(1)
    Linux-ps命令
    Linux-tcpdump命令
    转载-测试新人培训方法之目标法
  • 原文地址:https://www.cnblogs.com/youmi/p/5806449.html
Copyright © 2020-2023  润新知