diff --git a/.travis.yml b/.travis.yml index 014e4eb..99bece3 100644 --- a/.travis.yml +++ b/.travis.yml @@ -4,9 +4,9 @@ os: - linux - osx julia: - - 0.7 - 1.0 - 1.1 + - 1.2 - nightly notifications: email: false diff --git a/Project.toml b/Project.toml new file mode 100644 index 0000000..8ccb4bf --- /dev/null +++ b/Project.toml @@ -0,0 +1,19 @@ +name = "WignerSymbols" +uuid = "9f57e263-0b3d-5e2e-b1be-24f2bb48858b" +authors = ["Jutho Haegeman"] +version = "1.0.0" + +[deps] +HalfIntegers = "f0d1745a-41c9-11e9-1dd9-e5d34d218721" +Primes = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae" + +[compat] +julia = "1" + +[extras] +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" + +[targets] +test = ["Test", "LinearAlgebra", "Random"] diff --git a/REQUIRE b/REQUIRE deleted file mode 100644 index 61680a6..0000000 --- a/REQUIRE +++ /dev/null @@ -1,2 +0,0 @@ -julia 0.7 -Primes diff --git a/appveyor.yml b/appveyor.yml index d2c4499..b0b3e53 100644 --- a/appveyor.yml +++ b/appveyor.yml @@ -1,7 +1,8 @@ environment: matrix: - - julia_version: 0.7 - julia_version: 1 + - julia_version: 1.1 + - julia_version: 1.2 - julia_version: nightly platform: diff --git a/src/WignerSymbols.jl b/src/WignerSymbols.jl index eeb6d3b..59221a8 100644 --- a/src/WignerSymbols.jl +++ b/src/WignerSymbols.jl @@ -3,8 +3,8 @@ module WignerSymbols export δ, Δ, clebschgordan, wigner3j, wigner6j, racahV, racahW, HalfInteger using Base.GMP.MPZ +using HalfIntegers -include("halfinteger.jl") include("primefactorization.jl") const Wigner3j = Dict{Tuple{UInt,UInt,UInt,Int,Int},Tuple{Rational{BigInt},Rational{BigInt}}}() @@ -29,13 +29,7 @@ end Checks the triangle conditions `j₃ <= j₁ + j₂`, `j₁ <= j₂ + j₃` and `j₂ <= j₃ + j₁`. """ -function δ(j₁, j₂, j₃) - j₃ <= j₁ + j₂ || return false - j₁ <= j₂ + j₃ || return false - j₂ <= j₃ + j₁ || return false - isinteger(j₁+j₂+j₃) || return false - return true -end +δ(j₁, j₂, j₃) = (j₃ <= j₁ + j₂) && (j₁ <= j₂ + j₃) && (j₂ <= j₃ + j₁) && isinteger(j₁+j₂+j₃) # triangle coefficient """ @@ -51,7 +45,7 @@ throws a `DomainError` if the `jᵢ`s are not (half)integer Δ(j₁, j₂, j₃) = Δ(Float64, j₁, j₂, j₃) function Δ(T::Type{<:AbstractFloat}, j₁, j₂, j₃) for jᵢ in (j₁, j₂, j₃) - (ishalfinteger(jᵢ) && jᵢ >= 0) || throw(DomainError("invalid jᵢ", jᵢ)) + (ishalfinteger(jᵢ) && jᵢ >= zero(jᵢ)) || throw(DomainError("invalid jᵢ", jᵢ)) end if !δ(j₁, j₂, j₃) return zero(T) @@ -109,7 +103,8 @@ function wigner3j(T::Type{<:AbstractFloat}, j₁, j₂, j₃, m₁, m₂, m₃ = Wigner3j[(β₁, β₂, β₃, α₁, α₂)] = (r,s) end - return sgn*sqrt(convert(T, r.num)/convert(T, r.den))*(convert(T, s.num)/convert(T, s.den)) + sn, sd, rn, rd = convert.(T, (s.num, s.den, r.num, r.den)) + return sgn*(sn/sd)*sqrt(rn/rd) end """ @@ -121,7 +116,8 @@ as a type `T` floating point number. By default, `T = Float64` and `m₃ = m₁+ Returns `zero(T)` if the triangle condition `δ(j₁, j₂, j₃)` is not satisfied, but throws a `DomainError` if the `jᵢ`s and `mᵢ`s are not (half)integer or `abs(mᵢ) > jᵢ`. """ -clebschgordan(j₁, m₁, j₂, m₂, j₃, m₃ = m₁+m₂) = clebschgordan(Float64, j₁, m₁, j₂, m₂, j₃, m₃) +clebschgordan(j₁, m₁, j₂, m₂, j₃, m₃ = m₁+m₂) = + clebschgordan(Float64, j₁, m₁, j₂, m₂, j₃, m₃) function clebschgordan(T::Type{<:AbstractFloat}, j₁, m₁, j₂, m₂, j₃, m₃ = m₁+m₂) s = wigner3j(T, j₁, j₂, j₃, m₁, m₂, -m₃) iszero(s) && return s @@ -162,13 +158,13 @@ wigner6j(j₁, j₂, j₃, j₄, j₅, j₆) = wigner6j(Float64, j₁, j₂, j function wigner6j(T::Type{<:AbstractFloat}, j₁, j₂, j₃, j₄, j₅, j₆) # check validity of `jᵢ`s for jᵢ in (j₁, j₂, j₃, j₄, j₅, j₆) - (ishalfinteger(jᵢ) && jᵢ >= 0) || throw(DomainError("invalid jᵢ", jᵢ)) + (ishalfinteger(jᵢ) && jᵢ >= zero(jᵢ)) || throw(DomainError("invalid jᵢ", jᵢ)) end - α̂₁ = map(converthalfinteger, (j₁, j₂, j₃)) - α̂₂ = map(converthalfinteger, (j₁, j₆, j₅)) - α̂₃ = map(converthalfinteger, (j₂, j₄, j₆)) - α̂₄ = map(converthalfinteger, (j₃, j₄, j₅)) + α̂₁ = (j₁, j₂, j₃) + α̂₂ = (j₁, j₆, j₅) + α̂₃ = (j₂, j₄, j₆) + α̂₄ = (j₃, j₄, j₅) # check triangle conditions if !(δ(α̂₁...) && δ(α̂₂...) && δ(α̂₃...) && δ(α̂₄...)) @@ -206,7 +202,8 @@ function wigner6j(T::Type{<:AbstractFloat}, j₁, j₂, j₃, j₄, j₅, j₆) Wigner6j[(β₁, β₂, β₃, α₁, α₂, α₃)] = (r, s) end - return sqrt(convert(T, r.num)/convert(T, r.den))*(convert(T, s.num)/convert(T, s.den)) + sn, sd, rn, rd = convert.(T, (s.num, s.den, r.num, r.den)) + return (sn/sd)*sqrt(rn/rd) end """ diff --git a/src/halfinteger.jl b/src/halfinteger.jl deleted file mode 100644 index 0065c80..0000000 --- a/src/halfinteger.jl +++ /dev/null @@ -1,190 +0,0 @@ -# HalfInteger - -""" - struct HalfInteger <: Real - -Represents half-integer values. - ---- - -HalfInteger(numerator::Integer, denominator::Integer) - -Constructs a `HalfInteger` object as a rational number from the given integer numerator -and denominator values. - -# Examples - -```jldoctest -julia> HalfInteger(1, 2) -1/2 - -julia> HalfInteger(-2, 1) --2 -``` -""" -struct HalfInteger <: Real - numerator::Int # with an implicit denominator of 2 - - function HalfInteger(num::Integer, den::Integer) - (den == 2) && return new(num) - (den == 1) && return new(2*num) - (den == 0) && throw(ArgumentError("Denominator can not be zero.")) - # If non-trivial, we'll see if we can reduce it down to a half-integer - numerator, r = divrem(2*num, den) - if r == 0 - return new(numerator) - else - throw(ArgumentError("$num // $den is not a half-integer value.")) - end - end -end - -""" - HalfInteger(x::Real) - -Attempts to create a `HalfInteger` out of the real number `x`. Throws an `InexactError` if -`x` can not be represented as a half-integer value. - -# Examples - -```jldoctest -julia> HalfInteger(3) -3 - -julia> HalfInteger(1.5) -3/2 -``` -""" -HalfInteger(x::Real) = convert(HalfInteger, x) - -Base.promote_rule(::Type{HalfInteger}, ::Type{<:Integer}) = HalfInteger -Base.promote_rule(::Type{HalfInteger}, T::Type{<:Rational}) = T -Base.promote_rule(::Type{HalfInteger}, T::Type{<:Real}) = T - -Base.convert(::Type{HalfInteger}, n::Integer) = HalfInteger(2*n, 2) -function Base.convert(::Type{HalfInteger}, r::Rational) - if r.den == 1 - return HalfInteger(2*r.num, 2) - elseif r.den == 2 - return HalfInteger(r.num, 2) - else - throw(InexactError(:HalfInteger, HalfInteger, r)) - end -end -function Base.convert(::Type{HalfInteger}, r::Real) - num = 2*r - if isinteger(num) - return HalfInteger(convert(Int, num), 2) - else - throw(InexactError(:HalfInteger, HalfInteger, r)) - end -end -Base.convert(T::Type{<:Integer}, s::HalfInteger) = iseven(s.numerator) ? convert(T, s.numerator>>1) : throw(InexactError(Symbol(T), T, s)) -Base.convert(T::Type{<:Rational}, s::HalfInteger) = convert(T, s.numerator//2) -Base.convert(T::Type{<:AbstractFloat}, s::HalfInteger) = convert(T, s.numerator) / T(2) -Base.convert(::Type{HalfInteger}, s::HalfInteger) = s - -# Arithmetic - -Base.:+(a::HalfInteger, b::HalfInteger) = HalfInteger(a.numerator+b.numerator, 2) -Base.:-(a::HalfInteger, b::HalfInteger) = HalfInteger(a.numerator-b.numerator, 2) -Base.:-(a::HalfInteger) = HalfInteger(-a.numerator, 2) -Base.:*(a::Integer, b::HalfInteger) = HalfInteger(a * b.numerator, 2) -Base.:*(a::HalfInteger, b::Integer) = b * a -Base.:<=(a::HalfInteger, b::HalfInteger) = a.numerator <= b.numerator -Base.:<(a::HalfInteger, b::HalfInteger) = a.numerator < b.numerator -Base.one(::Type{HalfInteger}) = HalfInteger(2, 2) -Base.zero(::Type{HalfInteger}) = HalfInteger(0, 2) - -Base.floor(x::HalfInteger) = isinteger(x) ? x : x - HalfInteger(1, 2) -Base.floor(::Type{T}, x::HalfInteger) where T <: Integer = convert(T, floor(x)) - -Base.ceil(x::HalfInteger) = isinteger(x) ? x : x + HalfInteger(1, 2) -Base.ceil(::Type{T}, x::HalfInteger) where T <: Integer = convert(T, ceil(x)) - -# Hashing - -function Base.hash(a::HalfInteger, h::UInt) - iseven(a.numerator) && return hash(a.numerator>>1, h) - num, den = a.numerator, 2 - den = 1 - pow = -1 - if abs(num) < 9007199254740992 - return hash(ldexp(Float64(num),pow), h) - end - h = Base.hash_integer(den, h) - h = Base.hash_integer(pow, h) - h = Base.hash_integer(num, h) - return h -end - -# Parsing and printing - -""" - parse(HalfInteger, s) - -Parses the string `s` into the corresponding `HalfInteger`-value. String can either be a -number or a fraction of the form `n/2`. -""" -function Base.parse(::Type{HalfInteger}, s::AbstractString) - if in('/', s) - num, den = split(s, '/'; limit=2) - parse(Int, den) == 2 || - throw(ArgumentError("Denominator not 2 in HalfInteger string '$s'.")) - HalfInteger(parse(Int, num), 2) - elseif !isempty(strip(s)) - HalfInteger(parse(Int, s)) - else - throw(ArgumentError("input string is empty or only contains whitespace")) - end -end - -Base.show(io::IO, x::HalfInteger) = - print(io, iseven(x.numerator) ? "$(div(x.numerator, 2))" : "$(x.numerator)/2") - -# Other methods - -Base.isinteger(a::HalfInteger) = iseven(a.numerator) -ishalfinteger(a::HalfInteger) = true -ishalfinteger(a::Integer) = true -ishalfinteger(a::Rational) = a.den == 1 || a.den == 2 -ishalfinteger(a::Real) = isinteger(2*a) - -converthalfinteger(a::Number) = convert(HalfInteger, a) - -Base.numerator(a::HalfInteger) = iseven(a.numerator) ? div(a.numerator, 2) : a.numerator -Base.denominator(a::HalfInteger) = iseven(a.numerator) ? 1 : 2 - -# Range of HalfIntegers - -""" - struct HalfIntegerRange <: AbstractVector{HalfInteger} - -A range of `HalfInteger` values from `start` to `stop`, spaced by `1`. The `a:b` syntax -where both `a` and `b` are `HalfInteger`s can also be use to construct this range. -""" -struct HalfIntegerRange <: AbstractVector{HalfInteger} - start :: HalfInteger - stop :: HalfInteger - - function HalfIntegerRange(start::HalfInteger, stop::HalfInteger) - (start <= stop) || - throw(ArgumentError("Second argument must be greater or equal to the first.")) - return new(start, stop) - end -end -Base.iterate(it::HalfIntegerRange) = (it.start, it.start + 1) -Base.iterate(it::HalfIntegerRange, s) = (s <= it.stop) ? (s, s+1) : nothing -Base.length(it::HalfIntegerRange) = floor(Int, it.stop - it.start) + 1 -Base.size(it::HalfIntegerRange) = (length(it),) -function Base.getindex(it::HalfIntegerRange, i::Integer) - 1 <= i <= length(it) || throw(BoundsError(it, i)) - it.start + i - 1 -end - -""" - (:)(i::HalfInteger, j::HalfInteger) - -Constructs a `HalfIntegerRange` out of two `HalfInteger` values. -""" -Base.:(:)(i::HalfInteger, j::HalfInteger) = HalfIntegerRange(i, j) diff --git a/test/halfinteger.jl b/test/halfinteger.jl deleted file mode 100644 index 2198f9c..0000000 --- a/test/halfinteger.jl +++ /dev/null @@ -1,194 +0,0 @@ -using Test -using WignerSymbols: HalfInteger, ishalfinteger, HalfIntegerRange - -@testset "HalfInteger" begin - @testset "HalfInteger type" begin - # HalfInteger constructors - @test HalfInteger(1, 2).numerator == 1 - @test HalfInteger(1, 1).numerator == 2 - @test HalfInteger(0, 1).numerator == 0 - @test HalfInteger(0, 2).numerator == 0 - @test HalfInteger(0, 5).numerator == 0 - @test HalfInteger(10, 5).numerator == 4 - @test HalfInteger(21, 14).numerator == 3 - @test HalfInteger(-3, 2).numerator == -3 - @test HalfInteger(3, -2).numerator == -3 - @test HalfInteger(-3, -2).numerator == 3 - @test_throws ArgumentError HalfInteger(1, 0) - @test_throws ArgumentError HalfInteger(1, 3) - @test_throws ArgumentError HalfInteger(1, -3) - @test_throws ArgumentError HalfInteger(-5, 3) - @test_throws ArgumentError HalfInteger(-1000, -999) - - # convert methods - @test convert(HalfInteger, 2) === HalfInteger(2, 1) - @test convert(HalfInteger, 1//2) === HalfInteger(1, 2) - @test convert(HalfInteger, 1.5) === HalfInteger(3, 2) - @test_throws InexactError convert(HalfInteger, 1//3) - @test_throws InexactError convert(HalfInteger, 0.6) - @test convert(HalfInteger, 2) === HalfInteger(2, 1) - @test convert(HalfInteger, 1//2) === HalfInteger(1, 2) - @test convert(HalfInteger, 1.5) === HalfInteger(3, 2) - - @test convert(Integer, HalfInteger(2, 1)) === 2 - @test_throws InexactError convert(Integer, HalfInteger(1, 2)) - @test convert(Float64, HalfInteger(3, 2)) isa Float64 - @test convert(Float32, HalfInteger(3, 2)) isa Float32 - @test convert(Float64, HalfInteger(3, 2)) == 1.5 - @test convert(Real, HalfInteger(3, 2)) === HalfInteger(3, 2) - - # single-argument constructor - @test HalfInteger(0) == HalfInteger(0, 2) - @test HalfInteger(1) == HalfInteger(1, 1) - @test HalfInteger(2) == HalfInteger(2, 1) - @test HalfInteger(-30) == HalfInteger(-60, 2) - @test HalfInteger(0//2) == HalfInteger(0, 1) - @test HalfInteger(1//2) == HalfInteger(1, 2) - @test HalfInteger(-5//2) == HalfInteger(-5, 2) - end - - a = HalfInteger(2) - b = HalfInteger(3, 2) - - @testset "HalfInteger arithmetic" begin - @test a + b == 2 + 3//2 - @test a - b == 2 - 3//2 - @test zero(a) == 0 - @test one(a) == 1 - @test a > b - @test b < a - @test b <= a - @test a >= b - @test a == a - @test a != b - @test 2 * HalfInteger(0) == HalfInteger(0) - @test 2 * HalfInteger(1, 2) == HalfInteger(1) - @test HalfInteger(1) * 2 == HalfInteger(2) - @test 2 * a == HalfInteger(4) - @test (-1) * b == HalfInteger(-3//2) - - @test floor(HalfInteger(0)) === HalfInteger(0) - @test floor(HalfInteger(-1)) === HalfInteger(-1) - @test floor(HalfInteger(1, 2)) === HalfInteger(0) - @test floor(HalfInteger(-1, 2)) === HalfInteger(-1) - @test floor(Int, HalfInteger(0)) === 0 - @test floor(Int, HalfInteger(1, 2)) === 0 - @test floor(Int32, HalfInteger(-5, 2)) === Int32(-3) - @test floor(Int32, HalfInteger(5)) === Int32(5) - - @test ceil(HalfInteger(0)) === HalfInteger(0) - @test ceil(HalfInteger(-1)) === HalfInteger(-1) - @test ceil(HalfInteger(1, 2)) === HalfInteger(1) - @test ceil(HalfInteger(-1, 2)) === HalfInteger(0) - @test ceil(Int, HalfInteger(0)) === 0 - @test ceil(Int, HalfInteger(1, 2)) === 1 - @test ceil(Int32, HalfInteger(-5, 2)) === Int32(-2) - @test ceil(Int32, HalfInteger(5)) === Int32(5) - - for n in -98:7:98 - halfint, rat = HalfInteger(n, 2), n // 2 - @test halfint == rat - @test halfint == HalfInteger(n / 2) - iseven(n) && @test halfint == HalfInteger(div(n, 2)) - @test ceil(halfint) == ceil(rat) - @test floor(halfint) == floor(rat) - end - end - - @testset "Parsing and printing" begin - @test string(HalfInteger(0)) == "0" - @test string(HalfInteger(1)) == "1" - @test string(HalfInteger(-1)) == "-1" - @test string(HalfInteger(1, 2)) == "1/2" - @test string(HalfInteger(-3, 2)) == "-3/2" - - @test parse(HalfInteger, "0") == HalfInteger(0) - @test parse(HalfInteger, "1") == HalfInteger(1) - @test parse(HalfInteger, "210938") == HalfInteger(210938) - @test parse(HalfInteger, "-15") == HalfInteger(-15) - @test parse(HalfInteger, "1/2") == HalfInteger(1//2) - @test parse(HalfInteger, "-3/2") == HalfInteger(-3//2) - @test_throws ArgumentError parse(HalfInteger, "") - @test_throws ArgumentError parse(HalfInteger, "-50/100") - @test_throws ArgumentError parse(HalfInteger, "1/3") - end - - @testset "HalfInteger hashing" begin - @test hash(a) == hash(2) - @test hash(b) == hash(1.5) - end - - @testset "Other HalfInteger methods" begin - @test isinteger(HalfInteger(0)) - @test isinteger(HalfInteger(1)) - @test !isinteger(HalfInteger(1, 2)) - - @test ishalfinteger(1) - @test ishalfinteger(1.0) - @test ishalfinteger(-0.5) - @test ishalfinteger(HalfInteger(0)) - @test ishalfinteger(HalfInteger(1, 2)) - @test ishalfinteger(1//1) - @test ishalfinteger(1//2) - @test !ishalfinteger(0.3) - @test !ishalfinteger(-5//7) - - @test numerator(HalfInteger(0)) == 0 - @test numerator(HalfInteger(1, 2)) == 1 - @test numerator(HalfInteger(1)) == 1 - @test numerator(HalfInteger(-3, 2)) == -3 - - @test denominator(HalfInteger(0)) == 1 - @test denominator(HalfInteger(1, 2)) == 2 - @test denominator(HalfInteger(1)) == 1 - @test denominator(HalfInteger(-3, 2)) == 2 - end - - @testset "HalfIntegerRange" begin - hi(x) = HalfInteger(x) - - @test length(HalfIntegerRange(hi(0), hi(0))) == 1 - @test length(HalfIntegerRange(hi(0), hi(2))) == 3 - let hirange = HalfIntegerRange(hi(-1//2), hi(1//2)) - @test length(hirange) == 2 - @test size(hirange) == (2,) - @test collect(hirange) == [hi(-1//2), hi(1//2)] - end - let hirange = HalfIntegerRange(hi(0), hi(1//2)) - @test length(hirange) == 1 - @test size(hirange) == (1,) - @test collect(hirange) == [hi(0)] - end - let hirange = HalfIntegerRange(hi(1//2), hi(3)) - @test length(hirange) == 3 - @test size(hirange) == (3,) - @test collect(hirange) == [hi(1//2), hi(3//2), hi(5//2)] - end - - @test hi(5):hi(7) == HalfIntegerRange(hi(5), hi(7)) - @test hi(-1//2):hi(1//2) == HalfIntegerRange(hi(-1//2), hi(1//2)) - - @test collect(hi(0) : hi(2)) == [hi(0), hi(1), hi(2)] - @test collect(hi(-3//2) : hi(1//2)) == [hi(-3//2), hi(-1//2), hi(1//2)] - - let hirange = hi(-3//2):hi(0) - @test length(hirange) == 2 - @test size(hirange) == (2,) - @test collect(hirange) == [hi(-3//2), hi(-1//2)] - end - - @test hi(1//2) ∈ hi(-1//2) : hi(1//2) - @test 1 ∈ hi(0) : hi(2) - @test 1//2 ∈ hi(-1//2) : hi(7//2) - @test !(hi(1//2) ∈ hi(0) : hi(1)) - @test !(1//2 ∈ hi(-1) : hi(7)) - - r = hi(-3//2) : hi(3//2) - @test r[1] == hi(-3//2) - @test r[2] == hi(-1//2) - @test r[3] == hi(1//2) - @test r[4] == hi(3//2) - @test_throws BoundsError r[0] - @test_throws BoundsError r[5] - end -end diff --git a/test/runtests.jl b/test/runtests.jl index 58de51f..5f2ac8f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,8 +1,9 @@ using Test using WignerSymbols using LinearAlgebra +using Random -include("halfinteger.jl") +Random.seed!(1234) smalljlist = 0:1//2:10 largejlist = 0:1//2:1000 @@ -99,8 +100,9 @@ end 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) + 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 @@ -121,15 +123,15 @@ end 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) - for (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 + 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)