How can I get more exact decimal values in Python?

Question

from math import sqrt a=1e-8 b=10 c=1e-8 x1 = ((-b)-sqrt((b**2)-(4*a*c)))/(2*a) x2 = ((-b)+sqrt((b**2)-(4*a*c)))/(2*a) print 'x1 = {}'.format(x1) print 'x2 = {}'.format(x2) print (4*a*c) print (sqrt(b**2-4*a*c)) print b**2 print 2*a When I run the program, this returns: x1 = -1e+09 x2 = 0.0 4e-16 10.0 100.0 2e-08 What I need is for x2 to be equal to -1e-9. The problem seems to be with the sqrt((b**2)-(4*a*c)) as it gives 10 as a result, obviously because 4*(10^-8)*(10^-8) is almost equal to 0, and is considered 0 by python. This results in: sqrt((b**2)-(4*a*c)) = sqrt(b**2) = sqrt(10**2) = 10 Any help would be greatly appreciated

Answer 1

Given:

a = 1e-8

b = 10

c = 1e-8

x1 = ((-b) - sqrt(b**2 - 4*a*c)) / (2*a) (1st formula)

x2 = ((-b) + sqrt(b**2 - 4*a*c)) / (2*a) (2nd formula)

The problem is not how many decimals Python is printing. The problem is that the 2nd formula subtracts two numbers that are almost equal, and when that happens almost every correct digit gets cancelled out. This is called catastrophic cancellation.

Lets follow what happens inside the 2nd formula:

b**2 = 100

4*a*c = 4e-16

b**2 - 4*a*c = 100 - 4e-16 Subtract

= 100.0 A float only carries about 16 significant digits, so the 4e-16 falls off the end

sqrt(100.0) = 10.0 Take the square root

(-b) + 10.0 = -10 + 10.0 Plug into the top of the 2nd formula

= 0.0 Every correct digit is now gone

_______________________________

x2 = 0.0 / 2e-08 = 0.0

Notice the 1st formula never has this trouble, because (-b) and -sqrt(...) carry the same sign, so nothing cancels:

x1 = (-10 - 10.0) / 2e-08

= -20 / 2e-08 Simplify the top

= -1e+09 This value is already correct

So the fix is to compute the root that does not cancel, then get the other one from the fact that for a*x**2 + b*x + c = 0 the two roots multiply to c/a.

Firstly, we can write:

x1 * x2 = c/a Product of the roots

x2 = c / (a * x1) Divide by a*x1 on both sides to solve for x2

Now lets plug in the numbers:

x2 = 1e-8 / (1e-8 * -1e+09) Plug in a = 1e-8, c = 1e-8, x1 = -1e+09

= 1e-8 / -10.0 Simplify the bottom

= -1e-09

In code:

from math import sqrt

a = 1e-8

b = 10.0

c = 1e-8

d = sqrt(b*b - 4*a*c)

if b >= 0:

x1 = ((-b) - d) / (2*a)

else:

x1 = ((-b) + d) / (2*a)

x2 = c / (a * x1)

print 'x1 = {}'.format(x1)

print 'x2 = {}'.format(x2)

which gives:

x1 = -1000000000.0

x2 = -1e-09

The if statement is there so that the square root is always added to (-b) with the same sign, no matter whether b comes in positive or negative. That way the cancelling never happens on either root.

If you ever do need genuinely more digits and not just the correct answer, you can switch to the decimal module and raise the precision:

from decimal import Decimal, getcontext

getcontext().prec = 50

a = Decimal('1e-8')

but that is slower, and for this problem the rearrangement above is enough on its own.

x2 = -1e-09

Answer 2

When you are using floating point numbers in python, the best answer you can expect is double precision arithmetic, which means about 15 digits of precision.

In this case, the problem is that when you subtract that mall number from 100 you get a number very close to 100 and taking the square root gives a number very close to 10. However the difference from 10 is not displayed exactly correctly because of the limited precision. So yes, you have correctly identified the problem.

The best way to handle that situation is to detect when b² is much greater than 4ac, and to then use an approximation for -b + sqrt(b² -4ac) in those situations. (only when b is positive. If b is negative, then the addition is the problem)

When b is positive, and b² is much greater than 4ac, then -b + sqrt(b^2 - 4ac) is approximately equal to -2ac/b. You might check, for example that absolute value of (b^2/(4ac) is greater than about 10 when b is positive. If so, then use the approximation, otherwise calculate normally.

If you use the approximation I gave, you will get x = -c/b for those problem situations. That is exactly what you want.. The other root with the minus sign - b - sqrt (b² - 4ac) can be calculated normally when b is positive.

And for the cases where b is negative, then use the approximation with the negative square root.