diff --git a/src/WignerSymbols.jl b/src/WignerSymbols.jl index f6ec63d..22fff82 100644 --- a/src/WignerSymbols.jl +++ b/src/WignerSymbols.jl @@ -1,9 +1,9 @@ module WignerSymbols -export wigner3j, wigner6j +export δ, Δ, clebschgordan, wigner3j, wigner6j include("primefactorization.jl") -const Wigner3j = Dict{NTuple{5,UInt},Tuple{Rational{BigInt},Rational{BigInt}}}() +const Wigner3j = Dict{Tuple{UInt,UInt,UInt,Int,Int},Tuple{Rational{BigInt},Rational{BigInt}}}() const Wigner6j = Dict{NTuple{6,UInt},Tuple{Rational{BigInt},Rational{BigInt}}}() clebschgordan(j₁, m₁, j₂, m₂, j₃, m₃ = m₁+m₂) = clebschgordan(Float64, j₁, m₁, j₂, m₂, j₃, m₃) @@ -21,19 +21,35 @@ wigner3j(j₁, j₂, j₃, m₁, m₂, m₃ = -m₂-m₃) = wigner3j(Float64, j function wigner3j(T::Type{<:AbstractFloat}, j₁, j₂, j₃, m₁, m₂, m₃ = -m₂-m₃) # check angular momenta and triangle condition if !(δ(j₁, j₂, j₃) && ϵ(j₁, m₁) && ϵ(j₂, m₂) && ϵ(j₃, m₃)) - throw(DomainError()) + throw(DomainError((j₁, j₂, j₃, m₁, m₂, m₃))) end iszero(m₁+m₂+m₃) || return zero(T) # we reorder such that j₁ >= j₂ >= j₃ and m₁ >= 0 or m₁ == 0 && m₂ >= 0 j₁, j₂, j₃, m₁, m₂, m₃, sgn = reorder3j(j₁, j₂, j₃, m₁, m₂, m₃) + # do we also want to use Regge symmetries? + α₁ = convert(Int, j₂ - m₁ - j₃ ) # can be negative + α₂ = convert(Int, j₁ + m₂ - j₃ ) # can be negative + β₁ = convert(Int, j₁ + j₂ - j₃ ) + β₂ = convert(Int, j₁ - m₁ ) + β₃ = convert(Int, j₂ + m₂ ) + + # extra sign in definition + sgn = isodd(j₁ - j₃ - m₃) ? -sgn : sgn # dictionary lookup or compute - if haskey(Wigner3j, (j₁, j₂, j₃, m₁, m₂)) - r, s = Wigner3j[(j₁, j₂, j₃, m₁, m₂)] + if haskey(Wigner3j, (β₁, β₂, β₃, α₁, α₂)) + r, s = Wigner3j[(β₁, β₂, β₃, α₁, α₂)] else - r, s = compute3j(j₁, j₂, j₃, m₁, m₂) - Wigner3j[(j₁, j₂, j₃, m₁, m₂)] = (r,s) + s1 = Δ²(j₁, j₂, j₃) + s2 = prod(map(primefactorial, (β₂, β₁ - α₁, β₁ - α₂, β₃, β₃ - α₁, β₂ - α₂))) + + snum, rnum = splitsquare(s1.num * s2) + sden, rden = splitsquare(s1.den) + s = convert(BigInt, snum) // convert(BigInt, sden) + r = convert(BigInt, rnum) // convert(BigInt, rden) + s *= compute3jseries(β₁, β₂, β₃, α₁, α₂) + Wigner3j[(β₁, β₂, β₃, α₁, α₂)] = (r,s) end return sgn*sqrt(convert(T, r))*convert(T, s) @@ -48,7 +64,8 @@ function wigner6j(T::Type{<:AbstractFloat}, j₁, j₂, j₃, j₄, j₅, j₆) # check triangle conditions if !(δ(α̂₁...) && δ(α̂₂...) && δ(α̂₃...) && δ(α̂₄...)) - throw(DomainError()) + return zero(T) + # throw(DomainError()) end # reduce α₁ = convert(UInt, +(α̂₁...)) @@ -100,6 +117,15 @@ function δ(j₁, j₂, j₃) return true end +# triangle coefficient +function Δ(T::Type{<:AbstractFloat}, j₁, j₂, j₃) + if !δ(j₁, j₂, j₃) + throw(DomainError()) + end + v = Δ²(j₁, j₂, j₃) + return sqrt(convert(T, convert(BigInt, v.num) // convert(BigInt, v.den))) +end + # squared triangle coefficient function Δ²(j₁, j₂, j₃) # also checks the triangle conditions by converting to unsigned integer: @@ -113,11 +139,11 @@ end # reorder parameters determining the 3j symbol to canonical order: # j₁ >= j₂ >= j₃ and m₁ >= 0 or m₁ == 0 && m₂ >= 0 -function reorder3j(j₁, j₂, j₃, m₁, m₂, m₃, sgn = one(UInt8)) +function reorder3j(j₁, j₂, j₃, m₁, m₂, m₃, sign = one(UInt8)) if j₁ < j₂ - return reorder3j(j₂, j₁, j₃, m₁, m₂, m₃, -sign) + return reorder3j(j₂, j₁, j₃, m₂, m₁, m₃, -sign) elseif j₂ < j₃ - return reorder3j(j₁, j₃, j₂, m₁, m₂, m₃, -sign) + return reorder3j(j₁, j₃, j₂, m₁, m₃, m₂, -sign) elseif m₁ < zero(m₁) return reorder3j(j₁, j₂, j₃, -m₁, -m₂, -m₃, -sign) elseif iszero(m₁) && m₂ < zero(m₂) @@ -149,37 +175,42 @@ function reorder6j(β₁, β₂, β₃, α₁, α₂, α₃, α₄) end end -function compute3j(j₁, j₂, j₃, m₁, m₂) - m₃ = -m₁ - m₂ - - α₁ = convert(UInt, j₂ - m₁ - j₃ ) - α₂ = convert(UInt, j₁ + m₂ - j₃ ) - β₁ = convert(UInt, j₁ + j₂ - j₃ ) - β₂ = convert(UInt, j₁ - m₁ ) - β₃ = convert(UInt, j₂ + m₂ ) - - krange = max(α₁,α₂,zero(UInt)):min(β₁,β₂,β₃) - -end - - - -# compute the sum appearing in the 6j symbol -function compute6jseries(β₁, β₂, β₃, α₁, α₂, α₃, α₄) - krange = max(α₁,α₂,α₃,α₄):min(β₁,β₂,β₃) - - nums = Vector{typeof(snum)}(length(krange)) - dens = Vector{typeof(snum)}(length(krange)) +# compute the sum appearing in the 3j symbol +function compute3jseries(β₁, β₂, β₃, α₁, α₂) + krange = max(α₁, α₂, zero(α₁)):min(β₁, β₂, β₃) + T = PrimeFactorization{eltype(eltype(factorialtable))} + nums = Vector{T}(length(krange)) + dens = Vector{T}(length(krange)) for (i, k) in enumerate(krange) - num = iseven(k) ? primefactorial(k+1) : -primefactorial(k+1) - den = primefactorial(k-α₁)*primefactorial(k-α₂)*primefactorial(k-α₃)* - primefactorial(k-α₄)*primefactorial(β₁-k)*primefactorial(β₂-k)*primefactorial(β₃-k) - nums[i], dens[i] = divgcd(num, den) + num = iseven(k) ? primefactorial(1) : -primefactorial(1) + den = primefactorial(k)*primefactorial(k-α₁)*primefactorial(k-α₂)* + primefactorial(β₁-k)*primefactorial(β₂-k)*primefactorial(β₃-k) + nums[i], dens[i] = divgcd!(num, den) end den = commondenominator!(nums, dens) totalnum = sumlist!(nums) totalden = convert(BigInt, PrimeFactorization(den)) + return totalnum//totalden +end + +# compute the sum appearing in the 6j symbol +function compute6jseries(β₁, β₂, β₃, α₁, α₂, α₃, α₄) + krange = max(α₁, α₂, α₃, α₄):min(β₁, β₂, β₃) + T = PrimeFactorization{eltype(eltype(factorialtable))} + + nums = Vector{T}(length(krange)) + dens = Vector{T}(length(krange)) + for (i, k) in enumerate(krange) + num = iseven(k) ? primefactorial(k+1) : -primefactorial(k+1) + den = primefactorial(k-α₁)*primefactorial(k-α₂)*primefactorial(k-α₃)* + primefactorial(k-α₄)*primefactorial(β₁-k)*primefactorial(β₂-k)*primefactorial(β₃-k) + nums[i], dens[i] = divgcd!(num, den) + end + den = commondenominator!(nums, dens) + totalnum = sumlist!(nums) + totalden = convert(BigInt, PrimeFactorization(den)) + return totalnum//totalden end end # module diff --git a/src/primefactorization.jl b/src/primefactorization.jl index 79916c5..19a4693 100644 --- a/src/primefactorization.jl +++ b/src/primefactorization.jl @@ -57,6 +57,7 @@ end PrimeFactorization(powers::Vector{T}) where {T<:Unsigned} = PrimeFactorization{T}(powers, one(Int8)) Base.copy(a::PrimeFactorization) = PrimeFactorization(copy(a.powers), a.sign) +Base.:(==)(a::P, b::P) where {P<:PrimeFactorization} = a.powers == b.powers && a.sign == b.sign # define our own factor function, returning an instance of PrimeFactorization function primefactor(n::Integer) @@ -81,11 +82,11 @@ function primefactor(n::Integer) end push!(factortable, powers) end - @inbounds return PrimeFactorization(factortable[n], sn) + @inbounds return PrimeFactorization(copy(factortable[n]), sn) end function primefactorial(n::Integer) - n < 0 && return DomainError(n) + n < 0 && throw(DomainError(n)) m = length(factorialtable)-1 @inbounds while m < n prevfactorial = factorialtable[m+1] @@ -100,7 +101,7 @@ function primefactorial(n::Integer) end push!(factorialtable, powers) end - @inbounds return PrimeFactorization(factorialtable[n+1]) + @inbounds return PrimeFactorization(copy(factorialtable[n+1])) end # Conversion from PrimeFactorization to `BigInt` using `bigprime` @@ -134,7 +135,7 @@ function Base.gcd(a::PrimeFactorization{T}, b::PrimeFactorization{T}) where {T} elseif b.sign ==0 return a else - return Primefactorization(_vmin!(copy(a.powers), b.powers)) + return PrimeFactorization(_vmin!(copy(a.powers), b.powers)) end end function Base.lcm(a::PrimeFactorization{T}, b::PrimeFactorization{T}) where {T} @@ -143,33 +144,37 @@ function Base.lcm(a::PrimeFactorization{T}, b::PrimeFactorization{T}) where {T} elseif b.sign ==0 return b else - return Primefactorization(_vmax!(copy(a.powers), b.powers)) + return PrimeFactorization(_vmax!(copy(a.powers), b.powers)) end end -function Base.divgcd(a::PrimeFactorization{T}, b::PrimeFactorization{T}) where {T} +Base.divgcd(a::PrimeFactorization{T}, b::PrimeFactorization{T}) where {T} = divgcd!(copy(a), copy(b)) +function divgcd!(a::PrimeFactorization{T}, b::PrimeFactorization{T}) where {T} af, bf = a.powers, b.powers - - ag = copy(af) - bg = copy(bf) - for k = 1:min(length(ag), length(bg)) - gk = min(ag[k], bg[k]) - ag[k] -= gk - bg[k] -= gk + for k = 1:min(length(af), length(bf)) + gk = min(af[k], bf[k]) + af[k] -= gk + bf[k] -= gk end - while length(ag) > 0 && iszero(last(ag)) - pop!(ag) + while length(af) > 0 && iszero(last(af)) + pop!(af) end - while length(bg) > 0 && iszero(last(bg)) - pop!(bg) + while length(bf) > 0 && iszero(last(bf)) + pop!(bf) end - return PrimeFactorization{T}(ag, a.sign), PrimeFactorization{T}(bg, b.sign) + return a, b end # split `a::PrimeFactorization` into a square `s` and a remainder `r`, such that # `a = s^2 * r` and the powers in the prime factorization of `r` are zero or one function splitsquare(a::PrimeFactorization) r = PrimeFactorization(map(p->convert(UInt8, isodd(p)), a.powers), a.sign) + while length(r.powers) > 0 && iszero(last(r.powers)) + pop!(r.powers) + end s = PrimeFactorization(map(p->(p>>1), a.powers)) + while length(s.powers) > 0 && iszero(last(s.powers)) + pop!(s.powers) + end return s, r end @@ -191,8 +196,8 @@ end # auxiliary function to compute sums of a list of PrimeFactorizations as quickly as possible function sumlist!(list::Vector{<:PrimeFactorization}, ind = 1:length(list)) - # depending on length, sum smaller parts first - g = copy(list[ind[1]]) + # first compute gcd to take out common factors + g = PrimeFactorization(copy(list[ind[1]].powers)) for k in ind _vmin!(g.powers, list[k].powers) end @@ -204,8 +209,7 @@ function sumlist!(list::Vector{<:PrimeFactorization}, ind = 1:length(list)) l = L >> 1 s = sumlist!(list, first(ind)+(0:l-1)) + sumlist!(list, first(ind)+(l:L-1)) else - # first compute gcd to take out common factors - + # do sum s = big(0) for k in ind Base.GMP.MPZ.add!(s, convert(BigInt, list[k]))