Dear all, I am not sure if I misunderstood the Polynomial.roots() method, but I get strange behavior. The roots are way outside of the expected values. That is not the case with np.roots(). What am I doing wrong? Code: import numpy as np def roots1(coef, x, y): for i in range(y.size): coefN = coef.copy() coefN[0] -= y[i] print(np.roots(coefN[::-1]) - x[i]) def roots2(coef, x, y): for i in range(y.size): print((np.polynomial.polynomial.Polynomial(coef) - y[i]).roots() - x[i]) x = np.arange(4) * 0.1 alpha, beta, gamma = 0.16, 1, 1e-15 coef = np.array([alpha, beta, gamma]) y = alpha + beta*x + gamma*x**2 print("roots") roots1(coef, x, y) print('Polynomial') roots2(coef, x, y) Outputs: roots [-1.e+15 0.e+00] [-1.e+15 0.e+00] [-1.00000000e+15 2.77555756e-17] [-1.00000000e+15 5.55111512e-17] Polynomial [-1.e+15 0.e+00] [-1.0e+15 2.5e-02] [-1.0e+15 -7.5e-02] [-1.e+15 -5.e-02] Best Regards, Cédric Hannotier