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 ",num++,DP(n)); 46 } 47 return 0; 48 }