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 }

View Code

UVa11762