473,249 Members | 1,805 Online
Bytes | Software Development & Data Engineering Community
Post Job

Home Posts Topics Members FAQ

Join Bytes to post your question to a community of 473,249 software developers and data experts.

Trouble getting loop to break

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 =
41082209310966896423914890844331787613887964701399 5774.295143146627078225759757325948668733624698494 2
k = 427
exp(123.45) =
41082209310966896423914890844331787613887964701399 5774.295143146627078225759757325948668733624698494 2
e**123.45 =
41082209310966896423914890844331787613887964701399 5774.295143146627078225759757325948668733624698494 2
"""
================================================== ==

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

Nov 20 '07 #1
10 1154
At 12:45 AM 11/20/2007, Dennis Lee Bieber wrote:
>On Mon, 19 Nov 2007 23:41:02 -0800, Dick Moores <rd*@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.
If you run the program you'll see exactly that, if I understand you
correctly. <http://python.pastebin.com/f2f06fd76shows the full
output for a precision of 50 and x == 5.

Dick
Nov 20 '07 #2
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.

--Nathan Davis

On Nov 20, 1:41 am, Dick Moores <r...@rcblue.comwrote:
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 =
41082209310966896423914890844331787613887964701399 5774.295143146627078225759757325948668733624698494 2
k = 427
exp(123.45) =
41082209310966896423914890844331787613887964701399 5774.295143146627078225759757325948668733624698494 2
e**123.45 =
41082209310966896423914890844331787613887964701399 5774.295143146627078225759757325948668733624698494 2
"""
================================================== ==

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
Nov 20 '07 #3
At 10:42 AM 11/20/2007, da*********@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.
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 '07 #4
On Nov 20, 2007 10:00 PM, Dick Moores <rd*@rcblue.comwrote:
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 '07 #5
On Tue, 20 Nov 2007 10:42:48 -0800, da*********@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 20 '07 #6
At 01:32 PM 11/20/2007, Fredrik Johansson wrote:
>On Nov 20, 2007 10:00 PM, Dick Moores <rd*@rcblue.comwrote:
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
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 '07 #7
On Nov 25, 7:00 pm, Dick Moores <r...@rcblue.comwrote:
At 01:32 PM 11/20/2007, Fredrik Johansson wrote:


On Nov 20, 2007 10:00 PM, Dick Moores <r...@rcblue.comwrote:
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

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>
(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 '07 #8
On Nov 25, 2007 9:00 AM, Dick Moores <rd*@rcblue.comwrote:
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 edited the pastebin code, see: http://py77.python.pastebin.com/m6b2b34b7

Fredrik
Nov 25 '07 #9
At 03:26 AM 11/25/2007, Fredrik Johansson wrote:
>On Nov 25, 2007 9:00 AM, Dick Moores <rd*@rcblue.comwrote:
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 edited the pastebin code, see: http://py77.python.pastebin.com/m6b2b34b7

Fredrik
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 '07 #10
On Nov 25, 2007 2:47 PM, Dick Moores <rd*@rcblue.comwrote:
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
Nov 25 '07 #11

This thread has been closed and replies have been disabled. Please start a new discussion.

Similar topics

14
by: Ankit Aneja | last post by:
The code of classes given below is for server to which clients connect i want to get ip address of client which has connected pls help how can i get //listen class public class listen {
11
by: Snowbleach | last post by:
Hello, I am having trouble with the following console application. I am using delegates/event to make a clock raise an alarm when a certain amount of time has been reached. The following compiles,...
3
by: bmerlover | last post by:
I believe my problem lies inside the while loop. When I click the play button on the gui app, it goes inside the while loop, reads the file and calls the necessary function to do what it needs to do....
0
by: bmerlover | last post by:
This code makes sense to me, I'm just having trouble trying to understand why it doesn't work correctly. This is a GUI APP. When the Play button is Clicked, the play_Click(System::Object * sender,...
0
by: jthep | last post by:
Hi, I'm trying to get user input for a record book but I'm having trouble as I think I'm not making the getline function read the buffer correctly. I have the following variables declared in...
2
by: AceX | last post by:
Hi, I'm trying to use PHP to to a lot of repetetive work for me on a site. I've gotten most of it working fine - EXCEPT: the PHP is generating a table. It's looping through a CSV file until it...
9
Steel546
by: Steel546 | last post by:
This program is used to calculate GPA. I'm having trouble actually getting an output. Alright, I KNOW that I don't have an output statement, but I don't know where to put it... heh. Netbeans keeps...
5
matheussousuke
by: matheussousuke | last post by:
Hi guys, good morning. I've just get this script for converting mysql tables from wordpress, and I want to use it in my server, but no with wordpress, with oscommerce, a friend of mine told me a...
5
matheussousuke
by: matheussousuke | last post by:
Hello, I'm using tiny MCE plugin on my oscommerce and it is inserting my website URL when I use insert image function in the emails. The goal is: Make it send the email with the URL...
2
isladogs
by: isladogs | last post by:
The next Access Europe meeting will be on Wednesday 7 Feb 2024 starting at 18:00 UK time (6PM UTC) and finishing at about 19:30 (7.30PM). In this month's session, the creator of the excellent VBE...
0
by: MeoLessi9 | last post by:
I have VirtualBox installed on Windows 11 and now I would like to install Kali on a virtual machine. However, on the official website, I see two options: "Installer images" and "Virtual machines"....
0
by: DolphinDB | last post by:
The formulas of 101 quantitative trading alphas used by WorldQuant were presented in the paper 101 Formulaic Alphas. However, some formulas are complex, leading to challenges in calculation. Take...
0
by: Aftab Ahmad | last post by:
Hello Experts! I have written a code in MS Access for a cmd called "WhatsApp Message" to open WhatsApp using that very code but the problem is that it gives a popup message everytime I clicked on...
0
by: ryjfgjl | last post by:
ExcelToDatabase: batch import excel into database automatically...
0
isladogs
by: isladogs | last post by:
The next Access Europe meeting will be on Wednesday 6 Mar 2024 starting at 18:00 UK time (6PM UTC) and finishing at about 19:15 (7.15PM). In this month's session, we are pleased to welcome back...
0
by: marcoviolo | last post by:
Dear all, I would like to implement on my worksheet an vlookup dynamic , that consider a change of pivot excel via win32com, from an external excel (without open it) and save the new file into a...
0
isladogs
by: isladogs | last post by:
The next Access Europe meeting will be on Wednesday 6 Mar 2024 starting at 18:00 UK time (6PM UTC) and finishing at about 19:15 (7.15PM). In this month's session, we are pleased to welcome back...
0
by: Vimpel783 | last post by:
Hello! Guys, I found this code on the Internet, but I need to modify it a little. It works well, the problem is this: Data is sent from only one cell, in this case B5, but it is necessary that data...

By using Bytes.com and it's services, you agree to our Privacy Policy and Terms of Use.

To disable or enable advertisements and analytics tracking please visit the manage ads & tracking page.