[FOI2020WC模拟]看上去很简单

题目

题目描述
道生一，一生二，二生三，三生万物。由此可知，$g_i=3g_{i-1}$ 。

天之道，损有余而补不足。故而真正的 $g_i=3g_{i-1}-g_{i-2}$ 。

世事如棋局局新。$g_0=\alpha ,g_1=\beta$，此皆变者也。

子子孙孙无穷尽也，而值加增，何苦不知 $f_{n,k}=f_{g_n,k-1}$ 之值也！

天下有道，丘不与易也。$f_{n,0}=n$，此不变者也。

数据范围与提示
数据组数 $T\le 10^3$，而 $n\le 10^9$，好在 $k\le 10^2$ 。

为了方便输出，只需要求出模 $p$ 意义下的值。$p\le 10^9$，每组数据都不一样。

思路

好像 $g_n$ 不能直接取模 $p$，因为要求 $g_{g_n}$ 的值……
$$\begin{gathered} G(x)-a-bx=(3x-x^2)\cdot G(x)-3ax \\ (x^2-3x+1)\cdot G(x)=(b-3a)x+a \end{gathered}$$

我们有理由相信 $g_n$ 会循环，因为只有 $p^2$ 种相邻元素的取值。而一旦这个重复，后面就都一样了。可是怎么找啊？

考虑 $a,b$ 的系数分别计算。对于 $b$ 的系数，相当于初值 $a=0,b=1$，即
$$\begin{gathered} (x^2-3x+1)\cdot G_b(x)=x \\ b{g}_n=\frac{1}{\sqrt{5}}\left[{\left(\frac{3+\sqrt{5}}{2}\right)}^n-{\left(\frac{3-\sqrt{5}}{2}\right)}^n\right] \end{gathered}$$

当 $n=p$ 时，$\left(\genfrac{}{}{0ex}{}{n}{i}\right)$ 只有 i ∈ { 0 , p } i\in\{0,p\} i∈{0,p} 可以模 $p$ 不为零，于是只剩 $\frac{(\sqrt{5}{)}^p}{2^{p-1}}$ 。然后外面有个 $\frac{1}{\sqrt{5}}$，总体就是 5 p − 1 2 2 p − 1 ≡ 5 p − 1 2 ∈ { − 1 , 1 } {5^{p-1\over 2}\over 2^{p-1}}\equiv 5^{p-1\over 2}\in\{-1,1\} 2p−152p−1​​≡52p−1​∈{−1,1} 。（特殊情况也找到了，就是 $p=5$ 嘛。）

若 $n=p+1$ 呢？也不难计算出，其值为 $\frac{3}{2}(1+5^{\frac{p-1}{2}})$ 。那么我们可以分类讨论一下：

• 若 $5^{\frac{p-1}{2}}\equiv -1$，那么 $b{g}_{n+1}=0,b{g}_{n+2}=1$，于是它的周期是 $n+1$ 。
• 若 $5^{\frac{p-1}{2}}\equiv 1$，那么 $b{g}_{n+1}=3$，解方程 $3b{g}_n-b{g}_{n-1}=b{g}_{n+1}$ 可知 $b{g}_{n-1}=0$，于是存在一个 $n-1$ 的周期。

现在 $b$ 的周期得证了。我们再来看 $a$ 。它相当于初值 $a=1,b=0$，即
$$\begin{gathered} (x^2-3x+1)\cdot G_a(x)=-3x+1 \\ a{g}_n=\frac{3\sqrt{5}-5}{10}{\left(\frac{3+\sqrt{5}}{2}\right)}^n+\frac{15-3\sqrt{5}}{10}{\left(\frac{3-\sqrt{5}}{2}\right)}^n \end{gathered}$$

类似的，最后都归结到 $5^{\frac{p-1}{2}}$ 上面去。总之一句话：最小循环节要么是 $p-1$ 的因数，要么是 $p+1$ 的因数。虽然考试的时候我是直接打表发现的规律。

当然，如果你更细心一点，你会发现，若 p   m o d   5 ∈ { 1 , 4 } p\bmod 5\in\{1,4\} pmod5∈{1,4} 那么 $5^{\frac{p-1}{2}}\equiv 1$，否则 $5^{\frac{p-1}{2}}\equiv -1$ 。因为 $5^{\frac{p-1}{2}}=\left(\frac{5}{p}\right)=\left(\frac{p}{5}\right)$ 嘛（这是勒让德符号，用了二次互反律），而模 $5$ 的二次剩余仅有 $1,4$ 两个。不知道这个也没关系，暴力检查一下$\sout{p+1}$是不是周期也可以判断。

那么 $p^2$ 呢？打表发现，它是 $p\cdot T_p$ 或者 $T_p$ 。具体的证明差不多还是用卢卡斯定理。推广开来，$T_{p^k}\mid p^k\cdot T_p$ 。我们有理由相信，它是读者自证不难的。

