[多校赛20210406]迫害 DJ
题目
定义数列 \(\{g_n\}\) ：
对于 \(k\in \mathbb N,n\in \mathbb Z\) ，定义 \(f_{n,k}\) 为：
现在给定 \(a,b,n,k,p\) ，请求出 \(f_{n,k}\bmod p\) 的值。
数据范围：对于 \(100\%\) 的数据，满足 \(1\le n,k\le 10^9, 0\le a,b<p,1\le k\le 100\) ，且每个数据点中有最多 1000 组询问。
分析
显然可以发现题目要求的就是：
一样的套娃。但是如果直接递归显然是没法做的，我们必须缩减下标：
这里有一个误区，需要注意：
此处的 \(\{g_n\}\) 显然有通项 \(g_n=c_0\phi^n+c_1\hat\phi^n\) ，其中 \(\phi\) 和 \(\hat\phi\) 为特征方程 \(x^2-3x+1=0\) 的两个根，不难解出 \(\phi=\frac{3+\sqrt 5}{2},\hat\phi=\frac{3-\sqrt5}{2}\) 。
那么可能会想到扩展欧拉定理，但是很不幸不可行。这是因为该定理要求底数在模 \(p\) 的整数域中。如果 \(5\) 不是模 \(p\) 的二次剩余就是错误的！
事实上，基于 \(g\) 类似于 Fibonacci 数列的结构，我们可以猜想它在 \(\bmod P(P>1)\) 是存在循环节的！ 设 \(C(P)\) 为模 \(P\) 意义下的循环节，那么 \(g_{g_n}\bmod P=g_{g_n\bmod {C(P)}}\bmod 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\right)\left(\frac P 5\right)= -1^{\frac{p-1}{2}\cdot\frac{5-1}{2}}\) ，可知此时 \(\left(\frac{P}{5}\right)\equiv 1\pmod P\) ，那么可以得出 \(P\equiv \pm 1\pmod 5\) 。
此时 \(\phi,\hat\phi\) 都是有定义的，因此在模 \(P\) 意义下推导：
\[\begin{aligned} \phi^P&=\phi\\ \hat\phi^P&=\hat\phi \end{aligned} \]因而可以得出 \(g_{P-1}=g_0\Rightarrow C(P)=P-1\) 。
-
如果 5 不是 \(P\) 的二次剩余，可以此时 \(\left(\frac P 5\right)\equiv -1\pmod P\) ，那么可以得出 \(P\equiv 2\pmod 5\) 或 \(P\equiv 3\pmod 5\)。
此时 \(\phi\) 是没有定义的，但是我们可以将 \(\sqrt 5\) 扩入模 \(P\) 整数域中：
\[\begin{aligned} \phi^P &=\frac{(3+\sqrt 5)^P}{2^P}\\ &=\frac{3^P+\left(\sqrt 5\right)^P}{2^P}\\ &=\frac{3+5^{\frac{P-1}{2}}\times \sqrt 5}{2}\\ &=\frac{3-\sqrt 5}2=\hat\phi \end{aligned} \]从而有：
\[\begin{aligned} \phi^{2P+2} &=(\phi^P)^2\phi^2\\ &=\hat\phi^2\phi^2\\ &=(\phi\times \hat\phi)^2=1 \end{aligned} \]同理也可以得出 \(\hat\phi^{2P+2}=1\) 。因此得出 \(g_{2P+2}=g_0\Rightarrow C(P)=2P+2\) 。
-
-
合数：
这里不展开讲，只记一下结论吧：
设 \(n=p_1^{e_1}\times p_2^{e_2}\times \dots\times p_k^{e_k}\) ，那么 \(C(n)=\operatorname{lcm}_{1\le j\le 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( '\n' );
}
return 0;
}
