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
313
13
Sample Output
Case 1: 0.0000000000
Case 2: 2.0000000000
Case 3: 6.0000000000
题意：
给出一个整数N，每次在不超过N的素数中随机选择一个P，如果P是N的因数，则把N变成N/P，否则N不变。问平均情况下需要多少次选择才能把N变成1.
分析：
本题可以看成一个可以随机转移的状态机，其随机过程是马尔可夫过程，从每个状态出发的个转移状态的概率之和为1。
设f(i)表示当前数为i时接下来需要选择的次数期望，则根据期望的线性以及全期望公式可以为每个状态列出一个方程，设不超过x的素数有p(x)个，其中有g(x)个是x的因数，则：
f(x) = 1 + f(x) * (1 – g(x) / p(x)) + sum{f(x / y) * (1 / p(x)) | y是x的素因子}
边界为f(1) = 0，移项整理得f(x) = (sum{f(x / y) | y是x的素因子} + p(x)) / g(x)
注意到x / y < x，所以我们采用记忆化搜索的方式计算f(x)。
1 #include <cstdio>
2 #include <iostream>
3 #include <cstring>
4 #include <cmath>
5 using namespace std;
6 #define MAX 1000010
7 #define MAXP 700000
8 int vis[MAX],prime[MAXP],n;
9 double dp[MAX];
10 bool flag[MAX];
11 int primes_num;
12 void sieve(int n){
13 int m = (int)sqrt(n + 0.5);
14 memset(vis,0,sizeof(vis));
15 for(int i = 2 ; i <= m ; ++i)if(!vis[i])
16 for(int j = i * i ; j <= n ; j += i) vis[j] = 1;
17 }
18 int gen_primes(int n){
19 sieve(n);
20 int c = 0;
21 for(int i = 2 ; i <= n ; ++i)if(!vis[i]) prime[c++] = i;
22 return c;
23 }
24 double DP(int x){
25 if(x == 1) return 0;
26 if(flag[x]) return dp[x];
27 flag[x] = true;
28 double &ans = dp[x];
29 int g = 0,p = 0;
30 ans = 0;
31 for(int i = 0 ; i < primes_num && prime[i] <= x ; ++i){
32 p++;
33 if(x % prime[i] == 0) g++,ans += DP(x / prime[i]);
34 }
35 ans = (ans + p) / g;
36 return ans;
37 }
38 int main(){
39 primes_num = gen_primes(1000002);
40 memset(flag,false,sizeof(flag));
41 int T,num = 1;
42 scanf("%d",&T);
43 while(T--){
44 scanf("%d",&n);
45 printf("Case %d: %.10lf\n",num++,DP(n));
46 }
47 return 0;
48 }