首页 > 其他分享 >常数比较小码量不大的 MTT(4次FFT)/任意模数多项式乘法

常数比较小码量不大的 MTT(4次FFT)/任意模数多项式乘法

时间:2023-01-05 09:23:10浏览次数:59  
标签:15 小码量 Comp FFT MTT int len reg

先根据 Prean 的题解 写出一个常数较小的 5 次 FFT 写法。

inline ll get(const double x){return (ll(x+0.5))%mod;}
inline void MTT(const int *A,const int *B,int *C,int n,int m){
	static Comp F[maxn],G[maxn],T[maxn];len=1;while(len<=n+m)len<<=1;init(len);
	memset(F+n+1,0,(len-n-1)*sizeof(Comp)),memset(G+n+1,0,(len-n-1)*sizeof(Comp));
	memset(T+m+1,0,(len-m-1)*sizeof(Comp));
	for(reg i=0;i<=n;++i)F[i]=Comp(A[i]&32767,0),G[i]=Comp(A[i]>>15,0);
	for(reg i=0;i<=m;++i)T[i]=Comp(B[i]&32767,B[i]>>15);
	FFT(F,len),FFT(G,len),FFT(T,len);
	for(reg i=0;i<len;++i)F[i]*=T[i];
	for(reg i=0;i<len;++i)G[i]*=T[i];
	IFFT(F,len),IFFT(G,len);
	for(reg i=0;i<len;++i)C[i]=(get(F[i].real)+(get(F[i].imag+G[i].real)<<15)+(get(G[i].imag)<<30))%mod;
}


又发现

于是根据 Kewth 的题解,写出了 4 次 FFT 版本。

inline void FFT2(Comp *a,Comp *b,int n){
	for(reg i=0;i<n;++i)a[i].imag=b[i].real;
	FFT(a,n);
	for(reg i=0;i<n;++i)b[i]=conj(a[i?n-i:0]);
	for(reg i=0;i<n;++i){
		Comp p=a[i],q=b[i];
		a[i]=(p+q)*0.5,b[i]=(q-p)*0.5*Comp(0,1);
	}
}
int mod;
inline ll get(const double x){return (ll(x+0.5))%mod;}
inline void MTT(const int *A,const int *B,int *C,int n,int m){
	static Comp F[maxn],G[maxn],T[maxn];len=1;while(len<=n+m)len<<=1;init(len);
	memset(F+n+1,0,(len-n-1)*sizeof(Comp)),memset(G+n+1,0,(len-n-1)*sizeof(Comp));
	memset(T+m+1,0,(len-m-1)*sizeof(Comp));
	for(reg i=0;i<=n;++i)F[i]=Comp(A[i]&32767,0),G[i]=Comp(A[i]>>15,0);
	for(reg i=0;i<=m;++i)T[i]=Comp(B[i]&32767,B[i]>>15);
	FFT2(F,G,len),FFT(T,len);
	for(reg i=0;i<len;++i)F[i]*=T[i];
	for(reg i=0;i<len;++i)G[i]*=T[i];
	IFFT(F,len),IFFT(G,len);
	for(reg i=0;i<len;++i)C[i]=(get(F[i].real)+(get(F[i].imag+G[i].real)<<15)+(get(G[i].imag)<<30))%mod;
}

优化一下常数,最终卡进了最优第二页(最底端)

