Trouble getting loop to break

[Fredrik Johansson]

On Nov 20, 2007 8:41 AM, Dick Moores <rdm at rcblue.com> 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.

Fredrik