题目
定义数列 ({g_n}) ：
对于 (kin mathbb N,nin mathbb Z) ，定义 (f_{n,k}) 为：
现在给定 (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_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 的分析过程：
-
小质数：直接打表可得 (C(2)=3,C(3)=8,C(5)=10) ；
-
大质数：
-
如果 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) 。
-
如果 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) 。
-
-
合数：
这里不展开讲，只记一下结论吧：
设 (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) 层，求出每一层的模数，再倒过来用矩阵加速求每一层的值。
小结：
- 注意区分数列循环节和指数循环节！
- 这类分析循环节的方法，通过此题已经证明是比较通用的了。
代码
#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;
}