Repository: sympy/sympy

Issue:
diophantine returns a non-integer solution
The following code:

    import sympy
    from sympy.solvers.diophantine import diophantine
    sympy.simplify(diophantine('-48 - 2 * m * (3*m - 1) + n * (3*n - 1)',
                               param=sympy.Symbol('t')).pop()[1].replace('t', 1))

returns -71816005847/3 on 0.7.6

Discussion:

Patch:
diff --git a/AUTHORS b/AUTHORS
index 07e0a52f59..47a4f03692 100644
--- a/AUTHORS
+++ b/AUTHORS
@@ -387,10 +387,5 @@ harshil goel <harshil158@gmail.com>
 Lokesh Sharma <lokeshhsharma@gmail.com>
 Sartaj Singh <singhsartaj94@gmail.com>
 Chris Swierczewski <cswiercz@gmail.com>
-Sumith <sumith1896@gmail.com>
 Konstantin Togoi <togoi.konstantin@outlook.com>
 Param Singh <paramsingh258@gmail.com>
-Philippe Bouafia <philippe.bouafia@ensea.fr>
-Lucas Jones <lucas@lucasjones.co.uk>
-Peter Schmidt <peter@peterjs.com>
-Jiaxing Liang <liangjiaxing57@gmail.com>
diff --git a/doc/src/aboutus.rst b/doc/src/aboutus.rst
index 1f02eff625..66cdf0e007 100644
--- a/doc/src/aboutus.rst
+++ b/doc/src/aboutus.rst
@@ -392,13 +392,8 @@ want to be mentioned here, so see our repository history for a full list).
 #. Lokesh Sharma: reorder setup.py file imports to correct NameError
 #. Sartaj Singh: use \left \| instead of \lvert for latex Abs
 #. Chris Swierczewski: RootOf.evalf watches for root on interval boundary
-#. Sumith: implemented columnspace of a matrix
 #. Konstantin Togoi: added truth_table and binary digit utilities to boolalg
 #. Param Singh: added missing units from mks system
-#. Philippe Bouafia: Fixed a typo.
-#. Lucas Jones: Add a test in integrals
-#. Peter Schmidt: Removed duplicate tests
-#. Jiaxing Liang: Fix infinite recursion issue in test_f
 
 Up-to-date list in the order of the first contribution is given in the `AUTHORS
 <https://github.com/sympy/sympy/blob/master/AUTHORS>`_ file.
diff --git a/doc/src/modules/integrals/integrals.rst b/doc/src/modules/integrals/integrals.rst
index a5fb8cf471..2eb43ac901 100644
--- a/doc/src/modules/integrals/integrals.rst
+++ b/doc/src/modules/integrals/integrals.rst
@@ -93,6 +93,13 @@ API reference
 
 .. autofunction:: sympy.integrals.integrate
 .. autofunction:: sympy.integrals.line_integrate
+.. autofunction:: sympy.integrals.deltaintegrate
+.. autofunction:: sympy.integrals.ratint
+.. autofunction:: sympy.integrals.heurisch
+.. autofunction:: sympy.integrals.heurisch_wrapper
+.. autofunction:: sympy.integrals.trigintegrate
+.. autofunction:: sympy.integrals.manualintegrate
+.. autofunction:: sympy.integrals.manualintegrate.integral_steps
 
 The class `Integral` represents an unevaluated integral and has some methods that help in the integration of an expression.
 
diff --git a/doc/src/modules/solvers/diophantine.rst b/doc/src/modules/solvers/diophantine.rst
index 102046b987..0129975280 100644
--- a/doc/src/modules/solvers/diophantine.rst
+++ b/doc/src/modules/solvers/diophantine.rst
@@ -164,14 +164,17 @@ solutions. Consider the below cases where `\Delta = 8`.
 
 >>> diophantine(x**2 - 4*x*y + 2*y**2 - 3*x + 7*y - 5)
 set()
->>> from sympy import expand
+>>> from sympy import sqrt
 >>> n = symbols("n", integer=True)
 >>> s = diophantine(x**2 -  2*y**2 - 2*x - 4*y, n)
