{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\newcommand{\\stirling}[2]{\\begin{Bmatrix} #1\\\\#2\\end{Bmatrix}}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Surjections"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pour tous ensembles finis $E$ et $F$, notons $\\mathcal S(E,F)$ l'ensemble des surjections de $E$ sur $F$. Le cardinal de $\\mathcal S(E,F)$ ne dépend que des cardinaux $n$ et $p$ de $E$ et $F$. Notons-le $S^n_p$.\n",
    "\n",
    "Pour tous $n, p\\ge 1$,\n",
    "\n",
    "$$S^n_p=p\\left(S^{n-1}_{p - 1} + S^{n - 1}_p\\right)$$\n",
    "\n",
    "De plus, \n",
    "\n",
    "- $S^0_0=1$\n",
    "- pour tout $p\\ge 1$, $S^0_p=0$\n",
    "- pour tout $n\\ge 1$, $S^n_0=0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def matrice(p, q):\n",
    "    A = p * [None]\n",
    "    for i in range(p):\n",
    "        A[i] = q * [0]\n",
    "    return A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `surjections` renvoie la matrice des $S^n_p$ pour $0\\le n, p\\le N$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def surjections(N):\n",
    "    s = matrice(N + 1, N + 1)\n",
    "    s[0][0] = 1\n",
    "    for n in range(1, N + 1):\n",
    "        for p in range(1, n + 1):\n",
    "            s[n][p] = p * (s[n - 1][p - 1] + s[n - 1][p])\n",
    "    return s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Complexité de `surjections` : $O(N^2)$ opérations arithmétiques."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 6\n",
    "surj = surjections(N)\n",
    "for n in range(N + 1):\n",
    "    for p in range(N + 1):\n",
    "        print('%10d' % surj[n][p], end='')\n",
    "    print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "surj = surjections(456)\n",
    "surj[456][123]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(str(surj[456][123]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Partitions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 Dénombrements de partitions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soit $\\stirling n p$ le nombre de partitions d'un ensemble de cardinal $n$ en $p$ sous ensembles.\n",
    "\n",
    "Pour tous $n, p\\ge 1$,\n",
    "\n",
    "$$\\stirling n p=\\stirling{n-1}{p - 1} + p\\stirling{n - 1}p$$\n",
    "\n",
    "De plus,\n",
    "\n",
    "- $\\stirling0 0=1$\n",
    "- pour tout $p\\ge 1$, $\\stirling0 p=0$\n",
    "- pour tout $n\\ge 1$, $\\stirling n 0=0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def partitions(N):\n",
    "    s = matrice(N + 1, N + 1)\n",
    "    s[0][0] = 1\n",
    "    for n in range(1, N + 1):\n",
    "        for p in range(1, n + 1):\n",
    "            s[n][p] = s[n - 1][p - 1] + p * s[n - 1][p]\n",
    "    return s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Complexité de `partitions` : $O(N^2)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 7\n",
    "surj = partitions(N)\n",
    "for n in range(N + 1):\n",
    "    for p in range(N + 1):\n",
    "        print('%10d' % surj[n][p], end='')\n",
    "    print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 939\n",
    "p = 393\n",
    "partitions(n)[n][p]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_ % (10 ** 9)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Nombres de Bell"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Le $n$ième nombre de Bell est le nombre de partitions d'un ensemble de cardinal $n$.\n",
    "\n",
    "$$B_n=\\sum_{p=0}^n\\stirling n p$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bell(n):\n",
    "    s = partitions(n)[-1]\n",
    "    b = 0\n",
    "    for p in range(n + 1):\n",
    "        b += s[p]\n",
    "    return b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Complexité de `bell` : $O(n^2)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for n in range(15):\n",
    "    print(n, bell(n))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(bell(1000))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(str(bell(933)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Un autre calcul de $B_n$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On a $B_0=1$ et pour tout $n\\in\\mathbb N$,\n",
    "\n",
    "$$B_{n+1}=\\sum_{k=0}^n\\binom n kB_k$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `bell2` renvoie la liste $[B_0,\\ldots,B_N]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bell2(N):\n",
    "    cbs = [1]\n",
    "    bs = (N + 1) * [0]\n",
    "    bs[0] = 1\n",
    "    for n in range(1, N + 1):\n",
    "        for k in range(n):\n",
    "            bs[n] += cbs[k] * bs[k]\n",
    "        cbs1 = (n + 1) * [0]\n",
    "        cbs1[0] = 1\n",
    "        for k in range(1, n):\n",
    "            cbs1[k] = cbs[k - 1] + cbs[k]\n",
    "        cbs1[n] = 1\n",
    "        cbs = cbs1\n",
    "    return bs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Complexité de `bell2` : $O(N^2)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bell2(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bell2(1000)[1000]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Une valeur approchée de $B_n$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La *Formule de Dobinski* donne $B_n$ comme somme d'une série.\n",
    "\n",
    "$$B_n=\\frac 1 e \\sum_{k=0}^\\infty \\frac{k^n}{k!}$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fact(n):\n",
    "    p = 1\n",
    "    for k in range(1, n + 1):\n",
    "        p *= k\n",
    "    return p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `approx_bell` renvoie la $N$ème somme partielle de la série."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def approx_bell(n, N=20):\n",
    "    s = 0\n",
    "    for k in range(N + 1):\n",
    "        s = s + k ** n / fact(k)\n",
    "    return s / math.e"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quelle valeur de $N$ prendre pour avoir une bonne approximation ? That is the question. Tout à fait empiriquement, pour $n=30$ et $N=34$, tous les chiffres de l'approximation flottante de $B_n$ sont exacts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 34\n",
    "n = 30\n",
    "print(approx_bell(n, N - 1))\n",
    "print(approx_bell(n, N))\n",
    "print(approx_bell(n, N + 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (8, 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ns = range(31)\n",
    "bs = [bell(n) for n in ns]\n",
    "plt.semilogy(ns, bs, 'k', lw=1)\n",
    "bs1 = [approx_bell(n) for n in ns]\n",
    "plt.semilogy(ns, bs1, 'or', ms=4)\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}