diff --git a/Project.toml b/Project.toml
index c269e58..1b9c39d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -8,8 +8,14 @@ Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 Oscar = "f1435218-dba5-11e9-1e4d-f1a5fab5fc13"
 
+[extras]
+AbstractAlgebra = "c3fe647b-3220-5bb0-a1ea-a7954cac585d"
+Hecke = "3e1990a7-5d81-5526-99ce-9ba3ff248f21"
+
 [compat]
+AbstractAlgebra = "0.44.12"
 Compat = "4.11"
+Hecke = "0.35.13"
 JSON = "^0.20, ^0.21"
 Oscar = "1.3"
 julia = "1.6"
diff --git a/docs/make.jl b/docs/make.jl
index fa02f5f..9bed9bc 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -33,6 +33,7 @@ makedocs(;
     "modify.md",
     "cyclo.md",
     "book.md",
+    "printing.md",
     "unexported.md",
     "references.md",
   ],
diff --git a/docs/src/printing.md b/docs/src/printing.md
new file mode 100644
index 0000000..b279057
--- /dev/null
+++ b/docs/src/printing.md
@@ -0,0 +1,13 @@
+```@meta
+CurrentModule = GenericCharacterTables
+DocTestSetup = :(using GenericCharacterTables, Oscar)
+```
+
+# Cyclotomic factorization
+
+Sometimes it can be interesting to see which cyclotomic polynomial divides a given
+polynomial representing a group or centralizer order.
+
+```@docs
+factor_cyclotomic
+```
diff --git a/src/CyclotomicFac.jl b/src/CyclotomicFac.jl
new file mode 100644
index 0000000..80ace23
--- /dev/null
+++ b/src/CyclotomicFac.jl
@@ -0,0 +1,111 @@
+@doc raw"""
+    CyclotomicFac{T <: PolyRingElem}
+
+The type representing a factorization into cyclotomic polynomials.
+"""
+struct CyclotomicFac{T <: PolyRingElem}
+  unit::T
+  fac::Dict{T, Int}
+  cyclo_fac::Dict{Int, Int}
+  function CyclotomicFac(f::Fac{T}) where T <: PolyRingElem
+    fac = Dict{T, Int}()
+    cyclo_fac = Dict{Int, Int}()
+    for (fact, expo) in f
+      b, n = is_cyclotomic_polynomial_with_data(fact)
+      if b
+        cyclo_fac[n] = expo
+      else
+        fac[fact] = expo
+      end
+    end
+    return new{T}(unit(f), fac, cyclo_fac)
+  end
+end
+
+function expressify(x::CyclotomicFac; context=nothing)
+  if length(x.fac) == length(x.cyclo_fac) == 0
+    return expressify(x.unit, context=context)
+  end
+  prod = Expr(:call, :*)
+  for (fact, expo) in x.fac
+    ep = expressify(fact, context=context)
+    if isone(expo)
+      push!(prod.args, ep)
+    else
+      push!(prod.args, Expr(:call, :^, ep, expo))
+    end
+  end
+  for (fact, expo) in x.cyclo_fac
+    ep = Symbol("Φ_$(fact)")
+    if isone(expo)
+      push!(prod.args, ep)
+    else
+      push!(prod.args, Expr(:call, :^, ep, expo))
+    end
+  end
+  return prod
+end
+
+@enable_all_show_via_expressify CyclotomicFac
+
+@doc raw"""
+    factor_cyclotomic(p::PolyRingElem)
+
+Return a factorization of `p` into cyclotomic polynomials. The main difference
+to `factor` is the way the result gets printed.
+# Examples
+```jldoctest
+julia> g=green_function_table("GL6");
+
+julia> factor_cyclotomic(order(g))
+q^15*Φ_5*Φ_4*Φ_6*Φ_2^3*Φ_3^2*Φ_1^6
+```
+"""
+factor_cyclotomic(p::PolyRingElem) = CyclotomicFac(factor(p))
+
+@doc raw"""
+    factor_cyclotomic(p::Generic.UnivPoly)
+
+Return a factorization of `p` into cyclotomic polynomials. The main difference
+to `factor` is the way the result gets printed.
+# Examples
+```jldoctest
+julia> g=generic_character_table("GL3");
+
+julia> factor_cyclotomic(order(g))
+q^3*Φ_2*Φ_3*Φ_1^3
+```
+"""
+factor_cyclotomic(p::Generic.UnivPoly) = CyclotomicFac(factor(to_univariate(p)))
+
+@doc raw"""
+    evaluate(a::CyclotomicFac)
+
+Multiply out the factorization into a single element.
+# Examples
+```jldoctest
+julia> Qx, x = QQ["x"];
+
+julia> f = x^16 - x^15 - x^14 + 2*x^11 - x^8 - x^7 + x^6;
+
+julia> fac = factor_cyclotomic(f)
+x^6*Φ_4*Φ_2^2*Φ_3*Φ_1^4
+
+julia> evaluate(fac)
+x^16 - x^15 - x^14 + 2*x^11 - x^8 - x^7 + x^6
+```
+"""
+function evaluate(a::CyclotomicFac)
+   r = a.unit
+   for (p, e) in Iterators.flatten((a.fac, Iterators.map(x -> cyclotomic_polynomial(first(x), parent(a.unit)) => last(x), a.cyclo_fac)))
+      r *= p^e
+   end
+   return r
+end
+
+@doc raw"""
+    length(a::CyclotomicFac)
+
+Return the number of factors of $a$, not including the unit.
+"""
+length(a::CyclotomicFac) = length(a.fac) + length(a.cyclo_fac)
diff --git a/src/Exports.jl b/src/Exports.jl
index 25fa5a1..5fe29fd 100644
--- a/src/Exports.jl
+++ b/src/Exports.jl
@@ -5,6 +5,7 @@ export conjugacy_class_type
 export conjugacy_class_type_index
 export degree
 export evaluate
+export factor_cyclotomic
 export generic_character_table
 export generic_cyclotomic_ring
 export green_function_table
diff --git a/src/GenericCharacterTables.jl b/src/GenericCharacterTables.jl
index 7c27f42..fee175b 100644
--- a/src/GenericCharacterTables.jl
+++ b/src/GenericCharacterTables.jl
@@ -36,5 +36,6 @@ include("Iteration.jl")
 include("Congruence.jl")
 include("Exports.jl")
 include("Deprecations.jl")
+include("CyclotomicFac.jl")
 
 end