updates

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₂))
+    if haskey(Wigner3j, (β₁, β₂, β₃, α₁, α₂))
-        r, s = Wigner3j[(j₁, j₂, j₃, m₁, m₂)]
+        r, s = Wigner3j[(β₁, β₂, β₃, α₁, α₂)]
     else
-        r, s = compute3j(j₁, j₂, j₃, m₁, m₂)
+        s1 = Δ²(j₁, j₂, j₃)
-        Wigner3j[(j₁, j₂, j₃, m₁, m₂)] = (r,s)
+        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₂)
+# compute the sum appearing in the 3j symbol
-    m₃ = -m₁ - m₂
+function compute3jseries(β₁, β₂, β₃, α₁, α₂)
-
+    krange = max(α₁, α₂, zero(α₁)):min(β₁, β₂, β₃)
-    α₁ = convert(UInt, j₂ - m₁ - j₃ )
+    T = PrimeFactorization{eltype(eltype(factorialtable))}
-    α₂ = 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))
 
+    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)
+        num = iseven(k) ? primefactorial(1) : -primefactorial(1)
-        den = primefactorial(k-α₁)*primefactorial(k-α₂)*primefactorial(k-α₃)*
+        den = primefactorial(k)*primefactorial(k-α₁)*primefactorial(k-α₂)*
-            primefactorial(k-α₄)*primefactorial(β₁-k)*primefactorial(β₂-k)*primefactorial(β₃-k)
+            primefactorial(β₁-k)*primefactorial(β₂-k)*primefactorial(β₃-k)
-        nums[i], dens[i] = divgcd(num, den)
+        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

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
-
+    for k = 1:min(length(af), length(bf))
-    ag = copy(af)
+        gk = min(af[k], bf[k])
-    bg = copy(bf)
+        af[k] -= gk
-    for k = 1:min(length(ag), length(bg))
+        bf[k] -= gk
-        gk = min(ag[k], bg[k])
-        ag[k] -= gk
-        bg[k] -= gk
     end
-    while length(ag) > 0 && iszero(last(ag))
+    while length(af) > 0 && iszero(last(af))
-        pop!(ag)
+        pop!(af)
     end
-    while length(bg) > 0 && iszero(last(bg))
+    while length(bf) > 0 && iszero(last(bf))
-        pop!(bg)
+        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
+    # first compute gcd to take out common factors
-    g = copy(list[ind[1]])
+    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]))