hdu3802:Ipad,IPhone 之快速矩阵幂+二分快速幂

题目地址:http://acm.hdu.edu.cn/showproblem.php?pid=3802

此题需要用到以下定理:

1)a^n=a^(n % phi(p))(mod p)(p为素数),

2)以及一些如何从特征方程求出递推公式的定理:

An+2=pA(n+1)+qAn,

得起特征方程:x^2=px+q,

若该方程有相异根，则An=c1*s1^n+c2*s2^n,

若该方程有相同根，则An=(c1+nc2)*s1^n。

设zn=x^n+y^n;

x=(sqrt(a)+sqrt(b))^2;y=(sqrt(a)-sqrt(b))^2;

得到p=x+y=2(a+b),q=-x*y=-(a-b)^2;

则

An+2=2(a+b)A(n+1)-(a-b)^2An;

可由矩阵求其解。

当然还需具备一些其他的知识点！

好了，废话不多说:

给出我参考别人的代码再自己写下的代码:

#include<iostream>
#include<algorithm>
#include<cstdio>

using namespace std;
#define bigt __int64
#define SIZE 2

bigt orgf[1][2];
bigt orgz[1][2];
bigt ansf[1][2];
bigt ansz[1][2];
bigt mod;
typedef struct {
    bigt a[SIZE][SIZE];
}mat;
bigt mod2;
mat unit, fibm, zm;
bigt modgod(bigt kk, bigt mod){
    while (kk<0){
        kk += mod * 1000000;
    }
    kk %= mod;
    return kk;
}
bigt FastM(bigt a, bigt b, bigt p){
    bigt d = 1;
    while (b){
        if (b & 1) {
            d = (d*a) % p;
        }
        b >>= 1;
        a = (a*a) % p;
    }
    return d;
}
mat matMul(mat a, mat b, bigt modd){
    mat c;
    int i, j, k;
    for (i = 0; i < SIZE; i++){

        for (j = 0; j < SIZE; j++){
            c.a[i][j] = 0;
            for (k = 0; k < SIZE; k++){
                c.a[i][j] += a.a[i][k] * b.a[k][j];
                c.a[i][j] = modgod(c.a[i][j], modd);
            }

        }
    }
    return c;
}

void init(bigt a, bigt b){
    int i, j;
    for (i = 0; i < SIZE; i++){
        for (j = 0; j < SIZE; j++){
            unit.a[i][j] = (i == j);
            fibm.a[i][j] = 1;
        }
    }
    fibm.a[0][0] = 0;
    zm.a[0][0] = 0;
    zm.a[0][1] = modgod(-(a - b)*(a - b), mod);
    zm.a[1][0] = 1;
    zm.a[1][1] = (2 * (a + b)) % mod;
    orgz[0][0] = 2;
    orgz[0][1] = (2 * (a + b)) % mod;
    orgf[0][0] = 0;
    orgf[0][1] = 1;
}
void computeFib(bigt a[][2], bigt k, mat fibm){
    mat mm = unit;

    while (k){
        if (k & 1){
            mm = matMul(mm, fibm, mod2);
        }
        k >>= 1;
        fibm = matMul(fibm, fibm, mod2);
    }
    int i, j, t;
    for (i = 0; i < 1; i++){
        for (j = 0; j < 2; j++){
            a[i][j] = 0;
            for (t = 0; t < 2; t++){
                a[i][j] += orgf[i][t] * mm.a[t][j];
                a[i][j] = modgod(a[i][j], mod2);
            }
            //a[i][j] %= mod;
            //	printf("%I64d ",a[i][j]);
        }
        //printf("\n");
    }
}
void computeZm(bigt a[][2], bigt k, mat zm){
    mat mm = unit;
    while (k){
        if (k & 1){
            mm = matMul(mm, zm, mod);

        }
        k >>= 1;
        zm = matMul(zm, zm, mod);
    }
    int i, j, t;
    for (i = 0; i < 1; i++){
        for (j = 0; j < 2; j++){
            a[i][j] = 0;
            for (t = 0; t < 2; t++){
                a[i][j] += orgz[i][t] * mm.a[t][j];
                a[i][j] = modgod(a[i][j], mod);
            }
            //a[i][j] %= mod;
        }
    }
}
int main(){
    int T;
    scanf("%d", &T);
    while (T--){
        bigt a, b, n;
        scanf("%I64d%I64d%I64d%I64d", &a, &b, &n, &mod);
        //n--;
        //	printf("%I64d %I64d %I64d %I64d\n",a,b,n,mod);
        mod2 = mod - 1;
        init(a, b);
        computeFib(ansf, n, fibm);
        bigt fib = ansf[0][1];
        //	printf("fib%I64d ",fib);
        //	fib=fib%mod2;
        computeZm(ansz, fib, zm);
        bigt z = ansz[0][0] % mod;
        bigt p1 = FastM(a, mod / 2, mod) % mod;
        bigt p2 = FastM(b, mod / 2, mod) % mod;
        if (p1 == -1 || p2 == -1) { puts("0"); continue; }
        p1++;
        p2++;
        p1 = (p1*p2) % mod;
            bigt ans = (p1*z) % mod;
        printf("%I64d\n", ans);
    }
    return 0;
}