Polonomial division

algorithms/maths/polynomial.py

@@ -4,7 +4,7 @@
 from typing import Dict, Union, Set, Iterable
 from numbers import Rational
 from functools import reduce
-
+from math import inf
 
 class Monomial:
     """
@@ -283,6 +283,56 @@ def __str__(self) -> str:
             result += temp
         return result + ')'
 
+    def degree(self):
+        """
+        Get the degree of the monomial (sum of the exponents)
+        """
+        return sum(self.variables.values())
+
+    def degree_with_respect_to(self, other):
+        """
+        Same as degree(), except if this monomial shares no variables with other, then return -inf.
+        """
+        for other_variable, other_degree in other.variables.items():
+            if other_variable not in self.variables:
+                return -inf
+            self_degree = self.variables[other_variable]
+            if self_degree < other_degree:
+                return -inf
+
+        return self.degree()
+
+    def graded_reverse_lexicographic_order(self):
+        """
+        Sorts monomials first by their degree, and then by the reverse order of their variables.
+
+        For example, monomials would be sorted as follows:
+
+            x^2 > xy > y^2 > xz > yz > z^2
+
+        Reference: https://en.wikipedia.org/wiki/Monomial_order#Graded_reverse_lexicographic_order
+        """
+        return (self.degree(), self.variables_flattened_reverse())
+
+    def variables_flattened_reverse(self):
+        """
+        Given a monomial, returns the variables names as a list (keeping the
+        exponent as the number of repetitions).
+
+        For example, the monomial:
+            x^3 y^2 z^4
+        would turn into:
+            xxxyyzzzz
+        and then reversed:
+            zzzzyyxxx
+        """
+        names = []
+        for variable, degree in self.variables.items():
+            names += [variable] * degree
+        names.sort()
+        names.reverse()
+        return names
+
 
 class Polynomial:
     """
@@ -429,40 +479,90 @@ def __mul__(self, other: Union[int, float, Fraction, Monomial]):
             raise ValueError('Can only multiple int, float, Fraction, Monomials, or Polynomials with Polynomials.')
 
     # def __floordiv__(self, other: Union[int, float, Fraction, Monomial, Polynomial]) -> Polynomial:
-    def __floordiv__(self, other: Union[int, float, Fraction, Monomial]):
+    def __floordiv__(self, other: Union[int, float, Fraction, Monomial, 'Polynomial']):
         """
         For Polynomials, floordiv is the same
         as truediv.
         """
+        if isinstance(other, Polynomial):
+            if len(other.all_monomials()) == 1:
+                monomial = other.all_monomials().pop()
+                return self.__truediv__(monomial)
+            quotient, remainder = self.polynomial_division(other)
+            return quotient
         return self.__truediv__(other)
 
     # def __truediv__(self, other: Union[int, float, Fraction, Monomial, Polynomial]) -> Polynomial:
-    def __truediv__(self, other: Union[int, float, Fraction, Monomial]):
+    def __truediv__(self, other: Union[int, float, Fraction, Monomial, 'Polynomial']):
         """
         For Polynomials, only division by a monomial
         is defined.
-
-        TODO: Implement polynomial / polynomial.
         """
         if isinstance(other, int) or isinstance(other, float) or isinstance(other, Fraction):
             return self.__truediv__( Monomial({}, other) )
         elif isinstance(other, Monomial):
             poly_temp = reduce(lambda acc, val: acc + val, map(lambda x: x / other, [z for z in self.all_monomials()]), Polynomial([Monomial({}, 0)]))
             return poly_temp
         elif isinstance(other, Polynomial):
-            if Monomial({}, 0) in other.all_monomials():
-                if len(other.all_monomials()) == 2:
-                    temp_set = {x for x in other.all_monomials() if x != Monomial({}, 0)}
-                    only = temp_set.pop()
-                    return self.__truediv__(only)
-            elif len(other.all_monomials()) == 1:
-                temp_set = {x for x in other.all_monomials()}
-                only = temp_set.pop()
-                return self.__truediv__(only)
+            if len(other.all_monomials()) == 1:
+                monomial = other.all_monomials().pop()
+                return self.__truediv__(monomial)
+            quotient, remainder = self.polynomial_division(other)
+            assert len(remainder.all_monomials()) == 0, "polynomial division yielded non-zero remainder"
+            return quotient
 
         raise ValueError('Can only divide a polynomial by an int, float, Fraction, or a Monomial.')
+
+    def polynomial_division(self, other):
+        """
+        Perform polynomial division.
+        Order of division is performed following this rating:
+        1. Dividend term shares variable with divisor term of higest degree.
+        2. Degree
+        ie: if divisor is xy^2z^3 then a term of the dividend is divisible iff it is a produt of
+        xy^2z^2 and it is of the same or a (total) higher degree.
+        """
+
+        # Check that all monomial degrees are non-negative
+        monomials = self.all_monomials() | other.all_monomials()
+        if any((degree <= 0 for monomial in monomials for degree in monomial.variables.values())):
+            raise ValueError("cannot divide polynomials with negative or zero exponent")
+
+        # Implementation using standard long division using a Gröbner basis
+        # Reference: https://en.wikipedia.org/wiki/Polynomial_long_division
+        # Reference: https://en.wikipedia.org/wiki/Gr%C3%B6bner_basis#Reduction
+        # Reference: https://en.wikipedia.org/wiki/Monomial_order#Graded_reverse_lexicographic_order
+        remainder = self
+        quotient = Polynomial([])
+        while len(remainder.all_monomials()) > 0:
+            # Find the monomial with highest degree in the divisor
+            other_max = max(other.all_monomials(), key=lambda m: m.graded_reverse_lexicographic_order())
+
+            # Grab the term/monomial (with the highest degree) from the dividend that shares a
+            # variable with the highest degree term/monomial from the divisor.
+            dividend_max = max(remainder.all_monomials(), key = lambda term: term.degree_with_respect_to(other_max))
+
+            # Continue until the remainder's degree cannot be reduced further, or no common variable between the dividend and divisor remains
+            if dividend_max.degree_with_respect_to(other_max) == -inf:
+                break
+
+            # Find the factor required to reduce the remainders degree by one (or more)
+            factor = dividend_max / other_max
+            remainder -= other * factor
+            quotient += factor
+
+        return quotient, remainder
 
-        return
+    @staticmethod
+    def check_only_one_variable(monomial: Monomial) -> bool:
+        count = 0
+        for variable, degree in monomial.variables.items():
+            if degree == 0:
+                continue
+            count += 1
+
+        # Can either contain zero or one variables
+        return count <= 1
 
     # def clone(self) -> Polynomial:
     def clone(self):
@@ -531,4 +631,3 @@ def __str__(self) -> str:
         the polynomial.
         """
         return ' + '.join(str(m) for m in self.all_monomials() if m.coeff != Fraction(0, 1))
-