# Just for fun ;)
# Usage - if you have recurrent equation u(n) + a_k-1 u(n-1) + ... + a_0 u(n-k) = 0
# and condition that u(0) = b_0, u(1) = b_1, ... , u(k-1) = b_k-1
# then you just put the coefficients 1, a_k-1, a_k-2, ... , a_0 to initial_recurrent_coeffs list
# and put b_0, b_1, ... , b_k-1 to first_values list

# For me, I have an equation 1*u(n) + 2*u(n-1) + 1*u(n-2) = 0, and u(0) = -1, u(1) = -1,
# so I put 1, 2, 1 into initial_recurrent_coeffs and -1, -1 into first_values

initial_recurrent_coeffs = [1, 2, 1]
# Add first k-1 values of recurrent sequence
first_values = [-1, -1]

import sympy
from sympy import Symbol
from sympy import symbols
from sympy import poly
from sympy import Idx, Indexed, Matrix
from sympy.solvers import solve_linear_system

def solve_recurrent(coeffs):
    from sympy.abc import x
    ex = 0
    for i in xrange(len(coeffs)):
        ex += coeffs[i]*x**(len(coeffs) - 1 - i)
    return ex

def roots_counts(polynom):
    roots = polynom.all_roots()
    counts = {}
    for i in xrange(len(roots)):
        if counts.has_key(roots[i]):
            counts[roots[i]] += 1
        else:
            counts[roots[i]] = 1
    return counts

def solve_linear(counts, first_values):
    pow_values = counts.values()

    # Warning! I didn't check that str(pow_values[i]) is right choice for cc making 
    cc = [list(symbols('C'+str(i)+':'+str(i+1)+':'+str(pow_values[i]))) for i in xrange(len(pow_values))]
    
    counts_list = list(counts)
    print 'counts list: ', counts_list
    
    under_matrix = [[] for i in xrange(sum(counts.values()))]
    
    for power_of_roots, ex in enumerate(under_matrix):
        for i_idx, root_value in enumerate(counts):
            for j_idx in xrange(counts[root_value]):
                under_matrix[power_of_roots].append(power_of_roots**j_idx * root_value**power_of_roots)

    for i, item in enumerate(under_matrix):
        item.append(first_values[i])
    linear_system = Matrix(under_matrix)

    all_variables = [elem for item in cc for elem in item]
    
    solution = solve_linear_system(linear_system, *all_variables)
    
    # Build recurrent function
    result_expression = 0
    
    x = Symbol('x')
    
    print cc
    
    root_values = list(counts.keys())
    for root_idx, root_coeffs in enumerate(cc):
        subresult_expression = 0
        for coeff_idx, coeff_symbol in enumerate(root_coeffs):
            subresult_expression += solution[cc[root_idx][coeff_idx]] * x**coeff_idx
        subresult_expression *= root_values[root_idx]**x
        result_expression += subresult_expression

    print result_expression
    print sympy.simplify(result_expression)

    return sympy.lambdify(x, result_expression)


def naive_recursive_equation(n):
    if n < 0: return 0
    if n == 0: return -2
    if n == 1: return 1
    return 3*naive_recursive_equation(n-1) - 2*naive_recursive_equation(n-2)

def naive_recursive_equation_b(n):
    if n < 0: return 0
    if n == 0: return -1
    if n == 1: return -1
    return -2*naive_recursive_equation_b(n-1) - naive_recursive_equation_b(n-2)

# Make polynom from recurrent coefficients
characteristic_polynom = sympy.poly_from_expr(solve_recurrent(initial_recurrent_coeffs))[0]

# find all roots of polynom with powers of each root
counts = roots_counts(characteristic_polynom)

# solve recurrent equations and give all Cij, that exactly describes solution
recurrent_generator = solve_linear(counts, first_values)

for i in xrange(20):
    print recurrent_generator(i), naive_recursive_equation_b(i)


# Nice job, my friend 8-)