In each turn he randomly chooses a prime number less than or equal to D. If D is divisible by the prime number then he divides D by the prime number to obtain new D. Otherwise he keeps the old D. He repeats this procedure until D becomes 1. What is the expected number of moves required for N to become 1.
[We say that an integer is said to be prime if its divisible by exactly two different integers. So, 1 is not a prime, by definition. List of first few primes are 2, 3, 5, 7, 11, …]
Input
Input will start with an integer T (T <= 1000), which indicates the number of test cases. Each of the next T lines will contain one integer N (1 <= N <= 1000000).
Output
For each test case output a single line giving the case number followed by the expected number of turn required. Errors up to 1e-6 will be accepted.
Sample Input
3
1
3
13
Output for Sample Input
Case 1: 0.0000000000
Case 2: 2.0000000000
Case 3: 6.0000000000
解析
这是我写出的第一道概率DP,用的记忆化搜索
先推公式,f(i)表示从i变化到1的期望步数,按题意有:
f(1)=0
f(2)=1 + f(1)
f(3)=1 + 1/2 f(3) + 1/2 f(1) ==> 1/2 f(3)=1 + 1/2 f(1) ==> f(3)=2 + f(1)
f(4)=1 + 1/2 f(2) + 1/2 f(4) ==> 1/2 f(4)=1 + 1/2 f(2) ==> f(4)=2 + f(2)
f(5)=1 + 2/3 f(5) + 1/3 f(1) ==> 1/3 f(5)=1 + 1/3 f(1) ==> f(5)=3 + f(1)
f(6)=1 + 1/3 f(3) + 1/3 f(2) + 1/3 f(6) ==> 2/3 f(6)=1 + 1/3 f(3) + 1/3 f(2) ==> f(6)=[3+f(3)+f(2) ] / 2
我用bfprm[i]表示i前面有多少个包括i自己的素数,这里我有一个小小的优化,通过刚才的推导我发现对于素数a有 f(a)=bfprm[a]
对于其他的合数有f(x)=1 + [1 - x的因子数/bfprm[x] ] * f(x) + f(x/p1) / bfprm[x] + f(x/p2) / bfprm[x]+...+ f(x/pi) / bfprm[x]
其中p1 p2 ...pi 都表示x的因数
参阅《算法竞赛入门经典训练指南(第一版)》P143
#include<cstdio>
#include<cstring>
#include<cmath>
using namespace std;
double f[1000100];
int bfprm[1000100],prm[79000];
bool vis[1000100];
bool is_prime(int a)
{
int e=(int)sqrt(a);
for(int i=2;i<=e;i++)
if(a%i==0) return 0;
return 1;
}
double dp(int a)
{
if(vis[a]) return f[a];
int factor=0;
for(int j=1;j<=bfprm[a];j++)
if(a%prm[j] == 0)
{
factor++;
f[a]+=dp(a/prm[j]);
}
f[a]+=(double)bfprm[a];
f[a]/=factor;
vis[a]=1;
return f[a];
}
int main()
{
freopen("B.in","r",stdin);
freopen("B.out","w",stdout);
memset(f,0,sizeof(f));
memset(vis,0,sizeof(vis));
f[1]=0.0000000000;vis[1]=1;
int cnt=0;
for(int i=2;i<1000010;i++)
{
if(is_prime(i))
{
prm[++cnt]=i;
f[i]=(double)cnt;
vis[i]=1;
}
bfprm[i]=cnt;
}
int T;
scanf("%d",&T);
for(int i=1;i<=T;i++)
{
printf("Case %d: ",i);
int n; scanf("%d",&n);
printf("%.10lf",dp(n));
printf("\n");
}
return 0;
}