Plotting polynomials roots

NOTE: See follow up to this question here

I created a simple python script to plot quadratic, cubic and quartic polynomials with integer coefficients between -4 and 4. It uses numpy to find the roots for the polynomials and matplotlib for the actual plotting of the points. Because I want to plot all possible polynomial with integer coefficients between -4 and 4 I figured I needed to do a recursive function with a loop. However if I'm not mistaken this will give a time complexity of $$\O(n!)\$$.

import matplotlib

# Use Qt4Agg because otherwise matplotlib crashes on my computer running Arch Linux
matplotlib.use('Qt4Agg')
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import numpy as np

min_degree = 2
max_degree = 3
colours = []

def polynomial_scatter(ax, min_x, max_x, min_degree, max_degree, cur_degree, coeff):
global colours

if cur_degree <= max_degree:
for c in range(min_x, max_x):
new_coeff = coeff + [c]
polynomial_scatter(ax, min_x, max_x, min_degree, max_degree, cur_degree + 1, new_coeff)
if cur_degree >= min_degree:
roots = np.roots(new_coeff)

# check if no possible solution exists
if len(roots) == 0:
continue

# the real part of the roots will be our x values, the imaginary parts will be our y values
x, y = zip(*[(t.real, t.imag) for t in roots])

# scatter the roots for this polynomial with a colour corresponding to the degree of the polynomial
ax.scatter(x, y, c=colours[cur_degree - max_degree], alpha=.5, s=1)

if __name__ == "__main__":
f, ax = plt.subplots()
colours = cm.rainbow(np.linspace(0, 1, max_degree - min_degree + 1))
polynomial_scatter(ax, -4, 4, min_degree, max_degree, 0, [])
plt.savefig("plot.svg", format="svg")
plt.show()

Here are some pictures showing the result:

if I set max_degree to 5 it takes between 7-15 hours on my computer which is very long. Perhaps there is a way to implement multithreading.

How can I improve the performance of the polynomial_scatter algorithm?

• I think there are some symmetries in the coefficient space you're exploring. If all coefficients are multiplied by a constant, the roots stay the same, right? This would mean that any solution (except for the trivial solution) is plotted at least twice; e.g. after negating coefficients. So you could choose to loop over range(0,max_x) for the last coefficient. Also I think that coeff[0]==0 implies a lower degree, correct? If so, the function should either call itself or calculate roots.
– mvds
Jun 2, 2017 at 0:11
• Expanding on the previous comment: if the last coefficient is 0, you're actually solving a lower degree polynomial (ignoring the trivial solution). That implies that you could also loop over range(1,max_x). This should be 55% faster already.
– mvds
Jun 2, 2017 at 0:14
• @mvds that is an interesting optimization. Not sure how to implement that into a recursive function however. Jun 2, 2017 at 8:06
• in pseudo code: range(cur_degree==max_degree?1:min_x, max_x)
– mvds
Jun 2, 2017 at 15:06
• @mvds I might implement that, but I switched to C++ and have another performance problem now. Here is the new post: codereview.stackexchange.com/questions/164790/… Jun 2, 2017 at 15:21