work on special functions

JuliaApproximation · Aug 26, 2017 · 12a498c · 12a498c
1 parent 5251e69
commit 12a498c
Show file tree

Hide file tree

Showing 5 changed files with 17 additions and 7 deletions.
diff --git a/src/HypergeometricFunctions/HypergeometricFunctions.jl b/src/HypergeometricFunctions/HypergeometricFunctions.jl
@@ -5,10 +5,10 @@ This module calculates (generalized) hypergeometric functions:
 """
 module HypergeometricFunctions
     using DualNumbers
-    import ApproxFun: @clenshaw, real, eps, reverseorientation
-    import FastTransforms: pochhammer
-    import SingularIntegralEquations: undirected
-    
+import ApproxFun: @clenshaw, real, eps, reverseorientation
+import FastTransforms: pochhammer
+import SingularIntegralEquations: undirected, Directed
+
 export _₂F₁, _₃F₂, mFn
 
 include("Gauss.jl")

diff --git a/src/HypergeometricFunctions/specialfunctions.jl b/src/HypergeometricFunctions/specialfunctions.jl
@@ -44,6 +44,7 @@ expnlog1psinhatanhsqrt(n,x::Union{T,Dual{T}}) where {T<:Real} = x == 0 ? one(x)
 expm1nlog1p(n,x) = x == 0 ? one(x) : expm1(n*log1p(x))/(n*x)
 
 speciallog(x) = x == 0 ? one(x) : x > 0 ? (s = sqrt(x); 3(atanh(s)-s)/s^3) : (s = sqrt(-x); 3(s-atan(s))/s^3)
+speciallog(x::Directed) = (s = sqrt(-x); 3(s-atan(s))/s^3)
 function speciallog(x::Float64)
     if x > 0.2
         s = sqrt(x)

diff --git a/src/SingularIntegralEquations.jl b/src/SingularIntegralEquations.jl
@@ -94,8 +94,8 @@ end
 Base.log(x::Directed{true}) = log(-x.x) - π*im
 Base.log(x::Directed{false}) = log(-x.x) + π*im
 Base.log1p(x::Directed) = log(1+x)
-Base.sqrt(x::Directed{true}) = -im*sqrt(-x.x)
-Base.sqrt(x::Directed{false}) = im*sqrt(-x.x)
+Base.sqrt(x::Directed{true}) = x.x ≥ 0 ? sqrt(complex(x.x)) : -im*sqrt(-x.x)
+Base.sqrt(x::Directed{false}) = x.x ≥ 0 ? sqrt(complex(x.x)) : im*sqrt(-x.x)
 ^(x::Directed{true},a::Integer) = x.x^a
 ^(x::Directed{false},a::Integer) = x.x^a
 ^(x::Directed{true},a::Number) = exp(-a*π*im)*(-x.x)^a

diff --git a/src/stieltjesmoment.jl b/src/stieltjesmoment.jl
@@ -21,7 +21,7 @@ stieltjesmoment(S::JacobiQ,z) = stieltjesjacobimoment(S.a,S.b,z)
 
 normalization(n::Int,α::Real,β::Real) = 2^(α+β)*gamma(n+α+1)*gamma(n+β+1)/gamma(2n+α+β+2)
 stieltjesjacobimoment(α::Real,β::Real,n::Int,z) =
-    (x = 2./(1-z);normalization(n,α,β)*(-x)^(n+1)*_₂F₁(n+1,n+α+1,2n+α+β+2,x))
+    (x = 2./(1-z);normalization(n,α,β)*undirected(-x)^(n+1)*_₂F₁(n+1,n+α+1,2n+α+β+2,x))
 stieltjesjacobimoment(α::Real,β::Real,z) = stieltjesjacobimoment(α,β,0,z)
 
 

diff --git a/test/runtests.jl b/test/runtests.jl
@@ -35,6 +35,9 @@ x=Fun()
 @test sum(logabs(x-2.0)) ≈ logabslegendremoment(2.0)
 @test sum(logabs(x-(2.0+im))) ≈ logabslegendremoment(2.0+im)
 
+@test sqrt(Directed{true}(-0.1)) ≈ sqrt(-0.1-0.0im)
+@test sqrt(Directed{false}(-0.1)) ≈ sqrt(-0.1+0.0im)
+
 @test stieltjesmoment(Legendre(),0,0.1+0im) ≈ stieltjesmoment(Legendre(),0,Directed{true}(0.1))
 @test stieltjesmoment(Legendre(),0,0.1-0im) ≈ stieltjesmoment(Legendre(),0,Directed{false}(0.1))
 @test stieltjesjacobimoment(0.,0.,1,0.1+0im) ≈ stieltjesjacobimoment(0.,0.,1,0.1*⁺)
@@ -43,6 +46,12 @@ x=Fun()
 @test (stieltjesjacobimoment(0,0,0,0.1+0im)-stieltjesjacobimoment(0,0,0,0.1-0im))/(-2π*im) ≈ 1.0
 @test (stieltjesjacobimoment(0,0,1,0.1+0im)-stieltjesjacobimoment(0,0,1,0.1-0im))/(-2π*im) ≈ 0.1
 
+@test stieltjesjacobimoment(0.5,0,0,Directed{true}(0.1)) ≈ stieltjesjacobimoment(0.5,0,0,0.1+0im)
+@test stieltjesjacobimoment(0.5,0,0,Directed{false}(0.1)) ≈ stieltjesjacobimoment(0.5,0,0,0.1-0.0im)
+
+
+
+
 println("Hilbert test")
 include("HilbertTest.jl")