Brute-force Algorithm_矩阵快速幂&&欧拉公式*****

Problem Description

Professor Brute is not good at algorithm design. Once he was asked to solve a path finding problem. He worked on it for several days and finally came up with the following algorithm:

Any fool but Brute knows that the function “funny” will be called too many times. Brute wants to investigate the number of times the function will be called, but he is too lazy to do it.

Now your task is to calculate how many times the function “funny” will be called, for the given a, b and n. Because the answer may be too large, you should output the answer module by P.

 

Input

There are multiple test cases. The first line of the input contains an integer T, meaning the number of the test cases.

For each test cases, there are four integers a, b, P and n in a single line.
You can assume that 1≤n≤1000000000, 1≤P≤1000000, 0≤a, b<1000000.

 

Output

For each test case, output the answer with case number in a single line.

 

Sample Input

3 3 4 10 3 4 5 13 5 3 2 19 100

 

Sample Output

Case #1: 2 Case #2: 11 Case #3: 12

【参考】http://www.cnblogs.com/183zyz/archive/2012/05/11/2495401.html

【题意、思路】f[n]=f[n-1]*f[n-2];

f的前面几项可以罗列出来：

a^1*b^0,a^0*b^1,a^1*b^1,a^1*b^2,a^2*b^3....

可以发现a的指数和b的指数均类似于斐波那契数列。

用矩阵的快速幂可以很快的求出第n项a和b的指数分别是多少。

但是这个指数会非常大，存不下来，需要对一个数去模。

这里需要用到一个公式:

A^B%C=A^( B%Phi[C] + Phi[C] )%C   (B>=Phi[C])

Phi[C]表示不大于C的数中与C互质的数的个数，可以用欧拉函数来求：

找到C的所有素因子。

Phi[C]=C*(1-1/q1)*(1-1/q2)*(1-1/q3)*....*(1-1-qk);

q1,q2,q3...qk表示C的素因子。

到这里基本上就能解决了。

#include<iostream>
#include<stdio.h>
#include<string.h>
using namespace std;
const int N=1000005;
long long int a,b,p,n;
int vis[N],prime[N],K;
long long int phi;
struct mat
{
    long long int a[2][2];
};
void init()
{
    memset(vis,0,sizeof(vis));
    for(int i=2;i<=1000;i++)
    {
        if(vis[i]==0)
        {
            for(int j=i+i;j<=1000000;j+=i)
            {
                vis[j]=1;
            }
        }
    }
    K=0;
    for(int j=2;j<=1000000;j++)
        if(vis[j]==0) prime[++K]=j;
}
mat power(mat a1,mat b1)//矩阵乘法
{
    mat c;
    for(int i=0;i<=1;i++)
    {
        for(int j=0;j<=1;j++)
        {
            c.a[i][j]=0;
            for(int k=0;k<=1;k++)
            {
                c.a[i][j]+=a1.a[i][k]*b1.a[k][j];
                if(c.a[i][j]>phi)
                    c.a[i][j]=c.a[i][j]%phi+phi;
            }
        }
    }
    return c;
}
mat mod(mat x,long long int k)
{
    mat ans;
    ans.a[0][0]=1;
    ans.a[0][1]=0;
    ans.a[1][0]=0;
    ans.a[1][1]=1;
    while(k)
    {
        if(k%2) ans=power(ans,x);
        x=power(x,x);
        k/=2;
    }
    return ans;
}
long long int mod1(long long int a,long long int k)
{
    long long int ans;
    ans=1;
    while(k)
    {
        if(k%2)
        {
            ans=ans*a;
            ans%=p;
        }
        a=a*a;
        a%=p;
        k/=2;
    }
    return ans%p;
}
int main()
{
    int t,cas=1;
    long long int aa,bb,num,num2,num1;
    mat now,ans;
    init();
    scanf("%d",&t);
    while(t--)
    {
        scanf("%lld%lld%lld%lld",&a,&b,&p,&n);
        printf("Case #%d: ",cas++);

        if(n==1) {printf("%lld
",a%p);continue;}
        else if(n==2) {printf("%lld
",b%p);continue;}
        else if(n==3) {printf("%lld
",a*b%p);continue;}
        if(p==1) {printf("0
");continue;}
        else
        {
            now.a[0][0]=0;
            now.a[0][1]=1;
            now.a[1][0]=1;
            now.a[1][1]=1;
//  A^B % C = A ^ ( B % phi[C] + phi[C] ) %C  ( B >= phi[C] ) ,phi[C]表示与C互质的数的个数
            phi=1;
            num=p;
            for(int i=1; i<=K; i++)
            {
                if(prime[i]>p) break;
                if(p%prime[i]==0)
                {
                    phi*=(prime[i]-1);
                    num/=prime[i];
                }

            }
//phi[C]=C*(1-1/p1)*(1-1/p2)*...*(1-1/pt);p1,p2,...pt表示C的素因子
            phi*=num;//phi表示phi[c]
            ans=mod(now,n-3);//求指数

            num1=ans.a[1][1];//a的指数
            num2=ans.a[0][1]+ans.a[1][1];//b的指数
            if(num2>phi) num2=num2%phi+phi;

            aa=mod1(a,num1);//a^num1%p
            bb=mod1(b,num2);//b^num2%p
            printf("%lld
",aa*bb%p);
        }

    }
    return 0;
}