程式師世界 >> 編程語言 >> C語言 >> C++ >> C++入門知識 >> hdu 3221 Brute-force Algorithm(快速冪取模，矩陣快速冪求fib)

hdu 3221 Brute-force Algorithm(快速冪取模，矩陣快速冪求fib)

編輯：C++入門知識

hdu 3221 Brute-force Algorithm(快速冪取模，矩陣快速冪求fib)

一晚上搞出來這麼一道題。。Mark。

給出這麼一個程序，問funny函數調用了多少次。

我們定義數組為所求：f[1] = a，f[2] = b, f[3] = f[2]*f[3]......f[n] = f[n-1]*f[n-2]。對應的值表示也可為a^1*b^0%p，a^0*b^1%p，a^1*b^1%p,.....a^fib[n-3]*b^fib[n-2]%p。即a,b的指數從n=3以後與fib數列一樣。

因為n很大，fib[n]也想當大。a^fib[n]%p可以利用a^fib[n]%p = a^(fib[n]%phi[p]+phi[p])%p進行降冪，條件時fib[n]>=phi[p]。求fib[n]%phi[p]可以構造矩陣，利用矩陣快速冪求fib[n]%phi[p]。

#include 
#include 
#include
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#define LL long long
#define _LL __int64
#define eps 1e-12
#define PI acos(-1.0)
#define C 240
#define S 20
using namespace std;
const int maxn = 110;

struct matrix
{
    LL mat[3][3];
    void init()
    {
        memset(mat,0,sizeof(mat));
        for(int i = 0; i < 2; i++)
            mat[i][i] = 1;
    }
} m;

LL a,b,p,n,phi_p;
LL fib[10000000];

//phi[p]
LL Eular(LL num)
{
    LL res = num;
    for(int i = 2; i*i <= num; i++)
    {
        if(num%i == 0)
        {
            res -= res/i;
            while(num%i == 0)
                num /= i;
        }
    }
    if(num > 1)
        res -= res/num;
    return res;
}
//矩陣相乘
matrix mul_matrix(matrix x, matrix y)
{
    matrix ans;
    memset(ans.mat,0,sizeof(ans.mat));
    for(int i = 0; i < 2; i++)
    {
        for(int k = 0; k < 2; k++)
        {
            if(x.mat[i][k] == 0) continue;
            for(int j = 0; j < 2; j++)
            {
                ans.mat[i][j] = (ans.mat[i][j] + x.mat[i][k]*y.mat[k][j])%phi_p;
            }
        }
    }
    return ans;
}
//a^t%phi_p
LL pow_matrix(LL t)
{
    matrix a,b;
    a.mat[0][0] = a.mat[0][1] = a.mat[1][0] = 1;
    a.mat[1][1] = 0;
    b.init();
    while(t)
    {
        if(t&1)
            b = mul_matrix(a,b);
        a = mul_matrix(a,a);
        t >>= 1;
    }
    return b.mat[0][0];
}
//a^t%p
LL pow(LL a, LL t)
{
    LL res = 1;
    a %= p;
    while(t)
    {
        if(t&1)
            res = res*a%p;
        a = a*a%p;
        t >>= 1;
    }
    return res;
}
//a^fib[t]%p轉化為a^(fib[t]%phi[p]+phi[p])%p,fib[t] >= phi[p]。
LL solve(LL a, LL t)
{
    fib[0] = 1;
    fib[1] = 1;
    int i;
    for(i = 2; i <= t; i++)
    {
        fib[i] = fib[i-1] + fib[i-2];
        if(fib[i] >= phi_p)
            break;
    }
    if(i <= t) //當滿足條件fib[t] >= phi[p]時，進行降冪
    {
        LL c = pow_matrix(t) + phi_p;
        return pow(a,c);
    }
    else
        return pow(a,fib[t]);
}

int main()
{
    int test;
    scanf(%d,&test);
    for(int item = 1; item <= test; item++)
    {
        scanf(%lld %lld %lld %lld,&a,&b,&p,&n);
        printf(Case #%d: ,item);
        if(n == 1)
        {
            printf(%lld
,a%p);
            continue;
        }
        if(n == 2)
        {
            printf(%lld
,b%p);
            continue;
        }
        if(n == 3)
        {
            printf(%lld
,a*b%p);
            continue;
        }
        if(p == 1)
        {
            printf(0
);
            continue;
        }
        phi_p = Eular(p);
        LL res = solve(a,n-3)*solve(b,n-2)%p;
        printf(%lld
,res);
    }
    return 0;
}