読者です 読者をやめる 読者になる 読者になる

fortran66のブログ

fortran について書きます。

Python で 3j-symbol

Python は確かに楽に有理数・多倍長計算ができる。

フォーマット文にあたるものがよくわからなくて整形出力できない・・・IPython では TeX 表記もできる様だから、数学記号的に出力させられるのではないかとも思うのだが、よく分からんw

In [9]: chk3(4,4,4, 4, 0, -4)
0.104297703129124 sqrt(2002)/429 ({2: 1, 7: 1}, {11: 1, 3: 2, 13: 1})

In [10]: chk3(4,4,2, 1, -1, 0)
0.144400037573994 17*sqrt(385)/2310 ({17: 2}, {11: 1, 2: 2, 3: 2, 5: 1, 7: 1})

In [11]: chk3(3.5,3,1.5, 1.5, -1, -0.5)
0 0 ({0: 1}, {})

余りチェックしていない・・・

pylab 環境下で動作。

import sympy
import numpy

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

def threeJ(j1, j2, j3, m1, m2, m3):
    j1 = sympy.Rational(j1)
    j2 = sympy.Rational(j2)
    j3 = sympy.Rational(j3)
    m1 = sympy.Rational(m1)
    m2 = sympy.Rational(m2)
    m3 = sympy.Rational(m3)
    if m1 + m2 + m3 <> 0: return 0
    if j1 + j2 < j3: return 0
    if j2 + j3 < j1: return 0
    if j3 + j1 < j2: return 0
    J = j1 + j2 + j3 + 1
    jj3 =  j1 + j2 - j3
    jj2 =  j1 - j2 + j3 
    jj1 = -j1 + j2 + j3
    jm1 = j1 - m1
    jp1 = j1 + m1
    jm2 = j2 - m2
    jp2 = j2 + m2
    jm3 = j3 - m3
    jp3 = j3 + m3
    phase  = (-1)**(j1 - j2 - m3)
    factor2 = sympy.Rational(fa(jj3) * fa(jj2) * fa(jj1) *
              fa(jp1) * fa(jm1) * fa(jp2) * fa(jm2) * fa(jp3) * fa(jm3), fa(J))
    factor = sympy.sqrt(factor2)
#
    jm231 = j2 - j3 - m1
    jm132 = j1 - j3 + m2
    kmin = max(0, jm231, jm132)
    kmax = min(jj3, jm1, jp2)
    s = 0
    for k in arange(kmin, kmax + 1):
        s = s + sympy.Rational( (-1)**k, fa(k) * fa(jj3 - k) * fa(jm1 - k) * fa(jp2 - k) * fa(k - jm231) * fa(k - jm132) ) 
    return phase * factor * s

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(fa(2 * k1) * fa(2 * k2) * fa(2 * k3), fa(J2 + 1))
    t = sympy.Rational(fa(J), fact(k1) * fa(k2) * fa(k3))
    return (-1)**J * sympy.sqrt(s) * t

def dicsort(d):
    return sorted(d.items(), key=lambda x: x[0])

def pr(i):
    a = sympy.factorint((i * i).p)
    b = sympy.factorint((i * i).q)
    return a, b #dicsort(a), dicsort(b)
    
def chk0(j1, j2, j3):
    x = threeJ0(j1, j2, j3)
    print x.evalf(), x, pr(x)

def chk3(j1, j2, j3, m1, m2, m3):
    x = threeJ(j1, j2, j3, m1, m2, m3)
    print x.evalf(), x, pr(x)

def chk03(j1, j2, j3):
    print chk0(j1, j2, j3)
    print chk3( j1, j2, j3, 0, 0, 0)
    return