mirror of
https://github.com/tgorordo/WignerSymbols.jl.git
synced 2026-06-05 15:42:15 -07:00
9j caching to canonical lexicographical order w/ regge sym
This commit is contained in:
parent
8e3d78114e
commit
e21df477fc
GPG key ID:
0CBED22BB0D94490
2 changed files with 76 additions and 46 deletions
|
|
@ -269,19 +269,6 @@ function wigner9j(T::Type{<:Real}, j₁, j₂, j₃, j₄, j₅, j₆, j₇, j
|
|||
return _wigner9j(T, HalfInteger.((j₁, j₂, j₃, j₄, j₅, j₆, j₇, j₈, j₉))...)
|
||||
end
|
||||
|
||||
const _perms9j = [([1, 2, 3], [1, 2, 3]) ([1, 2, 3], [1, 3, 2]) ([1, 2, 3], [2, 1, 3]) ([1, 2, 3], [2, 3, 1]) ([1, 2, 3], [3, 1, 2]) ([1, 2, 3], [3, 2, 1]);
|
||||
([1, 3, 2], [1, 2, 3]) ([1, 3, 2], [1, 3, 2]) ([1, 3, 2], [2, 1, 3]) ([1, 3, 2], [2, 3, 1]) ([1, 3, 2], [3, 1, 2]) ([1, 3, 2], [3, 2, 1]);
|
||||
([2, 1, 3], [1, 2, 3]) ([2, 1, 3], [1, 3, 2]) ([2, 1, 3], [2, 1, 3]) ([2, 1, 3], [2, 3, 1]) ([2, 1, 3], [3, 1, 2]) ([2, 1, 3], [3, 2, 1]);
|
||||
([2, 3, 1], [1, 2, 3]) ([2, 3, 1], [1, 3, 2]) ([2, 3, 1], [2, 1, 3]) ([2, 3, 1], [2, 3, 1]) ([2, 3, 1], [3, 1, 2]) ([2, 3, 1], [3, 2, 1]);
|
||||
([3, 1, 2], [1, 2, 3]) ([3, 1, 2], [1, 3, 2]) ([3, 1, 2], [2, 1, 3]) ([3, 1, 2], [2, 3, 1]) ([3, 1, 2], [3, 1, 2]) ([3, 1, 2], [3, 2, 1]);
|
||||
([3, 2, 1], [1, 2, 3]) ([3, 2, 1], [1, 3, 2]) ([3, 2, 1], [2, 1, 3]) ([3, 2, 1], [2, 3, 1]) ([3, 2, 1], [3, 1, 2]) ([3, 2, 1], [3, 2, 1])]
|
||||
|
||||
const _signs9j = [1 -1 -1 1 1 -1;
|
||||
-1 1 1 -1 -1 1;
|
||||
-1 1 1 -1 -1 1;
|
||||
1 -1 -1 1 1 -1;
|
||||
1 -1 -1 1 1 -1;
|
||||
-1 1 1 -1 -1 1]
|
||||
|
||||
function _wigner9j(T::Type{<:Real}, j₁::HalfInteger, j₂::HalfInteger, j₃::HalfInteger,
|
||||
j₄::HalfInteger, j₅::HalfInteger, j₆::HalfInteger,
|
||||
|
|
@ -292,27 +279,23 @@ function _wigner9j(T::Type{<:Real}, j₁::HalfInteger, j₂::HalfInteger, j₃::
|
|||
return zero(T)
|
||||
end
|
||||
|
||||
# dictionary lookup, check all 72 permutations
|
||||
k = [j₁ j₂ j₃; j₄ j₅ j₆; j₇ j₈ j₉]
|
||||
for (p, m) in zip(_perms9j, _signs9j)
|
||||
kk = Tuple(reshape(k[p...], 9))
|
||||
kkT = Tuple(reshape(transpose(k[p...]), 9))
|
||||
if haskey(Wigner9j, kk)
|
||||
r, s = Wigner9j[kk]
|
||||
return m * _convert(T, s) * convert(T, signedroot(r))
|
||||
elseif haskey(Wigner9j, kkT)
|
||||
r, s = Wigner9j[kkT]
|
||||
return m * _convert(T, s) * convert(T, signedroot(r))
|
||||
end
|
||||
# dictionary lookup by regge symmetries
|
||||
k, m = reorder9j(j₁, j₂, j₃, j₄, j₅, j₆, j₇, j₈, j₉)
|
||||
if haskey(Wigner9j, k)
|
||||
r, s = Wigner9j[k]
|
||||
return m * _convert(T, s) * convert(T, signedroot(r))
|
||||
end
|
||||
|
||||
# canonicalized entries
|
||||
j₁c, j₂c, j₃c, j₄c, j₅c, j₆c, j₇c, j₈c, j₉c = k
|
||||
|
||||
# order irrelevant: product remains the same
|
||||
n₁, d₁ = Δ²(j₁, j₂, j₃)
|
||||
n₂, d₂ = Δ²(j₄, j₅, j₆)
|
||||
n₃, d₃ = Δ²(j₇, j₈, j₉)
|
||||
n₄, d₄ = Δ²(j₁, j₄, j₇)
|
||||
n₅, d₅ = Δ²(j₂, j₅, j₈)
|
||||
n₆, d₆ = Δ²(j₃, j₆, j₉)
|
||||
n₁, d₁ = Δ²(j₁c, j₂c, j₃c)
|
||||
n₂, d₂ = Δ²(j₄c, j₅c, j₆c)
|
||||
n₃, d₃ = Δ²(j₇c, j₈c, j₉c)
|
||||
n₄, d₄ = Δ²(j₁c, j₄c, j₇c)
|
||||
n₅, d₅ = Δ²(j₂c, j₅c, j₈c)
|
||||
n₆, d₆ = Δ²(j₃c, j₆c, j₉c)
|
||||
|
||||
snum, rnum = splitsquare(n₁ * n₂ * n₃ * n₄ * n₅ * n₆)
|
||||
sden, rden = splitsquare(d₁ * d₂ * d₃ * d₄ * d₅ * d₆)
|
||||
|
|
@ -321,7 +304,7 @@ function _wigner9j(T::Type{<:Real}, j₁::HalfInteger, j₂::HalfInteger, j₃::
|
|||
rr = Base.unsafe_rational(convert(BigInt, rnum), convert(BigInt, rden))
|
||||
|
||||
snum, sden = divgcd!(snum, sden)
|
||||
snum2 = compute9jseries(j₁, j₂, j₃, j₄, j₅, j₆, j₇, j₈, j₉)
|
||||
snum2 = compute9jseries(j₁c, j₂c, j₃c, j₄c, j₅c, j₆c, j₇c, j₈c, j₉c)
|
||||
for n = 1:length(sden.powers)
|
||||
p = bigprime(n)
|
||||
while sden.powers[n] > 0
|
||||
|
|
@ -336,8 +319,8 @@ function _wigner9j(T::Type{<:Real}, j₁::HalfInteger, j₂::HalfInteger, j₃::
|
|||
end
|
||||
s = Base.unsafe_rational(convert(BigInt, snum)*snum2, convert(BigInt, sden))
|
||||
|
||||
Wigner9j[(j₁, j₂, j₃, j₄, j₅, j₆, j₇, j₈, j₉)] = (rr, s)
|
||||
return _convert(T, s) * convert(T, signedroot(rr))
|
||||
Wigner9j[(j₁c, j₂c, j₃c, j₄c, j₅c, j₆c, j₇c, j₈c, j₉c)] = (rr, s)
|
||||
return m * _convert(T, s) * convert(T, signedroot(rr))
|
||||
end
|
||||
|
||||
# COMPUTATIONAL ROUTINES
|
||||
|
|
@ -392,6 +375,49 @@ function reorder6j(β₁, β₂, β₃, α₁, α₂, α₃, α₄)
|
|||
end
|
||||
end
|
||||
|
||||
const _perms = [[1, 2, 3], [1, 3, 2], [2, 1, 3], [2, 3, 1], [3, 1, 2], [3, 2, 1]]
|
||||
const _perms9j = [(r, c) for r in _perms, c in _perms]
|
||||
_permparity(p) = count(i < j && p[i] > p[j] for i in 1:length(p), j in 1:length(p)) % 2
|
||||
|
||||
#const _perms9j = [([1, 2, 3], [1, 2, 3]) ([1, 2, 3], [1, 3, 2]) ([1, 2, 3], [2, 1, 3]) ([1, 2, 3], [2, 3, 1]) ([1, 2, 3], [3, 1, 2]) ([1, 2, 3], [3, 2, 1]);
|
||||
# ([1, 3, 2], [1, 2, 3]) ([1, 3, 2], [1, 3, 2]) ([1, 3, 2], [2, 1, 3]) ([1, 3, 2], [2, 3, 1]) ([1, 3, 2], [3, 1, 2]) ([1, 3, 2], [3, 2, 1]);
|
||||
# ([2, 1, 3], [1, 2, 3]) ([2, 1, 3], [1, 3, 2]) ([2, 1, 3], [2, 1, 3]) ([2, 1, 3], [2, 3, 1]) ([2, 1, 3], [3, 1, 2]) ([2, 1, 3], [3, 2, 1]);
|
||||
# ([2, 3, 1], [1, 2, 3]) ([2, 3, 1], [1, 3, 2]) ([2, 3, 1], [2, 1, 3]) ([2, 3, 1], [2, 3, 1]) ([2, 3, 1], [3, 1, 2]) ([2, 3, 1], [3, 2, 1]);
|
||||
# ([3, 1, 2], [1, 2, 3]) ([3, 1, 2], [1, 3, 2]) ([3, 1, 2], [2, 1, 3]) ([3, 1, 2], [2, 3, 1]) ([3, 1, 2], [3, 1, 2]) ([3, 1, 2], [3, 2, 1]);
|
||||
# ([3, 2, 1], [1, 2, 3]) ([3, 2, 1], [1, 3, 2]) ([3, 2, 1], [2, 1, 3]) ([3, 2, 1], [2, 3, 1]) ([3, 2, 1], [3, 1, 2]) ([3, 2, 1], [3, 2, 1])]
|
||||
|
||||
# reorder parameters determining the 9j symbol to canonical (lexicographic) order with regge factor
|
||||
function reorder9j(j₁, j₂, j₃, j₄, j₅, j₆, j₇, j₈, j₉)
|
||||
|
||||
J = [j₁ j₂ j₃; j₄ j₅ j₆; j₇ j₈ j₉]
|
||||
|
||||
k = (j₁, j₂, j₃, j₄, j₅, j₆, j₇, j₈, j₉)
|
||||
m = 1
|
||||
|
||||
rg = (-1)^sum(k)
|
||||
|
||||
for p in _perms9j
|
||||
JJ = J[p...]
|
||||
kk = Tuple(reshape(JJ, 9))
|
||||
kkT = Tuple(reshape(JJ', 9))
|
||||
|
||||
s = iseven(_permparity(p[1]) + _permparity(p[2])) ? 1 : rg
|
||||
|
||||
if kk < kkT
|
||||
if kk < k
|
||||
k = kk
|
||||
m = s
|
||||
end
|
||||
else # kk > kkT
|
||||
if kkT < k
|
||||
k = kkT
|
||||
m = s
|
||||
end
|
||||
end
|
||||
end
|
||||
return k, m
|
||||
end
|
||||
|
||||
# compute the sum appearing in the 3j symbol
|
||||
function compute3jseries(β₁, β₂, β₃, α₁, α₂)
|
||||
krange = max(α₁, α₂, zero(α₁)):min(β₁, β₂, β₃)
|
||||
|
|
|
|||
|
|
@ -208,21 +208,25 @@ end
|
|||
end
|
||||
end
|
||||
|
||||
|
||||
@threads for i = 1:N
|
||||
@testset "wigner9j: relation to sum over 6j products, thread $i" begin
|
||||
for k = 1:10_000
|
||||
let (j1, j2, j3, j4, j5, j6, j7, j8, j9) = rand(smalljlist, 9)
|
||||
@test wigner9j(j1, j2, j3,
|
||||
j4, j5, j6,
|
||||
j7, j8, j9) ≈ sum(largejlist) do x # lazy choice for range of this sum, but good enough
|
||||
(iseven(2x) ? (2x + 1) : -(2x + 1)) *
|
||||
wigner6j(j1, j4, j7,
|
||||
j8, j9, x ) *
|
||||
wigner6j(j2, j5, j8,
|
||||
j4, x , j6) *
|
||||
wigner6j(j3, j6, j9,
|
||||
x , j1, j2)
|
||||
end
|
||||
debug = false
|
||||
js = rand(smalljlist, 9)
|
||||
|
||||
debug |= !(@test(wigner9j(js...) ≈ sum(largejlist) do x # lazy choice for range of this sum, but good enough
|
||||
(iseven(2x) ? (2x + 1) : -(2x + 1)) *
|
||||
wigner6j(js[1], js[4], js[7],
|
||||
js[8], js[9], x ) *
|
||||
wigner6j(js[2], js[5], js[8],
|
||||
js[4], x , js[6]) *
|
||||
wigner6j(js[3], js[6], js[9],
|
||||
x , js[1], js[2])
|
||||
end) isa Test.Pass)
|
||||
|
||||
if debug
|
||||
@warn "Failed on input" js
|
||||
end
|
||||
end
|
||||
end
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue