TopCoder SRM527 Div1Hard - Algorithmer’s note

問題

ある国のコインは、次のような単純な規則に従っている。

それぞれのコインは、正の整数の価値を持つ。
価値1のコインが存在する。
任意のコインの価値は、それより価値の低いコインの価値の整数倍である。

ある金額を払う時、これらのコインを使って何通りの払い方があるか。

解法

金額が大きすぎてDPができない。しかし基本精神はDP。
簡単のため、コインの金額が1, 3, 6の場合を考える。
すると、DPで解く場合、

A_i = B_i = C_i = 0 (i < 0)
A₀ = B₀ = C₀ = 1
A_i = A_i-1 (i >= 1)
B_i = A_i + B_i-3 (i >= 1)
C_i = B_i + C_i-6 (i >= 1)

という線形漸化式を考え、C_iが答えになることがわかる。
この漸化式は変形して、後の2式を

B_i = A_i-1 + B_i-3 (i >= 1)
C_i = A_i-1 + B_i-3 + C_i-6 (i >= 1)

としておく。（そのほうが後々都合が良い）
ここで、(A_i-1, B_i-1, C_i-1) -> (A_i, B_i, C_i) という線形変換があるととってもうれしい。しかし誰がどう考えてもそんな線形変換は存在しない。B_iの決定にはB_i-3が必要だからだ。
だが、iが「十分小さければ」B_i-3の項は0になるので無視できる。具体的には、i<=2ならばよい。
i=3の場合は、まずこの変換を3回行う（変換を複数回行うのは、行列の累乗なので簡単である）。すると、(A_i-3, B_i-3, C_i-3) -> (A_i, B_i, C_i)の変換「もどき」がわかる。「もどき」というのは、B_i-3が抜けているという意味である。ここにB_i-3の項を補完してやればよい。
すると、まだC_i-6の項を考えていないが、これも同様にすると補完できる。最終的に、(A_i-6, B_i-6, C_i-6) -> (A_i, B_i, C_i) の線形変換を表す行列までわかる。
結果、わかっている線形変換は、

(A_i-1, B_i-1, C_i-1) -> (A_i, B_i, C_i) (i<=2のとき用)
(A_i-3, B_i-3, C_i-3) -> (A_i, B_i, C_i) (i<=5のとき用)
(A_i-6, B_i-6, C_i-6) -> (A_i, B_i, C_i) (iはいつでもよい)

である。この3つを適切に組み合わせると、どんな場合でもC_iが求められる。

コード

#include <vector>
using namespace std;

typedef long long i64;
const int MOD = 1000003;

int N;

struct matrix
{
	i64 val[40][40];
};

matrix mul(const matrix& A, const matrix& B)
{
	matrix ret;
	for(int i=0;i<N;i++){
		for(int j=0;j<N;j++){
			ret.val[i][j] = 0;
			for(int k=0;k<N;k++){
				ret.val[i][j] += A.val[i][k] * B.val[k][j];
			}
			ret.val[i][j] %= MOD;
		}
	}
	return ret;
}

matrix pow(const matrix& A, i64 p)
{
	if(p==1) return A;
	matrix tmp = pow(A, p/2);
	tmp = mul(tmp, tmp);
	if(p%2==1) tmp = mul(tmp, A);
	return tmp;
}

class P8XCoinChange
{
	i64 ret[40];
	
	void linear(const matrix& A)
	{
		i64 r2[40];
		for(int i=0;i<N;i++){
			r2[i] = 0;
			for(int j=0;j<N;j++) r2[i] += A.val[i][j] * ret[j];
			r2[i] %= MOD;
		}
		for(int i=0;i<N;i++) ret[i] = r2[i];
	}
	
public:
	int solve(i64 coins, vector<i64> val)
	{
		N = val.size();
		
		matrix op;
		for(int i=0;i<N;i++){
			ret[i] = 1;
			for(int j=0;j<N;j++) op.val[i][j] = 0;
		}
		
		for(int i=0;i<N;i++){
			for(int j=i;j<N;j++) op.val[j][i]++;
			
			i64 times = (i==N-1)?(coins / val[i]):((coins % val[i+1]) / val[i]);
			if(times) linear(pow(op, times));
			
			if(i!=N-1) op = pow(op, val[i+1] / val[i]);
			
		}
		return (int)ret[N-1];
		
	}
};