HDU 4686 Arc of Dream（矩阵快速幂）

Arc of Dream

Time Limit: 2000/2000 MS (Java/Others) Memory Limit: 65535/65535 K (Java/Others)

Total Submission(s): 2992 Accepted Submission(s): 921



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

1
1 2 3
4 5 6
2
1 2 3
4 5 6
3
1 2 3
4 5 6

Sample Output

4
134
1902

Author

Zejun Wu (watashi)

Source

2013 Multi-University Training Contest 9

题目大意：求

，已知a[0]=A0，a[i]=a[i-1]*AX+AY；b[0]=B0，b[i]=b[i-1]*BX+BY。

解题思路：

a[i]=a[i-1]*AX+AY；

b[i]=b[i-1]*BX+BY；

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

AoD(i)=AoD(i-1)+a[i]*b[i]；

则a
，b
，a
*b
，AoD(n)都存在递推公式，那么可以尝试构造矩阵，用来表示递推关系。

a[i] = AX*a[i-1]
+ 0*b[i-1] + 0*a[i-1]*b[i-1] + 0*AoD(i-1) + AY*1；

b[i] = 0*a[i-1] + BX*b[i-1]
+ 0*a[i-1]*b[i-1] + 0*AoD(i-1) + BY*1；

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

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

代码如下：

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <ctime>
#include <iostream>
#include <algorithm>
#include <string>
#include <vector>
#include <deque>
#include <list>
#include <set>
#include <map>
#include <stack>
#include <queue>
#include <numeric>
#include <iomanip>
#include <bitset>
#include <sstream>
#include <fstream>
#include <limits.h>
#define debug "output for debug\n"
#define pi (acos(-1.0))
#define eps (1e-6)
#define inf (1<<28)
#define sqr(x) (x) * (x)
#define mod 1000000007
using namespace std;
typedef long long ll;
typedef unsigned long long ULL;
#define MAX 15
ll n;
struct Matrix
{
ll a[MAX][MAX];
Matrix()
{
memset(a,0,sizeof(a));
}
};
Matrix operator *(Matrix a,Matrix b)
{
Matrix c;
memset(c.a,0,sizeof(c.a));
for(ll k=0;k<5;k++)
{
for(ll i=0;i<5;i++)
{
if(a.a[i][k]==0)
continue;
for(ll j=0;j<5;j++)
{
if(b.a[k][j]==0)
continue;
c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j])%mod;
c.a[i][j]%=mod;
}
}
}
return c;
}
Matrix operator ^(Matrix a,ll k)
{
Matrix c;
for(ll i=0;i<5;i++)
c.a[i][i]=1;
while(k)
{
if(k&1)
c=c*a;
a=a*a;
k>>=1;
}
return c;
}
int main()
{
ll i,j,k,t;
ll A0,AX,AY,B0,BX,BY;
Matrix a,b,c;
while(~scanf("%I64d%I64d%I64d%I64d%I64d%I64d%I64d",&n,&A0,&AX,&AY,&B0,&BX,&BY))
{
if(n==0)
{
printf("0\n");
continue;
}
a.a[0][0]=A0;a.a[0][1]=B0;a.a[0][2]=A0*B0%mod;a.a[0][3]=A0*B0%mod;a.a[0][4]=1;
b.a[0][0]=AX;b.a[0][1]=0; b.a[0][2]=AX*BY%mod;b.a[0][3]=AX*BY%mod;b.a[0][4]=0;
b.a[1][0]=0; b.a[1][1]=BX;b.a[1][2]=AY*BX%mod;b.a[1][3]=AY*BX%mod;b.a[1][4]=0;
b.a[2][0]=0; b.a[2][1]=0; b.a[2][2]=AX*BX%mod;b.a[2][3]=AX*BX%mod;b.a[2][4]=0;
b.a[3][0]=0; b.a[3][1]=0; b.a[3][2]=0;        b.a[3][3]=1;        b.a[3][4]=0;
b.a[4][0]=AY;b.a[4][1]=BY;b.a[4][2]=AY*BY%mod;b.a[4][3]=AY*BY%mod;b.a[4][4]=1;
c=b^(n-1);
ll ans=0;
for(i=0;i<5;i++)
ans=(ans+a.a[0][i]*c.a[i][3])%mod;
printf("%I64d\n",ans%mod);
}
return 0;
}