gamma_pole_series.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Symbolically computing the series of $\\Gamma(-n+\\varepsilon)$ with https://functions.wolfram.com/GammaBetaErf/Gamma/06/01/05/01/0006/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import *\n",
    "from sympy.abc import n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def c(i):\n",
    "    if i % 2 == 0:\n",
    "        return 0\n",
    "    k = i // 2\n",
    "    r = Number(-1)**k * 2 *(Number(2)**(2*k + 1) - 1) * bernoulli(2*k+2) * pi**(2*k + 1) / factorial(2*k + 2)\n",
    "    return r\n",
    "\n",
    "def b(k):\n",
    "    return Derivative(gamma(n+1), n, k) / (factorial(n) * factorial(k))\n",
    "\n",
    "def p(j, k):\n",
    "    if k == 0:\n",
    "        return 1\n",
    "    su = 0\n",
    "    for m in range(1, k+1):\n",
    "        su += (j*m + m - k) * b(m) * p(j, k-m)\n",
    "    return su/k"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def term(k):\n",
    "    ss = 0\n",
    "    for j in range(0, k+2):\n",
    "        su = 0\n",
    "        for r in range(0, j+1):\n",
    "            su += Number(-1)**(j+r) * binomial(j, r) * p(r, j) / (r + 1)\n",
    "        ss += (j + 1) * su * c(k-j)\n",
    "    return pi * ss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{\\operatorname{polygamma}^{5}{\\left(0,n + 1 \\right)}}{120} - \\frac{\\operatorname{polygamma}^{3}{\\left(0,n + 1 \\right)} \\operatorname{polygamma}{\\left(1,n + 1 \\right)}}{12} + \\frac{\\pi^{2} \\operatorname{polygamma}^{3}{\\left(0,n + 1 \\right)}}{36} + \\frac{\\operatorname{polygamma}^{2}{\\left(0,n + 1 \\right)} \\operatorname{polygamma}{\\left(2,n + 1 \\right)}}{12} + \\frac{\\operatorname{polygamma}{\\left(0,n + 1 \\right)} \\operatorname{polygamma}^{2}{\\left(1,n + 1 \\right)}}{8} - \\frac{\\pi^{2} \\operatorname{polygamma}{\\left(0,n + 1 \\right)} \\operatorname{polygamma}{\\left(1,n + 1 \\right)}}{12} - \\frac{\\operatorname{polygamma}{\\left(0,n + 1 \\right)} \\operatorname{polygamma}{\\left(3,n + 1 \\right)}}{24} + \\frac{7 \\pi^{4} \\operatorname{polygamma}{\\left(0,n + 1 \\right)}}{360} - \\frac{\\operatorname{polygamma}{\\left(1,n + 1 \\right)} \\operatorname{polygamma}{\\left(2,n + 1 \\right)}}{12} + \\frac{\\pi^{2} \\operatorname{polygamma}{\\left(2,n + 1 \\right)}}{36} + \\frac{\\operatorname{polygamma}{\\left(4,n + 1 \\right)}}{120}$"
      ],
      "text/plain": [
       "polygamma(0, n + 1)**5/120 - polygamma(0, n + 1)**3*polygamma(1, n + 1)/12 + pi**2*polygamma(0, n + 1)**3/36 + polygamma(0, n + 1)**2*polygamma(2, n + 1)/12 + polygamma(0, n + 1)*polygamma(1, n + 1)**2/8 - pi**2*polygamma(0, n + 1)*polygamma(1, n + 1)/12 - polygamma(0, n + 1)*polygamma(3, n + 1)/24 + 7*pi**4*polygamma(0, n + 1)/360 - polygamma(1, n + 1)*polygamma(2, n + 1)/12 + pi**2*polygamma(2, n + 1)/36 + polygamma(4, n + 1)/120"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# one needs to simplify twice, naturally\n",
    "expr = term(4).simplify().simplify()\n",
    "expr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{3 \\psi_{0}^{5} + \\psi_{0}^{3} \\left(- 30 \\psi_{1} + 10 \\pi^{2}\\right) + 30 \\psi_{0}^{2} \\psi_{2} + \\psi_{0} \\left(45 \\psi_{1}^{2} - 30 \\pi^{2} \\psi_{1} - 15 \\psi_{3} + 7 \\pi^{4}\\right) - 30 \\psi_{1} \\psi_{2} + 10 \\pi^{2} \\psi_{2} + 3 \\psi_{4}}{360}$"
      ],
      "text/plain": [
       "(3*\\psi_0**5 + \\psi_0**3*(-30*\\psi_1 + 10*pi**2) + 30*\\psi_0**2*\\psi_2 + \\psi_0*(45*\\psi_1**2 - 30*pi**2*\\psi_1 - 15*\\psi_3 + 7*pi**4) - 30*\\psi_1*\\psi_2 + 10*pi**2*\\psi_2 + 3*\\psi_4)/360"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "subsdict = {polygamma(i, n+1): Symbol(\"\\\\psi_%d\"%i) for i in range(0, 10)}\n",
    "expr.subs(subsdict).factor(Symbol(\"\\\\psi_0\"))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.9.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}