->>> x_n, y_n = s.pop()
->>> expand(x_n)
--(-2*sqrt(2) + 3)**n/2 + sqrt(2)*(-2*sqrt(2) + 3)**n/2 - sqrt(2)*(2*sqrt(2) + 3)**n/2 - (2*sqrt(2) + 3)**n/2 + 1
->>> expand(y_n)
--sqrt(2)*(-2*sqrt(2) + 3)**n/4 + (-2*sqrt(2) + 3)**n/2 + sqrt(2)*(2*sqrt(2) + 3)**n/4 + (2*sqrt(2) + 3)**n/2 - 1
+>>> x_1, y_1 = s.pop()
+>>> x_2, y_2 = s.pop()
+>>> x_n = -(-2*sqrt(2) + 3)**n/2 + sqrt(2)*(-2*sqrt(2) + 3)**n/2 - sqrt(2)*(2*sqrt(2) + 3)**n/2 - (2*sqrt(2) + 3)**n/2 + 1
+>>> x_1 == x_n or x_2 == x_n
+True
+>>> y_n = -sqrt(2)*(-2*sqrt(2) + 3)**n/4 + (-2*sqrt(2) + 3)**n/2 + sqrt(2)*(2*sqrt(2) + 3)**n/4 + (2*sqrt(2) + 3)**n/2 - 1
+>>> y_1 == y_n or y_2 == y_n
+True
 
 Here `n` is an integer. Although x_n and y_n may not look like
 integers, substituting in specific values for n (and simplifying) shows that they
diff --git a/examples/galgebra/physics_check_latex.py b/examples/galgebra/physics_check_latex.py
index 3b69622d9c..e7dda69f59 100755
--- a/examples/galgebra/physics_check_latex.py
+++ b/examples/galgebra/physics_check_latex.py
@@ -28,7 +28,7 @@ def Maxwells_Equations_in_Geometric_Calculus():
     print('\\text{Electromagnetic Field Bi-Vector\\;\\;} F = E+IB =', F)
     print('%\\text{Four Current Density\\;\\;} J =', J)
     gradF = grad*F
-    print('#Geometric Derivative of Electromagnetic Field Bi-Vector')
+    print('#Geometric Derivative of Electomagnetic Field Bi-Vector')
     gradF.Fmt(3, 'grad*F')
 
     print('#Maxwell Equations')
diff --git a/examples/intermediate/coupled_cluster.py b/examples/intermediate/coupled_cluster.py
index eebda39f4c..5f5ab28b92 100755
--- a/examples/intermediate/coupled_cluster.py
+++ b/examples/intermediate/coupled_cluster.py
@@ -88,7 +88,7 @@ def main():
     comm4 = evaluate_deltas(comm4)
     comm4 = substitute_dummies(comm4)
 
-    print("construct Hausdorff expansion...")
+    print("construct Hausdoff expansion...")
     eq = H + comm1 + comm2/2 + comm3/6 + comm4/24
     eq = eq.expand()
     eq = evaluate_deltas(eq)
diff --git a/sympy/core/expr.py b/sympy/core/expr.py
index 0de02ff424..14142321ad 100644
--- a/sympy/core/expr.py
+++ b/sympy/core/expr.py
@@ -682,6 +682,10 @@ def equals(self, other, failing_expression=False):
             return diff
         return None
 
+    def _eval_is_composite(self):
+        if self.is_integer and self.is_positive and self.is_prime is False:
+            return True
+
     def _eval_is_positive(self):
         from sympy.polys.numberfields import minimal_polynomial
         from sympy.polys.polyerrors import NotAlgebraic
diff --git a/sympy/core/exprtools.py b/sympy/core/exprtools.py
index 654b7f3bf4..7e786424f3 100644
--- a/sympy/core/exprtools.py
+++ b/sympy/core/exprtools.py
@@ -72,7 +72,7 @@ def _monotonic_sign(self):
         rv = _monotonic_sign(-self)
         return rv if rv is None else -rv
 
-    if not self.is_Add and self.as_numer_denom()[1].is_number:
+    if self.is_Symbol:
         s = self
         if s.is_prime:
             if s.is_odd:
diff --git a/sympy/core/facts.py b/sympy/core/facts.py
index 4cc36ac997..ef4af6b5ad 100644
--- a/sympy/core/facts.py
+++ b/sympy/core/facts.py
@@ -186,7 +186,7 @@ def apply_beta_to_alpha_route(alpha_implications, beta_rules):
                     ximpls.add(bimpl)
 
                     # we introduced new implication - now we have to restore
-                    # completeness of the whole set.
+                    # completness of the whole set.
                     bimpl_impl = x_impl.get(bimpl)
                     if bimpl_impl is not None:
                         ximpls |= bimpl_impl[0]
