diff --git a/.travis.yml b/.travis.yml index 99bece3..2763dc7 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,5 +1,6 @@ # Documentation: http://docs.travis-ci.com/user/languages/julia/ language: julia + os: - linux - osx @@ -7,16 +8,29 @@ julia: - 1.0 - 1.1 - 1.2 + - 1.3 - nightly +# env: +# - JULIA_NUM_THREADS=1 +# - JULIA_NUM_THREADS=4 + notifications: email: false matrix: - allow_failures: - - julia: nightly + allow_failures: + - julia: nightly -## uncomment the following lines to override the default test script +codecov: true +coveralls: true -after_success: - # push coverage results to Codecov and Coveralls - - julia -e 'using Pkg; Pkg.add("Coverage"); using Coverage; Codecov.submit(Codecov.process_folder()); Coveralls.submit(Coveralls.process_folder())' +# 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 8ccb4bf..0eaa8c6 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ authors = ["Jutho Haegeman"] version = "1.0.0" [deps] +RationalRoots = "308eb6b3-cc68-5ff3-9e97-c3c4da4fa681" HalfIntegers = "f0d1745a-41c9-11e9-1dd9-e5d34d218721" Primes = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae" diff --git a/appveyor.yml b/appveyor.yml index b0b3e53..cf5896a 100644 --- a/appveyor.yml +++ b/appveyor.yml @@ -3,6 +3,7 @@ environment: - julia_version: 1 - julia_version: 1.1 - julia_version: 1.2 + - julia_version: 1.3 - julia_version: nightly platform: diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl index 012fe65..c8692b6 100644 --- a/benchmark/benchmarks.jl +++ b/benchmark/benchmarks.jl @@ -24,19 +24,19 @@ function compute3jmax(jmax) m₁ = 0 M₂ = min(j₂, j₃-m₁) for m₂ = M₂:(-1):0 - wigner3j(j₁,j₂,j₃,0, m₂) + wigner3j(Float64, j₁,j₂,j₃,0, m₂) end for m₁ = 1:j₁ M₂ = min(j₂, j₃-m₁) for m₂ = -M₂:M₂ - wigner3j(j₁, j₂, j₃, m₁, m₂) + wigner3j(Float64, j₁, j₂, j₃, m₁, m₂) end end else for m₁ = 1//2:j₁ M₂ = min(j₂, j₃-m₁) for m₂ = -M₂:M₂ - wigner3j(j₁, j₂, j₃, m₁, m₂) + wigner3j(Float64, j₁, j₂, j₃, m₁, m₂) end end end @@ -55,12 +55,12 @@ function computeinteger3jmax(jmax) m₁ = 0 M₂ = min(j₂, j₃-m₁) for m₂ = M₂:(-1):0 - wigner3j(j₁,j₂,j₃,0, m₂) + wigner3j(Float64, j₁,j₂,j₃,0, m₂) end for m₁ = 1:j₁ M₂ = min(j₂, j₃-m₁) for m₂ = -M₂:M₂ - wigner3j(j₁, j₂, j₃, m₁, m₂) + wigner3j(Float64, j₁, j₂, j₃, m₁, m₂) end end end @@ -76,7 +76,7 @@ function compute6jmax(jmax) for j₄ = 0:1//2:jmax for j₅ = abs(j₃-j₄):min(jmax,j₃+j₄) for j₆ = min(abs(j₂-j₄),abs(j₁-j₅)):min(j₂+j₄,j₁+j₅,jmax) - wigner6j(j₁, j₂, j₃, j₄, j₅, j₆) + wigner6j(Float64, j₁, j₂, j₃, j₄, j₅, j₆) end end end @@ -95,7 +95,7 @@ function computeinteger6jmax(jmax) for j₄ = 0:jmax for j₅ = abs(j₃-j₄):min(jmax,j₃+j₄) for j₆ = min(abs(j₂-j₄),abs(j₁-j₅)):min(j₂+j₄,j₁+j₅,jmax) - wigner6j(j₁, j₂, j₃, j₄, j₅, j₆) + wigner6j(Float64, j₁, j₂, j₃, j₄, j₅, j₆) end end end diff --git a/src/WignerSymbols.jl b/src/WignerSymbols.jl index 21cab7f..3cb0cb9 100644 --- a/src/WignerSymbols.jl +++ b/src/WignerSymbols.jl @@ -4,14 +4,15 @@ export δ, Δ, clebschgordan, wigner3j, wigner6j, racahV, racahW, HalfInteger using Base.GMP.MPZ using HalfIntegers +using RationalRoots +const RRBig = RationalRoot{BigInt} +import RationalRoots: _convert include("primefactorization.jl") 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 FASTCUTOFF = convert(BigInt, typemax(Int)) - function __init__() global bigone, bigprimetable, Wigner3j, Wigner6j bigone[] = big(1) @@ -44,8 +45,8 @@ as a type `T` floating point number. Returns `zero(T)` if the triangle condition `δ(j₁, j₂, j₃)` is not satisfied, but 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₃) +Δ(j₁, j₂, j₃) = Δ(RRBig, j₁, j₂, j₃) +function Δ(T::Type{<:Real}, j₁, j₂, j₃) for jᵢ in (j₁, j₂, j₃) (ishalfinteger(jᵢ) && jᵢ >= zero(jᵢ)) || throw(DomainError("invalid jᵢ", jᵢ)) end @@ -53,7 +54,7 @@ function Δ(T::Type{<:AbstractFloat}, j₁, j₂, j₃) return zero(T) end n, d = Δ²(j₁, j₂, j₃) - return sqrt(convert(T, convert(BigInt, n) // convert(BigInt, d))) + return convert(T, signedroot(RationalRoot{BigInt}, n//d)) end """ @@ -67,8 +68,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ᵢ`. """ -wigner3j(j₁, j₂, j₃, m₁, m₂, m₃ = -m₁-m₂) = wigner3j(Float64, j₁, j₂, j₃, m₁, m₂, m₃) -function wigner3j(T::Type{<:AbstractFloat}, j₁, j₂, j₃, m₁, m₂, m₃ = -m₁-m₂) +wigner3j(j₁, j₂, j₃, m₁, m₂, m₃ = -m₁-m₂) = wigner3j(RRBig, j₁, j₂, j₃, m₁, m₂, m₃) +function wigner3j(T::Type{<:Real}, j₁, j₂, j₃, m₁, m₂, m₃ = -m₁-m₂) # check angular momenta for (jᵢ,mᵢ) in ((j₁, m₁), (j₂, m₂), (j₃, m₃)) ϵ(jᵢ, mᵢ) || throw(DomainError((jᵢ, mᵢ), "invalid combination (jᵢ, mᵢ)")) @@ -104,13 +105,7 @@ function wigner3j(T::Type{<:AbstractFloat}, j₁, j₂, j₃, m₁, m₂, m₃ = s *= compute3jseries(β₁, β₂, β₃, α₁, α₂) Wigner3j[(β₁, β₂, β₃, α₁, α₂)] = (r,s) end - - if T != BigFloat && all((s.num, s.den, r.num, r.den) .< FASTCUTOFF) - sn, sd, rn, rd = convert.(T, (s.num, s.den, r.num, r.den)) - return sgn*(sn/sd)*sqrt(rn/rd) - else - return convert(T, sgn*s*sqrt(r)) - end + return _convert(T, sgn*s)*convert(T, signedroot(r)) end """ @@ -123,11 +118,11 @@ Returns `zero(T)` if the triangle condition `δ(j₁, j₂, j₃)` is not satisf 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₃) -function clebschgordan(T::Type{<:AbstractFloat}, j₁, m₁, j₂, m₂, j₃, m₃ = m₁+m₂) + clebschgordan(RRBig, j₁, m₁, j₂, m₂, j₃, m₃) +function clebschgordan(T::Type{<:Real}, j₁, m₁, j₂, m₂, j₃, m₃ = m₁+m₂) s = wigner3j(T, j₁, j₂, j₃, m₁, m₂, -m₃) iszero(s) && return s - s *= sqrt(convert(T, j₃+j₃+one(j₃))) + s *= convert(T, signedroot(RRBig, j₃+j₃+one(j₃))) return isodd(convert(Int,j₁ - j₂ + m₃)) ? -s : s end @@ -142,8 +137,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ᵢ`. """ -racahV(j₁, j₂, j₃, m₁, m₂, m₃ = -m₁-m₂) = racahV(Float64, j₁, j₂, j₃, m₁, m₂, m₃) -function racahV(T::Type{<:AbstractFloat}, j₁, j₂, j₃, m₁, m₂, m₃ = -m₁-m₂) +racahV(j₁, j₂, j₃, m₁, m₂, m₃ = -m₁-m₂) = racahV(RRBig, j₁, j₂, j₃, m₁, m₂, m₃) +function racahV(T::Type{<:Real}, j₁, j₂, j₃, m₁, m₂, m₃ = -m₁-m₂) s = wigner3j(T, j₁, j₂, j₃, m₁, m₂, m₃) return isodd(convert(Int, -j₁ + j₂ + j₃)) ? -s : s end @@ -160,8 +155,8 @@ Returns `zero(T)` if any of triangle conditions `δ(j₁, j₂, j₃)`, `δ(j₁ `δ(j₂, j₄, j₆)`, `δ(j₃, j₄, j₅)` are not satisfied, but throws a `DomainError` if the `jᵢ`s are not integer or halfinteger. """ -wigner6j(j₁, j₂, j₃, j₄, j₅, j₆) = wigner6j(Float64, j₁, j₂, j₃, j₄, j₅, j₆) -function wigner6j(T::Type{<:AbstractFloat}, j₁, j₂, j₃, j₄, j₅, j₆) +wigner6j(j₁, j₂, j₃, j₄, j₅, j₆) = wigner6j(RRBig, j₁, j₂, j₃, j₄, j₅, j₆) +function wigner6j(T::Type{<:Real}, j₁, j₂, j₃, j₄, j₅, j₆) # check validity of `jᵢ`s for jᵢ in (j₁, j₂, j₃, j₄, j₅, j₆) (ishalfinteger(jᵢ) && jᵢ >= zero(jᵢ)) || throw(DomainError("invalid jᵢ", jᵢ)) @@ -201,19 +196,15 @@ function wigner6j(T::Type{<:AbstractFloat}, 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) s *= compute6jseries(β₁, β₂, β₃, α₁, α₂, α₃, α₄) Wigner6j[(β₁, β₂, β₃, α₁, α₂, α₃)] = (r, s) end - - if T != BigFloat && all((s.num, s.den, r.num, r.den) .< FASTCUTOFF) - sn, sd, rn, rd = convert.(T, (s.num, s.den, r.num, r.den)) - return (sn/sd)*sqrt(rn/rd) - else - return convert(T, s*sqrt(r)) - end + return _convert(T, s) * convert(T, signedroot(r)) end """ @@ -227,8 +218,8 @@ Returns `zero(T)` if any of triangle conditions `δ(j₁, j₂, J₁₂)`, `δ(j `δ(j₁, J₂₃, J)`, `δ(J₁₂, j₃, J)` are not satisfied, but throws a `DomainError` if the `jᵢ`s and `J`s are not integer or halfinteger. """ -racahW(j₁, j₂, J, j₃, J₁₂, J₂₃) = racahW(Float64, j₁, j₂, J, j₃, J₁₂, J₂₃) -function racahW(T::Type{<:AbstractFloat}, j₁, j₂, J, j₃, J₁₂, J₂₃) +racahW(j₁, j₂, J, j₃, J₁₂, J₂₃) = racahW(RRBig, j₁, j₂, J, j₃, J₁₂, J₂₃) +function racahW(T::Type{<:Real}, j₁, j₂, J, j₃, J₁₂, J₂₃) s = wigner6j(T, j₁, j₂, J₁₂, j₃, J, J₂₃) if !iszero(s) && isodd(convert(Int, j₁ + j₂ + j₃ + J)) return -s diff --git a/src/primefactorization.jl b/src/primefactorization.jl index 556dafa..a173e8f 100644 --- a/src/primefactorization.jl +++ b/src/primefactorization.jl @@ -116,13 +116,13 @@ Base.promote_rule(P::Type{<:PrimeFactorization},::Type{BigInt}) = BigInt Base.convert(P::Type{<:PrimeFactorization}, n::Integer) = convert(P, primefactor(n)) function Base.convert(::Type{BigInt}, a::PrimeFactorization) - A = big(a.sign) + A = one(BigInt) for (n, e) in enumerate(a.powers) if !iszero(e) MPZ.mul!(A, bigprime(n, e)) end end - return A + return a.sign < 0 ? MPZ.neg!(A) : A end Base.convert(::Type{PrimeFactorization{U}}, a::PrimeFactorization{U}) where {U<:Unsigned} = a