多项式桶1-NTT（快速数论变换以及任意模数形式）以及FFT的一点杂谈

FFT杂谈
NTT
任意模数NTT
分治FFT
- 多项式求逆
附原根表
参考资料

FFT杂谈

因为之前的我已经不想去动他了。QAQ

其实FFT并不难学，大纲是这样：

前置知识
分治
蝴蝶变换
学完了，没了

之前那个前面不是是抄你谷多项式奆佬的，后面的各种优化没什么实际意义还疯狂掉精度，基本在实战中不会去碰他，所以那篇文章我可以说是弃飞了。

FFT的前置知识并不多，向量和单位根，向量相乘的两个性质一个用相似，一个用两点距离公式爆算都能证明，也就没什么难度。

会用到的性质：

向量相乘模相乘
向量相乘辐角相加

单位根的定义百度一下，你就知道。

当然，单位根用到的性质在下面列出来了：

\(w_{n}^0=1,w_{n}^k=w_{n}^{k\mod n}\)
{\(x∈C\)|\(x=w_{n}^i(i∈[0,n-1],i∈Z)\)}的元素个数是\(n-1\)个（其实意思就是\(w_{n}^i\)互不相等）
\(w_{n}^k=w_{2n}^{2k}\)
\(w_{n}^k=-w_{n}^{k+\frac{n}{2}}(n≡0\mod2)\)
\(\sum\limits_{i=0}^{n-1}(w_{n}^{j-k})^i=n*[j=k]\)（\([A]\)表示如果\(A\)条件成立，则\([A]=1\)，否则为\(0\)）
\(w_{n}^k*w_{n}^t=w_{n}^{k+t}\)

只要用到这几个性质，就可以完美的推导出FFT。

这里也不展开，主要就是奇偶分治。

至于蝴蝶迭代，这个我之前的博客也就，在这里就不展开了，代码以前也是有的，我看了一下，我之前的代码长得还行(⊙﹏⊙)。

但是如果这里就这么一点东西，那我还做个P的杂谈啊(╯‵□′)╯︵┻━┻

首先，关于快速傅里叶变换，你可以觉得其很神奇，因为他巧妙的利用了单位根的性质，之前那个IDFT的推导，也很神奇，好像一气呵成似的。

当然，这里介绍一种不同的理解方式（听说卷积跟线性代数有关，但是那个我不会，就不展开了）

其实你会发现，其实DFT和IDFT就是一个\(1*n\)的矩阵乘一个\(n*n\)的矩阵的过程，所以理论上只要\(IDFT\)和\(DFT\)的矩阵互逆，这样就可以完美的相互抵消。

当然，DFT中的矩阵叫范德蒙德矩阵。

图片来自参考链接的奆佬博客Orz中
在这里插入图片描述
所以，这看似运气十足的FFT背后也有其科学的数学道理。

NTT

全称：快速速论变换

请注意，本文中对于一个多项式模一个整数，只是对每一个系数取模而已。

优点：

全是整数，不用担心掉精度的问题。
可以计算特定模数的卷积。
因为不用计算复数，所以乘法少了不少，貌似常数加快了b（￣▽￣）d　。

缺点：

不能计算小数，过程涉及到取模，所以一旦数域超过了模数，就难以计算。
模数必须是\(k*2^t+1\)的形式，且必须是一个质数。
由于模运算比较慢，所以常数经常不如FFT。

你是不是以为NTT十分的妙？

不，其实就是FFT的过程然后把单位根换了一个更加神奇的东西。

上面也列了，我们会利用到单位根的一些性质，但是其实只要我们能在整数中找到一个同样能满足上述性质的东西，就可以完美的换成整数啦。

当然，在一般的整数域我们很难找到满足条件的数字。

但是如果套上一个模，如果这个模有原根，则有个与原根相关的东西刚好可以满足上面。

设模为\(p\)，\(p\)是形式为\(k*2^t+1\)形式的质数，同时原根为\(g\)。

（原根定义：在模\(p\)的时候，如果对于\(n\)，\(φ(p)\)是最小的正整数 \(t\) 满足 \(n^t≡1\)，则 \(n\) 为\(p的原根\)， \(φ(p)\)表示小于\(p\)的正整数中与 \(p\) 互质的数字个数）

