数论变换NTT

0. 前言

我们在学Fast Fourier Transforms的时候就会发现输出栏有res[i]=(unsigned long)(a[i].real()/limit+.5)

这里需要加上0.5以保证输出精度
输出精度是怎么产生的呢?
我们用复数运算,这样复数的实部和虚部都需要使用双精度浮点数来运算。
但是浮点数的存储方式和整数不太一样,而且有时候三角函数值出来没准还是无限不循环小数。
误差就这样产生了。
再者我们也有某个神奇的结论:

整形的运算比浮点型运算快5~20倍

所以我们就想,能不能用一些神奇的算法,使得用整数也可以算多项式呢?
我们回想一下为什么多项式里面我们要用复数。
因为我们要用单位根,单位根是复数。
为什么要用单位根?因为(ωn)n=1,ωnk+n2=ωnk
满足这两个性质,我们就可以优化Fourier Transforms
我们看看怎么使用整数也可以满足这两条性质

  1. 其中第一个性质是用来解决二分时的取值问题的,(ldq大佬:也就是要满足循环卷积),我们只需要保证我们能够快速去到我们需要的值就可以了。
  2. 至于第二个性质,确实在保证第一个性质下第二个性质很难满足。但是我们退一步想,我们需要xi=xi+n2是因为我们在f(x)=f0(x2)+xf1(x2)里面使用二分,也就是要保证xi2=xi+n22。所以并不是说一定要相反数,我们只需要保证这两个数平方相等就可以了。
    至于平方相等,我们很容易想到的是12=12(废话),单位根也正是利用了这个性质:
    (ωni)2=(ωni)2ωnn=(ωni+n2)2,这里有ωnn=1,才能满足这条性质。

所以,我们就是要找一个n次方变成1的整数集合X,并且集合中的数满足xX,!yX,yx,x2=y2
这么看来在整数内幂为1的就数只有一个。。。
貌似不太对劲。
没关系,我们看看大部分的题目

答案对998244353取模

这就好办了,在取模意义下大把大把的有,更何况998244353还是一个质数
所以我们只需要找出在模数意义下的X集合就可以了


(没见过前言讲这么长的吧?)

原根

那么整数,模意义下我们知道了很多数的幂都可以是1
但是我们回过头看Fast Fourier Transforms

然后他们又二分,又二分。。。
一路二分下去,最后的时间复杂度就降到了O(nlog2n)
需要满足的条件就是当前多项式的次数必须是奇数,除非你递归到了边界:0
所以我们就要求对于所有的次数必须都是偶数,也就是一开始多项式的次数必须是2的整数次幂。

这么说我们这个幂指数n必须是2的正整数次幂
我们就要找出一个g,使得gn1(mod998244353),n=2k,kZ
那么这个g就可以代替Fast Fourier Transforms里面的n次单位根ωn
这个g就叫做原根。
在观看了这篇文章之后相信我对原根有了很深的体会。
首先数学概念里有一个东西叫做
阶就是说,有两个互质的数a,p>1,满足an1(modp)的最小正整数n就叫做ap意义下的阶,记作δp(a)=n


然后我们再来了解一个欧拉函数φ(a),表示的是小于a且与a互质的数的个数。
比如说小于10且和10互质的数有1,3,7,94个,那么φ(10)=4
至于这个欧拉函数的求法,很简单也很粗暴。
a的质因数的集合是pi,那么φ(a)=a×(11pi)
其实很好理解,不和a互质就是pi的倍数,那么小于a里面是pi的倍数的有npi个,那么剩下的n(11pi)个就不是pi的倍数,就有可能和n互质。
那么我们对于剩下的质因数也执行这样的操作,我们就可以找齐所有不是pi的倍数的数,就是和n互质的数了。


原根,就是满足δm(g)=φ(m)的正整数a叫做m的一个原根。
很明显原根g满足gφ(m)1(modm)
那么我们就可以发现,m有原根,当且仅当m=2,4,pa,2pa,aN,p为奇素数。
至于为什么。。。我也不知道。欢迎评论区补充。
然后m的原根会满足gigj(modm)(0i<j<φ(m)),很明显这么说的意思就是g0gφ(m)1(modm)
这个还可以证明,用反证法

0i<j<φ(m),gigj(modm),就会有gji1(modm),且ji0
但是i<j<φ(m)=δm(g),他们的差ji是不可能达到δm(g)的,而阶δm(g)的定义就是满足gn1(modm)最小正整数n
那么i,j,0<ji<δm(g)都不能满足gji1(modm)
故原命题成立,gigj(modm)(0i<j<φ(m))

