求和谐数的比例

uk702

mathe 发表于 2021-1-8 17:32
前面弄错了，30个数才到0.228800,40个数到0.229411,50个数到0.229762,60个数到0.230061

我用 julia 写了程序楞算，概率至少 > 0.232323，并且不是随着 n 增加而严格递增
n = 10000000, t = 0.2322919
n = 11000000, t = 0.2322950909090909
n = 12000000, t = 0.23229758333333334
n = 13000000, t = 0.23229561538461538
n = 14000000, t = 0.2322962142857143
n = 15000000, t = 0.23229826666666667
n = 16000000, t = 0.232299125
n = 17000000, t = 0.23229911764705882
n = 18000000, t = 0.23230083333333335
n = 19000000, t = 0.23230373684210526
n = 20000000, t = 0.23230505
n = 21000000, t = 0.23230633333333334
n = 22000000, t = 0.23230577272727274
n = 23000000, t = 0.2323062608695652
n = 24000000, t = 0.232308625
n = 25000000, t = 0.23230984
n = 26000000, t = 0.23230926923076922
n = 27000000, t = 0.23230948148148148
n = 28000000, t = 0.23231017857142858

1. using Primes
   
2. 
   
3. function factors(n)
   
4.     f = [one(n)]
   
5.     for (p,e) in factor(n)
   
6.         f = reduce(vcat, [f*p^j for j in 1:e], init=f)
   
7.     end
   
8.     return length(f) == 1 ? [one(n), n] : sort!(f)
   
9. end
   
10. 
    
11. t = 0
    
12. i = 1
    
13. while true
    
14.     a = factors(i)
    
15.     for u=1:length(a)
    
16.         for v=u+1:length(a)
    
17.             if a[u]+a[v] in a
    
18.                 # println((i, a[u], a[v]))
    
19.                 global t = t+1
    
20.                 @goto _next
    
21.             end   
    
22.         end   
    
23.     end
    
24. 
    
25. @label _next
    
26. global i=i+1
    
27. if i % 1000000 == 0 println("n = $i, t = $(t/i)") end
    
28. end
    

复制代码