则我们可以把\(w_{n}^i\)替换成\(g^{\frac{p-1}{n}*i}\)（假设\(n|p-1\)），用\(g_{n}^i\)简化符号。

当然，其是否满足上面的式子，我们来看一下：

\(g_{n}^0≡1,g_{n}^i≡g_{n}^{i\mod n}\)
证明：\(g_{n}^0≡g^0=1\)，将 \(i\) 表示成 \(kn+t(0≤t<n)\) 的形式，所以\(g_{n}^i≡g_{n}^{kn+t}≡g^{\frac{p-1}{n}*(kn+t)}≡g^{\frac{p-1}{n}*kn}*g^{\frac{p-1}{n}*t}≡g^{\frac{p-1}{n}*t}≡g_{n}^t\)，证毕，同理，因为都等于\(g_{n}^t\)，所以\(g_{n}^i≡g_{n}^{i±n}\)
{\(x∈Z\)|\(x=g_{n}^i(i∈[0,n-1],i∈Z)\)}的元素个数是\(n-1\)个（意思就是\(g_{n}^i\)互不相等）
证明：假设\(g_{n}^i≡g_{n}^j(i>j)\)
\(g^{\frac{p-1}{n}*i}≡g^{\frac{p-1}{n}}*j\)
\((g^{\frac{p-1}{n}*(i-j)}-1)*g^{\frac{p-1}{n}*j}≡0\)
那么\(p|(g^{\frac{p-1}{n}*(i-j)}-1)*g^{\frac{p-1}{n}*j}\)，但是如果\(p∤g\)，则\(p∤p^k(k∈Z^{*})\)
所以\(p|(g^{\frac{p-1}{n}*(i-j)}-1)\)
所以\(g^{\frac{p-1}{n}*(i-j)}≡1\)，这显然是不成立的，因为\(\frac{p-1}{n}*(i-j)<p-1=φ(p)\)
\(g_{n}^i≡g_{2n}^{2i}\)
证明：\(g_{n}^i≡g^{\frac{p-1}{n}*i}≡g^{\frac{p-1}{2n}*2i}\)。
\(g_{n}^i≡-g_{n}^{i+\frac{n}{2}}(2|n)\)
证明：\(g_{n}^{i+\frac{n}{2}}≡g^{\frac{p-1}{n}*(i+\frac{n}{2})}≡g^{\frac{p-1}{n}*i}*g^{\frac{p-1}{n}*\frac{n}{2}}≡g_{n}^i*g^{\frac{p-1}{2}}=-g_{n}^i\)
当然，至于为什么\(g^{\frac{p-1}{2}}≡-1\)，我们还需要证明一个东西：
对于模质数\(p\)，\(x^2≡1\)的解只有两个\(1,-1\)（或者说是其的同余）。
\(x^2≡1\)
则\((x+1)(x-1)≡0\)，那么\(p|(x+1)\)或者\(p|(x-1)\)，所以\(x≡p+1≡1\)或\(x=p-1≡-1\)。
但是因为\(g\)是原根，所以\(g^{\frac{p-1}{2}}\)只能是\(-1\)。
\(g_{n}^{i}*g_{n}^j≡g^{i+j}_{n}\)。
5,6点换了一个顺序，海星。
证明：展开即可，没什么好说的
\(\sum\limits_{i=0}^{n-1}(g_{n}^{j-k})^i≡n*[j=k]\)
证明：其实完全没有什么必要证明的，反正这个东西只要用相同的性质套到FFT中便可以证明了：
先求\(\sum\limits_{i=0}^{n-1}(g_{n}^j)i\)
当\(j≠0\)时
\(\sum\limits_{i=0}^{n-1}(g_{n}^j)i≡\frac{g_{n}^{nj}-1}{g_{n}^j-1}≡0\)。
但是，如果\(j=0\)
那么显然\(\sum\limits_{i=0}^{n-1}(g_{n}^j)i≡n\)
那么什么时候\(j=0\)呢？也就是上述式子中的\(j=k\)时成立，证毕。

