Repository: sympy/sympy

Issue:
diophantine misses trivial solution
```py
In [6]: diophantine((3*(x**2 + y**2 + z**2) - 14*(x*y + y*z + z*x)))
Out[6]: set([])
```

See http://math.stackexchange.com/questions/1133888/hard-number-theory-problem. The only solution is `(0, 0, 0)`. Nothing in the docstring indicates that it should skip trivial solutions. 

Discussion:
asmeurer: Don't know what your aims are here, but I have found a [couple of issues](https://github.com/sympy/sympy/issues?utf8=%E2%9C%93&q=is%3Aopen+is%3Aissue+author%3Aasmeurer+diophantine) I've found with diophantine in the past that don't appear to be addressed here. 
smichr: > I have found a couple of issues

Those (and others) have now been addressed.

Patch:
diff --git a/doc/src/modules/solvers/diophantine.rst b/doc/src/modules/solvers/diophantine.rst
index 743a28f1d4..cb07f93588 100644
--- a/doc/src/modules/solvers/diophantine.rst
+++ b/doc/src/modules/solvers/diophantine.rst
@@ -45,8 +45,10 @@ structured in the following manner.
     - :py:meth:`~sympy.solvers.diophantine.diop_linear`
     - :py:meth:`~sympy.solvers.diophantine.diop_quadratic`
     - :py:meth:`~sympy.solvers.diophantine.diop_ternary_quadratic`
+    - :py:meth:`~sympy.solvers.diophantine.diop_ternary_quadratic_normal`
     - :py:meth:`~sympy.solvers.diophantine.diop_general_pythagorean`
     - :py:meth:`~sympy.solvers.diophantine.diop_general_sum_of_squares`
+    - :py:meth:`~sympy.solvers.diophantine.diop_general_sum_of_even_powers`
 
   - :py:meth:`~sympy.solvers.diophantine.merge_solution`
 
@@ -135,7 +137,7 @@ For linear Diophantine equations, the customized parameter is the prefix used
 for each free variable in the solution. Consider the following example:
 
 >>> diop_solve(2*x + 3*y - 5*z + 7, m)
-(m_0, -9*m_0 - 5*m_1 - 14, -5*m_0 - 3*m_1 - 7)
+(m_0, m_0 + 5*m_1 - 14, m_0 + 3*m_1 - 7)
 
 In the solution above, m_0 and m_1 are independent free variables.
 
@@ -228,7 +230,7 @@ We can solve an equation of the form `X^2 - DY^2 = N` by passing `D` and `N` to
 >>> diop_DN(5, 920)
 []
 
-Unfortunately, our equation does not have solutions.
+Unfortunately, our equation has no solution.
 
 Now let's turn to homogeneous ternary quadratic equations. These equations are
 of the form `ax^2 + by^2 + cz^2 + dxy + eyz + fzx = 0`. These type of equations
@@ -236,11 +238,11 @@ either have infinitely many solutions or no solutions (except the obvious
 solution (0, 0, 0))
 
 >>> diophantine(3*x**2 + 4*y**2 - 5*z**2 + 4*x*y + 6*y*z + 7*z*x)
-set()
+set([(0, 0, 0)])
 >>> diophantine(3*x**2 + 4*y**2 - 5*z**2 + 4*x*y - 7*y*z + 7*z*x)
 set([(-16*p**2 + 28*p*q + 20*q**2, 3*p**2 + 38*p*q - 25*q**2, 4*p**2 - 24*p*q + 68*q**2)])
 
-If you are only interested about a base solution rather than the parameterized
+If you are only interested in a base solution rather than the parameterized
 general solution (to be more precise, one of the general solutions), you can
 use :py:meth:`~sympy.solvers.diophantine.diop_ternary_quadratic`.
 
@@ -272,15 +274,41 @@ also be solved using the Diophantine module.
 set([(70*t1**2 + 70*t2**2 + 70*t3**2 + 70*t4**2 - 70*t5**2, 105*t1*t5, 420*t2*t5, 60*t3*t5, 210*t4*t5, 42*t1**2 + 42*t2**2 + 42*t3**2 + 42*t4**2 + 42*t5**2)])
 
 function :py:meth:`~sympy.solvers.diophantine.diop_general_pythagorean` can
-also be called directly to solve the same equation. This is true about the
-general sum of squares too. Either you can call
+also be called directly to solve the same equation. Either you can call
 :py:meth:`~sympy.solvers.diophantine.diop_general_pythagorean` or use the high
-level API.
-
->>> diophantine(a**2 + b**2 + c**2 + d**2 + e**2 + f**2 - 112)
-set([(8, 4, 4, 4, 0, 0)])
-
-If you want to get a more thorough idea about the the Diophantine module please
+level API. For the general sum of squares, this is also true, but one advantage
+of calling :py:meth:`~sympy.solvers.diophantine.diop_general_sum_of_squares` is that
+you can control how many solutions are returned.
+
+>>> from sympy.solvers.diophantine import diop_general_sum_of_squares
+>>> eq = a**2 + b**2 + c**2 + d**2 - 18
+>>> diophantine(eq)
+set([(0, 0, 3, 3), (0, 1, 1, 4), (1, 2, 2, 3)])
+>>> diop_general_sum_of_squares(eq, 2)
+set([(0, 0, 3, 3), (1, 2, 2, 3)])
+
+The :py:meth:`~sympy.solvers.diophantine.sum_of_squares` routine will
+providean iterator that returns solutions and one may control whether
+the solutions contain zeros or not (and the solutions not containing
+zeros are returned first):
+
+>>> from sympy.solvers.diophantine import sum_of_squares
+>>> sos = sum_of_squares(18, 4, zeros=True)
+>>> next(sos)
+(1, 2, 2, 3)
+>>> next(sos)
+(0, 0, 3, 3)
+
+
+Simple Eqyptian fractions can be found with the Diophantine module, too.
+For example, here are the ways that one might represent 1/2 as a sum of two
+unit fractions:
+
+>>> from sympy import Eq, S
+>>> diophantine(Eq(1/x + 1/y, S(1)/2))
+set([(-2, 1), (1, -2), (3, 6), (4, 4), (6, 3)])
+
+To get a more thorough understanding of the Diophantine module, please
 refer to the following blog.
 
 http://thilinaatsympy.wordpress.com/
@@ -301,21 +329,28 @@ References
 User Functions
 --------------
 
-These are functions that are imported into the global namespace with ``from
-sympy import *``. These functions are intended for use by ordinary users of SymPy.
+This functions is imported into the global namespace
+with ``from sympy import *``:
 
 :func:`diophantine`
 ^^^^^^^^^^^^^^^^^^^
 .. autofunction:: sympy.solvers.diophantine.diophantine
 
-:func:`diop_solve`
-^^^^^^^^^^^^^^^^^^
-.. autofunction:: sympy.solvers.diophantine.diop_solve
+And this function is imported with ``from sympy.solvers.diophantine import *``:
 
 :func:`classify_diop`
 ^^^^^^^^^^^^^^^^^^^^^
 .. autofunction:: sympy.solvers.diophantine.classify_diop
 
+Internal Functions
+------------------
+
+These functions are intended for internal use in the Diophantine module.
+
+:func:`diop_solve`
+^^^^^^^^^^^^^^^^^^
+.. autofunction:: sympy.solvers.diophantine.diop_solve
+
 :func:`diop_linear`
 ^^^^^^^^^^^^^^^^^^^
 .. autofunction:: sympy.solvers.diophantine.diop_linear
@@ -368,6 +403,10 @@ sympy import *``. These functions are intended for use by ordinary users of SymP
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 .. autofunction:: sympy.solvers.diophantine.diop_general_sum_of_squares
 
+:func:`diop_general_sum_of_even_powers`
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+.. autofunction:: sympy.solvers.diophantine.diop_general_sum_of_even_powers
+
 :func:`partition`
 ^^^^^^^^^^^^^^^^^
 .. autofunction:: sympy.solvers.diophantine.partition
@@ -380,9 +419,13 @@ sympy import *``. These functions are intended for use by ordinary users of SymP
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^
 .. autofunction:: sympy.solvers.diophantine.sum_of_four_squares
 
-Internal Functions
-------------------
-These functions are intended for the internal use in Diophantine module.
+:func:`sum_of_powers`
+^^^^^^^^^^^^^^^^^^^^^
+.. autofunction:: sympy.solvers.diophantine.sum_of_powers
+
+:func:`sum_of_squares`
+^^^^^^^^^^^^^^^^^^^^^^
+.. autofunction:: sympy.solvers.diophantine.sum_of_squares
 
 :obj:`merge_solution`
 ^^^^^^^^^^^^^^^^^^^^^
@@ -392,10 +435,6 @@ These functions are intended for the internal use in Diophantine module.
 ^^^^^^^^^^^^^^^^
 .. autofunction:: sympy.solvers.diophantine.divisible
 
-:obj:`extended_euclid`
-^^^^^^^^^^^^^^^^^^^^^^
-.. autofunction:: sympy.solvers.diophantine.extended_euclid
-
 :obj:`PQa`
 ^^^^^^^^^^
 .. autofunction:: sympy.solvers.diophantine.PQa
@@ -404,10 +443,6 @@ These functions are intended for the internal use in Diophantine module.
 ^^^^^^^^^^^^^^^^^
 .. autofunction:: sympy.solvers.diophantine.equivalent
 
-:obj:`simplified`
-^^^^^^^^^^^^^^^^^
-.. autofunction:: sympy.solvers.diophantine.simplified
-
 :obj:`parametrize_ternary_quadratic`
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 .. autofunction:: sympy.solvers.diophantine.parametrize_ternary_quadratic
@@ -429,16 +464,12 @@ These functions are intended for the internal use in Diophantine module.
 .. autofunction:: sympy.solvers.diophantine.holzer
 
 :obj:`prime_as_sum_of_two_squares`
-^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 .. autofunction:: sympy.solvers.diophantine.prime_as_sum_of_two_squares
 
-:obj:`pairwise_prime`
-^^^^^^^^^^^^^^^^^^^^^
-.. autofunction:: sympy.solvers.diophantine.pairwise_prime
-
-:obj:`make_prime`
+:obj:`sqf_normal`
 ^^^^^^^^^^^^^^^^^
-.. autofunction:: sympy.solvers.diophantine.make_prime
+.. autofunction:: sympy.solvers.diophantine.sqf_normal
 
 :obj:`reconstruct`
 ^^^^^^^^^^^^^^^^^^
diff --git a/sympy/core/add.py b/sympy/core/add.py
index 21caa0bc22..b70c51c12f 100644
--- a/sympy/core/add.py
+++ b/sympy/core/add.py
@@ -333,17 +333,13 @@ def as_coeff_add(self, *deps):
             return coeff, notrat + self.args[1:]
         return S.Zero, self.args
 
-    def as_coeff_Add(self):
+    def as_coeff_Add(self, rational=False):
         """Efficiently extract the coefficient of a summation. """
         coeff, args = self.args[0], self.args[1:]
 
-        if coeff.is_Number:
-            if len(args) == 1:
-                return coeff, args[0]
-            else:
-                return coeff, self._new_rawargs(*args)
-        else:
-            return S.Zero, self
+        if coeff.is_Number and not rational or coeff.is_Rational:
+            return coeff, self._new_rawargs(*args)
+        return S.Zero, self
 
     # Note, we intentionally do not implement Add.as_coeff_mul().  Rather, we
     # let Expr.as_coeff_mul() just always return (S.One, self) for an Add.  See
diff --git a/sympy/core/expr.py b/sympy/core/expr.py
index 3c81fb4ed1..1da1255cf6 100644
--- a/sympy/core/expr.py
+++ b/sympy/core/expr.py
@@ -2832,7 +2832,7 @@ def as_coeff_Mul(self, rational=False):
         """Efficiently extract the coefficient of a product. """
         return S.One, self
 
-    def as_coeff_Add(self):
+    def as_coeff_Add(self, rational=False):
         """Efficiently extract the coefficient of a summation. """
         return S.Zero, self
 
diff --git a/sympy/core/numbers.py b/sympy/core/numbers.py
index ca0b2dd8f3..aeff883ef7 100644
--- a/sympy/core/numbers.py
+++ b/sympy/core/numbers.py
@@ -549,9 +549,11 @@ def as_coeff_Mul(self, rational=False):
             return S.One, self
         return (self, S.One) if self else (S.One, self)
 
-    def as_coeff_Add(self):
+    def as_coeff_Add(self, rational=False):
         """Efficiently extract the coefficient of a summation. """
-        return self, S.Zero
+        if not rational:
+            return self, S.Zero
+        return S.Zero, self
 
     def gcd(self, other):
         """Compute GCD of `self` and `other`. """
@@ -1671,6 +1673,14 @@ def as_content_primitive(self, radical=False, clear=True):
             return -self, S.NegativeOne
         return S.One, self
 
+    def as_coeff_Mul(self, rational=False):
+        """Efficiently extract the coefficient of a product. """
+        return self, S.One
+
+    def as_coeff_Add(self, rational=False):
+        """Efficiently extract the coefficient of a summation. """
+        return self, S.Zero
+
 
 # int -> Integer
 _intcache = {}
@@ -2276,6 +2286,10 @@ def __nonzero__(self):
 
     __bool__ = __nonzero__
 
+    def as_coeff_Mul(self, rational=False):  # XXX this routine should be deleted
+        """Efficiently extract the coefficient of a summation. """
+        return S.One, self
+
 
 class One(with_metaclass(Singleton, IntegerConstant)):
     """The number one.
diff --git a/sympy/ntheory/continued_fraction.py b/sympy/ntheory/continued_fraction.py
index dd94a59f10..07f86f2720 100644
--- a/sympy/ntheory/continued_fraction.py
+++ b/sympy/ntheory/continued_fraction.py
@@ -1,6 +1,7 @@
 from sympy.core.numbers import Integer, Rational
 
 
+
 def continued_fraction_periodic(p, q, d=0):
     r"""
     Find the periodic continued fraction expansion of a quadratic irrational.
diff --git a/sympy/ntheory/factor_.py b/sympy/ntheory/factor_.py
index 9242120dfd..2138c24262 100644
--- a/sympy/ntheory/factor_.py
+++ b/sympy/ntheory/factor_.py
@@ -1669,7 +1669,7 @@ def core(n, t=2):
     See Also
     ========
 
-    factorint
+    factorint, sympy.solvers.diophantine.square_factor
     """
 
     n = as_int(n)
diff --git a/sympy/printing/llvmjitcode.py b/sympy/printing/llvmjitcode.py
index 51e0e65c41..f1e092f8ad 100644
--- a/sympy/printing/llvmjitcode.py
+++ b/sympy/printing/llvmjitcode.py
@@ -34,10 +34,6 @@ def __init__(self, module, builder, fn, *args, **kwargs):
         self.builder = builder
         self.fn = fn
         self.ext_fn = {}  # keep track of wrappers to external functions
-        self.tmp_var = {}
-
-    def _add_tmp_var(self, name, value):
-        self.tmp_var[name] = value
 
     def _print_Number(self, n, **kwargs):
         return ll.Constant(self.fp_type, float(n))
@@ -46,13 +42,8 @@ def _print_Integer(self, expr):
         return ll.Constant(self.fp_type, float(expr.p))
 
     def _print_Symbol(self, s):
-        val = self.tmp_var.get(s)
-        if not val:
-            # look up parameter with name s
-            val = self.func_arg_map.get(s)
-        if not val:
-            raise LookupError("Symbol not found: %s" % s)
-        return val
+        # look up parameter with name s
+        return self.func_arg_map.get(s)
 
     def _print_Pow(self, expr):
         base0 = self._print(expr.base)
@@ -122,13 +113,7 @@ def _print_Indexed(self, expr):
         return value
 
     def _print_Symbol(self, s):
-        val = self.tmp_var.get(s)
-        if val:
-            return val
-
         array, idx = self.func_arg_map.get(s, [None, 0])
-        if not array:
-            raise LookupError("Symbol not found: %s" % s)
         array_ptr = self.builder.gep(array, [ll.Constant(ll.IntType(32), idx)])
         fp_array_ptr = self.builder.bitcast(array_ptr,
                                             ll.PointerType(self.fp_type))
@@ -165,8 +150,6 @@ def _from_ctype(self, ctype):
             return ll.PointerType(self.fp_type)
         if ctype == ctypes.c_void_p:
             return ll.PointerType(ll.IntType(32))
-        if ctype == ctypes.py_object:
-            return ll.PointerType(ll.IntType(32))
 
         print("Unhandled ctype = %s" % str(ctype))
 
@@ -201,65 +184,12 @@ def _create_function(self, expr):
         lj = LLVMJitPrinter(self.module, builder, self.fn,
                             func_arg_map=self.param_dict)
 
-        ret = self._convert_expr(lj, expr)
-        lj.builder.ret(self._wrap_return(lj, ret))
+        ret = lj._print(expr)
+        lj.builder.ret(ret)
 
         strmod = str(self.module)
         return strmod
 
-    def _wrap_return(self, lj, vals):
-        # Return a single double if there is one return value,
-        #  else return a tuple of doubles.
-
-        # Don't wrap return value in this case
-        if self.signature.ret_type == ctypes.c_double:
-            return vals[0]
-
-        # Use this instead of a real PyObject*
-        void_ptr = ll.PointerType(ll.IntType(32))
-
-        # Create a wrapped double: PyObject* PyFloat_FromDouble(double v)
-        wrap_type = ll.FunctionType(void_ptr, [self.fp_type])
-        wrap_fn = ll.Function(lj.module, wrap_type, "PyFloat_FromDouble")
-
-        wrapped_vals = [lj.builder.call(wrap_fn, [v]) for v in vals]
-        if len(vals) == 1:
-            final_val = wrapped_vals[0]
-        else:
-            # Create a tuple: PyObject* PyTuple_Pack(Py_ssize_t n, ...)
-
-            # This should be Py_ssize_t
-            tuple_arg_types = [ll.IntType(32)]
-
-            tuple_arg_types.extend([void_ptr]*len(vals))
-            tuple_type = ll.FunctionType(void_ptr, tuple_arg_types)
-            tuple_fn = ll.Function(lj.module, tuple_type, "PyTuple_Pack")
-
-            tuple_args = [ll.Constant(ll.IntType(32), len(wrapped_vals))]
-            tuple_args.extend(wrapped_vals)
-
-            final_val = lj.builder.call(tuple_fn, tuple_args)
-
-        return final_val
-
-    def _convert_expr(self, lj, expr):
-        try:
-            # Match CSE return data structure.
-            if len(expr) == 2:
-                tmp_exprs = expr[0]
-                final_exprs = expr[1]
-                if len(final_exprs) != 1 and self.signature.ret_type == ctypes.c_double:
-                    raise NotImplementedError("Return of multiple expressions not supported for this callback")
-                for name, e in tmp_exprs:
-                    val = lj._print(e)
-                    lj._add_tmp_var(name, val)
-        except TypeError:
-            final_exprs = [expr]
-
-        vals = [lj._print(e) for e in final_exprs]
-
-        return vals
-
     def _compile_function(self, strmod):
         global exe_engines
         llmod = llvm.parse_assembly(strmod)
@@ -308,18 +238,13 @@ def _create_function(self, expr):
         lj = LLVMJitCallbackPrinter(self.module, builder, self.fn,
                                     func_arg_map=self.param_dict)
 
-        ret = self._convert_expr(lj, expr)
+        ret = lj._print(expr)
 
         if self.signature.ret_arg:
-            output_fp_ptr = builder.bitcast(self.fn.args[self.signature.ret_arg],
-                                            ll.PointerType(self.fp_type))
-            for i, val in enumerate(ret):
-                index = ll.Constant(ll.IntType(32), i)
-                output_array_ptr = builder.gep(output_fp_ptr, [index])
-                builder.store(val, output_array_ptr)
+            builder.store(ret, self.fn.args[self.signature.ret_arg])
             builder.ret(ll.Constant(ll.IntType(32), 0))  # return success
         else:
-            lj.builder.ret(self._wrap_return(lj, ret))
+            lj.builder.ret(ret)
 
         strmod = str(self.module)
         return strmod
@@ -370,7 +295,7 @@ def llvm_callable(args, expr, callback_type=None):
         Arguments to the generated function.  Usually the free symbols in
         the expression.  Currently each one is assumed to convert to
         a double precision scalar.
-    expr : Expr, or (Replacements, Expr) as returned from 'cse'
+    expr : Expr
         Expression to compile.
     callback_type : string
         Create function with signature appropriate to use as a callback.
@@ -417,29 +342,12 @@ def llvm_callable(args, expr, callback_type=None):
     The second ('scipy.integrate.test') is only useful for directly calling
     the function using ctypes variables. It will not pass the signature checks
     for scipy.integrate.
-
-    The return value from the cse module can also be compiled.  This
-    can improve the performance of the compiled function.  If multiple
-    expressions are given to cse, the compiled function returns a tuple.
-    The 'cubature' callback handles multiple expressions (set `fdim`
-    to match in the integration call.)
-    >>> import sympy.printing.llvmjitcode as jit
-    >>> from sympy import cse, exp
-    >>> from sympy.abc import x,y
-    >>> e1 = x*x + y*y
-    >>> e2 = 4*(x*x + y*y) + 8.0
-    >>> after_cse = cse([e1,e2])
-    >>> after_cse
-    ([(x0, x**2), (x1, y**2)], [x0 + x1, 4*x0 + 4*x1 + 8.0])
-    >>> j1 = jit.llvm_callable([x,y], after_cse)
-    >>> j1(1.0, 2.0)
-    (5.0, 28.0)
     '''
 
     if not llvmlite:
         raise ImportError("llvmlite is required for llvmjitcode")
 
-    signature = CodeSignature(ctypes.py_object)
+    signature = CodeSignature(ctypes.c_double)
 
     arg_ctypes = []
     if callback_type is None:
@@ -447,7 +355,6 @@ def llvm_callable(args, expr, callback_type=None):
             arg_ctype = ctypes.c_double
             arg_ctypes.append(arg_ctype)
     elif callback_type == 'scipy.integrate' or callback_type == 'scipy.integrate.test':
-        signature.ret_type = ctypes.c_double
         arg_ctypes = [ctypes.c_int, ctypes.POINTER(ctypes.c_double)]
         arg_ctypes_formal = [ctypes.c_int, ctypes.c_double]
         signature.input_arg = 1
diff --git a/sympy/simplify/fu.py b/sympy/simplify/fu.py
index 3eadde52d1..e74b1a0c54 100644
--- a/sympy/simplify/fu.py
+++ b/sympy/simplify/fu.py
@@ -176,11 +176,8 @@
 
 References
 ==========
-
-Fu, Hongguang, Xiuqin Zhong, and Zhenbing Zeng. "Automated and readable
-simplification of trigonometric expressions." Mathematical and computer
-modelling 44.11 (2006): 1169-1177.
-http://rfdz.ph-noe.ac.at/fileadmin/Mathematik_Uploads/ACDCA/DESTIME2006/DES_contribs/Fu/simplification.pdf
+http://rfdz.ph-noe.ac.at/fileadmin/Mathematik_Uploads/ACDCA/
+DESTIME2006/DES_contribs/Fu/simplification.pdf
 
 http://www.sosmath.com/trig/Trig5/trig5/pdf/pdf.html gives a formula sheet.
 
diff --git a/sympy/simplify/simplify.py b/sympy/simplify/simplify.py
index e93fdf8246..a8413332a6 100644
--- a/sympy/simplify/simplify.py
+++ b/sympy/simplify/simplify.py
@@ -1288,14 +1288,16 @@ def clear_coefficients(expr, rhs=S.Zero):
     """
     was = None
     free = expr.free_symbols
-    while was != expr:
+    if expr.is_Rational:
+        return (S.Zero, rhs - expr)
+    while expr and was != expr:
         was = expr
         m, expr = (
             expr.as_content_primitive()
             if free else
-            factor_terms(expr).as_coeff_Mul())
+            factor_terms(expr).as_coeff_Mul(rational=True))
         rhs /= m
-        c, expr = expr.as_coeff_Add()
+        c, expr = expr.as_coeff_Add(rational=True)
         rhs -= c
     expr = signsimp(expr, evaluate = False)
     if _coeff_isneg(expr):
diff --git a/sympy/solvers/__init__.py b/sympy/solvers/__init__.py
index ba51543f56..b5e7f6355f 100644
--- a/sympy/solvers/__init__.py
+++ b/sympy/solvers/__init__.py
@@ -12,6 +12,8 @@
     solve_undetermined_coeffs, nsolve, solve_linear, checksol, \
     det_quick, inv_quick
 
+from .diophantine import diophantine
+
 from .recurr import rsolve, rsolve_poly, rsolve_ratio, rsolve_hyper
 
 from .ode import checkodesol, classify_ode, dsolve, \
diff --git a/sympy/solvers/diophantine.py b/sympy/solvers/diophantine.py
index f85bf48485..856e09deae 100644
--- a/sympy/solvers/diophantine.py
+++ b/sympy/solvers/diophantine.py
@@ -1,47 +1,121 @@
 from __future__ import print_function, division
 
-from sympy import (Add, ceiling, divisors, factor_list, factorint, floor, igcd,
-    ilcm, Integer, integer_nthroot, isprime, Matrix, Mul, nextprime,
-    perfect_power, Poly, S, sign, solve, sqrt, Subs, Symbol, symbols, sympify,
-    Wild)
-
+from sympy.core.add import Add
+from sympy.core.compatibility import as_int, is_sequence, range
+from sympy.core.exprtools import factor_terms
 from sympy.core.function import _mexpand
-from sympy.simplify.radsimp import rad_rationalize
-from sympy.utilities import default_sort_key, numbered_symbols
-from sympy.core.numbers import igcdex
-from sympy.ntheory.residue_ntheory import sqrt_mod
-from sympy.core.compatibility import range
+from sympy.core.mul import Mul
+from sympy.core.numbers import Rational
+from sympy.core.numbers import igcdex, ilcm, igcd
+from sympy.core.power import integer_nthroot
 from sympy.core.relational import Eq
+from sympy.core.singleton import S
+from sympy.core.symbol import Symbol, symbols
+from sympy.functions.elementary.complexes import sign
+from sympy.functions.elementary.integers import floor
+from sympy.functions.elementary.miscellaneous import sqrt
+from sympy.matrices.dense import MutableDenseMatrix as Matrix
+from sympy.ntheory.factor_ import (
+    divisors, factorint, multiplicity, perfect_power)
+from sympy.ntheory.generate import nextprime
+from sympy.ntheory.primetest import is_square, isprime
+from sympy.ntheory.residue_ntheory import sqrt_mod
+from sympy.polys.polyerrors import GeneratorsNeeded
+from sympy.polys.polytools import Poly, factor_list
+from sympy.simplify.simplify import signsimp
 from sympy.solvers.solvers import check_assumptions
