/*
** This program was inspired by the article "The Impossible Dream:
** Computing e to 116,000 places with a Personal Computer" by
** Stephen Wozniak, published in Byte Magazine Vol 6, Issue 6
** (June 1981, pp392).
**
** Essentially, the following program implements a partial expansion
** of the equation: e = 1 + (1/1!) + (1/2!) + (1/3!) + ... + (1/n!)
** by inverting the equation, such that it reads...
** e = (((((...(((((1/n) + 1) / (n-1)) + 1) ... / 3) + 1) / 2) + 1) / 1) + 1
** 
** I originally wrote this program in Z80 Assembler code, to run under
** CP/M (Dec 1986), and subsequently have recoded the program in
** 370 Assembler to run under MVS and VSE, and C to run under MSDOS
** and Linux
**
** Lew Pitcher, 1999
*/
#include <stdio.h>

#define	I	1		/* number of bytes in integer portion	*/
#define	X	2078		/* number of bytes in fraction portion	*/
#define	N	1776		/* number of divide iterations		*/
#define	P	5000		/* number of multiply iterations	*/

char e[I+X];			/* holds fixed point binary value of e	*/

int main(void)
{
	void divide(char *, unsigned, unsigned),
	     multiply(char *, unsigned, unsigned);

	unsigned n,		/* divisor for calculation loop		*/
		 p;		/* precision for print loop		*/

	e[0] = 1;		/* start the series at 1.anything	*/
	for (n = N; n > 0; --n)
	{
		divide(e,n,I+X);
		e[0] += 1;
	}

	fputs(" e == ",stdout); fputc(e[0]+'0',stdout); fputc('.',stdout);
	for (p = 0; p < P; ++p)
	{
		e[0] = 0;
		multiply(e,10,I+X);
		fputc(e[0]+'0',stdout);
		if (p%5 == 4)
			if (p%50 == 49)
				fputs("\n\t",stdout);
			else
				fputc(' ',stdout);
	}
	fputc('\n',stdout);

	return 0;
}

void divide(char *dvd, unsigned dvr, unsigned prc)
{
	char *Wk_q;
	unsigned bit, rmdr;

	rmdr = 0;
	for (Wk_q = dvd; Wk_q < dvd+prc; ++Wk_q)
	{
		for(bit = 0; bit < 8; ++bit)
		{
			rmdr *= 2;
			if (*Wk_q & 0200) rmdr += 1;
			*Wk_q <<= 1;
			if (rmdr >= dvr)
			{
				*Wk_q += 1;
				rmdr -= dvr;
			}
		}
	}
}

void multiply(char *mcd, unsigned mpy, unsigned prc)
{
	char *Wk_r;
	unsigned bit, ovflo;

	ovflo = 0;
	for (Wk_r = mcd+prc-1; Wk_r >= mcd; --Wk_r)
	{
		for (bit = 8; bit > 0; --bit)
		{
			if (*Wk_r & 1) ovflo += mpy;
			*Wk_r >>= 1; *Wk_r &= 0177;
			if (ovflo & 1) *Wk_r = *Wk_r | 0200;
			ovflo /= 2;
		}
	}
}