Repository: sympy/sympy
Pull request: Add class DiophantineEquationType and move classification logic
#### References to other Issues or PRs
#18493

Closes #18370.

#### Brief description of what is fixed or changed
Introduce new class `DiophantineEquationType`. The class has various attributes for properties of the equation such as `homogeneous` and `total_degree`.

Move classification logic out of `classify_diop` and into individual equation type classes.

Generate `diop_known` from the list of equation type classes.

#### Release Notes

<!-- BEGIN RELEASE NOTES -->
NO ENTRY
<!-- END RELEASE NOTES -->

Problem statement:
Confusing test cases in test_classify_diop
We have these two test cases:

```
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')
```

In what way are these equations "inhomogeneous"?

Hints from issue discussion:
~~Unfortunately it would not be the first time that dubious test results made their way into the test suite...~~

~~I'd say the tests you cited are bugs in `classify_diop`. (Or at least badly chosen default names for some classes of equations that diophantine is not currently capable of solving.) IMO it should be fixed.~~
hey how x*y+z & x*y + z + w + x**2 are homogeneous? I am considering this defination "A polynomial is homogeneous if all its terms have the same degree".
Oops, I retract my previous comment. I was probably thinking of linear equations where homogeneity reduces to the presence or absence of a constant term.
is this issue still valid?

`homogeneous` clearly is referring to the presence of a constant term here:

https://github.com/sympy/sympy/blob/9f98339b3344de9109f4d2780faf9822fb110848/sympy/solvers/diophantine.py#L497

Here the code checks that `homogeneous` is true and then immediately returns an inhomogeneous equation type:

https://github.com/sympy/sympy/blob/9f98339b3344de9109f4d2780faf9822fb110848/sympy/solvers/diophantine.py#L507-L509
In general,we look for the coefficient only in the case of linear equation which look like this `a1x1+a2x2+...+anxn=b` but  ` xy+z & xy + z + w + x**2 are polynomial` in which the only thing that matters is the degree should be same for all the term.
This is why we look for constant to be zero in linear equation 
`A homogeneous linear equation, say, 2x+3y+5z=0, is really special because it has an important property: if (x,y,z) is a solution then (λx,λy,λz) is also a solution for arbitrary constant factor λ. If you replace 0 in the right hand side of the equation by a non-zero constant, and get, say, 2x+3y+5z=1, this property is lost.`
@ethankward :
> `homogeneous` clearly is referring to the presence of a constant term here:
> 
> https://github.com/sympy/sympy/blob/9f98339b3344de9109f4d2780faf9822fb110848/sympy/solvers/diophantine.py#L497

Logically that boolean variable expresses `possibly_homogeneous` and not (definitely) `homogeneous`.

> Here the code checks that `homogeneous` is true and then immediately returns an inhomogeneous equation type:
> 
> https://github.com/sympy/sympy/blob/9f98339b3344de9109f4d2780faf9822fb110848/sympy/solvers/diophantine.py#L507-L509

As I understand it, the intersection of `set(coeff)` and `set(var)` is non-empty (hence the `if`-condition true) exactly when a "naked" element of `var` appears in the expression, i.e. a term of degree 1. But the total degree in this case is 2, so in order to have a properly homogeneous polynomial all terms would have to be of degree 2, which is not the case. Hence the equation is indeed inhomogeneous.

The naming of the `homogeneous` variable is indeed unfortunate... But the logic of that excerpt seems valid.
The ode module also makes use of both meanings of "homogeneous" and calls the uniform degree version "homogenous order" https://docs.sympy.org/1.5.1/modules/solvers/ode.html#homogeneous-order
We can also add new things such as whether polynomial is symmetric or not?
I am not sure that whether sympy contain a module or method for `Elementary Symmetric Polynomials`
If not then I think we can add that functionality too. @asmeurer 

`Elementary Symmetric Polynomials` = https://en.wikipedia.org/wiki/Elementary_symmetric_polynomial

Commit message: introduce DiophantineEquationType; simplify classify_diop; update diop_known

Patch to review:
diff --git a/sympy/solvers/diophantine/diophantine.py b/sympy/solvers/diophantine/diophantine.py
index be94c6fcd30d..efa1dfa32740 100644
--- a/sympy/solvers/diophantine/diophantine.py
+++ b/sympy/solvers/diophantine/diophantine.py
@@ -36,22 +36,6 @@
 __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"}
