Gregory’s series

Another short ALGOL 68 program that I wrote last year while trying to get to grips with the language again, and as before this was based on a existing PASCAL program that I’d written previously.

The original program used Gregory’s Series to find Pi but this converges very slowly and while it is useful for using up spare clock cycles it isn’t much good for finding Pi.

After doing a bit of research I found an alternative related series that converges much faster and modified the original program to use that instead.

PROGRAM series

{ Program to calculate Pi - 05 Aug 12

  Copyright (C) 2012  - MEJT

  This program is free software: you can redistribute it and/or modify it
  under the terms of the GNU General Public License as published by the Free
  Software Foundation, either version 3 of the License, or (at your option)
  any later version.

  This program is distributed in the hope that it will be useful, but 
  WITHOUT ANY WARRANTY; without even the implied warranty of ERCHANTABILITY
  or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  for more details. 

  You should have received a copy of the GNU General Public License along
  with this program.  If not, see http://www.gnu.org/licenses/ 

  Uses a Gregory's series expansion to find Pi, an example of which is the
  Leibnitz Formula.

             Pi        1     1     1     1
            ---- = 1 - -  +  -  -  -  +  -  - ...
              4         3     5     7     9

  This series is very very slow to converge so the following series
  expansion has been used.

          1        Pi      1            1                   1
  atan(-------)  = -- = ------- - --------------- + ----------------- - ...
       sqrt(3)      6   sqrt(3)   3 * 3 * sqrt(3)   5 * 3^2 * sqrt(3)

  By expanding this and simplifying we get a series that converges much
  faster and appears to yeld aproximatly 16 decimal places in just 31
  iterations. 

          Pi              1         1           1
   ------------- = 1 - ----- + --------- - ------------- + ...
    2 * sqrt(3)        3 * 3   5 * 3 * 3   7 * 3 * 3 * 3

  When executed this program prints successive approximations for Pi.

    1                   1  3.4641016151377546  
    2                   9  3.0792014356780041  
    3                  45  3.1561814715699542  
    4                 189  3.1378528915956804  
    5                 729  3.1426047456630847  
    6                2673  3.1413087854628836  
    7                9477  3.1416743126988377  
    8               32805  3.1415687159417843  
    9              111537  3.1415997738115058  
   10              373977  3.1415905109380801  
   11             1240029  3.1415933045030816  
   12             4074381  3.1415924542876464  
   13            13286025  3.1415927150203798  
   14            43046721  3.1415926345473139  
   15           138706101  3.1415926595217137  
   16           444816117  3.1415926517339976  
   17          1420541793  3.1415926541725754  
   18          4519905705  3.1415926534061652  
   19         14334558093  3.1415926536478261  
   20         45328197213  3.1415926535714034  
   21        142958160441  3.1415926535956350  
   22        449795187729  3.1415926535879335  
   23       1412147682405  3.1415926535903866  
   24       4424729404869  3.1415926535896037  
   25      13839047287569  3.1415926535898540  
   26      43211719081593  3.1415926535897738  
   27     134718888901437  3.1415926535897996  
   28     419407861674285  3.1415926535897913  
   29    1303977169932777  3.1415926535897939  
   30    4049192264528097  3.1415926535897931  
   31   12559359057773589  3.1415926535897933  
   32   38913423965888661  3.1415926535897932 

}

BEGIN

   REAL term := 1.0, factor;
   REAL number := 1.0, sum := 0.0;
   INT count := 0, limit := 31;

   factor := 2 * sqrt (3.0);

   WHILE (count <= limit) DO
      term := (2 * count + 1) * 3.0 ^ (count);
      IF ODD count THEN
         sum := sum - (1.0 / term)
      ELSE
         sum := sum + (1.0 / term)
      FI;
      count +:= 1;
      write ((fixed (count, -3, 0), "  ", fixed (term, -18, 0), "  "));
      write ((fixed (sum * factor, -18, 16), "  ",  newline))
   OD

END

FINISH
This entry was posted in Programming and tagged . Bookmark the permalink.

Leave a Reply

Please log in using one of these methods to post your comment:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s