real_root_isolation.py

import typing
import sys
import sympy as sp
import numpy as np
from sympy.parsing.sympy_parser import parse_expr

x, y = sp.symbols('x y')


result = []


def real_root_isolation(p, c, d):

    n = p.degree()
    pI = (y + 1)**n * p.subs({x: (c * y + d) / (y + 1)})
    pI_coeffs = sp.Poly.from_expr(pI).simplify().coeffs()
    mu = count_sign_changes(pI_coeffs)

    print("intermediate result:", result, "\n")
    print("c:", c, "d:", d)
    print("coefficients:", pI_coeffs)

    if mu == 0:
        result.append((c, d, "-> keine Nullstelle"))
        return
    if mu == 1:
        result.append((c, d))
        return

    m = (c + d) / 2
    p_at_m = p.evalf(subs={x: m})
    if p_at_m == 0:
        result.append(m)

    real_root_isolation(p, c, m)
    real_root_isolation(p, m, d)


def count_sign_changes(a):
    sign_changes = 0
    last_sign = np.sign(a[0])

    for i in range(1, len(a)):
        this_sign = np.sign(a[i])
        if this_sign != last_sign:
            sign_changes += 1
            last_sign = this_sign
    return sign_changes


def calculate_b(p):
    coeffs = p.coeffs()[::-1]

    last = coeffs[len(coeffs) - 1]

    l = []
    for i in range(0, len(coeffs)):
        new_element = coeffs[i] / last
        new_element = abs(new_element) # entspricht nicht der Beschreibung auf dem Übungsblatt, aber funktioniert so.
        l.append(new_element)

    print(l)
    return 1 + max(l)


print("Enter Polynomial:")
in_string = sys.stdin.readline()
plnm = sp.Poly.from_expr(parse_expr(in_string))

b = calculate_b(plnm)
print("b", b)

real_root_isolation(plnm, -b, b)
print(result)