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

# The investigation below suggests that the found polynomial families are appell sequences <=> are have translation properties <=> have interesting derivative properties 
# 
# 1. https://en.wikipedia.org/wiki/Bernoulli_polynomials#Differences_and_derivatives
# 2. https://en.wikipedia.org/wiki/Appell_sequence
# 
# Q: are $p_m(n, k)$ linearly independent? Sol: show $p_m(0, k) = k^m$. Then this uses the fact the $x^m$ are linearly independent polynomials. Simple using the recurrence definition. 
# 
# Q: given a polynomial $p$ that is BH, is it a linear comb. of $\{p_m\}$? Given any polynomial in (n,k), is it a linear combination of $\{p_m\}$? Scratch: the bernoulli polynomials have an inversion formula, take a look at that. 
# 
# 

# In[1]:


from sympy import *


# In[66]:


x = symbols('x')
k = symbols('k')
a = symbols('a')
n = symbols('n')


# In[161]:


import functools


def bh(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)) * bh(i)(n, k)
            return running_sum
        
        def fm(n, k):
            return k * bh(m-1)(n,k) + n * R(n, k)
        
        return fm
    
bh_ = lambda m: (lambda x: bh(m)(n, x)) # fixing n


# In[167]:


def f(m):
    def fm(n, k):
        return (k + n) ** (m)
    return fm

f_ = lambda m: (lambda x: f(m)(n, x)) # fixing n


# In[169]:


for s in [f_, bh_]:
    assert expand(s(2)(x-2)) == expand(s(2)(x) - 2 * s(1)(x) * 2 + s(0)(x) * 4)
    assert expand(s(3)(x+1)) == expand(s(3)(x) + 3 * s(2)(x) + 3 * s(1)(x) + s(0)(x))


# In[181]:


# satifies appell seq. 
assert expand(diff(bh_(3)(x), x)) == expand(3 * bh_(2)(x))
assert expand(diff(bh_(4)(x), x)) == expand(4 * bh_(3)(x))
assert expand(diff(bh_(5)(x), x)) == expand(5 * bh_(4)(x))


# In[228]:


collect( (expand(bh(4)(1, k+1) + bh(4)(1, k-1)))/2, n)


# In[230]:


expand(f(3)(1, k + 1) + f(3)(1, k - 1)) 


# In[203]:


expand(2 *binomial(3,0)*bh(3)(n, k) + 2 *binomial(3,1)*bh(1)(n, k))


# In[231]:


expand(bh(4)(0, k))


# In[234]:


expand(bh(4)(n+1, k+1) + bh(3)(n+1, k+1) +  bh(4)(n+1, k-1) + bh(3)(n+1, k-1) )


# In[237]:


expand(2 * bh(4)(n, k) + 2*bh(3)(n, k))


# In[238]:


expand(bh(4)(n+1, k+1) + bh(3)(n+1, k+1))


# In[240]:


expand(bh(4)(n+1, k-1) + bh(3)(n+1, k-1))


# In[242]:


expand(bh(4)(n, k) )


# In[243]:


expand(bh(6)(n, k))


# In[246]:


for m in range(7):
    print(m, expand(bh(m)(n, k)))


# In[252]:


expand(bh(6)(n, k))


# In[ ]: