hdu 4686 Arc of Dream（矩阵快速幂乘法）

Problem Description

An Arc of Dream is a curve defined by following function:

where
a0 = A0
ai = ai-1*AX+AY
b0 = B0
bi = bi-1*BX+BY
What is the value of AoD(N) modulo 1,000,000,007?

Input

There are multiple test cases. Process to the End of File.
Each test case contains 7 nonnegative integers as follows:
N
A0 AX AY
B0 BX BY
N is no more than 1018, and all the other integers are no more than 2×109.

Output

For each test case, output AoD(N) modulo 1,000,000,007.

Sample Input

Sample Output

4 
134 
1902

Author

Zejun Wu (watashi)

Source

2013 Multi-University Training Contest 9

因为：a[i]*b[i]=(a[i-1]*AX+AY)*(b[i-1]*BX+BY)

=(a[i-1]*b[i-1]*AX*BX+a[i-1]*AX*BY+b[i-1]*BX*AY+AY*BY)

构造矩阵：

| 1 0 0 0 0 |

| AX*BY AX 0 AX*BY 0 |

{AoD(n-1),a[i-1],b[i-1],a[i-1]*b[i-1],1}* | BX*AY 0 BX BX*AY 0 | ={AoD(n),a[i],b[i],a[i]*b[i],1}

| AX*BX 0 0 AX*BX 0 |

| AY*BY AY BY AY*BY 1 |

另外注意：

if(n==0){//这个判断条件很重要，没有就会超时
printf("0 ");
continue;
}

  1 #pragma comment(linker, "/STACK:1024000000,1024000000")
  2 #include<iostream>
  3 #include<cstdio>
  4 #include<cstring>
  5 #include<cmath>
  6 #include<math.h>
  7 #include<algorithm>
  8 #include<queue>
  9 #include<set>
 10 #include<bitset>
 11 #include<map>
 12 #include<vector>
 13 #include<stdlib.h>
 14 #include <stack>
 15 using namespace std;
 16 #define PI acos(-1.0)
 17 #define max(a,b) (a) > (b) ? (a) : (b)
 18 #define min(a,b) (a) < (b) ? (a) : (b)
 19 #define ll long long
 20 #define eps 1e-10
 21 #define MOD 1000000007
 22 #define N 1000000
 23 #define inf 1e12
 24 ll n;
 25 ll A0,Ax,Ay,B0,Bx,By;
 26 struct Matrix{
 27    ll mp[5][5];
 28 };
 29 Matrix Mul(Matrix a,Matrix b){
 30    Matrix res;
 31    for(ll i=0;i<5;i++){
 32       for(ll j=0;j<5;j++){
 33          res.mp[i][j]=0;
 34          for(ll k=0;k<5;k++){
 35             res.mp[i][j]=(res.mp[i][j]+(a.mp[i][k]*b.mp[k][j])%MOD+MOD)%MOD;
 36          }
 37       }
 38    }
 39    return res;
 40 }
 41 Matrix fastm(Matrix a,ll b){
 42    Matrix res;
 43    memset(res.mp,0,sizeof(res.mp));
 44    for(ll i=0;i<5;i++){
 45       res.mp[i][i]=1;
 46    }
 47    while(b){
 48       if(b&1){
 49          res=Mul(res,a);
 50       }
 51       a=Mul(a,a);
 52       b>>=1;
 53    }
 54    return res;
 55 }
 56 int main()
 57 {
 58    while(scanf("%I64d",&n)==1){
 59       scanf("%I64d%I64d%I64d%I64d%I64d%I64d",&A0,&Ax,&Ay,&B0,&Bx,&By);
 60       
 61       if(n==0){//这个判断条件很重要，没有就会超时
 62          printf("0
");
 63          continue;
 64       }
 65       
 66       
 67       ll a0=A0;
 68       ll b0=B0;
 69 
 70       Matrix tmp;
 71       memset(tmp.mp,0,sizeof(tmp.mp));
 72       tmp.mp[0][0]=1%MOD;
 73       tmp.mp[1][0]=Ax*By%MOD;
 74       tmp.mp[1][1]=Ax%MOD;
 75       tmp.mp[1][3]=Ax*By%MOD;
 76       tmp.mp[2][0]=Bx*Ay%MOD;
 77       tmp.mp[2][2]=Bx%MOD;
 78       tmp.mp[2][3]=Bx*Ay%MOD;
 79       tmp.mp[3][0]=Ax*Bx%MOD;
 80       tmp.mp[3][3]=Ax*Bx%MOD;
 81       tmp.mp[4][0]=Ay*By%MOD;
 82       tmp.mp[4][1]=Ay%MOD;
 83       tmp.mp[4][2]=By%MOD;
 84       tmp.mp[4][3]=Ay*By%MOD;
 85       tmp.mp[4][4]=1%MOD;
 86 
 87       Matrix cnt=fastm(tmp,n-1);
 88 
 89       Matrix g;
 90       memset(g.mp,0,sizeof(g.mp));
 91       g.mp[0][0]=a0*b0%MOD;
 92       g.mp[0][1]=a0%MOD;
 93       g.mp[0][2]=b0%MOD;
 94       g.mp[0][3]=a0*b0%MOD;
 95       g.mp[0][4]=1%MOD;
 96       Matrix ans=Mul(g,cnt);
 97       printf("%I64d
",ans.mp[0][0]);
 98    }
 99     return 0;
100 }

View Code