How does NumPy solve nth (5 and higher) degree polynomials?

There is a function in NumPy that solves any polynomial with given coefficient (numpy.roots()). So how does NumPy solve it if there is no formula for 5th and higher degree polynomials?

I know about Newton's method but I wonder how exactly NumPy applies it.

I tried finding information about it in the NumPy documentation and other sources but I did not find anything about that function.

So, the answer is in numpy documentation. But since that was the opportunity for me to play with it, I put in an answer an experiment illustrating how it can be done.

Let's say we want to solve zeros of
X**5 - 7*X**4 + 15*X**3 - 5*X**2 - 16*X + 12

I obtained this polynomial by doing

import sympy
X=sympy.symbols("X")
sympy.expand((X-1)*(X-2)**2*(X-3)*(X+1))

So, spoiler alert, we expect to find -1,1,2,3 as zeros, with 2 being a double one.

For that, we can build a companion matrix

import numpy as np
M=np.array([[0,0,0,0,-12], [1,0,0,0,+16], [0,1,0,0,+5], [0,0,1,0,-15], [0,0,0,1,+7]])
#array([[  0,   0,   0,   0, -12],
#       [  1,   0,   0,   0,  16],
#       [  0,   1,   0,   0,   5],
#       [  0,   0,   1,   0, -15],
#       [  0,   0,   0,   1,   7]])

The eigenvalues of this matrix are the zeros of the polynomial (since the characteristic polynomial of this matrix is the one we are interested in).

We can check this

np.linalg.eigvals(M)
#array([-1.        ,  1.        ,  1.99999992,  2.00000008,  3.        ])

So, that shifted the question:

how those eigenvalues are computed?

Since, at school, the way we learn to compute eigenvalues is... by solving zeros of characteristic polynomial!

So, obviously here that can't be the solution, or else it wouldn't help a lot.

There are several method.

For example, we could use gram-schmidt to compute both Q and R, such as M=QR

def QR(M):
    Q=np.zeros_like(M)
    R=np.zeros_like(M)
    for i in range(len(M)):
        ei=M[:,i].copy()
        for j in range(i):
            R[j,i] = Q[:,j]@M[:,i]
            ei -= R[j,i]/(Q[:,j]@Q[:,j]) * Q[:,j]
        R[i,i] = np.sqrt(ei@ei)
        Q[:,i] = ei / R[i,i]
    return Q,R

Then the iterative algorithm to compute the eigenvalues

A=M.copy()
for i in range(100): # That's an arbitrary number of iterations. In real life
        # we check for convergence. But we aren't in real life: in real life we
        # we don't code a QR decomposition in pure python :D
    Q,R=QR(A)
    A=R@Q


np.diag(A)
#array([ 3.        ,  2.02134245,  1.97865755, -1.01636249,  1.01636249])

Result are the eigenvalues of M, that is the zeros of polynomial.