ipython_nbs/efficiency/greatest-common-divisor.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sebatian Raschka \n",
      "last updated: 2016-06-02 \n",
      "\n",
      "CPython 3.5.1\n",
      "IPython 4.2.0\n"
     ]
    }
   ],
   "source": [
    "%load_ext watermark\n",
    "%watermark -a 'Sebatian Raschka' -u -d -v"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Greatest Common Divisor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, the *greatest common divisor* (GCD) is the largest natural number $d$ that divides $a$ and $b$ in a fraction $\\frac{a}{b}$ without a remainder. \n",
    "\n",
    "For example, the GCD of the fraction $\\frac{6}{9}$ is 3: $$\\frac{6/3}{9/3} = \\frac{2}{3}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let us start with the \"intuitive,\" yet naive, brute-force implementation. Here, we simply iterate through all integers from 1 to max(a, b) to find the largest common divisor of the fraction $\\frac{a}{b}$ or equivalently $\\frac{a}{b}$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "In: 1/1, Out: 1\n",
      "In: 1/2, Out: 1\n",
      "In: 3/9, Out: 3\n",
      "In: 12/24, Out: 12\n",
      "In: 12/26, Out: 2\n",
      "In: 26/12, Out: 2\n",
      "In: 13/17, Out: 1\n"
     ]
    }
   ],
   "source": [
    "def naive_gcd(a, b):\n",
    "    gcd = 0\n",
    "    if a < b:\n",
    "        n = a\n",
    "    else:\n",
    "        n = a\n",
    "    for d in range(1, n + 1):\n",
    "        if not a % d and not b % d:\n",
    "            gcd = d\n",
    "    return gcd\n",
    "\n",
    "print('In: 1/1,', 'Out:', naive_gcd(1, 1))\n",
    "print('In: 1/2,', 'Out:', naive_gcd(1, 2))\n",
    "print('In: 3/9,', 'Out:', naive_gcd(3, 9))\n",
    "print('In: 12/24,', 'Out:', naive_gcd(12, 24))\n",
    "print('In: 12/26,', 'Out:', naive_gcd(12, 26))\n",
    "print('In: 26/12,', 'Out:', naive_gcd(26, 12))\n",
    "print('In: 13/17,', 'Out:', naive_gcd(13, 17))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Euclidean Algorithm\n",
    "\n",
    "The Greek mathematician Euclid described this algorithm approx. 300 BC in \n",
    "It is still being used in the field of number theory and conseqently cryptography to reduce fractions to their simplest form. For additional, interesting details and the proof by Gabriel Lamé in 1844, please take a look at the excellent [Wikipedia](https://en.wikipedia.org/wiki/Euclidean_algorithm) page.\n",
    "\n",
    "The idea behind the more efficient version, the Euclidean division (in contrast to the substraction-based) approach is the following:\n",
    "\n",
    "Given that we want to compute `gcd(a, b)`, the greatest common divisor of the fraction{a}{b}, we first compute the remainder of the division $\\frac{a}{b}$; we call this remainder a'. Then, we compute `gcd(a', b)`. We repeat this procedure in recursive manner until b=0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "In: 1/1, Out: 1\n",
      "In: 1/2, Out: 1\n",
      "In: 3/9, Out: 3\n",
      "In: 12/24, Out: 12\n",
      "In: 12/26, Out: 2\n",
      "In: 26/12, Out: 2\n",
      "In: 13/17, Out: 1\n"
     ]
    }
   ],
   "source": [
    "def eucl_gcd_recurse(a, b):\n",
    "    if not b:\n",
    "        return a\n",
    "    else:\n",
    "        return eucl_gcd_recurse(b, a % b)\n",
    "    \n",
    "print('In: 1/1,', 'Out:', naive_gcd(1, 1))\n",
    "print('In: 1/2,', 'Out:', naive_gcd(1, 2))\n",
    "print('In: 3/9,', 'Out:', naive_gcd(3, 9))\n",
    "print('In: 12/24,', 'Out:', naive_gcd(12, 24))\n",
    "print('In: 12/26,', 'Out:', naive_gcd(12, 26))\n",
    "print('In: 26/12,', 'Out:', naive_gcd(26, 12))\n",
    "print('In: 13/17,', 'Out:', naive_gcd(13, 17))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Euclidean GCD algorithm will reduce either of the number at least by half at each step (see the [Wikipedia](https://en.wikipedia.org/wiki/Euclidean_algorithm) page for details). Thus, the time complexity of this algorithm is \n",
    "\n",
    "$O(log_2 b) + O(log_2 a) = O(log_2 n),$\n",
    "\n",
    "where $n = max(a, b)$. In contrast, our previous, naive implementation has an upper bound of $O(n)$.\n",
    "\n",
    "Since Python is \"notoriously bad\" at recursion, let us implement a dynamic version of this algorithm. (One problem of recursion via Python is the limited stack size, and the other one is that tail recursion optimization not implemented.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "In: 1/1, Out: 1\n",
      "In: 1/2, Out: 1\n",
      "In: 3/9, Out: 3\n",
      "In: 12/24, Out: 12\n",
      "In: 12/26, Out: 2\n",
      "In: 26/12, Out: 2\n",
      "In: 13/17, Out: 1\n"
     ]
    }
   ],
   "source": [
    "def eucl_gcd_dynamic(a, b):\n",
    "    while b:\n",
    "       tmp = b \n",
    "       b = a % b \n",
    "       a = tmp \n",
    "    return a\n",
    "\n",
    "print('In: 1/1,', 'Out:', naive_gcd(1, 1))\n",
    "print('In: 1/2,', 'Out:', naive_gcd(1, 2))\n",
    "print('In: 3/9,', 'Out:', naive_gcd(3, 9))\n",
    "print('In: 12/24,', 'Out:', naive_gcd(12, 24))\n",
    "print('In: 12/26,', 'Out:', naive_gcd(12, 26))\n",
    "print('In: 26/12,', 'Out:', naive_gcd(26, 12))\n",
    "print('In: 13/17,', 'Out:', naive_gcd(13, 17))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given an arbitrary fraction $\\frac{a}{b}$, let us use the `%timeit` module for a quick comparison:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "a = 12313432\n",
    "b = 34234232342"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5 loops, best of 3: 1.78 s per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit -r 3 -n 5 naive_gcd(a, b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5 loops, best of 3: 3.23 µs per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit -r 3 -n 5 eucl_gcd_recurse(a, b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5 loops, best of 3: 2.21 µs per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit -r 3 -n 5 eucl_gcd_dynamic(a, b)"
   ]
  }
 ],
 "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}