2018-01-22 hihoCoder 1388 FFT 快速傅里叶变换

题目描述:

Profess X is an expert in signal processing. He has a device which can send a particular 1 second signal repeatedly. The signal is A0 ...An-1 under n Hz sampling.
One day, the device fell on the ground accidentally. Profess X wanted to check whether the device can still work properly. So he ran another n Hz sampling to the fallen device and got B0 ...Bn-1.
To compare two periodic signals, Profess X define the DIFFERENCE of signal A and B as follow:

You may assume that two signals are the same if their DIFFERENCE is small enough.
Profess X is too busy to calculate this value. So the calculation is on you.

输入:
The first line contains a single integer T, indicating the number of test cases.
In each test case, the first line contains an integer n. The second line contains n integers, A0 ...An-1. The third line contains n integers, B0 ... Bn-1.
T≤40 including several small test cases and no more than 4 large test cases.
For small test cases, 0<n≤6⋅103.
For large test cases, 0<n≤6⋅104.
For all test cases, 0≤Ai,Bi<220.

输出:
For each test case, print the answer in a single line.

样例输入:
2
9
3 0 1 4 1 5 9 2 6
5 3 5 8 9 7 9 3 2
5
1 2 3 4 5
2 3 4 5 1

样例输出:
80
0

题目大意, 呃.....

给你举个例子吧, 比如给你两组数列A和b, 都是n个数, 这里我们简单点描述, n取4, 编号从0开始, 到n-1结束.

然后有以下几种组合, 问你这四种组合里面最小的组合的值是多少.

组合0: ( A0 - B0 )^2 + ( A1 - B1 )^2 + ( A2 - B2 )^2 + ( A3 - B3 )^2 (对应原题面的k=0)
组合1: ( A0 - B1 )^2 + ( A1 - B2 )^2 + ( A2 - B3 )^2 + ( A3 - B0 )^2 (对应原题面的k=1)
组合2: ( A0 - B2 )^2 + ( A1 - B3 )^2 + ( A2 - B0 )^2 + ( A3 - B1 )^2 (对应原题面的k=2)
组合3: ( A0 - B3 )^2 + ( A1 - B0 )^2 + ( A2 - B1 )^2 + ( A3 - B2 )^2 (对应原题面的k=3, 也就是最大的k=n-1=3)

如果n取任意数的话, 用形象一点的描述的话, 就是, 固定住A的各项, 然后不断地 "滚动" B的各项.

就像上面举的例子, 一开始B的各项的顺序是0123, 滚动一下变成1230, 在滚动一下变成2301, 再滚动一下变成3012, 再滚就滚回0123了.

现在我们以5为例讲一下这个题是怎么跟FFT扯上关系的吧.

n=5, 可以 "滚" 5次,那就是5个组合.

组合0: ( A0 - B0 )^2 + ( A1 - B1 )^2 + ( A2 - B2 )^2 + ( A3 - B3 )^2 + ( A4 - B4 )^2 (对应原题面的k=0)
组合1: ( A0 - B1 )^2 + ( A1 - B2 )^2 + ( A2 - B3 )^2 + ( A3 - B4 )^2 + ( A3 - B0 )^2 (对应原题面的k=1)
组合2: ( A0 - B2 )^2 + ( A1 - B3 )^2 + ( A2 - B4 )^2 + ( A3 - B0 )^2 + ( A3 - B1 )^2 (对应原题面的k=2)
组合3: ( A0 - B3 )^2 + ( A1 - B4 )^2 + ( A2 - B0 )^2 + ( A3 - B1 )^2 + ( A3 - B2 )^2 (对应原题面的k=3)
组合4: ( A0 - B4 )^2 + ( A1 - B0 )^2 + ( A2 - B1 )^2 + ( A3 - B2 )^2 + ( A3 - B3 )^2 (对应原题面的k=4, 也就是最大的k=n-1=4)

如果打开的的话, 那么会发现这五种组合都有共同的10项, 那就是

A0^2 + A1^2 + A2^2 + A3^2 + A4^2 + B0^2 + B1^2 + B2^2 + B3^2 + B4^2

那这个时候就相当于比较这些东西谁比较大了(原本的式子减去公共的10项再取个负就是下面的东西)

组合0' :  A0*B0 *2 + A1*B1 *2 + A2*B2 *2 + A3*B3 *2 + A4*B4 *2
组合1' : A0*B1 *2 + A1*B2 *2 + A2*B3 *2 + A3*B4 *2 + A4*B0 *2
组合2' : A0*B2 *2 + A1*B3 *2 + A2*B4 *2 + A3*B0 *2 + A4*B1 *2
组合3' : A0*B3 *2 + A1*B4 *2 + A2*B0 *2 + A3*B1 *2 + A4*B2 *2
组合4' : A0*B4 *2 + A1*B0 *2 + A2*B1 *2 + A3*B2 *2 + A4*B3 *2

这些东西怎么算呢?

上图!(主要是电脑的公式编辑器我不怎么会用, 所以就干脆手写了)(加了滤镜把一些有的没的的东西去掉了)


啊对了, 我字丑, 下面的十六个字是 项的次数 对应系数 项的次数 对应系数