现在，我们知道 $g(n)\mathrm{m}\mathrm{o}\mathrm{d}p$ 循环节是 $T_p$，那么我们求 $g[g(n)]\mathrm{m}\mathrm{o}\mathrm{d}p$ 只需要求 $g(n)\mathrm{m}\mathrm{o}\mathrm{d}{T}_p$ 即可，这是一个递归的过程。

每次递归都要 $\mathcal{O}(\sqrt{p})$ 做质因数分解吗？其实也可以存下来。但是并没有快多少。

而且，根据上面的过程可知，最大的质因数如果次数为 $1$，那么每次都至少变小一半。很快 $p$ 就会塌缩成只有 $2,3,5$ 构成。最后会塌缩到 $2^2\times 3\times 5^t(t\in {\mathbb{N}}^{+})$ 或 $2^2\times 3$ 。

哪怕是 $2\times 3\times 5\times \cdots \times 41=304,250,263,527,210$ 这样的庞然大物，它也只能变成 $10,835,137,422,950,400$，大概是原来的 $35$ 倍而已……但是 $\mathrm{l}\mathrm{c}\mathrm{m}$ 的过程中，丢失的质因数是非常多的。一方面，最大的质因数只能提供一个 $\frac{p+1}{p}$，然后就会立刻被分解（因为 $2\mid (p+1)$，至少要缩小到原来的一半）。另一方面，$2\mid p$ 时，$\mathrm{l}\mathrm{c}\mathrm{m}$ 过程中百分之百会丢失至少一个 $2$，所以整个数字也在缩小到原来的一半。而 $2\ne p$ 肯定有 $2\mid p^{\prime }$，杠杠的！

总结起来就是：变大是不稳定且困难的，变小是容易且稳定的。我们有理由相信，$p$ 不会大的离谱。

接下来，这道题就变成了玄学优化题。请洒潘江，各倾陆海云尔！

代码

最有用的代码是这份打表代码。它帮助我找到了结论。

	int p = 25; // whatever you want
	for(int i=2,g0=0,g1=1; true; ++i){
		g0 = (3ll*g1-g0+p)%p; swap(g0,g1);
		if(g0 == 0 && g1 == 1){
			printf("period = %d",i-1);
			break;
		}
	}

$AC$ 代码反而挺无趣的……

#include <cstdio>
#include <iostream>
#include <cstring>
#include <vector>
#include <algorithm>
#include <map>
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)
# define __int128 long long
inline int readint(){
	int a = 0; char c = getchar(), f = 1;
	for(; c<'0'||c>'9'; c=getchar())
		if(c == '-') f = -f;
	for(; '0'<=c&&c<='9'; c=getchar())
		a = (a<<3)+(a<<1)+(c^48);
	return a*f;
}
inline void writeint(__int128 x){
	if(x > 9) writeint(x/10);
	putchar((x-x/10*10)^48);
}

const int MaxN = 1000002;
vector< int > primes;
int least[MaxN], next_one[MaxN], num[MaxN];
bool isPrime[MaxN]; int period[MaxN];
void sieve(int n = MaxN-1){
	memset(isPrime+2,1,n-1);
	for(int i=2,len=0; i<=n; ++i){
		if(isPrime[i]){
			++ len, primes.push_back(i);
			least[i] = i, num[i] = 1;
			next_one[i] = 1; // end
		}
		for(int j=0; j<len&&primes[j]<=n/i; ++j){
			isPrime[i*primes[j]] = false;
			least[i*primes[j]] = primes[j];
			if(i%primes[j] == 0){
				num[i*primes[j]] = num[i]+1;
				next_one[i*primes[j]] = next_one[i];
				break; // fenjie can be unordered
			}
			num[i*primes[j]] = 1;
			next_one[i*primes[j]] = i;
		}
	}
}

__int128 Mod; // vary
inline __int128 mul(__int128 a,__int128 b){
	a = (a*b-Mod*(__int128)((long double)a/Mod*b+1e-10))%Mod;
	return (a < 0) ? a+Mod : a; // avoid negative
}
struct Matrix{
	__int128 a[2][2];
	Matrix operator * (const Matrix &b) const {
		Matrix c; // no need to clear
		rep(i,0,1) rep(k,0,1) // manual
			c.a[i][k] = (mul(a[i][0],b.a[0][k])
				+mul(a[i][1],b.a[1][k]))%Mod;
		return c;
	}
};
Matrix qkpow(Matrix b,__int128 q){
	Matrix a = b; -- q;
	for(; q; q>>=1,b=b*b)
		if(q&1) a = a*b;
	return a;
}
__int128 qkpow(__int128 b,int q){
	__int128 a = 1;
	for(; q; q>>=1,b=mul(b,b))
		if(q&1) a = mul(a,b);
	return a;
}

