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
不为失败找借口,只为成功找方法
原文地址:https://www.cnblogs.com/youmi/p/5806449.html