快速傅里叶变换(fft)————数学与算法的魔术
nk8J
链接:牛客网-该题目为付费比赛题目,请购买后查看
来源:牛客网
In this problem, you are asked to solve an easy task -- AAA multiplies BBB.
Let U={a+b−2∣a,b∈Z}U = \{a + b \sqrt{-2} | a, b \in \mathbb{Z}\}U={a+b−2∣a,b∈Z} be the whole set where multiplications will be taken. Obviously, the product of any two elements in UUU is also an element in UUU.
Each element x∈Ux \in Ux∈U can be represented in base −2\sqrt{-2}−2 as a 01-string of length nnn in the form of (cn−1cn−2…c0)−2(c_{n-1} c_{n-2} \ldots c_0)_{\sqrt{-2}}(cn−1cn−2…c0)−2, where x=∑i=0n−1ci−2ix = \sum_{i = 0}^{n - 1}{c_i \sqrt{-2}^i}x=∑i=0n−1ci−2i, ci∈{0,1}c_i \in \{0, 1\}ci∈{0,1} (i=0,1,…,n−1i = 0, 1, \ldots, n - 1i=0,1,…,n−1) and when n>1n > 1n>1, cn−1=1c_{n - 1} = 1cn−1=1. For example, −1+−2=−20+−2+−22=(111)−2-1 + \sqrt{-2} = \sqrt{-2}^0 + \sqrt{-2} + \sqrt{-2}^2 = (111)_{\sqrt{-2}}−1+−2=−20+−2+−22=(111)−2 and 2−4−2=(4−2)+(−8+4)−2=−22+−24+−25+−27=(10110100)−22 - 4 \sqrt{-2} = (4 - 2) + (-8 + 4) \sqrt{-2} = \sqrt{-2}^2 + \sqrt{-2}^4 + \sqrt{-2}^5 + \sqrt{-2}^7 = (10110100)_{\sqrt{-2}}2−4−2=(4−2)+(−8+4)−2=−22+−24+−25+−27=(10110100)−2.
Your task is to calculate the product of two elements AAA and BBB of UUU in base −2\sqrt{-2}−2 representation.
输入描述:
The first line contains an integer TTT (1≤T≤1051 \leq T \leq 10^51≤T≤105), indicating the number of test cases. Then follow TTT test cases. For each test case: The only line contains two 01-strings AAA and BBB (∣A∣,∣B∣≥1,∣A∣+∣B∣≤2×105|A|, |B| \geq 1, |A| + |B| \leq 2 \times 10^5∣A∣,∣B∣≥1,∣A∣+∣B∣≤2×105). It is guaranteed that the sum of (∣A∣+∣B∣)(|A| + |B|)(∣A∣+∣B∣) over TTT test cases does not exceed 2×1062 \times 10^62×106.
输出描述:
For each test case, output a 01-string in one line, representing the product of AAA and BBB in base −2\sqrt{-2}−2.
示例1
输入
复制5 0 1 10 10 10 11 101 101 110 110
5 0 1 10 10 10 11 101 101 110 110
输出
复制0 100 110 1 10110100
0 100 110 1 10110100
说明
For the last test case of the example, (−2+−2)2=2−4−2(-2 + \sqrt{-2})^2 = 2 - 4 \sqrt{-2}(−2+−2)2=2−4−2.
觉得有点复杂?这其实是一道fft的裸题,要优美的解决这个问题,首先我们要知道--fft是什么?
我们先来看一个问题:给定两个形如f(x)=Ax^2+Bx+C的多项式,求他们的乘积;
在数学上,我们可以直接通过乘法分配律解决,而在计算机里,我们尝试把它映射到两个list里面,这就叫多项式的系数表示法,两个d阶的多项式相乘能够得到2d阶的多项式,用分配律解决的时间复杂度在O(n^2);
现在,让我们尝试改进这个算法
对于一个一次函数f(x)=kx+b;除了已知k,b之外,还有什么能够确定这条直线呢---两点确定一条直线
那么高阶多项式呢?d阶多项式能否由d+1个点确定?
答案是肯定的,我们不妨尝试证明一下:
对于f(x)=p0+p1*x+p2*x^2.....
我们将d+1个点代入,得到d+1个方程
p(x0)=p0+p1x0+p2x0^2+......
p(x1)=p0+p1x1+p2x1^2+.......
把它写成矩阵-向量形式就是
|p(x0)| |1 x0 x0^2......x0^d| |p0|
|p(x1)|=|1 x1 x1^2......x1^d|*|p1|
|p(x2)| |1 x2 x2^2......x2^d| |p2|
deg1 2 3
只要x0~xd不重复,deg2可逆,此为范德蒙德行列式,同时p(x0)....唯一,所以多项式系数唯一
(这里用到了克拉默法则和范德蒙德行列式的性质,学过线性代数的读者应该感到熟悉)
所以得证,故我们尝试使用值表示法替代系数表示法;
对于C(x)=A(x)*B(x),我们假设A,B均为2阶多项式,则C为4阶多项式,我们需要在A,B中各找5个点相乘得到C的5个点
发现了吗,运用多项式乘法,O(d^2)的时间复杂度被我们降低到了O(d)!
所以问题变成了
1.计算两个多项式在2*d+1上的值
2.将它们相乘
3.将值重新变成系数表示
让我们一步一步来;
首先,如何把系数表示转化成值表示呢?
给定d阶多项式f(x)=p0+p1*x+p2*x^2....和n个点,计算f在这些点上的值,n>=d+1;
最简单的方式是,直接找n个点带入,但。。。
每个点的计算都是O(d)的,这么算的时间复杂度是O(nd)>O(d*d)!
这还不如不优化?
这时候你可能已经想到了,怎么能够用更少的点得出更多的值呢?或者说,在什么情况下知道了一个点的值,就相当于知道了另一个点的值呢?
偶/奇函数!
那么原函数里面有偶函数吗,我们将偶数和奇数次幂的项整合得到
f(x)=(p0+p2*x^2+...)+x*(p1+p3*x^2+...)
我们把前面的函数命名为pe(even),后面的命名为po(odd),同时用x^2代替x,就是
f(x)=pe(x^2)+x*po(x^2);
这样,对于一对相反数,我们只需计算其中一个的值就行,同时,如果我们将两个部分看成x^2的函数而不是x的,它的幂次就被降到了2
我们尝试把它推广,对于n阶多项式,我们将它拆成奇偶两部分的n/2-1次,而对于这两部分,这又是一个和开始一样的求值问题,只不过所求点变成了开始的平方;
正好,我们开始取的相反数,平方后它们的个数吧也变成了n/2-1个!
熟悉算法竞赛的读者想必已经注意到了,这不就是递归吗!
我们只需要递归的求pe和po的值,然后回代,就可以得到原多项式的n个值。
如果这一切都能走通,那么我们的时间复杂度仅为0(nlogn),因为两个子问题都是原问题规模的一半,而且我们只需要线性时间来计算原多项式的值。
但是等等,聪明的你可能已经发现了问题,在每一次递归中,我们都假设了相反数对,而x^2的值一定是正数-------推论不成立?
那么我们能不能把每次递归用的数设置成相反数呢?
正在阅读这段文字的你可能已经想到了,有且只有一种可能--把它们都设置成复数!
这也是fft最为天才之处/
同时我们需要专门寻找一些复数,使得它们平方后仍然成双成对的出现。
那么我们该怎么取呢?这可能有点难
举个例子,对于P(x)=x^3+x^2-x+1,我们需要4个点来求值,它们是一对一对的相反数,我们将它们设置成(x1,-x1,x2,-x2),它们平方后为(x1^2,x2^2),x1^2=-x2^2,再次平方得到x1^4,现在我们可以随意取它了,我们不妨设它是1,则很容易想到原来的4个数分别是(1,-1,i,-i)
现在,我们把问题推广到d阶,我们取n>d且为2的幂次,这些点就是1的n个n次方根
我们来尝试解释这个推论,
我们知道,1的n次方根可以被解释为复平面上沿着单位圆等距排布的一系列点,任一相邻两点的夹角为2pi/n,从而这些点的最简便的表示方法就是用欧拉公式写作复指数。
我们定义w=exp(2pi*i/n),那么不同的取值实际上就代表了不同的单位复n次方根
w^0=1;
w^1=e^(2pii/n)
w^2=e^(2pi*2i)/n
所以当我们说我们想在1的n次方根求值时,我们实际上是在1,w,...,w(n-1)上求值
那么不同的取值实际上就代表了不同的单位复n次方根;
所以回到原来的问题:为啥n次方根就可以用于递归过程呢?有两个主要原因
1.我们要求所取的点是正负成对的, 这里w^(i+n/2)=-W正好是这么一对;
2.当我们平方,并且考虑对氏阶奇/偶函数时,我们得到的正好是n/2个n/2 次方根,它们也是正负配对的,完美适合递归;我们再平方,也可以得到n/4个n/4次方根,适于下一层的递归,最终直到我们只剩一个点
于是,我们终于可以开始认识fft了:
FFT算法的输入是多项式的n个系数p0,...,pn-1,其中n是2的幂次(当然不是幂次也行,这里以幂次为例),我们取W=en(2pi i/n)为1的n次方根,递归的底层 情况是n=1,此时我们只在一个点求多项式的值,而递归的主要命令,就是把多项式拆成奇/偶函数部分,两部分即调用函数两次,此时这两个子多项式函数都是n/2阶的,所以对应的求值点将是1的n/2次方根,我们假设递归调用.得到的两部分函数值为ye和yo,我们把这两部分函数值结合一下,计算出原多项式函数在对应的点的函数值,我们之前看到了,结合起来的核心想法就是利用正负对,不过这里我们要为n次方根 稍微改一改我们的表达式,我们知道xj此时应该写为wN, 同时,我们之前也观察到,这里的正负点对是-w^jiw^i+n/2),更棒的是,ye和yo上第个元素,正好对应了Pe(w(2j))和Polw^(2))
==》
p(w^j)=ye[j]+w^j*yo[j];
p(w^(j+n/2))=ye[j]-w^j*yo[j];
j 1-n/2-1;
以下是伪代码的实现过程:
定义函数fft(P);
n=len(p);
//递归的底层,当n==1时输出p;
w=e^(2pi*i/n);
//将多项式拆解成奇偶两个部分;
//pe,po;
//对这两部分分别调用fft;
ye=fft(pe);
yo=fft(po);
//将这两部分拼接起来;
//初始化最后的输出list为0;
memset(y,0,sizeof y);
//根据刚才的公式计算
for(j,0,n/2-1){
y[j]=ye[j]+w^j*yo[j];
y[j+n/2]=ye[j]-w^j*yo[j];
}
//最后返回y;
return y;
现在,让我们回头看看开始的那个问题,对于两个多项式相乘,我们现在可以用FFT高效地把系数表示转换成值表示,现在只剩下反过来了:把值表示转换为系数表示。我们一般称这个过程为插值,实际上求值和插值是紧密联系的。回忆我们之前可以把求值表示为矩阵-向量乘积,我们用系数向量乘以点对应的(范德蒙德)矩阵,从而得到对应点的值。因为FFT中xk是对应的第k个n次方根,n-1 所以我们可以吧这个矩阵乘法改写
|p(w^0)| |1 1 1 | |p0|
|p(w^1)|=|1 w w^2 | * |p1|
|p(w^2)| |1 w^2 w^4| |p2|
这个矩阵称为离散傅立叶变换(DFT) 矩阵
实际上FFT一般被描述为快速计算DFT和向量乘法的算法
我们刚才说的在n次方根处求值,实际上就是DFT矩阵和向量乘法的特例,也正因此我们可以用FFT,利用这种解释,插值也变得非常易于理解;
==》
不就是把这个DFT矩阵求了个逆么!
插值实际上就是给了给定位置的函数 值,计算函数的系数,写作矩阵就是DFT矩阵的逆和值表示向量相乘;
实际上我们只需要把原来的W换成w^(-1) 即可
这表明我们几乎可以在插值问题上套用FFT
逆FFT算法的输入是多项式在n次单位根处的函数值,输出多项式的系数,即反向操作FFT算法,我们前面已经看到了,我们需要做一个DFT矩阵的逆和向量的乘法,求DFT矩阵的逆,我们只需把原来矩阵中的w换成1/n*w^(-1)
妙就妙在,这意味着逆FFT实际上和FFT 的操作一模一样,只是右边的向量变成
了值表示向量,而且w换成了1/ n*w(-1)
就这样,我们已经得到了可以计算插值的逆FFT算法
只需要在
w=e^(2pi*i/n);
处,修改
w=e^(-2pi*i/n)
然后在外面×1/n就可以了!
让我们回到牛客的L,我们可以很轻松的敲出如下代码:
#include<bits/stdc++.h>
using namespace std;
typedef double db;
typedef long long ll;
//
int n,m;
string a1,b1;
const db pi=acos(-1.0);
const int N=100;
void solve(){
cin>>a1>>b1;
n=a1.size();
m=b1.size();
vector<db>a(n),b(m);
for(int i=0;i<n;i++){
a[i]=(a1[n-i-1]=='1')?1.0:0.0;
}
for(int i=0;i<m;i++){
b[i]=(b1[m-i-1]=='1')?1.0:0.0;
}
int tot=n+m+1;
int len=1;
while(len<tot){
len<<=1;
}
vector<complex<db>>fa(len,0.0);
vector<complex<db>>fb(len,0.0);
vector<ll>c(len+100,0);
int r=-1;
for(int i=0;i<n;i++){
fa[i]=complex<db>(a[i],0);
}
for(int i=0;i<m;i++){
fb[i]=complex<db>(b[i],0);
}
auto fft=[&](auto &&self,vector<complex<db>>&a,int o)->void{
if(a.size()==0) return ;
int n=a.size();
vector<int>rev(n);
int bt=0;
while((1<<bt)<n)bt++;
for(int i=0;i<n;i++){
rev[i]=(rev[i>>1]>>1|((i&1)<<(bt-1)));
if(i<rev[i]) swap(a[i],a[rev[i]]);
}
for(int i=1;i<n;i<<=1){
complex<db>w(cos(pi/i),-o*sin(pi/i));
for(int j=0;j<n;j+=i<<1){
complex<db>ww(1,0);
for(int k=0;k<i;k++,ww*=w){
complex<db>x=a[k+j],y=a[i+j+k]*ww;
a[j+k]=x+y;
a[i+j+k]=x-y;
}
}
}
if(o==-1){
for(int i=0;i<n;i++){
a[i]/=n;
}
}
};
fft(fft,fa,1);
fft(fft,fb,1);
for(int i=0;i<len;i++){
fa[i]=fa[i]*fb[i];
}
fft(fft,fa,-1);
for(int i=0;i<tot;i++){
c[i]=llround(fa[i].real());
}
for(int k=0;k<len+N;k++){
if(c[k]==0)continue;
ll val=c[k];
ll rtf=val%2;
if(rtf<0){
rtf+=2;
}
ll tt=(val-rtf)/2;
c[k]=rtf;
c[k+2]-=tt;
}
for(int i=len+N-1;i>=0;i--){
if(c[i]!=0){
r=i;
break;
}
}
if(r==-1)cout<<0<<"\n";
else{
for(int i=r;i>=0;i--){
cout<<char(c[i]+'0');
}
cout<<"\n";
}
}
int main(){
ios::sync_with_stdio(false);
cin.tie(0);
int t;
cin>>t;
while(t--){
solve();
}
}
MINNGYAAA,8.10 23:56;
更多推荐



所有评论(0)