comparing to "zero" is problematic
try changing the first few line to
aa=abs(a)
ab=abs(b)
ac=abs(c)
tol=max((aa,ab,ac))/1099511627776 #/2**40
if (aa<tol and ab<tol): # Case for handling Liner Equation
return np.array([(-d * 1.0) / c]) # Returning linear root as numpy array.
elif (aa < tol):
for a=-5.32282969291e-15 b=0.43709356404 c=0.819917378017 d=-0.199301067523
goes from
complex_roots: [ 8.211676669305905e+13 +0.j
-8.867187500000000e-01+353649.2408007157j
-8.867187500000000e-01-353649.2408007157j]
to
complex_roots: [ 0.2177888401916245 -2.093628352534504 ]
when numpy.roots says
complex_roots: [ 8.2116766693059141e+13 -2.0936283525344557e+00 2.1778884019162431e-01]
this came up in the context of a highly specialized pchip spline implementation I was working on