FFT

最新推荐文章于 2023-02-02 16:48:41 发布

mwhybmj

最新推荐文章于 2023-02-02 16:48:41 发布

阅读量387

点赞数

分类专栏： mathmatic

本文链接：https://blog.csdn.net/tigercoder/article/details/71512905

版权

mathmatic 专栏收录该内容

1 篇文章 0 订阅

订阅专栏

http://202.121.199.212/JudgeOnline/problem.php?id=1966

https://vjudge.net/problem/UVA-12298

FFT:

卷积：A[] = A[0] , A[1], A[2]……

B[] = B[0], B[1], B[2] ……

A与B的卷积就是C[]

C[0] = A[0] * B[0]

C[1] = A[0] * B[1] + A[1] * B[0]

C[2] = A[0] * B[2] + A[1] * B[1] + A[2] *B[0]

……

那，当我们需要算出两个数组的卷积是就需要用到O(n^2)的时间，但当数据量较大时，我们需要采取更快的方法，这就是fft的作用。

原理：

我们可以将A和B都看成多项式，即：A[] = a₀* x⁰ + a₁ * x¹+ a₂ * x² +…

B也同理。这样求卷积这么一个过程就是在求x的指数相同的系数的和。如a₀*x⁰*b_n*xⁿ=a₀*b_n*xⁿ, 而 a₁*x¹*b₁*x^n-1= a₁*b₁*xⁿ …

FFT的作用是将多项式表达式转换为点值表达式。

这里要补充一个定理：

一个n次多项式可以由（n+1）个点来确定，如1+x²+3*x,需要3个点即可确定。

所以A可以由若干个点确定，B也可以由若干个点确定。那么如果我们可以找到满足C方程的点，即可逆向求出C的多项式表达式。

那么如何找到满足C方程的点呢？由一个定理可知：假如满足A的点为(x1,y1),(x2,y2)…

满足B的点为(a1,b1),(a2,b2)…那么满足C的点就为A*B，即(x1*a1,y1*b1) , (x2*a2,y2*b2)…

这样只需要O(n)的时间即可求出确定C所需要的点。然后我们再逆FFT，即可求出C的多项式表达式，其系数就是对应的卷积。

len开多大合适呢？要对几个数组做fft就开到2的多少次方（小于等于）

#include<iostream>
#include<cstdio>
#include<string.h>
#include<math.h>
#include<string>
#include<map>
#include<set>
#include<vector>
#include<algorithm>
#include<queue>
#include<iomanip>
using namespace std;
const int INF = 0x3f3f3f3f;
const int NINF = 0xc0c0c0c0;
const int maxn = 300100;
const double PI = acos(-1.0);
int  suma[maxn];

int n;

struct complex
{
    double x,y;
    complex(double _x=0.0,double _y=0.0)
    {
        x = _x;
        y = _y;
    }
    complex operator -(const complex &b) const
    {
        return complex(x-b.x,y-b.y);
    }
    complex operator +(const complex &b) const
    {
        return complex(x+b.x,y+b.y);
    }
    complex operator *(const complex &b) const
    {
        return complex(x*b.x-y*b.y,x*b.y+y*b.x);
    }
};
complex a[maxn],b[maxn],c[maxn];

void change(complex y[],int len)
{
    int i,j,k;
    for(i=1,j=len/2;i<len-1;i++){
        if(i<j) swap(y[i],y[j]);
        k = len/2;
        while(j >= k)
        {
            j -= k;
            k /= 2;
        }
        if(j < k) j+=k;
    }
}

void fft(complex y[], int len, int on)
{
    change(y,len);
    for(int h=2;h<=len;h<<=1)
    {
        complex wn(cos(-on*2*PI/h),sin(-on*2*PI/h));
        for(int j=0;j<len;j+=h){
            complex w(1,0);
            for(int k=j;k<j+h/2;k++){
                complex u = y[k];
                complex t = w*y[k+h/2];
                y[k] = u+t;
                y[k+h/2] = u-t;
                w = w*wn;
            }
        }
    }
    if(on == -1){
        for(int i=0;i<len;i++){
            y[i].x /= len;
        }
    }
}
int main()
{
    while(cin >> n){
        for(int i=0;i<n;i++){
            cin >> a[i].x;
        }
        for(int i=0;i<n;i++){
            cin >> b[i].x;
        }
        int l = 1;
        while(l < 2*n) l <<= 1;
        for(int i=n;i<l;i++){                   //attention!   initialize the array everytime!
            a[i] = complex(0.0,0.0);
            b[i] = complex(0.0,0.0);
        }
        for(int i=0;i<n;i++){
            a[i] = complex(a[i].x,0.0);
            b[i] = complex(b[i].x,0.0);
        }
        fft(a,l,1);
        fft(b,l,1);
        for(int i=0;i<l;i++){
            c[i] = a[i]*b[i];                   // maybe excceed the int
        }
        fft(c,l,-1);
        for(int i=0;i<l;i++){                   // jingdu
            suma[i] = (int)(c[i].x+0.5);
        }
        for(int i=0;i<n;i++){
            if(i == 0) cout << suma[i];
            else cout << ' ' << suma[i];
        }
        cout << '\n';
    }
}

