多项式-FFT-NTT-数论-学习笔记
i207M
2018-11-20 19:09:38
**数组下标01搞清楚**
关于
```cpp
for(; fn<=(n<<1); fn<<=1);
```
那个等号是否加，情况是这样，如果输入n，给你n+1个数，等号必须加；给你n个数，可加可不加。
**在n不同时，r也不同！**
**$\pi$的值一定要用double存！**
------
## FFT
像我这种蒟蒻写一篇介绍FFT的博客有些超出我的能力范围了，就放一个很棒的博客的链接，我就是和他学的。
[FFT-学习笔记](https://www.luogu.org/blog/command-block/fft-xue-xi-bi-ji)
实现？
为了减小常数，使用迭代法。首先我们要找到每个元素最终的位置，其实这就是它二进制表示下的倒序，可以通过**蝴蝶变换**得到。
然后就是大力推式子了，具体的就看代码吧。
```cpp
void dft(cpx f[],int op)
{
for(ri i=0; i<n; ++i) if(i<r[i]) swap(f[i],f[r[i]]);
for(ri p=2; p<=n; p<<=1)
{
int len=p>>1;
cpx nx=cpx(cos(pi/len),op*sin(pi/len));
for(ri k=0; k<n; k+=p)
{
cpx w=cpx(1,0);
for(ri i=k; i<k+len; ++i)
{
cpx t=w*f[i+len];
f[i+len]=f[i]-t;
f[i]=f[i]+t;
w=w*nx;
}
}
}
}
```
我们一般写的FFT是处理下标和为定值的情况，即$c[k]=\sum a_ib_{k-i}$，如果要求下标的差为定值的情况$c[k]=\sum a_ib_{i-k}$，只需要把一个数组反过来，这样得到的答案数组再反过来即为所求。
------
## NTT
就是模意义下的FFT，可以在已知最大数不会爆模数的情况下，代替FFT的浮点运算来加快。
唯一的变化：单位根是什么？设md的原根为g，则$w=g^{(md-1)/p}$，这要求这个模数要有原根，且$md-1$有很多质因数2；
```cpp
void dft(int f[],int op)
{
for(ri i=0; i<n; ++i) if(i<r[i]) swap(f[i],f[r[i]]);
for(ri p=2; p<=n; p<<=1)
{
int len=p>>1,nx=qpow(op==1?3:332748118,(md-1)/p);
for(ri k=0; k<n; k+=p)
{
int w=1;
for(ri i=k; i<k+len; ++i)
{
int t=(LL)f[i+len]*w%md;
f[i+len]=(f[i]-t+md)%md;
f[i]=(f[i]+t)%md;
w=(LL)w*nx%md;
}
}
}
}
```
------
## 任意模数NTT-MTT
实质就是NTT+CRT合并。常用的NTT模数为
```cpp
const int md[3]= {469762049,998244353,1004535809};
```
对于模板题来说，最大的数字为$10^{18+5=23}$，比long long还要大，所以我们选择三个乘积大于$10^{24}$的质数作为模数即可。
如何合并？这三个质数乘起来也爆long long了（~~高精？~~）
首先用中国剩余定理合并前两个，得到：
$\left\{\begin{aligned}x \equiv A \pmod M \\x \equiv a_3 \pmod {m_3} \\\end{aligned}\right.$
然后设$ans = kM + A$
$k \equiv (a_3 - A)M^{-1} \pmod {m_3}$
回代求值即可。
注意：
1. CRT合并时，逆元是对于m不是M的。
2. CRT合并要用快速乘。
3. k的值是在模m3意义下的。
------
**求逆、开根的时候，建议把数组的长度设为2的整次幂，这样不会出现奇怪的问题**
## 多项式求逆
[感谢miskcoo](http://blog.miskcoo.com/2015/05/polynomial-inverse)
为什么会有无限项？因为多项式中含有未知数，它们相乘怎么也不会不含未知数，所以这个逆元实际上是无限项之和为1，于是我们就要求出$\mod x^n$的结果。
**由这个过程可以看出，一个多项式有没有逆元完全取决于其常数项是否有逆元**
------
## 多项式除法
~~舒服！一遍A！~~
[感谢miskcoo](http://blog.miskcoo.com/2015/05/polynomial-division)
为了方便，将式子写为$A(x)=B(x)C(x)+D(x)$
核心思想：通过**系数反转（${A}^R$）**，消除$D(x)$的影响，算出C，回代得到D。
这个值得放一下代码+注释：
```cpp
void dft(int n,int f[],int op)
{
for(ri i=0; i<n; ++i) r[i]=(r[i>>1]>>1)|((i&1)?n>>1:0);
for(ri i=0; i<n; ++i) if(i<r[i]) swap(f[i],f[r[i]]);
for(ri p=2; p<=n; p<<=1)
{
int len=p>>1,nx=qpow(op==1?3:332748118,(md-1)/p);
for(ri k=0; k<n; k+=p)
{
int w=1;
for(ri i=k; i<k+len; ++i)
{
int t=(LL)w*f[i+len]%md;
f[i+len]=(f[i]-t+md)%md;
f[i]=(f[i]+t)%md;
w=(LL)w*nx%md;
}
}
}
if(op==-1)
{
int inv=qpow(n,md-2);
for(ri i=0; i<n; ++i) f[i]=(LL)f[i]*inv%md;
}
}
void inverse(int n,int a[],int b[])
{
if(n==1)
{
b[0]=qpow(a[0],md-2);
return;
}
inverse((n+1)>>1,a,b); // 递归
int fn=1;
for(; fn<(n<<1); fn<<=1);
for(ri i=0; i<n; ++i) tmpc[i]=a[i];
for(ri i=n; i<fn; ++i) tmpc[i]=0; // 记得将高位清零
dft(fn,tmpc,1), dft(fn,b,1);
for(ri i=0; i<fn; ++i) b[i]=b[i]*(2-(LL)tmpc[i]*b[i]%md+md)%md;
dft(fn,b,-1);
for(ri i=n; i<fn; ++i) b[i]=0; // 记得将高位清零
}
void divi(int n,int a[],int m,int b[],int c[],int d[])
{
int t=n-m+1,fn=1;
for(; fn<(t<<1); fn<<=1);
for(ri i=0; i<t; ++i) tmpa[i]=a[n-i],tmpb[i]=b[m-i]; // 反转后存入数组
inverse(t,tmpb,c);
dft(fn,tmpa,1), dft(fn,c,1);
for(ri i=0; i<fn; ++i) c[i]=(LL)c[i]*tmpa[i]%md; // 乘逆元
dft(fn,c,-1);
reverse(c,c+t);
for(ri i=t; i<fn; ++i) c[i]=0; // 得到C后高位清零
for(fn=1; fn<(m<<1); fn<<=1);
for(ri i=0; i<fn; ++i) tmpa[i]=tmpb[i]=0;
for(ri i=0; i<m; ++i) tmpb[i]=b[i],tmpa[i]=c[i];
dft(fn,tmpb,1), dft(fn,tmpa,1);
for(ri i=0; i<fn; ++i) tmpa[i]=(LL)tmpa[i]*tmpb[i]%md; // 计算 B*C
dft(fn,tmpa,-1);
for(ri i=0; i<m; ++i) d[i]=(a[i]-tmpa[i]+md)%md;
}
```
------
~~多项式开根、指数、对数：未完待续~~
## 多项式开根
~~手推出来了好开心~~
据说常数是DFT的14倍...
```cpp
void Sqrt(int n,int a[],int b[])
{
if(n==1)
{
b[0]=1;
return;
}
Sqrt(n>>1,a,b);
inv(n,b,tmpb);
int m=pre(n<<1);
for(ri i=0; i<n; ++i) tmpa[i]=a[i];
dft(tmpa,m,1); dft(tmpb,m,1); dft(b,m,1);
for(ri i=0; i<m; ++i) b[i]=((LL)tmpa[i]*tmpb[i]%md+b[i])%md*inv2%md;
dft(b,m,-1);
for(ri i=n; i<m; ++i) b[i]=0;
for(ri i=0; i<m; ++i) tmpa[i]=tmpb[i]=0;
}
```
------------------
~~本来不打算去学对数的，但是看起来很简单就学了~~
## 多项式对数
$lnA(x)=\int (lnA(x))'=\int\frac{A'(x)}{A(x)}$
因为多项式很容易求导、积分，所以这种方法比较妙。
$b_i=(i+1)a_{i+1}$
$c_i=\frac {a_{i-1}}i$
关于取对数时要保证$a_0=1$的原因：当x=0时，$\ln A(0)=\ln 1=0$，我们求出的答案B(x)恰好常数项为0，x=0时$B(x)=0$，符合题意。
否则不定积分的常数项$C=\ln a_0 \mod p$，在模意义下取自然对数...反正我不会。
~~指数需要牛顿迭代法，先缓缓~~
------
## cdq分治FFT
$f[i]=\sum_{j=1}^if[i-j]g[j]$
f这个多项式没有完全给你，只能一边算一边求。像这种问题，可以使用——cdq分治！
每次统计左区间对右区间的贡献。
把$f[l...mid]$拿出来，放到$a[0...(mid-l)]$里。
把$g[0...r-l]$拿出来，放到$b$。
把它们卷起来，结果放到$a$。
把$a[(mid-l+1)...(r-l)]$对应加到$f[(mid+1)...r]$上。
```cpp
int a[N*4],b[N*4];
int f[N],g[N];
void fft(int l,int mid,int R)
{
int len=R-l+1,m=1;
for(ri i=l; i<=mid; ++i) a[i-l]=f[i];
for(ri i=0; i<len; ++i) b[i]=g[i];
len+=mid-l+1;
for(; m<=len; m<<=1);
for(ri i=1,hb=m>>1; i<m; ++i) r[i]=(r[i>>1]>>1)|((i&1)?hb:0);
dft(a,m,1), dft(b,m,1);
for(ri i=0; i<m; ++i) a[i]=(LL)a[i]*b[i]%md;
dft(a,m,-1);
for(ri i=mid+1; i<=R; ++i) (f[i]+=a[i-l])%=md;
for(ri i=0; i<m; ++i) a[i]=b[i]=0;
}
void cdq(int l,int r)
{
if(l>=r) return;
int mid=(l+r)/2;
cdq(l,mid);
fft(l,mid,r);
cdq(mid+1,r);
}
```
------------------
## 多项式快速幂
求一个次数为n的多项式的k次方，保证$nk<=10^6$
当然可以使用传说中的（真·）分治FFT来搞，复杂度$O(N\log N),N=nk$。
但是我们还有更加优美的方法：
直接把原多项式dft成点值表达，然后把每个位置变成$k$次方，再idft回去即可。记得要把范围扩大。
求多项式快速幂模另一个多项式的结果，就想快速幂一样，步步取模即可。写起来挺恶心。
------
## 多项式牛顿迭代
### 牛顿迭代
牛顿迭代的精度每迭代一次，精度可以近似看作翻倍。
*进入正题*。
我们已知多项式$G(x)$，要求出$F(x)$使得$G(F(x))=0 \mod x^n$。
好，接下来我们实际应用多项式牛顿迭代来解决多项式开方：
有一个疑问：为什么在求导的时候，忽略了$A(x)$？原因是我们求导的时候，只把$F(x)$看成$x$，不含$F(x)$的想都看成常数。
你问我为什么不用多项式牛顿迭代解决多项式求逆？因为它需要多项式求逆啊！
**注意：** 多项式牛顿迭代不一定能完美拟合所有函数！
----------------
## 多项式指数函数
$\exp$常数项是1
```cpp
void Inv(int n,int a[],int b[])
{
if(n==1)
{
b[0]=1;
return;
}
Inv(n>>1,a,b);
static int c[N];
for(ri i=0; i<n; ++i) c[i]=a[i];
int m=pre(n<<1);
dft(c,m,1); dft(b,m,1);
for(ri i=0; i<m; ++i) b[i]=(2-(LL)b[i]*c[i]%md+md)*b[i]%md;
dft(b,m,-1);
for(ri i=n; i<m; ++i) b[i]=0;
for(ri i=0; i<m; ++i) c[i]=0;
}
int iv[N];
void Ln(int n,int a[],int b[])
{
static int c[N];
for(ri i=1; i<n; ++i) b[i-1]=(LL)a[i]*i%md;
Inv(n,a,c);
int m=pre(n<<1);
dft(b,m,1); dft(c,m,1);
for(ri i=0; i<m; ++i) b[i]=(LL)b[i]*c[i]%md;
dft(b,m,-1);
for(ri i=n-1; i>=1; --i) b[i]=(LL)b[i-1]*iv[i]%md;
b[0]=0;
for(ri i=n; i<m; ++i) b[i]=0;
for(ri i=0; i<m; ++i) c[i]=0;
}
void Exp(int n,int a[],int b[])
{
if(n==1)
{
b[0]=1;
return;
}
Exp(n>>1,a,b);
static int c[N];
Ln(n,b,c);
for(ri i=0; i<n; ++i) c[i]=(a[i]-c[i]+md)%md;
(++c[0])%=md;
int m=pre(n<<1);
dft(b,m,1); dft(c,m,1);
for(ri i=0; i<m; ++i) b[i]=(LL)b[i]*c[i]%md;
dft(b,m,-1);
for(ri i=n; i<m; ++i) b[i]=0;
for(ri i=0; i<m; ++i) c[i]=0;
}
int n;
int a[N],b[N];
signed main()
{
#ifdef M207
freopen("in.in","r",stdin);
// freopen("out.out","w",stdout);
#endif
in(n);
for(ri i=0; i<n; ++i) in(a[i]);
int m=1; for(; m<n; m<<=1);
for(ri i=1; i<m; ++i) iv[i]=qpow(i,md-2);
Exp(m,a,b);
print(b,n);
return 0;
}
```
**通过ln和exp解决最低次项不为1的多项式快速幂：除以最低次项**
多项式$\ln,\exp$主要是用来解决排列问题（多重集合的排列），对于图的问题，可用来切换联通与不联通。eg.设$g$为n个点的无向图的个数，$f$为n个点的无向连通图的个数，则$g=\exp f,f=\ln g$
例题：城市规划
我们来证明一下上个式子（by Mr_Spade）：
可以发现一个无向图就是由若干个无向连通图组成的，因此我们可以枚举组成这个无向图的连通块个数计算贡献。假设我们当前枚举的是$k$个，它的贡献就为：
$\large \frac{f^k}{k!}$
$f^k$很好理解，为什么要除以$k!$呢？因为我们计算的$g_i$是要除以$i!$的，相当于要消除标号的影响，而$f^k$所乘的$k$个元素是有先后顺序之分的，因此需要除以全排列的数目$k!$。那么$f$和$g$的关系式就为：
$\large g=\sum \limits_{k=0}^{\infty}\frac{f^k}{k!}$
这个式子是不是有点眼熟呢？没错！右边那项就是$e^f$的麦克劳林级数。
其实我觉得这个题解的解释就挺好，我再解释一下。首先，指数生成函数相乘，系数就是多重集合的排列。而$f^k$就是把n拆成k个集合的情况。但是我们的每个多项式是不同的，所以每一项都多算了$k!$次。
如果你还是不理解，可以尝试一种暴力理解方式：
~~考虑一个整数划分会被多项式乘出来多少次？是**多重集合排列**次，设我们关注第n项，即为$\frac{n!}{a!b!...}$，好，我们先无脑除以前面的式子，现在的系数是$\frac{a!b!...}{n!}$，但是，这样就去重完了吗？考虑$4=2+2$这种拆分，我们在多重集合排列时认为这两个2是不同的，但其实它们是相同的，所以要除以$2!$。回到原来的式子，我们要除以$a!b!...$，鼓捣一番后发现，真的只剩下分母的$n!$了。~~
不要看上面的暴力解释了，我也不知道我在说啥。
我们设f是不一定联通的方案数，g是联通的方案数，那么有
$$f(n)=\sum_{i=1}^n {n-1 \choose i-1}f(n-i)g(i)$$
移项之后为
$$\frac{f(n)}{(n-1)!}=\sum_{i=1}^n \frac{g(i)}{(i-1)!}\frac{f(n-i)}{(n-i)!}$$
也就是
$$F'=G'F$$
~~接下来开始伪证~~
若$F=e^G$，则$F'=G'F$
所以得证。
----------------
## 多项式三角函数
目前来看没啥用，所以先mark着。[by Semiwaker](https://blog.csdn.net/semiwaker/article/details/73251486)
求$sin A(x),cos A(x)$
对于NTT也可以的！[Link](https://www.luogu.org/blog/user7035/duo-xiang-shi-zong-jie)
[应用](https://en.wikipedia.org/wiki/Alternating_permutation)
---------------------
## 多项式多点求值
~~hu~调了半天发现竟然是对负数取模出错了...~~
我们的老朋友，分治。我们考虑将问题分治成两个小问题。我们就简单粗暴的把询问分成前一半和后一半，如果这个时候，多项式的系数也能减半，就好了。我们构造多项式
$\begin{eqnarray*} P^{[0]}(x) &=& \prod_{i=0}^{\lfloor \frac{n}{2} \rfloor} (x-x_i) \end{eqnarray*}$
这样的多项式，在我们要求的地方值为0，所以不会影响最终答案。我们可以任意减去它——也就是多项式取模！
```cpp
const int rt=1;
#define ls (x<<1)
#define rs (x<<1|1)
#define gm int mid((l+r)/2)
int BUK[N*50],*buk=BUK;
int f[N],a[N];
int *pol[N];
int *init(int sz)
{
int *res=buk; buk+=sz;
return res;
}
void build(int x,int l,int r)
{
if(l==r)
{
pol[x]=init(2);
pol[x][0]=md-a[l],pol[x][1]=1;
return;
}
gm;
build(ls,l,mid); build(rs,mid+1,r);
pol[x]=init(r-l+2);
Mul(mid-l+2,pol[ls],pol[rs],pol[x]);
}
void del(int sz)
{
while(sz--) *--buk=0;
}
int ans[N];
void solve(int x,int l,int r,int *fuc)
{
if(r-l<=100)
{
for(ri i=l; i<=r; ++i)
{
int &s=ans[i];
for(ri j=r-l; j>=0; --j)
s=((LL)s*a[i]+fuc[j])%md;
}
return;
}
gm;
int *p=init(mid-l+1);
Div(r-l+1,fuc,mid-l+2,pol[ls],p);
solve(ls,l,mid,p);
Div(r-l+1,fuc,r-mid+1,pol[rs],p);
solve(rs,mid+1,r,p);
del(r-mid);
}
int n,m;
signed main()
{
#ifdef M207
freopen("in.in","r",stdin);
// freopen("out.out","w",stdout);
#endif
in(n,m); ++n;
for(ri i=0; i<n; ++i) in(f[i]);
for(ri i=1; i<=m; ++i) in(a[i]);
int m2=1; for(; m2<m; m2<<=1);
swap(m,m2);
build(rt,1,m);
if(n>m) Div(n,f,m+1,pol[rt],f);
solve(rt,1,m,f);
for(ri i=1; i<=m2; ++i) out(ans[i]);
return 0;
}
```
----------------------
## 多项式快速插值
[by xyz32768](https://blog.csdn.net/xyz32768/article/details/82832467)
```cpp
int n;
int X[N],Y[N];
int *tf[N];
void solve(int x,int l,int r)
{
if(l==r)
{
tf[x]=init(1);
tf[x][0]=(LL)Y[l]*qpow(gv[l],md-2)%md;
return;
}
gm;
solve(ls,l,mid); solve(rs,mid+1,r);
tf[x]=init(r-l+1);
static int tmp[N];
int t=Mul(mid-l+1,tf[ls],r-mid+1,pol[rs],tmp);
for(ri i=0; i<t; ++i) tf[x][i]=tmp[i],tmp[i]=0;
t=Mul(r-mid,tf[rs],mid-l+2,pol[ls],tmp);
for(ri i=0; i<t; ++i) (tf[x][i]+=tmp[i])%=md,tmp[i]=0;
}
int dg[N];
signed main()
{
#ifdef M207
freopen("in.in","r",stdin);
// freopen("out.out","w",stdout);
#endif
in(n);
for(ri i=1; i<=n; ++i)
{
in(X[i],Y[i]);
a[i]=X[i];
}
build(rt,1,n);
for(ri i=1; i<=n; ++i) dg[i-1]=(LL)pol[rt][i]*i%md;
calc(rt,1,n,dg);
solve(rt,1,n);
print(tf[rt],n);
return 0;
}
```
~~实用的~~多项式科技基本完结了吧...