NTT笔记和多项式全家桶

1.点值表示法

点值表示法是多项式的另一种表示方法，多项式一般是用表达式表示，对于一个n次的多项式我们可以代入n+1个点来确定这个多项式。
假设A(x)=a0+a1x+a2x²+a3x³…+anxⁿ
这n+1个点确定了这个多项式A（x）
换句话说A这个数唯一地对应了这n+1个点的集合，根据这n+1个点的集合也能反推出A。这种表示方法叫点值表示法

2.原根

定义：
假设一个数g是P的原根，那么gⁱ mod P的结果两两不同，且有 1<g<P，0<i<P，归根到底就是g^(P-1) ≡ 1 (mod P)当且仅当指数为P-1的时候成立.(这里P是素数)。

阶的定义
设m>1,gcd(a,m)=1,使得
成立的最小的r，称为a对模m的阶。
性质：
1.对于正整数m，只有当m=2，4，q^a,
2q^a 时m才有原根（q为奇素数，a≥1）
2.然后呢对于m的原根g，g^i(mod p)（0<=i<=p-2）的值两两不相同，且在[0,p-2]内
常见模数的原根
998244353，1004535809的原根为3
为什么原根能用
FFT中使用的单位根，可以被替换成原根，原根也具有相类似的性质，我们设g为p的原根,令g(N)=g^((p-1)/N),要求p-1能被N整除
性质1 g(N)^N≡1(mod p)
性质2 g(N)^N/2≡-1(mod p)
性质3 g(2N)^2*k≡ g(N)^K(mod p)
性质4 g(N)^k+N/2≡-g(N)^k(mod p）
可见单位根有的性质原根也有，所以我们可以用原根替换单位根来解决计算过程中有模数的问题

3.蝴蝶变换

如果使用递归实现的话会有很大的常数，遇到有些卡常题会TLE，是因为不断复制数组和每次都要排序而导致的，所以我们可以一开始就把最终的位置算出来，然后进行倍增合并，这样常数省下来了，且代码也更好写一点，这个合并的过程就叫做蝴蝶变换
那么如何快速得到每个数最后在那个位置呢？
打印几个结果后，发现是有规律的，规律就是每个数最后位置的下标等于当前下表的二进制的翻转
列如6的二进制为(110),翻转后为(011),所以下标为6的数最后会在下标为3的位置上，这样我们就可以O(n)的得到最终的数列了

void change(ll y[],ll len)
{int i,j,k;//cout<<i<<" "<<j<<endl;;for(int i=1,j=len/2; i<len-1; i++){if(i<j)swap(y[i],y[j]);int k=len/2;while(j>=k){j-=k;k=k/2;}j+=k;}
}

4.NTT实现(O(nlogn))

所以只需要把FFT中的单位根换成原根就可以了,然后计算全部在取模的环境下就可以了
DFT
IDFT

a是原根，a^-ik是代表在mod M下的逆元
模板题
代码实现

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const long double PI = acos(-1.0);
const ll p=998244353;//模数
const ll G=3;//原根
typedef long long ll;
typedef long double db;
ll ksm(ll a,ll b)
{ll ans=1;while(b){if(b&1)ans=ans*a;a=a*a;b=b>>1;}return ans;
}
void change(ll y[],ll len)
{int i,j,k;//cout<<i<<" "<<j<<endl;;for(int i=1,j=len/2; i<len-1; i++){if(i<j)swap(y[i],y[j]);int k=len/2;while(j>=k){j-=k;k=k/2;}j+=k;}
}
void ntt(ll y[], int len, int on)
{change(y, len);for (int h = 2; h <= len; h <<= 1){ll wn=ksm(G,(p-1)/h);if(on==-1){wn=ksm(wn,p-2);}for (int j = 0; j < len; j += h){ll w=1;for (int k = j; k < j + h / 2; k++){ll u = y[k]%p;ll t = w * y[k + h / 2]%p;y[k] = (u + t)%p;y[k + h / 2] = (u - t+p)%p;w = (w * wn)%p;}}}if (on == -1){ll len_inv = ksm(len,p-2);//N的逆元（N是limit, 指的是2的整数幂）for(int i = 0; i < len; ++ i)y[i] = (y[i] * len_inv) % p;//NTT还是要除以n的，但是这里把除换成逆元了，inv就是n在模p意义下的逆元}
}
ll a[4000010],b[4000010];
void prepare(int len1, int len2)
{int c=1;while(c<(max(len1,len2))*2)c=c<<1;// cout<<c<<endl;for(int i=len1; i<c; i++)a[i]=0;for(int i=len2; i<c; i++)b[i]=0;ntt(a,c,1);ntt(b,c,1);//cout<<1<<endl;for(int i=0; i<c; i++)a[i]=a[i]*b[i]%p;ntt(a,c,-1);for(int i=0; i<len1+len2-1; i++){cout<<a[i]<<" ";}
}
int main()
{int n,m;cin>>n>>m;for(int i=0; i<=n; i++){int x;cin>>x;a[i]=x;}for(int i=0; i<=m; i++){int x;cin>>x;b[i]=x;}prepare(n+1,m+1);
}

5.多项式求逆(倍增 O(nlogn))

定义：
对于多项式A(x),如果存在A(x)B(x)≡1(mod xⁿ),那么称B(x)是在mod xⁿ下A(x)的逆元
这里的模xⁿ是指舍弃所有x的次数大于等于n的项，当然也可以理解为所有x的次数大于等于n的项的系数赋为零。
性质1 设有整数m，1<=m<=n，则A(x)B(x)≡(mod x^m)，在mod xⁿ的基础上舍去（m到n）的项不会影响结论
证明
我们设B(x)为A(x)在mod xⁿ的意义下的逆元，C(x)为A(x)在mod x^(m),（m为n除2向上取整）
由
A(x)^-1≡B(x)(mod xⁿ)
A(x)^-1≡C(x)(mod x^m)
推出
B(x)≡C(x)(mod x^m)
即(B(x)-C(x))≡0(mod x^m)
(B(x)-C(x))²≡0(mod xⁿ)
得
B(x)²+C(x)²-2B(x)C(x)≡0(mod xⁿ)
两边同时乘以A(x)
A(x)B(x)B(x)+A(x)C(x)C(x)-2A(x)B(x)C(x)≡0(mod xⁿ)
根据定义A(x)B(x)≡1(mod xⁿ)可得
B(x)+A(x)C(x)C(x)-2C(x)≡0(mod xⁿ)
即B(x)=C(x)(2-A(x)C(x))(mod xⁿ)
所以可以通过求出C(x)，然后倍增求出B(x)

题目链接
代码实现

void ntt_inv(ll *a,ll *b,ll len)
{if(len==1){b[0]=inv(a[0])%p;return ;}ntt_inv(a,b,(len+1)>>1);ll cnt=1;while(cnt<=len*2)cnt<<=1;for(int i=0;i<len;i++)x[i]=a[i];for(int i=len;i<cnt;i++)x[i]=0;ntt(b,cnt,1),ntt(x,cnt,1);for(int i=0;i<cnt;i++)b[i]=(2-b[i]*x[i]%p+p)%p*b[i]%p;ntt(b,cnt,-1);for(int i=len;i<cnt;i++)b[i]=0;
}

6.多项式开方(O(nlog²n))

定义：
设有多项式A(x)，求B(x)使得A(x)≡B(x)²(mod xⁿ)
证明
和多项式求逆相似,假设我们已经求出G(x)
1.A(x)≡G(x)²(mod x^m)(m为n除2向上取整)
2.A(x)-G(X)²≡0(mod x^m)
两边同时平方
3.(A(x)-G(x)²)²≡0(mod xⁿ)
两边同时加上4G(x)²A(x)(mod xⁿ)
4.(A(x)+G(x)²)²≡4G(x)²A(x)(mod xⁿ)
5.(A(x)+G(x)²)²/(4G(x)²)≡A(x)(mod xⁿ)
6.((A(x)+G(x)²)/(2G(x)))²≡A(x)
根据定义
B(x)=(A(x)+G(x)²)/(2G(x))
题目
代码实现
分子部分直接倍增递推求解，分母部分用多项式求逆得出

void ntt_sqrt(ll *a,ll *b,ll len)
{if(len==1){b[0]=1;return ;}ntt_sqrt(a,b,(len+1)>>1);ll cnt=1;while(cnt<=len*2)cnt<<=1;for(int i=0;i<len;i++)c[i]=a[i];for(int i=len;i<cnt;i++)c[i]=0;for(int i=0;i<cnt;i++)d[i]=0;ntt_inv(b,d,len);ntt(c,cnt,1);//A[i]ntt(b,cnt,1);//B[i]ntt(d,cnt,1);//B[i]的乘法逆元//公式B[i]=(A[i]+B[i]^2)/(2*B[i])=(c[i]*d[i]+b[i])*inv_2for(int i=0;i<cnt;i++){b[i]=(c[i]*d[i]%p+b[i])%p*inv_2%p;}ntt(b,cnt,-1);for(int i=len;i<cnt;i++)b[i]=0;
}

7.多项式指数（B(x)=e^A(x)）

证明
实现要用到牛顿迭代的公式

设G(x)=e^F(x),直接求解不好接，我们可以两边取对数
即lnG(x)=F(x）
即lnG(x)-F(x)=0,
然后可以用牛顿迭代求零点的方法求，设f(G(x))=lnG(x)-F(x)
f(G(x))'=1/G(x)
带入牛顿迭代的公式，过程如下图

然后递推求解即可
代码实现

void ntt_exp(ll *a,ll *b,ll len)
{if(len==1){b[0]=1;return ;}ntt_exp(a,b,(len+1)>>1);ll cnt=1;while(cnt<len*2)cnt<<=1;for(int i=0;i<cnt;i++)f[i]=g[i]=0;for(int i=0;i<len;i++)f[i]=a[i];ntt_ln(b,g,c,d,len);//求ln//b[i]=G[i]//g[i]=lnG[i]//f[i]=F[i];ntt(b,cnt,1);ntt(g,cnt,1);ntt(f,cnt,1);for(int i=0;i<cnt;i++)b[i]=(1-g[i]+f[i]+p)%p*b[i]%p;ntt(b,cnt,-1);for(int i=len;i<cnt;i++)b[i]=0;
}

8.多项式对数（B(x)=lnA(x)）

定义：
设有多项式A(X),求B(x)使得B(x)≡lnA(x)(mod xⁿ)
因为无法直接求对数，我们可以先求导，在积分求出目标多项式
即(lnA(x))求导为(A’(x)/A(x))
然后积分即可得到B(x)
要用到的求导公式
(x^a)'=ax^a-1
要用的积分公式
积分(x^a)=1/(a+1)*x^a+1
代码实现

void ntt_ln(ll *a,ll *b,ll *c,ll len)
{ll cnt=1;while(cnt<len*2)cnt<<=1;ntt_qd(a,c,cnt);ntt_inv(a,b,cnt);for(int i=len;i<cnt;i++)c[i]=0;for(int i=len;i<cnt;i++)b[i]=0;ntt_mul(c,b,cnt);ntt_jf(c,b,cnt);
}

9.多项式快速幂

B(x)=A^k(X)
两边同时取对数
lnB(x)=klnA(x)
所以B(x)=e^klnA(x)
使用exp函数和ln函数即可
代码实现

void ntt_ksm(ll *a,ll *b,ll len,ll k)
{ll cnt=1;while(cnt<len*2)cnt<<=1;ntt_ln(a,h,j,len);for(int i=0;i<len;i++)h[i]=h[i]*k%p;memset(b,0,sizeof(b));ntt_exp(h,b,len);
}
int main()
{char s;ll k=0;scanf("%lld%s",&n,s);while(isdigit(s=getchar()))k=(k*10+(s^48))%p;for(int i=0; i<=n-1; i++){cin>>a[i];}// ntt_sqrt(a,b,c,d,n);// ntt_inv(a,b,c,n);// ntt_ln(a,b,c,d,n);ntt_ksm(a,b,n,k);for(int i=0;i<n;i++){cout<<b[i]<<" ";}cout<<endl;
}

10. 多项式全家桶

代码

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const long double PI = acos(-1.0);
const ll p=998244353;//模数
const ll G=3;//原根
const ll Gi=332748118;//逆元
const ll inv_2=499122177;
typedef long long ll;
typedef long double db;
int n,m;
ll ksm(ll a,ll b)
{ll ans=1;while(b){if(b&1)ans=(ans*a)%p;a=(a*a)%p;b=b>>1;}return ans%p;
}
ll inv(ll x){return ksm(x,p-2);}
void change(ll y[],ll len)
{int i,j,k;//cout<<i<<" "<<j<<endl;;for(int i=1,j=len/2; i<len-1; i++){if(i<j)swap(y[i],y[j]);int k=len/2;while(j>=k){j-=k;k=k/2;}j+=k;}
}
void ntt(ll *y, int len, int on)
{change(y, len);// for(int i=0;i<len;i++)cout<<y[i]<<" ";// cout<<endl;for (int h = 2; h <= len; h <<= 1){ll wn=ksm(G,(p-1)/h);if(on==-1){wn=inv(wn);}for (int j = 0; j < len; j += h){ll w=1;for (int k = j; k < j + h / 2; k++){ll u = y[k];ll t = w * y[k + h / 2]%p;y[k] = (u + t)%p;y[k + h / 2] = (u - t+p)%p;w = (w * wn)%p;}}}if (on == -1){ll len_inv = ksm(len,p-2);//N的逆元（N是limit, 指的是2的整数幂）for(int i = 0; i < len; ++ i){// cout<<y[i]<<" ";y[i] = (y[i] * len_inv) % p;//NTT还是要除以n的，但是这里把除换成逆元了，inv就是n在模p意义下的逆元}// cout<<endl;}
}
const ll maxn=4000010;
ll a[maxn],b[maxn];
ll x[maxn];
ll c[maxn],d[maxn],e[maxn];
ll f[maxn],g[maxn];
ll h[maxn],j[maxn];
void ntt_inv(ll *a,ll *b,ll len)
{if(len==1){b[0]=inv(a[0])%p;return ;}ntt_inv(a,b,(len+1)>>1);ll cnt=1;while(cnt<=len*2)cnt<<=1;for(int i=0;i<len;i++)x[i]=a[i];for(int i=len;i<cnt;i++)x[i]=0;ntt(b,cnt,1),ntt(x,cnt,1);for(int i=0;i<cnt;i++)b[i]=(2-b[i]*x[i]%p+p)%p*b[i]%p;ntt(b,cnt,-1);for(int i=len;i<cnt;i++)b[i]=0;
}
void ntt_sqrt(ll *a,ll *b,ll len)
{if(len==1){b[0]=1;return ;}ntt_sqrt(a,b,(len+1)>>1);ll cnt=1;while(cnt<=len*2)cnt<<=1;for(int i=0;i<len;i++)c[i]=a[i];for(int i=len;i<cnt;i++)c[i]=0;for(int i=0;i<cnt;i++)d[i]=0;ntt_inv(b,d,len);ntt(c,cnt,1);//A[i]ntt(b,cnt,1);//B[i]ntt(d,cnt,1);//B[i]的乘法逆元//公式B[i]=(A[i]+B[i]^2)/(2*B[i])=(c[i]*d[i]+b[i])*inv_2for(int i=0;i<cnt;i++){b[i]=(c[i]*d[i]%p+b[i])%p*inv_2%p;}ntt(b,cnt,-1);for(int i=len;i<cnt;i++)b[i]=0;
}
void ntt_qd(ll *a,ll *b,ll len)
{for(int i=1;i<len;i++){b[i-1]=i*a[i]%p;}b[len-1]=0;
}
void ntt_jf(ll *a,ll *b,ll len)
{for(int i=1;i<len;i++){b[i]=a[i-1]*inv(i)%p;}b[0]=0;
}
void ntt_mul(ll *a,ll *b,ll len)
{ll cnt=len;ntt(a,cnt,1);ntt(b,cnt,1);for(int i=0;i<cnt;i++){a[i]=a[i]*b[i]%p;}ntt(a,cnt,-1);
}
void ntt_ln(ll *a,ll *b,ll *c,ll len)
{ll cnt=1;while(cnt<len*2)cnt<<=1;ntt_qd(a,c,cnt);ntt_inv(a,b,cnt);for(int i=len;i<cnt;i++)c[i]=0;for(int i=len;i<cnt;i++)b[i]=0;ntt_mul(c,b,cnt);ntt_jf(c,b,cnt);
}
void ntt_exp(ll *a,ll *b,ll len)
{if(len==1){b[0]=1;return ;}ntt_exp(a,b,(len+1)>>1);ll cnt=1;while(cnt<len*2)cnt<<=1;for(int i=0;i<cnt;i++)f[i]=g[i]=0;for(int i=0;i<len;i++)f[i]=a[i];ntt_ln(b,g,c,len);//求ln//b[i]=G[i]//g[i]=lnG[i]//f[i]=F[i];ntt(b,cnt,1);ntt(g,cnt,1);ntt(f,cnt,1);for(int i=0;i<cnt;i++)b[i]=(1-g[i]+f[i]+p)%p*b[i]%p;ntt(b,cnt,-1);for(int i=len;i<cnt;i++)b[i]=0;
}
void ntt_ksm(ll *a,ll *b,ll len,ll k)
{ll cnt=1;while(cnt<len*2)cnt<<=1;ntt_ln(a,h,j,len);for(int i=0;i<len;i++)h[i]=h[i]*k%p;memset(b,0,sizeof(b));ntt_exp(h,b,len);
}
int main()
{/*char s;ll k=0;scanf("%lld%s",&n,s);while(isdigit(s=getchar()))k=(k*10+(s^48))%p;*/cin>>n;for(int i=0; i<=n-1; i++){cin>>a[i];}// ntt_sqrt(a,b,c,d,n);// ntt_inv(a,b,c,n);// ntt_ln(a,b,c,d,n);ntt_ksm(a,b,n,k);for(int i=0;i<n;i++){cout<<b[i]<<" ";}cout<<endl;
}