Factorising univariate polynomials

UniPolynomialFactorising.py

"""
Author: Isa (Seniatical)
License: GNU General Public License (GNU GPL)

An algorithmic approach to factorising univariate polynomials,
by combining methods such as the factor theorm to solve factors in an efficient manor.

Finding factors:
~~~~~~~~~~~~~~~~
This algorithm finds factors of all univariate polynomials which are factorisable, i.e. have roots.
The method this algorithm uses is better suited for cubics and other expressions which continue,
for quadratics and linear expressions, it is more simpler and faster to solve directly.

It is far simpler and faster to solve for roots by finding multiples from the first or last term,
as these are the result of multiplication, not addition.

This algorithm uses the last term, which is a real number with no variable.
By listing all of its factors we can repeatdly test whether any factors can divise the polynomial with no remainder,
if any do we thus have a factor of the polynomial.
We then repeat this step until we reach a linear or quadratic expression which returns our final factors.

If the last term is negative, we use the absolute value to find its factors and adapt the theorm later in the algorithm.

Example Steps:
1> IN: x^3 - 17x^2 + 82x - 120
2> LAST_TERM = 120
3> FACTORS = {1, 2, 3, 4, 5, 6, 40, 8, 10, 12, 15, 20, 120, 24, 60, 30}
4> IN.DIVIDE(..., N) -> IN.DIVIDE(-3)
5> OUT : x^2 -14x + 40
6> FACTOR OUT [QUADRATIC]
-4, -10
7> FACTORS = [-3, -4, -10]

For the full pseudocode example of the algorithm see below!
Full example written in python below.


Testing if polynomial factors are completely negative:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
If a polynomials factors are completely negative then we can test by reading every next coefficient in the polynomial,
if they are all negative then we can say all factors are negative.

Quadratic proof
(x - 3)(x - 4) -> x^2 - 4x - 3x + 12
               -> x^2 - 7x + 12
                      ^^^^ 

Cubic Proof
(x - 3)(x - 4)(x - 10) -> x^3 - 17x^2 + 82x - 120
                              ^^^^^^^       ^^^^^

Quartic proof
(x - 2)(x - 3)(x - 4)(x - 10) -> x^4 - 19x^3 + 116x^2 - 284x + 240
                                     ^^^^^^^          ^^^^^^


Example Test:
```
function factors_negative(poly[N...])
    is_odd = false

    for coefficient in poly
        if is_odd == false
            is_odd = true
            continue
        endif

        if coefficient >= 0
            return false
        endif

        is_odd = false
    endfor

    return true

endfunction
```

this test will indicicate whether we need to negate the multiple of P[N - 1] or not when testing the factor theorm

---------------------------------------------------------------------------------------------------------------

Example algorithm structure:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~
```
function factors(poly[N...])
    if poly.degree == 1
        return factor_linear(poly)
    else if poly.degree == 2
        return factor_quadratic(poly)
    endif

    right_most = poly[-1]
    is_negative = right_most < 0
    if is_negative:
        right_most = right_most.absolute()
    right_multiples = factorize(right_most)
    factors = []

    factors_neg = factors_negative(poly)

    n_poly = poly
    for multiple in right_multiples
        if poly.sub(-multiple) == 0 and not factors_neg
            factors.push(multiple)
            n_poly = n_poly.divide(multiple)
            endloop
        else if poly.sub(multiple) == 0 and factors_neg
            factors.push(-multiple)
            n_poly = n_poly.divide(factors_neg)
            endloop
        endif
    endfor
    
    if factors.length == poly.degree or factors.length == 0
        return factors
    endif
    
    factors.push_many(n_poly.factors())
    return factors

endfunction
```

References:
Polynomial Division: https://www.rosettacode.org/wiki/Polynomial_long_division
Factoring real numbers: https://stackoverflow.com/a/6800214

****
NOTE
****
This approach is subject to tweaks and changes to improve its accuracy and efficiency.
"""
from __future__ import annotations
from typing import *

from functools import reduce


