diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml new file mode 100644 index 0000000..c8ce604 --- /dev/null +++ b/.github/workflows/CompatHelper.yml @@ -0,0 +1,25 @@ +name: CompatHelper +on: + schedule: + - cron: 0 0 * * * + workflow_dispatch: +jobs: + CompatHelper: + runs-on: ubuntu-latest + steps: + - name: "Install CompatHelper" + run: | + import Pkg + name = "CompatHelper" + uuid = "aa819f21-2bde-4658-8897-bab36330d9b7" + version = "2" + Pkg.add(; name, uuid, version) + shell: julia --color=yes {0} + - name: "Run CompatHelper" + run: | + import CompatHelper + CompatHelper.main() + shell: julia --color=yes {0} + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + COMPATHELPER_PRIV: ${{ secrets.COMPATHELPER_PRIV }} diff --git a/.github/workflows/ci-julia-nightly.yml b/.github/workflows/ci-julia-nightly.yml new file mode 100644 index 0000000..37d6001 --- /dev/null +++ b/.github/workflows/ci-julia-nightly.yml @@ -0,0 +1,33 @@ +name: CI (Julia nightly) +on: + - push + - pull_request +jobs: + test: + name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + version: + - 'nightly' + os: + - ubuntu-latest + - macOS-latest + - windows-latest + arch: + - x64 + steps: + - uses: actions/checkout@v2 + - uses: julia-actions/setup-julia@v1 + with: + version: ${{ matrix.version }} + arch: ${{ matrix.arch }} + - uses: julia-actions/julia-buildpkg@latest + - uses: julia-actions/julia-runtest@latest + env: + JULIA_NUM_THREADS: 2 + - uses: julia-actions/julia-processcoverage@v1 + - uses: codecov/codecov-action@v1 + with: + file: lcov.info diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml new file mode 100644 index 0000000..d72335f --- /dev/null +++ b/.github/workflows/ci.yml @@ -0,0 +1,34 @@ +name: CI +on: + - push + - pull_request +jobs: + test: + name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + version: + - '1.0' + - '1' + os: + - ubuntu-latest + - macOS-latest + - windows-latest + arch: + - x64 + steps: + - uses: actions/checkout@v2 + - uses: julia-actions/setup-julia@v1 + with: + version: ${{ matrix.version }} + arch: ${{ matrix.arch }} + - uses: julia-actions/julia-buildpkg@latest + - uses: julia-actions/julia-runtest@latest + env: + JULIA_NUM_THREADS: 2 + - uses: julia-actions/julia-processcoverage@v1 + - uses: codecov/codecov-action@v1 + with: + file: lcov.info diff --git a/.travis.yml b/.travis.yml deleted file mode 100644 index c37bba4..0000000 --- a/.travis.yml +++ /dev/null @@ -1,34 +0,0 @@ -# Documentation: http://docs.travis-ci.com/user/languages/julia/ -language: julia - -os: - - linux - - osx -julia: - - 1.0 - - 1 - - nightly -# env: -# - JULIA_NUM_THREADS=1 -# - JULIA_NUM_THREADS=4 - -notifications: - email: false - -matrix: - allow_failures: - - julia: nightly - -codecov: true -coveralls: true - -# jobs: -# include: -# - stage: "Documentation" -# julia: 1.1 -# os: linux -# script: -# - julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); -# Pkg.instantiate()' -# - julia --project=docs/ docs/make.jl -# after_success: skip diff --git a/Project.toml b/Project.toml index 712e195..eafb1f7 100644 --- a/Project.toml +++ b/Project.toml @@ -1,17 +1,19 @@ name = "WignerSymbols" uuid = "9f57e263-0b3d-5e2e-b1be-24f2bb48858b" authors = ["Jutho Haegeman"] -version = "1.1.0" +version = "2.0.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" [compat] -RationalRoots = "0.1" -HalfIntegers = "1.0" -Primes = "0.4" +RationalRoots = "0.1 - 1" +HalfIntegers = "1" +Primes = "0.4 - 1" +LRUCache = "1.3" julia = "1" [extras] diff --git a/appveyor.yml b/appveyor.yml deleted file mode 100644 index 5150cd5..0000000 --- a/appveyor.yml +++ /dev/null @@ -1,43 +0,0 @@ -environment: - matrix: - - julia_version: 1.0 - - julia_version: 1 - - julia_version: nightly - -platform: - - x86 # 32-bit - - x64 # 64-bit - -# # Uncomment the following lines to allow failures on nightly julia -# # (tests will run but not make your overall status red) -matrix: - allow_failures: - - julia_version: nightly - -branches: - only: - - master - - /release-.*/ - -notifications: - - provider: Email - on_build_success: false - on_build_failure: false - on_build_status_changed: false - -install: - - ps: iex ((new-object net.webclient).DownloadString("https://raw.githubusercontent.com/JuliaCI/Appveyor.jl/version-1/bin/install.ps1")) - -build_script: - - echo "%JL_BUILD_SCRIPT%" - - C:\julia\bin\julia -e "%JL_BUILD_SCRIPT%" - -test_script: - - echo "%JL_TEST_SCRIPT%" - - C:\julia\bin\julia -e "%JL_TEST_SCRIPT%" - -# # Uncomment to support code coverage upload. Should only be enabled for packages -# # which would have coverage gaps without running on Windows -# on_success: -# - echo "%JL_CODECOV_SCRIPT%" -# - C:\julia\bin\julia -e "%JL_CODECOV_SCRIPT%" diff --git a/src/WignerSymbols.jl b/src/WignerSymbols.jl index 5aefe57..41279c8 100644 --- a/src/WignerSymbols.jl +++ b/src/WignerSymbols.jl @@ -2,16 +2,24 @@ __precompile__(true) module WignerSymbols export δ, Δ, clebschgordan, wigner3j, wigner6j, racahV, racahW, HalfInteger -using Base.GMP.MPZ using HalfIntegers using RationalRoots +using LRUCache const RRBig = RationalRoot{BigInt} import RationalRoots: _convert +include("growinglist.jl") include("primefactorization.jl") +convert(BigInt, primefactorial(401)) # trigger compilation and generate some fixed data -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}}}() +const Key3j = Tuple{UInt,UInt,UInt,Int,Int} +const Key6j = NTuple{6,UInt} + +# const Wigner3j = Dict{Key3j,Tuple{Rational{BigInt},Rational{BigInt}}}() +# const Wigner6j = Dict{Key6j,Tuple{Rational{BigInt},Rational{BigInt}}}() +# +const Wigner3j = LRU{Key3j,Tuple{Rational{BigInt},Rational{BigInt}}}(; maxsize = 10^6) +const Wigner6j = LRU{Key6j,Tuple{Rational{BigInt},Rational{BigInt}}}(; maxsize = 10^6) # check integerness and correctness of (j,m) angular momentum ϵ(j, m) = (abs(m) <= j && ishalfinteger(j) && isinteger(j-m) && isinteger(j+m)) @@ -44,7 +52,8 @@ function Δ(T::Type{<:Real}, j₁, j₂, j₃) return zero(T) end n, d = Δ²(j₁, j₂, j₃) - return convert(T, signedroot(RationalRoot{BigInt}, n//d)) + r = Base.unsafe_rational(n, d) + return convert(T, signedroot(RationalRoot{BigInt}, r)) end """ @@ -64,6 +73,11 @@ function wigner3j(T::Type{<:Real}, j₁, j₂, j₃, m₁, m₂, m₃ = -m₁-m for (jᵢ,mᵢ) in ((j₁, m₁), (j₂, m₂), (j₃, m₃)) ϵ(jᵢ, mᵢ) || throw(DomainError((jᵢ, mᵢ), "invalid combination (jᵢ, mᵢ)")) end + return _wigner3j(T, HalfInteger.((j₁, j₂, j₃, m₁, m₂, m₃))...) +end + +function _wigner3j(T::Type{<:Real}, j₁::HalfInteger, j₂::HalfInteger, j₃::HalfInteger, + m₁::HalfInteger, m₂::HalfInteger, m₃::HalfInteger) # check triangle condition and m₁+m₂+m₃ == 0 if !δ(j₁, j₂, j₃) || !iszero(m₁+m₂+m₃) return zero(T) @@ -74,9 +88,9 @@ function wigner3j(T::Type{<:Real}, j₁, j₂, j₃, m₁, m₂, m₃ = -m₁-m # TODO: 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₂ ) + β₁ = convert(UInt, j₁ + j₂ - j₃ ) + β₂ = convert(UInt, j₁ - m₁ ) + β₃ = convert(UInt, j₂ + m₂ ) # extra sign in definition: α₁ - α₂ = j₁ + m₂ - j₂ + m₁ = j₁ - j₂ + m₃ sgn = isodd(α₁ - α₂) ? -sgn : sgn @@ -90,8 +104,10 @@ function wigner3j(T::Type{<:Real}, j₁, j₂, j₃, m₁, m₂, m₃ = -m₁-m snum, rnum = splitsquare(s1n*s2n) sden, rden = splitsquare(s1d) - s = convert(BigInt, snum) // convert(BigInt, sden) - r = convert(BigInt, rnum) // convert(BigInt, rden) + 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 *= compute3jseries(β₁, β₂, β₃, α₁, α₂) Wigner3j[(β₁, β₂, β₃, α₁, α₂)] = (r,s) end @@ -151,7 +167,11 @@ function wigner6j(T::Type{<:Real}, j₁, j₂, j₃, j₄, j₅, j₆) for jᵢ in (j₁, j₂, j₃, j₄, j₅, j₆) (ishalfinteger(jᵢ) && jᵢ >= zero(jᵢ)) || throw(DomainError("invalid jᵢ", jᵢ)) end + return _wigner6j(T, HalfInteger.((j₁, j₂, j₃, j₄, j₅, j₆))...) +end +function _wigner6j(T::Type{<:Real}, j₁::HalfInteger, j₂::HalfInteger, j₃::HalfInteger, + j₄::HalfInteger, j₅::HalfInteger, j₆::HalfInteger) α̂₁ = (j₁, j₂, j₃) α̂₂ = (j₁, j₆, j₅) α̂₃ = (j₂, j₄, j₆) @@ -186,10 +206,10 @@ function wigner6j(T::Type{<:Real}, j₁, j₂, j₃, j₄, j₅, j₆) snum, rnum = splitsquare(n₁ * n₂ * n₃ * n₄) sden, rden = splitsquare(d₁ * d₂ * d₃ * d₄) - snu, sden = divgcd!(snum, sden) - rnu, rden = divgcd!(rnum, rden) - s = convert(BigInt, snum) // convert(BigInt, sden) - r = convert(BigInt, rnum) // convert(BigInt, rden) + 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 *= compute6jseries(β₁, β₂, β₃, α₁, α₂, α₃, α₄) Wigner6j[(β₁, β₂, β₃, α₁, α₂, α₃)] = (r, s) @@ -223,12 +243,13 @@ end # squared triangle coefficient function Δ²(j₁, j₂, j₃) # also checks the triangle conditions by converting to unsigned integer: - n1 = primefactorial( convert(UInt, + j₁ + j₂ - j₃) ) + n1 = copy(primefactorial( convert(UInt, + j₁ + j₂ - j₃) )) n2 = primefactorial( convert(UInt, + j₁ - j₂ + j₃) ) n3 = primefactorial( convert(UInt, - j₁ + j₂ + j₃) ) - d = primefactorial( convert(UInt, j₁ + j₂ + j₃ + 1) ) + num = mul!(mul!(n1, n2), n3) + den = copy(primefactorial( convert(UInt, j₁ + j₂ + j₃ + 1) )) # result - return (n1*n2*n3), d + return divgcd!(num, den) end # reorder parameters determining the 3j symbol to canonical order: @@ -278,14 +299,30 @@ function compute3jseries(β₁, β₂, β₃, α₁, α₂) dens = Vector{T}(undef, length(krange)) for (i, k) in enumerate(krange) num = iseven(k) ? one(T) : -one(T) - den = primefactorial(k)*primefactorial(k-α₁)*primefactorial(k-α₂)* - primefactorial(β₁-k)*primefactorial(β₂-k)*primefactorial(β₃-k) - nums[i], dens[i] = divgcd!(num, den) + den = copy(primefactorial(k)) + den = mul!(mul!(den, primefactorial(k-α₁)), primefactorial(k-α₂)) + den = mul!(mul!(mul!(den, primefactorial(β₁-k)), + primefactorial(β₂-k)), + primefactorial(β₃-k)) + nums[i], dens[i] = num, den end den = commondenominator!(nums, dens) totalnum = sumlist!(nums) totalden = convert(BigInt, den) - return totalnum//totalden + for n = 1:length(den.powers) + p = bigprime(n) + while den.powers[n] > 0 + q, r = divrem(totalnum, p) + if iszero(r) + totalnum = q + den.powers[n] -= 1 + else + break + end + end + end + totalden = convert(BigInt, den) + return Base.unsafe_rational(totalnum, totalden) end # compute the sum appearing in the 6j symbol @@ -296,15 +333,32 @@ function compute6jseries(β₁, β₂, β₃, α₁, α₂, α₃, α₄) nums = Vector{T}(undef, length(krange)) dens = Vector{T}(undef, 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) + num = iseven(k) ? copy(primefactorial(k+1)) : neg!(copy(primefactorial(k+1))) + den = copy(primefactorial(k-α₁)) + den = mul!(mul!(mul!(den, primefactorial(k-α₂)), + primefactorial(k-α₃)), + primefactorial(k-α₄)) + den = mul!(mul!(mul!(den, primefactorial(β₁-k)), + primefactorial(β₂-k)), + primefactorial(β₃-k)) nums[i], dens[i] = divgcd!(num, den) end den = commondenominator!(nums, dens) totalnum = sumlist!(nums) + for n = 1:length(den.powers) + p = bigprime(n) + while den.powers[n] > 0 + q, r = divrem(totalnum, p) + if iszero(r) + totalnum = q + den.powers[n] -= 1 + else + break + end + end + end totalden = convert(BigInt, den) - return totalnum//totalden + return Base.unsafe_rational(totalnum, totalden) end end # module diff --git a/src/growinglist.jl b/src/growinglist.jl new file mode 100644 index 0000000..52a2c6a --- /dev/null +++ b/src/growinglist.jl @@ -0,0 +1,154 @@ +using Base.Threads: Atomic, SpinLock + +# ListSegment represents a segment from a GrowingList; it has a list `data` to hold the elements, filled up to `currentlength`, and possibly a reference to the next segment, if it is not the final segment. +mutable struct ListSegment{T} + data::Vector{T} + currentlength::Int + next::Base.RefValue{ListSegment{T}} +end +ListSegment{T}(data::Vector{T}, currentlength::Int) where T = + ListSegment{T}(data, currentlength, Ref{ListSegment{T}}()) + +# getindex, assumes that index is inbounds, traverses the linked list +function _unsafe_getindex(l::ListSegment, i::Int) + if i <= l.currentlength + getindex(l.data, i) + else + _unsafe_getindex(l.next[], i - l.currentlength) + end +end + +# get or push a new element at the end; in itself not thread safe, should be protected by the lock in the parent GrowingList +function _unsafe_get!(l::ListSegment{T}, n::Int, default, newlength) where T + N = length(l.data) + if n > N + if isassigned(l.next) + return _unsafe_get!(l.next[], n - N, default, newlength) + else + newsegment = Vector{T}(undef, newlength) + newsegment[1] = default + l.next = Ref(ListSegment{T}(newsegment, 1)) + return default + end + else + if n <= l.currentlength + @inbounds return getindex(l.data, n) + else + @assert n == l.currentlength+1 + l.data[n] = default + l.currentlength += 1 + return default + end + end +end + +""" + GrowingList{T} <: AbstractVector{T} + + GrowingList([iter,]; sizehint = max(16, length(iter)), growthfactor = 2.) + +A thread safe vector / list data structure where new elements can be added at the back. +Once an element is set, it cannot be changed or removed. This ensures thread safe +`getindex` of that element without requiring a lock. The `length` of a `GrowingList` +instance can also be probed without a lock, but the return value will be a lower bound, +i.e. the list can already have increased in length at the same time. + +New elements can be added using the syntax +`get!(l::GrowingList, i::Int, value)` +`get!(value_generator::Callable, l::GrowingList, i::Int)` + +where the new element `value` or `value_generator()` will only be added if `i` is +`length(l)+1`. If multiple tasks or threads try to `get!` the same index `i`, only one of +them will actually be adding that element. The `value` or `value_generator()` produced by +the different threads should be the same to avoid unpredictable results. + +The list is grown by adding new segments using a linked list data structure. This guarantees that existing data does never have to move in memory, which is required in order to make `getindex` threadsafe without lock. +""" +mutable struct GrowingList{T} <: AbstractVector{T} + first::ListSegment{T} + totallength::Atomic{Int} + growthfactor::Float64 + lock::SpinLock + function GrowingList{T}(iter; + sizehint = max(16, length(iter)), + growthfactor = 2.) where {T} + firstsegment = Vector{T}(undef, sizehint) + i = 0 + next = iterate(iter) + @inbounds while i < sizehint && next !== nothing + i += 1 + val, state = next + firstsegment[i] = val + next = iterate(iter, state) + end + first = ListSegment{T}(firstsegment, i) + while next !== nothing + i += 1 + val, state = next + _unsafe_getindex(first, i, val, ceil(Int, (i-1)*growthfactor)) + next = iterate(iter, state) + end + return new{T}(first, Atomic{Int}(i), growthfactor, SpinLock()) + end +end +GrowingList(v::Vector{T}; sizehint = max(16, length(v)), growthfactor = 2.) where {T} = + GrowingList{T}(v; sizehint = sizehint, growthfactor = growthfactor) + +GrowingList{T}(; sizehint = 16, growthfactor = 2.) where {T} = + GrowingList{T}((); sizehint = sizehint, growthfactor = growthfactor) + +GrowingList(; sizehint = 16, growthfactor = 2.) = + GrowingList{Any}((); sizehint = sizehint, growthfactor = growthfactor) + +Base.length(l::GrowingList) = l.totallength[] +Base.size(l::GrowingList) = (length(l),) + +@inline function Base.getindex(l::GrowingList, n::Int) + @boundscheck checkbounds(l, n) + return _unsafe_getindex(l.first, n) +end + +function Base.get!(l::GrowingList, n::Int, default) + if n <= l.totallength[] + return _unsafe_getindex(l.first, n) + else + lock(l.lock) + len = length(l) + if n <= len # try again, maybe already ok now + unlock(l.lock) + return _unsafe_getindex(l.first, n) + elseif n == len+1 + _unsafe_get!(l.first, n, default, ceil(Int, (l.growthfactor-1)*len)) + Base.Threads.atomic_add!(l.totallength, 1) + unlock(l.lock) + return default + else + @show Base.Threads.threadid(), l.totallength[], n + unlock(l.lock) + throw(ArgumentError("can only insert new element at next index: $(len+1)")) + end + end +end + +function Base.get!(default::Base.Callable, l::GrowingList, n::Int) + if n <= l.totallength[] + return _unsafe_getindex(l.first, n) + else + v = default() + lock(l.lock) + len = l.totallength[] + if n <= len # try again, maybe already ok now + unlock(l.lock) + return _unsafe_getindex(l.first, n) + elseif n == len+1 + _unsafe_get!(l.first, n, v, ceil(Int, (l.growthfactor-1)*len)) + Base.Threads.atomic_add!(l.totallength, 1) + unlock(l.lock) + return v + else + @show Base.Threads.threadid(), l.totallength[], n + unlock(l.lock) + throw(ArgumentError("can only insert new element at next index: $(len+1)")) + end + end +end diff --git a/src/primefactorization.jl b/src/primefactorization.jl index a173e8f..9a7b00a 100644 --- a/src/primefactorization.jl +++ b/src/primefactorization.jl @@ -1,15 +1,15 @@ using Primes: isprime import Base.divgcd -const primetable = - [2,3,5] -const factortable = - [UInt8[], UInt8[1], UInt8[0,1], UInt8[2], UInt8[0,0,1]] -const factorialtable = - [UInt32[], UInt32[], UInt32[1], UInt32[1,1], UInt32[3,1], UInt32[3,1,1]] -const bigprimetable = - [[big(2)], [big(3)], [big(5)]] -const bigone = Ref{BigInt}(big(1)) +using Base.GMP.MPZ + +const primetable = GrowingList([2, 3]; sizehint = 256) +const factortable = GrowingList([UInt8[], UInt8[1], UInt8[0,1]]; sizehint = 1024) +const factorialtable = GrowingList([UInt32[], UInt32[1], UInt32[1,1]]; sizehint = 1024) +const bigprimetable = GrowingList([GrowingList([big(2)]; sizehint = 512), + GrowingList([big(3)]; sizehint = 256)]; + sizehint = 256) +const bigone = big(1) # Make a prime iterator struct PrimeIterator @@ -22,51 +22,71 @@ Base.eltype(::PrimeIterator) = Int # Get the `n`th prime; store all primes up to the `n`th if not yet available function prime(n::Int) - p = last(primetable) - while length(primetable) < n + k = min(length(primetable), length(bigprimetable)) + p = primetable[k] + while k < n p = p + 2 while !isprime(p) p += 2 end - push!(primetable, p) - push!(bigprimetable, [big(p)]) + k += 1 + # these lines do not get but set new elements; provided no other task did so earlier + get!(primetable, k, p) + get!(bigprimetable, k, GrowingList([big(p)]; sizehint = 256)) + k = min(length(primetable), length(bigprimetable)) + # other threads might have inserted additional entries, + # make sure they are finished with both primetable and bigprimetable end - @inbounds return primetable[n] + return primetable[n] end Base.iterate(::PrimeIterator, n = 1) = prime(n), n+1 # get primes and their powers as `BigInt`, also cache all results function bigprime(n::Integer, e::Integer=1) - e == 0 && return bigone[] + e == 0 && return bigone p = prime(n) # triggers computation of prime(n) if necessary - @inbounds l = length(bigprimetable[n]) + powerlist = bigprimetable[n] + l = length(powerlist) @inbounds while l < e # compute next prime power as approximate square of existing results - k = (l+1)>>1 - push!(bigprimetable[n], bigprimetable[n][k]*bigprimetable[n][l+1-k]) l += 1 + k = l>>1 + get!(powerlist, l, powerlist[k]*powerlist[l-k]) + l = length(powerlist) # other threads might have inserted more powers end - @inbounds return bigprimetable[n][e] + @inbounds return powerlist[e] end # A custom `Integer` subtype to store an integer as its prime factorization -struct PrimeFactorization{U<:Unsigned} <: Integer +# mutable to allow in place update of sign +mutable struct PrimeFactorization{U<:Unsigned} <: Integer powers::Vector{U} sign::Int8 + PrimeFactorization{U}(powers::Vector, sign = one(Int8)) where {U<:Unsigned} = + new{U}(convert(Vector{U}, powers), sign) +end +# convenience constructor: normalizes powers to have last entry nonzero +PrimeFactorization(powers::Vector{U}, sign = one(Int8)) where {U<:Unsigned} = + PrimeFactorization{U}(_normalize_powers!(powers), sign) + +function _normalize_powers!(v::Vector{<:Integer}) + i = findlast(!iszero, v) + l = ifelse(i === nothing, 0, i) + l < length(v) && resize!(v, l) + return v end -PrimeFactorization(powers::Vector{U}) where {U<:Unsigned} = - PrimeFactorization{U}(powers, one(Int8)) # define our own factor function, returning an instance of PrimeFactorization function primefactor(n::Integer) - iszero(n) && return PrimeFactorization(UInt8[], zero(Int8)) + iszero(n) && return PrimeFactorization{UInt8}(UInt8[], zero(Int8)) sn = n < 0 ? -one(Int8) : one(Int8) n = abs(n) m = length(factortable) while m < abs(n) m += 1 - powers = UInt8[] # should be sufficient for all integers up to 2^255 + powers = UInt8[] + # should be sufficient for all integers up to 2^255 a = m for p in primes() f = 0 @@ -79,16 +99,18 @@ function primefactor(n::Integer) push!(powers, f) a == 1 && break end - push!(factortable, powers) + get!(factortable, m, powers) + m = length(factortable) # other threads may have inserted other entries end - @inbounds return PrimeFactorization(copy(factortable[n]), sn) + @inbounds return PrimeFactorization{UInt8}(factortable[n], sn) end function primefactorial(n::Integer) - n < 0 && throw(DomainError(n)) - m = length(factorialtable)-1 + n < 0 && throw(DomainError(n,"primefactorial only works for non-negative numbers")) + n <= 1 && return PrimeFactorization{UInt32}(UInt32[], one(Int8)) + m = length(factorialtable) @inbounds while m < n - prevfactorial = factorialtable[m+1] + prevfactorial = factorialtable[m] m += 1 f = primefactor(m).powers powers = copy(prevfactorial) @@ -98,36 +120,48 @@ function primefactorial(n::Integer) for k = 1:length(f) powers[k] += f[k] end - push!(factorialtable, powers) + get!(factorialtable, m, powers) + m = length(factorialtable) # other threads may have inserted other entries end - @inbounds return PrimeFactorization(copy(factorialtable[n+1])) + @inbounds return PrimeFactorization{UInt32}(factorialtable[n]) end # Methods for PrimeFactorization: Base.copy(a::PrimeFactorization) = PrimeFactorization(copy(a.powers), a.sign) +function Base.copy!(c::PrimeFactorization, a::PrimeFactorization) + c.sign = a.sign + copy!(c.powers, a.powers) + return c +end Base.one(::Type{PrimeFactorization{U}}) where {U<:Unsigned} = - PrimeFactorization(Vector{U}()) + PrimeFactorization{U}(Vector{U}(), one(Int8)) Base.zero(::Type{PrimeFactorization{U}}) where {U<:Unsigned} = - PrimeFactorization(Vector{U}(), zero(Int8)) + PrimeFactorization{U}(Vector{U}(), zero(Int8)) + +one!(c::PrimeFactorization) = (c.sign = one(Int8); empty!(c.powers); return c) +zero!(c::PrimeFactorization) = (c.sign = zero(Int8); empty!(c.powers); return c) Base.promote_rule(P::Type{<:PrimeFactorization},::Type{<:Integer}) = P -Base.promote_rule(P::Type{<:PrimeFactorization},::Type{BigInt}) = BigInt +Base.promote_rule(::Type{<:PrimeFactorization},::Type{BigInt}) = BigInt +Base.promote_rule(::Type{PrimeFactorization{U1}}, + ::Type{PrimeFactorization{U2}}) where {U1<:Unsigned, U2<:Unsigned} = PrimeFactorization{promote_type(U1, U2)} Base.convert(P::Type{<:PrimeFactorization}, n::Integer) = convert(P, primefactor(n)) -function Base.convert(::Type{BigInt}, a::PrimeFactorization) - A = one(BigInt) +function _convert!(x::BigInt, a::PrimeFactorization) + MPZ.set!(x, bigone) for (n, e) in enumerate(a.powers) if !iszero(e) - MPZ.mul!(A, bigprime(n, e)) + MPZ.mul!(x, bigprime(n, e)) end end - return a.sign < 0 ? MPZ.neg!(A) : A + return a.sign < 0 ? MPZ.neg!(x) : x end +Base.convert(::Type{BigInt}, a::PrimeFactorization) = _convert!(one(BigInt), a) Base.convert(::Type{PrimeFactorization{U}}, a::PrimeFactorization{U}) where {U<:Unsigned} = a Base.convert(::Type{PrimeFactorization{U}}, a::PrimeFactorization) where {U<:Unsigned} = - PrimeFactorization(convert(Vector{U}, a.powers), a.sign) + PrimeFactorization{U}(convert(Vector{U}, a.powers), a.sign) Base.:(==)(a::PrimeFactorization, b::PrimeFactorization) = a.powers == b.powers && a.sign == b.sign @@ -138,7 +172,8 @@ function Base.:<(a::PrimeFactorization, b::PrimeFactorization) return <(-b, -a) else ag, bg = divgcd(a, b) - if length(ag.powers) <= length(bg.powers) && + ag == bg && return false + if length(ag.powers) <= length(bg.powers) all(k->ag.powers[k] la + throw(DivideError()) + end + lc = la + if lb == lc + while lc > 0 && a.powers[lc] == b.powers[lc] + lc -= 1 + end + end + lc == length(c.powers) || resize!(c.powers, lc) + @inbounds for k = 1:min(lb, lc) + if b.powers[k] > a.powers[k] + throw(DivideError()) + end + c.powers[k] = a.powers[k] - b.powers[k] + end + if c !== a + @inbounds for k = lb+1:lc + c.powers[k] = a.powers[k] + end + end end + return c end -function Base.lcm(a::PrimeFactorization{T}, b::PrimeFactorization{T}) where {T} +function gcd!(c::PrimeFactorization, a::PrimeFactorization, b::PrimeFactorization) if a.sign == 0 - return a + copy!(c.powers, b.powers) elseif b.sign ==0 - return b + copy!(c.powers, a.powers) else - return PrimeFactorization(_vmax!(copy(a.powers), b.powers)) + c.sign = one(Int8) + la = length(a.powers) + lb = length(b.powers) + lc = min(la, lb) + lc === length(c.powers) || resize!(c.powers, lc) + @inbounds for k = 1:lc + c.powers[k] = min(a.powers[k], b.powers[k]) + end end + c.sign = one(Int8) + return c +end +function lcm!(c::PrimeFactorization, a::PrimeFactorization, b::PrimeFactorization) + if a.sign == 0 || b.sign == 0 + return zero!(c) + else + c.sign = one(Int8) + la = length(a.powers) + lb = length(b.powers) + lc = max(la, lb) + lc === length(c.powers) || resize!(c.powers, lc) + @inbounds for k = 1:min(la,lb) + c.powers[k] = max(a.powers[k], b.powers[k]) + end + if c !== a + @inbounds for k = lb+1:la + c.powers[k] = a.powers[k] + end + end + @inbounds for k = la+1:lb + c.powers[k] = b.powers[k] + end + end + c.sign = one(Int8) + return c end -Base.divgcd(a::PrimeFactorization, b::PrimeFactorization) = divgcd!(copy(a), copy(b)) function divgcd!(a::PrimeFactorization, b::PrimeFactorization) af, bf = a.powers, b.powers for k = 1:min(length(af), length(bf)) @@ -187,26 +294,50 @@ function divgcd!(a::PrimeFactorization, b::PrimeFactorization) af[k] -= gk bf[k] -= gk end - while length(af) > 0 && iszero(last(af)) - pop!(af) - end - while length(bf) > 0 && iszero(last(bf)) - pop!(bf) - end + _normalize_powers!(a.powers) + _normalize_powers!(b.powers) return a, b end +mul!(a::PrimeFactorization, b::PrimeFactorization) = mul!(a, a, b) +divexact!(a::PrimeFactorization, b::PrimeFactorization) = divexact!(a, a, b) +gcd!(a::PrimeFactorization, b::PrimeFactorization) = gcd!(a, a, b) +lcm!(a::PrimeFactorization, b::PrimeFactorization) = lcm!(a, a, b) + +Base.:-(a::PrimeFactorization) = neg!(copy(a)) +function Base.:*(a::PrimeFactorization, b::PrimeFactorization) + P = promote_type(typeof(a), typeof(b)) + if length(a.powers) >= length(b.powers) + return typeof(a) == P ? mul!(copy(a), b) : mul!(convert(P, a), b) + else + return typeof(b) == P ? mul!(copy(b), a) : mul!(convert(P, b), a) + end +end +function Base.lcm(a::PrimeFactorization, b::PrimeFactorization) + P = promote_type(typeof(a), typeof(b)) + if length(a.powers) >= length(b.powers) + return typeof(a) == P ? lcm!(copy(a), b) : lcm!(convert(P, a), b) + else + return typeof(b) == P ? lcm!(copy(b), a) : lcm!(convert(P, b), a) + end +end +function Base.gcd(a::PrimeFactorization, b::PrimeFactorization) + P = promote_type(typeof(a), typeof(b)) + if length(a.powers) <= length(b.powers) + return typeof(a) == P ? lcm!(copy(a), b) : lcm!(convert(P, a), b) + else + return typeof(b) == P ? lcm!(copy(b), a) : lcm!(convert(P, b), a) + end +end +Base.divgcd(a::PrimeFactorization, b::PrimeFactorization) = divgcd!(copy(a), copy(b)) +# no promotion necessary, should be smaller than a +divexact(a::PrimeFactorization, b::PrimeFactorization) = divexact!(copy(a), b) + # 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 @@ -215,13 +346,13 @@ end function commondenominator!(nums::Vector{P}, dens::Vector{P}) where {P<:PrimeFactorization} isempty(nums) && return one(P) # accumulate lcm of denominator - den = PrimeFactorization(copy(dens[1].powers)) + den = copy(dens[1]) for i = 2:length(dens) - _vmax!(den.powers, dens[i].powers) + lcm!(den, dens[i]) end # rescale numerators for i = 1:length(nums) - _vsub!(_vadd!(nums[i].powers, den.powers), dens[i].powers) + divexact!(mul!(nums[i], den), dens[i]) end return den end @@ -229,68 +360,17 @@ end # auxiliary function to compute sums of a list of PrimeFactorizations as quickly as possible function sumlist!(list::Vector{<:PrimeFactorization}, ind = 1:length(list)) # first compute gcd to take out common factors - g = PrimeFactorization(copy(list[ind[1]].powers)) + g = copy(list[ind[1]]) for k in ind - _vmin!(g.powers, list[k].powers) + gcd!(g, list[k]) end for k in ind - _vsub!(list[k].powers, g.powers) + divexact!(list[k], g) end - L = length(ind) - if L > 32 - l = L >> 1 - s = sumlist!(list, first(ind).+(0:l-1)) + sumlist!(list, first(ind).+(l:L-1)) - else - # do sum - s = big(0) - for k in ind - MPZ.add!(s, convert(BigInt, list[k])) - end + s = big(0) + i = big(1) + for p in list + MPZ.add!(s, _convert!(i, p)) end - return MPZ.mul!(s, convert(BigInt, g)) -end - -# Mutating vector methods that also grow and shrink as required -function _vmin!(af::Vector{U}, bf::Vector{U}) where {U<:Unsigned} - while length(af) > length(bf) - pop!(af) - end - @inbounds for k = 1:length(af) - af[k] = min(af[k], bf[k]) - end - while length(af) > 0 && iszero(last(af)) - pop!(af) - end - return af -end -function _vmax!(af::Vector{U}, bf::Vector{U}) where {U<:Unsigned} - while length(bf) > length(af) - push!(af, zero(U)) - end - @inbounds for k = 1:length(bf) - af[k] = max(af[k], bf[k]) - end - return af -end -function _vadd!(af::Vector{U}, bf::Vector{U}) where {U<:Unsigned} - while length(bf) > length(af) - push!(af, zero(U)) - end - @inbounds for k = 1:length(bf) - af[k] = +(af[k], bf[k]) - end - return af -end -function _vsub!(af::Vector{U}, bf::Vector{U}) where {U<:Unsigned} - if length(bf) > length(af) - throw(OverflowError()) - end - @inbounds for k = 1:length(bf) - bf[k] > af[k] && throw(OverflowError()) - af[k] -= bf[k] - end - while length(af) > 0 && iszero(last(af)) - pop!(af) - end - return af + return MPZ.mul!(s, _convert!(i, g)) end diff --git a/test/runtests.jl b/test/runtests.jl index 2cad8ab..af787d1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,165 +2,203 @@ using Test using WignerSymbols using LinearAlgebra using Random - +using Base.Threads +N = Base.Threads.nthreads() Random.seed!(1234) smalljlist = 0:1//2:10 largejlist = 0:1//2:1000 -@testset "triangle coefficient" begin - for j1 in smalljlist, j2 in smalljlist - for j3 = abs(j1-j2):(j1+j2) - @test Δ(j1,j2,j3) ≈ sqrt(factorial(big(Int(j1+j2-j3)))* - factorial(big(Int(j1-j2+j3)))*factorial(big(Int(j2+j3-j1)))/ - factorial(big(Int(j1+j2+j3+1)))) +@threads for i = 1:N + @testset "triangle coefficient, thread $i" begin + for k = i:N:length(smalljlist) + j1 = smalljlist[k] + for j2 in smalljlist + for j3 = abs(j1-j2):(j1+j2) + @test Δ(j1,j2,j3) ≈ sqrt(factorial(big(Int(j1+j2-j3)))* + factorial(big(Int(j1-j2+j3)))* + factorial(big(Int(j2+j3-j1)))/ + factorial(big(Int(j1+j2+j3+1)))) + end + end end end end # test 3j: #-------- -@testset "clebschgordan: test orthogonality" begin - for j1 in smalljlist, j2 in smalljlist - d1::Int = 2*j1+1 - d2::Int = 2*j2+1 - M = zeros(Float64, (d1*d2, d1*d2)) - ind1 = 1 - for m1 in -j1:j1, m2 in -j2:j2 - ind2 = 1 - @inbounds for j3 in abs(j1-j2):(j1+j2), m3 in -j3:j3 - M[ind1,ind2] = clebschgordan(j1,m1,j2,m2,j3,m3) - ind2 += 1 +@threads for i = 1:N + @testset "clebschgordan: orthogonality relations, thread $i" begin + for k = i:N:length(smalljlist) + j1 = smalljlist[k] + for j2 in smalljlist + d1::Int = 2*j1+1 + d2::Int = 2*j2+1 + M = zeros(Float64, (d1*d2, d1*d2)) + ind1 = 1 + for m1 in -j1:j1, m2 in -j2:j2 + ind2 = 1 + @inbounds for j3 in abs(j1-j2):(j1+j2), m3 in -j3:j3 + M[ind1,ind2] = clebschgordan(j1,m1,j2,m2,j3,m3) + ind2 += 1 + end + ind1 += 1 + end + @test M'*M ≈ one(M) end - ind1 += 1 end - @test M'*M ≈ one(M) end end # test recurrence relations: Phys Rev E 57, 7274 (1998) -@testset "wigner3j: test recurrence relations" begin - for k = 1:10 - j2 = convert(BigFloat, rand(0:1//2:1000)) - j3 = convert(BigFloat, rand(0:1//2:1000)) - m2 = convert(BigFloat, rand(-j2:j2)) - m3 = convert(BigFloat, rand(-j3:j3)) +@threads for i = 1:N + @testset "wigner3j: recurrence relations, thread $i" begin + for k = 1:div(8,N) + j2 = convert(BigFloat, rand(largejlist)) + j3 = convert(BigFloat, rand(largejlist)) + m2 = convert(BigFloat, rand(-j2:j2)) + m3 = convert(BigFloat, rand(-j3:j3)) - for j in max(abs(j2-j3),abs(m2+m3))+1:(j2+j3)-1 - X = j*sqrt(((j+1)^2-(j2-j3)^2)*((j2+j3+1)^2-(j+1)^2)*((j+1)^2-(m2+m3)^2)) - Y = (2*j+1)*((m2+m3)*(j2*(j2+1)-j3*(j3+1)) - (m2-m3)*j*(j+1)) - Z = (j+1)*sqrt((j^2-(j2-j3)^2)*((j2+j3+1)^2-j^2)*(j^2-(m2+m3)^2)) - tol = 10*max(abs(X),abs(Y),abs(Z))*eps(BigFloat) - @test (X*wigner3j(BigFloat,j+1,j2,j3,-m2-m3,m2,m3) + - Z*wigner3j(BigFloat,j-1,j2,j3,-m2-m3,m2,m3)) ≈ - (-Y*wigner3j(BigFloat,j,j2,j3,-m2-m3,m2,m3)) atol=tol + for j in max(abs(j2-j3),abs(m2+m3))+1:(j2+j3)-1 + X = j*sqrt(((j+1)^2-(j2-j3)^2)*((j2+j3+1)^2-(j+1)^2)*((j+1)^2-(m2+m3)^2)) + Y = (2*j+1)*((m2+m3)*(j2*(j2+1)-j3*(j3+1)) - (m2-m3)*j*(j+1)) + Z = (j+1)*sqrt((j^2-(j2-j3)^2)*((j2+j3+1)^2-j^2)*(j^2-(m2+m3)^2)) + tol = 10*max(abs(X),abs(Y),abs(Z))*eps(BigFloat) + @test (X*wigner3j(BigFloat,j+1,j2,j3,-m2-m3,m2,m3) + + Z*wigner3j(BigFloat,j-1,j2,j3,-m2-m3,m2,m3)) ≈ + (-Y*wigner3j(BigFloat,j,j2,j3,-m2-m3,m2,m3)) atol=tol + end end end end -@testset "wigner3j: test orthogonality relations" begin - # equivalent to Clebsch-Gordan orthogonality, now test using Float32 - for j1 in smalljlist, j2 in smalljlist - d1::Int = 2*j1+1 - d2::Int = 2*j2+1 - M = zeros(Float32, (d1*d2, d1*d2)) - ind2 = 1 - for m1 in -j1:j1, m2 in -j2:j2 - ind1 = 1 - @inbounds for j3 in abs(j1-j2):(j1+j2), m3 in -j3:j3 - d3::Int = 2*j3+1 - M[ind1,ind2] += sqrt(d3) * wigner3j(Float32, j1, j2, j3, m1, m2, m3) - ind1 += 1 +@threads for i = 1:N + @testset "wigner3j: orthogonality relations, thread $i" begin + # equivalent to Clebsch-Gordan orthogonality, now test using Float32 + for k = i:N:length(smalljlist) + j1 = smalljlist[k] + for j2 in smalljlist + d1::Int = 2*j1+1 + d2::Int = 2*j2+1 + M = zeros(Float32, (d1*d2, d1*d2)) + ind2 = 1 + for m1 in -j1:j1, m2 in -j2:j2 + ind1 = 1 + @inbounds for j3 in abs(j1-j2):(j1+j2), m3 in -j3:j3 + d3::Int = 2*j3+1 + M[ind1,ind2] += sqrt(d3) * wigner3j(Float32, j1, j2, j3, m1, m2, m3) + ind1 += 1 + end + ind2 += 1 + end + @test M'*M ≈ one(M) # orthogonality relation type 1 + @test M*M' ≈ one(M) # orthogonality relation type 2 end - ind2 += 1 end - @test M'*M ≈ one(M) # orthogonality relation type 1 - @test M*M' ≈ one(M) # orthogonality relation type 2 end end # test 6j #---------- -@testset "wigner6j: test orthogonality" begin - for j1 in smalljlist, j2 in smalljlist, j4 in smalljlist - for j5 in max(abs(j1-j2-j4),abs(j1-j2+j4),abs(j1+j2-j4)):(j1+j2+j4) - j6range = max(abs(j2-j4),abs(j1-j5)):min((j2+j4),(j1+j5)) - j3range = max(abs(j1-j2),abs(j4-j5)):min((j1+j2),(j4+j5)) - @test length(j6range) == length(j3range) - M = zeros(Float64, (length(j3range), length(j6range))) - for (k2,j6) in enumerate(j6range) - for (k1,j3) in enumerate(j3range) - M[k1,k2] = sqrt(2*j3+1)*sqrt(2*j6+1)*wigner6j(j1,j2,j3,j4,j5,j6) +@threads for i = 1:N + @testset "wigner6j: orthogonality relations, thread $i" begin + for k = i:N:length(smalljlist) + j1 = smalljlist[k] + for j2 in smalljlist, j4 in smalljlist + for j5 in max(abs(j1-j2-j4),abs(j1-j2+j4),abs(j1+j2-j4)):(j1+j2+j4) + j6range = max(abs(j2-j4),abs(j1-j5)):min((j2+j4),(j1+j5)) + j3range = max(abs(j1-j2),abs(j4-j5)):min((j1+j2),(j4+j5)) + @test length(j6range) == length(j3range) + M = zeros(Float64, (length(j3range), length(j6range))) + for (k2,j6) in enumerate(j6range) + for (k1,j3) in enumerate(j3range) + M[k1,k2] = sqrt(2*j3+1)*sqrt(2*j6+1)*wigner6j(j1,j2,j3,j4,j5,j6) + end + end + @test M'*M ≈ one(M) end end - @test M'*M ≈ one(M) end end end -@testset "wigner6j: test special cases" begin - for j1 in smalljlist, j2 in smalljlist - j6 = 0 - j4 = j2 - j5 = j1 - for j3 in abs(j1-j2):(j1+j2) - @test wigner6j(j1,j2,j3,j4,j5,j6) ≈ (-1)^(j1+j2+j3)/sqrt((2*j1+1)*(2*j2+1)) +@threads for i = 1:N + @testset "wigner6j: special cases, thread $i" begin + for k = i:N:length(smalljlist) + j1 = smalljlist[k] + for j2 in smalljlist + j6 = 0 + j4 = j2 + j5 = j1 + for j3 in abs(j1-j2):(j1+j2) + @test wigner6j(j1,j2,j3,j4,j5,j6) ≈ + (-1)^(j1+j2+j3)/sqrt((2*j1+1)*(2*j2+1)) + end + end end end end -@testset "wigner6j: test recurrence relation" begin - for k = 1:10 - j2 = convert(BigFloat,rand(largejlist)) - j3 = convert(BigFloat,rand(largejlist)) - l1 = convert(BigFloat,rand(largejlist)) - l2 = convert(BigFloat,rand(abs(l1-j3):(l1+j3))) - l3 = convert(BigFloat,rand(abs(l1-j2):min(l1+j2))) +@threads for i = 1:N + @testset "wigner6j: recurrence relation, thread $i" begin + for k = 1:div(8,N) + j2 = convert(BigFloat,rand(largejlist)) + j3 = convert(BigFloat,rand(largejlist)) + l1 = convert(BigFloat,rand(largejlist)) + l2 = convert(BigFloat,rand(abs(l1-j3):(l1+j3))) + l3 = convert(BigFloat,rand(abs(l1-j2):min(l1+j2))) - for j in intersect(abs(j2-j3):(j2+j3), abs(l2-l3):(l2+l3)) - X = j * sqrt( ((j+1)^2-(j2-j3)^2) * ((j2+j3+1)^2-(j+1)^2) * - ((j+1)^2-(l2-l3)^2) * ((l2+l3+1)^2 - (j+1)^2) ) - Y = (2*j+1) * ( j*(j+1)*( -j*(j+1) + j2*(j2+1) + j3*(j3+1) - 2*l1*(l1+1)) + - l2*(l2+1)*( j*(j+1) + j2*(j2+1) - j3*(j3+1) ) + - l3*(l3+1)*( j*(j+1) - j2*(j2+1) + j3*(j3+1) ) ) - Z = (j+1) * sqrt( (j^2-(j2-j3)^2) * ((j2+j3+1)^2-j^2) * - (j^2-(l2-l3)^2) * ((l2+l3+1)^2-j^2) ) - tol = 10 * max(abs(X), abs(Y), abs(Z)) * eps(BigFloat) - @test (X*wigner6j(BigFloat, j+1, j2, j3, l1, l2, l3) + - Z*wigner6j(BigFloat, j-1, j2, j3, l1, l2, l3)) ≈ - (-Y*wigner6j(BigFloat, j, j2, j3, l1, l2, l3)) atol=tol + for j in intersect(abs(j2-j3):(j2+j3), abs(l2-l3):(l2+l3)) + X = j * sqrt( ((j+1)^2-(j2-j3)^2) * ((j2+j3+1)^2-(j+1)^2) * + ((j+1)^2-(l2-l3)^2) * ((l2+l3+1)^2 - (j+1)^2) ) + Y = (2*j+1) * ( j*(j+1)*( -j*(j+1) + j2*(j2+1) + j3*(j3+1) - 2*l1*(l1+1)) + + l2*(l2+1)*( j*(j+1) + j2*(j2+1) - j3*(j3+1) ) + + l3*(l3+1)*( j*(j+1) - j2*(j2+1) + j3*(j3+1) ) ) + Z = (j+1) * sqrt( (j^2-(j2-j3)^2) * ((j2+j3+1)^2-j^2) * + (j^2-(l2-l3)^2) * ((l2+l3+1)^2-j^2) ) + tol = 10 * max(abs(X), abs(Y), abs(Z)) * eps(BigFloat) + @test (X*wigner6j(BigFloat, j+1, j2, j3, l1, l2, l3) + + Z*wigner6j(BigFloat, j-1, j2, j3, l1, l2, l3)) ≈ + (-Y*wigner6j(BigFloat, j, j2, j3, l1, l2, l3)) atol=tol + end end end end -@testset "test recoupling relation between 3j/clebschgordan and 6j/racahW symbols" begin - smallerjlist = 0:1//2:5 - for j1 in smallerjlist, j2 in smallerjlist, j3 in smallerjlist - m1range = -j1:j1 - m2range = -j2:j2 - m3range = -j3:j3 - V1 = Array{Float64}(undef, length(m1range),length(m2range),length(m3range)) - V2 = Array{Float64}(undef, length(m1range),length(m2range),length(m3range)) - for J in max(abs(j1-j2-j3),abs(j1-j2+j3),abs(j1+j2-j3)):(j1+j2+j3) - J12range = max(abs(j1-j2),abs(J-j3)):min((j1+j2),(J+j3)) - J23range = max(abs(j2-j3),abs(j1-J)):min((j2+j3),(j1+J)) - for J12 in J12range, J23 in J23range - M = rand(-J:J) # only test for one instance of M in -J:J, should be independent of M anyway - fill!(V1,0) - fill!(V2,0) - for (k1,m1) in enumerate(m1range), (k2,m2) in enumerate(m2range) - abs(m1+m2)<=J12 || continue - for (k3,m3) in enumerate(m3range) - abs(m2+m3)<=J23 || continue - m1+m2+m3==M || continue - V1[k1,k2,k3] = clebschgordan(j1,m1,j2,m2,J12) * - clebschgordan(J12,m1+m2,j3,m3,J) - V2[k1,k2,k3] = clebschgordan(j2,m2,j3,m3,J23) * - clebschgordan(j1,m1,J23,m2+m3,J) +@threads for i = 1:N + @testset "recoupling relation between 3j/CG and 6j/racahW symbols, thread $i" begin + smallerjlist = 0:1//2:5 + for k = i:N:length(smallerjlist) + j1 = smallerjlist[k] + for j2 in smallerjlist, j3 in smallerjlist + m1range = -j1:j1 + m2range = -j2:j2 + m3range = -j3:j3 + V1 = Array{Float64}(undef, length(m1range),length(m2range),length(m3range)) + V2 = Array{Float64}(undef, length(m1range),length(m2range),length(m3range)) + for J in max(abs(j1-j2-j3),abs(j1-j2+j3),abs(j1+j2-j3)):(j1+j2+j3) + J12range = max(abs(j1-j2),abs(J-j3)):min((j1+j2),(J+j3)) + J23range = max(abs(j2-j3),abs(j1-J)):min((j2+j3),(j1+J)) + for J12 in J12range, J23 in J23range + M = rand(-J:J) # only test for one instance of M in -J:J + # should be independent of M anyway + fill!(V1,0) + fill!(V2,0) + for (k1,m1) in enumerate(m1range), (k2,m2) in enumerate(m2range) + abs(m1+m2)<=J12 || continue + for (k3,m3) in enumerate(m3range) + abs(m2+m3)<=J23 || continue + m1+m2+m3==M || continue + V1[k1,k2,k3] = clebschgordan(j1,m1,j2,m2,J12) * + clebschgordan(J12,m1+m2,j3,m3,J) + V2[k1,k2,k3] = clebschgordan(j2,m2,j3,m3,J23) * + clebschgordan(j1,m1,J23,m2+m3,J) + end + end + @test racahW(j1,j2,J,j3,J12,J23) ≈ + dot(V2,V1)/sqrt((2*J12+1)*(2*J23+1)) atol=10*eps(Float64) end end - @test racahW(j1,j2,J,j3,J12,J23) ≈ - dot(V2,V1)/sqrt((2*J12+1)*(2*J23+1)) atol=10*eps(Float64) end end end