#!/usr/koeki/bin/ruby
# coding: euc-jp


include Math

def factorial(k)
  if k == 1
    1
  else
    k * factorial(k-1)
  end
end

def lnfac(k)
  j = 1
  s = 0
  while j <= k 
    s += log(j)
    j += 1
  end
  s
end

def stirling(k)
#  0.5 * log(2 * PI) + 0.5 * log(k) + 
  k * (log(k) - 1)
end

def stirling1(k)
0.5 * log(k) + k * (log(k) - 1)
end

def stirling2(k)
0.5 * log(2 * PI) + 0.5 * log(k) + k * (log(k) - 1)
end


printf("%7s %10s\t%8s %14s\n","m!","計算結果","近似式","誤差[%]")
for d in 1 .. 20
  m = 5 * d
  printf("%6d! %15e\t%15e %6.1f\n", 
         m, 
         factorial(m),
         (sqrt(2 * PI * m)) * (m / E) ** m,
         (factorial(m) - (sqrt(2 * PI * m)) * (m / E) ** m)/factorial(m).to_f * 100)
end
puts""
printf("%7s %11s\t%9s %15s\n","m!","オーダー","近似式","誤差[%]")
for d in 1 .. 20
  m = 5 * d
  printf("%6d! %15e\t%15e %6.1f\n", 
         m, 
         lnfac(m),
         stirling(m),
         (lnfac(m) - stirling(m))/lnfac(m) * 100
         ) 
end

puts""
printf("%7s %11s\t%9s %15s\n","m!","オーダー","定数を除いた項","誤差[%]")
for d in 1 .. 20
  m = 5 * d
  printf("%6d! %15e\t%15e %6.1f\n", 
         m, 
         lnfac(m),
         stirling1(m),
         (lnfac(m) - stirling1(m))/lnfac(m) * 100
         ) 
end
puts""
printf("%7s %11s\t%9s %15s\n","m!","オーダー","完全な式","誤差[%]")
for d in 1 .. 20
  m = 5 * d
  printf("%6d! %15e\t%15e %6.1f\n", 
         m, 
         lnfac(m),
         stirling2(m),
         (lnfac(m) - stirling2(m))/lnfac(m) * 100
         ) 
end