2018.09.28【HDU5564】Clarke and digits（数位DP）（矩阵快速幂）

传送门

---

解析：

哦，这可怕的数据范围。
但是仔细想一想，只限制了位数的话，能不能给我们什么特殊性质

思路：

首先要知道动态规划暴力 做法。
设$dp[i][j][k]$表示处理到第$i$位，前$i$位%7余数是$j$，当前位是$k$的方案数，于是状态转移就是$$dp[i+1][(j\ast 10+x)mod7][x]+=dp[i][j][k],1\le k,x\le 9$$

而最后的答案就是$\sum dp[r][0][0-9]-\sum dp[l-1][0][0-9]$
然而这个范围太大了，根本没有办法直接转移。

再观察一下，发现这个转移方程以$i$分为多个阶段，不考虑$i$的话，每个阶段的转移是一模一样的。

所以我们可以构造矩阵快速幂。

然而现在除掉$i$后状态是四维，怎么构造矩阵？

状态压缩啊。
我们考虑第三维状态只有10种情况
而第二维状态只有7种情况
也就是说：我们可以将第二维*10+第三维，来表示所有不同的70种情况。

我们还要统计前缀和。
我们可以在矩阵的最后一列10个位置加上1，最后一行最后一列的位置设为1，那么答案矩阵的这个位置就是我们要的前缀和。

---

代码：

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define re register
#define gc getchar
#define pc putchar
#define cs const

inline
int getint(){
	re int num;
	re int c;
	while(!isdigit(c=gc()));num=c^48;
	while(isdigit(c=gc()))num=(num<<1)+(num<<3)+(c^48);
	return num;
}

cs ll mod=1000000007;
cs int N=72;
struct matrix{
	ll a[N][N];
	matrix(int sign=0){
		memset(a,0,sizeof a);
		if(sign)for(int re i=0;i<N;++i)a[i][i]=sign;
	}
	ll *const operator[](int offset){
		return a[offset];
	}
	matrix friend operator*(cs matrix&A,cs matrix&B){
		matrix C(0);
		for(int re i=0;i<N;++i)
		for(int re j=0;j<N;++j)
		for(int re k=0;k<N;++k){
			C[i][j]=(C[i][j]+A.a[i][k]*B.a[k][j]%mod)%mod;
		}
		return C;
	}
}A,B,resl,resr;

matrix quickpow(matrix a,ll b){
	matrix ans(1);
	while(b){
		if(b&1)ans=ans*a;
		a=a*a;
		b>>=1;
	}
	return ans;
}

int sign(int i,int j){
	return i*10+j;
}

int T,K;
ll L,R;
signed main(){
	for(int re i=1;i<10;++i)A[0][sign(i%7,i)]=1;
	T=getint();
	while(T--){
		L=getint(),R=getint();
		K=getint();
		B=matrix(0);
		for(int re i=0;i<7;++i)
		for(int re j=0;j<10;++j)
		for(int re k=0;k<10;++k){
			if(j+k==K)continue;
			B[sign(i,j)][sign((i*10+k)%7,k)]=1;
		}
		for(int re i=0;i<10;++i)B[i][N-1]=1;
		B[N-1][N-1]=1;
		resr=A*quickpow(B,R);
		resl=A*quickpow(B,L-1);
		cout<<((resr[0][N-1]-resl[0][N-1])%mod+mod)%mod<<"\n";
	}
	return 0;
}