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
a ₀ = A0
a _i = a _i-1*AX+AY
b ₀ = B0
b _i = b _i-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 10 ¹⁸, and all the other integers are no more than 2×10 ⁹.

 

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[n]，b[n]，a[n]*b[n]，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;
}