3次方程式を解く

import numpy as np

def third_order_eq(a, b, c, d):

	term1 = (((-27*a**2*d + 9*a*b*c - 2*b**3)**2 + 4*(3*a*c - b**2)**3)**(1/2) - 27*a**2*d + 9*a*b*c - 2*b**3)**(1/3)
	term2 = 3*2**(1/3)*a
	term3 = 3*a*c - b**2
	term4 = b/(3*a)
	term5 = -1/(6*2**(1/3)*a)
	term8 = 3*2**(2/3)*a
	term6 = complex(1, -3**(1/2))
	term7 = complex(1, 3**(1/2))
	
	sol1 = term1/term2 - (2**(1/3)*term3)/(3*a*term1) - term4
	sol2 = term5*term6*term1 + (term7*term3)/(term8*term1) - term4
	sol3 = term5*term7*term1 + (term6*term3)/(term8*term1) - term4
	
	return sol1, sol2, sol3
	
if __name__ == '__main__':
	
	#3次方程式 係数
	coeff2 = [123, 456, 789, 101]

	x1, x2, x3 = np.roots(coeff2)
	print(x1)
	print(x2)
	print(x3)
	
	x1, x2, x3 = (third_order_eq(coeff2[0], coeff2[1], coeff2[2], coeff2[3]))
	print(x1)
	print(x2)
	print(x3)
	
	if x1.imag < 1e-10 and x1.imag > -1e-10:
		print('x1 = %.10f' % x1.real)
	if x2.imag < 1e-10 and x2.imag > -1e-10:
		print('x2 = %.10f' % x2.real)
	if x3.imag < 1e-10 and x3.imag > -1e-10:
		print('x3 = %.10f' % x3.real)