use RationalRoots representation

.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:
+ allow_failures:
-  - julia: nightly
+ - julia: nightly
 
-## uncomment the following lines to override the default test script
+codecov: true
+coveralls: true
 
-after_success:
+# jobs:
-  # push coverage results to Codecov and Coveralls
+#   include:
-  - julia -e 'using Pkg; Pkg.add("Coverage"); using Coverage; Codecov.submit(Codecov.process_folder()); Coveralls.submit(Coveralls.process_folder())'
+#     - 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

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"
 

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:

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

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₃)
+Δ(j₁, j₂, j₃) = Δ(RRBig, j₁, j₂, j₃)
-function Δ(T::Type{<:AbstractFloat}, 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₃)
+wigner3j(j₁, j₂, j₃, m₁, m₂, m₃ = -m₁-m₂) = wigner3j(RRBig, j₁, j₂, j₃, m₁, m₂, m₃)
-function wigner3j(T::Type{<:AbstractFloat}, j₁, j₂, j₃, m₁, m₂, 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
-
+    return _convert(T, sgn*s)*convert(T, signedroot(r))
-    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
 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₃)
+    clebschgordan(RRBig, j₁, m₁, j₂, m₂, j₃, m₃)
-function clebschgordan(T::Type{<:AbstractFloat}, j₁, m₁, j₂, m₂, j₃, m₃ = m₁+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₃)
+racahV(j₁, j₂, j₃, m₁, m₂, m₃ = -m₁-m₂) = racahV(RRBig, j₁, j₂, j₃, m₁, m₂, m₃)
-function racahV(T::Type{<:AbstractFloat}, j₁, j₂, j₃, m₁, m₂, 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₆)
+wigner6j(j₁, j₂, j₃, j₄, j₅, j₆) = wigner6j(RRBig, j₁, j₂, j₃, j₄, j₅, j₆)
-function wigner6j(T::Type{<:AbstractFloat}, 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
-
+    return _convert(T, s) * convert(T, signedroot(r))
-    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
 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₂₃)
+racahW(j₁, j₂, J, j₃, J₁₂, J₂₃) = racahW(RRBig, j₁, j₂, J, j₃, J₁₂, J₂₃)
-function racahW(T::Type{<:AbstractFloat}, 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

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