polynomials

# Modesti Dennis VR429663
# Polinomi.py

from typing import List
real = float
polynomial = List[real]


def degree(pol: polynomial) -> int:
    if pol == []:
        return -1
    last_zero = 0
    for (n, i) in enumerate(pol):
        if i != 0:
            last_zero = n
    return last_zero


def extend(pol: polynomial, deg: int) -> polynomial:
    return pol + [0 for _ in range(deg - (len(pol) - 1))]


def reduction(pol: polynomial) -> polynomial:
    return pol[:degree(pol) + 1]


def sum(a: polynomial, b: polynomial) -> polynomial:
    deg_max = max(degree(a), degree(b))
    pol_a = extend(a, deg_max)
    pol_b = extend(b, deg_max)
    pol_sum = []
    for (i, j) in zip(pol_a, pol_b):
        pol_sum.append(i + j)
    return reduction(pol_sum)


def negative(pol: polynomial) -> polynomial:
    return reduction([-i for i in pol])


def subtraction(a: polynomial, b: polynomial) -> polynomial:
    return sum(a, negative(b))


def scalar_product(pol: polynomial, mul: real) -> polynomial:
    return reduction([i * mul for i in pol])


# Returns factorial(n) / factorial(d), n > d
def multi_factorial(n: int, d: int = 1) -> int:
    val = 1
    while n > d:
        val *= n
        n -= 1
    return val


def derivative(pol: polynomial, amount: int = 1) -> polynomial:
    return reduction(
        [multi_factorial(n + amount, n) * i for (n, i) in enumerate(pol[amount:])]
        )


def apply(pol: polynomial, value: real) -> real:
    val = 0
    for (n, i) in enumerate(pol):
        val += value**n * i
    return val


def is_equal(a: polynomial, b: polynomial) -> bool:
    return reduction(a) == reduction(b)


def increase_deg(pol: polynomial, deg: int = 1) -> polynomial:
    return [0 for _ in range(deg)] + pol


def product(a: polynomial, b: polynomial) -> polynomial:
    pol = []
    for (n, i) in enumerate(a):
        pol = sum(pol, scalar_product(increase_deg(b, n), i))
    return pol


# Returns the sign of a*b
def multi_sign(a: real, b: real) -> int:
    return [1, -1][(a < 0) ^ (b < 0)]


def find_zero(
    pol: polynomial,
    low_extreme: real = -1e4,
    high_extreme: real = 1e4,
    epsilon: real = 0.001
) -> real:
    if pol == [0]:
        return 0
    if pol == []:
        return None
    if low_extreme == high_extreme:
        return None

    f_le = apply(pol, low_extreme)
    f_he = apply(pol, high_extreme)
    m = (low_extreme + high_extreme) / 2

    if abs(f_le) < epsilon:
        return low_extreme

    if abs(f_he) < epsilon:
        return high_extreme

    if multi_sign(f_le, f_he) == 1:
        m = find_zero(derivative(pol), low_extreme, high_extreme, epsilon)
        if m == None:
            return None
        if multi_sign(f_le, apply(pol, m)) == 1:
            return None

    while abs(low_extreme - m) > epsilon:
        f_m = apply(pol, m)
        if multi_sign(f_le, f_m) == -1:
            high_extreme = m
        else:
            low_extreme = m
            f_le = f_m
        m = (low_extreme + high_extreme) / 2
    return m


def create_monomial(val: real, deg: int) -> polynomial:
    return [0 for _ in range(deg)] + [val]


def divide(pol: polynomial, d: polynomial) -> (polynomial, polynomial):
    dividend = reduction(pol)
    divisor = reduction(d)
    if divisor == [0] or divisor == []:
        return None
    quotient = []
    while degree(dividend) >= degree(divisor):
        divide_coeff = dividend[-1] / divisor[-1]
        bottom_step = create_monomial(
            divide_coeff, degree(dividend) - degree(divisor))
        dividend = subtraction(dividend, product(bottom_step, divisor))
        quotient = sum(quotient, bottom_step)
    return (quotient, dividend)