{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gauss-Laguerre Quadratur für die Gamma-Funktion\n",
    "\n",
    "$$\n",
    "    \\Gamma(z)\n",
    "    = \n",
    "    \\int_0^\\infty t^{z-1}e^{-t}dt\n",
    "$$\n",
    "\n",
    "$$\n",
    "    \\int_0^\\infty f(x) e^{-x} dx \n",
    "    \\approx \n",
    "    \\sum_{i=1}^{N} f(x_i) w_i\n",
    "    \\qquad\\text{ wobei }\n",
    "    w_i = \\frac{x_i}{(n+1)^2 [L_{n+1}(x_i)]^2}\n",
    "$$\n",
    "und $x_i$ sind Nullstellen des Laguerre Polynoms $L_n(x)$\n",
    "\n",
    "Der Fehler ist gegeben als\n",
    "\n",
    "$$\n",
    "    E \n",
    "    =\n",
    "    \\frac{(n!)^2}{(2n)!} f^{(2n)}(\\xi) \n",
    "    = \n",
    "    \\frac{(-2n + z)_{2n}}{(z-m)_m} \\frac{(n!)^2}{(2n)!} \\xi^{z + m - 2n - 1}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from cmath import exp, pi, sin, sqrt\n",
    "import scipy.special\n",
    "\n",
    "EPSILON = 1e-07\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lanczos_p = [\n",
    "    676.5203681218851,\n",
    "    -1259.1392167224028,\n",
    "    771.32342877765313,\n",
    "    -176.61502916214059,\n",
    "    12.507343278686905,\n",
    "    -0.13857109526572012,\n",
    "    9.9843695780195716e-6,\n",
    "    1.5056327351493116e-7,\n",
    "]\n",
    "\n",
    "\n",
    "def drop_imag(z):\n",
    "    if abs(z.imag) <= EPSILON:\n",
    "        z = z.real\n",
    "    return z\n",
    "\n",
    "\n",
    "def lanczos_gamma(z):\n",
    "    z = complex(z)\n",
    "    if z.real < 0.5:\n",
    "        y = pi / (sin(pi * z) * lanczos_gamma(1 - z))  # Reflection formula\n",
    "    else:\n",
    "        z -= 1\n",
    "        x = 0.99999999999980993\n",
    "        for (i, pval) in enumerate(lanczos_p):\n",
    "            x += pval / (z + i + 1)\n",
    "        t = z + len(lanczos_p) - 0.5\n",
    "        y = sqrt(2 * pi) * t ** (z + 0.5) * exp(-t) * x\n",
    "    return drop_imag(y)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zeros, weights = np.polynomial.laguerre.laggauss(8)\n",
    "# zeros = np.array(\n",
    "#     [\n",
    "#         1.70279632305101000e-1,\n",
    "#         9.03701776799379912e-1,\n",
    "#         2.25108662986613069e0,\n",
    "#         4.26670017028765879e0,\n",
    "#         7.04590540239346570e0,\n",
    "#         1.07585160101809952e1,\n",
    "#         1.57406786412780046e1,\n",
    "#         2.28631317368892641e1,\n",
    "#     ]\n",
    "# )\n",
    "\n",
    "# weights = np.array(\n",
    "#     [\n",
    "#         3.69188589341637530e-1,\n",
    "#         4.18786780814342956e-1,\n",
    "#         1.75794986637171806e-1,\n",
    "#         3.33434922612156515e-2,\n",
    "#         2.79453623522567252e-3,\n",
    "#         9.07650877335821310e-5,\n",
    "#         8.48574671627253154e-7,\n",
    "#         1.04800117487151038e-9,\n",
    "#     ]\n",
    "# )\n",
    "\n",
    "\n",
    "def pochhammer(z, n):\n",
    "    return np.prod(z + np.arange(n))\n",
    "\n",
    "\n",
    "def find_shift(z, target):\n",
    "    factor = 1.0\n",
    "    steps = int(np.floor(target - np.real(z)))\n",
    "    zs = z + steps\n",
    "    if steps > 0:\n",
    "        factor = 1 / pochhammer(z, steps)\n",
    "    elif steps < 0:\n",
    "        factor = pochhammer(zs, -steps)\n",
    "    return zs, factor\n",
    "\n",
    "\n",
    "def laguerre_gamma(z, x, w, target=11):\n",
    "    # res = 0.0\n",
    "    z = complex(z)\n",
    "    # if z.real < 1e-3:\n",
    "    #     res = pi / (\n",
    "    #         sin(pi * z) * laguerre_gamma(1 - z, x, w, target)\n",
    "    #     )  # Reflection formula\n",
    "    # else:\n",
    "        # z_shifted, correction_factor = find_shift(z, target)\n",
    "        # res = np.sum(x ** (z_shifted - 1) * w)\n",
    "        # res *= correction_factor\n",
    "    z_shifted, correction_factor = find_shift(z, target)\n",
    "    res = np.sum(x ** (z_shifted - 1) * w)\n",
    "    res *= correction_factor\n",
    "    res = drop_imag(res)\n",
    "    return res\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def eval_laguerre(x, target=12):\n",
    "    return np.array([laguerre_gamma(xi, zeros, weights, target) for xi in x])\n",
    "\n",
    "\n",
    "def eval_lanczos(x):\n",
    "    return np.array([lanczos_gamma(xi) for xi in x])\n",
    "\n",
    "\n",
    "def eval_mean_laguerre(x, targets):\n",
    "    return np.mean([eval_laguerre(x, target) for target in targets], 0)\n",
    "\n",
    "\n",
    "def calc_rel_error(x, y):\n",
    "    return (y - x) / x\n",
    "\n",
    "\n",
    "def evaluate(x, target=12):\n",
    "    lanczos_gammas = eval_lanczos(x)\n",
    "    laguerre_gammas = eval_laguerre(x, target)\n",
    "    rel_error = calc_rel_error(lanczos_gammas, laguerre_gammas)\n",
    "    return rel_error\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Test with real values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Empirische Tests zeigen:\n",
    "- $n=4 \\Rightarrow m=6$\n",
    "- $n=5 \\Rightarrow m=7$ oder $m=8$\n",
    "- $n=6 \\Rightarrow m=9$\n",
    "- $n=7 \\Rightarrow m=10$\n",
    "- $n=8 \\Rightarrow m=11$ oder $m=12$\n",
    "- $n=9 \\Rightarrow m=13$\n",
    "- $n=10 \\Rightarrow m=14$\n",
    "- $n=11 \\Rightarrow m=15$ oder $m=16$\n",
    "- $n=12 \\Rightarrow m=17$\n",
    "- $n=13 \\Rightarrow m=18 \\Rightarrow $ Beginnt numerisch instabil zu werden \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zeros, weights = np.polynomial.laguerre.laggauss(8)\n",
    "targets = np.arange(9, 14)\n",
    "mean_targets = ((9, 10),)\n",
    "x = np.linspace(EPSILON, 1 - EPSILON, 101)\n",
    "_, axs = plt.subplots(\n",
    "    2, sharex=True, clear=True, constrained_layout=True, figsize=(12, 12)\n",
    ")\n",
    "\n",
    "lanczos = eval_lanczos(x)\n",
    "for mean_target in mean_targets:\n",
    "    vals = eval_mean_laguerre(x, mean_target)\n",
    "    rel_error_mean = calc_rel_error(lanczos, vals)\n",
    "    axs[0].plot(x, rel_error_mean, label=mean_target)\n",
    "    axs[1].semilogy(x, np.abs(rel_error_mean), label=mean_target)\n",
    "\n",
    "mins = []\n",
    "maxs = []\n",
    "for target in targets:\n",
    "    rel_error = evaluate(x, target)\n",
    "    mins.append(np.min(np.abs(rel_error[(0.05 <= x) & (x <= 0.95)])))\n",
    "    maxs.append(np.max(np.abs(rel_error)))\n",
    "    axs[0].plot(x, rel_error, label=target)\n",
    "    axs[1].semilogy(x, np.abs(rel_error), label=target)\n",
    "# axs[0].set_ylim(*(np.array([-1, 1]) * 3.5e-8))\n",
    "\n",
    "axs[0].set_xlim(x[0], x[-1])\n",
    "axs[1].set_ylim(np.min(mins), 1.04*np.max(maxs))\n",
    "for ax in axs:\n",
    "    ax.legend()\n",
    "    ax.grid(which=\"both\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "targets = (16, 17)\n",
    "xmax = 15\n",
    "x = np.linspace(-xmax + EPSILON, xmax - EPSILON, 1000)\n",
    "\n",
    "mean_lag = eval_mean_laguerre(x, targets)\n",
    "lanczos = eval_lanczos(x)\n",
    "rel_error = calc_rel_error(lanczos, mean_lag)\n",
    "rel_error_simple = evaluate(x, targets[-1])\n",
    "# rel_error = evaluate(x, target)\n",
    "\n",
    "_, axs = plt.subplots(\n",
    "    2, sharex=True, clear=True, constrained_layout=True, figsize=(12, 12)\n",
    ")\n",
    "axs[0].plot(x, rel_error, label=targets)\n",
    "axs[1].semilogy(x, np.abs(rel_error), label=targets)\n",
    "axs[0].plot(x, rel_error_simple, label=targets[-1])\n",
    "axs[1].semilogy(x, np.abs(rel_error_simple), label=targets[-1])\n",
    "axs[0].set_xlim(x[0], x[-1])\n",
    "# axs[0].set_ylim(*(np.array([-1, 1]) * 4.2e-8))\n",
    "# axs[1].set_ylim(1e-10, 5e-8)\n",
    "for ax in axs:\n",
    "    ax.legend()\n",
    "\n",
    "x2 = np.linspace(-5 + EPSILON, 5, 4001)\n",
    "_, ax = plt.subplots(constrained_layout=True, figsize=(8, 6))\n",
    "ax.plot(x2, eval_mean_laguerre(x2, targets))\n",
    "ax.set_xlim(x2[0], x2[-1])\n",
    "ax.set_ylim(-7.5, 25)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Test with complex values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "targets = (16, 17)\n",
    "vals = np.linspace(-5 + EPSILON, 5, 100)\n",
    "x, y = np.meshgrid(vals, vals)\n",
    "mesh = x + 1j * y\n",
    "input = mesh.flatten()\n",
    "\n",
    "mean_lag = eval_mean_laguerre(input, targets).reshape(mesh.shape)\n",
    "lanczos = eval_lanczos(input).reshape(mesh.shape)\n",
    "rel_error = np.abs(calc_rel_error(lanczos, mean_lag))\n",
    "\n",
    "lag = eval_laguerre(input, targets[-1]).reshape(mesh.shape)\n",
    "rel_error_simple = np.abs(calc_rel_error(lanczos, lag))\n",
    "# rel_error = evaluate(x, target)\n",
    "\n",
    "fig, axs = plt.subplots(\n",
    "    2,\n",
    "    2,\n",
    "    sharex=True,\n",
    "    sharey=True,\n",
    "    clear=True,\n",
    "    constrained_layout=True,\n",
    "    figsize=(12, 10),\n",
    ")\n",
    "_c = axs[0, 1].pcolormesh(x, y, np.log10(np.abs(lanczos - mean_lag)), shading=\"gouraud\")\n",
    "_c = axs[0, 0].pcolormesh(x, y, np.log10(np.abs(lanczos - lag)), shading=\"gouraud\")\n",
    "fig.colorbar(_c, ax=axs[0, :])\n",
    "_c = axs[1, 1].pcolormesh(x, y, np.log10(rel_error), shading=\"gouraud\")\n",
    "_c = axs[1, 0].pcolormesh(x, y, np.log10(rel_error_simple), shading=\"gouraud\")\n",
    "fig.colorbar(_c, ax=axs[1, :])\n",
    "_ = axs[0, 0].set_title(\"Absolute Error\")\n",
    "_ = axs[1, 0].set_title(\"Relative Error\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z = 0.5\n",
    "ns = [4, 5, 5, 6, 7, 8, 8, 9, 10, 11, 11, 12]  # np.arange(4, 13)\n",
    "ms = np.arange(6, 18)\n",
    "xi = np.logspace(0, 2, 201)[:, None]\n",
    "lanczos = eval_lanczos([z])[0]\n",
    "\n",
    "_, ax = plt.subplots(clear=True, constrained_layout=True, figsize=(12, 8))\n",
    "ax.grid(1)\n",
    "for n, m in zip(ns, ms):\n",
    "    zeros, weights = np.polynomial.laguerre.laggauss(n)\n",
    "    c = scipy.special.factorial(n) ** 2 / scipy.special.factorial(2 * n)\n",
    "    e = np.abs(\n",
    "        scipy.special.poch(z - 2 * n, 2 * n)\n",
    "        / scipy.special.poch(z - m, m)\n",
    "        * c\n",
    "        * xi ** (z - 2 * n + m - 1)\n",
    "    )\n",
    "    ez = np.sum(\n",
    "        scipy.special.poch(z - 2 * n, 2 * n)\n",
    "        / scipy.special.poch(z - m, m)\n",
    "        * c\n",
    "        * zeros[:, None] ** (z - 2 * n + m - 1),\n",
    "        0,\n",
    "    )\n",
    "    lag = eval_laguerre([z], m)[0]\n",
    "    err = np.abs(lanczos - lag)\n",
    "    # print(m+z,ez)\n",
    "    # for zi,ezi in zip(z[0], ez):\n",
    "    #     print(f\"{m+zi}: {ezi}\")\n",
    "    # ax.semilogy(xi, e, color=color)\n",
    "    lines = ax.loglog(xi, e, label=str(n))\n",
    "    ax.axhline(err, color=lines[0].get_color())\n",
    "    # ax.set_xticks(np.arange(xi[-1] + 1))\n",
    "    # ax.set_ylim(1e-8, 1e5)\n",
    "_ = ax.legend()\n",
    "# _ = ax.legend([f\"z={zi}\" for zi in z[0]])\n",
    "# _ = [ax.axvline(x) for x in zeros]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bests = []\n",
    "N = 200\n",
    "step = 1 / (N - 1)\n",
    "a = 11 / 8\n",
    "b = 1 / 2\n",
    "x = np.linspace(step, 1 - step, N + 1)\n",
    "ns = np.arange(2, 13)\n",
    "for n in ns:\n",
    "    zeros, weights = np.polynomial.laguerre.laggauss(n)\n",
    "    est = np.ceil(b + a * n)\n",
    "    targets = np.arange(max(est - 2, 0), est + 3)\n",
    "    rel_errors = np.stack([np.abs(evaluate(x, target)) for target in targets], -1)\n",
    "    best = np.argmin(rel_errors, -1) + targets[0]\n",
    "    bests.append(best)\n",
    "bests = np.stack(bests, 0)\n",
    "\n",
    "fig, ax = plt.subplots(clear=True, constrained_layout=True, figsize=(5, 3))\n",
    "v = ax.imshow(bests, cmap=\"inferno\", aspect=\"auto\")\n",
    "plt.colorbar(v, ax=ax, label=r'$m$')\n",
    "ticks = np.arange(0, N + 1, 10)\n",
    "ax.set_xlim(0, 1)\n",
    "ax.set_xticks(ticks, [f\"{v:.2f}\" for v in ticks / N])\n",
    "ax.set_xticks(np.arange(N + 1), minor=True)\n",
    "ax.set_yticks(np.arange(len(ns)), ns)\n",
    "ax.set_xlabel(r\"$z$\")\n",
    "ax.set_ylabel(r\"$n$\")\n",
    "# for best in bests:\n",
    "#     print(\", \".join([f\"{int(b):2d}\" for b in best]))\n",
    "# print(np.unique(bests, return_counts=True))\n",
    "\n",
    "targets = np.mean(bests, -1)\n",
    "intercept, bias = np.polyfit(ns, targets, 1)\n",
    "_, axs2 = plt.subplots(2, sharex=True, clear=True, constrained_layout=True)\n",
    "xl = np.array([1, ns[-1] + 1])\n",
    "axs2[0].plot(ns, intercept * ns + bias)\n",
    "axs2[0].plot(ns, targets, \"x\")\n",
    "axs2[1].plot(ns, ((intercept * ns + bias) - targets), \"-x\")\n",
    "print(np.mean(bests, -1))\n",
    "print(f\"Intercept={intercept:.6g}, Bias={bias:.6g}\")\n",
    "\n",
    "\n",
    "predicts = np.ceil(intercept * ns[:, None] + bias - x)\n",
    "print(np.sum(np.abs(bests-predicts)))\n",
    "# for best in predicts:\n",
    "#     print(\", \".join([f\"{int(b):2d}\" for b in best]))\n"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "767d51c1340bd893661ea55ea3124f6de3c7a262a8b4abca0554b478b1e2ff90"
  },
  "kernelspec": {
   "display_name": "Python 3.8.10 64-bit",
   "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.8.10"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}