矩阵快速幂
要点
- 两个矩阵做乘法的时候，累积到一个全 0 的矩阵上。
- 在做快速幂的时候需要设为单位矩阵，相当于 1。
- 矩阵中出现负数，记得要及时取正数。
习题
2019牛客第五场 B. generator 1 （矩阵快速幂）
链接：https://ac.nowcoder.com/acm/contest/885/B
思路：矩阵快速幂，先找到转移矩阵
[
a
b
1
0
]
×
[
f
n
−
1
f
n
−
2
]
=
[
f
n
f
n
−
1
]
\begin{bmatrix} a & b \\ 1 & 0 \end{bmatrix} \times \begin{bmatrix} f_{n-1}\\ f_{n-2} \end{bmatrix}= \begin{bmatrix} f_{n}\\ f_{n-1} \end{bmatrix}
[a1b0]×[fn−1fn−2]=[fnfn−1]
变形之后， a [ 2 ] [ 1 ] × f 1 + a [ 2 ] [ 2 ] × f 0 a[2][1]\times f_1+a[2][2]\times f_0 a[2][1]×f1+a[2][2]×f0 就是答案
[ a b 1 0 ] n × [ f 1 f 0 ] = [ f n + 1 f n ] \begin{bmatrix} a&b\\ 1&0 \end{bmatrix}^n \times \begin{bmatrix} f_1\\ f_0 \end{bmatrix}= \begin{bmatrix} f_{n+1}\\ f_n \end{bmatrix} [a1b0]n×[f1f0]=[fn+1fn]
但是n会很大 1 0 ( 1 0 6 ) 10^{(10^6)} 10(106)，可以采用 10 进制倍增
#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int maxn=1e5+5;
int x0,x1,a,b,mod;
string n;
struct Node
{
int a[3][3];
} A[5],res;
Node mul(Node A,Node B)
{
Node res;
memset(res.a,0,sizeof(res.a));
for(int i=1; i<=2; ++i)
for(int j=1; j<=2; ++j)
for(int k=1; k<=2; ++k)
res.a[i][j]=(1ll*res.a[i][j]+1ll*A.a[i][k]*B.a[k][j]%mod)%mod;
return res;
}
int main()
{
cin>>x0>>x1>>a>>b>>n>>mod;
A[0].a[1][1]=a,A[0].a[1][2]=b;
A[0].a[2][1]=1,A[0].a[2][2]=0;
res.a[1][1]=1,res.a[1][2]=0;
res.a[2][1]=0,res.a[2][2]=1;
int len=n.size();
for(int i=len-1; i>=0; --i)
{
A[1]=mul(A[0],A[0]);
A[2]=mul(A[1],A[1]);
A[3]=mul(A[2],A[2]);
int d=n[i]-'0';
for(int j=0; j<4; ++j)
if(d>>j&1) res=mul(res,A[j]);
A[0]=mul(A[1],A[3]);
}
int ans=(1ll*res.a[2][1]*x1%mod+1ll*res.a[2][2]*x0%mod)%mod;
printf("%d\n",ans);
return 0;
}
递推例题
题意：给定前 k k k 项 a 1 , … , a k a_1,\dots,a_k a1,…,ak 以及 x 1 , … , x k ， 其 中 x i = ∑ j = 1 k a j x i − j x_1,\dots,x_k，其中 x_i=\sum_{j=1}^k a_jx_{i-j} x1,…,xk，其中xi=∑j=1kajxi−j，求 x n x_n xn 。 ( 1 ≤ k ≤ 100 , 1 ≤ n ≤ 1 0 18 ) (1\le k \le 100, 1\le n \le 10^{18}) (1≤k≤100,1≤n≤1018)
思路：可以矩阵递推，复杂度 k 3 l o g 2 ( m ) k^3log_2(m) k3log2(m)，也可以施展BM，复杂度 k 2 l o g 2 ( m ) k^2log_2(m) k2log2(m)。
- 利用矩阵递推的时候，由公式得出的只有一项，其他的都是现有的项，通过补 0 补 1 得到。下面以 k = 5 为例：
[ 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 a 1 a 2 a 3 a 4 a 5 ] × [ x n − 5 x n − 4 x n − 3 x n − 2 x n − 1 ] = [ x n − 4 x n − 3 x n − 2 x n − 1 x n ] \begin{bmatrix} 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ a_1 & a_2 & a_3 & a_4 & a_5 \end{bmatrix} \times \begin{bmatrix} x_{n-5}\\ x_{n-4}\\ x_{n-3}\\ x_{n-2}\\ x_{n-1} \end{bmatrix}= \begin{bmatrix} x_{n-4}\\ x_{n-3}\\ x_{n-2}\\ x_{n-1}\\ x_{n} \end{bmatrix} ⎣⎢⎢⎢⎢⎡0000a11000a20100a30010a40001a5⎦⎥⎥⎥⎥⎤×⎣⎢⎢⎢⎢⎡xn−5xn−4xn−3xn−2xn−1⎦⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎡xn−4xn−3xn−2xn−1xn⎦⎥⎥⎥⎥⎤
只有 x n x_n xn 是公式得到的，其他的项都是通过 0 、1 得到的。对构造的矩阵快速幂即可。
#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int maxn=55,mod=1e9+7;
int t,n,k,a[maxn],x[maxn];
struct Matrix
{
int a[maxn][maxn];
};
Matrix mul(Matrix a,Matrix b,int len)
{
Matrix res;
memset(res.a,0,sizeof(res.a));
for(int i=1; i<=len; ++i)
for(int j=1; j<=len; ++j)
for(int k=1; k<=len; ++k)
res.a[i][j]=(1ll*res.a[i][j]+1ll*a.a[i][k]*b.a[k][j]%mod)%mod;
return res;
}
Matrix qpow(Matrix b,int n,int len)
{
Matrix res;
memset(res.a,0,sizeof(res.a));
for(int i=1; i<=len; ++i) res.a[i][i]=1;
while(n)
{
if(n&1) res=mul(res,b,len);
b=mul(b,b,len);
n>>=1;
}
return res;
}
int main()
{
scanf("%d",&t);
while(t--)
{
scanf("%d",&n);
for(int i=1; i<=n; ++i) scanf("%d",&x[i]);
for(int i=1; i<=n; ++i) scanf("%d",&a[i]);
scanf("%d",&k);
Matrix res;
memset(res.a,0,sizeof(res.a));
for(int i=1; i<=n-1; ++i) res.a[i][i+1]=1;
for(int i=1; i<=n; ++i) res.a[n][i]=a[i];
res=qpow(res,k,n);
int ans=0;
for(int i=1; i<=n; ++i) ans=(1ll*ans+1ll*res.a[n][i]*x[i]%mod)%mod;
printf("%d\n",ans);
}
return 0;
}
Matrix Power Series POJ - 3233
链接：http://poj.org/problem?id=3233
题意：给定一个 n × n n \times n n×n 的矩阵A，求 s k = A + A 2 + ⋯ + A k s_k = A+A^2+\dots +A^k sk=A+A2+⋯+Ak
思路：递推一下， s k = s k − 1 + A k s_k=s_{k-1}+A^k sk=sk−1+Ak
[ E E O A ] × [ S k − 1 A k ] = [ S k A k + 1 ] \begin{bmatrix} E&E\\ O&A \end{bmatrix} \times \begin{bmatrix} S_{k-1}\\ A^k \end{bmatrix}= \begin{bmatrix} S_k\\ A^{k+1} \end{bmatrix} [EOEA]×[Sk−1Ak]=[SkAk+1]
[ E E O A ] k − 1 × [ S 1 A 2 ] = [ S k A k + 1 ] \begin{bmatrix} E&E\\ O&A \end{bmatrix} ^{k-1} \times \begin{bmatrix} S_1\\ A^2 \end{bmatrix}= \begin{bmatrix} S_k\\ A^{k+1} \end{bmatrix} [EOEA]k−1×[S1A2]=[SkAk+1]
矩阵套矩阵写错了，直接用矩阵乘矩阵好了，还是挺麻烦的。
- 以后可以这样考虑，矩阵快速幂涉及 2 n × 2 n 2n\times 2n 2n×2n，最后一步涉及 2 n × 2 n 2n\times 2n 2n×2n 矩阵和 2 n × n 2n \times n 2n×n 的矩阵。处理好这两个矩阵乘法就可以了。
#include <cstdio>
#include <cstring>
#define ll long long
using namespace std;
const int MAXN=125;
ll n,k,mod;
struct Matrix
{
ll a[MAXN][MAXN];
int n;
void setSize(int nn)
{
n=nn;
}
void init()
{
memset(a,0,sizeof(a));
for(int i=0; i<n; ++i)
a[i][i]=1;
}
Matrix operator*(const Matrix& b) const
{
Matrix res;
res.setSize(n);
memset(res.a,0,sizeof(res.a));
for(int i=0; i<n; ++i)
for(int j=0; j<n; ++j)
for(int k=0; k<n; ++k)
res.a[i][j]=(res.a[i][j]+a[i][k]*b.a[k][j]%mod+mod)%mod;
return res;
}
};
Matrix qpow(Matrix b,int n,int size)
{
Matrix res;
res.setSize(size);
res.init();
while(n)
{
if(n&1) res=res*b;
b=b*b;
n>>=1;
}
return res;
}
Matrix Mul(Matrix a,Matrix b)
{
Matrix res;
memset(res.a,0,sizeof(res.a));
for(int i=0; i<2*n; ++i)
for(int j=0; j<n; ++j)
for(int k=0; k<2*n; ++k)
res.a[i][j]=(res.a[i][j]+a.a[i][k]*b.a[k][j]%mod)%mod;
return res;
}
int main()
{
scanf("%lld%lld%lld",&n,&k,&mod);
Matrix A,A2,A3;
A.setSize(n);
for(int i=0; i<n; ++i)
for(int j=0; j<n; ++j)
scanf("%lld",&A.a[i][j]),A.a[i][j]%=mod;
A2=A*A;
A3.setSize(2*n);
A3=A;
for(int i=0; i<n; ++i)
for(int j=0; j<n; ++j)
A3.a[i+n][j]=A2.a[i][j];
Matrix res;
res.setSize(2*n);
memset(res.a,0,sizeof(res.a));
for(int i=0; i<n; ++i)
res.a[i][i]=1;
for(int i=0; i<n; ++i)
res.a[i][i+n]=1;
for(int i=0; i<n; ++i)
for(int j=0; j<n; ++j)
res.a[i+n][j+n]=A.a[i][j];
res=qpow(res,k-1,2*n);
Matrix ans=Mul(res,A3);
for(int i=0; i<n; ++i)
for(int j=0; j<n; ++j)
printf("%lld%c",ans.a[i][j],j==n-1?'\n':' ');
return 0;
}
HDU - 2276 Kiki & Little Kiki 2
链接：http://acm.hdu.edu.cn/showproblem.php?pid=2276
题意：给定 n 盏灯顺时针排成一个环，每秒钟灯的状态会发生变化，如果每盏灯左边的灯是亮的，那么当前灯的状态就会翻转。问 t 秒钟后，每盏灯的状态。
思路： A i A_i Ai 这盏灯的状态由上一层的 A i − 1 A_{i-1} Ai−1 和 A i A_i Ai 共同决定。
| A i − 1 A_{i-1} Ai−1 | A i A_i Ai | 当前 |
|---|---|---|
| 0 | 0 | 0 |
| 0 | 1 | 1 |
| 1 | 0 | 1 |
| 1 | 1 | 0 |
由此可以得到规律：
A
i
=
(
A
i
−
1
+
A
i
)
%
2
A_i=(A_{i-1}+A_i) \%2
Ai=(Ai−1+Ai)%2 或者是
A
i
=
A
i
−
1
x
o
r
A
i
A_i=A_{i-1} \ xor \ A_i
Ai=Ai−1 xor Ai
这样就可以得到递推：
[ 1 0 0 1 1 1 0 0 0 1 1 0 0 0 1 1 ] × [ a 1 a 2 a 3 a 4 ] = [ a 1 x o r a 4 a 2 x o r a 1 a 3 x o r a 2 a 4 x o r a 3 ] \begin{bmatrix} 1 & 0 & 0 & 1\\ 1 & 1 & 0 & 0\\ 0 & 1 & 1 & 0\\ 0 & 0 & 1 & 1 \end{bmatrix} \times \begin{bmatrix} a_1\\ a_2\\ a_3\\ a_4 \end{bmatrix} = \begin{bmatrix} a_1 \ xor \ a_4\\ a_2 \ xor \ a_1\\ a_3 \ xor \ a_2\\ a_4 \ xor \ a_3 \end{bmatrix} ⎣⎢⎢⎡1100011000111001⎦⎥⎥⎤×⎣⎢⎢⎡a1a2a3a4⎦⎥⎥⎤=⎣⎢⎢⎡a1 xor a4a2 xor a1a3 xor a2a4 xor a3⎦⎥⎥⎤
在矩阵乘法中，异或表现为对结果模 2 ，也可以用 & 1 代替
#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int maxn=100+10;
int n,a[maxn];
string s;
struct Matrix
{
int a[maxn][maxn];
};
Matrix mul(Matrix a,Matrix b,int len)
{
Matrix res;
memset(res.a,0,sizeof(res.a));
for(int i=0; i<len; ++i)
for(int j=0; j<len; ++j)
for(int k=0; k<len; ++k)
res.a[i][j]=(res.a[i][j]+a.a[i][k]*b.a[k][j])&1;
return res;
}
Matrix mul2(Matrix a,Matrix b,int len)
{
Matrix res;
memset(res.a,0,sizeof(res.a));
for(int i=0; i<len; ++i)
for(int j=0; j<1; ++j)
for(int k=0; k<len; ++k)
res.a[i][j]=(res.a[i][j]+a.a[i][k]*b.a[k][j])&1;
return res;
}
Matrix qpow(Matrix b,int n,int len)
{
Matrix res;
memset(res.a,0,sizeof(res.a));
for(int i=0; i<len; ++i) res.a[i][i]=1;
while(n>0)
{
if(n&1) res=mul(res,b,len);
b=mul(b,b,len);
n>>=1;
}
return res;
}
int main()
{
while(~scanf("%d",&n))
{
cin>>s;
int len=s.size();
for(int i=0; i<len; ++i)
a[i]=s[i]-'0';
Matrix res;
memset(res.a,0,sizeof(res.a));
for(int i=0; i<len; ++i)
res.a[i][(i-1+len)%len]=res.a[i][i]=1;
res=qpow(res,n,len);
Matrix b;
memset(b.a,0,sizeof(b.a));
for(int i=0; i<len; ++i) b.a[i][0]=a[i];
Matrix ans=mul2(res,b,len);
for(int i=0; i<len; ++i)
printf("%d",ans.a[i][0]);
puts("");
}
return 0;
}
HDU - 5015 233 Matrix
链接：http://acm.hdu.edu.cn/showproblem.php?pid=5015
题意：给出矩阵的第 0 行
233
,
2333
,
…
,
2333
…
3
233,2333,\dots, 2333\dots 3
233,2333,…,2333…3（
a
0
,
1
,
a
0
,
2
…
a
0
,
m
a_{0,1},a_{0,2}\dots a_{0,m}
a0,1,a0,2…a0,m），然后给出第 0 列的 n 个数（
a
1
,
0
,
a
2
,
0
,
…
,
a
n
,
0
a_{1,0},a_{2,0},\dots,a_{n,0}
a1,0,a2,0,…,an,0）。在矩阵中
a
i
,
j
=
a
i
−
1
,
j
+
a
i
,
j
−
1
a_{i,j}=a_{i-1,j}+a_{i,j-1}
ai,j=ai−1,j+ai,j−1，问
a
n
,
m
a_{n,m}
an,m 是多少
思路：
假设给定 n =3 ，m=3 时，可构造如下矩阵
[
1
0
0
0
0
1
10
0
0
0
0
1
1
0
0
0
1
1
1
0
0
1
1
1
1
]
×
[
3
a
0
,
m
a
1
,
m
−
1
a
2
,
m
−
1
a
3
,
m
−
1
]
=
[
3
a
0
,
m
+
1
a
1
,
m
a
2
,
m
a
3
,
m
]
\begin{bmatrix} 1 & 0 & 0 & 0 &0\\ 1 & 10 & 0 & 0 & 0\\ 0 & 1 & 1 & 0 & 0\\ 0 & 1 & 1 & 1 & 0\\ 0 & 1 & 1 & 1 & 1\\ \end{bmatrix} \times \begin{bmatrix} 3\\ a_{0,m}\\ a_{1,m-1}\\ a_{2,m-1}\\ a_{3,m-1}\\ \end{bmatrix} = \begin{bmatrix} 3\\ a_{0,m+1}\\ a_{1,m}\\ a_{2,m}\\ a_{3,m}\\ \end{bmatrix}
⎣⎢⎢⎢⎢⎡11000010111001110001100001⎦⎥⎥⎥⎥⎤×⎣⎢⎢⎢⎢⎡3a0,ma1,m−1a2,m−1a3,m−1⎦⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎡3a0,m+1a1,ma2,ma3,m⎦⎥⎥⎥⎥⎤
#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int maxn=100+10,mod=10000007;
int n,m;
ll a[15][2];
struct Matrix
{
ll a[maxn][maxn];
Matrix() {memset(a,0,sizeof(a));}
};
Matrix mul(Matrix a,Matrix b,int len)
{
Matrix res;
for(int i=1; i<=len; ++i)
for(int j=1; j<=len; ++j)
for(int k=1; k<=len; ++k)
res.a[i][j]=(res.a[i][j]+1ll*a.a[i][k]*b.a[k][j]%mod)%mod;
return res;
}
Matrix qpow(Matrix b,int n,int len)
{
Matrix res;
for(int i=1; i<=len; ++i) res.a[i][i]=1;
while(n>0)
{
if(n&1) res=mul(res,b,len);
b=mul(b,b,len);
n>>=1;
}
return res;
}
int main()
{
while(~scanf("%d%d",&n,&m))
{
a[1][1]=3;
a[2][1]=233;
for(int i=1; i<=n; ++i)
scanf("%lld",&a[i+2][1]);
Matrix res;
res.a[1][1]=1,res.a[2][1]=1,res.a[2][2]=10;
for(int i=3; i<=n+2; ++i) res.a[i][2]=1;
for(int i=3; i<=n+2; ++i)
for(int j=3; j<=i; ++j)
res.a[i][j]=1;
res=qpow(res,m,n+2);
ll ans=0;
for(int i=1; i<=n+2; ++i)
ans=(ans+1ll*res.a[n+2][i]*a[i][1]%mod)%mod;
printf("%lld\n",ans);
}
return 0;
}
CF1182 E. Product Oriented Recurrence
链接：https://codeforces.com/contest/1182/problem/E
题意： f ( x ) = c 2 x − 6 ⋅ f ( x − 1 ) ⋅ f ( x − 2 ) ⋅ f ( x − 3 ) f(x)=c^{2x−6}⋅f(x−1)⋅f(x−2)⋅f(x−3) f(x)=c2x−6⋅f(x−1)⋅f(x−2)⋅f(x−3)，求 f ( n ) m o d ( 1 e 9 + 7 ) f(n)\ mod (1e9+7) f(n) mod(1e9+7)
思路：矩阵快速幂可以解决线性递推关系。虽然我们不能直接求
f
(
n
)
f(n)
f(n)的递推关系，但我们可以把题目转变成求有多少个
f
(
1
)
,
f
(
2
)
,
f
(
3
)
,
c
f(1),f(2),f(3),c
f(1),f(2),f(3),c，也就是说可以把题意转化成求
f
(
1
)
,
f
(
2
)
,
f
(
3
)
,
c
f(1),f(2),f(3),c
f(1),f(2),f(3),c 的指数。设
f
1
(
n
)
,
f
2
(
n
)
,
f
3
(
n
)
f_1(n),f_2(n),f_3(n)
f1(n),f2(n),f3(n)分别表示
f
(
n
)
f(n)
f(n)中
f
(
1
)
,
f
(
2
)
,
f
(
3
)
f(1),f(2),f(3)
f(1),f(2),f(3)的个数。因此可以得到线性关系
f
1
(
n
)
=
f
1
(
n
−
1
)
+
f
1
(
n
−
2
)
+
f
1
(
n
−
3
)
f_1(n)=f_1(n-1)+f_1(n-2)+f_1(n-3)
f1(n)=f1(n−1)+f1(n−2)+f1(n−3)从而得到转移矩阵
[
1
1
1
1
0
0
0
1
0
]
∗
[
f
1
(
n
−
1
)
f
1
(
n
−
2
)
f
1
(
n
−
3
)
]
=
[
f
1
(
n
)
f
1
(
n
−
1
)
f
1
(
n
−
2
)
]
\begin{bmatrix} 1&1&1\\ 1&0&0\\ 0&1&0 \end{bmatrix} *\begin{bmatrix} f_1(n-1)\\ f_1(n-2)\\ f_1(n-3) \end{bmatrix}= \begin{bmatrix} f_1(n)\\ f_1(n-1)\\ f_1(n-2) \end{bmatrix}
⎣⎡110101100⎦⎤∗⎣⎡f1(n−1)f1(n−2)f1(n−3)⎦⎤=⎣⎡f1(n)f1(n−1)f1(n−2)⎦⎤
变形之后，因为
n
≥
4
，
并
且
f
1
(
3
)
=
0
,
f
1
(
2
)
=
0
,
f
1
(
1
)
=
1
n\ge 4，并且 f_1(3)=0,f_1(2)=0,f_1(1)=1
n≥4，并且f1(3)=0,f1(2)=0,f1(1)=1，因此
f
1
(
n
)
的
指
数
就
是
f_1(n)的指数就是
f1(n)的指数就是
T
n
−
3
[
1
]
[
3
]
T^{n-3}[1][3]
Tn−3[1][3]
[
1
1
1
1
0
0
0
1
0
]
n
−
3
∗
[
f
1
(
3
)
f
1
(
2
)
f
1
(
1
)
]
=
[
f
1
(
n
)
f
1
(
n
−
1
)
f
1
(
n
−
2
)
]
\begin{bmatrix} 1&1&1\\ 1&0&0\\ 0&1&0 \end{bmatrix}^{n-3} *\begin{bmatrix} f_1(3)\\ f_1(2)\\ f_1(1) \end{bmatrix}= \begin{bmatrix} f_1(n)\\ f_1(n-1)\\ f_1(n-2) \end{bmatrix}
⎣⎡110101100⎦⎤n−3∗⎣⎡f1(3)f1(2)f1(1)⎦⎤=⎣⎡f1(n)f1(n−1)f1(n−2)⎦⎤
同理可得
f
2
(
n
)
和
f
3
(
n
)
的
值
f_2(n)和f_3(n)的值
f2(n)和f3(n)的值
而c的指数同样可以得到线性关系
g
(
n
)
=
g
(
n
−
1
)
+
g
(
n
−
2
)
+
g
(
n
−
3
)
+
2
n
−
6
g(n)=g(n-1)+g(n-2)+g(n-3)+2n-6
g(n)=g(n−1)+g(n−2)+g(n−3)+2n−6转移矩阵为
[
1
1
1
2
−
4
1
0
0
0
0
0
1
0
0
0
0
0
0
1
1
0
0
0
0
1
]
∗
[
g
(
n
−
1
)
g
(
n
−
2
)
g
(
n
−
3
)
n
−
1
1
]
=
[
g
(
n
)
g
(
n
−
1
)
g
(
n
−
2
)
n
1
]
\begin{bmatrix} 1&1&1&2&-4\\ 1&0&0&0&0\\ 0&1&0&0&0\\ 0&0&0&1&1\\ 0&0&0&0&1 \end{bmatrix} * \begin{bmatrix} g(n-1)\\ g(n-2)\\ g(n-3)\\ n-1\\ 1 \end{bmatrix}= \begin{bmatrix} g(n)\\ g(n-1)\\ g(n-2)\\ n\\ 1 \end{bmatrix}
⎣⎢⎢⎢⎢⎡11000101001000020010−40011⎦⎥⎥⎥⎥⎤∗⎣⎢⎢⎢⎢⎡g(n−1)g(n−2)g(n−3)n−11⎦⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎡g(n)g(n−1)g(n−2)n1⎦⎥⎥⎥⎥⎤
我们最终的所求就是
f
(
n
)
=
c
a
1
⋅
f
(
1
)
a
2
⋅
f
(
2
)
a
3
⋅
f
(
3
)
a
4
m
o
d
1
e
9
+
7
f(n)=c^{a_1} \cdot f(1)^{a_2} \cdot f(2)^{a_3} \cdot f(3)^{a_4} mod 1e9+7
f(n)=ca1⋅f(1)a2⋅f(2)a3⋅f(3)a4mod1e9+7
根据欧拉降幂公式可知，指数是关于
φ
(
p
)
\varphi(p)
φ(p)循环的，所以在求它们的指数的时候，它们的模数应该是
φ
(
p
)
=
1
e
9
+
6
\varphi(p)=1e9+6
φ(p)=1e9+6
#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int maxn=1e5+5,mod=1e9+6,Mod=1e9+7;
ll n,f1,f2,f3,c;
struct Node
{
ll a[6][6];
} g;
Node mul(Node A,Node B,int len)
{
Node res;
memset(res.a,0,sizeof(res.a));
for(int i=1; i<=len; ++i)
for(int j=1; j<=len; ++j)
for(int k=1; k<=len; ++k)
res.a[i][j]=(1ll*res.a[i][j]+1ll*A.a[i][k]*B.a[k][j]%mod+mod)%mod;
return res;
}
Node qpow(Node b,ll n,int len)
{
Node res;
memset(res.a,0,sizeof(res.a));
for(int i=1; i<=len; ++i) res.a[i][i]=1;
while(n)
{
if(n&1) res=mul(res,b,len);
b=mul(b,b,len);
n>>=1;
}
return res;
}
int qpow(int b,int n,int mod)
{
int res=1;
while(n)
{
if(n&1) res=1ll*res*b%mod;
b=1ll*b*b%mod;
n>>=1;
}
return res;
}
int main()
{
cin>>n>>f1>>f2>>f3>>c;
Node A,B;
memset(A.a,0,sizeof(A.a));
A.a[1][1]=1,A.a[1][2]=1,A.a[1][3]=1;
A.a[2][1]=1,A.a[2][2]=0,A.a[2][3]=0;
A.a[3][1]=0,A.a[3][2]=1,A.a[3][3]=0;
B=qpow(A,n-3,3);
int exp1=B.a[1][3],exp2=B.a[1][2],exp3=B.a[1][1];
int f[10]= {0,f3,f2,f1,4,1};
g.a[1][1]=1, g.a[1][2]=1, g.a[1][3]=1, g.a[1][4]=2, g.a[1][5]=-6;
g.a[2][1]=1, g.a[2][2]=0, g.a[2][3]=0, g.a[2][4]=0, g.a[2][5]=0;
g.a[3][1]=0, g.a[3][2]=1, g.a[3][3]=0, g.a[3][4]=0, g.a[3][5]=0;
g.a[4][1]=0, g.a[4][2]=0, g.a[4][3]=0, g.a[4][4]=1, g.a[4][5]=1;
g.a[5][1]=0, g.a[5][2]=0, g.a[5][4]=0, g.a[5][4]=0, g.a[5][5]=1;
g=qpow(g,n-3,5);
int expc=(1ll*g.a[1][4]*4+g.a[1][5]+mod)%mod;
f1=qpow(f1,exp1,Mod);
f2=qpow(f2,exp2,Mod);
f3=qpow(f3,exp3,Mod);
int ans=1ll*qpow(c,expc,Mod)*f1%Mod*f2%Mod*f3%Mod;
printf("%d\n",ans);
return 0;
}
扫一扫
987
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
