已知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;}新聞熱點
疑難解答
