题意:给一个整数N,每次可以在不超过N的素数中等概率随机选择一个P,如果P是N的约数,N变为N/P,否则N不变。问N变到1的期望次数是多少。
思路:首先筛法打表,把范围内的素数都记录下来。我们可以用动态规划记忆化搜索的方法求解(递推会超时)。显然,N=1的时候,结果为0(边界)。根据期望的定义,可以列出方程。其中p(j)是选择概率,i取遍所有不大于j的素数。
E(j)+=p(j)*(E(i)+1) 当j可以被i整除时
E(j)+=p(j)*(E(j)+1) 当j不能被i整除时
最后求解这个方程算E(j)就是了。
这题告诉我们有时候记忆化搜索比递推是要快的,因为只求解了需要求解的状态,而不是全部的状态。另外,白书上说,这个随机过程叫马尔可夫过程。
#include <iostream>
#include <stdio.h>
#include <cmath>
#include <algorithm>
#include <iomanip>
#include <cstdlib>
#include <string>
#include <string.h>
#include <vector>
#include <queue>
#include <stack>
#include <map>
#include <assert.h>
#include <set>
#include <ctype.h>
#define ll long long
#define max3(a,b,c) max(a,max(b,c))
using namespace std;
#define maxn 10000001
bool vis[maxn];
int prime[maxn];
int cnt[maxn];
double dp[maxn];
double fun(int n){
if(!n)return 0.0;
if(dp[n]>0.000000001)return dp[n];
double p=1.0/cnt[n];
int k=0;
double a=0.0;
for(int j=0;j<cnt[n];j++){
if((n%prime[j])==0){
a+=(fun(n/prime[j])+1)*p;
k++;
}else{
}
}
if(!k){
dp[n]=a;
return a;
}
double t=(k+0.0)/cnt[n];
dp[n]=(a+1.0-t)/t;
return dp[n];
}
int main(){
memset(vis,0,sizeof(vis));
memset(cnt,0,sizeof(cnt));
vis[0]=vis[1]=1;
for(int i=2;i<sqrt(maxn);i++){
for(int j=i*i;j<maxn;j+=i){
vis[j]=1;
}
}
int k=0;
for(int i=1;i<maxn;i++){
if(!vis[i]){
prime[k++]=i;
cnt[i]=k;
}else{
cnt[i]=cnt[i-1];
}
}
int t;
cin>>t;
int cas=0;
while(t--){
cas++;
int n;
cin>>n;
double ans=fun(n);
printf("Case %d: %.10lf\n",cas,ans);
}
return 0;
}