mirror of
https://github.com/tgorordo/WignerSymbols.jl.git
synced 2026-06-13 02:02:14 -07:00
add MPZ for 0.6 compatibility; update REQUIRE, add testsets
This commit is contained in:
Jutho Haegeman
parent
cb09b427c9
commit
258246f854
5 changed files with 249 additions and 114 deletions
|
|
@ -1,9 +1,9 @@
|
|||
# WignerSymbols
|
||||
|
||||
[](https://travis-ci.org/jutho/WignerSymbols.jl)
|
||||
[](https://travis-ci.org/Jutho/WignerSymbols.jl)
|
||||
[](LICENSE.md)
|
||||
[](https://coveralls.io/github/Jutho/WignerSymbols.jl?branch=master)
|
||||
[](http://codecov.io/github/Jutho/WignerSymbols.jl?branch=master)
|
||||
[](http://codecov.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.
|
||||
|
||||
|
|
|
|||
2
REQUIRE
2
REQUIRE
|
|
@ -1,2 +1,2 @@
|
|||
julia 0.7-
|
||||
julia 0.6
|
||||
Primes
|
||||
|
|
|
|||
116
src/mpz.jl
Normal file
116
src/mpz.jl
Normal file
|
|
@ -0,0 +1,116 @@
|
|||
# taken from Julia 0.7-Dev base library
|
||||
|
||||
module MPZ
|
||||
# wrapping of libgmp functions
|
||||
# - "output parameters" are labeled x, y, z, and are returned when appropriate
|
||||
# - constant input parameters are labeled a, b, c
|
||||
# - a method modifying its input has a "!" appendend to its name, according to Julia's conventions
|
||||
# - some convenient methods are added (in addition to the pure MPZ ones), e.g. `add(a, b) = add!(BigInt(), a, b)`
|
||||
# and `add!(x, a) = add!(x, x, a)`.
|
||||
using Base.GMP: BigInt, Limb
|
||||
|
||||
const mpz_t = Ref{BigInt}
|
||||
const bitcnt_t = Culong
|
||||
|
||||
gmpz(op::Symbol) = (Symbol(:__gmpz_, op), :libgmp)
|
||||
|
||||
init!(x::BigInt) = (ccall((:__gmpz_init, :libgmp), Void, (mpz_t,), x); x)
|
||||
init2!(x::BigInt, a) = (ccall((:__gmpz_init2, :libgmp), Void, (mpz_t, bitcnt_t), x, a); x)
|
||||
|
||||
realloc2!(x, a) = (ccall((:__gmpz_realloc2, :libgmp), Void, (mpz_t, bitcnt_t), x, a); x)
|
||||
realloc2(a) = realloc2!(BigInt(), a)
|
||||
|
||||
sizeinbase(a::BigInt, b) = Int(ccall((:__gmpz_sizeinbase, :libgmp), Csize_t, (mpz_t, Cint), a, b))
|
||||
|
||||
for op in (:add, :sub, :mul, :fdiv_q, :tdiv_q, :fdiv_r, :tdiv_r, :gcd, :lcm, :and, :ior, :xor)
|
||||
op! = Symbol(op, :!)
|
||||
@eval begin
|
||||
$op!(x::BigInt, a::BigInt, b::BigInt) = (ccall($(gmpz(op)), Void, (mpz_t, mpz_t, mpz_t), x, a, b); x)
|
||||
$op(a::BigInt, b::BigInt) = $op!(BigInt(), a, b)
|
||||
$op!(x::BigInt, b::BigInt) = $op!(x, x, b)
|
||||
end
|
||||
end
|
||||
|
||||
invert!(x::BigInt, a::BigInt, b::BigInt) =
|
||||
ccall((:__gmpz_invert, :libgmp), Cint, (mpz_t, mpz_t, mpz_t), x, a, b)
|
||||
invert(a::BigInt, b::BigInt) = invert!(BigInt(), a, b)
|
||||
invert!(x::BigInt, b::BigInt) = invert!(x, x, b)
|
||||
|
||||
for op in (:add_ui, :sub_ui, :mul_ui, :mul_2exp, :fdiv_q_2exp, :pow_ui, :bin_ui)
|
||||
op! = Symbol(op, :!)
|
||||
@eval begin
|
||||
$op!(x::BigInt, a::BigInt, b) = (ccall($(gmpz(op)), Void, (mpz_t, mpz_t, Culong), x, a, b); x)
|
||||
$op(a::BigInt, b) = $op!(BigInt(), a, b)
|
||||
$op!(x::BigInt, b) = $op!(x, x, b)
|
||||
end
|
||||
end
|
||||
|
||||
ui_sub!(x::BigInt, a, b::BigInt) = (ccall((:__gmpz_ui_sub, :libgmp), Void, (mpz_t, Culong, mpz_t), x, a, b); x)
|
||||
ui_sub(a, b::BigInt) = ui_sub!(BigInt(), a, b)
|
||||
|
||||
for op in (:scan1, :scan0)
|
||||
@eval $op(a::BigInt, b) = Int(ccall($(gmpz(op)), Culong, (mpz_t, Culong), a, b))
|
||||
end
|
||||
|
||||
mul_si!(x::BigInt, a::BigInt, b) = (ccall((:__gmpz_mul_si, :libgmp), Void, (mpz_t, mpz_t, Clong), x, a, b); x)
|
||||
mul_si(a::BigInt, b) = mul_si!(BigInt(), a, b)
|
||||
mul_si!(x::BigInt, b) = mul_si!(x, x, b)
|
||||
|
||||
for op in (:neg, :com, :sqrt, :set)
|
||||
op! = Symbol(op, :!)
|
||||
@eval begin
|
||||
$op!(x::BigInt, a::BigInt) = (ccall($(gmpz(op)), Void, (mpz_t, mpz_t), x, a); x)
|
||||
$op(a::BigInt) = $op!(BigInt(), a)
|
||||
end
|
||||
op == :set && continue # MPZ.set!(x) would make no sense
|
||||
@eval $op!(x::BigInt) = $op!(x, x)
|
||||
end
|
||||
|
||||
for (op, T) in ((:fac_ui, Culong), (:set_ui, Culong), (:set_si, Clong), (:set_d, Cdouble))
|
||||
op! = Symbol(op, :!)
|
||||
@eval begin
|
||||
$op!(x::BigInt, a) = (ccall($(gmpz(op)), Void, (mpz_t, $T), x, a); x)
|
||||
$op(a) = $op!(BigInt(), a)
|
||||
end
|
||||
end
|
||||
|
||||
popcount(a::BigInt) = Int(ccall((:__gmpz_popcount, :libgmp), Culong, (mpz_t,), a))
|
||||
|
||||
mpn_popcount(d::Ptr{Limb}, s::Integer) = Int(ccall((:__gmpn_popcount, :libgmp), Culong, (Ptr{Limb}, Csize_t), d, s))
|
||||
mpn_popcount(a::BigInt) = mpn_popcount(a.d, abs(a.size))
|
||||
|
||||
function tdiv_qr!(x::BigInt, y::BigInt, a::BigInt, b::BigInt)
|
||||
ccall((:__gmpz_tdiv_qr, :libgmp), Void, (mpz_t, mpz_t, mpz_t, mpz_t), x, y, a, b)
|
||||
x, y
|
||||
end
|
||||
tdiv_qr(a::BigInt, b::BigInt) = tdiv_qr!(BigInt(), BigInt(), a, b)
|
||||
|
||||
powm!(x::BigInt, a::BigInt, b::BigInt, c::BigInt) =
|
||||
(ccall((:__gmpz_powm, :libgmp), Void, (mpz_t, mpz_t, mpz_t, mpz_t), x, a, b, c); x)
|
||||
powm(a::BigInt, b::BigInt, c::BigInt) = powm!(BigInt(), a, b, c)
|
||||
powm!(x::BigInt, b::BigInt, c::BigInt) = powm!(x, x, b, c)
|
||||
|
||||
function gcdext!(x::BigInt, y::BigInt, z::BigInt, a::BigInt, b::BigInt)
|
||||
ccall((:__gmpz_gcdext, :libgmp), Void, (mpz_t, mpz_t, mpz_t, mpz_t, mpz_t), x, y, z, a, b)
|
||||
x, y, z
|
||||
end
|
||||
gcdext(a::BigInt, b::BigInt) = gcdext!(BigInt(), BigInt(), BigInt(), a, b)
|
||||
|
||||
cmp(a::BigInt, b::BigInt) = Int(ccall((:__gmpz_cmp, :libgmp), Cint, (mpz_t, mpz_t), a, b))
|
||||
cmp_si(a::BigInt, b) = Int(ccall((:__gmpz_cmp_si, :libgmp), Cint, (mpz_t, Clong), a, b))
|
||||
cmp_ui(a::BigInt, b) = Int(ccall((:__gmpz_cmp_ui, :libgmp), Cint, (mpz_t, Culong), a, b))
|
||||
cmp_d(a::BigInt, b) = Int(ccall((:__gmpz_cmp_d, :libgmp), Cint, (mpz_t, Cdouble), a, b))
|
||||
|
||||
mpn_cmp(a::Ptr{Limb}, b::Ptr{Limb}, c) = ccall((:__gmpn_cmp, :libgmp), Cint, (Ptr{Limb}, Ptr{Limb}, Clong), a, b, c)
|
||||
mpn_cmp(a::BigInt, b::BigInt, c) = mpn_cmp(a.d, b.d, c)
|
||||
|
||||
get_str!(x, a, b::BigInt) = (ccall((:__gmpz_get_str,:libgmp), Ptr{Cchar}, (Ptr{Cchar}, Cint, mpz_t), x, a, b); x)
|
||||
set_str!(x::BigInt, a, b) = Int(ccall((:__gmpz_set_str, :libgmp), Cint, (mpz_t, Ptr{UInt8}, Cint), x, a, b))
|
||||
get_d(a::BigInt) = ccall((:__gmpz_get_d, :libgmp), Cdouble, (mpz_t,), a)
|
||||
|
||||
limbs_write!(x::BigInt, a) = ccall((:__gmpz_limbs_write, :libgmp), Ptr{Limb}, (mpz_t, Clong), x, a)
|
||||
limbs_finish!(x::BigInt, a) = ccall((:__gmpz_limbs_finish, :libgmp), Void, (mpz_t, Clong), x, a)
|
||||
import!(x::BigInt, a, b, c, d, e, f) = ccall((:__gmpz_import, :libgmp), Void,
|
||||
(mpz_t, Csize_t, Cint, Csize_t, Cint, Csize_t, Ptr{Void}), x, a, b, c, d, e, f)
|
||||
|
||||
end # module MPZ
|
||||
|
|
@ -1,6 +1,13 @@
|
|||
using Primes.isprime
|
||||
import Base.divgcd
|
||||
|
||||
if VERSION <= v"0.7.0-DEV.262"
|
||||
include("mpz.jl")
|
||||
using .MPZ
|
||||
else
|
||||
using Base.GMP.MPZ
|
||||
end
|
||||
|
||||
|
||||
const primetable = [2,3,5]
|
||||
const factortable = [UInt8[], UInt8[1], UInt8[0,1], UInt8[2], UInt8[0,0,1]]
|
||||
|
|
@ -115,7 +122,7 @@ function Base.convert(::Type{BigInt}, a::PrimeFactorization)
|
|||
A = big(a.sign)
|
||||
for (n, e) in enumerate(a.powers)
|
||||
if !iszero(e)
|
||||
Base.GMP.MPZ.mul!(A, bigprime(n, e))
|
||||
MPZ.mul!(A, bigprime(n, e))
|
||||
end
|
||||
end
|
||||
return A
|
||||
|
|
@ -236,10 +243,10 @@ function sumlist!(list::Vector{<:PrimeFactorization}, ind = 1:length(list))
|
|||
# do sum
|
||||
s = big(0)
|
||||
for k in ind
|
||||
Base.GMP.MPZ.add!(s, convert(BigInt, list[k]))
|
||||
MPZ.add!(s, convert(BigInt, list[k]))
|
||||
end
|
||||
end
|
||||
return Base.GMP.MPZ.mul!(s, convert(BigInt, g))
|
||||
return MPZ.mul!(s, convert(BigInt, g))
|
||||
end
|
||||
|
||||
# Mutating vector methods that also grow and shrink as required
|
||||
|
|
|
|||
|
|
@ -1,20 +1,26 @@
|
|||
if VERSION < v"0.7.0-DEV.2005"
|
||||
const Test = Base.Test
|
||||
end
|
||||
|
||||
using Test
|
||||
using WignerSymbols
|
||||
using Base.Test
|
||||
|
||||
smalljlist = 0:1//2:10
|
||||
largejlist = 0:1//2:1000
|
||||
# test triangle coefficient
|
||||
for j1 in smalljlist, j2 in smalljlist
|
||||
|
||||
@testset "triangle coefficient" begin
|
||||
for j1 in smalljlist, j2 in smalljlist
|
||||
for j3 = abs(j1-j2):(j1+j2)
|
||||
@test Δ(j1,j2,j3) ≈ sqrt(factorial(float(j1+j2-j3))*factorial(float(j1-j2+j3))*
|
||||
factorial(float(j2+j3-j1))/factorial(float(j1+j2+j3+1)))
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
# test 3j:
|
||||
#--------
|
||||
# test cg orthogonality
|
||||
for j1 in smalljlist, j2 in smalljlist
|
||||
@testset "clebschgordan: test orthogonality" begin
|
||||
for j1 in smalljlist, j2 in smalljlist
|
||||
d1 = 2*j1+1
|
||||
d2 = 2*j2+1
|
||||
M = zeros(d1*d2, d1*d2)
|
||||
|
|
@ -30,10 +36,12 @@ for j1 in smalljlist, j2 in smalljlist
|
|||
ind1 += 1
|
||||
end
|
||||
@test M'*M ≈ one(M)
|
||||
end
|
||||
end
|
||||
|
||||
# test recurrence relations: Phys Rev E 57, 7274 (1998)
|
||||
for k = 1:10
|
||||
@testset "wigner3j: test recurrence relations" begin
|
||||
for k = 1:10
|
||||
j2 = convert(BigFloat,rand(0:1//2:1000))
|
||||
j3 = convert(BigFloat,rand(0:1//2:1000))
|
||||
m2 = convert(BigFloat,rand(-j2:j2))
|
||||
|
|
@ -46,12 +54,13 @@ for k = 1:10
|
|||
tol = 10*max(abs(X),abs(Y),abs(Z))*eps(BigFloat)
|
||||
@test (X*wigner3j(BigFloat,j+1,j2,j3,-m2-m3,m2,m3) + Z*wigner3j(BigFloat,j-1,j2,j3,-m2-m3,m2,m3))≈(-Y*wigner3j(BigFloat,j,j2,j3,-m2-m3,m2,m3)) atol=tol
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
# test 6j
|
||||
#----------
|
||||
# test orthogonality
|
||||
for j1 in smalljlist, j2 in smalljlist, j4 in smalljlist
|
||||
@testset "wigner6j: test orthogonality" begin
|
||||
for j1 in smalljlist, j2 in smalljlist, j4 in smalljlist
|
||||
for j5 in max(abs(j1-j2-j4),abs(j1-j2+j4),abs(j1+j2-j4)):(j1+j2+j4)
|
||||
j6range = max(abs(j2-j4),abs(j1-j5)):min((j2+j4),(j1+j5))
|
||||
j3range = max(abs(j1-j2),abs(j4-j5)):min((j1+j2),(j4+j5))
|
||||
|
|
@ -64,20 +73,22 @@ for j1 in smalljlist, j2 in smalljlist, j4 in smalljlist
|
|||
end
|
||||
@test M'*M ≈ one(M)
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
# test special case
|
||||
for j1 in smalljlist, j2 in smalljlist
|
||||
@testset "wigner6j: test special cases" begin
|
||||
for j1 in smalljlist, j2 in smalljlist
|
||||
j6 = 0
|
||||
j4 = j2
|
||||
j5 = j1
|
||||
for j3 in abs(j1-j2):(j1+j2)
|
||||
@test wigner6j(j1,j2,j3,j4,j5,j6) ≈ (-1)^(j1+j2+j3)/sqrt((2*j1+1)*(2*j2+1))
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
# test recurrence relations: Phys Rev E 57, 7274 (1998)
|
||||
for k = 1:10
|
||||
@testset "wigner6j: test recurrence relation" begin
|
||||
for k = 1:10
|
||||
j2 = convert(BigFloat,rand(largejlist))
|
||||
j3 = convert(BigFloat,rand(largejlist))
|
||||
l1 = convert(BigFloat,rand(largejlist))
|
||||
|
|
@ -93,12 +104,12 @@ for k = 1:10
|
|||
tol = 10*max(abs(X),abs(Y),abs(Z))*eps(BigFloat)
|
||||
@test (X*wigner6j(BigFloat,j+1,j2,j3,l1,l2,l3) + Z*wigner6j(BigFloat,j-1,j2,j3,l1,l2,l3))≈(-Y*wigner6j(BigFloat,j,j2,j3,l1,l2,l3)) atol=tol
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
# test recoupling, i.e. relation between 3j and 6j symbols (but use Clebsch-Gordan and RacahW)
|
||||
#---------------------------------------------------------------------------------------------
|
||||
smallerjlist = 0:1//2:5
|
||||
for j1 in smallerjlist, j2 in smallerjlist, j3 in smallerjlist
|
||||
@testset "test recoupling relation between 3j/clebschgordan and 6j/racahW symbols" begin
|
||||
smallerjlist = 0:1//2:5
|
||||
for j1 in smallerjlist, j2 in smallerjlist, j3 in smallerjlist
|
||||
m1range = -j1:j1
|
||||
m2range = -j2:j2
|
||||
m3range = -j3:j3
|
||||
|
|
@ -125,4 +136,5 @@ for j1 in smallerjlist, j2 in smallerjlist, j3 in smallerjlist
|
|||
@test racahW(j1,j2,J,j3,J12,J23) ≈ vecdot(V2,V1)/sqrt((2*J12+1)*(2*J23+1)) atol=10*eps(Float64)
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue