您的位置:首页 > 其它

中国(北方)大学生程序设计训练赛(第一周)(Problem E: Water Problem-矩阵快速幂)

2017-03-07 16:06 871 查看
已知f(1),f(2),n,

f(n+1)=f(n)+f(n-1)+sin(n*Pi/2),(n>=2)

求f(n)

矩阵快速幂,周期乘4个矩阵

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<algorithm>
#include<functional>
#include<iostream>
#include<cmath>
#include<cctype>
#include<ctime>
using namespace std;
#define For(i,n) for(int i=1;i<=n;i++)
#define Fork(i,k,n) for(int i=k;i<=n;i++)
#define ForkD(i,k,n) for(int i=n;i>=k;i--)
#define Rep(i,n) for(int i=0;i<n;i++)
#define ForD(i,n) for(int i=n;i;i--)
#define RepD(i,n) for(int i=n;i>=0;i--)
#define Forp(x) for(int p=pre[x];p;p=next[p])
#define Forpiter(x) for(int &p=iter[x];p;p=next[p])
#define Lson (x<<1)
#define Rson ((x<<1)+1)
#define MEM(a) memset(a,0,sizeof(a));
#define MEMI(a) memset(a,127,sizeof(a));
#define MEMi(a) memset(a,128,sizeof(a));
#define INF (2139062143)
#define F (1000000007)
#define eps (1e-3)
#define MAXN (5+10)
#define MAXM (5+10)
typedef long long ll;
ll mul(ll a,ll b){return (a*b)%F;}
ll add(ll a,ll b){return (a+b)%F;}
ll sub(ll a,ll b){return (a-b+llabs(a-b)/F*F+F)%F;}
void upd(ll &a,ll b){a=(a%F+b%F)%F;}

struct M
{
int n,m;
ll a[MAXN][MAXN];
M(int _n=0){n=m=_n;MEM(a);}
M(int _n,int _m){n=_n,m=_m;MEM(a);}
void mem (int _n=0){n=m=_n;MEM(a);}
void mem (int _n,int _m){n=_n,m=_m;MEM(a);}
void fibM(int t) {
n=m=3; MEM(a)
a[1][1]=a[1][2]=a[2][1]=a[3][3]=1;;
a[1][3]=t;
}
friend M operator*(M a,M b)
{
M c(a.n,b.m);
For(k,a.m)
For(i,a.n)
For(j,b.m)
c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j])%F;
return c;
}
friend M operator+(M a,M b)
{
For(i,a.n)
For(j,a.m)
a.a[i][j]=(a.a[i][j]+b.a[i][j])%F;
return a;
}
void make_I(int _n)
{
n=m=_n; MEM(a)
For(i,n) a[i][i]=1;
}
}A,B,C,D,I,ans,E;
M pow2(M a,ll b)
{
M c;c.make_I(a.n);
static bool a2[1000000];
int n=0;while (b) a2[++n]=b&1,b>>=1;
For(i,n)
{
if (a2[i]) c=c*a;
a=a*a;
}
return c;
}
int main()
{
// freopen("E.in","r",stdin);

A.fibM(1);
B.fibM(0);
C.fibM(-1);
D=((A*B)*C)*B;
ll f1,f2,n;
while(scanf("%lld%lld%lld",&f1,&f2,&n)==3) {
if (n==1) {
printf("%lld\n",f1);
} else if (n==2) {
printf("%lld\n",f2);
} else {
ll c=n-2;
ll p=c/4;
I.make_I(3);
if (p) I=pow2(D,p);
c%=4;
if (c) {
I=B*I;c--;
}if (c) {
I=C*I;c--;
}if (c) {
I=B*I;c--;
}
ll ans=0;
upd(ans,mul(I.a[1][1],f2));
upd(ans,mul(I.a[1][2],f1));
upd(ans,I.a[1][3]);
ans=sub(ans,0);
printf("%lld\n",ans);
}

}
return 0;
}
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: 
相关文章推荐