diff --git a/mathics/builtin/numbers/linalg.py b/mathics/builtin/numbers/linalg.py
index 1a9da9992..d670c4727 100644
--- a/mathics/builtin/numbers/linalg.py
+++ b/mathics/builtin/numbers/linalg.py
@@ -124,10 +124,6 @@ class Eigenvalues(Builtin):
 
     >> Eigenvalues[{{7, 1}, {-4, 3}}]
      = {5, 5}
-
-    #> Eigenvalues[{{1, 0}, {0}}]
-     : Argument {{1, 0}, {0}} at position 1 is not a non-empty rectangular matrix.
-     = Eigenvalues[{{1, 0}, {0}}]
     """
 
     messages = {
@@ -221,9 +217,6 @@ class Eigenvectors(Builtin):
     >> Eigenvectors[{{0.1, 0.2}, {0.8, 0.5}}]
      = ...
     ### = {{-0.355518, -1.15048}, {-0.62896, 0.777438}}
-
-    #> Eigenvectors[{{-2, 1, -1}, {-3, 2, 1}, {-1, 1, 0}}]
-     = {{1, 7, 3}, {1, 1, 0}, {0, 0, 0}}
     """
 
     messages = {
@@ -365,18 +358,6 @@ class LeastSquares(Builtin):
     >> LeastSquares[{{1, 1, 1}, {1, 1, 2}}, {1, 3}]
      : Solving for underdetermined system not implemented.
      = LeastSquares[{{1, 1, 1}, {1, 1, 2}}, {1, 3}]
-
-    ## Inconsistent system - ideally we'd print a different message
-    #> LeastSquares[{{1, 1, 1}, {1, 1, 1}}, {1, 0}]
-     : Solving for underdetermined system not implemented.
-     = LeastSquares[{{1, 1, 1}, {1, 1, 1}}, {1, 0}]
-
-    #> LeastSquares[{1, {2}}, {1, 2}]
-     : Argument {1, {2}} at position 1 is not a non-empty rectangular matrix.
-     = LeastSquares[{1, {2}}, {1, 2}]
-    #> LeastSquares[{{1, 2}, {3, 4}}, {1, {2}}]
-     : Argument {1, {2}} at position 2 is not a non-empty rectangular matrix.
-     = LeastSquares[{{1, 2}, {3, 4}}, {1, {2}}]
     """
 
     messages = {
@@ -510,13 +491,6 @@ class LinearSolve(Builtin):
     >> LinearSolve[{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, {1, -2, 3}]
      : Linear equation encountered that has no solution.
      = LinearSolve[{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, {1, -2, 3}]
-
-    #> LinearSolve[{1, {2}}, {1, 2}]
-     : Argument {1, {2}} at position 1 is not a non-empty rectangular matrix.
-     = LinearSolve[{1, {2}}, {1, 2}]
-    #> LinearSolve[{{1, 2}, {3, 4}}, {1, {2}}]
-     : Argument {1, {2}} at position 2 is not a non-empty rectangular matrix.
-     = LinearSolve[{{1, 2}, {3, 4}}, {1, {2}}]
     """
 
     messages = {
@@ -582,13 +556,6 @@ class MatrixExp(Builtin):
 
     >> MatrixExp[{{1.5, 0.5}, {0.5, 2.0}}]
      = {{5.16266, 3.02952}, {3.02952, 8.19218}}
-
-    #> MatrixExp[{{a, 0}, {0, b}}]
-     = {{E ^ a, 0}, {0, E ^ b}}
-
-    #> MatrixExp[{{1, 0}, {0}}]
-     : Argument {{1, 0}, {0}} at position 1 is not a non-empty rectangular matrix.
-     = MatrixExp[{{1, 0}, {0}}]
     """
 
     messages = {
@@ -628,13 +595,6 @@ class MatrixPower(Builtin):
 
     >> MatrixPower[{{1, 2}, {2, 5}}, -3]
      = {{169, -70}, {-70, 29}}
-
-    #> MatrixPower[{{0, x}, {0, 0}}, n]
-     = MatrixPower[{{0, x}, {0, 0}}, n]
-
-    #> MatrixPower[{{1, 0}, {0}}, 2]
-     : Argument {{1, 0}, {0}} at position 1 is not a non-empty rectangular matrix.
-     = MatrixPower[{{1, 0}, {0}}, 2]
     """
 
     messages = {
@@ -681,10 +641,6 @@ class MatrixRank(Builtin):
      = 3
     >> MatrixRank[{{a, b}, {3 a, 3 b}}]
      = 1
-
-    #> MatrixRank[{{1, 0}, {0}}]
-     : Argument {{1, 0}, {0}} at position 1 is not a non-empty rectangular matrix.
-     = MatrixRank[{{1, 0}, {0}}]
     """
 
     messages = {
@@ -721,10 +677,6 @@ class NullSpace(Builtin):
      = {}
     >> MatrixRank[A]
      = 3
-
-    #> NullSpace[{1, {2}}]
-     : Argument {1, {2}} at position 1 is not a non-empty rectangular matrix.
-     = NullSpace[{1, {2}}]
     """
 
     messages = {
@@ -764,10 +716,6 @@ class PseudoInverse(Builtin):
 
     >> PseudoInverse[{{1.0, 2.5}, {2.5, 1.0}}]
      = {{-0.190476, 0.47619}, {0.47619, -0.190476}}
-
-    #> PseudoInverse[{1, {2}}]
-    : Argument {1, {2}} at position 1 is not a non-empty rectangular matrix.
-    = PseudoInverse[{1, {2}}]
     """
 
     messages = {
@@ -798,10 +746,6 @@ class QRDecomposition(Builtin):
 
     >> QRDecomposition[{{1, 2}, {3, 4}, {5, 6}}]
      = {{{Sqrt[35] / 35, 3 Sqrt[35] / 35, Sqrt[35] / 7}, {13 Sqrt[210] / 210, 2 Sqrt[210] / 105, -Sqrt[210] / 42}}, {{Sqrt[35], 44 Sqrt[35] / 35}, {0, 2 Sqrt[210] / 35}}}
-
-    #> QRDecomposition[{1, {2}}]
-     : Argument {1, {2}} at position 1 is not a non-empty rectangular matrix.
-     = QRDecomposition[{1, {2}}]
     """
 
     messages = {
@@ -844,10 +788,6 @@ class RowReduce(Builtin):
      . 0   1   2
      .
      . 0   0   0
-
-    #> RowReduce[{{1, 0}, {0}}]
-     : Argument {{1, 0}, {0}} at position 1 is not a non-empty rectangular matrix.
-     = RowReduce[{{1, 0}, {0}}]
     """
 
     messages = {
@@ -881,15 +821,6 @@ class SingularValueDecomposition(Builtin):
 
     >> SingularValueDecomposition[{{1.5, 2.0}, {2.5, 3.0}}]
      = {{{0.538954, 0.842335}, {0.842335, -0.538954}}, {{4.63555, 0.}, {0., 0.107862}}, {{0.628678, 0.777666}, {-0.777666, 0.628678}}}
-
-
-    #> SingularValueDecomposition[{{3/2, 2}, {5/2, 3}}]
-     : Symbolic SVD is not implemented, performing numerically.
-     = {{{0.538954, 0.842335}, {0.842335, -0.538954}}, {{4.63555, 0.}, {0., 0.107862}}, {{0.628678, 0.777666}, {-0.777666, 0.628678}}}
-
-    #> SingularValueDecomposition[{1, {2}}]
-     : Argument {1, {2}} at position 1 is not a non-empty rectangular matrix.
-     = SingularValueDecomposition[{1, {2}}]
     """
 
     # Sympy lacks symbolic SVD
diff --git a/mathics/builtin/numbers/numbertheory.py b/mathics/builtin/numbers/numbertheory.py
index 85687c3db..9e9972238 100644
--- a/mathics/builtin/numbers/numbertheory.py
+++ b/mathics/builtin/numbers/numbertheory.py
@@ -91,14 +91,6 @@ class Divisors(Builtin):
      = {1, 2, 4, 8, 11, 16, 22, 32, 44, 64, 88, 176, 352, 704}
     >> Divisors[{87, 106, 202, 305}]
      = {{1, 3, 29, 87}, {1, 2, 53, 106}, {1, 2, 101, 202}, {1, 5, 61, 305}}
-    #> Divisors[0]
-     = Divisors[0]
-    #> Divisors[{-206, -502, -1702, 9}]
-     = {{1, 2, 103, 206}, {1, 2, 251, 502}, {1, 2, 23, 37, 46, 74, 851, 1702}, {1, 3, 9}}
-    #> Length[Divisors[1000*369]]
-     = 96
-    #> Length[Divisors[305*176*369*100]]
-     = 672
     """
 
     # TODO: support GaussianIntegers
@@ -275,21 +267,6 @@ class FractionalPart(Builtin):
 
     >> FractionalPart[-5.25]
      = -0.25
-
-    #> FractionalPart[b]
-     = FractionalPart[b]
-
-    #> FractionalPart[{-2.4, -2.5, -3.0}]
-     = {-0.4, -0.5, 0.}
-
-    #> FractionalPart[14/32]
-     = 7 / 16
-
-    #> FractionalPart[4/(1 + 3 I)]
-     = 2 / 5 - I / 5
-
-    #> FractionalPart[Pi^20]
-     = -8769956796 + Pi ^ 20
     """
 
     attributes = A_LISTABLE | A_NUMERIC_FUNCTION | A_READ_PROTECTED | A_PROTECTED
@@ -370,47 +347,6 @@ class MantissaExponent(Builtin):
 
     >> MantissaExponent[10, b]
      = MantissaExponent[10, b]
-
-    #> MantissaExponent[E, Pi]
-     = {E / Pi, 1}
-
-    #> MantissaExponent[Pi, Pi]
-     = {1 / Pi, 2}
-
-    #> MantissaExponent[5/2 + 3, Pi]
-     = {11 / (2 Pi ^ 2), 2}
-
-    #> MantissaExponent[b]
-     = MantissaExponent[b]
-
-    #> MantissaExponent[17, E]
-     = {17 / E ^ 3, 3}
-
-    #> MantissaExponent[17., E]
-     = {0.84638, 3}
-
-    #> MantissaExponent[Exp[Pi], 2]
-     = {E ^ Pi / 32, 5}
-
-    #> MantissaExponent[3 + 2 I, 2]
-     : The value 3 + 2 I is not a real number
-     = MantissaExponent[3 + 2 I, 2]
-
-    #> MantissaExponent[25, 0.4]
-     : Base 0.4 is not a real number greater than 1.
-     = MantissaExponent[25, 0.4]
-
-    #> MantissaExponent[0.0000124]
-     = {0.124, -4}
-
-    #> MantissaExponent[0.0000124, 2]
-     = {0.812646, -16}
-
-    #> MantissaExponent[0]
-     = {0, 0}
-
-    #> MantissaExponent[0, 2]
-     = {0, 0}
     """
 
     attributes = A_LISTABLE | A_PROTECTED
@@ -674,9 +610,6 @@ class PrimePowerQ(Builtin):
 
     >> PrimePowerQ[371293]
      = True
-
-    #> PrimePowerQ[1]
-     = False
     """
 
     attributes = A_LISTABLE | A_PROTECTED | A_READ_PROTECTED
@@ -687,19 +620,19 @@ class PrimePowerQ(Builtin):
 
     # TODO: GaussianIntegers option
     """
-    #> PrimePowerQ[5, GaussianIntegers -> True]
+    ##> PrimePowerQ[5, GaussianIntegers -> True]
      = False
     """
 
     # TODO: Complex args
     """
-    #> PrimePowerQ[{3 + I, 3 - 2 I, 3 + 4 I, 9 + 7 I}]
+    ##> PrimePowerQ[{3 + I, 3 - 2 I, 3 + 4 I, 9 + 7 I}]
      = {False, True, True, False}
     """
 
     # TODO: Gaussian rationals
     """
-    #> PrimePowerQ[2/125 - 11 I/125]
+    ##> PrimePowerQ[2/125 - 11 I/125]
      = True
     """
 
@@ -740,12 +673,6 @@ class RandomPrime(Builtin):
 
     >> RandomPrime[{10,30}, {2,5}]
      = ...
-
-    #> RandomPrime[{10,12}, {2,2}]
-     = {{11, 11}, {11, 11}}
-
-    #> RandomPrime[2, {3,2}]
-     = {{2, 2}, {2, 2}, {2, 2}}
     """
 
     messages = {
diff --git a/mathics/builtin/specialfns/bessel.py b/mathics/builtin/specialfns/bessel.py
index e2e2455da..ebe20fea2 100644
--- a/mathics/builtin/specialfns/bessel.py
+++ b/mathics/builtin/specialfns/bessel.py
@@ -24,7 +24,6 @@
 
 
 class _Bessel(MPMathFunction):
-
     attributes = A_LISTABLE | A_NUMERIC_FUNCTION | A_PROTECTED | A_READ_PROTECTED
 
     nargs = {2}
@@ -109,18 +108,6 @@ class AiryAiZero(Builtin):
 
     >> N[AiryAiZero[1]]
      = -2.33811
-
-    #> AiryAiZero[1]
-     = AiryAiZero[1]
-
-    #> AiryAiZero[1.]
-     = AiryAiZero[1.]
-
-    #> AiryAi[AiryAiZero[1]]
-     = 0
-
-    #> N[AiryAiZero[2], 100]
-     = -4.087949444130970616636988701457391060224764699108529754984160876025121946836047394331169160758270562
     """
 
     # TODO: 'AiryAiZero[$k$, $x0$]' - $k$th zero less than x0
@@ -235,18 +222,6 @@ class AiryBiZero(Builtin):
 
     >> N[AiryBiZero[1]]
      = -1.17371
-
-    #> AiryBiZero[1]
-     = AiryBiZero[1]
-
-    #> AiryBiZero[1.]
-     = AiryBiZero[1.]
-
-    #> AiryBi[AiryBiZero[1]]
-     = 0
-
-    #> N[AiryBiZero[2], 100]
-     = -3.271093302836352715680228240166413806300935969100284801485032396261130864238742879252000673830055014
     """
 
     # TODO: 'AiryBiZero[$k$, $x0$]' - $k$th zero less than x0
@@ -380,9 +355,6 @@ class BesselJ(_Bessel):
     >> BesselJ[0, 5.2]
      = -0.11029
 
-    #> BesselJ[2.5, 1]
-     = 0.0494968
-
     >> D[BesselJ[n, z], z]
      = -BesselJ[1 + n, z] / 2 + BesselJ[-1 + n, z] / 2
 
diff --git a/mathics/builtin/specialfns/gamma.py b/mathics/builtin/specialfns/gamma.py
index 33746fd11..7369e0ba6 100644
--- a/mathics/builtin/specialfns/gamma.py
+++ b/mathics/builtin/specialfns/gamma.py
@@ -155,8 +155,6 @@ class Factorial(PostfixOperator, MPMathFunction):
     >> !a! //FullForm
      = Not[Factorial[a]]
 
-    #> 0!
-     = 1
     """
 
     attributes = A_NUMERIC_FUNCTION | A_PROTECTED
@@ -301,22 +299,6 @@ class Gamma(MPMathMultiFunction):
     Both 'Gamma' and 'Factorial' functions are continuous:
     >> Plot[{Gamma[x], x!}, {x, 0, 4}]
      = -Graphics-
-
-    ## Issue 203
-    #> N[Gamma[24/10], 100]
-     = 1.242169344504305404913070252268300492431517240992022966055507541481863694148882652446155342679460339
-    #> N[N[Gamma[24/10],100]/N[Gamma[14/10],100],100]
-     = 1.400000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
-    #> % // Precision
-     = 100.
-
-    #> Gamma[1.*^20]
-     : Overflow occurred in computation.
-     = Overflow[]
-
-    ## Needs mpmath support for lowergamma
-    #> Gamma[1., 2.]
-     = Gamma[1., 2.]
     """
 
     mpmath_names = {
diff --git a/mathics/builtin/specialfns/orthogonal.py b/mathics/builtin/specialfns/orthogonal.py
index 464b148f9..bd2cb002e 100644
--- a/mathics/builtin/specialfns/orthogonal.py
+++ b/mathics/builtin/specialfns/orthogonal.py
@@ -269,9 +269,6 @@ class SphericalHarmonicY(MPMathFunction):
     ## Results depend on sympy version
     >> SphericalHarmonicY[3, 1, theta, phi]
      = ...
-
-    #> SphericalHarmonicY[1,1,x,y]
-     = -Sqrt[6] E ^ (I y) Sin[x] / (4 Sqrt[Pi])
     """
 
     nargs = {4}
diff --git a/test/builtin/numbers/test_linalg.py b/test/builtin/numbers/test_linalg.py
index 73a914082..3494e65e2 100644
--- a/test/builtin/numbers/test_linalg.py
+++ b/test/builtin/numbers/test_linalg.py
@@ -88,3 +88,139 @@ def test_inverse(str_expr, str_expected, fail_msg, warnings):
     check_evaluation(
         str_expr, str_expected, failure_message="", expected_messages=warnings
     )
+
+
+@pytest.mark.parametrize(
+    ("str_expr", "msgs", "str_expected", "fail_msg"),
+    [
+        (
+            "Eigenvalues[{{1, 0}, {0}}]",
+            (
+                "Argument {{1, 0}, {0}} at position 1 is not a non-empty rectangular matrix.",
+            ),
+            "Eigenvalues[{{1, 0}, {0}}]",
+            None,
+        ),
+        (
+            "Eigenvectors[{{-2, 1, -1}, {-3, 2, 1}, {-1, 1, 0}}]",
+            None,
+            "{{1, 7, 3}, {1, 1, 0}, {0, 0, 0}}",
+            None,
+        ),
+        ## Inconsistent system - ideally we'd print a different message
+        (
+            "LeastSquares[{{1, 1, 1}, {1, 1, 1}}, {1, 0}]",
+            ("Solving for underdetermined system not implemented.",),
+            "LeastSquares[{{1, 1, 1}, {1, 1, 1}}, {1, 0}]",
+            None,
+        ),
+        (
+            "LeastSquares[{1, {2}}, {1, 2}]",
+            ("Argument {1, {2}} at position 1 is not a non-empty rectangular matrix.",),
+            "LeastSquares[{1, {2}}, {1, 2}]",
+            None,
+        ),
+        (
+            "LeastSquares[{{1, 2}, {3, 4}}, {1, {2}}]",
+            ("Argument {1, {2}} at position 2 is not a non-empty rectangular matrix.",),
+            "LeastSquares[{{1, 2}, {3, 4}}, {1, {2}}]",
+            None,
+        ),
+        (
+            "LinearSolve[{1, {2}}, {1, 2}]",
+            ("Argument {1, {2}} at position 1 is not a non-empty rectangular matrix.",),
+            "LinearSolve[{1, {2}}, {1, 2}]",
+            None,
+        ),
+        (
+            "LinearSolve[{{1, 2}, {3, 4}}, {1, {2}}]",
+            ("Argument {1, {2}} at position 2 is not a non-empty rectangular matrix.",),
+            "LinearSolve[{{1, 2}, {3, 4}}, {1, {2}}]",
+            None,
+        ),
+        ("MatrixExp[{{a, 0}, {0, b}}]", None, "{{E ^ a, 0}, {0, E ^ b}}", None),
+        (
+            "MatrixExp[{{1, 0}, {0}}]",
+            (
+                "Argument {{1, 0}, {0}} at position 1 is not a non-empty rectangular matrix.",
+            ),
+            "MatrixExp[{{1, 0}, {0}}]",
+            None,
+        ),
+        (
+            "MatrixPower[{{0, x}, {0, 0}}, n]",
+            None,
+            "MatrixPower[{{0, x}, {0, 0}}, n]",
+            None,
+        ),
+        (
+            "MatrixPower[{{1, 0}, {0}}, 2]",
+            (
+                "Argument {{1, 0}, {0}} at position 1 is not a non-empty rectangular matrix.",
+            ),
+            "MatrixPower[{{1, 0}, {0}}, 2]",
+            None,
+        ),
+        (
+            "MatrixRank[{{1, 0}, {0}}]",
+            (
+                "Argument {{1, 0}, {0}} at position 1 is not a non-empty rectangular matrix.",
+            ),
+            "MatrixRank[{{1, 0}, {0}}]",
+            None,
+        ),
+        (
+            "NullSpace[{1, {2}}]",
+            ("Argument {1, {2}} at position 1 is not a non-empty rectangular matrix.",),
+            "NullSpace[{1, {2}}]",
+            None,
+        ),
+        (
+            "PseudoInverse[{1, {2}}]",
+            ("Argument {1, {2}} at position 1 is not a non-empty rectangular matrix.",),
+            "PseudoInverse[{1, {2}}]",
+            None,
+        ),
+        (
+            "QRDecomposition[{1, {2}}]",
+            ("Argument {1, {2}} at position 1 is not a non-empty rectangular matrix.",),
+            "QRDecomposition[{1, {2}}]",
+            None,
+        ),
+        (
+            "RowReduce[{{1, 0}, {0}}]",
+            (
+                "Argument {{1, 0}, {0}} at position 1 is not a non-empty rectangular matrix.",
+            ),
+            "RowReduce[{{1, 0}, {0}}]",
+            None,
+        ),
+        (
+            "SingularValueDecomposition[{{3/2, 2}, {5/2, 3}}]",
+            ("Symbolic SVD is not implemented, performing numerically.",),
+            (
+                "{{{0.538954, 0.842335}, {0.842335, -0.538954}}, "
+                "{{4.63555, 0.}, {0., 0.107862}}, "
+                "{{0.628678, 0.777666}, {-0.777666, 0.628678}}}"
+            ),
+            None,
+        ),
+        (
+            "SingularValueDecomposition[{1, {2}}]",
+            ("Argument {1, {2}} at position 1 is not a non-empty rectangular matrix.",),
+            "SingularValueDecomposition[{1, {2}}]",
+            None,
+        ),
+    ],
+)
+def test_private_doctests_linalg(str_expr, msgs, str_expected, fail_msg):
+    """ """
+    check_evaluation(
+        str_expr,
+        str_expected,
+        to_string_expr=True,
+        to_string_expected=True,
+        hold_expected=True,
+        failure_message=fail_msg,
+        expected_messages=msgs,
+    )
diff --git a/test/builtin/numbers/test_numbertheory.py b/test/builtin/numbers/test_numbertheory.py
new file mode 100644
index 000000000..8e5e149b4
--- /dev/null
+++ b/test/builtin/numbers/test_numbertheory.py
@@ -0,0 +1,70 @@
+# -*- coding: utf-8 -*-
+"""
+Unit tests for mathics.builtins.numbers.numbertheory
+"""
+from test.helper import check_evaluation
+
+import pytest
+
+
+@pytest.mark.parametrize(
+    ("str_expr", "msgs", "str_expected", "fail_msg"),
+    [
+        ("Divisors[0]", None, "Divisors[0]", None),
+        (
+            "Divisors[{-206, -502, -1702, 9}]",
+            None,
+            (
+                "{{1, 2, 103, 206}, "
+                "{1, 2, 251, 502}, "
+                "{1, 2, 23, 37, 46, 74, 851, 1702}, "
+                "{1, 3, 9}}"
+            ),
+            None,
+        ),
+        ("Length[Divisors[1000*369]]", None, "96", None),
+        ("Length[Divisors[305*176*369*100]]", None, "672", None),
+        ("FractionalPart[b]", None, "FractionalPart[b]", None),
+        ("FractionalPart[{-2.4, -2.5, -3.0}]", None, "{-0.4, -0.5, 0.}", None),
+        ("FractionalPart[14/32]", None, "7 / 16", None),
+        ("FractionalPart[4/(1 + 3 I)]", None, "2 / 5 - I / 5", None),
+        ("FractionalPart[Pi^20]", None, "-8769956796 + Pi ^ 20", None),
+        ("MantissaExponent[E, Pi]", None, "{E / Pi, 1}", None),
+        ("MantissaExponent[Pi, Pi]", None, "{1 / Pi, 2}", None),
+        ("MantissaExponent[5/2 + 3, Pi]", None, "{11 / (2 Pi ^ 2), 2}", None),
+        ("MantissaExponent[b]", None, "MantissaExponent[b]", None),
+        ("MantissaExponent[17, E]", None, "{17 / E ^ 3, 3}", None),
+        ("MantissaExponent[17., E]", None, "{0.84638, 3}", None),
+        ("MantissaExponent[Exp[Pi], 2]", None, "{E ^ Pi / 32, 5}", None),
+        (
+            "MantissaExponent[3 + 2 I, 2]",
+            ("The value 3 + 2 I is not a real number",),
+            "MantissaExponent[3 + 2 I, 2]",
+            None,
+        ),
+        (
+            "MantissaExponent[25, 0.4]",
+            ("Base 0.4 is not a real number greater than 1.",),
+            "MantissaExponent[25, 0.4]",
+            None,
+        ),
+        ("MantissaExponent[0.0000124]", None, "{0.124, -4}", None),
+        ("MantissaExponent[0.0000124, 2]", None, "{0.812646, -16}", None),
+        ("MantissaExponent[0]", None, "{0, 0}", None),
+        ("MantissaExponent[0, 2]", None, "{0, 0}", None),
+        ("PrimePowerQ[1]", None, "False", None),
+        ("RandomPrime[{10,12}, {2,2}]", None, "{{11, 11}, {11, 11}}", None),
+        ("RandomPrime[2, {3,2}]", None, "{{2, 2}, {2, 2}, {2, 2}}", None),
+    ],
+)
+def test_private_doctests_numbertheory(str_expr, msgs, str_expected, fail_msg):
+    """ """
+    check_evaluation(
+        str_expr,
+        str_expected,
+        to_string_expr=True,
+        to_string_expected=True,
+        hold_expected=True,
+        failure_message=fail_msg,
+        expected_messages=msgs,
+    )
diff --git a/test/builtin/specialfns/test_bessel.py b/test/builtin/specialfns/test_bessel.py
index b6201d76e..6b7296763 100644
--- a/test/builtin/specialfns/test_bessel.py
+++ b/test/builtin/specialfns/test_bessel.py
@@ -1,6 +1,7 @@
 # -*- coding: utf-8 -*-
 """
-Unit tests for mathics.builtins.arithmetic.bessel
+Unit tests for mathics.builtins.specialfns.bessel and
+mathics.builtins.specialfns.orthogonal
 """
 from test.helper import check_evaluation
 
@@ -30,3 +31,64 @@ def test_add(str_expr, str_expected, assert_failure_msg):
     check_evaluation(
         str_expr, str_expected, hold_expected=True, failure_message=assert_failure_msg
     )
+
+
+@pytest.mark.parametrize(
+    ("str_expr", "msgs", "str_expected", "fail_msg"),
+    [
+        ("AiryAiZero[1]", None, "AiryAiZero[1]", None),
+        ("AiryAiZero[1.]", None, "AiryAiZero[1.]", None),
+        ("AiryAi[AiryAiZero[1]]", None, "0", None),
+        (
+            "N[AiryAiZero[2], 100]",
+            None,
+            "-4.087949444130970616636988701457391060224764699108529754984160876025121946836047394331169160758270562",
+            None,
+        ),
+        ("AiryBiZero[1]", None, "AiryBiZero[1]", None),
+        ("AiryBiZero[1.]", None, "AiryBiZero[1.]", None),
+        ("AiryBi[AiryBiZero[1]]", None, "0", None),
+        (
+            "N[AiryBiZero[2], 100]",
+            None,
+            "-3.271093302836352715680228240166413806300935969100284801485032396261130864238742879252000673830055014",
+            None,
+        ),
+        ("BesselJ[2.5, 1]", None, "0.0494968", None),
+    ],
+)
+def test_private_doctests_bessel(str_expr, msgs, str_expected, fail_msg):
+    """ """
+    check_evaluation(
+        str_expr,
+        str_expected,
+        to_string_expr=True,
+        to_string_expected=True,
+        hold_expected=True,
+        failure_message=fail_msg,
+        expected_messages=msgs,
+    )
+
+
+@pytest.mark.parametrize(
+    ("str_expr", "msgs", "str_expected", "fail_msg"),
+    [
+        (
+            "SphericalHarmonicY[1,1,x,y]",
+            None,
+            "-Sqrt[6] E ^ (I y) Sin[x] / (4 Sqrt[Pi])",
+            None,
+        ),
+    ],
+)
+def test_private_doctests_orthogonal(str_expr, msgs, str_expected, fail_msg):
+    """ """
+    check_evaluation(
+        str_expr,
+        str_expected,
+        to_string_expr=True,
+        to_string_expected=True,
+        hold_expected=True,
+        failure_message=fail_msg,
+        expected_messages=msgs,
+    )
diff --git a/test/builtin/specialfns/test_gamma.py b/test/builtin/specialfns/test_gamma.py
new file mode 100644
index 000000000..b5d1d4148
--- /dev/null
+++ b/test/builtin/specialfns/test_gamma.py
@@ -0,0 +1,46 @@
+# -*- coding: utf-8 -*-
+"""
+Unit tests for mathics.builtins.specialfns.gamma
+"""
+from test.helper import check_evaluation
+
+import pytest
+
+
+@pytest.mark.parametrize(
+    ("str_expr", "msgs", "str_expected", "fail_msg"),
+    [
+        ("0!", None, "1", None),
+        (
+            "N[Gamma[24/10], 100]",
+            None,
+            "1.242169344504305404913070252268300492431517240992022966055507541481863694148882652446155342679460339",
+            "Issue 203",
+        ),
+        (
+            "res=N[N[Gamma[24/10],100]/N[Gamma[14/10],100],100]",
+            None,
+            "1.400000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000",
+            "Issue 203",
+        ),
+        ("res // Precision", None, "100.", None),
+        (
+            "Gamma[1.*^20]",
+            ("Overflow occurred in computation.",),
+            "Overflow[]",
+            "Overflow",
+        ),
+        ("Gamma[1., 2.]", None, "Gamma[1., 2.]", "needs mpmath for lowergamma"),
+    ],
+)
+def test_private_doctests_gamma(str_expr, msgs, str_expected, fail_msg):
+    """ """
+    check_evaluation(
+        str_expr,
+        str_expected,
+        to_string_expr=True,
+        to_string_expected=True,
+        hold_expected=True,
+        failure_message=fail_msg,
+        expected_messages=msgs,
+    )