HDU 6327 Random Sequence（动态规划）

题目链接：Random Sequence

题意

有一个长度为 $n$ 的序列 $a_1,a_2,\cdots,a_n~(0\leq a_i\leq m)$ 和长度为 $m$ 的序列 $v_1,v_2,\cdots,v_m$，如果 $a_i=0$，则 $a_i$ 的值可能为 $1$ 到 $m$ 之间的任意值，取得任意值的概率为 $\frac{1}{m}$，求公式：

$$\prod_{i=1}^{n-3}v_{\gcd(a_i,a_{i+1},a_{i+2},a_{i+3})}$$ 的期望值。

输入

第一行为一个整数 $T~(1\leq T\leq10)$，接下去有 $T$ 组数据，每组数据第一行为两个整数 $n,m~(4\leq n\leq100,1\leq m\leq100)$，第二行为 $n$ 个整数 $a_1,a_2,\cdots,a_n~(0\leq a_i\leq m)$，第三行为 $m$ 个整数 $v_1,v_2,\cdots,v_m~(1\leq v_i\leq10^9)$。

输出

输出所求公式的期望值对 $10^9+7$ 取模后的结果，若期望值无法用整数表示，则先将期望值化简为 $\frac{A}{B}$ 后，输出 $C~(0\leq C<10^9+7)$ 的值，其中 $C$ 满足 $A=B\times C(\mod 10^9+7)$。

样例

输入
2 6 8 4 8 8 4 6 5 10 20 30 40 50 60 70 80 4 3 0 0 0 0 3 2 4
输出
8000 3

题解

由于任意两个 $v$ 之间的选择是相互独立的，所以 $E(v_1v_2)=E(v_1)E(v_2)$。
如果定义状态 $dp[i][x][y][z]$ 表示到第 $i$ 位 $a_i=x,a_{i-1}=y,a_{i-2}=z$ 时 $\prod_{j=4}^iv_{\gcd(a_j,a_{j-1},a_{j-2},a_{j-3})}$ 的所有合法情况的和，则可以通过枚举第 $a_{i+1}$ 位的所有合法值 $xx$ 计算 $dp[i+1][xx][x][y]$ 的值，其递推公式为：

$$dp[i+1][xx][x][y]=\sum_{x,y,z\in[1,m]}dp[i][x][y][z]\times v_{\gcd(xx,x,y,z)}$$ 注意当 $a_i$ 的值非 $0$ 时 $x$ 的值只能取 $a_i$，只有 $a_i$ 的值为 $0$ 时 $x$ 才可以从 $1$ 到 $m$ 枚举所有值，对 $xx,y,z$ 同理。这种递推式的时间复杂度为 $O(nm^4)$，必然超时。
我们可以定义状态 $dp[i][x][y][z]$ 表示到第 $i$ 位 $a_i=x,\gcd(x,a_{i-1})=y,\gcd(x,a_{i-1},a_{i-2})=z$ 时所有合法情况下 $\prod_{j=4}^iv_{\gcd(a_j,a_{j-1},a_{j-2},a_{j-3})}$ 的和，由公约数关系可以得到， $x,y,z$ 必须满足 $z$ 能整除 $y$， $y$ 整除 $x$，于是有下列递推式：
$$dp[i+1][xx][\gcd(xx,x)][\gcd(xx,y)]=\sum_{x,y,z\in[1,m],z|y,y|x}dp[i][x][y][z]\times v_{\gcd(xx,z)}$$ 其中 $xx,x,y,z$ 在 $a_i$ 非 $0$ 时的取值限制同上，$a|b$ 表示 $a$ 能整除 $b$，在 $z|y,y|x$ 的限制下， $x,y,z$ 所有合法的状态可以减少到 $1471$ 种，因此时间复杂度为 $O(1471\times nm)$，其中 $\gcd$ 的计算可以 $O(m^2\log m)$ 通过预处理将复杂度降低到 $O(1)$。
最后是 $dp$ 的初始值的确定：
$$dp[3][x][\gcd(x,y)][\gcd(x,y,z)]=\sum_{x,y,z\in[1,m]}1$$ 其中 $x,y,z$ 的取值限制同上， $dp$ 初始值即表示取前 $3$ 位为 $x,y,z$ 的合法状态数，系数为 $1$，初始化方式为暴力，时间复杂度为 $O(m^3)$。
最后将 $\sum_{x,y,z\in[1,m],z|y,y|x}dp[n][x][y][z]$ 乘上总方案数对 $10^9+7$ 的逆元就是答案，总方案数为 $m^{cnt}$，其中 $cnt$ 为 $a$ 序列中 $0$ 的个数。

