diff --git a/LICENSE.md b/LICENSE.md index 77b3272..3f79163 100644 --- a/LICENSE.md +++ b/LICENSE.md @@ -1,17 +1,17 @@ The WignerSymbols.jl package is licensed under the MIT "Expat" License: -> Copyright (c) 2017: Jutho. -> +> Copyright (c) 2017: Jutho Haegeman. +> > Permission is hereby granted, free of charge, to any person obtaining a copy > of this software and associated documentation files (the "Software"), to deal > in the Software without restriction, including without limitation the rights > to use, copy, modify, merge, publish, distribute, sublicense, and/or sell > copies of the Software, and to permit persons to whom the Software is > furnished to do so, subject to the following conditions: -> +> > The above copyright notice and this permission notice shall be included in all > copies or substantial portions of the Software. -> +> > THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR > IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, > FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE @@ -19,4 +19,4 @@ The WignerSymbols.jl package is licensed under the MIT "Expat" License: > LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, > OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE > SOFTWARE. -> +> diff --git a/README.md b/README.md index dfe18c4..270c620 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,53 @@ # WignerSymbols -[![Build Status](https://travis-ci.org/jutho/WignerSymbols.jl.svg?branch=master)](https://travis-ci.org/jutho/WignerSymbols.jl) +[![Build Status](https://travis-ci.org/Jutho/WignerSymbols.jl.svg?branch=master)](https://travis-ci.org/jutho/WignerSymbols.jl) +[![License](http://img.shields.io/badge/license-MIT-brightgreen.svg?style=flat)](LICENSE.md) +[![Coverage Status](https://coveralls.io/repos/Jutho/WignerSymbols.jl/badge.svg?branch=master&service=github)](https://coveralls.io/github/Jutho/WignerSymbols.jl?branch=master) +[![codecov.io](http://codecov.io/github/jutho/WignerSymbols.jl/coverage.svg?branch=master)](http://codecov.io/github/Jutho/WignerSymbols.jl?branch=master) -[![Coverage Status](https://coveralls.io/repos/jutho/WignerSymbols.jl/badge.svg?branch=master&service=github)](https://coveralls.io/github/jutho/WignerSymbols.jl?branch=master) +Compute Wigner's 3j and 6j symbols, and related quantities such as Clebsch-Gordan coefficients and Racah's symbols. -[![codecov.io](http://codecov.io/github/jutho/WignerSymbols.jl/coverage.svg?branch=master)](http://codecov.io/github/jutho/WignerSymbols.jl?branch=master) +## Requirements +This requires a recent master edition of Julia (i.e. v0.7.0-DEV), because it depends on some changes in `Base.GMP`. In particular, it uses the mutating functions for reducing allocation overhead while working with `BigInts` (namely [JuliaLang/julia#21654](https://github.com/JuliaLang/julia/pull/21654)). It also depends on `Primes.jl` for generating prime numbers. + +## Installation +Until it is registered, install via `Pkg.clone("https://github.com/Jutho/WignerSymbols.jl.git")`. + +## 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₁) -> ::T` +* `wigner6j(T::Type{<:AbstractFloat} = Float64, j₁, j₂, j₃, j₄, j₅, j₆) -> ::T` +* `clebschgordan(T::Type{<:AbstractFloat} = Float64, j₁, j₂, j₃, m₁, m₂, m₃ = m₁+m₂) -> ::T` +* `racahV(T::Type{<:AbstractFloat} = Float64, j₁, j₂, j₃, m₁, m₂, m₃ = -m₁-m₂) -> ::T` +* `racahW(T::Type{<:AbstractFloat} = Float64, j₁, j₂, J, j₃, J₁₂, J₂₃) -> ::T` +* `δ(j₁, j₂, j₃) -> ::Bool` +* `Δ(T::Type{<:AbstractFloat} = Float64, j₁, j₂, j₃) -> ::T` + +## 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](https://doi.org/10.1137/15M1021908) ([arXiv:1504.08329](https://arxiv.org/abs/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 `BigInt`s) 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](https://doi.org/10.1137/S1064827503422932) + +for caching the computed 3j and 6j symbols. + +# Benchmark: + +## Todo: +* Wigner 9-j symbols, as explained in [1] and based on + + [3] [L. Wei, New formula for 9-j symbols and their direct calculation, Computers in Physics, 12 (1998), 632–634.](ftp://gravity.physics.uwa.edu.au/pub/Clebsch-Gordan/Papers/Wei9j.pdf) + +* Convenient syntax to get the exact results in the `√(r) * s` format, but in such a way that it can be used by + the Julia type system and can be converted afterwards. diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl new file mode 100644 index 0000000..012fe65 --- /dev/null +++ b/benchmark/benchmarks.jl @@ -0,0 +1,106 @@ +# some random benchmarks; for j=50000 you need more than 8GB of memory +using WignerSymbols + +@time wigner3j(BigFloat, 15, 30, 40, 2, 2, -4) +@time wigner3j(BigFloat, 200, 200, 200, -10, 60, -50) +@time wigner3j(BigFloat, 50000, 50000, 50000, 1000, -6000, 5000) + +j = 8 +@time wigner6j(BigFloat, j, j, j, j, j, j) +j = 200 +@time wigner6j(BigFloat, j, j, j, j, j, j) +j = 600 +@time wigner6j(BigFloat, j, j, j, j, j, j) +j = 10000 +@time wigner6j(BigFloat, j, j, j, j, j, j) +j = 50000 +@time wigner6j(BigFloat, j, j, j, j, j, j) + +function compute3jmax(jmax) + for j₁ = 0:1//2:jmax + for j₂ = 0:1//2:j₁ + for j₃ = abs(j₁-j₂):j₂ + if isinteger(j₁) + m₁ = 0 + M₂ = min(j₂, j₃-m₁) + for m₂ = M₂:(-1):0 + wigner3j(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₂) + end + end + else + for m₁ = 1//2:j₁ + M₂ = min(j₂, j₃-m₁) + for m₂ = -M₂:M₂ + wigner3j(j₁, j₂, j₃, m₁, m₂) + end + end + end + end + end + end + return nothing +end + +@time compute3jmax(10) + +function computeinteger3jmax(jmax) + for j₁ = 0:jmax + for j₂ = 0:j₁ + for j₃ = abs(j₁-j₂):j₂ + m₁ = 0 + M₂ = min(j₂, j₃-m₁) + for m₂ = M₂:(-1):0 + wigner3j(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₂) + end + end + end + end + end + return nothing +end + +function compute6jmax(jmax) + for j₁ = 0:1//2:jmax + for j₂ = 0:1//2:j₁ + for j₃ = abs(j₁-j₂):j₂ + 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₆) + end + end + end + end + end + end + return nothing +end +@time compute6jmax(10) + + +function computeinteger6jmax(jmax) + for j₁ = 0:jmax + for j₂ = 0:j₁ + for j₃ = abs(j₁-j₂):j₂ + 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₆) + end + end + end + end + end + end + return nothing +end