Doublefloats and hypergeometric functions (#92)

* doublefloats and lower incomplete gamma function

* conversion of doublefloats to arbfloat and arbreal

Author: chillenb

Project.toml

@@ -17,6 +17,7 @@ Readables = "0d4725de-cd7c-5e44-8a85-a48caeef9fa5"
 GMP_jll = "781609d7-10c4-51f6-84f2-b8444358ff6d"
 MPFR_jll = "3a97d323-0669-5f0c-9066-3539efd106a3"
 FLINT_jll = "e134572f-a0d5-539d-bddf-3cad8db41a82"
+DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
 
 [compat]
 FLINT_jll = "301.300 - 304.700"

src/ArbNumerics.jl

@@ -62,6 +62,7 @@ export ArbNumber,
        hypergeometric_0F1, hypergeometric_1F1, hypergeometric_2F1,
        regular_hypergeometric_0F1, regular_hypergeometric_1F1, regular_hypergeometric_2F1,
        F₀₁,  F₁₁, F₂₁, regularF₀₁,  regularF₁₁, regularF₂₁,
+       hypergeometric_U, hypergeometric_gamma_lower, regular_hypergeometric_gamma_lower, further_regular_hypergeometric_gamma_lower,
        elliptic_k, elliptic_e, elliptic_pi, elliptic_f,
        elliptic_k2, elliptic_e2, elliptic_pi2, elliptic_f2, # modulus^2
        elliptic_rf, elliptic_rg, elliptic_rj,
@@ -127,6 +128,8 @@ import SpecialFunctions: gamma, lgamma, lfact, digamma, invdigamma, polygamma, t
      besselj, besselj0, besselj1, bessely, bessely0, bessely1, besseli, besselk,
      eta, zeta
 
+import DoubleFloats: DoubleFloat
+
 using FLINT_jll
 
 macro libarb(function_name)

src/float/hypergeometric.jl

@@ -169,25 +169,76 @@ function regular_hypergeometric_2F1(a::ArbReal{P}, b::ArbReal{P}, c::ArbReal{P},
     return result
 end
 
+"""
+    hypergeometric_gamma_lower(s, z)
+
+lower incomplete gamma function, ``small_gamma(s, z) = (z^s / s) * ₁F₁(s, s+1, -z)``
+"""
+function hypergeometric_gamma_lower(s::ArbReal{P}, z::ArbReal{P}) where {P}
+    result = ArbReal{P}()
+    regularized = Cint(0)
+    ccall(@libarb(arb_hypgeom_gamma_lower), Cvoid, (Ref{ArbReal}, Ref{ArbReal}, Ref{ArbReal}, Cint, Clong),
+         result, s, z, regularized, P)
+    return result
+end
+
+"""
+    regular_hypergeometric_gamma_lower(s, z)
+
+regularized lower incomplete gamma function, ``P(s, z) = small_gamma(s, z) / gamma(s)``
+"""
+function regular_hypergeometric_gamma_lower(s::ArbReal{P}, z::ArbReal{P}) where {P}
+    result = ArbReal{P}()
+    regularized = Cint(1)
+    ccall(@libarb(arb_hypgeom_gamma_lower), Cvoid, (Ref{ArbReal}, Ref{ArbReal}, Ref{ArbReal}, Cint, Clong),
+         result, s, z, regularized, P)
+    return result
+end
+
+"""
+    further_regular_hypergeometric_gamma_lower(s, z)
+
+'further' regularized lower incomplete gamma function, ``small_gamma_star(s, z) = z^-s * P(s, z)``
+"""
+function further_regular_hypergeometric_gamma_lower(s::ArbReal{P}, z::ArbReal{P}) where {P}
+    result = ArbReal{P}()
+    regularized = Cint(2)
+    ccall(@libarb(arb_hypgeom_gamma_lower), Cvoid, (Ref{ArbReal}, Ref{ArbReal}, Ref{ArbReal}, Cint, Clong),
+         result, s, z, regularized, P)
+    return result
+end
+
+
+
+
 # ArbFloat
 
 hypergeometric_0F1(a::ArbFloat{P}, z::ArbFloat{P}) where {P} =
-    ArbFloat{P}(hypgeom0f1(ArbReal{P}(a), ArbReal{P}(z)))
+    ArbFloat{P}(hypergeometric_0F1(ArbReal{P}(a), ArbReal{P}(z)))
 
 hypergeometric_0F1_regularized(a::ArbFloat{P}, z::ArbFloat{P}) where {P} =
-    ArbFloat{P}(hypgeom0f1reg(ArbReal{P}(a), ArbReal{P}(z)))
+    ArbFloat{P}(hypergeometric_0F1_regularized(ArbReal{P}(a), ArbReal{P}(z)))
 
-hypergeometric_U(a::ArbFloat{P}, z::ArbFloat{P}) where {P} =
-    ArbFloat{P}(hypgeomu(ArbReal{P}(a), ArbReal{P}(z)))
+hypergeometric_U(a::ArbFloat{P}, b::ArbFloat{P}, z::ArbFloat{P}) where {P} =
+    ArbFloat{P}(hypergeometric_U(ArbReal{P}(a), ArbReal{P}(b), ArbReal{P}(z)))
 
 hypergeometric_1F1(a::ArbFloat{P}, b::ArbFloat{P}, z::ArbFloat{P}) where {P} =
-    ArbFloat{P}(hypgeom1f1(ArbReal{P}(a), ArbReal{P}(b), ArbReal{P}(z)))
+    ArbFloat{P}(hypergeometric_1F1(ArbReal{P}(a), ArbReal{P}(b), ArbReal{P}(z)))
 
 hypergeometric_1F1_regularized(a::ArbFloat{P}, b::ArbFloat{P}, z::ArbFloat{P}) where {P} =
-    ArbFloat{P}(hypgeom1f1reg(ArbReal{P}(a), ArbReal{P}(b), ArbReal{P}(z)))
+    ArbFloat{P}(hypergeometric_1F1_regularized(ArbReal{P}(a), ArbReal{P}(b), ArbReal{P}(z)))
 
 hypergeometric_1F2(a::ArbFloat{P}, b::ArbFloat{P}, c::ArbFloat{P}, z::ArbFloat{P}) where {P} =
-    ArbFloat{P}(hypgeom1f2(ArbReal{P}(a), ArbReal{P}(b), ArbReal{P}(c), ArbReal{P}(z)))
+    ArbFloat{P}(hypergeometric_1F2(ArbReal{P}(a), ArbReal{P}(b), ArbReal{P}(c), ArbReal{P}(z)))
 
 hypergeometric_1F2_regularized(a::ArbFloat{P}, b::ArbFloat{P}, c::ArbFloat{P}, z::ArbFloat{P}) where {P} =
-    ArbFloat{P}(hypgeom1f2reg(ArbReal{P}(a), ArbReal{P}(b), ArbReal{P}(c), ArbReal{P}(z)))
+    ArbFloat{P}(hypergeometric_1F2_regularized(ArbReal{P}(a), ArbReal{P}(b), ArbReal{P}(c), ArbReal{P}(z)))
+
+hypergeometric_gamma_lower(s::ArbFloat{P}, z::ArbFloat{P}) where {P} =
+    ArbFloat{P}(hypergeometric_gamma_lower(ArbReal{P}(s), ArbReal{P}(z)))
+
+regular_hypergeometric_gamma_lower(s::ArbFloat{P}, z::ArbFloat{P}) where {P} =
+    ArbFloat{P}(hypergeometric_gamma_lower(ArbReal{P}(s), ArbReal{P}(z)))
+
+further_regular_hypergeometric_gamma_lower(s::ArbFloat{P}, z::ArbFloat{P}) where {P} =
+    ArbFloat{P}(further_hypergeometric_gamma_lower(ArbReal{P}(s), ArbReal{P}(z)))

src/libarb/ArbFloat.jl

@@ -107,6 +107,8 @@ end
 
 ArbFloat{P}(x::T) where {P,T<:Integer} = ArbFloat{P}(BigFloat(x, precision=P + 2))
 
+ArbFloat{P}(x::DoubleFloat{T}) where {P,T<:IEEEFloat} = ArbFloat{P}(x.hi) + ArbFloat{P}(x.lo)
+
 function ArbFloat{P}(x::Base.IEEEFloat) where {P}
     z = ArbFloat{P}()
     ccall(@libarb(arf_set_d), Cvoid, (Ref{ArbFloat}, Cdouble), z, x)

src/libarb/ArbReal.jl

@@ -60,6 +60,9 @@ function ArbReal{P}(x::USlong) where {P}
 end
 ArbReal{P}(x::T) where {P, T<:Integer} = ArbReal{P}(BigFloat(x, precision=P + 2))
 
+ArbReal{P}(x::DoubleFloat{T}) where {P,T<:IEEEFloat} = ArbReal{P}(x.hi) + ArbReal{P}(x.lo)
+
+
 function ArbReal{P}(x::Base.IEEEFloat) where {P}
     z = ArbReal{P}()
     ccall(@libarb(arb_set_d), Cvoid, (Ref{ArbReal}, Cdouble), z, x)