有没有看出点门道?

组合0' 正好是x^4的系数和x^9的系数之和 恰巧9=5+4
组合1' 正好是x^3的系数和x^8的系数之和 恰巧8=5+3
组合2' 正好是x^2的系数和x^7的系数之和 恰巧7=5+2
组合3' 正好是x^1的系数和x^6的系数之和 恰巧6=5+1
组合4' 正好是x^0的系数和x^5的系数之和 恰巧5=5+0

还不会找规律我觉得你不适合写这个.....

而图里面的最上面的多项式相乘, 可以用FFT进行优化, 找到之后生成的生成函数系数数组就是图上面的东西, 挨个儿看过去哪个x^i的系数和x^(n+i)的系数之和最大就行, 那我们也找到了最大的组合(n-1-i)', 也就找到了最小的组合(n-1-i).

比如我发现x^2的系数和x^(n+2)的系数之和最大, 那就是组合(n-1-2)'是最大的, 那么组合(n-1-2)就是最小的.

行了 , 比如我找到了某个组合是最大的, 但是, 我要告诫你, 由于某些原因, 如果你直接把找到的最大的组合某撇的值取负然后加上A B两个数组各项平方之和的当结果的话, 你一定会WA, 因为精度不够.

那重算就行啦, 不知道你有没有意识到我刻意设置组合标号的含义---"滚"B数组的次数, 这个时候不难了吧, 比如我想求滚了num次之后的组合对应的值, A[i]要减掉的是B[(i+num)%n]  这个模可以理解成从末尾又滚回去了.

然后求n下和, 就是结果

代码如下:

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <algorithm>
#include <iostream>
using namespace std;
const double PI=acos(-1);//这个是一定要有的, 因为fft有应用到三角函数,因此头文件math.h也要包含进去
const int maxx=300000;

struct complex
{
    double a,b;//a表示实部,b表示虚部
    complex(double aa=0.0,double bb=0.0)
    {
        a=aa;
        b=bb;
    }
    complex operator +(const complex &e)
    {//重载运算符+用以直接进行虚数相加,下同
        return complex(a+e.a,b+e.b);
    }
    complex operator -(const complex &e)
    {//重载运算符-
        return complex(a-e.a,b-e.b);
    }
    complex operator *(const complex &e)
    {//重载运算符*
        return complex(a*e.a-b*e.b,a*e.b+b*e.a);
    }
};

void change(complex *y,long long int len)
{//change这个词还挺常用于给自己写的功能函数重名的,要注意一下,不过有可能以后更新模板就换掉了
    long long 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>>=1;
        }
        if(j<k)
        {
            j+=k;
        }
    }
}
void fft(complex *y,long long int len,int on)
{
    change(y, len);
    int i,j,k;
    for(i=2;i<=len;i<<=1)
    {
        complex wn(cos(-on*2*PI/i),sin(-on*2*PI/i));
        for(j=0;j<len;j+=i)
        {
            complex w(1,0);
            for(k=j;k<j+i/2;k++)
            {
                complex u=y[k],t=w*y[k+i/2];
                y[k]=u+t;
                y[k+i/2]=u-t;
                w=w*wn;
            }
        }
    }
    if(on==-1)//还原回原来的多项式
        for(i=0;i<len;i++)
            y[i].a/=len;
}
complex a[maxx],b[maxx],ans[maxx];
long long int A[100000],B[100000],ab[maxx];
int main()
{
    long long int T,i,n,len_a,len,num,res;
    scanf("%lld",&T);
    while(T--)
    {
        memset(a,0,sizeof(a));
        memset(b,0,sizeof(b));
        memset(A,0,sizeof(A));
        memset(B,0,sizeof(B));
        memset(ab,0,sizeof(ab));
        memset(ans,0,sizeof(ans));
        scanf("%lld",&n);
        for(i=0;i<n;i++)
            scanf("%lld",&A[i]);
        for(i=0;i<n;i++)
            scanf("%lld",&B[n-1-i]);//注意这里是倒置哦!
        len=1;
        len_a=n;
        while(len<2*len_a)
            len=len<<1;
        for(i=0;i<n;i++)
        {
            a[i].a=A[i];
            b[i].a=B[i];
        }
        fft(a,len,1);
        fft(b,len,1);
        for(i=0;i<len;i++)
            ans[i]=a[i]*b[i];
        fft(ans,len,-1);
        for(i=0;i<len;i++)
            ab[i]=(long long int)(ans[i].a+0.5);
        res=0;
        num=0;
        for(i=0;i<n;i++)
        {
            if(res<(ab[i]+ab[i+n]))//为什么这么写参见上面的解释
            {
                res=ab[i]+ab[i+n];
                num=n-1-i;
            }
        }
        res=0;
        for(i=0;i<n;i++)
            res+=(A[i]-B[n-1-((i+num)%n)])*(A[i]-B[n-1-((i+num)%n)]);//因为是倒置输入的,取B[(i+num)%n]实际上在代码里面得写B[n-1-((i+num)%n)]
        printf("%lld\n",res);
    }
    return 0;
}


评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值