??

FFT

装13神器

Posted by JU on August 19, 2018

LuoguP3803

#pragma GCC optimize("Ofast,no-stack-protector")
#include <bits/stdc++.h>
#define pr(p) printf("%d ",p)
const int oo=2139063143;
const int N=3010000;
const int mod=1000000007;
using namespace std;
typedef long long LL;
typedef double db;
inline void sc (int &x)
{
    x=0; static int p; p=1; static char c; c=getchar();
    while (!isdigit(c)) { if (c=='-') p=-1; c=getchar(); }
    while ( isdigit(c)) { x=(x<<1)+(x<<3)+(c-48); c=getchar(); }
    x*=p;
}
#define PI 3.1415926535897932384626
typedef complex <db> cp;
int rev[N];
void FFT (cp *a,int n,int kind)
{
    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)
    {
        cp wn=exp (cp(0,kind*PI/i));
        for (int j=0; j< n; j+=i<<1)
        {
            cp wnk(1,0);
            for (int k=j; k< i+j; k++)
            {
                cp x=a[k],y=wnk*a[k+i];//remember to multiply wnk
                a[k]=x+y;
                a[k+i]=x-y;
                wnk*=wn;
            }
        }
    }
    if (kind==-1)
        for (int i=0; i< n; i++)
            a[i]/=n;
}
int n,m;
int init ()
{
    int s=2,bit=1;
    for (; (1<<bit)<=n+m; bit++)
        s<<=1;
    for (int i=0; i< (1<<bit); i++)
        rev[i]=(rev[i>>1]>>1) | ((i&1)<<(bit-1));
    return s;
}
cp a[N],b[N];
int ans[N];
int main()
{
    //freopen (".in","r",stdin);
    //freopen (".out","w",stdout);
    sc(n),sc(m); ++n,++m;
    int sb;
    for (int i=0; i< n; i++)
        sc(sb),a[n-i-1]=sb;
    for (int i=0; i< m; i++)
        sc(sb),b[m-i-1]=sb;
    int s=init ();
    FFT (a,s,1),FFT (b,s,1);
    for (int i=0; i< s; i++)
        a[i]*=b[i];
    FFT (a,s,-1);
    for (int i=0; i< s; i++)
        ans[i]=(int)(a[i].real ()+0.5);
    int len=n+m;
    while (!ans[len]&&len>=0) --len;
    if (len==-1) pr(0);
    for (int i=len; i>=0; i--)
        pr(ans[i]);
    return 0;
}

预处理单位根

有时候用来避免爆精度。
还是Luogu3803（不知为何预处理之后反而变得奇慢无比）

const double PI=acos (-1);
typedef complex <db> cp;
int rev[N];
cp w[N];
void FFT (cp *a,int n,int kind)
{
    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; j< n; j+=i<<1)
        {
            int cur=0;
            for (int k=j; k< i+j; k++)
            {
                cp x=a[k],y=w[cur]*a[k+i];
                a[k]=x+y;
                a[k+i]=x-y;
                cur=(cur+kind*n/(i<<1)+n)&(n-1);
            }
        }
    if (kind==-1)
        for (int i=0; i< n; i++)
            a[i]/=n;
}
int n,m;
int init ()
{
    int bit=1,s=2;
    for (; (1<<bit)<=n+m; bit++)
        s<<=1;
    for (int i=0; i< (1<<bit); i++)
        rev[i]=(rev[i>>1]>>1) | ((i&1)<<(bit-1));
    for (int i=0; i< s; i++)
        w[i]=cp(cos(2*PI*i/s),sin(2*PI*i/s));
    return s;
}

求卷积

个人理解，卷积大概是：

\[\sum_{i=0}^{n} f[i]g[n-i]\]

栗:LuoguP3338

\[Ei=\sum_{i>j}^{}{q_j \over (i-j)^2}-\sum_{i<j}^{}{q_j \over (j-i)^2}\] \[Ei=\sum_{j=1}^{i-1}{q_j \over (i-j)^2}-\sum_{j=i+1}^{n}{q_j \over (j-i)^2}\] \[Ei=\sum_{j=1}^{i-1}{a[j]b[i-j]}-\sum_{j=i+1}^{n}{a[j]b[j-i]}\] \[a[i]=q_i,b[i]={1 \over i^2}\]

第一个式子已经化成了卷积的形式，
对于第二个式子可以换一下元 。
$a’$表示$a$反向后的值。

\[\sum_{j=1}^{n-i}{a'[i]b[j]}\] 
#include <bits/stdc++.h>
#define pr(p) printf("%.3lf\n",p)
const int oo=2139063143;
const int N=1010000;
const int mod=1000000007;
using namespace std;
typedef long long LL;
typedef double db;
inline void sc (int &x)
{
    x=0; static int p; p=1; static char c; c=getchar();
    while (!isdigit(c)) { if (c=='-') p=-1; c=getchar(); }
    while ( isdigit(c)) { x=(x<<1)+(x<<3)+(c-48); c=getchar(); }
    x*=p;
}
#define PI 3.1415926535897932384626
typedef complex <db> cp;
int rev[N];
void FFT (cp *a,int n,int kind)
{
    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)
    {
        cp wn=exp (cp(0,kind*PI/i));
        for (int j=0; j< n; j+=i<<1)
        {
            cp wnk(1,0);
            for (int k=j; k< i+j; k++)
            {
                cp x=a[k],y=wnk*a[k+i];
                a[k]=x+y;
                a[k+i]=x-y;
                wnk*=wn;
            }
        }
    }
    if (kind==-1)
        for (int i=0; i< n; i++)
            a[i]/=n;
}
int n;
int init ()
{
    int s=2,bit=1;
    for (; (1<<bit)<=(n<<1); bit++)
        s<<=1;
    for (int i=0; i< (1<<bit); i++)
        rev[i]=(rev[i>>1]>>1) | ((i&1)<<(bit-1));
    return s;
}
cp a[N],b[N];
db q[N],ans[N];
int main()
{
    //freopen (".in","r",stdin);
    //freopen (".out","w",stdout);
    sc(n);
    int s=init ();

    for (int i=1; i<=n; i++)
        scanf ("%lf",&q[i]),a[i]=q[i];
    for (int i=1; i<=n; i++)
        b[i]=1.0/i/i;
    FFT (a,s,1),FFT (b,s,1);
    for (int i=0; i< s; i++)
        a[i]*=b[i];
    FFT (a,s,-1);
    for (int i=1; i<=n; i++)
        ans[i]=a[i].real ();

    memset (a,0,sizeof (a));
    for (int i=1; i<=n; i++)
        a[i]=q[n-i+1];
    FFT (a,s,1);
    for (int i=0; i< s; i++)
        a[i]*=b[i];
    FFT (a,s,-1);
    for (int i=1; i<=n; i++)
        ans[i]-=a[n-i+1].real ();
    for (int i=1; i<=n; i++)
        pr(ans[i]);
    return 0;
}
CC 原创文章采用CC BY-NC-SA 4.0协议进行许可，转载请注明：“转载自：FFT”
CATALOG
    FRIENDS