Matrix Q, S; int A, B;
__int128 calc(__int128 n,const __int128 &p){
	if(n <= 1) return n ? B : A;
	else Mod = p; // setup
	/* [ g_{i-1},g_i ] * S = [ g_i,g_{i+1} ] */ ;
	S.a[0][0] = 0, S.a[1][0] = 1;
	S.a[0][1] = p-1, S.a[1][1] = 3;
	S = qkpow(S,n-1); // calculate it
	Q.a[0][0] = A%p, Q.a[0][1] = B%p;
	Q.a[1][0] = Q.a[1][1] = 0; Q = Q*S;
	return Q.a[0][1];
}

bool __switch_of_period;
__int128 getPeriod(const __int128 &p){
	if(p == 5) return 10; // the only exception
	if(p%5 == 1 || p%5 == 4) return p-1;
	else return p+1; // only two cases
	if(__switch_of_period){ // enabled
		int rnk = lower_bound(primes.begin(),
			primes.end(),p)-primes.begin();
		if((rnk^primes.size()) && primes[rnk] == p)
			return period[rnk]; // solved
	}
	__int128 gp_1 = calc(p-1,p), gp = calc(p,p);
	if(gp_1 == A%p && gp == B%p) return p-1;
	else return p+1; // believe in conclusion
}

bool miller_rabbin(const int &p){
	if(p <= 7) return p == 2 or
		p == 3 or p == 5 or p == 7;
	return (Mod = p) != 1
		and qkpow(2,p-1) == 1
		and qkpow(3,p-1) == 1
		and qkpow(5,p-1) == 1
		and qkpow(7,p-1) == 1;
}

map< int,int > cur; vector<int> sy;
typedef map< int,int >::iterator ddq;
void getFenJie(int v){
	if(miller_rabbin(v)) return void(++ cur[v]);
	for(int i=0,t; primes[i]<=v/primes[i]&&v>=MaxN; ++i)
		if(v%primes[i] == 0){
			for(t=0; !(v%primes[i]); ++t)
				v /= primes[i];
			cur[primes[i]] += t;
		}
	if(v >= MaxN) return void(++ cur[v]);
	for(; v!=1; v=next_one[v])
		cur[least[v]] += num[v];
}
__int128 hardCalc(__int128 p){
	if(miller_rabbin(p)) return getPeriod(p);
	__int128 res = 1, v; // make lcm right
	for(int i=0; primes[i]<=p/primes[i]&&p>=MaxN; ++i)
		if(p%primes[i] == 0){
			for(v=1; !(p%primes[i]); v*=primes[i])
				p /= primes[i];
			v = v/primes[i]*getPeriod(primes[i]);
			res = res/__gcd(res,v)*v;
		}
	if(p >= MaxN) return res*getPeriod(p);
	Mod = p; // whatever... big is enough
	for(; p!=1; p=next_one[p]){
		v = qkpow(least[p],num[p]-1);
		v *= getPeriod(least[p]);
		res = res/__gcd(res,v)*v;
	}
	return res;
}

const int Threshold = MaxN;
unordered_map< __int128,__int128 > path;
__int128 getG(int n,int k,__int128 p){
	sy.clear(); if(!k) return n%p; // exit
	if(path.count(p)) // recorded
		return calc(getG(n,k-1,path[p]),p);
	__int128 phi = p; ddq it;
	if(!cur.empty() && cur.rbegin()->first < Threshold){
		phi = hardCalc(p); goto FUCK_IT_ALL;
	}
	for(it=cur.begin(); it!=cur.end(); ){
		if(it->first != 5) // 5 won't fade
			phi /= it->first, -- it->second;
		sy.push_back(it->first); // record
		if(!it->second) cur.erase(it ++);
		else ++ it; // move ahead
	}
	for(int i=0,zxy=sy.size(); i<zxy; ++i){
		__int128 v = getPeriod(sy[i]);
		v /= __gcd(phi,v); phi *= v; // lcm
		getFenJie(v); // add it
	}
FUCK_IT_ALL:
	path[p] = phi; // memorize
	return calc(getG(n,k-1,phi),p);
}

int main(){
	freopen("hakugai.in","r",stdin);
	freopen("hakugai.out","w",stdout);
	sieve(); int len = primes.size();
	A = 0, B = 1; // whatever value
	rep(i,0,len-1) // prepare period
		period[i] = getPeriod(primes[i]);
	__switch_of_period = true;
	for(int T=readint(); T; --T){
		A = readint(), B = readint();
		int n = readint(), k = readint();
		int p = readint(); // stupid!!!
		cur.clear(); getFenJie(p);
		writeint(getG(n,k,p)), putchar('\n');
	}
	return 0;
}
/*
2
161412041 518937363 163320992 2 657844984
33527462 107256621 442389778 1 311584514

Ans: 78179392 217296843
*/