Mathematics : Finding Pi using Spigot algorithm

Irrational numbers are interesting bit of mathematics. Those numbers can't be written in fractions. Their decimal goes on like forever without repeating. Pi, e, Golden Ratio are some of well known irrational numbers. And Pi is the subject of interest here. "Pi is an infinite, non-repeating decimal -- meaning that every possible number combination exists somewhere in pi. Converted into ASCII text, somewhere in that infinite string of digits is the name of every person you will ever love, the date, time, and manner of your death, and the answers to all great questions of the universe." (Copied that one from internet :D )
If its true that irrational numbers goes on forever then building one of the Pi computation engine would be a tough thing. Then a bright light of hope from far googleeast revealed the almighty Spigot Algorithm. Spigot Algorithm is simple and only needs integer arithmetics to compute the value of Pi. You can check wiki about the algorithm here . Just follow the links on references. I read one from Rabinowitz, Stanley; Wagon, Stan (1995). The algorithm presented on the paper uses a mixed base representation of PI to calculate its value.
Firstly I used python code as protype to compute first 1000 digits. As people never use python for its speed, I shifted to C++. Using C++ code, it took near about 0.4s in my laptop(i3 dual core) to compute first 1000 digits which was happily far better than python code. Then I tried 10000 digits, the time was near about 18 seconds. Then 100000 digits. This one made me wait quite long time. Waiting was tough. While my laptop was busy calculating digits, I grabbed a sheet of paper and roughly calculated the time. Considering 2Ghz processor and total loop counts, 27 minutes was my conclusion. And I was pretty much sure that it would take more than that. To my surprise, the program stopped with time 27.906 minutes :D. However the program got final 3 digits wrong. A quick google showed it was a perfect valid condition for a spigot algorithm. They inherit such problems in general and are less likely to be used for calculating higher digits.
And here I present my code snippet I had built. First a quick python prototype.
def reform():
    i = -1;
    while(pirepo[i] >= 10):
        q = (int) (pirepo[i]/10)
        pirepo[i] %= 10;
        i -= 1
        pirepo[i] += q
    return
 
 
digits = 1000
bucketsize = (int)(10 * digits / 3) + 1
bucket = [2 for i in xrange(bucketsize)]
 
pirepo = []
 
for i in xrange(digits):
    bucket = [(bucket[j] * 10) for j in xrange(bucketsize)]
 
    quotent = 0
    for j in xrange(bucketsize):
        index = bucketsize - j - 1
        bucket[index] = ( bucket[index] + quotent * (index + 1))
        quotent = (int)(bucket[index] / (2 * index + 1))
        bucket[index] = bucket[index] % (2 * index + 1)
 
    bucket[0] = quotent % 10
    output = (int)(quotent / 10)
 
    pirepo.append(output)
 
    reform()
 
print pirepo
And implementation on C++
#include 
class PI{
private:
    const static int digits = 100000;
    const static int bucketSize = (10*digits/3) + 1;
 
    int piStore[digits];
    int bucket[bucketSize];
 
public:
    void calculate(){
        for(int i = 0; i < bucketSize ; i++)
            bucket[i] = 2;
 
        int q = 0;
        for (int i = 0; i < digits ; i++){
            q = 0;
            for (int j = bucketSize-1; j >= 0 ; j--){
                bucket[j] *= 10;
                bucket[j] = bucket[j] + (q * (j + 1));
                q = bucket[j] / (2 * j + 1);
                bucket[j] = bucket[j] % (2 * j + 1);
            }
            piStore[i] = q/10;
            bucket[0] = q%10;
 
        }
 
        revalidate();
 
        for ( int i = 0 ; i < digits ; i++)
            std::cout << piStore[i];
    }
    void revalidate(){
        int q;
        int index;
        for (int i = 0 ; i < digits ; i++){
            q = 0;
            index = 0;
            while (piStore[i-index]>9){
                q = piStore[i-index]/10;
                piStore[i-index] %= 10;
                index++;
                piStore[i-index] += q;
            }
        }
    }
};
int main(){
    PI piObj;
    piObj.calculate();
    return 0;
}
The code recipes provided here are my versions. There are far better versions available on internet, you can google them or you can simply try to build your own version. I suggest building your own version, its more fun. My journey to find Pi seems to end here. Happy with my 100,000 digits. Thanks to my friend who talked to me about the Pi. And finally, Happy coding to all of you .

No comments :

Post a Comment