数论变换NTT
0. 前言
我们在学Fast Fourier Transforms的时候就会发现输出栏有res[i]=(unsigned long)(a[i].real()/limit+.5)
这里需要加上\(0.5\)以保证输出精度
输出精度是怎么产生的呢？
我们用复数运算，这样复数的实部和虚部都需要使用双精度浮点数来运算。
但是浮点数的存储方式和整数不太一样，而且有时候三角函数值出来没准还是无限不循环小数。
误差就这样产生了。
再者我们也有某个神奇的结论：
整形的运算比浮点型运算快5~20倍
所以我们就想，能不能用一些神奇的算法，使得用整数也可以算多项式呢？
我们回想一下为什么多项式里面我们要用复数。
因为我们要用单位根，单位根是复数。
为什么要用单位根？因为\((\omega_n)^n=1,\omega_n^{k+\frac n2}=-\omega_n^k\)
满足这两个性质，我们就可以优化Fourier Transforms
我们看看怎么使用整数也可以满足这两条性质
- 其中第一个性质是用来解决二分时的取值问题的，(ldq大佬:也就是要满足循环卷积)，我们只需要保证我们能够快速去到我们需要的值就可以了。
- 至于第二个性质，确实在保证第一个性质下第二个性质很难满足。但是我们退一步想，我们需要\(x_i=-x_{i+\frac n2}\)是因为我们在\(f(x)=f_0(x^2)+x\cdot f_1(x^2)\)里面使用二分，也就是要保证\(x_i^2=x_{i+\frac n2}^2\)。所以并不是说一定要相反数，我们只需要保证这两个数平方相等就可以了。
至于平方相等，我们很容易想到的是\(1^2=1^2\)(废话)，单位根也正是利用了这个性质：
\((\omega_n^i)^2=(\omega_n^i)^2\omega_n^n=(\omega_n^{i+\frac n2})^2\),这里有\(\omega_n^n=1\)，才能满足这条性质。
所以，我们就是要找一个\(n\)次方变成\(1\)的整数集合\(X\)，并且集合中的数满足\(\forall x\in X,\exist! y\in X, y\neq x, x^2=y^2\)
这么看来在整数内幂为1的就数只有一个。。。
貌似不太对劲。
没关系，我们看看大部分的题目
答案对\(998244353\)取模
这就好办了，在取模意义下大把大把的有，更何况\(998244353\)还是一个质数
所以我们只需要找出在模数意义下的\(X\)集合就可以了
（没见过前言讲这么长的吧？）
原根
那么正整数，模意义下我们知道了很多数的幂都可以是1
但是我们回过头看Fast Fourier Transforms
然后他们又二分，又二分。。。
一路二分下去，最后的时间复杂度就降到了\(O(n\log2n)\)
需要满足的条件就是当前多项式的次数必须是奇数，除非你递归到了边界：\(0\)次
所以我们就要求对于所有的次数必须都是偶数，也就是一开始多项式的次数必须是\(2\)的整数次幂。
这么说我们这个幂指数\(n\)必须是2的正整数次幂
我们就要找出一个\(g\),使得\(g^n \equiv 1 \pmod {998244353}, n=2^k,k\in \Z\)
那么这个\(g\)就可以代替Fast Fourier Transforms里面的\(n\)次单位根\(\omega_n\)了
这个\(g\)就叫做原根。
在观看了这篇文章之后相信我对原根有了很深的体会。
首先数学概念里有一个东西叫做阶
阶就是说，有两个互质的数\(a,p>1\),满足\(a^n\equiv 1 \pmod p\)的最小正整数\(n\)就叫做\(a\)模\(p\)意义下的阶，记作\(\delta_p(a)=n\)
然后我们再来了解一个欧拉函数\(\varphi(a)\)，表示的是小于\(a\)且与\(a\)互质的数的个数。
比如说小于\(10\)且和10互质的数有\(1,3,7,9\)共\(4\)个，那么\(\varphi(10)=4\)
至于这个欧拉函数的求法，很简单也很粗暴。
设\(a\)的质因数的集合是\({p_i}\)，那么\(\varphi(a)=a\times\prod(1-\frac 1{p_i})\)
其实很好理解，不和\(a\)互质就是\(p_i\)的倍数，那么小于\(a\)里面是\(p_i\)的倍数的有\(\frac n{p_i}\)个，那么剩下的\(n(1-\frac 1{p_i})\)个就不是\(p_i\)的倍数，就有可能和\(n\)互质。
那么我们对于剩下的质因数也执行这样的操作，我们就可以找齐所有不是\(p_i\)的倍数的数，就是和\(n\)互质的数了。
原根，就是满足\(\delta_m(g)=\varphi(m)\)的正整数\(a\)叫做\(m\)的一个原根。
很明显原根\(g\)满足\(g^{\varphi(m)}\equiv 1 \pmod m\)
那么我们就可以发现，\(m\)有原根，当且仅当\(m=2,4,p^a,2p^a,a\in \N*, p\)为奇素数。
至于为什么。。。我也不知道。欢迎评论区补充。
然后\(m\)的原根会满足\(g^i \not\equiv g^j \pmod m(0 \leq i < j < \varphi(m))\),很明显这么说的意思就是\(g^0\equiv g^{\varphi(m)}\equiv 1 \pmod m\)
这个还可以证明，用反证法
设\(\exist 0\leq i<j<\varphi(m), g^i\equiv g^j \pmod m\),就会有\(g^{j-i}\equiv 1 \pmod m\),且\(j-i\neq 0\)
但是\(i<j<\varphi(m)=\delta_m(g)\),他们的差\(j-i\)是不可能达到\(\delta_m(g)\)的，而阶\(\delta_m(g)\)的定义就是满足\(g^n\equiv 1 \pmod m\)的最小正整数\(n\)。
那么\(\forall i,j, 0<j-i<\delta_m(g)\)都不能满足\(g^{j-i}\equiv 1 \pmod m\)
故原命题成立，\(g^i \not\equiv g^j \pmod m(0 \leq i < j < \varphi(m))\)
现在设\(p\)为奇素数，那么有\(\varphi(p)=p-1\)
那么我们根据原根的性质就会有对于\(p\)的原根\(g\)满足\(g^i,i\in(0,p-1)\)互不相同，
但是有\(g^0\equiv g^{p-1}\equiv 1 \pmod p\),这个就和单位根\(\omega_n^0=\omega_n^n=1\)的性质很相似
那么我们就仿照一下就会有
然后我们再发现给的模数\(998244353=0b111011100000000000000000000001\)
代入\(p\)，有\(p-1=119*2^{23}\),恰好又满足了我们的\(n\)需要是2的正整数次幂（反正你也取不到\(2^{23}\)题目最多给到就是\(2^{21}\)所以是可以跑的）
现在到了最后一个问题：怎么求原根。
那很明显，我们需要满足的第一个条件就是\(g^{\varphi(m)}\equiv 1 \pmod m\)
但是只满足这个条件的不一定是原根。这只是一个必要条件。
我们当初说的是\(\delta_m(g)=\varphi(m)\),这个\(\delta_m(g)\)是最小的\(n\)使得\(g^n\equiv 1 \pmod m\)
如果我们的\(\varphi(m)\)不是最小的，那么他就不是\(\delta_m(g)\),这个\(g\)就不是原根了
举个例子，我们求\(7\)的原根，有\(\varphi(7)=6,2^6\equiv 1 \pmod 7\),但是其实\(2^3\equiv 1 \pmod 7\)
所以\(\delta_7(2)=3\neq \varphi(7)=6\),2不是7的原根
那么我们怎么排除这种情况呢？这里用到了裴蜀定理：
裴蜀定理是说明了对任何整数 \(a,b\)和它们的最大公约数\(\gcd(a,b)=d\) ，关于未知数\(x,y\)的线性的丢番图方程（称为裴蜀等式）。
丢番图方程又是什么呢？就是不定方程，他的求解只在整数范围内进行，最后这个限制使得丢番图方程求解与实数范围方程求解有根本的不同。
裴蜀定理说明了两个命题
\( A:\forall x,y \in \Z, \gcd(a,b)|(ax+by)\\ B:\exist x,y\in \Z, (ax+by)=\gcd(a,b) \)
它的一个重要推论是：a,b互质的充分必要条件是存在整数x,y使ax+by=1。
他的证明是这样的：（其实百度百科下面有）：
设\(d=\gcd(a,b)\)，则\(d|a,d|b\)，有\(d|ax+by(x,y\in\Z)\),这很简单，命题\(A\)得证。
接下来一大段都是证明命题\(B\)的
设\(s=ax_0+by_0\)为\(ax+by\)的所有取值中的最小正值，那么再令\(\displaystyle q=\left\lfloor\frac as \right\rfloor, r=a\%s=a-qs=a-q(ax_0+by_0)=a(1-qx_0)+b(-qy_0)\)
那么我们知道这个\(r\)也是\(a,b\)的线性组合，如果我们令\(x_1=1-qx_0, y_1=-qy_0\),那么就有\(\gcd(a,b))|ax_1+by_1, d|r\)
\(\because r=a\%s\quad\therefore 0\leq r<s\)
又因为\(s\)应该是\(ax+by\)的最小正值，那么\(r\)想要是所有\(ax+by\)中比\(s\)小的非负值，就只能\(a=b=0\)的时候\(r=0\)了
\(\therefore r=a\%s=0, s|a\)
同理，如果我们设\(\displaystyle q=\left\lfloor\frac bs \right\rfloor, r=b\%s\)就可以证明\(s|b\)，这里不再赘述。
所以\(s|a,s|b\),\(s\)就是\(a,b\)的公约数，而\(d\)是他们的最大公约数，所以有\(d\geq s\)
然后我们也知道对于这个最大公约数，我们一定有\(d|a,d|b\),然后就有\(d|(ax_0+by_0), d|s\)
所以就可以推断出\(d\leq s\)
\(\because d\leq s, d\geq s\quad\therefore d=s\)
所以当\(x=x_0,y=y_0\)时，我们就有\(ax+by=\gcd(a,b)\)，所以\(\exist x,y\in\Z, ax+by=\gcd(a,b)\)，命题\(B\)得证。
说完了裴蜀定理我们就来看看怎么排除不是原根的情况。
首先我们满足了\(g^{\varphi(m)}\equiv 1 \pmod m\),然后我们假设\(\exist k\in[0,\varphi(m)), g^k \equiv 1 \pmod m\),这样就使得\(g\)不是\(m\)的原根
根据裴蜀定理，\(\exist x,(-y)\in \Z, kx+\varphi(m)\cdot(-y)=\gcd(k,\varphi(m))\)
即\(kx=\varphi(m)y+\gcd(k,\varphi(m))\)
同时我们有\(\begin{cases} g^k\equiv 1 \pmod m,& g^{kx}\equiv 1^x \pmod m\\ g^{\varphi(m)}\equiv 1\pmod m,& g^{\varphi(m)y}\equiv 1^y\pmod m \end{cases}\)
那么很简单了，我们可以推断出\(g^k\equiv g^{kx}\pmod m\)，代入\(kx=\varphi(m)y+\gcd(k,\varphi(m))\)
有\(\displaystyle g^{k}\equiv g^{\varphi(m)y+\gcd(k,\varphi(m))}\equiv g^{\varphi(m)y}\cdot g^{\gcd(k,\varphi(m))}\equiv 1\cdot g^{\gcd(k,\varphi(m))}\pmod m\)
只取恒等式最左边和最右边就是\(g^k\equiv g^{\gcd(k,\varphi(m))}\pmod m\),同时\(g^k\equiv 1 \pmod m\)那么\(g^{\gcd(k,\varphi(m))}\equiv 1\pmod m\)
又因为\(k<\varphi(m)\),所以\(\gcd(k,\varphi(m))\leq k < \varphi(m)\),这个保证了\(\gcd(k,\varphi(m))\neq\varphi(m)\)
同时我们有\(\gcd(k,\varphi(m))|\varphi(m)\)。
我们知道妨碍\(g\)不能成为\(m\)的原根主要就是因为有不知道哪个数\(k<\varphi(m),g^k\equiv 1\pmod m\)
那么我们根据上面的等价变换（都是等价的吧上面的推理）
我们可以知道，如果存在这样一个\(k\)使得\(g^k\equiv 1\pmod m\)，那么一定存在\(k'=\gcd(k,\varphi(m))\)也能使\(g^{k'}\equiv 1\pmod m\)
而这个\(k'|\varphi(m)\)，所以我们只需要验证\(\varphi(m)\)所有的因数\(n\)看他能不能使\(g^n\equiv 1\pmod m\),如果可以，那么这个\(g\)就不是\(m\)的原根。
但是要验证所有的因数好麻烦啊。。。
不慌我们有一个很快的方法。
如果合数\(n=p_1p_2\)是\(\varphi(m)\)的因数，那么\(p_1,p_2\)肯定都是\(\varphi(m)\)的因数，只要\(g^n\not\equiv 1\pmod m\),那么\(g^{p_1}, g^{p_2}\)在模\(m\)意义下肯定也不会和\(1\)同余，（反证法：）如果同余那么\((g^{p_1})^{p_2},(g^{p_2})^{p_1}\)之中的一个就是1为底的指数式，这样出来的结果在模\(m\)意义下肯定也是和1同余的，但是他们两个其实和\(g^n\)是等价的，这就违背了\(g^n\not\equiv 1\pmod m\),所以我们只需要把他所有的质因数都乘起来再验证这个数就可以了。
那么。。。这个数不就是\(m\)本身吗。。。
我们看看是哪里出了问题：我们验证\(n\)是\(\varphi(m)\)的因数，以及\(g^n\not \equiv 1\pmod m\),这个大前提是因为\(n\)有可能会妨碍\(g\)成为原根。那么这个\(n\)就一定满足\(n<\varphi(m)\)，如果等于了还妨碍什么？
我们把所有质因数都乘起来那么就直接等于\(m\)了，所以不可行。那么我们可以先把所有质因数都乘起来，然后抠掉一个，再判断，这样就可以知道有没有一个数妨碍\(g\)成为原根了
如果\(\forall p_i, g^{\frac{\varphi(m)}{p_i}}\not\equiv1\pmod m\),其中\(p_i\)是\(m\)的质因数，那么我们就可以知道这个\(g\)就是\(m\)的一个原根。
对，我说的是一个原根。
原根可以有\(0\)个或者\(\varphi(\varphi(m))\)个
这个……是因为在\(\{\{g^0,g^1,g^2,\dots,g^{\varphi(m)-1}\},\otimes\}\)其中\(\otimes\)是模意义下的乘法。在这个群中，如果是有原根的存在，那么\(g^1\)必定是原根（因为他的指数是最小的正数），那么他就是一个循环群。其实原根就是这个循环群的生成元，那么这个\(\varphi(m)\)阶循环群的生成元数目\(\{g^r|,0<r\leq \varphi(m), r\perp \varphi(m)\}\),那么这个原根集合的元素个数就是\(\varphi(\varphi(m))\)了
所以我们知道\(3\)是\(998344353\)的原根，同样的\(3^3\)也是\(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\)的次数的奇偶性进行分类。
在FFT过程中我们代入\(\omega_n^i\)，为了保证二分之后的取值能顺利进行，所以我们取的\(n\)就是当前的次数\(limit\)
那么我们的NTT很明显，在第一层递归的时候要取\(n\)个点，第二层就要取\(\frac n2\)个点，第三层要取\(\frac n4\)个点。。。
那这个很容易实现，如果最小的单位根是\(g\)，那么第一层我们代入\(g'=g^{\varphi(m)/n}\pmod m\),第二层代入\(g'=(g^{\varphi(m)/n})^2 \pmod m\)
那么不管对于哪一层的\(g'\)和\(n\)，我们都能保证\((g')^n\equiv 1 \pmod m\)
Inverse Number Theoretic Transforms
数论逆变换。
其实和IFFT是一样的。
我们系数转点值的时候其实是一个矩阵乘向量的过程
然后我们的点值向量再乘起来，得到一个点值向量。
然后这个点值向量就是卷积之后的函数\(h=f\cdot g\)的系数转化而来的向量。
那么我们这个点值其实也是一个变量矩阵乘一个系数向量的结果。
然后我们把矩阵"除"过去，就会得到
注意到\(n\)是\(2\)的正整数次幂，所以我们根据大量的草稿纸可以得出，这个矩阵的逆元就是
小心前面还有\(\frac 1n\)
所以原来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;
}
