0 views

Skip to first unread message

Nov 20, 2007, 2:41:02 AM11/20/07

to pytho...@python.org

I'm writing a demo of the infinite series

x**0/0! + x**1/1! + x**2/2! + x**3/3! + ... = e**x (x is non-negative)

It works OK for many x, but for many the loop doesn't break. Is there

a way to get it to break where I want it to, i.e., when the sum

equals the limit as closely as the precision allows?

Here's what I have:

======= series_xToN_OverFactorialN.py ==========================

#!/usr/bin/env python

#coding=utf-8

# series_xToN_OverFactorialN.py limit is e**x from p.63 in The

Pleasures of Pi,e

from mpmath import mpf, e, exp, factorial

import math

import time

precision = 100

mpf.dps = precision

n = mpf(0)

x = mpf(raw_input("Enter a non-negative int or float: "))

term = 1

sum = 0

limit = e**x

k = 0

while True:

k += 1

term = x**n/factorial(n)

sum += term

print " sum = %s k = %d" % (sum, k)

print "exp(%s) = %s" % (x, exp(x))

print " e**%s = %s" % (x, e**x)

print

if sum >= limit:

print "math.e**%s = %f" % (x, math.e**x)

print "last term = %s" % term

break

time.sleep(0.2)

n += 1

"""

Output for x == mpf(123.45):

sum =

410822093109668964239148908443317876138879647013995774.2951431466270782257597573259486687336246984942

k = 427

exp(123.45) =

410822093109668964239148908443317876138879647013995774.2951431466270782257597573259486687336246984942

e**123.45 =

410822093109668964239148908443317876138879647013995774.2951431466270782257597573259486687336246984942

"""

====================================================

This is also on the web at <http://python.pastebin.com/f1a5b9e03>.

Examples of problem x's: 10, 20, 30, 40, 100, 101

Examples of OK x's: 0.2, 5, 10.1, 11, 33.3, 123.45

Thanks,

Dick Moores

Message has been deleted

Nov 20, 2007, 4:45:02 AM11/20/07

to pytho...@python.org

At 12:45 AM 11/20/2007, Dennis Lee Bieber wrote:

>On Mon, 19 Nov 2007 23:41:02 -0800, Dick Moores <r...@rcblue.com>

>declaimed the following in comp.lang.python:

>

> > a way to get it to break where I want it to, i.e., when the sum

> > equals the limit as closely as the precision allows?

> >

>

> > if sum >= limit:

>

> Well, since it ISN'T a case of testing for an absolute equivalence

>with floats...

>

> Perhaps putting a "print sum, limit" before that point would reveal

>what type of values you are encountering.

>On Mon, 19 Nov 2007 23:41:02 -0800, Dick Moores <r...@rcblue.com>

>declaimed the following in comp.lang.python:

>

> > a way to get it to break where I want it to, i.e., when the sum

> > equals the limit as closely as the precision allows?

> >

>

>

> Well, since it ISN'T a case of testing for an absolute equivalence

>with floats...

>

> Perhaps putting a "print sum, limit" before that point would reveal

>what type of values you are encountering.

If you run the program you'll see exactly that, if I understand you

correctly. <http://python.pastebin.com/f2f06fd76> shows the full

output for a precision of 50 and x == 5.

Dick

Nov 20, 2007, 6:53:50 AM11/20/07

to Dick Moores, pytho...@python.org

On Nov 20, 2007 8:41 AM, Dick Moores <r...@rcblue.com> wrote:

> I'm writing a demo of the infinite series

>

> x**0/0! + x**1/1! + x**2/2! + x**3/3! + ... = e**x (x is non-negative)

>

> It works OK for many x, but for many the loop doesn't break. Is there

> I'm writing a demo of the infinite series

>

> x**0/0! + x**1/1! + x**2/2! + x**3/3! + ... = e**x (x is non-negative)

>

> It works OK for many x, but for many the loop doesn't break. Is there

> a way to get it to break where I want it to, i.e., when the sum

> equals the limit as closely as the precision allows?

>

> equals the limit as closely as the precision allows?

