#!/usr/local/bin/perl

# ----------------------------------------------------------------
# John Kerl
# kerl.john.r@gmail.com
# 2005-01-11
#
# This is an example of numerical series estimation.
#
# The Maclaurin series for exp(x) is sum from k = 0 to infinity of x^k/k!.  For
# x = 1.0, one expects e, i.e. 2.71828182845905 (to fourteen places, recalling
# that Perl does double-precision floating-point).
#
# Also, recall the identity pi^2/6 = sum from k = 1 to infinity of 1/k^2.  One
# expects pi = 3.14159265358979 (to fourteen places).
#
# This program prints:
#
#    0   1.00000000000000               ----------------
#    1   2.00000000000000 7.18282e-01   2.44948974278318 6.92103e-01 
#    2   2.50000000000000 2.18282e-01   2.73861278752583 4.02980e-01 
#    3   2.66666666666667 5.16152e-02   2.85773803324704 2.83855e-01 
#    4   2.70833333333333 9.94850e-03   2.92261298612503 2.18980e-01 
#    5   2.71666666666667 1.61516e-03   2.96338770103857 1.78205e-01 
#    6   2.71805555555556 2.26273e-04   2.99137649474842 1.50216e-01 
#    7   2.71825396825397 2.78602e-05   3.01177394784621 1.29819e-01 
#    8   2.71827876984127 3.05862e-06   3.02729785665784 1.14295e-01 
#    9   2.71828152557319 3.02886e-07   3.03950758956105 1.02085e-01 
#   10   2.71828180114638 2.73127e-08   3.04936163598207 9.22310e-02 
#   11   2.71828182619849 2.26055e-09   3.05748150670756 8.41111e-02 
#   12   2.71828182828617 1.72876e-10   3.06428781783393 7.73048e-02 
#   13   2.71828182844676 1.22857e-11   3.07007537189322 7.15173e-02 
#   14   2.71828182845823 8.14904e-13   3.07505691557136 6.65357e-02 
#   15   2.71828182845899 5.01821e-14   3.07938982603209 6.22028e-02 
#   16   2.71828182845904 2.22045e-15   3.08319302033945 5.83996e-02 
#   17   2.71828182845905 4.44089e-16   3.08655802575371 5.50346e-02 
#   18   2.71828182845905 4.44089e-16   3.08955643496978 5.20362e-02 
#   19   2.71828182845905 4.44089e-16   3.09224505230070 4.93476e-02 
#   20   2.71828182845905 4.44089e-16   3.09466952411370 4.69231e-02 
#   21   2.71828182845905 4.44089e-16   3.09686694994393 4.47257e-02 
#   22   2.71828182845905 4.44089e-16   3.09886779322221 4.27249e-02 
#   23   2.71828182845905 4.44089e-16   3.10069730139518 4.08954e-02 
#   24   2.71828182845905 4.44089e-16   3.10237657635981 3.92161e-02 
#   25   2.71828182845905 4.44089e-16   3.10392339170058 3.76693e-02 
#   26   2.71828182845905 4.44089e-16   3.10535282394625 3.62398e-02 
#   27   2.71828182845905 4.44089e-16   3.10667774541646 3.49149e-02 
#   28   2.71828182845905 4.44089e-16   3.10790921281339 3.36834e-02 
#   29   2.71828182845905 4.44089e-16   3.10905677641032 3.25359e-02 
#   30   2.71828182845905 4.44089e-16   3.11012872814126 3.14639e-02 
#   31   2.71828182845905 4.44089e-16   3.11113230222817 3.04604e-02 
#   32   2.71828182845905 4.44089e-16   3.11207383861109 2.95188e-02 
#   33   2.71828182845905 4.44089e-16   3.11295891698571 2.86337e-02 
#   34   2.71828182845905 4.44089e-16   3.11379246743573 2.78002e-02 
#   35   2.71828182845905 4.44089e-16   3.11457886229313 2.70138e-02 
#   36   2.71828182845905 4.44089e-16   3.11532199284004 2.62707e-02 
#   37   2.71828182845905 4.44089e-16   3.11602533369232 2.55673e-02 
#   38   2.71828182845905 4.44089e-16   3.11669199711265 2.49007e-02 
#   39   2.71828182845905 4.44089e-16   3.11732477904398 2.42679e-02 
#
# This illustrates the rapid convergence of the first series (due to the
# factorial in the denominator) vs. the slow convergence of the second.
# ----------------------------------------------------------------

$kmax = 20;
$kmax = $ARGV[0] if (@ARGV);
$x    = 1.0;

# The exp series starts at k=0, but the pi^2/6 series starts at k=1.
# So, initialize the exp sum to already have the value for k=0.
$expsum = 1.0;
$kfact  = 1;
$xpower = 1.0;

$pisum = 0.0;

$pi = 4.0 * atan2(1.0, 1.0);
$e  = exp(1.0);

$k = 0;
printf "%4d %18.14f               ----------------\n", $k, $expsum, $pisum;

for ($k = 1; $k < $kmax; $k++) {
	$kfact  *= $k;
	$xpower *= $x;
	$expterm = $xpower / $kfact;
	$expsum += $expterm;

	$piterm = 1.0 / $k / $k;
	$pisum += $piterm;

	$approxpi = sqrt(6.0 * $pisum);

	printf "%4d %18.14f %9.5e %18.14f %9.5e \n",
		$k,
		$expsum, abs($expsum - $e),
		sqrt(6.0*$pisum), abs($approxpi - $pi);
}