diff --git a/Project.toml b/Project.toml
index 643ea0f47..e4bf95ef4 100644
--- a/Project.toml
+++ b/Project.toml
@@ -9,6 +9,8 @@ ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 DensityInterface = "b429d917-457f-4dbc-8f4c-0cc954292b1d"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
+IrrationalConstants = "92d709cd-6900-40b7-9082-c6be49f344b6"
 PDMats = "90014a1f-27ba-587c-ab20-58faa44d9150"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
@@ -43,6 +45,8 @@ FiniteDifferences = "0.12"
 ForwardDiff = "0.10, 1"
 JSON = "0.21"
 LinearAlgebra = "<0.0.1, 1"
+LogExpFunctions = "0.3.29"
+IrrationalConstants = "0.2.6"
 OffsetArrays = "1"
 PDMats = "0.11.35"
 Printf = "<0.0.1, 1"
diff --git a/src/Distributions.jl b/src/Distributions.jl
index fcf699434..b0d540fdb 100644
--- a/src/Distributions.jl
+++ b/src/Distributions.jl
@@ -25,8 +25,9 @@ import StatsBase: kurtosis, skewness, entropy, mode, modes,
 import PDMats: dim, PDMat, invquad
 
 using SpecialFunctions
+import LogExpFunctions: logabssinh, log1psq
 using Base.MathConstants: eulergamma
-
+import IrrationalConstants: invπ
 import AliasTables
 
 export
diff --git a/src/univariate/continuous/cosine.jl b/src/univariate/continuous/cosine.jl
index 4d11a9ddd..33a6527e7 100644
--- a/src/univariate/continuous/cosine.jl
+++ b/src/univariate/continuous/cosine.jl
@@ -59,37 +59,66 @@ kurtosis(d::Cosine{T}) where {T<:Real} = 6*(90-T(π)^4) / (5*(T(π)^2-6)^2)
 
 #### Evaluation
 
-function pdf(d::Cosine{T}, x::Real) where T<:Real
+function pdf(d::Cosine{T}, x::S) where {T<:Real, S<:Real}
     if insupport(d, x)
         z = (x - d.μ) / d.σ
         return (1 + cospi(z)) / (2d.σ)
     else
-        return zero(T)
+        return zero(promote_type(T, S))
     end
 end
 
-logpdf(d::Cosine, x::Real) = log(pdf(d, x))
-
-function cdf(d::Cosine{T}, x::Real) where T<:Real
-    if x < d.μ - d.σ
-        return zero(T)
-    end
-    if x > d.μ + d.σ
-        return one(T)
+function logpdf(d::Cosine{T}, x::S) where {T<:Real, S<:Real}
+    if insupport(d, x)
+        z = (x - d.μ) / d.σ
+        return log1p(cospi(z)) - log(2d.σ)
+    else
+        return typemin(promote_type(T, S))
     end
+end
+
+function cdf(d::Cosine{T}, x::S) where {T<:Real, S<:Real}
+    W = promote_type(T, S)
+
+    x < d.μ - d.σ && return zero(W)
+    x > d.μ + d.σ && return one(W)
+
     z = (x - d.μ) / d.σ
     (1 + z + sinpi(z) * invπ) / 2
 end
 
-function ccdf(d::Cosine{T}, x::Real) where T<:Real
-    if x < d.μ - d.σ
-        return one(T)
-    end
-    if x > d.μ + d.σ
-        return zero(T)
-    end
+function ccdf(d::Cosine{T}, x::S) where {T<:Real, S<:Real}
+    W = promote_type(T, S)
+
+    x < d.μ - d.σ && return one(W)
+    x > d.μ + d.σ && return zero(W)
+
     nz = (d.μ - x) / d.σ
     (1 + nz + sinpi(nz) * invπ) / 2
 end
 
 quantile(d::Cosine, p::Real) = quantile_bisect(d, p)
+
+function mgf(d::Cosine, t::Real)
+    σt, μt = d.σ * t, d.μ * t
+    z = iszero(σt) ? one(float(σt)) : sinh(σt)/σt
+    return exp(μt) * (z / (1 + (invπ * σt)^2))
+end
+
+function cgf(d::Cosine, t::Real)
+    σt, μt = d.σ * t, d.μ * t
+    z = iszero(σt) ? zero(float(σt)) : logabssinh(σt) - log(σt)
+    return μt + z - log1psq(invπ * σt)
+end
+
+function cf(d::Cosine, t::Real)
+    σt, μt = d.σ * t, d.μ * t
+    abs(σt) ≈ π && return cis(μt) / 2
+    z = iszero(σt) ? one(float(σt)) : sin(σt)/σt
+    return cis(μt) * z / (1 - (invπ * σt)^2)
+end
+
+#### Affine transformations
+
+Base.:+(d::Cosine, a::Real) = Cosine(d.μ + a, d.σ)
+Base.:*(c::Real, d::Cosine) = Cosine(c * d.μ, abs(c) * d.σ)
diff --git a/test/univariate/continuous/cosine.jl b/test/univariate/continuous/cosine.jl
new file mode 100644
index 000000000..e05b86da6
--- /dev/null
+++ b/test/univariate/continuous/cosine.jl
@@ -0,0 +1,13 @@
+@testset "Cosine" begin
+    @testset "mgf, cgf, cf" begin
+        @test cf(Cosine(1, 1), 1) ≈ cis(1) * sin(1)/(1 - π^-2)
+        @test mgf(Cosine(-1, 2), -1) ≈ exp(1) * sinh(-2) / (-2*(1 + 4π^-2))
+        @test cgf(Cosine(-1, 2), 2) ≈ -2 + log(sinh(4)) - log(4) - log(1 + 16π^-2)
+    end
+
+    @testset "affine transformations" begin
+        @test -Cosine(1, 1) == Cosine(-1, 1)
+        @test 3Cosine(1, 2) == Cosine(3, 6)
+        @test Cosine(0, 2) + 42 == Cosine(42, 2)
+    end
+end
\ No newline at end of file