然后就没有什么了，直接把代码中的单位根换成\(g\)，就可以求解了，但是，仅适用于整数卷积。

当然，需要注意的是：

\(g_{n}^i\)中\(n\)必须整除\(p-1\)，而\(FFT\)中\(n=2^i(i∈[1,r])\)，所以\(p=k*2^{t}+1\)就是这么一个道理。
同样的，对于补充到\(2\)的次幂的多项式长度\(len'=2^r\)，我们要求\(r<=t\)，所以对于原来的多项式长度\(len\)，我们要求\(len≤2^t\)，所以\(p\)的\(2^t\)要尽可能的大，那是最好的。
其次，值域的上届减下界\(+1\)要小于\(p\)，因为模\(p\)最怕的就是在值域中有两个数字同余。
对于\(g_{n}^i(i<0)\)，此时\(g_{n}^i≡\frac{1}{g^{\frac{p-1}{n}*(-i)}}≡(\frac{1}{g})^{\frac{p-1}{n}*(-i)}\)，所以只要求出\(g\)的逆就可以计算\(g_{n}^i(i<0)\)了，当然，\(g\)的逆也是一个原根，证明如下：
\(g'^{i}≡(g^i)'\)（其中\(x'\)表示\(x\)的逆）
而\(1\)的逆是其本身，且因为一个数不存在两个不同余的逆，所以很容易知道\(g'\)也是原根。
在IDFT中应当用\(g_{n}^{-i}\)替换\(w_{n}^i\)，应该不会只有我一个憨憨用\(-g_{n}^i\)去代替吧QAQ。
在最后输出除长度\(n\)时，需要注意的是，不能直接除\(p\)，即使值域中的每个数字都小于\(p\)，因为虽然值域的长度小于\(p\)，但是不代表要求值域的长度\(*n\)以后还小于\(p\)，所以还是乘 \(n'\) 的更加保险，当然，如果你这值域长度小的可怜，当我没说。
没什么好提示了，还不快去🐎代码？

当然，这里提一下常用的一个模数：\(998244353=119*2^{23}+1\)，\(2^{23}=8388608\)，他的原根是\(3\)。

一般模数是这个就不用担心长度不够的问题了，因为如果长度还大于\(8388608\)，那多半是你做法的问题。

模板题

时间复杂度：\(O(n\log n)\)

#include<cstdio>
#include<cstring>
#include<cstdlib> 
#define  N  3100000
using  namespace  std;
typedef  long  long  LL;
const  LL  mod=998244353;
inline  LL  ksm(LL  x,LL y)
{
	LL ans=1;
	while(y)
	{
		if(y&1)ans=(ans*x)%mod;
		x=x*x%mod;y>>=1;
	}
	return ans;
}
template<class  T>
inline void zwap(T &x,T &y){x^=y^=x^=y;}
LL  n,m;
LL  limit,l;
LL  a[N],b[N],r[N]/*蝴蝶变换*/;//用来相乘的多项式
LL  g[N],g1=3,g2=ksm(3,mod-2),gg;
void NTT(LL *A,LL type)
{
	type==1?gg=g1:gg=g2;
	for(int  i=1;i<limit;i++)
	{
		if(i<r[i])zwap(A[i],A[r[i]]);
	}
	for(LL  lon=1;lon<limit;lon<<=1)//表示由长度为2^i合并到2^i+1 
	{
		LL  llon=lon<<1;
		g[0]=1;g[1]=ksm(gg,(mod-1)/llon);
		for(LL  j=2;j<lon;j++)g[j]=g[j-1]*g[1]%mod;
		for(LL  L=0;L<limit;L+=llon)
		{
			LL  mid=L+lon;
			for(LL  k=0;k<lon;k++)
			{
				LL  x=A[k+L],y=A[k+mid]*g[k]%mod;
				A[k+L]=(x+y)%mod;A[k+mid]=(x-y+mod)%mod;
			}
		}
	}
}
int  main()
{
	scanf("%lld%lld",&n,&m);
	for(LL  i=0;i<=n;i++)scanf("%lld",&a[i]);
	for(LL  i=0;i<=m;i++)scanf("%lld",&b[i]);
	limit=1;l=0;
	while(limit<=n+m)limit<<=1,l++;
	r[0]=0;for(int i=1;i<limit;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
	NTT(a,1);NTT(b,1);
	for(LL  i=0;i<limit;i++)a[i]=a[i]*b[i]%mod;
	NTT(a,-1);
	LL  ed=n+m,fuck=ksm(limit,mod-2);
	for(LL  i=0;i<=ed;i++)printf("%lld ",a[i]*fuck%mod);//我之前就犯了第6点错误QAQ
	printf("\n");
	return  0;
}

任意模数NTT

前排说明，没有代码，没有代码！

模板题

当然，这道题目用MTT也是可以轻松过掉的

当模数为任意数字的时候，我们便不再考虑用使用题目中的模数来跑NTT了，而应该考虑一种更加神仙的东西，即：中国剩余定理合并。

即我们考虑用大质数作为模数，扩大值域，然后直接计算出其准确的值再进行模，讲真其实MTT又何尝不是呢？

考虑三个大质数\(a,b,c\)，通过中国剩余定理我们知道，我们可以把值域\([0,a-1],[0,b-1],[0,c-1]\)合并成\([0,abc-1]\)。

当然，对于\(a,b,c\)的选取，要注意\(a^2,b^2,c^2\)都要在long long范围内，否则同样会爆掉，这可就功亏一篑了啊。

当然，对于\(a,b,c\)的选取，有个经典组合：998244353,1004535809,469762049，而且这三个模数的原根都是\(3\)，也方便写代码。

然后考虑一下这道题目的值域为多少：

大概是在\(10^9*10^9*10^5=10^{23}\)级别。

所以如果三个大质数的值域在\(10^{24}\)级别是绝对不会爆炸的。

而三个经典质数的值域可以去到：

\(4.7106432275×10^{26}\)

完全没有任何担忧。

至于合并，其实在合并的时候也有技巧可以不去使用快速乘，又是一个黑科技呢。

首先，已知：
\(x≡x_1\mod a\)
\(x≡x_2\mod b\)
\(x≡x_3\mod c\)

前两个直接合并，当然，合并方法可以是CRT合并，可以是EXGCD合并，当然，也可以用下面的奇技淫巧。

\(x≡x_1+k_1a\mod b\)
\(x_2≡x_1+k_1a\mod b\)

所以\(k1≡\frac{x2-x1}{a}\mod b\)

求逆元即可，轻松合并，时间复杂度都是一样的，别想着偷懒。

假定\(x_4≡x_1+k_1a\mod ab\)

可以发现，\(x_4\)在\(\mod a\)和\(\mod b\)的时候都满足对应的要求。

当然，第二次就是\(ab\)和\(c\)的合并。

同理：

\(x≡x_4+k_2ab\mod c\)
\(x_3≡x_4+k_2ab\mod c\)

所以\(k2≡\frac{x3-x4}{ab}\mod c\)

然后同理便可以求出\(k2\)，所以\(x=x_{4}+k_2ab\)，代回去检验也会发现没有任何的问题。

而且，可以发现\(x∈[0,abc),x_{4}∈[0,ab),k_2∈[0,c)\)，除了\(x\)都没有爆long long，可以完美的模其给定的素数\(p\)完成计算，太漂亮。

看完了还不赶紧A了模板题？

🐎代码去！

我就不🐎了。

分治FFT

还是没有🐎代码。

讲真我觉得这玩意实际上就是FFT放在了分治上，甚至你用NTT放在分治上有时候效果都比这个好。

其实就是分治，然后计算左边对右边的贡献的时候用一下FFT统计贡献就行了。

时间复杂度：\(T(n)=n\log{n}+2T(n)\)

这个复杂度啊，那不就是\(O(n\log^2{n})\)

上模板，🐎！

模板

当然，在这里顺便提一下这个模板的另外一个做法。

我的牛顿迭代做法要求复合函数，时间复杂度爆炸了啊QAQ。

多项式求逆

转载自https://www.luogu.com.cn/blog/cain12/solution-p4721

在这里插入图片描述
时间复杂度：\(O(n\log{n})\)

#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
#define md 998244353LL
#define inv(x) power(x,md-2)
#define max_n 1000000
using namespace std;

int n,N,i;

long long num[max_n],ans[max_n];

//const double pi=acos(-1);
//
//struct complex
//{
//  double x,y;
//  complex operator +(complex a)
//  {
//      return (complex){x+a.x,y+a.y};
//  }
//  complex operator -(complex a)
//  {
//      return (complex){x-a.x,y-a.y};
//  }
//  complex operator *(complex a)
//  {
//      return (complex){x*a.x-y*a.y,x*a.y+y*a.x};
//  }
//  complex operator /(double a)
//  {
//      return (complex){x/a,y/a};
//  }
//};

long long power(long long a,long long n)
{
    long long ans=1;
    while(n)
    {
        if(n&1)
          ans=ans*a%md;
        a=a*a%md;
        n>>=1;
    }
    return ans;
}

int re[max_n];
void build_re(int N)
{
    int x,now=N>>1,len=0;
    while(now)
      len++,now>>=1;
    for(int i=0;i<N;i++)
    {
        x=i;now=0;
        for(int j=0;j<len;j++)
          now=(now<<1)|(x&1),x>>=1;
        re[i]=now;
    }
}

//void FFT(complex a[],int len,int opt)
//{
//  for(int i=0;i<len;i++)
//    if(i<re[i])
//      swap(a[i],a[re[i]]);
//  complex wn,w,x;
//  for(int i=2;i<=len;i<<=1)
//    for(int j=(wn=(complex){cos(pi*2/i),opt*sin(pi*2/i)},0);j<len;j+=i)
//      for(int k=(w=(complex){1,0},0);k<i>>1;k++,w=w*wn)
//        x=w*a[j+k+(i>>1)],a[j+k+(i>>1)]=a[j+k]-x,a[j+k]=a[j+k]+x;
//  if(opt==-1)
//    for(int i=0;i<len;i++)
//      a[i]=a[i]/len;
//}
void NTT(long long a[],int len,int opt)
{
    for(int i=0;i<len;i++)
      if(i<re[i])
        swap(a[i],a[re[i]]);
    long long wn,w,x;
    for(int i=2;i<=len;i<<=1)
      for(int j=(wn=power(3,(md-1)/i),(opt==-1?wn=power(wn,md-2):0),0);j<len;j+=i)
        for(int k=(w=1,0);k<i>>1;k++,w=w*wn%md)
          x=a[j+k+(i>>1)]*w%md,a[j+k+(i>>1)]=(a[j+k]-x+md)%md,a[j+k]=(a[j+k]+x)%md;
    if(opt==-1)
    {
        long long inv=power(len,md-2);
        for(int i=0;i<len;i++)
          a[i]=a[i]*inv%md;
    }
}
//void poly_mul(long long a[],long long b[],long long tar[],int len)
//{
//  static long long A[4*max_n],B[4*max_n];
//  memcpy(A,a,sizeof(long long)*len);memcpy(B,b,sizeof(long long)*len);
//  build_re(len);
//  NTT(A,len,1);NTT(B,len,1);
//  for(int i=0;i<len;i++)
//    tar[i]=A[i]*B[i]%md;
//  NTT(tar,len,-1);
//}
void poly_inv(long long a[],long long tar[],int len)
{
    static long long now[4*max_n],ret[4*max_n];
    memset(now,0,sizeof(long long)*len);memset(ret,0,sizeof(long long)*len);
    ret[0]=inv(a[0]);
    for(int i=2;i<=len;i<<=1)
    {
        build_re(i<<1);
        memcpy(now,a,sizeof(long long)*i);
        NTT(ret,i<<1,1);NTT(now,i<<1,1);
        for(int j=0;j<i<<1;j++)
          ret[j]=ret[j]*(2LL-now[j]*ret[j]%md+md)%md;
        NTT(ret,i<<1,-1);
        memset(ret+i,0,sizeof(long long)*i);
    }
    memcpy(tar,ret,sizeof(long long)*len);
}
//void poly_dir(long long a[],long long tar[],int len)
//{
//  for(int i=0;i<len-1;i++)
//    tar[i]=a[i+1]*(i+1)%md;
//}
//void poly_int(long long a[],long long tar[],int len)
//{
//  for(int i=len;i>=1;i--)
//    tar[i]=a[i-1]*inv(i)%md;
//  tar[0]=0;
//}
//void poly_ln(long long a[],long long tar[],int len)
//{
//  static long long now[4*max_n],ret[4*max_n];
//  memset(now,0,sizeof(long long)*len);memset(ret,0,sizeof(long long)*len);
//  poly_inv(a,ret,len);
//  poly_dir(a,now,len);
//  poly_mul(now,ret,ret,len<<1);
//  poly_int(ret,ret,len);
//  memcpy(tar,ret,sizeof(long long)*len);
//}
//void poly_exp(long long a[],long long tar[],int len)
//{
//  static long long now[4*max_n],ret[4*max_n],ln[4*max_n];
//  memset(now,0,sizeof(long long)*len);memset(ret,0,sizeof(long long)*len);memset(ln,0,sizeof(long long)*len);
//  ret[0]=1;
//  for(int i=2;i<=len;i<<=1)
//  {
//      poly_ln(ret,now,i);
//      for(int j=0;j<i;j++)
//        now[j]=(a[j]-now[j]+(j==0)+md)%md;
//      poly_mul(now,ret,ret,i<<1);
//  }
//  memcpy(tar,ret,sizeof(long long)*len);
//}

char Getchar()
{
    return getchar();
    static char buff[1000000],*p,*end=p;
    if(p==end)
      end=buff+fread(p=buff,1,1000000,stdin);
    return *(p++);
}
template<typename T>void read(T &x)
{
    static char rc;
    x=0;rc=Getchar();
    while(!isdigit(rc))
      rc=Getchar();
    while(isdigit(rc))
      x=x*10+rc-'0',rc=Getchar();
}

int main()
{
    read(n);num[0]=1;
    for(i=1;i<n;i++)
      read(num[i]),num[i]=-num[i];
    N=1;
    while(N<n)
      N<<=1;
    poly_inv(num,ans,N);
    for(i=0;i<n;i++)
      cout<<ans[i]<<' ';
    return 0;
}

附原根表

转自https://www.cnblogs.com/Guess2/p/8422205.html

常用素数：
P = 1004535809  ====>  pr = 3
P = 998244353  =====>  pr = 3

//（g 是mod(r*2^k+1)的原根）

素数  r  k  g
3   1   1   2
5   1   2   2
17  1   4   3
97  3   5   5
193 3   6   5
257 1   8   3
7681    15  9   17
12289   3   12  11
40961   5   13  3
65537   1   16  3
786433  3   18  10
5767169 11  19  3
7340033 7   20  3
23068673    11  21  3
104857601   25  22  3
167772161   5   25  3
469762049   7   26  3
1004535809  479 21  3
2013265921  15  27  31
2281701377  17  27  3
3221225473  3   30  5
75161927681 35  31  3
77309411329 9   33  7
206158430209    3   36  22
2061584302081   15  37  7
2748779069441   5   39  3
6597069766657   3   41  5
39582418599937  9   42  5
79164837199873  9   43  5
263882790666241 15  44  7
1231453023109121    35  45  3
1337006139375617    19  46  3
3799912185593857    27  47  5
4222124650659841    15  48  19
7881299347898369    7   50  6
31525197391593473   7   52  3
180143985094819841  5   55  6
1945555039024054273 27  56  5
4179340454199820289 29  57  3

可以根据自己的癖好选择

参考资料

这个师兄真的写的很不错啊QAQ，直接看懂了啊（虽然FFT是在你谷学的）：https://blog.csdn.net/a_forever_dream/article/details/106520376

奆佬博客Orz：https://www.luogu.com.cn/blog/command-block/wei-yun-suan-juan-ji-yu-ji-kuo-zhan

你谷题解

posted @ 2021-02-24 21:26 敌敌畏58 阅读(313) 评论(0) 编辑收藏举报

会员力量，点亮园子希望

刷新页面返回顶部

敌敌畏58