HDU5667 - Sequence 矩阵快速幂 + 费马小定理

f[1] = 1

f[2] = a ^ b

f[n]\quad =\quad { a }^{ b }\quad *\quad { f[n-1] }^{ c }\quad *\quad f[n-2]

其实不是很好的去想到取log的 

两边同时取log

\log _{ a }{ f[n]\quad =\quad \log _{ a }{ ({ a }^{ b }\quad *\quad { f[n-1] }^{ c }\quad *\quad f[n-2]) } } \\ \qquad \qquad \quad \quad =\quad b\quad +\quad c\log _{ a }{ f[n-1]\quad + } \log _{ a }{ f[n-1] } \quad \qquad \\ \\ F[n]\quad =\quad \log _{ a }{ f[n] } \qquad \\ F[n]\quad =\quad b\quad +\quad cF[n-1]\quad +\quad F[n-2]\\ F[3]\quad =\quad b\quad +\quad cF[2]\quad +\quad F[1]\\ \begin{pmatrix} F[2] & F[1] & b \end{pmatrix}\begin{pmatrix} c & 1 & 0 \\ 1 & 0 & 0 \\ 1 & 0 & 1 \end{pmatrix}\quad =\quad \begin{pmatrix} F[3] & F[2] & b \end{pmatrix}\\

然后 F[2] = b  F[1] = 0

则 f[n] = a^(F[n]) % p  

费马小定理 : 

① 判断素数,对于大素数的判定,Miller-Rabin 素数判定
②求解逆元 ,设a模p的逆元为x,则a*x≡1(mod p) ,(a,p)=1;由费马小定理可以知道x=a^(p-2)
③对于计算ab(modp)ab(modp) 可简化
            对于素数p,任取跟他互素的数a,有a^(p-1)(mod p)=1
            所以任取b,有a^b%p=a^(b%(p-1))(%p)从而简化运算。

本文来自 杨美人 的博客 ,全文地址请点击:https://blog.csdn.net/qq_40679299/article/details/80596406?utm_source=copy

 

然后 f[n] = a^(F[n] mod (p - 1) ) mod p

代码:

//
//  HDU5667 - Sequence.cpp
//  数论
//
//  Created by Terry on 2018/9/29.
//  Copyright © 2018年 Terry. All rights reserved.
//

#include <stdio.h>
typedef long long LL;
long long read(){
    long long x = 0, f = 1;
    char ch=getchar();
    while(ch < '0' || ch > '9'){if(ch=='-') f = -1; ch = getchar();}
    while(ch >= '0' && ch <= '9'){x = x * 10 + ch - '0'; ch = getchar();}
    return x * f;
}
const int maxn = 3;
LL mod;
struct Matrix{
    long long int m[maxn][maxn];
}unit;
Matrix operator * (Matrix a, Matrix b){
    Matrix ret;
    LL x;
    int n = 3;
    for(int i = 0; i < n; ++i){
        for(int j = 0; j < n; ++j){
            x = 0;
            for(int k = 0; k < n; ++k){
                x += ((LL)a.m[i][k] * b.m[k][j]) % mod;
            }
            ret.m[i][j] = x % mod;
        }
    }
    return ret;
}
void init_unit(){
    // 单位矩阵
    for(int i = 0; i < maxn; ++i){
        unit.m[i][i] = 1;
    }
}
Matrix pow_mat(Matrix a, LL n){
    Matrix ans = unit;
    while (n) {
        if(n & 1){
            ans = ans * a;
        }
        a = a * a;
        n >>= 1;
    }
    return ans;
}
LL multi(LL a, LL b, LL Mod){
    // a * b % Mod
    LL ret = 0;
    while(b > 0){
        if(b & 1){
            ret = (ret + a) % Mod;
        }
        b >>= 1;
        a = (a << 1) % Mod;
    }
    return ret;
}

LL quick_mod(LL a, LL b, LL Mod){
    LL ans = 1;
    while(b){
        if(b & 1){
            ans = multi(ans, a, Mod);
            b--;
        }
        b /= 2;
        a = multi(a, a, Mod);
    }
    return ans;
}
int main(){
    init_unit();
    LL T = read();
    while (T--) {
        LL n = read();
        LL a = read();
        LL b = read();
        LL c = read();
        mod = read();
        if(n == 1){
            printf("1\n");
        }
        else if(n == 2){
            printf("%lld\n", quick_mod(a, b, mod));
        }
        else{
            if(a % mod == 0){
                printf("0\n");
                continue;
            }
            Matrix ans, A;
            A.m[0][0] = c; A.m[0][1] = 1; A.m[0][2] = 0;
            A.m[1][0] = 1; A.m[1][1] = 0; A.m[1][2] = 0;
            A.m[2][0] = 1; A.m[2][1] = 0; A.m[2][2] = 1;

            ans.m[0][0] = b; ans.m[0][1] = 0; ans.m[0][2] = b;

            mod--; // 费马小定理
            A = pow_mat(A, n - 2);

            ans = ans * A;

            mod++;

            printf("%lld\n", quick_mod(a, ans.m[0][0]%(mod - 1), mod));
        }
    }
    return 0;
}

 

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值