zoukankan      html  css  js  c++  java
  • 任意模数FFT

    任意模数FFT

    这是一个神奇的魔法,但是和往常一样,在这之前,先

    ( exttt{orz} color{orange}{ exttt{matthew99}})

    问题描述

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

    系数对 p 取模(2 le p le 10^9+9)

    拆系数FFT

    我们考虑令(M)(sqrt{p}),那么我们可以将原本的多项式拆成4个。

    (F(x)=A(x)*M+B(x))

    (G(X)=C(X)*M+D(x))

    然后(A),(B),(C),(D)随便乘一下就可以求出答案。

    这样需要的(FFT)次数为7次。

    合并FFT

    我们思考一下有没有什么优化的方法呢?

    (P(x)=A(x)+iB(x))(Q(x)=A(x)-iB(x))

    (F_p[k])(F_q[k])分别表示(P)(Q)做了(DFT)后的(k)项系数。

    推一推式子可以发现(F_q[k]=conj(F_p[2L-k]))

    那么我们只需要1遍(DFT)就可以求出这两个东西。

    此时我们把模意义下的多项式乘法转换成了复数意义下的多项式乘法,就只需要按照(FFT)的过程实现。

    把实部与虚部合并即可。

    注意此时(FFT)的精度十分爆炸,所以用(long double)并预处理单位根比较好。

    /*
      mail: mleautomaton@foxmail.com
      author: MLEAutoMaton
      This Code is made by MLEAutoMaton
    */
    #include<stdio.h>
    #include<stdlib.h>
    #include<string.h>
    #include<math.h>
    #include<algorithm>
    #include<queue>
    #include<set>
    #include<map>
    #include<iostream>
    using namespace std;
    #define ll long long
    #define re register
    #define file(a) freopen(a".in","r",stdin);freopen(a".out","w",stdout)
    inline int gi()
    {
    	int f=1,sum=0;char ch=getchar();
    	while(ch>'9' || ch<'0'){if(ch=='-')f=-1;ch=getchar();}
    	while(ch>='0' && ch<='9'){sum=(sum<<3)+(sum<<1)+ch-'0';ch=getchar();}
    	return f*sum;
    }
    const int N=1000010;const long double Pi=acos(-1.0);
    int n,m,f[N],g[N],Mod;
    struct node
    {
    	long double x,y;
    	node operator+(const node &b)const{return (node){x+b.x,y+b.y};}
    	node operator-(const node &b)const{return (node){x-b.x,y-b.y};}
    	node operator*(const node &b)const{return (node){x*b.x-y*b.y,x*b.y+y*b.x};}
    };
    int r[N],limit,ans[N];
    void fft(node *a,int opt)
    {
    	for(int i=0;i<limit;i++)
    		if(i<r[i])swap(a[i],a[r[i]]);
    	for(int mid=1;mid<limit;mid<<=1)
    	{
    		node Root=(node){cos(Pi/mid),opt*sin(Pi/mid)};
    		for(int R=mid<<1,i=0;i<limit;i+=R)
    		{
    			node W=(node){1,0};
    			for(int j=0;j<mid;j++,W=W*Root)
    			{
    				node X=a[i+j],Y=W*a[mid+i+j];
    				a[i+j]=X+Y;a[mid+i+j]=X-Y;
    			}
    		}
    	}
    }
    node a[N],b[N],c[N],d[N];
    void mtt()
    {
    	int LIM=(1<<15)-1;
    	for(int i=0;i<=n;i++)f[i]=(f[i]+Mod)%Mod;for(int i=0;i<=m;i++)g[i]=(g[i]+Mod)%Mod;
    	for(int i=0;i<=n;i++)a[i]=(node){f[i]&LIM,f[i]>>15};
    	for(int i=0;i<=m;i++)b[i]=(node){g[i]&LIM,g[i]>>15};
    	fft(a,1);fft(b,1);
    	for(int i=0;i<limit;i++)c[i]=a[i],d[i]=b[i];
    	for(int i=0;i<limit;i++)
    	{
    		int j=(limit-i)&(limit-1);
    		a[i]=(node){0.5*(c[i].x+c[j].x),0.5*(c[i].y-c[j].y)}*d[i];
    		b[i]=(node){0.5*(c[i].y+c[j].y),0.5*(c[j].x-c[i].x)}*d[i];
    	}
    	fft(a,-1);fft(b,-1);
    	for(int i=0;i<limit;i++)
    	{
    		ll ia,ib,ic,id;
    		ia=(ll)(a[i].x/limit+0.5)%Mod;
    		ib=(ll)(a[i].y/limit+0.5)%Mod;
    		ic=(ll)(b[i].x/limit+0.5)%Mod;
    		id=(ll)(b[i].y/limit+0.5)%Mod;
    		ans[i]=((ia%Mod+((ib+ic)%Mod<<15))%Mod+(id%Mod<<30))%Mod;
    	}
    }
    int main()
    {
    #ifndef ONLINE_JUDGE
    	freopen("in.in","r",stdin);
    #endif
    	n=gi();m=gi();Mod=gi();
    	for(int i=0;i<=n;i++)f[i]=gi();for(int i=0;i<=m;i++)g[i]=gi();
    	limit=1;int l=0;
    	while(limit<=(n+m))limit<<=1,l++;
    	for(int i=0;i<limit;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    	mtt();
    	for(int i=0;i<=n+m;i++)printf("%d ",(ans[i]+Mod)%Mod);
    	return 0;
    }
    

    题目

    任意模数NTT

    Shell Necklace

  • 相关阅读:
    小白开学Asp.Net Core 《一》
    小白开学Asp.Net Core 开篇
    分享一个.NET平台开源免费跨平台的大数据分析框架.NET for Apache Spark
    微信支付退款中发现的一个问题
    发布mvc遇到的HTTP错误 403.14-Forbidden解决办法
    English--元音
    开发工具--浅谈Git
    开发工具--搭建python环境
    开发工具--PyCharm
    English--介词省略句型与总结
  • 原文地址:https://www.cnblogs.com/mleautomaton/p/11428886.html
Copyright © 2011-2022 走看看