现在设p为奇素数,那么有φ(p)=p1
那么我们根据原根的性质就会有对于p的原根g满足gi,i(0,p1)互不相同,
但是有g0gp11(modp),这个就和单位根ωn0=ωnn=1的性质很相似
那么我们就仿照一下就会有

(gi+p12)2(gi)2(gp12)2(gi)2(modp)

然后我们再发现给的模数998244353=0b111011100000000000000000000001
代入p,有p1=119223,恰好又满足了我们的n需要是2的正整数次幂(反正你也取不到223题目最多给到就是221所以是可以跑的)
现在到了最后一个问题:怎么求原根。
那很明显,我们需要满足的第一个条件就是gφ(m)1(modm)
但是只满足这个条件的不一定是原根。这只是一个必要条件。
我们当初说的是δm(g)=φ(m),这个δm(g)是最小的n使得gn1(modm)
如果我们的φ(m)不是最小的,那么他就不是δm(g),这个g就不是原根了
举个例子,我们求7的原根,有φ(7)=6,261(mod7),但是其实231(mod7)
所以δ7(2)=3φ(7)=6,2不是7的原根
那么我们怎么排除这种情况呢?这里用到了裴蜀定理

裴蜀定理是说明了对任何整数 a,b和它们的最大公约数gcd(a,b)=d ,关于未知数x,y的线性的丢番图方程(称为裴蜀等式)。
丢番图方程又是什么呢?就是不定方程,他的求解只在整数范围内进行,最后这个限制使得丢番图方程求解与实数范围方程求解有根本的不同。
裴蜀定理说明了两个命题
A:x,yZ,gcd(a,b)|(ax+by)B:x,yZ,(ax+by)=gcd(a,b)
它的一个重要推论是:a,b互质的充分必要条件是存在整数x,y使ax+by=1。
他的证明是这样的:(其实百度百科下面有):
d=gcd(a,b),则d|a,d|b,有d|ax+by(x,yZ),这很简单,命题A得证。
接下来一大段都是证明命题B
s=ax0+by0ax+by的所有取值中的最小正值,那么再令q=as,r=a%s=aqs=aq(ax0+by0)=a(1qx0)+b(qy0)
那么我们知道这个r也是a,b的线性组合,如果我们令x1=1qx0,y1=qy0,那么就有gcd(a,b))|ax1+by1,d|r
r=a%s0r<s
又因为s应该是ax+by的最小正值,那么r想要是所有ax+by中比s小的非负值,就只能a=b=0的时候r=0
r=a%s=0,s|a
同理,如果我们设q=bs,r=b%s就可以证明s|b,这里不再赘述。
所以s|a,s|b,s就是a,b的公约数,而d是他们的最大公约数,所以有ds
然后我们也知道对于这个最大公约数,我们一定有d|a,d|b,然后就有d|(ax0+by0),d|s
所以就可以推断出ds
ds,dsd=s
所以当x=x0,y=y0时,我们就有ax+by=gcd(a,b),所以x,yZ,ax+by=gcd(a,b),命题B得证。

说完了裴蜀定理我们就来看看怎么排除不是原根的情况。

