new Gaunt method for 32-bit word sizes

Author: MikaelSlevinsky

src/gaunt.jl

@@ -14,7 +14,7 @@ This is a Julia implementation of the stable recurrence described in:
 
 Y.-l. Xu, Fast evaluation of Gaunt coefficients: recursive approach, *J. Comp. Appl. Math.*, **85**:53–65, 1997.
 """
-function gaunt{T}(::Type{T},m::Int,n::Int,μ::Int,ν::Int;normalized::Bool=false)
+function gaunt{T}(::Type{T},m::Integer,n::Integer,μ::Integer,ν::Integer;normalized::Bool=false)
     if normalized
         normalizedgaunt(T,m,n,μ,ν)
     else
@@ -24,15 +24,17 @@ end
 doc"""
 Calculates the Gaunt coefficients in 64-bit floating-point arithmetic.
 """
-gaunt(m::Int,n::Int,μ::Int,ν::Int;kwds...) = gaunt(Float64,m,n,μ,ν;kwds...)
+gaunt(m::Integer,n::Integer,μ::Integer,ν::Integer;kwds...) = gaunt(Float64,m,n,μ,ν;kwds...)
 
-function normalization{T}(::Type{T},m::Int,n::Int,μ::Int,ν::Int)
+gaunt{T}(::Type{T},m::Int32,n::Int32,μ::Int32,ν::Int32;normalized::Bool=false) = gaunt(T,Int64(m),Int64(n),Int64(μ),Int64(ν);normalized=normalized)
+
+function normalization{T}(::Type{T},m::Integer,n::Integer,μ::Integer,ν::Integer)
     pochhammer(n+one(T),n)*pochhammer(ν+one(T),ν)/pochhammer(n+ν+one(T),n+ν)*gamma(n+ν-m-μ+one(T))/gamma(n-m+one(T))/gamma(ν-μ+one(T))
 end
 
-normalization(::Type{Float64},m::Int,n::Int,μ::Int,ν::Int) = normalization1(Float64,n,ν)*normalization2(Float64,n-m,ν-μ)
+normalization(::Type{Float64},m::Integer,n::Integer,μ::Integer,ν::Integer) = normalization1(Float64,n,ν)*normalization2(Float64,n-m,ν-μ)
 
-function normalization1(::Type{Float64},n::Int,ν::Int)
+function normalization1(::Type{Float64},n::Integer,ν::Integer)
     if n ≥ 8
         if ν ≥ 8
             return exp((n+0.5)*log1p(n/(n+1))+(ν+0.5)*log1p(ν/(ν+1))+(n+ν+0.5)*log1p(-(n+ν)/(2n+2ν+1))+n*log1p(-2ν/(2n+2ν+1))+ν*log1p(-2n/(2n+2ν+1)))*stirlingseries(2n+1.0)*stirlingseries(2ν+1.0)*stirlingseries(n+ν+1.0)/stirlingseries(n+1.0)/stirlingseries(ν+1.0)/stirlingseries(2n+2ν+1.0)
@@ -46,7 +48,7 @@ function normalization1(::Type{Float64},n::Int,ν::Int)
     end
 end
 
-function normalization2(::Type{Float64},nm::Int,νμ::Int)
+function normalization2(::Type{Float64},nm::Integer,νμ::Integer)
     if nm ≥ 8
         if νμ ≥ 8
             return edivsqrt2pi*exp((nm+0.5)*log1p(νμ/(nm+1))+(νμ+0.5)*log1p(nm/(νμ+1)))/sqrt(nm+νμ+1.0)*stirlingseries(nm+νμ+1.0)/stirlingseries(nm+1.0)/stirlingseries(νμ+1.0)
@@ -60,7 +62,7 @@ function normalization2(::Type{Float64},nm::Int,νμ::Int)
     end
 end
 
-function normalizedgaunt{T}(::Type{T},m::Int,n::Int,μ::Int,ν::Int)
+function normalizedgaunt{T}(::Type{T},m::Integer,n::Integer,μ::Integer,ν::Integer)
     qmax = min(n,ν,(n+ν-abs(m+μ))÷2)
     a = Vector{T}(qmax+1)
     a[1] = one(T)
@@ -106,15 +108,15 @@ function normalizedgaunt{T}(::Type{T},m::Int,n::Int,μ::Int,ν::Int)
     a
 end
 
-function secondinitialcondition{T}(::Type{T},m::Int,n::Int,μ::Int,ν::Int)
+function secondinitialcondition{T}(::Type{T},m::Integer,n::Integer,μ::Integer,ν::Integer)
     n₄ = n+ν-m-μ
     mn = m-n
     μν = μ-ν
     temp = 2n+2ν-one(T)
     return (temp-2)/2*(1-temp/n₄/(n₄-1)*(mn*(mn+one(T))/(2n-1)+μν*(μν+one(T))/(2ν-1)))
 end
 
-function thirdinitialcondition{T}(::Type{T},m::Int,n::Int,μ::Int,ν::Int)
+function thirdinitialcondition{T}(::Type{T},m::Integer,n::Integer,μ::Integer,ν::Integer)
     n₄ = n+ν-m-μ
     mn = m-n
     μν = μ-ν
@@ -124,14 +126,14 @@ function thirdinitialcondition{T}(::Type{T},m::Int,n::Int,μ::Int,ν::Int)
     return temp*(temp-6)/4*( (temp-2)/n₄/(n₄-1)*temp2 + one(T)/2 )
 end
 
-α{T}(::Type{T},n::Int,ν::Int,p::Int) = (p^2-(n+ν+1)^2)*(p^2-(n-ν)^2)/(4p^2-one(T))
-A{T}(::Type{T},m::Int,n::Int,μ::Int,ν::Int,p::Int) = p*(p-one(T))*(m-μ)-(m+μ)*(n-ν)*(n+ν+one(T))
+α{T}(::Type{T},n::Integer,ν::Integer,p::Integer) = (p^2-(n+ν+1)^2)*(p^2-(n-ν)^2)/(4p^2-one(T))
+A{T}(::Type{T},m::Integer,n::Integer,μ::Integer,ν::Integer,p::Integer) = p*(p-one(T))*(m-μ)-(m+μ)*(n-ν)*(n+ν+one(T))
 
-c₀{T}(::Type{T},m::Int,n::Int,μ::Int,ν::Int,p::Int,p₁::Int) = (p+2)*(p+3)*(p₁+1)*(p₁+2)*A(T,m,n,μ,ν,p+4)*α(T,n,ν,p+1)
-c₁{T}(::Type{T},m::Int,n::Int,μ::Int,ν::Int,p::Int,p₁::Int,p₂::Int) = A(T,m,n,μ,ν,p+2)*A(T,m,n,μ,ν,p+3)*A(T,m,n,μ,ν,p+4) + (p+1)*(p+3)*(p₁+2)*(p₂+2)*A(T,m,n,μ,ν,p+4)*α(T,n,ν,p+2) + (p+2)*(p+4)*(p₁+3)*(p₂+3)*A(T,m,n,μ,ν,p+2)*α(T,n,ν,p+3)
-c₂{T}(::Type{T},m::Int,n::Int,μ::Int,ν::Int,p::Int,p₂::Int) = -(p+2)*(p+3)*(p₂+3)*(p₂+4)*A(T,m,n,μ,ν,p+2)*α(T,n,ν,p+4)
+c₀{T}(::Type{T},m::Integer,n::Integer,μ::Integer,ν::Integer,p::Integer,p₁::Integer) = (p+2)*(p+3)*(p₁+1)*(p₁+2)*A(T,m,n,μ,ν,p+4)*α(T,n,ν,p+1)
+c₁{T}(::Type{T},m::Integer,n::Integer,μ::Integer,ν::Integer,p::Integer,p₁::Integer,p₂::Integer) = A(T,m,n,μ,ν,p+2)*A(T,m,n,μ,ν,p+3)*A(T,m,n,μ,ν,p+4) + (p+1)*(p+3)*(p₁+2)*(p₂+2)*A(T,m,n,μ,ν,p+4)*α(T,n,ν,p+2) + (p+2)*(p+4)*(p₁+3)*(p₂+3)*A(T,m,n,μ,ν,p+2)*α(T,n,ν,p+3)
+c₂{T}(::Type{T},m::Integer,n::Integer,μ::Integer,ν::Integer,p::Integer,p₂::Integer) = -(p+2)*(p+3)*(p₂+3)*(p₂+4)*A(T,m,n,μ,ν,p+2)*α(T,n,ν,p+4)
 
-d₀{T}(::Type{T},m::Int,n::Int,μ::Int,ν::Int,p::Int,p₁::Int) = (p+2)*(p+3)*(p+5)*(p₁+2)*(p₁+4)*A(T,m,n,μ,ν,p+6)*α(T,n,ν,p+1)
-d₁{T}(::Type{T},m::Int,n::Int,μ::Int,ν::Int,p::Int,p₁::Int,p₂::Int) = (p+5)*(p₁+4)*A(T,m,n,μ,ν,p+6)*( A(T,m,n,μ,ν,p+2)*A(T,m,n,μ,ν,p+3) + (p+1)*(p+3)*(p₁+2)*(p₂+2)*α(T,n,ν,p+2) )
-d₂{T}(::Type{T},m::Int,n::Int,μ::Int,ν::Int,p::Int,p₁::Int,p₂::Int) = (p+2)*(p₂+3)*A(T,m,n,μ,ν,p+2)*( A(T,m,n,μ,ν,p+5)*A(T,m,n,μ,ν,p+6) + (p+4)*(p+6)*(p₁+5)*(p₂+5)*α(T,n,ν,p+5) )
-d₃{T}(::Type{T},m::Int,n::Int,μ::Int,ν::Int,p::Int,p₂::Int) = -(p+2)*(p+4)*(p+5)*(p₂+3)*(p₂+5)*(p₂+6)*A(T,m,n,μ,ν,p+2)*α(T,n,ν,p+6)
+d₀{T}(::Type{T},m::Integer,n::Integer,μ::Integer,ν::Integer,p::Integer,p₁::Integer) = (p+2)*(p+3)*(p+5)*(p₁+2)*(p₁+4)*A(T,m,n,μ,ν,p+6)*α(T,n,ν,p+1)
+d₁{T}(::Type{T},m::Integer,n::Integer,μ::Integer,ν::Integer,p::Integer,p₁::Integer,p₂::Integer) = (p+5)*(p₁+4)*A(T,m,n,μ,ν,p+6)*( A(T,m,n,μ,ν,p+2)*A(T,m,n,μ,ν,p+3) + (p+1)*(p+3)*(p₁+2)*(p₂+2)*α(T,n,ν,p+2) )
+d₂{T}(::Type{T},m::Integer,n::Integer,μ::Integer,ν::Integer,p::Integer,p₁::Integer,p₂::Integer) = (p+2)*(p₂+3)*A(T,m,n,μ,ν,p+2)*( A(T,m,n,μ,ν,p+5)*A(T,m,n,μ,ν,p+6) + (p+4)*(p+6)*(p₁+5)*(p₂+5)*α(T,n,ν,p+5) )
+d₃{T}(::Type{T},m::Integer,n::Integer,μ::Integer,ν::Integer,p::Integer,p₂::Integer) = -(p+2)*(p+4)*(p+5)*(p₂+3)*(p₂+5)*(p₂+6)*A(T,m,n,μ,ν,p+2)*α(T,n,ν,p+6)

test/gaunttests.jl

@@ -3,21 +3,19 @@ using FastTransforms, Base.Test
 import FastTransforms: δ
 
 @testset "Gaunt coefficients" begin
-    if Sys.WORD_SIZE == 64
-        println("Testing Table 2 of Y.-l. Xu, JCAM 85:53–65, 1997.")
-        for (m,n) in ((0,2),(1,2),(1,8),(6,8),(3,18),
-                      (10,18),(5,25),(-23,25),(2,40),(-35,40),
-                      (28,62),(-42,62),(1,99),(90,99),(10,120),
-                      (80,120),(23,150),(88,150))
-            @test norm(gaunt(m,n,-m,n)[end]./(big(-1.0)^m/(2n+1))-1, Inf) < 400eps()
-        end
-        println("Testing Table 3 of Y.-l. Xu, JCAM 85:53–65, 1997.")
-        for (m,n,μ,ν) in ((0,1,0,5),(0,5,0,10),(0,9,0,10),(0,10,0,12),
-                          (0,11,0,15),(0,12,0,20),(0,20,0,45),(0,40,0,80),
-                          (0,45,0,100),(3,5,-3,6),(4,9,-4,15),(-8,18,8,23),
-                          (-10,20,10,30),(5,25,-5,45),(15,50,-15,60),(-28,68,28,75),
-                          (32,78,-32,88),(45,82,-45,100))
-            @test norm(sum(gaunt(m,n,μ,ν))-δ(m,0), Inf) < 15000eps()
-        end
+    println("Testing Table 2 of Y.-l. Xu, JCAM 85:53–65, 1997.")
+    for (m,n) in ((0,2),(1,2),(1,8),(6,8),(3,18),
+                  (10,18),(5,25),(-23,25),(2,40),(-35,40),
+                  (28,62),(-42,62),(1,99),(90,99),(10,120),
+                  (80,120),(23,150),(88,150))
+        @test norm(gaunt(m,n,-m,n)[end]./(big(-1.0)^m/(2n+1))-1, Inf) < 400eps()
+    end
+    println("Testing Table 3 of Y.-l. Xu, JCAM 85:53–65, 1997.")
+    for (m,n,μ,ν) in ((0,1,0,5),(0,5,0,10),(0,9,0,10),(0,10,0,12),
+                      (0,11,0,15),(0,12,0,20),(0,20,0,45),(0,40,0,80),
+                      (0,45,0,100),(3,5,-3,6),(4,9,-4,15),(-8,18,8,23),
+                      (-10,20,10,30),(5,25,-5,45),(15,50,-15,60),(-28,68,28,75),
+                      (32,78,-32,88),(45,82,-45,100))
+        @test norm(sum(gaunt(m,n,μ,ν))-δ(m,0), Inf) < 15000eps()
     end
 end