diff --git a/sympy/functions/elementary/exponential.py b/sympy/functions/elementary/exponential.py
index 01559e77d8..5bc5af8089 100644
--- a/sympy/functions/elementary/exponential.py
+++ b/sympy/functions/elementary/exponential.py
@@ -301,7 +301,7 @@ def taylor_term(n, x, *previous_terms):
             p = previous_terms[-1]
             if p is not None:
                 return p * x / n
-        return x**n/factorial(n)
+        return x**n/factorial()(n)
 
     def as_real_imag(self, deep=True, **hints):
         """
@@ -687,9 +687,6 @@ def _eval_is_positive(self):
     def _eval_is_zero(self):
         return (self.args[0] - 1).is_zero
 
-    def _eval_is_nonnegative(self):
-        return (self.args[0] - 1).is_nonnegative
-
     def _eval_nseries(self, x, n, logx):
         # NOTE Please see the comment at the beginning of this file, labelled
         #      IMPORTANT.
diff --git a/sympy/functions/elementary/trigonometric.py b/sympy/functions/elementary/trigonometric.py
index df0117467b..4ad9068349 100644
--- a/sympy/functions/elementary/trigonometric.py
+++ b/sympy/functions/elementary/trigonometric.py
@@ -1942,9 +1942,6 @@ def _eval_is_rational(self):
     def _eval_is_positive(self):
         return self.args[0].is_positive
 
-    def _eval_is_nonnegative(self):
-        return self.args[0].is_nonnegative
-
     @classmethod
     def eval(cls, arg):
         if arg.is_Number:
diff --git a/sympy/integrals/__init__.py b/sympy/integrals/__init__.py
index 139dd599a3..46c0dcd31a 100644
--- a/sympy/integrals/__init__.py
+++ b/sympy/integrals/__init__.py
@@ -9,7 +9,11 @@
     >>> integrate(sin(x),x)
     -cos(x)
 """
+from .deltafunctions import deltaintegrate
+from .heurisch import heurisch, heurisch_wrapper
 from .integrals import integrate, Integral, line_integrate
+from .manualintegrate import manualintegrate, integral_steps
+from .rationaltools import ratint
 from .transforms import (mellin_transform, inverse_mellin_transform,
                         MellinTransform, InverseMellinTransform,
                         laplace_transform, inverse_laplace_transform,
@@ -22,3 +26,4 @@
                         CosineTransform, InverseCosineTransform,
                         hankel_transform, inverse_hankel_transform,
                         HankelTransform, InverseHankelTransform)
+from .trigonometry import trigintegrate
diff --git a/sympy/matrices/expressions/matexpr.py b/sympy/matrices/expressions/matexpr.py
index b07b3750ad..01b7420823 100644
--- a/sympy/matrices/expressions/matexpr.py
+++ b/sympy/matrices/expressions/matexpr.py
@@ -322,14 +322,6 @@ class MatrixElement(Expr):
     j = property(lambda self: self.args[2])
     _diff_wrt = True
 
