daruma3940の日記

理解や文章に間違い等あればどんなことでもご指摘お願いします

(2^n n!)^{-1/2}のオーバーフローを起こさない計算法

  (2^n n!)^{-1/2}
こういう計算よくありますよね

これをpythonで愚直に計算してみましょう

import numpy as np
import math

nx = 100

ans = (((2 ** nx) * (math.factorial(nx))) ** (-0.5))

print(ans)

かっこが多すぎる糞コードですね

実行結果:
9.193870653018146e-95

nx=100の場合はちゃんと計算できているみたいですがnx=200にするとどうでしょうか

 ans = (((2 ** nx) * (math.factorial(nx))) ** (-0.5))
OverflowError: int too large to convert to float

あららオーバーフローしちゃいましたね

ここで
  a^x = exp( xln(a) )
を用いて工夫してみましょう
  a \equiv 2^n n!
とすると
  ln(a) \Rightarrow ln(2^n n!) =nln(2)+\sum_n ln( n)

  a^{-1/2} \Rightarrow (2^n n!)^{-1/2} =exp(- \frac{1}{2}(nln(2)+\sum_n ln( n) ) )

import numpy as np
import math

nx = 100

ans = (((2 ** nx) * (math.factorial(nx))) ** (-0.5))

def lnfact(n):
    sum = 0
    for i in np.arange(1, n+1, 1):
        sum += np.log(i)
    return sum


ans2 = np.exp( (-0.5) * (nx * np.log(2) + lnfact(nx) )  )
print(ans, ans2)

実行結果
9.193870653018146e-95 9.19387065302e-95

一致してますね
n=200の場合も計算してみましょう

import numpy as np
import math

nx = 200

#ans = (((2 ** nx) * (math.factorial(nx))) ** (-0.5))

def lnfact(n):
    sum = 0
    for i in np.arange(1, n+1, 1):
        sum += np.log(i)
    return sum


ans2 = np.exp( (-0.5) * (nx * np.log(2) + lnfact(nx) )  )
#print(ans, ans2)
print(ans2)

実行結果
2.80902829347e-218

計算が出来ましたやったね!


(研究が行き詰っているからって過去に学んだ内容をブログに書くのに逃げるな)

追記

  (2^n n!)^{-1/2}
が必要な場合っていうのは
調和振動子 - Wikipedia
の量子的な調和振動子のエネルギー固有状態とかです

次回
daruma3940.hatenablog.com