material/computer-arithmetics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computer Arithmetics and Round-off Methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the ideal mathematical world, operations like $1+2=3$, $4\\times 3 = 12$, $(\\sqrt{2})^2 = 2$ are unambiguously defined, however, when one is representing numbers in a computer, this is no longer true. The main reason of this is the so-called *finite arithmetic*, what is the way in which a computer performs basic operations. Some features of *finite arithmetic* are stated below:\n",
    "\n",
    "- Only integer and rational numbers can be exactly represented.\n",
    "- The elements of the set in which arithmetic is performed is necessarily finite.\n",
    "- Any arithmetic operation between two or more numbers of this set should be another element of the set.\n",
    "- Non-representable numbers like irrational numbers are approximated to the closest rational number within the defined set.\n",
    "- Extremely large numbers produce overflows and extremely small numbers produce underflows, which are taken as null.\n",
    "- Operations over non-representable numbers are not exact.\n",
    "\n",
    "In spite of this, defining adequately the set of elements in which our computer will operate, round-off methods can be systematically neglected, yielding correct results within reasonable error margins. In some pathological cases, when massive iterations are required, these errors must be taken into account more seriously.\n",
    "\n",
    "- - - \n",
    "\n",
    "- [Binary machine numbers](#Binary-machine-numbers)\n",
    "    - [Single-precision numbers](#Single-precision-numbers)\n",
    "    - [Double-precision numbers](#Double-precision-numbers)\n",
    "- [Finite Arithmetic](#Finite-Arithmetic)\n",
    "    - [Addition](#Addition)\n",
    "    - [Multiplication](#Multiplication)\n",
    "\n",
    "- - - "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Binary machine numbers    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As everyone knows, the base of the modern computation is the binary numbers. The binary base or base-2 numeral system is the simplest one among the existing numeral bases. As every electronic devices are based on logic circuits (circuits operating with [logic gates](#LogicGates)), the implementation of a binary base is straightforward, besides, any other numeral system can be reduced to a binary representation.\n",
    "\n",
    "![LogicGates](http://www.ee.surrey.ac.uk/Projects/CAL/digital-logic/gatesfunc/graphics/symtab.gif)\n",
    "\n",
    "According to the standard [IEEE 754-2008](http://en.wikipedia.org/wiki/IEEE_floating_point), representation of real numbers can be done in several ways, [single-precision](http://en.wikipedia.org/wiki/Single-precision_floating-point_format) and [double precision](http://en.wikipedia.org/wiki/Double-precision_floating-point_format) are the most used ones."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Single-precision numbers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Single-precision numbers are used when one does not need very accurate results and/or need to save memory. These numbers are represented by a **32-bits** (Binary digIT) lenght binary number, where the real number is stored following the next rules:\n",
    "\n",
    "![32-bits](http://upload.wikimedia.org/wikipedia/commons/thumb/d/d2/Float_example.svg/590px-Float_example.svg.png)\n",
    "\n",
    "1. The fist digit (called *s*) indicates the sign of the number (s=0 means a positive number, s=1 a negative one).\n",
    "2. The next 8 bits represent the exponent of the number.\n",
    "3. The last 23 bits represent the fractional part of the number.\n",
    "\n",
    "The formula for recovering the real number is then given by:\n",
    "\n",
    "$$r = (-1)^s\\times \\left( 1 + \\sum_{i=1}^{23}b_{23-i}2^{-i} \\right)\\times 2^{e-127}$$\n",
    "\n",
    "where $s$ is the sign, $b_{23-i}$ the fraction bits and $e$ is given by:\n",
    "\n",
    "$$e = \\sum_{i=0}^7 b_{23+i}2^i$$\n",
    "\n",
    "Next, it is shown a little routine for calculating the value of the represented 32-bits number"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def number32( binary ):\n",
    "    #Inverting binary string\n",
    "    binary = binary[::-1]\n",
    "    #Decimal part\n",
    "    dec = 1\n",
    "    for i in xrange(1,24):\n",
    "        dec += int(binary[23-i])*2**-i\n",
    "    #Exponent part\n",
    "    exp = 0\n",
    "    for i in xrange(0,8):\n",
    "        exp += int(binary[23+i])*2**i\n",
    "    #Total number\n",
    "    number = (-1)**int(binary[31])*2**(exp-127)*dec\n",
    "    return number"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.15625"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "number32( \"00111110001000000000000000000000\" )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Single-precision system can represent real numbers within the interval $\\pm 10^{-38} \\cdots 10^{38}$, with $7-8$ decimal digits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "Decimal digits contributions for single precision number\n",
      "1.19209289551e-07 3.0517578125e-05 0.03125 \n",
      "\n",
      "Largest and smallest exponent for single precision number\n",
      "3.40282366921e+38 5.87747175411e-39 \n",
      "\n"
     ]
    }
   ],
   "source": [
    "#Decimal digits \n",
    "print \"\\n\"\n",
    "print \"Decimal digits contributions for single precision number\"\n",
    "print 2**-23., 2**-15., 2**-5. , \"\\n\"\n",
    "\n",
    "#Largest and smallest exponent  \n",
    "suma = 0 \n",
    "for i in xrange(0,8):\n",
    "    suma += 2**i\n",
    "print \"Largest and smallest exponent for single precision number\"    \n",
    "print 2**(suma-127.), 2**(-127.),\"\\n\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Double-precision numbers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Double-precision numbers are used when high accuracy is required. These numbers are represented by a **64-bits** (Binary digIT) lenght binary number, where the real number is stored following the next rules:\n",
    "\n",
    "![64-bits](http://upload.wikimedia.org/wikipedia/commons/thumb/a/a9/IEEE_754_Double_Floating_Point_Format.svg/618px-IEEE_754_Double_Floating_Point_Format.svg.png)\n",
    "\n",
    "1. The fist digit (called *s*) indicates the sign of the number (s=0 means a positive number, s=1 a negative one).\n",
    "2. The next 11 bits represent the exponent of the number.\n",
    "3. The last  bits represent the fractional part of the number.\n",
    "\n",
    "The formula for recovering the real number is then given by:\n",
    "\n",
    "$$r = (-1)^s\\times \\left( 1 + \\sum_{i=1}^{52}b_{52-i}2^{-i} \\right)\\times 2^{e-1023}$$\n",
    "\n",
    "where $s$ is the sign, $b_{23-i}$ the fraction bits and $e$ is given by:\n",
    "\n",
    "$$e = \\sum_{i=0}^{10} b_{52+i}2^i$$\n",
    "\n",
    "Double-precision system can represent real numbers within the interval $\\pm 10^{-308} \\cdots 10^{308}$, with $16-17$ decimal digits."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###**ACTIVITY**\n",
    "\n",
    "**1.** Write a python script that calculates the double precision number represented by a 64-bits binary.\n",
    "\n",
    "    \n",
    "**2.** What is the number represented by:\n",
    "\n",
    "0 10000000011 1011100100001111111111111111111111111111111111111111\n",
    "\n",
    "<font color='white'>\n",
    "    **ANSWER:**  27.56640625"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Finite Arithmetic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The most basic arithmetic operations are addition and multiplication. Further operations such as subtraction, division and power are secondary as they can be reached by iteratively use the latter ones."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Addition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As mentioned before, arithmetic operations are not exact in a computer due to the inherent limitations in number representing. Even when adding two already approximate numbers, say a single-precision couple of numbers, the result may not be a representable number, being necessary to apply approximation rules."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.999755859375\n"
     ]
    }
   ],
   "source": [
    "N = 9\n",
    "x = 0\n",
    "for i in xrange(N):\n",
    "    x += np.float16(1.0/N)\n",
    "print x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the sucessive application of rounded-off numbers produces a final result less precise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5/7 0.714286\n",
      "1/3 0.333333\n",
      "1.04762 1.04761904762\n",
      "Error: 5.67663283046e-08\n"
     ]
    }
   ],
   "source": [
    "print \"5/7\", np.float32(5/7.)\n",
    "print \"1/3\", np.float32(1/3.)\n",
    "print np.float32(5/7.+1/3.), 22/21.\n",
    "print \"Error:\", np.float32(5/7.+1/3.)-22/21."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although the **float16** or half-float precision is standard according to the IEEE 754-2008, many devices do not support it well."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multiplication"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For multiplication it is applied the same round-off rules as the addition, however, be aware that multiplicative errors propagate more quickly than additive errors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.99580530418 0.71436\n"
     ]
    }
   ],
   "source": [
    "N = 20\n",
    "x = 1\n",
    "for i in xrange(N):\n",
    "    x *= np.float16(2.0**(1.0/N))\n",
    "print x, np.float16(5/7.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The final result has an error at the third decimal digit, one more than the case of addition."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**ACTIVITY**\n",
    "\n",
    "Find the error associated to the finite representation in the next operations \n",
    "\n",
    "\n",
    "\n",
    "$$\n",
    "x-u, \\frac{x-u}{w}, (x-u)*v, u+v \n",
    "$$\n",
    "\n",
    "considering the values \n",
    "\n",
    "$$\n",
    "x = \\frac{5}{7}, y = \\frac{1}{3}, u = 0.71425\n",
    "$$\n",
    "\n",
    "\n",
    "\n",
    "$$\n",
    "v = 0.98765\\times 10^5, w = 0.111111\\times 10^{-4}\n",
    "$$\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}