数学研发论坛

 找回密码
 欢迎注册
12
返回列表 发新帖
楼主: lsr314

[提问] 求和谐数的比例

[复制链接]
发表于 2021-1-9 11:58:20 | 显示全部楼层
本帖最后由 uk702 于 2021-1-9 12:06 编辑
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. function factors(n)
  3.     f = [one(n)]
  4.     for (p,e) in factor(n)
  5.         f = reduce(vcat, [f*p^j for j in 1:e], init=f)
  6.     end
  7.     return length(f) == 1 ? [one(n), n] : sort!(f)
  8. end

  9. t = 0
  10. i = 1
  11. while true
  12.     a = factors(i)
  13.     for u=1:length(a)
  14.         for v=u+1:length(a)
  15.             if a[u]+a[v] in a
  16.                 # println((i, a[u], a[v]))
  17.                 global t = t+1
  18.                 @goto _next
  19.             end   
  20.         end   
  21.     end

  22. @label _next
  23. global i=i+1
  24. if i % 1000000 == 0 println("n = $i, t = $(t/i)") end
  25. end
复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2021-2-27 15:41 , Processed in 0.048750 second(s), 14 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表