メモ帳 とりあえず定義式通りに。計算量が少ない変形もあるんですが・・・
import sympy import numpy def fa(n): if n < 0: return 0 else: return math.factorial(n) 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 triangle(a, b, c): if a + b < c: return 0 if b + c < a: return 0 if c + a < b: return 0 a = sympy.Rational(a) b = sympy.Rational(b) c = sympy.Rational(c) s1 = a + b - c s2 = a - b + c s3 = -a + b + c s4 = a + b + c + 1 t = sympy.Rational( fa(s1) * fa(s2) * fa(s3), fa(s4) ) return sympy.sqrt(t) def sixJ(j1, j2, j3, l1, l2, l3): j1 = sympy.Rational(j1) j2 = sympy.Rational(j2) j3 = sympy.Rational(j3) l1 = sympy.Rational(l1) l2 = sympy.Rational(l2) l3 = sympy.Rational(l3) factor = triangle(j1, j2, j3) * triangle(j1, l2, l3) * triangle(l1, j2, l3) * triangle(l1, l2, j3) J = j1 + j2 + j3 jl1 = j1 + l2 + l3 jl2 = l1 + j2 + l3 jl3 = l1 + l2 + j3 jl12 = j1 + j2 + l1 + l2 jl23 = j2 + j3 + l2 + l3 jl31 = j3 + j1 + l3 + l1 kmin = max(J, jl1, jl2, jl3) kmax = min(jl12, jl23, jl31) s = 0 for k in arange(kmin, kmax + 1): s = s + sympy.Rational( (-1)**k * fa(k + 1), fa(k - J) * fa(k - jl1) * fa(k - jl2) * fa(k - jl3) * fa(jl12 -k) * fa(jl23 - k) * fa(jl31 - k) ) return factor * s def chk6(j1, j2, j3, l1, l2, l3): x = sixJ(j1, j2, j3, l1, l2, l3) print x.evalf(), x, pr(x)
