{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e935c4e5-efe0-4630-b245-842ad4512b69",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Représentation des nombres"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0eeeed8f-1558-4035-8711-a097786330a9",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# importing libraries used in this book\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01087337-5e01-4deb-a8c8-83fe147d63f3",
   "metadata": {},
   "source": [
    "## Représentation de nombres : Exemple $e$\n",
    "\n",
    "L'expression suivante du nombre d'Euler $e$ est connue :\n",
    "$$e = \\lim_{n\\to\\infty}\\left( 1 + \\frac{1}{n} \\right)^{n}.$$\n",
    " On s'attend donc à ce que $e_n = \\left(1 + \\frac{1}{n} \\right)^{n}$ donne des approximations de plus en plus bonnes de $e$\n",
    " En arithmétique exacte, c'est effectivement le cas. Sur l'ordinateur, la suite \\emph{calculée} $\\hat e_n$ se comporte tout à fait différemment :\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d3ef1dd-bc69-41fd-b6e6-1061b0f57b94",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Approximation of e, numpy package needed to include e\n",
    "import numpy as np\n",
    "print('10^i \\t\\t e_n \\t\\t\\t e_n - e')\n",
    "for i in range(1,16):\n",
    "    n = 10.0 ** i; en = (1 + 1/n) ** n\n",
    "    print('10^%2d \\t %20.15f \\t %20.15f' % (i,en,en-np.e))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f962ce28-6a55-4207-93d9-eaa5a7849bcd",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Représentation de nombres : Exemple $e^x$\n",
    "\n",
    "La série de Taylor pour la fonction exponentielle converge pour tout $x \\in \\mathbb R$ :\n",
    "$$e^x = \\sum\\limits_{k=0}^{\\infty}\\frac{x^{k}}{k!} = \n",
    "        1 + x + \\frac{x^{2}}{2} + \\frac{x^{3}}{6} + \\frac{x^{4}}{24} + \\dots.$$\n",
    "L'ordinateur peut seulement calculer la série partielle\n",
    "$$s_i(x) = \\sum\\limits_{k=0}^{i}\\frac{x^{k}}{k!}.$$\n",
    "Le reste de Taylor est donné par\n",
    "$e^x - s_i(x)  = \\frac{e^\\xi x^{i+1}}{(i+1)!}\n",
    "\\text{ pour un } \\xi \\in \\mathbb R \\text{ avec }0<|\\xi|<|x|$$\n",
    "\n",
    "Si on choisi une tolerance $\\mathsf{tol}>0$ et $i$ tel que\n",
    "$|x|^{i+1} / (i+1)! \\le \\mathsf{tol} \\cdot s_i(x)$, pour $x<0$ on obtient\n",
    "$$|e^x - s_i(x)| \\le   \\frac{|x|^{i+1}}{(i+1)!}  \\le \\mathsf{tol} \\cdot s_i(x) \\approx \\mathsf{tol} \\cdot e^x.$$\n",
    "Au même temps,  **l'erreur rélative** $|e^x - s_i(x)| / |e^x|$ est bornée par  $\\mathsf{tol}$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4a7c17a-00fa-47c8-9f01-26bf50abdb73",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def expeval(x, tol):\n",
    "    #Approximation of e^x\n",
    "    s = 1; k = 1\n",
    "    term = 1\n",
    "    while (abs(term)>tol*abs(s)):\n",
    "        term = term * x / k\n",
    "        s = s + term ; k = k + 1\n",
    "    return s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c226fa49-0d5f-4458-aa38-97be980ab158",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# reproduisez ici le tableau\n",
    "# $x$ & Computed $\\hat s_i(x)$ & $\\exp(x)$ & $\\frac{|\\exp(x)-\\hat s_i(x)|}{\\exp(x)}$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be14a2b4-f962-4fa2-9f5b-27fc6d640d81",
   "metadata": {
    "tags": []
   },
   "source": [
    "## IEEE 754 in Python\n",
    "En Python, toutes les opérations sur les nombres réels sont exécutées par défaut en double précision.\n",
    "Les variables en simple précision sont générées par la commande ```numpy.float32()```."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e5ad2fc5-856c-4ab7-950e-21401b65924e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.float_info.min             # 2.2251e-308"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "def89b00-2173-4f20-b94a-f5eb7cfeb106",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sys.float_info.max             # 1.7977e+308"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17167230-0a86-4b90-9317-6ed5d1ca45a8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "1 / 0                          # Divide by zero error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03aeeee4-074f-4f34-b9e3-5bff9a70c445",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "3 * float('inf')               # np.inf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eddbb58b-ed38-4627-a595-62b55435243a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "-1 / 0                         # Divide by zero error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da4157f4-d940-42f1-983e-bc6e0ea25afb",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "0 / 0                          # Divide by zero error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7a75c76-1132-487a-8181-f92be0492246",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "float('inf') - float('inf')    # nan"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76ad9faf-2be2-4b23-923c-c636cdc8cfd3",
   "metadata": {},
   "source": [
    "De manière assez surprenante, et non en accord avec le standard IEEE 754 https://wusun.name/blog/2017-12-18-python-zerodiv/, Python renvoit une erreur lorsqu'une division par zéro se produit, au lieu de retourner un $\\infty$ signé. Pour éviter cela, on peut utiliser ``float64`` de ``numpy``. Par exemple, ``1/np.float64(0)`` renvoie ``inf`` (comme il se doit).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e27340b1-213d-4b4f-bef7-6ecec40c37c5",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Arrondis\n",
    "\n",
    "### Sommes\n",
    "\n",
    "Comme exemple, regardons ce qu'il se pase avec de ```float16``` et une somme de nombre de plus en plus petits.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68662282-a72e-4791-8d3b-f5d16f7d0ba9",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# see what happens with n = 2, 10, 100, 1000, 10000\n",
    "n = 2 # COMPLETE HERE\n",
    "\n",
    "# t and y are 64 bit float (double) numbers\n",
    "t = np.sin(np.linspace(0,1,n))\n",
    "y = np.linspace(0,1,n) \n",
    "\n",
    "y[1::2] = -1 * y[1::2] # odd numbers set to negatives\n",
    "\n",
    "# reduced accuracy\n",
    "x= np.float16(y) \n",
    "r= np.float16(t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b24996b9-c181-4492-950f-e42e7ce7a141",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Use the sum function\n",
    "print(sum(y*t))\n",
    "print(sum(x*r))\n",
    "print(np.float16(sum(x*r)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "752a87e2-a28d-4377-8ff8-ecb8a2dd2336",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Use the sum function, reverse order\n",
    "# To reverse the order of an array, use y[::-1]\n",
    "\n",
    "# COMPLETE HERE\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f4907c6-3166-49d3-b800-52473aada8df",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# loop for the sum from the first to last, in float16\n",
    "s = np.float16(0);\n",
    "for k in range(n) :\n",
    "    s = s + x[k]*r[k]\n",
    "print(s)\n",
    "\n",
    "# loop for the sum from the first to last\n",
    "s = np.float16(0);\n",
    "for k in range(n) :\n",
    "    s = s + x[n-k-1]*r[n-k-1]\n",
    "print(s)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4fee6e0b-91b5-49a3-816f-a5856452a720",
   "metadata": {},
   "source": [
    "## Cancellation, exemple\n",
    "\n",
    "Considérons\n",
    "$$f(x) = \\frac{1- \\cos(x)}{x^2} \\;, \\quad x = b^{-k}, k=5,...,9.$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2830e01-006e-4c0b-bd73-541b54db4df7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "del x\n",
    "def f (x):\n",
    "    return (1 - np.cos(x))/ (x*x)\n",
    "\n",
    "# Try with b = 2,3,4,10\n",
    "b = 10\n",
    "\n",
    "k=np.array(range(6,20))\n",
    "x = 1/b**k\n",
    "\n",
    "# Uncomment the lines below\n",
    "# print('k \\t', k)\n",
    "# print('x=10^-k ',x)\n",
    "# print('f(x) \\t', f(x))\n",
    "\n",
    "## Plot a loglog graph with x and f(x)\n",
    "# USE plt.loglog(n,errN, 'b')\n",
    "\n",
    "# Labels\n",
    "plt.xlabel('x'); plt.ylabel('f(x)'); plt.title('f(x) with respect to x')\n",
    "plt.legend(['error'])\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e17b172b-a30a-4f19-8b06-d8305180ee4d",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Exemple, différences finies\n",
    "\n",
    "Considérons $f(x) = e^x$, $x_0=0$ et \n",
    "$$f^\\prime(x_0) = \\lim_{h\\to 0} \\frac{f(x_0+h) - f(x_0)}{h}.$$\n",
    "\n",
    "Attente : L'approximation de $f^\\prime(x_0)$ par un quotient de différences finies tend à s'améliorer lorsque $h$ s'approche de $0$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8b79693e-5af0-4732-8481-3e1d7ecf4786",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "\n",
    "def f(x) :\n",
    "    return np.exp(x)\n",
    "def df(x) :\n",
    "    return np.exp(x)\n",
    "x0 = 0\n",
    "\n",
    "# prenez h = 10{-k} avec k=0,1,...,16\n",
    "k = np.array(range(0,17))\n",
    "h = 1/10**k\n",
    "\n",
    "# calculez la difference finie (f(x+h) - f(x))/h pour tout les h\n",
    "fprime = (f(h) - f(0))/h\n",
    "\n",
    "# calculez l'erreur avec la solution exacte errH = |df(x0) - ..|\n",
    "# COMPLETEZ ICI\n",
    "\n",
    "# dessinez l'erreur en log-log \n",
    "\n",
    "# Labels\n",
    "plt.xlabel('h'); plt.ylabel('ErrH'); plt.title('Error in finite difference')\n",
    "plt.legend(['error'])\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a45e6948-bd26-4ce0-bf6a-01aad2c77e4e",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Exemple, le nombre $e$\n",
    "\n",
    "$$e = \\lim_{n\\to \\infty} \\big( 1 + \\frac1n)^n$$\n",
    "\n",
    "Attente : plus $n$ est grand, plus $e_n = \\big( 1 + \\frac1n)^n$ s'approche de  $e$ \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4546b369-675b-4e76-8a0b-e3c70928be4a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def limN(n) :\n",
    "    return (1+1/n)**n\n",
    "\n",
    "# prenez n = 10{k} avec k=0,1,...,16\n",
    "k = np.array(range(0,17))\n",
    "n = 10**k\n",
    "\n",
    "# calculez (1+1/n)**n\n",
    "approxE = limN(n)\n",
    "\n",
    "# calculez l'erreur avec la solution exacte errN = |e - ..|\n",
    "# COPLETEZ ICI    \n",
    "\n",
    "\n",
    "# dessinez l'erreur en log-log \n",
    "\n",
    "# Labels\n",
    "plt.xlabel('n'); plt.ylabel('ErrN'); plt.title('Error in finite difference')\n",
    "plt.legend(['error'])\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e5e9797-8b87-4857-8873-64e8ca3e25f4",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}