CF1626F

,

Solution

观察到 $k$ 很小，我们可以对值域按照模 $lcm(1,…,k)$ 进行分类。同一类的数进行相同操作，减掉的数的大小都是相同的。

令 $M=lcm(1,…,k)$

根据期望的线性性，我们可以对每类数分别计算。

令 $f[i][j]$ 为进行了 $i$ 轮，选到 $j$ 这个数的概率。
$$
\begin{cases}f[i+1][j]+=\dfrac{n-1}{n}\times f[i][j]\f[i+1][j-j\bmod i]+=\dfrac{1}{n}\times f[i][j]\end{cases}
$$
累加答案，$\sum\limits_{1\le i\le k}\sum\limits_{0\le j< M} f[i][j]\times j$

注意到最后一轮，模 $k$ 之后不再产生贡献。为了减少常数，$M=lca(1,…,k-1)$ 即可。

code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
#include<bits/stdc++.h>
#define int long long
#define ull unsigned long long
#define maxn 13000005
#define put() putchar('\n')
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
using namespace std;
inline void read(int &x){
    int f=1;x=0;char c=getchar();
    while (c<'0'||c>'9') {if (c=='-') f=-1;c=getchar();}
    while (c>='0'&&c<='9') {x=x*10+c-'0';c=getchar();}
    x*=f;
}
namespace Debug{
	Tp void _debug(char* f,Ty t){cerr<<f<<'='<<t<<endl;}
	Ts void _debug(char* f,Ty x,Ar... y){while(*f!=',') cerr<<*f++;cerr<<'='<<x<<",";_debug(f+1,y...);}
	Tp ostream& operator<<(ostream& os,vector<Ty>& V){os<<"[";for(auto& vv:V) os<<vv<<",";os<<"]";return os;}
	#define gdb(...) _debug((char*)#__VA_ARGS__,__VA_ARGS__)
}using namespace Debug;
const int mod=998244353;
int base=720720*17;
int n,a[maxn],A,B,P,k,Invn;
int nums[maxn],f[2][1000005];
inline int power(int x,int y) {
	int ans=1;
	while (y) {
		if (y&1) ans=ans*x%mod;
		x=x*x%mod;y>>=1;
	}
	return ans;
}
int sum,ans,fn[20];
signed main(void){
	int i,j;
	read(n);read(a[1]);read(A);read(B);read(k);read(P);
	for (base=1,fn[0]=1,i=1;i<k;i++) fn[i]=fn[i-1]*n%mod,base=base/__gcd(i,base)*i;
	Invn=power(n,mod-2);
	for (i=2;i<=n;i++) {
		a[i]=(a[i-1]*A+B)%P;
		sum=(sum+a[i]-a[i]%base)%mod;
		f[1][a[i]%base]=(f[1][a[i]%base]+Invn)%mod;
	}
	sum=(sum+a[1]-a[1]%base)*Invn%mod*k%mod;
	ans=sum;
	f[1][a[1]%base]=(f[1][a[1]%base]+Invn)%mod;
	for (i=1;i<=k;i++) {
		int now=i&1,pre=now^1;
		memset(f[pre],0,sizeof(f[pre]));
		for (j=0;j<=base;j++) {
			f[pre][j]=(f[pre][j]+f[now][j]*(1-Invn+mod))%mod;
			f[pre][j-j%i]=(f[pre][j-j%i]+f[now][j]*Invn)%mod;
			ans=(ans+f[now][j]*j%mod)%mod;
		}
	}
	printf("%lld\n",(ans*power(n,k))%mod);
	return 0;
}