数论变换NTT
0. 前言
我们在学Fast Fourier Transforms的时候就会发现输出栏有res[i]=(unsigned long)(a[i].real()/limit+.5)
这里需要加上以保证输出精度
输出精度是怎么产生的呢?
我们用复数运算,这样复数的实部和虚部都需要使用双精度浮点数来运算。
但是浮点数的存储方式和整数不太一样,而且有时候三角函数值出来没准还是无限不循环小数。
误差就这样产生了。
再者我们也有某个神奇的结论:
整形的运算比浮点型运算快5~20倍
所以我们就想,能不能用一些神奇的算法,使得用整数也可以算多项式呢?
我们回想一下为什么多项式里面我们要用复数。
因为我们要用单位根,单位根是复数。
为什么要用单位根?因为
满足这两个性质,我们就可以优化Fourier Transforms
我们看看怎么使用整数也可以满足这两条性质
- 其中第一个性质是用来解决二分时的取值问题的,(ldq大佬:也就是要满足循环卷积),我们只需要保证我们能够快速去到我们需要的值就可以了。
- 至于第二个性质,确实在保证第一个性质下第二个性质很难满足。但是我们退一步想,我们需要是因为我们在里面使用二分,也就是要保证。所以并不是说一定要相反数,我们只需要保证这两个数平方相等就可以了。
至于平方相等,我们很容易想到的是(废话),单位根也正是利用了这个性质:
,这里有,才能满足这条性质。
所以,我们就是要找一个次方变成的整数集合,并且集合中的数满足
这么看来在整数内幂为1的就数只有一个。。。
貌似不太对劲。
没关系,我们看看大部分的题目
答案对取模
这就好办了,在取模意义下大把大把的有,更何况还是一个质数
所以我们只需要找出在模数意义下的集合就可以了
(没见过前言讲这么长的吧?)
原根
那么正整数,模意义下我们知道了很多数的幂都可以是1
但是我们回过头看Fast Fourier Transforms
然后他们又二分,又二分。。。
一路二分下去,最后的时间复杂度就降到了
需要满足的条件就是当前多项式的次数必须是奇数,除非你递归到了边界:次
所以我们就要求对于所有的次数必须都是偶数,也就是一开始多项式的次数必须是的整数次幂。
这么说我们这个幂指数必须是2的正整数次幂
我们就要找出一个,使得
那么这个就可以代替Fast Fourier Transforms里面的次单位根了
这个就叫做原根。
在观看了这篇文章之后相信我对原根有了很深的体会。
首先数学概念里有一个东西叫做阶
阶就是说,有两个互质的数,满足的最小正整数就叫做模意义下的阶,记作
然后我们再来了解一个欧拉函数,表示的是小于且与互质的数的个数。
比如说小于且和10互质的数有共个,那么
至于这个欧拉函数的求法,很简单也很粗暴。
设的质因数的集合是,那么
其实很好理解,不和互质就是的倍数,那么小于里面是的倍数的有个,那么剩下的个就不是的倍数,就有可能和互质。
那么我们对于剩下的质因数也执行这样的操作,我们就可以找齐所有不是的倍数的数,就是和互质的数了。
原根,就是满足的正整数叫做的一个原根。
很明显原根满足
那么我们就可以发现,有原根,当且仅当为奇素数。
至于为什么。。。我也不知道。欢迎评论区补充。
然后的原根会满足,很明显这么说的意思就是
这个还可以证明,用反证法
设,就会有,且
但是,他们的差是不可能达到的,而阶的定义就是满足的最小正整数。
那么都不能满足
故原命题成立,
现在设为奇素数,那么有
那么我们根据原根的性质就会有对于的原根满足互不相同,
但是有,这个就和单位根的性质很相似
那么我们就仿照一下就会有
然后我们再发现给的模数
代入,有,恰好又满足了我们的需要是2的正整数次幂(反正你也取不到题目最多给到就是所以是可以跑的)
现在到了最后一个问题:怎么求原根。
那很明显,我们需要满足的第一个条件就是
但是只满足这个条件的不一定是原根。这只是一个必要条件。
我们当初说的是,这个是最小的使得
如果我们的不是最小的,那么他就不是,这个就不是原根了
举个例子,我们求的原根,有,但是其实
所以,2不是7的原根
那么我们怎么排除这种情况呢?这里用到了裴蜀定理:
裴蜀定理是说明了对任何整数 和它们的最大公约数 ,关于未知数的线性的丢番图方程(称为裴蜀等式)。
丢番图方程又是什么呢?就是不定方程,他的求解只在整数范围内进行,最后这个限制使得丢番图方程求解与实数范围方程求解有根本的不同。
裴蜀定理说明了两个命题
它的一个重要推论是:a,b互质的充分必要条件是存在整数x,y使ax+by=1。
他的证明是这样的:(其实百度百科下面有):
设,则,有,这很简单,命题得证。
接下来一大段都是证明命题的
设为的所有取值中的最小正值,那么再令
那么我们知道这个也是的线性组合,如果我们令,那么就有
又因为应该是的最小正值,那么想要是所有中比小的非负值,就只能的时候了
同理,如果我们设就可以证明,这里不再赘述。
所以,就是的公约数,而是他们的最大公约数,所以有
然后我们也知道对于这个最大公约数,我们一定有,然后就有
所以就可以推断出
所以当时,我们就有,所以,命题得证。
说完了裴蜀定理我们就来看看怎么排除不是原根的情况。
首先我们满足了,然后我们假设,这样就使得不是的原根
根据裴蜀定理,
即
同时我们有
那么很简单了,我们可以推断出,代入
有
只取恒等式最左边和最右边就是,同时那么
又因为,所以,这个保证了
同时我们有。
我们知道妨碍不能成为的原根主要就是因为有不知道哪个数
那么我们根据上面的等价变换(都是等价的吧上面的推理)
我们可以知道,如果存在这样一个使得,那么一定存在也能使
而这个,所以我们只需要验证所有的因数看他能不能使,如果可以,那么这个就不是的原根。
但是要验证所有的因数好麻烦啊。。。
不慌我们有一个很快的方法。
如果合数是的因数,那么肯定都是的因数,只要,那么在模意义下肯定也不会和同余,(反证法:)如果同余那么之中的一个就是1为底的指数式,这样出来的结果在模意义下肯定也是和1同余的,但是他们两个其实和是等价的,这就违背了,所以我们只需要把他所有的质因数都乘起来再验证这个数就可以了。
那么。。。这个数不就是本身吗。。。
我们看看是哪里出了问题:我们验证是的因数,以及,这个大前提是因为有可能会妨碍成为原根。那么这个就一定满足,如果等于了还妨碍什么?
我们把所有质因数都乘起来那么就直接等于了,所以不可行。那么我们可以先把所有质因数都乘起来,然后抠掉一个,再判断,这样就可以知道有没有一个数妨碍成为原根了
如果,其中是的质因数,那么我们就可以知道这个就是的一个原根。
对,我说的是一个原根。
原根可以有个或者个
这个……是因为在其中是模意义下的乘法。在这个群中,如果是有原根的存在,那么必定是原根(因为他的指数是最小的正数),那么他就是一个循环群。其实原根就是这个循环群的生成元,那么这个阶循环群的生成元数目,那么这个原根集合的元素个数就是了
所以我们知道是的原根,同样的也是的原根。
(写到这里我都快忘了我们原来是要求原根用来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;
}
然后我们跑出来的一个原根是
对于我们要跑NTT的程序来说这个运算可以在编译期通过另一个程序完成,叫做编译器常量,这样求原根的时间复杂度在求NTT的时候就只会是
(很快啊)
那么我们就可以很愉快地写NTT了
Number Theoretic Transforms
模仿Fast Fourier Transforms
我们把他的表达式按照的下标奇偶性,或者的次数的奇偶性进行分类。
在FFT过程中我们代入,为了保证二分之后的取值能顺利进行,所以我们取的就是当前的次数
那么我们的NTT很明显,在第一层递归的时候要取个点,第二层就要取个点,第三层要取个点。。。
那这个很容易实现,如果最小的单位根是,那么第一层我们代入,第二层代入
那么不管对于哪一层的和,我们都能保证
Inverse Number Theoretic Transforms
数论逆变换。
其实和IFFT是一样的。
我们系数转点值的时候其实是一个矩阵乘向量的过程
然后我们的点值向量再乘起来,得到一个点值向量。
然后这个点值向量就是卷积之后的函数的系数转化而来的向量。
那么我们这个点值其实也是一个变量矩阵乘一个系数向量的结果。
然后我们把矩阵"除"过去,就会得到
注意到是的正整数次幂,所以我们根据大量的草稿纸可以得出,这个矩阵的逆元就是
小心前面还有
所以原来NTT我们是正着取原根的正整数次幂,现在我们INTT要反着取原根的负整数次幂,其实就是原根的逆元的正整数次幂。
所以我们修改一下,把逆元代入计算就可以得到INTT了。
是不是很像FFT和IFFT?
是的,NTT只是在模意义下用整数优化了FFT,其他细节。。。和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;
}
【推荐】国内首个AI IDE,深度理解中文开发场景,立即下载体验Trae
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 开发者必知的日志记录最佳实践
· 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工具