-    def doit(self, **kwargs):
-        deep = kwargs.get('deep', True)
-        if deep:
-            args = [arg.doit(**kwargs) for arg in self.args]
-        else:
-            args = self.args
-        return args[0][args[1], args[2]]
-
 
 class MatrixSymbol(MatrixExpr):
     """Symbolic representation of a Matrix object
diff --git a/sympy/physics/quantum/qubit.py b/sympy/physics/quantum/qubit.py
index 216f9d556f..58ba60d675 100644
--- a/sympy/physics/quantum/qubit.py
+++ b/sympy/physics/quantum/qubit.py
@@ -481,7 +481,7 @@ def matrix_to_density(mat):
 
 
 def qubit_to_matrix(qubit, format='sympy'):
-    """Converts an Add/Mul of Qubit objects into it's matrix representation
+    """Coverts an Add/Mul of Qubit objects into it's matrix representation
 
     This function is the inverse of ``matrix_to_qubit`` and is a shorthand
     for ``represent(qubit)``.
diff --git a/sympy/solvers/diophantine.py b/sympy/solvers/diophantine.py
index e5378624a8..e7e5479381 100644
--- a/sympy/solvers/diophantine.py
+++ b/sympy/solvers/diophantine.py
@@ -715,6 +715,11 @@ def _diop_quadratic(var, coeff, t):
             # In this case equation can be transformed into a Pell equation
             #n = symbols("n", integer=True)
 
+            fund_solns = solns_pell
+            solns_pell = set(fund_solns)
+            for X, Y in fund_solns:
+                solns_pell.add((-X, -Y))
+
             a = diop_DN(D, 1)
             T = a[0][0]
             U = a[0][1]
@@ -736,50 +741,33 @@ def _diop_quadratic(var, coeff, t):
                     l.add((x_n, y_n))
 
             else:
-                L = ilcm(S(P[0]).q, ilcm(S(P[1]).q, ilcm(S(P[2]).q, ilcm(S(P[3]).q, ilcm(S(Q[0]).q, S(Q[1]).q)))))
+                L = ilcm(S(P[0]).q, ilcm(S(P[1]).q, ilcm(S(P[2]).q,
+                         ilcm(S(P[3]).q, ilcm(S(Q[0]).q, S(Q[1]).q)))))
 
-                k = 0
-                done = False
+                k = 1
 
                 T_k = T
                 U_k = U
 
-                while not done:
-                    k = k + 1
-
-                    if (T_k - 1) % L == 0 and U_k % L == 0:
-                        done = True
+                while (T_k - 1) % L != 0 or U_k % L != 0:
                     T_k, U_k = T_k*T + D*U_k*U, T_k*U + U_k*T
+                    k += 1
 
-                for soln in solns_pell:
-
-                    X = soln[0]
-                    Y = soln[1]
-
-                    done = False
+                for X, Y in solns_pell:
 
                     for i in range(k):
-
-                        X_1 = X*T + D*U*Y
-                        Y_1 = X*U + Y*T
-
-                        x = (P*Matrix([X_1, Y_1]) + Q)[0]
-                        y = (P*Matrix([X_1, Y_1]) + Q)[1]
-
-                        if is_solution_quad(var, coeff, x, y):
-                            done = True
-
-
-                            x_n = S( (X_1 + sqrt(D)*Y_1)*(T + sqrt(D)*U)**(t*L) + (X_1 - sqrt(D)*Y_1)*(T - sqrt(D)*U)**(t*L) )/ 2
-                            y_n = S( (X_1 + sqrt(D)*Y_1)*(T + sqrt(D)*U)**(t*L) - (X_1 - sqrt(D)*Y_1)*(T - sqrt(D)*U)**(t*L) )/ (2*sqrt(D))
-
-                            x_n = _mexpand(x_n)
-                            y_n = _mexpand(y_n)
-                            x_n, y_n = (P*Matrix([x_n, y_n]) + Q)[0], (P*Matrix([x_n, y_n]) + Q)[1]
-                            l.add((x_n, y_n))
-
-                        if done:
-                            break
+                        Z = P*Matrix([X, Y]) + Q
+                        x, y = Z[0], Z[1]
+
+                        if isinstance(x, Integer) and isinstance(y, Integer):
+                            Xt = S((X + sqrt(D)*Y)*(T_k + sqrt(D)*U_k)**t +
+                                  (X - sqrt(D)*Y)*(T_k - sqrt(D)*U_k)**t)/ 2
+                            Yt = S((X + sqrt(D)*Y)*(T_k + sqrt(D)*U_k)**t -
+                                  (X - sqrt(D)*Y)*(T_k - sqrt(D)*U_k)**t)/ (2*sqrt(D))
+                            Zt = P*Matrix([Xt, Yt]) + Q
+                            l.add((Zt[0], Zt[1]))
+
+                        X, Y = X*T + D*U*Y, X*U + Y*T
 
 
     return l


Test patch:
diff --git a/bin/coverage_doctest.py b/bin/coverage_doctest.py
index dfe02413e7..02e4c77198 100755
--- a/bin/coverage_doctest.py
+++ b/bin/coverage_doctest.py
@@ -396,7 +396,7 @@ def coverage(module_path, verbose=False, no_color=False, sphinx=True):
     contained. It then goes through each of the classes/functions to get
     the docstring and doctest coverage of the module. """
 
-    # Import the package and find members
+    # Import the package and find membmers
     m = None
     try:
         __import__(module_path)
diff --git a/sympy/core/tests/test_assumptions.py b/sympy/core/tests/test_assumptions.py
index 4e29f4e8ca..3c62a807f0 100644
--- a/sympy/core/tests/test_assumptions.py
+++ b/sympy/core/tests/test_assumptions.py
@@ -1,4 +1,4 @@
-from sympy import I, sqrt, log, exp, sin, asin, factorial
+from sympy import I, sqrt, log, exp, sin, asin
 from sympy.core import Symbol, S, Rational, Integer, Dummy, Wild, Pow
 from sympy.core.facts import InconsistentAssumptions
 from sympy import simplify
