diff --git a/.gitignore b/.gitignore index 8c960ec..312718a 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,4 @@ +Manifest.toml *.jl.cov *.jl.*.cov *.jl.mem diff --git a/Project.toml b/Project.toml index 2e945cf..dedcec9 100644 --- a/Project.toml +++ b/Project.toml @@ -1,25 +1,25 @@ name = "WignerSymbols" uuid = "9f57e263-0b3d-5e2e-b1be-24f2bb48858b" authors = ["Jutho Haegeman"] -version = "2.0.0" +version = "2.1.0" [deps] -RationalRoots = "308eb6b3-cc68-5ff3-9e97-c3c4da4fa681" HalfIntegers = "f0d1745a-41c9-11e9-1dd9-e5d34d218721" -Primes = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae" LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637" +Primes = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae" +RationalRoots = "308eb6b3-cc68-5ff3-9e97-c3c4da4fa681" [compat] -RationalRoots = "0.1 - 0.9" HalfIntegers = "1" -Primes = "0.4 - 0.9" LRUCache = "1.3" +Primes = "0.4 - 0.9" +RationalRoots = "0.1 - 0.9" julia = "1.5" [extras] -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] test = ["Test", "LinearAlgebra", "Random"] diff --git a/README.md b/README.md index b6de15e..8286e22 100644 --- a/README.md +++ b/README.md @@ -52,6 +52,7 @@ While the following function signatures are probably self-explanatory, you can q for them in the Julia REPL to get further details. * `wigner3j(T::Type{<:Real} = RationalRoot{BigInt}, j₁, j₂, j₃, m₁, m₂, m₃ = -m₂-m₁) -> ::T` * `wigner6j(T::Type{<:Real} = RationalRoot{BigInt}, j₁, j₂, j₃, j₄, j₅, j₆) -> ::T` +* `wigner9j(T::Type{<:Real} = RationalRoot{BigInt}, j₁, j₂, j₃, j₄, j₅, j₆, j₇, j₈, j₉) -> ::T` * `clebschgordan(T::Type{<:Real} = RationalRoot{BigInt}, j₁, m₁, j₂, m₂, j₃, m₃ = m₁+m₂) -> ::T` * `racahV(T::Type{<:Real} = RationalRoot{BigInt}, j₁, j₂, j₃, m₁, m₂, m₃ = -m₁-m₂) -> ::T` * `racahW(T::Type{<:Real} = RationalRoot{BigInt}, j₁, j₂, J, j₃, J₁₂, J₂₃) -> ::T` @@ -74,7 +75,7 @@ Largely based on reading the paper (but not the code): with some additional modifications to further improve efficiency for large `j` (angular momenta quantum numbers). -In particular, 3j and 6j symbols are computed exactly, in the format `√(r) * s` where `r` +In particular, 3j, 6j, and 9j symbols are computed exactly, in the format `√(r) * s` where `r` and `s` are exactly computed as `Rational{BigInt}`, using an intermediate representation based on prime number factorization. This exact representation is captured by the `RationalRoot` type. For further calculations, these values probably need to be converted @@ -89,9 +90,8 @@ Also uses ideas from [2] [J. Rasch and A. C. H. Yu, SIAM Journal on Scientific Compututing 25 (2003), 1416–1428](https://doi.org/10.1137/S1064827503422932) -for caching the computed 3j and 6j symbols. +for caching the computed 3j and 6j symbols. Wigner 9-j symbols implemented based on -## Todo -* Wigner 9-j symbols, as explained in [1] and based on +[3] [L. Wei, New formula for 9-j symbols and their direct calculation, Computers in Physics, 12 (1998), 632–634.](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.481.5946&rep=rep1&type=pdf) - [3] [L. Wei, New formula for 9-j symbols and their direct calculation, Computers in Physics, 12 (1998), 632–634.](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.481.5946&rep=rep1&type=pdf) +with binomials additionally optimized using the methods described in [1]. diff --git a/src/WignerSymbols.jl b/src/WignerSymbols.jl index f7d820a..f00d5b5 100644 --- a/src/WignerSymbols.jl +++ b/src/WignerSymbols.jl @@ -1,5 +1,5 @@ module WignerSymbols -export δ, Δ, clebschgordan, wigner3j, wigner6j, racahV, racahW, HalfInteger +export δ, Δ, clebschgordan, wigner3j, wigner6j, wigner9j, racahV, racahW, HalfInteger using HalfIntegers using RationalRoots @@ -14,9 +14,11 @@ convert(BigInt, primefactorial(401)) # trigger compilation and generate some fix const Key3j = Tuple{UInt,UInt,UInt,Int,Int} const Key6j = NTuple{6,UInt} +const Key9j = NTuple{9,HalfInteger} const Wigner3j = LRU{Key3j,Tuple{Rational{BigInt},Rational{BigInt}}}(; maxsize = 10^6) const Wigner6j = LRU{Key6j,Tuple{Rational{BigInt},Rational{BigInt}}}(; maxsize = 10^6) +const Wigner9j = LRU{Key9j,Tuple{Rational{BigInt},Rational{BigInt}}}(; maxsize = 10^6) function set_buffer3j_size!(; maxsize) resize!(Wigner3j; maxsize = maxsize) @@ -24,6 +26,9 @@ end function set_buffer6j_size!(; maxsize) resize!(Wigner6j; maxsize = maxsize) end +function set_buffer9j_size!(; maxsize) + resize!(Wigner9j; maxsize = maxsize) +end # check integerness and correctness of (j,m) angular momentum ϵ(j, m) = (abs(m) <= j && ishalfinteger(j) && isinteger(j-m) && isinteger(j+m)) @@ -242,6 +247,99 @@ function racahW(T::Type{<:Real}, j₁, j₂, J, j₃, J₁₂, J₂₃) end end +""" + wigner9j(T::Type{<:Real} = RationalRoot{BigInt}, j₁, j₂, j₃, j₄, j₅, j₆, j₇, j₈, j₉) -> ::T + +Compute the value of the Wigner-9j symbol +⎧ j₁ j₂ j₃ ⎫ +⎨ j₄ j₅ j₆ ⎬ +⎩ j₇ j₈ j₉ ⎭ +as a type `T` real number. By default `T = RationalRoot{BigInt}`. + +Returns `zero(T)` if any of the triangle conditions TODO +are not satisfied, but throws a `DomainError` if +the `jᵢ`s are not integer or halfinteger. +""" +wigner9j(j₁, j₂, j₃, j₄, j₅, j₆, j₇, j₈, j₉) = wigner9j(RRBig, j₁, j₂, j₃, j₄, j₅, j₆, j₇, j₈, j₉) +function wigner9j(T::Type{<:Real}, j₁, j₂, j₃, j₄, j₅, j₆, j₇, j₈, j₉) + # check validity of `jᵢ`s + for jᵢ in (j₁, j₂, j₃, j₄, j₅, j₆, j₇, j₈, j₉) + (ishalfinteger(jᵢ) && jᵢ >= zero(jᵢ)) || throw(DomainError("invalid jᵢ", jᵢ)) + end + 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, + j₇::HalfInteger, j₈::HalfInteger, j₉::HalfInteger) + + # check triangle conditions + if !(δ(j₁, j₂, j₃) && δ(j₄, j₅, j₆) && δ(j₇, j₈, j₉) && δ(j₁, j₄, j₇) && δ(j₂, j₅, j₈) && δ(j₃, j₆, 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 + end + + # 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₉) + + snum, rnum = splitsquare(n₁ * n₂ * n₃ * n₄ * n₅ * n₆) + sden, rden = splitsquare(d₁ * d₂ * d₃ * d₄ * d₅ * d₆) + + rnum, rden = divgcd!(rnum, rden) + 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₉) + for n = 1:length(sden.powers) + p = bigprime(n) + while sden.powers[n] > 0 + q, r = divrem(snum2, p) + if iszero(r) + snum2 = q + sden.powers[n] -= 1 + else + break + end + end + 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)) +end + # COMPUTATIONAL ROUTINES #------------------------ # squared triangle coefficient @@ -312,7 +410,6 @@ function compute3jseries(β₁, β₂, β₃, α₁, α₂) end den = commondenominator!(nums, dens) totalnum = sumlist!(nums) - totalden = convert(BigInt, den) for n = 1:length(den.powers) p = bigprime(n) while den.powers[n] > 0 @@ -365,13 +462,84 @@ function compute6jseries(β₁, β₂, β₃, α₁, α₂, α₃, α₄) return Base.unsafe_rational(totalnum, totalden) end +# compute the sum appearing in the 9j symbol +function compute9jseries(a, b, c, d, e, f, g, h, j) + I1 = max(abs(h - d), abs(b - f), abs(a - j)) + I2 = min(h + d, b + f, a + j) + krange = I1:I2 + + s = sum(krange) do k + p = iseven(2k) ? big(2k + 1) : -big(2k + 1) + + b₁ = let (m₁, m₂, m₃, m₄, m₅, m₆) = (a, b, c, f, j, k) + α₁ = convert(BigInt, m₁ + m₅ - m₆) + α₂ = convert(BigInt, m₁ - m₅ + m₆) + α₃ = convert(BigInt, -m₁ + m₅ + m₆) + β₁ = convert(BigInt, m₁ + m₅ + m₆) + β₂ = convert(BigInt, m₂ + m₄ + m₆) + β₃ = convert(BigInt, m₃ + m₄ + m₅) + β₀ = convert(BigInt, m₁ + m₂ + m₃) + wei9jbracket(α₁, α₂, α₃, β₁, β₂, β₃, β₀) + end + + b₂ = let (m₁, m₂, m₃, m₄, m₅, m₆) = (f, d, e, h, b, k) + α₁ = convert(BigInt, m₁ + m₅ - m₆) + α₂ = convert(BigInt, m₁ - m₅ + m₆) + α₃ = convert(BigInt, -m₁ + m₅ + m₆) + β₁ = convert(BigInt, m₁ + m₅ + m₆) + β₂ = convert(BigInt, m₂ + m₄ + m₆) + β₃ = convert(BigInt, m₃ + m₄ + m₅) + β₀ = convert(BigInt, m₁ + m₂ + m₃) + wei9jbracket(α₁, α₂, α₃, β₁, β₂, β₃, β₀) + end + + b₃ = let (m₁, m₂, m₃, m₄, m₅, m₆) = (h, j, g, a, d, k) + α₁ = convert(BigInt, m₁ + m₅ - m₆) + α₂ = convert(BigInt, m₁ - m₅ + m₆) + α₃ = convert(BigInt, -m₁ + m₅ + m₆) + β₁ = convert(BigInt, m₁ + m₅ + m₆) + β₂ = convert(BigInt, m₂ + m₄ + m₆) + β₃ = convert(BigInt, m₃ + m₄ + m₅) + β₀ = convert(BigInt, m₁ + m₂ + m₃) + wei9jbracket(α₁, α₂, α₃, β₁, β₂, β₃, β₀) + end + + p * b₁ * b₂ * b₃ + end + + return convert(BigInt, s) +end + +# Wei square bracket terms appearing the 9j series +function wei9jbracket(α₁::BigInt, α₂::BigInt, α₃::BigInt, β₁::BigInt, β₂::BigInt, β₃::BigInt, β₀::BigInt) + p = max(β₁, β₂, β₃, β₀) + q = min(α₁ + β₂, α₂ + β₃, α₃ + β₀) + lrange = p:q + T = PrimeFactorization{eltype(eltype(factorialtable))} + + terms = Vector{T}(undef, length(lrange)) + for (j, l) in enumerate(lrange) + b₁ = iseven(l) ? copy(primebinomial(big(l + 1), l - β₁)) : neg!(copy(primebinomial(big(l + 1), l - β₁))) + b₂ = primebinomial(α₁, l - β₂) + b₃ = primebinomial(α₂, l - β₃) + b₄ = primebinomial(α₃, l - β₀) + + terms[j] = mul!(mul!(mul!(b₁, b₂), b₃), b₄) + end + total = sumlist!(terms) + return total +end + function _precompile_() @assert precompile(wigner3j, (Type{Float64}, Int, Int, Int, Int, Int, Int)) @assert precompile(wigner6j, (Type{Float64}, Int, Int, Int, Int, Int, Int)) + @assert precompile(wigner9j, (Type{Float64}, Int, Int, Int, Int, Int, Int, Int, Int, Int)) @assert precompile(wigner3j, (Type{BigFloat}, Rational{Int}, Rational{Int}, Rational{Int}, Rational{Int}, Rational{Int}, Rational{Int})) @assert precompile(wigner6j, (Type{BigFloat}, Rational{Int}, Rational{Int}, Rational{Int}, Rational{Int}, Rational{Int}, Rational{Int})) + @assert precompile(wigner9j, (Type{BigFloat}, Rational{Int}, Rational{Int}, Rational{Int}, Rational{Int}, Rational{Int}, Rational{Int}, Rational{Int}, Rational{Int}, Rational{Int})) @assert precompile(wigner3j, (Type{RRBig}, HalfInt, HalfInt, HalfInt, HalfInt, HalfInt, HalfInt)) @assert precompile(wigner6j, (Type{RRBig}, HalfInt, HalfInt, HalfInt, HalfInt, HalfInt, HalfInt)) + @assert precompile(wigner9j, (Type{RRBig}, HalfInt, HalfInt, HalfInt, HalfInt, HalfInt, HalfInt, HalfInt, HalfInt, HalfInt)) end _precompile_() diff --git a/src/primefactorization.jl b/src/primefactorization.jl index f1fab0e..0697f23 100644 --- a/src/primefactorization.jl +++ b/src/primefactorization.jl @@ -384,3 +384,10 @@ function sumlist!(list::Vector{<:PrimeFactorization}, ind = 1:length(list)) end return MPZ.mul!(s, _convert!(i, g)) end + +function primebinomial(n::BigInt, k::BigInt) + num = copy(primefactorial(n)) + divexact!(num, primefactorial(k)) + return divexact!(num, primefactorial(n-k)) +end + diff --git a/test/runtests.jl b/test/runtests.jl index 5fd30c2..f583fbc 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -207,3 +207,23 @@ end 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 + end + end + end +end