Unraveling of a bad factorial implementation

bad_factorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e0166afc-86c1-48b0-8d0d-547e11fc6df3",
   "metadata": {},
   "source": [
    "# Original code\n",
    "\n",
    "Taken from [this Reddit post](https://www.reddit.com/r/programminghorror/comments/1g6ea5q/god_i_love_python/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "9c7b3653-3108-4903-8549-6ba319b1ea2e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3628800"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def Factorial(val=-1, data=[]):\n",
    "    if val > 0:\n",
    "        data.clear()\n",
    "        data.append(val)\n",
    "        return Factorial()\n",
    "    pain = data.pop()\n",
    "    if pain == 1:\n",
    "        return pain\n",
    "    data.append(pain-1)\n",
    "    return Factorial()*pain\n",
    "\n",
    "\n",
    "Factorial(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7283b032-1870-4e98-bd24-dd91034fba9c",
   "metadata": {},
   "source": [
    "# Unraveling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a7aaabb8-d18b-4d4e-ab56-d28a40ec337d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3628800"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "class FactorialCalculator:\n",
    "    '''Compute the factorial of n recursively by storing the parameter\n",
    "    to the recursive application in a register variable and calling a\n",
    "    recurse method without parameters.\n",
    "\n",
    "    Silly, but this is how the algorithm above operates.\n",
    "    '''\n",
    "\n",
    "    # We don't need a list here, only one optional element\n",
    "    # The original code abused the fact that Python stores\n",
    "    # default-parameter lists by reference so the default\n",
    "    # list acted as a register\n",
    "    register: int | None\n",
    "\n",
    "    def __init__(self) -> None:\n",
    "        self.register = None\n",
    "\n",
    "    def recurse(self) -> int:\n",
    "        n = self.register\n",
    "\n",
    "        if n is None or n < 2:  # May be 0, 1 or None; the original only handles 1\n",
    "            return 1\n",
    "\n",
    "        self.register -= 1\n",
    "        return n * self.recurse()\n",
    "\n",
    "    def calculate(self, n: int) -> int:\n",
    "        if n < 0:\n",
    "            raise ValueError('What do you expect?')\n",
    "\n",
    "        self.register = n\n",
    "        return self.recurse()\n",
    "\n",
    "\n",
    "FactorialCalculator().calculate(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "abf7f9ee-0d74-4658-ae6b-a726ab2832ee",
   "metadata": {},
   "source": [
    "# Back to normal recursion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "bd244714-bdec-4a93-ba0d-430d5ba7e595",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3628800"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def factorial(n: int) -> int:\n",
    "    '''Look, ma, no need for registers. So no need for persistent storage,\n",
    "    which reduces the algorithm above to the well-known single-function implementation.'''\n",
    "\n",
    "    if n < 0:\n",
    "        raise ValueError('What do you expect?')\n",
    "\n",
    "    if n < 2:\n",
    "        return 1\n",
    "\n",
    "    return n * factorial(n - 1)\n",
    "\n",
    "\n",
    "factorial(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96ddea15-8706-420a-a827-8176b449a5f6",
   "metadata": {},
   "source": [
    "# __Bonus__: No recursion\n",
    "\n",
    "Although factorials are commonly presented as a recursive problem that can be optimized in various ways, there is a simple numerical solution (modulo the inherent complexity of numerical mathematics, especially with floating-point numbers involved)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "c76f5d22-5fc3-45ef-b332-574c29206f12",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3628800"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import math\n",
    "\n",
    "\n",
    "def factorial_via_gamma(n: int) -> int:\n",
    "    if n < 0:\n",
    "        raise ValueError('What do you expect?')\n",
    "\n",
    "    return int(math.gamma(n + 1))\n",
    "\n",
    "\n",
    "factorial_via_gamma(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9bd390ad-91f3-4eb4-8c74-5c99993912a8",
   "metadata": {},
   "source": [
    "If the above looks like cheating, here is a naive Monte-Carlo integration algorithm. For $10!$ it is much more computation-heavy, but the point is clear."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "eeb2e3bc-40d7-424b-80b4-7fb36c1c47b8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3628800"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.stats.qmc import Sobol\n",
    "\n",
    "\n",
    "def factorial_via_gamma_monte_carlo(n: int, iterations: int = 2 ** 15) -> int:\n",
    "    if n < 0:\n",
    "        raise ValueError('What do you expect?')\n",
    "\n",
    "    x = Sobol(1).random(iterations)\n",
    "    # We use the transform t ↦ 1 / (1 + t) on the conventional integrand\n",
    "    t = (1 - x) / x\n",
    "    y = np.pow(t, 10) * np.exp(-t) * np.pow(x, -2)\n",
    "\n",
    "    return int(y.mean())\n",
    "\n",
    "\n",
    "factorial_via_gamma_monte_carlo(10)"
   ]
  }
 ],
 "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.13.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}