{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "8d70673e",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "b40aa144",
   "metadata": {},
   "outputs": [],
   "source": [
    "n , k = symbols('n, k')\n",
    "c6 = Function('c6')\n",
    "c5 = Function('c5')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "c6ca9089",
   "metadata": {},
   "outputs": [],
   "source": [
    "f0 = lambda n, k: 1\n",
    "f1 = lambda n, k: k"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "46443345",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle k^{2} - n$"
      ],
      "text/plain": [
       "k**2 - n"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f2 = lambda n, k: k * f1(n, k) - 1 * n * f0(n, k)\n",
    "expand(f2(n, k))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "fcd1fc66",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle k^{3} - 3 k n$"
      ],
      "text/plain": [
       "k**3 - 3*k*n"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f3 = lambda n, k:  k * f2(n, k) - 2 * n * f1(n, k)\n",
    "expand(f3(n, k))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "7eac2eae",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle k^{4} - 6 k^{2} n + 3 n^{2} + 2 n$"
      ],
      "text/plain": [
       "k**4 - 6*k**2*n + 3*n**2 + 2*n"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f4 = lambda n, k:  k * f3(n, k) - 3 * n * f2(n, k) + 2 * n\n",
    "expand(f4(n, k))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "28d52775",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle k^{4} - 6 k^{2} n + 3 n^{2} + 2 n$"
      ],
      "text/plain": [
       "k**4 - 6*k**2*n + 3*n**2 + 2*n"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expand((f4(n+1, k+1) + f4(n+1, k-1))/2 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "f29cda28",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle k^{5} - 10 k^{3} n + 15 k n^{2} + 10 k n$"
      ],
      "text/plain": [
       "k**5 - 10*k**3*n + 15*k*n**2 + 10*k*n"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f5 = lambda n, k:  k * f4(n, k) - 4 * n * f3(n, k) + 8 * n * k\n",
    "expand(f5(n, k))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "b9d28d1f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle k^{5} - 10 k^{3} n + 15 k n^{2} + 10 k n$"
      ],
      "text/plain": [
       "k**5 - 10*k**3*n + 15*k*n**2 + 10*k*n"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expand((f5(n+1, k+1) + f5(n+1, k-1))/2 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "cd683968",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle k^{6} - 15 k^{4} n + 45 k^{2} n^{2} + 30 k^{2} n - 15 n^{3} - 30 n^{2} - 16 n$"
      ],
      "text/plain": [
       "k**6 - 15*k**4*n + 45*k**2*n**2 + 30*k**2*n - 15*n**3 - 30*n**2 - 16*n"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f6 = lambda n, k:  k * f5(n, k) - 5 * n * f4(n, k) + 20 * n * k ** 2 - 16 * n - 20 * n ** 2\n",
    "expand(f6(n, k))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "54d584a4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 0$"
      ],
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expand((f6(n+1, k+1) + f6(n+1, k-1))/2 ) - expand(f6(n, k))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "642a6707",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 4 n \\left(5 k^{2} - 5 n - 4\\right)$"
      ],
      "text/plain": [
       "4*n*(5*k**2 - 5*n - 4)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "simplify(20 * n * k ** 2 - 36 * n - 20 * n ** 2 + 20 * n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "33bf0b2b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 4 n \\left(5 k^{2} - 5 n - 4\\right)$"
      ],
      "text/plain": [
       "4*n*(5*k**2 - 5*n - 4)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "2 * n * f0(n,k)\n",
    "8 * n * f1(n, k)\n",
    "simplify(20 * n *(f2(n, k)) - 16 * n * f0(n,k))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "953338e8",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "a652f30e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle k^{6} - 15 k^{4} n + 45 k^{2} n^{2} + 30 k^{2} n - 15 n^{3} - 30 n^{2} - 16 n$"
      ],
      "text/plain": [
       "k**6 - 15*k**4*n + 45*k**2*n**2 + 30*k**2*n - 15*n**3 - 30*n**2 - 16*n"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f7 = lambda n, k:  k * f6(n, k) - 6 * n * f5(n, k) + 40 * n * f3(n, k) - 96 * n * f1(n, k)\n",
    "expand(f6(n, k))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "adacc2eb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 0$"
      ],
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expand((f7(n+1, k+1) + f7(n+1, k-1))/2 ) - expand(f7(n, k))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "53da835d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 8 k n \\left(5 k^{2} - 15 n - 12\\right)$"
      ],
      "text/plain": [
       "8*k*n*(5*k**2 - 15*n - 12)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "simplify(40 * n * f3(n, k) - 96 * n * f1(n, k))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "c8e80d85",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 0$"
      ],
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "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)\n",
    "expand((f8(n+1, k+1) + f8(n+1, k-1))/2 ) - expand(f8(n, k))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "225c4f0e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 0$"
      ],
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "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)\n",
    "expand((f9(n+1, k+1) + f9(n+1, k-1))/2 ) - expand(f9(n, k))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "2aeec120",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 0$"
      ],
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "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)\n",
    "expand((f10(n+1, k+1) + f10(n+1, k-1))/2 ) - expand(f10(n, k))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "33f4b228",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 0$"
      ],
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "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)\n",
    "expand((f11(n+1, k+1) + f11(n+1, k-1))/2 ) - expand(f11(n, k))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "0f627091",
   "metadata": {},
   "outputs": [],
   "source": [
    "import functools\n",
    "\n",
    "\n",
    "def f(m):\n",
    "    if m < 0:\n",
    "        return lambda n, k: 0\n",
    "    if m == 0:\n",
    "        return lambda n, k: 1\n",
    "    elif m == 1:\n",
    "        return lambda n, k: k\n",
    "    else:\n",
    "        \n",
    "        def alternating_tangent(n):\n",
    "            # https://mathworld.wolfram.com/TangentNumber.html \n",
    "            return 2 ** (2*n) * (2**(2*n)-1) * Abs(bernoulli(2*n)) / (2*n)\n",
    "        \n",
    "        def R(n, k):\n",
    "            running_sum = 0\n",
    "            for j in range(1, round(m/2) + 1):\n",
    "                i = m - 2*j\n",
    "                running_sum += (-1)**(j) * alternating_tangent(j) * binomial(m-1, (2 * j - 1)) * f(i)(n, k)\n",
    "            return running_sum\n",
    "        \n",
    "        def fm(n, k):\n",
    "            return k * f(m-1)(n,k) + n * R(n, k)\n",
    "        \n",
    "        \n",
    "        return fm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "id": "1f48d7c9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "k^{7} - 21 k^{5} n + 105 k^{3} n^{2} + 70 k^{3} n - 105 k n^{3} - 210 k n^{2} - 112 k n\n"
     ]
    }
   ],
   "source": [
    "print(latex(expand(f(8)(n,k))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "2c21c9ca",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle -9600$"
      ],
      "text/plain": [
       "-9600"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f(6)(6,6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 167,
   "id": "fb2d7ade",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle k^{4} - 6 k^{2} n - 6 k^{2} + 3 n^{2} + 8 n + 5$"
      ],
      "text/plain": [
       "k**4 - 6*k**2*n - 6*k**2 + 3*n**2 + 8*n + 5"
      ]
     },
     "execution_count": 167,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expand(f(4)(n+1,k))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 168,
   "id": "e99d75c6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle k^{4} - 6 k^{2} n - 6 k^{2} + 3 n^{2} + 8 n + 5$"
      ],
      "text/plain": [
       "k**4 - 6*k**2*n - 6*k**2 + 3*n**2 + 8*n + 5"
      ]
     },
     "execution_count": 168,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expand(f(4)(n,k) - 6 * f(2)(n, k) + 5*f(0)(n, k))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 193,
   "id": "1087f30d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle k^{5} - 10 k^{3} n - 10 k^{3} + 15 k n^{2} + 40 k n + 25 k$"
      ],
      "text/plain": [
       "k**5 - 10*k**3*n - 10*k**3 + 15*k*n**2 + 40*k*n + 25*k"
      ]
     },
     "execution_count": 193,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expand(f(5)(n+1,k))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 194,
   "id": "12557e5a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle k^{5} - 10 k^{3} n - 10 k^{3} + 15 k n^{2} + 40 k n + 25 k$"
      ],
      "text/plain": [
       "k**5 - 10*k**3*n - 10*k**3 + 15*k*n**2 + 40*k*n + 25*k"
      ]
     },
     "execution_count": 194,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expand(f(5)(n,k) -10 * f(3)(n,k) + 25 * f(1)(n,k))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 195,
   "id": "74ad8942",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle k^{6} - 15 k^{4} n - 15 k^{4} + 45 k^{2} n^{2} + 120 k^{2} n + 75 k^{2} - 15 n^{3} - 75 n^{2} - 121 n - 61$"
      ],
      "text/plain": [
       "k**6 - 15*k**4*n - 15*k**4 + 45*k**2*n**2 + 120*k**2*n + 75*k**2 - 15*n**3 - 75*n**2 - 121*n - 61"
      ]
     },
     "execution_count": 195,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expand(f(6)(n+1,k))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 206,
   "id": "11cc0ff2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle k^{6} - 15 k^{4} n - 15 k^{4} + 45 k^{2} n^{2} + 120 k^{2} n + 75 k^{2} - 15 n^{3} - 75 n^{2} - 121 n - 61$"
      ],
      "text/plain": [
       "k**6 - 15*k**4*n - 15*k**4 + 45*k**2*n**2 + 120*k**2*n + 75*k**2 - 15*n**3 - 75*n**2 - 121*n - 61"
      ]
     },
     "execution_count": 206,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expand(f(6)(n,k) - 15* f(4)(n,k) + 75 * f(2)(n, k) - 61 * f(0)(n,k))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 208,
   "id": "d1725992",
   "metadata": {},
   "outputs": [],
   "source": [
    "# the above coefficients look like http://www.luschny.de/math/seq/SwissKnifePolynomials.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 220,
   "id": "a7d5b904",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - 8 k^{3} + 24 k n - 8 k$"
      ],
      "text/plain": [
       "-8*k**3 + 24*k*n - 8*k"
      ]
     },
     "execution_count": 220,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expand(f(4)(n,k-1) - f(4)(n,k+1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 218,
   "id": "5d7ac620",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle k^{4} - 4 k^{3} - 6 k^{2} n + 6 k^{2} + 12 k n - 4 k + 3 n^{2} - 4 n + 1$"
      ],
      "text/plain": [
       "k**4 - 4*k**3 - 6*k**2*n + 6*k**2 + 12*k*n - 4*k + 3*n**2 - 4*n + 1"
      ]
     },
     "execution_count": 218,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "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))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e11bf6b",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.0rc2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}