hdu 4686 Arc of Dream

题目链接： 

　　http://acm.hdu.edu.cn/showproblem.php?pid=4686

题目描述：

　　题目上说的很清楚。

解题思路：

　　就是递推，用快速幂+构造矩阵解决，因为n的取值范围好大，好全。

/*
    |1 0      0   0   0 |
    |1 ax*bx  0   0   0 |
    |0 ax*by  ax  0   0 |*|s[n] f[n] a[n] b[n] 1|=|s[n+1] f[n+1] a[n+1] b[n+1] 1|
    |0 bx*ay  0   bx  0 |
    |0 by*ay  ay  by  1 |
    
    a[n+1]=a[n]*ax+ay;
    b[n+1]=b[n]*bx+by;
    f[n]=a[n]*b[n];
    s[n+1]=s[n]+f[n];
*/

#include <iostream>
#include <cstdlib>
#include <cstdio>
#include <algorithm>
#include <vector>
#include <queue>
#include <cstring>
using namespace std;

const int maxn = 5;
#define LL long long
#define mod 1000000007
struct mat
{
    LL p[maxn][maxn];
};

mat mul (mat a, mat b);
mat pow (LL n, mat a, mat b);
int main ()
{
    LL n, a, ax, ay, b, bx, by;
    mat A, B;
    while (scanf ("%lld", &n) != EOF)
    {
        memset (A.p, 0, sizeof(A.p));
        memset (B.p, 0, sizeof(B.p));
        scanf ("%lld %lld %lld %lld %lld %lld", &a, &ax, &ay, &b, &bx, &by);
        a %= mod, b %= mod, ax %= mod, bx %= mod, ay %= mod, by %= mod;
        B.p[0][0] = 0;//贴这个代码，我什么也不想说明，只想表明矩阵相乘取余要全面，要细心，（wa哭了）
        B.p[0][1] = a*b%mod;//心好累，在上面取完模还不够，在这里也要取模
        A.p[1][1] = ax*bx%mod;
        A.p[2][1] = ax*by%mod;
        A.p[3][1] = ay*bx%mod;
        A.p[4][1] = ay*by%mod;
        B.p[0][2] = a;
        B.p[0][3] = b;
        B.p[0][4] = 1;
        A.p[0][0] = A.p[1][0] = 1;
        
        A.p[2][2] = ax;
        A.p[3][3] = bx;
        A.p[4][4] = 1;
        A.p[4][2] = ay;
        A.p[4][3] = by;
        B = pow (n, A, B);
        printf ("%lld
", B.p[0][0]);
    }
    return 0;
}

mat mul (mat a, mat b)
{
    int i, j, k;
    mat c;
    memset (c.p, 0, sizeof(c.p));

    for (i=0; i<5; i++)
        for (j=0; j<5; j++)
        {
            for (k=0; k<5; k++)
                c.p[i][j] = (c.p[i][j] + a.p[i][k] * b.p[k][j]) % mod;
        }
    return c;
}
mat pow (LL n, mat a, mat b)
{
    while (n)
    {
        if (n % 2)
        {
            b = mul (b, a);
        }
        a = mul (a, a);
        n /= 2;
    }
    return b;
}