Conference rusure::math

    P.S.
    I CAN NOT USE ANY SYSTEM SERVICES ONLY PASCAL.

I have a copy of a paper by A. H. J. Sale, describing an ALGOL procedure 
for this problem, dated 1968. I think it's from an IEEE publication, or
possibly Computer Journal - unfortunately the paper does not have the 
journal name or exact date on it! Send me your SnailMail address for a 
hardcopy.

Lynn Yarbrough

    Re .0:
    
    Keep more than one digit in an array element.  (You will need to write
    a routine to output array elements with leading zeroes.)  Select a
    large base but leave yourself room to handle carrying.  You only need
    two operations:  divide (set remainder to zero, for each element
    left-to-right:  add remainder to element, divide by divisor, put
    quotient back in element, keep remainder for next element), and add
    (set carry to zero, for each element right-to-left: add two elements
    and carry, if over base set carry to one and subtract base otherwise
    set carry to zero, return result to element).
    
    You could also keep track of which elements have become zero, to avoid
    dividing or adding zeroes unnecessarily.
    
    
    				-- edp

program e(input, output);

const
    n = 450;

var
    a : array[1..n] of integer;
    i, j, q : integer;

begin

    a[1] := 2;
    for i := 2 to n do a[i] := 1;

    for i := 1 to 1000 do begin
	q := 0;
	for j := n downto 1 do begin
	    a[j] := a[j]*10+q;
	    q := a[j] div j;
	    a[j] := a[j] rem j;
	    end;
    write(q:1);
    if (i mod 64) = 0 then writeln;
    end;

end.

    Also, see note 92 for a similar algorithm that computes digits of pi.

Re .4:

Is that Sale's algorithm, an implementation of edp's "BIGNUM" suggestion,
or something else entirely?
				/AHM/THX

The algorithm in .4 appears to be the same as Sale's.

Perhaps of help is the following:

e = 1 + 1
	--------------------------------------
	2 + 1
	    --------------------------------------
	    1 + 1
		--------------------------------------
		1 + 1
		    --------------------------------------
		    4 + 1
			-------------------------------------
			1 + 1
		    	    ------------------------------------
		    	    1 + 1
				-----------------------------------
				6 + 1
				    ----------------------------------
				    1 + 1
					---------------------------------
					1 + 1
					    -------------------------------
					    8 + 1
						------------------------------
						1 + 1
						    ----------------------------
						    1 + 1
							------------------------
							10 . . .

/Eric

One thing that's always worried me about published "1000 digits of pi"
and "1000 digits of e" etc. is, are they correct ?

Without careful checking, it seems to me quite likely that computer
roundoff errors can cause unexpected glitches in the output.

Any of you that claim to have working programs that produce 1000 digits
of e, perhaps you can program in a check to verify it.

For instance, use some totally different method to take the natural
log of the result and see if you get something that is closer to 1
than if you take the natural log of, say, all 2000 values produced by
taking your answer and modifying each digit in turn to one value too
high and one value too low.

/Eric

re .9:  Most of the algorithms I've seen use integer arithmetic, so roundoff
errors don't factor into the algorithm at all - only overflows.  The algorithm
is chosen so that overflows are impossible, carrys from word to word are done
before the word range of the machine's integers are exceeded.  Another way to
think of this is the calculations are performed using base n numbers, where n
is large and dependant on the word size of the machine.

Similarly, our decimal counting system can represent any number to any degree
of accuracy given enough digits, even though the range of each digit is
restricted to the range 0-9. (equivalent to a machine whose "word" is between
3 and 4 bits long!)

-Mike

I don't care how much you want to flail your arms and shout and jump up
and down and guarantee that your huge numbers are accurate.

I STILL don't trust them until I see a verification check by an independent
method.

Back in engineering school, in order to get credit on an exam problem,
we had to include a "check".  For instance, if the problem was to calculate
I by dividing V/R, we could get full credit by multiplying I*R to see if
we got V again.

/Eric

    Re .11:
    
program e(input, output);

const
    n = 450;

var
    a : array[1..n] of integer;
    i, j, q : integer;                

