Added 'limits of the real roots'

Author: cfgnunes

README.md

@@ -57,6 +57,7 @@ python3 main.py
 
 - Briot-Ruffini method
 - Newton's Divided-Difference method
+- Limits of the real roots
 
 ### Numerical differentiation
 

main.py

@@ -252,6 +252,20 @@ def run_example_neville():
     print(f"\tq =\n{q}")
 
 
+@print_docstring
+def run_example_root_limits():
+    """Run an example 'Polynomials: Root limits'."""
+    c = np.array([1, 2, -13, -14, 24])
+
+    print("Input:")
+    print(f"\tc = {c}")
+
+    limits = polynomials.root_limits(c)
+
+    print("Output:")
+    print(f"\tlimits = {limits}")
+
+
 @print_docstring
 def run_example_briot_ruffini():
     """Run an example 'Polynomials: Briot-Ruffini'."""
@@ -656,6 +670,7 @@ def main():
     run_example_newton()
     run_example_lagrange()
     run_example_neville()
+    run_example_root_limits()
     run_example_briot_ruffini()
     run_example_newton_divided_difference()
     run_example_derivative_backward_difference()

polynomials.py

@@ -1,5 +1,7 @@
 """Methods for polynomials."""
 
+import math
+
 import numpy as np
 
 
@@ -66,3 +68,76 @@ def newton_divided_difference(x, y):
     print("")
 
     return [f]
+
+
+def root_limits(c):
+    """Find the limits of the real roots of a polynomial equation.
+
+    Using Lagrange's Theorem, whose proof is given by Demidovich and Maron.
+
+    Args:
+        c: polynomial coefficients.
+
+    Returns:
+        lim: lower and upper limits of positive and negative roots.
+    """
+    lim = np.zeros(4)
+    n = len(c) - 1
+    c = np.concatenate(([0], c))
+    c = np.concatenate((c, [0]))
+
+    if c[1] == 0:
+        raise ValueError("The coefficient c[1] is null.")
+
+    t = n + 1
+    c[t + 1] = 0
+
+    # If c[t+1] is null, then the polynomial is deflated.
+    while True:
+        if c[t] != 0:
+            break
+        t -= 1
+
+    # Compute the four limits of real roots.
+    for i in range(0, 4):
+        if i in (1, 3):
+            # Inversion of the order of the coefficients.
+            for j in range(1, t // 2 + 1):
+                c[j], c[t - j + 1] = c[t - j + 1], c[j]
+        else:
+            if i == 2:
+                # Reinversion of the order and exchange
+                # of signs of the coefficients.
+                for j in range(1, t // 2 + 1):
+                    c[j], c[t - j + 1] = c[t - j + 1], c[j]
+                for j in range(t - 1, 0, -2):
+                    c[j] = -c[j]
+
+        # If c[1] is negative, then all coefficients are swapped.
+        if c[1] < 0:
+            for j in range(1, t + 1):
+                c[j] = -c[j]
+
+        # Calculation of 'k', the largest index of the negative coefficients.
+        k = 2
+        while True:
+            if c[k] < 0 or k > t:
+                break
+            k += 1
+
+        # Calculation of 'b', the largest negative coefficient in modulus.
+        if k <= t:
+            b = 0
+            for j in range(2, t + 1):
+                if c[j] < 0 and math.fabs(c[j]) > b:
+                    b = math.fabs(c[j])
+
+            # Limit of positive roots of 'P(x) = 0' and auxiliary equations.
+            lim[i] = 1 + (b / c[1]) ** (1 / (k - 1))
+        else:
+            lim[i] = 10 ** 100
+
+    # Limit of positive and negative roots of 'P(x) = 0'.
+    lim[0], lim[1], lim[2], lim[3] = 1 / lim[1], lim[0], -lim[2], -1 / lim[3]
+
+    return lim