At 03:53 AM 11/20/2007, Fredrik Johansson wrote:
>On Nov 20, 2007 8:41 AM, Dick Moores <[EMAIL PROTECTED]> 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
>
>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.
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
--
http://mail.python.org/mailman/listinfo/python-list
