Velocity Reviews - Computer Hardware Reviews

Velocity Reviews > Newsgroups > Programming > Python > Trouble getting loop to break

Reply
Thread Tools

Trouble getting loop to break

 
 
Dick Moores
Guest
Posts: n/a
 
      11-20-2007
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

 
Reply With Quote
 
 
 
 
Dennis Lee Bieber
Guest
Posts: n/a
 
      11-20-2007
On Mon, 19 Nov 2007 23:41:02 -0800, Dick Moores <(E-Mail Removed)>
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.
--
Wulfraed Dennis Lee Bieber KD6MOG
http://www.velocityreviews.com/forums/(E-Mail Removed) (E-Mail Removed)
HTTP://wlfraed.home.netcom.com/
(Bestiaria Support Staff: (E-Mail Removed))
HTTP://www.bestiaria.com/
 
Reply With Quote
 
 
 
 
Dick Moores
Guest
Posts: n/a
 
      11-20-2007
At 12:45 AM 11/20/2007, Dennis Lee Bieber wrote:
>On Mon, 19 Nov 2007 23:41:02 -0800, Dick Moores <(E-Mail Removed)>
>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/f2f06fd76> shows the full
output for a precision of 50 and x == 5.

Dick


 
Reply With Quote
 
davisn90210@gmail.com
Guest
Posts: n/a
 
      11-20-2007
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 <(E-Mail Removed)> 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
> 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


 
Reply With Quote
 
Dick Moores
Guest
Posts: n/a
 
      11-20-2007
At 10:42 AM 11/20/2007, (E-Mail Removed) 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


 
Reply With Quote
 
Fredrik Johansson
Guest
Posts: n/a
 
      11-20-2007
On Nov 20, 2007 10:00 PM, Dick Moores <(E-Mail Removed)> wrote:
> 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
 
Reply With Quote
 
Steven D'Aprano
Guest
Posts: n/a
 
      11-20-2007
On Tue, 20 Nov 2007 10:42:48 -0800, (E-Mail Removed) 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.
 
Reply With Quote
 
Dick Moores
Guest
Posts: n/a
 
      11-25-2007
At 01:32 PM 11/20/2007, Fredrik Johansson wrote:
>On Nov 20, 2007 10:00 PM, Dick Moores <(E-Mail Removed)> wrote:
> > 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




 
Reply With Quote
 
John Machin
Guest
Posts: n/a
 
      11-25-2007
On Nov 25, 7:00 pm, Dick Moores <(E-Mail Removed)> wrote:
> At 01:32 PM 11/20/2007, Fredrik Johansson wrote:
>
>
>
>
>
>
>
> >On Nov 20, 2007 10:00 PM, Dick Moores <(E-Mail Removed)> wrote:
> > > 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.
 
Reply With Quote
 
Fredrik Johansson
Guest
Posts: n/a
 
      11-25-2007
On Nov 25, 2007 9:00 AM, Dick Moores <(E-Mail Removed)> 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>
> 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
 
Reply With Quote
 
 
 
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are Off


Similar Threads
Thread Thread Starter Forum Replies Last Post
Triple nested loop python (While loop insde of for loop inside ofwhile loop) Isaac Won Python 9 03-04-2013 10:08 AM
`if (!p ? i++ : 0) break;' == `if (!p){ i++; break;}' ? lovecreatesbea...@gmail.com C Programming 12 04-14-2008 07:59 AM
Trouble getting past the loop. mmoski Java 5 02-03-2007 04:30 AM
How to break a while loop inside a switch statement? howachen@gmail.com Java 16 07-11-2006 10:23 AM
Getting a loop to activate a loop above it Byte Python 4 03-24-2006 03:04 AM



Advertisments