{ "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": {} } ] }