fortran66のブログ

fortran について書きます。

Python で有理数計算

よくわからなくて苦しいw ライブラリの関係も謎すぎる、どの程度信じていいのか?w

ThreeJ symbol の m=0 の時でチェック。

In [190]: chk(10,6,6)
-0.0884757071422 -6*sqrt(2028117)/96577 ({2: 2, 3: 3, 7: 1}, {17: 1, 19: 1, 13: 1, 23: 1})

Python ソース

import sympy
import numpy

def fact(n):
    if n < 0:
        return 0
    else:
        return math.factorial(n)


def threeJ0(j1, j2, j3):
    J2 = j1 + j2 + j3 
    if J2 % 2 == 1:
        return 0
    J = J2 / 2
    k1 = J - j1
    k2 = J - j2
    k3 = J - j3
    s = sympy.Rational(fact(2 * k1) * fact(2 * k2) * fact(2 * k3), fact(J2 + 1))
    t = sympy.Rational(fact(J), fact(k1) * fact(k2) * fact(k3))
    return (-1)**J * sympy.sqrt(s) * t

def pr(i):
    a = sympy.factorint((i * i).p)
    b = sympy.factorint((i * i).q)
    return a, b

def chk(j1, j2, j3):
    x = threeJ0(j1, j2, j3)
    print x.evalf(), x, pr(x)

値を持たない時の扱いは放置プレイで。