【BZOJ2194】快速傅里叶之二,FFT和一点奇怪的想法

传送门(权限题)
题意:给定数列a,b,求 ck=i=kn1aibik , n105
思路:
(偶然发现一点新东西的奇妙感,如果神牛们早就知道了就无视掉我这篇胡扯的博文吧)
下面的n默认为是2的正整数次幂
传统FFT模板题,网上的传统做法是把a或b翻转过来,化成卷积形式然后直接做就可以了
但是我在做这道题的时候并没有想到这种方法,而是一直在考虑在多项式方面的问题,从一般感觉上来说, (ik)+i=k ,那怎么化成负的呢?
我们知道一般形式为 A(x)=a0+a1x+a2x2+...+an1xn1 ,称其为x的一个次数界为n的多项式;对于多项式乘法,设结果为c,相乘的两个多项式为a,b; ai,bi,ci 分别表示a,b,c次数为i的项,那么 ci=j=0iajbij ,这是一个卷积形式,这也是FFT是处理卷积问题的利器的原因。回到原来的问题上,我们从多项式的项次数上考虑,那么我们要求的 ck 不就是”a的次数为i的项乘b的次数为-(i-k)的项”的累和?
抱着试一试与乱搞的心态,我写了一个这样的东西
A(x)=a0+a1x+a2x2+...+an1xn1
B(x)=b0+b1x1+b2x2+...+bn1xn1
C(x)=A(x)B(x1)
那么多项式 C(x) 的0~n-1次项的系数就是答案
处理 B(x) 的FFT的时候,像逆DFT一样,把单位复数根 ωn 变成 ω1n 就好了
但是 C(x) 里有负指数项,正确答案在哪里呢?
抱着弃疗的心态输出了 c0,c1..cn1 ,发现样例过了,对拍过了,然后就A了
静下心来仔细思考,我们在算出A,B的DFT,再做完向量乘法时,实际上的C是这样的
ci=0<j<n,0<j+i<najbj+i
C(x)=c(n1)x(n1)+c(n2)x(n2)+..+c0+c1x+c2x2+...cn1xn1
但在做逆DFT时,因为逆DFT是插值法,求出来的形式一定与多项式的一般形式相符合而不会出现负指数,而插值正好是2n次单位复数根的k(k=0~2n-1)次方,把负号移到它们头上,所求出来的项的次数是在 mod 2n 意义下的,也就是这些负指数项正好加上2n(我们做FFT时默认A,B,C次数为2n,不过它们实际最高次数并不一定是2n)
上面可能说的很模糊,举例来说,我们在 C(x) 的原式中所谓的第 j cjxj ,因为 x=ω02n,ω12n..ω2n12n ,所以这里的 cjxj=cjx2nj ,显然 2nj>0
说白了,我们对 A(x)B(x1) 的2n次单位复数根处插值时,得到的多项式实际是
C(x)=c0+c1x+c2x2+...+cn1xn1+c(n1)xn+1+...+c2x2n2+c1x2n1
不得不说,这样做的基础就是n次单位复数根在乘法下构成n阶群以及插值法确定多项式的唯一性
(上面这句话是我瞎扯淡的,请无视)
代码:

#include<cstdio>
#include<complex>
using namespace std;
int n;
int pos[1<<18],bit[1<<18];
const double pi=acos(-1.0);
typedef complex<double> Co;
Co a[1<<18],b[1<<18];
void FFT(Co a[],bool tp,int len)
{
    for (int i=0;i<len;++i)
        if (pos[i]>i) swap(a[pos[i]],a[i]);
    for (int i=2;i<=len;i=i<<1)
    {
        Co w(cos(2*pi/i),(tp?-1:1)*sin(2*pi/i)),wn,p,q;
        for (int j=0;j<len;j+=i)
        {
            wn=Co(1,0);
            for (int k=0;k<i/2;++k)
                p=a[k+j],
                q=wn*a[k+j+i/2],
                a[k+j]=p+q,
                a[k+j+i/2]=p-q,
                wn=w*wn;
        }
    }
}
main()
{
    scanf("%d",&n);
    for (int x,y,i=0;i<n;++i)
        scanf("%d%d",&x,&y),
        a[i]=Co(x,0),b[i]=Co(y,0); 
    int wei=log2(n*2)+1,len=1<<wei;
    for (int t=1,i=0;i<wei;++i,t+=t) bit[t]=i;
    for (int i=1;i<len;++i) pos[i]=pos[i-(i&-i)]|(1<<wei-bit[i&-i]-1);
    FFT(a,0,len);
    FFT(b,1,len);
    for (int i=0;i<len;++i) a[i]=a[i]*b[i];
    FFT(a,1,len);
    for (int i=0;i<n;++i) printf("%d\n",(int)(a[i].real()/len+0.5));
}
  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值