調和数

数学物語テトラちゃんとハーモニック・ナンバーに出てきた、調和数のグラフ用のデータを作ります。
まずは、Floatを使って、ナイーブに。

def harmonic(n)
  s = 0.0
  (1..n).each do |k|
    s += 1.0 / k
  end
  s
end

def print_harmonic(n)
  print "H(#{n}) = #{harmonic(n)}, "
  print "H(#{n})-log(#{n} + 1) = #{harmonic(n) - Math.log(n + 1)}"
  print "\n"
end

(1..6).each do |n|
  print_harmonic(n)
end
print "\n"

(1..6001).step(1000) do |n|
  print_harmonic(n)
end
print "\n"

(1..60001).step(10000) do |n|
  print_harmonic(n)
end
print "\n"

実行結果です。

H(1) = 1.0, H(1)-log(1 + 1) = 0.306852819440055
H(2) = 1.5, H(2)-log(2 + 1) = 0.40138771133189
H(3) = 1.83333333333333, H(3)-log(3 + 1) = 0.447038972213443
H(4) = 2.08333333333333, H(4)-log(4 + 1) = 0.473895420899233
H(5) = 2.28333333333333, H(5)-log(5 + 1) = 0.491573864105278
H(6) = 2.45, H(6)-log(6 + 1) = 0.504089850944687

H(1) = 1.0, H(1)-log(1 + 1) = 0.306852819440055
H(1001) = 7.48646986154934, H(1001)-log(1001 + 1) = 0.576716579904534
H(2001) = 8.17886785373522, H(2001)-log(2001 + 1) = 0.576965893860056
H(3001) = 8.58408311221842, H(3001)-log(3001 + 1) = 0.577049100025008
H(4001) = 8.87164023731082, H(4001)-log(4001 + 1) = 0.57709072216714
H(5001) = 9.0947088129924, H(5001)-log(5001 + 1) = 0.577115701554838
H(6001) = 9.27698038302374, H(6001)-log(6001 + 1) = 0.577132357023428

H(1) = 1.0, H(1)-log(1 + 1) = 0.306852819440055
H(10001) = 9.78770602604535, H(10001)-log(10001 + 1) = 0.577165674066498
H(20001) = 10.4807782147294, H(20001)-log(20001 + 1) = 0.577190667192985
H(30001) = 10.8862183243422, H(30001)-log(30001 + 1) = 0.577198999253353
H(40001) = 11.1738878973205, H(40001)-log(40001 + 1) = 0.577203165474423
H(50001) = 11.3970239488785, H(50001)-log(50001 + 1) = 0.577205665268208
H(60001) = 11.5793405058048, H(60001)-log(60001 + 1) = 0.577207331822818

次に、BigDecimalを使って。logの有効桁数は10にしました(適当)。

require 'bigdecimal'
require "bigdecimal/math"
require 'bigdecimal/util'
include BigMath

def harmonic(n)
  s = BigDecimal("0.0")
  (1..n).each do |k|
    s += BigDecimal("1.0") / BigDecimal(k.to_s)
  end
  s
end

def print_harmonic(n)
  print "H(#{n}) = #{harmonic(n)}, "
  print "H(#{n})-log(#{n} + 1) = #{harmonic(n) - BigMath.log(BigDecimal((n + 1).to_s), 10)}"
  print "\n"
end

(1..6).each do |n|
  print_harmonic(n)
end
print "\n"

(1..6001).step(1000) do |n|
  print_harmonic(n)
end
print "\n"

(1..60001).step(10000) do |n|
  print_harmonic(n)
end
print "\n"

実行結果です。

H(1) = 0.1E1, H(1)-log(1 + 1) = 0.306852819440054690582767887411396223645312E0
H(2) = 0.15E1, H(2)-log(2 + 1) = 0.401387711331890308604754774273176076453984E0
H(3) = 0.18333333333333333E1, H(3)-log(3 + 1) = 0.447038972213442681165535776124198204664712E0
H(4) = 0.20833333333333333E1, H(4)-log(4 + 1) = 0.473895420899232925399240652011540608564226E0
H(5) = 0.22833333333333333E1, H(5)-log(5 + 1) = 0.491573864105278299187522651996827171163856E0
H(6) = 0.245E1, H(6)-log(6 + 1) = 0.50408985094468669489464725551705015967369E0

H(1) = 0.1E1, H(1)-log(1 + 1) = 0.306852819440054690582767887411396223645312E0
H(1001) = 0.74864698615493461E1, H(1001)-log(1001 + 1) = 0.57671657990453599192781827686703277069014E0
H(2001) = 0.81788678537352194E1, H(2001)-log(2001 + 1) = 0.57696589386005350536207258596906317034144E0
H(3001) = 0.85840831122184352E1, H(3001)-log(3001 + 1) = 0.57704910002502793640094159265678243970808E0
H(4001) = 0.88716402373108505E1, H(4001)-log(4001 + 1) = 0.57709072216717178119768641380977576329254E0
H(5001) = 0.90947088129924366E1, H(5001)-log(5001 + 1) = 0.57711570155487223796485957412710820589396E0
H(6001) = 0.92769803830237625E1, H(6001)-log(6001 + 1) = 0.5771323570234500759406694327874107902215E0

H(1) = 0.1E1, H(1)-log(1 + 1) = 0.306852819440054690582767887411396223645312E0
H(10001) = 0.978770602604538530001E1, H(10001)-log(10001 + 1) = 0.57716567406653629720786861251604258541684E0
H(20001) = 0.1048077821472945570242135E2, H(20001)-log(20001 + 1) = 0.5771906671929943485978231877599297030072E0
H(30001) = 0.10886218324342161756084465E2, H(30001)-log(30001 + 1) = 0.57719899925332612368022214899526109094976E0
H(40001) = 0.11173887897320541642696126E2, H(40001)-log(40001 + 1) = 0.57720316547442662268733509118070030868322E0
H(50001) = 0.11397023948878493773713902E2, H(50001)-log(50001 + 1) = 0.57720566526818933035036277169505698709286E0
H(60001) = 0.11579340505804852437056151E2, H(60001)-log(60001 + 1) = 0.57720733182282457702662164497600665290402E0