@@ -479,7 +479,7 @@ def test_composite():
     assert S(17).is_composite is False
     assert S(4).is_composite is True
     x = Dummy(integer=True, positive=True, prime=False)
-    assert x.is_composite is None # x could be 1
+    assert x.is_composite
     assert (x + 1).is_composite is None
 
 
@@ -900,8 +900,3 @@ def test_issues_8632_8633_8638_8675_8992():
     assert ((p + 2)**3 - S.Half).is_positive
     n = Dummy(negative=True)
     assert (n - 3).is_nonpositive
-
-def test_issue_9115():
-    n = Dummy('n', integer=True, nonnegative=True)
-    assert (factorial(n) >= 1) == True
-    assert (factorial(n) < 1) == False
diff --git a/sympy/functions/elementary/tests/test_exponential.py b/sympy/functions/elementary/tests/test_exponential.py
index a2057365d8..61991ec99d 100644
--- a/sympy/functions/elementary/tests/test_exponential.py
+++ b/sympy/functions/elementary/tests/test_exponential.py
@@ -1,5 +1,5 @@
 from sympy import (
-    symbols, log, ln, Float, nan, oo, zoo, I, pi, E, exp, Symbol,
+    symbols, log, Float, nan, oo, zoo, I, pi, E, exp, Symbol,
     LambertW, sqrt, Rational, expand_log, S, sign, conjugate,
     sin, cos, sinh, cosh, tanh, exp_polar, re, Function, simplify)
 
@@ -131,10 +131,6 @@ def test_exp_leading_term():
     assert exp(1/x).as_leading_term(x) == exp(1/x)
     assert exp(2 + x).as_leading_term(x) == exp(2)
 
-def test_exp_taylor_term():
-    x = symbols('x')
-    assert exp(x).taylor_term(1, x) == x
-    assert exp(x).taylor_term(3, x) == x**3/6
 
 def test_log_values():
     assert log(nan) == nan
@@ -444,7 +440,6 @@ def test_log_product():
     expr = log(Product(-2, (n, 0, 4)))
     assert simplify(expr) == expr
 
-
 def test_issue_8866():
     x = Symbol('x')
     assert simplify(log(x, 10, evaluate=False)) == simplify(log(x, 10))
@@ -457,10 +452,3 @@ def test_issue_8866():
     b2 = log(exp(y), exp(5), evaluate=False)
     assert simplify(log(l1, b1)) == simplify(log(l2, b2))
     assert expand_log(log(l1, b1)) == expand_log(log(l2, b2))
-
-
-def test_issue_9116():
-    n = Symbol('n', positive=True, integer=True)
-
-    assert ln(n).is_nonnegative is True
-    assert log(n).is_nonnegative is True
diff --git a/sympy/functions/elementary/tests/test_trigonometric.py b/sympy/functions/elementary/tests/test_trigonometric.py
index 5bc62b70d9..f94d1c341f 100644
--- a/sympy/functions/elementary/tests/test_trigonometric.py
+++ b/sympy/functions/elementary/tests/test_trigonometric.py
@@ -1390,8 +1390,3 @@ def test_issue_8653():
     assert sin(n).is_irrational is None
     assert cos(n).is_irrational is None
     assert tan(n).is_irrational is None
-
-
-def test_issue_9157():
-    n = Symbol('n', integer=True, positive=True)
-    atan(n - 1).is_nonnegative is True
diff --git a/sympy/integrals/tests/test_integrals.py b/sympy/integrals/tests/test_integrals.py
index 3b59ab1713..b34db5b37d 100644
--- a/sympy/integrals/tests/test_integrals.py
+++ b/sympy/integrals/tests/test_integrals.py
@@ -1087,9 +1087,3 @@ def test_issue_8901():
     assert integrate(sinh(1.0*x)) == 1.0*cosh(1.0*x)
     assert integrate(tanh(1.0*x)) == 1.0*x - 1.0*log(tanh(1.0*x) + 1)
     assert integrate(tanh(x)) == x - log(tanh(x) + 1)
