{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Lab session: Feb 10 - 14"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Some basic number theory"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Sieve of Eratosthenes"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The problem is the following. Input a number $n$. Then find all the prime numbers between $1$ and $n$ ($n$ inclusive) using the Sieve of Eratosthenes. For those, who have forgotten or don't know what Sieve of Eratosthenes is, here is a short description :\n",
      "\n",
      "1. Make a list of numbers $2, 3, 4, ... n$.\n",
      "2. Set p = 2 and mark out every 2nd number after that; i.e., mark 4, 6, ...\n",
      "3. Set p to be the next unmarked number, here 3.\n",
      "4. Mark out every 3rd number after 3; i.e. mark 6, 9, ...\n",
      "5. Set p to be the next unmarked no and so on."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Remark : The following solution uses a bit more memomry than what is absolutely needed. I'll give a modified solution after this."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# We shall replace the number by tuples (n, mark) where mark is true or false.\n",
      "# Since we have to mark numbers, we write a function\n",
      "def mark(list_of_marked_pairs, k) :\n",
      "    # We know that the list will start from 2 and go upto n. So k should be in\n",
      "    # the k-2'th place.\n",
      "    if k <= 2 or k >= n+2 :\n",
      "        print \"Number out of range : Cannot mark.\"\n",
      "    else :\n",
      "        list_of_marked_pairs[k-2] = (k, True)\n",
      "\n",
      "def sieve(n) :\n",
      "    primes = [(k, False) for k in range(2, n+1)]\n",
      "    # (n, True) means n is marked\n",
      "    # Uncomment to see step by step evaluation :\n",
      "    # print primes\n",
      "    \n",
      "    p = 2 \n",
      "    # p is the present prime\n",
      "    while p <= n :\n",
      "        for x in range(2*p, n+1, p) :\n",
      "            mark(primes, x)\n",
      "        # Uncomment to see step by step evaluation : \n",
      "        # print primes\n",
      "        \n",
      "        # Now find the next unmarked number after p and mark that as p\n",
      "        k = p+1\n",
      "        while k <= n and primes[k-2][1] == True :\n",
      "            k += 1\n",
      "        p = k\n",
      "    # Now the primes are those numbers which have false associated with them\n",
      "    list_of_primes = []\n",
      "    for n, m in primes :\n",
      "        if m == False :\n",
      "            list_of_primes.append(n)\n",
      "    return list_of_primes\n",
      "\n",
      "n = input(\"Enter a number greather than 2 : \")\n",
      "print sieve(n)\n",
      "    \n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "name": "stdout",
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Enter a number greather than 2 : 100\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#Another solution\n",
      "def sieve_eff(n) :\n",
      "    prime_marks = [False for i in range(2, n+1)]\n",
      "    # print prime_marks\n",
      "    # Note prime_marks[k] corresponds to the numbe k + 2.\n",
      "    #p = present prime\n",
      "    p = 2\n",
      "    while p <= n :\n",
      "        for i in range(2*p, n+1, p) :\n",
      "            prime_marks[i-2] = True\n",
      "        # The next p :\n",
      "        k = p+1\n",
      "        while k <= n and prime_marks[k-2] == True :\n",
      "            k += 1\n",
      "        p = k\n",
      "    # print prime_marks\n",
      "    \n",
      "    # Now collect the primes :\n",
      "    primes = []\n",
      "    for n, val in enumerate(prime_marks) :\n",
      "        if val == False :\n",
      "            primes.append(n+2)\n",
      "    return primes\n",
      "        \n",
      "n = input(\"Input a number greater than 2 : \")\n",
      "print sieve_eff(n)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "name": "stdout",
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Input a number greater than 2 : 100\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Pythagorean triples"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Given an integer $N$ find all triples $(a, b, c)$ where $a, b, c$ are integers lying between $1$ and $N$ such that $a^2 + b^2 = c^2$. In this first part of the problem we consider $(3, 4, 5)$ and $(4, 3, 5)$ to be different triples."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def pythagorean_triples(N) :\n",
      "    pt = [] # pt = pythagorean triples\n",
      "    for a in range(1, N+1) :\n",
      "        for b in range(1, N+1) :\n",
      "            for c in range(1, N+1) :\n",
      "                if a*a + b*b - c*c == 0 :\n",
      "                    pt.append((a, b, c))\n",
      "    return pt\n",
      "\n",
      "N = input(\"Find Pythagorean triples less than or equal to N = ? \")\n",
      "for i, j, k in pythagorean_triples(N):\n",
      "    print \"%4d^2 + %4d^2 = %4d^2\" %(i, j, k)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "name": "stdout",
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Find Pythagorean triples less than or equal to N = ? 100\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "   3^2 +    4^2 =    5^2\n",
        "   4^2 +    3^2 =    5^2\n",
        "   5^2 +   12^2 =   13^2\n",
        "   6^2 +    8^2 =   10^2\n",
        "   7^2 +   24^2 =   25^2\n",
        "   8^2 +    6^2 =   10^2\n",
        "   8^2 +   15^2 =   17^2\n",
        "   9^2 +   12^2 =   15^2\n",
        "   9^2 +   40^2 =   41^2\n",
        "  10^2 +   24^2 =   26^2\n",
        "  11^2 +   60^2 =   61^2\n",
        "  12^2 +    5^2 =   13^2\n",
        "  12^2 +    9^2 =   15^2\n",
        "  12^2 +   16^2 =   20^2\n",
        "  12^2 +   35^2 =   37^2\n",
        "  13^2 +   84^2 =   85^2\n",
        "  14^2 +   48^2 =   50^2\n",
        "  15^2 +    8^2 =   17^2\n",
        "  15^2 +   20^2 =   25^2\n",
        "  15^2 +   36^2 =   39^2\n",
        "  16^2 +   12^2 =   20^2\n",
        "  16^2 +   30^2 =   34^2\n",
        "  16^2 +   63^2 =   65^2\n",
        "  18^2 +   24^2 =   30^2\n",
        "  18^2 +   80^2 =   82^2\n",
        "  20^2 +   15^2 =   25^2\n",
        "  20^2 +   21^2 =   29^2\n",
        "  20^2 +   48^2 =   52^2\n",
        "  21^2 +   20^2 =   29^2\n",
        "  21^2 +   28^2 =   35^2\n",
        "  21^2 +   72^2 =   75^2\n",
        "  24^2 +    7^2 =   25^2\n",
        "  24^2 +   10^2 =   26^2\n",
        "  24^2 +   18^2 =   30^2\n",
        "  24^2 +   32^2 =   40^2\n",
        "  24^2 +   45^2 =   51^2\n",
        "  24^2 +   70^2 =   74^2\n",
        "  25^2 +   60^2 =   65^2\n",
        "  27^2 +   36^2 =   45^2\n",
        "  28^2 +   21^2 =   35^2\n",
        "  28^2 +   45^2 =   53^2\n",
        "  28^2 +   96^2 =  100^2\n",
        "  30^2 +   16^2 =   34^2\n",
        "  30^2 +   40^2 =   50^2\n",
        "  30^2 +   72^2 =   78^2\n",
        "  32^2 +   24^2 =   40^2\n",
        "  32^2 +   60^2 =   68^2\n",
        "  33^2 +   44^2 =   55^2\n",
        "  33^2 +   56^2 =   65^2\n",
        "  35^2 +   12^2 =   37^2\n",
        "  35^2 +   84^2 =   91^2\n",
        "  36^2 +   15^2 =   39^2\n",
        "  36^2 +   27^2 =   45^2\n",
        "  36^2 +   48^2 =   60^2\n",
        "  36^2 +   77^2 =   85^2\n",
        "  39^2 +   52^2 =   65^2\n",
        "  39^2 +   80^2 =   89^2\n",
        "  40^2 +    9^2 =   41^2\n",
        "  40^2 +   30^2 =   50^2\n",
        "  40^2 +   42^2 =   58^2\n",
        "  40^2 +   75^2 =   85^2\n",
        "  42^2 +   40^2 =   58^2\n",
        "  42^2 +   56^2 =   70^2\n",
        "  44^2 +   33^2 =   55^2\n",
        "  45^2 +   24^2 =   51^2\n",
        "  45^2 +   28^2 =   53^2\n",
        "  45^2 +   60^2 =   75^2\n",
        "  48^2 +   14^2 =   50^2\n",
        "  48^2 +   20^2 =   52^2\n",
        "  48^2 +   36^2 =   60^2\n",
        "  48^2 +   55^2 =   73^2\n",
        "  48^2 +   64^2 =   80^2\n",
        "  51^2 +   68^2 =   85^2\n",
        "  52^2 +   39^2 =   65^2\n",
        "  54^2 +   72^2 =   90^2\n",
        "  55^2 +   48^2 =   73^2\n",
        "  56^2 +   33^2 =   65^2\n",
        "  56^2 +   42^2 =   70^2\n",
        "  57^2 +   76^2 =   95^2\n",
        "  60^2 +   11^2 =   61^2\n",
        "  60^2 +   25^2 =   65^2\n",
        "  60^2 +   32^2 =   68^2\n",
        "  60^2 +   45^2 =   75^2\n",
        "  60^2 +   63^2 =   87^2\n",
        "  60^2 +   80^2 =  100^2\n",
        "  63^2 +   16^2 =   65^2\n",
        "  63^2 +   60^2 =   87^2\n",
        "  64^2 +   48^2 =   80^2\n",
        "  65^2 +   72^2 =   97^2\n",
        "  68^2 +   51^2 =   85^2\n",
        "  70^2 +   24^2 =   74^2\n",
        "  72^2 +   21^2 =   75^2\n",
        "  72^2 +   30^2 =   78^2\n",
        "  72^2 +   54^2 =   90^2\n",
        "  72^2 +   65^2 =   97^2\n",
        "  75^2 +   40^2 =   85^2\n",
        "  76^2 +   57^2 =   95^2\n",
        "  77^2 +   36^2 =   85^2\n",
        "  80^2 +   18^2 =   82^2\n",
        "  80^2 +   39^2 =   89^2\n",
        "  80^2 +   60^2 =  100^2\n",
        "  84^2 +   13^2 =   85^2\n",
        "  84^2 +   35^2 =   91^2\n",
        "  96^2 +   28^2 =  100^2\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now try to print the list where permutations are considered to be the same triple."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Note that c is the largest. We can get rid of permutations just by imposing a <= b\n",
      "\n",
      "def pythagorean_triples_2(N) :\n",
      "    pt = []\n",
      "    for a in range(1, N+1) :\n",
      "        for b in range(a, N+1) :\n",
      "            for c in range(b, N+1) :\n",
      "                # c cannot be less than b!\n",
      "                if a*a + b*b - c*c == 0 :\n",
      "                    pt.append((a, b, c))\n",
      "    return pt\n",
      "\n",
      "N = input(\"N = \")\n",
      "for a, b, c in pythagorean_triples_2(N) :\n",
      "    print \"%4d^2 + %4d^2 = %4d^2\" % (a, b, c)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "name": "stdout",
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "N = 100\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "   3^2 +    4^2 =    5^2\n",
        "   5^2 +   12^2 =   13^2\n",
        "   6^2 +    8^2 =   10^2\n",
        "   7^2 +   24^2 =   25^2\n",
        "   8^2 +   15^2 =   17^2\n",
        "   9^2 +   12^2 =   15^2\n",
        "   9^2 +   40^2 =   41^2\n",
        "  10^2 +   24^2 =   26^2\n",
        "  11^2 +   60^2 =   61^2\n",
        "  12^2 +   16^2 =   20^2\n",
        "  12^2 +   35^2 =   37^2\n",
        "  13^2 +   84^2 =   85^2\n",
        "  14^2 +   48^2 =   50^2\n",
        "  15^2 +   20^2 =   25^2\n",
        "  15^2 +   36^2 =   39^2\n",
        "  16^2 +   30^2 =   34^2\n",
        "  16^2 +   63^2 =   65^2\n",
        "  18^2 +   24^2 =   30^2\n",
        "  18^2 +   80^2 =   82^2\n",
        "  20^2 +   21^2 =   29^2\n",
        "  20^2 +   48^2 =   52^2\n",
        "  21^2 +   28^2 =   35^2\n",
        "  21^2 +   72^2 =   75^2\n",
        "  24^2 +   32^2 =   40^2\n",
        "  24^2 +   45^2 =   51^2\n",
        "  24^2 +   70^2 =   74^2\n",
        "  25^2 +   60^2 =   65^2\n",
        "  27^2 +   36^2 =   45^2\n",
        "  28^2 +   45^2 =   53^2\n",
        "  28^2 +   96^2 =  100^2\n",
        "  30^2 +   40^2 =   50^2\n",
        "  30^2 +   72^2 =   78^2\n",
        "  32^2 +   60^2 =   68^2\n",
        "  33^2 +   44^2 =   55^2\n",
        "  33^2 +   56^2 =   65^2\n",
        "  35^2 +   84^2 =   91^2\n",
        "  36^2 +   48^2 =   60^2\n",
        "  36^2 +   77^2 =   85^2\n",
        "  39^2 +   52^2 =   65^2\n",
        "  39^2 +   80^2 =   89^2\n",
        "  40^2 +   42^2 =   58^2\n",
        "  40^2 +   75^2 =   85^2\n",
        "  42^2 +   56^2 =   70^2\n",
        "  45^2 +   60^2 =   75^2\n",
        "  48^2 +   55^2 =   73^2\n",
        "  48^2 +   64^2 =   80^2\n",
        "  51^2 +   68^2 =   85^2\n",
        "  54^2 +   72^2 =   90^2\n",
        "  57^2 +   76^2 =   95^2\n",
        "  60^2 +   63^2 =   87^2\n",
        "  60^2 +   80^2 =  100^2\n",
        "  65^2 +   72^2 =   97^2\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now notice that if $(a, b, c)$ is a pythagorean triple, then so is $(na, nb, nc)$ for any $n$. Can you write code so that if some triple is obtained from another one by multiplying by some integer, then it is not printed in the list?"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Hint : Suppose $g$ is the g.c.d. of $a$ and $b$. Will $g$ divide $c$? If so, will $(a/g, b/g, c/g)$ be a Pythagorean triple?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Answers to both the questions above are yes.\n",
      "# This gives the hint that we have to only look at those triples where gcd(a,b) = 1\n",
      "\n",
      "def gcd(m, n) :\n",
      "    q = m / n\n",
      "    r = m % n\n",
      "    while r > 0 :\n",
      "        m = n\n",
      "        n = r\n",
      "        q = m / n\n",
      "        r = m % n\n",
      "    return n\n",
      "\n",
      "def gcd2(m, n) :\n",
      "    while n >0 :\n",
      "        if m < n :\n",
      "            k = m\n",
      "            m = n\n",
      "            n = k\n",
      "        else :\n",
      "            m = m - n\n",
      "    return m\n",
      "        \n",
      "def pythagorean_triples_3(N) :\n",
      "    pt = []\n",
      "    for a in range(1, N+1) :\n",
      "        for b in range(a+1, N+1) :\n",
      "            if gcd2(a, b) == 1 :\n",
      "                for c in range(b+1, N+1) :\n",
      "                    if a*a + b*b - c*c == 0 :\n",
      "                        pt.append((a, b, c))\n",
      "    return pt\n",
      "\n",
      "N = input(\"N = \")\n",
      "for a, b, c in pythagorean_triples_3(N) :\n",
      "    print \"%4d^2 + %4d^2 = %4d^2\" % (a, b, c)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "name": "stdout",
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "N = 100\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "   3^2 +    4^2 =    5^2\n",
        "   5^2 +   12^2 =   13^2\n",
        "   7^2 +   24^2 =   25^2\n",
        "   8^2 +   15^2 =   17^2\n",
        "   9^2 +   40^2 =   41^2\n",
        "  11^2 +   60^2 =   61^2\n",
        "  12^2 +   35^2 =   37^2\n",
        "  13^2 +   84^2 =   85^2\n",
        "  16^2 +   63^2 =   65^2\n",
        "  20^2 +   21^2 =   29^2\n",
        "  28^2 +   45^2 =   53^2\n",
        "  33^2 +   56^2 =   65^2\n",
        "  36^2 +   77^2 =   85^2\n",
        "  39^2 +   80^2 =   89^2\n",
        "  48^2 +   55^2 =   73^2\n",
        "  65^2 +   72^2 =   97^2\n"
       ]
      }
     ],
     "prompt_number": 5
    }
   ],
   "metadata": {}
  }
 ]
}