{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "98e12805-2379-47fa-ada5-a53406f8210c",
   "metadata": {},
   "source": [
    "# Factorials and rounding errors"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2443558c-27fe-4bf5-8692-ff50a3ac7b24",
   "metadata": {},
   "source": [
    "Credit: Alex Gezerlis, \"Numerical Methods in Physics with Python\" (Cambridge University Press 2023)\n",
    "\n",
    "Consider the integrals   \n",
    "$f(n) = \\int_0^1 dx x^n e^{-x}$  \n",
    "so that  \n",
    "$f(0) = \\int_0^1 dx e^{-x} = 1-e^{-1}$.     \n",
    "You can find larger $n$ results doing partial integration,  \n",
    "$f(n) = \\int_0^1 dx x^n e^{-x} = -\\int_0^1 dx x^n \\frac{d e^{-x}}{dx} \n",
    "=  [-x^n e^{-x}]_0^1 + n \\int_0^1 dx x^{n-1} e^{-x}\n",
    "=  -e^{-1} + n \\int_0^1 dx x^{n-1} e^{-x}\n",
    "= -e^{-1} + f(n-1)$,    \n",
    "so we found a recursive relations  \n",
    "$f(n) = -e^{-1} + n*f(n-1)$,  \n",
    "and we know $f(0)$. Why not compute the integrals numerically:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "c2d9900e-9d42-4daa-8485-0eae8b9a5346",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 0.6321205588285577\n",
      "1 0.26424111765711533\n",
      "2 0.16060279414278833\n",
      "3 0.11392894125692266\n",
      "4 0.08783632385624829\n",
      "5 0.07130217810979911\n",
      "6 0.059933627487352314\n",
      "7 0.05165595124002387\n",
      "8 0.045368168748748605\n",
      "9 0.04043407756729511\n",
      "10 0.03646133450150879\n",
      "11 0.033195238345154365\n",
      "12 0.03046341897041005\n",
      "13 0.028145005443888316\n",
      "14 0.026150635042994086\n",
      "15 0.024380084473468955\n",
      "16 0.022201910404060943\n",
      "17 0.009553035697593693\n",
      "18 -0.19592479861475587\n",
      "19 -4.090450614851804\n",
      "20 -82.17689173820752\n"
     ]
    }
   ],
   "source": [
    "from numpy import exp\n",
    "def f(n):\n",
    "    if n==0:\n",
    "        return 1-exp(-1)\n",
    "    else:\n",
    "        return -exp(-1) + n*f(n-1)\n",
    "\n",
    "# test\n",
    "up_recursion_results={}\n",
    "for n in range(21):\n",
    "    fn = f(n)\n",
    "    print(n,fn)\n",
    "    up_recursion_results[n] = fn\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d953640-b034-456d-bec1-818c1434990e",
   "metadata": {},
   "source": [
    "<font color=\"red\"> Started fine for small $n$, but $n=20$ can't be correct! </font>    \n",
    "The problem is that error accumulates during iteration. The value of $f(0)$ is correct up to machine precision,\n",
    "which is about $10^{-16}$, double precision accuracy commonly used. You can see it by trying "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "3a237fd7-b88b-44d6-bd8e-3c92c752b76f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "False\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "x = 1.0\n",
    "print(x + 1e-15 == x)\n",
    "print(x + 1e-16 == x)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99f885f8-717f-469c-ad3d-0fec0942a3f3",
   "metadata": {},
   "source": [
    "The recursion adds to $-\\exp(-1) \\approx -0.36787944117144233$ values that are\n",
    "multiplied by $n$.  \n",
    "This is how error propagates:  \n",
    "$f(0) = (1-0.36787944117144233) \\pm 10^{-16}  = 0.6321205588285577 \\pm 10^{-16}$\n",
    "$f(1) =  0.26424111765711533 \\pm 10^{-16} + 1*0.6321205588285577 \\pm (1*10^{-16})$   \n",
    "$f(2) = -0.36787944117144233 \\pm 10^{-16} + 2*0.26424111765711533 \\pm (1*2*10^{-16})$  \n",
    "$...$   \n",
    "$f(20) = $ correct value $\\pm (1*2*3*...*20)*10^{-16}$  \n",
    "$~~~~~~~= $ 0.0183504676972562 $\\pm 20!*10^{-16}$  \n",
    "$~~~~~~~= $ 0.0183504676972562 $\\pm 243.290200817664$.  \n",
    "The error is inacceptable. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f2ad280-c46f-432d-b351-0fae189cf8cd",
   "metadata": {},
   "source": [
    "Try solving the relation this way:  \n",
    "$f(n) = -e^{-1} + n*f(n-1) \\Leftrightarrow  f(n-1) = \\frac{e^{-1} + f(n)}{n}$.   \n",
    "This is a reverse recursion, starting from some large-$n$ value, say, $f(40)$, and working down.  \n",
    "Now you say that \"but I don't have any idea what $f(40)$ is!\"  \n",
    "The beauty is that you don't *need* to know it accurately!   \n",
    "If you guess it to be  \n",
    "$f(40) = 0.1$,   \n",
    "you actually guess it to be   \n",
    "$f(40) = 0.1 \\pm error$.  \n",
    "The recursion downward divides the $error$ first by $40$, then by $40*39$ and so on.  \n",
    "Once you reach $f(20)$, your final error is $error/(40*39*38*...*21) = error/(40!/20!) \\approx error/(3.354*10^{29})$, which is negligible!  \n",
    "Code this to another recursion: \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "d5c3ff4a-579d-418e-bdd5-a5f9ad0e60b2",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "39 0.011696986029286057\n",
      "38 0.009732728902582779\n",
      "37 0.009937162370369082\n",
      "36 0.010211259555184092\n",
      "35 0.010502519464628511\n",
      "34 0.010810913161030595\n",
      "33 0.01113795159801391\n",
      "32 0.011485375538468371\n",
      "31 0.011855150522184709\n",
      "30 0.012249502957858937\n",
      "29 0.012670964804310042\n",
      "28 0.013122427792267324\n",
      "27 0.013607209605846772\n",
      "26 0.014129135213973671\n",
      "25 0.014692637553285232\n",
      "24 0.015302883148989102\n",
      "23 0.015965930180017976\n",
      "22 0.016688929189193926\n",
      "21 0.01748038047093801\n",
      "20 0.018350467697256206\n",
      "19 0.01931149544343493\n",
      "18 0.020378470348151434\n",
      "17 0.021569883973310763\n",
      "16 0.0229087838320443\n",
      "15 0.024424264062717915\n",
      "14 0.026153580348944015\n",
      "13 0.028145215822884737\n",
      "12 0.030463435153409775\n",
      "11 0.03319523969373767\n",
      "10 0.03646133462410727\n",
      "9 0.04043407757955496\n",
      "8 0.045368168750110814\n",
      "7 0.05165595124019414\n",
      "6 0.059933627487376635\n",
      "5 0.07130217810980316\n",
      "4 0.0878363238562491\n",
      "3 0.11392894125692285\n",
      "2 0.1606027941427884\n",
      "1 0.2642411176571154\n",
      "0 0.6321205588285577\n"
     ]
    }
   ],
   "source": [
    "def f_down(n):\n",
    "    if n==40:\n",
    "        return 0.1\n",
    "    else:\n",
    "        return (exp(-1) + f_down(n+1))/(n+1)\n",
    "\n",
    "# test \n",
    "for n in reversed(range(40)):\n",
    "    print(n, f_down(n))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b19d545-ee55-4e33-87f8-2526ee3d13fe",
   "metadata": {},
   "source": [
    "The accurate value was $f(20) = 0.0183504676972562$, so the deviation is just   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "51fca24d-b4dd-4c04-99c0-620718630258",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.018350467697256206\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "-6.938893903907228e-18"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print(f_down(20))\n",
    "0.0183504676972562 - f_down(20)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "53613fea-ba6e-4300-a70e-e9bc5a6dc5fe",
   "metadata": {},
   "source": [
    "Change the guessed value of $f(40)$, and you get a surprising result!"
   ]
  }
 ],
 "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.12.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}