printpi3.sd7

(********************************************************************)
(*                                                                  *)
(*  printpi3.sd7  Print 1000 digits of PI with Newtons formula      *)
(*  Copyright (C) 2006  Thomas Mertes                               *)
(*                                                                  *)
(*  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 2 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  *)
(*  MERCHANTABILITY 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, write to the           *)
(*  Free Software Foundation, Inc., 51 Franklin Street,             *)
(*  Fifth Floor, Boston, MA  02110-1301, USA.                       *)
(*                                                                  *)
(********************************************************************)


$ include "seed7_05.s7i";
  include "stdio.s7i";

# Newtons formula for PI is:
# PI / 2 = sum_n_from_0_to_inf(n! / (2 * n + 1)!!)

# This can be written as:
# PI / 2 = 1 + 1/3 * (1 + 2/5 * (1 + 3/7 * (1 + 4/9 * (1 + ... ))))

# This algorithm puts 2 * 1000 on the right side and computes everything from inside out.

const integer: SCALE is 10000;
const integer: MAXARR is 3500;
const integer: ARRINIT is 2000;

const proc: main is func
  local
    var integer: i is 0;
    var integer: j is 0;
    var integer: carry is 0;
    var array integer: arr is MAXARR times ARRINIT;
    var integer: sum is 0;
  begin
    for i range MAXARR downto 1 step 14 do
      sum := 0;
      for j range i downto 1 do
        sum := sum*j + SCALE*arr[j];
        arr[j] := sum rem pred(j*2);
        sum := sum div pred(j*2);
      end for;
      write(carry + sum div SCALE lpad0 4);
      carry := sum rem SCALE;
    end for;
    writeln;
  end func;