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−1​cn−2​…c0​)−2​​, where x=∑i=0n−1ci−2ix = \sum_{i = 0}^{n - 1}{c_i \sqrt{-2}^i}x=∑i=0n−1​ci​−2​i, 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​=−2​0+−2​+−2​2=(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​=−2​2+−2​4+−2​5+−2​7=(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;


 

Logo

助力广东及东莞地区开发者,代码托管、在线学习与竞赛、技术交流与分享、资源共享、职业发展,成为松山湖开发者首选的工作与学习平台

更多推荐