UVa11762

Posted 2020-08-04 Yan_Bin

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