单位根和单位根反演
本文作者为 JustinRochester。
单位根和单位根反演
引入
根据上一篇的内容，我们了解到复数可以表示为 \(z=re^{i\theta}\) 的形式，其中 \(r,\theta\) 分别表示复数的模长和幅角。
当 \(r>1\) 时，复数 \(z\) 反复乘自己将使得自己的模长越来越大；而当 \(r<1\) 时，将使得模长越来越小。太大或太小的数都会对精度产生一定的影响。
为了避免这个情况，我们意识到：当 \(r=1\) 时，复数 \(z\) 不论怎么取乘幂都不会增加自己的模长，这具有一定的意义。
单位根
我们令 \(\omega_n(n\in N_+)\) 表示复平面上的 \(n\) 次单位根，即 \(x^n=1\) 的解。因此，\(\omega_n^n=1\) 。
然而，很显然，这样的解一共有 \(n\) 个，分别为 \(\exp(i\cdot {2\pi\over n}\cdot 0), \exp(i\cdot {2\pi\over n}\cdot 1), \exp(i\cdot {2\pi\over n}\cdot 2),\cdots ,\exp(i\cdot {2\pi\over n}\cdot (n-1))\) 。为了方便，我们取非零最小幅角的 \(\exp(i\cdot {2\pi\over n})\) 来当作 \(\omega_n\) 的值。
于是，有 \(\omega_n=e^{2\pi i\over n}=\cos {2\pi\over n} + i\sin {2\pi\over n}\) 。
如果我们作图就会发现，一个复数 \(z\) 乘上单位根 \(\omega_n\) ，其实就是将它绕原点逆时针旋转 \(1\over n\) 个圆周。
单位根的性质
- \(\omega_n^0, \omega_n^1, \cdots, \omega_n^{n-1}\) 互不相同。
我们不妨假设 \(0\leq p,q<n\) 且 \(\omega_n^p=\omega_n^q\) 。
于是 \(\omega_n^{p-q}=1=(\omega_n^n)^k=\omega_n^{kn}, k\in Z\) 从而得到 \(p-q=kn\) 。
而由于 \(0\leq p=q+kn<n\wedge 0\leq q<n\) ，因此有且仅有解 \(k=0\) ，故 \(p=q\) 。
- \(\omega_n^k=\omega_n^{k\bmod n}\) 。
\(\omega_n^k=(\omega_n^n)^{\lfloor{k\over n}\rfloor}\cdot \omega_n^{k\bmod n}=1^{\lfloor{k\over n}\rfloor}\cdot \omega_n^{k\bmod n}=\omega_n^{k\bmod n}\) 。
一个重要的推论是 \(\omega_n^{-k}=\omega_n^{n-k}(k>0)\) 。
- \(\omega_{tn}^{tk}=\omega_n^k\)
\(\omega_{tn}^{tk}=\exp({2\pi i\over tn}\cdot tk)=\exp({2\pi i\over n}\cdot k)=\omega_n^k\) 。
一个比较重要的推论是 \(\omega_{2n}^{2k}=\omega_n^k\) （有的地方称它为折半引理）。
- \(\omega_n^{n\over 2}=-1\) 。
\(\omega_n^{n\over 2}=\exp({2\pi i\over n}\cdot {n\over 2})=\exp(\pi i)=-1\) 。
或是由 \(3\) 的结论可以得到 \(\omega_n^{n\over 2}=\omega_{2n}^n=\omega_2^1=-1\) 。
- \(\omega_n^{k+n\over 2}=-\omega_n^k\) 。
由 \(4\) 得 \(\omega_n^{n\over 2}=-1\) ，故：
\(\omega_n^{k+n\over 2}=\omega_n^{n\over 2}\cdot \omega_n^k=-\omega_n^k\) 。
- \(\overline {\omega_n^k}=\omega_n^{-k}\) 。
\(\overline {\omega_n^k}=\overline {\cos{{2\pi\over n}\cdot k} + i\sin{{2\pi\over n}\cdot k}} = \cos{{2\pi\over n}\cdot k}-i\sin {{2\pi\over n}\cdot k}=\cos{{2\pi\over n}\cdot (-k)}+i\sin{{2\pi\over n}\cdot (-k)}=\omega_n^{-k}\)
- \(n\) 的本原单位根有 \(\boldsymbol \varphi(n)\) 个。
本原单位根指满足 \(x^n=1\) 且 \(x^t\neq 1(1\leq t<n\wedge t\in Z)\) 的解。
由于满足 \(x^n=1\) ，故解的选项只有 \(\omega_n^0,\omega_n^1,\cdots,\omega_n^{n-1}\) 这 \(n\) 个。
其次，若 \(\omega_n^k\) 满足 \(x^t=1\) 则 \(\omega_n^{kt}=1\) ，于是 \(n\mid kt\) 。
\(n\mid kt\) 等价于 \({n\over \gcd(n, k)}\mid {k\over \gcd(n,k)}\cdot t\) ，而 \(\gcd({n\over \gcd(n,k)}, {k\over \gcd(n,k)})=1\) 故 \({n\over \gcd(n,k)}\mid t\) 。
因此，令 \(\omega_n^k\) 满足 \(x^t=1\) 的最小 \(t\) 为 \(t_{min}={n\over \gcd(n,k)}\) 。
当 \(\gcd(n,k)>1\) 时，\(0\leq t_{min}<n\) 有解，故 \(\omega_n^k\) 不为本原单位根；否则 \(t_{min}=n\geq n\)，故 \(\omega_n^k\) 为本原单位根。
因此，本原单位根的个数为 \(\displaystyle \sum_{k=0}^{n-1}[\gcd(n,k)=1]=\sum_{k=1}^n[\gcd(n,k)=1]=\boldsymbol \varphi(n)\) 。
单位根反演
我们尝试将 \(x^n=1\) 的某个解 \(\omega_n^k(k\in Z)\) 的 \(0\) 到 \(n-1\) 次幂都加起来，就能得到一个式子：
\(\displaystyle \sum_{i=0}^{n-1}(\omega_n^k)^i\)
这个式子是等比数列求和的形式，我们代入公式即可得到 \(\displaystyle \sum_{i=0}^{n-1}(\omega_n^k)^i={1-(\omega_n^k)^n\over 1-\omega_n^k}={1-\omega_n^{nk}\over 1-\omega_n^k}={1-1\over 1-\omega_n^k}=0\)
但是，这个答案并不是正确的，问题在于等比数列求和的公式只适用于 \(q\neq 1\) 的情况，而 \(n\mid k\) 时 \(\omega_n^k=1\) ，此时 \(\displaystyle \sum_{i=0}^{n-1}(\omega_n^k)^i=\sum_{i=0}^{n-1}1^i=n\) 。
因此，我们综合上述两式得到 \(\displaystyle \sum_{i=0}^{n-1}(\omega_n^k)^i=n\cdot [n\mid k]\) ，其中 \([P]\) 为艾弗森括号，当 \(P\) 命题为真时值为 \(1\) ，否则为 \(0\) 。
我们整理式子即可得到单位根反演的式子：\(\displaystyle [n\mid k]={1\over n}\sum_{i=0}^{n-1}\omega_n^{ki}\)
模意义下的单位根
实数的单位根由于涉及到三角函数，具有精度误差。而很多时候，我们的题目需要将答案对质数 \(P\) 取模，而 \(n\mid (P-1)\) 。
此时，数论中有个概念叫做原根：原根 \(g\) 满足 \(g^k\equiv 1\pmod P\) 的最小正整数解为 \(\boldsymbol \varphi(P)\) 。
根据数论知识，质数必然存在原根 \(g\leq O(P^{1\over 4})\) 。
于是我们只需要枚举出 \(P-1\) 的所有质因子 \(p_i\) ，枚举数值 \(x\) 检验是否满足 \(\forall i\to x^{P-1\over p_i}\not\equiv 1\pmod P\) 。复杂度为 \(O(P^{1\over 4}\cdot \omega(P-1))\) 。
其中 \(\omega(P-1)\) 表示 \(P-1\) 中互异的质因子数，是 \(o(\log (P-1))\) 的。
已知原根了怎么求模意义下的单位根呢？由于 \(g^{P-1}\equiv 1\pmod P\) ，和 \(\omega_n^n=1\) 结构一致；而 \(n\mid (P-1)\) ，因此直接取 \(\omega_n\equiv g^{P-1\over n}\pmod P\) 即可得到了模意义下的单位根。
例题
大意
给定一个长为 \(L\) 的字符串，包含前 \(k(k\geq 2)\) 个字母，共有 \(k^L\) 个串。
对于每一对 \(0\leq i,j<n\) 的二元组 \(\left<i,j\right>\) ，找出有多少个字符串包含 \(p\) 个字母 'a' 和 \(q\) 个字母 'b'，使得 \(p\equiv i\pmod n, q\equiv j\pmod n\) 。
输出矩阵 \(A_{i,j}\) 的每个元素分别表示二元组 \(\left<i,j\right>\) 的答案。
\(2\leq k\leq 26, 1\leq L\leq 10^{18}, 1\leq n\leq 500\) ，答案对 \(P=10^9+7\) 取模，保证 \(n\mid (P-1)\) 。
分析
首先分析 \(A_{i,j}\) 的结果：
我们枚举字母 'a' 和 'b' 的出现次数和位置：
\(\begin{aligned} &A_{i,j}\\ =&\sum_{t=0}^L\dbinom L t(k-2)^{L-t}\sum_{p+q=t}\dbinom t {p,q}[p\equiv i\pmod n][q\equiv j\pmod n]&\text{(we define that }0^0=1\text{)}\\ =&\sum_{t=0}^L\dbinom L t(k-2)^{L-t}\sum_{p+q=t}\dbinom t {p,q}[n\mid (p-i)][n\mid (q-j)]\\ =&\sum_{t=0}^L\dbinom L t(k-2)^{L-t}\sum_{p+q=t}\dbinom t {p,q}{1\over n}\sum_{x=0}^{n-1}\omega_n^{(p-i)x}{1\over n}\sum_{y=0}^{n-1}\omega_n^{(q-j)y}\\ =&{1\over n^2}\sum_{x=0}^{n-1}\sum_{y=0}^{n-1}\omega_n^{-ix-jy}\sum_{t=0}^L\dbinom L t(k-2)^{L-t}\sum_{p+q=t}\dbinom t {p,q}(\omega_n^x)^p\cdot (\omega_n^y)^q\\ =&{1\over n^2}\sum_{x=0}^{n-1}\sum_{y=0}^{n-1}\omega_n^{-ix-jy}\sum_{t=0}^L\dbinom L t(k-2)^{L-t}(\omega_n^x+\omega_n^y)^t\\ =&{1\over n^2}\sum_{x=0}^{n-1}\sum_{y=0}^{n-1}\omega_n^{-ix-jy}(k-2+\omega_n^x+\omega_n^y)^L\\ \end{aligned}\)
当然，若我们暴力计算，单项的复杂度为 \(O(n^2\log L)\) ，总复杂度为 \(O(n^4\log L)\) ，是不可接受的。
首先是我们可以预处理 \(\displaystyle W_{x,y}={(k-2+\omega_n^x+\omega_n^y)^L\over n^2}\) 的结果，复杂度为 \(O(n^2\log L)\) ；
后续我们注意到 \(\displaystyle A_{i,j}=\sum_{x=0}^{n-1}\sum_{y=0}^{n-1}\omega_n^{-ix-jy}W_{x,y}=\sum_{x=0}^{n-1}\omega_n^{-ix}(\sum_{y=0}^{n-1}\omega_n^{-jy}W_{x,y})\)
而 \(\displaystyle \sum_{y=0}^{n-1}\omega_n^{-jy}W_{x,y}\) 是仅有参数 \(x\) 的，因此我们记 \(\displaystyle Y_x=\sum_{y=0}^{n-1}\omega_n^{-jy}W_{x,y}\) ；计算出所有 \(Y_x\) 的复杂度为 \(O(n^2)\) 。
最后，我们由 \(\displaystyle A_{i,j}=\sum_{x=0}^{n-1}\omega_n^{-ix}Y_x\) 可以 \(O(n)\) 时间内求出一个解，总复杂度降低为 \(O(n^3)\) 。
因此该算法的总复杂度为 \(O(n^2\log L)+O(n^2)+O(n^3)=O(n^3)\) 。
模意义下的单位根可以打表获得具体数值，也可以放在程序开头进行计算。
代码
代码是 以前 写的，可能和这一份分析的命名有些差异，但思路是一致的
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef pair<ll, ll> pii;
typedef double db;
#define fi first
#define se second
const int MOD=1e9+9, g=13;
ll k, l, n, omega[512], powit[512][512], t[512][512];
inline ll fpow(ll a,ll x) { ll ans=1; for(;x;x>>=1,a=a*a%MOD) if(x&1) ans=ans*a%MOD; return ans; }
inline int w(int val) { return (val%=n)>=0?omega[val]:omega[n+val]; }
inline int add(int a, int b) { return (a+=b)>=MOD?a-MOD:a; }
inline void pre(){
cin>>k>>l>>n;
omega[0]=1;
omega[1]=fpow(g, (MOD-1)/n);
for(int i=2;i<n;++i) omega[i]=omega[i-1]*omega[1]%MOD;
for(int p=0;p<n;++p)
for(int q=0;q<n;++q)
powit[p][q]=fpow(k-2+w(-p)+w(-q), l);
for(int a=0;a<n;++a)
for(int j=0;j<n;++j){
int res=0;
for(int b=0;b<n;++b)
res=add(res, w(b*j)*powit[a][b]%MOD);
t[a][j]=res;
}
}
int main(){
ios::sync_with_stdio(0);
cin.tie(0); cout.tie(0);
pre();
ll C=fpow(n, MOD+MOD-4);
for(int i=0;i<n;++i){
for(int j=0;j<n;++j){
int res=0;
for(int a=0;a<n;++a)
res=add(res, w(a*i)*t[a][j]%MOD);
cout<<res*C%MOD<<" ";
}
cout<<"\n";
}
cout.flush();
return 0;
}
