zoukankan      html  css  js  c++  java
  • 任意模数FFT(拆系数FFT)

    作用

    求多项式卷积对 (p) 取模的结果,(p) 不一定能写成 (2^k+1) 的形式。

    做法

    菜鸡不会三模数只好拆系数

    将两个多项式拆开乘就不会乘爆,最后加一下取个模就行了。

    把多项式 ((A 和 B)) 的系数拆成 (a_{i}=A_{i}>>15 , b_{i}=A_{i}& 32767 ,c_{i}=B_{i}>>15 , d_{i}=B_{i}& 32767).

    然后分别乘一下, (F=a*c , G=a*b+c*d , H=b*d).

    最后答案就是 : (F<<30+G<<15+H) 对于每一项先取模再右移,可以保证不爆。

    板子题

    #include<cmath>
    #include<cstdio>
    #include<cstring>
    #include<iostream>
    #include<algorithm>
    #define N 300005
    #define LL long long 
    #define LB long double 
    using namespace std;
    
    int mod;
    int n,m,bit,tot,rev[N];
    const LB pi=3.1415926535897932384626433832795028841;
    struct complex
    {
    	LB x,y;
    	complex operator+ (const complex &t) const 
    	{
    		return (complex){x+t.x,y+t.y};
    	}
    	complex operator- (const complex &t) const
    	{
    		return (complex){x-t.x,y-t.y};
    	}
    	complex operator* (const complex &t) const 
    	{
    		return (complex){x*t.x-y*t.y,x*t.y+y*t.x};
    	}
    }a[N],b[N],c[N],d[N],f[N],g[N],h[N];
    
    inline int qr()
    {
    	char a=0;int w=1,x=0;
    	while(a<'0'||a>'9'){if(a=='-')w=-1;a=getchar();}
    	while(a<='9'&&a>='0'){x=(x<<3)+(x<<1)+(a%48);a=getchar();}
    	return x*w;
    }
    
    void init_rev(int len)
    {
    	while((1<<bit)<len)
    		bit++;
    	tot=(1<<bit);
    	for(register int i=0;i<tot;i++)
    		rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
    }
    
    inline void FFT(complex *a,int opt)
    {
    	for(register int i=0;i<tot;i++)
    		if(i<rev[i])
    			swap(a[i],a[rev[i]]);
    	for(register int mid=1;mid<tot;mid<<=1)
    	{
    		complex w1=(complex){cos(pi/(LB)mid),(LB)opt*sin(pi/(LB)mid)};
    		for(register int i=0;i<tot;i+=mid<<1)
    		{
    			complex wk=(complex){1,0};
    			for(register int j=0;j<mid;j++)
    			{
    				complex x=a[i+j];
    				complex y=wk*a[i+j+mid];
    				a[i+j]=x+y;
    				a[i+j+mid]=x-y;
    				wk=wk*w1;
    			}
    		}
    	}
    	if(opt==-1)
    	{
    		LB op=1.00000/(LB)tot;
    		for(register int i=0;i<tot;i++)
    			a[i].x=(LL)(a[i].x*op+0.5);
    	}
    }
    
    int main()
    {
    	n=qr();
    	m=qr();
    	mod=qr();
    	init_rev(m+n+1);
    	for(register int i=0;i<=n;i++)
    	{
    		int x=qr();
    		a[i].x=x>>15;
    		b[i].x=x&32767;
    	}
    	for(register int i=0;i<=m;i++)
    	{
    		int x=qr();
    		c[i].x=x>>15;
    		d[i].x=x&32767;
    	}
    	FFT(a,1);
    	FFT(b,1);
    	FFT(c,1);
    	FFT(d,1);
    	for(register int i=0;i<tot;i++)
    	{
    		f[i]=a[i]*c[i];
    		g[i]=a[i]*d[i]+b[i]*c[i];
    		h[i]=b[i]*d[i];
    	}
    	FFT(f,-1);
    	FFT(g,-1);
    	FFT(h,-1);
    	for(register int i=0;i<n+m+1;i++)
    	{
    		LL op1=((LL)(f[i].x)%mod)<<30ll;
    		LL op2=((LL)(g[i].x)%mod)<<15ll;
    		LL op3=(LL)(h[i].x)%mod;
    		printf("%lld ",(op1+op2+op3)%mod);
    	}
    	printf("
    ");
    	return 0;
    }
    
  • 相关阅读:
    sprint2第五天任务完成情况
    sprint2第四天任务完成情况
    sprint2第三天任务完成情况
    spark编程基础1
    git基本命令
    自定义bean对象实现序列化接口(Writable)
    HDFS 2.X新特性
    win10-idea连接hdfs集群
    centos6-yum源失效问题
    hadoop-源码编译
  • 原文地址:https://www.cnblogs.com/isonder/p/14412627.html
Copyright © 2011-2022 走看看