Floating Point Arithmetic.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Floating Point Arithmetic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* Finite computing and memory capacity\n",
    "* Need to find a way to represent floating point numbers\n",
    "* IEEE 754 defines common standard on modern computers to represent floating\n",
    "point numbers.\n",
    "* Two most commonly used types: `IEEE double precision` and `IEEE single precision`.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Floating point types in Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The Numpy module defines convenient ways to query properties of floating point numbers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np # import the numpy extension module\n",
    "                   # and call it np as short form"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The data type names in `Numpy` for floating point types are:\n",
    "\n",
    "* `IEEE double precision`: np.float, np.double, np.float64\n",
    "* `IEEE single precision`: np.single, np.float32\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Let us query some properties of these numbers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "double_precision_info = np.finfo(np.float64)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The biggest and smallest (by absolute value) normalized floating point numbers are"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.7976931348623157e+308"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "double_precision_info.max"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.2250738585072014e-308"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "double_precision_info.tiny"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Floating point numbers can not be arbitrarily close to each other. There is a smallest relative distance, which we define shortly. It is given as"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.220446049250313e-16"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "double_precision_info.eps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "This leads to a limited precision of floating point numbers. The approximate relative precision is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "15"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "double_precision_info.precision"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Task: What are the values for single precision arithmetic?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[3.4028235e+38, 1.1754944e-38, 1.1920929e-07, 6]\n"
     ]
    }
   ],
   "source": [
    "single_precision_info = np.finfo(np.float32)\n",
    "print([single_precision_info.max,\n",
    "       single_precision_info.tiny,\n",
    "       single_precision_info.eps,\n",
    "       single_precision_info.precision])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Definition of floating point numbers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The set of floating point numbers is defined as follows:\n",
    "\n",
    "$$\n",
    "\\mathcal{F} = \\left\\{(-1)^s\\cdot b^e \\cdot \\frac{m}{b^{p-1}} :\\right.\n",
    "\\left. s = 0,1; e_{min}\\leq e \\leq e_{max}; b^{p-1}\\leq m\\leq b^{p}-1\\right\\}.\n",
    "$$\n",
    "\n",
    "* `IEEE double precision`: $e_{min} = -1022, e_{max} = 1023, p = 53$\n",
    "* `IEEE single precision`: $e_{min} = -126, e_{max} = 127, p = 24$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Typically, b = 2. For Mantissa have:\n",
    "\n",
    "$$\n",
    "\\frac{m}{b^{p-1}} = 1, 1+b^{1-p}, 1+2b^{1-p}, \\dots, b-b^{1-p}\n",
    "$$\n",
    "\n",
    "$\\rightarrow$ Distance of neighbouring floats is $2^e b^{1-p}$.\n",
    "\n",
    "$\\epsilon_{rel} = b^{1-p}$ is smallest number such that\n",
    "$$\n",
    "1 + \\epsilon_{rel} \\neq 1.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0000000000000002"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1 + double_precision_info.eps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1 + .25 * double_precision_info.eps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Approximating numbers in floating point arithmetic\n",
    "Let $x\\in\\mathbb{R}$, $b^{e_{min}}\\leq x < b^{e_{max}+1}$. Define $\\epsilon_{mach}:=\\frac{1}{2}b^{1-p}$\n",
    "\n",
    "There exists $x'\\in \\mathcal{F}$ such that $|x-x'|\\leq\\epsilon_{mach}|x|$.\n",
    "\n",
    "**$\\epsilon_{mach}$ is relative distance to the next floating point number in $\\mathcal{F}$.**\n",
    "\n",
    "Define the projection\n",
    "$$fl:fl(x)\\rightarrow x',$$\n",
    "where $x'$ is the closest floating point number in $\\mathcal{F}$. \n",
    "\n",
    "It follows that $fl(x)=x*(1+\\epsilon)$ for some $|\\epsilon|\\leq \\epsilon_{mach}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Fundamental Axiom of Floating Point Arithmetic**\n",
    "Define $x\\odot y = fl(x \\cdot y)$, where $\\cdot$ is one of $+,-,\\times,\\div$. Then for all $x,y\\in\\mathcal{F}$ there exists $\\epsilon$ with $|\\epsilon| \\leq \\epsilon_{mach}$ such that\n",
    "$$ x\\odot y = (x \\cdot y)(1+\\epsilon).$$\n",
    "\n",
    "Most modern architectures guarantee this property."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Special symbols in floating point arithmetic\n",
    "\n",
    "In addition to numbers several other important symbols are defined in the floating point standard. The most important are:\n",
    "\n",
    "* NaN: Not a number\n",
    "* $\\pm$inf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "nan\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/nbuser/anaconda3_501/lib/python3.6/site-packages/ipykernel/__main__.py:3: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  app.launch_new_instance()\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "a = np.inf\n",
    "b = np.float64(0) / 0\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python [conda env:dev]",
   "language": "python",
   "name": "conda-env-dev-py"
  },
  "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}