过题代码

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <climits>
#include <cstring>
#include <string>
#include <vector>
#include <list>
#include <queue>
#include <stack>
#include <map>
#include <set>
#include <bitset>
#include <algorithm>
#include <functional>
#include <iomanip>
using namespace std;

#define LL long long
const int maxn = 101;
const int MOD = 1000000007;
int T, n, m, ans, cnt, inv;
int a[maxn], v[maxn];
int Gcd[maxn][maxn], dp[maxn][maxn][maxn][maxn];

int add(int a, int b) {
    a += b;
    if(a >= MOD) {
        return a - MOD;
    }
    if(a < 0) {
        return a + MOD;
    }
    return a;
}

int gcd(int x, int y) {
    return y == 0? x: gcd(y, x % y);
}

void prepare_gcd() {
    for(int i = 1; i < maxn; ++i) {
        Gcd[0][i] = Gcd[i][0] = i;
        for(int j = 1; j < maxn; ++j) {
            Gcd[i][j] = gcd(i, j);
        }
    }
}

void Init() {
    for(int i = 3; i <= n; ++i) {
        for(int z = 1; z <= m; ++z) {
            for(int y = z; y <= m; y += z) {
                for(int x = y; x <= m; x += y) {
                    dp[i][x][y][z] = 0;
                }
            }
        }
    }
    for(int x = 1; x <= m; ++x) {
        if(a[3] != 0 && x != a[3]) {
            continue;
        }
        for(int y = 1; y <= m; ++y) {
            if(a[2] != 0 && a[2] != y) {
                continue;
            }
            for(int z = 1; z <= m; ++z) {
                if(a[1] != 0 && z != a[1]) {
                    continue;
                }
                ++dp[3][x][Gcd[x][y]][Gcd[Gcd[x][y]][z]];
            }
        }
    }
}

LL fast_pow(LL res, LL n) {
    LL ans;
    for(ans = 1; n != 0; n >>= 1) {
        if((n & 1) == 1) {
            ans = (ans * res) % MOD;
        }
        res = (res * res) % MOD;
    }
    return ans;
}

int main() {
    #ifdef Dmaxiya
    freopen("test.txt", "r", stdin);
//    freopen("test1.out", "w", stdout);
    #endif // Dmaxiya
    ios::sync_with_stdio(false);

    prepare_gcd();
    scanf("%d", &T);
    while(T--) {
        cnt = 0;
        scanf("%d%d", &n, &m);
        for(int i = 1; i <= n; ++i) {
            scanf("%d", &a[i]);
            if(a[i] == 0) {
                ++cnt;
            }
        }
        for(int i = 1; i <= m; ++i) {
            scanf("%d", &v[i]);
        }
        Init();
        for(int i = 3; i < n; ++i) {
            for(int z = 1; z <= m; ++z) {
                for(int y = z; y <= m; y += z) {
                    for(int x = y; x <= m; x += y) {
                        if(dp[i][x][y][z] == 0) {
                            continue;
                        }
                        for(int xx = 1; xx <= m; ++xx) {
                            if(a[i + 1] != 0 && a[i + 1] != xx) {
                                continue;
                            }
                            dp[i + 1][xx][Gcd[xx][x]][Gcd[xx][y]] = add(dp[i + 1][xx][Gcd[xx][x]][Gcd[xx][y]], (LL)dp[i][x][y][z] * v[Gcd[xx][z]] % MOD);
                        }
                    }
                }
            }
        }
        ans = 0;
        for(int k = 1; k <= m; ++k) {
            for(int j = k; j <= m; j += k) {
                for(int i = j; i <= m; i += j) {
                    ans = add(ans, dp[n][i][j][k]);
                }
            }
        }
        ans = (LL)ans * fast_pow(fast_pow(m, cnt), MOD - 2) % MOD;
        printf("%d\n", ans);
    }

    return 0;
}