NYOJ 301 遞推求值【矩陣快速冪取模】
阿新 • • 發佈：2018-12-24
遞推求值
時間限制：1000 ms | 記憶體限制：65535 KB 難度：4- 描述
-
給你一個遞推公式：
f(x)=a*f(x-2)+b*f(x-1)+c
並給你f(1),f(2)的值，請求出f(n)的值，由於f(n)的值可能過大，求出f(n)對1000007取模後的值。
注意:-1對3取模後等於2
分析：由於n的值比較大，所以常規方法肯定會超時。根據遞推式求第n個表示式的值時，通常用矩陣乘法來做。
本題要構造兩個矩陣，其中一個為矩陣A，作為初始矩陣f2 0 0
f1 0 0
1 0 0
另一個為矩陣B
b a c
1 0 0
0 0 1
因為F(2)和F(1)是已知的，當n>=3時，每次都乘以矩陣B，就能推出下一個矩陣。而矩陣的第一行第一列的元素就是所求的結果。
所以利用矩陣快速冪能夠快速準確地求出結果。
//B^(n-2)採用矩陣快速冪計算
AC程式碼：
#include <iostream> #include <cstdio> #include <cstring> using namespace std; typedef long long LL; const int maxn=100+5; const int mod=1000007; const int N=3; struct Mat { LL mat[N][N]; }; Mat mul(Mat a,Mat b) { Mat ans; for(int i=0; i<N; i++) { for(int j=0; j<N; j++) { ans.mat[i][j]=0; for(int k=0; k<N; k++) { ans.mat[i][j]+=a.mat[i][k]*b.mat[k][j]; ans.mat[i][j]%=mod; } } } return ans; } Mat quickPow(Mat a,int n) { Mat ans= { 1,0,0, 0,1,0, 0,0,1 }; while(n) { if(n&1) ans=mul(ans,a); n>>=1; a=mul(a,a); } return ans; } int main() { int T; cin>>T; int a,b,c,n; int f[2]; Mat tmp,ans; Mat e= { 0,0,0, 1,0,0, 0,0,1 }; tmp.mat[2][0]=1; while(T--) { cin>>f[0]>>f[1]>>a>>b>>c>>n; if(n==1||n==2) { cout<<((f[n-1]+mod)%mod)<<endl; continue; } e.mat[0][0]=b; e.mat[0][1]=a; e.mat[0][2]=c; tmp.mat[0][0]=f[1]; tmp.mat[1][0]=f[0]; // ans=mul(quickPow(e,n-2),tmp); cout<<(ans.mat[0][0]+mod)%mod<<endl; } return 0; }
NYOJ排行第一程式碼：
執行號：548458 執行人:張燚
#include<stdio.h>
long long m,n,a,b,c,k;
long long mol=1000007;
long long x1,x2,x3;
long long y1,y2,y3;
long long z1,z2,z3;
void f()
{
long long h=k;
if(h==1)
{
printf("%lld\n",m);
return;
}
if(h==2)
{
printf("%lld\n",n);
return;
}
long long r1=1,r2=0,r3=0;
long long s1=0,s2=1,s3=0;
long long t1=0,t2=0,t3=1;
h-=2;
while(h>=1)
{
int o1,o2,o3,p1,p2,p3,q1,q2,q3;
if(h%2==1)
{
o1=(r1*x1+s1*x2+t1*x3)%mol;
p1=(r1*y1+s1*y2+t1*y3)%mol;
q1=(r1*z1+s1*z2+t1*z3)%mol;
o2=(r2*x1+s2*x2+t2*x3)%mol;
p2=(r2*y1+s2*y2+t2*y3)%mol;
q2=(r2*z1+s2*z2+t2*z3)%mol;
o3=(r3*x1+s3*x2+t3*x3)%mol;
p3=(r3*y1+s3*y2+t3*y3)%mol;
q3=(r3*z1+s3*z2+t3*z3)%mol;
r1=o1,r2=o2,r3=o3;
s1=p1,s2=p2,s3=p3;
t1=q1,t2=q2,t3=q3;
}
o1=(x1*x1+y1*x2+z1*x3)%mol;
p1=(x1*y1+y1*y2+z1*y3)%mol;
q1=(x1*z1+y1*z2+z1*z3)%mol;
o2=(x2*x1+y2*x2+z2*x3)%mol;
p2=(x2*y1+y2*y2+z2*y3)%mol;
q2=(x2*z1+y2*z2+z2*z3)%mol;
o3=(x3*x1+y3*x2+z3*x3)%mol;
p3=(x3*y1+y3*y2+z3*y3)%mol;
q3=(x3*z1+y3*z2+z3*z3)%mol;
x1=o1,x2=o2,x3=o3;
y1=p1,y2=p2,y3=p3;
z1=q1,z2=q2,z3=q3;
h/=2;
}
n=(m*r2+n*s2+t2)%mol;
n=(n+mol)%mol;
printf("%lld\n",n);
}
int main()
{
int N;
scanf("%d",&N);
while(N--)
{
scanf("%lld%lld%lld%lld%lld%lld",&m,&n,&a,&b,&c,&k);
x1=0,y1=1,z1=0;
x2=a,y2=b,z2=c;
x3=0,y3=0,z3=1;
f();
}
}
標程程式碼：
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<iostream>
using namespace std;
#define CLR(arr,val) memset(arr,val,sizeof(arr))
const int MAX=100;
char str[MAX+2];
struct Matrix;
int mod=1000007;
struct Matrix
{
Matrix() {}
void init(int r,int c) //全零矩陣
{
row=r;
col=c;
CLR(data,0);
}
void init(int size)
{
row=size;
col=size;
CLR(data,0);
for(int i=0; i!=size; i++) data[i][i]=1;
}
int data[MAX][MAX];
int row,col;
const Matrix& pow(int m);
};
Matrix ret,temp,pret;
const Matrix& operator*(const Matrix& m1,const Matrix &m2)
{
ret.init(m1.row,m2.col);
for(int i=0; i!=m1.row; i++)
for(int j=0; j!=m1.col; j++)
if(m1.data[i][j]) for(int k=0; k!=m2.col; k++)
if(m2.data[j][k]) ret.data[i][k]=(ret.data[i][k]+((long long)m1.data[i][j]*m2.data[j][k]))%mod;
return ret;
}
const Matrix& Matrix::pow(int m)
{
temp=*this;
pret.init(row);
while(m)
{
if(m&1)
pret=temp*pret;
temp=temp*temp;
m>>=1;
}
return pret;
}
Matrix m1,m2;
int main()
{
int t,k;
m1.init(3,3);
m1.data[1][0]=1;
m1.data[2][2]=1;
m2.init(1,3);
m2.data[0][2]=1;
scanf("%d",&t);
while(t--)
{
scanf("%d%d%d%d%d%d",&m2.data[0][0],&m2.data[0][1],&m1.data[0][1],&m1.data[1][1],&m1.data[2][1],&k);
if(k==1) printf("%d\n",m2.data[0][0]);
else printf("%d\n",((m2*m1.pow(k-2)).data[0][1]+mod)%mod);
}
}
AC程式碼3：
#include <stdio.h>
#include <string.h>
#define N 3
#define mod 1000007
#define ll long long
struct Matrix{
ll mat[N][N];//乘積時(可能接近mod^2)超出int,
};
struct Matrix mul(struct Matrix a,struct Matrix b)
{
int i,j,k;
struct Matrix res;
for(i=0;i<N;i++)
{
for(j=0;j<N;j++)
{
res.mat[i][j] = 0;
for(k=0;k<N;k++)
{
res.mat[i][j] += a.mat[i][k]*b.mat[k][j];
res.mat[i][j] %= mod;//必須每次取餘
}
}
}
return res;
}
struct Matrix mul_matrix(struct Matrix b,int n)
{
struct Matrix res = {
1,0,0,
0,1,0,
0,0,1
};//單位陣
while(n)
{
if(n&1)
res = mul(res,b);
n >>= 1;
b = mul(b,b);
}
return res;
}
int main()
{
int T,f[2],a,b,c,n;
struct Matrix tmp,res;
struct Matrix tmp_matrix = {//中間矩陣
0,0,0,//每組測試將此三處的值更新為b,a,c
1,0,0,
0,0,1
};
memset(tmp.mat,0,sizeof(tmp.mat));//tmp值只有(0,0),(1,0)需要填入f(2),f(1)的值,其餘不變
tmp.mat[2][0] = 1;
scanf("%d",&T);
while(T--)
{
scanf("%d%d%d%d%d%d",&f[0],&f[1],&a,&b,&c,&n);
if(n == 1 || n == 2)
printf("%d\n",(f[n-1]+mod)%mod);//注意對負值的處理
else{
tmp_matrix.mat[0][0] = b;
tmp_matrix.mat[0][1] = a;
tmp_matrix.mat[0][2] = c;
tmp.mat[0][0] = f[1];
tmp.mat[1][0] = f[0];
//實際上最終我們只需進行矩陣和向量(f(2),f(1),1)'的運算,但是因為在進行矩陣快速冪時我們已經定義了矩陣乘積
//所以不妨借用,最後取(0,0)的值即可
res = mul(mul_matrix(tmp_matrix,n-2),tmp);
printf("%lld\n",(res.mat[0][0]+mod)%mod);
}
}
return 0;
} 尊重原創，轉載請註明出處：