#!/usr/bin/env python
# coding: utf-8

# In[1]:


from sympy import *


# In[2]:


n , k = symbols('n, k')
c6 = Function('c6')
c5 = Function('c5')


# In[3]:


f0 = lambda n, k: 1
f1 = lambda n, k: k


# In[4]:


f2 = lambda n, k: k * f1(n, k) - 1 * n * f0(n, k)
expand(f2(n, k))


# In[5]:


f3 = lambda n, k:  k * f2(n, k) - 2 * n * f1(n, k)
expand(f3(n, k))


# In[6]:


f4 = lambda n, k:  k * f3(n, k) - 3 * n * f2(n, k) + 2 * n
expand(f4(n, k))


# In[7]:


expand((f4(n+1, k+1) + f4(n+1, k-1))/2 )


# In[8]:


f5 = lambda n, k:  k * f4(n, k) - 4 * n * f3(n, k) + 8 * n * k
expand(f5(n, k))


# In[9]:


expand((f5(n+1, k+1) + f5(n+1, k-1))/2 )


# In[10]:


f6 = lambda n, k:  k * f5(n, k) - 5 * n * f4(n, k) + 20 * n * k ** 2 - 16 * n - 20 * n ** 2
expand(f6(n, k))


# In[11]:


expand((f6(n+1, k+1) + f6(n+1, k-1))/2 ) - expand(f6(n, k))


# In[12]:


simplify(20 * n * k ** 2 - 36 * n - 20 * n ** 2 + 20 * n)


# In[13]:


2 * n * f0(n,k)
8 * n * f1(n, k)
simplify(20 * n *(f2(n, k)) - 16 * n * f0(n,k))


# In[ ]:





# In[14]:


f7 = lambda n, k:  k * f6(n, k) - 6 * n * f5(n, k) + 40 * n * f3(n, k) - 96 * n * f1(n, k)
expand(f6(n, k))


# In[15]:


expand((f7(n+1, k+1) + f7(n+1, k-1))/2 ) - expand(f7(n, k))


# In[16]:


simplify(40 * n * f3(n, k) - 96 * n * f1(n, k))


# In[17]:


f8 = lambda n, k:  k * f7(n, k) - 7 * n * f6(n, k) + 70 * n * f4(n, k) - 336 * n * f2(n, k) + 272 * n * f0(n, k)
expand((f8(n+1, k+1) + f8(n+1, k-1))/2 ) - expand(f8(n, k))


# In[18]:


f9 = lambda n, k:  k * f8(n, k) - 8 * n * f7(n, k) + 112 * n * f5(n, k) - 896 * n * f3(n, k) + 2176 * n * f1(n, k)
expand((f9(n+1, k+1) + f9(n+1, k-1))/2 ) - expand(f9(n, k))


# In[19]:


f10 = lambda n, k:  k * f9(n, k) - 9 * n * f8(n, k) + 168 * n * f6(n, k) - 2016 * n * f4(n, k) + 9792 * n * f2(n, k) - 7936 * n * f0(n, k)
expand((f10(n+1, k+1) + f10(n+1, k-1))/2 ) - expand(f10(n, k))


# In[20]:


f11 = lambda n, k:  k * f10(n, k) - 10 * n * f9(n, k) + 240 * n * f7(n, k) - 4032 * n * f5(n, k) + 32640 * n * f3(n, k) - 79360 * n * f1(n, k)
expand((f11(n+1, k+1) + f11(n+1, k-1))/2 ) - expand(f11(n, k))


# In[21]:


import functools


def f(m):
    if m < 0:
        return lambda n, k: 0
    if m == 0:
        return lambda n, k: 1
    elif m == 1:
        return lambda n, k: k
    else:
        
        def alternating_tangent(n):
            # https://mathworld.wolfram.com/TangentNumber.html 
            return 2 ** (2*n) * (2**(2*n)-1) * Abs(bernoulli(2*n)) / (2*n)
        
        def R(n, k):
            running_sum = 0
            for j in range(1, round(m/2) + 1):
                i = m - 2*j
                running_sum += (-1)**(j) * alternating_tangent(j) * binomial(m-1, (2 * j - 1)) * f(i)(n, k)
            return running_sum
        
        def fm(n, k):
            return k * f(m-1)(n,k) + n * R(n, k)
        
        
        return fm


# In[126]:


print(latex(expand(f(8)(n,k))))


# In[23]:


f(6)(6,6)


# In[167]:


expand(f(4)(n+1,k))


# In[168]:


expand(f(4)(n,k) - 6 * f(2)(n, k) + 5*f(0)(n, k))


# In[193]:


expand(f(5)(n+1,k))


# In[194]:


expand(f(5)(n,k) -10 * f(3)(n,k) + 25 * f(1)(n,k))


# In[195]:


expand(f(6)(n+1,k))


# In[206]:


expand(f(6)(n,k) - 15* f(4)(n,k) + 75 * f(2)(n, k) - 61 * f(0)(n,k))


# In[208]:


# the above coefficients look like http://www.luschny.de/math/seq/SwissKnifePolynomials.html


# In[220]:


expand(f(4)(n,k-1) - f(4)(n,k+1))


# In[218]:


expand(f(4)(n,k) - 4 * f(3)(n, k) + 6*f(2)(n, k) - 4*f(1)(n, k) + 1 * f(0)(n, k))


# In[ ]: