R Programming: How to Get Euler's Number

R programming: How do I get Euler's number?

The R expression

exp(1)

represents e, and

exp(2)

represents e^2.

This works because exp is the exponentiation function with base e.

How to get Euler–Mascheroni's constant in R?

Mathmatically we can write the constant of interest as the negative of the derivative of the gamma function evaluated at 1. R has the derivative of the gamma function as digamma so it's just a matter of plugging this in.

-digamma(1)
#[1] 0.5772157

Calculating the Euler's Number gives a limited number of decimals on C#

A fairly simple improvement is to perform the addition starting from the smallest terms first, so they are similar in size, and the processor can perform the addition without as much loss of precision.

Note: factorial(50) is 50 * factorial(49), and the last two terms are 1/factorial(49) + 1/factorial(50). Apply some algebra to get 1/factorial(49) + 1.0/50/factorial(49), which is (1 + 1.0/50) / factorial(49)

Say you only calculate the numerator, and keep track of what factorial appears in the denominator without calculating it. Now you have two very nice properties:

1. You never have to calculate numbers that overflow (like factorial(i)) and
2. Whatever rounding error there is in the update equation, doesn't matter to the final answer, because it's an error in a term that's going to get divided some more and become even smaller

That leads to the following code:

double accumulator = 1.0;

for( int i = 50; i > 0; --i )
{
    accumulator = 1.0 + accumulator / i;
}

Demo: https://rextester.com/FEZB44588

---

Extension to use .NET's BigInteger class allows us to have many more digits of precision

BigInteger scale = BigInteger.Pow(10, 60);
BigInteger accumulator = scale;

for( int i = 75; i > 0; --i )
{
    accumulator = scale + accumulator / i;
}

result (insert the decimal point):

2.718281828459045235360287471352662497757247093699959574966967

first 50 decimal places from Wikipedia:

2.71828182845904523536028747135266249775724709369995...

Note that Wikipedia's verbiage is slightly wrong, this is not the value rounded to 50 decimal places, these are the first 50 decimal digits of a sequence that continues

Euler's number with stop condition

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.

Approximate Euler's Number in script

Firstly I would rewrite your function for finding the factorial to the simpler

function n = factorial(n)
    n = prod(1:n);
end

The loop in your question is unnecessary as you never use the loop variable i. I wouldn't use this function though for my solution as it can be quite slow because you have to compute redundant information at each loop iteration.

If you still want to use a for loop you need to rewrite it to

function f = factorial(n)
    f = 1; % 0 factorial
    for i = 1:n
        f = f * i;
    end
end

---

You can use the natural logarithm and the rules of logs to determine a very accurate value of e with which you can compare against. The value of e you can check against is given by x^(1 / log(x)) where x can be any positive real number other than 1, like 2. We can see this in

                                                                

Now how do we check that we've computed a value of e to 10 decimal places of accuracy. Well given that b from above is a very accurate representation of e we can compare against it to determine when we've reached an accurate solution

x = 2; % Any positive number other than 1
c = x^(1 / log(x));
...
if (abs(e - c) < 1e-10)
    break;
end

In my solution e is the approximate value I've computed with the infinite sum. Note: the absolute value is taken to prevent false positives when e - c is a negative number.

---

Now, an efficient method for computing the infinite sum. We can exploit how the factorial is computed to not have to compute it during each iteration hugely improving the efficiency. Firstly, we need a sum variable, e in my case, to keep track of our approximate solution. Then we need another variable to keep track of the factorial, f in my case. As 0 is a funny case we'll start off with it

e = 0;
f = 1; % 0 factorial

e = e + 1 / f;

and now we have the first element in our infinite sum. Next we can use the infinite sum to compute a more accurate approximate to e. The factorial can be updated during each iteration with f = f * n; leading to

for n = 1:inf
    f = f * n; % Compute new factorial
    e = e + 1 / f; % Infinite sum
    ...
end

Now putting that altogether produces

x = 2; % Any positive number other than 1
c = x^(1 / log(x));

e = 0;
f = 1; % 0 factorial

e = e + 1 / f;

for n = 1:inf
    f = f * n; % Compute new factorial
    e = e + 1 / f; % Infinite sum
    if (abs(e - c) < 1e-10)
        break;
    end
end

Calculating Euler's constant with OCaml

As I recall, scheme has a "numerical tower" that tries to keep you from losing accuracy in numeric computations.

OCaml doesn't have any fancy automatic handling of numbers. Your code is using type int, which is a fixed size integer (31 or 63 bits in the usual implementations). Thus, your expression 1 / factorial accuracy will be 0 in almost all cases and values for factorial accuracy will be unrepresentable for all but the smallest values. The value of factorial 100 will be 0 because it is a multiple of 2^63:

# let rec fact n = if n < 2 then 1 else n * fact (n - 1);;
val fact : int -> int = <fun>
# fact 100;;
- : int = 0

There are no floats in your code, and hence there couldn't possibly be any mixing. But OCaml doesn't allow mixing in the first place. It's a strongly typed language where int and float are two different types.

Here is your code converted to use floats:

let rec factorial n =
    if n = 0.0 then 1.0 else n *. factorial (n -. 1.0) 

let rec find_e accuracy sum =
    if accuracy = 0.0 then
        sum +. 1.0
    else
        find_e (accuracy -. 1.0) (sum +. 1.0 /. factorial accuracy)

let result = find_e 100.0 0.0

If I copy/paste this into the OCaml REPL (also known as "toplevel") I see this:

val factorial : float -> float = <fun>
val find_e : float -> float -> float = <fun>
val result : float = 2.71828182845904509

As a side comment, the equality comparison operator in OCaml is =. Don't use == for comparisons. It's a completely different operator:

# 1.0 == 1.0;;
- : bool = false

---

Related Topics