mirror of
https://github.com/tgorordo/WignerSymbols.jl.git
synced 2026-06-13 02:02:14 -07:00
Switch to using HalfIntegers (#6)
* swich to using HalfIntegers, add project.toml, bump version, update CI * add Random and add seed to avoid unlikely test failure
This commit is contained in:
Jutho
GitHub
parent
6ddbd340b0
commit
a5ee2d5fc6
8 changed files with 50 additions and 417 deletions
|
|
@ -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
|
||||
|
||||
"""
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
Loading…
Add table
Add a link
Reference in a new issue