题目链接：HDU4686
Arc of Dream
Time Limit: 2000/2000 MS (Java/Others) Memory Limit: 65535/65535 K (Java/Others)
Total Submission(s): 5059 Accepted Submission(s): 1584
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
分析：ai*bi=(ai-1 *ax+ay)(bi-1 *bx+by)
=(ai-1 * bi-1 ax*bx)+(ai-1 *ax*by)+(bi-1 *bx*ay)+(ay*by)
设p=ax*bx, q=ax*by, r=ay*bx, s=ay*by
所以ai*bi=p(ai-1 bi-1)+q(ai-1)+r(bi-1)+s*
虽然可以用递推来求出每一项，但是n太大了，直接求绝对会超时的。
设f(n)=an*bn, a(n)=an, b(n)=bn
s(n)=sum(ai*bi),i=0,1,…n
则f(i)=p*f(i-1)+q*a(i-1)+r*b(i-1)+s
这是一个递推式，对于任何一个递推式，我们都可以用矩阵法来优化，加快速度求出第n项或前n项和。
我们可以构造一个5*5的矩阵A，使得
*【f(n-1),a(n-1),b(n-1),1,s(n-2)】*A=【f(n),a(n),b(n),1,s(n-1)】
=【p*f(n-1)+q*a(n-1)+r*b(n-1)+s, a(n-1)ax+ay, b(n-1)*bx+by, 1, s(n-2)+f(n-1)】
所以我们容易得出矩阵A：
axbx 0 0 0 1
axby ax 0 0 0
aybx 0 bx 0 0
ayay ay by 1 0
0 0 0 0 1
由【f(1), a(1) ,b(1), 1, s(0)】*A = 【f(2), a(2), b(2), 1, s(1)】
以此类推得，【f(1), a(1) ,b(1), 1, s(0)】*A^(n-1) = 【f(n), a(n), b(n), 1, s(n-1)】
这样就可以快速的求出s(n-1)了，
其中f(1)=a1*b1, a(1)=a0*ax+ay,
b(1)=b0*bx+by, s(0)=a0*b0
接下来就是矩阵快速幂了。
注意：n==0时，直接输出0，不然会死循环TLE的，还有就是要用long long,也要记得mod
#include<iostream>
#include<algorithm>
#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<queue>
#include<map>
#include<set>
#include<stack>
#include<bitset>
#include<numeric>
#include<vector>
#include<string>
#include<iterator>
#include<cstring>
#include<ctime>
#include<functional>
#define INF 0x3f3f3f3f
#define ms(a,b) memset(a,b,sizeof(a))
#define pi 3.14159265358979
#define mod 1000000007
#define lson l,m,rt<<1
#define rson m+1,r,rt<<1|1
using namespace std;
typedef pair<int, int> P;
typedef long long ll;
typedef unsigned long long ull;
const int maxn = 50;
struct Matrix {
ll a[6][6];
}ori, A, temp, res, ans;
int n;
Matrix mul(Matrix x, Matrix y)
{
ms(temp.a, 0);
for (int i = 1; i <= n; i++)
{
for (int j = 1; j <= n; j++)
{
ll sum = 0;
for (int k = 1; k <= n; k++)
{
sum = (sum + x.a[i][k] * y.a[k][j] % mod) % mod;
}
temp.a[i][j] = sum;
}
}
return temp;
}
void quickpow(ll k)
{
ms(res.a, 0);
for (int i = 1; i <= n; i++) res.a[i][i] = 1;
while (k)
{
if (k & 1) res = mul(res, A);
A = mul(A, A);
k >>= 1;
}
}
int main()
{
ll N, a0, ax, ay, b0, bx, by;
while (~scanf("%lld %lld %lld %lld %lld %lld %lld", &N, &a0, &ax, &ay, &b0, &bx, &by))
{
if (!N)
{
puts("0");
continue;
}
ll a1 = (a0*ax%mod + ay) % mod, b1 = (b0*bx%mod + by) % mod, f1 = (a1*b1) % mod, s0 = (a0*b0) % mod;
n = 5;
ms(ori.a, 0);
ori.a[1][1] = f1;
ori.a[1][2] = a1;
ori.a[1][3] = b1;
ori.a[1][4] = 1;
ori.a[1][5] = s0;
memset(A.a, 0, sizeof(A.a));
A.a[1][1] = (ax*bx) % mod;
A.a[1][5] = 1;
A.a[2][1] = (ax*by) % mod;
A.a[2][2] = ax%mod;
A.a[3][1] = (ay*bx) % mod;
A.a[3][3] = bx%mod;
A.a[4][1] = (ay*by) % mod;
A.a[4][2] = ay%mod;
A.a[4][3] = by%mod;
A.a[4][4] = 1;
A.a[5][5] = 1;
quickpow(N - 1);
ans = mul(ori, res);
printf("%lld
", ans.a[1][5]);
}
return 0;
}