def factors(n):
    ''' See REFERENCES:Factorising real numbers, for a better explanation '''
    return set(reduce(list.__add__, 
                ([i, n//i] for i in range(1, int(n**0.5) + 1) if not n % i)))


def factor_linear(*args) -> float:
    ''' Factors a linear equation, takes 2 arguments **m** and **c** '''
    m, c = args
    return (-c) / m


def factor_quadratic(*args) -> float:
    ''' Solves a quadratic equation, takes 2 arguments **a**, **b** and **c** '''
    a, b, c = args
    s1 = ((b * -1) + (((b ** 2) - (4 * a * c)) ** 0.5)) / (2 * a)
    s2 = ((b * -1) - (((b ** 2) - (4 * a * c)) ** 0.5)) / (2 * a)

    ## return whole factor
    if a > 1:
        if s1 - int(s1) != 0:
            s1 *= a
        elif s2 - int(s2) != 0:
            s2 *= a

    return s1, s2


def degree(poly):
    ''' Returns the degree of a polynomial '''
    while poly and poly[-1] == 0:
        poly.pop()
    return len(poly)-1


class Polynomial(object):
    '''
    An object that represents a polynomial.
    
    ..code-block:: python
    
        >>> Polynomial(3, 4)
        Polynomial[3, 4] # y = 3x + 4
        >>> Polynomial(1, 6, 9) # y = x^2 + 6x + 9
        ...
    
    Parameters
    ----------
    *coefficients: *float
        Coefficients for each term in the polynomials,
        degree of polynomial is automatically calculated after.
    '''
    coefficients: List[float]
    '''coefficients of polynomial'''
    degree: int
    '''degree of polynomial'''

    def __init__(self, *coefficients) -> None:
        self.coefficients = list(coefficients)
    
    @property
    def degree(self) -> int:
        return degree(self.coefficients)

    def copy(self) -> Polynomial:
        '''Returns a shallow copy of the polynomial'''
        return Polynomial(*self.coefficients.copy())

    def sub(self, X: float) -> float:
        '''
        SEE REFERENCES:Polynomial Division for more info on how this functions works.
        
        Divides a polynomial by another polynomial,
        for example (x^2 + 6x + 9) / (x + 3) -> Polynomial.divide(Polynomial(1, 6, 9), Polynomial(1, 3)).
        
        Returns `Q`, `R` or the quotient and remainder.
        -> Returns Q/R as polynomial if possible
        -> Returns Q/R as a float/None
        . Varies on results.
        '''
        return sum(
            (X**power) * coefficient
            for power, coefficient in enumerate(self.coefficients[::-1])
        )

    @classmethod
    def divide(cls, left: Polynomial, right: Polynomial) -> List[List[float], List[float]]:
        '''
        Divides a polynomial by another polynomial,
        for example (x^2 + 6x + 9) / (x + 3) -> Polynomial.divide(Polynomial(1, 6, 9), Polynomial(1, 3))
        '''
        left_d = left.coefficients.copy()
        right_d = right.coefficients.copy()

        if right.degree < 0:
            raise ZeroDivisionError()
        dL = left.degree
        dR = right.degree
        
        if dL >= dR:
            q = [0] * dL
            while dL >= dR:
                d = ([0] * (dL - dR)) + right_d
                multiple = q[dL - dR] = left_d[-1] / float(right_d[-1])
                d = [coefficient * multiple for coefficient in d]
                left_d = [cL - cd for cL, cd in zip(left_d, d)]
                dL = degree(left_d)
        else:
            q = [0]

        r = left_d
        q = Polynomial(*q) if degree(q) > 0 else (q or [None])[0]
        r = Polynomial(*r) if degree(r) > 0 else (r or [None])[0]

        return q, r

    

    def factors_negative(self) -> bool:
        ''' 
        Tests whether current polynomials factors are all negatives,
        by skipping a term per iteration and checking if the coefficients are all negative.
        If not then one or more factors may be positive
        '''
        is_odd = False
        
        for coefficient in self.coefficients:
            if not is_odd:
                is_odd = True
                continue
            
            if coefficient >= 0:
                return False
            is_odd = False
        
        return True


    def factors(self) -> List[float]:
        '''
        Calculates the factors of a standard 1D univariate polynomial
        '''
        match (self.degree):
            case (1):
                return [-factor_linear(*self.coefficients)]
            case (2):
                # Calculates root, we need factor
                r1, r2 = factor_quadratic(*self.coefficients)
                return [(r1 * -1), (r2 * -1)]

        right_most = self.coefficients[-1]
        if right_most < 0:
            right_most = abs(right_most)
        right_multiples = factors(right_most)
        factors_negative = self.factors_negative()

        F = []
        nPoly = self
        ## Test factors using factor theorm
        for multiple in right_multiples:
            if self.sub(-multiple) == 0 and not factors_negative:
                F.append(multiple)
                nPoly, R = Polynomial.divide(nPoly, Polynomial(1, multiple))
                break
            elif factors_negative and self.sub(multiple) == 0:
                F.append(-multiple)
                nPoly, R = Polynomial.divide(nPoly, Polynomial(1, -multiple))
                break
        
        ## Breaking point, e.g. cannot reduce no more
        if nPoly == self or not nPoly:
            return F
        F.extend(nPoly.factors())
        return F
            

    def __repr__(self) -> str:
        return f'Polynomial{self.coefficients}'
    
    def __eq__(self, other: Polynomial) -> bool:
        return set(self.coefficients) == set(other.coefficients)


if __name__ == '__main__':
    print('Testing Division')
    print(Polynomial.divide(Polynomial(1, 6, 9), Polynomial(1, 3)))
    print(Polynomial.divide(Polynomial(1, 6, 9), Polynomial(1, 3)) == (Polynomial(1, 3), None))
    print()

    print('Calculating factors')
    # (x + 3)(x + 4)(x + 1)
    P = Polynomial(1, 8, 19, 12)
    print(P, '=', P.factors())

    # (x + 3)(x - 4)(x + 3)
    P = Polynomial(1, 2, -15, -36)
    print(P, '=', P.factors())

    # (x + 3)(3x - 4)(x + 3)
    P = Polynomial(3, 14, 3, -36)
    print(P, '=', P.factors())

    # (x + 2)(x + 3)(x + 4)(x + 10)
    P = Polynomial(1, 15, 48, -44, -240)
    print(P, '=', P.factors())

    # (x - 2)(x - 3)(x - 4)(x - 10)
    P = Polynomial(1, -19, 116, -284, 240)
    print(P, '=', P.factors())

    # (x - 4)(x - 3)(x - 2.5)
    P = Polynomial(1, -9.5, 29.5, -30)
    print(P, '=', P.factors())

Did some refactoring for you:

from __future__ import annotations
from typing import *

from functools import reduce


def factors(n):
    """ See REFERENCES:Factorising real numbers, for a better explanation """
    return set(reduce(list.__add__,
                      ([i, n // i] for i in range(1, int(n ** 0.5) + 1) if not n % i)))


def factor_linear(*args) -> float:
    """ Factors a linear equation, takes 2 arguments **m** and **c** """
    m, c = args
    return (-c) / m


def factor_quadratic(*args) -> Tuple[float | Any, float | Any]:
    """ Solves a quadratic equation, takes 2 arguments **a**, **b** and **c** """
    a, b, c = args
    s1 = ((b * -1) + (((b ** 2) - (4 * a * c)) ** 0.5)) / (2 * a)
    s2 = ((b * -1) - (((b ** 2) - (4 * a * c)) ** 0.5)) / (2 * a)

    ## return whole factor
    if a > 1:
        if s1 - int(s1) != 0:
            s1 *= a
        elif s2 - int(s2) != 0:
            s2 *= a

    return s1, s2


def degree(poly):
    """ Returns the degree of a polynomial """
    while poly and poly[-1] == 0:
        poly.pop()
    return len(poly) - 1


class Polynomial(object):
    """
    An object that represents a polynomial.

    ..code-block:: python

        >>> Polynomial(3, 4)
        Polynomial[3, 4] # y = 3x + 4
        >>> Polynomial(1, 6, 9) # y = x^2 + 6x + 9
        ...

    Parameters
    ----------
    *coefficients: *float
        Coefficients for each term in the polynomials,
        degree of polynomial is automatically calculated after.
    """
    coefficients: List[float]
    '''coefficients of polynomial'''
    degree: int
    '''degree of polynomial'''

    def __init__(self, *coefficients) -> None:
        self.coefficients = list(coefficients)

    @property
    def degree(self) -> int:
        return degree(self.coefficients)

    def copy(self) -> Polynomial:
        """Returns a shallow copy of the polynomial"""
        return Polynomial(*self.coefficients.copy())

    def sub(self, X: float) -> float:
        """
        SEE REFERENCES:Polynomial Division for more info on how this functions works.

        Divides a polynomial by another polynomial,
        for example (x^2 + 6x + 9) / (x + 3) -> Polynomial.divide(Polynomial(1, 6, 9), Polynomial(1, 3)).

        Returns `Q`, `R` or the quotient and remainder.
        -> Returns Q/R as polynomial if possible
        -> Returns Q/R as a float/None
        . Varies on results.
        """
        return sum(
            (X**power) * coefficient
            for power, coefficient in enumerate(self.coefficients[::-1])
        )

    @classmethod
    def divide(cls, left: Polynomial, right: Polynomial) -> Tuple[Polynomial | int, Polynomial | float]:
        """
        Divides a polynomial by another polynomial,
        for example (x^2 + 6x + 9) / (x + 3) -> Polynomial.divide(Polynomial(1, 6, 9), Polynomial(1, 3))
        """
        left_d = left.coefficients.copy()
        right_d = right.coefficients.copy()

        if right.degree < 0:
            raise ZeroDivisionError()
        dL = left.degree
        dR = right.degree

        if dL >= dR:
            q = [0] * dL
            while dL >= dR:
                d = ([0] * (dL - dR)) + right_d
                multiple = q[dL - dR] = left_d[-1] / float(right_d[-1])
                d = [coefficient * multiple for coefficient in d]
                left_d = [cL - cd for cL, cd in zip(left_d, d)]
                dL = degree(left_d)
        else:
            q = [0]
        r = left_d
        q = Polynomial(*q) if degree(q) > 0 else (q or [None])[0]
        r = Polynomial(*r) if degree(r) > 0 else (r or [None])[0]
        return q, r

    def factors_negative(self) -> bool:
        """
        Tests whether current polynomials factors are all negatives,
        by skipping a term per iteration and checking if the coefficients are all negative.
        If not then one or more factors may be positive
        """
        is_odd = False

        for coefficient in self.coefficients:
            if not is_odd:
                is_odd = True
                continue

            if coefficient >= 0:
                return False
            is_odd = False

        return True

    def factors(self) -> List[float]:
        """
        Calculates the factors of a standard 1D univariate polynomial
        """
        match self.degree:
            case (1):
                return [-factor_linear(*self.coefficients)]
            case (2):
                # Calculates root, we need factor
                r1, r2 = factor_quadratic(*self.coefficients)
                return [(r1 * -1), (r2 * -1)]

        RIGHT_MOST = self.coefficients[-1]
        if RIGHT_MOST < 0:
            RIGHT_MOST = abs(RIGHT_MOST)
        RIGHT_MULTIPLES = factors(RIGHT_MOST)
        FACTORS_NEGATIVE = self.factors_negative()

        F = []
        nPoly = self
        ## Test factors using factor theorm
        for multiple in RIGHT_MULTIPLES:
            if self.sub(-multiple) == 0 and not FACTORS_NEGATIVE:
                F.append(multiple)
                nPoly, R = Polynomial.divide(nPoly, Polynomial(1, multiple))
                break
            elif FACTORS_NEGATIVE and self.sub(multiple) == 0:
                F.append(-multiple)
                nPoly, R = Polynomial.divide(nPoly, Polynomial(1, -multiple))
                break

        ## Breaking point, e.g. cannot reduce no more
        if nPoly == self or not nPoly:
            return F
        F.extend(nPoly.factors())
        return F

    def __repr__(self) -> str:
        return f'Polynomial{self.coefficients}'

    def __eq__(self, other: Polynomial) -> bool:
        return set(self.coefficients) == set(other.coefficients)

Refactored some variable names and took changed snippets from @cherriae 's code