任意模数NTT（拆系数FFT）

原创已于 2023-07-10 21:17:55 修改 · 470 阅读

2 ·

CC 4.0 BY-SA版权

文章标签：

#算法

于 2023-07-10 19:39:24 首次发布

算法同时被 2 个专栏收录

44 篇文章

订阅专栏

信竞数学

10 篇文章

订阅专栏

文章介绍了在不满足特定模数条件时如何使用拆系数FFT进行任意模数的快速傅里叶变换(FFT)以实现多项式乘法。这种方法将系数拆分为高位和低位，然后通过多次FFT运算和组合来完成计算，具有较好的理解和实现性，时间复杂度为O(nlogn)。文章还提供了一个洛谷P4245题目作为示例，展示了如何应用该算法。

前置知识： $FFT$

介绍

一般的NTT的模数 $P$ 都要满足 $P=r\times 2^k+1$ 。但如果不满足这个条件，那就不能直接用NTT了，需要用到任意模数NTT。

可以了解一下三模数NTT。

下面就来介绍一种任意模数NTT——拆系数FFT。

拆系数FFT

对于多项式 $f (x)$ 的系数 $A$ ，可以将 $A$ 分为高位 $A_1$ 和低位 $A_0$ ，比如令 $A=A_1\times 2^{15}+A_0$ 。那么，可以将 $A_1$ 作为系数构造多项式 $f_1(x)$ ，将 $A_0$ 作为系数构造多项式 $f_0(x)$ 。同理，将多项式 $g (x)$ 的系数 $B$ 分为高位 $B_1$ 和低位 $B_0$ ，得到 $g_1(x)$ 和 $g_0(x)$ 。

那么， $f(x)g(x)=f_1(x)g_1(x)\times 2^{30}+(f_1(x)g_0(x)+f_0(x)g_1(x))\times 2^{15}+f_0(x)g_0(x)$

用 $FFT$ 来处理，正变换需要四次。本来逆变换也需要四次，但因为 $f_1(x)g_0(x)$ 和 $f_0(x)g_1(x)$ 是同阶的，可以相加之后在进行逆变换，所以逆变换只需要三次。

在得出三个部分之后，在乘各个部分的阶数的时候，

整个过程中，总共有 $7$ 次 $D FT$ 或 $I D FT$ 。

时间复杂度为 $O(n\log n)$ 。

相对于三模数NTT说，拆系数FFT好理解、好实现、跑得快，但有时候精度可能会出问题，而且如果结果多项式中的系数较大，long double存不下的话，也不能用，在做题时要有所注意。

例题

洛谷P4245 【模板】任意模数多项式乘法

#include<bits/stdc++.h>
using namespace std;
const int N=500000;
const long double pi=acos(-1.0);
int n,m,l=1;
long long p,f[N+5],g[N+5],w1[N+5],w2[N+5],w3[N+5],ans[N+5];
struct cp{
	long double a,b;
	cp operator +(const cp ax)const{
		return (cp){a+ax.a,b+ax.b};
	}
	cp operator -(const cp ax)const{
		return (cp){a-ax.a,b-ax.b};
	}
	cp operator *(const cp ax)const{
		return (cp){a*ax.a-b*ax.b,b*ax.a+a*ax.b};
	}
}w,wn,v1[N+5],v2[N+5],v3[N+5],v4[N+5],h1[N+5],h2[N+5],h3[N+5];
void ch(cp *a){
	for(int i=1,j=l/2,k;i<l-1;i++){
		if(i<j) swap(a[i],a[j]);
		k=l/2;
		while(j>=k){
			j-=k;k>>=1;
		}
		j+=k;
	}
}
void fft(cp *a,int fl){
	ch(a);
	for(int i=2;i<=l;i<<=1){
		wn=(cp){cos(fl*2*pi/i),sin(fl*2*pi/i)};
		for(int j=0;j<l;j+=i){
			w=(cp){1,0};
			for(int k=j;k<j+i/2;k++,w=w*wn){
				cp t=a[k],u=w*a[k+i/2];
				a[k]=t+u;
				a[k+i/2]=t-u;
			}
		}
	}
}
int main()
{
	scanf("%d%d%lld",&n,&m,&p);
	for(int i=0;i<=n;i++){
		scanf("%lld",&f[i]);f[i]%=p;
	}
	for(int i=0;i<=m;i++){
		scanf("%lld",&g[i]);g[i]%=p;
	}
	while(l<n+m+1) l<<=1;
	for(int i=0;i<l;i++){
		v1[i].a=f[i]>>15;
		v2[i].a=f[i]&((1<<15)-1);
		v3[i].a=g[i]>>15;
		v4[i].a=g[i]&((1<<15)-1);
	}
	fft(v1,1);fft(v2,1);fft(v3,1);fft(v4,1);
	for(int i=0;i<l;i++){
		h1[i]=v1[i]*v3[i];
		h2[i]=v1[i]*v4[i]+v2[i]*v3[i];
		h3[i]=v2[i]*v4[i];
	}
	fft(h1,-1);fft(h2,-1);fft(h3,-1);
	for(int i=0;i<l;i++){
		w1[i]=(long long)(h1[i].a/l+0.5)%p;
		w2[i]=(long long)(h2[i].a/l+0.5)%p;
		w3[i]=(long long)(h3[i].a/l+0.5)%p;
		ans[i]=(w1[i]*(1<<30)%p+w2[i]*(1<<15)%p+w3[i])%p;
	}
	for(int i=0;i<=n+m;i++){
		printf("%lld ",ans[i]);
	}
	return 0;
}