6ddbd340b0
As HalfInteger <: Real, there should be no reason to convert anything in this situation. It happens because the convert method resorts to Float64 as an intermediate value. To still get conversion to floats, we can just dispatch on AbstractFloat instead. However, it should be better to convert the numerator to T first and then divide, so that we would not use a potentially lower precision intermediate value. This solves the problem where calling sum on an vector of HalfIntegers yields a floating point value, even though there is no reason to convert in the summation: julia> sum([HalfInteger(1//2), HalfInteger(3//2)]) 2.0 This is because there is an implicit convert(::Real) in the Base.add_sum function. With this patch the sum call correctly yields a HalfInteger. It also updates the tests related to HalfInteger convert methods: - Make sure that the convert tests also check types - Add a few tests for converting out of HalfInteger |
||
|---|---|---|
| benchmark | ||
| src | ||
| test | ||
| .gitignore | ||
| .travis.yml | ||
| appveyor.yml | ||
| LICENSE.md | ||
| README.md | ||
| REQUIRE | ||
WignerSymbols
Compute Wigner's 3j and 6j symbols, and related quantities such as Clebsch-Gordan coefficients and Racah's symbols.
Requirements
Latest version is compatible with Julia v0.7 only, but older versions can be installed on Julia v0.6.
Installation
Install with the new package manager via ]add WignerSymbols or
using Pkg
Pkg.add("WignerSymbols")
Available functions
While the following function signatures are probably self-explanatory, you can query help for them in the Julia REPL to get further details.
wigner3j(T::Type{<:AbstractFloat} = Float64, j₁, j₂, j₃, m₁, m₂, m₃ = -m₂-m₁) -> ::Twigner6j(T::Type{<:AbstractFloat} = Float64, j₁, j₂, j₃, j₄, j₅, j₆) -> ::Tclebschgordan(T::Type{<:AbstractFloat} = Float64, j₁, m₁, j₂, m₂, j₃, m₃ = m₁+m₂) -> ::TracahV(T::Type{<:AbstractFloat} = Float64, j₁, j₂, j₃, m₁, m₂, m₃ = -m₁-m₂) -> ::TracahW(T::Type{<:AbstractFloat} = Float64, j₁, j₂, J, j₃, J₁₂, J₂₃) -> ::Tδ(j₁, j₂, j₃) -> ::BoolΔ(T::Type{<:AbstractFloat} = Float64, j₁, j₂, j₃) -> ::T
The package also defines the HalfInteger type that can be used to represent half-integer values. Construct if as HalfInteger(a::Real) or HalfInteger(numerator::Integer, denominator::Integer). Furthermore, the range operator a:b can be used to create ranges of HalfInteger values (a HalfIntegerRange).
Implementation
Largely based on reading the paper (but not the code):
[1] H. T. Johansson and C. Forssén, SIAM Journal on Scientific Compututing 38 (2016) 376-384 (arXiv:1504.08329)
with some additional modifications to further improve efficiency for large j (angular momenta quantum numbers).
In particular, 3j and 6j symbols are computed exactly, in the format √(r) * s where r and s are exactly computed as Rational{BigInt},
using an intermediate representation based on prime number factorization. As a consequence thereof, all of the above functions can be called
requesting BigFloat precision for the result. There is currently no convenient syntax for obtaining r and s directly (see TODO).
Most intermediate calculations (prime factorizations of numbers and their factorials, conversion between prime powers and BigInts) are cached to improve the efficiency, but this can result in large use of memory when querying Wigner symbols for large values of j.
Also uses ideas from
[2] J. Rasch and A. C. H. Yu, SIAM Journal on Scientific Compututing 25 (2003), 1416–1428
for caching the computed 3j and 6j symbols.
Todo
-
Wigner 9-j symbols, as explained in [1] and based on
-
Convenient syntax to get the exact results in the
√(r) * sformat, but in such a way that it can be used by the Julia type system and can be converted afterwards.