diff --git a/docs/src/mpolynomial.md b/docs/src/mpolynomial.md
index d1b44e3286..447c7c8e5e 100644
--- a/docs/src/mpolynomial.md
+++ b/docs/src/mpolynomial.md
@@ -321,6 +321,10 @@ The following functionality is also provided for all multivariate polynomials.
 is_univariate(::MPolyRing{T}) where T <: RingElement
 ```
 
+```@docs
+var_indices(p::MPolyRingElem{T}) where T <: RingElement
+```
+
 ```@docs
 vars(p::MPolyRingElem{T}) where T <: RingElement
 ```
diff --git a/src/MPoly.jl b/src/MPoly.jl
index 96a668f193..5579082933 100644
--- a/src/MPoly.jl
+++ b/src/MPoly.jl
@@ -89,6 +89,19 @@ Return the number of variables in `R`.
 """
 number_of_generators(R::MPolyRing) = number_of_variables(R)
 
+@doc raw"""
+    var_indices(p::MPolyRingElem{T}) where {T <: RingElement}
+
+Return the indices of the variables actually occurring in $p$.
+"""
+function var_indices(p::MPolyRingElem{T}) where {T <: RingElement}
+   isused = zeros(Int, nvars(parent(p)))
+   for v in exponent_vectors(p)
+      isused .|= v  # accumulate by OR-ing the exponent vectors
+   end
+   return findall(!iszero, isused)
+end
+
 @doc raw"""
     vars(p::MPolyRingElem{T}) where {T <: RingElement}
 
@@ -96,12 +109,7 @@ Return the variables actually occurring in $p$.
 """
 function vars(p::MPolyRingElem{T}) where {T <: RingElement}
    R = parent(p)
-   n = nvars(R)
-   isused = zeros(Int, n)
-   for v in exponent_vectors(p)
-      isused .|= v  # accumulate by OR-ing the exponent vectors
-   end
-   return [R[i] for i in 1:n if isused[i] != 0]
+   return [gen(R, i) for i in var_indices(p)]
 end
 
 @doc raw"""
diff --git a/src/exports.jl b/src/exports.jl
index 049c9c7dbf..5dd005a49f 100644
--- a/src/exports.jl
+++ b/src/exports.jl
@@ -598,6 +598,7 @@ export use_karamul
 export valuation
 export var
 export var_index
+export var_indices
 export vars
 export vector_space
 export vector_space_dim
diff --git a/src/generic/MPoly.jl b/src/generic/MPoly.jl
index f8e37b72ae..713f90426a 100644
--- a/src/generic/MPoly.jl
+++ b/src/generic/MPoly.jl
@@ -77,18 +77,17 @@ function gen(a::MPolyRing{T}, i::Int) where {T <: RingElement}
    return a([one(base_ring(a))], exps)
 end
 
-function vars(p::MPoly{T}) where {T <: RingElement}
-   vars_in_p = Vector{MPoly{T}}(undef, 0)
-   n = nvars(p.parent)
+function var_indices(p::MPoly)
+   vars_in_p = Int[]
+   n = nvars(parent(p))
    exps = p.exps
-   gen_list = gens(p.parent)
    for j = 1:n
       for i = 1:length(p)
          if exps[j, i] > 0
             if p.parent.ord == :degrevlex
-               push!(vars_in_p, gen_list[j])
+               push!(vars_in_p, j)
             else
-               push!(vars_in_p, gen_list[n - j + 1])
+               push!(vars_in_p, n - j + 1)
             end
             break
          end
@@ -97,7 +96,12 @@ function vars(p::MPoly{T}) where {T <: RingElement}
    if p.parent.ord != :degrevlex
       vars_in_p = reverse(vars_in_p)
    end
-   return(vars_in_p)
+   return vars_in_p
+end
+
+function vars(p::MPoly{T}) where {T <: RingElement}
+   R = parent(p)
+   return MPoly{T}[gen(R, i) for i in var_indices(p)]
 end
 
 @doc raw"""
diff --git a/src/generic/imports.jl b/src/generic/imports.jl
index 928fd83b8f..28957330bb 100644
--- a/src/generic/imports.jl
+++ b/src/generic/imports.jl
@@ -230,6 +230,7 @@ import ..AbstractAlgebra: use_karamul
 import ..AbstractAlgebra: valuation
 import ..AbstractAlgebra: var
 import ..AbstractAlgebra: var_index
+import ..AbstractAlgebra: var_indices
 import ..AbstractAlgebra: vars
 import ..AbstractAlgebra: vector_space_dim
 import ..AbstractAlgebra: zero!