-
-
 class DiophantineSolutionSet(set):
     """
     Container for a set of solutions to a particular diophantine equation.
@@ -152,6 +136,223 @@ def __call__(self, *args):
         return self.subs(zip(self.parameters, args))
 
 
+class DiophantineEquationType:
+    """
+    Internal representation of a particular diophantine equation type.
+
+    Parameters
+    ==========
+
+    equation
+        The diophantine equation that is being solved.
+    free_symbols: list (optional)
+        The symbols being solved for.
+
+    Attributes
+    ==========
+
+    total_degree
+        The maximum of the degrees of all terms in the equation
+    homogeneous
+        Does the equation contain a term of degree 0
+    homogeneous_order
+        Does the equation contain any coefficient that is in the symbols being solved for
+    dimension
+        The number of symbols being solved for
+    """
+    name = None
+
+    def __init__(self, equation, free_symbols=None):
+        try:
+            if free_symbols is None:
+                free_symbols = list(equation.free_symbols)
+            assert free_symbols
+        except (AttributeError, AssertionError):
+            raise ValueError('equation should have 1 or more free symbols')
+
+        free_symbols.sort(key=default_sort_key)
+        self.free_symbols = free_symbols
+        self.equation = equation.expand(force=True)
+        self.coeff = self.equation.as_coefficients_dict()
+        if not all(_is_int(c) for c in self.coeff.values()):
+            raise TypeError("Coefficients should be Integers")
+
+        self.total_degree = Poly(self.equation).total_degree()
+        self.homogeneous = 1 not in self.coeff
+        self.homogeneous_order = not (set(self.coeff) & set(self.free_symbols))
+        self.dimension = len(self.free_symbols)
+
+    def matches(self):
+        """
+        Determine whether the given equation can be matched to the particular equation type.
+        """
+        return False
+
+
+class Univariate(DiophantineEquationType):
+    name = 'univariate'
+
+    def matches(self):
+        return self.dimension == 1
+
+
+class Linear(DiophantineEquationType):
+    name = 'linear'
+
+    def matches(self):
+        return self.total_degree == 1
+
+
+class BinaryQuadratic(DiophantineEquationType):
+    name = 'binary_quadratic'
+
+    def matches(self):
+        return self.total_degree == 2 and self.dimension == 2
+
+
+class InhomogeneousTernaryQuadratic(DiophantineEquationType):
+    name = 'inhomogeneous_ternary_quadratic'
+
+    def matches(self):
+        if not (self.total_degree == 2 and self.dimension == 3):
+            return False
+        if not self.homogeneous:
+            return False
+        return not self.homogeneous_order
+
+
+class HomogeneousTernaryQuadraticNormal(DiophantineEquationType):
+    name = 'homogeneous_ternary_quadratic_normal'
+
+    def matches(self):
+        if not (self.total_degree == 2 and self.dimension == 3):
+            return False
+        if not self.homogeneous:
+            return False
+        if not self.homogeneous_order:
+            return False
+
+        nonzero = [k for k in self.coeff if self.coeff[k]]
+        return len(nonzero) == 3 and all(i**2 in nonzero for i in self.free_symbols)
+
+
+class HomogeneousTernaryQuadratic(DiophantineEquationType):
+    name = 'homogeneous_ternary_quadratic'
+
+    def matches(self):
+        if not (self.total_degree == 2 and self.dimension == 3):
+            return False
+        if not self.homogeneous:
+            return False
+        if not self.homogeneous_order:
+            return False
+
+        nonzero = [k for k in self.coeff if self.coeff[k]]
+        return not (len(nonzero) == 3 and all(i**2 in nonzero for i in self.free_symbols))
+
+
+class InhomogeneousGeneralQuadratic(DiophantineEquationType):
+    name = 'inhomogeneous_general_quadratic'
+
+    def matches(self):
+        if not (self.total_degree == 2 and self.dimension >= 3):
+            return False
+        if not self.homogeneous_order:
+            return True
+        else:
+            # there may be Pow keys like x**2 or Mul keys like x*y
+            if any(k.is_Mul for k in self.coeff): # cross terms
+                return not self.homogeneous
+        return False
+
+
+class HomogeneousGeneralQuadratic(DiophantineEquationType):
+    name = 'homogeneous_general_quadratic'
+
+    def matches(self):
+        if not (self.total_degree == 2 and self.dimension >= 3):
+            return False
+        if not self.homogeneous_order:
+            return False
+        else:
+            # there may be Pow keys like x**2 or Mul keys like x*y
+            if any(k.is_Mul for k in self.coeff): # cross terms
+                return self.homogeneous
+        return False
+
+
+class GeneralSumOfSquares(DiophantineEquationType):
+    name = 'general_sum_of_squares'
+
+    def matches(self):
+        if not (self.total_degree == 2 and self.dimension >= 3):
+            return False
+        if not self.homogeneous_order:
+            return False
+        if any(k.is_Mul for k in self.coeff):
+            return False
+        return all(self.coeff[k] == 1 for k in self.coeff if k != 1)
+
+
+class GeneralPythagorean(DiophantineEquationType):
+    name = 'general_pythagorean'
+
+    def matches(self):
+        if not (self.total_degree == 2 and self.dimension >= 3):
+            return False
+        if not self.homogeneous_order:
+            return False
+        if any(k.is_Mul for k in self.coeff):
+            return False
+        if all(self.coeff[k] == 1 for k in self.coeff if k != 1):
+            return False
+        if not all(is_square(abs(self.coeff[k])) for k in self.coeff):
+            return False
+        # all but one has the same sign
+        # e.g. 4*x**2 + y**2 - 4*z**2
+        return abs(sum(sign(self.coeff[k]) for k in self.coeff)) == self.dimension - 2
+
+
+class CubicThue(DiophantineEquationType):
+    name = 'cubic_thue'
+
+    def matches(self):
+        return self.total_degree == 3 and self.dimension == 2
+
+
+class GeneralSumOfEvenPowers(DiophantineEquationType):
+    name = 'general_sum_of_even_powers'
+
+    def matches(self):
+        if not self.total_degree > 3:
+            return False
+        if self.total_degree % 2 != 0:
+            return False
+        if not all(k.is_Pow and k.exp == self.total_degree for k in self.coeff if k != 1):
+            return False
+        return all(self.coeff[k] == 1 for k in self.coeff if k != 1)
+
+
+# these types are known (but not necessarily handled)
+# note that order is important here (in the current solver state)
+all_diop_classes = [
+    Linear,
+    Univariate,
+    BinaryQuadratic,
+    HomogeneousTernaryQuadratic,
+    HomogeneousGeneralQuadratic,
+    HomogeneousTernaryQuadraticNormal,
+    InhomogeneousTernaryQuadratic,
+    InhomogeneousGeneralQuadratic,
+    GeneralPythagorean,
+    GeneralSumOfEvenPowers,
+    GeneralSumOfSquares,
+    CubicThue
+]
+
+diop_known = set([diop_class.name for diop_class in all_diop_classes])
+
+
 def _is_int(i):
     try:
         as_int(i)
@@ -322,12 +523,12 @@ def diophantine(eq, param=symbols("t", integer=True), syms=None,
         if permute:
             len_var = len(v)
             permute_signs_for = [
-                'general_sum_of_squares',
-                'general_sum_of_even_powers']
+                GeneralSumOfSquares.name,
+                GeneralSumOfEvenPowers.name]
             permute_signs_check = [
-                'homogeneous_ternary_quadratic',
-                'homogeneous_ternary_quadratic_normal',
-                'binary_quadratic']
+                HomogeneousTernaryQuadratic.name,
+                HomogeneousTernaryQuadraticNormal.name,
+                BinaryQuadratic.name]
             if t in permute_signs_for:
                 do_permute_signs_var = True
             elif t in permute_signs_check:
@@ -413,17 +614,17 @@ def diophantine(eq, param=symbols("t", integer=True), syms=None,
         solution = diop_solve(base, param)
 
         if eq_type in [
-                "linear",
-                "homogeneous_ternary_quadratic",
-                "homogeneous_ternary_quadratic_normal",
-                "general_pythagorean"]:
+                Linear.name,
+                HomogeneousTernaryQuadratic.name,
+                HomogeneousTernaryQuadraticNormal.name,
+                GeneralPythagorean.name]:
             sols.add(merge_solution(var, var_t, solution))
 
         elif eq_type in [
-                "binary_quadratic",
-                "general_sum_of_squares",
-                "general_sum_of_even_powers",
-                "univariate"]:
+                BinaryQuadratic.name,
+                GeneralSumOfSquares.name,
+                GeneralSumOfEvenPowers.name,
+                Univariate.name]:
             for sol in solution:
                 sols.add(merge_solution(var, var_t, sol))
 
@@ -536,28 +737,28 @@ def diop_solve(eq, param=symbols("t", integer=True)):
     """
     var, coeff, eq_type = classify_diop(eq, _dict=False)
 
