gistfile1.txt

{
 "metadata": {
  "name": "Calculating pi"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The Indian mathematician, Ramanujan did most of his works on number series. Several of these series was about calculating the value of $\\pi$.\n",
      "\n",
      "The following series to calculate $\\pi$ is one of his works, which was done back in 1914:\n",
      "$\\frac{1}{\\pi} = \\frac{2\\sqrt{2}}{9801} \\sum_{k=0}^{\\infty} \\frac{(4k)! (1103+26390k)}{(k!)^4 396^{4k}}$\n",
      "\n",
      "Even today, his series are the base of the fastest method to calculate $\\pi$.\n",
      "\n",
      "We can translate his series into code. Each \"step\" in the series will produce approximately 8 new digits of $\\pi$.\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "First we define some constants and helper functions needed. We see at each step the factor is constant $\\frac{2 \\sqrt{2}}{9801}$. Also Python doesn't have a built-in factorial function."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from decimal import *\n",
      "\n",
      "getcontext().prec = 1000 # make decimal precise to 1k points.\n",
      "\n",
      "factorial = lambda n: reduce((lambda a, b: a*b), range(1, n+1)) if n > 0 else 1\n",
      "\n",
      "factor = Decimal(2) * Decimal(2).sqrt() / Decimal(9801)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": "*"
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we're ready to define the function to calculate $\\pi$ based on the series. We do not continue $\\infty$, but set a precision given by how many digits we wish to calculate, and allow the function to run until we have reached the desired precision."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "step_val = lambda k: factor * (factorial(Decimal(4*k)) * Decimal(1103 + 26390*k)) /\n",
      "                     (factorial(Decimal(k))**4 * Decimal(396**(4*k)))\n",
      "\n",
      "def pi():\n",
      "    digits_val = Decimal(1) / 1000\n",
      "    total      = Decimal(0)\n",
      "    k          = 0\n",
      "    while True:\n",
      "        total      += step_val(k)\n",
      "        if abs(step_val) < digits_val: # check if current precision is below wanted digits.\n",
      "            break\n",
      "        k += 1\n",
      "        \n",
      "    return Decimal(1) / total"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": "*"
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can test the method by looking at historical values for $\\pi$. For example Aitken recited $\\pi$ to the first 1000 digits as follows:\n",
      "\n",
      "<pre>\n",
      "3.14159265358979323846264338327950288419716939937510\n",
      "  58209749445923078164062862089986280348253421170679\n",
      "  82148086513282306647093844609550582231725359408128\n",
      "  48111745028410270193852110555964462294895493038196\n",
      "  44288109756659334461284756482337867831652712019091\n",
      "  45648566923460348610454326648213393607260249141273\n",
      "  72458700660631558817488152092096282925409171536436\n",
      "  78925903600113305305488204665213841469519415116094\n",
      "  33057270365759591953092186117381932611793105118548\n",
      "  07446237996274956735188575272489122793818301194912\n",
      "  98336733624406566430860213949463952247371907021798\n",
      "  60943702770539217176293176752384674818467669405132\n",
      "  00056812714526356082778577134275778960917363717872\n",
      "  14684409012249534301465495853710507922796892589235\n",
      "  42019956112129021960864034418159813629774771309960\n",
      "  51870721134999999837297804995105973173281609631859\n",
      "  50244594553469083026425223082533446850352619311881\n",
      "  71010003137838752886587533208381420617177669147303\n",
      "  59825349042875546873115956286388235378759375195778\n",
      "  18577805321712268066130019278766111959092164201989\n",
      "</pre>"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pi_aitken = Decimal(\"\"\"\n",
      "3.14159265358979323846264338327950288419716939937510\n",
      "  58209749445923078164062862089986280348253421170679\n",
      "  82148086513282306647093844609550582231725359408128\n",
      "  48111745028410270193852110555964462294895493038196\n",
      "  44288109756659334461284756482337867831652712019091\n",
      "  45648566923460348610454326648213393607260249141273\n",
      "  72458700660631558817488152092096282925409171536436\n",
      "  78925903600113305305488204665213841469519415116094\n",
      "  33057270365759591953092186117381932611793105118548\n",
      "  07446237996274956735188575272489122793818301194912\n",
      "  98336733624406566430860213949463952247371907021798\n",
      "  60943702770539217176293176752384674818467669405132\n",
      "  00056812714526356082778577134275778960917363717872\n",
      "  14684409012249534301465495853710507922796892589235\n",
      "  42019956112129021960864034418159813629774771309960\n",
      "  51870721134999999837297804995105973173281609631859\n",
      "  50244594553469083026425223082533446850352619311881\n",
      "  71010003137838752886587533208381420617177669147303\n",
      "  59825349042875546873115956286388235378759375195778\n",
      "  18577805321712268066130019278766111959092164201989\"\"\".replace(\"\\n\", \"\").replace(\" \",\"\"))\n",
      "\n",
      "pi_r = str(pi())\n",
      "\n",
      "for x in [ pi_r[x:x+50] for x in range(0, len(pi_r), 50) ]:\n",
      "    print x\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "3.141592653589793877998905826306013094216645029322\n",
        "84887917396379150578440064851167469327839384943948\n",
        "18144977362409840921649870099469289482700906055700\n",
        "08232448569070565225117261929671542266072657406722\n",
        "90709808328625814224710211298682485180629048757430\n",
        "98898800845798127138825309404173040482622867066073\n",
        "75938877516032550654677185287604746228636591216480\n",
        "12735103486547664740694156241119551876454392787753\n",
        "21505579554979836645433665279782657168281439378387\n",
        "28277223685045670035214900432839230508193216393433\n",
        "51314936535921545202929797110682358958754390440264\n",
        "03274443324585783524364308738628065910884486435913\n",
        "33438579491981872958490965113844396301375105180910\n",
        "30421401751866665015216405974399352597351913912127\n",
        "84344167910354391500166890711450850452292129140834\n",
        "63490624556734326918499821494683541355534470994519\n",
        "19713559759727285669336702201068084247758003830617\n",
        "21463488377300493137375051032516374542525216125277\n",
        "94273688664200248278275264156995095320711353610503\n",
        "40798331126329185941155506107015675450606088823007\n",
        "1\n"
       ]
      }
     ],
     "prompt_number": 25
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "False\n"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}