[多校赛20210406]迫害 DJ

题目

定义数列 ({g_n}) ：

[g_n= egin{cases} a&n=0\ b&n=1\ 3g_{n-1}-g_{n-2}&n>1 end{cases} ]

对于 (kin mathbb N,nin mathbb Z) ，定义 (f_{n,k}) 为：

[f_{n,k}= egin{cases} n&k=0\ f_{g_n,k}&k>0 end{cases} ]

现在给定 (a,b,n,k,p) ，请求出 (f_{n,k}mod p) 的值。

数据范围：对于 (100\%) 的数据，满足 (1le n,kle 10^9, 0le a,b<p,1le kle 100) ，且每个数据点中有最多 1000 组询问。

分析

显然可以发现题目要求的就是：

[g_{g_{dots_{g_n}}} ]

一样的套娃。但是如果直接递归显然是没法做的，我们必须缩减下标：

这里有一个误区，需要注意：

此处的 ({g_n}) 显然有通项 (g_n=c_0phi^n+c_1hatphi^n) ，其中 (phi) 和 (hatphi) 为特征方程 (x^2-3x+1=0) 的两个根，不难解出 (phi=frac{3+sqrt 5}{2},hatphi=frac{3-sqrt5}{2}) 。

那么可能会想到扩展欧拉定理，但是很不幸不可行。这是因为该定理要求底数在模 (p) 的整数域中。如果 (5) 不是模 (p) 的二次剩余就是错误的！

事实上，基于 (g) 类似于 Fibonacci 数列的结构，我们可以猜想它在 (mod P(P>1)) 是存在循环节的！ 设 (C(P)) 为模 (P) 意义下的循环节，那么 (g_{g_n}mod P=g_{g_nmod {C(P)}}mod P) 。

将 ((g_n,g_{n+1})) 看作一个有序数对，合法的数对只有 (P^2) 个。

并且可以说明 ({g_n}) 是纯循环的，因为在循环节中 (g_n) 的前一个元素唯一确定。

类似于 Fibonacci 的分析过程：

1. 小质数：直接打表可得 (C(2)=3,C(3)=8,C(5)=10) ；

2. 大质数：

   1. 如果 5 是 (P) 的二次剩余，根据二次互反律 (left(frac 5 P ight)left(frac P 5 ight)= -1^{frac{p-1}{2}cdotfrac{5-1}{2}}) ，可知此时 (left(frac{P}{5} ight)equiv 1pmod P) ，那么可以得出 (Pequiv pm 1pmod 5) 。

      此时 (phi,hatphi) 都是有定义的，因此在模 (P) 意义下推导：

      [egin{aligned} phi^P&=phi\ hatphi^P&=hatphi end{aligned} ]

      因而可以得出 (g_{P-1}=g_0Rightarrow C(P)=P-1) 。

   2. 如果 5 不是 (P) 的二次剩余，可以此时 (left(frac P 5 ight)equiv -1pmod P) ，那么可以得出 (Pequiv 2pmod 5) 或 (Pequiv 3pmod 5)。

      此时 (phi) 是没有定义的，但是我们可以将 (sqrt 5) 扩入模 (P) 整数域中：

      [egin{aligned} phi^P &=frac{(3+sqrt 5)^P}{2^P}\ &=frac{3^P+left(sqrt 5 ight)^P}{2^P}\ &=frac{3+5^{frac{P-1}{2}} imes sqrt 5}{2}\ &=frac{3-sqrt 5}2=hatphi end{aligned} ]

      从而有：

      [egin{aligned} phi^{2P+2} &=(phi^P)^2phi^2\ &=hatphi^2phi^2\ &=(phi imes hatphi)^2=1 end{aligned} ]

      同理也可以得出 (hatphi^{2P+2}=1) 。因此得出 (g_{2P+2}=g_0Rightarrow C(P)=2P+2) 。

3. 合数：

   这里不展开讲，只记一下结论吧：

   设 (n=p_1^{e_1} imes p_2^{e_2} imes dots imes p_k^{e_k}) ，那么 (C(n)=operatorname{lcm}_{1le jle k}{p_j^{e_j-1}C(p_j)}) ，其中 (p_1,p_2,dots,p_k) 均为质数。

每次对 (P) 进行分解低于 (O(sqrt P)) ，而 (P) 不会超过 long long ，因此可以递归 (k) 层，求出每一层的模数，再倒过来用矩阵加速求每一层的值。

小结：

1. 注意区分数列循环节和指数循环节！
2. 这类分析循环节的方法，通过此题已经证明是比较通用的了。

代码

#include <cstdio>
#include <cstring>
#include <iostream>

#define rep( i, a, b ) for( int i = (a) ; i <= (b) ; i ++ )
#define per( i, a, b ) for( int i = (a) ; i >= (b) ; i -- )

typedef long long LL;

const int MAXN = 1e6 + 5, MAXLOG = 105;