#include<bits/stdc++.h>
#define EL puts("Elaina")
#define reg register int
typedef long long ll;
using namespace std;
namespace IO{
	const int siz=1<<18;char buf[siz],*p1,*p2;
	inline char getc(){return p1==p2&&(p2=(p1=buf)+fread(buf,1,siz,stdin),p1==p2)?EOF:*p1++;}
	inline int read(){
		int x=0;char ch=getc();
		while(!isdigit(ch))ch=getc();
		while(isdigit(ch))x=x*10+(ch^48),ch=getc();
		return x;
	}
}using IO::read;
const int maxn=4e5+3;const double Pi=acos(-1);
struct Comp{
	double real,imag;
	Comp(){}
	Comp(double a,double b){real=a,imag=b;}
	inline Comp operator +(const Comp &a)const{return Comp(real+a.real,imag+a.imag);}
	inline Comp operator -(const Comp &a)const{return Comp(real-a.real,imag-a.imag);}
	inline Comp operator *(const Comp &a)const{return Comp(real*a.real-imag*a.imag,imag*a.real+real*a.imag);}
	inline Comp operator *(const double &a)const{return Comp(real*a,imag*a);}
	inline void operator+=(const Comp &a){real+=a.real,imag+=a.imag;}
	inline void operator*=(const Comp &a){(*this)=(*this)*a;}
}w[maxn];
int rev[maxn],len;
inline void init(int len){
	int mid=len>>1;w[mid]=Comp(1,0);
	for(reg i=0;i<len;++i){rev[i]=rev[i>>1]>>1;if(i&1)rev[i]|=mid;}
	for(reg i=1;i<mid;++i)w[i+mid]=Comp(cos(i*2*Pi/len),sin(i*2*Pi/len));
	for(reg i=mid-1;i>0;--i)w[i]=w[i<<1];
}
inline void FFT(Comp *a,int n){
	Comp t;for(reg i=0;i<n;++i)if(i<rev[i])swap(a[i],a[rev[i]]);
	for(reg u=2,v=1;u<=n;v=u,u<<=1)for(reg L=0;L<n;L+=u)
		for(reg p=L,x=v;p<L+v;++p,++x)t=w[x]*a[p+v],a[p+v]=a[p]-t,a[p]+=t;
}
inline void IFFT(Comp *a,int n){
	FFT(a,n),reverse(a+1,a+n);for(reg i=0;i<n;++i)a[i].real/=n,a[i].imag/=n;
}
int mod;
inline ll get(const double x){return (ll(x+0.5))%mod;}
inline void MTT(const int *A,const int *B,int *C,int n,int m){
	static Comp F[maxn],G[maxn],T[maxn];len=1;while(len<=n+m)len<<=1;init(len);
	memset(F+n+1,0,(len-n-1)*sizeof(Comp)),memset(T+m+1,0,(len-m-1)*sizeof(Comp));
	for(reg i=0;i<=n;++i)F[i]=Comp(A[i]&32767,A[i]>>15);
	for(reg i=0;i<=m;++i)T[i]=Comp(B[i]&32767,B[i]>>15);
	FFT(F,len),FFT(T,len),G[0]=F[0];
	for(reg i=1;i<len;++i)G[i]=F[len-i];
	for(reg i=0;i<len;++i){
		Comp p=F[i],q=G[i];
		F[i]=Comp(p.real+q.real,p.imag-q.imag)*T[i]*0.5;
		G[i]=Comp(p.imag+q.imag,q.real-p.real)*T[i]*0.5;
	}
	IFFT(F,len),IFFT(G,len);
	for(reg i=0;i<len;++i)C[i]=(get(F[i].real)+(get(F[i].imag+G[i].real)<<15)+(get(G[i].imag)<<30))%mod;
}
int n,m,F[maxn],G[maxn];
inline void MyDearMoments(){
	n=read(),m=read(),mod=read();
	for(reg i=0;i<=n;++i)F[i]=read();
	for(reg i=0;i<=m;++i)G[i]=read();
	MTT(F,G,F,n,m);
	for(reg i=0;i<=n+m;++i)printf("%d ",F[i]);puts("");
}
int main(){return MyDearMoments(),0;}

标签:15,小码量,Comp,FFT,MTT,int,len,reg
From: https://www.cnblogs.com/muel-imj/p/17026555.html

相关文章

  • 通过二维FFT变换对比加入窗函数之后的图像频谱和相位
    目录一、理论基础1.1二维FFT变换1.2窗函数二、核心程序三、测试结果一、理论基础1.1二维FFT变换以下公式定义m×n矩阵X的离散傅里叶变换 Y:    i......
  • 快速傅里叶变换(FFT)(下)
    关于学习FFT算法的资料个人最推荐的还是算法导论上的第30章(第三版),多项式与快速傅里叶变换,基础知识都讲得很全面。FFT算法基本概念:FFT(FastFourierTransformation)即......
  • FFT
    最近刚学了FFT,给我这个蒟蒻的脑细胞干废了,写篇笔记浅浅记录一下先讲一下卷积,卷积就是求\(\int_{-\infty}^{\infty}f(\tau)g(x-\tau)d\tau\),这玩意的用处有一个是多项式乘......
  • 手撕fft算法--fft原理和源码解析
    一前言 在音频信号处理中,fft变换是一个无法绕过过去的存在。借着一次算法出来的机会,把fft熟悉一下不为过啊。 二问题 这里,其实是由一个问题驱动的,那......
  • FFT入门——学习笔记
    FFT入门给一个非常好的入门视频:快速傅里叶变换复数与单位根定义:\(i^2=-1\)为虚数单位,我们称形如\(a+bi(a,b\inR)\)的数为复数。我们可以用复数在复平面上表示点\((0,......
  • 选带快速傅立叶变换ZOOM-FFT的matlab实现
    up目录一、理论基础二、核心程序三、测试结果一、理论基础zoom-fft是一个信号传输过程,其中输入信号被向下混频到基带,然后被抽取,然后被传递到标准FFT。ZOOM-FFT称为细......
  • Windows7下git配置difftool
    GIT是一个代码版本控制工具,是软件开发团队中必不可少的一类工具,类似的工具还有像SVN,CVS等;在此之前我一直使用的SVN,因为SVN在windows下有很好的客户端【小乌龟】,使用起来简单......
  • 利用FFT计算非平稳随机信号的WVD分布
    up目录一、理论基础二、核心程序三、测试结果一、理论基础fft:快速傅里叶变换(fastFouriertransform),即利用计算机计算离散傅里叶变换(DFT)的高效、快速计算方法......
  • matlab_fft函数c语言实现
    前言最近工作移植PPG算法,将MATLAB上代码移植到嵌入式设备上去。因为心率算法利用FFT实现会较为简单,所以又重新了解了一下大学里学的FFT,并写了C语言实现MATLAB的FFT接口的......
  • quartus中使用FFT IP核
     一、准备工作  首先需要把需要的器材准备好,我使用的是quartus18.0,并且要使用IP核被破解的版本,不然无法使用其中的FFT和NCO,一定要注意,quartus对于版本非常敏感,一定要......