diff --git a/README.md b/README.md index d12d284..672d493 100644 --- a/README.md +++ b/README.md @@ -1,9 +1,9 @@ # 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) +[![codecov.io](http://codecov.io/github/Jutho/WignerSymbols.jl/coverage.svg?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. diff --git a/REQUIRE b/REQUIRE index cd4a598..48df863 100644 --- a/REQUIRE +++ b/REQUIRE @@ -1,2 +1,2 @@ -julia 0.7- +julia 0.6 Primes diff --git a/src/mpz.jl b/src/mpz.jl new file mode 100644 index 0000000..a448819 --- /dev/null +++ b/src/mpz.jl @@ -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 diff --git a/src/primefactorization.jl b/src/primefactorization.jl index e478a9c..e0145b5 100644 --- a/src/primefactorization.jl +++ b/src/primefactorization.jl @@ -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 diff --git a/test/runtests.jl b/test/runtests.jl index 4ff1178..70ec8ca 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,128 +1,140 @@ +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 - 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))) + +@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 - d1 = 2*j1+1 - d2 = 2*j2+1 - M = zeros(d1*d2, d1*d2) - ind1 = 1 - for m1 in -j1:j1, m2 in -j2:j2 - ind2 = 1 - for j3 in abs(j1-j2):(j1+j2) - for m3 in -j3:j3 - M[ind1,ind2] = clebschgordan(j1,m1,j2,m2,j3,m3) - ind2 += 1 - end - end - ind1 += 1 - end - @test M'*M ≈ one(M) -end - -# test recurrence relations: Phys Rev E 57, 7274 (1998) -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)) - m3 = convert(BigFloat,rand(-j3:j3)) - - for j in max(abs(j2-j3),abs(m2+m3))+1:(j2+j3)-1 - X = j*sqrt(((j+1)^2-(j2-j3)^2)*((j2+j3+1)^2-(j+1)^2)*((j+1)^2-(m2+m3)^2)) - Y = (2*j+1)*((m2+m3)*(j2*(j2+1)-j3*(j3+1)) - (m2-m3)*j*(j+1)) - Z = (j+1)*sqrt((j^2-(j2-j3)^2)*((j2+j3+1)^2-j^2)*(j^2-(m2+m3)^2)) - 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 - -# test 6j -#---------- -# test orthogonality -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)) - @test length(j6range) == length(j3range) - M = zeros(length(j3range),length(j6range)) - for (k2,j6) in enumerate(j6range) - for (k1,j3) in enumerate(j3range) - M[k1,k2] = sqrt(2*j3+1)*sqrt(2*j6+1)*wigner6j(j1,j2,j3,j4,j5,j6) +@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) + ind1 = 1 + for m1 in -j1:j1, m2 in -j2:j2 + ind2 = 1 + for j3 in abs(j1-j2):(j1+j2) + for m3 in -j3:j3 + M[ind1,ind2] = clebschgordan(j1,m1,j2,m2,j3,m3) + ind2 += 1 + end end + ind1 += 1 end @test M'*M ≈ one(M) end end -# test special case -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 - # test recurrence relations: Phys Rev E 57, 7274 (1998) -for k = 1:10 - j2 = convert(BigFloat,rand(largejlist)) - j3 = convert(BigFloat,rand(largejlist)) - l1 = convert(BigFloat,rand(largejlist)) - l2 = convert(BigFloat,rand(abs(l1-j3):(l1+j3))) - l3 = convert(BigFloat,rand(abs(l1-j2):min(l1+j2))) +@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)) + m3 = convert(BigFloat,rand(-j3:j3)) - for j in intersect(abs(j2-j3):(j2+j3), abs(l2-l3):(l2+l3)) - X = j*sqrt(((j+1)^2-(j2-j3)^2)*((j2+j3+1)^2-(j+1)^2)*((j+1)^2-(l2-l3)^2)*((l2+l3+1)^2 - (j+1)^2)) - Y = (2*j+1)*( j*(j+1)*( -j*(j+1) + j2*(j2+1) + j3*(j3+1) - 2*l1*(l1+1)) + - l2*(l2+1)*( j*(j+1) + j2*(j2+1) - j3*(j3+1) ) + - l3*(l3+1)*( j*(j+1) - j2*(j2+1) + j3*(j3+1) ) ) - Z = (j+1)*sqrt((j^2-(j2-j3)^2)*((j2+j3+1)^2-j^2)*(j^2-(l2-l3)^2)*((l2+l3+1)^2 - j^2)) - 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 - -# 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 - m1range = -j1:j1 - m2range = -j2:j2 - m3range = -j3:j3 - V1 = Array{Float64}(length(m1range),length(m2range),length(m3range)) - V2 = Array{Float64}(length(m1range),length(m2range),length(m3range)) - for J in max(abs(j1-j2-j3),abs(j1-j2+j3),abs(j1+j2-j3)):(j1+j2+j3) - J12range = max(abs(j1-j2),abs(J-j3)):min((j1+j2),(J+j3)) - J23range = max(abs(j2-j3),abs(j1-J)):min((j2+j3),(j1+J)) - for J12 in J12range, J23 in J23range - M = J # only test for J, should be independent - fill!(V1,0) - fill!(V2,0) - for (k1,m1) in enumerate(m1range) - for (k2,m2) in enumerate(m2range) - abs(m1+m2)<=J12 || continue - for (k3,m3) in enumerate(m3range) - abs(m2+m3)<=J23 || continue - m1+m2+m3==M || continue - V1[k1,k2,k3] = clebschgordan(j1,m1,j2,m2,J12)*clebschgordan(J12,m1+m2,j3,m3,J) - V2[k1,k2,k3] = clebschgordan(j2,m2,j3,m3,J23)*clebschgordan(j1,m1,J23,m2+m3,J) - end - end - end - @test racahW(j1,j2,J,j3,J12,J23) ≈ vecdot(V2,V1)/sqrt((2*J12+1)*(2*J23+1)) atol=10*eps(Float64) + for j in max(abs(j2-j3),abs(m2+m3))+1:(j2+j3)-1 + X = j*sqrt(((j+1)^2-(j2-j3)^2)*((j2+j3+1)^2-(j+1)^2)*((j+1)^2-(m2+m3)^2)) + Y = (2*j+1)*((m2+m3)*(j2*(j2+1)-j3*(j3+1)) - (m2-m3)*j*(j+1)) + Z = (j+1)*sqrt((j^2-(j2-j3)^2)*((j2+j3+1)^2-j^2)*(j^2-(m2+m3)^2)) + 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 +#---------- +@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)) + @test length(j6range) == length(j3range) + M = zeros(length(j3range),length(j6range)) + for (k2,j6) in enumerate(j6range) + for (k1,j3) in enumerate(j3range) + M[k1,k2] = sqrt(2*j3+1)*sqrt(2*j6+1)*wigner6j(j1,j2,j3,j4,j5,j6) + end + end + @test M'*M ≈ one(M) + end + end +end + +@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 + +@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)) + l2 = convert(BigFloat,rand(abs(l1-j3):(l1+j3))) + l3 = convert(BigFloat,rand(abs(l1-j2):min(l1+j2))) + + for j in intersect(abs(j2-j3):(j2+j3), abs(l2-l3):(l2+l3)) + X = j*sqrt(((j+1)^2-(j2-j3)^2)*((j2+j3+1)^2-(j+1)^2)*((j+1)^2-(l2-l3)^2)*((l2+l3+1)^2 - (j+1)^2)) + Y = (2*j+1)*( j*(j+1)*( -j*(j+1) + j2*(j2+1) + j3*(j3+1) - 2*l1*(l1+1)) + + l2*(l2+1)*( j*(j+1) + j2*(j2+1) - j3*(j3+1) ) + + l3*(l3+1)*( j*(j+1) - j2*(j2+1) + j3*(j3+1) ) ) + Z = (j+1)*sqrt((j^2-(j2-j3)^2)*((j2+j3+1)^2-j^2)*(j^2-(l2-l3)^2)*((l2+l3+1)^2 - j^2)) + 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 + +@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 + V1 = Array{Float64}(length(m1range),length(m2range),length(m3range)) + V2 = Array{Float64}(length(m1range),length(m2range),length(m3range)) + for J in max(abs(j1-j2-j3),abs(j1-j2+j3),abs(j1+j2-j3)):(j1+j2+j3) + J12range = max(abs(j1-j2),abs(J-j3)):min((j1+j2),(J+j3)) + J23range = max(abs(j2-j3),abs(j1-J)):min((j2+j3),(j1+J)) + for J12 in J12range, J23 in J23range + M = J # only test for J, should be independent + fill!(V1,0) + fill!(V2,0) + for (k1,m1) in enumerate(m1range) + for (k2,m2) in enumerate(m2range) + abs(m1+m2)<=J12 || continue + for (k3,m3) in enumerate(m3range) + abs(m2+m3)<=J23 || continue + m1+m2+m3==M || continue + V1[k1,k2,k3] = clebschgordan(j1,m1,j2,m2,J12)*clebschgordan(J12,m1+m2,j3,m3,J) + V2[k1,k2,k3] = clebschgordan(j2,m2,j3,m3,J23)*clebschgordan(j1,m1,J23,m2+m3,J) + end + end + end + @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