Bug in Decimal??

[Mark Dickinson]

<pleasedontspam <at> isp.com> writes:

> I've tested on all platforms I know of and confirmed it. The wrong digit
> occurs in the middle of the number.  Propagation error would have a bad digit
> near the end, and garbage after that. Here there's a perfect sequence of
> numbers, but with one single digit changed in the middle of the number.  No
> error propagation in a series expansion can do that.

I can see how it might be surprising if you don't think about it too hard, but
I'm afraid that you're wrong here: error propagation is *exactly* what's
causing the effects you're seeing.

Here's another way of looking at it: if you truncate the Taylor series about 0
for (1 + x) / (1 - x) to k (>= 1) terms, you get the polynomial (1 + x - 2x^k)
/ (1 - x).  For example, taking k to be 3, we're getting (1 + x - 2x^3) / (1 -
x).  Given that the particular value of x you're testing with has the form
10**<negative>, rounding your intermediate result to the working precision has
exactly the effect of truncating the series at some k.

Now you can compute and compare (by hand, via Wolfram alpha, or however you
like) the Taylor series expansions for log((1 + x) / (1 - x)) and log((1 + x -
2x^3) / (1 - x)).  For the first you'll see:

  2x + 2/3 x^3  + 2/5 x^5 - 2/7 x^7 + 2/9 x^9 - ...

and for the second you'll get:

  2x - 4/3 x^3 + 2 x^4 - 8/5 x^5 + 16/7 x^7 - ...

The difference between the two series is:

  -2x^3 + 2x^4 - 2x^5 + 2x^7 - 4x^8 + ...

So again with x a small power of 10, you're going to see a single-digit error
from the -2x^3 term, and another single-digit error further along from
the 2x^3 term, and so on.

Here's a simpler example of the same phenomenon.  Note how the error propagation
leads to a single incorrect digit in the *middle* of the digit string.

Python 3.4.0 (default, Mar 25 2014, 11:07:05) 
[GCC 4.2.1 Compatible Apple LLVM 5.1 (clang-503.0.38)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> from decimal import Decimal
>>> x = Decimal('1e-15')
>>> y = (1 - 2 * x) / (1 - x)
>>> 2 * x + (y - 1) * (1 - x)  # Mathematically, expect to get 'x' back.
Decimal('1.000000000000001000000000000E-15')
>>> x
Decimal('1E-15')

-- 
Mark