diff --git a/Project.toml b/Project.toml index 2e945cf..02e6f6f 100644 --- a/Project.toml +++ b/Project.toml @@ -1,25 +1,26 @@ 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" +Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" 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..1a7b209 100644 --- a/README.md +++ b/README.md @@ -89,9 +89,6 @@ 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 symbol implementation is 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) diff --git a/src/WignerSymbols.jl b/src/WignerSymbols.jl index f7d820a..29cd56e 100644 --- a/src/WignerSymbols.jl +++ b/src/WignerSymbols.jl @@ -1,11 +1,12 @@ module WignerSymbols -export δ, Δ, clebschgordan, wigner3j, wigner6j, racahV, racahW, HalfInteger +export δ, Δ, clebschgordan, wigner3j, wigner6j, wigner9j, racahV, racahW, HalfInteger using HalfIntegers using RationalRoots using LRUCache const RRBig = RationalRoot{BigInt} import RationalRoots: _convert +using Combinatorics include("growinglist.jl") include("bigint.jl") # additional GMP BigInt functionality not wrapped in Base.GMP.MPZ @@ -14,9 +15,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 +27,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 +248,75 @@ 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 + +_perms9j = [(i, j) for i in Combinatorics.permutations([1, 2, 3]), + j in Combinatorics.permutations([1, 2, 3])] + +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 in _perms9j + kk = tuple(reshape(k[p...], 9)) + kkT = tuple(reshape(transpose(k[p...]), 9)) + if haskey(Wigner9j, kk) + r, s = Wigner9j[kk] + return _convert(T, s) * convert(T, signedroot(r)) + elseif haskey(Wigner9j, kkT) + r, s = Wigner9j[kkT] + return _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₆) + snum, sden = divgcd!(snum, sden) + rnum, rden = divgcd!(rnum, rden) + s = Base.unsafe_rational(convert(BigInt, snum), convert(BigInt, sden)) + r = Base.unsafe_rational(convert(BigInt, rnum), convert(BigInt, rden)) + s *= compute9jseries(j₁, j₂, j₃, j₄, j₅, j₆, j₇, j₈, j₉) + + Wigner9j[(j₁, j₂, j₃, j₄, j₅, j₆, j₇, j₈, j₉)] = (r, s) + return _convert(T, s) * convert(T, signedroot(r)) +end + + # COMPUTATIONAL ROUTINES #------------------------ # squared triangle coefficient @@ -312,7 +387,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 +439,80 @@ 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 + + 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 +end + +export wei9jbracket +# 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..1e7b5ba 100644 --- a/src/primefactorization.jl +++ b/src/primefactorization.jl @@ -384,3 +384,31 @@ function sumlist!(list::Vector{<:PrimeFactorization}, ind = 1:length(list)) end return MPZ.mul!(s, _convert!(i, g)) end + +#= +# A cached binomial implementation +bcache = LRU{Tuple{BigInt, BigInt}, PrimeFactorization}(; maxsize=10^6) +function primebinomial(n::BigInt, k::BigInt) + T = PrimeFactorization{eltype(eltype(factorialtable))} + if k == 0 + return one(T) + end # guard + if haskey(bcache, (n, k)) + return bcache[(n, k)] + else + den = primefactor(k) + num = mul!(copy(primefactor(n + 1 - k)), primebinomial(n, k - 1)) + res = divexact!(num, den) + bcache[(n, k)] = res + return res + end +end +=# + +function primebinomial(n::BigInt, k::BigInt) + num = copy(primefactorial(n)) + den = mul!(copy(primefactorial(k)), primefactorial(n - k)) + return divexact!(num, den) +end + + diff --git a/test/runtests.jl b/test/runtests.jl index 5fd30c2..43da50c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -13,6 +13,7 @@ Random.seed!(1234) smalljlist = 0:1//2:10 largejlist = 0:1//2:1000 + @threads for i = 1:N @testset "triangle coefficient, thread $i" begin for k = i:N:length(smalljlist) @@ -207,3 +208,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 + @testset 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