-    if eq_type == "linear":
+    if eq_type == Linear.name:
         return diop_linear(eq, param)
 
-    elif eq_type == "binary_quadratic":
+    elif eq_type == BinaryQuadratic.name:
         return diop_quadratic(eq, param)
 
-    elif eq_type == "homogeneous_ternary_quadratic":
+    elif eq_type == HomogeneousTernaryQuadratic.name:
         return diop_ternary_quadratic(eq, parameterize=True)
 
-    elif eq_type == "homogeneous_ternary_quadratic_normal":
+    elif eq_type == HomogeneousTernaryQuadraticNormal.name:
         return diop_ternary_quadratic_normal(eq, parameterize=True)
 
-    elif eq_type == "general_pythagorean":
+    elif eq_type == GeneralPythagorean.name:
         return diop_general_pythagorean(eq, param)
 
-    elif eq_type == "univariate":
+    elif eq_type == Univariate.name:
         return diop_univariate(eq)
 
-    elif eq_type == "general_sum_of_squares":
+    elif eq_type == GeneralSumOfSquares.name:
         return diop_general_sum_of_squares(eq, limit=S.Infinity)
 
-    elif eq_type == "general_sum_of_even_powers":
+    elif eq_type == GeneralSumOfEvenPowers.name:
         return diop_general_sum_of_even_powers(eq, limit=S.Infinity)
 
     if eq_type is not None and eq_type not in diop_known:
