Trouble getting loop to break

Discussion in 'Python' started by Dick Moores, Nov 20, 2007.

  1. Dick Moores

    Dick Moores Guest

    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
    Dick Moores, Nov 20, 2007
    #1
    1. Advertising

  2. 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 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://wlfraed.home.netcom.com/
    (Bestiaria Support Staff: )
    HTTP://www.bestiaria.com/
    Dennis Lee Bieber, Nov 20, 2007
    #2
    1. Advertising

  3. Dick Moores

    Dick Moores Guest

    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 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
    Dick Moores, Nov 20, 2007
    #3
  4. Dick Moores

    Guest

    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 <> 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 =
    > 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
    , Nov 20, 2007
    #4
  5. Dick Moores

    Dick Moores Guest

    At 10:42 AM 11/20/2007, 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
    Dick Moores, Nov 20, 2007
    #5
  6. On Nov 20, 2007 10:00 PM, Dick Moores <> 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 Johansson, Nov 20, 2007
    #6
  7. On Tue, 20 Nov 2007 10:42:48 -0800, 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.
    Steven D'Aprano, Nov 20, 2007
    #7
  8. Dick Moores

    Dick Moores Guest

    At 01:32 PM 11/20/2007, Fredrik Johansson wrote:
    >On Nov 20, 2007 10:00 PM, Dick Moores <> 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
    Dick Moores, Nov 25, 2007
    #8
  9. Dick Moores

    John Machin Guest

    On Nov 25, 7:00 pm, Dick Moores <> wrote:
    > At 01:32 PM 11/20/2007, Fredrik Johansson wrote:
    >
    >
    >
    >
    >
    >
    >
    > >On Nov 20, 2007 10:00 PM, Dick Moores <> 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.
    John Machin, Nov 25, 2007
    #9
  10. On Nov 25, 2007 9:00 AM, Dick Moores <> 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
    Fredrik Johansson, Nov 25, 2007
    #10
  11. Dick Moores

    Dick Moores Guest

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


    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
    Dick Moores, Nov 25, 2007
    #11
  12. On Nov 25, 2007 2:47 PM, Dick Moores <> 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."


    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
    Fredrik Johansson, Nov 25, 2007
    #12
    1. Advertising

Want to reply to this thread or ask your own question?

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.
Similar Threads
  1. Replies:
    16
    Views:
    3,574
    steve
    Jul 11, 2006
  2. Byte
    Replies:
    4
    Views:
    413
  3. mmoski
    Replies:
    5
    Views:
    328
    Luc The Perverse
    Feb 3, 2007
  4. Replies:
    12
    Views:
    954
  5. Isaac Won
    Replies:
    9
    Views:
    364
    Ulrich Eckhardt
    Mar 4, 2013
Loading...

Share This Page