UVA 11762 Race to 1

Dilu have learned a new thing about integers, which is - any positive integer greater than 1 can be divided by at least one prime number less than or equal to that number. So, he is now playing with this property. He selects a number N. And he calls this D.
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;}