HDU 4686（矩阵快速幂）

题目链接：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;
}