+from sympy.solvers.solveset import solveset_real
+from sympy.utilities import default_sort_key, numbered_symbols
+from sympy.utilities.misc import filldedent
+
+
+# these are imported with 'from sympy.solvers.diophantine import *
+__all__ = ['diophantine', 'classify_diop']
+
+
+# these types are known (but not necessarily handled)
+diop_known = {
+    "binary_quadratic",
+    "cubic_thue",
+    "general_pythagorean",
+    "general_sum_of_even_powers",
+    "general_sum_of_squares",
+    "homogeneous_general_quadratic",
+    "homogeneous_ternary_quadratic",
+    "homogeneous_ternary_quadratic_normal",
+    "inhomogeneous_general_quadratic",
+    "inhomogeneous_ternary_quadratic",
+    "linear",
+    "univariate"}
+
+
+def _is_int(i):
+    try:
+        as_int(i)
+        return True
+    except ValueError:
+        pass
+
+
+def _sorted_tuple(*i):
+    return tuple(sorted(i))
+
 
-__all__ = ['base_solution_linear', 'classify_diop', 'cornacchia', 'descent',
-           'diop_bf_DN', 'diop_DN', 'diop_general_pythagorean',
-           'diop_general_sum_of_squares', 'diop_linear', 'diop_quadratic',
-           'diop_solve', 'diop_ternary_quadratic', 'diophantine', 'find_DN',
-           'partition', 'square_factor', 'sum_of_four_squares',
-           'sum_of_three_squares', 'transformation_to_DN']
+def _remove_gcd(*x):
+    try:
+        g = igcd(*x)
+        return tuple([i//g for i in x])
+    except ValueError:
+        return x
+    except TypeError:
+        raise TypeError('_remove_gcd(a,b,c) or _remove_gcd(*container)')
 
-def diophantine(eq, param=symbols("t", integer=True)):
+
+def _rational_pq(a, b):
+    # return `(numer, denom)` for a/b; sign in numer and gcd removed
+    return _remove_gcd(sign(b)*a, abs(b))
+
+
+def _nint_or_floor(p, q):
+    # return nearest int to p/q; in case of tie return floor(p/q)
+    w, r = divmod(p, q)
+    if abs(r) <= abs(q)//2:
+        return w
+    return w + 1
+
+
+def _odd(i):
+    return i % 2 != 0
+
+
+def _even(i):
+    return i % 2 == 0
+
+
+def diophantine(eq, param=symbols("t", integer=True), syms=None):
     """
     Simplify the solution procedure of diophantine equation ``eq`` by
     converting it into a product of terms which should equal zero.
 
     For example, when solving, `x^2 - y^2 = 0` this is treated as
-    `(x + y)(x - y) = 0` and `x+y = 0` and `x-y = 0` are solved independently
-    and combined. Each term is solved by calling ``diop_solve()``.
+    `(x + y)(x - y) = 0` and `x + y = 0` and `x - y = 0` are solved
+    independently and combined. Each term is solved by calling
+    ``diop_solve()``.
 
-    Output of ``diophantine()`` is a set of tuples. Each tuple represents a
-    solution of the input equation. In a tuple, solution for each variable is
-    listed according to the alphabetic order of input variables. i.e. if we have
-    an equation with two variables `a` and `b`, first element of the tuple will
-    give the solution for `a` and the second element will give the solution for
-    `b`.
+    Output of ``diophantine()`` is a set of tuples. The elements of the
+    tuple are the solutions for each variable in the the equation and
+    are arranged according to the alphabetic ordering of the variables.
+    e.g. For an equation with two variables, `a` and `b`, the first
+    element of the tuple is the solution for `a` and the second for `b`.
 
     Usage
     =====
 
-    ``diophantine(eq, t)``: Solve the diophantine equation ``eq``.
-    ``t`` is the parameter to be used by ``diop_solve()``.
+    ``diophantine(eq, t, syms, factor)``: Solve the diophantine
+    equation ``eq``.
+    ``t`` is the optional parameter to be used by ``diop_solve()``.
+    ``syms`` is an optional list of symbols which determines the
+    order of the elements in the returned tuple.
 
     Details
     =======
@@ -55,12 +129,12 @@ def diophantine(eq, param=symbols("t", integer=True)):
     >>> from sympy.solvers.diophantine import diophantine
     >>> from sympy.abc import x, y, z
     >>> diophantine(x**2 - y**2)
-    set([(-t_0, -t_0), (t_0, -t_0)])
+    set([(t_0, -t_0), (t_0, t_0)])
 
-    #>>> diophantine(x*(2*x + 3*y - z))
-    #set([(0, n1, n2), (3*t - z, -2*t + z, z)])
-    #>>> diophantine(x**2 + 3*x*y + 4*x)
-    #set([(0, n1), (3*t - 4, -t)])
+    >>> diophantine(x*(2*x + 3*y - z))
+    set([(0, n1, n2), (t_0, t_1, 2*t_0 + 3*t_1)])
+    >>> diophantine(x**2 + 3*x*y + 4*x)
+    set([(0, n1), (3*t_0 - 4, -t_0)])
 
     See Also
     ========
@@ -70,36 +144,85 @@ def diophantine(eq, param=symbols("t", integer=True)):
     if isinstance(eq, Eq):
         eq = eq.lhs - eq.rhs
 
-    eq = Poly(eq).as_expr()
-    if not eq.is_polynomial() or eq.is_number:
-        raise TypeError("Equation input format not supported")
-
-    var = list(eq.expand(force=True).free_symbols)
-    var.sort(key=default_sort_key)
-
-    terms = factor_list(eq)[1]
+    try:
+        var = list(eq.expand(force=True).free_symbols)
+        var.sort(key=default_sort_key)
+        if syms:
+            if not is_sequence(syms):
+                raise TypeError(
+                    'syms should be given as a sequence, e.g. a list')
+            syms = [i for i in syms if i in var]
+            if syms != var:
+                map = dict(zip(syms, range(len(syms))))
+                return set([tuple([t[map[i]] for i in var])
+                    for t in diophantine(eq, param)])
+        n, d = eq.as_numer_denom()
+        if not n.free_symbols:
+            return set()
+        if d.free_symbols:
+            dsol = diophantine(d)
+            good = diophantine(n) - dsol
+            return set([s for s in good if _mexpand(d.subs(zip(var, s)))])
+        else:
+            eq = n
+        eq = factor_terms(eq)
+        assert not eq.is_number
+        eq = eq.as_independent(*var, as_Add=False)[1]
+        p = Poly(eq)
+        assert not any(g.is_number for g in p.gens)
+        eq = p.as_expr()
+        assert eq.is_polynomial()
+    except (GeneratorsNeeded, AssertionError, AttributeError):
+        raise TypeError(filldedent('''
+    Equation should be a polynomial with Rational coefficients.'''))
+
+    try:
+        # if we know that factoring should not be attempted, skip
+        # the factoring step
+        v, c, t = classify_diop(eq)
+        if t  == 'general_sum_of_squares':
+            # trying to factor such expressions will sometimes hang
+            terms = [(eq, 1)]
+        else:
+            raise TypeError
+    except (TypeError, NotImplementedError):
+        terms = factor_list(eq)[1]
 
     sols = set([])
 
     for term in terms:
 
-        base = term[0]
-
-        var_t, jnk, eq_type = classify_diop(base)
-        if not var_t:
-            continue
+        base, _ = term
+        var_t, _, eq_type = classify_diop(base, _dict=False)
+        _, base = signsimp(base, evaluate=False).as_coeff_Mul()
         solution = diop_solve(base, param)
 
-        if eq_type in ["linear", "homogeneous_ternary_quadratic", "general_pythagorean"]:
-            if merge_solution(var, var_t, solution) != ():
-                sols.add(merge_solution(var, var_t, solution))
-
-        elif eq_type in ["binary_quadratic",  "general_sum_of_squares", "univariate"]:
+        if eq_type in [
+                "linear",
+                "homogeneous_ternary_quadratic",
+                "homogeneous_ternary_quadratic_normal",
+                "general_pythagorean"]:
+            sols.add(merge_solution(var, var_t, solution))
+
+        elif eq_type in [
+                "binary_quadratic",
+                "general_sum_of_squares",
+                "general_sum_of_even_powers",
+                "univariate"]:
             for sol in solution:
-                if merge_solution(var, var_t, sol) != ():
-                    sols.add(merge_solution(var, var_t, sol))
+                sols.add(merge_solution(var, var_t, sol))
 
-    return sols
+        else:
+            raise NotImplementedError('unhandled type: %s' % eq_type)
+
+    # remove null merge results
+    if () in sols:
+        sols.remove(())
+    null = tuple([0]*len(var))
+    # if there is no solution, return trivial solution
+    if not sols and eq.subs(zip(var, null)) is S.Zero:
+        sols.add(null)
+    return set([S(i) for i in sols])
 
 
 def merge_solution(var, var_t, solution):
@@ -108,39 +231,39 @@ def merge_solution(var, var_t, solution):
     equations.
 
     For example when solving the equation `(x - y)(x^2 + y^2 - z^2) = 0`,
-    solutions for each of the equations `x-y = 0` and `x^2 + y^2 - z^2` are
+    solutions for each of the equations `x - y = 0` and `x^2 + y^2 - z^2` are
     found independently. Solutions for `x - y = 0` are `(x, y) = (t, t)`. But
     we should introduce a value for z when we output the solution for the
     original equation. This function converts `(t, t)` into `(t, t, n_{1})`
     where `n_{1}` is an integer parameter.
     """
-    l = []
+    sol = []
 
     if None in solution:
         return ()
 
     solution = iter(solution)
-    params = numbered_symbols("n", Integer=True, start=1)
+    params = numbered_symbols("n", integer=True, start=1)
     for v in var:
         if v in var_t:
-            l.append(next(solution))
+            sol.append(next(solution))
         else:
-            l.append(next(params))
+            sol.append(next(params))
 
-    for val, symb in zip(l, var):
+    for val, symb in zip(sol, var):
         if check_assumptions(val, **symb.assumptions0) is False:
             return tuple()
 
-    return tuple(l)
+    return tuple(sol)
 
 
 def diop_solve(eq, param=symbols("t", integer=True)):
     """
     Solves the diophantine equation ``eq``.
 
-    Similar to ``diophantine()`` but doesn't try to factor ``eq`` as latter
-    does. Uses ``classify_diop()`` to determine the type of the eqaution and
-    calls the appropriate solver function.
+    Unlike ``diophantine()``, factoring of ``eq`` is not attempted. Uses
+    ``classify_diop()`` to determine the type of the equation and calls
+    the appropriate solver function.
 
     Usage
     =====
@@ -161,19 +284,19 @@ def diop_solve(eq, param=symbols("t", integer=True)):
     >>> from sympy.abc import x, y, z, w
     >>> diop_solve(2*x + 3*y - 5)
     (3*t_0 - 5, -2*t_0 + 5)
-    >>> diop_solve(4*x + 3*y -4*z + 5)
-    (t_0, -4*t_1 + 5, t_0 - 3*t_1 + 5)
-    >>> diop_solve(x + 3*y - 4*z + w -6)
-    (t_0, t_0 + t_1, -2*t_0 - 3*t_1 - 4*t_2 - 6, -t_0 - 2*t_1 - 3*t_2 - 6)
+    >>> diop_solve(4*x + 3*y - 4*z + 5)
+    (t_0, 8*t_0 + 4*t_1 + 5, 7*t_0 + 3*t_1 + 5)
+    >>> diop_solve(x + 3*y - 4*z + w - 6)
+    (t_0, t_0 + t_1, 6*t_0 + 5*t_1 + 4*t_2 - 6, 5*t_0 + 4*t_1 + 3*t_2 - 6)
     >>> diop_solve(x**2 + y**2 - 5)
-    set([(-2, -1), (-2, 1), (2, -1), (2, 1)])
+    set([(-1, -2), (-1, 2), (1, -2), (1, 2)])
 
     See Also
     ========
 
     diophantine()
     """
-    var, coeff, eq_type = classify_diop(eq)
+    var, coeff, eq_type = classify_diop(eq, _dict=False)
 
     if eq_type == "linear":
         return _diop_linear(var, coeff, param)
@@ -183,167 +306,158 @@ def diop_solve(eq, param=symbols("t", integer=True)):
 
     elif eq_type == "homogeneous_ternary_quadratic":
         x_0, y_0, z_0 = _diop_ternary_quadratic(var, coeff)
-        return _parametrize_ternary_quadratic((x_0, y_0, z_0), var, coeff)
+        return _parametrize_ternary_quadratic(
+            (x_0, y_0, z_0), var, coeff)
+
+    elif eq_type == "homogeneous_ternary_quadratic_normal":
+        x_0, y_0, z_0 = _diop_ternary_quadratic_normal(var, coeff)
+        return _parametrize_ternary_quadratic(
+            (x_0, y_0, z_0), var, coeff)
 
     elif eq_type == "general_pythagorean":
         return _diop_general_pythagorean(var, coeff, param)
 
     elif eq_type == "univariate":
-        l = solve(eq)
-        s = set([])
-
-        for soln in l:
-            if isinstance(soln, Integer):
-                s.add((soln,))
-        return s
+        return set([(int(i),) for i in solveset_real(
+            eq, var[0]).intersect(S.Integers)])
 
     elif eq_type == "general_sum_of_squares":
-        return _diop_general_sum_of_squares(var, coeff)
-
-
-def classify_diop(eq):
-    """
-    Helper routine used by diop_solve() to find the type of the ``eq`` etc.
-
-    Returns a tuple containing the type of the diophantine equation along with
-    the variables(free symbols) and their coefficients. Variables are returned
-    as a list and coefficients are returned as a dict with the key being the
-    respective term and the constant term is keyed to Integer(1). Type is an
-    element in the set {"linear", "binary_quadratic", "general_pythagorean",
-    "homogeneous_ternary_quadratic", "univariate", "general_sum_of_squares"}
-
-    Usage
-    =====
-
-    ``classify_diop(eq)``: Return variables, coefficients and type of the
-    ``eq``.
-
-    Details
-    =======
-
-    ``eq`` should be an expression which is assumed to be zero.
+        return _diop_general_sum_of_squares(var, -int(coeff[1]), limit=S.Infinity)
+
+    elif eq_type == "general_sum_of_even_powers":
+        for k in coeff.keys():
+            if k.is_Pow and coeff[k]:
+                p = k.exp
+        return _diop_general_sum_of_even_powers(var, p, -int(coeff[1]), limit=S.Infinity)
+
+    if eq_type is not None and eq_type not in diop_known:
+            raise ValueError(filldedent('''
+    Alhough this type of equation was identified, it is not yet
+    handled. It should, however, be listed in `diop_known` at the
+    top of this file. Developers should see comments at the end of
+    `classify_diop`.
+            '''))  # pragma: no cover
+    else:
+        raise NotImplementedError(
+            'No solver has been written for %s.' % eq_type)
 
-    Examples
-    ========
 
-    >>> from sympy.solvers.diophantine import classify_diop
-    >>> from sympy.abc import x, y, z, w, t
-    >>> classify_diop(4*x + 6*y - 4)
-    ([x, y], {1: -4, x: 4, y: 6}, 'linear')
-    >>> classify_diop(x + 3*y -4*z + 5)
-    ([x, y, z], {1: 5, x: 1, y: 3, z: -4}, 'linear')
-    >>> classify_diop(x**2 + y**2 - x*y + x + 5)
-    ([x, y], {1: 5, x: 1, x**2: 1, y: 0, y**2: 1, x*y: -1}, 'binary_quadratic')
-    """
+def classify_diop(eq, _dict=True):
+    # docstring supplied externally
+    try:
+        var = list(eq.free_symbols)
+        assert var
+    except (AttributeError, AssertionError):
+        raise ValueError('equation should have 1 or more free symbols')
+    var.sort(key=default_sort_key)
     eq = eq.expand(force=True)
-
     coeff = eq.as_coefficients_dict()
-    diop_type = None
-
-    var = []
-    if isinstance(eq, Symbol):
-        var.append(eq)
-        coeff[eq] = Integer(1)
-    elif isinstance(eq, Mul) and Poly(eq).total_degree() == 1:
-        var.append(eq.as_two_terms()[1])
-        coeff[eq.as_two_terms()[1]] = Integer(eq.as_two_terms()[0])
-    else:
-        var = list(eq.free_symbols)
-        var.sort(key=default_sort_key)
-        coeff = dict([reversed(t.as_independent(*var)) for t in eq.args])
-
-    for c in coeff:
-        if not isinstance(coeff[c], Integer):
-            raise TypeError("Coefficients should be Integers")
+    if not all(_is_int(c) for c in coeff.values()):
+        raise TypeError("Coefficients should be Integers")
 
-    if Poly(eq).total_degree() == 1:
+    diop_type = None
+    total_degree = Poly(eq).total_degree()
+    homogeneous = 1 not in coeff
+    if total_degree == 1:
         diop_type = "linear"
+
     elif len(var) == 1:
         diop_type = "univariate"
-    elif Poly(eq).total_degree() == 2 and len(var) == 2:
+
+    elif total_degree == 2 and len(var) == 2:
         diop_type = "binary_quadratic"
-        x, y = var[:2]
 
-        if isinstance(eq, Mul):
-            coeff = {x**2: 0, x*y: eq.args[0], y**2: 0, x: 0, y: 0, Integer(1): 0}
+    elif total_degree == 2 and len(var) == 3 and homogeneous:
+        if set(coeff) & set(var):
+            diop_type = "inhomogeneous_ternary_quadratic"
         else:
-            for term in [x**2, y**2, x*y, x, y, Integer(1)]:
-                if term not in coeff.keys():
-                    coeff[term] = Integer(0)
-
-    elif Poly(eq).total_degree() == 2 and len(var) == 3 and Integer(1) not in coeff.keys():
-        for v in var:
-            if v in coeff.keys():
-                diop_type = "inhomogeneous_ternary_quadratic"
-                break
-        else:
-            diop_type = "homogeneous_ternary_quadratic"
+            nonzero = [k for k in coeff if coeff[k]]
+            if len(nonzero) == 3 and all(i**2 in nonzero for i in var):
+                diop_type = "homogeneous_ternary_quadratic_normal"
+            else:
+                diop_type = "homogeneous_ternary_quadratic"
 
-            x, y, z = var[:3]
+    elif total_degree == 2 and len(var) >= 3:
+        if set(coeff) & set(var):
+            diop_type = "inhomogeneous_general_quadratic"
+        else:
+            # there may be Pow keys like x**2 or Mul keys like x*y
+            if any(k.is_Mul for k in coeff):  # cross terms
+                if not homogeneous:
+                    diop_type = "inhomogeneous_general_quadratic"
+                else:
+                    diop_type = "homogeneous_general_quadratic"
+            else:  # all squares: x**2 + y**2 + ... + constant
+                if all(coeff[k] == 1 for k in coeff if k != 1):
+                    diop_type = "general_sum_of_squares"
+                elif all(is_square(abs(coeff[k])) for k in coeff):
+                    if abs(sum(sign(coeff[k]) for k in coeff)) == \
+                            len(var) - 2:
+                        # all but one has the same sign
+                        # e.g. 4*x**2 + y**2 - 4*z**2
+                        diop_type = "general_pythagorean"
 
-            for term in [x**2, y**2, z**2, x*y, y*z, x*z]:
-                if term not in coeff.keys():
-                    coeff[term] = Integer(0)
+    elif total_degree == 3 and len(var) == 2:
+        diop_type = "cubic_thue"
 
-    elif Poly(eq).degree() == 2 and len(var) >= 3:
+    elif (total_degree > 3 and total_degree % 2 == 0 and
+            all(k.is_Pow for k in coeff if k != 1)):
+        if all(coeff[k] == 1 for k in coeff if k != 1):
+            diop_type = 'general_sum_of_even_powers'
 
-        for v in var:
-            if v in coeff.keys():
-                diop_type = "inhomogeneous_general_quadratic"
-                break
+    if diop_type is not None:
+        return var, dict(coeff) if _dict else coeff, diop_type
 
-        else:
-            if Integer(1) in coeff.keys():
-                constant_term = True
-            else:
-                constant_term = False
+    # new diop type instructions
+    # --------------------------
+    # if this error raises and the equation *can* be classified,
+    #  * it should be identified in the if-block above
+    #  * the type should be added to the diop_known
+    # if a solver can be written for it,
+    #  * a dedicated handler should be written (e.g. diop_linear)
+    #  * it should be passed to that handler in diop_solve
+    raise NotImplementedError(filldedent('''
+        This equation is not yet recognized or else has not been
+        simplified sufficiently to put it in a form recognized by
+        diop_classify().'''))
 
-            non_square_degree_2_terms = False
-            for v in var:
-                for u in var:
-                    if u != v and u*v in coeff.keys():
-                        non_square_degree_2_terms = True
-                        break
-                if non_square_degree_2_terms:
-                    break
 
-            if constant_term and non_square_degree_2_terms:
-                diop_type = "inhomogeneous_general_quadratic"
+classify_diop.func_doc = '''
+    Helper routine used by diop_solve() to find information about ``eq``.
 
-            elif constant_term and not non_square_degree_2_terms:
-                for v in var:
-                    if coeff[v**2] != 1:
-                        break
-                else:
-                    diop_type = "general_sum_of_squares"
+    Returns a tuple containing the type of the diophantine equation
+    along with the variables (free symbols) and their coefficients.
+    Variables are returned as a list and coefficients are returned
+    as a dict with the key being the respective term and the constant
+    term is keyed to 1. The type is one of the following:
 
-            elif not constant_term and non_square_degree_2_terms:
-                diop_type = "homogeneous_general_quadratic"
+    * %s
 
-            else:
-                coeff_sign_sum = 0
+    Usage
+    =====
 
-                for v in var:
-                    if not isinstance(sqrt(abs(Integer(coeff[v**2]))), Integer):
-                        break
-                    coeff_sign_sum = coeff_sign_sum + sign(coeff[v**2])
-                else:
-                    if abs(coeff_sign_sum) == len(var) - 2 and not constant_term:
-                        diop_type = "general_pythagorean"
+    ``classify_diop(eq)``: Return variables, coefficients and type of the
+    ``eq``.
 
-    elif Poly(eq).total_degree() == 3 and len(var) == 2:
+    Details
+    =======
 
-        x, y = var[:2]
-        diop_type = "cubic_thue"
+    ``eq`` should be an expression which is assumed to be zero.
+    ``_dict`` is for internal use: when True (default) a dict is returned,
+    otherwise a defaultdict which supplies 0 for missing keys is returned.
 
-        for term in [x**3, x**2*y, x*y**2, y**3, Integer(1)]:
-            if term not in coeff.keys():
-                coeff[term] == Integer(0)
+    Examples
+    ========
 
-    if diop_type is not None:
-        return var, coeff, diop_type
-    else:
-        raise NotImplementedError("Still not implemented")
+    >>> from sympy.solvers.diophantine import classify_diop
+    >>> from sympy.abc import x, y, z, w, t
+    >>> classify_diop(4*x + 6*y - 4)
+    ([x, y], {1: -4, x: 4, y: 6}, 'linear')
+    >>> classify_diop(x + 3*y -4*z + 5)
+    ([x, y, z], {1: 5, x: 1, y: 3, z: -4}, 'linear')
+    >>> classify_diop(x**2 + y**2 - x*y + x + 5)
+    ([x, y], {1: 5, x: 1, x**2: 1, y**2: 1, x*y: -1}, 'binary_quadratic')
+    ''' % ('\n    * '.join(sorted(diop_known)))
 
 
 def diop_linear(eq, param=symbols("t", integer=True)):
@@ -372,14 +486,13 @@ def diop_linear(eq, param=symbols("t", integer=True)):
 
     >>> from sympy.solvers.diophantine import diop_linear
     >>> from sympy.abc import x, y, z, t
-    >>> from sympy import Integer
-    >>> diop_linear(2*x - 3*y - 5) #solves equation 2*x - 3*y -5 = 0
-    (-3*t_0 - 5, -2*t_0 - 5)
+    >>> diop_linear(2*x - 3*y - 5) # solves equation 2*x - 3*y - 5 == 0
+    (3*t_0 - 5, 2*t_0 - 5)
 
     Here x = -3*t_0 - 5 and y = -2*t_0 - 5
 
     >>> diop_linear(2*x - 3*y - 4*z -3)
-    (t_0, -6*t_0 - 4*t_1 + 3, 5*t_0 + 3*t_1 - 3)
+    (t_0, 2*t_0 + 4*t_1 + 3, -t_0 - 3*t_1 - 3)
 
     See Also
     ========
@@ -387,7 +500,8 @@ def diop_linear(eq, param=symbols("t", integer=True)):
     diop_quadratic(), diop_ternary_quadratic(), diop_general_pythagorean(),
     diop_general_sum_of_squares()
     """
-    var, coeff, diop_type = classify_diop(eq)
+    from sympy.core.function import count_ops
+    var, coeff, diop_type = classify_diop(eq, _dict=False)
 
     if diop_type == "linear":
         return _diop_linear(var, coeff, param)
@@ -401,32 +515,29 @@ def _diop_linear(var, coeff, param):
 
     Note that no solution exists if gcd(a_0, ..., a_n) doesn't divide c.
     """
-    if len(var) == 0:
-        return None
 
-    if Integer(1) in coeff:
-        #coeff[] is negated because input is of the form: ax + by + c ==  0
-        #                                 but is used as: ax + by     == -c
-        c = -coeff[Integer(1)]
+    if 1 in coeff:
+        # negate coeff[] because input is of the form: ax + by + c ==  0
+        #                              but is used as: ax + by     == -c
+        c = -coeff[1]
     else:
         c = 0
 
     # Some solutions will have multiple free variables in their solutions.
-    params = [str(param) + "_" + str(i) for i in range(len(var))]
-    params = [symbols(p, integer=True) for p in params]
+    if param is None:
+        params = [symbols('t')]*len(var)
+    else:
+        temp = str(param) + "_%i"
+        params = [symbols(temp % i, integer=True) for i in range(len(var))]
 
     if len(var) == 1:
-        if coeff[var[0]] == 0:
-            if c == 0:
-                return tuple([params[0]])
-            else:
-                return tuple([None])
-        elif divisible(c, coeff[var[0]]):
-            return tuple([c/coeff[var[0]]])
+        q, r = divmod(c, coeff[var[0]])
+        if not r:
+            return (q,)
         else:
-            return tuple([None])
+            return (None,)
 
-    """
+    '''
     base_solution_linear() can solve diophantine equations of the form:
 
     a*x + b*y == c
@@ -466,7 +577,7 @@ def _diop_linear(var, coeff, param):
 
     The arrays A and B are the arrays of integers used for
     'a' and 'b' in each of the n-1 bivariate equations we solve.
-    """
+    '''
 
     A = [coeff[v] for v in var]
     B = []
@@ -481,7 +592,7 @@ def _diop_linear(var, coeff, param):
             B.insert(0, gcd)
     B.append(A[-1])
 
-    """
+    '''
     Consider the trivariate linear equation:
 
     4*x_0 + 6*x_1 + 3*x_2 == 2
@@ -532,64 +643,51 @@ def _diop_linear(var, coeff, param):
 
     This method is generalised for many variables, below.
 
-    """
+    '''
 
     solutions = []
-    no_solution = tuple([None] * len(var))
     for i in range(len(B)):
-        tot_x, tot_y = 0, 0
+        tot_x, tot_y = [], []
 
-        if isinstance(c, Add):
-            # example: 5 + t_0 + 3*t_1
-            args = c.args
-        else: # c is a Mul, a Symbol, or an Integer
-            args = [c]
-
-        for j in range(len(args)):
-            if isinstance(args[j], Mul):
+        for j, arg in enumerate(Add.make_args(c)):
+            if arg.is_Integer:
+                # example: 5 -> k = 5
+                k, p = arg, S.One
+                pnew = params[0]
+            else:  # arg is a Mul or Symbol
                 # example: 3*t_1 -> k = 3
-                k = args[j].as_two_terms()[0]
-                param_index = params.index(args[j].as_two_terms()[1]) + 1
-            elif isinstance(args[j], Symbol):
                 # example: t_0 -> k = 1
-                k = 1
-                param_index = params.index(args[j]) + 1
-            else: #args[j] is an Integer
-                # example: 5 -> k = 5
-                k = args[j]
-                param_index = 0
+                k, p = arg.as_coeff_Mul()
+                pnew = params[params.index(p) + 1]
 
-            sol_x, sol_y = base_solution_linear(k, A[i], B[i], params[param_index])
-            if isinstance(args[j], Mul) or isinstance(args[j], Symbol):
-                if isinstance(sol_x, Add):
-                    sol_x = sol_x.args[0]*params[param_index - 1] + sol_x.args[1]
-                elif isinstance(sol_x, Integer):
-                    sol_x = sol_x*params[param_index - 1]
-
-                if isinstance(sol_y, Add):
-                    sol_y = sol_y.args[0]*params[param_index - 1] + sol_y.args[1]
-                elif isinstance(sol_y, Integer):
-                    sol_y = sol_y*params[param_index - 1]
+            sol = sol_x, sol_y = base_solution_linear(k, A[i], B[i], pnew)
 
+            if p is S.One:
+                if None in sol:
+                    return tuple([None]*len(var))
             else:
-                if sol_x is None or sol_y is None:
-                    return no_solution
-
-            tot_x += sol_x
-            tot_y += sol_y
+                # convert a + b*pnew -> a*p + b*pnew
+                if isinstance(sol_x, Add):
+                    sol_x = sol_x.args[0]*p + sol_x.args[1]
+                if isinstance(sol_y, Add):
+                    sol_y = sol_y.args[0]*p + sol_y.args[1]
 
-        solutions.append(tot_x)
-        c = tot_y
+            tot_x.append(sol_x)
+            tot_y.append(sol_y)
 
-    solutions.append(tot_y)
+        solutions.append(Add(*tot_x))
+        c = Add(*tot_y)
 
+    solutions.append(c)
+    if param is None:
+        # just keep the additive constant (i.e. replace t with 0)
+        solutions = [i.as_coeff_Add()[0] for i in solutions]
     return tuple(solutions)
 
 
 def base_solution_linear(c, a, b, t=None):
     """
-    Return the base solution for a linear diophantine equation with two
-    variables.
+    Return the base solution for the linear equation, `ax + by = c`.
 
     Used by ``diop_linear()`` to find the base solution of a linear
     Diophantine equation. If ``t`` is given then the parametrized solution is
@@ -615,70 +713,37 @@ def base_solution_linear(c, a, b, t=None):
     >>> base_solution_linear(0, 5, 7, t) # equation 5*x + 7*y = 0
     (7*t, -5*t)
     """
-    d = igcd(a, igcd(b, c))
-    a = a // d
-    b = b // d
-    c = c // d
+    a, b, c = _remove_gcd(a, b, c)
 
     if c == 0:
-        if t != None:
+        if t is not None:
+            if b < 0:
+                t = -t
             return (b*t , -a*t)
         else:
-            return (S.Zero, S.Zero)
+            return (0, 0)
     else:
-        x0, y0, d = extended_euclid(int(abs(a)), int(abs(b)))
+        x0, y0, d = igcdex(abs(a), abs(b))
 
-        x0 = x0 * sign(a)
-        y0 = y0 * sign(b)
+        x0 *= sign(a)
+        y0 *= sign(b)
 
         if divisible(c, d):
-            if t != None:
+            if t is not None:
+                if b < 0:
+                    t = -t
                 return (c*x0 + b*t, c*y0 - a*t)
             else:
-                return (Integer(c*x0), Integer(c*y0))
+                return (c*x0, c*y0)
         else:
             return (None, None)
 
 
-def extended_euclid(a, b):
-    """
-    For given ``a``, ``b`` returns a tuple containing integers `x`, `y` and `d`
-    such that `ax + by = d`. Here `d = gcd(a, b)`.
-
-    Usage
-    =====
-
-    ``extended_euclid(a, b)``: returns `x`, `y` and `\gcd(a, b)`.
-
-    Details
-    =======
-
-    ``a`` Any instance of Integer.
-    ``b`` Any instance of Integer.
-
-    Examples
-    ========
-
-    >>> from sympy.solvers.diophantine import extended_euclid
-    >>> extended_euclid(4, 6)
-    (-1, 1, 2)
-    >>> extended_euclid(3, 5)
-    (2, -1, 1)
-    """
-    if b == 0:
-        return (1, 0, a)
-
-    x0, y0, d = extended_euclid(b, a%b)
-    x, y = y0, x0 - (a//b) * y0
-
-    return x, y, d
-
-
 def divisible(a, b):
     """
     Returns `True` if ``a`` is divisible by ``b`` and `False` otherwise.
     """
-    return igcd(int(a), int(b)) == abs(int(b))
+    return not a % b
 
 
 def diop_quadratic(eq, param=symbols("t", integer=True)):
@@ -713,7 +778,7 @@ def diop_quadratic(eq, param=symbols("t", integer=True)):
     References
     ==========
 
-    .. [1] Methods to solve Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0,[online],
+    .. [1] Methods to solve Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0, [online],
           Available: http://www.alpertron.com.ar/METHODS.HTM
     .. [2] Solving the equation ax^2+ bxy + cy^2 + dx + ey + f= 0, [online],
           Available: http://www.jpr2718.org/ax2p.pdf
@@ -724,7 +789,7 @@ def diop_quadratic(eq, param=symbols("t", integer=True)):
     diop_linear(), diop_ternary_quadratic(), diop_general_sum_of_squares(),
     diop_general_pythagorean()
     """
-    var, coeff, diop_type = classify_diop(eq)
+    var, coeff, diop_type = classify_diop(eq, _dict=False)
 
     if diop_type == "binary_quadratic":
         return _diop_quadratic(var, coeff, param)
@@ -732,134 +797,121 @@ def diop_quadratic(eq, param=symbols("t", integer=True)):
 
 def _diop_quadratic(var, coeff, t):
 
-    x, y = var[:2]
-
-    for term in [x**2, y**2, x*y, x, y, Integer(1)]:
-        if term not in coeff.keys():
-            coeff[term] = Integer(0)
+    x, y = var
 
     A = coeff[x**2]
     B = coeff[x*y]
     C = coeff[y**2]
     D = coeff[x]
     E = coeff[y]
-    F = coeff[Integer(1)]
+    F = coeff[1]
 
-    d = igcd(A, igcd(B, igcd(C, igcd(D, igcd(E, F)))))
-    A = A // d
-    B = B // d
-    C = C // d
-    D = D // d
-    E = E // d
-    F = F // d
+    A, B, C, D, E, F = [as_int(i) for i in _remove_gcd(A, B, C, D, E, F)]
 
-    # (1) Linear case: A = B = C = 0 ==> considered under linear diophantine equations
-
-    # (2) Simple-Hyperbolic case:A = C = 0, B != 0
+    # (1) Simple-Hyperbolic case: A = C = 0, B != 0
     # In this case equation can be converted to (Bx + E)(By + D) = DE - BF
     # We consider two cases; DE - BF = 0 and DE - BF != 0
     # More details, http://www.alpertron.com.ar/METHODS.HTM#SHyperb
 
-    l = set([])
-
+    sol = set([])
+    discr = B**2 - 4*A*C
     if A == 0 and C == 0 and B != 0:
 
         if D*E - B*F == 0:
-            if divisible(int(E), int(B)):
-                l.add((-E/B, t))
-            if divisible(int(D), int(B)):
-                l.add((t, -D/B))
-
+            q, r = divmod(E, B)
+            if not r:
+                sol.add((-q, t))
+            q, r = divmod(D, B)
+            if not r:
+                sol.add((t, -q))
         else:
             div = divisors(D*E - B*F)
             div = div + [-term for term in div]
-
             for d in div:
-                if divisible(int(d - E), int(B)):
-                    x0  = (d - E) // B
-                    if divisible(int(D*E - B*F), int(d)):
-                        if divisible(int((D*E - B*F)// d - D), int(B)):
-                            y0 = ((D*E - B*F) // d - D) // B
-                            l.add((x0, y0))
-
-    # (3) Parabolic case: B**2 - 4*A*C = 0
+                x0, r = divmod(d - E, B)
+                if not r:
+                    q, r = divmod(D*E - B*F, d)
+                    if not r:
+                        y0, r = divmod(q - D, B)
+                        if not r:
+                            sol.add((x0, y0))
+
+    # (2) Parabolic case: B**2 - 4*A*C = 0
     # There are two subcases to be considered in this case.
     # sqrt(c)D - sqrt(a)E = 0 and sqrt(c)D - sqrt(a)E != 0
     # More Details, http://www.alpertron.com.ar/METHODS.HTM#Parabol
 
-    elif B**2 - 4*A*C == 0:
+    elif discr == 0:
 
         if A == 0:
             s = _diop_quadratic([y, x], coeff, t)
             for soln in s:
-                l.add((soln[1], soln[0]))
+                sol.add((soln[1], soln[0]))
 
         else:
-            g = igcd(A, C)
-            g = abs(g) * sign(A)
+            g = sign(A)*igcd(A, C)
             a = A // g
             b = B // g
             c = C // g
             e = sign(B/A)
 
 
-            if e*sqrt(c)*D - sqrt(a)*E == 0:
+            _c = e*sqrt(c)*D - sqrt(a)*E
+            if _c == 0:
                 z = symbols("z", real=True)
-                roots = solve(sqrt(a)*g*z**2 + D*z + sqrt(a)*F)
+                eq = sqrt(a)*g*z**2 + D*z + sqrt(a)*F
+                roots = solveset_real(eq, z).intersect(S.Integers)
                 for root in roots:
-                    if isinstance(root, Integer):
-                        l.add((diop_solve(sqrt(a)*x + e*sqrt(c)*y - root)[0], diop_solve(sqrt(a)*x + e*sqrt(c)*y - root)[1]))
+                    ans = diop_solve(sqrt(a)*x + e*sqrt(c)*y - root)
+                    sol.add((ans[0], ans[1]))
 
-            elif isinstance(e*sqrt(c)*D - sqrt(a)*E, Integer):
-                solve_x = lambda u: e*sqrt(c)*g*(sqrt(a)*E - e*sqrt(c)*D)*t**2 - (E + 2*e*sqrt(c)*g*u)*t\
-                    - (e*sqrt(c)*g*u**2 + E*u + e*sqrt(c)*F) // (e*sqrt(c)*D - sqrt(a)*E)
+            elif _is_int(_c):
+                solve_x = lambda u: -e*sqrt(c)*g*_c*t**2 - (E + 2*e*sqrt(c)*g*u)*t\
+                    - (e*sqrt(c)*g*u**2 + E*u + e*sqrt(c)*F) // _c
 
-                solve_y = lambda u: sqrt(a)*g*(e*sqrt(c)*D - sqrt(a)*E)*t**2 + (D + 2*sqrt(a)*g*u)*t \
-                    + (sqrt(a)*g*u**2 + D*u + sqrt(a)*F) // (e*sqrt(c)*D - sqrt(a)*E)
+                solve_y = lambda u: sqrt(a)*g*_c*t**2 + (D + 2*sqrt(a)*g*u)*t \
+                    + (sqrt(a)*g*u**2 + D*u + sqrt(a)*F) // _c
 
-                for z0 in range(0, abs(e*sqrt(c)*D - sqrt(a)*E)):
-                    if divisible(sqrt(a)*g*z0**2 + D*z0 + sqrt(a)*F, e*sqrt(c)*D - sqrt(a)*E):
-                        l.add((solve_x(z0), solve_y(z0)))
+                for z0 in range(0, abs(_c)):
+                    if divisible(
+                            sqrt(a)*g*z0**2 + D*z0 + sqrt(a)*F,
+                            _c):
+                        sol.add((solve_x(z0), solve_y(z0)))
 
-    # (4) Method used when B**2 - 4*A*C is a square, is descibed in p. 6 of the below paper
+    # (3) Method used when B**2 - 4*A*C is a square, is described in p. 6 of the below paper
     # by John P. Robertson.
     # http://www.jpr2718.org/ax2p.pdf
 
-    elif isinstance(sqrt(B**2 - 4*A*C), Integer):
+    elif is_square(discr):
         if A != 0:
-            r = sqrt(B**2 - 4*A*C)
+            r = sqrt(discr)
             u, v = symbols("u, v", integer=True)
-            eq = _mexpand(4*A*r*u*v + 4*A*D*(B*v + r*u + r*v - B*u) + 2*A*4*A*E*(u - v) + 4*A*r*4*A*F)
+            eq = _mexpand(
+                4*A*r*u*v + 4*A*D*(B*v + r*u + r*v - B*u) +
+                2*A*4*A*E*(u - v) + 4*A*r*4*A*F)
 
-            sol = diop_solve(eq, t)
-            sol = list(sol)
+            solution = diop_solve(eq, t)
 
-            for solution in sol:
-                s0 = solution[0]
-                t0 = solution[1]
+            for s0, t0 in solution:
 
-                x_0 = S(B*t0 + r*s0 + r*t0 - B*s0)/(4*A*r)
+                num = B*t0 + r*s0 + r*t0 - B*s0
+                x_0 = S(num)/(4*A*r)
                 y_0 = S(s0 - t0)/(2*r)
-
                 if isinstance(s0, Symbol) or isinstance(t0, Symbol):
                     if check_param(x_0, y_0, 4*A*r, t) != (None, None):
-                        l.add((check_param(x_0, y_0, 4*A*r, t)[0], check_param(x_0, y_0, 4*A*r, t)[1]))
+                        ans = check_param(x_0, y_0, 4*A*r, t)
+                        sol.add((ans[0], ans[1]))
+                elif x_0.is_Integer and y_0.is_Integer:
+                    if is_solution_quad(var, coeff, x_0, y_0):
+                        sol.add((x_0, y_0))
 
-                elif divisible(B*t0 + r*s0 + r*t0 - B*s0, 4*A*r):
-                    if divisible(s0 - t0, 2*r):
-                        if is_solution_quad(var, coeff, x_0, y_0):
-                            l.add((x_0, y_0))
         else:
-            _var = var
-            _var[0], _var[1] = _var[1], _var[0] # Interchange x and y
-            s = _diop_quadratic(_var, coeff, t)
-
-            while len(s) > 0:
-                sol = s.pop()
-                l.add((sol[1], sol[0]))
+            s = _diop_quadratic(var[::-1], coeff, t)  # Interchange x and y
+            while s:                                  #         |
+                sol.add(s.pop()[::-1])  # and solution <--------+
 
 
-    # (5) B**2 - 4*A*C > 0 and B**2 - 4*A*C not a square or B**2 - 4*A*C < 0
+    # (4) B**2 - 4*A*C > 0 and B**2 - 4*A*C not a square or B**2 - 4*A*C < 0
 
     else:
 
@@ -871,42 +923,34 @@ def _diop_quadratic(var, coeff, t):
             for solution in solns_pell:
                 for X_i in [-solution[0], solution[0]]:
                     for Y_i in [-solution[1], solution[1]]:
-                        x_i, y_i = (P*Matrix([X_i, Y_i]) + Q)[0], (P*Matrix([X_i, Y_i]) + Q)[1]
-                        if isinstance(x_i, Integer) and isinstance(y_i, Integer):
-                            l.add((x_i, y_i))
+                        s = P*Matrix([X_i, Y_i]) + Q
+                        try:
+                            sol.add(tuple([as_int(_) for _ in s]))
+                        except ValueError:
+                            pass
 
         else:
             # 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 = set(solns_pell)
+            for X, Y in list(solns_pell):
                 solns_pell.add((-X, -Y))
 
             a = diop_DN(D, 1)
             T = a[0][0]
             U = a[0][1]
 
-            if (isinstance(P[0], Integer) and isinstance(P[1], Integer) and isinstance(P[2], Integer)
-                and isinstance(P[3], Integer) and isinstance(Q[0], Integer) and isinstance(Q[1], Integer)):
-
-                for sol in solns_pell:
-
-                    r = sol[0]
-                    s = sol[1]
-                    x_n = S((r + s*sqrt(D))*(T + U*sqrt(D))**t + (r - s*sqrt(D))*(T - U*sqrt(D))**t)/2
-                    y_n = S((r + s*sqrt(D))*(T + U*sqrt(D))**t - (r - s*sqrt(D))*(T - U*sqrt(D))**t)/(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 all(_is_int(_) for _ in P[:4] + Q[:2]):
+                for r, s in solns_pell:
+                    _a = (r + s*sqrt(D))*(T + U*sqrt(D))**t
+                    _b = (r - s*sqrt(D))*(T - U*sqrt(D))**t
+                    x_n = _mexpand(S(_a + _b)/2)
+                    y_n = _mexpand(S(_a - _b)/(2*sqrt(D)))
+                    s = P*Matrix([x_n, y_n]) + Q
+                    sol.add(tuple(s))
 
             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(*[_.q for _ in P[:4] + Q[:2]])
 
                 k = 1
 
@@ -920,21 +964,18 @@ def _diop_quadratic(var, coeff, t):
                 for X, Y in solns_pell:
 
                     for i in range(k):
-                        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]))
+                        if all(_is_int(_) for _ in P*Matrix([X, Y]) + Q):
+                            _a = (X + sqrt(D)*Y)*(T_k + sqrt(D)*U_k)**t
+                            _b = (X - sqrt(D)*Y)*(T_k - sqrt(D)*U_k)**t
+                            Xt = S(_a + _b)/2
+                            Yt = S(_a - _b)/(2*sqrt(D))
+                            s = P*Matrix([Xt, Yt]) + Q
+                            sol.add(tuple(s))
 
                         X, Y = X*T + D*U*Y, X*U + Y*T
 
 
-    return l
+    return sol
 
 
 def is_solution_quad(var, coeff, u, v):
@@ -945,24 +986,22 @@ def is_solution_quad(var, coeff, u, v):
 
     Not intended for use by normal users.
     """
-    x, y = var[:2]
-
-    eq = x**2*coeff[x**2] + x*y*coeff[x*y] + y**2*coeff[y**2] + x*coeff[x] + y*coeff[y] + coeff[Integer(1)]
-
-    return _mexpand(Subs(eq, (x, y), (u, v)).doit()) == 0
+    reps = dict(zip(var, (u, v)))
+    eq = Add(*[j*i.xreplace(reps) for i, j in coeff.items()])
+    return _mexpand(eq) == 0
 
 
 def diop_DN(D, N, t=symbols("t", integer=True)):
     """
     Solves the equation `x^2 - Dy^2 = N`.
 
-    Mainly concerned in the case `D > 0, D` is not a perfect square, which is
-    the same as generalized Pell equation. To solve the generalized Pell
-    equation this function Uses LMM algorithm. Refer [1]_ for more details on
-    the algorithm.
-    Returns one solution for each class of the solutions. Other solutions of
-    the class can be constructed according to the values of ``D`` and ``N``.
-    Returns a list containing the solution tuples `(x, y)`.
+    Mainly concerned with the case `D > 0, D` is not a perfect square,
+    which is the same as the generalized Pell equation. The LMM
+    algorithm [1]_ is used to solve this equation.
+
+    Returns one solution tuple, (`x, y)` for each class of the solutions.
+    Other solutions of the class can be constructed according to the
+    values of ``D`` and ``N``.
 
     Usage
     =====
@@ -1005,52 +1044,55 @@ def diop_DN(D, N, t=symbols("t", integer=True)):
     """
     if D < 0:
         if N == 0:
-            return [(S.Zero, S.Zero)]
+            return [(0, 0)]
         elif N < 0:
             return []
         elif N > 0:
-            d = divisors(square_factor(N))
             sol = []
-
-            for divisor in d:
-                sols = cornacchia(1, -D, N // divisor**2)
+            for d in divisors(square_factor(N)):
+                sols = cornacchia(1, -D, N // d**2)
                 if sols:
                     for x, y in sols:
-                        sol.append((divisor*x, divisor*y))
+                        sol.append((d*x, d*y))
 
             return sol
 
     elif D == 0:
-        if N < 0 or not isinstance(sqrt(N), Integer):
+        if N < 0:
             return []
         if N == 0:
-            return [(S.Zero, t)]
-        if isinstance(sqrt(N), Integer):
-            return [(sqrt(N), t)]
-
-    else: # D > 0
-        if isinstance(sqrt(D), Integer):
-            r = sqrt(D)
+            return [(0, t)]
+        sN, _exact = integer_nthroot(N, 2)
+        if _exact:
+            return [(sN, t)]
+        else:
+            return []
 
+    else:  # D > 0
+        sD, _exact = integer_nthroot(D, 2)
+        if _exact:
             if N == 0:
-                return [(r*t, t)]
+                return [(sD*t, t)]
             else:
                 sol = []
 
-                for y in range(floor(sign(N)*(N - 1)/(2*r)) + 1):
-                    if isinstance(sqrt(D*y**2 + N), Integer):
-                        sol.append((sqrt(D*y**2 + N), y))
+                for y in range(floor(sign(N)*(N - 1)/(2*sD)) + 1):
+                    try:
+                        sq, _exact = integer_nthroot(D*y**2 + N, 2)
+                    except ValueError:
+                        _exact = False
+                    if _exact:
+                        sol.append((sq, y))
 
                 return sol
         else:
             if N == 0:
-                return [(S.Zero, S.Zero)]
+                return [(0, 0)]
 
             elif abs(N) == 1:
 
                 pqa = PQa(0, 1, D)
-                a_0 = floor(sqrt(D))
-                l = 0
+                j = 0
                 G = []
                 B = []
 
@@ -1060,29 +1102,29 @@ def diop_DN(D, N, t=symbols("t", integer=True)):
                     G.append(i[5])
                     B.append(i[4])
 
-                    if l != 0 and a == 2*a_0:
+                    if j != 0 and a == 2*sD:
                         break
-                    l = l + 1
+                    j = j + 1
 
-                if l % 2 == 1:
+                if _odd(j):
 
                     if N == -1:
-                        x = G[l-1]
-                        y = B[l-1]
+                        x = G[j - 1]
+                        y = B[j - 1]
                     else:
-                        count = l
-                        while count < 2*l - 1:
+                        count = j
+                        while count < 2*j - 1:
                             i = next(pqa)
                             G.append(i[5])
                             B.append(i[4])
-                            count = count + 1
+                            count += 1
 
                         x = G[count]
                         y = B[count]
                 else:
                     if N == 1:
-                        x = G[l-1]
-                        y = B[l-1]
+                        x = G[j - 1]
+                        y = B[j - 1]
                     else:
                         return []
 
@@ -1100,18 +1142,17 @@ def diop_DN(D, N, t=symbols("t", integer=True)):
 
                 for f in fs:
                     m = N // f**2
-                    zs = sqrt_mod(D, abs(m), True)
 
+                    zs = sqrt_mod(D, abs(m), all_roots=True)
                     zs = [i for i in zs if i <= abs(m) // 2 ]
+
                     if abs(m) != 2:
-                        zs = zs + [-i for i in zs]
-                        if S.Zero in zs:
-                            zs.remove(S.Zero) # Remove duplicate zero
+                        zs = zs + [-i for i in zs if i]  # omit dupl 0
 
                     for z in zs:
 
                         pqa = PQa(z, abs(m), D)
-                        l = 0
+                        j = 0
                         G = []
                         B = []
 
@@ -1121,9 +1162,9 @@ def diop_DN(D, N, t=symbols("t", integer=True)):
                             G.append(i[5])
                             B.append(i[4])
 
-                            if l != 0 and abs(i[1]) == 1:
-                                r = G[l-1]
-                                s = B[l-1]
+                            if j != 0 and abs(i[1]) == 1:
+                                r = G[j-1]
+                                s = B[j-1]
 
                                 if r**2 - D*s**2 == m:
                                     sol.append((f*r, f*s))
@@ -1134,8 +1175,8 @@ def diop_DN(D, N, t=symbols("t", integer=True)):
 
                                 break
 
-                            l = l + 1
-                            if l == length(z, abs(m), D):
+                            j = j + 1
+                            if j == length(z, abs(m), D):
                                 break
 
                 return sol
@@ -1148,8 +1189,8 @@ def cornacchia(a, b, m):
     Uses the algorithm due to Cornacchia. The method only finds primitive
     solutions, i.e. ones with `\gcd(x, y) = 1`. So this method can't be used to
     find the solutions of `x^2 + y^2 = 20` since the only solution to former is
-    `(x,y) = (4, 2)` and it is not primitive. When ` a = b = 1`, only the
-    solutions with `x \geq y` are found. For more details, see the References.
+    `(x, y) = (4, 2)` and it is not primitive. When `a = b`, only the
+    solutions with `x \leq y` are found. For more details, see the References.
 
     Examples
     ========
@@ -1158,7 +1199,7 @@ def cornacchia(a, b, m):
     >>> cornacchia(2, 3, 35) # equation 2x**2 + 3y**2 = 35
     set([(2, 3), (4, 1)])
     >>> cornacchia(1, 1, 25) # equation x**2 + y**2 = 25
-    set([(4, 3)])
+    set([(3, 4)])
 
     References
     ===========
@@ -1167,18 +1208,18 @@ def cornacchia(a, b, m):
     .. [2] Solving the diophantine equation ax**2 + by**2 = m by Cornacchia's
         method, [online], Available:
         http://www.numbertheory.org/php/cornacchia.html
+
+    See Also
+    ========
+    sympy.utilities.iterables.signed_permutations
     """
-    sols = set([])
+    sols = set()
 
     a1 = igcdex(a, m)[0]
-    v = sqrt_mod(-b*a1, m, True)
-
-    if v is None:
+    v = sqrt_mod(-b*a1, m, all_roots=True)
+    if not v:
         return None
 
-    if not isinstance(v, list):
-        v = [v]
-
     for t in v:
         if t < m // 2:
             continue
@@ -1194,8 +1235,10 @@ def cornacchia(a, b, m):
 
         if m1 % b == 0:
             m1 = m1 // b
-            if isinstance(sqrt(m1), Integer):
-                s = sqrt(m1)
+            s, _exact = integer_nthroot(m1, 2)
+            if _exact:
+                if a == b and r > s:
+                    r, s = s, r
                 sols.add((int(r), int(s)))
 
     return sols
@@ -1303,6 +1346,9 @@ def diop_bf_DN(D, N, t=symbols("t", integer=True)):
     .. [1] Solving the generalized Pell equation x**2 - D*y**2 = N, John P.
         Robertson, July 31, 2004, Page 15. http://www.jpr2718.org/pell.pdf
     """
+    D = as_int(D)
+    N = as_int(N)
+
     sol = []
     a = diop_DN(D, 1)
     u = a[0][0]
@@ -1314,27 +1360,33 @@ def diop_bf_DN(D, N, t=symbols("t", integer=True)):
 
     elif N > 1:
         L1 = 0
-        L2 = floor(sqrt(S(N*(u - 1))/(2*D))) + 1
+        L2 = integer_nthroot(int(N*(u - 1)/(2*D)), 2)[0] + 1
 
     elif N < -1:
-        L1 = ceiling(sqrt(S(-N)/D))
-        L2 = floor(sqrt(S(-N*(u + 1))/(2*D))) + 1
+        L1, _exact = integer_nthroot(-int(N/D), 2)
+        if not _exact:
+            L1 += 1
+        L2 = integer_nthroot(-int(N*(u + 1)/(2*D)), 2)[0] + 1
 
-    else:
+    else:  # N = 0
         if D < 0:
-            return [(S.Zero, S.Zero)]
+            return [(0, 0)]
         elif D == 0:
-            return [(S.Zero, t)]
+            return [(0, t)]
         else:
-            if isinstance(sqrt(D), Integer):
-                return [(sqrt(D)*t, t), (-sqrt(D)*t, t)]
+            sD, _exact = integer_nthroot(D, 2)
+            if _exact:
+                return [(sD*t, t), (-sD*t, t)]
             else:
-                return [(S.Zero, S.Zero)]
+                return [(0, 0)]
 
 
     for y in range(L1, L2):
-        if isinstance(sqrt(N + D*y**2), Integer):
-            x = sqrt(N + D*y**2)
+        try:
+            x, _exact = integer_nthroot(N + D*y**2, 2)
+        except ValueError:
+            _exact = False
+        if _exact:
             sol.append((x, y))
             if not equivalent(x, y, -x, y, D, N):
                 sol.append((-x, y))
@@ -1384,7 +1436,7 @@ def length(P, Q, D):
     continued fraction representation of `\\frac{P + \sqrt{D}}{Q}`.
 
     It is important to remember that this does NOT return the length of the
-    periodic part but the addition of the legths of the two parts as mentioned
+    periodic part but the sum of the legths of the two parts as mentioned
     above.
 
     Usage
@@ -1406,35 +1458,21 @@ def length(P, Q, D):
     >>> length(-2 , 4, 5) # (-2 + sqrt(5))/4
     3
     >>> length(-5, 4, 17) # (-5 + sqrt(17))/4
-    4
-    """
-    x = P + sqrt(D)
-    y = Q
-
-    x = sympify(x)
-    v, res = [], []
-    q = x/y
-
-    if q < 0:
-        v.append(q)
-        res.append(floor(q))
-        q = q - floor(q)
-        num, den = rad_rationalize(1, q)
-        q = num / den
-
-    while 1:
-        v.append(q)
-        a = int(q)
-        res.append(a)
-
-        if q == a:
-            return len(res)
-
-        num, den = rad_rationalize(1,(q - a))
-        q = num / den
+    5
 
-        if q in v:
-            return len(res)
+    See Also
+    ========
+    sympy.ntheory.continued_fraction.continued_fraction_periodic
+    """
+    from sympy.ntheory.continued_fraction import continued_fraction_periodic
+    v = continued_fraction_periodic(P, Q, D)
+    if type(v[-1]) is list:
+        rpt = len(v[-1])
+        nonrpt = len(v) - 1
+    else:
+        rpt = 0
+        nonrpt = len(v)
+    return rpt + nonrpt
 
 
 def transformation_to_DN(eq):
@@ -1477,7 +1515,7 @@ def transformation_to_DN(eq):
     will give an equation of the form `x^2 - Dy^2 = N`.
 
     >>> from sympy.abc import X, Y
-    >>> from sympy import Matrix, simplify, Subs
+    >>> from sympy import Matrix, simplify
     >>> u = (A*Matrix([X, Y]) + B)[0] # Transformation for x
     >>> u
     X/26 + 3*Y/26 - 6/13
@@ -1488,7 +1526,7 @@ def transformation_to_DN(eq):
     Next we will substitute these formulas for `x` and `y` and do
     ``simplify()``.
 
-    >>> eq = simplify(Subs(x**2 - 3*x*y - y**2 - 2*y + 1, (x, y), (u, v)).doit())
+    >>> eq = simplify((x**2 - 3*x*y - y**2 - 2*y + 1).subs(zip((x, y), (u, v))))
     >>> eq
     X**2/676 - Y**2/52 + 17/13
 
@@ -1513,65 +1551,52 @@ def transformation_to_DN(eq):
            http://www.jpr2718.org/ax2p.pdf
     """
 
-
-    var, coeff, diop_type = classify_diop(eq)
+    var, coeff, diop_type = classify_diop(eq, _dict=False)
     if diop_type == "binary_quadratic":
         return _transformation_to_DN(var, coeff)
 
 
 def _transformation_to_DN(var, coeff):
 
-    x, y = var[:2]
+    x, y = var
 
     a = coeff[x**2]
     b = coeff[x*y]
     c = coeff[y**2]
     d = coeff[x]
     e = coeff[y]
-    f = coeff[Integer(1)]
+    f = coeff[1]
 
-    g = igcd(a, igcd(b, igcd(c, igcd(d, igcd(e, f)))))
-    a = a // g
-    b = b // g
-    c = c // g
-    d = d // g
-    e = e // g
-    f = f // g
+    a, b, c, d, e, f = [as_int(i) for i in _remove_gcd(a, b, c, d, e, f)]
 
     X, Y = symbols("X, Y", integer=True)
 
-    if b != Integer(0):
-        B = (S(2*a)/b).p
-        C = (S(2*a)/b).q
-        A = (S(a)/B**2).p
-        T = (S(a)/B**2).q
+    if b:
+        B, C = _rational_pq(2*a, b)
+        A, T = _rational_pq(a, B**2)
 
         # eq_1 = A*B*X**2 + B*(c*T - A*C**2)*Y**2 + d*T*X + (B*e*T - d*T*C)*Y + f*T*B
-        coeff = {X**2: A*B, X*Y: 0, Y**2: B*(c*T - A*C**2), X: d*T, Y: B*e*T - d*T*C, Integer(1): f*T*B}
+        coeff = {X**2: A*B, X*Y: 0, Y**2: B*(c*T - A*C**2), X: d*T, Y: B*e*T - d*T*C, 1: f*T*B}
         A_0, B_0 = _transformation_to_DN([X, Y], coeff)
         return Matrix(2, 2, [S(1)/B, -S(C)/B, 0, 1])*A_0, Matrix(2, 2, [S(1)/B, -S(C)/B, 0, 1])*B_0
 
     else:
-        if d != Integer(0):
-            B = (S(2*a)/d).p
-            C = (S(2*a)/d).q
-            A = (S(a)/B**2).p
-            T = (S(a)/B**2).q
+        if d:
+            B, C = _rational_pq(2*a, d)
+            A, T = _rational_pq(a, B**2)
 
             # eq_2 = A*X**2 + c*T*Y**2 + e*T*Y + f*T - A*C**2
-            coeff = {X**2: A, X*Y: 0, Y**2: c*T, X: 0, Y: e*T, Integer(1): f*T - A*C**2}
+            coeff = {X**2: A, X*Y: 0, Y**2: c*T, X: 0, Y: e*T, 1: f*T - A*C**2}
             A_0, B_0 = _transformation_to_DN([X, Y], coeff)
             return Matrix(2, 2, [S(1)/B, 0, 0, 1])*A_0, Matrix(2, 2, [S(1)/B, 0, 0, 1])*B_0 + Matrix([-S(C)/B, 0])
 
         else:
-            if e != Integer(0):
-                B = (S(2*c)/e).p
-                C = (S(2*c)/e).q
-                A = (S(c)/B**2).p
-                T = (S(c)/B**2).q
+            if e:
+                B, C = _rational_pq(2*c, e)
+                A, T = _rational_pq(c, B**2)
 
                 # eq_3 = a*T*X**2 + A*Y**2 + f*T - A*C**2
-                coeff = {X**2: a*T, X*Y: 0, Y**2: A, X: 0, Y: 0, Integer(1): f*T - A*C**2}
+                coeff = {X**2: a*T, X*Y: 0, Y**2: A, X: 0, Y: 0, 1: f*T - A*C**2}
                 A_0, B_0 = _transformation_to_DN([X, Y], coeff)
                 return Matrix(2, 2, [1, 0, 0, S(1)/B])*A_0, Matrix(2, 2, [1, 0, 0, S(1)/B])*B_0 + Matrix([0, -S(C)/B])
 
@@ -1621,81 +1646,54 @@ def find_DN(eq):
            John P.Robertson, May 8, 2003, Page 7 - 11.
            http://www.jpr2718.org/ax2p.pdf
     """
-    var, coeff, diop_type = classify_diop(eq)
+    var, coeff, diop_type = classify_diop(eq, _dict=False)
     if diop_type == "binary_quadratic":
         return _find_DN(var, coeff)
 
 
 def _find_DN(var, coeff):
 
-    x, y = var[:2]
+    x, y = var
     X, Y = symbols("X, Y", integer=True)
-    A , B = _transformation_to_DN(var, coeff)
+    A, B = _transformation_to_DN(var, coeff)
 
     u = (A*Matrix([X, Y]) + B)[0]
     v = (A*Matrix([X, Y]) + B)[1]
-    eq = x**2*coeff[x**2] + x*y*coeff[x*y] + y**2*coeff[y**2] + x*coeff[x] + y*coeff[y] + coeff[Integer(1)]
-
-    simplified = _mexpand(Subs(eq, (x, y), (u, v)).doit())
+    eq = x**2*coeff[x**2] + x*y*coeff[x*y] + y**2*coeff[y**2] + x*coeff[x] + y*coeff[y] + coeff[1]
 
-    coeff = dict([reversed(t.as_independent(*[X, Y])) for t in simplified.args])
+    simplified = _mexpand(eq.subs(zip((x, y), (u, v))))
 
-    for term in [X**2, Y**2, Integer(1)]:
-        if term not in coeff.keys():
-            coeff[term] = Integer(0)
+    coeff = simplified.as_coefficients_dict()
 
-    return -coeff[Y**2]/coeff[X**2], -coeff[Integer(1)]/coeff[X**2]
+    return -coeff[Y**2]/coeff[X**2], -coeff[1]/coeff[X**2]
 
 
 def check_param(x, y, a, t):
     """
-    Check if there is a number modulo ``a`` such that ``x`` and ``y`` are both
-    integers. If exist, then find a parametric representation for ``x`` and
-    ``y``.
+    If there is a number modulo ``a`` such that ``x`` and ``y`` are both
+    integers, then return a parametric representation for ``x`` and ``y``
+    else return (None, None).
 
     Here ``x`` and ``y`` are functions of ``t``.
     """
-    k, m, n = symbols("k, m, n", integer=True)
-    p = Wild("p", exclude=[k])
-    q = Wild("q", exclude=[k])
-    ok = False
-
-    for i in range(a):
-
-        z_x = _mexpand(Subs(x, t, a*k + i).doit()).match(p*k + q)
-        z_y = _mexpand(Subs(y, t, a*k + i).doit()).match(p*k + q)
-
-        if (isinstance(z_x[p], Integer) and isinstance(z_x[q], Integer) and
-            isinstance(z_y[p], Integer) and isinstance(z_y[q], Integer)):
-            ok = True
-            break
-
-    if ok == True:
-
-        x_param = x.match(p*t + q)
-        y_param = y.match(p*t + q)
-
-        if x_param[p] == 0 or y_param[p] == 0:
-            if x_param[p] == 0:
-                l1, junk = Poly(y).clear_denoms()
-            else:
-                l1 = 1
+    from sympy.simplify.simplify import clear_coefficients
 
-            if y_param[p] == 0:
-                l2, junk = Poly(x).clear_denoms()
-            else:
-                l2 = 1
+    if x.is_number and not x.is_Integer:
+        return (None, None)
 
-            return x*ilcm(l1, l2), y*ilcm(l1, l2)
+    if y.is_number and not y.is_Integer:
+        return (None, None)
 
-        eq = S(m - x_param[q])/x_param[p] - S(n - y_param[q])/y_param[p]
+    m, n = symbols("m, n", integer=True)
+    c, p = (m*x + n*y).as_content_primitive()
+    if a % c.q:
+        return (None, None)
 
-        lcm_denom, junk = Poly(eq).clear_denoms()
-        eq = eq * lcm_denom
+    # clear_coefficients(mx + b, R)[1] -> (R - b)/m
+    eq = clear_coefficients(x, m)[1] - clear_coefficients(y, n)[1]
+    junk, eq = eq.as_content_primitive()
 
-        return diop_solve(eq, t)[0], diop_solve(eq, t)[1]
-    else:
-        return (None, None)
+    return diop_solve(eq, t)
 
 
 def diop_ternary_quadratic(eq):
@@ -1732,18 +1730,18 @@ def diop_ternary_quadratic(eq):
     >>> diop_ternary_quadratic(x**2 - 49*y**2 - z**2 + 13*z*y -8*x*y)
     (9, 1, 5)
     """
-    var, coeff, diop_type = classify_diop(eq)
+    var, coeff, diop_type = classify_diop(eq, _dict=False)
 
-    if diop_type == "homogeneous_ternary_quadratic":
+    if diop_type in (
+            "homogeneous_ternary_quadratic",
+            "homogeneous_ternary_quadratic_normal"):
         return _diop_ternary_quadratic(var, coeff)
 
 
 def _diop_ternary_quadratic(_var, coeff):
 
-    x, y, z = _var[:3]
-
-    var = [x]*3
-    var[0], var[1], var[2] = _var[0], _var[1], _var[2]
+    x, y, z = _var
+    var = [x, y, z]
 
     # Equations of the form B*x*y + C*z*x + E*y*z = 0 and At least two of the
     # coefficients A, B, C are non-zero.
@@ -1754,8 +1752,8 @@ def _diop_ternary_quadratic(_var, coeff):
     # using methods for binary quadratic diophantine equations. Let's select the
     # solution which minimizes |x| + |z|
 
-    if coeff[x**2] == 0 and coeff[y**2] == 0 and coeff[z**2] == 0:
-        if coeff[x*z] != 0:
+    if not any(coeff[i**2] for i in var):
+        if coeff[x*z]:
             sols = diophantine(coeff[x*y]*x + coeff[y*z]*z - x*z)
             s = sols.pop()
             min_sum = abs(s[0]) + abs(s[1])
@@ -1771,7 +1769,7 @@ def _diop_ternary_quadratic(_var, coeff):
             var[0], var[1] = _var[1], _var[0]
             y_0, x_0, z_0 = _diop_ternary_quadratic(var, coeff)
 
-        return simplified(x_0, y_0, z_0)
+        return _remove_gcd(x_0, y_0, z_0)
 
     if coeff[x**2] == 0:
         # If the coefficient of x is zero change the variables
@@ -1784,7 +1782,7 @@ def _diop_ternary_quadratic(_var, coeff):
             y_0, x_0, z_0 = _diop_ternary_quadratic(var, coeff)
 
     else:
-        if coeff[x*y] != 0 or coeff[x*z] != 0:
+        if coeff[x*y] or coeff[x*z]:
         # Apply the transformation x --> X - (B*y + C*z)/(2*A)
             A = coeff[x**2]
             B = coeff[x*y]
@@ -1804,11 +1802,11 @@ def _diop_ternary_quadratic(_var, coeff):
 
             X_0, y_0, z_0 = _diop_ternary_quadratic(var, _coeff)
 
-            if X_0 == None:
+            if X_0 is None:
                 return (None, None, None)
 
-            l = (S(B*y_0 + C*z_0)/(2*A)).q
-            x_0, y_0, z_0 = X_0*l - (S(B*y_0 + C*z_0)/(2*A)).p, y_0*l, z_0*l
+            p, q = _rational_pq(B*y_0 + C*z_0, 2*A)
+            x_0, y_0, z_0 = X_0*q - p, y_0*q, z_0*q
 
         elif coeff[z*y] != 0:
             if coeff[y**2] == 0:
@@ -1817,8 +1815,7 @@ def _diop_ternary_quadratic(_var, coeff):
                     A = coeff[x**2]
                     E = coeff[y*z]
 
-                    b = (S(-E)/A).p
-                    a = (S(-E)/A).q
+                    b, a = _rational_pq(-E, A)
 
                     x_0, y_0, z_0 = b, a, b
 
@@ -1836,32 +1833,44 @@ def _diop_ternary_quadratic(_var, coeff):
             # Ax**2 + D*y**2 + F*z**2 = 0, C may be zero
             x_0, y_0, z_0 = _diop_ternary_quadratic_normal(var, coeff)
 
-    return simplified(x_0, y_0, z_0)
+    return _remove_gcd(x_0, y_0, z_0)
 
 
 def transformation_to_normal(eq):
     """
-    Returns the transformation Matrix from general ternary quadratic equation
-    `eq` to normal form.
-
-    General form of the ternary quadratic equation is `ax^2 + by^2 cz^2 + dxy +
-    eyz + fxz`. This function returns a 3X3 transformation Matrix which
-    transforms the former equation to the form `ax^2 + by^2 + cz^2 = 0`. This
-    is not used in solving ternary quadratics. Only implemented for the sake
-    of completeness.
+    Returns the transformation Matrix that converts a general ternary
+    quadratic equation `eq` (`ax^2 + by^2 + cz^2 + dxy + eyz + fxz`)
+    to a form without cross terms: `ax^2 + by^2 + cz^2 = 0`. This is
+    not used in solving ternary quadratics; it is only implemented for
+    the sake of completeness.
     """
-    var, coeff, diop_type = classify_diop(eq)
+    var, coeff, diop_type = classify_diop(eq, _dict=False)
 
-    if diop_type == "homogeneous_ternary_quadratic":
+    if diop_type in (
+            "homogeneous_ternary_quadratic",
+            "homogeneous_ternary_quadratic_normal"):
         return _transformation_to_normal(var, coeff)
 
 
 def _transformation_to_normal(var, coeff):
 
-    _var = [var[0]]*3
-    _var[1], _var[2] = var[1], var[2]
-
-    x, y, z = var[:3]
+    _var = list(var)  # copy
+    x, y, z = var
+
+    if not any(coeff[i**2] for i in var):
+        # https://math.stackexchange.com/questions/448051/transform-quadratic-ternary-form-to-normal-form/448065#448065
+        a = coeff[x*y]
+        b = coeff[y*z]
+        c = coeff[x*z]
+        swap = False
+        if not a:  # b can't be 0 or else there aren't 3 vars
+            swap = True
+            a, b = b, a
+        T = Matrix(((1, 1, -b/a), (1, -1, -c/a), (0, 0, 1)))
+        if swap:
+            T.row_swap(0, 1)
+            T.col_swap(0, 1)
+        return T
 
     if coeff[x**2] == 0:
         # If the coefficient of x is zero change the variables
@@ -1879,65 +1888,52 @@ def _transformation_to_normal(var, coeff):
             T.col_swap(0, 1)
             return T
 
-    else:
-        # Apply the transformation x --> X - (B*Y + C*Z)/(2*A)
-        if coeff[x*y] != 0 or coeff[x*z] != 0:
-            A = coeff[x**2]
-            B = coeff[x*y]
-            C = coeff[x*z]
-            D = coeff[y**2]
-            E = coeff[y*z]
-            F = coeff[z**2]
-
-            _coeff = dict()
+    # Apply the transformation x --> X - (B*Y + C*Z)/(2*A)
+    if coeff[x*y] != 0 or coeff[x*z] != 0:
+        A = coeff[x**2]
+        B = coeff[x*y]
+        C = coeff[x*z]
+        D = coeff[y**2]
+        E = coeff[y*z]
+        F = coeff[z**2]
 
-            _coeff[x**2] = 4*A**2
-            _coeff[y**2] = 4*A*D - B**2
-            _coeff[z**2] = 4*A*F - C**2
-            _coeff[y*z] = 4*A*E - 2*B*C
-            _coeff[x*y] = 0
-            _coeff[x*z] = 0
+        _coeff = dict()
 
-            T_0 = _transformation_to_normal(_var, _coeff)
-            return Matrix(3, 3, [1, S(-B)/(2*A), S(-C)/(2*A), 0, 1, 0, 0, 0, 1]) * T_0
+        _coeff[x**2] = 4*A**2
+        _coeff[y**2] = 4*A*D - B**2
+        _coeff[z**2] = 4*A*F - C**2
+        _coeff[y*z] = 4*A*E - 2*B*C
+        _coeff[x*y] = 0
+        _coeff[x*z] = 0
 
-        elif coeff[y*z] != 0:
-            if coeff[y**2] == 0:
-                if coeff[z**2] == 0:
-                    # Equations of the form A*x**2 + E*yz = 0.
-                    # Apply transformation y -> Y + Z ans z -> Y - Z
-                    return Matrix(3, 3, [1, 0, 0, 0, 1, 1, 0, 1, -1])
+        T_0 = _transformation_to_normal(_var, _coeff)
+        return Matrix(3, 3, [1, S(-B)/(2*A), S(-C)/(2*A), 0, 1, 0, 0, 0, 1])*T_0
 
-                else:
-                    # Ax**2 + E*y*z + F*z**2  = 0
-                    _var[0], _var[2] = var[2], var[0]
-                    T = _transformtion_to_normal(_var, coeff)
-                    T.row_swap(0, 2)
-                    T.col_swap(0, 2)
-                    return T
+    elif coeff[y*z] != 0:
+        if coeff[y**2] == 0:
+            if coeff[z**2] == 0:
+                # Equations of the form A*x**2 + E*yz = 0.
+                # Apply transformation y -> Y + Z ans z -> Y - Z
+                return Matrix(3, 3, [1, 0, 0, 0, 1, 1, 0, 1, -1])
 
             else:
-                # A*x**2 + D*y**2 + E*y*z + F*z**2 = 0, F may be zero
-                _var[0], _var[1] = var[1], var[0]
+                # Ax**2 + E*y*z + F*z**2  = 0
+                _var[0], _var[2] = var[2], var[0]
                 T = _transformation_to_normal(_var, coeff)
-                T.row_swap(0, 1)
-                T.col_swap(0, 1)
+                T.row_swap(0, 2)
+                T.col_swap(0, 2)
                 return T
 
         else:
-            return Matrix(3, 3, [1, 0, 0, 0, 1, 0, 0, 0, 1])
-
-
-def simplified(x, y, z):
-    """
-    Simplify the solution `(x, y, z)`.
-    """
-    if x == None or y == None or z == None:
-        return (x, y, z)
-
-    g = igcd(x, igcd(y, z))
+            # A*x**2 + D*y**2 + E*y*z + F*z**2 = 0, F may be zero
+            _var[0], _var[1] = var[1], var[0]
+            T = _transformation_to_normal(_var, coeff)
+            T.row_swap(0, 1)
+            T.col_swap(0, 1)
+            return T
 
-    return x // g, y // g, z // g
+    else:
+        return Matrix.eye(3)
 
 
 def parametrize_ternary_quadratic(eq):
@@ -1969,41 +1965,48 @@ def parametrize_ternary_quadratic(eq):
            Press, Cambridge, 1998.
 
     """
-    var, coeff, diop_type = classify_diop(eq)
+    var, coeff, diop_type = classify_diop(eq, _dict=False)
 
-    if diop_type == "homogeneous_ternary_quadratic":
+    if diop_type in (
+            "homogeneous_ternary_quadratic",
+            "homogeneous_ternary_quadratic_normal"):
         x_0, y_0, z_0 = _diop_ternary_quadratic(var, coeff)
-        return _parametrize_ternary_quadratic((x_0, y_0, z_0), var, coeff)
+        return _parametrize_ternary_quadratic(
+            (x_0, y_0, z_0), var, coeff)
 
 
 def _parametrize_ternary_quadratic(solution, _var, coeff):
+    # called for a*x**2 + b*y**2 + c*z**2 + d*x*y + e*y*z + f*x*z = 0
+    assert 1 not in coeff
 
-    x, y, z = _var[:3]
+    x, y, z = _var
 
-    x_0, y_0, z_0 = solution[:3]
+    x_0, y_0, z_0 = solution
 
-    v = [x]*3
-    v[0], v[1], v[2] = _var[0], _var[1], _var[2]
+    v = list(_var)  # copy
 
-    if x_0 == None:
+    if x_0 is None:
+        return (None, None, None)
+
+    if solution.count(0) >= 2:
+        # if there are 2 zeros the the equation reduces
+        # to k*X**2 == 0 where X is x, y, or z so X must
+        # be zero, too. So there is only the trivial
+        # solution.
         return (None, None, None)
 
     if x_0 == 0:
-        if y_0 == 0:
-            v[0], v[2] = v[2], v[0]
-            z_p, y_p, x_p = _parametrize_ternary_quadratic((z_0, y_0, x_0), v, coeff)
-            return x_p, y_p, z_p
-        else:
-            v[0], v[1] = v[1], v[0]
-            y_p, x_p, z_p = _parametrize_ternary_quadratic((y_0, x_0, z_0), v, coeff)
-            return x_p, y_p, z_p
+        v[0], v[1] = v[1], v[0]
+        y_p, x_p, z_p = _parametrize_ternary_quadratic(
+            (y_0, x_0, z_0), v, coeff)
+        return x_p, y_p, z_p
 
-    x, y, z = v[:3]
+    x, y, z = v
     r, p, q = symbols("r, p, q", integer=True)
 
-    eq = x**2*coeff[x**2] + y**2*coeff[y**2] + z**2*coeff[z**2] + x*y*coeff[x*y] + y*z*coeff[y*z] + z*x*coeff[z*x]
-    eq_1 = Subs(eq, (x, y, z), (r*x_0, r*y_0 + p, r*z_0 + q)).doit()
-    eq_1 = _mexpand(eq_1)
+    eq = sum(k*v for k, v in coeff.items())
+    eq_1 = _mexpand(eq.subs(zip(
+        (x, y, z), (r*x_0, r*y_0 + p, r*z_0 + q))))
     A, B = eq_1.as_independent(r, as_Add=True)
 
 
@@ -2042,38 +2045,29 @@ def diop_ternary_quadratic_normal(eq):
     >>> diop_ternary_quadratic_normal(34*x**2 - 3*y**2 - 301*z**2)
     (4, 9, 1)
     """
-    var, coeff, diop_type = classify_diop(eq)
-
-    if diop_type == "homogeneous_ternary_quadratic":
+    var, coeff, diop_type = classify_diop(eq, _dict=False)
+    if diop_type == "homogeneous_ternary_quadratic_normal":
         return _diop_ternary_quadratic_normal(var, coeff)
 
 
 def _diop_ternary_quadratic_normal(var, coeff):
 
-    x, y, z = var[:3]
+    x, y, z = var
 
     a = coeff[x**2]
     b = coeff[y**2]
     c = coeff[z**2]
-
-    if a*b*c == 0:
-        raise ValueError("Try factoring out you equation or using diophantine()")
-
-    g = igcd(a, igcd(b, c))
-
-    a = a // g
-    b = b // g
-    c = c // g
-
-    a_0 = square_factor(a)
-    b_0 = square_factor(b)
-    c_0 = square_factor(c)
-
-    a_1 = a // a_0**2
-    b_1 = b // b_0**2
-    c_1 = c // c_0**2
-
-    a_2, b_2, c_2 = pairwise_prime(a_1, b_1, c_1)
+    try:
+        assert len([k for k in coeff if coeff[k]]) == 3
+        assert all(coeff[i**2] for i in var)
+    except AssertionError:
+        raise ValueError(filldedent('''
+    coeff dict is not consistent with assumption of this routine:
+    coefficients should be those of an expression in the form
+    a*x**2 + b*y**2 + c*z**2 where a*b*c != 0.'''))
+
+    (sqf_of_a, sqf_of_b, sqf_of_c), (a_1, b_1, c_1), (a_2, b_2, c_2) = \
+        sqf_normal(a, b, c, steps=True)
 
     A = -a_2*c_2
     B = -b_2*c_2
@@ -2082,20 +2076,19 @@ def _diop_ternary_quadratic_normal(var, coeff):
     if A < 0 and B < 0:
         return (None, None, None)
 
-    if (sqrt_mod(-b_2*c_2, a_2) == None or sqrt_mod(-c_2*a_2, b_2) == None or
-        sqrt_mod(-a_2*b_2, c_2) == None):
+    if (
+            sqrt_mod(-b_2*c_2, a_2) is None or
+            sqrt_mod(-c_2*a_2, b_2) is None or
+            sqrt_mod(-a_2*b_2, c_2) is None):
         return (None, None, None)
 
     z_0, x_0, y_0 = descent(A, B)
 
-    if divisible(z_0, c_2) == True:
-        z_0 = z_0 // abs(c_2)
-    else:
-        x_0 = x_0*(S(z_0)/c_2).q
-        y_0 = y_0*(S(z_0)/c_2).q
-        z_0 = (S(z_0)/c_2).p
+    z_0, q = _rational_pq(z_0, abs(c_2))
+    x_0 *= q
+    y_0 *= q
 
-    x_0, y_0, z_0 = simplified(x_0, y_0, z_0)
+    x_0, y_0, z_0 = _remove_gcd(x_0, y_0, z_0)
 
     # Holzer reduction
     if sign(a) == sign(b):
@@ -2109,126 +2102,118 @@ def _diop_ternary_quadratic_normal(var, coeff):
     y_0 = reconstruct(a_1, c_1, y_0)
     z_0 = reconstruct(a_1, b_1, z_0)
 
-    l = ilcm(a_0, ilcm(b_0, c_0))
+    sq_lcm = ilcm(sqf_of_a, sqf_of_b, sqf_of_c)
 
-    x_0 = abs(x_0*l//a_0)
-    y_0 = abs(y_0*l//b_0)
-    z_0 = abs(z_0*l//c_0)
+    x_0 = abs(x_0*sq_lcm//sqf_of_a)
+    y_0 = abs(y_0*sq_lcm//sqf_of_b)
+    z_0 = abs(z_0*sq_lcm//sqf_of_c)
 
-    return simplified(x_0, y_0, z_0)
+    return _remove_gcd(x_0, y_0, z_0)
 
 
-def square_factor(a):
+def sqf_normal(a, b, c, steps=False):
     """
-    Returns an integer `c` s.t. `a = c^2k, \ c,k \in Z`. Here `k` is square
-    free.
+    Return `a', b', c'`, the coefficients of the square-free normal
+    form of `ax^2 + by^2 + cz^2 = 0`, where `a', b', c'` are pairwise
+    prime.  If `steps` is True then also return three tuples:
+    `sq`, `sqf`, and `(a', b', c')` where `sq` contains the square
+    factors of `a`, `b` and `c` after removing the `gcd(a, b, c)`;
+    `sqf` contains the values of `a`, `b` and `c` after removing
+    both the `gcd(a, b, c)` and the square factors.
+
+    The solutions for `ax^2 + by^2 + cz^2 = 0` can be
+    recovered from the solutions of `a'x^2 + b'y^2 + c'z^2 = 0`.
 
     Examples
     ========
 
-    >>> from sympy.solvers.diophantine import square_factor
-    >>> square_factor(24)
-    2
-    >>> square_factor(36)
-    6
-    >>> square_factor(1)
-    1
-    """
-    f = factorint(abs(a))
-    c = 1
-
-    for p, e in f.items():
-        c = c * p**(e//2)
-
-    return c
-
-
-def pairwise_prime(a, b, c):
-    """
-    Transform `ax^2 + by^2 + cz^2 = 0` into an equivalent equation
-    `a'x^2 + b'y^2 + c'z^2 = 0` where `a', b', c'` are pairwise relatively
-    prime.
+    >>> from sympy.solvers.diophantine import sqf_normal
+    >>> sqf_normal(2 * 3**2 * 5, 2 * 5 * 11, 2 * 7**2 * 11)
+    (11, 1, 5)
+    >>> sqf_normal(2 * 3**2 * 5, 2 * 5 * 11, 2 * 7**2 * 11, True)
+    ((3, 1, 7), (5, 55, 11), (11, 1, 5))
 
-    Returns a tuple containing `a', b', c'`. `\gcd(a, b, c)` should equal `1`
-    for this to work. The solutions for `ax^2 + by^2 + cz^2 = 0` can be
-    recovered from the solutions of `a'x^2 + b'y^2 + c'z^2 = 0`.
+    References
+    ==========
 
-    Examples
-    ========
+    .. [1] Legendre's Theorem, Legrange's Descent,
+           http://public.csusm.edu/aitken_html/notes/legendre.pdf
 
-    >>> from sympy.solvers.diophantine import pairwise_prime
-    >>> pairwise_prime(6, 15, 10)
-    (5, 2, 3)
 
     See Also
     ========
 
-    make_prime(), reocnstruct()
-    """
-    a, b, c = make_prime(a, b, c)
-    b, c, a = make_prime(b, c, a)
-    c, a, b = make_prime(c, a, b)
-
-    return a, b, c
+    reconstruct()
+    """
+    ABC = A, B, C = _remove_gcd(a, b, c)
+    sq = tuple(square_factor(i) for i in ABC)
+    sqf = A, B, C = tuple([i//j**2 for i,j in zip(ABC, sq)])
+    pc = igcd(A, B)
+    A /= pc
+    B /= pc
+    pa = igcd(B, C)
+    B /= pa
+    C /= pa
+    pb = igcd(A, C)
+    A /= pb
+    B /= pb
+
+    A *= pa
+    B *= pb
+    C *= pc
+
+    if steps:
+        return (sq, sqf, (A, B, C))
+    else:
+        return A, B, C
 
 
-def make_prime(a, b, c):
+def square_factor(a):
     """
-    Transform the equation `ax^2 + by^2 + cz^2 = 0` to an equivalent equation
-    `a'x^2 + b'y^2 + c'z^2 = 0` with `\gcd(a', b') = 1`.
-
-    Returns a tuple `(a', b', c')` which satisfies above conditions. Note that
-    in the returned tuple `\gcd(a', c')` and `\gcd(b', c')` can take any value.
+    Returns an integer `c` s.t. `a = c^2k, \ c,k \in Z`. Here `k` is square
+    free. `a` can be given as an integer or a dictionary of factors.
 
     Examples
     ========
 
-    >>> from sympy.solvers.diophantine import make_prime
-    >>> make_prime(4, 2, 7)
-    (2, 1, 14)
+    >>> from sympy.solvers.diophantine import square_factor
+    >>> square_factor(24)
+    2
+    >>> square_factor(-36*3)
+    6
+    >>> square_factor(1)
+    1
+    >>> square_factor({3: 2, 2: 1, -1: 1})  # -18
+    3
 
     See Also
     ========
-
-    pairwaise_prime(), reconstruct()
+    sympy.ntheory.factor_.core
     """
-    g = igcd(a, b)
-
-    if g != 1:
-        f = factorint(g)
-        for p, e in f.items():
-            a = a // p**e
-            b = b // p**e
+    f = a if isinstance(a, dict) else factorint(a)
+    return Mul(*[p**(e//2) for p, e in f.items()])
 
-            if e % 2 == 1:
-                c = p*c
 
-    return a, b, c
-
-
-def reconstruct(a, b, z):
+def reconstruct(A, B, z):
     """
     Reconstruct the `z` value of an equivalent solution of `ax^2 + by^2 + cz^2`
-    from the `z` value of a solution of a transformed version of the above
-    equation.
+    from the `z` value of a solution of the square-free normal form of the
+    equation, `a'*x^2 + b'*y^2 + c'*z^2`, where `a'`, `b'` and `c'` are square
+    free and `gcd(a', b', c') == 1`.
     """
-    g = igcd(a, b)
-
-    if g != 1:
-        f = factorint(g)
-        for p, e in f.items():
-            if e %2 == 0:
-                z = z*p**(e//2)
-            else:
-                z = z*p**((e//2)+1)
-
+    f = factorint(igcd(A, B))
+    for p, e in f.items():
+        if e != 1:
+            raise ValueError('a and b should be square-free')
+        z *= p
     return z
 
 
 def ldescent(A, B):
     """
-    Uses Lagrange's method to find a non trivial solution to
-    `w^2 = Ax^2 + By^2`.
+    Return a non-trivial solution to `w^2 = Ax^2 + By^2` using
+    Lagrange's method; return None if there is no such solution.
+    .
 
     Here, `A \\neq 0` and `B \\neq 0` and `A` and `B` are square free. Output a
     tuple `(w_0, x_0, y_0)` which is a solution to the above equation.
@@ -2255,17 +2240,21 @@ def ldescent(A, B):
            London Mathematical Society Student Texts 41, Cambridge University
            Press, Cambridge, 1998.
     .. [2] Efficient Solution of Rational Conices, J. E. Cremona and D. Rusin,
-           Mathematics of Computation, Volume 00, Number 0.
+           [online], Available:
+           http://eprints.nottingham.ac.uk/60/1/kvxefz87.pdf
     """
     if abs(A) > abs(B):
         w, y, x = ldescent(B, A)
         return w, x, y
 
     if A == 1:
-        return (S.One, S.One, 0)
+        return (1, 1, 0)
 
     if B == 1:
-        return (S.One, 0, S.One)
+        return (1, 0, 1)
+
+    if B == -1:  # and A == -1
+        return
 
     r = sqrt_mod(A, B)
 
@@ -2279,27 +2268,24 @@ def ldescent(A, B):
         B_0 = None
 
         for i in div:
-            if isinstance(sqrt(abs(Q) // i), Integer):
-                B_0, d = sign(Q)*i, sqrt(abs(Q) // i)
+            sQ, _exact = integer_nthroot(abs(Q) // i, 2)
+            if _exact:
+                B_0, d = sign(Q)*i, sQ
                 break
 
-    if B_0 != None:
+    if B_0 is not None:
         W, X, Y = ldescent(A, B_0)
-        return simplified((-A*X + r*W), (r*X - W), Y*(B_0*d))
-    # In this module Descent will always be called with inputs which have solutions.
+        return _remove_gcd((-A*X + r*W), (r*X - W), Y*(B_0*d))
 
 
 def descent(A, B):
     """
-    Lagrange's `descent()` with lattice-reduction to find solutions to
-    `x^2 = Ay^2 + Bz^2`.
+    Returns a non-trivial solution, (x, y, z), to `x^2 = Ay^2 + Bz^2`
+    using Lagrange's descent method with lattice-reduction. `A` and `B`
+    are assumed to be valid for such a solution to exist.
 
-    Here `A` and `B` should be square free and pairwise prime. Always should be
-    called with suitable ``A`` and ``B`` so that the above equation has
-    solutions.
-
-    This is more faster than the normal Lagrange's descent algorithm because
-    the gaussian reduction is used.
+    This is faster than the normal Lagrange's descent algorithm because
+    the Gaussian reduction is used.
 
     Examples
     ========
@@ -2327,8 +2313,6 @@ def descent(A, B):
         return (1, 0, 1)
     if A == 1:
         return (1, 1, 0)
-    if B == -1:
-        return (None, None, None)
     if B == -A:
         return (0, 1, 1)
     if B == A:
@@ -2344,7 +2328,7 @@ def descent(A, B):
 
     x_1, z_1, y_1 = descent(A, t_1)
 
-    return simplified(x_0*x_1 + A*z_0*z_1, z_0*x_1 + x_0*z_1, t_1*t_2*y_1)
+    return _remove_gcd(x_0*x_1 + A*z_0*z_1, z_0*x_1 + x_0*z_1, t_1*t_2*y_1)
 
 
 def gaussian_reduce(w, a, b):
@@ -2361,9 +2345,9 @@ def gaussian_reduce(w, a, b):
     ==========
 
     .. [1] Gaussian lattice Reduction [online]. Available:
-        http://home.ie.cuhk.edu.hk/~wkshum/wordpress/?p=404
+           http://home.ie.cuhk.edu.hk/~wkshum/wordpress/?p=404
     .. [2] Efficient Solution of Rational Conices, J. E. Cremona and D. Rusin,
-        Mathematics of Computation, Volume 00, Number 0.
+           Mathematics of Computation, Volume 00, Number 0.
     """
     u = (0, 1)
     v = (1, 0)
@@ -2394,8 +2378,8 @@ def dot(u, v, w, a, b):
     `v = (v_{1}, v_{2})` which is defined in order to reduce solution of
     the congruence equation `X^2 - aZ^2 \equiv 0 \ (mod \ b)`.
     """
-    u_1, u_2 = u[:2]
-    v_1, v_2 = v[:2]
+    u_1, u_2 = u
+    v_1, v_2 = v
     return (w*u_1 + b*u_2)*(w*v_1 + b*v_2) + abs(a)*u_1*v_1
 
 
@@ -2405,39 +2389,72 @@ def norm(u, w, a, b):
     defined by `u \cdot v = (wu_{1} + bu_{2})(w*v_{1} + bv_{2}) + |a|*u_{1}*v_{1}`
     where `u = (u_{1}, u_{2})` and `v = (v_{1}, v_{2})`.
     """
-    u_1, u_2 = u[:2]
+    u_1, u_2 = u
     return sqrt(dot((u_1, u_2), (u_1, u_2), w, a, b))
 
 
-def holzer(x_0, y_0, z_0, a, b, c):
-    """
-    Simplify the solution `(x_{0}, y_{0}, z_{0})` of the equation
-    `ax^2 + by^2 = cz^2` with `a, b, c > 0` and `z_{0}^2 \geq \mid ab \mid` to
-    a new reduced solution `(x, y, z)` such that `z^2 \leq \mid ab \mid`.
+def holzer(x, y, z, a, b, c):
     """
-    while z_0 > sqrt(a*b):
+    Simplify the solution `(x, y, z)` of the equation
+    `ax^2 + by^2 = cz^2` with `a, b, c > 0` and `z^2 \geq \mid ab \mid` to
+    a new reduced solution `(x', y', z')` such that `z'^2 \leq \mid ab \mid`.
 
-        if c % 2 == 0:
-            k = c // 2
-            u_0, v_0 = base_solution_linear(k, y_0, -x_0)
+    The algorithm is an interpretation of Mordell's reduction as described
+    on page 8 of Cremona and Rusin's paper [1]_ and the work of Mordell in
+    reference [2]_.
 
-        else:
-            k = 2*c
-            u_0, v_0 = base_solution_linear(c, y_0, -x_0)
-
-        w = -(a*u_0*x_0 + b*v_0*y_0) // (c*z_0)
+    References
+    ==========
 
-        if c % 2 == 1:
-            if w % 2 != (a*u_0 + b*v_0) % 2:
-                w = w + 1
+    .. [1] Efficient Solution of Rational Conices, J. E. Cremona and D. Rusin,
+           Mathematics of Computation, Volume 00, Number 0.
+    .. [2] Diophantine Equations, L. J. Mordell, page 48.
 
-        x = (x_0*(a*u_0**2 + b*v_0**2 + c*w**2) - 2*u_0*(a*u_0*x_0 + b*v_0*y_0 + c*w*z_0)) // k
-        y = (y_0*(a*u_0**2 + b*v_0**2 + c*w**2) - 2*v_0*(a*u_0*x_0 + b*v_0*y_0 + c*w*z_0)) // k
-        z = (z_0*(a*u_0**2 + b*v_0**2 + c*w**2) - 2*w*(a*u_0*x_0 + b*v_0*y_0 + c*w*z_0)) // k
+    """
 
+    if _odd(c):
+        k = 2*c
+    else:
+        k = c//2
+
+    small = a*b*c
+    step = 0
+    while True:
+        t1, t2, t3 = a*x**2, b*y**2, c*z**2
+        # check that it's a solution
+        if t1 + t2 != t3:
+            if step == 0:
+                raise ValueError('bad starting solution')
+            break
         x_0, y_0, z_0 = x, y, z
+        if max(t1, t2, t3) <= small:
+            # Holzer condition
+            break
 
-    return x_0, y_0, z_0
+        uv = u, v = base_solution_linear(k, y_0, -x_0)
+        if None in uv:
+            break
+
+        p, q = -(a*u*x_0 + b*v*y_0), c*z_0
+        r = Rational(p, q)
+        if _even(c):
+            w = _nint_or_floor(p, q)
+            assert abs(w - r) <= S.Half
+        else:
+            w = p//q  # floor
+            if _odd(a*u + b*v + c*w):
+                w += 1
+            assert abs(w - r) <= S.One
+
+        A = (a*u**2 + b*v**2 + c*w**2)
+        B = (a*u*x_0 + b*v*y_0 + c*w*z_0)
+        x = Rational(x_0*A - 2*u*B, k)
+        y = Rational(y_0*A - 2*v*B, k)
+        z = Rational(z_0*A - 2*w*B, k)
+        assert all(i.is_Integer for i in (x, y, z))
+        step += 1
+
+    return tuple([int(i) for i in (x_0, y_0, z_0)])
 
 
 def diop_general_pythagorean(eq, param=symbols("m", integer=True)):
@@ -2465,7 +2482,7 @@ def diop_general_pythagorean(eq, param=symbols("m", integer=True)):
     >>> diop_general_pythagorean(9*a**2 - 4*b**2 + 16*c**2 + 25*d**2 + e**2)
     (10*m1**2  + 10*m2**2  + 10*m3**2 - 10*m4**2, 15*m1**2  + 15*m2**2  + 15*m3**2  + 15*m4**2, 15*m1*m4, 12*m2*m4, 60*m3*m4)
     """
-    var, coeff, diop_type  = classify_diop(eq)
+    var, coeff, diop_type  = classify_diop(eq, _dict=False)
 
     if diop_type == "general_pythagorean":
         return _diop_general_pythagorean(var, coeff, param)
@@ -2475,7 +2492,7 @@ def _diop_general_pythagorean(var, coeff, t):
 
     if sign(coeff[var[0]**2]) + sign(coeff[var[1]**2]) + sign(coeff[var[2]**2]) < 0:
         for key in coeff.keys():
-            coeff[key] = coeff[key] * -1
+            coeff[key] = -coeff[key]
 
     n = len(var)
     index = 0
@@ -2484,26 +2501,19 @@ def _diop_general_pythagorean(var, coeff, t):
         if sign(coeff[v**2]) == -1:
             index = i
 
-    m = symbols(str(t) + "1:" + str(n), integer=True)
-    l = []
-    ith = 0
-
-    for m_i in m:
-        ith = ith + m_i**2
-
-    l.append(ith - 2*m[n - 2]**2)
-
-    for i in range(n - 2):
-        l.append(2*m[i]*m[n-2])
-
-    sol = l[:index] + [ith] + l[index:]
+    m = symbols('%s1:%i' % (t, n), integer=True)
+    ith = sum(m_i**2 for m_i in m)
+    L = [ith - 2*m[n - 2]**2]
+    L.extend([2*m[i]*m[n-2] for i in range(n - 2)])
+    sol = L[:index] + [ith] + L[index:]
 
     lcm = 1
     for i, v in enumerate(var):
         if i == index or (index > 0 and i == 0) or (index == 0 and i == 1):
             lcm = ilcm(lcm, sqrt(abs(coeff[v**2])))
         else:
-            lcm = ilcm(lcm, sqrt(coeff[v**2]) if sqrt(coeff[v**2]) % 2 else sqrt(coeff[v**2]) // 2)
+            s = sqrt(coeff[v**2])
+            lcm = ilcm(lcm, s if _odd(s) else s//2)
 
     for i, v in enumerate(var):
         sol[i] = (lcm*sol[i]) / sqrt(abs(coeff[v**2]))
@@ -2515,17 +2525,14 @@ def diop_general_sum_of_squares(eq, limit=1):
     """
     Solves the equation `x_{1}^2 + x_{2}^2 + . . . + x_{n}^2 - k = 0`.
 
-    Returns at most ``limit`` number of solutions. Currently there is no way to
-    set ``limit`` using higher level API's like ``diophantine()`` or
-    ``diop_solve()`` but that will be fixed soon.
+    Returns at most ``limit`` number of solutions.
 
     Usage
     =====
 
     ``general_sum_of_squares(eq, limit)`` : Here ``eq`` is an expression which
     is assumed to be zero. Also, ``eq`` should be in the form,
-    `x_{1}^2 + x_{2}^2 + . . . + x_{n}^2 - k = 0`. At most ``limit`` number of
-    solutions are returned.
+    `x_{1}^2 + x_{2}^2 + . . . + x_{n}^2 - k = 0`.
 
     Details
     =======
@@ -2539,60 +2546,109 @@ def diop_general_sum_of_squares(eq, limit=1):
     >>> from sympy.solvers.diophantine import diop_general_sum_of_squares
     >>> from sympy.abc import a, b, c, d, e, f
     >>> diop_general_sum_of_squares(a**2 + b**2 + c**2 + d**2 + e**2 - 2345)
-    set([(0, 48, 5, 4, 0)])
+    set([(2, 4, 4, 10, 47)])
 
     Reference
     =========
 
-    .. [1] Representing an Integer as a sum of three squares, [online],
+    .. [1] Representing an integer as a sum of three squares, [online],
         Available:
         http://www.proofwiki.org/wiki/Integer_as_Sum_of_Three_Squares
     """
-    var, coeff, diop_type = classify_diop(eq)
+    var, coeff, diop_type = classify_diop(eq, _dict=False)
 
     if diop_type == "general_sum_of_squares":
-        return _diop_general_sum_of_squares(var, coeff, limit)
-
+        return _diop_general_sum_of_squares(var, -coeff[1], limit)
 
-def _diop_general_sum_of_squares(var, coeff, limit=1):
 
+def _diop_general_sum_of_squares(var, k, limit=1):
+    # solves Eq(sum(i**2 for i in var), k)
     n = len(var)
-    k = -int(coeff[Integer(1)])
-    s = set([])
+    if n < 3:
+        raise ValueError('n must be greater than 2')
 
-    if k < 0:
-        return set([])
+    s = set()
 
-    if n == 3:
-        s.add(sum_of_three_squares(k))
-    elif n == 4:
-        s.add(sum_of_four_squares(k))
-    else:
+    if k < 0 or limit < 1:
+        return s
+
+    sign = [-1 if x.is_nonpositive else 1 for x in var]
+    negs = sign.count(-1) != 0
 
-        m = n // 4
-        f = partition(k, m, True)
+    took = 0
+    for t in sum_of_squares(k, n, zeros=True):
+        if negs:
+            s.add(tuple([sign[i]*j for i, j in enumerate(t)]))
+        else:
+            s.add(t)
+        took += 1
+        if took == limit:
+            break
+    return s
 
-        for j in range(limit):
 
-            soln = []
-            try:
-                l = next(f)
-            except StopIteration:
-                break
+def diop_general_sum_of_even_powers(eq, limit=1):
+    """
+    Solves the equation `x_{1}^e + x_{2}^e + . . . + x_{n}^e - k = 0`
+    where `e` is an even, integer power.
+
+    Returns at most ``limit`` number of solutions.
+
+    Usage
+    =====
+
+    ``general_sum_of_even_powers(eq, limit)`` : Here ``eq`` is an expression which
+    is assumed to be zero. Also, ``eq`` should be in the form,
+    `x_{1}^e + x_{2}^e + . . . + x_{n}^e - k = 0`.
+
+    Examples
+    ========
+
+    >>> from sympy.solvers.diophantine import diop_general_sum_of_even_powers
+    >>> from sympy.abc import a, b
+    >>> diop_general_sum_of_even_powers(a**4 + b**4 - (2**4 + 3**4))
+    set([(2, 3)])
+
+    See Also
+    ========
+    power_representation()
+    """
+    var, coeff, diop_type = classify_diop(eq, _dict=False)
+
+    if diop_type == "general_sum_of_even_powers":
+        for k in coeff.keys():
+            if k.is_Pow and coeff[k]:
+                p = k.exp
+        return _diop_general_sum_of_even_powers(var, p, -coeff[1], limit)
 
-            for n_i in l:
-                a, b, c, d = sum_of_four_squares(n_i)
-                soln = soln + [a, b, c, d]
 
-            soln = soln + [0] * (n % 4)
+def _diop_general_sum_of_even_powers(var, p, n, limit=1):
+    # solves Eq(sum(i**2 for i in var), n)
+    k = len(var)
 
-            s.add(tuple(soln))
+    s = set()
 
+    if n < 0 or limit < 1:
+        return s
+
+    sign = [-1 if x.is_nonpositive else 1 for x in var]
+    negs = sign.count(-1) != 0
+
+    took = 0
+    for t in power_representation(n, p, k):
+        if negs:
+            s.add(tuple([sign[i]*j for i, j in enumerate(t)]))
+        else:
+            s.add(t)
+        took += 1
+        if took == limit:
+            break
     return s
 
 
-## Functions below this comment can be more suitably grouped under an Additive number theory module
-## rather than the Diophantine equation module.
+## Functions below this comment can be more suitably grouped under
+## an Additive number theory module rather than the Diophantine
+## equation module.
 
 
 def partition(n, k=None, zeros=False):
@@ -2600,12 +2656,11 @@ def partition(n, k=None, zeros=False):
     Returns a generator that can be used to generate partitions of an integer
     `n`.
 
-    A partition of `n` is a set of positive integers which add upto `n`. For
-    example, partitions of 3 are 3 , 1 + 2, 1 + 1+ 1. A partition is returned
+    A partition of `n` is a set of positive integers which add up to `n`. For
+    example, partitions of 3 are 3, 1 + 2, 1 + 1 + 1. A partition is returned
     as a tuple. If ``k`` equals None, then all possible partitions are returned
     irrespective of their size, otherwise only the partitions of size ``k`` are
-    returned. If there are no partions of `n` with size `k` then an empty tuple
-    is returned. If the ``zero`` parameter is set to True then a suitable
+    returned. If the ``zero`` parameter is set to True then a suitable
     number of zeros are added at the end of every partition of size less than
     ``k``.
 
@@ -2631,101 +2686,51 @@ def partition(n, k=None, zeros=False):
     (1, 1, 1, 2)
     >>> g = partition(5, 3)
     >>> next(g)
-    (3, 1, 1)
+    (1, 1, 3)
     >>> next(g)
-    (2, 2, 1)
-
-    Reference
-    =========
+    (1, 2, 2)
+    >>> g = partition(5, 3, zeros=True)
+    >>> next(g)
+    (0, 0, 5)
 
-    .. [1] Generating Integer Partitions, [online],
-        Available: http://jeromekelleher.net/partitions.php
     """
-    if n < 1:
-        yield tuple()
-
-    if k is not None:
-        if k < 1:
-            yield tuple()
-
-        elif k > n:
-            if zeros:
-                for i in range(1, n):
-                    for t in partition(n, i):
-                        yield (t,) + (0,) * (k - i)
-            else:
-                yield tuple()
-
-        else:
-            a = [1 for i in range(k)]
-            a[0] = n - k + 1
-
-            yield tuple(a)
-
-            i = 1
-            while a[0] >= n // k + 1:
-                j = 0
-
-                while j < i and j + 1 < k:
-                    a[j] = a[j] - 1
-                    a[j + 1] = a[j + 1] + 1
-
-                    yield tuple(a)
-
-                    j = j + 1
-
-                i = i + 1
-
-            if zeros:
-                for m in range(1, k):
-                    for a in partition(n, m):
-                        yield tuple(a) + (0,) * (k - m)
-
+    from sympy.utilities.iterables import ordered_partitions
+    if not zeros or k is None:
+        for i in ordered_partitions(n, k):
+            yield tuple(i)
     else:
-        a = [0 for i in range(n + 1)]
-        l = 1
-        y = n - 1
-
-        while l != 0:
-            x = a[l - 1] + 1
-            l -= 1
-
-            while 2*x <= y:
-                a[l] = x
-                y -= x
-                l += 1
-
-            m = l + 1
-            while x <= y:
-                a[l] = x
-                a[m] = y
-                yield tuple(a[:l + 2])
-                x += 1
-                y -= 1
-
-            a[l] = x + y
-            y = x + y - 1
-            yield tuple(a[:l + 1])
+        for m in range(1, k + 1):
+            for i in ordered_partitions(n, m):
+                i = tuple(i)
+                yield (0,)*(k - len(i)) + i
 
 
 def prime_as_sum_of_two_squares(p):
     """
-    Represent a prime `p` which is congruent to 1 mod 4, as a sum of two
-    squares.
+    Represent a prime `p` as a unique sum of two squares; this can
+    only be done if the prime is congruent to 1 mod 4.
 
     Examples
     ========
 
     >>> from sympy.solvers.diophantine import prime_as_sum_of_two_squares
+    >>> prime_as_sum_of_two_squares(7)  # can't be done
     >>> prime_as_sum_of_two_squares(5)
-    (2, 1)
+    (1, 2)
 
     Reference
     =========
 
     .. [1] Representing a number as a sum of four squares, [online],
-        Available: http://www.schorn.ch/howto.html
+        Available: http://schorn.ch/lagrange.html
+
+    See Also
+    ========
+    sum_of_squares()
     """
+    if not p % 4 == 1:
+        return
+
     if p % 8 == 5:
         b = 2
     else:
@@ -2740,7 +2745,7 @@ def prime_as_sum_of_two_squares(p):
     while b**2 > p:
         a, b = b, a % b
 
-    return (b, a % b)
+    return (int(a % b), int(b))  # convert from long
 
 
 def sum_of_three_squares(n):
@@ -2748,7 +2753,7 @@ def sum_of_three_squares(n):
     Returns a 3-tuple `(a, b, c)` such that `a^2 + b^2 + c^2 = n` and
     `a, b, c \geq 0`.
 
-    Returns (None, None, None) if `n = 4^a(8m + 7)` for some `a, m \in Z`. See
+    Returns None if `n = 4^a(8m + 7)` for some `a, m \in Z`. See
     [1]_ for more details.
 
     Usage
@@ -2761,13 +2766,17 @@ def sum_of_three_squares(n):
 
     >>> from sympy.solvers.diophantine import sum_of_three_squares
     >>> sum_of_three_squares(44542)
-    (207, 37, 18)
+    (18, 37, 207)
 
     References
     ==========
 
     .. [1] Representing a number as a sum of three squares, [online],
-        Available: http://www.schorn.ch/howto.html
+        Available: http://schorn.ch/lagrange.html
+
+    See Also
+    ========
+    sum_of_squares()
     """
     special = {1:(1, 0, 0), 2:(1, 1, 0), 3:(1, 1, 1), 10: (1, 3, 0), 34: (3, 3, 4), 58:(3, 7, 0),
         85:(6, 7, 0), 130:(3, 11, 0), 214:(3, 6, 13), 226:(8, 9, 9), 370:(8, 9, 15),
@@ -2779,47 +2788,46 @@ def sum_of_three_squares(n):
     if n == 0:
         return (0, 0, 0)
 
-    while n % 4 == 0:
-        v = v + 1
-        n = n // 4
+    v = multiplicity(4, n)
+    n //= 4**v
 
     if n % 8 == 7:
-        return (None, None, None)
+        return
 
     if n in special.keys():
         x, y, z = special[n]
-        return (2**v*x, 2**v*y, 2**v*z)
+        return _sorted_tuple(2**v*x, 2**v*y, 2**v*z)
 
-    l = int(sqrt(n))
+    s, _exact = integer_nthroot(n, 2)
 
-    if n == l**2:
-        return (2**v*l, 0, 0)
+    if _exact:
+        return (2**v*s, 0, 0)
 
     x = None
 
     if n % 8 == 3:
-        l = l if l % 2 else l - 1
+        s = s if _odd(s) else s - 1
 
-        for i in range(l, -1, -2):
+        for i in range(s, -1, -2):
             if isprime((n - i**2) // 2):
                 x = i
                 break
 
         y, z = prime_as_sum_of_two_squares((n - x**2) // 2)
-        return (2**v*x, 2**v*(y + z), 2**v*abs(y - z))
+        return _sorted_tuple(2**v*x, 2**v*(y + z), 2**v*abs(y - z))
 
     if n % 8 == 2 or n % 8 == 6:
-        l = l if l % 2 else l - 1
+        s = s if _odd(s) else s - 1
     else:
-        l = l - 1 if l % 2 else l
+        s = s - 1 if _odd(s) else s
 
-    for i in range(l, -1, -2):
+    for i in range(s, -1, -2):
         if isprime(n - i**2):
             x = i
             break
 
     y, z = prime_as_sum_of_two_squares(n - x**2)
-    return (2**v*x, 2**v*y, 2**v*z)
+    return _sorted_tuple(2**v*x, 2**v*y, 2**v*z)
 
 
 def sum_of_four_squares(n):
@@ -2838,23 +2846,25 @@ def sum_of_four_squares(n):
 
     >>> from sympy.solvers.diophantine import sum_of_four_squares
     >>> sum_of_four_squares(3456)
-    (8, 48, 32, 8)
+    (8, 8, 32, 48)
     >>> sum_of_four_squares(1294585930293)
-    (0, 1137796, 2161, 1234)
+    (0, 1234, 2161, 1137796)
 
     References
     ==========
 
     .. [1] Representing a number as a sum of four squares, [online],
-        Available: http://www.schorn.ch/howto.html
+        Available: http://schorn.ch/lagrange.html
+
+    See Also
+    ========
+    sum_of_squares()
     """
     if n == 0:
         return (0, 0, 0, 0)
 
-    v = 0
-    while n % 4 == 0:
-        v = v + 1
-        n = n // 4
+    v = multiplicity(4, n)
+    n //= 4**v
 
     if n % 8 == 7:
         d = 2
@@ -2867,70 +2877,218 @@ def sum_of_four_squares(n):
 
     x, y, z = sum_of_three_squares(n)
 
-    return (2**v*d, 2**v*x, 2**v*y, 2**v*z)
+    return _sorted_tuple(2**v*d, 2**v*x, 2**v*y, 2**v*z)
 
 
 def power_representation(n, p, k, zeros=False):
     """
-    Returns a generator for finding k-tuples `(n_{1}, n_{2}, . . . n_{k})` such
-    that `n = n_{1}^p + n_{2}^p + . . . n_{k}^p`.
-
-    Here `n` is a non-negative integer. StopIteration exception is raised after
-    all the solutions are generated, so should always be used within a try-
-    catch block.
+    Returns a generator for finding k-tuples of integers,
+    `(n_{1}, n_{2}, . . . n_{k})`, such that
+    `n = n_{1}^p + n_{2}^p + . . . n_{k}^p`.
 
     Usage
     =====
 
-    ``power_representation(n, p, k, zeros)``: Represent number ``n`` as a sum
-    of ``k``, ``p``th powers. If ``zeros`` is true, then the solutions will
-    contain zeros.
+    ``power_representation(n, p, k, zeros)``: Represent non-negative number
+    ``n`` as a sum of ``k`` ``p``th powers. If ``zeros`` is true, then the
+    solutions is allowed to contain zeros.
 
     Examples
     ========
 
     >>> from sympy.solvers.diophantine import power_representation
-    >>> f = power_representation(1729, 3, 2) # Represent 1729 as a sum of two cubes
+
+    Represent 1729 as a sum of two cubes:
+
+    >>> f = power_representation(1729, 3, 2)
     >>> next(f)
-    (12, 1)
+    (1, 12)
     >>> next(f)
-    (10, 9)
+    (9, 10)
+
+    If the flag `zeros` is True, the solution may contain tuples with
+    zeros; any such solutions will be generated after the solutions
+    without zeros:
+
+    >>> list(power_representation(125, 2, 3, zeros=True))
+    [(3, 4, 10), (5, 6, 8), (0, 2, 11), (0, 5, 10)]
+
+    For even `p` the `permute_sign` function can be used to get all
+    signed values:
+
+    >>> from sympy.utilities.iterables import permute_signs
+    >>> list(permute_signs((1, 12)))
+    [(1, 12), (-1, 12), (1, -12), (-1, -12)]
+
+    All possible signed permutations can also be obtained:
+
+    >>> from sympy.utilities.iterables import signed_permutations
+    >>> list(signed_permutations((1, 12)))
+    [(1, 12), (-1, 12), (1, -12), (-1, -12), (12, 1), (-12, 1), (12, -1), (-12, -1)]
     """
-    if p < 1 or k < 1 or n < 1:
-        raise ValueError("Expected: n > 0 and k >= 1 and p >= 1")
+    n, p, k = [as_int(i) for i in (n, p, k)]
+
+    if n < 0:
+        if p % 2:
+            for t in power_representation(-n, p, k, zeros):
+                yield tuple(-i for i in t)
+        return
+
+    if p < 1 or k < 1:
+        raise ValueError(filldedent('''
+    Expecting positive integers for `(p, k)`, but got `(%s, %s)`'''
+    % (p, k)))
+
+    if n == 0:
+        if zeros:
+            yield (0,)*k
+        return
 
     if k == 1:
-        if perfect_power(n):
-            yield (perfect_power(n)[0],)
+        if p == 1:
+            yield (n,)
         else:
-            yield tuple()
-
-    elif p == 1:
-        for t in partition(n, k, zeros):
+            be = perfect_power(n)
+            if be:
+                b, e = be
+                d, r = divmod(e, p)
+                if not r:
+                    yield (b**d,)
+        return
+
+    if p == 1:
+        for t in partition(n, k, zeros=zeros):
+            yield t
+        return
+
+    if p == 2:
+        feasible = _can_do_sum_of_squares(n, k)
+        if not feasible:
+            return
+        if not zeros and n > 33 and k >= 5 and k <= n and n - k in (
+                13, 10, 7, 5, 4, 2, 1):
+            '''Todd G. Will, "When Is n^2 a Sum of k Squares?", [online].
+                Available: https://www.maa.org/sites/default/files/Will-MMz-201037918.pdf'''
+            return
+        if feasible is 1:  # it's prime and k == 2
+            yield prime_as_sum_of_two_squares(n)
+            return
+
+    if k == 2 and p > 2:
+        be = perfect_power(n)
+        if be and be[1] % p == 0:
+            return  # Fermat: a**n + b**n = c**n has no solution for n > 2
+
+    if n >= k:
+        a = integer_nthroot(n - (k - 1), p)[0]
+        for t in pow_rep_recursive(a, k, n, [], p):
             yield t
 
-    else:
-        l = []
+    if zeros:
         a = integer_nthroot(n, p)[0]
+        for i in range(1, k):
+            for t in pow_rep_recursive(a, i, n, [], p):
+                yield (0,) * (k - i) + t
 
-        for t in pow_rep_recursive(a, k, n, [], p):
-                yield t
 
-        if zeros:
-            for i in range(2, k):
-                for t in pow_rep_recursive(a, i, n, [], p):
-                    yield t + (0,) * (k - i)
+sum_of_powers = power_representation
 
 
 def pow_rep_recursive(n_i, k, n_remaining, terms, p):
 
     if k == 0 and n_remaining == 0:
-        yield tuple(terms)
+        yield _sorted_tuple(*terms)
     else:
         if n_i >= 1 and k > 0 and n_remaining >= 0:
             if n_i**p <= n_remaining:
                 for t in pow_rep_recursive(n_i, k - 1, n_remaining - n_i**p, terms + [n_i], p):
-                    yield t
+                    yield _sorted_tuple(*t)
 
             for t in pow_rep_recursive(n_i - 1, k, n_remaining, terms, p):
-                yield t
+                yield _sorted_tuple(*t)
+
+
+def sum_of_squares(n, k, zeros=False):
+    """Return a generator that yields the k-tuples of nonnegative
+    values, the squares of which sum to n. If zeros is False (default)
+    then the solution will not contain zeros. The nonnegative
+    elements of a tuple are sorted.
+
+    * If k == 1 and n is square, (n,) is returned.
+
+    * If k == 2 then n can only be written as a sum of squares if
+      every prime in the factorization of n that has the form
+      4*k + 3 has an even multiplicity. If n is prime then
+      it can only be written as a sum of two squares if it is
+      in the form 4*k + 1.
+
+    * if k == 3 then n can be written as a sum of squares if it does
+      not have the form 4**m*(8*k + 7).
+
+    * all integers can be written as the sum of 4 squares.
+
+    * if k > 4 then n can be partitioned and each partition can
+      be written as a sum of 4 squares; if n is not evenly divisible
+      by 4 then n can be written as a sum of squares only if the
+      an additional partition can be written as as sum of squares.
+      For example, if k = 6 then n is partitioned into two parts,
+      the first being written as a sum of 4 squares and the second
+      being written as a sum of 2 squares -- which can only be
+      done if the contition above for k = 2 can be met, so this will
+      automatically reject certain partitions of n.
+
+    Examples
+    ========
+
+    >>> from sympy.solvers.diophantine import sum_of_squares
+    >>> list(sum_of_squares(25, 2))
+    [(3, 4)]
+    >>> list(sum_of_squares(25, 2, True))
+    [(3, 4), (0, 5)]
+    >>> list(sum_of_squares(25, 4))
+    [(1, 2, 2, 4)]
+
+    See Also
+    ========
+    sympy.utilities.iterables.signed_permutations
+    """
+    for t in power_representation(n, 2, k, zeros):
+        yield t
+
+
+def _can_do_sum_of_squares(n, k):
+    """Return True if n can be written as the sum of k squares,
+    False if it can't, or 1 if k == 2 and n is prime (in which
+    case it *can* be written as a sum of two squares). A False
+    is returned only if it can't be written as k-squares, even
+    if 0s are allowed.
+    """
+    if k < 1:
+        return False
+    if n < 0:
+        return False
+    if n == 0:
+        return True
+    if k == 1:
+        return is_square(n)
+    if k == 2:
+        if n in (1, 2):
+            return True
+        if isprime(n):
+            if n % 4 == 1:
+                return 1  # signal that it was prime
+            return False
+        else:
+            f = factorint(n)
+            for p, m in f.items():
+                # we can proceed iff no prime factor in the form 4*k + 3
+                # has an odd multiplicity
+                if (p % 4 == 3) and m % 2:
+                    return False
+            return True
+    if k == 3:
+        if (n//4**multiplicity(4, n)) % 8 == 7:
+            return False
+    # every number can be written as a sum of 4 squares; for k > 4 partitions
+    # can be 0
+    return True
diff --git a/sympy/tensor/array/arrayop.py b/sympy/tensor/array/arrayop.py
index ce53884225..b9737f924c 100644
--- a/sympy/tensor/array/arrayop.py
+++ b/sympy/tensor/array/arrayop.py
@@ -123,39 +123,16 @@ def tensorcontraction(array, *contraction_axes):
         cum_shape[rank - i - 1] = _cumul
         _cumul *= int(array.shape[rank - i - 1])
 
-    # DEFINITION: by absolute position it is meant the position along the one
-    # dimensional array containing all the tensor components.
-
-    # Possible future work on this module: move computation of absolute
-    # positions to a class method.
-
-    # Determine absolute positions of the uncontracted indices:
     remaining_indices = [[cum_shape[i]*j for j in range(array.shape[i])]
                          for i in range(rank) if i not in taken_dims]
 
-    # Determine absolute positions of the contracted indices:
-    summed_deltas = []
-    for axes_group in contraction_axes:
-        lidx = []
-        for js in range(array.shape[axes_group[0]]):
-            lidx.append(sum([cum_shape[ig] * js for ig in axes_group]))
-        summed_deltas.append(lidx)
-
-    # Compute the contracted array:
-    #
-    # 1. external for loops on all uncontracted indices.
-    #    Uncontracted indices are determined by the combinatorial product of
-    #    the absolute positions of the remaining indices.
-    # 2. internal loop on all contracted indices.
-    #    It sum the values of the absolute contracted index and the absolute
-    #    uncontracted index for the external loop.
     contracted_array = []
     for icontrib in itertools.product(*remaining_indices):
-        index_base_position = sum(icontrib)
+        i = sum(icontrib)
         isum = S.Zero
-        for sum_to_index in itertools.product(*summed_deltas):
-            isum += array[index_base_position + sum(sum_to_index)]
-
+        for axes_group in contraction_axes:
+            for js in range(array.shape[axes_group[0]]):
+                isum += array[i + sum([cum_shape[ig]*js for ig in axes_group])]
         contracted_array.append(isum)
 
     if len(remaining_indices) == 0:
diff --git a/sympy/utilities/iterables.py b/sympy/utilities/iterables.py
index 1ee65bf5d8..98aecd3f98 100644
--- a/sympy/utilities/iterables.py
+++ b/sympy/utilities/iterables.py
@@ -1309,7 +1309,7 @@ def multiset_partitions(multiset, m=None):
 
 
 def partitions(n, m=None, k=None, size=False):
-    """Generate all partitions of integer n (>= 0).
+    """Generate all partitions of positive integer, n.
 
     Parameters
     ==========
@@ -1342,7 +1342,8 @@ def partitions(n, m=None, k=None, size=False):
     {1: 6}
 
     The maximum number of parts in the partition (the sum of the values in
-    the returned dict) are limited with m:
+    the returned dict) are limited with m (default value, None, gives
+    partitions from 1 through n):
 
     >>> for p in partitions(6, m=2):  # doctest: +SKIP
     ...     print(p)
@@ -1377,20 +1378,34 @@ def partitions(n, m=None, k=None, size=False):
     sympy.combinatorics.partitions.IntegerPartition
 
     """
-    if n < 0:
-        raise ValueError("n must be >= 0")
-    if m == 0:
-        raise ValueError("m must be > 0")
-    m = min(m or n, n)
-    if m < 1:
-        raise ValueError("maximum numbers in partition, m, must be > 0")
-    k = min(k or n, n)
-    if k < 1:
-        raise ValueError("maximum value in partition, k, must be > 0")
+    if (
+            n <= 0 or
+            m is not None and m < 1 or
+            k is not None and k < 1 or
+            m and k and m*k < n):
+        # the empty set is the only way to handle these inputs
+        # and returning {} to represent it is consistent with
+        # the counting convention, e.g. nT(0) == 1.
+        if size:
+            yield 0, {}
+        else:
+            yield {}
+        return
+
+    if m is None:
+        m = n
+    else:
+        m = min(m, n)
 
-    if m*k < n:
+    if n == 0:
+        if size:
+            yield 1, {0: 1}
+        else:
+            yield {0: 1}
         return
 
+    k = min(k or n, n)
+
     n, m, k = as_int(n), as_int(m), as_int(k)
     q, r = divmod(n, k)
     ms = {k: q}
@@ -1445,6 +1460,140 @@ def partitions(n, m=None, k=None, size=False):
             yield ms
 
 
+def ordered_partitions(n, m=None, sort=True):
+    """Generates ordered partitions of integer ``n``.
+
+    Parameters
+    ==========
+
+    ``m`` : integer (default gives partitions of all sizes) else only
+        those with size m. In addition, if ``m`` is not None then
+        partitions are generated *in place* (see examples).
+    ``sort`` : bool (default True) controls whether partitions are
+        returned in sorted order when ``m`` is not None; when False,
+        the partitions are returned as fast as possible with elements
+        sorted, but when m|n the partitions will not be in
+        ascending lexicographical order.
+
+    Examples
+    ========
+
+    >>> from sympy.utilities.iterables import ordered_partitions
+
+    All partitions of 5 in ascending lexicographical:
+
+    >>> for p in ordered_partitions(5):
+    ...     print(p)
+    [1, 1, 1, 1, 1]
+    [1, 1, 1, 2]
+    [1, 1, 3]
+    [1, 2, 2]
+    [1, 4]
+    [2, 3]
+    [5]
+
+    Only partitions of 5 with two parts:
+
+    >>> for p in ordered_partitions(5, 2):
+    ...     print(p)
+    [1, 4]
+    [2, 3]
+
+    When ``m`` is given, a given list objects will be used more than
+    once for speed reasons so you will not see the correct partitions
+    unless you make a copy of each as it is generated:
+
+    >>> [p for p in ordered_partitions(7, 3)]
+    [[1, 1, 1], [1, 1, 1], [1, 1, 1], [2, 2, 2]]
+    >>> [list(p) for p in ordered_partitions(7, 3)]
+    [[1, 1, 5], [1, 2, 4], [1, 3, 3], [2, 2, 3]]
+
+    When ``n`` is a multiple of ``m``, the elements are still sorted
+    but the partitions themselves will be *unordered* if sort is False;
+    the default is to return them in ascending lexicographical order.
+
+    >>> for p in ordered_partitions(6, 2):
+    ...     print(p)
+    [1, 5]
+    [2, 4]
+    [3, 3]
+
+    But if speed is more important than ordering, sort can be set to
+    False:
+
+    >>> for p in ordered_partitions(6, 2, sort=False):
+    ...     print(p)
+    [1, 5]
+    [3, 3]
+    [2, 4]
+
+    References
+    ==========
+
+    .. [1] Generating Integer Partitions, [online],
+        Available: http://jeromekelleher.net/generating-integer-partitions.html
+    .. [2] Jerome Kelleher and Barry O'Sullivan, "Generating All
+        Partitions: A Comparison Of Two Encodings", [online],
+        Available: http://arxiv.org/pdf/0909.2331v2.pdf
+    """
+    if n < 1 or m is not None and m < 1:
+        # the empty set is the only way to handle these inputs
+        # and returning {} to represent it is consistent with
+        # the counting convention, e.g. nT(0) == 1.
+        yield []
+        return
+
+    if m is None:
+        # The list `a`'s leading elements contain the partition in which
+        # y is the biggest element and x is either the same as y or the
+        # 2nd largest element; v and w are adjacent element indices
+        # to which x and y are being assigned, respectively.
+        a = [1]*n
+        y = -1
+        v = n
+        while v > 0:
+            v -= 1
+            x = a[v] + 1
+            while y >= 2 * x:
+                a[v] = x
+                y -= x
+                v += 1
+            w = v + 1
+            while x <= y:
+                a[v] = x
+                a[w] = y
+                yield a[:w + 1]
+                x += 1
+                y -= 1
+            a[v] = x + y
+            y = a[v] - 1
+            yield a[:w]
+    elif m == 1:
+        yield [n]
+    elif n == m:
+        yield [1]*n
+    else:
+        # recursively generate partitions of size m
+        for b in range(1, n//m + 1):
+            a = [b]*m
+            x = n - b*m
+            if not x:
+                if sort:
+                    yield a
+            elif not sort and x <= m:
+                for ax in ordered_partitions(x, sort=False):
+                    mi = len(ax)
+                    a[-mi:] = [i + b for i in ax]
+                    yield a
+                    a[-mi:] = [b]*mi
+            else:
+                for mi in range(1, m):
+                    for ax in ordered_partitions(x, mi, sort=True):
+                        a[-mi:] = [i + b for i in ax]
+                        yield a
+                        a[-mi:] = [b]*mi
+
+
 def binary_partitions(n):
     """
     Generates the binary partition of n.
@@ -2084,3 +2233,38 @@ def partition(lista, bins):
     else:
         raise ValueError(
             'ordered must be one of 00, 01, 10 or 11, not %s' % ordered)
+
+
+def permute_signs(t):
+    """Return iterator in which the signs of non-zero elements
+    of t are permuted.
+
+    Examples
+    ========
+
+    >>> from sympy.utilities.iterables import permute_signs
+    >>> list(permute_signs((0, 1, 2)))
+    [(0, 1, 2), (0, -1, 2), (0, 1, -2), (0, -1, -2)]
+    """
+    for signs in cartes(*[(1, -1)]*(len(t) - t.count(0))):
+        signs = list(signs)
+        yield type(t)([i*signs.pop() if i else i for i in t])
+
+
+def signed_permutations(t):
+    """Return iterator in which the signs of non-zero elements
+    of t and the order of the elements are permuted.
+
+    Examples
+    ========
+
+    >>> from sympy.utilities.iterables import signed_permutations
+    >>> list(signed_permutations((0, 1, 2)))
+    [(0, 1, 2), (0, -1, 2), (0, 1, -2), (0, -1, -2), (0, 2, 1),
+    (0, -2, 1), (0, 2, -1), (0, -2, -1), (1, 0, 2), (-1, 0, 2),
+    (1, 0, -2), (-1, 0, -2), (1, 2, 0), (-1, 2, 0), (1, -2, 0),
+    (-1, -2, 0), (2, 0, 1), (-2, 0, 1), (2, 0, -1), (-2, 0, -1),
+    (2, 1, 0), (-2, 1, 0), (2, -1, 0), (-2, -1, 0)]
+    """
+    return (type(t)(i) for j in permutations(t)
+        for i in permute_signs(j))


Test patch:
diff --git a/sympy/core/tests/test_expr.py b/sympy/core/tests/test_expr.py
index 2f3247f44a..2ad2d33f5b 100644
--- a/sympy/core/tests/test_expr.py
+++ b/sympy/core/tests/test_expr.py
@@ -1305,6 +1305,7 @@ def test_as_coeff_Add():
     assert (Integer(3) + x).as_coeff_Add() == (Integer(3), x)
     assert (Rational(3, 4) + x).as_coeff_Add() == (Rational(3, 4), x)
     assert (Float(5.0) + x).as_coeff_Add() == (Float(5.0), x)
+    assert (Float(5.0) + x).as_coeff_Add(rational=True) == (0, Float(5.0) + x)
 
     assert (Integer(3) + x + y).as_coeff_Add() == (Integer(3), x + y)
     assert (Rational(3, 4) + x + y).as_coeff_Add() == (Rational(3, 4), x + y)
diff --git a/sympy/printing/tests/test_llvmjit.py b/sympy/printing/tests/test_llvmjit.py
index 5bb9134855..59d9c63259 100644
--- a/sympy/printing/tests/test_llvmjit.py
+++ b/sympy/printing/tests/test_llvmjit.py
@@ -1,6 +1,5 @@
 
 from sympy.external import import_module
-from sympy.utilities.pytest import raises
 import ctypes
 
 
@@ -134,92 +133,11 @@ def test_callback_alt_two():
     assert isclose(jit_res, res)
 
 
-def test_multiple_statements():
-    # Match return from CSE
-    e = [[(b, 4.0*a)], [b + 5]]
-    f = g.llvm_callable([a], e)
-    b_val = e[0][0][1].subs({a: 1.5})
-    res = float(e[1][0].subs({b: b_val}).evalf())
-    jit_res = f(1.5)
-    assert isclose(jit_res, res)
-
-    f_callback = g.llvm_callable([a], e, callback_type='scipy.integrate.test')
-    m = ctypes.c_int(1)
-    array_type = ctypes.c_double * 1
-    array = array_type(1.5)
-    jit_callback_res = f_callback(m, array)
-    assert isclose(jit_callback_res, res)
-
-
-def test_cse():
-    e = a*a + b*b + sympy.exp(-a*a - b*b)
-    e2 = sympy.cse(e)
-    f = g.llvm_callable([a, b], e2)
-    res = float(e.subs({a: 2.3, b: 0.1}).evalf())
-    jit_res = f(2.3, 0.1)
-
-    assert isclose(jit_res, res)
-
-
-def eval_cse(e, sub_dict):
-    tmp_dict = dict()
-    for tmp_name, tmp_expr in e[0]:
-        e2 = tmp_expr.subs(sub_dict)
-        e3 = e2.subs(tmp_dict)
-        tmp_dict[tmp_name] = e3
-    return [e.subs(sub_dict).subs(tmp_dict) for e in e[1]]
-
-
-def test_cse_multiple():
-    e1 = a*a
-    e2 = a*a + b*b
-    e3 = sympy.cse([e1, e2])
-
-    raises(NotImplementedError,
-           lambda: g.llvm_callable([a, b], e3, callback_type='scipy.integrate'))
-
-    f = g.llvm_callable([a, b], e3)
-    jit_res = f(0.1, 1.5)
-    assert len(jit_res) == 2
-    res = eval_cse(e3, {a: 0.1, b: 1.5})
-    assert isclose(res[0], jit_res[0])
-    assert isclose(res[1], jit_res[1])
-
-
-def test_callback_cubature_multiple():
-    e1 = a*a
-    e2 = a*a + b*b
-    e3 = sympy.cse([e1, e2, 4*e2])
-    f = g.llvm_callable([a, b], e3, callback_type='cubature')
-
-    # Number of input variables
-    ndim = 2
-    # Number of output expression values
-    outdim = 3
-
-    m = ctypes.c_int(ndim)
-    fdim = ctypes.c_int(outdim)
-    array_type = ctypes.c_double * ndim
-    out_array_type = ctypes.c_double * outdim
-    inp = {a: 0.2, b: 1.5}
-    array = array_type(inp[a], inp[b])
-    out_array = out_array_type()
-    jit_ret = f(m, array, None, fdim, out_array)
-
-    assert jit_ret == 0
-
-    res = eval_cse(e3, inp)
-
-    assert isclose(out_array[0], res[0])
-    assert isclose(out_array[1], res[1])
-    assert isclose(out_array[2], res[2])
-
-
-def test_symbol_not_found():
-    e = a*a + b
-    raises(LookupError, lambda: g.llvm_callable([a], e))
-
-
 def test_bad_callback():
     e = a
-    raises(ValueError, lambda: g.llvm_callable([a], e, callback_type='bad_callback'))
+    try:
+        g.llvm_callable([a], e, callback_type='bad_callback')
+    except ValueError:
+        pass
+    else:
+        assert False, "Should raise exception with unknown callback"
diff --git a/sympy/simplify/tests/test_simplify.py b/sympy/simplify/tests/test_simplify.py
index 63f3a209c6..ddb48b6988 100644
--- a/sympy/simplify/tests/test_simplify.py
+++ b/sympy/simplify/tests/test_simplify.py
@@ -644,3 +644,7 @@ def test_clear_coefficients():
     assert clear_coefficients(4*y*(6*x + 3) - 2, x) == (y*(2*x + 1), x/12 + S(1)/6)
     assert clear_coefficients(sqrt(2) - 2) == (sqrt(2), 2)
     assert clear_coefficients(4*sqrt(2) - 2) == (sqrt(2), S.Half)
+    assert clear_coefficients(S(3), x) == (0, x - 3)
+    assert clear_coefficients(S.Infinity, x) == (S.Infinity, x)
+    assert clear_coefficients(-S.Pi, x) == (S.Pi, -x)
+    assert clear_coefficients(2 - S.Pi/3, x) == (pi, -3*x + 6)
diff --git a/sympy/solvers/tests/test_diophantine.py b/sympy/solvers/tests/test_diophantine.py
index 3617763686..95ed44ad4b 100644
--- a/sympy/solvers/tests/test_diophantine.py
+++ b/sympy/solvers/tests/test_diophantine.py
@@ -1,27 +1,65 @@
+from sympy import (Add, factor_list, igcd, Matrix, Mul, S, simplify,
+    Symbol, symbols, Eq, pi, factorint, oo, powsimp)
+from sympy.core.function import _mexpand
+from sympy.core.compatibility import range
+from sympy.functions.elementary.trigonometric import sin
 from sympy.solvers.diophantine import (descent, diop_bf_DN, diop_DN,
     diop_solve, diophantine, divisible, equivalent, find_DN, ldescent, length,
-    pairwise_prime, partition, power_representation,
+    reconstruct, partition, power_representation,
     prime_as_sum_of_two_squares, square_factor, sum_of_four_squares,
-    sum_of_three_squares, transformation_to_DN, transformation_to_normal)
+    sum_of_three_squares, transformation_to_DN, transformation_to_normal,
+    classify_diop, base_solution_linear, cornacchia, sqf_normal,
+    diop_ternary_quadratic_normal, _diop_ternary_quadratic_normal,
+    gaussian_reduce, holzer,diop_general_pythagorean,
+    _diop_general_sum_of_squares, _nint_or_floor, _odd, _even,
+    _remove_gcd, check_param, parametrize_ternary_quadratic,
+    diop_ternary_quadratic, diop_linear, diop_quadratic,
+    diop_general_sum_of_squares, sum_of_powers, sum_of_squares,
+    diop_general_sum_of_even_powers, _can_do_sum_of_squares)
+from sympy.utilities import default_sort_key
 
-from sympy import (Add, factor_list, igcd, Integer, Matrix, Mul, S, simplify,
-    Subs, Symbol, symbols)
+from sympy.utilities.pytest import slow, raises, XFAIL
+
+a, b, c, d, p, q, x, y, z, w, t, u, v, X, Y, Z = symbols(
+    "a, b, c, d, p, q, x, y, z, w, t, u, v, X, Y, Z", integer=True)
+t_0, t_1, t_2, t_3, t_4, t_5, t_6 = symbols("t_:7", integer=True)
+m1, m2, m3 = symbols('m1:4', integer=True)
+n1 = symbols('n1', integer=True)
 
-from sympy.core.function import _mexpand
-from sympy.core.compatibility import range
-from sympy.functions.elementary.trigonometric import sin
-from sympy.utilities.pytest import slow, raises
-from sympy.utilities import default_sort_key
 
-x, y, z, w, t, X, Y, Z = symbols("x, y, z, w, t, X, Y, Z", integer=True)
-t_0, t_1, t_2, t_3, t_4, t_5, t_6 = symbols("t_0, t_1, t_2, t_3, t_4, t_5, t_6", integer=True)
+def diop_simplify(eq):
+    return _mexpand(powsimp(_mexpand(eq)))
+
 
 def test_input_format():
     raises(TypeError, lambda: diophantine(sin(x)))
+    raises(TypeError, lambda: diophantine(3))
+    raises(TypeError, lambda: diophantine(x/pi - 3))
+
 
 def test_univariate():
-    assert diop_solve((x - 1)*(x - 2)**2) == set([(Integer(1),), (Integer(2),)])
-    assert diop_solve((x - 1)*(x - 2)) == set([(Integer(1),), (Integer(2),)])
+    assert diop_solve((x - 1)*(x - 2)**2) == set([(1,), (2,)])
+    assert diop_solve((x - 1)*(x - 2)) == set([(1,), (2,)])
+
+
+def test_classify_diop():
+    raises(TypeError, lambda: classify_diop(x**2/3 - 1))
+    raises(ValueError, lambda: classify_diop(1))
+    raises(NotImplementedError, lambda: classify_diop(w*x*y*z - 1))
+    assert classify_diop(14*x**2 + 15*x - 42) == (
+        [x], {1: -42, x: 15, x**2: 14}, 'univariate')
+    assert classify_diop(x*y + z) == (
+        [x, y, z], {x*y: 1, z: 1}, 'inhomogeneous_ternary_quadratic')
+    assert classify_diop(x*y + z + w + x**2) == (
+        [w, x, y, z], {x*y: 1, w: 1, x**2: 1, z: 1}, 'inhomogeneous_general_quadratic')
+    assert classify_diop(x*y + x*z + x**2 + 1) == (
+        [x, y, z], {x*y: 1, x*z: 1, x**2: 1, 1: 1}, 'inhomogeneous_general_quadratic')
+    assert classify_diop(x*y + z + w + 42) == (
+        [w, x, y, z], {x*y: 1, w: 1, 1: 42, z: 1}, 'inhomogeneous_general_quadratic')
+    assert classify_diop(x*y + z*w) == (
+        [w, x, y, z], {x*y: 1, w*z: 1}, 'homogeneous_general_quadratic')
+    assert classify_diop(x*y**2 + 1) == (
+        [x, y], {x*y**2: 1, 1: 1}, 'cubic_thue')
 
 
 def test_linear():
@@ -36,33 +74,37 @@ def test_linear():
     assert diop_solve(y + x - 0) == (t_0, -t_0)
     assert diop_solve(0*x - y - 5) == (-5,)
     assert diop_solve(3*y + 2*x - 5) == (3*t_0 - 5, -2*t_0 + 5)
-    assert diop_solve(2*x - 3*y - 5) == (-3*t_0 - 5, -2*t_0 - 5)
-    assert diop_solve(-2*x - 3*y - 5) == (-3*t_0 + 5, 2*t_0 - 5)
+    assert diop_solve(2*x - 3*y - 5) == (3*t_0 - 5, 2*t_0 - 5)
+    assert diop_solve(-2*x - 3*y - 5) == (3*t_0 + 5, -2*t_0 - 5)
     assert diop_solve(7*x + 5*y) == (5*t_0, -7*t_0)
     assert diop_solve(2*x + 4*y) == (2*t_0, -t_0)
     assert diop_solve(4*x + 6*y - 4) == (3*t_0 - 2, -2*t_0 + 2)
     assert diop_solve(4*x + 6*y - 3) == (None, None)
-    assert diop_solve(0*x + 3*y - 4*z + 5) == (-4*t_0 + 5, -3*t_0 + 5)
-    assert diop_solve(4*x + 3*y - 4*z + 5) == (t_0, -4*t_1 + 5, t_0 - 3*t_1 + 5)
+    assert diop_solve(0*x + 3*y - 4*z + 5) == (4*t_0 + 5, 3*t_0 + 5)
+    assert diop_solve(4*x + 3*y - 4*z + 5) == (t_0, 8*t_0 + 4*t_1 + 5, 7*t_0 + 3*t_1 + 5)
+    assert diop_solve(4*x + 3*y - 4*z + 5, None) == (0, 5, 5)
     assert diop_solve(4*x + 2*y + 8*z - 5) == (None, None, None)
-    assert diop_solve(5*x + 7*y - 2*z - 6) == (t_0, -7*t_0 - 2*t_1 + 6, -22*t_0 - 7*t_1 + 18)
+    assert diop_solve(5*x + 7*y - 2*z - 6) == (t_0, -3*t_0 + 2*t_1 + 6, -8*t_0 + 7*t_1 + 18)
     assert diop_solve(3*x - 6*y + 12*z - 9) == (2*t_0 + 3, t_0 + 2*t_1, t_1)
-    assert diop_solve(6*w + 9*x + 20*y - z) == (t_0, t_1, -t_1 - t_2, 6*t_0 - 11*t_1 - 20*t_2)
+    assert diop_solve(6*w + 9*x + 20*y - z) == (t_0, t_1, t_1 + t_2, 6*t_0 + 29*t_1 + 20*t_2)
+
+    # to ignore constant factors, use diophantine
+    raises(TypeError, lambda: diop_solve(x/2))
 
 
 def test_quadratic_simple_hyperbolic_case():
     # Simple Hyperbolic case: A = C = 0 and B != 0
     assert diop_solve(3*x*y + 34*x - 12*y + 1) == \
-        set([(-Integer(133), -Integer(11)), (Integer(5), -Integer(57))])
+        set([(-133, -11), (5, -57)])
     assert diop_solve(6*x*y + 2*x + 3*y + 1) == set([])
-    assert diop_solve(-13*x*y + 2*x - 4*y - 54) == set([(Integer(27), Integer(0))])
-    assert diop_solve(-27*x*y - 30*x - 12*y - 54) == set([(-Integer(14), -Integer(1))])
-    assert diop_solve(2*x*y + 5*x + 56*y + 7) == set([(-Integer(161), -Integer(3)),\
-        (-Integer(47),-Integer(6)), (-Integer(35), -Integer(12)), (-Integer(29), -Integer(69)),\
-        (-Integer(27), Integer(64)), (-Integer(21), Integer(7)),(-Integer(9), Integer(1)),\
-        (Integer(105), -Integer(2))])
+    assert diop_solve(-13*x*y + 2*x - 4*y - 54) == set([(27, 0)])
+    assert diop_solve(-27*x*y - 30*x - 12*y - 54) == set([(-14, -1)])
+    assert diop_solve(2*x*y + 5*x + 56*y + 7) == set([(-161, -3),\
+        (-47,-6), (-35, -12), (-29, -69),\
+        (-27, 64), (-21, 7),(-9, 1),\
+        (105, -2)])
     assert diop_solve(6*x*y + 9*x + 2*y + 3) == set([])
-    assert diop_solve(x*y + x + y + 1) == set([(-Integer(1), t), (t, -Integer(1))])
+    assert diop_solve(x*y + x + y + 1) == set([(-1, t), (t, -1)])
     assert diophantine(48*x*y)
 
 
@@ -71,12 +113,12 @@ def test_quadratic_elliptical_case():
     # Two test cases highlighted require lot of memory due to quadratic_congruence() method.
     # This above method should be replaced by Pernici's square_mod() method when his PR gets merged.
 
-    #assert diop_solve(42*x**2 + 8*x*y + 15*y**2 + 23*x + 17*y - 4915) == set([(-Integer(11), -Integer(1))])
+    #assert diop_solve(42*x**2 + 8*x*y + 15*y**2 + 23*x + 17*y - 4915) == set([(-11, -1)])
     assert diop_solve(4*x**2 + 3*y**2 + 5*x - 11*y + 12) == set([])
-    assert diop_solve(x**2 + y**2 + 2*x + 2*y + 2) == set([(-Integer(1), -Integer(1))])
-    #assert diop_solve(15*x**2 - 9*x*y + 14*y**2 - 23*x - 14*y - 4950) == set([(-Integer(15), Integer(6))])
+    assert diop_solve(x**2 + y**2 + 2*x + 2*y + 2) == set([(-1, -1)])
+    #assert diop_solve(15*x**2 - 9*x*y + 14*y**2 - 23*x - 14*y - 4950) == set([(-15, 6)])
     assert diop_solve(10*x**2 + 12*x*y + 12*y**2 - 34) == \
-        set([(Integer(1), -Integer(2)), (-Integer(1), -Integer(1)),(Integer(1), Integer(1)), (-Integer(1), Integer(2))])
+        set([(1, -2), (-1, -1),(1, 1), (-1, 2)])
 
 
 def test_quadratic_parabolic_case():
@@ -115,7 +157,10 @@ def test_quadratic_non_perfect_square():
 
 
 def test_issue_9106():
-    assert check_integrality(-48 - 2*x*(3*x - 1) + y*(3*y - 1))
+    eq = -48 - 2*x*(3*x - 1) + y*(3*y - 1)
+    v = (x, y)
+    for sol in diophantine(eq):
+        assert not diop_simplify(eq.xreplace(dict(zip(v, sol))))
 
 
 @slow
@@ -151,8 +196,8 @@ def test_DN():
 
     # When equation is x**2 + y**2 = N
     # Solutions are interchangeable
-    assert diop_DN(-1, 5) == [(2, 1)]
-    assert diop_DN(-1, 169) == [(12, 5), (0, 13)]
+    assert diop_DN(-1, 5) == [(1, 2)]
+    assert diop_DN(-1, 169) == [(5, 12), (0, 13)]
 
     # D > 0 and D is not a square
 
@@ -184,8 +229,8 @@ def test_DN():
 
     assert diop_DN(13, 27) == [(220, 61), (40, 11), (768, 213), (12, 3)]
     assert set(diop_DN(157, 12)) == \
-    set([(Integer(13), Integer(1)), (Integer(10663), Integer(851)), (Integer(579160), Integer(46222)), \
-        (Integer(483790960),Integer(38610722)), (Integer(26277068347), Integer(2097138361)), (Integer(21950079635497), Integer(1751807067011))])
+    set([(13, 1), (10663, 851), (579160, 46222), \
+        (483790960,38610722), (26277068347, 2097138361), (21950079635497, 1751807067011)])
     assert diop_DN(13, 25) == [(3245, 900)]
     assert diop_DN(192, 18) == []
     assert diop_DN(23, 13) == [(-6, 1), (6, 1)]
@@ -199,6 +244,9 @@ def test_DN():
 
     assert diop_DN(123, -23) == [(-10, 1), (10, 1)]
 
+    assert diop_DN(0, 0, t) == [(0, t)]
+    assert diop_DN(0, -1, t) == []
+
 
 def test_bf_pell():
     assert diop_bf_DN(13, -4) == [(3, 1), (-3, 1), (36, 10)]
@@ -208,15 +256,22 @@ def test_bf_pell():
     assert diop_bf_DN(89, -8) == [(9, 1), (-9, 1)]
     assert diop_bf_DN(21257, -1) == [(13913102721304, 95427381109)]
     assert diop_bf_DN(340, -4) == [(756, 41)]
+    assert diop_bf_DN(-1, 0, t) == [(0, 0)]
+    assert diop_bf_DN(0, 0, t) == [(0, t)]
+    assert diop_bf_DN(4, 0, t) == [(2*t, t), (-2*t, t)]
+    assert diop_bf_DN(3, 0, t) == [(0, 0)]
+    assert diop_bf_DN(1, -2, t) == []
 
 
 def test_length():
+    assert length(2, 1, 0) == 1
     assert length(-2, 4, 5) == 3
-    assert length(-5, 4, 17) == 4
+    assert length(-5, 4, 17) == 5
     assert length(0, 4, 13) == 6
-    assert length(-31, 8, 613) == 67
+    assert length(-31, 8, 613) == 69
     assert length(7, 13, 11) == 23
     assert length(-40, 5, 23) == 4
+    assert length(1, 6, 4) == 2
 
 
 def is_pell_transformation_ok(eq):
@@ -230,7 +285,7 @@ def is_pell_transformation_ok(eq):
     A, B = transformation_to_DN(eq)
     u = (A*Matrix([X, Y]) + B)[0]
     v = (A*Matrix([X, Y]) + B)[1]
-    simplified = _mexpand(Subs(eq, (x, y), (u, v)).doit())
+    simplified = diop_simplify(eq.subs(zip((x, y), (u, v))))
 
     coeff = dict([reversed(t.as_independent(*[X, Y])) for t in simplified.args])
 
@@ -238,12 +293,13 @@ def is_pell_transformation_ok(eq):
         if term in coeff.keys():
             return False
 
-    for term in [X**2, Y**2, Integer(1)]:
+    for term in [X**2, Y**2, 1]:
         if term not in coeff.keys():
-            coeff[term] = Integer(0)
+            coeff[term] = 0
 
     if coeff[X**2] != 0:
-        return isinstance(S(coeff[Y**2])/coeff[X**2], Integer) and isinstance(S(coeff[Integer(1)])/coeff[X**2], Integer)
+        return divisible(coeff[Y**2], coeff[X**2]) and \
+        divisible(coeff[1], coeff[X**2])
 
     return True
 
@@ -276,6 +332,7 @@ def test_ldescent():
     for a, b in u:
         w, x, y = ldescent(a, b)
         assert a*x**2 + b*y**2 == w**2
+    assert ldescent(-1, -1) is None
 
 
 def test_diop_ternary_quadratic_normal():
@@ -296,7 +353,7 @@ def test_diop_ternary_quadratic_normal():
 def is_normal_transformation_ok(eq):
     A = transformation_to_normal(eq)
     X, Y, Z = A*Matrix([x, y, z])
-    simplified = _mexpand(Subs(eq, (x, y, z), (X, Y, Z)).doit())
+    simplified = diop_simplify(eq.subs(zip((x, y, z), (X, Y, Z))))
 
     coeff = dict([reversed(t.as_independent(*[X, Y, Z])) for t in simplified.args])
     for term in [X*Y, Y*Z, X*Z]:
@@ -314,6 +371,9 @@ def test_transformation_to_normal():
     assert is_normal_transformation_ok(x**2 + 23*x*y - 34*y*z + 12*x*z)
     assert is_normal_transformation_ok(z**2 + 34*x*y - 23*y*z + x*z)
     assert is_normal_transformation_ok(x**2 + y**2 + z**2 - x*y - y*z - x*z)
+    assert is_normal_transformation_ok(x**2 + 2*y*z + 3*z**2)
+    assert is_normal_transformation_ok(x*y + 2*x*z + 3*y*z)
+    assert is_normal_transformation_ok(2*x*z + 3*y*z)
 
 
 def test_diop_ternary_quadratic():
@@ -333,15 +393,23 @@ def test_diop_ternary_quadratic():
     assert check_solutions(x**2 + 3*y**2 + z**2 - 13*x*y - 16*y*z + 12*x*z)
     assert check_solutions(x*y - 7*y*z + 13*x*z)
 
-
-def test_pairwise_prime():
-    assert pairwise_prime(6, 10, 15) == (5, 3, 2)
-    assert pairwise_prime(2, 3, 5) == (2, 3, 5)
-    assert pairwise_prime(1, 4, 7) == (1, 4, 7)
-    assert pairwise_prime(4, 6, 5) == (1, 6, 5)
-    assert pairwise_prime(6, 10, -15) == (5, 3, -2)
-    assert pairwise_prime(-6, -10, -15) == (-5, -3, -2)
-    assert pairwise_prime(4, -6, -5) == (1, -6, -5)
+    assert diop_ternary_quadratic_normal(x**2 + y**2 + z**2) == (None, None, None)
+    assert diop_ternary_quadratic_normal(x**2 + y**2) is None
+    raises(ValueError, lambda:
+        _diop_ternary_quadratic_normal((x, y, z),
+        {x*y: 1, x**2: 2, y**2: 3, z**2: 0}))
+    eq = -2*x*y - 6*x*z + 7*y**2 - 3*y*z + 4*z**2
+    assert diop_ternary_quadratic(eq) == (7, 2, 0)
+    assert diop_ternary_quadratic_normal(4*x**2 + 5*y**2 - z**2) == \
+        (1, 0, 2)
+    assert diop_ternary_quadratic(x*y + 2*y*z) == \
+        (-2, 0, n1)
+    eq = -5*x*y - 8*x*z - 3*y*z + 8*z**2
+    assert parametrize_ternary_quadratic(eq) == \
+        (64*p**2 - 24*p*q, -64*p*q + 64*q**2, 40*p*q)
+    # this cannot be tested with diophantine because it will
+    # factor into a product
+    assert diop_solve(x*y + 2*y*z) == (-4*p*q, -2*n1*p**2 + 2*p**2, 2*p*q)
 
 
 def test_square_factor():
@@ -355,6 +423,7 @@ def test_square_factor():
     assert square_factor(52) == 2
     assert square_factor(49) == 7
     assert square_factor(392) == 14
+    assert square_factor(factorint(-12)) == 2
 
 
 def test_parametrize_ternary_quadratic():
@@ -385,6 +454,13 @@ def test_descent():
     for a, b in u:
         w, x, y = descent(a, b)
         assert a*x**2 + b*y**2 == w**2
+    # the docstring warns against bad input, so these are expected results
+    # - can't both be negative
+    raises(TypeError, lambda: descent(-1, -3))
+    # A can't be zero unless B != 1
+    raises(ZeroDivisionError, lambda: descent(0, 3))
+    # supposed to be square-free
+    raises(TypeError, lambda: descent(4, 3))
 
 
 def test_diophantine():
@@ -405,6 +481,54 @@ def test_diophantine():
     assert check_solutions(y**2 - 7*x*y + 4*y*z)
     assert check_solutions(x**2 - 2*x + 1)
 
+    assert diophantine(x - y) == diophantine(Eq(x, y))
+    assert diophantine(3*x*pi - 2*y*pi) == set([(2*t_0, 3*t_0)])
+    assert diophantine(x**2 + y**2 + z**2 - 14) == set([(1, 2, 3)])
+    assert diophantine(x**2 + 15*x/14 - 3) == set()
+    # test issue 11049
+    eq = 92*x**2 - 99*y**2 - z**2
+    coeff = eq.as_coefficients_dict()
+    assert _diop_ternary_quadratic_normal((x, y, z), coeff) == \
+        (9, 7, 51)
+    assert diophantine(eq) == set([(
+        891*p**2 + 9*q**2, -693*p**2 - 102*p*q + 7*q**2,
+        5049*p**2 - 1386*p*q - 51*q**2)])
+    eq = 2*x**2 + 2*y**2 - z**2
+    coeff = eq.as_coefficients_dict()
+    assert _diop_ternary_quadratic_normal((x, y, z), coeff) == \
+        (1, 1, 2)
+    assert diophantine(eq) == set([(
+        2*p**2 - q**2, -2*p**2 + 4*p*q - q**2,
+        4*p**2 - 4*p*q + 2*q**2)])
+    eq = 411*x**2+57*y**2-221*z**2
+    coeff = eq.as_coefficients_dict()
+    assert _diop_ternary_quadratic_normal((x, y, z), coeff) == \
+        (2021, 2645, 3066)
+    assert diophantine(eq) == \
+        set([(115197*p**2 - 446641*q**2, -150765*p**2 + 1355172*p*q -
+        584545*q**2, 174762*p**2 - 301530*p*q + 677586*q**2)])
+    eq = 573*x**2+267*y**2-984*z**2
+    coeff = eq.as_coefficients_dict()
+    assert _diop_ternary_quadratic_normal((x, y, z), coeff) == \
+        (49, 233, 127)
+    assert diophantine(eq) == \
+        set([(4361*p**2 - 16072*q**2, -20737*p**2 + 83312*p*q - 76424*q**2,
+        11303*p**2 - 41474*p*q + 41656*q**2)])
+    # this produces factors during reconstruction
+    eq = x**2 + 3*y**2 - 12*z**2
+    coeff = eq.as_coefficients_dict()
+    assert _diop_ternary_quadratic_normal((x, y, z), coeff) == \
+        (0, 2, 1)
+    assert diophantine(eq) == \
+        set([(24*p*q, 2*p**2 - 24*q**2, p**2 + 12*q**2)])
+    # solvers have not been written for every type
+    raises(NotImplementedError, lambda: diophantine(x*y**2 + 1))
+
+    # rational expressions
+    assert diophantine(1/x) == set()
+    assert diophantine(1/x + 1/y - S.Half)
+    set([(6, 3), (-2, 1), (4, 4), (1, -2), (3, 6)])
+
 
 def test_general_pythagorean():
     from sympy.abc import a, b, c, d, e
@@ -418,51 +542,56 @@ def test_general_pythagorean():
     assert check_solutions(16*a**2 - b**2 + 9*c**2 + d**2 + 25*e**2)
 
 
-def test_diop_general_sum_of_squares():
-    from sympy.abc import a, b, c, d, e, f, g, h, i
-
-    assert check_solutions(a**2 + b**2 + c**2 - 5)
-    assert check_solutions(a**2 + b**2 + c**2 - 57)
-    assert check_solutions(a**2 + b**2 + c**2 - 349560695)
-    assert check_solutions(a**2 + b**2 + c**2 + d**2 - 304)
-    assert check_solutions(a**2 + b**2 + c**2 + d**2 - 23345)
-    assert check_solutions(a**2 + b**2 + c**2 + d**2 - 23345494)
-    assert check_solutions(a**2 + b**2 + c**2 + d**2 + e**2 - 1344545)
-    assert check_solutions(a**2 + b**2 + c**2 + d**2 + e**2 + f**2 - 6933949)
-    assert check_solutions(a**2 + b**2 + c**2 + d**2 + e**2 + f**2 + g**2 - 753934)
-    assert check_solutions(a**2 + b**2 + c**2 + d**2 + e**2 + f**2 + g**2 + h**2 - 5)
-    assert check_solutions(a**2 + b**2 + c**2 + d**2 + e**2 + f**2 + g**2 + h**2 + i**2 - 693940)
-
-
-def test_partition():
-    tests = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
-
-    for test in tests:
-        f = partition(test)
-        while True:
-            try:
-                l = next(f)
-            except StopIteration:
-                break
-
-    tests_k = [8, 10]
-
-    for test in tests_k:
-        for k in range(8):
-            f = partition(test, k)
-
-            while True:
-                try:
-                    l = next(f)
-                    assert len(l) == k
-                except StopIteration:
-                    break
+def test_diop_general_sum_of_squares_quick():
+    for i in range(3, 10):
+        assert check_solutions(sum(i**2 for i in symbols(':%i' % i)) - i)
+    raises(ValueError, lambda: _diop_general_sum_of_squares((x, y), 2))
+    assert _diop_general_sum_of_squares((x, y, z), -2) == set()
+    eq = x**2 + y**2 + z**2 - (1 + 4 + 9)
+    assert diop_general_sum_of_squares(eq) == \
+        set([(1, 2, 3)])
+    eq = u**2 + v**2 + x**2 + y**2 + z**2 - 1313
+    assert len(diop_general_sum_of_squares(eq, 3)) == 3
+    # issue 11016
+    var = symbols(':5') + (symbols('6', negative=True),)
+    eq = Add(*[i**2 for i in var]) - 112
+    assert diophantine(eq) == set(
+        [(0, 1, 1, 5, 6, -7), (1, 1, 1, 3, 6, -8), (2, 3, 3, 4,
+        5, -7), (0, 1, 1, 1, 3, -10), (0, 0, 4, 4, 4, -8), (1, 2, 3,
+        3, 5, -8), (0, 1, 2, 3, 7, -7), (2, 2, 4, 4, 6, -6), (1, 1,
+        3, 4, 6, -7), (0, 2, 3, 3, 3, -9), (0, 0, 2, 2, 2, -10), (1,
+        1, 2, 3, 4, -9), (0, 1, 1, 2, 5, -9), (0, 0, 2, 6, 6, -6),
+        (1, 3, 4, 5, 5, -6), (0, 2, 2, 2, 6, -8), (0, 3, 3, 3, 6,
+        -7), (0, 2, 3, 5, 5, -7), (0, 1, 5, 5, 5, -6)])
+    # handle negated squares with signsimp
+    assert diophantine(12 - x**2 - y**2 - z**2) == set([(2, 2, 2)])
+    # diophantine handles simplification, so classify_diop should
+    # not have to look for additional patterns that are removed
+    # by diophantine
+    eq = a**2 + b**2 + c**2 + d**2 - 4
+    raises(NotImplementedError, lambda: classify_diop(-eq))
+
+
+def test_diop_partition():
+    for n in [8, 10]:
+        for k in range(1, 8):
+            for p in partition(n, k):
+                assert len(p) == k
+    assert [p for p in partition(3, 5)] == []
+    assert [list(p) for p in partition(3, 5, 1)] == [
+        [0, 0, 0, 0, 3], [0, 0, 0, 1, 2], [0, 0, 1, 1, 1]]
+    assert list(partition(0)) == [()]
+    assert list(partition(1, 0)) == [()]
+    assert [list(i) for i in partition(3)] == [[1, 1, 1], [1, 2], [3]]
 
 
 def test_prime_as_sum_of_two_squares():
     for i in [5, 13, 17, 29, 37, 41, 2341, 3557, 34841, 64601]:
         a, b = prime_as_sum_of_two_squares(i)
         assert a**2 + b**2 == i
+    assert prime_as_sum_of_two_squares(7) is None
+    ans = prime_as_sum_of_two_squares(800029)
+    assert ans == (450, 773) and type(ans[0]) is int
 
 
 def test_sum_of_three_squares():
@@ -471,17 +600,25 @@ def test_sum_of_three_squares():
         a, b, c = sum_of_three_squares(i)
         assert a**2 + b**2 + c**2 == i
 
-    assert sum_of_three_squares(7) == (None, None, None)
-    assert sum_of_three_squares((4**5)*15) == (None, None, None)
+    assert sum_of_three_squares(7) is None
+    assert sum_of_three_squares((4**5)*15) is None
+    assert sum_of_three_squares(25) == (5, 0, 0)
+    assert sum_of_three_squares(4) == (0, 0, 2)
 
 
 def test_sum_of_four_squares():
     from random import randint
 
-    for i in range(10):
-        n = randint(1, 100000000000000)
-        a, b, c, d = sum_of_four_squares(n)
-        assert a**2 + b**2 + c**2 + d**2 == n
+    # this should never fail
+    n = randint(1, 100000000000000)
+    assert sum(i**2 for i in sum_of_four_squares(n)) == n
+
+    assert sum_of_four_squares(0) == (0, 0, 0, 0)
+    assert sum_of_four_squares(14) == (0, 1, 2, 3)
+    assert sum_of_four_squares(15) == (1, 1, 2, 3)
+    assert sum_of_four_squares(18) == (1, 2, 2, 3)
+    assert sum_of_four_squares(19) == (0, 1, 3, 3)
+    assert sum_of_four_squares(48) == (0, 4, 4, 4)
 
 
 def test_power_representation():
@@ -505,6 +642,29 @@ def test_power_representation():
             except StopIteration:
                 break
 
+    assert list(power_representation(20, 2, 4, True)) == \
+        [(1, 1, 3, 3), (0, 0, 2, 4)]
+    raises(ValueError, lambda: list(power_representation(1.2, 2, 2)))
+    raises(ValueError, lambda: list(power_representation(2, 0, 2)))
+    raises(ValueError, lambda: list(power_representation(2, 2, 0)))
+    assert list(power_representation(-1, 2, 2)) == []
+    assert list(power_representation(1, 1, 1)) == [(1,)]
+    assert list(power_representation(3, 2, 1)) == []
+    assert list(power_representation(4, 2, 1)) == [(2,)]
+    assert list(power_representation(3**4, 4, 6, zeros=True)) == \
+        [(1, 2, 2, 2, 2, 2), (0, 0, 0, 0, 0, 3)]
+    assert list(power_representation(3**4, 4, 5, zeros=False)) == []
+    assert list(power_representation(-2, 3, 2)) == [(-1, -1)]
+    assert list(power_representation(-2, 4, 2)) == []
+    assert list(power_representation(0, 3, 2, True)) == [(0, 0)]
+    assert list(power_representation(0, 3, 2, False)) == []
+    # when we are dealing with squares, do feasibility checks
+    assert len(list(power_representation(4**10*(8*10 + 7), 2, 3))) == 0
+    # there will be a recursion error if these aren't recognized
+    big = 2**30
+    for i in [13, 10, 7, 5, 4, 2, 1]:
+        assert list(sum_of_powers(big, 2, big - i)) == []
+
 
 def test_assumptions():
     """
@@ -529,44 +689,182 @@ def check_solutions(eq):
     """
     s = diophantine(eq)
 
-    terms = factor_list(eq)[1]
+    factors = Mul.make_args(eq)
 
     var = list(eq.free_symbols)
     var.sort(key=default_sort_key)
 
-    okay = True
-
-    while len(s) and okay:
+    while s:
         solution = s.pop()
-
-        okay = False
-
-        for term in terms:
-            subeq = term[0]
-
-            if simplify(_mexpand(Subs(subeq, var, solution).doit())) == 0:
-                okay = True
+        for f in factors:
+            if diop_simplify(f.subs(zip(var, solution))) == 0:
                 break
+        else:
+            return False
+    return True
 
-    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
+def test_diopcoverage():
+    eq = (2*x + y + 1)**2
+    assert diop_solve(eq) == set([(t_0, -2*t_0 - 1)])
+    eq = 2*x**2 + 6*x*y + 12*x + 4*y**2 + 18*y + 18
+    assert diop_solve(eq) == set([(t_0, -t_0 - 3), (2*t_0 - 3, -t_0)])
+    assert diop_quadratic(x + y**2 - 3) == set([(-t**2 + 3, -t)])
+
+    assert diop_linear(x + y - 3) == (t_0, 3 - t_0)
+
+    assert base_solution_linear(0, 1, 2, t=None) == (0, 0)
+    ans = (3*t - 1, -2*t + 1)
+    assert base_solution_linear(4, 8, 12, t) == ans
+    assert base_solution_linear(4, 8, 12, t=None) == tuple(_.subs(t, 0) for _ in ans)
+
+    assert cornacchia(1, 1, 20) is None
+    assert cornacchia(1, 1, 5) == set([(1, 2)])
+    assert cornacchia(1, 2, 17) == set([(3, 2)])
+
+    raises(ValueError, lambda: reconstruct(4, 20, 1))
+
+    assert gaussian_reduce(4, 1, 3) == (1, 1)
+    eq = -w**2 - x**2 - y**2 + z**2
+
+    assert diop_general_pythagorean(eq) == \
+        diop_general_pythagorean(-eq) == \
+            (m1**2 + m2**2 - m3**2, 2*m1*m3,
+            2*m2*m3, m1**2 + m2**2 + m3**2)
+
+    assert check_param(S(3) + x/3, S(4) + x/2, S(2), x) == (None, None)
+    assert check_param(S(3)/2, S(4) + x, S(2), x) == (None, None)
+    assert check_param(S(4) + x, S(3)/2, S(2), x) == (None, None)
+
+    assert _nint_or_floor(16, 10) == 2
+    assert _odd(1) == (not _even(1)) == True
+    assert _odd(0) == (not _even(0)) == False
+    assert _remove_gcd(2, 4, 6) == (1, 2, 3)
+    raises(TypeError, lambda: _remove_gcd((2, 4, 6)))
+    assert sqf_normal(2 * 3**2 * 5, 2 * 5 * 11, 2 * 7**2 * 11)  == \
+        (11, 1, 5)
+
+    # it's ok if these pass some day when the solvers are implemented
+    raises(NotImplementedError, lambda: diophantine(x**2 + y**2 + x*y + 2*y*z - 12))
+    raises(NotImplementedError, lambda: diophantine(x**3 + y**2))
+
+
+def test_holzer():
+    # if the input is good, don't let it diverge in holzer()
+    # (but see test_fail_holzer below)
+    assert holzer(2, 7, 13, 4, 79, 23) == (2, 7, 13)
+
+    # None in uv condition met; solution is not Holzer reduced
+    # so this will hopefully change but is here for coverage
+    assert holzer(2, 6, 2, 1, 1, 10) == (2, 6, 2)
+
+    raises(ValueError, lambda: holzer(2, 7, 14, 4, 79, 23))
+
+
+@XFAIL
+def test_fail_holzer():
+    eq = lambda x, y, z: a*x**2 + b*y**2 - c*z**2
+    a, b, c = 4, 79, 23
+    x, y, z = xyz = 26, 1, 11
+    X, Y, Z = ans = 2, 7, 13
+    assert eq(*xyz) == 0
+    assert eq(*ans) == 0
+    assert max(a*x**2, b*y**2, c*z**2) <= a*b*c
+    assert max(a*X**2, b*Y**2, c*Z**2) <= a*b*c
+    h = holzer(x, y, z, a, b, c)
+    assert h == ans  # it would be nice to get the smaller soln
+
+
+def test_issue_9539():
+    assert diophantine(6*w + 9*y + 20*x - z) == \
+        set([(t_0, t_1, t_1 + t_2, 6*t_0 + 29*t_1 + 9*t_2)])
+
+
+def test_issue_8943():
+    assert diophantine(
+        (3*(x**2 + y**2 + z**2) - 14*(x*y + y*z + z*x))) == \
+        set([(0, 0, 0)])
+
+
+def test_diop_sum_of_even_powers():
+    eq = x**4 + y**4 + z**4 - 2673
+    assert diop_solve(eq) == set([(3, 6, 6), (2, 4, 7)])
+    assert diop_general_sum_of_even_powers(eq, 2) == set(
+        [(3, 6, 6), (2, 4, 7)])
+    raises(NotImplementedError, lambda: diop_general_sum_of_even_powers(-eq, 2))
+    neg = symbols('neg', negative=True)
+    eq = x**4 + y**4 + neg**4 - 2673
+    assert diop_general_sum_of_even_powers(eq) == set([(-2, 4, 7)])
+    assert diophantine(x**4 + y**4 + 2) == set()
+    assert diop_general_sum_of_even_powers(x**4 + y**4 - 2, limit=0) == set()
+
+
+def test_sum_of_squares_powers():
+    tru = set([
+    (0, 0, 1, 1, 11), (0, 0, 5, 7, 7), (0, 1, 3, 7, 8), (0, 1, 4, 5, 9),
+    (0, 3, 4, 7, 7), (0, 3, 5, 5, 8), (1, 1, 2, 6, 9), (1, 1, 6, 6, 7),
+    (1, 2, 3, 3, 10), (1, 3, 4, 4, 9), (1, 5, 5, 6, 6), (2, 2, 3, 5, 9),
+    (2, 3, 5, 6, 7), (3, 3, 4, 5, 8)])
+    eq = u**2 + v**2 + x**2 + y**2 + z**2 - 123
+    ans = diop_general_sum_of_squares(eq, oo)  # allow oo to be used
+    assert len(ans) == 14
+
+    raises(ValueError, lambda: list(sum_of_squares(10, -1)))
+    assert list(sum_of_squares(-10, 2)) == []
+    assert list(sum_of_squares(2, 3)) == []
+    assert list(sum_of_squares(0, 3, True)) == [(0, 0, 0)]
+    assert list(sum_of_squares(0, 3)) == []
+    assert list(sum_of_squares(4, 1)) == [(2,)]
+    assert list(sum_of_squares(5, 1)) == []
+    assert list(sum_of_squares(50, 2)) == [(1, 7), (5, 5)]
+    assert list(sum_of_squares(11, 5, True)) == [
+        (1, 1, 1, 2, 2), (0, 0, 1, 1, 3)]
+    assert list(sum_of_squares(8, 8)) == [(1, 1, 1, 1, 1, 1, 1, 1)]
+
+    assert [len(list(sum_of_squares(i, 5, True))) for i in range(30)] == [
+        1, 1, 1, 1, 2,
+        2, 1, 1, 2, 2,
+        2, 2, 2, 3, 2,
+        1, 3, 3, 3, 3,
+        4, 3, 3, 2, 2,
+        4, 4, 4, 4, 5]
+    assert [len(list(sum_of_squares(i, 5))) for i in range(30)] == [
+        0, 0, 0, 0, 0,
+        1, 0, 0, 1, 0,
+        0, 1, 0, 1, 1,
+        0, 1, 1, 0, 1,
+        2, 1, 1, 1, 1,
+        1, 1, 1, 1, 3]
+    for i in range(30):
+        s1 = set(sum_of_squares(i, 5, True))
+        assert not s1 or all(sum(j**2 for j in t) == i for t in s1)
+        s2 = set(sum_of_squares(i, 5))
+        assert all(sum(j**2 for j in t) == i for t in s2)
+
+    raises(ValueError, lambda: list(sum_of_powers(2, -1, 1)))
+    raises(ValueError, lambda: list(sum_of_powers(2, 1, -1)))
+    assert list(sum_of_powers(-2, 3, 2)) == [(-1, -1)]
+    assert list(sum_of_powers(-2, 4, 2)) == []
+    assert list(sum_of_powers(2, 1, 1)) == [(2,)]
+    assert list(sum_of_powers(2, 1, 3, True)) == [(0, 0, 2), (0, 1, 1)]
+    assert list(sum_of_powers(5, 1, 2, True)) == [(0, 5), (1, 4), (2, 3)]
+    assert list(sum_of_powers(6, 2, 2)) == []
+    assert list(sum_of_powers(3**5, 3, 1)) == []
+    assert list(sum_of_powers(3**6, 3, 1)) == [(9,)] and (9**3 == 3**6)
+    assert list(sum_of_powers(2**1000, 5, 2)) == []
+
+
+def test__can_do_sum_of_squares():
+    assert _can_do_sum_of_squares(3, -1) is False
+    assert _can_do_sum_of_squares(-3, 1) is False
+    assert _can_do_sum_of_squares(0, 1)
+    assert _can_do_sum_of_squares(4, 1)
+    assert _can_do_sum_of_squares(1, 2)
+    assert _can_do_sum_of_squares(2, 2)
+    assert _can_do_sum_of_squares(3, 2) is False
+
+
+def test_issue_9538():
+    eq = x - 3*y + 2
+    assert diophantine(eq, syms=[y,x]) == set([(t_0, 3*t_0 - 2)])
+    raises(TypeError, lambda: diophantine(eq, syms=set([y,x])))
diff --git a/sympy/tensor/array/tests/test_arrayop.py b/sympy/tensor/array/tests/test_arrayop.py
index 7992120eb3..bd63fb881f 100644
--- a/sympy/tensor/array/tests/test_arrayop.py
+++ b/sympy/tensor/array/tests/test_arrayop.py
@@ -65,22 +65,3 @@ def test_derivative_by_array():
     assert derive_by_array([[x, y], [z, t]], [x, y]) == Array([[[1, 0], [0, 0]], [[0, 1], [0, 0]]])
     assert derive_by_array([[x, y], [z, t]], [[x, y], [z, t]]) == Array([[[[1, 0], [0, 0]], [[0, 1], [0, 0]]],
                                                                          [[[0, 0], [1, 0]], [[0, 0], [0, 1]]]])
-
-
-def test_issue_emerged_while_discussing_10972():
-    ua = Array([-1,0])
-    Fa = Array([[0, 1], [-1, 0]])
-    po = tensorproduct(Fa, ua, Fa, ua)
-    assert tensorcontraction(po, (1, 2), (4, 5)) == Array([[0, 0], [0, 1]])
-
-    sa = symbols('a0:144')
-    po = Array(sa, [2, 2, 3, 3, 2, 2])
-    assert tensorcontraction(po, (0, 1), (2, 3), (4, 5)) == sa[0] + sa[108] + sa[111] + sa[124] + sa[127] + sa[140] + sa[143] + sa[16] + sa[19] + sa[3] + sa[32] + sa[35]
-    assert tensorcontraction(po, (0, 1, 4, 5), (2, 3)) == sa[0] + sa[111] + sa[127] + sa[143] + sa[16] + sa[32]
-    assert tensorcontraction(po, (0, 1), (4, 5)) == Array([[sa[0] + sa[108] + sa[111] + sa[3], sa[112] + sa[115] + sa[4] + sa[7],
-                                                             sa[11] + sa[116] + sa[119] + sa[8]], [sa[12] + sa[120] + sa[123] + sa[15],
-                                                             sa[124] + sa[127] + sa[16] + sa[19], sa[128] + sa[131] + sa[20] + sa[23]],
-                                                            [sa[132] + sa[135] + sa[24] + sa[27], sa[136] + sa[139] + sa[28] + sa[31],
-                                                             sa[140] + sa[143] + sa[32] + sa[35]]])
-    assert tensorcontraction(po, (0, 1), (2, 3)) == Array([[sa[0] + sa[108] + sa[124] + sa[140] + sa[16] + sa[32], sa[1] + sa[109] + sa[125] + sa[141] + sa[17] + sa[33]],
-                                                           [sa[110] + sa[126] + sa[142] + sa[18] + sa[2] + sa[34], sa[111] + sa[127] + sa[143] + sa[19] + sa[3] + sa[35]]])
diff --git a/sympy/utilities/tests/test_iterables.py b/sympy/utilities/tests/test_iterables.py
index 4b0a944abb..ddb406c5e6 100644
--- a/sympy/utilities/tests/test_iterables.py
+++ b/sympy/utilities/tests/test_iterables.py
@@ -15,7 +15,7 @@
     multiset_permutations, necklaces, numbered_symbols, ordered, partitions,
     permutations, postfixes, postorder_traversal, prefixes, reshape,
     rotate_left, rotate_right, runs, sift, subsets, take, topological_sort,
-    unflatten, uniq, variations)
+    unflatten, uniq, variations, ordered_partitions)
 from sympy.utilities.enumerative import (
     factoring_visitor, multiset_partitions_taocp )
 
@@ -390,6 +390,15 @@ def test():
 
 
 def test_partitions():
+    ans = [[{}], [(0, {})]]
+    for i in range(2):
+        assert list(partitions(0, size=i)) == ans[i]
+        assert list(partitions(1, 0, size=i)) == ans[i]
+        assert list(partitions(6, 2, 2, size=i)) == ans[i]
+        assert list(partitions(6, 2, None, size=i)) != ans[i]
+        assert list(partitions(6, None, 2, size=i)) != ans[i]
+        assert list(partitions(6, 2, 0, size=i)) == ans[i]
+
     assert [p.copy() for p in partitions(6, k=2)] == [
         {2: 3}, {1: 2, 2: 2}, {1: 4, 2: 1}, {1: 6}]
 
@@ -397,8 +406,6 @@ def test_partitions():
         {3: 2}, {1: 1, 2: 1, 3: 1}, {1: 3, 3: 1}, {2: 3}, {1: 2, 2: 2},
         {1: 4, 2: 1}, {1: 6}]
 
-    assert [p.copy() for p in partitions(6, k=2, m=2)] == []
-
     assert [p.copy() for p in partitions(8, k=4, m=3)] == [
         {4: 2}, {1: 1, 3: 1, 4: 1}, {2: 2, 4: 1}, {2: 1, 3: 2}] == [
         i.copy() for i in partitions(8, k=4, m=3) if all(k <= 4 for k in i)
@@ -411,7 +418,6 @@ def test_partitions():
         {1: 1, 3: 1}, {2: 2}, {1: 2, 2: 1}, {1: 4}] == [
         i.copy() for i in partitions(4) if all(k <= 3 for k in i)]
 
-    raises(ValueError, lambda: list(partitions(3, 0)))
 
     # Consistency check on output of _partitions and RGS_unrank.
     # This provides a sanity test on both routines.  Also verifies that
@@ -422,7 +428,7 @@ def test_partitions():
         i  = 0
         for m, q  in _set_partitions(n):
             assert  q == RGS_unrank(i, n)
-            i = i+1
+            i += 1
         assert i == RGS_enum(n)
 
 def test_binary_partitions():
@@ -707,3 +713,16 @@ def test__partition():
         ['b', 'e'], ['a', 'c'], ['d']]
     output = (3, [1, 0, 1, 2, 0])
     assert _partition('abcde', *output) == [['b', 'e'], ['a', 'c'], ['d']]
+
+
+def test_ordered_partitions():
+    from sympy.functions.combinatorial.numbers import nT
+    f = ordered_partitions
+    assert list(f(0, 1)) == [[]]
+    assert list(f(1, 0)) == [[]]
+    for i in range(1, 7):
+        for j in [None] + list(range(1, i)):
+            assert (
+                sum(1 for p in f(i, j, 1)) ==
+                sum(1 for p in f(i, j, 0)) ==
+                nT(i, j))