template<typename _T>
void read( _T &x )
{
	x = 0; char s = getchar(); int f = 1;
	while( s < '0' || '9' < s ) { f = 1; if( s == '-' ) f = -1; s = getchar(); }
	while( '0' <= s && s <= '9' ) { x = ( x << 3 ) + ( x << 1 ) + ( s - '0' ), s = getchar(); }
	x *= f;
}

template<typename _T>
void write( _T x )
{
	if( x < 0 ) putchar( '-' ), x = - x;
	if( 9 < x ) write( x / 10 );
	putchar( x % 10 + '0' );
}

template<typename _T>
_T MIN( const _T a, const _T b )
{
	return a < b ? a : b;
}

inline LL Mul( LL x, LL v ); 
inline LL Sub( LL x, LL v );
inline LL Add( LL x, LL v );

struct Matrix
{
	LL mat[2][2] = {};
	
	Matrix& operator *= ( const Matrix &b ) { return *this = *this * b; }
	
	Matrix operator * ( const Matrix &b ) const
	{
		Matrix ret;
		for( int i = 0 ; i < 2 ; i ++ )
			for( int k = 0 ; k < 2 ; k ++ )
				if( mat[i][k] ) for( int j = 0 ; j < 2 ; j ++ )
					ret.mat[i][j] = Add( ret.mat[i][j], Mul( mat[i][k], b.mat[k][j] ) );
		return ret;
	}
};

Matrix T;

LL MOD[105];

int prime[MAXN], pn;
bool isPrime[MAXN];

LL mod;
LL A, B, N, K, P;

inline LL Sub( LL x, LL v ) { return ( x -= v ) < 0 ? x + mod : x; }
inline LL Add( LL x, LL v ) { return ( x += v ) >= mod ? x - mod : x; }
inline LL Gcd( const LL a, const LL b ) { return b ? Gcd( b, a % b ) : a; }
inline LL LCM( const LL a, const LL b ) { return a / Gcd( a, b ) * b; }
inline LL Mul( LL a, LL b ) { return ( a * b - ( LL )( ( long double ) a / mod * b ) * mod + mod ) % mod; }

void Init()
{
	T.mat[0][0] = 3, T.mat[0][1] = mod - 1;
	T.mat[1][0] = 1, T.mat[1][1] = 0;
}

void EulerSieve( const int n = MAXN - 5 )
{
	for( int i = 2 ; i <= n ; i ++ )
	{
		if( ! isPrime[i] ) prime[++ pn] = i;
		for( int j = 1 ; j <= pn && 1ll * i * prime[j] <= n ; j ++ )
		{
			isPrime[i * prime[j]] = true;
			if( ! ( i % prime[j] ) ) break;
		}
	}
}

inline Matrix Qkpow( Matrix b, LL idx )
{
	Matrix ret;
	ret.mat[0][0] = ret.mat[1][1] = 1;
	while( idx )
	{
		if( idx & 1 ) ret *= b;
		b *= b, idx >>= 1;
	}
	return ret;
}

inline LL G( const LL n, const LL m ) 
{ 
	if( n == 0 ) return A;
	if( n == 1 ) return B;
	mod = m, Init(); 
	Matrix tmp = Qkpow( T, n );
	return Add( Mul( tmp.mat[1][0], B ), Mul( tmp.mat[1][1], A ) );
}

LL GetCir( LL p )
{
	LL ret = 1, tmp;
	for( int i = 1 ; i <= pn && 1ll * prime[i] * prime[i] <= p ; i ++ )
		if( ! ( p % prime[i] ) )
		{
			p /= prime[i];
			if( prime[i] == 2 ) tmp = 3;
			if( prime[i] == 3 ) tmp = 8;
			if( prime[i] == 5 ) tmp = 10;
			if( prime[i] % 5 == 1 || prime[i] % 5 == 4 ) tmp = prime[i] - 1;
			if( prime[i] % 5 == 2 || prime[i] % 5 == 3 ) tmp = 2 * prime[i] + 2;
			while( ! ( p % prime[i] ) ) tmp *= prime[i], p /= prime[i];
			ret = LCM ( ret, tmp );
		}
	if( p > 1 ) 
	{
		if( p == 2 ) tmp = 3;
		if( p == 3 ) tmp = 8;
		if( p == 5 ) tmp = 10;
		if( p % 5 == 1 || p % 5 == 4 ) tmp = p - 1;
		if( p % 5 == 2 || p % 5 == 3 ) tmp = 2 * p + 2;
		ret = LCM( ret, tmp );
	}
	return ret;
}

signed main()
{
	freopen( "hakugai.in", "r", stdin );
	freopen( "hakugai.out", "w", stdout );
	int T;
	read( T );
	EulerSieve();
	while( T -- )
	{
		read( A ), read( B ), read( N ), read( K ), read( P );
		MOD[K] = P; LL res = N;
		per( i, K - 1, 0 ) MOD[i] = GetCir( MOD[i + 1] );
		rep( i, 1, K ) 
			res = G( res, MOD[i] );
		write( res ), putchar( '
' );
	}
	return 0;
}