{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Les coefficients binomiaux\n",
    "\n",
    "Marc Lorenzi\n",
    "\n",
    "27 octobre 2020"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Pari** : Calculer $\\binom {2 \\text{ millions}}{1\\text{ million}}$ en moins d'une minute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "import random\n",
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (10, 6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1 Qu'allons nous faire ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nous allons dans ce notebook nous intéresser aux **coefficients binomiaux**. Pour tous entiers $n,k\\in\\mathbb N$, nous allons définir un nombre $\\binom n k\\in\\mathbb N$. Nous allons étudier quelques propriétés bien connues de ces nombres, et décrire quelques algorithmes permettant de les calculer.\n",
    "\n",
    "Dans la dernière section nous obtiendrons un algorithme très efficace permettant de factoriser les coefficients binomiaux. De cet algorithme de factorisation nous déduirons facilement un algorithme de calcul des coefficients binomiaux. Celui-ci nous permattra de calculer en un temps très raisonnable le coefficient $\\binom {2000000}{1000000}$, qui est un nombre d'environ 600000 chiffres."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 Compteurs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Définissons une classe `Compteur`. Un objet de cette classe est un *compteur*. On peut l'incrémenter ou le décrémenter (de 1 par défaut) par les méthodes `incr` et `decr`, ou lui donner une valeur précise (0 par défaut) par la méthode `reset`. On peut également obtenir la valeur du compteur par la méthode `val`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Compteur:\n",
    "    \n",
    "    def __init__(self, val=0): self.val = val\n",
    "        \n",
    "    def incr(self, step=1): self.val = self.val + step\n",
    "    def decr(self, step=1): self.val = self.val - step\n",
    "    def reset(self, val=0): self.val = val\n",
    "    def val(self): return self.val\n",
    "        \n",
    "    def __str__(self): return str(self.val)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Définissons un compteur global `CPTR`. Celui-ci nous servira dans quelques unes de nos expériences futures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "CPTR = Compteur()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Faisons un petit test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "CPTR.reset()\n",
    "for k in range(10000):\n",
    "    r = random.randint(0,1)\n",
    "    if r == 0: CPTR.incr()\n",
    "    else: CPTR.decr()\n",
    "print(CPTR)\n",
    "CPTR.reset()\n",
    "print(CPTR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Coefficients binomiaux"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nous allons dans ce qui suit parler de cardinaux. Voici ce qu'il y a à savoir :\n",
    "\n",
    "- Un ensemble $A$ est dit **fini** lorsqu'il existe $n\\in\\mathbb N$ tel que $A$ soit en bijection avec $[|1,n|]$. Un tel entier $n$ est alors unique. On l'appelle le cardinal de $A$ et on le note $|A|$.\n",
    "- Deux ensembles finis $A$ et $B$ ont le même cardinal si et seulement si il existe une bijection de $A$ sur $B$.\n",
    "- Une réunion d'ensembles finis **disjoints** est encore un ensemble fini, et son cardinal est la somme de leurs cardinaux.\n",
    "- Des « évidences », comme par exemple le fait que le cardinal d'un sous-ensemble est plus petit que le cardinal de l'ensemble.\n",
    "\n",
    "Tout ceci sera détaillé dans le cours de dénombrement."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 Parties d'un ensemble"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.1.1 Parties quelconques"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Notation.** Soit $E$ un ensemble fini. On note $\\mathcal P(E)$ l'ensemble des parties de $E$. \n",
    "\n",
    "**Proposition.** Soit $E$ un ensemble fini de cardinal $n$. On a \n",
    "\n",
    "> $$|\\mathcal P(E)|=2^n$$\n",
    "\n",
    "**Démonstration.** on procède par récurrence sur $n$.\n",
    "\n",
    "- Le seul ensemble de cardinal 0 est $\\emptyset$, qui possède $2^0=1$ partie : lui-même.\n",
    "- Soit $n\\in\\mathbb N$. Supposons que pour tout ensemble fini $E$ de cardinal $n$, on a $|\\mathcal P(E)|=2^n$. Soit $E$ un ensemble fini de cardinal $n+1$. Donnons-nous $a\\in E$ et écrivons $E=E'\\cup \\{a\\}$ où $E'$ est de cardinal $n$. On a alors\n",
    "\n",
    "$$\\mathcal P(E)=\\mathcal P'(E)\\cup \\mathcal P''(E)$$\n",
    "\n",
    "où\n",
    "\n",
    "$$\\mathcal P'(E)=\\{X\\in\\mathcal P(E), a\\in X\\}$$\n",
    "\n",
    "$$\\mathcal P''(E)=\\{X\\in\\mathcal P(E), a\\not\\in X\\}=\\mathcal P(E')$$\n",
    "\n",
    "Clairement, l'application $X\\mapsto X\\cup\\{a\\}$ est une bijection de $\\mathcal P(E')$ sur $\\mathcal P'(E)$. Ainsi, $|\\mathcal P(E')|=|\\mathcal P'(E)|$. De là,\n",
    "\n",
    "$$|\\mathcal P(E)|=|\\mathcal P'(E)|+|\\mathcal P''(E)|=2|\\mathcal P(E')|=2^{n+1}$$\n",
    "\n",
    "en appliquant l'hypothèse de récurrence à $E'$, qui est de cardinal $n$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.1.2 Parties de cardinal donné"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Notation.** Soit $E$ un ensemble fini. Pour tout $k\\in\\mathbb N$, on note $\\mathcal P_k(E)$ l'ensemble des parties de $E$ de cardinal $k$. On note également\n",
    "\n",
    "$$\\binom n k = |\\mathcal P_k(E)|$$\n",
    "\n",
    "Les entiers $\\binom n k$ sont appelés les **coefficients binomiaux** (lire « $k$ parmi $n$ », c'est le « nombre de façons » de « choisir » $k$ objets parmi $n$ objets). Ils sont définis pour $n,k\\in\\mathbb N$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.1.3 Deux formules bien connues"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Proposition.** Soit $n\\in\\mathbb N$. On a\n",
    "\n",
    "> $$\\sum_{k=0}^n\\binom n k=2^n$$\n",
    "\n",
    "**Démonstration.** Soit $E$ un ensemble de cardinal $n$. On a\n",
    "\n",
    "$$\\mathcal P(E)=\\bigcup_{k=0}^n \\mathcal P_k(E)$$\n",
    "\n",
    "et cette réunion est disjointe. Ainsi,\n",
    "\n",
    "$$2^n=|\\mathcal P(E)|=\\sum_{k=0}^n |\\mathcal P_k(E)|=\\sum_{k=0}^n \\binom n k$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Proposition.** Soit $n\\in\\mathbb N$. Soit $0\\le k\\le n$. On a\n",
    "\n",
    "> $$\\binom n k = \\binom n {n-k}$$\n",
    "\n",
    "**Démonstration.** Soit $E$ un ensemble de cardinal $n$. L'application $A\\mapsto E\\setminus A$ est une bijection de $\\mathcal P_k(E)$ sur $\\mathcal P_{n-k}(E)$. Ces deux ensembles ont donc le même cardinal."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Une première fonction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notre première fonction de calcul des coefficients binomiaux recherche toutes les parties à $k$ éléments d'un ensemble à $n$ éléments. Puis elle compte combien elle a trouvé de parties.\n",
    "\n",
    "Comment créer toutes les parties à $k$ éléments d'un ensemble $E$ de cardinal $n$ ? \n",
    "\n",
    "- Si $k < 0$, il n'y a aucune telle partie.\n",
    "- Si $k = 0$, il y a une seule partie, $\\emptyset$.\n",
    "- Si $k\\ge 1$ et $n=0$ il n'y a aucune partie.\n",
    "- Si $k\\ge 1$ et $n\\ge 1$, prenons $a\\in E$ et écrivons \n",
    "\n",
    "$$E=E'\\cup\\{a\\}$$\n",
    "\n",
    "où $E'$ est de cardinal $n-1$. Posons\n",
    "\n",
    "$$\\mathcal P'_k(E)=\\{X\\in\\mathcal P_k(E), a\\in X\\}$$\n",
    "\n",
    "et\n",
    "\n",
    "$$\\mathcal P''_k(E)=\\{X\\in\\mathcal P_k(E), a\\not\\in X\\}$$\n",
    "\n",
    "Ces deux ensembles sont disjoints, et\n",
    "\n",
    "$$\\mathcal P_k(E)=\\mathcal P'_k(E)\\cup \\mathcal P''_k(E)$$\n",
    "\n",
    "La fonction `parties` prend en paramètres un ensemble $E$ représenté par une liste Python et un entier $k\\in\\mathbb Z$. Elle renvoie la liste des parties de $E$ de cardinal $k$. Dans le cas intéressant ($n,k\\ge 1$), cette fonction utilise ce que nous venons de dire, en prenant $a=E[0]$ et $E'=E\\setminus\\{a\\}=E[1:]$.\n",
    "\n",
    "**Remarque.** Nous comptons, grâce au compteur `CPTR` le nombre de concaténations de listes effectuées par la fonction. Ainsi, à l'avant dernière ligne, `CPTR` est incrémenté de 2 car à la ligne suivante on effectue deux concaténations de listes (deux signes $+$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def parties(E, k):\n",
    "    if k == 0: return [[]]\n",
    "    else:\n",
    "        n = len(E)\n",
    "        if n == 0: return []\n",
    "        else:\n",
    "            a = E[0]\n",
    "            E1 = E[1:]\n",
    "            P1 = parties(E1, k - 1)\n",
    "            P2 = parties(E1, k)\n",
    "            CPTR.incr(2)\n",
    "            return P2 + [[a] + X for X in P1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(parties([1, 2, 3, 4, 5], 3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voici notre première fonction de calcul des coefficients binomiaux."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def binomial0(n, k):\n",
    "    return len(parties(list(range(n)), k))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Testons, en affichant le **triangle de Pascal** (la fonction `fmt` nous servira pour afficher de façon jolie des listes d'entiers. Nous verrons plus loin comment nous en servir)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fmt(n, d=6):\n",
    "    return n * ('%' + str(d) + 'd')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 11\n",
    "for n in range(N):\n",
    "    print(fmt(n + 1, 4) % tuple([binomial0(n, k) for k in range(n + 1)]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bien entendu, cette fonction est terriblement inefficace, en temps comme en espace."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "CPTR.reset()\n",
    "print(binomial0(20, 10))\n",
    "print('CPTR =', CPTR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Le calcul de $\\binom {20} {10}$ nécessite environ 1 million de concaténations de listes ! De plus, la fonction crée 184756 listes d'entiers de longueur 10, ce qui nécessite le stockage de 1847560 entiers.\n",
    "\n",
    "Il va nous falloir faire mieux que cela. Première idée, nous n'avons pas besoin de **chercher** les parties d'un ensemble, mais seulement de les **compter**. Mais avant, faisons une petite digression."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3 La formule du binôme"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Si les coefficients binomiaux s'appellent ainsi, c'est parce qu'ils apparaissent dans la **formule du binôme**. On se propose de montrer cette formule uniquement avec ce qui précède dans le notebook ... c'est à dire la définition des coefficients binomiaux.\n",
    "\n",
    "Donnons-nous $2n$ réels (ou complexes, ou autres $(*)$) $a_1,\\ldots,a_n,b_1,\\ldots,b_n$ et développons le produit\n",
    "\n",
    "$$P = (a_1+b_1)(a_2+b_2)\\ldots(a_n+b_n)=F_1F_2\\ldots F_n$$\n",
    "\n",
    "$P$ est un produit de sommes. Développer $P$, c'est écrire $P$ comme une somme de produits. Chacun des termes de la somme développée est un produit de $n$ facteurs. Ces facteurs sont obtenus en choisissant, pour chacun des facteurs $F_i$ de $P$, l'un des termes $a_i$ ou $b_i$. On fait cela pour tous les choix possibles. En gros,\n",
    "\n",
    "$$P=\\sum_{\\text{tous les choix}}x_1\\dots x_n$$\n",
    "\n",
    "où $x_i$ vaut $a_i$ ou $b_i$. Comment préciser ces choix ? Il suffit de dire quels sont les facteurs de $P$ dans lesquels on choisit $a_i$ (pour les autres, on choisit $b_i$). Un tel choix de facteurs est caractérisé par la donnée d'une partie $A$ de l'ensemble $E=[|1,n|]$ : $A$ est l'ensemble des numéros $i$ des facteurs $F_i$ pour lesquels on choisit $a_i$. On somme ensuite pour toutes les parties $A$ possibles. On a ainsi\n",
    "\n",
    "$$P=\\sum_{A\\in\\mathcal P(E)}\\prod_{i\\in A}a_i\\prod_{i\\not\\in A}b_i$$\n",
    "\n",
    "Prenons maintenant tous les $a_i$ égaux à un réel $a$, et tous les $b_i$ égaux à un réel $b$. On a\n",
    "\n",
    "$$P=\\sum_{A\\in\\mathcal P(E)}\\prod_{i\\in A}a\\prod_{i\\not\\in A}b=\\sum_{A\\in\\mathcal P(E)}a^{|A|}b^{|E\\setminus A|}=\\sum_{A\\in\\mathcal P(E)}a^{|A|}b^{n-|A|}$$\n",
    "\n",
    "Réordonnons cette somme en regroupant les ensembles $A$ de même cardinal $k$, puis en sommant sur toutes les valeurs de $k$ (c'est à dire $0,1,\\ldots,n$). Il vient\n",
    "\n",
    "$$\\begin{array}{lll}\n",
    "P&=&\\sum_{k=0}^n\\sum_{A\\in\\mathcal P_k(E)}a^{|A|}b^{n-|A|}\\\\\n",
    "&=&\\sum_{k=0}^n\\sum_{A\\in\\mathcal P_k(E)}a^{k}b^{n-k}\\\\\n",
    "&=&\\sum_{k=0}^n\\left(\\sum_{A\\in\\mathcal P_k(E)} 1\\right)a^{k}b^{n-k}\\\\\n",
    "\\end{array}$$\n",
    "\n",
    "Mais $\\sum_{A\\in\\mathcal P_k(E)} 1=|\\mathcal P_k(E)|=\\binom n k$. Ainsi,\n",
    "\n",
    "$$P=\\sum_{k=0}^n\\binom n ka^{k}b^{n-k}$$\n",
    "\n",
    "Nous venons de montrer la formule du binôme.\n",
    "\n",
    "**Théorème.** Soient $a,b\\in\\mathbb R$. Soit $n\\in\\mathbb N$. On a\n",
    "\n",
    "$$(a+b)^n=\\sum_{k=0}^n\\binom n ka^{k}b^{n-k}$$\n",
    "\n",
    "$(*)$ *De façon très générale, la formule du binôme est vraie lorsque $a,b\\in \\mathbb A$ où $\\mathbb A$ est un **anneau**, et que de plus $ab=ba$.* "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.4 Une relation de récurrence"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Proposition.** Pour tous $n, k\\in\\mathbb N^*$,\n",
    "\n",
    "> $$\\binom n k = \\binom {n-1}{k-1}+\\binom {n-1}k$$\n",
    "\n",
    "**Démonstration.** Reprenons les notations vues dans le paragraphe 2.2. Prenons $a\\in E$ et écrivons \n",
    "\n",
    "$$E=E'\\cup\\{a\\}$$\n",
    "\n",
    "où $E'$ est de cardinal $n-1$. Posons\n",
    "\n",
    "$$\\mathcal P'_k(E)=\\{X\\in\\mathcal P_k(E), a\\in X\\}$$\n",
    "\n",
    "et\n",
    "\n",
    "$$\\mathcal P''_k(E)=\\{X\\in\\mathcal P_k(E), a\\not\\in X\\}$$\n",
    "\n",
    "Ces deux ensembles sont disjoints, et\n",
    "\n",
    "$$\\mathcal P_k(E)=\\mathcal P'_k(E)\\cup \\mathcal P''_k(E)$$\n",
    "\n",
    "Faisons les remarques suivantes :\n",
    "\n",
    "- $\\mathcal P''_k(E)=\\mathcal P_k(E')$, donc $|\\mathcal P''_k(E)|=|\\mathcal P_k(E')|=\\binom {n-1} k$.\n",
    "- $\\mathcal P'_k(E)=\\{X\\cup\\{a\\}, X\\in\\mathcal P_{k-1}(E')\\}$. L'application $X\\mapsto X\\cup\\{a\\}$ est une bijection de $\\mathcal P_{k-1}(E')$ sur $\\mathcal P'_k(E)$, donc $|\\mathcal P'_k(E)|=|\\mathcal P_{k-1}(E')|=\\binom {n-1} {k-1}$.\n",
    "\n",
    "De là,\n",
    "\n",
    "$$\\binom n k = |\\mathcal P_k(E)|=|\\mathcal P'_k(E)|+|\\mathcal P''_k(E)|=\\binom{n-1}{k-1}+\\binom{n-1}{k}$$\n",
    "\n",
    "Cette relation permet d'écrire une fonction qui calcule les coefficients binomiaux. Appelons-la `binomial1`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def binomial1(n, k):\n",
    "    if k == 0: return 1\n",
    "    elif n == 0: return 0\n",
    "    else:\n",
    "        CPTR.incr()\n",
    "        return binomial1(n - 1, k - 1) + binomial1(n - 1, k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 11\n",
    "for n in range(N):\n",
    "    print(fmt(n + 1, 4) % tuple([binomial1(n, k) for k in range(n + 1)]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sans surprise, on retrouve la même table que celle que nous avions obtenue avec `binomial0`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.5 Complexité de la fonction `binomial1`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pour tout entiers naturels $n$ et $k$, notons $A_{n,k}$ le nombre d'additions effectuées par `binomial1` lors du calcul de $\\binom n k$. On a\n",
    "\n",
    "- Pour tout $n\\ge 0$, $A_{n,0}=0$.\n",
    "- Pour tout $k > 0$, $A_{0,k}=0$.\n",
    "- Pour tout $n\\ge 1$ et tout $k\\ge 1$, $A_{n,k}=A_{n-1,k-1}+A_{n-1,k}+1$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def A(n, k):\n",
    "    if k == 0: return 0\n",
    "    elif n == 0: return 0\n",
    "    else: return A(n - 1, k - 1) + A(n - 1, k) + 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voici les premières valeurs de $A_{n,k}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 11\n",
    "for n in range(N):\n",
    "    print(fmt(n + 1, 5) % tuple([A(n, k) for k in range(n + 1)]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Histoire de nous rassurer, contrôlons par exemple le nombre d'additions effectuées pour le calcul de $\\binom {10} 7$. Le tableau nous dit 967. Et dans la réalité ? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "CPTR.reset()\n",
    "x = binomial1(10, 7)\n",
    "print(CPTR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nous voici rassurés."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Peut-on obtenir une formule explicite pour $A_{n,k}$ ? Oui, plus ou moins. Posons $B_{n,k}=A_{n,k}+1$. On a alors\n",
    "\n",
    "- Pour tout $n\\ge 0$, $B_{n,0}=1$.\n",
    "- Pour tout $k > 0$, $B_{0,k}=1$.\n",
    "- Pour tout $n\\ge 1$ et tout $k\\ge 1$, $B_{n,k}=B_{n-1,k-1}+B_{n-1,k}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 11\n",
    "for n in range(N):\n",
    "    print(fmt(n + 1, 5) % tuple([A(n, k) + 1 for k in range(n + 1)]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Proposition.** On a pour tous entiers $n,k\\ge 0$\n",
    "\n",
    "$$B_{n,k}=\\sum_{j=0}^{k} \\binom n j$$\n",
    "\n",
    "**Démonstration.** On fait une récurrence sur $n$. C'est clair pour $n=0$ puisque dans ce cas $B_{n,k}=1$. Soit $n\\in\\mathbb N^*$. Supposons que pour tout $k\\in\\mathbb N$, $B_{n-1,k}=\\sum_{j=0}^{k} \\binom {n-1} j$. Montrons que ceci est encore vrai pour $n$. \n",
    "\n",
    "Pour $k=0$, la propriété est vraie puisque $B_{n,0}=1$. Soit $k\\in\\mathbb N^*$. On a\n",
    "\n",
    "$$\\begin{array}{lll}\n",
    "B_{n,k}&=&B_{n-1,k-1}+B_{n-1,k}\\\\\n",
    "&=&\\sum_{j=0}^{k-1} \\binom {n-1} j+\\sum_{j=0}^{k} \\binom {n-1} j\\\\\n",
    "&=&\\sum_{j=1}^{k} \\binom {n-1} {j-1}+\\sum_{j=0}^{k} \\binom {n-1} j\\\\\n",
    "&=&\\binom {n-1} 0+\\sum_{j=1}^{k} \\left(\\binom {n-1} {j-1}+\\binom {n-1} j\\right)\\\\\n",
    "&=&1+\\sum_{j=1}^{k} \\binom {n} j\\\\\n",
    "&=&\\sum_{j=0}^{k} \\binom {n} j\\\\\n",
    "\\end{array}$$\n",
    "\n",
    "**Corollaire.** On a pour tous entiers $n,k\\ge 0$\n",
    "\n",
    "$$A_{n,k}=\\sum_{j=1}^{k} \\binom n j$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La valeur obtenue pour $A_{n,k}$ n'est pas du tout rassurante. Par exemple, dans le cas où $k=n$, on obtient $A_{n,n}=2^n-1$. Il faut $2^n-1$ additions pour calculer $\\binom n n$ avec notre fonction `binomial1`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "CPTR.reset()\n",
    "x = binomial1(20, 20)\n",
    "print(CPTR)\n",
    "print(2 ** 20 - 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Inutile de dire que ce n'est pas avec une telle fonction que nous calculerons des coefficients binomiaux où $n=2000000$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. La taille des coefficients binomiaux"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Avant de nous plonger dans de nouveaux algorithmes, posons nous quelques questions sur les coefficients binomiaux. Pour un entier $n$ fixé, quel est le comportement de la suite $\\left(\\binom n k\\right)_{0\\le k\\le n}$ ? Quelle est la valeur maximale de cette suite ? Peut-on obtenir un équivalent simple de ce maximum ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 Représentation graphique"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `plot_binom` prend en paramètre un entier $n$. Elle trace le graphe des coefficients binomiaux $\\binom n k$ pour $0\\le k\\le n$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_binom(n):\n",
    "    xs = list(range(n + 1))\n",
    "    ys = [binomial1(n, k) for k in range(n + 1)]\n",
    "    zs = [1]\n",
    "    for k in range(n + 1):\n",
    "        plt.fill([k - 0.5, k + 0.5, k + 0.5, k - 0.5, k - 0.5], [0, 0, ys[k], ys[k], 0], 'r')\n",
    "        plt.plot([k - 0.5, k + 0.5, k + 0.5, k - 0.5, k - 0.5], [0, 0, ys[k], ys[k], 0], 'k')\n",
    "    plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Traçons pour $n=14$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_binom(14)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Et aussi pour $n=15$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_binom(15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ça monte puis ça descend. Et le maximum est au milieu, avec une légère différence selon que $n$ est pair (max en un unique point) ou impair (max en deux points)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Confirmations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soient $n,k\\in\\mathbb N$, $n\\ge 1$, $0\\le k \\le n-1$. On a \n",
    "\n",
    "$$\\begin{array}{ccc}\n",
    "\\binom n {k+1} - \\binom n k &=& \\frac{n!}{(k+1)!(n-1-k)!}-\\frac{n!}{k!(n-k)!}\\\\\n",
    "&=& \\frac{n!}{(k+1)!(n-k)!}\\left(n-2k-1\\right)\\\\\n",
    "\\end{array}$$\n",
    "\n",
    "Cette quantité est du signe de $n-2k-1$. Discutons selon la parité de $n$.\n",
    "\n",
    "- Cas 1, $n=2p$ où $p\\ge 1$. On a $n-2k-1> 0$ si et seulement si $k< p-\\frac 1 2$ c'est à dire, puisque $k$ est entier, $k\\le p-1=\\frac n 2 - 1$. Ainsi, la suite $\\left(\\binom n k\\right)_{0\\le k\\le n}$ croît strictement pour $0\\le k\\le p$, passe par un maximum pour $k=p$, puis décroît strictement.\n",
    "\n",
    "- Cas 2, $n=2p+1$ où $p\\ge 0$. On a $n-2k-1> 0$ si et seulement si $k< p$. Ainsi, la suite $\\left(\\binom n k\\right)_{0\\le k\\le n}$ croît strictement pour $0\\le k\\le p$, prend deux valeurs égales pour $k=p$ et $k=p+1$, puis décroît strictement.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3 La valeur du maximum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Que « vaut » le coefficient maximal ? Regardons le cas où $n=2p$ est pair (le cas $n$ impair est similaire). Le coefficient maximal est obtenu pour $k=p$, et vaut $\\binom {2p}p$. Calculons un équivalent de ce coefficient lorsque $p$ tend vers l'infini, en utilisant la formule de Stirling. Cette formule nous dit que\n",
    "\n",
    "$$n!\\sim \\left(\\frac n e\\right)^n\\sqrt{2\\pi n}$$\n",
    "\n",
    "On a donc\n",
    "\n",
    "$$(2p)!\\sim \\left(\\frac{2p}e\\right)^{2p}\\sqrt{2\\pi 2p}$$\n",
    "\n",
    "et\n",
    "\n",
    "$$p!^2\\sim \\left(\\frac{p}e\\right)^{2p}{2\\pi p}$$\n",
    "\n",
    "Ainsi, après simplifications,\n",
    "\n",
    "$$\\binom{2p}p=\\frac{(2p)!}{p!^2}\\sim\\frac{2^{2p}}{\\sqrt{\\pi p}}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def equiv_max(n):\n",
    "    p = n / 2\n",
    "    return 2 ** n / math.sqrt(math.pi * n / 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(equiv_max(100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Les coefficients binomiaux peuvent donc être énormes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Une méthode un peu moins naïve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1 Une formule pour les coefficients binomiaux"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Proposition.** Pour tous entiers $n,k\\in\\mathbb N$ tels que $0\\le k\\le n$, on a\n",
    "\n",
    "> $$\\binom n k = \\frac{n!}{k!(n-k)!}$$\n",
    "\n",
    "**Démonstration.** On procède par récurrence sur $n$.\n",
    "\n",
    "- $n=0$ impose $k=0$ et la formule est clairement vraie dans ce cas, puisque $\\binom 0 0=1$ et $0!=1$.\n",
    "- Soit $n\\ge 1$. Supposons la propriété vraie pour l'entier $n-1$. Soit $0\\le k\\le n$. Supposons tout d'abord $k\\ne 0$ et $k\\ne n$. On a alors\n",
    "\n",
    "$$\\begin{array}{lll}\n",
    "\\binom n k &=&\\binom{n-1}{k-1}+\\binom{n-1}k\\\\\n",
    "&=&\\frac{(n-1)!}{(k-1)!(n-k)!}+\\frac{(n-1)!}{k!(n-1-k)!}\\\\\n",
    "&=&\\frac{(n-1)!}{k!(n-k)!}\\left(k + (n-k)\\right)\\\\\n",
    "&=&\\frac{n!}{k!(n-k)!}\n",
    "\\end{array}$$\n",
    "\n",
    "Si $k=0$, ou $k=n$, la formule est clairement vraie. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.2 Fonction Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Posons, pour $n\\ge 0$ et $k\\ge 1$,\n",
    "\n",
    "$$n^{\\underline k}=n(n-1)\\ldots (n-k+1)$$\n",
    "\n",
    "Posons également $n^{\\underline 0}=1$. Nous appellerons $n^{\\underline k}$ la $k$ième **puissance descendante** de $n$. On a pour tous entiers $n,k\\ge 0$ tels que $0\\le k\\le n$,\n",
    "\n",
    "$$\\binom n k=\\frac{n!}{k!(n-k)!}=\\frac{n^{\\underline k}}{k!}=\\frac{n^{\\underline k}}{k^{\\underline k}}$$\n",
    "\n",
    "Le calcul des coefficients binomiaux se ramène donc à celui des puissances descendantes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def power_down(n, k):\n",
    "    p = 1\n",
    "    for j in range(k): p = p * (n - j)\n",
    "    return p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def binomial2(n, k):\n",
    "    if k < 0 or k > n: return 0\n",
    "    else:\n",
    "        return power_down(n, k) // power_down(k, k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 11\n",
    "for n in range(N):\n",
    "    print(fmt(n + 1, 4) % tuple([binomial2(n, k) for k in range(n + 1)]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Le calcul de `power_down(n, k)` demande $k$ multiplications. Celui de $\\binom n k$ par cette méthode demande donc $2k$ multiplications et une division. Évidemment, ces multiplications peuvent être très coûteuses, car les factorielles deviennent énormes lorsque $n$ augmente.\n",
    "\n",
    "Notre fonction `binomial2` est nettement plus efficace que `binomial1`. Nous pouvons maintenant calculer des coefficients binomiaux de taille plus importante."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rappelez-vous notre estimation de $\\binom {100}{50}$, qui était $1.01\\times 10^{29}$. Quelle est la valeur exacte de ce coefficient ? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = binomial2(100, 50)\n",
    "print(x)\n",
    "print(float(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "binomial2(2000, 1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Regardons le temps mis par `binomial2` pour calculer de « grands » coefficients binomiaux. Pour la seconde cellule ci-dessous un peu de patience est nécessaire."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t1 = time.time()\n",
    "x = binomial2(100000, 50000)\n",
    "t2 = time.time()\n",
    "print(t2 - t1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t1 = time.time()\n",
    "x = binomial2(200000, 100000)\n",
    "t2 = time.time()\n",
    "print(t2 - t1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "J'ai tenté `binomial2(2000000, 1000000)` et arrêté le calcul après une dizaine de minutes. Voyez un peu plus loin pour une estimation du temps nécessaire à ce calcul avec la fonction `binomial2`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3 Temps réel de calcul"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Le calcul de $\\binom{2n}{n}$ nécessite donc $4n$ multiplications et une division. Pour $n=10^6$, il faut donc effectuer 4 millions de multiplications (et une division). Mais nous faisons des multiplications de **grands** nombres !\n",
    "\n",
    "Un algorithme naïf pour multiplier un nombre de $m$ chiffres par un nombre de $n$ chiffres demande (sans entrer dans les détails) de l'ordre de $m\\times n$ opérations. Voir à ce sujet votre cours de CM2.\n",
    "\n",
    "Quel est le nombre de chiffres en base 2 de l'entier $n$ ? Soit $K$ ce nombre de chiffres. On a\n",
    "\n",
    "$$2^{K-1}\\le n < 2^K$$\n",
    "\n",
    "Passant au logarithme en base 2, et ajoutant 1, il vient\n",
    "\n",
    "$$K\\le \\log_2 n + 1< K + 1$$\n",
    "\n",
    "Ainsi,\n",
    "\n",
    "$$K = \\lfloor \\log_2(n)\\rfloor + 1$$\n",
    "\n",
    "La fonction `nb_chiffres` renvoie le nombre de chiffres de l'écriture de $n$ en base 2. Elle renvoie 0 si $n$ est nul."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def nb_chiffres(n):\n",
    "    if n == 0: return 0\n",
    "    else:\n",
    "        return math.floor(math.log(n, 2)) + 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Réécrivons `power_down` en incrémentant `CPTR` de la valeur adéquate à chaque multiplication."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def power_down2(n, k):\n",
    "    p = 1\n",
    "    for j in range(k): \n",
    "        CPTR.incr(nb_chiffres(p) * nb_chiffres(n - j))\n",
    "        p = p * (n - j)\n",
    "    return p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def binomial3(n, k):\n",
    "    if k < 0 or k > n: return 0\n",
    "    else:\n",
    "        CPTR.incr() # pour la division\n",
    "        return power_down2(n, k) // power_down2(k, k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_binomial3(N):\n",
    "    s = []\n",
    "    ns = list(range(N, 11 * N, N))\n",
    "    for n in ns:\n",
    "        CPTR.reset()\n",
    "        x = binomial3(2 * n, n)\n",
    "        s.append(CPTR.val)\n",
    "    plt.loglog(ns, s, '-ok')\n",
    "    plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Le tracé ci-dessous est en coordonnées logarithmiques."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_binomial3(1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plus réalistement, nous pouvons tracer la courbe des « vrais » temps d'exécution. Voici ce que l'on obtient (cela dépend évidemment de la machine utilisée)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_binomial3bis(N):\n",
    "    s = []\n",
    "    ns = list(range(N, 11 * N, N))\n",
    "    for n in ns:\n",
    "        t1 = time.time()\n",
    "        x = binomial3(2 * n, n)\n",
    "        t2 = time.time()\n",
    "        s.append(t2 - t1)\n",
    "    plt.loglog(ns, s, '-ok')\n",
    "    plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_binomial3bis(1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On constate que la courbe est approximativement une droite. Si nous appelons $T_n$ le temps mis pour calculer $\\binom {2n}n$ on a donc à peu près\n",
    "\n",
    "$$\\ln T_n = c + \\alpha \\log n$$\n",
    "\n",
    "où $c$ et $\\alpha$ sont deux constantes. $\\alpha$ est la pente de la droite."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t1 = time.time()\n",
    "na = 50000\n",
    "x = binomial2(2 * na, na)\n",
    "t2 = time.time()\n",
    "ta = t2 - t1\n",
    "\n",
    "t1 = time.time()\n",
    "nb = 100000\n",
    "x = binomial2(2 * nb, nb)\n",
    "t2 = time.time()\n",
    "\n",
    "tb = t2 - t1\n",
    "\n",
    "alpha = (math.log(tb) - math.log(ta)) / (math.log(nb) - math.log(na))\n",
    "print(alpha)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Obtenir la constante $c$ est maintenant facile."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c = math.log(tb) - alpha * math.log(nb)\n",
    "print(c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Maintenant, $T_n = e^c n^\\alpha$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def temps_binomial2(n):\n",
    "    k = math.exp(c)\n",
    "    return k * n ** alpha"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vérifions la pertinence de tout cela avec des valeurs de $n$ pour lesquelles nous avons déjà mesuré le temps de calcul."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "temps_binomial2(50000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "temps_binomial2(100000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nous retrouvons à peu près les temps que nous avions mesuré. Notre modèle n'est donc pas trop mauvais. Et si $n=10^6$ ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "temps_binomial2(1000000) / 60"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Il faudrait laisser tourner notre machine un peu plus de 20 minutes. Je vous laisse tenter l'expérience ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Factorisation des coefficients binomiaux"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.1 Valuation $p$-adique d'un entier"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notons $\\mathcal P$ l'ensemble des nombres premiers. Pour tout $n\\in\\mathbb N^*$ et tout nombre premier $p$, posons \n",
    "\n",
    "$$\\nu_p(n)=\\max\\{k\\in\\mathbb N, p^k\\mathbin |n\\}$$\n",
    "\n",
    "$\\nu_p(n)$ est appelé la **valuation $p$-adique** de $n$. C'est la plus grande puissance de $p$ que l'on puisse factoriser dans l'entier $n$. Tout entier naturel $n\\ge 1$ peut s'écrire sous la forme\n",
    "\n",
    "$$n=\\prod_{p\\in\\mathcal P}p^{\\nu_p(n)}$$\n",
    "\n",
    "Le produit ci-dessus est en réalité fini, puisque si $p$ est assez grand, $p$ ne divise pas $n$ et donc $\\nu_p(n)=0$.\n",
    "\n",
    "On montre facilement les propriétés suivantes :\n",
    "\n",
    "- Pour tous entiers $m, n\\ge 1$, $\\nu_p(mn)=\\nu_p(m)+\\nu_p(n)$.\n",
    "- Pour tous entiers $m, n\\ge 1$ tels que $m$ divise $n$, $\\nu_p(\\frac n m)=\\nu_p(n)-\\nu_p(m)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.2 Valuation $p$-adique des coefficients binomiaux"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La clé de ce qui va suivre est la proposition suivante :"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Proposition.** Soient $n\\in\\mathbb N$ et $p$ un nombre premier. On a\n",
    "\n",
    "$$\\nu_p(n!)=\\sum_{i=1}^\\infty\\left\\lfloor \\frac n {p^i}\\right\\rfloor$$\n",
    "\n",
    "**Démonstration.** Posons $E=[|1,n|]$. On a $E=\\bigcup_{k\\ge 0} E_k$ où $E_k=\\{x\\in E, \\nu_p(x)=k\\}$. On a alors \n",
    "\n",
    "$$\\begin{array}{lll}\n",
    "\\nu_p(n!)&=&\\sum_{x\\in E} \\nu_p(x)\\\\\n",
    "&=&\\sum_{k\\ge 0}\\sum_{x\\in E_k}\\nu_p(x)\\\\\n",
    "&=&\\sum_{k\\ge 0}k|E_k|\n",
    "\\end{array}$$\n",
    "\n",
    "Considérons maintenant, pour tout $k\\ge 0$, l'ensemble $F_k=\\{x\\in E, p^k|x\\}$. On a \n",
    "\n",
    "- $F_{k+1}\\subset F_k$ \n",
    "- $E_k=F_k\\setminus F_{k+1}$. \n",
    "\n",
    "Ainsi $|E_k|=|F_k|-|F_{k+1}|$. De là,\n",
    "\n",
    "$$\\begin{array}{lll}\n",
    "\\nu_p(n!)&=&\\sum_{k\\ge 0}k(|F_k|-|F_{k+1}|)\\\\\n",
    "&=&\\sum_{k\\ge 0}k|F_k|-\\sum_{k\\ge0}k|F_{k+1}|\\\\\n",
    "&=&\\sum_{k\\ge 0}k|F_k|-\\sum_{k\\ge 1}(k-1)|F_{k}|\\\\\n",
    "&=&0|F_0|+\\sum_{k\\ge 1}(k-(k-1))|F_k|\\\\\n",
    "&=&\\sum_{k\\ge 1}|F_k|\n",
    "\\end{array}$$\n",
    "\n",
    "Pour terminer on remarque que, pour tout $k\\ge 1$, $F_k=\\{p^k,2p^k,\\ldots,\\alpha p^k\\}$ où $\\alpha p^k\\le n < (\\alpha+1)p^k$. Le cardinal de $F_k$ est donc égal à $\\alpha =\\left\\lfloor \\frac n {p^k}\\right\\rfloor$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `val_fact` ci-dessous renvoie la valuation $p$-adique de $n!$. Elle se contente d'appliquer la formule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def val_fact(n, p):\n",
    "    s = 0\n",
    "    while n != 0:\n",
    "        n = n // p\n",
    "        s = s + n\n",
    "    return s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quelle est la valuation dyadique de $1000!$ ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "val_fact(1000, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$1000!$ est divisible par $2^{994}$, mais pas par $2^{995}$.\n",
    "\n",
    "Un corollaire immédiat du théorème précédent nous donne la valuation des coefficients binomiaux."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Corollaire.** Soient $n, k\\in\\mathbb N$ et $p$ un nombre premier. On a\n",
    "\n",
    "$$\\nu_p\\left(\\binom n k\\right)=\\sum_{i=1}^\\infty \\left(\\left\\lfloor \\frac n {p^i}\\right\\rfloor-\\left\\lfloor \\frac k {p^i}\\right\\rfloor - \\left\\lfloor \\frac {n-k} {p^i}\\right\\rfloor\\right)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nous allons voir que l'on peut déduire de cette formule des relations de récurrence permettant de calculer la valuation $p$-adique de $\\binom n k$ récursivement."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.3 Une relation de récurrence"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Notation.** Pour tout entier $n$ et tout nombre premier $p$, notons $\\widehat n=\\left\\lfloor \\frac n p\\right\\rfloor$. L'entier $\\widehat n$ est le quotient de la division euclidienne de $n$ par $p$.\n",
    "\n",
    "\n",
    "Soit $n\\in\\mathbb N$. Soit $p$ un nombre premier. Posons $n=p\\widehat n+c_n$ où $0\\le c_n<p$. On a pour tout $i\\ge 1$,\n",
    "\n",
    "$$\\left\\lfloor \\frac n {p^i}\\right\\rfloor=\\left\\lfloor \\frac {\\widehat n} {p^{i-1}} + \\frac {c_n} {p^i}\\right\\rfloor=\\left\\lfloor \\frac {\\widehat n} {p^{i-1}}\\right\\rfloor$$\n",
    "\n",
    "De là,\n",
    "\n",
    "$$\\nu_p(n!)=\\sum_{i\\ge 1}\\left\\lfloor \\frac n {p^i}\\right\\rfloor=\\widehat n+\\sum_{i\\ge 2}\\left\\lfloor \\frac {\\widehat n} {p^{i-1}}\\right\\rfloor=\\widehat n+\\nu_p(\\widehat n!)$$\n",
    "\n",
    "Ainsi, on a pour tous entiers $n,k$ tels que $0\\le k\\le n$,\n",
    "\n",
    "$$\\nu_p \\binom n k =\\widehat n- \\widehat k-\\widehat{n-k}+\\nu_p(\\widehat n!)-\\nu_p(\\widehat k!)-\\nu_p(\\widehat{n-k}!)$$\n",
    "\n",
    "Remarquons que $n-k=p(\\widehat n-\\widehat k)+(c_n-c_k)$. Plusieurs cas se présentent alors.\n",
    "\n",
    "- Cas 0 : $k=0$. Dans ce cas, $\\nu_p \\binom n k = 0$.\n",
    "- Cas 1 : $c_n\\ge c_k$. Alors, $\\widehat {n-k}=\\widehat n-\\widehat k$ et donc $\\nu_p \\binom n k = \\nu_p \\binom {\\widehat n}{\\widehat k}$.\n",
    "- Cas 2 : $c_n < c_k$. On a alors $n-k=p(\\widehat n-\\widehat k-1)+(c_n-c_k+p)$. Comme $0\\le c_n-c_k+p < p$, on a donc $\\widehat {n-k}=\\widehat n-\\widehat k-1$. Posons $\\widehat n=p\\widehat {\\widehat n}+b_n$ où $0\\le b_n<p$, et de même pour $\\widehat k$.\n",
    "\n",
    "    - Cas 2.1 : $b_n\\ne b_k$. $\\widehat n-\\widehat k$ n'est pas un multiple de $p$, et donc $\\nu_p(\\widehat n-\\widehat k-1)!)=\\nu_p((\\widehat n-\\widehat k)!)$. Dans ce cas,\n",
    "    \n",
    "    $$\\nu_p\\binom n k = \\nu_p \\binom{\\widehat n}{\\widehat k} + 1$$\n",
    "    \n",
    "    - Cas 2.2 : $b_n=b_k\\ne 0$. $\\widehat n$ n'est pas un multiple de $p$, donc $\\nu_p(\\widehat n!)=\\nu_p((\\widehat n-1)!)$. Dans ce cas,\n",
    "    \n",
    "    $$\\nu_p\\binom n k = \\nu_p \\binom{\\widehat n-1}{\\widehat k} + 1$$\n",
    "    \n",
    "    - Cas 2.3 : $b_n=b_k= 0$. $\\widehat k$ est un multiple de $p$, donc $\\widehat k+1$ n'en est pas un. Ainsi, $\\nu_p(\\widehat k!)=\\nu_p((\\widehat k+1)!)$. Dans ce cas,\n",
    "    \n",
    "    $$\\nu_p\\binom n k = \\nu_p \\binom{\\widehat n}{\\widehat k+1} + 1$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nous sommes donc ramenés, pour calculer $\\nu_p \\binom n k$, au calcul de la valuation $p$ adique d'un coefficient binomial avec des paramètres strictement plus petits que $n$ et $k$ (ce n'est pas tout à fait vrai, voir le paragraphe 5.5)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.4 La fonction Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pour tout entier $n$ et tout nombre premier $p$, posons \n",
    "\n",
    "$$n=a_n p^2+b_np+c_n$$ \n",
    "\n",
    "où $a_n=\\widehat{\\widehat n}\\in\\mathbb N$ et $b_n, c_n\\in[|0,p-1|]$. On a $c_n = n \\bmod p$ et $b_n= (\\frac 1 p (n - c_n)) \\bmod p$. \n",
    "\n",
    "La fonction `decomp` ci-dessous prend un entier $n$ en paramètre et renvoie le triplet $(a_n,b_n,c_n)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def decomp(n, p):\n",
    "    c = n % p\n",
    "    n1 = n // p\n",
    "    b = n1 % p\n",
    "    a = n1 // p\n",
    "    return (a, b, c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `val_binom` utilise les relations que nous avons montrées dans le paragraphe précédent pour calculer $\\nu_p \\binom n k$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def val_binom(n, k, p):\n",
    "    if k == 0: return 0\n",
    "    else:\n",
    "        an, bn, cn = decomp(n, p)\n",
    "        ak, bk, ck = decomp(k, p)\n",
    "        if cn >= ck: return val_binom(n // p, k // p, p)\n",
    "        elif bn != bk: return val_binom(n // p, k // p, p) + 1\n",
    "        elif bn != 0: return val_binom(n // p - 1, k // p, p) + 1\n",
    "        else: return val_binom(n // p, k // p + 1, p) + 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Remarque.** Tout ce que nous pouvons dire pour l'instant c'est que **SI** `val_binom(n, k, p)` termine, alors il renvoie le bon résultat. À ceux qui trouveraient cette phrase bizarre, je propose la définition universelle suivante : `def f(x): return f(x)`. Eh bien, si l'appel `f(x)` termine, il renvoie bien $f(x)$. Mais il ne termine pas 😀.\n",
    "\n",
    "\n",
    "Avant de prouver la terminaison de `val_binom`, prenons un exemple pour nous rassurer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(fmt(11) % tuple([binomial2(10, k) for k in range(11)]))\n",
    "print(fmt(11) % tuple([val_binom(10, k, 2) for k in range(11)]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.5 La terminaison et la complexité de `val_binom`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notons $C_{n,k}$ le nombre d'appels récursifs effectués par `val_binom(n, k, p)`, en convenant que $C_{n,k}=\\infty$  si l'appel ne termine pas. Quelles sont les équations vérifiées par $C_{n,k}$ ? Il suffit de renprendre le code de `val_binom`. Avec les notations de cette fonction (et en convenant que $\\infty+1=\\infty$),\n",
    "\n",
    "- $C_{n,0}=0$.\n",
    "- Si $k\\ne 0$ et $c_n\\ge c_k$, $C_{n,k}=C_{\\widehat n,\\widehat k}+1$.\n",
    "- Si $k\\ne 0$, $c_n < c_k$ et $b_n\\ne b_k$, $C_{n,k}=C_{\\widehat n,\\widehat k}+1$.\n",
    "- Si $k\\ne 0$, $c_n < c_k$ et $b_n= b_k\\ne 0$, $C_{n,k}=C_{\\widehat n - 1,\\widehat k}+1$.\n",
    "- Si $k\\ne 0$, $c_n < c_k$ et $b_n= b_k= 0$, $C_{n,k}=C_{\\widehat n - 1,\\widehat k}+1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def complex_val_binom(n, k, p):\n",
    "    if k == 0: return 0\n",
    "    elif k > n: return 0\n",
    "    else:\n",
    "        an, bn, cn = decomp(n, p)\n",
    "        ak, bk, ck = decomp(k, p)\n",
    "        if cn >= ck: return complex_val_binom(n // p, k // p, p) + 1\n",
    "        elif bn != bk: return complex_val_binom(n // p, k // p, p) + 1\n",
    "        elif bn != 0: return complex_val_binom(n // p - 1, k // p, p) + 1\n",
    "        else: return complex_val_binom(n // p, k // p + 1, p) + 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 30\n",
    "p = 2\n",
    "for n in range(1, N):\n",
    "    print(fmt(N, 2) % tuple([complex_val_binom(n, k, p) for k in range(N)]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plus « graphiquement » :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (10, 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_complex(N, p):\n",
    "    m = (N + 1) * [None]\n",
    "    for i in range(N + 1): m[i] = (N + 1) * [None]\n",
    "    for n in range(N + 1):\n",
    "        for k in range(N + 1):\n",
    "            m[n][k] = complex_val_binom(n, k, p)\n",
    "    plt.imshow(m, interpolation='none', cmap='hot')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_complex(511, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tout ceci est très joli mais il est temps de se mettre au travail."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Proposition.** Pour tout entier $n\\ge 1$ et tout entier $k\\le n$, $C_{n,k}\\ne\\infty$ et\n",
    "\n",
    "$$C_{n,k}\\le \\log_p n+ 1$$\n",
    "\n",
    "**Démonstration.** On fait une récurrence forte sur $n$.\n",
    "\n",
    "- Pour $n=1$ c'est évident.\n",
    "- Soit $n>1$. Supposons la propriété vérifiée pour $1\\le n'<n$. \n",
    "\n",
    "$$C_{n,k}= C_{\\widehat n-\\varepsilon,k'}+1$$\n",
    "\n",
    "où $\\varepsilon$ vaut 0 ou 1 et, **admettons-le provisoirement**, $k'\\le \\widehat n-\\varepsilon$. Clairement, $\\widehat n-\\varepsilon<n$. Par l'hypothèse de récurrence, $C_{\\widehat n-\\varepsilon,k'}\\ne\\infty$. Deux cas se présentent.\n",
    "\n",
    "- Si $\\widehat n\\ge 1$, on applique l'hypothèse de récurrence :\n",
    "\n",
    "$$C_{n,k}\\le \\log_p \\widehat n + 1 + 1\\le \\log_p\\frac n p + 2 = \\log_p n + 1$$ \n",
    "\n",
    "- Si $\\widehat n=0$, nous avons $0\\le k\\le n<p$. Avec les notations de la fonction `val_binom`, on a $c_n=n$, $c_k=k$. Ainsi, `val_binom` s'appelle récursivement sur `val_binom(0, 0, p)` qui renvoie imméfiatement 0. Il y a donc 1 appel récursif. Ainsi, $C_{n,k}=1\\le \\log_p n + 1$ puisque $\\log_p n<1$.\n",
    "\n",
    "*La terminaison et la complexité de `val_binom` seront donc prouvées dès que nous aurons montré la proposition suivante.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Lemme.** Soient $n,k\\in\\mathbb N$ vérifiant $k\\le n$. Si `val_binom(n, k, p)` se rappelle récursivement, c'est avec deux paramètres $n',k'\\in\\mathbb N$ tels que $k'\\le n'$.\n",
    "\n",
    "**Démonstration.** Un certain nombre de cas sont à considérer :\n",
    "\n",
    "1. $k\\ne 0$ et $c_n\\ge c_k$ : $n'=\\widehat n$, $k'=\\widehat k$.\n",
    "2. $k\\ne 0$, $c_n < c_k$ et $b_n\\ne b_k$ : $n'=\\widehat n$, $k'=\\widehat k$.\n",
    "3. $k\\ne 0$, $c_n < c_k$ et $b_n= b_k\\ne 0$ : $n'=\\widehat n - 1$, $k'=\\widehat k$.\n",
    "4. $k\\ne 0$, $c_n < c_k$ et $b_n= b_k= 0$ : $n'=\\widehat n$, $k'=\\widehat k + 1$.\n",
    "\n",
    "Puisque $k\\le n$, on a $\\frac k p\\le \\frac n p$. Par croissance de la partie entière on obtient $\\widehat k\\le \\widehat n$. Les cas 1 et 2 sont donc réglés.\n",
    "\n",
    "Peut-on avoir $\\widehat k = \\widehat n$ ? Supposons que cela soit le cas, appelons $a$ la valeur commune de ces deux nombres. On a alors\n",
    "\n",
    "$$a\\le \\frac k p<a+1 \\text{ et } a \\le \\frac n p < a + 1$$\n",
    "\n",
    "On en tire facilement que $n-k<p$. Si l'on est dans le cas 3 ou le cas 4, alors\n",
    "\n",
    "$$n=a_np^2+b_np+c_n\\text{ et }k=a_kp^2+b_kp+c_k$$\n",
    "\n",
    "On a $b_n=b_k$, donc\n",
    "\n",
    "$$n-k=(a_n-a_k)p^2+(b_n-b_k)p+c_n-c_k=(a_n-a_k)p^2+c_n-c_k$$\n",
    "\n",
    "Comme $c_n<c_k$ et $n-k\\ge 0$, on a nécessairement $a_n-a_k\\ge 1$. Ainsi,\n",
    "\n",
    "$$n-k=(a_n-a_k)p^2+c_n-c_k\\ge p^2+c_n-c_k\\ge p^2-p=p(p-1)\\ge p$$\n",
    "\n",
    "On n'a donc pas $n-k<p$. Ainsi, $\\widehat k < \\widehat n$, et donc \n",
    "\n",
    "- Pour le cas 3 : $k'=\\widehat k \\le \\widehat n - 1 = n'$.\n",
    "- Pour le cas 4 : $k'=\\widehat k +1 \\le \\widehat n = n'$.\n",
    "\n",
    "Les cas 3 et 4 sont donc aussi réglés."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**La fonction `val_binom` est donc de complexité logarithmique en $n$, c'est à dire linéaire en le nombre de chiffres de $n$ (en base $p$). On pouvait difficilement espérer mieux !**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.6 Factorisation des coefficients binomiaux"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nous voici maintenant capables de factoriser les coefficients binomiaux. Pour factoriser $\\binom n k$, effectuer les opérations suivantes :\n",
    "\n",
    "- Remarquer que si un nombre premier $p$ divise $n!$, alors il divise un entier inférieur ou égal à $n$. Ainsi, $p\\le n$.\n",
    "- Pour tout nombre premier $p\\le n$, déterminer $\\nu_p \\binom n k$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voici le code Python. La fonction `est_premier` teste naïvement si $p$ est un nombre premier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def est_premier(p):\n",
    "    k = 2\n",
    "    while k * k <= p and p % k != 0: k = k + 1\n",
    "    return k * k > p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print([p for p in range(2, 100) if est_premier(p)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `factor_binomial` renvoie la factorisation de $\\binom n k$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def factor_binomial(n, k):\n",
    "    ps = [p for p in range(2, n + 1) if est_premier(p)]\n",
    "    s = []\n",
    "    for p in ps:\n",
    "        m = val_binom(n, k, p)\n",
    "        if m != 0: s.append((p, m))\n",
    "    return s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "L'exemple ci-dessous est pris dans l'article donné en référence à la fin du notebook. Cet article date de 1987, l'auteur avait alors obtenu la factorisation en 1/2 seconde sur un PC équipé d'un processeur 8086 à 8 MHz avec un programme écrit en langage Pascal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(factor_binomial(1000, 353))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notre machine donne le résultat en zéro seconde, mais nous n'avons aucun mérite. C'est juste qu'elle est équipée d'un processeur à 4 coeurs tournant à 2.93 GHz 😀."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.7 Une fonction de calcul rapide des coefficients binomiaux"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nous voici en possession d'un nouvel algorithme permettant de calculer les coefficients binomiaux. En effet,\n",
    "\n",
    "$$\\binom n k = \\prod_{p\\in\\mathcal P, p\\le n}p^{\\nu_p \\binom n k}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def binomial4(n, k):\n",
    "    ps = factor_binomial(n, k)\n",
    "    b = 1\n",
    "    for p, m in ps:\n",
    "        b = b * p ** m\n",
    "    return b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 11\n",
    "for n in range(N):\n",
    "    print(fmt(n + 1, 4) % tuple([binomial4(n, k) for k in range(n + 1)]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(binomial4(1000, 353))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ici encore, résultat en zéro seconde. L'auteur de l'article avait dû patienter 8 secondes pour effectuer le produit."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.8 Complexité"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comparons les temps d'exécution de `binomial2` (en noir) et `binomial4` (en rouge)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_binomial4(N):\n",
    "    s2 = []\n",
    "    s4 = []\n",
    "    ns = list(range(N, 11 * N, N))\n",
    "    for n in ns:\n",
    "        t1 = time.time()\n",
    "        x = binomial4(2 * n, n)\n",
    "        t2 = time.time()\n",
    "        s4.append(t2 - t1)\n",
    "        t1 = time.time()\n",
    "        x = binomial2(2 * n, n)\n",
    "        t2 = time.time()\n",
    "        s2.append(t2 - t1)\n",
    "    plt.plot(ns, s2, '-ok')\n",
    "    plt.plot(ns, s4, '-or')\n",
    "    plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_binomial4(2000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pour de petites valeurs de $n$, `binomial2`est le vainqueur. Mais en gros pour $n\\ge 8000$, c'est `binomial4`qui est le meilleur. \n",
    "\n",
    "Reprenons maintenant pour `binomial4` les tests que nous avions faits pour `binomial2`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_binomial4bis(N):\n",
    "    s = []\n",
    "    ns = list(range(N, 11 * N, N))\n",
    "    for n in ns:\n",
    "        t1 = time.time()\n",
    "        x = binomial4(2 * n, n)\n",
    "        t2 = time.time()\n",
    "        s.append(t2 - t1)\n",
    "    plt.loglog(ns, s, '-ok')\n",
    "    plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_binomial4bis(3000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La courbe est quasiment une droite. Si nous appelons $T_n$ le temps mis pour calculer $\\binom {2n}n$ on a donc à peu près\n",
    "\n",
    "$$\\ln T_n = d + \\beta \\log n$$\n",
    "\n",
    "où $d$ et $\\beta$ sont deux constantes. $\\beta$ est la pente de la droite."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t1 = time.time()\n",
    "na = 50000\n",
    "x = binomial4(2 * na, na)\n",
    "t2 = time.time()\n",
    "ta = t2 - t1\n",
    "\n",
    "t1 = time.time()\n",
    "nb = 100000\n",
    "x = binomial4(2 * nb, nb)\n",
    "t2 = time.time()\n",
    "\n",
    "tb = t2 - t1\n",
    "\n",
    "beta = (math.log(tb) - math.log(ta)) / (math.log(nb) - math.log(na))\n",
    "print(beta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rappelons-nous que pour la fonction `binomial2` nous avions une pente supérieure à 2. Obtenir la constante $d$ est maintenant facile."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = math.log(tb) - beta * math.log(nb)\n",
    "print(d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Maintenant, $T_n = e^d n^\\beta$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def temps_binomial4(n):\n",
    "    k = math.exp(d)\n",
    "    return k * n ** beta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "temps_binomial4(1000000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cette estimation nous rend confiants. Nous devrions pouvoir enfin calculer $\\binom {2000000}{1000000}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t1 = time.time()\n",
    "x = binomial4(2000000, 1000000)\n",
    "t2 = time.time()\n",
    "print('temps : %fs'% (t2 - t1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "C'est un peu plus que le temps prévu, mais nous avons gagné notre pari."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Références"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- P. Goetgheluck, Computing Binomial Coefficients - The American Mathematical Monthly, Vol. 94, No. 4 (Apr. 1987), pp. 360-365."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}