Fun with the Fellowship(矩阵快速幂)

Problem Description

Legolas, Gimli and Aragorn have just fought an army of orcs and now want to relax for the day. Their way of having fun is to play an archery contest with Legolas trying to hit any object that the other two point to. However, Aragorn soon realises that Legolas is too good for this game. He adds a twist to the game. Aragorn will now point to 4 targets and assign some number of points (a,b,c,d) to those 4 targets - the points will be unique for each target. Gimli will then give Legolas a number n and Legolas has to tell how many ways he can hit any of those 4 targets one after the other in any order so that the total points scored will be n. Now Legolas is a genius in archery but he's stumped by this question. Help Legolas!

You have to print the answer modulo 1000000007 (10^9+7).

Input Format

There is only one line in input with 5 space-separated integers - a b c d n

The last 3 test cases are for Extra Credit.

Constraints

1 <= a,b,c,d <= 10 and a,b,c,d are unique 1 <= n <= 10^5 For extra credit: 1 <= n <= 10^18

Output Format

Print the number of ways Legolas can get to the target number modulo 1000000007 (10^9+7). If it is impossible, print 0.

Sample Input 0

2 3 5 7 8

Sample Output 0

6

Explanation 0

The 6 ways of scoring a total of 8 are (2 + 2 + 2 + 2), (2 + 3 + 3), (3 + 2 + 3), (3 + 3 + 2), (3 + 5), and (5 + 3).

分析：

当n <= 10^5时，可以用一维dp解决，类似于多重背包，状态转移方程为：

dp[i] = dp[i-a] + dp[i-b] + dp[i-c] + dp[i-d]

当n达到10^18级别的时候，很明显可以根据一维dp的状态转移方程构造出系数矩阵，利用矩阵快速幂解决。

考虑到1<=a,b,c,d<=10，所以系数矩阵应该是10x10的，构造方法同一般的矩阵快速幂一样。

这个就是根据样例构造的系数矩阵，只需要在a,b,c,d对应的位置置为1，其余全为0即可。

其中，A就是上面构造的那个系数矩阵，接下来要做的就是求出f[1..10]和A^(n-10)。

下面给出代码：

#include<stdio.h>
#include<algorithm>
#define N 1000005
#define mod 1000000007
#define N 1000005
#include<string.h>
using namespace std;
typedef long long ll;

struct matrix
{
    ll mat[10][10];
}base,ans;

int dp[N];
ll n;
int a,b,c,d;
int f[11],w[4];

void init()
{
    sort(w, w+4);
    memset(f, 0, sizeof(f));
    a = w[0],b = w[1],c = w[2],d = w[3];
    f[a] = f[b] = f[c] = f[d] = 1;
    for(int i = 1; i <= 10; ++i)
        for(int j = 0; j < 4; ++j)
            if(i > w[j])
                f[i] += f[i-w[j]];
    memset(base.mat, 0, sizeof(base.mat));
    for(int i = 0; i < 4; ++i)
        base.mat[0][w[i]-1] = 1;
    for(int i = 1; i < 10; ++i)
        base.mat[i][i-1] = 1;
    memset(ans.mat, 0, sizeof(ans.mat));
    for(int i = 0; i < 10; ++i)
        ans.mat[i][i] = 1;
}

matrix matrixMultiply(matrix a, matrix b)
{
    matrix c;
    for(int i = 0; i < 10; ++i)
        for(int j = 0; j < 10; ++j)
        {
            c.mat[i][j] = 0;
            for(int k = 0; k < 10; ++k)
                c.mat[i][j] = (c.mat[i][j] + a.mat[i][k] * b.mat[k][j]) % mod;
        }
    return c;
}

void quickPow()
{
    ll m = n - 10;
    while(m)
    {
        if(m & 1)
            ans = matrixMultiply(ans, base);
        base = matrixMultiply(base, base);
        m >>= 1;
    }
}

int main()
{
    while(~scanf("%d%d%d%d%lld", &w[0],&w[1],&w[2],&w[3],&n))
    {
        init();
        if(n <= 10)
            printf("%d\n", f[n]);
        else
        {
            quickPow();
            ll ret = 0;
            for(int i = 0; i < 10; ++i)
                ret = (ret + ans.mat[0][i] * f[10-i]) % mod;
            printf("%lld\n", ret);
        }
    }
    return 0;
}

其中，init()函数中，初始化f[1..10]的过程就是求解小范围n的过程。只需要把10改成n即可。