公式f(x) = 1+f(x)*(1-g(x)/p(x))+ f(x/y)/p(x)
p(x)是不超过x的素数的个数
g(x)是p(x)个素数中是x因子的个数
然后移项(和我上次做的那一题差不多)
f(x)=(f(x/y)+p(x))/g(x)
y是x的因子 要枚举y
#include <cstdio> #include <cstring> const int maxn = 1000010; int vis[maxn]; int prime[maxn]; int prime_cnt; double f[maxn]; void sieve() { for(int i = 2; i <= 1000; i++) if(!vis[i]) for(int j = i*i; j <= 1000000; j += i) vis[j] = 1; } void get_primes() { sieve(); prime_cnt = 0; for(int i = 2; i <= 1000000; i++) if(!vis[i]) prime[prime_cnt++] = i; } double dp(int x) { if(x == 1) return 0.0; if(vis[x]) return f[x]; vis[x] = 1; double& ans = f[x]; ans = 0; int g = 0, p = 0; for(int i = 0; i < prime_cnt && prime[i] <= x; i++) { p++; if(x % prime[i] == 0) { g++; ans += dp(x / prime[i]); } } ans = (ans + p) / g; return ans; } int main() { get_primes(); memset(vis, 0, sizeof(vis)); int cas = 1; int T; scanf("%d", &T); while(T--) { int n; scanf("%d", &n); printf("Case %d: ", cas++); printf("%.10f\n", dp(n)); } return 0; }
原文:http://blog.csdn.net/u011686226/article/details/18939643