hdu3306：Another kind of Fibonacci

A(0)=A(1)=1，A(i)=X*A(i-1)+Y*A(i-2)，求S(n)=A(0)^2+A(1)^2+A(2)^2+A(3)^2+……+A(n)^2。

这个矩阵有点毒。。

 1 #include<stdio.h>
 2 #include<string.h>
 3 #include<algorithm>
 4 #include<math.h>
 5 //#include<iostream>
 6 using namespace std;
 7
 8 #define LL long long
 9 LL n,x,y;
10 typedef LL mat[6][6];
11 mat a,f,ans;
12 const int mod=10007;
13 void copy(mat &a,mat b)
14 {
15     for (int i=1;i<=4;i++)
16         for (int j=1;j<=4;j++)
17             a[i][j]=b[i][j];
18 }
19 void mul(mat a,mat b,mat &ans)
20 {
21     mat t;memset(t,0,sizeof(t));
22     for (int i=1;i<=4;i++)
23         for (int j=1;j<=4;j++)
24             for (int k=1;k<=4;k++)
25                 t[i][j]=(t[i][j]+a[i][k]*b[k][j]%mod)%mod;
26     copy(ans,t);
27 }
28 void init(mat &a)
29 {
30     memset(a,0,sizeof(a));
31     for (int i=1;i<=4;i++) a[i][i]=1;
32 }
33 void pow(mat a,LL b,mat &ans)
34 {
35     mat tmp,last;copy(tmp,a);init(last);
36     while (b)
37     {
38         if (b&1) mul(last,tmp,last);
39         mul(tmp,tmp,tmp);
40         b>>=1;
41     }
42     copy(ans,last);
43 }
44 int main()
45 {
46     while (~scanf("%lld%lld%lld",&n,&x,&y))
47     {
48         memset(a,0,sizeof(a));
49         a[1][1]=a[1][2]=a[3][2]=1;
50         a[2][2]=x*x%mod;
51         a[2][3]=y*y%mod;
52         a[2][4]=2*x*y%mod;
53         a[4][2]=x;a[4][4]=y;
54         f[1][1]=f[2][1]=f[3][1]=f[4][1]=1;
55         pow(a,n,a);
56         mul(a,f,ans);
57         printf("%lld\n",ans[1][1]);
58     }
59     return 0;
60 }