diff --git a/src/SpecialFunctions.jl b/src/SpecialFunctions.jl
index c4fe6120..2869a1cd 100644
--- a/src/SpecialFunctions.jl
+++ b/src/SpecialFunctions.jl
@@ -78,6 +78,7 @@ export
     expintx,
     sinint,
     cosint,
+    cosint2,
     lbinomial
 
 include("bessel.jl")
diff --git a/src/chainrules.jl b/src/chainrules.jl
index 76a8f07d..52b0b44f 100644
--- a/src/chainrules.jl
+++ b/src/chainrules.jl
@@ -193,6 +193,7 @@ ChainRulesCore.@scalar_rule(
 ChainRulesCore.@scalar_rule(expinti(x), exp(x) / x)
 ChainRulesCore.@scalar_rule(sinint(x), sinc(invπ * x))
 ChainRulesCore.@scalar_rule(cosint(x), cos(x) / x)
+ChainRulesCore.@scalar_rule(cosint2(x),  iszero(x) ? zero(x) : (one(x)-cos(x))/x )
 
 # elliptic integrals
 ChainRulesCore.@scalar_rule(
diff --git a/src/sincosint.jl b/src/sincosint.jl
index 783ec6b0..9d8fe2c1 100644
--- a/src/sincosint.jl
+++ b/src/sincosint.jl
@@ -216,7 +216,18 @@ function cosint(x::Float64)
     end
 end
 
-for f in (:sinint, :cosint)
+function cosint2(x::Float64)
+    y = abs(x)
+    if y < 0.5
+        z = y*y
+        return z * @horner(z, 1/4, -1/96, 1/4320, -1/322560, 1/36288000, 1/5748019200)
+    else
+        return MathConstants.eulergamma + log(y) - cosint(y)
+    end
+end
+
+
+for f in (:sinint, :cosint, :cosint2)
     @eval begin
         ($f)(x::Float32)        = Float32(($f)(Float64(x)))
         ($f)(x::Float16)        = Float16(($f)(Float64(x)))
diff --git a/test/sincosint.jl b/test/sincosint.jl
index 6251111c..4e8d872d 100644
--- a/test/sincosint.jl
+++ b/test/sincosint.jl
@@ -23,16 +23,22 @@
 
     @test sinint(Float16(1)) == Float16(sinint(1))
     @test cosint(Float16(1)) == Float16(cosint(1))
+    @test cosint2(Float16(1)) == Float16(cosint2(1))
     @test sinint(Float32(1)) == Float32(sinint(1))
     @test cosint(Float32(1)) == Float32(cosint(1))
+    @test cosint2(Float32(1)) == Float32(cosint2(1))
 
     @test sinint(1//2) == sinint(0.5)
     @test cosint(1//2) == cosint(0.5)
+    @test cosint2(1//2) == cosint2(0.5)
 
     @test sinint(1e300)     ≅ π/2
     @test cosint(1e300)     ≅ -8.17881912115908554103E-301
     @test sinint(1e-300)    ≅ 1.0E-300
     @test cosint(1e-300)    ≅ -690.1983122333121
+    
+    @test cosint2(0.0) == 0.0
+    @test cosint2(0.05) ≅ 0.0006249348994501105 # Computed with Mathematica: SetPrecision[EulerGamma + Log[1/20] - CosIntegral[1/20], 20]
 
     @test sinint(Inf) == π/2
     @test cosint(Inf) == 0.0
@@ -41,10 +47,12 @@
 
     @test_throws ErrorException sinint(big(1.0))
     @test_throws ErrorException cosint(big(1.0))
+    @test_throws ErrorException cosint2(big(1.0))
     @test_throws DomainError cosint(-1.0)
     @test_throws DomainError cosint(Float32(-1.0))
     @test_throws DomainError cosint(Float16(-1.0))
     @test_throws DomainError cosint(-1//2)
     @test_throws MethodError sinint(Complex(1))
     @test_throws MethodError cosint(Complex(1))
+    @test_throws MethodError cosint2(Complex(1))
 end