矩阵快速幂
矩阵模板：
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mod = 1e9 + 7;
const int NUM = 7; // 最大存储空间
int N; //实际进行运算的方阵的维度N
struct mat{
int a[NUM][NUM];
void clear(){
memset(a,0,sizeof(a));
}
void init(){ //初始化成单位方阵
memset(a,0,sizeof(a));
for(int i = 0; i < N; ++i) a[i][i] = 1;
}
void show() const{ //输出矩阵
for(int i = 0; i < N; ++i){
for(int j = 0; j < N; ++j){
printf("%d%c",a[i][j],(j == N - 1) ? '\n' : ' ');
}
}
}
mat operator * (const mat &b) const{
mat c;
for(int i = 0; i < N; ++i){
for(int j = 0; j < N; ++j){
c.a[i][j] = 0;
for(int k = 0; k < N; ++k){
c.a[i][j] = (c.a[i][j] + 1LL * a[i][k] * b.a[k][j] % mod) % mod; //注意中间过程精度
}
}
}
return c;
}
mat operator + (const mat &b) const{
mat c;
for(int i = 0; i < N; ++i){
for(int j = 0;j < N; ++j){
c.a[i][j] = ((ll)a[i][j] + b.a[i][j]) % mod;
}
}
return c;
}
};
mat fast_pow(mat a,ll b){ //求矩阵a ^ n幂次方之和,注意幂次，int 还是 long long
mat res;
res.init();
while(b > 0){
if(b & 1) res = res * a;
a = a * a;
b >>= 1;
}
return res;
}
mat pow_sum(mat a,ll n){ //(a+a^2+a^3....+a^n)%mod 矩阵的幂和,也需要注意幂次，int or long long
if(n == 1) return a;
ll mid = n >> 1;
mat pre = pow_sum(a,mid);
mat ans = pre + pre * fast_pow(a,mid);
if(n & 1)
ans = ans + fast_pow(a,mid);
return ans;
}
int det(mat a){ // 求方阵的行列式
int ans = 1,sign = 0;
for(int i = 0; i < N; ++i){
for(int j = i + 1; j < N; ++j){
int x = i, y = j;
while(a.a[y][i]){
int t = a.a[x][i] / a.a[y][i];
for(int k = i; k < N; ++k){
a.a[x][k] = (a.a[x][k] - a.a[y][k] * t) % mod;
}
swap(x,y);
}
if(x != i){
for(int k = 0; k < N; ++k){
swap(a.a[x][k],a.a[i][k]);
}
sign ^= 1;
}
}
if(a.a[i][i] == 0){
ans = 0;
return ans;
}else{
ans = ans * a.a[i][i] % mod;
}
}
if(sign) ans *= -1;
if(ans < 0) ans = (ans % mod + mod) % mod;
return ans;
}
int main(){
mat a;
int n;
while(cin >> n){
N = n;
for(int i = 0; i < n; ++i){
for(int j = 0; j < n; ++j){
cin >> a.a[i][j];
}
}
a.show();
cout << det(a) << endl;
mat tmp = fast_pow(a, 2);
tmp.show();
mat b = pow_sum(a,2);
b.show();
}
return 0;
}
矩阵快速幂主要用于递推：
第一步推出递推关系很重要，然后是构建变换矩阵：
先举个简单栗子：
再说下经典的例题：
cdoe：
#include<iostream>
#include<algorithm>
#include<cstdio>
#include<cstring>
using namespace std;
typedef long long ll;
const int mod = 10000;
const int NUM = 2;
int N = 2;
struct mat{
int a[NUM][NUM];
void clear(){
memset(a,0,sizeof(a));
}
void init(){
memset(a,0,sizeof(a));
for(int i = 0; i < N; ++i) a[i][i] = 1;
}
void show() const{
for(int i = 0; i < N; ++i){
for(int j = 0; j < N; ++j){
printf("%d%c",a[i][j],(j == N - 1) ? '\n' : ' ');
}
}
}
mat operator * (const mat &b) const{
mat c;
for(int i = 0; i < N; ++i){
for(int j = 0; j < N; ++j){
c.a[i][j] = 0;
for(int k = 0; k < N; ++k){
c.a[i][j] = (c.a[i][j] + 1LL * a[i][k] * b.a[k][j] % mod) % mod; //注意中间过程精度
}
}
}
return c;
}
mat operator + (const mat &b) const{
mat c;
for(int i = 0; i < N; ++i){
for(int j = 0;j < N; ++j){
c.a[i][j] = ((ll)a[i][j] + b.a[i][j]) % mod;
}
}
return c;
}
};
mat fast_pow(mat a,ll b){
mat res;
res.init();
while(b > 0){
if(b & 1) res = res * a;
a = a * a;
b >>= 1;
}
return res;
}
mat pow_sum(mat a,ll n){
if(n == 1) return a;
ll mid = n >> 1;
mat pre = pow_sum(a,mid);
mat ans = pre + pre * fast_pow(a,mid);
if(n & 1)
ans = ans + fast_pow(a,mid);
return ans;
}
int det(mat a){
int ans = 1,sign = 0;
for(int i = 0; i < N; ++i){
for(int j = i + 1; j < N; ++j){
int x = i, y = j;
while(a.a[y][i]){
int t = a.a[x][i] / a.a[y][i];
for(int k = i; k < N; ++k){
a.a[x][k] = (a.a[x][k] - a.a[y][k] * t) % mod;
}
swap(x,y);
}
if(x != i){
for(int k = 0; k < N; ++k){
swap(a.a[x][k],a.a[i][k]);
}
sign ^= 1;
}
}
if(a.a[i][i] == 0){
ans = 0;
return ans;
}else{
ans = ans * a.a[i][i] % mod;
}
}
if(sign) ans *= -1;
if(ans < 0) ans = (ans % mod + mod) % mod;
return ans;
}
int main(){
mat a;
a.a[0][0] = 1, a.a[0][1] = 1;
a.a[1][0] = 1, a.a[1][1] = 0;
ll n;
int f[]={0,1,1,2,3};
while(cin >> n){
if(n == -1) break ;
if(n <= 2){
cout << f[n] << endl;
continue ;
}
mat F;
F = fast_pow(a,n - 2);
printf("%d\n",(F.a[0][0] + F.a[0][1]) % mod);
}
return 0;
}
举例说：
n ^ k = (n - 1 + 1) ^ k二项式展开获得系数Cm n;
举两道例题：
数列计算
问题描述
Cai 给了如下递归定义的序列：
Fn = 3Fn-1 +2 Fn-2 +n * 3 ^n，（n>1)，F1=2， F0 =1
他要求你对给定一个整数n，求Fn 是多少？
输入
多组数据，每组一行，其上有一个整数n(0<n<=10000000000000)。
输出
每组输出一行，输出Fn (mod 100000007)
输入样例
3
1000
输出样例
55
93320452
思路：
最后注意：n加多1，因为板子里面A ^ 0 是返回的是单位矩阵，但推式子的时候认为我们认为 A ^ 0 = A，所以要多加一次幂。
AC code：
#include<bits/stdc++.h>
using namespace std;
const int NUM = 4;
const int mod = 100000007;
typedef long long ll;
int N = 4;
struct mat{
int a[NUM][NUM];
void clear(){
memset(a,0,sizeof(a));
}
void init(){
memset(a,0,sizeof(a));
for(int i = 0; i < N; ++i) a[i][i] = 1;
}
void show(){
for(int i = 0; i < N; ++i){
for(int j = 0; j < N; ++j){
printf("%d%c",a[i][j], (j == N - 1) ? '\n' : ' ');
}
}
}
mat operator * (const mat &b) const{
mat c;
for(int i = 0; i < N; ++i){
for(int j = 0; j < N; ++j){
c.a[i][j] = 0;
for(int k = 0; k < N; ++k){
c.a[i][j] = (c.a[i][j] + 1LL * a[i][k] * b.a[k][j] % mod) % mod;
}
}
}
return c;
}
};
mat fast_pow(mat a,ll b){
mat res;
res.init();
while(b > 0){
if(b & 1) res = res * a;
a = a * a;
b >>= 1;
}
return res;
}
int main(){
mat A;
A.clear();
A.a[0][0] = 3, A.a[0][1] = 2, A.a[0][2] = 3, A.a[0][3] = 3;
A.a[1][0] = 1;
A.a[2][2] = 3, A.a[2][3] = 3;
A.a[3][3] = 3;
ll n;
while(cin >> n){
if(n == 0){
cout << 1 <<endl;
continue ;
}
if(n == 1){
cout << 2 <<endl;
continue ;
}
n += 1; //注意加多1，因为板子里面A^0是单位矩阵，但推式子的时候认为A^0=A，所以要多加一次幂
mat F;
F = fast_pow(A,n - 2);
int ans = 0;
ans = (ans + F.a[0][0] * 2LL) % mod;
ans = (ans + F.a[0][1] * 1LL) % mod;
ans = (ans + F.a[0][2] * 3LL) % mod;
ans = (ans + F.a[0][3] * 3LL) % mod;
cout << ans << endl;
}
return 0;
}
[HDU5950 Recursive sequence]
题意：
求 f(n) = f(n−1) + 2 f∗(n−2) + n^4，其中 f(1) = a, f(2) = b;
思路：
情形3情况，还是套路题。
AC code：
#include<bits/stdc++.h>
using namespace std;
const int NUM = 8;
const long long mod = 2147493647;
typedef long long ll;
int N = 7;
struct mat{
ll a[NUM][NUM];
void init(){
memset(a,0,sizeof(a));
for(int i = 0; i < N; ++i) a[i][i] = 1;
}
mat operator * (const mat &b) const{
mat c;
for(int i = 0; i < N; ++i){
for(int j = 0; j < N; ++j){
c.a[i][j] = 0;
for(int k = 0; k < N; ++k){
c.a[i][j] = (c.a[i][j] + a[i][k] * b.a[k][j] % mod) % mod;
}
}
}
return c;
}
};
mat fast_pow(mat a,ll b){
mat res;
res.init();
while(b > 0){
if(b & 1) res = res * a;
a = a * a;
b >>= 1;
}
return res;
}
int main(){
mat ma;
memset(ma.a,0,sizeof(ma.a));
ma.a[0][0] = 1, ma.a[0][1] = 2, ma.a[0][2] = 1,ma.a[0][3] = 4,ma.a[0][4] = 6;
ma.a[0][5] = 4, ma.a[0][6] = 1;
ma.a[1][0] = 1;
ma.a[2][2] = 1, ma.a[2][3] = 4, ma.a[2][4] = 6, ma.a[2][5] = 4, ma.a[2][6] = 1;
ma.a[3][3] = 1, ma.a[3][4] = 3, ma.a[3][5] = 3, ma.a[3][6] = 1;
ma.a[4][4] = 1, ma.a[4][5] = 2, ma.a[4][6] = 1;
ma.a[5][5] = 1, ma.a[5][6] = 1;
ma.a[6][6] = 1;
int T;
ll n,a,b;
cin >> T;
while(T--){
cin >> n >> a >> b;
if(n == 1){
cout << a % mod << endl;
continue ;
}
if(n == 2){
cout << b % mod << endl;
continue ;
}
mat tmp;
tmp = fast_pow(ma,n - 2);
ll ans = 0;
ans = (ans + tmp.a[0][0] * b) % mod;
ans = (ans + tmp.a[0][1] * a) % mod;
ans = (ans + tmp.a[0][2] * 16) % mod;
ans = (ans + tmp.a[0][3] * 8) % mod;
ans = (ans + tmp.a[0][4] * 4) % mod;
ans = (ans + tmp.a[0][5] * 2) % mod;
ans = (ans + tmp.a[0][6]) % mod;
cout << ans << endl;
}
return 0;
}
488
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
