Ipad,IPhone（矩陣求遞推項+歐拉定理）

Ipad,IPhone

Time Limit: 2000/1000 MS (Java/Others) Memory Limit: 32768/32768 K (Java/Others)

Total Submission(s): 100 Accepted Submission(s): 46

Problem Description In ACM_DIY, there is one master called 「Lost」. As we know he is a 「-2Dai」, which means he has a lot of money.    Well, Lost use Ipad and IPhone to reward the ones who solve the following problem.    In this problem, we define F( n ) as :    Then Lost denote a function G(a,b,n,p) as Here a, b, n, p are all positive integer! If you could tell Lost the value of G(a,b,n,p) , then you will get one Ipad and one IPhone!

Input The first line is one integer T indicates the number of the test cases. (T <= 100) Then for every case, only one line containing 4 positive integers a, b, n and p. (1 ≤a, b, n, p≤2*10 9 , p is an odd prime number and a,b < p.)

Output Output one line,the value of the G(a,b,n,p) .

Sample Input 4 2 3 1 10007 2 3 2 10007 2 3 3 10007 2 3 4 10007

Sample Output 40 392 3880 9941

Author AekdyCoin

Recommend notonlysuccess

/*
題意：給出如圖所示的式子，讓你求出值

初步思路：斐波那契數列後推一位，直接求是不能夠的，由於幾十項就會爆，因此用到矩陣求斐波那契數列，式子的前半部分能夠直接用快速冪求得
    如今有一個問題，就是X^t mod q是否是等於 X^(t mod q) 先爆一發試試

#錯誤：上面的想法顯然是錯誤的！

#改進：後邊的比較麻煩,後半部分用矩陣求遞推項，這幾天看這個好麻煩啊，轉專業沒學線性代數，構造矩陣就麻煩死了。
    首先令    Xn=(sqrt(a)+sqrt(b))^(2*fn)
            Yn=(sqrt(a)-sqrt(b))^(2*fn)
    而後得    X1=a+b+2*sqrt(ab)
            Y1=a+b-2*sqrt(ab)
        得    X1+Y1=2*(a+b)
            X1*Y1=(a-b)^2
由韋達定理
求特徵方程：u^2-2*(a+b)*u+(a-b)^2=0;

得特徵方程    Zn^2=2*(a+b)*Zn-1 - (a-b)^2*Zn-2

由此構造矩陣 Zn=(Z0,Z1)*{ 0 , -(a-b)^2 }^n
                        { 1 , 2*(a+b)  }    
                        
*/
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
ll mod;
/******************快速冪*************************/
ll my_pow(ll a,ll b){
    ll ans = 1;
    while(b){
        if (b%2)
            ans = (ans*a)%mod;
        b/=2;
        a=(a*a)%mod;
    }
    return ans;
}
/******************快速冪*************************/

/*******************矩陣快速冪**********************/
struct Mar{
    ll mat[2][2];
};
Mar E={1,0,0,1},unit2={0,1,1,1};
Mar mul(Mar a,Mar b){
    ll fumod(ll);
    Mar ans={0,0,0,0};
    for (int i = 0; i < 2; ++i)
        for (int j = 0; j < 2; ++j)
            for (int k = 0; k <2; ++k){
                ans.mat[i][j]+=(a.mat[i][k]*b.mat[k][j])%mod;
                ans.mat[i][j]%=mod;
        }
    return ans;
}
Mar pow2(Mar a,ll n){
    Mar ans = E;
    while(n){
        if (n%2)
            ans = mul(ans,a);
        n/=2;
        a=mul(a,a);
    }
    return ans;
}
/*******************矩陣快速冪**********************/

ll fib(ll n){//求斐波那契數列
    Mar uu={1,1,1,1};
    Mar p=pow2(unit2,n);
    p=mul(uu,p);
    return p.mat[0][0];
}
int main(){
    //freopen("in.txt","r",stdin);
    ll T,a,b,n,p;
    scanf("%lld",&T);
    while(T--){
        scanf("%lld%lld%lld%lld",&a,&b,&n,&p);
        mod = p;
        //前邊的兩項
        ll temp1=(my_pow(a,(p-1)/2)+1)%mod;
        ll temp2=(my_pow(b,(p-1)/2)+1)%mod;
        
        if (!temp1 || !temp2){cout << "0" << endl;continue;}
        --mod;
        
        ll fn=fib(n);
        ++mod;
        
        //矩陣求後邊的兩項
        Mar unit={2,2*(a+b)%mod,1,1};
        Mar tmp1={0,-(a-b)*(a-b)%mod,1,2*(a+b)%mod};
        Mar p=pow2(tmp1,fn);
        p=mul(unit,p);
        ll pn=p.mat[0][0];
        ll ans = pn%mod*temp1%mod*temp2%mod;
        while(ans < 0)ans+=mod;
        printf("%lld\n",ans);
    }
    return 0;
}