diff --git a/mathics/builtin/numbers/algebra.py b/mathics/builtin/numbers/algebra.py
index 1b366531dd..1ef7bafa55 100644
--- a/mathics/builtin/numbers/algebra.py
+++ b/mathics/builtin/numbers/algebra.py
@@ -58,6 +58,7 @@
     coeff_power,
     coefficient,
     eval_Apart,
+    eval_PolynomialMod,
     expand_polynomial,
     find_all_vars,
     get_exponents_sorted,
@@ -1387,6 +1388,41 @@ def eval(self, expr, evaluation: Evaluation, options: dict):
         return eval_Numerator(expr, evaluation, algebraic_options)
 
 
+class PolynomialMod(Builtin):
+    """
+    <url>:WMA link:
+    https://reference.wolfram.com/language/ref/PolynomialMod.html</url>
+
+    <dl>
+      <dt>'PolynomialMod'[$poly$, $m$]
+      <dd>for a positive integer $m$, gives $poly$ expanded, with all its \
+          coefficients reduced modulo $m$.
+    </dl>
+
+    >> PolynomialMod[3 x ^ 2 + 2 x + 1, 2]
+     = 1 + x ^ 2
+    >> PolynomialMod[3 x ^ 3 + 21 x ^ 2 - 7 x + 55, 9]
+     = 1 + 2 x + 3 x ^ 2 + 3 x ^ 3
+    >> PolynomialMod[1234 (x + y) ^ 7 - (x - y) ^ 5, 7]
+     = 6 x ^ 5 + 2 x ^ 7 + 5 x ^ 4 y + 4 x ^ 3 y ^ 2 + 3 x ^ 2 y ^ 3 + 2 x y ^ 4 + y ^ 5 + 2 y ^ 7
+
+    Parts that are not polynomials are treated as variables:
+    >> PolynomialMod[3 f[1] + 2 f[2] - f[3], 2]
+     = f[1] + f[3]
+
+    Reduction modulo a polynomial is not supported yet.
+    """
+
+    summary_text = "reduce the coefficients of a polynomial modulo an integer"
+
+    def eval(self, poly, m: Integer, evaluation: Evaluation):
+        "PolynomialMod[poly_, m_Integer]"
+        modulus = m.value
+        if modulus <= 0:
+            return None
+        return eval_PolynomialMod(poly, modulus)
+
+
 class PolynomialQ(Builtin):
     """
     <url>:Polynomial:
diff --git a/mathics/eval/numbers/algebra/polynomial.py b/mathics/eval/numbers/algebra/polynomial.py
index 126dc9ff2a..5b9f436efc 100644
--- a/mathics/eval/numbers/algebra/polynomial.py
+++ b/mathics/eval/numbers/algebra/polynomial.py
@@ -705,3 +705,50 @@ def get_exponents_sorted(expr, var) -> list:
             expos = [term.as_coeff_exponent(x)[1] for mul in muls for term in mul]
             result.add(from_sympy(tracing.run_sympy(sympy.Max, *[e for e in expos])))
     return sorted(result)
+
+
+def eval_PolynomialMod(poly: BaseElement, modulus: int) -> Optional[BaseElement]:
+    """
+    Evaluation routine for PolynomialMod[poly, m] with a positive integer m:
+    expand poly and reduce each integer coefficient to the range 0..m-1.
+
+    As in WMA, parts that are not polynomials, like f[1] or Sin[x], are
+    treated as variables. If a coefficient is not an integer (a Rational,
+    Real or Complex), None is returned and the expression is left
+    unevaluated.
+    """
+    expanded = expand_polynomial(poly, True, False, modulus=None, trig=False)
+    if expanded.has_form(SymbolList, None):
+        return None
+    terms = expanded.elements if expanded.has_form(SymbolPlus, None) else (expanded,)
+    reduced = []
+    for term in terms:
+        rest: Optional[BaseElement]
+        if isinstance(term, Number):
+            coefficient, rest = term, None
+        elif term.has_form(SymbolTimes, 2, None) and isinstance(
+            term.elements[0], Number
+        ):
+            coefficient = term.elements[0]
+            factors = term.elements[1:]
+            rest = (
+                factors[0] if len(factors) == 1 else Expression(SymbolTimes, *factors)
+            )
+        else:
+            coefficient, rest = Integer1, term
+        if not isinstance(coefficient, Integer):
+            return None
+        value = coefficient.value % modulus
+        if value == 0:
+            continue
+        if rest is None:
+            reduced.append(Integer(value))
+        elif value == 1:
+            reduced.append(rest)
+        else:
+            reduced.append(Expression(SymbolTimes, Integer(value), rest))
+    if not reduced:
+        return Integer0
+    if len(reduced) == 1:
+        return reduced[0]
+    return Expression(SymbolPlus, *reduced)
diff --git a/test/builtin/numbers/test_algebra.py b/test/builtin/numbers/test_algebra.py
index aebbb8e5dd..c40f654888 100644
--- a/test/builtin/numbers/test_algebra.py
+++ b/test/builtin/numbers/test_algebra.py
@@ -3,6 +3,7 @@
 Unit tests for mathics.builtins.numbers.algebra and
 mathics.builtins.numbers.integer
 """
+
 from test.helper import check_arg_counts, check_evaluation
 
 import pytest
@@ -388,6 +389,27 @@ def test_factor_terms_list():
         check_evaluation(str_expr, str_expected)
 
 
+def test_polynomialmod():
+    for str_expr, str_expected in [
+        ("PolynomialMod[3 x ^ 2 + 2 x + 1, 2]", "1 + x ^ 2"),
+        ("PolynomialMod[3 (x + 1) ^ 2, 2]", "1 + x ^ 2"),
+        ("PolynomialMod[-x, 2]", "x"),
+        ("PolynomialMod[-7, 3]", "2"),
+        ("PolynomialMod[2 x, 2]", "0"),
+        # Parts that are not polynomials are variables.
+        ("PolynomialMod[3 f[1] + 2 f[2], 2]", "f[1]"),
+        ("PolynomialMod[g[a, b], 2]", "g[a, b]"),
+        ("PolynomialMod[3 Sin[x] + 4 Sqrt[2] y, 3]", "Sqrt[2] y"),
+        # Not supported: left unevaluated rather than guessed.
+        ("PolynomialMod[x / 2, 3]", "PolynomialMod[x / 2, 3]"),
+        ("PolynomialMod[2.5 x, 2]", "PolynomialMod[2.5 x, 2]"),
+        ("PolynomialMod[x, -2]", "PolynomialMod[x, -2]"),
+        ("PolynomialMod[x ^ 2 + 1, x + 1]", "PolynomialMod[1 + x ^ 2, 1 + x]"),
+        ("PolynomialMod[{3 x}, 2]", "PolynomialMod[{3 x}, 2]"),
+    ]:
+        check_evaluation(str_expr, str_expected)
+
+
 def test_polynomialq():
     for str_expr, str_expected in [
         ("PolynomialQ[1/x[1]^2]", "False"),