Link
首先有一个结论是\(f^p=\epsilon\)，因此\(f=g^{k^{-1}\bmod p}\)。
我们现在将问题转化为给定\(f\)，求\(g=f^k\)。
设\(F(z)=\sum\limits_{n\ge1}\frac{f(n)}{n^z},G(z)=\sum\limits_{n\ge1}\frac{g(n)}{n^z}\)。
定义\(F'(z)=\sum\limits_{n\ge1}\frac{f(n)\ln n}{n^z}\)，即\(f'(n)=f(n)\ln n\)。
我们知道\(kF'G=FG'\)，对比左右两边\([n^z]\)可以得到

\[k\sum\limits_{d|n}f'(d)g(\frac nd)=\sum\limits_{d|n}f(d)g'(\frac nd) \]

我们知道\(g(1)=1,g'(1)=0\)，因此我们可以在\(O(n\log n)\)的时间内递推出\(g\)。
有一个问题是\(\ln\)在\(\mathbb N_+\mapsto\mathbb{F_p}\)意义下没有定义，那么我们换用同样具有完全加性的\(\Omega\)即可。

我们是在\(\mathbb N_+\mapsto\mathbb{F_p}\)下进行操作，因此我们不得不舍弃导数原本的定义。
而事实上我们定义了一个线性变换\(T\)，它与导数一样满足线性性，以及\(T(fg)=T(f)g+T(g)f\)。
我们只需要用到这几条性质即可推出结论\(kF'G=FG'\)。
同样的，对于\(\ln\)，我们在推导的过程中实际上只运用到了它的完全加性，因此用\(\Omega\)替代是没有问题的。

h#include<cctype>
#include<cstdio>
using i64=long long;
const int N=1000007,P=998244353;
char ibuf[1<<24|1],*iS=ibuf;
int Omega[N],f[N],g[N],df[N];
int pow(int a,int b){int r=1;for(;b;b>>=1,a=1ll*a*a%P)if(b&1)r=1ll*a*r%P;return r;}
int read(){int x=0;while(isspace(*iS))++iS;while(isdigit(*iS))(x*=10)+=*iS++&15;return x;}
int main()
{
    fread(ibuf,1,1<<24,stdin);
    int n=read(),k=pow(read(),P-2);
    for(int i=1;i<=n;++i) f[i]=read();
    for(int i=2;i<=n;++i) if(!Omega[i]) for(i64 j=i;j<=n;j*=i) for(int k=1;j*k<=n;++k) ++Omega[j*k];
    for(int i=1;i<=n;++i) df[i]=1ll*f[i]*Omega[i]%P*k%P;
    for(int i=1;i<=n;++i)
    {
	g[i]=i==1? 1:1ll*g[i]*pow(Omega[i],P-2)%P;
	for(int j=2;j<=n/i;++j) g[i*j]=(g[i*j]+1ll*df[j]*g[i])%P;
	for(int j=2;j<=n/i;++j) g[i*j]=(g[i*j]+P-1ll*f[j]*g[i]%P*Omega[i]%P)%P;
    }
    for(int i=1;i<=n;++i) printf("%d ",g[i]);
}