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.