-
-
-def test_issue_7130():
-    i, L, a, b = symbols('i L a b')
-    integrand = (cos(pi*i*x/L)**2 / (a + b*x)).rewrite(exp)
-    assert x not in integrate(integrand, (x, 0, L)).free_symbols
diff --git a/sympy/matrices/expressions/tests/test_matrix_exprs.py b/sympy/matrices/expressions/tests/test_matrix_exprs.py
index d1f0dcbb23..77a8f1f9ed 100644
--- a/sympy/matrices/expressions/tests/test_matrix_exprs.py
+++ b/sympy/matrices/expressions/tests/test_matrix_exprs.py
@@ -207,9 +207,3 @@ def test_single_indexing():
 
 def test_MatrixElement_diff():
     assert (A[3, 0]*A[0, 0]).diff(A[0, 0]) == A[3, 0]
-
-
-def test_MatrixElement_doit():
-    u = MatrixSymbol('u', 2, 1)
-    v = ImmutableMatrix([3, 5])
-    assert u[0, 0].subs(u, v).doit() == v[0, 0]
diff --git a/sympy/series/tests/test_limits.py b/sympy/series/tests/test_limits.py
index bd160ac9fa..16dd4efeae 100644
--- a/sympy/series/tests/test_limits.py
+++ b/sympy/series/tests/test_limits.py
@@ -37,6 +37,9 @@ def test_basic1():
     assert limit(gamma(1/x + 3), x, oo) == 2
     assert limit(S.NaN, x, -oo) == S.NaN
     assert limit(Order(2)*x, x, S.NaN) == S.NaN
+    assert limit(gamma(1/x + 3), x, oo) == 2
+    assert limit(S.NaN, x, -oo) == S.NaN
+    assert limit(Order(2)*x, x, S.NaN) == S.NaN
     assert limit(1/(x - 1), x, 1, dir="+") == oo
     assert limit(1/(x - 1), x, 1, dir="-") == -oo
     assert limit(1/(5 - x)**3, x, 5, dir="+") == -oo
diff --git a/sympy/series/tests/test_residues.py b/sympy/series/tests/test_residues.py
index 23db6f4fee..676f603677 100644
--- a/sympy/series/tests/test_residues.py
+++ b/sympy/series/tests/test_residues.py
@@ -26,8 +26,9 @@ def test_basic2():
 
 
 def _test_f():
+    # FIXME: we get infinite recursion here:
     f = Function("f")
-    assert residue(f(x)/x**5, x, 0) == f(x).diff(x, 4).subs(x, 0)/24
+    assert residue(f(x)/x**5, x, 0) == f.diff(x, 4)/24
 
 
 def test_functions():
diff --git a/sympy/solvers/tests/test_diophantine.py b/sympy/solvers/tests/test_diophantine.py
index 056124908e..8aea7c93fd 100644
--- a/sympy/solvers/tests/test_diophantine.py
+++ b/sympy/solvers/tests/test_diophantine.py
@@ -118,12 +118,15 @@ def test_quadratic_non_perfect_square():
     assert check_solutions(x**2 - x*y - y**2 - 3*y)
     assert check_solutions(x**2 - 9*y**2 - 2*x - 6*y)
 
+def test_issue_9106():
+    assert check_integrality(-48 - 2*x*(3*x - 1) + y*(3*y - 1))
 
 @slow
 def test_quadratic_non_perfect_slow():
 
     assert check_solutions(8*x**2 + 10*x*y - 2*y**2 - 32*x - 13*y - 23)
-    assert check_solutions(5*x**2 - 13*x*y + y**2 - 4*x - 4*y - 15)
+    # This leads to very large numbers.
+    # assert check_solutions(5*x**2 - 13*x*y + y**2 - 4*x - 4*y - 15)
     assert check_solutions(-3*x**2 - 2*x*y + 7*y**2 - 5*x - 7)
     assert check_solutions(-4 - x + 4*x**2 - y - 3*x*y - 4*y**2)
     assert check_solutions(1 + 2*x + 2*x**2 + 2*y + x*y - 2*y**2)
@@ -579,3 +582,23 @@ def check_solutions(eq):
                 break
 
     return okay
+
+def check_integrality(eq):
+    """
+    Check that the solutions returned by diophantine() are integers.
+    This should be seldom needed except for general quadratic
+    equations which are solved with rational transformations.
+    """
+    def _check_values(x):
+        """ Check a number of values. """
+        for i in range(-4, 4):
+            if not isinstance(simplify(x.subs(t, i)), Integer):
+                return False
+        return True
+
+    for soln in diophantine(eq, param=t):
+        for x in soln:
+            if not _check_values(x):
+                return False
+
+    return True