Approximating e - Getting as many precision digits as possible

Question

I have been trying to approximate e using series representation to get as many precision digits as possible using the code below, but no matter how many terms I compute, the number of precision digits seems to remain the same. ie:

2.71828198432922363281250000000000000000000000000000

Is it my approach that's wrong? Here is the code:

  1 #include <stdio.h>
  2 #include <iostream>
  3 #include <math.h>
  4 using namespace std;
  5 
  6 float factorial (float a)
  7 {
  8         if (a > 1)
  9         {
 10                 return (a * factorial (a-1));
 11         } else
 12         {
 13                 return 1;
 14         }
 15 }
 16 
 17 int main()
 18 {
 19         float sum  = 0;
 20         int range=100000;
 21 
 22         for (int i=0; i<=range;i++)
 23         {
 24                 sum += pow(-1,i)/factorial(i);
 25         }
 26         sum = pow(sum,-1);
 27         printf("%4.50f\n", sum);
 28 } 

Answer 1

To get more exact digits, you should write your on data class which store more digit, say, 1000 digits. The hardest part is to wirte the +, -, *, / operations.

If what you want is just to experiment with the math formula, you can choose another language, such as Python. It has data types like Decimal, Fraction that can do more precise calculating.

I love math so I do write a python script to test the formula:

from decimal import Decimal, getcontext
prec = 100
getcontext().prec = prec

fac = Decimal(1)/2
sum = Decimal(0)
eps = Decimal(0.1)
eps = eps**prec

realE = '2.71828182845904523536028747135266249775724709369995'

i = 3
while 1:
    ds = fac - fac/i
    sum += ds
    if ds < eps: break
    fac /= i * (i+1)
    i += 2

print 'After %d iteration:' % i
print realE
print str(1/sum)

Here's the result:

After 71 iteration:
2.71828182845904523536028747135266249775724709369995
2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166429

Answer 2

You are hitting the limit where the number added is much smaller than the sum, and because floats are basically rational numbers, this is cut off. Here is a good read about the subtleties of floating point numbers

Example:

 12345.123
+    0.0001
------------
 12345.123

if only the first 8 digits are saved in a number.

A easy fix would be to iterate from range to 0, so that you start your sum with the small numbers, and to keep track of the lost digits. As an example:

sum0 = 12345.123
b0   =     0.0001
sum1 = sum0 + b0 # 12345.123
diff1 = (sum1 - sum0) - b0 # 0.0001
# the result you want is sum1 + diff1

# keep on iterating and summing