begin
    
    {(0) Let maximum = 450.}
    a[1] := 2;
    for i := 2 to n do a[i] := 1;
    {(1) Therefore a[i] <= j+1 for 0 < i <= n.}
    
    for i := 1 to 1000 do begin
	q := 0;
        {(2) Therefore q <= 2*maximum.}
	for j := n downto 1 do begin
            {(3) From (2) or (7), q <= 2*maximum.}
            {(3.5) From (1) or (8), a[j] <= j+1.}
	    a[j] := a[j]*10+q;
            {(4) Therefore a[j] <= 10*(j+1)+2*maximum.}
	    q := a[j] div j;
            {(5) From (4), q <= 10*(j+1)/j+2*maximum/j.}
    	    {(6) From (5), if j >= 2, q <= 10*(j+1)/j+maximum.}
            {(7) 10*(j+1)/2 is at most 15, which is less than maximum,
                 so q <= maximum+maximum, so q <= 2*maximum if j >=2.}
            {Note that (3) is arrived at only from (2) or from (7) when
             j was greater than or equal to 2 at (7).}
	    a[j] := a[j] rem j;
    	    {(8) Therefore a[j] < j, so a[j] <= j+1.}
	    end;
    write(q:1);
    if (i mod 64) = 0 then writeln;
    end;

end.
    
    It can now be seen that no number, i, j, n, q, a[i], or intermediate
    expression ever exceeds 10*(j+1)+2*maximum for the greatest j, 450,
    which is 5410, well less than the maximum integer of VAX computers.
    Inspection also shows no number becomes negative.


				-- edp

	People who calculate pi to millions of decimal places usually
	do an independent check before publishing the digits.  That's
	because there may have been some computer errors during the
	computation, or some errors in the algorithm or the program
	that implements the algorithm.  Typically, a careful analysis
	of the algorithm is done to determine the effect of accumulated
	round-off errors.

	Let's consider the calculation of e with the following algorithm.
	It's based on the series e = 1 + 1/1! + 1/2! + 1/3! + ...

	x := 10**1000;		{ a multi-precision package is used }
	y := 10**1000;
	for n := 1 to 450 do
	    begin
	    y := y div n;	{ the result is truncated to an integer }
				{ y is roughly 1/n! * 10**1000 }
	    x := x + y;		{ exact; x accumulates the sum }
	    end;
	{ Now x is roughly e * 10**1000 }

	The algorithm is fairly straight-forward, and the multi-precision
	package itself has been extensively tested.

	The calculated value of 1/n! * 10**1000 is always too small (for n > 2).
	We assert that for n > 2, 1/n! * 10**1000 - 1 < y(n) < 1/n! * 10**1000.
	This is easily demonstrated for n=3.  For induction, we must show that

	     1/n! * 10**1000 - 1 < y(n) < 1/n! * 10**1000

	implies

	    (1/n!*10**1000)/(n+1)-1 <= [y(n)/(n+1)] < (1/n!*10**1000)/(n+1)

	where [x] indicates the floor of x.

	The right-hand inequality is easy:

	    [y(n)/(n+1)] <= y(n)/(n+1) < (1/n!*10**1000)/(n+1)

	The left-hand inequality is pretty simple.  We first note that

	    x <= y  implies  x/b - (b-1)/b <= [y/b], for positive integers y,b.

	So we have:

	    1/n! * 10**1000 - 1 <= y(n), (given)
	    (1/n! * 10**1000 - 1)/(n+1) - n/(n+1) <= [y(n)/(n+1)]
	    1/(n+1)! * 10**1000 - 1 <= [y(n+1)]


	This tells us that y is indeed close to 1/n!, and gives an estimate
	of the error.  It's easy to show that y is always >= 0, and so when
	n >= 450, the inequality indicates y must be zero (450! > 10**1000).

	Each iteration of the loop may accumulate errors in x.  Since each
	value of y has an error of 1 or less, the accumulated error in x
	must be less than 450 (the values for n=1 and n=2 are computed exactly).
	After the computation, we know that

		x <= e * 10**1000 <= x + 450

	We still haven't determined the number of digits that are correct.
	This can be determined by calculating x+450, and comparing the digits
	with x; where they are the same (before the first difference), we are
	assured that digit is correct in the infinite decimal expansion of e.