首页
登录 | 注册

HDU 6327 2018HDU多校赛 第三场 Random Sequence(概率期望dp+数论)

Problem I. Random Sequence

Time Limit: 2000/1000 MS (Java/Others)    Memory Limit: 524288/524288 K (Java/Others)
Total Submission(s): 200    Accepted Submission(s): 69


 

Problem Description

There is a positive integer sequence a1,a2,...,an with some unknown positions, denoted by 0. Little Q will replace each 0 by a random integer within the range [1,m] equiprobably. After that, he will calculate the value of this sequence using the following formula :

∏i=1n−3v[gcd(ai,ai+1,ai+2,ai+3)]


Little Q is wondering what is the expected value of this sequence. Please write a program to calculate the expected value.

 

 

Input

The first line of the input contains an integer T(1≤T≤10) , denoting the number of test cases.
In each test case, there are 2 integers n,m(4≤n≤100,1≤m≤100) in the first line, denoting the length of the sequence and the bound of each number.
In the second line, there are n integers a1,a2,...,an(0≤ai≤m) , denoting the sequence.
In the third line, there are m integers v1,v2,...vm(1≤vi≤109) , denoting the array v .

 

 

Output

For each test case, print a single line containing an integer, denoting the expected value. If the answer is AB , please print C(0≤C<109+7) where AC×B(mod109+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

 

 

Source

2018 Multi-University Training Contest 3

 

 

大致题意:给你一个数列,里面的数字要么是0,要么是1~m中的一个数字。你需要随机的把这些0替换成1~m中的任意一个,然后再计算着整个数列的权值。权值定义为HDU 6327 2018HDU多校赛 第三场 Random Sequence(概率期望dp+数论)
,现在问你最后权值的期望。

由于这个v数组是输入的,也就是说没有什么性质,所以这个权值的表达式也不能用数论的知识推出什么结论。因此我们考虑暴力的dp。我们令dp[i][j][k][l]表示当前计算到第i位,且a[i]取j,a[i-1]取k,a[i-2]取l时的期望的分子。那么,从这个状态出发,枚举a[i+1]的数值,就看求出着四个数字的gcd,对应转移过去即可。如此一来,时间复杂度就是HDU 6327 2018HDU多校赛 第三场 Random Sequence(概率期望dp+数论)
。显然是不可以接受的。

接下来,我们考虑进行优化。我们注意到,如果a[i]取2,a[i-1]取4和6都是一样的,也就是说,每一维我们根本没有必要从1..m的枚举,我们考虑枚举对应的gcd即可。由此,我们重新定义dp数组,dp[i][j][k][l]表示当前计算到第i位,且a[i]取j,gcd(a[i],a[i-1])取k,gcd(k,a[i-2])取l时的期望的分子。但是这么表示,枚举起来还是很复杂,于是我们考虑提前预处理所有合法的状态,然后把每一个状态对应一个数字,如此数组就变成一个二维的dp[i][j]表示当前处理到第i个数字,然后最近三个数字的状态是状态j的期望的分子。

然后转移目标就是下一位,然后对应修改三个gcd的值,转移的量就是四个数字的gcd与dp[i][j]的乘积。这样总的复杂度就是O(N*M*tot+M^3) 其中tot表示合法的状态数目,M^3位初始化的复杂度。最后的结论是,合法状态数不超过1471个。具体实现见代码:

#include<bits/stdc++.h>
#define INF 0x3f3f3f3f
#define IO ios::sync_with_stdio(0);cin.tie(0);cout.tie(0)
#define mod 1000000007
#define N 110
#define LL long long


using namespace std;

int dp[N][1500],gcd[N][N],id[N][N][N];
struct node{int x,y,z;} st[1500];
int n,m,tot,v[N],a[N],inv[N];

void init()
{
    tot=0; inv[1]=1;
    for(int i=2;i<N;i++)
        inv[i]=1LL*(mod-mod/i)*inv[mod%i]%mod;
    for(int i=1;i<=100;i++)
        for(int j=1;j<=i;j++)
            gcd[i][j]=gcd[j][i]=__gcd(i,j);
    for(int i=1;i<=100;i++)
        for(int j=1;j<=i;j++)
            if (i%j==0)
            for(int k=1;k<=j;k++)
            {
                if (j%k) continue;
                st[++tot]=node{i,j,k};
                id[i][j][k]=tot;
            }
}

int main()
{
    IO; init();
    int T; cin>>T;
    while(T--)
    {
        cin>>n>>m;
        memset(dp,0,sizeof(dp));
        for(int i=1;i<=n;i++) cin>>a[i];
        for(int i=1;i<=m;i++) cin>>v[i];
        for(int i=1;i<=m;i++)
            if (!a[3]||a[3]==i)
            for(int j=1;j<=m;j++)
                if (!a[2]||a[2]==j)
                for(int k=1;k<=m;k++)
                {
                    if (a[1]&&a[1]!=k) continue;
                    int y=gcd[i][j],z=gcd[y][k];
                    int nxt=id[i][y][z];
                    dp[3][nxt]++;
                }
        for(int j=3;j<n;j++)
        {
            for(int i=1;i<=tot;i++)
            {
                if (!dp[j][i]) continue;
                int x=st[i].x,y=st[i].y,z=st[i].z;
                for(int k=1;k<=m;k++)
                {
                    if (a[j+1]&&k!=a[j+1]) continue;
                    int yy=gcd[k][x],zz=gcd[k][y];
                    int nxt=id[k][yy][zz];
                    dp[j+1][nxt]+=1LL*v[gcd[k][z]]*dp[j][i]%mod;
                    if (dp[j+1][nxt]>=mod) dp[j+1][nxt]-=mod;
                }
            }
        }
        LL ans=0;
        for(int i=1;i<=tot;i++)
        {
            ans+=dp[n][i];
            if (ans>=mod) ans-=mod;
        }
        for(int i=1;i<=n;i++)
            if (!a[i]) ans=ans*inv[m]%mod;
        cout<<ans<<endl;
    }
    return 0;
}

 



2020 jeepxie.net webmaster#jeepxie.net
10 q. 0.009 s.
京ICP备10005923号