设：A^N=A^(k*x+y)。即：N=k*x+y，x=N/k，y=N%k。由于N为10^9，所以，k取33333左右就可以了，这样x和y的取值都不超过33333了。
则快速幂变成了：A^N mod P=A^(k*x+y) mod P=(A^(k*x) * A^y )mod P=(A^(k*x) mod P * (A^y) mod P) mod P。由于A，k，P都是定值，则A^N mod P 的值只取决于x和y。分别用dpx和dpy来记录A^(k*x) mod P和A^y mod P，则又可得到：
A^N=(dpx *dpy)mod P。
下面是该如何求dpx和dpy呢？
显然：A^y mod P=A*A^(y-1) mod P=(A mod P *(A^(y-1)) mod P) mod P。即：dpy=(A mod P * dp(y-1)) mod P。
同样的：A^(k*x) mod P=(A^k mod P * A^(k*(x-1)) mod P) mod P。即：dpx=(A^k mod P *dp(x-1)) mod P。
#include <stdio.h> #include <string.h> #include <algorithm> using namespace std; long long f[1000009]; long long dpx[33335],dpy[33335]; long long t, A, K, a, b, m, P,n,cnt=1; /* 万万没想到log(n)都能优化，我服！！ */ void init() { dpx[0]=dpy[0]=1; for(int i=1;i<=33333;i++) dpy[i]=((A%P)*dpy[i-1])%P; dpx[1]=dpy[33333]; for(int j=2;j<=33333;j++) dpx[j]=(dpx[1]*dpx[j-1])%P; } long long solve() { long long ans=0; for(int i=1;i<=n;i++) { ans+=(dpx[f[i]/33333]*dpy[f[i]3333])%P; ans%=P; } return ans; } int main() { scanf("%lld",&t); while(t--) { scanf("%lld%lld%lld%lld%lld%lld%lld", &n, &A, &K, &a, &b, &m, &P); f[1]=K; for(int i=2;i<=n;i++) { f[i]=(a*f[i-1]+b)%m; } init(); printf("Case #%lld: %lld\n",cnt++,solve()); } return 0; }
