@@ -6,6 +6,10 @@ const BACKWARD = false
66const sqrtpi = 1.772453850905516027298
77const edivsqrt2pi = 1.084437551419227546612
88
9+ half (x:: Number ) = oftype (x,0.5 )
10+ half {T<:Number} (:: Type{T} ) = convert (T,0.5 )
11+ two (x:: Number ) = oftype (x,2 )
12+ two {T<:Number} (:: Type{T} ) = convert (T,2 )
913
1014"""
1115Pochhammer symbol (x)_n = Γ(x+n)/Γ(x) for the rising factorial.
@@ -120,7 +124,7 @@ Anαβ{T<:Integer}(n::AbstractMatrix{T},α::Number,β::Number) = [ Anαβ(n[i,j]
120124#
121125# I. Bogaert and B. Michiels and J. Fostier, 𝒪(1) computation of Legendre polynomials and Gauss--Legendre nodes and weights for parallel computing, SIAM J. Sci. Comput., 34:C83--C101, 2012.
122126#
123- Cx (x:: Number ) = exp (lgamma (x+ 1 / 2 ) - lgamma (x+ 1 ))
127+ Cx (x:: Number ) = exp (lgamma (x+ half (x)) - lgamma (x+ one (x) ))
124128function Cx (x:: Float64 )
125129 if x > 9.84475
126130 xp = x+ 0.25
@@ -132,12 +136,12 @@ end
132136@vectorize_1arg Number Cx
133137
134138Cnλ (n:: Integer ,λ:: Float64 ) = 2 ^ λ/ sqrtpi* Cx (n+ λ)
135- Cnλ (n:: Integer ,λ:: Number ) = 2 ^ λ/ sqrt (convert ( typeof (λ) ,π))* Cx (n+ λ)
139+ Cnλ (n:: Integer ,λ:: Number ) = 2 ^ λ/ sqrt (oftype (λ ,π))* Cx (n+ λ)
136140function Cnλ {T<:Integer} (n:: UnitRange{T} ,λ:: Number )
137141 ret = Vector {typeof(λ)} (length (n))
138142 ret[1 ] = Cnλ (first (n),λ)
139143 for i= 2 : length (n)
140- ret[i] = (n[i]+ λ- 1 / 2 )/ (n[i]+ λ)* ret[i- 1 ]
144+ ret[i] = (n[i]+ λ- half (λ) )/ (n[i]+ λ)* ret[i- 1 ]
141145 end
142146 ret
143147end
199203function absf (α:: Number ,β:: Number ,m:: Int ,θ:: Number )
200204 ret = zero (θ)
201205 for l= 0 : m
202- ret += pochhammer (1 / 2 + α,l)* pochhammer (1 / 2 - α,l)* pochhammer (1 / 2 + β,m- l)* pochhammer (1 / 2 - β,m- l)/ factorial (l)/ factorial (m- l)/ sinpi (θ/ 2 )^ (l+ α+ 1 / 2 ) / cospi (θ/ 2 )^ (m- l+ β+ 1 / 2 )
206+ ret += pochhammer (half (α) + α,l)* pochhammer (half (α) - α,l)* pochhammer (half (β) + β,m- l)* pochhammer (half (β) - β,m- l)/ factorial (l)/ factorial (m- l)/ sinpi (θ/ 2 )^ (l+ α+ half (α)) / cospi (θ/ 2 )^ (m- l+ β+ half (β) )
203207 end
204208 ret
205209end
@@ -213,10 +217,10 @@ function absf{T<:Number}(α::Number,β::Number,m::Int,θ::AbstractArray{T,1})
213217 ret = zero (θ)
214218 cfs = zeros (T,m+ 1 )
215219 for l= 0 : m
216- @inbounds cfs[l+ 1 ] = pochhammer (1 / 2 + α,l)* pochhammer (1 / 2 - α,l)* pochhammer (1 / 2 + β,m- l)* pochhammer (1 / 2 - β,m- l)/ factorial (l)/ factorial (m- l)
220+ @inbounds cfs[l+ 1 ] = pochhammer (half (α) + α,l)* pochhammer (half (α) - α,l)* pochhammer (half (β) + β,m- l)* pochhammer (half (β) - β,m- l)/ factorial (l)/ factorial (m- l)
217221 end
218222 @inbounds for i= 1 : length (θ),l= 0 : m
219- ret[i] += cfs[l+ 1 ]/ sinpi (θ[i]/ 2 )^ (l+ α+ 1 / 2 ) / cospi (θ[i]/ 2 )^ (m- l+ β+ 1 / 2 )
223+ ret[i] += cfs[l+ 1 ]/ sinpi (θ[i]/ 2 )^ (l+ α+ half (α)) / cospi (θ[i]/ 2 )^ (m- l+ β+ half (β) )
220224 end
221225 ret
222226end
@@ -245,12 +249,12 @@ end
245249function compute_umvm! {T<:AbstractFloat} (um:: Vector{T} ,vm:: Vector{T} ,cfs:: Matrix{T} ,α:: T ,β:: T ,tempcos:: Vector{T} ,tempsin:: Vector{T} ,tempcosβsinα:: Vector{T} ,m:: Int ,θ:: Vector{T} ,ir:: UnitRange{Int64} )
246250 @inbounds for i in ir
247251 temp = inv (tempcos[i]^ m* tempcosβsinα[i])
248- ϑ = (α+ 1 / 2 )/ 2 - (α+ β+ m+ 1 )* θ[i]/ 2
252+ ϑ = (α+ half (α) )/ 2 - (α+ β+ m+ 1 )* θ[i]/ 2
249253 um[i] = cfs[m+ 1 ,1 ]* cospi (ϑ)* temp
250254 vm[i] = cfs[m+ 1 ,1 ]* sinpi (ϑ)* temp
251255 @inbounds for l= 1 : m
252256 temp *= tempcos[i]/ tempsin[i]
253- ϑ = (α+ l+ 1 / 2 )/ 2 - (α+ β+ m+ 1 )* θ[i]/ 2
257+ ϑ = (α+ l+ half (α) )/ 2 - (α+ β+ m+ 1 )* θ[i]/ 2
254258 um[i] += cfs[m+ 1 ,l+ 1 ]* cospi (ϑ)* temp
255259 vm[i] += cfs[m+ 1 ,l+ 1 ]* sinpi (ϑ)* temp
256260 end
260264function compute_umvm! {T<:AbstractFloat} (um:: Vector{T} ,vm:: Vector{T} ,λ:: T ,tempsin:: Vector{T} ,tempsinλ:: Vector{T} ,m:: Int ,θ:: Vector{T} ,ir:: UnitRange{Int64} )
261265 @inbounds for i in ir
262266 temp = inv (tempsin[i]^ m* tempsinλ[i])
263- ϑ = (m+ λ)* (1 / 2 - θ[i])
267+ ϑ = (m+ λ)* (half (T) - θ[i])
264268 um[i] = cospi (ϑ)* temp
265269 vm[i] = sinpi (ϑ)* temp
266270 end
325329function init_cfs {T<:AbstractFloat} (α:: T ,β:: T ,M:: Int )
326330 cfs = zeros (T,M+ 1 ,M+ 1 )
327331 @inbounds for m= 0 : M,l= 0 : m
328- cfs[m+ 1 ,l+ 1 ] = pochhammer (1 / 2 + α,l)* pochhammer (1 / 2 - α,l)* pochhammer (1 / 2 + β,m- l)* pochhammer (1 / 2 - β,m- l)/ factorial (l)/ factorial (m- l)
332+ cfs[m+ 1 ,l+ 1 ] = pochhammer (half (α) + α,l)* pochhammer (half (α) - α,l)* pochhammer (half (β) + β,m- l)* pochhammer (half (β) - β,m- l)/ factorial (l)/ factorial (m- l)
329333 end
330334 cfs
331335end
0 commit comments