Computing e

June 07, 2012
Home

I saw a job posting on Hacker News with an interesting problem which promptly sucked me in.

Send a bash program that calculates Euler’s number to N decimal places (where N is specified on the command line) to email.

I thought it looked like an interesting problem so I opened up Vim and started working on it in my language of choice, C, with the intention of porting it to Bash after I got it working.

The first step was determining how to calculate e programmatically, which brought me to Wikipedia. The infinite series at the top of the article seemed like a good place to start so I implemented it quickly and hackily.

#include <stdio.h>
int main(int argc, const char *argv[])
{
	double e = 0;
	double e_prev = 0;

	double n = 1;
	int i = 2;
	int places;
	int temp;

	if (argc < 2) {
		printf("Usage: %s n\nn - number of decimal places\n", argv[0]);
		return -1;
	}
	places = atoi(argv[1]);

	printf("2.");
	n = n*(i++);
	e = e+(1.0/n);
	while (places) {
		e_prev = e;
		n = n*(i++);
		e = e+(1.0/n);
		// if digit in tenths place is the same, output and shift left
		while(places && (int)(e*10) == (int)(e_prev*10)) {
			temp = (int)(e*10);
			printf("%d", temp);
			e *= 10;
			e_prev *= 10;
			n = n / 10;
			e -= temp;
			e_prev -= temp;
			places--;
		}
	}
	printf("\n");
	return 0;
}

That worked for the first 11 digits (2.71828182845), but failed on the 12th (2.718281828458 instead of 2.78281828459)

I thought that the issue would be the denominator of the infinite series being extremely large and overflowing, but that’s not the case. Before that happens, the double holding e loses precision after too many divide-by-10s.

I did some Googling (actually DuckDuckGo-ing, but that doesn’t have the same ring to it) and found a website describing how one person wrote a program to do exactly what I wanted. It referenced a paper on an finding digits of pi. I read through the paper and concluded that I am bad at reading and understanding academic papers, especially since the code snippets are written in Haskell, a language that I don’t know. I was able to read the blog post though and implement it in C (after spending way too much time figuring out the matrix math).

#include <stdio.h>
int main(int argc, const char *argv[])
{
	double e_low = 0;
	double e_high = 1;
	double m [2][2];
	int temp;
	
	int max;
	int i=3;

	if (argc < 2) {
		printf("Usage: %s n\nn - number of decimal places\n", argv[0]);
		return -1;
	}
	max = atoi(argv[1]);

	m[0][0] = .5;
	m[0][1] = 1;
	m[1][0] = 0;
	m[1][1] = 1;

	printf("2.");
	while (max--) {
		// keep multiplying terms until we converge on a digit
		while (((int)(e_low*10)) != ((int)(e_high*10))) {
			m[0][1] = m[0][0] + m[0][1];
			m[0][0] = m[0][0] * (1.0/i);
			i++;
			//m[1][0] = 0;
			//m[1][1] = 1;

			e_low = m[0][0]*1 + m[0][1];
			e_high = m[0][0]*2 + m[0][1];
		}
		temp = (int)(e_low*10);
		e_low = e_low*10;
		e_low = e_low-temp;
		e_high = e_high*10;
		e_high = e_high-temp;
		m[0][0] *= 10;
		m[0][1] *= 10;
		m[0][1] -= temp;
		printf("%d", temp%10);
		
	}
	printf("\n");
	return 0;
}

This algorithm works using a different infinite sequence, but it still has the same precision errors. After 15 digits, the first error occurs (2.7182818284590453 instead of 2.7182818284590452). I spent a long time going over the code and couldn’t find any errors that were causing my program to fail while it seemed to work for the guy who’s blog I was reading. I came to the conclusion that the author must have used BigInteger or another library.

I still really wanted to solve the problem, so I decided to give it one last shot. I decided to go back to the original problem statement, and I wrote it in Bash. It isn’t written using only Bash, but at least it gives me some closure. It works for up to 1,000,000 digits of e.

#!/bin/bash

if [[ ( $# < 1 ) || !( $1 =~ ^[0-9]+$ ) ]]; then
	echo "Usage: $0 n"
	echo "   n - number of decimal places of e to output"
	exit -1
fi

# output the first 2 digits ("2.")
n=$((2 + $1))

# grab first 1 million digits from internet, remove header and newlines, print
# first $n characters
curl -s http://apod.nasa.gov/htmltest/gifcity/e.1mil | tail -c +1298 | \
	tr -d '\n' | head -c +$n
echo

If I have time, and motivation, I may revisit this problem. The next step I’d try would be to use a fixed point library, such as the GNU MPFR.