numerics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How to Break Math\n",
    "## Numerical Computing with Floating-Point Numbers\n",
    "\n",
    "Real numbers may have an infinite number of decimal digits: for example, $\\frac{1}{3} = 0.333\\dots$, and $\\pi = 3.14159\\dots$. However, computers do not have an infinite amount of storage, and for efficiency, most computer programs represent real numbers using floating-point values.\n",
    "\n",
    "Floating point is similar to scientific notation. For example, the speed of light can be written simply as $299792458\\,ms^{-1}$ or in scientific notation, $2.99792458 \\times 10^8\\,ms^{-1}$. A number in scientific notation has three parts: the *significand* or *mantissa*, the *base*, and the *exponent*:\n",
    "\\begin{equation}\n",
    "299792458 = \\underbrace{2.99792458}_{\\text{significand}} \\times {\\underbrace{10}_{\\text{base}}}^{\\hspace{-0.6em}\\overbrace{8}^{\\text{exponent}}}\n",
    "\\end{equation}\n",
    "\n",
    "By convention, the significand has one non-zero digit before the point (unless the number is exactly zero).\n",
    "\n",
    "In binary floating point, the base is always 2, and the signficand is made up entirely of the digits 0 and 1. For example, the number 42 can be represented as a binary integer as 101010, or in binary floating point:\n",
    "\\begin{equation}\n",
    "42 = 1.01010 \\times 2^5\n",
    "\\end{equation}\n",
    "\n",
    "Computers use the [IEEE-754 standard floating point format](https://en.wikipedia.org/wiki/IEEE_754), which defines standard \"single\" and \"double\" precision floating-point numbers. The bits which represent an IEEE-754 value are divided into three fields: the sign bit, the exponent, and the significand.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "caption": "Image source: Fresheneesz (Wikimedia Commons) CC-ShareAlike 3.0"
   },
   "source": [
    "![IEEE-754 Single-Precision Floating Point format](images/IEEE754_single.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A sign bit of `0` represents a positive number, while a sign bit of `1` represents a negative number. \n",
    "\n",
    "The leading '1' in the significand is not stored, so the value `0b0100...` represents the signficand '1.0100...'. \n",
    "\n",
    "The exponent is *biased* - it is read as an integer value, from which the value of the bias is subtracted. For example, in single precision the bias is 127, so the exponent `0b10000100` (decimal 132) is read as (132-127) = 5, and the exponent `0b01111100` (decimal 124) is read as (124-127) = -3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function print_bitstring(a::Float32)\n",
    "    s = bitstring(a)\n",
    "    printstyled(s[1], color=:blue)\n",
    "    print(\" \")\n",
    "    printstyled(s[2:9], color=:green)\n",
    "    print(\" \")\n",
    "    printstyled(s[10:32], color=:red)\n",
    "    println()\n",
    "end\n",
    "\n",
    "function print_bitstring(a::Float64)\n",
    "    # TODO complete - print the bits of the double-precision floating-point number a, separated into fields\n",
    "end\n",
    "\n",
    "function exponent_ieee754(a::Float32)\n",
    "    s = bitstring(a)\n",
    "    e = parse(Int32,s[2:9],base=2)\n",
    "    return e - (2^7-1)\n",
    "end\n",
    "\n",
    "function exponent_ieee754(a::Float64)\n",
    "    # TODO complete - extract exponent field from the double-precision floating point number a\n",
    "end\n",
    "    \n",
    "x = Float32(42)\n",
    "print(x,\" = \")\n",
    "print_bitstring(x)\n",
    "println(\"exponent = \", exponent_ieee754(x))\n",
    "\n",
    "y = Float32(0.15625)\n",
    "print(y,\" = \")\n",
    "print_bitstring(y)\n",
    "println(\"exponent = \", exponent_ieee754(y))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Rounding Error\n",
    "\n",
    "Rounding error occurs whenever an input value or the result of a calculation cannot be represented exactly in the chosen floating-point format, for example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print_bitstring(0.1)\n",
    "print_bitstring(Float64(pi))\n",
    "print_bitstring(1.0/3.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cancellation\n",
    "\n",
    "Cancellation occurs when two close-together values are subtracted, leading to a large relative error in the result. For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "epsilon = 1e-8\n",
    "x = (1+epsilon)\n",
    "print(\"x \")\n",
    "print_bitstring(x)\n",
    "y = (1.0)\n",
    "print(\"y \")\n",
    "print_bitstring(y)\n",
    "difference = epsilon\n",
    "z = x-y\n",
    "print(\"z \")\n",
    "print_bitstring(z)\n",
    "println(\"Correct result is \", difference, \" but got \", z, \" (relative error \", (difference-z)/difference, \")\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rump's Example\n",
    "\n",
    "Siegfried Rump created a [nasty example](http://www.ti3.tu-harburg.de/paper/rump/Ru88a.pdf) of floating point error. He asks us to compute the value of the function:\n",
    "\n",
    "$333.75 b^6 + a^2 (11 a^2 b^2 – b^6 – 121 b^4 – 2) + 5.5 b^8 + \\frac{a}{(2b)}$\n",
    "\n",
    "for $a = 77617.0$ and $b = 33096.0$ . We can define a Julia function to compute the value using floating-point numbers of type `T`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function rump(T)\n",
    "  a = T(77617.0)\n",
    "  b = T(33096.0)\n",
    "  a2 = a*a\n",
    "  b2 = b*b\n",
    "  b4 = b2*b2\n",
    "  b6 = b4*b2\n",
    "  b8 = b4*b4\n",
    "  return T(333.75)*b6 + a2 * (11*a2*b2 - b6 - 121*b4 - 2) + T(5.5)*b8 + (a/(2*b))\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compute that function in single-precision floating point:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rump(Float32)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To check the correctness, let's compute it again in double-precision:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rump(Float64)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Oh dear! Our results aren't even remotely the same - they're out by more than eight orders of magnitude. \n",
    "\n",
    "Let's try increasing the precision still further. Julia supports arbitrary-precision arithmetic using [GNU MPFR](https://www.mpfr.org/) with a special type `BigFloat`. We can use this to create a function that will compute Rump's example using a chosen number of bits of precision in the significand:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function rump(numbits::Int)\n",
    "  Base.MPFR.setprecision(numbits)\n",
    "  return rump(BigFloat)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's replicate our previous results using BigFloat. First single-precision, which in [IEEE-754](https://en.wikipedia.org/wiki/IEEE_754) has 24 bits in the significand:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rump(24)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, double precision, which has 53 bits in the significand:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rump(53)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's try computing in quadruple precision, which in IEEE-754 has 113 bits in the significand:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rump(113)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This looks even worse - now even the sign of the result has changed!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i = 113:130\n",
    "    println(\"Precision \", i, \" result is \", rump(i))\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rump gives the exact answer as $\\frac{-54767}{66192} = -0.8273960599\\dots$\n",
    "\n",
    "Why do we need such high precision to correctly compute Rump's example? Use the `bitstring_ieee754(Float64)` function you created earlier to print the representation of some of the subexpressions that are computed in the function, to see if you can understand why the error is so great."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function rump_print(T)\n",
    "  a = T(77617.0)\n",
    "  b = T(33096.0)\n",
    "  a2 = a*a\n",
    "  b2 = b*b\n",
    "  b4 = b2*b2\n",
    "  b6 = b4*b2\n",
    "  b8 = b4*b4\n",
    "  # TODO print some of the subexpressions computed below\n",
    "  intermediate_value1 = # ?\n",
    "  intermediate_value2 = # ?\n",
    "  print_bitstring(intermediate_value1)\n",
    "  print_bitstring(intermediate_value2)\n",
    "  return T(333.75)*b6 + a2 * (11*a2*b2 - b6 - 121*b4 - 2) + T(5.5)*b8 + (a/(2*b))\n",
    "end\n",
    "\n",
    "result = rump_print(Float64)\n",
    "print_bitstring(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Patriot vs. Scud\n",
    "\n",
    "The [Patriot missile battery](http://www-users.math.umn.edu/~arnold//disasters/patriot.html) counted time using a clock which ticked ten times each second. To compute the expected location of an incoming target, the number of ticks was divided by 10 to produce a fixed-point value representing the number of seconds since the system started up. We can simulate this using a single-precision floating point value.\n",
    "\n",
    "Modify the code below to compute the actual time in seconds given the current value of the clock. Why is `current_time_seconds` different? Hint: try printing out the bitstring representation of `tick_length`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Printf\n",
    "const TICKS_PER_SECOND = 10\n",
    "const SCUD_MISSILE_SPEED_METRES_PER_SECOND = 1676\n",
    "\n",
    "function on_the_hour(ticks)\n",
    "  return (ticks % (TICKS_PER_SECOND * 60 * 60) == 0)\n",
    "end\n",
    "\n",
    "Base.MPFR.setprecision(20)\n",
    "tick_internal = BigFloat(1 / TICKS_PER_SECOND)\n",
    "tick_length = Float32(tick_internal)\n",
    "\n",
    "current_ticks = 0\n",
    "current_time_seconds = 0.0\n",
    "one_hundred_hours = 100 * 60 * 60 * TICKS_PER_SECOND\n",
    "for i=0:one_hundred_hours\n",
    "    current_ticks += 1\n",
    "    current_time_seconds += tick_length\n",
    "    if (on_the_hour(current_ticks))\n",
    "        actual_time_seconds = # TODO compute correct time in seconds\n",
    "        h = actual_time_seconds / 3600\n",
    "        error_time_seconds = actual_time_seconds - current_time_seconds\n",
    "        @printf(\"hour = %3d; error (time) = %.2e s; error (target position) %.2e m\\n\", h, \n",
    "            error_time_seconds, error_time_seconds*SCUD_MISSILE_SPEED_METRES_PER_SECOND)\n",
    "    end\n",
    "end"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.4.0",
   "language": "julia",
   "name": "julia-1.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.4.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}