Calculating Pi to arbitrary number of digits

Background

Introduction

I always wanted to calculate π to arbitrary number of decimal places (100, 1000 etc.) just for the fun. Recently the google search led me to Chris Hills' excellent article on the very same subject. It took me some time to understand his code and what he is saying about it through the faq. I am using some ( well, most ) of his code. But there were pieces of code which I had independently thought about ( believe me ). I chose to write the article for my own clarity and thought that someone else could be interested in it as well.

Theory

Calculating π to any arbitrarily large numbers involves two major issues.
  • Every operating system has a finite width to represent integer. On 32 bit operating systems, the integer is often represented by 4 bytes i.e. 32 bits. So large integers cannot be represented without the loss of accuracy. If I want to represent '97616742900186289193' accurately, I cannot use the built-in integer data type because signed integers on 32 bit OS can only represent integers upto 2^32/2 - 1. So the first piece of puzzle is to implement my own integer data type which can hold arbitrarily large strings of digits. Having my own integer data type means I have to implement my own arithmetic operations such as +, - and /. It is very refreshing to implement these operators !
  • The second and more interesting but not so obvious piece is selection of number system ( base ). I started with usual decimal number system i.e. base 10 but Chris' article explained why choosing larger number system is key to the implementation. Large numbers can be represented compactly if the base is big. So if the base is 100, all the number from 0 to 99 would have unique representation and 11 in this system would be 101. The project code configures the base depending upon the number of digits desired beyond decimal of π. Because the number system is not decimal (10), I had to think of the arithmetic operators hard §
  • .

Mathematical equations

I used two formulae to compute π. One is Machin's formula and other is Euler-Gregory-Madhava series. Machin's formula converges fast while as Euler-Gregory-Madhava series is extremely slow to converge ( I will explain why in a minute ).

Machin's formula is given by,

(1)
π
4
   
=      4 * arctan 1
5
    -   arctan 1
239
 
4 * Σ (-1)k
(2k+1)*(52k+1)
   -  Σ (-1)k
(2k+1)*(2392k+1)







The Euler-Gregory-Madhava series is given by,

(2)
π
4
   
=      Σ (-1)k
2k+1
 
1     -   1
3
    +   1
5
    -   1
7
    +  ...




In both the above equations (1) and (2), k takes values from 0 to ∞.
As you can see from equation (1), the denominator gets bigger within short range of terms which means the result of the division gets smaller and smaller which means the terms start contributing to the later part of π. Compare equation (2) with equation (1). The denominator in equation (2) increases extremely slowly so the result of the division still has significant impact on the digits closer to decimal place. So one must continue to compute much larger range of terms in equation (2). This is the reason Euler-Gregory-Madhava series converges so slowly.

Reference

Calculating pi to one million decimal places, Chris Hills,
URL : http://www.codeproject.com/KB/recipes/CRHpi.aspx


†   URL: http://milan.milanovic.org/math/english/pi/machin.html
§   The reason why I enjoyed doing this exercise so much is I would try very hard on my own to figure out the concepts and the implementation. Only when I was absolutely stuck I referred to Chris' article. I conveniently customized the number system in power of 10 because then the conversion from this number system to 10 at the time of printing value of π becomes trivial.


Downloads

Comments

  • There are no comments yet. Be the first to comment!

Leave a Comment
  • Your email address will not be published. All fields are required.

Top White Papers and Webcasts

  • With JRebel, developers get to see their code changes immediately, fine-tune their code with incremental changes, debug, explore and deploy their code with ease (both locally and remotely), and ultimately spend more time coding instead of waiting for the dreaded application redeploy to finish. Every time a developer tests a code change it takes minutes to build and deploy the application. JRebel keeps the app server running at all times, so testing is instantaneous and interactive.

  • Ever-increasing workloads and the challenge of containing costs leave companies conflicted by the need for increased processing capacity while limiting physical expansion. Migration to HP's new generation of increased-density rack-and-blade servers can address growing demands for compute capacity while reducing costly sprawl. Sponsored by: HP and Intel® Xeon® processors Intel, the Intel logo, and Xeon Inside are trademarks of Intel Corporation in the U.S. and/or other countries. HP is the sponsor …

Most Popular Programming Stories

More for Developers

Latest Developer Headlines

RSS Feeds