Last active
October 28, 2023 15:16
-
-
Save Seniatical/ebb3d897c2db18f7339cfd80e574c1de to your computer and use it in GitHub Desktop.
Factorising univariate polynomials
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
""" | |
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()) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Refactored some variable names and took changed snippets from @cherriae 's code