洛谷 P4245 【模板】MTT

题目描述

给定 2 个多项式 F(x),G(x) ，请求出 F(x)∗G(x) 。

系数对 p 取模，且不保证 p 可以分解成 p=a⋅2^k+1 之形式。

输入输出格式

输入格式：

输入共 3 行。
第一行 3 个整数 n,m,p ，分别表示 F(x),G(x) 的次数以及模数 p 。
第二行为 n+1 个整数， 第 i 个整数 a_i 表示 F(x) 的 i-1 次项的系数。
第三行为 m+1 个整数， 第 i 个整数 b_i 表示 G(x) 的 i-1 次项的系数。

输出格式：

输出 n+m+1 个整数， 第 i 个整数 c_i 表示 F(x) * G(x) 的 i-1 次项的系数。

输入输出样例

输入样例#1：

5 8 28
19 32 0 182 99 95
77 54 15 3 98 66 21 20 38

输出样例#1：

7 18 25 19 5 13 12 2 9 22 5 27 6 26

说明

$1 \leq n \leq 10^5, 0 \leq a_i, b_i \leq 10^9, 2 \leq p \leq 10^9 + 9,1≤n≤10^ 5$

Solution

• 第一道拆系数FFT，被神奇错误卡了很久很久……

• 把系数拆成 $kM+b$ ，这样乘就不会有精度问题了。

• 然而问题是这个M要取2的整数次幂（即 $M=2^{15}$），不然会有迷之错误。

• 有高手知道可以指导我一下……

Code

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<cctype>
using namespace std;
typedef long long LL;
const int N=1e5+5,pp=1<<15;
const double Pi=acos(-1.0);
struct comp
{
    double r,i;
    comp(){}
    comp(double rr,double ii){r=rr,i=ii;}
    friend comp operator +(comp x,comp y)
    {
        return comp(x.r+y.r,x.i+y.i);
    }
    friend comp operator -(comp x,comp y)
    {
        return comp(x.r-y.r,x.i-y.i);
    }
    friend comp operator *(comp x,comp y)
    {
        return comp(x.r*y.r-x.i*y.i,x.r*y.i+x.i*y.r);
    }
};
int n,m,p,len;
int rev[N<<2],f[N<<1];
comp a1[N<<2],b1[N<<2],a2[N<<2],b2[N<<2];
comp c1[N<<2],c2[N<<2],c3[N<<2],w[N<<2];
inline int read()
{
    int X=0,w=0; char ch=0;
    while(!isdigit(ch)) w|=ch=='-',ch=getchar();
    while(isdigit(ch)) X=(X<<1)+(X<<3)+(ch^48),ch=getchar();
    return w?-X:X;
}
inline void write(int x)
{
    if(x>9) write(x/10);
    putchar(x%10+'0');
}
inline void FFT(comp *y,int ff)
{
    for(int i=0;i<len;i++)
        if(i<rev[i]) swap(y[i],y[rev[i]]);
    for(int h=2;h<=len;h<<=1)
        for(int i=0,mid=h/2,num=len/h;i<len;i+=h)
        {
            for(int k=i,cnt=0;k<i+mid;k++,cnt++)
            {
                comp u=y[k],t=w[ff>0?num*cnt:len-num*cnt]*y[k+mid];
                y[k]=u+t;
                y[k+mid]=u-t;
            }
        }
    if(ff==-1) for(int i=0;i<len;i++) y[i].r/=len;
}
int main()
{
    n=read()+1,m=read()+1,p=read();
    for(int i=0;i<n;i++)
    {
        int x=read();
        a1[i].r=x/pp;
        b1[i].r=x%pp;
    }
    for(int i=0;i<m;i++)
    {
        int x=read();
        a2[i].r=x/pp;
        b2[i].r=x%pp;
    }
    int l=0;
    for(len=1;len<n+m;len<<=1) l++;
    for(int i=0;i<len;i++) rev[i]=rev[i>>1]>>1|(i&1)<<l-1;
    for(int i=0;i<=len;i++) w[i]=comp(cos(2*Pi*i/len),sin(2*Pi*i/len));
    FFT(a1,1),FFT(b1,1);
    FFT(a2,1),FFT(b2,1);
    for(int i=0;i<len;i++)
    {
        c1[i]=a1[i]*a2[i];
        c2[i]=a1[i]*b2[i]+a2[i]*b1[i];
        c3[i]=b1[i]*b2[i];
    }
    FFT(c1,-1),FFT(c2,-1),FFT(c3,-1);
    for(int i=0;i<n+m-1;i++)
    {
        f[i]=(LL)(c1[i].r+0.5)%p*pp%p*pp%p;
        f[i]=(f[i]+(LL)(c2[i].r+0.5)%p*pp%p)%p;
        f[i]=(f[i]+(LL)(c3[i].r+0.5)%p)%p;
    }
    for(int i=0;i<n+m-1;i++) write(f[i]),putchar(' ');
    return 0;
}