diff --git a/src/sage/combinat/sf/kfpoly.py b/src/sage/combinat/sf/kfpoly.py
index 01470a24b5a..a5f09c7e0e7 100644
--- a/src/sage/combinat/sf/kfpoly.py
+++ b/src/sage/combinat/sf/kfpoly.py
@@ -24,6 +24,8 @@
 
 from sage.combinat.partition import _Partitions
 from sage.combinat.partitions import ZS1_iterator
+from sage.combinat.skew_partition import SkewPartitions
+from sage.combinat.skew_tableau import SemistandardSkewTableaux
 from sage.rings.polynomial.polynomial_ring import polygen
 from sage.rings.integer_ring import ZZ
 
@@ -31,17 +33,27 @@
 def KostkaFoulkesPolynomial(mu, nu, t=None):
     r"""
     Return the Kostka-Foulkes polynomial `K_{\mu, \nu}(t)`.
+    Here, `\mu` is a partition or a skew partition, whereas
+    `\nu` is a partition of the same size.
+
+    The Kostka-Foulkes polynomial is defined to be the sum
+    of the monomials `t^{\operatorname{charge}(T)}` over all
+    semistandard tableaux `T` of shape `\lambda / \mu``,
+    where `\operatorname{charge}(T)` denotes the charge
+    of the reading word of `T`
+    (see :meth:`sage.combinat.words.finite_word.FiniteWord_class.charge`).
 
     INPUT:
 
-    - ``mu``, ``nu`` -- partitions
+    - ``mu`` -- partition or skew partition
+    - ``nu`` -- partition
     - ``t`` -- an optional parameter (default: ``None``)
 
     OUTPUT:
 
-    - the Koskta-Foulkes polynomial indexed by partitions ``mu`` and ``nu`` and
+    - the Kostka-Foulkes polynomial indexed by ``mu`` and ``nu`` and
       evaluated at the parameter ``t``.  If ``t`` is ``None`` the resulting
-      polynomial is in the polynomial ring `\ZZ['t']`.
+      polynomial is in the polynomial ring `\ZZ[t]`.
 
     EXAMPLES::
 
@@ -56,13 +68,15 @@ def KostkaFoulkesPolynomial(mu, nu, t=None):
         sage: q = PolynomialRing(QQ,'q').gen()
         sage: KostkaFoulkesPolynomial([2,2],[2,1,1],q)
         q
+        sage: KostkaFoulkesPolynomial([[3,2],[1]],[2,2],q)
+        q + 1
 
     TESTS::
 
         sage: KostkaFoulkesPolynomial([2,4],[2,2])
         Traceback (most recent call last):
         ...
-        ValueError: mu must be a partition
+        ValueError: mu must be a partition or a skew partition
         sage: KostkaFoulkesPolynomial([2,2],[2,4])
         Traceback (most recent call last):
         ...
@@ -73,7 +87,9 @@ def KostkaFoulkesPolynomial(mu, nu, t=None):
         ValueError: mu and nu must be partitions of the same size
     """
     if mu not in _Partitions:
-        raise ValueError("mu must be a partition")
+        if mu in SkewPartitions():
+            return kfpoly_skew(mu, nu, t)
+        raise ValueError("mu must be a partition or a skew partition")
     if nu not in _Partitions:
         raise ValueError("nu must be a partition")
     if sum(mu) != sum(nu):
@@ -94,7 +110,7 @@ def kfpoly(mu, nu, t=None):
 
     OUTPUT:
 
-    - the Koskta-Foulkes polynomial indexed by partitions ``mu`` and ``nu`` and
+    - the Kostka-Foulkes polynomial indexed by partitions ``mu`` and ``nu`` and
       evaluated at the parameter ``t``.  If ``t`` is ``None`` the resulting polynomial
       is in the polynomial ring `\ZZ['t']`.
 
@@ -127,6 +143,47 @@ def kfpoly(mu, nu, t=None):
 
     return sum(f(rg) for rg in riggings(mu))
 
+def kfpoly_skew(lamu, nu, t=None):
+    r"""
+    Return the Kostka-Foulkes polynomial `K_{\lambda / \mu, \nu}(t)`
+    by summing `t^{\operatorname{charge}(T)}` over all semistandard
+    tableaux `T` of shape `\lambda / \mu``.
+
+    INPUT:
+
+    - ``lamu`` -- skew partition `\lambda / \mu`
+    - ``nu`` -- partition `\nu`
+    - ``t`` -- an optional parameter (default: ``None``)
+
+    OUTPUT:
+
+    - the Kostka-Foulkes polynomial indexed by ``mu`` and ``nu`` and
+      evaluated at the parameter ``t``.  If ``t`` is ``None`` the
+      resulting polynomial is in the polynomial ring `\ZZ['t']`.
+
+    EXAMPLES::
+
+        sage: from sage.combinat.sf.kfpoly import kfpoly_skew
+        sage: kfpoly_skew([[3,3], [1,1]], [2,1,1])
+        t
+        sage: kfpoly_skew([[3,3], [1,1]], [2,1,1], 5)
+        5
+        sage: kfpoly_skew([[5], [1]], [2,2])
+        t^2
+
+    TESTS::
+
+        sage: from sage.combinat.sf.kfpoly import kfpoly, kfpoly_skew
+        sage: all(kfpoly_skew(SkewPartition([b,[]]), c) == kfpoly(b,c)
+        ....:     for n in range(6) for b in Partitions(n)
+        ....:     for c in Partitions(n))
+        True
+    """
+    if t is None:
+        t = polygen(ZZ, 't')
+
+    return t.parent().sum(t ** (T.to_word().charge())
+                          for T in SemistandardSkewTableaux(lamu, nu))
 
 def schur_to_hl(mu, t=None):
     r"""