Description
给出一长度为
n
n
的序列,如果
ai=0
a
i
=
0
则
ai
a
i
等概率取
[1,m]
[
1
,
m
]
中任意一个整数,给出权值序列
v1,...,vm
v
1
,
.
.
.
,
v
m
,求
的期望
Input
第一行一整数 T T 表示用例组数,每组用例首先输入两个整数,之后输入 n n 个整数,最后输入 m m 个整数
(1≤T≤10,4≤n≤100,1≤m≤100,0≤ai≤m,1≤vi≤109) ( 1 ≤ T ≤ 10 , 4 ≤ n ≤ 100 , 1 ≤ m ≤ 100 , 0 ≤ a i ≤ m , 1 ≤ v i ≤ 10 9 )
Output
输出期望值,结果模 109+7 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 d p ,以 dp[i][x][y][z] d p [ i ] [ x ] [ y ] [ z ] 表示 ai−2=x,ai−1=y,ai=z a i − 2 = x , a i − 1 = y , a i = z 时的答案,每次枚举 ai+1 a i + 1 的值转移即可,这样的时间复杂度为 O(nm4) O ( n m 4 ) ,注意到实际有用的不是 ai−2,ai−1 a i − 2 , a i − 1 的值,而是 gcd(ai−2,ai−1,ai) g c d ( a i − 2 , a i − 1 , a i ) 和 gcd(ai−1,ai) g c d ( a i − 1 , a i ) 的值,那么以 dp[i][x][y][z] d p [ i ] [ x ] [ y ] [ z ] 表示 gcd(ai−2,ai−1,ai)=x,gcd(ai−1,ai)=y,ai=z g c d ( a i − 2 , a i − 1 , a i ) = x , g c d ( a i − 1 , a i ) = y , a i = z 时的答案,那么显然需要 x|y|z≤m x | y | z ≤ m ,这样的三元组在 m=100 m = 100 时只有 res=1471 r e s = 1471 个,进而可以直接转移
给所有三元组 (x,y,z) ( x , y , z ) 编号 id(x,y,z) i d ( x , y , z ) ,以 dp[i][j] d p [ i ] [ j ] 表示 gcd(ai−2,ai−1,ai),gcd(ai−1,ai),ai g c d ( a i − 2 , a i − 1 , a i ) , g c d ( a i − 1 , a i ) , a i 为第 j j 个三元组时的答案,并预处理表示第 i i 个三元组加上后转移到的三元组编号, val(i,j) v a l ( i , j ) 表示第 i i 个三元组加上的 gcd g c d 值,进而有转移
dp[i][suf(j,k)]+=vval(j,k)⋅dp[i−1][j],k∈[li,ri] d p [ i ] [ s u f ( j , k ) ] + = v v a l ( j , k ) ⋅ d p [ i − 1 ] [ j ] , k ∈ [ l i , r i ]
其中 [li,ri] [ l i , r i ] 表示 ai a i 的取值区间,答案即为 ∑dp[n][i]mnum ∑ d p [ n ] [ i ] m n u m ,其中 num n u m 为 ai=0 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;
}