diff --git a/NEWS.md b/NEWS.md
index 65dbd86392592..e0a5a8073d6c7 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -300,6 +300,8 @@ Deprecated or removed
     full path if you need access to executables or libraries in the `JULIA_HOME` directory, e.g.
     `joinpath(JULIA_HOME, "7z.exe")` for `7z.exe` ([#21540]).
 
+  * `expm` has been deprecated in favor of `exp` ([#23233]).
+
   * Calling `union` with no arguments is deprecated; construct an empty set with an appropriate
     element type using `Set{T}()` instead ([#23144]).
 
@@ -1165,3 +1167,4 @@ Command-line option changes
 [#23157]: https://github.com/JuliaLang/julia/issues/23157
 [#23187]: https://github.com/JuliaLang/julia/issues/23187
 [#23207]: https://github.com/JuliaLang/julia/issues/23207
+[#23233]: https://github.com/JuliaLang/julia/issues/23233
diff --git a/base/deprecated.jl b/base/deprecated.jl
index a9a5c93e269e5..4f89a34e45113 100644
--- a/base/deprecated.jl
+++ b/base/deprecated.jl
@@ -220,7 +220,7 @@ for f in (:sin, :sinh, :sind, :asin, :asinh, :asind,
         :tan, :tanh, :tand, :atan, :atanh, :atand,
         :sinpi, :cosc, :ceil, :floor, :trunc, :round,
         :log1p, :expm1, :abs, :abs2,
-        :log, :log2, :log10, :exp, :exp2, :exp10, :sinc, :cospi,
+        :log, :log2, :log10, :exp2, :exp10, :sinc, :cospi,
         :cos, :cosh, :cosd, :acos, :acosd,
         :cot, :coth, :cotd, :acot, :acotd,
         :sec, :sech, :secd, :asech,
@@ -251,7 +251,7 @@ for f in (
         # base/special/gamma.jl
         :gamma, :lfact,
         # base/math.jl
-        :cbrt, :sinh, :cosh, :tanh, :atan, :asinh, :exp, :exp2,
+        :cbrt, :sinh, :cosh, :tanh, :atan, :asinh, :exp2,
         :expm1, :exp10, :sin, :cos, :tan, :asin, :acos, :acosh, :atanh,
         #=:log,=# :log2, :log10, :lgamma, #=:log1p,=# :sqrt,
         # base/floatfuncs.jl
@@ -1653,6 +1653,10 @@ function SymTridiagonal(dv::AbstractVector{T}, ev::AbstractVector{S}) where {T,S
     SymTridiagonal(convert(Vector{R}, dv), convert(Vector{R}, ev))
 end
 
+# deprecate expm in favor of exp
+@deprecate expm! exp!
+@deprecate expm exp
+
 # PR #23092
 @eval LibGit2 begin
     function prompt(msg::AbstractString; default::AbstractString="", password::Bool=false)
diff --git a/base/exports.jl b/base/exports.jl
index eecd75f36d46e..adcbaee1615e5 100644
--- a/base/exports.jl
+++ b/base/exports.jl
@@ -579,7 +579,6 @@ export
     eigvals,
     eigvals!,
     eigvecs,
-    expm,
     eye,
     factorize,
     givens,
diff --git a/base/linalg/dense.jl b/base/linalg/dense.jl
index bde22d2c2d58a..68193231f6880 100644
--- a/base/linalg/dense.jl
+++ b/base/linalg/dense.jl
@@ -420,12 +420,12 @@ function (^)(A::AbstractMatrix{T}, p::Real) where T
     # Otherwise, use Schur decomposition
     return schurpow(A, p)
 end
-(^)(A::AbstractMatrix, p::Number) = expm(p*logm(A))
+(^)(A::AbstractMatrix, p::Number) = exp(p*logm(A))
 
 # Matrix exponential
 
 """
-    expm(A)
+    exp(A::AbstractMatrix)
 
 Compute the matrix exponential of `A`, defined by
 
@@ -445,22 +445,21 @@ julia> A = eye(2, 2)
  1.0  0.0
  0.0  1.0
 
-julia> expm(A)
+julia> exp(A)
 2×2 Array{Float64,2}:
  2.71828  0.0
  0.0      2.71828
 ```
 """
-expm(A::StridedMatrix{<:BlasFloat}) = expm!(copy(A))
-expm(A::StridedMatrix{<:Integer}) = expm!(float(A))
-expm(x::Number) = exp(x)
+exp(A::StridedMatrix{<:BlasFloat}) = exp!(copy(A))
+exp(A::StridedMatrix{<:Integer}) = exp!(float(A))
 
 ## Destructive matrix exponential using algorithm from Higham, 2008,
 ## "Functions of Matrices: Theory and Computation", SIAM
-function expm!(A::StridedMatrix{T}) where T<:BlasFloat
+function exp!(A::StridedMatrix{T}) where T<:BlasFloat
     n = checksquare(A)
     if ishermitian(A)
-        return full(expm(Hermitian(A)))
+        return full(exp(Hermitian(A)))
     end
     ilo, ihi, scale = LAPACK.gebal!('B', A)    # modifies A
     nA   = norm(A, 1)
diff --git a/base/linalg/diagonal.jl b/base/linalg/diagonal.jl
index 439e409958c52..c98087b915802 100644
--- a/base/linalg/diagonal.jl
+++ b/base/linalg/diagonal.jl
@@ -326,8 +326,7 @@ end
 eye(::Type{Diagonal{T}}, n::Int) where {T} = Diagonal(ones(T,n))
 
 # Matrix functions
-expm(D::Diagonal) = Diagonal(exp.(D.diag))
-expm(D::Diagonal{<:AbstractMatrix}) = Diagonal(expm.(D.diag))
+exp(D::Diagonal) = Diagonal(exp.(D.diag))
 logm(D::Diagonal) = Diagonal(log.(D.diag))
 logm(D::Diagonal{<:AbstractMatrix}) = Diagonal(logm.(D.diag))
 sqrtm(D::Diagonal) = Diagonal(sqrt.(D.diag))
diff --git a/base/linalg/linalg.jl b/base/linalg/linalg.jl
index a510f0efc8b51..50f1e9a1aaebe 100644
--- a/base/linalg/linalg.jl
+++ b/base/linalg/linalg.jl
@@ -6,7 +6,7 @@ import Base: \, /, *, ^, +, -, ==
 import Base: A_mul_Bt, At_ldiv_Bt, A_rdiv_Bc, At_ldiv_B, Ac_mul_Bc, A_mul_Bc, Ac_mul_B,
     Ac_ldiv_B, Ac_ldiv_Bc, At_mul_Bt, A_rdiv_Bt, At_mul_B
 import Base: USE_BLAS64, abs, big, broadcast, ceil, conj, convert, copy, copy!,
-    ctranspose, eltype, eye, findmax, findmin, fill!, floor, full, getindex,
+    ctranspose, eltype, exp, eye, findmax, findmin, fill!, floor, full, getindex,
     hcat, imag, indices, inv, isapprox, isone, IndexStyle, kron, length, map,
     ndims, oneunit, parent, power_by_squaring, print_matrix, promote_rule, real, round,
     setindex!, show, similar, size, transpose, trunc, typed_hcat
@@ -80,7 +80,6 @@ export
     eigvals,
     eigvals!,
     eigvecs,
-    expm,
     eye,
     factorize,
     givens,
diff --git a/base/linalg/symmetric.jl b/base/linalg/symmetric.jl
index 8d98709daaf7f..95c8b824e684e 100644
--- a/base/linalg/symmetric.jl
+++ b/base/linalg/symmetric.jl
@@ -589,11 +589,11 @@ function ^(A::Hermitian{T}, p::Real) where T
     end
 end
 
-function expm(A::Symmetric)
+function exp(A::Symmetric)
     F = eigfact(A)
     return Symmetric((F.vectors * Diagonal(exp.(F.values))) * F.vectors')
 end
-function expm(A::Hermitian{T}) where T
+function exp(A::Hermitian{T}) where T
     n = checksquare(A)
     F = eigfact(A)
     retmat = (F.vectors * Diagonal(exp.(F.values))) * F.vectors'
diff --git a/doc/src/manual/linear-algebra.md b/doc/src/manual/linear-algebra.md
index 8bb0e2101a8d6..3a0fb0fb4e6cc 100644
--- a/doc/src/manual/linear-algebra.md
+++ b/doc/src/manual/linear-algebra.md
@@ -177,8 +177,8 @@ as well as whether hooks to various optimized methods for them in LAPACK are ava
 
 | Matrix type               | `+` | `-` | `*` | `\` | Other functions with optimized methods                              |
 |:------------------------- |:--- |:--- |:--- |:--- |:------------------------------------------------------------------- |
-| [`Symmetric`](@ref)       |     |     |     | MV  | [`inv()`](@ref), [`sqrtm()`](@ref), [`expm()`](@ref)                |
-| [`Hermitian`](@ref)       |     |     |     | MV  | [`inv()`](@ref), [`sqrtm()`](@ref), [`expm()`](@ref)                |
+| [`Symmetric`](@ref)       |     |     |     | MV  | [`inv()`](@ref), [`sqrtm()`](@ref), [`exp()`](@ref)                |
+| [`Hermitian`](@ref)       |     |     |     | MV  | [`inv()`](@ref), [`sqrtm()`](@ref), [`exp()`](@ref)                |
 | [`UpperTriangular`](@ref) |     |     | MV  | MV  | [`inv()`](@ref), [`det()`](@ref)                                    |
 | [`LowerTriangular`](@ref) |     |     | MV  | MV  | [`inv()`](@ref), [`det()`](@ref)                                    |
 | [`SymTridiagonal`](@ref)  | M   | M   | MS  | MV  | [`eigmax()`](@ref), [`eigmin()`](@ref)                              |
diff --git a/test/linalg/dense.jl b/test/linalg/dense.jl
index 7eca457de6f2b..98a8e7ac2de43 100644
--- a/test/linalg/dense.jl
+++ b/test/linalg/dense.jl
@@ -139,7 +139,6 @@ bimg  = randn(n,2)/2
         A = zeros(eltya,1,1)
         A[1,1] = α
         @test diagm(α) == A # Test behavior of `diagm` when passed a scalar
-        @test expm(α) == exp(α) # `expm` should behave like `exp` with scalar argument
     end
 
     @testset "Factorize" begin
@@ -416,7 +415,7 @@ end
         eA1 = convert(Matrix{elty}, [147.866622446369 127.781085523181  127.781085523182;
                                      183.765138646367 183.765138646366  163.679601723179;
                                      71.797032399996  91.8825693231832 111.968106246371]')
-        @test expm(A1) ≈ eA1
+        @test exp(A1) ≈ eA1
 
         A2  = convert(Matrix{elty},
                       [29.87942128909879    0.7815750847907159 -2.289519314033932;
@@ -426,21 +425,21 @@ end
                       [  5496313853692458.0 -18231880972009236.0 -30475770808580460.0;
                        -18231880972009252.0  60605228702221920.0 101291842930249760.0;
                        -30475770808580480.0 101291842930249728.0 169294411240851968.0])
-        @test expm(A2) ≈ eA2
+        @test exp(A2) ≈ eA2
 
         A3  = convert(Matrix{elty}, [-131 19 18;-390 56 54;-387 57 52])
         eA3 = convert(Matrix{elty}, [-1.50964415879218 -5.6325707998812  -4.934938326092;
                                      0.367879439109187 1.47151775849686  1.10363831732856;
                                      0.135335281175235 0.406005843524598 0.541341126763207]')
-        @test expm(A3) ≈ eA3
+        @test exp(A3) ≈ eA3
 
         A4 = convert(Matrix{elty}, [0.25 0.25; 0 0])
         eA4 = convert(Matrix{elty}, [1.2840254166877416 0.2840254166877415; 0 1])
-        @test expm(A4) ≈ eA4
+        @test exp(A4) ≈ eA4
 
         A5 = convert(Matrix{elty}, [0 0.02; 0 0])
         eA5 = convert(Matrix{elty}, [1 0.02; 0 1])
-        @test expm(A5) ≈ eA5
+        @test exp(A5) ≈ eA5
 
         # Hessenberg
         @test hessfact(A1)[:H] ≈ convert(Matrix{elty},
@@ -454,20 +453,20 @@ end
                                      1/3 1/4 1/5 1/6;
                                      1/4 1/5 1/6 1/7;
                                      1/5 1/6 1/7 1/8])
-        @test expm(logm(A4)) ≈ A4
+        @test exp(logm(A4)) ≈ A4
 
         A5  = convert(Matrix{elty}, [1 1 0 1; 0 1 1 0; 0 0 1 1; 1 0 0 1])
-        @test expm(logm(A5)) ≈ A5
+        @test exp(logm(A5)) ≈ A5
 
         A6  = convert(Matrix{elty}, [-5 2 0 0 ; 1/2 -7 3 0; 0 1/3 -9 4; 0 0 1/4 -11])
-        @test expm(logm(A6)) ≈ A6
+        @test exp(logm(A6)) ≈ A6
 
         A7  = convert(Matrix{elty}, [1 0 0 1e-8; 0 1 0 0; 0 0 1 0; 0 0 0 1])
-        @test expm(logm(A7)) ≈ A7
+        @test exp(logm(A7)) ≈ A7
     end
 
     A8 = 100 * [-1+1im 0 0 1e-8; 0 1 0 0; 0 0 1 0; 0 0 0 1]
-    @test expm(logm(A8)) ≈ A8
+    @test exp(logm(A8)) ≈ A8
 end
 
 @testset "issue 5116" begin
@@ -476,19 +475,19 @@ end
            0.006540706968939  -0.999786072879326   0.0   0.0
            0.0                 0.0                 1.0   0.0
            0.013081413937878  -3.999572145758650   0.0   1.0]
-    @test expm(A9) ≈ eA9
+    @test exp(A9) ≈ eA9
 
     A10  = [ 0. 0. 0. 0. ; 0. 0. -im 0.; 0. im 0. 0.; 0. 0. 0. 0.]
     eA10 = [ 1.0+0.0im   0.0+0.0im                 0.0+0.0im                0.0+0.0im
             0.0+0.0im   1.543080634815244+0.0im   0.0-1.175201193643801im  0.0+0.0im
             0.0+0.0im   0.0+1.175201193643801im   1.543080634815243+0.0im  0.0+0.0im
             0.0+0.0im   0.0+0.0im                 0.0+0.0im                1.0+0.0im]
-    @test expm(A10) ≈ eA10
+    @test exp(A10) ≈ eA10
 end
 
 @testset "Additional matrix logarithm tests" for elty in (Float64, Complex{Float64})
     A11 = convert(Matrix{elty}, [3 2; -5 -3])
-    @test expm(logm(A11)) ≈ A11
+    @test exp(logm(A11)) ≈ A11
 
     A12 = convert(Matrix{elty}, [1 -1; 1 -1])
     @test typeof(logm(A12)) == Array{Complex{Float64}, 2}
@@ -498,7 +497,7 @@ end
                                     0.2310490602 1.295566591 0.2651438179;
                                     0.2310490602 0.1969543025 1.363756107])
     @test logm(A1) ≈ logmA1
-    @test expm(logm(A1)) ≈ A1
+    @test exp(logm(A1)) ≈ A1
 
     A4  = convert(Matrix{elty}, [1/2 1/3 1/4 1/5+eps();
                                  1/3 1/4 1/5 1/6;
@@ -509,7 +508,7 @@ end
                                     0.4462766564 2.994142974 -7.351095988 3.318413247;
                                     0.2414170219 0.5865285289 3.318413247 -5.444632124])
     @test logm(A4) ≈ logmA4
-    @test expm(logm(A4)) ≈ A4
+    @test exp(logm(A4)) ≈ A4
 end
 
 @testset "issue #7181" begin
@@ -615,7 +614,7 @@ end
 
 @testset "test ops on Numbers for $elty" for elty in [Float32,Float64,Complex64,Complex128]
     a = rand(elty)
-    @test expm(a) == exp(a)
+    @test exp(a) == exp(a)
     @test isposdef(one(elty))
     @test sqrtm(a) == sqrt(a)
     @test logm(a) ≈ log(a)
diff --git a/test/linalg/diagonal.jl b/test/linalg/diagonal.jl
index b1e8df0bd3ffb..faf800f3653a7 100644
--- a/test/linalg/diagonal.jl
+++ b/test/linalg/diagonal.jl
@@ -60,7 +60,7 @@ srand(1)
             @test func(D) ≈ func(DM) atol=n^2*eps(relty)*(1+(elty<:Complex))
         end
         if relty <: BlasFloat
-            for func in (expm,)
+            for func in (exp,)
                 @test func(D) ≈ func(DM) atol=n^3*eps(relty)
             end
             @test logm(Diagonal(abs.(D.diag))) ≈ logm(abs.(DM)) atol=n^3*eps(relty)
@@ -381,7 +381,7 @@ end
     @test ishermitian(Dherm) == true
     @test ishermitian(Dsym) == false
 
-    @test expm(D) == Diagonal([expm([1 2; 3 4]), expm([1 2; 3 4])])
+    @test exp(D) == Diagonal([exp([1 2; 3 4]), exp([1 2; 3 4])])
     @test logm(D) == Diagonal([logm([1 2; 3 4]), logm([1 2; 3 4])])
     @test sqrtm(D) == Diagonal([sqrtm([1 2; 3 4]), sqrtm([1 2; 3 4])])
 end
diff --git a/test/linalg/lapack.jl b/test/linalg/lapack.jl
index 7d5d92b5204db..56838b1410f2a 100644
--- a/test/linalg/lapack.jl
+++ b/test/linalg/lapack.jl
@@ -613,7 +613,7 @@ end
 
 # Issue 13976
 let A = [NaN 0.0 NaN; 0 0 0; NaN 0 NaN]
-    @test_throws ArgumentError expm(A)
+    @test_throws ArgumentError exp(A)
 end
 
 # Issue 14065 (and 14220)
diff --git a/test/linalg/symmetric.jl b/test/linalg/symmetric.jl
index 54cfb2753af27..b0a2734f03518 100644
--- a/test/linalg/symmetric.jl
+++ b/test/linalg/symmetric.jl
@@ -12,16 +12,16 @@ end
 @testset "Hermitian matrix exponential/log" begin
     A1 = randn(4,4) + im*randn(4,4)
     A2 = A1 + A1'
-    @test expm(A2) ≈ expm(Hermitian(A2))
+    @test exp(A2) ≈ exp(Hermitian(A2))
     @test logm(A2) ≈ logm(Hermitian(A2))
     A3 = A1 * A1' # posdef
-    @test expm(A3) ≈ expm(Hermitian(A3))
+    @test exp(A3) ≈ exp(Hermitian(A3))
     @test logm(A3) ≈ logm(Hermitian(A3))
 
     A1 = randn(4,4)
     A3 = A1 * A1'
     A4 = A1 + A1.'
-    @test expm(A4) ≈ expm(Symmetric(A4))
+    @test exp(A4) ≈ exp(Symmetric(A4))
     @test logm(A3) ≈ logm(Symmetric(A3))
     @test logm(A3) ≈ logm(Hermitian(A3))
 end
diff --git a/test/linalg/triangular.jl b/test/linalg/triangular.jl
index f6dbe41552cc0..05643f243df4d 100644
--- a/test/linalg/triangular.jl
+++ b/test/linalg/triangular.jl
@@ -174,9 +174,9 @@ for elty1 in (Float32, Float64, BigFloat, Complex64, Complex128, Complex{BigFloa
             @test B == viewA1.'
         end
 
-        #expm/logm
+        #exp/logm
         if (elty1 == Float64 || elty1 == Complex128) && (t1 == UpperTriangular || t1 == LowerTriangular)
-            @test expm(full(logm(A1))) ≈ full(A1)
+            @test exp(full(logm(A1))) ≈ full(A1)
         end
 
         # scale