>

Hi,

Checking that sum >= e**x will generally not work, because e**x might

have been rounded up while the sum might repeatedly be rounding down.

If this happens, no matter how many terms you add, the sum will never

reach the limit (one of the curses of finite-precision arithmetic).

One solution is to use directed rounding. First compute the limit with

downward rounding:

mpf.round_down()

limit = e**x

mpf.round_default()

Then compute every term in the sum with upward rounding:

mpf.round_down()

fac = factorial(n)

mpf.round_up()

term = x**n / fac

sum += term

(Note that the factorial should be rounded down to obtain upward

rounding in the term, since you're taking its reciprocal.)

This should guarantee that the sum eventually exceeds the limit.

As a simpler, less rigorous alternative, instead of checking if sum >=

limit, check (for example) whether

abs(sum - limit) / limit <= mpf(10)**(-precision+3)

i.e., if the sum is within 3 digits of the limit. This is the usual

way to test for numerical equality of floating-point numbers.

Fredrik

Nov 20, 2007, 1:35:45 PM11/20/07

to pytho...@python.org

I tried out both ways, and found that the second one best suited my

purposes. Please see the 2 highlighted lines in

<http://python.pastebin.com/fcc23b10>

Note that to break the loop I found that this does the job:

if abs(sum - limit) / limit <= mpf(10)**(-precision+1): # I

changed your +3 to +1

break

Fredrik, thanks VERY much for your terrific instruction!

Dick

Nov 20, 2007, 1:42:48 PM11/20/07

to

Instead of comparing sum to the "known" value of e**x, why not test

for convergence? I.e., if sum == last_sum: break. Seems like that

would be more robust (you don't need to know the answer to computer

the answer), since it seems like it should converge.

for convergence? I.e., if sum == last_sum: break. Seems like that

would be more robust (you don't need to know the answer to computer

the answer), since it seems like it should converge.

--Nathan Davis

Nov 20, 2007, 4:00:04 PM11/20/07

to pytho...@python.org

At 10:42 AM 11/20/2007, davis...@gmail.com wrote:

>Instead of comparing sum to the "known" value of e**x, why not test

>for convergence? I.e., if sum == last_sum: break. Seems like that

>would be more robust (you don't need to know the answer to computer

>the answer), since it seems like it should converge.

>Instead of comparing sum to the "known" value of e**x, why not test

>for convergence? I.e., if sum == last_sum: break. Seems like that

>would be more robust (you don't need to know the answer to computer

>the answer), since it seems like it should converge.

Yes! And believe it or not I did that before seeing your post. Works

well. See <http://python.pastebin.com/f7c37186a>

And also with the amazing Chudnovsky algorithm for pi. See

<http://python.pastebin.com/f4410f3dc>

Thanks,

Dick

Nov 20, 2007, 4:32:55 PM11/20/07

to Dick Moores, pytho...@python.org

On Nov 20, 2007 10:00 PM, Dick Moores <r...@rcblue.com> wrote:

> And also with the amazing Chudnovsky algorithm for pi. See

> <http://python.pastebin.com/f4410f3dc>

> And also with the amazing Chudnovsky algorithm for pi. See

> <http://python.pastebin.com/f4410f3dc>

Nice! I'd like to suggest two improvements for speed.

First, the Chudnovsky algorithm uses lots of factorials, and it's

rather inefficient to call mpmath's factorial function from scratch

each time. You could instead write a custom factorial function that

only uses multiplications and caches results, something like this:

cached_factorials = [mpf(1)]

def f(n):

n = int(n)

if n < len(cached_factorials):

return cached_factorials[n]

p = cached_factorials[-1]

for i in range(len(cached_factorials), n+1):

p *= i

cached_factorials.append(p)

return p

(In some future version of mpmath, the factorial function might be

optimized so that you won't have to do this.)

Second, to avoid unnecessary work, factor out the fractional power of

640320 that occurs in each term. That is, change the "denom =" line to

denom = (f(3*k) * ((f(k))**3) * (640320**(3*k)))

and then multiply it back in at the end:

print 1/(12*sum/640320**(mpf(3)/2))

With these changes, the time to compute 1,000 digits drops to only 0.05 seconds!

Further improvements are possible.

Fredrik

Nov 20, 2007, 5:23:34 PM11/20/07

to

On Tue, 20 Nov 2007 10:42:48 -0800, davis...@gmail.com wrote:

> Instead of comparing sum to the "known" value of e**x, why not test for

> convergence? I.e., if sum == last_sum: break. Seems like that would be

> more robust (you don't need to know the answer to computer the answer),

> since it seems like it should converge.

That will only work if you know that the terms in your sequence are

monotonically decreasing: that is, if each term is smaller than the last.

It also assumes that the terms decrease reasonably rapidly, and you want

the full precision available in floats. Are you sure your algorithm is

that precise? It's one thing to produce 16 decimal places in your result,

but if only the first 12 of them are meaningful, and the last four are

inaccurate, you might as well not bother with the extra work.

--

Steven.

Nov 25, 2007, 3:00:42 AM11/25/07

to pytho...@python.org

Fredrik,

I'm afraid I'm just too dumb to see how to use your first suggestion

of cached_factorials. Where do I put it and def()? Could you show me,

even on-line, what to do? <http://py77.python.pastebin.com/f48e4151c>

You (or anyone) can submit an amendment to my code using the textbox.

I did make the denom change, and see that it does improve the speed a bit.

Thanks,

Dick

Nov 25, 2007, 5:55:38 AM11/25/07

to

(1) Replace line 8 (from mpmath import factorial as f) with Fredrik's

code. (2) Test it to see that it gave the same answers as before ...

[time passes]

Wow! [w/o psyco] Pi to 1000 decimal places takes 13 seconds with

original code and 0.2 seconds with Fredrik's suggestion. 2000: 99

seconds -> 0.5 seconds.

Nov 25, 2007, 6:26:05 AM11/25/07

to Dick Moores, pytho...@python.org

On Nov 25, 2007 9:00 AM, Dick Moores <r...@rcblue.com> wrote:

> Fredrik,

>

> I'm afraid I'm just too dumb to see how to use your first suggestion

> of cached_factorials. Where do I put it and def()? Could you show me,

> even on-line, what to do? <http://py77.python.pastebin.com/f48e4151c>

> Fredrik,

>

> I'm afraid I'm just too dumb to see how to use your first suggestion

> of cached_factorials. Where do I put it and def()? Could you show me,

> even on-line, what to do? <http://py77.python.pastebin.com/f48e4151c>

> You (or anyone) can submit an amendment to my code using the textbox.

>

> I did make the denom change, and see that it does improve the speed a bit.

>

> I did make the denom change, and see that it does improve the speed a bit.

I edited the pastebin code, see: http://py77.python.pastebin.com/m6b2b34b7

Fredrik

Nov 25, 2007, 8:47:06 AM11/25/07

to pytho...@python.org

Wow. your f() is ingenious, Frederik. Thanks very much.

Any more tricks up your sleeve? You did say a post or so ago,

"Further improvements are possible."

Dick

Nov 25, 2007, 1:32:34 PM11/25/07

to Dick Moores, pytho...@python.org

On Nov 25, 2007 2:47 PM, Dick Moores <r...@rcblue.com> wrote:

> Wow. your f() is ingenious, Frederik. Thanks very much.

>

> Any more tricks up your sleeve? You did say a post or so ago,

> "Further improvements are possible."

> Wow. your f() is ingenious, Frederik. Thanks very much.

>

> Any more tricks up your sleeve? You did say a post or so ago,

> "Further improvements are possible."

The next improvement would be to find a recurrence formula for the

terms instead of computing them directly. So for example, if summing

over c[n] = x**n / n!, instead of computing c[n] = x**n / factorial(n)

for each n, you'd compute c[0] and then just do c[n] = c[n-1] * x / n

for each of the remaining terms.

Fredrik

Reply all

Reply to author

Forward

0 new messages

Search

Clear search

Close search

Google apps

Main menu