[模板] MTT实现之FFT拆系数

51 篇文章 1 订阅
12 篇文章 0 订阅

题目描述:

求两个多项式的卷积,系数对P取模,不保证P可以分解成 P=a2k+1 P = a ⋅ 2 k + 1

题目分析:

P不保证分解成那个形式,那么我们就不可以用NTT了.
如何解决这个问题呢?
拆系数+FFT~
基本上就是找个模数
然后拆成4个多项式,分别是第一个多项式 / 模数 %模数 第二个多项式 / 模数 %模数 的值
然后两两相乘,结果 乘上对应项数的模数次方然后%要取的模数
这样的话是做4次DFT,4次IDFT
然而单位复根的精度貌似会出问题
可以预处理,懒得话可以开long double

题目链接:

Luogu 4245

Ac 代码:

#include<cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <algorithm>
#include <iostream>
#define ll long long 
#define double long double
const double PI=std::acos(-1);
const int maxm=5e5+100,M=32768;
struct complex{
    double real,imag;
    complex(){};
    complex(double _real,double _imag):real(_real),imag(_imag){}
};
inline complex operator + (complex x,complex y){return (complex){x.real+y.real,x.imag+y.imag};}
inline complex operator - (complex x,complex y){return (complex){x.real-y.real,x.imag-y.imag};}
inline complex operator * (complex x,complex y){return (complex){x.real*y.real-x.imag*y.imag,x.real*y.imag+y.real*x.imag};}
int F[maxm],G[maxm],rev[maxm],l;
inline void FFT(complex *a,int n,int f)
{
    for(int i=0;i<n;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
    for(int i=0;i<n;i++) if(i<rev[i]) std::swap(a[i],a[rev[i]]);
    for(int i=1;i<n;i<<=1)
    {
        complex wn=(complex){std::cos(PI/i),f*std::sin(PI/i)};
        for(int j=0;j<n;j+=(i<<1))
        {
            complex w=(complex){1,0};
            for(int k=0;k<i;k++,w=w*wn)
            {
                complex x=a[j+k];
                complex y=a[i+j+k]*w;
                a[j+k]=x+y;
                a[i+j+k]=x-y;
            }
        }
    }
    if(f==1) return;
    for(int i=0;i<n;i++) a[i].real=a[i].real/(double)n;
}
int n,m,p;
ll ans[maxm];
complex a[maxm],b[maxm],c[maxm],d[maxm],e[maxm],f[maxm],g[maxm],h[maxm];
inline void MTT()
{
    n+=m;
    for(m=1;m<=n;m<<=1) l++;
    for(int i=0;i<n;i++) 
    {
        a[i].real=F[i]/M,b[i].real=F[i]%M;
        c[i].real=G[i]/M,d[i].real=G[i]%M;
    }
    FFT(a,m,1),FFT(b,m,1),FFT(c,m,1),FFT(d,m,1);
    for(int i=0;i<m;i++)
    {
        e[i]=a[i]*c[i],f[i]=a[i]*d[i];
        g[i]=b[i]*c[i],h[i]=b[i]*d[i];

    }
    FFT(e,m,-1),FFT(f,m,-1),FFT(g,m,-1),FFT(h,m,-1);
    for(int i=0;i<m;i++)
    {
        ans[i]=(ll)(round(e[i].real))%p*M%p*M%p;
        ans[i]+=(ll)(round(f[i].real))%p*M%p;
        ans[i]+=(ll)(round(g[i].real))%p*M%p;
        ans[i]+=(ll)(round(h[i].real))%p;
        ans[i]%=p;
    }
}
int main()
{
    scanf("%d%d%d",&n,&m,&p);
    for(int i=0;i<=n;i++) scanf("%d",&F[i]);
    for(int i=0;i<=m;i++) scanf("%d",&G[i]);
    MTT();
    for(int i=0;i<=n;i++) printf("%lld ",ans[i]);
    return 0;
}
  • 0
    点赞
  • 3
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值