HDU 6327 Problem I. Random Sequence（dp)

Description

给出一长度为$n$的序列$a_1,...,a_n$，如果$a_i=0$则$a_i$等概率取$[1,m]$中任意一个整数，给出权值序列$v_1,...,v_m$，求

$$\sum\limits_{i=1}^{n-3}v[gcd(a_i,a_{i+1},a_{i+2},a_{i+3})]$$
的期望

Input

第一行一整数$T$表示用例组数，每组用例首先输入两个整数$n,m$，之后输入$n$个整数$a_1,...,a_n$，最后输入$m$个整数$v_1,...,v_m$

$(1\le T\le 10,4\le n\le 100,1\le m\le 100,0\le a_i\le m,1\le v_i\le 10^9)$

Output

输出期望值，结果模$10^9+7$

Sample Input

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

Sample Output

8000
3

Solution

考虑朴素$dp$，以$dp[i][x][y][z]$表示$a_{i-2}=x,a_{i-1}=y,a_i=z$时的答案，每次枚举$a_{i+1}$的值转移即可，这样的时间复杂度为$O(nm^4)$，注意到实际有用的不是$a_{i-2},a_{i-1}$的值，而是$gcd(a_{i-2},a_{i-1},a_i)$和$gcd(a_{i-1},a_i)$的值，那么以$dp[i][x][y][z]$表示$gcd(a_{i-2},a_{i-1},a_i)=x,gcd(a_{i-1},a_i)=y,a_i=z$时的答案，那么显然需要$x|y|z\le m$，这样的三元组在$m=100$时只有$res=1471$个，进而可以直接转移

给所有三元组$(x,y,z)$编号$id(x,y,z)$，以$dp[i][j]$表示$gcd(a_{i-2},a_{i-1},a_i),gcd(a_{i-1},a_i),a_i$为第$j$个三元组时的答案，并预处理$suf(i,j)$表示第$i$个三元组加上$j$后转移到的三元组编号，$val(i,j)$表示第$i$个三元组加上$j$的$gcd$值，进而有转移

$dp[i][suf(j,k)]+=v_{val(j,k)}\cdot dp[i-1][j],k\in [l_i,r_i]$

其中$[l_i,r_i]$表示$a_i$的取值区间，答案即为$\frac{\sum dp[n][i]}{m^{num}}$，其中$num$为$a_i=0$的个数

Code

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
#define mod 1000000007
int mul(int x,int y)
{
    ll z=1ll*x*y;
    return z-z/mod*mod;
}
int add(int x,int y)
{
    x+=y;
    if(x>=mod)x-=mod;
    return x;
}
int Pow(int x,int y)
{
    int ans=1;
    while(y)
    {
        if(y&1)ans=(ll)ans*x%mod;
        x=(ll)x*x%mod;
        y>>=1;
    }
    return ans;
}
int T,n,m,a[105],v[105];
int res,dp[105][1500],id[105][105][105],gcd[105][105],suf[1500][105],val[1500][105];
void init(int n=100)
{
    for(int i=1;i<=n;i++)gcd[i][0]=gcd[0][i]=i;
    for(int i=1;i<=n;i++)
    {
        gcd[i][i]=i;
        for(int j=1;j<i;j++)
            gcd[i][j]=gcd[j][i]=gcd[j][i%j];
    }
    res=0;
    for(int x=1;x<=n;x++)
        for(int y=x;y<=n;y+=x)
            for(int z=y;z<=n;z+=y)
                id[x][y][z]=res++;
    for(int x=1;x<=n;x++)
        for(int y=x;y<=n;y+=x)
            for(int z=y;z<=n;z+=y)
                for(int w=1;w<=n;w++)
                {
                    suf[id[x][y][z]][w]=id[gcd[y][w]][gcd[z][w]][w];
                    val[id[x][y][z]][w]=gcd[x][w];
                }
}
int main()
{
    init();
    scanf("%d",&T);
    while(T--)
    {
        scanf("%d%d",&n,&m);
        for(int i=1;i<=n;i++)scanf("%d",&a[i]);
        for(int i=1;i<=m;i++)scanf("%d",&v[i]);
        memset(dp,0,sizeof(dp));
        int cnt=id[m][m][m]+1;
        for(int i=(a[1]?a[1]:1);i<=(a[1]?a[1]:m);i++)
            for(int j=(a[2]?a[2]:1);j<=(a[2]?a[2]:m);j++)
                for(int k=(a[3]?a[3]:1);k<=(a[3]?a[3]:m);k++)
                    dp[3][id[gcd[i][gcd[j][k]]][gcd[j][k]][k]]++;
        for(int i=4;i<=n;i++)
            for(int j=0;j<cnt;j++)
                if(dp[i-1][j])
                    for(int k=(a[i]?a[i]:1);k<=(a[i]?a[i]:m);k++)
                        dp[i][suf[j][k]]=add(dp[i][suf[j][k]],mul(v[val[j][k]],dp[i-1][j]));
        int ans=0;
        for(int i=0;i<cnt;i++)ans=add(ans,dp[n][i]);
        int num=0;
        for(int i=1;i<=n;i++)
            if(!a[i])num++;
        ans=mul(ans,Pow(Pow(m,mod-2),num));
        printf("%d\n",ans); 
    }
    return 0;
}