Merge branch 'AndreasFMueller:master' into master

13 files changed, 622 insertions, 495 deletions

diff --git a/buch/papers/laguerre/scripts/estimates.py b/buch/papers/laguerre/scripts/estimates.py
new file mode 100644
index 0000000..1acd7f7
--- /dev/null
+++ b/buch/papers/laguerre/scripts/estimates.py
@@ -0,0 +1,49 @@
+if __name__ == "__main__":
+    import matplotlib as mpl
+    import matplotlib.pyplot as plt
+    import numpy as np
+
+    import gamma_approx as ga
+    import targets
+
+    mpl.rcParams.update(
+        {
+            "mathtext.fontset": "stix",
+            "font.family": "serif",
+            "font.serif": "TeX Gyre Termes",
+        }
+    )
+    
+    N = 200
+    ns = np.arange(1, 13)
+    step = 1 / (N - 1)
+    x = np.linspace(step, 1 - step, N + 1)
+
+    bests = targets.find_best_loc(N, ns=ns)
+    mean_m = np.mean(bests, -1)
+
+    intercept, bias = np.polyfit(ns, mean_m, 1)
+    fig, axs = plt.subplots(
+        2, num=1, sharex=True, clear=True, constrained_layout=True, figsize=(4.5, 3.6)
+    )
+    xl = np.array([ns[0] - 0.5, ns[-1] + 0.5])
+    axs[0].plot(xl, intercept * xl + bias, label=r"$\hat{m}$")
+    axs[0].plot(ns, mean_m, "x", label=r"$\overline{m}$")
+    axs[1].plot(
+        ns, ((intercept * ns + bias) - mean_m), "-x", label=r"$\hat{m} - \overline{m}$"
+    )
+    axs[0].set_xlim(*xl)
+    axs[0].set_xticks(ns)
+    axs[0].set_yticks(np.arange(np.floor(mean_m[0]), np.ceil(mean_m[-1]) + 0.1, 2))
+    # axs[0].set_title("Schätzung von Mittelwert")
+    # axs[1].set_title("Fehler")
+    axs[-1].set_xlabel(r"$n$")
+    for ax in axs:
+        ax.grid(1)
+        ax.legend()
+    # fig.savefig(f"{ga.img_path}/estimates.pgf")
+    fig.savefig(f"{ga.img_path}/estimates.pdf")
+
+    print(f"Intercept={intercept:.6g}, Bias={bias:.6g}")
+    predicts = np.ceil(intercept * ns[:, None] + bias - np.real(x))
+    print(f"Error: {np.mean(np.abs(bests - predicts))}")
diff --git a/buch/papers/laguerre/scripts/gamma_approx.ipynb b/buch/papers/laguerre/scripts/gamma_approx.ipynb
deleted file mode 100644
index 44f3abd..0000000
--- a/buch/papers/laguerre/scripts/gamma_approx.ipynb
+++ /dev/null
@@ -1,395 +0,0 @@
-{
- "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",
-    "    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(12)\n",
-    "targets = np.arange(16, 21)\n",
-    "mean_targets = ((16, 17),)\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.1 <= x) & (x <= 0.9)])))\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"
-   ]
-  }
- ],
- "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
-}
diff --git a/buch/papers/laguerre/scripts/gamma_approx.py b/buch/papers/laguerre/scripts/gamma_approx.py
new file mode 100644
index 0000000..5b09e59
--- /dev/null
+++ b/buch/papers/laguerre/scripts/gamma_approx.py
@@ -0,0 +1,116 @@
+from pathlib import Path
+
+import numpy as np
+import scipy.special
+
+EPSILON = 1e-7
+root = str(Path(__file__).parent)
+img_path = f"{root}/../images"
+fontsize = "medium"
+cmap = "plasma"
+
+
+def _prep_zeros_and_weights(x, w, n):
+    if x is None or w is None:
+        return np.polynomial.laguerre.laggauss(n)
+    return x, w
+
+
+def drop_imag(z):
+    if abs(z.imag) <= EPSILON:
+        z = z.real
+    return z
+
+
+def pochhammer(z, n):
+    return np.prod(z + np.arange(n))
+
+
+def find_optimal_shift(z, n):
+    mhat = 1.34093 * n + 0.854093
+    steps = int(np.floor(mhat - np.real(z)))
+    return steps
+
+
+def get_shifting_factor(z, steps):
+    factor = 1.0
+    if steps > 0:
+        factor = 1 / pochhammer(z, steps)
+    elif steps < 0:
+        factor = pochhammer(z + steps, -steps)
+    return factor
+
+
+def laguerre_gamma_shifted(z, x=None, w=None, n=8, target=11):
+    x, w = _prep_zeros_and_weights(x, w, n)
+    n = len(x)
+
+    z += 0j
+    steps = int(np.floor(target - np.real(z)))
+    z_shifted = z + steps
+    correction_factor = get_shifting_factor(z, steps)
+
+    res = np.sum(x ** (z_shifted - 1) * w)
+    res *= correction_factor
+    res = drop_imag(res)
+    return res
+
+
+def laguerre_gamma_opt_shifted(z, x=None, w=None, n=8):
+    if z == 0.0:
+        return np.infty
+    x, w = _prep_zeros_and_weights(x, w, n)
+    n = len(x)
+
+    z += 0j
+    opt_shift = find_optimal_shift(z, n)
+    correction_factor = get_shifting_factor(z, opt_shift)
+    z_shifted = z + opt_shift
+
+    res = np.sum(x ** (z_shifted - 1) * w)
+    res *= correction_factor
+    res = drop_imag(res)
+    return res
+
+
+def laguerre_gamma_simple(z, x=None, w=None, n=8):
+    if z == 0.0:
+        return np.infty
+    x, w = _prep_zeros_and_weights(x, w, n)
+    z += 0j
+    res = np.sum(x ** (z - 1) * w)
+    res = drop_imag(res)
+    return res
+
+
+def laguerre_gamma_mirror(z, x=None, w=None, n=8):
+    if z == 0.0:
+        return np.infty
+    x, w = _prep_zeros_and_weights(x, w, n)
+    z += 0j
+    if z.real < 1e-3:
+        return np.pi / (
+            np.sin(np.pi * z) * laguerre_gamma_simple(1 - z, x, w)
+        )  # Reflection formula
+    return laguerre_gamma_simple(z, x, w)
+
+
+def eval_laguerre_gamma(z, x=None, w=None, n=8, func="simple", **kwargs):    
+    x, w = _prep_zeros_and_weights(x, w, n)    
+    if func == "simple":
+        f = laguerre_gamma_simple
+    elif func == "mirror":
+        f = laguerre_gamma_mirror
+    elif func == "optimal_shifted":
+        f = laguerre_gamma_opt_shifted
+    else:
+        f = laguerre_gamma_shifted
+    return np.array([f(zi, x, w, n, **kwargs) for zi in z])
+
+
+def calc_rel_error(x, y):
+    mask = np.abs(x) != np.infty
+    rel_error = np.zeros_like(y)
+    rel_error[mask] = (y[mask] - x[mask]) / x[mask]
+    rel_error[~mask] = 0.0
+    return rel_error
diff --git a/buch/papers/laguerre/scripts/integrand.py b/buch/papers/laguerre/scripts/integrand.py
new file mode 100644
index 0000000..e970721
--- /dev/null
+++ b/buch/papers/laguerre/scripts/integrand.py
@@ -0,0 +1,42 @@
+#!/usr/bin/env python3
+# -*- coding:utf-8 -*-
+"""Plot for integrand of gamma function with shifting terms."""
+
+if __name__ == "__main__":
+    import os
+    from pathlib import Path
+
+    import matplotlib as mpl
+    import matplotlib.pyplot as plt
+    import numpy as np
+
+    mpl.rcParams.update(
+        {
+            "mathtext.fontset": "stix",
+            "font.family": "serif",
+            "font.serif": "TeX Gyre Termes",
+        }
+    )
+    
+    EPSILON = 1e-12
+    xlims = np.array([-3, 3])
+
+    root = str(Path(__file__).parent)
+    img_path = f"{root}/../images"
+    os.makedirs(img_path, exist_ok=True)
+
+    t = np.logspace(*xlims, 1001)[:, None]
+
+    z = np.array([-4.5, -2, -1, -0.5, 0.0, 0.5, 1, 2, 4.5])
+    r = t ** z
+
+    fig, ax = plt.subplots(num=1, clear=True, constrained_layout=True, figsize=(4, 2.4))
+    ax.semilogx(t, r)
+    ax.set_xlim(*(10.0 ** xlims))
+    ax.set_ylim(1e-3, 40)
+    ax.set_xlabel(r"$x$")
+    # ax.set_ylabel(r"$x^z$")
+    ax.grid(1, "both")
+    labels = [f"$z={zi: 3.1f}$" for zi in np.squeeze(z)]
+    ax.legend(labels, ncol=2, loc="upper left", fontsize="small")
+    fig.savefig(f"{img_path}/integrand.pdf")
diff --git a/buch/papers/laguerre/scripts/integrand_exp.py b/buch/papers/laguerre/scripts/integrand_exp.py
new file mode 100644
index 0000000..e649b26
--- /dev/null
+++ b/buch/papers/laguerre/scripts/integrand_exp.py
@@ -0,0 +1,46 @@
+#!/usr/bin/env python3
+# -*- coding:utf-8 -*-
+"""Plot for integrand of gamma function with shifting terms."""
+
+if __name__ == "__main__":
+    import os
+    from pathlib import Path
+
+    import matplotlib as mpl
+    import matplotlib.pyplot as plt
+    import numpy as np
+    
+    mpl.rcParams.update(
+        {
+            "mathtext.fontset": "stix",
+            "font.family": "serif",
+            "font.serif": "TeX Gyre Termes",
+        }
+    )
+
+    EPSILON = 1e-12
+    xlims = np.array([-3, 3])
+
+    root = str(Path(__file__).parent)
+    img_path = f"{root}/../images"
+    os.makedirs(img_path, exist_ok=True)
+
+    t = np.logspace(*xlims, 1001)[:, None]
+
+    z = np.array([-1, -0.5, 0.0, 0.5, 1, 2, 3, 4, 4.5])
+    e = np.exp(-t)
+    r = t ** z * e
+
+    fig, ax = plt.subplots(num=2, clear=True, constrained_layout=True, figsize=(4, 2.4))
+    ax.semilogx(t, r)
+    # ax.plot(t,np.exp(-t))
+    ax.set_xlim(10 ** (-2), 20)
+    ax.set_ylim(1e-3, 10)
+    ax.set_xlabel(r"$x$")
+    # ax.set_ylabel(r"$x^z e^{-x}$")
+    ax.grid(1, "both")
+    labels = [f"$z={zi: 3.1f}$" for zi in np.squeeze(z)]
+    ax.legend(labels, ncol=2, loc="upper left", fontsize="small")
+    # fig.savefig(f"{img_path}/integrand_exp.pgf")
+    fig.savefig(f"{img_path}/integrand_exp.pdf")
+    # plt.show()
diff --git a/buch/papers/laguerre/scripts/laguerre_plot.py b/buch/papers/laguerre/scripts/laguerre_plot.py
deleted file mode 100644
index b9088d0..0000000
--- a/buch/papers/laguerre/scripts/laguerre_plot.py
+++ /dev/null
@@ -1,100 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding:utf-8 -*-
-"""Some plots for Laguerre Polynomials."""
-
-import os
-from pathlib import Path
-
-import matplotlib.pyplot as plt
-import numpy as np
-import scipy.special as ss
-
-
-def get_ticks(start, end, step=1):
-    ticks = np.arange(start, end, step)
-    return ticks[ticks != 0]
-
-
-N = 1000
-step = 5
-t = np.linspace(-1.05, 10.5, N)[:, None]
-root = str(Path(__file__).parent)
-img_path = f"{root}/../images"
-os.makedirs(img_path, exist_ok=True)
-
-
-# fig = plt.figure(num=1, clear=True, tight_layout=True, figsize=(5.5, 3.7))
-# ax = fig.add_subplot(axes_class=AxesZero)
-fig, ax = plt.subplots(num=1, clear=True, constrained_layout=True, figsize=(6, 4))
-for n in np.arange(0, 8):
-    k = np.arange(0, n + 1)[None]
-    L = np.sum((-1) ** k * ss.binom(n, k) / ss.factorial(k) * t ** k, -1)
-    ax.plot(t, L, label=f"n={n}")
-
-ax.set_xticks(get_ticks(int(t[0]), t[-1]), minor=True)
-ax.set_xticks(get_ticks(0, t[-1], step))
-ax.set_xlim(t[0], t[-1] + 0.1 * (t[1] - t[0]))
-ax.set_xlabel(r"$x$", x=1.0, labelpad=-10, rotation=0, fontsize="large")
-
-ylim = 13
-ax.set_yticks(np.arange(-ylim, ylim), minor=True)
-ax.set_yticks(np.arange(-step * (ylim // step), ylim, step))
-ax.set_ylim(-ylim, ylim)
-ax.set_ylabel(r"$y$", y=0.95, labelpad=-18, rotation=0, fontsize="large")
-
-ax.legend(ncol=2, loc=(0.125, 0.01), fontsize="large")
-
-# set the x-spine
-ax.spines[["left", "bottom"]].set_position("zero")
-ax.spines[["right", "top"]].set_visible(False)
-ax.xaxis.set_ticks_position("bottom")
-hlx = 0.4
-dx = t[-1, 0] - t[0, 0]
-dy = 2 * ylim
-hly = dy / dx * hlx
-dps = fig.dpi_scale_trans.inverted()
-bbox = ax.get_window_extent().transformed(dps)
-width, height = bbox.width, bbox.height
-
-# manual arrowhead width and length
-hw = 1.0 / 60.0 * dy
-hl = 1.0 / 30.0 * dx
-lw = 0.5  # axis line width
-ohg = 0.0  # arrow overhang
-
-# compute matching arrowhead length and width
-yhw = hw / dy * dx * height / width
-yhl = hl / dx * dy * width / height
-
-# draw x and y axis
-ax.arrow(
-    t[-1, 0] - hl,
-    0,
-    hl,
-    0.0,
-    fc="k",
-    ec="k",
-    lw=lw,
-    head_width=hw,
-    head_length=hl,
-    overhang=ohg,
-    length_includes_head=True,
-    clip_on=False,
-)
-
-ax.arrow(
-    0,
-    ylim - yhl,
-    0.0,
-    yhl,
-    fc="k",
-    ec="k",
-    lw=lw,
-    head_width=yhw,
-    head_length=yhl,
-    overhang=ohg,
-    length_includes_head=True,
-    clip_on=False,
-)
-
-fig.savefig(f"{img_path}/laguerre_polynomes.pdf")
diff --git a/buch/papers/laguerre/scripts/laguerre_poly.py b/buch/papers/laguerre/scripts/laguerre_poly.py
new file mode 100644
index 0000000..05db5d3
--- /dev/null
+++ b/buch/papers/laguerre/scripts/laguerre_poly.py
@@ -0,0 +1,108 @@
+import numpy as np
+
+
+def get_ticks(start, end, step=1):
+    ticks = np.arange(start, end, step)
+    return ticks[ticks != 0]
+
+
+if __name__ == "__main__":
+    import os
+    from pathlib import Path
+
+    import matplotlib as mpl
+    import matplotlib.pyplot as plt
+    import scipy.special as ss
+    
+    mpl.rcParams.update(
+        {
+            "mathtext.fontset": "stix",
+            "font.family": "serif",
+            "font.serif": "TeX Gyre Termes",
+        }
+    )
+
+    N = 1000
+    step = 5
+    t = np.linspace(-1.05, 10.5, N)[:, None]
+    root = str(Path(__file__).parent)
+    img_path = f"{root}/../images"
+    os.makedirs(img_path, exist_ok=True)
+
+    # fig = plt.figure(num=1, clear=True, tight_layout=True, figsize=(5.5, 3.7))
+    # ax = fig.add_subplot(axes_class=AxesZero)
+    fig, ax = plt.subplots(num=1, clear=True, constrained_layout=True, figsize=(6, 4))
+    for n in np.arange(0, 8):
+        k = np.arange(0, n + 1)[None]
+        L = np.sum((-1) ** k * ss.binom(n, k) / ss.factorial(k) * t ** k, -1)
+        ax.plot(t, L, label=f"$n={n}$")
+
+    ax.set_xticks(get_ticks(int(t[0]), t[-1]), minor=True)
+    ax.set_xticks(get_ticks(0, t[-1], step))
+    ax.set_xlim(t[0], t[-1] + 0.1 * (t[1] - t[0]))
+    ax.set_xlabel(r"$x$", x=1.0, labelpad=-10, rotation=0, fontsize="large")
+
+    ylim = 13
+    ax.set_yticks(get_ticks(-ylim, ylim), minor=True)
+    ax.set_yticks(get_ticks(-step * (ylim // step), ylim, step))
+    ax.set_ylim(-ylim, ylim)
+    ax.set_ylabel(r"$y$", y=0.95, labelpad=-14, rotation=0, fontsize="large")
+
+    ax.legend(ncol=2, loc=(0.125, 0.01), fontsize="large")
+
+    # set the x-spine
+    ax.spines[["left", "bottom"]].set_position("zero")
+    ax.spines[["right", "top"]].set_visible(False)
+    ax.xaxis.set_ticks_position("bottom")
+    hlx = 0.4
+    dx = t[-1, 0] - t[0, 0]
+    dy = 2 * ylim
+    hly = dy / dx * hlx
+    dps = fig.dpi_scale_trans.inverted()
+    bbox = ax.get_window_extent().transformed(dps)
+    width, height = bbox.width, bbox.height
+
+    # manual arrowhead width and length
+    hw = 1.0 / 60.0 * dy
+    hl = 1.0 / 30.0 * dx
+    lw = 0.5  # axis line width
+    ohg = 0.0  # arrow overhang
+
+    # compute matching arrowhead length and width
+    yhw = hw / dy * dx * height / width
+    yhl = hl / dx * dy * width / height
+
+    # draw x and y axis
+    ax.arrow(
+        t[-1, 0] - hl,
+        0,
+        hl,
+        0.0,
+        fc="k",
+        ec="k",
+        lw=lw,
+        head_width=hw,
+        head_length=hl,
+        overhang=ohg,
+        length_includes_head=True,
+        clip_on=False,
+    )
+
+    ax.arrow(
+        0,
+        ylim - yhl,
+        0.0,
+        yhl,
+        fc="k",
+        ec="k",
+        lw=lw,
+        head_width=yhw,
+        head_length=yhl,
+        overhang=ohg,
+        length_includes_head=True,
+        clip_on=False,
+    )
+
+    # fig.savefig(f"{img_path}/laguerre_poly.pgf")
+    fig.savefig(f"{img_path}/laguerre_poly.pdf")
+    # plt.show()
diff --git a/buch/papers/laguerre/scripts/rel_error_complex.py b/buch/papers/laguerre/scripts/rel_error_complex.py
new file mode 100644
index 0000000..4a714fa
--- /dev/null
+++ b/buch/papers/laguerre/scripts/rel_error_complex.py
@@ -0,0 +1,43 @@
+if __name__ == "__main__":
+    import matplotlib as mpl
+    import matplotlib.pyplot as plt
+    import numpy as np
+    import scipy.special
+
+    import gamma_approx as ga
+
+    mpl.rcParams.update(
+        {
+            "mathtext.fontset": "stix",
+            "font.family": "serif",
+            "font.serif": "TeX Gyre Termes",
+        }
+    )
+    
+    xmax = 4
+    vals = np.linspace(-xmax + ga.EPSILON, xmax, 100)
+    x, y = np.meshgrid(vals, vals)
+    mesh = x + 1j * y
+    input = mesh.flatten()
+
+    lanczos = scipy.special.gamma(mesh)
+    lag = ga.eval_laguerre_gamma(input, n=8, func="optimal_shifted").reshape(mesh.shape)
+    rel_error = np.abs(ga.calc_rel_error(lanczos, lag))
+
+    fig, ax = plt.subplots(clear=True, constrained_layout=True, figsize=(3.5, 2.1))
+    _c = ax.pcolormesh(
+        x, y, rel_error, shading="gouraud", cmap=ga.cmap, norm=mpl.colors.LogNorm()
+    )
+    cbr = fig.colorbar(_c, ax=ax)
+    cbr.minorticks_off()
+    # ax.set_title("Relative Error")
+    ax.set_xlabel("Re($z$)")
+    ax.set_ylabel("Im($z$)")
+    minor_ticks = np.arange(-xmax, xmax + ga.EPSILON)
+    ticks = np.arange(-xmax, xmax + ga.EPSILON, 2)
+    ax.set_xticks(ticks)
+    ax.set_xticks(minor_ticks, minor=True)
+    ax.set_yticks(ticks)
+    ax.set_yticks(minor_ticks, minor=True)
+    fig.savefig(f"{ga.img_path}/rel_error_complex.pdf")
+    # plt.show()
diff --git a/buch/papers/laguerre/scripts/rel_error_mirror.py b/buch/papers/laguerre/scripts/rel_error_mirror.py
new file mode 100644
index 0000000..7348d5e
--- /dev/null
+++ b/buch/papers/laguerre/scripts/rel_error_mirror.py
@@ -0,0 +1,38 @@
+if __name__ == "__main__":
+    import matplotlib as mpl
+    import matplotlib.pyplot as plt
+    import numpy as np
+    import scipy.special
+
+    import gamma_approx as ga
+    
+    mpl.rcParams.update(
+        {
+            "mathtext.fontset": "stix",
+            "font.family": "serif",
+            "font.serif": "TeX Gyre Termes",
+        }
+    )
+
+    xmin = -15
+    xmax = 15
+    ns = np.arange(2, 12, 2)
+    ylim = np.array([-11, 1])
+    x = np.linspace(xmin + ga.EPSILON, xmax - ga.EPSILON, 400)
+    gamma = scipy.special.gamma(x)
+    fig, ax = plt.subplots(num=1, clear=True, constrained_layout=True, figsize=(5, 2.5))
+    for n in ns:
+        gamma_lag = ga.eval_laguerre_gamma(x, n=n, func="mirror")
+        rel_err = ga.calc_rel_error(gamma, gamma_lag)
+        ax.semilogy(x, np.abs(rel_err), label=f"$n={n}$")
+    ax.set_xlim(x[0], x[-1])
+    ax.set_ylim(*(10.0 ** ylim))
+    ax.set_xticks(np.arange(xmin, xmax + ga.EPSILON, 5))
+    ax.set_xticks(np.arange(xmin, xmax), minor=True)
+    ax.set_yticks(10.0 ** np.arange(*ylim, 2))
+    ax.set_xlabel(r"$z$")
+    # ax.set_ylabel("Relativer Fehler")
+    ax.legend(ncol=1, loc="upper left", fontsize=ga.fontsize)
+    ax.grid(1, "both")
+    # fig.savefig(f"{ga.img_path}/rel_error_mirror.pgf")
+    fig.savefig(f"{ga.img_path}/rel_error_mirror.pdf")
diff --git a/buch/papers/laguerre/scripts/rel_error_range.py b/buch/papers/laguerre/scripts/rel_error_range.py
new file mode 100644
index 0000000..ece3b6d
--- /dev/null
+++ b/buch/papers/laguerre/scripts/rel_error_range.py
@@ -0,0 +1,41 @@
+if __name__ == "__main__":
+    import matplotlib as mpl
+    import matplotlib.pyplot as plt
+    import numpy as np
+    import scipy.special
+
+    import gamma_approx as ga
+    
+    mpl.rcParams.update(
+        {
+            "mathtext.fontset": "stix",
+            "font.family": "serif",
+            "font.serif": "TeX Gyre Termes",
+        }
+    )
+    N = 1201
+    xmax = 6
+    xmin = -xmax
+    ns = np.arange(2, 12, 2)
+    ylim = np.array([-11, -1.2])
+
+    x = np.linspace(xmin + ga.EPSILON, xmax - ga.EPSILON, N)
+    gamma = scipy.special.gamma(x)
+    fig, ax = plt.subplots(num=1, clear=True, constrained_layout=True, figsize=(5, 2))
+    for n in ns:
+        gamma_lag = ga.eval_laguerre_gamma(x, n=n, func="optimal_shifted")
+        rel_err = ga.calc_rel_error(gamma, gamma_lag)
+        ax.semilogy(x, np.abs(rel_err), label=f"$n={n}$")
+    ax.set_xlim(x[0], x[-1])
+    ax.set_ylim(*(10.0 ** ylim))
+    ax.set_xticks(np.arange(xmin, xmax + ga.EPSILON, 2))
+    ax.set_xticks(np.arange(xmin, xmax + ga.EPSILON), minor=True)
+    ax.set_yticks(10.0 ** np.arange(*ylim, 2))
+    ax.set_yticks(10.0 ** np.arange(*ylim, 1), "", minor=True)
+    ax.set_xlabel(r"$z$")
+    # ax.set_ylabel("Relativer Fehler")
+    ax.legend(ncol=1, loc="upper left", fontsize=ga.fontsize)
+    ax.grid(1, "both")
+    # fig.savefig(f"{ga.img_path}/rel_error_range.pgf")
+    fig.savefig(f"{ga.img_path}/rel_error_range.pdf")
+    # plt.show()
diff --git a/buch/papers/laguerre/scripts/rel_error_shifted.py b/buch/papers/laguerre/scripts/rel_error_shifted.py
new file mode 100644
index 0000000..f53c89b
--- /dev/null
+++ b/buch/papers/laguerre/scripts/rel_error_shifted.py
@@ -0,0 +1,40 @@
+if __name__ == "__main__":
+    import matplotlib as mpl
+    import matplotlib.pyplot as plt
+    import numpy as np
+    import scipy.special
+
+    import gamma_approx as ga
+
+    mpl.rcParams.update(
+        {
+            "mathtext.fontset": "stix",
+            "font.family": "serif",
+            "font.serif": "TeX Gyre Termes",
+        }
+    )
+    n = 8  # order of Laguerre polynomial
+    N = 200  # number of points in interval
+
+    step = 1 / (N - 1)
+    x = np.linspace(step, 1 - step, N + 1)
+    targets = np.arange(10, 14)
+    gamma = scipy.special.gamma(x)
+    fig, ax = plt.subplots(num=1, clear=True, constrained_layout=True, figsize=(5, 2.1))
+    for target in targets:
+        gamma_lag = ga.eval_laguerre_gamma(x, target=target, n=n, func="shifted")
+        rel_error = np.abs(ga.calc_rel_error(gamma, gamma_lag))
+        ax.semilogy(x, rel_error, label=f"$m={target}$", linewidth=3)
+    gamma_lgo = ga.eval_laguerre_gamma(x, n=n, func="optimal_shifted")
+    rel_error = np.abs(ga.calc_rel_error(gamma, gamma_lgo))
+    ax.semilogy(x, rel_error, "m", linestyle=":", label="$m^*$", linewidth=3)
+    ax.set_xlim(x[0], x[-1])
+    ax.set_ylim(5e-9, 5e-8)
+    ax.set_xlabel(r"$z$")
+    ax.set_xticks(np.linspace(0, 1, 6))
+    ax.set_xticks(np.linspace(0, 1, 11), minor=True)
+    ax.grid(1, "both")
+    ax.legend(ncol=1, fontsize=ga.fontsize)
+    # fig.savefig(f"{ga.img_path}/rel_error_shifted.pgf")
+    fig.savefig(f"{ga.img_path}/rel_error_shifted.pdf")
+    # plt.show()
diff --git a/buch/papers/laguerre/scripts/rel_error_simple.py b/buch/papers/laguerre/scripts/rel_error_simple.py
new file mode 100644
index 0000000..e1ea36a
--- /dev/null
+++ b/buch/papers/laguerre/scripts/rel_error_simple.py
@@ -0,0 +1,41 @@
+if __name__ == "__main__":
+    import matplotlib as mpl
+    import matplotlib.pyplot as plt
+    import numpy as np
+    import scipy.special
+
+    import gamma_approx as ga
+
+    # mpl.rc("text", usetex=True)
+    mpl.rcParams.update(
+        {
+            "mathtext.fontset": "stix",
+            "font.family": "serif",
+            "font.serif": "TeX Gyre Termes",
+        }
+    )
+    # mpl.rcParams.update({"font.family": "serif", "font.serif": "TeX Gyre Termes"})
+
+    # Simple / naive
+    xmin = -5
+    xmax = 25
+    ns = np.arange(2, 12, 2)
+    ylim = np.array([-11, 6])
+    x = np.linspace(xmin + ga.EPSILON, xmax - ga.EPSILON, 400)
+    gamma = scipy.special.gamma(x)
+    fig, ax = plt.subplots(num=1, clear=True, constrained_layout=True, figsize=(5, 2.5))
+    for n in ns:
+        gamma_lag = ga.eval_laguerre_gamma(x, n=n)
+        rel_err = ga.calc_rel_error(gamma, gamma_lag)
+        ax.semilogy(x, np.abs(rel_err), label=f"$n={n}$")
+    ax.set_xlim(x[0], x[-1])
+    ax.set_ylim(*(10.0 ** ylim))
+    ax.set_xticks(np.arange(xmin, xmax + ga.EPSILON, 5))
+    ax.set_xticks(np.arange(xmin, xmax), minor=True)
+    ax.set_yticks(10.0 ** np.arange(*ylim, 2))
+    ax.set_xlabel(r"$z$")
+    # ax.set_ylabel("Relativer Fehler")
+    ax.legend(ncol=3, fontsize=ga.fontsize)
+    ax.grid(1, "both")
+    # fig.savefig(f"{ga.img_path}/rel_error_simple.pgf")
+    fig.savefig(f"{ga.img_path}/rel_error_simple.pdf")
diff --git a/buch/papers/laguerre/scripts/targets.py b/buch/papers/laguerre/scripts/targets.py
new file mode 100644
index 0000000..69f94ba
--- /dev/null
+++ b/buch/papers/laguerre/scripts/targets.py
@@ -0,0 +1,58 @@
+import numpy as np
+import scipy.special
+
+import gamma_approx as ga
+
+
+def find_best_loc(N=200, a=1.375, b=0.5, ns=None):
+    if ns is None:
+        ns = np.arange(2, 13)
+    bests = []
+    step = 1 / (N - 1)
+    x = np.linspace(step, 1 - step, N + 1)
+    gamma = scipy.special.gamma(x)
+    for n in ns:
+        zeros, weights = np.polynomial.laguerre.laggauss(n)
+        est = np.ceil(b + a * n)
+        targets = np.arange(max(est - 2, 0), est + 3)
+        rel_error = []
+        for target in targets:
+            gamma_lag = ga.eval_laguerre_gamma(x, target=target, x=zeros, w=weights, func="shifted")
+            rel_error.append(np.abs(ga.calc_rel_error(gamma, gamma_lag)))
+        rel_error = np.stack(rel_error, -1)
+        best = np.argmin(rel_error, -1) + targets[0]
+        bests.append(best)
+    return np.stack(bests, 0)
+
+
+if __name__ == "__main__":
+    import matplotlib as mpl
+    import matplotlib.pyplot as plt
+
+    mpl.rcParams.update(
+        {
+            "mathtext.fontset": "stix",
+            "font.family": "serif",
+            "font.serif": "TeX Gyre Termes",
+        }
+    )
+    
+    N = 200
+    ns = np.arange(1, 13)
+    
+    bests = find_best_loc(N, ns=ns)
+
+    fig, ax = plt.subplots(num=1, clear=True, constrained_layout=True, figsize=(3.5, 2.1))
+    v = ax.imshow(bests, cmap=ga.cmap, aspect="auto", interpolation="nearest")
+    plt.colorbar(v, ax=ax, label=r"$m^*$")
+    ticks = np.arange(0, N + 1, N // 5)
+    ax.set_xlim(0, 1)
+    ax.set_xticks(ticks)
+    ax.set_xticklabels([f"{v:.2f}" for v in ticks / N])
+    ax.set_xticks(np.arange(0, N + 1, N // 20), minor=True)
+    ax.set_yticks(np.arange(len(ns)))
+    ax.set_yticklabels(ns)
+    ax.set_xlabel(r"$z$")
+    ax.set_ylabel(r"$n$")
+    # fig.savefig(f"{ga.img_path}/targets.pgf")
+    fig.savefig(f"{ga.img_path}/targets.pdf")