hdu 5564 Clarke and digits (数位dp + 矩阵快速幂优化)

@(K ACMer)

---

题意:
求长度在$[l, r]$之间,且相邻数和不为$k$的能被7整除的数的个数?

---

知识补充:
大数取模:对一个大数取k模,有这样的性质:

$$(k * 10 + a)\ \%\ mod\ = \ (k \ \%\ mod * 10+ a)\ \%\ mod$$ 所以对大数只需要从左到右的往复执行$乘10取模,加一位数再取模$操作即可.
数位DP:数位DP是指给给定一个数的长度的区间$[l, r]$问你在这个区间内满足某个条件的数的个数.通常就需要设置状态来进行递推求解.也通常是用$[0, r]$的满足条件的数的个数减去$[0,l- 1]$满足条件的数的个数.

---

分析:
对于本题,我们可以看出是典型的数位dp.设计状态:$dp[i][j][k]$为长度为i最后一位数位j且取余7为k的数的个数.这样有这样的状态转移方程:

$$dp[i][l][(k * 10 + l) \% 7]\qquad += \qquad dp[i][j][k]$$ 但是注意到这里的r的范围非常大为 $1e9$,显然需要 $O(logn)$的递推才行.
我们可以把每一个状态 $a[j][k](结尾数位j,取余7为k的数的个数)$用一个矩阵表示,然后从这个状态转移到下一个状态的操作就是乘以一个转移矩阵 $x$,即下一个状态为$a[j][k] * x$.所以我们需要跟据转移方程的特点构建这个转移矩阵.
其中利用了一个特性,若转移矩阵的 $x[i][j] = 1$就表示下一个状态的矩阵的第j列会加上原来矩阵的第i列.
我们把状态的二维通过压缩变成一维数组,然后就可以得到: $$x.v[compressor(j, k)][compressor(l, (k * 10 + l) \% 7)] = 1, (j + l != k);$$ 用这个来构建转移矩阵.
然后发现求最终状态dp[r],其实就是用初始状态(长度为1的数)dp[1]来乘上 $x^{r- 1}$,这里就转化成矩阵的幂次了,就可以用矩阵快速幂来加速为 $O(logn)$的递推了!

---

代码:

#include <iostream>
#include <cstdio>
#include <cstring>
#include <set>
#include <map>
#include <stack>
#include <vector>
#include <string>
#include <queue>
#include <cstdlib>
#include <cmath>
#include <algorithm>
using namespace std;
typedef pair<int, int> pii;
typedef unsigned long long ull;
typedef long long ll;
typedef vector<int> vi;
#define xx first
#define yy second
#define rep(i, a, n) for(int i = a; i < n; i++)
#define vep(c) for(decltype((c).begin()) it = (c).begin(); it != (c).end(); it++)
int mod = int(1e9) + 7, INF = 0x3fffffff, maxn = 1e5 - 1;

ll st[71], ed[71];

int compressor(int i, int j)
{
    return i * 7 + j;
}

class matrix
{
public:
    ll v[71][71];
    matrix(void) {
        memset(v, 0, sizeof(v));
    }
    matrix operator*(matrix& a) {
        matrix ans, b = *this;
        rep(i, 0, 71) {
            rep(j, 0, 71) {
                rep(k, 0, 71) {
                    ans.v[i][j] = (ans.v[i][j] + b.v[i][k] * a.v[k][j]) % mod;
                }
            }
        }
        return ans;
    }
    matrix operator^(int b) {
        matrix ans, a = *this;
        if (b < 0) return ans;
        rep(i, 0, 71) ans.v[i][i] = 1;
        while (b) {
            if (b & 1) ans = ans * a;
            a = a * a;
            b >>= 1;
        }
        return ans;
    }
};

int main(void)
{
    int T;
    scanf("%d", &T);
    rep(i, 1, 10) st[compressor(i, i % 7)]++;
    while (T--) {
        int l, r, m;
        scanf("%d%d%d", &l, &r, &m);
        matrix x;
        rep (l, 0, 10) {
            rep (j, 0, 10) {
                if (l + j == m) continue;
                rep (k, 0, 7) {
                    x.v[compressor(j, k)][compressor(l, (k * 10 + l) % 7)] = 1;
                }
            }
        }
        rep (i, 0, 10) x.v[compressor(i, 0)][70] = 1;
        x.v[70][70] = 1;

        matrix x1 = x ^ (r - 1);

        memset(ed, 0, sizeof(ed));
        rep (i, 0, 71) {
            rep (j, 0, 71) {
                ed[i] += st[j] * x1.v[j][i] % mod;
            }
        }
        ll ansr = ed[70];
        rep (i, 0, 10) ansr = (ansr + ed[compressor(i, 0)]) % mod;

        matrix x2 = x ^ (l - 2);
        memset(ed, 0, sizeof(ed));
        rep (i, 0, 71) {
            rep (j, 0, 71) {
                ed[i] += st[j] * x2.v[j][i] % mod;
            }
        }
        ll ansl = ed[70];
        rep (i, 0, 10) ansl = (ansl + ed[compressor(i, 0)]) % mod;
        printf("%lld\n", (ansr - ansl + mod) % mod);
    }
    return 0;
}