BeginPackage["ISP`Blocks`"]

Needs["NumericalMath`ComputerArithmetic`"]

BlockFast::usage = "BlockFast[n] computes the extension possible with a stack of n blocks"

BlockFloat::usage = "BlockFloat[n] computes the extension possible with a stack of n blocks by summing n floating-point numbers"

BlockPrecision::usage = "BlockPrecision[n, p] computes the extension possible with a stack of n blocks by summing n p-digit floating-point numbers"

BlockRat::usage = "BlockRat[n] computes the extension possible with a stack of n blocks by summing n rational numbers"

Begin["Private`"]

BlockFast = Function[n, (1./2.)*(N[EulerGamma] + PolyGamma[0., 1.+n])]

BlockFloat = Function[n, Module[{i,sum},
				i=1; sum=0.;
 				While[i<=n, sum=sum+1./(2.*i); i=i+1];
				sum]]

BlockRat = Function[n, Module[{i,sum},
				i=1; sum=0;
				While[i<=n, sum=sum+1/(2*i); i=i+1];
				sum]]

BlockPrecision = Function[{n,p}, Module[{i,sum},
			 Which [p <= 10,
				SetArithmetic[p, MixedMode->True];
				sum=ComputerNumber[0];
 				If[n>0,
                                   i=1;
                                   While[i<=n, sum=sum+1/(ComputerNumber[2] *
                                                          ComputerNumber[i]);
                                            i=i+1],
                                   i=-n;
                                   While[i>=1, sum=sum+1/(ComputerNumber[2] *
       							  ComputerNumber[i]);
                                            i=i-1]];
				sum,

			        p > 16,
				sum = 0;
 				If[n>0,
                                   i=1;
                                   While[i<=n, sum=sum+1/N[2*i,p];
                                            i=i+1],
                                   i=-n;
                                   While[i>=1, sum=sum+1/N[2*i,p];
                                            i=i-1]];
				sum,

			        True,
				sum = 0;
 				If[n>0,
                                   i=1;
                                   While[i<=n, sum=sum+1/(2.*i);
                                            i=i+1],
                                   i=-n;
                                   While[i>=1, sum=sum+1/(2.*i);
                                            i=i-1]];
				sum]]]

End[]

EndPackage[]