首先我们满足了gφ(m)1(modm),然后我们假设k[0,φ(m)),gk1(modm),这样就使得g不是m的原根
根据裴蜀定理,x,(y)Z,kx+φ(m)(y)=gcd(k,φ(m))
kx=φ(m)y+gcd(k,φ(m))
同时我们有{gk1(modm),gkx1x(modm)gφ(m)1(modm),gφ(m)y1y(modm)
那么很简单了,我们可以推断出gkgkx(modm),代入kx=φ(m)y+gcd(k,φ(m))
gkgφ(m)y+gcd(k,φ(m))gφ(m)yggcd(k,φ(m))1ggcd(k,φ(m))(modm)
只取恒等式最左边和最右边就是gkggcd(k,φ(m))(modm),同时gk1(modm)那么ggcd(k,φ(m))1(modm)
又因为k<φ(m),所以gcd(k,φ(m))k<φ(m),这个保证了gcd(k,φ(m))φ(m)
同时我们有gcd(k,φ(m))|φ(m)
我们知道妨碍g不能成为m的原根主要就是因为有不知道哪个数k<φ(m),gk1(modm)
那么我们根据上面的等价变换(都是等价的吧上面的推理)
我们可以知道,如果存在这样一个k使得gk1(modm),那么一定存在k=gcd(k,φ(m))也能使gk1(modm)
而这个k|φ(m),所以我们只需要验证φ(m)所有的因数n看他能不能使gn1(modm),如果可以,那么这个g就不是m的原根。
但是要验证所有的因数好麻烦啊。。。
不慌我们有一个很快的方法。
如果合数n=p1p2φ(m)的因数,那么p1,p2肯定都是φ(m)的因数,只要gn1(modm),那么gp1,gp2在模m意义下肯定也不会和1同余,(反证法:)如果同余那么(gp1)p2,(gp2)p1之中的一个就是1为底的指数式,这样出来的结果在模m意义下肯定也是和1同余的,但是他们两个其实和gn是等价的,这就违背了gn1(modm),所以我们只需要把他所有的质因数都乘起来再验证这个数就可以了。
那么。。。这个数不就是m本身吗。。。
我们看看是哪里出了问题:我们验证nφ(m)的因数,以及gn1(modm),这个大前提是因为n有可能会妨碍g成为原根。那么这个n就一定满足n<φ(m),如果等于了还妨碍什么?
我们把所有质因数都乘起来那么就直接等于m了,所以不可行。那么我们可以先把所有质因数都乘起来,然后抠掉一个,再判断,这样就可以知道有没有一个数妨碍g成为原根了
如果pi,gφ(m)pi1(modm),其中pim的质因数,那么我们就可以知道这个g就是m的一个原根。

对,我说的是一个原根。
原根可以有0个或者φ(φ(m))
这个……是因为在{{g0,g1,g2,,gφ(m)1},}其中是模意义下的乘法。在这个群中,如果是有原根的存在,那么g1必定是原根(因为他的指数是最小的正数),那么他就是一个循环群。其实原根就是这个循环群的生成元,那么这个φ(m)阶循环群的生成元数目{gr|,0<rφ(m),rφ(m)},那么这个原根集合的元素个数就是φ(φ(m))
所以我们知道3998344353的原根,同样的33也是998244353的原根。

(写到这里我都快忘了我们原来是要求原根用来Number Theoretic Transforms的了)
好了我们回来。
先写出一个代码求原根。
那么肯定需要分解质因数啦。
大概就应该是这样

View Code
#include <iostream>
#include <cmath> // for sqrt
#include <vector> // for vector
#include <utility> // for pair

template <typename T>
const T fastpow(const T x,size_t k,const T mod)
{
    if(k==0) return 1;
    if(k==1) return x;
    register T res=fastpow(x,k>>1,mod);
    res=(res*res)%mod;
    if(k&1) res=(res*x)%mod;
    return res;
}

template <typename T>
inline std::vector<std::pair<T,size_t> > primes(const T n)
{
    register T x=n;
    std::pair<T,size_t> tmp;
    std::vector<std::pair<T,size_t> > vec;
    for(register T i=2;i*i<=x;++i)
    {
        if(x%i==0)
        {
            tmp.first=i;
            tmp.second=0;
            do {x/=i; tmp.second++;}
            while(x%i==0);
            vec.push_back(tmp);
        }
    }
    tmp.first=x; tmp.second=1;
    if(x!=1) vec.push_back(tmp);
    return vec;
}


template <typename T>
inline T phi(const T n)
{
    std::vector<std::pair<T,size_t> > vec=primes(n);
    register T x=n;
    for(register typename std::vector<std::pair<T,size_t> >::iterator it=vec.begin();it!=vec.end();++it) x-=x/it->first;
    return x;
}

template <typename T>
std::ostream& operator << (std::ostream& stream, std::vector<std::pair<T,size_t> > p)
{
    for(register typename std::vector<std::pair<T,size_t> >::iterator it=p.begin();it!=p.end();++it) stream<<it->first<<' '<<it->second<<std::endl;
    return stream;
}

template <typename T>
inline T is_Primitive_Root(const T i, std::vector<std::pair<T,size_t> >vec,const T PHI,const T mod)
{
    for(register typename std::vector<std::pair<T,size_t> >::iterator it=vec.begin();it!=vec.end();++it)
    {
        if(fastpow(i,PHI/it->first,mod)==1) return it->first;
    }
    return 0;
}

template <typename T>
inline T Primitive_Root(const T n)
{
    const T PHI=phi(n);
    std::vector<std::pair<T,size_t> > vec=primes(PHI);
    for(register T i=2;i<=PHI;++i)
    {
        if(fastpow(i,PHI,n)==1)
        {
            if(is_Primitive_Root(i,vec,PHI,n)==0) return i;
        }
    }
    return 0;
}

int main()
{
    register unsigned long long n,opt=0,a,b;
    register long long unsigned PHI=phi(b);
    do
    {
        switch(opt)
        {
            case 1:
                std::cin>>n;
                std::cout<<Primitive_Root(n)<<std::endl;
                break;
            case 2:
                std::cin>>a>>b;
                PHI=phi(b);
                if(fastpow(a,PHI,b)!=1) std::cout<<-1<<std::endl; else
                std::cout<<is_Primitive_Root(a,primes(PHI),PHI,b)<<std::endl;
                break;
            default:
                break;
        }
        std::cout<<std::endl<<std::endl<<"input opt to ask your question:\n1: input n and output the Primitive Root of n\n2: input a and b to check if a is the Primitive Root of b(if it's true, it will output 0, or it will output n which makes a^{\\varphi(a)/i}\\equiv 1 \\pmod b, or output -1 meaning a^{\\varphi(a)} \\not\\equiv 1 \\pmod b\n"<<std::endl;
    }while(std::cin>>opt);
    return 0;
}

然后我们跑出来998244353的一个原根是3
对于我们要跑NTT的程序来说这个运算可以在编译期通过另一个程序完成,叫做编译器常量,这样求原根的时间复杂度在求NTT的时候就只会是O(1)
(很快啊)
那么我们就可以很愉快地写NTT了

Number Theoretic Transforms

模仿Fast Fourier Transforms
我们把他的表达式按照a的下标奇偶性,或者x的次数的奇偶性进行分类。

f0(x)=a0x0+a2x1+a4x2+f1(x)=a1x0+a3x1+a5x2+f(x)=f0(x2)+xf1(x2)

在FFT过程中我们代入ωni,为了保证二分之后的取值能顺利进行,所以我们取的n就是当前的次数limit
那么我们的NTT很明显,在第一层递归的时候要取n个点,第二层就要取n2个点,第三层要取n4个点。。。
那这个很容易实现,如果最小的单位根是g,那么第一层我们代入g=gφ(m)/n(modm),第二层代入g=(gφ(m)/n)2(modm)
那么不管对于哪一层的gn,我们都能保证(g)n1(modm)

Inverse Number Theoretic Transforms

数论逆变换。
其实和IFFT是一样的。
我们系数转点值的时候其实是一个矩阵乘向量的过程

[h(ωn0)=y0h(ωn1)=y1h(ωn2)=y2h(ωnn)=yn]=[(ωn0)0(ωn0)1(ωn0)2(ωn0)n1(ωn1)0(ωn1)1(ωn1)2(ωn1)n1(ωn2)0(ωn2)1(ωn2)2(ωn2)n1(ωnn1)0(ωnn1)1(ωnn1)2(ωnn1)n1][a0a1a2an1]

然后我们的点值向量再乘起来,得到一个点值向量。
然后这个点值向量就是卷积之后的函数h=fg的系数转化而来的向量。
那么我们这个点值其实也是一个变量矩阵乘一个系数向量的结果。

[h(g0)=y0h(g1)=y1h(g2)=y2h(gn)=yn]=[(g0)0(g0)1(g0)2(g0)n1(g1)0(g1)1(g1)2(g1)n1(g2)0(g2)1(g2)2(g2)n1(gn1)0(gn1)1(gn1)2(gn1)n1][a0a1a2an1]

然后我们把矩阵"除"过去,就会得到

[a0a1a2an1]=[(g0)0(g0)1(g0)2(g0)n1(g1)0(g1)1(g1)2(g1)n1(g2)0(g2)1(g2)2(g2)n1(gn1)0(gn1)1(gn1)2(gn1)n1]1[y0y1y2yn]

注意到n2的正整数次幂,所以我们根据大量的草稿纸可以得出,这个矩阵的逆元就是

[(g0)0(g0)1(g0)2(g0)n1(g1)0(g1)1(g1)2(g1)n1(g2)0(g2)1(g2)2(g2)n1(gn1)0(gn1)1(gn1)2(gn1)n1]1=1n[(g0)0(g0)1(g0)2(g0)(n1)(g1)0(g1)1(g1)2(g1)(n1)(g2)0(g2)1(g2)2(g2)(n1)(gn1)0(gn1)1(gn1)2(gn1)(n1)]

小心前面还有1n
所以原来NTT我们是正着取原根的正整数次幂,现在我们INTT要反着取原根的负整数次幂,其实就是原根的逆元的正整数次幂。
所以我们修改一下,把逆元代入计算就可以得到INTT了。
是不是很像FFT和IFFT?
是的,NTT只是在模意义下用整数优化了FFT,其他细节。。。和FFT一样。

P3803 【模板】多项式乘法(FFT)

View Code
#include <cstdio>
#include <cstring>

namespace number
{

inline unsigned fastpow(unsigned a,unsigned k,const unsigned& mod)
{
	unsigned res=1;
	while(k)
	{
		if(k&1) res=unsigned((long long unsigned)res*a%mod);
		a=unsigned((long long unsigned)a*a%mod);
		k>>=1;
	}
	return res%mod;
}

const unsigned mod=998244353, PHI=998244352, G=3, invG=332748118;
struct number
{
	unsigned dat[2097154],len;
	static unsigned B[2097154];
	inline void read(const unsigned n)
	{
		len=n;
		for(register unsigned i=0;i<=len;++i) scanf("%u",&dat[i]);
		return;
	}
	number()
	{
		len=0;
		memset(dat,0,sizeof(dat));
	}
	number(const number& b)
	{
		len=b.len;
		memcpy(dat,b.dat,sizeof(dat));
	}
	inline number& operator =(long long unsigned n)
	{
		len=-1;
		while(n)
		{
			dat[++len]=n%10;
			n/=10;
		}
		for(register unsigned i=0;(i<<1)<len;++i)
		{
			register unsigned &a=dat[i], &b=dat[len-i];
			a^=b^=a^=b;
		}
		return *this;
	}
	inline void write(const char ch=' ') const
	{
		for(register unsigned i=0;i<=len;++i) printf("%d%c",dat[i],ch);
	}
};
unsigned number::B[2097154];
inline void NTT(unsigned a[],unsigned limit,const unsigned flag=1)
{
	for(register unsigned i=1;i<limit;++i)
	{
		/// 蝴蝶变换
		if(i<number::B[i]) a[i]^=a[number::B[i]]^=a[i]^=a[number::B[i]]; // 一种神奇的交换整数的方法。
	}
	for(register unsigned mid=1;mid<limit;mid<<=1)
	{
		// 用原根代替单位根
		register const unsigned g=fastpow(flag==1? G:invG,PHI/(mid<<1),mod);
		for(register unsigned R=mid<<1,j=0;j<limit;j+=R)
		{
			register unsigned rt=1;
			for(register unsigned opt=0;opt<mid;++opt,rt=unsigned((long long unsigned)rt*g%mod))
			{
				register const unsigned x=a[j+opt],y=unsigned((long long unsigned)rt*a[j+opt+mid]%mod);
				a[j+opt]=(x+y)%mod; a[j+opt+mid]=(x-y+mod)%mod;
			}
		}
	}
}

inline number operator *(const number& _a,const number& _b)
{
	number a=_a,b=_b;
	register unsigned k=0,limit;
	register const unsigned n=a.len,m=b.len;
	while(n+m>=(1u<<++k));
	limit=1<<k;
	memset(number::B,0,sizeof(number::B));
	for(register unsigned i=1;i<=limit;++i) number::B[i]=(number::B[i>>1]>>1)|((i&1)<<(k-1));
	NTT(a.dat,limit);
	NTT(b.dat,limit);
	for(unsigned i=0;i<=limit;++i) a.dat[i]=int((long long)a.dat[i]*b.dat[i]%mod);
	NTT(a.dat,limit,-1);
	int inv=fastpow(limit,mod-2,mod);
	for(register unsigned i=0;i<=n+m;++i) a.dat[i]=int((long long)a.dat[i]*inv%mod);
	a.len=n+m;
	return a;
}


}
unsigned n,m;
number::number a,b;

int main()
{
	scanf("%u%u",&n,&m);
	a.read(n); b.read(m);
	a=a*b;
	a.write();
	return 0;
}
posted @   IdanSuce  阅读(231)  评论(0编辑  收藏  举报
编辑推荐:
· 开发者必知的日志记录最佳实践
· SQL Server 2025 AI相关能力初探
· Linux系列:如何用 C#调用 C方法造成内存泄露
· AI与.NET技术实操系列(二):开始使用ML.NET
· 记一次.NET内存居高不下排查解决与启示
阅读排行:
· 阿里最新开源QwQ-32B,效果媲美deepseek-r1满血版,部署成本又又又降低了!
· 开源Multi-agent AI智能体框架aevatar.ai,欢迎大家贡献代码
· Manus重磅发布:全球首款通用AI代理技术深度解析与实战指南
· 被坑几百块钱后,我竟然真的恢复了删除的微信聊天记录!
· AI技术革命,工作效率10个最佳AI工具
点击右上角即可分享
微信分享提示