Feature: type-generic function interfaces (#4)

src/Cl.jl

@@ -173,10 +173,10 @@ function cl_series(n::Integer, x::Float64)
 end
 
 """
-    cl(n::Integer, x::Float64)::Float64
+    cl(n::Integer, x::Real)::Real
 
 Returns the value of the Clausen function ``\\operatorname{Cl}_n(x)``
-for integers ``n > 0`` and a real angle ``x`` of type `Float64`.  This
+for integers ``n > 0`` and a real angle ``x`` of type `Real`.  This
 function is defined as
 
 ```math
@@ -186,8 +186,8 @@ function is defined as
 \\end{aligned}
 ```
 
-Note: ``\\operatorname{Cl}_1(x)`` is not defined for ``x=2n\\pi`` with
-``n\\in\\mathbb{Z}``.
+Note: ``\\operatorname{Cl}_1(x)`` is not defined for ``x=2k\\pi`` with
+``k\\in\\mathbb{Z}``.
 
 The implementation follows the approach presented in [Jiming Wu,
 Xiaoping Zhang, Dongjie Liu, "An efficient calculation of the Clausen
@@ -204,7 +204,13 @@ julia> cl(10, 1.0)
 0.8423605391686301
 ```
 """
-function cl(n::Integer, x::Float64)::Float64
+cl(n::Integer, x::Real) = _cl(n, float(x))
+
+_cl(n::Integer, x::Float16) = oftype(x, _cl(n, Float32(x)))
+
+_cl(n::Integer, x::Float32) = oftype(x, _cl(n, Float64(x)))
+
+function _cl(n::Integer, x::Float64)::Float64
     n < 1  && throw(DomainError(n, "cl(n,x) undefined for n < 1"))
     n == 1 && return cl1(x)
     n == 2 && return cl2(x)

src/Cl1.jl

@@ -1,15 +1,15 @@
 """
-    cl1(x::Float64)::Float64
+    cl1(x::Real)::Real
 
 Returns the value of the Clausen function ``\\operatorname{Cl}_1(x)``
-for a real angle ``x`` of type `Float64`.  This function is defined as
+for a real angle ``x`` of type `Real`.  This function is defined as
 
 ```math
 \\operatorname{Cl}_1(x) = \\Re[\\operatorname{Li}_1(e^{ix})] = \\Re[-\\log(1 - e^{ix})]
 ```
 
-Note: ``\\operatorname{Cl}_1(x)`` is not defined for ``x=2n\\pi`` with
-``n\\in\\mathbb{Z}``.
+Note: ``\\operatorname{Cl}_1(x)`` is not defined for ``x=2k\\pi`` with
+``k\\in\\mathbb{Z}``.
 
 Author: Alexander Voigt
 
@@ -21,12 +21,12 @@ julia> cl1(1.0)
 0.04201950582536895
 ```
 """
-function cl1(x::Float64)::Float64
+function cl1(x::Real)::Real
     x = range_reduce_odd(x)
 
     if x == zero(x)
         return Inf
     end
 
-    -log(2.0*sin(0.5*x))
+    -log(2*sin(one(x)/2*x))
 end

src/Cl2.jl

@@ -1,8 +1,8 @@
 """
-    cl2(x::Float64)::Float64
+    cl2(x::Real)::Real
 
 Returns the value of the Clausen function ``\\operatorname{Cl}_2(x)``
-for a real angle ``x`` of type `Float64`.  The Clausen function is
+for a real angle ``x`` of type `Real`.  The Clausen function is
 defined as
 
 ```math
@@ -19,7 +19,13 @@ julia> cl2(1.0)
 1.0139591323607684
 ```
 """
-function cl2(x::Float64)::Float64
+cl2(x::Real) = _cl2(float(x))
+
+_cl2(x::Float16) = oftype(x, _cl2(Float32(x)))
+
+_cl2(x::Float32) = oftype(x, _cl2(Float64(x)))
+
+function _cl2(x::Float64)::Float64
     pi28 = pi*pi/8.0
 
     (x, sgn) = range_reduce_even(x)

src/Cl3.jl

@@ -1,8 +1,8 @@
 """
-    cl3(x::Float64)::Float64
+    cl3(x::Real)::Real
 
 Returns the value of the Clausen function ``\\operatorname{Cl}_3(x)``
-for a real angle ``x`` of type `Float64`.  This function is defined as
+for a real angle ``x`` of type `Real`.  This function is defined as
 
 ```math
 \\operatorname{Cl}_3(x) = \\Re[\\operatorname{Li}_3(e^{ix})] = \\sum_{k=1}^\\infty \\frac{\\cos(kx)}{k^3}
@@ -18,7 +18,13 @@ julia> cl3(1.0)
 0.44857300728001737
 ```
 """
-function cl3(x::Float64)::Float64
+cl3(x::Real) = _cl3(float(x))
+
+_cl3(x::Float16) = oftype(x, _cl3(Float32(x)))
+
+_cl3(x::Float32) = oftype(x, _cl3(Float64(x)))
+
+function _cl3(x::Float64)::Float64
     zeta3 = 1.2020569031595943
     pi28 = pi*pi/8.0
 

src/Cl4.jl

@@ -1,8 +1,8 @@
 """
-    cl4(x::Float64)::Float64
+    cl4(x::Real)::Real
 
 Returns the value of the Clausen function ``\\operatorname{Cl}_4(x)``
-for a real angle ``x`` of type `Float64`.  This function is defined as
+for a real angle ``x`` of type `Real`.  This function is defined as
 
 ```math
 \\operatorname{Cl}_4(x) = \\Im[\\operatorname{Li}_4(e^{ix})] = \\sum_{k=1}^\\infty \\frac{\\sin(kx)}{k^4}
@@ -18,7 +18,13 @@ julia> cl4(1.0)
 0.8958052386793799
 ```
 """
-function cl4(x::Float64)::Float64
+cl4(x::Real) = _cl4(float(x))
+
+_cl4(x::Float16) = oftype(x, _cl4(Float32(x)))
+
+_cl4(x::Float32) = oftype(x, _cl4(Float64(x)))
+
+function _cl4(x::Float64)::Float64
     zeta3 = 1.2020569031595943
     pi28 = pi*pi/8.0
 

src/Cl5.jl

@@ -1,8 +1,8 @@
 """
-    cl5(x::Float64)::Float64
+    cl5(x::Real)::Real
 
 Returns the value of the Clausen function ``\\operatorname{Cl}_5(x)``
-for a real angle ``x`` of type `Float64`.  This function is defined as
+for a real angle ``x`` of type `Real`.  This function is defined as
 
 ```math
 \\operatorname{Cl}_5(x) = \\Re[\\operatorname{Li}_5(e^{ix})] = \\sum_{k=1}^\\infty \\frac{\\cos(kx)}{k^5}
@@ -18,7 +18,13 @@ julia> cl5(1.0)
 0.5228208076420943
 ```
 """
-function cl5(x::Float64)::Float64
+cl5(x::Real) = _cl5(float(x))
+
+_cl5(x::Float16) = oftype(x, _cl5(Float32(x)))
+
+_cl5(x::Float32) = oftype(x, _cl5(Float64(x)))
+
+function _cl5(x::Float64)::Float64
     zeta5 = 1.0369277551433699
     pi28 = pi*pi/8.0
 

src/Cl6.jl

@@ -1,8 +1,8 @@
 """
-    cl6(x::Float64)::Float64
+    cl6(x::Real)::Real
 
 Returns the value of the Clausen function ``\\operatorname{Cl}_6(x)``
-for a real angle ``x`` of type `Float64`.  This function is defined as
+for a real angle ``x`` of type `Real`.  This function is defined as
 
 ```math
 \\operatorname{Cl}_6(x) = \\Im[\\operatorname{Li}_6(e^{ix})] = \\sum_{k=1}^\\infty \\frac{\\sin(kx)}{k^6}
@@ -18,7 +18,13 @@ julia> cl6(1.0)
 0.855629273183937
 ```
 """
-function cl6(x::Float64)::Float64
+cl6(x::Real) = _cl6(float(x))
+
+_cl6(x::Float16) = oftype(x, _cl6(Float32(x)))
+
+_cl6(x::Float32) = oftype(x, _cl6(Float64(x)))
+
+function _cl6(x::Float64)::Float64
     zeta3 = 1.2020569031595943
     pi28 = pi*pi/8.0
 

src/Sl.jl

@@ -30,11 +30,11 @@ function sl_series(n::Integer, x::Float64)
 end
 
 """
-    sl(n::Integer, x::Float64)::Float64
+    sl(n::Integer, x::Real)::Real
 
 Returns the value of the Glaisher-Clausen function
 ``\\operatorname{Sl}_n(x)`` for integers ``n > 0`` and a real angle
-``x`` of type `Float64`.  This function is defined as
+``x`` of type `Real`.  This function is defined as
 
 ```math
 \\begin{aligned}
@@ -56,7 +56,13 @@ julia> sl(10, 1.0)
 0.5398785706335891
 ```
 """
-function sl(n::Integer, x::Float64)::Float64
+sl(n::Integer, x::Real) = _sl(n, float(x))
+
+_sl(n::Integer, x::Float16) = oftype(x, _sl(n, Float32(x)))
+
+_sl(n::Integer, x::Float32) = oftype(x, _sl(n, Float64(x)))
+
+function _sl(n::Integer, x::Float64)::Float64
     n < 1 && throw(DomainError(n, "sl(n,x) undefined for n < 1"))
 
     (x, sgn) = range_reduce(n + 1, x)

src/range_reduction.jl

@@ -1,11 +1,11 @@
 # returns range-reduced x in [0,pi] for odd n
-function range_reduce_odd(x::Float64)
+function range_reduce_odd(x::Real)
     if x < zero(x)
         x = -x
     end
 
-    if x >= 2.0*pi
-        x = mod2pi(x)
+    if x >= 2pi
+        x = mod(x, 2pi)
     end
 
     if x > pi
@@ -18,16 +18,16 @@ function range_reduce_odd(x::Float64)
 end
 
 # returns (x, sign) with range-reduced x in [0,pi] for even n
-function range_reduce_even(x::Float64)
+function range_reduce_even(x::Real)
     sgn = one(x)
 
     if x < zero(x)
         x = -x
         sgn = -one(x)
     end
 
-    if x >= 2.0*pi
-        x = mod2pi(x)
+    if x >= 2pi
+        x = mod(x, 2pi)
     end
 
     if x > pi
@@ -41,7 +41,7 @@ function range_reduce_even(x::Float64)
 end
 
 # returns (x, sign) with range-reduced x in [0,pi]
-function range_reduce(n::Integer, x::Float64)
+function range_reduce(n::Integer, x::Real)
     if iseven(n)
         range_reduce_even(x)
     else

test/Cl.jl

@@ -13,10 +13,19 @@
 
             @test cl == clm
             @test cl ≈ expected atol=1e-14
+            @test ClausenFunctions.cl(n, Float16(x)) ≈ Float16(expected) atol=30*eps(Float16) rtol=30*eps(Float16)
+            @test ClausenFunctions.cl(n, Float32(x)) ≈ Float32(expected) atol=30*eps(Float32) rtol=30*eps(Float32)
         end
 
     end
 
+    @test ClausenFunctions.cl(1, 1//2) ≈ 0.70358563513784466 atol=1e-14
+    @test ClausenFunctions.cl(2, 1//2) ≈ 0.84831187770367927 atol=1e-14
+    @test ClausenFunctions.cl(3, 1//2) ≈ 0.92769631047023043 atol=1e-14
+    @test ClausenFunctions.cl(4, 1//2) ≈ 0.54837172654589549 atol=1e-14
+    @test ClausenFunctions.cl(5, 1//2) ≈ 0.89390286951083851 atol=1e-14
+    @test ClausenFunctions.cl(6, 1//2) ≈ 0.49419627977618802 atol=1e-14
+
     @test_throws DomainError ClausenFunctions.cl(0, 1.0)
     @test_throws DomainError ClausenFunctions.cl(-1, 1.0)
     @test_throws DomainError ClausenFunctions.cl(-2, 1.0)

test/Cl1.jl

@@ -12,5 +12,6 @@
         @test ClausenFunctions.cl1(x + 2.0*pi) ≈ expected atol=1e-12
     end
 
+    @test ClausenFunctions.cl1(1//2) ≈ 0.70358563513784466 atol=1e-14
     @test ClausenFunctions.cl1(0.0) == Inf
 end

test/Cl2.jl

@@ -10,5 +10,11 @@
         @test ClausenFunctions.cl2(-x) ≈ -expected atol=1e-14
         @test ClausenFunctions.cl2(x - 2.0*pi) ≈ expected atol=1e-13
         @test ClausenFunctions.cl2(x + 2.0*pi) ≈ expected atol=1e-13
+
+        @test ClausenFunctions.cl2(Float16(x)) ≈ Float16(expected) atol=30*eps(Float16) rtol=30*eps(Float16)
+        @test ClausenFunctions.cl2(Float32(x)) ≈ Float32(expected) atol=30*eps(Float32) rtol=30*eps(Float32)
     end
+
+    @test ClausenFunctions.cl2(pi/2) ≈ 0.915965594177219015054603514932384110 atol=1e-14
+    @test ClausenFunctions.cl2(1//2) ≈ 0.84831187770367927 atol=1e-14
 end

test/Cl3.jl

@@ -10,5 +10,10 @@
         @test ClausenFunctions.cl3(-x) ≈ expected atol=1e-14
         @test ClausenFunctions.cl3(x - 2.0*pi) ≈ expected atol=1e-13
         @test ClausenFunctions.cl3(x + 2.0*pi) ≈ expected atol=1e-13
+
+        @test ClausenFunctions.cl3(Float16(x)) ≈ Float16(expected) atol=30*eps(Float16) rtol=30*eps(Float16)
+        @test ClausenFunctions.cl3(Float32(x)) ≈ Float32(expected) atol=30*eps(Float32) rtol=30*eps(Float32)
     end
+
+    @test ClausenFunctions.cl3(1//2) ≈ 0.92769631047023043 atol=1e-14
 end

test/Cl4.jl

@@ -10,5 +10,10 @@
         @test ClausenFunctions.cl4(-x) ≈ -expected atol=1e-14
         @test ClausenFunctions.cl4(x - 2.0*pi) ≈ expected atol=1e-13
         @test ClausenFunctions.cl4(x + 2.0*pi) ≈ expected atol=1e-13
+
+        @test ClausenFunctions.cl4(Float16(x)) ≈ Float16(expected) atol=30*eps(Float16) rtol=30*eps(Float16)
+        @test ClausenFunctions.cl4(Float32(x)) ≈ Float32(expected) atol=30*eps(Float32) rtol=30*eps(Float32)
     end
+
+    @test ClausenFunctions.cl4(1//2) ≈ 0.54837172654589549 atol=1e-14
 end

test/Cl5.jl

@@ -10,5 +10,10 @@
         @test ClausenFunctions.cl5(-x) ≈ expected atol=1e-14
         @test ClausenFunctions.cl5(x - 2.0*pi) ≈ expected atol=1e-13
         @test ClausenFunctions.cl5(x + 2.0*pi) ≈ expected atol=1e-13
+
+        @test ClausenFunctions.cl5(Float16(x)) ≈ Float16(expected) atol=30*eps(Float16) rtol=30*eps(Float16)
+        @test ClausenFunctions.cl5(Float32(x)) ≈ Float32(expected) atol=30*eps(Float32) rtol=30*eps(Float32)
     end
+
+    @test ClausenFunctions.cl5(1//2) ≈ 0.89390286951083851 atol=1e-14
 end

test/Cl6.jl

@@ -10,5 +10,10 @@
         @test ClausenFunctions.cl6(-x) ≈ -expected atol=1e-14
         @test ClausenFunctions.cl6(x - 2.0*pi) ≈ expected atol=1e-13
         @test ClausenFunctions.cl6(x + 2.0*pi) ≈ expected atol=1e-13
+
+        @test ClausenFunctions.cl6(Float16(x)) ≈ Float16(expected) atol=30*eps(Float16) rtol=30*eps(Float16)
+        @test ClausenFunctions.cl6(Float32(x)) ≈ Float32(expected) atol=30*eps(Float32) rtol=30*eps(Float32)
     end
+
+    @test ClausenFunctions.cl6(1//2) ≈ 0.49419627977618802 atol=1e-14
 end

test/Sl.jl

@@ -26,9 +26,18 @@
 
             @test ClausenFunctions.sl(n, x) == (-1)^n*ClausenFunctions.sl(n, -x)
             @test ClausenFunctions.sl(n, x) ≈ expected atol=1e-14
+            @test ClausenFunctions.sl(n, Float16(x)) ≈ Float16(expected) atol=30*eps(Float16) rtol=30*eps(Float16)
+            @test ClausenFunctions.sl(n, Float32(x)) ≈ Float32(expected) atol=30*eps(Float32) rtol=30*eps(Float32)
         end
     end
 
+    @test ClausenFunctions.sl(1, 1//2) ≈ 1.3207963267948966 atol=1e-14
+    @test ClausenFunctions.sl(2, 1//2) ≈ 0.92203590345077813 atol=1e-14
+    @test ClausenFunctions.sl(3, 1//2) ≈ 0.63653415924141781 atol=1e-14
+    @test ClausenFunctions.sl(4, 1//2) ≈ 0.90812931549667023 atol=1e-14
+    @test ClausenFunctions.sl(5, 1//2) ≈ 0.51085256423059275 atol=1e-14
+    @test ClausenFunctions.sl(6, 1//2) ≈ 0.88593812938731573 atol=1e-14
+
     @test_throws DomainError ClausenFunctions.sl(0, 1.0)
     @test_throws DomainError ClausenFunctions.sl(-1, 1.0)
     @test_throws DomainError ClausenFunctions.sl(-2, 1.0)

test/data/Sl1.txt

@@ -1,3 +1,4 @@
+0.0000000000000000000000000000000000000000      0.0
 0.06283185307179586476925286766559005768394	1.539380400258998686846695257806956413257
 0.1256637061435917295385057353311801153679	1.507964473723100754462068823974161384415
 0.1884955592153875943077586029967701730518	1.476548547187202822077442390141366355573
@@ -97,4 +98,3 @@
 6.094689747964198882617528163562235595343	-1.476548547187202822077442390141366355573
 6.157521601035994747386781031227825653026	-1.507964473723100754462068823974161384415
 6.22035345410779061215603389889341571071	-1.539380400258998686846695257806956413257
-6.283185307179586476925286766559005768394	-1.570796326794896619231321691639751442099