Python Decimal type precision error

Question

I have a quite puzzling issue that I suspect has to do with scientific notation and decimal precision. Here is part of my code:

    def atan(x):
        # Calculate arctan(1/x)
        x = Decimal(x)
        current_value = Decimal(0)
        divisor = 1
        x_squared = x * x
        current_term = 1 / x


        while True:
            current_value += current_term

            divisor += 2
            current_term = (-current_term / x_squared) / divisor
            print(current_term)

            # The issue
            if current_term == Decimal(0):
                break

        return current_value

    print(atan(5))

This is based on the formula atan(1/x) = 1/x - 1/(3x^3) + 1/(5x^5) - ...

However, I discovered that current_term, which gets smaller every loop iteration, is going into values like 4E-80000. Since I've set my decimal precision getcontext().prec to 20, current term should not even support these values. I think somehow current_term is not of decimal type but of scientific notation/float type, but python tells me it's still decimal type.

The correct value for arctan(1/5) is about 0.1973955. I get a value of 0.1973545, which is wrong starting at the 5th digit. Even if I manually break the loop the value is still wrong for some reason. Any help fixing this issue is appreciated.

I think your algorithm itself is wrong - you never change x_squared in your loop. — Tony Suffolk 66
– Tony Suffolk 66, Commented Jul 26, 2014 at 5:12
@TonySuffolk66 x_squared is supposed to be constant; current_term is divided by x_squared every loop — qwr
– qwr, Commented Jul 26, 2014 at 5:16
Why do you think 4E-80000 is not a valid value for a precision of 20? It has just one significant digit — the 4. Precision doesn't determine the total number of digits after the decimal point. If it would then Decimal wouldn't be floating point values but fixed point values. — BlackJack
– BlackJack, Commented Jul 26, 2014 at 5:26
@BlackJack That is a very good point I didn't think of. Even then, though, 4E-80000==Decimal(0) should return True (I tested it) — qwr
– qwr, Commented Jul 26, 2014 at 5:28

Raymond Hettinger · Accepted Answer · 2014-07-26 05:34:08Z

3

You code doesn't match the formula. It got a little too tricky with inferring one term from the next ;-) The 1/(5x^5) term isn't a multiple of 1/(3x^3) term.

Here is code that models the formula directly:

from decimal import Decimal

def atan_recip(x):
    # Calculate arctan(1/x)
    x = Decimal(x)

    total = Decimal(0)
    sign = 1
    for i in range(1, 35, 2):
        total += sign / (i * x ** i)
        sign = -sign
        print(total)

atan_recip(5)

The output is what you expected:

0.2
0.1973333333333333333333333333
0.1973973333333333333333333333
0.1973955047619047619047619047
0.1973955616507936507936507936
0.1973955597889754689754689754
0.1973955598519908535908535908
0.1973955598498063202575202575
0.1973955598498834214339908457
0.1973955598498806620234645299
0.1973955598498807618878454823
0.1973955598498807582406246127
0.1973955598498807583748423407
0.1973955598498807583698713137
0.1973955598498807583700564416
0.1973955598498807583700495142
0.1973955598498807583700497745

edited Jul 26, 2014 at 5:34

answered Jul 26, 2014 at 5:29

Raymond Hettinger

229k67 gold badges405 silver badges504 bronze badges

Sign up to request clarification or add additional context in comments.

5 Comments

qwr Over a year ago

Your code seems to work, but I don't see where my code doesn't match the formula. Also, I still can't figure out why current_term==Decimal(0) doesn't evaluate to True.

Raymond Hettinger Over a year ago

The issue with your code doesn't "has to do with scientific notation and decimal precision". The algebra is a bit off. The ratio of 5x^5 to 3x^3 isn't x^2, it is 5/3 x^2.

qwr Over a year ago

Ah, I see now why my formula was incorrect. However, shouldn't current_term==Decimal(0) evaluate to True in the loop?

Raymond Hettinger Over a year ago

Nope, the term gets smaller and smaller but doesn't hit zero (otherwise it wouldn't be an approximation, it would be exact). The decimal value keeps 28 places of precision, not 28 places to the right of the decimal point (notice the exponents in the your print out of the current_term.

qwr Over a year ago

Thank you very much for your help, I had some confusion over basic algebra and fixed point vs floating point.

Collectives™ on Stack Overflow

Python Decimal type precision error

1 Answer 1

5 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

5 Comments

Your Answer

Sign up or log in

Post as a guest

Related