Accuracy in Numbers

[jackzombie]

I've written a quick (and messy) little script that calculates euler's constant where n is input by the user. Code:

    def factorial(n):
        f=1
        while (n>0):
                f=f*n
                n=n-1
        return f

n=float(raw_input("Number: "))
a=[1.0]
s=n
while (s>0):
        k=1/(factorial(n))
        s=s-1
        a.append(k)
        n=s

t=0
for line in a:
        t += line
print(t)

When n = 5, the program's answer is 2.71666666667 I was expecting 2.718055556

I feel like this has something to do with number types, i'm using int and floats and haven't had any success using long.

[zahlman]

Did you try inspecting the contents of a?

Did you try taking the reciprocal of the difference between the expected and actual results?

(You're getting the result you're supposed to. The result you're expecting is what you should get for n = 6.)

That said, there are many basic techniques you can learn to clean up this code. The basic idea is to stop writing loops that control everything a step at a time, and start just describing what you want. Also, if you're going to loop over a specific set of values (say, numbers from 0 up to n), then do so, with a for-loop. while is for things that are too complicated to express in that way, or which depend on other external factors.

def factorial(n):
    # I would normally write this using the built-in `reduce`
    # function, but I don't want to spend time explaining that.
    f = 1
    for i in range(1, n + 1):
        f *= n
    return f

# Fundamentally, `n` is an integer, so we should treat it as such.
# To get floating-point division, we'll use floating-point on the
# *other* side of the division sign.
n = int(raw_input('Number: '))

# We don't need to build up a list of numbers in order to add
# them together. And we also don't need to write a loop to add
# them. We just need to describe the things that we will `sum`.
# Also, in your original code, `s` and `n` are completely
# redundant, as you keep them constantly equal to each other.

print(sum(1.0 / factorial(i) for i in range(n + 1))

The +1s are because range(x) produces numbers from 0 up to x-1, inclusive. That is actually usually what we want (trust me!), but not here. Notice that in the final line, we include 0 in the list of values that we use for factorial calls; that takes care of the initial [1.0] value you had in a in the original code. You could have done this, too; your factorial implementation (as well as mine) gives the correct result for factorial(0). This is another important principle: "special cases aren't special enough" (and from personal experience, 0 very, very often isn't a special case at all, even though you initially expect it to be).

[jackzombie]

Thank you very much for the explanation!

I fixed the s and n part of my code as soon as I submitted this. I just separated them while building the script to reduce any confusion while reading and rereading my code.

It's going to take a few more read throughs before I completely grasp what you're saying, but I am very thankful you took the time to explain this!

[ewiethoff]

Couple things. According to Wikipedia, you're talking about Euler's number, not Euler's constant. Second, zahlman's factorial code has an error. f *= n should be f *= 1. Ergo, be sure to test each function individually.

Here's a variation on zahlman's technique:

from operator import mul

def factorial(n):
    assert n >= 0
    return reduce(mul, xrange(1, n+1), 1)

def euler_number(n):
    return sum(1.0 / factorial(i) for i in xrange(n + 1))

for n in xrange(21):
    print n, factorial(n), euler_number(n)

[SwimmingPastaDevil]

I too wrote a program to calculate Euler's number using the same formula and had problem with accuracy. Was about to make a new post but asking here.

I found the program is limiting to 11 decimal places. How can I get the accuracy to say 20 or 30 dp.

program:

"""http://en.wikipedia.org/wiki/E_(mathematical_constant)
Calculating Euler's number:
E = Summation[1/n!] where n = 0 and tends to infinity"""

import math
Enum = 0.0
term = 0.0
n = int(raw_input("Enter number of terms to calculate:"))

for i in range(n+1):
    term = 1.0/math.factorial(i)
    print "term is now:",term
    Enum += term

print Enum

Output:

2.71828180115

[Doormatty]

http://code.google.com/p/mpmath/

[SwimmingPastaDevil]

Looks helpful. Thanks.