[EC-Final2019]狄利克雷 k 次根

OneInDark

已于 2022-01-24 15:57:22 修改

阅读量275

点赞数

分类专栏： C++ 数学 # 多项式/生成函数文章标签：我不会做数学

于 2022-01-24 08:44:37 首次发布

本文链接：https://blog.csdn.net/qq_42101694/article/details/122660689

版权

C++ 同时被 3 个专栏收录

547 篇文章 3 订阅

订阅专栏

数学

280 篇文章 1 订阅

订阅专栏

多项式/生成函数

51 篇文章 1 订阅

订阅专栏

题目

传送门 to LOJ

思路

最强的一点在于：指数也可以对 $p = 998244353$ 取模。这真的非常神奇。

利用讲过的狄利克雷生成函数，我们似乎只需要求出 $\exp\frac{\ln F(s)}{k}$ 。但是 $k^{-1}$ 是什么呢？

神奇的事情发生了：对于形式幂级数 $F(x)^p\equiv F(x^p)\pmod{p}\;(p\in\text{Prime})$ ，而对于狄利克雷生成函数有 $F(s)^p\equiv F(sp)\pmod{p}$ 。从积性函数的角度，就是 $f (n)$ 移动到了 $f(n^p)$ 去。

证明方法就有很多了，~~包括暴力打表验证~~。比如模仿形式幂级数的推导，有
$[n^{-s}]F(s)=\sum_{\prod q^{t_q}=n}\frac{p!\prod a_i^{t_i}}{\prod(t_i!)}$

只有 $t_v=p$ 时，右侧的阶乘之比可能模 $p$ 不为零。此时 $n=v^p$ ，系数是 $f(v)^p\equiv f(v)\pmod{p}$ ，也就是说 $f(v^p)$ 的系数是原来的 $f (v)$ ，证毕。

从贝尔级数的角度，就更简单了，由 $F_p(x)^q\equiv F_p(x^q)\pmod{q}\;(q\in\text{Prime})$ 一步到位。但是这就只适用于积性函数了（“贝尔生成函数” 难以推导）。

而一般情况下 $2^p\gg n$ ，所以除 $1$ 以外的最低项 $2^p$ 也已经超出有用的范围了。于是我们可以粗略地理解为，一个函数的狄利克雷卷积 $p$ 次幂只有 $1$ 处有值 $f (1)$ 。而积性函数为底数时， $f (1) = 1$ ，所以指数可以对 $p$ 取模。

如果只需要证上面这玩意儿， $\sf EI$ 给了一个更高级的证明：由于 $\ln$ 和 $\exp$ 算子都是 线性变换，故所有运算应当在 $\Bbb F_p$ 下进行，证毕。

于是我们只需要计算 $F(s)^{k^{-1}}$ 。显然可以用 $\exp\ln$ 组合技。但是还有一种更简单的推导方法（并且同样适用于形式幂级数）：由于导数是良定义的， $F(s)^k]'=kF(s)^{k-1}F'(s)$ 。设 $G(s)=F(s)^k$ ，两边同时求导后乘 $F (s)$ 得
$G^{'} (s) F (s) = k G (s) F^{'} (s)$

显然它也可以做到 $\mathcal O(n\ln n)$ ，并且看上去更好写。不过，作为狄利克雷生成函数的第一道例题，我还是老老实实使用了 $\exp$ 和 $\ln$ 以测试其正确性。

代码

#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
#include <cctype>
using namespace std;
# define rep(i,a,b) for(int i=(a); i<=(b); ++i)
# define drep(i,a,b) for(int i=(a); i>=(b); --i)
typedef long long llong;
inline int readint(){
	int a = 0, c = getchar(), f = 1;
	for(; !isdigit(c); c=getchar())
		if(c == '-') f = -f;
	for(; isdigit(c); c=getchar())
		a = (a<<3)+(a<<1)+(c^48);
	return a*f;
}

const int MOD = 998244353;
inline llong qkpow(llong b,int q){
	llong a = 1;
	for(; q; q>>=1,b=b*b%MOD)
		if(q&1) a = a*b%MOD;
	return a;
}

const int MAXN = 1000001;
bool isPrime[MAXN];
int primes[MAXN], primes_size, coe[MAXN];
void sieve(int n = MAXN-1){
	memset(isPrime+2,true,n-1);
	rep(i,2,n){
		if(isPrime[i]) primes[primes_size ++] = i, coe[i] = 1;
		for(int j=0; j!=primes_size&&primes[j]<=n/i; ++j){
			isPrime[i*primes[j]] = false;
			coe[i*primes[j]] = coe[i]+1;
			if(i%primes[j] == 0) break;
		}
	}
}

void divide(int f[],const int g[],int n){
	const llong inv = qkpow(g[1],MOD-2);
	for(int i=1; i<=n; ++i){
		f[i] = int(inv*f[i]%MOD);
		rep(j,2,n/i) f[i*j] = int((f[i*j]-
			llong(f[i])*g[j]%MOD+MOD)%MOD);
	}
}
inline void jifen(int f[],int n){
	rep(i,1,n) f[i] = int(qkpow(coe[i],MOD-2)*f[i]%MOD);
}
void getLn(const int f[],int g[],int n){
	rep(i,1,n) g[i] = int(llong(f[i])*coe[i]%MOD);
	divide(g,f,n); jifen(g,n);
}

void exponent(const int f[],int g[],int n){
	memset(g+1,0,n<<2); rep(i,g[1]=1,n){
		if(coe[i]) g[i] = int(qkpow(coe[i],MOD-2)*g[i]%MOD);
		rep(j,2,n/i) g[i*j] = int((g[i*j]+
			llong(f[j])*coe[j]%MOD*g[i])%MOD);
	}
}

int f[MAXN], g[MAXN];
int main(){
	sieve(); // actually just need n
	int n = readint(), k = readint();
	rep(i,1,n) f[i] = readint();
	getLn(f,g,n); const llong invk = qkpow(k,MOD-2);
	rep(i,1,n) g[i] = int(invk*g[i]%MOD);
	exponent(g,f,n); rep(i,1,n) printf("%d ",f[i]);
	return 0;
}