这几天做了一些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;
}