# Functions for Linear Problems

# ----------------------------------------------------------------------------------------------------------------------------- #

# for a parameter list 'P', and the minimal polynomial of a 'd'-th matrix 'Rd', and a 'd'-the vector 'Ud'
# find 'Vd' or 'KERd' such that 'Rd * Ud = Vd + Kd' where 'Kd' is an element of the kernerl of 'Rd'
# if 'vector_opt == 0' then compute only 'Kd' 

def linear_problem_solver(P, Rd, Ud, Cd, vector_opt=0, lin_log_opt=0):
    # a preparation
    if lin_log_opt != 0:
       TIME0 = cputime(subprocesses=True)
    C = Cd.coefficients(sparse=False)
    x = parent(Cd).gen()
    d = len(C)
    j = floor(sqrt(float(d)))
    # the vactor 'Ud' of 'L[2](PHI[d-1]) + ... + L[d-1](PHI[2])'
    if C[0] != 0:
       c0 = C[0]
       C.pop(0)
       Vd = expanded_horner_method(P, V, -add([(C[i]/c0)*(x**(i)) for i in range(d-1)]), Rd, Ud, j)
       if lin_log_opt != 0:
          TIME1 = cputime(subprocesses=True)
          print(lin_log_opt+1, '-th Linear Algebra Problem (Vector):', TIME1-TIME0, 'seconds')
          print(lin_log_opt+1, '-th Linear Algebra Problem (Vector):', Vd)
       Kd = 0
       if lin_log_opt != 0:
          TIME2 = cputime(subprocesses=True)
          print(lin_log_opt+1, '-th Linear Algebra Problem (Kernel):', TIME2-TIME1, 'seconds')
          print(lin_log_opt+1, '-th Linear Algebra Problem (Kernel):', Kd)
    else:
       c1 = C[1]
       C.pop(0)
       if vector_opt == 0:
          Vd = expanded_horner_method(P, V, -add([(C[i+1]/c1)*(x**i) for i in range(d-2)]), Rd, Ud, j)
          if lin_log_opt != 0:
             TIME1 = cputime(subprocesses=True)
             print(lin_log_opt+1, '-th Linear Algebra Problem (Vector):', TIME1-TIME0, 'seconds')
             print(lin_log_opt+1, '-th Linear Algebra Problem (Vector):', Vd)
          Kd = expanded_horner_method(P, V, -add([(C[i]/c1)*(x**i) for i in range(d-1)]), Rd, Ud, j)
          if lin_log_opt != 0:
             TIME2 = cputime(subprocesses=True)
             print(lin_log_opt+1, '-th Linear Algebra Problem (Kernel):', TIME2-TIME1, 'seconds')
             print(lin_log_opt+1, '-th Linear Algebra Problem (Kernel):', Kd)
       else:
          Vd = 0
          if lin_log_opt != 0:
             TIME1 = cputime(subprocesses=True)
             print(lin_log_opt+1, '-th Linear Algebra Problem (Vector):', TIME1-TIME0, 'seconds')
             print(lin_log_opt+1, '-th Linear Algebra Problem (Vector):', Vd)
          Kd = expanded_horner_method(P, V, -add([(C[i]/c1)*(x**(i)) for i in range(d-1)]), Rd, Ud, j)
          if lin_log_opt != 0:
             TIME2 = cputime(subprocesses=True)
             print(lin_log_opt+1, '-th Linear Algebra Problem (Kernel):', TIME2-TIME1, 'seconds')
             print(lin_log_opt+1, '-th Linear Algebra Problem (Kernel):', Kd)
    return [Vd, Kd]

# ----------------------------------------------------------------------------------------------------------------------------- #

# for a parameter list 'P', and a variable list 'V',
#     and a polynomial 'G', and a matrix 'Rd', and a vector 'Ud', and a natural number 'j'
# generate the vector 'Fd(Rd * Ud)' by using the expaned horner's method
# where 'Fd(x) = an * (x**n) + ... + a0 = bk(x) * (x**(d*k)) + ... + b0(x)'

def expanded_horner_method(P, V, Fd, Rd, Ud, j):
    # 'k+1' lists of the '0,...,j'-th, 'j+1,...,2j'-th, ..., 'k*j+1,...,(k+1)*j'-th, '(k+1)*j+1, ..., n'-th coefficients of 'Fd'
    n = Fd.degree()
    k = floor(n/j)
    Cj = Fd.coefficients(sparse=False)
    Cj = [Cj[i*j:i*j+j] for i in range(k+1)]
    PF = parent(Cj[0][0])
    # the list of '(Rd**(i+1))*Ud' for i in range(j), the 'j'-th power 'Rd**j' of 'Rd'
    Rd_p = [Ud] + [(Rd**(i+1))*Ud for i in range(j-1)]
    Rd_j = Rd**(j)
    # the polynomial preresentaion of 'b0(Rd*Ud), ..., Bk(Rd*Ud)'
    Pd = [sum([cj[i]*Rd_p[i] for i in range(len(cj))]) for cj in Cj]
    # the sum of '(Rd**(d*k))*(bk(Rd*Ud)) + ... + b0(Rd*Ud)'
    Hj = Pd[-1]
    Pd = Pd[:-1]
    Hj = recursive_expanded_horner_method(Hj, Ud, Pd, Rd_j)
    Hj = [PF((SR(h[0]).numerator()/(SR(h[0]).denominator()))) for h in Hj]
    Hj = [h.numerator()*((h.denominator()).denominator())/(h.denominator()).numerator() for h in Hj]
    return Hj

def recursive_expanded_horner_method(Hj, Ud, Pd, Rd_p):
    if len(Pd) != 0:
       Hj = Rd_p*Hj + Pd[-1]
       Pd = Pd[:-1]
       return(recursive_expanded_horner_method(Hj, Ud, Pd, Rd_p))
    else:
       return(Hj)

# ----------------------------------------------------------------------------------------------------------------------------- #

def routh_hurwitz_determinant(f):
    d = f.degree()
    l = list(reversed(f.list()))
    H = matrix([[l[2*i-j+1] if 2*i-j+1 >= 0 and 2*i-j+1 <= d else 0 for j in range(d)] for i in range(d)])
    return [l[-1], [H[0:k,0:k].det() for k in range(1,d)]]

# ----------------------------------------------------------------------------------------------------------------------------- #