{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Partitions d'un entier\n",
    "\n",
    "Marc Lorenzi - 21 août 2018"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from math import *\n",
    "import sys\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (8, 8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dans le film _The Man Who Knew Infinity_, le mathématicien S. Ramanujan fait le pari avec P.A. Macmahon qu'il sera capable en quelques mois de calculer $P(200)$, le nombre de partitions de l'entier 200. La scène est probablement inventée. Ce qui est vrai, en revanche, c'est que Percy Alexander Macmahon utilisa une formule due à Euler (nous la verrons dans ce notebook) pour calculer les valeurs de $P(n)$ jusqu'à $n=200$. Cela se passait au tout début du XXème siècle. Et nous, pouvons-nous, avec Python, faire ce même calcul en moins de 1 seconde ? Pouvons-nous calculer $P(100000)$ en moins d'une minute ? Et au fait, c'est quoi une partition ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1 C'est quoi une partition ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Définition__ : Soit $n$ un entier naturel non nul. Une __partition__ de $n$ est une suite (finie) d'entiers naturels non nuls $(k_1,k_2,\\ldots,k_p)$ (où $p\\ge 1$) telle que $\\sum_{i=1}^pk_i=n$. Toute permutation d'une telle suite étant encore une partition de $n$, nous supposerons sans perte de généralité que $k_1\\ge k_2\\ge\\ldots\\ge k_p$. \n",
    "\n",
    "__Conventions__ : \n",
    "- Nous conviendrons que l'entier 0 possède une unique partition : la suite vide.\n",
    "- Nous conviendrons aussi que les entiers négatifs n'ont pas de partition."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Notation__ : Pour tout $n\\in\\mathbb Z$, on note $P(n)$ le nombre de partitions de l'entier $n$.\n",
    "\n",
    "__Cas particuliers__ : On a donc $P(0)=1$ et pour tout $n<0$, $P(n)=0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Exemples__ : \n",
    "- 1 n'a qu'une partition qui est $(1)$. On a $P(1)=1$.\n",
    "- Les partitions de 2 sont $(2)$ et $(1,1)$. On a $P(2)=2$.\n",
    "- Les partitions de 3 sont $(3)$, $(2,1)$ et $(1,1,1)$. On a $P(3)=3$.\n",
    "- Les partitions de 4 sont $(4)$, $(3,1)$, $(2,2)$, $(2,1,1)$ et $(1,1,1,1)$. On a $P(4)=5$.\n",
    "- Les partitions de 5 sont $(5)$, $(4,1)$, $(3,1,1)$, $(2,1,1,1)$, $(1,1,1,1,1)$, $(2,2,1)$ et $(3,2)$. On a $P(5)=7$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Exercice__ : Déterminer les partitions de 6."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 Qu'allons-nous faire ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La théorie des partitions est une théorie très riche. Il y a énormément de choses à dire sur le sujet. Nous allons rester très modestes : nous allons uniquement dans ce qui suit décrire trois algorithmes de calcul de $P(n)$. \n",
    "\n",
    "- Le premier est naïf et terriblement inefficace, il permet d'envisager le calcul de $P(50)$, et pas beaucoup plus.\n",
    "- Le second est une adaptation du premier. Il permet d'aller jusqu'à environ $P(1000)$. \n",
    "- Le troisième utilise une formule due à Euler. Il nous donnera sereinement $P(100000)$.\n",
    "\n",
    "Nous n'atteindrons pas $P(1000000)$, mais nous donnerons à la fin du notebook un équivalent de $P(n)$ lorsque $n$ tend vers l'infini. Cela nous permettra d'avoir une idée de la _taille_ de ce nombre."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Calculer toutes les partitions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comment déterminer toutes les partitions de l'entier $n$ ? Étant donnée une partition $(k_1,k_2,\\ldots,k_p)$ de $n$, la suite $(k_1,k_2,\\ldots,k_{p-1})$ est une partition de l'entier $n-k_p$. Mais, attention,  cette partition de $n-k_p$ a une propriété supplémentaire : les nombres de cette partition sont supérieurs à $k_p$, puisque nous avons convenu d'écrire une partition dans l'ordre décroissant. \n",
    "\n",
    "Nous allons donc écrire une fonction `part(n, k)` qui renvoie la liste des partitions de $n$ dont tous les éléments sont supérieurs ou égaux à $k$. Bien entendu, le calcul des partitions de $n$ se fera en évaluant `part(n, 1)`.\n",
    "\n",
    "- Si $k > n$ c'est évident : il n'y a aucune telle partition.\n",
    "- Sinon, pour tous les entiers $j$ compris entre $k$ et $n-k$, on calcule les partitions de $n-j$ dont les éléments sont supérieurs ou égaux à $j$ (hypothèse de décroissance des suites) et à $k$ (uniquement des entiers supérieurs à $k$). On rajoute $j$ à la fin de ces partitions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def part(n, k=1):\n",
    "    if k > n: return []\n",
    "    else:\n",
    "        s = [[n]]\n",
    "        for j in range(k, n - k + 1):\n",
    "            s1 = part(n - j, max(j, k))\n",
    "            s1 = [t + [j] for t in s1]\n",
    "            s = s + s1\n",
    "        return s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "part(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(part(10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(part(10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notre fonction `part` a évidemment une complexité épouvantable. Nous verrons un peu plus loin dans ce notebook que le nombre de partitions de $n$ croît \"exponentiellement\" avec $n$. Toute fonction renvoyant les partitions de $n$ a donc une complexité exponentielle. Cessons donc de nous intéresser au calcul des partitions, et intéressons-nous uniquement au __nombre__ de partitions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Nombre de partitions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 Première tentative"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comment calculer le nombre de partitions d'un entier $n$ ? On peut évidemment calculer la liste de toutes les partitions, puis déterminer la longueur de cette liste comme nous l'avons fait pour $n=10$ un peu plus haut. Cet algorithme gaspille un espace énorme : on n'a pas besoin de stocker les partitions, on veut juste les compter. Adaptons donc la fonction `part` ci-dessus pour qu'elle compte au lieu de stocker ... la fonction `npart1` ci-dessous prend deux paramètres $n$ et $k$. Elle renvoie le nombre de partitions de $n$ dont tous les éléments sont supérieurs ou égaux à $k$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def npart1(n, k=1):\n",
    "    if k > n: return 0\n",
    "    else:\n",
    "        s = 1\n",
    "        for j in range(k, n - k + 1):\n",
    "            s = s + npart1(n - j, max(j, k))\n",
    "        return s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "npart1(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voici les premières valeurs de $P(n)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for k in range(21):\n",
    "    print('%5d%5d' % (k, npart1(k)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "npart1(50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La complexité de `npart1` est exponentielle en $n$. Inutile, donc, d'essayer avec celle-ci d'obtenir le nombre de partitions de, mettons, 100. Peut-on faire mieux ? Oui, heureusement."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Exercice__ : quel est le plus grand entier $n$ pour lequel (sur votre machine) la fonction `npart1` renvoie le nombre de partitions en moins d'une minute ? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Une fonction de complexité polynomiale cubique en $n$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Un inconvénient majeur de la fonction `npart` est qu'elle fait des appels récursifs un grand nombre de fois avec les mêmes valeurs du paramètre. Bref, elle calcule, recalcule et rerecalcule tout le temps la même chose. Pour nous en convaincre, la fonction `testnpart1 ci-dessous` prend en paramètres deux entiers $n$ et $k$ et une matrice $n$. Après exécution de `testnpart1(n, m, k)`, $m_{ij}$ est égal au nombre de fois que `npart1(i, j)` a été calculé lors de l'appel `npart1(n, k)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def testnpart1(n, m, k=1):\n",
    "    if k > n: return 0\n",
    "    else:\n",
    "        m[n][k] += 1\n",
    "        s = 1\n",
    "        for j in range(k, n - k + 1):\n",
    "            s = s + testnpart1(n - j, m, max(j, k))\n",
    "        return s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 40\n",
    "m = (n + 1) * [None]\n",
    "for k in range(n + 1):\n",
    "    m[k] = (n + 1) * [0]\n",
    "testnpart1(n, m)\n",
    "\n",
    "# Afficher joliment la partie intéressante de m\n",
    "sys.stdout.write('    ')\n",
    "for j in range(n // 2 + 1):\n",
    "    sys.stdout.write('%4d' % (j))\n",
    "print()\n",
    "for i in range(n + 1):\n",
    "    sys.stdout.write('%4d' % (i))\n",
    "    for j in range(n // 2 + 1):\n",
    "        if m[i][j] == 0: sys.stdout.write('%4s' % u'\\u2022')\n",
    "        else: sys.stdout.write('%4d' % (m[i][j]))\n",
    "    print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Par exemple, lors de l'appel `npart1(40)`, `npart1(10, 6)` est calculé 532 fois. Essayez avec $n=50$, certains calculs sont effectués ... des milliers de fois !"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " Pour éviter cela il suffit de stocker tous les résultats dans une matrice $P$. Pour $1\\le i\\le n$ et $1\\le j\\le n$, $P_{ij}$ est le nombre de partitions de l'entier $i$ ne contenant que des entiers supérieurs ou égaux à $j$. On peut même éliminer tout recours à la récursivité."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def npart2(n, mat=False):\n",
    "    P = (n + 1) * [None]\n",
    "    for k in range(n + 1):\n",
    "        P[k] = (n + 1) * [0]\n",
    "    for j in range(n + 1): P[0][j] = 1\n",
    "    for i in range(1, n + 1):\n",
    "        for j in range(1, i + 1):\n",
    "            P[i][j] = 1\n",
    "            for k in range(j, i - j + 1):\n",
    "                P[i][j] = P[i][j] + P[i - k][max(k, j)]\n",
    "    if mat: return P\n",
    "    else: return P[n][1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "npart2(10, mat=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "npart2(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quelle est la complexité de notre fonction `npart2` ? Sans entrer dans les détails, la triple boucle `for` en $i,j,k$ annonce une complexité en $\\mathcal O(n^3)$. Nous avons donc maintenant une complexité polynomiale en $n$. Bien entendu l'exposant $n^3$ nous dit qu'il faudra rester raisonnables. Mais le calcul de $P(100)$ est maintenant immédiat."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(npart2(100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Et $P(1000)$ ? Eh bien, $1000^3=10^9$, c'est donc faisable en étant un tout petit peu patient (une quarantaine de secondes sur ma machine)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(npart2(1000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Le calcul de $P(10000)$ prendrait $10^3=1000$ fois plus de temps. Il faut trouver mieux que notre algorithme cubique."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3 Peut-on faire mieux ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Oui, on peut faire mieux. La formule ci-dessous est non triviale, elle a été énoncée par Euler.\n",
    "\n",
    "__Théorème__ : Pour tout $n\\ge 1$, on a $P(n)=\\sum_{k=1}^\\infty (-1)^{k-1}[P(n-\\frac 1 2 k(3k-1))+P(n-\\frac 1 2 k(3k+1)]$\n",
    "\n",
    "__Remarque__ : La somme ci-dessus est en réalité finie puisque son terme général est nul dès que $\\frac 1 2 k(3k-1)>n$. Cela nous donne environ $2\\sqrt{\\frac{2n}3}$ termes, ce qui est donc bien mieux que les sommes que nous avons considérées jusqu'à présent. Gardez cette valeur dans un coin de votre esprit, en mathématiques tout fait partie d'un grand plan :-). \n",
    "\n",
    "__Démonstration__ : Je ne ferai pas la démonstration de la formule d'Euler (encore une !!!) mais je ne résiste pas à la tentation d'évoquer sa provenance. On peut montrer que pour tout réel $x$ tel que $|x|<1$, on a la merveilleuse égalité\n",
    "\n",
    "$$\\sum_{n=0}^\\infty P(n)x^n=\\frac 1{(1-x)(1-x^2)(1-x^3)...}=\\prod_{k=1}^\\infty \\frac 1 {1-x^k}$$\n",
    "\n",
    "La série du membre de gauche est la __série génératrice__ de la suite $(P(n))_{n\\ge 0}$. Le membre de droite est une fonction de $x$, écrite sous forme d'un __produit infini__. Cette fonction de $x$ contient en essence les valeurs de $P(n)$ pour __tous__ les entiers $n$. En développant le produit infini $\\prod_{k=1}^\\infty (1-x^k)$, puis en multipliant le résultat obtenu par $\\sum_{n=0}^\\infty P(n)x^n$, on obtient, en vertu de la merveilleuse égalité ... 1. La formule qui nous intéresse en découle."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voici maintenant la fonction `npart3`. Elle prend un entier $n$ en paramètre. Puis elle calcule $P(i)$ pour $i=1,2,\\ldots, n$ en utilisant la formule d'Euler. Nous suivons donc la méthode employée par Macmahon pour calculer $P(200)$ à la main. Les valeurs des $P(i)$ sont stockées au fur et à mesure dans une liste $P$.\n",
    "\n",
    "Je n'entrerai pas dans les détails, mais une écriture soigneuse de l'algorithme permet le calcul des $P(k)$, $1\\le k\\le n$ en faisant $O(n^2)$ opérations élémentaires. On peut vérifier que la fonction `npart3` effectue $O(n^{\\frac 3 2})$ additions sur des entiers. Bien évidemment, les entiers en question peuvent être très grands et une simple addition peut devenir coûteuse. Une estimation plus fine de la taille des entiers mis en jeu donnerait (ou pas :-)) la complexité $O(n^2)$ espérée. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def npart3(n, mat=False):\n",
    "    P = (n+1) * [0]\n",
    "    P[0] = 1\n",
    "    for i in range(1, n + 1): \n",
    "        k = 1\n",
    "        s = 0\n",
    "        while True:\n",
    "            u = i - k * (3 * k - 1) // 2\n",
    "            v = i - k * (3 * k + 1) // 2\n",
    "            if u < 0 and v < 0: break\n",
    "            if u >= 0: s = s + (-1) ** (k - 1) * P[u]\n",
    "            if v >= 0: s = s + (-1) ** (k - 1) * P[v]\n",
    "            k = k + 1\n",
    "        P[i] = s\n",
    "    if mat: return P\n",
    "    else: return P[n]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Maintenant, le calcul de $P(1000)$ est immédiat."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "npart3(1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Le calcul de $P(10000)$ est quasi-immédiat."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "npart3(10000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Et $P(100000)$ ? Aussi, mais soyez un petit peu patients si vous voulez avoir la réponse (environ 50 secondes sur ma machine)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "npart3(100000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Exercice__ : Faites comme Macmahon et calculez _à la main_ $P(k)$ pour $k=1,2,3,\\ldots$ avec l'algorithme utilisé par `npart3`. Où en êtes-vous au bout d'une heure ? Dans le film cité plus haut, il est dit que Macmahon a obtenu $P(100)$ en un mois :-). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.4 La \"complexité expérimentale\" de `npart3`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nous avons dit un peu plus haut que la fonction `npart3` effectuait $O(n^{\\frac 3 2})$ additions pour calculer les valeurs de $P(k)$, $k=1,\\ldots,n$. Nous avons également remarqué que, ces additions opérant sur de grands entiers, la complexité réelle de notre fonction était plutôt en $O(n^2)$. Qu'en est-il en réalité ?\n",
    "\n",
    "Calculons le temps mis par notre machine pour effectuer ces calculs, ceci pour $n=1000, 2000,\\ldots,30000$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import timeit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Attention, l'évaluation de la cellule ci-dessous prend du temps. Allez prendre un café."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tms = [timeit.timeit('npart3(1000*%d)' % (n), setup='from __main__ import npart3', number=1) for n in range(1, 31)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(tms)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Imaginons maintenant, en appelant $t_n$ le temps mis pour calculer $P(1),\\ldots,P(n)$, que $t_n\\sim Cn^\\alpha$ où $C$ et $\\alpha$ sont deux réels strictement positifs ? $C$ dépend évidemment de la machine choisie, mais que vaut $\\alpha$ ? Remarquons que $\\ln t_n \\sim \\ln C + \\alpha\\ln n\\sim \\alpha\\ln n$, le graphe de $t_n$ en fonction de $n$ ressemblera donc fort à une droite."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot([log(k) for k in range(1, 31)], [log(t) for t in tms])\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Effectivement, cela ressemble fort à une droite. Quelle est sa pente, $\\alpha$ ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(log(tms[29])-log(tms[10])) / (log(30000) - log(11000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On obtient $\\alpha\\simeq \\frac 3 2$. Surprenant ... que peut-on en conclure ? Que les effets \"grands entiers\" ne se sont pas encore fait sentir ? Si nous étions courageux, nous calculerions les temps mis par notre fonction pour des entiers $n$ plus grands. Il y a fort à parier que notre \"constante\" $\\alpha$ serait plus grande que $\\frac 3 2$ ...\n",
    "\n",
    "Laissons tout cela au conditionnel, c'est l'heure de l'apéro :-)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.5 Peut-on faire encore mieux ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Oui, par exemple en utilisant des transformées de Fourier, mais ceci est une autre histoire que je n'aborderai pas ici. On peut calculer les $P(k)$, $1\\le k\\le n$ avec une complexité en $O(n^{\\frac 3 2}\\log^2 n)$ opérations élémentaires. Pour $n=1000000$ et le logarithme décimal, la quantité dans le grand $O$ est égale à 36 milliards, un nombre d'opérations grand mais pas inatteignable. On peut également paralléliser les calculs, c'est à dire utiliser plusieurs machines qui effectuent chacune une partie des calculs, puis recombiner le tout. De telles techniques permettent d'atteindre la valeur de $P(n)$ pour des valeurs de $n$ de l'ordre ... du milliard."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Un équivalent de $P(n)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1 L'équivalent de Hardy et Ramanujan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "J'ai annoncé au début de ce notebook que $P(n)$ croît \"exponentiellement\" avec $n$. Qu'en est-il exactement ? On possède un équivalent de $P(n)$ lorsque $n$ tend vers l'infini. Ce résultat est dû à G.H. Hardy et S. Ramanujan. Le voici sans preuve :"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Théorème__ : $P(n)\\sim \\frac 1 {4 n\\sqrt 3}\\exp(\\pi \\sqrt{\\frac{2n}3})$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Remarque__ : Ce n'est pas la première fois que nous rencontrons $\\sqrt{\\frac{2n}3}$ ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Remarque__ : Cet équivalent est le premier terme d'un développement asymptotique de $P(n)$ (Hardy, Ramanujan, Rademacher) dont le reste tend vers 0 lorsque le nombre de termes du développement tend vers l'infini. On dispose de plus d'un majorant de l'erreur commise. En utilisant ce développement asymptotique on peut donc calculer $P(n)$ de façon exacte puisque ce nombre est ... un entier ! Des algorithmes efficaces de calcul de $P(n)$ ont été développés par cette technique."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Exercice__ : Que signifient le G.H. et le S. de G.H. Hardy et S. Ramanujan ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `approx_rama` calcule cet équivalent. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def approx_rama(n):\n",
    "    return 1 / (4 * n * sqrt(3)) * exp(pi * sqrt(2 * n / 3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Prenons par exemple $n=200$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "approx_rama(200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "float(npart3(200))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "approx_rama(200) / npart3(200)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "L'approximation de Hardy et Ramanujan nous donne donc $P(200)$ avec une erreur de l'ordre de 3%, ce qui est remarquable. Si vous ne trouvez pas cela remarquable, refaite l'exercice \"sur les pas de Macmahon\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Traçons sur un même graphe $P(n)$ et son approximation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 100\n",
    "rg = range(1, n + 1)\n",
    "P = npart3(n, mat=True)\n",
    "Q = [approx_rama(k) for k in rg]\n",
    "plt.semilogy(rg, P[1:],'k')\n",
    "plt.semilogy(rg, Q, 'r')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Les deux courbes sont quasi-superposées. Autre vision : le tracé de l'erreur relative."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 100\n",
    "rg = range(2, n + 1)\n",
    "P = npart3(n, mat=True)\n",
    "Q = [approx_rama(k) / P[k] for k in rg]\n",
    "plt.plot(rg, Q, 'k')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.2 La taille de $P(n)$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soit $N$ un entier non nul. Le nombre de chiffres en base 10 de l'entier $n$ est $\\lfloor\\log n\\rfloor + 1$, où $\\log$ désigne le logarithme en base 10."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def taille(n): return floor(log(n) / log(10)) + 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "taille(12345)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Il est donc immédiat d'avoir le nombre de chiffres de $P(n)$ lorsqu'on peut calculer $P(n)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "taille(npart3(10000))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def approx_taille(n):\n",
    "    return floor((pi * sqrt(2 * n / 3) - log(4 * n * sqrt(3))) / log(10)) + 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "approx_taille(10000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quel est (environ, ou peut-être exactement, qui sait ?) le nombre de chiffres de $P(1000000)$ ? De $P(1000000000)$ ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "approx_taille(10 ** 6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "approx_taille(10 ** 9)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notons $\\varphi(n)=\\frac 1 {4 n\\sqrt 3}\\exp(\\pi \\sqrt{\\frac{2n}3})$. On a $\\frac {P(n)}{\\varphi(n)}\\to 1$ lorsque $n$ tend vers l'infini. En passant aux logarithmes, il vient donc que $\\log P(n)-\\log\\varphi(n)\\to 0$ lorsque $n$ tend vers l'infini. La fonction `approx_taille` renvoie donc la valeur à peu près exacte (à 1 près) du nombre de chiffres de $P(n)$, pour tout $n$ assez grand."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 500\n",
    "rg = range(2, n + 1)\n",
    "P = npart3(n, mat=True)\n",
    "Q = [taille(approx_rama(k)) - taille(P[k]) for k in rg]\n",
    "plt.plot(rg, Q, 'k')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La taille de `approx_rama(n)` serait-elle la même que celle de $P(n)$ à 1 près, pour __tout__ n ? Pourquoi ces pics par ci par là ? Il y en a un, par exemple, pour $n=60$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "npart3(60)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "approx_rama(60)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ah ben oui. $P(60)$ est tout juste inférieur à une puissance de 10. Je n'irai pas plus loin, bien conscient de n'avoir fait qu'effleurer un sujet dont l'étendue est monumentale."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Exercice__ : il y a 7 pics sur le graphe ci-dessous. Le premier d'entre eux se situe à $n=60$. Trouvez la position exacte des 6 autres."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Références"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- The Man who Knew Infinity (2015). Réalisateur : Matt Brown. Acteurs : Dev Patel, Jeremy Irons.\n",
    "- An Introduction to the Theory of Numbers, G.H. Hardy et E.M. Wright (6ème édition, 2008). Chapitre XIX : Partitions.\n",
    "\n",
    "Et, bien entendu, une recherche sur Internet avec les mots clés \"integer partition\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}