矩阵快速幂的应用 hdu5451

发表于 |

上次已经写过了矩阵快速幂，这次再写的原因是因为此题hdu5451用了很多的黑科技，因此还是值得记录的。

问题简述如下:
$$ y = (5+ 2\sqrt{6})^{1+2^x}$$
where $0 \leq x < 2^{32} $ and a prime number $p(p<46337)$,calculate $r = \floor{y} \;mod\; p$.

实际上我们需要计算的是
$$ r = (5+ 2\sqrt{6})^{1+2^x} + (5 - 2\sqrt{6})^{1+2^x} -1 $$
设 $a_n = (5+ 2\sqrt{6})^{1+2^x} + (5 - 2\sqrt{6})^{1+2^x} $ 则我们有
$$ a_0=2,a_1=10,a_n = a_{n-1} +a_{n-2}$$
最终答案即为
$$ (a_{1+2^n} - 1) \;mod\; M $$
另一方面，由递推关系式子，我们知道
$$ \left( \begin{matrix} a_{n} \\ a_{n-1} \end{matrix} \right) = A \left( \begin{matrix} a_{n-1} \\ a_{n-2} \end{matrix} \right) = A^{n-1} \left( \begin{matrix} a_1 \\ a_0 \end{matrix} \right) $$
因此
$$ \left( \begin{matrix} r \\ * \end{matrix} \right) = A^{2^x} \left( \begin{matrix} 10 \\ 2 \end{matrix} \right) $$
由于我们需要的是最终结果模 $p$.
而 $GL(n,p)$可看成 $n$ 阶方针到自身的可逆变换全体。推理易知道
$$ |GL(n,p)| = \prod_{i=0}^{n-1} (p^n-p^i) $$
对于此题中 $ mp = |GL(2,p)| = (p^2-1)(p^2-p) $，由群的Lagrange定理知 $A^{|GL(2,p)|}=I_2$。因此 当矩阵$A$可逆时，我们只需计算$A^{2^x \mod mp}$ 用快速幂解决即可，但是这里有一个问题就是 mp 有可能很大，导致快速幂乘法的时候会溢出，于是乘法用快速加法来实现，具体见代码
其实对于 $n=2$ 的情形，周期为 $p^2-1$。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
//#pragma comment(linker,"/STACK:10240000,10240000")
#include <cstring>
#include <cstdio>
#include <cstdlib>
#include <ctime>
#include <cmath>
#include <iostream>
#include <algorithm>
#include <vector>
#include <map>
#include <set>
#include <queue>
#include <stack>
#include <string>
#include <bitset>
using namespace std;
typedef long long LL;
typedef pair<int,int> PII;
typedef pair<LL,LL> PLL;
#define clr(a,b) memset(a,b,sizeof(a))
#define MP make_pair
#define PB push_back
#define lrt rt<<1
#define rrt rt<<1|1
#define lson l,m,lrt
#define rson m+1,r,rrt
/*------ Welcome to visit blog of dna049: http://dna049.com ------*/
LL mod;
class Matrix{
public:
    const static int N = 2; //col
    LL a[N][N];
    Matrix(LL a00=0,LL a01=0,LL a10=0,LL a11=0){
        a[0][0]=a00;a[0][1]=a01;
        a[1][0]=a10;a[1][1]=a11;
    }
    Matrix operator+(const Matrix& A)const{
        Matrix R;
        for(int i=0;i<N;++i){
            for(int j=0;j<N;++j){
                R.a[i][j]=a[i][j]+A.a[i][j];
                if(R.a[i][j]>=mod)   R.a[i][j]-=mod;
            }
        }
        return R;
    }
    Matrix operator*(const Matrix& A)const{
        Matrix R;
        for(int i=0;i<N;++i){
            for(int k=0;k<N;++k){
                for(int j=0;j<N;++j){
                    R.a[i][j] = (R.a[i][j]+a[i][k]*A.a[k][j])%mod;
                }
            }
        }
        return R;
    }
};
Matrix pow(Matrix A,LL n){
    Matrix R(1,0,0,1);
    while(n){
        if(n&1) R=R*A;
        n>>=1;  A=A*A;
    }
    return R;
}
LL mul(LL x,LL n,LL p){
    LL r = 0;
    while(n){
        if(n&1) r+=x;
        n>>=1;  x<<=1;
        if(r>=p)    r-=p;
        if(x>=p)    x-=p;
    }
    return r;
}
LL pow_mod(LL x,LL n,LL p){
    LL r = 1;
    while(n){
        if(n&1) r=mul(r,x,p);
        n>>=1;  x=mul(x,x,p);
    }
    return r;
}
int main(){
//    freopen("/Users/dna049/Desktop/AC/in","r",stdin);
//    freopen("/Users/dna049/Desktop/AC/out","w",stdout);
    int T;
    scanf("%d",&T);
    for(int t=1;t<=T;++t){
        unsigned int x;
        cin>>x>>mod;
        LL p = mod;
        p = (p*p-1)*(p*p-p);
        Matrix A(10,-1,1,0);
        Matrix R = pow(A,pow_mod(2,x,p));
        LL r = (R.a[0][0]*10+R.a[0][1]*2-1)%mod;
        if(r<0) r+=mod;
        cout<<"Case #"<<t<<": "<<r<<endl;
    }
    return 0;
}