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;
}