這幾天做了一些Fibnacci的題目，基本思想是構造矩陣，對于求余的時候，模上的數比較小的時候可以用循環節來做。
1、http://acm.hdu.edu.cn/showproblem.php?pid=1021
求fib(n) mod 3是否為0，構造矩陣，代碼略。
2、http://acm.hdu.edu.cn/showproblem.php?pid=1568
求fib的前4位數，當n比較小的時候，直接求；當n比較大的時候，用Fibnacci數的通項公式，((1-sqrt(5))/2)^n -> 0可以直接略去，然后使用一般的求一個數的前幾位數的方法求出即可，代碼略。
3、http://acm.hdu.edu.cn/showproblem.php?pid=3306
定義s(n) = fib(0)^2+fib(1)^2+......+fib(n)^2，求s(n)%10007 ，構造一個4階的矩陣，代碼略。
4、http://acm.hdu.edu.cn/showproblem.php?pid=2814
用循環節來做，要想清楚，我做這題的時候糾結了一兩天。。AC大牛的題，厲害！
5、http://acm.tzc.edu.cn/acmhome/problemdetail.do?&method=showdetail&id=3012
這道題是第三題的進化版，求 a^s(x) mod n，當時馬上想到了歐拉定理，但是發現 a,n不一定互素，決定使用離散對數，感覺還是有問題，topsky是用離散對數過的。我糾結了一個禮拜，想了很多種方法，感覺都有問題，最后問了zsut 的yygy，他告訴我一個重要結論，對于一般的情況下（a，n不互素） a^phi( n ) % n = a^( k*phi( n ) ) % n 。有了這個結論，這道題就簡單了:(1)當s( x ) >= n 時，進行化簡 s( x ) = s( x ) % phi( n ) + phi( n )，然后用快速冪。(2)當s( x ) < n 時，直接計算。 代碼如下：

#include <stdio.h>
#include <iostream>
#include <algorithm>
#include <queue>
#include <map>
#include <vector>
#include <cmath>
#define MAXN 100005
using namespace std;
typedef __int64 LL;
LL p[MAXN];
bool vis[MAXN];
int len;
LL s[100],fib[100];
void prime( )
{
    LL i,j,k;
    len=0;
    memset( vis, false, sizeof( vis ) );
    for( i=2, k=4; i < MAXN; ++i, k+=i+i-1 )
    {
        if( !vis[i] )
        {
            p[len++]=i;
            if( k < MAXN )for( j=k; j < MAXN; j+=i )vis[j]=true;
        }
    }
}

void init( )
{
    prime( );
    fib[1] = 1, fib[2] = 2;
    s[1] = 1, s[2] = 5;
    for( int i = 3; i <= 20; i++ )
    {
        fib[i] = fib[i-1] + fib[i-2];
        s[i] = s[i-1] + fib[i] * fib[i];
    }
}

LL euler( LL n )     //在刷出素數表的前提下的歐拉函數
{
    LL phi=n;
    for( LL i = 0; (LL)p[i] * p[i] <= n; i++ )
    {
        if( n%p[i]==0 )
        {
            while( !( n % p[i] ) ) n /= p[i];
            phi = phi / p[i] * ( p[i] - 1 );
        }
    }
    if( n > 1 ) phi = phi / n * ( n - 1 );
    return phi;
}

LL pow_mod( LL a, LL x, LL n )
{
    LL ans = 1, d = a % n;
    while( x )
    {
        if( x & 1 ) ans = ( ans * d ) % n;
        d = ( d * d ) % n;
        x >>= 1;
    }
    return ans;
}

void matrix_mul( LL a[][4], LL b[][4], LL p )
{
    int i,j,k;
    LL c[4][4];
    memset( c, 0, sizeof( c ) );
    for( i=0; i<4; i++ )
        for( j=0; j<4; j++ )
            for( k=0; k<4; k++ )
                c[i][j]+=a[i][k]*b[k][j];
    for( i=0; i<4; i++ )
        for( j=0; j<4; j++ )
            a[i][j]=c[i][j]%p;
}

LL Sx( LL n, LL p )
{
    int i;
    LL I[4][4],a[4][4]; 
    memset( I, 0, sizeof( I ) );
    for( i=0; i<4; i++ )
        I[i][i]=1;
    a[0][0]=1,a[0][1]=1,a[0][2]=0,a[0][3]=0;
    a[1][0]=0,a[1][1]=1,a[1][2]=1,a[1][3]=2;
    a[2][0]=0,a[2][1]=1,a[2][2]=0,a[2][3]=0;
    a[3][0]=0,a[3][1]=1,a[3][2]=0,a[3][3]=1;
    while( n )
    {
        if( n & 1 ) matrix_mul( I, a, p );
        matrix_mul( a, a, p );
        n>>=1;
    }
    return (I[0][0]+I[0][1]*4+I[0][2]+I[0][3]*2)%p;
}


int main( int argc, char *argv[] )
{
    LL a, x, n, PHI, ans, N;
    init( );
    while( 1 )
    {
        scanf("%I64d%I64d%I64d",&a,&x,&n);
        if( a == 0 && x == 0 && n == 0 )break;
        PHI = euler( n );
        if( x>=1 && x<=20 ) ans = pow_mod( a, s[x], n );
        else 
        {
            N = Sx( x - 1, PHI ) + PHI;
            ans = pow_mod( a, N, n );
        }
        printf("%I64d\n",ans);
    }
    return 0;
}