Euler's number with stop condition

Question

original outdated code: Write an algorithm that compute the Euler's number until

My professor from Algorithms course gave me the following homework:

Write a C/C++ program that calculates the value of the Euler's number (e) with a given accuracy of eps > 0. Hint: The number e = 1 + 1/1! +1/2! + ... + 1 / n! + ... = 2.7172 ... can be calculated as the sum of elements of the sequence x_0, x_1, x_2, ..., where x_0 = 1, x_1 = 1+ 1/1 !, x_2 = 1 + 1/1! +1/2 !, ..., the summation continues as long as the condition |x_(i+1) - x_i| >= eps is valid.

As he further explained, eps is the precision of the algorithm. For example, the precision could be 1/100 |x_(i + 1) - x_i| = absolute value of ( x_(i+1) - x_i )

Currently, my program looks in the following way:

#include<iostream>
#include<cstdlib>
#include<math.h>

// Euler's number

using namespace std;

double factorial(double n)
{
    double result = 1;
    for(double i = 1; i <= n; i++)
    {
        result = result*i;
    }
    return result;
}

int main()
{
    long double euler = 2;

    long double counter = 2;
    long double epsilon = 1.0/1000;
    long double moduloDifference;

    do
    {
        euler+=  1 / factorial(counter);
        counter++;
        moduloDifference = (euler + 1 / factorial(counter+1) - euler);
    } while(moduloDifference >= epsilon);
    printf("%.35Lf ", euler );
    return 0;
}

Issues:

1. It seems my epsilon value does not work properly. It is supposed to control the precision. For example, when I wish precision of 5 digits, I initialize it to 1.0/10000, and it outputs 3 digits before they get truncated after 8 (.7180).
2. When I use long double data type, and epsilon = 1/10000, my epsilon gets the value 0, and my program runs infinitely. Yet, if change the data type from long double to double, it works. Why epsilon becomes 0 when using long double data type?
3. How can I optimize the algorithm of finding Euler's number? I know, I can rid off the function and calculate the Euler's value on the fly, but after each attempt to do that, I receive other errors.

Answer 1

One problem with computing Euler's constant this way is pretty simple: you're starting with some fairly large numbers, but since the denominator in each term is N!, the amount added by each successive term shrinks very quickly. Using naive summation, you quickly reach a point where the value you're adding is small enough that it no longer affects the sum.

In the specific case of Euler's constant, since the numbers constantly decrease, one way we can deal with them quite a bit better is to compute and store all the terms, then add them up in reverse order.

Another possibility that's more general is to use Kahan's summation algorithm instead. This keeps track of a running error while it's doing the summation, and takes the current error into account as it's adding each successive term.

For example, I've rewritten your code to use Kahan summation to compute to (approximately) the limit of precision of a typical (80-bit) long double:

#include<iostream>
#include<cstdlib>
#include<math.h>
#include <vector>
#include <iomanip>
#include <limits>

// Euler's number

using namespace std;

long double factorial(long double n)
{
    long double result = 1.0L;
    for(int i = 1; i <= n; i++)
    {
        result = result*i;
    }
    return result;
}

template <class InIt>
typename std::iterator_traits<InIt>::value_type accumulate(InIt begin, InIt end) {
    typedef typename std::iterator_traits<InIt>::value_type real;
    real sum = real();
    real running_error = real();

    for ( ; begin != end; ++begin) {
        real difference = *begin - running_error;
        real temp = sum + difference;
        running_error = (temp - sum) - difference;
        sum = temp;
    }
    return sum;
}

int main()
{  
  std::vector<long double> terms;
  long double epsilon = 1e-19;

  long double i = 0;
  double term;
 
  for (int i=0; (term=1.0L/factorial(i)) >= epsilon; i++)
    terms.push_back(term);

  int width = std::numeric_limits<long double>::digits10;

  std::cout << std::setw(width) << std::setprecision(width) << accumulate(terms.begin(), terms.end()) << "\n";
}

Result: 2.71828182845904522

In fairness, I should actually add that I haven't checked what happens with your code using naive summation--it's possible the problem you're seeing is from some other source. On the other hand, this does fit fairly well with a type of situation where Kahan summation stands at least a reasonable chance of improving results.