多项式基础操作全家桶
多项式入门教程
基础概念
形如 \(F(x)=a_1+a_2x+a_3x^2+...+a_nx^n\) 的东西叫多项式。
然后，这个 \(n\) 可以是无穷大的。
其中上面的 \(a_i\) 称为多项式的第 \(i\) 项系数。
多项式乘法
即各项相乘，假设有两个多项式 \(F(x)=\sum_{i=0}^na_ix^i\) 和 \(G(x)=\sum_{i=0}^nb_ix^i\)
那么：
然后我们发现这个式子里有一个部分很有意思，我们专门拿出来：
设两多项式的乘积的第 \(i\) 项系数为 \(c_i\) 那么有：
我们称这个为加法卷积。如果把那个加号换成别的就变成其他的什么卷积了。
那么我们如何计算多项式乘法？你发现可以容易的在 \(O(n^2)\) 内计算。但是这是不优的，我们接下来介绍一种 \(O(n\log n)\) 内计算多项式乘法的算法，先了解一些前置知识。
复数
什么叫复数？形如 \(a+bi\) 的叫复数。其中 \(i^2=-1\)，\(a\) 称为一个复数的实部，\(b\) 称为一个复数的虚部。\(|a+bi|\) 表示其模长。其共轭复数为 \(a-bi\)，对于复数 \(z\) ，其共轭复数记作 \(\overline{z}\)。
然后其四则运算为：
然后模长 \(|z|=|a+b_i|=\sqrt{z\overline{z} }=\sqrt{a^2+b^2}\)
\(a+bi\not=0\) 时存在逆元 \(\frac{1}{a+bi}=\frac{a-bi}{|a+b_i|}\)
复数几何意义
我们把 \(x\) 坐标表示实部， \(y\) 坐标表示虚部的平面称为复平面。可以发现是数轴的扩展。
可以把复数看做一个在复平面上的向量 \((a,b)\)，记其幅角为 \(\theta\) ，模长为 \(r\) 那么有 $a=r\cos{\theta} $ 且 \(b=r\sin{\theta}\)，那么其实复数还可表示为 \(r(\cos{\theta}+i\sin{\theta})\) 。
然后你会惊奇的发现复数相乘相当于模长相乘，幅角相加。
单位根
我们称满足 \(z^n=1\) 的复数 \(z\) 为 \(n\) 次单位根。
首先由于模长相乘，幅角相加的乘法，我们知道 \(n\) 次单位根模长为 \(1\)，\(n\theta=2k\pi\)。因此不同的 \(n\) 次单位根总共有 \(n\) 个，它们将单位圆 \(n\) 等分。
一般我们用 \(w^i_n\) 表示第 \(i+1\) 个 \(n\) 次单位根，即 \(\theta_n^i=\frac{2i\pi}{n}\) 。
一些性质
性质1
性质2
角度所占比相同，模长都为 \(1\) ，因此固然相同。
性质3
单位根可以看做 \(\frac{1}{i}\) 圆周。
性质4
旋转了一个半圆，因此相反。
性质5
证明：
当 \(a=0\) 时，等价于 \(n\) 个第 \(1\) 个 \(n\) 次单位根相加，就是 \(n\) 个 \(1\) 相加，于是得到 \(n\)。
当 \(a\not=0\)，
先了解一下等比数列求和公式
对于 \(a,aq,aq^2,...,aq^n\) 求 \(S_n=a+aq+aq^2+...+aq^n\)
其等价于 \(af(q),f(x)=1+x+x^2+...+x^n\)
我们可以知道 \(f(x)=xf(x)+1-x^{n+1}\) 即 \(f(x)=\frac{1-x^{n+1} }{1-x}\)
那么 \(S_n=af(q)=a\frac{1-q^{n+1} }{1-q}\)
好，那么把 \(w_n^0\) 代入 \(a\) ，\(w_n^a\) 代入 \(q\) 得到
Q.E.D.
离散傅里叶变换
简称 DFT (Discrete Fourier Transform)。
其逆变换称为 IDFT 。
所以实际上
我们知道， \(n\) 个方程可以确定 \(n\) 个未知数。也就是说对于一个 \(n\) 次多项式，我们只需要 \(n+1\) 个点值就可以求出多项式的系数。点值是什么？将 \(x\) 的取值代入多项式得到的值。
然后经过简单验证可以发现，\(F,G\) 在 \(x\) 处的点值的乘积，等于 \(F\times G\) 在 \(x\) 处的点值。
因为我们发现多项式在点值上相乘更快，所以我们考虑一种系数表示转点值表示再转回系数表示的做法。
而这个东西的作用和点值表示有一些相同的性质，以及一个更为优秀的性质。因为傅里叶变换的本质是把单位根代入多项式，DFT 本质是一种特殊的 插值，而 IDFT 也是一种特殊的 多点求值。
回到正题，傅里叶变换。
假设对于一个多项式 \(f\) 我们知道其在 \(0\) 到 \(n-1\) 处的系数。
那么定义其傅里叶变换 \(F[f]\) 满足：（两式的单位根的次方的正负号可以交换）
先证明一下这两个互为逆运算，即证明：
证明：
怎么证明？全部展开！
同理另一个得证。
Q.E.D.
然后实际上，你发现：
卷积定理
然后有：
证明：
也就是傅里叶变换的先后是无关的，所以可以用傅里叶变换的乘积进行逆傅里叶变换得到乘积。
快速傅里叶变换
简称 FFT (Fast Fourier Transform)，用来快速计算傅里叶变换。
为方便处理，我们把 \(n\) 取到大于等于自身的最小的 \(2\) 的幂次。（不用担心会计算错误，因为不存在的高位系数都是 \(0\)）。
记 \(m=\frac{n}{2}\)
设：
然后有 \(F[f](k)=F_0[f](k)+F_1[f](k)\)
观察一下可以发现，\(F_0[f](k+m)=F_0[f](k)\) 且 \(F_1[f](k+m)=w_n^{m}F_1[f](k)\) 且 \(w^{k+m}_n=-w^{k}_n\) 。（三个结论暴力拆解易证）
所以 \(F_0\) 和 \(F_1\) 都只计算前 \(m\) 项，这个分治处理，边界是 \(m-1=1\)，此时两式都是原多项式。
然后复杂度 \(T(n)=2T(\frac{n}{2})+O(n)\) 于是总复杂度 \(O(n\log n)\)
然后再观察一下 \(F\) 和 \(F^-1\):
如果令 \(w'_n=w^{-1}_n\) ：
发现和 \(F\) 长得其实巨像无比，所以我们只需要将 \(w'_n\) 作为单位根进行傅里叶变换就可以实现逆傅里叶变换。
代码实现
由于考虑到代码常数的问题，我们就不考虑实现递归版了。因为这个理解和实现起来都简单。这个递归版的实现可以看 这里
为了保证一个能够过题的常数，我们考虑以下的一些优化。
别的没什么好说的，我们重点说一下 FFT 的二进制翻转优化。
先给出结论，递归到最底层时，位置为 \(i\) 处理的是在原序列位置为 \(rev(i)\) 的地方，其中 \(rev(i)\) 表示 \(i\) 的二进制翻转。
我们直接把序列放到这个位置便是对的，因为我们只在乎递归边界的值的大小，中间的位置我们不需要在意。
那么为什么递归的最底层满足位置的二进制翻转？
证明：
考虑归纳法证明。
当 \(n=1\)，有 \(rev(0)=0\) 成立。
当 \(n=2\)，有 \(rev(0)=0,rev(1)=1\) 成立。
当 \(n=4\)， 有 \(rev(0)=0,rev(1)=2,rev(2)=1,rev(3)=3\) 成立。
当 \(n=8\)，有：
0 1 2 3 4 5 6 7 第1层
0 2 4 6|1 3 5 7 第2层
0 4|2 6|1 5|3 7 第3层
0|4|2|6|1|5|3|7 第4层
成立。
假设结论对于 \(n=2^k\) 成立，证明对于 \(n=2^{k+1}\) 也成立。
有
左边的位置全部除去二进制下最低位的 \(0\) 得到
我们知道，对于 \(n=2^k\) 结论成立，于是左边重新加回最低位 \(0\) ，发现此操作前后二进制翻转不变，因此左半边结论成立。
右半边除去最低位 \(1\) 得到相同的式子。
由于对于 \(n=2^k\) 成立，右边加回最低位 \(1\) ，发现操作前后二进制翻转相差 \(2^k\)，减去左半边下标 \(2^k\) 正好对应右半边的二进制翻转。于是右半边结论成立。
综上结论成立。Q.E.D.
那么我们直接把这位的多项式的值换到二进制翻转的位置计算即可。至于怎么线性求二进制翻转？看 这里
可能你有点疑惑，为什么递归最底层的值不需要操作可以直接用？你代入 \(F[f](0)\) 你发现 \(w^{0}=1\) ，直接就是原本的多项式。
#include<bits/stdc++.h>
#define ll long long
#define db double
#define file(a) freopen(#a".in","r",stdin),freopen(#a".out","w",stdout)
#define Better_IO true
namespace IO{
#if Better_IO
char buf[(1<<20)+3],*p1(buf),*p2(buf);
const int lim=1<<20;
inline char gc(){
if(p1==p2) p2=(p1=buf)+fread(buf,1,lim,stdin);
return p1==p2?EOF:*p1++;
}
#define pc putchar
#else
#define gc getchar
#define pc putchar
#endif
inline bool blank(char c){
return c=='\t' or c=='\n' or c=='\r' or c==' ' or c==EOF;
}
inline void gs(char *s){
char ch=gc();
while(blank(ch) ) ch=gc();
while(!blank(ch) ) {*s++=ch;ch=gc();}
*s=0;
}
inline void ps(char *s){
while(*s!=0) {pc(*s++);}
}
template<class T>
inline void read(T &s){
s=0;char ch=gc();bool f=0;
while(ch<'0'||'9'<ch) {if(ch=='-') f=1;ch=gc();}
while('0'<=ch&&ch<='9') {s=s*10+(ch^48);ch=gc();}
if(ch=='.'){
db p=0.1;ch=gc();
while('0'<=ch&&ch<='9') {s=s+p*(ch^48);ch=gc();p*=0.1;}
}
s=f?-s:s;
}
template<class T,class ...A>
inline void read(T &s,A &...a){
read(s);read(a...);
}
};
using IO::gs;
using IO::ps;
using IO::read;
using IO::gc;
const int N=(1<<20)+3;
const db pi=acos(-1);
struct Complex{
db a,b;
Complex(){};
Complex(db _a,db _b){
a=_a;b=_b;
}
inline friend Complex operator + (Complex x,Complex y){
return {x.a+y.a,x.b+y.b};
}
inline friend Complex operator - (Complex x,Complex y){
return {x.a-y.a,x.b-y.b};
}
inline friend Complex operator * (Complex x,Complex y){
return {x.a*y.a-x.b*y.b,x.a*y.b+x.b*y.a};
}
};
int n,m;
Complex F[N<<1],G[N<<1],w[N];
int rev[N<<1];
inline void getrev(){
for(int i=0;i<n;++i){
rev[i]=(rev[i>>1]>>1)|((i&1)?(n>>1):0);
}
}
inline void FFT(Complex *f,bool I){
for(int i=0;i<n;++i){
if(rev[i]<i)
std::swap(f[i],f[rev[i] ]);
}
for(int p=2;p<=n;p<<=1){
int len=p>>1;
w[0]=Complex(1,0);w[1]=Complex(cos(2*pi/p),sin(2*pi/p) );
if(I) w[1].b*=-1;
for(int i=2;i<len;++i){
w[i]=w[i-1]*w[1];
}
for(int k=0;k<n;k+=p){
for(int l=k;l<k+len;++l){
Complex tmp=f[l+len]*w[l-k];
f[l+len]=f[l]-tmp;
f[l]=f[l]+tmp;
}
}
}
}
int main(){
file(a);
read(n,m);
for(int i=0;i<=n;++i){
read(F[i].a);
}
for(int i=0;i<=m;++i){
read(G[i].a);
}
for(m+=n,n=1;n<=m;n<<=1);
getrev();
FFT(F,0);FFT(G,0);
for(int i=0;i<n;++i) F[i]=F[i]*G[i];
FFT(F,1);
for(int i=0;i<=m;++i){
printf("%d ",(int)(F[i].a/n+0.5) );
}
return 0;
}
快速数论变换
简称 NTT (Number-theoretic transform)。
很多时候计数问题都是在取模意义下的，而 FFT 使用浮点数，不单单是巨慢无比，更可怕的是基本上用不了了，怎么办？于是考虑选择一个模意义下的玩意来代替单位根。最终大佬们找到了 原根。然后就形成了 NTT。
原根与阶
阶
定义：\(a\) 在模 \(p\) 意义下的阶为一个最小的整数 \(k\) 使得 \(a^k\equiv 1\pmod p\)。其中 \(gcd(a,p)=1\)。记作 \(\delta_p(a)\)。
本质上是幂次的最小循环节。
性质1
根据欧拉定理，我们发现一个数 \(a\) 的阶会小于等于 \(\varphi(a)\)。
性质2
若 \(a^n\equiv 1\pmod m\) , 则 \(\delta_m(a)|n\)
证明：
令 \(n=\delta_m(a)q+r\quad (0\le r<\delta_m(a) )\)
当 \(r>0\) 有 \(a^r\equiv a^r(a^{\delta_m(a)})^q\equiv 1\pmod m\)
然后就是有一个比 \(\delta_m(a)\) 更小的满足阶的性质的东西 \(r\) ，于是矛盾。
于是 \(r=0\) ，即 \(\delta_m(a)|n\) ，得证，Q.E.D.
性质3
对于 \(a,b,m\) ，且 \(gcd(a,m)=gcd(b,m)=1\) 则：
\(\delta_m(ab)=\delta_m(a)\delta_m(b)\) 的充分必要条件是 \(gcd(\delta_m(a),\delta_m(b) )=1\)
证明：
必要性证明：
由阶的定义可得
由 性质2 可知
由 \(\delta_m(ab)=\delta_m(a)\delta_m(b)\) 得：
由公式 \(ab=lcm(a,b)\gcd(a,b)\) 可知：
\(gcd(\delta_m(a),\delta_m(b) )=1\)
Q.E.D.
充分性证明：
可得 \(\delta(a)|\delta(ab)\delta_m(b)\) 再结合 \(gcd(\delta_m(a),\delta_m(b) )=1\) 得到 \(\delta_m(a)|\delta_m(ab)\)
同理可得： \(\delta_m(b)|\delta_m(ab)\)
因此 \(\delta_m(a)\delta_m(b)|\delta_m(ab)\)
然后类似地可以得到 \(\delta_m(ab)|\delta_m(a)\delta_m(b)\)
性质4
对于 \(\gcd(a,m)=1\)，有:
证明：
由：
可得 \(\delta_m(a)|k\delta_m(a^k)\) ，则 \(\frac{\delta_m(a) }{\gcd(\delta_m(a),k)}|\delta_m(a^k)\)
再由：
于是有 \(\delta_m(a^k)|\frac{\delta_m(a) }{\gcd(\delta_m(a),k)}\)
然后把刚才这个式子与上面得到的式子结合，可得：
Q.E.D.
原根
定义：一个数 \(a\) 在模 \(p\) 意义下满足 \(\delta_p(a)=\varphi(p)\) ，则称之为 \(p\) 的原根。
证明 性质1 前，需要证明两个定理：
拉格朗日定理
对于素数 \(p\)，对于模 \(p\) 意义下的整系数多项式 \(f(x)\):
其同余方程 \(f(x)\equiv 0 \pmod p\) 在模 \(p\) 意义下至多有 \(n\) 个不同解。
证明：
当 \(n=0\) 时，有 \(p\not|a_0\) ，因此 \(f(x)\equiv 0\pmod p\) 是无解的，定理对于 \(n=0\) 成立。
假设定理对于 \(\deg f<n\) 成立，考虑证明定理对于 \(\deg f=n\) 成立。
反证法，假设存在 \(n+1\) 个符合条件的 \(x_0,x_1,...,x_n\) 。
现在设 \(g(x)\) 满足 \(f(x)-f(x_0)=(x-x_0)g(x)\) ，因此可知 \(g(x)\) 在模 \(p\) 意义下最多是 \(n-1\) 次多项式，又因为对于 \(i\in[0,n]\) 都满足 \(f(x_i)\equiv 0 \pmod p\) ，于是有：
对于 \(i\in[1,n]\) 有，
又因为 \(x_i\not \equiv x_0\pmod p\) ，得到 \(g(x)\equiv 0 \pmod p\) ，此时 \(g(x)\equiv 0\pmod p\) 这个同余方程有 \(n\) 个根，与归纳法的假设冲突，矛盾。
因此定理对 \(\deg f=n\) 也成立，得证。
Q.E.D.
定理1
对于 \(a,b,p\) 满足 \(gcd(a,p)=gcd(b,p)=1\) ,存在 \(c\in Z\) 使得 \(\delta_p(c)=lcm(\delta_p(a),\delta_p(b))\)
证明：
对于 \(\delta_p(a),\delta_p(b)\) 质因数分解，得到：
令 \(k_1\) 为所有 \(n_i\le m_i\) 的 \(p^{n_i}\) 的乘积，\(k_2\) 为所有 \(n_i>m_i\) 的 \(p_i^{m_i}\) 的乘积。
然后取 \(x,y\) 使得 \(\delta_p(a)=k_1x,\delta_P(b)=k_2y\)。可以发现 \(\gcd(x,y)=1\) 且 \(lcm(\delta_p(a),\delta_p(b) )=xy\) 。
然后由 阶的性质4 ，\(\delta_p(a^{k_1})=\frac{\delta_p(a)}{\gcd(\delta_p(a),k_1)}=x,\delta_p(b^{k_2})=y\) 。
再由 阶的性质3 可得：（条件：\(\gcd(\delta_p(a),\delta_p(b))=1\) 且 \(\gcd(a,p)=\gcd(b,p)=1\)）
可得 \(c=a^{k_1}b^{k_2}\) 。
Q.E.D.
性质1
对于任意奇素数（非 \(2\) 素数）\(p\) ，其存在原根。
证明：
由上面的 定理1 ，可以发现存在一个 \(g\) 使得 \(\delta_p(g)=lcm(\delta_p(1),...,\delta_p(p-1) )\)
于是，对于任意 \(i\in[1,p-1]\) 有 \(\delta_p(i)|\delta_p(g)\) 即 \(i^{\delta_p(g)}\equiv 1\pmod p\)。也就是说，这些 \(i\) 全部都是同余方程 \(x^{\delta_p(g)}\equiv 1\pmod p\) 的解。
有 拉格朗日定理 可知，\(\delta_p(g)\le p-1\) ，同时由费马小定理，有 \(\delta_p(g)\ge p-1\) ，可得 \(\delta_p(g)=p-1=\varphi(p)\) 。
所以 \(g\) 是模 \(p\) 意义下的原根。
Q.E.D.
性质2
对于奇素数 \(p\) 和 \(a\in N^*\) ，\(p^a\) 有原根。
证明：
考虑先证明 \(p^2\) 的原根存在。
记 \(p\) 的一个原根为 \(g\) ，其在模 \(p^2\) 的阶为 \(d\) ，则有 \(d|\varphi(p^2)\) （欧拉定理+阶的性质2），再由 \(g^d\equiv 1\pmod p\) 于是有 \(g^d\equiv 1\pmod{p^2}\)，因此 \(\varphi(p)|d\)。
由于 \(\varphi(p)=(p-1)|d\) 且 \(d|p(p-1)=\varphi(p^2)\) ，\(p\) 为奇素数，可知 \(d=p-1\) 或 \(d=p(p-1)\) 。
当 \(d=p-1\) ，则 \(g\) 为 \(p^2\) 的原根。当 \(d=p(p-1)\) 时，考虑证明 \((g+p)\) 是 \(p^2\) 的原根。
首先 \((g+p)\) 为 \(p\) 的原根。同上可知 \(g+p\) 在 \(p^2\) 的阶为 \((p-1)\) 或 \(p(p-1)\) 。由 \(g^{p-1}\equiv 1\pmod{p^2}\) 结合二项式定理可得：
可知 \((g+p)^{p-1}\not \equiv 1\pmod{p^2}\) ,则 \((g+p)^{p(p-1)}\equiv 1\pmod{p^2}\)。
然后可以归纳证明：
假设 \(g\) 是 \(p^a\) 的原根 ，\(d\) 是 \(g\) 对于模 \(p^{a+1}\) 的阶，考虑证明 \(g\) 是 \(p^{a+1}\) 。于是类似的可以得到 \(d=p^{a-1}(p-1)\) 或 \(d=p^a(p-1)\)。
由 欧拉定理: \(g^{p^{a-2}(p-1) }\equiv 1\pmod {p^{a-1}}\) ，即 \(g^{p^{a-2}(p-1) }=1+kp^{a-1}\)
由于 \(g\) 为 \(p^a\) 的原根,\(d>p^{a-2}(p-1)\)，有 \(g^{p^{a-2}(p-1)}\not \equiv 1 \pmod{p^a}\)
于是又有 \(p^a\not |kp^{a-1}\) 即 \(\gcd(k,p)=1\)
然后二项式定理展开，有：
于是 \(d=p^a(p-1)\) , \(p^{a+1}\) 存在原根 \(g\)。
Q.E.D.
性质3
对于奇素数 \(p\) 和 \(a\in N^*\) ，\(2p^a\) 有原根。
证明：
假设对于 \(p^a\) 存在原根 \(g\)。首先我们知道对于 \(2p^a\) ，它的原根必然是奇数的。而对于 \(p^a\) 的原根 \(g\) ，以及 \((g+p^a)\) ，其中至少有一个是奇数，我们记这个奇数为 \(b\)。
设 \(b\) 在模 \(2p^a\) 的意义下的阶为 \(d\) ，于是有 \(d|\varphi(2p^a)=\varphi(p^a)\)（阶的性质2）且 \(\varphi(p^a)|d\)（阶的性质2），即 \(d=\varphi(p^a)=\varphi(2p^a)\)。
所以 \(b\) 为模 \(2p^a\) 的原根。
Q.E.D.
性质4
对于 \(p\not = 2,4,p^c,2p^c\) ，其不存在原根。
证明：
分三种情况来讨论：
1.对于 \(p=2^a,a\ge 3\)
2.对于 \(p\) 为 \(2^a\) 与一个奇素数的乘积且 \(a\ge 2\)
3.对于 \(p\) 为 \(2^a\) 与多个奇素数的乘积且 \(a\ge 1\)
我们先处理第 \(2,3\) 种情况。这两种情况其实可以概括成：
对于 \(p=rs\) ，有 \(r,s>2\) 且 \(\gcd(r,s)=1\)：
由 \(\gcd(r,s)=1\) 可知 \(\varphi(p)=\varphi(r)\varphi(s)\) ，然后由 \(r,s>2\) 得到 \(4|\varphi(p)\) （因为此时两者的欧拉函数都为偶数）。
记 \(d=\frac{1}{2}\varphi(p)\) ，那么 \(\varphi(r),\varphi(s)\) 都是 \(d\) 的因子。于是对于任意 \(\gcd(a,p)=1\)，有 \(a^d\equiv 1\pmod r\) 和 \(a^d\equiv 1\pmod s\) ，则有 \(a^d\equiv 1\pmod p\) 。发现 \(d<\varphi(p)\)，于是没有原根，
而第 \(1\) 种情况：
由于互质的限制，我们只需考虑奇数。
因此对于任意奇数 \(b=2k+1\) ，有：
解释一下：
其中第二步，由于从第 \(4\) 项开始，都具有 \(2^{a-1}(2^{a-2}-2)\) 的组合，具有 \(2^a\) 的因子，于是可以删去。
其中第三步 \(k+(2^{a-2}-1)k^2\) 是偶数。（因为两项奇偶性相同）
于是 \(a^{2^{a-1}}\not \equiv 1\pmod{2^a}\)，没有原根。
Q.E.D.
性质5:原根存在定理
形如 \(1,2,4,p^c,2p^c\) （\(p\) 为奇质数）的数才存在原根。
证明：
对于 \(p=1,2,4\) 的情况，分别有原根 \(1,1,3\)。
对于 \(p^c\) 我们在 性质2 中已经证明，对于 \(2p^c\) 我们在 性质3 中得到证明。
对于 \(p\not=2,4,p^c,2p^c\) 的，我们在 性质4 中证明。
综上，该定理成立。
Q.E.D.
性质6:原根的数量级
可以证明，如果 \(m\) 有原根，那么其最小原根 \(g\le m^{0.25}\) ，证明略（我不道啊）。
所以找最小原根是可以暴力找的。怎么找呢？
性质7:原根判定定理
对于 \(m\ge 3,\gcd(a,m)=1\) ，则 \(a\) 是模 \(m\) 意义下的原根的充要条件是：对于 \(\varphi(m)\) 的每个质因数 \(p\)，都有 \(a^{\frac{\varphi(m)}{p}}\not \equiv 1\pmod m\)。
证明：
考虑反证，假设存在一个非模 \(m\) 的原根 \(a\) 使得存在一个质因数 \(p\) 满足 \(a^{\frac{\varphi(m)}{p} }\not\equiv 1\pmod m\) 。
于是存在 \(t<\varphi(m)\) 使得 \(a^t\equiv 1\pmod m\)。
由 裴蜀定理 可知，存在 \(x,y\) 使得 \(ty=\varphi(m)x+\gcd(\varphi(m),t)\)，于是有 \(a^{ty}\equiv 1\equiv a^{\varphi(m)x+\gcd(\varphi(m),t)} \equiv a^{\gcd(t,\varphi(m))} \pmod m\) ，中间的转换根据欧拉定理。
则 \(\gcd(t,\varphi(m))|\varphi(m)\) ，且 \(\gcd(t,\varphi(m)) \le t<\varphi(m)\)。
然后可知存在 \(\varphi(m)\) 的质因数 \(p\) 使得 \(\gcd(t,\varphi(m) )|\frac{\varphi(m)}{p}\) 于是有 \(a^{\frac{\varphi(m)}{p}}\equiv 1\pmod m\) ，矛盾。
Q.E.D.
求出原根
有了上面的性质的支撑，我们就可以计算原根了。我们直接暴力枚举然后用原根判定定理检验，其他原根可由这个原根推出。具体来说，\(g^i,i\in[0,\varphi(p)]\) 在模 \(p\) 意义下互不相同，我们只需要检验其阶是否为 \(\varphi(p)\) 即可（凭借 阶的性质4 可以转化为 \(i\) 与 \(\varphi(p)\) 是否互质，其揭露了原根循环群的本质）。
复杂度 \(O(p^{0.25}\varphi(p) \log \varphi(p)+\sqrt p)\) 。
给出参考代码：
#include<bits/stdc++.h>
#define ll long long
#define db double
#define file(a) freopen(#a".in","r",stdin),freopen(#a".out","w",stdout)
#define sky fflush(stdout)
#define Better_IO true
namespace IO{
#if Better_IO==true
char buf[(1<<20)+3],*p1(buf),*p2(buf);
const int lim=1<<20;
inline char gc(){
if(p1==p2) p2=(p1=buf)+fread(buf,1,lim,stdin);
return p1==p2?EOF:*p1++;
}
#define pc putchar
#else
#define gc getchar
#define pc putchar
#endif
inline bool blank(const char &c){
return c==' ' or c=='\n' or c=='\t' or c=='\r' or c==EOF;
}
inline void gs(char *s){
char ch=gc();
while(blank(ch) ) {ch=gc();}
while(!blank(ch) ) {*s++=ch;ch=gc();}
*s=0;
}
inline void gs(std::string &s){
s="";char ch=gc();s+='#';
while(blank(ch) ) {ch=gc();}
while(!blank(ch) ) {s+=ch;ch=gc();}
}
inline void ps(char *s){
while(*s!=0) pc(*s++);
}
inline void ps(const std::string &s){
for(auto it:s)
if(it!='#') pc(it);
}
template<class T>
inline void read(T &s){
s=0;char ch=gc();bool f=0;
while(ch<'0'||'9'<ch) {if(ch=='-') f=1;ch=gc();}
while('0'<=ch&&ch<='9') {s=s*10+(ch^48);ch=gc();}
if(ch=='.'){
db p=0.1;ch=gc();
while('0'<=ch&&ch<='9') {s=s+p*(ch^48);p*=0.1;ch=gc();}
}
s=f?-s:s;
}
template<class T,class ...A>
inline void read(T &s,A &...a){
read(s);read(a...);
}
};
using IO::read;
using IO::gs;
using IO::ps;
const int N=2e6+3;
int n,d;
inline int inc(int x,int y){
return (x+=y)>=n?x-n:x;
}
inline int dec(int x,int y){
return (x-=y)<0?x+n:x;
}
inline int mul(int x,int y){
return 1ll*x*y%n;
}
inline int qpow(int a,int b){
int ans=1;
while(b){
if(b&1) ans=mul(ans,a);
a=mul(a,a);
b>>=1;
}
return ans;
}
std::vector<int>pri;
db phi;
inline void get_phi(int x){
phi=x;
for(auto y:pri){
phi*=(y-1);phi/=y;
}
}
inline bool check(int x){
for(auto y:pri){
if(qpow(x,phi/y)==1)
return 0;
}
return 1;
}
int gcd(int x,int y){
if(y==0) return x;
return gcd(y,x%y);
}
inline void answer(int x){
static int ans[N],top;
top=0;
for(int i=0;i<phi;++i){
if(gcd(i,phi)==1){
ans[++top]=qpow(x,i);
}
}
std::sort(ans+1,ans+1+top);
printf("%d\n",top);
for(int i=d;i<=top;i+=d){
printf("%d ",ans[i]);
}
}
int main(){
file(a);
int T;read(T);
while(T--){
read(n,d);
if(n==2){
printf("1\n");
if(d==1) printf("1\n");
else printf("\n");
continue;
}
int x=n;
pri.clear();
for(int i=2;i*i<=x;++i){
if(x%i==0){
while(x%i==0) x/=i;
pri.push_back(i);
}
}
if(x!=1){
pri.push_back(x);
}
get_phi(n);
x=phi;
pri.clear();
for(int i=2;i*i<=x;++i){
if(x%i==0){
while(x%i==0) x/=i;
pri.push_back(i);
}
}
if(x!=1){
pri.push_back(x);
}
bool f=0;
for(int i=1;i<=n;++i){
if(check(i) and gcd(i,n)==1){
answer(i);f=1;
break;
}
}
if(!f) printf("0\n");
pc('\n');
}
return 0;
}
再论快速数论变换
回到 NTT 本身，考虑使用原根来代替单位根。
首先，我们看看在 FFT 中使用了那些单位根的性质：
1.\(\omega_n^k=(\omega_n^1)^k\)
因此只需求出一个单位根 \(\omega_n^1\)。
2.\(\omega_n^i,i\in[0,n)\) 互不相同
保证可以得到 \(n\) 个不同点值。
3.\(\omega_n^k=\omega_n^{k\bmod n}\)
有了这条就然后需要推出：
展现出一个循环，或者说圆。
4.\(\omega_{na}^{ka}=\omega_n^k\)
等分性质。
我们看看怎么使得原根满足这几个性质。
其实是考虑一个模意义下阶为 \(n\) 的一个元素 \(g\)。
发现对于模素数 \(p\)，其原根 \(g\) 的阶为 \(p-1\)。
于是有 \(\delta_p(g^k)=\frac{\delta_m(g)}{\gcd(k,\delta_p(g) )}=\frac{p-1}{\gcd(k,p-1)}\)
我们发现，当 \(n|(p-1)\) 时， \(\delta_p(g^{\frac{p-1}{n} })=n\) ，此时 \(g^{\frac{p-1}{n}}\) 满足条件 1,2,3。
至于为什么满足 3 的推理：
由 \((g^{\frac{n}{2} })^2=g^n\equiv 1\pmod p\)
由 \(g^i,i\in[1,n]\) 互不相同，得到 \(g^{\frac{n}{2}}\equiv -1\pmod p\)
得 \(-g^{k+\frac{n}{2} }\equiv g^k\pmod p\) 。、
Q.E.D.
至于另一个推理则同理单位根的证明了。
4 也是易证的: \((g^{\frac{p-1}{kn}})^k=g^{\frac{p-1}{n}}\) 。
于是 \(n|(p-1)\) 时，\(g^{\frac{n}{2} }\) 是可以代替单位根存在的。
我们发现 \(p=998244353\) 时， \(p-1=2^{23}\times 7\times 17\)。而我们的 \(n\) 是 \(2\) 的幂次，于是总是整除的。我们可以把这种模数称为 NTT模数。
接下来直接替换即可。顺带一提, \(998244353\) 的一个原根为 \(3\)。
至于 IDFT 是由其逆元来弄的，把原根换成其逆元即可。
代码实现
#include<bits/stdc++.h>
#define ll long long
#define db double
#define file(a) freopen(#a".in","r",stdin),freopen(#a".out","w",stdout)
#define sky fflush(stdout)
#define Better_IO true
namespace IO{
#if Better_IO==true
char buf[(1<<20)+3],*p1(buf),*p2(buf);
const int lim=1<<20;
inline char gc(){
if(p1==p2) p2=(p1=buf)+fread(buf,1,lim,stdin);
return p1==p2?EOF:*p1++;
}
#define pc putchar
#else
#define gc getchar
#define pc putchar
#endif
inline bool blank(const char &c){
return c==' ' or c=='\n' or c=='\t' or c=='\r' or c==EOF;
}
inline void gs(char *s){
char ch=gc();
while(blank(ch) ) {ch=gc();}
while(!blank(ch) ) {*s++=ch;ch=gc();}
*s=0;
}
inline void gs(std::string &s){
s="";char ch=gc();s+='#';
while(blank(ch) ) {ch=gc();}
while(!blank(ch) ) {s+=ch;ch=gc();}
}
inline void ps(char *s){
while(*s!=0) pc(*s++);
}
inline void ps(const std::string &s){
for(auto it:s)
if(it!='#') pc(it);
}
template<class T>
inline void read(T &s){
s=0;char ch=gc();bool f=0;
while(ch<'0'||'9'<ch) {if(ch=='-') f=1;ch=gc();}
while('0'<=ch&&ch<='9') {s=s*10+(ch^48);ch=gc();}
if(ch=='.'){
db p=0.1;ch=gc();
while('0'<=ch&&ch<='9') {s=s+p*(ch^48);p*=0.1;ch=gc();}
}
s=f?-s:s;
}
template<class T,class ...A>
inline void read(T &s,A &...a){
read(s);read(a...);
}
};
using IO::read;
using IO::gs;
using IO::ps;
const int N=(1<<20)+3;
const int mods=998244353;
inline int inc(int x,int y){
return (x+=y)>=mods?x-mods:x;
}
inline int dec(int x,int y){
return (x-=y)<0?x+mods:x;
}
inline int mul(int x,int y){
return 1ll*x*y%mods;
}
inline int qpow(int a,int b){
int res=1;
while(b){
if(b&1) res=mul(res,a);
a=mul(a,a);
b>>=1;
}
return res;
}
int n,m;
int g=3,invg=qpow(g,mods-2);
int F[N<<1],G[N<<1],w[N];
int rev[N<<1];
inline void get_rev(){
for(int i=1;i<n;++i){
rev[i]=(rev[i>>1]>>1)|((i&1)?(n>>1):0);
}
}
inline void NTT(int *f,bool I){
for(int i=0;i<n;++i){
if(i<rev[i]) std::swap(f[i],f[rev[i] ]);
}
for(int p=2;p<=n;p<<=1){
int len=p>>1;
w[0]=1;w[1]=qpow(I?invg:g,(mods-1)/p);
for(int i=2;i<len;++i){
w[i]=mul(w[i-1],w[1]);
}
for(int k=0;k<n;k+=p){
for(int l=k;l<k+len;++l){
int tmp=mul(f[l+len],w[l-k]);
f[l+len]=dec(f[l],tmp);
f[l]=inc(f[l],tmp);
}
}
}
}
int main(){
file(a);
read(n,m);
for(int i=0;i<=n;++i){
read(F[i]);
}
for(int i=0;i<=m;++i){
read(G[i]);
}
for(m+=n,n=1;n<=m;n<<=1);
get_rev();
NTT(F,0);NTT(G,0);
for(int i=0;i<n;++i) F[i]=mul(F[i],G[i]);
NTT(F,1);
int invn=qpow(n,mods-2);
for(int i=0;i<=m;++i){
printf("%d ",mul(F[i],invn) );
}
return 0;
}
多项式插值
拉格朗日插值
对于 \(n\) 次多项式 \(f(x)\) ，我们知道其 \(n+1\) 个点值 \((x_i,y_i)\)，于是有：
由性质 \(f(x)\equiv f(a) \pmod{(x-a)}\) 结合中国剩余定理易得。
计算的复杂度为 \(O(n^2)\)。
如果点值连续，可通过小转换得到：
通过预处理可以 \(O(n)\) 计算。
多项式快速插值
快速沃尔什变换
快速沃尔什变换, 即 FWT(Fast Walsh Transform).
其意义是位运算卷积，也就是将我们熟知的加法卷积的加法换为位运算。即：
其中我们用 \(\oplus\) 来代指某种位运算。
我们期望的做法同样是将多项式转为点值进行运算后在转回来。
考虑构造一组 \(FWT[f(x)]\) 的变化，形如：
其中 \(c(i,j)\) 表示式子的系数。
我们知道， \(A\times B=C\) ，则 \(FWT[A]\times FWT[B]=FWT[C]\) 。
则：
得到 \(c(i,j)c(i,k)=c(i,j,j\oplus k)\) 。然后不难发现 \(c(i,x)\) 可以把 \(x\) 按照当前二进制操作拆分然后相乘。
我们考虑如下的做法：
发现左边没有 \(n\) 的最高位 \(0\) ，而右边有。
记 \(x'\) 表示 \(x\) 在二进制上去除最高位，\(x_i\) 表示 \(x\) 二进制上从高到底的第 \(i\) 位。
得到：
发现一直这么拆是 \(O(\log n)\) ，合并两式复杂度为 \(O(n)\) ，总复杂度 \(O(n\log n)\)。
发现 \(c\) 是个矩阵，我们对其求逆然后再 FWT ，就可以做到 IFWT 。
然后我们要用的 \(c\) 矩阵只有 \(2\times 2\)。
接下来为 \(3\) 种位运算构造矩阵。构造方法是把式子全列出来，然后手动枚举找可行情况，过于繁琐，于是直接给出构造和细节。
Or
或
一般采用第二种，因为其同时具有子集卷积的性质。也是记忆技巧。
其逆矩阵为：
And
或
采用第二种，其逆矩阵为：
记忆技巧，为 \(Or\) 矩阵的转置。
具有超集卷积的性质。
Xor
或
仍然采用第二种。
记忆技巧，\(c(i,j)=(-1)^{i\&j}\)。
其逆矩阵为：
#include<bits/stdc++.h>
#define ll long long
#define db double
#define file(a) freopen(#a".in","r",stdin),freopen(#a".out","w",stdout)
#define sky fflush(stdout)
#define gc getchar
#define pc putchar
namespace IO{
template<class T>
inline void read(T &s){
s=0;char ch=gc();bool f=0;
while(ch<'0'||'9'<ch) {if(ch=='-') f=1;ch=gc();}
while('0'<=ch&&ch<='9') {s=s*10+(ch^48);ch=gc();}
if(ch=='.'){
T p=0.1;ch=gc();
while('0'<=ch&&ch<='9') {s=s+p*(ch^48);p/=10;ch=gc();}
}
s=f?-s:s;
}
template<class T,class ...A>
inline void read(T &s,A &...a){
read(s);read(a...);
}
};
using IO::read;
const int mod=998244353;
inline int inc(int x,int y){
return (x+=y)>=mod?x-mod:x;
}
inline int dec(int x,int y){
return (x-=y)<0?x+mod:x;
}
inline int mul(int x,int y){
return 1ll*x*y%mod;
}
inline int qpow(int a,int b){
int res=1;
while(b){
if(b&1) res=mul(res,a);
a=mul(a,a);
b>>=1;
}
return res;
}
const int N=(1<<17)+3;
int _n,n;
int a[N],b[N],ans1[N],ans2[N],ans3[N];
int oa[N],ob[N];
const int inv2=qpow(2,mod-2);
const int
Cor[2][2]={{1,0},{1,1} },ICor[2][2]={{1,0},{mod-1,1} },
Cand[2][2]={{1,1},{0,1} },ICand[2][2]={{1,mod-1},{0,1} },
Cxor[2][2]={{1,1},{1,mod-1} },ICxor[2][2]={{inv2,inv2},{inv2,dec(mod,inv2)} };
inline void FWT(int *f,const int c[2][2]){
for(int p=2;p<=n;p<<=1){
int len=(p>>1);
for(int k=0;k<n;k+=p){
for(int l=k;l<k+len;++l){
int tl=f[l],tr=f[l+len];
f[l]=inc(mul(c[0][0],tl),mul(c[0][1],tr) );
f[l+len]=inc(mul(c[1][0],tl),mul(c[1][1],tr) );
}
}
}
}
int main(){
file(a);
read(_n);n=(1<<_n);
for(int i=0;i<n;++i){
read(a[i]);
oa[i]=a[i];
}
for(int i=0;i<n;++i){
read(b[i]);
ob[i]=b[i];
}
FWT(a,Cor);FWT(b,Cor);
for(int i=0;i<n;++i) ans1[i]=mul(a[i],b[i]);
FWT(ans1,ICor);
for(int i=0;i<n;++i){
a[i]=oa[i];
b[i]=ob[i];
}
FWT(a,Cand);FWT(b,Cand);
for(int i=0;i<n;++i) ans2[i]=mul(a[i],b[i]);
FWT(ans2,ICand);
for(int i=0;i<n;++i){
a[i]=oa[i];
b[i]=ob[i];
}
FWT(a,Cxor);FWT(b,Cxor);
for(int i=0;i<n;++i) ans3[i]=mul(a[i],b[i]);
FWT(ans3,ICxor);
for(int i=0;i<n;++i){
printf("%d ",ans1[i]);
}pc('\n');
for(int i=0;i<n;++i){
printf("%d ",ans2[i]);
}pc('\n');
for(int i=0;i<n;++i){
printf("%d ",ans3[i]);
}pc('\n');
return 0;
}
FWT本质的探索
把数看成很多 \(01\) 维，那么 \(Or\) 运算等价于每一维度取 \(max\) ， \(And\) 运算等价于每一维取 \(min\)，而 \(Xor\) 就是每维加起来 \(\bmod 2\)。
然后考虑扩展到 \(k\) 进制数，那么 \(Or\) 运算等价于每一维度取 \(max\) ， \(And\) 运算等价于每一维取 \(min\)，而 \(Xor\) 就是每维加起来 \(\bmod k\)。
Or
考虑 \(Or\) 运算，发现矩阵里面只会有 \(0,1\) 。然后考虑同一行，对于 \(y_1<y_2\) 有 \(c(x,y_1)c(x,y_2)=c(x,max\{y_1,y_2\})=c(x,y_2)\) 。也就是说，同一行内的构成必定是先 \(1\) 后 \(0\) 。然后为了保证矩阵可逆，其需要满秩，那么每行互不相同即可。我们采用如下矩阵：（假设 \(k=4\)）
及其逆：
就是每个 \(1\) 下面带个 \(-1\) 。
两个矩阵可以看做先（对各维子集）做前缀和，再做差分。
套用 FWT ，复杂度 \(O(nk\log_k n)\)。
直接计算高维前缀和复杂度 \(O(n\log_k n)\)（每次 \(k\) 进制下把一位挖掉），下同。
And
和 \(Or\) 正好相反，矩阵也恰好是为 \(Or\) 矩阵的转置。
及其逆：
两个矩阵可以看做先（对各维超集）做后缀和，再做差分。