@@ -574,69 +775,17 @@ def diop_solve(eq, param=symbols("t", integer=True)):
 
 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()
-    if not all(_is_int(c) for c in coeff.values()):
-        raise TypeError("Coefficients should be Integers")
 
+    matched = False
     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 total_degree == 2 and len(var) == 2:
-        diop_type = "binary_quadratic"
-
-    elif total_degree == 2 and len(var) == 3 and homogeneous:
-        if set(coeff) & set(var):
-            diop_type = "inhomogeneous_ternary_quadratic"
-        else:
-            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"
+    for diop_class in all_diop_classes:
+        diop_type = diop_class(eq)
+        if diop_type.matches():
+            matched = True
+            break
 
-    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"
-
-    elif total_degree == 3 and len(var) == 2:
-        diop_type = "cubic_thue"
-
-    elif (total_degree > 3 and total_degree % 2 == 0 and
-            all(k.is_Pow and k.exp == total_degree 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'
-
-    if diop_type is not None:
-        return var, dict(coeff) if _dict else coeff, diop_type
+    if matched:
+        return diop_type.free_symbols, dict(diop_type.coeff) if _dict else diop_type.coeff, diop_type.name
 
     # new diop type instructions
     # --------------------------
@@ -733,7 +882,7 @@ def diop_linear(eq, param=symbols("t", integer=True)):
     """
     var, coeff, diop_type = classify_diop(eq, _dict=False)
 
-    if diop_type == "linear":
+    if diop_type == Linear.name:
         result = _diop_linear(var, coeff, param)
 
         if param is None:
@@ -1011,7 +1160,7 @@ def diop_univariate(eq):
     """
     var, coeff, diop_type = classify_diop(eq, _dict=False)
 
-    if diop_type == "univariate":
+    if diop_type == Univariate.name:
         return set([(int(i),) for i in solveset_real(
             eq, var[0]).intersect(S.Integers)])
 
@@ -1068,7 +1217,7 @@ def diop_quadratic(eq, param=symbols("t", integer=True)):
     """
     var, coeff, diop_type = classify_diop(eq, _dict=False)
 
-    if diop_type == "binary_quadratic":
+    if diop_type == BinaryQuadratic.name:
         return set(_diop_quadratic(var, coeff, param))
 
 
@@ -1928,7 +2077,7 @@ def transformation_to_DN(eq):
     """
 
     var, coeff, diop_type = classify_diop(eq, _dict=False)
-    if diop_type == "binary_quadratic":
+    if diop_type == BinaryQuadratic.name:
         return _transformation_to_DN(var, coeff)
 
 
@@ -2023,7 +2172,7 @@ def find_DN(eq):
            http://www.jpr2718.org/ax2p.pdf
     """
     var, coeff, diop_type = classify_diop(eq, _dict=False)
-    if diop_type == "binary_quadratic":
+    if diop_type == BinaryQuadratic.name:
         return _find_DN(var, coeff)
 
 
@@ -2468,7 +2617,7 @@ def diop_ternary_quadratic_normal(eq, parameterize=False):
     (4, 9, 1)
     """
     var, coeff, diop_type = classify_diop(eq, _dict=False)
-    if diop_type == "homogeneous_ternary_quadratic_normal":
+    if diop_type == HomogeneousTernaryQuadraticNormal.name:
         sol = _diop_ternary_quadratic_normal(var, coeff)
         if len(sol) > 0:
             x_0, y_0, z_0 = list(sol)[0]
@@ -2916,7 +3065,7 @@ def diop_general_pythagorean(eq, param=symbols("m", integer=True)):
     """
     var, coeff, diop_type  = classify_diop(eq, _dict=False)
 
-    if diop_type == "general_pythagorean":
+    if diop_type == GeneralPythagorean.name:
         return list(_diop_general_pythagorean(var, coeff, param))[0]
 
 
@@ -2991,7 +3140,7 @@ def diop_general_sum_of_squares(eq, limit=1):
     """
     var, coeff, diop_type = classify_diop(eq, _dict=False)
 
-    if diop_type == "general_sum_of_squares":
+    if diop_type == GeneralSumOfSquares.name:
         return set(_diop_general_sum_of_squares(var, -int(coeff[1]), limit))
 
 
@@ -3050,7 +3199,7 @@ def diop_general_sum_of_even_powers(eq, limit=1):
     """
     var, coeff, diop_type = classify_diop(eq, _dict=False)
 
-    if diop_type == "general_sum_of_even_powers":
+    if diop_type == GeneralSumOfEvenPowers.name:
         for k in coeff.keys():
             if k.is_Pow and coeff[k]:
                 p = k.exp
diff --git a/sympy/solvers/diophantine/tests/test_diophantine.py b/sympy/solvers/diophantine/tests/test_diophantine.py
index d4c31b8e5b1b..8e5e8cc8a717 100644
--- a/sympy/solvers/diophantine/tests/test_diophantine.py
+++ b/sympy/solvers/diophantine/tests/test_diophantine.py
@@ -70,6 +70,8 @@ def test_classify_diop():
         [x, y], {x*y**2: 1, 1: 1}, 'cubic_thue')
     assert classify_diop(x**4 + y**4 + z**4 - (1 + 16 + 81)) == (
         [x, y, z], {1: -98, x**4: 1, z**4: 1, y**4: 1}, 'general_sum_of_even_powers')
+    assert classify_diop(x**2 + y**2 + z**2) == (
+        [x, y, z], {x**2: 1, y**2: 1, z**2: 1}, 'homogeneous_ternary_quadratic_normal')
 
 
 def test_linear():


Code review comments:
diff_hunk: @@ -152,6 +136,223 @@ def __call__(self, *args):
         return self.subs(zip(self.parameters, args))
 
 
+class DiophantineEquationType:
+    """
+    Internal representation of a particular diophantine equation type.
+
+    Parameters
+    ==========
+
+    equation
+        The diophantine equation that is being solved.
+    free_symbols: list (optional)
+        The symbols being solved for.
+
+    Attributes
+    ==========
+
+    total_degree
+        The maximum of the degrees of all terms in the equation
+    homogeneous
+        Does the equation contain a term of degree 0
+    homogeneous_order
+        Does the equation contain any coefficient that is in the symbols being solved for
+    dimension
+        The number of symbols being solved for
+    """
+    name = None
+
+    def __init__(self, equation, free_symbols=None):
+        try:
+            if free_symbols is None:
+                free_symbols = list(equation.free_symbols)
+            assert free_symbols
+        except (AttributeError, AssertionError):
+            raise ValueError('equation should have 1 or more free symbols')
path: sympy/solvers/diophantine/diophantine.py
text: @user1:
I realise that you're just moving the code around here but I really don't like catching either `AttributeError` or `AssertionError` as these can hide bugs. Is it possible to remove them? Why is `AttributeError` being caught anyway?

@author:
It's raised with the test `raises(ValueError, lambda: classify_diop(1))`. Probably it would be better to just sympify the argument.