Python Forum

Hello forum. I am currently writing a program to determine the roots of a cubic equation. However I have come across some trouble with my script.

import math
def sols(a, b, c, d ):
    disc = 18*a*b*c*d - 4*(b**3)*d + (b**2)*(c**2) - 4*a*(c**3) - 27*(a**2)*(d**2)
    x = b**2.0 - 3.0*a*c
    y = 2.0*(b**3.0) - 9.0*a*b*c + 27.0*(a**2.0)*d
    alpha = complex(-0.5, 0.5*math.sqrt(3))
    if disc >= 0:
        C = ((y + complex (0, math.sqrt (27.*(a**2)*disc)))/2)**(1.0/3)
        sol1 = -(1/3*a)*(b + C + x/C)
        sol2 = -(1/3*a)*(b + C*alpha + x/(C*alpha))
        sol3 = -(1/3*a)*(b + C*(alpha**2)) + x/(C*(alpha**2))

This is my program so far. The problem I have is that,  for the values a = 2, b = 3, c= -11, d = -6 :

In[11]: C

Out[12]: (7.5+4.330127018922192j)


In[13]: x/C

Out[14]: (7.500000000000002-4.330127018922193j)

In order to get sol1 I must add C and x/C. I want this to equal (15.00000000000002 + 0j)

Instead I get:

In[26]: C + x/C

Out[27]: (15.000000000000002-8.881784197001252e-16j)

I dont understand why I have gotten this result. I believe the complex parts should cancel each other. I am aware the last decimal point of C and x/C differ and i believe this is the problem. However the two complex parts should cancel, so I am not sure why I have this problem of the differing decimal value.

Anyone have any idea? 

Many thanks.

Also if you cannot already tell, I am very much a python newbie. Apologies if I have not made myself clear.

It's a floating point error. Note that the final complex part ends with e-16. That means it's time 10^-16. So the complex part isn't -8.88, it's -0.0000000000000000888. Floating point doesn't exactly capture all real numbers, and this can lead to really small error like the one you are seeing. When using floating point, you generally don't test for equality of two numbers, you test for the difference between the two numbers being less than some very small value.

You could format the numbers for better output:

#python3
c = 15.000000000000002-8.881784197001252e-16j

print(c)
print("{0:.16f}".format(c))

==>
(15.000000000000002-8.881784197001252e-16j)
15.0000000000000018-0.0000000000000009j