440,235 Members | 1,008 Online Need help? Post your question and get tips & solutions from a community of 440,235 IT Pros & Developers. It's quick & easy.

# Trouble getting loop to break

 P: n/a 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 . 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 Replies

 P: n/a At 12:45 AM 11/20/2007, Dennis Lee Bieber wrote: >On Mon, 19 Nov 2007 23:41:02 -0800, Dick Moores 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 equivalencewith floats... Perhaps putting a "print sum, limit" before that point would revealwhat type of values you are encountering. If you run the program you'll see exactly that, if I understand you correctly.

 P: n/a 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 = 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 . 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

 P: n/a At 10:42 AM 11/20/2007, da*********@gmail.com wrote: >Instead of comparing sum to the "known" value of e**x, why not testfor convergence? I.e., if sum == last_sum: break. Seems like thatwould be more robust (you don't need to know the answer to computerthe answer), since it seems like it should converge. Yes! And believe it or not I did that before seeing your post. Works well. See And also with the amazing Chudnovsky algorithm for pi. See Thanks, Dick Nov 20 '07 #4

 P: n/a On Nov 20, 2007 10:00 PM, Dick Moores 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

 P: n/a 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

 P: n/a At 01:32 PM 11/20/2007, Fredrik Johansson wrote: >On Nov 20, 2007 10:00 PM, Dick Moores Nice! I'd like to suggest two improvements for speed.First, the Chudnovsky algorithm uses lots of factorials, and it'srather inefficient to call mpmath's factorial function from scratcheach time. You could instead write a custom factorial function thatonly 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 beoptimized so that you won't have to do this.)Second, to avoid unnecessary work, factor out the fractional power of640320 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 only0.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? 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

 P: n/a On Nov 25, 7:00 pm, Dick Moores 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? (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

 P: n/a On Nov 25, 2007 9:00 AM, Dick Moores 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

 P: n/a At 03:26 AM 11/25/2007, Fredrik Johansson wrote: >On Nov 25, 2007 9:00 AM, Dick Moores 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/m6b2b34b7Fredrik 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

 P: n/a On Nov 25, 2007 2:47 PM, Dick Moores

### This discussion thread is closed

Replies have been disabled for this discussion. 