LightOJ 1038 Race to 1 Again（概率DP）

Race to 1 Again

Rimi learned a new thing about integers, which is - any positive integer greater than 1 can be divided by its divisors. 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 divisor of D (1 to D). Then he divides D by the number to obtain new D. He repeats this procedure until D becomes 1. What is the expected number of moves required for N to become 1.

Input
Input starts with an integer T (≤ 10000), denoting the number of test cases.
Each case begins with an integer N (1 ≤ N ≤ 105).

Output
For each case of input you have to print the case number and the expected value. Errors less than 10-6 will be ignored.

Sample Input

3
1
2
50

Sample Output

Case 1: 0
Case 2: 2.00
Case 3: 3.0333333333

题意
对于一个数$n$，求把它变回1，所需要的除的次数的期望

题解
概率DP。
令$dp[i]$表示需要次数为$i$的期望
状态转移方程：
$$\begin{array}{rl}dp[i]& =\frac{{\sum }_{d|i}dp[d]}{num[i]}+1\\ dp[i]num[i]& =\sum _{d|i}dp[d]+num[i]\\ dp[i](num[i]-1)& =\sum _{d|i,d≠i}dp[d]+num[i]\\ dp[i]& =\frac{{\sum }_{d|i,d≠i}dp[d]+num[i]}{num[i]-1}\end{array}$$
其中$d|i$表示$d$是$i$的因子，$num[i]$表示$i$的因子个数

求解时，对每个$i(i>1)$都通过公式求出其$dp[i]$,并对$num[i]+1$（$i$本身也是因子），然后对$i$的倍数$j$都进行$dp[j]+=dp[i]$并且$num[j]+1$,
代码

#include <iostream>
#include <algorithm>
#include <cstdio>
#include <queue>
#include <cmath>
#include <string>
#include <cstring>
#include <map>
#include <set>
#include <cmath>

using namespace std;
#define me(x,y) memset(x,y,sizeof x)
#define MIN(x,y) x < y ? x : y
#define MAX(x,y) x > y ? x : y
typedef long long ll;
typedef unsigned long long ull;
const int maxn = 1e5;
const ll INF = 0x3f3f3f3f;
const int MOD = 998244353;
const double eps = 1e-06;

double dp[maxn+10];
int num[maxn+10];
void get_dp(){
    me(dp,0);
    me(num,0);
    dp[1] = 0;
    for(int i = 1; i <= maxn; ++i){
        num[i]++;
        if(i>1)dp[i] = (dp[i]+num[i])/(num[i]-1);
        for(int j = i*2; j <= maxn; j += i){
            dp[j] += dp[i];
            num[j]++;
        }
    }
}

int main(){
    int t;
    get_dp();
    cin>>t;
    for(int ca = 1; ca <= t; ca++){
        int n;
        cin>>n;
        printf("Case %d: %lf\n",ca,dp[n]);
    }
    return 0;
}


/*


*/