题目：http://www.lydsy.com/JudgeOnline/problem.php?id=2432
我们考虑一下所求数列mod k意义下的情况（样例）：
1 1 2 3 5 0 5 5 3 0 3 3 6 2 …
发现每个-1的地方都是0，每两个这种相邻地方的中间是一个类FIB数列，然后我们可以拆成一堆小数列的组合：
1 1 2 3 5 0
5 5 3 0
3 3 6 2 …我们发现每个小数列只要开头确定了，那么这个数列就确定了，下一个数列的开头也确定了，然后小数列的开头只有k中情况，所以如果小数列数量大于k，一定会出现循环，然后对于一个小数列，就是一个斐波那契数列乘上一个常数a，设Fi a = 1 ( mod k ，那么Fi就是a的乘法逆元，如果Fi不存在，那么就这个小数列会一直持续下去，直接矩阵算掉，否则用Fi来查i，然后我就卡了，最后跑去VFK的BLOG找了个神奇的定理，模G意义下的斐波那契数列循环节不会超过6G,那么就相当好做了先预处理出一定长度的FIB数列，算出Fi最早出现的位置，然后每次查一下就好了，然后一个数列一个数列的跳，如果循环了就把循环用矩阵跳掉，这样似乎就可以O(k log k)了QAQ。
代码（注意处理小数列的循环的时候，我们维护模p意义下的数要用矩阵快速幂去算，因为这个卡了好久额QAQ）：
#include <cstdio>
#include <algorithm>
#include <cstring>
using namespace std ;
#define REP( i , l , r ) for ( int i = l ; i <= r ; ++ i )
#define rep( i , x ) for ( int i = 0 ; i ++ < x ; )
#define DOWN( i , r , l ) for ( int i = r ; i >= l ; -- i )
#define Rep( i , x ) for ( int i = 0 ; i < x ; ++ i )
const int maxk = 1010000 ;
typedef long long ll ;
struct mat {
int n , m , a[ 4 ][ 4 ] ;
mat( ) {
n = m = 0 ;
memset( a , 0 , sizeof( a ) ) ;
}
inline void Init( int _n ) {
n = m = _n ;
memset( a , 0 , sizeof( a ) ) ;
rep( i , n ) a[ i ][ i ] = 1 ;
}
} e ;
inline mat mul( const mat &x , const mat &y , int mod ) {
mat ret ;
ret.n = x.n , ret.m = y.m ;
rep( i , x.n ) rep( j , y.m ) rep( k , x.m ) {
ret.a[ i ][ j ] = ( ret.a[ i ][ j ] + ( ll ) x.a[ i ][ k ] * y.a[ k ][ j ] % mod ) % mod ;
}
return ret ;
}
inline mat power( mat x , ll cnt , int mod ) {
mat ret ; ret.Init( x.n ) ;
for ( ; cnt ; cnt >>= 1 ) {
if ( cnt & 1 ) ret = mul( ret , x , mod ) ;
x = mul( x , x , mod ) ;
}
return ret ;
}
inline void Init_mat( ) {
e.n = e.m = 2 , e.a[ 1 ][ 1 ] = 0 , e.a[ 1 ][ 2 ] = e.a[ 2 ][ 1 ] = e.a[ 2 ][ 2 ] = 1 ;
}
inline void getf( int f1 , int f2 , ll cnt , int mod , int &x , int &y ) {
if ( cnt == 1 ) {
y = f1 ; return ;
}
if ( cnt == 2 ) {
x = f1 , y = f2 ; return ;
}
mat rec = power( e , cnt - ll( 2 ) , mod ) , ret ;
ret.n = 2 , ret.m = 1 , ret.a[ 1 ][ 1 ] = f1 % mod , ret.a[ 2 ][ 1 ] = f2 % mod ;
ret = mul( rec , ret , mod ) ;
x = ret.a[ 1 ][ 1 ] , y = ret.a[ 2 ][ 1 ] ;
}
inline int gcd( int x , int y ) {
if ( x < y ) swap( x , y ) ;
for ( int z ; y ; z = y , y = x % y , x = z ) ;
return x ;
}
void exgcd( int a , int b , int &x , int &y ) {
if ( ! b ) {
x = 1 , y = 0 ; return ;
}
int tx , ty ; exgcd( b , a % b , tx , ty ) ;
x = ty , y = tx - ( a / b ) * ty ;
}
inline int Inv( int a , int mod ) {
if ( gcd( a , mod ) != 1 ) return 0 ;
int x , y ;
if ( a > mod ) exgcd( a , mod , x , y ) ; else exgcd( mod , a , y , x ) ;
if ( x > 0 ) x %= mod ;
if ( x < 0 ) x += mod * ( ( - x ) / mod + ( ( - x ) % mod > 0 ) ) ;
return x ;
}
int pos[ maxk ] , k , p , f1 , f2 ;
ll n ;
int last[ maxk ] , cnt = 0 ;
ll sum[ maxk ] ;
int main( ) {
Init_mat( ) ;
memset( pos , 0 , sizeof( pos ) ) ;
scanf( "%lld%d%d" , &n , &k , &p ) ;
f1 = f2 = 1 ;
for ( int t , i = 3 , j = 6 * k ; i <= j ; ++ i ) {
t = ( f1 + f2 ) % k ;
f1 = f2 ; f2 = t ;
if ( ! pos[ t ] ) pos[ t ] = i ;
}
memset( last , 0 , sizeof( last ) ) ;
sum[ cnt = 0 ] = 0 ;
if ( n <= 2 ) {
printf( "1\n" ) ; return 0 ;
}
for ( int a = 1 , b1 = 1 , b2 = 1 , len , inv , x , y ; ; ) {
inv = Inv( a , k ) ;
if ( ! inv || ! pos[ inv ] || n <= pos[ inv ] ) {
getf( b1 , b2 , n , p , x , y ) ;
if ( y % k == 1 ) y = ( y - 1 + p ) % p ;
printf( "%d\n" , y ) ;
return 0 ;
}
len = pos[ inv ] ;
++ cnt ;
sum[ cnt ] = sum[ cnt - 1 ] + ll( len ) , last[ a ] = cnt ;
getf( a , a , len - 1 , k , x , y ) ; a = y ;
getf( b1 , b2 , len , p , x , y ) ; y = ( y - 1 + p ) % p ;
b1 = ( x + y ) % p ; b2 = ( y + b1 ) % p ;
n -= len ;
if ( last[ a ] && n > sum[ cnt ] - sum[ last[ a ] - 1 ] ) {
mat Ex , Ea , Eb , Et ;
Ea.n = Ea.m = Eb.n = Eb.m = 3 ;
Ea.a[ 1 ][ 1 ] = 0 , Ea.a[ 1 ][ 2 ] = 1 , Ea.a[ 1 ][ 3 ] = 0 ;
Ea.a[ 2 ][ 1 ] = 1 , Ea.a[ 2 ][ 2 ] = 1 , Ea.a[ 2 ][ 3 ] = 0 ;
Ea.a[ 3 ][ 1 ] = 0 , Ea.a[ 3 ][ 2 ] = 0 , Ea.a[ 3 ][ 3 ] = 1 ;
Eb.a[ 1 ][ 1 ] = 1 , Eb.a[ 1 ][ 2 ] = 1 , Eb.a[ 1 ][ 3 ] = p - 1 ;
Eb.a[ 2 ][ 1 ] = 1 , Eb.a[ 2 ][ 2 ] = 2 , Eb.a[ 2 ][ 3 ] = p - 2 ;
Eb.a[ 3 ][ 1 ] = 0 , Eb.a[ 3 ][ 2 ] = 0 , Eb.a[ 3 ][ 3 ] = 1 ;
Ex.Init( 3 ) ;
REP( i , last[ a ] , cnt ) {
Ex = mul( power( Ea , sum[ i ] - sum[ i - 1 ] - 2 , p ) , Ex , p ) ;
Ex = mul( Eb , Ex , p ) ;
}
ll sz = sum[ cnt ] - sum[ last[ a ] - 1 ] , t ;
if ( n % sz ) {
t = n / sz ;
n %= sz ;
} else {
t = n / sz - 1 ;
n = sz ;
}
Et.n = 3 , Et.m = 1 , Et.a[ 1 ][ 1 ] = b1 , Et.a[ 2 ][ 1 ] = b2 , Et.a[ 3 ][ 1 ] = 1 ;
Et = mul( power( Ex , t , p ) , Et , p ) ;
b1 = Et.a[ 1 ][ 1 ] , b2 = Et.a[ 2 ][ 1 ] ;
}
}
return 0 ;
}
