单位根反演 小猪佩奇学数学

题面
一个比较纯粹的推式子题目，第一步先将整除拆开：
${\sum }_{i=0}^n\left(\genfrac{}{}{0ex}{}{n}{i}\right)p^i\frac{1}{k}(i-i\mathrm{m}\mathrm{o}\mathrm{d}k)$
拆开括号并枚举余数：
$\frac{1}{k}({\sum }_{i=0}^n(\genfrac{}{}{0ex}{}{n}{i})p^{i}i-(\genfrac{}{}{0ex}{}{n}{i})p^i{\sum }_{j=0}^{k-1}[(i-j)\mathrm{m}\mathrm{o}\mathrm{d}k=0]j)$
利用组合数公式以及单位根反演得：
$\frac{1}{k}({\sum }_{i=1}^n(\genfrac{}{}{0ex}{}{n-1}{i-1})n{p}^i-(\genfrac{}{}{0ex}{}{n}{i})p^i{\sum }_{j=0}^{k-1}\frac{j}{k}{\sum }_{t=0}^{k-1}{\omega }_k^{t(i-j)})$
根据二项式定理整理：
$\frac{1}{k}(np(p+1{)}^{n-1}-\frac{1}{k}{\sum }_{j=0}^{k-1}j{\sum }_{t=0}^{k-1}{\omega }_k^{-tj}{\sum }_{i=0}^n\left(\genfrac{}{}{0ex}{}{n}{i}\right)p^i{\omega }_k^{ti})$
再用二项式定理：
$\frac{1}{k}(np(p+1{)}^{n-1}-\frac{1}{k}{\sum }_{j=0}^{k-1}{\sum }_{t=0}^{k-1}j{\omega }_k^{-tj}(p{\omega }_k^t+1{)}^n)$
根据这个式子的结构可以用任意模长DFT来求解，但又由于是等差乘等比的结构，可以使用特殊的方法：
$\frac{1}{k}(np(p+1{)}^{n-1}-\frac{1}{k}{\sum }_{j=0}^{k-1}{\sum }_{t=0}^{k-1}(p{\omega }_k^t+1{)}^n{\sum }_{i=0}^{j-1}{\omega }_k^{-tj})$
交换求和顺序：
$\frac{1}{k}(np(p+1{)}^{n-1}-\frac{1}{k}{\sum }_{t=0}^{k-1}(p{\omega }_k^t+1{)}^n{\sum }_{i=0}^{k-2}{\sum }_{j=i+1}^{k-1}{\omega }_k^{-tj})$
利用等比数列求和：
$\frac{1}{k}(np(p+1{)}^{n-1}-\frac{1}{k}{\sum }_{t=0}^{k-1}(p{\omega }_k^t+1{)}^n{\sum }_{i=0}^{k-2}\frac{{\omega }_k^{-tk}-{\omega }_k^{-t(i+1)}}{{\omega }_k^{-t}-1})$
整理式子，得到：
$\frac{1}{k}(np(p+1{)}^{n-1}-\frac{1}{k}{\sum }_{t=0}^{k-1}\frac{(p{\omega }_k^t+1{)}^n}{{\omega }_k^{-t}-1}{\sum }_{i=0}^{k-1}(1-{\omega }_k^{-t(i+1)}))$
根据单位根的性质：
$\frac{1}{k}(np(p+1{)}^{n-1}-\frac{1}{k}{\sum }_{t=0}^{k-1}\frac{(p{\omega }_k^t+1{)}^n}{{\omega }_k^{-t}-1}k)$
最终得到：
$\frac{1}{k}(np(p+1{)}^{n-1}-{\sum }_{i=0}^{k-1}\frac{(p{\omega }_k^i+1{)}^n}{{\omega }_k^{-i}-1})$
这里的式子当 $i$ 取$0$时是错误的，中途不能使用等比数列求和公式。实际答案应该为：
$\frac{1}{k}(np(p+1{)}^{n-1}-\frac{k-1}{2}(p+1{)}^n-{\sum }_{i=1}^{k-1}\frac{(p{\omega }_k^i+1{)}^n}{{\omega }_k^{-i}-1})$
时间复杂度 $O(k{\log }_{2}n)$，空间复杂度 $O(1)$

#include<iostream>
using namespace std;
#define R register int
#define L long long
#define I inline
#define P 998244353
I int PowMod(int x,int y){
	int res=1;
	while(y!=0){
		if((y&1)==1){
			res=(L)res*x%P;
		}
		y>>=1;
		x=(L)x*x%P;
	}
	return res;
}
int main(){
	int n,p,k,w,pw=1,ans=0;
	cin>>n>>p>>k;
	w=PowMod(3,(P-1)/k);
	for(R i=1;i!=k;i++){
		pw=(L)pw*w%P;
		ans=((L)PowMod(((L)p*pw+1)%P,n)*PowMod(PowMod(pw,P-2)+P-1,P-2)+ans)%P;
	}
	ans=(499122177ll*(k-1)%P*PowMod(p+1,n)+ans)%P;
	ans=((L)n*p%P*PowMod(p+1,n-1)-ans+P)%P;
	cout<<(L)ans*PowMod(k,P-2)%P;
	return 0;
}