#include<iostream>
#include<cstdio>
#include<string.h>
#include<math.h>
#include<string>
#include<map>
#include<set>
#include<vector>
#include<algorithm>
#include<queue>
#include<iomanip>
using namespace std;
const int INF = 0x3f3f3f3f;
const int NINF = 0xc0c0c0c0;
const long double PI = acos(-1.0);
struct complex
{
    long double x,y;
    complex (long double _x = 0.0 , long double _y = 0.0)
    {
        x = _x;
        y = _y;
    }
    complex operator -(const complex &b)const
    {
        return complex(x-b.x,y-b.y);
    }
    complex operator +(const complex &b)const
    {
        return complex(x+b.x,y+b.y);
    }
    complex operator *(const complex &b)const
    {
        return complex(x*b.x-y*b.y,x*b.y+y*b.x);
    }
};

void change(complex y[],int len)
{
    int i,j,k;
    for(i=1,j=len/2;i<len-1;i++)
    {
        if(i<j) swap(y[i],y[j]);
        k = len/2;
        while(j>=k){
            j -= k;
            k /= 2;
        }
        if(j<k) j+= k;
        
    }
}

void fft(complex y[],int len,int on)
{
    change(y,len);
    for(int h=2;h<=len;h<<=1){
        complex wn(cos(-on*2*PI/h),sin(-on*2*PI/h));
        for(int j=0;j<len;j+=h){
            complex w(1,0);
            for(int k=j;k<j+h/2;k++){
                complex u = y[k];
                complex t = w*y[k+h/2];
                y[k] = u+t;
                y[k+h/2] = u-t;
                w = w*wn;
            }
        }
    }
    if(on == -1){
        for(int i=0;i< len;i++){
            y[i].x /= len;
        }
    }
}
const int maxn = 50000*20;
complex S[maxn],H[maxn],C[maxn],D[maxn];
long long  suma[maxn];
bool prime[50010];
int main()
{
    int a,b,c;
    prime[0]=1;
    prime[1]=1;
    for(long long i=2;i<50010;i++){
        if(prime[i] == 0){
            for(long long  j = i; j*i < 50010;j++){
                prime[j*i] = 1;
            }
        }
    }
    while(scanf("%d%d%d",&a,&b,&c)){
        if(a == 0 && b==0 && c ==0){
            break;
        }
        memset(S,0,sizeof(S));
        memset(H,0,sizeof(H));
        memset(C,0,sizeof(C));
        memset(D,0,sizeof(D));
        int len = 1;
        while(len <= 8*b) len <<= 1;        //attention!    
        for(int i=4;i<=b;i++){
            if(prime[i] == 1){
                S[i] = complex(1,0);
                H[i] = complex(1,0);
                C[i] = complex(1,0);
                D[i] = complex(1,0);
            }
            else{
                S[i] = complex(0,0);
                H[i] = complex(0,0);
                C[i] = complex(0,0);
                D[i] = complex(0,0);
            }
        }
        for(int i=b+1;i<len;i++){
            S[i] = complex(0,0);
            H[i] = complex(0,0);
            C[i] = complex(0,0);
            D[i] = complex(0,0);
        }
        for(int i=0;i<c;i++){
            int sumt;
            char temp;
            scanf("%d%c",&sumt,&temp);
            if(temp == 'S'){
                S[sumt] = complex(0,0);
            }
            else if(temp == 'H'){
                H[sumt] = complex(0,0);
            }
            else if(temp == 'C'){
                C[sumt] = complex(0,0);
            }
            else D[sumt] = complex(0,0);
        }
        
        fft(S,len,1);
        fft(H,len,1);
        fft(C,len,1);
        fft(D,len,1);
        for(int i=0;i<len;i++){
            S[i] = S[i]*H[i]*C[i]*D[i];
        }
        fft(S,len,-1);
        for(int i=0;i<len;i++){
            suma[i] = (long long)(S[i].x+0.5);
        }
        for(int i=a;i<=b;i++){
            //if(i<16){
             //   cout << 0 << endl;
              //  continue;
           // }
            printf("%lld\n",suma[i]);
        }     
        printf("\n");
        
    }
}

mwhybmj

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
FFT

http://202.121.199.212/JudgeOnline/problem.php?id=1966#include#include#include#include#include#include#include#include#include#include#includeusing namespace std;const int IN
复制链接

扫一扫

专栏目录